diff --git a/dunereco/LowEUtils/LowEUtils.cc b/dunereco/LowEUtils/LowEUtils.cc index 0ad2adb4..4ce6c9bd 100644 --- a/dunereco/LowEUtils/LowEUtils.cc +++ b/dunereco/LowEUtils/LowEUtils.cc @@ -22,21 +22,36 @@ 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 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 + 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") - fAdjOpFlashMinPEAttenuationStrength(p.get("AdjOpFlashMinPEAttenuationStrength", 10)), // Strength of the asymptotic attenuation for minimum PE cut (in powers of 10) + fAdjOpFlashMaxPEAttenuate(p.get("AdjOpFlashMaxPEAttenuate", "flat")), // Type of attenuation for maximum PE cut ("light_map", "asymptotic", "linear" or "flat") + 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") + 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 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; @@ -171,9 +186,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 @@ -323,7 +344,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 @@ -425,6 +449,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 @@ -476,16 +512,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 @@ -513,6 +558,7 @@ namespace lowe return RecoY; } + //...................................................... void LowEUtils::FillClusterVariables( std::vector>> Clusters, std::vector> &ClNHits, @@ -520,6 +566,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(); @@ -542,6 +595,7 @@ namespace lowe } // End of loop over clusters } // End of loop over planes } + void LowEUtils::FillClusterVariables( std::set SignalTrackIDs, std::vector>> Clusters, @@ -558,6 +612,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(); @@ -631,25 +701,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 = {{}, {}, {}}; @@ -886,13 +976,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 = {{}, {}, {}}; @@ -1018,7 +1116,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) { @@ -1036,7 +1139,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. @@ -1256,12 +1374,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 @@ -1314,84 +1442,82 @@ 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(); 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) { @@ -1400,6 +1526,7 @@ namespace lowe matchedFlashIndex = i; matchedRecoX = clusterX; } + } } if (debug) { if (matchedFlashIndex != -1) { @@ -1421,20 +1548,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; } @@ -1457,5 +1610,160 @@ namespace lowe else { return true; } } return true; - } // SelectPDSFlashPE + } + + + 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") { + 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; } + } + 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; + } + + void LowEUtils::GetLightMapParameters( + const std::string &LightMapType, + const float &ClusterCharge, + float &a, + float &b, + float &c) + /* + 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; + 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 bool &IsFirstFlash, + 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 + - 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 * (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 + 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 { + 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 6b92fed6..b77ffe9e 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, @@ -230,6 +234,35 @@ 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); + + void GetLightMapParameters( + const std::string &LightMapType, + const float &ClusterCharge, + float &a, + float &b, + float &c); + + bool SelectPDSFlash( + const bool &IsFirstFlash, + 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 @@ -250,14 +283,22 @@ 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 fAdjOpFlashMinPECut; // Minimum photoelectrons for adjacent flash 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 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 + 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 float fFlashMatchByPELightMapExponent; // Exponent for PE weighting in light map flash matching std::unique_ptr producer; // Pointer to the ProducerUtils instance }; } // namespace lowe diff --git a/dunereco/LowEUtils/fcl/LightMapData.fcl b/dunereco/LowEUtils/fcl/LightMapData.fcl new file mode 100644 index 00000000..ead8cd1a --- /dev/null +++ b/dunereco/LowEUtils/fcl/LightMapData.fcl @@ -0,0 +1,39 @@ +BEGIN_PROLOG + +light_map_data_dune10kt_1x2x6_AdjOpFlashMinPELightMap: [ + ["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.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.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", [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", [-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]] + ] + +light_map_data_dunevd10kt_1x8x14_3view_30deg_AdjOpFlashPELightMap: [ + ["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 diff --git a/dunereco/LowEUtils/fcl/SolarEvent.fcl b/dunereco/LowEUtils/fcl/SolarEvent.fcl index e2509efb..9502e986 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,14 @@ 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. + FlashMatchByPELightMapExponent: 1 # Exponent for PE light map matching. Debug: true } @@ -37,11 +38,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