DQARhopi.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 "VertexFit/IVertexDbSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/ISvcLocator.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 "McTruth/McParticle.h"
 
#include "TH1F.h"
#include "TMath.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"
#include "GaudiKernel/ITHistSvc.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 "DQARhopiAlg/DQARhopi.h"
 
#include "VertexFit/KinematicFit.h"
#include "VertexFit/VertexFit.h"
#include "ParticleID/ParticleID.h"
 
#include <vector>
using namespace Event;
//const double twopi = 6.2831853;
//const double pi = 3.1415927;
const double mpi = 0.13957;
const double mk  = 0.493677;
const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
//const double velc = 29.9792458;  tof_path unit in cm.
const double velc = 299.792458;   // tof path unit in mm
typedef std::vector<int> Vint;
typedef std::vector<HepLorentzVector> Vp4;
 
const HepLorentzVector ecms(0.034,0,0,3.097);
const Hep3Vector u_cms = -ecms.boostVector();
 
 
/////////////////////////////////////////////////////////////////////////////
 
DQARhopi::DQARhopi(const std::string& name, ISvcLocator* pSvcLocator) :
  Algorithm(name, pSvcLocator) {
  
  //Declare the properties  
  declareProperty("Vr0cut", m_vr0cut=1.0);
  declareProperty("Vz0cut", m_vz0cut=5.0);
  declareProperty("Vctcut", m_cthcut=0.93);
  declareProperty("EnergyThreshold", m_energyThreshold=0.05);
  declareProperty("GammaAngCut", m_gammaAngCut=25.0);
  declareProperty("Test4C", m_test4C = 1);
  declareProperty("Test5C", m_test5C = 1);
  declareProperty("CheckDedx", m_checkDedx = 1);
  declareProperty("CheckTof",  m_checkTof = 1);
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode DQARhopi::initialize(){
  MsgStream log(msgSvc(), name());
 
  log << MSG::INFO << "in initialize()" << endmsg;
 
  Ncut0=Ncut1=Ncut2=Ncut3=Ncut4=Ncut5=Ncut6=Ncut7=Ncut8=Ncut9=Ncut10=0;
  
  StatusCode status;
 
    NTuplePtr nt4(ntupleSvc(), "DQAFILE/Rhopi");
    if ( nt4 ) m_tuple4 = nt4;
    else {
      m_tuple4 = ntupleSvc()->book ("DQAFILE/Rhopi", CLID_ColumnWiseTuple, "ks N-Tuple example");
      if ( m_tuple4 )    {
        status = m_tuple4->addItem ("run",    m_run);
        status = m_tuple4->addItem ("rec",    m_rec);     
        status = m_tuple4->addItem ("nch",    m_nch);
        status = m_tuple4->addItem ("nneu",   m_nneu);
    status = m_tuple4->addItem ("chi1",   m_chi1);
    status = m_tuple4->addItem ("mpi0",   m_mpi0);
        status = m_tuple4->addItem ("mprho0", m_prho0);
        status = m_tuple4->addItem ("mprhop", m_prhop);
        status = m_tuple4->addItem ("mprhom", m_prhom);
        status = m_tuple4->addItem ("mgood",  m_good);
        status = m_tuple4->addItem ("mgam",   m_gam);
        status = m_tuple4->addItem ("mpip",   m_pip);
        status = m_tuple4->addItem ("mpim",   m_pim);
        status = m_tuple4->addItem ("m2gam",  m_2gam);
        status = m_tuple4->addItem ("ngch",     m_ngch,  0, 10);
        status = m_tuple4->addItem ("nggneu",   m_nggneu,0, 10);
        status = m_tuple4->addItem ("moutpi0",m_outpi0);
        status = m_tuple4->addItem ("cosang", m_cosang);
        status = m_tuple4->addItem ("moutpip",m_outpip);
        status = m_tuple4->addItem ("moutpim",m_outpim);
        status = m_tuple4->addItem ("menpip", m_enpip);
        status = m_tuple4->addItem ("mdcpip", m_dcpip );
        status = m_tuple4->addItem ("menpim", m_enpim );
        status = m_tuple4->addItem ("mdcpim", m_dcpim );   
        status = m_tuple4->addItem ("mpipf",  m_pipf);        
        status = m_tuple4->addItem ("mpimf",  m_pimf);        
        status = m_tuple4->addItem ("mpi0f",  m_pi0f);
 
        status = m_tuple4->addItem ("mpmax",   m_pmax);
        status = m_tuple4->addItem ("mppx",    m_ppx);
        status = m_tuple4->addItem ("mppy",    m_ppy);
        status = m_tuple4->addItem ("mppz",    m_ppz);
        status = m_tuple4->addItem ("mcosthep",m_costhep);
        status = m_tuple4->addItem ("mppxkal", m_ppxkal);
        status = m_tuple4->addItem ("mppykal", m_ppykal);
        status = m_tuple4->addItem ("mppzkal", m_ppzkal);
        status = m_tuple4->addItem ("mmpx",    m_mpx);
        status = m_tuple4->addItem ("mmpy",    m_mpy);
        status = m_tuple4->addItem ("mmpz",    m_mpz);
        status = m_tuple4->addItem ("mcosthem",m_costhem);
        status = m_tuple4->addItem ("mmpxkal", m_mpxkal);
        status = m_tuple4->addItem ("mmpykal", m_mpykal);
        status = m_tuple4->addItem ("mmpzkal", m_mpzkal);
 
        status = m_tuple4->addItem ("mvxpin",  m_vxpin);
        status = m_tuple4->addItem ("mvypin",  m_vypin);
        status = m_tuple4->addItem ("mvzpin",  m_vzpin);
        status = m_tuple4->addItem ("mvrpin",  m_vrpin);
        status = m_tuple4->addItem ("mcosthepin",m_costhepin);
        status = m_tuple4->addItem ("mvxmin",  m_vxmin);
        status = m_tuple4->addItem ("mvymin",  m_vymin);
        status = m_tuple4->addItem ("mvzmin",  m_vzmin);
        status = m_tuple4->addItem ("mvrmin",  m_vrmin);
        status = m_tuple4->addItem ("mcosthemin",m_costhemin);
        status = m_tuple4->addItem ("mvxp",    m_vxp);
        status = m_tuple4->addItem ("mvyp",    m_vyp);
        status = m_tuple4->addItem ("mvzp",    m_vzp);
        status = m_tuple4->addItem ("mvrp",    m_vrp);
        status = m_tuple4->addItem ("mvxm",    m_vxm);
        status = m_tuple4->addItem ("mvym",    m_vym);
        status = m_tuple4->addItem ("mvzm",    m_vzm);
        status = m_tuple4->addItem ("mvrm",    m_vrm);
        status = m_tuple4->addItem ("nangecc",m_nangecc,0,10);
        status = m_tuple4->addIndexedItem ("mdthec", m_nangecc, m_dthec);
        status = m_tuple4->addIndexedItem ("mdphic", m_nangecc, m_dphic);
        status = m_tuple4->addIndexedItem ("mdangc", m_nangecc, m_dangc);
        status = m_tuple4->addIndexedItem ("mspippim", m_nangecc, m_mspippim);
 
        status = m_tuple4->addItem ("angsg",   dangsg);
        status = m_tuple4->addItem ("thesg",   dthesg);
        status = m_tuple4->addItem ("phisg",   dphisg);
        status = m_tuple4->addItem ("cosgth1", cosgth1);
        status = m_tuple4->addItem ("cosgth2", cosgth2);
 
        status = m_tuple4->addItem ("mchi5",   m_chi5);
        status = m_tuple4->addItem ("mkpi0",   m_kpi0);
        status = m_tuple4->addItem ("mkpkm",   m_kpkm);
        status = m_tuple4->addItem ("mkpp0",   m_kpp0);
        status = m_tuple4->addItem ("mkmp0",   m_kmp0);            
        status = m_tuple4->addItem ("mpgam2pi1",m_pgam2pi1);
        status = m_tuple4->addItem ("mpgam2pi2",m_pgam2pi2);
        status = m_tuple4->addItem ("cosva1",   cosva1);
        status = m_tuple4->addItem ("cosva2",   cosva2);
        status = m_tuple4->addItem ("laypi1",   m_laypi1);
        status = m_tuple4->addItem ("hit1",     m_hit1);
        status = m_tuple4->addItem ("laypi2",   m_laypi2);
        status = m_tuple4->addItem ("hit2",     m_hit2);
        status = m_tuple4->addItem ("anglepm",  m_anglepm);
 
        status = m_tuple4->addIndexedItem ("mptrk", m_ngch, m_ptrk);
        status = m_tuple4->addIndexedItem ("chie",   m_ngch, m_chie);
        status = m_tuple4->addIndexedItem ("chimu",  m_ngch,m_chimu);
        status = m_tuple4->addIndexedItem ("chipi",  m_ngch,m_chipi);
        status = m_tuple4->addIndexedItem ("chik",   m_ngch,m_chik);
        status = m_tuple4->addIndexedItem ("chip",   m_ngch,m_chip);
        status = m_tuple4->addIndexedItem ("probPH", m_ngch,m_probPH);
        status = m_tuple4->addIndexedItem ("normPH", m_ngch,m_normPH);
        status = m_tuple4->addIndexedItem ("ghit",   m_ngch,m_ghit);
        status = m_tuple4->addIndexedItem ("thit",   m_ngch,m_thit);
 
        status = m_tuple4->addIndexedItem ("ptot_etof", m_ngch,m_ptot_etof);
        status = m_tuple4->addIndexedItem ("cntr_etof", m_ngch,m_cntr_etof);
        status = m_tuple4->addIndexedItem ("te_etof",   m_ngch,m_te_etof);
        status = m_tuple4->addIndexedItem ("tmu_etof",  m_ngch,m_tmu_etof);
        status = m_tuple4->addIndexedItem ("tpi_etof",  m_ngch,m_tpi_etof);
        status = m_tuple4->addIndexedItem ("tk_etof",   m_ngch,m_tk_etof);
        status = m_tuple4->addIndexedItem ("tp_etof",   m_ngch,m_tp_etof);
        status = m_tuple4->addIndexedItem ("ph_etof",   m_ngch,m_ph_etof);
        status = m_tuple4->addIndexedItem ("rhit_etof", m_ngch,m_rhit_etof);
        status = m_tuple4->addIndexedItem ("qual_etof", m_ngch,m_qual_etof);
        status = m_tuple4->addIndexedItem ("ec_toff_e", m_ngch,m_ec_toff_e);
        status = m_tuple4->addIndexedItem ("ec_toff_mu",m_ngch,m_ec_toff_mu);
        status = m_tuple4->addIndexedItem ("ec_toff_pi",m_ngch,m_ec_toff_pi);
        status = m_tuple4->addIndexedItem ("ec_toff_k", m_ngch,m_ec_toff_k);
        status = m_tuple4->addIndexedItem ("ec_toff_p", m_ngch,m_ec_toff_p);
        status = m_tuple4->addIndexedItem ("ec_tsig_e", m_ngch,m_ec_tsig_e);
        status = m_tuple4->addIndexedItem ("ec_tsig_mu",m_ngch,m_ec_tsig_mu);
        status = m_tuple4->addIndexedItem ("ec_tsig_pi",m_ngch,m_ec_tsig_pi);
        status = m_tuple4->addIndexedItem ("ec_tsig_k", m_ngch,m_ec_tsig_k);
        status = m_tuple4->addIndexedItem ("ec_tsig_p", m_ngch,m_ec_tsig_p);
        status = m_tuple4->addIndexedItem ("ec_tof",    m_ngch,m_ec_tof);
        status = m_tuple4->addIndexedItem ("ptot_btof1",m_ngch,m_ptot_btof1);
        status = m_tuple4->addIndexedItem ("cntr_btof1",m_ngch,m_cntr_btof1);
        status = m_tuple4->addIndexedItem ("te_btof1",  m_ngch,m_te_btof1);
        status = m_tuple4->addIndexedItem ("tmu_btof1", m_ngch,m_tmu_btof1);
        status = m_tuple4->addIndexedItem ("tpi_btof1", m_ngch,m_tpi_btof1);
        status = m_tuple4->addIndexedItem ("tk_btof1",  m_ngch,m_tk_btof1);
        status = m_tuple4->addIndexedItem ("tp_btof1",  m_ngch,m_tp_btof1);
        status = m_tuple4->addIndexedItem ("ph_btof1",  m_ngch,m_ph_btof1);
        status = m_tuple4->addIndexedItem ("zhit_btof1",m_ngch,m_zhit_btof1);
        status = m_tuple4->addIndexedItem ("qual_btof1",m_ngch,m_qual_btof1);
        status = m_tuple4->addIndexedItem ("b1_toff_e", m_ngch,m_b1_toff_e);
        status = m_tuple4->addIndexedItem ("b1_toff_mu",m_ngch,m_b1_toff_mu);
        status = m_tuple4->addIndexedItem ("b1_toff_pi",m_ngch,m_b1_toff_pi);
        status = m_tuple4->addIndexedItem ("b1_toff_k", m_ngch,m_b1_toff_k);
        status = m_tuple4->addIndexedItem ("b1_toff_p", m_ngch,m_b1_toff_p);
        status = m_tuple4->addIndexedItem ("b1_tsig_e", m_ngch,m_b1_tsig_e);
        status = m_tuple4->addIndexedItem ("b1_tsig_mu",m_ngch,m_b1_tsig_mu);
        status = m_tuple4->addIndexedItem ("b1_tsig_pi",m_ngch,m_b1_tsig_pi);
        status = m_tuple4->addIndexedItem ("b1_tsig_k", m_ngch,m_b1_tsig_k);
        status = m_tuple4->addIndexedItem ("b1_tsig_p", m_ngch,m_b1_tsig_p);
        status = m_tuple4->addIndexedItem ("b1_tof",    m_ngch,m_b1_tof);                                                    
 
        status = m_tuple4->addIndexedItem ("mdedx_pid", m_ngch,m_dedx_pid);
        status = m_tuple4->addIndexedItem ("mtof1_pid", m_ngch,m_tof1_pid);
        status = m_tuple4->addIndexedItem ("mtof2_pid", m_ngch,m_tof2_pid);
        status = m_tuple4->addIndexedItem ("mprob_pid", m_ngch,m_prob_pid);
        status = m_tuple4->addIndexedItem ("mptrk_pid", m_ngch,m_ptrk_pid);
        status = m_tuple4->addIndexedItem ("mcost_pid", m_ngch,m_cost_pid);
        status = m_tuple4->addItem ("mpnp",   m_pnp);
        status = m_tuple4->addItem ("mpnm",   m_pnm);
 
        status = m_tuple4->addIndexedItem ("numHits",      m_nggneu,m_numHits);    // Total number of hits
        status = m_tuple4->addIndexedItem ("secondmoment", m_nggneu,m_secondmoment);
        status = m_tuple4->addIndexedItem ("mx",           m_nggneu,m_x);       //  Shower coordinates and errors
        status = m_tuple4->addIndexedItem ("my",           m_nggneu,m_y);
        status = m_tuple4->addIndexedItem ("mz",           m_nggneu,m_z);
        status = m_tuple4->addIndexedItem ("cosemc",       m_nggneu,m_cosemc);   // Shower Counter angles and errors
        status = m_tuple4->addIndexedItem ("phiemc",       m_nggneu,m_phiemc);
        status = m_tuple4->addIndexedItem ("energy",       m_nggneu,m_energy);  // Total energy observed in Emc
        status = m_tuple4->addIndexedItem ("eseed",        m_nggneu,m_eSeed);
        status = m_tuple4->addIndexedItem ("me9",          m_nggneu,m_e3x3); 
        status = m_tuple4->addIndexedItem ("me25",         m_nggneu,m_e5x5); 
        status = m_tuple4->addIndexedItem ("mlat",         m_nggneu,m_lat);
        status = m_tuple4->addIndexedItem ("ma20",         m_nggneu,m_a20);
        status = m_tuple4->addIndexedItem ("ma42",         m_nggneu,m_a42);  
 
      }
      else    { 
    log << MSG::ERROR << "    Cannot book N-tuple:" << long(m_tuple4) << endmsg;
    return StatusCode::FAILURE;
      }
    }
 
  if(service("THistSvc", m_thsvc).isFailure()) {
    log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
    return StatusCode::FAILURE;
  }
 
  // "DQAHist" is fixed
  // "Rhopi" is control sample name, just as ntuple case.
 
  TH1F* mrho0 = new TH1F("mrho0", "mass for #rho^{0}->#pi^{+}#pi^{-}", 160, 0.0, 3.2);
  if(m_thsvc->regHist("/DQAHist/Rhopi/mrho0", mrho0).isFailure()) {
      log << MSG::ERROR << "Couldn't register mrho0" << endreq;
  }
 
  TH1F* mrhop = new TH1F("mrhop", "mass for #rho^{+}->#pi^{+}#pi^{0}", 160, 0.0,3.2);
  if(m_thsvc->regHist("/DQAHist/Rhopi/mrhop", mrhop).isFailure()) {
      log << MSG::ERROR << "Couldn't register mrhop" << endreq;
  }
 
  TH1F* mrhom = new TH1F("mrhom", "mass for #rho^{-}->#pi^{-}#pi^{0}", 160, 0.0, 3.2);
  if(m_thsvc->regHist("/DQAHist/Rhopi/mrhom", mrhom).isFailure()) {
      log << MSG::ERROR << "Couldn't register mrhom" << endreq;
  }
 
 
  TH1F* mpi0 = new TH1F("mpi0", "mass for #pi^{0}->#gamma #gamma", 50,0.08, 0.18);
  if(m_thsvc->regHist("/DQAHist/Rhopi/mpi0", mpi0).isFailure()) {
      log << MSG::ERROR << "Couldn't register mpi0" << endreq;
  }
 
  //
  //--------end of book--------
  //
 
  log << MSG::INFO << "successfully return from initialize()" <<endmsg;
  return StatusCode::SUCCESS;
 
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode DQARhopi::execute() {
  
//  std::cout << "execute()" << std::endl;
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in execute()" << endreq;
 
    setFilterPassed(false);
 
  SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
  int run=eventHeader->runNumber();
  int event=eventHeader->eventNumber();
  log << MSG::DEBUG <<"run, evtnum = "
      << run << " , "
      << event <<endreq;
 
  m_run = eventHeader->runNumber();  
  m_rec = eventHeader->eventNumber();
 
//  cout<<"event "<<event<<endl;
  Ncut0++;
  SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
    log << MSG::INFO << "get event tag OK" << endreq;
    log << MSG::DEBUG <<"ncharg, nneu, tottks = " 
      << evtRecEvent->totalCharged() << " , "
      << evtRecEvent->totalNeutral() << " , "
      << evtRecEvent->totalTracks() <<endreq;
 
  m_nch  = evtRecEvent->totalCharged();
  m_nneu = evtRecEvent->totalNeutral();
 
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
 
  //
  // check x0, y0, z0, r0
  // suggest cut: |z0|<5 && r0<1
  //
  Vint iGood, ipip, ipim, ipnp,ipnm;
  iGood.clear();
  ipip.clear();
  ipim.clear();
  ipnp.clear();
  ipnm.clear();
  Vp4 ppip, ppim , ppnp, ppnm;
  ppip.clear();
  ppim.clear();
  ppnp.clear();
  ppnm.clear();
 
  Hep3Vector xorigin(0,0,0);
 
 
  //if (m_reader.isRunNumberValid(runNo)) {
   IVertexDbSvc*  vtxsvc;
  Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
  if(vtxsvc->isVertexValid()){ 
  double* dbv = vtxsvc->PrimaryVertex();
  double*  vv = vtxsvc->SigmaPrimaryVertex();
    xorigin.setX(dbv[0]);
    xorigin.setY(dbv[1]);
    xorigin.setZ(dbv[2]);
  }
 
  int nCharge = 0;
  for(int i = 0; i < evtRecEvent->totalCharged(); i++){
    EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
    if(!(*itTrk)->isMdcTrackValid()) continue;
    if (!(*itTrk)->isMdcKalTrackValid()) continue;
    RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
    RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();    
 
    double pch =mdcTrk->p();
    double x0  =mdcTrk->x();
    double y0  =mdcTrk->y();
    double z0  =mdcTrk->z();
    double phi0=mdcTrk->helix(1);
    double xv=xorigin.x();
    double yv=xorigin.y();
    double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0));
    double m_vx0 = x0;
    double m_vy0 = y0;
    double m_vz0 = z0-xorigin.z();
    double m_vr0 = Rxy;
    double m_Vct=cos(mdcTrk->theta());
//    m_tuple1->write();
    if(fabs(m_vz0) >= m_vz0cut) continue;
    if(m_vr0 >= m_vr0cut) continue;
    if(fabs(m_Vct)>=m_cthcut) continue;
 
    iGood.push_back((*itTrk)->trackId());
    nCharge += mdcTrk->charge();
 
  }
  
  //
  // Finish Good Charged Track Selection
  //
  int nGood = iGood.size();
  log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
  if((nGood != 2)||(nCharge!=0)){
    return StatusCode::SUCCESS;
  }
  Ncut1++;
 
  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;
 
      if(fabs(thed) < fabs(dthe)) dthe = thed;
      if(fabs(phid) < fabs(dphi)) dphi = phid;
      if(angd < dang) dang = angd;
    }
    if(dang>=200) continue;
    double eraw = emcTrk->energy();
    double theta1 = emcTrk->theta();
    double time = emcTrk->time();
    dthe = dthe * 180 / (CLHEP::pi);
    dphi = dphi * 180 / (CLHEP::pi);
    dang = dang * 180 / (CLHEP::pi);
    double m_dthe = dthe;
    double m_dphi = dphi;
    double m_dang = dang;
    double m_eraw = eraw;
