Nonfiducial filter for ECAL pn

Biasing/include/Biasing/NonFiducialFilter.h

@@ -0,0 +1,74 @@
+/**
+ * @file NonFiducialFilter.h
+ * @class NonFiducialFilter
+ * @brief Class defining a UserActionPlugin that allows a user to filter out
+ *        non-fiducial events,  i.e.when  the electron is not contained
+ *        in the ECAL and so can act like the signal
+ * @author David Jiang, UCSB
+ * @author Tamas Almos Vami, UCSB
+ */
+
+#ifndef BIASING_NONFIDUCIALFILTER_H
+#define BIASING_NONFIDUCIALFILTER_H
+
+//----------------//
+//   C++ StdLib   //
+//----------------//
+#include <algorithm>
+
+/*~~~~~~~~~~~~~*/
+/*   SimCore   */
+/*~~~~~~~~~~~~~*/
+#include "SimCore/UserAction.h"
+
+/*~~~~~~~~~~~~~~~*/
+/*   Framework   */
+/*~~~~~~~~~~~~~~~*/
+#include "Framework/Configure/Parameters.h"
+#include "Framework/EventProcessor.h"
+
+namespace biasing {
+
+/**
+ * User action that allows a user to filter out events that are non-fiducial,
+ * i.e. the electron is not contained in the ECAL and so can act like the signal
+ */
+class NonFiducialFilter : public simcore::UserAction {
+ public:
+  /// Constructor
+  NonFiducialFilter(const std::string& name,
+                    framework::config::Parameters& parameters);
+
+  /// Destructor
+  virtual ~NonFiducialFilter() = default;
+
+  /**
+   * Implement the stepping action which performs the target volume biasing.
+   * @param step The Geant4 step.
+   */
+  void stepping(const G4Step* step) final override;
+
+  /**
+   * Method called at the end of every event.
+   * @param event Geant4 event object.
+   */
+  void EndOfEventAction(const G4Event*) final override;
+
+  /// Retrieve the type of actions this class defines
+  std::vector<simcore::TYPE> getTypes() final override {
+    return {simcore::TYPE::EVENT, simcore::TYPE::STACKING,
+            simcore::TYPE::STEPPING};
+  }
+
+ private:
+  /// Recoil electron threshold.
+  double recoil_max_p_{1500};  // MeV
+  /// If turned on, this aborts fiducial events.
+  bool abort_fiducial_{true};
+  /// Enable logging
+  enableLogging("NonFiducialFilter")
+
+};  // NonFiducialFilter
+}  // namespace biasing
+
+#endif  // BIASING_NONFIDUCIALFILTER_H

Biasing/python/ecal.py

@@ -75,6 +75,42 @@ def photo_nuclear( detector, generator ) :
 
     return sim
 
