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