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