ATLAS Offline Software
TruthDressingTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // TruthDressingTool.cxx
7 // Author: Kevin Finelli (kevin.finelli@cern.ch)
8 // Create dressed (i.e. including FSR photons) 4-vectors of truth objects
9 
17 
19 
20 #include "fastjet/PseudoJet.hh"
21 #include "fastjet/ClusterSequence.hh"
22 #include <vector>
23 #include <string>
24 #include <algorithm>
25 #include <memory>
26 namespace {
27  static const SG::ConstAccessor<unsigned int> acc_origin("Classification");
28 }
29 // Constructor
31  const std::string& n,
32  const IInterface* p )
33  : base_class(t,n,p)
34 {
35 }
36 
37 // Destructor
39 }
40 
41 // Athena initialize and finalize
43 {
44  // Initialise handle keys
45  ATH_CHECK(m_particlesKey.initialize());
46  ATH_CHECK(m_dressParticlesKey.initialize());
47  m_decorator_eKey = m_dressParticlesKey.key() + ".e_dressed";
48  ATH_CHECK(m_decorator_eKey.initialize());
49  m_decorator_ptKey = m_dressParticlesKey.key() + ".pt_dressed";
50  ATH_CHECK(m_decorator_ptKey.initialize());
51  m_decorator_etaKey= m_dressParticlesKey.key() + ".eta_dressed";
52  ATH_CHECK(m_decorator_etaKey.initialize());
53  m_decorator_phiKey= m_dressParticlesKey.key() + ".phi_dressed";
54  ATH_CHECK(m_decorator_phiKey.initialize());
55  m_decorator_pt_visKey = m_dressParticlesKey.key() + ".pt_vis_dressed";
56  ATH_CHECK(m_decorator_pt_visKey.initialize());
57  m_decorator_eta_visKey= m_dressParticlesKey.key() + ".eta_vis_dressed";
58  ATH_CHECK(m_decorator_eta_visKey.initialize());
59  m_decorator_phi_visKey= m_dressParticlesKey.key() + ".phi_vis_dressed";
60  ATH_CHECK(m_decorator_phi_visKey.initialize());
61  m_decorator_m_visKey= m_dressParticlesKey.key() + ".m_vis_dressed";
62  ATH_CHECK(m_decorator_m_visKey.initialize());
63  m_decorator_nphotonKey = m_dressParticlesKey.key() + ".nPhotons_dressed";
64  ATH_CHECK(m_decorator_nphotonKey.initialize());
65  if (!m_decorationName.empty()) {m_decorationKey = m_particlesKey.key()+"."+m_decorationName;}
66  else {m_decorationKey = m_particlesKey.key()+".unusedPhotonDecoration";}
67  ATH_CHECK(m_decorationKey.initialize());
68  m_truthClassKey = m_dressParticlesKey.key() + "." + SG::AuxTypeRegistry::instance().getName(acc_origin.auxid());
69  ATH_CHECK(m_truthClassKey.initialize());
70 
71  // ensure we are not mixing truth taus with other truth particles
72  if (m_listOfPIDs.size() > 1) {
73  if (std::find(m_listOfPIDs.begin(), m_listOfPIDs.end(), 15) != m_listOfPIDs.end()) {
74  ATH_MSG_ERROR("Truth taus must be dressed separately from other truth particles");
75  return StatusCode::FAILURE;
76  }
77  }
78 
79  return StatusCode::SUCCESS;
80 }
81 
82 // Function to do dressing, implements interface in IAugmentationTool
84 {
85  // Get the event context
86  const EventContext& ctx = Gaudi::Hive::currentContext();
87 
88  // Retrieve the truth collections
89  SG::ReadHandle<xAOD::TruthParticleContainer> truthParticles(m_particlesKey,ctx);
90  if (!truthParticles.isValid()) {
91  ATH_MSG_ERROR("Couldn't retrieve TruthParticle collection with name " << m_particlesKey);
92  return StatusCode::FAILURE;
93  }
94 
95  SG::ReadHandle<xAOD::TruthParticleContainer> dressTruthParticles(m_dressParticlesKey,ctx);
96  if (!dressTruthParticles.isValid()) {
97  ATH_MSG_ERROR("Couldn't retrieve TruthParticle collection with name " << m_dressParticlesKey);
98  return StatusCode::FAILURE;
99  }
100 
101  // Decorators
102  SG::WriteDecorHandle< xAOD::TruthParticleContainer,float > decorator_e(m_decorator_eKey, ctx);
103  SG::WriteDecorHandle< xAOD::TruthParticleContainer,float > decorator_pt(m_decorator_ptKey, ctx);
104  SG::WriteDecorHandle< xAOD::TruthParticleContainer,float > decorator_eta(m_decorator_etaKey, ctx);
105  SG::WriteDecorHandle< xAOD::TruthParticleContainer,float > decorator_phi(m_decorator_phiKey, ctx);
106  // for truth taus, use 'vis' in the decoration name to avoid prompt/visible tau momentum ambiguity
107  // use (pt,eta,phi,m) for taus, for consistency with other TauAnalysisTools decorations
108  SG::WriteDecorHandle< xAOD::TruthParticleContainer,float > decorator_pt_vis(m_decorator_pt_visKey, ctx);
109  SG::WriteDecorHandle< xAOD::TruthParticleContainer,float > decorator_eta_vis(m_decorator_eta_visKey, ctx);
110  SG::WriteDecorHandle< xAOD::TruthParticleContainer,float > decorator_phi_vis(m_decorator_phi_visKey, ctx);
111  SG::WriteDecorHandle< xAOD::TruthParticleContainer,float > decorator_m_vis(m_decorator_m_visKey, ctx);
112  SG::WriteDecorHandle< xAOD::TruthParticleContainer,int > decorator_nphoton(m_decorator_nphotonKey, ctx);
113  SG::WriteDecorHandle< xAOD::TruthParticleContainer,char > dressDec(m_decorationKey, ctx);
114  // If we want to decorate, then we need to decorate everything with false to begin with
115  if (!m_decorationName.empty()){
116  for (const auto * particle : *truthParticles){
117  dressDec(*particle) = 0;
118  }
119  } // We are using the decoration
120 
121  // accessors for truth tau visible momentum
122  static const SG::ConstAccessor<double> pt_visAcc("pt_vis");
123  static const SG::ConstAccessor<double> eta_visAcc("eta_vis");
124  static const SG::ConstAccessor<double> phi_visAcc("phi_vis");
125  static const SG::ConstAccessor<double> mvisAcc("m_vis");
126 
127  //get struct of helper functions
129 
130  std::vector<const xAOD::TruthParticle*> listOfParticlesToDress;
131  std::vector<xAOD::TruthParticle::FourMom_t> listOfDressedP4;
132  std::vector<xAOD::TruthParticle::FourMom_t> listOfBareP4;
133 
134  if (m_listOfPIDs.size()==1 && std::abs(m_listOfPIDs[0])==15) {
135  // when dressing truth taus, it is assumed that the truth tau container has been built beforehand and is used as input
136  for (auto *pItr : *dressTruthParticles) {
137  if (!pItr->isTau()) {
138  ATH_MSG_ERROR("Input particles should be truth taus.");
139  return StatusCode::FAILURE;
140  }
141  if (!pt_visAcc.isAvailable(*pItr) || !eta_visAcc.isAvailable(*pItr) || !phi_visAcc.isAvailable(*pItr) || !mvisAcc.isAvailable(*pItr)) {
142  ATH_MSG_ERROR("Visible momentum not available for truth taus, cannot perform dressing!");
143  return StatusCode::FAILURE;
144  }
145 
146  listOfParticlesToDress.push_back(pItr);
147 
148  // we dresss the visible 4-momentum
150  bare_part.SetPtEtaPhiM(pt_visAcc(*pItr), eta_visAcc(*pItr), phi_visAcc(*pItr), mvisAcc(*pItr));
151  listOfDressedP4.push_back(bare_part);
152  }
153  // in the cone-based approach, make a copy of bare P4 to avoid recomputing it
154  if (!m_useAntiKt) {
155  listOfBareP4 = listOfDressedP4;
156  }
157  }
158  else {
159  // non-prompt particles are still included here to ensure all particles
160  // will get the decoration; however further down only the prompt particles
161  // are actually dressed depending on the value of m_useLeptonsFromHadrons
162  decayHelper.constructListOfFinalParticles(dressTruthParticles.ptr(), listOfParticlesToDress, m_listOfPIDs, true);
163 
164  for (const auto* part : listOfParticlesToDress) {
165  listOfDressedP4.push_back(part->p4());
166  }
167  }
168 
169  std::vector<int> dressedParticlesNPhot(listOfParticlesToDress.size(), 0);
170 
171  //fill the photon list
172  std::vector<const xAOD::TruthParticle*> photonsFSRList;
173  std::vector<int> photonPID{22};
174  const bool pass = decayHelper.constructListOfFinalParticles(truthParticles.ptr(), photonsFSRList,
175  photonPID, m_usePhotonsFromHadrons);
176  if (!pass) {
177  ATH_MSG_WARNING("Cannot construct the list of final state particles "<<m_truthClassKey.fullKey());
178  }
179 
180  // Do dR-based photon dressing (default)
181  if (!m_useAntiKt){
182  //loop over photons, uniquely associate each to nearest bare particle
183  for (const auto* phot : photonsFSRList) {
184  double dRmin = m_coneSize;
185  int idx = -1;
186 
187  for (size_t i = 0; i < listOfParticlesToDress.size(); ++i) {
188  if (!m_useLeptonsFromHadrons) {
189  if (!acc_origin.isAvailable(*listOfParticlesToDress[i])) {
190  ATH_MSG_WARNING("MCTruthClassifier "<<m_truthClassKey.fullKey() <<" not available, cannot apply notFromHadron veto!");
191  }
192  unsigned int result = acc_origin(*listOfParticlesToDress[i]);
193  const bool isPrompt = MCTruthPartClassifier::isPrompt(result, true);
194  if (!isPrompt) continue;
195  }
197  if (listOfParticlesToDress[i]->isTau()) {
198  bare_part = listOfBareP4[i];
199  }
200  else {
201  bare_part = listOfParticlesToDress[i]->p4();
202  }
203 
204  double dR = bare_part.DeltaR(phot->p4());
205  if (dR < dRmin) {
206  dRmin = dR;
207  idx = i;
208  }
209  }
210 
211  if(idx > -1) {
212  listOfDressedP4[idx] += phot->p4();
213  dressedParticlesNPhot[idx]++;
214  if (!m_decorationName.empty()){
215  dressDec(*phot) = 1;
216  }
217  }
218  }
219 
220  //loop over particles and add decorations
221  for (size_t i = 0; i < listOfParticlesToDress.size(); ++i) {
222  const xAOD::TruthParticle* part = listOfParticlesToDress[i];
223  const xAOD::TruthParticle::FourMom_t& dressedVec = listOfDressedP4[i];
224 
225  if(part->isTau()) {
226  decorator_pt_vis(*part) = dressedVec.Pt();
227  decorator_eta_vis(*part) = dressedVec.Eta();
228  decorator_phi_vis(*part) = dressedVec.Phi();
229  decorator_m_vis(*part) = dressedVec.M();
230  }
231  else {
232  decorator_e(*part) = dressedVec.E();
233  decorator_pt(*part) = dressedVec.Pt();
234  decorator_eta(*part) = dressedVec.Eta();
235  decorator_phi(*part) = dressedVec.Phi();
236  }
237  decorator_nphoton(*part) = dressedParticlesNPhot[i];
238  }
239  } // end of the dR matching part
240 
241  //build the anti-kt jet list
242  if (m_useAntiKt) {
243  std::vector<fastjet::PseudoJet> sorted_jets;
244  std::vector<fastjet::PseudoJet> fj_particles;
245  for (const auto* part : listOfParticlesToDress) {
246  if(part->isTau()) {
247  const xAOD::TruthParticle::FourMom_t& tauvis = listOfDressedP4[part->index()];
248  fj_particles.emplace_back(tauvis.Px(), tauvis.Py(), tauvis.Pz(), tauvis.E());
249  }
250  else {
251  fj_particles.emplace_back(part->px(), part->py(), part->pz(), part->e());
252  }
253 
254  fj_particles.back().set_user_index(HepMC::uniqueID(part));
255  }
256  for (const auto* part : photonsFSRList) {
257  fj_particles.emplace_back(part->px(), part->py(), part->pz(), part->e());
258  fj_particles.back().set_user_index(HepMC::uniqueID(part));
259  }
260 
261  //run the clustering
262  fastjet::JetAlgorithm alg=fastjet::antikt_algorithm;
263  const fastjet::JetDefinition jet_def(alg, m_coneSize);
264  fastjet::ClusterSequence cseq(fj_particles, jet_def);
265  sorted_jets = sorted_by_pt(cseq.inclusive_jets(0));
266  //associate clustered jets back to bare particles
267  std::vector<int> photon_uniqueIDs(50);
268  photon_uniqueIDs.clear();
269  for (const auto* part : listOfParticlesToDress) {
270  //loop over fastjet pseudojets and associate one with this particle
271  bool found=false;
272  auto pjItr=sorted_jets.begin();
273  while(!found && pjItr!=sorted_jets.end()) {
274  std::vector<fastjet::PseudoJet> constituents = pjItr->constituents();
275  for(const auto& constit : constituents) {
276  if (HepMC::uniqueID(part)==constit.user_index()) {
277 
278  // shall we count the number of photons among the pseudojet constituents
279  // to decorate leptons with the number of dressing photons found by the anti-kt algorithm?
280 
281  if(part->isTau()) {
282  decorator_pt_vis(*part) = pjItr->pt();
283  decorator_eta_vis(*part) = pjItr->pseudorapidity();
284  decorator_phi_vis(*part) = pjItr->phi_std(); //returns phi in [-pi,pi]
285  decorator_m_vis(*part) = pjItr->m();
286  }
287  else {
288  decorator_e(*part) = pjItr->e();
289  decorator_pt(*part) = pjItr->pt();
290  decorator_eta(*part) = pjItr->pseudorapidity();
291  decorator_phi(*part) = pjItr->phi_std(); //returns phi in [-pi,pi]
292  }
293  found=true;
294  break;
295  } // Found the matching unique ID
296  } // Loop over the jet constituents
297  if (found){
298  for(const auto& constit : constituents) {
299  photon_uniqueIDs.push_back(constit.user_index());
300  } // Loop over the constituents
301  } // Found one of the key leptons in this jet
302  ++pjItr;
303  }
304  if (!found) {
305  if(part->isTau()) {
306  decorator_pt_vis(*part) = 0.;
307  decorator_eta_vis(*part) = 0.;
308  decorator_phi_vis(*part) = 0.;
309  decorator_m_vis(*part) = 0.;
310  }
311  else {
312  decorator_e(*part) = 0;
313  decorator_pt(*part) = 0;
314  decorator_eta(*part) = 0;
315  decorator_phi(*part) = 0;
316  }
317  ATH_MSG_WARNING("Bare particle not found in constituents ");
318  }
319  }
320  // Check if we wanted to decorate photons used for dressing
321  if (!m_decorationName.empty()){
322  //loop over photons, uniquely associate each to nearest bare particle
323  for (const auto* phot : photonsFSRList ) {
324  bool found = std::find(photon_uniqueIDs.begin(), photon_uniqueIDs.end(), HepMC::uniqueID(phot)) != photon_uniqueIDs.end();
325  if (found) {
326  dressDec(*phot) = 1;
327  }
328  } // End of loop over photons
329  } // End of decoration of photons used in dressing
330  } // End of anti-kT dressing
331 
332  return StatusCode::SUCCESS;
333 }
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
DerivationFramework::TruthDressingTool::initialize
StatusCode initialize()
Definition: TruthDressingTool.cxx:42
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:79
get_generator_info.result
result
Definition: get_generator_info.py:21
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
DerivationFramework::TruthDressingTool::addBranches
virtual StatusCode addBranches() const
Definition: TruthDressingTool.cxx:83
SG::AuxTypeRegistry::instance
static AuxTypeRegistry & instance()
Return the singleton registry instance.
Definition: AuxTypeRegistry.cxx:639
SG::ReadHandle< xAOD::TruthParticleContainer >
SG::AuxTypeRegistry::getName
std::string getName(SG::auxid_t auxid) const
Return the name of an aux data item.
Definition: AuxTypeRegistry.cxx:881
SG::ConstAccessor< unsigned int >
MCTruthClassifier.h
DerivationFramework::DecayGraphHelper::constructListOfFinalParticles
bool constructListOfFinalParticles(const xAOD::TruthParticleContainer *allParticles, std::vector< const xAOD::TruthParticle * > &selectedlist, const std::vector< int > &pdgId, bool allowFromHadron=false, bool chargedOnly=false) const
Definition: DecayGraphHelper.h:218
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
DecorKeyHelpers.h
Some common helper functions used by decoration handles.
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
jet::ClusterSequence
fastjet::ClusterSequence ClusterSequence
Definition: ClusterSequence.h:21
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:729
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
xAOD::TruthParticle_v1
Class describing a truth particle in the MC record.
Definition: TruthParticle_v1.h:37
SG::WriteDecorHandle
Handle class for adding a decoration to an object.
Definition: StoreGate/StoreGate/WriteDecorHandle.h:100
WriteDecorHandle.h
Handle class for adding a decoration to an object.
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:116
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TruthDressingTool.h
DerivationFramework::TruthDressingTool::TruthDressingTool
TruthDressingTool(const std::string &t, const std::string &n, const IInterface *p)
Definition: TruthDressingTool.cxx:30
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
isTau
bool isTau(const T &p)
Definition: AtlasPID.h:205
MagicNumbers.h
DerivationFramework::DecayGraphHelper
Definition: DecayGraphHelper.h:26
RegSelToolConfig.alg
alg
Definition: RegSelToolConfig.py:332
checkTriggerxAOD.found
found
Definition: checkTriggerxAOD.py:328
xAOD::JetAlgorithmType::antikt_algorithm
@ antikt_algorithm
Definition: JetContainerInfo.h:33
SG::ReadHandle::ptr
const_pointer_type ptr()
Dereference the pointer.
DerivationFramework::TruthDressingTool::~TruthDressingTool
~TruthDressingTool()
Definition: TruthDressingTool.cxx:38
MCTruthPartClassifier::isPrompt
int isPrompt(const unsigned int classify, bool allow_prompt_tau_decays=true)
Definition: TruthClassifiers.h:180
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
ReadDecorHandle.h
Handle class for reading a decoration on an object.
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
TruthEventContainer.h
xAOD::TruthParticle_v1::FourMom_t
IParticle::FourMom_t FourMom_t
Definition of the 4-momentum type.
Definition: TruthParticle_v1.h:142