#include "PedAlgoRollingMean.h"

Inheritance diagram for pmtana::PedAlgoRollingMean:

Public Member Functions
	PedAlgoRollingMean (const std::string name="PedRollingMean")
	Default constructor. More...

	PedAlgoRollingMean (const fhicl::ParameterSet &pset, const std::string name="PedRollingMean")
	Alternative ctor. More...

virtual	~PedAlgoRollingMean ()
	Default destructor. More...

const std::string &	Name () const
	Name getter. More...

bool	Evaluate (const pmtana::Waveform_t &wf)
	Method to compute a pedestal. More...

double	Mean (size_t i) const
	Getter of the pedestal mean value. More...

const pmtana::PedestalMean_t &	Mean () const
	Getter of the pedestal mean value. More...

double	Sigma (size_t i) const
	Getter of the pedestal standard deviation. More...

const pmtana::PedestalSigma_t &	Sigma () const
	Getter of the pedestal standard deviation. More...

Protected Member Functions
bool	ComputePedestal (const pmtana::Waveform_t &wf, pmtana::PedestalMean_t &mean_v, pmtana::PedestalSigma_t &sigma_v)
	Method to compute a pedestal of the input waveform using "nsample" ADC samples from "start" index. More...

virtual bool	ComputePedestal (const ::pmtana::Waveform_t &wf, pmtana::PedestalMean_t &mean_v, pmtana::PedestalSigma_t &sigma_v)=0

Private Attributes
size_t	_sample_size

float	_max_sigma

float	_ped_range_max

float	_ped_range_min

double	_threshold

double	_diff_threshold

double	_diff_adc_count

int	_n_presamples

Detailed Description

A class that calculates pedestal mean & standard deviation (here and elsewhere called as "RMS").

Definition at line 30 of file PedAlgoRollingMean.h.

Constructor & Destructor Documentation

pmtana::PedAlgoRollingMean::PedAlgoRollingMean ( const std::string name = "PedRollingMean" )

Default constructor.

Definition at line 18 of file PedAlgoRollingMean.cxx.

     : PMTPedestalBase(name)
       //*****************************************************************
   {
     srand(static_cast<unsigned int>(time(0)));
   }

pmtana::PedAlgoRollingMean::PedAlgoRollingMean	(	const fhicl::ParameterSet &	pset,
		const std::string	name = `"PedRollingMean"`
	)

Alternative ctor.

Definition at line 26 of file PedAlgoRollingMean.cxx.

References _diff_adc_count, _diff_threshold, _max_sigma, _n_presamples, _ped_range_max, _ped_range_min, _sample_size, _threshold, and fhicl::ParameterSet::get().

     : PMTPedestalBase(name)
       //############################################################
   {
 
     _sample_size    = pset.get<size_t>("SampleSize");
     _max_sigma      = pset.get<float> ("MaxSigma");
     _ped_range_max  = pset.get<float> ("PedRangeMax");
     _ped_range_min  = pset.get<float> ("PedRangeMin");
 
     // _range          = pset.get<int>   ("RandomRange");
     // _divisions      = pset.get<double>("RandomRangeDivisions");
     _threshold      = pset.get<double>("Threshold");
     _diff_threshold = pset.get<double>("DiffBetweenGapsThreshold");
     _diff_adc_count = pset.get<double>("DiffADCCounts");
 
     _n_presamples   = pset.get<int>("NPrePostSamples");
     
     //_random_shift   = pset.get<double>("RandomRangeShift");
     // Random seed number generator
     //srand(static_cast<unsigned int>(time(0)));
   }

pmtana::PedAlgoRollingMean::~PedAlgoRollingMean ( )

virtual

Default destructor.

Definition at line 52 of file PedAlgoRollingMean.cxx.

54 {}

Member Function Documentation

bool pmtana::PedAlgoRollingMean::ComputePedestal	(	const pmtana::Waveform_t &	wf,
		pmtana::PedestalMean_t &	mean_v,
		pmtana::PedestalSigma_t &	sigma_v
	)

protected

Method to compute a pedestal of the input waveform using "nsample" ADC samples from "start" index.

Definition at line 57 of file PedAlgoRollingMean.cxx.

