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