initial copy of MidShowerNuclear but re-worked for Muon-Conv

initial testing looks like it is functional (i.e. debug printouts when
uncommented show that events are only being kept when muon pairs are
generated at above the input threshold energy), will need to check the
physical distributions after generating an actual sample.

Biasing/include/Biasing/MidShowerDiMuonBkgdFilter.h

@@ -0,0 +1,142 @@
+#ifndef BIASING_MIDSHOWERDIMUONBKGDFILTER_H_
+#define BIASING_MIDSHOWERDIMUONBKGDFILTER_H_
+
+/*~~~~~~~~~~~~*/
+/*   SimCore  */
+/*~~~~~~~~~~~~*/
+#include "SimCore/UserAction.h"
+
+/// Forward declaration of virtual process class
+class G4VProcess;
+
+namespace biasing {
+
+/**
+ * @class MidShowerDiMuonBkgdFilter
+ *
+ * The basic premis of this filter is to add up all of the
+ * energy "lost" to muons created within the calorimeters.
+ * When the PartialEnergySorter has run out of "high" energy 
+ * particles to process (when NewStage is called)
+ * we check if the running total is high enough to keep the event.
+ *
+ * @see PartialEnergySorter
+ * Here we assume that the partial energy sorter is being run in sequence with
+ * this filter.
+ */
+class MidShowerDiMuonBkgdFilter : public simcore::UserAction {
+ public:
+  /**
+   * Class constructor.
+   *
+   * Retrieve the necessary configuration parameters
+   */
+  MidShowerDiMuonBkgdFilter(const std::string& name,
+                             framework::config::Parameters& parameters);
+
+  /**
+   * Class destructor.
+   */
+  ~MidShowerDiMuonBkgdFilter() {}
+
+  /**
+   * Get the types of actions this class can do
+   *
+   * @return list of action types this class does
+   */
+  std::vector<simcore::TYPE> getTypes() final override {
+    return {simcore::TYPE::STACKING, simcore::TYPE::STEPPING,
+            simcore::TYPE::EVENT};
+  }
+
+  /**
+   * Reset the total energy going to the muons.
+   *
+   * @param[in] event not used
+   */
+  void BeginOfEventAction(const G4Event* event) final override;
+
+  /**
+   * We follow the simulation along each step and check
+   * if any secondaries of the input process were created.
+   *
+   * If they were created, we add the change in energy
+   * to the total energy "lost" to the input process.
+   *
+   * @note Only includes the steps that are within the
+   * 'CalorimeterRegion' in the search for interesting
+   * products.
+   *
+   * @see isOutsideCalorimeterRegion
+   * @see anyCreatedViaProcess
+   *
+   * @param[in] step current G4Step
+   */
+  void stepping(const G4Step* step) final override;
+
+  /**
+   * When using the PartialEnergySorter,
+   * the *first* time that a new stage begins is when
+   * all particles are now below the threshold.
+   *
+   * We take this opportunity to make sure that
+   * enough energy has gone to the products of
+   * the input process.
+   *
+   * @see PartialEnergySort::NewStage
+   * @see AbortEvent
+   */
+  void NewStage() final override;
+
+ private:
+  /**
+   * Checks if the passed step is outside of the CalorimeterRegion.
+   *
+   * @param[in] step const G4Step to check
+   * @returns true if passed step is outside of the CalorimeterRegion.
+   */
+  bool isOutsideCalorimeterRegion(const G4Step* step) const;
+
+  /**
+   * Helper to save the passed track
+   *
+   * @note Assumes user track information has already been created
+   * for the input track.
+   *
+   * @param[in] track G4Track to persist into output
+   */
+  void save(const G4Track* track) const;
+
+  /**
+   * Helper to abort an event with a message
+   *
+   * Tells the RunManger to abort the current event
+   * after displaying the input message.
+   *
+   * @param[in] reason reason for aborting the event
+   */
+  void AbortEvent(const std::string& reason) const;
+
+ private:
+  /**
+   * Minimum energy [MeV] that the process products need to have
+   * to keep the event.
+   *
+   * Also used by PartialEnergySorter to determine
+   * which tracks should be processed first.
+   *
+   * Parameter Name: 'threshold'
+   */
+  double threshold_;
+
+  /**
+   * Total energy gone to the process in the current event
+   *
+   * Reset to 0. in BeginOfEventAction
+   */
+  double total_process_energy_{0.};
+
+};  // MidShowerDiMuonBkgdFilter
+}  // namespace biasing
+
+#endif  // BIASING_MIDSHOWERDIMUONBKGDFILTER_H__

Biasing/src/Biasing/MidShowerDiMuonBkgdFilter.cxx