References _diff_adc_count, _diff_threshold, _max_sigma, _n_presamples, _ped_range_max, _sample_size, _threshold, pmtana::BinnedMaxOccurrence(), pmtana::edge_aware_mean(), pmtana::mean(), and pmtana::std().

   {
 
     // parameters
     if(wf.size()<=(_sample_size * 2))
       return false;
 
     const size_t window_size = _sample_size*2;
 
     // middle mean
     for(size_t i=0; i< wf.size(); ++i) {
 
       mean_v[i]  = 0;
       sigma_v[i] = 0;
 
       if( i < _sample_size || i >= (wf.size() - _sample_size) ) continue;
 
       mean_v[i]  = mean (wf,i - _sample_size,window_size);
       sigma_v[i] = std  (wf,mean_v[i],i - _sample_size,window_size);
     }
 
     // front mean
     for(size_t i=0; i<_sample_size; ++i) {
 
       mean_v[i]  = mean_v [_sample_size];
       sigma_v[i] = sigma_v[_sample_size];
       
     }
     // tail mean
     for(size_t i=(wf.size() - _sample_size); i<wf.size(); ++i) {
 
       mean_v[i]  = mean_v [wf.size() - _sample_size -1];
       sigma_v[i] = sigma_v[wf.size() - _sample_size -1];
 
     }
     
     float  best_sigma = 1.1e9;
     size_t best_sigma_index = 0;
     size_t num_good_adc = 0;
 
     for(size_t i=0; i<sigma_v.size(); ++i) {
       // Only consider adcs which mean is in the allowed range
       auto const& mean  = mean_v[i];
       
       if( mean < _ped_range_min || mean > _ped_range_max ) continue;
 
       auto const& sigma = sigma_v[i];
       if(sigma < best_sigma) {
         best_sigma = sigma;
         best_sigma_index = i;
       }
 
       if(sigma < _max_sigma) num_good_adc += 1;
     }
 
 
     if( num_good_adc < 1 ) {
       std::cerr << "\033[93m<<" << __FUNCTION__ << ">>\033[00m Could not find good pedestal at all..." << std::endl;
       return false;
     }
 
     // If not enough # of good mean indices, use the best guess within this waveform
     if(best_sigma > _max_sigma || num_good_adc < 3) {
       for(size_t i=0; i<mean_v.size(); ++i) {
         mean_v[i]  = mean_v.at  ( best_sigma_index );
         sigma_v[i] = sigma_v.at ( best_sigma_index );
       }
 
       return true;
     }
 
     
     // Else do extrapolation, or random seed depending on what we find...
 
     unsigned nbins = 1000;  
 
     const auto mode_mean  = BinnedMaxOccurrence(mean_v ,nbins);
     const auto mode_sigma = BinnedMaxOccurrence(sigma_v,nbins);
 
     //auto mode_mean  = BinnedMaxTH1D(mean_v ,nbins);
     //auto mode_sigma = BinnedMaxTH1D(sigma_v,nbins);
 
     //std::cout<<mode_mean<<" +/- "<<mode_sigma<<std::endl;
 
     _diff_threshold *= mode_sigma;
     
     double diff_cutoff = _diff_threshold < _diff_adc_count ? _diff_adc_count : _diff_threshold;
 
     int last_good_index = -1;
 
     for(size_t i=0; i < wf.size(); ++i) {
       
       auto const mean  = mean_v[i];
       auto const sigma = sigma_v[i];
 
       // if(sigma <= _max_sigma && mean < _ped_range_max && mean > _ped_range_min) {
       // not sure if this works well for basline that is "linear" seen by David K
       
       if(sigma <= _threshold * mode_sigma && fabs(mean - mode_mean) <= _threshold * mode_sigma) {
 
         if(last_good_index<0) {
           last_good_index = (int)i;
           continue;
         }
 
         if( ( last_good_index + 1 ) < (int) i ) {
 
           auto diff = fabs(mean_v.at(last_good_index) - mean);
 
           if ( diff >  diff_cutoff) {
             //this should become generic interpolation function, for now lets leave.
             float slope = (mean - mean_v.at(last_good_index)) / (float(i - last_good_index));
             
             for(size_t j = last_good_index + 1; j < i; ++j) {
               mean_v.at(j)  = slope * ( float(j - last_good_index) ) + mean_v.at(last_good_index);
               sigma_v.at(j) = mode_sigma;
             }
           }
           else {
             //difference roughly the same lets fill with constant baseline
             
             auto presample_mean  = edge_aware_mean(wf,last_good_index - _n_presamples, last_good_index);
             auto postsample_mean = edge_aware_mean(wf,i, _n_presamples);
             
             auto diff_pre  = fabs(presample_mean  - mode_mean);
             auto diff_post = fabs(postsample_mean - mode_mean);
 
             auto constant_mean = diff_pre <= diff_post ? presample_mean : postsample_mean;
             
             for(size_t j = last_good_index + 1; j < i; ++j) {
               //mean_v.at(j)  = floor( mean_v.at(last_good_index) ) + _random_shift + (double) ( rand() % _range) / _divisions;
               mean_v.at(j)  = constant_mean;
               sigma_v.at(j) = mode_sigma;
             }
           }
         }
         last_good_index = (int)i;
       }
     }
     
 
     // Next do extrapolation to the first and end (if needed)
     // vic: for now we leave this i'm not sure if this really needs
     //      to be tuned until I can make unit test
     // update: yes this needs work...
 
     if(sigma_v.front() > mode_sigma) {
 
       int first_index  = -1;
       int second_index = -1;
 
       for(size_t i=0; i < wf.size(); ++i) {
         if( sigma_v.at(i) < mode_sigma ) {
           if( first_index < 0 ) first_index = (int)i;
           else if( second_index < 0 ) {
             second_index = (int)i;
             break;
           }
         }
       }
       
       if(first_index < 0 || second_index < 0) {
         std::cerr <<"\033[93m<<" << __FUNCTION__ << ">>\033[00m Could not find good pedestal for CDF"
                   << "\n"
                   << "first_index:  " << first_index << "\n"
                   << "second_index: " << second_index << "\n"
                   << "If one of these is less than zero, this means there is pulse\n"
                   << "on first sample and baseline never went back down. Returning false here.";
         return false;
       }
 
 
       auto diff = fabs(mean_v.at(first_index) - mean_v.at(second_index));
       
       if ( diff >  diff_cutoff) {
         
         float slope = (mean_v.at(second_index) - mean_v.at(first_index)) / (float(second_index - first_index));
         
         for(int j=0; j < first_index; ++j) {
           mean_v.at(j)  = mean_v.at(first_index) - slope * (first_index - j);
           sigma_v.at(j) = _max_sigma;
         }
         
       } else {
 
         auto postsample_mean = edge_aware_mean(wf,first_index, first_index + _n_presamples);
 
         for(int j=0; j < first_index; ++j) {
           //mean_v.at(j)  = floor( mean_v.at(second_index) ) + _random_shift + (double) ( rand() % _range) / _divisions;
           mean_v.at(j)  = postsample_mean;
           sigma_v.at(j) = mode_sigma;
         }
       }
       
     }
     
     
     if(sigma_v.back() > mode_sigma) {
 
       int first_index  = -1;
       int second_index = -1;
       
       for(int i = wf.size()-1; i >= 0; --i) {
         if(sigma_v.at(i) < mode_sigma) {
           if( second_index < 0 ) second_index = (int)i;
           else if( first_index < 0 ) {
             first_index = (int)i;
             break;
           }
         }
       }
       
       if(first_index < 0 || second_index < 0) {
         std::cerr <<"\033[93m<<" << __FUNCTION__ << ">>\033[00m Could not find good pedestal for CDF"
                   << "\n"
                   << "first_index:  " << first_index << "\n"
                   << "second_index: " << second_index << "\n"
                   << "If one of these is less than zero, this means there is pulse\n"
                   << "on the last sample and baseline never went back down. Returning false here.";
         return false;
       }
             
       
       auto diff = fabs(mean_v.at(second_index) - mean_v.at(first_index) );
       
       if ( diff >  diff_cutoff) {
     
         float slope = (mean_v.at(second_index) - mean_v.at(first_index)) / (float(second_index - first_index));
         for(int j = second_index+1; j < int(wf.size()); ++j) {
           mean_v.at(j)  = mean_v.at(second_index) + slope * (j-second_index);
           sigma_v.at(j) = _max_sigma;
         }
 
       }
       else {
         
         auto presample_mean  = edge_aware_mean(wf,first_index - _n_presamples, second_index);
         
         for(int j = second_index+1; j < int(wf.size()); ++j) {
           //mean_v.at(j)  = floor( mean_v.at(first_index) ) + _random_shift + (double) ( rand() % _range) / _divisions;
           mean_v.at(j)  = presample_mean;
           sigma_v.at(j) = mode_sigma;
         }
       }
     }
 
     
     return true;
 
   }

