EventPreSelectAlg Class Reference

#include <EventPreSelectAlg.h>

Inheritance diagram for EventPreSelectAlg:

Public Member Functions
StatusCode	initialize () override
StatusCode	execute () override
StatusCode	finalize () override

Detailed Description

Definition at line 11 of file EventPreSelectAlg.h.

Member Function Documentation

execute()

StatusCode EventPreSelectAlg::execute ( )

override

Definition at line 209 of file EventPreSelectAlg.cxx.

                                      {
  info() << "In execute" << endmsg;
 
  bool isBarrelBhabha   = false;
  bool isEndcapBhabha   = false;
  bool isBarrelDimu     = false;
  bool isEndcapDimu     = false;
  bool isHadron         = false;
  bool isBarrelDiphoton = false;
  bool isEndcapDiphoton = false;
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  if ( !eventHeader )
  {
    error() << "Could not retrieve event header" << endmsg;
    return StatusCode::FAILURE;
  }
 
  int run   = eventHeader->runNumber();
  int event = eventHeader->eventNumber();
  m_events++;
 
  SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
  if ( !evtRecEvent )
  {
    error() << "Could not retrieve EvtRecEvent" << endmsg;
    return StatusCode::FAILURE;
  }
 
  debug() << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
          << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
 
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
  if ( !evtRecTrkCol )
  {
    error() << "Could not retrieve EvtRecTrackCol" << endmsg;
    return StatusCode::FAILURE;
  }
 
  if ( evtRecEvent->totalTracks() != evtRecTrkCol->size() ) return StatusCode::SUCCESS;
 
  // -------- Good Charged Track Selection
  std::vector<EvtRecTrack*> goodTracks;
 
  int nCharge = 0;
  for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
  {
    EvtRecTrack* recTrk = ( *evtRecTrkCol )[i];
    if ( !recTrk->isMdcTrackValid() ) continue;
 
    RecMdcTrack* mdcTrk = recTrk->mdcTrack();
    double       vx0    = mdcTrk->x();
    double       vy0    = mdcTrk->y();
    double       vz0    = mdcTrk->z();
    double       vr0    = mdcTrk->r();
    double       theta0 = mdcTrk->theta();
    double       p0     = mdcTrk->p();
    double       pt0    = mdcTrk->pxy();
 
    if ( m_output )
    {
      m_vx0    = vx0;
      m_vy0    = vy0;
      m_vz0    = vz0;
      m_vr0    = vr0;
      m_theta0 = theta0;
      m_p0     = p0;
      m_pt0    = pt0;
      m_tuple1->write().ignore();
    }
 
    if ( fabs( vz0 ) >= m_vz0cut ) continue;
    if ( vr0 >= m_vr0cut ) continue;
    if ( pt0 >= m_pt0HighCut ) continue;
    if ( pt0 <= m_pt0LowCut ) continue;
 
    goodTracks.push_back( recTrk );
    nCharge += mdcTrk->charge();
  }
  int nGood = goodTracks.size();
 
  // -------- Good Photon Selection
  std::vector<EvtRecTrack*> goodGammas;
  for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
  {
    EvtRecTrack* recTrk = ( *evtRecTrkCol )[i];
    if ( !recTrk->isEmcShowerValid() ) continue;
 
    RecEmcShower* emcTrk = recTrk->emcShower();
    double        eraw   = emcTrk->energy();
    if ( m_output )
    {
      m_eraw = eraw;
      m_tuple2->write().ignore();
    }
 
    if ( eraw < m_energyThreshold ) continue;
    goodGammas.push_back( recTrk );
  }
  int nGam = goodGammas.size();
 
  // -------- Assign 4-momentum to each charged track
  Vp4 ppip{ 0 }, ppim{ 0 };
 
  double echarge = 0.; // total energy of charged track
  double ptot    = 0.; // total momentum of charged track
  double etot    = 0.; // total energy in MDC and EMC
  double eemc    = 0.; // total energy in EMC
  double pp      = 0.;
  double pm      = 0.;
 
  for ( const auto& goodTrk : goodTracks )
  {
    if ( goodTrk->isMdcTrackValid() )
    {
      // 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 If
      // use this algorithm in reconstruction job, mdcKalTrk->p()=0!!! I don't
      // know why.
 
      RecMdcTrack* mdcTrk = goodTrk->mdcTrack();
 
      ptot += mdcTrk->p();
 
      HepLorentzVector ptrk;
      ptrk.setPx( mdcTrk->px() );
      ptrk.setPy( mdcTrk->py() );
      ptrk.setPz( mdcTrk->pz() );
      double p3 = ptrk.mag();
      ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
      ptrk = ptrk.boost( -0.011, 0, 0 ); // boost to cms
 
      echarge += ptrk.e();
      etot += ptrk.e();
 
      if ( mdcTrk->charge() > 0 )
      {
        ppip.push_back( ptrk );
        pp = ptrk.rho();
      }
      else
      {
        ppim.push_back( ptrk );
        pm = ptrk.rho();
      }
    }
 
    if ( goodTrk->isEmcShowerValid() )
    {
      RecEmcShower* emcTrk = goodTrk->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
 
      eemc += ptrk.e();
      etot += ptrk.e();
    }
  }
 
  // -------- Assign 4-momentum to each photon
  Vp4 pGam{ 0 };
 
  double eneu = 0; // total energy of neutral track
  for ( const auto& goodGamTrk : goodGammas )
  {
    RecEmcShower*    emcTrk = goodGamTrk->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 );
    eneu += ptrk.e();
    etot += ptrk.e();
    eemc += ptrk.e();
  }
 
  double esum = echarge + eneu;
 
  // -------- Use EMC shower information only
  double maxE    = 0.;
  double secE    = 0.;
  double maxThe  = 999.;
  double maxPhi  = 999.;
  double secThe  = 999.;
  double secPhi  = 999.;
  int    npart   = 999.;
  double dphi    = 999.;
  double dthe    = 999.;
  int    mdcHit1 = 0.;
  int    mdcHit2 = 0.;
 
  SmartDataPtr<RecEmcShowerCol> aShowerCol( eventSvc(), "/Event/Recon/RecEmcShowerCol" );
  if ( !aShowerCol )
  {
    warning() << "Could not find RecEmcShowerCol" << endmsg;
    return StatusCode::SUCCESS;
  }
 
  int                       ishower = 0;
  RecEmcShowerCol::iterator iShowerCol;
  for ( iShowerCol = aShowerCol->begin(); iShowerCol != aShowerCol->end(); iShowerCol++ )
  {
 
    if ( ishower == 0 )
    {
      maxE   = ( *iShowerCol )->e5x5();
      maxThe = ( *iShowerCol )->theta();
      maxPhi = ( *iShowerCol )->phi();
      npart  = ( *iShowerCol )->module();
    }
    else if ( ishower == 1 )
    {
      secE   = ( *iShowerCol )->e5x5();
      secThe = ( *iShowerCol )->theta();
      secPhi = ( *iShowerCol )->phi();
    }
    else if ( ishower == 2 ) { break; }
 
    ishower++;
  }
 
  if ( aShowerCol->size() >= 2 )
  {
    dphi = ( fabs( maxPhi - secPhi ) - CLHEP::pi ) * 180. / CLHEP::pi;
    dthe = ( fabs( maxThe + secThe ) - CLHEP::pi ) * 180. / CLHEP::pi;
 
    double phi1 = maxPhi < 0 ? maxPhi + CLHEP::twopi : maxPhi;
    double phi2 = secPhi < 0 ? secPhi + CLHEP::twopi : secPhi;
 
    // Define sector (phi11,phi12) and (phi21,phi22)
    double phi11 = min( phi1, phi2 );
    double phi22 = max( phi1, phi2 );
    double phi12 = ( phi11 + phi22 - CLHEP::pi ) * 0.5;
    double phi21 = ( phi11 + phi22 + CLHEP::pi ) * 0.5;
    if ( phi12 < 0. ) phi12 += CLHEP::twopi;
    if ( phi21 > CLHEP::twopi ) phi21 -= CLHEP::twopi;
 
    SmartDataPtr<MdcDigiCol> mdcDigiCol( evtSvc(), "/Event/Digi/MdcDigiCol" );
    if ( !mdcDigiCol )
    {
      error() << "Could not find MdcDigiCol!" << endmsg;
      return StatusCode::FAILURE;
    }
    int hitnum = mdcDigiCol->size();
    for ( int i = 0; i < hitnum; i++ )
    {
      MdcDigi* digi   = dynamic_cast<MdcDigi*>( mdcDigiCol->containedObject( i ) );
      double   time   = digi->getTimeChannel();
      double   charge = digi->getChargeChannel();
 
      if ( time == 0x7FFFFFFF || charge == 0x7FFFFFFF ) continue;
 
      Identifier   id( digi->identify() );
      unsigned int iphi   = MdcID::wire( id );
      unsigned int ilayer = MdcID::layer( id );
 
      if ( ilayer >= 43 ) error() << "MDC(" << ilayer << "," << iphi << ")" << endmsg;
 
      double phi = CLHEP::twopi * iphi / idmax[ilayer];
 
      if ( WhetherSector( phi, phi11, phi12 ) ) mdcHit1++;
      else if ( WhetherSector( phi, phi21, phi22 ) ) mdcHit2++;
    }
  }
 
  // If it's bhabha, return.
  // if(nGood==2 && pp+pm>3.4) return StatusCode::SUCCESS;
 
  // -------- Select each event
  // Select bhabha
  if ( eemc > m_bhabhaEmcECut && maxE > m_bhabhaMaxECut && secE > m_bhabhaSecECut &&
       abs( dthe ) < m_bhabhaDTheCut &&
       ( ( dphi > m_bhabhaDPhiCut1 && dphi < m_bhabhaDPhiCut2 ) ||
         ( dphi > m_bhabhaDPhiCut3 && dphi < m_bhabhaDPhiCut4 ) ) )
  {
    if ( npart == 1 && mdcHit1 > m_bhabhaMdcHitCutB && mdcHit2 > m_bhabhaMdcHitCutB )
    {
      isBarrelBhabha = true;
      m_barrelBhabhaNumber++;
    }
    else if ( ( npart == 0 || npart == 2 ) && mdcHit1 > m_bhabhaMdcHitCutE &&
              mdcHit2 > m_bhabhaMdcHitCutE )
    {
      isEndcapBhabha = true;
      m_endcapBhabhaNumber++;
    }
  }
 
  // Select dimu
  if ( m_selectDimu )
  {
    if ( m_dimuTool->IsDimu() == 1 )
    {
      isBarrelDimu = true;
      m_barrelDimuNumber++;
    }
    else if ( m_dimuTool->IsDimu() == 2 )
    {
      isEndcapDimu = true;
      m_endcapDimuNumber++;
    }
  }
 
  // Select hadron
  if ( ( nGood >= 1 && esum > m_hadronChaECut ) || ( nGood == 0 && esum > m_hadronNeuECut ) )
  {
    isHadron = true;
    m_hadronNumber++;
  }
 
  // Select diphoton
  if ( nGood == 0 && eemc > m_diphotonEmcECut && secE > m_diphotonSecECut &&
       fabs( dthe ) < m_diphotonDTheCut && dphi > m_diphotonDPhiCut1 &&
       dphi < m_diphotonDPhiCut2 )
  {
    if ( npart == 1 )
    {
      isBarrelDiphoton = true;
      m_barrelDiphotonNumber++;
    }
    else if ( npart == 0 || npart == 2 )
    {
      isEndcapDiphoton = true;
      m_endcapDiphotonNumber++;
    }
  }
 
  // -------- Write to root file
  if ( m_selectBhabha && isBarrelBhabha )
    m_preSelTools["SelectBarrelBhabha"]->write().ignore();
  if ( m_selectBhabha && isEndcapBhabha )
    m_preSelTools["SelectEndcapBhabha"]->write().ignore();
 
  if ( m_selectDimu && isBarrelDimu ) m_preSelTools["SelectBarrelDimu"]->write().ignore();
  if ( m_selectDimu && isEndcapDimu ) m_preSelTools["SelectEndcapDimu"]->write().ignore();
 
  if ( m_selectHadron && isHadron ) m_preSelTools["SelectHadron"]->write().ignore();
 
  if ( m_selectDiphoton && isBarrelDiphoton )
    m_preSelTools["SelectBarrelDiphoton"]->write().ignore();
  if ( m_selectDiphoton && isEndcapDiphoton )
    m_preSelTools["SelectEndcapDiphoton"]->write().ignore();
 
  if ( m_output )
  {
    m_runnb   = run;
    m_evtnb   = event;
    m_esum    = esum;
    m_eemc    = eemc;
    m_etot    = etot;
    m_nCharge = nCharge;
    m_nGood   = nGood;
    m_nGam    = nGam;
    m_ptot    = ptot;
    m_pp      = pp;
    m_pm      = pm;
    m_maxE    = maxE;
    m_secE    = secE;
    m_dThe    = dthe;
    m_dPhi    = dphi;
    m_mdcHit1 = mdcHit1;
    m_mdcHit2 = mdcHit2;
    m_tuple0->write().ignore();
  }
 
  return StatusCode::SUCCESS;
}

