308 {
309
310
311
312 MsgStream log(
msgSvc(), name() );
313 log << MSG::INFO << "in execute()" << endmsg;
314 m_cout_all++;
315 StatusCode sc = StatusCode::SUCCESS;
316
317 setFilterPassed( false );
318
319 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
320 if ( !eventHeader )
321 {
322 log << MSG::ERROR << "EventHeader not found" << endmsg;
323 return StatusCode::SUCCESS;
324 }
325 int run( eventHeader->runNumber() );
326 int event( eventHeader->eventNumber() );
327 if ( event % 1000 == 0 ) cout << "run: " << run << " event: " << event << endl;
328
329 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), "/Event/EvtRec/EvtRecEvent" );
330 if ( !evtRecEvent )
331 {
332 log << MSG::ERROR << "EvtRecEvent not found" << endmsg;
333 return StatusCode::SUCCESS;
334 }
335 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
336 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
337
338 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), "/Event/EvtRec/EvtRecTrackCol" );
339 if ( !evtRecTrkCol )
340 {
341 log << MSG::ERROR << "EvtRecTrackCol" << endmsg;
342 return StatusCode::SUCCESS;
343 }
344
345 if ( m_trigger_flag )
346 {
348 if ( !trigData )
349 {
350 log << MSG::FATAL << "Could not find Trigger Data for physics analysis" << endmsg;
351 return StatusCode::FAILURE;
352 }
353
354 log << MSG::DEBUG << "Trigger conditions: " << endmsg;
355 for ( int i = 0; i < 48; i++ )
356 {
357 log << MSG::DEBUG << "i:" << i << " name:" << trigData->getTrigCondName( i )
358 << " cond:" << trigData->getTrigCondition( i ) << endmsg;
359 }
360
361 int m_trig_tot( 0 ), m_trig_which( -1 );
362 if ( m_eventrate )
363 {
364 for ( int j = 0; j < 16; j++ )
365 {
366 if ( trigData->getTrigChannel( j ) )
367 {
368 m_trig_tot++;
369 m_trig_which = j + 1;
370 }
371 }
372 if ( m_trig_tot == 1 && m_trig_which == m_chan_det )
m_cout_everat++;
373 return sc;
374 }
375 }
376
378 if ( evtRecEvent->totalCharged() < 3 || evtRecTrkCol->size() < 3 ||
379 evtRecEvent->totalTracks() > 99 || evtRecTrkCol->size() > 99 )
380 return StatusCode::SUCCESS;
382
383
385 iGood.clear();
386 int m_num[4] = { 0, 0, 0, 0 };
387 int nCharge = 0;
388 m_pion_matched = 0;
389 m_lep_matched = 0;
390 HepLorentzVector m_lv_pionp, m_lv_pionm, m_lv_lepp, m_lv_lepm, m_lv_ele, m_lv_pos, m_lv_mum,
391 m_lv_mup;
392
393 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
394 {
396 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
397 RecMdcKalTrack* mdcTrk = ( *itTrk )->mdcKalTrack();
398
403 if ( fabs( m_vz0 ) >= m_vz0cut ) continue;
404 if ( m_vr0 >= m_vr0cut ) continue;
405 iGood.push_back( i );
406 nCharge += mdcTrk->
charge();
407 if ( mdcTrk->
p() < 1.0 )
408 {
409 if ( ( *itTrk )->isEmcShowerValid() ) m_pion_matched++;
410 }
411 else
412 {
413 if ( ( *itTrk )->isEmcShowerValid() ) m_lep_matched++;
414 }
415
416 if ( mdcTrk->
charge() > 0 )
417 {
418 if ( mdcTrk->
p() < 1.0 )
419 {
421
422
423
424 m_lv_pionp = mdcTrk->
p4(
xmass[2] );
425 m_num[0]++;
426 }
427 else
428 {
430 m_lv_pos = mdcTrk->
p4(
xmass[0] );
432 m_lv_mup = mdcTrk->
p4(
xmass[1] );
433 m_num[2]++;
434 }
435 }
436 else
437 {
438 if ( mdcTrk->
p() < 1.0 )
439 {
441 m_lv_pionm = mdcTrk->
p4(
xmass[2] );
442 m_num[1]++;
443 }
444 else
445 {
447 m_lv_ele = mdcTrk->
p4(
xmass[0] );
449 m_lv_mum = mdcTrk->
p4(
xmass[1] );
450 m_num[3]++;
451 }
452 }
453 }
454
455 int nGood = iGood.size();
456 log << MSG::DEBUG << "With KalmanTrack, ngood, totcharge = " << nGood << " , " << nCharge
457 << endmsg;
458 if ( nGood < 3 || nGood > 4 ) return sc;
460
461 m_ep_ratio = 0;
462 for ( int i = 0; i < evtRecEvent->totalTracks(); i++ )
463 {
465 if ( !( *itTrk )->isEmcShowerValid() ) continue;
466 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
467 m_ep_ratio += emcTrk->
energy();
468 }
469
470 if ( m_ep_ratio > m_distin_emuon )
471 {
472 m_lv_lepp = m_lv_pos;
473 m_lv_lepm = m_lv_ele;
474 }
475 else
476 {
477 m_lv_lepp = m_lv_mup;
478 m_lv_lepm = m_lv_mum;
479 }
480
481 HepLorentzVector m_lv_lab( 0.04, 0, 0, 3.686 );
482 if ( nGood == 4 )
483 {
484 if ( nCharge ) return sc;
485 m_event_flag = 4;
486 }
487 else
488 {
489 if ( m_num[0] > 1 || m_num[1] > 1 || m_num[2] > 1 || m_num[3] > 1 ) return sc;
490 if ( m_num[0] == 0 )
491 {
492 if ( nCharge != -1 ) return StatusCode::SUCCESS;
493 m_lv_pionp = m_lv_lab - m_lv_pionm - m_lv_lepp - m_lv_lepm;
494 if ( m_lv_pionp.vect().cosTheta() > m_cosThetaCut ) return StatusCode::SUCCESS;
495 m_event_flag = 0;
496 }
497 if ( m_num[1] == 0 )
498 {
499 if ( nCharge != 1 ) return StatusCode::SUCCESS;
500 m_lv_pionm = m_lv_lab - m_lv_pionp - m_lv_lepp - m_lv_lepm;
501 if ( m_lv_pionm.vect().cosTheta() > m_cosThetaCut ) return StatusCode::SUCCESS;
502 m_event_flag = 1;
503 }
504 if ( m_num[2] == 0 )
505 {
506 if ( nCharge != -1 ) return StatusCode::SUCCESS;
507 m_lv_lepp = m_lv_lab - m_lv_pionp - m_lv_pionm - m_lv_lepm;
508 if ( m_lv_lepp.vect().cosTheta() > m_cosThetaCut ) return StatusCode::SUCCESS;
509 m_event_flag = 2;
510 }
511 if ( m_num[3] == 0 )
512 {
513 if ( nCharge != 1 ) return StatusCode::SUCCESS;
514 m_lv_lepm = m_lv_lab - m_lv_pionp - m_lv_pionm - m_lv_lepp;
515 if ( m_lv_lepm.vect().cosTheta() > m_cosThetaCut ) return StatusCode::SUCCESS;
516 m_event_flag = 3;
517 }
518 }
520
521
522
523 HepLorentzVector m_lv_recoil, m_lv_jpsi;
524 m_lv_recoil = m_lv_lab - m_lv_pionp - m_lv_pionm;
525 m_lv_jpsi = m_lv_lepp + m_lv_lepm;
526
527 m_mass_twopi = ( m_lv_pionp + m_lv_pionm ).m();
528 m_mass_recoil = m_lv_recoil.m();
529 m_mass_jpsi = m_lv_jpsi.m();
530
531
532 if ( m_mass_recoil < 3.05 || m_mass_recoil > 3.15 ) return sc;
533 if ( m_mass_jpsi < 3.0 || m_mass_jpsi > 3.2 ) return sc;
535
536 HepLorentzVector m_ttm( m_lv_jpsi + m_lv_pionp + m_lv_pionm );
537 if ( m_ttm.m() > 4 || m_ttm.m() < 3 ) return sc;
538
539
540 m_pipi_dang = m_lv_pionp.vect().cosTheta( m_lv_pionm.vect() );
541
542 m_mom_pionp = m_lv_pionp.vect().mag();
543 m_mom_pionm = m_lv_pionm.vect().mag();
544 m_mom_lepp = m_lv_lepp.vect().mag();
545 m_mom_lepm = m_lv_lepm.vect().mag();
546 m_trans_ratio_lepp = m_lv_lepp.vect().perp() / m_lv_lepp.vect().mag();
547 m_trans_ratio_lepm = m_lv_lepm.vect().perp() / m_lv_lepm.vect().mag();
548 m_trans_ratio_pionp = m_lv_pionp.vect().perp() / m_lv_pionp.vect().mag();
549 m_trans_ratio_pionm = m_lv_pionm.vect().perp() / m_lv_pionm.vect().mag();
550
551 Hep3Vector m_boost_jpsi( m_lv_recoil.boostVector() );
552 HepLorentzVector m_lv_cms_lepp( boostOf( m_lv_lepp, -m_boost_jpsi ) );
553 HepLorentzVector m_lv_cms_lepm( boostOf( m_lv_lepm, -m_boost_jpsi ) );
554 m_cms_lepm = m_lv_cms_lepm.vect().mag();
555 m_cms_lepp = m_lv_cms_lepp.vect().mag();
556 log << MSG::DEBUG << "jpsi four momentum in cms " << m_lv_cms_lepp + m_lv_cms_lepm << endmsg;
557
558 m_inv_mass = m_ttm.m();
559 m_tot_e = m_ttm.e();
560 m_tot_px = m_ttm.px();
561 m_tot_py = m_ttm.py();
562 m_tot_pz = m_ttm.pz();
563 m_run = run;
564 m_event = event;
565 HepLorentzVector m_lv_book( 0, 0, 0, 0 );
566 for ( m_index = 0; m_index < 4; m_index++ )
567 {
568 switch ( m_index )
569 {
570 case 0: m_lv_book = m_lv_pionp; break;
571 case 1: m_lv_book = m_lv_pionm; break;
572 case 2: m_lv_book = m_lv_lepp; break;
573 case 3: m_lv_book = m_lv_lepm; break;
574 default: m_lv_book.setE( 2008 );
575 }
576 m_cos_theta[m_index] = m_lv_book.vect().cosTheta();
577 m_phi[m_index] = m_lv_book.vect().phi();
578 m_four_mom[m_index][0] = m_lv_book.e();
579 m_four_mom[m_index][1] = m_lv_book.px();
580 m_four_mom[m_index][2] = m_lv_book.py();
581 m_four_mom[m_index][3] = m_lv_book.pz();
582 }
583
584 if ( m_subsample_flag ) setFilterPassed( true );
585 else if ( m_mass_recoil > 3.08 && m_mass_recoil < 3.12 && m_mass_jpsi > 3.0 &&
586 m_mass_jpsi < 3.2 && m_cms_lepp < 1.7 && m_cms_lepp > 1.4 && m_cms_lepm < 1.7 &&
587 m_cms_lepm > 1.4 && m_event_flag == 4 && m_pipi_dang < m_pipi_dang_cut )
588 setFilterPassed( true );
589
590
591
592 SmartDataPtr<Event::McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
593 if ( m_run < 0 )
594 {
595 int m_numParticle( 0 ), m_true_pid( 0 );
596 if ( !mcParticleCol )
597 {
598 log << MSG::ERROR << "Could not retrieve McParticelCol" << endmsg;
599 return StatusCode::FAILURE;
600 }
601 else
602 {
603 bool psipDecay( false );
604 int rootIndex( -1 );
605 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
606 for ( ; iter_mc != mcParticleCol->end(); iter_mc++ )
607 {
608 if ( ( *iter_mc )->primaryParticle() ) continue;
609 if ( !( *iter_mc )->decayFromGenerator() ) continue;
610
611 if ( ( *iter_mc )->particleProperty() == 100443 )
612 {
613 psipDecay = true;
614 rootIndex = ( *iter_mc )->trackIndex();
615 }
616 if ( !psipDecay ) continue;
617 int mcidx = ( ( *iter_mc )->mother() ).trackIndex() - rootIndex;
618 int pdgid = ( *iter_mc )->particleProperty();
619 m_pdgid[m_numParticle] = pdgid;
620 m_motheridx[m_numParticle] = mcidx;
621 m_numParticle++;
622
623
624 if ( ( *iter_mc )->particleProperty() == 211 )
625 m_true_pionp = ( *iter_mc )->initialFourMomentum().vect().mag();
626 if ( ( *iter_mc )->particleProperty() == -211 )
627 m_true_pionm = ( *iter_mc )->initialFourMomentum().vect().mag();
628 }
629 m_idxmc = m_numParticle;
630 }
631 }
632
633 m_tuple1->write();
634 m_tuple8->write();
635
636
638 iGam.clear();
639 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
640 {
642 if ( !( *itTrk )->isEmcShowerValid() ) continue;
643 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
644 Hep3Vector emcpos( emcTrk->
x(), emcTrk->
y(), emcTrk->
z() );
645
646 double dthe = 200.;
647 double dphi = 200.;
648 double dang = 200.;
649 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
650 {
652 if ( !( *jtTrk )->isExtTrackValid() ) continue;
653 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
656
657 double angd = extpos.angle( emcpos );
658 double thed = extpos.theta() - emcpos.theta();
659 double phid = extpos.deltaPhi( emcpos );
660 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi +
pi, CLHEP::twopi ) - CLHEP::pi;
661 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi +
pi, CLHEP::twopi ) - CLHEP::pi;
662
663 if ( fabs( thed ) < fabs( dthe ) ) dthe = thed;
664 if ( fabs( phid ) < fabs( dphi ) ) dphi = phid;
665 if ( angd < dang ) dang = angd;
666 }
667 if ( dang >= 200 ) continue;
668 double eraw = emcTrk->
energy();
669 dthe = dthe * 180 / ( CLHEP::pi );
670 dphi = dphi * 180 / ( CLHEP::pi );
671 dang = dang * 180 / ( CLHEP::pi );
672 m_dthe = dthe;
673 m_dphi = dphi;
674 m_dang = dang;
675 m_eraw = eraw;
676
677
678
679 iGam.push_back( ( *itTrk )->trackId() );
680 }
681
682 m_nGam = iGam.size();
683 log << MSG::DEBUG << "num Good Photon " << m_nGam << " , " << evtRecEvent->totalNeutral()
684 << endmsg;
685 m_tuple2->write();
686
687
688
689
690 if ( m_checkDedx )
691 {
692 int m_dedx_cout( 0 );
693 for ( int i = 0; i < nGood; i++ )
694 {
696 if ( !( *itTrk )->isMdcDedxValid() ) continue;
697 RecMdcKalTrack* mdcTrk = ( *itTrk )->mdcKalTrack();
698 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
699
700 m_ptrk = mdcTrk->
p();
701 m_chie = dedxTrk->
chiE();
702 m_chimu = dedxTrk->
chiMu();
703 m_chipi = dedxTrk->
chiPi();
704 m_chik = dedxTrk->
chiK();
705 m_chip = dedxTrk->
chiP();
708 m_probPH = dedxTrk->
probPH();
709 m_normPH = dedxTrk->
normPH();
710
711 m_tuple3->write();
712 }
713 }
714
715
716
717
718 if ( m_checkTof )
719 {
720 int m_endcap_cout( 0 ), m_layer1_cout( 0 ), m_layer2_cout( 0 );
721 for ( int i = 0; i < nGood; i++ )
722 {
724 if ( !( *itTrk )->isTofTrackValid() ) continue;
725
726 RecMdcKalTrack* mdcTrk = ( *itTrk )->mdcKalTrack();
727 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
728
729 double ptrk = mdcTrk->
p();
730
731 for ( SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
732 iter_tof != tofTrkCol.end(); iter_tof++ )
733 {
734 TofHitStatus* status = new TofHitStatus;
735 status->
setStatus( ( *iter_tof )->status() );
737 {
739 if ( status->
layer() != 0 )
continue;
740 double path = ( *iter_tof )->path();
741 double tof = ( *iter_tof )->tof();
742 double ph = ( *iter_tof )->ph();
743 double rhit = ( *iter_tof )->zrhit();
744 double qual = 0.0 + ( *iter_tof )->quality();
745 double cntr = 0.0 + ( *iter_tof )->tofID();
746 double texp[5];
747 for ( int j = 0; j < 5; j++ )
748 {
749 double gb =
xmass[j] / ptrk;
750 double beta = sqrt( 1 + gb * gb );
751 texp[j] = path * beta /
velc;
752 }
753 m_cntr_etof = cntr;
754 m_ptot_etof = ptrk;
755 m_path_etof = path;
756 m_ph_etof = ph;
757 m_rhit_etof = rhit;
758 m_qual_etof = qual;
759 m_tof_etof = tof;
760 m_te_etof = tof - texp[0];
761 m_tmu_etof = tof - texp[1];
762 m_tpi_etof = tof - texp[2];
763 m_tk_etof = tof - texp[3];
764 m_tp_etof = tof - texp[4];
765
766 m_tuple4->write();
767 }
768 else
769 {
771 if ( status->
layer() == 1 )
772 {
773 double path = ( *iter_tof )->path();
774 double tof = ( *iter_tof )->tof();
775 double ph = ( *iter_tof )->ph();
776 double rhit = ( *iter_tof )->zrhit();
777 double qual = 0.0 + ( *iter_tof )->quality();
778 double cntr = 0.0 + ( *iter_tof )->tofID();
779 double texp[5];
780 for ( int j = 0; j < 5; j++ )
781 {
782 double gb =
xmass[j] / ptrk;
783 double beta = sqrt( 1 + gb * gb );
784 texp[j] = path * beta /
velc;
785 }
786
787 m_cntr_btof1 = cntr;
788 m_ptot_btof1 = ptrk;
789 m_path_btof1 = path;
790 m_ph_btof1 = ph;
791 m_zhit_btof1 = rhit;
792 m_qual_btof1 = qual;
793 m_tof_btof1 = tof;
794 m_te_btof1 = tof - texp[0];
795 m_tmu_btof1 = tof - texp[1];
796 m_tpi_btof1 = tof - texp[2];
797 m_tk_btof1 = tof - texp[3];
798 m_tp_btof1 = tof - texp[4];
799
800 m_tuple5->write();
801 }
802
803 if ( status->
layer() == 2 )
804 {
805 double path = ( *iter_tof )->path();
806 double tof = ( *iter_tof )->tof();
807 double ph = ( *iter_tof )->ph();
808 double rhit = ( *iter_tof )->zrhit();
809 double qual = 0.0 + ( *iter_tof )->quality();
810 double cntr = 0.0 + ( *iter_tof )->tofID();
811 double texp[5];
812 for ( int j = 0; j < 5; j++ )
813 {
814 double gb =
xmass[j] / ptrk;
815 double beta = sqrt( 1 + gb * gb );
816 texp[j] = path * beta /
velc;
817 }
818
819 m_cntr_btof2 = cntr;
820 m_ptot_btof2 = ptrk;
821 m_path_btof2 = path;
822 m_ph_btof2 = ph;
823 m_zhit_btof2 = rhit;
824 m_qual_btof2 = qual;
825 m_tof_btof2 = tof;
826 m_te_btof2 = tof - texp[0];
827 m_tmu_btof2 = tof - texp[1];
828 m_tpi_btof2 = tof - texp[2];
829 m_tk_btof2 = tof - texp[3];
830 m_tp_btof2 = tof - texp[4];
831
832 m_tuple6->write();
833 }
834 }
835
836 delete status;
837 }
838 }
839 }
840
841 return StatusCode::SUCCESS;
842}
EvtRecTrackCol::iterator EvtRecTrackIterator
static long m_cout_col(0)
static long m_cout_everat(0)
static long m_cout_charge(0)
static long m_cout_nGood(0)
static long m_cout_recoil(0)
static long m_cout_mom(0)
const Hep3Vector emcPosition() const
const int emcVolumeNumber() const
static void setPidType(PidType pidType)
const HepLorentzVector p4() const
unsigned int layer() const
void setStatus(unsigned int status)
_EXTERN_ std::string TrigData