1 // HepMCInterface.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2013 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 // Author: Mikhail Kirsanov, Mikhail.Kirsanov@cern.ch
7 // Function definitions (not found in the header) for the I_Pythia8 class,
8 // which converts a PYTHIA event record to the standard HepMC format.
10 #include "HepMCInterface.h"
11 #include "HepMC/GenEvent.h"
15 //==========================================================================
17 // Main method for conversion from PYTHIA event to HepMC event.
18 // Read one event from Pythia8 and fill GenEvent,
19 // and return T/F = success/failure.
21 bool I_Pythia8::fill_next_event( Pythia8::Event& pyev, GenEvent* evt,
22 int ievnum, Pythia8::Info* pyinfo, Pythia8::Settings* pyset) {
24 // 1. Error if no event passed.
26 std::cerr << "I_Pythia8::fill_next_event error - passed null event."
31 // Event number counter.
33 evt->set_event_number(ievnum);
34 m_internal_event_number = ievnum;
36 evt->set_event_number(m_internal_event_number);
37 ++m_internal_event_number;
40 // Conversion factors from Pythia units GeV and mm to HepMC ones.
41 double momFac = HepMC::Units::conversion_factor(HepMC::Units::GEV,
42 evt->momentum_unit());
43 double lenFac = HepMC::Units::conversion_factor(HepMC::Units::MM,
46 // 2. Create a particle instance for each entry and fill a map, and
47 // a vector which maps from the particle index to the GenParticle address.
48 std::vector<GenParticle*> hepevt_particles( pyev.size() );
49 for (int i = 1; i < pyev.size(); ++i) {
52 hepevt_particles[i] = new GenParticle(
53 FourVector( momFac * pyev[i].px(), momFac * pyev[i].py(),
54 momFac * pyev[i].pz(), momFac * pyev[i].e() ),
55 pyev[i].id(), pyev.statusHepMC(i) );
56 hepevt_particles[i]->suggest_barcode(i);
57 hepevt_particles[i]->set_generated_mass( momFac * pyev[i].m() );
59 // Colour flow uses index 1 and 2.
60 int colType = pyev[i].colType();
61 if (colType == 1 || colType == 2)
62 hepevt_particles[i]->set_flow(1, pyev[i].col());
63 if (colType == -1 || colType == 2)
64 hepevt_particles[i]->set_flow(2, pyev[i].acol());
67 // Here we assume that the first two particles in the list
68 // are the incoming beam particles.
69 evt->set_beam_particles( hepevt_particles[1], hepevt_particles[2] );
71 // 3. Loop over particles AGAIN, this time creating vertices.
72 // We build the production vertex for each entry in hepevt.
73 // The HEPEVT pointers are bi-directional, so gives decay vertices as well.
74 for (int i = 1; i < pyev.size(); ++i) {
75 GenParticle *p = hepevt_particles[i];
77 // 3a. Search to see if a production vertex already exists.
78 std::vector<int> mothers = pyev.motherList(i);
79 unsigned int imother = 0;
80 int mother = -1; // note that in Pythia8 there is a particle number 0!
81 if ( !mothers.empty() ) mother = mothers[imother];
82 GenVertex* prod_vtx = p->production_vertex();
83 while ( !prod_vtx && mother > 0 ) {
84 prod_vtx = hepevt_particles[mother]->end_vertex();
85 if ( prod_vtx ) prod_vtx->add_particle_out( p );
86 mother = ( ++imother < mothers.size() ) ? mothers[imother] : -1;
89 // 3b. If no suitable production vertex exists - and the particle has
90 // at least one mother or position information to store - make one.
91 FourVector prod_pos( lenFac * pyev[i].xProd(), lenFac * pyev[i].yProd(),
92 lenFac * pyev[i].zProd(), lenFac * pyev[i].tProd() );
93 if ( !prod_vtx && ( mothers.size() > 0 || prod_pos != FourVector() ) ) {
94 prod_vtx = new GenVertex();
95 prod_vtx->add_particle_out( p );
96 evt->add_vertex( prod_vtx );
99 // 3c. If prod_vtx doesn't already have position specified, fill it.
100 if ( prod_vtx && prod_vtx->position() == FourVector() )
101 prod_vtx->set_position( prod_pos );
103 // 3d. loop over mothers to make sure their end_vertices are consistent.
106 if ( !mothers.empty() ) mother = mothers[imother];
107 while ( prod_vtx && mother > 0 ) {
109 // If end vertex of the mother isn't specified, do it now.
110 if ( !hepevt_particles[mother]->end_vertex() ) {
111 prod_vtx->add_particle_in( hepevt_particles[mother] );
113 // Problem scenario: the mother already has a decay vertex which
114 // differs from the daughter's production vertex. This means there is
115 // internal inconsistency in the HEPEVT event record. Print an error.
116 // Note: we could provide a fix by joining the two vertices with a
117 // dummy particle if the problem arises often.
118 } else if (hepevt_particles[mother]->end_vertex() != prod_vtx ) {
119 if ( m_print_inconsistency ) std::cerr
120 << "HepMC::I_Pythia8: inconsistent mother/daugher "
121 << "information in Pythia8 event " << std::endl
122 << "i = " << i << " mother = " << mother
123 << "\n This warning can be turned off with the "
124 << "I_Pythia8::print_inconsistency switch." << std::endl;
127 // End of vertex-setting loops.
128 mother = ( ++imother < mothers.size() ) ? mothers[imother] : -1;
132 // If hadronization switched on then no final coloured particles.
133 bool doHadr = (pyset == 0) ? m_free_parton_warnings
134 : pyset->flag("HadronLevel:all") && pyset->flag("HadronLevel:Hadronize");
136 // 4. Check for particles which come from nowhere, i.e. are without
137 // mothers or daughters. These need to be attached to a vertex, or else
138 // they will never become part of the event.
139 for (int i = 1; i < pyev.size(); ++i) {
140 if ( !hepevt_particles[i]->end_vertex() &&
141 !hepevt_particles[i]->production_vertex() ) {
142 std::cerr << "hanging particle " << i << std::endl;
143 GenVertex* prod_vtx = new GenVertex();
144 prod_vtx->add_particle_out( hepevt_particles[i] );
145 evt->add_vertex( prod_vtx );
148 // Also check for free partons (= gluons and quarks; not diquarks?).
149 if ( doHadr && m_free_parton_warnings ) {
150 if ( hepevt_particles[i]->pdg_id() == 21 &&
151 !hepevt_particles[i]->end_vertex() ) {
152 std::cerr << "gluon without end vertex " << i << std::endl;
153 if ( m_crash_on_problem ) exit(1);
155 if ( abs(hepevt_particles[i]->pdg_id()) <= 6 &&
156 !hepevt_particles[i]->end_vertex() ) {
157 std::cerr << "quark without end vertex " << i << std::endl;
158 if ( m_crash_on_problem ) exit(1);
163 // 5. Store PDF, weight, cross section and other event information.
164 // Flavours of incoming partons.
165 if (m_store_pdf && pyinfo != 0) {
166 int id1pdf = pyinfo->id1pdf();
167 int id2pdf = pyinfo->id2pdf();
168 if ( m_convert_gluon_to_0 ) {
169 if (id1pdf == 21) id1pdf = 0;
170 if (id2pdf == 21) id2pdf = 0;
173 // Store PDF information.
174 evt->set_pdf_info( PdfInfo( id1pdf, id2pdf, pyinfo->x1pdf(),
175 pyinfo->x2pdf(), pyinfo->QFac(), pyinfo->pdf1(), pyinfo->pdf2() ) );
178 // Store process code, scale, alpha_em, alpha_s.
179 if (m_store_proc && pyinfo != 0) {
180 evt->set_signal_process_id( pyinfo->code() );
181 evt->set_event_scale( pyinfo->QRen() );
182 if (evt->alphaQED() <= 0) evt->set_alphaQED( pyinfo->alphaEM() );
183 if (evt->alphaQCD() <= 0) evt->set_alphaQCD( pyinfo->alphaS() );
186 // Store cross-section information in pb and event weight. The latter is
187 // usually dimensionless, but in units of pb for Les Houches strategies +-4.
188 if (m_store_xsec && pyinfo != 0) {
189 HepMC::GenCrossSection xsec;
190 xsec.set_cross_section( pyinfo->sigmaGen() * 1e9, pyinfo->sigmaErr() * 1e9);
191 evt->set_cross_section(xsec);
192 evt->weights().push_back( pyinfo->weight() );
200 //==========================================================================
202 } // end namespace HepMC