From 7f44e94300f1cff97fb2fd5f3b5511be69e7fb2b Mon Sep 17 00:00:00 2001 From: mantheys Date: Wed, 17 Dec 2025 12:08:46 +0100 Subject: [PATCH 1/9] ADD MAX PE CUT FOR FLASH MATCHING SELECTION CRITERIA --- dunereco/LowEUtils/LowEUtils.cc | 60 ++++++++++++++++++++++++++++++++- dunereco/LowEUtils/LowEUtils.h | 17 ++++++++++ 2 files changed, 76 insertions(+), 1 deletion(-) diff --git a/dunereco/LowEUtils/LowEUtils.cc b/dunereco/LowEUtils/LowEUtils.cc index 0ad2adb4..44d561a2 100644 --- a/dunereco/LowEUtils/LowEUtils.cc +++ b/dunereco/LowEUtils/LowEUtils.cc @@ -23,13 +23,18 @@ namespace lowe fAdjOpFlashY(p.get("AdjOpFlashY", 100.0)), // Y coordinate for flash projection [cm] fAdjOpFlashZ(p.get("AdjOpFlashZ", 100.0)), // Z coordinate for flash projection [cm] fAdjOpFlashMinPECut(p.get("AdjOpFlashMinPECut", 0.0)), // Minimum PE for flash selection + fAdjOpFlashMaxPECut(p.get("AdjOpFlashMaxPECut", 1e10)), // Maximum PE for flash selection fAdjOpFlashMaxPERatioCut(p.get("AdjOpFlashMaxPERatioCut", 1.0)), // Maximum PE ratio for flash selection fAdjOpFlashMembraneProjection(p.get("AdjOpFlashMembraneProjection", false)), // Whether to project flashes onto the membrane fAdjOpFlashEndCapProjection(p.get("AdjOpFlashEndCapProjection", false)), // Whether to project flashes onto the end cap fAdjOpFlashMinPEAttenuation(p.get("AdjOpFlashMinPEAttenuation", 1)), // PE cut attenuation over drift time [us] to compensate for light attenuation + fAdjOpFlashMaxPEAttenuation(p.get("AdjOpFlashMaxPEAttenuation", 1)), // PE cut attenuation over drift time [us] to compensate for light attenuation fAdjOpFlashMinPEAttenuate(p.get("AdjOpFlashMinPEAttenuate", "flat")), // Type of attenuation for minimum PE cut ("light_map", "asymptotic", "linear" or "flat") + fAdjOpFlashMaxPEAttenuate(p.get("AdjOpFlashMaxPEAttenuate", "flat")), // Type of attenuation for maximum PE cut ("light_map", "asymptotic", "linear" or "flat") fAdjOpFlashMinPEAttenuationStrength(p.get("AdjOpFlashMinPEAttenuationStrength", 10)), // Strength of the asymptotic attenuation for minimum PE cut (in powers of 10) + fAdjOpFlashMaxPEAttenuationStrength(p.get("AdjOpFlashMaxPEAttenuationStrength", 10)), // Strength of the asymptotic attenuation for maximum PE cut (in powers of 10) fAdjOpFlashMinPELightMap(p.get>>>("AdjOpFlashMinPELightMap", {})), // Light map file and histogram name for light map attenuation + fAdjOpFlashMaxPELightMap(p.get>>>("AdjOpFlashMaxPELightMap", {})), // Light map file and histogram name for light map attenuation producer(new ProducerUtils(p)) { // Initialize the LowEUtils instance @@ -1421,6 +1426,16 @@ namespace lowe const float &MatchedDriftTime, const float &ClusterCharge, const float &OpFlashPE) + { + return CutPDSFlashMinPE(TPCDriftTime, MatchedDriftTime, ClusterCharge, OpFlashPE) && + CutPDSFlashMaxPE(TPCDriftTime, MatchedDriftTime, ClusterCharge, OpFlashPE); + } // SelectPDSFlashPE + + bool LowEUtils::CutPDSFlashMinPE( + const float &TPCDriftTime, + const float &MatchedDriftTime, + const float &ClusterCharge, + const float &OpFlashPE) { if (fAdjOpFlashMinPEAttenuate == "light_map") { double a, b, c; @@ -1457,5 +1472,48 @@ namespace lowe else { return true; } } return true; - } // SelectPDSFlashPE + } + + bool LowEUtils::CutPDSFlashMaxPE( + const float &TPCDriftTime, + const float &MatchedDriftTime, + const float &ClusterCharge, + const float &OpFlashPE) + { + if (fAdjOpFlashMaxPEAttenuate == "light_map") { + double a, b, c; + for (int i = 0; i < int(fAdjOpFlashMaxPELightMap.size()); i++) { + if (fAdjOpFlashMaxPELightMap[i].first == "Amplitude") { + a = fAdjOpFlashMaxPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMaxPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMaxPELightMap[i].second[2]; + } + else if (fAdjOpFlashMaxPELightMap[i].first == "Attenuation") { + b = fAdjOpFlashMaxPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMaxPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMaxPELightMap[i].second[2]; + } + else if (fAdjOpFlashMaxPELightMap[i].first == "Correction") { + c = fAdjOpFlashMaxPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMaxPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMaxPELightMap[i].second[2]; + } + } + double MaxPE = pow(10, a - a * b * MatchedDriftTime / TPCDriftTime + c * pow(MatchedDriftTime / TPCDriftTime, 2)); + if (OpFlashPE > MaxPE) { return false; } + else { return true; } + } + else if (fAdjOpFlashMaxPEAttenuate == "asymptotic") { + if (OpFlashPE > fAdjOpFlashMaxPECut - fAdjOpFlashMaxPEAttenuation * (1-pow(pow(10,fAdjOpFlashMaxPEAttenuationStrength), -MatchedDriftTime / TPCDriftTime)) * fAdjOpFlashMaxPECut) { return false; } + else { return true; } + } + else if (fAdjOpFlashMaxPEAttenuate == "linear") { + if (OpFlashPE > fAdjOpFlashMaxPECut - fAdjOpFlashMaxPEAttenuation * pow(MatchedDriftTime / TPCDriftTime, 2) * fAdjOpFlashMaxPECut) { return false; } + else { return true; } + } + else if (fAdjOpFlashMaxPEAttenuate == "flat") { + if (OpFlashPE > fAdjOpFlashMaxPECut) { return false; } + else { return true; } + } + else { + ProducerUtils::PrintInColor("Warning: Unknown fAdjOpFlashMaxPEAttenuate option " + fAdjOpFlashMaxPEAttenuate + ", using flat cut", ProducerUtils::GetColor("red"), "Warning"); + if (OpFlashPE > fAdjOpFlashMaxPECut) { return false; } + else { return true; } + } + return true; + } } // namespace lowe \ No newline at end of file diff --git a/dunereco/LowEUtils/LowEUtils.h b/dunereco/LowEUtils/LowEUtils.h index 6b92fed6..4aeb7322 100644 --- a/dunereco/LowEUtils/LowEUtils.h +++ b/dunereco/LowEUtils/LowEUtils.h @@ -230,6 +230,18 @@ namespace lowe const float &ClusterCharge, const float &OpFlashPE); + bool CutPDSFlashMinPE( + const float &TPCDriftTime, + const float &MatchedDriftTime, + const float &ClusterCharge, + const float &OpFlashPE); + + bool CutPDSFlashMaxPE( + const float &TPCDriftTime, + const float &MatchedDriftTime, + const float &ClusterCharge, + const float &OpFlashPE); + // Declare member data here. private: // From fhicl configuration @@ -251,13 +263,18 @@ namespace lowe const double fAdjOpFlashY; // Maximum Y distance for adjacent flash matching [cm] const double fAdjOpFlashZ; // Maximum Z distance for adjacent flash matching [cm] const double fAdjOpFlashMinPECut; // Minimum photoelectrons for adjacent flash + const double fAdjOpFlashMaxPECut; // Maximum photoelectrons for adjacent flash const double fAdjOpFlashMaxPERatioCut; // Maximum photoelectrons ratio for adjacent flash const bool fAdjOpFlashMembraneProjection; // If true, project the TPC reco onto the membrane const bool fAdjOpFlashEndCapProjection; // If true, project the TPC reco onto the end cap const double fAdjOpFlashMinPEAttenuation; // Attenuation factor for minimum PE cut based on drift time [us] + const double fAdjOpFlashMaxPEAttenuation; // Attenuation factor for maximum PE cut based on drift time [us] const std::string fAdjOpFlashMinPEAttenuate; // Type of attenuation for minimum PE cut ("light_map", "asymptotic", "linear" or "flat") + const std::string fAdjOpFlashMaxPEAttenuate; // Type of attenuation for maximum PE cut ("light_map", "asymptotic", "linear" or "flat") const int fAdjOpFlashMinPEAttenuationStrength; // Strength of the asymptotic attenuation for minimum PE cut (in powers of 10) + const int fAdjOpFlashMaxPEAttenuationStrength; // Strength of the asymptotic attenuation for maximum PE cut (in powers of 10) const std::vector>> fAdjOpFlashMinPELightMap; // Light map file and histogram name for light map attenuation + const std::vector>> fAdjOpFlashMaxPELightMap; // Light map file and histogram name for light map attenuation std::unique_ptr producer; // Pointer to the ProducerUtils instance }; } // namespace lowe From 4aaa954096f85da696db12aeded599feaf67fe78 Mon Sep 17 00:00:00 2001 From: mantheys Date: Sun, 4 Jan 2026 03:37:41 -0600 Subject: [PATCH 2/9] UPDATE LIGHTMAP DATA AND MOVE TO SEPARATE FHICL --- dunereco/LowEUtils/fcl/LightMapData.fcl | 39 +++++++++++++++++++++++++ dunereco/LowEUtils/fcl/SolarEvent.fcl | 26 ++++++++--------- 2 files changed, 51 insertions(+), 14 deletions(-) create mode 100644 dunereco/LowEUtils/fcl/LightMapData.fcl diff --git a/dunereco/LowEUtils/fcl/LightMapData.fcl b/dunereco/LowEUtils/fcl/LightMapData.fcl new file mode 100644 index 00000000..90f4974a --- /dev/null +++ b/dunereco/LowEUtils/fcl/LightMapData.fcl @@ -0,0 +1,39 @@ +BEGIN_PROLOG + +light_map_data_dune10kt_1x2x6_AdjOpFlashMinPELightMap: [ + ["Amplitude", [-1.2723e-07, 7.1570e-04, 7.1178e-01]], + ["Attenuation", [1.4256e-07, -4.7368e-04, 8.9290e-01]], + ["Correction", [-4.0734e-08, 3.7639e-04, -1.0007e+00]] + ] + +light_map_data_dune10kt_1x2x6_AdjOpFlashMaxPELightMap: [ + ["Amplitude", [-1.1814e-07, 5.7685e-04, 2.8689e+00]], + ["Attenuation", [3.4830e-08, -1.3746e-04, 3.9619e-01]], + ["Correction", [6.1369e-08, -2.1512e-04, 1.2846e-01]] + ] + +light_map_data_dune10kt_1x2x6_AdjOpFlashPELightMap: [ + ["Amplitude", [-1.2784e-07, 6.5108e-04, 2.3573e+00]], + ["Attenuation", [2.5445e-08, -1.2782e-04, 4.4181e-01]], + ["Correction", [-1.1812e-08, -1.6604e-05, -6.5942e-02]] + ] + +light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMinPELightMap: [ + ["Amplitude", [2.5276e-08, 2.9590e-04, -7.2231e-02]], + ["Attenuation", [-1.1732e-06, 8.1172e-03, -1.5276e+01]], + ["Correction", [7.6294e-07, -2.8500e-03, 1.9332e+00]] + ] + +light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMaxPELightMap: [ + ["Amplitude", [-8.8197e-08, 5.2212e-04, 2.1964e+00]], + ["Attenuation", [4.8050e-08, -9.5364e-05, -5.7437e-01]], + ["Correction", [3.7871e-07, -1.2702e-03, 7.0594e-01]] + ] + +light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashPELightMap: [ + ["Amplitude", [-1.2053e-07, 7.1549e-04, 7.2371e-01]], + ["Attenuation", [-7.9112e-08, 4.6725e-04, -1.7473e+00]], + ["Correction", [2.5340e-07, -1.0071e-03, 5.7146e-01]] + ] + +END_PROLOG \ No newline at end of file diff --git a/dunereco/LowEUtils/fcl/SolarEvent.fcl b/dunereco/LowEUtils/fcl/SolarEvent.fcl index e2509efb..1b45e950 100644 --- a/dunereco/LowEUtils/fcl/SolarEvent.fcl +++ b/dunereco/LowEUtils/fcl/SolarEvent.fcl @@ -1,3 +1,5 @@ +#include LightMapData.fcl + BEGIN_PROLOG solar_event_dune10kt_1x2x6: @@ -21,15 +23,13 @@ solar_event_dune10kt_1x2x6: AdjOpFlashMinNHitCut: 3 # Cut on the minimum number of OpHits in the OpFlash. AdjOpFlashMembraneProjection: true # If true, the OpFlash matching is projected on the membrane planes for VD. AdjOpFlashEndCapProjection: false # If true, the OpFlash matching is projected on the end cap planes for VD. - AdjOpFlashMinPEAttenuation: 1.0 # Attenuation factor over full drift time. - AdjOpFlashMinPEAttenuate: "flat" # Type of attenuation: "flat", "linear", "asymptotic" or "light_map" - AdjOpFlashMinPEAttenuationStrength: 10 # Strength of asymptotic attenuation (in powers of 10). - AdjOpFlashMinPELightMap: [ - ["Amplitude", [-1.1011e-7, 5.0670e-4, 1.6480e0]], - ["Attenuation", [1.7061e-7, -7.2303e-4, 1.0623e0]], - ["Correction", [3.1135e-7, -1.3050e-3, 7.6113e-1]] - ] # Quadratic fit parameters for attenuation correction based in logarithm light map. - + AdjOpFlashMinPEAttenuation: 1.0 # Attenuation factor over full drift time. + AdjOpFlashMinPEAttenuate: "flat" # Type of attenuation: "flat", "linear", "asymptotic" or "light_map" + AdjOpFlashMinPEAttenuationStrength: 10 # Strength of asymptotic attenuation (in powers of 10). + AdjOpFlashMinPELightMap: @local::light_map_data_dune10kt_1x2x6_AdjOpFlashMinPELightMap + AdjOpFlashMaxPELightMap: @local::light_map_data_dune10kt_1x2x6_AdjOpFlashMaxPELightMap + AdjOpFlashPELightMap: @local::light_map_data_dune10kt_1x2x6_AdjOpFlashPELightMap + FlashMatchBy: "light_map" # Match flashes by residual. Alternative is to match by MaxFlashPE. Debug: true } @@ -37,11 +37,9 @@ solar_event_dunevd10kt_1x8x14_3view_30deg: @local::solar_event_dune10kt_1x2x6 solar_event_dunevd10kt_1x8x14_3view_30deg.AdjOpFlashX: 140 solar_event_dunevd10kt_1x8x14_3view_30deg.AdjOpFlashY: 140 solar_event_dunevd10kt_1x8x14_3view_30deg.AdjOpFlashZ: 140 -solar_event_dunevd10kt_1x8x14_3view_30deg.AdjOpFlashMinPELightMap: [ - ["Amplitude", [1.7578e-7, -4.9589e-4, 1.2519e0]], - ["Attenuation", [4.9682e-7, -2.5242e-3, 1.2541e0]], - ["Correction", [1.0516e-7, -1.1637e-3, 1.5118e0]] -] +solar_event_dunevd10kt_1x8x14_3view_30deg.AdjOpFlashMinPELightMap: @local::light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMinPELightMap +solar_event_dunevd10kt_1x8x14_3view_30deg.AdjOpFlashMaxPELightMap: @local::light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMaxPELightMap +solar_event_dunevd10kt_1x8x14_3view_30deg.AdjOpFlashPELightMap: @local::light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashPELightMap solar_event_dunevd10kt_1x8x6_3view_30deg: @local::solar_event_dunevd10kt_1x8x14_3view_30deg From 4be0824e403d324924619841a6a958802bf0d34b Mon Sep 17 00:00:00 2001 From: mantheys Date: Sun, 4 Jan 2026 03:38:25 -0600 Subject: [PATCH 3/9] REFORMAT AND EXPAND BY ADDING LIGHT MAP MATCHING FUNCTIONS --- dunereco/LowEUtils/LowEUtils.cc | 412 +++++++++++++++++++++++++------- dunereco/LowEUtils/LowEUtils.h | 52 ++-- 2 files changed, 360 insertions(+), 104 deletions(-) diff --git a/dunereco/LowEUtils/LowEUtils.cc b/dunereco/LowEUtils/LowEUtils.cc index 44d561a2..10b6c152 100644 --- a/dunereco/LowEUtils/LowEUtils.cc +++ b/dunereco/LowEUtils/LowEUtils.cc @@ -22,26 +22,35 @@ namespace lowe fAdjOpFlashX(p.get("AdjOpFlashX", 100.0)), // X coordinate for flash projection [cm] fAdjOpFlashY(p.get("AdjOpFlashY", 100.0)), // Y coordinate for flash projection [cm] fAdjOpFlashZ(p.get("AdjOpFlashZ", 100.0)), // Z coordinate for flash projection [cm] - fAdjOpFlashMinPECut(p.get("AdjOpFlashMinPECut", 0.0)), // Minimum PE for flash selection - fAdjOpFlashMaxPECut(p.get("AdjOpFlashMaxPECut", 1e10)), // Maximum PE for flash selection fAdjOpFlashMaxPERatioCut(p.get("AdjOpFlashMaxPERatioCut", 1.0)), // Maximum PE ratio for flash selection + fAdjOpFlashMinPECut(p.get("AdjOpFlashMinPECut", 0.0)), // Minimum PE for flash selection + fAdjOpFlashMaxPECut(p.get("AdjOpFlashMaxPECut", 1e9)), // Maximum PE for flash selection fAdjOpFlashMembraneProjection(p.get("AdjOpFlashMembraneProjection", false)), // Whether to project flashes onto the membrane fAdjOpFlashEndCapProjection(p.get("AdjOpFlashEndCapProjection", false)), // Whether to project flashes onto the end cap - fAdjOpFlashMinPEAttenuation(p.get("AdjOpFlashMinPEAttenuation", 1)), // PE cut attenuation over drift time [us] to compensate for light attenuation - fAdjOpFlashMaxPEAttenuation(p.get("AdjOpFlashMaxPEAttenuation", 1)), // PE cut attenuation over drift time [us] to compensate for light attenuation + fAdjOpFlashMinPEAttenuation(p.get("AdjOpFlashMinPEAttenuation", 0.0)), // PE cut attenuation over drift time [us] to compensate for light attenuation + fAdjOpFlashMaxPEAttenuation(p.get("AdjOpFlashMaxPEAttenuation", 0.0)), // PE cut attenuation over drift time [us] to compensate for light attenuation fAdjOpFlashMinPEAttenuate(p.get("AdjOpFlashMinPEAttenuate", "flat")), // Type of attenuation for minimum PE cut ("light_map", "asymptotic", "linear" or "flat") fAdjOpFlashMaxPEAttenuate(p.get("AdjOpFlashMaxPEAttenuate", "flat")), // Type of attenuation for maximum PE cut ("light_map", "asymptotic", "linear" or "flat") - fAdjOpFlashMinPEAttenuationStrength(p.get("AdjOpFlashMinPEAttenuationStrength", 10)), // Strength of the asymptotic attenuation for minimum PE cut (in powers of 10) - fAdjOpFlashMaxPEAttenuationStrength(p.get("AdjOpFlashMaxPEAttenuationStrength", 10)), // Strength of the asymptotic attenuation for maximum PE cut (in powers of 10) + fAdjOpFlashMinPEAttenuationStrength(p.get("AdjOpFlashMinPEAttenuationStrength", 0)), // Strength of the asymptotic attenuation for minimum PE cut (in powers of 10) + fAdjOpFlashMaxPEAttenuationStrength(p.get("AdjOpFlashMaxPEAttenuationStrength", 0)), // Strength of the asymptotic attenuation for maximum PE cut (in powers of 10) fAdjOpFlashMinPELightMap(p.get>>>("AdjOpFlashMinPELightMap", {})), // Light map file and histogram name for light map attenuation fAdjOpFlashMaxPELightMap(p.get>>>("AdjOpFlashMaxPELightMap", {})), // Light map file and histogram name for light map attenuation + fAdjOpFlashPELightMap(p.get>>>("AdjOpFlashPELightMap", {})), // Light map file and histogram name for PE attenuation + fFlashMatchBy(p.get("FlashMatchBy", "maximum")), // Method to match flashes ("maximum" or "light_map") producer(new ProducerUtils(p)) { // Initialize the LowEUtils instance producer->PrintInColor("LowEUtils initialized with parameters from FHiCL configuration.", ProducerUtils::GetColor("green"), "Debug"); } - void LowEUtils::MakeClusterVector(std::vector &ClusterVec, std::vector>> &Clusters, art::Event const &evt) + //...................................................... + void LowEUtils::MakeClusterVector( + std::vector &ClusterVec, + std::vector>> &Clusters, + art::Event const &evt) + /* + Create a vector of RawPerPlaneCluster from a vector of clusters of hits + */ { mf::LogDebug("LowEUtils") << "Charge variable set to " << fClusterChargeVariable; int ID = 0; @@ -176,9 +185,15 @@ namespace lowe return; } - void LowEUtils::CalcAdjHits(std::vector MyVec, std::vector> &Clusters, TH1I *MyHist, TH1F *ADCIntHist, art::Event const &evt, bool debug) + //...................................................... + void LowEUtils::CalcAdjHits( + std::vector MyVec, + std::vector> &Clusters, + TH1I *MyHist, TH1F *ADCIntHist, + art::Event const &evt, + bool debug) /* - Find adjacent hits in time and space: + Find adjacent hits in time and space. Obsolete (kept for backwards compatibility), use the other overload! - MyVec is the vector of hits to be clustered - Clusters is the vector of clusters - MyHist is the histogram to be filled with the number of hits in each cluster @@ -328,7 +343,10 @@ namespace lowe return; } - void LowEUtils::CalcAdjHits(std::vector> MyVec, std::vector>> &Clusters, std::vector> &ClusterIdx, art::Event const &evt) + void LowEUtils::CalcAdjHits(std::vector> MyVec, + std::vector>> &Clusters, + std::vector> &ClusterIdx, + art::Event const &evt) /* Find adjacent hits in time and space: - MyVec is the vector of hits to be clustered @@ -430,6 +448,18 @@ namespace lowe std::vector &Dir, detinfo::DetectorClocksData ClockData) /* + Fill hit-level vectors for a given cluster of hits + - Cluster: vector of hits in the cluster + - TPC: output vector of TPC numbers + - Channel: output vector of channel numbers + - Charge: output vector of hit charges + - Time: output vector of hit times + - SigmaTime: output cluster time standard deviation + - DeltaTime: output cluster time span + - Y: output vector of hit Y positions + - Z: output vector of hit Z positions + - Dir: output vector of hit directions (dZ/dY) + - ClockData: detector clock data for time conversion */ { // --- Clear the vectors @@ -481,16 +511,25 @@ namespace lowe //...................................................... std::vector LowEUtils::ComputeRecoY( - int Event, - std::vector &HIndTPC, - std::vector &Z, - std::vector &Time, - std::vector &IndZ, - std::vector &IndY, - std::vector &IndT, - std::vector &IndDir, - bool debug) + int Event, + std::vector &HIndTPC, + std::vector &Z, + std::vector &Time, + std::vector &IndZ, + std::vector &IndY, + std::vector &IndT, + std::vector &IndDir, + bool debug) /* + Compute the reconstructed Y positions for a set of hits given the induction plane cluster reference points + - Event: event number (for debugging purposes) + - HIndTPC: vector of TPC numbers for the hits + - Z: vector of Z positions for the hits + - Time: vector of times for the hits + - IndZ: vector of Z reference positions for the induction plane cluster + - IndY: vector of Y reference positions for the induction plane cluster + - IndT: vector of time reference positions for the induction plane cluster + - IndDir: vector of direction cosines (dZ/dY) for the induction plane cluster */ { // Create the interpolator @@ -518,6 +557,7 @@ namespace lowe return RecoY; } + //...................................................... void LowEUtils::FillClusterVariables( std::vector>> Clusters, std::vector> &ClNHits, @@ -525,6 +565,13 @@ namespace lowe std::vector> &ClCharge, art::Event const &evt, bool debug) + /* + Fill basic cluster variables: number of hits, time, charge + - Clusters: input vector of clusters of hits + - ClNHits: output vector of number of hits per cluster + - ClT: output vector of cluster times + - ClCharge: output vector of cluster charges + */ { detinfo::DetectorClocksData const &clockData = art::ServiceHandle()->DataFor(evt); auto TickPeriod = clockData.TPCClock().TickPeriod(); @@ -547,6 +594,7 @@ namespace lowe } // End of loop over clusters } // End of loop over planes } + void LowEUtils::FillClusterVariables( std::set SignalTrackIDs, std::vector>> Clusters, @@ -563,6 +611,22 @@ namespace lowe std::vector> &ClCompleteness, detinfo::DetectorClocksData const &clockData, bool debug) + /* + Fill detailed cluster variables: main track ID, number of hits, TPC, channel, time, Y, Z, direction, charge, purity, completeness + - SignalTrackIDs: set of signal track IDs for purity/completeness calculation + - Clusters: input vector of clusters of hits + - ClMainID: output vector of main track IDs per cluster + - ClNHits: output vector of number of hits per cluster + - ClTPC: output vector of TPC numbers per cluster + - ClChannel: output vector of channel numbers per cluster + - ClT: output vector of cluster times + - ClY: output vector of cluster Y positions + - ClZ: output vector of cluster Z positions + - ClDir: output vector of cluster directions (dZ/dY) + - ClCharge: output vector of cluster charges + - ClPurity: output vector of cluster purities + - ClCompleteness: output vector of cluster completenesses + */ { art::ServiceHandle bt_serv; geo::WireReadoutGeom const &wireReadout = art::ServiceHandle()->Get(); @@ -636,25 +700,45 @@ namespace lowe } // End of loop over planes } + //...................................................... void LowEUtils::MatchClusters( - std::set SignalTrackIDs, - std::vector> &MatchedClustersIdx, - std::vector>> &MatchedClusters, - std::vector> ClustersIdx, - std::vector>> Clusters, - std::vector> &ClMainID, - std::vector> &ClNHits, - std::vector> &ClTPC, - std::vector> &ClChannel, - std::vector> &ClT, - std::vector> &ClY, - std::vector> &ClZ, - std::vector> &ClDir, - std::vector> &ClCharge, - std::vector> &ClPurity, - std::vector> &ClCompleteness, - detinfo::DetectorClocksData const &clockData, - bool debug) + std::set SignalTrackIDs, + std::vector> &MatchedClustersIdx, + std::vector>> &MatchedClusters, + std::vector> ClustersIdx, + std::vector>> Clusters, + std::vector> &ClMainID, + std::vector> &ClNHits, + std::vector> &ClTPC, + std::vector> &ClChannel, + std::vector> &ClT, + std::vector> &ClY, + std::vector> &ClZ, + std::vector> &ClDir, + std::vector> &ClCharge, + std::vector> &ClPurity, + std::vector> &ClCompleteness, + detinfo::DetectorClocksData const &clockData, + bool debug) + /* + Match clusters across the three planes based on time and other criteria + - SignalTrackIDs: set of signal track IDs for purity/completeness calculation + - MatchedClustersIdx: output vector of matched cluster indices + - MatchedClusters: output vector of matched clusters of hits + - ClustersIdx: input vector of cluster indices + - Clusters: input vector of clusters of hits + - ClMainID: input vector of main track IDs per cluster + - ClNHits: input vector of number of hits per cluster + - ClTPC: input vector of TPC numbers per cluster + - ClChannel: input vector of channel numbers per cluster + - ClT: input vector of cluster times + - ClY: input vector of cluster Y positions + - ClZ: input vector of cluster Z positions + - ClDir: input vector of cluster directions (dZ/dY) + - ClCharge: input vector of cluster charges + - ClPurity: input vector of cluster purities + - ClCompleteness: input vector of cluster completenesses + */ { LowEUtils::FillClusterVariables(SignalTrackIDs, Clusters, ClMainID, ClNHits, ClTPC, ClChannel, ClT, ClY, ClZ, ClDir, ClCharge, ClPurity, ClCompleteness, clockData, debug); std::vector> MatchedClMainID = {{}, {}, {}}, MatchedClNHits = {{}, {}, {}}, MatchedClTPC = {{}, {}, {}}, MatchedClChannel = {{}, {}, {}}; @@ -891,13 +975,21 @@ namespace lowe } void LowEUtils::MatchClusters( - std::vector>> &MatchedClusters, - std::vector>> Clusters, - std::vector> &ClNHits, - std::vector> &ClT, - std::vector> &ClCharge, - art::Event const &evt, - bool debug) + std::vector>> &MatchedClusters, + std::vector>> Clusters, + std::vector> &ClNHits, + std::vector> &ClT, + std::vector> &ClCharge, + art::Event const &evt, + bool debug) + /* + Match clusters across the three planes based on time and other criteria + - MatchedClusters: output vector of matched clusters of hits + - Clusters: input vector of clusters of hits + - ClNHits: output vector of number of hits per cluster + - ClT: output vector of cluster times + - ClCharge: output vector of cluster charges + */ { LowEUtils::FillClusterVariables(Clusters, ClNHits, ClT, ClCharge, evt, debug); std::vector> MatchedClNHits = {{}, {}, {}}; @@ -1023,7 +1115,12 @@ namespace lowe } //...................................................... - double LowEUtils::STD(const std::vector& Vec) + double LowEUtils::STD( + const std::vector& Vec) + /* + Calculate the standard deviation of a vector of doubles + - Vec: input vector of doubles + */ { if (Vec.size() == 0) { @@ -1041,7 +1138,22 @@ namespace lowe } //...................................................... - void LowEUtils::FindPrimaryClusters(const std::vector> &SolarClusterVector, std::vector &EventCandidateFound, std::vector>> &EventCandidateVector, std::vector> &EventCandidateIdx, const detinfo::DetectorClocksData &clockData, const art::Event &evt) + void LowEUtils::FindPrimaryClusters( + const std::vector> &SolarClusterVector, + std::vector &EventCandidateFound, + std::vector>> &EventCandidateVector, + std::vector> &EventCandidateIdx, + const detinfo::DetectorClocksData &clockData, + const art::Event &evt) + /* + Find primary clusters and group adjacent clusters into event candidates + - SolarClusterVector: input vector of LowECluster pointers + - EventCandidateFound: output vector of booleans indicating if an event candidate was found for each cluster + - EventCandidateVector: output vector of vectors of LowECluster pointers for each event candidate + - EventCandidateIdx: output vector of vectors of indices for each event candidate + - clockData: detector clocks data + - evt: art event + */ { // This is the low energy primary cluster algorithm. It groups all input clusters into event candidates by finding the primary clusters (charge > adjacent clusters up to distance fAdjClusterRad). // The algorithm outputs vectors of clusters where the first entry is the primary cluster and the rest are the corresponding adjacent clusters. @@ -1261,12 +1373,22 @@ namespace lowe return; } + int LowEUtils::MatchPDSFlash( const std::vector> &SolarClusterVector, const std::vector> &PDSFlashes, const detinfo::DetectorClocksData &clockData, const art::Event &evt, bool debug) + /* + Match the main cluster in SolarClusterVector to the best matching PDS flash in PDSFlashes. + Returns the index of the matched flash, or -1 if no match is found. + - SolarClusterVector: vector of LowECluster pointers to match + - PDSFlashes: vector of OpFlash pointers to match against + - clockData: detector clocks data for time conversions + - evt: art event for geometry and detector properties + - debug: if true, print debug information + */ { std::string sFlashMatching = "LowEUtils::MatchPDSFlash " + ProducerUtils::str(SolarClusterVector.size()) + " clusters and " + ProducerUtils::str(PDSFlashes.size()) + " flashes found in the event\n"; // get drift properties @@ -1319,24 +1441,22 @@ namespace lowe // Loop through the flashes to find the best match for (int i = 0; i < int(PDSFlashes.size()); i++) { - // Reset cluster X for each flash - const auto &flash = PDSFlashes[i]; - double flashPE = flash->TotalPE(); - // double flashR = -1e6; - // Check if the flash time is within the acceptable range - double flashTime = flash->Time(); - double dT = clusterTime - flashTime; // Time difference between the cluster and the flash - // If dT is bigger that drift time, skip this flash - if (dT > driftTime || dT < 0) { - continue; // Skip this flash if it's too far in time or in the future - } - if (int(flash->PEs().size()) < fAdjOpFlashMinNHitCut || *std::max_element(flash->PEs().begin(), flash->PEs().end()) / flash->TotalPE() > fAdjOpFlashMaxPERatioCut) { - continue; // Skip flashes with insufficient hits or too much concentration in one XA - } - if ( SelectPDSFlashPE(driftTime, dT, clusterCharge, flashPE) == false ) { - continue; // Skip flashes that don't pass the PE selection - } - + // Reset cluster X for each flash + const auto &flash = PDSFlashes[i]; + double flashPE = flash->TotalPE(); + // double flashR = -1e6; + // Check if the flash time is within the acceptable range + double flashTime = flash->Time(); + double dT = clusterTime - flashTime; // Time difference between the cluster and the flash + // If dT is bigger that drift time, skip this flash + if (dT > driftTime || dT < 0) { + continue; // Skip this flash if it's too far in time or in the future + } + if (int(flash->PEs().size()) < fAdjOpFlashMinNHitCut || *std::max_element(flash->PEs().begin(), flash->PEs().end()) / flash->TotalPE() > fAdjOpFlashMaxPERatioCut) { + continue; // Skip flashes with insufficient hits or too much concentration in one XA + } + if ( SelectPDSFlashPE(driftTime, dT, clusterCharge, flashPE) ) + { // Calculate the distance in space double flashY = flash->YCenter(); double flashZ = flash->ZCenter(); @@ -1396,7 +1516,7 @@ namespace lowe clusterX = driftLength / 2 - clusterX; // Convert to the VD geometry coordinate system } else { - continue; // Unknown geometry, skip this matching + continue; // Unknown geometry, skip this matching } if (flashPE > matchedFlashPE || matchedFlashIndex == -1) { @@ -1405,6 +1525,7 @@ namespace lowe matchedFlashIndex = i; matchedRecoX = clusterX; } + } } if (debug) { if (matchedFlashIndex != -1) { @@ -1426,30 +1547,46 @@ namespace lowe const float &MatchedDriftTime, const float &ClusterCharge, const float &OpFlashPE) + /* + Select PDS flash based on PE cuts + - TPCDriftTime: total drift time in the TPC + - MatchedDriftTime: drift time corresponding to the matched flash + - ClusterCharge: charge of the cluster + - OpFlashPE: total PE of the flash + */ { return CutPDSFlashMinPE(TPCDriftTime, MatchedDriftTime, ClusterCharge, OpFlashPE) && CutPDSFlashMaxPE(TPCDriftTime, MatchedDriftTime, ClusterCharge, OpFlashPE); } // SelectPDSFlashPE + bool LowEUtils::CutPDSFlashMinPE( const float &TPCDriftTime, const float &MatchedDriftTime, const float &ClusterCharge, const float &OpFlashPE) + /* + Apply minimum PE cut to PDS flash + - TPCDriftTime: total drift time in the TPC + - MatchedDriftTime: drift time corresponding to the matched flash + - ClusterCharge: charge of the cluster + - OpFlashPE: total PE of the flash + */ { if (fAdjOpFlashMinPEAttenuate == "light_map") { - double a, b, c; - for (int i = 0; i < int(fAdjOpFlashMinPELightMap.size()); i++) { - if (fAdjOpFlashMinPELightMap[i].first == "Amplitude") { - a = fAdjOpFlashMinPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMinPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMinPELightMap[i].second[2]; - } - else if (fAdjOpFlashMinPELightMap[i].first == "Attenuation") { - b = fAdjOpFlashMinPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMinPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMinPELightMap[i].second[2]; - } - else if (fAdjOpFlashMinPELightMap[i].first == "Correction") { - c = fAdjOpFlashMinPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMinPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMinPELightMap[i].second[2]; - } - } + float a, b, c; + GetLightMapParameters("min", ClusterCharge, a, b, c); + // for (int i = 0; i < int(fAdjOpFlashMinPELightMap.size()); i++) { + // if (fAdjOpFlashMinPELightMap[i].first == "Amplitude") { + // a = fAdjOpFlashMinPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMinPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMinPELightMap[i].second[2]; + // } + // else if (fAdjOpFlashMinPELightMap[i].first == "Attenuation") { + // b = fAdjOpFlashMinPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMinPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMinPELightMap[i].second[2]; + // } + // else if (fAdjOpFlashMinPELightMap[i].first == "Correction") { + // c = fAdjOpFlashMinPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMinPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMinPELightMap[i].second[2]; + // } + // } double MinPE = pow(10, a - a * b * MatchedDriftTime / TPCDriftTime + c * pow(MatchedDriftTime / TPCDriftTime, 2)); if (OpFlashPE < MinPE) { return false; } else { return true; } @@ -1474,25 +1611,34 @@ namespace lowe return true; } + bool LowEUtils::CutPDSFlashMaxPE( const float &TPCDriftTime, const float &MatchedDriftTime, const float &ClusterCharge, const float &OpFlashPE) + /* + Apply maximum PE cut to PDS flash + - TPCDriftTime: total drift time in the TPC + - MatchedDriftTime: drift time corresponding to the matched flash + - ClusterCharge: charge of the cluster + - OpFlashPE: total PE of the flash + */ { if (fAdjOpFlashMaxPEAttenuate == "light_map") { - double a, b, c; - for (int i = 0; i < int(fAdjOpFlashMaxPELightMap.size()); i++) { - if (fAdjOpFlashMaxPELightMap[i].first == "Amplitude") { - a = fAdjOpFlashMaxPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMaxPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMaxPELightMap[i].second[2]; - } - else if (fAdjOpFlashMaxPELightMap[i].first == "Attenuation") { - b = fAdjOpFlashMaxPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMaxPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMaxPELightMap[i].second[2]; - } - else if (fAdjOpFlashMaxPELightMap[i].first == "Correction") { - c = fAdjOpFlashMaxPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMaxPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMaxPELightMap[i].second[2]; - } - } + float a, b, c; + GetLightMapParameters("max", ClusterCharge, a, b, c); + // for (int i = 0; i < int(fAdjOpFlashMaxPELightMap.size()); i++) { + // if (fAdjOpFlashMaxPELightMap[i].first == "Amplitude") { + // a = fAdjOpFlashMaxPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMaxPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMaxPELightMap[i].second[2]; + // } + // else if (fAdjOpFlashMaxPELightMap[i].first == "Attenuation") { + // b = fAdjOpFlashMaxPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMaxPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMaxPELightMap[i].second[2]; + // } + // else if (fAdjOpFlashMaxPELightMap[i].first == "Correction") { + // c = fAdjOpFlashMaxPELightMap[i].second[0] * pow(ClusterCharge, 2) + fAdjOpFlashMaxPELightMap[i].second[1] * ClusterCharge + fAdjOpFlashMaxPELightMap[i].second[2]; + // } + // } double MaxPE = pow(10, a - a * b * MatchedDriftTime / TPCDriftTime + c * pow(MatchedDriftTime / TPCDriftTime, 2)); if (OpFlashPE > MaxPE) { return false; } else { return true; } @@ -1516,4 +1662,92 @@ namespace lowe } return true; } + + void LowEUtils::GetLightMapParameters( + const std::string &LightMapType, + const float &ClusterCharge, + float &a, + float &b, + float &c) + /* + Store light map parameters for flash selection + - LightMapType: type of light map ("min", "max", "med") + - a: Amplitude parameters + - b: Attenuation parameters + - c: Correction parameters + */ + { + std::vector>> SelectedLightMap; + if (LightMapType == "min") { + SelectedLightMap = fAdjOpFlashMinPELightMap; + } + else if (LightMapType == "max") { + SelectedLightMap = fAdjOpFlashMaxPELightMap; + } + else { + if (LightMapType != "med") { + producer->PrintInColor("Warning: Unknown light map type " + LightMapType + ", using med", ProducerUtils::GetColor("red"), "Warning"); + } + SelectedLightMap = fAdjOpFlashPELightMap; + } + for (int i = 0; i < int(SelectedLightMap.size()); i++) { + if (SelectedLightMap[i].first == "Amplitude") { + a = SelectedLightMap[i].second[0] * pow(ClusterCharge, 2) + SelectedLightMap[i].second[1] * ClusterCharge + SelectedLightMap[i].second[2]; + } + else if (SelectedLightMap[i].first == "Attenuation") { + b = SelectedLightMap[i].second[0] * pow(ClusterCharge, 2) + SelectedLightMap[i].second[1] * ClusterCharge + SelectedLightMap[i].second[2]; + } + else if (SelectedLightMap[i].first == "Correction") { + c = SelectedLightMap[i].second[0] * pow(ClusterCharge, 2) + SelectedLightMap[i].second[1] * ClusterCharge + SelectedLightMap[i].second[2]; + } + } + return; + } + + bool LowEUtils::SelectPDSFlash( + const float &TPCDriftTime, + const float &ClusterTime, + const float &ClusterCharge, + const float &RefOpFlashTime, + const float &RefOpFlashPE, + const float &OpFlashTime, + const float &OpFlashPE) + /* + Select PDS flash based on time and PE. There are 2 options: "maximum" or "light_map" + - TPCDriftTime: total drift time in the TPC + - ClusterTime: time of the cluster + - ClusterCharge: charge of the cluster + - RefOpFlashTime: time of the reference flash + - RefOpFlashTime: time of the reference flash + - OpFlashPE: total PE of the flash + - OpFlashPE: total PE of the flash + */ + { + if (fFlashMatchBy == "light_map") { + float a, b, c; + GetLightMapParameters("med", ClusterCharge, a, b, c); + double RefPE = pow(10, a - a * b * std::abs(ClusterTime - RefOpFlashTime) / TPCDriftTime + c * pow(std::abs(ClusterTime - RefOpFlashTime) / TPCDriftTime, 2)); + double SelectedPE = pow(10, a - a * b * std::abs(ClusterTime - OpFlashTime) / TPCDriftTime + c * pow(std::abs(ClusterTime - OpFlashTime) / TPCDriftTime, 2)); + // If the flash PE is closer to the selected PE, than the reference flash PE, select it + if (std::abs(OpFlashPE - SelectedPE) < std::abs(RefOpFlashPE - RefPE)) { + return true; + } + else { + return false; + } + } + else { + // If the flashMatchBy is not "maximum" print a warning + if (fFlashMatchBy != "maximum") { + ProducerUtils::PrintInColor("Warning: Unknown fFlashMatchBy option " + fFlashMatchBy + ", using maximum PE selection", ProducerUtils::GetColor("red"), "Warning"); + } + // If the flash PE is bigger than the reference flash PE, select it + if (OpFlashPE > RefOpFlashPE) { + return true; + } + else { + return false; + } + } + } } // namespace lowe \ No newline at end of file diff --git a/dunereco/LowEUtils/LowEUtils.h b/dunereco/LowEUtils/LowEUtils.h index 4aeb7322..3a8e504d 100644 --- a/dunereco/LowEUtils/LowEUtils.h +++ b/dunereco/LowEUtils/LowEUtils.h @@ -113,11 +113,24 @@ namespace lowe std::vector ClusterIdxVec; }; - void CalcAdjHits(std::vector MyVec, std::vector> &Clusters, TH1I *MyHist, TH1F *ADCIntHist, art::Event const &evt, bool debug); + void CalcAdjHits( + std::vector MyVec, + std::vector> &Clusters, + TH1I *MyHist, + TH1F *ADCIntHist, + art::Event const &evt, + bool debug); - void CalcAdjHits(std::vector> MyVec, std::vector>> &Clusters, std::vector> &ClusterIdx, art::Event const &evt); + void CalcAdjHits( + std::vector> MyVec, + std::vector>> &Clusters, + std::vector> &ClusterIdx, + art::Event const &evt); - void MakeClusterVector(std::vector &ClusterVec, std::vector>> &Clusters, art::Event const &evt); + void MakeClusterVector( + std::vector &ClusterVec, + std::vector>> &Clusters, + art::Event const &evt); void FillClusterVariables( std::vector>> Clusters, @@ -197,7 +210,8 @@ namespace lowe std::vector &IndDir, bool debug = false); - double STD(const std::vector& Vec); + double STD( + const std::vector& Vec); void FindPrimaryClusters( const std::vector> &SolarClusterVector, @@ -207,16 +221,6 @@ namespace lowe const detinfo::DetectorClocksData &clockData, const art::Event &evt); - /** - * @brief MatchPDSFlash matches a vector of LowEClusters with a vector of PDS flashes. - * @param SolarClusterVector Vector of LowEClusters to match. - * @param PDSFlashes Vector of PDS flashes to match against. - * @param clockData DetectorClocksData for time calculations. - * @param debug If true, enables debug output. - * @return Returns the index of the matched flash, or -1 if no match is found. - * @note This function uses the time and charge of the clusters to find the best matching flash. - * It assumes that the clusters are already sorted by time. - */ int MatchPDSFlash( const std::vector> &SolarClusterVector, const std::vector> &PDSFlashes, @@ -242,6 +246,22 @@ namespace lowe const float &ClusterCharge, const float &OpFlashPE); + void GetLightMapParameters( + const std::string &LightMapType, + const float &ClusterCharge, + float &a, + float &b, + float &c); + + bool SelectPDSFlash( + const float &TPCDriftTime, + const float &ClusterTime, + const float &ClusterCharge, + const float &RefOpFlashTime, + const float &RefOpFlashPe, + const float &OpFlashTime, + const float &OpFlashPE); + // Declare member data here. private: // From fhicl configuration @@ -262,9 +282,9 @@ namespace lowe const double fAdjOpFlashX; // Maximum X distance for adjacent flash matching [cm] const double fAdjOpFlashY; // Maximum Y distance for adjacent flash matching [cm] const double fAdjOpFlashZ; // Maximum Z distance for adjacent flash matching [cm] + const double fAdjOpFlashMaxPERatioCut; // Maximum photoelectrons ratio for adjacent flash const double fAdjOpFlashMinPECut; // Minimum photoelectrons for adjacent flash const double fAdjOpFlashMaxPECut; // Maximum photoelectrons for adjacent flash - const double fAdjOpFlashMaxPERatioCut; // Maximum photoelectrons ratio for adjacent flash const bool fAdjOpFlashMembraneProjection; // If true, project the TPC reco onto the membrane const bool fAdjOpFlashEndCapProjection; // If true, project the TPC reco onto the end cap const double fAdjOpFlashMinPEAttenuation; // Attenuation factor for minimum PE cut based on drift time [us] @@ -275,6 +295,8 @@ namespace lowe const int fAdjOpFlashMaxPEAttenuationStrength; // Strength of the asymptotic attenuation for maximum PE cut (in powers of 10) const std::vector>> fAdjOpFlashMinPELightMap; // Light map file and histogram name for light map attenuation const std::vector>> fAdjOpFlashMaxPELightMap; // Light map file and histogram name for light map attenuation + const std::vector>> fAdjOpFlashPELightMap; // Light map file and histogram name for PE attenuation + const std::string fFlashMatchBy; // Method to match flashes ("maximum" or "light_map") std::unique_ptr producer; // Pointer to the ProducerUtils instance }; } // namespace lowe From ce507e39f99c7a3e7160075bf2a40424802e0f42 Mon Sep 17 00:00:00 2001 From: mantheys Date: Fri, 9 Jan 2026 14:01:19 -0600 Subject: [PATCH 4/9] UPDATE MAPS --- dunereco/LowEUtils/fcl/LightMapData.fcl | 36 ++++++++++++------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/dunereco/LowEUtils/fcl/LightMapData.fcl b/dunereco/LowEUtils/fcl/LightMapData.fcl index 90f4974a..793d6290 100644 --- a/dunereco/LowEUtils/fcl/LightMapData.fcl +++ b/dunereco/LowEUtils/fcl/LightMapData.fcl @@ -1,39 +1,39 @@ BEGIN_PROLOG light_map_data_dune10kt_1x2x6_AdjOpFlashMinPELightMap: [ - ["Amplitude", [-1.2723e-07, 7.1570e-04, 7.1178e-01]], - ["Attenuation", [1.4256e-07, -4.7368e-04, 8.9290e-01]], - ["Correction", [-4.0734e-08, 3.7639e-04, -1.0007e+00]] + ["Amplitude", [-1.2043e-07, 6.8624e-04, 7.3629e-01]], + ["Attenuation", [1.3162e-07, -6.0140e-04, 1.3823e+00]], + ["Correction", [-3.7917e-08, 1.6326e-04, -4.3063e-01]] ] light_map_data_dune10kt_1x2x6_AdjOpFlashMaxPELightMap: [ - ["Amplitude", [-1.1814e-07, 5.7685e-04, 2.8689e+00]], - ["Attenuation", [3.4830e-08, -1.3746e-04, 3.9619e-01]], - ["Correction", [6.1369e-08, -2.1512e-04, 1.2846e-01]] + ["Amplitude", [-1.1345e-07, 5.6632e-04, 2.8691e+00]], + ["Attenuation", [2.5484e-08, -1.2190e-04, 4.8002e-01]], + ["Correction", [4.7424e-08, -2.1438e-04, 2.7895e-01]] ] light_map_data_dune10kt_1x2x6_AdjOpFlashPELightMap: [ - ["Amplitude", [-1.2784e-07, 6.5108e-04, 2.3573e+00]], - ["Attenuation", [2.5445e-08, -1.2782e-04, 4.4181e-01]], - ["Correction", [-1.1812e-08, -1.6604e-05, -6.5942e-02]] + ["Amplitude", [-1.2095e-07, 6.4042e-04, 2.3421e+00]], + ["Attenuation", [1.0674e-08, -9.3331e-05, 5.2349e-01]], + ["Correction", [-4.1823e-08, 6.1425e-05, 2.3361e-02]] ] light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMinPELightMap: [ - ["Amplitude", [2.5276e-08, 2.9590e-04, -7.2231e-02]], - ["Attenuation", [-1.1732e-06, 8.1172e-03, -1.5276e+01]], - ["Correction", [7.6294e-07, -2.8500e-03, 1.9332e+00]] + ["Amplitude", [5.2745e-08, 1.9580e-04, -4.6346e-02]], + ["Attenuation", [4.0968e-07, 3.3389e-03, -1.2823e+01]], + ["Correction", [8.3383e-07, -3.1493e-03, 2.2614e+00]] ] light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMaxPELightMap: [ - ["Amplitude", [-8.8197e-08, 5.2212e-04, 2.1964e+00]], - ["Attenuation", [4.8050e-08, -9.5364e-05, -5.7437e-01]], - ["Correction", [3.7871e-07, -1.2702e-03, 7.0594e-01]] + ["Amplitude", [8.3383e-07, -3.1493e-03, 2.2614e+00]], + ["Attenuation", [5.9746e-08, -1.5539e-04, -5.1042e-01]], + ["Correction", [3.3565e-07, -1.1742e-03, 7.2596e-01]] ] light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashPELightMap: [ - ["Amplitude", [-1.2053e-07, 7.1549e-04, 7.2371e-01]], - ["Attenuation", [-7.9112e-08, 4.6725e-04, -1.7473e+00]], - ["Correction", [2.5340e-07, -1.0071e-03, 5.7146e-01]] + ["Amplitude", [-1.2489e-07, 7.2077e-04, 6.7834e-01]], + ["Attenuation", [-1.4162e-08, 2.2639e-04, -1.5893e+00]], + ["Correction", [2.4396e-07, -1.0174e-03, 6.4482e-01]] ] END_PROLOG \ No newline at end of file From 10f9ec96d5aaae278aed915dcbeb5fd3609578d3 Mon Sep 17 00:00:00 2001 From: mantheys Date: Fri, 9 Jan 2026 14:01:45 -0600 Subject: [PATCH 5/9] SAMLL IMPROVEMENTS --- dunereco/LowEUtils/LowEUtils.cc | 118 ++++++++++++++++++-------------- dunereco/LowEUtils/LowEUtils.h | 4 +- 2 files changed, 70 insertions(+), 52 deletions(-) diff --git a/dunereco/LowEUtils/LowEUtils.cc b/dunereco/LowEUtils/LowEUtils.cc index 10b6c152..0519575d 100644 --- a/dunereco/LowEUtils/LowEUtils.cc +++ b/dunereco/LowEUtils/LowEUtils.cc @@ -37,6 +37,7 @@ namespace lowe fAdjOpFlashMaxPELightMap(p.get>>>("AdjOpFlashMaxPELightMap", {})), // Light map file and histogram name for light map attenuation fAdjOpFlashPELightMap(p.get>>>("AdjOpFlashPELightMap", {})), // Light map file and histogram name for PE attenuation fFlashMatchBy(p.get("FlashMatchBy", "maximum")), // Method to match flashes ("maximum" or "light_map") + fFlashMatchByPELightMapExponent(p.get("FlashMatchByPELightMapExponent", 1)), // 0 implies matching against absolute PE error, 1 against relative PE error, 2 adds an additional weight to higher PE flashes. producer(new ProducerUtils(p)) { // Initialize the LowEUtils instance @@ -1463,60 +1464,60 @@ namespace lowe producer->ComputeDistanceX(clusterX, clusterTime, flashTime, driftLength, driftTime); if (fGeometry == "HD") { - // For DUNE 10kt geometry, we have different projections based on the plane - // Change the sign of clusterX for the collection plane - if (flash->Frame() == 0) // Collection plane - { - clusterX = -clusterX; // Convert to the collection plane coordinate system - } + // For DUNE 10kt geometry, we have different projections based on the plane + // Change the sign of clusterX for the collection plane + if (flash->Frame() == 0) // Collection plane + { + clusterX = -clusterX; // Convert to the collection plane coordinate system + } + if (pow(clusterY - flashY, 2) / pow(fAdjOpFlashY, 2) + pow(clusterZ - flashZ, 2) / pow(fAdjOpFlashZ, 2) > 1) + { + continue; + } + // flashR = sqrt(pow(clusterY - flashY, 2) + pow(clusterZ - flashZ, 2)); + } + else if (fGeometry == "VD") { + // Convert clusterX to the VD geometry [-driftLength/2, driftLength/2] + if (flash->Frame() == 0) {// Cathode flashes if (pow(clusterY - flashY, 2) / pow(fAdjOpFlashY, 2) + pow(clusterZ - flashZ, 2) / pow(fAdjOpFlashZ, 2) > 1) { continue; } // flashR = sqrt(pow(clusterY - flashY, 2) + pow(clusterZ - flashZ, 2)); - } - else if (fGeometry == "VD") { - // Convert clusterX to the VD geometry [-driftLength/2, driftLength/2] - if (flash->Frame() == 0) {// Cathode flashes - if (pow(clusterY - flashY, 2) / pow(fAdjOpFlashY, 2) + pow(clusterZ - flashZ, 2) / pow(fAdjOpFlashZ, 2) > 1) - { - continue; - } - // flashR = sqrt(pow(clusterY - flashY, 2) + pow(clusterZ - flashZ, 2)); } else if (flash->Frame() == 1 || flash->Frame() == 2) {// Membrane flashes - if (fAdjOpFlashMembraneProjection) { - if (clusterY * flashY < 0) {continue;} - if (pow(clusterX, 2) / pow(fAdjOpFlashX, 2) + pow(clusterZ - flashZ, 2) / pow(fAdjOpFlashZ, 2) > 1) { - continue; - } - } - else { - if (pow(clusterX, 2) / pow(fAdjOpFlashX, 2) + pow(clusterY - flashY, 2) / pow(fAdjOpFlashY, 2) + pow(clusterZ - flashZ, 2) / pow(fAdjOpFlashZ, 2) > 1) { - continue; - } - } - // flashR = sqrt(pow(clusterX, 2) + pow(clusterZ - flashZ, 2)); + if (fAdjOpFlashMembraneProjection) { + if (clusterY * flashY < 0) {continue;} + if (pow(clusterX, 2) / pow(fAdjOpFlashX, 2) + pow(clusterZ - flashZ, 2) / pow(fAdjOpFlashZ, 2) > 1) { + continue; + } + } + else { + if (pow(clusterX, 2) / pow(fAdjOpFlashX, 2) + pow(clusterY - flashY, 2) / pow(fAdjOpFlashY, 2) + pow(clusterZ - flashZ, 2) / pow(fAdjOpFlashZ, 2) > 1) { + continue; + } + } + // flashR = sqrt(pow(clusterX, 2) + pow(clusterZ - flashZ, 2)); } else if (flash->Frame() == 3 || flash->Frame() == 4) {// End-Cap flashes - if (fAdjOpFlashEndCapProjection) { - if (clusterZ < fidVolZ/2 && flash->Frame() == 3) {continue;} // Skip if the cluster is in the bottom half and the flash is in the top end-cap - if (clusterZ > fidVolZ/2 && flash->Frame() == 4) {continue;} // Skip if the cluster is in the top half and the flash is in the bottom end-cap - if (pow(clusterX, 2) / pow(fAdjOpFlashX, 2) + pow(clusterY - flashY, 2) / pow(fAdjOpFlashY, 2) > 1) { - continue; - } - } - else { - if (pow(clusterX, 2) / pow(fAdjOpFlashX, 2) + pow(clusterY - flashY, 2) / pow(fAdjOpFlashY, 2) + pow(clusterZ - flashZ, 2) / pow(fAdjOpFlashZ, 2) > 1) { - continue; - } - } - // flashR = sqrt(pow(clusterX, 2) + pow(clusterY - flashY, 2) + pow(clusterZ - flashZ, 2)); + if (fAdjOpFlashEndCapProjection) { + if (clusterZ < fidVolZ/2 && flash->Frame() == 3) {continue;} // Skip if the cluster is in the bottom half and the flash is in the top end-cap + if (clusterZ > fidVolZ/2 && flash->Frame() == 4) {continue;} // Skip if the cluster is in the top half and the flash is in the bottom end-cap + if (pow(clusterX, 2) / pow(fAdjOpFlashX, 2) + pow(clusterY - flashY, 2) / pow(fAdjOpFlashY, 2) > 1) { + continue; + } + } + else { + if (pow(clusterX, 2) / pow(fAdjOpFlashX, 2) + pow(clusterY - flashY, 2) / pow(fAdjOpFlashY, 2) + pow(clusterZ - flashZ, 2) / pow(fAdjOpFlashZ, 2) > 1) { + continue; + } + } + // flashR = sqrt(pow(clusterX, 2) + pow(clusterY - flashY, 2) + pow(clusterZ - flashZ, 2)); } clusterX = driftLength / 2 - clusterX; // Convert to the VD geometry coordinate system } else { - continue; // Unknown geometry, skip this matching + continue; // Unknown geometry, skip this matching } if (flashPE > matchedFlashPE || matchedFlashIndex == -1) { @@ -1670,11 +1671,17 @@ namespace lowe float &b, float &c) /* - Store light map parameters for flash selection - - LightMapType: type of light map ("min", "max", "med") - - a: Amplitude parameters - - b: Attenuation parameters - - c: Correction parameters + Retrieve light map parameters for flash selection + - LightMapType: type of light map to use ("min", "max", "med"), selecting which set + of polynomial coefficients to apply. + - ClusterCharge: reconstructed charge of the cluster used to evaluate the light map + parameter polynomials. + - a: amplitude parameter (output), computed from the "Amplitude" light map polynomial + evaluated at ClusterCharge and used to scale the predicted flash amplitude. + - b: attenuation parameter (output), computed from the "Attenuation" light map + polynomial evaluated at ClusterCharge and used to describe light attenuation. + - c: correction parameter (output), computed from the "Correction" light map + polynomial evaluated at ClusterCharge and used as an additional correction factor. */ { std::vector>> SelectedLightMap; @@ -1705,6 +1712,7 @@ namespace lowe } bool LowEUtils::SelectPDSFlash( + const bool &IsFirstFlash, const float &TPCDriftTime, const float &ClusterTime, const float &ClusterCharge, @@ -1718,18 +1726,26 @@ namespace lowe - ClusterTime: time of the cluster - ClusterCharge: charge of the cluster - RefOpFlashTime: time of the reference flash - - RefOpFlashTime: time of the reference flash - - OpFlashPE: total PE of the flash + - RefOpFlashPE: total PE of the reference flash + - OpFlashTime: time of the flash to evaluate - OpFlashPE: total PE of the flash */ { + // If this is the first flash, always select it + if (IsFirstFlash) { + return true; + } if (fFlashMatchBy == "light_map") { float a, b, c; GetLightMapParameters("med", ClusterCharge, a, b, c); - double RefPE = pow(10, a - a * b * std::abs(ClusterTime - RefOpFlashTime) / TPCDriftTime + c * pow(std::abs(ClusterTime - RefOpFlashTime) / TPCDriftTime, 2)); - double SelectedPE = pow(10, a - a * b * std::abs(ClusterTime - OpFlashTime) / TPCDriftTime + c * pow(std::abs(ClusterTime - OpFlashTime) / TPCDriftTime, 2)); + double RefPE = pow(10, a - a * b * (ClusterTime - RefOpFlashTime) / TPCDriftTime + c * pow((ClusterTime - RefOpFlashTime) / TPCDriftTime, 2)); + double SelectedPE = pow(10, a - a * b * (ClusterTime - OpFlashTime) / TPCDriftTime + c * pow((ClusterTime - OpFlashTime) / TPCDriftTime, 2)); // If the flash PE is closer to the selected PE, than the reference flash PE, select it - if (std::abs(OpFlashPE - SelectedPE) < std::abs(RefOpFlashPE - RefPE)) { + float SelectedError = std::abs((OpFlashPE - SelectedPE) / pow(OpFlashPE, fFlashMatchByPELightMapExponent)); + float RefError = std::abs((RefOpFlashPE - RefPE) / pow(RefOpFlashPE, fFlashMatchByPELightMapExponent)); + if (SelectedError < RefError) { + // std::cout << "Selecting flash with PE " << OpFlashPE << " closer to predicted PE " << SelectedPE << " than reference flash PE " << RefOpFlashPE << " predicted PE " << RefPE << std::endl; + // std::cout << "The differences are " << SelectedError << " and " << RefError << std::endl; return true; } else { diff --git a/dunereco/LowEUtils/LowEUtils.h b/dunereco/LowEUtils/LowEUtils.h index 3a8e504d..fbd8fe5d 100644 --- a/dunereco/LowEUtils/LowEUtils.h +++ b/dunereco/LowEUtils/LowEUtils.h @@ -254,6 +254,7 @@ namespace lowe float &c); bool SelectPDSFlash( + const bool &IsFirstFlash, const float &TPCDriftTime, const float &ClusterTime, const float &ClusterCharge, @@ -296,7 +297,8 @@ namespace lowe const std::vector>> fAdjOpFlashMinPELightMap; // Light map file and histogram name for light map attenuation const std::vector>> fAdjOpFlashMaxPELightMap; // Light map file and histogram name for light map attenuation const std::vector>> fAdjOpFlashPELightMap; // Light map file and histogram name for PE attenuation - const std::string fFlashMatchBy; // Method to match flashes ("maximum" or "light_map") + const std::string fFlashMatchBy; // Method to match flashes ("maximum" or "light_map") + const int fFlashMatchByPELightMapExponent; // Exponent for PE weighting in light map flash matching std::unique_ptr producer; // Pointer to the ProducerUtils instance }; } // namespace lowe From 71d933212dd3ffd9268f8987050a7b81f02638b8 Mon Sep 17 00:00:00 2001 From: mantheys Date: Fri, 9 Jan 2026 14:01:53 -0600 Subject: [PATCH 6/9] ADD DEFAULTS --- dunereco/LowEUtils/fcl/SolarEvent.fcl | 1 + 1 file changed, 1 insertion(+) diff --git a/dunereco/LowEUtils/fcl/SolarEvent.fcl b/dunereco/LowEUtils/fcl/SolarEvent.fcl index 1b45e950..f4709219 100644 --- a/dunereco/LowEUtils/fcl/SolarEvent.fcl +++ b/dunereco/LowEUtils/fcl/SolarEvent.fcl @@ -30,6 +30,7 @@ solar_event_dune10kt_1x2x6: AdjOpFlashMaxPELightMap: @local::light_map_data_dune10kt_1x2x6_AdjOpFlashMaxPELightMap AdjOpFlashPELightMap: @local::light_map_data_dune10kt_1x2x6_AdjOpFlashPELightMap FlashMatchBy: "light_map" # Match flashes by residual. Alternative is to match by MaxFlashPE. + FlashMatchByPELightMapExponent: 1 # Exponent for PE light map matching. Debug: true } From 60a0474891a05ed1ea82c86e8656a6c1bc6b23e1 Mon Sep 17 00:00:00 2001 From: mantheys Date: Sat, 10 Jan 2026 07:13:19 -0600 Subject: [PATCH 7/9] FIX TYPO --- dunereco/LowEUtils/fcl/SolarEvent.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dunereco/LowEUtils/fcl/SolarEvent.fcl b/dunereco/LowEUtils/fcl/SolarEvent.fcl index f4709219..9502e986 100644 --- a/dunereco/LowEUtils/fcl/SolarEvent.fcl +++ b/dunereco/LowEUtils/fcl/SolarEvent.fcl @@ -1,4 +1,4 @@ -#include LightMapData.fcl +#include "LightMapData.fcl" BEGIN_PROLOG From 433fe4bf653663c81b792379312727e567c22de1 Mon Sep 17 00:00:00 2001 From: mantheys Date: Sat, 10 Jan 2026 07:23:50 -0600 Subject: [PATCH 8/9] GENERALIZE VARIABLE TYPE --- dunereco/LowEUtils/LowEUtils.cc | 2 +- dunereco/LowEUtils/LowEUtils.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/dunereco/LowEUtils/LowEUtils.cc b/dunereco/LowEUtils/LowEUtils.cc index 0519575d..4ce6c9bd 100644 --- a/dunereco/LowEUtils/LowEUtils.cc +++ b/dunereco/LowEUtils/LowEUtils.cc @@ -37,7 +37,7 @@ namespace lowe fAdjOpFlashMaxPELightMap(p.get>>>("AdjOpFlashMaxPELightMap", {})), // Light map file and histogram name for light map attenuation fAdjOpFlashPELightMap(p.get>>>("AdjOpFlashPELightMap", {})), // Light map file and histogram name for PE attenuation fFlashMatchBy(p.get("FlashMatchBy", "maximum")), // Method to match flashes ("maximum" or "light_map") - fFlashMatchByPELightMapExponent(p.get("FlashMatchByPELightMapExponent", 1)), // 0 implies matching against absolute PE error, 1 against relative PE error, 2 adds an additional weight to higher PE flashes. + fFlashMatchByPELightMapExponent(p.get("FlashMatchByPELightMapExponent", 1)), // 0 implies matching against absolute PE error, 1 against relative PE error, 2 adds an additional weight to higher PE flashes. producer(new ProducerUtils(p)) { // Initialize the LowEUtils instance diff --git a/dunereco/LowEUtils/LowEUtils.h b/dunereco/LowEUtils/LowEUtils.h index fbd8fe5d..b77ffe9e 100644 --- a/dunereco/LowEUtils/LowEUtils.h +++ b/dunereco/LowEUtils/LowEUtils.h @@ -298,7 +298,7 @@ namespace lowe const std::vector>> fAdjOpFlashMaxPELightMap; // Light map file and histogram name for light map attenuation const std::vector>> fAdjOpFlashPELightMap; // Light map file and histogram name for PE attenuation const std::string fFlashMatchBy; // Method to match flashes ("maximum" or "light_map") - const int fFlashMatchByPELightMapExponent; // Exponent for PE weighting in light map flash matching + const float fFlashMatchByPELightMapExponent; // Exponent for PE weighting in light map flash matching std::unique_ptr producer; // Pointer to the ProducerUtils instance }; } // namespace lowe From d8d82a4850eb9adcb1b6401a86cc0fa12ccf3b6c Mon Sep 17 00:00:00 2001 From: mantheys Date: Wed, 21 Jan 2026 00:28:59 +0100 Subject: [PATCH 9/9] CORRECT WRONG PARAMETERS --- dunereco/LowEUtils/fcl/LightMapData.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dunereco/LowEUtils/fcl/LightMapData.fcl b/dunereco/LowEUtils/fcl/LightMapData.fcl index 793d6290..ead8cd1a 100644 --- a/dunereco/LowEUtils/fcl/LightMapData.fcl +++ b/dunereco/LowEUtils/fcl/LightMapData.fcl @@ -25,7 +25,7 @@ light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMinPELightMap: [ ] light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashMaxPELightMap: [ - ["Amplitude", [8.3383e-07, -3.1493e-03, 2.2614e+00]], + ["Amplitude", [-1.1320e-07, 5.9340e-04, 2.1169e+00]], ["Attenuation", [5.9746e-08, -1.5539e-04, -5.1042e-01]], ["Correction", [3.3565e-07, -1.1742e-03, 7.2596e-01]] ]