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-2019 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 "HepMC3/HEPEVT_Wrapper.h"
18 #include "HepMC3/GenEvent.h"
19 #include "HepMC3/GenParticle.h"
20 #include "HepMC3/GenVertex.h"
21 #include <algorithm>
22 #include <set>
23 #include <vector>
24 namespace HepMC3
25 {
26 
27 struct HEPEVT* hepevtptr;
28 
29 
30 /** @brief comparison of two particles */
32 {
33  /** @brief comparison of two particles */
34  bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx ) const
35  {
36  /* Cannot use id as it could be different*/
37  if (lx->pid() !=rx->pid()) return (lx->pid() < rx->pid());
38  if (lx->status() !=rx->status()) return (lx->status() < rx->status());
39  /*Hopefully it will reach this point not too ofter.*/
40  return (lx->momentum().e()<rx->momentum().e());
41 
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 TODO: implement "Trust daughters"*/
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 );
245  HEPEVT_Wrapper::set_parents( i, 0, 0 );
246  }
247  HEPEVT_Wrapper::set_children( i, 0, 0 );
248  }
249  return true;
250 }
251 }
252 #endif
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Definition of class HEPEVT_Wrapper.
Generic 4-vector.
Definition: FourVector.h:35
double t() const
Time component of position/displacement.
Definition: FourVector.h:101
double x() const
x-component of position/displacement
Definition: FourVector.h:80
double y() const
y-component of position/displacement
Definition: FourVector.h:87
double z() const
z-component of position/displacement
Definition: FourVector.h:94
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:267
int event_number() const
Get event number.
Definition: GenEvent.h:135
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:137
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
static double pz(const int &index)
Get Z momentum.
static void set_mass(const int &index, double mass)
Set mass.
static double py(const int &index)
Get Y momentum.
static void set_momentum(const int &index, const double &px, const double &py, const double &pz, const double &e)
Set 4-momentum.
static void set_number_entries(const int &noentries)
Set number of entries.
static double m(const int &index)
Get generated mass.
static bool GenEvent_to_HEPEVT(const GenEvent *evt)
Convert GenEvent to HEPEVT.
static void set_id(const int &index, const int &id)
Set PDG particle id.
static int max_number_entries()
Block size.
static void set_status(const int &index, const int &status)
Set status code.
static double y(const int &index)
Get Y Production vertex.
static double t(const int &index)
Get production time.
static int last_parent(const int &index)
Get index of last mother.
static double e(const int &index)
Get Energy.
static double z(const int &index)
Get Z Production vertex.
static double x(const int &index)
Get X Production vertex.
static int event_number()
Get event number.
static int status(const int &index)
Get status code.
static bool HEPEVT_to_GenEvent(GenEvent *evt)
Convert HEPEVT to GenEvent.
static void set_parents(const int &index, const int &firstparent, const int &lastparent)
Set parents.
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 double px(const int &index)
Get X momentum.
static int first_parent(const int &index)
Get index of 1st mother.
static void set_children(const int &index, const int &firstchild, const int &lastchild)
Set children.
static int number_entries()
Get number of entries.
static int id(const int &index)
Get PDG particle id.
static void set_event_number(const int &evtno)
Set event number.
HepMC3 main namespace.
void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map< ConstGenVertexPtr, int > &pathl)
Calculates the path to the top (beam) particles.
Fortran common block HEPEVT.
comparison of two particles
bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
comparison of two particles
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....