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]);
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Definition of class HEPEVT_Wrapper.
double t() const
Time component of position/displacement.
double x() const
x-component of position/displacement
double y() const
y-component of position/displacement
double z() const
z-component of position/displacement
Stores event-related information.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
int event_number() const
Get event number.
void set_event_number(const int &num)
Set event number.
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
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.
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....