diff --git a/sbncode/DetSim/CMakeLists.txt b/sbncode/DetSim/CMakeLists.txt index 2a81d6a62..113aa375a 100644 --- a/sbncode/DetSim/CMakeLists.txt +++ b/sbncode/DetSim/CMakeLists.txt @@ -32,6 +32,9 @@ set( MODULE_LIBRARIES ) cet_build_plugin(AdjustSimForTrigger art::module LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(FilterSimEnergyDeposits art::module LIBRARIES ${MODULE_LIBRARIES}) + +add_subdirectory(fcl) #install_headers() install_fhicl() diff --git a/sbncode/DetSim/FilterSimEnergyDeposits_module.cc b/sbncode/DetSim/FilterSimEnergyDeposits_module.cc new file mode 100644 index 000000000..802887c91 --- /dev/null +++ b/sbncode/DetSim/FilterSimEnergyDeposits_module.cc @@ -0,0 +1,125 @@ +/** + * @file sbncode/DetSim/FilterSimEnergyDeposits_module.cc + * @date October 17, 2024 + * @author Jacob Zettlemoyer + */ + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "larcore/Geometry/Geometry.h" +#include "larcorealg/Geometry/GeometryCore.h" +#include "larcore/CoreUtils/ServiceUtil.h" +#include "larcorealg/Geometry/BoxBoundedGeo.h" + +#include "lardataobj/Simulation/SimEnergyDeposit.h" + +#include + +class FilterSimEnergyDeposits; + + +class FilterSimEnergyDeposits : public art::EDProducer { +public: + explicit FilterSimEnergyDeposits(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + FilterSimEnergyDeposits(FilterSimEnergyDeposits const&) = delete; + FilterSimEnergyDeposits(FilterSimEnergyDeposits&&) = delete; + FilterSimEnergyDeposits& operator=(FilterSimEnergyDeposits const&) = delete; + FilterSimEnergyDeposits& operator=(FilterSimEnergyDeposits&&) = delete; + // Required functions. + void produce(art::Event& e) override; + +private: + + // Declare member data here. + + art::InputTag fInitSimEnergyDepositLabel; + geo::BoxBoundedGeo fBox; + static constexpr auto kModuleName = "FilterSimEnergyDeposit"; + double ShiftX(double z) const; + double fA; + double fB; + double fC; + +}; + + +FilterSimEnergyDeposits::FilterSimEnergyDeposits(fhicl::ParameterSet const& p) + : EDProducer{p} + , fInitSimEnergyDepositLabel{p.get("InitSimEnergyDepositLabel", "undefined")} + , fBox{p.get("P1_X"), p.get("P2_X"), + p.get("P1_Y"), p.get("P2_Y"), + p.get("P1_Z"), p.get("P2_Z")} + , fA{p.get("A")} + , fB{p.get("B")} + , fC{p.get("C")} +{ + // Call appropriate produces<>() functions here. + // Call appropriate consumes<>() for any products to be retrieved by this module. + produces>(); +} + + +double FilterSimEnergyDeposits::ShiftX(double z) const +{ + //known as a "Hill equation" from fitting distribution of angle corrected tracks looking at the deflection as it gets closer to z = 0 + double a = fA; + double b = fB; + double c = fC; + double num = a*std::pow(z, b); + double denom = c + std::pow(z, b); + return num/denom; +} + +void FilterSimEnergyDeposits::produce(art::Event& e) +{ + auto const& simEDeps = + e.getProduct>(fInitSimEnergyDepositLabel); + auto pSimEDeps = std::make_unique>(); + pSimEDeps->reserve(simEDeps.size()); + + for(auto const& sedep : simEDeps) + { + if(fBox.ContainsPosition(sedep.Start())) + continue; + art::ServiceHandle geom; + geo::TPCGeo const* TPC = geom->PositionToTPCptr(sedep.MidPoint()); + if(!TPC) { + mf::LogVerbatim(kModuleName) << "TPC ID not found! Not performing a shift!"; + } + const int numphotons = sedep.NumPhotons(); + const int numelectrons = sedep.NumElectrons(); + const double syratio = sedep.ScintYieldRatio(); + const double energy = sedep.Energy(); + geo::Point_t start = sedep.Start(); + if (TPC) TPC->DriftPoint(start, ShiftX(std::abs(sedep.Start().Z()))); + geo::Point_t end = sedep.End(); + if (TPC) TPC->DriftPoint(end, ShiftX(std::abs(sedep.End().Z()))); + const double startT = sedep.StartT(); + const double endT = sedep.EndT(); + const int thisID = sedep.TrackID(); + const int thisPDG = sedep.PdgCode(); + const int origID = sedep.OrigTrackID(); + pSimEDeps->push_back(sim::SimEnergyDeposit(numphotons, + numelectrons, + syratio, + energy, + start, + end, + startT, + endT, + thisID, + thisPDG, + origID)); + } + e.put(std::move(pSimEDeps)); +} + +DEFINE_ART_MODULE(FilterSimEnergyDeposits) diff --git a/sbncode/DetSim/fcl/CMakeLists.txt b/sbncode/DetSim/fcl/CMakeLists.txt new file mode 100644 index 000000000..eca2a52d6 --- /dev/null +++ b/sbncode/DetSim/fcl/CMakeLists.txt @@ -0,0 +1 @@ +install_fhicl() \ No newline at end of file diff --git a/sbncode/DetSim/fcl/icarus_simedepfilter.fcl b/sbncode/DetSim/fcl/icarus_simedepfilter.fcl new file mode 100644 index 000000000..124a50981 --- /dev/null +++ b/sbncode/DetSim/fcl/icarus_simedepfilter.fcl @@ -0,0 +1,19 @@ +BEGIN_PROLOG + +simedepfilter_ind1gap: { + module_type: "FilterSimEnergyDeposits" + InitSimEnergyDepositLabel: "ionization" + P1_X: -400 #cm, create a box that spans the full x,y range. Not super-precise. + P1_Y: -200 #cm, create a box that spans the full x,y range. Not super-precise. + P1_Z: -2 #cm, one edge of the filtered gap based on the 37 mm total size of the wire gap + P2_X: 400 #cm, create a box that spans the full x,y range. Not super-precise. + P2_Y: 200 #cm, create a box that spans the full x,y range. Not super-precise. + P2_Z: 2 #cm, other edge of the filtered gap based on the 37 mm total size of the wire gap + #From SBN doc-db 38701: + A: -0.431172 + B: -2.52724 + C: 0.22043 +} + + +END_PROLOG \ No newline at end of file