finalize()

StatusCode EventPreSelectAlg::finalize ( )

override

Definition at line 591 of file EventPreSelectAlg.cxx.

                                       {
  info() << "In finalize()" << endmsg;
 
  if ( m_selectDimu && m_output ) { m_dimuTool->Print(); }
 
  always() << "Total events: " << m_events << endmsg;
 
  if ( m_selectBhabha )
  {
    always() << "Selected barrel bhabha: " << m_barrelBhabhaNumber << endmsg;
    always() << "Selected endcap bhabha: " << m_endcapBhabhaNumber << endmsg;
  }
 
  if ( m_selectDimu )
  {
    always() << "Selected barrel dimu: " << m_barrelDimuNumber << endmsg;
    always() << "Selected endcap dimu: " << m_endcapDimuNumber << endmsg;
  }
 
  if ( m_selectHadron ) { always() << "Selected hadron: " << m_hadronNumber << endmsg; }
 
  if ( m_selectDiphoton )
  {
    always() << "Selected barrel diphoton: " << m_barrelDiphotonNumber << endmsg;
    always() << "Selected endcap diphoton: " << m_endcapDiphotonNumber << endmsg;
  }
 
  return StatusCode::SUCCESS;
}

initialize()

StatusCode EventPreSelectAlg::initialize ( )

