BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
Gam4pikp Class Reference

#include <Gam4pikp.h>

+ Inheritance diagram for Gam4pikp:

Public Member Functions

 Gam4pikp (const std::string &name, ISvcLocator *pSvcLocator)
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 

Detailed Description

Definition at line 9 of file Gam4pikp.h.

Constructor & Destructor Documentation

◆ Gam4pikp()

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

Definition at line 62 of file Gam4pikp.cxx.

62 :
63 Algorithm(name, pSvcLocator) {
64
65 //Declare the properties
66 declareProperty("skim4pi", m_skim4pi=false);
67 declareProperty("skim4k", m_skim4k=false);
68 declareProperty("skim2pi2pr", m_skim2pi2pr=false);
69 declareProperty("rootput", m_rootput=false);
70 declareProperty("Ecms", m_ecms=3.6862);
71 declareProperty("Vr0cut", m_vr0cut=1.0);
72 declareProperty("Vz0cut", m_vz0cut=5.0);
73 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
74 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
75 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
76 declareProperty("GammaDangCut", m_gammadangcut=20.0);
77// declareProperty("Test4C", m_test4C = 1);
78// declareProperty("Test5C", m_test5C = 1);
79// declareProperty("CheckDedx", m_checkDedx = 1);
80// declareProperty("CheckTof", m_checkTof = 1);
81 }

Member Function Documentation

◆ execute()

StatusCode Gam4pikp::execute ( )

Definition at line 450 of file Gam4pikp.cxx.

