rhopi.cxx

Go to the documentation of this file.

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/PropertyMgr.h"
 
#include "EventModel/EventModel.h"
#include "EventModel/Event.h"
 
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
#include "DstEvent/TofHitStatus.h"
#include "EventModel/EventHeader.h"
 
#include "TMath.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IHistogramSvc.h"
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Vector/TwoVector.h"
using CLHEP::Hep3Vector;
using CLHEP::Hep2Vector;
using CLHEP::HepLorentzVector;
#include "CLHEP/Geometry/Point3D.h"
#ifndef ENABLE_BACKWARDS_COMPATIBILITY
   typedef HepGeom::Point3D<double> HepPoint3D;
#endif
 
#include "VertexFit/KinematicFit.h"
#include "VertexFit/VertexFit.h"
#include "ParticleID/ParticleID.h"
 
#include "ValidJpsiRhopi/rhopi.h"
 
#include <vector>
 
const double mpi0c = 0.1349766;
const double mpi = 0.13957;
const double ecms = 3.097;
const double mk = 0.493677;
const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
//                          e         mu         pi      K         p
 
typedef std::vector<int> Vint;
typedef std::vector<HepLorentzVector> Vp4;
 
/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//   J/psi-->rho               +        pi                                 //
//             |-->pi pi                                                   //
//                                                                         //
//                   Final : (2Gamma + 2pi)            xugf 2008.09        //
///////////////////////////////////////////////////////////////////////////// 
DECLARE_COMPONENT(rhopi)
rhopi::rhopi(const std::string& name, ISvcLocator* pSvcLocator) :
  Algorithm(name, pSvcLocator) {
  
  m_tuple4 = 0;
  m_tuple5 = 0;
 
  for(int i = 0; i < 12; i++) m_pass[i] = 0;
 
//Declare the properties
  declareProperty("Vr0cut", m_vr0cut=5.0);
  declareProperty("Vz0cut", m_vz0cut=15.0);
  declareProperty("EnergyThreshold", m_energyThreshold=0.03);
  declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
  declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode rhopi::initialize(){
  MsgStream log(msgSvc(), name());
 
  log << MSG::INFO << "in initialize()" << endmsg;
  
  StatusCode status;
 
     status = service("THistSvc", m_thistsvc);
     if(status.isFailure() ){
          log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endreq;
          return status;
     }
     m_pi_vx = new TH1F( "pi_vx", "pi_vx", 100,-1.0,1.0 );
     status = m_thistsvc->regHist("/VAL/PHY/pi_vx", m_pi_vx);
     m_pi_vy = new TH1F( "pi_vy", "pi_vy", 100,-1.0,1.0 );
     status = m_thistsvc->regHist("/VAL/PHY/pi_vy", m_pi_vy);
     m_pi_vz = new TH1F( "pi_vz", "pi_vz", 100,-6.0,6.0 );
     status = m_thistsvc->regHist("/VAL/PHY/pi_vz", m_pi_vz);
     m_pip_mom = new TH1F( "pip_mom", "pip_moment", 100,0.0,1.6 );
     status = m_thistsvc->regHist("/VAL/PHY/pip_mom", m_pip_mom);
     m_pim_mom = new TH1F( "pim_mom", "pim_moment", 100,0.0,1.6 );
     status = m_thistsvc->regHist("/VAL/PHY/pim_mom", m_pim_mom);
     m_rhop_mass = new TH1F( "rhop_mass", "rhop_mass", 100,0.0,1.5 );
     status = m_thistsvc->regHist("/VAL/PHY/rhop_mass", m_rhop_mass);
     m_rhom_mass = new TH1F( "rhom_mass", "rhom_mass", 100,0.0,1.5 );
     status = m_thistsvc->regHist("/VAL/PHY/rhom_mass", m_rhom_mass);
     m_rho0_mass = new TH1F( "rho0_mass", "rho0_mass", 100,0.0,1.5 );
     status = m_thistsvc->regHist("/VAL/PHY/rho0_mass", m_rho0_mass);
     m_pi0_mass = new TH1F( "pi0_mass", "pi0_mass", 100,0.0,0.5 );
     status = m_thistsvc->regHist("/VAL/PHY/pi0_mass", m_pi0_mass);
     m_chisq_4c = new TH1F( "chisq_4c", "chisq_4c", 100,0.0,150. );
     status = m_thistsvc->regHist("/VAL/PHY/chisq_4c", m_chisq_4c);
     m_chisq_5c = new TH1F( "chisq_5c", "chisq_5c", 100,0.0,150. );
     status = m_thistsvc->regHist("/VAL/PHY/chisq_5c", m_chisq_5c);
     
     m_cos_pip = new TH1F( "cos_pip", "cos_pip", 100, -1.0,1.0);
     status = m_thistsvc->regHist("/VAL/PHY/cos_pip", m_cos_pip);
     m_cos_pim = new TH1F( "cos_pim", "cos_pim", 100, -1.0,1.0);
     status = m_thistsvc->regHist("/VAL/PHY/cos_pim", m_cos_pim);
 
     m_chipi_dedx = new TH1F( "chipi_dedx", "chipi_dedx", 60,-5.0,10. );
     status = m_thistsvc->regHist("/VAL/PHY/chipi_dedx", m_chipi_dedx);
     m_chie_dedx = new TH1F( "chie_dedx", "chie_dedx", 60,-15.0,5. );
     status = m_thistsvc->regHist("/VAL/PHY/chie_dedx", m_chie_dedx);
     m_chimu_dedx = new TH1F( "chimu_dedx", "chimu_dedx", 60,-5.0,10. );
     status = m_thistsvc->regHist("/VAL/PHY/chimu_dedx", m_chimu_dedx);
     m_chik_dedx = new TH1F( "chik_dedx", "chik_dedx", 100,-20.0,10. );
     status = m_thistsvc->regHist("/VAL/PHY/chik_dedx", m_chik_dedx);
     m_chip_dedx = new TH1F( "chip_dedx", "chip_dedx", 100,-20.0,5. );
     status = m_thistsvc->regHist("/VAL/PHY/chip_dedx", m_chip_dedx);
 
    NTuplePtr nt4(ntupleSvc(), "FILE1/h4");
    if ( nt4 ) m_tuple4 = nt4;
    else {
      m_tuple4 = ntupleSvc()->book ("FILE1/h4", CLID_ColumnWiseTuple, "4gam6pi Ntuple");
      if ( m_tuple4 )    {
    status = m_tuple4->addItem ("ngam",           m_ngoodneu);
    status = m_tuple4->addItem ("npip",           m_npip);
    status = m_tuple4->addItem ("npim",           m_npim);
    status = m_tuple4->addItem ("ngoodch",        m_ngoodch);
    status = m_tuple4->addItem ("chisq4c",        m_chisq4c);
    status = m_tuple4->addItem ("ppi04c",         m_ppi0);
    status = m_tuple4->addItem ("mpi04c",         m_mpi0);
    status = m_tuple4->addItem ("chisq5c",        m_chisq5c);
    status = m_tuple4->addItem ("ppi05c",         m_ppi0fit);
    status = m_tuple4->addItem ("mpi05c",         m_mpi0fit);
    status = m_tuple4->addItem ("g1inpi0the",     m_g1inpi0the);
    status = m_tuple4->addItem ("g2inpi0the",     m_g2inpi0the);
    status = m_tuple4->addItem ("theta2pi",       m_theta2pi);
    status = m_tuple4->addItem ("ppip",           m_ppip);
    status = m_tuple4->addItem ("ppim",           m_ppim);
    status = m_tuple4->addItem ("p2pi",           m_p2pi);
    status = m_tuple4->addItem ("m2pi",           m_m2pi);
    status = m_tuple4->addItem ("ppip0",          m_ppip0);
    status = m_tuple4->addItem ("mpip0",          m_mpip0);
    status = m_tuple4->addItem ("ppim0",          m_ppim0);
    status = m_tuple4->addItem ("mpim0",          m_mpim0);
    status = m_tuple4->addItem ("eneumiss",       m_eneumiss);
    status = m_tuple4->addItem ("pneumiss",       m_pneumiss);
    status = m_tuple4->addItem ("mneumiss",       m_mneumiss);
//  status = m_tuple4->addIndexedItem ("kal_cos",  m_ngoodch , m_kal_cos);
      }
      else    { 
    log << MSG::ERROR << "    Cannot book N-tuple:" << long(m_tuple4) << endmsg;
    return StatusCode::FAILURE;
      }
    }
 
    NTuplePtr nt5(ntupleSvc(), "FILE1/h5");
    if ( nt5 ) m_tuple5 = nt5;
    else {
      m_tuple5 = ntupleSvc()->book ("FILE1/h5", CLID_ColumnWiseTuple, "4gam6pi Ntuple");
      if ( m_tuple5 )    {
    status = m_tuple5->addItem ("m2graw",         m_m2graw);
    status = m_tuple5->addItem ("emiss",          m_emiss);
    status = m_tuple5->addItem ("pmiss",          m_pmiss);
    status = m_tuple5->addItem ("mmiss",          m_mmiss);
    status = m_tuple5->addItem ("prho0",          m_prho0);
    status = m_tuple5->addItem ("mrho0",          m_mrho0);
    status = m_tuple5->addItem ("pmrho0",         m_pmrho0);
    status = m_tuple5->addItem ("mmrho0",         m_mmrho0);
    status = m_tuple5->addItem ("prhop",          m_prhop);
    status = m_tuple5->addItem ("mrhop",          m_mrhop);
    status = m_tuple5->addItem ("pmrhom",         m_pmrhom);
    status = m_tuple5->addItem ("mmrhom",         m_mmrhom);
    status = m_tuple5->addItem ("prhom",          m_prhom);
    status = m_tuple5->addItem ("ppipraw",        m_ppipraw);
    status = m_tuple5->addItem ("mrhom",          m_mrhom);
 
      }
      else    { 
    log << MSG::ERROR << "    Cannot book N-tuple:" << long(m_tuple4) << endmsg;
    return StatusCode::FAILURE;
      }
    }
 
 
  //
  //--------end of book--------
  //
 
  log << MSG::INFO << "successfully return from initialize()" <<endmsg;
  return StatusCode::SUCCESS;
 
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode rhopi::execute() {
  
  StatusCode sc = StatusCode::SUCCESS;
 
  MsgStream log(msgSvc(), name());
//  log << MSG::INFO << "in execute()" << endreq;
 
  m_pass[0] += 1;
 
  SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
  //SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::Recon::EvtRecEvent);
  SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
  
//  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::Recon::EvtRecTrackCol);
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(),  EventModel::EvtRec::EvtRecTrackCol);
  //
  // check x0, y0, z0, r0
  // suggest cut: |z0|<10 && r0<5
  //
 
  Vint ipip, ipim, iGood;
  iGood.clear();
  ipip.clear();
  ipim.clear();
 
  Vp4 ppip, ppim;
  ppip.clear();
  ppim.clear();
 
  int TotCharge = 0;
  for(int i = 0; i < evtRecEvent->totalCharged(); i++){
    EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
    if(!(*itTrk)->isMdcTrackValid()) continue;
    RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
 
     m_pi_vx->Fill(mdcTrk->x());
     m_pi_vy->Fill(mdcTrk->y());
     m_pi_vz->Fill(mdcTrk->z());
 
    if(fabs(mdcTrk->z()) >= m_vz0cut) continue;
    if(mdcTrk->r() >= m_vr0cut) continue;
    iGood.push_back((*itTrk)->trackId());
    TotCharge += mdcTrk->charge();
  }
  //
  // Finish Good Charged Track Selection
  //
  int nGood = iGood.size();
    m_ngoodch = nGood;
 
  Vint iGam;
  iGam.clear();
  for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
    EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
    if(!(*itTrk)->isEmcShowerValid()) continue;
    RecEmcShower *emcTrk = (*itTrk)->emcShower();
    Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
    // find the nearest charged track
    double dthe = 200.;
    double dphi = 200.;
    double dang = 200.;
    for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
      EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
      if(!(*jtTrk)->isExtTrackValid()) continue;
      RecExtTrack *extTrk = (*jtTrk)->extTrack();
      if(extTrk->emcVolumeNumber() == -1) continue;
      Hep3Vector extpos = extTrk->emcPosition();
      //      double ctht = extpos.cosTheta(emcpos);
      double angd = extpos.angle(emcpos);
      double thed = extpos.theta() - emcpos.theta();
      double phid = extpos.deltaPhi(emcpos);
      thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
      phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
      angd = fmod(angd, CLHEP::twopi);
 
      if(angd < dang){
              dang = angd;
              dthe = thed;
              dphi = phid;
        }
    }
 
    if(dang >= 200.) continue;
 
    double eraw = emcTrk->energy();
    dthe = dthe * 180 / (CLHEP::pi);
    dphi = dphi * 180 / (CLHEP::pi);
    dang = dang * 180 / (CLHEP::pi);
    if(eraw < m_energyThreshold) continue;
//    if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
    //
    // good photon cut will be set here
    //
    iGam.push_back((*itTrk)->trackId());
  }
  //
  // Finish Good Photon Selection
  //
  int nGam = iGam.size();
    m_ngoodneu = nGam;
  //
  // Assign 4-momentum to each photon
  //
 
       HepLorentzVector ptcharg, ptneu, ptchargp, ptchargm;
       Vp4 pGam;
       pGam.clear();
  for(int i = 0; i < nGam; i++) {
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
    RecEmcShower* emcTrk = (*itTrk)->emcShower();
    double eraw = emcTrk->energy();
    double phi = emcTrk->phi();
    double the = emcTrk->theta();
    HepLorentzVector ptrk;
    ptrk.setPx(eraw*sin(the)*cos(phi));
    ptrk.setPy(eraw*sin(the)*sin(phi));
    ptrk.setPz(eraw*cos(the));
    ptrk.setE(eraw);
    pGam.push_back(ptrk);
  }
 