+def nonfiducial_photo_nuclear( detector, generator ) :
+
+    # Instantiate the simulator. 
+    sim = simulator.simulator("photo-nuclear")
+
+    # Set the path to the detector to use.
+    #   the second parameter says we want to include scoring planes
+    sim.setDetector( detector , True )
+
+    # Set run parameters
+    sim.description = "ECal photo-nuclear, xsec bias 450"
+    sim.beamSpotSmear = [20., 80., 0.] #mm
+
+    sim.generators.append( generator )
+
+    # Enable and configure the biasing
+    sim.biasing_operators = [ bias_operators.PhotoNuclear('ecal',450.,2500.,only_children_of_primary = True) ]
+
+    # the following filters are in a library that needs to be included
+    includeBiasing.library()
+
+    # Configure the sequence in which user actions should be called.
+    sim.actions.extend([
+            filters.TaggerVetoFilter(),
+            # Only consider events where a hard brem occurs
+            filters.TargetBremFilter(),
+            # Only considers events that are Non-Fiducial (Doesn't enter an ECal volume)
+            filters.NonFiducialFilter(),
+            # Only consider events where a PN reaction happens in the ECal
+            filters.EcalProcessFilter(),     
+            # Tag all photo-nuclear tracks to persist them to the event.
+            util.TrackProcessFilter.photo_nuclear()
+    ])
+
+    return sim
+
 def gamma_mumu(detector, generator) :
     """Example configuration for biasing gamma to mu+ mu- conversions in the ecal.
 

Biasing/python/filters.py

@@ -39,6 +39,26 @@ def __init__(self,recoil_max_p = 1500.,brem_min_e = 2500.) :
         self.brem_min_energy_threshold = brem_min_e
         self.kill_recoil_track = False
 
+class NonFiducialFilter(simcfg.UserAction):
+    """ Configuration for rejecting events that are fiducial.
+
+    Parameters
+    ----------
+    recoil_max_p : float
+        Maximum momentum the recoil electron can have [MeV]
+    abort_fiducial: bool
+        If true, aborts fiducial events. Otherwise, the fiducial events will be only tagged
+    """
+
+    def __init__(self,recoil_max_momentum = 1500., abort_fid_event = True) :
+        super().__init__("nonfiducial_filter", "biasing::NonFiducialFilter")
+
+        from LDMX.Biasing import include
+        include.library()
+
+        self.recoil_max_p = recoil_max_momentum
+        self.abort_fiducial = abort_fid_event
+
 class EcalProcessFilter(simcfg.UserAction):
     """ Configuration for filtering events that don't see a hard brem undergo a photo-nuclear reaction in the ECal. 
 

Biasing/src/Biasing/NonFiducialFilter.cxx

@@ -0,0 +1,110 @@
+#include "Biasing/NonFiducialFilter.h"
+
+/*~~~~~~~~~~~~*/
+/*   Geant4   */
+/*~~~~~~~~~~~~*/
+#include "G4EventManager.hh"
+#include "G4RunManager.hh"
+#include "G4Step.hh"
+#include "G4String.hh"
+#include "G4Track.hh"
+
+/*~~~~~~~~~~~~~*/
+/*   SimCore   */
+/*~~~~~~~~~~~~~*/
+#include "SimCore/UserEventInformation.h"
+#include "SimCore/UserTrackInformation.h"
+
+namespace biasing {
+
+bool nonFiducial_ = false;
+
+NonFiducialFilter::NonFiducialFilter(const std::string& name,
+                                     framework::config::Parameters& parameters)
+    : simcore::UserAction(name, parameters) {
+  recoil_max_p_ = parameters.getParameter<double>("recoil_max_p");
+  abort_fiducial_ = parameters.getParameter<bool>("abort_fiducial");
+}
+
+void NonFiducialFilter::stepping(const G4Step* step) {
+  // Get the track associated with this step.
+  auto track{step->GetTrack()};
+
+  // Get the PDG ID of the track and make sure it's an electron.
+  if (auto pdgID{track->GetParticleDefinition()->GetPDGEncoding()};
+      pdgID != 11) {
+    return;
+  }
+
+  // Only process the primary electron track
+  int parentID{step->GetTrack()->GetParentID()};
+  if (parentID != 0) {
+    return;
+  }
+
+  // Check in which volume the electron is currently
+  auto volume{track->GetVolume()->GetLogicalVolume()
+                  ? track->GetVolume()->GetLogicalVolume()->GetName()
+                  : "undefined"};
+
+  // Check if the track is tagged.
+  auto electronCheck{simcore::UserTrackInformation::get(track)};
+  if (electronCheck->isRecoilElectron() == true) {
+    if (track->GetMomentum().mag() > recoil_max_p_) {
+      // Kill the track if its momemntum is too high
+      track->SetTrackStatus(fKillTrackAndSecondaries);
+      G4RunManager::GetRunManager()->AbortEvent();
+      ldmx_log(debug) << " Recoil track momentum is too high, expected to be "
+                         "fiducial, exiting\n";
+      return;
+    }
+    // Check if the track ever enters the ECal. If it does, kill the track and
+    // abort the event.
+    auto isInEcal{((volume.contains("Si") || volume.contains("W") ||
+                    volume.contains("PCB") || volume.contains("strongback") ||
+                    volume.contains("Glue") || volume.contains("CFMix") ||
+                    volume.contains("Al") || volume.contains("C")) &&
+                   volume.contains("volume")) ||
+                  (volume.contains("nohole_motherboard"))};
+
+    // isInEcal should be taken from
+    // simcore::logical_volume_tests::isInEcal(volume) but for now it's under
+    // its own namespace so I cannot reach it here see issue
+    // https://github.com/LDMX-Software/ldmx-sw/issues/1286
+    if (abort_fiducial_ && isInEcal) {
+      track->SetTrackStatus(fKillTrackAndSecondaries);
+      G4RunManager::GetRunManager()->AbortEvent();
+      ldmx_log(debug) << ">> This event is fiducial, exiting";
+      nonFiducial_ = false;
+      return;
+    }
+    // I comment the following debug out since it would print per step and it's
+    // hard to read but it could be otherwise useful if somebody wants to do a
+    // step-by-step debugging ldmx_log(debug) << "  >> In this step this is
+    // non-fiducial, keeping it so far";
+    nonFiducial_ = true;
+    return;
+  } else {
+    // Check if the particle enters the recoil tracker.
+    if (volume.compareTo("recoil") == 0) {
+      /* Tag the tracks that:
+       1) Have a recoil electron
+       2) Enter/Exit the Target */
+      auto trackInfo{simcore::UserTrackInformation::get(track)};
+      trackInfo->tagRecoilElectron();  // tag the target recoil electron
+      ldmx_log(debug) << "  >> This track is the recoil electron, tagging it";
+      return;
+    }
+  }
+}
+
+void NonFiducialFilter::EndOfEventAction(const G4Event*) {
+  if (nonFiducial_) {
+    ldmx_log(debug) << "  >> This event is non-fiducial in ECAL, keeping it";
+  } else {
+    ldmx_log(debug) << ">> This event is fiducial, exiting";
+  }
+}
+}  // namespace biasing
+
+DECLARE_ACTION(biasing, NonFiducialFilter)

