diff --git a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc index 1515339f..7a519418 100644 --- a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc +++ b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc @@ -5,10 +5,14 @@ using namespace producer; namespace solar { AdjOpHitsUtils::AdjOpHitsUtils(fhicl::ParameterSet const &p) - : fOpHitTimeVariable(p.get("OpHitTimeVariable", "PeakTime")), // Variable to use for time sorting ("StartTime" or "PeakTime") + : fOpWaveformLabel(p.get("OpWaveformLabel", "opdetwaveform")), // Label for OpDetWaveform collection. + fOpHitLabel(p.get("OpHitLabel", "ophit")), // Label for OpHit collection. + fOpHitTimeVariable(p.get("OpHitTimeVariable", "PeakTime")), // Variable to use for time sorting ("StartTime" or "PeakTime"). + fOpHitTime2us(p.get("OpHitTime2us", false)), // Conversion factor from OpHit time units to microseconds. Default factor is TickPeriod() from DetectorClocksService. fOpFlashAlgoNHit(p.get("OpFlashAlgoNHit", 3)), // Minimum number of OpHits in a cluster to consider it for flash creation. - fOpFlashAlgoMinTime(p.get("OpFlashAlgoMinTime", 0.010)), // Negative time window to look for adj. OpHits. Default for HD 10 ns [0.6 tick] - fOpFlashAlgoMaxTime(p.get("OpFlashAlgoMaxTime", 0.016)), // Positive time window to look for adj. OpHits. Default for HD 16 ns [1 tick] + fOpFlashAlgoMinTime(p.get("OpFlashAlgoMinTime", 0.32)), // Negative time window to look for adj. OpHits in [us]. + fOpFlashAlgoMaxTime(p.get("OpFlashAlgoMaxTime", 0.96)), // Positive time window to look for adj. OpHits in [us]. + fOpFlashAlgoWeightedTime(p.get("OpFlashAlgoWeightedTime", false)), // Whether to use weighted time of trigger time for flash time reference. fOpFlashAlgoRad(p.get("OpFlashAlgoRad", 300.0)), // Radius to look for adj. OpHits in [cm] units. fOpFlashAlgoPE(p.get("OpFlashAlgoPE", 1.5)), // Minimum PE of OpHit to consider it for flash creation. fOpFlashAlgoTriggerPE(p.get("OpFlashAlgoTriggerPE", 1.5)), // Minimum PE of OpHit to consider it as a trigger for flash creation. @@ -20,34 +24,36 @@ namespace solar fXAStartCapZ(p.get("XAStartCapZ", -96.5)) { } + + void AdjOpHitsUtils::MakeFlashVector(std::vector &FlashVec, std::vector>> &OpHitClusters, art::Event const &evt) { auto const clockData = art::ServiceHandle()->DataFor(evt); auto TickPeriod = clockData.OpticalClock().TickPeriod(); - if (fOpFlashAlgoHotVertexThld > 1) - { + if (fOpFlashAlgoHotVertexThld > 1) { ProducerUtils::PrintInColor("Hot vertex threshold must be between 0 and 1", ProducerUtils::GetColor("red"), "Error"); return; } + for (std::vector> Cluster : OpHitClusters) { - if (!Cluster.empty()) + if (Cluster.empty()) { - if (fOpHitTimeVariable == "StartTime") - std::stable_sort(Cluster.begin(), Cluster.end(), [](art::Ptr a, art::Ptr b) - { return a->StartTime() < b->StartTime(); }); - else // Default to PeakTime - std::stable_sort(Cluster.begin(), Cluster.end(), [](art::Ptr a, art::Ptr b) - { return a->PeakTime() < b->PeakTime(); }); + continue; } + int Plane = GetOpHitPlane(Cluster[0], 0.1); + int MaxIdx = 0; + int Idx = 0; int NHit = 0; - double Time = -1e6; + double TimeMax = -1e6; double TimeWidth = 0; + double TimeWeighted = -1e6; double TimeSum = 0; + double Amplitude = 0; double PE = 0; double MaxPE = 0; - std::vector PEperOpDet; + std::vector PEperOpDet = {}; double FastToTotal = 1; double X = 0; double Y = 0; @@ -59,30 +65,62 @@ namespace solar double YSum = 0; double ZSum = 0; double STD = 0; + std::vector MainOpWaveform = {}; // Make vector for main OpWaveform with 1000 entries (max waveform size) + bool MainOpWaveformValid = false; + float MainOpWaveformTime = -1e6; + + std::vector OpHitWvfValid = {}; + std::vector OpHitWvfTime = {}; + std::vector> OpHitWvfIntVector = {}; + GetOpHitSignal(Cluster, OpHitWvfIntVector, OpHitWvfTime, OpHitWvfValid, evt); // Get OpHit waveforms // Compute total number of PE and MaxPE. for (art::Ptr PDSHit : Cluster) { NHit++; + double thisTime = -1e6; + double thisPE = PDSHit->PE(); + double thisAmp = PDSHit->Amplitude(); auto thisPlane = GetOpHitPlane(PDSHit, 0.1); + if (thisPlane != Plane) { ProducerUtils::PrintInColor("OpHit in cluster not in same plane: CH " + ProducerUtils::str(PDSHit->OpChannel()) + " Plane " + ProducerUtils::str(thisPlane) + " Expected Plane " + ProducerUtils::str(Plane), ProducerUtils::GetColor("red"), "Error"); Plane = -1; // Set plane to -1 if hits in cluster are not in the same plane } - PE += PDSHit->PE(); - if (PDSHit->PE() > MaxPE) - MaxPE = PDSHit->PE(); + if (fOpHitTimeVariable == "StartTime") { + thisTime = PDSHit->StartTime(); // Use StartTime + } + else { + thisTime = PDSHit->PeakTime(); // Default to PeakTime + } + if (fOpHitTime2us) { + thisTime *= TickPeriod; // Convert to microseconds + } + + PE += thisPE; + TimeSum += thisTime * thisPE; + + if (thisPE > MaxPE) { + Amplitude = thisAmp; + MaxPE = thisPE; + MaxIdx = Idx; + TimeMax = thisTime; + } + + PEperOpDet.push_back(thisPE); + Idx++; + } - PEperOpDet.push_back(PDSHit->PE()); - if (fOpHitTimeVariable == "StartTime") - TimeSum += PDSHit->StartTime() * TickPeriod * PDSHit->PE(); - else // Default to PeakTime - TimeSum += PDSHit->PeakTime() * TickPeriod * PDSHit->PE(); + if (OpHitWvfValid[MaxIdx]) { + MainOpWaveform = OpHitWvfIntVector[MaxIdx]; + MainOpWaveformTime = OpHitWvfTime[MaxIdx]; + MainOpWaveformValid = true; } float HotPE = 0; - Time = TimeSum / PE; + TimeWeighted = TimeSum / PE; + // Compute flash center from weighted average of "hottest" ophits. for (art::Ptr PDSHit : Cluster) { @@ -94,6 +132,7 @@ namespace solar HotPE += PDSHit->PE(); } } + X = XSum / HotPE; Y = YSum / HotPE; Z = ZSum / HotPE; @@ -104,12 +143,18 @@ namespace solar std::vector varXZ; for (art::Ptr PDSHit : Cluster) { + float ThisOpHitTime = -1e6; auto OpHitXYZ = wireReadout.OpDetGeoFromOpChannel(PDSHit->OpChannel()).GetCenter(); + if (fOpHitTimeVariable == "StartTime") - TimeWidth += (PDSHit->StartTime() * TickPeriod - Time) * (PDSHit->StartTime() * TickPeriod - Time); + ThisOpHitTime = PDSHit->StartTime(); else // Default to PeakTime - TimeWidth += (PDSHit->PeakTime() * TickPeriod - Time) * (PDSHit->PeakTime() * TickPeriod - Time); + ThisOpHitTime = PDSHit->PeakTime(); + if (fOpHitTime2us) { + ThisOpHitTime *= TickPeriod; + } + TimeWidth += (TimeMax - ThisOpHitTime) * (TimeMax - ThisOpHitTime); XWidth += (OpHitXYZ.X() - X) * (OpHitXYZ.X() - X); YWidth += (OpHitXYZ.Y() - Y) * (OpHitXYZ.Y() - Y); ZWidth += (OpHitXYZ.Z() - Z) * (OpHitXYZ.Z() - Z); @@ -128,21 +173,48 @@ namespace solar // Compute FastToTotal according to the #PEs arriving within the first 10% of the time window wrt the total #PEs for (art::Ptr PDSHit : Cluster) { - auto thisTime = 0.0; + auto thisTime = -1e6; if (fOpHitTimeVariable == "StartTime") - thisTime = PDSHit->StartTime() * TickPeriod; + thisTime = PDSHit->StartTime(); else // Default to PeakTime - thisTime = PDSHit->PeakTime() * TickPeriod; + thisTime = PDSHit->PeakTime(); + if (fOpHitTime2us) { + thisTime *= TickPeriod; + } - if (thisTime < Time + TimeWidth / 10) + // Check if thisTime is within 10% of the time window in positive and negative direction + if (std::abs(thisTime - TimeMax) <= 0.05 * TimeWidth) FastToTotal += PDSHit->PE(); } FastToTotal /= PE; - FlashVec.push_back(FlashInfo{Plane, NHit, Time, TimeWidth, PE, MaxPE, PEperOpDet, FastToTotal, X, Y, Z, XWidth, YWidth, ZWidth, STD}); + FlashVec.push_back(FlashInfo{ + Plane, + NHit, + MaxPE, + TimeMax, + Amplitude, + TimeWidth, + TimeWeighted, + PE, + PEperOpDet, + FastToTotal, + X, + Y, + Z, + XWidth, + YWidth, + ZWidth, + STD, + MainOpWaveform, + MainOpWaveformTime, + MainOpWaveformValid + } + ); } return; } + float AdjOpHitsUtils::GetOpFlashPlaneSTD(const int Plane, const std::vector varXY, const std::vector varYZ, const std::vector varXZ) { std::vector var; @@ -185,6 +257,7 @@ namespace solar return varstd; } + void AdjOpHitsUtils::CalcAdjOpHits(const std::vector> &OpHitVector, std::vector>> &OpHitClusters, std::vector> &OpHitClusterIdx, art::Event const &evt) { // This is the low energy flash (hit clustering) algorithm aka flip-flop. It is based on the idea that the hits in the same flash are close in time and space and follow a 1/r² signal decay. @@ -222,14 +295,21 @@ namespace solar const auto &hit = OpHitVector[*it]; if (hit->PE() < fOpFlashAlgoTriggerPE) continue; - + float hitTime = -1e6; + if (fOpHitTimeVariable == "StartTime") + hitTime = hit->StartTime(); + else // Default to PeakTime + hitTime = hit->PeakTime(); + if (fOpHitTime2us) { + hitTime *= TickPeriod; + } bool main_hit = true; // If a trigger hit is found, start a new cluster with the hits around it that are within the time and radius range ClusteredHits[*it] = true; AdjHitIdx.push_back(*it); AdjHitVec.push_back(hit); - sOpHitClustering += "Trigger hit found: PE " + ProducerUtils::str(hit->PE()) + " CH " + ProducerUtils::str(hit->OpChannel()) + " Time " + ProducerUtils::str(hit->PeakTime() * TickPeriod) + "\n"; + sOpHitClustering += "Trigger hit found: PE " + ProducerUtils::str(hit->PE()) + " CH " + ProducerUtils::str(hit->OpChannel()) + " Time " + ProducerUtils::str(hitTime) + "\n"; int refHit1 = hit->OpChannel(); auto ref1 = wireReadout.OpDetGeoFromOpChannel(refHit1).GetCenter(); @@ -255,9 +335,13 @@ namespace solar thisHitTime = hit->PeakTime(); thisAdjHitTime = adjHit->PeakTime(); } + if (fOpHitTime2us) { + thisHitTime *= TickPeriod; + thisAdjHitTime *= TickPeriod; + } - if (std::abs(thisAdjHitTime - thisHitTime) * TickPeriod > fOpFlashAlgoMaxTime) { - sOpHitClustering += "Breaking time loop at dT " + ProducerUtils::str(std::abs(thisAdjHitTime - thisHitTime) * TickPeriod) + " us\n"; + if (std::abs(thisAdjHitTime - thisHitTime) > fOpFlashAlgoMaxTime) { + sOpHitClustering += "Breaking time loop at dT " + ProducerUtils::str(std::abs(thisAdjHitTime - thisHitTime)) + " us\n"; break; } @@ -271,7 +355,7 @@ namespace solar // If hit has already been clustered, skip if (ClusteredHits[*it2] == true && !fOpFlashAlgoHitDuplicates) { - sOpHitClustering += "Skipping already clustered hit: CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(adjHit->PeakTime() * TickPeriod) + "\n"; + sOpHitClustering += "Skipping already clustered hit: CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(thisAdjHitTime) + "\n"; continue; } @@ -281,20 +365,27 @@ namespace solar { if (adjHit->PE() > hit->PE()) { - sOpHitClustering += "¡¡¡Hit with PE > TriggerPE found: PE " + ProducerUtils::str(adjHit->PE()) + " CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(adjHit->PeakTime() * TickPeriod) + "\n"; + sOpHitClustering += "¡¡¡Hit with PE > TriggerPE found: PE " + ProducerUtils::str(adjHit->PE()) + " CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(thisAdjHitTime) + "\n"; ClusteredHits[*it] = false; main_hit = false; // Reset the ClusteredHits values for the hits that have been added to the cluster for (auto it3 = AdjHitIdx.begin(); it3 != AdjHitIdx.end(); ++it3) { - sOpHitClustering += "---Removing hit: CH " + ProducerUtils::str(OpHitVector[*it3]->OpChannel()) + " Time " + ProducerUtils::str(OpHitVector[*it3]->PeakTime() * TickPeriod) + "\n"; + auto refTime = OpHitVector[*it3]->PeakTime(); + auto refChannel = OpHitVector[*it3]->OpChannel(); + if (fOpHitTimeVariable == "StartTime") + refTime = OpHitVector[*it3]->StartTime(); + if (fOpHitTime2us) { + refTime *= TickPeriod; + } + sOpHitClustering += "---Removing hit: CH " + ProducerUtils::str(refChannel) + " Time " + ProducerUtils::str(refTime) + "\n"; ClusteredHits[*it3] = false; } break; } else { - sOpHitClustering += "+++Adding hit: PE " + ProducerUtils::str(adjHit->PE()) + " CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(adjHit->PeakTime() * TickPeriod) + "\n"; + sOpHitClustering += "+++Adding hit: PE " + ProducerUtils::str(adjHit->PE()) + " CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(thisAdjHitTime) + "\n"; AdjHitVec.push_back(adjHit); AdjHitIdx.push_back(*it2); ClusteredHits[*it2] = true; @@ -330,8 +421,12 @@ namespace solar thisHitTime = hit->PeakTime(); thisAdjHitTime = adjHit->PeakTime(); } + if (fOpHitTime2us) { + thisHitTime *= TickPeriod; + thisAdjHitTime *= TickPeriod; + } - if (std::abs(thisHitTime - thisAdjHitTime) * TickPeriod > fOpFlashAlgoMinTime) + if (std::abs(thisHitTime - thisAdjHitTime) > fOpFlashAlgoMinTime) break; if (adjHit->PE() < fOpFlashAlgoPE) @@ -344,7 +439,7 @@ namespace solar // if hit has already been clustered, skip if (ClusteredHits[*it4] == true && !fOpFlashAlgoHitDuplicates) { - sOpHitClustering += "Skipping already clustered hit: CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(adjHit->PeakTime() * TickPeriod) + "\n"; + sOpHitClustering += "Skipping already clustered hit: CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(thisAdjHitTime) + "\n"; continue; } @@ -354,20 +449,27 @@ namespace solar { if (adjHit->PE() > hit->PE()) { - sOpHitClustering += "¡¡¡Hit with PE > TriggerPE found: PE " + ProducerUtils::str(adjHit->PE()) + " CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(adjHit->PeakTime() * TickPeriod) + "\n"; + sOpHitClustering += "¡¡¡Hit with PE > TriggerPE found: PE " + ProducerUtils::str(adjHit->PE()) + " CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(thisAdjHitTime) + "\n"; ClusteredHits[*it] = false; main_hit = false; // Reset the ClusteredHits values for the hits that have been added to the cluster for (auto it5 = AdjHitIdx.begin(); it5 != AdjHitIdx.end(); ++it5) { - sOpHitClustering += "---Removing hit: CH " + ProducerUtils::str(OpHitVector[*it5]->OpChannel()) + " Time " + ProducerUtils::str(OpHitVector[*it5]->PeakTime() * TickPeriod) + "\n"; + auto refTime = OpHitVector[*it5]->PeakTime(); + auto refChannel = OpHitVector[*it5]->OpChannel(); + if (fOpHitTimeVariable == "StartTime") + refTime = OpHitVector[*it5]->StartTime(); + if (fOpHitTime2us) { + refTime *= TickPeriod; + } + sOpHitClustering += "---Removing hit: CH " + ProducerUtils::str(refChannel) + " Time " + ProducerUtils::str(refTime) + "\n"; ClusteredHits[*it5] = false; } break; } else { - sOpHitClustering += "Adding hit: PE " + ProducerUtils::str(adjHit->PE()) + " CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(adjHit->PeakTime() * TickPeriod) + "\n"; + sOpHitClustering += "Adding hit: PE " + ProducerUtils::str(adjHit->PE()) + " CH " + ProducerUtils::str(adjHit->OpChannel()) + " Time " + ProducerUtils::str(thisAdjHitTime) + "\n"; AdjHitVec.push_back(adjHit); AdjHitIdx.push_back(*it4); ClusteredHits[*it4] = true; @@ -390,6 +492,7 @@ namespace solar return; } + void AdjOpHitsUtils::FlashMatchResidual(float &Residual, std::vector> Hits, double x, double y, double z) { if (Hits.size() == 0) @@ -451,6 +554,7 @@ namespace solar return; } + int AdjOpHitsUtils::GetOpHitPlane(const art::Ptr &hit, float buffer) { std::string geoName = geom->DetectorName(); @@ -483,9 +587,9 @@ namespace solar } } - // Define a function and map to get the plane of the hit. If geometry is VD the number of planes is 5 (cathode, leftmembrane, rightmembrane, startcap, finalcap). If geometry is HD the number of planes is 2 (leftdrift, rightdrift). + std::map AdjOpHitsUtils::GetOpHitPlaneMap(const std::vector> &OpHitVector) - { + { // Define a function and map to get the plane of the hit. If geometry is VD the number of planes is 5 (cathode, leftmembrane, rightmembrane, startcap, finalcap). If geometry is HD the number of planes is 2 (leftdrift, rightdrift). std::map OpHitPlane; for (const auto &hit : OpHitVector) { @@ -494,9 +598,9 @@ namespace solar return OpHitPlane; } - // Define a function to check if the adjhit is in the same plane as the reference hit using the OpHitPlane map + bool AdjOpHitsUtils::CheckOpHitPlane(std::map OpHitPlane, int refHit, int adjHit) - { + { // Define a function to check if the adjhit is in the same plane as the reference hit using the OpHitPlane map // Check that both hits are in the map if (OpHitPlane.find(refHit) == OpHitPlane.end() || OpHitPlane.find(adjHit) == OpHitPlane.end()) { ProducerUtils::PrintInColor("Hit not found in OpHitPlane map: refHit " + ProducerUtils::str(refHit) + " adjHit " + ProducerUtils::str(adjHit), ProducerUtils::GetColor("red"), "Error"); @@ -509,4 +613,157 @@ namespace solar return false; } + + void AdjOpHitsUtils::GetOpHitWaveforms(const std::vector> &OpHitVector, std::vector> &OpHitWvfVector, std::vector &OpHitWvfValid, art::Event const &evt) + { // Define a function to get the OpHit waveforms that correspond to each OpHit in the input vector + // Get the OpDetWaveform handle + art::Handle> opDetWvfHandle; + evt.getByLabel(fOpWaveformLabel, opDetWvfHandle); + auto const clockData = art::ServiceHandle()->DataFor(evt); + if (!opDetWvfHandle.isValid()) + { + ProducerUtils::PrintInColor("Invalid OpDetWaveform handle", ProducerUtils::GetColor("red"), "Error"); + for (size_t i = 0; i < OpHitVector.size(); i++) + { + OpHitWvfVector.push_back(art::Ptr()); // Fill with a null pointer to keep the indices consistent + OpHitWvfValid.push_back(false); + } + return; + } + + // Loop over the OpHit vector and get the corresponding OpDetWaveform + for (const auto &hit : OpHitVector) + { + bool found = false; + float hitTime = -1e6; + unsigned int opChannel = hit->OpChannel(); + auto TickPeriod = clockData.OpticalClock().TickPeriod(); + + if (fOpHitTimeVariable == "StartTime") + hitTime = hit->StartTime(); + else // Default to PeakTime + hitTime = hit->PeakTime(); + if (fOpHitTime2us) { + hitTime *= TickPeriod; + } + + for (size_t i = 0; i < opDetWvfHandle->size(); i++) + { // Match by channel number and amplitude + if ((*opDetWvfHandle)[i].ChannelNumber() == opChannel && hitTime >= (*opDetWvfHandle)[i].TimeStamp() && hitTime <= (*opDetWvfHandle)[i].TimeStamp() + (*opDetWvfHandle)[i].Waveform().size() * TickPeriod) + { + ProducerUtils::PrintInColor("Matching OpDetWaveform found for OpHit channel " + ProducerUtils::str(int(opChannel)), ProducerUtils::GetColor("green"), "Debug"); + art::Ptr wvfPtr(opDetWvfHandle, i); + OpHitWvfVector.push_back(wvfPtr); + OpHitWvfValid.push_back(true); + found = true; + break; + } + } + if (!found) + { + ProducerUtils::PrintInColor("No matching OpDetWaveform found for OpHit channel " + ProducerUtils::str(int(opChannel)), ProducerUtils::GetColor("red"), "Debug"); + OpHitWvfVector.push_back(art::Ptr()); // Fill with a null pointer to keep the indices consistent + OpHitWvfValid.push_back(false); + } + } + return; + } + + + void AdjOpHitsUtils::GetOpHitWaveforms(const std::vector> &OpHitVector, std::vector> &OpHitWvfVector, std::vector &OpHitWvfValid, art::Event const &evt) + { // Repeat the function but for waveforms of type std::vector + // Get the OpWaveform handle + art::Handle> opWvfHandle; + evt.getByLabel(fOpWaveformLabel, opWvfHandle); + auto const clockData = art::ServiceHandle()->DataFor(evt); + if (!opWvfHandle.isValid()) + { + ProducerUtils::PrintInColor("Invalid OpWaveform handle", ProducerUtils::GetColor("red"), "Error"); + return; + } + + // Loop over the OpHit vector and get the corresponding OpWaveform + for (const auto &hit : OpHitVector) + { + bool found = false; + float hitTime = -1e6; + unsigned int opChannel = hit->OpChannel(); + auto TickPeriod = clockData.OpticalClock().TickPeriod(); + + if (fOpHitTimeVariable == "StartTime") + hitTime = hit->StartTime(); + else // Default to PeakTime + hitTime = hit->PeakTime(); + if (fOpHitTime2us) { + hitTime *= TickPeriod; + } + + for (size_t i = 0; i < opWvfHandle->size(); i++) + { // Match by channel number and amplitude + if ((*opWvfHandle)[i].Channel() == opChannel && hitTime >= (*opWvfHandle)[i].TimeStamp() && hitTime <= (*opWvfHandle)[i].TimeStamp() + (*opWvfHandle)[i].Signal().size() * TickPeriod) + { + ProducerUtils::PrintInColor("Matching OpWaveform found for OpHit channel " + ProducerUtils::str(int(opChannel)), ProducerUtils::GetColor("green"), "Debug"); + art::Ptr wvfPtr(opWvfHandle, i); + OpHitWvfVector.push_back(wvfPtr); + OpHitWvfValid.push_back(true); + found = true; + break; + } + } + + if (!found) { + ProducerUtils::PrintInColor("No matching OpWaveform found for OpHit channel " + ProducerUtils::str(int(opChannel)), ProducerUtils::GetColor("red"), "Debug"); + OpHitWvfVector.push_back(art::Ptr()); + OpHitWvfValid.push_back(false); + } + } + return; + } + + + void AdjOpHitsUtils::GetOpHitSignal(const std::vector> &OpHitVector, std::vector> &OpHitWvfIntVector, std::vector &OpHitWvfTime, std::vector &OpHitWvfValid, art::Event const &evt) + { + auto geoName = geom->DetectorName(); + if (geoName.find("dune10kt") != std::string::npos) { + std::vector> OpHitWvfVector; + GetOpHitWaveforms(OpHitVector, OpHitWvfVector, OpHitWvfValid, evt); + // Get entries in OpHitWvfvector and save to OpHitWvfIntVector as integers + for (int i = 0; i < int(OpHitWvfVector.size()); i++) + { + art::Ptr wvf = OpHitWvfVector[i]; + if (OpHitWvfValid[i]) { + std::vector wvfInt; + for (auto adc : wvf->Signal()) { + wvfInt.push_back(int(adc)); + } + OpHitWvfIntVector.push_back(wvfInt); + OpHitWvfTime.push_back(wvf->TimeStamp()); + } + else { + OpHitWvfIntVector.push_back(std::vector{}); + OpHitWvfTime.push_back(-1e6); + } + } + } + else if (geoName.find("dunevd10kt") != std::string::npos) { + std::vector> OpHitWvfVector; + GetOpHitWaveforms(OpHitVector, OpHitWvfVector, OpHitWvfValid, evt); + for (int i = 0; i < int(OpHitWvfVector.size()); i++) + { + art::Ptr wvf = OpHitWvfVector[i]; + if (OpHitWvfValid[i]) { + std::vector wvfInt; + for (auto adc : wvf->Waveform()) { + wvfInt.push_back(int(adc)); + } + OpHitWvfIntVector.push_back(wvfInt); + OpHitWvfTime.push_back(wvf->TimeStamp()); + } + else { + OpHitWvfIntVector.push_back(std::vector{}); + OpHitWvfTime.push_back(-1e6); + } + } + } + } } // namespace solar \ No newline at end of file diff --git a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h index eda7f8d9..0739e73d 100644 --- a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h +++ b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h @@ -24,6 +24,8 @@ #include "lardata/Utilities/AssociationUtil.h" #include "larcorealg/Geometry/GeometryCore.h" #include "larcorealg/Geometry/OpDetGeo.h" +#include "lardataobj/RawData/OpDetWaveform.h" +#include "lardataobj/RecoBase/OpWaveform.h" #include "lardataobj/RecoBase/OpHit.h" #include "lardataobj/RecoBase/OpFlash.h" #include "messagefacility/MessageLogger/MessageLogger.h" @@ -44,10 +46,12 @@ namespace solar { int Plane; int NHit; - double Time; + double MainOpHitPE; + double MainOpHitTime; + double MainOpHitAmplitude; double TimeWidth; + double TimeWeighted; double PE; - double MaxPE; std::vector PEperOpDet; double FastToTotal; double X; @@ -57,6 +61,9 @@ namespace solar double YWidth; double ZWidth; double STD; + std::vector MainOpWaveform; + float MainOpWaveformTime; + bool MainOpWaveformValid; }; struct OpFlashes { @@ -78,17 +85,22 @@ namespace solar int GetOpHitPlane(const art::Ptr &hit, float buffer); std::map GetOpHitPlaneMap(const std::vector> &OpHitVector); bool CheckOpHitPlane(std::map OpHitPlane, int refHit1, int refHit2); - // void CalcCentroid(std::vector> Hits, double x, double y, double z); - // double GaussianPDF(double x, double mean, double sigma); + void GetOpHitWaveforms(const std::vector> &OpHitVector, std::vector> &OpDetWvfVector, std::vector &OpHitWvfValid, art::Event const &evt); + void GetOpHitWaveforms(const std::vector> &OpHitVector, std::vector> &OpDetWvfVector, std::vector &OpHitWvfValid, art::Event const &evt); + void GetOpHitSignal(const std::vector> &OpHitVector, std::vector> &OpDetWvfIntVector, std::vector &OpHitWvfTime, std::vector &OpHitWvfValid, art::Event const &evt); private: art::ServiceHandle geom; geo::WireReadoutGeom const& wireReadout = art::ServiceHandle()->Get(); // From fhicl configuration + const std::string fOpWaveformLabel; + const std::string fOpHitLabel; const std::string fOpHitTimeVariable; + const bool fOpHitTime2us; const int fOpFlashAlgoNHit; const float fOpFlashAlgoMinTime; const float fOpFlashAlgoMaxTime; + const float fOpFlashAlgoWeightedTime; const float fOpFlashAlgoRad; const float fOpFlashAlgoPE; const float fOpFlashAlgoTriggerPE; diff --git a/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc b/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc index edcdf081..3f6e58ef 100644 --- a/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc +++ b/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc @@ -64,11 +64,13 @@ namespace solar private: // The parameters we'll read from the .fcl file. art::ServiceHandle geo; + std::string fOpWaveformLabel; // Input tag for OpDetWaveform collection std::string fOpHitLabel; // Input tag for OpHit collection std::string fOpHitTimeVariable; int fOpFlashAlgoNHit; float fOpFlashAlgoMinTime; float fOpFlashAlgoMaxTime; + bool fOpFlashAlgoWeightedTime; float fOpFlashAlgoRad; float fOpFlashAlgoPE; float fOpFlashAlgoTriggerPE; @@ -96,8 +98,8 @@ namespace solar fOpHitLabel(p.get("OpHitLabel", "ophitspe")), fOpHitTimeVariable(p.get("OpHitTimeVariable", "PeakTime")), fOpFlashAlgoNHit(p.get("OpFlashAlgoNHit", 3)), - fOpFlashAlgoMinTime(p.get("OpFlashAlgoMinTime", 0.010)), // 10 ns [0.6 tick] - fOpFlashAlgoMaxTime(p.get("OpFlashAlgoMaxTime", 0.016)), // 16 ns [1 tick] + fOpFlashAlgoMinTime(p.get("OpFlashAlgoMinTime", 0.32)), // [20 PDS tick] + fOpFlashAlgoMaxTime(p.get("OpFlashAlgoMaxTime", 0.96)), // [60 PDS tick] fOpFlashAlgoRad(p.get("OpFlashAlgoRad", 500.0)), fOpFlashAlgoPE(p.get("OpFlashAlgoPE", 1.5)), fOpFlashAlgoTriggerPE(p.get("OpFlashAlgoTriggerPE", 1.5)), @@ -205,7 +207,13 @@ namespace solar double OpFlashdY = OpFlash.YWidth; double OpFlashZ = OpFlash.Z; double OpFlashdZ = OpFlash.ZWidth; - double OpFlashT = OpFlash.Time - fOpFlashTimeOffset; // Convert to time to us happens in MakeFlashVector + double OpFlashT = -1e6; + if (fOpFlashAlgoWeightedTime) { + OpFlashT = OpFlash.TimeWeighted - fOpFlashTimeOffset; // Convert to time to us happens in MakeFlashVector + } + else { + OpFlashT = OpFlash.MainOpHitTime - fOpFlashTimeOffset; // Convert to time to us happens in MakeFlashVector + } double OpFlashdT = OpFlash.TimeWidth; // Convert to time to us happens in MakeFlashVector std::vector OpFlashPEs = OpFlash.PEperOpDet; @@ -234,7 +242,7 @@ namespace solar if (InBeamFrame) OnBeamTime = 1; - oflashes.emplace_back(OpFlashT, OpFlashdT, OpFlashT, Plane, + oflashes.emplace_back(OpFlashT, OpFlashdT, OpFlashT, Plane, // Using Frame to save the OpFlash Plane OpFlashPEs, OnPlane, OnBeamTime, OpFlashFast2Tot, OpFlashX, OpFlashdX, OpFlashY, OpFlashdY, OpFlashZ, OpFlashdZ); } diff --git a/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl b/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl index fba0355e..2f586eb4 100644 --- a/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl +++ b/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl @@ -3,10 +3,12 @@ BEGIN_PROLOG solar_opflash_dune10kt_1x2x6: { - OpHitLabel: "ophitspe" # The label for the process which ran the OpHit + OpWaveformLabel: "opdec" # The label for the process which ran the OpDigi + OpHitLabel: "ophitspe" # The label for the process which ran the OpHit module_type: "SolarOpFlash" - OpHitTimeVariable: "PeakTime" # Time variable to use for OpHit selection: "PeakTime" or "StartTime". - OpFlashTimeOffset: 18.1 # Time offset to be applied to the OpFlash time in [us] units. + OpHitTimeVariable: "PeakTime" # Time variable to use for OpHit selection: "PeakTime" or "StartTime". + OpHitTime2us: false # If true, the OpHit time is in [ticks] and will be converted to [us]. + OpFlashTimeOffset: 18.1 # Time offset to be applied to the OpFlash time in [us] units. XACathodeX: 0 # X position of the VD cathode XAs in [cm]. XAMembraneY: 0 # Y position of the VD membrane XAs in [cm]. @@ -14,8 +16,8 @@ solar_opflash_dune10kt_1x2x6: XAStartCapZ: 0 # Z position of the VD start cap XAs in [cm]. OpFlashAlgoNHit: 3 # Min number of hits to consider a flash. Change to > 3 for bkg productions. - OpFlashAlgoMinTime: 0.005 # Negative time window to look for adj. OpHits in [us] units. - OpFlashAlgoMaxTime: 0.016 # Positive time window to look for adj. OpHits in [us] units. + OpFlashAlgoMinTime: 0.32 # Negative time window to look for adj. OpHits in [us] units. + OpFlashAlgoMaxTime: 0.96 # Positive time window to look for adj. OpHits in [us] units. OpFlashAlgoRad: 300 # Distance to look for adj. OpHits in [cm] units. OpFlashAlgoPE: 1.5 # PE threshold to look for adj. OpHits. OpFlashAlgoTriggerPE: 1.5 # PE threshold to trigger an OpFlash. @@ -26,6 +28,7 @@ solar_opflash_dune10kt_1x2x6: } solar_opflash_dunevd10kt_1x8x14_3view_30deg: @local::solar_opflash_dune10kt_1x2x6 +solar_opflash_dunevd10kt_1x8x14_3view_30deg.OpWaveformLabel: "opdigi10ppm" solar_opflash_dunevd10kt_1x8x14_3view_30deg.OpHitLabel: "ophit10ppm" solar_opflash_dunevd10kt_1x8x14_3view_30deg.OpHitTimeVariable: "StartTime" solar_opflash_dunevd10kt_1x8x14_3view_30deg.OpFlashTimeOffset: 0