File indexing completed on 2024-04-06 12:14:06
0001
0002
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
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
0043
0044
0045
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
0060
0061 if (trust_beam_particles()) {
0062 evt->set_beam_particles(hepevt_particle[1], hepevt_particle[2]);
0063 }
0064
0065
0066
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
0071
0072
0073
0074
0075
0076 if (m_trust_mothers_before_daughters || m_trust_both_mothers_and_daughters) {
0077 build_production_vertex(i, hepevt_particle, evt);
0078 }
0079
0080
0081
0082 if (!m_trust_mothers_before_daughters || m_trust_both_mothers_and_daughters) {
0083 build_end_vertex(i, hepevt_particle, evt);
0084 }
0085 }
0086
0087
0088
0089
0090
0091
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
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
0116
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
0127
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
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
0155
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
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
0184 if (++mother > EPOS_Wrapper::last_parent(i))
0185 mother = 0;
0186 }
0187
0188
0189
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
0197 if (prod_vtx && prod_vtx->position() == FourVector(0, 0, 0, 0)) {
0198 prod_vtx->set_position(prod_pos);
0199 }
0200
0201
0202 mother = EPOS_Wrapper::first_parent(i);
0203 while (prod_vtx && mother > 0) {
0204 if (!hepevt_particle[mother]->end_vertex()) {
0205
0206 prod_vtx->add_particle_in(hepevt_particle[mother]);
0207 } else if (hepevt_particle[mother]->end_vertex() != prod_vtx) {
0208
0209
0210
0211
0212
0213
0214
0215
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
0231 HepMC::GenParticle* p = hepevt_particle[i];
0232
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
0243
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
0250
0251
0252 daughter = EPOS_Wrapper::first_child(i);
0253 while (end_vtx && daughter > 0) {
0254 if (!hepevt_particle[daughter]->production_vertex()) {
0255
0256 end_vtx->add_particle_out(hepevt_particle[daughter]);
0257
0258
0259 if (end_vtx->position() == FourVector(0, 0, 0, 0)) {
0260
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
0271
0272
0273
0274
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
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 }