BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DQA_MUC.cxx
Go to the documentation of this file.
1#include <vector>
2
3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/SmartDataPtr.h"
5#include "GaudiKernel/PropertyMgr.h"
6#include "GaudiKernel/Bootstrap.h"
7#include "GaudiKernel/Algorithm.h"
8
9#include "GaudiKernel/INTupleSvc.h"
10#include "GaudiKernel/NTuple.h"
11#include "GaudiKernel/ITHistSvc.h"
12
13#include "CLHEP/Vector/ThreeVector.h"
14#include "CLHEP/Vector/LorentzVector.h"
15
17#include "EventModel/Event.h"
18
23
25#include "MucRawEvent/MucDigi.h"
29
30#include "DQAEvent/DQAEvent.h"
31#include "DQA_MUC/DQA_MUC.h"
32
33using CLHEP::Hep3Vector;
34/////////////////////////////////////////////////////////////////////////////
35DECLARE_COMPONENT(DQA_MUC)
36DQA_MUC::DQA_MUC(const std::string& name, ISvcLocator* pSvcLocator) :
37 Algorithm(name, pSvcLocator) {
38
39 //Declare the properties
40 declareProperty("EffWindow", m_effWindow = 4);
41 declareProperty("DstFileOnly", m_dstFileOnly = true);
42}
43
44// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
46 MsgStream log(msgSvc(), name());
47
48 log << MSG::INFO << "in initialize()" << endmsg;
49 StatusCode status;
50
51 // Call Histogram service
52 if(service("THistSvc", m_thsvc).isFailure()) {
53 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
54 return StatusCode::FAILURE;
55 }
56
57 char name[60];
58 char title[60];
59
60 // Resolution histograms
61 for( int i=0; i<B_LAY_NUM; i++ )
62 {
63 sprintf( name, "BrResDist_All_L%d", i );
64 sprintf( title, "Barrel spacial resolution in L%d",i );
65 m_hBrResDist[i][0] = new TH1F(name,title, 200, -100, 100 );
66 sprintf( name, "/DQAHist/MUC/BrResDist_All_L%d", i );
67 if(m_thsvc->regHist(name, m_hBrResDist[i][0]).isFailure())
68 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
69
70 sprintf( name, "BrResDist_Dimu_L%d", i );
71 m_hBrResDist[i][1] = new TH1F(name,title, 200, -100, 100 );
72 sprintf( name, "/DQAHist/MUC/BrResDist_Dimu_L%d", i );
73 if(m_thsvc->regHist(name, m_hBrResDist[i][1]).isFailure())
74 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
75 }
76
77 for( int i=0; i<E_LAY_NUM; i++ )
78 {
79 sprintf( name, "EcResDist_All_L%d", i );
80 sprintf( title, "Endcap spacial resolution in L%d",i );
81 m_hEcResDist[i][0] = new TH1F(name,title, 200, -100, 100 );
82 sprintf( name, "/DQAHist/MUC/EcResDist_All_L%d", i );
83 if(m_thsvc->regHist(name, m_hEcResDist[i][0]).isFailure())
84 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
85
86 sprintf( name, "EcResDist_Dimu_L%d", i );
87 m_hEcResDist[i][1] = new TH1F(name,title, 200, -100, 100 );
88 sprintf( name, "/DQAHist/MUC/EcResDist_Dimu_L%d", i );
89 if(m_thsvc->regHist(name, m_hEcResDist[i][1]).isFailure())
90 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
91 }
92
93 //m_hBrRes[0] = new TH1F("BrRes_All", "BrResAll", B_LAY_NUM+1, -0.5, B_LAY_NUM+0.5 );
94 //m_hBrRes[0]->GetXaxis()->SetTitle("Layer id");
95 //m_hBrRes[0]->GetXaxis()->CenterTitle();
96 //m_hBrRes[0]->SetStats(0);
97 //m_hBrRes[1] = new TH1F("BrRes_Dimu", "BrResDimu", B_LAY_NUM+1, -0.5, B_LAY_NUM+0.5 );
98 //m_hBrRes[1]->GetXaxis()->SetTitle("Layer id");
99 //m_hBrRes[1]->GetXaxis()->CenterTitle();
100 //m_hBrRes[1]->SetStats(0);
101 //
102 //m_hEcRes[0] = new TH1F("EcRes_All", "EcResAll", E_LAY_NUM+1, -0.5, E_LAY_NUM+0.5 );
103 //m_hEcRes[0]->GetXaxis()->SetTitle("Layer id");
104 //m_hEcRes[0]->GetXaxis()->CenterTitle();
105 //m_hEcRes[0]->SetStats(0);
106 //m_hEcRes[1] = new TH1F("EcRes_Dimu", "EcResDimu", E_LAY_NUM+1, -0.5, E_LAY_NUM+0.5 );
107 //m_hEcRes[1]->GetXaxis()->SetTitle("Layer id");
108 //m_hEcRes[1]->GetXaxis()->CenterTitle();
109 //m_hEcRes[1]->SetStats(0);
110
111 //if(m_thsvc->regHist("/DQAHist/MUC/BrRes_All", m_hBrRes[0]).isFailure())
112 //{ log << MSG::ERROR << "Couldn't register BrRes_All" << endreq; }
113 //if(m_thsvc->regHist("/DQAHist/MUC/BrRes_Dimu", m_hBrRes[1]).isFailure())
114 //{ log << MSG::ERROR << "Couldn't register BrRes_Dimu" << endreq; }
115 //if(m_thsvc->regHist("/DQAHist/MUC/EcRes_All", m_hEcRes[0]).isFailure())
116 //{ log << MSG::ERROR << "Couldn't register EcRes_All" << endreq; }
117 //if(m_thsvc->regHist("/DQAHist/MUC/EcRes_Dimu", m_hEcRes[1]).isFailure())
118 //{ log << MSG::ERROR << "Couldn't register BrRes_Dimu" << endreq; }
119
120 // Efficiency histograms, j=0: numberator, j=1; denominator
121 for(int i=0; i<TAGN; i++)
122 {
123 for(int j=0; j<2; j++)
124 {
125 for(int k=0; k<LVLN; k++)
126 {
127 sprintf( name, "%sEff_%s_%s", LNAME[k], ENAME[i], HTYPE[j] );
128 sprintf( title, "%s efficiency",LNAME[k] );
129
130 if(k==0) {
131 m_hEff[k][i][j] = new TH1F(name,title, LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
132 m_hEff[k][i][j]->GetXaxis()->SetTitle("Layer id");
133 }else if(k==1) {
134 m_hEff[k][i][j] = new TH1F(name,title, BOX_MAX+1, -0.5, BOX_MAX+0.5);
135 m_hEff[k][i][j]->GetXaxis()->SetTitle("Box id [EE:0~31, BR:32~103, WE:104~135]");
136 }else {
137 m_hEff[k][i][j] = new TH1F(name,title, STRIP_MAX+1, -0.5, STRIP_MAX+0.5);
138 m_hEff[k][i][j]->GetXaxis()->SetTitle("Strip id");
139 }
140
141 m_hEff[k][i][j]->GetXaxis()->CenterTitle();
142 //m_hEff[k][i][j]->SetStats(0);
143
144 sprintf( name, "/DQAHist/MUC/%sEff_%s_%s", LNAME[k], ENAME[i], HTYPE[j] );
145 if(m_thsvc->regHist(name, m_hEff[k][i][j]).isFailure())
146 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
147 }
148
149 sprintf( name, "%sNosRatio_%s_%s", LNAME[1], ENAME[i], HTYPE[j] );
150 sprintf( title, "%s noise ratio",LNAME[1] );
151 m_hNosRatio[i][j] = new TH1F(name, title, BOX_MAX+1, -0.5, BOX_MAX+0.5);
152 m_hNosRatio[i][j]->GetXaxis()->SetTitle("Box id [EE:0~31, BR:32~103, WE:104~135]");
153 m_hNosRatio[i][j]->GetXaxis()->CenterTitle();
154 //m_hNosRatio[i][j]->GetYaxis()->SetRangeUser(0., 0.3);
155 //m_hNosRatio[i][j]->SetStats(0);
156
157 sprintf( name, "/DQAHist/MUC/%sNosRatio_%s_%s", LNAME[1], ENAME[i], HTYPE[j] );
158 if(m_thsvc->regHist(name, m_hNosRatio[i][j]).isFailure())
159 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
160 }
161
162 // Costheta and phi
163 sprintf( name, "Costheta_%s", ENAME[i]);
164 sprintf( title, "Costheta");
165 m_hCostheta[i] = new TH1F(name, title, 360, -1.0, 1.0);
166 m_hCostheta[i]->GetXaxis()->SetTitle("cos#theta");
167 m_hCostheta[i]->GetXaxis()->CenterTitle();
168
169 sprintf( name, "Phi_%s", ENAME[i]);
170 sprintf( title, "Phi");
171 m_hPhi[i] = new TH1F(name, title, 360, -PI, PI);
172 m_hPhi[i]->GetXaxis()->SetTitle("#phi");
173 m_hPhi[i]->GetXaxis()->CenterTitle();
174
175 sprintf( name, "/DQAHist/MUC/Costheta_%s", ENAME[i] );
176 if(m_thsvc->regHist(name, m_hCostheta[i]).isFailure())
177 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
178 sprintf( name, "/DQAHist/MUC/Phi_%s", ENAME[i] );
179 if(m_thsvc->regHist(name, m_hPhi[i]).isFailure())
180 { log << MSG::ERROR << "Couldn't register Phi_All" << endreq; }
181 }
182
183 //m_hStripEff[0] = new TH1F("StripEff_All", "Strip efficiency", 101, -0.005, 1.005 );
184 //m_hStripEff[0]->GetXaxis()->SetTitle("Efficiency");
185 //m_hStripEff[0]->GetXaxis()->CenterTitle();
186 //m_hStripEff[1] = new TH1F("StripEff_Dimu", "Strip efficiency", 101, -0.005, 1.005 );
187 //m_hStripEff[1]->GetXaxis()->SetTitle("Efficiency");
188 //m_hStripEff[1]->GetXaxis()->CenterTitle();
189
190 // Parameters
191 for( int i=0; i<PART_MAX; i++ )
192 for( int j=0; j<SEGMENT_MAX; j++ )
193 for( int k=0; k<LAYER_MAX; k++ )
194 for( int n=0; n<STRIP_INBOX_MAX; n++ )
195 {
196 m_recordAll[i][j][k][n][0] = 0;
197 m_recordAll[i][j][k][n][1] = 0;
198 m_recordAll[i][j][k][n][2] = 0;
199 m_recordAll[i][j][k][n][3] = 0;
200 m_recordAll[i][j][k][n][4] = 0;
201
202 m_recordDimu[i][j][k][n][0] = 0;
203 m_recordDimu[i][j][k][n][1] = 0;
204 m_recordDimu[i][j][k][n][2] = 0;
205 m_recordDimu[i][j][k][n][3] = 0;
206 m_recordDimu[i][j][k][n][4] = 0;
207 }
208
209 m_ptrMucMark = new MucMark(0,0,0,0);
210 m_ptrIdTr = new MucIdTransform();
211
212 m_digiCol.resize(0);
213 m_expHitCol.resize(0);
214 m_effHitCol.resize(0);
215 m_nosHitCol.resize(0);
216 m_clusterCol.resize(0);
217
218 for( int i=0; i<PART_MAX; i++ )
219 for( int j=0; j<SEGMENT_MAX; j++ ) {
220 m_segDigiCol[i][j].resize(0);
221 }
222
223 log << MSG::INFO << "Initialize done!" <<endmsg;
224 return StatusCode::SUCCESS;
225}
226
227// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
228StatusCode DQA_MUC::execute() {
229
230 MsgStream log(msgSvc(), name());
231 log << MSG::INFO << "in execute()" << endreq;
232
233 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
234 m_run = eventHeader->runNumber();
235 m_event = eventHeader->eventNumber();
236 log << MSG::DEBUG <<"Run " << m_run << "\tEvent " << m_event <<endreq;
237
238 SmartDataPtr<DQAEvent::DQAEvent> dqaevt(eventSvc(), "/Event/DQATag");
239 if ( dqaevt ) {
240 log << MSG::INFO << "success get DQAEvent" << endreq;
241 } else {
242 log << MSG::ERROR << "Error accessing DQAEvent" << endreq;
243 return StatusCode::FAILURE;
244 }
245 log << MSG::DEBUG << "DQA event tag = " << dqaevt->EventTag() << endreq;
246
247 int part, segment, layer, strip;
248 char name[100];
249 bool isDimu = false;
250 if ( dqaevt->Dimu() ) isDimu = true;
251 log << MSG::INFO << "DQADimuTag:\t" << dqaevt->Dimu() << endreq;
252
253 //---> Retrieve MUC digi
254 if(!m_dstFileOnly)
255 {
256 log << MSG::INFO << "Retrieve digis" << endreq;
257 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
258 if(!mucDigiCol) {
259 log << MSG::FATAL << "Could not find MUC digi" << endreq;
260 return( StatusCode::FAILURE);
261 }
262
263 Identifier mucId;
264 MucDigiCol::iterator digiIter = mucDigiCol->begin();
265 int eventDigi = 0;
266 for ( int digiId =0; digiIter != mucDigiCol->end(); digiIter++, digiId++ )
267 {
268 mucId = (*digiIter)->identify();
269 part = MucID::barrel_ec(mucId); segment = MucID::segment(mucId);
270 layer = MucID::layer(mucId); strip = MucID::channel(mucId);
271 eventDigi ++;
272
273 // Add digi
274 MucMark *aMark = new MucMark( part, segment, layer, strip );
275 m_digiCol.push_back( aMark );
276 m_segDigiCol[part][segment].push_back( aMark );
277 m_recordAll[part][segment][layer][strip][0] ++;
278 if( isDimu ) m_recordDimu[part][segment][layer][strip][0] ++;
279 }
280 log << MSG::INFO << "Total digis in this event: " << eventDigi << endreq;
281 if( eventDigi > 500) cout <<"Run:\t"<< m_run <<"\tEvent:\t"<< m_event << "\tdigits inflation:\t" << eventDigi << endreq;
282 }
283 // Search cluster in digis
284 m_clusterCol = (*m_ptrMucMark).CreateClusterCol(1, m_digiCol );
285 log << MSG::INFO << "Total clusters in this event: " << m_clusterCol.size() << endreq;
286 //<--- End retrieve digis
287
288 //---> Retrieve rec tracks
289 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
290 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc(),"/Event/Recon/RecMucTrackCol");
291 if (!mucTrackCol) {
292// log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
293 log << MSG::WARNING << "Could not find RecMucTrackCol" << endreq; // until BOSS 7.0.3, if not MUC track, no RecMucTrackCol either.
294 return( StatusCode::SUCCESS);
295 }
296
297 RecMucTrackCol *aRecMucTrackCol = mucTrackCol;
298 if (aRecMucTrackCol->size() < 1) {
299 log << MSG::INFO << "No MUC tracks in this event" << endreq;
300 return StatusCode::SUCCESS;
301 }
302 log << MSG::INFO << "Total tracks of this event: " << aRecMucTrackCol->size() << endreq;
303
304 int trackHitNum, expectedHitNum, segNum, lastLayerBR, lastLayerEC;
305 int layerPassNum[3], passMax[TRACK_SEG_MAX][2];
306 bool firedLay[TRACK_SEG_MAX][LAYER_MAX];
307 double costheta, phi;
308 TH1* htmp(0);
309 trackHitNum = expectedHitNum = segNum = lastLayerBR = lastLayerEC = 0;
310 layerPassNum[0] = layerPassNum[1] = layerPassNum[2] = 0;
311 for( int segi=0; segi<TRACK_SEG_MAX; segi++ ) {
312 passMax[segi][0] = passMax[segi][1] = 0;
313 for( int layi=0; layi<LAYER_MAX; layi++ ) firedLay[segi][layi] = 0;
314 }
315
316 vector<MucRecHit*> mucRawHitCol;
317 vector<MucRecHit*> mucExpHitCol;
318 RecMucTrackCol::iterator trackIter = mucTrackCol->begin();
319 for (int trackId = 0; trackIter != mucTrackCol->end(); trackIter++, trackId++)
320 {
321 trackHitNum = (*trackIter)->numHits();
322
323 if( trackHitNum == 0 ) {
324 log << MSG::INFO << "Track " << trackId << " no hits" << endreq;
325 continue;
326 }
327
328 lastLayerBR = (*trackIter)->brLastLayer();
329 lastLayerEC = (*trackIter)->ecLastLayer();
330 // First fit position in MUC
331 CLHEP::Hep3Vector a3Vector((*trackIter)->xPos(),(*trackIter)->yPos(),(*trackIter)->zPos());
332 costheta = a3Vector.cosTheta();
333 phi = a3Vector.phi();
334 if(m_thsvc->getHist("/DQAHist/MUC/Costheta_All", htmp).isSuccess()) {
335 htmp->Fill( costheta );
336 } else {
337 log << MSG::ERROR << "Fail to retrieve Costheta_All" << endreq;
338 }
339 if(m_thsvc->getHist("/DQAHist/MUC/Phi_All", htmp).isSuccess()) {
340 htmp->Fill( phi );
341 } else {
342 log << MSG::ERROR << "Fail to retrieve Phi_All" << endreq;
343 }
344
345 if( isDimu )
346 {
347 if(m_thsvc->getHist("/DQAHist/MUC/Costheta_Dimu", htmp).isSuccess()) {
348 htmp->Fill( costheta );
349 } else {
350 log << MSG::ERROR << "Fail to retrieve Costheta_Dimu" << endreq;
351 }
352 if(m_thsvc->getHist("/DQAHist/MUC/Phi_Dimu", htmp).isSuccess()) {
353 htmp->Fill( phi );
354 } else {
355 log << MSG::ERROR << "Fail to retrieve Phi_Dimu" << endreq;
356 }
357 }
358 log << MSG::INFO << "Fill costheta and phi:\t" << costheta << "\t" << phi << endreq;
359
360 MucRecHit* pMucRawHit;
361 MucRecHit* pMucExpHit;
362
363
364 // Expected hits in this rec track
365 mucExpHitCol = (*trackIter)->GetExpectedHits();
366 //mucExpHitCol = (*trackIter)->getExpHits();
367 expectedHitNum += mucExpHitCol.size();
368 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
369 {
370 pMucRawHit = mucExpHitCol[ hitId ];
371 part = pMucRawHit->Part(); segment = pMucRawHit->Seg();
372 layer = pMucRawHit->Gap(); strip = pMucRawHit->Strip();
373
374 MucMark* currentMark = new MucMark( part, segment, layer, strip );
375 m_expHitCol.push_back( currentMark );
376
377 // Judge efficiency hit
378 int isInPos = -1;
379 bool isInEffWindow = false;
380 isInPos = currentMark->IsInCol( m_segDigiCol[part][segment] );
381
382 // Avoid bias in outer layers caused by low momentum tracks
383 if( part == BRID && (layer-lastLayerBR>1) ) continue;
384 if( part != BRID && (layer-lastLayerEC>1) ) continue;
385
386 // Avoid bias in both sides of the innermost layer of Barrel
387 if( part==BRID && layer==0 && (strip<2 || strip>45) )
388 {
389 if( isInPos != -1) // expHit is fired
390 {
391 m_recordAll[part][segment][layer][strip][2] ++; // Efficiency hit number
392 m_recordAll[part][segment][layer][strip][1] ++; // Rec track number
393 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
394
395 if(isDimu) {
396 m_recordDimu[part][segment][layer][strip][2] ++; // Efficiency hit number
397 m_recordDimu[part][segment][layer][strip][1] ++; // Rec track number
398 }
399 }
400 else {
401 m_recordAll[part][segment][layer][strip][1] ++;
402 if(isDimu) m_recordDimu[part][segment][layer][strip][1] ++;
403 }
404 continue; // Judge next hit
405 }
406 // Eff calibration
407 if( isInPos != -1 ) // expHit is fired
408 {
409 m_recordAll[part][segment][layer][strip][2] ++; // Efficiency hit number
410 m_recordAll[part][segment][layer][strip][1] ++; // Rec track number
411 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
412
413 if(isDimu) {
414 m_recordDimu[part][segment][layer][strip][2] ++; // Efficiency hit number
415 m_recordDimu[part][segment][layer][strip][1] ++; // Rec track number
416 }
417
418 continue; // Judge next hit
419 }
420 else for(int tempStrip=0, hiti=-m_effWindow; hiti<=m_effWindow; hiti++ )
421 {
422 if( hiti == 0 ) continue;
423 tempStrip = strip + hiti;
424 if( tempStrip < 0 || tempStrip > m_ptrIdTr->GetStripMax(part,segment,layer) ) continue;
425
426 isInPos = m_ptrMucMark->IsInCol( part, segment, layer, tempStrip, m_segDigiCol[part][segment] );
427 if( isInPos != -1 )
428 {
429 m_recordAll[part][segment][layer][tempStrip][2] ++; // Efficiency hit number
430 m_recordAll[part][segment][layer][tempStrip][1] ++; // Rec track number
431 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
432
433 if(isDimu) {
434 m_recordDimu[part][segment][layer][tempStrip][2] ++; // Efficiency hit number
435 m_recordDimu[part][segment][layer][tempStrip][1] ++; // Rec track number
436 }
437
438 isInEffWindow = true;
439 }
440 } // End else
441
442 if( isInEffWindow ) { continue; } // Judge next hit
443 else { // A hit should be fired but not fired and not in the EffWindow
444 m_recordAll[part][segment][layer][strip][1] ++; // Rec track number
445 if(isDimu) m_recordDimu[part][segment][layer][strip][1] ++;
446 }
447
448 } // End expected hits
449
450 // Fill residual, and for the other way of eff calculation
451 log << MSG::INFO << "Fill residual" << endreq;
452 if(!m_dstFileOnly)
453 {
454 vector<float> m_lineResCol = (*trackIter)->getDistHits();
455 // Digis belong to this rec track
456 mucRawHitCol = (*trackIter)->GetHits(); // Get hit collection of a track
457 //mucRawHitCol = (*trackIter)->getVecHits(); // Get hit collection of a track
458
459 if( trackHitNum > 4 && m_lineResCol[0] != -99) // track is good for res
460 {
461 // Fill res histograms
462 bool firedFlag[PART_MAX][LAYER_MAX][2];
463 for(int iprt=0; iprt<PART_MAX; iprt++)
464 for(int jlay=0; jlay<LAYER_MAX; jlay++)
465 firedFlag[iprt][jlay][0] = firedFlag[iprt][jlay][1] = false;
466
467 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
468 {
469 pMucExpHit = mucExpHitCol[ hitId ];
470 part = pMucExpHit->Part(); segment = pMucExpHit->Seg(); layer = pMucExpHit->Gap();
471 firedFlag[part][layer][0] = true;
472 }
473
474 for(unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
475 {
476 pMucRawHit = mucRawHitCol[ hitId ];
477 part = pMucRawHit->Part(); segment = pMucRawHit->Seg(); layer = pMucRawHit->Gap();
478
479 // if exp is true and fired is true, and not filled yet
480 if( firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false )
481 {
482 if( part == BRID )
483 {
484 sprintf( name, "/DQAHist/MUC/BrResDist_All_L%d", layer );
485 if(m_thsvc->getHist(name, htmp).isSuccess()) {
486 htmp->Fill( m_lineResCol[hitId] );
487 } else {
488 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
489 }
490 //m_hBrResDist[layer][0]->Fill( m_lineResCol[hitId] );
491
492 if( isDimu )
493 {
494 sprintf( name, "/DQAHist/MUC/BrResDist_Dimu_L%d", layer );
495 if(m_thsvc->getHist(name, htmp).isSuccess()) {
496 htmp->Fill( m_lineResCol[hitId] );
497 } else {
498 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
499 }
500 //m_hBrResDist[layer][1]->Fill( m_lineResCol[hitId] );
501 }
502 } else {
503 sprintf( name, "/DQAHist/MUC/EcResDist_All_L%d", layer );
504 if(m_thsvc->getHist(name, htmp).isSuccess()) {
505 htmp->Fill( m_lineResCol[hitId] );
506 } else {
507 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
508 }
509 //m_hEcResDist[layer][0]->Fill( m_lineResCol[hitId] );
510
511 if( isDimu )
512 {
513 sprintf( name, "/DQAHist/MUC/EcResDist_Dimu_L%d", layer );
514 if(m_thsvc->getHist(name, htmp).isSuccess()) {
515 htmp->Fill( m_lineResCol[hitId] );
516 } else {
517 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
518 }
519 //m_hEcResDist[layer][1]->Fill( m_lineResCol[hitId] );
520 }
521 }
522 }
523
524 firedFlag[part][layer][1] = true;
525 }
526 } // End fill residual, if track is good for res
527 }
528 mucRawHitCol.clear();
529 mucExpHitCol.clear();
530
531 } // End read all tracks
532 //<--- End retrieve rec tracks
533
534 //---> Searching inc/noise hits
535 log << MSG::INFO << "Searching inc/noise hits" << endreq;
536 bool isNosHit;
537 bool hasEffHit;
538 for( unsigned int i=0; i < m_digiCol.size(); i++ )
539 {
540 isNosHit = true;
541
542 if( m_digiCol[i]->IsInCol( m_effHitCol )!=-1 ) continue; // digi in effHitCol
543 else
544 {
545 for( unsigned int j=0; j < m_clusterCol.size(); j++ )
546 {
547 hasEffHit = false;
548 for( unsigned int k=0; k<m_clusterCol[j].size(); k++)
549 {
550 if( m_clusterCol[j][k]->IsInCol(m_effHitCol) != -1)
551 { // Clusters have efficiency hit
552 hasEffHit = true;
553 break; // Out a cluster
554 }
555 }
556
557 if ( hasEffHit && (m_digiCol[i]->IsInCol( m_clusterCol[j] ) != -1) ) {
558 isNosHit = false;
559 break; // Out cluster collection
560 }
561 } // End cluster col
562
563 if( isNosHit ) {
564 part = (*m_digiCol[i]).Part(); segment = (*m_digiCol[i]).Segment();
565 layer = (*m_digiCol[i]).Layer(); strip = (*m_digiCol[i]).Strip();
566
567 //log << MSG::INFO << "Noise hit:\t" << part<<"\t"<<segment<<"\t"<<layer<<"\t"<<strip << endreq;
568 m_recordAll[part][segment][layer][strip][3] ++;
569 if( isDimu ) m_recordDimu[part][segment][layer][strip][3] ++;
570 m_nosHitCol.push_back(m_digiCol[i]);
571 }
572 }// End else
573 } // End digi collection
574 //<--- End searching inc/noise hits
575
576 DQA_MUC::FillHistograms( isDimu );
577
578 // Reset mark collection
579 for( unsigned int i=0; i<m_digiCol.size(); i++ ) {
580 if( m_digiCol[i] != NULL ) delete m_digiCol[i];
581 }
582
583 for( unsigned int i=0; i<m_expHitCol.size(); i++ ) {
584 if( m_expHitCol[i] != NULL ) delete m_expHitCol[i];
585 }
586
587 if( m_effHitCol.size() != 0 ) m_effHitCol.clear();
588 if( m_expHitCol.size() != 0 ) m_expHitCol.clear();
589 if( m_nosHitCol.size() != 0 ) m_nosHitCol.clear();
590 if( m_digiCol.size() != 0 ) m_digiCol.clear();
591 if( m_clusterCol.size() != 0 ) m_clusterCol.clear();
592
593 for( int i=0; i<PART_MAX; i++ ) {
594 for( int j=0; j<SEGMENT_MAX; j++ ) {
595 if( m_segDigiCol[i][j].size() != 0 ) m_segDigiCol[i][j].clear();
596 }
597 }
598
599 return StatusCode::SUCCESS;
600}
601
602// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
603StatusCode DQA_MUC::finalize() {
604
605 MsgStream log(msgSvc(), name());
606 log << MSG::INFO << "in finalize()" << endmsg;
607
608/*
609 TH1 *htmp(0);
610
611 // Fit spacial resolution
612 log << MSG::INFO << "Spacial resolution fitting" << endreq;
613 double resSigma, resSigmaErr;
614 resSigma = resSigmaErr = 0.;
615
616 for( int i=0; i<B_LAY_NUM; i++ )
617 {
618 m_hBrResDist[i][0]->Fit("gaus");
619 resSigma = m_hBrResDist[i][0]->GetFunction("gaus")->GetParameter("Sigma");
620 resSigmaErr = m_hBrResDist[i][0]->GetFunction("gaus")->GetParError(2);
621 if(m_thsvc->getHist("/DQAHist/MUC/BrRes_All", htmp).isSuccess()) {
622 htmp->Fill( i, resSigma );
623 htmp->SetBinError( i+1, resSigmaErr );
624 } else {
625 log << MSG::ERROR << "Fail to retrieve BrRes_All" << endreq;
626 }
627 //log << MSG::INFO << "Barrle layer:\t" << i << "\t" << resSigma << "\t" << resSigmaErr << endreq;
628
629 m_hBrResDist[i][1]->Fit("gaus");
630 resSigma = m_hBrResDist[i][1]->GetFunction("gaus")->GetParameter("Sigma");
631 resSigmaErr = m_hBrResDist[i][1]->GetFunction("gaus")->GetParError(2);
632 if(m_thsvc->getHist("/DQAHist/MUC/BrRes_Dimu", htmp).isSuccess()) {
633 htmp->Fill( i, resSigma );
634 htmp->SetBinError( i+1, resSigmaErr );
635 } else {
636 log << MSG::ERROR << "Fail to retrieve BrRes_Dimu" << endreq;
637 }
638
639 }
640
641 for( int i=0; i<E_LAY_NUM; i++ )
642 {
643 m_hEcResDist[i][0]->Fit("gaus");
644 resSigma = m_hEcResDist[i][0]->GetFunction("gaus")->GetParameter("Sigma");
645 resSigmaErr = m_hEcResDist[i][0]->GetFunction("gaus")->GetParError(2);
646 if(m_thsvc->getHist("/DQAHist/MUC/EcRes_All", htmp).isSuccess()) {
647 htmp->Fill( i, resSigma );
648 htmp->SetBinError( i+1, resSigmaErr );
649 } else {
650 log << MSG::ERROR << "Fail to retrieve EcRes_All" << endreq;
651 }
652 //log << MSG::INFO << "Endcap layer:\t" << i << "\t" << resSigma << "\t" << resSigmaErr << endreq;
653
654 m_hEcResDist[i][1]->Fit("gaus");
655 resSigma = m_hEcResDist[i][1]->GetFunction("gaus")->GetParameter("Sigma");
656 resSigmaErr = m_hEcResDist[i][0]->GetFunction("gaus")->GetParError(2);
657 if(m_thsvc->getHist("/DQAHist/MUC/EcRes_Dimu", htmp).isSuccess()) {
658 htmp->Fill( i, resSigma );
659 htmp->SetBinError( i+1, resSigmaErr );
660 } else {
661 log << MSG::ERROR << "Fail to retrieve EcRes_Dimu" << endreq;
662 }
663 }
664
665
666 // Efficiency and noise ratio
667 int part, segment, layer, stripMax;
668 part = segment = layer = stripMax = 0;
669 double digi[TAGN][3], effHit[TAGN][3], trackNum[TAGN][3], nosHit[TAGN][3];
670 double eff[TAGN], effErr[TAGN], nos[TAGN], nosErr[TAGN];
671 for(int i=0; i<TAGN; i++) {
672 for( int j=0; j<3; j++) {
673 digi[i][j] = effHit[i][j] = trackNum[i][j] = nosHit[i][j] = 0.;
674 }
675
676 eff[i] = effErr[i] = nos[i] = nosErr[i] = 0.;
677 }
678
679 // Layer efficiency
680 for( int i=0; i<LAYER_MAX; i++ )
681 {
682 effHit[0][0] = trackNum[0][0] = 0;
683 effHit[1][0] = trackNum[1][0] = 0;
684
685 for( int j=0; j<PART_MAX; j++ ) {
686 for( int k=0; k<SEGMENT_MAX; k++) {
687 stripMax = m_ptrIdTr->GetStripMax( j, k, i );
688 for( int n=0; n<stripMax; n++ ) {
689 trackNum[0][0] += m_recordAll[j][k][i][n][1];
690 effHit[0][0] += m_recordAll[j][k][i][n][2];
691
692 trackNum[1][0] += m_recordDimu[j][k][i][n][1];
693 effHit[1][0] += m_recordDimu[j][k][i][n][2];
694 }
695 }
696 }
697
698 if( trackNum[0][0] == 0 ) {
699 eff[0] = effErr[0] = 0.0;
700 } else {
701 eff[0] = ( (double)effHit[0][0] )/trackNum[0][0];
702 effErr[0] = sqrt( eff[0]*(1-eff[0])/trackNum[0][0] );
703 }
704 if(m_thsvc->getHist("/DQAHist/MUC/LayerEff_All", htmp).isSuccess()) {
705 htmp->Fill( i, eff[0] );
706 htmp->SetBinError( i+1, effErr[0] );
707 } else {
708 log << MSG::ERROR << "Fail to retrieve LayerEff_All" << endreq;
709 }
710 log << MSG::INFO << "LayerEff:\t" << i << "\t" << eff[0] << "\t" << effErr[0] << endreq;
711
712 if( trackNum[1][0] == 0 ) {
713 eff[1] = effErr[1] = 0.0;
714 } else {
715 eff[1] = ( (double)effHit[1][0] )/trackNum[1][0];
716 effErr[1] = sqrt( eff[1]*(1-eff[1])/trackNum[1][0] );
717 }
718 if(m_thsvc->getHist("/DQAHist/MUC/LayerEff_Dimu", htmp).isSuccess()) {
719 htmp->Fill( i, eff[1] );
720 htmp->SetBinError( i+1, effErr[1] );
721 } else {
722 log << MSG::ERROR << "Fail to retrieve LayerEff_Dimu" << endreq;
723 }
724
725 } // End layer
726
727 // Box efficiency and noise ratio, strip efficiency
728 for( int i=0; i<BOX_MAX; i++ )
729 {
730 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
731 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
732
733 digi[0][1] = effHit[0][1] = trackNum[0][1] = nosHit[0][1] = 0;
734 digi[1][1] = effHit[1][1] = trackNum[1][1] = nosHit[1][1] = 0;
735 for( int j=0; j<stripMax; j++ )
736 {
737 // For box
738 digi[0][1] += m_recordAll[part][segment][layer][j][0];
739 trackNum[0][1] += m_recordAll[part][segment][layer][j][1];
740 effHit[0][1] += m_recordAll[part][segment][layer][j][2];
741 nosHit[0][1] += m_recordAll[part][segment][layer][j][3];
742
743 digi[1][1] += m_recordDimu[part][segment][layer][j][0];
744 trackNum[1][1] += m_recordDimu[part][segment][layer][j][1];
745 effHit[1][1] += m_recordDimu[part][segment][layer][j][2];
746 nosHit[1][1] += m_recordDimu[part][segment][layer][j][3];
747
748 // For strip
749 trackNum[0][2] = m_recordAll[part][segment][layer][j][1];
750 effHit[0][2] = m_recordAll[part][segment][layer][j][2];
751 trackNum[1][2] = m_recordDimu[part][segment][layer][j][1];
752 effHit[1][2] = m_recordDimu[part][segment][layer][j][2];
753
754 if( trackNum[0][2] == 0 ) {
755 eff[0] = 0.0;
756 } else {
757 eff[0] = ( (double)effHit[0][2] )/trackNum[0][2];
758 }
759 if(m_thsvc->getHist("/DQAHist/MUC/StripEff_All", htmp).isSuccess()) {
760 htmp->Fill( eff[0] );
761 } else {
762 log << MSG::ERROR << "Fail to retrieve StripEff_All" << endreq;
763 }
764
765 if( trackNum[1][2] == 0 ) {
766 eff[1] = 0.0;
767 } else {
768 eff[1] = ( (double)effHit[1][2] )/trackNum[1][2];
769 }
770 if(m_thsvc->getHist("/DQAHist/MUC/StripEff_Dimu", htmp).isSuccess()) {
771 htmp->Fill( eff[1] );
772 } else {
773 log << MSG::ERROR << "Fail to retrieve StripEff_Dimu" << endreq;
774 }
775 } // End StripMax
776
777 // Box efficiency
778 if( trackNum[0][1] == 0 ) {
779 eff[0] = effErr[0] = 0.0;
780 } else {
781 eff[0] = ( (double)effHit[0][1] )/trackNum[0][1];
782 effErr[0] = sqrt( eff[0]*(1-eff[0])/trackNum[0][1] );
783 }
784 if(m_thsvc->getHist("/DQAHist/MUC/BoxEff_All", htmp).isSuccess()) {
785 htmp->Fill( i, eff[0] );
786 htmp->SetBinError( i+1, effErr[0] );
787 } else {
788 log << MSG::ERROR << "Fail to retrieve BoxEff_All" << endreq;
789 }
790 log << MSG::INFO << "BoxEff:\t" << i << "\t" << eff[0] << "\t" << effErr[0] << endreq;
791
792 if( trackNum[1][1] == 0 ) {
793 eff[1] = effErr[1] = 0.0;
794 } else {
795 eff[1] = ( (double)effHit[1][1] )/trackNum[1][1];
796 effErr[1] = sqrt( eff[1]*(1-eff[1])/trackNum[1][1] );
797 }
798 if(m_thsvc->getHist("/DQAHist/MUC/BoxEff_Dimu", htmp).isSuccess()) {
799 htmp->Fill( i, eff[1] );
800 htmp->SetBinError( i+1, effErr[1] );
801 } else {
802 log << MSG::ERROR << "Fail to retrieve BoxEff_Dimu" << endreq;
803 }
804
805 // Box noise ratio
806 if( digi[0][1] == 0 ) {
807 nos[0] = nosErr[0] = 0.0;
808 } else {
809 nos[0] = ( (double)nosHit[0][1] )/digi[0][1];
810 nosErr[0] = sqrt( nos[0]*(1-nos[0])/digi[0][1] );
811 }
812 if(m_thsvc->getHist("/DQAHist/MUC/BoxNosRatio_All", htmp).isSuccess()) {
813 htmp->Fill( i, nos[0] );
814 htmp->SetBinError( i+1, nosErr[0] );
815 } else {
816 log << MSG::ERROR << "Fail to retrieve BoxNosRatio_All" << endreq;
817 }
818 log << MSG::INFO << "BoxNosRatio:\t" << i << "\t" << nos[0] << "\t" << nosErr[0] << endreq;
819
820 if( digi[1][1] == 0 ) {
821 nos[1] = nosErr[1] = 0.0;
822 } else {
823 nos[1] = ( (double)nosHit[1][1] )/digi[1][1];
824 nosErr[1] = sqrt( nos[1]*(1-nos[1])/digi[1][1] );
825 }
826 if(m_thsvc->getHist("/DQAHist/MUC/BoxNosRatio_Dimu", htmp).isSuccess()) {
827 htmp->Fill( i, nos[1] );
828 htmp->SetBinError( i+1, nosErr[1] );
829 } else {
830 log << MSG::ERROR << "Fail to retrieve BoxNosRatio_Dimu" << endreq;
831 }
832 }// End BOX_MAX
833*/
834
835 log << MSG::INFO << "Finalize successfully! " << endmsg;
836
837 return StatusCode::SUCCESS;
838}
839
840StatusCode DQA_MUC::FillHistograms(bool isDimu) {
841
842 MsgStream log(msgSvc(), name());
843 log << MSG::INFO << "Filling histograms" << endmsg;
844 char name[100];
845 int part, segment, layer, strip, boxId, strId;
846 part = segment = layer = strip = boxId = strId = 0;
847 TH1 *hEff[2][2][3];
848
849 for(int i=0; i<TAGN; i++) {
850 for(int j=0; j<2; j++) {
851 for(int k=0; k<LVLN; k++)
852 {
853 sprintf( name, "/DQAHist/MUC/%sEff_%s_%s", LNAME[k], ENAME[i], HTYPE[j] );
854 if(!m_thsvc->getHist(name, hEff[i][j][k]).isSuccess())
855 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
856 }
857 }
858 }
859
860 // Numberator of eff
861 for(int i=0; i<m_effHitCol.size(); i++)
862 {
863 part = m_effHitCol[i]->Part();
864 segment = m_effHitCol[i]->Segment();
865 layer = m_effHitCol[i]->Layer();
866 strip = m_effHitCol[i]->Strip();
867 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
868 strId = m_ptrIdTr->GetStripId(part, segment, layer, strip);
869 hEff[0][0][0]->Fill(layer);
870 hEff[0][0][1]->Fill(boxId);
871 hEff[0][0][2]->Fill(strId);
872
873 if( isDimu ) {
874 hEff[1][0][0]->Fill(layer);
875 hEff[1][0][1]->Fill(boxId);
876 hEff[1][0][2]->Fill(strId);
877 }
878 }
879
880 // Denominator of eff
881 for(int i=0; i<m_expHitCol.size(); i++)
882 {
883 part = m_expHitCol[i]->Part();
884 segment = m_expHitCol[i]->Segment();
885 layer = m_expHitCol[i]->Layer();
886 strip = m_expHitCol[i]->Strip();
887 //cout<<"ExpHit:\t"<<part<<"\t"<<segment<<"\t"<<layer<<"\t"<<strip<<endl;
888 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
889 strId = m_ptrIdTr->GetStripId(part, segment, layer, strip);
890 hEff[0][1][0]->Fill(layer);
891 hEff[0][1][1]->Fill(boxId);
892 hEff[0][1][2]->Fill(strId);
893
894 if( isDimu ) {
895 hEff[1][1][0]->Fill(layer);
896 hEff[1][1][1]->Fill(boxId);
897 hEff[1][1][2]->Fill(strId);
898 }
899
900 }
901
902 TH1 *hNosRatio[2][2];
903 for(int i=0; i<TAGN; i++) {
904 for(int j=0; j<2; j++) {
905 sprintf( name, "/DQAHist/MUC/%sNosRatio_%s_%s", LNAME[1], ENAME[i], HTYPE[j] );
906 if(!m_thsvc->getHist(name, hNosRatio[i][j]).isSuccess())
907 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
908 }
909 }
910
911 // Numberator of box noise ratio
912 for(int i=0; i<m_nosHitCol.size(); i++)
913 {
914 part = m_nosHitCol[i]->Part();
915 segment = m_nosHitCol[i]->Segment();
916 layer = m_nosHitCol[i]->Layer();
917 strip = m_nosHitCol[i]->Strip();
918 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
919 strId = m_ptrIdTr->GetStripId(part, segment, layer, strip);
920 hNosRatio[0][0]->Fill(boxId);
921
922 if( isDimu ) hNosRatio[1][0]->Fill(boxId);
923 }
924
925 // Denominator of box noise ratio
926 for(int i=0; i<m_digiCol.size(); i++)
927 {
928 part = m_digiCol[i]->Part();
929 segment = m_digiCol[i]->Segment();
930 layer = m_digiCol[i]->Layer();
931 strip = m_digiCol[i]->Strip();
932 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
933 hNosRatio[0][1]->Fill(boxId);
934
935 if( isDimu ) hNosRatio[1][1]->Fill(boxId);
936 }
937
938 return StatusCode::SUCCESS;
939}
940
941//END
942
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
const char ENAME[TAGN][10]
Definition: DQA_MUC.h:18
const char HTYPE[2][10]
Definition: DQA_MUC.h:19
const char LNAME[LVLN][10]
Definition: DQA_MUC.h:20
const int TAGN
Definition: DQA_MUC.h:16
const int LVLN
Definition: DQA_MUC.h:17
const Int_t n
titledef title[20]
const int STRIP_MAX
Definition: MucMappingAlg.h:24
const double PI
Definition: PipiJpsi.cxx:55
ObjectVector< RecMucTrack > RecMucTrackCol
Definition: RecMucTrack.h:394
IMessageSvc * msgSvc()
#define Segment
Definition: TTrackBase.h:31
#define NULL
StatusCode execute()
Definition: DQA_MUC.cxx:228
StatusCode initialize()
Definition: DQA_MUC.cxx:45
StatusCode finalize()
Definition: DQA_MUC.cxx:603
static int barrel_ec(const Identifier &id)
Values of different levels.
Definition: MucID.cxx:41
static int layer(const Identifier &id)
Definition: MucID.cxx:61
static int channel(const Identifier &id)
Definition: MucID.cxx:71
static int segment(const Identifier &id)
Definition: MucID.cxx:51
int GetStripId(int part, int segment, int layer, int subid)
int GetStripMax(int part, int segment, int layer)
int GetBoxId(int part, int segment, int layer)
int IsInCol(int part, int segment, int layer, int strip, mark_col &aCol)
Definition: MucMark.cxx:100
int Seg() const
Get Seg.
Definition: MucRecHit.h:74
int Part() const
Get Part.
Definition: MucRecHit.h:71
int Gap() const
Get Gap.
Definition: MucRecHit.h:77
int Strip() const
Get Strip.
Definition: MucRecHit.h:80
_EXTERN_ std::string EvtRecEvent
Definition: EventModel.h:116
float costheta