6 #ifndef Pythia8_Pythia8ToHepMC3_H
7 #define Pythia8_Pythia8ToHepMC3_H
10 #include "Pythia8/Pythia.h"
24 m_print_inconsistency(
true),
25 m_free_parton_warnings(
true),
26 m_crash_on_problem(
false),
27 m_convert_gluon_to_0(
false),
31 m_store_weights(
true) {}
36 bool fill_next_event( Pythia8::Pythia& pythia,
GenEvent* evt,
int ievnum = -1 )
38 return fill_next_event( pythia.event, evt, ievnum, &pythia.info, &pythia.settings);
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)
47 bool fill_next_event( Pythia8::Event& pyev,
GenEvent* evt,
48 int ievnum = -1, Pythia8::Info* pyinfo = 0,
49 Pythia8::Settings* pyset = 0)
55 std::cerr <<
"Pythia8ToHepMC3::fill_next_event error - passed null event."
63 m_internal_event_number = ievnum;
67 ++m_internal_event_number;
73 std::vector<GenParticlePtr> hepevt_particles;
74 hepevt_particles.reserve( pyev.size() );
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() )
81 hepevt_particles[i]->set_generated_mass( pyev[i].m() );
85 std::vector<GenVertexPtr> vertex_cache;
86 std::vector<GenParticlePtr> beam_particles;
87 for(
int i=0; i<pyev.size(); ++i) {
89 std::vector<int> mothers = pyev[i].motherList();
92 GenVertexPtr prod_vtx = hepevt_particles[mothers[0]]->end_vertex();
95 prod_vtx = std::make_shared<GenVertex>();
96 vertex_cache.push_back(prod_vtx);
98 for(
unsigned int j=0; j<mothers.size(); ++j) {
99 prod_vtx->add_particle_in( hepevt_particles[mothers[j]] );
102 FourVector prod_pos( pyev[i].xProd(), pyev[i].yProd(),pyev[i].zProd(), pyev[i].tProd() );
105 if(!prod_pos.
is_zero() && prod_vtx->position().is_zero()) prod_vtx->set_position( prod_pos );
107 prod_vtx->add_particle_out( hepevt_particles[i] );
109 else beam_particles.push_back(hepevt_particles[i]);
113 evt->
reserve( hepevt_particles.size(), vertex_cache.size() );
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);
122 for(
int i=0; i<pyev.size(); ++i) {
125 int colType = pyev[i].colType();
126 if (colType == -1 ||colType == 1 || colType == 2)
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));
137 bool doHadr = (pyset == 0) ? m_free_parton_warnings : pyset->flag(
"HadronLevel:all") && pyset->flag(
"HadronLevel:Hadronize");
142 for (
int i = 1; i < pyev.size(); ++i) {
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] );
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);
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);
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;
177 GenPdfInfoPtr pdfinfo = std::make_shared<GenPdfInfo>();
178 pdfinfo->set(id1pdf, id2pdf, pyinfo->x1pdf(), pyinfo->x2pdf(), pyinfo->QFac(), pyinfo->pdf1(), pyinfo->pdf2() );
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()));
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);
200 if (m_store_weights && pyinfo != 0) {
202 for (
int iweight=0; iweight < pyinfo->nWeights(); ++iweight) {
203 evt->
weights().push_back(pyinfo->weight(iweight));
212 bool print_inconsistency()
const {
213 return m_print_inconsistency;
215 bool free_parton_warnings()
const {
216 return m_free_parton_warnings;
218 bool crash_on_problem()
const {
219 return m_crash_on_problem;
221 bool convert_gluon_to_0()
const {
222 return m_convert_gluon_to_0;
224 bool store_pdf()
const {
227 bool store_proc()
const {
230 bool store_xsec()
const {
233 bool store_weights()
const {
234 return m_store_weights;
238 void set_print_inconsistency(
bool b =
true) {
239 m_print_inconsistency = b;
241 void set_free_parton_warnings(
bool b =
true) {
242 m_free_parton_warnings = b;
244 void set_crash_on_problem(
bool b =
false) {
245 m_crash_on_problem = b;
247 void set_convert_gluon_to_0(
bool b =
false) {
248 m_convert_gluon_to_0 = b;
250 void set_store_pdf(
bool b =
true) {
253 void set_store_proc(
bool b =
true) {
256 void set_store_xsec(
bool b =
true) {
259 void set_store_weights(
bool b =
true) {
266 virtual bool fill_next_event(
GenEvent* ) {
269 virtual void write_event(
const GenEvent* ) {}
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;
283 bool m_store_weights;
void set_cross_section(GenCrossSectionPtr cs)
Set cross-section information.
void add_vertex(GenVertexPtr v)
Add vertex.
Definition of class GenParticle.
Definition of class GenVertex.
bool is_zero() const
Check if the length of this vertex is zero.
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
Stores event-related information.
void set_pdf_info(GenPdfInfoPtr pi)
Set PDF information.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition of class GenEvent.
Feature< Feature_type > abs(const Feature< Feature_type > &input)
Obtain the absolute value of a Feature. This works as you'd expect. If foo is a valid Feature...
void set_event_number(const int &num)
Set event number.