Add new filters

Biasing/include/Biasing/PhotoNuclearTopologyFilters.h

@@ -0,0 +1,124 @@
+#ifndef PHOTONUCLEARTOPOLOGYFILTERS_H
+#define PHOTONUCLEARTOPOLOGYFILTERS_H
+
+/*~~~~~~~~~~~~~*/
+/*   SimCore   */
+/*~~~~~~~~~~~~~*/
+#include "SimCore/UserAction.h"
+#include "SimCore/UserTrackInformation.h"
+/*~~~~~~~~~~~~~~~*/
+/*   Framework   */
+/*~~~~~~~~~~~~~~~*/
+#include <G4RunManager.hh>
+#include <G4Step.hh>
+
+#include "Framework/Configure/Parameters.h"
+// Forward declaration
+
+namespace biasing {
+
+/**
+ * Abstract base class for a user action used to filter out photo-nuclear events
+ * that don't match the topology the user is interested in studying. Derived
+ * classes implement the event rejection
+ *
+ * Similar to the PhotoNuclearProductsFilter, this user action will only process
+ * steps whose associated track has been tagged as a "PN Gamma". This tag is
+ * currently only set in ECalProcessFilter and needs to be placed in the
+ * UserAction pipeline before this class.
+ *
+ */
+class PhotoNuclearTopologyFilter : public simcore::UserAction {
+ public:
+  /**
+   * Constructor
+   *
+   * @param[in] name The name of this class instance.
+   * @param[in] parameters The parameters used to configure this class.
+   */
+  PhotoNuclearTopologyFilter(const std::string& name,
+                             framework::config::Parameters& parameters);
+
+  /// Destructor
+  ~PhotoNuclearTopologyFilter() = default;
+
+  /**
+   * Callback that allows a user to take some actions at the end of
+   * a step.
+   *
+   * @param[in] step The Geant4 step containing transient information
+   *      about the step taken by a track.
+   */
+  void stepping(const G4Step* step) override;
+
+  virtual bool rejectEvent(const std::vector<G4Track*>& secondaries) const = 0;
+
+  /// Retrieve the type of actions this class defines
+  std::vector<simcore::TYPE> getTypes() final override {
+    return {simcore::TYPE::STEPPING};
+  }
+
+ protected:
+  /**
+   *  Check if the PDG code corresponds to a light ion nucleus.
+   *
+   *  Nuclear PDG codes are given by ±10LZZZAAAI So to find the atomic
+   *  number, we first divide by 10 (to lose the I-component) and then
+   *  take the modulo with 1000.
+   *
+   *  TODO: Repeated code from SimCore, could probably live elsewhere.
+   *
+   */
+  constexpr bool isLightIon(const int pdgCode) const {
+    //
+    // TODO: Is the < check necessary?
+    if (pdgCode > 1000000000 && pdgCode < 10000000000) {
+      // Check if the atomic number is less than or equal to 4
+      return ((pdgCode / 10) % 1000) <= 4;
+    }
+    return false;
+  }
+
+  /**
+   * Whether or not to include a particular particle type in any counting.
+   * Unless \ref count_light_ions_ is set, we don't count anything with a
+   * nuclear PDG code. This is consistent with the counting behaviour used in
+   * the PhotoNuclearDQM.
+   *
+   * If \ref count_light_ions_ is set, we also match PDG codes for nuclei with
+   * atomic number < 4. 
+   *
+     TODO: Repeated code from SimCore, could probably live elsewhere.
+   * @see isLightIon
+   *
+   */
+  constexpr bool skipCountingParticle(const int pdgcode) const {
+    return !(pdgcode < 10000 || count_light_ions_ && isLightIon(pdgcode));
+  }
+
+  constexpr bool isNeutron(const int pdgID) const { return pdgID == 2112; }
+
+  bool count_light_ions_;
+  double hard_particle_threshold_;
+};  // PhotoNuclearTopologyFilter
+
+class SingleNeutronFilter : public PhotoNuclearTopologyFilter {
+ public:
+  SingleNeutronFilter(const std::string& name,
+                      framework::config::Parameters& parameters)
+      : PhotoNuclearTopologyFilter{name, parameters} {}
+
+  bool rejectEvent(const std::vector<G4Track*>& secondaries) const override;
+};
+
+class NothingHardFilter : public PhotoNuclearTopologyFilter {
+ public:
+  NothingHardFilter(const std::string& name,
+                    framework::config::Parameters& parameters)
+      : PhotoNuclearTopologyFilter{name, parameters} {}
+
+  bool rejectEvent(const std::vector<G4Track*>& secondaries) const override;
+};
+}  // namespace biasing
+
+#endif /* PHOTONUCLEARTOPOLOGYFILTERS_H */

Biasing/python/particle_filter.py

@@ -43,3 +43,69 @@ def kaon() :
                 321  # K^+
         ]
         return particle_filter
