BOSS 7.0.4
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 log << MSG::WARNING << "Could not find RecMucTrackCol" << endreq; // until BOSS 7.0.3, if not MUC track, no RecMucTrackCol either.
291 return( StatusCode::SUCCESS);
292 }
293
294 RecMucTrackCol *aRecMucTrackCol = mucTrackCol;
295 if (aRecMucTrackCol->size() < 1) {
296 log << MSG::INFO << "No MUC tracks in this event" << endreq;
297 return StatusCode::SUCCESS;
298 }
299 log << MSG::INFO << "Total tracks of this event: " << aRecMucTrackCol->size() << endreq;
300
301 int trackHitNum, expectedHitNum, segNum, lastLayerBR, lastLayerEC;
302 int layerPassNum[3], passMax[TRACK_SEG_MAX][2];
303 bool firedLay[TRACK_SEG_MAX][LAYER_MAX];
304 double costheta, phi;
305 TH1* htmp(0);
306 trackHitNum = expectedHitNum = segNum = lastLayerBR = lastLayerEC = 0;
307 layerPassNum[0] = layerPassNum[1] = layerPassNum[2] = 0;
308 for( int segi=0; segi<TRACK_SEG_MAX; segi++ ) {
309 passMax[segi][0] = passMax[segi][1] = 0;
310 for( int layi=0; layi<LAYER_MAX; layi++ ) firedLay[segi][layi] = 0;
311 }
312
313 vector<MucRecHit*> mucRawHitCol;
314 vector<MucRecHit*> mucExpHitCol;
315 RecMucTrackCol::iterator trackIter = mucTrackCol->begin();
316 for (int trackId = 0; trackIter != mucTrackCol->end(); trackIter++, trackId++)
317 {
318 trackHitNum = (*trackIter)->numHits();
319
320 if( trackHitNum == 0 ) {
321 log << MSG::INFO << "Track " << trackId << " no hits" << endreq;
322 continue;
323 }
324
325 lastLayerBR = (*trackIter)->brLastLayer();
326 lastLayerEC = (*trackIter)->ecLastLayer();
327 // First fit position in MUC
328 CLHEP::Hep3Vector a3Vector((*trackIter)->xPos(),(*trackIter)->yPos(),(*trackIter)->zPos());
329 costheta = a3Vector.cosTheta();
330 phi = a3Vector.phi();
331 if(m_thsvc->getHist("/DQAHist/MUC/Costheta_All", htmp).isSuccess()) {
332 htmp->Fill( costheta );
333 } else {
334 log << MSG::ERROR << "Fail to retrieve Costheta_All" << endreq;
335 }
336 if(m_thsvc->getHist("/DQAHist/MUC/Phi_All", htmp).isSuccess()) {
337 htmp->Fill( phi );
338 } else {
339 log << MSG::ERROR << "Fail to retrieve Phi_All" << endreq;
340 }
341
342 if( isDimu )
343 {
344 if(m_thsvc->getHist("/DQAHist/MUC/Costheta_Dimu", htmp).isSuccess()) {
345 htmp->Fill( costheta );
346 } else {
347 log << MSG::ERROR << "Fail to retrieve Costheta_Dimu" << endreq;
348 }
349 if(m_thsvc->getHist("/DQAHist/MUC/Phi_Dimu", htmp).isSuccess()) {
350 htmp->Fill( phi );
351 } else {
352 log << MSG::ERROR << "Fail to retrieve Phi_Dimu" << endreq;
353 }
354 }
355 log << MSG::INFO << "Fill costheta and phi:\t" << costheta << "\t" << phi << endreq;
356
357 MucRecHit* pMucRawHit;
358 MucRecHit* pMucExpHit;
359
360
361 // Expected hits in this rec track
362 mucExpHitCol = (*trackIter)->GetExpectedHits();
363 //mucExpHitCol = (*trackIter)->getExpHits();
364 expectedHitNum += mucExpHitCol.size();
365 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
366 {
367 pMucRawHit = mucExpHitCol[ hitId ];
368 part = pMucRawHit->Part(); segment = pMucRawHit->Seg();
369 layer = pMucRawHit->Gap(); strip = pMucRawHit->Strip();
370
371 MucMark* currentMark = new MucMark( part, segment, layer, strip );
372 m_expHitCol.push_back( currentMark );
373
374 // Judge efficiency hit
375 int isInPos = -1;
376 bool isInEffWindow = false;
377 isInPos = currentMark->IsInCol( m_segDigiCol[part][segment] );
378
379 // Avoid bias in outer layers caused by low momentum tracks
380 if( part == BRID && (layer-lastLayerBR>1) ) continue;
381 if( part != BRID && (layer-lastLayerEC>1) ) continue;
382
383 // Avoid bias in both sides of the innermost layer of Barrel
384 if( part==BRID && layer==0 && (strip<2 || strip>45) )
385 {
386 if( isInPos != -1) // expHit is fired
387 {
388 m_recordAll[part][segment][layer][strip][2] ++; // Efficiency hit number
389 m_recordAll[part][segment][layer][strip][1] ++; // Rec track number
390 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
391
392 if(isDimu) {
393 m_recordDimu[part][segment][layer][strip][2] ++; // Efficiency hit number
394 m_recordDimu[part][segment][layer][strip][1] ++; // Rec track number
395 }
396 }
397 else {
398 m_recordAll[part][segment][layer][strip][1] ++;
399 if(isDimu) m_recordDimu[part][segment][layer][strip][1] ++;
400 }
401 continue; // Judge next hit
402 }
403 // Eff calibration
404 if( isInPos != -1 ) // expHit is fired
405 {
406 m_recordAll[part][segment][layer][strip][2] ++; // Efficiency hit number
407 m_recordAll[part][segment][layer][strip][1] ++; // Rec track number
408 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
409
410 if(isDimu) {
411 m_recordDimu[part][segment][layer][strip][2] ++; // Efficiency hit number
412 m_recordDimu[part][segment][layer][strip][1] ++; // Rec track number
413 }
414
415 continue; // Judge next hit
416 }
417 else for(int tempStrip=0, hiti=-m_effWindow; hiti<=m_effWindow; hiti++ )
418 {
419 if( hiti == 0 ) continue;
420 tempStrip = strip + hiti;
421 if( tempStrip < 0 || tempStrip > m_ptrIdTr->GetStripMax(part,segment,layer) ) continue;
422
423 isInPos = m_ptrMucMark->IsInCol( part, segment, layer, tempStrip, m_segDigiCol[part][segment] );
424 if( isInPos != -1 )
425 {
426 m_recordAll[part][segment][layer][tempStrip][2] ++; // Efficiency hit number
427 m_recordAll[part][segment][layer][tempStrip][1] ++; // Rec track number
428 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
429
430 if(isDimu) {
431 m_recordDimu[part][segment][layer][tempStrip][2] ++; // Efficiency hit number
432 m_recordDimu[part][segment][layer][tempStrip][1] ++; // Rec track number
433 }
434
435 isInEffWindow = true;
436 }
437 } // End else
438
439 if( isInEffWindow ) { continue; } // Judge next hit
440 else { // A hit should be fired but not fired and not in the EffWindow
441 m_recordAll[part][segment][layer][strip][1] ++; // Rec track number
442 if(isDimu) m_recordDimu[part][segment][layer][strip][1] ++;
443 }
444
445 } // End expected hits
446
447 // Fill residual, and for the other way of eff calculation
448 log << MSG::INFO << "Fill residual" << endreq;
449 vector<float> m_lineResCol = (*trackIter)->getDistHits();
450 // Digis belong to this rec track
451 mucRawHitCol = (*trackIter)->GetHits(); // Get hit collection of a track
452 //mucRawHitCol = (*trackIter)->getVecHits(); // Get hit collection of a track
453
454 if( trackHitNum > 4 && m_lineResCol[0] != -99) // track is good for res
455 {
456 // Fill res histograms
457 bool firedFlag[PART_MAX][LAYER_MAX][2];
458 for(int iprt=0; iprt<PART_MAX; iprt++)
459 for(int jlay=0; jlay<LAYER_MAX; jlay++)
460 firedFlag[iprt][jlay][0] = firedFlag[iprt][jlay][1] = false;
461
462 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
463 {
464 pMucExpHit = mucExpHitCol[ hitId ];
465 part = pMucExpHit->Part(); segment = pMucExpHit->Seg(); layer = pMucExpHit->Gap();
466 firedFlag[part][layer][0] = true;
467 }
468
469 for(unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
470 {
471 pMucRawHit = mucRawHitCol[ hitId ];
472 part = pMucRawHit->Part(); segment = pMucRawHit->Seg(); layer = pMucRawHit->Gap();
473
474 // if exp is true and fired is true, and not filled yet
475 if( firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false )
476 {
477 if( part == BRID )
478 {
479 sprintf( name, "/DQAHist/MUC/BrResDist_All_L%d", layer );
480 if(m_thsvc->getHist(name, htmp).isSuccess()) {
481 htmp->Fill( m_lineResCol[hitId] );
482 } else {
483 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
484 }
485 //m_hBrResDist[layer][0]->Fill( m_lineResCol[hitId] );
486
487 if( isDimu )
488 {
489 sprintf( name, "/DQAHist/MUC/BrResDist_Dimu_L%d", layer );
490 if(m_thsvc->getHist(name, htmp).isSuccess()) {
491 htmp->Fill( m_lineResCol[hitId] );
492 } else {
493 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
494 }
495 //m_hBrResDist[layer][1]->Fill( m_lineResCol[hitId] );
496 }
497 } else {
498 sprintf( name, "/DQAHist/MUC/EcResDist_All_L%d", layer );
499 if(m_thsvc->getHist(name, htmp).isSuccess()) {
500 htmp->Fill( m_lineResCol[hitId] );
501 } else {
502 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
503 }
504 //m_hEcResDist[layer][0]->Fill( m_lineResCol[hitId] );
505
506 if( isDimu )
507 {
508 sprintf( name, "/DQAHist/MUC/EcResDist_Dimu_L%d", layer );
509 if(m_thsvc->getHist(name, htmp).isSuccess()) {
510 htmp->Fill( m_lineResCol[hitId] );
511 } else {
512 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
513 }
514 //m_hEcResDist[layer][1]->Fill( m_lineResCol[hitId] );
515 }
516 }
517 }
518
519 firedFlag[part][layer][1] = true;
520 }
521 } // End fill residual, if track is good for res
522
523 mucRawHitCol.clear();
524 mucExpHitCol.clear();
525
526 } // End read all tracks
527 //<--- End retrieve rec tracks
528
529 //---> Searching inc/noise hits
530 log << MSG::INFO << "Searching inc/noise hits" << endreq;
531 bool isNosHit;
532 bool hasEffHit;
533 for( unsigned int i=0; i < m_digiCol.size(); i++ )
534 {
535 isNosHit = true;
536
537 if( m_digiCol[i]->IsInCol( m_effHitCol )!=-1 ) continue; // digi in effHitCol
538 else
539 {
540 for( unsigned int j=0; j < m_clusterCol.size(); j++ )
541 {
542 hasEffHit = false;
543 for( unsigned int k=0; k<m_clusterCol[j].size(); k++)
544 {
545 if( m_clusterCol[j][k]->IsInCol(m_effHitCol) != -1)
546 { // Clusters have efficiency hit
547 hasEffHit = true;
548 break; // Out a cluster
549 }
550 }
551
552 if ( hasEffHit && (m_digiCol[i]->IsInCol( m_clusterCol[j] ) != -1) ) {
553 isNosHit = false;
554 break; // Out cluster collection
555 }
556 } // End cluster col
557
558 if( isNosHit ) {
559 part = (*m_digiCol[i]).Part(); segment = (*m_digiCol[i]).Segment();
560 layer = (*m_digiCol[i]).Layer(); strip = (*m_digiCol[i]).Strip();
561
562 //log << MSG::INFO << "Noise hit:\t" << part<<"\t"<<segment<<"\t"<<layer<<"\t"<<strip << endreq;
563 m_recordAll[part][segment][layer][strip][3] ++;
564 if( isDimu ) m_recordDimu[part][segment][layer][strip][3] ++;
565 m_nosHitCol.push_back(m_digiCol[i]);
566 }
567 }// End else
568 } // End digi collection
569 //<--- End searching inc/noise hits
570
571 DQA_MUC::FillHistograms( isDimu );
572
573 // Reset mark collection
574 for( unsigned int i=0; i<m_digiCol.size(); i++ ) {
575 if( m_digiCol[i] != NULL ) delete m_digiCol[i];
576 }
577
578 for( unsigned int i=0; i<m_expHitCol.size(); i++ ) {
579 if( m_expHitCol[i] != NULL ) delete m_expHitCol[i];
580 }
581
582 if( m_effHitCol.size() != 0 ) m_effHitCol.clear();
583 if( m_expHitCol.size() != 0 ) m_expHitCol.clear();
584 if( m_nosHitCol.size() != 0 ) m_nosHitCol.clear();
585 if( m_digiCol.size() != 0 ) m_digiCol.clear();
586 if( m_clusterCol.size() != 0 ) m_clusterCol.clear();
587
588 for( int i=0; i<PART_MAX; i++ ) {
589 for( int j=0; j<SEGMENT_MAX; j++ ) {
590 if( m_segDigiCol[i][j].size() != 0 ) m_segDigiCol[i][j].clear();
591 }
592 }
593
594 return StatusCode::SUCCESS;
595}
596
597// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
598StatusCode DQA_MUC::finalize() {
599
600 MsgStream log(msgSvc(), name());
601 log << MSG::INFO << "in finalize()" << endmsg;
602
603/*
604 TH1 *htmp(0);
605
606 // Fit spacial resolution
607 log << MSG::INFO << "Spacial resolution fitting" << endreq;
608 double resSigma, resSigmaErr;
609 resSigma = resSigmaErr = 0.;
610
611 for( int i=0; i<B_LAY_NUM; i++ )
612 {
613 m_hBrResDist[i][0]->Fit("gaus");
614 resSigma = m_hBrResDist[i][0]->GetFunction("gaus")->GetParameter("Sigma");
615 resSigmaErr = m_hBrResDist[i][0]->GetFunction("gaus")->GetParError(2);
616 if(m_thsvc->getHist("/DQAHist/MUC/BrRes_All", htmp).isSuccess()) {
617 htmp->Fill( i, resSigma );
618 htmp->SetBinError( i+1, resSigmaErr );
619 } else {
620 log << MSG::ERROR << "Fail to retrieve BrRes_All" << endreq;
621 }
622 //log << MSG::INFO << "Barrle layer:\t" << i << "\t" << resSigma << "\t" << resSigmaErr << endreq;
623
624 m_hBrResDist[i][1]->Fit("gaus");
625 resSigma = m_hBrResDist[i][1]->GetFunction("gaus")->GetParameter("Sigma");
626 resSigmaErr = m_hBrResDist[i][1]->GetFunction("gaus")->GetParError(2);
627 if(m_thsvc->getHist("/DQAHist/MUC/BrRes_Dimu", htmp).isSuccess()) {
628 htmp->Fill( i, resSigma );
629 htmp->SetBinError( i+1, resSigmaErr );
630 } else {
631 log << MSG::ERROR << "Fail to retrieve BrRes_Dimu" << endreq;
632 }
633
634 }
635
636 for( int i=0; i<E_LAY_NUM; i++ )
637 {
638 m_hEcResDist[i][0]->Fit("gaus");
639 resSigma = m_hEcResDist[i][0]->GetFunction("gaus")->GetParameter("Sigma");
640 resSigmaErr = m_hEcResDist[i][0]->GetFunction("gaus")->GetParError(2);
641 if(m_thsvc->getHist("/DQAHist/MUC/EcRes_All", htmp).isSuccess()) {
642 htmp->Fill( i, resSigma );
643 htmp->SetBinError( i+1, resSigmaErr );
644 } else {
645 log << MSG::ERROR << "Fail to retrieve EcRes_All" << endreq;
646 }
647 //log << MSG::INFO << "Endcap layer:\t" << i << "\t" << resSigma << "\t" << resSigmaErr << endreq;
648
649 m_hEcResDist[i][1]->Fit("gaus");
650 resSigma = m_hEcResDist[i][1]->GetFunction("gaus")->GetParameter("Sigma");
651 resSigmaErr = m_hEcResDist[i][0]->GetFunction("gaus")->GetParError(2);
652 if(m_thsvc->getHist("/DQAHist/MUC/EcRes_Dimu", htmp).isSuccess()) {
653 htmp->Fill( i, resSigma );
654 htmp->SetBinError( i+1, resSigmaErr );
655 } else {
656 log << MSG::ERROR << "Fail to retrieve EcRes_Dimu" << endreq;
657 }
658 }
659
660
661 // Efficiency and noise ratio
662 int part, segment, layer, stripMax;
663 part = segment = layer = stripMax = 0;
664 double digi[TAGN][3], effHit[TAGN][3], trackNum[TAGN][3], nosHit[TAGN][3];
665 double eff[TAGN], effErr[TAGN], nos[TAGN], nosErr[TAGN];
666 for(int i=0; i<TAGN; i++) {
667 for( int j=0; j<3; j++) {
668 digi[i][j] = effHit[i][j] = trackNum[i][j] = nosHit[i][j] = 0.;
669 }
670
671 eff[i] = effErr[i] = nos[i] = nosErr[i] = 0.;
672 }
673
674 // Layer efficiency
675 for( int i=0; i<LAYER_MAX; i++ )
676 {
677 effHit[0][0] = trackNum[0][0] = 0;
678 effHit[1][0] = trackNum[1][0] = 0;
679
680 for( int j=0; j<PART_MAX; j++ ) {
681 for( int k=0; k<SEGMENT_MAX; k++) {
682 stripMax = m_ptrIdTr->GetStripMax( j, k, i );
683 for( int n=0; n<stripMax; n++ ) {
684 trackNum[0][0] += m_recordAll[j][k][i][n][1];
685 effHit[0][0] += m_recordAll[j][k][i][n][2];
686
687 trackNum[1][0] += m_recordDimu[j][k][i][n][1];
688 effHit[1][0] += m_recordDimu[j][k][i][n][2];
689 }
690 }
691 }
692
693 if( trackNum[0][0] == 0 ) {
694 eff[0] = effErr[0] = 0.0;
695 } else {
696 eff[0] = ( (double)effHit[0][0] )/trackNum[0][0];
697 effErr[0] = sqrt( eff[0]*(1-eff[0])/trackNum[0][0] );
698 }
699 if(m_thsvc->getHist("/DQAHist/MUC/LayerEff_All", htmp).isSuccess()) {
700 htmp->Fill( i, eff[0] );
701 htmp->SetBinError( i+1, effErr[0] );
702 } else {
703 log << MSG::ERROR << "Fail to retrieve LayerEff_All" << endreq;
704 }
705 log << MSG::INFO << "LayerEff:\t" << i << "\t" << eff[0] << "\t" << effErr[0] << endreq;
706
707 if( trackNum[1][0] == 0 ) {
708 eff[1] = effErr[1] = 0.0;
709 } else {
710 eff[1] = ( (double)effHit[1][0] )/trackNum[1][0];
711 effErr[1] = sqrt( eff[1]*(1-eff[1])/trackNum[1][0] );
712 }
713 if(m_thsvc->getHist("/DQAHist/MUC/LayerEff_Dimu", htmp).isSuccess()) {
714 htmp->Fill( i, eff[1] );
715 htmp->SetBinError( i+1, effErr[1] );
716 } else {
717 log << MSG::ERROR << "Fail to retrieve LayerEff_Dimu" << endreq;
718 }
719
720 } // End layer
721
722 // Box efficiency and noise ratio, strip efficiency
723 for( int i=0; i<BOX_MAX; i++ )
724 {
725 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
726 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
727
728 digi[0][1] = effHit[0][1] = trackNum[0][1] = nosHit[0][1] = 0;
729 digi[1][1] = effHit[1][1] = trackNum[1][1] = nosHit[1][1] = 0;
730 for( int j=0; j<stripMax; j++ )
731 {
732 // For box
733 digi[0][1] += m_recordAll[part][segment][layer][j][0];
734 trackNum[0][1] += m_recordAll[part][segment][layer][j][1];
735 effHit[0][1] += m_recordAll[part][segment][layer][j][2];
736 nosHit[0][1] += m_recordAll[part][segment][layer][j][3];
737
738 digi[1][1] += m_recordDimu[part][segment][layer][j][0];
739 trackNum[1][1] += m_recordDimu[part][segment][layer][j][1];
740 effHit[1][1] += m_recordDimu[part][segment][layer][j][2];
741 nosHit[1][1] += m_recordDimu[part][segment][layer][j][3];
742
743 // For strip
744 trackNum[0][2] = m_recordAll[part][segment][layer][j][1];
745 effHit[0][2] = m_recordAll[part][segment][layer][j][2];
746 trackNum[1][2] = m_recordDimu[part][segment][layer][j][1];
747 effHit[1][2] = m_recordDimu[part][segment][layer][j][2];
748
749 if( trackNum[0][2] == 0 ) {
750 eff[0] = 0.0;
751 } else {
752 eff[0] = ( (double)effHit[0][2] )/trackNum[0][2];
753 }
754 if(m_thsvc->getHist("/DQAHist/MUC/StripEff_All", htmp).isSuccess()) {
755 htmp->Fill( eff[0] );
756 } else {
757 log << MSG::ERROR << "Fail to retrieve StripEff_All" << endreq;
758 }
759
760 if( trackNum[1][2] == 0 ) {
761 eff[1] = 0.0;
762 } else {
763 eff[1] = ( (double)effHit[1][2] )/trackNum[1][2];
764 }
765 if(m_thsvc->getHist("/DQAHist/MUC/StripEff_Dimu", htmp).isSuccess()) {
766 htmp->Fill( eff[1] );
767 } else {
768 log << MSG::ERROR << "Fail to retrieve StripEff_Dimu" << endreq;
769 }
770 } // End StripMax
771
772 // Box efficiency
773 if( trackNum[0][1] == 0 ) {
774 eff[0] = effErr[0] = 0.0;
775 } else {
776 eff[0] = ( (double)effHit[0][1] )/trackNum[0][1];
777 effErr[0] = sqrt( eff[0]*(1-eff[0])/trackNum[0][1] );
778 }
779 if(m_thsvc->getHist("/DQAHist/MUC/BoxEff_All", htmp).isSuccess()) {
780 htmp->Fill( i, eff[0] );
781 htmp->SetBinError( i+1, effErr[0] );
782 } else {
783 log << MSG::ERROR << "Fail to retrieve BoxEff_All" << endreq;
784 }
785 log << MSG::INFO << "BoxEff:\t" << i << "\t" << eff[0] << "\t" << effErr[0] << endreq;
786
787 if( trackNum[1][1] == 0 ) {
788 eff[1] = effErr[1] = 0.0;
789 } else {
790 eff[1] = ( (double)effHit[1][1] )/trackNum[1][1];
791 effErr[1] = sqrt( eff[1]*(1-eff[1])/trackNum[1][1] );
792 }
793 if(m_thsvc->getHist("/DQAHist/MUC/BoxEff_Dimu", htmp).isSuccess()) {
794 htmp->Fill( i, eff[1] );
795 htmp->SetBinError( i+1, effErr[1] );
796 } else {
797 log << MSG::ERROR << "Fail to retrieve BoxEff_Dimu" << endreq;
798 }
799
800 // Box noise ratio
801 if( digi[0][1] == 0 ) {
802 nos[0] = nosErr[0] = 0.0;
803 } else {
804 nos[0] = ( (double)nosHit[0][1] )/digi[0][1];
805 nosErr[0] = sqrt( nos[0]*(1-nos[0])/digi[0][1] );
806 }
807 if(m_thsvc->getHist("/DQAHist/MUC/BoxNosRatio_All", htmp).isSuccess()) {
808 htmp->Fill( i, nos[0] );
809 htmp->SetBinError( i+1, nosErr[0] );
810 } else {
811 log << MSG::ERROR << "Fail to retrieve BoxNosRatio_All" << endreq;
812 }
813 log << MSG::INFO << "BoxNosRatio:\t" << i << "\t" << nos[0] << "\t" << nosErr[0] << endreq;
814
815 if( digi[1][1] == 0 ) {
816 nos[1] = nosErr[1] = 0.0;
817 } else {
818 nos[1] = ( (double)nosHit[1][1] )/digi[1][1];
819 nosErr[1] = sqrt( nos[1]*(1-nos[1])/digi[1][1] );
820 }
821 if(m_thsvc->getHist("/DQAHist/MUC/BoxNosRatio_Dimu", htmp).isSuccess()) {
822 htmp->Fill( i, nos[1] );
823 htmp->SetBinError( i+1, nosErr[1] );
824 } else {
825 log << MSG::ERROR << "Fail to retrieve BoxNosRatio_Dimu" << endreq;
826 }
827 }// End BOX_MAX
828*/
829
830 log << MSG::INFO << "Finalize successfully! " << endmsg;
831
832 return StatusCode::SUCCESS;
833}
834
835StatusCode DQA_MUC::FillHistograms(bool isDimu) {
836
837 MsgStream log(msgSvc(), name());
838 log << MSG::INFO << "Filling histograms" << endmsg;
839 char name[100];
840 int part, segment, layer, strip, boxId, strId;
841 part = segment = layer = strip = boxId = strId = 0;
842 TH1 *hEff[2][2][3];
843
844 for(int i=0; i<TAGN; i++) {
845 for(int j=0; j<2; j++) {
846 for(int k=0; k<LVLN; k++)
847 {
848 sprintf( name, "/DQAHist/MUC/%sEff_%s_%s", LNAME[k], ENAME[i], HTYPE[j] );
849 if(!m_thsvc->getHist(name, hEff[i][j][k]).isSuccess())
850 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
851 }
852 }
853 }
854
855 // Numberator of eff
856 for(int i=0; i<m_effHitCol.size(); i++)
857 {
858 part = m_effHitCol[i]->Part();
859 segment = m_effHitCol[i]->Segment();
860 layer = m_effHitCol[i]->Layer();
861 strip = m_effHitCol[i]->Strip();
862 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
863 strId = m_ptrIdTr->GetStripId(part, segment, layer, strip);
864 hEff[0][0][0]->Fill(layer);
865 hEff[0][0][1]->Fill(boxId);
866 hEff[0][0][2]->Fill(strId);
867
868 if( isDimu ) {
869 hEff[1][0][0]->Fill(layer);
870 hEff[1][0][1]->Fill(boxId);
871 hEff[1][0][2]->Fill(strId);
872 }
873 }
874
875 // Denominator of eff
876 for(int i=0; i<m_expHitCol.size(); i++)
877 {
878 part = m_expHitCol[i]->Part();
879 segment = m_expHitCol[i]->Segment();
880 layer = m_expHitCol[i]->Layer();
881 strip = m_expHitCol[i]->Strip();
882 //cout<<"ExpHit:\t"<<part<<"\t"<<segment<<"\t"<<layer<<"\t"<<strip<<endl;
883 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
884 strId = m_ptrIdTr->GetStripId(part, segment, layer, strip);
885 hEff[0][1][0]->Fill(layer);
886 hEff[0][1][1]->Fill(boxId);
887 hEff[0][1][2]->Fill(strId);
888
889 if( isDimu ) {
890 hEff[1][1][0]->Fill(layer);
891 hEff[1][1][1]->Fill(boxId);
892 hEff[1][1][2]->Fill(strId);
893 }
894
895 }
896
897 TH1 *hNosRatio[2][2];
898 for(int i=0; i<TAGN; i++) {
899 for(int j=0; j<2; j++) {
900 sprintf( name, "/DQAHist/MUC/%sNosRatio_%s_%s", LNAME[1], ENAME[i], HTYPE[j] );
901 if(!m_thsvc->getHist(name, hNosRatio[i][j]).isSuccess())
902 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
903 }
904 }
905
906 // Numberator of box noise ratio
907 for(int i=0; i<m_nosHitCol.size(); i++)
908 {
909 part = m_nosHitCol[i]->Part();
910 segment = m_nosHitCol[i]->Segment();
911 layer = m_nosHitCol[i]->Layer();
912 strip = m_nosHitCol[i]->Strip();
913 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
914 strId = m_ptrIdTr->GetStripId(part, segment, layer, strip);
915 hNosRatio[0][0]->Fill(boxId);
916
917 if( isDimu ) hNosRatio[1][0]->Fill(boxId);
918 }
919
920 // Denominator of box noise ratio
921 for(int i=0; i<m_digiCol.size(); i++)
922 {
923 part = m_digiCol[i]->Part();
924 segment = m_digiCol[i]->Segment();
925 layer = m_digiCol[i]->Layer();
926 strip = m_digiCol[i]->Strip();
927 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
928 hNosRatio[0][1]->Fill(boxId);
929
930 if( isDimu ) hNosRatio[1][1]->Fill(boxId);
931 }
932
933 return StatusCode::SUCCESS;
934}
935
936//END
937
n
const Int_t n
Definition: DataBase/tau_mode.c:65
title
titledef title[20]
Definition: EvtPhokharaDef.hh:58
ENAME
const char ENAME[TAGN][10]
Definition: InstallArea/include/DQA_MUC/DQA_MUC/DQA_MUC.h:18
HTYPE
const char HTYPE[2][10]
Definition: InstallArea/include/DQA_MUC/DQA_MUC/DQA_MUC.h:19
LNAME
const char LNAME[LVLN][10]
Definition: InstallArea/include/DQA_MUC/DQA_MUC/DQA_MUC.h:20
TAGN
const int TAGN
Definition: InstallArea/include/DQA_MUC/DQA_MUC/DQA_MUC.h:16
LVLN
const int LVLN
Definition: InstallArea/include/DQA_MUC/DQA_MUC/DQA_MUC.h:17
RecMucTrackCol
ObjectVector< RecMucTrack > RecMucTrackCol
Definition: InstallArea/include/MucRecEvent/MucRecEvent/RecMucTrack.h:394
PI
const double PI
Definition: PipiJpsi.cxx:55
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:598
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
MucIdTransform::GetBoxId
int GetBoxId(int part, int segment, int layer)
Definition: Muc/MucCalibAlg/MucCalibAlg-00-02-16/src/MucIdTransform.cxx:51
MucMark::IsInCol
int IsInCol(int part, int segment, int layer, int strip, mark_col &aCol)
Definition: MucMark.cxx:100
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)