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