|
| 1 | +/** |
| 2 | + * @file sbncode/DetSim/FilterSimEnergyDeposits_module.cc |
| 3 | + * @date October 17, 2024 |
| 4 | + * @author Jacob Zettlemoyer |
| 5 | + */ |
| 6 | + |
| 7 | +#include "art/Framework/Core/EDProducer.h" |
| 8 | +#include "art/Framework/Core/ModuleMacros.h" |
| 9 | +#include "art/Framework/Principal/Event.h" |
| 10 | +#include "canvas/Utilities/InputTag.h" |
| 11 | +#include "fhiclcpp/ParameterSet.h" |
| 12 | +#include "messagefacility/MessageLogger/MessageLogger.h" |
| 13 | +#include "larcore/Geometry/Geometry.h" |
| 14 | +#include "larcore/Geometry/WireReadout.h" |
| 15 | +#include "larcorealg/Geometry/GeometryCore.h" |
| 16 | +#include "larcorealg/Geometry/WireReadoutGeom.h" |
| 17 | +#include "larcore/CoreUtils/ServiceUtil.h" |
| 18 | +#include "larcorealg/Geometry/BoxBoundedGeo.h" |
| 19 | + |
| 20 | +#include "lardataobj/Simulation/SimEnergyDeposit.h" |
| 21 | + |
| 22 | +#include <memory> |
| 23 | + |
| 24 | +class FilterSimEnergyDeposits; |
| 25 | + |
| 26 | + |
| 27 | +class FilterSimEnergyDeposits : public art::EDProducer { |
| 28 | +public: |
| 29 | + explicit FilterSimEnergyDeposits(fhicl::ParameterSet const& p); |
| 30 | + // The compiler-generated destructor is fine for non-base |
| 31 | + // classes without bare pointers or other resource use. |
| 32 | + |
| 33 | + // Plugins should not be copied or assigned. |
| 34 | + FilterSimEnergyDeposits(FilterSimEnergyDeposits const&) = delete; |
| 35 | + FilterSimEnergyDeposits(FilterSimEnergyDeposits&&) = delete; |
| 36 | + FilterSimEnergyDeposits& operator=(FilterSimEnergyDeposits const&) = delete; |
| 37 | + FilterSimEnergyDeposits& operator=(FilterSimEnergyDeposits&&) = delete; |
| 38 | + // Required functions. |
| 39 | + void produce(art::Event& e) override; |
| 40 | + |
| 41 | +private: |
| 42 | + |
| 43 | + // Declare member data here. |
| 44 | + |
| 45 | + art::InputTag fInitSimEnergyDepositLabel; |
| 46 | + geo::BoxBoundedGeo fBox; |
| 47 | + static constexpr auto kModuleName = "FilterSimEnergyDeposit"; |
| 48 | + double ShiftX(double z) const; |
| 49 | + double fA; |
| 50 | + double fB; |
| 51 | + double fC; |
| 52 | + |
| 53 | +}; |
| 54 | + |
| 55 | + |
| 56 | +FilterSimEnergyDeposits::FilterSimEnergyDeposits(fhicl::ParameterSet const& p) |
| 57 | + : EDProducer{p} |
| 58 | + , fInitSimEnergyDepositLabel{p. get<art::InputTag>( "InitSimEnergyDepositLabel", "undefined")} |
| 59 | + , fBox{p.get<double>("P1_X"), p.get<double>("P2_X"), |
| 60 | + p.get<double>("P1_Y"), p.get<double>("P2_Y"), |
| 61 | + p.get<double>("P1_Z"), p.get<double>("P2_Z")} |
| 62 | + , fA{p.get<double>("A")} |
| 63 | + , fB{p.get<double>("B")} |
| 64 | + , fC{p.get<double>("C")} |
| 65 | +{ |
| 66 | + // Call appropriate produces<>() functions here. |
| 67 | + // Call appropriate consumes<>() for any products to be retrieved by this module. |
| 68 | + produces<std::vector<sim::SimEnergyDeposit>>(); |
| 69 | +} |
| 70 | + |
| 71 | + |
| 72 | +double FilterSimEnergyDeposits::ShiftX(double z) const |
| 73 | +{ |
| 74 | + //known as a "Hill equation" from fitting distribution of angle corrected tracks looking at the deflection as it gets closer to z = 0 |
| 75 | + double a = fA; |
| 76 | + double b = fB; |
| 77 | + double c = fC; |
| 78 | + double num = a*std::pow(z, b); |
| 79 | + double denom = c + std::pow(z, b); |
| 80 | + return num/denom; |
| 81 | +} |
| 82 | + |
| 83 | +void FilterSimEnergyDeposits::produce(art::Event& e) |
| 84 | +{ |
| 85 | + auto const& simEDeps = |
| 86 | + e.getProduct<std::vector<sim::SimEnergyDeposit>>(fInitSimEnergyDepositLabel); |
| 87 | + auto pSimEDeps = std::make_unique<std::vector<sim::SimEnergyDeposit>>(); |
| 88 | + pSimEDeps->reserve(simEDeps.size()); |
| 89 | + |
| 90 | + for(auto const& sedep : simEDeps) |
| 91 | + { |
| 92 | + if(fBox.ContainsPosition(sedep.Start())) |
| 93 | + continue; |
| 94 | + art::ServiceHandle<geo::Geometry const> geom; |
| 95 | + geo::TPCGeo const* TPC = geom->PositionToTPCptr(sedep.MidPoint()); |
| 96 | + if(!TPC) { |
| 97 | + mf::LogVerbatim(kModuleName) << "TPC ID not found! Not performing a shift!"; |
| 98 | + } |
| 99 | + const int numphotons = sedep.NumPhotons(); |
| 100 | + const int numelectrons = sedep.NumElectrons(); |
| 101 | + const double syratio = sedep.ScintYieldRatio(); |
| 102 | + const double energy = sedep.Energy(); |
| 103 | + geo::Point_t start = sedep.Start(); |
| 104 | + geo::Point_t end = sedep.End(); |
| 105 | + if (TPC) { |
| 106 | + auto const& wireReadout = art::ServiceHandle<geo::WireReadout const>()->Get(); |
| 107 | + auto const& referencePlane = wireReadout.FirstPlane(TPC->ID()); |
| 108 | + referencePlane.DriftPoint(start, ShiftX(std::abs(sedep.Start().Z()))); |
| 109 | + referencePlane.DriftPoint(end, ShiftX(std::abs(sedep.End().Z()))); |
| 110 | + } |
| 111 | + const double startT = sedep.StartT(); |
| 112 | + const double endT = sedep.EndT(); |
| 113 | + const int thisID = sedep.TrackID(); |
| 114 | + const int thisPDG = sedep.PdgCode(); |
| 115 | + const int origID = sedep.OrigTrackID(); |
| 116 | + pSimEDeps->push_back(sim::SimEnergyDeposit(numphotons, |
| 117 | + numelectrons, |
| 118 | + syratio, |
| 119 | + energy, |
| 120 | + start, |
| 121 | + end, |
| 122 | + startT, |
| 123 | + endT, |
| 124 | + thisID, |
| 125 | + thisPDG, |
| 126 | + origID)); |
| 127 | + } |
| 128 | + e.put(std::move(pSimEDeps)); |
| 129 | +} |
| 130 | + |
| 131 | +DEFINE_ART_MODULE(FilterSimEnergyDeposits) |
0 commit comments