ExponentialFitter.h

Go to the documentation of this file.

#ifndef _utl_ExponentialFitter_h_
#define _utl_ExponentialFitter_h_

#include <utl/ExponentialFitData.h>
#include <utl/AugerException.h>


namespace utl {

  template<class Histogram>
  class ExponentialFitter {
  public:
    ExponentialFitter(const Histogram& histogram,
                      const int startBin, const int stopBin)
      : fStartBin(startBin), fStopBin(stopBin), fHistogram(histogram)
    {
      if (fStartBin < 0 || fStopBin >= histogram.GetNBins())
        throw OutOfBoundException("start, stop bins out of bound");
      if (fStopBin - fStartBin < 2)
        throw OutOfBoundException("too few data points to fit exponential");
    }

    ExponentialFitter(const Histogram& histogram,
                      const double xStart, const double xStop)
      : fStartBin(histogram.GetBinIndex(xStart)),
        fStopBin(histogram.GetBinIndex(xStop)),
        fHistogram(histogram)
    {
      if (fStartBin < 0)
        throw OutOfBoundException("range start error");
      if (fStopBin - fStartBin < 2)
        throw OutOfBoundException("range order mismatch");
    }

    double
    GetSlope(double& slopeError, const double offset = 0)
    {
      // assumes sigma(y) = sqrt(y) without offset
      // z = log(y)
      fSum1 = 0;
      fSumx = 0;
      fSumx2 = 0;
      fSumxz = 0;
      fSumz = 0;
      fN = 0;
      for (int i = fStartBin; i < fStopBin; ++i) {
        const double counts = fHistogram.GetBin(i);
        if (counts > 0) {
          const double width = fHistogram.GetBinSize(i);
          const double y = (counts - offset) / width;
          if (y > 0) {
            const double invSigmaz2 = counts;
            const double z = log(y);
            const double x = fHistogram.GetBinCenter(i);
            fSum1 += invSigmaz2;
            const double s2x = invSigmaz2 * x;
            fSumx += s2x;
            fSumx2 += s2x * x;
            fSumxz += s2x * z;
            fSumz += invSigmaz2 * z;
            ++fN;
          }
        }
      }

      if (fN > 1) {
        fDet = fSum1*fSumx2 - Sqr(fSumx);
        if (fDet) {
          slopeError = fSum1 / fDet;
          return (fSum1*fSumxz - fSumx*fSumz) / fDet;
        }
      }

      fDet = slopeError = 0;
      return 0;
    }

    void
    GetFit(double& amplitude, double& amplitudeError,
           double& slope, double& slopeError,
           const double offset = 0)
    {
      slope = GetSlope(slopeError, offset);
      if (fDet) {
        amplitude = exp((fSumx2*fSumz - fSumx*fSumxz) / fDet);
        amplitudeError = amplitude * (fSumx2 / fDet);
      } else {
        amplitude = 0;
        amplitudeError = 0;
      }
    }

    void
    GetFit(double& amplitude, double& amplitudeError,
           double& slope, double& slopeError,
           double& chi2, int& ndof,
           const double offset = 0)
    {
      GetFit(amplitude, amplitudeError, slope, slopeError, offset);
      if (fDet) {
        const double linOffset = log(amplitude);
        chi2 = Sqr(linOffset)*fSum1 + Sqr(slope)*fSumx2 + 2*(slope*linOffset*fSumx - linOffset*fSumz - slope*fSumxz);
        ndof = fN - 2;
      } else {
        chi2 = 0;
        ndof = 0;
      }
    }

    void
    GetFit(ExponentialFitData& ef, const double offset = 0)
    {
      ef.fStart = fHistogram.GetBinCenter(fStartBin);
      ef.fStop = fHistogram.GetBinCenter(fStopBin);
      GetFit(ef.fAmplitude, ef.fAmplitudeError, ef.fSlope, ef.fSlopeError, ef.fChi2, ef.fNdof, offset);
      ef.fOffset = offset;
    }

  private:
    const int fStartBin;
    const int fStopBin;
    const Histogram& fHistogram;

    double fSum1;
    double fSumx;
    double fSumx2;
    double fSumxz;
    double fSumz;
    double fDet;
    int fN;
  };


  template<class Histogram, typename T>
  inline
  ExponentialFitter<Histogram>
  MakeExponentialFitter(const Histogram& histogram,
                        const T& start, const T& stop)
  {
    return ExponentialFitter<Histogram>(histogram, start, stop);
  }


  template<class Histogram, typename T>
  inline
  ExponentialFitter<Histogram>
  MakeExponentialFitter(const Histogram& histogram,
                        const std::pair<T, T>& startStop)
  {
    return ExponentialFitter<Histogram>(histogram, startStop.first, startStop.second);
  }

}

#endif