diff --git a/CMakeLists.txt b/CMakeLists.txt index 3de19365a..f8a3014c1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,7 +16,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR) find_package(cetmodules 3.20.00 REQUIRED) -project(sbncode VERSION 09.88.00.04 LANGUAGES CXX) +project(sbncode VERSION 09.89.01 LANGUAGES CXX) message(STATUS "\n\n ========================== ${PROJECT_NAME} ==========================") diff --git a/sbncode/DetSim/AdjustSimForTrigger_module.cc b/sbncode/DetSim/AdjustSimForTrigger_module.cc index 195615fee..14c7cabe2 100644 --- a/sbncode/DetSim/AdjustSimForTrigger_module.cc +++ b/sbncode/DetSim/AdjustSimForTrigger_module.cc @@ -17,9 +17,9 @@ #include "fhiclcpp/ParameterSet.h" #include "messagefacility/MessageLogger/MessageLogger.h" -// #include "lardataalg/DetectorInfo/DetectorClocksData.h" #include "lardataobj/RawData/OpDetWaveform.h" #include "lardataobj/RawData/TriggerData.h" +#include "lardataobj/Simulation/AuxDetSimChannel.h" #include "lardataobj/Simulation/BeamGateInfo.h" #include "lardataobj/Simulation/SimEnergyDeposit.h" #include "lardataobj/Simulation/SimPhotons.h" @@ -45,44 +45,52 @@ class AdjustSimForTrigger : public art::EDProducer { void produce(art::Event& e) override; private: - art::InputTag fInputSimEnergyDepositLabel; - art::InputTag fInputSimPhotonsLabel; art::InputTag fInputTriggerLabel; - art::InputTag fInputWaveformLabel; - art::InputTag fInputBeamGateInfoLabel; - double fAdditionalOffset; + art::InputTag fInitAuxDetSimChannelLabel; + art::InputTag fInitBeamGateInfoLabel; + art::InputTag fInitSimEnergyDepositLabel; + art::InputTag fInitSimPhotonsLabel; + art::InputTag fInitWaveformLabel; + bool fShiftAuxDetIDEs; + bool fShiftBeamGateInfo; bool fShiftSimEnergyDeposits; bool fShiftSimPhotons; bool fShiftWaveforms; - bool fShiftBeamGateInfo; + double fAdditionalOffset; + static constexpr auto& kModuleName = "AdjustSimForTrigger"; }; AdjustSimForTrigger::AdjustSimForTrigger(fhicl::ParameterSet const& p) : EDProducer{p} - , fInputSimEnergyDepositLabel{p.get("InputSimEnergyDepositLabel", "")} - , fInputSimPhotonsLabel{p.get("InputSimPhotonsLabel", "")} - , fInputTriggerLabel{p.get("InputTriggerLabel", "")} - , fInputWaveformLabel(p.get("InputWaveformLabel", "")) - , fInputBeamGateInfoLabel{p.get("InputBeamGateInfoLabel", "")} - , fAdditionalOffset{p.get("AdditionalOffset")} + , fInputTriggerLabel{p.get("InputTriggerLabel", "undefined")} + , fInitAuxDetSimChannelLabel(p.get("InitAuxDetSimChannelLabel", "undefined")) + , fInitBeamGateInfoLabel{p.get("InitBeamGateInfoLabel", "undefined")} + , fInitSimEnergyDepositLabel{p.get("InitSimEnergyDepositLabel", "undefined")} + , fInitSimPhotonsLabel{p.get("InitSimPhotonsLabel", "undefined")} + , fInitWaveformLabel(p.get("InitWaveformLabel", "undefined")) + , fShiftAuxDetIDEs{p.get("ShiftAuxDetIDEs", false)} + , fShiftBeamGateInfo{p.get("ShiftBeamGateInfo", false)} , fShiftSimEnergyDeposits{p.get("ShiftSimEnergyDeposits", false)} , fShiftSimPhotons{p.get("ShiftSimPhotons", false)} , fShiftWaveforms{p.get("ShiftWaveforms", false)} - , fShiftBeamGateInfo{p.get("ShiftBeamGateInfo", false)} + , fAdditionalOffset{p.get("AdditionalOffset", 0.)} { - if (!(fShiftSimEnergyDeposits | fShiftSimPhotons | fShiftWaveforms | fShiftBeamGateInfo)) { + if (!(fShiftSimEnergyDeposits || fShiftSimPhotons || fShiftWaveforms || fShiftAuxDetIDEs || + fShiftBeamGateInfo)) { throw art::Exception(art::errors::EventProcessorFailure) - << "NO SHIFTS ENABLED!\n" - << "SHIFTING SIMENERGYDEPOSITS? " << fShiftSimEnergyDeposits << '\n' - << "SHIFTING SIMPHOTONS? " << fShiftSimPhotons << '\n' - << "SHIFTING OPDETWAVEFORMS? " << fShiftWaveforms << '\n' - << "SHIFTING BEAMGATEINFO? " << fShiftBeamGateInfo << '\n'; + << kModuleName << ": NO SHIFTS ENABLED!\n"; } + mf::LogInfo(kModuleName) << std::boolalpha << "SHIFTING AUXDETIDES? " << fShiftAuxDetIDEs << '\n' + << "SHIFTING BEAMGATEINFO? " << fShiftBeamGateInfo << '\n' + << "SHIFTING SIMENERGYDEPOSITS? " << fShiftSimEnergyDeposits << '\n' + << "SHIFTING SIMPHOTONS? " << fShiftSimPhotons << '\n' + << "SHIFTING OPDETWAVEFORMS? " << fShiftWaveforms; + if (fShiftAuxDetIDEs) { produces>(); } + if (fShiftBeamGateInfo) { produces>(); } if (fShiftSimEnergyDeposits) { produces>(); } if (fShiftSimPhotons) { produces>(); } if (fShiftWaveforms) { produces>(); } - if (fShiftBeamGateInfo) { produces>(); } } void AdjustSimForTrigger::produce(art::Event& e) @@ -92,10 +100,11 @@ void AdjustSimForTrigger::produce(art::Event& e) if (triggers.size() != 1) { if (triggers.empty()) { - throw art::Exception(art::errors::EventProcessorFailure) << "NO TRIGGER IDENTIFIED!\n"; + throw art::Exception(art::errors::EventProcessorFailure) + << kModuleName << ": NO TRIGGER IDENTIFIED!"; } throw art::Exception(art::errors::EventProcessorFailure) - << "MORE THAN ONE TRIGGER IN EVENT... why?\n"; + << kModuleName << ": MORE THAN ONE TRIGGER IN EVENT... why?"; } // Assuming there is a trigger, get time shift @@ -111,13 +120,59 @@ void AdjustSimForTrigger::produce(art::Event& e) hasValidTriggerTime ? clock_data.TriggerTime() - trigger.TriggerTime() + fAdditionalOffset : 0.; const double timeShiftForTrigger_ns = 1000. * timeShiftForTrigger_us; - mf::LogInfo("AdjustSimForTrigger") - << "FOR THIS EVENT THE TIME SHIFT BEING ASSUMED IS " << timeShiftForTrigger_ns << " ns ...\n"; + mf::LogInfo(kModuleName) << "FOR THIS EVENT THE TIME SHIFT BEING ASSUMED IS " + << timeShiftForTrigger_ns << " ns ..."; + + // Loop over the sim::AuxDetIDE and shift time BACK by the TRIGGER + if (fShiftAuxDetIDEs) { + auto const& simChannels = + e.getProduct>(fInitAuxDetSimChannelLabel); + auto pSimChannels = std::make_unique>(); + + pSimChannels->reserve(simChannels.size()); + + for (auto const& simChannel : simChannels) { + std::vector shiftedAuxDetIDEs = simChannel.AuxDetIDEs(); + for (auto& auxDetIDE : shiftedAuxDetIDEs) { + auxDetIDE.entryT += timeShiftForTrigger_ns; + auxDetIDE.exitT += timeShiftForTrigger_ns; + } + pSimChannels->emplace_back(sim::AuxDetSimChannel( + simChannel.AuxDetID(), shiftedAuxDetIDEs, simChannel.AuxDetSensitiveID())); + } + e.put(std::move(pSimChannels)); + } + + // Repeat for sim::BeamGateInfo + if (fShiftBeamGateInfo) { + auto const& beamGates = e.getProduct>(fInitBeamGateInfoLabel); + + if (beamGates.size() != 1) { + if (beamGates.empty()) { + throw art::Exception(art::errors::EventProcessorFailure) + << kModuleName << ": THERE IS NO BEAM GATE INFO!\n"; + } + throw art::Exception(art::errors::EventProcessorFailure) + << kModuleName << ": MORE THAN ONE BEAM GATE?\n"; + } + + const auto& beamGate = beamGates[0]; + + const double shiftedBeamGateStart = beamGate.Start() + timeShiftForTrigger_ns; + const double gateWidth = beamGate.Width(); + const sim::BeamType_t beam = beamGate.BeamType(); + + const sim::BeamGateInfo shiftedBeamGate(shiftedBeamGateStart, gateWidth, beam); + + auto pBeamGateInfos = std::make_unique>(1, shiftedBeamGate); + + e.put(std::move(pBeamGateInfos)); + } - // Loop over the SimEnergyDeposit objects and shift time BACK by the TRIGGER + // Repeat for sim::SimEnergyDeposit if (fShiftSimEnergyDeposits) { auto const& simEDeps = - e.getProduct>(fInputSimEnergyDepositLabel); + e.getProduct>(fInitSimEnergyDepositLabel); auto pSimEDeps = std::make_unique>(); pSimEDeps->reserve(simEDeps.size()); @@ -136,24 +191,24 @@ void AdjustSimForTrigger::produce(art::Event& e) const int thisPDG = inSimEDep.PdgCode(); const int origID = inSimEDep.OrigTrackID(); - pSimEDeps->push_back(sim::SimEnergyDeposit(numphotons, - numelectrons, - syratio, - energy, - start, - end, - startT, - endT, - thisID, - thisPDG, - origID)); + pSimEDeps->emplace_back(sim::SimEnergyDeposit(numphotons, + numelectrons, + syratio, + energy, + start, + end, + startT, + endT, + thisID, + thisPDG, + origID)); } e.put(std::move(pSimEDeps)); } // Repeat for sim::SimPhotons if (fShiftSimPhotons) { - auto const& simPhotons = e.getProduct>(fInputSimPhotonsLabel); + auto const& simPhotons = e.getProduct>(fInitSimPhotonsLabel); auto pSimPhotonss = std::make_unique>(simPhotons); for (auto& photons : *pSimPhotonss) { @@ -166,7 +221,7 @@ void AdjustSimForTrigger::produce(art::Event& e) // Repeat for raw::OpDetWaveform if (fShiftWaveforms) { - auto const& waveforms = e.getProduct>(fInputWaveformLabel); + auto const& waveforms = e.getProduct>(fInitWaveformLabel); auto pWaveforms = std::make_unique>(waveforms); for (auto& waveform : *pWaveforms) { @@ -174,30 +229,6 @@ void AdjustSimForTrigger::produce(art::Event& e) } e.put(std::move(pWaveforms)); } - - // Repeat for sim::BeamGateInfo - if (fShiftBeamGateInfo) { - auto const& beamGates = e.getProduct>(fInputBeamGateInfoLabel); - - if (beamGates.size() != 1) { - if (beamGates.empty()) { - throw art::Exception(art::errors::EventProcessorFailure) << "THERE IS NO BEAM GATE INFO!\n"; - } - throw art::Exception(art::errors::EventProcessorFailure) << "MORE THAN ONE BEAM GATE?\n"; - } - - const auto& beamGate = beamGates[0]; - - const double shiftedBeamGateStart = beamGate.Start() + timeShiftForTrigger_ns; - const double gateWidth = beamGate.Width(); - const sim::BeamType_t beam = beamGate.BeamType(); - - const sim::BeamGateInfo shiftedBeamGate(shiftedBeamGateStart, gateWidth, beam); - - auto pBeamGateInfos = std::make_unique>(1, shiftedBeamGate); - - e.put(std::move(pBeamGateInfos)); - } } DEFINE_ART_MODULE(AdjustSimForTrigger) diff --git a/sbncode/DetSim/CMakeLists.txt b/sbncode/DetSim/CMakeLists.txt index 2a81d6a62..0bf4e263d 100644 --- a/sbncode/DetSim/CMakeLists.txt +++ b/sbncode/DetSim/CMakeLists.txt @@ -28,10 +28,14 @@ set( MODULE_LIBRARIES ROOT::Geom ROOT::XMLIO ROOT::Gdml + ROOT::Core FFTW3::FFTW3 ) 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..c70075f07 --- /dev/null +++ b/sbncode/DetSim/FilterSimEnergyDeposits_module.cc @@ -0,0 +1,166 @@ +//////////////////////////////////////////////////////////////////////// +// Class: FilterSimEnergyDeposits +// Plugin Type: producer (Unknown Unknown) +// File: FilterSimEnergyDeposits_module.cc +// +// Generated at Thu Oct 17 20:16:12 2024 by Jacob Zettlemoyer using cetskelgen +// from cetlib version 3.18.02. +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.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 "lardata/DetectorInfoServices/LArPropertiesService.h" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" + +#include "lardataobj/Simulation/SimEnergyDeposit.h" + +#include "TVector3.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; + std::string fOutSimEnergyDepositLabel; + double fP1_x; + double fP1_y; + double fP1_z; + double fP2_x; + double fP2_y; + double fP2_z; + static constexpr auto& kModuleName = "FilterSimEnergyDeposit"; + + bool isWithinVolume(TVector3 p1, TVector3 p2, TVector3 sedep); + double ShiftX(double z); + +}; + + +FilterSimEnergyDeposits::FilterSimEnergyDeposits(fhicl::ParameterSet const& p) + : EDProducer{p} + , fInitSimEnergyDepositLabel{p.get("InitSimEnergyDepositLabel", "undefined")} + , fOutSimEnergyDepositLabel{p.get("OutSimEnergyDepositLabel", "undefined")} + , fP1_x{p.get("P1_X", 0.0)} + , fP1_y{p.get("P1_Y", 0.0)} + , fP1_z{p.get("P1_Z", 0.0)} + , fP2_x{p.get("P2_X", 0.0)} + , fP2_y{p.get("P2_Y", 0.0)} + , fP2_z{p.get("P2_Z", 0.0)} +{ + // Call appropriate produces<>() functions here. + // Call appropriate consumes<>() for any products to be retrieved by this module. + produces>(); +} + +bool FilterSimEnergyDeposits::isWithinVolume(TVector3 p_min, TVector3 p_max, TVector3 edep) +{ + bool ingap = false; + if(edep.X() >= p_min.X() && edep.X() <= p_max.X() && + edep.Y() >= p_min.Y() && edep.Y() <= p_max.Y() && + edep.Z() >= p_min.Z() && edep.Z() <= p_max.Z()) + { + ingap = true; + } + return ingap; +} + +double FilterSimEnergyDeposits::ShiftX(double z) +{ + //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 = -0.431172; + double b = -2.52724; + double c = 0.22043; + 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()); + + TVector3 p1(fP1_x,fP1_y,fP1_z); + TVector3 p2(fP2_x,fP2_y,fP2_z); + + for(auto const& sedep : simEDeps) + { + TVector3 dep(sedep.StartX(), sedep.StartY(), sedep.StartZ()); + if(!(isWithinVolume(p1, p2, dep))) + { + art::ServiceHandle fGeometry; + geo::TPCID tpcid = fGeometry->PositionToTPCID(sedep.MidPoint()); + short int sign = 1; + //if not in TPC, don't do anything to be sure + if(bool(tpcid)) + { + const geo::TPCGeo& tpcGeo = fGeometry->TPC(tpcid); + sign = -1*tpcGeo.DetectDriftDirection(); //this gives us direction towards the anode, I want to shift towards the cathode so take the inverse of what I get + } + else + { + std::cout << "TPC ID not found! Not performing a shift!" << std::endl; + sign = 0; + } + const int numphotons = sedep.NumPhotons(); + const int numelectrons = sedep.NumElectrons(); + const double syratio = sedep.ScintYieldRatio(); + const double energy = sedep.Energy(); + //need to shift in -x for x positions > cathode position and +x for x positions < cathode position in each TPC + double shift_start_x = sedep.Start().X() + sign*ShiftX(std::abs(sedep.Start().Z())); + const geo::Point_t start = {shift_start_x, sedep.Start().Y(), sedep.Start().Z()}; + double shift_end_x = sedep.End().X() + sign*ShiftX(std::abs(sedep.End().Z())); + const geo::Point_t end = {shift_end_x, sedep.End().Y(), 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->emplace_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..e95961ee1 --- /dev/null +++ b/sbncode/DetSim/fcl/icarus_simedepfilter.fcl @@ -0,0 +1,43 @@ +BEGIN_PROLOG + +simedepfilter_ind1gap: { + module_type: "FilterSimEnergyDeposits" + InitSimEnergyDepositLabel: "ionization" + OutSimEnergyDepositLabel: "filteredSED" + P1_X: -400 + P1_Y: -200 + P1_Z: -2 + P2_X: 400 + P2_Y: 200 + P2_Z: 2 +} + +simedepfilter_cathodewest: { + module_type: "FilterSimEnergyDeposits" + InitSimEnergyDepositLabel: "ionization" + OutSimEnergyDepositLabel: "filteredSED" + P1_X: 192.2 + P1_Y: -200 + P1_Z: -900 + P2_X: 232.2 + P2_Y: 200 + P2_Z: 900 +} + +simedepfilter_cathodeeast: { + module_type: "FilterSimEnergyDeposits" + InitSimEnergyDepositLabel: "ionization" + OutSimEnergyDepositLabel: "filteredSED" + P1_X: -232.2 + P1_Y: -200 + P1_Z: -900 + P2_X: -198.2 + P2_Y: 200 + P2_Z: 900 +} + +filter_ind1gap: @local::simedepfilter_ind1gap +filter_cathodewest: @local::simedepfilter_cathodewest +filter_cathodeeast: @local::simedepfilter_cathodeeast + +END_PROLOG \ No newline at end of file diff --git a/sbncode/EventGenerator/MultiPart/multipartvertex_sbnd.fcl b/sbncode/EventGenerator/MultiPart/multipartvertex_sbnd.fcl deleted file mode 100644 index 8c03a1476..000000000 --- a/sbncode/EventGenerator/MultiPart/multipartvertex_sbnd.fcl +++ /dev/null @@ -1,43 +0,0 @@ -BEGIN_PROLOG - -MultiPartVertex : { - module_type : "MultiPartVertex" - DebugMode : 0 - G4Time : 1700 - G4TimeJitter : 0 - TPCRange : [[0,0],[0,1]] - XRange : [10] - YRange : [10] - ZRange : [10] - MultiMax : 6 - MultiMin : 1 - ParticleParameter: { - PDGCode : [ [11,13], [22], [211,-211], [2212]] - MinMulti : [ 0, 0, 0, 0] - MaxMulti : [ 1, 4, 3, 3] - ProbWeight : [ 3, 2, 3, 3] - KERange : [ [40,1500], [40,1500], [40,1500], [50,700]] - } -} - -MultiPartRain : { - module_type : "MultiPartRain" - DebugMode : 0 - G4Time : 0 - G4TimeJitter : 0 - TPCRange : [[0,0],[0,1]] - XRange : [10] - YRange : [10] - ZRange : [10] - MultiMax : 6 - MultiMin : 1 - ParticleParameter: { - PDGCode : [ [11,13], [22], [211,-211], [2212]] - MinMulti : [ 0, 0, 0, 0] - MaxMulti : [ 1, 4, 3, 3] - ProbWeight : [ 3, 2, 3, 3] - KERange : [ [40,1500], [40,1500], [40,1500], [50,700]] - } -} - -END_PROLOG diff --git a/ups/product_deps b/ups/product_deps index 9e73859bb..199bb72ff 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -254,7 +254,7 @@ libdir fq_dir lib product version qual flags genie_xsec v3_04_00 - larcv2 v2_2_6 - -larsoft v09_88_00 - +larsoft v09_89_01 - sbnanaobj v09_23_00_01 - sbndaq_artdaq_core v1_09_00of1 - sbndata v01_07 -