219 {
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
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();
265 eventDigi++;
266
267
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
280 m_clusterCol = ( *m_ptrMucMark ).CreateClusterCol( 1, m_digiCol );
281 log << MSG::INFO << "Total clusters in this event: " << m_clusterCol.size() << endmsg;
282
283
284
286 SmartDataPtr<RecMucTrackCol> mucTrackCol( eventSvc(), "/Event/Recon/RecMucTrackCol" );
287 if ( !mucTrackCol )
288 {
289
290 log << MSG::WARNING << "Could not find RecMucTrackCol"
291 << endmsg;
292 return ( StatusCode::SUCCESS );
293 }
294
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
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
357 mucExpHitCol = ( *trackIter )->GetExpectedHits();
358
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
372 int isInPos = -1;
373 bool isInEffWindow = false;
374 isInPos = currentMark->
IsInCol( m_segDigiCol[part][segment] );
375
376
377 if ( part ==
BRID && ( layer - lastLayerBR > 1 ) )
continue;
378 if ( part !=
BRID && ( layer - lastLayerEC > 1 ) )
continue;
379
380
381 if ( part ==
BRID && layer == 0 && ( strip < 2 || strip > 45 ) )
382 {
383 if ( isInPos != -1 )
384 {
385 m_recordAll[part][segment][layer][strip][2]++;
386 m_recordAll[part][segment][layer][strip][1]++;
387 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
388
389 if ( isDimu )
390 {
391 m_recordDimu[part][segment][layer][strip][2]++;
392 m_recordDimu[part][segment][layer][strip][1]++;
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;
401 }
402
403 if ( isInPos != -1 )
404 {
405 m_recordAll[part][segment][layer][strip][2]++;
406 m_recordAll[part][segment][layer][strip][1]++;
407 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
408
409 if ( isDimu )
410 {
411 m_recordDimu[part][segment][layer][strip][2]++;
412 m_recordDimu[part][segment][layer][strip][1]++;
413 }
414
415 continue;
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]++;
430 m_recordAll[part][segment][layer][tempStrip][1]++;
431 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
432
433 if ( isDimu )
434 {
435 m_recordDimu[part][segment][layer][tempStrip][2]++;
436 m_recordDimu[part][segment][layer][tempStrip][1]++;
437 }
438
439 isInEffWindow = true;
440 }
441 }
442
443 if ( isInEffWindow ) { continue; }
444 else
445 {
446 m_recordAll[part][segment][layer][strip][1]++;
447 if ( isDimu ) m_recordDimu[part][segment][layer][strip][1]++;
448 }
449
450 }
451
452
453 log << MSG::INFO << "Fill residual" << endmsg;
454 if ( !m_dstFileOnly )
455 {
456 vector<float> m_lineResCol = ( *trackIter )->getDistHits();
457
458 mucRawHitCol = ( *trackIter )->GetHits();
459
460
461 if ( trackHitNum > 4 && m_lineResCol[0] != -99 )
462 {
463
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
486 if ( firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false )
487 {
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
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
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
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
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" << 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;
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 {
559 isNosHit = false;
560 break;
561 }
562 }
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
572
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 }
578 }
579
580
581 DQA_MUC::FillHistograms( isDimu ).ignore();
582
583
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}
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