HepMC_helper.cxx

Go to the documentation of this file.

#include "AtlfastUtils/HepMC_helper/HepMC_helper.h"
#include "AtlfastEvent/ParticleCodes.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "HepMC/GenEvent.h"
#include "AtlfastEvent/ChargeService.h"

#include <cassert>
namespace HepMC_helper {
  using std::vector;
  using std::abs;
  //****************************************************************
  //*                    All                                       *
  //****************************************************************
  bool All::operator()( const Particle* const) const {
    return true;
  }
  bool All::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 
  IMCselector* All::create() const {return new All();}
  //****************************************************************
  //*                    IsFinalState                              *
  //***********************final*****************************************
  IsFinalState::IsFinalState()
    : m_igs() {}
  
  bool IsFinalState::operator()( const Particle* const p ) const {
    if ( !m_igs(p) ) return false;
    if ( m_chargeService.operator()(p) == -999. ){
      std::cout << "HepMC_helper::IsFinalState ignoring this one, particle will NOT be fast-simulated" << std::endl;
      return false;
    }    
    return true;
  }
  
  bool IsFinalState::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 
  
  IMCselector* IsFinalState::create() const {return new IsFinalState(*this);}

  //****************************************************************
  //*                    IsStatusxx                              *
  //***********************final**************************************
  bool IsStatusxx::operator()( const Particle* const p ) const {
    const int fsCode = 3;

    int pStatus1 = p->status();
    int pStatus2 = pStatus1%100;
    
    if(m_stat != fsCode){return ( pStatus2 == m_stat)? true:false;}
    
    //This IsStatus has been initialised with the "usual" final state code.
    
    //Apply herwig history if the generator is Herwig. 

    int generCode  = p->parent_event()->signal_process_id()/1000000;
    
    //Herwig
    if (generCode == 2){
      return ((pStatus1 >= 101 && pStatus1 <= 103) ||
              (pStatus1 >= 120 && pStatus1 <= 125)) ? true:false;
    }

    //All other generators.....
    return ( pStatus2 == fsCode)? true:false;

  }
  bool IsStatusxx::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 
  IMCselector* IsStatusxx::create() const {return new IsStatusxx(*this);}
  //****************************************************************
  //*                    IsFromHardScatter                           *
  //*****************************************************************
  bool IsFromHardScatter::operator()( const Particle* const /*p*/ ) const {

    // This just returns true for the moment
    // Expect to be able to query GenEvent for status in future
    // p->parent_event()
    
    return true;

  }
  bool IsFromHardScatter::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 
  IMCselector* IsFromHardScatter::create() const {return new IsFromHardScatter(*this);}
  //****************************************************************
  //* RMS                  NCutter                                 *
  //****************************************************************
  bool NCutter::operator()( const Particle* const p ) const {
    for (vector<IMCselector*>::const_iterator i = m_selectors.begin();
         i != m_selectors.end(); i++) { 
      if ( !(*i)->operator()(p) ) return false;
    }
    return true;
  }
  bool NCutter::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 
  IMCselector* NCutter::create() const {return new NCutter(*this);}
  //****************************************************************
  //*                    SelectType                                *
  //****************************************************************
  SelectType::SelectType(vector<int> types):m_requiredTypes(types){}
  SelectType::SelectType(int type){m_requiredTypes.push_back(type);}
  SelectType::SelectType(const SelectType& rhs):
    IMCselector(),
    m_requiredTypes(rhs.m_requiredTypes){}

  bool SelectType::operator() ( const Particle* const p ) const {
    vector<int>::const_iterator itype = m_requiredTypes.begin();
    for(; itype < m_requiredTypes.end(); ++itype ){ 
      if( abs(p->pdg_id()) == (*itype) ){
        return true ;
      }
    }
    return false ;
  }
  bool SelectType::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 
  IMCselector* SelectType::create() const {return new SelectType(*this);}
  
  //****************************************************************
  //*                    RejectType                                *
  //****************************************************************
  
  RejectType::RejectType(int type): m_selectType(type){}
  RejectType::RejectType(vector<int> types): m_selectType(types){}
  RejectType::RejectType(const RejectType& rhs):
    IMCselector(), 
    m_selectType(rhs.m_selectType){}

  bool RejectType::operator() ( const Particle* const p ) const{
      return !m_selectType(p);}
  bool RejectType::operator()( const Particle& p ) const {
      return this->operator()(&p);} 
  IMCselector* RejectType::create() const {return new RejectType(*this);}

  //****************************************************************
  //*                     IsCharged                                *
  //****************************************************************


  IsCharged::IsCharged(){}

