JetTruthParticleSelectorTool.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
#include "JetSimTools/JetTruthParticleSelectorTool.h"
#include "xAODTruth/TruthParticleContainer.h"
 
#include <string>
 
#include <map>
 
#include <functional>
 
 
#include "xAODTruth/TruthVertex.h"
 
#include "JetUtils/JetDistances.h"
 
#include "TruthUtils/MagicNumbers.h"
#include "TruthUtils/HepMCHelpers.h"
 
namespace {
 
  bool isWZDecay(const xAOD::TruthParticle* p) {
    int pdg_id = std::abs(p->pdgId() );
    int mom_pdg_id = pdg_id;
    const xAOD::TruthVertex* vprod = p->prodVtx();
    const xAOD::TruthVertex*oldVprod = vprod;
    // Ascend decay chain looking for when actual decay occurs (not jsut evolution of particle)
    while (pdg_id == mom_pdg_id) {
      const xAOD::TruthParticle* mother = vprod->incomingParticle(0) ;
      if (mother) {
    mom_pdg_id = abs(mother->pdgId());
      } else break;
      if (pdg_id != mom_pdg_id) break;
        
      // Protect against strange infinite reference in sherpa samples
      if(!mother->hasProdVtx()) break;
      if(oldVprod==mother->prodVtx())  break;
      oldVprod = vprod;        
      vprod = mother->prodVtx();
      if (!vprod || vprod->nIncomingParticles() == 0) break;
    } // End while loop
 
      // The following vertex-based identification of W/Z's is needed for SHERPA
      // samples where the W/Z particle is not explicitly in the particle record.
      // At this point if we have a valid vertex, it should be a true decay vertex.
      // If it is a W or Z then two of those decay products should be lepton/neutrino
    int nDecay=0;
    // Prompt W/Z's should come from a vertex with more than one incoming particle
    // This suppresses fave Z's from conversions
    if (vprod && vprod->nIncomingParticles() >1)
      {
    for( const auto& elTruth : vprod->outgoingParticleLinks() ){
      if (MC::isSMLepton(*elTruth)) nDecay++;
    }
      }
 
    bool isWZ = (nDecay==2);
    return (mom_pdg_id == 23 || mom_pdg_id == 24 || isWZ);
  }
 
 
  // This part is copied and adapted from PhysicsAnalysis/D3PDMaker/TruthD3PDMaker/trunk/src/TruthJetFilterTool.cxx
  bool isWZDecay(const xAOD::TruthParticle* p, std::vector<const xAOD::TruthParticle*>& wzLeptons, double photonCone2) {
    int pdg_id = abs(p->pdgId() );
 
    if(! p->hasProdVtx() ) return false;
    const xAOD::TruthVertex* vprod = p->prodVtx();
         
    if (vprod && vprod->nIncomingParticles() > 0) {
      if (isWZDecay(p)) { // paricle decends from W or Z
        //  Save lepton reference for comparison to FSR photons (only for muons and electrons)
        if(MC::isStable(p) && (  (pdg_id==11) || (pdg_id==13) ) ) {
      wzLeptons.push_back(p);
    }
 
        // Only exclude photons within deltaR of leptons (if m_photonCone<0, exclude all photons)
        if( pdg_id == 22 && photonCone2>0)
          {
            for( const auto *lep: wzLeptons) {
              double deltaR2 = jet::JetDistances::deltaR2(*p, *lep);
              if( deltaR2 < photonCone2 ) {
        return true;
          }
            }
          }
        else // not a photon so exclude
          return true;
      }
    }
    return false;
  }
 
 
  // This part is copied and adapted from PhysicsAnalysis/D3PDMaker/TruthD3PDMaker/trunk/src/TruthJetFilterTool.cxx
  bool isLeptonFromTau(const xAOD::TruthParticle* part){
 
    int pdg = part->pdgId();
 
    if(!MC::isSMLepton(pdg)) return false; // all leptons including tau.
    if(!part->hasProdVtx()) return false;
    const xAOD::TruthVertex* prod = part->prodVtx();
    if(!prod) return false; // no parent.
    
    for( const auto& elParent : prod->incomingParticleLinks() ){
      int parentId = (*elParent)->pdgId();
      if(abs(parentId) == 15 && isWZDecay(*elParent)) {
        return true; // Has tau parent
      }
 
      if(parentId == pdg) { // Same particle just a different MC status
        // Go up the generator record until a tau is found or not.
        // Note that this requires a connected *lepton* chain, while calling
        //  isFromTau would allow leptons from hadrons from taus
        if(isLeptonFromTau(*elParent)) {
          return true;
        }
      }
    }
    return false;
  }
 
 
}
 
 
JetTruthParticleSelectorTool::JetTruthParticleSelectorTool(const std::string& n) :
  asg::AsgTool(n),
  m_min_pt(-1),
  m_max_pt(-1),
  m_maxAbsEta(-1),
  m_minEta(-5),
  m_maxEta( 5),
  m_wzPhotonCone(0.1),
  m_useOnlyCharged(false),
  m_listPDGofStables(false),
  m_selectionModeName("StableNoMuonNoNu"),
  m_selectionMode()
{
 
  declareInterface<JetTruthParticleSelectorTool>(this);
  declareProperty("min_pT",                    m_min_pt);
  declareProperty("max_pT",                    m_max_pt);
  declareProperty("max_absEta",                m_maxAbsEta);
  declareProperty("minEta",                    m_minEta);
  declareProperty("maxEta",                    m_maxEta);
  declareProperty("WZModePhotonCone", m_wzPhotonCone);
 
  declareProperty("useChargedParticlesOnly",   m_useOnlyCharged);
  declareProperty("MakeListOfStableParticles", m_listPDGofStables);
  declareProperty("SelectionMode",     m_selectionModeName);
}
 