//    m_tuple2->write();
    double fc_theta = fabs(cos(theta1)); 
    if ( fc_theta < 0.80) {     // Barrel EMC
        if (eraw <= m_energyThreshold/2) continue;
    }
    else if  ( fc_theta >0.86 && fc_theta < 0.92 ) { // Endcap EMC
        if (eraw <= m_energyThreshold) continue;
    }
    else continue;
 
    if (time < 0 || time > 14) continue;
    if(dang < m_gammaAngCut) continue;
    //
    // good photon cut will be set here
    //
    iGam.push_back((*itTrk)->trackId());
  }
  
  //
  // Finish Good Photon Selection
  //
  int nGam = iGam.size();
 
  log << MSG::DEBUG << "num Good Photon " << nGam  << " , " <<evtRecEvent->totalNeutral()<<endreq;
  if(nGam<2){
    return StatusCode::SUCCESS;
  }
  Ncut2++;
 
 
  //
  // Assign 4-momentum to each photon
  // 
 
  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);
 
//    ptrk = ptrk.boost(-0.011,0,0);// boost to cms
 
    pGam.push_back(ptrk);
  }
 
  log << MSG::DEBUG << "liuf1 Good Photon " <<endreq;
 
 
  for(int i = 0; i < nGood; i++) {//for rhopi without PID
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
    if(!(*itTrk)->isMdcTrackValid()) continue;
    if (!(*itTrk)->isMdcKalTrackValid()) continue;
    RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack(); 
    RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
    RecMdcKalTrack::setPidType  (RecMdcKalTrack::pion);
 
    if(mdcTrk->charge() >0) {
      ipip.push_back(iGood[i]);
      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));
      ppip.push_back(ptrk);
    } else {
      ipim.push_back(iGood[i]);
      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));
      ppim.push_back(ptrk);
    }
  }// without PID
 
  int npip = ipip.size();
  int npim = ipim.size();
  log << MSG::DEBUG << "num of pion "<< ipip.size()<<","<<ipim.size()<<endreq;
  Ncut3++;
 
 
 
  //***********************************************//
  // the angle between the two charged tracks      //
  //***********************************************//
 
    int langcc=0;
    double dthec = 200.; 
    double dphic = 200.; 
    double dangc = 200.; 
  for(int i=0; i < npip; i++) {
    EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + ipip[i] ;
    RecMdcTrack* mdcTrk1 = (*itTrk)->mdcTrack();
    RecMdcKalTrack* mdcKalTrk1 = (*itTrk)->mdcKalTrack();
    Hep3Vector emcpos(ppip[i].px(), ppip[i].py(), ppip[i].pz());
   
    for(int j = 0; j < npim; j++) {
      EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + ipim[j];
      RecMdcTrack* mdcTrk2 = (*jtTrk)->mdcTrack();
      RecMdcKalTrack* mdcKalTrk2 = (*jtTrk)->mdcKalTrack();
 
      HepLorentzVector pippim=ppip[i]+ppim[j];
 
      Hep3Vector extpos(ppim[j].px(), ppim[j].py(), ppim[j].pz());
       
      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;
 
      m_dthec[langcc] = thed * 180 / (CLHEP::pi);
      m_dphic[langcc] = phid * 180 / (CLHEP::pi);
      m_dangc[langcc] = angd * 180 / (CLHEP::pi);
      m_mspippim[langcc] =pippim.m();
      langcc++;
    }
   }
      m_nangecc=langcc;
 
  //
  // Loop each gamma pair, check ppi0 and pTot 
  //
 
  double m_m2gg,m_momentpi0;
  HepLorentzVector pTot,p2g;
 
  log << MSG::DEBUG << "liuf2 Good Photon " <<langcc<<endreq;
    //******************************************************//
    //   asign the momentum of MDC and KALFIT               //
    //   the deposite energy of EMC  for charged tracks     //
    //******************************************************// 
    double m_momentpip,m_momentpim,extmot1,extmot2;
    double emcTg1=0.0;
    double emcTg2=0.0;
    double nlaypi1=0;
    double nhit1=0;
    double nlaypi2=0;
    double nhit2=0;
 
      EvtRecTrackIterator  itTrk11 = evtRecTrkCol->begin() + ipip[0];
      RecMdcTrack* mdcTrk11 = (*itTrk11)->mdcTrack();
      RecMdcKalTrack* mdcKalTrk11 = (*itTrk11)->mdcKalTrack();
      RecEmcShower* emcTrk11 = (*itTrk11)->emcShower();
      RecMucTrack *mucTrk11=(*itTrk11)->mucTrack();
      double phi01=mdcTrk11->helix(1);
 
      EvtRecTrackIterator  itTrk12 = evtRecTrkCol->begin() + ipim[0];
      RecMdcTrack* mdcTrk12 = (*itTrk12)->mdcTrack();
      RecMdcKalTrack* mdcKalTrk12 = (*itTrk12)->mdcKalTrack();
      RecEmcShower* emcTrk12 = (*itTrk12)->emcShower();
      RecMucTrack *mucTrk12=(*itTrk12)->mucTrack();
      double phi02=mdcTrk12->helix(1);
 
    m_vxpin = mdcTrk11->x();
    m_vypin = mdcTrk11->y();
    m_vzpin = mdcTrk11->z()-xorigin.z();
    m_vrpin = fabs((mdcTrk11->x()-xorigin.x())*cos(phi01)+(mdcTrk11->y()-xorigin.y())*sin(phi01));
    m_costhepin =cos(mdcTrk11->theta());
 
    m_momentpip=mdcTrk11->p();
    m_ppx  =mdcTrk11->px();
    m_ppy  =mdcTrk11->py();
    m_ppz  =mdcTrk11->pz();
 
    m_vxp = mdcKalTrk11->x();
    m_vyp = mdcKalTrk11->y();
    m_vzp = mdcKalTrk11->z()-xorigin.z();
    m_vrp = fabs((mdcKalTrk11->x()-xorigin.x())*cos(phi01)+(mdcKalTrk11->y()-xorigin.y())*sin(phi01));
    m_costhep =cos(mdcKalTrk11->theta());
 
