FdDoubleBumpFinderKG/ProfileFitFunctions.h

Go to the documentation of this file.

#ifndef _FdDoubleBumpFinderKG_ProfileFitFunctions_h_
#define _FdDoubleBumpFinderKG_ProfileFitFunctions_h_

#include <cmath>


namespace profileFit {

  using namespace std;

  enum EFitFunctionType {
    eGaisserHillas,
    eGaisserHillasUnconstrained,
    eDoubleGaisserHillas,
    eGaussInAge,
    eDoubleGaussInAge
  };


  enum EGHFunctionPar {
    eGHX0 = 0,
    eGHXmax,
    eGHnMax,
    eGHLambda,
    eNGHpars
  };


  inline
  double
  GHFunction(const double X,
             const double x0, const double lambda,
             const double xMax, const double nMax)
  {
    if (X <= x0)
      return 0;

    const double tmp = xMax - x0;

    return nMax * pow((X-x0)/tmp, tmp/lambda) *
      exp((xMax - X) / lambda);
  }


  enum EDoubleGHFunctionPar {
    eDGHX01 = 0,
    eDGHXmax1,
    eDGHnMax1,
    eDGHLambda1,
    eDGHX02,
    eDGHXmax2,
    eDGHnMax2,
    eDGHLambda2,
    eNDGHpars
  };


  inline
  double
  DoubleGHFunction(const double X,
                   const double x01, const double lambda1,
                   const double xMax1, const double nMax1,
                   const double x02, const double lambda2,
                   const double xMax2, const double nMax2)
  {
    return GHFunction(X, x01, lambda1, xMax1, nMax1) +
           GHFunction(X, x02, lambda2, xMax2, nMax2);

  }


  enum EGaussInAgeFunctionPar {
    eGiAXmax = 0,
    eGiAnMax,
    eGiASigma1,
    eGiASigma2,
    eNGiApars
  };


  inline
  double
  GaussInAgeFunction(const double X,
                     const double nMax,
                     const double xMax,
                     const double sigma1,
                     const double sigma2)
  {
    const double age = 3*X / (X + 2*xMax);
    const double sigma = (age < 1) ? sigma1 : sigma2;
    const double z = (age - 1) / sigma;
    return nMax * exp(-0.5*z*z);
  }

}


#endif