virtual bool pmtana::PMTPedestalBase::ComputePedestal	(	const ::pmtana::Waveform_t &	wf,
		pmtana::PedestalMean_t &	mean_v,
		pmtana::PedestalSigma_t &	sigma_v
	)

protectedpure virtualinherited

Method to compute pedestal: mean and sigma array should be filled per ADC. The length of each array is guaranteed to be same.

Referenced by pmtana::PMTPedestalBase::Evaluate().

bool pmtana::PMTPedestalBase::Evaluate ( const pmtana::Waveform_t & wf )

inherited

Method to compute a pedestal.

Definition at line 33 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_mean_v, pmtana::PMTPedestalBase::_sigma_v, and pmtana::PMTPedestalBase::ComputePedestal().

Referenced by pmtana::PedAlgoUB::ComputePedestal(), and pmtana::PulseRecoManager::Reconstruct().

   {
     _mean_v.resize(wf.size(),0);
     _sigma_v.resize(wf.size(),0);
 
     for(size_t i=0; i<wf.size(); ++i)
       _mean_v[i] = _sigma_v[i] = 0;
 
     const bool res = ComputePedestal(wf, _mean_v, _sigma_v);
 
     if(wf.size() != _mean_v.size())
       throw OpticalRecoException("Internal error: computed pedestal mean array length changed!");
     if(wf.size() != _sigma_v.size())
       throw OpticalRecoException("Internal error: computed pedestal sigma array length changed!");
 
     return res;
   }