JetTruthParticleSelectorTool::~JetTruthParticleSelectorTool()
= default;
 
void JetTruthParticleSelectorTool::setupEvent() {
  m_wzLeptons.clear();
}
  
StatusCode
JetTruthParticleSelectorTool::initialize()
{
 
 
  if(     m_selectionModeName == "StableNoMuonNoNu") m_selectionMode = StableNoMuonNoNu ;
  else if(m_selectionModeName == "MuonOnly") m_selectionMode = MuonOnly ;
  else if(m_selectionModeName == "NuOnly") m_selectionMode = NuOnly ;
  else if(m_selectionModeName == "NoWZDecay") m_selectionMode = NoWZDecay ;
  else {
    ATH_MSG_ERROR(" Unknown truth selection mode "<< m_selectionMode);
    return StatusCode::FAILURE;
  }
 
  m_pdgList.clear();
  m_avP.clear();
  m_av2P.clear();
  m_avPt.clear();
  m_av2Pt.clear();
  m_avEta.clear();
  m_av2Eta.clear();
  m_avPhi.clear();
  m_av2Phi.clear();
  m_avM.clear();
  m_av2M.clear();
  
  return StatusCode::SUCCESS;
}
 
 
bool JetTruthParticleSelectorTool::passKinematics(const xAOD::TruthParticle* truthPart) const{
  auto truth4Mom = truthPart->p4();
  // some safety checks ... both will crash the eta calculation
  if ( truth4Mom.E() <= 0 )                                return false;
  if ( truth4Mom.P() - truth4Mom.Pz() <= 0 )                     return false;
  
  double pt  = truth4Mom.Pt();
  if ( m_min_pt > 0 && pt < m_min_pt )               return false;
  if ( m_max_pt > 0 && pt > m_max_pt )               return false;
  double eta = truth4Mom.Eta();
  if ( m_maxAbsEta > 0 && fabs(eta) > m_maxAbsEta )  return false;
  if ( eta < m_minEta )                              return false;
  if ( eta > m_maxEta )                              return false;
  return true;
}
 
