12#include "CLHEP/Vector/LorentzVector.h"
13#include "CLHEP/Vector/ThreeVector.h"
14#include "CLHEP/Vector/LorentzVector.h"
15#include "CLHEP/Vector/TwoVector.h"
16#include "CLHEP/Geometry/Point3D.h"
17using CLHEP::Hep3Vector;
18using CLHEP::Hep2Vector;
19using CLHEP::HepLorentzVector;
21#include "GaudiKernel/IHistogramSvc.h"
22#include "GaudiKernel/MsgStream.h"
23#include "GaudiKernel/AlgFactory.h"
24#include "GaudiKernel/ISvcLocator.h"
25#include "GaudiKernel/SmartDataPtr.h"
26#include "GaudiKernel/SmartDataLocator.h"
27#include "GaudiKernel/IDataProviderSvc.h"
28#include "GaudiKernel/Bootstrap.h"
29#include "GaudiKernel/PropertyMgr.h"
30#include "GaudiKernel/INTupleSvc.h"
31#include "GaudiKernel/NTuple.h"
32#include "GaudiKernel/Algorithm.h"
67 MsgStream log(
msgSvc,
"MucCalibMgr");
68 log << MSG::INFO <<
"-------Initializing------- " << endreq;
71 log << MSG::INFO <<
"Initialize Gaudi service" << endreq;
72 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
79 log << MSG::INFO <<
"Initialize parameters" << endreq;
81 m_vs = jobInfo[0]; m_hv = jobInfo[1]; m_th = jobInfo[2];
83 if( jobInfo[3] <= 0 ) m_er = TRIGGER_RATE;
84 else m_er = jobInfo[3];
86 if( jobInfo[4] <= 0 || jobInfo[4] >5 ) m_tag = 0;
87 else m_tag = jobInfo[4];
89 if( configInfo[0] < 0 || configInfo[0] > 2 )
92 m_recMode = configInfo[0];
94 m_usePad = configInfo[1];
96 if( configInfo[2] < 0 || configInfo[2] > STRIP_INBOX_MAX )
97 m_effWindow = EFF_WINDOW;
98 else m_effWindow = configInfo[2];
100 if( configInfo[3]<0 || configInfo[3] >4 )
101 m_clusterMode = DEFAULT_BUILD_MODE;
102 else m_clusterMode = configInfo[3];
104 m_clusterSave = configInfo[4];
105 m_checkEvent = configInfo[5];
106 m_dimuSelect = configInfo[6];
107 m_dimuOnly = configInfo[7];
109 m_outputFile = outputFile;
112 if( m_clusterMode ==1 && m_clusterSave == 1 )
113 m_fdata =
new ofstream(
"MucCluster.dat", ios::out);
115 m_fStartRun = m_fEndRun = 0;
116 m_currentRun = m_currentEvent = 0;
120 for(
int j=0; j<SEGMENT_MAX; j++ )
121 for(
int k=0; k<LAYER_MAX; k++ )
122 for(
int n=0;
n<STRIP_INBOX_MAX;
n++ )
124 m_record[i][j][k][
n][0] = 0;
125 m_record[i][j][k][
n][1] = 0;
126 m_record[i][j][k][
n][2] = 0;
127 m_record[i][j][k][
n][3] = 0;
128 m_record[i][j][k][
n][4] = 0;
131 for(
int i=0; i<LAYER_MAX; i++ )
133 m_layerResults[0][i] = 0;
134 m_layerResults[1][i] = 0;
135 m_layerResults[2][i] = 0;
136 m_layerResults[3][i] = 0;
137 m_layerResults[4][i] = 0;
140 for(
int i=0; i<BOX_MAX; i++ )
142 m_boxResults[0][i] = 0;
143 m_boxResults[1][i] = 0;
144 m_boxResults[2][i] = 0;
145 m_boxResults[3][i] = 0;
146 m_boxResults[4][i] = 0;
151 m_stripResults[0][i] = 0;
152 m_stripResults[1][i] = 0;
153 m_stripResults[2][i] = 0;
154 m_stripResults[3][i] = 0;
160 m_ptrMucMark =
new MucMark(0,0,0,0);
164 m_calHitCol.resize(0);
165 m_expHitCol.resize(0);
166 m_effHitCol.resize(0);
167 m_nosHitCol.resize(0);
168 m_clusterCol.resize(0);
171 for(
int j=0; j<SEGMENT_MAX; j++ ) {
172 m_segDigiCol[i][j].resize(0);
186 log << MSG::INFO <<
"Initialize done!" << endreq;
203 if(m_ptrMucMark)
delete m_ptrMucMark;
204 if(m_ptrIdTr)
delete m_ptrIdTr;
210 m_clusterCol.clear();
213 for(
int j=0; j<SEGMENT_MAX; j++ ) {
214 m_segDigiCol[i][j].clear();
217 if(m_record)
delete []m_record;
218 if(m_layerResults)
delete []m_layerResults;
219 if(m_boxResults)
delete []m_boxResults;
220 if(m_stripResults)
delete []m_stripResults;
226 if(m_hHitMapBarrel_Lay)
delete []m_hHitMapBarrel_Lay;
227 if(m_hHitMapEndcap_Lay)
delete []m_hHitMapEndcap_Lay;
228 if(m_hHitMapBarrel_Seg)
delete []m_hHitMapBarrel_Seg;
229 if(m_hHitMapEndcap_Seg)
delete []m_hHitMapEndcap_Seg;
231 if(m_hHitVsEvent)
delete m_hHitVsEvent;
232 if(m_hTrackDistance)
delete m_hTrackDistance;
233 if(m_hTrackPosPhiDiff)
delete m_hTrackPosPhiDiff;
234 if(m_hTrackPosThetaDiff)
delete m_hTrackPosThetaDiff;
235 if(m_hTrackMomPhiDiff)
delete m_hTrackMomPhiDiff;
236 if(m_hTrackMomThetaDiff)
delete m_hTrackMomThetaDiff;
237 if(m_hDimuTracksPosDiff)
delete m_hDimuTracksPosDiff;
238 if(m_hDimuTracksMomDiff)
delete m_hDimuTracksMomDiff;
239 if(m_hPhiCosTheta)
delete m_hPhiCosTheta;
241 return StatusCode::SUCCESS;
247 if(m_hBrLayerFire)
delete m_hBrLayerFire;
248 if(m_hEcLayerFire)
delete m_hEcLayerFire;
249 if(m_hLayerFire)
delete m_hLayerFire;
250 if(m_hLayerExpHit)
delete m_hLayerExpHit;
251 if(m_hLayerExpHit)
delete m_hLayerEffHit;
252 if(m_hLayerNosHit)
delete m_hLayerNosHit;
253 if(m_hLayerEff)
delete m_hLayerEff;
254 if(m_hLayerNosRatio)
delete m_hLayerNosRatio;
255 if(m_hLayerArea)
delete m_hLayerArea;
256 if(m_hLayerNos)
delete m_hLayerNos;
257 if(m_hLayerCnt)
delete m_hLayerCnt;
259 return StatusCode::SUCCESS;
265 if(m_hBoxFire)
delete m_hBoxFire;
266 if(m_hBoxExpHit)
delete m_hBoxExpHit;
267 if(m_hBoxEffHit)
delete m_hBoxEffHit;
268 if(m_hBoxNosHit)
delete m_hBoxNosHit;
269 if(m_hBoxEff)
delete m_hBoxEff;
270 if(m_hBoxNosRatio)
delete m_hBoxNosRatio;
271 if(m_hBoxArea)
delete m_hBoxArea;
272 if(m_hBoxNos)
delete m_hBoxNos;
273 if(m_hBoxCnt)
delete m_hBoxCnt;
275 return StatusCode::SUCCESS;
281 for(
int i=0; i< BOX_MAX; i++ )
283 if(m_hStripFireMap[i])
delete m_hStripFireMap[i];
284 if(m_hStripExpHitMap[i])
delete m_hStripExpHitMap[i];
285 if(m_hStripEffHitMap[i])
delete m_hStripEffHitMap[i];
286 if(m_hStripNosHitMap[i])
delete m_hStripNosHitMap[i];
287 if(m_hStripEffMap[i])
delete m_hStripEffMap[i];
288 if(m_hStripNosRatioMap[i])
delete m_hStripNosRatioMap[i];
291 if(m_hStripFire)
delete m_hStripFire;
292 if(m_hStripExpHit)
delete m_hStripExpHit;
293 if(m_hStripEffHit)
delete m_hStripEffHit;
294 if(m_hStripNosHit)
delete m_hStripNosHit;
295 if(m_hStripEff)
delete m_hStripEff;
296 if(m_hStripNosRatio)
delete m_hStripNosRatio;
297 if(m_hStripArea)
delete m_hStripArea;
298 if(m_hStripNos)
delete m_hStripNos;
299 if(m_hStripCnt)
delete m_hStripCnt;
301 return StatusCode::SUCCESS;
312 if(m_h2DExpMap[i][j][k])
delete m_h2DExpMap[i][j][k];
313 if(m_h2DHitMap[i][j][k])
delete m_h2DHitMap[i][j][k];
314 if(m_h2DEffMap[i][j][k])
delete m_h2DEffMap[i][j][k];
317 return StatusCode::SUCCESS;
323 for(
int i=0; i<LAYER_MAX; i++ ) {
324 if(m_hLayerCluster[i])
delete m_hLayerCluster[i];
327 for(
int i=0; i<BOX_MAX; i++ ) {
328 if(m_hBoxCluster[i])
delete m_hBoxCluster[i];
331 if(m_hLayerClusterCmp)
delete m_hLayerClusterCmp;
332 if(m_hBoxClusterCmp)
delete m_hBoxClusterCmp;
334 return StatusCode::SUCCESS;
341 if(m_hBarrelResDist[i])
delete m_hBarrelResDist[i];
344 if(m_hEndcapResDist[i])
delete m_hEndcapResDist[i];
347 if(m_hBarrelResComp[0])
delete m_hBarrelResComp[0];
348 if(m_hBarrelResComp[1])
delete m_hBarrelResComp[1];
349 if(m_hEndcapResComp[0])
delete m_hEndcapResComp[0];
350 if(m_hEndcapResComp[1])
delete m_hEndcapResComp[1];
352 return StatusCode::SUCCESS;
359 if(m_tJobLog)
delete m_tJobLog;
360 if(m_tStatLog)
delete m_tStatLog;
361 if(m_tLayConst)
delete m_tLayConst;
362 if(m_tBoxConst)
delete m_tBoxConst;
363 if(m_tStrConst)
delete m_tStrConst;
365 return StatusCode::SUCCESS;
376 MsgStream log(
msgSvc,
"MucCalibMgr");
377 log << MSG::INFO <<
"Initialize NTuples" << endreq;
380 Gaudi::svcLocator() -> service(
"NTupleSvc",
ntupleSvc);
385 NTuplePtr nt1(
ntupleSvc,
"FILE450/EventLog");
386 if ( nt1 ) { m_eventLogTuple = nt1; }
389 m_eventLogTuple =
ntupleSvc->book (
"FILE450/EventLog", CLID_RowWiseTuple,
"MucCalibConst N-Tuple");
390 if ( m_eventLogTuple )
392 sc = m_eventLogTuple->addItem (
"event_id", m_ntEventId);
393 sc = m_eventLogTuple->addItem (
"event_tag", m_ntEventTag);
394 sc = m_eventLogTuple->addItem (
"start_time", m_ntEsTime);
395 sc = m_eventLogTuple->addItem (
"digi_num", m_ntDigiNum);
396 sc = m_eventLogTuple->addItem (
"track_num", m_ntTrackNum);
397 sc = m_eventLogTuple->addItem (
"exphit_num", m_ntExpHitNum);
398 sc = m_eventLogTuple->addItem (
"effhit_num", m_ntEffHitNum);
399 sc = m_eventLogTuple->addItem (
"noshit_num", m_ntNosHitNum);
400 sc = m_eventLogTuple->addItem (
"cluster_num", m_ntClusterNum);
401 sc = m_eventLogTuple->addItem (
"event_time", m_ntEventTime);
404 log << MSG::ERROR <<
"Cannot book N-tuple:" << long(m_eventLogTuple) << endmsg;
405 return StatusCode::FAILURE;
411 NTuplePtr nt2(
ntupleSvc,
"FILE450/MdcTrkInfo");
412 if ( nt2 ) { m_mdcTrkInfoTuple = nt2; }
415 m_mdcTrkInfoTuple =
ntupleSvc->book (
"FILE450/MdcTrkInfo", CLID_RowWiseTuple,
"MucCalibConst N-Tuple");
416 if ( m_mdcTrkInfoTuple )
418 sc = m_mdcTrkInfoTuple->addItem (
"charge", m_charge);
419 sc = m_mdcTrkInfoTuple->addItem (
"mdcpx", m_mdcpx);
420 sc = m_mdcTrkInfoTuple->addItem (
"mdcpy", m_mdcpy);
421 sc = m_mdcTrkInfoTuple->addItem (
"mdcpz", m_mdcpz);
422 sc = m_mdcTrkInfoTuple->addItem (
"mdcpt", m_mdcpt);
423 sc = m_mdcTrkInfoTuple->addItem (
"mdcpp", m_mdcpp);
424 sc = m_mdcTrkInfoTuple->addItem (
"mdcphi", m_mdcphi);
425 sc = m_mdcTrkInfoTuple->addItem (
"mdctheta", m_mdctheta);
428 log << MSG::ERROR <<
"Cannot book N-tuple:" << long(m_mdcTrkInfoTuple) << endmsg;
429 return StatusCode::FAILURE;
434 NTuplePtr nt3(
ntupleSvc,
"FILE450/TrackInfo");
435 if ( nt3 ) { m_trackInfoTuple = nt3; }
438 m_trackInfoTuple =
ntupleSvc->book (
"FILE450/TrackInfo", CLID_RowWiseTuple,
"MucCalibConst N-Tuple");
439 if ( m_trackInfoTuple )
441 sc = m_trackInfoTuple->addItem (
"track_event", m_ntTrackEvent);
442 sc = m_trackInfoTuple->addItem (
"track_tag", m_ntTrackTag);
443 sc = m_trackInfoTuple->addItem (
"track_hits", m_ntTrackHits);
444 sc = m_trackInfoTuple->addItem (
"segment_fly", m_ntTrackSegFly);
445 sc = m_trackInfoTuple->addItem (
"layer_fly_a", m_ntTrackLayFlyA);
446 sc = m_trackInfoTuple->addItem (
"layer_fly_b", m_ntTrackLayFlyB);
447 sc = m_trackInfoTuple->addItem (
"layer_fly_c", m_ntTrackLayFlyC);
448 sc = m_trackInfoTuple->addItem (
"rec_mode", m_trkRecMode);
449 sc = m_trackInfoTuple->addItem (
"chi2", m_chi2);
450 sc = m_trackInfoTuple->addItem (
"px", m_px);
451 sc = m_trackInfoTuple->addItem (
"py", m_py);
452 sc = m_trackInfoTuple->addItem (
"pz", m_pz);
453 sc = m_trackInfoTuple->addItem (
"pt", m_pt);
454 sc = m_trackInfoTuple->addItem (
"pp", m_pp);
455 sc = m_trackInfoTuple->addItem (
"r", m_r );
456 sc = m_trackInfoTuple->addItem (
"costheta", m_cosTheta);
457 sc = m_trackInfoTuple->addItem (
"theta", m_theta);
458 sc = m_trackInfoTuple->addItem (
"phi", m_phi);
459 sc = m_trackInfoTuple->addItem (
"depth", m_depth);
460 sc = m_trackInfoTuple->addItem (
"br_last_lay", m_brLastLayer);
461 sc = m_trackInfoTuple->addItem (
"ec_last_lay", m_ecLastLayer);
462 sc = m_trackInfoTuple->addItem (
"total_hits", m_totalHits);
463 sc = m_trackInfoTuple->addItem (
"fired_layers", m_totalLayers);
464 sc = m_trackInfoTuple->addItem (
"maxhits_in_layer", m_maxHitsInLayer);
467 log << MSG::ERROR <<
"Cannot book N-tuple:" << long(m_trackInfoTuple) << endmsg;
468 return StatusCode::FAILURE;
473 NTuplePtr nt4(
ntupleSvc,
"FILE450/TrackDiff");
474 if ( nt4 ) { m_trackDiffTuple = nt4; }
477 m_trackDiffTuple =
ntupleSvc->book (
"FILE450/TrackDiff", CLID_RowWiseTuple,
"MucCalibConst N-Tuple");
478 if ( m_trackDiffTuple ) {
479 sc = m_trackDiffTuple->addItem (
"dimu_tag", m_ntDimuTag);
480 sc = m_trackDiffTuple->addItem (
"pos_phi_diff", m_ntPosPhiDiff);
481 sc = m_trackDiffTuple->addItem (
"pos_theta_diff", m_ntPosThetaDiff);
482 sc = m_trackDiffTuple->addItem (
"mom_phi_diff", m_ntMomPhiDiff);
483 sc = m_trackDiffTuple->addItem (
"mom_theta_diff", m_ntMomThetaDiff);
486 log << MSG::ERROR <<
"Cannot book N-tuple:" << long(m_trackDiffTuple) << endmsg;
487 return StatusCode::FAILURE;
492 NTuplePtr nt5(
ntupleSvc,
"FILE450/ClusterSize");
493 if ( nt5 ) { m_clusterSizeTuple = nt5; }
496 m_clusterSizeTuple =
ntupleSvc->book (
"FILE450/ClusterSize", CLID_RowWiseTuple,
"MucCalibConst N-Tuple");
497 if ( m_clusterSizeTuple ) {
498 sc = m_clusterSizeTuple->addItem (
"cluster_size", m_ntClusterSize);
501 log << MSG::ERROR <<
"Cannot book N-tuple:" << long(m_clusterSizeTuple) << endmsg;
502 return StatusCode::FAILURE;
507 NTuplePtr nt6(
ntupleSvc,
"FILE450/EffWindow");
508 if ( nt6 ) { m_effWindowTuple = nt6; }
511 m_effWindowTuple =
ntupleSvc->book (
"FILE450/EffWindow", CLID_RowWiseTuple,
"MucCalibConst N-Tuple");
512 if ( m_effWindowTuple ) {
513 sc = m_effWindowTuple->addItem (
"hit_window", m_ntEffWindow);
516 log << MSG::ERROR <<
"Cannot book N-tuple:" << long(m_effWindowTuple) << endmsg;
517 return StatusCode::FAILURE;
521 NTuplePtr nt7(
ntupleSvc,
"FILE450/ResInfo");
522 if ( nt7 ) { m_resInfoTuple = nt7; }
525 m_resInfoTuple =
ntupleSvc->book (
"FILE450/ResInfo", CLID_RowWiseTuple,
"MucCalibConst N-Tuple");
526 if ( m_resInfoTuple ) {
527 sc = m_resInfoTuple->addItem (
"line_res", m_lineRes);
528 sc = m_resInfoTuple->addItem (
"quad_res", m_quadRes);
529 sc = m_resInfoTuple->addItem (
"extr_res", m_extrRes);
530 sc = m_resInfoTuple->addItem (
"res_part", m_resPart);
531 sc = m_resInfoTuple->addItem (
"res_segment", m_resSegment);
532 sc = m_resInfoTuple->addItem (
"res_layer", m_resLayer);
533 sc = m_resInfoTuple->addItem (
"res_fired", m_resFired);
534 sc = m_resInfoTuple->addItem (
"res_mode", m_resMode);
537 log << MSG::ERROR <<
"Cannot book N-tuple:" << long(m_resInfoTuple) << endmsg;
538 return StatusCode::FAILURE;
572 return StatusCode::SUCCESS;
578 MsgStream log(
msgSvc,
"MucCalibMgr");
579 log << MSG::INFO <<
"Initialize Histograms" << endreq;
598 return StatusCode::SUCCESS;
612 if(i%2 != 0 ) { stripMax = 96*7 + 112; }
613 else { stripMax = 48*8; }
614 sprintf( name,
"HitMapBarrel_Lay_L%d", i );
616 m_hHitMapBarrel_Lay[i] =
new TH1F(name,
title, stripMax, 0, stripMax);
621 sprintf( name,
"HitMapEastEndcap_L%d", i );
623 m_hHitMapEndcap_Lay[0][i] =
new TH1F(name,
title, stripMax, 0, stripMax);
625 sprintf( name,
"HitMapWestEndcap_L%d", i );
627 m_hHitMapEndcap_Lay[1][i] =
new TH1F(name,
title, stripMax, 0, stripMax);
635 sprintf( name,
"HitMapBarrel_Seg_S%d", i );
637 if(i==2 ) { stripMax = 48*5 + 112*4; }
638 else { stripMax = 48*5 + 96*4; }
639 m_hHitMapBarrel_Seg[i] =
new TH1F(name,
title, stripMax, 0, stripMax);
643 sprintf( name,
"HitMapEastEndcap_S%d", i );
644 sprintf(
title,
"Hit map in east-endcap segment %d", i );
646 m_hHitMapEndcap_Seg[0][i] =
new TH1F(name,
title, stripMax, 0, stripMax);
648 sprintf( name,
"HitMapWestEndcap_S%d", i );
649 sprintf(
title,
"Hit map in west-endcap segment %d", i );
650 m_hHitMapEndcap_Seg[1][i] =
new TH1F(name,
title, stripMax, 0, stripMax);
656 m_hHitVsEvent =
new TH1F(
"HitVsEvent",
"Hit VS event",10000,0,10000);
657 m_hHitVsEvent->GetXaxis()->SetTitle(
"Event NO.");
658 m_hHitVsEvent->GetXaxis()->CenterTitle();
659 m_hHitVsEvent->GetYaxis()->SetTitle(
"Hits");
660 m_hHitVsEvent->GetYaxis()->CenterTitle();
662 m_hTrackDistance =
new TH1F(
"TrackDistance",
"Track distance", 1000,-500,500);
663 m_hTrackDistance->GetXaxis()->SetTitle(
"Distance of fit pos and hit pos on 1st layer [cm]");
664 m_hTrackDistance->GetXaxis()->CenterTitle();
666 m_hTrackPosPhiDiff =
new TH1F(
"TrackPosPhiDiff",
"#Delta#phi of tracks pos", 720,-2*
PI,2*
PI);
667 m_hTrackPosPhiDiff->GetXaxis()->SetTitle(
"#Delta#phi [rad]");
668 m_hTrackPosPhiDiff->GetXaxis()->CenterTitle();
670 m_hTrackPosThetaDiff =
new TH1F(
"TrackPosThetaDiff",
"#Delta#theta of tracks pos", 720,-2*
PI,2*
PI);
671 m_hTrackPosThetaDiff->GetXaxis()->SetTitle(
"#Delta#theta [rad]");
672 m_hTrackPosThetaDiff->GetXaxis()->CenterTitle();
674 m_hTrackMomPhiDiff =
new TH1F(
"TrackMomPhiDiff",
"#Delta#phi of tracks mom", 720,-2*
PI,2*
PI);
675 m_hTrackMomPhiDiff->GetXaxis()->SetTitle(
"#Delta#phi [rad]");
676 m_hTrackMomPhiDiff->GetXaxis()->CenterTitle();
678 m_hTrackMomThetaDiff =
new TH1F(
"TrackMomThetaDiff",
"#Delta#theta of tracks mom", 720,-2*
PI,2*
PI);
679 m_hTrackMomThetaDiff->GetXaxis()->SetTitle(
"#Delta#theta [rad]");
680 m_hTrackMomThetaDiff->GetXaxis()->CenterTitle();
683 m_hDimuTracksPosDiff =
new TH2F(
"DimuTracksPosDiff",
"#Delta#phi VS #Delta#theta of dimu tracks pos", 720, -
PI,
PI, 720, -
PI,
PI);
684 m_hDimuTracksPosDiff->GetXaxis()->SetTitle(
"#Delta#theta");
685 m_hDimuTracksPosDiff->GetXaxis()->CenterTitle();
686 m_hDimuTracksPosDiff->GetYaxis()->SetTitle(
"#Delta#phi");
687 m_hDimuTracksPosDiff->GetYaxis()->CenterTitle();
689 m_hDimuTracksMomDiff =
new TH2F(
"DimuTracksMomDiff",
"#Delta#phi VS #Delta#theta of dimu tracks mom", 720, -
PI,
PI, 720, -
PI,
PI);
690 m_hDimuTracksMomDiff->GetXaxis()->SetTitle(
"#Delta#theta");
691 m_hDimuTracksMomDiff->GetXaxis()->CenterTitle();
692 m_hDimuTracksMomDiff->GetYaxis()->SetTitle(
"#Delta#phi");
693 m_hDimuTracksMomDiff->GetYaxis()->CenterTitle();
695 m_hPhiCosTheta =
new TH2F(
"PhiVsCosTheta",
"#phi VS cos(#theta)", 720, -1, 1, 720, -
PI,
PI);
696 m_hPhiCosTheta->GetXaxis()->SetTitle(
"cos(#theta)");
697 m_hPhiCosTheta->GetXaxis()->CenterTitle();
698 m_hPhiCosTheta->GetYaxis()->SetTitle(
"#phi");
699 m_hPhiCosTheta->GetYaxis()->CenterTitle();
701 return StatusCode::SUCCESS;
707 m_hBrLayerFire =
new TH1F(
"BrLayerFire",
"Fires per layer in Barrel", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
708 m_hEcLayerFire =
new TH1F(
"EcLayerFire",
"Fires per layer in Endcap", 2*LAYER_MAX-1, -LAYER_MAX+1, LAYER_MAX-0.5);
710 m_hLayerFire =
new TH1F(
"LayerFire",
"Fires per layer", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
711 m_hLayerExpHit =
new TH1F(
"LayerExpHit",
"Exp hits per layer", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
712 m_hLayerEffHit =
new TH1F(
"LayerEffHit",
"Eff hits per layer", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
713 m_hLayerNosHit =
new TH1F(
"LayerNosHit",
"Nos hits per layer", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
714 m_hLayerEff =
new TH1F(
"LayerEff",
"Layer efficiency", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
715 m_hLayerNosRatio=
new TH1F(
"LayerNosRatio",
"Layer noise hit ratio", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
716 m_hLayerArea =
new TH1F(
"LayerArea",
"Layer area", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
717 m_hLayerNos =
new TH1F(
"LayerNos",
"Layer noise", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
718 m_hLayerCnt =
new TH1F(
"LayerCnt",
"Layer count", LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
720 return StatusCode::SUCCESS;
726 m_hBoxFire =
new TH1F(
"BoxFire",
"Fires per box", BOX_MAX+1, -0.5, BOX_MAX+0.5);
727 m_hBoxExpHit =
new TH1F(
"BoxExpHit",
"Exp hits per box", BOX_MAX+1, -0.5, BOX_MAX+0.5);
728 m_hBoxEffHit =
new TH1F(
"BoxEffHit",
"Eff hits per box", BOX_MAX+1, -0.5, BOX_MAX+0.5);
729 m_hBoxNosHit =
new TH1F(
"BoxNosHit",
"Nos hits per box", BOX_MAX+1, -0.5, BOX_MAX+0.5);
730 m_hBoxEff =
new TH1F(
"BoxEff",
"Box efficiency", BOX_MAX+1, -0.5, BOX_MAX+0.5);
731 m_hBoxNosRatio =
new TH1F(
"BoxNosRatio",
"Box noise hit ratio", BOX_MAX+1, -0.5, BOX_MAX+0.5);
732 m_hBoxArea =
new TH1F(
"BoxArea",
"Box area", BOX_MAX+1, -0.5, BOX_MAX+0.5);
733 m_hBoxNos =
new TH1F(
"BoxNos",
"Box noise", BOX_MAX+1, -0.5, BOX_MAX+0.5);
734 m_hBoxCnt =
new TH1F(
"BoxCnt",
"Box count", BOX_MAX+1, -0.5, BOX_MAX+0.5);
736 return StatusCode::SUCCESS;
744 int part, segment, layer, stripMax;
745 part = segment = layer = stripMax = 0;
747 for(
int i=0; i<BOX_MAX; i++ )
749 m_ptrIdTr->
SetBoxPos( i, &part, &segment, &layer );
750 stripMax = m_ptrIdTr->
GetStripMax( part, segment, layer );
752 sprintf( name,
"StripFireMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
753 sprintf(
title,
"Fires per strip in P%d_S%d_L%d Box%d",part, segment, layer, i );
754 m_hStripFireMap[i] =
new TH1F(name,
title, stripMax, 0, stripMax);
756 sprintf( name,
"StripExpHitMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
757 sprintf(
title,
"Exp hits per strip in P%d_S%d_L%d Box%d",part, segment, layer, i );
758 m_hStripExpHitMap[i] =
new TH1F(name,
title, stripMax, 0, stripMax);
760 sprintf( name,
"StripEffHitMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
761 sprintf(
title,
"Eff hits per strip in P%d_S%d_L%d Box%d",part, segment, layer, i );
762 m_hStripEffHitMap[i] =
new TH1F(name,
title, stripMax, 0, stripMax);
764 sprintf( name,
"StripNosHitMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
765 sprintf(
title,
"Inc hits per strip in P%d_S%d_L%d Box%d",part, segment, layer, i );
766 m_hStripNosHitMap[i] =
new TH1F(name,
title, stripMax, 0, stripMax);
768 sprintf( name,
"StripEffMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
769 sprintf(
title,
"Strip efficiency in P%d_S%d_L%d Box%d",part, segment, layer, i );
770 m_hStripEffMap[i] =
new TH1F(name,
title, stripMax, 0, stripMax);
772 sprintf( name,
"StripNosRatioMap_P%d_S%d_L%d_Box%d", part, segment, layer, i );
773 sprintf(
title,
"Strip noise hit ratio in P%d_S%d_L%d Box%d",part, segment, layer, i );
774 m_hStripNosRatioMap[i] =
new TH1F(name,
title, stripMax, 0, stripMax);
777 m_hStripFire =
new TH1F(
"StripFire",
"Fires per strip",
STRIP_MAX+1, -0.5,
STRIP_MAX+0.5 );
778 m_hStripExpHit =
new TH1F(
"StripExpHit",
"Exp hit per strip",
STRIP_MAX+1, -0.5,
STRIP_MAX+0.5 );
779 m_hStripEffHit =
new TH1F(
"StripEffHit",
"Eff hit per strip",
STRIP_MAX+1, -0.5,
STRIP_MAX+0.5 );
780 m_hStripNosHit =
new TH1F(
"StripNoshit",
"Nos hit per strip",
STRIP_MAX+1, -0.5,
STRIP_MAX+0.5 );
781 m_hStripEff =
new TH1F(
"StripEff",
"Strip efficiency",
STRIP_MAX+1, -0.5,
STRIP_MAX+0.5 );
782 m_hStripNosRatio=
new TH1F(
"StripNosRatio",
"Strip noise hit ratio",
STRIP_MAX+1, -0.5,
STRIP_MAX+0.5 );
783 m_hStripArea =
new TH1F(
"StripArea",
"Strip area",
STRIP_MAX+1, -0.5,
STRIP_MAX+0.5 );
787 return StatusCode::SUCCESS;
794 int stripMax, boxID, stripID, xBin, yBin;
795 stripMax = boxID = stripID = xBin = yBin = 0;
796 double xStart, xEnd, yStart, yEnd, sWidth;
797 xStart = xEnd = yStart = yEnd = sWidth = 0.;
799 m_histArray =
new TObjArray();
809 xStart = -2000, xEnd = 2000;
810 sWidth = B_STR_DST[k];
811 xBin = int((xEnd - xStart)/sWidth);
812 yStart = -B_BOX_WT[k]/2-100, yEnd = B_BOX_WT[k]/2+100;
813 yBin = int((B_BOX_WT[k]+200)/sWidth);
817 xStart = yStart = -1250;
820 xBin = yBin = int((xEnd - xStart)/sWidth);
823 sprintf(name,
"ExpMap2D_P%d_S%d_L%d_Box%d", i, j, k, boxID);
824 m_h2DExpMap[i][j][k] =
new TH2F(name, name, xBin, xStart, xEnd, yBin, yStart, yEnd);
825 sprintf(name,
"HitMap2D_P%d_S%d_L%d_Box%d", i, j, k, boxID);
826 m_h2DHitMap[i][j][k] =
new TH2F(name, name, xBin, xStart, xEnd, yBin, yStart, yEnd);
827 sprintf(name,
"EffMap2D_P%d_S%d_L%d_Box%d", i, j, k, boxID);
828 m_h2DEffMap[i][j][k] =
new TH2F(name, name, xBin, xStart, xEnd, yBin, yStart, yEnd);
832 m_histArray->Add(m_h2DEffMap[i][j][k]);
835 return StatusCode::SUCCESS;
843 int part, segment, layer, stripMax;
844 part = segment = layer = stripMax = 0;
846 for(
int i=0; i<LAYER_MAX; i++ )
848 sprintf( name,
"LayerCluster_L%d", i );
850 m_hLayerCluster[i] =
new TH1F(name,
title, CLUSTER_RANGE, 0, CLUSTER_RANGE );
853 for(
int i=0; i<BOX_MAX; i++ )
855 m_ptrIdTr->
SetBoxPos( i, &part, &segment, &layer );
856 stripMax = m_ptrIdTr->
GetStripMax( part, segment, layer );
857 sprintf( name,
"BoxCluster_P%d_S%d_L%d_Box%d", part, segment, layer, i );
858 sprintf(
title,
"Clusters in P%d_S%d_L%d Box%d", part, segment, layer, i );
859 m_hBoxCluster[i] =
new TH1F(name,
title, CLUSTER_RANGE, 0, CLUSTER_RANGE );
862 m_hLayerClusterCmp =
new TH1F(
"LayerCluster",
"LayerCluster", LAYER_MAX+1, -0.5, LAYER_MAX+0.5 );
863 m_hBoxClusterCmp =
new TH1F(
"BoxCluster",
"BoxCluster", BOX_MAX+1, -0.5, BOX_MAX+0.5 );
865 return StatusCode::SUCCESS;
876 sprintf( name,
"BarrelRes_L%d", i );
878 m_hBarrelResDist[i] =
new TH1F(name,
title, 200, -100, 100 );
883 sprintf( name,
"EndcapRes_L%d", i );
885 m_hEndcapResDist[i] =
new TH1F(name,
title, 200, -100, 100 );
888 m_hBarrelResComp[0] =
new TH1F(
"BarrelResMean",
"BarrelResMean",
B_LAY_NUM+1, -0.5,
B_LAY_NUM+0.5 );
889 m_hBarrelResComp[1] =
new TH1F(
"BarrelResSigma",
"BarrelResSigma",
B_LAY_NUM+1, -0.5,
B_LAY_NUM+0.5 );
890 m_hEndcapResComp[0] =
new TH1F(
"EndcapResMean",
"EndcapResMean",
E_LAY_NUM+1, -0.5,
E_LAY_NUM+0.5 );
891 m_hEndcapResComp[1] =
new TH1F(
"EndcapResSigma",
"EndcapResSigma",
E_LAY_NUM+1, -0.5,
E_LAY_NUM+0.5 );
893 return StatusCode::SUCCESS;
899 MsgStream log(
msgSvc,
"MucCalibMgr");
900 log << MSG::INFO <<
"Initialize Trees" << endreq;
903 m_tJobLog =
new TTree(
"JobLog",
"Job information");
904 m_tJobLog->Branch(
"version", &m_vs,
"m_vs/D");
905 m_tJobLog->Branch(
"high_voltage", &m_hv,
"m_hv/D");
906 m_tJobLog->Branch(
"threshold", &m_th,
"m_th/D");
907 m_tJobLog->Branch(
"event_rate", &m_er,
"m_er/D");
908 m_tJobLog->Branch(
"input_tag", &m_tag,
"m_tag/D");
909 m_tJobLog->Branch(
"rec_mode", &m_recMode,
"m_recMode/I");
910 m_tJobLog->Branch(
"use_pad", &m_usePad,
"m_usePad/I");
911 m_tJobLog->Branch(
"eff_window", &m_effWindow,
"m_effWindow/I");
912 m_tJobLog->Branch(
"cluster_mode", &m_clusterMode,
"m_clusterMode/I");
913 m_tJobLog->Branch(
"check_event", &m_checkEvent,
"m_checkEvent/I");
914 m_tJobLog->Branch(
"dimu_select", &m_dimuSelect,
"m_dimuSelect/I");
915 m_tJobLog->Branch(
"dimu_only", &m_dimuOnly,
"m_dimuOnly/I");
917 m_tJobLog->Branch(
"run_start", &m_fStartRun,
"m_fStartRun/I");
918 m_tJobLog->Branch(
"run_end", &m_fEndRun,
"m_fEndRun/I");
919 m_tJobLog->Branch(
"daq_time", &m_fTotalDAQTime,
"m_fTotalDAQTime/D");
920 m_tJobLog->Branch(
"job_time", &m_fTotalJobTime,
"m_fTotalJobTime/D");
921 m_tJobLog->Branch(
"calib_layer", &m_fCalibLayerNum,
"m_fCalibLayerNum/D");
922 m_tJobLog->Branch(
"calib_box", &m_fCalibBoxNum,
"m_fCalibBoxNum/D");
923 m_tJobLog->Branch(
"calib_strip", &m_fCalibStripNum,
"m_fCalibStripNum/D");
924 m_tJobLog->Branch(
"total_event", &m_fTotalEvent,
"m_fTotalEvent/D");
925 m_tJobLog->Branch(
"total_digi", &m_fTotalDigi,
"m_fTotalDigi/D");
926 m_tJobLog->Branch(
"total_exphit", &m_fTotalExpHit,
"m_fTotalExpHit/D");
927 m_tJobLog->Branch(
"total_effhit", &m_fTotalEffHit,
"m_fTotalEffHit/D");
928 m_tJobLog->Branch(
"total_noshit", &m_fTotalNosHit,
"m_fTotalNosHit/D");
929 m_tJobLog->Branch(
"total_cluster", &m_fTotalClstNum,
"m_fTotalClstNum/D");
930 m_tJobLog->Branch(
"total_strip_area",&m_fTotalStripArea,
"m_fTotalStripArea/D");
931 m_tJobLog->Branch(
"layer_coverage", &m_fLayerCoverage,
"m_fLayerCoverage/D");
932 m_tJobLog->Branch(
"box_coverage", &m_fBoxCoverage,
"m_fBoxCoverage/D");
933 m_tJobLog->Branch(
"strip_coverage", &m_fStripCoverage,
"m_fStripCoverage/D");
936 m_tStatLog =
new TTree(
"StatLog",
"Statistic results");
939 m_tLayConst =
new TTree(
"LayConst",
"layer constants");
940 m_tLayConst->Branch(
"layer_id", &m_fLayerId,
"m_fLayerId/D");
941 m_tLayConst->Branch(
"layer_eff", &m_fLayerEff,
"m_fLayerEff/D");
942 m_tLayConst->Branch(
"layer_efferr", &m_fLayerEffErr,
"m_fLayerEffErr/D");
943 m_tLayConst->Branch(
"layer_nosratio", &m_fLayerNosRatio,
"m_fLayerNosRatio/D");
944 m_tLayConst->Branch(
"layer_digi", &m_fLayerDigi,
"m_fLayerDigi/D");
945 m_tLayConst->Branch(
"layer_noise", &m_fLayerNos,
"m_fLayerNos/D");
946 m_tLayConst->Branch(
"layer_cnt", &m_fLayerCnt,
"m_fLayerCnt/D");
947 m_tLayConst->Branch(
"layer_exphit", &m_fLayerExpHit,
"m_fLayerExpHit/D");
948 m_tLayConst->Branch(
"layer_effHit", &m_fLayerEffHit,
"m_fLayerEffHit/D");
949 m_tLayConst->Branch(
"layer_nosHit", &m_fLayerNosHit,
"m_fLayerNosHit/D");
950 m_tLayConst->Branch(
"layer_cluster", &m_fLayerCluster,
"m_fLayerCluster/D");
951 m_tLayConst->Branch(
"layer_trknum", &m_fLayerTrkNum,
"m_fLayerTrkNum/D");
954 m_tBoxConst =
new TTree(
"BoxConst",
"box constants");
955 m_tBoxConst->Branch(
"box_id", &m_fBoxId,
"m_fBoxId/D");
956 m_tBoxConst->Branch(
"box_part", &m_fBoxPart,
"m_fBoxPart/D");
957 m_tBoxConst->Branch(
"box_segment", &m_fBoxSegment,
"m_fBoxSegment/D");
958 m_tBoxConst->Branch(
"box_layer", &m_fBoxLayer,
"m_fBoxLayer/D");
959 m_tBoxConst->Branch(
"box_eff", &m_fBoxEff,
"m_fBoxEff/D");
960 m_tBoxConst->Branch(
"box_efferr", &m_fBoxEffErr,
"m_fBoxEffErr/D");
961 m_tBoxConst->Branch(
"box_nosratio", &m_fBoxNosRatio,
"m_fBoxNosRatio/D");
962 m_tBoxConst->Branch(
"box_digi", &m_fBoxDigi,
"m_fBoxDigi/D");
963 m_tBoxConst->Branch(
"box_noise", &m_fBoxNos,
"m_fBoxNos/D");
964 m_tBoxConst->Branch(
"box_cnt", &m_fBoxCnt,
"m_fBoxCnt/D");
965 m_tBoxConst->Branch(
"box_exphit", &m_fBoxExpHit,
"m_fBoxExpHit/D");
966 m_tBoxConst->Branch(
"box_effhit", &m_fBoxEffHit,
"m_fBoxEffHit/D");
967 m_tBoxConst->Branch(
"box_noshit", &m_fBoxNosHit,
"m_fBoxNosHit/D");
968 m_tBoxConst->Branch(
"box_cluster", &m_fBoxCluster,
"m_fBoxCluster/D");
969 m_tBoxConst->Branch(
"box_trknum", &m_fBoxTrkNum,
"m_fBoxTrkNum/D");
972 m_tStrConst =
new TTree(
"StrConst",
"strip constants");
973 m_tStrConst->Branch(
"strip_id", &m_fStripId,
"m_fStripId/D");
974 m_tStrConst->Branch(
"strip_part", &m_fStripPart,
"m_fStripPart/D");
975 m_tStrConst->Branch(
"strip_segment", &m_fStripSegment,
"m_fStripSegment/D");
976 m_tStrConst->Branch(
"strip_layer", &m_fStripLayer,
"m_fStripLayer/D");
977 m_tStrConst->Branch(
"strip_eff", &m_fStripEff,
"m_fStripEff/D");
978 m_tStrConst->Branch(
"strip_efferr", &m_fStripEffErr,
"m_fStripEffErr/D");
979 m_tStrConst->Branch(
"strip_nosratio", &m_fStripNosRatio,
"m_fStripNosRatio/D");
980 m_tStrConst->Branch(
"strip_digi", &m_fStripDigi,
"m_fStripDigi/D");
981 m_tStrConst->Branch(
"strip_noise", &m_fStripNos,
"m_fStripNos/D");
982 m_tStrConst->Branch(
"strip_cnt", &m_fStripCnt,
"m_fStripCnt/D");
983 m_tStrConst->Branch(
"strip_effhit", &m_fStripEffHit,
"m_fStripEffHit/D");
984 m_tStrConst->Branch(
"strip_exphit", &m_fStripExpHit,
"m_fStripExpHit/D");
985 m_tStrConst->Branch(
"strip_noshit", &m_fStripNosHit,
"m_fStripNosHit/D");
986 m_tStrConst->Branch(
"strip_trknum", &m_fStripTrkNum,
"m_fStripTrkNum/D");
988 return StatusCode::SUCCESS;
994 int stripMax, boxID, stripID;
995 stripMax = boxID = stripID = 0;
996 double totalArea = 0;
1003 boxID = m_ptrIdTr->
GetBoxId(i, j, k);
1004 m_boxResults[3][boxID] = aBox->
GetArea();
1005 m_layerResults[3][k] += aBox->
GetArea();
1009 for(
int m=0; m<stripMax; m++ )
1013 m_stripResults[3][stripID] = aStrip->
GetArea();
1014 totalArea += m_stripResults[3][stripID];
1019 for(
int i=0; i<LAYER_MAX; i++ ) m_hLayerArea->Fill(i, m_layerResults[3][i]);
1020 for(
int i=0; i<BOX_MAX; i++ ) m_hBoxArea->Fill(i, m_boxResults[3][i]);
1021 for(
int i=0; i<
STRIP_MAX; i++ ) m_hStripArea->Fill(i, m_stripResults[3][i]);
1023 m_fTotalStripArea = totalArea;
1025 return StatusCode::SUCCESS;
1035 MsgStream log(
msgSvc,
"MucCalibMgr");
1036 log << MSG::INFO <<
"Dimu select" << endreq;
1037 Gaudi::svcLocator()->service(
"EventDataSvc",
eventSvc);
1039 if( m_tag > 0) { m_eventTag = (int)m_tag;
return ( StatusCode::SUCCESS ); }
1042 SmartDataPtr<Event::EventHeader> eventHeader(
eventSvc,
"/Event/EventHeader");
1044 log << MSG::FATAL <<
"Could not find event header" << endreq;
1045 return( StatusCode::FAILURE );
1049 SmartDataPtr<RecMdcTrackCol> mdcTrackCol(
eventSvc,
"/Event/Recon/RecMdcTrackCol");
1051 log << MSG::FATAL <<
"Could not find mdc tracks" << endreq;
1052 return( StatusCode::FAILURE);
1054 log << MSG::INFO << mdcTrackCol->size() << endreq;
1055 if( mdcTrackCol->size() != 2 )
return ( StatusCode::FAILURE );
1057 log << MSG::INFO <<
"r1:\t" << (*mdcTrackCol)[0]->r() <<
"\tz1:" << (*mdcTrackCol)[0]->z() << endreq;
1058 log << MSG::INFO <<
"r2:\t" << (*mdcTrackCol)[1]->r() <<
"\tz2:" << (*mdcTrackCol)[1]->z() << endreq;
1059 log << MSG::INFO <<
"p1:\t" << (*mdcTrackCol)[0]->p() <<
"\tp2:" << (*mdcTrackCol)[1]->p() << endreq;
1061 bool mdcFlag1 = (*mdcTrackCol)[0]->charge() + (*mdcTrackCol)[1]->charge() == 0;
1062 bool mdcFlag2 = fabs((*mdcTrackCol)[0]->r())<=1 && fabs((*mdcTrackCol)[0]->z())<3;
1063 bool mdcFlag3 = fabs((*mdcTrackCol)[1]->r())<=1 && fabs((*mdcTrackCol)[1]->z())<3;
1064 bool mdcFlag4 = (*mdcTrackCol)[0]->p()<2 && (*mdcTrackCol)[1]->p()<2;
1065 bool mdcFlag5 = fabs( (*mdcTrackCol)[0]->p() - (*mdcTrackCol)[1]->p() )/( (*mdcTrackCol)[0]->p() + (*mdcTrackCol)[1]->p() ) < 0.5;
1066 bool mdcPass =
false;
1067 if( mdcFlag1 && mdcFlag2 && mdcFlag3 && mdcFlag4 && mdcFlag5 ) mdcPass =
true;
1070 SmartDataPtr<RecTofTrackCol> tofTrackCol(
eventSvc,
"/Event/Recon/RecTofTrackCol");
1072 log << MSG::FATAL <<
"Could not find RecTofTrackCol!" << endreq;
1073 return( StatusCode::FAILURE);
1075 log << MSG::INFO <<
"TOF tracks:\t" << tofTrackCol->size() << endreq;
1080 if(tofTrackCol->size() > 7 && tofTrackCol->size() < 21)
1083 for(
unsigned int itof = 0; itof < tofTrackCol->size(); itof++)
1086 status->
setStatus((*tofTrackCol)[itof]->status());
1088 if( !(status->
is_cluster()) ) {
delete status;
continue; }
1089 if(goodtof==0) t1 = (*tofTrackCol)[itof]->tof();
1090 if(goodtof==1) t2 = (*tofTrackCol)[itof]->tof();
1096 bool tofPass =
false;
1097 if( fabs( t1-t2 ) < 4) tofPass =
true;
1100 SmartDataPtr<RecEmcShowerCol> emcShowerCol(
eventSvc,
"/Event/Recon/RecEmcShowerCol");
1101 if (!emcShowerCol) {
1102 log << MSG::FATAL <<
"Could not find RecEmcShowerCol!" << endreq;
1103 return( StatusCode::FAILURE);
1105 log << MSG::INFO <<
"EMC showers:\t" << emcShowerCol->size() << endreq;
1107 if( emcShowerCol->size() < 2 )
return ( StatusCode::SUCCESS );
1112 e1 = (*emcShowerCol)[0]->energy();
e2 = (*emcShowerCol)[1]->energy();
1113 theta1 = (*emcShowerCol)[0]->theta();
theta2 = (*emcShowerCol)[1]->theta();
1114 phi1 = (*emcShowerCol)[0]->phi();
phi2 = (*emcShowerCol)[1]->phi();
1115 part = (*emcShowerCol)[0]->module();
1117 log << MSG::INFO <<
"e1:\t" <<
e1 <<
"\te2:\t" <<
e2 << endreq;
1118 log << MSG::INFO <<
"theta1:\t" <<
theta1 <<
"\ttheta2:\t" <<
theta2 << endreq;
1119 log << MSG::INFO <<
"phi1:\t" <<
phi1 <<
"\tphi2:\t" <<
phi2 << endreq;
1121 bool emcFlag1 =
e1 < 1.0 &&
e1 > 0.1;
1122 bool emcFlag2 =
e2 < 1.0 &&
e2 > 0.1;
1124 bool emcFlag4 = fabs(
phi1 -
phi2) -
PI < 0.4;
1126 bool emcPass =
false;
1127 if( emcFlag1 && emcFlag2 && emcFlag3 && emcFlag4 ) emcPass =
true;
1130 SmartDataPtr<MucDigiCol> mucDigiCol(
eventSvc,
"/Event/Digi/MucDigiCol");
1132 log << MSG::FATAL <<
"Could not find MUC digi" << endreq;
1133 return( StatusCode::FAILURE);
1135 SmartDataPtr<RecMucTrackCol> mucTrackCol(
eventSvc,
"/Event/Recon/RecMucTrackCol");
1137 log << MSG::FATAL <<
"Could not find RecMucTrackCol" << endreq;
1138 return( StatusCode::FAILURE);
1140 log << MSG::INFO <<
"digi:\t" << mucDigiCol->size() <<
"tracks:\t" << mucTrackCol->size() << endreq;
1142 bool mucFlag1 = mucDigiCol->size()>6;
1143 bool mucFlag2 = mucTrackCol->size()>1;
1144 bool mucPass =
false;
1145 if( mucFlag1 && mucFlag2 ) mucPass =
true;
1147 if( mdcPass && tofPass && emcPass && mucPass ) m_eventTag = 1;
1148 else m_eventTag = (int)m_tag;
1150 return( StatusCode::SUCCESS );
1156 MsgStream log(
msgSvc,
"MucCalibMgr");
1157 log << MSG::INFO <<
"Read event" << endreq;
1159 Gaudi::svcLocator()->service(
"EventDataSvc",
eventSvc);
1160 m_evtBegin = clock();
1163 SmartDataPtr<Event::EventHeader> eventHeader(
eventSvc,
"/Event/EventHeader");
1165 log << MSG::FATAL <<
"Could not find event header" << endreq;
1166 return( StatusCode::FAILURE );
1169 m_currentRun = eventHeader->runNumber();
1170 m_currentEvent = eventHeader->eventNumber();
1171 if( m_fStartRun == 0 ) m_fStartRun = m_currentRun;
1172 m_fEndRun = m_currentRun;
1175 log << MSG::INFO <<
"Run [ " << m_currentRun <<
" ]\tEvent [ " << m_currentEvent <<
" ]" << endreq;
1176 if( ((
long)m_fTotalEvent)%2000 == 0 ) cout << m_fTotalEvent <<
"\tdone!" << endl;
1179 if( m_dimuSelect ) {
1181 log << MSG::INFO <<
"Event tag:\t" << m_eventTag << endreq;
1182 if( m_dimuOnly && m_eventTag != 1 )
return( StatusCode::FAILURE );
1186 log << MSG::INFO <<
"Retrieve digis" << endreq;
1188 SmartDataPtr<MucDigiCol> mucDigiCol(
eventSvc,
"/Event/Digi/MucDigiCol");
1190 log << MSG::FATAL <<
"Could not find MUC digi" << endreq;
1191 return( StatusCode::FAILURE);
1194 int part, segment, layer, strip, pad;
1195 part = segment = layer = strip = pad = 0;
1196 double padX, padY, padZ;
1197 padX = padY = padZ = 0.;
1201 MucDigiCol::iterator digiIter = mucDigiCol->begin();
1203 for (
int digiId =0; digiIter != mucDigiCol->end(); digiIter++, digiId++ )
1205 mucId = (*digiIter)->identify();
1211 log << MSG::DEBUG <<
"[" << part <<
"\t" << segment <<
"\t" << layer <<
"\t" << strip <<
"]\t" ;
1212 if( (digiId+1)%8 == 0 ) log << MSG::DEBUG << endreq;
1216 if(
abs(part)>=
PART_MAX ||
abs(segment)>=SEGMENT_MAX ||
abs(layer)>=LAYER_MAX ||
abs(strip)>=STRIP_INBOX_MAX) {
1217 log << MSG::ERROR << endreq <<
"Digi IDs slop over!" << endreq;
1223 m_digiCol.push_back( aMark );
1224 m_segDigiCol[part][segment].push_back( aMark );
1226 log << MSG::DEBUG << endreq;
1227 log << MSG::INFO <<
"Total digits of this event: " << eventDigi << endreq;
1228 if( eventDigi > 200) log << MSG::ERROR <<
"Event: " << m_currentEvent <<
"\tdigits sharply rise:\t" << eventDigi << endreq;
1246 int clusterNum, bigClusterNum, clusterSize;
1247 clusterNum = bigClusterNum = clusterSize = 0;
1248 if( m_clusterMode ) {
1249 log << MSG::INFO <<
"Searching clusters" << endreq;
1250 m_clusterCol = (*m_ptrMucMark).CreateClusterCol(m_clusterMode, m_digiCol );
1253 for(
unsigned int i=0; i<m_clusterCol.size(); i++ )
1255 clusterSize = m_clusterCol[i].size();
1257 if( clusterSize > CLUSTER_ALARM )
1259 log << MSG::WARNING <<
"Big cluster:" << endreq;
1260 part = (*m_clusterCol[i][0]).Part();
1261 segment = (*m_clusterCol[i][0]).
Segment();
1262 layer = (*m_clusterCol[i][0]).Layer();
1264 if( m_clusterSave ) (*m_fdata) <<
"Event:\t" << m_currentEvent <<
"\tbig cluster " << bigClusterNum << endl;
1266 for(
int j=0; j<clusterSize; j++ )
1268 strip = (*m_clusterCol[i][j]).Strip();
1269 log << MSG::WARNING <<
"[" << part <<
"\t" << segment <<
"\t" << layer <<
"\t" << strip <<
"]\t";
1270 if( (j+1)%8 == 0 ) log << MSG::WARNING << endreq;
1271 if( m_clusterSave ) (*m_fdata) << part <<
"\t" << segment <<
"\t" << layer <<
"\t" << strip << endl;
1273 log << MSG::WARNING << endreq;
1276 else if( clusterSize > 1 )
1278 log << MSG::DEBUG <<
"cluster: " << clusterNum << endreq;
1279 clusterNum ++, m_fTotalClstNum ++;
1280 part = (*m_clusterCol[i][0]).Part();
1281 segment = (*m_clusterCol[i][0]).
Segment();
1282 layer = (*m_clusterCol[i][0]).Layer();
1283 for(
int j=0; j<clusterSize; j++ )
1285 strip = (*m_clusterCol[i][j]).Strip();
1286 log << MSG::DEBUG <<
"[" << part <<
"\t" << segment <<
"\t" << layer <<
"\t" << strip <<
"]\t";
1287 if( (j+1)%8 == 0 ) log << MSG::DEBUG << endreq;
1289 log << MSG::DEBUG << endreq;
1293 if( m_clusterMode) log << MSG::INFO <<
"Total clusters in this event: " << clusterNum << endreq;
1294 else log << MSG::INFO <<
"Clusters not built" << endreq;
1298 log << MSG::INFO <<
"Retrieve tracks" << endreq;
1300 SmartDataPtr<RecMdcTrackCol> mdcTrackCol(
eventSvc,
"/Event/Recon/RecMdcTrackCol");
1302 log << MSG::FATAL <<
"Could not find mdc tracks" << endreq;
1303 return( StatusCode::FAILURE);
1306 RecMdcTrackCol::iterator mdctrkIter = mdcTrackCol->begin();
1307 for (; mdctrkIter != mdcTrackCol->end(); mdctrkIter++)
1309 m_charge = (*mdctrkIter)->charge();
1310 m_mdcpx = (*mdctrkIter)->px();
1311 m_mdcpy = (*mdctrkIter)->py();
1312 m_mdcpz = (*mdctrkIter)->pz();
1313 m_mdcpt = (*mdctrkIter)->pxy();
1314 m_mdcpp = (*mdctrkIter)->p();
1315 m_mdcphi = (*mdctrkIter)->phi();
1316 m_mdctheta = (*mdctrkIter)->theta();
1317 m_mdcTrkInfoTuple->write();
1321 SmartDataPtr<RecMucTrackCol> mucTrackCol(
eventSvc,
"/Event/Recon/RecMucTrackCol");
1323 log << MSG::FATAL <<
"Could not find RecMucTrackCol" << endreq;
1324 return( StatusCode::FAILURE);
1328 if (aRecMucTrackCol->size() < 1) {
1329 log << MSG::INFO <<
"No MUC tracks in this event" << endreq;
1330 return StatusCode::SUCCESS;
1332 log << MSG::INFO <<
"Total tracks of this event: " << aRecMucTrackCol->size() << endreq;
1339 SmartDataPtr<RecEsTimeCol> aRecEsTimeCol(
eventSvc,
"/Event/Recon/RecEsTimeCol");
1340 if( ! aRecEsTimeCol ){
1341 log << MSG::ERROR <<
"Could not find RecEsTimeCol" << endreq;
1342 return StatusCode::FAILURE;
1344 RecEsTimeCol::iterator iter_evt = aRecEsTimeCol->begin();
1347 m_ntEsTime = (*iter_evt)->getStat();
1348 if( (*iter_evt)->getStat() != 211 ) {
1349 log << MSG::WARNING <<
"Event time not by TOF, skip!" << endreq;
1350 return StatusCode::SUCCESS;
1357 phiDiff = thetaDiff = 0.;
1358 if( aRecMucTrackCol->size()==2 && (*aRecMucTrackCol)[0]->GetTotalHits() > 4 && (*aRecMucTrackCol)[1]->GetTotalHits() > 4 )
1361 phi1 = (*aRecMucTrackCol)[0]->getMucPos().phi();
phi2 = (*aRecMucTrackCol)[1]->getMucPos().phi();
1366 theta1 = (*aRecMucTrackCol)[0]->getMucPos().theta();
theta2 = (*aRecMucTrackCol)[1]->getMucPos().theta();
1368 m_hTrackPosPhiDiff->Fill( phiDiff );
1369 m_hTrackPosThetaDiff->Fill( thetaDiff );
1370 m_hDimuTracksPosDiff->Fill( thetaDiff, phiDiff );
1371 m_ntPosPhiDiff = phiDiff;
1372 m_ntPosThetaDiff = thetaDiff;
1374 log << MSG::INFO <<
"PosPhiDiff:\t" << phiDiff <<
"\tPosThetaDiff:\t" << thetaDiff << endreq;
1377 phi1 = (*aRecMucTrackCol)[0]->getMucMomentum().phi();
phi2 = (*aRecMucTrackCol)[1]->getMucMomentum().phi();
1382 theta1 = (*aRecMucTrackCol)[0]->getMucMomentum().theta();
theta2 = (*aRecMucTrackCol)[1]->getMucMomentum().theta();
1385 m_hTrackMomPhiDiff->Fill( phiDiff );
1386 m_hTrackMomThetaDiff->Fill( thetaDiff );
1387 m_hDimuTracksMomDiff->Fill( thetaDiff, phiDiff );
1388 m_ntMomPhiDiff = phiDiff;
1389 m_ntMomThetaDiff = thetaDiff;
1391 log << MSG::INFO <<
"MomPhiDiff:\t" << phiDiff <<
"\tMomThetaDiff:\t" << thetaDiff << endreq;
1392 m_ntDimuTag = m_eventTag;
1393 m_trackDiffTuple->write();
1397 RecMucTrackCol::iterator trackIter = mucTrackCol->begin();
1398 int trackHitNum, rawHitNum, expectedHitNum, segNum, trkRecMode, lastLayerBR, lastLayerEC;
1399 int layerPassNum[3], passMax[TRACK_SEG_MAX][2];
1400 bool firedLay[TRACK_SEG_MAX][LAYER_MAX];
1401 bool seedList[
PART_MAX][LAYER_MAX];
1402 trackHitNum = rawHitNum = expectedHitNum = segNum = trkRecMode = lastLayerBR = lastLayerEC = 0;
1403 layerPassNum[0] = layerPassNum[1] = layerPassNum[2] = 0;
1404 for(
int segi=0; segi<TRACK_SEG_MAX; segi++ ) {
1405 passMax[segi][0] = passMax[segi][1] = 0;
1406 for(
int layi=0; layi<LAYER_MAX; layi++ ) firedLay[segi][layi] = 0;
1413 vector<int> mucRawHitCol;
1414 vector<int> mucExpHitCol;
1416 for (
int trackId = 0; trackIter != mucTrackCol->end(); trackIter++, trackId++)
1421 trackHitNum = (*trackIter)->numHits();
1422 log << MSG::DEBUG <<
"Track: " << trackId <<
" Hits: " << trackHitNum << endreq;
1424 if( trackHitNum == 0 ) {
1425 log << MSG::INFO <<
"Track " << trackId <<
" no hits" << endreq;
1429 m_ntTrackHits = trackHitNum;
1431 m_trkRecMode = trkRecMode = (*trackIter)->GetRecMode();
1432 m_chi2 = (*trackIter)->chi2();
1433 m_px = (*trackIter)->getMucMomentum().x();
1434 m_py = (*trackIter)->getMucMomentum().y();
1435 m_pz = (*trackIter)->getMucMomentum().z();
1436 m_pt = sqrt(m_px*m_px + m_py*m_py);
1437 m_pp = sqrt(m_px*m_px + m_py*m_py + m_pz*m_pz);
1440 m_r = (*trackIter)->getMucPos().mag();
1441 m_cosTheta = (*trackIter)->getMucPos().cosTheta();
1442 m_theta = (*trackIter)->getMucPos().theta();
1443 m_phi = (*trackIter)->getMucPos().phi();
1444 m_depth = (*trackIter)->depth();
1445 m_brLastLayer = lastLayerBR = (*trackIter)->brLastLayer();
1446 m_ecLastLayer = lastLayerEC = (*trackIter)->ecLastLayer();
1447 m_totalHits = (*trackIter)->numHits();
1448 m_totalLayers = (*trackIter)->numLayers();
1449 m_maxHitsInLayer = (*trackIter)->maxHitsInLayer();
1452 m_hPhiCosTheta->Fill(m_cosTheta, m_phi);
1453 log << MSG::INFO <<
"Fill track info" << endreq;
1457 if( m_calHitCol.size() != 0 ) m_calHitCol.clear();
1460 log << MSG::DEBUG <<
"Reconstruction hits(digis in a track): " << endreq;
1462 mucRawHitCol =(*trackIter)->getVecHits();
1463 rawHitNum += mucRawHitCol.size();
1472 for(
unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
1480 mucId = mucRawHitCol[hitId];
1487 log << MSG::DEBUG <<
"[" << part <<
"\t" << segment <<
"\t" << layer <<
"\t" << strip <<
"]\t";
1492 m_calHitCol.push_back( aMark );
1498 if(hitId == 0) { trkSeg[segNum].push_back( aMark ); segNum ++; }
1501 log << MSG::DEBUG <<
"segNum: " << segNum << endreq;
1502 bool notInSeg =
true;
1503 for(
int segi=0; segi<segNum; segi++ )
1507 trkSeg[segi].push_back( aMark );
1513 if( notInSeg ==
true )
1515 trkSeg[segNum].push_back( aMark );
1517 if( segNum > TRACK_SEG_MAX ) {
1518 log << MSG::ERROR <<
"Track segment overflow: " << segNum << endreq;
1524 log << MSG::DEBUG << endreq;
1527 layerPassNum[0] = layerPassNum[1] = layerPassNum[2] = 0;
1528 for(
int segi=0; segi<segNum; segi++ )
1531 passMax[segi][0] = passMax[segi][1] = trkSeg[segi][0]->Layer();
1532 for(
unsigned int hiti=1; hiti<trkSeg[segi].size(); hiti++ )
1534 if( trkSeg[segi][hiti]->Layer() < passMax[segi][0] )
1535 passMax[segi][0] = trkSeg[segi][hiti]->Layer();
1536 if( trkSeg[segi][hiti]->Layer() > passMax[segi][1] )
1537 passMax[segi][1] = trkSeg[segi][hiti]->Layer();
1538 firedLay[segi][trkSeg[segi][hiti]->Layer()] = 1;
1541 for(
int layi=0; layi<LAYER_MAX; layi++ ) {
1542 if( firedLay[segi][layi] ) tmpLayNum ++;
1545 if( segi == 0 ) layerPassNum[0] += passMax[segi][1] + 1;
1546 else layerPassNum[0] += (passMax[segi][1] - passMax[segi][0] + 1);
1548 layerPassNum[1] += (passMax[segi][1] - passMax[segi][0] + 1);
1549 layerPassNum[2] += tmpLayNum;
1551 trkSeg[segi].clear();
1553 m_ntTrackEvent = m_currentEvent;
1554 m_ntTrackTag = m_eventTag;
1555 m_ntTrackSegFly = segNum;
1556 m_ntTrackLayFlyA = layerPassNum[0];
1557 m_ntTrackLayFlyB = layerPassNum[1];
1558 m_ntTrackLayFlyC = layerPassNum[2];
1559 m_trackInfoTuple->write();
1560 log << MSG::INFO <<
"Track\t" << trackId <<
"\tsegment(s):\t" << segNum
1561 <<
"\tlayer passed:\t" << layerPassNum[0] <<
"\t" << layerPassNum[1] <<
"\t" << layerPassNum[2] << endreq;
1566 log << MSG::DEBUG <<
"Fitting hits(expected hits in a track): " << endreq;
1568 mucExpHitCol = (*trackIter)->getExpHits();
1570 expectedHitNum += mucExpHitCol.size();
1572 for(
unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
1578 mucId =mucExpHitCol[hitId];
1601 MucMark* currentMark =
new MucMark( part, segment, layer, strip );
1602 m_expHitCol.push_back( currentMark );
1608 bool isInEffWindow =
false;
1609 isInPos = currentMark->
IsInCol( m_segDigiCol[part][segment] );
1612 if( part ==
BRID && (layer-lastLayerBR>1) )
continue;
1613 if( part !=
BRID && (layer-lastLayerEC>1) )
continue;
1616 if( part==
BRID && layer==0 && (strip<2 || strip>45) )
1620 m_record[part][segment][layer][strip][2] ++;
1621 m_record[part][segment][layer][strip][1] ++;
1622 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
1624 if( m_usePad != 0 ) {
1625 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1626 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
1631 m_record[part][segment][layer][strip][1] ++;
1632 if( m_usePad != 0 ) m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1640 m_record[part][segment][layer][strip][2] ++;
1641 m_record[part][segment][layer][strip][1] ++;
1642 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
1644 if( m_usePad != 0 ) {
1645 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1646 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
1651 else for(
int tempStrip=0, hiti=-m_effWindow; hiti<=m_effWindow; hiti++ )
1653 if( hiti == 0 )
continue;
1654 tempStrip = strip + hiti;
1655 if( tempStrip < 0 || tempStrip > m_ptrIdTr->
GetStripMax(part,segment,layer) )
continue;
1657 isInPos = m_ptrMucMark->
IsInCol( part, segment, layer, tempStrip, m_segDigiCol[part][segment] );
1660 m_record[part][segment][layer][tempStrip][2] ++;
1661 m_record[part][segment][layer][tempStrip][1] ++;
1662 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
1664 if( m_usePad != 0 ) {
1665 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1666 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
1669 m_ntEffWindow = hiti;
1670 m_effWindowTuple->write();
1671 isInEffWindow =
true;
1676 if( isInEffWindow ) {
continue; }
1678 m_record[part][segment][layer][strip][1] ++;
1679 if( m_usePad != 0 ) m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1685 log << MSG::INFO <<
"Fill residual" << endreq;
1686 vector<float> m_lineResCol = (*trackIter)->getDistHits();
1687 vector<float> m_quadResCol = (*trackIter)->getQuadDistHits();
1688 vector<float> m_extrResCol = (*trackIter)->getExtDistHits();
1689 int mode = (*trackIter)->GetRecMode();
1691 for(
unsigned int nres = 0; nres < m_lineResCol.size(); nres++ )
1692 if( fabs(m_lineResCol[nres])>resMax ) resMax = fabs(m_lineResCol[nres]);
1694 log << MSG::INFO <<
"Good track for res" << endreq;
1695 if( trackHitNum > 4 && m_lineResCol[0] != -99)
1698 bool firedFlag[
PART_MAX][LAYER_MAX][2];
1699 for(
int iprt=0; iprt<
PART_MAX; iprt++)
1700 for(
int jlay=0; jlay<LAYER_MAX; jlay++)
1701 firedFlag[iprt][jlay][0] = firedFlag[iprt][jlay][1] =
false;
1703 for(
unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
1708 mucId = mucExpHitCol[hitId];
1714 firedFlag[part][layer][0] =
true;
1717 log << MSG::INFO <<
"Fit res" << endreq;
1718 for(
unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
1723 mucId = mucRawHitCol[hitId];
1728 if( part ==
BRID ) m_hBarrelResDist[layer]->Fill( m_lineResCol[hitId] );
1729 else m_hEndcapResDist[layer]->Fill( m_lineResCol[hitId] );
1732 if( firedFlag[part][layer][0] ==
true && firedFlag[part][layer][1] ==
false )
1735 m_resSegment = segment;
1737 m_lineRes = m_lineResCol[hitId];
1744 m_resInfoTuple->write();
1748 firedFlag[part][layer][1] =
true;
1751 log << MSG::INFO <<
"Exp res" << endreq;
1752 for(
unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
1757 mucId = mucExpHitCol[hitId];
1762 if(firedFlag[part][layer][0] ==
true && firedFlag[part][layer][1] ==
false)
1765 m_resSegment = segment;
1772 m_resInfoTuple->write();
1778 mucRawHitCol.clear();
1779 mucExpHitCol.clear();
1785 m_ntTrackNum = mucTrackCol->size();
1787 m_fTotalEffHit += rawHitNum;
1788 log << MSG::INFO <<
"Total hits in this event, raw: " << rawHitNum <<
"\texpected: " << expectedHitNum << endreq;
1792 log << MSG::INFO <<
"Searching inc/noise hits" << endreq;
1795 for(
unsigned int i=0; i < m_digiCol.size(); i++ )
1799 if( m_digiCol[i]->IsInCol( m_effHitCol ) !=-1)
continue;
1802 for(
unsigned int j=0; j < m_clusterCol.size(); j++ )
1805 for(
unsigned int k=0; k<m_clusterCol[j].size(); k++)
1807 if( m_clusterCol[j][k]->IsInCol(m_effHitCol) != -1)
1814 if ( hasEffHit && (m_digiCol[i]->IsInCol( m_clusterCol[j] ) != -1) ) {
1821 m_nosHitCol.push_back( m_digiCol[i] );
1827 return StatusCode::SUCCESS;
1833 MsgStream log(
msgSvc,
"MucCalibMgr");
1834 log << MSG::INFO <<
"Check event" << endreq;
1838 log << MSG::INFO <<
"Searching over digi" << endreq;
1840 for(
unsigned int i=0; i < m_digiCol.size(); i++ )
1841 for(
unsigned int j=i+1; j < m_digiCol.size(); j++ ) {
1842 if( (*m_digiCol[i]) == (*m_digiCol[j]) ) overDigi ++;
1846 log << MSG::ERROR << overDigi <<
" over digi found in digis!" << endreq;
1850 log << MSG::INFO <<
"Searching over hit" << endreq;
1852 for(
unsigned int i=0; i < m_expHitCol.size(); i++ )
1853 for(
unsigned int j=i+1; j < m_expHitCol.size(); j++ ) {
1854 if( (*m_expHitCol[i]) == (*m_expHitCol[j]) ) overHit ++;
1858 log << MSG::WARNING << overHit <<
" over hits found in rec hits!" << endreq;
1861 log << MSG::INFO <<
"Searching new hit" << endreq;
1864 for(
unsigned int i=0; i < m_expHitCol.size(); i++ )
1866 num = m_expHitCol[i]->NumInCol( m_digiCol );
1869 log << MSG::ERROR <<
"Exp hit not in digis!"
1870 <<
"prt: " << (*m_expHitCol[i]).Part()
1871 <<
"\tseg: " << (*m_expHitCol[i]).
Segment()
1872 <<
"\tlay: " << (*m_expHitCol[i]).Layer()
1873 <<
"\tstr: " << (*m_expHitCol[i]).Strip() << endreq;
1880 log << MSG::WARNING << newHit <<
" new hit(s) found in rec hits!" << endreq;
1882 return StatusCode::SUCCESS;
1888 MsgStream log(
msgSvc,
"MucCalibMgr");
1889 log << MSG::INFO <<
"Fill event" << endreq;
1891 int part, segment, layer, strip, size;
1892 part = segment = layer = strip = size = 0;
1895 log << MSG::INFO <<
"Fill digis" << endreq;
1896 for(
unsigned int i=0; i<m_digiCol.size(); i++ )
1898 part = (*m_digiCol[i]).Part();
1899 segment = (*m_digiCol[i]).
Segment();
1900 layer = (*m_digiCol[i]).Layer();
1901 strip = (*m_digiCol[i]).Strip();
1903 FillDigi( part, segment, layer, strip );
1904 m_record[part][segment][layer][strip][0] ++;
1909 log << MSG::INFO <<
"Fill rec hits" << endreq;
1910 for(
unsigned int i=0; i<m_expHitCol.size(); i++ )
1912 part = (*m_expHitCol[i]).Part();
1913 segment = (*m_expHitCol[i]).
Segment();
1914 layer = (*m_expHitCol[i]).Layer();
1915 strip = (*m_expHitCol[i]).Strip();
1918 m_record[part][segment][layer][strip][4] ++;
1923 log << MSG::INFO <<
"Fill eff hits" << endreq;
1924 for(
unsigned int i=0; i<m_effHitCol.size(); i++ )
1926 part = (*m_effHitCol[i]).Part();
1927 segment = (*m_effHitCol[i]).
Segment();
1928 layer = (*m_effHitCol[i]).Layer();
1929 strip = (*m_effHitCol[i]).Strip();
1936 log << MSG::INFO <<
"Fill clusters" << endreq;
1937 for(
unsigned int i=0; i<m_clusterCol.size(); i++ )
1939 size = m_clusterCol[i].size();
1941 if( size > CLUSTER_CUT )
1943 part = (*m_clusterCol[i][0]).Part();
1944 segment = (*m_clusterCol[i][0]).
Segment();
1945 layer = (*m_clusterCol[i][0]).Layer();
1948 m_ntClusterSize = size;
1949 m_clusterSizeTuple->write();
1954 log << MSG::INFO <<
"Fill inc/noise hits" << endreq;
1955 for(
unsigned int i=0; i<m_nosHitCol.size(); i++ )
1957 part = (*m_nosHitCol[i]).Part();
1958 segment = (*m_nosHitCol[i]).
Segment();
1959 layer = (*m_nosHitCol[i]).Layer();
1960 strip = (*m_nosHitCol[i]).Strip();
1963 m_record[part][segment][layer][strip][3] ++;
1967 m_ntEventId = m_currentEvent;
1968 m_ntEventTag = m_eventTag;
1969 m_ntDigiNum = m_digiCol.size();
1970 m_ntExpHitNum = m_expHitCol.size();
1971 m_ntEffHitNum = m_effHitCol.size();
1972 m_ntNosHitNum = m_nosHitCol.size();
1973 m_ntClusterNum = m_clusterCol.size();
1976 for(
unsigned int i=0; i<m_digiCol.size(); i++ ) {
1977 if( m_digiCol[i] !=
NULL )
delete m_digiCol[i];
1980 for(
unsigned int i=0; i<m_expHitCol.size(); i++ ) {
1981 if( m_expHitCol[i] !=
NULL )
delete m_expHitCol[i];
1989 for(
unsigned int i=0; i<m_calHitCol.size(); i++ ) {
1990 if( m_calHitCol[i] !=
NULL )
delete m_calHitCol[i];
1993 if( m_digiCol.size() != 0 ) m_digiCol.clear();
1994 if( m_expHitCol.size() != 0 ) m_expHitCol.clear();
1995 if( m_calHitCol.size() != 0 ) m_calHitCol.clear();
1996 if( m_effHitCol.size() != 0 ) m_effHitCol.clear();
1997 if( m_nosHitCol.size() != 0 ) m_nosHitCol.clear();
1998 if( m_clusterCol.size() != 0 ) m_clusterCol.clear();
2002 for(
int j=0; j<SEGMENT_MAX; j++ ) {
2003 if( m_segDigiCol[i][j].size() != 0 ) m_segDigiCol[i][j].clear();
2009 m_ntEventTime = (double(m_evtEnd - m_evtBegin))/CLOCKS_PER_SEC;
2010 log << MSG::INFO <<
"Event time:\t" << m_ntEventTime <<
"s" << endreq;
2011 m_eventLogTuple->write();
2013 return StatusCode::SUCCESS;
2021 if( (
int)m_tag || m_dimuOnly==0 || (m_dimuOnly==1&&m_eventTag==1) )
2026 if( layer%2 == 0 ) {
2027 if( segment > 4 ) tmpId = segment*48 + (47 - strip);
2028 else tmpId = segment*48 + strip;
2030 else if( segment < 3 ) tmpId = segment*96 + strip;
2031 else tmpId = (segment-1)*96 + 112 + strip;
2033 m_hHitMapBarrel_Lay[layer]->Fill(tmpId);
2036 if( segment != B_TOP ) {
2038 tmpId = 48*((layer+1)/2) + 96*(layer/2) + ( ((layer+1)%2)*48 + (layer%2)*96 - strip - 1 );
2039 else tmpId = 48*((layer+1)/2) + 96*(layer/2) + strip;
2041 else tmpId = 48*((layer+1)/2) + 112*(layer/2) + strip;
2043 m_hHitMapBarrel_Seg[segment]->Fill(tmpId);
2045 else if( part ==
EEID )
2048 tmpId = segment*64 + strip;
2049 m_hHitMapEndcap_Lay[0][layer]->Fill(tmpId);
2051 tmpId = layer*64 + strip;
2052 m_hHitMapEndcap_Seg[0][segment]->Fill(tmpId);
2054 else if( part == EWID )
2057 tmpId = segment*64 + strip;
2058 m_hHitMapEndcap_Lay[1][layer]->Fill(tmpId);
2060 tmpId = layer*64 + strip;
2061 m_hHitMapEndcap_Seg[1][segment]->Fill(tmpId);
2066 int boxId = m_ptrIdTr->
GetBoxId( part, segment, layer );
2067 int strId = m_ptrIdTr->
GetStripId( part, segment, layer, strip );
2069 m_hStripFireMap[boxId]->Fill( strip );
2070 m_hStripFire->Fill( strId );
2071 m_hBoxFire->Fill( boxId );
2072 m_hLayerFire->Fill( layer );
2073 if( part==
BRID ) m_hBrLayerFire->Fill( layer );
2074 else if( part==
EEID ) m_hEcLayerFire->Fill( layer+1 );
2075 else m_hEcLayerFire->Fill( -(layer+1) );
2077 return StatusCode::SUCCESS;
2082 int boxId = m_ptrIdTr->
GetBoxId( part, segment, layer );
2083 int strId = m_ptrIdTr->
GetStripId( part, segment, layer, strip );
2085 m_hStripExpHitMap[boxId]->Fill( strip );
2086 m_hStripExpHit->Fill( strId );
2087 m_hBoxExpHit->Fill( boxId );
2088 m_hLayerExpHit->Fill( layer );
2090 return StatusCode::SUCCESS;
2095 int boxId = m_ptrIdTr->
GetBoxId( part, segment, layer );
2096 int strId = m_ptrIdTr->
GetStripId( part, segment, layer, strip );
2098 m_hStripEffHitMap[boxId]->Fill( strip );
2099 m_hStripEffHit->Fill( strId );
2100 m_hBoxEffHit->Fill( boxId );
2101 m_hLayerEffHit->Fill( layer );
2103 return StatusCode::SUCCESS;
2108 int boxId = m_ptrIdTr->
GetBoxId( part, segment, layer );
2109 int strId = m_ptrIdTr->
GetStripId( part, segment, layer, strip );
2111 m_hStripNosHitMap[boxId]->Fill( strip );
2112 m_hStripNosHit->Fill( strId );
2113 m_hBoxNosHit->Fill( boxId );
2114 m_hLayerNosHit->Fill( layer );
2116 return StatusCode::SUCCESS;
2121 int boxId = m_ptrIdTr->
GetBoxId( part, segment, layer );
2123 m_hLayerCluster[layer]->Fill( size );
2124 m_hBoxCluster[boxId]->Fill( size );
2126 return StatusCode::SUCCESS;
2136 MsgStream log(
msgSvc,
"MucCalibMgr");
2137 log<< MSG::INFO <<
"Start efficiency and noise calibration" << endreq;
2140 m_fTotalDAQTime = m_fTotalEvent/m_er;
2141 log<< MSG::INFO <<
"Total run time:\t" << m_fTotalDAQTime << endreq;
2147 if( m_usePad != 0 )
PadEff();
2149 return StatusCode::SUCCESS;
2154 MsgStream log(
msgSvc,
"MucCalibMgr");
2155 log<< MSG::INFO <<
"Analyse layer efficiency and noise" << endreq;
2157 int part, segment, layer, stripMax;
2158 part = segment = layer = stripMax = 0;
2160 double digi, effHit, trackNum, nosHit, recHit;
2161 double eff, effErr, noise, nosRatio, nosRatioErr, cnt, cluster;
2162 eff = effErr = noise = nosRatio = nosRatioErr = cnt = cluster = 0.;
2164 for(
int i=0; i<LAYER_MAX; i++ )
2166 digi = effHit = trackNum = nosHit = recHit = 0;
2170 for(
int k=0; k<SEGMENT_MAX; k++)
2173 for(
int n=0;
n<stripMax;
n++ )
2175 digi += m_record[j][k][i][
n][0];
2176 trackNum += m_record[j][k][i][
n][1];
2177 effHit += m_record[j][k][i][
n][2];
2178 nosHit += m_record[j][k][i][
n][3];
2179 recHit += m_record[j][k][i][
n][4];
2185 if( trackNum == 0 ) {
2192 eff = ( (double)effHit )/trackNum;
2193 effErr = sqrt( eff*(1-eff)/trackNum );
2194 m_fCalibLayerNum ++;
2198 if( m_layerResults[3][i] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
2199 noise = DEFAULT_NOS_VALUE;
2201 if( m_recMode == 2 ) noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW*m_layerResults[3][i]);
2202 else noise = (double)nosHit/(m_fTotalDAQTime * m_layerResults[3][i]);
2206 nosRatio = ( (double)nosHit )/digi;
2207 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
2210 nosRatio = DEFAULT_INC_VALUE;
2215 if( m_recMode == 2 )
2216 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_layerResults[3][i]);
2218 cnt = (double)digi/(m_fTotalDAQTime * m_layerResults[3][i]);
2221 cluster = m_hLayerCluster[i]->GetMean();
2223 m_layerResults[0][ i ] = eff;
2224 m_layerResults[1][ i ] = effErr;
2225 m_layerResults[2][ i ] = noise;
2226 m_layerResults[4][ i ] = cluster;
2227 m_layerResults[5][ i ] = trackNum;
2230 m_hLayerEff->Fill( i, eff );
2231 m_hLayerEff->SetBinError( i+1, effErr );
2232 m_hLayerNosRatio->Fill( i, nosRatio );
2233 m_hLayerNosRatio->SetBinError( i+1, nosRatioErr );
2234 m_hLayerNos->Fill( i, noise );
2235 m_hLayerCnt->Fill( i, cnt );
2238 m_hLayerClusterCmp->Fill( i, cluster );
2243 m_fLayerEffErr = effErr;
2244 m_fLayerTrkNum = trackNum;
2245 m_fLayerExpHit = recHit;
2246 m_fLayerEffHit = effHit;
2247 m_fLayerNosRatio = nosRatio;
2248 m_fLayerDigi = digi;
2249 m_fLayerNosHit = nosHit;
2250 m_fLayerNos = noise;
2252 m_tLayConst->Fill();
2255 m_fLayerCluster = cluster;
2259 return StatusCode::SUCCESS;
2264 MsgStream log(
msgSvc,
"MucCalibMgr");
2265 log<< MSG::INFO <<
"Analyse box efficiency and noise" << endreq;
2267 int part, segment, layer, stripMax;
2268 part = segment = layer = stripMax = 0;
2270 double digi, effHit, trackNum, nosHit, recHit;
2271 double eff, effErr, noise, nosRatio, nosRatioErr, cnt, cluster;
2272 eff = effErr = noise = nosRatio = nosRatioErr = cnt = cluster = 0.;
2274 for(
int i=0; i<BOX_MAX; i++ )
2276 m_ptrIdTr->
SetBoxPos( i, &part, &segment, &layer );
2277 stripMax = m_ptrIdTr->
GetStripMax( part, segment, layer );
2279 digi = effHit = trackNum = nosHit = recHit = 0;
2280 for(
int j=0; j<stripMax; j++ )
2282 digi += m_record[part][segment][layer][j][0];
2283 trackNum += m_record[part][segment][layer][j][1];
2284 effHit += m_record[part][segment][layer][j][2];
2285 nosHit += m_record[part][segment][layer][j][3];
2286 recHit += m_record[part][segment][layer][j][4];
2290 if( trackNum == 0 ) {
2297 eff = ( (double)effHit )/trackNum;
2298 effErr = sqrt( eff*(1-eff)/trackNum );
2303 if( m_boxResults[3][i] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
2304 noise = DEFAULT_NOS_VALUE;
2306 if( m_recMode == 2 )
2307 noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW * m_boxResults[3][i]);
2309 noise = (double)nosHit/(m_fTotalDAQTime * m_boxResults[3][i]);
2313 nosRatio = ( (double)nosHit )/digi;
2314 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
2317 nosRatio = DEFAULT_INC_VALUE;
2322 if( m_recMode == 2 )
2323 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_boxResults[3][i]);
2325 cnt = (double)digi/(m_fTotalDAQTime * m_boxResults[3][i]);
2328 cluster = m_hBoxCluster[i]->GetMean();
2330 m_boxResults[0][ i ] = eff;
2331 m_boxResults[1][ i ] = effErr;
2332 m_boxResults[2][ i ] = noise;
2333 m_boxResults[4][ i ] = cluster;
2334 m_boxResults[5][ i ] = trackNum;
2337 m_hBoxEff->Fill( i, eff );
2338 m_hBoxEff->SetBinError( i+1, effErr );
2339 m_hBoxNosRatio->Fill( i, nosRatio );
2340 m_hBoxNosRatio->SetBinError( i+1, nosRatioErr );
2341 m_hBoxNos->Fill( i, noise );
2342 m_hBoxCnt->Fill( i, cnt );
2345 m_hBoxClusterCmp->Fill( i, cluster );
2350 m_fBoxSegment = segment;
2351 m_fBoxLayer = layer;
2353 m_fBoxEffErr = effErr;
2354 m_fBoxTrkNum = trackNum;
2355 m_fBoxExpHit = recHit;
2356 m_fBoxEffHit = effHit;
2357 m_fBoxNosRatio = nosRatio;
2359 m_fBoxNosHit = nosHit;
2362 m_tBoxConst->Fill();
2365 m_fBoxCluster = cluster;
2369 return StatusCode::SUCCESS;
2374 MsgStream log(
msgSvc,
"MucCalibMgr");
2375 log<< MSG::INFO <<
"Analyse strip efficiency and noise" << endreq;
2377 int part, segment, layer, stripMax;
2378 part = segment = layer = stripMax = 0;
2380 double digi, effHit, trackNum, nosHit, recHit;
2381 double eff, effErr, noise, nosRatio, nosRatioErr, cnt, cluster;
2382 eff = effErr = noise = nosRatio = nosRatioErr = cnt = cluster = 0.;
2384 for(
int i=0; i<BOX_MAX; i++ )
2386 m_ptrIdTr->
SetBoxPos( i, &part, &segment, &layer );
2387 stripMax = m_ptrIdTr->
GetStripMax( part, segment, layer );
2389 for(
int j=0; j<stripMax; j++ )
2391 digi = m_record[part][segment][layer][j][0];
2392 trackNum = m_record[part][segment][layer][j][1];
2393 effHit = m_record[part][segment][layer][j][2];
2394 nosHit = m_record[part][segment][layer][j][3];
2395 recHit = m_record[part][segment][layer][j][4];
2397 int stripId = m_ptrIdTr->
GetStripId( part, segment, layer, j );
2400 if( trackNum == 0 ) {
2407 eff = ( (double)effHit )/trackNum;
2408 effErr = sqrt( eff*(1-eff)/trackNum );
2409 m_fCalibStripNum ++;
2413 if( m_stripResults[3][stripId] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
2414 noise = DEFAULT_NOS_VALUE;
2416 if( m_recMode == 2 )
2417 noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW * m_stripResults[3][stripId]);
2419 noise = (double)nosHit/(m_fTotalDAQTime * m_stripResults[3][stripId]);
2423 nosRatio = ( (double)nosHit )/digi;
2424 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
2427 nosRatio = DEFAULT_INC_VALUE;
2432 if( m_recMode == 2 )
2433 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_stripResults[3][stripId]);
2435 cnt = (double)digi/(m_fTotalDAQTime * m_stripResults[3][stripId]);
2438 m_stripResults[0][ stripId ] = eff;
2439 m_stripResults[1][ stripId ] = effErr;
2440 m_stripResults[2][ stripId ] = noise;
2441 m_stripResults[5][ stripId ] = trackNum;
2444 m_hStripEffMap[i]->Fill( j, eff );
2445 m_hStripEffMap[i]->SetBinError( j+1, effErr );
2446 m_hStripEff->Fill( stripId, eff );
2447 m_hStripEff->SetBinError( stripId+1, effErr );
2448 m_hStripNosRatioMap[i]->Fill( j, nosRatio );
2449 m_hStripNosRatioMap[i]->SetBinError( j+1, nosRatioErr );
2450 m_hStripNosRatio->Fill( stripId, nosRatio );
2451 m_hStripNosRatio->SetBinError( stripId+1, nosRatioErr );
2452 m_hStripNos->Fill( stripId, noise );
2453 m_hStripCnt->Fill( stripId, cnt );
2456 m_fStripId = stripId;
2457 m_fStripPart = part;
2458 m_fStripSegment = segment;
2459 m_fStripLayer = layer;
2461 m_fStripEffErr = effErr;
2462 m_fStripNosRatio = nosRatio;
2463 m_fStripDigi = digi;
2464 m_fStripNos = noise;
2466 m_fStripEffHit = effHit;
2467 m_fStripExpHit = recHit;
2468 m_fStripNosHit = nosHit;
2469 m_fStripTrkNum = trackNum;
2470 m_tStrConst->Fill();
2475 return StatusCode::SUCCESS;
2480 MsgStream log(
msgSvc,
"MucCalibMgr");
2482 int xBinStart, xBinEnd, yBinStart, yBinEnd;
2483 double expHit, firedHit, eff;
2484 eff = expHit = firedHit = 0.;
2490 xBinStart = m_h2DExpMap[i][j][k]->GetXaxis()->GetFirst();
2491 xBinEnd = m_h2DExpMap[i][j][k]->GetXaxis()->GetLast() + 1;
2492 yBinStart = m_h2DExpMap[i][j][k]->GetYaxis()->GetFirst();
2493 yBinEnd = m_h2DExpMap[i][j][k]->GetYaxis()->GetLast() + 1;
2495 for(
int xi = xBinStart; xi<xBinEnd; xi++ )
2496 for(
int yi = yBinStart; yi<yBinEnd; yi++ )
2498 expHit = m_h2DExpMap[i][j][k]->GetBinContent(xi, yi);
2499 firedHit = m_h2DHitMap[i][j][k]->GetBinContent(xi, yi);
2501 if( expHit !=0 ) eff = firedHit / expHit;
2505 cout<<
"Eff error:\t["<<i<<
"\t"<<j<<
"\t"<<k<<
",\t"<<xi<<
"\t"<<yi<<
"]:\t"
2506 <<eff<<
"\t,"<<firedHit<<
" "<<expHit<<endl;
2508 m_h2DEffMap[i][j][k]->SetBinContent(xi, yi, eff);
2511 return StatusCode::SUCCESS;
2516 MsgStream log(
msgSvc,
"MucCalibMgr");
2518 return StatusCode::SUCCESS;
2523 MsgStream log(
msgSvc,
"MucCalibMgr");
2524 log << MSG::INFO <<
"Analyse spacial resolution" << endreq;
2526 double resMean, resMeanErr, resSigma, resSigmaErr;
2527 resMean = resMeanErr = resSigma = resSigmaErr = 0.;
2531 m_hBarrelResDist[i]->Fit(
"gaus");
2532 resMean = m_hBarrelResDist[i]->GetFunction(
"gaus")->GetParameter(
"Mean");
2533 resSigma = m_hBarrelResDist[i]->GetFunction(
"gaus")->GetParameter(
"Sigma");
2534 resMeanErr = m_hBarrelResDist[i]->GetFunction(
"gaus")->GetParError(1);
2535 resSigmaErr = m_hBarrelResDist[i]->GetFunction(
"gaus")->GetParError(2);
2537 m_hBarrelResComp[0]->Fill( i, resMean );
2538 m_hBarrelResComp[0]->SetBinError( i+1, resMeanErr );
2539 m_hBarrelResComp[1]->Fill( i, resSigma );
2540 m_hBarrelResComp[1]->SetBinError( i+1, resSigmaErr );
2545 m_hEndcapResDist[i]->Fit(
"gaus");
2546 resMean = m_hEndcapResDist[i]->GetFunction(
"gaus")->GetParameter(
"Mean");
2547 resSigma = m_hEndcapResDist[i]->GetFunction(
"gaus")->GetParameter(
"Sigma");
2548 resMeanErr = m_hEndcapResDist[i]->GetFunction(
"gaus")->GetParError(1);
2549 resSigmaErr = m_hEndcapResDist[i]->GetFunction(
"gaus")->GetParError(2);
2551 m_hEndcapResComp[0]->Fill( i, resMean );
2552 m_hEndcapResComp[0]->SetBinError( i+1, resMeanErr );
2553 m_hEndcapResComp[1]->Fill( i, resSigma );
2554 m_hEndcapResComp[1]->SetBinError( i+1, resSigmaErr );
2557 return StatusCode::SUCCESS;
2562 MsgStream log(
msgSvc,
"MucCalibMgr");
2563 log << MSG::INFO <<
"Save calibration constants" << endreq;
2567 double layerXY[2][LAYER_MAX];
2568 double layerEXY[2][LAYER_MAX];
2569 for(
int i=0; i<LAYER_MAX; i++ )
2573 if( m_layerResults[5][i] >= 100*TRACK_THRESHOLD ) {
2574 layerXY[1][i] = m_layerResults[0][i];
2575 layerEXY[1][i] = m_layerResults[1][i];
2582 m_geLayerEff =
new TGraphErrors(LAYER_MAX, layerXY[0], layerXY[1], layerEXY[0], layerEXY[1]);
2583 m_geLayerEff->SetMarkerStyle(25);
2584 m_geLayerEff->SetMarkerSize(0.5);
2585 m_cv[0] =
new TCanvas(
"GoodLayerEff",
"Layer efficiency", 50, 50, 800, 600);
2586 m_cv[0]->SetFillColor(0);
2587 m_cv[0]->SetBorderMode(0);
2588 m_geLayerEff->Draw(
"AP");
2592 double boxXY[2][BOX_MAX];
2593 double boxEXY[2][BOX_MAX];
2594 for(
int i=0; i<BOX_MAX; i++ )
2598 if( m_boxResults[5][i] >= 10*TRACK_THRESHOLD ) {
2599 boxXY[1][i] = m_boxResults[0][i];
2600 boxEXY[1][i] = m_boxResults[1][i];
2607 m_geBoxEff =
new TGraphErrors(BOX_MAX, boxXY[0], boxXY[1], boxEXY[0], boxEXY[1]);
2608 m_geBoxEff->SetMarkerStyle(25);
2609 m_geBoxEff->SetMarkerSize(0.5);
2610 m_cv[1] =
new TCanvas(
"GoodBoxEff",
"Box efficiency", 75, 75, 800, 600);
2611 m_cv[1]->SetFillColor(0);
2612 m_cv[1]->SetBorderMode(0);
2613 m_geBoxEff->Draw(
"AP");
2623 if( m_stripResults[5][i] >= TRACK_THRESHOLD ) {
2624 stripXY[1][i] = m_stripResults[0][i];
2625 stripEXY[1][i] = m_stripResults[1][i];
2632 m_geStripEff =
new TGraphErrors(
STRIP_MAX, stripXY[0], stripXY[1], stripEXY[0], stripEXY[1]);
2633 m_geStripEff->SetMarkerStyle(25);
2634 m_geStripEff->SetMarkerSize(0.5);
2635 m_cv[2] =
new TCanvas(
"GoodStripEff",
"Strip efficiency", 100, 100, 800, 600);
2636 m_cv[2]->SetFillColor(0);
2637 m_cv[2]->SetBorderMode(0);
2638 m_geStripEff->Draw(
"AP");
2644 m_hHitMapBarrel_Lay[i]->Write();
2647 m_hHitMapEndcap_Lay[0][i]->Write();
2648 m_hHitMapEndcap_Lay[1][i]->Write();
2654 m_hHitMapBarrel_Seg[i]->Write();
2657 m_hHitMapEndcap_Seg[0][i]->Write();
2658 m_hHitMapEndcap_Seg[1][i]->Write();
2661 m_hTrackDistance->Fit(
"gaus");
2662 m_hTrackPosPhiDiff->Fit(
"gaus");
2663 m_hTrackPosThetaDiff->Fit(
"gaus");
2664 m_hTrackMomPhiDiff->Fit(
"gaus");
2665 m_hTrackMomThetaDiff->Fit(
"gaus");
2667 m_hTrackDistance->Write();
2668 m_hTrackPosPhiDiff->Write();
2669 m_hTrackPosThetaDiff->Write();
2670 m_hTrackMomPhiDiff->Write();
2671 m_hTrackMomThetaDiff->Write();
2676 m_hBarrelResComp[0]->Write();
2677 m_hBarrelResComp[1]->Write();
2678 m_hEndcapResComp[0]->Write();
2679 m_hEndcapResComp[1]->Write();
2681 if( m_usePad != 0 ) m_histArray->Write();
2683 for(
int i=0; i<BOX_MAX; i++ )
2685 m_hStripFireMap[ i ]->Write();
2686 m_hStripExpHitMap[ i ]->Write();
2687 m_hStripEffHitMap[ i ]->Write();
2688 m_hStripNosHitMap[ i ]->Write();
2689 m_hStripEffMap[ i ]->Write();
2690 m_hStripNosRatioMap[ i ]->Write();
2692 m_hStripFire->Write();
2693 m_hStripExpHit->Write();
2694 m_hStripEffHit->Write();
2695 m_hStripNosHit->Write();
2696 m_hStripEff->Write();
2697 m_hStripArea->Write();
2698 m_hStripNos->Write();
2699 m_hStripNosRatio->Write();
2700 log << MSG::INFO <<
"Save LV2 histograms done!" << endreq;
2702 m_hBoxFire->Write();
2703 m_hBoxExpHit->Write();
2704 m_hBoxEffHit->Write();
2705 m_hBoxNosHit->Write();
2707 m_hBoxArea->Write();
2709 m_hBoxNosRatio->Write();
2710 log << MSG::INFO <<
"Save LV1 histograms done!" << endreq;
2712 m_hBrLayerFire->Write();
2713 m_hEcLayerFire->Write();
2714 m_hLayerFire->Write();
2715 m_hLayerExpHit->Write();
2716 m_hLayerEffHit->Write();
2717 m_hLayerNosHit->Write();
2718 m_hLayerEff->Write();
2719 m_hLayerArea->Write();
2720 m_hLayerNos->Write();
2721 m_hLayerNosRatio->Write();
2723 for(
int i=0; i<LAYER_MAX; i++ ) m_hLayerCluster[i]->
Write();
2724 for(
int i=0; i<BOX_MAX; i++ ) m_hBoxCluster[i]->
Write();
2725 m_hLayerClusterCmp->Write();
2726 m_hBoxClusterCmp->Write();
2728 log << MSG::INFO <<
"Save histograms done!" << endreq;
2731 m_fLayerCoverage = 100*(double)m_fCalibLayerNum/LAYER_MAX;
2732 m_fBoxCoverage = 100*(double)m_fCalibBoxNum/BOX_MAX;
2733 m_fStripCoverage = 100*(double)m_fCalibStripNum/
STRIP_MAX;
2735 long digi_num, trk_num, eff_hit, nos_hit, exp_hit;
2736 m_tStatLog->Branch(
"digi_num", &digi_num,
"digi_num/I");
2737 m_tStatLog->Branch(
"trk_num", &trk_num,
"trk_num/I");
2738 m_tStatLog->Branch(
"eff_hit", &eff_hit,
"eff_hit/I");
2739 m_tStatLog->Branch(
"nos_hit", &nos_hit,
"nos_hit/I");
2740 m_tStatLog->Branch(
"exp_hit", &exp_hit,
"exp_hit/I");
2748 for(
int n=0;
n<stripMax;
n++ )
2750 digi_num = m_record[i][j][k][
n][0];
2751 trk_num = m_record[i][j][k][
n][1];
2752 eff_hit = m_record[i][j][k][
n][2];
2753 nos_hit = m_record[i][j][k][
n][3];
2754 exp_hit = m_record[i][j][k][
n][4];
2759 m_jobFinish = clock();
2760 m_fTotalJobTime = (double)(m_jobFinish - m_jobStart)/CLOCKS_PER_SEC;
2764 m_tStatLog->Write();
2765 m_tLayConst->Write();
2766 m_tBoxConst->Write();
2767 m_tStrConst->Write();
2770 if( m_fdata !=
NULL ) m_fdata->close();
2772 return StatusCode::SUCCESS;
2777 MsgStream log(
msgSvc,
"MucCalibMgr");
2778 log << MSG::INFO << endreq;
2779 log << MSG::INFO <<
"Total events : " << m_fTotalEvent << endreq;
2780 log << MSG::INFO <<
"Total digis : " << m_fTotalDigi << endreq;
2781 log << MSG::INFO <<
"Total rec hits : " << m_fTotalExpHit << endreq;
2782 log << MSG::INFO <<
"Total eff hits : " << m_fTotalEffHit << endreq;
2783 log << MSG::INFO <<
"Total inc hits : " << m_fTotalNosHit << endreq;
2784 log << MSG::INFO <<
"Total clusters : " << m_fTotalClstNum<< endreq;
2786 log << MSG::INFO <<
"Strip calibrated percentage: "
2787 << 100*(double)m_fCalibStripNum/
STRIP_MAX <<
"%" << endreq;
2788 log << MSG::INFO <<
"Box calibrated percentage: "
2789 << 100*(double)m_fCalibBoxNum/BOX_MAX <<
"%" << endreq;
2790 log << MSG::INFO <<
"Layer calibrated percentage: "
2791 << 100*(double)m_fCalibLayerNum/LAYER_MAX <<
"%" << endreq;
2793 log << MSG::INFO <<
"Calibration run successfully" << endreq << endreq;
2795 return StatusCode::SUCCESS;
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
vector< MucMark * > mark_col
ObjectVector< RecMucTrack > RecMucTrackCol
StatusCode ClearHistoLV2()
StatusCode InitResHisto()
StatusCode FillEffHit(int part, int segment, int layer, int strip)
StatusCode InitOnlineHisto()
StatusCode Init2DEffHisto()
StatusCode FillExpHit(int part, int segment, int layer, int strip)
StatusCode ClearResHisto()
StatusCode ClearHistoLV0()
StatusCode ClearClusterHisto()
StatusCode AnalyseCluster()
StatusCode ClearHistoLV1()
StatusCode InitHistoLV1()
StatusCode InitHistoLV2()
StatusCode InitClusterHisto()
StatusCode FillCluster(int part, int segment, int layer, int size)
StatusCode ClearOnlineHisto()
StatusCode EffAndNoiseLV2()
MucCalibMgr(std::vector< double > jobInfo, std::vector< int > configInfo, std::string outputFile)
StatusCode InitHistoLV0()
IDataProviderSvc * eventSvc
StatusCode EffAndNoiseLV1()
StatusCode FillNosHit(int part, int segment, int layer, int strip)
StatusCode EffAndNoiseLV0()
StatusCode InitConstTree()
StatusCode ClearConstTree()
StatusCode AnalyseEffAndNoise()
StatusCode Clear2DEffHisto()
StatusCode FillDigi(int part, int segment, int layer, int strip)
static int barrel_ec(const Identifier &id)
Values of different levels.
static int layer(const Identifier &id)
static int channel(const Identifier &id)
static int segment(const Identifier &id)
bool IsInSegWith(MucMark &other)
int IsInCol(int part, int segment, int layer, int strip, mark_col &aCol)
void setStatus(unsigned int status)