diff --git a/Biasing/include/Biasing/PhotoNuclearTopologyFilters.h b/Biasing/include/Biasing/PhotoNuclearTopologyFilters.h new file mode 100644 index 000000000..1c5abf248 --- /dev/null +++ b/Biasing/include/Biasing/PhotoNuclearTopologyFilters.h @@ -0,0 +1,124 @@ +#ifndef PHOTONUCLEARTOPOLOGYFILTERS_H +#define PHOTONUCLEARTOPOLOGYFILTERS_H + +/*~~~~~~~~~~~~~*/ +/* SimCore */ +/*~~~~~~~~~~~~~*/ +#include "SimCore/UserAction.h" +#include "SimCore/UserTrackInformation.h" +/*~~~~~~~~~~~~~~~*/ +/* Framework */ +/*~~~~~~~~~~~~~~~*/ +#include +#include + +#include "Framework/Configure/Parameters.h" +// Forward declaration + +namespace biasing { + +/** + * Abstract base class for a user action used to filter out photo-nuclear events + * that don't match the topology the user is interested in studying. Derived + * classes implement the event rejection + * + * Similar to the PhotoNuclearProductsFilter, this user action will only process + * steps whose associated track has been tagged as a "PN Gamma". This tag is + * currently only set in ECalProcessFilter and needs to be placed in the + * UserAction pipeline before this class. + * + */ +class PhotoNuclearTopologyFilter : public simcore::UserAction { + public: + /** + * Constructor + * + * @param[in] name The name of this class instance. + * @param[in] parameters The parameters used to configure this class. + */ + PhotoNuclearTopologyFilter(const std::string& name, + framework::config::Parameters& parameters); + + /// Destructor + ~PhotoNuclearTopologyFilter() = default; + + /** + * Callback that allows a user to take some actions at the end of + * a step. + * + * @param[in] step The Geant4 step containing transient information + * about the step taken by a track. + */ + void stepping(const G4Step* step) override; + + virtual bool rejectEvent(const std::vector& secondaries) const = 0; + + /// Retrieve the type of actions this class defines + std::vector getTypes() final override { + return {simcore::TYPE::STEPPING}; + } + + protected: + /** + * Check if the PDG code corresponds to a light ion nucleus. + * + * Nuclear PDG codes are given by ±10LZZZAAAI So to find the atomic + * number, we first divide by 10 (to lose the I-component) and then + * take the modulo with 1000. + * + * TODO: Repeated code from SimCore, could probably live elsewhere. + * + */ + constexpr bool isLightIon(const int pdgCode) const { + // + // TODO: Is the < check necessary? + if (pdgCode > 1000000000 && pdgCode < 10000000000) { + // Check if the atomic number is less than or equal to 4 + return ((pdgCode / 10) % 1000) <= 4; + } + return false; + } + + /** + * Whether or not to include a particular particle type in any counting. + * Unless \ref count_light_ions_ is set, we don't count anything with a + * nuclear PDG code. This is consistent with the counting behaviour used in + * the PhotoNuclearDQM. + * + * If \ref count_light_ions_ is set, we also match PDG codes for nuclei with + * atomic number < 4.  + * + TODO: Repeated code from SimCore, could probably live elsewhere. + * @see isLightIon + * + */ + constexpr bool skipCountingParticle(const int pdgcode) const { + return !(pdgcode < 10000 || count_light_ions_ && isLightIon(pdgcode)); + } + + constexpr bool isNeutron(const int pdgID) const { return pdgID == 2112; } + + bool count_light_ions_; + double hard_particle_threshold_; +}; // PhotoNuclearTopologyFilter + +class SingleNeutronFilter : public PhotoNuclearTopologyFilter { + public: + SingleNeutronFilter(const std::string& name, + framework::config::Parameters& parameters) + : PhotoNuclearTopologyFilter{name, parameters} {} + + bool rejectEvent(const std::vector& secondaries) const override; +}; + +class NothingHardFilter : public PhotoNuclearTopologyFilter { + public: + NothingHardFilter(const std::string& name, + framework::config::Parameters& parameters) + : PhotoNuclearTopologyFilter{name, parameters} {} + + bool rejectEvent(const std::vector& secondaries) const override; +}; +} // namespace biasing + +#endif /* PHOTONUCLEARTOPOLOGYFILTERS_H */ diff --git a/Biasing/python/particle_filter.py b/Biasing/python/particle_filter.py index 6fea4593b..300dad13e 100644 --- a/Biasing/python/particle_filter.py +++ b/Biasing/python/particle_filter.py @@ -43,3 +43,69 @@ def kaon() : 321 # K^+ ] return particle_filter + +class PhotoNuclearTopologyFilter(simcfg.UserAction): + """Configuration for keeping events with a PN reaction producing a + particular event topology. + + Parameters + ---------- + name : str + Name for this filter + filter_class: + Name of the class that implements this filter. Should correspond to the + invocation of DECLARE_ACTION in PhotoNuclearTopologyFilters.cxx + + Attributes + ---------- + hard_particle_threshold: float + The kinetic energy threshold required to count a particle as "hard" + + count_light_ions: bool + Whether or not light ions (e.g. deutrons) should be considered when + applying the filter. Setting this to False can produce misleading + results such as classifying events with very high kinetic energy light + ions as "Nothing Hard" but it is the current default choice in the + PhotoNuclearDQM. + + """ + def __init__(self, name, filter_class): + super().__init__(name, filter_class) + from LDMX.Biasing import include + include.library() + self.count_light_ions = True + + + def SingleNeutronFilter(): + """Configuration for keeping photonuclear events where a single neutron + carries most of the kinetic energy from the interaction + + + Returns + ------- + Instance of configured PhotoNuclearTopologyFilter + + """ + filter = PhotoNuclearTopologyFilter(name='SingleNeutron', + filter_class='biasing::SingleNeutronFilter' ) + filter.hard_particle_threshold = 200. + return filter + + def NothingHardFilter(): + """Configuration for keeping photonuclear events no particles received + significant kinetic energy. + + Returns + ------- + Instance of configured PhotoNuclearTopologyFilter + + """ + filter = PhotoNuclearTopologyFilter(name='NothingHard', + filter_class='biasing::NothingHardFilter') + filter.hard_particle_threshold = 200. + return filter + + + + + diff --git a/Biasing/src/Biasing/PhotoNuclearTopologyFilters.cxx b/Biasing/src/Biasing/PhotoNuclearTopologyFilters.cxx new file mode 100644 index 000000000..22ef7c4d5 --- /dev/null +++ b/Biasing/src/Biasing/PhotoNuclearTopologyFilters.cxx @@ -0,0 +1,79 @@ +#include "Biasing/PhotoNuclearTopologyFilters.h" + +namespace biasing { + +bool NothingHardFilter::rejectEvent( + const std::vector& secondaries) const { + for (const auto& secondary : secondaries) { + // Get the PDG ID of the track + const auto pdgID{ + std::abs(secondary->GetParticleDefinition()->GetPDGEncoding())}; + if (skipCountingParticle(pdgID)) { + continue; + } + auto energy{secondary->GetKineticEnergy()}; + if (energy > hard_particle_threshold_) { + return true; + } + } + return false; +} +bool SingleNeutronFilter::rejectEvent( + const std::vector& secondaries) const { + int hard_particles{0}; + int hard_neutrons{0}; + for (const auto& secondary : secondaries) { + // Get the PDG ID of the track + const auto pdgID{ + std::abs(secondary->GetParticleDefinition()->GetPDGEncoding())}; + if (skipCountingParticle(pdgID)) { + continue; + } + auto energy{secondary->GetKineticEnergy()}; + if (energy > hard_particle_threshold_) { + hard_particles++; + if (isNeutron(pdgID)) { + hard_neutrons++; + } + } + } + auto reject{hard_particles != hard_neutrons || hard_particles != 1}; + return reject; +} + +PhotoNuclearTopologyFilter::PhotoNuclearTopologyFilter( + const std::string& name, framework::config::Parameters& parameters) + : UserAction{name, parameters}, + hard_particle_threshold_{ + parameters.getParameter("hard_particle_threshold")}, + count_light_ions_{parameters.getParameter("count_light_ions")} {} + +void PhotoNuclearTopologyFilter::stepping(const G4Step* step) { + // Get the track associated with this step. + auto track{step->GetTrack()}; + + // Get the track info and check if this track has been tagged as the + // photon that underwent a photo-nuclear reaction. Only those tracks + // tagged as PN photos will be processed. The track is currently only + // tagged by the UserAction ECalProcessFilter which needs to be run + // before this UserAction. + auto trackInfo{simcore::UserTrackInformation::get(track)}; + if ((trackInfo != nullptr) && !trackInfo->isPNGamma()) return; + + // Get the PN photon daughters. + auto secondaries{step->GetSecondary()}; + + if (rejectEvent(*secondaries)) { + track->SetTrackStatus(fKillTrackAndSecondaries); + G4RunManager::GetRunManager()->AbortEvent(); + } + + // Once the PN gamma has been procesed, untag it so its not reprocessed + // again. + trackInfo->tagPNGamma(false); +} + +} // namespace biasing + +DECLARE_ACTION(biasing, NothingHardFilter) +DECLARE_ACTION(biasing, SingleNeutronFilter)