]>
Commit | Line | Data |
---|---|---|
0ca57c2f | 1 | // from Andy Buckley |
2 | ||
3 | #include "HepMC/GenEvent.h" | |
4 | ||
5 | void filterEvent(HepMC::GenEvent* ge) { | |
6 | // We treat beam particles a bit specially | |
7 | const std::pair<HepMC::GenParticle*, HepMC::GenParticle*> beams = ge->beam_particles(); | |
8 | ||
9 | // Make list of non-physical particle entries | |
10 | std::vector<HepMC::GenParticle*> unphys_particles; | |
11 | for (HepMC::GenEvent::particle_const_iterator pi = ge->particles_begin(); | |
12 | pi != ge->particles_end(); ++pi) { | |
13 | // Beam particles might not have status = 4, but we want them anyway | |
14 | if (beams.first == *pi || beams.second == *pi) continue; | |
15 | // Filter by status | |
16 | const int status = (*pi)->status(); | |
17 | if (status != 1 && status != 2 && status != 4) { | |
18 | unphys_particles.push_back(*pi); | |
19 | } | |
20 | } | |
21 | ||
22 | // Remove each unphysical particle from the list | |
23 | while (unphys_particles.size()) { | |
24 | HepMC::GenParticle* gp = unphys_particles.back(); | |
25 | ||
26 | // Get start and end vertices | |
27 | HepMC::GenVertex* vstart = gp->production_vertex(); | |
28 | HepMC::GenVertex* vend = gp->end_vertex(); | |
29 | ||
30 | if (vend == vstart) { | |
31 | // Deal with loops | |
32 | vstart->remove_particle(gp); | |
33 | } else { | |
34 | ||
35 | // Connect decay particles from end vertex to start vertex | |
36 | /// @todo Have to build a list, since the GV::add_particle_out method modifies the end vertex! | |
37 | if (vend && vend->particles_out_size()) { | |
38 | std::vector<HepMC::GenParticle*> end_particles; | |
39 | for (HepMC::GenVertex::particles_out_const_iterator gpe = vend->particles_out_const_begin(); | |
40 | gpe != vend->particles_out_const_end(); ++gpe) { | |
41 | end_particles.push_back(*gpe); | |
42 | } | |
43 | // Reset production vertices of child particles to bypass unphysical particle | |
44 | for (std::vector<HepMC::GenParticle*>::const_iterator gpe = end_particles.begin(); | |
45 | gpe != end_particles.end(); ++gpe) { | |
46 | //std::cout << vstart << ", " << vend << std::endl; | |
47 | if (vstart) vstart->add_particle_out(*gpe); | |
48 | } | |
49 | } else { | |
50 | // If null end_vertex... stable unphysical particle? | |
51 | } | |
52 | ||
53 | // Delete unphysical particle and its orphaned end vertex | |
54 | delete vend; | |
55 | if (vstart) { | |
56 | delete vstart->remove_particle(gp); | |
57 | }// else { | |
58 | /// @todo Why does this cause an error? | |
59 | // delete gp; | |
60 | //} | |
61 | } | |
62 | ||
63 | // Remove deleted particle from list | |
64 | unphys_particles.pop_back(); | |
65 | //std::cout << unphys_particles.size() << std::endl; | |
66 | } | |
67 | ||
68 | // Delete any orphaned vertices | |
69 | std::vector<HepMC::GenVertex*> orphaned_vtxs; | |
70 | for (HepMC::GenEvent::vertex_const_iterator vi = ge->vertices_begin(); | |
71 | vi != ge->vertices_end(); ++vi) { | |
72 | if ((*vi)->particles_in_size() == 0 && (*vi)->particles_out_size() == 0) { | |
73 | orphaned_vtxs.push_back(*vi); | |
74 | } | |
75 | } | |
76 | while (orphaned_vtxs.size()) { | |
77 | delete orphaned_vtxs.back(); | |
78 | orphaned_vtxs.pop_back(); | |
79 | } | |
80 | } | |
81 | ||
82 |