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/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 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 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)