ATLAS Offline Software
xAODtoHepMCTool.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 #ifndef XAOD_STANDALONE
8 #endif
10 
12  : asg::AsgTool(name)
13 {
14 }
15 
17 {
18  ATH_MSG_INFO("Initializing xAODtoHepMCTool " << name() << "...");
19  ATH_MSG_INFO("SignalOnly = " << m_signalOnly);
20  return StatusCode::SUCCESS;
21 }
22 
24 {
25  ATH_MSG_INFO("==============================================================");
26  ATH_MSG_INFO("========== xAOD -> HepMC Tool :: Run Summary ==========");
27  ATH_MSG_INFO("==============================================================");
28  if (m_noProdVtx)
29  {
30  ATH_MSG_INFO("Number of events that have a missing production vertex = " << m_noProdVtx);
31  }
32  else
33  {
34  ATH_MSG_INFO("No missing production vertices");
35  }
36  if (m_badBeams)
37  {
38  ATH_MSG_INFO("Number events with improperly defined beamparticles = " << m_badBeams);
39  ATH_MSG_INFO("Inconsistencies with the beam particles storage in the event record!");
40  ATH_MSG_INFO("... check xAOD::TruthEvent record ...");
41  }
42  else
43  {
44  ATH_MSG_INFO("No events with undefined beams.");
45  }
46  ATH_MSG_INFO("==============================================================");
47  ATH_MSG_INFO("=================== End Run Summary ===================");
48  ATH_MSG_INFO("==============================================================");
49  return StatusCode::SUCCESS;
50 }
51 
52 std::vector<HepMC::GenEvent> xAODtoHepMCTool ::getHepMCEvents(const xAOD::TruthEventContainer *xTruthEventContainer, const xAOD::EventInfo *eventInfo) const
53 {
54  ++m_evtCount;
55  bool doPrint = m_evtCount < m_maxCount;
56  // CREATE HEPMC EVENT COLLECTION
57  std::vector<HepMC::GenEvent> mcEventCollection;
58  // Loop over xAOD truth container
59  // Signal event is always first, followed by pileup events
60  ATH_MSG_DEBUG("Start event loop");
61  for (const xAOD::TruthEvent *xAODEvent : *xTruthEventContainer)
62  {
63  if (doPrint)
64  ATH_MSG_DEBUG("XXX Printing xAOD Event");
65  if (doPrint)
66  printxAODEvent(xAODEvent, eventInfo);
67  // Create GenEvent for each xAOD truth event
68  ATH_MSG_DEBUG("Create new GenEvent");
69  HepMC::GenEvent&& hepmcEvent = createHepMCEvent(xAODEvent, eventInfo);
70  #ifdef HEPMC3
71  std::shared_ptr<HepMC3::GenRunInfo> runinfo = std::make_shared<HepMC3::GenRunInfo>(*(hepmcEvent.run_info().get()));
72  #endif
73  //possibly print info, before moving the object
74  if (doPrint){
75  ATH_MSG_DEBUG("XXX Printing HepMC Event");
76  HepMC::Print::line(std::cout, hepmcEvent);
77  }
78  // Insert into McEventCollection
79  mcEventCollection.push_back(std::move(hepmcEvent));
80  #ifdef HEPMC3
81  mcEventCollection[mcEventCollection.size()-1].set_run_info(runinfo);
82  #endif
83  // Quit if signal only
84  if (m_signalOnly)
85  break;
86  }
87  return mcEventCollection;
88 }
89 
90 HepMC::GenEvent xAODtoHepMCTool::createHepMCEvent(const xAOD::TruthEvent *xEvt, const xAOD::EventInfo *eventInfo) const
91 {
92 
94  HepMC::GenEvent genEvt;
95  #ifdef HEPMC3
96  genEvt.set_units(HepMC3::Units::MEV, HepMC3::Units::MM);
97  #endif
98 
99  long long int evtNum = eventInfo->eventNumber();
100  genEvt.set_event_number(evtNum);
101  ATH_MSG_DEBUG("Start createHepMCEvent for event " << evtNum);
102 
103  //Weights
104  #ifndef XAOD_STANDALONE
105  const std::vector<float> weights = xEvt->weights();
106  std::map<std::string, int> weightNameMap;
107  if (AthAnalysisHelper::retrieveMetadata("/Generation/Parameters","HepMCWeightNames", weightNameMap).isFailure()) {
108  ATH_MSG_DEBUG("Couldn't find meta-data for weight names.");
109  }
110  #ifdef HEPMC3
111  std::shared_ptr<HepMC3::GenRunInfo> runinfo = std::make_shared<HepMC3::GenRunInfo>();
112  genEvt.set_run_info(runinfo);
113  std::vector<std::string> wnames;
114  wnames.reserve(weights.size());
115  for (int idx = 0; idx < int(weights.size()); ++idx) {
116  for (const auto& it : weightNameMap) {
117  if (it.second == idx) {
118  wnames.push_back(it.first);
119  break;
120  }
121  }
122  }
123  genEvt.run_info()->set_weight_names(wnames);
124  for ( std::vector<float>::const_iterator wgt = weights.begin(); wgt != weights.end(); ++wgt ) {
125  genEvt.weights().push_back(*wgt);
126  }
127  #else
128  if (weightNameMap.size()) {
129  HepMC::WeightContainer& wc = genEvt.weights();
130  wc.clear();
131  for (int idx = 0; idx < int(weights.size()); ++idx) {
132  for (const auto& it : weightNameMap) {
133  if (it.second == idx) {
134  wc[ it.first ] = weights[idx];
135  break;
136  }
137  }
138  }
139  }
140  else {
141  for ( std::vector<float>::const_iterator wgt = weights.begin(); wgt != weights.end(); ++wgt ) {
142  genEvt.weights().push_back(*wgt);
143  }
144  }
145  #endif
146  #endif
147 
148  // PARTICLES AND VERTICES
149  // Map of existing vertices - needed for the tree linking
150  std::map<const xAOD::TruthVertex *, HepMC::GenVertexPtr> vertexMap;
151 
152  // Loop over all of the particles in the event, call particle builder
153  for (auto tlink : xEvt->truthParticleLinks()) {
154  if (!tlink.isValid()) { continue; }
155  const xAOD::TruthParticle *xPart = *tlink;
156 
157  // sanity check
158  if (xPart == nullptr) {
159  ATH_MSG_WARNING("xAOD TruthParticle is equal to NULL. This should not happen!");
160  continue;
161  }
162 
163  if (!xPart->hasProdVtx() && !xPart->hasDecayVtx()) {
164  ATH_MSG_WARNING("xAOD particle with no vertices, uniqueID = " << HepMC::uniqueID(xPart));
165  continue;
166  }
167 
168  // skip particles which are Geant4 secondaries
169  if (HepMC::is_simulation_particle(xPart)) { continue; }
170 
171  // Create GenParticle
172  // presumably the GenEvent takes ownership of this, but creating a unique_ptr here as that will only happen if there's an associated vertex
173 #ifdef HEPMC3
174  auto hepmcParticle = createHepMCParticle(xPart);
175 #else
176  std::unique_ptr<HepMC::GenParticle> hepmcParticle(createHepMCParticle(xPart));
177 #endif
178 
179  // Get the production and decay vertices
180  if (xPart->hasProdVtx()) {
181  const xAOD::TruthVertex *xAODProdVtx = xPart->prodVtx();
182  // skip production vertices which are Geant4 secondaries
183  if (HepMC::is_simulation_vertex(xAODProdVtx)) { continue; }
184  bool prodVtxSeenBefore(false); // is this new?
185  auto hepmcProdVtx = vertexHelper(xAODProdVtx, vertexMap, prodVtxSeenBefore);
186  // Set the decay/production links
187 #ifdef HEPMC3
188  hepmcProdVtx->add_particle_out(hepmcParticle);
189 #else
190  hepmcProdVtx->add_particle_out(hepmcParticle.get());
191 #endif
192  // Insert into Event
193  if (!prodVtxSeenBefore) {
194  genEvt.add_vertex(hepmcProdVtx);
195  }
196  }
197  else {
198  ATH_MSG_VERBOSE("No production vertex found for particle " << HepMC::uniqueID(xPart));
199  }
200 
201  if (xPart->hasDecayVtx()) {
202  const xAOD::TruthVertex *xAODDecayVtx = xPart->decayVtx();
203  // skip decay vertices which are Geant4 secondaries
204  if (HepMC::is_simulation_vertex(xAODDecayVtx)) {
206 #ifndef HEPMC3
207  if (xPart->hasProdVtx()) { (void)hepmcParticle.release(); }
208 #endif
209  continue;
210  }
211  bool decayVtxSeenBefore(false); // is this new?
212  auto hepmcDecayVtx = vertexHelper(xAODDecayVtx, vertexMap, decayVtxSeenBefore);
213  // Set the decay/production links
214 #ifdef HEPMC3
215  hepmcDecayVtx->add_particle_in(hepmcParticle);
216 #else
217  hepmcDecayVtx->add_particle_in(hepmcParticle.get());
218 #endif
219  // Insert into Event
220  if (!decayVtxSeenBefore) {
221  genEvt.add_vertex(hepmcDecayVtx);
222  }
223  }
224 #ifndef HEPMC3
225  (void)hepmcParticle.release();
226 #endif
227 
228  } // end of particle loop
229 
230  return genEvt;
231 }
232 
233 // Helper to check whether a vertex exists or not using a map;
234 // calls createHepMCVertex if not
236  std::map<const xAOD::TruthVertex *, HepMC::GenVertexPtr> &vertexMap,
237  bool &seenBefore) const
238 {
239 
240  HepMC::GenVertexPtr hepmcVertex;
242  vMapItr = vertexMap.find(xaodVertex);
243  // Vertex seen before?
244  if (vMapItr != vertexMap.end()) {
245  // YES: use the HepMC::Vertex already in the map
246  hepmcVertex = (*vMapItr).second;
247  seenBefore = true;
248  }
249  else {
250  // NO: create a new HepMC::Vertex
251  vertexMap[xaodVertex] = createHepMCVertex(xaodVertex);
252  hepmcVertex = vertexMap[xaodVertex];
253  seenBefore = false;
254  }
255  return hepmcVertex;
256 }
257 
258 // Create the HepMC GenParticle
260 {
261  ATH_MSG_VERBOSE("Creating GenParticle for uniqueID " << HepMC::uniqueID(particle));
262  const HepMC::FourVector fourVec(m_momFac * particle->px(), m_momFac * particle->py(), m_momFac * particle->pz(), m_momFac * particle->e());
263  auto hepmcParticle = HepMC::newGenParticlePtr(fourVec, particle->pdgId(), particle->status());
264  hepmcParticle->set_generated_mass(m_momFac * particle->m());
265  return hepmcParticle;
266 }
267 
268 // Create the HepMC GenVertex
270 {
271  ATH_MSG_VERBOSE("Creating GenVertex for uniqueID " << HepMC::uniqueID(vertex));
272  HepMC::FourVector prod_pos(m_lenFac * vertex->x(), m_lenFac * vertex->y(), m_lenFac * vertex->z(), m_lenFac * vertex->t());
273  auto genVertex = HepMC::newGenVertexPtr(prod_pos);
274  return genVertex;
275 }
276 
277 // Print xAODTruth Event. The printout is particle oriented, unlike the HepMC particle/vertex printout. Regenerated Geant (i.e. ones surving an interaction) and pileup particles (e.g. suppressed pileup barcodes), with barcode > 1000000 are omitted. Thus we output only the original particles created in Geant4.
278 
280 {
281 
282  std::vector<int> uidPars;
283  std::vector<int> uidKids;
284 
285  long long int evtNum = eventInfo->eventNumber();
286 
287  std::cout << "======================================================================================" << std::endl;
288  std::cout << "xAODTruth Event " << evtNum << std::endl;
289  std::cout << " UniqueID PDG Id Status px(GeV) py(GeV) pz(GeV) E(GeV) Parent: Decay" << std::endl;
290  std::cout << " -----------------------------------------------------------------------------------" << std::endl;
291 
292  int nPart = event->nTruthParticles();
293  for (int i = 0; i < nPart; ++i)
294  {
295  const xAOD::TruthParticle *part = event->truthParticle(i);
296  if (part == nullptr) continue;
297  if (HepMC::generations(part) > 0) continue;
298  int id = part->pdgId();
299  if (id != 25) continue;
300  int stat = part->status();
301  float px = part->px() / 1000.;
302  float py = part->py() / 1000.;
303  float pz = part->pz() / 1000.;
304  float e = part->e() / 1000.;
305  uidPars.clear();
306  uidKids.clear();
307 
308  if (part->hasProdVtx())
309  {
310  const xAOD::TruthVertex *pvtx = part->prodVtx();
311  if (pvtx)
312  uidPars.push_back(HepMC::uniqueID(pvtx));
313  }
314 
315  if (part->hasDecayVtx())
316  {
317  const xAOD::TruthVertex *dvtx = part->decayVtx();
318  if (dvtx)
319  uidKids.push_back(HepMC::uniqueID(dvtx));
320  }
321 
322  std::cout << std::setw(10) << HepMC::uniqueID(part) << std::setw(12) << id
323  << std::setw(8) << stat
324  << std::setprecision(2) << std::fixed
325  << std::setw(10) << px << std::setw(10) << py
326  << std::setw(10) << pz << std::setw(10) << e << " ";
327  std::cout << "P: ";
328  for (unsigned int k = 0; k < uidPars.size(); ++k)
329  {
330  std::cout << uidPars[k] << " ";
331  }
332  std::cout << " D: ";
333  for (unsigned int k = 0; k < uidKids.size(); ++k)
334  {
335  std::cout << uidKids[k] << " ";
336  }
337  std::cout << std::endl;
338  }
339  std::cout << "======================================================================================" << std::endl;
340 }
HepMC::GenVertexPtr
HepMC::GenVertex * GenVertexPtr
Definition: GenVertex.h:59
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:79
xAODtoHepMCTool::m_signalOnly
Gaudi::Property< bool > m_signalOnly
Definition: xAODtoHepMCTool.h:49
xAODtoHepMCTool::xAODtoHepMCTool
xAODtoHepMCTool(const std::string &name)
Definition: xAODtoHepMCTool.cxx:11
test_pyathena.px
px
Definition: test_pyathena.py:18
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
xAOD::EventInfo_v1::eventNumber
uint64_t eventNumber() const
The current event's event number.
HepMC::GenParticlePtr
GenParticle * GenParticlePtr
Definition: GenParticle.h:37
skel.it
it
Definition: skel.GENtoEVGEN.py:407
asg
Definition: DataHandleTestTool.h:28
MM
@ MM
Definition: RegSelEnums.h:38
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
HepMC::Print::line
void line(std::ostream &os, const GenEvent &e)
Definition: GenEvent.h:676
DeMoUpdate.runinfo
def runinfo
Definition: DeMoUpdate.py:542
xAODtoHepMCTool::m_badBeams
int m_badBeams
Definition: xAODtoHepMCTool.h:60
xAODtoHepMCTool::printxAODEvent
void printxAODEvent(const xAOD::TruthEvent *event, const xAOD::EventInfo *eventInfo) const
Definition: xAODtoHepMCTool.cxx:279
xAODtoHepMCTool::m_evtCount
std::atomic< int > m_evtCount
Counters.
Definition: xAODtoHepMCTool.h:57
xAODtoHepMCTool::m_noProdVtx
int m_noProdVtx
Definition: xAODtoHepMCTool.h:59
xAODtoHepMCTool.h
xAODtoHepMCTool::m_maxCount
Gaudi::Property< int > m_maxCount
Definition: xAODtoHepMCTool.h:52
xAODtoHepMCTool::createHepMCEvent
HepMC::GenEvent createHepMCEvent(const xAOD::TruthEvent *xEvt, const xAOD::EventInfo *eventInfo) const
Definition: xAODtoHepMCTool.cxx:90
xAOD::TruthParticle_v1::hasDecayVtx
bool hasDecayVtx() const
Check for a decay vertex on this particle.
xAODtoHepMCTool::finalize
StatusCode finalize() override
Definition: xAODtoHepMCTool.cxx:23
event
POOL::TEvent event(POOL::TEvent::kClassAccess)
HepMC::is_simulation_particle
bool is_simulation_particle(const T &p)
Method to establish if a particle (or barcode) was created during the simulation (TODO update to be s...
Definition: MagicNumbers.h:354
HepMC::newGenVertexPtr
GenVertexPtr newGenVertexPtr(const HepMC::FourVector &pos=HepMC::FourVector(0.0, 0.0, 0.0, 0.0), const int i=0)
Definition: GenVertex.h:64
lumiFormat.i
int i
Definition: lumiFormat.py:85
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
xAOD::TruthParticle_v1
Class describing a truth particle in the MC record.
Definition: TruthParticle_v1.h:37
xAODtoHepMCTool::m_lenFac
Gaudi::Property< float > m_lenFac
Definition: xAODtoHepMCTool.h:46
xAOD::TruthParticle_v1::hasProdVtx
bool hasProdVtx() const
Check for a production vertex on this particle.
Definition: TruthParticle_v1.cxx:69
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:116
Amg::pz
@ pz
Definition: GeoPrimitives.h:40
xAOD::TruthEvent_v1
Class describing a signal truth event in the MC record.
Definition: TruthEvent_v1.h:35
HepMC::is_simulation_vertex
bool is_simulation_vertex(const T &v)
Method to establish if the vertex was created during simulation (TODO migrate to be based on status).
Definition: MagicNumbers.h:360
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
beamspotman.stat
stat
Definition: beamspotman.py:264
xAOD::TruthParticle_v1::decayVtx
const TruthVertex_v1 * decayVtx() const
The decay vertex of this particle.
xAOD::TruthEvent_v1::weights
const std::vector< float > & weights() const
Const access to the weights vector.
xAOD::TruthParticle_v1::prodVtx
const TruthVertex_v1 * prodVtx() const
The production vertex of this particle.
Definition: TruthParticle_v1.cxx:75
xAOD::TruthVertex_v1
Class describing a truth vertex in the MC record.
Definition: TruthVertex_v1.h:37
xAODtoHepMCTool::createHepMCParticle
HepMC::GenParticlePtr createHepMCParticle(const xAOD::TruthParticle *) const
Definition: xAODtoHepMCTool.cxx:259
Amg::py
@ py
Definition: GeoPrimitives.h:39
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
MagicNumbers.h
weights
Definition: herwig7_interface.h:38
xAOD::EventInfo_v1
Class describing the basic event information.
Definition: EventInfo_v1.h:43
AthAnalysisHelper.h
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
HepMC::newGenParticlePtr
GenParticlePtr newGenParticlePtr(const HepMC::FourVector &mom=HepMC::FourVector(0.0, 0.0, 0.0, 0.0), int pid=0, int status=0)
Definition: GenParticle.h:39
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
xAODtoHepMCTool::initialize
virtual StatusCode initialize() override
Dummy implementation of the initialisation function.
Definition: xAODtoHepMCTool.cxx:16
xAODtoHepMCTool::vertexHelper
HepMC::GenVertexPtr vertexHelper(const xAOD::TruthVertex *, std::map< const xAOD::TruthVertex *, HepMC::GenVertexPtr > &, bool &) const
Definition: xAODtoHepMCTool.cxx:235
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
xAODtoHepMCTool::createHepMCVertex
HepMC::GenVertexPtr createHepMCVertex(const xAOD::TruthVertex *) const
Definition: xAODtoHepMCTool.cxx:269
xAODtoHepMCTool::getHepMCEvents
std::vector< HepMC::GenEvent > getHepMCEvents(const xAOD::TruthEventContainer *xTruthEventContainer, const xAOD::EventInfo *eventInfo) const override
Definition: xAODtoHepMCTool.cxx:52
xAODtoHepMCTool::m_momFac
Gaudi::Property< float > m_momFac
Input container key (job property)
Definition: xAODtoHepMCTool.h:43
xAOD::TruthEventBase_v1::truthParticleLinks
const TruthParticleLinks_t & truthParticleLinks() const
Get all the truth particles.
HepMC::generations
int generations(const T &p)
Method to return how many interactions a particle has undergone during simulation (TODO migrate to be...
Definition: MagicNumbers.h:357
fitman.k
k
Definition: fitman.py:528
AthAnalysisHelper::retrieveMetadata
static std::string retrieveMetadata(const std::string &folder, const std::string &key, const ServiceHandle< StoreGateSvc > &inputMetaStore)
method that always returns as a string you can use from, e.g, pyROOT with evt = ROOT....
Definition: AthAnalysisHelper.h:254