CopyTruthJetParticles.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#include "ParticleJetTools/CopyTruthJetParticles.h"
#include "xAODTruth/TruthVertex.h"
#include "TruthUtils/HepMCHelpers.h"
#include "TruthUtils/MagicNumbers.h"
 
#include "xAODTruth/TruthParticleContainer.h"
#include "AthContainers/ConstDataVector.h"
#include "AsgDataHandles/ReadHandle.h"
#include "AsgDataHandles/WriteHandle.h"
#include "AsgMessaging/Check.h"
#include "AsgTools/CurrentContext.h"
#include "CxxUtils/checker_macros.h"
 
#include <mutex>          // std::call_once, std::once_flag
 
#ifndef XAOD_STANDALONE
// Usage of metadata is for now only possible in Athena...
//#include "CoralBase/AttributeListException.h"
#include "AthAnalysisBaseComps/AthAnalysisHelper.h"
#endif
 
// For std::find in comesFrom()
#include <algorithm>
 
using namespace std;
using namespace MCTruthPartClassifier;
 
CopyTruthJetParticles::CopyTruthJetParticles(const std::string& name)
  : AsgTool(name) {}
 
StatusCode CopyTruthJetParticles::initialize() {
 
  ATH_CHECK(m_truthParticleKey.initialize());
  ATH_CHECK(m_outTruthPartKey.initialize());
  ATH_CHECK(m_dressingNames.initialize());
 
  // ATLASRECTS-8290: this is for backward compatability, remove eventually
  if (m_use_barcode) m_uid = SG::ConstAccessor<int>("barcode");
 
  return StatusCode::SUCCESS;
}
 
 
 
bool CopyTruthJetParticles::classifyJetInput(const xAOD::TruthParticle* tp,
                                             std::vector<const xAOD::TruthParticle*>& promptLeptons,
                                             std::map<const xAOD::TruthParticle*,unsigned int>& tc_results) const {
 
  // Needed for the dressed photon decorations
  const EventContext& ctx = Gaudi::Hive::currentContext();
  
  // Check if this thing is a candidate to be in a truth jet
  //  First block is largely copied from isGenStable, which works on HepMC only
  if (HepMC::is_simulation_particle(tp)) return false; // Particle is from G4
  int pdgid = tp->pdgId();
  if (MC::isZeroEnergyPhoton(tp)) return false; // Work around for an old generator bug
 
  // -- changed for dark jet clustering -- //
  if ( !MC::isStable(tp) && !m_includeDark ) return false; // dark hadrons will not be status 1
  // ----------------------------------- //
  
  // Easy classifiers by PDG ID
  if(MC::isNeutrino(pdgid)) {
    if (!m_includeNu) return false;
  } else {
    if (!m_includeBSMNonInt && !MC::isInteracting(pdgid)) return false;
  }
  if (!m_includeMu && MC::isMuon(pdgid)) return false;
 
  // Already built a list of prompt leptons, just use it here
  if (!m_includePromptLeptons && std::find(promptLeptons.begin(),promptLeptons.end(),tp)!=promptLeptons.end()){
    ATH_MSG_VERBOSE("Veto prompt lepton (" << pdgid << ") with pt " << tp->pt());
    return false;
  }
 
  // Extra catch.  If we aren't supposed to include prompt leptons, we aren't supposed to include prompt neutrinos
  unsigned int tc_res = getTCresult(tp, tc_results);
  if (!m_includePromptLeptons && MC::isNeutrino(pdgid) && MCTruthPartClassifier::isPrompt(tc_res)) {
    return false;
  }
 
  // -- added for dark jet clustering -- //
  // new classifiers to account for dark particles
  // for dark jets: ignore SM particles; include only "stable" dark hadrons
  if (!m_includeSM && !MC::isHiddenValley(tp)) return false;
  if (m_includeDark) {
    if (abs(tp->pdgId()) <= 4900101) return false; // ignore Xd, qd, gd
    if (tp->hasDecayVtx()) {
       size_t good_hadrons = 0;
       for (size_t p = 0; p < tp->end_vertex()->nOutgoingParticles(); ++p) {
           if (!MC::isHiddenValley(tp->child(p))) {
                good_hadrons++; 
           }
       }
       if (good_hadrons == 0) return false; // ignore "non-stable" dark hadrons (decaying to dark sector) -- "stable" if decaying to SM
    }
  }
  // for SM jets: ignore dark particles - probably unnecessary bc of status requirement above
  if (!m_includeDark && MC::isHiddenValley(tp)) return false;
  // ----------------------------------- //
 
  if (!m_includePromptPhotons && MC::isPhoton(pdgid) && tp->hasProdVtx()){
    //ParticleOrigin orig = getPartOrigin(tp, originMap);
    //if (orig==Higgs || orig==HiggsMSSM) return false;
    if (MCTruthPartClassifier::isPrompt(tc_res))  return false;
  }
 
  // If we want to remove photons via the dressing decoration
  if (!m_dressingNames.empty()){
    // Accessor for the dressing decoration above
    bool foundDressDec{false};
    for(const auto &decName : m_dressingNames){
      SG::ReadDecorHandle<xAOD::TruthParticleContainer, char> dressAcc(decName, ctx);
      if (MC::isPhoton(pdgid) && dressAcc(*tp)) foundDressDec = true;
    }
    if (foundDressDec) return false;
  } // End of removal via dressing decoration
 
  // Pseudo-rapidity cut
  if(std::abs(tp->eta())>m_maxAbsEta) return false;
 
  // Vetoes of specific origins.  Not fast, but if no particles are specified should not execute
  if (m_vetoPDG_IDs.size()>0){
    std::vector<int> used_vertices;
    for (int anID : m_vetoPDG_IDs){
      used_vertices.clear();
      if (comesFrom(tp,anID,used_vertices)) return false;
    }
  }
 
  // Made it! 
  return true;
}
 
 
unsigned int CopyTruthJetParticles::getTCresult(const xAOD::TruthParticle* tp,
                                                std::map<const xAOD::TruthParticle*,unsigned int>& tc_results) const
{
  if(tc_results.find(tp) == tc_results.end()) {
    const unsigned int result = tp ? std::get<0>(MCTruthPartClassifier::defOrigOfParticle(tp)) : 0;
    tc_results[tp] = result;
  }
  return tc_results[tp];
}
 
int CopyTruthJetParticles::execute() const {
 
  // retrieve barcode Offset for this event from metadata.
  // We'd need a cleaner solution where this offset is set only at 
  // each new file, but this requires some tool interface which does 
  // not exist in RootCore yet. 
  // So we use the less disruptive solution in Athena for now...
 
  // the function used below is 
   // std::call_once(metaDataFlag,basicMetaDataCheck(), this);
  //  std::call_once(metaDataFlag,this->basicMetaDataCheck());
  // the syntax is explained in http://stackoverflow.com/questions/23197333/why-is-this-pointer-needed-when-calling-stdcall-once
  // this syntax requires the call_once function to receive the object the function is called upon
  //": these are all non-static member functions , they act on objects 
  // they need the "this" pointer which always point to the object the function is working on
  //http://www.learncpp.com/cpp-tutorial/812-static-member-functions/ 
 
  std::vector<const xAOD::TruthParticle*> promptLeptons;
  promptLeptons.reserve(10);
 
  // Retrieve the xAOD truth objects
  auto truthParticles = SG::makeHandle(m_truthParticleKey);
  if ( !truthParticles.isValid() ) {
    ATH_MSG_ERROR("Failed to retrieve truth particle container " << m_truthParticleKey.key());
    return 1;
  }
  
  // Classify particles for tagging and add to the TruthParticleContainer
  std::unique_ptr<ConstDataVector<xAOD::TruthParticleContainer> > ptruth(new ConstDataVector<xAOD::TruthParticleContainer>(SG::VIEW_ELEMENTS));
  std::map<const xAOD::TruthParticle*,unsigned int> tc_results;
  tc_results.clear();
  size_t numCopied = 0;
 
  for (const xAOD::TruthParticle* tp : *truthParticles) {
    if(!tp) continue;
    if (tp->pt() < m_ptmin)
        continue;
    // Cannot use the truth helper functions; they're written for HepMC
    // Last two switches only apply if the thing is a lepton and not a tau
    if ((MC::isElectron(tp) || MC::isMuon(tp)) && tp->hasProdVtx()) {
      // If this is a prompt, generator stable lepton, then we can use it
      if(MC::isStable(tp) && !HepMC::is_simulation_particle(tp) && MCTruthPartClassifier::isPrompt(getTCresult(tp, tc_results))) {
        promptLeptons.push_back(tp);
      }
    }
  }
 
  for (const xAOD::TruthParticle* tp : *truthParticles) {
    if(!tp) continue;
    if (tp->pt() < m_ptmin)
        continue;
 
    if (classifyJetInput(tp, promptLeptons, tc_results)) {
      ptruth->push_back(tp);
      numCopied += 1;
    }
  }
 
  ATH_MSG_DEBUG("Copied " << numCopied << " truth particles into " << m_outTruthPartKey.key() << " TruthParticle container");
 
  // record
  auto truthParticles_out = SG::makeHandle(m_outTruthPartKey);
  ATH_MSG_DEBUG("Recorded truth particle collection " << m_outTruthPartKey.key());
  // notify
  if (!truthParticles_out.put(std::move(ptruth))) {
    ATH_MSG_ERROR("Unable to write new TruthParticleContainer to event store: " 
                  << m_outTruthPartKey.key());
  } else {
    ATH_MSG_DEBUG("Created new TruthParticleContainer in event store: " 
                  << m_outTruthPartKey.key());
  }
 
  return 0;
}
 
 
bool CopyTruthJetParticles::comesFrom( const xAOD::TruthParticle* tp, const int pdgID, std::vector<int>& used_vertices ) const {
  // If it's not a particle, then it doesn't come from something...
  if (!tp) return false;
  // If it doesn't have a production vertex or has no parents, it doesn't come from much of anything
  if (!tp->prodVtx() || tp->nParents()==0) return false;
  // If we have seen it before, then skip this production vertex
  // ATLASRECTS-8290: this should be replaced with ->uid()
  if (std::find(used_vertices.begin(),used_vertices.end(), m_uid(*tp->prodVtx()))!=used_vertices.end()) return false;
  // Add the production vertex to our used list
  // ATLASRECTS-8290: this should be replaced with ->uid()
  used_vertices.push_back( m_uid(*tp->prodVtx()) );
  // Loop over the parents
  for (size_t par=0;par<tp->nParents();++par){
    // Check for null pointers in case of skimming
    if (!tp->parent(par)) continue;
    // Check for a match
    if (tp->parent(par)->absPdgId()==pdgID) return true;
    // Recurse on this parent
    if (comesFrom(tp->parent(par), pdgID, used_vertices)) return true;
  }
  // No hits -- all done with the checks!
  return false;
}