//    extmot1=mdcKalTrk11->p(); // only has value when read data from dst file
    m_ppxkal   =mdcKalTrk11->px();
    m_ppykal   =mdcKalTrk11->py();
    m_ppzkal   =mdcKalTrk11->pz();
    extmot1 = sqrt(m_ppxkal*m_ppxkal + m_ppykal*m_ppykal + m_ppzkal*m_ppzkal);
 
    m_vxmin = mdcTrk12->x();  
    m_vymin = mdcTrk12->y();  
    m_vzmin = mdcTrk12->z();  
    m_vrmin = fabs((mdcTrk12->x()-xorigin.x())*cos(phi02)+(mdcTrk12->y()-xorigin.y())*sin(phi02));
    m_costhemin =cos(mdcTrk12->theta());
 
    m_momentpim=mdcTrk12->p();
    m_mpx  =mdcTrk12->px();
    m_mpy  =mdcTrk12->py();
    m_mpz  =mdcTrk12->pz();
 
    m_vxm = mdcKalTrk12->x();
    m_vym = mdcKalTrk12->y();
    m_vzm = mdcKalTrk12->z();  
    m_vrm = fabs((mdcKalTrk12->x()-xorigin.x())*cos(phi02)+(mdcKalTrk12->y()-xorigin.y())*sin(phi02));
    m_costhem =cos(mdcKalTrk12->theta());
 
//    extmot2    =mdcKalTrk12->p();
    m_mpxkal   =mdcKalTrk12->px();
    m_mpykal   =mdcKalTrk12->py();
    m_mpzkal   =mdcKalTrk12->pz();
    extmot2 = sqrt(m_mpxkal*m_mpxkal + m_mpykal*m_mpykal + m_mpzkal*m_mpzkal);
 
    if((*itTrk11)->isEmcShowerValid()){
    emcTg1 = emcTrk11->energy();
    }
    if((*itTrk12)->isEmcShowerValid()){
    emcTg2 = emcTrk12->energy();
    }
    if((*itTrk11)->isMucTrackValid()){
    nlaypi1=mucTrk11->numLayers();
    nhit1  =mucTrk11->numHits();
    }
    if((*itTrk12)->isMucTrackValid()){
    nlaypi2=mucTrk12->numLayers();
    nhit2  =mucTrk12->numHits();  
    }
 
    m_laypi1=nlaypi1;
    m_hit1  =nhit1;  
    m_laypi2=nlaypi2;
    m_hit2  =nhit2;  
 
  log << MSG::DEBUG << "liuf3 Good Photon " << ipim[0] <<endreq;
 
  RecMdcKalTrack *pipTrk = (*(evtRecTrkCol->begin()+ipip[0]))->mdcKalTrack();
  RecMdcKalTrack *pimTrk = (*(evtRecTrkCol->begin()+ipim[0]))->mdcKalTrack();
 
  log << MSG::DEBUG << "liuf4 Good Photon "  <<endreq;
 
  WTrackParameter wvpipTrk, wvpimTrk,wvkpTrk, wvkmTrk;
  wvpipTrk = WTrackParameter(mpi, pipTrk->getZHelix(), pipTrk->getZError());
  wvpimTrk = WTrackParameter(mpi, pimTrk->getZHelix(), pimTrk->getZError());
 
  wvkpTrk = WTrackParameter(mk, pipTrk->getZHelixK(), pipTrk->getZErrorK());//kaon
  wvkmTrk = WTrackParameter(mk, pimTrk->getZHelixK(), pimTrk->getZErrorK());//kaon
 
