311 {
312
313
314
315 MsgStream log(
msgSvc(), name());
316 log << MSG::INFO << "in execute()" << endreq;
317
318 if( m_writeDst) m_subalg8->execute();
319 if( m_writeRec) m_subalg9->execute();
320
321 m_isBarrelBhabha = false;
322 m_isEndcapBhabha = false;
323 m_isBarrelDimu = false;
324 m_isEndcapDimu = false;
325 m_isHadron = false;
326 m_isBarrelDiphoton = false;
327 m_isEndcapDiphoton = false;
328
329 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
330 if(!eventHeader)
331 {
332 cout<<" eventHeader "<<endl;
333 return StatusCode::FAILURE;
334 }
335
336 int run=eventHeader->runNumber();
337 int event=eventHeader->eventNumber();
338
339 if(m_events%1000==0) cout<< m_events << " -------- run,event: "<<run<<","<<event<<endl;
340 m_events++;
341
343 if(!evtRecEvent ) {
344 cout<<" evtRecEvent "<<endl;
345 return StatusCode::FAILURE;
346 }
347
348 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
349 << evtRecEvent->totalCharged() << " , "
350 << evtRecEvent->totalNeutral() << " , "
351 << evtRecEvent->totalTracks() <<endreq;
353 if(!evtRecTrkCol){
354 cout<<" evtRecTrkCol "<<endl;
355 return StatusCode::FAILURE;
356 }
357
358 if(evtRecEvent->totalTracks()!=evtRecTrkCol->size()) return StatusCode::SUCCESS;
359
360
362 iGood.clear();
363 int nCharge = 0;
364 for(int i = 0;i < evtRecEvent->totalCharged(); i++)
365 {
366
368 if(!(*itTrk)->isMdcTrackValid()) continue;
370 double vx0 = mdcTrk->
x();
371 double vy0 = mdcTrk->
y();
372 double vz0 = mdcTrk->
z();
373 double vr0 = mdcTrk->
r();
374 double theta0 = mdcTrk->
theta();
375 double p0 = mdcTrk->
p();
376 double pt0 = mdcTrk->
pxy();
377
378 if(m_output) {
379 m_vx0 = vx0;
380 m_vy0 = vy0;
381 m_vz0 = vz0;
382 m_vr0 = vr0;
383 m_theta0 = theta0;
384 m_p0 = p0;
385 m_pt0 = pt0;
386 m_tuple1->write();
387 }
388
389 if(fabs(vz0) >= m_vz0cut) continue;
390 if(vr0 >= m_vr0cut) continue;
391 if(pt0 >= m_pt0HighCut) continue;
392 if(pt0 <= m_pt0LowCut) continue;
393
394 iGood.push_back((*itTrk)->trackId());
395 nCharge += mdcTrk->
charge();
396 }
397 int nGood = iGood.size();
398
399
401 iGam.clear();
402 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
403
405 if(!(*itTrk)->isEmcShowerValid()) continue;
407 double eraw = emcTrk->
energy();
408 if(m_output) {
409 m_eraw = eraw;
410 m_tuple2->write();
411 }
412 if(eraw < m_energyThreshold) continue;
413 iGam.push_back((*itTrk)->trackId());
414 }
415 int nGam = iGam.size();
416
417
419 ipip.clear();
420 ipim.clear();
422 ppip.clear();
423 ppim.clear();
424
425
426 double echarge = 0.;
427 double ptot = 0.;
429 double eemc = 0.;
430 double pp = 0.;
431 double pm = 0.;
432
433 for(int i = 0; i < nGood; i++) {
435
436 if((*itTrk)->isMdcTrackValid()) {
437
438
439
440
441
443
445
446 HepLorentzVector ptrk;
447 ptrk.setPx(mdcTrk->
px());
448 ptrk.setPy(mdcTrk->
py());
449 ptrk.setPz(mdcTrk->
pz());
450 double p3 = ptrk.mag();
451 ptrk.setE(sqrt(p3*p3+
mpi*
mpi));
452 ptrk = ptrk.boost(-0.011,0,0);
453
454 echarge += ptrk.e();
456
458 ppip.push_back(ptrk);
459 pp = ptrk.rho();
460 } else {
461 ppim.push_back(ptrk);
462 pm = ptrk.rho();
463 }
464
465 }
466
467 if((*itTrk)->isEmcShowerValid()) {
468
470 double eraw = emcTrk->
energy();
471 double phi = emcTrk->
phi();
472 double the = emcTrk->
theta();
473
474 HepLorentzVector ptrk;
475 ptrk.setPx(eraw*
sin(the)*
cos(phi));
476 ptrk.setPy(eraw*
sin(the)*
sin(phi));
477 ptrk.setPz(eraw*
cos(the));
478 ptrk.setE(eraw);
479 ptrk = ptrk.boost(-0.011,0,0);
480
481 eemc += ptrk.e();
483
484 }
485
486 }
487
488
489
490
492 pGam.clear();
493 double eneu=0;
494 for(int i = 0; i < nGam; i++) {
497 double eraw = emcTrk->
energy();
498 double phi = emcTrk->
phi();
499 double the = emcTrk->
theta();
500 HepLorentzVector ptrk;
501 ptrk.setPx(eraw*
sin(the)*
cos(phi));
502 ptrk.setPy(eraw*
sin(the)*
sin(phi));
503 ptrk.setPz(eraw*
cos(the));
504 ptrk.setE(eraw);
505 ptrk = ptrk.boost(-0.011,0,0);
506 pGam.push_back(ptrk);
507 eneu += ptrk.e();
509 eemc += ptrk.e();
510 }
511
512 double esum = echarge + eneu;
513
514
515 double maxE = 0.;
516 double secE = 0.;
517 double maxThe = 999.;
518 double maxPhi = 999.;
519 double secThe = 999.;
520 double secPhi = 999.;
521 int npart = 999.;
522 double dphi = 999.;
523 double dthe = 999.;
524 int mdcHit1 = 0.;
525 int mdcHit2 = 0.;
526
527 SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
528 if (!aShowerCol) {
529 log << MSG::WARNING << "Could not find RecEmcShowerCol" << endreq;
530 return( StatusCode::SUCCESS);
531 }
532
533 int ishower = 0;
534 RecEmcShowerCol::iterator iShowerCol;
535 for(iShowerCol=aShowerCol->begin();
536 iShowerCol!=aShowerCol->end();
537 iShowerCol++) {
538
539 if(ishower == 0) {
540 maxE = (*iShowerCol)->e5x5();
541 maxThe = (*iShowerCol)->theta();
542 maxPhi = (*iShowerCol)->phi();
543 npart = (*iShowerCol)->module();
544 } else if(ishower == 1) {
545 secE = (*iShowerCol)->e5x5();
546 secThe = (*iShowerCol)->theta();
547 secPhi = (*iShowerCol)->phi();
548 } else if(ishower == 2) {
549 break;
550 }
551
552 ishower++;
553
554 }
555
556 if(aShowerCol->size() >= 2) {
557
558 dphi = (fabs(maxPhi-secPhi)-CLHEP::pi)*180./CLHEP::pi;
559 dthe = (fabs(maxThe+secThe)-CLHEP::pi)*180./CLHEP::pi;
560
561 double phi1 = maxPhi<0 ? maxPhi+CLHEP::twopi : maxPhi;
562 double phi2 = secPhi<0 ? secPhi+CLHEP::twopi : secPhi;
563
564
567 double phi12=(phi11+phi22-CLHEP::pi)*0.5;
568 double phi21=(phi11+phi22+CLHEP::pi)*0.5;
569 if(phi12<0.) phi12 += CLHEP::twopi;
570 if(phi21>CLHEP::twopi) phi21 -= CLHEP::twopi;
571
572 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc(), "/Event/Digi/MdcDigiCol");
573 if (!mdcDigiCol) {
574 log << MSG::FATAL << "Could not find MdcDigiCol!" << endreq;
575 return StatusCode::FAILURE;
576 }
577 int hitnum = mdcDigiCol->size();
578 for (int i = 0;i< hitnum;i++ ) {
579 MdcDigi* digi=
dynamic_cast<MdcDigi*
>(mdcDigiCol->containedObject(i));
582 if (
time == 0x7FFFFFFF || charge == 0x7FFFFFFF)
continue;
586 if(ilayer>=43)
587 log << MSG::ERROR << "MDC(" << ilayer <<","<<iphi <<")"<<endreq;
588 double phi=CLHEP::twopi*iphi/idmax[ilayer];
591
592
593 }
594 }
595
596
597
598
599
600
601 if( eemc>m_bhabhaEmcECut && maxE>m_bhabhaMaxECut && secE>m_bhabhaSecECut
602 &&
abs(dthe)<m_bhabhaDTheCut &&
603 ( (dphi>m_bhabhaDPhiCut1&&dphi<m_bhabhaDPhiCut2) ||
604 (dphi>m_bhabhaDPhiCut3&&dphi<m_bhabhaDPhiCut4) ) ) {
605 if( npart==1 && mdcHit1>m_bhabhaMdcHitCutB && mdcHit2>m_bhabhaMdcHitCutB ) {
606 m_isBarrelBhabha = true;
607 m_barrelBhabhaNumber++;
608 } else if( ( npart==0 || npart==2 )
609 && mdcHit1>m_bhabhaMdcHitCutE && mdcHit2>m_bhabhaMdcHitCutE ) {
610 m_isEndcapBhabha = true;
611 m_endcapBhabhaNumber++;
612 }
613 }
614
615
616
617
618
619
620
621
622
623
624
625
626
627 if(m_selectDimu) {
628 if( m_dimuAlg->
IsDimu() == 1 ) {
629 m_isBarrelDimu = true;
630 m_barrelDimuNumber++;
631 }
else if(m_dimuAlg->
IsDimu() == 2 ) {
632 m_isEndcapDimu = true;
633 m_endcapDimuNumber++;
634 }
635 }
636
637
638 if( (nGood>=1 && esum>m_hadronChaECut)
639 || (nGood==0 && esum>m_hadronNeuECut) ) {
640 m_isHadron = true;
641 m_hadronNumber++;
642 }
643
644
645 if( nGood==0 && eemc>m_diphotonEmcECut && secE>m_diphotonSecECut
646 && fabs(dthe)<m_diphotonDTheCut
647 && dphi>m_diphotonDPhiCut1 && dphi<m_diphotonDPhiCut2 ) {
648 if( npart==1 ) {
649 m_isBarrelDiphoton = true;
650 m_barrelDiphotonNumber++;
651 } else if( npart==0 || npart==2 ) {
652 m_isEndcapDiphoton = true;
653 m_endcapDiphotonNumber++;
654 }
655 }
656
657
658 if( m_selectBhabha && m_isBarrelBhabha ) m_subalg1->execute();
659 if( m_selectBhabha && m_isEndcapBhabha ) m_subalg2->execute();
660 if( m_selectDimu && m_isBarrelDimu ) m_subalg3->execute();
661 if( m_selectDimu && m_isEndcapDimu ) m_subalg4->execute();
662 if( m_selectHadron && m_isHadron ) m_subalg5->execute();
663 if( m_selectDiphoton && m_isBarrelDiphoton ) m_subalg6->execute();
664 if( m_selectDiphoton && m_isEndcapDiphoton ) m_subalg7->execute();
665
666
667 if(m_output) {
668 m_runnb = run;
669 m_evtnb = event;
670 m_esum = esum;
671 m_eemc = eemc;
673 m_nCharge = nCharge;
674 m_nGood = nGood;
675 m_nGam = nGam;
676 m_ptot = ptot;
677 m_pp = pp;
678 m_pm = pm;
679 m_maxE = maxE;
680 m_secE = secE;
681 m_dThe = dthe;
682 m_dPhi = dphi;
683 m_mdcHit1 = mdcHit1;
684 m_mdcHit2 = mdcHit2;
685 m_tuple0->write();
686 }
687
688 return StatusCode::SUCCESS;
689
690}
double sin(const BesAngle a)
double cos(const BesAngle a)
double abs(const EvtComplex &c)
EvtRecTrackCol::iterator EvtRecTrackIterator
std::vector< HepLorentzVector > Vp4
const double theta() const
bool WhetherSector(double ph, double ph1, double ph2)
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
static int wire(const Identifier &id)
virtual Identifier identify() const
unsigned int getChargeChannel() const
unsigned int getTimeChannel() const
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol