diff --git a/duneopdet/OpticalDetector/WaveformPreProcessing.fcl b/duneopdet/OpticalDetector/WaveformPreProcessing.fcl index 961c8d26..7ba8202d 100644 --- a/duneopdet/OpticalDetector/WaveformPreProcessing.fcl +++ b/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 diff --git a/duneopdet/OpticalDetector/WaveformPreProcessing_module.cc b/duneopdet/OpticalDetector/WaveformPreProcessing_module.cc index 92b8e81e..8bf4cf34 100644 --- a/duneopdet/OpticalDetector/WaveformPreProcessing_module.cc +++ b/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 wf, Float_t fDynamicRangeSaturation, size_t fMaxTicksSat); void BaselineExtractor(std::vector& wf); - void DenoisingAlgo(std::vector& waveform, double lambda); + void DentCorrection(std::vector &wf, Float_t fDynamicRangeSaturation, size_t MaxTicksDent, int channel, double MaxDentDepth); + void DenoisingAlgo(std::vector& waveform, double lambda, Float_t fDynamicRangeSaturation); bool TV1D_denoise(std::vector& waveform, std::vector& outwaveform, const double lambda); std::vector ComputeMovingAverage(const std::vector& data, int n); double median_func(std::vector 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 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("InputModule"); fApplyDenoising = p.get("ApplyDenoising"); - fIsPDVD = p.get("IsPDVD"); + fIgnoreChannels = p.get >("IgnoreChannels"); fLambda = p.get("Lambda"); //parameter for the denoising algorithm fMaxTicks = p.get("MaxTicks"); //maximum number of ticks in the waveform + fMaxTicksDent = p.get("MaxTicksDent"); //maximum number of ticks for dent formation + fMaxTicksSat = p.get("MaxTicksSat"); //maximum number of saturated ticks in the waveform to apply baseline removal fSecondBaselineSub = p.get("SecondBaselineSub"); //the mode of lowest SecondBaslineSub of the signal for second baseline estimate fFirstBaselineSub = p.get("FirstBaselineSub"); //around 3 times the expected large signals + fMaxDentDepth = p.get("MaxDentDepth"); //Maximum percentage of falling from saturation to be considered a laser failure + fDynamicRangeSaturation = p.get("DynamicRangeSaturation"); + fRemoveBaselineFluctuation = p.get("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 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 &wf){ std::vector signal_base, waveform_full_bs; std::vector temp = wf; - //MA, baseline estimate, basline subtraction + //MA, baseline estimate, baseline subtraction std::vector signalsma = ComputeMovingAverage(temp, 4); std::vector basev = EstimateBaselineOpeningCentered(signalsma, fFirstBaselineSub, 1); for (size_t j=0; j &waveform, double lambda){ + void WaveformPreProcessing::DentCorrection(std::vector &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]saturation + } + + bool WaveformPreProcessing::CheckSaturation(std::vector 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 &waveform, double lambda, Float_t fDynamicRangeSaturation){ int length = waveform.size(); std::vector 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] WaveformPreProcessing::VectorOfDoublesToVectorOfShorts (std::vector const& vectorOfDoubles) { - // Don't bother to round properly, it's faster this way - return std::vector(vectorOfDoubles.begin(), vectorOfDoubles.end()); + // Rounding properly + std::vector result; + result.reserve(vectorOfDoubles.size()); + for (double val : vectorOfDoubles) { + result.push_back(static_cast(std::round(val))); + } + return result; } }//namespace opdet