228 {
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
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
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();
271 eventDigi ++;
272
273
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
284 m_clusterCol = (*m_ptrMucMark).CreateClusterCol(1, m_digiCol );
285 log << MSG::INFO << "Total clusters in this event: " << m_clusterCol.size() << endreq;
286
287
288
290 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc(),"/Event/Recon/RecMucTrackCol");
291 if (!mucTrackCol) {
292
293 log << MSG::WARNING << "Could not find RecMucTrackCol" << endreq;
294 return( StatusCode::SUCCESS);
295 }
296
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];
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
331 CLHEP::Hep3Vector a3Vector((*trackIter)->xPos(),(*trackIter)->yPos(),(*trackIter)->zPos());
333 phi = a3Vector.phi();
334 if(m_thsvc->getHist("/DQAHist/MUC/Costheta_All", htmp).isSuccess()) {
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()) {
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
362
363
364
365 mucExpHitCol = (*trackIter)->GetExpectedHits();
366
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
378 int isInPos = -1;
379 bool isInEffWindow = false;
380 isInPos = currentMark->
IsInCol( m_segDigiCol[part][segment] );
381
382
383 if( part == BRID && (layer-lastLayerBR>1) ) continue;
384 if( part != BRID && (layer-lastLayerEC>1) ) continue;
385
386
387 if( part==BRID && layer==0 && (strip<2 || strip>45) )
388 {
389 if( isInPos != -1)
390 {
391 m_recordAll[part][segment][layer][strip][2] ++;
392 m_recordAll[part][segment][layer][strip][1] ++;
393 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
394
395 if(isDimu) {
396 m_recordDimu[part][segment][layer][strip][2] ++;
397 m_recordDimu[part][segment][layer][strip][1] ++;
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;
405 }
406
407 if( isInPos != -1 )
408 {
409 m_recordAll[part][segment][layer][strip][2] ++;
410 m_recordAll[part][segment][layer][strip][1] ++;
411 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
412
413 if(isDimu) {
414 m_recordDimu[part][segment][layer][strip][2] ++;
415 m_recordDimu[part][segment][layer][strip][1] ++;
416 }
417
418 continue;
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] ++;
430 m_recordAll[part][segment][layer][tempStrip][1] ++;
431 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
432
433 if(isDimu) {
434 m_recordDimu[part][segment][layer][tempStrip][2] ++;
435 m_recordDimu[part][segment][layer][tempStrip][1] ++;
436 }
437
438 isInEffWindow = true;
439 }
440 }
441
442 if( isInEffWindow ) { continue; }
443 else {
444 m_recordAll[part][segment][layer][strip][1] ++;
445 if(isDimu) m_recordDimu[part][segment][layer][strip][1] ++;
446 }
447
448 }
449
450
451 log << MSG::INFO << "Fill residual" << endreq;
452 if(!m_dstFileOnly)
453 {
454 vector<float> m_lineResCol = (*trackIter)->getDistHits();
455
456 mucRawHitCol = (*trackIter)->GetHits();
457
458
459 if( trackHitNum > 4 && m_lineResCol[0] != -99)
460 {
461
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
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
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
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
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
520 }
521 }
522 }
523
524 firedFlag[part][layer][1] = true;
525 }
526 }
527 }
528 mucRawHitCol.clear();
529 mucExpHitCol.clear();
530
531 }
532
533
534
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;
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 {
552 hasEffHit = true;
553 break;
554 }
555 }
556
557 if ( hasEffHit && (m_digiCol[i]->IsInCol( m_clusterCol[j] ) != -1) ) {
558 isNosHit = false;
559 break;
560 }
561 }
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
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 }
573 }
574
575
576 DQA_MUC::FillHistograms( isDimu );
577
578
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
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}
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)
ObjectVector< RecMucTrack > RecMucTrackCol
static int barrel_ec(const Identifier &id)
Values of different levels.
static int layer(const Identifier &id)
static int channel(const Identifier &id)
static int segment(const Identifier &id)
int IsInCol(int part, int segment, int layer, int strip, mark_col &aCol)
int Part() const
Get Part.
int Strip() const
Get Strip.
_EXTERN_ std::string EvtRecEvent