ReggeGribovPartonMCInterface/src/IO_EPOS.cc

0001
0002 // EPOS IO class
0003
0004 #include "GeneratorInterface/ReggeGribovPartonMCInterface/interface/IO_EPOS.h"
0005 #include "HepMC/GenEvent.h"
0006 #include <cstdio>  // needed for formatted output using sprintf
0007
0008 using namespace HepMC;
0009
0010 namespace EPOS {
0011
0012   unsigned int EPOS_Wrapper::s_sizeof_int = 4;
0013   unsigned int EPOS_Wrapper::s_sizeof_real = sizeof(double);
0014   unsigned int EPOS_Wrapper::s_max_number_entries = 99900;
0015
0016   IO_EPOS::IO_EPOS()
0017       : m_trust_mothers_before_daughters(true),
0018         m_trust_both_mothers_and_daughters(false),
0019         m_print_inconsistency_errors(true),
0020         m_trust_beam_particles(true),
0021         m_skip_nucl_frag(false) {}
0022
0023   IO_EPOS::~IO_EPOS() {}
0024
0025   void IO_EPOS::print(std::ostream& ostr) const {
0026     ostr << "IO_EPOS: reads an event from the FORTRAN EPOS g "
0027          << "common block. \n"
0028          << " trust_mothers_before_daughters = " << m_trust_mothers_before_daughters
0029          << " trust_both_mothers_and_daughters = " << m_trust_both_mothers_and_daughters
0030          << ", print_inconsistency_errors = " << m_print_inconsistency_errors << std::endl;
0031   }
0032
0033   bool IO_EPOS::fill_next_event(HepMC::GenEvent* evt) {
0034     //
0035     // 1. test that evt pointer is not null and set event number
0036     if (!evt) {
0037       std::cerr << "IO_EPOS::fill_next_event error - passed null event." << std::endl;
0038       return false;
0039     }
0040     evt->set_event_number(EPOS_Wrapper::event_number());
0041     //
0042     // 2. create a particle instance for each EPOS entry and fill a map
0043     //    create a vector which maps from the EPOS particle index to the
0044     //    GenParticle address
0045     //    (+1 in size accounts for hepevt_particle[0] which is unfilled)
0046     std::vector<HepMC::GenParticle*> hepevt_particle(EPOS_Wrapper::number_entries() + 1);
0047     hepevt_particle[0] = nullptr;
0048     for (int i1 = 1; i1 <= EPOS_Wrapper::number_entries(); ++i1) {
0049       hepevt_particle[i1] = build_particle(i1);
0050     }
0051
0052     HepMC::GenVertex* primaryVertex = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), 0);
0053     evt->add_vertex(primaryVertex);
0054     if (!evt->signal_process_vertex())
0055       evt->set_signal_process_vertex(primaryVertex);
0056
0057     std::set<HepMC::GenVertex*> new_vertices;
0058     //
0059     // Here we assume that the first two particles in the list
0060     // are the incoming beam particles.
0061     if (trust_beam_particles()) {
0062       evt->set_beam_particles(hepevt_particle[1], hepevt_particle[2]);
0063     }
0064     //
0065     // 3.+4. loop over EPOS particles AGAIN, this time creating vertices
0066     //MODIFICATION FROM HEPMC!! skipping nuclear fragments in event if option is set
0067     for (int i = 1; i <= EPOS_Wrapper::number_entries(); ++i) {
0068       if (m_skip_nucl_frag && abs(hepevt_particle[i]->pdg_id()) >= 1000000000)
0069         continue;
0070       // We go through and build EITHER the production or decay
0071       // vertex for each entry in hepevt, depending on the switch
0072       // m_trust_mothers_before_daughters (new 2001-02-28)
0073       // Note: since the EPOS pointers are bi-directional, it is
0074       //
0075       // 3. Build the production_vertex (if necessary)
0076       if (m_trust_mothers_before_daughters || m_trust_both_mothers_and_daughters) {
0077         build_production_vertex(i, hepevt_particle, evt);
0078       }
0079       //
0080       // 4. Build the end_vertex (if necessary)
0081       //    Identical steps as for production vertex
0082       if (!m_trust_mothers_before_daughters || m_trust_both_mothers_and_daughters) {
0083         build_end_vertex(i, hepevt_particle, evt);
0084       }
0085     }
0086     // 5.             01.02.2000
0087     // handle the case of particles in EPOS which come from nowhere -
0088     //  i.e. particles without mothers or daughters.
0089     //  These particles need to be attached to a vertex, or else they
0090     //  will never become part of the event. check for this situation
0091     //MODIFICATION FROM HEPMC!! skipping nuclear fragments in event if option is set
0092     for (int i3 = 1; i3 <= EPOS_Wrapper::number_entries(); ++i3) {
0093       if (m_skip_nucl_frag && abs(hepevt_particle[i3]->pdg_id()) >= 1000000000)
0094         continue;
0095       if (!hepevt_particle[i3]->end_vertex() && !hepevt_particle[i3]->production_vertex()) {
0096         HepMC::GenVertex* prod_vtx = new GenVertex();
0097         prod_vtx->add_particle_out(hepevt_particle[i3]);
0098         evt->add_vertex(prod_vtx);
0099       }
0100     }
0101     return true;
0102   }
0103
0104   void IO_EPOS::write_event(const GenEvent* evt) {
0105     //
0106     if (!evt)
0107       return;
0108     //
0109     // map all particles onto a unique index
0110     std::vector<HepMC::GenParticle*> index_to_particle(EPOS_Wrapper::max_number_entries() + 1);
0111     index_to_particle[0] = nullptr;
0112     std::map<HepMC::GenParticle*, int> particle_to_index;
0113     int particle_counter = 0;
0114     for (HepMC::GenEvent::vertex_const_iterator v = evt->vertices_begin(); v != evt->vertices_end(); ++v) {
0115       // all "mothers" or particles_in are kept adjacent in the list
0116       // so that the mother indices in hepevt can be filled properly
0117       for (HepMC::GenVertex::particles_in_const_iterator p1 = (*v)->particles_in_const_begin();
0118            p1 != (*v)->particles_in_const_end();
0119            ++p1) {
0120         ++particle_counter;
0121         if (particle_counter > EPOS_Wrapper::max_number_entries())
0122           break;
0123         index_to_particle[particle_counter] = *p1;
0124         particle_to_index[*p1] = particle_counter;
0125       }
0126       // daughters are entered only if they aren't a mother of
0127       // another vtx
0128       for (HepMC::GenVertex::particles_out_const_iterator p2 = (*v)->particles_out_const_begin();
0129            p2 != (*v)->particles_out_const_end();
0130            ++p2) {
0131         if (!(*p2)->end_vertex()) {
0132           ++particle_counter;
0133           if (particle_counter > EPOS_Wrapper::max_number_entries()) {
0134             break;
0135           }
0136           index_to_particle[particle_counter] = *p2;
0137           particle_to_index[*p2] = particle_counter;
0138         }
0139       }
0140     }
0141     if (particle_counter > EPOS_Wrapper::max_number_entries()) {
0142       particle_counter = EPOS_Wrapper::max_number_entries();
0143     }
0144     //
0145     // fill the EPOS event record
0146     EPOS_Wrapper::set_event_number(evt->event_number());
0147     EPOS_Wrapper::set_number_entries(particle_counter);
0148     for (int i = 1; i <= particle_counter; ++i) {
0149       EPOS_Wrapper::set_status(i, index_to_particle[i]->status());
0150       EPOS_Wrapper::set_id(i, index_to_particle[i]->pdg_id());
0151       FourVector m = index_to_particle[i]->momentum();
0152       EPOS_Wrapper::set_momentum(i, m.px(), m.py(), m.pz(), m.e());
0153       EPOS_Wrapper::set_mass(i, index_to_particle[i]->generatedMass());
0154       // there should ALWAYS be particles in any vertex, but some generators
0155       // are making non-kosher HepMC events
0156       if (index_to_particle[i]->production_vertex() && index_to_particle[i]->production_vertex()->particles_in_size()) {
0157         FourVector p = index_to_particle[i]->production_vertex()->position();
0158         EPOS_Wrapper::set_position(i, p.x(), p.y(), p.z(), p.t());
0159         int num_mothers = index_to_particle[i]->production_vertex()->particles_in_size();
0160         int first_mother =
0161             find_in_map(particle_to_index, *(index_to_particle[i]->production_vertex()->particles_in_const_begin()));
0162         int last_mother = first_mother + num_mothers - 1;
0163         if (first_mother == 0)
0164           last_mother = 0;
0165         EPOS_Wrapper::set_parents(i, first_mother, last_mother);
0166       } else {
0167         EPOS_Wrapper::set_position(i, 0, 0, 0, 0);
0168         EPOS_Wrapper::set_parents(i, 0, 0);
0169       }
0170       EPOS_Wrapper::set_children(i, 0, 0);
0171     }
0172   }
0173
0174   void IO_EPOS::build_production_vertex(int i, std::vector<HepMC::GenParticle*>& hepevt_particle, GenEvent* evt) {
0175     HepMC::GenParticle* p = hepevt_particle[i];
0176     // a. search to see if a production vertex already exists
0177     int mother = EPOS_Wrapper::first_parent(i);
0178     HepMC::GenVertex* prod_vtx = p->production_vertex();
0179     while (!prod_vtx && mother > 0) {
0180       prod_vtx = hepevt_particle[mother]->end_vertex();
0181       if (prod_vtx)
0182         prod_vtx->add_particle_out(p);
0183       // increment mother for next iteration
0184       if (++mother > EPOS_Wrapper::last_parent(i))
0185         mother = 0;
0186     }
0187     // b. if no suitable production vertex exists - and the particle
0188     // has atleast one mother or position information to store -
0189     // make one
0190     HepMC::FourVector prod_pos(EPOS_Wrapper::x(i), EPOS_Wrapper::y(i), EPOS_Wrapper::z(i), EPOS_Wrapper::t(i));
0191     if (!prod_vtx && (EPOS_Wrapper::number_parents(i) > 0 || prod_pos != FourVector(0, 0, 0, 0))) {
0192       prod_vtx = new HepMC::GenVertex();
0193       prod_vtx->add_particle_out(p);
0194       evt->add_vertex(prod_vtx);
0195     }
0196     // c. if prod_vtx doesn't already have position specified, fill it
0197     if (prod_vtx && prod_vtx->position() == FourVector(0, 0, 0, 0)) {
0198       prod_vtx->set_position(prod_pos);
0199     }
0200     // d. loop over mothers to make sure their end_vertices are
0201     //     consistent
0202     mother = EPOS_Wrapper::first_parent(i);
0203     while (prod_vtx && mother > 0) {
0204       if (!hepevt_particle[mother]->end_vertex()) {
0205         // if end vertex of the mother isn't specified, do it now
0206         prod_vtx->add_particle_in(hepevt_particle[mother]);
0207       } else if (hepevt_particle[mother]->end_vertex() != prod_vtx) {
0208         // problem scenario --- the mother already has a decay
0209         // vertex which differs from the daughter's produciton
0210         // vertex. This means there is internal
0211         // inconsistency in the EPOS event record. Print an
0212         // error
0213         // Note: we could provide a fix by joining the two
0214         //       vertices with a dummy particle if the problem
0215         //       arrises often with any particular generator.
0216         if (m_print_inconsistency_errors)
0217           std::cerr << "HepMC::IO_EPOS: inconsistent mother/daugher "
0218                     << "information in EPOS event " << EPOS_Wrapper::event_number() << ". \n I recommend you try "
0219                     << "inspecting the event first with "
0220                     << "\n\tEPOS_Wrapper::check_hepevt_consistency()"
0221                     << "\n This warning can be turned off with the "
0222                     << "IO_EPOS::print_inconsistency_errors switch." << std::endl;
0223       }
0224       if (++mother > EPOS_Wrapper::last_parent(i))
0225         mother = 0;
0226     }
0227   }
0228
0229   void IO_EPOS::build_end_vertex(int i, std::vector<HepMC::GenParticle*>& hepevt_particle, GenEvent* evt) {
0230     //    Identical steps as for build_production_vertex
0231     HepMC::GenParticle* p = hepevt_particle[i];
0232     // a.
0233     int daughter = EPOS_Wrapper::first_child(i);
0234     HepMC::GenVertex* end_vtx = p->end_vertex();
0235     while (!end_vtx && daughter > 0) {
0236       end_vtx = hepevt_particle[daughter]->production_vertex();
0237       if (end_vtx)
0238         end_vtx->add_particle_in(p);
0239       if (++daughter > EPOS_Wrapper::last_child(i))
0240         daughter = 0;
0241     }
0242     // b. (different from 3c. because EPOS particle can not know its
0243     //        decay position )
0244     if (!end_vtx && EPOS_Wrapper::number_children(i) > 0) {
0245       end_vtx = new GenVertex();
0246       end_vtx->add_particle_in(p);
0247       evt->add_vertex(end_vtx);
0248     }
0249     // c+d. loop over daughters to make sure their production vertices
0250     //    point back to the current vertex.
0251     //    We get the vertex position from the daughter as well.
0252     daughter = EPOS_Wrapper::first_child(i);
0253     while (end_vtx && daughter > 0) {
0254       if (!hepevt_particle[daughter]->production_vertex()) {
0255         // if end vertex of the mother isn't specified, do it now
0256         end_vtx->add_particle_out(hepevt_particle[daughter]);
0257         //
0258         // 2001-03-29 M.Dobbs, fill vertex the position.
0259         if (end_vtx->position() == FourVector(0, 0, 0, 0)) {
0260           // again mm to cm conversion
0261           FourVector prod_pos(EPOS_Wrapper::x(daughter),
0262                               EPOS_Wrapper::y(daughter),
0263                               EPOS_Wrapper::z(daughter),
0264                               EPOS_Wrapper::t(daughter));
0265           if (prod_pos != FourVector(0, 0, 0, 0)) {
0266             end_vtx->set_position(prod_pos);
0267           }
0268         }
0269       } else if (hepevt_particle[daughter]->production_vertex() != end_vtx) {
0270         // problem scenario --- the daughter already has a prod
0271         // vertex which differs from the mother's end
0272         // vertex. This means there is internal
0273         // inconsistency in the EPOS event record. Print an
0274         // error
0275         if (m_print_inconsistency_errors)
0276           std::cerr << "HepMC::IO_EPOS: inconsistent mother/daugher "
0277                     << "information in EPOS event " << EPOS_Wrapper::event_number() << ". \n I recommend you try "
0278                     << "inspecting the event first with "
0279                     << "\n\tEPOS_Wrapper::check_hepevt_consistency()"
0280                     << "\n This warning can be turned off with the "
0281                     << "IO_EPOS::print_inconsistency_errors switch." << std::endl;
0282       }
0283       if (++daughter > EPOS_Wrapper::last_child(i))
0284         daughter = 0;
0285     }
0286     if (!p->end_vertex() && !p->production_vertex()) {
0287       // Added 2001-11-04, to try and handle Isajet problems.
0288       build_production_vertex(i, hepevt_particle, evt);
0289     }
0290   }
0291
0292   HepMC::GenParticle* IO_EPOS::build_particle(int index) {
0293     //
0294     HepMC::GenParticle* p = new GenParticle(
0295         FourVector(EPOS_Wrapper::px(index), EPOS_Wrapper::py(index), EPOS_Wrapper::pz(index), EPOS_Wrapper::e(index)),
0296         EPOS_Wrapper::id(index),
0297         EPOS_Wrapper::status(index));
0298     p->setGeneratedMass(EPOS_Wrapper::m(index));
0299     p->suggest_barcode(index);
0300     return p;
0301   }
0302
0303   int IO_EPOS::find_in_map(const std::map<HepMC::GenParticle*, int>& m, HepMC::GenParticle* p) const {
0304     std::map<HepMC::GenParticle*, int>::const_iterator iter = m.find(p);
0305     if (iter == m.end())
0306       return 0;
0307     return iter->second;
0308   }
0309
0310 }  // namespace EPOS