PedAlgoRmsSlider.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//
//  PedAlgoRmsSlider source
//
////////////////////////////////////////////////////////////////////////

#include "PedAlgoRmsSlider.h"
#include "OpticalRecoException.h"
#include "UtilFunc.h"
#include "fhiclcpp/ParameterSet.h"

#include <iostream>
#include <fstream>
#include <numeric>

namespace pmtana{

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

  //**************************************************************************
  PedAlgoRmsSlider::PedAlgoRmsSlider(const fhicl::ParameterSet &pset,
                                             const std::string name)
    : 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;
    }

  }

  //*******************************************
  void PedAlgoRmsSlider::PrintInfo()
  //*******************************************
  {
    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;
  }

  //****************************************************************************
  double PedAlgoRmsSlider::CalcMean(const std::vector<double>& wf, size_t start, size_t nsample)
  //****************************************************************************
  {
    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 PedAlgoRmsSlider::CalcStd(const std::vector<double>& wf, const double ped_mean, size_t start, size_t nsample)
  //****************************************************************************
  {
    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 PedAlgoRmsSlider::ComputePedestal( const pmtana::Waveform_t& wf,
                                            pmtana::PedestalMean_t&   mean_v,
                                            pmtana::PedestalSigma_t&  sigma_v)
  //****************************************************************************
  {

    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;

  }





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

    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;

  }
}