From 528e2998ecc060c954806f72608676a1c1acfdba Mon Sep 17 00:00:00 2001 From: Yashwanth Bezawada Date: Tue, 21 Jan 2025 12:13:36 -0800 Subject: [PATCH] New waveform analysis class to interpolate waveforms using tapered sinc 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 | 16 +++ ratdb/IO.ratdb | 2 +- src/cmd/src/ProcBlockManager.cc | 4 +- src/daq/CMakeLists.txt | 1 + src/daq/include/RAT/WaveformAnalysisSinc.hh | 102 ++++++++++++++++ src/daq/src/WaveformAnalysisSinc.cc | 125 ++++++++++++++++++++ 6 files changed, 248 insertions(+), 2 deletions(-) create mode 100644 src/daq/include/RAT/WaveformAnalysisSinc.hh create mode 100644 src/daq/src/WaveformAnalysisSinc.cc diff --git a/ratdb/DIGITIZER_ANALYSIS.ratdb b/ratdb/DIGITIZER_ANALYSIS.ratdb index 8578d6fd..b2c07a3b 100644 --- a/ratdb/DIGITIZER_ANALYSIS.ratdb +++ b/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 } \ No newline at end of file diff --git a/ratdb/IO.ratdb b/ratdb/IO.ratdb index 62b8057f..f24627a9 100644 --- a/ratdb/IO.ratdb +++ b/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"] } diff --git a/src/cmd/src/ProcBlockManager.cc b/src/cmd/src/ProcBlockManager.cc index 7a3eb5a6..bc9943d6 100644 --- a/src/cmd/src/ProcBlockManager.cc +++ b/src/cmd/src/ProcBlockManager.cc @@ -29,6 +29,7 @@ #include #include #include +#include #include namespace RAT { @@ -94,8 +95,9 @@ ProcBlockManager::ProcBlockManager(ProcBlock *theMainBlock) { AppendProcessor(); AppendProcessor(); AppendProcessor(); - AppendProcessor(); AppendProcessor(); + AppendProcessor(); + AppendProcessor(); // Misc AppendProcessor(); AppendProcessor(); diff --git a/src/daq/CMakeLists.txt b/src/daq/CMakeLists.txt index 697c8d40..70c5f8f4 100644 --- a/src/daq/CMakeLists.txt +++ b/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 diff --git a/src/daq/include/RAT/WaveformAnalysisSinc.hh b/src/daq/include/RAT/WaveformAnalysisSinc.hh new file mode 100644 index 00000000..cabfb667 --- /dev/null +++ b/src/daq/include/RAT/WaveformAnalysisSinc.hh @@ -0,0 +1,102 @@ +//////////////////////////////////////////////////////////////////// +/// \class RAT::WaveformAnalysisSinc +/// +/// \brief Interpolate waveforms by convolution using sinc kernel +/// +/// \author Yashwanth Bezawada +/// \author Tanner Kaptanoglu +/// \author Ravi Pitelka +/// \author James Shen +/// +/// 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 + +#include +#include +#include +#include +#include +#include +#include + +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 ConvolveWaveform(const std::vector &wfm); + void InterpolateWaveform(const std::vector &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 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(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 &digitWfm) override; +}; + +} // namespace RAT + +#endif diff --git a/src/daq/src/WaveformAnalysisSinc.cc b/src/daq/src/WaveformAnalysisSinc.cc new file mode 100644 index 00000000..7636200d --- /dev/null +++ b/src/daq/src/WaveformAnalysisSinc.cc @@ -0,0 +1,125 @@ +#include +#include +#include + +#include +#include + +#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& 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 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 WaveformAnalysisSinc::ConvolveWaveform(const std::vector& wfm) { + int wl = wfm.size(); // waveform length + int numNewPoints = (wl - 1) * fNumInterpPoints + 1; + + std::vector 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& 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 sampleWfm; + int sampleLow = static_cast(floor(bf / fTimeStep)); + int sampleHigh = static_cast(floor(tf / fTimeStep)); + for (int j = sampleLow; j < sampleHigh + 1; j++) { + sampleWfm.push_back(voltWfm[j]); + } + + // Interpolated waveform + std::vector interpWfm = ConvolveWaveform(sampleWfm); + std::pair 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