6#ifndef Pythia8_Pythia8ToHepMC3_H
7#define Pythia8_Pythia8ToHepMC3_H
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.")
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."
14#include "Pythia8/Pythia.h"
28 m_print_inconsistency(
true),
29 m_free_parton_warnings(
true),
30 m_crash_on_problem(
false),
31 m_convert_gluon_to_0(
false),
35 m_store_weights(
true) {}
40 bool fill_next_event( Pythia8::Pythia& pythia,
GenEvent* evt,
int ievnum = -1 )
42 return fill_next_event( pythia.event, evt, ievnum, &pythia.info, &pythia.settings);
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)
51 bool fill_next_event( Pythia8::Event& pyev,
GenEvent* evt,
52 int ievnum = -1, Pythia8::Info* pyinfo = 0,
53 Pythia8::Settings* pyset = 0)
59 std::cerr <<
"Pythia8ToHepMC3::fill_next_event error - passed null event."
67 m_internal_event_number = ievnum;
71 ++m_internal_event_number;
77 std::vector<GenParticlePtr> hepevt_particles;
78 hepevt_particles.reserve( pyev.size() );
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() )
85 hepevt_particles[i]->set_generated_mass( pyev[i].m() );
89 std::vector<GenVertexPtr> vertex_cache;
90 std::vector<GenParticlePtr> beam_particles;
91 for(
int i=0; i<pyev.size(); ++i) {
93 std::vector<int> mothers = pyev[i].motherList();
96 GenVertexPtr prod_vtx = hepevt_particles[mothers[0]]->end_vertex();
99 prod_vtx = std::make_shared<GenVertex>();
100 vertex_cache.push_back(prod_vtx);
102 for(
unsigned int j=0; j<mothers.size(); ++j) {
103 prod_vtx->add_particle_in( hepevt_particles[mothers[j]] );
106 FourVector prod_pos( pyev[i].xProd(), pyev[i].yProd(),pyev[i].zProd(), pyev[i].tProd() );
109 if(!prod_pos.
is_zero() && prod_vtx->position().is_zero()) prod_vtx->set_position( prod_pos );
111 prod_vtx->add_particle_out( hepevt_particles[i] );
113 else beam_particles.push_back(hepevt_particles[i]);
117 evt->
reserve( hepevt_particles.size(), vertex_cache.size() );
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);
126 for(
int i=0; i<pyev.size(); ++i) {
129 int colType = pyev[i].colType();
130 if (colType == -1 ||colType == 1 || colType == 2)
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));
141 bool doHadr = (pyset == 0) ? m_free_parton_warnings : pyset->flag(
"HadronLevel:all") && pyset->flag(
"HadronLevel:Hadronize");
146 for (
int i = 1; i < pyev.size(); ++i) {
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] );
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);
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);
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;
181 GenPdfInfoPtr pdfinfo = std::make_shared<GenPdfInfo>();
182 pdfinfo->set(id1pdf, id2pdf, pyinfo->x1pdf(), pyinfo->x2pdf(), pyinfo->QFac(), pyinfo->pdf1(), pyinfo->pdf2() );
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()));
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);
204 if (m_store_weights && pyinfo != 0) {
206 for (
int iweight=0; iweight < pyinfo->nWeights(); ++iweight) {
207 evt->
weights().push_back(pyinfo->weight(iweight));
216 bool print_inconsistency()
const {
217 return m_print_inconsistency;
219 bool free_parton_warnings()
const {
220 return m_free_parton_warnings;
222 bool crash_on_problem()
const {
223 return m_crash_on_problem;
225 bool convert_gluon_to_0()
const {
226 return m_convert_gluon_to_0;
228 bool store_pdf()
const {
231 bool store_proc()
const {
234 bool store_xsec()
const {
237 bool store_weights()
const {
238 return m_store_weights;
242 void set_print_inconsistency(
bool b =
true) {
243 m_print_inconsistency = b;
245 void set_free_parton_warnings(
bool b =
true) {
246 m_free_parton_warnings = b;
248 void set_crash_on_problem(
bool b =
false) {
249 m_crash_on_problem = b;
251 void set_convert_gluon_to_0(
bool b =
false) {
252 m_convert_gluon_to_0 = b;
254 void set_store_pdf(
bool b =
true) {
257 void set_store_proc(
bool b =
true) {
260 void set_store_xsec(
bool b =
true) {
263 void set_store_weights(
bool b =
true) {
270 virtual bool fill_next_event(
GenEvent* ) {
273 virtual void write_event(
const GenEvent* ) {}
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;
287 bool m_store_weights;
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
bool is_zero() const
Check if the length of this vertex is zero.
Stores event-related information.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
void add_vertex(GenVertexPtr v)
Add vertex.
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
void set_event_number(const int &num)
Set event number.
void set_cross_section(GenCrossSectionPtr cs)
Set cross-section information.
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
const std::vector< double > & weights() const
Get event weight values as a vector.
void set_pdf_info(GenPdfInfoPtr pi)
Set PDF information.
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.