QuadraticFitter.h

Go to the documentation of this file.

#ifndef _utl_QuadraticFitter_h_
#define _utl_QuadraticFitter_h_

#include <cmath>

#include <utl/Math.h>
#include <utl/Histogram.h>
#include <utl/AugerException.h>
#include <utl/QuadraticFitData.h>


namespace utl {

  class PoissonianErrors {
  public:
    static double GetError(const double& value)
    { return std::sqrt(value); }

    static double GetError2(const double& value)
    { return value; }
  };


  class FlatErrors {
  public:
    static double GetError(const double& /*value*/) { return 1; }
    static double GetError2(const double& /*value*/) { return 1; }
  };


  template<class Histogram, class ErrorPolicy /*= utl::PoissonianErrors*/>
  class QuadraticFitter {
  public:
    QuadraticFitter(const Histogram& h)
      : fStartBin(0), fStopBin(h.GetNBins()), fHistogram(h) { }

    QuadraticFitter(const Histogram& h,
                    const int startBin, const int stopBin)
      : fStartBin(startBin), fStopBin(stopBin), fHistogram(h)
    {
      if (startBin < 0 || stopBin - startBin < 3 || stopBin >= int(h.GetNBins()))
        throw utl::OutOfBoundException("start, stop bin out of bounds of histogram");
    }

    QuadraticFitter(const Histogram& h,
                    const double xStart, const double xStop)
      : fHistogram(h)
    {
      int start = h.GetBinIndex(xStart);
      if (start == Histogram::eUnderflow)
        start = 0;
      fStartBin = start;
      int stop = h.GetBinIndex(xStop);
      if (stop == Histogram::eOverflow)
        stop = h.GetNBins() - 1;
      fStopBin = stop;
      if (start < 0 || stop < 0 || stop - start < 3)
        throw utl::OutOfBoundException("illegal start or stop, or not enough bins to fit");
    }

    double
    GetExtremePosition()
    {
      MakeSums();
      NeedsCoeficient<1>();
      NeedsCoeficient<2>();
      return fMeanX - fK[1] / (2 * fK[2]);
    }

    double
    GetExtremePosition(double& posError, int& ndof)
    {
      const double pos = GetExtremePosition(); // delivers: 1, 2

      const double sigma2k1 = Sqr(fX2Sum) - fX0Sum * fX4Sum; // /det
      const double sigma2k2 = Sqr(fX1Sum) - fX0Sum * fX2Sum; // /det
      const double covk1k2  = fX0Sum * fX3Sum - fX1Sum * fX2Sum; // /det
      const double k22 = Sqr(fK[2]); // /det^2

      NeedsDeterminant();

      posError =
        0.5 * sqrt(fDet*(sigma2k1 + (pos - fMeanX)*covk1k2 + sigma2k2*Sqr(fK[1])/k22) / k22);

      ndof = fN - 3;

      return pos;
    }

    double
    GetExtremePosition(double& posError, double& chi2, int& ndof)
    {
      const double pos = GetExtremePosition(posError, ndof); // delivers: det, 1, 2

      NeedsCoeficient<0>();

      chi2 = ((Sqr(fK[0]) * fX0Sum + 2*fK[0]*fK[1] * fX1Sum + (Sqr(fK[1]) + 2*fK[0]*fK[2]) * fX2Sum
               + 2*fK[1]*fK[2] * fX3Sum + Sqr(fK[2]) * fX4Sum) / Sqr(fDet)
              - 2 * (fK[0] * fX0YSum + fK[1] * fX1YSum + fK[2] * fX2YSum) / fDet + fX0Y2Sum);

      return pos;
    }

    double
    GetExtremePosition(double& posError, double& chi2, int& ndof, double polynomialCoeff[3])
    {
      const double pos = GetExtremePosition(posError, chi2, ndof); // delivers: det, 0, 1, 2
      // this gives expansion around the fMeanX
      for (int i = 0; i < 3; ++i)
        polynomialCoeff[i] = fK[i] / fDet;
      // polynomial coefficient are exported as expansion around the position of the extreme
      polynomialCoeff[0] +=
        fMeanY -
        polynomialCoeff[1] * fMeanX +
        polynomialCoeff[2] * (Sqr(fMeanX) - Sqr(pos));
      polynomialCoeff[1] = 0;
      // polynomialCoeff[2] is the same
      return pos;
    }

    void
    GetFitData(QuadraticFitData& fit)
    {
      fit.fStart = fHistogram.GetBinCenter(fStartBin);
      fit.fStop = fHistogram.GetBinCenter(fStopBin);
      fit.fExtremePosition =
        GetExtremePosition(fit.fExtremePositionError, fit.fChi2, fit.fNdof, fit.fCoefficients);
    }

