HepMC3 event record library
HEPEVT_Wrapper.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2020 The HepMC collaboration (see AUTHORS for details)
5 //
6 /**
7  * @file HEPEVT_Wrapper.cc
8  * @brief Implementation of conversion functions for HEPEVT block
9  ***********************************************************************
10  * Some parts from HepMC2 library
11  * Matt.Dobbs@Cern.CH, January 2000
12  * HEPEVT IO class
13  ***********************************************************************
14  *
15  */
16 #ifndef HEPEVT_WRAPPER_HEADER_ONLY
17 #include <algorithm>
18 #include <set>
19 #include <vector>
20 
21 #include "HepMC3/HEPEVT_Wrapper.h"
22 #include "HepMC3/GenEvent.h"
23 #include "HepMC3/GenParticle.h"
24 #include "HepMC3/GenVertex.h"
25 namespace HepMC3
26 {
27 
28 struct HEPEVT* hepevtptr;
29 
30 
31 /** @brief comparison of two particles */
33 {
34  /** @brief comparison of two particles */
35  bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
36  {
37  /* Cannot use id as it could be different.*/
38  if (lx->pid() != rx->pid()) return (lx->pid() < rx->pid());
39  if (lx->status() != rx->status()) return (lx->status() < rx->status());
40  /*Hopefully it will reach this point not too often.*/
41  return (lx->momentum().e() < rx->momentum().e());
42  }
43 };
44 /** @brief Order vertices with equal paths. */
46 {
47  /** @brief Order vertices with equal paths. If the paths are equal, order in other quantities.
48  * We cannot use id, as it can be assigned in different way*/
49  bool operator()(const std::pair<ConstGenVertexPtr, int>& lx, const std::pair<ConstGenVertexPtr, int>& rx) const
50  {
51  if (lx.second != rx.second) return (lx.second < rx.second);
52  if (lx.first->particles_in().size() != rx.first->particles_in().size()) return (lx.first->particles_in().size() < rx.first->particles_in().size());
53  if (lx.first->particles_out().size() != rx.first->particles_out().size()) return (lx.first->particles_out().size() < rx.first->particles_out().size());
54  /* The code below is usefull mainly for debug. Assures strong ordering.*/
55  std::vector<int> lx_id_in;
56  std::vector<int> rx_id_in;
57  for (ConstGenParticlePtr pp: lx.first->particles_in()) lx_id_in.push_back(pp->pid());
58  for (ConstGenParticlePtr pp: rx.first->particles_in()) rx_id_in.push_back(pp->pid());
59  std::sort(lx_id_in.begin(), lx_id_in.end());
60  std::sort(rx_id_in.begin(), rx_id_in.end());
61  for (unsigned int i = 0; i < lx_id_in.size(); i++) if (lx_id_in[i] != rx_id_in[i]) return (lx_id_in[i] < rx_id_in[i]);
62 
63  std::vector<int> lx_id_out;
64  std::vector<int> rx_id_out;
65  for (ConstGenParticlePtr pp: lx.first->particles_in()) lx_id_out.push_back(pp->pid());
66  for (ConstGenParticlePtr pp: rx.first->particles_in()) rx_id_out.push_back(pp->pid());
67  std::sort(lx_id_out.begin(), lx_id_out.end());
68  std::sort(rx_id_out.begin(), rx_id_out.end());
69  for (unsigned int i = 0; i < lx_id_out.size(); i++) if (lx_id_out[i] != rx_id_out[i]) return (lx_id_out[i] < rx_id_out[i]);
70 
71  std::vector<double> lx_mom_in;
72  std::vector<double> rx_mom_in;
73  for (ConstGenParticlePtr pp: lx.first->particles_in()) lx_mom_in.push_back(pp->momentum().e());
74  for (ConstGenParticlePtr pp: rx.first->particles_in()) rx_mom_in.push_back(pp->momentum().e());
75  std::sort(lx_mom_in.begin(), lx_mom_in.end());
76  std::sort(rx_mom_in.begin(), rx_mom_in.end());
77  for (unsigned int i = 0; i < lx_mom_in.size(); i++) if (lx_mom_in[i] != rx_mom_in[i]) return (lx_mom_in[i] < rx_mom_in[i]);
78 
79  std::vector<double> lx_mom_out;
80  std::vector<double> rx_mom_out;
81  for (ConstGenParticlePtr pp: lx.first->particles_in()) lx_mom_out.push_back(pp->momentum().e());
82  for (ConstGenParticlePtr pp: rx.first->particles_in()) rx_mom_out.push_back(pp->momentum().e());
83  std::sort(lx_mom_out.begin(), lx_mom_out.end());
84  std::sort(rx_mom_out.begin(), rx_mom_out.end());
85  for (unsigned int i = 0; i < lx_mom_out.size(); i++) if (lx_mom_out[i] != rx_mom_out[i]) return (lx_mom_out[i] < rx_mom_out[i]);
86  /* The code above is usefull mainly for debug*/
87 
88  return (lx.first < lx.first); /*This is random. This should never happen*/
89  }
90 };
91 /** @brief Calculates the path to the top (beam) particles */
92 void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map<ConstGenVertexPtr, int>& pathl)
93 {
94  int p = 0;
95  for (ConstGenParticlePtr pp: v->particles_in()) {
96  ConstGenVertexPtr v2 = pp->production_vertex();
97  if (v2 == v) continue; //LOOP! THIS SHOULD NEVER HAPPEN FOR A PROPER EVENT!
98  if (!v2) p = std::max(p, 1);
99  else
100  {if (pathl.find(v2) == pathl.end()) calculate_longest_path_to_top(v2, pathl); p = std::max(p, pathl[v2]+1);}
101  }
102  pathl[v] = p;
103  return;
104 }
105 
106 
108 {
109  if ( !evt ) { std::cerr << "IO_HEPEVT::fill_next_event error - passed null event." << std::endl; return false;}
111  std::map<GenParticlePtr, int > hepevt_particles;
112  std::map<int, GenParticlePtr > particles_index;
113  std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > hepevt_vertices;
114  std::map<int, GenVertexPtr > vertex_index;
115  for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); i++ )
116  {
117  GenParticlePtr p = std::make_shared<GenParticle>();
119  p->set_status(HEPEVT_Wrapper::status(i));
120  p->set_pid(HEPEVT_Wrapper::id(i)); //Confusing!
121  p->set_generated_mass(HEPEVT_Wrapper::m(i));
122  hepevt_particles[p] = i;
123  particles_index[i] = p;
124  GenVertexPtr v = std::make_shared<GenVertex>();
126  v->add_particle_out(p);
127  std::set<int> in;
128  std::set<int> out;
129  out.insert(i);
130  vertex_index[i] = v;
131  hepevt_vertices[v] = std::pair<std::set<int>, std::set<int> >(in, out);
132  }
133  /* The part above is always correct as it is a raw information without any topology.*/
134 
135  /* In this way we trust mother information. The "Trust daughters" is not implemented.*/
136  for (std::map<GenParticlePtr, int >::iterator it1 = hepevt_particles.begin(); it1 != hepevt_particles.end(); ++it1)
137  for (std::map<GenParticlePtr, int >::iterator it2 = hepevt_particles.begin(); it2 != hepevt_particles.end(); ++it2)
138  if (HEPEVT_Wrapper::first_parent(it2->second) <= it1->second && it1->second <= HEPEVT_Wrapper::last_parent(it2->second)) /*I'm you father, Luck!*/
139  hepevt_vertices[it2->first->production_vertex()].first.insert(it1->second);
140  /* Now all incoming sets are correct for all vertices. But we have duplicates.*/
141 
142  /* Disconnect all particles from the vertices*/
143  for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); i++ ) vertex_index[i]->remove_particle_out(particles_index[i]);
144 
145  /*Fill container with vertices with unique sets of incoming particles. Merge the outcoming particle sets.*/
146  std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > final_vertices_map;
147  for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator vs = hepevt_vertices.begin(); vs != hepevt_vertices.end(); ++vs)
148  {
149  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*/
150  std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator v2;
151  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;}
152  if (v2 == final_vertices_map.end()) final_vertices_map.insert(*vs);
153  }
154 
155  std::vector<GenParticlePtr> final_particles;
156  std::set<int> used;
157  for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >:: iterator it = final_vertices_map.begin(); it != final_vertices_map.end(); ++it)
158  {
159  GenVertexPtr v = it->first;
160  std::set<int> in = it->second.first;
161  std::set<int> out = it->second.second;
162  used.insert(in.begin(), in.end());
163  used.insert(out.begin(), out.end());
164  for (std::set<int>::iterator el = in.begin(); el != in.end(); ++el) v->add_particle_in(particles_index[*el]);
165  if (in.size() !=0 ) for (std::set<int>::iterator el = out.begin(); el != out.end(); ++el) v->add_particle_out(particles_index[*el]);
166  }
167  for (std::set<int>::iterator el = used.begin(); el != used.end(); ++el) final_particles.push_back(particles_index[*el]);
168  /* One can put here a check on the number of particles/vertices*/
169 
170  evt->add_tree(final_particles);
171 
172  return true;
173 }
174 
175 
177 {
178  /// This writes an event out to the HEPEVT common block. The daughters
179  /// field is NOT filled, because it is possible to contruct graphs
180  /// for which the mothers and daughters cannot both be make sequential.
181  /// This is consistent with how pythia fills HEPEVT (daughters are not
182  /// necessarily filled properly) and how IO_HEPEVT reads HEPEVT.
183  //
184  if ( !evt ) return false;
185 
186  /*AV Sorting the vertices by the lengths of their longest incoming paths assures the mothers will not go before the daughters*/
187  /* Calculate all paths*/
188  std::map<ConstGenVertexPtr, int> longest_paths;
189  for (ConstGenVertexPtr v: evt->vertices()) calculate_longest_path_to_top(v, longest_paths);
190  /* Sort paths*/
191  std::vector<std::pair<ConstGenVertexPtr, int> > sorted_paths;
192  copy(longest_paths.begin(), longest_paths.end(), std::back_inserter(sorted_paths));
193  sort(sorted_paths.begin(), sorted_paths.end(), pair_GenVertexPtr_int_greater());
194 
195  std::vector<ConstGenParticlePtr> sorted_particles;
196  std::vector<ConstGenParticlePtr> stable_particles;
197  /*For a valid "Trust mothers" HEPEVT record we must keep mothers together*/
198  for (std::pair<ConstGenVertexPtr, int> it: sorted_paths)
199  {
200  std::vector<ConstGenParticlePtr> Q = it.first->particles_in();
201  sort(Q.begin(), Q.end(), GenParticlePtr_greater_order());
202  copy(Q.begin(), Q.end(), std::back_inserter(sorted_particles));
203  /*For each vertex put all outgoing particles w/o end vertex. Ordering of particles to produces reproduceable record*/
204  for (ConstGenParticlePtr pp: it.first->particles_out())
205  if (!(pp->end_vertex())) stable_particles.push_back(pp);
206  }
207  sort(stable_particles.begin(), stable_particles.end(), GenParticlePtr_greater_order());
208  copy(stable_particles.begin(), stable_particles.end(), std::back_inserter(sorted_particles));
209 
210  int particle_counter = std::min(int(sorted_particles.size()), HEPEVT_Wrapper::max_number_entries());
211  /* fill the HEPEVT event record (MD code)*/
213  HEPEVT_Wrapper::set_number_entries(particle_counter);
214  for ( int i = 1; i <= particle_counter; ++i )
215  {
216  HEPEVT_Wrapper::set_status(i, sorted_particles[i-1]->status());
217  HEPEVT_Wrapper::set_id(i, sorted_particles[i-1]->pid());
218  FourVector m = sorted_particles[i-1]->momentum();
219  HEPEVT_Wrapper::set_momentum(i, m.px(), m.py(), m.pz(), m.e());
220  HEPEVT_Wrapper::set_mass(i, sorted_particles[i-1]->generated_mass());
221  // there should ALWAYS be particles in any vertex, but some generators
222  // are making non-kosher HepMC events
223  if ( sorted_particles[i-1]->production_vertex() &&
224  sorted_particles[i-1]->production_vertex()->particles_in().size())
225  {
226  FourVector p = sorted_particles[i-1]->production_vertex()->position();
227  HEPEVT_Wrapper::set_position(i, p.x(), p.y(), p.z(), p.t() );
228  std::vector<int> mothers;
229  mothers.clear();
230 
231  for (ConstGenParticlePtr it: sorted_particles[i-1]->production_vertex()->particles_in())
232  for ( int j = 1; j <= particle_counter; ++j )
233  if (sorted_particles[j-1] == (it))
234  mothers.push_back(j);
235  sort(mothers.begin(), mothers.end());
236  if (mothers.size() == 0)
237  mothers.push_back(0);
238  if (mothers.size() == 1) mothers.push_back(mothers[0]);
239 
240  HEPEVT_Wrapper::set_parents(i, mothers.front(), mothers.back());
241  }
242  else
243  {
244  HEPEVT_Wrapper::set_position(i, 0, 0, 0, 0);
246  }
248  }
249  return true;
250 }
251 }
252 #endif
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
static double m(const int &index)
Get generated mass.
Order vertices with equal paths.
static double z(const int &index)
Get Z Production vertex.
double t() const
Time component of position/displacement.
Definition: FourVector.h:102
static double y(const int &index)
Get Y Production vertex.
Definition of class GenParticle.
static double t(const int &index)
Get production time.
static int status(const int &index)
Get status code.
static int event_number()
Get event number.
static double x(const int &index)
Get X Production vertex.
Definition of class GenVertex.
static void set_mass(const int &index, double mass)
Set mass.
static double pz(const int &index)
Get Z momentum.
static void set_children(const int &index, const int &firstchild, const int &lastchild)
Set children.
static double e(const int &index)
Get Energy.
double z() const
z-component of position/displacement
Definition: FourVector.h:95
struct HEPEVT * hepevtptr
This is a pointer to HEPEVT common block.
double x() const
x-component of position/displacement
Definition: FourVector.h:81
int event_number() const
Get event number.
Definition: GenEvent.h:136
static void set_position(const int &index, const double &x, const double &y, const double &z, const double &t)
Set position in time-space.
static void set_parents(const int &index, const int &firstparent, const int &lastparent)
Set parents.
static void set_id(const int &index, const int &id)
Set PDG particle id.
double e() const
Energy component of momentum.
Definition: FourVector.h:131
Stores event-related information.
Definition: GenEvent.h:41
comparison of two particles
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.
Generic 4-vector.
Definition: FourVector.h:36
double px() const
x-component of momentum
Definition: FourVector.h:110
static double py(const int &index)
Get Y momentum.
Fortran common block HEPEVT.
static void set_number_entries(const int &noentries)
Set number of entries.
static bool HEPEVT_to_GenEvent(GenEvent *evt)
Convert HEPEVT to GenEvent.
void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map< ConstGenVertexPtr, int > &pathl)
Calculates the path to the top (beam) particles.
static int id(const int &index)
Get PDG particle id.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:265
static void set_event_number(const int &evtno)
Set event number.
static bool GenEvent_to_HEPEVT(const GenEvent *evt)
Convert GenEvent to HEPEVT.
double y() const
y-component of position/displacement
Definition: FourVector.h:88
bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
comparison of two particles
static int first_parent(const int &index)
Get index of 1st mother.
static double px(const int &index)
Get X momentum.
double pz() const
z-component of momentum
Definition: FourVector.h:124
static void set_status(const int &index, const int &status)
Set status code.
Definition of class GenEvent.
static int last_parent(const int &index)
Get index of last mother.
double py() const
y-component of momentum
Definition: FourVector.h:117
static void set_momentum(const int &index, const double &px, const double &py, const double &pz, const double &e)
Set 4-momentum.
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:138
Definition of class HEPEVT_Wrapper.
static int number_entries()
Get number of entries.
static int max_number_entries()
Block size.