From 8a715545cabc016ebb8f88eff37217131ae78cdf Mon Sep 17 00:00:00 2001
From: mantheys <sergio.manthey@ciemat.es>
Date: Tue, 18 Nov 2025 12:43:04 +0100
Subject: [PATCH 1/3] ADD OPHIT TO WAVEFORM ASSOCIATION FUNCTIONS AND ADD MAIN
 WAVEFORM INFO TO OPFLASH

---
 duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc      | 212 ++++++++++++++++--
 duneopdet/LowEPDSUtils/AdjOpHitsUtils.h       |  12 +-
 duneopdet/LowEPDSUtils/SolarOpFlash_module.cc |  14 +-
 duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl   |   2 +
 4 files changed, 212 insertions(+), 28 deletions(-)

diff --git a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc
index 1515339f..c6824909 100644
--- a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc
+++ b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc
@@ -5,10 +5,13 @@ using namespace producer;
 namespace solar
 {
   AdjOpHitsUtils::AdjOpHitsUtils(fhicl::ParameterSet const &p)
-      : fOpHitTimeVariable(p.get<std::string>("OpHitTimeVariable", "PeakTime")), // Variable to use for time sorting ("StartTime" or "PeakTime")
+      : fOpWaveformLabel(p.get<std::string>("OpWaveformLabel", "opdetwaveform")), // Label for OpDetWaveform collection
+        fOpHitLabel(p.get<std::string>("OpHitLabel", "ophit")), // Label for OpHit collection
+        fOpHitTimeVariable(p.get<std::string>("OpHitTimeVariable", "PeakTime")), // Variable to use for time sorting ("StartTime" or "PeakTime")
         fOpFlashAlgoNHit(p.get<int>("OpFlashAlgoNHit", 3)),       // Minimum number of OpHits in a cluster to consider it for flash creation.
-        fOpFlashAlgoMinTime(p.get<float>("OpFlashAlgoMinTime", 0.010)), // Negative time window to look for adj. OpHits. Default for HD 10 ns [0.6 tick]
+        fOpFlashAlgoMinTime(p.get<float>("OpFlashAlgoMinTime", 0.08)), // Negative time window to look for adj. OpHits. Default for HD 8 ns [0.5 tick]
         fOpFlashAlgoMaxTime(p.get<float>("OpFlashAlgoMaxTime", 0.016)), // Positive time window to look for adj. OpHits. Default for HD 16 ns [1 tick]
+        fOpFlashAlgoWeightedTime(p.get<bool>("OpFlashAlgoWeightedTime", false)), // Whether to use weighted time of trigger time for flash time reference.
         fOpFlashAlgoRad(p.get<float>("OpFlashAlgoRad", 300.0)),     // Radius to look for adj. OpHits in [cm] units.
         fOpFlashAlgoPE(p.get<float>("OpFlashAlgoPE", 1.5)),         // Minimum PE of OpHit to consider it for flash creation.
         fOpFlashAlgoTriggerPE(p.get<float>("OpFlashAlgoTriggerPE", 1.5)), // Minimum PE of OpHit to consider it as a trigger for flash creation.
@@ -20,15 +23,17 @@ namespace solar
         fXAStartCapZ(p.get<float>("XAStartCapZ", -96.5))
   {
   }
+  
+  
   void AdjOpHitsUtils::MakeFlashVector(std::vector<FlashInfo> &FlashVec, std::vector<std::vector<art::Ptr<recob::OpHit>>> &OpHitClusters, art::Event const &evt)
   {
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->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<art::Ptr<recob::OpHit>> Cluster : OpHitClusters)
     {
       if (!Cluster.empty())
@@ -40,10 +45,14 @@ namespace solar
           std::stable_sort(Cluster.begin(), Cluster.end(), [](art::Ptr<recob::OpHit> a, art::Ptr<recob::OpHit> b)
                           { return a->PeakTime() < b->PeakTime(); });
       }
+
       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 PE = 0;
       double MaxPE = 0;
@@ -59,30 +68,52 @@ namespace solar
       double YSum = 0;
       double ZSum = 0;
       double STD = 0;
+      std::vector<int> MainOpWaveform = {}; // Make vector for main OpWaveform with 1000 entries (max waveform size)
+
+      std::vector<bool> OpHitWvfValid = {};
+      std::vector<std::vector<int>> OpHitWvfIntVector = {};
+      GetOpHitSignal(Cluster, OpHitWvfIntVector, OpHitWvfValid, evt); // Get OpHit waveforms
 
       // Compute total number of PE and MaxPE.
       for (art::Ptr<recob::OpHit> PDSHit : Cluster)
       {
         NHit++;
+        double thisTime = -1e6;
+        double thisPE = PDSHit->PE();
         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() * TickPeriod; // Use StartTime
+        }
+        else {
+          thisTime = PDSHit->PeakTime() * TickPeriod; // Default to PeakTime
+        }
+        
+        PE += thisPE;
+        TimeSum += thisTime * thisPE;
+        
+        if (thisPE > MaxPE) {
+          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];
       }
 
       float HotPE = 0;
