16 #ifndef HEPEVT_WRAPPER_HEADER_ONLY
34 bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx )
const
37 if (lx->pid() !=rx->pid())
return (lx->pid() < rx->pid());
38 if (lx->status() !=rx->status())
return (lx->status() < rx->status());
40 return (lx->momentum().e()<rx->momentum().e());
49 bool operator()(
const std::pair<ConstGenVertexPtr,int>& lx,
const std::pair<ConstGenVertexPtr,int>& rx )
const
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());
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]);
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]);
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]);
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]);
88 return (lx.first<lx.first);
95 for(ConstGenParticlePtr pp: v->particles_in()) {
96 ConstGenVertexPtr v2 = pp->production_vertex();
98 if (!v2) p=std::max(p,1);
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;
117 GenParticlePtr p=std::make_shared<GenParticle>();
122 hepevt_particles[p]=i;
123 particles_index[i]=p;
124 GenVertexPtr v=std::make_shared<GenVertex>();
126 v->add_particle_out(p);
131 hepevt_vertices[v]=std::pair<std::set<int>,std::set<int> >(in,out);
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)
139 hepevt_vertices[it2->first->production_vertex()].first.insert(it1->second);
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)
149 if ((final_vertices_map.size()==0)||(vs->second.first.size()==0&&vs->second.second.size()!=0)) { final_vertices_map.insert(*vs);
continue; }
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);
155 std::vector<GenParticlePtr> final_particles;
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)
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]);
167 for (std::set<int>::iterator el=used.begin(); el!=used.end(); ++el) final_particles.push_back(particles_index[*el]);
184 if ( !evt )
return false;
188 std::map<ConstGenVertexPtr,int> longest_paths;
191 std::vector<std::pair<ConstGenVertexPtr,int> > sorted_paths;
192 copy(longest_paths.begin(),longest_paths.end(),std::back_inserter(sorted_paths));
195 std::vector<ConstGenParticlePtr> sorted_particles;
196 std::vector<ConstGenParticlePtr> stable_particles;
198 for (std::pair<ConstGenVertexPtr,int> it: sorted_paths)
200 std::vector<ConstGenParticlePtr> Q = it.first->particles_in();
202 copy(Q.begin(),Q.end(),std::back_inserter(sorted_particles));
204 for(ConstGenParticlePtr pp: it.first->particles_out())
205 if(!(pp->end_vertex())) stable_particles.push_back(pp);
208 copy(stable_particles.begin(),stable_particles.end(),std::back_inserter(sorted_particles));
214 for (
int i = 1; i <= particle_counter; ++i )
223 if ( sorted_particles[i-1]->production_vertex() &&
224 sorted_particles[i-1]->production_vertex()->particles_in().size())
226 FourVector p = sorted_particles[i-1]->production_vertex()->position();
228 std::vector<int> mothers;
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]);
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
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.
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
double x() const
x-component of position/displacement
int event_number() const
Get event number.
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.
Stores event-related information.
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.
double px() const
x-component of momentum
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.
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
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
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
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 of class HEPEVT_Wrapper.
static int number_entries()
Get number of entries.
static int max_number_entries()
Block size.