Line data    Source code

       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           0 :     const std::pair<HepMC::GenParticle*, HepMC::GenParticle*> beams = ge->beam_particles();
       8             : 
       9             :     // Make list of non-physical particle entries
      10           0 :     std::vector<HepMC::GenParticle*> unphys_particles;
      11           0 :     for (HepMC::GenEvent::particle_const_iterator pi = ge->particles_begin(); 
      12           0 :          pi != ge->particles_end(); ++pi) {
      13             :       // Beam particles might not have status = 4, but we want them anyway
      14           0 :       if (beams.first == *pi || beams.second == *pi) continue;
      15             :       // Filter by status
      16           0 :       const int status = (*pi)->status();
      17           0 :       if (status != 1 && status != 2 && status != 4) {
      18           0 :         unphys_particles.push_back(*pi);
      19             :       }      
      20           0 :     }
      21             : 
      22             :     // Remove each unphysical particle from the list
      23           0 :     while (unphys_particles.size()) {
      24           0 :       HepMC::GenParticle* gp = unphys_particles.back();
      25             : 
      26             :       // Get start and end vertices
      27           0 :       HepMC::GenVertex* vstart = gp->production_vertex();
      28           0 :       HepMC::GenVertex* vend = gp->end_vertex();
      29             : 
      30           0 :       if (vend == vstart) {
      31             :         // Deal with loops
      32           0 :         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           0 :         if (vend && vend->particles_out_size()) {
      38           0 :           std::vector<HepMC::GenParticle*> end_particles;
      39           0 :           for (HepMC::GenVertex::particles_out_const_iterator gpe = vend->particles_out_const_begin(); 
      40           0 :                gpe != vend->particles_out_const_end(); ++gpe) {
      41           0 :             end_particles.push_back(*gpe);
      42             :           }
      43             :           // Reset production vertices of child particles to bypass unphysical particle
      44           0 :           for (std::vector<HepMC::GenParticle*>::const_iterator gpe = end_particles.begin();
      45           0 :                gpe != end_particles.end(); ++gpe) {
      46             :             //std::cout << vstart << ", " << vend << std::endl;
      47           0 :             if (vstart) vstart->add_particle_out(*gpe);
      48             :           }
      49           0 :         } else {
      50             :           // If null end_vertex... stable unphysical particle?
      51             :         }
      52             : 
      53             :         // Delete unphysical particle and its orphaned end vertex
      54           0 :         delete vend;
      55           0 :         if (vstart) {
      56           0 :           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           0 :       unphys_particles.pop_back();
      65             :       //std::cout << unphys_particles.size() << std::endl;
      66             :     }
      67             : 
      68             :     // Delete any orphaned vertices
      69           0 :     std::vector<HepMC::GenVertex*> orphaned_vtxs;
      70           0 :     for (HepMC::GenEvent::vertex_const_iterator vi = ge->vertices_begin(); 
      71           0 :          vi != ge->vertices_end(); ++vi) {
      72           0 :       if ((*vi)->particles_in_size() == 0 && (*vi)->particles_out_size() == 0) {
      73           0 :         orphaned_vtxs.push_back(*vi);
      74             :       }
      75             :     }
      76           0 :     while (orphaned_vtxs.size()) {
      77           0 :       delete orphaned_vtxs.back();
      78           0 :       orphaned_vtxs.pop_back();
      79             :     }
      80           0 :   }
      81             : 
      82             :