450 {
451
452// std::cout << "execute()" << std::endl;
453 setFilterPassed(false);
454
455
456 MsgStream log(msgSvc(), name());
457 log << MSG::INFO << "in execute()" << endreq;
458
459 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
460 int run=eventHeader->runNumber();
461 int event=eventHeader->eventNumber();
462 log << MSG::DEBUG <<"run, evtnum = "
463 << run << " , "
464 << event <<endreq;
465// cout<<"run "<<run<<endl;
466// cout<<"event "<<event<<endl;
467if(m_rootput)
468{
469 m_run=run;
470 m_rec=event;
471}
472 Ncut[0]++;
473 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
474 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
475 << evtRecEvent->totalCharged() << " , "
476 << evtRecEvent->totalNeutral() << " , "
477 << evtRecEvent->totalTracks() <<endreq;
478
479 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
480
481 if(m_rootput)
482{
483 Gam4pikp::InitVar();
484}
485
486 // for MC topo analysis
487if(m_rootput)
488{
489 if (eventHeader->runNumber()<0)
490 {
491 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
492 int m_numParticle = 0;
493 if (!mcParticleCol)
494 {
495// std::cout << "Could not retrieve McParticelCol" << std::endl;
496 return StatusCode::FAILURE;
497 }
498 else
499 {
500 bool psipDecay = false;
501 int rootIndex = -1;
502 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
503 for (; iter_mc != mcParticleCol->end(); iter_mc++)
504 {
505 if ((*iter_mc)->primaryParticle()) continue;
506 if (!(*iter_mc)->decayFromGenerator()) continue;
507 if ((*iter_mc)->particleProperty()==100443)
508 {
509 psipDecay = true;
510 rootIndex = (*iter_mc)->trackIndex();
511 }
512 if (!psipDecay) continue;
513 int mcidx = ((*iter_mc)->mother()).trackIndex() - rootIndex;
514 int pdgid = (*iter_mc)->particleProperty();
515 m_pdgid[m_numParticle] = pdgid;
516 m_motheridx[m_numParticle] = mcidx;
517 m_numParticle += 1;
518 }
519 }
520 m_idxmc = m_numParticle;
521 }
522}
523
524 Vint iGood, ipip, ipim;
525 iGood.clear();
526 ipip.clear();
527 ipim.clear();
528 Vp4 ppip, ppim;
529 ppip.clear();
530 ppim.clear();
531 int nCharge = 0;
532 int ngdcgx=0;
533 int ncgx=0;
534 Vdouble pTrkCh;
535 pTrkCh.clear();
536 Hep3Vector v(0,0,0);
537 Hep3Vector vv(0,0,0);
538
539 IVertexDbSvc* vtxsvc;
540 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
541 if(vtxsvc->isVertexValid())
542 {
543 double* vertex = vtxsvc->PrimaryVertex();
544 double* vertexsigma = vtxsvc->SigmaPrimaryVertex();
545 v.setX(vertex[0]);
546 v.setY(vertex[1]);
547 v.setZ(vertex[2]);
548 vv.setX(vertexsigma[0]);
549 vv.setY(vertexsigma[1]);
550 vv.setZ(vertexsigma[2]);
551 }
552
553
554 if(run<0)
555 {
556 v[0]=0.;
557 v[1]=0.;
558 v[2]=0.;
559 vv[0]=0.;
560 vv[1]=0.;
561 vv[2]=0.;
562
563 }
564
565
566 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
567 if(i>=evtRecTrkCol->size()) break;
568 ncgx++;
569 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
570 if(!(*itTrk)->isMdcTrackValid()) continue;
571 if (!(*itTrk)->isMdcKalTrackValid()) continue;
572 ngdcgx++;
573 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
574
575 HepVector a = mdcTrk->helix();
576 HepSymMatrix Ea = mdcTrk->err();
577 HepPoint3D pivot(0.,0.,0.);
578 HepPoint3D IP(v[0],v[1],v[2]);
579 VFHelix helixp(pivot,a,Ea);
580 helixp.pivot(IP);
581 HepVector vec = helixp.a();
582 pTrkCh.push_back(mdcTrk->p()*mdcTrk->charge());
583 iGood.push_back((*itTrk)->trackId());
584 nCharge += mdcTrk->charge();
585 }
586
587
588 int nGood = iGood.size();
589 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
590 if((nGood < 4))
591 {
592 return StatusCode::SUCCESS;
593 }
594
595 Ncut[1]++;
596 Vint iGam;
597 int ngamx=0;
598 int ngdgamx=0;
599 iGam.clear();
600 Vint iGamnofit;
601 iGamnofit.clear();
602 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
603 ngamx++;
604 if(i>=evtRecTrkCol->size()) break;
605 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
606 if(!(*itTrk)->isEmcShowerValid()) continue;
607 RecEmcShower *emcTrk = (*itTrk)->emcShower();
608 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
609 // find the nearest charged track
610// double dthe = 200.;
611// double dphi = 200.;
612// double dang = 200.;
613 double dthe = 200.;
614 double dphi = 200.;
615 double dang = 200.;
616
617 ngdgamx++;
618 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
619 if(j>=evtRecTrkCol->size()) break;
620 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
621 if(!(*jtTrk)->isExtTrackValid()) continue;
622 RecExtTrack *extTrk = (*jtTrk)->extTrack();
623 if(extTrk->emcVolumeNumber() == -1) continue;
624 Hep3Vector extpos = extTrk->emcPosition();
625 // double ctht = extpos.cosTheta(emcpos);
626 double angd = extpos.angle(emcpos);
627 double thed = extpos.theta() - emcpos.theta();
628 double phid = extpos.deltaPhi(emcpos);
629 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
630 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
631
632// if(fabs(thed) < fabs(dthe)) dthe = thed;
633// if(fabs(phid) < fabs(dphi)) dphi = phid;
634// if(angd < dang) dang = angd;
635
636 if(angd < dang){ dang=angd;
637 dthe=thed;
638 dphi=phid;
639 }
640
641 }
642 // if(dang>=200) continue;
643 double eraw = emcTrk->energy();
644 dthe = dthe * 180 / (CLHEP::pi);
645 dphi = dphi * 180 / (CLHEP::pi);
646 dang = dang * 180 / (CLHEP::pi);
647// dthert = dthert * 180 / (CLHEP::pi);
648// dphirt = dphirt * 180 / (CLHEP::pi);
649// dangrt = dangrt * 180 / (CLHEP::pi);
650 // if(eraw < m_energyThreshold) continue;
651 iGam.push_back((*itTrk)->trackId());
652 if(eraw < m_energyThreshold) continue;
653 if(dang < 20.0) continue;
654 iGamnofit.push_back((*itTrk)->trackId());
655
656 }
657 // Finish Good Photon Selection
658 //
659 int nGam = iGam.size();
660 int nGamnofit = iGamnofit.size();
661
662 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
663 if(nGam<1){
664 return StatusCode::SUCCESS;
665 }
666 Ncut[2]++;
667
668 if(nGood>20||nGam>30) return StatusCode::SUCCESS;
669
670if(m_rootput)
671{
672m_idxmdc = nGood;
673m_idxemc=nGam;
674}
675//
676 Vp4 pGam;
677 pGam.clear();
678 for(int i = 0; i < nGam; i++) {
679 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
680 RecEmcShower* emcTrk = (*itTrk)->emcShower();
681 double eraw = emcTrk->energy();
682 double phi = emcTrk->phi();
683 double the = emcTrk->theta();
684 HepLorentzVector ptrk;
685 ptrk.setPx(eraw*sin(the)*cos(phi));
686 ptrk.setPy(eraw*sin(the)*sin(phi));
687 ptrk.setPz(eraw*cos(the));
688 ptrk.setE(eraw);
689
690 // ptrk = ptrk.boost(-0.011,0,0);// boost to cms
691
692 pGam.push_back(ptrk);
693 }
694 //
695 // Assign 4-momentum to each charged track
696 //
697 ParticleID *pid = ParticleID::instance();
698
699
700 for(int i = 0; i < nGood; i++) {
701 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
702 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
703
704 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();//After ParticleID, use RecMdcKalTrack substitute RecMdcTrack
705 RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);//PID can set to electron, muon, pion, kaon and proton;The default setting is pion
706
707 if(mdcKalTrk->charge() >0) {
708 ipip.push_back(iGood[i]);
709 HepLorentzVector ptrk;
710 ptrk.setPx(mdcKalTrk->px());
711 ptrk.setPy(mdcKalTrk->py());
712 ptrk.setPz(mdcKalTrk->pz());
713 double p3 = ptrk.mag();
714 ptrk.setE(sqrt(p3*p3+mpi*mpi));
715
716
717 ppip.push_back(ptrk);
718 } else {
719 ipim.push_back(iGood[i]);
720 HepLorentzVector ptrk;
721 ptrk.setPx(mdcKalTrk->px());
722 ptrk.setPy(mdcKalTrk->py());
723 ptrk.setPz(mdcKalTrk->pz());
724 double p3 = ptrk.mag();
725 ptrk.setE(sqrt(p3*p3+mpi*mpi));
726
727 // ptrk = ptrk.boost(-0.011,0,0);//boost to cms
728
729 ppim.push_back(ptrk);
730 }
731 }
732 int npip = ipip.size();
733 int npim = ipim.size();
734 // if(npip*npim != 4) return SUCCESS;
735 if((npip < 2)||(npim < 2)) return SUCCESS;
736
737 Ncut[3]++;
738
739 //
740 // Loop each gamma pair, check ppi0 and pTot
741 //
742
743 int icgp1=-1;
744 int icgp2=-1;
745 int icgm1=-1;
746 int icgm2=-1;
747 int igam=-1;
748 double chisqtrk = 9999.;
749 WTrackParameter wcgp1;
750 WTrackParameter wcgp2;
751 WTrackParameter wcgm1;
752 WTrackParameter wcgm2;
753 int ipip1js=-1;
754 int ipip2js=-1;
755 int ipim1js=-1;
756 int ipim2js=-1;
757 double mcompall=9999;
758 double mppreclst=9999;
759 double meelst=9999;;
760 double mchic2kpilst=9999;
761 double chis4c2kpilst=9999;
762 int type2kpilst=0;
763 double dtpr2kpilst[4]={9999,9999,9999,9999};
764
765 double mc1=mpi,mc2=mpi,mc3=mpi,mc4=mpi;
766 double mchic01=0.0;
767 //maqm HepLorentzVector pgam=(0,0,0,0);
768 HepLorentzVector pgam(0,0,0,0);
769 HepPoint3D xorigin(0.,0.,0.);
770 xorigin[0]=9999.;
771 xorigin[1]=9999.;
772 xorigin[2]=9999.;
773 HepSymMatrix xem(3,0);
774
775 int cl4pi=0;
776 int clajs=0;
777 HepLorentzVector p4psipjs(0.011*m_ecms, 0.0, 0.0, m_ecms);
778 double psipBetajs = (p4psipjs.vect()).mag()/(p4psipjs.e());
779 HepLorentzVector p4psip(0.011*m_ecms, 0.0, 0.0, m_ecms);
780 double psipBeta = (p4psip.vect()).mag()/(p4psip.e()); // beta = P/E
781
782 for(int ii = 0; ii < npip-1; ii++) {
783 RecMdcKalTrack *pip1Trk = (*(evtRecTrkCol->begin()+ipip[ii]))->mdcKalTrack();
784 for(int ij = ii+1; ij < npip; ij++) {
785 RecMdcKalTrack *pip2Trk = (*(evtRecTrkCol->begin()+ipip[ij]))->mdcKalTrack();
786 for(int ik = 0; ik < npim-1; ik++) {
787 RecMdcKalTrack *pim1Trk = (*(evtRecTrkCol->begin()+ipim[ik]))->mdcKalTrack();
788 for(int il = ik+1; il < npim; il++) {
789 RecMdcKalTrack *pim2Trk = (*(evtRecTrkCol->begin()+ipim[il]))->mdcKalTrack();
790 double squar[3]={9999.,9999.,9999.};
791 double squarkpi[6]={9999.,9999.,9999.,9999.,9999.,9999.};
792 WTrackParameter wvpip1Trk, wvpim1Trk, wvpip2Trk, wvpim2Trk;
793
794 Vint iGoodjs;
795 iGoodjs.clear();
796 Vdouble pTrkjs;
797 pTrkjs.clear();
798
799
800
801 pTrkjs.push_back(pip1Trk->p()*pip1Trk->charge());
802 pTrkjs.push_back(pip2Trk->p()*pip2Trk->charge());
803 pTrkjs.push_back(pim1Trk->p()*pim1Trk->charge());
804 pTrkjs.push_back(pim2Trk->p()*pim2Trk->charge());
805 iGoodjs.push_back(ipip[ii]);
806 iGoodjs.push_back(ipip[ij]);
807 iGoodjs.push_back(ipim[ik]);
808 iGoodjs.push_back(ipim[il]);
809
810 Gam4pikp::BubbleSort(pTrkjs, iGoodjs);
811 Vint jGoodjs;
812 jGoodjs.clear();
813 Vp4 p4chTrkjs;
814 p4chTrkjs.clear();
815 Vint i4cpip1js, i4cpip2js, i4cpim1js, i4cpim2js;
816 i4cpip1js.clear();
817 i4cpip2js.clear();
818 i4cpim1js.clear();
819 i4cpim2js.clear();
820 i4cpip1js.push_back(iGoodjs[2]);
821 i4cpip2js.push_back(iGoodjs[3]);
822 i4cpim1js.push_back(iGoodjs[1]);
823 i4cpim2js.push_back(iGoodjs[0]);
824 jGoodjs.push_back(i4cpip1js[0]);
825 jGoodjs.push_back(i4cpim1js[0]);
826 jGoodjs.push_back(i4cpip2js[0]);
827 jGoodjs.push_back(i4cpim2js[0]);
828
829 for (int i = 0; i < 4; i++)
830 {
831 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGoodjs[i];
832 if (!(*itTrk)->isMdcTrackValid()) continue;
833 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
834 if (!(*itTrk)->isMdcKalTrackValid()) continue;
835 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
836 double ptrk;
837 HepLorentzVector p4trk;
838 if (i<2)
839 {
840 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
841 ptrk = mdcKalTrk->p();
842 p4trk.setPx(mdcKalTrk->px());
843 p4trk.setPy(mdcKalTrk->py());
844 p4trk.setPz(mdcKalTrk->pz());
845 p4trk.setE(sqrt(ptrk*ptrk+mpi*mpi));
846 }
847 else{
848 mdcKalTrk->setPidType(RecMdcKalTrack::electron);
849 ptrk = mdcKalTrk->p();
850 p4trk.setPx(mdcKalTrk->px());
851 p4trk.setPy(mdcKalTrk->py());
852 p4trk.setPz(mdcKalTrk->pz());
853 p4trk.setE(sqrt(ptrk*ptrk+xmass[0]*xmass[0]));
854 }
855 p4trk.boost(-1.0*psipBetajs, 0.0, 0.0);
856 p4chTrkjs.push_back(p4trk);
857 }
858 p4psipjs.boost(-1.0*psipBetajs, 0.0, 0.0);
859
860 HepLorentzVector p4pipijs = p4chTrkjs[0] + p4chTrkjs[1];
861 HepLorentzVector p4eejs = p4chTrkjs[2] + p4chTrkjs[3];
862 HepLorentzVector p4pipirecjs = p4psipjs - p4pipijs;
863 double mpprecjs=p4pipirecjs.m();
864 double mpipijs=p4pipijs.m();
865 double meejs=p4eejs.m();
866 double mcomp=sqrt((mpprecjs-3.097)*(mpprecjs-3.097)+(meejs-3.097)*(meejs-3.097));
867 if(mcomp<mcompall)
868 {
869 mcompall=mcomp;
870 ipip1js=i4cpip1js[0];
871 ipim1js=i4cpim1js[0];
872 ipip2js=i4cpip2js[0];
873 ipim2js=i4cpim2js[0];
874 mppreclst=mpprecjs;
875 meelst=meejs;
876 }
877
878if(m_rootput)
879{
880 m_mpprecall=mppreclst;
881 m_meeall=meelst;
882}
883 }
884
885 }
886 }
887 }
888
889 //{
890// HepLorentzVector p4psip(0.011*m_ecms, 0.0, 0.0, m_ecms);
891// double psipBeta = (p4psip.vect()).mag()/(p4psip.e()); // beta = P/E
892 Vint iGood4c;
893 iGood4c.clear();
894 Vdouble pTrk4c;
895 pTrk4c.clear();
896 Vint jGood;
897 jGood.clear();
898 Vint jsGood;
899 jsGood.clear();
900
901 if(mcompall>9997)
902 return StatusCode::SUCCESS;
903
904 jsGood.push_back(ipip1js);
905 jsGood.push_back(ipim1js);
906 jsGood.push_back(ipip2js);
907 jsGood.push_back(ipim2js);
908
909 for(int i = 0; i < evtRecEvent->totalCharged(); i++)
910 {
911 if(i>=evtRecTrkCol->size()) break;
912 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
913 if(!(*itTrk)->isMdcTrackValid()) continue;
914 if (!(*itTrk)->isMdcKalTrackValid()) continue;
915 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
916 if((i!=ipip1js)&&(i!=ipim1js)&&(i!=ipip2js)&&(i!=ipim2js))
917 {
918 jsGood.push_back(i);
919 }
920 }
921
922 int njsGood=jsGood.size();
923 //ParticleID *pid = ParticleID::instance();
924if(m_rootput)
925{
926 for (int i = 0; i < njsGood; i++)
927 {
928 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jsGood[i];
929 if (!(*itTrk)->isMdcTrackValid()) continue;
930 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
931 if (!(*itTrk)->isMdcKalTrackValid()) continue;
932 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
933 double ptrk;
934 ptrk = mdcKalTrk->p();
935 m_x0js[i] = mdcTrk->x();
936 m_y0js[i] = mdcTrk->y();
937 m_z0js[i] = mdcTrk->z();
938 m_r0js[i] = mdcTrk->r();
939 m_ppmdcjs[i] = mdcTrk->p();
940 m_pxmdcjs[i] = mdcTrk->px();
941 m_pymdcjs[i] = mdcTrk->py();
942 m_pzmdcjs[i] = mdcTrk->pz();
943 m_ppkaljs[i] = mdcKalTrk->p();
944 Hep3Vector p3jsi(mdcTrk->px(),mdcTrk->py(),mdcTrk->pz());
945 if(njsGood>4){
946 EvtRecTrackIterator itTrk5 = evtRecTrkCol->begin() + jsGood[4];
947 RecMdcTrack *mdcTrk5 = (*itTrk5)->mdcTrack();
948 Hep3Vector p3js5(mdcTrk5->px(),mdcTrk5->py(),mdcTrk5->pz());
949 m_angjs5[i]=p3jsi.angle(p3js5);
950 m_nearjs5[i]=p3jsi.howNear(p3js5);
951 }
952 if(njsGood>5){
953 EvtRecTrackIterator itTrk6 = evtRecTrkCol->begin() + jsGood[5];
954 RecMdcTrack *mdcTrk6 = (*itTrk6)->mdcTrack();
955 Hep3Vector p3js6(mdcTrk6->px(),mdcTrk6->py(),mdcTrk6->pz());
956 m_angjs6[i]=p3jsi.angle(p3js6);
957 m_nearjs6[i]=p3jsi.howNear(p3js6);
958 }
959
960
961 m_ptmdcjs[i] = mdcTrk->pxy();
962 m_ptkaljs[i] = mdcKalTrk->pxy();
963 double x0=mdcTrk->x();
964 double y0=mdcTrk->y();
965 double z0=mdcTrk->z();
966 double phi0=mdcTrk->helix(1);
967 double xv=v[0];
968 double yv=v[1];
969 double zv=v[2];
970 double Rxy=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
971 double Rz=z0-zv;
972 m_Rxyjs[i]=Rxy;
973 m_Rzjs[i]=Rz;
974 HepVector a = mdcTrk->helix();
975 HepSymMatrix Ea = mdcTrk->err();
976 HepPoint3D pivot(0.,0.,0.);
977 HepPoint3D IP(v[0],v[1],v[2]);
978 VFHelix helixp(pivot,a,Ea);
979 helixp.pivot(IP);
980 HepVector vec = helixp.a();
981 m_Rnxyjs[i]=vec[0];
982 m_Rnzjs[i]=vec[3];
983 m_phinjs[i]=vec[1];
984 //tof
985 if ((*itTrk)->isTofTrackValid())
986 { m_bjtofvaljs[i]=1;
987 m_cotof1js[i]=1;
988 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
989 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
990 for (;iter_tof != tofTrkCol.end(); iter_tof++ ){
991 TofHitStatus *status = new TofHitStatus;
992 status->setStatus((*iter_tof)->status());
993 if(status->is_cluster()){
994 m_bjtofjs[i]=1;
995 m_counterjs[i] = status->is_counter();
996 m_barreljs[i] = status->is_barrel();
997 m_layertofjs[i] = status->layer();
998 m_readoutjs[i] = status->is_readout();
999 m_clusterjs[i] = status->is_cluster();
1000 m_cotof2js[i]=2;
1001 m_betajs[i]= (*iter_tof)->beta();
1002 m_tofjs[i] = (*iter_tof)->tof();
1003 m_tofpathjs[i] = (*iter_tof)->path();
1004 m_zhitjs[i]= (*iter_tof)->zrhit();
1005 m_texejs[i] = (*iter_tof)->texpElectron();
1006 m_texmujs[i] = (*iter_tof)->texpMuon();
1007 m_texpijs[i] = (*iter_tof)->texpPion();
1008 m_texkjs[i] = (*iter_tof)->texpKaon();
1009 m_texprjs[i] = (*iter_tof)->texpProton();
1010 m_dtejs[i] = m_tofjs[i] - m_texejs[i];
1011 m_dtmujs[i] = m_tofjs[i] - m_texmujs[i];
1012 m_dtpijs[i] = m_tofjs[i] - m_texpijs[i];
1013 m_dtkjs[i] = m_tofjs[i] - m_texkjs[i];
1014 m_dtprjs[i] = m_tofjs[i] - m_texprjs[i];
1015
1016 m_sigmaetofjs[i]=(*iter_tof)->sigma(0);
1017 m_sigmamutofjs[i]=(*iter_tof)->sigma(1);
1018 m_sigmapitofjs[i]=(*iter_tof)->sigma(2);
1019 m_sigmaktofjs[i]=(*iter_tof)->sigma(3);
1020 m_sigmaprtofjs[i]=(*iter_tof)->sigma(4);
1021 m_t0tofjs[i]=(*iter_tof)->t0();
1022 m_errt0tofjs[i]=(*iter_tof)->errt0();
1023
1024 m_tofIDjs[i]=(*iter_tof)->tofID();
1025 m_clusterIDjs[i]=(*iter_tof)->tofTrackID();
1026 }
1027 delete status;
1028 }
1029 }
1030 //dedx
1031 if ((*itTrk)->isMdcDedxValid())
1032 {
1033
1034 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
1035 m_chiejs[i] = dedxTrk->chiE();
1036 m_chimujs[i] = dedxTrk->chiMu();
1037 m_chipijs[i] = dedxTrk->chiPi();
1038 m_chikjs[i] = dedxTrk->chiK();
1039 m_chipjs[i] = dedxTrk->chiP();
1040 m_ghitjs[i] = dedxTrk->numGoodHits();
1041 m_thitjs[i] = dedxTrk->numTotalHits();
1042 m_probphjs[i] = dedxTrk->probPH();
1043 m_normphjs[i] = dedxTrk->normPH();
1044 }
1045
1046
1047
1048 //emc
1049 if ( (*itTrk)->isEmcShowerValid() )
1050 {
1051 RecEmcShower* emcTrk = (*itTrk)->emcShower();
1052 m_bjemcjs[i]=1;
1053 m_emcjs[i] = emcTrk->energy();
1054 m_evpjs[i] = emcTrk->energy()/ptrk;
1055 m_timecgjs[i] = emcTrk->time();
1056 // totalEnergy += emcTrk->energy();
1057 }
1058
1059 //muc
1060 if ( (*itTrk)->isMucTrackValid() )
1061 {
1062 RecMucTrack* mucTrk = (*itTrk)->mucTrack();
1063 double dpthp = mucTrk->depth()/25.0;//why?
1064 // m_depthmucjs[i] = dpthp<10.0 ? dpthp : 10.0;
1065 m_bjmucjs[i]=1;
1066 m_depthmucjs[i]=mucTrk->depth();
1067 m_layermucjs[i] = mucTrk->numLayers();
1068 }
1069
1070 m_cosmdcjs[i] = cos(mdcTrk->theta());
1071 m_phimdcjs[i] = mdcTrk->phi();
1072 //pid
1073 pid->init();
1074 pid->setMethod(pid->methodProbability());
1075 pid->setChiMinCut(4);
1076 pid->setRecTrack(*itTrk);
1077 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
1078 pid->identify(pid->onlyPion() | pid->onlyKaon() | pid->onlyProton());
1079 pid->identify(pid->onlyElectron() | pid->onlyMuon());
1080 pid->calculate();
1081 if(!(pid->IsPidInfoValid())) continue;
1082 // m_costpidjs[i] = cos(mdcTrk->theta());
1083 m_dedxpidjs[i] = pid->chiDedx(2);
1084 m_tof1pidjs[i] = pid->chiTof1(2);
1085 m_tof2pidjs[i] = pid->chiTof2(2);
1086 m_probejs[i] = pid->probElectron();
1087 m_probmujs[i] = pid->probMuon();
1088 m_probpijs[i] = pid->probPion();
1089 m_probkjs[i] = pid->probKaon();
1090 m_probprjs[i] = pid->probProton();
1091
1092
1093
1094 }
1095
1096}
1097
1098 Vint jGam2kpi;
1099 jGam2kpi.clear();
1100
1101 Vint iGood2kpi, ipip2kpi, ipim2kpi;
1102 iGood2kpi.clear();
1103 ipip2kpi.clear();
1104 ipim2kpi.clear();
1105 Vp4 ppip2kpi, ppim2kpi;
1106 ppip2kpi.clear();
1107 ppim2kpi.clear();
1108
1109 Vint ipipnofit, ipimnofit, ikpnofit, ikmnofit, ipropnofit, ipromnofit;
1110 ipipnofit.clear();
1111 ipimnofit.clear();
1112 ikpnofit.clear();
1113 ikmnofit.clear();
1114 ipropnofit.clear();
1115 ipromnofit.clear();
1116 Vp4 ppipnofit, ppimnofit, pkpnofit, pkmnofit, ppropnofit, ppromnofit;
1117 ppipnofit.clear();
1118 ppimnofit.clear();
1119 pkpnofit.clear();
1120 pkmnofit.clear();
1121 ppropnofit.clear();
1122 ppromnofit.clear();
1123
1124 Vdouble p3pip2kpi;
1125 p3pip2kpi.clear();
1126 Vdouble p3pim2kpi;
1127 p3pim2kpi.clear();
1128
1129 Vint itrak2kpi;
1130 itrak2kpi.clear();
1131 int cls2kpi;
1132 double chis2kpi=9999.;
1133 double m4piall=9999.;
1134 double m4kall=9999.;
1135 double mgam4piall=9999.;
1136 double mgam4kall=9999.;
1137 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
1138 if(i>=evtRecTrkCol->size()) break;
1139 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
1140 if(!(*itTrk)->isMdcTrackValid()) continue;
1141 if (!(*itTrk)->isMdcKalTrackValid()) continue;
1142 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
1143 double z02kpi = mdcTrk->z();
1144 double r02kpi = mdcTrk->r();
1145 HepVector a = mdcTrk->helix();
1146 HepSymMatrix Ea = mdcTrk->err();
1147 HepPoint3D pivot(0.,0.,0.);
1148 HepPoint3D IP(v[0],v[1],v[2]);
1149 VFHelix helixp(pivot,a,Ea);
1150 helixp.pivot(IP);
1151 HepVector vec = helixp.a();
1152 double Rnxy02kpi=vec[0];
1153 double Rnz02kpi=vec[3];
1154 // if(fabs(Rnxy02kpi>1.0)) continue;
1155 // if(fabs(Rnz02kpi>10.0)) continue;
1156 iGood2kpi.push_back((*itTrk)->trackId());
1157 }
1158 int nGood2kpi = iGood2kpi.size();
1159 if(nGood2kpi==4)
1160 {
1161 for(int i = 0; i < nGood2kpi; i++) {
1162 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood2kpi[i];
1163 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
1164 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
1165 if(mdcKalTrk->charge() >0) {
1166 ipip2kpi.push_back(iGood2kpi[i]);
1167 p3pip2kpi.push_back(mdcKalTrk->p());
1168 HepLorentzVector ptrk;
1169 ptrk.setPx(mdcKalTrk->px());
1170 ptrk.setPy(mdcKalTrk->py());
1171 ptrk.setPz(mdcKalTrk->pz());
1172 double p3 = ptrk.mag();
1173 ptrk.setE(sqrt(p3*p3+mpi*mpi));
1174 ppip2kpi.push_back(ptrk);
1175 } else {
1176 ipim2kpi.push_back(iGood2kpi[i]);
1177 p3pim2kpi.push_back(mdcKalTrk->p());
1178 HepLorentzVector ptrk;
1179 ptrk.setPx(mdcKalTrk->px());
1180 ptrk.setPy(mdcKalTrk->py());
1181 ptrk.setPz(mdcKalTrk->pz());
1182 double p3 = ptrk.mag();
1183 ptrk.setE(sqrt(p3*p3+mpi*mpi));
1184 ppim2kpi.push_back(ptrk);
1185 }
1186 }
1187 int npip2kpi = ipip2kpi.size();
1188 int npim2kpi = ipim2kpi.size();
1189
1190
1191
1192 ParticleID *pid = ParticleID::instance();
1193
1194if(m_rootput)
1195{
1196 m_cy2kpi=6;
1197}
1198
1199 if((npip2kpi == 2)&&(npim2kpi == 2))
1200 {
1201 //if(nGamnofit>=1)
1202
1203
1204Ncut[4]++;
1205
1206 HepPoint3D xorigin2kpi(0.,0.,0.);
1207 xorigin2kpi[0]=9999.;
1208 xorigin2kpi[1]=9999.;
1209 xorigin2kpi[2]=9999.;
1210 HepSymMatrix xem2kpi(3,0);
1211
1212 Gam4pikp::BubbleSort(p3pip2kpi, ipip2kpi);
1213 Gam4pikp::BubbleSort(p3pim2kpi, ipim2kpi);
1214
1215 Mcut[1]++;
1216 RecMdcKalTrack *pip12kpiTrk = (*(evtRecTrkCol->begin()+ipip2kpi[0]))->mdcKalTrack();
1217 RecMdcKalTrack *pip22kpiTrk = (*(evtRecTrkCol->begin()+ipip2kpi[1]))->mdcKalTrack();
1218 RecMdcKalTrack *pim12kpiTrk = (*(evtRecTrkCol->begin()+ipim2kpi[0]))->mdcKalTrack();
1219 RecMdcKalTrack *pim22kpiTrk = (*(evtRecTrkCol->begin()+ipim2kpi[1]))->mdcKalTrack();
1220 double squar2kpi[10]={9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.};
1221 double mc12kpi,mc22kpi,mc32kpi,mc42kpi;
1222 // double chis2kpi=9999.;
1223 WTrackParameter wcgp12kpi;
1224 WTrackParameter wcgp22kpi;
1225 WTrackParameter wcgm12kpi;
1226 WTrackParameter wcgm22kpi;
1227 int icgp12kpi=-1;
1228 int icgp22kpi=-1;
1229 int icgm12kpi=-1;
1230 int icgm22kpi=-1;
1231 int igam2kpi=-1;
1232 double m2kpi=9999;
1233
1234 int n20=0;
1235 int n30=0;
1236 WTrackParameter wvpip12kpiTrk, wvpim12kpiTrk, wvpip22kpiTrk, wvpim22kpiTrk;
1237 for(int k=0;k<6;k++)
1238 {
1239 if(k==0){mc12kpi=mpi;mc22kpi=mpi;mc32kpi=mpi;mc42kpi=mpi;}
1240 if(k==1){mc12kpi=mk;mc22kpi=mk;mc32kpi=mk;mc42kpi=mk;}
1241 if(k==2){mc12kpi=mpi;mc22kpi=mpro;mc32kpi=mpi;mc42kpi=mpro;}
1242 if(k==3){mc12kpi=mpro;mc22kpi=mpi;mc32kpi=mpro;mc42kpi=mpi;}
1243 if(k==4){mc12kpi=mpi;mc22kpi=mpro;mc32kpi=mpro;mc42kpi=mpi; }
1244 if(k==5){mc12kpi=mpro;mc22kpi=mpi;mc32kpi=mpi;mc42kpi=mpro;}
1245
1246
1247
1248 wvpip12kpiTrk = WTrackParameter(mc12kpi, pip12kpiTrk->getZHelix(), pip12kpiTrk->getZError());
1249 wvpip22kpiTrk = WTrackParameter(mc22kpi, pip22kpiTrk->getZHelix(), pip22kpiTrk->getZError());
1250 wvpim12kpiTrk = WTrackParameter(mc32kpi, pim12kpiTrk->getZHelix(), pim12kpiTrk->getZError());
1251 wvpim22kpiTrk = WTrackParameter(mc42kpi, pim22kpiTrk->getZHelix(), pim22kpiTrk->getZError());
1252 HepPoint3D vx(0., 0., 0.);
1253 HepSymMatrix Evx(3, 0);
1254 double bx = 1E+6;
1255 double by = 1E+6;
1256 double bz = 1E+6;
1257 Evx[0][0] = bx*bx;
1258 Evx[1][1] = by*by;
1259 Evx[2][2] = bz*bz;
1260 VertexParameter vxpar;
1261 vxpar.setVx(vx);
1262 vxpar.setEvx(Evx);
1263 VertexFit* vtxfit = VertexFit::instance();
1264 vtxfit->init();
1265 vtxfit->AddTrack(0, wvpip12kpiTrk);
1266 vtxfit->AddTrack(1, wvpim12kpiTrk);
1267 vtxfit->AddTrack(2, wvpip22kpiTrk);
1268 vtxfit->AddTrack(3, wvpim22kpiTrk);
1269 vtxfit->AddVertex(0, vxpar,0, 1, 2, 3);
1270 if(!vtxfit->Fit(0)) continue;
1271 vtxfit->Swim(0);
1272 WTrackParameter wpip12kpi = vtxfit->wtrk(0);
1273 WTrackParameter wpim12kpi = vtxfit->wtrk(1);
1274 WTrackParameter wpip22kpi = vtxfit->wtrk(2);
1275 WTrackParameter wpim22kpi = vtxfit->wtrk(3);
1276 WTrackParameter wpip12kpiys = vtxfit->wtrk(0);
1277 WTrackParameter wpim12kpiys = vtxfit->wtrk(1);
1278 WTrackParameter wpip22kpiys = vtxfit->wtrk(2);
1279 WTrackParameter wpim22kpiys = vtxfit->wtrk(3);
1280 xorigin2kpi = vtxfit->vx(0);
1281 xem2kpi = vtxfit->Evx(0);
1282 KalmanKinematicFit * kmfit = KalmanKinematicFit::instance();
1283
1284 HepLorentzVector ecms(0.040547,0,0,3.68632);
1285
1286 double chisq = 9999.;
1287 int ig1 = -1;
1288 for(int i = 0; i < nGam; i++) {
1289 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
1290 kmfit->init();
1291 kmfit->setBeamPosition(xorigin2kpi);
1292 kmfit->setVBeamPosition(xem2kpi);
1293 kmfit->AddTrack(0, wpip12kpi);
1294 kmfit->AddTrack(1, wpim12kpi);
1295 kmfit->AddTrack(2, wpip22kpi);
1296 kmfit->AddTrack(3, wpim22kpi);
1297 kmfit->AddTrack(4, 0.0, g1Trk);
1298 kmfit->AddFourMomentum(0, p4psip);
1299 // if(!kmfit->Fit(0)) continue;
1300 bool oksq = kmfit->Fit();
1301 if(oksq) {
1302 double chi2 = kmfit->chisq();
1303 if(chi2 < chisq) {
1304 chisq = chi2;
1305 squar2kpi[k]=chi2;
1306 ig1 = iGam[i];
1307
1308 }
1309 }
1310 }
1311 if(squar2kpi[k]<200&&squar2kpi[k]<chis2kpi)
1312 {
1313// m_comcs2kpi[k]=squar2kpi[k];
1314 chis2kpi=squar2kpi[k];
1315 if(squar2kpi[k]<20) n20=n20+1;
1316 if(squar2kpi[k]<30) n30=n30+1;
1317
1318 icgp12kpi=ipip2kpi[0];
1319 icgp22kpi=ipip2kpi[1];
1320 icgm12kpi=ipim2kpi[0];
1321 icgm22kpi=ipim2kpi[1];
1322 igam2kpi=ig1;
1323 wcgp12kpi=wpip12kpiys;
1324 wcgp22kpi=wpip22kpiys;
1325 wcgm12kpi=wpim12kpiys;
1326 wcgm22kpi=wpim22kpiys;
1327 cls2kpi=k;
1328
1329 if(k==0){
1330 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
1331 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
1332 }
1333
1334 if(k==1){
1335 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
1336 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
1337 }
1338
1339
1340 if(k==2){
1341 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
1342 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
1343 }
1344
1345 if(k==3){
1346 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
1347 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
1348 }
1349
1350 if(k==4){
1351 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[1]);
1352 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[0]);
1353 }
1354
1355 if(k==5){
1356 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[0]);
1357 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[1]);
1358 }
1359
1360
1361
1362
1363
1364 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+igam2kpi))->emcShower();
1365 kmfit->init();
1366 kmfit->setBeamPosition(xorigin2kpi);
1367 kmfit->setVBeamPosition(xem2kpi);
1368 kmfit->AddTrack(0, wpip12kpi);
1369 kmfit->AddTrack(1, wpim12kpi);
1370 kmfit->AddTrack(2, wpip22kpi);
1371 kmfit->AddTrack(3, wpim22kpi);
1372 kmfit->AddTrack(4, 0.0, g1Trk);
1373 kmfit->AddFourMomentum(0, p4psip);
1374 bool oksq = kmfit->Fit();
1375 if(oksq){
1376 HepLorentzVector pchic2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3);
1377 HepLorentzVector ppsip2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3) + kmfit->pfit(4);
1378 mchic2kpilst=pchic2kpi.m();
1379 chis4c2kpilst=kmfit->chisq();
1380if(m_rootput)
1381{
1382 m_mchic2kpi = pchic2kpi.m();
1383 m_chis2kpi = kmfit->chisq();
1384 m_mpsip2kpi=ppsip2kpi.m();
1385}
1386
1387
1388
1389 }
1390
1391
1392 }
1393 }
1394
1395
1396 if(chis2kpi<200)
1397 {
1398Ncut[5]++;
1399 if(m_rootput)
1400 {
1401 m_ncy20=n20;
1402 m_ncy30=n30;
1403 m_cla2kpi=cls2kpi;
1404 }
1405 type2kpilst=cls2kpi;
1406
1407 Vp4 p2kpi;
1408 p2kpi.clear();
1409
1410 for (int i = 0; i < 4; i++)
1411 {
1412 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + itrak2kpi[i];
1413 if (!(*itTrk)->isMdcTrackValid()) continue;
1414 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
1415 if (!(*itTrk)->isMdcKalTrackValid()) continue;
1416 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
1417 double ptrk2kpi;
1418 HepLorentzVector p4trk2kpi;
1419 if(cls2kpi==1)
1420 {
1421 mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
1422 ptrk2kpi = mdcKalTrk->p();
1423 p4trk2kpi.setPx(mdcKalTrk->px());
1424 p4trk2kpi.setPy(mdcKalTrk->py());
1425 p4trk2kpi.setPz(mdcKalTrk->pz());
1426 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mk*mk));
1427 p2kpi.push_back(p4trk2kpi);
1428 }
1429
1430 if(cls2kpi==2)
1431 {
1432 if (i<2)
1433 {
1434 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
1435 ptrk2kpi = mdcKalTrk->p();
1436 p4trk2kpi.setPx(mdcKalTrk->px());
1437 p4trk2kpi.setPy(mdcKalTrk->py());
1438 p4trk2kpi.setPz(mdcKalTrk->pz());
1439 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mpi*mpi));
1440 p2kpi.push_back(p4trk2kpi);
1441
1442 }
1443 else{
1444 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
1445 ptrk2kpi = mdcKalTrk->p();
1446 p4trk2kpi.setPx(mdcKalTrk->px());
1447 p4trk2kpi.setPy(mdcKalTrk->py());
1448 p4trk2kpi.setPz(mdcKalTrk->pz());
1449 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mpro*mpro));
1450 p2kpi.push_back(p4trk2kpi);
1451
1452 }
1453
1454 }
1455 if(cls2kpi!=1&&cls2kpi!=2)
1456 {
1457 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
1458 ptrk2kpi = mdcKalTrk->p();
1459 p4trk2kpi.setPx(mdcKalTrk->px());
1460 p4trk2kpi.setPy(mdcKalTrk->py());
1461 p4trk2kpi.setPz(mdcKalTrk->pz());
1462 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mpi*mpi));
1463 p2kpi.push_back(p4trk2kpi);
1464 }
1465
1466
1467
1468
1469 }
1470
1471 for (int i = 0; i < 4; i++)
1472 {
1473 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + itrak2kpi[i];
1474 if (!(*itTrk)->isMdcTrackValid()) continue;
1475 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
1476 if (!(*itTrk)->isMdcKalTrackValid()) continue;
1477 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
1478 if ((*itTrk)->isTofTrackValid())
1479 {
1480 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
1481 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1482 for (;iter_tof != tofTrkCol.end(); iter_tof++ ){
1483 TofHitStatus *status = new TofHitStatus;
1484 status->setStatus((*iter_tof)->status());
1485 if(status->is_cluster()){ dtpr2kpilst[i]=((*iter_tof)->tof()-(*iter_tof)->texpProton());
1486 }
1487 delete status;
1488 }
1489 }
1490 }
1491
1492
1493 /*
1494 HepLorentzVector ecmsb(0.040547,0,0,3.68632);
1495 double chisq = 9999.;
1496 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+igam2kpi))->emcShower();
1497 KalmanKinematicFit * kmfit = KalmanKinematicFit::instance();
1498 kmfit->init();
1499 kmfit->setBeamPosition(xorigin2kpi);
1500 kmfit->setVBeamPosition(xem2kpi);
1501 kmfit->AddTrack(0, wcgp12kpi);
1502 kmfit->AddTrack(1, wcgp22kpi);
1503 kmfit->AddTrack(2, wcgm12kpi);
1504 kmfit->AddTrack(3, wcgm22kpi);
1505 kmfit->AddTrack(4, 0.0, g1Trk);
1506 kmfit->AddFourMomentum(0, p4psip);
1507 bool oksq = kmfit->Fit();
1508 if(oksq) {
1509 Mcut[3]++;
1510 HepLorentzVector pchic4c2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3);
1511 HepLorentzVector ppsip4c2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3) + kmfit->pfit(4);
1512
1513 m_mchic4c2kpi = pchic4c2kpi.m();
1514// m2kpi=m_mchic4c2kpi;
1515 m_chis4c2kpi = kmfit->chisq();
1516 m_mpsip4c2kpi=ppsip4c2kpi.m();
1517
1518
1519
1520
1521 }
1522
1523*/
1524 Mcut[2]++;
1525if(m_rootput)
1526{
1527 for (int i = 0; i < 4; i++)
1528 {
1529 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + itrak2kpi[i];
1530 if (!(*itTrk)->isMdcTrackValid()) continue;
1531 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
1532 if (!(*itTrk)->isMdcKalTrackValid()) continue;
1533 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
1534 m_ppmdc2kpi[i]=mdcTrk->p();
1535 m_pxmdc2kpi[i]=mdcTrk->px();
1536 m_pymdc2kpi[i]=mdcTrk->py();
1537 m_pzmdc2kpi[i]=mdcTrk->pz();
1538 m_ppkal2kpi[i]=mdcKalTrk->p();
1539 m_charge2kpi[i]=mdcKalTrk->charge();
1540 double ptrk;
1541 ptrk=mdcKalTrk->p();
1542
1543 if (eventHeader->runNumber()<0)
1544 {double mcall=9999;
1545 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
1546 int m_numParticle = 0;
1547 if (!mcParticleCol)
1548 {
1549// std::cout << "Could not retrieve McParticelCol" << std::endl;
1550 return StatusCode::FAILURE;
1551 }
1552 else
1553 {
1554 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
1555 for (; iter_mc != mcParticleCol->end(); iter_mc++)
1556 { double commc;
1557 int pdgcode = (*iter_mc)->particleProperty();
1558 float px,py,pz;
1559 px=(*iter_mc)->initialFourMomentum().x();
1560 py=(*iter_mc)->initialFourMomentum().y();
1561 pz=(*iter_mc)->initialFourMomentum().z();
1562
1563 commc=ptrk*ptrk-px*px-py*py-pz*pz;
1564 if(fabs(commc)<fabs(mcall))
1565 {
1566 mcall=commc;
1567 m_pdg[i]=pdgcode;
1568 m_cbmc[i]=commc;
1569 }
1570
1571 }
1572
1573 }
1574
1575 }
1576
1577 if ((*itTrk)->isTofTrackValid())
1578 { m_bjtofval2kpi[i]=1;
1579 m_cotof12kpi[i]=1;
1580 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
1581 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1582 for (;iter_tof != tofTrkCol.end(); iter_tof++ ){
1583 TofHitStatus *status = new TofHitStatus;
1584 status->setStatus((*iter_tof)->status());
1585
1586 if(status->is_cluster()){
1587
1588
1589 m_bjtof2kpi[i]=1;
1590 m_counter2kpi[i] = status->is_counter();
1591 m_barrel2kpi[i] = status->is_barrel();
1592 m_layertof2kpi[i] = status->layer();
1593 m_readout2kpi[i] = status->is_readout();
1594 m_cluster2kpi[i] = status->is_cluster();
1595 m_cotof22kpi[i]=2;
1596 m_beta2kpi[i]= (*iter_tof)->beta();
1597 m_tof2kpi[i] = (*iter_tof)->tof();
1598 m_tofpath2kpi[i] = (*iter_tof)->path();
1599 m_zhit2kpi[i]= (*iter_tof)->zrhit();
1600 m_texe2kpi[i] = (*iter_tof)->texpElectron();
1601 m_texmu2kpi[i] = (*iter_tof)->texpMuon();
1602 m_texpi2kpi[i] = (*iter_tof)->texpPion();
1603 m_texk2kpi[i] = (*iter_tof)->texpKaon();
1604 m_texpr2kpi[i] = (*iter_tof)->texpProton();
1605 m_dte2kpi[i] = m_tof2kpi[i] - m_texe2kpi[i];
1606 m_dtmu2kpi[i] = m_tof2kpi[i] - m_texmu2kpi[i];
1607 m_dtpi2kpi[i] = m_tof2kpi[i] - m_texpi2kpi[i];
1608 m_dtk2kpi[i] = m_tof2kpi[i] - m_texk2kpi[i];
1609 m_dtpr2kpi[i] = m_tof2kpi[i] - m_texpr2kpi[i];
1610 m_tofID2kpi[i]=(*iter_tof)->tofID();
1611 m_clusterID2kpi[i]=(*iter_tof)->tofTrackID();
1612 m_sigmaetof2kpi[i]=(*iter_tof)->sigma(0);
1613 m_sigmamutof2kpi[i]=(*iter_tof)->sigma(1);
1614 m_sigmapitof2kpi[i]=(*iter_tof)->sigma(2);
1615 m_sigmaktof2kpi[i]=(*iter_tof)->sigma(3);
1616 m_sigmaprtof2kpi[i]=(*iter_tof)->sigma(4);
1617 m_t0tof2kpi[i]=(*iter_tof)->t0();
1618 m_errt0tof2kpi[i]=(*iter_tof)->errt0();
1619
1620 }
1621 delete status;
1622 }
1623 }
1624
1625
1626
1627
1628
1629
1630 if ((*itTrk)->isMdcDedxValid())
1631 {
1632 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
1633 m_chie2kpi[i] = dedxTrk->chiE();
1634 m_chimu2kpi[i] = dedxTrk->chiMu();
1635 m_chipi2kpi[i] = dedxTrk->chiPi();
1636 m_chik2kpi[i] = dedxTrk->chiK();
1637 m_chip2kpi[i] = dedxTrk->chiP();
1638 m_ghit2kpi[i] = dedxTrk->numGoodHits();
1639 m_thit2kpi[i] = dedxTrk->numTotalHits();
1640 m_probph2kpi[i] = dedxTrk->probPH();
1641 m_normph2kpi[i] = dedxTrk->normPH();
1642 }
1643
1644 //emc
1645 if ( (*itTrk)->isEmcShowerValid() )
1646 {
1647 RecEmcShower* emcTrk = (*itTrk)->emcShower();
1648 m_bjemc2kpi[i]=1;
1649 m_emc2kpi[i] = emcTrk->energy();
1650 m_evp2kpi[i] = emcTrk->energy()/ptrk;
1651 m_timecg2kpi[i] = emcTrk->time();
1652 }
1653
1654
1655 //muc
1656 if ( (*itTrk)->isMucTrackValid() )
1657 {
1658 RecMucTrack* mucTrk = (*itTrk)->mucTrack();
1659 double dpthp = mucTrk->depth()/25.0;//why?
1660 m_bjmuc2kpi[i]=1;
1661 m_depthmuc2kpi[i]=mucTrk->depth();
1662 m_layermuc2kpi[i] = mucTrk->numLayers();
1663 }
1664
1665 m_cosmdc2kpi[i] = cos(mdcTrk->theta());
1666 m_phimdc2kpi[i] = mdcTrk->phi();
1667
1668 m_pidnum2kpi[i]=(*itTrk)->partId();
1669
1670 if(m_skim4pi)
1671 {
1672 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
1673 {
1674 if(i<4) (*itTrk)->setPartId(3);
1675 }
1676
1677if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==1)
1678 {
1679 if(i<4) (*itTrk)->setPartId(4);
1680 }
1681
1682
1683if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==2)
1684 {
1685 if(i==0||i==1) (*itTrk)->setPartId(3);
1686 if(i==2||i==3) (*itTrk)->setPartId(5);
1687 }
1688
1689 }
1690 //pid
1691 ParticleID *pid = ParticleID::instance();
1692 pid->init();
1693 pid->setMethod(pid->methodProbability());
1694 pid->setChiMinCut(4);
1695 pid->setRecTrack(*itTrk);
1696 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
1697 pid->identify(pid->onlyPion() | pid->onlyKaon() | pid->onlyProton());
1698 pid->identify(pid->onlyElectron() | pid->onlyMuon());
1699 pid->calculate();
1700 if(!(pid->IsPidInfoValid())) continue;
1701 m_costpid2kpi[i] = cos(mdcTrk->theta());
1702
1703
1704 m_probe2kpi[i] = pid->probElectron();
1705 m_probmu2kpi[i] = pid->probMuon();
1706 m_probpi2kpi[i] = pid->probPion();
1707 m_probk2kpi[i] = pid->probKaon();
1708 m_probpr2kpi[i] = pid->probProton();
1709
1710 m_chipidxpid2kpi[i] = pid->chiDedx(2);
1711 m_chipitof1pid2kpi[i] = pid->chiTof1(2);
1712 m_chipitof2pid2kpi[i] = pid->chiTof2(2);
1713 m_chipitofpid2kpi[i]=pid->chiTof(2);
1714 m_chipitofepid2kpi[i]=pid->chiTofE(2);
1715 m_chipitofqpid2kpi[i]=pid->chiTofQ(2);
1716 m_probpidxpid2kpi[i]=pid->probDedx(2);
1717 m_probpitofpid2kpi[i]=pid->probTof(2);
1718
1719 m_chikdxpid2kpi[i] = pid->chiDedx(3);
1720 m_chiktof1pid2kpi[i] = pid->chiTof1(3);
1721 m_chiktof2pid2kpi[i] = pid->chiTof2(3);
1722 m_chiktofpid2kpi[i]=pid->chiTof(3);
1723 m_chiktofepid2kpi[i]=pid->chiTofE(3);
1724 m_chiktofqpid2kpi[i]=pid->chiTofQ(3);
1725 m_probkdxpid2kpi[i]=pid->probDedx(3);
1726 m_probktofpid2kpi[i]=pid->probTof(3);
1727
1728 m_chiprdxpid2kpi[i] = pid->chiDedx(4);
1729 m_chiprtof1pid2kpi[i] = pid->chiTof1(4);
1730 m_chiprtof2pid2kpi[i] = pid->chiTof2(4);
1731 m_chiprtofpid2kpi[i]=pid->chiTof(4);
1732 m_chiprtofepid2kpi[i]=pid->chiTofE(4);
1733 m_chiprtofqpid2kpi[i]=pid->chiTofQ(4);
1734 m_probprdxpid2kpi[i]=pid->probDedx(4);
1735 m_probprtofpid2kpi[i]=pid->probTof(4);
1736
1737
1738 if ((*itTrk)->isTofTrackValid())
1739 {
1740 SmartRefVector<RecTofTrack> dstTofTrackCol = (*itTrk)->tofTrack();
1741 SmartRefVector<RecTofTrack>::iterator iter_tof = dstTofTrackCol.begin();
1742 for( ; iter_tof!=dstTofTrackCol.end(); iter_tof++ )
1743 {
1744
1745 m_trk_pidtype=0;
1746 m_trk_type=0;
1747 if((((m_mchic2kpi>3.39&&m_mchic2kpi<3.44)||(m_mchic2kpi>3.5&&m_mchic2kpi<3.57))&&m_chis2kpi<20&&m_mpprecall<3.06&&m_energy[0]<0.3&&m_energy[0]>0.1&&m_cla2kpi==0&&m_probpi2kpi[i]>m_probk2kpi[i]&&m_probpi2kpi[i]>m_probpr2kpi[i]))
1748 m_trk_pidtype=1;
1749
1750 if((((m_mchic2kpi>3.39&&m_mchic2kpi<3.44)||(m_mchic2kpi>3.5&&m_mchic2kpi<3.57))&&m_chis2kpi<20&&m_mpprecall<3.06&&m_energy[0]<0.3&&m_energy[0]>0.1&&m_cla2kpi==1&&m_probk2kpi[i]>m_probpi2kpi[i]&&m_probk2kpi[i]>m_probpr2kpi[i]))
1751 m_trk_pidtype=2;
1752
1753 if((((m_mchic2kpi>3.39&&m_mchic2kpi<3.44)||(m_mchic2kpi>3.5&&m_mchic2kpi<3.57))&&m_chis2kpi<20&&m_mpprecall<3.06&&m_energy[0]<0.3&&m_energy[0]>0.1&&m_cla2kpi==2&&m_probpr2kpi[i]>m_probpi2kpi[i]&&m_probpr2kpi[i]>m_probk2kpi[i]&&(i==2||i==3)))
1754 m_trk_pidtype=3;
1755
1756 if((((m_mchic2kpi>3.39&&m_mchic2kpi<3.44)||(m_mchic2kpi>3.5&&m_mchic2kpi<3.57))&&m_chis2kpi<20&&m_mpprecall<3.06&&m_energy[0]<0.3&&m_energy[0]>0.1&&m_cla2kpi==0))
1757 { RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);
1758 m_trk_type=1;
1759 }
1760 if((((m_mchic2kpi>3.39&&m_mchic2kpi<3.44)||(m_mchic2kpi>3.5&&m_mchic2kpi<3.57))&&m_chis2kpi<20&&m_mpprecall<3.06&&m_energy[0]<0.3&&m_energy[0]>0.1&&m_cla2kpi==1))
1761 {RecMdcKalTrack::setPidType (RecMdcKalTrack::kaon);
1762 m_trk_type=2;
1763 }
1764 if((((m_mchic2kpi>3.39&&m_mchic2kpi<3.44)||(m_mchic2kpi>3.5&&m_mchic2kpi<3.57))&&m_chis2kpi<20&&m_mpprecall<3.06&&m_energy[0]<0.3&&m_energy[0]>0.1&&m_cla2kpi==2&&(i==2||i==3)))
1765 {
1766 RecMdcKalTrack::setPidType (RecMdcKalTrack::proton);
1767 m_trk_type=3;
1768
1769 }
1770 m_trk_trackid = (*iter_tof)->trackID();
1771 m_trk_tofid = (*iter_tof)->tofID();
1772 TofHitStatus* hitStatus = new TofHitStatus;
1773 hitStatus->setStatus( (*iter_tof)->status() );
1774
1775 m_trk_raw = hitStatus->is_raw();
1776 m_trk_readout = hitStatus->is_readout();
1777 m_trk_counter = hitStatus->is_counter();
1778 m_trk_cluster = hitStatus->is_cluster();
1779 m_trk_barrel = hitStatus->is_barrel();
1780 m_trk_east = hitStatus->is_east();
1781 m_trk_layer = hitStatus->layer();
1782 delete hitStatus;
1783 m_trk_path = (*iter_tof)->path();
1784 m_trk_zrhit = (*iter_tof)->zrhit();
1785 m_trk_ph = (*iter_tof)->ph();
1786 m_trk_tof = (*iter_tof)->tof();
1787 m_trk_beta = (*iter_tof)->beta();
1788 m_trk_texpe = (*iter_tof)->texpElectron();
1789 m_trk_texpmu = (*iter_tof)->texpMuon();
1790 m_trk_texppi = (*iter_tof)->texpPion();
1791 m_trk_texpk = (*iter_tof)->texpKaon();
1792 m_trk_texpp = (*iter_tof)->texpProton();
1793
1794 m_trk_offe = (*iter_tof)->toffsetElectron();
1795 m_trk_offmu = (*iter_tof)->toffsetMuon();
1796 m_trk_offpi = (*iter_tof)->toffsetPion();
1797 m_trk_offk = (*iter_tof)->toffsetKaon();
1798 m_trk_offp = (*iter_tof)->toffsetProton();
1799 m_trk_quality = (*iter_tof)->quality();
1800 m_trk_ppmdc=mdcTrk->p();
1801 m_trk_ppkal=mdcKalTrk->p();
1802 m_trk_cosmdc=cos(mdcTrk->theta());
1803 m_trk_phimdc=mdcTrk->phi();
1804 m_trk_coskal=cos(mdcKalTrk->theta());
1805 m_trk_phikal=mdcKalTrk->phi();
1806
1807 if(i==0||i==2)
1808 m_trk_charge=1;
1809 if(i==2||i==3)
1810 m_trk_charge=2;
1811 trk_tuple->write();
1812 }
1813 }
1814 }
1815}
1816
1817 }
1818// Ncut[10]++;
1819
1820if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
1821{
1822 jGam2kpi.push_back(igam2kpi);
1823}
1824
1825 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
1826 if(i>=evtRecTrkCol->size()) break;
1827 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
1828 if(!(*itTrk)->isEmcShowerValid()) continue;
1829if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
1830 {
1831 if(i!=igam2kpi) jGam2kpi.push_back((*itTrk)->trackId());
1832 }
1833 else{
1834 jGam2kpi.push_back((*itTrk)->trackId());
1835 }
1836
1837 }
1838
1839int ngam2kpi=jGam2kpi.size();
1840
1841if(m_rootput)
1842{
1843for(int i = 0; i< ngam2kpi; i++) {
1844 if(i>=evtRecTrkCol->size()) break;
1845 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + jGam2kpi[i];
1846 if(!(*itTrk)->isEmcShowerValid()) continue;
1847 RecEmcShower *emcTrk = (*(evtRecTrkCol->begin()+jGam2kpi[i]))->emcShower();
1848 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
1849 double dthe = 200.;
1850 double dphi = 200.;
1851 double dang = 200.;
1852// double dthert = 200.;
1853// double dphirt = 200.;
1854// double dangrt = 200.;
1855 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
1856 if(j>=evtRecTrkCol->size()) break;
1857 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
1858 if(!(*jtTrk)->isExtTrackValid()) continue;
1859 RecExtTrack *extTrk = (*jtTrk)->extTrack();
1860 if(extTrk->emcVolumeNumber() == -1) continue;
1861 Hep3Vector extpos = extTrk->emcPosition();
1862 double angd = extpos.angle(emcpos);
1863 double thed = extpos.theta() - emcpos.theta();
1864 double phid = extpos.deltaPhi(emcpos);
1865 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
1866 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
1867// if(fabs(thed) < fabs(dthe)) dthe = thed;
1868// if(fabs(phid) < fabs(dphi)) dphi = phid;
1869// if(angd < dang) dang = angd;
1870 if(angd < dang) { dang = angd;
1871 dthe = thed;
1872 dphi = phid;
1873 }
1874
1875 }
1876 if(dang>=200) continue;
1877 double eraw = emcTrk->energy();
1878 dthe = dthe * 180 / (CLHEP::pi);
1879 dphi = dphi * 180 / (CLHEP::pi);
1880 dang = dang * 180 / (CLHEP::pi);
1881// dthert = dthert * 180 / (CLHEP::pi);
1882// dphirt = dphirt * 180 / (CLHEP::pi);
1883// dangrt = dangrt * 180 / (CLHEP::pi);
1884 m_numHits[i]= emcTrk->numHits();
1885 m_secondmoment[i] = emcTrk->secondMoment();
1886 m_latmoment[i] = emcTrk->latMoment();
1887 m_timegm[i] = emcTrk->time();
1888 m_cellId[i]=emcTrk->cellId();
1889 m_module[i]=emcTrk->module();
1890 m_a20Moment[i]=emcTrk->a20Moment();
1891 m_a42Moment[i]=emcTrk->a42Moment();
1892 m_getShowerId[i]=emcTrk->getShowerId();
1893 m_getClusterId[i]=emcTrk->getClusterId();
1894 m_getEAll[i]=emcTrk->getEAll();
1895 m_x[i]= emcTrk->x();
1896 m_y[i]= emcTrk->y();
1897 m_z[i]= emcTrk->z();
1898 m_cosemc[i] = cos(emcTrk->theta());
1899 m_phiemc[i] = emcTrk->phi();
1900 m_energy[i] = emcTrk->energy();
1901 m_eSeed[i] = emcTrk->eSeed();
1902 m_e3x3[i] = emcTrk->e3x3();
1903 m_e5x5[i] = emcTrk->e5x5();
1904 m_dang4c[i] = dang;
1905 m_dthe4c[i] = dthe;
1906 m_dphi4c[i] = dphi;
1907// m_dang4crt[i] = dangrt;
1908// m_dthe4crt[i] = dthert;
1909// m_dphi4crt[i] = dphirt;
1910
1911}
1912}
1913}
1914
1915}
1916
1917
1918
1919
1920
1921 //4pi
1922 if(m_skim4pi)
1923 {
1924
1925if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
1926 {m_subalg1->execute(); // save gam4pi data
1927 Ncut[6]++;
1928 }
1929 }
1930
1931 //4k
1932 if(m_skim4k)
1933 {
1934
1935if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==1&&dtpr2kpilst[0]<-0.4&&dtpr2kpilst[1]<-0.4&&dtpr2kpilst[2]<-0.4&&dtpr2kpilst[3]<-0.4)
1936 {m_subalg2->execute(); // save gam4k data
1937 Ncut[7]++;
1938 }
1939 }
1940
1941 if(m_skim2pi2pr)
1942 {
1943if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==2)
1944 {m_subalg3->execute(); // save gam 2(pi p) data
1945 // pi+ pi- with low momentum and p pbar with high momentum.
1946 Ncut[8]++;
1947 }
1948 }
1949
1950
1951
1952//cout<<"chis4c2kpilst="<<chis4c2kpilst<<endl;
1953if(m_rootput)
1954{
1955
1956if((mppreclst<3.06&&chis4c2kpilst<40)||((meelst>3.06&&meelst<3.12)&&(fabs(mppreclst-3.097)<0.01)))
1957 { Ncut[9]++;
1958 m_tuple1->write();
1959
1960 }
1961}
1962
1963 return StatusCode::SUCCESS;
1964 }
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition EvtRecTrack.h:111
Vp4
std::vector< HepLorentzVector > Vp4
Definition Gam4pikp.cxx:53
xmass
const double xmass[5]
Definition Gam4pikp.cxx:50
mpro
const double mpro
Definition Gam4pikp.cxx:49
mk
const double mk
Definition Gam4pikp.cxx:48
Vdouble
std::vector< double > Vdouble
Definition Gam4pikp.cxx:54
mpi
const double mpi
Definition Gam4pikp.cxx:47
Vint
std::vector< int > Vint
Definition Gam4pikp.cxx:52
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition KarLud.h:35
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
DstEmcShower::cellId
int cellId() const
Definition DstEmcShower.h:32
DstEmcShower::latMoment
double latMoment() const
Definition DstEmcShower.h:52
DstEmcShower::a42Moment
double a42Moment() const
Definition DstEmcShower.h:54
DstEmcShower::eSeed
double eSeed() const
Definition DstEmcShower.h:47
DstEmcShower::theta
double theta() const
Definition DstEmcShower.h:38
DstEmcShower::module
int module() const
Definition DstEmcShower.h:33
DstEmcShower::e3x3
double e3x3() const
Definition DstEmcShower.h:48
DstEmcShower::phi
double phi() const
Definition DstEmcShower.h:39
DstEmcShower::secondMoment
double secondMoment() const
Definition DstEmcShower.h:51
DstEmcShower::x
double x() const
Definition DstEmcShower.h:35
DstEmcShower::e5x5
double e5x5() const
Definition DstEmcShower.h:49
DstEmcShower::time
double time() const
Definition DstEmcShower.h:50
DstEmcShower::z
double z() const
Definition DstEmcShower.h:37
DstEmcShower::numHits
int numHits() const
Definition DstEmcShower.h:30
DstEmcShower::a20Moment
double a20Moment() const
Definition DstEmcShower.h:53
DstEmcShower::energy
double energy() const
Definition DstEmcShower.h:45
DstEmcShower::y
double y() const
Definition DstEmcShower.h:36
DstExtTrack::emcPosition
const Hep3Vector emcPosition() const
Definition DstExtTrack.h:126
DstExtTrack::emcVolumeNumber
const int emcVolumeNumber() const
Definition DstExtTrack.h:132
DstMdcDedx::probPH
double probPH() const
Definition DstMdcDedx.h:66
DstMdcDedx::chiE
double chiE() const
Definition DstMdcDedx.h:59
DstMdcDedx::numTotalHits
int numTotalHits() const
Definition DstMdcDedx.h:65
DstMdcDedx::numGoodHits
int numGoodHits() const
Definition DstMdcDedx.h:64
DstMdcDedx::normPH
double normPH() const
Definition DstMdcDedx.h:67
DstMdcDedx::chiPi
double chiPi() const
Definition DstMdcDedx.h:61
DstMdcDedx::chiK
double chiK() const
Definition DstMdcDedx.h:62
DstMdcDedx::chiMu
double chiMu() const
Definition DstMdcDedx.h:60
DstMdcDedx::chiP
double chiP() const
Definition DstMdcDedx.h:63
DstMdcKalTrack::theta
const double theta() const
Definition DstMdcKalTrack.h:70
DstMdcKalTrack::phi
const double phi() const
Definition DstMdcKalTrack.h:71
DstMdcKalTrack::px
const double px() const
Definition DstMdcKalTrack.h:67
DstMdcKalTrack::setPidType
static void setPidType(PidType pidType)
Definition DstMdcKalTrack.h:59
DstMdcKalTrack::pz
const double pz() const
Definition DstMdcKalTrack.h:69
DstMdcKalTrack::py
const double py() const
Definition DstMdcKalTrack.h:68
DstMdcKalTrack::p
const double p() const
Definition DstMdcKalTrack.h:82
DstMdcKalTrack::charge
const int charge() const
Definition DstMdcKalTrack.h:65
DstMdcKalTrack::pxy
const double pxy() const
Definition DstMdcKalTrack.h:66
DstMdcTrack::theta
const double theta() const
Definition DstMdcTrack.h:59
DstMdcTrack::r
const double r() const
Definition DstMdcTrack.h:64
DstMdcTrack::py
const double py() const
Definition DstMdcTrack.h:56
DstMdcTrack::err
const HepSymMatrix err() const
Definition DstMdcTrack.cxx:155
DstMdcTrack::charge
const int charge() const
Definition DstMdcTrack.h:53
DstMdcTrack::px
const double px() const
Definition DstMdcTrack.h:55
DstMdcTrack::phi
const double phi() const
Definition DstMdcTrack.h:60
DstMdcTrack::pz
const double pz() const
Definition DstMdcTrack.h:57
DstMdcTrack::pxy
const double pxy() const
Definition DstMdcTrack.h:54
DstMdcTrack::helix
const HepVector helix() const
......
Definition DstMdcTrack.cxx:151
DstMdcTrack::z
const double z() const
Definition DstMdcTrack.h:63
DstMdcTrack::p
const double p() const
Definition DstMdcTrack.h:58
DstMdcTrack::y
const double y() const
Definition DstMdcTrack.h:62
DstMdcTrack::x
const double x() const
Definition DstMdcTrack.h:61
DstMucTrack::depth
double depth() const
Definition DstMucTrack.h:45
DstMucTrack::numLayers
int numLayers() const
Definition DstMucTrack.h:42
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
KalmanKinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition KalmanKinematicFit.cxx:353
KalmanKinematicFit::pfit
HepLorentzVector pfit(int n) const
Definition KalmanKinematicFit.h:154
KalmanKinematicFit::instance
static KalmanKinematicFit * instance()
Definition KalmanKinematicFit.cxx:20
ParticleIDBase::useTof2
int useTof2() const
Definition ParticleIDBase.h:67
ParticleIDBase::onlyProton
int onlyProton() const
Definition ParticleIDBase.h:84
ParticleIDBase::methodProbability
int methodProbability() const
Definition ParticleIDBase.h:88
ParticleIDBase::useDedx
int useDedx() const
Definition ParticleIDBase.h:64
ParticleIDBase::onlyMuon
int onlyMuon() const
Definition ParticleIDBase.h:81
ParticleIDBase::onlyKaon
int onlyKaon() const
Definition ParticleIDBase.h:83
ParticleIDBase::onlyElectron
int onlyElectron() const
Definition ParticleIDBase.h:80
ParticleIDBase::onlyPion
int onlyPion() const
Definition ParticleIDBase.h:82
ParticleIDBase::useTof1
int useTof1() const
Definition ParticleIDBase.h:66
ParticleIDBase::setChiMinCut
void setChiMinCut(const double chi=4)
Definition ParticleIDBase.h:39
ParticleIDBase::setRecTrack
void setRecTrack(EvtRecTrack *trk)
Definition ParticleIDBase.h:37
ParticleID::probKaon
double probKaon() const
Definition ParticleID.h:124
ParticleID::setMethod
void setMethod(const int method)
Definition ParticleID.h:94
ParticleID::probTof
double probTof(int n) const
ParticleID::chiTof2
double chiTof2(int n) const
ParticleID::chiTofE
double chiTofE(int n) const
ParticleID::identify
void identify(const int pidcase)
Definition ParticleID.h:103
ParticleID::probMuon
double probMuon() const
Definition ParticleID.h:122
ParticleID::probElectron
double probElectron() const
Definition ParticleID.h:121
ParticleID::usePidSys
void usePidSys(const int pidsys)
Definition ParticleID.h:97
ParticleID::chiTofQ
double chiTofQ(int n) const
ParticleID::instance
static ParticleID * instance()
Definition ParticleID.cxx:22
ParticleID::IsPidInfoValid
bool IsPidInfoValid() const
ParticleID::probDedx
double probDedx(int n) const
ParticleID::probPion
double probPion() const
Definition ParticleID.h:123
ParticleID::chiTof1
double chiTof1(int n) const
ParticleID::calculate
void calculate()
Definition ParticleID.cxx:97
ParticleID::init
void init()
Definition ParticleID.cxx:27
ParticleID::probProton
double probProton() const
Definition ParticleID.h:125
ParticleID::chiTof
double chiTof(int n) const
ParticleID::chiDedx
double chiDedx(int n) const
RecEmcShower::getClusterId
RecEmcID getClusterId() const
Definition RecEmcShower.h:59
RecEmcShower::getShowerId
RecEmcID getShowerId() const
Definition RecEmcShower.h:55
RecEmcShower::getEAll
RecEmcEnergy getEAll() const
Definition RecEmcShower.h:92
RecMdcKalTrack::getZHelix
const HepVector & getZHelix() const
Definition RecMdcKalTrack.h:66
RecMdcKalTrack::getZError
const HepSymMatrix & getZError() const
Definition RecMdcKalTrack.h:67
TofHitStatus::is_barrel
bool is_barrel() const
Definition TofHitStatus.h:26
TofHitStatus::layer
unsigned int layer() const
Definition TofHitStatus.h:28
TofHitStatus::is_cluster
bool is_cluster() const
Definition TofHitStatus.h:25
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition TofHitStatus.cxx:134
TofHitStatus::is_counter
bool is_counter() const
Definition TofHitStatus.h:24
TofHitStatus::is_east
bool is_east() const
Definition TofHitStatus.h:27
TofHitStatus::is_readout
bool is_readout() const
Definition TofHitStatus.h:23
TofHitStatus::is_raw
bool is_raw() const
Definition TofHitStatus.h:22
TrackPool::setVBeamPosition
void setVBeamPosition(const HepSymMatrix VBeamPosition)
Definition TrackPool.h:144
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:22
TrackPool::setBeamPosition
void setBeamPosition(const HepPoint3D BeamPosition)
Definition TrackPool.h:143
VertexFit::wtrk
WTrackParameter wtrk(int n) const
Definition VertexFit.h:79
VertexFit::init
void init()
Definition VertexFit.cxx:29
VertexFit::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFit.cxx:89
VertexFit::Evx
HepSymMatrix Evx(int n) const
Definition VertexFit.h:87
VertexFit::vx
HepPoint3D vx(int n) const
Definition VertexFit.h:85
VertexFit::instance
static VertexFit * instance()
Definition VertexFit.cxx:15
VertexFit::Swim
void Swim(int n)
Definition VertexFit.h:59
VertexFit::Fit
bool Fit()
Definition VertexFit.cxx:301
VertexParameter::setEvx
void setEvx(const HepSymMatrix &eVx)
Definition VertexParameter.h:29
VertexParameter::setVx
void setVx(const HepPoint3D &vx)
Definition VertexParameter.h:54
ecms
const double ecms
Definition inclkstar.cxx:42
EventModel::EvtRec::EvtRecEvent
_EXTERN_ std::string EvtRecEvent
Definition EventModel.h:116
EventModel::EvtRec::EvtRecTrackCol
_EXTERN_ std::string EvtRecTrackCol
Definition EventModel.h:117
ptrk
float ptrk
Definition prepare_data.cxx:28
p3
double double * p3
Definition qcdloop1.h:76
pi
const float pi
Definition vector3.h:133

