[571] | 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 |
|
---|