diff --git a/src/cmd/src/ProcBlockManager.cc b/src/cmd/src/ProcBlockManager.cc index a161847b..9bf763dd 100644 --- a/src/cmd/src/ProcBlockManager.cc +++ b/src/cmd/src/ProcBlockManager.cc @@ -28,6 +28,7 @@ #include #include #include +#include #include namespace RAT { @@ -94,6 +95,7 @@ ProcBlockManager::ProcBlockManager(ProcBlock *theMainBlock) { AppendProcessor(); AppendProcessor(); AppendProcessor(); + AppendProcessor(); // Misc AppendProcessor(); AppendProcessor(); diff --git a/src/daq/CMakeLists.txt b/src/daq/CMakeLists.txt index 6a258266..a6263689 100644 --- a/src/daq/CMakeLists.txt +++ b/src/daq/CMakeLists.txt @@ -20,6 +20,7 @@ add_library(daq OBJECT src/WaveformPrep.cc src/WaveformAnalyzerBase.cc src/WaveformAnalysisLognormal.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..5f691185 --- /dev/null +++ b/src/daq/include/RAT/WaveformAnalysisSinc.hh @@ -0,0 +1,70 @@ +//////////////////////////////////////////////////////////////////// +/// \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. +//////////////////////////////////////////////////////////////////// +#ifndef __RAT_WaveformAnalysisSinc__ +#define __RAT_WaveformAnalysisSinc__ + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace RAT { + +class WaveformAnalysisSinc : public WaveformAnalyzerBase { + public: + WaveformAnalysisSinc() : WaveformAnalysisSinc(""){}; + WaveformAnalysisSinc(std::string config_name) : WaveformAnalyzerBase("WaveformAnalysisSinc", config_name) { + Configure(config_name); + }; + virtual ~WaveformAnalysisSinc(){}; + void Configure(const std::string &config_name) override; + virtual void SetD(std::string param, double value) override; + + // Interpolate the digitized waveform by convolution using a sinc kernel + std::vector convolve_wfm(const std::vector &wfm, const std::vector &kernel); + void InterpolateWaveform(const std::vector &voltWfm); + + protected: + // Digitizer settings + DBLinkPtr fDigit; + + // Analysis constants + double fFitWindowLow; + double fFitWindowHigh; + + // Coming from WaveformPrep + double fDigitTimeInWindow; + + // Fit variables + double fFitTime; + double fFitCharge; + double fFitPeak; + + 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..bde8816b --- /dev/null +++ b/src/daq/src/WaveformAnalysisSinc.cc @@ -0,0 +1,110 @@ +#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"); + } catch (DBNotFoundError) { + RAT::Log::Die("WaveformAnalysisSinc: Unable to find analysis parameters."); + } +} + +void WaveformAnalysisSinc::SetD(std::string param, double value) { + if (param == "fit_window_low") { + fFitWindowLow = value; + } else if (param == "fit_window_high") { + fFitWindowHigh = 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(); + + if (std::isinf(fDigitTimeInWindow) || fDigitTimeInWindow < 0) { + return; + } + + InterpolateWaveform(voltWfm); + + DS::WaveformAnalysisResult* fit_result = digitpmt->GetOrCreateWaveformAnalysisResult("Sinc"); + fit_result->AddPE(fFitTime, fFitCharge, {{"peak", fFitPeak}}); +} + +std::vector WaveformAnalysisSinc::convolve_wfm(const std::vector& wfm, + const std::vector& kernel) { + int n = wfm.size(); + int m = kernel.size(); + std::vector result(n + m - 1, 0.0); // Size of the output + // Perform the convolution + for (int i = 0; i < n; ++i) { + for (int j = 0; j < m; ++j) { + result[i + j] += wfm[i] * kernel[j]; + } + } + return result; +} + +void WaveformAnalysisSinc::InterpolateWaveform(const std::vector& voltWfm) { + // Fit range + double bf = fDigitTimeInWindow - fFitWindowLow; + double tf = fDigitTimeInWindow + fFitWindowHigh; + + // Check the fit range is within the digitizer window + bf = (bf > 0) ? bf : 0; + tf = (tf > voltWfm.size() * fTimeStep) ? voltWfm.size() * fTimeStep : tf; + + // Get samples values within the fit range + std::vector fit_wfm; + int sample_low = static_cast(floor(bf / fTimeStep)); + int sample_high = static_cast(floor(tf / fTimeStep)); + sample_high = std::min(static_cast(voltWfm.size()) - 1, sample_high); + for (int j = sample_low; j < sample_high + 1; j++) { + fit_wfm.push_back(voltWfm[j]); + } + + // Calculating tapered sinc kernel + std::vector tsinc_kernel; + int N = 8; // number of interpolated points per data point + double T = 30.0; // tapering constant + for (int k = -48; k < 49; k++) { + double val = k * 3.1415 / (float)N; + double sinc = sin(val); + if (k == 0) + sinc = 1.0; + else + sinc /= val; + tsinc_kernel.push_back(sinc * exp(-pow(k / T, 2))); + } + + // Interpolated waveform + std::vector interp_wfm = convolve_wfm(fit_wfm, tsinc_kernel); + std::pair peakSampleVolt = WaveformUtil::FindHighestPeak(interp_wfm); + + fFitPeak = peakSampleVolt.second; + fFitTime = bf + peakSampleVolt.first * fTimeStep / (float)N; + fFitCharge = 0; + for (auto& vlt : interp_wfm) { + fFitCharge += WaveformUtil::VoltagetoCharge(vlt, fTimeStep / (float)N, fTermOhms); // in pC + } +} + +} // namespace RAT