BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
MucCalibMgr.cxx
Go to the documentation of this file.
1//------------------------------------------------------------------------------|
2// [File ]: MucCalibMgr.cxx |
3// [Brief ]: Manager of MUC calibration algrithom |
4// [Author]: Xie Yuguang, <[email protected]> |
5// [Date ]: Mar 28, 2006 |
6// [Log ]: See ChangLog |
7//------------------------------------------------------------------------------|
8#include <cstdio>
9#include <iostream>
10#include <fstream>
11
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;
20
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"
33
35#include "Identifier/MucID.h"
39
40#include "McTruth/McKine.h"
41#include "McTruth/McParticle.h"
42#include "McTruth/MucMcHit.h"
43
48
49//#include "ReconEvent/ReconEvent.h"
50#include "MucRawEvent/MucDigi.h"
54
57#include "MucCalibAlg/MucMark.h"
60
61using namespace std;
62using namespace Event;
63
64//------------------------------------------- Constructor -----------------------------------------
65MucCalibMgr::MucCalibMgr( std::vector<double> jobInfo, std::vector<int> configInfo, std::string outputFile )
66{
67 MsgStream log(msgSvc, "MucCalibMgr");
68 log << MSG::INFO << "-------Initializing------- " << endreq;
69
70 // init Gaudi service
71 log << MSG::INFO << "Initialize Gaudi service" << endreq;
72 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
74 eventSvc = NULL;
75
76 m_jobStart = clock();
77
78 // init parameters
79 log << MSG::INFO << "Initialize parameters" << endreq;
80
81 m_vs = jobInfo[0]; m_hv = jobInfo[1]; m_th = jobInfo[2];
82
83 if( jobInfo[3] <= 0 ) m_er = TRIGGER_RATE; // 4000Hz
84 else m_er = jobInfo[3];
85
86 if( jobInfo[4] <= 0 || jobInfo[4] >5 ) m_tag = 0;
87 else m_tag = jobInfo[4];
88
89 if( configInfo[0] < 0 || configInfo[0] > 2 )
90 m_recMode = 0;
91 else
92 m_recMode = configInfo[0];
93
94 m_usePad = configInfo[1];
95
96 if( configInfo[2] < 0 || configInfo[2] > STRIP_INBOX_MAX )
97 m_effWindow = EFF_WINDOW; // 4 strip-width
98 else m_effWindow = configInfo[2];
99
100 if( configInfo[3]<0 || configInfo[3] >4 )
101 m_clusterMode = DEFAULT_BUILD_MODE;
102 else m_clusterMode = configInfo[3];
103
104 m_clusterSave = configInfo[4];
105 m_checkEvent = configInfo[5];
106 m_dimuSelect = configInfo[6];
107 m_dimuOnly = configInfo[7];
108
109 m_outputFile = outputFile;
110
111 m_fdata = NULL;
112 if( m_clusterMode ==1 && m_clusterSave == 1 )
113 m_fdata = new ofstream("MucCluster.dat", ios::out);
114
115 m_fStartRun = m_fEndRun = 0;
116 m_currentRun = m_currentEvent = 0;
117 m_eventTag = 0;
118
119 for( int i=0; i<PART_MAX; i++ )
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++ )
123 {
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;
129 }
130
131 for( int i=0; i<LAYER_MAX; i++ )
132 {
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;
138 }
139
140 for( int i=0; i<BOX_MAX; i++ )
141 {
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;
147 }
148
149 for( int i=0; i<STRIP_MAX; i++ )
150 {
151 m_stripResults[0][i] = 0;
152 m_stripResults[1][i] = 0;
153 m_stripResults[2][i] = 0;
154 m_stripResults[3][i] = 0;
155 // strip no cluster
156 }
157
158 m_fTotalDAQTime = 0;
159
160 m_ptrMucMark = new MucMark(0,0,0,0);
161 m_ptrIdTr = new MucIdTransform();
162
163 m_digiCol.resize(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);
169
170 for( int i=0; i<PART_MAX; i++ )
171 for( int j=0; j<SEGMENT_MAX; j++ ) {
172 m_segDigiCol[i][j].resize(0);
173 }
174
175 // init Gaudi Ntuple
176 InitNtuple();
177
178 // init ROOT histograms
179 InitHisto();
180
181 // init Tree
182 InitConstTree();
183
184 // init area of units
185 InitArea();
186 log << MSG::INFO << "Initialize done!" << endreq;
187}
188
189//----------------------------------------------- Destructor --------------------------------------
191{
193
197 if( m_usePad != 0 ) Clear2DEffHisto();
198
202
203 if(m_ptrMucMark) delete m_ptrMucMark;
204 if(m_ptrIdTr) delete m_ptrIdTr;
205
206 m_digiCol.clear();
207 m_calHitCol.clear();
208 m_effHitCol.clear();
209 m_nosHitCol.clear();
210 m_clusterCol.clear();
211
212 for( int i=0; i<PART_MAX; i++ )
213 for( int j=0; j<SEGMENT_MAX; j++ ) {
214 m_segDigiCol[i][j].clear();
215 }
216
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;
221}
222
223// Clear online histograms
225{
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;
230
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;
240
241 return StatusCode::SUCCESS;
242}
243
244// Clear level 0 histograms
246{
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;
258
259 return StatusCode::SUCCESS;
260}
261
262// Clear level 1 histograms
264{
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;
274
275 return StatusCode::SUCCESS;
276}
277
278// Clear level 2 histograms
280{
281 for( int i=0; i< BOX_MAX; i++ )
282 {
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];
289 }
290
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;
300
301 return StatusCode::SUCCESS;
302}
303
304// Clear 2D eff histogrames
306{
307
308 for( int i=0; i<PART_MAX; i++ )
309 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
310 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
311 {
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];
315 }
316
317 return StatusCode::SUCCESS;
318}
319
320// Clear cluster histograms
322{
323 for( int i=0; i<LAYER_MAX; i++ ) {
324 if(m_hLayerCluster[i]) delete m_hLayerCluster[i];
325 }
326
327 for( int i=0; i<BOX_MAX; i++ ) {
328 if(m_hBoxCluster[i]) delete m_hBoxCluster[i];
329 }
330
331 if(m_hLayerClusterCmp) delete m_hLayerClusterCmp;
332 if(m_hBoxClusterCmp) delete m_hBoxClusterCmp;
333
334 return StatusCode::SUCCESS;
335}
336
337// Clear resolution histograms
339{
340 for( int i=0; i<B_LAY_NUM; i++ ) {
341 if(m_hBarrelResDist[i]) delete m_hBarrelResDist[i];
342 }
343 for( int i=0; i<E_LAY_NUM; i++ ) {
344 if(m_hEndcapResDist[i]) delete m_hEndcapResDist[i];
345 }
346
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];
351
352 return StatusCode::SUCCESS;
353}
354
355// Clear Tree
357{
358
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;
364
365 return StatusCode::SUCCESS;
366}
367
368//-------------------------------------------------------------------------------------------------
369//-------------------------------- Member functions for initializing ------------------------------
370//-------------------------------------------------------------------------------------------------
371
372
373// init Ntuples
375{
376 MsgStream log(msgSvc, "MucCalibMgr");
377 log << MSG::INFO << "Initialize NTuples" << endreq;
378
379 // Book ntuple
380 Gaudi::svcLocator() -> service("NTupleSvc", ntupleSvc);
381
382 StatusCode sc;
383
384 // event log
385 NTuplePtr nt1(ntupleSvc, "FILE450/EventLog");
386 if ( nt1 ) { m_eventLogTuple = nt1; }
387 else
388 {
389 m_eventLogTuple = ntupleSvc->book ("FILE450/EventLog", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
390 if ( m_eventLogTuple )
391 {
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);
402 }
403 else {
404 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_eventLogTuple) << endmsg;
405 return StatusCode::FAILURE;
406 }
407 }
408
409 // track info
410
411 NTuplePtr nt2(ntupleSvc, "FILE450/MdcTrkInfo");
412 if ( nt2 ) { m_mdcTrkInfoTuple = nt2; }
413 else
414 {
415 m_mdcTrkInfoTuple = ntupleSvc->book ("FILE450/MdcTrkInfo", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
416 if ( m_mdcTrkInfoTuple )
417 {
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);
426 }
427 else {
428 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_mdcTrkInfoTuple) << endmsg;
429 return StatusCode::FAILURE;
430 }
431 }
432
433
434 NTuplePtr nt3(ntupleSvc, "FILE450/TrackInfo");
435 if ( nt3 ) { m_trackInfoTuple = nt3; }
436 else
437 {
438 m_trackInfoTuple = ntupleSvc->book ("FILE450/TrackInfo", CLID_RowWiseTuple, "MucCalibConst N-Tuple");
439 if ( m_trackInfoTuple )
440 {
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);
465 }
466 else {
467 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_trackInfoTuple) << endmsg;
468 return StatusCode::FAILURE;
469 }
470 }
471
472 // track collinearity
473 NTuplePtr nt4(ntupleSvc, "FILE450/TrackDiff");
474 if ( nt4 ) { m_trackDiffTuple = nt4; }
475 else
476 {
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);
484 }
485 else {
486 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_trackDiffTuple) << endmsg;
487 return StatusCode::FAILURE;
488 }
489 }
490
491 // cluster size
492 NTuplePtr nt5(ntupleSvc, "FILE450/ClusterSize");
493 if ( nt5 ) { m_clusterSizeTuple = nt5; }
494 else
495 {
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);
499 }
500 else {
501 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_clusterSizeTuple) << endmsg;
502 return StatusCode::FAILURE;
503 }
504 }
505
506 // eff window
507 NTuplePtr nt6(ntupleSvc, "FILE450/EffWindow");
508 if ( nt6 ) { m_effWindowTuple = nt6; }
509 else
510 {
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);
514 }
515 else {
516 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_effWindowTuple) << endmsg;
517 return StatusCode::FAILURE;
518 }
519 }
520 // res info
521 NTuplePtr nt7(ntupleSvc, "FILE450/ResInfo");
522 if ( nt7 ) { m_resInfoTuple = nt7; }
523 else
524 {
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);
535 }
536 else {
537 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_resInfoTuple) << endmsg;
538 return StatusCode::FAILURE;
539 }
540 }
541
542/*
543 NTuplePtr nt7(ntupleSvc, "FILE450/ResInfo");
544 if ( nt7 ) { m_resInfoTuple = nt7; }
545 else
546 {
547 m_resInfoTuple = ntupleSvc->book ("FILE450/ResInfo", CLID_ColumnWiseTuple, "MucCalibConst N-Tuple");
548 if ( m_resInfoTuple ) {
549 sc = m_resInfoTuple->addItem ("exp_num", m_nExpNum, 0, 30);
550 sc = m_resInfoTuple->addIndexedItem ("res", m_nExpNum, m_res);
551 sc = m_resInfoTuple->addIndexedItem ("res_part", m_nExpNum, m_resPart);
552 sc = m_resInfoTuple->addIndexedItem ("res_segment", m_nExpNum, m_resSegment);
553 sc = m_resInfoTuple->addIndexedItem ("res_layer", m_nExpNum, m_resLayer);
554 sc = m_resInfoTuple->addIndexedItem ("res_fired", m_nExpNum, m_resFired);
555 }
556 else {
557 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_resInfoTuple) << endmsg;
558 return StatusCode::FAILURE;
559 }
560 }
561
562 for(int i=0; i<24; i++ )
563 {
564 m_res[i] = 211;
565 m_resPart[i] = -1;
566 m_resSegment[i] = -1;
567 m_resLayer[i] = -1;
568 m_resFired[i] = false;
569 }
570*/
571
572 return StatusCode::SUCCESS;
573}
574
575// Initialize histogram maps of fired strips and reconstructed hits
577{
578 MsgStream log(msgSvc, "MucCalibMgr");
579 log << MSG::INFO << "Initialize Histograms" << endreq;
580
581 // Init online monitor histo
583
584 // Init eff map and so on
585 InitHistoLV2();
586 InitHistoLV1();
587 InitHistoLV0();
588
589 // Init 2D eff map
590 if( m_usePad != 0 ) Init2DEffHisto();
591
592 // Init cluster histo
594
595 // Init spacial resolution histo
596 InitResHisto();
597
598 return StatusCode::SUCCESS;
599}
600
601// Init histograme for online monitor
603{
604 char name[60];
605 char title[60];
606 int stripMax = 0;
607
608 // Init hit map
609 for( int i=0; i<B_LAY_NUM; i++ )
610 {
611 // According to layer
612 if(i%2 != 0 ) { stripMax = 96*7 + 112; }
613 else { stripMax = 48*8; }
614 sprintf( name, "HitMapBarrel_Lay_L%d", i );
615 sprintf( title, "Hit map in barrel layer %d", i );
616 m_hHitMapBarrel_Lay[i] = new TH1F(name,title, stripMax, 0, stripMax);
617
618 if( i<E_LAY_NUM )
619 {
620 stripMax = 64*4; // 256
621 sprintf( name, "HitMapEastEndcap_L%d", i );
622 sprintf( title, "Hit map in east-endcap layer %d", i );
623 m_hHitMapEndcap_Lay[0][i] = new TH1F(name,title, stripMax, 0, stripMax);
624
625 sprintf( name, "HitMapWestEndcap_L%d", i );
626 sprintf( title, "Hit map in west-endcap layer %d", i );
627 m_hHitMapEndcap_Lay[1][i] = new TH1F(name,title, stripMax, 0, stripMax);
628 }
629
630 }
631
632 for( int i=0; i<B_SEG_NUM; i++ )
633 {
634 // According to segment
635 sprintf( name, "HitMapBarrel_Seg_S%d", i );
636 sprintf( title, "Hit map in barrel segment %d", i );
637 if(i==2 ) { stripMax = 48*5 + 112*4; } // 688
638 else { stripMax = 48*5 + 96*4; } // 624
639 m_hHitMapBarrel_Seg[i] = new TH1F(name,title, stripMax, 0, stripMax);
640
641 if( i<E_SEG_NUM )
642 {
643 sprintf( name, "HitMapEastEndcap_S%d", i );
644 sprintf( title, "Hit map in east-endcap segment %d", i );
645 stripMax = 64*8; // 512
646 m_hHitMapEndcap_Seg[0][i] = new TH1F(name,title, stripMax, 0, stripMax);
647
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);
651 }
652 }
653
654 // Init histograms for online monitor
655 // 1D histograms
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();
661
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();
665
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();
669
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();
673
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();
677
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();
681
682 // 2D histograms
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();
688
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();
694
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();
700
701 return StatusCode::SUCCESS;
702}
703
704// Init histograme for calibration level 0
706{
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);
709
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);
719
720 return StatusCode::SUCCESS;
721}
722
723// Init histograme for calibration level 1
725{
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);
735
736 return StatusCode::SUCCESS;
737}
738
739// Init histograme for calibration level 2
741{
742 char name[60];
743 char title[60];
744 int part, segment, layer, stripMax;
745 part = segment = layer = stripMax = 0;
746
747 for( int i=0; i<BOX_MAX; i++ )
748 {
749 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
750 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
751
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);
755
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);
759
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);
763
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);
767
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);
771
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);
775 }
776
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 );
784 m_hStripNos = new TH1F("StripNos", "Strip noise", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
785 m_hStripCnt = new TH1F("StripCnt", "Strip count", STRIP_MAX+1, -0.5, STRIP_MAX+0.5 );
786
787 return StatusCode::SUCCESS;
788}
789
790// Init 2D eff histograme
792{
793 char name[60];
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.;
798
799 m_histArray = new TObjArray();
800
801 for( int i=0; i<PART_MAX; i++ )
802 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
803 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
804 {
805 boxID = m_ptrIdTr->GetBoxId(i,j,k);
806
807 if( i==BRID )
808 {
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);
814 }
815 else
816 {
817 xStart = yStart = -1250;
818 xEnd = yEnd = 1250;
819 sWidth = E_STR_DST;
820 xBin = yBin = int((xEnd - xStart)/sWidth);
821 }
822
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);
829
830 //m_histArray->Add(m_h2DExpMap[i][j][k]);
831 //m_histArray->Add(m_h2DHitMap[i][j][k]);
832 m_histArray->Add(m_h2DEffMap[i][j][k]);
833 }
834
835 return StatusCode::SUCCESS;
836}
837
838// Init histograme for cluster
840{
841 char name[60];
842 char title[60];
843 int part, segment, layer, stripMax;
844 part = segment = layer = stripMax = 0;
845
846 for( int i=0; i<LAYER_MAX; i++ )
847 {
848 sprintf( name, "LayerCluster_L%d", i );
849 sprintf( title, "Clusters in L%d",i );
850 m_hLayerCluster[i] = new TH1F(name,title, CLUSTER_RANGE, 0, CLUSTER_RANGE );
851 }
852
853 for( int i=0; i<BOX_MAX; i++ )
854 {
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 );
860 }
861
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 );
864
865 return StatusCode::SUCCESS;
866}
867
868// Init histograms for spacial resolution
870{
871 char name[60];
872 char title[60];
873
874 for( int i=0; i<B_LAY_NUM; i++ )
875 {
876 sprintf( name, "BarrelRes_L%d", i );
877 sprintf( title, "Barrel spacial resolution in L%d",i );
878 m_hBarrelResDist[i] = new TH1F(name,title, 200, -100, 100 );
879 }
880
881 for( int i=0; i<E_LAY_NUM; i++ )
882 {
883 sprintf( name, "EndcapRes_L%d", i );
884 sprintf( title, "Endcap spacial resolution in L%d",i );
885 m_hEndcapResDist[i] = new TH1F(name,title, 200, -100, 100 );
886 }
887
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 );
892
893 return StatusCode::SUCCESS;
894}
895
896// Init Tree
898{
899 MsgStream log(msgSvc, "MucCalibMgr");
900 log << MSG::INFO << "Initialize Trees" << endreq;
901
902 // job log tree
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");
916
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");
934
935 // statistic log tree
936 m_tStatLog = new TTree("StatLog","Statistic results");
937
938 // layer constants tree, level 0
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");
952
953 // box constants tree, level 1
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");
970
971 // strip constants tree, level 2
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");
987
988 return StatusCode::SUCCESS;
989}
990
991// Init area of units
993{
994 int stripMax, boxID, stripID;
995 stripMax = boxID = stripID = 0;
996 double totalArea = 0;
997
998 for( int i=0; i<PART_MAX; i++ )
999 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
1000 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
1001 {
1002 MucBoxCal* aBox = new MucBoxCal( i, j, k );
1003 boxID = m_ptrIdTr->GetBoxId(i, j, k);
1004 m_boxResults[3][boxID] = aBox->GetArea();
1005 m_layerResults[3][k] += aBox->GetArea();
1006 delete aBox;
1007
1008 stripMax = m_ptrIdTr->GetStripMax( i, j, k );
1009 for( int m=0; m<stripMax; m++ )
1010 {
1011 MucStripCal* aStrip = new MucStripCal( i, j, k, m );
1012 stripID =m_ptrIdTr->GetStripId(i, j, k, m );
1013 m_stripResults[3][stripID] = aStrip->GetArea();
1014 totalArea += m_stripResults[3][stripID];
1015 delete aStrip;
1016 }
1017 }
1018
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]);
1022
1023 m_fTotalStripArea = totalArea;
1024
1025 return StatusCode::SUCCESS;
1026}
1027
1028//-------------------------------------------------------------------------------------------------
1029//-------------------------------- Member functions for executing --------------------------------
1030//-------------------------------------------------------------------------------------------------
1031
1032// Dimu select
1034{
1035 MsgStream log(msgSvc, "MucCalibMgr");
1036 log << MSG::INFO << "Dimu select" << endreq;
1037 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
1038
1039 if( m_tag > 0) { m_eventTag = (int)m_tag; return ( StatusCode::SUCCESS ); }
1040
1041 m_eventTag = 0;
1042 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
1043 if(!eventHeader) {
1044 log << MSG::FATAL << "Could not find event header" << endreq;
1045 return( StatusCode::FAILURE );
1046 }
1047
1048 // Select by MDC Info
1049 SmartDataPtr<RecMdcTrackCol> mdcTrackCol(eventSvc,"/Event/Recon/RecMdcTrackCol");
1050 if(!mdcTrackCol) {
1051 log << MSG::FATAL << "Could not find mdc tracks" << endreq;
1052 return( StatusCode::FAILURE);
1053 }
1054 log << MSG::INFO << mdcTrackCol->size() << endreq;
1055 if( mdcTrackCol->size() != 2 ) return ( StatusCode::FAILURE );
1056
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;
1060
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;
1068
1069 // Select by TOF Info
1070 SmartDataPtr<RecTofTrackCol> tofTrackCol(eventSvc,"/Event/Recon/RecTofTrackCol");
1071 if(!tofTrackCol) {
1072 log << MSG::FATAL << "Could not find RecTofTrackCol!" << endreq;
1073 return( StatusCode::FAILURE);
1074 }
1075 log << MSG::INFO << "TOF tracks:\t" << tofTrackCol->size() << endreq;
1076
1077 double t1, t2;
1078 t1 = 0., t2 = 1000;
1079 // if(tofTrackCol->size() < 8 && tofTrackCol->size() > 20)
1080 if(tofTrackCol->size() > 7 && tofTrackCol->size() < 21)
1081 {
1082 int goodtof = 0;
1083 for(unsigned int itof = 0; itof < tofTrackCol->size(); itof++)
1084 {
1085 TofHitStatus *status = new TofHitStatus;
1086 status->setStatus((*tofTrackCol)[itof]->status());
1087
1088 if( !(status->is_cluster()) ) { delete status; continue; }
1089 if(goodtof==0) t1 = (*tofTrackCol)[itof]->tof();
1090 if(goodtof==1) t2 = (*tofTrackCol)[itof]->tof();
1091
1092 goodtof++;
1093 delete status;
1094 }
1095 }
1096 bool tofPass = false;
1097 if( fabs( t1-t2 ) < 4) tofPass = true; // ns
1098
1099 // Select by EMC Info
1100 SmartDataPtr<RecEmcShowerCol> emcShowerCol(eventSvc,"/Event/Recon/RecEmcShowerCol");
1101 if (!emcShowerCol) {
1102 log << MSG::FATAL << "Could not find RecEmcShowerCol!" << endreq;
1103 return( StatusCode::FAILURE);
1104 }
1105 log << MSG::INFO << "EMC showers:\t" << emcShowerCol->size() << endreq;
1106
1107 if( emcShowerCol->size() < 2 ) return ( StatusCode::SUCCESS );
1108
1109 double e1, e2, theta1, theta2, phi1, phi2;
1110 int part;
1111
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();
1116
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;
1120
1121 bool emcFlag1 = e1 < 1.0 && e1 > 0.1;
1122 bool emcFlag2 = e2 < 1.0 && e2 > 0.1;
1123 bool emcFlag3 = fabs(theta1 + theta2 - PI) < 0.05;
1124 bool emcFlag4 = fabs(phi1 - phi2) - PI < 0.4;
1125
1126 bool emcPass = false;
1127 if( emcFlag1 && emcFlag2 && emcFlag3 && emcFlag4 ) emcPass = true;
1128
1129 // Select by MUC Info
1130 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc,"/Event/Digi/MucDigiCol");
1131 if(!mucDigiCol) {
1132 log << MSG::FATAL << "Could not find MUC digi" << endreq;
1133 return( StatusCode::FAILURE);
1134 }
1135 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc,"/Event/Recon/RecMucTrackCol");
1136 if (!mucTrackCol) {
1137 log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
1138 return( StatusCode::FAILURE);
1139 }
1140 log << MSG::INFO << "digi:\t" << mucDigiCol->size() << "tracks:\t" << mucTrackCol->size() << endreq;
1141
1142 bool mucFlag1 = mucDigiCol->size()>6;
1143 bool mucFlag2 = mucTrackCol->size()>1;
1144 bool mucPass = false;
1145 if( mucFlag1 && mucFlag2 ) mucPass = true;
1146
1147 if( mdcPass && tofPass && emcPass && mucPass ) m_eventTag = 1;
1148 else m_eventTag = (int)m_tag;
1149
1150 return( StatusCode::SUCCESS );
1151}
1152
1153// Read event info from TDS
1155{
1156 MsgStream log(msgSvc, "MucCalibMgr");
1157 log << MSG::INFO << "Read event" << endreq;
1158
1159 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
1160 m_evtBegin = clock();
1161
1162 // Check event head
1163 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
1164 if(!eventHeader) {
1165 log << MSG::FATAL << "Could not find event header" << endreq;
1166 return( StatusCode::FAILURE );
1167 }
1168
1169 m_currentRun = eventHeader->runNumber();
1170 m_currentEvent = eventHeader->eventNumber();
1171 if( m_fStartRun == 0 ) m_fStartRun = m_currentRun;
1172 m_fEndRun = m_currentRun;
1173 m_fTotalEvent++;
1174
1175 log << MSG::INFO << "Run [ " << m_currentRun << " ]\tEvent [ " << m_currentEvent << " ]" << endreq;
1176 if( ((long)m_fTotalEvent)%2000 == 0 ) cout << m_fTotalEvent << "\tdone!" << endl;
1177
1178 // Select dimu
1179 if( m_dimuSelect ) {
1181 log << MSG::INFO << "Event tag:\t" << m_eventTag << endreq;
1182 if( m_dimuOnly && m_eventTag != 1 ) return( StatusCode::FAILURE );
1183 }
1184
1185 //---> Retrieve MUC digi
1186 log << MSG::INFO << "Retrieve digis" << endreq;
1187
1188 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc,"/Event/Digi/MucDigiCol");
1189 if(!mucDigiCol) {
1190 log << MSG::FATAL << "Could not find MUC digi" << endreq;
1191 return( StatusCode::FAILURE);
1192 }
1193
1194 int part, segment, layer, strip, pad;
1195 part = segment = layer = strip = pad = 0;
1196 double padX, padY, padZ;
1197 padX = padY = padZ = 0.;
1198 double resMax = 0.;
1199
1200 Identifier mucId;
1201 MucDigiCol::iterator digiIter = mucDigiCol->begin();
1202 int eventDigi = 0;
1203 for ( int digiId =0; digiIter != mucDigiCol->end(); digiIter++, digiId++ )
1204 {
1205 mucId = (*digiIter)->identify();
1206 part = MucID::barrel_ec(mucId);
1207 segment = MucID::segment(mucId);
1208 layer = MucID::layer(mucId);
1209 strip = MucID::channel(mucId);
1210
1211 log << MSG::DEBUG << "[" << part << "\t" << segment << "\t" << layer << "\t" << strip << "]\t" ;
1212 if( (digiId+1)%8 == 0 ) log << MSG::DEBUG << endreq;
1213
1214 eventDigi ++;
1215
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;
1218 continue;
1219 }
1220
1221 // Add digi
1222 MucMark *aMark = new MucMark( part, segment, layer, strip );
1223 m_digiCol.push_back( aMark );
1224 m_segDigiCol[part][segment].push_back( aMark );
1225 }
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;
1229
1230 /*
1231 if( (long)m_fTotalEvent%10000 == 0 ) {
1232 log << MSG::INFO << "Delete HitVsEvent." << endreq;
1233 if(m_hHitVsEvent != NULL) delete m_hHitVsEvent;
1234 m_hHitVsEvent = new TH1F("HitVsEvent","Hit VS event",10000,m_fTotalEvent,m_fTotalEvent+10000);
1235 m_hHitVsEvent->GetXaxis()->SetTitle("Event NO.");
1236 log << MSG::INFO << "Recreate HitVsEvent." << endreq;
1237 }
1238 log << MSG::INFO << "Fill HitVsEvent\t" << m_hHitVsEvent << "\t" << m_fTotalEvent << "\t" << m_currentEvent << endreq;
1239 //m_hHitVsEvent->Fill(m_currentEvent+m_currentEvent%10000, eventDigi);
1240 m_hHitVsEvent->Fill(m_fTotalEvent, eventDigi);
1241 log << MSG::INFO << "Fill HitVsEvent done." << endreq;
1242 */
1243
1244 // Search cluster in digis
1245 // Method detail in MucMark class
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 );
1251 }
1252
1253 for( unsigned int i=0; i<m_clusterCol.size(); i++ )
1254 {
1255 clusterSize = m_clusterCol[i].size();
1256 // real cluster, size >= 2
1257 if( clusterSize > CLUSTER_ALARM )
1258 {
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();
1263
1264 if( m_clusterSave ) (*m_fdata) << "Event:\t" << m_currentEvent << "\tbig cluster " << bigClusterNum << endl;
1265
1266 for( int j=0; j<clusterSize; j++ )
1267 {
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;
1272 }
1273 log << MSG::WARNING << endreq;
1274 bigClusterNum ++;
1275 }
1276 else if( clusterSize > 1 )
1277 {
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++ )
1284 {
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;
1288 }
1289 log << MSG::DEBUG << endreq;
1290 }
1291 } // End m_clusterCol.size()
1292
1293 if( m_clusterMode) log << MSG::INFO << "Total clusters in this event: " << clusterNum << endreq;
1294 else log << MSG::INFO << "Clusters not built" << endreq;
1295 //<--- End retrieve digis
1296
1297 //---> Retrieve rec tracks
1298 log << MSG::INFO << "Retrieve tracks" << endreq;
1299 // MDC tracks
1300 SmartDataPtr<RecMdcTrackCol> mdcTrackCol(eventSvc,"/Event/Recon/RecMdcTrackCol");
1301 if(!mdcTrackCol) {
1302 log << MSG::FATAL << "Could not find mdc tracks" << endreq;
1303 return( StatusCode::FAILURE);
1304 }
1305
1306 RecMdcTrackCol::iterator mdctrkIter = mdcTrackCol->begin();
1307 for (; mdctrkIter != mdcTrackCol->end(); mdctrkIter++)
1308 {
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();
1318 }
1319
1320 // MUC tracks
1321 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc,"/Event/Recon/RecMucTrackCol");
1322 if (!mucTrackCol) {
1323 log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
1324 return( StatusCode::FAILURE);
1325 }
1326
1327 RecMucTrackCol *aRecMucTrackCol = mucTrackCol;
1328 if (aRecMucTrackCol->size() < 1) {
1329 log << MSG::INFO << "No MUC tracks in this event" << endreq;
1330 return StatusCode::SUCCESS;
1331 }
1332 log << MSG::INFO << "Total tracks of this event: " << aRecMucTrackCol->size() << endreq;
1333
1334 // Get RecEsTimeCol
1335 //int esTimeflag;
1336 //if( m_recMode == 0 ) // only for ExtTrk
1337 if( 0 ) // only for ExtTrk
1338 {
1339 SmartDataPtr<RecEsTimeCol> aRecEsTimeCol(eventSvc,"/Event/Recon/RecEsTimeCol");
1340 if( ! aRecEsTimeCol ){
1341 log << MSG::ERROR << "Could not find RecEsTimeCol" << endreq;
1342 return StatusCode::FAILURE;
1343 }else{
1344 RecEsTimeCol::iterator iter_evt = aRecEsTimeCol->begin();
1345 //tes = (*iter_evt)->getTest();
1346 //esTimeflag = (*iter_evt)->getStat();
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;
1351 }
1352 }
1353 }
1354
1355 // Phi diff of two tracks( dimuon event )
1356 double phi1, phi2, phiDiff, theta1, theta2, thetaDiff;
1357 phiDiff = thetaDiff = 0.;
1358 if( aRecMucTrackCol->size()==2 && (*aRecMucTrackCol)[0]->GetTotalHits() > 4 && (*aRecMucTrackCol)[1]->GetTotalHits() > 4 )
1359 {
1360 // Pos phi diff
1361 phi1 = (*aRecMucTrackCol)[0]->getMucPos().phi(); phi2 = (*aRecMucTrackCol)[1]->getMucPos().phi();
1362 if( phi1 > 0 ) phiDiff = phi1 - phi2 - PI;
1363 else phiDiff = phi2 - phi1 - PI;
1364
1365 // Pos theta diff
1366 theta1 = (*aRecMucTrackCol)[0]->getMucPos().theta(); theta2 = (*aRecMucTrackCol)[1]->getMucPos().theta();
1367 thetaDiff = theta1 + theta2 - PI;
1368 m_hTrackPosPhiDiff->Fill( phiDiff );
1369 m_hTrackPosThetaDiff->Fill( thetaDiff );
1370 m_hDimuTracksPosDiff->Fill( thetaDiff, phiDiff );
1371 m_ntPosPhiDiff = phiDiff;
1372 m_ntPosThetaDiff = thetaDiff;
1373
1374 log << MSG::INFO << "PosPhiDiff:\t" << phiDiff << "\tPosThetaDiff:\t" << thetaDiff << endreq;
1375
1376 // Mom phi diff
1377 phi1 = (*aRecMucTrackCol)[0]->getMucMomentum().phi(); phi2 = (*aRecMucTrackCol)[1]->getMucMomentum().phi();
1378 if( phi1 > 0 ) phiDiff = phi1 - phi2 - PI;
1379 else phiDiff = phi2 - phi1 - PI;
1380
1381 // Mom theta diff
1382 theta1 = (*aRecMucTrackCol)[0]->getMucMomentum().theta(); theta2 = (*aRecMucTrackCol)[1]->getMucMomentum().theta();
1383 thetaDiff = theta1 + theta2 - PI;
1384
1385 m_hTrackMomPhiDiff->Fill( phiDiff );
1386 m_hTrackMomThetaDiff->Fill( thetaDiff );
1387 m_hDimuTracksMomDiff->Fill( thetaDiff, phiDiff );
1388 m_ntMomPhiDiff = phiDiff;
1389 m_ntMomThetaDiff = thetaDiff;
1390
1391 log << MSG::INFO << "MomPhiDiff:\t" << phiDiff << "\tMomThetaDiff:\t" << thetaDiff << endreq;
1392 m_ntDimuTag = m_eventTag;
1393 m_trackDiffTuple->write();
1394 }
1395
1396 // Retrieve tracks for calibration
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;
1407 }
1408
1409 mark_col trkSeg[TRACK_SEG_MAX];
1410
1411 // vector<MucRecHit*> mucRawHitCol;
1412 // vector<MucRecHit*> mucExpHitCol; //==========
1413 vector<int> mucRawHitCol;
1414 vector<int> mucExpHitCol;
1415
1416 for (int trackId = 0; trackIter != mucTrackCol->end(); trackIter++, trackId++)
1417 {
1418
1419//cout<<__LINE__<<"trackID "<<trackId<<endl;
1420 //trackHitNum = (*trackIter)->GetTotalHits();
1421 trackHitNum = (*trackIter)->numHits();
1422 log << MSG::DEBUG << "Track: " << trackId << " Hits: " << trackHitNum << endreq;
1423//cout<< "Track: " << trackId << " Hits: " << trackHitNum << endl;
1424 if( trackHitNum == 0 ) {
1425 log << MSG::INFO << "Track " << trackId << " no hits" << endreq;
1426 continue;
1427 }
1428
1429 m_ntTrackHits = trackHitNum;
1430
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);
1438
1439 // First fit position in MUC
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();
1450
1451
1452 m_hPhiCosTheta->Fill(m_cosTheta, m_phi);
1453 log << MSG::INFO << "Fill track info" << endreq;
1454
1455// MucRecHit* pMucRawHit;
1456// MucRecHit* pMucExpHit;
1457 if( m_calHitCol.size() != 0 ) m_calHitCol.clear(); // Fresh each track
1458
1459 // Digis belong to this rec track
1460 log << MSG::DEBUG << "Reconstruction hits(digis in a track): " << endreq;
1461 //mucRawHitCol = (*trackIter)->GetHits(); // Get hit collection of a track
1462 mucRawHitCol =(*trackIter)->getVecHits();
1463 rawHitNum += mucRawHitCol.size();
1464//cout<<__LINE__<<"rawHitNum "<<rawHitNum<<endl;
1465 segNum = 0;
1466/*
1467 if( trkRecMode == 3 ) { // By SlfTrk
1468 for(int iPart=0; iPart<PART_MAX; iPart++)
1469 for(int iLayer=0; iLayer<LAYER_MAX; iLayer++) seedList[iPart][iLayer] = false;
1470 }
1471*/
1472 for(unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
1473 {
1474 /* pMucRawHit = mucRawHitCol[ hitId ];
1475 part = pMucRawHit->Part();
1476 segment = pMucRawHit->Seg();
1477 layer = pMucRawHit->Gap();
1478 strip = pMucRawHit->Strip();
1479 */
1480 mucId = mucRawHitCol[hitId];
1481 int part = MucID::barrel_ec(mucId);
1482 int segment = MucID::segment(mucId);
1483 int layer = MucID::layer(mucId);
1484 int strip = MucID::channel(mucId);
1485
1486//cout<<__LINE__<<" mu p s l s "<<mucId<<" "<<part<<" "<<segment<<" "<<layer<<" "<<strip<<endl;
1487 log << MSG::DEBUG << "[" << part << "\t" << segment << "\t" << layer << "\t" << strip << "]\t";
1488 //if( (hitId+1)%8 == 0 ) log << MSG::DEBUG << endreq;
1489
1490 // Add hit
1491 MucMark *aMark = new MucMark( part, segment, layer, strip );
1492 m_calHitCol.push_back( aMark );
1493
1494 // Set seed flag
1495// if( trkRecMode == 3 ) seedList[part][layer] = pMucRawHit->HitIsSeed(); // By SlfTrk
1496
1497 // Find track segment
1498 if(hitId == 0) { trkSeg[segNum].push_back( aMark ); segNum ++; }
1499 else
1500 {
1501 log << MSG::DEBUG << "segNum: " << segNum << endreq;
1502 bool notInSeg = true;
1503 for( int segi=0; segi<segNum; segi++ )
1504 {
1505 if( aMark->IsInSegWith( *(trkSeg[segi][0]) ) )
1506 {
1507 trkSeg[segi].push_back( aMark );
1508 notInSeg = false;
1509 break;
1510 }
1511 }
1512 // new track seg
1513 if( notInSeg == true )
1514 {
1515 trkSeg[segNum].push_back( aMark );
1516 segNum ++;
1517 if( segNum > TRACK_SEG_MAX ) {
1518 log << MSG::ERROR << "Track segment overflow: " << segNum << endreq;
1519 break;
1520 }
1521 }
1522 } // End else
1523 } // End raw hits
1524 log << MSG::DEBUG << endreq;
1525
1526 // Find maximal layers passed in track segments
1527 layerPassNum[0] = layerPassNum[1] = layerPassNum[2] = 0;
1528 for( int segi=0; segi<segNum; segi++ )
1529 {
1530 int tmpLayNum = 0;
1531 passMax[segi][0] = passMax[segi][1] = trkSeg[segi][0]->Layer();
1532 for( unsigned int hiti=1; hiti<trkSeg[segi].size(); hiti++ )
1533 {
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;
1539 }
1540
1541 for( int layi=0; layi<LAYER_MAX; layi++ ) {
1542 if( firedLay[segi][layi] ) tmpLayNum ++;
1543 }
1544
1545 if( segi == 0 ) layerPassNum[0] += passMax[segi][1] + 1;
1546 else layerPassNum[0] += (passMax[segi][1] - passMax[segi][0] + 1);
1547
1548 layerPassNum[1] += (passMax[segi][1] - passMax[segi][0] + 1);
1549 layerPassNum[2] += tmpLayNum;
1550
1551 trkSeg[segi].clear();
1552 }
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;
1562 //if( layerPassNum[0]>B_LAY_NUM || layerPassNum[1]>B_LAY_NUM || layerPassNum[2]>B_LAY_NUM )
1563 // log << MSG::ERROR << "Over max layer:\t" << m_currentRun << "\t" << m_currentEvent << endreq;
1564
1565 // Expected hits in this rec track
1566 log << MSG::DEBUG << "Fitting hits(expected hits in a track): " << endreq;
1567 //mucExpHitCol = (*trackIter)->GetExpectedHits();
1568 mucExpHitCol = (*trackIter)->getExpHits();
1569
1570 expectedHitNum += mucExpHitCol.size();
1571//cout<<__LINE__<<"expectedHitNum "<<expectedHitNum<<endl;
1572 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
1573 {
1574// pMucRawHit = mucExpHitCol[ hitId ];
1575// part = pMucRawHit->Part(); segment = pMucRawHit->Seg();
1576// layer = pMucRawHit->Gap(); strip = pMucRawHit->Strip();
1577
1578 mucId =mucExpHitCol[hitId];
1579 int part = MucID::barrel_ec(mucId);
1580 int segment = MucID::segment(mucId);
1581 int layer = MucID::layer(mucId);
1582 int strip = MucID::channel(mucId);
1583//cout<<__LINE__<<part<<" "<<segment<<" "<<layer<<" "<<strip<<endl;
1584/*
1585 if( m_usePad != 0 )
1586 {
1587 pad = pMucRawHit->GetPadID(); padZ = pMucRawHit->GetIntersectZ();
1588 padX = pMucRawHit->GetIntersectY(); padY = pMucRawHit->GetIntersectX(); // Note: local coordinate
1589
1590 if( part != BRID )
1591 {
1592 if(segment == 1) { padX = -padX; }
1593 else if(segment == 2) { padX = -padX, padY = -padY; }
1594 else if(segment == 3) { padY = -padY; }
1595 }
1596 }
1597*/
1598 // Avoid bias in seed layers
1599 // if( seedList[part][layer] == true ) continue;
1600
1601 MucMark* currentMark = new MucMark( part, segment, layer, strip );
1602 m_expHitCol.push_back( currentMark );
1603 //log << MSG::DEBUG << "[" << part << "\t" << segment << "\t" << layer << "\t" << strip << "]\t";
1604 //if( (hitId+1)%8 == 0 ) log << MSG::DEBUG << endreq;
1605
1606 // Judge efficiency hit
1607 int isInPos = -1;
1608 bool isInEffWindow = false;
1609 isInPos = currentMark->IsInCol( m_segDigiCol[part][segment] );
1610
1611 // Avoid bias in outer layers caused by low momentum tracks
1612 if( part == BRID && (layer-lastLayerBR>1) ) continue;
1613 if( part != BRID && (layer-lastLayerEC>1) ) continue;
1614
1615 // Avoid bias in both sides of the innermost layer of Barrel
1616 if( part==BRID && layer==0 && (strip<2 || strip>45) )
1617 {
1618 if( isInPos != -1) // expHit is fired
1619 {
1620 m_record[part][segment][layer][strip][2] ++; // Efficiency hit number
1621 m_record[part][segment][layer][strip][1] ++; // Rec track number
1622 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
1623
1624 if( m_usePad != 0 ) {
1625 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1626 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
1627 }
1628
1629 }
1630 else {
1631 m_record[part][segment][layer][strip][1] ++;
1632 if( m_usePad != 0 ) m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1633 }
1634 continue; // Judge next hit
1635 }
1636
1637 // Eff calibration
1638 if( isInPos != -1 ) // expHit is fired
1639 {
1640 m_record[part][segment][layer][strip][2] ++; // Efficiency hit number
1641 m_record[part][segment][layer][strip][1] ++; // Rec track number
1642 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
1643
1644 if( m_usePad != 0 ) {
1645 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1646 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
1647 }
1648
1649 continue; // Judge next hit
1650 }
1651 else for(int tempStrip=0, hiti=-m_effWindow; hiti<=m_effWindow; hiti++ )
1652 {
1653 if( hiti == 0 ) continue;
1654 tempStrip = strip + hiti;
1655 if( tempStrip < 0 || tempStrip > m_ptrIdTr->GetStripMax(part,segment,layer) ) continue;
1656
1657 isInPos = m_ptrMucMark->IsInCol( part, segment, layer, tempStrip, m_segDigiCol[part][segment] );
1658 if( isInPos != -1 )
1659 {
1660 m_record[part][segment][layer][tempStrip][2] ++; // Efficiency hit number
1661 m_record[part][segment][layer][tempStrip][1] ++; // Rec track number
1662 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
1663
1664 if( m_usePad != 0 ) {
1665 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1666 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
1667 }
1668
1669 m_ntEffWindow = hiti;
1670 m_effWindowTuple->write();
1671 isInEffWindow = true;
1672 }
1673
1674 } // End else
1675
1676 if( isInEffWindow ) { continue; } // Judge next hit
1677 else { // A hit should be fired but not fired and not in the EffWindow
1678 m_record[part][segment][layer][strip][1] ++; // Rec track number
1679 if( m_usePad != 0 ) m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1680 }
1681
1682 } // End expected hits
1683
1684 // Fill residual, and for the other way of eff calculation
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();
1690
1691 for(unsigned int nres = 0; nres < m_lineResCol.size(); nres++ )
1692 if( fabs(m_lineResCol[nres])>resMax ) resMax = fabs(m_lineResCol[nres]);
1693
1694 log << MSG::INFO << "Good track for res" << endreq;
1695 if( trackHitNum > 4 && m_lineResCol[0] != -99) // track is good for res
1696 {
1697 // Fill res histograms
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;
1702
1703 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
1704 {
1705// pMucExpHit = mucExpHitCol[ hitId ];
1706// part = pMucExpHit->Part(); segment = pMucExpHit->Seg(); layer = pMucExpHit->Gap();
1707
1708 mucId = mucExpHitCol[hitId];
1709
1710 int part = MucID::barrel_ec(mucId);
1711 int segment = MucID::segment(mucId);
1712 int layer = MucID::layer(mucId);
1713//cout<<__LINE__<<" m p s l"<<mucId<<" "<<part<<" "<<segment<<" "<<layer<<endl;
1714 firedFlag[part][layer][0] = true;
1715 }
1716
1717 log << MSG::INFO << "Fit res" << endreq;
1718 for(unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
1719 {
1720 // pMucRawHit = mucRawHitCol[ hitId ];
1721 // part = pMucRawHit->Part(); segment = pMucRawHit->Seg(); layer = pMucRawHit->Gap();
1722
1723 mucId = mucRawHitCol[hitId];
1724 int part = MucID::barrel_ec(mucId);
1725 int segment = MucID::segment(mucId);
1726 int layer = MucID::layer(mucId);
1727//cout<<__LINE__<<" part "<< mucId<<" "<<part<<" "<<segment<<" "<<layer<<endl;
1728 if( part == BRID ) m_hBarrelResDist[layer]->Fill( m_lineResCol[hitId] );
1729 else m_hEndcapResDist[layer]->Fill( m_lineResCol[hitId] );
1730
1731 // if exp is true and fired is true, and not filled yet
1732 if( firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false )
1733 {
1734 m_resPart = part;
1735 m_resSegment = segment;
1736 m_resLayer = layer;
1737 m_lineRes = m_lineResCol[hitId];
1738 //m_quadRes = m_quadResCol[hitId];
1739 //m_extrRes = m_extrResCol[hitId];
1740 m_quadRes = 999;
1741 m_extrRes = 999;
1742 m_resFired = 1;
1743 m_resMode = mode;
1744 m_resInfoTuple->write();
1745//cout<<__LINE__<<" m_resInfoTuple->write "<<endl;
1746 }
1747
1748 firedFlag[part][layer][1] = true;
1749 }
1750
1751 log << MSG::INFO << "Exp res" << endreq;
1752 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
1753 {
1754// pMucExpHit = mucExpHitCol[ hitId ];
1755// part = pMucExpHit->Part(); segment = pMucExpHit->Seg(); layer = pMucExpHit->Gap();
1756
1757 mucId = mucExpHitCol[hitId];
1758 int part = MucID::barrel_ec(mucId);
1759 int segment = MucID::segment(mucId);
1760 int layer = MucID::layer(mucId);
1761//cout<<__LINE__<<mucId<<" "<<part<<" "<<segment<<" "<<layer<<endl;
1762 if(firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false)
1763 {
1764 m_resPart = part;
1765 m_resSegment = segment;
1766 m_resLayer = layer;
1767 m_lineRes = 1000;
1768 m_quadRes = 1000;
1769 m_extrRes = 1000;
1770 m_resFired = 0;
1771 m_resMode = mode;
1772 m_resInfoTuple->write();
1773 }
1774 }
1775
1776 } // End fill residual, if track is good for res
1777
1778 mucRawHitCol.clear();
1779 mucExpHitCol.clear();
1780
1781 } // End read all tracks
1782
1783 //if( resMax > 300 ) cout <<"Res too big!\t"<< m_fTotalEvent <<"\t"<< m_currentRun <<"\t"<< m_currentEvent <<"\t"<< resMax << endl;
1784
1785 m_ntTrackNum = mucTrackCol->size();
1786
1787 m_fTotalEffHit += rawHitNum;
1788 log << MSG::INFO << "Total hits in this event, raw: " << rawHitNum << "\texpected: " << expectedHitNum << endreq;
1789 //<--- End retrieve rec tracks
1790
1791 //---> Searching inc/noise hits
1792 log << MSG::INFO << "Searching inc/noise hits" << endreq;
1793 bool isNosHit;
1794 bool hasEffHit;
1795 for( unsigned int i=0; i < m_digiCol.size(); i++ )
1796 {
1797 isNosHit = true;
1798
1799 if( m_digiCol[i]->IsInCol( m_effHitCol ) !=-1) continue; // digi in effHitCol
1800 else
1801 {
1802 for( unsigned int j=0; j < m_clusterCol.size(); j++ )
1803 {
1804 hasEffHit = false;
1805 for( unsigned int k=0; k<m_clusterCol[j].size(); k++)
1806 {
1807 if( m_clusterCol[j][k]->IsInCol(m_effHitCol) != -1) // Clusters have efficiency hit
1808 {
1809 hasEffHit = true;
1810 break; // Out a cluster
1811 }
1812 }
1813
1814 if ( hasEffHit && (m_digiCol[i]->IsInCol( m_clusterCol[j] ) != -1) ) {
1815 isNosHit = false;
1816 break; // Out cluster collection
1817 }
1818 } // End cluster col
1819
1820 if( isNosHit ) {
1821 m_nosHitCol.push_back( m_digiCol[i] );
1822 m_fTotalNosHit ++;
1823 }
1824 }// End else
1825 } // End digi collection
1826
1827 return StatusCode::SUCCESS;
1828}
1829
1830// Check event
1832{
1833 MsgStream log(msgSvc, "MucCalibMgr");
1834 log << MSG::INFO << "Check event" << endreq;
1835
1836 // Find over digis in digis
1837 // Note that only one digi relates to one strip
1838 log << MSG::INFO << "Searching over digi" << endreq;
1839 int overDigi = 0;
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 ++;
1843 }
1844
1845 if( overDigi !=0 )
1846 log << MSG::ERROR << overDigi << " over digi found in digis!" << endreq;
1847
1848 // Find over hits in reconstruction hits
1849 // Note that two tracks may pass thought one strip
1850 log << MSG::INFO << "Searching over hit" << endreq;
1851 int overHit = 0;
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 ++;
1855 }
1856
1857 if( overHit != 0 )
1858 log << MSG::WARNING << overHit << " over hits found in rec hits!" << endreq;
1859
1860 // Find hits of reconstruction not in MUCdigis
1861 log << MSG::INFO << "Searching new hit" << endreq;
1862 int newHit = 0;
1863 int num = 0;
1864 for( unsigned int i=0; i < m_expHitCol.size(); i++ )
1865 {
1866 num = m_expHitCol[i]->NumInCol( m_digiCol );
1867 if( num == 0 )
1868 {
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;
1874
1875 newHit ++;
1876 }
1877 }
1878
1879 if( newHit != 0 )
1880 log << MSG::WARNING << newHit << " new hit(s) found in rec hits!" << endreq;
1881
1882 return StatusCode::SUCCESS;
1883}
1884
1885// Fill event
1887{
1888 MsgStream log(msgSvc, "MucCalibMgr");
1889 log << MSG::INFO << "Fill event" << endreq;
1890
1891 int part, segment, layer, strip, size;
1892 part = segment = layer = strip = size = 0;
1893
1894 // Fill digis
1895 log << MSG::INFO << "Fill digis" << endreq;
1896 for( unsigned int i=0; i<m_digiCol.size(); i++ )
1897 {
1898 part = (*m_digiCol[i]).Part();
1899 segment = (*m_digiCol[i]).Segment();
1900 layer = (*m_digiCol[i]).Layer();
1901 strip = (*m_digiCol[i]).Strip();
1902
1903 FillDigi( part, segment, layer, strip );
1904 m_record[part][segment][layer][strip][0] ++;
1905 m_fTotalDigi ++;
1906 }
1907
1908 // Fill rec hits
1909 log << MSG::INFO << "Fill rec hits" << endreq;
1910 for( unsigned int i=0; i<m_expHitCol.size(); i++ )
1911 {
1912 part = (*m_expHitCol[i]).Part();
1913 segment = (*m_expHitCol[i]).Segment();
1914 layer = (*m_expHitCol[i]).Layer();
1915 strip = (*m_expHitCol[i]).Strip();
1916
1917 FillExpHit( part, segment, layer, strip );
1918 m_record[part][segment][layer][strip][4] ++;
1919 m_fTotalExpHit ++;
1920 }
1921
1922 // Fill eff hits
1923 log << MSG::INFO << "Fill eff hits" << endreq;
1924 for( unsigned int i=0; i<m_effHitCol.size(); i++ )
1925 {
1926 part = (*m_effHitCol[i]).Part();
1927 segment = (*m_effHitCol[i]).Segment();
1928 layer = (*m_effHitCol[i]).Layer();
1929 strip = (*m_effHitCol[i]).Strip();
1930
1931 FillEffHit( part, segment, layer, strip );
1932 m_fTotalEffHit ++;
1933 }
1934
1935 // Fill clusters
1936 log << MSG::INFO << "Fill clusters" << endreq;
1937 for( unsigned int i=0; i<m_clusterCol.size(); i++ )
1938 {
1939 size = m_clusterCol[i].size();
1940 // may include the special cluster, size = 1
1941 if( size > CLUSTER_CUT )
1942 {
1943 part = (*m_clusterCol[i][0]).Part();
1944 segment = (*m_clusterCol[i][0]).Segment();
1945 layer = (*m_clusterCol[i][0]).Layer();
1946
1947 FillCluster( part, segment, layer, size );
1948 m_ntClusterSize = size;
1949 m_clusterSizeTuple->write();
1950 }
1951 }
1952
1953 // Fill inc/noise hits
1954 log << MSG::INFO << "Fill inc/noise hits" << endreq;
1955 for( unsigned int i=0; i<m_nosHitCol.size(); i++ )
1956 {
1957 part = (*m_nosHitCol[i]).Part();
1958 segment = (*m_nosHitCol[i]).Segment();
1959 layer = (*m_nosHitCol[i]).Layer();
1960 strip = (*m_nosHitCol[i]).Strip();
1961
1962 FillNosHit( part, segment, layer, strip );
1963 m_record[part][segment][layer][strip][3] ++;
1964 }
1965
1966 // Event log
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();
1974
1975 // Reset mark collection
1976 for( unsigned int i=0; i<m_digiCol.size(); i++ ) {
1977 if( m_digiCol[i] != NULL ) delete m_digiCol[i];
1978 }
1979
1980 for( unsigned int i=0; i<m_expHitCol.size(); i++ ) {
1981 if( m_expHitCol[i] != NULL ) delete m_expHitCol[i];
1982 }
1983/*
1984 for( unsigned int i=0; i<m_effHitCol.size(); i++ ) {
1985 if( m_effHitCol[i] != NULL ) delete m_effHitCol[i];
1986 }
1987 log << MSG::INFO << m_effHitCol.size() << endreq;
1988*/
1989 for( unsigned int i=0; i<m_calHitCol.size(); i++ ) {
1990 if( m_calHitCol[i] != NULL ) delete m_calHitCol[i];
1991 }
1992
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();
1999
2000
2001 for( int i=0; i<PART_MAX; i++ ) {
2002 for( int j=0; j<SEGMENT_MAX; j++ ) {
2003 if( m_segDigiCol[i][j].size() != 0 ) m_segDigiCol[i][j].clear();
2004 }
2005 }
2006
2007 m_evtEnd = clock();
2008
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();
2012
2013 return StatusCode::SUCCESS;
2014}
2015
2016StatusCode MucCalibMgr::FillDigi( int part, int segment, int layer, int strip )
2017{
2018 // Hit map for online check
2019 int tmpId;
2020
2021 if( (int)m_tag || m_dimuOnly==0 || (m_dimuOnly==1&&m_eventTag==1) )
2022 {
2023 if( part == BRID )
2024 {
2025 // According to layer
2026 if( layer%2 == 0 ) {
2027 if( segment > 4 ) tmpId = segment*48 + (47 - strip);
2028 else tmpId = segment*48 + strip;
2029 }
2030 else if( segment < 3 ) tmpId = segment*96 + strip;
2031 else tmpId = (segment-1)*96 + 112 + strip;
2032
2033 m_hHitMapBarrel_Lay[layer]->Fill(tmpId);
2034
2035 // According to segment
2036 if( segment != B_TOP ) {
2037 if( segment > 4 )
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;
2040 }
2041 else tmpId = 48*((layer+1)/2) + 112*(layer/2) + strip;
2042
2043 m_hHitMapBarrel_Seg[segment]->Fill(tmpId);
2044 }
2045 else if( part == EEID )
2046 {
2047 // According to layer
2048 tmpId = segment*64 + strip;
2049 m_hHitMapEndcap_Lay[0][layer]->Fill(tmpId);
2050 // According to segment
2051 tmpId = layer*64 + strip;
2052 m_hHitMapEndcap_Seg[0][segment]->Fill(tmpId);
2053 }
2054 else if( part == EWID )
2055 {
2056 // According to layer
2057 tmpId = segment*64 + strip;
2058 m_hHitMapEndcap_Lay[1][layer]->Fill(tmpId);
2059 // According to segment
2060 tmpId = layer*64 + strip;
2061 m_hHitMapEndcap_Seg[1][segment]->Fill(tmpId);
2062 }
2063 }
2064
2065 // Total units
2066 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2067 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
2068
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) );
2076
2077 return StatusCode::SUCCESS;
2078}
2079
2080StatusCode MucCalibMgr::FillExpHit( int part, int segment, int layer, int strip )
2081{
2082 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2083 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
2084
2085 m_hStripExpHitMap[boxId]->Fill( strip );
2086 m_hStripExpHit->Fill( strId );
2087 m_hBoxExpHit->Fill( boxId );
2088 m_hLayerExpHit->Fill( layer );
2089
2090 return StatusCode::SUCCESS;
2091}
2092
2093StatusCode MucCalibMgr::FillEffHit( int part, int segment, int layer, int strip )
2094{
2095 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2096 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
2097
2098 m_hStripEffHitMap[boxId]->Fill( strip );
2099 m_hStripEffHit->Fill( strId );
2100 m_hBoxEffHit->Fill( boxId );
2101 m_hLayerEffHit->Fill( layer );
2102
2103 return StatusCode::SUCCESS;
2104}
2105
2106StatusCode MucCalibMgr::FillNosHit( int part, int segment, int layer, int strip )
2107{
2108 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2109 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
2110
2111 m_hStripNosHitMap[boxId]->Fill( strip );
2112 m_hStripNosHit->Fill( strId );
2113 m_hBoxNosHit->Fill( boxId );
2114 m_hLayerNosHit->Fill( layer );
2115
2116 return StatusCode::SUCCESS;
2117}
2118
2119StatusCode MucCalibMgr::FillCluster( int part, int segment, int layer, int size )
2120{
2121 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2122
2123 m_hLayerCluster[layer]->Fill( size );
2124 m_hBoxCluster[boxId]->Fill( size );
2125
2126 return StatusCode::SUCCESS;
2127}
2128
2129
2130//-------------------------------------------------------------------------------------------------
2131//-------------------------------- Member functions for finalizing --------------------------------
2132//-------------------------------------------------------------------------------------------------
2133
2135{
2136 MsgStream log(msgSvc, "MucCalibMgr");
2137 log<< MSG::INFO << "Start efficiency and noise calibration" << endreq;
2138
2139 // Set time
2140 m_fTotalDAQTime = m_fTotalEvent/m_er;
2141 log<< MSG::INFO << "Total run time:\t" << m_fTotalDAQTime << endreq;
2142
2146
2147 if( m_usePad != 0 ) PadEff();
2148
2149 return StatusCode::SUCCESS;
2150}
2151
2153{
2154 MsgStream log(msgSvc, "MucCalibMgr");
2155 log<< MSG::INFO << "Analyse layer efficiency and noise" << endreq;
2156
2157 int part, segment, layer, stripMax;
2158 part = segment = layer = stripMax = 0;
2159
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.;
2163
2164 for( int i=0; i<LAYER_MAX; i++ )
2165 {
2166 digi = effHit = trackNum = nosHit = recHit = 0;
2167
2168 for( int j=0; j<PART_MAX; j++ )
2169 {
2170 for( int k=0; k<SEGMENT_MAX; k++)
2171 {
2172 stripMax = m_ptrIdTr->GetStripMax( j, k, i );
2173 for( int n=0; n<stripMax; n++ )
2174 {
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];
2180 }
2181 }
2182 }
2183
2184 // Efficiency
2185 if( trackNum == 0 ) {
2186 eff = effErr = 0.0;
2187 //eff = DEFAULT_EFF_VALUE;
2188 //effErr = DEFAULT_EFF_ERR;
2189 }
2190 else
2191 {
2192 eff = ( (double)effHit )/trackNum;
2193 effErr = sqrt( eff*(1-eff)/trackNum );
2194 m_fCalibLayerNum ++;
2195 }
2196
2197 // Noise
2198 if( m_layerResults[3][i] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
2199 noise = DEFAULT_NOS_VALUE;
2200 else {
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]);
2203 }
2204
2205 if( digi != 0 ) {
2206 nosRatio = ( (double)nosHit )/digi;
2207 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
2208 }
2209 else {
2210 nosRatio = DEFAULT_INC_VALUE;
2211 nosRatioErr = 0;
2212 }
2213
2214 // Counting rate
2215 if( m_recMode == 2 )
2216 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_layerResults[3][i]);
2217 else
2218 cnt = (double)digi/(m_fTotalDAQTime * m_layerResults[3][i]);
2219
2220 // Cluster
2221 cluster = m_hLayerCluster[i]->GetMean();
2222
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;
2228
2229 // Fill histogram
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 );
2236
2237 // Add cluster histogram
2238 m_hLayerClusterCmp->Fill( i, cluster );
2239
2240 // Fill tree
2241 m_fLayerId = i;
2242 m_fLayerEff = eff;
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;
2251 m_fLayerCnt = cnt;
2252 m_tLayConst->Fill();
2253
2254 // Add cluster ntuple
2255 m_fLayerCluster = cluster;
2256
2257 } // End layer
2258
2259 return StatusCode::SUCCESS;
2260}
2261
2263{
2264 MsgStream log(msgSvc, "MucCalibMgr");
2265 log<< MSG::INFO << "Analyse box efficiency and noise" << endreq;
2266
2267 int part, segment, layer, stripMax;
2268 part = segment = layer = stripMax = 0;
2269
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.;
2273
2274 for( int i=0; i<BOX_MAX; i++ )
2275 {
2276 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
2277 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
2278
2279 digi = effHit = trackNum = nosHit = recHit = 0;
2280 for( int j=0; j<stripMax; j++ )
2281 {
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];
2287 }
2288
2289 // Efficiency
2290 if( trackNum == 0 ) {
2291 eff = effErr = 0.0;
2292 //eff = DEFAULT_EFF_VALUE;
2293 //effErr = DEFAULT_EFF_ERR;
2294 }
2295 else
2296 {
2297 eff = ( (double)effHit )/trackNum;
2298 effErr = sqrt( eff*(1-eff)/trackNum );
2299 m_fCalibBoxNum ++;
2300 }
2301
2302 // Noise
2303 if( m_boxResults[3][i] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
2304 noise = DEFAULT_NOS_VALUE;
2305 else {
2306 if( m_recMode == 2 )
2307 noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW * m_boxResults[3][i]);
2308 else
2309 noise = (double)nosHit/(m_fTotalDAQTime * m_boxResults[3][i]);
2310 }
2311
2312 if( digi != 0 ) {
2313 nosRatio = ( (double)nosHit )/digi;
2314 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
2315 }
2316 else {
2317 nosRatio = DEFAULT_INC_VALUE;
2318 nosRatioErr = 0;
2319 }
2320
2321 // Counting rate
2322 if( m_recMode == 2 )
2323 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_boxResults[3][i]);
2324 else
2325 cnt = (double)digi/(m_fTotalDAQTime * m_boxResults[3][i]);
2326
2327 // Cluster
2328 cluster = m_hBoxCluster[i]->GetMean();
2329
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;
2335
2336 // Fill histograms
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 );
2343
2344 // add cluster histogram
2345 m_hBoxClusterCmp->Fill( i, cluster );
2346
2347 // Fill tree
2348 m_fBoxId = i;
2349 m_fBoxPart = part;
2350 m_fBoxSegment = segment;
2351 m_fBoxLayer = layer;
2352 m_fBoxEff = eff;
2353 m_fBoxEffErr = effErr;
2354 m_fBoxTrkNum = trackNum;
2355 m_fBoxExpHit = recHit;
2356 m_fBoxEffHit = effHit;
2357 m_fBoxNosRatio = nosRatio;
2358 m_fBoxDigi = digi;
2359 m_fBoxNosHit = nosHit;
2360 m_fBoxNos = noise;
2361 m_fBoxCnt = cnt;
2362 m_tBoxConst->Fill();
2363
2364 // add box cluster ntuple
2365 m_fBoxCluster = cluster;
2366
2367 }// End BOX_MAX
2368
2369 return StatusCode::SUCCESS;
2370}
2371
2373{
2374 MsgStream log(msgSvc, "MucCalibMgr");
2375 log<< MSG::INFO << "Analyse strip efficiency and noise" << endreq;
2376
2377 int part, segment, layer, stripMax;
2378 part = segment = layer = stripMax = 0;
2379
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.;
2383
2384 for( int i=0; i<BOX_MAX; i++ )
2385 {
2386 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
2387 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
2388
2389 for( int j=0; j<stripMax; j++ )
2390 {
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];
2396
2397 int stripId = m_ptrIdTr->GetStripId( part, segment, layer, j );
2398
2399 // Efficiency
2400 if( trackNum == 0 ) {
2401 eff = effErr = 0.0;
2402 //eff = DEFAULT_EFF_VALUE;
2403 //effErr = DEFAULT_EFF_ERR;
2404 }
2405 else
2406 {
2407 eff = ( (double)effHit )/trackNum;
2408 effErr = sqrt( eff*(1-eff)/trackNum );
2409 m_fCalibStripNum ++;
2410 }
2411
2412 // Noise
2413 if( m_stripResults[3][stripId] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
2414 noise = DEFAULT_NOS_VALUE;
2415 else {
2416 if( m_recMode == 2 )
2417 noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW * m_stripResults[3][stripId]);
2418 else
2419 noise = (double)nosHit/(m_fTotalDAQTime * m_stripResults[3][stripId]);
2420 }
2421
2422 if( digi != 0 ) {
2423 nosRatio = ( (double)nosHit )/digi;
2424 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
2425 }
2426 else {
2427 nosRatio = DEFAULT_INC_VALUE;
2428 nosRatioErr = 0.;
2429 }
2430
2431 // Counting rate
2432 if( m_recMode == 2 )
2433 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_stripResults[3][stripId]);
2434 else
2435 cnt = (double)digi/(m_fTotalDAQTime * m_stripResults[3][stripId]);
2436
2437 // Strip itself no clusters
2438 m_stripResults[0][ stripId ] = eff;
2439 m_stripResults[1][ stripId ] = effErr;
2440 m_stripResults[2][ stripId ] = noise;
2441 m_stripResults[5][ stripId ] = trackNum;
2442
2443 // Fill histotram
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 );
2454
2455 // Fill Tree
2456 m_fStripId = stripId;
2457 m_fStripPart = part;
2458 m_fStripSegment = segment;
2459 m_fStripLayer = layer;
2460 m_fStripEff = eff;
2461 m_fStripEffErr = effErr;
2462 m_fStripNosRatio = nosRatio;
2463 m_fStripDigi = digi;
2464 m_fStripNos = noise;
2465 m_fStripCnt = cnt;
2466 m_fStripEffHit = effHit;
2467 m_fStripExpHit = recHit;
2468 m_fStripNosHit = nosHit;
2469 m_fStripTrkNum = trackNum;
2470 m_tStrConst->Fill();
2471
2472 } // End stripMax
2473 } // End BOX_MAX
2474
2475 return StatusCode::SUCCESS;
2476}
2477
2479{
2480 MsgStream log(msgSvc, "MucCalibMgr");
2481
2482 int xBinStart, xBinEnd, yBinStart, yBinEnd;
2483 double expHit, firedHit, eff;
2484 eff = expHit = firedHit = 0.;
2485
2486 for( int i=0; i<PART_MAX; i++ )
2487 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
2488 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
2489 {
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;
2494
2495 for( int xi = xBinStart; xi<xBinEnd; xi++ )
2496 for( int yi = yBinStart; yi<yBinEnd; yi++ )
2497 {
2498 expHit = m_h2DExpMap[i][j][k]->GetBinContent(xi, yi);
2499 firedHit = m_h2DHitMap[i][j][k]->GetBinContent(xi, yi);
2500
2501 if( expHit !=0 ) eff = firedHit / expHit;
2502 else eff = 0;
2503
2504 if( eff>1.0 )
2505 cout<<"Eff error:\t["<<i<<"\t"<<j<<"\t"<<k<<",\t"<<xi<<"\t"<<yi<<"]:\t"
2506 <<eff<<"\t,"<<firedHit<<" "<<expHit<<endl;
2507
2508 m_h2DEffMap[i][j][k]->SetBinContent(xi, yi, eff);
2509 }
2510 }
2511 return StatusCode::SUCCESS;
2512}
2513
2515{
2516 MsgStream log(msgSvc, "MucCalibMgr");
2517
2518 return StatusCode::SUCCESS;
2519}
2520
2522{
2523 MsgStream log(msgSvc, "MucCalibMgr");
2524 log << MSG::INFO << "Analyse spacial resolution" << endreq;
2525
2526 double resMean, resMeanErr, resSigma, resSigmaErr;
2527 resMean = resMeanErr = resSigma = resSigmaErr = 0.;
2528
2529 for( int i=0; i<B_LAY_NUM; i++ )
2530 {
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);
2536
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 );
2541 }
2542
2543 for( int i=0; i<E_LAY_NUM; i++ )
2544 {
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);
2550
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 );
2555 }
2556
2557 return StatusCode::SUCCESS;
2558}
2559
2561{
2562 MsgStream log(msgSvc, "MucCalibMgr");
2563 log << MSG::INFO << "Save calibration constants" << endreq;
2564
2565 // Save calibrated eff in graphes
2566 // LV0
2567 double layerXY[2][LAYER_MAX];
2568 double layerEXY[2][LAYER_MAX];
2569 for( int i=0; i<LAYER_MAX; i++ )
2570 {
2571 layerXY[0][i] = i;
2572 layerEXY[0][i] = 0;
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];
2576 }
2577 else {
2578 layerXY[1][i] = 0;
2579 layerEXY[1][i] = 0;
2580 }
2581 }
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");
2589 m_cv[0]->Write();
2590
2591 // LV1
2592 double boxXY[2][BOX_MAX];
2593 double boxEXY[2][BOX_MAX];
2594 for( int i=0; i<BOX_MAX; i++ )
2595 {
2596 boxXY[0][i] = i;
2597 boxEXY[0][i] = 0;
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];
2601 }
2602 else {
2603 boxXY[1][i] = 0;
2604 boxEXY[1][i] = 0;
2605 }
2606 }
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");
2614 m_cv[1]->Write();
2615
2616 // LV2
2617 double stripXY[2][STRIP_MAX];
2618 double stripEXY[2][STRIP_MAX];
2619 for( int i=0; i<STRIP_MAX; i++ )
2620 {
2621 stripXY[0][i] = i;
2622 stripEXY[0][i] = 0;
2623 if( m_stripResults[5][i] >= TRACK_THRESHOLD ) {
2624 stripXY[1][i] = m_stripResults[0][i];
2625 stripEXY[1][i] = m_stripResults[1][i];
2626 }
2627 else {
2628 stripXY[1][i] = 0;
2629 stripEXY[1][i] = 0;
2630 }
2631 }
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");
2639 m_cv[2]->Write();
2640
2641 // Save histograms
2642 for(int i=0; i<B_LAY_NUM; i++ )
2643 {
2644 m_hHitMapBarrel_Lay[i]->Write();
2645
2646 if( i<E_LAY_NUM ) {
2647 m_hHitMapEndcap_Lay[0][i]->Write();
2648 m_hHitMapEndcap_Lay[1][i]->Write();
2649 }
2650 }
2651
2652 for(int i=0; i<B_SEG_NUM; i++)
2653 {
2654 m_hHitMapBarrel_Seg[i]->Write();
2655
2656 if( i< E_SEG_NUM ) {
2657 m_hHitMapEndcap_Seg[0][i]->Write();
2658 m_hHitMapEndcap_Seg[1][i]->Write();
2659 }
2660 }
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");
2666
2667 m_hTrackDistance->Write();
2668 m_hTrackPosPhiDiff->Write();
2669 m_hTrackPosThetaDiff->Write();
2670 m_hTrackMomPhiDiff->Write();
2671 m_hTrackMomThetaDiff->Write();
2672
2673 for( int i=0; i<B_LAY_NUM; i++ ) m_hBarrelResDist[i]->Write();
2674 for( int i=0; i<E_LAY_NUM; i++ ) m_hEndcapResDist[i]->Write();
2675
2676 m_hBarrelResComp[0]->Write();
2677 m_hBarrelResComp[1]->Write();
2678 m_hEndcapResComp[0]->Write();
2679 m_hEndcapResComp[1]->Write();
2680
2681 if( m_usePad != 0 ) m_histArray->Write();
2682
2683 for( int i=0; i<BOX_MAX; i++ )
2684 {
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();
2691 }
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;
2701
2702 m_hBoxFire->Write();
2703 m_hBoxExpHit->Write();
2704 m_hBoxEffHit->Write();
2705 m_hBoxNosHit->Write();
2706 m_hBoxEff->Write();
2707 m_hBoxArea->Write();
2708 m_hBoxNos->Write();
2709 m_hBoxNosRatio->Write();
2710 log << MSG::INFO << "Save LV1 histograms done!" << endreq;
2711
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();
2722
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();
2727
2728 log << MSG::INFO << "Save histograms done!" << endreq;
2729
2730 // Save trees
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;
2734
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");
2741
2742 int stripMax;
2743 for( int i=0; i<PART_MAX; i++ )
2744 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
2745 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
2746 {
2747 stripMax = m_ptrIdTr->GetStripMax( i, j, k );
2748 for( int n=0; n<stripMax; n++ )
2749 {
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];
2755 m_tStatLog->Fill();
2756 }
2757 }
2758
2759 m_jobFinish = clock();
2760 m_fTotalJobTime = (double)(m_jobFinish - m_jobStart)/CLOCKS_PER_SEC;
2761
2762 m_tJobLog->Fill();
2763 m_tJobLog->Write();
2764 m_tStatLog->Write();
2765 m_tLayConst->Write();
2766 m_tBoxConst->Write();
2767 m_tStrConst->Write();
2768
2769 // Close cluster output file
2770 if( m_fdata != NULL ) m_fdata->close();
2771
2772 return StatusCode::SUCCESS;
2773}
2774
2776{
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;
2785
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;
2792
2793 log << MSG::INFO << "Calibration run successfully" << endreq << endreq;
2794
2795 return StatusCode::SUCCESS;
2796}
2797
2798// END
curve Write()
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)
const Int_t n
Double_t phi2
Double_t phi1
Double_t e1
Double_t e2
titledef title[20]
const int STRIP_MAX
Definition: MucMappingAlg.h:24
vector< MucMark * > mark_col
Definition: MucMark.h:17
const double PI
Definition: PipiJpsi.cxx:55
ObjectVector< RecMucTrack > RecMucTrackCol
Definition: RecMucTrack.h:394
#define Segment
Definition: TTrackBase.h:31
@ theta2
Definition: TrkKalDeriv.h:24
@ theta1
Definition: TrkKalDeriv.h:24
#define NULL
StatusCode ClearHistoLV2()
StatusCode InitResHisto()
StatusCode FillEffHit(int part, int segment, int layer, int strip)
StatusCode DimuSelect()
StatusCode InitOnlineHisto()
StatusCode Init2DEffHisto()
StatusCode SaveConst()
StatusCode FillExpHit(int part, int segment, int layer, int strip)
StatusCode ClearResHisto()
StatusCode ClearHistoLV0()
StatusCode ClearClusterHisto()
StatusCode AnalyseCluster()
StatusCode ClearHistoLV1()
StatusCode PadEff()
StatusCode InitHistoLV1()
StatusCode InitHistoLV2()
StatusCode InitClusterHisto()
StatusCode FillCluster(int part, int segment, int layer, int size)
StatusCode ClearOnlineHisto()
StatusCode InitHisto()
StatusCode FillEvent()
StatusCode EffAndNoiseLV2()
INTupleSvc * ntupleSvc
Definition: MucCalibMgr.h:100
MucCalibMgr(std::vector< double > jobInfo, std::vector< int > configInfo, std::string outputFile)
Definition: MucCalibMgr.cxx:65
StatusCode InitHistoLV0()
IDataProviderSvc * eventSvc
Definition: MucCalibMgr.h:101
StatusCode EndRun()
StatusCode ReadEvent()
StatusCode EffAndNoiseLV1()
StatusCode AnalyseRes()
IMessageSvc * msgSvc
Definition: MucCalibMgr.h:99
StatusCode CheckEvent()
StatusCode FillNosHit(int part, int segment, int layer, int strip)
StatusCode EffAndNoiseLV0()
StatusCode InitConstTree()
StatusCode ClearConstTree()
StatusCode InitArea()
StatusCode AnalyseEffAndNoise()
StatusCode InitNtuple()
StatusCode Clear2DEffHisto()
StatusCode FillDigi(int part, int segment, int layer, int strip)
double GetArea()
static int barrel_ec(const Identifier &id)
Values of different levels.
Definition: MucID.cxx:41
static int layer(const Identifier &id)
Definition: MucID.cxx:61
static int channel(const Identifier &id)
Definition: MucID.cxx:71
static int segment(const Identifier &id)
Definition: MucID.cxx:51
int GetStripId(int part, int segment, int layer, int subid)
bool SetBoxPos(int boxid, int *part, int *segment, int *layer)
int GetStripMax(int part, int segment, int layer)
int GetBoxId(int part, int segment, int layer)
bool IsInSegWith(MucMark &other)
Definition: MucMark.cxx:144
int IsInCol(int part, int segment, int layer, int strip, mark_col &aCol)
Definition: MucMark.cxx:100
bool is_cluster() const
Definition: TofHitStatus.h:25
void setStatus(unsigned int status)
Definition: Event.h:21