//
// Dedx Check ================
//
    for(int i = 0; i < nGood; i++) {
      EvtRecTrackIterator  itTrk = evtRecTrkCol->begin() + iGood[i];
      if(!(*itTrk)->isMdcTrackValid()) continue;
      if(!(*itTrk)->isMdcDedxValid())continue;
      RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
      RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
      m_chie_dedx -> Fill(dedxTrk->chiE());
      m_chimu_dedx -> Fill(dedxTrk->chiMu());
      m_chipi_dedx -> Fill(dedxTrk->chiPi());
      m_chik_dedx -> Fill(dedxTrk->chiK());
      m_chip_dedx -> Fill(dedxTrk->chiP());
    }
 
  //
  // Assign 4-momentum to each charged track
  //
 
     int TotQ = 0;
 
  ParticleID *pid = ParticleID::instance();
  for(int i = 0; i < nGood; i++) {
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
    //    if(pid) delete pid;
    pid->init();
    pid->setMethod(pid->methodProbability());
    pid->setChiMinCut(4);
    pid->setRecTrack(*itTrk);
    pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID sub-system
    pid->identify(pid->onlyPion() | pid->onlyKaon());    // seperater Pion/Kaon
    pid->calculate();
    if(!(pid->IsPidInfoValid())) continue;
 
    if(!(*itTrk)->isMdcKalTrackValid()) continue ;
    RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
    RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
 
      TotQ += mdcKalTrk->charge();
 
//    RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
 
//    if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
    if(pid->probPion() < pid->probKaon()) continue;
      mdcKalTrk->setPidType(RecMdcKalTrack::pion);
      HepLorentzVector ptrk;
      ptrk.setPx(mdcKalTrk->px());
      ptrk.setPy(mdcKalTrk->py());
      ptrk.setPz(mdcKalTrk->pz());
      double p3 = ptrk.mag();
      ptrk.setE(sqrt(p3*p3+mpi*mpi));
    if(mdcTrk->charge() >0) {
      ipip.push_back(iGood[i]);
      ppip.push_back(ptrk);
    } else {
      ipim.push_back(iGood[i]);
      ppim.push_back(ptrk);
    }
  }
 
  int npip = ipip.size();
  int npim = ipim.size();
     m_npip = npip;
     m_npim = npim;
 
  if(nGood != 2 || TotCharge!=0) return sc;
  m_pass[1] += 1;
  if(nGam<2 ) return sc;
  m_pass[2] += 1;
 
//  if(npip != 1) return sc;
//  if(npip*npim != 1) return sc;
 
        HepLorentzVector BoostP(0.0,0.0,0.0,ecms);
           BoostP[0] = 2*sin(0.011)*(ecms/2);
        Hep3Vector u_Boost = -BoostP.boostVector();
 
  if(npip == 1) {
//
//---- Find the Best Pi0  -----
//
  HepLorentzVector p2g, ppi0;
  double delmpi = 999.;
  int ig11=-1, ig21=-1;
 
  for(int i = 0; i < nGam - 1; i++){
    for(int j = i+1; j < nGam; j++) {
      HepLorentzVector ppi0  = pGam[i] + pGam[j];
      if(fabs(ppi0.m()-mpi0c) > delmpi) continue;
      if(ppi0.m() > 0.2) continue;
      delmpi = fabs(ppi0.m()-mpi0c);
      ig11 = i;
      ig21 = j;
    }
  }
     if(ig11==-1 || ig21== -1) return sc;
     p2g = pGam[ig11] + pGam[ig21];
         m_m2graw = p2g.m();
 
         HepLorentzVector Ppip1 = ppip[0];
           Ppip1.boost(u_Boost);
           p2g.boost(u_Boost);
//
// ----- Pion Efficiency -----
//
//
//------ Missing pi- 
//
    HepLorentzVector Pmiss = -Ppip1 - p2g;
     m_pmiss = Pmiss.rho();
        double emiss = ecms - Ppip1.e() - p2g.e();
     m_emiss = emiss;
     m_mmiss = sqrt(fabs(emiss*emiss - Pmiss.rho()*Pmiss.rho()));
 
  HepLorentzVector Prho0, Prhop, Prhom, pmisspi, ptot;
     pmisspi.setPx(Pmiss.px());
     pmisspi.setPy(Pmiss.py());
     pmisspi.setPz(Pmiss.pz());
     pmisspi.setE(emiss);
 
      Prho0 = Ppip1 + pmisspi;
      Prhop = Ppip1 + p2g;
      Prhom = pmisspi + p2g;
      ptot = Ppip1 + pmisspi + p2g;
 
     m_pmrho0 = Prho0.rho();
     m_mmrho0 = Prho0.m();
//     m_prho0 = (Ppip1 + Ppim1).rho();
//     m_mrho0 = (Ppip1 + Ppim1).m();
     m_prhop = Prhop.rho();
     m_mrhop = Prhop.m();
//     m_prhom = (p2g + Ppim1).rho();
//     m_mrhom = (p2g + Ppim1).m();
     m_pmrhom = Prhom.rho();
     m_mmrhom = Prhom.m();
     m_ppipraw = Ppip1.rho();
 
           m_tuple5->write();
    }
 
  if(npip*npim != 1) return sc;
  m_pass[3] += 1;
 
         HepLorentzVector Ppim1 = ppim[0];
           Ppim1.boost(u_Boost);
         HepLorentzVector Ppip1 = ppip[0];
           Ppip1.boost(u_Boost);
 
  HepLorentzVector Pneumiss = -(Ppim1 + Ppip1);  
     m_pneumiss = Pneumiss.rho();
     double eneumiss = ecms - Ppim1.e() - Ppip1.e();
     m_eneumiss = eneumiss;
     m_mneumiss = sqrt(fabs(eneumiss*eneumiss - Pneumiss.rho()*Pneumiss.rho()));
 
