6 #ifndef HEPEVT_HELPERS_H 7 #define HEPEVT_HELPERS_H 13 #if defined(WIN32)&&(!defined(HEPMC3_NO_EXPORTS)) 14 #if defined(HepMC3_EXPORTS) 15 #define HEPMC3_EXPORT_API __declspec(dllexport) 17 #define HEPMC3_EXPORT_API __declspec(dllimport) 20 #define HEPMC3_EXPORT_API 38 template <
int max_particles,
typename momentum_type =
double>
47 momentum_type
phep [max_particles][5];
48 momentum_type
vhep [max_particles][4];
58 template<
typename momentum_type =
double>
75 bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx)
const;
82 bool operator()(
const std::pair<ConstGenVertexPtr, int>& lx,
const std::pair<ConstGenVertexPtr, int>& rx)
const;
92 if ( !evt ) { std::cerr <<
"HEPEVT_to_GenEvent_nonstatic - passed null event." << std::endl;
return false;}
94 std::map<GenParticlePtr, int > hepevt_particles;
95 std::map<int, GenParticlePtr > particles_index;
96 std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > hepevt_vertices;
97 std::map<int, GenVertexPtr > vertex_index;
99 ne=A->number_entries();
100 for (
int i = 1; i <= ne; i++ )
102 GenParticlePtr p = std::make_shared<GenParticle>();
103 p->set_momentum(
FourVector(A->px(i), A->py(i), A->pz(i), A->e(i)));
104 p->set_status(A->status(i));
105 p->set_pid(A->id(i));
106 p->set_generated_mass(A->m(i));
107 hepevt_particles[p] = i;
108 particles_index[i] = p;
109 GenVertexPtr v = std::make_shared<GenVertex>();
110 v->set_position(
FourVector(A->x(i), A->y(i), A->z(i), A->t(i)));
111 v->add_particle_out(p);
116 hepevt_vertices[v] = std::pair<std::set<int>, std::set<int> >(in, out);
121 for (std::map<GenParticlePtr, int >::iterator it1 = hepevt_particles.begin(); it1 != hepevt_particles.end(); ++it1)
122 for (std::map<GenParticlePtr, int >::iterator it2 = hepevt_particles.begin(); it2 != hepevt_particles.end(); ++it2) {
123 if (A->first_parent(it2->second) <= it1->second && it1->second <= A->last_parent(it2->second)) hepevt_vertices[it2->first->production_vertex()].first.insert(it1->second);
128 for (
int i = 1; i <= A->number_entries(); i++ ) vertex_index[i]->remove_particle_out(particles_index[i]);
131 std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > final_vertices_map;
132 for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator vs = hepevt_vertices.begin(); vs != hepevt_vertices.end(); ++vs)
134 if ((final_vertices_map.size() == 0) || (vs->second.first.size() == 0 && vs->second.second.size() != 0)) { final_vertices_map.insert(*vs);
continue; }
135 std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator v2;
136 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;}
137 if (v2 == final_vertices_map.end()) final_vertices_map.insert(*vs);
140 std::vector<GenParticlePtr> final_particles;
142 for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >:: iterator it = final_vertices_map.begin(); it != final_vertices_map.end(); ++it)
144 GenVertexPtr v = it->first;
145 std::set<int> in = it->second.first;
146 std::set<int> out = it->second.second;
147 used.insert(in.begin(), in.end());
148 used.insert(out.begin(), out.end());
149 for (std::set<int>::iterator el = in.begin(); el != in.end(); ++el) v->add_particle_in(particles_index[*el]);
150 if (in.size() !=0 )
for (std::set<int>::iterator el = out.begin(); el != out.end(); ++el) v->add_particle_out(particles_index[*el]);
152 for (std::set<int>::iterator el = used.begin(); el != used.end(); ++el) final_particles.push_back(particles_index[*el]);
168 if ( !evt )
return false;
172 std::map<ConstGenVertexPtr, int> longest_paths;
175 std::vector<std::pair<ConstGenVertexPtr, int> > sorted_paths;
176 std::copy(longest_paths.begin(), longest_paths.end(), std::back_inserter(sorted_paths));
179 std::vector<ConstGenParticlePtr> sorted_particles;
180 std::vector<ConstGenParticlePtr> stable_particles;
182 for (std::pair<ConstGenVertexPtr, int> it: sorted_paths)
184 std::vector<ConstGenParticlePtr> Q = it.first->particles_in();
186 std::copy(Q.begin(), Q.end(), std::back_inserter(sorted_particles));
188 for (ConstGenParticlePtr pp: it.first->particles_out())
189 if (!(pp->end_vertex())) stable_particles.push_back(pp);
192 std::copy(stable_particles.begin(), stable_particles.end(), std::back_inserter(sorted_particles));
194 int particle_counter;
195 particle_counter = std::min(
int(sorted_particles.size()), A->max_number_entries());
198 A->set_number_entries(particle_counter);
199 for (
int i = 1; i <= particle_counter; ++i )
201 A->set_status(i, sorted_particles[i-1]->status());
202 A->set_id(i, sorted_particles[i-1]->pid());
203 FourVector m = sorted_particles[i-1]->momentum();
204 A->set_momentum(i, m.
px(), m.
py(), m.
pz(), m.
e());
205 A->set_mass(i, sorted_particles[i-1]->generated_mass());
206 if ( sorted_particles[i-1]->production_vertex() &&
207 sorted_particles[i-1]->production_vertex()->particles_in().size())
209 FourVector p = sorted_particles[i-1]->production_vertex()->position();
210 A->set_position(i, p.
x(), p.
y(), p.
z(), p.
t() );
211 std::vector<int> mothers;
214 for (ConstGenParticlePtr it: sorted_particles[i-1]->production_vertex()->particles_in())
215 for (
int j = 1; j <= particle_counter; ++j )
216 if (sorted_particles[j-1] == (it))
217 mothers.push_back(j);
218 std::sort(mothers.begin(), mothers.end());
219 if (mothers.size() == 0)
220 mothers.push_back(0);
221 if (mothers.size() == 1) mothers.push_back(mothers[0]);
223 A->set_parents(i, mothers.front(), mothers.back());
227 A->set_position(i, 0, 0, 0, 0);
228 A->set_parents(i, 0, 0);
230 A->set_children(i, 0, 0);
240 if ( !evt ) { std::cerr <<
"HEPEVT_to_GenEvent_static - passed null event." << std::endl;
return false;}
242 std::map<GenParticlePtr, int > hepevt_particles;
243 std::map<int, GenParticlePtr > particles_index;
244 std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > hepevt_vertices;
245 std::map<int, GenVertexPtr > vertex_index;
247 ne=T::number_entries();
248 for (
int i = 1; i <= ne; i++ )
250 GenParticlePtr p = std::make_shared<GenParticle>();
251 p->set_momentum(
FourVector(T::px(i), T::py(i), T::pz(i), T::e(i)));
252 p->set_status(T::status(i));
253 p->set_pid(T::id(i));
254 p->set_generated_mass(T::m(i));
255 hepevt_particles[p] = i;
256 particles_index[i] = p;
257 GenVertexPtr v = std::make_shared<GenVertex>();
258 v->set_position(
FourVector(T::x(i), T::y(i), T::z(i), T::t(i)));
259 v->add_particle_out(p);
264 hepevt_vertices[v] = std::pair<std::set<int>, std::set<int> >(in, out);
269 for (std::map<GenParticlePtr, int >::iterator it1 = hepevt_particles.begin(); it1 != hepevt_particles.end(); ++it1)
270 for (std::map<GenParticlePtr, int >::iterator it2 = hepevt_particles.begin(); it2 != hepevt_particles.end(); ++it2) {
271 if (T::first_parent(it2->second) <= it1->second && it1->second <= T::last_parent(it2->second)) hepevt_vertices[it2->first->production_vertex()].first.insert(it1->second);
276 for (
int i = 1; i <= T::number_entries(); i++ ) vertex_index[i]->remove_particle_out(particles_index[i]);
279 std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > > final_vertices_map;
280 for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator vs = hepevt_vertices.begin(); vs != hepevt_vertices.end(); ++vs)
282 if ((final_vertices_map.size() == 0) || (vs->second.first.size() == 0 && vs->second.second.size() != 0)) { final_vertices_map.insert(*vs);
continue; }
283 std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >::iterator v2;
284 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;}
285 if (v2 == final_vertices_map.end()) final_vertices_map.insert(*vs);
288 std::vector<GenParticlePtr> final_particles;
290 for (std::map<GenVertexPtr, std::pair<std::set<int>, std::set<int> > >:: iterator it = final_vertices_map.begin(); it != final_vertices_map.end(); ++it)
292 GenVertexPtr v = it->first;
293 std::set<int> in = it->second.first;
294 std::set<int> out = it->second.second;
295 used.insert(in.begin(), in.end());
296 used.insert(out.begin(), out.end());
297 for (std::set<int>::iterator el = in.begin(); el != in.end(); ++el) v->add_particle_in(particles_index[*el]);
298 if (in.size() !=0 )
for (std::set<int>::iterator el = out.begin(); el != out.end(); ++el) v->add_particle_out(particles_index[*el]);
300 for (std::set<int>::iterator el = used.begin(); el != used.end(); ++el) final_particles.push_back(particles_index[*el]);
317 if ( !evt )
return false;
321 std::map<ConstGenVertexPtr, int> longest_paths;
324 std::vector<std::pair<ConstGenVertexPtr, int> > sorted_paths;
325 std::copy(longest_paths.begin(), longest_paths.end(), std::back_inserter(sorted_paths));
328 std::vector<ConstGenParticlePtr> sorted_particles;
329 std::vector<ConstGenParticlePtr> stable_particles;
331 for (std::pair<ConstGenVertexPtr, int> it: sorted_paths)
333 std::vector<ConstGenParticlePtr> Q = it.first->particles_in();
335 std::copy(Q.begin(), Q.end(), std::back_inserter(sorted_particles));
337 for (ConstGenParticlePtr pp: it.first->particles_out())
338 if (!(pp->end_vertex())) stable_particles.push_back(pp);
341 std::copy(stable_particles.begin(), stable_particles.end(), std::back_inserter(sorted_particles));
343 int particle_counter;
344 particle_counter = std::min(
int(sorted_particles.size()), T::max_number_entries());
347 T::set_number_entries(particle_counter);
348 for (
int i = 1; i <= particle_counter; ++i )
350 T::set_status(i, sorted_particles[i-1]->status());
351 T::set_id(i, sorted_particles[i-1]->pid());
352 FourVector m = sorted_particles[i-1]->momentum();
353 T::set_momentum(i, m.
px(), m.
py(), m.
pz(), m.
e());
354 T::set_mass(i, sorted_particles[i-1]->generated_mass());
355 if ( sorted_particles[i-1]->production_vertex() &&
356 sorted_particles[i-1]->production_vertex()->particles_in().size())
358 FourVector p = sorted_particles[i-1]->production_vertex()->position();
359 T::set_position(i, p.
x(), p.
y(), p.
z(), p.
t() );
360 std::vector<int> mothers;
363 for (ConstGenParticlePtr it: sorted_particles[i-1]->production_vertex()->particles_in())
364 for (
int j = 1; j <= particle_counter; ++j )
365 if (sorted_particles[j-1] == (it))
366 mothers.push_back(j);
367 std::sort(mothers.begin(), mothers.end());
368 if (mothers.size() == 0)
369 mothers.push_back(0);
370 if (mothers.size() == 1) mothers.push_back(mothers[0]);
372 T::set_parents(i, mothers.front(), mothers.back());
376 T::set_position(i, 0, 0, 0, 0);
377 T::set_parents(i, 0, 0);
379 T::set_children(i, 0, 0);
int event_number() const
Get event number.
int * isthep
Pointer to Status code.
Order vertices with equal paths.
int idhep[max_particles]
PDG ID.
int jdahep[max_particles][2]
Position of 1nd and 2nd (or last!) daughter.
int * nhep
Pointer to Number of entries in the event.
C structure representing Fortran common block HEPEVT T. Sjöstrand et al., "A proposed standard event...
comparison of two particles
double pz() const
z-component of momentum
Definition of class GenParticle.
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
bool HEPEVT_to_GenEvent_nonstatic(GenEvent *evt, T *A)
Converts HEPEVT into GenEvent.
momentum_type vhep[max_particles][4]
Time-space position: x, y, z, t.
Definition of class GenVertex.
C structure representing Fortran common block HEPEVT T. Sjöstrand et al., "A proposed standard event...
bool operator()(ConstGenParticlePtr lx, ConstGenParticlePtr rx) const
comparison of two particles
int nhep
Number of entries in the event.
momentum_type phep[max_particles][5]
Momentum: px, py, pz, e, m.
int jmohep[max_particles][2]
Position of 1st and 2nd (or last!) mother.
int * nevhep
Pointer to Event number.
double x() const
x-component of position/displacement
int isthep[max_particles]
Status code.
Stores event-related information.
bool HEPEVT_to_GenEvent_static(GenEvent *evt)
Converts HEPEVT into GenEvent.
double y() const
y-component of position/displacement
double t() const
Time component of position/displacement.
void calculate_longest_path_to_top(ConstGenVertexPtr v, std::map< ConstGenVertexPtr, int > &pathl)
Calculates the path to the top (beam) particles.
double py() const
y-component of momentum
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
momentum_type * phep
Pointer to momentum: px, py, pz, e, m.
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.
int * jdahep
Pointer to position of 1nd and 2nd (or last!) daughter.
double e() const
Energy component of momentum.
momentum_type * vhep
Pointer to time-space position: x, y, z, t.
bool GenEvent_to_HEPEVT_nonstatic(const GenEvent *evt, T *A)
Converts GenEvent into HEPEVT.
Definition of class GenEvent.
int * jmohep
Pointer to position of 1st and 2nd (or last!) mother.
double px() const
x-component of momentum
int * idhep
Pointer to PDG ID.
void set_event_number(const int &num)
Set event number.
bool GenEvent_to_HEPEVT_static(const GenEvent *evt)
Converts GenEvent into HEPEVT.
double z() const
z-component of position/displacement