BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
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 47 of file Gam4pikp.cxx.

48 : Algorithm( name, pSvcLocator ) {
49
50 // Declare the properties
51 declareProperty( "skim4pi", m_skim4pi = false );
52 declareProperty( "skim4k", m_skim4k = false );
53 declareProperty( "skim2pi2pr", m_skim2pi2pr = false );
54 declareProperty( "rootput", m_rootput = false );
55 declareProperty( "Ecms", m_ecms = 3.6862 );
56 declareProperty( "Vr0cut", m_vr0cut = 1.0 );
57 declareProperty( "Vz0cut", m_vz0cut = 5.0 );
58 declareProperty( "EnergyThreshold", m_energyThreshold = 0.05 );
59 declareProperty( "GammaPhiCut", m_gammaPhiCut = 20.0 );
60 declareProperty( "GammaThetaCut", m_gammaThetaCut = 20.0 );
61 declareProperty( "GammaDangCut", m_gammadangcut = 20.0 );
62 // declareProperty("Test4C", m_test4C = 1);
63 // declareProperty("Test5C", m_test5C = 1);
64 // declareProperty("CheckDedx", m_checkDedx = 1);
65 // declareProperty("CheckTof", m_checkTof = 1);
66}

Referenced by Gam4pikp().

Member Function Documentation

◆ execute()

StatusCode Gam4pikp::execute ( )

Definition at line 432 of file Gam4pikp.cxx.

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

◆ finalize()

StatusCode Gam4pikp::finalize ( )

Definition at line 2026 of file Gam4pikp.cxx.

2026 {
2027 // cout<<"total number: "<<Ncut[0]<<endl;
2028 // cout<<"nGood>=4, nCharge==0: "<<Ncut[1]<<endl;
2029 // cout<<"nGam>=1: "<<Ncut[2]<<endl;
2030 // cout<<"cgp>2 cgm>2: "<<Ncut[3]<<endl;
2031 // cout<<"2+ 2-: "<<Ncut[4]<<endl;
2032 // cout<<"all cg cycle, chisq<200: "<<Ncut[5]<<endl;
2033 // cout<<"4pi ok: "<<Ncut[6]<<endl;
2034 // cout<<"4k ok: "<<Ncut[7]<<endl;
2035 // cout<<"2pi 2p ok: "<<Ncut[8]<<endl;
2036 // cout<<"ntuple write: "<<Ncut[9]<<endl;
2037 MsgStream log( msgSvc(), name() );
2038 log << MSG::INFO << "in finalize()" << endmsg;
2039 return StatusCode::SUCCESS;
2040}

◆ initialize()

StatusCode Gam4pikp::initialize ( )

Definition at line 69 of file Gam4pikp.cxx.

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