From af974ddb24188eb07705ca7eb26587294fa754ed Mon Sep 17 00:00:00 2001
From: Markus Frank <markus.frank@cern.ch>
Date: Wed, 24 Sep 2014 16:35:09 +0000
Subject: [PATCH] More modular user particle handling with example
 Geant4TCUserParticleHandler

---
 DDG4/examples/CLICSidSimu.py                  |  9 +-
 DDG4/include/DDG4/Geant4Particle.h            |  3 +-
 DDG4/include/DDG4/Geant4ParticleHandler.h     |  7 +-
 DDG4/include/DDG4/Geant4UserParticleHandler.h | 42 +++++++-
 DDG4/plugins/Geant4TCUserParticleHandler.cpp  | 98 +++++++++++++++++++
 DDG4/python/DDG4.py                           | 14 +--
 DDG4/src/Geant4ParticleHandler.cpp            | 60 +++++++-----
 DDG4/src/Geant4UserParticleHandler.cpp        |  4 +
 8 files changed, 193 insertions(+), 44 deletions(-)
 create mode 100644 DDG4/plugins/Geant4TCUserParticleHandler.cpp

diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py
index ef58298f4..d37fd24dc 100644
--- a/DDG4/examples/CLICSidSimu.py
+++ b/DDG4/examples/CLICSidSimu.py
@@ -109,6 +109,7 @@ def run():
   gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/bbbb_3TeV.stdhep"
   #gen.Input = "LCIOFileReader|/home/frankm/SW/data/mcparticles_pi-_5GeV.slcio"
   gen.OutputLevel = 4 # generator_output_level
+  gen.MomentumScale = 0.1
   gen.Mask = 1
   kernel.generatorAction().adopt(gen)
   # Install vertex smearing for this interaction
@@ -153,12 +154,12 @@ def run():
   #part.SaveProcesses = ['conv','Decay']
   part.SaveProcesses = ['Decay']
   part.MinimalKineticEnergy = 100*MeV
-  part.Z_trackers = 165.70*cm  # EcalEndcap_zmin
-  part.R_trackers = 126.50*cm  # EcalBarrel_rmin
   part.OutputLevel = 5 # generator_output_level
   part.enableUI()