+
+class PhotoNuclearTopologyFilter(simcfg.UserAction):
+    """Configuration for keeping events with a PN reaction producing a
+    particular event topology.
+
+    Parameters
+    ----------
+    name : str
+        Name for this filter
+    filter_class:
+        Name of the class that implements this filter. Should correspond to the
+        invocation of DECLARE_ACTION in PhotoNuclearTopologyFilters.cxx
+
+    Attributes
+    ----------
+    hard_particle_threshold: float
+        The kinetic energy threshold required to count a particle as "hard"
+
+    count_light_ions: bool
+        Whether or not light ions (e.g. deutrons) should be considered when
+        applying the filter. Setting this to False can produce misleading
+        results such as classifying events with very high kinetic energy light
+        ions as "Nothing Hard" but it is the current default choice in the
+        PhotoNuclearDQM.
+
+    """
+    def __init__(self, name, filter_class):
+        super().__init__(name, filter_class)
+        from LDMX.Biasing import include
+        include.library()
+        self.count_light_ions = True
+
+
+    def SingleNeutronFilter():
+        """Configuration for keeping photonuclear events where a single neutron
+        carries most of the kinetic energy from the interaction
+
+
+        Returns
+        -------
+        Instance of configured PhotoNuclearTopologyFilter
+
+        """
+        filter = PhotoNuclearTopologyFilter(name='SingleNeutron',
+                                            filter_class='biasing::SingleNeutronFilter' )
+        filter.hard_particle_threshold = 200.
+        return filter
+
+    def NothingHardFilter():
+        """Configuration for keeping photonuclear events no particles received
+        significant kinetic energy.
+
+        Returns
+        -------
+        Instance of configured PhotoNuclearTopologyFilter
+
+        """
+        filter = PhotoNuclearTopologyFilter(name='NothingHard',
+                                            filter_class='biasing::NothingHardFilter')
+        filter.hard_particle_threshold = 200.
+        return filter
+
+
+
+
+

Biasing/src/Biasing/PhotoNuclearTopologyFilters.cxx

@@ -0,0 +1,79 @@
+#include "Biasing/PhotoNuclearTopologyFilters.h"
+
+namespace biasing {
+
+bool NothingHardFilter::rejectEvent(
+    const std::vector<G4Track*>& secondaries) const {
+  for (const auto& secondary : secondaries) {
+    // Get the PDG ID of the track
+    const auto pdgID{
+        std::abs(secondary->GetParticleDefinition()->GetPDGEncoding())};
+    if (skipCountingParticle(pdgID)) {
+      continue;
+    }
+    auto energy{secondary->GetKineticEnergy()};
+    if (energy > hard_particle_threshold_) {
+      return true;
+    }
+  }
+  return false;
+}
+bool SingleNeutronFilter::rejectEvent(
+    const std::vector<G4Track*>& secondaries) const {
+  int hard_particles{0};
+  int hard_neutrons{0};
+  for (const auto& secondary : secondaries) {
+    // Get the PDG ID of the track
+    const auto pdgID{
+        std::abs(secondary->GetParticleDefinition()->GetPDGEncoding())};
+    if (skipCountingParticle(pdgID)) {
+      continue;
+    }
+    auto energy{secondary->GetKineticEnergy()};
+    if (energy > hard_particle_threshold_) {
+      hard_particles++;
+      if (isNeutron(pdgID)) {
+        hard_neutrons++;
+      }
+    }
+  }
+  auto reject{hard_particles != hard_neutrons || hard_particles != 1};
+  return reject;
+}
+
+PhotoNuclearTopologyFilter::PhotoNuclearTopologyFilter(
+    const std::string& name, framework::config::Parameters& parameters)
+    : UserAction{name, parameters},
+      hard_particle_threshold_{
+          parameters.getParameter<double>("hard_particle_threshold")},
+      count_light_ions_{parameters.getParameter<bool>("count_light_ions")} {}
+
+void PhotoNuclearTopologyFilter::stepping(const G4Step* step) {
+  // Get the track associated with this step.
+  auto track{step->GetTrack()};
+
+  // Get the track info and check if this track has been tagged as the
+  // photon that underwent a photo-nuclear reaction. Only those tracks
+  // tagged as PN photos will be processed. The track is currently only
+  // tagged by the UserAction ECalProcessFilter which needs to be run
+  // before this UserAction.
+  auto trackInfo{simcore::UserTrackInformation::get(track)};
+  if ((trackInfo != nullptr) && !trackInfo->isPNGamma()) return;
+
+  // Get the PN photon daughters.
+  auto secondaries{step->GetSecondary()};
+
+  if (rejectEvent(*secondaries)) {
+    track->SetTrackStatus(fKillTrackAndSecondaries);
+    G4RunManager::GetRunManager()->AbortEvent();
+  }
+
+  // Once the PN gamma has been procesed, untag it so its not reprocessed
+  // again.
+  trackInfo->tagPNGamma(false);
+}
+
+}  // namespace biasing
+
+DECLARE_ACTION(biasing, NothingHardFilter)
+DECLARE_ACTION(biasing, SingleNeutronFilter)