/* Default is pion, for other particles:
  wvppTrk = WTrackParameter(mp, pipTrk->getZHelixP(), pipTrk->getZErrorP());//proton
  wvmupTrk = WTrackParameter(mmu, pipTrk->getZHelixMu(), pipTrk->getZErrorMu());//muon
  wvepTrk = WTrackParameter(me, pipTrk->getZHelixE(), pipTrk->getZErrorE());//electron
  wvkpTrk = WTrackParameter(mk, pipTrk->getZHelixK(), pipTrk->getZErrorK());//kaon
*/
  //
  //    Test vertex fit
  //
 
  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);
  
  //****************************************************//
  //    Test vertex fit                                 //
  //****************************************************//
 
  VertexFit* vtxfit = VertexFit::instance();
 
  //****************************************************//
  //   if the charged particle is kaon                  //
  //****************************************************//
    double chi5=9999.0;
    double jkpi0=-0.5;
    double jkpkm=0.0;
    double jkpp0=0.0;
    double jkmp0=0.0;  
  vtxfit->init();
  vtxfit->AddTrack(0,  wvkpTrk);
  vtxfit->AddTrack(1,  wvkmTrk);
  vtxfit->AddVertex(0, vxpar,0, 1);
  if(vtxfit->Fit(0)) {
  vtxfit->Swim(0);  
  WTrackParameter wkp = vtxfit->wtrk(0);
  WTrackParameter wkm = vtxfit->wtrk(1);
 
  KinematicFit * kmfit = KinematicFit::instance();
 
  //
  //  Apply Kinematic 4C fit
  // 
 
    double chisq = 9999.;
    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->setDynamicerror(1);
    kmfit->AddTrack(0, wkp);
    kmfit->AddTrack(1, wkm);
    kmfit->AddTrack(2, 0.0, g1Trk);
    kmfit->AddTrack(3, 0.0, g2Trk);
    kmfit->AddFourMomentum(0, ecms);
    bool oksq = kmfit->Fit();
    
    if(oksq) {
      double chi2 = kmfit->chisq();
      if(chi2 < chi5) {
    HepLorentzVector kpi0 = kmfit->pfit(2) + kmfit->pfit(3);
        HepLorentzVector kpkm = kmfit->pfit(0) + kmfit->pfit(1);
        HepLorentzVector kpp0 = kmfit->pfit(0) + kpi0;
        HepLorentzVector kmp0 = kmfit->pfit(1) + kpi0;
    chi5 = kmfit->chisq();
    jkpi0 = kpi0.m();
        jkpkm= kpkm.m();
        jkpp0= kpp0.m();
        jkmp0= kmp0.m();
      }
    }
      }
    }
    } // vetex is true
 
 
  //****************************************************//
  //   if the charged particles are pions for real      //
  //****************************************************//
 
  vtxfit->init();
  vtxfit->AddTrack(0,  wvpipTrk);
  vtxfit->AddTrack(1,  wvpimTrk);
  vtxfit->AddVertex(0, vxpar,0, 1);
  if(!vtxfit->Fit(0)) return SUCCESS;
  vtxfit->Swim(0);
  
  WTrackParameter wpip = vtxfit->wtrk(0);
  WTrackParameter wpim = vtxfit->wtrk(1);
 
  log << MSG::DEBUG << "liuf5 Good Photon " <<endreq;
 
  KinematicFit * kmfit = KinematicFit::instance();
 
  //
  //  Apply Kinematic 4C fit
  // 
 
    double chisq = 9999.;
    int ig1 = -1;
    int ig2 = -1;
    HepLorentzVector gg1,gg2;    
    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->setDynamicerror(1);
    kmfit->AddTrack(0, wpip);
    kmfit->AddTrack(1, wpim);
    kmfit->AddTrack(2, 0.0, g1Trk);
    kmfit->AddTrack(3, 0.0, g2Trk);
    kmfit->AddFourMomentum(0, ecms);
    bool oksq = kmfit->Fit();
    if(oksq) {
      double chi2 = kmfit->chisq();
      if(chi2 < chisq) {
        chisq = chi2;
        ig1 = iGam[i];
        ig2 = iGam[j];
            gg1 = pGam[i];
            gg2 = pGam[j];
      }
    }
      }
    }
    
      p2g = gg1 + gg2;
      m_pmax=gg1.e()>gg2.e()?gg1.e():gg2.e();
      m_m2gg = p2g.m();
      m_cosang=(gg1.e()-gg2.e())/p2g.rho();
      m_momentpi0=sqrt(p2g.px()*p2g.px()+p2g.py()*p2g.py()+p2g.pz()*p2g.pz());  
        log << MSG::DEBUG << " chisq for 4c " << chisq <<endreq; 
    if(chisq > 200) {
       return StatusCode::SUCCESS; 
      }
 
