pmtana::PedAlgoRmsSlider Class Reference

#include <PedAlgoRmsSlider.h>

Inheritance diagram for pmtana::PedAlgoRmsSlider:

Public Member Functions
	PedAlgoRmsSlider (const std::string name="PedRmsSlider")
	Default constructor. More...
	PedAlgoRmsSlider (const fhicl::ParameterSet &pset, const std::string name="PedRmsSlider")
	Alternative ctor. More...
void	PrintInfo ()
	Print settings. More...
Public Member Functions inherited from pmtana::PMTPedestalBase
	PMTPedestalBase (std::string name="noname")
	Default constructor. More...
virtual	~PMTPedestalBase ()
	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...
double	Sigma (size_t i) const
	Getter of the pedestal standard deviation. More...
const pmtana::PedestalMean_t &	Mean () const
	Getter of the pedestal mean value. 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...
Protected Member Functions inherited from pmtana::PMTPedestalBase
virtual bool	ComputePedestal (const ::pmtana::Waveform_t &wf, pmtana::PedestalMean_t &mean_v, pmtana::PedestalSigma_t &sigma_v)=0

Private Member Functions
double	CalcMean (const std::vector< double > &wf, size_t start, size_t nsample)
	Returns the mean of the elements of the vector from start to start+nsample. More...
double	CalcStd (const std::vector< double > &wf, const double ped_mean, size_t start, size_t nsample)
	Returns the std of the elements of the vector from start to start+nsample. More...
bool	CheckSanity (pmtana::PedestalMean_t &mean_v, pmtana::PedestalSigma_t &sigma_v)
	Checks the sanity of the estimated pedestal, returns false if not sane. More...

Private Attributes
size_t	_sample_size
	How many samples are used to calculate local rms and mean. More...
double	_threshold
	Threshold applied to local rms to claim a pulse. More...
float	_max_sigma
	Max sigma to consider adc as 'sane'. More...
float	_ped_range_max
	Max value of adc to consider adc as 'sane'. More...
float	_ped_range_min
	Min value of adc to consider adc as 'sane'. More...
bool	_verbose
	For debugging. More...
int	_n_wf_to_csvfile
	If greater than zero saves firsts waveforms with pedestal to csv file. More...
int	_wf_saved = 0
int	_num_presample
	number of ADCs to sample before the gap More...
int	_num_postsample
	number of ADCs to sample after the gap More...
std::ofstream	_csvfile

Detailed Description

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

Definition at line 32 of file PedAlgoRmsSlider.h.

Constructor & Destructor Documentation

pmtana::PedAlgoRmsSlider::PedAlgoRmsSlider

(

const std::string

name = "PedRmsSlider"

)

Default constructor.

Definition at line 19 of file PedAlgoRmsSlider.cxx.

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

pmtana::PedAlgoRmsSlider::PedAlgoRmsSlider	(	const fhicl::ParameterSet &	pset,
		const std::string	name = "PedRmsSlider"
	)

Alternative ctor.

Definition at line 27 of file PedAlgoRmsSlider.cxx.

    : PMTPedestalBase(name)
      //############################################################
  {

    _sample_size     = pset.get<size_t>("SampleSize",       7    );
    _threshold       = pset.get<double>("Threshold",        0.6  );
    _max_sigma       = pset.get<float> ("MaxSigma",         0.5  );
    _ped_range_max   = pset.get<float> ("PedRangeMax",      2150 );
    _ped_range_min   = pset.get<float> ("PedRangeMin",      100  );
    _num_presample   = pset.get<int>   ("NumPreSample",     0    );
    _num_postsample  = pset.get<int>   ("NumPostSample",    0    );
    _verbose         = pset.get<bool>  ("Verbose",          true );
    _n_wf_to_csvfile = pset.get<int>   ("NWaveformsToFile", 12   );

    if (_n_wf_to_csvfile > 0) {
      _csvfile.open ("wf_pedalgormsslider.csv", std::ofstream::out | std::ofstream::trunc);
      _csvfile << "n,time,wf,wf_ped_mean,wf_ped_rms" << std::endl;
    }

  }

Member Function Documentation

double pmtana::PedAlgoRmsSlider::CalcMean ( const std::vector< double > &  wf, size_t  start, size_t  nsample  )

private

Returns the mean of the elements of the vector from start to start+nsample.

