From 0d8ae7116527b336222550031de05634a4982d3b Mon Sep 17 00:00:00 2001
From: Andre Sailer <andre.philippe.sailer@cern.ch>
Date: Fri, 16 Feb 2018 10:50:38 +0100
Subject: [PATCH] DDG4::PrimaryHanlder: add property RejectPDGs

---
 DDG4/include/DDG4/Geant4PrimaryHandler.h |  5 +++++
 DDG4/python/DDG4.py                      |  1 +
 DDG4/src/Geant4InputHandling.cpp         | 14 ++++++++++----
 DDG4/src/Geant4PrimaryHandler.cpp        |  1 +
 4 files changed, 17 insertions(+), 4 deletions(-)

diff --git a/DDG4/include/DDG4/Geant4PrimaryHandler.h b/DDG4/include/DDG4/Geant4PrimaryHandler.h
index 53d13ae6c..1c43cc44b 100644
--- a/DDG4/include/DDG4/Geant4PrimaryHandler.h
+++ b/DDG4/include/DDG4/Geant4PrimaryHandler.h
@@ -39,6 +39,11 @@ namespace dd4hep {
       virtual ~Geant4PrimaryHandler();
       /// Event generation action callback
       virtual void operator()(G4Event* event);
+
+    public:
+      /// particles with these PDG IDs are not passed to geant for simulation
+      std::set<int> m_rejectPDGs={1,2,3,4,5,6,21,23,24};
+
     };
   }    // End namespace sim
 }      // End namespace dd4hep
diff --git a/DDG4/python/DDG4.py b/DDG4/python/DDG4.py
index a9a87471b..71ffbc349 100644
--- a/DDG4/python/DDG4.py
+++ b/DDG4/python/DDG4.py
@@ -704,6 +704,7 @@ class Geant4:
     # Finally generate Geant4 primaries
     if have_mctruth:
       gen = GeneratorAction(self.kernel(),"Geant4PrimaryHandler/PrimaryHandler")
+      gen.RejectPDGs="{1,2,3,4,5,6,21,23,24}"
       gen.enableUI()
       if output_level is not None:
         gen.OutputLevel = output_level
diff --git a/DDG4/src/Geant4InputHandling.cpp b/DDG4/src/Geant4InputHandling.cpp
index 7215c58e0..e0ae2dc4d 100644
--- a/DDG4/src/Geant4InputHandling.cpp
+++ b/DDG4/src/Geant4InputHandling.cpp
@@ -16,6 +16,7 @@
 #include "DDG4/Geant4Primary.h"
 #include "DDG4/Geant4Context.h"
 #include "DDG4/Geant4Action.h"
+#include "DDG4/Geant4PrimaryHandler.h"
 #include "CLHEP/Units/SystemOfUnits.h"
 #include "CLHEP/Units/PhysicalConstants.h"
 
@@ -358,6 +359,7 @@ static vector< pair<Geant4Particle*,G4PrimaryParticle*> >
 getRelevant(set<int>& visited,
             map<int,G4PrimaryParticle*>& prim,
             Geant4PrimaryInteraction::ParticleMap& pm,
+            const set<int>& rejectPDGs,
             const Geant4ParticleHandle p)
 {
   typedef vector< pair<Geant4Particle*,G4PrimaryParticle*> > Primaries;
@@ -381,7 +383,6 @@ getRelevant(set<int>& visited,
     double proper_time = fabs(dp->time-p->time) * me;
 
     // -- remove original --- if the pdgID is not known (generator "strings") or the particle is quark, gluon, Z, W, etc.
-    const std::set<int> rejectPDGs{1, 2, 3, 4, 5, 6, 21, 23, 24};
     if (p.definition() && rejectPDGs.count(abs(p->pdgID)) == 0) {
       map<int,G4PrimaryParticle*>::iterator ip4 = prim.find(p->id);
       G4PrimaryParticle* p4 = (ip4 == prim.end()) ? 0 : (*ip4).second;
@@ -392,7 +393,7 @@ getRelevant(set<int>& visited,
         Primaries daughters;
         for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i)  {
           if ( visited.find(*i) == visited.end() )  {
-            Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
+            Primaries tmp = getRelevant(visited,prim,pm,rejectPDGs,pm[*i]);
             daughters.insert(daughters.end(), tmp.begin(),tmp.end());
           }
         }
@@ -404,7 +405,7 @@ getRelevant(set<int>& visited,
     else  {
       for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i)  {
         if ( visited.find(*i) == visited.end() )  {
-          Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
+          Primaries tmp = getRelevant(visited,prim,pm,rejectPDGs,pm[*i]);
           res.insert(res.end(), tmp.begin(),tmp.end());
         }
       }
@@ -427,6 +428,11 @@ int dd4hep::sim::generatePrimaries(const Geant4Action* caller,
   map<int,G4PrimaryParticle*> prim;
   set<int> visited;
 
+  auto const* primHandler = dynamic_cast<const Geant4PrimaryHandler*>(caller);
+  auto const& rejectPDGs = primHandler ? primHandler->m_rejectPDGs : std::set<int>();
+
+  caller->debug("Rejecting PDGs: %s", [&rejectPDGs](){ std::stringstream str; for (int i: rejectPDGs) { str  << i << ", "; } return str;}().str().c_str());
+
   if ( interaction->locked )  {
     caller->abortRun("Locked interactions may not be used to generate primaries!",
                      "Cannot handle a native G4 primary record!");
@@ -449,7 +455,7 @@ int dd4hep::sim::generatePrimaries(const Geant4Action* caller,
 	    mask.set(G4PARTICLE_HAS_SECONDARIES);
 	  }
 	  if ( p->parents.size() == 0 )  {
-	    Primaries relevant = getRelevant(visited,prim,pm,p);
+	    Primaries relevant = getRelevant(visited,prim,pm,rejectPDGs,p);
 	    for(Primaries::const_iterator j=relevant.begin(); j!= relevant.end(); ++j)  {
 	      Geant4ParticleHandle r = (*j).first;
 	      G4PrimaryParticle* p4 = (*j).second;
diff --git a/DDG4/src/Geant4PrimaryHandler.cpp b/DDG4/src/Geant4PrimaryHandler.cpp
index 8b9cbe727..8c98a3e62 100644
--- a/DDG4/src/Geant4PrimaryHandler.cpp
+++ b/DDG4/src/Geant4PrimaryHandler.cpp
@@ -24,6 +24,7 @@ Geant4PrimaryHandler::Geant4PrimaryHandler(Geant4Context* ctxt, const std::strin
   : Geant4GeneratorAction(ctxt,nam)
 {
   InstanceCount::increment(this);
+  declareProperty("RejectPDGs", m_rejectPDGs);
 }
 
 /// Default destructor
-- 
GitLab