double pmtana::PMTPedestalBase::Mean ( size_t i ) const

inherited

Getter of the pedestal mean value.

Definition at line 53 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_mean_v, and ss.

Referenced by pmtana::PedAlgoUB::ComputePedestal(), and pmtana::PulseRecoManager::Reconstruct().

   {
     if(i > _mean_v.size()) {
       std::stringstream ss;
       ss << "Invalid index: no pedestal mean exist @ " << i;
       throw OpticalRecoException(ss.str());
     }
     return _mean_v[i];
   }

const PedestalMean_t & pmtana::PMTPedestalBase::Mean ( ) const

inherited

Getter of the pedestal mean value.

Definition at line 77 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_mean_v.

79 { return _mean_v; }

pmtana::PMTPedestalBase::_mean_v

pmtana::PedestalMean_t _mean_v

A variable holder for pedestal mean value.

Definition: PMTPedestalBase.h:72

const std::string & pmtana::PMTPedestalBase::Name ( ) const

inherited

Name getter.

Definition at line 28 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_name.

30 { return _name;}

pmtana::PMTPedestalBase::_name

std::string _name

Name.

Definition: PMTPedestalBase.h:69

double pmtana::PMTPedestalBase::Sigma ( size_t i ) const

inherited

Getter of the pedestal standard deviation.

Definition at line 65 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_sigma_v, and ss.

Referenced by pmtana::PedAlgoUB::ComputePedestal(), and pmtana::PulseRecoManager::Reconstruct().

   {
     if(i > _sigma_v.size()) {
       std::stringstream ss;
       ss << "Invalid index: no pedestal sigma exist @ " << i;
       throw OpticalRecoException(ss.str());
     }
     return _sigma_v[i];
   }

const PedestalSigma_t & pmtana::PMTPedestalBase::Sigma ( ) const

inherited

Getter of the pedestal standard deviation.

Definition at line 82 of file PMTPedestalBase.cxx.

References pmtana::PMTPedestalBase::_sigma_v.

84 { return _sigma_v; }

pmtana::PMTPedestalBase::_sigma_v

pmtana::PedestalSigma_t _sigma_v

A variable holder for pedestal standard deviation.

Definition: PMTPedestalBase.h:75

Member Data Documentation

double pmtana::PedAlgoRollingMean::_diff_adc_count

private

Definition at line 62 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

double pmtana::PedAlgoRollingMean::_diff_threshold

private

Definition at line 61 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

float pmtana::PedAlgoRollingMean::_max_sigma

private

Definition at line 54 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

int pmtana::PedAlgoRollingMean::_n_presamples

private

Definition at line 64 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

float pmtana::PedAlgoRollingMean::_ped_range_max

private

Definition at line 55 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

float pmtana::PedAlgoRollingMean::_ped_range_min

private

Definition at line 56 of file PedAlgoRollingMean.h.

Referenced by PedAlgoRollingMean().

size_t pmtana::PedAlgoRollingMean::_sample_size

private

Definition at line 53 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

double pmtana::PedAlgoRollingMean::_threshold

private

Definition at line 60 of file PedAlgoRollingMean.h.

Referenced by ComputePedestal(), and PedAlgoRollingMean().

The documentation for this class was generated from the following files:

larana/v06_21_01/source/larana/OpticalDetector/OpHitFinder/PedAlgoRollingMean.h
larana/v06_21_01/source/larana/OpticalDetector/OpHitFinder/PedAlgoRollingMean.cxx

Public Member Functions

Protected Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation