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
19namespace HepMC3 {
20
21
23
24public:
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
267private:
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
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
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Generic 4-vector.
Definition: FourVector.h:36
bool is_zero() const
Check if the length of this vertex is zero.
Definition: FourVector.h:193
Stores event-related information.
Definition: GenEvent.h:41
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:265
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:96
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
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:150
void set_cross_section(GenCrossSectionPtr cs)
Set cross-section information.
Definition: GenEvent.h:179
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Definition: GenEvent.cc:806
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:98
void set_pdf_info(GenPdfInfoPtr pi)
Set PDF information.
Definition: GenEvent.h:172
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Definition: GenEvent.cc:385
HepMC3 main namespace.