  bool IsCharged::operator() ( const Particle* const p )const{
    
    double charge = m_chargeService.operator()(p);
    if ( charge == 0 ) return false;
    if ( charge == -999. ){
      std::cout << "HepMC_helper::IsCharged ignoring this one" << std::endl;
      return false;
    }
    return true ;       
  }

  bool IsCharged::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 

  IMCselector* IsCharged::create() const {
    IMCselector* p = 0;
    try{
      p = new IsCharged(*this);
    }catch(std::string errMsg){
      throw errMsg;
    }
    return p;
  }

  //****************************************************************
  //*                     MCCuts                                   *
  //****************************************************************
  MCCuts::MCCuts(double ptMin, double etaMax): 
    m_ptMin(ptMin),m_etaMax(etaMax){}
  bool MCCuts::operator() ( const Particle* const p )const {
    if (p->momentum().perp() < m_ptMin) return false;
    if (abs(p->momentum().pseudoRapidity()) > m_etaMax) return false;
    return true;
  }

  bool MCCuts::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 
  IMCselector* MCCuts::create() const  {return new MCCuts(*this);}
  

  //****************************************************************
  //*                    Unseen                                    *
  //****************************************************************
  Unseen::Unseen( vector<int> requiredTypes,
                  double muonPtMax,
                  double muonEtaMax): 
    m_seltype( SelectType(requiredTypes) ),
    m_muonPtMax( muonPtMax ),
    m_muonEtaMax( muonEtaMax ){
  }
  
  Unseen::Unseen(const Unseen& src ): 
    IMCselector(),
    m_seltype( src.m_seltype ),
    m_muonPtMax( src.m_muonPtMax ),
    m_muonEtaMax( src.m_muonEtaMax ){ 
  }
  bool Unseen::operator() ( const Particle* const p )const{
    if(m_seltype(p)) return true;
    if( abs(p->pdg_id()) == 13){
      if(!m_ifs(p)) return false;
      if(abs(p->momentum().pseudoRapidity()) > m_muonEtaMax) return true;
      if(    p->momentum().perp()            < m_muonPtMax ) return true;
    }
    return false ;       
  }
  IMCselector* Unseen::create() const {return new Unseen(*this);}
  // The operator() method to determine acceptability of particle
  bool Unseen::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 
  //****************************************************************
  //*                    SelectJetTag                              *
  //****************************************************************
  
  SelectJetTag::SelectJetTag(int id, double pt, double eta):
    m_id(abs(id)), m_ptMin(pt), m_etaMax(eta) { }
  bool SelectJetTag::operator() ( const Particle* const p ) const{
    // check on kinematics
    if( abs(p->pdg_id())                    !=     m_id) return false ; 
    if( p->momentum().perp()                 <  m_ptMin) return false ; 
    if((abs(p->momentum().pseudoRapidity())) > m_etaMax) return false ; 

    // No documentary quark
    if (p->status() == 3)                                return false ;
    // No quark from B/D decays (Herwig)
    bool fromHadron        = false;
    HepMC::GenVertex *pvtx = p->production_vertex();
    if (pvtx) {
      HepMC::GenVertex::particle_iterator firstParent = pvtx->particles_begin(HepMC::ancestors);
      HepMC::GenVertex::particle_iterator endParent   = pvtx->particles_end(HepMC::ancestors);
      HepMC::GenVertex::particle_iterator thisParent  = firstParent;
      for(; thisParent != endParent; ++thisParent){
       int pdgAnces = (*thisParent)->pdg_id();
       if (isBBaryon(pdgAnces) || isBMeson(pdgAnces) || // B --> b,c
           isDBaryon(pdgAnces) || isDMeson(pdgAnces) )  // D --> c
         fromHadron = true;
       if (fromHadron) break;
      }
    }
    if (fromHadron) return false;

    if(p->end_vertex()){
      
      HepMC::GenVertex::particle_iterator firstChild = 
        p->end_vertex()->particles_begin(HepMC::children);
      
      HepMC::GenVertex::particle_iterator endChild = 
        p->end_vertex()->particles_end(HepMC::children);
      
      HepMC::GenVertex::particle_iterator thisChild = firstChild; 
      
      //Should be vetoing this tag if it has a child of the same pdg id
      // following kludge due to fact that parents are included in
      // the list of children!
      for(; thisChild!=endChild; ++thisChild){
        if(abs((*thisChild)->pdg_id()) == m_id) return false;
      }
    }
    return true ;       
  }
  IMCselector* SelectJetTag::create() const {return new SelectJetTag(*this);}
  bool SelectJetTag::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 

  //Taken from BSignalFilter
  bool SelectJetTag::isBBaryon(const int pID) const {
    // PdgID of B-baryon is of form ...xxx5xxx
    std::string idStr = longToStr( abs(pID) );
    char digit4 = idStr[ idStr.length() - 4 ];
    if( (digit4=='5') ) 
      return true;
    else 
      return false;
  }

