1 | // from Andy Buckley
|
---|
2 |
|
---|
3 | #include "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 |
|
---|