Biasing/test/ecal_pn_nonfiducial.py

@@ -0,0 +1,36 @@
+from LDMX.Framework import ldmxcfg
+p=ldmxcfg.Process("v14_nonfid")
+
+from LDMX.Ecal import EcalGeometry
+from LDMX.Hcal import HcalGeometry
+import LDMX.Ecal.ecal_hardcoded_conditions as ecal_conditions
+import LDMX.Hcal.hcal_hardcoded_conditions as hcal_conditions
+from LDMX.Biasing import ecal
+from LDMX.SimCore import generators
+
+mysim = ecal.nonfiducial_photo_nuclear('ldmx-det-v14-8gev', generators.single_8gev_e_upstream_tagger())
+mysim.description = "ECal Non-Fiducial Test Simulation"
+
+#from LDMX.Biasing import util
+#mysim.actions.append( util.StepPrinter(1) )
+
+import LDMX.Ecal.digi as ecal_digi
+import LDMX.Ecal.vetos as ecal_vetos
+
+p.outputFiles = [f'events_nonfiducial_test_production.root']
+p.histogramFile = f'hist_nonfiducial_test_production.root'
+
+p.maxTriesPerEvent = 10000
+p.maxEvents = 100
+p.run = 2
+p.logFrequency = 100
+#p.termLogLevel = 0
+
+p.sequence=[ mysim,
+        ecal_digi.EcalDigiProducer(),
+        ecal_digi.EcalRecProducer(),
+        ecal_vetos.EcalVetoProcessor()
+        ]
+
+from LDMX.DQM import dqm
+p.sequence.extend(dqm.ecal_dqm)