6 #ifndef Pythia8_Pythia8ToHepMC3_H
7 #define Pythia8_Pythia8ToHepMC3_H
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."
12 #include "Pythia8/Pythia.h"
26 m_print_inconsistency(
true),
27 m_free_parton_warnings(
true),
28 m_crash_on_problem(
false),
29 m_convert_gluon_to_0(
false),
33 m_store_weights(
true) {}
38 bool fill_next_event( Pythia8::Pythia& pythia,
GenEvent* evt,
int ievnum = -1 )
40 return fill_next_event( pythia.event, evt, ievnum, &pythia.info, &pythia.settings);
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)
49 bool fill_next_event( Pythia8::Event& pyev,
GenEvent* evt,
50 int ievnum = -1, Pythia8::Info* pyinfo = 0,
51 Pythia8::Settings* pyset = 0)
57 std::cerr <<
"Pythia8ToHepMC3::fill_next_event error - passed null event."
65 m_internal_event_number = ievnum;
69 ++m_internal_event_number;
75 std::vector<GenParticlePtr> hepevt_particles;
76 hepevt_particles.reserve( pyev.size() );
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() )
83 hepevt_particles[i]->set_generated_mass( pyev[i].m() );
87 std::vector<GenVertexPtr> vertex_cache;
88 std::vector<GenParticlePtr> beam_particles;
89 for(
int i=0; i<pyev.size(); ++i) {
91 std::vector<int> mothers = pyev[i].motherList();
94 GenVertexPtr prod_vtx = hepevt_particles[mothers[0]]->end_vertex();
97 prod_vtx = std::make_shared<GenVertex>();
98 vertex_cache.push_back(prod_vtx);
100 for(
unsigned int j=0; j<mothers.size(); ++j) {
101 prod_vtx->add_particle_in( hepevt_particles[mothers[j]] );
104 FourVector prod_pos( pyev[i].xProd(), pyev[i].yProd(),pyev[i].zProd(), pyev[i].tProd() );
107 if(!prod_pos.
is_zero() && prod_vtx->position().is_zero()) prod_vtx->set_position( prod_pos );
109 prod_vtx->add_particle_out( hepevt_particles[i] );
111 else beam_particles.push_back(hepevt_particles[i]);
115 evt->
reserve( hepevt_particles.size(), vertex_cache.size() );
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);
124 for(
int i=0; i<pyev.size(); ++i) {
127 int colType = pyev[i].colType();
128 if (colType == -1 ||colType == 1 || colType == 2)
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));
139 bool doHadr = (pyset == 0) ? m_free_parton_warnings : pyset->flag(
"HadronLevel:all") && pyset->flag(
"HadronLevel:Hadronize");
144 for (
int i = 1; i < pyev.size(); ++i) {
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] );
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);
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);
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;
179 GenPdfInfoPtr pdfinfo = std::make_shared<GenPdfInfo>();
180 pdfinfo->set(id1pdf, id2pdf, pyinfo->x1pdf(), pyinfo->x2pdf(), pyinfo->QFac(), pyinfo->pdf1(), pyinfo->pdf2() );
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()));
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);
202 if (m_store_weights && pyinfo != 0) {
204 for (
int iweight=0; iweight < pyinfo->nWeights(); ++iweight) {
205 evt->
weights().push_back(pyinfo->weight(iweight));
214 bool print_inconsistency()
const {
215 return m_print_inconsistency;
217 bool free_parton_warnings()
const {
218 return m_free_parton_warnings;
220 bool crash_on_problem()
const {
221 return m_crash_on_problem;
223 bool convert_gluon_to_0()
const {
224 return m_convert_gluon_to_0;
226 bool store_pdf()
const {
229 bool store_proc()
const {
232 bool store_xsec()
const {
235 bool store_weights()
const {
236 return m_store_weights;
240 void set_print_inconsistency(
bool b =
true) {
241 m_print_inconsistency = b;
243 void set_free_parton_warnings(
bool b =
true) {
244 m_free_parton_warnings = b;
246 void set_crash_on_problem(
bool b =
false) {
247 m_crash_on_problem = b;
249 void set_convert_gluon_to_0(
bool b =
false) {
250 m_convert_gluon_to_0 = b;
252 void set_store_pdf(
bool b =
true) {
255 void set_store_proc(
bool b =
true) {
258 void set_store_xsec(
bool b =
true) {
261 void set_store_weights(
bool b =
true) {
268 virtual bool fill_next_event(
GenEvent* ) {
271 virtual void write_event(
const GenEvent* ) {}
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;
285 bool m_store_weights;