ATLAS Offline Software
PrintMC.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // PrintMC - dump details of MC events
6 // Updated 29/10/2009 by Andy Buckley <andy.buckley@cern.ch>
7 // Updated 18.06.2020 by <andrii.verbytskyi@mpp.mpg.de>
8 #include "TruthIO/PrintMC.h"
9 #include "GeneratorObjects/McEventCollection.h"
10 
11 #include "HepPDT/ParticleData.hh"
12 #include "HepPDT/ParticleDataTable.hh"
13 #include "AtlasHepMC/GenEvent.h"
14 
15 
16 inline void drawLine(std::ostream& os) {
17  os << std::string(80, '_') << "\n" << std::endl;
18 }
19 
20 
21 PrintMC::PrintMC(const std::string& name, ISvcLocator* pSvcLocator)
22  : GenBase(name, pSvcLocator)
23 {
24  // Declare JO interface
25  declareProperty("VerboseOutput", m_VerboseOutput=true, "Print out detailed information or not");
26  declareProperty("PrintStyle", m_printsty="Barcode", "Type of info to print out: 'Barcode' or 'Vertex'");
27  declareProperty("VertexInfo", m_vertexinfo=false, "Print out vertex information");
28  declareProperty("TrustHepMC", m_trustHepMC=false, "Use the event numbers in the HepMC record to choose the events to print");
29  declareProperty("FirstEvent", m_firstEvt=1, "Lowest event number to print");
30  declareProperty("LastEvent", m_lastEvt=100000000, "Highest event number to print");
31 }
32 
33 
34 StatusCode PrintMC::initialize() {
35  CHECK(GenBase::initialize());
36  if(!m_trustHepMC) ATH_CHECK(m_evtInfoKey.initialize());
37 
38  // Check settings
39  if (m_printsty != "Barcode" && m_printsty != "Vertex") {
40  ATH_MSG_ERROR("Unknown PrintStyle: " << m_printsty << "! Choose Barcode/Vertex");
41  return StatusCode::FAILURE;
42  }
43 
44  // Check that last > first
45  if (m_lastEvt < m_firstEvt) m_lastEvt = m_firstEvt;
46 
47  return StatusCode::SUCCESS;
48 }
49 
50 
51 
53 StatusCode PrintMC::execute() {
54  // If output already turned off by passing last dumped event, just return
56  if (!m_VerboseOutput) return StatusCode::SUCCESS;
57 
58  // Loop over all events in McEventCollection
59 
60  for (const HepMC::GenEvent* evt : *events_const()) {
61 
62  // Get event number from HepMC
63  uint64_t evtnum = std::max(0,evt->event_number());
64  // Override with evtnum from Athena if enabled and functional
65  if (!m_trustHepMC) {
66  SG::ReadHandle<xAOD::EventInfo> evtInfo(m_evtInfoKey);
67  evtnum = evtInfo->eventNumber();
68  }
69 
70  // Return immediately if event number is outside the print range
71  if (evtnum < m_firstEvt) {
72  return StatusCode::SUCCESS;
73  }
74  if (evtnum > m_lastEvt) {
76  m_VerboseOutput = false;
77  return StatusCode::SUCCESS;
78  }
79 
80  // Printout header
81  ATH_MSG_INFO(">>> PrintMC from execute");
82 
83  // VERTEX format
85  if (m_VerboseOutput && m_printsty == "Vertex") {
86  HepMC::Print::line(std::cout, *evt); // standard HepMC dump of vertex list
87  }
88  // BARCODE format
90  else if (m_VerboseOutput && m_printsty == "Barcode") {
91  drawLine(std::cout);
92 
93  std::cout << "GenEvent: #" << evt->event_number()
94  << " ID=" << HepMC::signal_process_id(evt)
95  << " SignalProcessGenVertex Barcode: "
96  << ( HepMC::signal_process_vertex(evt) ? HepMC::barcode(HepMC::signal_process_vertex(evt)) : 0 ) << "\n";
97  std::cout << " Entries this event: " << evt->vertices_size() << " vertices, "
98  << evt->particles_size() << " particles.\n";
99 #ifdef HEPMC3
100  if (evt->heavy_ion()) {
101  std::cout << " HeavyIon: jatt=" << evt->heavy_ion()->Ncoll_hard
102  << " np=" << evt->heavy_ion()->Npart_proj
103  << " nt=" << evt->heavy_ion()->Npart_targ
104  << " ncoll=" << evt->heavy_ion()->Ncoll
105  << " specn=" << evt->heavy_ion()->spectator_neutrons
106  << " specp=" << evt->heavy_ion()->spectator_protons
107  << " n01=" << evt->heavy_ion()->N_Nwounded_collisions
108  << " n10=" << evt->heavy_ion()->Nwounded_N_collisions
109  << " n11=" << evt->heavy_ion()->Nwounded_Nwounded_collisions
110  << " impact=" << evt->heavy_ion()->impact_parameter
111  << " evplane=" << evt->heavy_ion()->event_plane_angle
112  << " ecc=" << evt->heavy_ion()->eccentricity
113  << " sigmaNNinel=" << evt->heavy_ion()->sigma_inel_NN
114  << std::endl;
115  }
116  else {
117  std::cout << "HeavyIon: EMPTY"
118  << std::endl;
119  }
120 
121 
122  // Weights
123  std::cout << " Weights(" << evt->weights().size() << ")=";
124  for ( auto wgt = evt->weights().begin();
125  wgt != evt->weights().end(); wgt++ ) { std::cout << *wgt << " "; }
126  std::cout << "\n";
127  std::cout << " EventScale " << (evt->attribute<HepMC3::DoubleAttribute>("event_scale")? evt->attribute<HepMC3::DoubleAttribute>("event_scale")->value():0.0)
128  << " [energy] \t alphaQCD=" << (evt->attribute<HepMC3::DoubleAttribute>("alphaQCD")? evt->attribute<HepMC3::DoubleAttribute>("alphaQCD")->value():0.0)
129  << "\t alphaQED=" << (evt->attribute<HepMC3::DoubleAttribute>("alphaQED")? evt->attribute<HepMC3::DoubleAttribute>("alphaQED")->value():0.0) << std::endl;
130 
131  if (evt->pdf_info()) {
132  std::cout << "PdfInfo: id1=" << evt->pdf_info()->parton_id[0]
133  << " id2=" << evt->pdf_info()->parton_id[1]
134  << " x1=" << evt->pdf_info()->x[0]
135  << " x2=" << evt->pdf_info()->x[1]
136  << " q=" << evt->pdf_info()->scale
137  << " xpdf1=" << evt->pdf_info()->pdf_id[0]
138  << " xpdf2=" << evt->pdf_info()->pdf_id[1]
139  << std::endl;
140  }
141  else {
142  std::cout << "PdfInfo: EMPTY"
143  << std::endl;
144  }
145 #else
146  if (evt->heavy_ion()) {
147  std::cout << " HeavyIon: jatt=" << evt->heavy_ion()->Ncoll_hard()
148  << " np=" << evt->heavy_ion()->Npart_proj()
149  << " nt=" << evt->heavy_ion()->Npart_targ()
150  << " ncoll=" << evt->heavy_ion()->Ncoll()
151  << " specn=" << evt->heavy_ion()->spectator_neutrons()
152  << " specp=" << evt->heavy_ion()->spectator_protons()
153  << " n01=" << evt->heavy_ion()->N_Nwounded_collisions()
154  << " n10=" << evt->heavy_ion()->Nwounded_N_collisions()
155  << " n11=" << evt->heavy_ion()->Nwounded_Nwounded_collisions()
156  << " impact=" << evt->heavy_ion()->impact_parameter()
157  << " evplane=" << evt->heavy_ion()->event_plane_angle()
158  << " ecc=" << evt->heavy_ion()->eccentricity()
159  << " sigmaNNinel=" << evt->heavy_ion()->sigma_inel_NN()
160  << std::endl;
161  }
162  else {
163  std::cout << "HeavyIon: EMPTY"
164  << std::endl;
165  }
166 
167 
168  // Weights
169  std::cout << " Weights(" << evt->weights().size() << ")=";
170  for (double w : evt->weights()) {
171  std::cout << w << " ";
172  }
173  std::cout << "\n";
174  std::cout << " EventScale " << evt->event_scale()
175  << " [energy] \t alphaQCD=" << evt->alphaQCD()
176  << "\t alphaQED=" << evt->alphaQED() << std::endl;
177 
178  if (evt->pdf_info()) {
179  std::cout << "PdfInfo: id1=" << evt->pdf_info()->id1()
180  << " id2=" << evt->pdf_info()->id2()
181  << " x1=" << evt->pdf_info()->x1()
182  << " x2=" << evt->pdf_info()->x2()
183  << " q=" << evt->pdf_info()->scalePDF()
184  << " xpdf1=" << evt->pdf_info()->pdf1()
185  << " xpdf2=" << evt->pdf_info()->pdf2()
186  << std::endl;
187  }
188  else {
189  std::cout << "PdfInfo: EMPTY"
190  << std::endl;
191  }
192 #endif
193 
194  // Print a legend to describe the particle info
195  char particle_legend[120];
196  sprintf( particle_legend," %9s %8s %-15s %4s %8s %8s (%9s,%9s,%9s,%9s,%9s)",
197  "Barcode","PDG ID","Name","Stat","ProdVtx","DecayVtx","Px","Py","Pz","E ","m");
198  std::cout << std::endl;
199  std::cout << " GenParticle Legend\n" << particle_legend << "\n";
200  if (m_vertexinfo) {
201  sprintf( particle_legend," %60s (%9s,%9s,%9s,%9s)"," ","Vx","Vy","Vz","Vct ");
202  std::cout << particle_legend << std::endl;
203  }
204  drawLine(std::cout);
205 
206  // Print all particles
207  for (auto p: *evt) {
208  int p_bcode = HepMC::barcode(p);
209  int p_pdg_id = p->pdg_id();
210  HepMC::FourVector mom=p->momentum();
211  int p_stat = p->status();
212  int p_prodvtx = p->production_vertex()?HepMC::barcode(p->production_vertex()):0;
213  int p_endvtx = p->end_vertex()?HepMC::barcode(p->end_vertex()):0;
214  HepMC::FourVector prodvtx=p->production_vertex()?p->production_vertex()->position():HepMC::FourVector(0.0,0.0,0.0,0.0);
215  // Access the PDG table to get the particle name (and mass?)
216  std::string sname;
217  double p_mass = p->generated_mass();
218  const HepPDT::ParticleData* ap = particleData(std::abs(p_pdg_id));
219  if (!ap) {
220  ATH_MSG_DEBUG("PID " << std::abs(p_pdg_id) << " is not in particle data table");
221  } else {
222  const double p_charge = ap->charge() * (p_pdg_id < 0 ? -1 : 1); // assuming that charged leptons are in the PDT...
223  // Build particle name string
224  sname = ap->name();
225  if (p_charge < 0) {
226  const size_t plusidx = sname.rfind("+");
227  if (plusidx != std::string::npos) {
228  sname.replace(plusidx, 1, "-");
229  }
230  }
231  }
232  // PDG table missing or exceptions
233  // quarks & gluons
234  if (p_pdg_id == 21) sname="g";
235  else if (p_pdg_id == 1) sname="d";
236  else if (p_pdg_id == -1) sname="d~";
237  else if (p_pdg_id == 2) sname="u";
238  else if (p_pdg_id == -2) sname="u~";
239  else if (p_pdg_id == 3) sname="s";
240  else if (p_pdg_id == -3) sname="s~";
241  else if (p_pdg_id == 4) sname="c";
242  else if (p_pdg_id == -4) sname="c~";
243  else if (p_pdg_id == 5) sname="b";
244  else if (p_pdg_id == -5) sname="b~";
245  else if (p_pdg_id == 6) sname="t";
246  else if (p_pdg_id == -6) sname="t~";
247  // shower-specific
248  else if (p_pdg_id == 91) sname="cluster";
249  else if (p_pdg_id == 92) sname="string";
250  else if (p_pdg_id == 9922212) sname="remn";
251  else if (p_pdg_id == 2101) sname="ud";
252  else if (p_pdg_id == 2203) sname="uu";
253 
254  // Calc mass if missing
256  if (p_mass == 0 && (p_stat == 2 || (p_stat != 1 && p_pdg_id != 22))) {
257  p_mass = mom.m();
258  }
259 
260  const char* p_name = sname.c_str() ;
261  char particle_entries[120];
262  sprintf(particle_entries, " %9i %8i %-15s %4i %8i %8i (%+9.3g,%+9.3g,%+9.3g,%+9.3g,%9.3g)",
263  p_bcode, p_pdg_id, p_name, p_stat, p_prodvtx, p_endvtx, mom.px(), mom.py(),mom.pz(), mom.e(), p_mass);
264  std::cout << particle_entries << "\n";
265  if (m_vertexinfo) {
266  sprintf(particle_entries," %60s (%+9.3g,%+9.3g,%+9.3g,%+9.3g)"," ",prodvtx.x(), prodvtx.y(),prodvtx.z(), prodvtx.t());
267  std::cout << particle_entries << "\n";
268  }
269  }
270 
271  }
272 
273  }
274  // End of execution for each event
275  return StatusCode::SUCCESS;
276 }
PrintMC::m_lastEvt
uint64_t m_lastEvt
Definition: PrintMC.h:30
GenEvent.h
max
#define max(a, b)
Definition: cfImp.cxx:41
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.
GenBase::events_const
const McEventCollection * events_const() const
Access the current event's McEventCollection (const)
Definition: GenBase.h:96
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
HepMC::signal_process_id
int signal_process_id(const GenEvent &e)
Definition: GenEvent.h:635
PrintMC::m_firstEvt
uint64_t m_firstEvt
Definition: PrintMC.h:29
PrintMC::m_printsty
std::string m_printsty
Definition: PrintMC.h:27
PrintMC::m_VerboseOutput
bool m_VerboseOutput
Definition: PrintMC.h:26
LArG4FSStartPointFilter.evt
evt
Definition: LArG4FSStartPointFilter.py:42
PrintMC::execute
virtual StatusCode execute() override
Definition: PrintMC.cxx:53
HepMC::Print::line
void line(std::ostream &os, const GenEvent &e)
Definition: GenEvent.h:676
python.AtlRunQueryParser.ap
ap
Definition: AtlRunQueryParser.py:826
PrintMC::m_evtInfoKey
SG::ReadHandleKey< xAOD::EventInfo > m_evtInfoKey
Definition: PrintMC.h:20
PrintMC.h
PrintMC::m_vertexinfo
bool m_vertexinfo
Definition: PrintMC.h:28
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
GenBase
Base class for common behaviour of MC truth algorithms.
Definition: GenBase.h:47
McEventCollection.h
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
HepMC::barcode
int barcode(const T *p)
Definition: Barcode.h:16
xAOD::uint64_t
uint64_t
Definition: EventInfo_v1.cxx:123
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
CHECK
#define CHECK(...)
Evaluate an expression and check for errors.
Definition: Control/AthenaKernel/AthenaKernel/errorcheck.h:422
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
ReadFromCoolCompare.os
os
Definition: ReadFromCoolCompare.py:231
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
drawLine
void drawLine(std::ostream &os)
Definition: PrintMC.cxx:16
PrintMC::m_trustHepMC
bool m_trustHepMC
Definition: PrintMC.h:31
PrintMC::initialize
virtual StatusCode initialize() override
Definition: PrintMC.cxx:34
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
GenBase::particleData
const HepPDT::ParticleData * particleData(int pid) const
Access an element in the particle data table.
Definition: GenBase.h:126
GenBase::initialize
virtual StatusCode initialize() override
Definition: GenBase.cxx:17
PrintMC::PrintMC
PrintMC(const std::string &name, ISvcLocator *pSvcLocator)
Definition: PrintMC.cxx:21
HepMC::signal_process_vertex
GenVertex * signal_process_vertex(const GenEvent *e)
Definition: GenEvent.h:625
Generated on Thu Nov 7 2024 21:24:12 for ATLAS Offline Software by   doxygen 1.8.18