  private:
    void
    MakeSums()
    {
      fN = fStopBin - fStartBin;

      fMeanX = fMeanY = 0;
      for (unsigned int i = fStartBin; i < fStopBin; ++i) {
        fMeanX += fHistogram.GetBinCenter(i);
        fMeanY += fHistogram.GetBinAverage(i);
      }
      fMeanX /= fN;
      fMeanY /= fN;

      fX0Sum = fX1Sum = fX2Sum = fX3Sum = fX4Sum =
        fX0YSum = fX1YSum = fX2YSum = fX0Y2Sum = 0;

      for (unsigned int i = fStartBin; i < fStopBin; ++i) {
        const double x = fHistogram.GetBinCenter(i) - fMeanX;
        const double raw = fHistogram.GetBin(i);
        const double invSigma2 = Sqr(fHistogram.GetBinSize(i)) / ErrorPolicy::GetError2(raw);
        const double y = fHistogram.GetBinAverage(i) - fMeanY;
        fX0Sum += invSigma2;
        const double xs = x * invSigma2;
        fX1Sum += xs;
        const double ys = y * invSigma2;
        fX0YSum += ys;
        fX0Y2Sum += y * ys;
        fX1YSum += x * ys;
        const double x2s = x * xs;
        fX2Sum += x2s;
        fX2YSum += y * x2s;
        const double x3s = x * x2s;
        fX3Sum += x3s;
        fX4Sum += x * x3s;
      }
    }

    void
    NeedsDeterminant()
    {
      if (!fDet)
        fDet =   Sqr(fX2Sum) * fX2Sum + fX0Sum * Sqr(fX3Sum)
               + Sqr(fX1Sum) * fX4Sum
               - fX2Sum * (2 * fX1Sum * fX3Sum + fX0Sum * fX4Sum);
    }

    template<unsigned int index>
    void
    NeedsCoeficient()
    {
      switch (index) {
      case 0:
        if (!fK[0])
          fK[0] =   Sqr(fX3Sum) * fX0YSum + fX4Sum * (fX1Sum * fX1YSum - fX2Sum * fX0YSum)
                  + Sqr(fX2Sum) * fX2YSum - fX3Sum * (fX2Sum * fX1YSum + fX1Sum * fX2YSum); // /det
        break;
      case 1:
        if (!fK[1])
          fK[1] =   fX4Sum * (fX1Sum * fX0YSum - fX0Sum * fX1YSum) + Sqr(fX2Sum) * fX1YSum
                  + fX0Sum * fX3Sum * fX2YSum - fX2Sum * (fX3Sum * fX0YSum + fX1Sum * fX2YSum); // /det
        break;
      case 2:
        if (!fK[2])
          fK[2] =   Sqr(fX2Sum) * fX0YSum + fX3Sum * (fX0Sum * fX1YSum - fX1Sum * fX0YSum)
                  + Sqr(fX1Sum) * fX2YSum - fX2Sum * (fX1Sum * fX1YSum + fX0Sum * fX2YSum); // /det
        break;
      }
    }

    unsigned int fStartBin = 0;
    unsigned int fStopBin = 0;
    const Histogram& fHistogram;

    double fMeanX = 0;
    double fMeanY = 0;
    double fX0Sum = 0;
    double fX1Sum = 0;
    double fX2Sum = 0;
    double fX3Sum = 0;
    double fX4Sum = 0;
    double fX0YSum = 0;
    double fX1YSum = 0;
    double fX2YSum = 0;
    double fX0Y2Sum = 0;
    int fN = 0;
    double fDet = 0;
    double fK[3] = { 0 };
  };


  // convinience makers
  template<class Histogram, class ErrorPolicy>
  inline
  QuadraticFitter<Histogram, ErrorPolicy>
  MakeQuadraticFitter(const Histogram& h, const ErrorPolicy /*ep*/)
  {
    return QuadraticFitter<Histogram, ErrorPolicy>(h);
  }

  template<class Histogram>
  inline
  QuadraticFitter<Histogram>
  MakeQuadraticFitter(const Histogram& h)
  {
    return QuadraticFitter<Histogram>(h);
  }


  template<class Histogram, typename Type, class ErrorPolicy>
  inline
  QuadraticFitter<Histogram, ErrorPolicy>
  MakeQuadraticFitter(const Histogram& h,
                      const Type xStart, const Type xStop,
                      [[maybe_unused]] const ErrorPolicy ep)
  {
    return QuadraticFitter<Histogram, ErrorPolicy>(h, xStart, xStop);
  }

  template<class Histogram, typename Type>
  inline
  QuadraticFitter<Histogram>
  MakeQuadraticFitter(const Histogram& h,
                      const Type xStart, const Type xStop)
  {
    return QuadraticFitter<Histogram>(h, xStart, xStop);
  }

}


#endif