TEvtGen/HepMC/filterEvent.cc

   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