UtilFunc.cxx

Go to the documentation of this file.

#include "UtilFunc.h"
#include "OpticalRecoException.h"
#include "larana/OpticalDetector/OpHitFinder/OpticalRecoTypes.h"

#include <algorithm>
#include <numeric>
#include <cmath>

#include "TH1D.h"

namespace pmtana {

  double mean(const std::vector<short>& 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) / ((double)nsample);

    return sum;
  }


  double edge_aware_mean(const std::vector<short>& wf, int start, int end) {

    auto m   = double{0.0};
    auto n_t = unsigned{0};

    for(int k = start; k < end; ++k) {
      if (k < 0 or k > (int)(wf.size()) - 1) continue;
      m += wf.at(k);
      ++n_t;
    }

    if( n_t > 0 ) m /= n_t;
    n_t = 0;

    return m;
  }

  double std(const std::vector<short>& 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;
  }

  double BinnedMaxOccurrence(const PedestalMean_t& mean_v,const size_t nbins)
  {
    if(nbins<1) throw OpticalRecoException("Cannot have 0 binning");

    auto res = std::minmax_element(std::begin(mean_v),std::end(mean_v));

    double bin_width = ((*res.second) - (*res.first)) / ((double)nbins);

    if(nbins==1 || bin_width == 0) return ((*res.first) + bin_width /2.);

    //std::cout<<"Min: "<<(*res.first)<<" Max: "<<(*res.second)<<" Width: "<<bin_width<<std::endl;

    // Construct array of nbins
    static std::vector<size_t> ctr_v(nbins,0);
    for(auto& v : ctr_v) v=0;
    for(auto const& v : mean_v) {

      size_t index = int((v - (*res.first))/bin_width);
      //std::cout<<"adc = "<<v<<" width = "<<bin_width<< " ... "
      //<<index<<" / "<<ctr_v.size()<<std::endl;

      ctr_v[index]++;

    }

    // Find max occurrence
    auto max_it = std::max_element(std::begin(ctr_v),std::end(ctr_v));

    // Get the mean of max-occurrence bins
    double mean_max_occurrence = 0;
    double num_occurrence = 0;
    for(size_t bin=0; bin<ctr_v.size(); ++bin) {

      if(ctr_v[bin] != (*max_it)) continue;

      mean_max_occurrence += ((*res.first) + bin_width / 2. + bin_width * bin);

      num_occurrence += 1.0;
    }

    return (mean_max_occurrence / num_occurrence);
  }


  // template<typename W>
  int sign(double val) {

    if (val > 0) return  1;
    if (val < 0) return -1;
    return 0;

  }

  double BinnedMaxTH1D(const std::vector<double>& v ,int bins){

    auto max_it = std::max_element(std::begin(v), std::end(v));
    auto min_it = std::min_element(std::begin(v), std::end(v));

    TH1D th("th",";;",bins,*min_it,*max_it);

    for (const auto & m : v) th.Fill(m);

    return th.GetXaxis()->GetBinCenter(th.GetMaximumBin());
  }


}