HepMC3 event record library
Pythia8ToHepMC3.h
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2021 The HepMC collaboration (see AUTHORS for details)
5 // Part of code was adopted from Pythia8-HepMC interface by Mikhail Kirsanov.
6 #ifndef Pythia8_Pythia8ToHepMC3_H
7 #define Pythia8_Pythia8ToHepMC3_H
8 #ifdef _MSC_VER
9 #pragma message("HepMC3 interface is available in the latest versions of Pythia8, see https://pythia.org. This interface will be removed in the future HepMC3 versions.")
10 #else
11 #warning "HepMC3 interface is available in the latest versions of Pythia8, see https://pythia.org. This interface will be removed in the future HepMC3 versions."
12 #endif
13 #include <vector>
14 #include "Pythia8/Pythia.h"
15 #include "HepMC3/GenVertex.h"
16 #include "HepMC3/GenParticle.h"
17 #include "HepMC3/GenEvent.h"
18 
19 namespace HepMC3 {
20 
21 
23 
24 public:
25 
26  // Constructor and destructor
27  Pythia8ToHepMC3(): m_internal_event_number(0),
28  m_print_inconsistency(true),
29  m_free_parton_warnings(true),
30  m_crash_on_problem(false),
31  m_convert_gluon_to_0(false),
32  m_store_pdf(true),
33  m_store_proc(true),
34  m_store_xsec(true),
35  m_store_weights(true) {}
36 
37  virtual ~Pythia8ToHepMC3() {}
38 
39  // The recommended method to convert Pythia events into HepMC ones
40  bool fill_next_event( Pythia8::Pythia& pythia, GenEvent* evt, int ievnum = -1 )
41  {
42  return fill_next_event( pythia.event, evt, ievnum, &pythia.info, &pythia.settings);
43  }
44 
45  // Alternative method to convert Pythia events into HepMC ones
46 #if defined(PYTHIA_VERSION_INTEGER) && (PYTHIA_VERSION_INTEGER>8299)
47  bool fill_next_event( Pythia8::Event& pyev, GenEvent* evt,
48  int ievnum = -1, const Pythia8::Info* pyinfo = 0,
49  Pythia8::Settings* pyset = 0)
50 #else
51  bool fill_next_event( Pythia8::Event& pyev, GenEvent* evt,
52  int ievnum = -1, Pythia8::Info* pyinfo = 0,
53  Pythia8::Settings* pyset = 0)
54 #endif
55  {
56 
57  // 1. Error if no event passed.
58  if (!evt) {
59  std::cerr << "Pythia8ToHepMC3::fill_next_event error - passed null event."
60  << std::endl;
61  return false;
62  }
63 
64  // Event number counter.
65  if ( ievnum >= 0 ) {
66  evt->set_event_number(ievnum);
67  m_internal_event_number = ievnum;
68  }
69  else {
70  evt->set_event_number(m_internal_event_number);
71  ++m_internal_event_number;
72  }
73 
74  evt->set_units(Units::GEV,Units::MM);
75 
76  // 2. Fill particle information
77  std::vector<GenParticlePtr> hepevt_particles;
78  hepevt_particles.reserve( pyev.size() );
79 
80  for(int i=0; i<pyev.size(); ++i) {
81  hepevt_particles.push_back( std::make_shared<GenParticle>( FourVector( pyev[i].px(), pyev[i].py(),
82  pyev[i].pz(), pyev[i].e() ),
83  pyev[i].id(), pyev[i].statusHepMC() )
84  );
85  hepevt_particles[i]->set_generated_mass( pyev[i].m() );
86  }
87 
88  // 3. Fill vertex information and find beam particles.
89  std::vector<GenVertexPtr> vertex_cache;
90  std::vector<GenParticlePtr> beam_particles;
91  for(int i=0; i<pyev.size(); ++i) {
92 
93  std::vector<int> mothers = pyev[i].motherList();
94 
95  if(mothers.size()) {
96  GenVertexPtr prod_vtx = hepevt_particles[mothers[0]]->end_vertex();
97 
98  if(!prod_vtx) {
99  prod_vtx = std::make_shared<GenVertex>();
100  vertex_cache.push_back(prod_vtx);
101 
102  for(unsigned int j=0; j<mothers.size(); ++j) {
103  prod_vtx->add_particle_in( hepevt_particles[mothers[j]] );
104  }
105  }
106  FourVector prod_pos( pyev[i].xProd(), pyev[i].yProd(),pyev[i].zProd(), pyev[i].tProd() );
107 
108  // Update vertex position if necessary
109  if(!prod_pos.is_zero() && prod_vtx->position().is_zero()) prod_vtx->set_position( prod_pos );
110 
111  prod_vtx->add_particle_out( hepevt_particles[i] );
112  }
113  else beam_particles.push_back(hepevt_particles[i]);
114  }
115 
116  // Reserve memory for the event
117  evt->reserve( hepevt_particles.size(), vertex_cache.size() );
118 
119  // Add particles and vertices in topological order
120  if (beam_particles.size()!=2) {
121  std::cerr << "There are " << beam_particles.size() <<"!=2 particles without mothers"<< std::endl;
122  if ( m_crash_on_problem ) exit(1);
123  }
124  evt->add_tree( beam_particles );
125  //Attributes should be set after adding the particles to event
126  for(int i=0; i<pyev.size(); ++i) {
127  /* TODO: Set polarization */
128  // Colour flow uses index 1 and 2.
129  int colType = pyev[i].colType();
130  if (colType == -1 ||colType == 1 || colType == 2)
131  {
132  int flow1=0, flow2=0;
133  if (colType == 1 || colType == 2) flow1=pyev[i].col();
134  if (colType == -1 || colType == 2) flow2=pyev[i].acol();
135  hepevt_particles[i]->add_attribute("flow1",std::make_shared<IntAttribute>(flow1));
136  hepevt_particles[i]->add_attribute("flow2",std::make_shared<IntAttribute>(flow2));
137  }
138  }
139 
140  // If hadronization switched on then no final coloured particles.
141  bool doHadr = (pyset == 0) ? m_free_parton_warnings : pyset->flag("HadronLevel:all") && pyset->flag("HadronLevel:Hadronize");
142 
143  // 4. Check for particles which come from nowhere, i.e. are without
144  // mothers or daughters. These need to be attached to a vertex, or else
145  // they will never become part of the event.
146  for (int i = 1; i < pyev.size(); ++i) {
147 
148  // Check for particles not added to the event
149  // NOTE: We have to check if this step makes any sense in HepMC event standard
150  if ( !hepevt_particles[i]->in_event() ) {
151  std::cerr << "hanging particle " << i << std::endl;
152  GenVertexPtr prod_vtx = std::make_shared<GenVertex>();
153  prod_vtx->add_particle_out( hepevt_particles[i] );
154  evt->add_vertex(prod_vtx);
155  }
156 
157  // Also check for free partons (= gluons and quarks; not diquarks?).
158  if ( doHadr && m_free_parton_warnings ) {
159  if ( hepevt_particles[i]->pid() == 21 && !hepevt_particles[i]->end_vertex() ) {
160  std::cerr << "gluon without end vertex " << i << std::endl;
161  if ( m_crash_on_problem ) exit(1);
162  }
163  if ( std::abs(hepevt_particles[i]->pid()) <= 6 && !hepevt_particles[i]->end_vertex() ) {
164  std::cerr << "quark without end vertex " << i << std::endl;
165  if ( m_crash_on_problem ) exit(1);
166  }
167  }
168  }
169 
170 
171  // 5. Store PDF, weight, cross section and other event information.
172  // Flavours of incoming partons.
173  if (m_store_pdf && pyinfo != 0) {
174  int id1pdf = pyinfo->id1pdf();
175  int id2pdf = pyinfo->id2pdf();
176  if ( m_convert_gluon_to_0 ) {
177  if (id1pdf == 21) id1pdf = 0;
178  if (id2pdf == 21) id2pdf = 0;
179  }
180 
181  GenPdfInfoPtr pdfinfo = std::make_shared<GenPdfInfo>();
182  pdfinfo->set(id1pdf, id2pdf, pyinfo->x1pdf(), pyinfo->x2pdf(), pyinfo->QFac(), pyinfo->pdf1(), pyinfo->pdf2() );
183  // Store PDF information.
184  evt->set_pdf_info( pdfinfo );
185  }
186 
187  // Store process code, scale, alpha_em, alpha_s.
188  if (m_store_proc && pyinfo != 0) {
189  evt->add_attribute("mpi",std::make_shared<IntAttribute>(pyinfo->nMPI()));
190  evt->add_attribute("signal_process_id",std::make_shared<IntAttribute>( pyinfo->code()));
191  evt->add_attribute("event_scale",std::make_shared<DoubleAttribute>(pyinfo->QRen()));
192  evt->add_attribute("alphaQCD",std::make_shared<DoubleAttribute>(pyinfo->alphaS()));
193  evt->add_attribute("alphaQED",std::make_shared<DoubleAttribute>(pyinfo->alphaEM()));
194  }
195 
196  // Store cross-section information in pb.
197  if (m_store_xsec && pyinfo != 0) {
198  GenCrossSectionPtr xsec = std::make_shared<GenCrossSection>();
199  xsec->set_cross_section( pyinfo->sigmaGen() * 1e9, pyinfo->sigmaErr() * 1e9);
200  evt->set_cross_section(xsec);
201  }
202 
203  // Store event weights.
204  if (m_store_weights && pyinfo != 0) {
205  evt->weights().clear();
206  for (int iweight=0; iweight < pyinfo->nWeights(); ++iweight) {
207  evt->weights().push_back(pyinfo->weight(iweight));
208  }
209  }
210 
211  // Done.
212  return true;
213  }
214 
215  // Read out values for some switches
216  bool print_inconsistency() const {
217  return m_print_inconsistency;
218  }
219  bool free_parton_warnings() const {
220  return m_free_parton_warnings;
221  }
222  bool crash_on_problem() const {
223  return m_crash_on_problem;
224  }
225  bool convert_gluon_to_0() const {
226  return m_convert_gluon_to_0;
227  }
228  bool store_pdf() const {
229  return m_store_pdf;
230  }
231  bool store_proc() const {
232  return m_store_proc;
233  }
234  bool store_xsec() const {
235  return m_store_xsec;
236  }
237  bool store_weights() const {
238  return m_store_weights;
239  }
240 
241  // Set values for some switches
242  void set_print_inconsistency(bool b = true) {
243  m_print_inconsistency = b;
244  }
245  void set_free_parton_warnings(bool b = true) {
246  m_free_parton_warnings = b;
247  }
248  void set_crash_on_problem(bool b = false) {
249  m_crash_on_problem = b;
250  }
251  void set_convert_gluon_to_0(bool b = false) {
252  m_convert_gluon_to_0 = b;
253  }
254  void set_store_pdf(bool b = true) {
255  m_store_pdf = b;
256  }
257  void set_store_proc(bool b = true) {
258  m_store_proc = b;
259  }
260  void set_store_xsec(bool b = true) {
261  m_store_xsec = b;
262  }
263  void set_store_weights(bool b = true) {
264  m_store_weights = b;
265  }
266 
267 private:
268 
269  // Following methods are not implemented for this class
270  virtual bool fill_next_event( GenEvent* ) {
271  return 0;
272  }
273  virtual void write_event( const GenEvent* ) {}
274 
275  // Use of copy constructor is not allowed
276  Pythia8ToHepMC3( const Pythia8ToHepMC3& ) {}
277 
278  // Data members
279  int m_internal_event_number;
280  bool m_print_inconsistency;
281  bool m_free_parton_warnings;
282  bool m_crash_on_problem;
283  bool m_convert_gluon_to_0;
284  bool m_store_pdf;
285  bool m_store_proc;
286  bool m_store_xsec;
287  bool m_store_weights;
288 };
289 
290 } // namespace HepMC3
291 #endif
void set_cross_section(GenCrossSectionPtr cs)
Set cross-section information.
Definition: GenEvent.h:179
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:96
Definition of class GenParticle.
Definition of class GenVertex.
bool is_zero() const
Check if the length of this vertex is zero.
Definition: FourVector.h:193
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Definition: GenEvent.cc:806
Stores event-related information.
Definition: GenEvent.h:41
Generic 4-vector.
Definition: FourVector.h:36
void set_pdf_info(GenPdfInfoPtr pi)
Set PDF information.
Definition: GenEvent.h:172
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:265
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Definition: GenEvent.cc:385
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
Definition: GenEvent.cc:391
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:98
Definition of class GenEvent.
Feature< Feature_type > abs(const Feature< Feature_type > &input)
Obtain the absolute value of a Feature. This works as you&#39;d expect. If foo is a valid Feature...
Definition: Feature.h:323
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:150