PYTHIA8/pythia8175/hepmcinterface/HepMCInterface.cc

   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.
   5
   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.
   9
  10 #include "HepMCInterface.h"
  11 #include "HepMC/GenEvent.h"
  12
  13 namespace HepMC {
  14
  15 //==========================================================================
  16
  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.
  20
  21 bool I_Pythia8::fill_next_event( Pythia8::Event& pyev, GenEvent* evt,
  22     int ievnum, Pythia8::Info* pyinfo, Pythia8::Settings* pyset) {
  23
  24   // 1. Error if no event passed.
  25   if (!evt) {
  26     std::cerr << "I_Pythia8::fill_next_event error - passed null event."
  27               << std::endl;
  28     return 0;
  29   }
  30
  31   // Event number counter.
  32   if ( ievnum >= 0 ) {
  33     evt->set_event_number(ievnum);
  34     m_internal_event_number = ievnum;
  35   } else {
  36     evt->set_event_number(m_internal_event_number);
  37     ++m_internal_event_number;
  38   }
  39
  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,
  44     evt->length_unit());
  45
  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) {
  50
  51     // Fill the particle.
  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() );
  58
  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());
  65   }
  66
  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] );
  70
  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];
  76
  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;
  87     }
  88
  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 );
  97     }
  98
  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 );
 102
 103     // 3d. loop over mothers to make sure their end_vertices are consistent.
 104     imother = 0;
 105     mother = -1;
 106     if ( !mothers.empty() ) mother = mothers[imother];
 107     while ( prod_vtx && mother > 0 ) {
 108
 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] );
 112
 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;
 125       }
 126
 127       // End of vertex-setting loops.
 128       mother = ( ++imother < mothers.size() ) ? mothers[imother] : -1;
 129     }
 130   }
 131
 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");
 135
 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 );
 146     }
 147
 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);
 154       }
 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);
 159       }
 160     }
 161   }
 162
 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;
 171     }
 172
 173     // Store PDF information.
 174     evt->set_pdf_info( PdfInfo( id1pdf, id2pdf, pyinfo->x1pdf(),
 175       pyinfo->x2pdf(), pyinfo->QFac(), pyinfo->pdf1(), pyinfo->pdf2() ) );
 176   }
 177
 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() );
 184   }
 185
 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() );
 193   }
 194
 195   // Done.
 196   return true;
 197
 198 }
 199
 200 //==========================================================================
 201
 202 } // end namespace HepMC