BOSS 7.0.7
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
34#include "Identifier/Identifier.h"
35#include "Identifier/MucID.h"
36#include "EventModel/EventHeader.h"
37#include "EventModel/EventModel.h"
38#include "EvTimeEvent/RecEsTime.h"
39
40#include "McTruth/McKine.h"
41#include "McTruth/McParticle.h"
42#include "McTruth/MucMcHit.h"
43
44#include "MdcRecEvent/RecMdcTrack.h"
45#include "EmcRecEventModel/RecEmcEventModel.h"
46#include "TofRecEvent/RecTofTrack.h"
47#include "DstEvent/TofHitStatus.h"
48
49//#include "ReconEvent/ReconEvent.h"
50#include "MucRawEvent/MucDigi.h"
51#include "MucRecEvent/MucRecHit.h"
52#include "MucRecEvent/RecMucTrack.h"
53#include "MucRecEvent/MucRecHitContainer.h"
54
55#include "MucCalibAlg/MucStructConst.h"
56#include "MucCalibAlg/MucCalibMgr.h"
57#include "MucCalibAlg/MucMark.h"
58#include "MucCalibAlg/MucBoxCal.h"
59#include "MucCalibAlg/MucStripCal.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);
73 ntupleSvc = NULL;
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 --------------------------------------
190MucCalibMgr::~MucCalibMgr()
191{
192 ClearOnlineHisto();
193
194 ClearHistoLV2();
195 ClearHistoLV1();
196 ClearHistoLV0();
197 if( m_usePad != 0 ) Clear2DEffHisto();
198
199 ClearClusterHisto();
200 ClearResHisto();
201 ClearConstTree();
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
224StatusCode MucCalibMgr::ClearOnlineHisto()
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
245StatusCode MucCalibMgr::ClearHistoLV0()
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
263StatusCode MucCalibMgr::ClearHistoLV1()
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
279StatusCode MucCalibMgr::ClearHistoLV2()
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
305StatusCode MucCalibMgr::Clear2DEffHisto()
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
321StatusCode MucCalibMgr::ClearClusterHisto()
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
338StatusCode MucCalibMgr::ClearResHisto()
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
356StatusCode MucCalibMgr::ClearConstTree()
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
374StatusCode MucCalibMgr::InitNtuple()
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
576StatusCode MucCalibMgr::InitHisto()
577{
578 MsgStream log(msgSvc, "MucCalibMgr");
579 log << MSG::INFO << "Initialize Histograms" << endreq;
580
581 // Init online monitor histo
582 InitOnlineHisto();
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
593 InitClusterHisto();
594
595 // Init spacial resolution histo
596 InitResHisto();
597
598 return StatusCode::SUCCESS;
599}
600
601// Init histograme for online monitor
602StatusCode MucCalibMgr::InitOnlineHisto()
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
705StatusCode MucCalibMgr::InitHistoLV0()
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
724StatusCode MucCalibMgr::InitHistoLV1()
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
740StatusCode MucCalibMgr::InitHistoLV2()
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
791StatusCode MucCalibMgr::Init2DEffHisto()
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
839StatusCode MucCalibMgr::InitClusterHisto()
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
869StatusCode MucCalibMgr::InitResHisto()
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
897StatusCode MucCalibMgr::InitConstTree()
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
992StatusCode MucCalibMgr::InitArea()
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
1033StatusCode MucCalibMgr::DimuSelect()
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
1154StatusCode MucCalibMgr::ReadEvent()
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 ) {
1180 MucCalibMgr::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 vector<MucRecHit*> mucRawHitCol;
1411 vector<MucRecHit*> mucExpHitCol;
1412
1413 for (int trackId = 0; trackIter != mucTrackCol->end(); trackIter++, trackId++)
1414 {
1415 trackHitNum = (*trackIter)->GetTotalHits();
1416 log << MSG::DEBUG << "Track: " << trackId << " Hits: " << trackHitNum << endreq;
1417
1418 if( trackHitNum == 0 ) {
1419 log << MSG::INFO << "Track " << trackId << " no hits" << endreq;
1420 continue;
1421 }
1422
1423 m_ntTrackHits = trackHitNum;
1424
1425 m_trkRecMode = trkRecMode = (*trackIter)->GetRecMode();
1426 m_chi2 = (*trackIter)->chi2();
1427 m_px = (*trackIter)->getMucMomentum().x();
1428 m_py = (*trackIter)->getMucMomentum().y();
1429 m_pz = (*trackIter)->getMucMomentum().z();
1430 m_pt = sqrt(m_px*m_px + m_py*m_py);
1431 m_pp = sqrt(m_px*m_px + m_py*m_py + m_pz*m_pz);
1432
1433 // First fit position in MUC
1434 m_r = (*trackIter)->getMucPos().mag();
1435 m_cosTheta = (*trackIter)->getMucPos().cosTheta();
1436 m_theta = (*trackIter)->getMucPos().theta();
1437 m_phi = (*trackIter)->getMucPos().phi();
1438 m_depth = (*trackIter)->depth();
1439 m_brLastLayer = lastLayerBR = (*trackIter)->brLastLayer();
1440 m_ecLastLayer = lastLayerEC = (*trackIter)->ecLastLayer();
1441 m_totalHits = (*trackIter)->numHits();
1442 m_totalLayers = (*trackIter)->numLayers();
1443 m_maxHitsInLayer = (*trackIter)->maxHitsInLayer();
1444
1445
1446 m_hPhiCosTheta->Fill(m_cosTheta, m_phi);
1447 log << MSG::INFO << "Fill track info" << endreq;
1448
1449 MucRecHit* pMucRawHit;
1450 MucRecHit* pMucExpHit;
1451 if( m_calHitCol.size() != 0 ) m_calHitCol.clear(); // Fresh each track
1452
1453 // Digis belong to this rec track
1454 log << MSG::DEBUG << "Reconstruction hits(digis in a track): " << endreq;
1455 mucRawHitCol = (*trackIter)->GetHits(); // Get hit collection of a track
1456 rawHitNum += mucRawHitCol.size();
1457
1458 segNum = 0;
1459 if( trkRecMode == 3 ) { // By SlfTrk
1460 for(int iPart=0; iPart<PART_MAX; iPart++)
1461 for(int iLayer=0; iLayer<LAYER_MAX; iLayer++) seedList[iPart][iLayer] = false;
1462 }
1463
1464 for(unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
1465 {
1466 pMucRawHit = mucRawHitCol[ hitId ];
1467 part = pMucRawHit->Part();
1468 segment = pMucRawHit->Seg();
1469 layer = pMucRawHit->Gap();
1470 strip = pMucRawHit->Strip();
1471
1472 log << MSG::DEBUG << "[" << part << "\t" << segment << "\t" << layer << "\t" << strip << "]\t";
1473 //if( (hitId+1)%8 == 0 ) log << MSG::DEBUG << endreq;
1474
1475 // Add hit
1476 MucMark *aMark = new MucMark( part, segment, layer, strip );
1477 m_calHitCol.push_back( aMark );
1478
1479 // Set seed flag
1480 if( trkRecMode == 3 ) seedList[part][layer] = pMucRawHit->HitIsSeed(); // By SlfTrk
1481
1482 // Find track segment
1483 if(hitId == 0) { trkSeg[segNum].push_back( aMark ); segNum ++; }
1484 else
1485 {
1486 log << MSG::DEBUG << "segNum: " << segNum << endreq;
1487 bool notInSeg = true;
1488 for( int segi=0; segi<segNum; segi++ )
1489 {
1490 if( aMark->IsInSegWith( *(trkSeg[segi][0]) ) )
1491 {
1492 trkSeg[segi].push_back( aMark );
1493 notInSeg = false;
1494 break;
1495 }
1496 }
1497 // new track seg
1498 if( notInSeg == true )
1499 {
1500 trkSeg[segNum].push_back( aMark );
1501 segNum ++;
1502 if( segNum > TRACK_SEG_MAX ) {
1503 log << MSG::ERROR << "Track segment overflow: " << segNum << endreq;
1504 break;
1505 }
1506 }
1507 } // End else
1508 } // End raw hits
1509 log << MSG::DEBUG << endreq;
1510
1511 // Find maximal layers passed in track segments
1512 layerPassNum[0] = layerPassNum[1] = layerPassNum[2] = 0;
1513 for( int segi=0; segi<segNum; segi++ )
1514 {
1515 int tmpLayNum = 0;
1516 passMax[segi][0] = passMax[segi][1] = trkSeg[segi][0]->Layer();
1517 for( unsigned int hiti=1; hiti<trkSeg[segi].size(); hiti++ )
1518 {
1519 if( trkSeg[segi][hiti]->Layer() < passMax[segi][0] )
1520 passMax[segi][0] = trkSeg[segi][hiti]->Layer();
1521 if( trkSeg[segi][hiti]->Layer() > passMax[segi][1] )
1522 passMax[segi][1] = trkSeg[segi][hiti]->Layer();
1523 firedLay[segi][trkSeg[segi][hiti]->Layer()] = 1;
1524 }
1525
1526 for( int layi=0; layi<LAYER_MAX; layi++ ) {
1527 if( firedLay[segi][layi] ) tmpLayNum ++;
1528 }
1529
1530 if( segi == 0 ) layerPassNum[0] += passMax[segi][1] + 1;
1531 else layerPassNum[0] += (passMax[segi][1] - passMax[segi][0] + 1);
1532
1533 layerPassNum[1] += (passMax[segi][1] - passMax[segi][0] + 1);
1534 layerPassNum[2] += tmpLayNum;
1535
1536 trkSeg[segi].clear();
1537 }
1538 m_ntTrackEvent = m_currentEvent;
1539 m_ntTrackTag = m_eventTag;
1540 m_ntTrackSegFly = segNum;
1541 m_ntTrackLayFlyA = layerPassNum[0];
1542 m_ntTrackLayFlyB = layerPassNum[1];
1543 m_ntTrackLayFlyC = layerPassNum[2];
1544 m_trackInfoTuple->write();
1545 log << MSG::INFO << "Track\t" << trackId << "\tsegment(s):\t" << segNum
1546 << "\tlayer passed:\t" << layerPassNum[0] <<"\t" << layerPassNum[1] << "\t" << layerPassNum[2] << endreq;
1547 //if( layerPassNum[0]>B_LAY_NUM || layerPassNum[1]>B_LAY_NUM || layerPassNum[2]>B_LAY_NUM )
1548 // log << MSG::ERROR << "Over max layer:\t" << m_currentRun << "\t" << m_currentEvent << endreq;
1549
1550 // Expected hits in this rec track
1551 log << MSG::DEBUG << "Fitting hits(expected hits in a track): " << endreq;
1552 mucExpHitCol = (*trackIter)->GetExpectedHits();
1553 expectedHitNum += mucExpHitCol.size();
1554 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
1555 {
1556 pMucRawHit = mucExpHitCol[ hitId ];
1557 part = pMucRawHit->Part(); segment = pMucRawHit->Seg();
1558 layer = pMucRawHit->Gap(); strip = pMucRawHit->Strip();
1559
1560 if( m_usePad != 0 )
1561 {
1562 pad = pMucRawHit->GetPadID(); padZ = pMucRawHit->GetIntersectZ();
1563 padX = pMucRawHit->GetIntersectY(); padY = pMucRawHit->GetIntersectX(); // Note: local coordinate
1564
1565 if( part != BRID )
1566 {
1567 if(segment == 1) { padX = -padX; }
1568 else if(segment == 2) { padX = -padX, padY = -padY; }
1569 else if(segment == 3) { padY = -padY; }
1570 }
1571 }
1572
1573 // Avoid bias in seed layers
1574 // if( seedList[part][layer] == true ) continue;
1575
1576 MucMark* currentMark = new MucMark( part, segment, layer, strip );
1577 m_expHitCol.push_back( currentMark );
1578 //log << MSG::DEBUG << "[" << part << "\t" << segment << "\t" << layer << "\t" << strip << "]\t";
1579 //if( (hitId+1)%8 == 0 ) log << MSG::DEBUG << endreq;
1580
1581 // Judge efficiency hit
1582 int isInPos = -1;
1583 bool isInEffWindow = false;
1584 isInPos = currentMark->IsInCol( m_segDigiCol[part][segment] );
1585
1586 // Avoid bias in outer layers caused by low momentum tracks
1587 if( part == BRID && (layer-lastLayerBR>1) ) continue;
1588 if( part != BRID && (layer-lastLayerEC>1) ) continue;
1589
1590 // Avoid bias in both sides of the innermost layer of Barrel
1591 if( part==BRID && layer==0 && (strip<2 || strip>45) )
1592 {
1593 if( isInPos != -1) // expHit is fired
1594 {
1595 m_record[part][segment][layer][strip][2] ++; // Efficiency hit number
1596 m_record[part][segment][layer][strip][1] ++; // Rec track number
1597 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
1598
1599 if( m_usePad != 0 ) {
1600 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1601 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
1602 }
1603
1604 }
1605 else {
1606 m_record[part][segment][layer][strip][1] ++;
1607 if( m_usePad != 0 ) m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1608 }
1609 continue; // Judge next hit
1610 }
1611
1612 // Eff calibration
1613 if( isInPos != -1 ) // expHit is fired
1614 {
1615 m_record[part][segment][layer][strip][2] ++; // Efficiency hit number
1616 m_record[part][segment][layer][strip][1] ++; // Rec track number
1617 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
1618
1619 if( m_usePad != 0 ) {
1620 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1621 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
1622 }
1623
1624 continue; // Judge next hit
1625 }
1626 else for(int tempStrip=0, hiti=-m_effWindow; hiti<=m_effWindow; hiti++ )
1627 {
1628 if( hiti == 0 ) continue;
1629 tempStrip = strip + hiti;
1630 if( tempStrip < 0 || tempStrip > m_ptrIdTr->GetStripMax(part,segment,layer) ) continue;
1631
1632 isInPos = m_ptrMucMark->IsInCol( part, segment, layer, tempStrip, m_segDigiCol[part][segment] );
1633 if( isInPos != -1 )
1634 {
1635 m_record[part][segment][layer][tempStrip][2] ++; // Efficiency hit number
1636 m_record[part][segment][layer][tempStrip][1] ++; // Rec track number
1637 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
1638
1639 if( m_usePad != 0 ) {
1640 m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1641 m_h2DHitMap[part][segment][layer]->Fill(padX, padY);
1642 }
1643
1644 m_ntEffWindow = hiti;
1645 m_effWindowTuple->write();
1646 isInEffWindow = true;
1647 }
1648
1649 } // End else
1650
1651 if( isInEffWindow ) { continue; } // Judge next hit
1652 else { // A hit should be fired but not fired and not in the EffWindow
1653 m_record[part][segment][layer][strip][1] ++; // Rec track number
1654 if( m_usePad != 0 ) m_h2DExpMap[part][segment][layer]->Fill(padX, padY);
1655 }
1656
1657 } // End expected hits
1658
1659 // Fill residual, and for the other way of eff calculation
1660 log << MSG::INFO << "Fill residual" << endreq;
1661 vector<float> m_lineResCol = (*trackIter)->getDistHits();
1662 vector<float> m_quadResCol = (*trackIter)->getQuadDistHits();
1663 vector<float> m_extrResCol = (*trackIter)->getExtDistHits();
1664 int mode = (*trackIter)->GetRecMode();
1665
1666 for(unsigned int nres = 0; nres < m_lineResCol.size(); nres++ )
1667 if( fabs(m_lineResCol[nres])>resMax ) resMax = fabs(m_lineResCol[nres]);
1668
1669 log << MSG::INFO << "Good track for res" << endreq;
1670 if( trackHitNum > 4 && m_lineResCol[0] != -99) // track is good for res
1671 {
1672 // Fill res histograms
1673 bool firedFlag[PART_MAX][LAYER_MAX][2];
1674 for(int iprt=0; iprt<PART_MAX; iprt++)
1675 for(int jlay=0; jlay<LAYER_MAX; jlay++)
1676 firedFlag[iprt][jlay][0] = firedFlag[iprt][jlay][1] = false;
1677
1678 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
1679 {
1680 pMucExpHit = mucExpHitCol[ hitId ];
1681 part = pMucExpHit->Part(); segment = pMucExpHit->Seg(); layer = pMucExpHit->Gap();
1682 firedFlag[part][layer][0] = true;
1683 }
1684
1685 log << MSG::INFO << "Fit res" << endreq;
1686 for(unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
1687 {
1688 pMucRawHit = mucRawHitCol[ hitId ];
1689 part = pMucRawHit->Part(); segment = pMucRawHit->Seg(); layer = pMucRawHit->Gap();
1690
1691 if( part == BRID ) m_hBarrelResDist[layer]->Fill( m_lineResCol[hitId] );
1692 else m_hEndcapResDist[layer]->Fill( m_lineResCol[hitId] );
1693
1694 // if exp is true and fired is true, and not filled yet
1695 if( firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false )
1696 {
1697 m_resPart = part;
1698 m_resSegment = segment;
1699 m_resLayer = layer;
1700 m_lineRes = m_lineResCol[hitId];
1701 m_quadRes = m_quadResCol[hitId];
1702 m_extrRes = m_extrResCol[hitId];
1703 m_resFired = 1;
1704 m_resMode = mode;
1705 m_resInfoTuple->write();
1706 }
1707
1708 firedFlag[part][layer][1] = true;
1709 }
1710
1711 log << MSG::INFO << "Exp res" << endreq;
1712 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
1713 {
1714 pMucExpHit = mucExpHitCol[ hitId ];
1715 part = pMucExpHit->Part(); segment = pMucExpHit->Seg(); layer = pMucExpHit->Gap();
1716
1717 if(firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false)
1718 {
1719 m_resPart = part;
1720 m_resSegment = segment;
1721 m_resLayer = layer;
1722 m_lineRes = 1000;
1723 m_quadRes = 1000;
1724 m_extrRes = 1000;
1725 m_resFired = 0;
1726 m_resMode = mode;
1727 m_resInfoTuple->write();
1728 }
1729 }
1730
1731 } // End fill residual, if track is good for res
1732
1733 mucRawHitCol.clear();
1734 mucExpHitCol.clear();
1735
1736 } // End read all tracks
1737
1738 if( resMax > 300 ) cout <<"Res too big!\t"<< m_fTotalEvent <<"\t"<< m_currentRun <<"\t"<< m_currentEvent <<"\t"<< resMax << endl;
1739
1740 m_ntTrackNum = mucTrackCol->size();
1741
1742 m_fTotalEffHit += rawHitNum;
1743 log << MSG::INFO << "Total hits in this event, raw: " << rawHitNum << "\texpected: " << expectedHitNum << endreq;
1744 //<--- End retrieve rec tracks
1745
1746 //---> Searching inc/noise hits
1747 log << MSG::INFO << "Searching inc/noise hits" << endreq;
1748 bool isNosHit;
1749 bool hasEffHit;
1750 for( unsigned int i=0; i < m_digiCol.size(); i++ )
1751 {
1752 isNosHit = true;
1753
1754 if( m_digiCol[i]->IsInCol( m_effHitCol ) !=-1) continue; // digi in effHitCol
1755 else
1756 {
1757 for( unsigned int j=0; j < m_clusterCol.size(); j++ )
1758 {
1759 hasEffHit = false;
1760 for( unsigned int k=0; k<m_clusterCol[j].size(); k++)
1761 {
1762 if( m_clusterCol[j][k]->IsInCol(m_effHitCol) != -1) // Clusters have efficiency hit
1763 {
1764 hasEffHit = true;
1765 break; // Out a cluster
1766 }
1767 }
1768
1769 if ( hasEffHit && (m_digiCol[i]->IsInCol( m_clusterCol[j] ) != -1) ) {
1770 isNosHit = false;
1771 break; // Out cluster collection
1772 }
1773 } // End cluster col
1774
1775 if( isNosHit ) {
1776 m_nosHitCol.push_back( m_digiCol[i] );
1777 m_fTotalNosHit ++;
1778 }
1779 }// End else
1780 } // End digi collection
1781
1782 return StatusCode::SUCCESS;
1783}
1784
1785// Check event
1786StatusCode MucCalibMgr::CheckEvent()
1787{
1788 MsgStream log(msgSvc, "MucCalibMgr");
1789 log << MSG::INFO << "Check event" << endreq;
1790
1791 // Find over digis in digis
1792 // Note that only one digi relates to one strip
1793 log << MSG::INFO << "Searching over digi" << endreq;
1794 int overDigi = 0;
1795 for( unsigned int i=0; i < m_digiCol.size(); i++ )
1796 for( unsigned int j=i+1; j < m_digiCol.size(); j++ ) {
1797 if( (*m_digiCol[i]) == (*m_digiCol[j]) ) overDigi ++;
1798 }
1799
1800 if( overDigi !=0 )
1801 log << MSG::ERROR << overDigi << " over digi found in digis!" << endreq;
1802
1803 // Find over hits in reconstruction hits
1804 // Note that two tracks may pass thought one strip
1805 log << MSG::INFO << "Searching over hit" << endreq;
1806 int overHit = 0;
1807 for( unsigned int i=0; i < m_expHitCol.size(); i++ )
1808 for( unsigned int j=i+1; j < m_expHitCol.size(); j++ ) {
1809 if( (*m_expHitCol[i]) == (*m_expHitCol[j]) ) overHit ++;
1810 }
1811
1812 if( overHit != 0 )
1813 log << MSG::WARNING << overHit << " over hits found in rec hits!" << endreq;
1814
1815 // Find hits of reconstruction not in MUCdigis
1816 log << MSG::INFO << "Searching new hit" << endreq;
1817 int newHit = 0;
1818 int num = 0;
1819 for( unsigned int i=0; i < m_expHitCol.size(); i++ )
1820 {
1821 num = m_expHitCol[i]->NumInCol( m_digiCol );
1822 if( num == 0 )
1823 {
1824 log << MSG::ERROR << "Exp hit not in digis!"
1825 << "prt: " << (*m_expHitCol[i]).Part()
1826 << "\tseg: " << (*m_expHitCol[i]).Segment()
1827 << "\tlay: " << (*m_expHitCol[i]).Layer()
1828 << "\tstr: " << (*m_expHitCol[i]).Strip() << endreq;
1829
1830 newHit ++;
1831 }
1832 }
1833
1834 if( newHit != 0 )
1835 log << MSG::WARNING << newHit << " new hit(s) found in rec hits!" << endreq;
1836
1837 return StatusCode::SUCCESS;
1838}
1839
1840// Fill event
1841StatusCode MucCalibMgr::FillEvent()
1842{
1843 MsgStream log(msgSvc, "MucCalibMgr");
1844 log << MSG::INFO << "Fill event" << endreq;
1845
1846 int part, segment, layer, strip, size;
1847 part = segment = layer = strip = size = 0;
1848
1849 // Fill digis
1850 log << MSG::INFO << "Fill digis" << endreq;
1851 for( unsigned int i=0; i<m_digiCol.size(); i++ )
1852 {
1853 part = (*m_digiCol[i]).Part();
1854 segment = (*m_digiCol[i]).Segment();
1855 layer = (*m_digiCol[i]).Layer();
1856 strip = (*m_digiCol[i]).Strip();
1857
1858 FillDigi( part, segment, layer, strip );
1859 m_record[part][segment][layer][strip][0] ++;
1860 m_fTotalDigi ++;
1861 }
1862
1863 // Fill rec hits
1864 log << MSG::INFO << "Fill rec hits" << endreq;
1865 for( unsigned int i=0; i<m_expHitCol.size(); i++ )
1866 {
1867 part = (*m_expHitCol[i]).Part();
1868 segment = (*m_expHitCol[i]).Segment();
1869 layer = (*m_expHitCol[i]).Layer();
1870 strip = (*m_expHitCol[i]).Strip();
1871
1872 FillExpHit( part, segment, layer, strip );
1873 m_record[part][segment][layer][strip][4] ++;
1874 m_fTotalExpHit ++;
1875 }
1876
1877 // Fill eff hits
1878 log << MSG::INFO << "Fill eff hits" << endreq;
1879 for( unsigned int i=0; i<m_effHitCol.size(); i++ )
1880 {
1881 part = (*m_effHitCol[i]).Part();
1882 segment = (*m_effHitCol[i]).Segment();
1883 layer = (*m_effHitCol[i]).Layer();
1884 strip = (*m_effHitCol[i]).Strip();
1885
1886 FillEffHit( part, segment, layer, strip );
1887 m_fTotalEffHit ++;
1888 }
1889
1890 // Fill clusters
1891 log << MSG::INFO << "Fill clusters" << endreq;
1892 for( unsigned int i=0; i<m_clusterCol.size(); i++ )
1893 {
1894 size = m_clusterCol[i].size();
1895 // may include the special cluster, size = 1
1896 if( size > CLUSTER_CUT )
1897 {
1898 part = (*m_clusterCol[i][0]).Part();
1899 segment = (*m_clusterCol[i][0]).Segment();
1900 layer = (*m_clusterCol[i][0]).Layer();
1901
1902 FillCluster( part, segment, layer, size );
1903 m_ntClusterSize = size;
1904 m_clusterSizeTuple->write();
1905 }
1906 }
1907
1908 // Fill inc/noise hits
1909 log << MSG::INFO << "Fill inc/noise hits" << endreq;
1910 for( unsigned int i=0; i<m_nosHitCol.size(); i++ )
1911 {
1912 part = (*m_nosHitCol[i]).Part();
1913 segment = (*m_nosHitCol[i]).Segment();
1914 layer = (*m_nosHitCol[i]).Layer();
1915 strip = (*m_nosHitCol[i]).Strip();
1916
1917 FillNosHit( part, segment, layer, strip );
1918 m_record[part][segment][layer][strip][3] ++;
1919 }
1920
1921 // Event log
1922 m_ntEventId = m_currentEvent;
1923 m_ntEventTag = m_eventTag;
1924 m_ntDigiNum = m_digiCol.size();
1925 m_ntExpHitNum = m_expHitCol.size();
1926 m_ntEffHitNum = m_effHitCol.size();
1927 m_ntNosHitNum = m_nosHitCol.size();
1928 m_ntClusterNum = m_clusterCol.size();
1929
1930 // Reset mark collection
1931 for( unsigned int i=0; i<m_digiCol.size(); i++ ) {
1932 if( m_digiCol[i] != NULL ) delete m_digiCol[i];
1933 }
1934
1935 for( unsigned int i=0; i<m_expHitCol.size(); i++ ) {
1936 if( m_expHitCol[i] != NULL ) delete m_expHitCol[i];
1937 }
1938/*
1939 for( unsigned int i=0; i<m_effHitCol.size(); i++ ) {
1940 if( m_effHitCol[i] != NULL ) delete m_effHitCol[i];
1941 }
1942 log << MSG::INFO << m_effHitCol.size() << endreq;
1943*/
1944 for( unsigned int i=0; i<m_calHitCol.size(); i++ ) {
1945 if( m_calHitCol[i] != NULL ) delete m_calHitCol[i];
1946 }
1947
1948 if( m_digiCol.size() != 0 ) m_digiCol.clear();
1949 if( m_expHitCol.size() != 0 ) m_expHitCol.clear();
1950 if( m_calHitCol.size() != 0 ) m_calHitCol.clear();
1951 if( m_effHitCol.size() != 0 ) m_effHitCol.clear();
1952 if( m_nosHitCol.size() != 0 ) m_nosHitCol.clear();
1953 if( m_clusterCol.size() != 0 ) m_clusterCol.clear();
1954
1955
1956 for( int i=0; i<PART_MAX; i++ ) {
1957 for( int j=0; j<SEGMENT_MAX; j++ ) {
1958 if( m_segDigiCol[i][j].size() != 0 ) m_segDigiCol[i][j].clear();
1959 }
1960 }
1961
1962 m_evtEnd = clock();
1963
1964 m_ntEventTime = (double(m_evtEnd - m_evtBegin))/CLOCKS_PER_SEC;
1965 log << MSG::INFO << "Event time:\t" << m_ntEventTime << "s" << endreq;
1966 m_eventLogTuple->write();
1967
1968 return StatusCode::SUCCESS;
1969}
1970
1971StatusCode MucCalibMgr::FillDigi( int part, int segment, int layer, int strip )
1972{
1973 // Hit map for online check
1974 int tmpId;
1975
1976 if( (int)m_tag || m_dimuOnly==0 || (m_dimuOnly==1&&m_eventTag==1) )
1977 {
1978 if( part == BRID )
1979 {
1980 // According to layer
1981 if( layer%2 == 0 ) {
1982 if( segment > 4 ) tmpId = segment*48 + (47 - strip);
1983 else tmpId = segment*48 + strip;
1984 }
1985 else if( segment < 3 ) tmpId = segment*96 + strip;
1986 else tmpId = (segment-1)*96 + 112 + strip;
1987
1988 m_hHitMapBarrel_Lay[layer]->Fill(tmpId);
1989
1990 // According to segment
1991 if( segment != B_TOP ) {
1992 if( segment > 4 )
1993 tmpId = 48*((layer+1)/2) + 96*(layer/2) + ( ((layer+1)%2)*48 + (layer%2)*96 - strip - 1 );
1994 else tmpId = 48*((layer+1)/2) + 96*(layer/2) + strip;
1995 }
1996 else tmpId = 48*((layer+1)/2) + 112*(layer/2) + strip;
1997
1998 m_hHitMapBarrel_Seg[segment]->Fill(tmpId);
1999 }
2000 else if( part == EEID )
2001 {
2002 // According to layer
2003 tmpId = segment*64 + strip;
2004 m_hHitMapEndcap_Lay[0][layer]->Fill(tmpId);
2005 // According to segment
2006 tmpId = layer*64 + strip;
2007 m_hHitMapEndcap_Seg[0][segment]->Fill(tmpId);
2008 }
2009 else if( part == EWID )
2010 {
2011 // According to layer
2012 tmpId = segment*64 + strip;
2013 m_hHitMapEndcap_Lay[1][layer]->Fill(tmpId);
2014 // According to segment
2015 tmpId = layer*64 + strip;
2016 m_hHitMapEndcap_Seg[1][segment]->Fill(tmpId);
2017 }
2018 }
2019
2020 // Total units
2021 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2022 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
2023
2024 m_hStripFireMap[boxId]->Fill( strip );
2025 m_hStripFire->Fill( strId );
2026 m_hBoxFire->Fill( boxId );
2027 m_hLayerFire->Fill( layer );
2028 if( part==BRID ) m_hBrLayerFire->Fill( layer );
2029 else if( part== EEID ) m_hEcLayerFire->Fill( layer+1 );
2030 else m_hEcLayerFire->Fill( -(layer+1) );
2031
2032 return StatusCode::SUCCESS;
2033}
2034
2035StatusCode MucCalibMgr::FillExpHit( int part, int segment, int layer, int strip )
2036{
2037 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2038 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
2039
2040 m_hStripExpHitMap[boxId]->Fill( strip );
2041 m_hStripExpHit->Fill( strId );
2042 m_hBoxExpHit->Fill( boxId );
2043 m_hLayerExpHit->Fill( layer );
2044
2045 return StatusCode::SUCCESS;
2046}
2047
2048StatusCode MucCalibMgr::FillEffHit( int part, int segment, int layer, int strip )
2049{
2050 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2051 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
2052
2053 m_hStripEffHitMap[boxId]->Fill( strip );
2054 m_hStripEffHit->Fill( strId );
2055 m_hBoxEffHit->Fill( boxId );
2056 m_hLayerEffHit->Fill( layer );
2057
2058 return StatusCode::SUCCESS;
2059}
2060
2061StatusCode MucCalibMgr::FillNosHit( int part, int segment, int layer, int strip )
2062{
2063 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2064 int strId = m_ptrIdTr->GetStripId( part, segment, layer, strip );
2065
2066 m_hStripNosHitMap[boxId]->Fill( strip );
2067 m_hStripNosHit->Fill( strId );
2068 m_hBoxNosHit->Fill( boxId );
2069 m_hLayerNosHit->Fill( layer );
2070
2071 return StatusCode::SUCCESS;
2072}
2073
2074StatusCode MucCalibMgr::FillCluster( int part, int segment, int layer, int size )
2075{
2076 int boxId = m_ptrIdTr->GetBoxId( part, segment, layer );
2077
2078 m_hLayerCluster[layer]->Fill( size );
2079 m_hBoxCluster[boxId]->Fill( size );
2080
2081 return StatusCode::SUCCESS;
2082}
2083
2084
2085//-------------------------------------------------------------------------------------------------
2086//-------------------------------- Member functions for finalizing --------------------------------
2087//-------------------------------------------------------------------------------------------------
2088
2089StatusCode MucCalibMgr::AnalyseEffAndNoise()
2090{
2091 MsgStream log(msgSvc, "MucCalibMgr");
2092 log<< MSG::INFO << "Start efficiency and noise calibration" << endreq;
2093
2094 // Set time
2095 m_fTotalDAQTime = m_fTotalEvent/m_er;
2096 log<< MSG::INFO << "Total run time:\t" << m_fTotalDAQTime << endreq;
2097
2098 EffAndNoiseLV2();
2099 EffAndNoiseLV1();
2100 EffAndNoiseLV0();
2101
2102 if( m_usePad != 0 ) PadEff();
2103
2104 return StatusCode::SUCCESS;
2105}
2106
2107StatusCode MucCalibMgr::EffAndNoiseLV0()
2108{
2109 MsgStream log(msgSvc, "MucCalibMgr");
2110 log<< MSG::INFO << "Analyse layer efficiency and noise" << endreq;
2111
2112 int part, segment, layer, stripMax;
2113 part = segment = layer = stripMax = 0;
2114
2115 double digi, effHit, trackNum, nosHit, recHit;
2116 double eff, effErr, noise, nosRatio, nosRatioErr, cnt, cluster;
2117 eff = effErr = noise = nosRatio = nosRatioErr = cnt = cluster = 0.;
2118
2119 for( int i=0; i<LAYER_MAX; i++ )
2120 {
2121 digi = effHit = trackNum = nosHit = recHit = 0;
2122
2123 for( int j=0; j<PART_MAX; j++ )
2124 {
2125 for( int k=0; k<SEGMENT_MAX; k++)
2126 {
2127 stripMax = m_ptrIdTr->GetStripMax( j, k, i );
2128 for( int n=0; n<stripMax; n++ )
2129 {
2130 digi += m_record[j][k][i][n][0];
2131 trackNum += m_record[j][k][i][n][1];
2132 effHit += m_record[j][k][i][n][2];
2133 nosHit += m_record[j][k][i][n][3];
2134 recHit += m_record[j][k][i][n][4];
2135 }
2136 }
2137 }
2138
2139 // Efficiency
2140 if( trackNum == 0 ) {
2141 eff = effErr = 0.0;
2142 //eff = DEFAULT_EFF_VALUE;
2143 //effErr = DEFAULT_EFF_ERR;
2144 }
2145 else
2146 {
2147 eff = ( (double)effHit )/trackNum;
2148 effErr = sqrt( eff*(1-eff)/trackNum );
2149 m_fCalibLayerNum ++;
2150 }
2151
2152 // Noise
2153 if( m_layerResults[3][i] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
2154 noise = DEFAULT_NOS_VALUE;
2155 else {
2156 if( m_recMode == 2 ) noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW*m_layerResults[3][i]);
2157 else noise = (double)nosHit/(m_fTotalDAQTime * m_layerResults[3][i]);
2158 }
2159
2160 if( digi != 0 ) {
2161 nosRatio = ( (double)nosHit )/digi;
2162 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
2163 }
2164 else {
2165 nosRatio = DEFAULT_INC_VALUE;
2166 nosRatioErr = 0;
2167 }
2168
2169 // Counting rate
2170 if( m_recMode == 2 )
2171 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_layerResults[3][i]);
2172 else
2173 cnt = (double)digi/(m_fTotalDAQTime * m_layerResults[3][i]);
2174
2175 // Cluster
2176 cluster = m_hLayerCluster[i]->GetMean();
2177
2178 m_layerResults[0][ i ] = eff;
2179 m_layerResults[1][ i ] = effErr;
2180 m_layerResults[2][ i ] = noise;
2181 m_layerResults[4][ i ] = cluster;
2182 m_layerResults[5][ i ] = trackNum;
2183
2184 // Fill histogram
2185 m_hLayerEff->Fill( i, eff );
2186 m_hLayerEff->SetBinError( i+1, effErr );
2187 m_hLayerNosRatio->Fill( i, nosRatio );
2188 m_hLayerNosRatio->SetBinError( i+1, nosRatioErr );
2189 m_hLayerNos->Fill( i, noise );
2190 m_hLayerCnt->Fill( i, cnt );
2191
2192 // Add cluster histogram
2193 m_hLayerClusterCmp->Fill( i, cluster );
2194
2195 // Fill tree
2196 m_fLayerId = i;
2197 m_fLayerEff = eff;
2198 m_fLayerEffErr = effErr;
2199 m_fLayerTrkNum = trackNum;
2200 m_fLayerExpHit = recHit;
2201 m_fLayerEffHit = effHit;
2202 m_fLayerNosRatio = nosRatio;
2203 m_fLayerDigi = digi;
2204 m_fLayerNosHit = nosHit;
2205 m_fLayerNos = noise;
2206 m_fLayerCnt = cnt;
2207 m_tLayConst->Fill();
2208
2209 // Add cluster ntuple
2210 m_fLayerCluster = cluster;
2211
2212 } // End layer
2213
2214 return StatusCode::SUCCESS;
2215}
2216
2217StatusCode MucCalibMgr::EffAndNoiseLV1()
2218{
2219 MsgStream log(msgSvc, "MucCalibMgr");
2220 log<< MSG::INFO << "Analyse box efficiency and noise" << endreq;
2221
2222 int part, segment, layer, stripMax;
2223 part = segment = layer = stripMax = 0;
2224
2225 double digi, effHit, trackNum, nosHit, recHit;
2226 double eff, effErr, noise, nosRatio, nosRatioErr, cnt, cluster;
2227 eff = effErr = noise = nosRatio = nosRatioErr = cnt = cluster = 0.;
2228
2229 for( int i=0; i<BOX_MAX; i++ )
2230 {
2231 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
2232 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
2233
2234 digi = effHit = trackNum = nosHit = recHit = 0;
2235 for( int j=0; j<stripMax; j++ )
2236 {
2237 digi += m_record[part][segment][layer][j][0];
2238 trackNum += m_record[part][segment][layer][j][1];
2239 effHit += m_record[part][segment][layer][j][2];
2240 nosHit += m_record[part][segment][layer][j][3];
2241 recHit += m_record[part][segment][layer][j][4];
2242 }
2243
2244 // Efficiency
2245 if( trackNum == 0 ) {
2246 eff = effErr = 0.0;
2247 //eff = DEFAULT_EFF_VALUE;
2248 //effErr = DEFAULT_EFF_ERR;
2249 }
2250 else
2251 {
2252 eff = ( (double)effHit )/trackNum;
2253 effErr = sqrt( eff*(1-eff)/trackNum );
2254 m_fCalibBoxNum ++;
2255 }
2256
2257 // Noise
2258 if( m_boxResults[3][i] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
2259 noise = DEFAULT_NOS_VALUE;
2260 else {
2261 if( m_recMode == 2 )
2262 noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW * m_boxResults[3][i]);
2263 else
2264 noise = (double)nosHit/(m_fTotalDAQTime * m_boxResults[3][i]);
2265 }
2266
2267 if( digi != 0 ) {
2268 nosRatio = ( (double)nosHit )/digi;
2269 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
2270 }
2271 else {
2272 nosRatio = DEFAULT_INC_VALUE;
2273 nosRatioErr = 0;
2274 }
2275
2276 // Counting rate
2277 if( m_recMode == 2 )
2278 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_boxResults[3][i]);
2279 else
2280 cnt = (double)digi/(m_fTotalDAQTime * m_boxResults[3][i]);
2281
2282 // Cluster
2283 cluster = m_hBoxCluster[i]->GetMean();
2284
2285 m_boxResults[0][ i ] = eff;
2286 m_boxResults[1][ i ] = effErr;
2287 m_boxResults[2][ i ] = noise;
2288 m_boxResults[4][ i ] = cluster;
2289 m_boxResults[5][ i ] = trackNum;
2290
2291 // Fill histograms
2292 m_hBoxEff->Fill( i, eff );
2293 m_hBoxEff->SetBinError( i+1, effErr );
2294 m_hBoxNosRatio->Fill( i, nosRatio );
2295 m_hBoxNosRatio->SetBinError( i+1, nosRatioErr );
2296 m_hBoxNos->Fill( i, noise );
2297 m_hBoxCnt->Fill( i, cnt );
2298
2299 // add cluster histogram
2300 m_hBoxClusterCmp->Fill( i, cluster );
2301
2302 // Fill tree
2303 m_fBoxId = i;
2304 m_fBoxPart = part;
2305 m_fBoxSegment = segment;
2306 m_fBoxLayer = layer;
2307 m_fBoxEff = eff;
2308 m_fBoxEffErr = effErr;
2309 m_fBoxTrkNum = trackNum;
2310 m_fBoxExpHit = recHit;
2311 m_fBoxEffHit = effHit;
2312 m_fBoxNosRatio = nosRatio;
2313 m_fBoxDigi = digi;
2314 m_fBoxNosHit = nosHit;
2315 m_fBoxNos = noise;
2316 m_fBoxCnt = cnt;
2317 m_tBoxConst->Fill();
2318
2319 // add box cluster ntuple
2320 m_fBoxCluster = cluster;
2321
2322 }// End BOX_MAX
2323
2324 return StatusCode::SUCCESS;
2325}
2326
2327StatusCode MucCalibMgr::EffAndNoiseLV2()
2328{
2329 MsgStream log(msgSvc, "MucCalibMgr");
2330 log<< MSG::INFO << "Analyse strip efficiency and noise" << endreq;
2331
2332 int part, segment, layer, stripMax;
2333 part = segment = layer = stripMax = 0;
2334
2335 double digi, effHit, trackNum, nosHit, recHit;
2336 double eff, effErr, noise, nosRatio, nosRatioErr, cnt, cluster;
2337 eff = effErr = noise = nosRatio = nosRatioErr = cnt = cluster = 0.;
2338
2339 for( int i=0; i<BOX_MAX; i++ )
2340 {
2341 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
2342 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
2343
2344 for( int j=0; j<stripMax; j++ )
2345 {
2346 digi = m_record[part][segment][layer][j][0];
2347 trackNum = m_record[part][segment][layer][j][1];
2348 effHit = m_record[part][segment][layer][j][2];
2349 nosHit = m_record[part][segment][layer][j][3];
2350 recHit = m_record[part][segment][layer][j][4];
2351
2352 int stripId = m_ptrIdTr->GetStripId( part, segment, layer, j );
2353
2354 // Efficiency
2355 if( trackNum == 0 ) {
2356 eff = effErr = 0.0;
2357 //eff = DEFAULT_EFF_VALUE;
2358 //effErr = DEFAULT_EFF_ERR;
2359 }
2360 else
2361 {
2362 eff = ( (double)effHit )/trackNum;
2363 effErr = sqrt( eff*(1-eff)/trackNum );
2364 m_fCalibStripNum ++;
2365 }
2366
2367 // Noise
2368 if( m_stripResults[3][stripId] < LIMIT_CUT || m_fTotalDAQTime < LIMIT_CUT )
2369 noise = DEFAULT_NOS_VALUE;
2370 else {
2371 if( m_recMode == 2 )
2372 noise = (double)nosHit/(m_fTotalEvent*TRIGGER_WINDOW * m_stripResults[3][stripId]);
2373 else
2374 noise = (double)nosHit/(m_fTotalDAQTime * m_stripResults[3][stripId]);
2375 }
2376
2377 if( digi != 0 ) {
2378 nosRatio = ( (double)nosHit )/digi;
2379 nosRatioErr = sqrt( nosRatio*(1-nosRatio)/digi );
2380 }
2381 else {
2382 nosRatio = DEFAULT_INC_VALUE;
2383 nosRatioErr = 0.;
2384 }
2385
2386 // Counting rate
2387 if( m_recMode == 2 )
2388 cnt = (double)digi/(m_fTotalEvent*TRIGGER_WINDOW * m_stripResults[3][stripId]);
2389 else
2390 cnt = (double)digi/(m_fTotalDAQTime * m_stripResults[3][stripId]);
2391
2392 // Strip itself no clusters
2393 m_stripResults[0][ stripId ] = eff;
2394 m_stripResults[1][ stripId ] = effErr;
2395 m_stripResults[2][ stripId ] = noise;
2396 m_stripResults[5][ stripId ] = trackNum;
2397
2398 // Fill histotram
2399 m_hStripEffMap[i]->Fill( j, eff );
2400 m_hStripEffMap[i]->SetBinError( j+1, effErr );
2401 m_hStripEff->Fill( stripId, eff );
2402 m_hStripEff->SetBinError( stripId+1, effErr );
2403 m_hStripNosRatioMap[i]->Fill( j, nosRatio );
2404 m_hStripNosRatioMap[i]->SetBinError( j+1, nosRatioErr );
2405 m_hStripNosRatio->Fill( stripId, nosRatio );
2406 m_hStripNosRatio->SetBinError( stripId+1, nosRatioErr );
2407 m_hStripNos->Fill( stripId, noise );
2408 m_hStripCnt->Fill( stripId, cnt );
2409
2410 // Fill Tree
2411 m_fStripId = stripId;
2412 m_fStripPart = part;
2413 m_fStripSegment = segment;
2414 m_fStripLayer = layer;
2415 m_fStripEff = eff;
2416 m_fStripEffErr = effErr;
2417 m_fStripNosRatio = nosRatio;
2418 m_fStripDigi = digi;
2419 m_fStripNos = noise;
2420 m_fStripCnt = cnt;
2421 m_fStripEffHit = effHit;
2422 m_fStripExpHit = recHit;
2423 m_fStripNosHit = nosHit;
2424 m_fStripTrkNum = trackNum;
2425 m_tStrConst->Fill();
2426
2427 } // End stripMax
2428 } // End BOX_MAX
2429
2430 return StatusCode::SUCCESS;
2431}
2432
2433StatusCode MucCalibMgr::PadEff()
2434{
2435 MsgStream log(msgSvc, "MucCalibMgr");
2436
2437 int xBinStart, xBinEnd, yBinStart, yBinEnd;
2438 double expHit, firedHit, eff;
2439 eff = expHit = firedHit = 0.;
2440
2441 for( int i=0; i<PART_MAX; i++ )
2442 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
2443 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
2444 {
2445 xBinStart = m_h2DExpMap[i][j][k]->GetXaxis()->GetFirst();
2446 xBinEnd = m_h2DExpMap[i][j][k]->GetXaxis()->GetLast() + 1;
2447 yBinStart = m_h2DExpMap[i][j][k]->GetYaxis()->GetFirst();
2448 yBinEnd = m_h2DExpMap[i][j][k]->GetYaxis()->GetLast() + 1;
2449
2450 for( int xi = xBinStart; xi<xBinEnd; xi++ )
2451 for( int yi = yBinStart; yi<yBinEnd; yi++ )
2452 {
2453 expHit = m_h2DExpMap[i][j][k]->GetBinContent(xi, yi);
2454 firedHit = m_h2DHitMap[i][j][k]->GetBinContent(xi, yi);
2455
2456 if( expHit !=0 ) eff = firedHit / expHit;
2457 else eff = 0;
2458
2459 if( eff>1.0 )
2460 cout<<"Eff error:\t["<<i<<"\t"<<j<<"\t"<<k<<",\t"<<xi<<"\t"<<yi<<"]:\t"
2461 <<eff<<"\t,"<<firedHit<<" "<<expHit<<endl;
2462
2463 m_h2DEffMap[i][j][k]->SetBinContent(xi, yi, eff);
2464 }
2465 }
2466 return StatusCode::SUCCESS;
2467}
2468
2469StatusCode MucCalibMgr::AnalyseCluster()
2470{
2471 MsgStream log(msgSvc, "MucCalibMgr");
2472
2473 return StatusCode::SUCCESS;
2474}
2475
2476StatusCode MucCalibMgr::AnalyseRes()
2477{
2478 MsgStream log(msgSvc, "MucCalibMgr");
2479 log << MSG::INFO << "Analyse spacial resolution" << endreq;
2480
2481 double resMean, resMeanErr, resSigma, resSigmaErr;
2482 resMean = resMeanErr = resSigma = resSigmaErr = 0.;
2483
2484 for( int i=0; i<B_LAY_NUM; i++ )
2485 {
2486 m_hBarrelResDist[i]->Fit("gaus");
2487 resMean = m_hBarrelResDist[i]->GetFunction("gaus")->GetParameter("Mean");
2488 resSigma = m_hBarrelResDist[i]->GetFunction("gaus")->GetParameter("Sigma");
2489 resMeanErr = m_hBarrelResDist[i]->GetFunction("gaus")->GetParError(1);
2490 resSigmaErr = m_hBarrelResDist[i]->GetFunction("gaus")->GetParError(2);
2491
2492 m_hBarrelResComp[0]->Fill( i, resMean );
2493 m_hBarrelResComp[0]->SetBinError( i+1, resMeanErr );
2494 m_hBarrelResComp[1]->Fill( i, resSigma );
2495 m_hBarrelResComp[1]->SetBinError( i+1, resSigmaErr );
2496 }
2497
2498 for( int i=0; i<E_LAY_NUM; i++ )
2499 {
2500 m_hEndcapResDist[i]->Fit("gaus");
2501 resMean = m_hEndcapResDist[i]->GetFunction("gaus")->GetParameter("Mean");
2502 resSigma = m_hEndcapResDist[i]->GetFunction("gaus")->GetParameter("Sigma");
2503 resMeanErr = m_hEndcapResDist[i]->GetFunction("gaus")->GetParError(1);
2504 resSigmaErr = m_hEndcapResDist[i]->GetFunction("gaus")->GetParError(2);
2505
2506 m_hEndcapResComp[0]->Fill( i, resMean );
2507 m_hEndcapResComp[0]->SetBinError( i+1, resMeanErr );
2508 m_hEndcapResComp[1]->Fill( i, resSigma );
2509 m_hEndcapResComp[1]->SetBinError( i+1, resSigmaErr );
2510 }
2511
2512 return StatusCode::SUCCESS;
2513}
2514
2515StatusCode MucCalibMgr::SaveConst()
2516{
2517 MsgStream log(msgSvc, "MucCalibMgr");
2518 log << MSG::INFO << "Save calibration constants" << endreq;
2519
2520 // Save calibrated eff in graphes
2521 // LV0
2522 double layerXY[2][LAYER_MAX];
2523 double layerEXY[2][LAYER_MAX];
2524 for( int i=0; i<LAYER_MAX; i++ )
2525 {
2526 layerXY[0][i] = i;
2527 layerEXY[0][i] = 0;
2528 if( m_layerResults[5][i] >= 100*TRACK_THRESHOLD ) {
2529 layerXY[1][i] = m_layerResults[0][i];
2530 layerEXY[1][i] = m_layerResults[1][i];
2531 }
2532 else {
2533 layerXY[1][i] = 0;
2534 layerEXY[1][i] = 0;
2535 }
2536 }
2537 m_geLayerEff = new TGraphErrors(LAYER_MAX, layerXY[0], layerXY[1], layerEXY[0], layerEXY[1]);
2538 m_geLayerEff->SetMarkerStyle(25);
2539 m_geLayerEff->SetMarkerSize(0.5);
2540 m_cv[0] = new TCanvas("GoodLayerEff","Layer efficiency", 50, 50, 800, 600);
2541 m_cv[0]->SetFillColor(0);
2542 m_cv[0]->SetBorderMode(0);
2543 m_geLayerEff->Draw("AP");
2544 m_cv[0]->Write();
2545
2546 // LV1
2547 double boxXY[2][BOX_MAX];
2548 double boxEXY[2][BOX_MAX];
2549 for( int i=0; i<BOX_MAX; i++ )
2550 {
2551 boxXY[0][i] = i;
2552 boxEXY[0][i] = 0;
2553 if( m_boxResults[5][i] >= 10*TRACK_THRESHOLD ) {
2554 boxXY[1][i] = m_boxResults[0][i];
2555 boxEXY[1][i] = m_boxResults[1][i];
2556 }
2557 else {
2558 boxXY[1][i] = 0;
2559 boxEXY[1][i] = 0;
2560 }
2561 }
2562 m_geBoxEff = new TGraphErrors(BOX_MAX, boxXY[0], boxXY[1], boxEXY[0], boxEXY[1]);
2563 m_geBoxEff->SetMarkerStyle(25);
2564 m_geBoxEff->SetMarkerSize(0.5);
2565 m_cv[1] = new TCanvas("GoodBoxEff","Box efficiency", 75, 75, 800, 600);
2566 m_cv[1]->SetFillColor(0);
2567 m_cv[1]->SetBorderMode(0);
2568 m_geBoxEff->Draw("AP");
2569 m_cv[1]->Write();
2570
2571 // LV2
2572 double stripXY[2][STRIP_MAX];
2573 double stripEXY[2][STRIP_MAX];
2574 for( int i=0; i<STRIP_MAX; i++ )
2575 {
2576 stripXY[0][i] = i;
2577 stripEXY[0][i] = 0;
2578 if( m_stripResults[5][i] >= TRACK_THRESHOLD ) {
2579 stripXY[1][i] = m_stripResults[0][i];
2580 stripEXY[1][i] = m_stripResults[1][i];
2581 }
2582 else {
2583 stripXY[1][i] = 0;
2584 stripEXY[1][i] = 0;
2585 }
2586 }
2587 m_geStripEff = new TGraphErrors(STRIP_MAX, stripXY[0], stripXY[1], stripEXY[0], stripEXY[1]);
2588 m_geStripEff->SetMarkerStyle(25);
2589 m_geStripEff->SetMarkerSize(0.5);
2590 m_cv[2] = new TCanvas("GoodStripEff","Strip efficiency", 100, 100, 800, 600);
2591 m_cv[2]->SetFillColor(0);
2592 m_cv[2]->SetBorderMode(0);
2593 m_geStripEff->Draw("AP");
2594 m_cv[2]->Write();
2595
2596 // Save histograms
2597 for(int i=0; i<B_LAY_NUM; i++ )
2598 {
2599 m_hHitMapBarrel_Lay[i]->Write();
2600
2601 if( i<E_LAY_NUM ) {
2602 m_hHitMapEndcap_Lay[0][i]->Write();
2603 m_hHitMapEndcap_Lay[1][i]->Write();
2604 }
2605 }
2606
2607 for(int i=0; i<B_SEG_NUM; i++)
2608 {
2609 m_hHitMapBarrel_Seg[i]->Write();
2610
2611 if( i< E_SEG_NUM ) {
2612 m_hHitMapEndcap_Seg[0][i]->Write();
2613 m_hHitMapEndcap_Seg[1][i]->Write();
2614 }
2615 }
2616 m_hTrackDistance->Fit("gaus");
2617 m_hTrackPosPhiDiff->Fit("gaus");
2618 m_hTrackPosThetaDiff->Fit("gaus");
2619 m_hTrackMomPhiDiff->Fit("gaus");
2620 m_hTrackMomThetaDiff->Fit("gaus");
2621
2622 m_hTrackDistance->Write();
2623 m_hTrackPosPhiDiff->Write();
2624 m_hTrackPosThetaDiff->Write();
2625 m_hTrackMomPhiDiff->Write();
2626 m_hTrackMomThetaDiff->Write();
2627
2628 for( int i=0; i<B_LAY_NUM; i++ ) m_hBarrelResDist[i]->Write();
2629 for( int i=0; i<E_LAY_NUM; i++ ) m_hEndcapResDist[i]->Write();
2630
2631 m_hBarrelResComp[0]->Write();
2632 m_hBarrelResComp[1]->Write();
2633 m_hEndcapResComp[0]->Write();
2634 m_hEndcapResComp[1]->Write();
2635
2636 if( m_usePad != 0 ) m_histArray->Write();
2637
2638 for( int i=0; i<BOX_MAX; i++ )
2639 {
2640 m_hStripFireMap[ i ]->Write();
2641 m_hStripExpHitMap[ i ]->Write();
2642 m_hStripEffHitMap[ i ]->Write();
2643 m_hStripNosHitMap[ i ]->Write();
2644 m_hStripEffMap[ i ]->Write();
2645 m_hStripNosRatioMap[ i ]->Write();
2646 }
2647 m_hStripFire->Write();
2648 m_hStripExpHit->Write();
2649 m_hStripEffHit->Write();
2650 m_hStripNosHit->Write();
2651 m_hStripEff->Write();
2652 m_hStripArea->Write();
2653 m_hStripNos->Write();
2654 m_hStripNosRatio->Write();
2655 log << MSG::INFO << "Save LV2 histograms done!" << endreq;
2656
2657 m_hBoxFire->Write();
2658 m_hBoxExpHit->Write();
2659 m_hBoxEffHit->Write();
2660 m_hBoxNosHit->Write();
2661 m_hBoxEff->Write();
2662 m_hBoxArea->Write();
2663 m_hBoxNos->Write();
2664 m_hBoxNosRatio->Write();
2665 log << MSG::INFO << "Save LV1 histograms done!" << endreq;
2666
2667 m_hBrLayerFire->Write();
2668 m_hEcLayerFire->Write();
2669 m_hLayerFire->Write();
2670 m_hLayerExpHit->Write();
2671 m_hLayerEffHit->Write();
2672 m_hLayerNosHit->Write();
2673 m_hLayerEff->Write();
2674 m_hLayerArea->Write();
2675 m_hLayerNos->Write();
2676 m_hLayerNosRatio->Write();
2677
2678 for( int i=0; i<LAYER_MAX; i++ ) m_hLayerCluster[i]->Write();
2679 for( int i=0; i<BOX_MAX; i++ ) m_hBoxCluster[i]->Write();
2680 m_hLayerClusterCmp->Write();
2681 m_hBoxClusterCmp->Write();
2682
2683 log << MSG::INFO << "Save histograms done!" << endreq;
2684
2685 // Save trees
2686 m_fLayerCoverage = 100*(double)m_fCalibLayerNum/LAYER_MAX;
2687 m_fBoxCoverage = 100*(double)m_fCalibBoxNum/BOX_MAX;
2688 m_fStripCoverage = 100*(double)m_fCalibStripNum/STRIP_MAX;
2689
2690 long digi_num, trk_num, eff_hit, nos_hit, exp_hit;
2691 m_tStatLog->Branch("digi_num", &digi_num, "digi_num/I");
2692 m_tStatLog->Branch("trk_num", &trk_num, "trk_num/I");
2693 m_tStatLog->Branch("eff_hit", &eff_hit, "eff_hit/I");
2694 m_tStatLog->Branch("nos_hit", &nos_hit, "nos_hit/I");
2695 m_tStatLog->Branch("exp_hit", &exp_hit, "exp_hit/I");
2696
2697 int stripMax;
2698 for( int i=0; i<PART_MAX; i++ )
2699 for( int j=0; j<((i==BRID)?B_SEG_NUM:E_SEG_NUM); j++ )
2700 for( int k=0; k<((i==BRID)?B_LAY_NUM:E_LAY_NUM); k++ )
2701 {
2702 stripMax = m_ptrIdTr->GetStripMax( i, j, k );
2703 for( int n=0; n<stripMax; n++ )
2704 {
2705 digi_num = m_record[i][j][k][n][0];
2706 trk_num = m_record[i][j][k][n][1];
2707 eff_hit = m_record[i][j][k][n][2];
2708 nos_hit = m_record[i][j][k][n][3];
2709 exp_hit = m_record[i][j][k][n][4];
2710 m_tStatLog->Fill();
2711 }
2712 }
2713
2714 m_jobFinish = clock();
2715 m_fTotalJobTime = (double)(m_jobFinish - m_jobStart)/CLOCKS_PER_SEC;
2716
2717 m_tJobLog->Fill();
2718 m_tJobLog->Write();
2719 m_tStatLog->Write();
2720 m_tLayConst->Write();
2721 m_tBoxConst->Write();
2722 m_tStrConst->Write();
2723
2724 // Close cluster output file
2725 if( m_fdata != NULL ) m_fdata->close();
2726
2727 return StatusCode::SUCCESS;
2728}
2729
2730StatusCode MucCalibMgr::EndRun()
2731{
2732 MsgStream log(msgSvc, "MucCalibMgr");
2733 log << MSG::INFO << endreq;
2734 log << MSG::INFO << "Total events : " << m_fTotalEvent << endreq;
2735 log << MSG::INFO << "Total digis : " << m_fTotalDigi << endreq;
2736 log << MSG::INFO << "Total rec hits : " << m_fTotalExpHit << endreq;
2737 log << MSG::INFO << "Total eff hits : " << m_fTotalEffHit << endreq;
2738 log << MSG::INFO << "Total inc hits : " << m_fTotalNosHit << endreq;
2739 log << MSG::INFO << "Total clusters : " << m_fTotalClstNum<< endreq;
2740
2741 log << MSG::INFO << "Strip calibrated percentage: "
2742 << 100*(double)m_fCalibStripNum/STRIP_MAX << "%" << endreq;
2743 log << MSG::INFO << "Box calibrated percentage: "
2744 << 100*(double)m_fCalibBoxNum/BOX_MAX << "%" << endreq;
2745 log << MSG::INFO << "Layer calibrated percentage: "
2746 << 100*(double)m_fCalibLayerNum/LAYER_MAX << "%" << endreq;
2747
2748 log << MSG::INFO << "Calibration run successfully" << endreq << endreq;
2749
2750 return StatusCode::SUCCESS;
2751}
2752
2753// END
Write
curve Write()
title
titledef title[20]
Definition: EvtPhokharaDef.hh:58
STRIP_MAX
const int STRIP_MAX
Definition: MucMappingAlg.h:24
mark_col
vector< MucMark * > mark_col
Definition: MucMark.h:17
PI
const double PI
Definition: PipiJpsi.cxx:55
RecMucTrackCol
ObjectVector< RecMucTrack > RecMucTrackCol
Definition: RecMucTrack.h:394
Segment
#define Segment
Definition: TTrackBase.h:31
theta2
@ theta2
Definition: TrkKalDeriv.h:24
theta1
@ theta1
Definition: TrkKalDeriv.h:24
MucCalibMgr::ClearHistoLV2
StatusCode ClearHistoLV2()
Definition: MucCalibMgr.cxx:279
MucCalibMgr::InitResHisto
StatusCode InitResHisto()
Definition: MucCalibMgr.cxx:869
MucCalibMgr::FillEffHit
StatusCode FillEffHit(int part, int segment, int layer, int strip)
Definition: MucCalibMgr.cxx:2048
MucCalibMgr::DimuSelect
StatusCode DimuSelect()
Definition: MucCalibMgr.cxx:1033
MucCalibMgr::InitOnlineHisto
StatusCode InitOnlineHisto()
Definition: MucCalibMgr.cxx:602
MucCalibMgr::Init2DEffHisto
StatusCode Init2DEffHisto()
Definition: MucCalibMgr.cxx:791
MucCalibMgr::SaveConst
StatusCode SaveConst()
Definition: MucCalibMgr.cxx:2515
MucCalibMgr::FillExpHit
StatusCode FillExpHit(int part, int segment, int layer, int strip)
Definition: MucCalibMgr.cxx:2035
MucCalibMgr::ClearResHisto
StatusCode ClearResHisto()
Definition: MucCalibMgr.cxx:338
MucCalibMgr::ClearHistoLV0
StatusCode ClearHistoLV0()
Definition: MucCalibMgr.cxx:245
MucCalibMgr::ClearClusterHisto
StatusCode ClearClusterHisto()
Definition: MucCalibMgr.cxx:321
MucCalibMgr::AnalyseCluster
StatusCode AnalyseCluster()
Definition: MucCalibMgr.cxx:2469
MucCalibMgr::ClearHistoLV1
StatusCode ClearHistoLV1()
Definition: MucCalibMgr.cxx:263
MucCalibMgr::PadEff
StatusCode PadEff()
Definition: MucCalibMgr.cxx:2433
MucCalibMgr::InitHistoLV1
StatusCode InitHistoLV1()
Definition: MucCalibMgr.cxx:724
MucCalibMgr::InitHistoLV2
StatusCode InitHistoLV2()
Definition: MucCalibMgr.cxx:740
MucCalibMgr::InitClusterHisto
StatusCode InitClusterHisto()
Definition: MucCalibMgr.cxx:839
MucCalibMgr::FillCluster
StatusCode FillCluster(int part, int segment, int layer, int size)
Definition: MucCalibMgr.cxx:2074
MucCalibMgr::ClearOnlineHisto
StatusCode ClearOnlineHisto()
Definition: MucCalibMgr.cxx:224
MucCalibMgr::InitHisto
StatusCode InitHisto()
Definition: MucCalibMgr.cxx:576
MucCalibMgr::FillEvent
StatusCode FillEvent()
Definition: MucCalibMgr.cxx:1841
MucCalibMgr::EffAndNoiseLV2
StatusCode EffAndNoiseLV2()
Definition: MucCalibMgr.cxx:2327
MucCalibMgr::ntupleSvc
INTupleSvc * ntupleSvc
Definition: MucCalibMgr.h:100
MucCalibMgr::MucCalibMgr
MucCalibMgr(std::vector< double > jobInfo, std::vector< int > configInfo, std::string outputFile)
Definition: MucCalibMgr.cxx:65
MucCalibMgr::InitHistoLV0
StatusCode InitHistoLV0()
Definition: MucCalibMgr.cxx:705
MucCalibMgr::eventSvc
IDataProviderSvc * eventSvc
Definition: MucCalibMgr.h:101
MucCalibMgr::EndRun
StatusCode EndRun()
Definition: MucCalibMgr.cxx:2730
MucCalibMgr::ReadEvent
StatusCode ReadEvent()
Definition: MucCalibMgr.cxx:1154
MucCalibMgr::EffAndNoiseLV1
StatusCode EffAndNoiseLV1()
Definition: MucCalibMgr.cxx:2217
MucCalibMgr::AnalyseRes
StatusCode AnalyseRes()
Definition: MucCalibMgr.cxx:2476
MucCalibMgr::msgSvc
IMessageSvc * msgSvc
Definition: MucCalibMgr.h:99
MucCalibMgr::CheckEvent
StatusCode CheckEvent()
Definition: MucCalibMgr.cxx:1786
MucCalibMgr::FillNosHit
StatusCode FillNosHit(int part, int segment, int layer, int strip)
Definition: MucCalibMgr.cxx:2061
MucCalibMgr::EffAndNoiseLV0
StatusCode EffAndNoiseLV0()
Definition: MucCalibMgr.cxx:2107
MucCalibMgr::InitConstTree
StatusCode InitConstTree()
Definition: MucCalibMgr.cxx:897
MucCalibMgr::ClearConstTree
StatusCode ClearConstTree()
Definition: MucCalibMgr.cxx:356
MucCalibMgr::InitArea
StatusCode InitArea()
Definition: MucCalibMgr.cxx:992
MucCalibMgr::AnalyseEffAndNoise
StatusCode AnalyseEffAndNoise()
Definition: MucCalibMgr.cxx:2089
MucCalibMgr::InitNtuple
StatusCode InitNtuple()
Definition: MucCalibMgr.cxx:374
MucCalibMgr::Clear2DEffHisto
StatusCode Clear2DEffHisto()
Definition: MucCalibMgr.cxx:305
MucCalibMgr::FillDigi
StatusCode FillDigi(int part, int segment, int layer, int strip)
Definition: MucCalibMgr.cxx:1971
MucEntityCal::GetArea
double GetArea()
Definition: MucEntityCal.cxx:95
MucID::barrel_ec
static int barrel_ec(const Identifier &id)
Values of different levels.
Definition: MucID.cxx:41
MucID::layer
static int layer(const Identifier &id)
Definition: MucID.cxx:61
MucID::channel
static int channel(const Identifier &id)
Definition: MucID.cxx:71
MucID::segment
static int segment(const Identifier &id)
Definition: MucID.cxx:51
MucIdTransform::GetStripId
int GetStripId(int part, int segment, int layer, int subid)
Definition: Muc/MucCalibAlg/MucCalibAlg-00-02-16/src/MucIdTransform.cxx:67
MucIdTransform::SetBoxPos
bool SetBoxPos(int boxid, int *part, int *segment, int *layer)
Definition: Muc/MucCalibAlg/MucCalibAlg-00-02-16/src/MucIdTransform.cxx:100
MucIdTransform::GetStripMax
int GetStripMax(int part, int segment, int layer)
Definition: Muc/MucCalibAlg/MucCalibAlg-00-02-16/src/MucIdTransform.cxx:38
MucIdTransform::GetBoxId
int GetBoxId(int part, int segment, int layer)
Definition: Muc/MucCalibAlg/MucCalibAlg-00-02-16/src/MucIdTransform.cxx:51
MucMark::IsInSegWith
bool IsInSegWith(MucMark &other)
Definition: MucMark.cxx:144
MucMark::IsInCol
int IsInCol(int part, int segment, int layer, int strip, mark_col &aCol)
Definition: MucMark.cxx:100
MucRecHit::GetIntersectZ
float GetIntersectZ() const
Definition: MucRecHit.h:110
MucRecHit::Seg
int Seg() const
Get Seg.
Definition: MucRecHit.h:74
MucRecHit::GetIntersectY
float GetIntersectY() const
Definition: MucRecHit.h:109
MucRecHit::GetIntersectX
float GetIntersectX() const
Definition: MucRecHit.h:108
MucRecHit::Part
int Part() const
Get Part.
Definition: MucRecHit.h:71
MucRecHit::Gap
int Gap() const
Get Gap.
Definition: MucRecHit.h:77
MucRecHit::Strip
int Strip() const
Get Strip.
Definition: MucRecHit.h:80
MucRecHit::GetPadID
int GetPadID() const
Definition: MucRecHit.h:107
MucRecHit::HitIsSeed
int HitIsSeed() const
Definition: MucRecHit.h:100
TofHitStatus::is_cluster
bool is_cluster() const
Definition: TofHitStatus.h:25
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:134
sprintf
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)
Event
Definition: Event.h:21