//
//  Let's play vertex fit here
//
 HepPoint3D vx(0., 0., 0.);
  HepSymMatrix Evx(3, 0);
  double bx = 1E+6;
  double by = 1E+6;
  double bz = 1E+6;
  Evx[0][0] = bx*bx;
  Evx[1][1] = by*by;
  Evx[2][2] = bz*bz;
 
  VertexParameter vxpar;
  vxpar.setVx(vx);
  vxpar.setEvx(Evx);
 
  VertexFit* vtxfit = VertexFit::instance();
  KinematicFit * kmfit = KinematicFit::instance();
 
  RecMdcKalTrack *pipTrk = (*(evtRecTrkCol->begin()+ipip[0]))->mdcKalTrack();
  pipTrk->setPidType(RecMdcKalTrack::pion);
  WTrackParameter wpip0(mpi, pipTrk->helix(), pipTrk->err());
 
  pipTrk->setPidType(RecMdcKalTrack::kaon);
  WTrackParameter wkp0(mk, pipTrk->helix(), pipTrk->err());
 
  RecMdcKalTrack *pimTrk = (*(evtRecTrkCol->begin()+ipim[0]))->mdcKalTrack();
  pimTrk->setPidType(RecMdcKalTrack::pion);
  WTrackParameter wpim0(mpi, pimTrk->helix(), pimTrk->err());
 
  pimTrk->setPidType(RecMdcKalTrack::kaon);
  WTrackParameter wkm0(mk, pimTrk->helix(), pimTrk->err());
 
 
  vtxfit->init();
  vtxfit->AddTrack(0, wpip0);
  vtxfit->AddTrack(1, wpim0);
  vtxfit->AddVertex(0, vxpar,0, 1);
  if(!vtxfit->Fit(0)) return sc;
 
 
 
  m_pass[4] += 1;
 
  WTrackParameter wpip = vtxfit->wtrk(0);
  WTrackParameter wpim = vtxfit->wtrk(1);
 
