Updating wvf pre-processing

duneopdet/OpticalDetector/WaveformPreProcessing.fcl

@@ -2,14 +2,19 @@ BEGIN_PROLOG
 
 protodunevd_wvfpreprocessing:
 {
-  module_type:          "WaveformPreProcessing"
-  InputModule:          "pdvddaphne:daq" // Input tag for OpDetWaveform collection (data)
-  ApplyDenoising:       true
-  IsPDVD:               true //Not to apply to PMTs
-  Lambda:               10.0
-  MaxTicks:             30000 //Maximum number of ticks in a waveform
-  SecondBaselineSub:    0.1 //the mode of lowest SecondBaslineSub of the signal for second baseline estimate
-  FirstBaselineSub:     900 //around 3 times the expected large signals
+  module_type:               "WaveformPreProcessing"
+  InputModule:               "pdvddaphne:daq" // Input tag for OpDetWaveform collection (data)
+  ApplyDenoising:            true
+  RemoveBaselineFluctuation: true
+  Lambda:                    7.0
+  MaxTicks:                  30000 //Maximum number of ticks in a waveform
+  MaxTicksDent:              1000 //Maximum number of ticks for "dent" formation
+  MaxTicksSat:               500  //Maximum number of saturated ticks to apply baseline fluctuation removal
+  MaxDentDepth:              0.8  //Maximum percentage of falling to be considered a laser failure
+  SecondBaselineSub:         0.1  //the mode of lowest SecondBaslineSub of the signal for second baseline estimate
+  FirstBaselineSub:          900  //around 3 times the expected large signals
+  DynamicRangeSaturation:    16383 # maximum ADC count in the waveform (=14 bits)
+  IgnoreChannels:            [2010, 2011, 2020, 2021, 2030, 2031, 2040, 2041, 3010, 3020, 3030, 3040, 3050, 3060, 3070, 3080, 3090, 3100, 3110, 3120, 3130, 3140, 3150, 3160, 3170, 3180, 3190, 3200, 3210, 3220, 3230, 3240] //All top membrane and PMTs
 }
 
 END_PROLOG

duneopdet/OpticalDetector/WaveformPreProcessing_module.cc

@@ -7,8 +7,10 @@
 // from cetpkgsupport v1_14_01.
 // This module works on data waveforms and can perform two tasks
 // 1. Remove baselines that fluctuate with time (code developed by A. Paudel)
+//    for waveforms which are not "too much" saturated
 // 2. Use a denoising algorithm (the same one used in PDSP) to make
-// waveforms smoother (fhicl adjustable)
+//    waveforms smoother 
+// (both fhicl adjustable)
 // This module produces a new set of OpDetWaveforms that should be
 // used as input to the OpHitFinder
 // TO DO: Implement a ROI finder to decrease the size of the optical
@@ -64,8 +66,10 @@ namespace opdet{
 
     // Required functions.
     void produce(art::Event & evt) override;
+    bool CheckSaturation(std::vector<double> wf, Float_t fDynamicRangeSaturation, size_t fMaxTicksSat);
     void BaselineExtractor(std::vector<double>& wf);
-    void DenoisingAlgo(std::vector<double>& waveform, double lambda);
+    void DentCorrection(std::vector<double> &wf, Float_t fDynamicRangeSaturation, size_t MaxTicksDent, int channel, double MaxDentDepth);
+    void DenoisingAlgo(std::vector<double>& waveform, double lambda, Float_t fDynamicRangeSaturation);
     bool TV1D_denoise(std::vector<double>& waveform, std::vector<double>& outwaveform, const double lambda);
     std::vector<double> ComputeMovingAverage(const std::vector<double>& data, int n);
     double median_func(std::vector<double> vec, double fraction);
@@ -81,11 +85,17 @@ namespace opdet{
     // The parameters we'll read from the .fcl file.
     std::string fInputModule; // Input tag for OpDetWaveform collection
     bool fApplyDenoising;
-    bool fIsPDVD;
+    bool fRemoveBaselineFluctuation;
+    std::vector<int> fIgnoreChannels;
+   // bool fIsPDVD; Causing warning errors in c14:prof
     double fLambda;
-    int fMaxTicks;
+    int fMaxTicks;        //Maximum number of ticks in waveform
+    size_t fMaxTicksDent;     //Maximum number of ticks for dent correction
+    size_t fMaxTicksSat;   //Maximum number of saturated ticks to apply baseline removal
+    Float_t fDynamicRangeSaturation;
     double fSecondBaselineSub;
     double fFirstBaselineSub;
+    double fMaxDentDepth; //Maximum percentage of falling from saturation to be considered a laser failure
   };
 }//namespace opdet
 
@@ -104,11 +114,16 @@ namespace opdet {
     // Indicate that the Input Module comes from .fcl
     fInputModule        = p.get<std::string>("InputModule"); 
     fApplyDenoising     = p.get<bool>("ApplyDenoising"); 
-    fIsPDVD             = p.get<bool>("IsPDVD"); 
+    fIgnoreChannels     = p.get<std::vector<int> >("IgnoreChannels");
     fLambda             = p.get<double>("Lambda"); //parameter for the denoising algorithm
     fMaxTicks           = p.get<int>("MaxTicks"); //maximum number of ticks in the waveform
+    fMaxTicksDent       = p.get<size_t>("MaxTicksDent"); //maximum number of ticks for dent formation
+    fMaxTicksSat        = p.get<size_t>("MaxTicksSat"); //maximum number of saturated ticks in the waveform to apply baseline removal
     fSecondBaselineSub  = p.get<double>("SecondBaselineSub"); //the mode of lowest SecondBaslineSub of the signal for second baseline estimate
     fFirstBaselineSub   = p.get<double>("FirstBaselineSub"); //around 3 times the expected large signals
+    fMaxDentDepth       = p.get<double>("MaxDentDepth");     //Maximum percentage of falling from saturation to be considered a laser failure
+    fDynamicRangeSaturation    = p.get<Float_t>("DynamicRangeSaturation");
+    fRemoveBaselineFluctuation = p.get<bool>("RemoveBaselineFluctuation"); 
 
     // This module produces (infrastructure piece)
     produces< std::vector< raw::OpDetWaveform > >();
@@ -128,10 +143,6 @@ namespace opdet {
     assert(wvfHandle.isValid());
 
     // Reserve a large enough array
-    // TODO: Commented out int totalsize to fix unused variable build error in clang.
-    //       Uncomment when implementing the full logic.
-    //int totalsize = 0;
-    //totalsize += wvfHandle->size();
     std::vector<double> fwaveform;
     fwaveform.reserve(fMaxTicks); 
 
@@ -147,13 +158,14 @@ namespace opdet {
         fwaveform[i] = wf[i];
       }
 
-      if (fIsPDVD && wf.ChannelNumber() > 3000 ) { //Excluding PMTs from the pre processing in PDVD
+      if(std::find(fIgnoreChannels.begin(), fIgnoreChannels.end(), wf.ChannelNumber()) != fIgnoreChannels.end()){
       }else{
-        if (fApplyDenoising){
+        DentCorrection(fwaveform, fDynamicRangeSaturation, fMaxTicksDent, wf.ChannelNumber(), fMaxDentDepth);
+        if(fRemoveBaselineFluctuation && !CheckSaturation(fwaveform, fDynamicRangeSaturation, fMaxTicksSat)){ BaselineExtractor(fwaveform);}
+        if(fApplyDenoising){
           double lambda = fLambda;
-          DenoisingAlgo(fwaveform, lambda);
+          DenoisingAlgo(fwaveform, lambda, fDynamicRangeSaturation);
         }
-        BaselineExtractor(fwaveform);
       }
       std::vector< short > waveformOfShorts = VectorOfDoublesToVectorOfShorts(fwaveform);
 
@@ -175,7 +187,7 @@ namespace opdet {
   void WaveformPreProcessing::BaselineExtractor(std::vector<double> &wf){
     std::vector<double> signal_base, waveform_full_bs;
     std::vector<double> temp = wf; 
-    //MA, baseline estimate, basline subtraction
+    //MA, baseline estimate, baseline subtraction
     std::vector<double> signalsma = ComputeMovingAverage(temp, 4);
     std::vector<double> basev = EstimateBaselineOpeningCentered(signalsma, fFirstBaselineSub, 1);
     for (size_t j=0; j<signalsma.size(); ++j){
@@ -188,7 +200,47 @@ namespace opdet {
     wf = waveform_full_bs;
   }
 
-  void WaveformPreProcessing::DenoisingAlgo(std::vector<double> &waveform, double lambda){
+  void WaveformPreProcessing::DentCorrection(std::vector<double> &wf, Float_t fDynamicRangeSaturation, size_t MaxTicksDent, int channel, double MaxDentDepth){
+  //Checks if the waveform is saturated and if it has a "dent" at the start of the saturation plateau (likely caused by a temporary laser failure). If yes, replace the dent by the maximum ADC value. 
+    Float_t Peak = *std::max_element(wf.begin(), wf.end()); //Find the highest adc value in the waveform
+    if(Peak >= fDynamicRangeSaturation){
+    // Iterate through the data to find if there is a "dent". We look for a pattern: High -> Low (Dent) -> High
+      for (size_t i = 1; i < wf.size(); i++) {
+        if (wf[i] < Peak) { // If the current point is below the plateau, check if it's flanked by saturated points. This may be a hole in the plateau, not a natural peak
+          if(wf[i - 1] >= Peak) {//check for saturation - dent - saturation
+            for (size_t j = i + 1; j < std::min(i + MaxTicksDent, wf.size()); j++) { // Look ahead to see if it saturates soon
+              if(wf[j] < MaxDentDepth*Peak) break; //There are two genuine saturated separated peaks within the interval
+              if (wf[j] >= Peak) {
+                wf[i] = Peak;
+                break;
+              }//end if
+            }// end for
+          }
+        }//end wf[i]<Peak
+      }//end for
+    }//end Peak>saturation
+  } 
+
+  bool WaveformPreProcessing::CheckSaturation(std::vector<double> wf, Float_t fDynamicRangeSaturation, size_t fMaxTicksSat){
+  //Checks if the waveform has a hit with a very large charge that prevents the baseline from returning to the original value
+    auto it = std::max_element(wf.begin(), wf.end()); //Find the highest adc value in the waveform
+    Float_t Peak = *it;
+    if(Peak < fDynamicRangeSaturation) return false; //wvf is not saturated
+    size_t consecutiveSatTicks = 0;
+    for(double sample : wf) {
+      if(sample >= fDynamicRangeSaturation){// We are inside a saturated region
+        consecutiveSatTicks++;
+        // Check if THIS specific peak has exceeded the limit
+        if (consecutiveSatTicks > fMaxTicksSat) return true; 
+      }else{ // Signal dropped below threshold; reset counter for the next peak
+        consecutiveSatTicks = 0;
+      }
+    }
+    // If we finished the entire waveform without the counter exceeding fMaxTicksSat
+    return false;
+  }
+
+  void WaveformPreProcessing::DenoisingAlgo(std::vector<double> &waveform, double lambda, Float_t fDynamicRangeSaturation){
     int length = waveform.size();
     std::vector<double> outwaveform; 
     outwaveform = waveform;  // copy
@@ -205,8 +257,8 @@ namespace opdet {
       }
     }
 
-    for(int i = 0; i < length; i++) { //replacing non-zero entries
-      if(outwaveform[i]) waveform[i] = outwaveform[i];
+    for(int i = 0; i < length; i++) { //replacing non-zero and non-saturated entries
+      if(outwaveform[i] && waveform[i]<fDynamicRangeSaturation) waveform[i] = outwaveform[i];
     }
   }
 
@@ -408,8 +460,13 @@ namespace opdet {
 
   std::vector<short> WaveformPreProcessing::VectorOfDoublesToVectorOfShorts (std::vector<double> const& vectorOfDoubles) 
   {
-    // Don't bother to round properly, it's faster this way
-    return std::vector<short>(vectorOfDoubles.begin(), vectorOfDoubles.end());
+    // Rounding properly
+    std::vector<short> result;
+    result.reserve(vectorOfDoubles.size());
+    for (double val : vectorOfDoubles) {
+      result.push_back(static_cast<short>(std::round(val)));
+    }
+    return result;
   }
 
 }//namespace opdet