227 {
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
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
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();
268 eventDigi ++;
269
270
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
281 m_clusterCol = (*m_ptrMucMark).CreateClusterCol(1, m_digiCol );
282 log << MSG::INFO << "Total clusters in this event: " << m_clusterCol.size() << endreq;
283
284
285
287 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc(),"/Event/Recon/RecMucTrackCol");
288 if (!mucTrackCol) {
289
290 log << MSG::WARNING << "Could not find RecMucTrackCol" << endreq;
291 return( StatusCode::SUCCESS);
292 }
293
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
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
359
360
361
362 mucExpHitCol = (*trackIter)->GetExpectedHits();
363
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
375 int isInPos = -1;
376 bool isInEffWindow = false;
377 isInPos = currentMark->
IsInCol( m_segDigiCol[part][segment] );
378
379
380 if( part == BRID && (layer-lastLayerBR>1) ) continue;
381 if( part != BRID && (layer-lastLayerEC>1) ) continue;
382
383
384 if( part==BRID && layer==0 && (strip<2 || strip>45) )
385 {
386 if( isInPos != -1)
387 {
388 m_recordAll[part][segment][layer][strip][2] ++;
389 m_recordAll[part][segment][layer][strip][1] ++;
390 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
391
392 if(isDimu) {
393 m_recordDimu[part][segment][layer][strip][2] ++;
394 m_recordDimu[part][segment][layer][strip][1] ++;
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;
402 }
403
404 if( isInPos != -1 )
405 {
406 m_recordAll[part][segment][layer][strip][2] ++;
407 m_recordAll[part][segment][layer][strip][1] ++;
408 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
409
410 if(isDimu) {
411 m_recordDimu[part][segment][layer][strip][2] ++;
412 m_recordDimu[part][segment][layer][strip][1] ++;
413 }
414
415 continue;
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] ++;
427 m_recordAll[part][segment][layer][tempStrip][1] ++;
428 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
429
430 if(isDimu) {
431 m_recordDimu[part][segment][layer][tempStrip][2] ++;
432 m_recordDimu[part][segment][layer][tempStrip][1] ++;
433 }
434
435 isInEffWindow = true;
436 }
437 }
438
439 if( isInEffWindow ) { continue; }
440 else {
441 m_recordAll[part][segment][layer][strip][1] ++;
442 if(isDimu) m_recordDimu[part][segment][layer][strip][1] ++;
443 }
444
445 }
446
447
448 log << MSG::INFO << "Fill residual" << endreq;
449 vector<float> m_lineResCol = (*trackIter)->getDistHits();
450
451 mucRawHitCol = (*trackIter)->GetHits();
452
453
454 if( trackHitNum > 4 && m_lineResCol[0] != -99)
455 {
456
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
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
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
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
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
515 }
516 }
517 }
518
519 firedFlag[part][layer][1] = true;
520 }
521 }
522
523 mucRawHitCol.clear();
524 mucExpHitCol.clear();
525
526 }
527
528
529
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;
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 {
547 hasEffHit = true;
548 break;
549 }
550 }
551
552 if ( hasEffHit && (m_digiCol[i]->IsInCol( m_clusterCol[j] ) != -1) ) {
553 isNosHit = false;
554 break;
555 }
556 }
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
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 }
568 }
569
570
571 DQA_MUC::FillHistograms( isDimu );
572
573
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
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}
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.
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)
_EXTERN_ std::string EvtRecEvent