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;
115 for (
int i = 1; i <= HEPEVT_Wrapper::number_entries(); i++ )
117 GenParticlePtr p=make_shared<GenParticle>();
118 p->set_momentum(
FourVector( HEPEVT_Wrapper::px(i), HEPEVT_Wrapper::py(i), HEPEVT_Wrapper::pz(i), HEPEVT_Wrapper::e(i) ));
119 p->set_status(HEPEVT_Wrapper::status(i));
120 p->set_pid(HEPEVT_Wrapper::id(i));
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>();
125 v->set_position(
FourVector( HEPEVT_Wrapper::x(i), HEPEVT_Wrapper::y(i), HEPEVT_Wrapper::z(i), HEPEVT_Wrapper::t(i)));
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)
138 if (HEPEVT_Wrapper::first_parent(it2->second)<=it1->second&&it1->second<=HEPEVT_Wrapper::last_parent(it2->second))
139 hepevt_vertices[it2->first->production_vertex()].first.insert(it1->second);
143 for (
int i = 1; i <= HEPEVT_Wrapper::number_entries(); i++ ) vertex_index[i]->remove_particle_out(particles_index[i]);
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));
210 int particle_counter=std::min(
int(sorted_particles.size()),HEPEVT_Wrapper::max_number_entries());
212 HEPEVT_Wrapper::set_event_number( evt->
event_number() );
213 HEPEVT_Wrapper::set_number_entries( particle_counter );
214 for (
int i = 1; i <= particle_counter; ++i )
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();
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]);