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