@@ -0,0 +1,117 @@
+#include "Biasing/MidShowerDiMuonBkgdFilter.h"
+
+#include "G4EventManager.hh"
+#include "G4RunManager.hh"
+#include "G4Step.hh"
+#include "G4Gamma.hh"
+#include "SimCore/UserTrackInformation.h"
+
+namespace biasing {
+
+MidShowerDiMuonBkgdFilter::MidShowerDiMuonBkgdFilter(const std::string& name,
+                                                       framework::config::Parameters& parameters)
+    : simcore::UserAction(name, parameters) {
+  threshold_ = parameters.getParameter<double>("threshold");
+  /* debug printout
+  std::cout
+      << "[ MidShowerDiMuonBkgdFilter ]: "
+      << "created instance " << name
+      << "with threshold " << threshold_ << " MeV"
+      << std::endl;
+   */
+}
+
+void MidShowerDiMuonBkgdFilter::BeginOfEventAction(const G4Event*) {
+  /* debug printout
+  std::cout
+      << "[ MidShowerDiMuonBkgdFilter ]: "
+      << "("
+      << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID()
+      << ") "
+      << "starting new simulation event." << std::endl;
+   */
+  total_process_energy_ = 0.;
+}
+
+void MidShowerDiMuonBkgdFilter::stepping(const G4Step* step) {
+  // skips steps that aren't done by photons (required for mu conv)
+  if (step->GetTrack()->GetParticleDefinition() != G4Gamma::Gamma()) return;
+  // skip steps that are outside the calorimeter region
+  if (isOutsideCalorimeterRegion(step)) return;
+  // check photon secondaries for muons
+  const G4TrackVector* secondaries{step->GetSecondary()};
+  if (secondaries == nullptr or secondaries->size() == 0) return;
+  // have left if secondaries is empty
+  bool found_muon{false};
+  for (const G4Track* secondary : *secondaries) {
+    if (abs(secondary->GetParticleDefinition()->GetPDGEncoding()) == 13) {
+      found_muon = true;
+      save(secondary);
+    }
+  }
+  if (not found_muon) return;
+  // there are interesting secondaries in this step
+  double pre_energy = step->GetPreStepPoint()->GetTotalEnergy();
+  double post_energy = step->GetPostStepPoint()->GetTotalEnergy();
+  total_process_energy_ += pre_energy - post_energy;
+  const G4Track* track = step->GetTrack();
+  /* debug printout
+  std::cout << "[ MidShowerDiMuonBkgdFilter ]: "
+            << "("
+            << G4EventManager::GetEventManager()
+                   ->GetConstCurrentEvent()
+                   ->GetEventID()
+            << ") "
+            << "Track " << track->GetParentID() << " created "
+            << track->GetTrackID() << " which went from " << pre_energy
+            << " MeV to " << post_energy << " MeV generating muons." 
+            << std::endl;
+   */
+  // make sure this track is saved
+  save(track);
+}
+
+void MidShowerDiMuonBkgdFilter::NewStage() {
+  /* debug printout
+  std::cout
+      << "[ MidShowerDiMuonBkgdFilter ]: "
+      << "("
+      << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID()
+      << ") " << total_process_energy_ << " MeV was muonic." << std::endl;
+   */
+  if (total_process_energy_ < threshold_)
+    AbortEvent("Not enough energy went to the input process.");
+  return;
+}
+
+bool MidShowerDiMuonBkgdFilter::isOutsideCalorimeterRegion(
+    const G4Step* step) const {
+  // the pointers in this chain are assumed to be always valid
+  auto reg{step->GetTrack()->GetVolume()->GetLogicalVolume()->GetRegion()};
+  if (reg) return (reg->GetName() != "CalorimeterRegion");
+  // region is nullptr ==> no region defined for current volume
+  //  ==> outside CalorimeterRegion
+  return true;
+}
+
+void MidShowerDiMuonBkgdFilter::save(const G4Track* track) const {
+  auto track_info{simcore::UserTrackInformation::get(track)};
+  track_info->setSaveFlag(true);
+  return;
+}
+
+void MidShowerDiMuonBkgdFilter::AbortEvent(const std::string& reason) const {
+  if (G4RunManager::GetRunManager()->GetVerboseLevel() > 1) {
+    std::cout << "[ MidShowerDiMuonBkgdFilter ]: "
+              << "("
+              << G4EventManager::GetEventManager()
+                     ->GetConstCurrentEvent()
+                     ->GetEventID()
+              << ") " << reason << " Aborting event." << std::endl;
+  }
+  G4RunManager::GetRunManager()->AbortEvent();
+  return;
+}
+}  // namespace biasing
+
+DECLARE_ACTION(biasing, MidShowerDiMuonBkgdFilter)