// select charge track and nneu track
  Vint jGood;
  jGood.clear();
  jGood.push_back(ipip[0]);  
  jGood.push_back(ipim[0]);
  m_ngch = jGood.size();
 
  Vint jGgam;
  jGgam.clear();
  jGgam.push_back(ig1);
  jGgam.push_back(ig2);
  m_nggneu=jGgam.size();
 
     HepLorentzVector ppip1=ppip[0];
     HepLorentzVector ppim1=ppim[0];
 
     HepLorentzVector Ppipboost = ppip1.boost(u_cms);
     HepLorentzVector Ppimboost = ppim1.boost(u_cms);
     Hep3Vector p3pi1 = Ppipboost.vect();          //pi1
     Hep3Vector p3pi2 = Ppimboost.vect();         //pi2   
     m_anglepm   = (p3pi1.angle(p3pi2))* 180 / (CLHEP::pi);     
 
//*******************************************************//
 
      RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+ig1))->emcShower();
      RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+ig2))->emcShower();
      kmfit->init();
//      kmfit->setDynamicerror(1);
      kmfit->AddTrack(0, wpip);
      kmfit->AddTrack(1, wpim);
      kmfit->AddTrack(2, 0.0, g1Trk);
      kmfit->AddTrack(3, 0.0, g2Trk);
      kmfit->AddFourMomentum(0, ecms);
      bool oksq = kmfit->Fit();
      if(!oksq) return SUCCESS;
 
    HepLorentzVector ppi0 = kmfit->pfit(2) + kmfit->pfit(3);
        HepLorentzVector prho0 = kmfit->pfit(0) + kmfit->pfit(1);
        HepLorentzVector prhop = kmfit->pfit(0) + ppi0;
        HepLorentzVector prhom = kmfit->pfit(1) + ppi0;
        HepLorentzVector pgam2pi1 = prho0 + kmfit->pfit(2);
        HepLorentzVector pgam2pi2 = prho0 + kmfit->pfit(3);
        HepLorentzVector pgam11 =kmfit->pfit(2);
        HepLorentzVector pgam12 =kmfit->pfit(3);
 
/*
        HepLorentzVector ppi0_a=ppi0;
        HepLorentzVector pgam11_a =pgam11;
        HepLorentzVector pgam12_a =pgam12;
 
        Hep3Vector pi0_cms = -ppi0_a.boostVector();
        HepLorentzVector pgam1 = pgam11_a.boost(pi0_cms);
        HepLorentzVector pgam2 = pgam12_a.boost(pi0_cms);
*/
    m_chi1 = kmfit->chisq();
    m_mpi0 = ppi0.m();
        m_prho0= prho0.m();
        m_prhop= prhop.m();
        m_prhom= prhom.m();
        m_good = nGood;
        m_gam  = nGam;
 
        m_pip  = npip;
        m_pim  = npim;
        m_2gam = m_m2gg;
        m_outpi0=m_momentpi0;
        m_outpip=m_momentpip;
        m_outpim=m_momentpim;
        m_enpip =emcTg1;
        m_dcpip =extmot1;
        m_enpim =emcTg2;
        m_dcpim =extmot2;
        m_pipf=kmfit->pfit(0).rho();
        m_pimf=kmfit->pfit(1).rho();
        m_pi0f=ppi0.rho();
 
        m_chi5=chi5;
        m_kpi0=jkpi0;
        m_kpkm=jkpkm;
        m_kpp0=jkpp0;
        m_kmp0=jkmp0;
        m_pgam2pi1=pgam2pi1.m();
        m_pgam2pi2=pgam2pi2.m();
        cosva1=kmfit->pfit(2).rho();
        cosva2=kmfit->pfit(3).rho();
//        m_pe1   =pe1;
//        m_pe2   =pe2;
 