◆ finalize()

StatusCode Gam4pikp::finalize ( )

Definition at line 1968 of file Gam4pikp.cxx.

1968 {
1969// cout<<"total number: "<<Ncut[0]<<endl;
1970// cout<<"nGood>=4, nCharge==0: "<<Ncut[1]<<endl;
1971// cout<<"nGam>=1: "<<Ncut[2]<<endl;
1972// cout<<"cgp>2 cgm>2: "<<Ncut[3]<<endl;
1973// cout<<"2+ 2-: "<<Ncut[4]<<endl;
1974// cout<<"all cg cycle, chisq<200: "<<Ncut[5]<<endl;
1975// cout<<"4pi ok: "<<Ncut[6]<<endl;
1976// cout<<"4k ok: "<<Ncut[7]<<endl;
1977// cout<<"2pi 2p ok: "<<Ncut[8]<<endl;
1978// cout<<"ntuple write: "<<Ncut[9]<<endl;
1979 MsgStream log(msgSvc(), name());
1980 log << MSG::INFO << "in finalize()" << endmsg;
1981 return StatusCode::SUCCESS;
1982 }

◆ initialize()

StatusCode Gam4pikp::initialize ( )

Definition at line 84 of file Gam4pikp.cxx.

84 {
85 MsgStream log(msgSvc(), name());
86
87 log << MSG::INFO << "in initialize()" << endmsg;
88 // setFilterPassed(false);
89 StatusCode status;
90
91 if(m_rootput)
92 {
93
94
95 NTuplePtr trk0(ntupleSvc(), "FILE1/toftrk");
96 if ( trk0 ) trk_tuple = trk0;
97 else { trk_tuple = ntupleSvc()->book ("FILE1/toftrk", CLID_ColumnWiseTuple, "ks N-Tuple example");
98 if ( trk_tuple ) {
99 status = trk_tuple->addItem ("trk_trackid", m_trk_trackid );
100 status = trk_tuple->addItem ("trk_tofid", m_trk_tofid );
101 status = trk_tuple->addItem ("trk_raw", m_trk_raw);
102 status = trk_tuple->addItem ("trk_readout", m_trk_readout);
103 status = trk_tuple->addItem ("trk_counter", m_trk_counter);
104 status = trk_tuple->addItem ("trk_cluster", m_trk_cluster );
105 status = trk_tuple->addItem ("trk_barrel", m_trk_barrel );
106 status = trk_tuple->addItem ("trk_east", m_trk_east );
107 status = trk_tuple->addItem ("trk_layer", m_trk_layer );
108 status = trk_tuple->addItem ("trk_path", m_trk_path );
109 status = trk_tuple->addItem ("trk_zrhit", m_trk_zrhit );
110 status = trk_tuple->addItem ("trk_ph", m_trk_ph );
111 status = trk_tuple->addItem ("trk_tof", m_trk_tof );
112 status = trk_tuple->addItem ("trk_beta", m_trk_beta );
113
114 status = trk_tuple->addItem ("trk_texpe", m_trk_texpe );
115 status = trk_tuple->addItem ("trk_texpmu", m_trk_texpmu );
116 status = trk_tuple->addItem ("trk_texppi", m_trk_texppi );
117 status = trk_tuple->addItem ("trk_texpk", m_trk_texpk );
118 status = trk_tuple->addItem ("trk_texpp", m_trk_texpp );
119 status = trk_tuple->addItem ("trk_offe", m_trk_offe );
120 status = trk_tuple->addItem ("trk_offmu", m_trk_offmu );
121 status = trk_tuple->addItem ("trk_offpi", m_trk_offpi );
122 status = trk_tuple->addItem ("trk_offk", m_trk_offk );
123 status = trk_tuple->addItem ("trk_offp", m_trk_offp );
124 status = trk_tuple->addItem ("trk_quality", m_trk_quality );
125 status = trk_tuple->addItem ("trk_type", m_trk_type );
126 status = trk_tuple->addItem ("trk_pidtype", m_trk_pidtype );
127 status = trk_tuple->addItem ("trk_ppmdc", m_trk_ppmdc );
128 status = trk_tuple->addItem ("trk_ppmdc", m_trk_ppmdc );
129 status = trk_tuple->addItem ("trk_ptmdc", m_trk_ptmdc );
130 status = trk_tuple->addItem ("trk_ppkal", m_trk_ppkal );
131 status = trk_tuple->addItem ("trk_ptkal", m_trk_ptkal );
132 status = trk_tuple->addItem ("trk_cosmdc", m_trk_cosmdc );
133 status = trk_tuple->addItem ("trk_phimdc", m_trk_phimdc );
134 status = trk_tuple->addItem ("trk_coskal", m_trk_coskal );
135 status = trk_tuple->addItem ("trk_phikal", m_trk_phikal );
136
137 status = trk_tuple->addItem ("trk_charge", m_trk_charge );
138
139 }
140 else {
141 log << MSG::ERROR << " Cannot book N-tuple:" << long(trk_tuple) << endmsg;
142 return StatusCode::FAILURE;
143 }
144 }
145
146
147
148
149 NTuplePtr nt1(ntupleSvc(), "FILE1/total4c");
150 if ( nt1 ) m_tuple1 = nt1;
151 else {
152 m_tuple1 = ntupleSvc()->book ("FILE1/total4c", CLID_ColumnWiseTuple, "ks N-Tuple example");
153 if ( m_tuple1 ) {
154
155 status = m_tuple1->addItem ("run", m_run );
156 status = m_tuple1->addItem ("rec", m_rec );
157 status = m_tuple1->addItem ("mpprecall", m_mpprecall );
158 status = m_tuple1->addItem ("meeall", m_meeall );
159 status = m_tuple1->addItem ("ncgjs", m_ncgjs );
160 status = m_tuple1->addItem ("cla2kpi", m_cla2kpi );
161 status = m_tuple1->addItem("indexmc", m_idxmc, 0, 100);
162 status = m_tuple1->addIndexedItem("pdgid", m_idxmc, m_pdgid);
163
164 status = m_tuple1->addIndexedItem("motheridx", m_idxmc, m_motheridx);
165 status = m_tuple1->addItem("indexmdc", m_idxmdc, 0, 5000);
166 status = m_tuple1->addIndexedItem ("x0js", m_idxmdc, m_x0js);
167 status = m_tuple1->addIndexedItem ("y0js", m_idxmdc, m_y0js);
168 status = m_tuple1->addIndexedItem ("z0js",m_idxmdc, m_z0js);
169 status = m_tuple1->addIndexedItem ("r0js",m_idxmdc, m_r0js);
170 status = m_tuple1->addIndexedItem ("Rxyjs",m_idxmdc, m_Rxyjs);
171 status = m_tuple1->addIndexedItem ("Rzjs",m_idxmdc, m_Rzjs);
172 status = m_tuple1->addIndexedItem ("Rnxyjs",m_idxmdc, m_Rnxyjs);
173 status = m_tuple1->addIndexedItem ("phinjs",m_idxmdc, m_phinjs);
174 status = m_tuple1->addIndexedItem ("Rnzjs",m_idxmdc, m_Rnzjs);
175 status = m_tuple1->addItem ("ncy20", m_ncy20);
176 status = m_tuple1->addItem ("ncy30", m_ncy30);
177 status = m_tuple1->addIndexedItem("angjs5", m_idxmdc, m_angjs5);
178 status = m_tuple1->addIndexedItem("nearjs5", m_idxmdc, m_nearjs5);
179 status = m_tuple1->addIndexedItem("angjs6", m_idxmdc, m_angjs6);
180 status = m_tuple1->addIndexedItem("nearjs6", m_idxmdc, m_nearjs6);
181 status = m_tuple1->addIndexedItem("ang4pi5", m_idxmdc, m_ang4pi5);
182 status = m_tuple1->addIndexedItem("near4pi5", m_idxmdc, m_near4pi5);
183 status = m_tuple1->addIndexedItem("ang4pi6", m_idxmdc, m_ang4pi6);
184 status = m_tuple1->addIndexedItem("near4pi6", m_idxmdc, m_near4pi6);
185 status = m_tuple1->addIndexedItem("ppmdcjs", m_idxmdc, m_ppmdcjs);
186 status = m_tuple1->addIndexedItem("pxmdcjs", m_idxmdc, m_pxmdcjs);
187 status = m_tuple1->addIndexedItem("pymdcjs", m_idxmdc, m_pymdcjs);
188 status = m_tuple1->addIndexedItem("pzmdcjs", m_idxmdc, m_pzmdcjs);
189 status = m_tuple1->addIndexedItem("ppkaljs", m_idxmdc, m_ppkaljs);
190 status = m_tuple1->addIndexedItem("ptmdcjs", m_idxmdc, m_ptmdcjs);
191 status = m_tuple1->addIndexedItem("ptkaljs", m_idxmdc, m_ptkaljs);
192 status = m_tuple1->addIndexedItem("ppmdc2kpi", m_idxmdc, m_ppmdc2kpi);
193 status = m_tuple1->addIndexedItem("pxmdc2kpi", m_idxmdc, m_pxmdc2kpi);
194 status = m_tuple1->addIndexedItem("pymdc2kpi", m_idxmdc, m_pymdc2kpi);
195 status = m_tuple1->addIndexedItem("pzmdc2kpi", m_idxmdc, m_pzmdc2kpi);
196 status = m_tuple1->addIndexedItem("ppkal2kpi", m_idxmdc, m_ppkal2kpi);
197 status = m_tuple1->addIndexedItem("ptmdc2kpi", m_idxmdc, m_ptmdc2kpi);
198 status = m_tuple1->addIndexedItem("charge2kpi", m_idxmdc, m_charge2kpi);
199 status = m_tuple1->addIndexedItem("ptkal2kpi", m_idxmdc, m_ptkal2kpi);
200 status = m_tuple1->addItem ("cy2pi", m_cy2kpi, 0, 100 );
201 status = m_tuple1->addIndexedItem("comcs2kpi", m_cy2kpi, m_comcs2kpi);
202 status = m_tuple1->addItem ("chiejs", m_idxmdc, m_chiejs);
203 status = m_tuple1->addItem ("chimujs", m_idxmdc, m_chimujs);
204 status = m_tuple1->addItem ("chipijs", m_idxmdc, m_chipijs);
205 status = m_tuple1->addItem ("chikjs", m_idxmdc, m_chikjs);
206 status = m_tuple1->addItem ("chipjs", m_idxmdc, m_chipjs);
207 status = m_tuple1->addItem ("ghitjs", m_idxmdc, m_ghitjs);
208 status = m_tuple1->addItem ("thitjs", m_idxmdc, m_thitjs);
209 status = m_tuple1->addIndexedItem("probphjs", m_idxmdc, m_probphjs);
210 status = m_tuple1->addIndexedItem("normphjs", m_idxmdc, m_normphjs);
211 status = m_tuple1->addItem ("pdg", m_idxmdc, m_pdg);
212 status = m_tuple1->addItem ("cbmc", m_idxmdc, m_cbmc);
213 status = m_tuple1->addIndexedItem("sigmaetof2kpi", m_idxmdc, m_sigmaetof2kpi);
214 status = m_tuple1->addIndexedItem("sigmamutof2kpi", m_idxmdc, m_sigmamutof2kpi);
215 status = m_tuple1->addIndexedItem("sigmapitof2kpi", m_idxmdc, m_sigmapitof2kpi);
216 status = m_tuple1->addIndexedItem("sigmaktof2kpi", m_idxmdc, m_sigmaktof2kpi);
217 status = m_tuple1->addIndexedItem("sigmaprtof2kpi", m_idxmdc, m_sigmaprtof2kpi);
218 status = m_tuple1->addIndexedItem("t0tof2kpi", m_idxmdc, m_t0tof2kpi);
219 status = m_tuple1->addIndexedItem("errt0tof2kpi", m_idxmdc, m_errt0tof2kpi);
220
221 status = m_tuple1->addItem ("chie2kpi", m_idxmdc, m_chie2kpi);
222 status = m_tuple1->addItem ("chimu2kpi", m_idxmdc, m_chimu2kpi);
223 status = m_tuple1->addItem ("chipi2kpi", m_idxmdc, m_chipi2kpi);
224 status = m_tuple1->addItem ("chik2kpi", m_idxmdc, m_chik2kpi);
225 status = m_tuple1->addItem ("chip2kpi", m_idxmdc, m_chip2kpi);
226 status = m_tuple1->addItem ("ghit2kpi", m_idxmdc, m_ghit2kpi);
227 status = m_tuple1->addItem ("thit2kpi", m_idxmdc, m_thit2kpi);
228 status = m_tuple1->addIndexedItem("probph2kpi", m_idxmdc, m_probph2kpi);
229 status = m_tuple1->addIndexedItem("normph2kpi", m_idxmdc, m_normph2kpi);
230 status = m_tuple1->addIndexedItem("pidnum2kpi", m_idxmdc, m_pidnum2kpi);
231 status = m_tuple1->addIndexedItem("bjmucjs", m_idxmdc, m_bjmucjs);
232 status = m_tuple1->addIndexedItem("bjmuc2kpi", m_idxmdc, m_bjmuc2kpi);
233 status = m_tuple1->addIndexedItem("bjemcjs", m_idxmdc, m_bjemcjs);
234 status = m_tuple1->addIndexedItem("bjemc2kpi", m_idxmdc, m_bjemc2kpi);
235 status = m_tuple1->addIndexedItem("bjtofjs", m_idxmdc, m_bjtofjs);
236 status = m_tuple1->addIndexedItem("bjtof2kpi", m_idxmdc, m_bjtof2kpi);
237 status = m_tuple1->addIndexedItem("bjtofvaljs", m_idxmdc, m_bjtofvaljs);
238 status = m_tuple1->addIndexedItem("bjtofval2kpi", m_idxmdc, m_bjtofval2kpi);
239
240 status = m_tuple1->addIndexedItem("emcjs", m_idxmdc, m_emcjs);
241 status = m_tuple1->addIndexedItem("evpjs", m_idxmdc, m_evpjs);
242 status = m_tuple1->addIndexedItem("timecgjs", m_idxmdc, m_timecgjs);
243 status = m_tuple1->addIndexedItem("depthjs", m_idxmdc, m_depthmucjs);
244 status = m_tuple1->addIndexedItem("layermucjs", m_idxmdc, m_layermucjs);
245
246 status = m_tuple1->addIndexedItem("emc2kpi", m_idxmdc, m_emc2kpi);
247 status = m_tuple1->addIndexedItem("evp2kpi", m_idxmdc, m_evp2kpi);
248 status = m_tuple1->addIndexedItem("timecg2kpi", m_idxmdc, m_timecg2kpi);
249 status = m_tuple1->addIndexedItem("depth2kpi", m_idxmdc, m_depthmuc2kpi);
250 status = m_tuple1->addIndexedItem("layermuc2kpi", m_idxmdc, m_layermuc2kpi);
251
252 status = m_tuple1->addIndexedItem("cotof1js", m_idxmdc, m_cotof1js);
253 status = m_tuple1->addIndexedItem("cotof2js", m_idxmdc, m_cotof2js);
254 status = m_tuple1->addIndexedItem("counterjs", m_idxmdc, m_counterjs);
255 status = m_tuple1->addIndexedItem("barreljs", m_idxmdc, m_barreljs);
256 status = m_tuple1->addIndexedItem("layertofjs", m_idxmdc, m_layertofjs);
257 status = m_tuple1->addIndexedItem("readoutjs", m_idxmdc, m_readoutjs);
258 status = m_tuple1->addIndexedItem("clusterjs", m_idxmdc, m_clusterjs);
259 status = m_tuple1->addIndexedItem("betajs", m_idxmdc, m_betajs);
260 status = m_tuple1->addIndexedItem("tofjs", m_idxmdc, m_tofjs);
261 status = m_tuple1->addIndexedItem("tofpathjs", m_idxmdc, m_tofpathjs);
262 status = m_tuple1->addIndexedItem("zhitjs", m_idxmdc, m_zhitjs);
263 status = m_tuple1->addIndexedItem("tofIDjs", m_idxmdc, m_tofIDjs);
264 status = m_tuple1->addIndexedItem("clusterIDjs", m_idxmdc, m_clusterIDjs);
265 status = m_tuple1->addIndexedItem("texejs", m_idxmdc, m_texejs);
266 status = m_tuple1->addIndexedItem("texmujs", m_idxmdc, m_texmujs);
267 status = m_tuple1->addIndexedItem("texpijs", m_idxmdc, m_texpijs);
268 status = m_tuple1->addIndexedItem("texkjs", m_idxmdc, m_texkjs);
269 status = m_tuple1->addIndexedItem("texprjs", m_idxmdc, m_texprjs);
270 status = m_tuple1->addIndexedItem("dtejs", m_idxmdc, m_dtejs);
271 status = m_tuple1->addIndexedItem("dtmujs", m_idxmdc, m_dtmujs);
272 status = m_tuple1->addIndexedItem("dtpijs", m_idxmdc, m_dtpijs);
273 status = m_tuple1->addIndexedItem("dtkjs", m_idxmdc, m_dtkjs);
274 status = m_tuple1->addIndexedItem("dtprjs", m_idxmdc, m_dtprjs);
275 status = m_tuple1->addIndexedItem("sigmaetofjs", m_idxmdc, m_sigmaetofjs);
276 status = m_tuple1->addIndexedItem("sigmamutofjs", m_idxmdc, m_sigmamutofjs);
277 status = m_tuple1->addIndexedItem("sigmapitofjs", m_idxmdc, m_sigmapitofjs);
278 status = m_tuple1->addIndexedItem("sigmaktofjs", m_idxmdc, m_sigmaktofjs);
279 status = m_tuple1->addIndexedItem("sigmaprtofjs", m_idxmdc, m_sigmaprtofjs);
280 status = m_tuple1->addIndexedItem("t0tofjs", m_idxmdc,m_t0tofjs);
281 status = m_tuple1->addIndexedItem("errt0tofjs", m_idxmdc,m_errt0tofjs);
282 status = m_tuple1->addIndexedItem("cotof12kpi", m_idxmdc, m_cotof12kpi);
283 status = m_tuple1->addIndexedItem("cotof22kpi", m_idxmdc, m_cotof22kpi);
284 status = m_tuple1->addIndexedItem("counter2kpi", m_idxmdc, m_counter2kpi);
285 status = m_tuple1->addIndexedItem("barrel2kpi", m_idxmdc, m_barrel2kpi);
286 status = m_tuple1->addIndexedItem("layertof2kpi", m_idxmdc, m_layertof2kpi);
287 status = m_tuple1->addIndexedItem("readout2kpi", m_idxmdc, m_readout2kpi);
288 status = m_tuple1->addIndexedItem("cluster2kpi", m_idxmdc, m_cluster2kpi);
289 status = m_tuple1->addIndexedItem("beta2kpi", m_idxmdc, m_beta2kpi);
290 status = m_tuple1->addIndexedItem("tof2kpi", m_idxmdc, m_tof2kpi);
291 status = m_tuple1->addIndexedItem("tofpath2kpi", m_idxmdc, m_tofpath2kpi);
292 status = m_tuple1->addIndexedItem("zhit2kpi", m_idxmdc, m_zhit2kpi);
293 status = m_tuple1->addIndexedItem("tofID2kpi", m_idxmdc, m_tofID2kpi);
294 status = m_tuple1->addIndexedItem("clusterID2kpi", m_idxmdc, m_clusterID2kpi);
295 status = m_tuple1->addIndexedItem("texe2kpi", m_idxmdc, m_texe2kpi);
296 status = m_tuple1->addIndexedItem("texmu2kpi", m_idxmdc, m_texmu2kpi);
297 status = m_tuple1->addIndexedItem("texpi2kpi", m_idxmdc, m_texpi2kpi);
298 status = m_tuple1->addIndexedItem("texk2kpi", m_idxmdc, m_texk2kpi);
299 status = m_tuple1->addIndexedItem("texpr2kpi", m_idxmdc, m_texpr2kpi);
300 status = m_tuple1->addIndexedItem("dte2kpi", m_idxmdc, m_dte2kpi);
301 status = m_tuple1->addIndexedItem("dtmu2kpi", m_idxmdc, m_dtmu2kpi);
302 status = m_tuple1->addIndexedItem("dtpi2kpi", m_idxmdc, m_dtpi2kpi);
303 status = m_tuple1->addIndexedItem("dtk2kpi", m_idxmdc, m_dtk2kpi);
304 status = m_tuple1->addIndexedItem("dtpr2kpi", m_idxmdc, m_dtpr2kpi);
305 status = m_tuple1->addIndexedItem("costpid2kpi", m_idxmdc, m_costpid2kpi);
306 status = m_tuple1->addIndexedItem("dedxpid2kpi", m_idxmdc, m_dedxpid2kpi);
307 status = m_tuple1->addIndexedItem("tof1pid2kpi", m_idxmdc, m_tof1pid2kpi);
308 status = m_tuple1->addIndexedItem("tof2pid2kpi", m_idxmdc, m_tof2pid2kpi);
309 status = m_tuple1->addIndexedItem("probe2kpi", m_idxmdc, m_probe2kpi);
310 status = m_tuple1->addIndexedItem("probmu2kpi", m_idxmdc, m_probmu2kpi);
311 status = m_tuple1->addIndexedItem("probpi2kpi", m_idxmdc, m_probpi2kpi);
312 status = m_tuple1->addIndexedItem("probk2kpi", m_idxmdc, m_probk2kpi);
313 status = m_tuple1->addIndexedItem("probpr2kpi", m_idxmdc, m_probpr2kpi);
314
315 status = m_tuple1->addIndexedItem("chipidxpid2kpi", m_idxmdc, m_chipidxpid2kpi);
316 status = m_tuple1->addIndexedItem("chipitof1pid2kpi", m_idxmdc, m_chipitof1pid2kpi);
317 status = m_tuple1->addIndexedItem("chipitof2pid2kpi", m_idxmdc, m_chipitof2pid2kpi);
318 status = m_tuple1->addIndexedItem("chipitofpid2kpi", m_idxmdc, m_chipitofpid2kpi);
319 status = m_tuple1->addIndexedItem("chipitofepid2kpi", m_idxmdc, m_chipitofepid2kpi);
320 status = m_tuple1->addIndexedItem("chipitofqpid2kpi", m_idxmdc, m_chipitofqpid2kpi);
321 status = m_tuple1->addIndexedItem("probpidxpid2kpi", m_idxmdc, m_probpidxpid2kpi);
322 status = m_tuple1->addIndexedItem("probpitofpid2kpi", m_idxmdc, m_probpitofpid2kpi);
323 status = m_tuple1->addIndexedItem("chikdxpid2kpi", m_idxmdc, m_chikdxpid2kpi);
324 status = m_tuple1->addIndexedItem("chiktof1pid2kpi", m_idxmdc, m_chiktof1pid2kpi);
325 status = m_tuple1->addIndexedItem("chiktof2pid2kpi", m_idxmdc, m_chiktof2pid2kpi);
326 status = m_tuple1->addIndexedItem("chiktofpid2kpi", m_idxmdc, m_chiktofpid2kpi);
327 status = m_tuple1->addIndexedItem("chiktofepid2kpi", m_idxmdc, m_chiktofepid2kpi);
328 status = m_tuple1->addIndexedItem("chiktofqpid2kpi", m_idxmdc, m_chiktofqpid2kpi);
329 status = m_tuple1->addIndexedItem("probkdxpid2kpi", m_idxmdc, m_probkdxpid2kpi);
330 status = m_tuple1->addIndexedItem("probktofpid2kpi", m_idxmdc, m_probktofpid2kpi);
331
332 status = m_tuple1->addIndexedItem("chiprdxpid2kpi", m_idxmdc, m_chiprdxpid2kpi);
333 status = m_tuple1->addIndexedItem("chiprtof1pid2kpi", m_idxmdc, m_chiprtof1pid2kpi);
334 status = m_tuple1->addIndexedItem("chiprtof2pid2kpi", m_idxmdc, m_chiprtof2pid2kpi);
335 status = m_tuple1->addIndexedItem("chiprtofpid2kpi", m_idxmdc, m_chiprtofpid2kpi);
336 status = m_tuple1->addIndexedItem("chiprtofepid2kpi", m_idxmdc, m_chiprtofepid2kpi);
337 status = m_tuple1->addIndexedItem("chiprtofqpid2kpi", m_idxmdc, m_chiprtofqpid2kpi);
338 status = m_tuple1->addIndexedItem("probprdxpid2kpi", m_idxmdc, m_probprdxpid2kpi);
339 status = m_tuple1->addIndexedItem("probprtofpid2kpi", m_idxmdc, m_probprtofpid2kpi);
340
341 status = m_tuple1->addIndexedItem("cosmdcjs", m_idxmdc, m_cosmdcjs);
342 status = m_tuple1->addIndexedItem("phimdcjs", m_idxmdc, m_phimdcjs);
343 status = m_tuple1->addIndexedItem("cosmdc2kpi", m_idxmdc, m_cosmdc2kpi);
344 status = m_tuple1->addIndexedItem("phimdc2kpi", m_idxmdc, m_phimdc2kpi);
345
346 status = m_tuple1->addIndexedItem("dedxpidjs", m_idxmdc, m_dedxpidjs);
347 status = m_tuple1->addIndexedItem("tof1pidjs", m_idxmdc, m_tof1pidjs);
348 status = m_tuple1->addIndexedItem("tof2pidjs", m_idxmdc, m_tof2pidjs);
349 status = m_tuple1->addIndexedItem("probejs", m_idxmdc, m_probejs);
350 status = m_tuple1->addIndexedItem("probmujs", m_idxmdc, m_probmujs);
351 status = m_tuple1->addIndexedItem("probpijs", m_idxmdc, m_probpijs);
352 status = m_tuple1->addIndexedItem("probkjs", m_idxmdc, m_probkjs);
353 status = m_tuple1->addIndexedItem("probprjs", m_idxmdc, m_probprjs);
354 status = m_tuple1->addItem ("mchic2kpi", m_mchic2kpi);
355 status = m_tuple1->addItem ("mpsip2kpi", m_mpsip2kpi);
356 status = m_tuple1->addItem ("chis2kpi", m_chis2kpi);
357 status = m_tuple1->addItem ("mchic4c2kpi", m_mchic4c2kpi);
358 status = m_tuple1->addItem ("mpsip4c2kpi", m_mpsip4c2kpi);
359 status = m_tuple1->addItem ("chis4c2kpi", m_chis4c2kpi);
360
361 status = m_tuple1->addItem("indexemc", m_idxemc, 0, 5000);
362 status = m_tuple1->addIndexedItem("numHits", m_idxemc, m_numHits);
363 status = m_tuple1->addIndexedItem("secmom", m_idxemc, m_secondmoment);
364 status = m_tuple1->addIndexedItem("latmom", m_idxemc, m_latmoment);
365 status = m_tuple1->addIndexedItem("timegm", m_idxemc, m_timegm);
366 status = m_tuple1->addIndexedItem("cellId", m_idxemc, m_cellId);
367 status = m_tuple1->addIndexedItem("module", m_idxemc, m_module);
368 status = m_tuple1->addIndexedItem("a20Moment", m_idxemc, m_a20Moment);
369 status = m_tuple1->addIndexedItem("a42Moment", m_idxemc, m_a42Moment);
370 status = m_tuple1->addIndexedItem("getEAll", m_idxemc, m_getEAll);
371 status = m_tuple1->addIndexedItem("getShowerId", m_idxemc, m_getShowerId);
372 status = m_tuple1->addIndexedItem("getClusterId", m_idxemc, m_getClusterId);
373 status = m_tuple1->addIndexedItem("x", m_idxemc, m_x);
374 status = m_tuple1->addIndexedItem("y", m_idxemc, m_y);
375 status = m_tuple1->addIndexedItem("z", m_idxemc, m_z);
376 status = m_tuple1->addIndexedItem("cosemc", m_idxemc, m_cosemc);
377 status = m_tuple1->addIndexedItem("phiemc", m_idxemc, m_phiemc);
378 status = m_tuple1->addIndexedItem("energy", m_idxemc, m_energy);
379 status = m_tuple1->addIndexedItem("e1", m_idxemc, m_eSeed);
380 status = m_tuple1->addIndexedItem("e9", m_idxemc, m_e3x3);
381 status = m_tuple1->addIndexedItem("e25", m_idxemc, m_e5x5);
382 status = m_tuple1->addIndexedItem("dang4c", m_idxemc, m_dang4c);
383 status = m_tuple1->addIndexedItem("dthe4c", m_idxemc, m_dthe4c);
384 status = m_tuple1->addIndexedItem("dphi4c", m_idxemc, m_dphi4c);
385 status = m_tuple1->addIndexedItem("dang4crt", m_idxemc, m_dang4crt);
386 status = m_tuple1->addIndexedItem("dthe4crt", m_idxemc, m_dthe4crt);
387 status = m_tuple1->addIndexedItem("dphi4crt", m_idxemc, m_dphi4crt);
388 status = m_tuple1->addIndexedItem("phtof", m_idxemc, 3, m_phgmtof,0.0,10000.0);
389 status = m_tuple1->addIndexedItem("phgmtof0", m_idxemc, m_phgmtof0);
390 status = m_tuple1->addIndexedItem("phgmtof1", m_idxemc, m_phgmtof1);
391 status = m_tuple1->addIndexedItem("phgmtof2", m_idxemc, m_phgmtof2);
392
393 }
394 else {
395 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
396 return StatusCode::FAILURE;
397 }
398 }
399}
400
401StatusCode sc;
402
403 if(m_skim4pi)
404 {
405 sc = createSubAlgorithm( "EventWriter", "Selectgam4pi", m_subalg1);
406 if( sc.isFailure() ) {
407 log << MSG::ERROR << "Error creating Sub-Algorithm Selectgam4pi" <<endreq;
408 return sc;
409 } else {
410 log << MSG::INFO << "Success creating Sub-Algorithm Selectgam4pi" <<endreq;
411 }
412 }
413
414
415 if(m_skim4k)
416 {
417 sc = createSubAlgorithm( "EventWriter", "Selectgam4k", m_subalg2);
418 if( sc.isFailure() ) {
419 log << MSG::ERROR << "Error creating Sub-Algorithm Selectgam4k" <<endreq;
420 return sc;
421 } else {
422 log << MSG::INFO << "Success creating Sub-Algorithm Selectgam4k" <<endreq;
423 }
424 }
425
426 if(m_skim2pi2pr)
427 {
428 sc = createSubAlgorithm( "EventWriter", "Selectgam2pi2pr", m_subalg3);
429 if( sc.isFailure() ) {
430 log << MSG::ERROR << "Error creating Sub-Algorithm Selectgam2pi2pr" <<endreq;
431 return sc;
432 } else {
433 log << MSG::INFO << "Success creating Sub-Algorithm Selectgam2pi2pr" <<endreq;
434 }
435 }
436
437
438
439
440 //
441 //--------end of book--------
442 //
443
444 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
445 return StatusCode::SUCCESS;
446
447}
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:55
The documentation for this class was generated from the following files: