ATLAS Offline Software
TruthDecayCollectionMaker.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // TruthDecayCollectionMaker.cxx
7 // Removes all truth particles/vertices which do not pass a user-defined cut
8 // Based on TruthCollectionMaker (but simpler)
9 
11 #include "StoreGate/ReadHandle.h"
12 #include "StoreGate/WriteHandle.h"
19 
21 
22 // STL includes
23 #include <vector>
24 #include <string>
25 // For a find in the vector
26 #include <algorithm>
27 
28 // Constructor
30  const std::string& n,
31  const IInterface* p)
32  : AthAlgTool(t,n,p)
33 {
34  declareInterface<DerivationFramework::IAugmentationTool>(this);
35 }
36 
37 // Destructor
39 }
40 
41 // Athena initialize
43 {
44  ATH_MSG_VERBOSE("initialize() ...");
45 
46  // Input truth particles
47  ATH_CHECK( m_particlesKey.initialize() );
48  ATH_MSG_INFO("Using " << m_particlesKey.key() << " as the input truth container key");
49 
50  // Output particle/vertex containers
51  if (m_collectionName.empty()) {
52  ATH_MSG_FATAL("No base name provided for the new truth particle/vertex containers");
53  return StatusCode::FAILURE;
54  } else {ATH_MSG_INFO("Base name for new truth particle/vertex containers: " << m_collectionName );}
55  m_outputParticlesKey = m_collectionName + "Particles";
56  ATH_CHECK(m_outputParticlesKey.initialize());
57  ATH_MSG_INFO("New truth particles container key: " << m_outputParticlesKey.key() );
58  m_outputVerticesKey = m_collectionName + "Vertices";
59  ATH_CHECK(m_outputVerticesKey.initialize());
60  ATH_MSG_INFO("New truth vertices container key: " << m_outputVerticesKey.key() );
61 
62  if (m_pdgIdsToKeep.empty() && !m_keepBHadrons && !m_keepCHadrons && !m_keepBSM) {
63  ATH_MSG_FATAL("No PDG IDs provided, not keeping b- or c-hadrons or BSM particles -- what do you want?");
64  return StatusCode::FAILURE;
65  }
66 
67  // Decorators
68  m_originDecoratorKey = m_outputParticlesKey.key()+".classifierParticleOrigin";
69  ATH_CHECK(m_originDecoratorKey.initialize());
70  m_typeDecoratorKey = m_outputParticlesKey.key()+".classifierParticleType";
71  ATH_CHECK(m_typeDecoratorKey.initialize());
72  m_outcomeDecoratorKey = m_outputParticlesKey.key()+".classifierParticleOutCome";
73  ATH_CHECK(m_outcomeDecoratorKey.initialize());
74  m_classificationDecoratorKey = m_outputParticlesKey.key()+".Classification";
75  ATH_CHECK(m_classificationDecoratorKey.initialize());
76  m_motherIDDecoratorKey = m_outputParticlesKey.key()+".motherID";
77  ATH_CHECK(m_motherIDDecoratorKey.initialize());
78  m_daughterIDDecoratorKey = m_outputParticlesKey.key()+".daughterID";
79  ATH_CHECK(m_daughterIDDecoratorKey.initialize());
80 
81  return StatusCode::SUCCESS;
82 }
83 
84 
85 // Selection and collection creation
87 {
88  // Event context for AthenaMT
89  const EventContext& ctx = Gaudi::Hive::currentContext();
90 
91  // Retrieve truth collections
92  SG::ReadHandle<xAOD::TruthParticleContainer> truthParticles(m_particlesKey,ctx);
93  if (!truthParticles.isValid()) {
94  ATH_MSG_ERROR("Couldn't retrieve TruthParticle collection with name " << m_particlesKey);
95  return StatusCode::FAILURE;
96  }
97 
98  // Create the new particle containers and WriteHandles
99  SG::WriteHandle<xAOD::TruthParticleContainer> newParticlesWriteHandle(m_outputParticlesKey, ctx);
100  ATH_CHECK(newParticlesWriteHandle.record(std::make_unique<xAOD::TruthParticleContainer>(),
101  std::make_unique<xAOD::TruthParticleAuxContainer>()));
102  ATH_MSG_DEBUG( "Recorded new TruthParticleContainer with key: " << (m_outputParticlesKey.key()));
103  // Create the new vertex containers and WriteHandles
104  SG::WriteHandle<xAOD::TruthVertexContainer> newVerticesWriteHandle(m_outputVerticesKey, ctx);
105  ATH_CHECK(newVerticesWriteHandle.record(std::make_unique<xAOD::TruthVertexContainer>(),
106  std::make_unique<xAOD::TruthVertexAuxContainer>()));
107  ATH_MSG_DEBUG( "Recorded new TruthVertexContainer with key: " << (m_outputVerticesKey.key()));
108 
109  // List of unique IDs for particles in our collection already. Because of the way we recurse,
110  // adding more particles as we go, there should be no need to add (or help from adding) the
111  // unique IDs of particles that we are *not* going to keep
112  std::vector<int> seen_particles;
113  // Go through that list of particles!
114  for (const auto * part : *truthParticles){
115  // If this passes my cuts, keep it
116  if (id_ok(*part)){
117  addTruthParticle( ctx, *part, newParticlesWriteHandle.ptr(), newVerticesWriteHandle.ptr(), seen_particles , m_generations );
118  }
119  } // Loop over the initial truth particle collection
120  return StatusCode::SUCCESS;
121 }
122 
124  const xAOD::TruthParticle& old_part,
125  xAOD::TruthParticleContainer* part_cont,
126  xAOD::TruthVertexContainer* vert_cont,
127  std::vector<int>& seen_particles,
128  const int generations) const {
129  // See if we've seen it - note, could also do this with a unary function on the container itself
130  if (std::find(seen_particles.begin(),seen_particles.end(),HepMC::uniqueID(&old_part))!=seen_particles.end()){
131  for (size_t p=0;p<part_cont->size();++p){
132  // Was it a hit?
133  const xAOD::TruthParticle *theParticle = (*part_cont)[p];
134  if (HepMC::is_same_particle(theParticle,&old_part)) return p;
135  } // Look through the old container
136  } // Found it in the old container
137  // Now we have seen it
138  seen_particles.push_back(HepMC::uniqueID(&old_part));
139  // Set up decorators
140  SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > originDecorator(m_originDecoratorKey, ctx);
141  SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > typeDecorator(m_typeDecoratorKey, ctx);
142  SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > outcomeDecorator(m_outcomeDecoratorKey, ctx);
143  SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > classificationDecorator(m_classificationDecoratorKey, ctx);
144  SG::WriteDecorHandle< xAOD::TruthParticleContainer, int > motherIDDecorator(m_motherIDDecoratorKey, ctx);
145  SG::WriteDecorHandle< xAOD::TruthParticleContainer, int > daughterIDDecorator(m_daughterIDDecoratorKey, ctx);
146  // Make a truth particle and add it to the container
147  xAOD::TruthParticle* xTruthParticle = new xAOD::TruthParticle();
148  part_cont->push_back( xTruthParticle );
149  // Fill with numerical content
150  xTruthParticle->setPdgId(old_part.pdgId());
151  xTruthParticle->setBarcode(HepMC::barcode(&old_part)); // FIXME barcode-based
152  xTruthParticle->setStatus(old_part.status());
153  xTruthParticle->setM(old_part.m());
154  xTruthParticle->setPx(old_part.px());
155  xTruthParticle->setPy(old_part.py());
156  xTruthParticle->setPz(old_part.pz());
157  xTruthParticle->setE(old_part.e());
158  // Copy over the polarization information if it's there
159  if (old_part.polarization().valid()){
162  }
163  // Make a link to this particle
164  int my_index = part_cont->size()-1;
165  ElementLink<xAOD::TruthParticleContainer> eltp(*part_cont, my_index);
166  // Decay vertex information
167  if (old_part.hasDecayVtx()) {
168  int vert_index = addTruthVertex( ctx, *old_part.decayVtx(), part_cont, vert_cont, seen_particles, generations);
169  ElementLink<xAOD::TruthVertexContainer> eltv( *vert_cont, vert_index );
170  xTruthParticle->setDecayVtxLink( eltv );
171  (*vert_cont)[vert_index]->addIncomingParticleLink( eltp );
172  }
173  // Copy over the decorations if they are available
174  static const SG::ConstAccessor<unsigned int> classifierParticleTypeAcc("classifierParticleType");
175  static const SG::ConstAccessor<unsigned int> classifierParticleOriginAcc("classifierParticleOrigin");
176  static const SG::ConstAccessor<unsigned int> classifierParticleOutComeAcc("classifierParticleOutCome");
177  static const SG::ConstAccessor<unsigned int> ClassificationAcc("Classification");
178 
179  typeDecorator(*xTruthParticle) = classifierParticleTypeAcc.withDefault (old_part, 0);
180  originDecorator(*xTruthParticle) = classifierParticleOriginAcc.withDefault (old_part, 0);
181  outcomeDecorator(*xTruthParticle) = classifierParticleOutComeAcc.withDefault (old_part, 0);
182  classificationDecorator(*xTruthParticle) = ClassificationAcc.withDefault (old_part, 0);
183 
184  // Return a link to this particle
185  return my_index;
186 }
187 
189  const xAOD::TruthVertex& old_vert,
190  xAOD::TruthParticleContainer* part_cont,
191  xAOD::TruthVertexContainer* vert_cont,
192  std::vector<int>& seen_particles,
193  const int generations) const {
194  // Make a new vertex and add it to the container
195  xAOD::TruthVertex* xTruthVertex = new xAOD::TruthVertex();
196  vert_cont->push_back( xTruthVertex );
197  // Get a link to this vertex -- will be used to set production vertices on all the next particles
198  int my_index = vert_cont->size()-1;
199  ElementLink<xAOD::TruthVertexContainer> eltv(*vert_cont, my_index);
200  // Set properties
201  xTruthVertex->setId(HepMC::status(old_vert));
202  xTruthVertex->setBarcode(HepMC::barcode(&old_vert)); // FIXME barcode-based
203  xTruthVertex->setX(old_vert.x());
204  xTruthVertex->setY(old_vert.y());
205  xTruthVertex->setZ(old_vert.z());
206  xTruthVertex->setT(old_vert.t());
207  // If we are done, then stop here
208  if (generations==0) return my_index;
209  // Add all the outgoing particles
210  for (size_t n=0;n<old_vert.nOutgoingParticles();++n){
211  if (!old_vert.outgoingParticle(n)) continue; // Just in case we removed some truth particles, e.g. G4 decays
212  if (m_rejectHadronChildren && old_vert.outgoingParticle(n)->isHadron()) { // Option to skip hadrons outright{
213  continue;
214  }
215  // Continue on the next generation; note that we only decrement the generation if this particle doesn't also pass our cuts
216  int part_index = addTruthParticle( ctx, *old_vert.outgoingParticle(n), part_cont, vert_cont, seen_particles,
217  generations-1+(id_ok(*old_vert.outgoingParticle(n))?1:0) );
218  ElementLink<xAOD::TruthParticleContainer> eltp( *part_cont, part_index);
219  xTruthVertex->addOutgoingParticleLink( eltp );
220  (*part_cont)[part_index]->setProdVtxLink( eltv );
221  }
222  // Return a link to this vertex
223  return my_index;
224 }
225 
227 {
228  // Check list of PDG IDs to keep
229  for (int id : m_pdgIdsToKeep){
230  if (part.absPdgId()==id){
231  return true;
232  } // Found a particle of interest!
233  } // Loop over the PDG IDs we want to keep
234  // Also check functions for B/C/BSM
235  return (m_keepBHadrons && part.isBottomHadron()) ||
236 
237  (m_keepCHadrons && part.isCharmHadron()) ||
238 
239  (m_keepBSM && part.isBSM());
240 }
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
xAOD::TruthVertex_v1::nOutgoingParticles
size_t nOutgoingParticles() const
Get the number of outgoing particles.
xAOD::TruthParticle_v1::setStatus
void setStatus(int value)
Set status code.
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
xAOD::TruthParticle_v1::pz
float pz() const
The z component of the particle's momentum.
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
xAOD::TruthParticle_v1::setE
void setE(float value)
Set the energy of the particle.
Definition: TruthParticle_v1.cxx:235
xAOD::TruthParticle_v1::setBarcode
void setBarcode(int value)
Set barcode.
TruthVertexContainer.h
TruthParticleContainer.h
DerivationFramework::TruthDecayCollectionMaker::TruthDecayCollectionMaker
TruthDecayCollectionMaker(const std::string &t, const std::string &n, const IInterface *p)
Definition: TruthDecayCollectionMaker.cxx:29
xAOD::TruthParticle_v1::setPolarizationParameter
bool setPolarizationParameter(float value, PolParam parameter)
Set method for polarization parameter values.
Definition: TruthParticle_v1.cxx:351
xAOD::TruthParticle_v1::Polarization::valid
bool valid() const
Check if the stored values are valid.
Definition: TruthParticle_v1.h:375
xAOD::TruthParticle_v1::polarization
Polarization polarization() const
Retrieve a full Polarization with a single call.
Definition: TruthParticle_v1.cxx:383
xAOD::TruthParticle_v1::px
float px() const
The x component of the particle's momentum.
DerivationFramework::TruthDecayCollectionMaker::initialize
StatusCode initialize()
Definition: TruthDecayCollectionMaker.cxx:42
xAOD::TruthVertex_v1::setT
void setT(float value)
Set the vertex time.
xAOD::TruthParticle_v1::setPx
void setPx(float value)
Set the x component of the particle's momentum.
xAOD::TruthParticle_v1::py
float py() const
The y component of the particle's momentum.
SG::ConstAccessor< unsigned int >
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
xAOD::TruthVertex_v1::addOutgoingParticleLink
void addOutgoingParticleLink(const TPLink_t &link)
Add one outgoing particle.
Definition: TruthVertex_v1.cxx:137
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
xAOD::TruthVertex_v1::setY
void setY(float value)
Set the y displacement of the vertex.
HepMC::is_same_particle
bool is_same_particle(const T1 &p1, const T2 &p2)
Method to establish if two particles in the GenEvent actually represent the same particle.
Definition: MagicNumbers.h:354
xAOD::TruthVertex_v1::y
float y() const
Vertex y displacement.
HepMC::BarcodeBased::generations
int generations(const T &p)
Method to return how many interactions a particle has undergone during simulation (only to be used in...
Definition: MagicNumbers.h:199
xAOD::TruthVertex_v1::t
float t() const
Vertex time.
TruthParticleAuxContainer.h
WriteHandle.h
Handle class for recording to StoreGate.
xAOD::TruthParticle_v1::hasDecayVtx
bool hasDecayVtx() const
Check for a decay vertex on this particle.
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
xAOD::TruthParticle_v1::setM
void setM(float value)
Also store the mass.
Definition: TruthParticle_v1.cxx:241
TruthDecayCollectionMaker.h
DerivationFramework::TruthDecayCollectionMaker::addTruthVertex
int addTruthVertex(const EventContext &, const xAOD::TruthVertex &old_vert, xAOD::TruthParticleContainer *part_cont, xAOD::TruthVertexContainer *vert_cont, std::vector< int > &seen_particles, const int generations=-1) const
Definition: TruthDecayCollectionMaker.cxx:188
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
xAOD::TruthParticle_v1::e
virtual double e() const override final
The total energy of the particle.
DerivationFramework::TruthDecayCollectionMaker::addTruthParticle
int addTruthParticle(const EventContext &ctx, const xAOD::TruthParticle &old_part, xAOD::TruthParticleContainer *part_cont, xAOD::TruthVertexContainer *vert_cont, std::vector< int > &seen_particles, const int generations=-1) const
Definition: TruthDecayCollectionMaker.cxx:123
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
SG::WriteHandle::ptr
pointer_type ptr()
Dereference the pointer.
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
HepMC::barcode
int barcode(const T *p)
Definition: Barcode.h:16
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.
xAOD::TruthParticle
TruthParticle_v1 TruthParticle
Typedef to implementation.
Definition: Event/xAOD/xAODTruth/xAODTruth/TruthParticle.h:15
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:109
xAOD::TruthParticle_v1::setPy
void setPy(float value)
Set the y component of the particle's momentum.
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
xAOD::TruthVertex
TruthVertex_v1 TruthVertex
Typedef to implementation.
Definition: TruthVertex.h:15
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
TruthVertexAuxContainer.h
xAOD::TruthParticle_v1::decayVtx
const TruthVertex_v1 * decayVtx() const
The decay vertex of this particle.
xAOD::TruthVertex_v1
Class describing a truth vertex in the MC record.
Definition: TruthVertex_v1.h:37
xAOD::TruthParticle_v1::setPdgId
void setPdgId(int pid)
Set PDG ID code.
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
DerivationFramework::TruthDecayCollectionMaker::id_ok
bool id_ok(const xAOD::TruthParticle &part) const
Definition: TruthDecayCollectionMaker.cxx:226
xAOD::TruthVertex_v1::setId
void setId(int value)
Obsolete function Set vertex ID code HepMC2 id == HepMC3 status, i.e.
xAOD::TruthParticle_v1::polarizationParameter
bool polarizationParameter(float &value, PolParam parameter) const
Accessor for polarization parameters.
Definition: TruthParticle_v1.cxx:328
xAOD::TruthVertex_v1::x
float x() const
Vertex x displacement.
xAOD::TruthVertex_v1::setZ
void setZ(float value)
Set the vertex's longitudinal distance from the origin.
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:76
xAOD::TruthParticle_v1::status
int status() const
Status code.
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
xAOD::TruthVertex_v1::setBarcode
void setBarcode(int value)
Set barcode.
xAOD::TruthParticle_v1::polarizationTheta
@ polarizationTheta
Polarization in ( )
Definition: TruthParticle_v1.h:319
xAOD::TruthParticle_v1::polarizationPhi
@ polarizationPhi
Polarization in ( )
Definition: TruthParticle_v1.h:318
xAOD::TruthVertex_v1::z
float z() const
Vertex longitudinal distance along the beam line form the origin.
xAOD::TruthParticle_v1::setPz
void setPz(float value)
Set the z component of the particle's momentum.
DerivationFramework::TruthDecayCollectionMaker::~TruthDecayCollectionMaker
~TruthDecayCollectionMaker()
Definition: TruthDecayCollectionMaker.cxx:38
xAOD::TruthParticle_v1::setDecayVtxLink
void setDecayVtxLink(const ElementLink< TruthVertexContainer > &link)
Set the decay vertex of the particle.
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
ReadHandle.h
Handle class for reading from StoreGate.
AthAlgTool
Definition: AthAlgTool.h:26
HepMC::status
int status(const T &p)
Definition: MagicNumbers.h:125
xAOD::TruthParticle_v1::pdgId
int pdgId() const
PDG ID code.
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
xAOD::TruthVertex_v1::outgoingParticle
const TruthParticle_v1 * outgoingParticle(size_t index) const
Get one of the outgoing particles.
Definition: TruthVertex_v1.cxx:119
xAOD::TruthParticle_v1::m
virtual double m() const override final
The mass of the particle.
HepMCHelpers.h
xAOD::TruthParticle_v1::isHadron
bool isHadron() const
Whether the particle is a hadron.
xAOD::TruthVertex_v1::setX
void setX(float value)
Set the x displacement of the vertex.
SG::ConstAccessor::withDefault
const_reference_type withDefault(const ELT &e, const T &deflt) const
Fetch the variable for one element, as a const reference, with a default.
DerivationFramework::TruthDecayCollectionMaker::addBranches
virtual StatusCode addBranches() const
Pass the thinning service
Definition: TruthDecayCollectionMaker.cxx:86