//
//          fill information of dedx and tof
//
 
  //
  // check dedx infomation
  //
 
    for(int ii = 0; ii < m_ngch; ii++) {
// dedx
      m_ptrk[ii] = 9999.0;
      m_chie[ii] = 9999.0;
      m_chimu[ii] = 9999.0;
      m_chipi[ii] = 9999.0;
      m_chik[ii] = 9999.0;
      m_chip[ii] = 9999.0;
      m_ghit[ii] = 9999.0;
      m_thit[ii] = 9999.0;
      m_probPH[ii] = 9999.0;
      m_normPH[ii] = 9999.0;
 
//endtof
          m_cntr_etof[ii] = 9999.0;
          m_ptot_etof[ii] = 9999.0;
          m_ph_etof[ii]   = 9999.0;
          m_rhit_etof[ii] = 9999.0;
          m_qual_etof[ii] = 9999.0;
          m_te_etof[ii]   = 9999.0;
          m_tmu_etof[ii]  = 9999.0;
          m_tpi_etof[ii]  = 9999.0;
          m_tk_etof[ii]   = 9999.0;
          m_tp_etof[ii]   = 9999.0;
          m_ec_tof[ii]    =  9999.0;
          m_ec_toff_e[ii]   = 9999.0;
          m_ec_toff_mu[ii]  = 9999.0;
          m_ec_toff_pi[ii]  = 9999.0;
          m_ec_toff_k[ii]   = 9999.0;
          m_ec_toff_p[ii]   = 9999.0;
          m_ec_tsig_e[ii]   = 9999.0;
          m_ec_tsig_mu[ii]  = 9999.0;
          m_ec_tsig_pi[ii]  = 9999.0;
          m_ec_tsig_k[ii]   = 9999.0;
          m_ec_tsig_p[ii]   = 9999.0;
 
// barrel tof
          m_cntr_btof1[ii]  = 9999.0;
          m_ptot_btof1[ii]  = 9999.0;
          m_ph_btof1[ii]    = 9999.0;
          m_zhit_btof1[ii]  = 9999.0;
          m_qual_btof1[ii]  = 9999.0;
          m_te_btof1[ii]    = 9999.0;
          m_tmu_btof1[ii]   = 9999.0;
          m_tpi_btof1[ii]   = 9999.0;
          m_tk_btof1[ii]    = 9999.0;
          m_tp_btof1[ii]    = 9999.0;
          m_b1_tof[ii]      = 9999.0;
          m_b1_toff_e[ii]   = 9999.0;
          m_b1_toff_mu[ii]  = 9999.0;
          m_b1_toff_pi[ii]  = 9999.0;
          m_b1_toff_k[ii]   = 9999.0;
          m_b1_toff_p[ii]   = 9999.0;
          m_b1_tsig_e[ii]   = 9999.0;
          m_b1_tsig_mu[ii]  = 9999.0;
          m_b1_tsig_pi[ii]  = 9999.0;
          m_b1_tsig_k[ii]   = 9999.0;
          m_b1_tsig_p[ii]   = 9999.0;
//pid
          m_dedx_pid[ii] = 9999.0;
          m_tof1_pid[ii] = 9999.0;
          m_tof2_pid[ii] = 9999.0;
          m_prob_pid[ii] = 9999.0; 
          m_ptrk_pid[ii] = 9999.0;  
          m_cost_pid[ii] = 9999.0;
       }
 
 
      int  indx0=0;
      for(int i = 0; i < m_ngch; i++) {
      EvtRecTrackIterator  itTrk = evtRecTrkCol->begin() + jGood[i];
      if(!(*itTrk)->isMdcTrackValid()) continue;
      if(!(*itTrk)->isMdcDedxValid())continue;
      RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
      RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
      m_ptrk[indx0] = mdcTrk->p();
 
      m_chie[indx0] = dedxTrk->chiE();
      m_chimu[indx0] = dedxTrk->chiMu();
      m_chipi[indx0] = dedxTrk->chiPi();
      m_chik[indx0] = dedxTrk->chiK();
      m_chip[indx0] = dedxTrk->chiP();
      m_ghit[indx0] = dedxTrk->numGoodHits();
      m_thit[indx0] = dedxTrk->numTotalHits();
      m_probPH[indx0] = dedxTrk->probPH();
      m_normPH[indx0] = dedxTrk->normPH();
      indx0++;
    }
 
 
  //
  // check TOF infomation
  //
 
 
     int indx1=0;
    for(int i = 0; i < m_ngch; i++) {
      EvtRecTrackIterator  itTrk = evtRecTrkCol->begin() + jGood[i];
      if(!(*itTrk)->isMdcTrackValid()) continue;
      if(!(*itTrk)->isTofTrackValid()) continue;
 
      RecMdcTrack * mdcTrk = (*itTrk)->mdcTrack();
      SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
 
      double ptrk = mdcTrk->p();
      SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
      for(;iter_tof != tofTrkCol.end(); iter_tof++ ) { 
        TofHitStatus *status = new TofHitStatus; 
        status->setStatus((*iter_tof)->status());
        if(!(status->is_barrel())){//endcap
          if( !(status->is_counter()) ) continue; // ? 
          if( status->layer()!=1 ) continue;//layer1
          double path=(*iter_tof)->path(); // ? 
          double tof  = (*iter_tof)->tof();
          double ph   = (*iter_tof)->ph();
          double rhit = (*iter_tof)->zrhit();
          double qual = 0.0 + (*iter_tof)->quality();
          double cntr = 0.0 + (*iter_tof)->tofID();
          double texp[5];
          double tsig[5];
          for(int j = 0; j < 5; j++) {//0 e, 1 mu, 2 pi, 3 K, 4 p
           texp[j] = (*iter_tof)->texp(j);
//           tsig[j] = (*iter_tof)->sigma(j);
//           toffset[j] = (*iter_tof)->offset(j);
          }
          m_cntr_etof[indx1]  = cntr;
          m_ptot_etof[indx1] = ptrk;
          m_ph_etof[indx1]   = ph;
          m_rhit_etof[indx1]  = rhit;
          m_qual_etof[indx1]  = qual;
          m_te_etof[indx1]    = tof - texp[0];
          m_tmu_etof[indx1]   = tof - texp[1];
          m_tpi_etof[indx1]   = tof - texp[2];
          m_tk_etof[indx1]    = tof - texp[3];
          m_tp_etof[indx1]    = tof - texp[4];
 
          m_ec_tof[indx1]      =  tof;
 
          m_ec_toff_e[indx1]   = (*iter_tof)->toffset(0);
          m_ec_toff_mu[indx1]  = (*iter_tof)->toffset(1);
          m_ec_toff_pi[indx1]  = (*iter_tof)->toffset(2);
          m_ec_toff_k[indx1]   = (*iter_tof)->toffset(3);
          m_ec_toff_p[indx1]   = (*iter_tof)->toffset(4);
 
          m_ec_tsig_e[indx1]   = (*iter_tof)->sigma(0);
          m_ec_tsig_mu[indx1]  = (*iter_tof)->sigma(1);
          m_ec_tsig_pi[indx1]  = (*iter_tof)->sigma(2);
          m_ec_tsig_k[indx1]   = (*iter_tof)->sigma(3);
          m_ec_tsig_p[indx1]   = (*iter_tof)->sigma(4);
 
        }
        else {//barrel
          if( !(status->is_cluster()) ) continue; // ? 
            double path=(*iter_tof)->path(); // ? 
            double tof  = (*iter_tof)->tof();
            double ph   = (*iter_tof)->ph();
            double rhit = (*iter_tof)->zrhit();
            double qual = 0.0 + (*iter_tof)->quality();
            double cntr = 0.0 + (*iter_tof)->tofID();
          double texp[5];
          for(int j = 0; j < 5; j++) {
            texp[j] = (*iter_tof)->texp(j);
            }
            m_cntr_btof1[indx1]  = cntr;
            m_ptot_btof1[indx1] = ptrk;
            m_ph_btof1[indx1]   = ph;
            m_zhit_btof1[indx1]  = rhit;
            m_qual_btof1[indx1]  = qual;
            m_te_btof1[indx1]    = tof - texp[0];
            m_tmu_btof1[indx1]   = tof - texp[1];
            m_tpi_btof1[indx1]   = tof - texp[2];
            m_tk_btof1[indx1]    = tof - texp[3];
            m_tp_btof1[indx1]    = tof - texp[4];
 
          m_b1_tof[indx1]      =  tof;
 
          m_b1_toff_e[indx1]   = (*iter_tof)->toffset(0);
          m_b1_toff_mu[indx1]  = (*iter_tof)->toffset(1);
          m_b1_toff_pi[indx1]  = (*iter_tof)->toffset(2);
          m_b1_toff_k[indx1]   = (*iter_tof)->toffset(3);
          m_b1_toff_p[indx1]   = (*iter_tof)->toffset(4);
 
          m_b1_tsig_e[indx1]   = (*iter_tof)->sigma(0);
          m_b1_tsig_mu[indx1]  = (*iter_tof)->sigma(1);
          m_b1_tsig_pi[indx1]  = (*iter_tof)->sigma(2);
          m_b1_tsig_k[indx1]   = (*iter_tof)->sigma(3);
          m_b1_tsig_p[indx1]   = (*iter_tof)->sigma(4);
 
          }
        delete status; 
      } 
      indx1++;
    } // loop all charged track
 
  //
  // Assign 4-momentum to each charged track
  //
    int indx2=0;
  ParticleID *pid = ParticleID::instance();
  for(int i = 0; i < m_ngch; i++) {
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i]; 
    //    if(pid) delete pid;
    pid->init();
    pid->setMethod(pid->methodProbability());
