include/HepMC3/HEPEVT_Helpers.h

0001 // -*- C++ -*-
0002 //
0003 // This file is part of HepMC
0004 // Copyright (C) 2014-2023 The HepMC collaboration (see AUTHORS for details)
0005 //
0006 #ifndef HEPEVT_HELPERS_H
0007 #define HEPEVT_HELPERS_H
0008 /**
0009  *  @file HEPEVT_Helpers.h
0010  *  @brief Helper functions used to manipulate with HEPEVT block
0011  *
0012  */
0013 #if defined(WIN32)&&(!defined(HEPMC3_NO_EXPORTS))
0014 #if defined(HepMC3_EXPORTS)
0015 #define HEPMC3_EXPORT_API __declspec(dllexport)
0016 #else
0017 #define HEPMC3_EXPORT_API __declspec(dllimport)
0018 #endif
0019 #else
0020 #define HEPMC3_EXPORT_API
0021 #endif
0022 #include <algorithm>
0023 #include <map>
0024
0025 #include "HepMC3/GenEvent.h"
0026 #include "HepMC3/GenParticle.h"
0027 #include "HepMC3/GenVertex.h"
0028 namespace HepMC3
0029 {
0030
0031 /** @struct HEPEVT_Templated
0032  *  @brief  C structure representing Fortran common block HEPEVT
0033  * T. Sjöstrand et al., "A proposed standard event record",
0034  *  in `Z physics at LEP 1', eds. G. Altarelli, R. Kleiss and C. Verzegnassi,
0035  * Geneva, Switzerland, September 4-5, 1989, CERN 89-08 (Geneva, 1989), Vol. 3, p. 327
0036  * Disk representation is given by Fortran WRITE/READ format.
0037  */
0038 template <int max_particles, typename momentum_type = double>
0039 struct HEPEVT_Templated
0040 {
0041     int        nevhep;             //!< Event number
0042     int        nhep;               //!< Number of entries in the event
0043     int        isthep[max_particles];     //!< Status code
0044     int        idhep [max_particles];     //!< PDG ID
0045     int        jmohep[max_particles][2];  //!< Position of 1st and 2nd (or last!) mother
0046     int        jdahep[max_particles][2];  //!< Position of 1nd and 2nd (or last!) daughter
0047     momentum_type phep  [max_particles][5];  //!< Momentum: px, py, pz, e, m
0048     momentum_type vhep  [max_particles][4];  //!< Time-space position: x, y, z, t
0049 };
0050
0051 /** @struct HEPEVT_Pointers
0052  *  @brief  C structure representing Fortran common block HEPEVT
0053  * T. Sjöstrand et al., "A proposed standard event record",
0054  *  in `Z physics at LEP 1', eds. G. Altarelli, R. Kleiss and C. Verzegnassi,
0055  * Geneva, Switzerland, September 4-5, 1989, CERN 89-08 (Geneva, 1989), Vol. 3, p. 327
0056  * Disk representation is given by Fortran WRITE/READ format.
0057  */
0058 template<typename momentum_type = double>
0059 struct HEPEVT_Pointers
0060 {
0061     int*        nevhep;             //!< Pointer to Event number
0062     int*        nhep;               //!< Pointer to Number of entries in the event
0063     int*        isthep;     //!< Pointer to Status code
0064     int*        idhep;      //!< Pointer to PDG ID
0065     int*        jmohep;     //!< Pointer to position of 1st and 2nd (or last!) mother
0066     int*        jdahep;     //!< Pointer to position of 1nd and 2nd (or last!) daughter
0067     momentum_type* phep;    //!< Pointer to momentum: px, py, pz, e, m
0068     momentum_type* vhep;    //!< Pointer to time-space position: x, y, z, t
0069 };
0070
0071 /** @brief comparison of two particles */
0072 struct GenParticlePtr_greater
0073 {
0074     /** @brief comparison of two particles */
0075     bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const;
0076 };
0077 /** @brief  Order vertices with equal paths. */
0078 struct pair_GenVertexPtr_int_greater
0079 {
0080     /** @brief  Order vertices with equal paths. If the paths are equal, order in other quantities.
0081      * We cannot use id, as it can be assigned in different way*/
0082     bool operator()(const std::pair<ConstGenVertexPtr, int>& lx, const std::pair<ConstGenVertexPtr, int>& rx) const;
0083 };
0084 /** @brief Calculates the path to the top (beam) particles */
0085 void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map<ConstGenVertexPtr, int>& pathl);
0086
0087
0088 /** @brief  Converts HEPEVT into GenEvent. */
0089 template <class T>
0090 bool HEPEVT_to_GenEvent_nonstatic(GenEvent* evt, T* A)
0091 {
0092     if ( !evt ) { std::cerr << "HEPEVT_to_GenEvent_nonstatic  - passed null event." << std::endl; return false;}
0093     evt->set_event_number(A->event_number());
0094     std::map<GenParticlePtr, int > hepevt_particles;
0095     std::map<int, GenParticlePtr > particles_index;
0096     std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > hepevt_vertices;
0097     std::map<int, GenVertexPtr > vertex_index;
0098     int ne=0;
0099     ne=A->number_entries();
0100     for ( int i = 1; i <= ne; i++ )
0101     {
0102         GenParticlePtr p = std::make_shared<GenParticle>();
0103         p->set_momentum(FourVector(A->px(i), A->py(i), A->pz(i), A->e(i)));
0104         p->set_status(A->status(i));
0105         p->set_pid(A->id(i)); //Confusing!
0106         p->set_generated_mass(A->m(i));
0107         hepevt_particles[p] = i;
0108         particles_index[i] = p;
0109         GenVertexPtr v = std::make_shared<GenVertex>();
0110         v->set_position(FourVector(A->x(i), A->y(i), A->z(i), A->t(i)));
0111         v->add_particle_out(p);
0112         std::set<int> in;
0113         std::set<int> out;
0114         out.insert(i);
0115         vertex_index[i] = v;
0116         hepevt_vertices[v] = std::pair<std::set<int>, std::set<int> >(in, out);
0117     }
0118     /* The part above is always correct as it is a raw information without any topology.*/
0119
0120     /* In this way we trust mother information. The "Trust daughters" is not implemented.*/
0121     for (std::map<GenParticlePtr, int >::iterator it1 = hepevt_particles.begin(); it1 != hepevt_particles.end(); ++it1)
0122         for (std::map<GenParticlePtr, int >::iterator it2 = hepevt_particles.begin(); it2 != hepevt_particles.end(); ++it2) {
0123             if   (A->first_parent(it2->second) <= it1->second && it1->second <= A->last_parent(it2->second)) hepevt_vertices[it2->first->production_vertex()].first.insert(it1->second);
0124         }
0125     /* Now all incoming sets are correct for all vertices. But we have duplicates.*/
0126
0127     /* Disconnect all particles from the vertices*/
0128     for ( int i = 1; i <= A->number_entries(); i++ ) vertex_index[i]->remove_particle_out(particles_index[i]);
0129
0130     /*Fill container with vertices with unique sets of incoming particles. Merge the outcoming particle sets.*/
0131     std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > final_vertices_map;
0132     for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator vs = hepevt_vertices.begin(); vs != hepevt_vertices.end(); ++vs)
0133     {
0134         if ((final_vertices_map.size() == 0) || (vs->second.first.size() == 0 && vs->second.second.size() != 0)) { final_vertices_map.insert(*vs);  continue; } /*Always insert particles out of nowhere*/
0135         std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator  v2;
0136         for (v2 = final_vertices_map.begin(); v2 != final_vertices_map.end(); ++v2) if (vs->second.first == v2->second.first) {v2->second.second.insert(vs->second.second.begin(), vs->second.second.end()); break;}
0137         if (v2 == final_vertices_map.end()) final_vertices_map.insert(*vs);
0138     }
0139
0140     std::vector<GenParticlePtr> final_particles;
0141     std::set<int> used;
0142     for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >:: iterator it = final_vertices_map.begin(); it != final_vertices_map.end(); ++it)
0143     {
0144         GenVertexPtr v = it->first;
0145         std::set<int> in = it->second.first;
0146         std::set<int> out = it->second.second;
0147         used.insert(in.begin(), in.end());
0148         used.insert(out.begin(), out.end());
0149         for (std::set<int>::iterator el = in.begin(); el != in.end(); ++el) v->add_particle_in(particles_index[*el]);
0150         if (in.size() !=0 ) for (std::set<int>::iterator el = out.begin(); el != out.end(); ++el) v->add_particle_out(particles_index[*el]);
0151     }
0152     for (std::set<int>::iterator el = used.begin(); el != used.end(); ++el) final_particles.push_back(particles_index[*el]);
0153     /* One can put here a check on the number of particles/vertices*/
0154
0155     evt->add_tree(final_particles);
0156
0157     return true;
0158 }
0159 /** @brief  Converts GenEvent into HEPEVT. */
0160 template <class T>
0161 bool GenEvent_to_HEPEVT_nonstatic(const GenEvent* evt, T* A)
0162 {
0163     /// This writes an event out to the HEPEVT common block. The daughters
0164     /// field is NOT filled, because it is possible to contruct graphs
0165     /// for which the mothers and daughters cannot both be make sequential.
0166     /// This is consistent with how pythia fills HEPEVT (daughters are not
0167     /// necessarily filled properly) and how IO_HEPEVT reads HEPEVT.
0168     if ( !evt ) return false;
0169
0170     /*AV Sorting the vertices by the lengths of their longest incoming paths assures the mothers will not go before the daughters*/
0171     /* Calculate all paths*/
0172     std::map<ConstGenVertexPtr, int> longest_paths;
0173     for (ConstGenVertexPtr v: evt->vertices()) calculate_longest_path_to_top(v, longest_paths);
0174     /* Sort paths*/
0175     std::vector<std::pair<ConstGenVertexPtr, int> > sorted_paths;
0176     std::copy(longest_paths.begin(), longest_paths.end(), std::back_inserter(sorted_paths));
0177     std::sort(sorted_paths.begin(), sorted_paths.end(), pair_GenVertexPtr_int_greater());
0178
0179     std::vector<ConstGenParticlePtr> sorted_particles;
0180     std::vector<ConstGenParticlePtr> stable_particles;
0181     /*For a valid "Trust mothers" HEPEVT record we must  keep mothers together*/
0182     for (std::pair<ConstGenVertexPtr, int> it: sorted_paths)
0183     {
0184         std::vector<ConstGenParticlePtr> Q = it.first->particles_in();
0185         std::sort(Q.begin(), Q.end(), GenParticlePtr_greater());
0186         std::copy(Q.begin(), Q.end(), std::back_inserter(sorted_particles));
0187         /*For each vertex put all outgoing particles w/o end vertex. Ordering of particles to produces reproduceable record*/
0188         for (ConstGenParticlePtr pp: it.first->particles_out())
0189             if (!(pp->end_vertex())) stable_particles.push_back(pp);
0190     }
0191     std::sort(stable_particles.begin(), stable_particles.end(), GenParticlePtr_greater());
0192     std::copy(stable_particles.begin(), stable_particles.end(), std::back_inserter(sorted_particles));
0193
0194     int particle_counter;
0195     particle_counter = std::min(int(sorted_particles.size()), A->max_number_entries());
0196     /* fill the HEPEVT event record (MD code)*/
0197     A->set_event_number(evt->event_number());
0198     A->set_number_entries(particle_counter);
0199     for ( int i = 1; i <= particle_counter; ++i )
0200     {
0201         A->set_status(i, sorted_particles[i-1]->status());
0202         A->set_id(i, sorted_particles[i-1]->pid());
0203         FourVector m = sorted_particles[i-1]->momentum();
0204         A->set_momentum(i, m.px(), m.py(), m.pz(), m.e());
0205         A->set_mass(i, sorted_particles[i-1]->generated_mass());
0206         if ( sorted_particles[i-1]->production_vertex()  &&
0207                 sorted_particles[i-1]->production_vertex()->particles_in().size())
0208         {
0209             FourVector p = sorted_particles[i-1]->production_vertex()->position();
0210             A->set_position(i, p.x(), p.y(), p.z(), p.t() );
0211             std::vector<int> mothers;
0212             mothers.clear();
0213
0214             for (ConstGenParticlePtr it: sorted_particles[i-1]->production_vertex()->particles_in())
0215                 for ( int j = 1; j <= particle_counter; ++j )
0216                     if (sorted_particles[j-1] == (it))
0217                         mothers.push_back(j);
0218             std::sort(mothers.begin(), mothers.end());
0219             if (mothers.size() == 0)
0220                 mothers.push_back(0);
0221             if (mothers.size() == 1) mothers.push_back(mothers[0]);
0222
0223             A->set_parents(i, mothers.front(), mothers.back());
0224         }
0225         else
0226         {
0227             A->set_position(i, 0, 0, 0, 0);
0228             A->set_parents(i, 0, 0);
0229         }
0230         A->set_children(i, 0, 0);
0231     }
0232     return true;
0233 }
0234
0235
0236 /** @brief  Converts  HEPEVT into GenEvent. */
0237 template <class T>
0238 bool HEPEVT_to_GenEvent_static(GenEvent* evt)
0239 {
0240     if ( !evt ) { std::cerr << "HEPEVT_to_GenEvent_static  - passed null event." << std::endl; return false;}
0241     evt->set_event_number(T::event_number());
0242     std::map<GenParticlePtr, int > hepevt_particles;
0243     std::map<int, GenParticlePtr > particles_index;
0244     std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > hepevt_vertices;
0245     std::map<int, GenVertexPtr > vertex_index;
0246     int ne=0;
0247     ne=T::number_entries();
0248     for ( int i = 1; i <= ne; i++ )
0249     {
0250         GenParticlePtr p = std::make_shared<GenParticle>();
0251         p->set_momentum(FourVector(T::px(i), T::py(i), T::pz(i), T::e(i)));
0252         p->set_status(T::status(i));
0253         p->set_pid(T::id(i)); //Confusing!
0254         p->set_generated_mass(T::m(i));
0255         hepevt_particles[p] = i;
0256         particles_index[i] = p;
0257         GenVertexPtr v = std::make_shared<GenVertex>();
0258         v->set_position(FourVector(T::x(i), T::y(i), T::z(i), T::t(i)));
0259         v->add_particle_out(p);
0260         std::set<int> in;
0261         std::set<int> out;
0262         out.insert(i);
0263         vertex_index[i] = v;
0264         hepevt_vertices[v] = std::pair<std::set<int>, std::set<int> >(in, out);
0265     }
0266     /* The part above is always correct as it is a raw information without any topology.*/
0267
0268     /* In this way we trust mother information. The "Trust daughters" is not implemented.*/
0269     for (std::map<GenParticlePtr, int >::iterator it1 = hepevt_particles.begin(); it1 != hepevt_particles.end(); ++it1)
0270         for (std::map<GenParticlePtr, int >::iterator it2 = hepevt_particles.begin(); it2 != hepevt_particles.end(); ++it2) {
0271             if   (T::first_parent(it2->second) <= it1->second && it1->second <= T::last_parent(it2->second)) hepevt_vertices[it2->first->production_vertex()].first.insert(it1->second);
0272         }
0273     /* Now all incoming sets are correct for all vertices. But we have duplicates.*/
0274
0275     /* Disconnect all particles from the vertices*/
0276     for ( int i = 1; i <= T::number_entries(); i++ ) vertex_index[i]->remove_particle_out(particles_index[i]);
0277
0278     /*Fill container with vertices with unique sets of incoming particles. Merge the outcoming particle sets.*/
0279     std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > final_vertices_map;
0280     for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator vs = hepevt_vertices.begin(); vs != hepevt_vertices.end(); ++vs)
0281     {
0282         if ((final_vertices_map.size() == 0) || (vs->second.first.size() == 0 && vs->second.second.size() != 0)) { final_vertices_map.insert(*vs);  continue; } /*Always insert particles out of nowhere*/
0283         std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator  v2;
0284         for (v2 = final_vertices_map.begin(); v2 != final_vertices_map.end(); ++v2) if (vs->second.first == v2->second.first) {v2->second.second.insert(vs->second.second.begin(), vs->second.second.end()); break;}
0285         if (v2 == final_vertices_map.end()) final_vertices_map.insert(*vs);
0286     }
0287
0288     std::vector<GenParticlePtr> final_particles;
0289     std::set<int> used;
0290     for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >:: iterator it = final_vertices_map.begin(); it != final_vertices_map.end(); ++it)
0291     {
0292         GenVertexPtr v = it->first;
0293         std::set<int> in = it->second.first;
0294         std::set<int> out = it->second.second;
0295         used.insert(in.begin(), in.end());
0296         used.insert(out.begin(), out.end());
0297         for (const auto&  el: in) v->add_particle_in(particles_index[el]);
0298         if (in.size() !=0 ) for (const auto&  el: out) v->add_particle_out(particles_index[el]);
0299     }
0300     for (const auto&  el: used) final_particles.emplace_back(particles_index[el]);
0301     /* One can put here a check on the number of particles/vertices*/
0302
0303     evt->add_tree(final_particles);
0304
0305     return true;
0306 }
0307
0308 /** @brief  Converts GenEvent into HEPEVT. */
0309 template <class T>
0310 bool GenEvent_to_HEPEVT_static(const GenEvent* evt)
0311 {
0312     /// This writes an event out to the HEPEVT common block. The daughters
0313     /// field is NOT filled, because it is possible to contruct graphs
0314     /// for which the mothers and daughters cannot both be make sequential.
0315     /// This is consistent with how pythia fills HEPEVT (daughters are not
0316     /// necessarily filled properly) and how IO_HEPEVT reads HEPEVT.
0317     if ( !evt ) return false;
0318
0319     /*AV Sorting the vertices by the lengths of their longest incoming paths assures the mothers will not go before the daughters*/
0320     /* Calculate all paths*/
0321     std::map<ConstGenVertexPtr, int> longest_paths;
0322     for (ConstGenVertexPtr v: evt->vertices()) calculate_longest_path_to_top(v, longest_paths);
0323     /* Sort paths*/
0324     std::vector<std::pair<ConstGenVertexPtr, int> > sorted_paths;
0325     std::copy(longest_paths.begin(), longest_paths.end(), std::back_inserter(sorted_paths));
0326     std::sort(sorted_paths.begin(), sorted_paths.end(), pair_GenVertexPtr_int_greater());
0327
0328     std::vector<ConstGenParticlePtr> sorted_particles;
0329     std::vector<ConstGenParticlePtr> stable_particles;
0330     /*For a valid "Trust mothers" HEPEVT record we must  keep mothers together*/
0331     for (std::pair<ConstGenVertexPtr, int> it: sorted_paths)
0332     {
0333         std::vector<ConstGenParticlePtr> Q = it.first->particles_in();
0334         std::sort(Q.begin(), Q.end(), GenParticlePtr_greater());
0335         std::copy(Q.begin(), Q.end(), std::back_inserter(sorted_particles));
0336         /*For each vertex put all outgoing particles w/o end vertex. Ordering of particles to produces reproduceable record*/
0337         for (ConstGenParticlePtr pp: it.first->particles_out())
0338             if (!(pp->end_vertex())) stable_particles.push_back(pp);
0339     }
0340     std::sort(stable_particles.begin(), stable_particles.end(), GenParticlePtr_greater());
0341     std::copy(stable_particles.begin(), stable_particles.end(), std::back_inserter(sorted_particles));
0342
0343     int particle_counter;
0344     particle_counter = std::min(int(sorted_particles.size()), T::max_number_entries());
0345     /* fill the HEPEVT event record (MD code)*/
0346     T::set_event_number(evt->event_number());
0347     T::set_number_entries(particle_counter);
0348     for ( int i = 1; i <= particle_counter; ++i )
0349     {
0350         T::set_status(i, sorted_particles[i-1]->status());
0351         T::set_id(i, sorted_particles[i-1]->pid());
0352         FourVector m = sorted_particles[i-1]->momentum();
0353         T::set_momentum(i, m.px(), m.py(), m.pz(), m.e());
0354         T::set_mass(i, sorted_particles[i-1]->generated_mass());
0355         if ( sorted_particles[i-1]->production_vertex()  &&
0356                 sorted_particles[i-1]->production_vertex()->particles_in().size())
0357         {
0358             FourVector p = sorted_particles[i-1]->production_vertex()->position();
0359             T::set_position(i, p.x(), p.y(), p.z(), p.t() );
0360             std::vector<int> mothers;
0361             mothers.clear();
0362
0363             for (ConstGenParticlePtr it: sorted_particles[i-1]->production_vertex()->particles_in())
0364                 for ( int j = 1; j <= particle_counter; ++j )
0365                     if (sorted_particles[j-1] == (it))
0366                         mothers.push_back(j);
0367             std::sort(mothers.begin(), mothers.end());
0368             if (mothers.size() == 0)
0369                 mothers.push_back(0);
0370             if (mothers.size() == 1) mothers.push_back(mothers[0]);
0371
0372             T::set_parents(i, mothers.front(), mothers.back());
0373         }
0374         else
0375         {
0376             T::set_position(i, 0, 0, 0, 0);
0377             T::set_parents(i, 0, 0);
0378         }
0379         T::set_children(i, 0, 0);
0380     }
0381     return true;
0382 }
0383
0384 }
0385 #endif