From e49c4307a149977c1b18b99d6a94b716ea37eb9a Mon Sep 17 00:00:00 2001 From: LassiW96 Date: Wed, 22 Jul 2026 10:14:21 -0500 Subject: [PATCH] Tests that use PE based waveform generator algorithm added. --- tests/UnitTest_SA.cxx | 109 +++----- tests/waveformGenerator.cxx | 537 +++++++++++++++++++++++------------- tests/waveformGenerator.h | 132 +++++---- 3 files changed, 470 insertions(+), 308 deletions(-) diff --git a/tests/UnitTest_SA.cxx b/tests/UnitTest_SA.cxx index 66e8cc52..9ca91bc4 100644 --- a/tests/UnitTest_SA.cxx +++ b/tests/UnitTest_SA.cxx @@ -1,4 +1,5 @@ #include "../Podd/FadcStandaloneAnalyzer.h" +#include "RtypesCore.h" #include "waveformGenerator.h" #include #include @@ -16,37 +17,33 @@ std::vector ToDoubleWave(const std::vector &wave) { FADCConfig GetTestConfig() { FADCConfig conf; conf.fNPedestalSamples = 4; - conf.fSampThreshold = 50.0; + conf.fSampThreshold = 25.0; conf.fNSB = 1; - conf.fNSA = 2; + conf.fNSA = 4; conf.fNSAT = 1; conf.fMaxNPulses = 3; return conf; } // ─── Tests ─────────────────────────────────────────────────────────────────── - -TEST_CASE("WaveformAnalyzer correctly processes fADC data", "[analyzer]") { - +// Waveform generator v2 tests +TEST_CASE( + "WaveformAnalyzer correctly processes waveforms generated by generator v2", + "[analyzer]") { FadcStandaloneAnalyzer analyzer; analyzer.setConfig(GetTestConfig()); FADCCalib calib = {0.2441, 4.0 / 50.0}; analyzer.setCalib(calib); - WaveformData wfData; - wfData.TotNSamples = 100; - wfData.InitPedSamples = 4; - wfData.Pedestal = 100; - wfData.PedRange = 4; - wfData.NPulses = 4; - wfData.PulseAmp = 300; - wfData.PulseWidth = 5; - wfData.PulseSeparation = 10; + const FADCConfig &conf = analyzer.getConfig(); + + // Keeping default values + WaveformParams wfParams; + DigitizerParams daqParams; waveformGenerator gen; - gen.setWaveformData(wfData); - gen.setAnalyzer(analyzer); // inject same config so singleWaveIntegral uses - // identical FADCConfig + gen.setDetectorParams(wfParams); + gen.setDigitizerParams(daqParams); SECTION("Empty waveform returns no pulses and doesn't crash") { std::vector empty_wave; @@ -61,69 +58,43 @@ TEST_CASE("WaveformAnalyzer correctly processes fADC data", "[analyzer]") { REQUIRE(results.fNSampPulses == 0); } - SECTION("Single pulse waveform is detected with a positive " - "pedestal-subtracted integral") { - std::vector int_wave = gen.singlePulseWaveform(); - std::vector wave = ToDoubleWave(int_wave); - - // --- FADCStandaloneAnalyzer result --- - auto results = analyzer.Analyze(wave, 1.0, 4.0); - - // At least one pulse must be detected. - REQUIRE(results.fNSampPulses != 0); - - // Pulse parameters using waveformGenerator class (for comparison) - gen.clearWaveformInfo(); - gen.calcPulseParameters(int_wave); - const WaveformInfo &wfInfo = gen.getWaveformInfo(); - REQUIRE(results.fSampPulseIntPedSub[0] == - Catch::Approx(static_cast(wfInfo.PedSubIntegral[0])) - .margin(wfData.PedRange)); - REQUIRE(results.fSampPulseAmp[0] == - Catch::Approx(static_cast(wfInfo.PulseAmp[0])) - .margin(wfData.PedRange)); - - // std::cout << "\nSingle Pulse:\nPulse Time: From Waveform generator " << - // wfInfo.PulseTime[0] << " from FADCStandaloneAnalyzer " << - // results.fSampPulseTime[0] << "\n"; - REQUIRE(results.fSampPulseTime[0] == - Catch::Approx(static_cast(wfInfo.PulseTime[0])) - .margin(wfData.PedRange)); - } - - SECTION("Multiple pulse waveform is detected with a positive " - "pedestal-subtracted integral") { - std::vector int_wave = gen.multiplePulseWaveform(); + SECTION("Generate exponential waveform") { + std::vector int_wave = gen.generateWaveform(); std::vector wave = ToDoubleWave(int_wave); // --- FADCStandaloneAnalyzer result --- auto results = analyzer.Analyze(wave, 1.0, 4.0); // Output shouldn't be empty - REQUIRE(results.fNSampPulses >= 2); - + REQUIRE(results.fNSampPulses >= 1); // Pulse parameters using waveformGenerator class (for comparison) - gen.clearWaveformInfo(); - gen.calcPulseParameters(int_wave); - const WaveformInfo &wfInfo = gen.getWaveformInfo(); - - // std::cout << "\nCalculated Integrals and times for Multiple Pulses:\n"; - for (int i = 0; i < (int)wfInfo.PedSubIntegral.size(); i++) { - // std::cout << "\tPulse " << i + 1 << ": " << wfInfo.PedSubIntegral[i] << - // std::endl; + // gen.clearTruthInfo(); + const TruthPulseInfo &wfInfo = gen.getTruthInfo(); + + Double_t int_margin = (daqParams.noise_sigma) * (conf.fNSA + conf.fNSB) * + daqParams.sample_rate; + Double_t amp_margin = daqParams.noise_sigma * 2; + Double_t time_margin = 64.0 * (wfParams.tts_sigma / daqParams.sample_rate); + + std::cout << "\nCalculated Integrals and times for Multiple Pulses:\n"; + for (int i = 0; i < (int)wfInfo.pulseIntegral.size(); i++) { + std::cout << "\tPulse Number: " << i + 1 << std::endl; + std::cout << "\tPulse Integral: From Waveform generator: " + << wfInfo.pulseIntegral[i] << " from FADCStandaloneAnalyzer " + << results.fSampPulseIntPedSub[i] << std::endl; REQUIRE(results.fSampPulseIntPedSub[i] == - Catch::Approx(static_cast(wfInfo.PedSubIntegral[i])) - .margin(wfData.PedRange)); + Catch::Approx(static_cast(wfInfo.pulseIntegral[i])) + .margin(int_margin)); REQUIRE(results.fSampPulseAmp[i] == - Catch::Approx(static_cast(wfInfo.PulseAmp[i])) - .margin(wfData.PedRange)); + Catch::Approx(static_cast(wfInfo.pulseAmplitude[i])) + .margin(amp_margin)); - // std::cout << "\tPulse Time: From Waveform generator " << - // wfInfo.PulseTime[i] << " from FADCStandaloneAnalyzer " << - // results.fSampPulseTime[i] << "\n\n"; + std::cout << "\tPulse Time: From Waveform generator " + << wfInfo.pulseTime[i] << " from FADCStandaloneAnalyzer " + << results.fSampPulseTime[i] << "\n\n"; REQUIRE(results.fSampPulseTime[i] == - Catch::Approx(static_cast(wfInfo.PulseTime[i])) - .margin(wfData.PedRange)); + Catch::Approx(static_cast(wfInfo.pulseTime[i])) + .margin(time_margin)); } } } \ No newline at end of file diff --git a/tests/waveformGenerator.cxx b/tests/waveformGenerator.cxx index 194a4859..badab668 100644 --- a/tests/waveformGenerator.cxx +++ b/tests/waveformGenerator.cxx @@ -1,229 +1,390 @@ #include "waveformGenerator.h" -#include "../Podd/FadcStandaloneAnalyzer.h" #include "RtypesCore.h" +#include +#include +#include #include waveformGenerator::waveformGenerator() : fRng(std::random_device{}()) {} std::vector waveformGenerator::noiseWaveform() const { - int min = fWaveformData.Pedestal - fWaveformData.PedRange / 2; - int max = fWaveformData.Pedestal + fWaveformData.PedRange / 2 + 1; - std::uniform_real_distribution rndm(0.0, 1.0); - std::vector wave; - for (int i = 0; i < fWaveformData.TotNSamples; i++) { + int min = fDaq.pedestal - fDaq.ped_sigma / 2.0; + int max = fDaq.pedestal + fDaq.ped_sigma / 2.0 + 1; + std::uniform_real_distribution rndm(0.0, 1.0); + std::vector wave; + for( int i = 0; i < fDaq.num_samples; i++ ) { wave.push_back(min + (max - min) * rndm(fRng)); - } - return wave; + } + return wave; } std::vector waveformGenerator::singlePulseWaveform() const { - int min_ped = fWaveformData.Pedestal - fWaveformData.PedRange / 2; - int max_ped = fWaveformData.Pedestal + fWaveformData.PedRange / 2 + 1; - std::uniform_real_distribution rndm(0.0, 1.0); - std::vector wave; - - int pulse_min = fWaveformData.InitPedSamples; - int pulse_max = fWaveformData.TotNSamples - 1; - int random_idx = pulse_min + int(pulse_max - pulse_min) * rndm(fRng); - for (int i = 0; i < fWaveformData.TotNSamples; i++) { - if (i >= random_idx && i < random_idx + fWaveformData.PulseWidth) { - int j = i - random_idx; - wave.push_back(min_ped + fWaveformData.PulseAmp / TMath::Power(2, j)); - } else { - wave.push_back(min_ped + (max_ped - min_ped) * rndm(fRng)); + int min_ped = fDaq.pedestal - fDaq.ped_sigma / 2.0; + int max_ped = fDaq.pedestal + fDaq.ped_sigma / 2.0 + 1; + std::uniform_real_distribution rndm(0.0, 1.0); + std::vector wave; + + int pulse_min = fDaq.init_ped_samples; + int pulse_max = fDaq.num_samples - 1; + int random_idx = pulse_min + int(pulse_max - pulse_min) * rndm(fRng); + for( int i = 0; i < fDaq.num_samples; i++ ) { + if(i >= random_idx && i < random_idx + fWf.pulse_width) { + int j = i - random_idx; + wave.push_back(min_ped + fWf.mean_gain/TMath::Power(2, j)); + } else { + wave.push_back(min_ped + (max_ped - min_ped) * rndm(fRng)); + } } - } - return wave; + return wave; } std::vector waveformGenerator::multiplePulseWaveform() const { - int min_ped = fWaveformData.Pedestal - fWaveformData.PedRange / 2; - int max_ped = fWaveformData.Pedestal + fWaveformData.PedRange / 2 + 1; - std::uniform_real_distribution rndm(0.0, 1.0); - std::vector wave; - std::vector> pulse_info; + int min_ped = fDaq.pedestal - fDaq.ped_sigma / 2.0; + int max_ped = fDaq.pedestal + fDaq.ped_sigma / 2.0 + 1; + std::uniform_real_distribution rndm(0.0, 1.0); + std::vector wave; + std::vector> pulse_info; - int num_pulses = 2 + int((fWaveformData.NPulses - 2) * rndm(fRng)); + int num_pulses = 2 + int((fDaq.max_n_pulses - 2) * rndm(fRng)); - // Maximum pulse space required - int max_space = (num_pulses * fWaveformData.PulseWidth) + - ((num_pulses - 1) * fWaveformData.PulseSeparation); + // Maximum pulse space required + int max_space = (num_pulses * fWf.pulse_width) + + ((num_pulses - 1) * fWf.pulse_separation); - // Maximum available space for pulse start - int max_first_start = fWaveformData.TotNSamples - max_space; + // Maximum available space for pulse start + int max_first_start = fDaq.num_samples - max_space; - if (max_first_start < fWaveformData.InitPedSamples) { - max_first_start = fWaveformData.InitPedSamples; - } + if(max_first_start < fDaq.init_ped_samples) { + max_first_start = fDaq.init_ped_samples; + } - // Current starting point - int current_start = - fWaveformData.InitPedSamples + - int((max_first_start - fWaveformData.InitPedSamples) * rndm(fRng)); + // Current starting point + int current_start = fDaq.init_ped_samples + + int((max_first_start - fDaq.init_ped_samples) * rndm(fRng)); - // Random amps and separations - for (int p = 0; p < num_pulses; p++) { - int min_amp = fWaveformData.PulseAmp / 2; - int current_amp = - min_amp + int((fWaveformData.PulseAmp - min_amp) * rndm(fRng)); + // Random amps and separations + for(int p = 0; p < num_pulses; p++) { + int min_amp = fWf.mean_gain / 2; + int current_amp = min_amp + int((fWf.mean_gain - min_amp) * rndm(fRng)); pulse_info.push_back({current_start, current_amp}); - int current_separation = fWaveformData.PulseSeparation / 2 + - int((fWaveformData.PulseSeparation - - fWaveformData.PulseSeparation / 2.) * - rndm(fRng)); - current_start += fWaveformData.PulseWidth + current_separation; - } + int current_separation = fWf.pulse_separation / 2 + + int((fWf.pulse_separation - fWf.pulse_separation / 2.) * rndm(fRng)); + current_start += fWf.pulse_width + current_separation; + } - // Build the waveform - for (int i = 0; i < fWaveformData.TotNSamples; i++) { + // Build the waveform + for( int i = 0; i < fDaq.num_samples; i++) { Bool_t in_pulse = false; Int_t pulse_relative_idx = 0; Int_t active_amp = 0; - for (auto &p : pulse_info) { - if (i >= p.first && i < p.first + fWaveformData.PulseWidth) { - in_pulse = true; - pulse_relative_idx = i - p.first; - active_amp = p.second; - break; - } + for(auto& p : pulse_info) { + if(i >= p.first && i < p.first + fWf.pulse_width) { + in_pulse = true; + pulse_relative_idx = i - p.first; + active_amp = p.second; + break; + } } - if (in_pulse) { - wave.push_back(min_ped + - active_amp / TMath::Power(2, pulse_relative_idx)); + if(in_pulse) { + wave.push_back(min_ped + active_amp/TMath::Power(2, pulse_relative_idx)); } else { - wave.push_back(min_ped + (max_ped - min_ped) * rndm(fRng)); + wave.push_back(min_ped + (max_ped - min_ped) * rndm(fRng)); + } + } + return wave; +} + +// PE based waveform generators +// Waveform generator v2 related functions - now the parameters are read from the object state +Double_t waveformGenerator::spePulse(Double_t t, Double_t t0, Double_t amp) const { + if(t < t0) return 0.0; + return amp * (std::exp(-(t - t0) / fWf.tau_f) - std::exp(-(t - t0) / fWf.tau_r)); +} + +Double_t waveformGenerator::squarePulse(Double_t t, Double_t t0, Double_t amp) const { + Double_t width_ns = fWf.pulse_width * fDaq.sample_rate; // Convert samples -> ns + if(t < t0 || t > t0 + width_ns) return 0.0; + return amp; +} + +Double_t waveformGenerator::trianglePulse(Double_t t, Double_t t0, Double_t amp) const { + Double_t width_ns = fWf.pulse_width * fDaq.sample_rate; + Int_t dt = fDaq.sample_rate; + t0 = round(t0 / dt) * dt; + Double_t half_w = width_ns / 2.0; + if(t < t0 - half_w || t > t0 + half_w) return 0.0; + + if(t < t0) { + return amp * (t - (t0 - half_w)) / (half_w); + } else { + return amp * ((t0 + half_w) - t) / (half_w); + } +} + +std::vector waveformGenerator::generatePEHits() const { + // Distributions for event and after pulse gen + std::uniform_int_distribution num_pulses_dist(1, fDaq.max_n_pulses); + std::uniform_real_distribution window_dist(10.0, fDaq.num_samples * fDaq.sample_rate - 50.0); + std::poisson_distribution pe_per_hit_dist(15); // Avg 15 primary PEs per hit + std::uniform_real_distribution ap_chance(0.0, 1.0); // Prob for after pulse + std::normal_distribution ap_time_dist(fWf.ap_delay_mean, fWf.ap_delay_sigma); // After pulse time dist + + std::vector hits; + Int_t num_pulses = num_pulses_dist(fRng); + + for (Int_t p = 0; p < num_pulses; p++) { + PEHit hit; + hit.hitTime = window_dist(fRng); + Int_t num_pes = pe_per_hit_dist(fRng); + + for (Int_t i = 0; i < num_pes; i++) { + hit.peTimes.push_back(hit.hitTime); + + // After pulse probability roll + // if(ap_chance(fRng) < fWf.ap_prob) { + // Double_t delayed_time = hit.hitTime + ap_time_dist(fRng); + // hit.peTimes.push_back(delayed_time); + // } + } + // Sort PE times within this hit + std::sort(hit.peTimes.begin(), hit.peTimes.end()); + hits.push_back(hit); } - } - return wave; + // Sort hits by arrival time + std::sort(hits.begin(), hits.end(), + [](const PEHit& a, const PEHit& b) { return a.hitTime < b.hitTime; }); + return hits; } -// std::vector waveformGenerator::calcPulseIntegral(const -// std::vector wave) const { -// FADCConfig fConfig = saAnalyzer.getConfig(); -// Int_t NPedSamps = fConfig.fNPedestalSamples; -// Int_t NSB = fConfig.fNSB; -// Int_t NSA = fConfig.fNSA; -// Int_t Threshold = fConfig.fSampThreshold; -// Int_t NSamples = (Int_t)wave.size(); - -// std::vector AvgIntegrals = {}; - -// // Pedestal: raw sum of first NPedSamps samples (matches analyzer's -// fSampPed) Int_t PedSum = 0; for (int i = 0; i < NPedSamps; i++) -// PedSum += wave[i]; -// Int_t AvgPed = PedSum / NPedSamps; - -// Bool_t PulseFound = false; - -// for(int i = 0; i < wave.size(); i++) { -// if(wave[i] - AvgPed > Threshold) { -// if(!PulseFound) { -// // Integration window: [FirstSamp - NSB, FirstSamp + NSA - -// 1], clamped to valid range Int_t lo = TMath::Max(i - NSB, 0); -// Int_t hi = TMath::Min(i + NSA - 1, NSamples - 1); -// // Pedestal-subtracted integral, using the same formula as -// the analyzer: -// // PedSubInt = RawInt - PedSum * (window_width / NPedSamps) -// // = RawInt - AvgPed * window_width -// Int_t WindowWidth = hi - lo + 1; -// Int_t AvgIntegral = 0; -// Int_t RawIntegral = 0; -// for (int j = lo; j <= hi; j++) { -// RawIntegral += wave[j]; -// } -// AvgIntegral = RawIntegral - AvgPed * WindowWidth; -// PulseFound = true; -// AvgIntegrals.push_back(AvgIntegral); -// } -// } else { -// PulseFound = false; -// } -// } -// return AvgIntegrals; -// } - -void waveformGenerator::calcPulseParameters(const std::vector wave) { - FADCConfig fConfig = saAnalyzer.getConfig(); - Int_t NPedSamps = fConfig.fNPedestalSamples; - Int_t NSB = fConfig.fNSB; - Int_t NSA = fConfig.fNSA; - Int_t Threshold = fConfig.fSampThreshold; - Int_t NSamples = (Int_t)wave.size(); - - std::vector AvgIntegrals = {}; - - // Pedestal: raw sum of first NPedSamps samples (matches analyzer's fSampPed) - Double_t PedSum = 0.0; - for (int i = 0; i < NPedSamps; i++) - PedSum += wave[i]; - // Use the same formula as the analyzer: - // AvgPed (float) = PedSum / NPedSamps — NOT integer division - Double_t AvgPed = PedSum / NPedSamps; - - Bool_t PulseFound = false; - - for (int i = 0; i < wave.size(); i++) { - if (wave[i] - AvgPed > Threshold) { - if (!PulseFound) { - // Integration window: [FirstSamp - NSB, FirstSamp + NSA - 1], clamped - // to valid range - Int_t lo = TMath::Max(i - NSB, 0); - Int_t hi = TMath::Min(i + NSA - 1, NSamples - 1); - // Pedestal-subtracted integral, using the same formula as the analyzer: - // PedSubInt = RawInt - PedSum * (window_width / NPedSamps) - // = RawInt - AvgPed * window_width - Int_t WindowWidth = hi - lo + 1; - Double_t AvgIntegral = 0.0; - Int_t RawIntegral = 0; - for (int j = lo; j <= hi; j++) { - RawIntegral += wave[j]; +std::vector waveformGenerator::generateWaveform() { + auto hits = generatePEHits(); + fTruthInfo = {}; // Auto-clear truth info + + // Distributions + std::normal_distribution tts_dist(0.0, fWf.tts_sigma); + std::normal_distribution gain_dist(fWf.mean_gain, fWf.gain_sigma); + std::normal_distribution noise_dist(0.0, fDaq.noise_sigma); + + std::vector analog_wf(fDaq.num_samples, 0.0); + + for(auto& hit : hits) { + // Per-hit analog contribution for truth extraction + std::vector hit_contrib(fDaq.num_samples, 0.0); + + for(Double_t pe_time : hit.peTimes) { + Double_t smeared_time = pe_time + tts_dist(fRng); + Double_t smeared_gain = gain_dist(fRng); + if(smeared_gain < 0.0) smeared_gain = 0.0; + + for(Int_t i = 0; i < fDaq.num_samples; i++) { + Double_t current_time = i * fDaq.sample_rate; + // Optimization + if(current_time > smeared_time - 5.0 && current_time < smeared_time + 50) { + Double_t val = spePulse(current_time, smeared_time, smeared_gain); + hit_contrib[i] += val; + analog_wf[i] += val; + } + } } - // Match analyzer: PedSubInt = RawInt - PedSum * (WindowWidth / - // NPedSamps) - AvgIntegral = RawIntegral - - PedSum * (static_cast(WindowWidth) / NPedSamps); - PulseFound = true; - fWaveformInfo.PedSubIntegral.push_back(AvgIntegral); - - Int_t pulse_amp = 0; - Int_t pulse_Time = 0; - - Int_t PeakBin = 0; - Double_t PeakVal = wave[i] - AvgPed; - for (Int_t nt = i + 1; nt < TMath::Min((i + NSA), int(NSamples)); - nt++) { - if ((wave[nt] - AvgPed) < PeakVal && PeakBin == 0) { - PeakBin = nt - 1; - } else if (PeakBin == 0) { - PeakVal = wave[nt] - AvgPed; - } + + // Extract truth features from per-hit contribution + // Find peak amplitude and bin + Double_t peakAmp = 0.0; + Int_t peakBin = 0; + for(Int_t i = 0; i < fDaq.num_samples; i++) { + if(hit_contrib[i] > peakAmp) { + peakAmp = hit_contrib[i]; + peakBin = i; + } } - if (PeakBin > 0) { - pulse_amp = wave[PeakBin]; - Int_t Time = i * 64; - Double_t VMid = (wave[PeakBin] - AvgPed) / 2.; - for (Int_t nt = TMath::Max(i - fConfig.fNSB, 0); - nt < TMath::Min(PeakBin, int(NSamples - 1)); nt++) { - if (VMid >= (wave[nt] - AvgPed) && VMid < (wave[nt + 1] - AvgPed)) { - Time = 64 * nt + - int(64 * (VMid - (wave[nt] - AvgPed)) / - ((wave[nt + 1] - AvgPed) - (wave[nt] - AvgPed))); + + // Integral: sum over non-negligible samples (>1% of peak), + // converted to ADC counts (pedestal-subtracted) + Double_t integral = 0.0; + Double_t threshold = peakAmp * 0.01; + for(Int_t i = 0; i < fDaq.num_samples; i++) { + if(hit_contrib[i] > threshold) { + integral += hit_contrib[i]; } - } - pulse_Time = Time; - } else { - pulse_amp = wave[i]; - pulse_Time = 64 * i; } + integral /= fDaq.v_lsb; // mV -> ADC counts - fWaveformInfo.PulseAmp.push_back(pulse_amp); - fWaveformInfo.PulseTime.push_back(pulse_Time); - } - } else { - PulseFound = false; + // Amplitude: convert to ADC with pedestal (matches fSampPulseAmp) + Double_t ampADC = peakAmp / fDaq.v_lsb + fDaq.pedestal; + + // Time: half-max crossing on rising edge, in 1/64th-sample units + // (matches analyzer's fSampPulseTime interpolation) + Double_t vMid = peakAmp / 2.0; + Double_t timeVal = peakBin * 64.0; // fallback: peak bin + for(Int_t i = 0; i < peakBin; i++) { + if(hit_contrib[i] <= vMid && hit_contrib[i + 1] > vMid) { + timeVal = 64.0 * i + 64.0 * (vMid - hit_contrib[i]) / + (hit_contrib[i + 1] - hit_contrib[i]); + break; + } + } + + fTruthInfo.pulseTime.push_back(timeVal); + fTruthInfo.pulseAmplitude.push_back(ampADC); + fTruthInfo.pulseIntegral.push_back(integral); + } + + // Add noise and digitize + std::vector adc_wf(fDaq.num_samples, 0); + for(Int_t i = 0; i < fDaq.num_samples; i++) { + Double_t noisy_voilage = analog_wf[i] + noise_dist(fRng); + Int_t adc_count = static_cast(std::floor(noisy_voilage/fDaq.v_lsb) + fDaq.pedestal); + adc_count = std::max(0, std::min(adc_count, fDaq.max_adc)); + adc_wf[i] = adc_count; + } + return adc_wf; +} + +std::vector waveformGenerator::squareWaveform() { + auto hits = generatePEHits(); + fTruthInfo = {}; // Auto-clear truth info + + std::vector square_wf(fDaq.num_samples, 0.0); + + for(auto& hit : hits) { + // Per-hit analog contribution for truth extraction + std::vector hit_contrib(fDaq.num_samples, 0.0); + + for(Double_t pe_time : hit.peTimes) { + for(Int_t i = 0; i < fDaq.num_samples; i++) { + Double_t current_time = i * fDaq.sample_rate; + // Optimization + if(current_time > pe_time - 5.0 && current_time < pe_time + 50) { + Double_t val = squarePulse(current_time, pe_time, fWf.mean_gain); + hit_contrib[i] += val; + square_wf[i] += val; + } + } + } + + // Extract truth features from per-hit contribution + // Find peak amplitude and bin + Double_t peakAmp = 0.0; + Int_t peakBin = 0; + for(Int_t i = 0; i < fDaq.num_samples; i++) { + if(hit_contrib[i] > peakAmp) { + peakAmp = hit_contrib[i]; + peakBin = i; + } + } + + // Integral: sum over non-negligible samples (>1% of peak), + // converted to ADC counts (pedestal-subtracted) + Double_t integral = 0.0; + Double_t threshold = peakAmp * 0.01; + for(Int_t i = 0; i < fDaq.num_samples; i++) { + if(hit_contrib[i] > threshold) { + integral += hit_contrib[i]; + } + } + integral /= fDaq.v_lsb; // mV -> ADC counts + + // Amplitude: convert to ADC with pedestal (matches fSampPulseAmp) + Double_t ampADC = peakAmp / fDaq.v_lsb + fDaq.pedestal; + + // Time: half-max crossing on rising edge, in 1/64th-sample units + Double_t vMid = peakAmp / 2.0; + Double_t timeVal = peakBin * 64.0; // fallback: peak bin + for(Int_t i = 0; i < peakBin; i++) { + if(hit_contrib[i] <= vMid && hit_contrib[i + 1] > vMid) { + timeVal = 64.0 * i + 64.0 * (vMid - hit_contrib[i]) / + (hit_contrib[i + 1] - hit_contrib[i]); + break; + } + } + + fTruthInfo.pulseTime.push_back(timeVal); + fTruthInfo.pulseAmplitude.push_back(ampADC); + fTruthInfo.pulseIntegral.push_back(integral); + } + + // Digitize + std::vector adc_wf(fDaq.num_samples, 0); + for(Int_t i = 0; i < fDaq.num_samples; i++) { + Int_t adc_count = static_cast(std::floor(square_wf[i]/fDaq.v_lsb) + fDaq.pedestal); + adc_count = std::max(0, std::min(adc_count, fDaq.max_adc)); + adc_wf[i] = adc_count; } - } + return adc_wf; } + +std::vector waveformGenerator::triangleWaveform() { + auto hits = generatePEHits(); + fTruthInfo = {}; // Auto-clear truth info + + std::vector triangle_wf(fDaq.num_samples, 0.0); + + for(auto& hit : hits) { + // Per-hit analog contribution for truth extraction + std::vector hit_contrib(fDaq.num_samples, 0.0); + + for(Double_t pe_time : hit.peTimes) { + for(Int_t i = 0; i < fDaq.num_samples; i++) { + Double_t current_time = i * fDaq.sample_rate; + // Optimization + if(current_time > pe_time - 5.0 && current_time < pe_time + 50) { + Double_t val = trianglePulse(current_time, pe_time, fWf.mean_gain); + hit_contrib[i] += val; + triangle_wf[i] += val; + } + } + } + + // Extract truth features from per-hit contribution + // Find peak amplitude and bin + Double_t peakAmp = 0.0; + Int_t peakBin = 0; + for(Int_t i = 0; i < fDaq.num_samples; i++) { + if(hit_contrib[i] > peakAmp) { + peakAmp = hit_contrib[i]; + peakBin = i; + } + } + + // Integral: sum over non-negligible samples (>1% of peak), + // converted to ADC counts (pedestal-subtracted) + Double_t integral = 0.0; + Double_t threshold = peakAmp * 0.01; + for(Int_t i = 0; i < fDaq.num_samples; i++) { + if(hit_contrib[i] > threshold) { + integral += hit_contrib[i]; + } + } + integral /= fDaq.v_lsb; // mV -> ADC counts + + // Amplitude: convert to ADC with pedestal (matches fSampPulseAmp) + Double_t ampADC = peakAmp / fDaq.v_lsb + fDaq.pedestal; + + // Time: half-max crossing on rising edge, in 1/64th-sample units + Double_t vMid = peakAmp / 2.0; + Double_t timeVal = peakBin * 64.0; // fallback: peak bin + for(Int_t i = 0; i < peakBin; i++) { + if(hit_contrib[i] <= vMid && hit_contrib[i + 1] > vMid) { + timeVal = 64.0 * i + 64.0 * (vMid - hit_contrib[i]) / + (hit_contrib[i + 1] - hit_contrib[i]); + break; + } + } + + fTruthInfo.pulseTime.push_back(timeVal); + fTruthInfo.pulseAmplitude.push_back(ampADC); + fTruthInfo.pulseIntegral.push_back(integral); + } + + // Digitize + std::vector adc_wf(fDaq.num_samples, 0); + for(Int_t i = 0; i < fDaq.num_samples; i++) { + Int_t adc_count = static_cast(std::floor(triangle_wf[i]/fDaq.v_lsb) + fDaq.pedestal); + adc_count = std::max(0, std::min(adc_count, fDaq.max_adc)); + adc_wf[i] = adc_count; + } + return adc_wf; +} \ No newline at end of file diff --git a/tests/waveformGenerator.h b/tests/waveformGenerator.h index 7a7411f1..6fde782a 100644 --- a/tests/waveformGenerator.h +++ b/tests/waveformGenerator.h @@ -1,69 +1,99 @@ #ifndef WAVEFORMGENERATOR_H #define WAVEFORMGENERATOR_H +#include +#include #include "Rtypes.h" #include "RtypesCore.h" #include "TMath.h" -#include -#include -#include "../Podd/FadcStandaloneAnalyzer.h" +struct DigitizerParams { + Double_t sample_rate = 4.0; // Sampling rate of the digitizer (250 MHz) + Int_t num_samples = 256; // #of samples in the sampling window + Double_t pedestal = 100.0; // Pedestal in ADC counts + Double_t noise_sigma = 1.5; // Electronic noise (mV) + Double_t v_lsb = 1.0; // Voltage resolution (mV per ADC count) + Int_t max_adc = 4095; // Maximum ADC value (12-bit) + Int_t max_n_pulses = 4; // Max # of pulses per window -struct WaveformData { - Int_t TotNSamples; - Int_t InitPedSamples; + // Params needed for the initial algorithm and constant shape pulses + Int_t init_ped_samples = 4; // # of samples needed to calculate initial pedestal + Int_t ped_sigma = 4; // Pedestal variation in number of samples +}; + +// Structs for the waveform generator algorithm v2 - keep the safety values +struct WaveformParams { + Double_t tau_r = 2.0; // Rise time (ns) + Double_t tau_f = 5.0; // Decay time (ns) + Double_t tts_sigma = 1.5; // Transit time spread (ns) + Double_t mean_gain = 10.0; // Mean amplitude per PE (mV) + Double_t gain_sigma = 1.0; // Gain variation (mV) - Int_t Pedestal; - Int_t PedRange; + // After pulse parameters + Double_t ap_prob = 0.03; // After pulse generation probability + Double_t ap_delay_mean = 80.0; // Mean delay for afterpulse (ns) + Double_t ap_delay_sigma = 5.0; // Spread in afterpulse delay (ns) - Int_t NPulses; - Int_t PulseAmp; - Int_t PulseWidth; - Int_t PulseSeparation; + // Parameters needed for the initial algorithm + Int_t pulse_width = 3; // Pulse width in samples + Int_t pulse_separation = 10; // Pulse separation in samples }; -struct WaveformInfo { - std::vector - PedSubIntegral; // Double to match analyzer's floating-point formula - std::vector PulseAmp; - std::vector PulseTime; +// Truth-level pulse info per hit, computed during waveform generation. +// One entry per hit group (not per PE). Values are converted to match +// the units of FADCStandaloneAnalyzer output for direct comparison. +struct TruthPulseInfo { + std::vector pulseIntegral; // Ped-subtracted integral in ADC counts (÷ v_lsb) + std::vector pulseAmplitude; // Peak in ADC counts (÷ v_lsb + pedestal) + std::vector pulseTime; // Half-max crossing in 1/64th-sample units }; class waveformGenerator { -public: - waveformGenerator(); // Seed the random number generator. New seed everytime - // run. - virtual ~waveformGenerator() = default; - - // Waveform Data Setters - void setWaveformData(const WaveformData &waveformData) { - fWaveformData = waveformData; - } - void setAnalyzer(const FadcStandaloneAnalyzer &analyzer) { - saAnalyzer = analyzer; - } - - // Waveform generators - std::vector noiseWaveform() const; - std::vector singlePulseWaveform() const; - std::vector multiplePulseWaveform() const; - - // Calc pulse parameters - // std::vector calcPulseIntegral(const std::vector wave) const; - void calcPulseParameters(const std::vector wave); - - // Accessors for pulse parameters computed by calcPulseParameters() - const WaveformInfo &getWaveformInfo() const { return fWaveformInfo; } - void clearWaveformInfo() { fWaveformInfo = {}; } - -private: - WaveformData fWaveformData; - WaveformInfo fWaveformInfo; - - FadcStandaloneAnalyzer saAnalyzer; - - // Per-instance RNG — thread-safe (no shared global state) - mutable std::mt19937 fRng; + public: + waveformGenerator(); // Seed the random number generator. New seed everytime run. + virtual ~waveformGenerator() = default; + + // NEW: Detector + Digitizer Setter + void setDetectorParams(const WaveformParams& params) {fWf = params;} + void setDigitizerParams(const DigitizerParams& params) {fDaq = params;} + + // Waveform generators + // Preliminary generators + std::vector noiseWaveform() const; + std::vector singlePulseWaveform() const; + std::vector multiplePulseWaveform() const; + + // PE based generators + std::vector generateWaveform(); + std::vector squareWaveform(); + std::vector triangleWaveform(); + + // Accessors for truth-level pulse info computed during generation + const TruthPulseInfo& getTruthInfo() const { return fTruthInfo; } + void clearTruthInfo() { fTruthInfo = {}; } + + private: + TruthPulseInfo fTruthInfo; + + // Members for waveform generator v2 + DigitizerParams fDaq; + WaveformParams fWf; + + // Per-instance RNG — thread-safe (no shared global state) + mutable std::mt19937 fRng; + + // Helpers for waveform generator v2 + // Waveform templates + Double_t spePulse(Double_t t, Double_t t0, Double_t amp) const; + Double_t squarePulse(Double_t t, Double_t t0, Double_t amp) const; + Double_t trianglePulse(Double_t t, Double_t t0, Double_t amp) const; + + // Hit-level PE grouping for structured truth extraction + struct PEHit { + Double_t hitTime; // Primary hit time (ns) + std::vector peTimes; // All PE times in this hit (incl. afterpulses) + }; + std::vector generatePEHits() const; }; #endif // WAVEFORMGENERATOR_H \ No newline at end of file