//    pid->setMethod(pid->methodLikelihood());  //for Likelihood Method  
 
    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->identify(pid->onlyPion());
    //  pid->identify(pid->onlyKaon());
    pid->calculate();
    if(!(pid->IsPidInfoValid())) continue;
    RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
 
    m_dedx_pid[indx2] = pid->chiDedx(2);
    m_tof1_pid[indx2] = pid->chiTof1(2);
    m_tof2_pid[indx2] = pid->chiTof2(2);
    m_prob_pid[indx2] = pid->probPion();
 
//    if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
//    if(pid->probPion() < 0.001) continue;
//    if(pid->pdf(2)<pid->pdf(3)) continue; //  for Likelihood Method(0=electron 1=muon 2=pion 3=kaon 4=proton) 
 
    RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();//After ParticleID, use RecMdcKalTrack substitute RecMdcTrack
    RecMdcKalTrack::setPidType  (RecMdcKalTrack::pion);//PID can set to electron, muon, pion, kaon and proton;The default setting is pion
 
//    m_ptrk_pid[indx2] = mdcKalTrk->p();
    m_cost_pid[indx2] = cos(mdcTrk->theta());
 
    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));
    //      ptrk = ptrk.boost(-0.011,0,0);//boost to cms
    m_ptrk_pid[indx2] = p3;
 
    if(mdcTrk->charge() >0 && pid->probPion() > pid->probKaon()) {
      ipnp.push_back(jGood[i]);
      ppnp.push_back(ptrk);
    }     //plus charge with with PID
      if(mdcTrk->charge() <0 && pid->probPion() > pid->probKaon()) {
      ipnm.push_back(jGood[i]);
      ppnm.push_back(ptrk);
    }     //minus charge with with PID
  }
  int npnp = ipnp.size();
  int npnm = ipnm.size();
 
  m_pnp  = npnp;
  m_pnm  = npnm;
 
  int iphoton = 0;
  for (int i=0; i<m_nggneu; i++)
  {
    EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + jGgam[i];
    if (!(*itTrk)->isEmcShowerValid()) continue; 
    RecEmcShower *emcTrk = (*itTrk)->emcShower();
    m_numHits[iphoton] = emcTrk->numHits();
    m_secondmoment[iphoton] = emcTrk->secondMoment();
    m_x[iphoton] = emcTrk->x();
    m_y[iphoton]= emcTrk->y();
    m_z[iphoton]= emcTrk->z();
    m_cosemc[iphoton]   = cos(emcTrk->theta());
    m_phiemc[iphoton]   = emcTrk->phi();      
    m_energy[iphoton]   = emcTrk->energy();
    m_eSeed[iphoton]    = emcTrk->eSeed(); 
    m_e3x3[iphoton]     = emcTrk->e3x3();  
    m_e5x5[iphoton]     = emcTrk->e5x5();  
    m_lat[iphoton]      = emcTrk->latMoment();
    m_a20[iphoton]      = emcTrk->a20Moment();
    m_a42[iphoton]      = emcTrk->a42Moment();
    iphoton++;
  }
    m_tuple4->write();
         Ncut4++;
 
      if(kmfit->chisq()>=chi5) return SUCCESS;
      if(pgam2pi2.m()<=1.0) return SUCCESS;
      if(pgam2pi1.m()<=1.0) return SUCCESS;
      if(nGam!=2)  return SUCCESS;
      if(emcTg1/extmot1>=0.8)  return SUCCESS;
      if(emcTg2/extmot2>=0.8)  return SUCCESS;
         Ncut5++;
 
    // DQA
    TH1 *h(0);
    if (m_thsvc->getHist("/DQAHist/Rhopi/mpi0", h).isSuccess()) {
        h->Fill(ppi0.m());
    } else {
        log << MSG::ERROR << "Couldn't retrieve mpi0" << endreq;
    }
 
      if(fabs(ppi0.m()-0.135)>=0.015) return SUCCESS;      
         Ncut6++;
 
    if (m_thsvc->getHist("/DQAHist/Rhopi/mrho0", h).isSuccess()) {
        h->Fill(prho0.m());
    } else {
        log << MSG::ERROR << "Couldn't retrieve mrho0" << endreq;
    }
 
 
    if (m_thsvc->getHist("/DQAHist/Rhopi/mrhop", h).isSuccess()) {
        h->Fill(prhop.m());
    } else {
        log << MSG::ERROR << "Couldn't retrieve mrhop" << endreq;
    }
 
    if (m_thsvc->getHist("/DQAHist/Rhopi/mrhom", h).isSuccess()) {
        h->Fill(prhom.m());
    } else {
        log << MSG::ERROR << "Couldn't retrieve mrhom" << endreq;
    }
 
 
 // Pid: 1 - electron, 2 - muon, 3 - pion, 4 - kaon, 5 - proton
 
  (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(3);
  (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(3);
 
  // Quality: defined by whether dE/dx or TOF is used to identify particle
  // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
  // 1 - only dE/dx (can be used for TOF calibration)
  // 2 - only TOF (can be used for dE/dx calibration)
  // 3 - Both dE/dx and TOF
 
  (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(0);
  (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(0);
 
         setFilterPassed(true);
 
  return StatusCode::SUCCESS;
}
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode DQARhopi::finalize() {
  cout<<"total number:         "<<Ncut0<<endl;
  cout<<"nGood==2, nCharge==0: "<<Ncut1<<endl;
  cout<<"nGam>=2:              "<<Ncut2<<endl;
  cout<<"Pass Pid:             "<<Ncut3<<endl;
  cout<<"Pass 4C:              "<<Ncut4<<endl;
  cout<<"final cut without pi0:"<<Ncut5<<endl;
  cout<<"final cut with pi0:   "<<Ncut6<<endl;
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in finalize()" << endmsg;
  return StatusCode::SUCCESS;
}