00001
00002
00003
00004
00005
00006
00007
00008 #include <cstdio>
00009 #include <iostream>
00010 #include <fstream>
00011
00012 #include "CLHEP/Vector/LorentzVector.h"
00013 #include "CLHEP/Vector/ThreeVector.h"
00014 #include "CLHEP/Vector/LorentzVector.h"
00015 #include "CLHEP/Vector/TwoVector.h"
00016 #include "CLHEP/Geometry/Point3D.h"
00017 using CLHEP::Hep3Vector;
00018 using CLHEP::Hep2Vector;
00019 using CLHEP::HepLorentzVector;
00020
00021 #include "GaudiKernel/IHistogramSvc.h"
00022 #include "GaudiKernel/MsgStream.h"
00023 #include "GaudiKernel/AlgFactory.h"
00024 #include "GaudiKernel/ISvcLocator.h"
00025 #include "GaudiKernel/SmartDataPtr.h"
00026 #include "GaudiKernel/SmartDataLocator.h"
00027 #include "GaudiKernel/IDataProviderSvc.h"
00028 #include "GaudiKernel/Bootstrap.h"
00029 #include "GaudiKernel/PropertyMgr.h"
00030 #include "GaudiKernel/INTupleSvc.h"
00031 #include "GaudiKernel/NTuple.h"
00032 #include "GaudiKernel/Algorithm.h"
00033
00034 #include "Identifier/Identifier.h"
00035 #include "Identifier/MucID.h"
00036 #include "EventModel/EventHeader.h"
00037 #include "EventModel/EventModel.h"
00038 #include "EvTimeEvent/RecEsTime.h"
00039
00040 #include "McTruth/McKine.h"
00041 #include "McTruth/McParticle.h"
00042 #include "McTruth/MucMcHit.h"
00043
00044 #include "MdcRecEvent/RecMdcTrack.h"
00045 #include "EmcRecEventModel/RecEmcEventModel.h"
00046 #include "TofRecEvent/RecTofTrack.h"
00047 #include "DstEvent/TofHitStatus.h"
00048
00049
00050 #include "MucRawEvent/MucDigi.h"
00051 #include "MucRecEvent/MucRecHit.h"
00052 #include "MucRecEvent/RecMucTrack.h"
00053 #include "MucRecEvent/MucRecHitContainer.h"
00054
00055 #include "MucCalibAlg/MucStructConst.h"
00056 #include "MucCalibAlg/MucCalibMgr.h"
00057 #include "MucCalibAlg/MucMark.h"
00058 #include "MucCalibAlg/MucBoxCal.h"
00059 #include "MucCalibAlg/MucStripCal.h"
00060
00061 using namespace std;
00062 using namespace Event;
00063
00064
00065 MucCalibMgr::MucCalibMgr( std::vector<double> jobInfo, std::vector<int> configInfo, std::string outputFile )
00066 {
00067 MsgStream log(msgSvc, "MucCalibMgr");
00068 log << MSG::INFO << "-------Initializing------- " << endreq;
00069
00070
00071 log << MSG::INFO << "Initialize Gaudi service" << endreq;
00072 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00073 ntupleSvc = NULL;
00074 eventSvc = NULL;
00075
00076 m_jobStart = clock();
00077
00078
00079 log << MSG::INFO << "Initialize parameters" << endreq;
00080
00081 m_vs = jobInfo[0]; m_hv = jobInfo[1]; m_th = jobInfo[2];
00082
00083 if( jobInfo[3] <= 0 ) m_er = TRIGGER_RATE;
00084 else m_er = jobInfo[3];
00085
00086 if( jobInfo[4] <= 0 || jobInfo[4] >5 ) m_tag = 0;
00087 else m_tag = jobInfo[4];
00088
00089 if( configInfo[0] < 0 || configInfo[0] > 2 )
00090 m_recMode = 0;
00091 else
00092 m_recMode = configInfo[0];
00093
00094 m_usePad = configInfo[1];
00095
00096 if( configInfo[2] < 0 || configInfo[2] > STRIP_INBOX_MAX )
00097 m_effWindow = EFF_WINDOW;
00098 else m_effWindow = configInfo[2];
00099
00100 if( configInfo[3]<0 || configInfo[3] >4 )
00101 m_clusterMode = DEFAULT_BUILD_MODE;
00102 else m_clusterMode = configInfo[3];
00103
00104 m_clusterSave = configInfo[4];
00105 m_checkEvent = configInfo[5];
00106 m_dimuSelect = configInfo[6];
00107 m_dimuOnly = configInfo[7];
00108
00109 m_outputFile = outputFile;
00110
00111 m_fdata = NULL;
00112 if( m_clusterMode ==1 && m_clusterSave == 1 )
00113 m_fdata = new ofstream("MucCluster.dat", ios::out);
00114
00115 m_fStartRun = m_fEndRun = 0;
00116 m_currentRun = m_currentEvent = 0;
00117 m_eventTag = 0;
00118
00119 for( int i=0; i<PART_MAX; i++ )
00120 for( int j=0; j<SEGMENT_MAX; j++ )
00121 for( int k=0; k<LAYER_MAX; k++ )
00122 for( int n=0; n<STRIP_INBOX_MAX; n++ )
00123 {
00124 m_record[i][j][k][n][0] = 0;
00125 m_record[i][j][k][n][1] = 0;
00126 m_record[i][j][k][n][2] = 0;
00127 m_record[i][j][k][n][3] = 0;
00128 m_record[i][j][k][n][4] = 0;
00129 }
00130
00131 for( int i=0; i<LAYER_MAX; i++ )
00132 {
00133 m_layerResults[0][i] = 0;
00134 m_layerResults[1][i] = 0;
00135 m_layerResults[2][i] = 0;
00136 m_layerResults[3][i] = 0;
00137 m_layerResults[4][i] = 0;
00138 }
00139
00140 for( int i=0; i<BOX_MAX; i++ )
00141 {
00142 m_boxResults[0][i] = 0;
00143 m_boxResults[1][i] = 0;
00144 m_boxResults[2][i] = 0;
00145 m_boxResults[3][i] = 0;
00146 m_boxResults[4][i] = 0;
00147 }
00148
00149 for( int i=0; i<STRIP_MAX; i++ )
00150 {
00151 m_stripResults[0][i] = 0;
00152 m_stripResults[1][i] = 0;
00153 m_stripResults[2][i] = 0;
00154 m_stripResults[3][i] = 0;
00155
00156 }
00157
00158 m_fTotalDAQTime = 0;
00159
00160 m_ptrMucMark = new MucMark(0,0,0,0);
00161 m_ptrIdTr = new MucIdTransform();
00162
00163 m_digiCol.resize(0);
00164 m_calHitCol.resize(0);
00165 m_expHitCol.resize(0);
00166 m_effHitCol.resize(0);
00167 m_nosHitCol.resize(0);
00168 m_clusterCol.resize(0);
00169
00170 for( int i=0; i<PART_MAX; i++ )
00171 for( int j=0; j<SEGMENT_MAX; j++ ) {
00172 m_segDigiCol[i][j].resize(0);
00173 }
00174
00175
00176 InitNtuple();
00177
00178
00179 InitHisto();
00180
00181
00182 InitConstTree();
00183
00184
00185 InitArea();
00186 log << MSG::INFO << "Initialize done!" << endreq;
00187 }
00188
00189
00190 MucCalibMgr::~MucCalibMgr()
00191 {
00192 ClearOnlineHisto();
00193
00194 ClearHistoLV2();
00195 ClearHistoLV1();
00196 ClearHistoLV0();
00197 if( m_usePad != 0 ) Clear2DEffHisto();
00198
00199 ClearClusterHisto();
00200 ClearResHisto();
00201 ClearConstTree();
00202
00203 if(m_ptrMucMark) delete m_ptrMucMark;
00204 if(m_ptrIdTr) delete m_ptrIdTr;
00205
00206 m_digiCol.clear();
00207 m_calHitCol.clear();
00208 m_effHitCol.clear();
00209 m_nosHitCol.clear();
00210 m_clusterCol.clear();
00211
00212 for( int i=0; i<PART_MAX; i++ )
00213 for( int j=0; j<SEGMENT_MAX; j++ ) {
00214 m_segDigiCol[i][j].clear();
00215 }
00216
00217 if(m_record) delete []m_record;
00218 if(m_layerResults) delete []m_layerResults;
00219 if(m_boxResults) delete []m_boxResults;
00220 if(m_stripResults) delete []m_stripResults;
00221 }
00222
00223
00224 StatusCode MucCalibMgr::ClearOnlineHisto()
00225 {
00226 if(m_hHitMapBarrel_Lay) delete []m_hHitMapBarrel_Lay;
00227 if(m_hHitMapEndcap_Lay) delete []m_hHitMapEndcap_Lay;
00228 if(m_hHitMapBarrel_Seg) delete []m_hHitMapBarrel_Seg;
00229 if(m_hHitMapEndcap_Seg) delete []m_hHitMapEndcap_Seg;
00230
00231 if(m_hHitVsEvent) delete m_hHitVsEvent;
00232 if(m_hTrackDistance) delete m_hTrackDistance;
00233 if(m_hTrackPosPhiDiff) delete m_hTrackPosPhiDiff;
00234 if(m_hTrackPosThetaDiff)delete m_hTrackPosThetaDiff;
00235 if(m_hTrackMomPhiDiff) delete m_hTrackMomPhiDiff;
00236 if(m_hTrackMomThetaDiff)delete m_hTrackMomThetaDiff;
00237 if(m_hDimuTracksPosDiff)delete m_hDimuTracksPosDiff;
00238 if(m_hDimuTracksMomDiff)delete m_hDimuTracksMomDiff;
00239 if(m_hPhiCosTheta) delete m_hPhiCosTheta;
00240
00241 return StatusCode::SUCCESS;
00242 }
00243
00244
00245 StatusCode MucCalibMgr::ClearHistoLV0()
00246 {
00247 if(m_hBrLayerFire) delete m_hBrLayerFire;
00248 if(m_hEcLayerFire) delete m_hEcLayerFire;
00249 if(m_hLayerFire) delete m_hLayerFire;
00250 if(m_hLayerExpHit) delete m_hLayerExpHit;
00251 if(m_hLayerExpHit) delete m_hLayerEffHit;
00252 if(m_hLayerNosHit) delete m_hLayerNosHit;
00253 if(m_hLayerEff) delete m_hLayerEff;
00254 if(m_hLayerNosRatio)delete m_hLayerNosRatio;
00255 if(m_hLayerArea) delete m_hLayerArea;
00256 if(m_hLayerNos) delete m_hLayerNos;
00257 if(m_hLayerCnt) delete m_hLayerCnt;
00258
00259 return StatusCode::SUCCESS;
00260 }
00261
00262
00263 StatusCode MucCalibMgr::ClearHistoLV1()
00264 {
00265 if(m_hBoxFire) delete m_hBoxFire;
00266 if(m_hBoxExpHit) delete m_hBoxExpHit;
00267 if(m_hBoxEffHit) delete m_hBoxEffHit;
00268 if(m_hBoxNosHit) delete m_hBoxNosHit;
00269 if(m_hBoxEff) delete m_hBoxEff;
00270 if(m_hBoxNosRatio)delete m_hBoxNosRatio;
00271 if(m_hBoxArea) delete m_hBoxArea;
00272 if(m_hBoxNos) delete m_hBoxNos;
00273 if(m_hBoxCnt) delete m_hBoxCnt;
00274
00275 return StatusCode::SUCCESS;
00276 }
00277
00278
00279 StatusCode MucCalibMgr::ClearHistoLV2()
00280 {
00281 for( int i=0; i< BOX_MAX; i++ )
00282 {
00283 if(m_hStripFireMap[i]) delete m_hStripFireMap[i];
00284 if(m_hStripExpHitMap[i]) delete m_hStripExpHitMap[i];
00285 if(m_hStripEffHitMap[i]) delete m_hStripEffHitMap[i];
00286 if(m_hStripNosHitMap[i]) delete m_hStripNosHitMap[i];
00287 if(m_hStripEffMap[i]) delete m_hStripEffMap[i];
00288 if(m_hStripNosRatioMap[i])delete m_hStripNosRatioMap[i];
00289 }
00290
00291 if(m_hStripFire) delete m_hStripFire;
00292 if(m_hStripExpHit) delete m_hStripExpHit;
00293 if(m_hStripEffHit) delete m_hStripEffHit;
00294 if(m_hStripNosHit) delete m_hStripNosHit;
00295 if(m_hStripEff) delete m_hStripEff;
00296 if(m_hStripNosRatio)delete m_hStripNosRatio;
00297 if(m_hStripArea) delete m_hStripArea;
00298 if(m_hStripNos) delete m_hStripNos;
00299 if(m_hStripCnt) delete m_hStripCnt;
00300
00301 return StatusCode::SUCCESS;
00302 }
00303
00304
00305 StatusCode MucCalibMgr::Clear2DEffHisto()
00306 {
00307
00308 for( int i=0; i<PART_MAX; i++ )
00309 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
00310 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
00311 {
00312 if(m_h2DExpMap[i][j][k]) delete m_h2DExpMap[i][j][k];
00313 if(m_h2DHitMap[i][j][k]) delete m_h2DHitMap[i][j][k];
00314 if(m_h2DEffMap[i][j][k]) delete m_h2DEffMap[i][j][k];
00315 }
00316
00317 return StatusCode::SUCCESS;
00318 }
00319
00320
00321 StatusCode MucCalibMgr::ClearClusterHisto()
00322 {
00323 for( int i=0; i<LAYER_MAX; i++ ) {
00324 if(m_hLayerCluster[i]) delete m_hLayerCluster[i];
00325 }
00326
00327 for( int i=0; i<BOX_MAX; i++ ) {
00328 if(m_hBoxCluster[i]) delete m_hBoxCluster[i];
00329 }
00330
00331 if(m_hLayerClusterCmp) delete m_hLayerClusterCmp;
00332 if(m_hBoxClusterCmp) delete m_hBoxClusterCmp;
00333
00334 return StatusCode::SUCCESS;
00335 }
00336
00337
00338 StatusCode MucCalibMgr::ClearResHisto()
00339 {
00340 for( int i=0; i<B_LAY_NUM; i++ ) {
00341 if(m_hBarrelResDist[i]) delete m_hBarrelResDist[i];
00342 }
00343 for( int i=0; i<E_LAY_NUM; i++ ) {
00344 if(m_hEndcapResDist[i]) delete m_hEndcapResDist[i];
00345 }
00346
00347 if(m_hBarrelResComp[0]) delete m_hBarrelResComp[0];
00348 if(m_hBarrelResComp[1]) delete m_hBarrelResComp[1];
00349 if(m_hEndcapResComp[0]) delete m_hEndcapResComp[0];
00350 if(m_hEndcapResComp[1]) delete m_hEndcapResComp[1];
00351
00352 return StatusCode::SUCCESS;
00353 }
00354
00355
00356 StatusCode MucCalibMgr::ClearConstTree()
00357 {
00358
00359 if(m_tJobLog) delete m_tJobLog;
00360 if(m_tStatLog) delete m_tStatLog;
00361 if(m_tLayConst) delete m_tLayConst;
00362 if(m_tBoxConst) delete m_tBoxConst;
00363 if(m_tStrConst) delete m_tStrConst;
00364
00365 return StatusCode::SUCCESS;
00366 }
00367
00368
00369
00370
00371
00372
00373
00374 StatusCode MucCalibMgr::InitNtuple()
00375 {
00376 MsgStream log(msgSvc, "MucCalibMgr");
00377 log << MSG::INFO << "Initialize NTuples" << endreq;
00378
00379
00380 Gaudi::svcLocator() -> service("NTupleSvc", ntupleSvc);
00381
00382 StatusCode sc;
00383
00384
00385 NTuplePtr nt1(ntupleSvc, "FILE450/EventLog");
00386 if ( nt1 ) { m_eventLogTuple = nt1; }
00387 else
00388 {
00389 m_eventLogTuple = ntupleSvc->book ("FILE450/EventLog", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
00390 if ( m_eventLogTuple )
00391 {
00392 sc = m_eventLogTuple->addItem ("event_id", m_ntEventId);
00393 sc = m_eventLogTuple->addItem ("event_tag", m_ntEventTag);
00394 sc = m_eventLogTuple->addItem ("start_time", m_ntEsTime);
00395 sc = m_eventLogTuple->addItem ("digi_num", m_ntDigiNum);
00396 sc = m_eventLogTuple->addItem ("track_num", m_ntTrackNum);
00397 sc = m_eventLogTuple->addItem ("exphit_num", m_ntExpHitNum);
00398 sc = m_eventLogTuple->addItem ("effhit_num", m_ntEffHitNum);
00399 sc = m_eventLogTuple->addItem ("noshit_num", m_ntNosHitNum);
00400 sc = m_eventLogTuple->addItem ("cluster_num", m_ntClusterNum);
00401 sc = m_eventLogTuple->addItem ("event_time", m_ntEventTime);
00402 }
00403 else {
00404 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_eventLogTuple) << endmsg;
00405 return StatusCode::FAILURE;
00406 }
00407 }
00408
00409
00410
00411 NTuplePtr nt2(ntupleSvc, "FILE450/MdcTrkInfo");
00412 if ( nt2 ) { m_mdcTrkInfoTuple = nt2; }
00413 else
00414 {
00415 m_mdcTrkInfoTuple = ntupleSvc->book ("FILE450/MdcTrkInfo", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
00416 if ( m_mdcTrkInfoTuple )
00417 {
00418 sc = m_mdcTrkInfoTuple->addItem ("charge", m_charge);
00419 sc = m_mdcTrkInfoTuple->addItem ("mdcpx", m_mdcpx);
00420 sc = m_mdcTrkInfoTuple->addItem ("mdcpy", m_mdcpy);
00421 sc = m_mdcTrkInfoTuple->addItem ("mdcpz", m_mdcpz);
00422 sc = m_mdcTrkInfoTuple->addItem ("mdcpt", m_mdcpt);
00423 sc = m_mdcTrkInfoTuple->addItem ("mdcpp", m_mdcpp);
00424 sc = m_mdcTrkInfoTuple->addItem ("mdcphi", m_mdcphi);
00425 sc = m_mdcTrkInfoTuple->addItem ("mdctheta", m_mdctheta);
00426 }
00427 else {
00428 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_mdcTrkInfoTuple) << endmsg;
00429 return StatusCode::FAILURE;
00430 }
00431 }
00432
00433
00434 NTuplePtr nt3(ntupleSvc, "FILE450/TrackInfo");
00435 if ( nt3 ) { m_trackInfoTuple = nt3; }
00436 else
00437 {
00438 m_trackInfoTuple = ntupleSvc->book ("FILE450/TrackInfo", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
00439 if ( m_trackInfoTuple )
00440 {
00441 sc = m_trackInfoTuple->addItem ("track_event", m_ntTrackEvent);
00442 sc = m_trackInfoTuple->addItem ("track_tag", m_ntTrackTag);
00443 sc = m_trackInfoTuple->addItem ("track_hits", m_ntTrackHits);
00444 sc = m_trackInfoTuple->addItem ("segment_fly", m_ntTrackSegFly);
00445 sc = m_trackInfoTuple->addItem ("layer_fly_a", m_ntTrackLayFlyA);
00446 sc = m_trackInfoTuple->addItem ("layer_fly_b", m_ntTrackLayFlyB);
00447 sc = m_trackInfoTuple->addItem ("layer_fly_c", m_ntTrackLayFlyC);
00448 sc = m_trackInfoTuple->addItem ("rec_mode", m_trkRecMode);
00449 sc = m_trackInfoTuple->addItem ("chi2", m_chi2);
00450 sc = m_trackInfoTuple->addItem ("px", m_px);
00451 sc = m_trackInfoTuple->addItem ("py", m_py);
00452 sc = m_trackInfoTuple->addItem ("pz", m_pz);
00453 sc = m_trackInfoTuple->addItem ("pt", m_pt);
00454 sc = m_trackInfoTuple->addItem ("pp", m_pp);
00455 sc = m_trackInfoTuple->addItem ("r", m_r );
00456 sc = m_trackInfoTuple->addItem ("costheta", m_cosTheta);
00457 sc = m_trackInfoTuple->addItem ("theta", m_theta);
00458 sc = m_trackInfoTuple->addItem ("phi", m_phi);
00459 sc = m_trackInfoTuple->addItem ("depth", m_depth);
00460 sc = m_trackInfoTuple->addItem ("br_last_lay", m_brLastLayer);
00461 sc = m_trackInfoTuple->addItem ("ec_last_lay", m_ecLastLayer);
00462 sc = m_trackInfoTuple->addItem ("total_hits", m_totalHits);
00463 sc = m_trackInfoTuple->addItem ("fired_layers", m_totalLayers);
00464 sc = m_trackInfoTuple->addItem ("maxhits_in_layer", m_maxHitsInLayer);
00465 }
00466 else {
00467 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_trackInfoTuple) << endmsg;
00468 return StatusCode::FAILURE;
00469 }
00470 }
00471
00472
00473 NTuplePtr nt4(ntupleSvc, "FILE450/TrackDiff");
00474 if ( nt4 ) { m_trackDiffTuple = nt4; }
00475 else
00476 {
00477 m_trackDiffTuple = ntupleSvc->book ("FILE450/TrackDiff", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
00478 if ( m_trackDiffTuple ) {
00479 sc = m_trackDiffTuple->addItem ("dimu_tag", m_ntDimuTag);
00480 sc = m_trackDiffTuple->addItem ("pos_phi_diff", m_ntPosPhiDiff);
00481 sc = m_trackDiffTuple->addItem ("pos_theta_diff", m_ntPosThetaDiff);
00482 sc = m_trackDiffTuple->addItem ("mom_phi_diff", m_ntMomPhiDiff);
00483 sc = m_trackDiffTuple->addItem ("mom_theta_diff", m_ntMomThetaDiff);
00484 }
00485 else {
00486 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_trackDiffTuple) << endmsg;
00487 return StatusCode::FAILURE;
00488 }
00489 }
00490
00491
00492 NTuplePtr nt5(ntupleSvc, "FILE450/ClusterSize");
00493 if ( nt5 ) { m_clusterSizeTuple = nt5; }
00494 else
00495 {
00496 m_clusterSizeTuple = ntupleSvc->book ("FILE450/ClusterSize", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
00497 if ( m_clusterSizeTuple ) {
00498 sc = m_clusterSizeTuple->addItem ("cluster_size", m_ntClusterSize);
00499 }
00500 else {
00501 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_clusterSizeTuple) << endmsg;
00502 return StatusCode::FAILURE;
00503 }
00504 }
00505
00506
00507 NTuplePtr nt6(ntupleSvc, "FILE450/EffWindow");
00508 if ( nt6 ) { m_effWindowTuple = nt6; }
00509 else
00510 {
00511 m_effWindowTuple = ntupleSvc->book ("FILE450/EffWindow", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
00512 if ( m_effWindowTuple ) {
00513 sc = m_effWindowTuple->addItem ("hit_window", m_ntEffWindow);
00514 }
00515 else {
00516 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_effWindowTuple) << endmsg;
00517 return StatusCode::FAILURE;
00518 }
00519 }
00520
00521 NTuplePtr nt7(ntupleSvc, "FILE450/ResInfo");
00522 if ( nt7 ) { m_resInfoTuple = nt7; }
00523 else
00524 {
00525 m_resInfoTuple = ntupleSvc->book ("FILE450/ResInfo", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
00526 if ( m_resInfoTuple ) {
00527 sc = m_resInfoTuple->addItem ("line_res", m_lineRes);
00528 sc = m_resInfoTuple->addItem ("quad_res", m_quadRes);
00529 sc = m_resInfoTuple->addItem ("extr_res", m_extrRes);
00530 sc = m_resInfoTuple->addItem ("res_part", m_resPart);
00531 sc = m_resInfoTuple->addItem ("res_segment", m_resSegment);
00532 sc = m_resInfoTuple->addItem ("res_layer", m_resLayer);
00533 sc = m_resInfoTuple->addItem ("res_fired", m_resFired);
00534 sc = m_resInfoTuple->addItem ("res_mode", m_resMode);
00535 }
00536 else {
00537 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_resInfoTuple) << endmsg;
00538 return StatusCode::FAILURE;
00539 }
00540 }
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572 return StatusCode::SUCCESS;
00573 }
00574
00575
00576 StatusCode MucCalibMgr::InitHisto()
00577 {
00578 MsgStream log(msgSvc, "MucCalibMgr");
00579 log << MSG::INFO << "Initialize Histograms" << endreq;
00580
00581
00582 InitOnlineHisto();
00583
00584
00585 InitHistoLV2();
00586 InitHistoLV1();
00587 InitHistoLV0();
00588
00589
00590 if( m_usePad != 0 ) Init2DEffHisto();
00591
00592
00593 InitClusterHisto();
00594
00595
00596 InitResHisto();
00597
00598 return StatusCode::SUCCESS;
00599 }
00600
00601
00602 StatusCode MucCalibMgr::InitOnlineHisto()
00603 {
00604 char name[60];
00605 char title[60];
00606 int stripMax = 0;
00607
00608
00609 for( int i=0; i<B_LAY_NUM; i++ )
00610 {
00611
00612 if(i%2 != 0 ) { stripMax = 96*7 + 112; }
00613 else { stripMax = 48*8; }
00614 sprintf( name, "HitMapBarrel_Lay_L%d", i );
00615 sprintf( title, "Hit map in barrel layer %d", i );
00616 m_hHitMapBarrel_Lay[i] = new TH1F(name,title, stripMax, 0, stripMax);
00617
00618 if( i<E_LAY_NUM )
00619 {
00620 stripMax = 64*4;
00621 sprintf( name, "HitMapEastEndcap_L%d", i );
00622 sprintf( title, "Hit map in east-endcap layer %d", i );
00623 m_hHitMapEndcap_Lay[0][i] = new TH1F(name,title, stripMax, 0, stripMax);
00624
00625 sprintf( name, "HitMapWestEndcap_L%d", i );
00626 sprintf( title, "Hit map in west-endcap layer %d", i );
00627 m_hHitMapEndcap_Lay[1][i] = new TH1F(name,title, stripMax, 0, stripMax);
00628 }
00629
00630 }
00631
00632 for( int i=0; i<B_SEG_NUM; i++ )
00633 {
00634
00635 sprintf( name, "HitMapBarrel_Seg_S%d", i );
00636 sprintf( title, "Hit map in barrel segment %d", i );
00637 if(i==2 ) { stripMax = 48*5 + 112*4; }
00638 else { stripMax = 48*5 + 96*4; }
00639 m_hHitMapBarrel_Seg[i] = new TH1F(name,title, stripMax, 0, stripMax);
00640
00641 if( i<E_SEG_NUM )
00642 {
00643 sprintf( name, "HitMapEastEndcap_S%d", i );
00644 sprintf( title, "Hit map in east-endcap segment %d", i );
00645 stripMax = 64*8;
00646 m_hHitMapEndcap_Seg[0][i] = new TH1F(name,title, stripMax, 0, stripMax);
00647
00648 sprintf( name, "HitMapWestEndcap_S%d", i );
00649 sprintf( title, "Hit map in west-endcap segment %d", i );
00650 m_hHitMapEndcap_Seg[1][i] = new TH1F(name,title, stripMax, 0, stripMax);
00651 }
00652 }
00653
00654
00655
00656 m_hHitVsEvent = new TH1F("HitVsEvent","Hit VS event",10000,0,10000);
00657 m_hHitVsEvent->GetXaxis()->SetTitle("Event NO.");
00658 m_hHitVsEvent->GetXaxis()->CenterTitle();
00659 m_hHitVsEvent->GetYaxis()->SetTitle("Hits");
00660 m_hHitVsEvent->GetYaxis()->CenterTitle();
00661
00662 m_hTrackDistance = new TH1F("TrackDistance","Track distance", 1000,-500,500);
00663 m_hTrackDistance->GetXaxis()->SetTitle("Distance of fit pos and hit pos on 1st layer [cm]");
00664 m_hTrackDistance->GetXaxis()->CenterTitle();
00665
00666 m_hTrackPosPhiDiff = new TH1F("TrackPosPhiDiff","#Delta#phi of tracks pos", 720,-2*PI,2*PI);
00667 m_hTrackPosPhiDiff->GetXaxis()->SetTitle("#Delta#phi [rad]");
00668 m_hTrackPosPhiDiff->GetXaxis()->CenterTitle();
00669
00670 m_hTrackPosThetaDiff = new TH1F("TrackPosThetaDiff","#Delta#theta of tracks pos", 720,-2*PI,2*PI);
00671 m_hTrackPosThetaDiff->GetXaxis()->SetTitle("#Delta#theta [rad]");
00672 m_hTrackPosThetaDiff->GetXaxis()->CenterTitle();
00673
00674 m_hTrackMomPhiDiff = new TH1F("TrackMomPhiDiff","#Delta#phi of tracks mom", 720,-2*PI,2*PI);
00675 m_hTrackMomPhiDiff->GetXaxis()->SetTitle("#Delta#phi [rad]");
00676 m_hTrackMomPhiDiff->GetXaxis()->CenterTitle();
00677
00678 m_hTrackMomThetaDiff = new TH1F("TrackMomThetaDiff","#Delta#theta of tracks mom", 720,-2*PI,2*PI);
00679 m_hTrackMomThetaDiff->GetXaxis()->SetTitle("#Delta#theta [rad]");
00680 m_hTrackMomThetaDiff->GetXaxis()->CenterTitle();
00681
00682
00683 m_hDimuTracksPosDiff = new TH2F("DimuTracksPosDiff", "#Delta#phi VS #Delta#theta of dimu tracks pos", 720, -PI, PI, 720, -PI, PI);
00684 m_hDimuTracksPosDiff->GetXaxis()->SetTitle("#Delta#theta");
00685 m_hDimuTracksPosDiff->GetXaxis()->CenterTitle();
00686 m_hDimuTracksPosDiff->GetYaxis()->SetTitle("#Delta#phi");
00687 m_hDimuTracksPosDiff->GetYaxis()->CenterTitle();
00688
00689 m_hDimuTracksMomDiff = new TH2F("DimuTracksMomDiff", "#Delta#phi VS #Delta#theta of dimu tracks mom", 720, -PI, PI, 720, -PI, PI);
00690 m_hDimuTracksMomDiff->GetXaxis()->SetTitle("#Delta#theta");
00691 m_hDimuTracksMomDiff->GetXaxis()->CenterTitle();
00692 m_hDimuTracksMomDiff->GetYaxis()->SetTitle("#Delta#phi");
00693 m_hDimuTracksMomDiff->GetYaxis()->CenterTitle();
00694
00695 m_hPhiCosTheta = new TH2F("PhiVsCosTheta", "#phi VS cos(#theta)", 720, -1, 1, 720, -PI, PI);
00696 m_hPhiCosTheta->GetXaxis()->SetTitle("cos(#theta)");
00697 m_hPhiCosTheta->GetXaxis()->CenterTitle();
00698 m_hPhiCosTheta->GetYaxis()->SetTitle("#phi");
00699 m_hPhiCosTheta->GetYaxis()->CenterTitle();
00700
00701 return StatusCode::SUCCESS;
00702 }
00703
00704
00705 StatusCode MucCalibMgr::InitHistoLV0()
00706 {
00707 m_hBrLayerFire = new TH1F("BrLayerFire","Fires per layer in Barrel", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00708 m_hEcLayerFire = new TH1F("EcLayerFire","Fires per layer in Endcap", 2*LAYER_MAX-1, -LAYER_MAX+1, LAYER_MAX-0.5);
00709
00710 m_hLayerFire = new TH1F("LayerFire","Fires per layer", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00711 m_hLayerExpHit = new TH1F("LayerExpHit","Exp hits per layer", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00712 m_hLayerEffHit = new TH1F("LayerEffHit","Eff hits per layer", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00713 m_hLayerNosHit = new TH1F("LayerNosHit","Nos hits per layer", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00714 m_hLayerEff = new TH1F("LayerEff","Layer efficiency", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00715 m_hLayerNosRatio= new TH1F("LayerNosRatio","Layer noise hit ratio", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00716 m_hLayerArea = new TH1F("LayerArea","Layer area", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00717 m_hLayerNos = new TH1F("LayerNos","Layer noise", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00718 m_hLayerCnt = new TH1F("LayerCnt","Layer count", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00719
00720 return StatusCode::SUCCESS;
00721 }
00722
00723
00724 StatusCode MucCalibMgr::InitHistoLV1()
00725 {
00726 m_hBoxFire = new TH1F("BoxFire","Fires per box", BOX_MAX+1, -0.5, BOX_MAX+0.5);
00727 m_hBoxExpHit = new TH1F("BoxExpHit","Exp hits per box", BOX_MAX+1, -0.5, BOX_MAX+0.5);
00728 m_hBoxEffHit = new TH1F("BoxEffHit","Eff hits per box", BOX_MAX+1, -0.5, BOX_MAX+0.5);
00729 m_hBoxNosHit = new TH1F("BoxNosHit","Nos hits per box", BOX_MAX+1, -0.5, BOX_MAX+0.5);
00730 m_hBoxEff = new TH1F("BoxEff","Box efficiency", BOX_MAX+1, -0.5, BOX_MAX+0.5);
00731 m_hBoxNosRatio = new TH1F("BoxNosRatio","Box noise hit ratio", BOX_MAX+1, -0.5, BOX_MAX+0.5);
00732 m_hBoxArea = new TH1F("BoxArea","Box area", BOX_MAX+1, -0.5, BOX_MAX+0.5);
00733 m_hBoxNos = new TH1F("BoxNos","Box noise", BOX_MAX+1, -0.5, BOX_MAX+0.5);
00734 m_hBoxCnt = new TH1F("BoxCnt","Box count", BOX_MAX+1, -0.5, BOX_MAX+0.5);
00735
00736 return StatusCode::SUCCESS;
00737 }
00738
00739
00740 StatusCode MucCalibMgr::InitHistoLV2()
00741 {
00742 char name[60];
00743 char title[60];
00744 int part, segment, layer, stripMax;
00745 part = segment = layer = stripMax = 0;
00746
00747 for( int i=0; i<BOX_MAX; i++ )
00748 {
00749 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
00750 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
00751
00752 sprintf( name, "StripFireMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
00753 sprintf( title, "Fires per strip in P%d_S%d_L%d Box%d",part, segment, layer, i );
00754 m_hStripFireMap[i] = new TH1F(name,title, stripMax, 0, stripMax);
00755
00756 sprintf( name, "StripExpHitMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
00757 sprintf( title, "Exp hits per strip in P%d_S%d_L%d Box%d",part, segment, layer, i );
00758 m_hStripExpHitMap[i] = new TH1F(name,title, stripMax, 0, stripMax);
00759
00760 sprintf( name, "StripEffHitMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
00761 sprintf( title, "Eff hits per strip in P%d_S%d_L%d Box%d",part, segment, layer, i );
00762 m_hStripEffHitMap[i] = new TH1F(name,title, stripMax, 0, stripMax);
00763
00764 sprintf( name, "StripNosHitMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
00765 sprintf( title, "Inc hits per strip in P%d_S%d_L%d Box%d",part, segment, layer, i );
00766 m_hStripNosHitMap[i] = new TH1F(name,title, stripMax, 0, stripMax);
00767
00768 sprintf( name, "StripEffMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
00769 sprintf( title, "Strip efficiency in P%d_S%d_L%d Box%d",part, segment, layer, i );
00770 m_hStripEffMap[i] = new TH1F(name,title, stripMax, 0, stripMax);
00771
00772 sprintf( name, "StripNosRatioMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
00773 sprintf( title, "Strip noise hit ratio in P%d_S%d_L%d Box%d",part, segment, layer, i );
00774 m_hStripNosRatioMap[i] = new TH1F(name,title, stripMax, 0, stripMax);
00775 }
00776
00777 m_hStripFire = new TH1F("StripFire", "Fires per strip", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
00778 m_hStripExpHit = new TH1F("StripExpHit", "Exp hit per strip", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
00779 m_hStripEffHit = new TH1F("StripEffHit", "Eff hit per strip", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
00780 m_hStripNosHit = new TH1F("StripNoshit", "Nos hit per strip", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
00781 m_hStripEff = new TH1F("StripEff", "Strip efficiency", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
00782 m_hStripNosRatio= new TH1F("StripNosRatio", "Strip noise hit ratio", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
00783 m_hStripArea = new TH1F("StripArea", "Strip area", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
00784 m_hStripNos = new TH1F("StripNos", "Strip noise", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
00785 m_hStripCnt = new TH1F("StripCnt", "Strip count", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
00786
00787 return StatusCode::SUCCESS;
00788 }
00789
00790
00791 StatusCode MucCalibMgr::Init2DEffHisto()
00792 {
00793 char name[60];
00794 int stripMax, boxID, stripID, xBin, yBin;
00795 stripMax = boxID = stripID = xBin = yBin = 0;
00796 double xStart, xEnd, yStart, yEnd, sWidth;
00797 xStart = xEnd = yStart = yEnd = sWidth = 0.;
00798
00799 m_histArray = new TObjArray();
00800
00801 for( int i=0; i<PART_MAX; i++ )
00802 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
00803 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
00804 {
00805 boxID = m_ptrIdTr->GetBoxId(i,j,k);
00806
00807 if( i==BRID )
00808 {
00809 xStart = -2000, xEnd = 2000;
00810 sWidth = B_STR_DST[k];
00811 xBin = int((xEnd - xStart)/sWidth);
00812 yStart = -B_BOX_WT[k]/2-100, yEnd = B_BOX_WT[k]/2+100;
00813 yBin = int((B_BOX_WT[k]+200)/sWidth);
00814 }
00815 else
00816 {
00817 xStart = yStart = -1250;
00818 xEnd = yEnd = 1250;
00819 sWidth = E_STR_DST;
00820 xBin = yBin = int((xEnd - xStart)/sWidth);
00821 }
00822
00823 sprintf(name, "ExpMap2D_P%d_S%d_L%d_Box%d", i, j, k, boxID);
00824 m_h2DExpMap[i][j][k] = new TH2F(name, name, xBin, xStart, xEnd, yBin, yStart, yEnd);
00825 sprintf(name, "HitMap2D_P%d_S%d_L%d_Box%d", i, j, k, boxID);
00826 m_h2DHitMap[i][j][k] = new TH2F(name, name, xBin, xStart, xEnd, yBin, yStart, yEnd);
00827 sprintf(name, "EffMap2D_P%d_S%d_L%d_Box%d", i, j, k, boxID);
00828 m_h2DEffMap[i][j][k] = new TH2F(name, name, xBin, xStart, xEnd, yBin, yStart, yEnd);
00829
00830
00831
00832 m_histArray->Add(m_h2DEffMap[i][j][k]);
00833 }
00834
00835 return StatusCode::SUCCESS;
00836 }
00837
00838
00839 StatusCode MucCalibMgr::InitClusterHisto()
00840 {
00841 char name[60];
00842 char title[60];
00843 int part, segment, layer, stripMax;
00844 part = segment = layer = stripMax = 0;
00845
00846 for( int i=0; i<LAYER_MAX; i++ )
00847 {
00848 sprintf( name, "LayerCluster_L%d", i );
00849 sprintf( title, "Clusters in L%d",i );
00850 m_hLayerCluster[i] = new TH1F(name,title, CLUSTER_RANGE, 0, CLUSTER_RANGE );
00851 }
00852
00853 for( int i=0; i<BOX_MAX; i++ )
00854 {
00855 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
00856 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
00857 sprintf( name, "BoxCluster_P%d_S%d_L%d_Box%d", part, segment, layer, i );
00858 sprintf( title, "Clusters in P%d_S%d_L%d Box%d", part, segment, layer, i );
00859 m_hBoxCluster[i] = new TH1F(name,title, CLUSTER_RANGE, 0, CLUSTER_RANGE );
00860 }
00861
00862 m_hLayerClusterCmp = new TH1F("LayerCluster", "LayerCluster", LAYER_MAX+1, -0.5, LAYER_MAX+0.5 );
00863 m_hBoxClusterCmp = new TH1F("BoxCluster", "BoxCluster", BOX_MAX+1, -0.5, BOX_MAX+0.5 );
00864
00865 return StatusCode::SUCCESS;
00866 }
00867
00868
00869 StatusCode MucCalibMgr::InitResHisto()
00870 {
00871 char name[60];
00872 char title[60];
00873
00874 for( int i=0; i<B_LAY_NUM; i++ )
00875 {
00876 sprintf( name, "BarrelRes_L%d", i );
00877 sprintf( title, "Barrel spacial resolution in L%d",i );
00878 m_hBarrelResDist[i] = new TH1F(name,title, 200, -100, 100 );
00879 }
00880
00881 for( int i=0; i<E_LAY_NUM; i++ )
00882 {
00883 sprintf( name, "EndcapRes_L%d", i );
00884 sprintf( title, "Endcap spacial resolution in L%d",i );
00885 m_hEndcapResDist[i] = new TH1F(name,title, 200, -100, 100 );
00886 }
00887
00888 m_hBarrelResComp[0] = new TH1F("BarrelResMean", "BarrelResMean", B_LAY_NUM+1, -0.5, B_LAY_NUM+0.5 );
00889 m_hBarrelResComp[1] = new TH1F("BarrelResSigma", "BarrelResSigma", B_LAY_NUM+1, -0.5, B_LAY_NUM+0.5 );
00890 m_hEndcapResComp[0] = new TH1F("EndcapResMean", "EndcapResMean", E_LAY_NUM+1, -0.5, E_LAY_NUM+0.5 );
00891 m_hEndcapResComp[1] = new TH1F("EndcapResSigma", "EndcapResSigma", E_LAY_NUM+1, -0.5, E_LAY_NUM+0.5 );
00892
00893 return StatusCode::SUCCESS;
00894 }
00895
00896
00897 StatusCode MucCalibMgr::InitConstTree()
00898 {
00899 MsgStream log(msgSvc, "MucCalibMgr");
00900 log << MSG::INFO << "Initialize Trees" << endreq;
00901
00902
00903 m_tJobLog = new TTree("JobLog","Job information");
00904 m_tJobLog->Branch("version", &m_vs, "m_vs/D");
00905 m_tJobLog->Branch("high_voltage", &m_hv, "m_hv/D");
00906 m_tJobLog->Branch("threshold", &m_th, "m_th/D");
00907 m_tJobLog->Branch("event_rate", &m_er, "m_er/D");
00908 m_tJobLog->Branch("input_tag", &m_tag, "m_tag/D");
00909 m_tJobLog->Branch("rec_mode", &m_recMode, "m_recMode/I");
00910 m_tJobLog->Branch("use_pad", &m_usePad, "m_usePad/I");
00911 m_tJobLog->Branch("eff_window", &m_effWindow, "m_effWindow/I");
00912 m_tJobLog->Branch("cluster_mode", &m_clusterMode, "m_clusterMode/I");
00913 m_tJobLog->Branch("check_event", &m_checkEvent, "m_checkEvent/I");
00914 m_tJobLog->Branch("dimu_select", &m_dimuSelect, "m_dimuSelect/I");
00915 m_tJobLog->Branch("dimu_only", &m_dimuOnly, "m_dimuOnly/I");
00916
00917 m_tJobLog->Branch("run_start", &m_fStartRun, "m_fStartRun/I");
00918 m_tJobLog->Branch("run_end", &m_fEndRun, "m_fEndRun/I");
00919 m_tJobLog->Branch("daq_time", &m_fTotalDAQTime, "m_fTotalDAQTime/D");
00920 m_tJobLog->Branch("job_time", &m_fTotalJobTime, "m_fTotalJobTime/D");
00921 m_tJobLog->Branch("calib_layer", &m_fCalibLayerNum, "m_fCalibLayerNum/D");
00922 m_tJobLog->Branch("calib_box", &m_fCalibBoxNum, "m_fCalibBoxNum/D");
00923 m_tJobLog->Branch("calib_strip", &m_fCalibStripNum, "m_fCalibStripNum/D");
00924 m_tJobLog->Branch("total_event", &m_fTotalEvent, "m_fTotalEvent/D");
00925 m_tJobLog->Branch("total_digi", &m_fTotalDigi, "m_fTotalDigi/D");
00926 m_tJobLog->Branch("total_exphit", &m_fTotalExpHit, "m_fTotalExpHit/D");
00927 m_tJobLog->Branch("total_effhit", &m_fTotalEffHit, "m_fTotalEffHit/D");
00928 m_tJobLog->Branch("total_noshit", &m_fTotalNosHit, "m_fTotalNosHit/D");
00929 m_tJobLog->Branch("total_cluster", &m_fTotalClstNum, "m_fTotalClstNum/D");
00930 m_tJobLog->Branch("total_strip_area",&m_fTotalStripArea,"m_fTotalStripArea/D");
00931 m_tJobLog->Branch("layer_coverage", &m_fLayerCoverage, "m_fLayerCoverage/D");
00932 m_tJobLog->Branch("box_coverage", &m_fBoxCoverage, "m_fBoxCoverage/D");
00933 m_tJobLog->Branch("strip_coverage", &m_fStripCoverage, "m_fStripCoverage/D");
00934
00935
00936 m_tStatLog = new TTree("StatLog","Statistic results");
00937
00938
00939 m_tLayConst = new TTree("LayConst","layer constants");
00940 m_tLayConst->Branch("layer_id", &m_fLayerId, "m_fLayerId/D");
00941 m_tLayConst->Branch("layer_eff", &m_fLayerEff, "m_fLayerEff/D");
00942 m_tLayConst->Branch("layer_efferr", &m_fLayerEffErr, "m_fLayerEffErr/D");
00943 m_tLayConst->Branch("layer_nosratio", &m_fLayerNosRatio,"m_fLayerNosRatio/D");
00944 m_tLayConst->Branch("layer_digi", &m_fLayerDigi, "m_fLayerDigi/D");
00945 m_tLayConst->Branch("layer_noise", &m_fLayerNos, "m_fLayerNos/D");
00946 m_tLayConst->Branch("layer_cnt", &m_fLayerCnt, "m_fLayerCnt/D");
00947 m_tLayConst->Branch("layer_exphit", &m_fLayerExpHit, "m_fLayerExpHit/D");
00948 m_tLayConst->Branch("layer_effHit", &m_fLayerEffHit, "m_fLayerEffHit/D");
00949 m_tLayConst->Branch("layer_nosHit", &m_fLayerNosHit, "m_fLayerNosHit/D");
00950 m_tLayConst->Branch("layer_cluster", &m_fLayerCluster, "m_fLayerCluster/D");
00951 m_tLayConst->Branch("layer_trknum", &m_fLayerTrkNum, "m_fLayerTrkNum/D");
00952
00953
00954 m_tBoxConst = new TTree("BoxConst","box constants");
00955 m_tBoxConst->Branch("box_id", &m_fBoxId, "m_fBoxId/D");
00956 m_tBoxConst->Branch("box_part", &m_fBoxPart, "m_fBoxPart/D");
00957 m_tBoxConst->Branch("box_segment", &m_fBoxSegment, "m_fBoxSegment/D");
00958 m_tBoxConst->Branch("box_layer", &m_fBoxLayer, "m_fBoxLayer/D");
00959 m_tBoxConst->Branch("box_eff", &m_fBoxEff, "m_fBoxEff/D");
00960 m_tBoxConst->Branch("box_efferr", &m_fBoxEffErr, "m_fBoxEffErr/D");
00961 m_tBoxConst->Branch("box_nosratio", &m_fBoxNosRatio,"m_fBoxNosRatio/D");
00962 m_tBoxConst->Branch("box_digi", &m_fBoxDigi, "m_fBoxDigi/D");
00963 m_tBoxConst->Branch("box_noise", &m_fBoxNos, "m_fBoxNos/D");
00964 m_tBoxConst->Branch("box_cnt", &m_fBoxCnt, "m_fBoxCnt/D");
00965 m_tBoxConst->Branch("box_exphit", &m_fBoxExpHit, "m_fBoxExpHit/D");
00966 m_tBoxConst->Branch("box_effhit", &m_fBoxEffHit, "m_fBoxEffHit/D");
00967 m_tBoxConst->Branch("box_noshit", &m_fBoxNosHit, "m_fBoxNosHit/D");
00968 m_tBoxConst->Branch("box_cluster", &m_fBoxCluster, "m_fBoxCluster/D");
00969 m_tBoxConst->Branch("box_trknum", &m_fBoxTrkNum, "m_fBoxTrkNum/D");
00970
00971
00972 m_tStrConst = new TTree("StrConst","strip constants");
00973 m_tStrConst->Branch("strip_id", &m_fStripId, "m_fStripId/D");
00974 m_tStrConst->Branch("strip_part", &m_fStripPart, "m_fStripPart/D");
00975 m_tStrConst->Branch("strip_segment", &m_fStripSegment, "m_fStripSegment/D");
00976 m_tStrConst->Branch("strip_layer", &m_fStripLayer, "m_fStripLayer/D");
00977 m_tStrConst->Branch("strip_eff", &m_fStripEff, "m_fStripEff/D");
00978 m_tStrConst->Branch("strip_efferr", &m_fStripEffErr, "m_fStripEffErr/D");
00979 m_tStrConst->Branch("strip_nosratio", &m_fStripNosRatio, "m_fStripNosRatio/D");
00980 m_tStrConst->Branch("strip_digi", &m_fStripDigi, "m_fStripDigi/D");
00981 m_tStrConst->Branch("strip_noise", &m_fStripNos, "m_fStripNos/D");
00982 m_tStrConst->Branch("strip_cnt", &m_fStripCnt, "m_fStripCnt/D");
00983 m_tStrConst->Branch("strip_effhit", &m_fStripEffHit, "m_fStripEffHit/D");
00984 m_tStrConst->Branch("strip_exphit", &m_fStripExpHit, "m_fStripExpHit/D");
00985 m_tStrConst->Branch("strip_noshit", &m_fStripNosHit, "m_fStripNosHit/D");
00986 m_tStrConst->Branch("strip_trknum", &m_fStripTrkNum, "m_fStripTrkNum/D");
00987
00988 return StatusCode::SUCCESS;
00989 }
00990
00991
00992 StatusCode MucCalibMgr::InitArea()
00993 {
00994 int stripMax, boxID, stripID;
00995 stripMax = boxID = stripID = 0;
00996 double totalArea = 0;
00997
00998 for( int i=0; i<PART_MAX; i++ )
00999 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
01000 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
01001 {
01002 MucBoxCal* aBox = new MucBoxCal( i, j, k );
01003 boxID = m_ptrIdTr->GetBoxId(i, j, k);
01004 m_boxResults[3][boxID] = aBox->GetArea();
01005 m_layerResults[3][k] += aBox->GetArea();
01006 delete aBox;
01007
01008 stripMax = m_ptrIdTr->GetStripMax( i, j, k );
01009 for( int m=0; m<stripMax; m++ )
01010 {
01011 MucStripCal* aStrip = new MucStripCal( i, j, k, m );
01012 stripID =m_ptrIdTr->GetStripId(i, j, k, m );
01013 m_stripResults[3][stripID] = aStrip->GetArea();
01014 totalArea += m_stripResults[3][stripID];
01015 delete aStrip;
01016 }
01017 }
01018
01019 for( int i=0; i<LAYER_MAX; i++ ) m_hLayerArea->Fill(i, m_layerResults[3][i]);
01020 for( int i=0; i<BOX_MAX; i++ ) m_hBoxArea->Fill(i, m_boxResults[3][i]);
01021 for( int i=0; i<STRIP_MAX; i++ ) m_hStripArea->Fill(i, m_stripResults[3][i]);
01022
01023 m_fTotalStripArea = totalArea;
01024
01025 return StatusCode::SUCCESS;
01026 }
01027
01028
01029
01030
01031
01032
01033 StatusCode MucCalibMgr::DimuSelect()
01034 {
01035 MsgStream log(msgSvc, "MucCalibMgr");
01036 log << MSG::INFO << "Dimu select" << endreq;
01037 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
01038
01039 if( m_tag > 0) { m_eventTag = (int)m_tag; return ( StatusCode::SUCCESS ); }
01040
01041 m_eventTag = 0;
01042 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
01043 if(!eventHeader) {
01044 log << MSG::FATAL << "Could not find event header" << endreq;
01045 return( StatusCode::FAILURE );
01046 }
01047
01048
01049 SmartDataPtr<RecMdcTrackCol> mdcTrackCol(eventSvc,"/Event/Recon/RecMdcTrackCol");
01050 if(!mdcTrackCol) {
01051 log << MSG::FATAL << "Could not find mdc tracks" << endreq;
01052 return( StatusCode::FAILURE);
01053 }
01054 log << MSG::INFO << mdcTrackCol->size() << endreq;
01055 if( mdcTrackCol->size() != 2 ) return ( StatusCode::FAILURE );
01056
01057 log << MSG::INFO << "r1:\t" << (*mdcTrackCol)[0]->r() << "\tz1:" << (*mdcTrackCol)[0]->z() << endreq;
01058 log << MSG::INFO << "r2:\t" << (*mdcTrackCol)[1]->r() << "\tz2:" << (*mdcTrackCol)[1]->z() << endreq;
01059 log << MSG::INFO << "p1:\t" << (*mdcTrackCol)[0]->p() << "\tp2:" << (*mdcTrackCol)[1]->p() << endreq;
01060
01061 bool mdcFlag1 = (*mdcTrackCol)[0]->charge() + (*mdcTrackCol)[1]->charge() == 0;
01062 bool mdcFlag2 = fabs((*mdcTrackCol)[0]->r())<=1 && fabs((*mdcTrackCol)[0]->z())<3;
01063 bool mdcFlag3 = fabs((*mdcTrackCol)[1]->r())<=1 && fabs((*mdcTrackCol)[1]->z())<3;
01064 bool mdcFlag4 = (*mdcTrackCol)[0]->p()<2 && (*mdcTrackCol)[1]->p()<2;
01065 bool mdcFlag5 = fabs( (*mdcTrackCol)[0]->p() - (*mdcTrackCol)[1]->p() )/( (*mdcTrackCol)[0]->p() + (*mdcTrackCol)[1]->p() ) < 0.5;
01066 bool mdcPass = false;
01067 if( mdcFlag1 && mdcFlag2 && mdcFlag3 && mdcFlag4 && mdcFlag5 ) mdcPass = true;
01068
01069
01070 SmartDataPtr<RecTofTrackCol> tofTrackCol(eventSvc,"/Event/Recon/RecTofTrackCol");
01071 if(!tofTrackCol) {
01072 log << MSG::FATAL << "Could not find RecTofTrackCol!" << endreq;
01073 return( StatusCode::FAILURE);
01074 }
01075 log << MSG::INFO << "TOF tracks:\t" << tofTrackCol->size() << endreq;
01076
01077 double t1, t2;
01078 t1 = 0., t2 = 1000;
01079
01080 if(tofTrackCol->size() > 7 && tofTrackCol->size() < 21)
01081 {
01082 int goodtof = 0;
01083 for(unsigned int itof = 0; itof < tofTrackCol->size(); itof++)
01084 {
01085 TofHitStatus *status = new TofHitStatus;
01086 status->setStatus((*tofTrackCol)[itof]->status());
01087
01088 if( !(status->is_cluster()) ) { delete status; continue; }
01089 if(goodtof==0) t1 = (*tofTrackCol)[itof]->tof();
01090 if(goodtof==1) t2 = (*tofTrackCol)[itof]->tof();
01091
01092 goodtof++;
01093 delete status;
01094 }
01095 }
01096 bool tofPass = false;
01097 if( fabs( t1-t2 ) < 4) tofPass = true;
01098
01099
01100 SmartDataPtr<RecEmcShowerCol> emcShowerCol(eventSvc,"/Event/Recon/RecEmcShowerCol");
01101 if (!emcShowerCol) {
01102 log << MSG::FATAL << "Could not find RecEmcShowerCol!" << endreq;
01103 return( StatusCode::FAILURE);
01104 }
01105 log << MSG::INFO << "EMC showers:\t" << emcShowerCol->size() << endreq;
01106
01107 if( emcShowerCol->size() < 2 ) return ( StatusCode::SUCCESS );
01108
01109 double e1, e2, theta1, theta2, phi1, phi2;
01110 int part;
01111
01112 e1 = (*emcShowerCol)[0]->energy(); e2 = (*emcShowerCol)[1]->energy();
01113 theta1 = (*emcShowerCol)[0]->theta(); theta2 = (*emcShowerCol)[1]->theta();
01114 phi1 = (*emcShowerCol)[0]->phi(); phi2 = (*emcShowerCol)[1]->phi();
01115 part = (*emcShowerCol)[0]->module();
01116
01117 log << MSG::INFO << "e1:\t" << e1 << "\te2:\t" << e2 << endreq;
01118 log << MSG::INFO << "theta1:\t" << theta1 << "\ttheta2:\t" << theta2 << endreq;
01119 log << MSG::INFO << "phi1:\t" << phi1 << "\tphi2:\t" << phi2 << endreq;
01120
01121 bool emcFlag1 = e1 < 1.0 && e1 > 0.1;
01122 bool emcFlag2 = e2 < 1.0 && e2 > 0.1;
01123 bool emcFlag3 = fabs(theta1 + theta2 - PI) < 0.05;
01124 bool emcFlag4 = fabs(phi1 - phi2) - PI < 0.4;
01125
01126 bool emcPass = false;
01127 if( emcFlag1 && emcFlag2 && emcFlag3 && emcFlag4 ) emcPass = true;
01128
01129
01130 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc,"/Event/Digi/MucDigiCol");
01131 if(!mucDigiCol) {
01132 log << MSG::FATAL << "Could not find MUC digi" << endreq;
01133 return( StatusCode::FAILURE);
01134 }
01135 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc,"/Event/Recon/RecMucTrackCol");
01136 if (!mucTrackCol) {
01137 log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
01138 return( StatusCode::FAILURE);
01139 }
01140 log << MSG::INFO << "digi:\t" << mucDigiCol->size() << "tracks:\t" << mucTrackCol->size() << endreq;
01141
01142 bool mucFlag1 = mucDigiCol->size()>6;
01143 bool mucFlag2 = mucTrackCol->size()>1;
01144 bool mucPass = false;
01145 if( mucFlag1 && mucFlag2 ) mucPass = true;
01146
01147 if( mdcPass && tofPass && emcPass && mucPass ) m_eventTag = 1;
01148 else m_eventTag = (int)m_tag;
01149
01150 return( StatusCode::SUCCESS );
01151 }
01152
01153
01154 StatusCode MucCalibMgr::ReadEvent()
01155 {
01156 MsgStream log(msgSvc, "MucCalibMgr");
01157 log << MSG::INFO << "Read event" << endreq;
01158
01159 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
01160 m_evtBegin = clock();
01161
01162
01163 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
01164 if(!eventHeader) {
01165 log << MSG::FATAL << "Could not find event header" << endreq;
01166 return( StatusCode::FAILURE );
01167 }
01168
01169 m_currentRun = eventHeader->runNumber();
01170 m_currentEvent = eventHeader->eventNumber();
01171 if( m_fStartRun == 0 ) m_fStartRun = m_currentRun;
01172 m_fEndRun = m_currentRun;
01173 m_fTotalEvent++;
01174
01175 log << MSG::INFO << "Run [ " << m_currentRun << " ]\tEvent [ " << m_currentEvent << " ]" << endreq;
01176 if( ((long)m_fTotalEvent)%2000 == 0 ) cout << m_fTotalEvent << "\tdone!" << endl;
01177
01178
01179 if( m_dimuSelect ) {
01180 MucCalibMgr::DimuSelect();
01181 log << MSG::INFO << "Event tag:\t" << m_eventTag << endreq;
01182 if( m_dimuOnly && m_eventTag != 1 ) return( StatusCode::FAILURE );
01183 }
01184
01185
01186 log << MSG::INFO << "Retrieve digis" << endreq;
01187
01188 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc,"/Event/Digi/MucDigiCol");
01189 if(!mucDigiCol) {
01190 log << MSG::FATAL << "Could not find MUC digi" << endreq;
01191 return( StatusCode::FAILURE);
01192 }
01193
01194 int part, segment, layer, strip, pad;
01195 part = segment = layer = strip = pad = 0;
01196 double padX, padY, padZ;
01197 padX = padY = padZ = 0.;
01198 double resMax = 0.;
01199
01200 Identifier mucId;
01201 MucDigiCol::iterator digiIter = mucDigiCol->begin();
01202 int eventDigi = 0;
01203 for ( int digiId =0; digiIter != mucDigiCol->end(); digiIter++, digiId++ )
01204 {
01205 mucId = (*digiIter)->identify();
01206 part = MucID::barrel_ec(mucId);
01207 segment = MucID::segment(mucId);
01208 layer = MucID::layer(mucId);
01209 strip = MucID::channel(mucId);
01210
01211 log << MSG::DEBUG << "[" << part << "\t" << segment << "\t" << layer << "\t" << strip << "]\t" ;
01212 if( (digiId+1)%8 == 0 ) log << MSG::DEBUG << endreq;
01213
01214 eventDigi ++;
01215
01216 if(abs(part)>=PART_MAX || abs(segment)>=SEGMENT_MAX || abs(layer)>=LAYER_MAX || abs(strip)>=STRIP_INBOX_MAX) {
01217 log << MSG::ERROR << endreq << "Digi IDs slop over!" << endreq;
01218 continue;
01219 }
01220
01221
01222 MucMark *aMark = new MucMark( part, segment, layer, strip );
01223 m_digiCol.push_back( aMark );
01224 m_segDigiCol[part][segment].push_back( aMark );
01225 }
01226 log << MSG::DEBUG << endreq;
01227 log << MSG::INFO << "Total digits of this event: " << eventDigi << endreq;
01228 if( eventDigi > 200) log << MSG::ERROR << "Event: " << m_currentEvent << "\tdigits sharply rise:\t" << eventDigi << endreq;
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246 int clusterNum, bigClusterNum, clusterSize;
01247 clusterNum = bigClusterNum = clusterSize = 0;
01248 if( m_clusterMode ) {
01249 log << MSG::INFO << "Searching clusters" << endreq;
01250 m_clusterCol = (*m_ptrMucMark).CreateClusterCol(m_clusterMode, m_digiCol );
01251 }
01252
01253 for( unsigned int i=0; i<m_clusterCol.size(); i++ )
01254 {
01255 clusterSize = m_clusterCol[i].size();
01256
01257 if( clusterSize > CLUSTER_ALARM )
01258 {
01259 log << MSG::WARNING << "Big cluster:" << endreq;
01260 part = (*m_clusterCol[i][0]).Part();
01261 segment = (*m_clusterCol[i][0]).Segment();
01262 layer = (*m_clusterCol[i][0]).Layer();
01263
01264 if( m_clusterSave ) (*m_fdata) << "Event:\t" << m_currentEvent << "\tbig cluster " << bigClusterNum << endl;
01265
01266 for( int j=0; j<clusterSize; j++ )
01267 {
01268 strip = (*m_clusterCol[i][j]).Strip();
01269 log << MSG::WARNING << "[" << part << "\t" << segment << "\t" << layer << "\t" << strip << "]\t";
01270 if( (j+1)%8 == 0 ) log << MSG::WARNING << endreq;
01271 if( m_clusterSave ) (*m_fdata) << part << "\t" << segment << "\t" << layer << "\t" << strip << endl;
01272 }
01273 log << MSG::WARNING << endreq;
01274 bigClusterNum ++;
01275 }
01276 else if( clusterSize > 1 )
01277 {
01278 log << MSG::DEBUG << "cluster: " << clusterNum << endreq;
01279 clusterNum ++, m_fTotalClstNum ++;
01280 part = (*m_clusterCol[i][0]).Part();
01281 segment = (*m_clusterCol[i][0]).Segment();
01282 layer = (*m_clusterCol[i][0]).Layer();
01283 for( int j=0; j<clusterSize; j++ )
01284 {
01285 strip = (*m_clusterCol[i][j]).Strip();
01286 log << MSG::DEBUG << "[" << part << "\t" << segment << "\t" << layer << "\t" << strip << "]\t";
01287 if( (j+1)%8 == 0 ) log << MSG::DEBUG << endreq;
01288 }
01289 log << MSG::DEBUG << endreq;
01290 }
01291 }
01292
01293 if( m_clusterMode) log << MSG::INFO << "Total clusters in this event: " << clusterNum << endreq;
01294 else log << MSG::INFO << "Clusters not built" << endreq;
01295
01296
01297
01298 log << MSG::INFO << "Retrieve tracks" << endreq;
01299
01300 SmartDataPtr<RecMdcTrackCol> mdcTrackCol(eventSvc,"/Event/Recon/RecMdcTrackCol");
01301 if(!mdcTrackCol) {
01302 log << MSG::FATAL << "Could not find mdc tracks" << endreq;
01303 return( StatusCode::FAILURE);
01304 }
01305
01306 RecMdcTrackCol::iterator mdctrkIter = mdcTrackCol->begin();
01307 for (; mdctrkIter != mdcTrackCol->end(); mdctrkIter++)
01308 {
01309 m_charge = (*mdctrkIter)->charge();
01310 m_mdcpx = (*mdctrkIter)->px();
01311 m_mdcpy = (*mdctrkIter)->py();
01312 m_mdcpz = (*mdctrkIter)->pz();
01313 m_mdcpt = (*mdctrkIter)->pxy();
01314 m_mdcpp = (*mdctrkIter)->p();
01315 m_mdcphi = (*mdctrkIter)->phi();
01316 m_mdctheta = (*mdctrkIter)->theta();
01317 m_mdcTrkInfoTuple->write();
01318 }
01319
01320
01321 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc,"/Event/Recon/RecMucTrackCol");
01322 if (!mucTrackCol) {
01323 log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
01324 return( StatusCode::FAILURE);
01325 }
01326
01327 RecMucTrackCol *aRecMucTrackCol = mucTrackCol;
01328 if (aRecMucTrackCol->size() < 1) {
01329 log << MSG::INFO << "No MUC tracks in this event" << endreq;
01330 return StatusCode::SUCCESS;
01331 }
01332 log << MSG::INFO << "Total tracks of this event: " << aRecMucTrackCol->size() << endreq;
01333
01334
01335
01336
01337 if( 0 )
01338 {
01339 SmartDataPtr<RecEsTimeCol> aRecEsTimeCol(eventSvc,"/Event/Recon/RecEsTimeCol");
01340 if( ! aRecEsTimeCol ){
01341 log << MSG::ERROR << "Could not find RecEsTimeCol" << endreq;
01342 return StatusCode::FAILURE;
01343 }else{
01344 RecEsTimeCol::iterator iter_evt = aRecEsTimeCol->begin();
01345
01346
01347 m_ntEsTime = (*iter_evt)->getStat();
01348 if( (*iter_evt)->getStat() != 211 ) {
01349 log << MSG::WARNING << "Event time not by TOF, skip!" << endreq;
01350 return StatusCode::SUCCESS;
01351 }
01352 }
01353 }
01354
01355
01356 double phi1, phi2, phiDiff, theta1, theta2, thetaDiff;
01357 phiDiff = thetaDiff = 0.;
01358 if( aRecMucTrackCol->size()==2 && (*aRecMucTrackCol)[0]->GetTotalHits() > 4 && (*aRecMucTrackCol)[1]->GetTotalHits() > 4 )
01359 {
01360
01361 phi1 = (*aRecMucTrackCol)[0]->getMucPos().phi(); phi2 = (*aRecMucTrackCol)[1]->getMucPos().phi();
01362 if( phi1 > 0 ) phiDiff = phi1 - phi2 - PI;
01363 else phiDiff = phi2 - phi1 - PI;
01364
01365
01366 theta1 = (*aRecMucTrackCol)[0]->getMucPos().theta(); theta2 = (*aRecMucTrackCol)[1]->getMucPos().theta();
01367 thetaDiff = theta1 + theta2 - PI;
01368 m_hTrackPosPhiDiff->Fill( phiDiff );
01369 m_hTrackPosThetaDiff->Fill( thetaDiff );
01370 m_hDimuTracksPosDiff->Fill( thetaDiff, phiDiff );
01371 m_ntPosPhiDiff = phiDiff;
01372 m_ntPosThetaDiff = thetaDiff;
01373
01374 log << MSG::INFO << "PosPhiDiff:\t" << phiDiff << "\tPosThetaDiff:\t" << thetaDiff << endreq;
01375
01376
01377 phi1 = (*aRecMucTrackCol)[0]->getMucMomentum().phi(); phi2 = (*aRecMucTrackCol)[1]->getMucMomentum().phi();
01378 if( phi1 > 0 ) phiDiff = phi1 - phi2 - PI;
01379 else phiDiff = phi2 - phi1 - PI;
01380
01381
01382 theta1 = (*aRecMucTrackCol)[0]->getMucMomentum().theta(); theta2 = (*aRecMucTrackCol)[1]->getMucMomentum().theta();
01383 thetaDiff = theta1 + theta2 - PI;
01384
01385 m_hTrackMomPhiDiff->Fill( phiDiff );
01386 m_hTrackMomThetaDiff->Fill( thetaDiff );
01387 m_hDimuTracksMomDiff->Fill( thetaDiff, phiDiff );
01388 m_ntMomPhiDiff = phiDiff;
01389 m_ntMomThetaDiff = thetaDiff;
01390
01391 log << MSG::INFO << "MomPhiDiff:\t" << phiDiff << "\tMomThetaDiff:\t" << thetaDiff << endreq;
01392 m_ntDimuTag = m_eventTag;
01393 m_trackDiffTuple->write();
01394 }
01395
01396
01397 RecMucTrackCol::iterator trackIter = mucTrackCol->begin();
01398 int trackHitNum, rawHitNum, expectedHitNum, segNum, trkRecMode, lastLayerBR, lastLayerEC;
01399 int layerPassNum[3], passMax[TRACK_SEG_MAX][2];
01400 bool firedLay[TRACK_SEG_MAX][LAYER_MAX];
01401 bool seedList[PART_MAX][LAYER_MAX];
01402 trackHitNum = rawHitNum = expectedHitNum = segNum = trkRecMode = lastLayerBR = lastLayerEC = 0;
01403 layerPassNum[0] = layerPassNum[1] = layerPassNum[2] = 0;
01404 for( int segi=0; segi<TRACK_SEG_MAX; segi++ ) {
01405 passMax[segi][0] = passMax[segi][1] = 0;
01406 for( int layi=0; layi<LAYER_MAX; layi++ ) firedLay[segi][layi] = 0;
01407 }
01408
01409 mark_col trkSeg[TRACK_SEG_MAX];
01410 vector<MucRecHit*> mucRawHitCol;
01411 vector<MucRecHit*> mucExpHitCol;
01412
01413 for (int trackId = 0; trackIter != mucTrackCol->end(); trackIter++, trackId++)
01414 {
01415 trackHitNum = (*trackIter)->GetTotalHits();
01416 log << MSG::DEBUG << "Track: " << trackId << " Hits: " << trackHitNum << endreq;
01417
01418 if( trackHitNum == 0 ) {
01419 log << MSG::INFO << "Track " << trackId << " no hits" << endreq;
01420 continue;
01421 }
01422
01423 m_ntTrackHits = trackHitNum;
01424
01425 m_trkRecMode = trkRecMode = (*trackIter)->GetRecMode();
01426 m_chi2 = (*trackIter)->chi2();
01427 m_px = (*trackIter)->getMucMomentum().x();
01428 m_py = (*trackIter)->getMucMomentum().y();
01429 m_pz = (*trackIter)->getMucMomentum().z();
01430 m_pt = sqrt(m_px*m_px + m_py*m_py);
01431 m_pp = sqrt(m_px*m_px + m_py*m_py + m_pz*m_pz);
01432
01433
01434 m_r = (*trackIter)->getMucPos().mag();
01435 m_cosTheta = (*trackIter)->getMucPos().cosTheta();
01436 m_theta = (*trackIter)->getMucPos().theta();
01437 m_phi = (*trackIter)->getMucPos().phi();
01438 m_depth = (*trackIter)->depth();
01439 m_brLastLayer = lastLayerBR = (*trackIter)->brLastLayer();
01440 m_ecLastLayer = lastLayerEC = (*trackIter)->ecLastLayer();
01441 m_totalHits = (*trackIter)->numHits();
01442 m_totalLayers = (*trackIter)->numLayers();
01443 m_maxHitsInLayer = (*trackIter)->maxHitsInLayer();
01444
01445
01446 m_hPhiCosTheta->Fill(m_cosTheta, m_phi);
01447 log << MSG::INFO << "Fill track info" << endreq;
01448
01449 MucRecHit* pMucRawHit;
01450 MucRecHit* pMucExpHit;
01451 if( m_calHitCol.size() != 0 ) m_calHitCol.clear();
01452
01453
01454 log << MSG::DEBUG << "Reconstruction hits(digis in a track): " << endreq;
01455 mucRawHitCol = (*trackIter)->GetHits();
01456 rawHitNum += mucRawHitCol.size();
01457
01458 segNum = 0;
01459 if( trkRecMode == 3 ) {
01460 for(int iPart=0; iPart<PART_MAX; iPart++)
01461 for(int iLayer=0; iLayer<LAYER_MAX; iLayer++) seedList[iPart][iLayer] = false;
01462 }
01463
01464 for(unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
01465 {
01466 pMucRawHit = mucRawHitCol[ hitId ];
01467 part = pMucRawHit->Part();
01468 segment = pMucRawHit->Seg();
01469 layer = pMucRawHit->Gap();
01470 strip = pMucRawHit->Strip();
01471
01472 log << MSG::DEBUG << "[" << part << "\t" << segment << "\t" << layer << "\t" << strip << "]\t";
01473
01474
01475
01476 MucMark *aMark = new MucMark( part, segment, layer, strip );
01477 m_calHitCol.push_back( aMark );
01478
01479
01480 if( trkRecMode == 3 ) seedList[part][layer] = pMucRawHit->HitIsSeed();
01481
01482
01483 if(hitId == 0) { trkSeg[segNum].push_back( aMark ); segNum ++; }
01484 else
01485 {
01486 log << MSG::DEBUG << "segNum: " << segNum << endreq;
01487 bool notInSeg = true;
01488 for( int segi=0; segi<segNum; segi++ )
01489 {
01490 if( aMark->IsInSegWith( *(trkSeg[segi][0]) ) )
01491 {
01492 trkSeg[segi].push_back( aMark );
01493 notInSeg = false;
01494 break;
01495 }
01496 }
01497
01498 if( notInSeg == true )
01499 {
01500 trkSeg[segNum].push_back( aMark );
01501 segNum ++;
01502 if( segNum > TRACK_SEG_MAX ) {
01503 log << MSG::ERROR << "Track segment overflow: " << segNum << endreq;
01504 break;
01505 }
01506 }
01507 }
01508 }
01509 log << MSG::DEBUG << endreq;
01510
01511
01512 layerPassNum[0] = layerPassNum[1] = layerPassNum[2] = 0;
01513 for( int segi=0; segi<segNum; segi++ )
01514 {
01515 int tmpLayNum = 0;
01516 passMax[segi][0] = passMax[segi][1] = trkSeg[segi][0]->Layer();
01517 for( unsigned int hiti=1; hiti<trkSeg[segi].size(); hiti++ )
01518 {
01519 if( trkSeg[segi][hiti]->Layer() < passMax[segi][0] )
01520 passMax[segi][0] = trkSeg[segi][hiti]->Layer();
01521 if( trkSeg[segi][hiti]->Layer() > passMax[segi][1] )
01522 passMax[segi][1] = trkSeg[segi][hiti]->Layer();
01523 firedLay[segi][trkSeg[segi][hiti]->Layer()] = 1;
01524 }
01525
01526 for( int layi=0; layi<LAYER_MAX; layi++ ) {
01527 if( firedLay[segi][layi] ) tmpLayNum ++;
01528 }
01529
01530 if( segi == 0 ) layerPassNum[0] += passMax[segi][1] + 1;
01531 else layerPassNum[0] += (passMax[segi][1] - passMax[segi][0] + 1);
01532
01533 layerPassNum[1] += (passMax[segi][1] - passMax[segi][0] + 1);
01534 layerPassNum[2] += tmpLayNum;
01535
01536 trkSeg[segi].clear();
01537 }
01538 m_ntTrackEvent = m_currentEvent;
01539 m_ntTrackTag = m_eventTag;
01540 m_ntTrackSegFly = segNum;
01541 m_ntTrackLayFlyA = layerPassNum[0];
01542 m_ntTrackLayFlyB = layerPassNum[1];
01543 m_ntTrackLayFlyC = layerPassNum[2];
01544 m_trackInfoTuple->write();
01545 log << MSG::INFO << "Track\t" << trackId << "\tsegment(s):\t" << segNum
01546 << "\tlayer passed:\t" << layerPassNum[0] <<"\t" << layerPassNum[1] << "\t" << layerPassNum[2] << endreq;
01547
01548
01549
01550
01551 log << MSG::DEBUG << "Fitting hits(expected hits in a track): " << endreq;
01552 mucExpHitCol = (*trackIter)->GetExpectedHits();
01553 expectedHitNum += mucExpHitCol.size();
01554 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
01555 {
01556 pMucRawHit = mucExpHitCol[ hitId ];
01557 part = pMucRawHit->Part(); segment = pMucRawHit->Seg();
01558 layer = pMucRawHit->Gap(); strip = pMucRawHit->Strip();
01559
01560 if( m_usePad != 0 )
01561 {
01562 pad = pMucRawHit->GetPadID(); padZ = pMucRawHit->GetIntersectZ();
01563 padX = pMucRawHit->GetIntersectY(); padY = pMucRawHit->GetIntersectX();
01564
01565 if( part != BRID )
01566 {
01567 if(segment == 1) { padX = -padX; }
01568 else if(segment == 2) { padX = -padX, padY = -padY; }
01569 else if(segment == 3) { padY = -padY; }
01570 }
01571 }
01572
01573
01574
01575
01576 MucMark* currentMark = new MucMark( part, segment, layer, strip );
01577 m_expHitCol.push_back( currentMark );
01578
01579
01580
01581
01582 int isInPos = -1;
01583 bool isInEffWindow = false;
01584 isInPos = currentMark->IsInCol( m_segDigiCol[part][segment] );
01585
01586
01587 if( part == BRID && (layer-lastLayerBR>1) ) continue;
01588 if( part != BRID && (layer-lastLayerEC>1) ) continue;
01589
01590
01591 if( part==BRID && layer==0 && (strip<2 || strip>45) )
01592 {
01593 if( isInPos != -1)
01594 {
01595 m_record[part][segment][layer][strip][2] ++;
01596 m_record[part][segment][layer][strip][1] ++;
01597 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
01598
01599 if( m_usePad != 0 ) {
01600 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
01601 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
01602 }
01603
01604 }
01605 else {
01606 m_record[part][segment][layer][strip][1] ++;
01607 if( m_usePad != 0 ) m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
01608 }
01609 continue;
01610 }
01611
01612
01613 if( isInPos != -1 )
01614 {
01615 m_record[part][segment][layer][strip][2] ++;
01616 m_record[part][segment][layer][strip][1] ++;
01617 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
01618
01619 if( m_usePad != 0 ) {
01620 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
01621 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
01622 }
01623
01624 continue;
01625 }
01626 else for(int tempStrip=0, hiti=-m_effWindow; hiti<=m_effWindow; hiti++ )
01627 {
01628 if( hiti == 0 ) continue;
01629 tempStrip = strip + hiti;
01630 if( tempStrip < 0 || tempStrip > m_ptrIdTr->GetStripMax(part,segment,layer) ) continue;
01631
01632 isInPos = m_ptrMucMark->IsInCol( part, segment, layer, tempStrip, m_segDigiCol[part][segment] );
01633 if( isInPos != -1 )
01634 {
01635 m_record[part][segment][layer][tempStrip][2] ++;
01636 m_record[part][segment][layer][tempStrip][1] ++;
01637 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
01638
01639 if( m_usePad != 0 ) {
01640 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
01641 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
01642 }
01643
01644 m_ntEffWindow = hiti;
01645 m_effWindowTuple->write();
01646 isInEffWindow = true;
01647 }
01648
01649 }
01650
01651 if( isInEffWindow ) { continue; }
01652 else {
01653 m_record[part][segment][layer][strip][1] ++;
01654 if( m_usePad != 0 ) m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
01655 }
01656
01657 }
01658
01659
01660 log << MSG::INFO << "Fill residual" << endreq;
01661 vector<float> m_lineResCol = (*trackIter)->getDistHits();
01662 vector<float> m_quadResCol = (*trackIter)->getQuadDistHits();
01663 vector<float> m_extrResCol = (*trackIter)->getExtDistHits();
01664 int mode = (*trackIter)->GetRecMode();
01665
01666 for(unsigned int nres = 0; nres < m_lineResCol.size(); nres++ )
01667 if( fabs(m_lineResCol[nres])>resMax ) resMax = fabs(m_lineResCol[nres]);
01668
01669 log << MSG::INFO << "Good track for res" << endreq;
01670 if( trackHitNum > 4 && m_lineResCol[0] != -99)
01671 {
01672
01673 bool firedFlag[PART_MAX][LAYER_MAX][2];
01674 for(int iprt=0; iprt<PART_MAX; iprt++)
01675 for(int jlay=0; jlay<LAYER_MAX; jlay++)
01676 firedFlag[iprt][jlay][0] = firedFlag[iprt][jlay][1] = false;
01677
01678 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
01679 {
01680 pMucExpHit = mucExpHitCol[ hitId ];
01681 part = pMucExpHit->Part(); segment = pMucExpHit->Seg(); layer = pMucExpHit->Gap();
01682 firedFlag[part][layer][0] = true;
01683 }
01684
01685 log << MSG::INFO << "Fit res" << endreq;
01686 for(unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
01687 {
01688 pMucRawHit = mucRawHitCol[ hitId ];
01689 part = pMucRawHit->Part(); segment = pMucRawHit->Seg(); layer = pMucRawHit->Gap();
01690
01691 if( part == BRID ) m_hBarrelResDist[layer]->Fill( m_lineResCol[hitId] );
01692 else m_hEndcapResDist[layer]->Fill( m_lineResCol[hitId] );
01693
01694
01695 if( firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false )
01696 {
01697 m_resPart = part;
01698 m_resSegment = segment;
01699 m_resLayer = layer;
01700 m_lineRes = m_lineResCol[hitId];
01701 m_quadRes = m_quadResCol[hitId];
01702 m_extrRes = m_extrResCol[hitId];
01703 m_resFired = 1;
01704 m_resMode = mode;
01705 m_resInfoTuple->write();
01706 }
01707
01708 firedFlag[part][layer][1] = true;
01709 }
01710
01711 log << MSG::INFO << "Exp res" << endreq;
01712 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
01713 {
01714 pMucExpHit = mucExpHitCol[ hitId ];
01715 part = pMucExpHit->Part(); segment = pMucExpHit->Seg(); layer = pMucExpHit->Gap();
01716
01717 if(firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false)
01718 {
01719 m_resPart = part;
01720 m_resSegment = segment;
01721 m_resLayer = layer;
01722 m_lineRes = 1000;
01723 m_quadRes = 1000;
01724 m_extrRes = 1000;
01725 m_resFired = 0;
01726 m_resMode = mode;
01727 m_resInfoTuple->write();
01728 }
01729 }
01730
01731 }
01732
01733 mucRawHitCol.clear();
01734 mucExpHitCol.clear();
01735
01736 }
01737
01738 if( resMax > 300 ) cout <<"Res too big!\t"<< m_fTotalEvent <<"\t"<< m_currentRun <<"\t"<< m_currentEvent <<"\t"<< resMax << endl;
01739
01740 m_ntTrackNum = mucTrackCol->size();
01741
01742 m_fTotalEffHit += rawHitNum;
01743 log << MSG::INFO << "Total hits in this event, raw: " << rawHitNum << "\texpected: " << expectedHitNum << endreq;
01744
01745
01746
01747 log << MSG::INFO << "Searching inc/noise hits" << endreq;
01748 bool isNosHit;
01749 bool hasEffHit;
01750 for( unsigned int i=0; i < m_digiCol.size(); i++ )
01751 {
01752 isNosHit = true;
01753
01754 if( m_digiCol[i]->IsInCol( m_effHitCol ) !=-1) continue;
01755 else
01756 {
01757 for( unsigned int j=0; j < m_clusterCol.size(); j++ )
01758 {
01759 hasEffHit = false;
01760 for( unsigned int k=0; k<m_clusterCol[j].size(); k++)
01761 {
01762 if( m_clusterCol[j][k]->IsInCol(m_effHitCol) != -1)
01763 {
01764 hasEffHit = true;
01765 break;
01766 }
01767 }
01768
01769 if ( hasEffHit && (m_digiCol[i]->IsInCol( m_clusterCol[j] ) != -1) ) {
01770 isNosHit = false;
01771 break;
01772 }
01773 }
01774
01775 if( isNosHit ) {
01776 m_nosHitCol.push_back( m_digiCol[i] );
01777 m_fTotalNosHit ++;
01778 }
01779 }
01780 }
01781
01782 return StatusCode::SUCCESS;
01783 }
01784
01785
01786 StatusCode MucCalibMgr::CheckEvent()
01787 {
01788 MsgStream log(msgSvc, "MucCalibMgr");
01789 log << MSG::INFO << "Check event" << endreq;
01790
01791
01792
01793 log << MSG::INFO << "Searching over digi" << endreq;
01794 int overDigi = 0;
01795 for( unsigned int i=0; i < m_digiCol.size(); i++ )
01796 for( unsigned int j=i+1; j < m_digiCol.size(); j++ ) {
01797 if( (*m_digiCol[i]) == (*m_digiCol[j]) ) overDigi ++;
01798 }
01799
01800 if( overDigi !=0 )
01801 log << MSG::ERROR << overDigi << " over digi found in digis!" << endreq;
01802
01803
01804
01805 log << MSG::INFO << "Searching over hit" << endreq;
01806 int overHit = 0;
01807 for( unsigned int i=0; i < m_expHitCol.size(); i++ )
01808 for( unsigned int j=i+1; j < m_expHitCol.size(); j++ ) {
01809 if( (*m_expHitCol[i]) == (*m_expHitCol[j]) ) overHit ++;
01810 }
01811
01812 if( overHit != 0 )
01813 log << MSG::WARNING << overHit << " over hits found in rec hits!" << endreq;
01814
01815
01816 log << MSG::INFO << "Searching new hit" << endreq;
01817 int newHit = 0;
01818 int num = 0;
01819 for( unsigned int i=0; i < m_expHitCol.size(); i++ )
01820 {
01821 num = m_expHitCol[i]->NumInCol( m_digiCol );
01822 if( num == 0 )
01823 {
01824 log << MSG::ERROR << "Exp hit not in digis!"
01825 << "prt: " << (*m_expHitCol[i]).Part()
01826 << "\tseg: " << (*m_expHitCol[i]).Segment()
01827 << "\tlay: " << (*m_expHitCol[i]).Layer()
01828 << "\tstr: " << (*m_expHitCol[i]).Strip() << endreq;
01829
01830 newHit ++;
01831 }
01832 }
01833
01834 if( newHit != 0 )
01835 log << MSG::WARNING << newHit << " new hit(s) found in rec hits!" << endreq;
01836
01837 return StatusCode::SUCCESS;
01838 }
01839
01840
01841 StatusCode MucCalibMgr::FillEvent()
01842 {
01843 MsgStream log(msgSvc, "MucCalibMgr");
01844 log << MSG::INFO << "Fill event" << endreq;
01845
01846 int part, segment, layer, strip, size;
01847 part = segment = layer = strip = size = 0;
01848
01849
01850 log << MSG::INFO << "Fill digis" << endreq;
01851 for( unsigned int i=0; i<m_digiCol.size(); i++ )
01852 {
01853 part = (*m_digiCol[i]).Part();
01854 segment = (*m_digiCol[i]).Segment();
01855 layer = (*m_digiCol[i]).Layer();
01856 strip = (*m_digiCol[i]).Strip();
01857
01858 FillDigi( part, segment, layer, strip );
01859 m_record[part][segment][layer][strip][0] ++;
01860 m_fTotalDigi ++;
01861 }
01862
01863
01864 log << MSG::INFO << "Fill rec hits" << endreq;
01865 for( unsigned int i=0; i<m_expHitCol.size(); i++ )
01866 {
01867 part = (*m_expHitCol[i]).Part();
01868 segment = (*m_expHitCol[i]).Segment();
01869 layer = (*m_expHitCol[i]).Layer();
01870 strip = (*m_expHitCol[i]).Strip();
01871
01872 FillExpHit( part, segment, layer, strip );
01873 m_record[part][segment][layer][strip][4] ++;
01874 m_fTotalExpHit ++;
01875 }
01876
01877
01878 log << MSG::INFO << "Fill eff hits" << endreq;
01879 for( unsigned int i=0; i<m_effHitCol.size(); i++ )
01880 {
01881 part = (*m_effHitCol[i]).Part();
01882 segment = (*m_effHitCol[i]).Segment();
01883 layer = (*m_effHitCol[i]).Layer();
01884 strip = (*m_effHitCol[i]).Strip();
01885
01886 FillEffHit( part, segment, layer, strip );
01887 m_fTotalEffHit ++;
01888 }
01889
01890
01891 log << MSG::INFO << "Fill clusters" << endreq;
01892 for( unsigned int i=0; i<m_clusterCol.size(); i++ )
01893 {
01894 size = m_clusterCol[i].size();
01895
01896 if( size > CLUSTER_CUT )
01897 {
01898 part = (*m_clusterCol[i][0]).Part();
01899 segment = (*m_clusterCol[i][0]).Segment();
01900 layer = (*m_clusterCol[i][0]).Layer();
01901
01902 FillCluster( part, segment, layer, size );
01903 m_ntClusterSize = size;
01904 m_clusterSizeTuple->write();
01905 }
01906 }
01907
01908
01909 log << MSG::INFO << "Fill inc/noise hits" << endreq;
01910 for( unsigned int i=0; i<m_nosHitCol.size(); i++ )
01911 {
01912 part = (*m_nosHitCol[i]).Part();
01913 segment = (*m_nosHitCol[i]).Segment();
01914 layer = (*m_nosHitCol[i]).Layer();
01915 strip = (*m_nosHitCol[i]).Strip();
01916
01917 FillNosHit( part, segment, layer, strip );
01918 m_record[part][segment][layer][strip][3] ++;
01919 }
01920
01921
01922 m_ntEventId = m_currentEvent;
01923 m_ntEventTag = m_eventTag;
01924 m_ntDigiNum = m_digiCol.size();
01925 m_ntExpHitNum = m_expHitCol.size();
01926 m_ntEffHitNum = m_effHitCol.size();
01927 m_ntNosHitNum = m_nosHitCol.size();
01928 m_ntClusterNum = m_clusterCol.size();
01929
01930
01931 for( unsigned int i=0; i<m_digiCol.size(); i++ ) {
01932 if( m_digiCol[i] != NULL ) delete m_digiCol[i];
01933 }
01934
01935 for( unsigned int i=0; i<m_expHitCol.size(); i++ ) {
01936 if( m_expHitCol[i] != NULL ) delete m_expHitCol[i];
01937 }
01938
01939
01940
01941
01942
01943
01944 for( unsigned int i=0; i<m_calHitCol.size(); i++ ) {
01945 if( m_calHitCol[i] != NULL ) delete m_calHitCol[i];
01946 }
01947
01948 if( m_digiCol.size() != 0 ) m_digiCol.clear();
01949 if( m_expHitCol.size() != 0 ) m_expHitCol.clear();
01950 if( m_calHitCol.size() != 0 ) m_calHitCol.clear();
01951 if( m_effHitCol.size() != 0 ) m_effHitCol.clear();
01952 if( m_nosHitCol.size() != 0 ) m_nosHitCol.clear();
01953 if( m_clusterCol.size() != 0 ) m_clusterCol.clear();
01954
01955
01956 for( int i=0; i<PART_MAX; i++ ) {
01957 for( int j=0; j<SEGMENT_MAX; j++ ) {
01958 if( m_segDigiCol[i][j].size() != 0 ) m_segDigiCol[i][j].clear();
01959 }
01960 }
01961
01962 m_evtEnd = clock();
01963
01964 m_ntEventTime = (double(m_evtEnd - m_evtBegin))/CLOCKS_PER_SEC;
01965 log << MSG::INFO << "Event time:\t" << m_ntEventTime << "s" << endreq;
01966 m_eventLogTuple->write();
01967
01968 return StatusCode::SUCCESS;
01969 }
01970
01971 StatusCode MucCalibMgr::FillDigi( int part, int segment, int layer, int strip )
01972 {
01973
01974 int tmpId;
01975
01976 if( (int)m_tag || m_dimuOnly==0 || (m_dimuOnly==1&&m_eventTag==1) )
01977 {
01978 if( part == BRID )
01979 {
01980
01981 if( layer%2 == 0 ) {
01982 if( segment > 4 ) tmpId = segment*48 + (47 - strip);
01983 else tmpId = segment*48 + strip;
01984 }
01985 else if( segment < 3 ) tmpId = segment*96 + strip;
01986 else tmpId = (segment-1)*96 + 112 + strip;
01987
01988 m_hHitMapBarrel_Lay[layer]->Fill(tmpId);
01989
01990
01991 if( segment != B_TOP ) {
01992 if( segment > 4 )
01993 tmpId = 48*((layer+1)/2) + 96*(layer/2) + ( ((layer+1)%2)*48 + (layer%2)*96 - strip - 1 );
01994 else tmpId = 48*((layer+1)/2) + 96*(layer/2) + strip;
01995 }
01996 else tmpId = 48*((layer+1)/2) + 112*(layer/2) + strip;
01997
01998 m_hHitMapBarrel_Seg[segment]->Fill(tmpId);
01999 }
02000 else if( part == EEID )
02001 {
02002
02003 tmpId = segment*64 + strip;
02004 m_hHitMapEndcap_Lay[0][layer]->Fill(tmpId);
02005
02006 tmpId = layer*64 + strip;
02007 m_hHitMapEndcap_Seg[0][segment]->Fill(tmpId);
02008 }
02009 else if( part == EWID )
02010 {
02011
02012 tmpId = segment*64 + strip;
02013 m_hHitMapEndcap_Lay[1][layer]->Fill(tmpId);
02014
02015 tmpId = layer*64 + strip;
02016 m_hHitMapEndcap_Seg[1][segment]->Fill(tmpId);
02017 }
02018 }
02019
02020
02021 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
02022 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
02023
02024 m_hStripFireMap[boxId]->Fill( strip );
02025 m_hStripFire->Fill( strId );
02026 m_hBoxFire->Fill( boxId );
02027 m_hLayerFire->Fill( layer );
02028 if( part==BRID ) m_hBrLayerFire->Fill( layer );
02029 else if( part== EEID ) m_hEcLayerFire->Fill( layer+1 );
02030 else m_hEcLayerFire->Fill( -(layer+1) );
02031
02032 return StatusCode::SUCCESS;
02033 }
02034
02035 StatusCode MucCalibMgr::FillExpHit( int part, int segment, int layer, int strip )
02036 {
02037 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
02038 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
02039
02040 m_hStripExpHitMap[boxId]->Fill( strip );
02041 m_hStripExpHit->Fill( strId );
02042 m_hBoxExpHit->Fill( boxId );
02043 m_hLayerExpHit->Fill( layer );
02044
02045 return StatusCode::SUCCESS;
02046 }
02047
02048 StatusCode MucCalibMgr::FillEffHit( int part, int segment, int layer, int strip )
02049 {
02050 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
02051 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
02052
02053 m_hStripEffHitMap[boxId]->Fill( strip );
02054 m_hStripEffHit->Fill( strId );
02055 m_hBoxEffHit->Fill( boxId );
02056 m_hLayerEffHit->Fill( layer );
02057
02058 return StatusCode::SUCCESS;
02059 }
02060
02061 StatusCode MucCalibMgr::FillNosHit( int part, int segment, int layer, int strip )
02062 {
02063 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
02064 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
02065
02066 m_hStripNosHitMap[boxId]->Fill( strip );
02067 m_hStripNosHit->Fill( strId );
02068 m_hBoxNosHit->Fill( boxId );
02069 m_hLayerNosHit->Fill( layer );
02070
02071 return StatusCode::SUCCESS;
02072 }
02073
02074 StatusCode MucCalibMgr::FillCluster( int part, int segment, int layer, int size )
02075 {
02076 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
02077
02078 m_hLayerCluster[layer]->Fill( size );
02079 m_hBoxCluster[boxId]->Fill( size );
02080
02081 return StatusCode::SUCCESS;
02082 }
02083
02084
02085
02086
02087
02088
02089 StatusCode MucCalibMgr::AnalyseEffAndNoise()
02090 {
02091 MsgStream log(msgSvc, "MucCalibMgr");
02092 log<< MSG::INFO << "Start efficiency and noise calibration" << endreq;
02093
02094
02095 m_fTotalDAQTime = m_fTotalEvent/m_er;
02096 log<< MSG::INFO << "Total run time:\t" << m_fTotalDAQTime << endreq;
02097
02098 EffAndNoiseLV2();
02099 EffAndNoiseLV1();
02100 EffAndNoiseLV0();
02101
02102 if( m_usePad != 0 ) PadEff();
02103
02104 return StatusCode::SUCCESS;
02105 }
02106
02107 StatusCode MucCalibMgr::EffAndNoiseLV0()
02108 {
02109 MsgStream log(msgSvc, "MucCalibMgr");
02110 log<< MSG::INFO << "Analyse layer efficiency and noise" << endreq;
02111
02112 int part, segment, layer, stripMax;
02113 part = segment = layer = stripMax = 0;
02114
02115 double digi, effHit, trackNum, nosHit, recHit;
02116 double eff, effErr, noise, nosRatio, nosRatioErr, cnt, cluster;
02117 eff = effErr = noise = nosRatio = nosRatioErr = cnt = cluster = 0.;
02118
02119 for( int i=0; i<LAYER_MAX; i++ )
02120 {
02121 digi = effHit = trackNum = nosHit = recHit = 0;
02122
02123 for( int j=0; j<PART_MAX; j++ )
02124 {
02125 for( int k=0; k<SEGMENT_MAX; k++)
02126 {
02127 stripMax = m_ptrIdTr->GetStripMax( j, k, i );
02128 for( int n=0; n<stripMax; n++ )
02129 {
02130 digi += m_record[j][k][i][n][0];
02131 trackNum += m_record[j][k][i][n][1];
02132 effHit += m_record[j][k][i][n][2];
02133 nosHit += m_record[j][k][i][n][3];
02134 recHit += m_record[j][k][i][n][4];
02135 }
02136 }
02137 }
02138
02139
02140 if( trackNum == 0 ) {
02141 eff = effErr = 0.0;
02142
02143
02144 }
02145 else
02146 {
02147 eff = ( (double)effHit )/trackNum;
02148 effErr = sqrt( eff*(1-eff)/trackNum );
02149 m_fCalibLayerNum ++;
02150 }
02151
02152
02153 if( m_layerResults[3][i] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
02154 noise = DEFAULT_NOS_VALUE;
02155 else {
02156 if( m_recMode == 2 ) noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW*m_layerResults[3][i]);
02157 else noise = (double)nosHit/(m_fTotalDAQTime * m_layerResults[3][i]);
02158 }
02159
02160 if( digi != 0 ) {
02161 nosRatio = ( (double)nosHit )/digi;
02162 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
02163 }
02164 else {
02165 nosRatio = DEFAULT_INC_VALUE;
02166 nosRatioErr = 0;
02167 }
02168
02169
02170 if( m_recMode == 2 )
02171 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_layerResults[3][i]);
02172 else
02173 cnt = (double)digi/(m_fTotalDAQTime * m_layerResults[3][i]);
02174
02175
02176 cluster = m_hLayerCluster[i]->GetMean();
02177
02178 m_layerResults[0][ i ] = eff;
02179 m_layerResults[1][ i ] = effErr;
02180 m_layerResults[2][ i ] = noise;
02181 m_layerResults[4][ i ] = cluster;
02182 m_layerResults[5][ i ] = trackNum;
02183
02184
02185 m_hLayerEff->Fill( i, eff );
02186 m_hLayerEff->SetBinError( i+1, effErr );
02187 m_hLayerNosRatio->Fill( i, nosRatio );
02188 m_hLayerNosRatio->SetBinError( i+1, nosRatioErr );
02189 m_hLayerNos->Fill( i, noise );
02190 m_hLayerCnt->Fill( i, cnt );
02191
02192
02193 m_hLayerClusterCmp->Fill( i, cluster );
02194
02195
02196 m_fLayerId = i;
02197 m_fLayerEff = eff;
02198 m_fLayerEffErr = effErr;
02199 m_fLayerTrkNum = trackNum;
02200 m_fLayerExpHit = recHit;
02201 m_fLayerEffHit = effHit;
02202 m_fLayerNosRatio = nosRatio;
02203 m_fLayerDigi = digi;
02204 m_fLayerNosHit = nosHit;
02205 m_fLayerNos = noise;
02206 m_fLayerCnt = cnt;
02207 m_tLayConst->Fill();
02208
02209
02210 m_fLayerCluster = cluster;
02211
02212 }
02213
02214 return StatusCode::SUCCESS;
02215 }
02216
02217 StatusCode MucCalibMgr::EffAndNoiseLV1()
02218 {
02219 MsgStream log(msgSvc, "MucCalibMgr");
02220 log<< MSG::INFO << "Analyse box efficiency and noise" << endreq;
02221
02222 int part, segment, layer, stripMax;
02223 part = segment = layer = stripMax = 0;
02224
02225 double digi, effHit, trackNum, nosHit, recHit;
02226 double eff, effErr, noise, nosRatio, nosRatioErr, cnt, cluster;
02227 eff = effErr = noise = nosRatio = nosRatioErr = cnt = cluster = 0.;
02228
02229 for( int i=0; i<BOX_MAX; i++ )
02230 {
02231 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
02232 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
02233
02234 digi = effHit = trackNum = nosHit = recHit = 0;
02235 for( int j=0; j<stripMax; j++ )
02236 {
02237 digi += m_record[part][segment][layer][j][0];
02238 trackNum += m_record[part][segment][layer][j][1];
02239 effHit += m_record[part][segment][layer][j][2];
02240 nosHit += m_record[part][segment][layer][j][3];
02241 recHit += m_record[part][segment][layer][j][4];
02242 }
02243
02244
02245 if( trackNum == 0 ) {
02246 eff = effErr = 0.0;
02247
02248
02249 }
02250 else
02251 {
02252 eff = ( (double)effHit )/trackNum;
02253 effErr = sqrt( eff*(1-eff)/trackNum );
02254 m_fCalibBoxNum ++;
02255 }
02256
02257
02258 if( m_boxResults[3][i] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
02259 noise = DEFAULT_NOS_VALUE;
02260 else {
02261 if( m_recMode == 2 )
02262 noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW * m_boxResults[3][i]);
02263 else
02264 noise = (double)nosHit/(m_fTotalDAQTime * m_boxResults[3][i]);
02265 }
02266
02267 if( digi != 0 ) {
02268 nosRatio = ( (double)nosHit )/digi;
02269 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
02270 }
02271 else {
02272 nosRatio = DEFAULT_INC_VALUE;
02273 nosRatioErr = 0;
02274 }
02275
02276
02277 if( m_recMode == 2 )
02278 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_boxResults[3][i]);
02279 else
02280 cnt = (double)digi/(m_fTotalDAQTime * m_boxResults[3][i]);
02281
02282
02283 cluster = m_hBoxCluster[i]->GetMean();
02284
02285 m_boxResults[0][ i ] = eff;
02286 m_boxResults[1][ i ] = effErr;
02287 m_boxResults[2][ i ] = noise;
02288 m_boxResults[4][ i ] = cluster;
02289 m_boxResults[5][ i ] = trackNum;
02290
02291
02292 m_hBoxEff->Fill( i, eff );
02293 m_hBoxEff->SetBinError( i+1, effErr );
02294 m_hBoxNosRatio->Fill( i, nosRatio );
02295 m_hBoxNosRatio->SetBinError( i+1, nosRatioErr );
02296 m_hBoxNos->Fill( i, noise );
02297 m_hBoxCnt->Fill( i, cnt );
02298
02299
02300 m_hBoxClusterCmp->Fill( i, cluster );
02301
02302
02303 m_fBoxId = i;
02304 m_fBoxPart = part;
02305 m_fBoxSegment = segment;
02306 m_fBoxLayer = layer;
02307 m_fBoxEff = eff;
02308 m_fBoxEffErr = effErr;
02309 m_fBoxTrkNum = trackNum;
02310 m_fBoxExpHit = recHit;
02311 m_fBoxEffHit = effHit;
02312 m_fBoxNosRatio = nosRatio;
02313 m_fBoxDigi = digi;
02314 m_fBoxNosHit = nosHit;
02315 m_fBoxNos = noise;
02316 m_fBoxCnt = cnt;
02317 m_tBoxConst->Fill();
02318
02319
02320 m_fBoxCluster = cluster;
02321
02322 }
02323
02324 return StatusCode::SUCCESS;
02325 }
02326
02327 StatusCode MucCalibMgr::EffAndNoiseLV2()
02328 {
02329 MsgStream log(msgSvc, "MucCalibMgr");
02330 log<< MSG::INFO << "Analyse strip efficiency and noise" << endreq;
02331
02332 int part, segment, layer, stripMax;
02333 part = segment = layer = stripMax = 0;
02334
02335 double digi, effHit, trackNum, nosHit, recHit;
02336 double eff, effErr, noise, nosRatio, nosRatioErr, cnt, cluster;
02337 eff = effErr = noise = nosRatio = nosRatioErr = cnt = cluster = 0.;
02338
02339 for( int i=0; i<BOX_MAX; i++ )
02340 {
02341 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
02342 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
02343
02344 for( int j=0; j<stripMax; j++ )
02345 {
02346 digi = m_record[part][segment][layer][j][0];
02347 trackNum = m_record[part][segment][layer][j][1];
02348 effHit = m_record[part][segment][layer][j][2];
02349 nosHit = m_record[part][segment][layer][j][3];
02350 recHit = m_record[part][segment][layer][j][4];
02351
02352 int stripId = m_ptrIdTr->GetStripId( part, segment, layer, j );
02353
02354
02355 if( trackNum == 0 ) {
02356 eff = effErr = 0.0;
02357
02358
02359 }
02360 else
02361 {
02362 eff = ( (double)effHit )/trackNum;
02363 effErr = sqrt( eff*(1-eff)/trackNum );
02364 m_fCalibStripNum ++;
02365 }
02366
02367
02368 if( m_stripResults[3][stripId] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
02369 noise = DEFAULT_NOS_VALUE;
02370 else {
02371 if( m_recMode == 2 )
02372 noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW * m_stripResults[3][stripId]);
02373 else
02374 noise = (double)nosHit/(m_fTotalDAQTime * m_stripResults[3][stripId]);
02375 }
02376
02377 if( digi != 0 ) {
02378 nosRatio = ( (double)nosHit )/digi;
02379 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
02380 }
02381 else {
02382 nosRatio = DEFAULT_INC_VALUE;
02383 nosRatioErr = 0.;
02384 }
02385
02386
02387 if( m_recMode == 2 )
02388 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_stripResults[3][stripId]);
02389 else
02390 cnt = (double)digi/(m_fTotalDAQTime * m_stripResults[3][stripId]);
02391
02392
02393 m_stripResults[0][ stripId ] = eff;
02394 m_stripResults[1][ stripId ] = effErr;
02395 m_stripResults[2][ stripId ] = noise;
02396 m_stripResults[5][ stripId ] = trackNum;
02397
02398
02399 m_hStripEffMap[i]->Fill( j, eff );
02400 m_hStripEffMap[i]->SetBinError( j+1, effErr );
02401 m_hStripEff->Fill( stripId, eff );
02402 m_hStripEff->SetBinError( stripId+1, effErr );
02403 m_hStripNosRatioMap[i]->Fill( j, nosRatio );
02404 m_hStripNosRatioMap[i]->SetBinError( j+1, nosRatioErr );
02405 m_hStripNosRatio->Fill( stripId, nosRatio );
02406 m_hStripNosRatio->SetBinError( stripId+1, nosRatioErr );
02407 m_hStripNos->Fill( stripId, noise );
02408 m_hStripCnt->Fill( stripId, cnt );
02409
02410
02411 m_fStripId = stripId;
02412 m_fStripPart = part;
02413 m_fStripSegment = segment;
02414 m_fStripLayer = layer;
02415 m_fStripEff = eff;
02416 m_fStripEffErr = effErr;
02417 m_fStripNosRatio = nosRatio;
02418 m_fStripDigi = digi;
02419 m_fStripNos = noise;
02420 m_fStripCnt = cnt;
02421 m_fStripEffHit = effHit;
02422 m_fStripExpHit = recHit;
02423 m_fStripNosHit = nosHit;
02424 m_fStripTrkNum = trackNum;
02425 m_tStrConst->Fill();
02426
02427 }
02428 }
02429
02430 return StatusCode::SUCCESS;
02431 }
02432
02433 StatusCode MucCalibMgr::PadEff()
02434 {
02435 MsgStream log(msgSvc, "MucCalibMgr");
02436
02437 int xBinStart, xBinEnd, yBinStart, yBinEnd;
02438 double expHit, firedHit, eff;
02439 eff = expHit = firedHit = 0.;
02440
02441 for( int i=0; i<PART_MAX; i++ )
02442 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
02443 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
02444 {
02445 xBinStart = m_h2DExpMap[i][j][k]->GetXaxis()->GetFirst();
02446 xBinEnd = m_h2DExpMap[i][j][k]->GetXaxis()->GetLast() + 1;
02447 yBinStart = m_h2DExpMap[i][j][k]->GetYaxis()->GetFirst();
02448 yBinEnd = m_h2DExpMap[i][j][k]->GetYaxis()->GetLast() + 1;
02449
02450 for( int xi = xBinStart; xi<xBinEnd; xi++ )
02451 for( int yi = yBinStart; yi<yBinEnd; yi++ )
02452 {
02453 expHit = m_h2DExpMap[i][j][k]->GetBinContent(xi, yi);
02454 firedHit = m_h2DHitMap[i][j][k]->GetBinContent(xi, yi);
02455
02456 if( expHit !=0 ) eff = firedHit / expHit;
02457 else eff = 0;
02458
02459 if( eff>1.0 )
02460 cout<<"Eff error:\t["<<i<<"\t"<<j<<"\t"<<k<<",\t"<<xi<<"\t"<<yi<<"]:\t"
02461 <<eff<<"\t,"<<firedHit<<" "<<expHit<<endl;
02462
02463 m_h2DEffMap[i][j][k]->SetBinContent(xi, yi, eff);
02464 }
02465 }
02466 return StatusCode::SUCCESS;
02467 }
02468
02469 StatusCode MucCalibMgr::AnalyseCluster()
02470 {
02471 MsgStream log(msgSvc, "MucCalibMgr");
02472
02473 return StatusCode::SUCCESS;
02474 }
02475
02476 StatusCode MucCalibMgr::AnalyseRes()
02477 {
02478 MsgStream log(msgSvc, "MucCalibMgr");
02479 log << MSG::INFO << "Analyse spacial resolution" << endreq;
02480
02481 double resMean, resMeanErr, resSigma, resSigmaErr;
02482 resMean = resMeanErr = resSigma = resSigmaErr = 0.;
02483
02484 for( int i=0; i<B_LAY_NUM; i++ )
02485 {
02486 m_hBarrelResDist[i]->Fit("gaus");
02487 resMean = m_hBarrelResDist[i]->GetFunction("gaus")->GetParameter("Mean");
02488 resSigma = m_hBarrelResDist[i]->GetFunction("gaus")->GetParameter("Sigma");
02489 resMeanErr = m_hBarrelResDist[i]->GetFunction("gaus")->GetParError(1);
02490 resSigmaErr = m_hBarrelResDist[i]->GetFunction("gaus")->GetParError(2);
02491
02492 m_hBarrelResComp[0]->Fill( i, resMean );
02493 m_hBarrelResComp[0]->SetBinError( i+1, resMeanErr );
02494 m_hBarrelResComp[1]->Fill( i, resSigma );
02495 m_hBarrelResComp[1]->SetBinError( i+1, resSigmaErr );
02496 }
02497
02498 for( int i=0; i<E_LAY_NUM; i++ )
02499 {
02500 m_hEndcapResDist[i]->Fit("gaus");
02501 resMean = m_hEndcapResDist[i]->GetFunction("gaus")->GetParameter("Mean");
02502 resSigma = m_hEndcapResDist[i]->GetFunction("gaus")->GetParameter("Sigma");
02503 resMeanErr = m_hEndcapResDist[i]->GetFunction("gaus")->GetParError(1);
02504 resSigmaErr = m_hEndcapResDist[i]->GetFunction("gaus")->GetParError(2);
02505
02506 m_hEndcapResComp[0]->Fill( i, resMean );
02507 m_hEndcapResComp[0]->SetBinError( i+1, resMeanErr );
02508 m_hEndcapResComp[1]->Fill( i, resSigma );
02509 m_hEndcapResComp[1]->SetBinError( i+1, resSigmaErr );
02510 }
02511
02512 return StatusCode::SUCCESS;
02513 }
02514
02515 StatusCode MucCalibMgr::SaveConst()
02516 {
02517 MsgStream log(msgSvc, "MucCalibMgr");
02518 log << MSG::INFO << "Save calibration constants" << endreq;
02519
02520
02521
02522 double layerXY[2][LAYER_MAX];
02523 double layerEXY[2][LAYER_MAX];
02524 for( int i=0; i<LAYER_MAX; i++ )
02525 {
02526 layerXY[0][i] = i;
02527 layerEXY[0][i] = 0;
02528 if( m_layerResults[5][i] >= 100*TRACK_THRESHOLD ) {
02529 layerXY[1][i] = m_layerResults[0][i];
02530 layerEXY[1][i] = m_layerResults[1][i];
02531 }
02532 else {
02533 layerXY[1][i] = 0;
02534 layerEXY[1][i] = 0;
02535 }
02536 }
02537 m_geLayerEff = new TGraphErrors(LAYER_MAX, layerXY[0], layerXY[1], layerEXY[0], layerEXY[1]);
02538 m_geLayerEff->SetMarkerStyle(25);
02539 m_geLayerEff->SetMarkerSize(0.5);
02540 m_cv[0] = new TCanvas("GoodLayerEff","Layer efficiency", 50, 50, 800, 600);
02541 m_cv[0]->SetFillColor(0);
02542 m_cv[0]->SetBorderMode(0);
02543 m_geLayerEff->Draw("AP");
02544 m_cv[0]->Write();
02545
02546
02547 double boxXY[2][BOX_MAX];
02548 double boxEXY[2][BOX_MAX];
02549 for( int i=0; i<BOX_MAX; i++ )
02550 {
02551 boxXY[0][i] = i;
02552 boxEXY[0][i] = 0;
02553 if( m_boxResults[5][i] >= 10*TRACK_THRESHOLD ) {
02554 boxXY[1][i] = m_boxResults[0][i];
02555 boxEXY[1][i] = m_boxResults[1][i];
02556 }
02557 else {
02558 boxXY[1][i] = 0;
02559 boxEXY[1][i] = 0;
02560 }
02561 }
02562 m_geBoxEff = new TGraphErrors(BOX_MAX, boxXY[0], boxXY[1], boxEXY[0], boxEXY[1]);
02563 m_geBoxEff->SetMarkerStyle(25);
02564 m_geBoxEff->SetMarkerSize(0.5);
02565 m_cv[1] = new TCanvas("GoodBoxEff","Box efficiency", 75, 75, 800, 600);
02566 m_cv[1]->SetFillColor(0);
02567 m_cv[1]->SetBorderMode(0);
02568 m_geBoxEff->Draw("AP");
02569 m_cv[1]->Write();
02570
02571
02572 double stripXY[2][STRIP_MAX];
02573 double stripEXY[2][STRIP_MAX];
02574 for( int i=0; i<STRIP_MAX; i++ )
02575 {
02576 stripXY[0][i] = i;
02577 stripEXY[0][i] = 0;
02578 if( m_stripResults[5][i] >= TRACK_THRESHOLD ) {
02579 stripXY[1][i] = m_stripResults[0][i];
02580 stripEXY[1][i] = m_stripResults[1][i];
02581 }
02582 else {
02583 stripXY[1][i] = 0;
02584 stripEXY[1][i] = 0;
02585 }
02586 }
02587 m_geStripEff = new TGraphErrors(STRIP_MAX, stripXY[0], stripXY[1], stripEXY[0], stripEXY[1]);
02588 m_geStripEff->SetMarkerStyle(25);
02589 m_geStripEff->SetMarkerSize(0.5);
02590 m_cv[2] = new TCanvas("GoodStripEff","Strip efficiency", 100, 100, 800, 600);
02591 m_cv[2]->SetFillColor(0);
02592 m_cv[2]->SetBorderMode(0);
02593 m_geStripEff->Draw("AP");
02594 m_cv[2]->Write();
02595
02596
02597 for(int i=0; i<B_LAY_NUM; i++ )
02598 {
02599 m_hHitMapBarrel_Lay[i]->Write();
02600
02601 if( i<E_LAY_NUM ) {
02602 m_hHitMapEndcap_Lay[0][i]->Write();
02603 m_hHitMapEndcap_Lay[1][i]->Write();
02604 }
02605 }
02606
02607 for(int i=0; i<B_SEG_NUM; i++)
02608 {
02609 m_hHitMapBarrel_Seg[i]->Write();
02610
02611 if( i< E_SEG_NUM ) {
02612 m_hHitMapEndcap_Seg[0][i]->Write();
02613 m_hHitMapEndcap_Seg[1][i]->Write();
02614 }
02615 }
02616 m_hTrackDistance->Fit("gaus");
02617 m_hTrackPosPhiDiff->Fit("gaus");
02618 m_hTrackPosThetaDiff->Fit("gaus");
02619 m_hTrackMomPhiDiff->Fit("gaus");
02620 m_hTrackMomThetaDiff->Fit("gaus");
02621
02622 m_hTrackDistance->Write();
02623 m_hTrackPosPhiDiff->Write();
02624 m_hTrackPosThetaDiff->Write();
02625 m_hTrackMomPhiDiff->Write();
02626 m_hTrackMomThetaDiff->Write();
02627
02628 for( int i=0; i<B_LAY_NUM; i++ ) m_hBarrelResDist[i]->Write();
02629 for( int i=0; i<E_LAY_NUM; i++ ) m_hEndcapResDist[i]->Write();
02630
02631 m_hBarrelResComp[0]->Write();
02632 m_hBarrelResComp[1]->Write();
02633 m_hEndcapResComp[0]->Write();
02634 m_hEndcapResComp[1]->Write();
02635
02636 if( m_usePad != 0 ) m_histArray->Write();
02637
02638 for( int i=0; i<BOX_MAX; i++ )
02639 {
02640 m_hStripFireMap[ i ]->Write();
02641 m_hStripExpHitMap[ i ]->Write();
02642 m_hStripEffHitMap[ i ]->Write();
02643 m_hStripNosHitMap[ i ]->Write();
02644 m_hStripEffMap[ i ]->Write();
02645 m_hStripNosRatioMap[ i ]->Write();
02646 }
02647 m_hStripFire->Write();
02648 m_hStripExpHit->Write();
02649 m_hStripEffHit->Write();
02650 m_hStripNosHit->Write();
02651 m_hStripEff->Write();
02652 m_hStripArea->Write();
02653 m_hStripNos->Write();
02654 m_hStripNosRatio->Write();
02655 log << MSG::INFO << "Save LV2 histograms done!" << endreq;
02656
02657 m_hBoxFire->Write();
02658 m_hBoxExpHit->Write();
02659 m_hBoxEffHit->Write();
02660 m_hBoxNosHit->Write();
02661 m_hBoxEff->Write();
02662 m_hBoxArea->Write();
02663 m_hBoxNos->Write();
02664 m_hBoxNosRatio->Write();
02665 log << MSG::INFO << "Save LV1 histograms done!" << endreq;
02666
02667 m_hBrLayerFire->Write();
02668 m_hEcLayerFire->Write();
02669 m_hLayerFire->Write();
02670 m_hLayerExpHit->Write();
02671 m_hLayerEffHit->Write();
02672 m_hLayerNosHit->Write();
02673 m_hLayerEff->Write();
02674 m_hLayerArea->Write();
02675 m_hLayerNos->Write();
02676 m_hLayerNosRatio->Write();
02677
02678 for( int i=0; i<LAYER_MAX; i++ ) m_hLayerCluster[i]->Write();
02679 for( int i=0; i<BOX_MAX; i++ ) m_hBoxCluster[i]->Write();
02680 m_hLayerClusterCmp->Write();
02681 m_hBoxClusterCmp->Write();
02682
02683 log << MSG::INFO << "Save histograms done!" << endreq;
02684
02685
02686 m_fLayerCoverage = 100*(double)m_fCalibLayerNum/LAYER_MAX;
02687 m_fBoxCoverage = 100*(double)m_fCalibBoxNum/BOX_MAX;
02688 m_fStripCoverage = 100*(double)m_fCalibStripNum/STRIP_MAX;
02689
02690 long digi_num, trk_num, eff_hit, nos_hit, exp_hit;
02691 m_tStatLog->Branch("digi_num", &digi_num, "digi_num/I");
02692 m_tStatLog->Branch("trk_num", &trk_num, "trk_num/I");
02693 m_tStatLog->Branch("eff_hit", &eff_hit, "eff_hit/I");
02694 m_tStatLog->Branch("nos_hit", &nos_hit, "nos_hit/I");
02695 m_tStatLog->Branch("exp_hit", &exp_hit, "exp_hit/I");
02696
02697 int stripMax;
02698 for( int i=0; i<PART_MAX; i++ )
02699 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
02700 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
02701 {
02702 stripMax = m_ptrIdTr->GetStripMax( i, j, k );
02703 for( int n=0; n<stripMax; n++ )
02704 {
02705 digi_num = m_record[i][j][k][n][0];
02706 trk_num = m_record[i][j][k][n][1];
02707 eff_hit = m_record[i][j][k][n][2];
02708 nos_hit = m_record[i][j][k][n][3];
02709 exp_hit = m_record[i][j][k][n][4];
02710 m_tStatLog->Fill();
02711 }
02712 }
02713
02714 m_jobFinish = clock();
02715 m_fTotalJobTime = (double)(m_jobFinish - m_jobStart)/CLOCKS_PER_SEC;
02716
02717 m_tJobLog->Fill();
02718 m_tJobLog->Write();
02719 m_tStatLog->Write();
02720 m_tLayConst->Write();
02721 m_tBoxConst->Write();
02722 m_tStrConst->Write();
02723
02724
02725 if( m_fdata != NULL ) m_fdata->close();
02726
02727 return StatusCode::SUCCESS;
02728 }
02729
02730 StatusCode MucCalibMgr::EndRun()
02731 {
02732 MsgStream log(msgSvc, "MucCalibMgr");
02733 log << MSG::INFO << endreq;
02734 log << MSG::INFO << "Total events : " << m_fTotalEvent << endreq;
02735 log << MSG::INFO << "Total digis : " << m_fTotalDigi << endreq;
02736 log << MSG::INFO << "Total rec hits : " << m_fTotalExpHit << endreq;
02737 log << MSG::INFO << "Total eff hits : " << m_fTotalEffHit << endreq;
02738 log << MSG::INFO << "Total inc hits : " << m_fTotalNosHit << endreq;
02739 log << MSG::INFO << "Total clusters : " << m_fTotalClstNum<< endreq;
02740
02741 log << MSG::INFO << "Strip calibrated percentage: "
02742 << 100*(double)m_fCalibStripNum/STRIP_MAX << "%" << endreq;
02743 log << MSG::INFO << "Box calibrated percentage: "
02744 << 100*(double)m_fCalibBoxNum/BOX_MAX << "%" << endreq;
02745 log << MSG::INFO << "Layer calibrated percentage: "
02746 << 100*(double)m_fCalibLayerNum/LAYER_MAX << "%" << endreq;
02747
02748 log << MSG::INFO << "Calibration run successfully" << endreq << endreq;
02749
02750 return StatusCode::SUCCESS;
02751 }
02752
02753
