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