-  #user = DDG4.Action(kernel,"Geant4UserParticleHandler/UserParticleHandler")
-  #part.adopt(user)
+  user = DDG4.Action(kernel,"Geant4TCUserParticleHandler/UserParticleHandler")
+  user.TrackingVolume_Zmax = 165.70*cm  # EcalEndcap_zmin
+  user.TrackingVolume_Rmax = 126.50*cm  # EcalBarrel_rmin
+  part.adopt(user)
 
   """
   rdr = DDG4.GeneratorAction(kernel,"LcioGeneratorAction/Reader")
diff --git a/DDG4/include/DDG4/Geant4Particle.h b/DDG4/include/DDG4/Geant4Particle.h
index 0dc6a874d..dc7d7de81 100644
--- a/DDG4/include/DDG4/Geant4Particle.h
+++ b/DDG4/include/DDG4/Geant4Particle.h
@@ -54,7 +54,8 @@ namespace DD4hep {
       G4PARTICLE_KEEP_PARENT = 1<<6,
       G4PARTICLE_CREATED_CALORIMETER_HIT = 1<<7,
       G4PARTICLE_CREATED_TRACKER_HIT = 1<<8,
-      G4PARTICLE_KEEP_ALWAYS = 1<<9,
+      G4PARTICLE_KEEP_USER = 1<<9,
+      G4PARTICLE_KEEP_ALWAYS = 1<<10,
 
       // Generator status for a given particles: bit 0...3 come from LCIO, rest is internal
       G4PARTICLE_GEN_EMPTY           = 1<<0,  // Empty line
diff --git a/DDG4/include/DDG4/Geant4ParticleHandler.h b/DDG4/include/DDG4/Geant4ParticleHandler.h
index ff291107d..f28154e9e 100644
--- a/DDG4/include/DDG4/Geant4ParticleHandler.h
+++ b/DDG4/include/DDG4/Geant4ParticleHandler.h
@@ -75,11 +75,8 @@ namespace DD4hep {
       /// Property: Energy cut below which particles are not collected, but assigned to the parent
       double m_kinEnergyCut;
 
-      double m_zTracker, m_rTracker;
-
-
       /// Global particle identifier. Obtained at the begin of the event.
-      int m_globalParticleID, m_initialParticleID;
+      int m_globalParticleID;
       /// User action pointer
       Geant4UserParticleHandler* m_userHandler;
       /// Primary map
@@ -113,7 +110,7 @@ namespace DD4hep {
       /// Default destructor
       virtual ~Geant4ParticleHandler();
       /// Adopt the user particle handler
-      bool adopt(Geant4UserParticleHandler* action);
+      bool adopt(Geant4Action* action);
       /// Event generation action callback
       virtual void operator()(G4Event* event);
       /// User stepping callback
diff --git a/DDG4/include/DDG4/Geant4UserParticleHandler.h b/DDG4/include/DDG4/Geant4UserParticleHandler.h
index 745cb3de8..3ae3f01eb 100644
--- a/DDG4/include/DDG4/Geant4UserParticleHandler.h
+++ b/DDG4/include/DDG4/Geant4UserParticleHandler.h
@@ -39,7 +39,7 @@ namespace DD4hep {
      *  Clients may inherit from this class and override the approriate methods
      *  to add additional information in form of a DataExtension object to the Particle.
      *
-     *  The dfault implementation is always empty!
+     *  The default implementation is always empty!
      *
      * @author  M.Frank
      * @version 1.0
@@ -58,13 +58,51 @@ namespace DD4hep {
       virtual void begin(const G4Event* event);
       /// Post-event action callback
       virtual void end(const G4Event* event);
+
       /// User stepping callback
+      /** Allow the user to intercept particle handling in the pre track action.
+       *  The default implementation is empty.
+       *
+       *  Note: The particle passed is a temporary and will be copied if kept.
+       */
       virtual void step(const G4Step* step, G4SteppingManager* mgr, Particle& particle);
-      /// Pre-track action callback
+
+      /// Pre-track action callback.
+      /** Allow the user to intercept particle handling in the pre track action.
+       *  e.g. attach relevant user information.
+       *  The default implementation is empty.
+       *
+       *  Note: The particle passed is a temporary and will be copied if kept.
+       */
       virtual void begin(const G4Track* track, Particle& particle);
+
       /// Post-track action callback
+      /** Allow the user to force the particle handling in the post track action
+       *  set the reason mask to NULL in order to drop the particle.
+       *  The parent's reasoning mask will be or'ed with the particle's mask
+       *  to preserve the MC truth for the hit creation.
+       *  The default implementation is empty.
+       *
+       *  Note: The particle passed is a temporary and will be copied if kept.
+       */
       virtual void end(const G4Track* track, Particle& particle);
+
+      /// Callback to be answered if the particle MUST be kept during recombination step
+      /** Allow the user to force the particle handling either by
+       *  or the reason mask with G4PARTICLE_KEEP_USER or
+       *  to set the reason mask to NULL in order to drop it.
+       *  The default implementation is empty.
+       *
+       *  Note: This may override all other decisions!
+       *        Default implementation is empty.
+       */
+      virtual void keepParticle(Particle& particle);
+
       /// Callback when parent should be combined
+      /** Called before a particle is removed from the final record.
+       *  Relevant particle properties of the parent may be updated.
+       *  The default implementation is empty.
+       */
       virtual void combine(Particle& to_be_deleted, Particle& remaining_parent);
     };
   }    // End namespace Simulation
diff --git a/DDG4/plugins/Geant4TCUserParticleHandler.cpp b/DDG4/plugins/Geant4TCUserParticleHandler.cpp
new file mode 100644
index 000000000..efa025a0f
--- /dev/null
+++ b/DDG4/plugins/Geant4TCUserParticleHandler.cpp
@@ -0,0 +1,98 @@
+// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+//  AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+//  Author     : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_GEANT4TCUSERPARTICLEHANDLER_H
+#define DD4HEP_DDG4_GEANT4TCUSERPARTICLEHANDLER_H
+
+// Framework include files
+#include "DDG4/Geant4UserParticleHandler.h"
+
+/*
+ *   DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+  /*
+   *   Simulation namespace declaration
+   */
+  namespace Simulation {
+
+    /** Geant4TCUserParticleHandler
+     *
+     *  Rejects to keep particles, which are created outside a tracking cylinder.
+     *
+     *  TC stands for TrackingCylinder ;-)
+     *
+     * @author  M.Frank
+     * @version 1.0
+     */
+    class Geant4TCUserParticleHandler : public Geant4UserParticleHandler  {
+      double m_zTracker, m_rTracker;
+    public:
+      /// Standard constructor
+      Geant4TCUserParticleHandler(Geant4Context* context, const std::string& nam);
+
+      /// Default destructor
+      virtual ~Geant4TCUserParticleHandler() {}
+
+      /// Post-track action callback
+      /** Allow the user to force the particle handling in the post track action
+       *  set the reason mask to NULL in order to drop the particle.
+       *  The parent's reasoning mask will be or'ed with the particle's mask
+       *  to preserve the MC truth for the hit creation.
+       *  The default implementation is empty.
+       *
+       *  Note: The particle passed is a temporary and will be copied if kept.
+       */
+      virtual void end(const G4Track* track, Particle& particle);
+
+    };
+  }    // End namespace Simulation
+}      // End namespace DD4hep
+
+#endif // DD4HEP_DDG4_GEANT4TCUSERPARTICLEHANDLER_H
+
+// $Id: Geant4Field.cpp 888 2013-11-14 15:54:56Z markus.frank@cern.ch $
+//====================================================================
+//  AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+//  Author     : M.Frank
+//
+//====================================================================
+// Framework include files
+//#include "DDG4/Geant4TCUserParticleHandler.h"
+#include "DDG4/Geant4Particle.h"
+#include "DDG4/Factories.h"
+
+using namespace DD4hep::Simulation;
+DECLARE_GEANT4ACTION(Geant4TCUserParticleHandler)
+
+/// Standard constructor
+Geant4TCUserParticleHandler::Geant4TCUserParticleHandler(Geant4Context* context, const std::string& nam)
+: Geant4UserParticleHandler(context,nam) 
+{
+  declareProperty("TrackingVolume_Zmax",m_zTracker=1e100);
+  declareProperty("TrackingVolume_Rmax",m_rTracker=1e100);
+}
+
+/// Post-track action callback
+void Geant4TCUserParticleHandler::end(const G4Track* /* track */, Particle& p)  {
+  double r_end  = sqrt(p.vex*p.vex + p.vey*p.vey);
+  double z_end  = fabs(p.vez);
+  // Keep particles with end-vertex inside tracking (including back-scattered particles)
+  if ( r_end <= m_rTracker && z_end <= m_zTracker )  {
+    return;
+  }
+  double r_prod = sqrt(p.vsx*p.vsx + p.vsy*p.vsy);
+  double z_prod = fabs(p.vsz);
+  // Otherwise: reject particles created outside the tracking volume.
+  if ( r_prod > m_rTracker || z_prod > m_zTracker )  {
+    p.reason = 0;
+  }
+}
diff --git a/DDG4/python/DDG4.py b/DDG4/python/DDG4.py
index f863a9222..b1e2a3a36 100644
--- a/DDG4/python/DDG4.py
+++ b/DDG4/python/DDG4.py
@@ -225,16 +225,16 @@ class Simple:
     self.kernel.generatorAction().add(gun)
     return gun
 
-  def setupUI(typ='csh',vis=False,ui=True):
+  def setupUI(self,typ='csh',vis=False,ui=True):
     # Configure UI
-    ui = Action(self.kernel,"Geant4UIManager/UI")
-    ui.HaveVIS = vis
-    ui.HaveUI = ui
-    ui.SessionType = typ
-    self.kernel.registerGlobalAction(ui)
+    ui_action = Action(self.kernel,"Geant4UIManager/UI")
+    ui_action.HaveVIS = vis
+    ui_action.HaveUI = ui
+    ui_action.SessionType = typ
+    self.kernel.registerGlobalAction(ui_action)
 
   def setupCshUI(self,typ='csh',vis=False,ui=True):
-    self.setupUI('csh',vis,ui)
+    self.setupUI(typ='csh',vis=vis,ui=ui)
 
   def setupROOTOutput(self,name,output):
     evt_root = EventAction(self.kernel,'Geant4Output2ROOT/'+name)
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index 1b7f1c3b0..8e9309789 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -69,8 +69,6 @@ Geant4ParticleHandler::Geant4ParticleHandler(Geant4Context* context, const strin
   declareProperty("KeepAllParticles",    m_keepAll = false);
   declareProperty("SaveProcesses",       m_processNames);
   declareProperty("MinimalKineticEnergy",m_kinEnergyCut = 100e0*MeV);
-  declareProperty("Z_trackers",m_zTracker=1e100);
-  declareProperty("R_trackers",m_rTracker=1e100);
   InstanceCount::increment(this);
 }
 
@@ -91,7 +89,7 @@ Geant4ParticleHandler& Geant4ParticleHandler::operator=(const Geant4ParticleHand
 }
 
 /// Adopt the user particle handler
-bool Geant4ParticleHandler::adopt(Geant4UserParticleHandler* action)    {
+bool Geant4ParticleHandler::adopt(Geant4Action* action)    {
   if ( action )   {
     if ( !m_userHandler )  {
       Geant4UserParticleHandler* h = dynamic_cast<Geant4UserParticleHandler*>(action);
@@ -289,11 +287,19 @@ void Geant4ParticleHandler::end(const G4Track* track)   {
   int g4_id = h.id();
   int track_reason = m_currTrack.reason;
   PropertyMask mask(m_currTrack.reason);
-
-  double r_prod = sqrt(ph->vsx*ph->vsx + ph->vsy*ph->vsy);
-  double z_prod = fabs(ph->vsz);
-  if ( r_prod > m_rTracker || z_prod > m_zTracker )  {
-    m_currTrack.reason = 0;
+  // Update vertex end point and final momentum
+  G4ThreeVector m = track->GetMomentum();
+  const G4ThreeVector& p = track->GetPosition();
+  ph->pex = m.x();
+  ph->pey = m.y();
+  ph->pez = m.z();
+  ph->vex = p.x();
+  ph->vey = p.y();
+  ph->vez = p.z();
+
+  /// Final update of the particle using the user handler
+  if ( m_userHandler )  {
+    m_userHandler->end(track, m_currTrack);
   }
 
   // These are candate tracks with a probability to be stored due to their properties:
@@ -304,20 +310,7 @@ void Geant4ParticleHandler::end(const G4Track* track)   {
   // - to be kept due to creator process 
   //
   if ( !mask.isNull() )   {
-    // Update vertex end point and final momentum
-    G4ThreeVector m = track->GetMomentum();
-    const G4ThreeVector& p = track->GetPosition();
-    ph->pex = m.x();
-    ph->pey = m.y();
-    ph->pez = m.z();
-    ph->vex = p.x();
-    ph->vey = p.y();
-    ph->vez = p.z();
     m_equivalentTracks[g4_id] = g4_id;
-    /// Final update of the particle using the user handler
-    if ( m_userHandler )  {
-      m_userHandler->begin(track, m_currTrack);
-    }
     ParticleMap::iterator ip = m_particleMap.find(g4_id);
     if ( mask.isSet(G4PARTICLE_PRIMARY) )   {
       ph.dump2(outputLevel()-1,name(),"Add Primary",h.id(),ip!=m_particleMap.end());
@@ -346,7 +339,7 @@ void Geant4ParticleHandler::end(const G4Track* track)   {
     if ( ip != m_particleMap.end() )
       (*ip).second->reason |= track_reason;
     else
-      ph.dump3(outputLevel()+3,name(),"No real parent present");
+      ph.dump3(outputLevel()+3,name(),"FATAL: No real particle parent present");
   }
 }
 
@@ -482,6 +475,15 @@ int Geant4ParticleHandler::recombineParents()  {
     int g4_id = (*i).first;
     Particle* p = (*i).second;
     set<int>& daughters = p->daughters;
+
+    // Allow the user to force the particle handling either by
+    // or the reason mask with G4PARTICLE_KEEP_USER or
+    // to set the reason mask to NULL in order to drop it.
+    // Note: This may override all other decisions!
+    if ( m_userHandler )  {
+      m_userHandler->keepParticle(*p);
+    }
+
     if ( daughters.size() == 0 )   {
       PropertyMask mask(p->reason);
       int  id             =  p->id;
@@ -495,11 +497,19 @@ int Geant4ParticleHandler::recombineParents()  {
       bool remove_me      = false;
 
       if ( id == break_trackID )   {  // Used for debugging to set break point
-	remove_me      = false;
+	remove_me = false;
       }
 
-      /// Primary particles MUST be kept!
-      if ( mask.isSet(G4PARTICLE_PRIMARY) )   {
+      if ( mask.isSet(G4PARTICLE_KEEP_USER) )  {
+	/// If user decides it must be kept, it MUST be kept!
+	mask.set(G4PARTICLE_KEEP_USER);
+	continue;
+      }
+      else if ( mask.isSet(G4PARTICLE_PRIMARY) )   {
+	/// Primary particles MUST be kept!
+	continue;
+      }
+      else if ( mask.isSet(G4PARTICLE_KEEP_ALWAYS) )   {
 	continue;
       }
       else if ( keep_parent )  {
diff --git a/DDG4/src/Geant4UserParticleHandler.cpp b/DDG4/src/Geant4UserParticleHandler.cpp
index c4d1b67a8..0ed51d41c 100644
--- a/DDG4/src/Geant4UserParticleHandler.cpp
+++ b/DDG4/src/Geant4UserParticleHandler.cpp
@@ -54,3 +54,7 @@ void Geant4UserParticleHandler::end(const G4Track* /* track */, Particle& /* par
 /// Callback when parent should be combined
 void Geant4UserParticleHandler::combine(Particle& /* to_be_deleted */, Particle& /* remaining_parent */)  {
 }
+
+/// Callback to be answered if the particle MUST be kept during recombination step
+void Geant4UserParticleHandler::keepParticle(Particle& /* particle */)   {
+}
-- 
GitLab