PYTHIA  8.312
HepMC3.h
1 // HepMC3.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 //
6 // Author: HepMC 3 Collaboration, hepmc-dev@.cern.ch
7 // Based on the HepMC2 interface by Mikhail Kirsanov, Mikhail.Kirsanov@cern.ch.
8 // Header file and function definitions for the Pythia8ToHepMC class,
9 // which converts a PYTHIA event record to the standard HepMC format.
10 
11 #ifndef Pythia8_HepMC3_H
12 #define Pythia8_HepMC3_H
13 
14 #ifdef Pythia8_HepMC2_H
15 #error Cannot include HepMC3.h if HepMC2.h has already been included.
16 #endif
17 
18 #include <vector>
19 #include "Pythia8/Pythia.h"
20 #include "Pythia8/HIInfo.h"
21 #include "HepMC3/GenVertex.h"
22 #include "HepMC3/GenParticle.h"
23 #include "HepMC3/GenEvent.h"
24 #include "HepMC3/WriterAscii.h"
25 #include "HepMC3/WriterAsciiHepMC2.h"
26 #include "HepMC3/GenHeavyIon.h"
27 #include "HepMC3/GenPdfInfo.h"
28 
29 namespace HepMC3 {
30 
31 //==========================================================================
32 
34 
35 public:
36 
37  // Constructor and destructor
38  Pythia8ToHepMC3(): m_internal_event_number(0), m_print_inconsistency(true),
39  m_free_parton_warnings(true), m_crash_on_problem(false),
40  m_convert_gluon_to_0(false), m_store_pdf(true), m_store_proc(true),
41  m_store_xsec(true), m_store_weights(true) {}
42  virtual ~Pythia8ToHepMC3() {}
43 
44  // The recommended method to convert Pythia events into HepMC3 ones.
45  bool fill_next_event( Pythia8::Pythia& pythia, GenEvent* evt,
46  int ievnum = -1 ) { return fill_next_event( pythia.event, evt,
47  ievnum, &pythia.info, &pythia.settings); }
48  bool fill_next_event( Pythia8::Pythia& pythia, GenEvent& evt) {
49  return fill_next_event( pythia, &evt); }
50 
51  // Alternative method to convert Pythia events into HepMC3 ones.
52  bool fill_next_event( Pythia8::Event& pyev, GenEvent&evt, int ievnum = -1,
53  const Pythia8::Info* pyinfo = 0, Pythia8::Settings* pyset = 0) {
54  return fill_next_event(pyev, &evt, ievnum, pyinfo, pyset); }
55  bool fill_next_event( Pythia8::Event& pyev, GenEvent* evt, int ievnum = -1,
56  const Pythia8::Info* pyinfo = 0, Pythia8::Settings* pyset = 0) {
57 
58  // 1. Error if no event passed.
59  if (!evt) {
60  std::cerr << "Pythia8ToHepMC3::fill_next_event error - "
61  << "passed null event." << std::endl;
62  return false;
63  }
64 
65  // Event number counter.
66  if ( ievnum >= 0 ) {
67  evt->set_event_number(ievnum);
68  m_internal_event_number = ievnum;
69  }
70  else {
71  evt->set_event_number(m_internal_event_number);
72  ++m_internal_event_number;
73  }
74 
75  // Set units to be GeV and mm, to agree with Pythia ones.
76  evt->set_units(Units::GEV,Units::MM);
77 
78  // 1a. If there is a HIInfo object fill info from that.
79  if ( pyinfo && pyinfo->hiInfo ) {
80  auto ion = make_shared<HepMC3::GenHeavyIon>();
81  ion->Ncoll_hard = pyinfo->hiInfo->nCollNDTot();
82  ion->Ncoll = pyinfo->hiInfo->nAbsProj() +
83  pyinfo->hiInfo->nDiffProj() +
84  pyinfo->hiInfo->nAbsTarg() +
85  pyinfo->hiInfo->nDiffTarg() -
86  pyinfo->hiInfo->nCollND() -
87  pyinfo->hiInfo->nCollDD();
88  ion->Npart_proj = pyinfo->hiInfo->nAbsProj() +
89  pyinfo->hiInfo->nDiffProj();
90  ion->Npart_targ = pyinfo->hiInfo->nAbsTarg() +
91  pyinfo->hiInfo->nDiffTarg();
92  ion->impact_parameter = pyinfo->hiInfo->b();
93  evt->set_heavy_ion(ion);
94  }
95 
96  // 2. Fill particle information.
97  std::vector<GenParticlePtr> hepevt_particles;
98  hepevt_particles.reserve( pyev.size() );
99  for(int i = 0; i < pyev.size(); ++i) {
100  hepevt_particles.push_back( std::make_shared<GenParticle>(
101  FourVector( pyev[i].px(), pyev[i].py(), pyev[i].pz(), pyev[i].e() ),
102  pyev[i].id(), pyev[i].statusHepMC() ) );
103  hepevt_particles[i]->set_generated_mass( pyev[i].m() );
104  }
105 
106  // 3. Fill vertex information.
107  std::vector<GenVertexPtr> vertex_cache;
108  for (int i = 3; i < pyev.size(); ++i) {
109  std::vector<int> mothers = pyev[i].motherList();
110  if (mothers.size()) {
111  GenVertexPtr prod_vtx = hepevt_particles[mothers[0]]->end_vertex();
112  if (!prod_vtx) {
113  prod_vtx = make_shared<GenVertex>();
114  vertex_cache.push_back(prod_vtx);
115  for (unsigned int j = 0; j < mothers.size(); ++j)
116  prod_vtx->add_particle_in( hepevt_particles[mothers[j]] );
117  }
118  FourVector prod_pos( pyev[i].xProd(), pyev[i].yProd(),pyev[i].zProd(),
119  pyev[i].tProd() );
120 
121  // Update vertex position if necessary.
122  if (!prod_pos.is_zero() && prod_vtx->position().is_zero())
123  prod_vtx->set_position( prod_pos );
124  prod_vtx->add_particle_out( hepevt_particles[i] );
125  }
126  }
127 
128  // Reserve memory for the event.
129  evt->reserve( hepevt_particles.size(), vertex_cache.size() );
130 
131  // Here we assume that the first two particles are the beam particles.
132  vector<GenParticlePtr> beam_particles;
133  beam_particles.push_back(hepevt_particles[1]);
134  beam_particles.push_back(hepevt_particles[2]);
135 
136  // Add particles and vertices in topological order.
137  evt->add_tree( beam_particles );
138  // Attributes should be set after adding the particles to event.
139  for (int i = 0; i < pyev.size(); ++i) {
140  /* TODO: Set polarization */
141  // Colour flow uses index 1 and 2.
142  int colType = pyev[i].colType();
143  if (colType == -1 ||colType == 1 || colType == 2) {
144  int flow1 = 0, flow2 = 0;
145  if (colType == 1 || colType == 2) flow1 = pyev[i].col();
146  if (colType == -1 || colType == 2) flow2 = pyev[i].acol();
147  hepevt_particles[i]->add_attribute("flow1",
148  make_shared<IntAttribute>(flow1));
149  hepevt_particles[i]->add_attribute("flow2",
150  make_shared<IntAttribute>(flow2));
151  }
152  }
153 
154  // If hadronization switched on then no final coloured particles.
155  bool doHadr = (pyset == 0) ? m_free_parton_warnings
156  : pyset->flag("HadronLevel:all") && pyset->flag("HadronLevel:Hadronize");
157 
158  // 4. Check for particles which come from nowhere, i.e. are without
159  // mothers or daughters. These need to be attached to a vertex, or else
160  // they will never become part of the event.
161  for (int i = 1; i < pyev.size(); ++i) {
162 
163  // Check for particles not added to the event.
164  // NOTE: We have to check if this step makes any sense in
165  // the HepMC event standard.
166  if ( !hepevt_particles[i] ) {
167  std::cerr << "hanging particle " << i << std::endl;
168  GenVertexPtr prod_vtx = make_shared<GenVertex>();
169  prod_vtx->add_particle_out( hepevt_particles[i] );
170  evt->add_vertex(prod_vtx);
171  }
172 
173  // Also check for free partons (= gluons and quarks; not diquarks?).
174  if ( doHadr && m_free_parton_warnings ) {
175  if ( hepevt_particles[i]->pid() == 21
176  && !hepevt_particles[i]->end_vertex() ) {
177  std::cerr << "gluon without end vertex " << i << std::endl;
178  if ( m_crash_on_problem ) exit(1);
179  }
180  if ( std::abs(hepevt_particles[i]->pid()) <= 6
181  && !hepevt_particles[i]->end_vertex() ) {
182  std::cerr << "quark without end vertex " << i << std::endl;
183  if ( m_crash_on_problem ) exit(1);
184  }
185  }
186  }
187 
188  // 5. Store PDF, weight, cross section and other event information.
189  // Flavours of incoming partons.
190  if (m_store_pdf && pyinfo != 0) {
191  int id1pdf = pyinfo->id1pdf();
192  int id2pdf = pyinfo->id2pdf();
193  if ( m_convert_gluon_to_0 ) {
194  if (id1pdf == 21) id1pdf = 0;
195  if (id2pdf == 21) id2pdf = 0;
196  }
197 
198  // Store PDF information.
199  GenPdfInfoPtr pdfinfo = make_shared<GenPdfInfo>();
200  pdfinfo->set(id1pdf, id2pdf, pyinfo->x1pdf(), pyinfo->x2pdf(),
201  pyinfo->QFac(), pyinfo->pdf1(), pyinfo->pdf2() );
202  evt->set_pdf_info( pdfinfo );
203  }
204 
205  // Store process code, scale, alpha_em, alpha_s.
206  if (m_store_proc && pyinfo != 0) {
207  evt->add_attribute("signal_process_id",
208  std::make_shared<IntAttribute>( pyinfo->code()));
209  evt->add_attribute("event_scale",
210  std::make_shared<DoubleAttribute>(pyinfo->QRen()));
211  evt->add_attribute("alphaQCD",
212  std::make_shared<DoubleAttribute>(pyinfo->alphaS()));
213  evt->add_attribute("alphaQED",
214  std::make_shared<DoubleAttribute>(pyinfo->alphaEM()));
215  }
216 
217  // Store event weights.
218  if (m_store_weights && pyinfo != 0) {
219  evt->weights().clear();
220  for (int iWeight = 0; iWeight < pyinfo->numberOfWeights(); ++iWeight)
221  evt->weights().push_back(pyinfo->weightValueByIndex(iWeight));
222  }
223 
224  // Store cross-section information in pb.
225  if (m_store_xsec && pyinfo != 0) {
226  // First set atribute to event, such that
227  // GenCrossSection::set_cross_section knows how many weights the
228  // event has and sets the number of cross sections accordingly.
229  GenCrossSectionPtr xsec = make_shared<GenCrossSection>();
230  evt->set_cross_section(xsec);
231  xsec->set_cross_section( pyinfo->sigmaGen() * 1e9,
232  pyinfo->sigmaErr() * 1e9);
233  // If multiweights with possibly different xsec, overwrite central value
234  vector<double> xsecVec = pyinfo->weightContainerPtr->getTotalXsec();
235  if (xsecVec.size() > 0) {
236  for (unsigned int iXsec = 0; iXsec < xsecVec.size(); ++iXsec) {
237  xsec->set_xsec(iXsec, xsecVec[iXsec]*1e9);
238  }
239  }
240  }
241 
242  // Done.
243  return true;
244  }
245 
246  // Read out values for some switches.
247  bool print_inconsistency() const { return m_print_inconsistency; }
248  bool free_parton_warnings() const { return m_free_parton_warnings; }
249  bool crash_on_problem() const { return m_crash_on_problem; }
250  bool convert_gluon_to_0() const { return m_convert_gluon_to_0; }
251  bool store_pdf() const { return m_store_pdf; }
252  bool store_proc() const { return m_store_proc; }
253  bool store_xsec() const { return m_store_xsec; }
254  bool store_weights() const { return m_store_weights; }
255 
256  // Set values for some switches.
257  void set_print_inconsistency(bool b = true) { m_print_inconsistency = b; }
258  void set_free_parton_warnings(bool b = true) { m_free_parton_warnings = b; }
259  void set_crash_on_problem(bool b = false) { m_crash_on_problem = b; }
260  void set_convert_gluon_to_0(bool b = false) { m_convert_gluon_to_0 = b; }
261  void set_store_pdf(bool b = true) { m_store_pdf = b; }
262  void set_store_proc(bool b = true) { m_store_proc = b; }
263  void set_store_xsec(bool b = true) { m_store_xsec = b; }
264  void set_store_weights(bool b = true) { m_store_weights = b; }
265 
266 private:
267 
268  // Following methods are not implemented for this class.
269  virtual bool fill_next_event( GenEvent* ) { return 0; }
270  virtual void write_event( const GenEvent* ) {}
271 
272  // Use of copy constructor is not allowed.
273  Pythia8ToHepMC3( const Pythia8ToHepMC3& ) {}
274 
275  // Data members.
276  int m_internal_event_number;
277  bool m_print_inconsistency, m_free_parton_warnings, m_crash_on_problem,
278  m_convert_gluon_to_0, m_store_pdf, m_store_proc, m_store_xsec,
279  m_store_weights;
280 
281 };
282 
283 //==========================================================================
284 
285 } // end namespace HepMC3
286 
287 namespace Pythia8 {
288 
289 //==========================================================================
290 // This a wrapper around HepMC::Pythia8ToHepMC in the Pythia8
291 // namespace that simplify the most common use cases. It stores the
292 // current GenEvent and output stream internally to avoid cluttering
293 // of user code. This class is also defined in HepMC2.h with the same
294 // signatures, and the user can therefore switch between HepMC version
295 // 2 and 3, by simply changing the include file.
296 class Pythia8ToHepMC : public HepMC3::Pythia8ToHepMC3 {
297 
298 public:
299 
300  // We can either have standard ascii output version 2 or three or
301  // none at all.
302  enum OutputType { none, ascii2, ascii3 };
303 
304  // Typedef for the version 3 specific classes used.
305  typedef HepMC3::GenEvent GenEvent;
306  typedef shared_ptr<GenEvent> EventPtr;
307  typedef HepMC3::Writer Writer;
308  typedef shared_ptr<Writer> WriterPtr;
309 
310  // The empty constructor does not creat an aoutput stream.
311  Pythia8ToHepMC() : runinfo(make_shared<HepMC3::GenRunInfo>()) {}
312 
313  // Construct an object with an internal output stream.
314  Pythia8ToHepMC(string filename, OutputType ft = ascii3)
315  : runinfo(make_shared<HepMC3::GenRunInfo>()) {
316  setNewFile(filename, ft);
317  }
318 
319  // Open a new external output stream.
320  bool setNewFile(string filename, OutputType ft = ascii3) {
321  switch ( ft ) {
322  case ascii3:
323  writerPtr = make_shared<HepMC3::WriterAscii>(filename);
324  break;
325  case ascii2:
326  writerPtr = make_shared<HepMC3::WriterAsciiHepMC2>(filename);
327  break;
328  case none:
329  break;
330  }
331  return writerPtr != nullptr;
332  }
333 
334  // Create a new GenEvent object and fill it with information from
335  // the given Pythia object.
336  bool fillNextEvent(Pythia & pythia) {
337  geneve = make_shared<HepMC3::GenEvent>(runinfo);
338  if (runinfo->weight_names().size() == 0)
339  setWeightNames(pythia.info.weightNameVector());
340  return fill_next_event(pythia, *geneve);
341  }
342 
343  // Write out the current GenEvent to the internal stream.
344  void writeEvent() {
345  writerPtr->write_event(*geneve);
346  }
347 
348  // Create a new GenEvent object and fill it with information from
349  // the given Pythia object and write it out directly to the
350  // internal stream.
351  bool writeNextEvent(Pythia & pythia) {
352  if ( !fillNextEvent(pythia) ) return false;
353  writeEvent();
354  return !writerPtr->failed();
355  }
356 
357  // Get a reference to the current GenEvent.
358  GenEvent & event() {
359  return *geneve;
360  }
361 
362  // Get a pointer to the current GenEvent.
363  EventPtr getEventPtr() {
364  return geneve;
365  }
366 
367  // Get a reference to the internal stream.
368  Writer & output() {
369  return *writerPtr;
370  }
371 
372  // Get a pointer to the internal stream.
373  WriterPtr outputPtr() {
374  return writerPtr;
375  }
376 
377  // Set cross section information in the current GenEvent.
378  void setXSec(double xsec, double xsecerr) {
379  auto xsecptr = geneve->cross_section();
380  if ( !xsecptr ) {
381  xsecptr = make_shared<HepMC3::GenCrossSection>();
382  geneve->set_cross_section(xsecptr);
383  }
384  xsecptr->set_cross_section(xsec, xsecerr);
385  }
386 
387  // Update all weights in the current GenEvent.
388  void setWeights(const vector<double> & wv) {
389  geneve->weights() = wv;
390  }
391 
392  // Set all weight names in the current run.
393  void setWeightNames(const vector<string> &wnv) {
394  runinfo->set_weight_names(wnv);
395  }
396 
397  // Update the PDF information in the current GenEvent
398  void setPdfInfo(int id1, int id2, double x1, double x2,
399  double scale, double xf1, double xf2,
400  int pdf1 = 0, int pdf2 = 0) {
401  auto pdf = make_shared<HepMC3::GenPdfInfo>();
402  pdf->set(id1, id2, x1, x2, scale, xf1, xf2, pdf1, pdf2);
403  geneve->set_pdf_info(pdf);
404  }
405 
406  // Add an additional attribute derived from HepMC3::Attribute
407  // to the current event.
408  template<class T>
409  void addAtribute(const string& name, T& attribute) {
410  shared_ptr<HepMC3::Attribute> att = make_shared<T>(attribute);
411  geneve->add_attribute(name, att);
412  }
413 
414  // Add an attribute of double type.
415  template<class T=double>
416  void addAttribute(const string& name, double& attribute) {
417  auto dAtt = HepMC3::DoubleAttribute(attribute);
418  shared_ptr<HepMC3::Attribute> att =
419  make_shared<HepMC3::DoubleAttribute>(dAtt);
420  geneve->add_attribute(name, att);
421  }
422 
423  // Remove an attribute from the current event.
424  void removeAttribute(const string& name) {
425  geneve->remove_attribute(name);
426  }
427 
428 private:
429 
430  // The current GenEvent
431  EventPtr geneve = nullptr;
432 
433  // The output stream.
434  WriterPtr writerPtr = nullptr;
435 
436  // The current run info.
437  shared_ptr<HepMC3::GenRunInfo> runinfo;
438 
439 };
440 
441 }
442 
443 #endif // end Pythia8_HepMC3_H
void addAttribute(const string &name, double &attribute)
Add an attribute of double type.
Definition: HepMC3.h:416
Definition: HepMC3.h:29
void setWeightNames(const vector< string > &wnv)
Set all weight names in the current run.
Definition: HepMC3.h:393
GenEvent & event()
Get a reference to the current GenEvent.
Definition: HepMC3.h:358
Definition: Info.h:45
The Event class holds all info on the generated event.
Definition: Event.h:453
void setPdfInfo(int id1, int id2, double x1, double x2, double scale, double xf1, double xf2, int pdf1=0, int pdf2=0)
Update the PDF information in the current GenEvent.
Definition: HepMC3.h:398
bool fill_next_event(Pythia8::Event &pyev, GenEvent *evt, int ievnum=-1, const Pythia8::Info *pyinfo=0, Pythia8::Settings *pyset=0)
Definition: HepMC3.h:55
bool setNewFile(string filename, OutputType ft=ascii3)
Open a new external output stream.
Definition: HepMC3.h:320
Settings settings
Settings: databases of flags/modes/parms/words to control run.
Definition: Pythia.h:332
void removeAttribute(const string &name)
Remove an attribute from the current event.
Definition: HepMC3.h:424
bool writeNextEvent(Pythia &pythia)
Definition: HepMC3.h:351
Definition: HepMC3.h:33
Event event
The event record for the complete event history.
Definition: Pythia.h:323
void setXSec(double xsec, double xsecerr)
Set cross section information in the current GenEvent.
Definition: HepMC3.h:378
Pythia8ToHepMC3()
Constructor and destructor.
Definition: HepMC3.h:38
void writeEvent()
Write out the current GenEvent to the internal stream.
Definition: HepMC3.h:344
Pythia8ToHepMC(string filename, OutputType ft=ascii3)
Construct an object with an internal output stream.
Definition: HepMC3.h:314
EventPtr getEventPtr()
Get a pointer to the current GenEvent.
Definition: HepMC3.h:363
void set_print_inconsistency(bool b=true)
Set values for some switches.
Definition: HepMC3.h:257
Pythia8ToHepMC()
The empty constructor does not creat an aoutput stream.
Definition: HepMC3.h:311
void addAtribute(const string &name, T &attribute)
Definition: HepMC3.h:409
int size() const
Event record size.
Definition: Event.h:504
HepMC3::GenEvent GenEvent
Typedef for the version 3 specific classes used.
Definition: HepMC3.h:305
WriterPtr outputPtr()
Get a pointer to the internal stream.
Definition: HepMC3.h:373
void setWeights(const vector< double > &wv)
Update all weights in the current GenEvent.
Definition: HepMC3.h:388
bool fill_next_event(Pythia8::Pythia &pythia, GenEvent *evt, int ievnum=-1)
The recommended method to convert Pythia events into HepMC3 ones.
Definition: HepMC3.h:45
bool fill_next_event(Pythia8::Event &pyev, GenEvent &evt, int ievnum=-1, const Pythia8::Info *pyinfo=0, Pythia8::Settings *pyset=0)
Alternative method to convert Pythia events into HepMC3 ones.
Definition: HepMC3.h:52
OutputType
We can either have standard ascii output or none at all.
Definition: HepMC2.h:413
bool print_inconsistency() const
Read out values for some switches.
Definition: HepMC3.h:247
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
The Pythia class contains the top-level routines to generate an event.
Definition: Pythia.h:71
bool fillNextEvent(Pythia &pythia)
Definition: HepMC3.h:336
const Info & info
Public information and statistic on the generation.
Definition: Pythia.h:326
Writer & output()
Get a reference to the internal stream.
Definition: HepMC3.h:368
Definition: Settings.h:195