From 8d6f86a783c67d54a6a819bec4a683b7f6bd8d40 Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Mon, 16 Oct 2023 13:21:45 -0500 Subject: [PATCH 1/3] 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/MidShowerDiMuonBkgdFilter.h | 142 ++++++++++++++++++ .../src/Biasing/MidShowerDiMuonBkgdFilter.cxx | 117 +++++++++++++++ 2 files changed, 259 insertions(+) create mode 100644 Biasing/include/Biasing/MidShowerDiMuonBkgdFilter.h create mode 100644 Biasing/src/Biasing/MidShowerDiMuonBkgdFilter.cxx diff --git a/Biasing/include/Biasing/MidShowerDiMuonBkgdFilter.h b/Biasing/include/Biasing/MidShowerDiMuonBkgdFilter.h new file mode 100644 index 000000000..8f6138dba --- /dev/null +++ b/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 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__ diff --git a/Biasing/src/Biasing/MidShowerDiMuonBkgdFilter.cxx b/Biasing/src/Biasing/MidShowerDiMuonBkgdFilter.cxx new file mode 100644 index 000000000..c4e96c3bb --- /dev/null +++ b/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("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) From a267c721fd50b82dfed7089916451ddaa09859b2 Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Mon, 16 Oct 2023 14:51:42 -0500 Subject: [PATCH 2/3] add config class for midshower dimuon filter --- Biasing/python/filters.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/Biasing/python/filters.py b/Biasing/python/filters.py index bebd6b087..f2b055f6a 100644 --- a/Biasing/python/filters.py +++ b/Biasing/python/filters.py @@ -182,3 +182,20 @@ def __init__(self,thresh) : include.library() self.threshold = thresh + +class MidShowerDiMuonBkgdFilter(simcfg.UserAction) : + """ Configuration used to reject events that don't have enough energy given to muon conversion + + Parameters + ---------- + thresh : float + Minimum energy [MeV] that needs to go into creating the muons + """ + + def __init__(self,thresh) : + super().__init__('midshower_dimuon_min_%d_MeV'%(thresh),'biasing::MidShowerDiMuonBkgdFilter') + + from LDMX.Biasing import include + include.library() + + self.threshold = thresh From 0fa56edf950fcb405d02fa6ecd47cba606fce53c Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Mon, 16 Oct 2023 14:52:40 -0500 Subject: [PATCH 3/3] add example of creating simulation for midshower dimuon --- Biasing/python/eat.py | 128 +++++++++++++++++++++++++++++++++++------- 1 file changed, 108 insertions(+), 20 deletions(-) diff --git a/Biasing/python/eat.py b/Biasing/python/eat.py index d16695e74..f7f3140e2 100644 --- a/Biasing/python/eat.py +++ b/Biasing/python/eat.py @@ -15,7 +15,7 @@ from LDMX.Biasing import filters from LDMX.Biasing import util -def midshower_nuclear( detector , bias_factor , bias_threshold , min_nuclear_energy ) : +def midshower_nuclear( detector , generator, bias_factor , bias_threshold , min_nuclear_energy ) : """Example configuration for producing mid-shower nuclear interactions in the ECal that fake a missing energy (ME) signal. @@ -28,6 +28,10 @@ def midshower_nuclear( detector , bias_factor , bias_threshold , min_nuclear_ene detector : str Name of the detector + generator : simcfg.PrimaryGenerator + Beam generator for this simulation which should be a ParticleGun + so we can configure the PrimaryToEcalFilter to select events where + beam electrons retain 87.5% of their energy. bias_factor : float Factor to multiply EN/PN xsecs by bias_threshold : float @@ -44,7 +48,14 @@ def midshower_nuclear( detector , bias_factor , bias_threshold , min_nuclear_ene ------- from LDMX.Biasing import eat - eat_pn_sim = eat.midshower_nuclear('ldmx-det-v12',1000.,1700.) + from LDMX.SimCore import generators + eat_dimuon = eat.midshower_nuclear( + 'ldmx-det-v14', + generators.single_e_4gev_upstream_tagger(), + 200, + 1500., + 2500. + ) """ @@ -60,8 +71,7 @@ def midshower_nuclear( detector , bias_factor , bias_threshold , min_nuclear_ene sim.description = "Biased Mid-Shower Nuclear Interactions ME Background" sim.beamSpotSmear = [20., 80., 0.] #mm - from LDMX.SimCore import generators - sim.generators = [ generators.single_4gev_e_upstream_tagger() ] + sim.generators = [ generator ] #Enable and configure the biasing from LDMX.SimCore import bias_operators @@ -72,8 +82,7 @@ def midshower_nuclear( detector , bias_factor , bias_threshold , min_nuclear_ene #Configure the sequence in which user actions should be called. sim.actions = [ - #get electron to ECal with 3500MeV - filters.PrimaryToEcalFilter(3500.), + filters.PrimaryToEcalFilter(0.875*generator.energy*1000), #Make sure all particles above 1500MeV are processed first util.PartialEnergySorter(bias_threshold), #Make sure a total of 1700MeV energy went PN in ECal @@ -82,7 +91,81 @@ def midshower_nuclear( detector , bias_factor , bias_threshold , min_nuclear_ene return sim -def dark_brem( ap_mass , lhe, detector ) : + +def midshower_dimuon( detector , generator, bias_factor , bias_threshold , min_dimuon_energy ) : + """Example configuration for producing mid-shower dimuon interactions in + the ECal that fake a missing energy (ME) signal. + + In this particular example, 4 GeV electrons are fired upstream of the + tagger tracker. Then simulation events are then put through a series of + filters. + + Parameters + ---------- + + detector : str + Name of the detector + generator : simcfg.PrimaryGenerator + Beam generator for this simulation which should be a ParticleGun + so we can configure the PrimaryToEcalFilter to select events where + beam electrons retain 87.5% of their energy. + bias_factor : float + Factor to multiply GammaToMuPair xsec by within the ECal + bias_threshold : float + Minium energy [MeV] to bias a photon + min_dimuon_energy : float + Minium total energy [MeV] that went to muons to keep event + + Returns + ------- + simulator : + configured for mid-shower dimuon interactions and far-upstream generator + + Example + ------- + + from LDMX.Biasing import eat + from LDMX.SimCore import generators + eat_dimuon = eat.midshower_dimuon( + 'ldmx-det-v14', + generators.single_e_4gev_upstream_tagger(), + 1e4, + 500., + 1000. + ) + + """ + + #Instantiate the simulator. + sim = simulator.simulator("eat-midshower-dimuon") + from LDMX.Ecal import EcalGeometry + + #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 = "Biased Mid-Shower DiMuon Interactions ME Background" + sim.beamSpotSmear = [20., 80., 0.] #mm + + from LDMX.SimCore import generators + sim.generators = [ generator ] + + #Enable and configure the biasing + from LDMX.SimCore import bias_operators + sim.biasing_operators = [ bias_operators.GammaToMuPair('ecal', bias_factor, bias_threshold) ] + + #Configure the sequence in which user actions should be called. + sim.actions = [ + filters.PrimaryToEcalFilter(0.875*generator.energy*1000), + util.PartialEnergySorter(bias_threshold), + filters.MidShowerDiMuonBkgdFilter(min_dimuon_energy), + ] + + return sim + + +def dark_brem( ap_mass , lhe, detector, generator) : """Example configuration for producing dark brem interactions in the ECal. This configures the simulator to fire a 4 GeV electron upstream of the @@ -95,9 +178,13 @@ def dark_brem( ap_mass , lhe, detector ) : ap_mass : float The mass of the A' in MeV. lhe : str - The path to the LHE file to use as vertices of the dark brem. + The path to the reference library to use as vertices of the dark brem. detector : str Path to the detector. + generator : simcfg.PrimaryGenerator + Beam generator for this simulation which should be a ParticleGun + so we can configure the PrimaryToEcalFilter to select events where + beam electrons retain 87.5% of their energy. Return ------ @@ -105,11 +192,14 @@ def dark_brem( ap_mass , lhe, detector ) : Example ------- - Here we use the example vertex library. This should not be used - for large (>50k) event samples. - from LDMX.SimCore import makePath - eat_ap_sim = eat.dark_brem(1000., makePath.makeLHEPath(1000.), 'ldmx-det-v12') + from LDMX.SimCore import generators + eat_ap_sim = eat.dark_brem( + 1000., + 'path/to/1GeV_mA_library', + 'ldmx-det-v14', + generators.single_e_4gev_upstream_tagger() + ) """ @@ -118,13 +208,13 @@ def dark_brem( ap_mass , lhe, detector ) : sim.description = "One e- fired far upstream with Dark Brem turned on and biased up in ECal" sim.setDetector( detector , True ) - sim.generators.append( generators.single_4gev_e_upstream_tagger() ) + sim.generators = [ generator ] sim.beamSpotSmear = [ 20., 80., 0. ] #mm #Activiate dark bremming with a certain A' mass and LHE library from LDMX.SimCore import dark_brem db_model = dark_brem.G4DarkBreMModel( lhe ) - db_model.threshold = 2. #GeV - minimum energy electron needs to have to dark brem + db_model.threshold = 0.5*generator.energy #GeV - minimum energy electron needs to have to dark brem db_model.epsilon = 0.01 #decrease epsilon from one to help with Geant4 biasing calculations sim.dark_brem.activate( ap_mass , db_model ) @@ -137,13 +227,11 @@ def dark_brem( ap_mass , lhe, detector ) : ) ] + beam_energy = generator.energy*1000 sim.actions = [ - #Make sure all particles above 2GeV are processed first - util.PartialEnergySorter(2000.), - #Abort events if the electron doesn't get to the ECal with 3.5GeV - filters.PrimaryToEcalFilter(3500.), - #Only keep events when a dark brem happens in the Ecal - filters.EcalDarkBremFilter(2000.) + util.PartialEnergySorter(0.5*beam_energy), + filters.PrimaryToEcalFilter(0.875*beam_energy), + filters.EcalDarkBremFilter(0.5*beam_energy) ] return sim