AlgoThreshold.cxx

Go to the documentation of this file.

//
//  AlgoThreshold source
//

#include "fhiclcpp/ParameterSet.h"

#include "AlgoThreshold.h"

namespace pmtana {

  //***************************************************************************
  AlgoThreshold::AlgoThreshold(const std::string name) : PMTPulseRecoBase(name)
  //***************************************************************************
  {
    //_adc_thres = 3;
    //_nsigma = 5;
    Reset();
  }

  //************************************************************
  AlgoThreshold::AlgoThreshold(const fhicl::ParameterSet& pset,
                               std::unique_ptr<pmtana::RiseTimeCalculatorBase> risetimecalculator,
                               //AlgoThreshold::AlgoThreshold(const ::fcllite::PSet &pset,
                               const std::string name)
    : PMTPulseRecoBase(name)
  //*******************************************************
  {

    _start_adc_thres = pset.get<double>("StartADCThreshold");
    _end_adc_thres = pset.get<double>("EndADCThreshold");

    //_nsigma = pset.get<double>("NSigmaThreshold");

    _nsigma_start = pset.get<double>("NSigmaThresholdStart");
    _nsigma_end = pset.get<double>("NSigmaThresholdEnd");

    _risetime_calc_ptr = std::move(risetimecalculator);

    Reset();
  }

  //***************************************************************
  void AlgoThreshold::Reset()
  //***************************************************************
  {
    PMTPulseRecoBase::Reset();
  }

  //***************************************************************
  bool AlgoThreshold::RecoPulse(const Waveform_t& wf,
                                const PedestalMean_t& mean_v,
                                const PedestalSigma_t& sigma_v)
  //***************************************************************
  {
    bool fire = false;

    double counter = 0;

    double ped_mean = mean_v.front();
    double ped_rms = sigma_v.front();

    //double threshold = ( _adc_thres > (_nsigma * ped_rms) ? _adc_thres : (_nsigma * ped_rms) );
    auto start_threshold =
      (_start_adc_thres > (_nsigma_start * ped_rms) ? _start_adc_thres : (_nsigma_start * ped_rms));
    auto end_threshold =
      (_end_adc_thres > (_nsigma_end * ped_rms) ? _end_adc_thres : (_nsigma_end * ped_rms));

    //    threshold += ped_mean

    start_threshold += ped_mean;
    end_threshold += ped_mean;

    Reset();

    for (auto const& value : wf) {

      if (!fire && ((double)value) >= start_threshold) {

        // Found a new pulse

        fire = true;

        _pulse.ped_mean = ped_mean;
        _pulse.ped_sigma = ped_rms;

        //vic: i move t_start back one, this helps with porch

        _pulse.t_start = counter - 1 > 0 ? counter - 1 : counter;
        //std::cout << "counter: " << counter << " tstart : " << _pulse.t_start << "\n";
      }

      if (fire && ((double)value) < end_threshold) {

        // Found the end of a pulse

        fire = false;

        //vic: i move t_start forward one, this helps with tail
        _pulse.t_end = counter < wf.size() ? counter : counter - 1;

        if (_risetime_calc_ptr)
          _pulse.t_rise = _risetime_calc_ptr->RiseTime(
            {wf.begin() + _pulse.t_start, wf.begin() + _pulse.t_end},
            {mean_v.begin() + _pulse.t_start, mean_v.begin() + _pulse.t_end},
            true);

        _pulse_v.push_back(_pulse);

        _pulse.reset_param();
      }

      //std::cout << "\tFire=" << fire << std::endl;

      if (fire) {

        // Add this adc count to the integral

        _pulse.area += ((double)value - (double)ped_mean);

        if (_pulse.peak < ((double)value - (double)ped_mean)) {

          // Found a new maximum

          _pulse.peak = ((double)value - (double)ped_mean);

          _pulse.t_max = counter;
        }
      }

      counter++;
    }

    if (fire) {

      // Take care of a pulse that did not finish within the readout window.

      fire = false;

      _pulse.t_end = counter - 1;

      if (_risetime_calc_ptr)
        _pulse.t_rise = _risetime_calc_ptr->RiseTime(
          {wf.begin() + _pulse.t_start, wf.begin() + _pulse.t_end},
          {mean_v.begin() + _pulse.t_start, mean_v.begin() + _pulse.t_end},
          true);

      _pulse_v.push_back(_pulse);

      _pulse.reset_param();
    }

    return true;
  }

}