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
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Generic 4-vector.
Definition: FourVector.h:36
double e() const
Energy component of momentum.
Definition: FourVector.h:131
double pz() const
z-component of momentum
Definition: FourVector.h:124
double t() const
Time component of position/displacement.
Definition: FourVector.h:102
double px() const
x-component of momentum
Definition: FourVector.h:110
double py() const
y-component of momentum
Definition: FourVector.h:117
double x() const
x-component of position/displacement
Definition: FourVector.h:81
double y() const
y-component of position/displacement
Definition: FourVector.h:88
double z() const
z-component of position/displacement
Definition: FourVector.h:95
Stores event-related information.
Definition: GenEvent.h:41
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:265
int event_number() const
Get event number.
Definition: GenEvent.h:148
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:150
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
HepMC3 main namespace.
bool GenEvent_to_HEPEVT_static(const GenEvent *evt)
Converts GenEvent into HEPEVT.
bool HEPEVT_to_GenEvent_static(GenEvent *evt)
Converts HEPEVT into GenEvent.
bool GenEvent_to_HEPEVT_nonstatic(const GenEvent *evt, T *A)
Converts GenEvent into HEPEVT.
void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map< ConstGenVertexPtr, int > &pathl)
Calculates the path to the top (beam) particles.
bool HEPEVT_to_GenEvent_nonstatic(GenEvent *evt, T *A)
Converts HEPEVT into GenEvent.
comparison of two particles
bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
comparison of two particles
C structure representing Fortran common block HEPEVT T. Sjöstrand et al., "A proposed standard event ...
int * jdahep
Pointer to position of 1nd and 2nd (or last!) daughter.
int * isthep
Pointer to Status code.
int * nevhep
Pointer to Event number.
int * idhep
Pointer to PDG ID.
momentum_type * vhep
Pointer to time-space position: x, y, z, t.
momentum_type * phep
Pointer to momentum: px, py, pz, e, m.
int * jmohep
Pointer to position of 1st and 2nd (or last!) mother.
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 ...
int jmohep[max_particles][2]
Position of 1st and 2nd (or last!) mother.
momentum_type vhep[max_particles][4]
Time-space position: x, y, z, t.
momentum_type phep[max_particles][5]
Momentum: px, py, pz, e, m.
int idhep[max_particles]
PDG ID.
int jdahep[max_particles][2]
Position of 1nd and 2nd (or last!) daughter.
int isthep[max_particles]
Status code.
int nevhep
Event number.
int nhep
Number of entries in the event.
Order vertices with equal paths.
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....