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
18#include "HepMC3/GenEvent.h"
19#include "HepMC3/GenParticle.h"
20#include "HepMC3/GenVertex.h"
21#include <algorithm>
22#include <set>
23#include <vector>
24namespace HepMC3
25{
26
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 */
92void 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=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=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 );
246 }
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:83
double x() const
x-component of position/displacement
Definition FourVector.h:62
double y() const
y-component of position/displacement
Definition FourVector.h:69
double z() const
z-component of position/displacement
Definition FourVector.h:76
Stores event-related information.
Definition GenEvent.h:42
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition GenEvent.cc:268
int event_number() const
Get event number.
Definition GenEvent.h:136
void set_event_number(const int &num)
Set event number.
Definition GenEvent.h:138
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition GenEvent.cc:44
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.
Definition ReaderGZ.h:28
void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map< ConstGenVertexPtr, int > &pathl)
Calculates the path to the top (beam) particles.
struct HEPEVT * hepevtptr
Pointer to external (e.g. in Pythia6) struct with HEPEVT.
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....