override

Definition at line 29 of file EventPreSelectAlg.cxx.

                                         {
  info() << "In initialize()" << endmsg;
 
  m_bhabhaEmcECut   = 0.7 * m_ecm;
  m_bhabhaMaxECut   = 0.3 * m_ecm;
  m_bhabhaSecECut   = 0.1 * m_ecm;
  m_dimuEHighCut    = 0.27 * m_ecm;
  m_dimuELowCut     = 0.027 * m_ecm;
  m_diphotonEmcECut = 0.7 * m_ecm;
  m_diphotonSecECut = 0.3 * m_ecm;
  m_hadronChaECut   = 0.3 * m_ecm;
  m_hadronNeuECut   = 0.3 * m_ecm;
 
  if ( m_ecm == 2.4 )
  {
    m_bhabhaEmcECut   = 0.7 * m_ecm;
    m_bhabhaMaxECut   = 0.15 * m_ecm;
    m_bhabhaSecECut   = 0.05 * m_ecm;
    m_diphotonEmcECut = 0.7 * m_ecm;
    m_diphotonSecECut = 0.2 * m_ecm;
  }
 
  info() << "Creating tools..." << endmsg;
 
  // retrieve IEvtWriterTool
  std::vector<std::string> tools_to_retrieve;
  if ( m_selectBhabha )
  {
    tools_to_retrieve.push_back( "SelectBarrelBhabha" );
    tools_to_retrieve.push_back( "SelectEndcapBhabha" );
  }
 
  if ( m_selectDimu )
  {
    tools_to_retrieve.push_back( "SelectBarrelDimu" );
    tools_to_retrieve.push_back( "SelectEndcapDimu" );
  }
 
  if ( m_selectDiphoton )
  {
    tools_to_retrieve.push_back( "SelectBarrelDiphoton" );
    tools_to_retrieve.push_back( "SelectEndcapDiphoton" );
  }
 
  if ( m_selectHadron ) { tools_to_retrieve.push_back( "SelectHadron" ); }
 
  for ( const auto& tool_name : tools_to_retrieve )
  {
    IEventWriterTool* tool = nullptr;
    StatusCode        sc = toolSvc()->retrieveTool( "EventWriterTool", tool_name, tool, this );
    if ( sc.isFailure() )
    {
      error() << "Error retrieving tool " << tool_name << endmsg;
      return sc;
    }
    else info() << "Successfully retrieved tool " << tool_name << endmsg;
 
    m_preSelTools[tool_name] = tool;
  }
 
  // retrieve IDimuPreSelectTool
  if ( m_selectDimu )
  {
    StatusCode sc = toolSvc()->retrieveTool( "DimuPreSelectTool", m_dimuTool, this );
    if ( sc.isFailure() )
    {
      error() << "Error retrieving tool DimuPreSelectTool" << endmsg;
      return sc;
    }
    else info() << "Successfully retrieved tool DimuPreSelectTool" << endmsg;
  }
 
  if ( m_output )
  {
    StatusCode status;
    NTuplePtr  nt0( ntupleSvc(), "FILE1/hadron" );
    if ( nt0 ) m_tuple0 = nt0;
    else
    {
      m_tuple0 = ntupleSvc()->book( "FILE1/hadron", CLID_ColumnWiseTuple, "N-Tuple example" );
      if ( m_tuple0 )
      {
        status = m_tuple0->addItem( "esum", m_esum );
        status = m_tuple0->addItem( "eemc", m_eemc );
        status = m_tuple0->addItem( "etot", m_etot );
        status = m_tuple0->addItem( "nGood", m_nGood );
        status = m_tuple0->addItem( "nCharge", m_nCharge );
        status = m_tuple0->addItem( "nGam", m_nGam );
        status = m_tuple0->addItem( "ptot", m_ptot );
        status = m_tuple0->addItem( "pp", m_pp );
        status = m_tuple0->addItem( "pm", m_pm );
        status = m_tuple0->addItem( "run", m_runnb );
        status = m_tuple0->addItem( "event", m_evtnb );
        status = m_tuple0->addItem( "maxE", m_maxE );
        status = m_tuple0->addItem( "secE", m_secE );
        status = m_tuple0->addItem( "dthe", m_dThe );
        status = m_tuple0->addItem( "dphi", m_dPhi );
        status = m_tuple0->addItem( "mdcHit1", m_mdcHit1 );
        status = m_tuple0->addItem( "mdcHit2", m_mdcHit2 );
      }
      else
      {
        error() << "    Cannot book N-tuple:" << long( m_tuple0 ) << endmsg;
        return StatusCode::FAILURE;
      }
    }
 
    NTuplePtr nt1( ntupleSvc(), "FILE1/vxyz" );
    if ( nt1 ) m_tuple1 = nt1;
    else
    {
      m_tuple1 = ntupleSvc()->book( "FILE1/vxyz", CLID_ColumnWiseTuple, "ks N-Tuple example" );
      if ( m_tuple1 )
      {
        status = m_tuple1->addItem( "vx0", m_vx0 );
        status = m_tuple1->addItem( "vy0", m_vy0 );
        status = m_tuple1->addItem( "vz0", m_vz0 );
        status = m_tuple1->addItem( "vr0", m_vr0 );
        status = m_tuple1->addItem( "theta0", m_theta0 );
        status = m_tuple1->addItem( "p0", m_p0 );
        status = m_tuple1->addItem( "pt0", m_pt0 );
      }
      else
      {
        error() << "    Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
        return StatusCode::FAILURE;
      }
    }
 
    NTuplePtr nt2( ntupleSvc(), "FILE1/photon" );
    if ( nt2 ) m_tuple2 = nt2;
    else
    {
      m_tuple2 =
          ntupleSvc()->book( "FILE1/photon", CLID_ColumnWiseTuple, "ks N-Tuple example" );
      if ( m_tuple2 )
      {
        status = m_tuple2->addItem( "dthe", m_dthe );
        status = m_tuple2->addItem( "dphi", m_dphi );
        status = m_tuple2->addItem( "dang", m_dang );
        status = m_tuple2->addItem( "eraw", m_eraw );
      }
      else
      {
        error() << "    Cannot book N-tuple:" << long( m_tuple2 ) << endmsg;
        return StatusCode::FAILURE;
      }
    }
 
    if ( m_selectDimu )
    {
      NTuplePtr nt3( ntupleSvc(), "FILE1/dimu" );
      if ( nt3 ) m_tuple3 = nt3;
      else
      {
        m_tuple3 =
            ntupleSvc()->book( "FILE1/dimu", CLID_ColumnWiseTuple, "ks N-Tuple example" );
        if ( m_tuple3 ) { m_dimuTool->BookNtuple( m_tuple3 ); }
        else
        {
          error() << "    Cannot book N-tuple:" << long( m_tuple3 ) << endmsg;
          return StatusCode::FAILURE;
        }
      }
    }
  }
 
  m_events               = 0;
  m_barrelBhabhaNumber   = 0;
  m_endcapBhabhaNumber   = 0;
  m_barrelDimuNumber     = 0;
  m_endcapDimuNumber     = 0;
  m_hadronNumber         = 0;
  m_barrelDiphotonNumber = 0;
  m_endcapDiphotonNumber = 0;
 
  info() << "successfully return from initialize()" << endmsg;
  return StatusCode::SUCCESS;
}

---

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

• EventPreSelectAlg.h
• EventPreSelectAlg.cxx