diff --git a/Biasing/include/Biasing/NonFiducialFilter.h b/Biasing/include/Biasing/NonFiducialFilter.h new file mode 100644 index 000000000..db976c445 --- /dev/null +++ b/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 + +/*~~~~~~~~~~~~~*/ +/* 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 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 diff --git a/Biasing/python/ecal.py b/Biasing/python/ecal.py index 3b9bb2c17..b222d476f 100644 --- a/Biasing/python/ecal.py +++ b/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. diff --git a/Biasing/python/filters.py b/Biasing/python/filters.py index f2b055f6a..1b5e670a2 100644 --- a/Biasing/python/filters.py +++ b/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. diff --git a/Biasing/src/Biasing/NonFiducialFilter.cxx b/Biasing/src/Biasing/NonFiducialFilter.cxx new file mode 100644 index 000000000..f92a56a8a --- /dev/null +++ b/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("recoil_max_p"); + abort_fiducial_ = parameters.getParameter("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) diff --git a/Biasing/test/ecal_pn_nonfiducial.py b/Biasing/test/ecal_pn_nonfiducial.py new file mode 100644 index 000000000..32e1d1b99 --- /dev/null +++ b/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)