Definition at line 62 of file PedAlgoRmsSlider.cxx.

  {
    if(!nsample) nsample = wf.size();
    if(start > wf.size() || (start+nsample) > wf.size())
      throw OpticalRecoException("Invalid start/end index!");

    double sum = std::accumulate(wf.begin()+start,wf.begin()+start+nsample,0.0);

    sum /= ((double)nsample);

    return sum;
  }

double pmtana::PedAlgoRmsSlider::CalcStd ( const std::vector< double > &  wf, const double  ped_mean, size_t  start, size_t  nsample  )

private

Returns the std of the elements of the vector from start to start+nsample.

Definition at line 77 of file PedAlgoRmsSlider.cxx.

  {
    if(!nsample) nsample = wf.size();
    if(start > wf.size() || (start+nsample) > wf.size())
      throw OpticalRecoException("Invalid start/end index!");

    double sigma = 0;

    for(size_t index=start; index < (start+nsample); ++index){
      sigma += pow( (wf[index] - ped_mean), 2 );
    }

    sigma = sqrt(sigma/((double)(nsample)));

    return sigma;
  }

bool pmtana::PedAlgoRmsSlider::CheckSanity ( pmtana::PedestalMean_t &  mean_v, pmtana::PedestalSigma_t &  sigma_v  )

private

Checks the sanity of the estimated pedestal, returns false if not sane.

Definition at line 411 of file PedAlgoRmsSlider.cxx.

  {

    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) {

       if(_verbose) {
         std::cout << "\033[93mPedAlgoRmsSlider\033[00m: Not enough number of good mean indices."
           << "Using the best guess within this waveform."
           << std::endl;
       }

      for(size_t i=0; i<mean_v.size(); ++i) {
        mean_v[i]  = mean_v  [ best_sigma_index ];
        sigma_v[i] = sigma_v [ best_sigma_index ];
      }
    }

    return true;

  }

