1237 add an option to print steps of the daughters of the particle of interest

Biasing/include/Biasing/PhotoNuclearProductsFilter.h

@@ -58,6 +58,7 @@ class PhotoNuclearProductsFilter : public simcore::UserAction {
  private:
   /// Container to hold the PDG IDs of products of interest
   std::vector<int> productsPdgID_;
+    double min_e;
 
 };  // PhotoNuclearProductsFilter
 

Biasing/include/Biasing/Utility/StepPrinter.h

@@ -21,27 +21,36 @@ class StepPrinter : public simcore::UserAction {
    * @param[in] parameters the parameters used to configure this
    *      UserAction.
    */
-  StepPrinter(const std::string& name, framework::config::Parameters& parameters);
+  StepPrinter(const std::string& name,
+              framework::config::Parameters& parameters);
 
   /// Destructor
-  ~StepPrinter();
+  virtual ~StepPrinter() = default;
 
   /**
    * Stepping action called when a step is taken during tracking of
    * a particle.
    *
    * @param[in] step Geant4 step
    */
-  void stepping(const G4Step* step) final override;
+  void stepping(const G4Step* step) override;
 
   /// Retrieve the type of actions this class defines
-  std::vector<simcore::TYPE> getTypes() final override {
+  std::vector<simcore::TYPE> getTypes() override {
     return {simcore::TYPE::STEPPING};
   }
 
  private:
   /// The track ID to filter on
   int trackID_{-9999};
+  std::string processName_{"UNDEFINED"};
+  int depth_{0};
+  std::unordered_map<int,int> trackParents_{};
+
+
+  /// Check if the given track is a descendent of the track we are interested in
+  /// up to a certain depth
+  bool isDescendent(const G4Track* track) const;
 
 };  // StepPrinter
 

Biasing/python/particle_filter.py

@@ -25,6 +25,7 @@ def __init__(self,name) :
         include.library()
 
         self.pdg_ids = [ ]
+        self.min_e = 0.
 
     def kaon() :
         """ Configuration for filtering photo-nuclear events whose products don't contain a kaon.

Biasing/python/util.py

@@ -30,10 +30,11 @@ class StepPrinter(BiasingUtilityAction) :
         Geant4 track ID to print each step of
     """
 
-    def __init__(self,track_id=1) :
+    def __init__(self,track_id=1, process_name='', depth=0) :
         super().__init__('print_steps_%s'%track_id,'StepPrinter')
-
+        self.process_name = process_name
         self.track_id = track_id
+        self.depth = depth
 
 class PartialEnergySorter(BiasingUtilityAction) :
     """Process particles such that all particles above

Biasing/src/Biasing/PhotoNuclearProductsFilter.cxx

@@ -23,6 +23,7 @@ PhotoNuclearProductsFilter::PhotoNuclearProductsFilter(
     const std::string& name, framework::config::Parameters& parameters)
     : simcore::UserAction(name, parameters) {
   productsPdgID_ = parameters.getParameter<std::vector<int> >("pdg_ids");
+  min_e = parameters.getParameter<double>("min_e");
 }
 
 PhotoNuclearProductsFilter::~PhotoNuclearProductsFilter() {}
@@ -54,7 +55,9 @@ void PhotoNuclearProductsFilter::stepping(const G4Step* step) {
     // Check if the PDG ID is in the list of products of interest
     if (std::find(productsPdgID_.begin(), productsPdgID_.end(), pdgID) !=
         productsPdgID_.end()) {
-      productFound = true;
+      if (secondary->GetKineticEnergy() > min_e) {
+        productFound = true;
+      }
       break;
     }
   }

Biasing/src/Biasing/Utility/StepPrinter.cxx

@@ -9,43 +9,94 @@
 namespace biasing {
 namespace utility {
 
-StepPrinter::StepPrinter(const std::string& name, framework::config::Parameters& parameters)
+StepPrinter::StepPrinter(const std::string& name,
+                         framework::config::Parameters& parameters)
     : simcore::UserAction(name, parameters) {
   trackID_ = parameters.getParameter<int>("track_id");
+  processName_ = parameters.getParameter<std::string>("process_name");
+  depth_ = parameters.getParameter<int>("depth");
 }
 
-StepPrinter::~StepPrinter() {}
-
+bool StepPrinter::isDescendent(const G4Track* track) const {
+  auto trackID{track->GetTrackID()};
+  int currentDepth{0};
+  int currentTrackID{trackID};
+  // Walk the tree until we either no longer have a parent or we reach the
+  // desired depth
+  while (currentDepth < depth_ &&
+         trackParents_.find(currentTrackID) != trackParents_.end()) {
+    // See if we have encountered the parent of the current track
+    //
+    // operator[] is not const, so we need to use at()
+    currentTrackID = trackParents_.at(currentTrackID);
+    if (currentTrackID == trackID_) {
+      // If one of the parents is the track of interest, we are done!
+      return true;
+    }
+    currentTrackID = track->GetParentID();
+    currentDepth++;
+  }
+  return false;
+}
 void StepPrinter::stepping(const G4Step* step) {
   // Get the track associated with this step
   auto track{step->GetTrack()};
 
-  if (auto trackID{track->GetTrackID()};
-      (trackID_ > 0) && (trackID != trackID_))
+  const auto trackID{track->GetTrackID()};
+  const auto parent{track->GetParentID()};
+  // Don't bother filling the map if we aren't going to use it
+  if (depth_ > 0) {
+    trackParents_[trackID] = parent;
+  }
+
+  auto process{track->GetCreatorProcess()};
+  std::string processName{process ? process->GetProcessName() : "Primary"};
+  // Unwrap biasing part of process name if present
+  if (processName.find("biasWrapper") != std::string::npos) {
+    std::size_t pos = processName.find_first_of("(") + 1;
+    processName = processName.substr(pos, processName.size() - pos - 1);
+  }
+
+  // This could be a negated condition, but it is easier to read this way
+  //
+  if (trackID == trackID_ || // We are the track of interest
+      isDescendent(track) || // We are a descendent of the track of interest
+      processName == processName_ // The parent process was the process of interest
+  ) {
+    // This is an interesting track -> Carry on processing
+  }  else {
     return;
-
+  }
   // Get the particle name.
-  auto particleName{track->GetParticleDefinition()->GetParticleName()};
+  const auto particleName{track->GetParticleDefinition()->GetParticleName()};
 
   // Get the energy of the particle
-  auto energy{step->GetPostStepPoint()->GetTotalEnergy()};
+  const auto energy{step->GetPostStepPoint()->GetTotalEnergy()};
 
   // Get the volume the particle is in.
-  auto volume{track->GetVolume()->GetName()};
+  auto volume{track->GetVolume()};
+  auto volumeName{volume->GetName()};
 
-  // Get the next volume
-  auto nextVolume{track->GetNextVolume()->GetName()};
+  // Get the next volume (can fail if current volume is WorldPV and next is
+  // outside the world)
+  auto nextVolume{track->GetNextVolume() ? track->GetNextVolume()->GetName()
+                                         : "undefined"};
 
   // Get the region
-  auto region{track->GetVolume()->GetLogicalVolume()->GetRegion()->GetName()};
+  auto regionName{volume->GetLogicalVolume()->GetRegion()->GetName()};
 
-  std::cout << " Step " << track->GetCurrentStepNumber() << " {"
+  std::cout << " Step " << track->GetCurrentStepNumber() << " ("
+            << track->GetParticleDefinition()->GetParticleName() << ") {"
             << " Energy: " << energy << " Track ID: " << track->GetTrackID()
-            << " Particle currently in: " << volume << " Region: " << region
-            << " Next volume: " << nextVolume
-            << " Weight: " << track->GetWeight() << " Children:";
-  for (auto const& track : *(step->GetSecondaryInCurrentStep()))
-    std::cout << " " << track->GetParticleDefinition()->GetPDGEncoding();
+            << " Particle currently in: " << volumeName
+            << " Region: " << regionName << " Next volume: " << nextVolume
+            << " Weight: " << track->GetWeight() << " Parent: " << parent
+            << " (" << processName << ") "
+            << " Children:";
+  for (auto const& child : *(step->GetSecondaryInCurrentStep())) {
+    std::cout << " (" << child->GetTotalEnergy()
+              << "): " << child->GetParticleDefinition()->GetPDGEncoding();
+  }
 
   std::cout << " }" << std::endl;
 }