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