//
//  Apply Kinematic 4C fit
// 
 
    double chisq4c = 9999.;
 
    int ig1 = -1; 
    int ig2 = -1;
 
    for(int i = 0; i < nGam-1; i++) {
     RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
      for(int j = i+1; j < nGam; j++) {
        RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
 
    kmfit->init();
    kmfit->AddTrack(0, wpip);
    kmfit->AddTrack(1, wpim);
        kmfit->AddTrack(2, 0.0, g1Trk);
        kmfit->AddTrack(3, 0.0, g2Trk);
    kmfit->AddFourMomentum(0, BoostP);
 
    if(!kmfit->Fit()) continue;
          double chi2 = kmfit->chisq();
          if(chi2 > chisq4c) continue;
            chisq4c = chi2;
            ig1 = i;
            ig2 = j;
        }    // gamma1
      }      // gamma2
 
 
  if(chisq4c > 500 || ig1 < 0) return sc; 
      RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[ig1]))->emcShower();
      RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[ig2]))->emcShower();
 
    kmfit->init();
    kmfit->AddTrack(0, wpip);
    kmfit->AddTrack(1, wpim);
        kmfit->AddTrack(2, 0.0, g1Trk);
        kmfit->AddTrack(3, 0.0, g2Trk);
    kmfit->AddFourMomentum(0, BoostP);
    if(!kmfit->Fit()) return sc;
          double chi = kmfit->chisq();
            m_chisq4c = chi;
       m_chisq_4c -> Fill(chi);
 
       HepLorentzVector p2gam = kmfit->pfit(2) + kmfit->pfit(3);
       m_ppi0=p2gam.rho();
       m_mpi0=p2gam.m();
 
       m_pi0_mass -> Fill(p2gam.m());
 
  m_pass[5] += 1;
 