bool pmtana::PedAlgoRmsSlider::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 96 of file PedAlgoRmsSlider.cxx.

  {

    if (_verbose)
      this->PrintInfo();

    if (wf.size() <= (_sample_size * 2))
      return false;

    // Prepare output
    mean_v.resize (wf.size(), 0);
    sigma_v.resize(wf.size(), 0);




    // **********
    // To start, set the pedestal equal to
    // the wf itself
    // **********

    pmtana::PedestalMean_t  mean_temp_v;
    mean_temp_v.resize( wf.size(), 0);

    for(size_t i=0; i< wf.size(); ++i) {
      mean_temp_v[i]  = wf[i];
      sigma_v[i] = 0;
    }





    // **********
    // Now look for rms variations
    // and change the mean and rms accordingly
    // **********
    int last_good_index = -1;
    double local_mean, local_rms;
    std::vector<double> local_mean_v(wf.size(),-1.);
    std::vector<double> local_sigma_v(wf.size(),-1.);

    for (size_t i = 0; i < wf.size() - _sample_size; i++) {

      local_mean = mean(wf, i, _sample_size);
      local_rms  = std(wf, local_mean, i, _sample_size);

      if(_verbose) std::cout << "\033[93mPedAlgoRmsSlider\033[00m: i " << i 
                             << "  local_mean: " << local_mean 
                             << "  local_rms: " << local_rms << std::endl;

      if (local_rms < _threshold) {

        local_mean_v[i] = local_mean;
        local_sigma_v[i] = local_rms;

        if(_verbose)
          std::cout << "\033[93mBelow threshold\033[00m: "
                    << "at i " << i
                    << " last good index was: " << last_good_index
                    << std::endl;

      }
    }

    // find the gaps (regions to be interpolated
    last_good_index = -1;
    std::vector<bool> ped_interapolated(wf.size(),false);
    for(size_t i=0; i < wf.size() - _sample_size; i++) {

      if(local_mean_v[i] > -0.1) {
        // good pedestal!

        if( ( last_good_index + 1 ) < (int)i ) {
          // finished the gap. try interpolation
          // 0) find where to start/end interpolation
          int start_tick  = last_good_index;
          int end_tick    = i;
          int start_bound = std::max(last_good_index - _num_presample, 0);
          int end_bound   = std::min(i + _num_postsample, (int)(wf.size()) - _sample_size);
          for(int j=start_tick; j>=start_bound; --j) {
            if(local_mean_v[j] < 0) continue;
            start_tick = j;
          }
          for(int j=end_tick; j<=end_bound; ++j) {
            if(local_mean_v[j] < 0) continue;
            end_tick = j;
          }

          //this should become generic interpolation function, for now lets leave.
          float slope = (local_mean_v[end_tick] - local_mean_v[start_tick]) / (float(end_tick - start_tick));

          for(int j = start_tick + 1; j < end_tick; ++j) {
            mean_temp_v[j]  = slope * ( float(j - start_tick) ) + local_mean_v[start_tick];
            // for sigma, put either the sigma in the region before the pulse or
            // after the pulse, depending on which one if != 0. If both are !=0 put the one after
            // the pulse (just default), if both are zero then put zero
            sigma_v[j] = (local_sigma_v[end_tick] != 0 ? local_sigma_v[end_tick] : local_sigma_v[start_tick]); // todo: fluctuate baseline
            ped_interapolated[j] = true;
          }
        }

        last_good_index = i;
      }
    }

    /*

    for (size_t i = 0; i < wf.size() - _sample_size; i++) {

      local_mean = mean(wf, i, _sample_size);
      local_rms  = std(wf, local_mean, i, _sample_size);

      if(_verbose) std::cout << "\033[93mPedAlgoRmsSlider\033[00m: i " << i << "  local_mean: " << local_mean << "  local_rms: " << local_rms << std::endl;

      if (local_rms < _threshold) {

        local_mean_v[i] = local_mean;
        local_sigma_v[i] = local_sigma;

        if(_verbose)
          std::cout << "\033[93mBelow threshold\033[00m: "
                    << "at i " << i
                    << " last good index was: " << last_good_index
                    << std::endl;

        if(last_good_index<0) {
          last_good_index = (int)i;
          last_local_mean = local_mean;
          last_local_rms  = local_rms;
          continue;
        }


        if( ( last_good_index + 1 ) < (int)i ) {

          //this should become generic interpolation function, for now lets leave.
          float slope = (local_mean - last_local_mean) / (float(i - last_good_index));

          for(size_t j = last_good_index + 1; j < i && j < wf.size(); ++j) {
            mean_temp_v.at(j)  = slope * ( float(j - last_good_index) ) + mean_temp_v.at(last_good_index);
            // for sigma, put either the sigma in the region before the pulse or
            // after the pulse, depending on which one if != 0. If both are !=0 put the one after
            // the pulse (just default), if both are zero then put zero
            sigma_v.at(j) = (local_rms != 0 ? local_rms : last_local_rms); // todo: fluctuate baseline
            ped_interapolated.at(j) = true;
          }
        }

        // record this mean & rms as good mean value
        last_good_index = i;
        last_local_mean = local_mean;
        last_local_rms  = local_rms;
        // if _num_postsample is specified, go back in time to look for it
        if(_num_postsample >0 && (i>_num_postsample)) {
          int loop_start = std::max(((int)i) - _num_postsample, 0);
          for(int j=loop_start; j>=0; --j) {
            if(local_mean_v[j] <0) continue;
            last_good_index = j;
            last_local_mean = local_mean_v[j];
            last_local_rms  = local_sigma_v[j];
            break;
          }
        }
      }


    }
     */



    // **********
    // Now look at special cases, if wf starts or
    // ends with a pulse
    // **********

    // At start

    bool end_found = false;

    local_mean = mean(wf, 0, _sample_size);
    local_rms  = std(wf, local_mean, 0, _sample_size);

    if (local_rms >= _threshold) {

      for (size_t i = 1; i < wf.size() - _sample_size; i++) {

        local_mean = mean(wf, i, _sample_size);
        local_rms  = std(wf, local_mean, i, _sample_size);

        if (local_rms < _threshold) {

          end_found = true;

          for (size_t j = 0; j < i; j++){
            mean_temp_v[j] = local_mean;
            sigma_v[j] = local_rms;
            ped_interapolated[j] = true;
          }
          break;
        }
      }

      if (!end_found) {
        std::cerr <<"\033[93m<<" << __FUNCTION__ << ">>\033[00m Could not find good pedestal for CDF"
                  << "There is pulse on first sample and baseline never went back down. Returning false here.";
        return false;
      }

    }


    // At end

    bool start_found = false;

    local_mean = mean(wf, wf.size()-1-_sample_size, _sample_size);
    local_rms  = std(wf, local_mean, wf.size()-1-_sample_size, _sample_size);

    if (local_rms >= _threshold) {

      size_t i = wf.size() - 1 - _sample_size;
      while (i-- > 0) {
        local_mean = mean(wf, i, _sample_size);
        local_rms  = std(wf, local_mean, i, _sample_size);

        if (local_rms < _threshold) {

          start_found = true;

          for (size_t j = wf.size()-1; j > i; j--){
            mean_temp_v[j] = local_mean;
            sigma_v[j] = local_rms;
            ped_interapolated[j] = true;
          }
          break;
        }
      }

      if (!start_found) {
        std::cerr <<"\033[93m<<" << __FUNCTION__ << ">>\033[00m Could not find good pedestal for CDF"
                  << "There is pulse on last sample and baseline never went back down. Returning false here.";
        return false;
      }

    }







    // **********
    // Now smooth it to estimate the final pedestal
    // **********

    const size_t window_size = _sample_size*2;

    // middle mean
    for(size_t i=0; i< mean_temp_v.size(); ++i) {

      if( i < _sample_size || i >= (wf.size() - _sample_size) ) continue;

      mean_v[i]  = this->CalcMean (mean_temp_v,i - _sample_size,window_size);
      if(!ped_interapolated[i]){
        sigma_v[i] = this->CalcStd  (mean_temp_v,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];
      if(!ped_interapolated[i]){
        sigma_v[i] = sigma_v[_sample_size];
      }
    }

    // tail mean
    for(size_t i=(mean_temp_v.size() - _sample_size); i<mean_temp_v.size(); ++i) {

      mean_v[i]  = mean_v [wf.size() - _sample_size -1];
      if(!ped_interapolated[i]){
        sigma_v[i] = sigma_v[wf.size() - _sample_size -1];
      }
    }




    // Save to file
    if (_wf_saved + 1 <= _n_wf_to_csvfile) {
      _wf_saved ++;
      for (size_t i = 0; i < wf.size(); i++) {
        _csvfile << _wf_saved-1 << "," << i << "," << wf[i] << "," << mean_v[i] << "," << sigma_v[i] << std::endl;
      }
    }


    bool is_sane = this->CheckSanity(mean_v, sigma_v);

    return is_sane;

  }

void pmtana::PedAlgoRmsSlider::PrintInfo

(

)

Print settings.

Definition at line 51 of file PedAlgoRmsSlider.cxx.

  {
    std::cout << "PedAlgoRmsSlider setting:"
              << "\n\t SampleSize:       " << _sample_size
              << "\n\t Threshold:        " << _threshold
              << "\n\t Verbose:          " << _verbose
              << "\n\t NWaveformsToFile: " << _n_wf_to_csvfile << std::endl;
  }

Member Data Documentation

std::ofstream pmtana::PedAlgoRmsSlider::_csvfile

private

Definition at line 67 of file PedAlgoRmsSlider.h.

float pmtana::PedAlgoRmsSlider::_max_sigma

private

Max sigma to consider adc as 'sane'.

Definition at line 58 of file PedAlgoRmsSlider.h.

int pmtana::PedAlgoRmsSlider::_n_wf_to_csvfile

private

If greater than zero saves firsts waveforms with pedestal to csv file.

Definition at line 63 of file PedAlgoRmsSlider.h.

int pmtana::PedAlgoRmsSlider::_num_postsample

private

number of ADCs to sample after the gap

Definition at line 66 of file PedAlgoRmsSlider.h.

int pmtana::PedAlgoRmsSlider::_num_presample

private

number of ADCs to sample before the gap

Definition at line 65 of file PedAlgoRmsSlider.h.

float pmtana::PedAlgoRmsSlider::_ped_range_max

private

Max value of adc to consider adc as 'sane'.

Definition at line 59 of file PedAlgoRmsSlider.h.

float pmtana::PedAlgoRmsSlider::_ped_range_min

private

Min value of adc to consider adc as 'sane'.

Definition at line 60 of file PedAlgoRmsSlider.h.

size_t pmtana::PedAlgoRmsSlider::_sample_size

private

How many samples are used to calculate local rms and mean.

Definition at line 55 of file PedAlgoRmsSlider.h.

double pmtana::PedAlgoRmsSlider::_threshold

private

Threshold applied to local rms to claim a pulse.

Definition at line 56 of file PedAlgoRmsSlider.h.

bool pmtana::PedAlgoRmsSlider::_verbose

private

For debugging.

Definition at line 62 of file PedAlgoRmsSlider.h.

int pmtana::PedAlgoRmsSlider::_wf_saved = 0

private

Definition at line 64 of file PedAlgoRmsSlider.h.

---

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

• larana/larana/OpticalDetector/OpHitFinder/PedAlgoRmsSlider.h
• larana/larana/OpticalDetector/OpHitFinder/PedAlgoRmsSlider.cxx