-      Time = TimeSum / PE;
+      TimeWeighted = TimeSum / PE;
+
       // Compute flash center from weighted average of "hottest" ophits.
       for (art::Ptr<recob::OpHit> PDSHit : Cluster)
       {
@@ -106,9 +137,9 @@ namespace solar
       {
         auto OpHitXYZ = wireReadout.OpDetGeoFromOpChannel(PDSHit->OpChannel()).GetCenter();
         if (fOpHitTimeVariable == "StartTime")
-          TimeWidth += (PDSHit->StartTime() * TickPeriod - Time) * (PDSHit->StartTime() * TickPeriod - Time);
+          TimeWidth += (TimeMax - PDSHit->StartTime() * TickPeriod) * (TimeMax - PDSHit->StartTime() * TickPeriod);
         else // Default to PeakTime
-          TimeWidth += (PDSHit->PeakTime() * TickPeriod - Time) * (PDSHit->PeakTime() * TickPeriod - Time);
+          TimeWidth += (TimeMax - PDSHit->PeakTime() * TickPeriod) * (TimeMax - PDSHit->PeakTime() * TickPeriod);
         
         XWidth += (OpHitXYZ.X() - X) * (OpHitXYZ.X() - X);
         YWidth += (OpHitXYZ.Y() - Y) * (OpHitXYZ.Y() - Y);
@@ -128,21 +159,23 @@ namespace solar
       // Compute FastToTotal according to the #PEs arriving within the first 10% of the time window wrt the total #PEs
       for (art::Ptr<recob::OpHit> PDSHit : Cluster)
       {
-        auto thisTime = 0.0;
+        auto thisTime = -1e6;
         if (fOpHitTimeVariable == "StartTime")
           thisTime = PDSHit->StartTime() * TickPeriod;
         else // Default to PeakTime
           thisTime = PDSHit->PeakTime() * 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, TimeMax, TimeWidth, TimeWeighted, PE, MaxPE, PEperOpDet, FastToTotal, X, Y, Z, XWidth, YWidth, ZWidth, STD, MainOpWaveform});
     }
     return;
   }
 
+
   float AdjOpHitsUtils::GetOpFlashPlaneSTD(const int Plane, const std::vector<float> varXY, const std::vector<float> varYZ, const std::vector<float> varXZ)
   {
     std::vector<float> var;
@@ -185,6 +218,7 @@ namespace solar
     return varstd;
   }
 
+
   void AdjOpHitsUtils::CalcAdjOpHits(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<std::vector<art::Ptr<recob::OpHit>>> &OpHitClusters, std::vector<std::vector<int>> &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.
@@ -390,6 +424,7 @@ namespace solar
     return;
   }
 
+
   void AdjOpHitsUtils::FlashMatchResidual(float &Residual, std::vector<art::Ptr<recob::OpHit>> Hits, double x, double y, double z)
   {
     if (Hits.size() == 0)
@@ -451,6 +486,7 @@ namespace solar
     return;
   }
 
+
   int AdjOpHitsUtils::GetOpHitPlane(const art::Ptr<recob::OpHit> &hit, float buffer)
   {
     std::string geoName = geom->DetectorName();
@@ -483,9 +519,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<int, int> AdjOpHitsUtils::GetOpHitPlaneMap(const std::vector<art::Ptr<recob::OpHit>> &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<int, int> OpHitPlane;
     for (const auto &hit : OpHitVector)
     {
@@ -494,9 +530,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<int, int> 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 +545,134 @@ namespace solar
       return false;
   }
 
+
+  void AdjOpHitsUtils::GetOpHitWaveforms(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<art::Ptr<raw::OpDetWaveform>> &OpHitWvfVector, std::vector<bool> &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<std::vector<raw::OpDetWaveform>> opDetWvfHandle;
+    evt.getByLabel(fOpWaveformLabel, opDetWvfHandle);
+    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->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<raw::OpDetWaveform>());  // 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;
+      unsigned int opChannel = hit->OpChannel();
+      float hitTime = hit->PeakTime() * clockData.OpticalClock().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() * clockData.OpticalClock().TickPeriod())
+        {
+          ProducerUtils::PrintInColor("Matching OpDetWaveform found for OpHit channel " + ProducerUtils::str(int(opChannel)), ProducerUtils::GetColor("green"), "Debug");
+          art::Ptr<raw::OpDetWaveform> 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<raw::OpDetWaveform>());  // Fill with a null pointer to keep the indices consistent
+        OpHitWvfValid.push_back(false);
+      }
+    }
+    return;
+  }
+
+
+  void AdjOpHitsUtils::GetOpHitWaveforms(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<art::Ptr<recob::OpWaveform>> &OpHitWvfVector, std::vector<bool> &OpHitWvfValid, art::Event const &evt)
+  { // Repeat the function but for waveforms of type std::vector<recob::OpWaveform>
+    // Get the OpWaveform handle
+    art::Handle<std::vector<recob::OpWaveform>> opWvfHandle;
+    evt.getByLabel(fOpWaveformLabel, opWvfHandle);
+    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->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;
+      unsigned int opChannel = hit->OpChannel();
+      float hitTime = hit->PeakTime() * clockData.OpticalClock().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() * clockData.OpticalClock().TickPeriod())
+        {
+          ProducerUtils::PrintInColor("Matching OpWaveform found for OpHit channel " + ProducerUtils::str(int(opChannel)), ProducerUtils::GetColor("green"), "Debug");
+          art::Ptr<recob::OpWaveform> 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<recob::OpWaveform>());
+        OpHitWvfValid.push_back(false);
+      }
+    }
+    return;
+  }
+
+
+  void AdjOpHitsUtils::GetOpHitSignal(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<std::vector<int>> &OpHitWvfIntVector, std::vector<bool> &OpHitWvfValid, art::Event const &evt)
+  {
+    auto geoName = geom->DetectorName();
+    if (geoName.find("dune10kt") != std::string::npos) {
+      std::vector<art::Ptr<recob::OpWaveform>> 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<recob::OpWaveform> wvf = OpHitWvfVector[i];
+        if (OpHitWvfValid[i]) {
+          std::vector<int> wvfInt;
+          for (auto adc : wvf->Signal()) {
+            wvfInt.push_back(int(adc));
+          }
+          OpHitWvfIntVector.push_back(wvfInt);
+        }
+        else {
+          OpHitWvfIntVector.push_back(std::vector<int>{});
+        }
+      } 
+    }
+    else if (geoName.find("dunevd10kt") != std::string::npos) {
+      std::vector<art::Ptr<raw::OpDetWaveform>> OpHitWvfVector;
+      GetOpHitWaveforms(OpHitVector, OpHitWvfVector, OpHitWvfValid, evt);
+      for (int i = 0; i < int(OpHitWvfVector.size()); i++)
+      {
+        art::Ptr<raw::OpDetWaveform> wvf = OpHitWvfVector[i];
+        if (OpHitWvfValid[i]) {
+          std::vector<int> wvfInt;
+          for (auto adc : wvf->Waveform()) {
+            wvfInt.push_back(int(adc));
+          }
+          OpHitWvfIntVector.push_back(wvfInt);
+        }
+        else {
+          OpHitWvfIntVector.push_back(std::vector<int>{});
+        }
+      }
+    }
+  }
+
 } // namespace solar
\ No newline at end of file
diff --git a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h
index eda7f8d9..d35630a8 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"
@@ -46,6 +48,7 @@ namespace solar
             int NHit;
             double Time;
             double TimeWidth;
+            double TimeWeighted;
             double PE;
             double MaxPE;
             std::vector<double> PEperOpDet;
@@ -57,6 +60,7 @@ namespace solar
             double YWidth;
             double ZWidth;
             double STD;
+            std::vector<int> MainOpWaveform;
         };
         struct OpFlashes
         {
@@ -78,17 +82,21 @@ namespace solar
         int GetOpHitPlane(const art::Ptr<recob::OpHit> &hit, float buffer);
         std::map<int, int> GetOpHitPlaneMap(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector);
         bool CheckOpHitPlane(std::map<int, int> OpHitPlane, int refHit1, int refHit2);
-        // void CalcCentroid(std::vector<art::Ptr<recob::OpHit>> Hits, double x, double y, double z);
-        // double GaussianPDF(double x, double mean, double sigma);
+        void GetOpHitWaveforms(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<art::Ptr<raw::OpDetWaveform>> &OpDetWvfVector, std::vector<bool> &OpHitWvfValid, art::Event const &evt);
+        void GetOpHitWaveforms(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<art::Ptr<recob::OpWaveform>> &OpDetWvfVector, std::vector<bool> &OpHitWvfValid, art::Event const &evt);
+        void GetOpHitSignal(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<std::vector<int>> &OpDetWvfIntVector, std::vector<bool> &OpHitWvfValid, art::Event const &evt);
 
     private:
         art::ServiceHandle<geo::Geometry> geom;
         geo::WireReadoutGeom const& wireReadout = art::ServiceHandle<geo::WireReadout>()->Get();
         // From fhicl configuration
+        const std::string fOpWaveformLabel;
+        const std::string fOpHitLabel;
         const std::string fOpHitTimeVariable;
         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..b294150b 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::Geometry> 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,7 +98,7 @@ namespace solar
         fOpHitLabel(p.get<std::string>("OpHitLabel", "ophitspe")),
         fOpHitTimeVariable(p.get<std::string>("OpHitTimeVariable", "PeakTime")),
         fOpFlashAlgoNHit(p.get<int>("OpFlashAlgoNHit", 3)),
-        fOpFlashAlgoMinTime(p.get<float>("OpFlashAlgoMinTime", 0.010)), // 10 ns [0.6 tick]
+        fOpFlashAlgoMinTime(p.get<float>("OpFlashAlgoMinTime", 0.008)), // 8 ns [0.5 tick]
         fOpFlashAlgoMaxTime(p.get<float>("OpFlashAlgoMaxTime", 0.016)), // 16 ns [1 tick]
         fOpFlashAlgoRad(p.get<float>("OpFlashAlgoRad", 500.0)),
         fOpFlashAlgoPE(p.get<float>("OpFlashAlgoPE", 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.Time - fOpFlashTimeOffset; // Convert to time to us happens in MakeFlashVector
+      }
       double OpFlashdT = OpFlash.TimeWidth; // Convert to time to us happens in MakeFlashVector
       std::vector<double> 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..8e19d6c9 100644
--- a/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl
+++ b/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl
@@ -3,6 +3,7 @@ BEGIN_PROLOG
 
 solar_opflash_dune10kt_1x2x6:
 {
+    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".
@@ -26,6 +27,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

From 9019af6b8381132e33ad2068dc1b7bb173f0a31f Mon Sep 17 00:00:00 2001
From: mantheys <sergio.manthey@ciemat.es>
Date: Wed, 19 Nov 2025 06:06:13 -0600
Subject: [PATCH 2/3] ADD NEW FLASH VARIABLES WITH WAVEFORM AND OPHIT
 BACKTRACKING

---
 duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc      | 59 +++++++++++++++----
 duneopdet/LowEPDSUtils/AdjOpHitsUtils.h       |  9 ++-
 duneopdet/LowEPDSUtils/SolarOpFlash_module.cc |  2 +-
 3 files changed, 55 insertions(+), 15 deletions(-)

diff --git a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc
index c6824909..baad1d81 100644
--- a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc
+++ b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc
@@ -36,15 +36,18 @@ namespace solar
     
     for (std::vector<art::Ptr<recob::OpHit>> Cluster : OpHitClusters)
     {
-      if (!Cluster.empty())
+      if (Cluster.empty())
       {
-        if (fOpHitTimeVariable == "StartTime")
-          std::stable_sort(Cluster.begin(), Cluster.end(), [](art::Ptr<recob::OpHit> a, art::Ptr<recob::OpHit> b)
-                           { return a->StartTime() < b->StartTime(); });
-        else // Default to PeakTime
-          std::stable_sort(Cluster.begin(), Cluster.end(), [](art::Ptr<recob::OpHit> a, art::Ptr<recob::OpHit> b)
-                          { return a->PeakTime() < b->PeakTime(); });
+        continue;
       }
+      
+      // Sorting already done in CalcAdjOpHits
+      // if (fOpHitTimeVariable == "StartTime")
+      //   std::stable_sort(Cluster.begin(), Cluster.end(), [](art::Ptr<recob::OpHit> a, art::Ptr<recob::OpHit> b)
+      //                    { return a->StartTime() < b->StartTime(); });
+      // else // Default to PeakTime
+      //   std::stable_sort(Cluster.begin(), Cluster.end(), [](art::Ptr<recob::OpHit> a, art::Ptr<recob::OpHit> b)
+      //                   { return a->PeakTime() < b->PeakTime(); });
 
       int Plane = GetOpHitPlane(Cluster[0], 0.1);
       int MaxIdx = 0;
@@ -54,9 +57,10 @@ namespace solar
       double TimeWidth = 0;
       double TimeWeighted = -1e6;
       double TimeSum = 0;
+      double Amplitude = 0;
       double PE = 0;
       double MaxPE = 0;
-      std::vector<double> PEperOpDet;
+      std::vector<double> PEperOpDet = {};
       double FastToTotal = 1;
       double X = 0;
       double Y = 0;
@@ -69,10 +73,13 @@ namespace solar
       double ZSum = 0;
       double STD = 0;
       std::vector<int> MainOpWaveform = {}; // Make vector for main OpWaveform with 1000 entries (max waveform size)
+      bool MainOpWaveformValid = false;
+      float MainOpWaveformTime = -1e6;
 
       std::vector<bool> OpHitWvfValid = {};
+      std::vector<float> OpHitWvfTime = {};
       std::vector<std::vector<int>> OpHitWvfIntVector = {};
-      GetOpHitSignal(Cluster, OpHitWvfIntVector, OpHitWvfValid, evt); // Get OpHit waveforms
+      GetOpHitSignal(Cluster, OpHitWvfIntVector, OpHitWvfTime, OpHitWvfValid, evt); // Get OpHit waveforms
 
       // Compute total number of PE and MaxPE.
       for (art::Ptr<recob::OpHit> PDSHit : Cluster)
@@ -80,6 +87,7 @@ namespace solar
         NHit++;
         double thisTime = -1e6;
         double thisPE = PDSHit->PE();
+        double thisAmp = PDSHit->Amplitude();
         auto thisPlane = GetOpHitPlane(PDSHit, 0.1);
         
         if (thisPlane != Plane) {
@@ -98,6 +106,7 @@ namespace solar
         TimeSum += thisTime * thisPE;
         
         if (thisPE > MaxPE) {
+          Amplitude = thisAmp;
           MaxPE = thisPE;
           MaxIdx = Idx;
           TimeMax = thisTime;
@@ -109,6 +118,8 @@ namespace solar
 
       if (OpHitWvfValid[MaxIdx]) {
         MainOpWaveform = OpHitWvfIntVector[MaxIdx];
+        MainOpWaveformTime = OpHitWvfTime[MaxIdx];
+        MainOpWaveformValid = true;
       }
 
       float HotPE = 0;
@@ -170,7 +181,29 @@ namespace solar
           FastToTotal += PDSHit->PE();
       }
       FastToTotal /= PE;
-      FlashVec.push_back(FlashInfo{Plane, NHit, TimeMax, TimeWidth, TimeWeighted, PE, MaxPE, PEperOpDet, FastToTotal, X, Y, Z, XWidth, YWidth, ZWidth, STD, MainOpWaveform});
+      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;
   }
@@ -633,7 +666,7 @@ namespace solar
   }
 
 
-  void AdjOpHitsUtils::GetOpHitSignal(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<std::vector<int>> &OpHitWvfIntVector, std::vector<bool> &OpHitWvfValid, art::Event const &evt)
+  void AdjOpHitsUtils::GetOpHitSignal(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<std::vector<int>> &OpHitWvfIntVector, std::vector<float> &OpHitWvfTime, std::vector<bool> &OpHitWvfValid, art::Event const &evt)
   {
     auto geoName = geom->DetectorName();
     if (geoName.find("dune10kt") != std::string::npos) {
@@ -649,9 +682,11 @@ namespace solar
             wvfInt.push_back(int(adc));
           }
           OpHitWvfIntVector.push_back(wvfInt);
+          OpHitWvfTime.push_back(wvf->TimeStamp());
         }
         else {
           OpHitWvfIntVector.push_back(std::vector<int>{});
+          OpHitWvfTime.push_back(-1e6);
         }
       } 
     }
@@ -667,9 +702,11 @@ namespace solar
             wvfInt.push_back(int(adc));
           }
           OpHitWvfIntVector.push_back(wvfInt);
+          OpHitWvfTime.push_back(wvf->TimeStamp());
         }
         else {
           OpHitWvfIntVector.push_back(std::vector<int>{});
+          OpHitWvfTime.push_back(-1e6);
         }
       }
     }
diff --git a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h
index d35630a8..e56fe8dc 100644
--- a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h
+++ b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h
@@ -46,11 +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<double> PEperOpDet;
             double FastToTotal;
             double X;
@@ -61,6 +62,8 @@ namespace solar
             double ZWidth;
             double STD;
             std::vector<int> MainOpWaveform;
+            float MainOpWaveformTime;
+            bool MainOpWaveformValid;
         };
         struct OpFlashes
         {
@@ -84,7 +87,7 @@ namespace solar
         bool CheckOpHitPlane(std::map<int, int> OpHitPlane, int refHit1, int refHit2);
         void GetOpHitWaveforms(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<art::Ptr<raw::OpDetWaveform>> &OpDetWvfVector, std::vector<bool> &OpHitWvfValid, art::Event const &evt);
         void GetOpHitWaveforms(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<art::Ptr<recob::OpWaveform>> &OpDetWvfVector, std::vector<bool> &OpHitWvfValid, art::Event const &evt);
-        void GetOpHitSignal(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<std::vector<int>> &OpDetWvfIntVector, std::vector<bool> &OpHitWvfValid, art::Event const &evt);
+        void GetOpHitSignal(const std::vector<art::Ptr<recob::OpHit>> &OpHitVector, std::vector<std::vector<int>> &OpDetWvfIntVector, std::vector<float> &OpHitWvfTime, std::vector<bool> &OpHitWvfValid, art::Event const &evt);
 
     private:
         art::ServiceHandle<geo::Geometry> geom;
diff --git a/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc b/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc
index b294150b..04d08535 100644
--- a/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc
+++ b/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc
@@ -212,7 +212,7 @@ namespace solar
         OpFlashT = OpFlash.TimeWeighted - fOpFlashTimeOffset; // Convert to time to us happens in MakeFlashVector
       }
       else {
-        OpFlashT = OpFlash.Time - fOpFlashTimeOffset; // Convert to time to us happens in MakeFlashVector
+        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<double> OpFlashPEs = OpFlash.PEperOpDet;

From 586e0d565100b931bab2f3ed7e9a8e4926d6a356 Mon Sep 17 00:00:00 2001
From: mantheys <sergio.manthey@ciemat.es>
Date: Fri, 21 Nov 2025 12:17:24 +0100
Subject: [PATCH 3/3] ADD OPTION TO CONVERT UNITS OF OPHITS IF THEY ARE IN
 TICKS

---
 duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc      | 128 +++++++++++++-----
 duneopdet/LowEPDSUtils/AdjOpHitsUtils.h       |   1 +
 duneopdet/LowEPDSUtils/SolarOpFlash_module.cc |   4 +-
 duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl   |  13 +-
 4 files changed, 101 insertions(+), 45 deletions(-)

diff --git a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc
index baad1d81..7a519418 100644
--- a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc
+++ b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.cc
@@ -5,12 +5,13 @@ using namespace producer;
 namespace solar
 {
   AdjOpHitsUtils::AdjOpHitsUtils(fhicl::ParameterSet const &p)
-      : fOpWaveformLabel(p.get<std::string>("OpWaveformLabel", "opdetwaveform")), // Label for OpDetWaveform collection
-        fOpHitLabel(p.get<std::string>("OpHitLabel", "ophit")), // Label for OpHit collection
-        fOpHitTimeVariable(p.get<std::string>("OpHitTimeVariable", "PeakTime")), // Variable to use for time sorting ("StartTime" or "PeakTime")
+      : fOpWaveformLabel(p.get<std::string>("OpWaveformLabel", "opdetwaveform")), // Label for OpDetWaveform collection.
+        fOpHitLabel(p.get<std::string>("OpHitLabel", "ophit")), // Label for OpHit collection.
+        fOpHitTimeVariable(p.get<std::string>("OpHitTimeVariable", "PeakTime")), // Variable to use for time sorting ("StartTime" or "PeakTime").
+        fOpHitTime2us(p.get<bool>("OpHitTime2us", false)), // Conversion factor from OpHit time units to microseconds. Default factor is TickPeriod() from DetectorClocksService.
         fOpFlashAlgoNHit(p.get<int>("OpFlashAlgoNHit", 3)),       // Minimum number of OpHits in a cluster to consider it for flash creation.
-        fOpFlashAlgoMinTime(p.get<float>("OpFlashAlgoMinTime", 0.08)), // Negative time window to look for adj. OpHits. Default for HD 8 ns [0.5 tick]
-        fOpFlashAlgoMaxTime(p.get<float>("OpFlashAlgoMaxTime", 0.016)), // Positive time window to look for adj. OpHits. Default for HD 16 ns [1 tick]
+        fOpFlashAlgoMinTime(p.get<float>("OpFlashAlgoMinTime", 0.32)), // Negative time window to look for adj. OpHits in [us].
+        fOpFlashAlgoMaxTime(p.get<float>("OpFlashAlgoMaxTime", 0.96)), // Positive time window to look for adj. OpHits in [us].
         fOpFlashAlgoWeightedTime(p.get<bool>("OpFlashAlgoWeightedTime", false)), // Whether to use weighted time of trigger time for flash time reference.
         fOpFlashAlgoRad(p.get<float>("OpFlashAlgoRad", 300.0)),     // Radius to look for adj. OpHits in [cm] units.
         fOpFlashAlgoPE(p.get<float>("OpFlashAlgoPE", 1.5)),         // Minimum PE of OpHit to consider it for flash creation.
@@ -40,14 +41,6 @@ namespace solar
       {
         continue;
       }
-      
-      // Sorting already done in CalcAdjOpHits
-      // if (fOpHitTimeVariable == "StartTime")
-      //   std::stable_sort(Cluster.begin(), Cluster.end(), [](art::Ptr<recob::OpHit> a, art::Ptr<recob::OpHit> b)
-      //                    { return a->StartTime() < b->StartTime(); });
-      // else // Default to PeakTime
-      //   std::stable_sort(Cluster.begin(), Cluster.end(), [](art::Ptr<recob::OpHit> a, art::Ptr<recob::OpHit> b)
-      //                   { return a->PeakTime() < b->PeakTime(); });
 
       int Plane = GetOpHitPlane(Cluster[0], 0.1);
       int MaxIdx = 0;
@@ -96,10 +89,13 @@ namespace solar
         }
 
         if (fOpHitTimeVariable == "StartTime") {
-          thisTime = PDSHit->StartTime() * TickPeriod; // Use StartTime
+          thisTime = PDSHit->StartTime(); // Use StartTime
         }
         else {
-          thisTime = PDSHit->PeakTime() * TickPeriod; // Default to PeakTime
+          thisTime = PDSHit->PeakTime(); // Default to PeakTime
+        }
+        if (fOpHitTime2us) {
+          thisTime *= TickPeriod; // Convert to microseconds
         }
         
         PE += thisPE;
@@ -136,6 +132,7 @@ namespace solar
           HotPE += PDSHit->PE();
         }
       }
+      
       X = XSum / HotPE;
       Y = YSum / HotPE;
       Z = ZSum / HotPE;
@@ -146,12 +143,18 @@ namespace solar
       std::vector<float> varXZ;
       for (art::Ptr<recob::OpHit> PDSHit : Cluster)
       {
+        float ThisOpHitTime = -1e6;
         auto OpHitXYZ = wireReadout.OpDetGeoFromOpChannel(PDSHit->OpChannel()).GetCenter();
+        
         if (fOpHitTimeVariable == "StartTime")
-          TimeWidth += (TimeMax - PDSHit->StartTime() * TickPeriod) * (TimeMax - PDSHit->StartTime() * TickPeriod);
+          ThisOpHitTime = PDSHit->StartTime();
         else // Default to PeakTime
-          TimeWidth += (TimeMax - PDSHit->PeakTime() * TickPeriod) * (TimeMax - PDSHit->PeakTime() * TickPeriod);
+          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);
@@ -172,9 +175,12 @@ namespace solar
       {
         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;
+        }
         
         // Check if thisTime is within 10% of the time window in positive and negative direction
         if (std::abs(thisTime - TimeMax) <= 0.05 * TimeWidth)
@@ -289,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();
@@ -322,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;
         }
         
@@ -338,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;
         }
 
@@ -348,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;
@@ -397,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)
@@ -411,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;
         }
 
@@ -421,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;
@@ -600,11 +635,21 @@ namespace solar
     for (const auto &hit : OpHitVector)
     {
       bool found = false;
+      float hitTime = -1e6;
       unsigned int opChannel = hit->OpChannel();
-      float hitTime = hit->PeakTime() * clockData.OpticalClock().TickPeriod();
+      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() * clockData.OpticalClock().TickPeriod())
+        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<raw::OpDetWaveform> wvfPtr(opDetWvfHandle, i);
@@ -641,11 +686,21 @@ namespace solar
     for (const auto &hit : OpHitVector)
     {
       bool found = false;
+      float hitTime = -1e6;
       unsigned int opChannel = hit->OpChannel();
-      float hitTime = hit->PeakTime() * clockData.OpticalClock().TickPeriod();
+      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() * clockData.OpticalClock().TickPeriod())
+        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<recob::OpWaveform> wvfPtr(opWvfHandle, i);
@@ -711,5 +766,4 @@ namespace solar
       }
     }
   }
-
 } // namespace solar
\ No newline at end of file
diff --git a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h
index e56fe8dc..0739e73d 100644
--- a/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h
+++ b/duneopdet/LowEPDSUtils/AdjOpHitsUtils.h
@@ -96,6 +96,7 @@ namespace solar
         const std::string fOpWaveformLabel;
         const std::string fOpHitLabel;
         const std::string fOpHitTimeVariable;
+        const bool fOpHitTime2us;
         const int fOpFlashAlgoNHit;
         const float fOpFlashAlgoMinTime;
         const float fOpFlashAlgoMaxTime;
diff --git a/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc b/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc
index 04d08535..3f6e58ef 100644
--- a/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc
+++ b/duneopdet/LowEPDSUtils/SolarOpFlash_module.cc
@@ -98,8 +98,8 @@ namespace solar
         fOpHitLabel(p.get<std::string>("OpHitLabel", "ophitspe")),
         fOpHitTimeVariable(p.get<std::string>("OpHitTimeVariable", "PeakTime")),
         fOpFlashAlgoNHit(p.get<int>("OpFlashAlgoNHit", 3)),
-        fOpFlashAlgoMinTime(p.get<float>("OpFlashAlgoMinTime", 0.008)), // 8 ns [0.5 tick]
-        fOpFlashAlgoMaxTime(p.get<float>("OpFlashAlgoMaxTime", 0.016)), // 16 ns [1 tick]
+        fOpFlashAlgoMinTime(p.get<double>("OpFlashAlgoMinTime", 0.32)), //  [20 PDS tick]
+        fOpFlashAlgoMaxTime(p.get<double>("OpFlashAlgoMaxTime", 0.96)), //  [60 PDS tick]
         fOpFlashAlgoRad(p.get<float>("OpFlashAlgoRad", 500.0)),
         fOpFlashAlgoPE(p.get<float>("OpFlashAlgoPE", 1.5)),
         fOpFlashAlgoTriggerPE(p.get<float>("OpFlashAlgoTriggerPE", 1.5)),
diff --git a/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl b/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl
index 8e19d6c9..2f586eb4 100644
--- a/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl
+++ b/duneopdet/LowEPDSUtils/fcl/SolarOpFlash.fcl
@@ -3,11 +3,12 @@ BEGIN_PROLOG
 
 solar_opflash_dune10kt_1x2x6:
 {
-    OpWaveformLabel:            "opdec"     # The label for the process which ran the OpDigi
-    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].
@@ -15,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.