HepMC3 event record library
HEPEVT_Helpers.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2021 The HepMC collaboration (see AUTHORS for details)
5 //
6 #ifndef HEPEVT_HELPERS_H
7 #define HEPEVT_HELPERS_H
8 /**
9  * @file HEPEVT_Helpers.h
10  * @brief Helper functions used to manipulate with HEPEVT block
11  *
12  */
13 #if defined(WIN32)&&(!defined(HEPMC3_NO_EXPORTS))
14 #if defined(HepMC3_EXPORTS)
15 #define HEPMC3_EXPORT_API __declspec(dllexport)
16 #else
17 #define HEPMC3_EXPORT_API __declspec(dllimport)
18 #endif
19 #else
20 #define HEPMC3_EXPORT_API
21 #endif
22 #include <algorithm>
23 #include <map>
24 
25 #include "HepMC3/GenEvent.h"
26 #include "HepMC3/GenParticle.h"
27 #include "HepMC3/GenVertex.h"
28 namespace HepMC3
29 {
30 
31 /** @struct HEPEVT_Templated
32  * @brief C structure representing Fortran common block HEPEVT
33  * T. Sjöstrand et al., "A proposed standard event record",
34  * in `Z physics at LEP 1', eds. G. Altarelli, R. Kleiss and C. Verzegnassi,
35  * Geneva, Switzerland, September 4-5, 1989, CERN 89-08 (Geneva, 1989), Vol. 3, p. 327
36  * Disk representation is given by Fortran WRITE/READ format.
37  */
38 template <int max_particles, typename momentum_type = double>
40 {
41  int nevhep; //!< Event number
42  int nhep; //!< Number of entries in the event
43  int isthep[max_particles]; //!< Status code
44  int idhep [max_particles]; //!< PDG ID
45  int jmohep[max_particles][2]; //!< Position of 1st and 2nd (or last!) mother
46  int jdahep[max_particles][2]; //!< Position of 1nd and 2nd (or last!) daughter
47  momentum_type phep [max_particles][5]; //!< Momentum: px, py, pz, e, m
48  momentum_type vhep [max_particles][4]; //!< Time-space position: x, y, z, t
49 };
50 
51 /** @struct HEPEVT_Pointers
52  * @brief C structure representing Fortran common block HEPEVT
53  * T. Sjöstrand et al., "A proposed standard event record",
54  * in `Z physics at LEP 1', eds. G. Altarelli, R. Kleiss and C. Verzegnassi,
55  * Geneva, Switzerland, September 4-5, 1989, CERN 89-08 (Geneva, 1989), Vol. 3, p. 327
56  * Disk representation is given by Fortran WRITE/READ format.
57  */
58 template<typename momentum_type = double>
60 {
61  int* nevhep; //!< Pointer to Event number
62  int* nhep; //!< Pointer to Number of entries in the event
63  int* isthep; //!< Pointer to Status code
64  int* idhep; //!< Pointer to PDG ID
65  int* jmohep; //!< Pointer to position of 1st and 2nd (or last!) mother
66  int* jdahep; //!< Pointer to position of 1nd and 2nd (or last!) daughter
67  momentum_type* phep; //!< Pointer to momentum: px, py, pz, e, m
68  momentum_type* vhep; //!< Pointer to time-space position: x, y, z, t
69 };
70 
71 /** @brief comparison of two particles */
73 {
74  /** @brief comparison of two particles */
75  bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const;
76 };
77 /** @brief Order vertices with equal paths. */
79 {
80  /** @brief Order vertices with equal paths. If the paths are equal, order in other quantities.
81  * We cannot use id, as it can be assigned in different way*/
82  bool operator()(const std::pair<ConstGenVertexPtr, int>& lx, const std::pair<ConstGenVertexPtr, int>& rx) const;
83 };
84 /** @brief Calculates the path to the top (beam) particles */
85 void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map<ConstGenVertexPtr, int>& pathl);
86 
87 
88 /** @brief Converts HEPEVT into GenEvent. */
89 template <class T>
91 {
92  if ( !evt ) { std::cerr << "HEPEVT_to_GenEvent_nonstatic - passed null event." << std::endl; return false;}
93  evt->set_event_number(A->event_number());
94  std::map<GenParticlePtr, int > hepevt_particles;
95  std::map<int, GenParticlePtr > particles_index;
96  std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > hepevt_vertices;
97  std::map<int, GenVertexPtr > vertex_index;
98  int ne=0;
99  ne=A->number_entries();
100  for ( int i = 1; i <= ne; i++ )
101  {
102  GenParticlePtr p = std::make_shared<GenParticle>();
103  p->set_momentum(FourVector(A->px(i), A->py(i), A->pz(i), A->e(i)));
104  p->set_status(A->status(i));
105  p->set_pid(A->id(i)); //Confusing!
106  p->set_generated_mass(A->m(i));
107  hepevt_particles[p] = i;
108  particles_index[i] = p;
109  GenVertexPtr v = std::make_shared<GenVertex>();
110  v->set_position(FourVector(A->x(i), A->y(i), A->z(i), A->t(i)));
111  v->add_particle_out(p);
112  std::set<int> in;
113  std::set<int> out;
114  out.insert(i);
115  vertex_index[i] = v;
116  hepevt_vertices[v] = std::pair<std::set<int>, std::set<int> >(in, out);
117  }
118  /* The part above is always correct as it is a raw information without any topology.*/
119 
120  /* In this way we trust mother information. The "Trust daughters" is not implemented.*/
121  for (std::map<GenParticlePtr, int >::iterator it1 = hepevt_particles.begin(); it1 != hepevt_particles.end(); ++it1)
122  for (std::map<GenParticlePtr, int >::iterator it2 = hepevt_particles.begin(); it2 != hepevt_particles.end(); ++it2) {
123  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);
124  }
125  /* Now all incoming sets are correct for all vertices. But we have duplicates.*/
126 
127  /* Disconnect all particles from the vertices*/
128  for ( int i = 1; i <= A->number_entries(); i++ ) vertex_index[i]->remove_particle_out(particles_index[i]);
129 
130  /*Fill container with vertices with unique sets of incoming particles. Merge the outcoming particle sets.*/
131  std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > final_vertices_map;
132  for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator vs = hepevt_vertices.begin(); vs != hepevt_vertices.end(); ++vs)
133  {
134  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*/
135  std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator v2;
136  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;}
137  if (v2 == final_vertices_map.end()) final_vertices_map.insert(*vs);
138  }
139 
140  std::vector<GenParticlePtr> final_particles;
141  std::set<int> used;
142  for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >:: iterator it = final_vertices_map.begin(); it != final_vertices_map.end(); ++it)
143  {
144  GenVertexPtr v = it->first;
145  std::set<int> in = it->second.first;
146  std::set<int> out = it->second.second;
147  used.insert(in.begin(), in.end());
148  used.insert(out.begin(), out.end());
149  for (std::set<int>::iterator el = in.begin(); el != in.end(); ++el) v->add_particle_in(particles_index[*el]);
150  if (in.size() !=0 ) for (std::set<int>::iterator el = out.begin(); el != out.end(); ++el) v->add_particle_out(particles_index[*el]);
151  }
152  for (std::set<int>::iterator el = used.begin(); el != used.end(); ++el) final_particles.push_back(particles_index[*el]);
153  /* One can put here a check on the number of particles/vertices*/
154 
155  evt->add_tree(final_particles);
156 
157  return true;
158 }
159 /** @brief Converts GenEvent into HEPEVT. */
160 template <class T>
162 {
163  /// This writes an event out to the HEPEVT common block. The daughters
164  /// field is NOT filled, because it is possible to contruct graphs
165  /// for which the mothers and daughters cannot both be make sequential.
166  /// This is consistent with how pythia fills HEPEVT (daughters are not
167  /// necessarily filled properly) and how IO_HEPEVT reads HEPEVT.
168  if ( !evt ) return false;
169 
170  /*AV Sorting the vertices by the lengths of their longest incoming paths assures the mothers will not go before the daughters*/
171  /* Calculate all paths*/
172  std::map<ConstGenVertexPtr, int> longest_paths;
173  for (ConstGenVertexPtr v: evt->vertices()) calculate_longest_path_to_top(v, longest_paths);
174  /* Sort paths*/
175  std::vector<std::pair<ConstGenVertexPtr, int> > sorted_paths;
176  std::copy(longest_paths.begin(), longest_paths.end(), std::back_inserter(sorted_paths));
177  std::sort(sorted_paths.begin(), sorted_paths.end(), pair_GenVertexPtr_int_greater());
178 
179  std::vector<ConstGenParticlePtr> sorted_particles;
180  std::vector<ConstGenParticlePtr> stable_particles;
181  /*For a valid "Trust mothers" HEPEVT record we must keep mothers together*/
182  for (std::pair<ConstGenVertexPtr, int> it: sorted_paths)
183  {
184  std::vector<ConstGenParticlePtr> Q = it.first->particles_in();
185  std::sort(Q.begin(), Q.end(), GenParticlePtr_greater());
186  std::copy(Q.begin(), Q.end(), std::back_inserter(sorted_particles));
187  /*For each vertex put all outgoing particles w/o end vertex. Ordering of particles to produces reproduceable record*/
188  for (ConstGenParticlePtr pp: it.first->particles_out())
189  if (!(pp->end_vertex())) stable_particles.push_back(pp);
190  }
191  std::sort(stable_particles.begin(), stable_particles.end(), GenParticlePtr_greater());
192  std::copy(stable_particles.begin(), stable_particles.end(), std::back_inserter(sorted_particles));
193 
194  int particle_counter;
195  particle_counter = std::min(int(sorted_particles.size()), A->max_number_entries());
196  /* fill the HEPEVT event record (MD code)*/
197  A->set_event_number(evt->event_number());
198  A->set_number_entries(particle_counter);
199  for ( int i = 1; i <= particle_counter; ++i )
200  {
201  A->set_status(i, sorted_particles[i-1]->status());
202  A->set_id(i, sorted_particles[i-1]->pid());
203  FourVector m = sorted_particles[i-1]->momentum();
204  A->set_momentum(i, m.px(), m.py(), m.pz(), m.e());
205  A->set_mass(i, sorted_particles[i-1]->generated_mass());
206  if ( sorted_particles[i-1]->production_vertex() &&
207  sorted_particles[i-1]->production_vertex()->particles_in().size())
208  {
209  FourVector p = sorted_particles[i-1]->production_vertex()->position();
210  A->set_position(i, p.x(), p.y(), p.z(), p.t() );
211  std::vector<int> mothers;
212  mothers.clear();
213 
214  for (ConstGenParticlePtr it: sorted_particles[i-1]->production_vertex()->particles_in())
215  for ( int j = 1; j <= particle_counter; ++j )
216  if (sorted_particles[j-1] == (it))
217  mothers.push_back(j);
218  std::sort(mothers.begin(), mothers.end());
219  if (mothers.size() == 0)
220  mothers.push_back(0);
221  if (mothers.size() == 1) mothers.push_back(mothers[0]);
222 
223  A->set_parents(i, mothers.front(), mothers.back());
224  }
225  else
226  {
227  A->set_position(i, 0, 0, 0, 0);
228  A->set_parents(i, 0, 0);
229  }
230  A->set_children(i, 0, 0);
231  }
232  return true;
233 }
234 
235 
236 /** @brief Converts HEPEVT into GenEvent. */
237 template <class T>
239 {
240  if ( !evt ) { std::cerr << "HEPEVT_to_GenEvent_static - passed null event." << std::endl; return false;}
241  evt->set_event_number(T::event_number());
242  std::map<GenParticlePtr, int > hepevt_particles;
243  std::map<int, GenParticlePtr > particles_index;
244  std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > hepevt_vertices;
245  std::map<int, GenVertexPtr > vertex_index;
246  int ne=0;
247  ne=T::number_entries();
248  for ( int i = 1; i <= ne; i++ )
249  {
250  GenParticlePtr p = std::make_shared<GenParticle>();
251  p->set_momentum(FourVector(T::px(i), T::py(i), T::pz(i), T::e(i)));
252  p->set_status(T::status(i));
253  p->set_pid(T::id(i)); //Confusing!
254  p->set_generated_mass(T::m(i));
255  hepevt_particles[p] = i;
256  particles_index[i] = p;
257  GenVertexPtr v = std::make_shared<GenVertex>();
258  v->set_position(FourVector(T::x(i), T::y(i), T::z(i), T::t(i)));
259  v->add_particle_out(p);
260  std::set<int> in;
261  std::set<int> out;
262  out.insert(i);
263  vertex_index[i] = v;
264  hepevt_vertices[v] = std::pair<std::set<int>, std::set<int> >(in, out);
265  }
266  /* The part above is always correct as it is a raw information without any topology.*/
267 
268  /* In this way we trust mother information. The "Trust daughters" is not implemented.*/
269  for (std::map<GenParticlePtr, int >::iterator it1 = hepevt_particles.begin(); it1 != hepevt_particles.end(); ++it1)
270  for (std::map<GenParticlePtr, int >::iterator it2 = hepevt_particles.begin(); it2 != hepevt_particles.end(); ++it2) {
271  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);
272  }
273  /* Now all incoming sets are correct for all vertices. But we have duplicates.*/
274 
275  /* Disconnect all particles from the vertices*/
276  for ( int i = 1; i <= T::number_entries(); i++ ) vertex_index[i]->remove_particle_out(particles_index[i]);
277 
278  /*Fill container with vertices with unique sets of incoming particles. Merge the outcoming particle sets.*/
279  std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > final_vertices_map;
280  for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator vs = hepevt_vertices.begin(); vs != hepevt_vertices.end(); ++vs)
281  {
282  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*/
283  std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator v2;
284  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;}
285  if (v2 == final_vertices_map.end()) final_vertices_map.insert(*vs);
286  }
287 
288  std::vector<GenParticlePtr> final_particles;
289  std::set<int> used;
290  for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >:: iterator it = final_vertices_map.begin(); it != final_vertices_map.end(); ++it)
291  {
292  GenVertexPtr v = it->first;
293  std::set<int> in = it->second.first;
294  std::set<int> out = it->second.second;
295  used.insert(in.begin(), in.end());
296  used.insert(out.begin(), out.end());
297  for (std::set<int>::iterator el = in.begin(); el != in.end(); ++el) v->add_particle_in(particles_index[*el]);
298  if (in.size() !=0 ) for (std::set<int>::iterator el = out.begin(); el != out.end(); ++el) v->add_particle_out(particles_index[*el]);
299  }
300  for (std::set<int>::iterator el = used.begin(); el != used.end(); ++el) final_particles.push_back(particles_index[*el]);
301  /* One can put here a check on the number of particles/vertices*/
302 
303  evt->add_tree(final_particles);
304 
305  return true;
306 }
307 
308 /** @brief Converts GenEvent into HEPEVT. */
309 template <class T>
311 {
312  /// This writes an event out to the HEPEVT common block. The daughters
313  /// field is NOT filled, because it is possible to contruct graphs
314  /// for which the mothers and daughters cannot both be make sequential.
315  /// This is consistent with how pythia fills HEPEVT (daughters are not
316  /// necessarily filled properly) and how IO_HEPEVT reads HEPEVT.
317  if ( !evt ) return false;
318 
319  /*AV Sorting the vertices by the lengths of their longest incoming paths assures the mothers will not go before the daughters*/
320  /* Calculate all paths*/
321  std::map<ConstGenVertexPtr, int> longest_paths;
322  for (ConstGenVertexPtr v: evt->vertices()) calculate_longest_path_to_top(v, longest_paths);
323  /* Sort paths*/
324  std::vector<std::pair<ConstGenVertexPtr, int> > sorted_paths;
325  std::copy(longest_paths.begin(), longest_paths.end(), std::back_inserter(sorted_paths));
326  std::sort(sorted_paths.begin(), sorted_paths.end(), pair_GenVertexPtr_int_greater());
327 
328  std::vector<ConstGenParticlePtr> sorted_particles;
329  std::vector<ConstGenParticlePtr> stable_particles;
330  /*For a valid "Trust mothers" HEPEVT record we must keep mothers together*/
331  for (std::pair<ConstGenVertexPtr, int> it: sorted_paths)
332  {
333  std::vector<ConstGenParticlePtr> Q = it.first->particles_in();
334  std::sort(Q.begin(), Q.end(), GenParticlePtr_greater());
335  std::copy(Q.begin(), Q.end(), std::back_inserter(sorted_particles));
336  /*For each vertex put all outgoing particles w/o end vertex. Ordering of particles to produces reproduceable record*/
337  for (ConstGenParticlePtr pp: it.first->particles_out())
338  if (!(pp->end_vertex())) stable_particles.push_back(pp);
339  }
340  std::sort(stable_particles.begin(), stable_particles.end(), GenParticlePtr_greater());
341  std::copy(stable_particles.begin(), stable_particles.end(), std::back_inserter(sorted_particles));
342 
343  int particle_counter;
344  particle_counter = std::min(int(sorted_particles.size()), T::max_number_entries());
345  /* fill the HEPEVT event record (MD code)*/
346  T::set_event_number(evt->event_number());
347  T::set_number_entries(particle_counter);
348  for ( int i = 1; i <= particle_counter; ++i )
349  {
350  T::set_status(i, sorted_particles[i-1]->status());
351  T::set_id(i, sorted_particles[i-1]->pid());
352  FourVector m = sorted_particles[i-1]->momentum();
353  T::set_momentum(i, m.px(), m.py(), m.pz(), m.e());
354  T::set_mass(i, sorted_particles[i-1]->generated_mass());
355  if ( sorted_particles[i-1]->production_vertex() &&
356  sorted_particles[i-1]->production_vertex()->particles_in().size())
357  {
358  FourVector p = sorted_particles[i-1]->production_vertex()->position();
359  T::set_position(i, p.x(), p.y(), p.z(), p.t() );
360  std::vector<int> mothers;
361  mothers.clear();
362 
363  for (ConstGenParticlePtr it: sorted_particles[i-1]->production_vertex()->particles_in())
364  for ( int j = 1; j <= particle_counter; ++j )
365  if (sorted_particles[j-1] == (it))
366  mothers.push_back(j);
367  std::sort(mothers.begin(), mothers.end());
368  if (mothers.size() == 0)
369  mothers.push_back(0);
370  if (mothers.size() == 1) mothers.push_back(mothers[0]);
371 
372  T::set_parents(i, mothers.front(), mothers.back());
373  }
374  else
375  {
376  T::set_position(i, 0, 0, 0, 0);
377  T::set_parents(i, 0, 0);
378  }
379  T::set_children(i, 0, 0);
380  }
381  return true;
382 }
383 
384 }
385 #endif
int event_number() const
Get event number.
Definition: GenEvent.h:148
int * isthep
Pointer to Status code.
Order vertices with equal paths.
int idhep[max_particles]
PDG ID.
int jdahep[max_particles][2]
Position of 1nd and 2nd (or last!) daughter.
HepMC3 main namespace.
int * nhep
Pointer to Number of entries in the event.
C structure representing Fortran common block HEPEVT T. Sjöstrand et al., "A proposed standard event...
comparison of two particles
double pz() const
z-component of momentum
Definition: FourVector.h:124
Definition of class GenParticle.
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
bool HEPEVT_to_GenEvent_nonstatic(GenEvent *evt, T *A)
Converts HEPEVT into GenEvent.
momentum_type vhep[max_particles][4]
Time-space position: x, y, z, t.
Definition of class GenVertex.
C structure representing Fortran common block HEPEVT T. Sjöstrand et al., "A proposed standard event...
bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
comparison of two particles
int nhep
Number of entries in the event.
momentum_type phep[max_particles][5]
Momentum: px, py, pz, e, m.
int jmohep[max_particles][2]
Position of 1st and 2nd (or last!) mother.
int * nevhep
Pointer to Event number.
double x() const
x-component of position/displacement
Definition: FourVector.h:81
int isthep[max_particles]
Status code.
Stores event-related information.
Definition: GenEvent.h:41
Generic 4-vector.
Definition: FourVector.h:36
bool HEPEVT_to_GenEvent_static(GenEvent *evt)
Converts HEPEVT into GenEvent.
double y() const
y-component of position/displacement
Definition: FourVector.h:88
double t() const
Time component of position/displacement.
Definition: FourVector.h:102
void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map< ConstGenVertexPtr, int > &pathl)
Calculates the path to the top (beam) particles.
double py() const
y-component of momentum
Definition: FourVector.h:117
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:265
momentum_type * phep
Pointer to momentum: px, py, pz, e, m.
bool operator()(const std::pair< ConstGenVertexPtr, int > &lx, const std::pair< ConstGenVertexPtr, int > &rx) const
Order vertices with equal paths. If the paths are equal, order in other quantities. We cannot use id, as it can be assigned in different way.
int * jdahep
Pointer to position of 1nd and 2nd (or last!) daughter.
int nevhep
Event number.
double e() const
Energy component of momentum.
Definition: FourVector.h:131
momentum_type * vhep
Pointer to time-space position: x, y, z, t.
bool GenEvent_to_HEPEVT_nonstatic(const GenEvent *evt, T *A)
Converts GenEvent into HEPEVT.
Definition of class GenEvent.
int * jmohep
Pointer to position of 1st and 2nd (or last!) mother.
double px() const
x-component of momentum
Definition: FourVector.h:110
int * idhep
Pointer to PDG ID.
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:150
bool GenEvent_to_HEPEVT_static(const GenEvent *evt)
Converts GenEvent into HEPEVT.
double z() const
z-component of position/displacement
Definition: FourVector.h:95