DQARhopi Class Reference

#include <DQARhopi.h>

Inheritance diagram for DQARhopi:

Public Member Functions
	DQARhopi (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()

Detailed Description

Definition at line 8 of file DQARhopi.h.

Constructor & Destructor Documentation

DQARhopi()

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

Definition at line 54 of file DQARhopi.cxx.

    : 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 );
}

Member Function Documentation

execute()

StatusCode DQARhopi::execute

(

)

Definition at line 299 of file DQARhopi.cxx.

                             {
 
  //  std::cout << "execute()" << std::endl;
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  setFilterPassed( false );
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  int                              run   = eventHeader->runNumber();
  int                              event = eventHeader->eventNumber();
  log << MSG::DEBUG << "run, evtnum = " << run << " , " << event << endmsg;
 
  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" << endmsg;
  log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
      << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
 
  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 ).ignore();
  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 << endmsg;
  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()
      << endmsg;
  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 " << endmsg;
 
  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() << endmsg;
  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 << endmsg;
  //******************************************************//
  //   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] << endmsg;
 
  RecMdcKalTrack* pipTrk = ( *( evtRecTrkCol->begin() + ipip[0] ) )->mdcKalTrack();
  RecMdcKalTrack* pimTrk = ( *( evtRecTrkCol->begin() + ipim[0] ) )->mdcKalTrack();
 
  log << MSG::DEBUG << "liuf4 Good Photon " << endmsg;
 
  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 StatusCode::SUCCESS;
  vtxfit->Swim( 0 );
 
  WTrackParameter wpip = vtxfit->wtrk( 0 );
  WTrackParameter wpim = vtxfit->wtrk( 1 );
 
  log << MSG::DEBUG << "liuf5 Good Photon " << endmsg;
 
  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 << endmsg;
  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 StatusCode::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().ignore();
  Ncut4++;
 
  if ( kmfit->chisq() >= chi5 ) return StatusCode::SUCCESS;
  if ( pgam2pi2.m() <= 1.0 ) return StatusCode::SUCCESS;
  if ( pgam2pi1.m() <= 1.0 ) return StatusCode::SUCCESS;
  if ( nGam != 2 ) return StatusCode::SUCCESS;
  if ( emcTg1 / extmot1 >= 0.8 ) return StatusCode::SUCCESS;
  if ( emcTg2 / extmot2 >= 0.8 ) return StatusCode::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" << endmsg; }
 
  if ( fabs( ppi0.m() - 0.135 ) >= 0.015 ) return StatusCode::SUCCESS;
  Ncut6++;
 
  if ( m_thsvc->getHist( "/DQAHist/Rhopi/mrho0", h ).isSuccess() ) { h->Fill( prho0.m() ); }
  else { log << MSG::ERROR << "Couldn't retrieve mrho0" << endmsg; }
 
  if ( m_thsvc->getHist( "/DQAHist/Rhopi/mrhop", h ).isSuccess() ) { h->Fill( prhop.m() ); }
  else { log << MSG::ERROR << "Couldn't retrieve mrhop" << endmsg; }
 
  if ( m_thsvc->getHist( "/DQAHist/Rhopi/mrhom", h ).isSuccess() ) { h->Fill( prhom.m() ); }
  else { log << MSG::ERROR << "Couldn't retrieve mrhom" << endmsg; }
 
  // 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;
}

finalize()

StatusCode DQARhopi::finalize

(

)

Definition at line 1269 of file DQARhopi.cxx.

                              {
  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;
}

initialize()

StatusCode DQARhopi::initialize

(

)

Definition at line 70 of file DQARhopi.cxx.

                                {
  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" << endmsg;
    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" << endmsg; }
 
  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" << endmsg; }
 
  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" << endmsg; }
 
  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" << endmsg; }
 
  //
  //--------end of book--------
  //
 
  log << MSG::INFO << "successfully return from initialize()" << endmsg;
  return StatusCode::SUCCESS;
}

---

The documentation for this class was generated from the following files:

• DQARhopi.h
• DQARhopi.cxx