  bool SelectJetTag::isBMeson(const int pID) const {
    // PdgID of B-meson is of form ...xxx05xx
    std::string idStr = longToStr( abs(pID) );
    char digit3 = idStr[ idStr.length() - 3 ];
    char digit4;
    if( idStr.length() < 4 ) { digit4 = '0'; }
    else { digit4 = idStr[ idStr.length() - 4 ]; };
    if( (digit4=='0') && (digit3=='5') ) 
      return true;
    else 
      return false;
  }
    bool SelectJetTag::isDBaryon(const int pID) const {
      // PdgID of D-baryon is of form ...xxx4xxx
      std::string idStr = longToStr( abs(pID) );
      char digit4 = idStr[ idStr.length() - 4 ];
      if( (digit4=='4') ) 
        return true;
      else 
        return false;
    }

    bool SelectJetTag::isDMeson(const int pID) const {
      // PdgID of D-meson is of form ...xxx04xx
      std::string idStr = longToStr( abs(pID) );
      char digit3 = idStr[ idStr.length() - 3 ];
      char digit4;
      if( idStr.length() < 4 ) { digit4 = '0'; }
      else { digit4 = idStr[ idStr.length() - 4 ]; };
      if( (digit4=='0') && (digit3=='4') ) 
        return true;
      else 
        return false;
    }

    std::string SelectJetTag::longToStr( const long n ) const {
      if (0==n) return "0"; 
      std::string str = "";
      for ( long m = n; m!=0; m/=10 )
        str = char( '0' + abs(m%10) ) + str;
      if ( n<0 )
        str = "-" + str;
      return str;
    }

  //****************************************************************
  //*                    SelectTauTag                              *
  //****************************************************************

  SelectTauTag::SelectTauTag(double pt, double eta):
    m_jetTagSelector(ParticleCodes::TAU, pt, eta){
  }
  bool SelectTauTag::operator() ( const Particle* const p ) const{
    if(p->end_vertex()){
      
      HepMC::GenVertex::particle_iterator firstChild = 
        p->end_vertex()->particles_begin(HepMC::children);
      
      HepMC::GenVertex::particle_iterator endChild = 
        p->end_vertex()->particles_end(HepMC::children);
      
      HepMC::GenVertex::particle_iterator thisChild = firstChild; 
      
      // Reject if the Tau does not decay hadronically
      for(; thisChild!=endChild; ++thisChild){
        if(abs((*thisChild)->pdg_id()) == 11 || 
           abs((*thisChild)->pdg_id()) == 13 || 
           abs((*thisChild)->pdg_id()) == 15 ) return false;
      }
    }
    return m_jetTagSelector(p);
  }
  
  IMCselector* SelectTauTag::create() const {return new SelectTauTag(*this);}
  
  bool SelectTauTag::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 

  //****************************************************************
  //*                    SelectZ0                                  *
  //****************************************************************

  bool SelectZ0::operator() ( const Particle* const p ) const{
    if(!m_z0Type(p))   return false;
    if(!m_kineCuts(p)) return false;
    return true;
  } 
  IMCselector* SelectZ0::create() const {return new SelectZ0(*this);}
  bool SelectZ0::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 


  //****************************************************************
  //*                    SelectPid                                 *
  //****************************************************************

  bool SelectPid::operator() ( const Particle* const p ) const{
    if(!m_selType(p))  return false;
    if(!m_kineCuts(p)) return false;
    return true;
  } 
  IMCselector* SelectPid::create() const {return new SelectPid(*this);}
  bool SelectPid::operator() ( const Particle& p ) const {
    return this->operator()(&p);
  } 

  //****************************************************************
  //*                    BFieldCutter                              *
  //****************************************************************
  
  bool BFieldCutter::operator()(const Particle* const a) const {
    if(m_chargeService->operator()(a) == 0.) return true;
    return (a->momentum().perp()<m_ptCut) ? false:true;
  }
   
  IMCselector* BFieldCutter::create() const {return new BFieldCutter(*this);}
  bool BFieldCutter::operator()( const Particle& p ) const {
    return this->operator()(&p);
  } 

  //****************************************************************
  //*                    AscendingEta                              *
  //****************************************************************

  bool AscendingEta::operator() ( const Particle* const a, 
                                  const Particle* const b  ) const{ 
    return( a->momentum().pseudoRapidity() < 
            b->momentum().pseudoRapidity() ); }
  bool AscendingEta::operator() ( const Particle& a, 
                                  const Particle& b  ) const{ 
    return( a.momentum().pseudoRapidity() < 
            b.momentum().pseudoRapidity() ); 
  }

}  // closing the namespace

---