BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
CalibEventSelect Class Reference

#include <CalibEventSelect.h>

Inheritance diagram for CalibEventSelect:

Public Member Functions

 CalibEventSelect (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode initialize ()
StatusCode execute ()
StatusCode finalize ()
bool WhetherSector (double ph, double ph1, double ph2)
void readbeamEfromDb (int runNo, double &beamE)

Detailed Description

Definition at line 15 of file CalibEventSelect.h.

Constructor & Destructor Documentation

◆ CalibEventSelect()

CalibEventSelect::CalibEventSelect ( const std::string & name,
ISvcLocator * pSvcLocator )

Definition at line 114 of file CalibEventSelect.cxx.

115 : Algorithm( name, pSvcLocator ) {
116 // Declare the properties
117
118 declareProperty( "Output", m_output = false );
119 declareProperty( "Display", m_display = false );
120 declareProperty( "PrintMonitor", m_printmonitor = false );
121 declareProperty( "SelectRadBhabha", m_selectRadBhabha = false );
122 declareProperty( "SelectGGEE", m_selectGGEE = false );
123 declareProperty( "SelectGG4Pi", m_selectGG4Pi = false );
124 declareProperty( "SelectRadBhabhaT", m_selectRadBhabhaT = false );
125 declareProperty( "SelectRhoPi", m_selectRhoPi = false );
126 declareProperty( "SelectKstarK", m_selectKstarK = false );
127 declareProperty( "SelectPPPiPi", m_selectPPPiPi = false );
128 declareProperty( "SelectRecoJpsi", m_selectRecoJpsi = false );
129 declareProperty( "SelectBhabha", m_selectBhabha = false );
130 declareProperty( "SelectDimu", m_selectDimu = false );
131 declareProperty( "SelectCosmic", m_selectCosmic = false );
132 declareProperty( "SelectGenProton", m_selectGenProton = false );
133 declareProperty( "SelectPsipRhoPi", m_selectPsipRhoPi = false );
134 declareProperty( "SelectPsipKstarK", m_selectPsipKstarK = false );
135 declareProperty( "SelectPsipPPPiPi", m_selectPsipPPPiPi = false );
136 declareProperty( "SelectPsippCand", m_selectPsippCand = false );
137
138 declareProperty( "BhabhaScale", m_radbhabha_scale = 1000 );
139 declareProperty( "RadBhabhaScale", m_bhabha_scale = 1000 );
140 declareProperty( "DimuScale", m_dimu_scale = 10 );
141 declareProperty( "CosmicScale", m_cosmic_scale = 3 );
142 declareProperty( "ProtonScale", m_proton_scale = 100 );
143
144 declareProperty( "CosmicLowp", m_cosmic_lowp = 1.0 );
145
146 declareProperty( "WriteDst", m_writeDst = false );
147 declareProperty( "WriteRec", m_writeRec = false );
148 declareProperty( "Ecm", m_ecm = 3.1 );
149 declareProperty( "Vr0cut", m_vr0cut = 1.0 );
150 declareProperty( "Vz0cut", m_vz0cut = 10.0 );
151 declareProperty( "Pt0HighCut", m_pt0HighCut = 5.0 );
152 declareProperty( "Pt0LowCut", m_pt0LowCut = 0.05 );
153 declareProperty( "EnergyThreshold", m_energyThreshold = 0.05 );
154 declareProperty( "GammaPhiCut", m_gammaPhiCut = 20.0 );
155 declareProperty( "GammaThetaCut", m_gammaThetaCut = 20.0 );
156}

Referenced by CalibEventSelect().

Member Function Documentation

◆ execute()

StatusCode CalibEventSelect::execute ( )

Definition at line 484 of file CalibEventSelect.cxx.

484 {
485
486 // setFilterPassed(false);
487
488 MsgStream log( msgSvc(), name() );
489 log << MSG::INFO << "in execute()" << endmsg;
490
491 if ( m_writeDst ) m_subalg1->execute();
492 if ( m_writeRec ) m_subalg2->execute();
493
494 m_isRadBhabha = false;
495 m_isGGEE = false;
496 m_isGG4Pi = false;
497 m_isRadBhabhaT = false;
498 m_isRhoPi = false;
499 m_isKstarK = false;
500 m_isRecoJpsi = false;
501 m_isPPPiPi = false;
502 m_isBhabha = false;
503 m_isDimu = false;
504 m_isCosmic = false;
505 m_isGenProton = false;
506 m_isPsipRhoPi = false;
507 m_isPsipKstarK = false;
508 m_isPsipPPPiPi = false;
509 m_isPsippCand = false;
510
511 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
512 if ( !eventHeader )
513 {
514 cout << " eventHeader " << endl;
515 return StatusCode::FAILURE;
516 }
517
518 int run = eventHeader->runNumber();
519 int event = eventHeader->eventNumber();
520
521 // get run by run ebeam
522 if ( m_run != run )
523 {
524 m_run = run;
525 double beamE = 0;
526 readbeamEfromDb( run, beamE );
527 cout << "new run is:" << m_run << endl;
528 cout << "beamE is:" << beamE << endl;
529 if ( beamE > 0 && beamE < 3 ) m_ecm = 2 * beamE;
530 }
531
532 if ( m_display && m_events % 1000 == 0 )
533 {
534 cout << m_events << " -------- run,event: " << run << "," << event << endl;
535 cout << "m_ecm=" << m_ecm << endl;
536 }
537 m_events++;
538
539 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
540 if ( !evtRecEvent )
541 {
542 cout << " evtRecEvent " << endl;
543 return StatusCode::FAILURE;
544 }
545
546 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
547 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
548 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
549 if ( !evtRecTrkCol )
550 {
551 cout << " evtRecTrkCol " << endl;
552 return StatusCode::FAILURE;
553 }
554
555 if ( evtRecEvent->totalTracks() != evtRecTrkCol->size() ) return StatusCode::SUCCESS;
556
557 // get pi0 reconstruction
558 SmartDataPtr<EvtRecPi0Col> recPi0Col( eventSvc(), "/Event/EvtRec/EvtRecPi0Col" );
559 if ( !recPi0Col )
560 {
561 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endmsg;
562 return StatusCode::FAILURE;
563 }
564
565 int nPi0 = recPi0Col->size();
566 EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
567 if ( nPi0 == 1 )
568 {
569 double mpi0 = ( *itpi0 )->unconMass();
570 if ( fabs( mpi0 - 0.135 ) > 0.015 ) nPi0 = 0;
571 }
572
573 /*
574 int nPi0=0;
575 for(EvtRecPi0Col::iterator it=itpi0; it!= recPi0Col->end(); it++){
576 double mpi0 = (*itpi0)->unconMass();
577 if ( fabs( mpi0 - 0.135 )<= 0.015 )
578 nPi0++;
579 }
580 */
581
582 // -------- Good Charged Track Selection
583 Vint iGood;
584 iGood.clear();
585 vector<int> iCP, iCN;
586 iCP.clear();
587 iCN.clear();
588
589 Vint iCosmicGood;
590 iCosmicGood.clear();
591
592 int nCharge = 0;
593 int nCosmicCharge = 0;
594
595 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
596 {
597 // if(i>=evtRecTrkCol->size()) break;
598 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
599
600 // if(!(*itTrk)->isMdcTrackValid()) continue;
601 // RecMdcTrack *mdcTrk =(*itTrk)->mdcTrack();
602 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
603 RecMdcKalTrack* mdcTrk = ( *itTrk )->mdcKalTrack();
604 // mdcTrk->setPidType(RecMdcKalTrack::electron);
605
606 double vx0 = mdcTrk->x();
607 double vy0 = mdcTrk->y();
608 double vz0 = mdcTrk->z();
609 double vr0 = mdcTrk->r();
610 double theta0 = mdcTrk->theta();
611 double p0 = P( mdcTrk );
612 double pt0 = sqrt( pow( Px( mdcTrk ), 2 ) + pow( Py( mdcTrk ), 2 ) );
613
614 if ( m_output )
615 {
616 m_vx0 = vx0;
617 m_vy0 = vy0;
618 m_vz0 = vz0;
619 m_vr0 = vr0;
620 m_theta0 = theta0;
621 m_p0 = p0;
622 m_pt0 = pt0;
623 m_tuple1->write();
624 }
625
626 // db
627
628 Hep3Vector xorigin( 0, 0, 0 );
629 IVertexDbSvc* vtxsvc;
630 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
631 if ( vtxsvc->isVertexValid() )
632 {
633 double* dbv = vtxsvc->PrimaryVertex();
634 double* vv = vtxsvc->SigmaPrimaryVertex();
635 xorigin.setX( dbv[0] );
636 xorigin.setY( dbv[1] );
637 xorigin.setZ( dbv[2] );
638 }
639 HepVector a = mdcTrk->getZHelix();
640 HepSymMatrix Ea = mdcTrk->getZError();
641 HepPoint3D point0( 0., 0., 0. ); // the initial point for MDC recosntruction
642 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
643 VFHelix helixip( point0, a, Ea );
644 helixip.pivot( IP );
645 HepVector vecipa = helixip.a();
646 double db = fabs( vecipa[0] );
647 double dz = vecipa[3];
648
649 if ( fabs( dz ) >= m_vz0cut ) continue;
650 if ( db >= m_vr0cut ) continue;
651
652 // cosmic tracks
653 if ( p0 > m_cosmic_lowp && p0 < 20 )
654 {
655 iCosmicGood.push_back( ( *itTrk )->trackId() );
656 nCosmicCharge += mdcTrk->charge();
657 }
658
659 if ( pt0 >= m_pt0HighCut ) continue;
660 if ( pt0 <= m_pt0LowCut ) continue;
661
662 iGood.push_back( ( *itTrk )->trackId() );
663 nCharge += mdcTrk->charge();
664 if ( mdcTrk->charge() > 0 ) iCP.push_back( ( *itTrk )->trackId() );
665 else if ( mdcTrk->charge() < 0 ) iCN.push_back( ( *itTrk )->trackId() );
666 }
667 int nGood = iGood.size();
668 int nCosmicGood = iCosmicGood.size();
669
670 // -------- Good Photon Selection
671 Vint iGam;
672 iGam.clear();
673 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
674 {
675 // if(i>=evtRecTrkCol->size()) break;
676 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
677 if ( !( *itTrk )->isEmcShowerValid() ) continue;
678 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
679 double eraw = emcTrk->energy();
680 double time = emcTrk->time();
681 double costh = cos( emcTrk->theta() );
682 if ( fabs( costh ) < 0.83 && eraw < 0.025 ) continue;
683 if ( fabs( costh ) >= 0.83 && eraw < 0.05 ) continue;
684 if ( time < 0 || time > 14 ) continue;
685
686 iGam.push_back( ( *itTrk )->trackId() );
687 }
688 int nGam = iGam.size();
689
690 // -------- Assign 4-momentum to each charged track
691 Vint ipip, ipim;
692 ipip.clear();
693 ipim.clear();
694 Vp4 ppip, ppim;
695 ppip.clear();
696 ppim.clear();
697
698 // cout<<"charged track:"<<endl;
699 double echarge = 0.; // total energy of charged track
700 double ptot = 0.; // total momentum of charged track
701 double etot = 0.; // total energy in MDC and EMC
702 double eemc = 0.; // total energy in EMC
703 double pp = 0.;
704 double pm = 0.;
705 double pmax = 0.0;
706 double psec = 0.0;
707 double eccmax = 0.0;
708 double eccsec = 0.0;
709 double phimax = 0.0;
710 double phisec = 0.0;
711 double eopmax = 0.0;
712 double eopsec = 0.0;
713 Hep3Vector Pt_charge( 0, 0, 0 );
714
715 for ( int i = 0; i < nGood; i++ )
716 {
717 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
718
719 double p3 = 0;
720 // if((*itTrk)->isMdcTrackValid()) {
721 // RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
722
723 if ( ( *itTrk )->isMdcKalTrackValid() )
724 {
725 RecMdcKalTrack* mdcTrk = ( *itTrk )->mdcKalTrack();
726 mdcTrk->setPidType( RecMdcKalTrack::electron );
727
728 ptot += P( mdcTrk );
729
730 // double phi= mdcTrk->phi();
731 double phi = Phi( mdcTrk );
732
733 HepLorentzVector ptrk;
734 ptrk.setPx( Px( mdcTrk ) );
735 ptrk.setPy( Py( mdcTrk ) );
736 ptrk.setPz( Pz( mdcTrk ) );
737 p3 = fabs( ptrk.mag() );
738
739 // cout<<"p3 before compa="<<p3<<endl;
740 // cout<<"pmax before compa ="<<pmax<<endl;
741 // cout<<"psec before compa ="<<psec<<endl;
742
743 Hep3Vector ptemp( Px( mdcTrk ), Py( mdcTrk ), 0 );
744 Pt_charge += ptemp;
745
746 double ecc = 0.0;
747 if ( ( *itTrk )->isEmcShowerValid() )
748 {
749 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
750 ecc = emcTrk->energy();
751 }
752
753 if ( p3 > pmax )
754 {
755 psec = pmax;
756 eccsec = eccmax;
757 phisec = phimax;
758 pmax = p3;
759 eccmax = ecc;
760 phimax = phi;
761 }
762 else if ( p3 < pmax && p3 > psec )
763 {
764 psec = p3;
765 eccsec = ecc;
766 phisec = phi;
767 }
768
769 // cout<<"p3 after compa="<<p3<<endl;
770 // cout<<"pmax after compa ="<<pmax<<endl;
771 // cout<<"psec after compa ="<<psec<<endl;
772
773 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
774 ptrk = ptrk.boost( -0.011, 0, 0 ); // boost to cms
775
776 echarge += ptrk.e();
777 etot += ptrk.e();
778
779 if ( mdcTrk->charge() > 0 )
780 {
781 ppip.push_back( ptrk );
782 pp = ptrk.rho();
783 }
784 else
785 {
786 ppim.push_back( ptrk );
787 pm = ptrk.rho();
788 }
789 }
790
791 if ( ( *itTrk )->isEmcShowerValid() )
792 {
793
794 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
795 double eraw = emcTrk->energy();
796 double phiemc = emcTrk->phi();
797 double the = emcTrk->theta();
798
799 HepLorentzVector pemc;
800 pemc.setPx( eraw * sin( the ) * cos( phiemc ) );
801 pemc.setPy( eraw * sin( the ) * sin( phiemc ) );
802 pemc.setPz( eraw * cos( the ) );
803 pemc.setE( eraw );
804 pemc = pemc.boost( -0.011, 0, 0 ); // boost to cms
805
806 // eemc += eraw;
807 etot += pemc.e();
808 }
809
810 } // end of looping over good charged track
811
812 if ( pmax != 0 ) eopmax = eccmax / pmax;
813 if ( psec != 0 ) eopsec = eccsec / psec;
814
815 eemc = eccmax + eccsec;
816
817 /*
818 if(nGood>1){
819 cout<<"pmax="<<pmax<<endl;
820 cout<<"psec="<<psec<<endl;
821 cout<<"eopmax="<<eopmax<<endl;
822 cout<<"dphi-180="<< (fabs(phimax-phisec)-CLHEP::pi)*180/CLHEP::pi <<endl;
823 }
824 */
825
826 // -------- Assign 4-momentum to each photon
827 // cout<<"neutral:"<<endl;
828 Vp4 pGam;
829 pGam.clear();
830 double eneu = 0; // total energy of neutral track
831 double egmax = 0;
832 Hep3Vector Pt_neu( 0, 0, 0 );
833
834 for ( int i = 0; i < nGam; i++ )
835 {
836 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
837 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
838 double eraw = emcTrk->energy();
839 double phi = emcTrk->phi();
840 double the = emcTrk->theta();
841 HepLorentzVector ptrk;
842 ptrk.setPx( eraw * sin( the ) * cos( phi ) );
843 ptrk.setPy( eraw * sin( the ) * sin( phi ) );
844 ptrk.setPz( eraw * cos( the ) );
845 ptrk.setE( eraw );
846 ptrk = ptrk.boost( -0.011, 0, 0 ); // boost to cms
847 pGam.push_back( ptrk );
848 eneu += ptrk.e();
849 etot += ptrk.e();
850 eemc += eraw;
851
852 Hep3Vector ptemp( eraw * sin( the ) * cos( phi ), eraw * sin( the ) * sin( phi ), 0 );
853 Pt_neu += ptemp;
854
855 if ( ptrk.e() > egmax ) egmax = ptrk.e();
856 }
857
858 double esum = echarge + eneu;
859 Hep3Vector Pt = Pt_charge + Pt_neu;
860
861 double phid = phimax - phisec;
862 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
863
864 // -------- Select each event
865
866 // select bhabha/dimu/cosmic events, need prescale
867
868 if ( nGood == 2 && nCharge == 0 && ( m_selectBhabha || m_selectDimu ) )
869 {
870
871 // bhabha
872 if ( m_events % m_bhabha_scale == 0 && m_selectBhabha && echarge / m_ecm > 0.8 &&
873 eopmax > 0.8 && eopsec > 0.8 )
874 {
875 m_isBhabha = true;
876 m_bhabhaNumber++;
877 }
878
879 // dimu
880 if ( iCP.size() == 1 && iCN.size() == 1 && m_events % m_dimu_scale == 0 && m_selectDimu &&
881 eemc / m_ecm < 0.3 )
882 {
883
884 EvtRecTrackIterator itp = evtRecTrkCol->begin() + iCP[0];
885 EvtRecTrackIterator itn = evtRecTrkCol->begin() + iCN[0];
886
887 /*
888 double time1=-99,depth1=-99;
889 double time2=-99,depth2=-99;
890 if( (*itp)->isTofTrackValid() ){
891 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
892 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
893
894 for(;iter_tof!=tofTrkCol.end();iter_tof++){
895 TofHitStatus* status =new TofHitStatus;
896 status->setStatus( (*iter_tof)->status() );
897 if(status->is_cluster()){
898 time1=(*iter_tof)->tof();
899 }
900 delete status;
901 }
902 }
903
904 if( (*itp)->isMucTrackValid() ){
905 RecMucTrack* mucTrk=(*itp)->mucTrack();
906 depth1=mucTrk->depth();
907 }
908
909 if( (*itn)->isTofTrackValid() ){
910 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
911 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
912
913 for(;iter_tof!=tofTrkCol.end();iter_tof++){
914 TofHitStatus* status =new TofHitStatus;
915 status->setStatus( (*iter_tof)->status() );
916 if(status->is_cluster()){
917 time2=(*iter_tof)->tof();
918 }
919 delete status;
920 }
921 }
922
923 if( (*itn)->isMucTrackValid() ){
924 RecMucTrack* mucTrk=(*itn)->mucTrack();
925 depth2=mucTrk->depth();
926 }
927
928
929 //dimu
930 //if( m_selectDimu && time1!=-99 && time2!=-99 && fabs(time1-time2)<5 && depth1>5 &&
931 depth2>5){ //tight, no endcap
932 // if( eopmax<0.3 && eopsec<0.3 && fabs( fabs(phid)-CLHEP::pi)*180.0/CLHEP::pi<6 &&
933 fabs(psec)>m_ecm/4.0 ){ //tight, no rad dimu
934 // if( eemc<0.3*m_ecm && (fabs(pmax)>0.45*m_ecm || fabs(psec)>0.45*m_ecm) ){
935 if( ((fabs(pmax)/0.5/m_ecm>0.15 && fabs(pmax)/0.5/m_ecm<.75) ||
936 (fabs(psec)/0.5/m_ecm>0.15 && fabs(psec)/0.5/m_ecm<.75)) && (eopmax<0.4 || eopsec<0.4)
937 && (depth1>=3 || depth2>=3)){
938 m_isDimu=true;
939 m_dimuNumber++;
940 }
941 */
942
943 double time1 = -99, depth1 = -99;
944 double time2 = -99, depth2 = -99;
945 double p1 = -99, p2 = -99;
946 double emc1 = 0, emc2 = 0;
947 if ( ( *itp )->isTofTrackValid() )
948 {
949 SmartRefVector<RecTofTrack> tofTrkCol = ( *itp )->tofTrack();
950 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
951
952 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
953 {
954 TofHitStatus* status = new TofHitStatus;
955 status->setStatus( ( *iter_tof )->status() );
956 if ( status->is_cluster() ) { time1 = ( *iter_tof )->tof(); }
957 delete status;
958 }
959 }
960
961 if ( ( *itp )->isMucTrackValid() )
962 {
963 RecMucTrack* mucTrk = ( *itp )->mucTrack();
964 depth1 = mucTrk->depth();
965 }
966
967 if ( ( *itp )->isMdcKalTrackValid() )
968 {
969 RecMdcKalTrack* mdcTrk = ( *itp )->mdcKalTrack();
970 mdcTrk->setPidType( RecMdcKalTrack::muon );
971 p1 = P( mdcTrk );
972 }
973
974 if ( ( *itp )->isEmcShowerValid() )
975 {
976 RecEmcShower* emcTrk = ( *itp )->emcShower();
977 emc1 = emcTrk->energy();
978 }
979
980 if ( ( *itn )->isTofTrackValid() )
981 {
982 SmartRefVector<RecTofTrack> tofTrkCol = ( *itn )->tofTrack();
983 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
984
985 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
986 {
987 TofHitStatus* status = new TofHitStatus;
988 status->setStatus( ( *iter_tof )->status() );
989 if ( status->is_cluster() ) { time2 = ( *iter_tof )->tof(); }
990 delete status;
991 }
992 }
993
994 if ( ( *itn )->isMucTrackValid() )
995 {
996 RecMucTrack* mucTrk = ( *itn )->mucTrack();
997 depth2 = mucTrk->depth();
998 }
999
1000 if ( ( *itn )->isMdcKalTrackValid() )
1001 {
1002 RecMdcKalTrack* mdcTrk = ( *itn )->mdcKalTrack();
1003 mdcTrk->setPidType( RecMdcKalTrack::muon );
1004 p2 = P( mdcTrk );
1005 }
1006
1007 if ( ( *itn )->isEmcShowerValid() )
1008 {
1009 RecEmcShower* emcTrk = ( *itn )->emcShower();
1010 emc2 = emcTrk->energy();
1011 }
1012
1013 double ebeam = m_ecm / 2.0;
1014 if ( fabs( time1 - time2 ) < 5 &&
1015 ( ( p1 / ebeam > 0.15 && p1 / ebeam < 0.75 ) ||
1016 ( p2 / ebeam > 0.15 && p2 / ebeam < 0.75 ) ) &&
1017 ( emc1 / p1 < 0.4 || emc2 / p2 < 0.4 ) && ( depth1 > 3 || depth2 > 3 ) &&
1018 emc1 > 0.15 && emc1 < 0.3 && emc2 > 0.15 && emc2 < 0.3 && emc1 / p1 < 0.8 &&
1019 emc2 / p2 < 0.8 &&
1020 ( ( depth1 > 80 * p1 - 60 || depth1 > 40 ) ||
1021 ( depth2 > 80 * p2 - 60 || depth2 > 40 ) ) )
1022 {
1023 m_isDimu = true;
1024 m_dimuNumber++;
1025 }
1026
1027 } // end of dimu
1028
1029 } // end of bhabha, dimu
1030
1031 if ( m_selectCosmic && nCosmicGood == 2 && eemc / m_ecm < 0.3 )
1032 {
1033
1034 EvtRecTrackIterator itp = evtRecTrkCol->begin() + iCosmicGood[0];
1035 EvtRecTrackIterator itn = evtRecTrkCol->begin() + iCosmicGood[1];
1036
1037 double time1 = -99, depth1 = -99;
1038 double time2 = -99, depth2 = -99;
1039 if ( ( *itp )->isTofTrackValid() )
1040 {
1041 SmartRefVector<RecTofTrack> tofTrkCol = ( *itp )->tofTrack();
1042 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
1043
1044 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
1045 {
1046 TofHitStatus* status = new TofHitStatus;
1047 status->setStatus( ( *iter_tof )->status() );
1048 if ( status->is_cluster() ) { time1 = ( *iter_tof )->tof(); }
1049 delete status;
1050 }
1051 }
1052
1053 if ( ( *itp )->isMucTrackValid() )
1054 {
1055 RecMucTrack* mucTrk = ( *itp )->mucTrack();
1056 depth1 = mucTrk->depth();
1057 }
1058
1059 if ( ( *itn )->isTofTrackValid() )
1060 {
1061 SmartRefVector<RecTofTrack> tofTrkCol = ( *itn )->tofTrack();
1062 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
1063
1064 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
1065 {
1066 TofHitStatus* status = new TofHitStatus;
1067 status->setStatus( ( *iter_tof )->status() );
1068 if ( status->is_cluster() ) { time2 = ( *iter_tof )->tof(); }
1069 delete status;
1070 }
1071 }
1072
1073 if ( ( *itn )->isMucTrackValid() )
1074 {
1075 RecMucTrack* mucTrk = ( *itn )->mucTrack();
1076 depth2 = mucTrk->depth();
1077 }
1078
1079 // cosmic
1080 // if( m_selectCosmic && time1!=-99 && time2!=-99 && fabs(time1-time2)>5 && depth1>5 &&
1081 // depth2>5){
1082 if ( m_selectCosmic && time1 != -99 && time2 != -99 && fabs( time1 - time2 ) > 5 )
1083 {
1084 m_isCosmic = true;
1085 m_cosmicNumber++;
1086 }
1087
1088 } // end of cosmic
1089
1090 // select generic protons
1091
1092 if ( m_events % m_proton_scale == 0 )
1093 {
1094
1095 bool protoncand = false;
1096 int ncharge = 0;
1097 int nproton = 0;
1098 for ( int i = 0; i < nGood; i++ )
1099 {
1100
1101 EvtRecTrackIterator iter = evtRecTrkCol->begin() + iGood[i];
1102 RecMdcKalTrack* mdcTrk = ( *iter )->mdcKalTrack();
1103 mdcTrk->setPidType( RecMdcKalTrack::proton );
1104
1105 double pp = P( mdcTrk );
1106 double dedx = -99;
1107 double chiP = -99;
1108
1109 if ( ( *iter )->isMdcDedxValid() )
1110 {
1111 RecMdcDedx* dedxTrk = ( *iter )->mdcDedx();
1112 dedx = dedxTrk->normPH();
1113 chiP = dedxTrk->chiP();
1114 }
1115
1116 if ( fabs( pp ) < 0.6 && dedx > 1.2 && fabs( chiP ) < 6 )
1117 {
1118 ncharge += mdcTrk->charge();
1119 nproton++;
1120 // protoncand=true;
1121 // break;
1122 }
1123 }
1124 if ( ( nproton == 1 && ncharge == -1 ) || ( nproton >= 2 && ncharge <= nproton - 2 ) )
1125 protoncand = true;
1126 if ( protoncand )
1127 {
1128 m_isGenProton = true;
1129 m_genProtonNumber++;
1130 }
1131
1132 } // end of generic proton
1133
1134 // fill monitoring histograms for rad bhabha
1135
1136 if ( m_printmonitor )
1137 {
1138 h_nGood->fill( nGood );
1139 h_nCharge->fill( nCharge );
1140 h_pmaxobeam->fill( pmax / ( m_ecm / 2.0 ) );
1141 h_eopmax->fill( eopmax );
1142 h_eopsec->fill( eopsec );
1143 h_deltap->fill( pmax - psec );
1144 h_esumoecm->fill( esum / m_ecm );
1145 h_egmax->fill( egmax );
1146 h_deltaphi->fill( phid );
1147 h_Pt->fill( Pt.mag() );
1148 }
1149
1150 // select radbhabha
1151 if ( nGood > 1 && pmax / ( m_ecm / 2.0 ) > 0.4 && eopmax > 0.5 && esum / m_ecm > 0.75 &&
1152 egmax > 0.08 && fabs( fabs( phid ) - CLHEP::pi ) * 180.0 / CLHEP::pi > 2.86 )
1153 {
1154 m_isRadBhabha = true;
1155 m_radBhabhaNumber++;
1156 }
1157 // select radbhabha tight
1158 if ( m_isRadBhabha )
1159 {
1160 // prescale high momentum tracks
1161 if ( nGood == 2 && nCharge == 0 && eemc / m_ecm >= 0.7 && eopmax >= 0.85 &&
1162 eopmax <= 1.15 && eopsec >= 0.85 && eopsec <= 1.15 )
1163 {
1164
1165 if ( fabs( fabs( pmax ) - m_ecm / 2.0 ) < 0.1 &&
1166 fabs( fabs( psec ) - m_ecm / 2.0 ) < 0.1 )
1167 {
1168 if ( m_events % 1000 == 0 )
1169 {
1170 m_isRadBhabhaT = true;
1171 m_radBhabhaTNumber++;
1172 }
1173 }
1174 else
1175 {
1176 m_isRadBhabhaT = true;
1177 m_radBhabhaTNumber++;
1178 }
1179 }
1180 }
1181
1182 // select ggee events, (exclude radee tight)
1183 // if(!m_isRadBhabhaT && nGood==2 && nCharge==0 && eopmax>=0.85 && eopsec>=0.85 &&
1184 // eemc/m_ecm<=0.8 && Pt.mag()<=0.2)
1185 if ( !m_isRadBhabhaT && nGood == 2 && nCharge == 0 && eopmax >= 0.85 && eopmax <= 1.15 &&
1186 eopsec >= 0.85 && eopsec <= 1.15 && eemc / m_ecm <= 0.8 && Pt.mag() <= 0.2 )
1187 {
1188 m_isGGEE = true;
1189 m_GGEENumber++;
1190 }
1191
1192 // select gg4pi events
1193 if ( m_selectGG4Pi && nGood == 4 && nCharge == 0 && pmax / ( m_ecm / 2.0 ) > 0.04 &&
1194 pmax / ( m_ecm / 2.0 ) < 0.9 && esum / m_ecm > 0.04 && esum / m_ecm < 0.75 &&
1195 Pt.mag() <= 0.2 )
1196 {
1197 m_isGG4Pi = true;
1198 m_GG4PiNumber++;
1199 }
1200
1201 // select rhopi/kstark events
1202 if ( ( m_selectRhoPi || m_selectKstarK ) && nGood == 2 && nCharge == 0 && nPi0 == 1 )
1203 {
1204
1205 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iGood[0];
1206 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iGood[1];
1207
1208 if ( ( *itone )->isMdcKalTrackValid() && ( *ittwo )->isMdcKalTrackValid() )
1209 {
1210
1211 RecMdcKalTrack* trk1 = ( *itone )->mdcKalTrack();
1212 RecMdcKalTrack* trk2 = ( *ittwo )->mdcKalTrack();
1213
1214 HepLorentzVector p4trk1;
1215 p4trk1.setPx( Px( trk1 ) );
1216 p4trk1.setPy( Py( trk1 ) );
1217 p4trk1.setPz( Pz( trk1 ) );
1218 p4trk1.setE( sqrt( pow( P( trk1 ), 2 ) + mkaon * mkaon ) );
1219
1220 HepLorentzVector p4trk2;
1221 p4trk2.setPx( Px( trk2 ) );
1222 p4trk2.setPy( Py( trk2 ) );
1223 p4trk2.setPz( Pz( trk2 ) );
1224 p4trk2.setE( sqrt( pow( P( trk2 ), 2 ) + mkaon * mkaon ) );
1225
1226 HepLorentzVector p4trk3;
1227 p4trk3.setPx( Px( trk1 ) );
1228 p4trk3.setPy( Py( trk1 ) );
1229 p4trk3.setPz( Pz( trk1 ) );
1230 p4trk3.setE( sqrt( pow( P( trk1 ), 2 ) + mpi * mpi ) );
1231
1232 HepLorentzVector p4trk4;
1233 p4trk4.setPx( Px( trk2 ) );
1234 p4trk4.setPy( Py( trk2 ) );
1235 p4trk4.setPz( Pz( trk2 ) );
1236 p4trk4.setE( sqrt( pow( P( trk2 ), 2 ) + mpi * mpi ) );
1237
1238 // EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
1239 itpi0 = recPi0Col->begin();
1240 HepLorentzVector p4gam1 = ( *itpi0 )->hiPfit();
1241 HepLorentzVector p4gam2 = ( *itpi0 )->loPfit();
1242 HepLorentzVector p4pi0 = p4gam1 + p4gam2;
1243
1244 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0; // kstark
1245 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0; // rhopi
1246
1247 double rhopimass = p4total2.m();
1248 double kstarkmass = p4total.m();
1249 if ( ( kstarkmass > 3.0 && kstarkmass < 3.15 ) ||
1250 ( rhopimass > 3.0 && rhopimass < 3.15 ) )
1251 {
1252
1253 // tight cuts
1254
1255 // remove bhabha background
1256 double eop1 = 0.0, eop2 = 0.0;
1257 if ( ( *itone )->isEmcShowerValid() )
1258 {
1259 RecEmcShower* emcTrk = ( *itone )->emcShower();
1260 double etrk = emcTrk->energy();
1261 // cout<<"trk1 p="<<P(trk1)<<endl;
1262 if ( P( trk1 ) != 0 )
1263 {
1264 eop1 = etrk / P( trk1 );
1265 // if( fabs(eop1)> 0.8 ) return StatusCode::SUCCESS;
1266 }
1267 }
1268
1269 if ( ( *ittwo )->isEmcShowerValid() )
1270 {
1271 RecEmcShower* emcTrk = ( *ittwo )->emcShower();
1272 double etrk = emcTrk->energy();
1273 // cout<<"trk2 p="<<P(trk2)<<endl;
1274 if ( P( trk2 ) != 0 )
1275 {
1276 eop2 = etrk / P( trk2 );
1277 // if( fabs(eop2)> 0.8 ) return StatusCode::SUCCESS;
1278 }
1279 }
1280
1281 if ( eop1 < 0.8 && eop2 < 0.8 )
1282 {
1283
1284 if ( rhopimass > 3.0 && rhopimass < 3.15 )
1285 {
1286 // kinematic fit
1287 // HepLorentzVector ecms(0.034,0,0,3.097);
1288 HepLorentzVector ecms( 0.034, 0, 0, m_ecm );
1289
1290 WTrackParameter wvpiTrk1, wvpiTrk2;
1291 wvpiTrk1 = WTrackParameter( mpi, trk1->getZHelix(), trk1->getZError() );
1292 wvpiTrk2 = WTrackParameter( mpi, trk2->getZHelix(), trk2->getZError() );
1293
1294 const EvtRecTrack* gTrk1 = ( *itpi0 )->hiEnGamma();
1295 const EvtRecTrack* gTrk2 = ( *itpi0 )->loEnGamma();
1296 RecEmcShower* gShr1 = const_cast<EvtRecTrack*>( gTrk1 )->emcShower();
1297 RecEmcShower* gShr2 = const_cast<EvtRecTrack*>( gTrk2 )->emcShower();
1298
1299 KinematicFit* kmfit = KinematicFit::instance();
1300 kmfit->init();
1301 kmfit->AddTrack( 0, wvpiTrk1 );
1302 kmfit->AddTrack( 1, wvpiTrk2 );
1303 kmfit->AddTrack( 2, 0.0, gShr1 );
1304 kmfit->AddTrack( 3, 0.0, gShr2 );
1305 kmfit->AddFourMomentum( 0, ecms );
1306
1307 bool oksq1 = kmfit->Fit();
1308 double chi1 = kmfit->chisq();
1309
1310 //
1311 if ( oksq1 && chi1 <= 60 )
1312 {
1313 m_isRhoPi = true;
1314 m_rhoPiNumber++;
1315 }
1316 } // end of selecting rhopi
1317
1318 if ( kstarkmass > 3.0 && kstarkmass < 3.15 )
1319 {
1320
1321 // kstark resonances
1322 double mkstarp = 0, mkstarm = 0;
1323 if ( trk1->charge() > 0 )
1324 {
1325 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
1326 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
1327 mkstarp = p4kstarp.m();
1328 mkstarm = p4kstarm.m();
1329 }
1330 else
1331 {
1332 HepLorentzVector p4kstarm = p4trk1 + p4pi0;
1333 HepLorentzVector p4kstarp = p4trk2 + p4pi0;
1334 mkstarp = p4kstarp.m();
1335 mkstarm = p4kstarm.m();
1336 }
1337
1338 if ( ( fabs( mkstarp - 0.89166 ) < 0.1 && fabs( mkstarm - 0.89166 ) > 0.1 ) ||
1339 ( fabs( mkstarm - 0.89166 ) < 0.1 && fabs( mkstarp - 0.89166 ) > 0.1 ) )
1340 {
1341 // kinematic fit
1342 // HepLorentzVector ecms(0.034,0,0,3.097);
1343 HepLorentzVector ecms( 0.034, 0, 0, m_ecm );
1344
1345 WTrackParameter wvpiTrk1, wvpiTrk2;
1346 wvpiTrk1 = WTrackParameter( mkaon, trk1->getZHelix(), trk1->getZError() );
1347 wvpiTrk2 = WTrackParameter( mkaon, trk2->getZHelix(), trk2->getZError() );
1348
1349 const EvtRecTrack* gTrk1 = ( *itpi0 )->hiEnGamma();
1350 const EvtRecTrack* gTrk2 = ( *itpi0 )->loEnGamma();
1351 RecEmcShower* gShr1 = const_cast<EvtRecTrack*>( gTrk1 )->emcShower();
1352 RecEmcShower* gShr2 = const_cast<EvtRecTrack*>( gTrk2 )->emcShower();
1353
1354 KinematicFit* kmfit = KinematicFit::instance();
1355 kmfit->init();
1356 kmfit->AddTrack( 0, wvpiTrk1 );
1357 kmfit->AddTrack( 1, wvpiTrk2 );
1358 kmfit->AddTrack( 2, 0.0, gShr1 );
1359 kmfit->AddTrack( 3, 0.0, gShr2 );
1360 kmfit->AddFourMomentum( 0, ecms );
1361
1362 bool oksq1 = kmfit->Fit();
1363 double chi1 = kmfit->chisq();
1364 //
1365
1366 if ( oksq1 && chi1 <= 60 )
1367 {
1368 m_isKstarK = true;
1369 m_kstarKNumber++;
1370 }
1371 }
1372
1373 } // end of selecting kstark
1374
1375 } // end of non di-electron
1376
1377 // end of tight cuts
1378 }
1379 }
1380
1381 } // end of selecting rhopi/kstark events
1382
1383 // select ppPiPi events
1384 if ( m_selectPPPiPi && nGood == 4 && nCharge == 0 && nPi0 == 0 )
1385 {
1386
1387 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iCP[0];
1388 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iCP[1];
1389 EvtRecTrackIterator itthr = evtRecTrkCol->begin() + iCN[0];
1390 EvtRecTrackIterator itfor = evtRecTrkCol->begin() + iCN[1];
1391 RecMdcKalTrack* trk1 = ( *itone )->mdcKalTrack();
1392 RecMdcKalTrack* trk2 = ( *ittwo )->mdcKalTrack();
1393 RecMdcKalTrack* trk3 = ( *itthr )->mdcKalTrack();
1394 RecMdcKalTrack* trk4 = ( *itfor )->mdcKalTrack();
1395
1396 HepLorentzVector p4trkpp;
1397 HepLorentzVector p4trkpm;
1398 HepLorentzVector p4trkpip;
1399 HepLorentzVector p4trkpim;
1400
1401 // hypothesis 1
1402
1403 trk1->setPidType( RecMdcKalTrack::proton );
1404 p4trkpp.setPx( Px( trk1 ) );
1405 p4trkpp.setPy( Py( trk1 ) );
1406 p4trkpp.setPz( Pz( trk1 ) );
1407 p4trkpp.setE( sqrt( pow( P( trk1 ), 2 ) + mproton * mproton ) );
1408
1409 trk3->setPidType( RecMdcKalTrack::proton );
1410 p4trkpm.setPx( Px( trk3 ) );
1411 p4trkpm.setPy( Py( trk3 ) );
1412 p4trkpm.setPz( Pz( trk3 ) );
1413 p4trkpm.setE( sqrt( pow( P( trk3 ), 2 ) + mproton * mproton ) );
1414
1415 trk2->setPidType( RecMdcKalTrack::pion );
1416 p4trkpip.setPx( Px( trk2 ) );
1417 p4trkpip.setPy( Py( trk2 ) );
1418 p4trkpip.setPz( Pz( trk2 ) );
1419 p4trkpip.setE( sqrt( pow( P( trk2 ), 2 ) + mpi * mpi ) );
1420
1421 trk4->setPidType( RecMdcKalTrack::pion );
1422 p4trkpim.setPx( Px( trk4 ) );
1423 p4trkpim.setPy( Py( trk4 ) );
1424 p4trkpim.setPz( Pz( trk4 ) );
1425 p4trkpim.setE( sqrt( pow( P( trk4 ), 2 ) + mpi * mpi ) );
1426
1427 double jpsimass1 = ( p4trkpp + p4trkpm + p4trkpip + p4trkpim ).m();
1428
1429 // hypothesis 2
1430
1431 trk1->setPidType( RecMdcKalTrack::proton );
1432 p4trkpp.setPx( Px( trk1 ) );
1433 p4trkpp.setPy( Py( trk1 ) );
1434 p4trkpp.setPz( Pz( trk1 ) );
1435 p4trkpp.setE( sqrt( pow( P( trk1 ), 2 ) + mproton * mproton ) );
1436
1437 trk4->setPidType( RecMdcKalTrack::proton );
1438 p4trkpm.setPx( Px( trk4 ) );
1439 p4trkpm.setPy( Py( trk4 ) );
1440 p4trkpm.setPz( Pz( trk4 ) );
1441 p4trkpm.setE( sqrt( pow( P( trk4 ), 2 ) + mproton * mproton ) );
1442
1443 trk2->setPidType( RecMdcKalTrack::pion );
1444 p4trkpip.setPx( Px( trk2 ) );
1445 p4trkpip.setPy( Py( trk2 ) );
1446 p4trkpip.setPz( Pz( trk2 ) );
1447 p4trkpip.setE( sqrt( pow( P( trk2 ), 2 ) + mpi * mpi ) );
1448
1449 trk3->setPidType( RecMdcKalTrack::pion );
1450 p4trkpim.setPx( Px( trk3 ) );
1451 p4trkpim.setPy( Py( trk3 ) );
1452 p4trkpim.setPz( Pz( trk3 ) );
1453 p4trkpim.setE( sqrt( pow( P( trk3 ), 2 ) + mpi * mpi ) );
1454
1455 double jpsimass2 = ( p4trkpp + p4trkpm + p4trkpip + p4trkpim ).m();
1456
1457 // hypothesis 3
1458
1459 trk2->setPidType( RecMdcKalTrack::proton );
1460 p4trkpp.setPx( Px( trk2 ) );
1461 p4trkpp.setPy( Py( trk2 ) );
1462 p4trkpp.setPz( Pz( trk2 ) );
1463 p4trkpp.setE( sqrt( pow( P( trk2 ), 2 ) + mproton * mproton ) );
1464
1465 trk3->setPidType( RecMdcKalTrack::proton );
1466 p4trkpm.setPx( Px( trk3 ) );
1467 p4trkpm.setPy( Py( trk3 ) );
1468 p4trkpm.setPz( Pz( trk3 ) );
1469 p4trkpm.setE( sqrt( pow( P( trk3 ), 2 ) + mproton * mproton ) );
1470
1471 trk1->setPidType( RecMdcKalTrack::pion );
1472 p4trkpip.setPx( Px( trk1 ) );
1473 p4trkpip.setPy( Py( trk1 ) );
1474 p4trkpip.setPz( Pz( trk1 ) );
1475 p4trkpip.setE( sqrt( pow( P( trk1 ), 2 ) + mpi * mpi ) );
1476
1477 trk4->setPidType( RecMdcKalTrack::pion );
1478 p4trkpim.setPx( Px( trk4 ) );
1479 p4trkpim.setPy( Py( trk4 ) );
1480 p4trkpim.setPz( Pz( trk4 ) );
1481 p4trkpim.setE( sqrt( pow( P( trk4 ), 2 ) + mpi * mpi ) );
1482
1483 double jpsimass3 = ( p4trkpp + p4trkpm + p4trkpip + p4trkpim ).m();
1484
1485 // hypothesis 4
1486
1487 trk2->setPidType( RecMdcKalTrack::proton );
1488 p4trkpp.setPx( Px( trk2 ) );
1489 p4trkpp.setPy( Py( trk2 ) );
1490 p4trkpp.setPz( Pz( trk2 ) );
1491 p4trkpp.setE( sqrt( pow( P( trk2 ), 2 ) + mproton * mproton ) );
1492
1493 trk4->setPidType( RecMdcKalTrack::proton );
1494 p4trkpm.setPx( Px( trk4 ) );
1495 p4trkpm.setPy( Py( trk4 ) );
1496 p4trkpm.setPz( Pz( trk4 ) );
1497 p4trkpm.setE( sqrt( pow( P( trk4 ), 2 ) + mproton * mproton ) );
1498
1499 trk1->setPidType( RecMdcKalTrack::pion );
1500 p4trkpip.setPx( Px( trk1 ) );
1501 p4trkpip.setPy( Py( trk1 ) );
1502 p4trkpip.setPz( Pz( trk1 ) );
1503 p4trkpip.setE( sqrt( pow( P( trk1 ), 2 ) + mpi * mpi ) );
1504
1505 trk3->setPidType( RecMdcKalTrack::pion );
1506 p4trkpim.setPx( Px( trk3 ) );
1507 p4trkpim.setPy( Py( trk3 ) );
1508 p4trkpim.setPz( Pz( trk3 ) );
1509 p4trkpim.setE( sqrt( pow( P( trk3 ), 2 ) + mpi * mpi ) );
1510
1511 double jpsimass4 = ( p4trkpp + p4trkpm + p4trkpip + p4trkpim ).m();
1512
1513 if ( fabs( jpsimass1 - 3.075 ) <= 0.075 || fabs( jpsimass2 - 3.075 ) <= 0.075 ||
1514 fabs( jpsimass3 - 3.075 ) <= 0.075 || fabs( jpsimass4 - 3.075 ) <= 0.075 )
1515 {
1516
1517 // tight cuts
1518
1519 // kinematic fits
1520 double chi1, chi2, chi3, chi4;
1521 int type = 0;
1522 double chimin = -999;
1523 HepLorentzVector ecms( 0.034, 0, 0, m_ecm );
1524
1525 WTrackParameter wvpiTrk1, wvpiTrk2, wvpiTrk3, wvpiTrk4;
1526
1527 {
1528 wvpiTrk1 = WTrackParameter( mproton, trk1->getZHelix(), trk1->getZError() );
1529 wvpiTrk2 = WTrackParameter( mpi, trk2->getZHelix(), trk2->getZError() );
1530 wvpiTrk3 = WTrackParameter( mproton, trk3->getZHelix(), trk3->getZError() );
1531 wvpiTrk4 = WTrackParameter( mpi, trk4->getZHelix(), trk4->getZError() );
1532
1533 KinematicFit* kmfit = KinematicFit::instance();
1534 kmfit->init();
1535 kmfit->AddTrack( 0, wvpiTrk1 );
1536 kmfit->AddTrack( 1, wvpiTrk2 );
1537 kmfit->AddTrack( 2, wvpiTrk3 );
1538 kmfit->AddTrack( 3, wvpiTrk4 );
1539 kmfit->AddFourMomentum( 0, ecms );
1540
1541 bool oksq1 = kmfit->Fit();
1542 chi1 = kmfit->chisq();
1543 if ( oksq1 )
1544 {
1545 chimin = chi1;
1546 type = 1;
1547 }
1548 //
1549 }
1550
1551 {
1552 wvpiTrk1 = WTrackParameter( mproton, trk1->getZHelix(), trk1->getZError() );
1553 wvpiTrk2 = WTrackParameter( mpi, trk2->getZHelix(), trk2->getZError() );
1554 wvpiTrk3 = WTrackParameter( mpi, trk3->getZHelix(), trk3->getZError() );
1555 wvpiTrk4 = WTrackParameter( mproton, trk4->getZHelix(), trk4->getZError() );
1556
1557 KinematicFit* kmfit = KinematicFit::instance();
1558 kmfit->init();
1559 kmfit->AddTrack( 0, wvpiTrk1 );
1560 kmfit->AddTrack( 1, wvpiTrk2 );
1561 kmfit->AddTrack( 2, wvpiTrk3 );
1562 kmfit->AddTrack( 3, wvpiTrk4 );
1563 kmfit->AddFourMomentum( 0, ecms );
1564
1565 bool oksq1 = kmfit->Fit();
1566 chi2 = kmfit->chisq();
1567 if ( oksq1 )
1568 {
1569 if ( type == 0 )
1570 {
1571 chimin = chi2;
1572 type = 2;
1573 }
1574 else if ( chi2 < chimin )
1575 {
1576 chimin = chi2;
1577 type = 2;
1578 }
1579 }
1580 //
1581 }
1582
1583 {
1584 wvpiTrk1 = WTrackParameter( mpi, trk1->getZHelix(), trk1->getZError() );
1585 wvpiTrk2 = WTrackParameter( mproton, trk2->getZHelix(), trk2->getZError() );
1586 wvpiTrk3 = WTrackParameter( mpi, trk3->getZHelix(), trk3->getZError() );
1587 wvpiTrk4 = WTrackParameter( mproton, trk4->getZHelix(), trk4->getZError() );
1588
1589 KinematicFit* kmfit = KinematicFit::instance();
1590 kmfit->init();
1591 kmfit->AddTrack( 0, wvpiTrk1 );
1592 kmfit->AddTrack( 1, wvpiTrk2 );
1593 kmfit->AddTrack( 2, wvpiTrk3 );
1594 kmfit->AddTrack( 3, wvpiTrk4 );
1595
1596 kmfit->AddFourMomentum( 0, ecms );
1597
1598 bool oksq1 = kmfit->Fit();
1599 chi3 = kmfit->chisq();
1600 if ( oksq1 )
1601 {
1602 if ( type == 0 )
1603 {
1604 chimin = chi3;
1605 type = 3;
1606 }
1607 else if ( chi3 < chimin )
1608 {
1609 chimin = chi3;
1610 type = 3;
1611 }
1612 }
1613
1614 //
1615 }
1616
1617 {
1618 wvpiTrk1 = WTrackParameter( mpi, trk1->getZHelix(), trk1->getZError() );
1619 wvpiTrk2 = WTrackParameter( mproton, trk2->getZHelix(), trk2->getZError() );
1620 wvpiTrk3 = WTrackParameter( mproton, trk3->getZHelix(), trk3->getZError() );
1621 wvpiTrk4 = WTrackParameter( mpi, trk4->getZHelix(), trk4->getZError() );
1622
1623 KinematicFit* kmfit = KinematicFit::instance();
1624 kmfit->init();
1625 kmfit->AddTrack( 0, wvpiTrk1 );
1626 kmfit->AddTrack( 1, wvpiTrk2 );
1627 kmfit->AddTrack( 2, wvpiTrk3 );
1628 kmfit->AddTrack( 3, wvpiTrk4 );
1629
1630 kmfit->AddFourMomentum( 0, ecms );
1631
1632 bool oksq1 = kmfit->Fit();
1633 chi4 = kmfit->chisq();
1634
1635 if ( oksq1 )
1636 {
1637 if ( type == 0 )
1638 {
1639 chimin = chi4;
1640 type = 4;
1641 }
1642 else if ( chi4 < chimin )
1643 {
1644 chimin = chi4;
1645 type = 4;
1646 }
1647 }
1648
1649 //
1650 }
1651
1652 if ( type != 0 && chimin < 100 )
1653 {
1654 m_isPPPiPi = true;
1655 m_ppPiPiNumber++;
1656 }
1657
1658 // end of tight cuts
1659
1660 } // end of loose cuts
1661
1662 } // end of selecting pppipi
1663
1664 // select recoil events
1665 // if( m_selectRecoJpsi && (nGood==4 || nGood==6) && nCharge==0 ){
1666 if ( ( m_selectPsipRhoPi || m_selectPsipKstarK || m_selectPsipPPPiPi ) &&
1667 ( nGood == 4 || nGood == 6 ) && nCharge == 0 )
1668 {
1669 EvtRecTrackIterator pione, pitwo;
1670 RecMdcKalTrack* pitrk1;
1671 RecMdcKalTrack* pitrk2;
1672
1673 double bestmass = 1.0;
1674 int pindex, nindex;
1675 vector<int> iJPsiP, iJPsiN;
1676 for ( int ip = 0; ip < iCP.size(); ip++ )
1677 {
1678 for ( int in = 0; in < iCN.size(); in++ )
1679 {
1680 pione = evtRecTrkCol->begin() + iCP[ip];
1681 pitwo = evtRecTrkCol->begin() + iCN[in];
1682 pitrk1 = ( *pione )->mdcKalTrack();
1683 pitrk2 = ( *pitwo )->mdcKalTrack();
1684 Hep3Vector p1( Px( pitrk1 ), Py( pitrk1 ), Pz( pitrk1 ) );
1685 Hep3Vector p2( Px( pitrk2 ), Py( pitrk2 ), Pz( pitrk2 ) );
1686 double E1 = sqrt( pow( P( pitrk1 ), 2 ) + mpi * mpi );
1687 double E2 = sqrt( pow( P( pitrk2 ), 2 ) + mpi * mpi );
1688 double recomass = sqrt( pow( 3.686 - E1 - E2, 2 ) - ( -( p1 + p2 ) ).mag2() );
1689 // if(fabs(recomass-3.686)< fabs(bestmass-3.686)){
1690 if ( fabs( recomass - 3.096 ) < fabs( bestmass - 3.096 ) )
1691 {
1692 bestmass = recomass;
1693 pindex = ip;
1694 nindex = in;
1695 }
1696 }
1697 }
1698
1699 // soft pions
1700 pione = evtRecTrkCol->begin() + iCP[pindex];
1701 pitwo = evtRecTrkCol->begin() + iCN[nindex];
1702 pitrk1 = ( *pione )->mdcKalTrack();
1703 pitrk2 = ( *pitwo )->mdcKalTrack();
1704
1705 // tracks from jpsi
1706 double jpsicharge = 0;
1707 for ( int ip = 0; ip < iCP.size(); ip++ )
1708 {
1709 if ( ip != pindex )
1710 {
1711 iJPsiP.push_back( iCP[ip] );
1712 EvtRecTrackIterator itertmp = evtRecTrkCol->begin() + iCP[ip];
1713 RecMdcKalTrack* trktmp = ( *itertmp )->mdcKalTrack();
1714 jpsicharge += trktmp->charge();
1715 }
1716 }
1717
1718 for ( int in = 0; in < iCN.size(); in++ )
1719 {
1720 if ( in != nindex )
1721 {
1722 iJPsiN.push_back( iCN[in] );
1723 EvtRecTrackIterator itertmp = evtRecTrkCol->begin() + iCN[in];
1724 RecMdcKalTrack* trktmp = ( *itertmp )->mdcKalTrack();
1725 jpsicharge += trktmp->charge();
1726 }
1727 }
1728
1729 HepLorentzVector ecms( 0.034, 0, 0, m_ecm );
1730
1731 // jpsi to 2 charged track and 1 pi0
1732 if ( ( m_selectPsipRhoPi || m_selectPsipKstarK ) &&
1733 ( bestmass > 3.075 && bestmass < 3.120 ) && nGood == 4 && jpsicharge == 0 &&
1734 nPi0 == 1 )
1735 {
1736
1737 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iJPsiP[0];
1738 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iJPsiN[0];
1739
1740 if ( ( *itone )->isMdcKalTrackValid() && ( *ittwo )->isMdcKalTrackValid() )
1741 {
1742
1743 RecMdcKalTrack* trk1 = ( *itone )->mdcKalTrack();
1744 RecMdcKalTrack* trk2 = ( *ittwo )->mdcKalTrack();
1745
1746 HepLorentzVector p4trk1;
1747 p4trk1.setPx( Px( trk1 ) );
1748 p4trk1.setPy( Py( trk1 ) );
1749 p4trk1.setPz( Pz( trk1 ) );
1750 p4trk1.setE( sqrt( pow( P( trk1 ), 2 ) + mkaon * mkaon ) );
1751
1752 HepLorentzVector p4trk2;
1753 p4trk2.setPx( Px( trk2 ) );
1754 p4trk2.setPy( Py( trk2 ) );
1755 p4trk2.setPz( Pz( trk2 ) );
1756 p4trk2.setE( sqrt( pow( P( trk2 ), 2 ) + mkaon * mkaon ) );
1757
1758 HepLorentzVector p4trk3;
1759 p4trk3.setPx( Px( trk1 ) );
1760 p4trk3.setPy( Py( trk1 ) );
1761 p4trk3.setPz( Pz( trk1 ) );
1762 p4trk3.setE( sqrt( pow( P( trk1 ), 2 ) + mpi * mpi ) );
1763
1764 HepLorentzVector p4trk4;
1765 p4trk4.setPx( Px( trk2 ) );
1766 p4trk4.setPy( Py( trk2 ) );
1767 p4trk4.setPz( Pz( trk2 ) );
1768 p4trk4.setE( sqrt( pow( P( trk2 ), 2 ) + mpi * mpi ) );
1769
1770 EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
1771 const EvtRecTrack* gTrk1 = ( *itpi0 )->hiEnGamma();
1772 const EvtRecTrack* gTrk2 = ( *itpi0 )->loEnGamma();
1773 RecEmcShower* gShr1 = const_cast<EvtRecTrack*>( gTrk1 )->emcShower();
1774 RecEmcShower* gShr2 = const_cast<EvtRecTrack*>( gTrk2 )->emcShower();
1775 RecEmcShower* gShr3 = const_cast<EvtRecTrack*>( gTrk1 )->emcShower();
1776 RecEmcShower* gShr4 = const_cast<EvtRecTrack*>( gTrk2 )->emcShower();
1777
1778 HepLorentzVector p4gam1 = ( *itpi0 )->hiPfit();
1779 HepLorentzVector p4gam2 = ( *itpi0 )->loPfit();
1780 HepLorentzVector p4pi0 = p4gam1 + p4gam2;
1781 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0;
1782 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0;
1783
1784 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
1785 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
1786 double mkstarp = p4kstarp.m();
1787 double mkstarm = p4kstarm.m();
1788
1789 if ( ( p4total.m() > 3.0 && p4total.m() < 3.15 ) ||
1790 ( p4total2.m() > 3.0 && p4total2.m() < 3.15 ) )
1791 {
1792 // m_isRecoJpsi = true;
1793 // m_recoJpsiNumber++;
1794
1795 // tighten cuts for rhopi and kstark
1796
1797 WTrackParameter wvpiTrk1, wvpiTrk2;
1798 wvpiTrk1 = WTrackParameter( mpi, trk1->getZHelix(), trk1->getZError() );
1799 wvpiTrk2 = WTrackParameter( mpi, trk2->getZHelix(), trk2->getZError() );
1800
1801 WTrackParameter wvkTrk1, wvkTrk2;
1802 wvkTrk1 = WTrackParameter( mkaon, trk1->getZHelixK(), trk1->getZErrorK() ); // kaon
1803 wvkTrk2 = WTrackParameter( mkaon, trk2->getZHelixK(), trk2->getZErrorK() ); // kaon
1804
1805 // soft pions
1806 WTrackParameter wvpiTrk3, wvpiTrk4;
1807 wvpiTrk3 = WTrackParameter( mpi, pitrk1->getZHelix(), pitrk1->getZError() );
1808 wvpiTrk4 = WTrackParameter( mpi, pitrk2->getZHelix(), pitrk2->getZError() );
1809
1810 if ( m_selectPsipRhoPi )
1811 {
1812 KinematicFit* kmfit = KinematicFit::instance();
1813 kmfit->init();
1814 kmfit->AddTrack( 0, wvpiTrk1 );
1815 kmfit->AddTrack( 1, wvpiTrk2 );
1816 kmfit->AddTrack( 2, 0.0, gShr1 );
1817 kmfit->AddTrack( 3, 0.0, gShr2 );
1818 kmfit->AddTrack( 4, wvpiTrk3 );
1819 kmfit->AddTrack( 5, wvpiTrk4 );
1820 kmfit->AddFourMomentum( 0, ecms );
1821
1822 bool oksq1 = kmfit->Fit();
1823 double chi1 = kmfit->chisq();
1824 //
1825
1826 if ( oksq1 && chi1 > 0 && chi1 < 100 )
1827 {
1828 m_isPsipRhoPi = true;
1829 m_psipRhoPiNumber++;
1830 }
1831 }
1832 if ( m_selectPsipKstarK )
1833 {
1834 KinematicFit* kmfit2 = KinematicFit::instance();
1835 kmfit2->init();
1836 kmfit2->AddTrack( 0, wvkTrk1 );
1837 kmfit2->AddTrack( 1, wvkTrk2 );
1838 kmfit2->AddTrack( 2, 0.0, gShr3 );
1839 kmfit2->AddTrack( 3, 0.0, gShr4 );
1840 kmfit2->AddTrack( 4, wvpiTrk3 );
1841 kmfit2->AddTrack( 5, wvpiTrk4 );
1842 kmfit2->AddFourMomentum( 0, ecms );
1843
1844 bool oksq2 = kmfit2->Fit();
1845 double chi2 = kmfit2->chisq();
1846
1847 if ( oksq2 && chi2 > 0 && chi2 < 200 &&
1848 ( ( fabs( mkstarp - 0.89166 ) < 0.1 && fabs( mkstarm - 0.89166 ) > 0.1 ) ||
1849 ( fabs( mkstarm - 0.89166 ) < 0.1 && fabs( mkstarp - 0.89166 ) > 0.1 ) ) )
1850 {
1851 m_isPsipKstarK = true;
1852 m_psipKstarKNumber++;
1853 }
1854 }
1855
1856 } // end of loose cuts
1857 }
1858
1859 } // recoil jpsi to 2tracks and 1 pi0
1860
1861 // jpsi to pppipi
1862 if ( m_selectPsipPPPiPi && ( bestmass > 3.075 && bestmass < 3.120 ) && nGood == 6 &&
1863 jpsicharge == 0 && nPi0 == 0 )
1864 {
1865
1866 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iJPsiP[0];
1867 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iJPsiP[1];
1868 EvtRecTrackIterator itthr = evtRecTrkCol->begin() + iJPsiN[0];
1869 EvtRecTrackIterator itfor = evtRecTrkCol->begin() + iJPsiN[1];
1870 RecMdcKalTrack* trk1 = ( *itone )->mdcKalTrack();
1871 RecMdcKalTrack* trk2 = ( *ittwo )->mdcKalTrack();
1872 RecMdcKalTrack* trk3 = ( *itthr )->mdcKalTrack();
1873 RecMdcKalTrack* trk4 = ( *itfor )->mdcKalTrack();
1874
1875 HepLorentzVector p4trkpp;
1876 HepLorentzVector p4trkpm;
1877 HepLorentzVector p4trkpip;
1878 HepLorentzVector p4trkpim;
1879
1880 // hypothesis 1
1881
1882 trk1->setPidType( RecMdcKalTrack::proton );
1883 p4trkpp.setPx( Px( trk1 ) );
1884 p4trkpp.setPy( Py( trk1 ) );
1885 p4trkpp.setPz( Pz( trk1 ) );
1886 p4trkpp.setE( sqrt( pow( P( trk1 ), 2 ) + mproton * mproton ) );
1887
1888 trk3->setPidType( RecMdcKalTrack::proton );
1889 p4trkpm.setPx( Px( trk3 ) );
1890 p4trkpm.setPy( Py( trk3 ) );
1891 p4trkpm.setPz( Pz( trk3 ) );
1892 p4trkpm.setE( sqrt( pow( P( trk3 ), 2 ) + mproton * mproton ) );
1893
1894 trk2->setPidType( RecMdcKalTrack::pion );
1895 p4trkpip.setPx( Px( trk2 ) );
1896 p4trkpip.setPy( Py( trk2 ) );
1897 p4trkpip.setPz( Pz( trk2 ) );
1898 p4trkpip.setE( sqrt( pow( P( trk2 ), 2 ) + mpi * mpi ) );
1899
1900 trk4->setPidType( RecMdcKalTrack::pion );
1901 p4trkpim.setPx( Px( trk4 ) );
1902 p4trkpim.setPy( Py( trk4 ) );
1903 p4trkpim.setPz( Pz( trk4 ) );
1904 p4trkpim.setE( sqrt( pow( P( trk4 ), 2 ) + mpi * mpi ) );
1905
1906 double jpsimass1 = ( p4trkpp + p4trkpm + p4trkpip + p4trkpim ).m();
1907
1908 // hypothesis 2
1909
1910 trk1->setPidType( RecMdcKalTrack::proton );
1911 p4trkpp.setPx( Px( trk1 ) );
1912 p4trkpp.setPy( Py( trk1 ) );
1913 p4trkpp.setPz( Pz( trk1 ) );
1914 p4trkpp.setE( sqrt( pow( P( trk1 ), 2 ) + mproton * mproton ) );
1915
1916 trk4->setPidType( RecMdcKalTrack::proton );
1917 p4trkpm.setPx( Px( trk4 ) );
1918 p4trkpm.setPy( Py( trk4 ) );
1919 p4trkpm.setPz( Pz( trk4 ) );
1920 p4trkpm.setE( sqrt( pow( P( trk4 ), 2 ) + mproton * mproton ) );
1921
1922 trk2->setPidType( RecMdcKalTrack::pion );
1923 p4trkpip.setPx( Px( trk2 ) );
1924 p4trkpip.setPy( Py( trk2 ) );
1925 p4trkpip.setPz( Pz( trk2 ) );
1926 p4trkpip.setE( sqrt( pow( P( trk2 ), 2 ) + mpi * mpi ) );
1927
1928 trk3->setPidType( RecMdcKalTrack::pion );
1929 p4trkpim.setPx( Px( trk3 ) );
1930 p4trkpim.setPy( Py( trk3 ) );
1931 p4trkpim.setPz( Pz( trk3 ) );
1932 p4trkpim.setE( sqrt( pow( P( trk3 ), 2 ) + mpi * mpi ) );
1933
1934 double jpsimass2 = ( p4trkpp + p4trkpm + p4trkpip + p4trkpim ).m();
1935
1936 // hypothesis 3
1937
1938 trk2->setPidType( RecMdcKalTrack::proton );
1939 p4trkpp.setPx( Px( trk2 ) );
1940 p4trkpp.setPy( Py( trk2 ) );
1941 p4trkpp.setPz( Pz( trk2 ) );
1942 p4trkpp.setE( sqrt( pow( P( trk2 ), 2 ) + mproton * mproton ) );
1943
1944 trk3->setPidType( RecMdcKalTrack::proton );
1945 p4trkpm.setPx( Px( trk3 ) );
1946 p4trkpm.setPy( Py( trk3 ) );
1947 p4trkpm.setPz( Pz( trk3 ) );
1948 p4trkpm.setE( sqrt( pow( P( trk3 ), 2 ) + mproton * mproton ) );
1949
1950 trk1->setPidType( RecMdcKalTrack::pion );
1951 p4trkpip.setPx( Px( trk1 ) );
1952 p4trkpip.setPy( Py( trk1 ) );
1953 p4trkpip.setPz( Pz( trk1 ) );
1954 p4trkpip.setE( sqrt( pow( P( trk1 ), 2 ) + mpi * mpi ) );
1955
1956 trk4->setPidType( RecMdcKalTrack::pion );
1957 p4trkpim.setPx( Px( trk4 ) );
1958 p4trkpim.setPy( Py( trk4 ) );
1959 p4trkpim.setPz( Pz( trk4 ) );
1960 p4trkpim.setE( sqrt( pow( P( trk4 ), 2 ) + mpi * mpi ) );
1961
1962 double jpsimass3 = ( p4trkpp + p4trkpm + p4trkpip + p4trkpim ).m();
1963
1964 // hypothesis 4
1965
1966 trk2->setPidType( RecMdcKalTrack::proton );
1967 p4trkpp.setPx( Px( trk2 ) );
1968 p4trkpp.setPy( Py( trk2 ) );
1969 p4trkpp.setPz( Pz( trk2 ) );
1970 p4trkpp.setE( sqrt( pow( P( trk2 ), 2 ) + mproton * mproton ) );
1971
1972 trk4->setPidType( RecMdcKalTrack::proton );
1973 p4trkpm.setPx( Px( trk4 ) );
1974 p4trkpm.setPy( Py( trk4 ) );
1975 p4trkpm.setPz( Pz( trk4 ) );
1976 p4trkpm.setE( sqrt( pow( P( trk4 ), 2 ) + mproton * mproton ) );
1977
1978 trk1->setPidType( RecMdcKalTrack::pion );
1979 p4trkpip.setPx( Px( trk1 ) );
1980 p4trkpip.setPy( Py( trk1 ) );
1981 p4trkpip.setPz( Pz( trk1 ) );
1982 p4trkpip.setE( sqrt( pow( P( trk1 ), 2 ) + mpi * mpi ) );
1983
1984 trk3->setPidType( RecMdcKalTrack::pion );
1985 p4trkpim.setPx( Px( trk3 ) );
1986 p4trkpim.setPy( Py( trk3 ) );
1987 p4trkpim.setPz( Pz( trk3 ) );
1988 p4trkpim.setE( sqrt( pow( P( trk3 ), 2 ) + mpi * mpi ) );
1989
1990 double jpsimass4 = ( p4trkpp + p4trkpm + p4trkpip + p4trkpim ).m();
1991
1992 if ( fabs( jpsimass1 - 3.075 ) <= 0.075 || fabs( jpsimass2 - 3.075 ) <= 0.075 ||
1993 fabs( jpsimass3 - 3.075 ) <= 0.075 || fabs( jpsimass4 - 3.075 ) <= 0.075 )
1994 {
1995
1996 // tighten cuts
1997 double chi1, chi2, chi3, chi4;
1998 int type = 0;
1999 double chimin = -999;
2000
2001 WTrackParameter wvpiTrk1, wvpiTrk2, wvpiTrk3, wvpiTrk4, wvpiTrk5, wvpiTrk6;
2002
2003 {
2004 wvpiTrk1 = WTrackParameter( mproton, trk1->getZHelix(), trk1->getZError() );
2005 wvpiTrk2 = WTrackParameter( mpi, trk2->getZHelix(), trk2->getZError() );
2006 wvpiTrk3 = WTrackParameter( mproton, trk3->getZHelix(), trk3->getZError() );
2007 wvpiTrk4 = WTrackParameter( mpi, trk4->getZHelix(), trk4->getZError() );
2008 wvpiTrk5 = WTrackParameter( mpi, pitrk1->getZHelix(), pitrk1->getZError() );
2009 wvpiTrk6 = WTrackParameter( mpi, pitrk2->getZHelix(), pitrk2->getZError() );
2010
2011 KinematicFit* kmfit = KinematicFit::instance();
2012 kmfit->init();
2013 kmfit->AddTrack( 0, wvpiTrk1 );
2014 kmfit->AddTrack( 1, wvpiTrk2 );
2015 kmfit->AddTrack( 2, wvpiTrk3 );
2016 kmfit->AddTrack( 3, wvpiTrk4 );
2017 kmfit->AddTrack( 4, wvpiTrk5 );
2018 kmfit->AddTrack( 5, wvpiTrk6 );
2019 kmfit->AddFourMomentum( 0, ecms );
2020
2021 bool oksq1 = kmfit->Fit();
2022 chi1 = kmfit->chisq();
2023 if ( oksq1 )
2024 {
2025 chimin = chi1;
2026 type = 1;
2027 }
2028 }
2029
2030 {
2031 wvpiTrk1 = WTrackParameter( mproton, trk1->getZHelix(), trk1->getZError() );
2032 wvpiTrk2 = WTrackParameter( mpi, trk2->getZHelix(), trk2->getZError() );
2033 wvpiTrk3 = WTrackParameter( mpi, trk3->getZHelix(), trk3->getZError() );
2034 wvpiTrk4 = WTrackParameter( mproton, trk4->getZHelix(), trk4->getZError() );
2035 wvpiTrk5 = WTrackParameter( mpi, pitrk1->getZHelix(), pitrk1->getZError() );
2036 wvpiTrk6 = WTrackParameter( mpi, pitrk2->getZHelix(), pitrk2->getZError() );
2037
2038 KinematicFit* kmfit = KinematicFit::instance();
2039 kmfit->init();
2040 kmfit->AddTrack( 0, wvpiTrk1 );
2041 kmfit->AddTrack( 1, wvpiTrk2 );
2042 kmfit->AddTrack( 2, wvpiTrk3 );
2043 kmfit->AddTrack( 3, wvpiTrk4 );
2044 kmfit->AddTrack( 4, wvpiTrk5 );
2045 kmfit->AddTrack( 5, wvpiTrk6 );
2046 kmfit->AddFourMomentum( 0, ecms );
2047
2048 bool oksq1 = kmfit->Fit();
2049 chi2 = kmfit->chisq();
2050 if ( oksq1 )
2051 {
2052 if ( type == 0 )
2053 {
2054 chimin = chi2;
2055 type = 2;
2056 }
2057 else if ( chi2 < chimin )
2058 {
2059 chimin = chi2;
2060 type = 2;
2061 }
2062 }
2063 //
2064 }
2065
2066 {
2067 wvpiTrk1 = WTrackParameter( mpi, trk1->getZHelix(), trk1->getZError() );
2068 wvpiTrk2 = WTrackParameter( mproton, trk2->getZHelix(), trk2->getZError() );
2069 wvpiTrk3 = WTrackParameter( mpi, trk3->getZHelix(), trk3->getZError() );
2070 wvpiTrk4 = WTrackParameter( mproton, trk4->getZHelix(), trk4->getZError() );
2071 wvpiTrk5 = WTrackParameter( mpi, pitrk1->getZHelix(), pitrk1->getZError() );
2072 wvpiTrk6 = WTrackParameter( mpi, pitrk2->getZHelix(), pitrk2->getZError() );
2073
2074 KinematicFit* kmfit = KinematicFit::instance();
2075 kmfit->init();
2076 kmfit->AddTrack( 0, wvpiTrk1 );
2077 kmfit->AddTrack( 1, wvpiTrk2 );
2078 kmfit->AddTrack( 2, wvpiTrk3 );
2079 kmfit->AddTrack( 3, wvpiTrk4 );
2080 kmfit->AddTrack( 4, wvpiTrk5 );
2081 kmfit->AddTrack( 5, wvpiTrk6 );
2082 kmfit->AddFourMomentum( 0, ecms );
2083
2084 bool oksq1 = kmfit->Fit();
2085 chi3 = kmfit->chisq();
2086 if ( oksq1 )
2087 {
2088 if ( type == 0 )
2089 {
2090 chimin = chi3;
2091 type = 3;
2092 }
2093 else if ( chi3 < chimin )
2094 {
2095 chimin = chi3;
2096 type = 3;
2097 }
2098 }
2099 // delete kmfit;
2100 }
2101
2102 {
2103 wvpiTrk1 = WTrackParameter( mpi, trk1->getZHelix(), trk1->getZError() );
2104 wvpiTrk2 = WTrackParameter( mproton, trk2->getZHelix(), trk2->getZError() );
2105 wvpiTrk3 = WTrackParameter( mproton, trk3->getZHelix(), trk3->getZError() );
2106 wvpiTrk4 = WTrackParameter( mpi, trk4->getZHelix(), trk4->getZError() );
2107 wvpiTrk5 = WTrackParameter( mpi, pitrk1->getZHelix(), pitrk1->getZError() );
2108 wvpiTrk6 = WTrackParameter( mpi, pitrk2->getZHelix(), pitrk2->getZError() );
2109
2110 KinematicFit* kmfit = KinematicFit::instance();
2111 kmfit->init();
2112 kmfit->AddTrack( 0, wvpiTrk1 );
2113 kmfit->AddTrack( 1, wvpiTrk2 );
2114 kmfit->AddTrack( 2, wvpiTrk3 );
2115 kmfit->AddTrack( 3, wvpiTrk4 );
2116 kmfit->AddTrack( 4, wvpiTrk5 );
2117 kmfit->AddTrack( 5, wvpiTrk6 );
2118 kmfit->AddFourMomentum( 0, ecms );
2119
2120 bool oksq1 = kmfit->Fit();
2121 chi4 = kmfit->chisq();
2122 if ( oksq1 )
2123 {
2124 if ( type == 0 )
2125 {
2126 chimin = chi4;
2127 type = 4;
2128 }
2129 else if ( chi4 < chimin )
2130 {
2131 chimin = chi4;
2132 type = 4;
2133 }
2134 }
2135 }
2136
2137 if ( chimin > 0 && chimin < 200 )
2138 {
2139 m_isPsipPPPiPi = true;
2140 m_psipPPPiPiNumber++;
2141 }
2142
2143 } // end of loose cuts
2144
2145 } // end of selecting recol jpsi to pppipi
2146
2147 } // end of recoil jpsi selection
2148
2149 // select psi'' events using dtaging
2150
2151 if ( m_selectPsippCand )
2152 {
2153
2154 SmartDataPtr<EvtRecDTagCol> evtRecDTagCol( eventSvc(), EventModel::EvtRec::EvtRecDTagCol );
2155 if ( !evtRecDTagCol )
2156 {
2157 log << MSG::FATAL << "Could not find EvtRecDTagCol" << endmsg;
2158 return StatusCode::FAILURE;
2159 }
2160 if ( evtRecDTagCol->size() > 0 )
2161 {
2162
2163 EvtRecDTagCol::iterator m_iterbegin = evtRecDTagCol->begin();
2164 EvtRecDTagCol::iterator m_iterend = evtRecDTagCol->end();
2165 for ( EvtRecDTagCol::iterator iter = m_iterbegin; iter != m_iterend; iter++ )
2166 {
2167
2168 if ( ( ( *iter )->decayMode() == EvtRecDTag::kD0toKPiPi0 &&
2169 fabs( ( *iter )->mBC() - 1.865 ) < 0.006 && ( *iter )->deltaE() > -0.05 &&
2170 ( *iter )->deltaE() < 0.03 ) ||
2171 ( ( *iter )->decayMode() == EvtRecDTag::kD0toKPiPi0Pi0 &&
2172 fabs( ( *iter )->mBC() - 1.865 ) < 0.006 && ( *iter )->deltaE() > -0.05 &&
2173 ( *iter )->deltaE() < 0.03 ) ||
2174 ( ( *iter )->decayMode() == EvtRecDTag::kD0toKPiPiPi &&
2175 fabs( ( *iter )->mBC() - 1.865 ) < 0.006 && ( *iter )->deltaE() > -0.03 &&
2176 ( *iter )->deltaE() < 0.03 ) ||
2177 ( ( *iter )->decayMode() == EvtRecDTag::kD0toKPiPiPiPi0 &&
2178 fabs( ( *iter )->mBC() - 1.865 ) < 0.006 && ( *iter )->deltaE() > -0.05 &&
2179 ( *iter )->deltaE() < 0.03 ) ||
2180 ( ( *iter )->decayMode() == EvtRecDTag::kD0toKsPiPi &&
2181 fabs( ( *iter )->mBC() - 1.865 ) < 0.006 && ( *iter )->deltaE() > -0.03 &&
2182 ( *iter )->deltaE() < 0.03 ) ||
2183 ( ( *iter )->decayMode() == EvtRecDTag::kD0toKsPiPiPi0 &&
2184 fabs( ( *iter )->mBC() - 1.865 ) < 0.006 && ( *iter )->deltaE() > -0.05 &&
2185 ( *iter )->deltaE() < 0.03 ) ||
2186 ( ( *iter )->decayMode() == EvtRecDTag::kDptoKPiPi &&
2187 fabs( ( *iter )->mBC() - 1.87 ) < 0.006 && ( *iter )->deltaE() > -0.03 &&
2188 ( *iter )->deltaE() < 0.03 ) ||
2189 ( ( *iter )->decayMode() == EvtRecDTag::kDptoKPiPiPi0 &&
2190 fabs( ( *iter )->mBC() - 1.87 ) < 0.006 && ( *iter )->deltaE() > -0.05 &&
2191 ( *iter )->deltaE() < 0.03 ) ||
2192 ( ( *iter )->decayMode() == EvtRecDTag::kDptoKsPiPi0 &&
2193 fabs( ( *iter )->mBC() - 1.87 ) < 0.006 && ( *iter )->deltaE() > -0.05 &&
2194 ( *iter )->deltaE() < 0.03 ) ||
2195 ( ( *iter )->decayMode() == EvtRecDTag::kDptoKsPiPiPi &&
2196 fabs( ( *iter )->mBC() - 1.87 ) < 0.006 && ( *iter )->deltaE() > -0.03 &&
2197 ( *iter )->deltaE() < 0.03 ) )
2198 {
2199 m_isPsippCand = true;
2200 m_psippCandNumber++;
2201 break;
2202 }
2203
2204 } // end of looping D modes
2205
2206 } // end of non-zero dtag list
2207
2208 } // end of selecting psi'' events
2209
2210 // -------- Write to root file
2211
2212 if ( m_selectRadBhabha && m_isRadBhabha ) m_subalg3->execute();
2213 if ( m_selectGGEE && m_isGGEE ) m_subalg4->execute();
2214 if ( m_selectGG4Pi && m_isGG4Pi ) m_subalg5->execute();
2215 if ( m_selectRadBhabhaT && m_isRadBhabhaT ) m_subalg6->execute();
2216 if ( m_selectRhoPi && m_isRhoPi ) m_subalg7->execute();
2217 if ( m_selectKstarK && m_isKstarK ) m_subalg8->execute();
2218 if ( m_selectPPPiPi && m_isPPPiPi ) m_subalg9->execute();
2219 if ( m_selectRecoJpsi && m_isRecoJpsi ) m_subalg10->execute();
2220 if ( m_selectBhabha && m_isBhabha ) m_subalg11->execute();
2221 if ( m_selectDimu && m_isDimu ) m_subalg12->execute();
2222 if ( m_selectCosmic && m_isCosmic ) m_subalg13->execute();
2223 if ( m_selectGenProton && m_isGenProton ) m_subalg14->execute();
2224 if ( m_selectPsipRhoPi && m_isPsipRhoPi ) m_subalg15->execute();
2225 if ( m_selectPsipKstarK && m_isPsipKstarK ) m_subalg16->execute();
2226 if ( m_selectPsipPPPiPi && m_isPsipPPPiPi ) m_subalg17->execute();
2227 if ( m_selectPsippCand && m_isPsippCand ) m_subalg18->execute();
2228
2229 // if(m_output) {
2230 // m_runnb = run;
2231 // m_evtnb = event;
2232 // m_esum = esum;
2233 // m_eemc = eemc;
2234 // m_etot = etot;
2235 // m_nCharge = nCharge;
2236 // m_nGood = nGood;
2237 // m_nGam = nGam;
2238 // m_ptot = ptot;
2239 // m_pp = pp;
2240 // m_pm = pm;
2241 // m_maxE = maxE;
2242 // m_secE = secE;
2243 // m_dThe = dthe;
2244 // m_dPhi = dphi;
2245 // m_mdcHit1 = mdcHit1;
2246 // m_mdcHit2 = mdcHit2;
2247 // m_tuple0->write();
2248 // }
2249
2250 return StatusCode::SUCCESS;
2251}
p2
double p2[4]
Definition BabayagaNLO.cxx:41
p1
double p1[4]
Definition BabayagaNLO.cxx:41
Px
double Px(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:56
Pz
double Pz(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:75
Phi
double Phi(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:96
P
double P(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:84
Py
double Py(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:65
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
mkaon
const double mkaon
Definition DSemilepAlg.cxx:35
etot
Double_t etot
Definition DataBase/tau_mode.c:8
time
Double_t time
Definition DataBase/tau_mode.c:7
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
mpi
double mpi
Definition EvtPhokharaDef.hh:17
pi
double pi
Definition EvtPhokharaDef.hh:17
ebeam
double ebeam
Definition EvtPhokharaDef.hh:18
mpi0
double mpi0
Definition EvtPhokharaDef.hh:17
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:83
Vp4
std::vector< HepLorentzVector > Vp4
Definition Gam4pikp.cxx:38
Vint
std::vector< int > Vint
Definition Gam4pikp.cxx:37
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
mproton
const double mproton
Definition PipiJpsi.cxx:49
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
CalibEventSelect::readbeamEfromDb
void readbeamEfromDb(int runNo, double &beamE)
Definition CalibEventSelect.cxx:2336
DstMdcKalTrack::setPx
void setPx(const double px, const int pid)
Definition Event/DstEvent/include/DstEvent/DstMdcKalTrack.h:99
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
KinematicFit::instance
static KinematicFit * instance()
Definition KinematicFit.cxx:18
KinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition KinematicFit.cxx:341
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition TofHitStatus.cxx:115
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:21
ecms
const double ecms
Definition inclkstar.cxx:26

◆ finalize()

StatusCode CalibEventSelect::finalize ( )

Definition at line 2254 of file CalibEventSelect.cxx.

2254 {
2255
2256 MsgStream log( msgSvc(), name() );
2257 log << MSG::INFO << "in finalize()" << endmsg;
2258
2259 cout << "Total events: " << m_events << endl;
2260
2261 if ( m_selectRadBhabha ) { cout << "Selected rad bhabha: " << m_radBhabhaNumber << endl; }
2262
2263 if ( m_selectGGEE ) { cout << "Selected ggee events: " << m_GGEENumber << endl; }
2264
2265 if ( m_selectGG4Pi ) { cout << "Selected gg4pi events: " << m_GG4PiNumber << endl; }
2266
2267 if ( m_selectRadBhabhaT )
2268 { cout << "Selected rad bhabha tight: " << m_radBhabhaTNumber << endl; }
2269
2270 if ( m_selectRhoPi ) { cout << "Selected rhopi events: " << m_rhoPiNumber << endl; }
2271
2272 if ( m_selectKstarK ) { cout << "Selected kstark events: " << m_kstarKNumber << endl; }
2273
2274 if ( m_selectPPPiPi ) { cout << "Selected pppipi events: " << m_ppPiPiNumber << endl; }
2275
2276 if ( m_selectRecoJpsi )
2277 { cout << "Selected recoil jsi events: " << m_recoJpsiNumber << endl; }
2278
2279 if ( m_selectBhabha ) { cout << "Selected bhabha events: " << m_bhabhaNumber << endl; }
2280
2281 if ( m_selectDimu ) { cout << "Selected dimu events: " << m_dimuNumber << endl; }
2282
2283 if ( m_selectCosmic ) { cout << "Selected cosmic events: " << m_cosmicNumber << endl; }
2284
2285 if ( m_selectGenProton )
2286 { cout << "Selected generic proton events: " << m_genProtonNumber << endl; }
2287
2288 if ( m_selectPsipRhoPi )
2289 { cout << "Selected recoil rhopi events: " << m_psipRhoPiNumber << endl; }
2290
2291 if ( m_selectPsipKstarK )
2292 { cout << "Selected recoil kstark events: " << m_psipKstarKNumber << endl; }
2293
2294 if ( m_selectPsipPPPiPi )
2295 { cout << "Selected recoil pppipi events: " << m_psipPPPiPiNumber << endl; }
2296
2297 if ( m_selectPsippCand )
2298 { cout << "Selected psi'' candi events: " << m_psippCandNumber << endl; }
2299
2300 return StatusCode::SUCCESS;
2301}

◆ initialize()

StatusCode CalibEventSelect::initialize ( )

Definition at line 159 of file CalibEventSelect.cxx.

159 {
160 MsgStream log( msgSvc(), name() );
161
162 log << MSG::INFO << "in initialize()" << endmsg;
163
164 m_run = 0;
165 // m_ecm=3.1;
166
167 h_nGood = histoSvc()->book( "1D/nGoodtracks", 1, "num of good chaged tracks", 20, 0, 20 );
168 h_nCharge = histoSvc()->book( "1D/nCharge", 1, "net charge", 20, -10, 10 );
169 h_pmaxobeam = histoSvc()->book( "1D/pmax", 1, "pmax over beam ratio", 100, 0, 3 );
170 h_eopmax = histoSvc()->book( "1D/eopmax", 1, "e over pmax ratio", 100, 0, 3 );
171 h_eopsec = histoSvc()->book( "1D/eopsec", 1, "e over psec ratio", 100, 0, 3 );
172 h_deltap = histoSvc()->book( "1D/deltap", 1, "pmax minus psec", 100, -3, 3 );
173 h_esumoecm =
174 histoSvc()->book( "1D/esumoverecm", 1, "total energy over ecm ratio", 100, 0, 3 );
175 h_egmax = histoSvc()->book( "1D/egmax", 1, "most energetic photon energy", 100, 0, 3 );
176 h_deltaphi = histoSvc()->book( "1D/deltaphi", 1, "phi_pmax minus phi_sec", 150, -4, 4 );
177 h_Pt = histoSvc()->book( "1D/Pt", 1, "total Transverse Momentum", 200, -1, 4 );
178
179 StatusCode sc;
180
181 log << MSG::INFO << "creating sub-algorithms...." << endmsg;
182
183 if ( m_writeDst )
184 {
185 sc = createSubAlgorithm( "EventWriter", "WriteDst", m_subalg1 );
186 if ( sc.isFailure() )
187 {
188 log << MSG::ERROR << "Error creating Sub-Algorithm WriteDst" << endmsg;
189 return sc;
190 }
191 else { log << MSG::INFO << "Success creating Sub-Algorithm WriteDst" << endmsg; }
192 }
193
194 if ( m_writeRec )
195 {
196 sc = createSubAlgorithm( "EventWriter", "WriteRec", m_subalg2 );
197 if ( sc.isFailure() )
198 {
199 log << MSG::ERROR << "Error creating Sub-Algorithm WriteRec" << endmsg;
200 return sc;
201 }
202 else { log << MSG::INFO << "Success creating Sub-Algorithm WriteRec" << endmsg; }
203 }
204
205 if ( m_selectRadBhabha )
206 {
207 sc = createSubAlgorithm( "EventWriter", "SelectRadBhabha", m_subalg3 );
208 if ( sc.isFailure() )
209 {
210 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRadBhabha" << endmsg;
211 return sc;
212 }
213 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectRadBhabha" << endmsg; }
214 }
215
216 if ( m_selectGGEE )
217 {
218 sc = createSubAlgorithm( "EventWriter", "SelectGGEE", m_subalg4 );
219 if ( sc.isFailure() )
220 {
221 log << MSG::ERROR << "Error creating Sub-Algorithm SelectGGEE" << endmsg;
222 return sc;
223 }
224 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectGGEE" << endmsg; }
225 }
226
227 if ( m_selectGG4Pi )
228 {
229 sc = createSubAlgorithm( "EventWriter", "SelectGG4Pi", m_subalg5 );
230 if ( sc.isFailure() )
231 {
232 log << MSG::ERROR << "Error creating Sub-Algorithm SelectGG4Pi" << endmsg;
233 return sc;
234 }
235 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectGG4Pi" << endmsg; }
236 }
237
238 if ( m_selectRadBhabhaT )
239 {
240 sc = createSubAlgorithm( "EventWriter", "SelectRadBhabhaT", m_subalg6 );
241 if ( sc.isFailure() )
242 {
243 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRadBhabhaT" << endmsg;
244 return sc;
245 }
246 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectRadBhabhaT" << endmsg; }
247 }
248
249 if ( m_selectRhoPi )
250 {
251 sc = createSubAlgorithm( "EventWriter", "SelectRhoPi", m_subalg7 );
252 if ( sc.isFailure() )
253 {
254 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRhoPi" << endmsg;
255 return sc;
256 }
257 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectRhoPi" << endmsg; }
258 }
259
260 if ( m_selectKstarK )
261 {
262 sc = createSubAlgorithm( "EventWriter", "SelectKstarK", m_subalg8 );
263 if ( sc.isFailure() )
264 {
265 log << MSG::ERROR << "Error creating Sub-Algorithm SelectKstarK" << endmsg;
266 return sc;
267 }
268 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectKstarK" << endmsg; }
269 }
270
271 if ( m_selectPPPiPi )
272 {
273 sc = createSubAlgorithm( "EventWriter", "SelectPPPiPi", m_subalg9 );
274 if ( sc.isFailure() )
275 {
276 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPPPiPi" << endmsg;
277 return sc;
278 }
279 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectPPPiPi" << endmsg; }
280 }
281
282 if ( m_selectRecoJpsi )
283 {
284 sc = createSubAlgorithm( "EventWriter", "SelectRecoJpsi", m_subalg10 );
285 if ( sc.isFailure() )
286 {
287 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRecoJpsi" << endmsg;
288 return sc;
289 }
290 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectRecoJpsi" << endmsg; }
291 }
292
293 if ( m_selectBhabha )
294 {
295 sc = createSubAlgorithm( "EventWriter", "SelectBhabha", m_subalg11 );
296 if ( sc.isFailure() )
297 {
298 log << MSG::ERROR << "Error creating Sub-Algorithm SelectBhabha" << endmsg;
299 return sc;
300 }
301 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectBhabha" << endmsg; }
302 }
303
304 if ( m_selectDimu )
305 {
306 sc = createSubAlgorithm( "EventWriter", "SelectDimu", m_subalg12 );
307 if ( sc.isFailure() )
308 {
309 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDimu" << endmsg;
310 return sc;
311 }
312 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectDimu" << endmsg; }
313 }
314
315 if ( m_selectCosmic )
316 {
317 sc = createSubAlgorithm( "EventWriter", "SelectCosmic", m_subalg13 );
318 if ( sc.isFailure() )
319 {
320 log << MSG::ERROR << "Error creating Sub-Algorithm SelectCosmic" << endmsg;
321 return sc;
322 }
323 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectCosmic" << endmsg; }
324 }
325
326 if ( m_selectGenProton )
327 {
328 sc = createSubAlgorithm( "EventWriter", "SelectGenProton", m_subalg14 );
329 if ( sc.isFailure() )
330 {
331 log << MSG::ERROR << "Error creating Sub-Algorithm SelectGenProton" << endmsg;
332 return sc;
333 }
334 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectGenProton" << endmsg; }
335 }
336
337 if ( m_selectPsipRhoPi )
338 {
339 sc = createSubAlgorithm( "EventWriter", "SelectPsipRhoPi", m_subalg15 );
340 if ( sc.isFailure() )
341 {
342 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsipRhoPi" << endmsg;
343 return sc;
344 }
345 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectPsipRhoPi" << endmsg; }
346 }
347
348 if ( m_selectPsipKstarK )
349 {
350 sc = createSubAlgorithm( "EventWriter", "SelectPsipKstarK", m_subalg16 );
351 if ( sc.isFailure() )
352 {
353 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsipKstarK" << endmsg;
354 return sc;
355 }
356 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectPsipKstarK" << endmsg; }
357 }
358
359 if ( m_selectPsipPPPiPi )
360 {
361 sc = createSubAlgorithm( "EventWriter", "SelectPsipPPPiPi", m_subalg17 );
362 if ( sc.isFailure() )
363 {
364 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsipPPPiPi" << endmsg;
365 return sc;
366 }
367 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectPsipPPPiPi" << endmsg; }
368 }
369
370 if ( m_selectPsippCand )
371 {
372 sc = createSubAlgorithm( "EventWriter", "SelectPsippCand", m_subalg18 );
373 if ( sc.isFailure() )
374 {
375 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsippCand" << endmsg;
376 return sc;
377 }
378 else { log << MSG::INFO << "Success creating Sub-Algorithm SelectPsippCand" << endmsg; }
379 }
380
381 if ( m_output )
382 {
383 StatusCode status;
384 NTuplePtr nt0( ntupleSvc(), "FILE1/hadron" );
385 if ( nt0 ) m_tuple0 = nt0;
386 else
387 {
388 m_tuple0 = ntupleSvc()->book( "FILE1/hadron", CLID_ColumnWiseTuple, "N-Tuple example" );
389 if ( m_tuple0 )
390 {
391 status = m_tuple0->addItem( "esum", m_esum );
392 status = m_tuple0->addItem( "eemc", m_eemc );
393 status = m_tuple0->addItem( "etot", m_etot );
394 status = m_tuple0->addItem( "nGood", m_nGood );
395 status = m_tuple0->addItem( "nCharge", m_nCharge );
396 status = m_tuple0->addItem( "nGam", m_nGam );
397 status = m_tuple0->addItem( "ptot", m_ptot );
398 status = m_tuple0->addItem( "pp", m_pp );
399 status = m_tuple0->addItem( "pm", m_pm );
400 status = m_tuple0->addItem( "run", m_runnb );
401 status = m_tuple0->addItem( "event", m_evtnb );
402 status = m_tuple0->addItem( "maxE", m_maxE );
403 status = m_tuple0->addItem( "secE", m_secE );
404 status = m_tuple0->addItem( "dthe", m_dthe );
405 status = m_tuple0->addItem( "dphi", m_dphi );
406 status = m_tuple0->addItem( "mdcHit1", m_mdcHit1 );
407 status = m_tuple0->addItem( "mdcHit2", m_mdcHit2 );
408 }
409 else
410 {
411 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple0 ) << endmsg;
412 return StatusCode::FAILURE;
413 }
414 }
415
416 NTuplePtr nt1( ntupleSvc(), "FILE1/vxyz" );
417 if ( nt1 ) m_tuple1 = nt1;
418 else
419 {
420 m_tuple1 = ntupleSvc()->book( "FILE1/vxyz", CLID_ColumnWiseTuple, "ks N-Tuple example" );
421 if ( m_tuple1 )
422 {
423 status = m_tuple1->addItem( "vx0", m_vx0 );
424 status = m_tuple1->addItem( "vy0", m_vy0 );
425 status = m_tuple1->addItem( "vz0", m_vz0 );
426 status = m_tuple1->addItem( "vr0", m_vr0 );
427 status = m_tuple1->addItem( "theta0", m_theta0 );
428 status = m_tuple1->addItem( "p0", m_p0 );
429 status = m_tuple1->addItem( "pt0", m_pt0 );
430 }
431 else
432 {
433 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
434 return StatusCode::FAILURE;
435 }
436 }
437
438 NTuplePtr nt2( ntupleSvc(), "FILE1/photon" );
439 if ( nt2 ) m_tuple2 = nt2;
440 else
441 {
442 m_tuple2 =
443 ntupleSvc()->book( "FILE1/photon", CLID_ColumnWiseTuple, "ks N-Tuple example" );
444 if ( m_tuple2 )
445 {
446 status = m_tuple2->addItem( "dthe", m_dthe );
447 status = m_tuple2->addItem( "dphi", m_dphi );
448 status = m_tuple2->addItem( "dang", m_dang );
449 status = m_tuple2->addItem( "eraw", m_eraw );
450 }
451 else
452 {
453 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple2 ) << endmsg;
454 return StatusCode::FAILURE;
455 }
456 }
457 }
458 //
459 //--------end of book--------
460 //
461
462 m_events = 0;
463 m_radBhabhaNumber = 0;
464 m_GGEENumber = 0;
465 m_GG4PiNumber = 0;
466 m_radBhabhaTNumber = 0;
467 m_rhoPiNumber = 0;
468 m_kstarKNumber = 0;
469 m_ppPiPiNumber = 0;
470 m_recoJpsiNumber = 0;
471 m_bhabhaNumber = 0;
472 m_dimuNumber = 0;
473 m_cosmicNumber = 0;
474 m_genProtonNumber = 0;
475 m_psipRhoPiNumber = 0;
476 m_psipKstarKNumber = 0;
477 m_psipPPPiPiNumber = 0;
478
479 log << MSG::INFO << "successfully return from initialize()" << endmsg;
480 return StatusCode::SUCCESS;
481}
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
histoSvc
IHistogramSvc * histoSvc()
Definition ServiceAccessor.h:23

◆ readbeamEfromDb()

void CalibEventSelect::readbeamEfromDb ( int runNo,
double & beamE )

Definition at line 2336 of file CalibEventSelect.cxx.

2336 {
2337
2338 const char host[] = "202.122.37.69";
2339 const char user[] = "guest";
2340 const char passwd[] = "guestpass";
2341 const char db[] = "run";
2342 unsigned int port_num = 3306;
2343
2344 MYSQL* mysql = mysql_init( NULL );
2345 mysql = mysql_real_connect( mysql, host, user, passwd, db, port_num,
2346 NULL, // socket
2347 0 ); // client_flag
2348
2349 if ( mysql == NULL )
2350 { fprintf( stderr, "can not open database: RunInfo for run %d lum\n", runNo ); }
2351
2352 char stmt[1024];
2353
2354 snprintf( stmt, 1024,
2355 "select BER_PRB, BPR_PRB "
2356 "from RunParams where run_number = %d",
2357 runNo );
2358 if ( mysql_real_query( mysql, stmt, strlen( stmt ) ) )
2359 { fprintf( stderr, "query error\n" ); }
2360
2361 MYSQL_RES* result_set = mysql_store_result( mysql );
2362 MYSQL_ROW row = mysql_fetch_row( result_set );
2363 if ( !row )
2364 {
2365 fprintf( stderr, "cannot find data for RunNo %d\n", runNo );
2366 return;
2367 }
2368
2369 double E_E = 0, E_P = 0;
2370 sscanf( row[0], "%lf", &E_E );
2371 sscanf( row[1], "%lf", &E_P );
2372
2373 beamE = ( E_E + E_P ) / 2.0;
2374}
runNo
int runNo
Definition DQA_TO_DB.cxx:13

Referenced by execute().

◆ WhetherSector()

bool CalibEventSelect::WhetherSector ( double ph,
double ph1,
double ph2 )

Definition at line 2303 of file CalibEventSelect.cxx.

2303 {
2304 double phi1 = min( ph1, ph2 );
2305 double phi2 = max( ph1, ph2 );
2306 double delta = 0.0610865; // 3.5*3.1415926/180.
2307 if ( ( phi2 - phi1 ) < CLHEP::pi )
2308 {
2309 phi1 -= delta;
2310 phi2 += delta;
2311 if ( phi1 < 0. ) phi1 += CLHEP::twopi;
2312 if ( phi2 > CLHEP::twopi ) phi2 -= CLHEP::twopi;
2313 double tmp1 = min( phi1, phi2 );
2314 double tmp2 = max( phi1, phi2 );
2315 phi1 = tmp1;
2316 phi2 = tmp2;
2317 }
2318 else
2319 {
2320 phi1 += delta;
2321 phi2 -= delta;
2322 }
2323 if ( ( phi2 - phi1 ) < CLHEP::pi )
2324 {
2325 if ( ph <= phi2 && ph >= phi1 ) return true;
2326 else return false;
2327 }
2328 else
2329 {
2330 if ( ph >= phi2 || ph <= phi1 ) return true;
2331 else return false;
2332 }
2333}
phi2
Double_t phi2
Definition DataBase/tau_mode.c:8
phi1
Double_t phi1
Definition DataBase/tau_mode.c:8
min
#define min(a, b)
Definition Eepipi/src/ee2eepp/basesv5.1/f2c.h:157
max
#define max(a, b)
Definition Eepipi/src/ee2eepp/basesv5.1/f2c.h:158
The documentation for this class was generated from the following files: