New waveform analysis class to interpolate waveforms using tapered si…

…nc kernel (#208)

* New waveform analysis class to interpolate waveforms using tapered sinc kernel.

* some changes

* some changes

* Updated tsinc interpolation code to correctly interpolate the waveform

* removed unnecessary cout line

* updates to the sinc interpolation code

ratdb/DIGITIZER_ANALYSIS.ratdb

@@ -46,4 +46,20 @@ fit_window_high: 15.0,
 gaussian_width: 1.5,
 gaussian_width_low: 0.5, // lower and upper limit for fitted sigma
 gaussian_width_high: 6.0,
+}
+
+{
+name: "DIGITIZER_ANALYSIS",
+index: "TSincInterpolation",
+valid_begin: [0,0],
+valid_end: [0,0],
+
+// Fit window, in ns
+fit_window_low: 10.0,
+fit_window_high: 15.0,
+
+// Interpolation variables
+num_interp_points: 8,
+tapering_constant: 30.0,
+num_sinc_lobes: 6
 }

ratdb/IO.ratdb

@@ -23,5 +23,5 @@ include_digitizerfits: true,
 include_digitizerwaveforms: false,
 include_untriggered_events: true,
 include_mchits: true,
-waveform_fitters: ["Lognormal", "Gaussian"]
+waveform_fitters: ["Lognormal", "Gaussian", "Sinc"]
 }

src/cmd/src/ProcBlockManager.cc

@@ -29,6 +29,7 @@
 #include <RAT/WaveformAnalysis.hh>
 #include <RAT/WaveformAnalysisGaussian.hh>
 #include <RAT/WaveformAnalysisLognormal.hh>
+#include <RAT/WaveformAnalysisSinc.hh>
 #include <RAT/WaveformPrep.hh>
 
 namespace RAT {
@@ -94,8 +95,9 @@ ProcBlockManager::ProcBlockManager(ProcBlock *theMainBlock) {
   AppendProcessor<TrueDAQProc>();
   AppendProcessor<WaveformAnalysis>();
   AppendProcessor<WaveformPrep>();
-  AppendProcessor<WaveformAnalysisLognormal>();
   AppendProcessor<WaveformAnalysisGaussian>();
+  AppendProcessor<WaveformAnalysisLognormal>();
+  AppendProcessor<WaveformAnalysisSinc>();
   // Misc
   AppendProcessor<CountProc>();
   AppendProcessor<PruneProc>();

src/daq/CMakeLists.txt

@@ -21,6 +21,7 @@ add_library(daq OBJECT
   src/WaveformAnalyzerBase.cc
   src/WaveformAnalysisLognormal.cc
   src/WaveformAnalysisGaussian.cc
+  src/WaveformAnalysisSinc.cc
   )
 
 # Set our include directories

src/daq/include/RAT/WaveformAnalysisSinc.hh

@@ -0,0 +1,102 @@
+////////////////////////////////////////////////////////////////////
+/// \class RAT::WaveformAnalysisSinc
+///
+/// \brief Interpolate waveforms by convolution using sinc kernel
+///
+/// \author Yashwanth Bezawada <[email protected]>
+/// \author Tanner Kaptanoglu <[email protected]>
+/// \author Ravi Pitelka <[email protected]>
+/// \author James Shen <[email protected]>
+///
+/// REVISION HISTORY:\n
+///     25 Nov 2024: Initial commit
+///
+/// \details
+/// This class provides full support for interpolating
+/// the waveforms by convolution using a tapered sinc (tsinc)
+/// kernel. Based on the implementation of the existing
+/// WaveformAnalysis class and WaveformAnalysisLognormal class.
+///
+/// Algorithm is from W. K. Warburton and W. Hennig,
+/// "New Algorithms for Improved Digital Pulse Arrival Timing With
+/// Sub-GSps ADCs," in IEEE Transactions on Nuclear Science, vol. 64,
+/// no. 12, pp. 2938-2950, Dec. 2017, doi: 10.1109/TNS.2017.2766074.
+///
+/// A window for interpolation is chosen around the peak of the
+/// digitized waveform. Interpolation is performed by convolving the
+/// samples within this window and the tsinc kernel. The algorithm
+/// for the discrete convolution code used here is similar to the
+/// one specified in the paper (optimized to make use of sinc being
+/// an even function).
+///
+////////////////////////////////////////////////////////////////////
+#ifndef __RAT_WaveformAnalysisSinc__
+#define __RAT_WaveformAnalysisSinc__
+
+#include <TObject.h>
+
+#include <RAT/DB.hh>
+#include <RAT/DS/DigitPMT.hh>
+#include <RAT/Digitizer.hh>
+#include <RAT/Processor.hh>
+#include <RAT/WaveformAnalyzerBase.hh>
+#include <cmath>
+#include <vector>
+
+namespace RAT {
+
+class WaveformAnalysisSinc : public WaveformAnalyzerBase {
+ public:
+  WaveformAnalysisSinc() : WaveformAnalysisSinc("TSincInterpolation"){};
+  WaveformAnalysisSinc(std::string config_name) : WaveformAnalyzerBase("WaveformAnalysisSinc", config_name) {
+    Configure(config_name);
+  };
+  virtual ~WaveformAnalysisSinc(){};
+  void Configure(const std::string &config_name) override;
+  virtual void SetI(std::string param, int value) override;
+  virtual void SetD(std::string param, double value) override;
+
+  // Interpolate the digitized waveform by convolution using a tapered sinc kernel
+  std::vector<double> ConvolveWaveform(const std::vector<double> &wfm);
+  void InterpolateWaveform(const std::vector<double> &voltWfm);
+
+ protected:
+  // Digitizer settings
+  DBLinkPtr fDigit;
+
+  // Analysis constants
+  double fFitWindowLow;
+  double fFitWindowHigh;
+  int fNumInterpPoints;   // number of interpolated points per data point (or the number of points in each lobe of
+                          // the tsinc kernel)
+  double fTaperingConst;  // tapering constant (determines how fast the kernel decays)
+  int fNumSincLobes;      // number of sinc lobes to be included in the kernel (determines the length of the kernel)
+
+  // Coming from WaveformPrep
+  double fDigitTimeInWindow;
+
+  // Fit variables
+  double fFitTime;
+  double fFitCharge;
+  double fFitPeak;
+
+  std::vector<double> tsinc_kernel;  // Stores the tsinc kernel
+  void calculateTSincKernel() {      // Calculates the tsinc kernel
+    for (int k = 0; k < fNumInterpPoints * fNumSincLobes + 1;
+         k++) {  // Storing only the positive k values since sinc is an even function
+      double val = k * M_PI / static_cast<float>(fNumInterpPoints);
+      double sinc = sin(val);
+      if (k == 0)
+        sinc = 1.0;
+      else
+        sinc /= val;
+      tsinc_kernel.push_back(sinc * exp(-pow(k / fTaperingConst, 2)));
+    }
+  }
+
+  void DoAnalysis(DS::DigitPMT *pmt, const std::vector<UShort_t> &digitWfm) override;
+};
+
+}  // namespace RAT
+
+#endif

src/daq/src/WaveformAnalysisSinc.cc

@@ -0,0 +1,125 @@
+#include <TF1.h>
+#include <TH1D.h>
+#include <TMath.h>
+
+#include <RAT/Log.hh>
+#include <RAT/WaveformAnalysisSinc.hh>
+
+#include "RAT/DS/DigitPMT.hh"
+#include "RAT/DS/WaveformAnalysisResult.hh"
+#include "RAT/WaveformUtil.hh"
+
+namespace RAT {
+
+void WaveformAnalysisSinc::Configure(const std::string& config_name) {
+  try {
+    fDigit = DB::Get()->GetLink("DIGITIZER_ANALYSIS", config_name);
+    fFitWindowLow = fDigit->GetD("fit_window_low");
+    fFitWindowHigh = fDigit->GetD("fit_window_high");
+    fNumInterpPoints = fDigit->GetI("num_interp_points");
+    fTaperingConst = fDigit->GetD("tapering_constant");
+    fNumSincLobes = fDigit->GetI("num_sinc_lobes");
+  } catch (DBNotFoundError) {
+    RAT::Log::Die("WaveformAnalysisSinc: Unable to find analysis parameters.");
+  }
+  calculateTSincKernel();
+}
+
+void WaveformAnalysisSinc::SetI(std::string param, int value) {
+  if (param == "num_interp_points") {
+    fNumInterpPoints = value;
+  } else if (param == "num_sinc_lobes") {
+    fNumSincLobes = value;
+  } else {
+    throw Processor::ParamUnknown(param);
+  }
+}
+
+void WaveformAnalysisSinc::SetD(std::string param, double value) {
+  if (param == "fit_window_low") {
+    fFitWindowLow = value;
+  } else if (param == "fit_window_high") {
+    fFitWindowHigh = value;
+  } else if (param == "tapering_constant") {
+    fTaperingConst = value;
+  } else {
+    throw Processor::ParamUnknown(param);
+  }
+}
+
+void WaveformAnalysisSinc::DoAnalysis(DS::DigitPMT* digitpmt, const std::vector<UShort_t>& digitWfm) {
+  double pedestal = digitpmt->GetPedestal();
+  if (pedestal == -9999) {
+    RAT::Log::Die("WaveformAnalysisSinc: Pedestal is invalid! Did you run WaveformPrep first?");
+  }
+  // Convert from ADC to mV
+  std::vector<double> voltWfm = WaveformUtil::ADCtoVoltage(digitWfm, fVoltageRes, pedestal = pedestal);
+  fDigitTimeInWindow = digitpmt->GetDigitizedTimeNoOffset();
+
+  // check to see if digitTime is well behaved
+  if (std::isinf(fDigitTimeInWindow) || fDigitTimeInWindow < 0 || fDigitTimeInWindow > voltWfm.size() * fTimeStep) {
+    return;
+  }
+
+  InterpolateWaveform(voltWfm);
+
+  DS::WaveformAnalysisResult* fit_result = digitpmt->GetOrCreateWaveformAnalysisResult("Sinc");
+  fit_result->AddPE(fFitTime, fFitCharge, {{"peak", fFitPeak}});
+}
+
+std::vector<double> WaveformAnalysisSinc::ConvolveWaveform(const std::vector<double>& wfm) {
+  int wl = wfm.size();  // waveform length
+  int numNewPoints = (wl - 1) * fNumInterpPoints + 1;
+
+  std::vector<double> newWfm(numNewPoints, 0.0);  // vector that stored the interpolated waveform
+  for (int j = 0; j < wl - 1; j++) {              // looping over the digitized waveform
+    for (int k = 0; k < fNumInterpPoints; k++)    // looping over the interpolated points
+    {
+      double interpValue = 0.;
+      for (int i = 0; i < fNumSincLobes; i++)  // looping over the tsinc kernel
+      {
+        if (j - i >= 0) {
+          interpValue += wfm[j - i] * tsinc_kernel[i * fNumInterpPoints + k];
+        }
+        if (j + i + 1 < wl) {
+          interpValue += wfm[j + 1 + i] * tsinc_kernel[(i + 1) * fNumInterpPoints - k];
+        }
+      }
+      newWfm[j * fNumInterpPoints + k] = interpValue;
+    }
+  }
+  newWfm[numNewPoints - 1] =
+      wfm[wl - 1];  // End point of the digitized waveform is the end point of the interpolated waveform
+  return newWfm;
+}
+
+void WaveformAnalysisSinc::InterpolateWaveform(const std::vector<double>& voltWfm) {
+  // Interpolation range
+  double bf = fDigitTimeInWindow - fFitWindowLow;
+  double tf = fDigitTimeInWindow + fFitWindowHigh;
+
+  // Check the interpolation range is within the digitizer window
+  bf = (bf > 0) ? bf : 0;
+  tf = (tf >= voltWfm.size() * fTimeStep) ? (voltWfm.size() * fTimeStep) - (fTimeStep / 2.0) : tf;
+
+  // Get samples values within the interpolation range
+  std::vector<double> sampleWfm;
+  int sampleLow = static_cast<int>(floor(bf / fTimeStep));
+  int sampleHigh = static_cast<int>(floor(tf / fTimeStep));
+  for (int j = sampleLow; j < sampleHigh + 1; j++) {
+    sampleWfm.push_back(voltWfm[j]);
+  }
+
+  // Interpolated waveform
+  std::vector<double> interpWfm = ConvolveWaveform(sampleWfm);
+  std::pair<int, double> peakSampleVolt = WaveformUtil::FindHighestPeak(interpWfm);
+
+  fFitPeak = peakSampleVolt.second;
+  fFitTime = ((sampleLow + 0.5) * fTimeStep) + peakSampleVolt.first * fTimeStep / fNumInterpPoints;
+  fFitCharge = 0;
+  for (const double& vlt : interpWfm) {
+    fFitCharge += WaveformUtil::VoltagetoCharge(vlt, fTimeStep / fNumInterpPoints, fTermOhms);  // in pC
+  }
+}
+
+}  // namespace RAT