//
//  Apply Kinematic 5C fit
// 
//
// 5c for J/psi-->pi+pi-pi0
//
    kmfit->init();
    kmfit->AddTrack(0, wpip);
    kmfit->AddTrack(1, wpim);
        kmfit->AddTrack(2, 0.0, g1Trk);
        kmfit->AddTrack(3, 0.0, g2Trk);
        kmfit->AddResonance(0,mpi0c,2,3);
    kmfit->AddFourMomentum(1, BoostP);
    if(!kmfit->Fit()) return sc;
          double chis = kmfit->chisq();
 
    m_pass[6] += 1;
//
// 5c for J/psi-->K+K-pi0
//
     double chisk5c = 9999.;  
 
      vtxfit->init();
      vtxfit->AddTrack(0, wkp0);
      vtxfit->AddTrack(1, wkm0);
      vtxfit->AddVertex(0, vxpar,0, 1);
 
      if(vtxfit->Fit(0)){
 
      WTrackParameter wkp = vtxfit->wtrk(0);
      WTrackParameter wkm = vtxfit->wtrk(1);
 
    kmfit->init();
    kmfit->AddTrack(0, wkp);
    kmfit->AddTrack(1, wkm);
        kmfit->AddTrack(2, 0.0, g1Trk);
        kmfit->AddTrack(3, 0.0, g2Trk);
        kmfit->AddResonance(0,mpi0c,2,3);
    kmfit->AddFourMomentum(1, BoostP);
    if(kmfit->Fit()) chisk5c = kmfit->chisq();
      }
      if(chisk5c<chis) return sc;
 
      m_pass[7] += 1;
      //
      // 5c Fit again for J/psi-->pi+pi-pi0
      //
      kmfit->init();
      kmfit->AddTrack(0, wpip);
      kmfit->AddTrack(1, wpim);
      kmfit->AddTrack(2, 0.0, g1Trk);
      kmfit->AddTrack(3, 0.0, g2Trk);
      kmfit->AddResonance(0,mpi0c,2,3);
      kmfit->AddFourMomentum(1, BoostP);
      if(!kmfit->Fit(0)) return sc;
      if(!kmfit->Fit()) return sc;
 
      m_pass[8] += 1;
      m_chisq5c = kmfit->chisq();
 
      m_chisq_5c -> Fill(kmfit->chisq());
 
      HepLorentzVector ppi1 = kmfit->pfit(0);
      HepLorentzVector ppi2 = kmfit->pfit(1);
      HepLorentzVector pgam1 = kmfit->pfit(2);
      HepLorentzVector pgam2 = kmfit->pfit(3);
      HepLorentzVector p2gam1 = kmfit->pfit(2) + kmfit->pfit(3);
      HepLorentzVector p2pi = kmfit->pfit(0) + kmfit->pfit(1);
      m_ppi0fit = p2gam1.rho();
      m_mpi0fit = p2gam1.m();
      m_ppip = ppi1.rho();
      m_ppim = ppi2.rho();
      m_p2pi = p2pi.rho();
      m_m2pi = p2pi.m();
      m_ppip0 = (ppi1+p2gam1).rho();
      m_mpip0 = (ppi1+p2gam1).m();
      m_ppim0 = (ppi2+p2gam1).rho();
      m_mpim0 = (ppi2+p2gam1).m();
 
      m_pip_mom -> Fill(ppi1.rho());
      m_pim_mom -> Fill(ppi2.rho());
      m_rho0_mass -> Fill(p2pi.m());
      m_rhop_mass -> Fill((ppi1+p2gam1).m());
      m_rhom_mass -> Fill((ppi2+p2gam1).m());
 
      //--------------Angle between two charged pions ---------
      double theta2pi= ppi1.px()*ppi2.px()+ppi1.py()*ppi2.py()+ppi1.pz()*ppi2.pz();
      double time2pi = ppi1.rho()*ppi2.rho();
      if(time2pi == 0.) return sc;
 
      m_pass[9] += 1;
 
      theta2pi = theta2pi/time2pi;
 
      //         cout <<"time2pi=====" << time2pi <<endl;
 
      if(theta2pi <= -1.) theta2pi=180.;
      else if(theta2pi >= 1.) theta2pi=0.;
      else if(fabs(theta2pi) < 1.) theta2pi = acos(theta2pi)*180./3.14159;
      m_theta2pi=theta2pi;
 
      //----- Theta angle of gamma in pi0 system --------
      HepLorentzVector pg1inpi0, pg2inpi0;
 
      //       HepLorentzVector betapi0 = p2gam1/(p2gam1.e());
      double betapi0 = p2gam1.beta();
      //    ==================================  
      Hep3Vector twogam(p2gam.px(),p2gam.py(),p2gam.pz());
      //       CLHEP::boostOf(pg1inpi0,gam1,betapi0);
      m_g1inpi0the = ((boostOf(pgam1,twogam/p2gam.rho(),betapi0)).theta())*180./3.14159;
      m_g2inpi0the = ((boostOf(pgam2,twogam/p2gam.rho(),betapi0)).theta())*180./3.14159;
      /*
      //
      // ----- Pion Efficiency -----
      //
      HepLorentzVector p2g = pGam[ig1] + pGam[ig2];
      m_m2graw = p2g.m();
      p2g.boost(u_Boost);
      //
      //------ Missing pi- 
      //
      HepLorentzVector Pmiss = -Ppip1 - p2g;
      m_pmiss = Pmiss.rho();
      double emiss = ecms - Ppip1.e() - p2g.e();
      m_emiss = emiss;
      m_mmiss = sqrt(fabs(emiss*emiss - Pmiss.rho()*Pmiss.rho()));
 
      HepLorentzVector Prho0, Prhop, Prhom, pmisspi, ptot;
      pmisspi.setPx(Pmiss.px());
      pmisspi.setPy(Pmiss.py());
      pmisspi.setPz(Pmiss.pz());
      pmisspi.setE(emiss);
 
      Prho0 = Ppip1 + pmisspi;
      Prhop = Ppip1 + p2g;
      Prhom = pmisspi + p2g;
      ptot = Ppip1 + pmisspi + p2g;
 
      m_pmrho0 = Prho0.rho();
      m_mmrho0 = Prho0.m();
      m_prho0 = (Ppip1 + Ppim1).rho();
      m_mrho0 = (Ppip1 + Ppim1).m();
      m_prhop = Prhop.rho();
      m_mrhop = Prhop.m();
      m_prhom = (p2g + Ppim1).rho();
      m_mrhom = (p2g + Ppim1).m();
      m_pmrhom = Prhom.rho();
      m_mmrhom = Prhom.m();
       */
      m_tuple4->write();
 
      return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode rhopi::finalize() {
 
    MsgStream log(msgSvc(), name());
    log << MSG::ALWAYS << "Total Entries :        " << m_pass[0] << endreq;
    log << MSG::ALWAYS << "Ncharge=2 Cut :        " << m_pass[1] << endreq;
    log << MSG::ALWAYS << "Ngam<2 cut :           " << m_pass[2] << endreq;
    log << MSG::ALWAYS << "Nch+=1,Nch-=1 cut :    " << m_pass[3] << endreq;
    log << MSG::ALWAYS << "Vertex Fit :           " << m_pass[4] << endreq;
    log << MSG::ALWAYS << "4c-Fit :               " << m_pass[5] << endreq;
    log << MSG::ALWAYS << "5c-Fit :               " << m_pass[6] << endreq;
    log << MSG::ALWAYS << "chi3pi<chikkpi cut :   " << m_pass[7] << endreq;
    log << MSG::ALWAYS << "5c-Fit Again:          " << m_pass[8] << endreq;
    log << MSG::ALWAYS << "Final :                " << m_pass[9] << endreq;
    return StatusCode::SUCCESS;
}