bool JetTruthParticleSelectorTool::selector(const xAOD::TruthParticle* truthPart) {
 
  bool result = false;
  if (MC::isGeantino(truthPart))  return false; // protect against unexpected geantinos
 
  switch( m_selectionMode) { // For now only 4 modes used in practice for jets...
                             // a switch statement is probably not optimal here...
  case StableNoMuonNoNu:
    result =(   MC::isStableOrSimDecayed(truthPart) && !HepMC::is_simulation_particle(truthPart) &&  !MC::isZeroEnergyPhoton(truthPart) &&
     passKinematics(truthPart) && !MC::isMuon(truthPart) && (!( MC::isStable(truthPart) && MC::isSpecialNonInteracting(truthPart))));
    break;
  case NuOnly:
    result= ( MC::isStableOrSimDecayed(truthPart) && !HepMC::is_simulation_particle(truthPart) &&  !MC::isZeroEnergyPhoton(truthPart) &&
     passKinematics(truthPart) &&  MC::isStable(truthPart) && MC::isSpecialNonInteracting(truthPart));
    break;
  case MuonOnly:
    result = MC::isMuon(truthPart);
    break;
  case NoWZDecay:
    result = (MC::isStableOrSimDecayed(truthPart) && !HepMC::is_simulation_particle(truthPart) &&  !MC::isZeroEnergyPhoton(truthPart) && 
               passKinematics(truthPart) &&
               !isWZDecay(truthPart, m_wzLeptons, m_wzPhotonCone*m_wzPhotonCone) &&
               !isLeptonFromTau(truthPart) );
    break;
  default:
    break;
  }
  
  if (m_listPDGofStables && result ) {
    int apId = abs(truthPart->pdgId());
    auto truth4Mom = truthPart->p4();
    double p = truth4Mom.P();
    double phi = truth4Mom.Phi();
    double pt  = truth4Mom.Pt();
    double eta  = truth4Mom.Eta();
    PDGList::iterator it = m_pdgList.find(apId);
    if ( it == m_pdgList.end() ) {
      m_pdgList[apId]=1;
      m_avP[apId]=p;
      m_av2P[apId]=p * p;
      m_avPt[apId]=pt;
      m_av2Pt[apId]=pt * pt;
      m_avEta[apId]=eta;
      m_av2Eta[apId]=eta * eta;
      m_avPhi[apId]=phi;
      m_av2Phi[apId]=phi * phi;
      m_avM[apId]=truthPart->m();
      m_av2M[apId]=truthPart->m() * truthPart->m();
    } else {
      it->second++;
      m_avP[apId]+=p;
      m_av2P[apId]+=p * p;
      m_avPt[apId]+=pt;
      m_av2Pt[apId]+=pt * pt;
      m_avEta[apId]+=eta;
      m_av2Eta[apId]+=eta * eta;
      m_avPhi[apId]+=phi;
      m_av2Phi[apId]+=phi * phi;
      m_avM[apId]+=truthPart->m();
      m_av2M[apId]+=truthPart->m() * truthPart->m();
    }
  }
 
  ATH_MSG_VERBOSE("Truth selection result: " << truthPart<<", result=" << result);
 
  return result;
}
 
StatusCode
JetTruthParticleSelectorTool::finalize()
{
  if (m_listPDGofStables)
    {
      msg(MSG::INFO) << "Counts of PDGs of all stable particles :" << endmsg;
      msg(MSG::INFO) << "========================================" << endmsg;
      msg(MSG::INFO) << "|    PDG    |  # particles  |" << endmsg;
      for ( PDGList::iterator it = m_pdgList.begin(); it != m_pdgList.end(); ++it )
        {
          msg(MSG::INFO) << "|"
                         << std::setw(10) << it->first << " |"
                         << std::setw(10) << it->second << " |"
                         << endmsg;
        }
      msg(MSG::INFO) << "|    PDG    |    <p>    |    rms    |    <pt>    |    rms    |    <eta>    |    rms    |" << endmsg;
      for ( PDGList::iterator it = m_pdgList.begin(); it != m_pdgList.end(); ++it )
        {
          int n=it->second;
          double p=m_avP[it->first]/n;
          double p2=std::sqrt(m_av2P[it->first]/n-
                         m_avP[it->first]*m_avP[it->first]/n/n);
          double pt=m_avPt[it->first]/n;
          double pt2=std::sqrt(m_av2Pt[it->first]/n-
                          m_avPt[it->first]*m_avPt[it->first]/n/n);
          double eta=m_avEta[it->first]/n;
          double eta2=std::sqrt(m_av2Eta[it->first]/n-
                           m_avEta[it->first]*m_avEta[it->first]/n/n);
          msg(MSG::INFO) << "|"
                         << std::setw(15) << it->first << " |"
                         << std::setw(15) << p << " |"
                         << std::setw(15) << p2 << " |"
                         << std::setw(15) << pt << " |"
                         << std::setw(15) << pt2 << " |"
                         << std::setw(15) << eta << " |"
                         << std::setw(15) << eta2 << " |"
                         << endmsg;
        }
      msg(MSG::INFO) << "|    PDG    |   <phi>   |    rms    |     <m>     |     rms    |" << endmsg;
      for ( PDGList::iterator it = m_pdgList.begin(); it != m_pdgList.end(); ++it )
        {
          int n=it->second;
          double phi=m_avPhi[it->first]/n;
          double phi2=std::sqrt(m_av2Phi[it->first]/n-
                           m_avPhi[it->first]*m_avPhi[it->first]/n/n);
          double m=m_avM[it->first]/n;
          double m2=0;
          if ( m_av2M[it->first]/n > m_avM[it->first]*m_avM[it->first]/n/n)
            m2=std::sqrt(m_av2M[it->first]/n-
                    m_avM[it->first]*m_avM[it->first]/n/n);
          msg(MSG::INFO) << "|"
                         << std::setw(15) << it->first << " |"
                         << std::setw(15) << phi << " |"
                         << std::setw(15) << phi2 << " |"
                         << std::setw(15) << m << " |"
                         << std::setw(15) << m2 << " |"
                         << endmsg;
        }
    }
  return StatusCode::SUCCESS;
}