NKGFermiLDF.h

Go to the documentation of this file.

#ifndef _LDFFinderKG_NKGFermiLDF_h_
#define _LDFFinderKG_NKGFermiLDF_h_


#include "VLDF.h"
#include <cmath>
#include <vector>
#include <utl/Math.h>
#include <utl/AugerUnits.h>


namespace LDFFinderKG {

  class NKGFermiLDF : public VLDF {

  public:
    NKGFermiLDF(const double refDistance) : VLDF(refDistance, 20, 2) { }

    double
    Value(const double r, const std::vector<double>& shape)
      const
    {
      const double rRef  = fReferenceDistance;
      const double beta  = shape[0];
      const double gamma = shape[1];
      const double mu = shape[2];
      const double tau = shape[3];

      return NKG(r, rRef, beta, gamma) * FermiSuppression(r, rRef, mu, tau);
    }

    double
    SecondDerivative(const double r, const std::vector<double>& shape)
      const
    {
      const double rRef = fReferenceDistance;
      const double beta = shape[0];
      const double gamma = shape[1];
      const double mu = shape[2];
      const double tau = shape[3];
      const double k700 = 700*utl::meter;

      const double expofunc = std::exp((r - mu)/tau);
      const double fermi = 1 / (expofunc + 1);
      const double fermi2 = utl::Sqr(fermi);
      const double nkgFirstDerivative =
        NKG(r, rRef, beta, gamma) * (2*beta*r + beta*k700 + gamma*r) / (r * (r + k700));
      return NKGSecondDerivative(r, rRef, beta, gamma) * fermi +
        nkgFirstDerivative * (expofunc/tau * fermi2) +
        NKG(r, rRef, beta, gamma) * fermi2 / (tau*tau) * expofunc * (2*fermi + 1);
    }

    std::vector<double>
    ShapeModel(const double cosTheta, const double showerSize)
      const
    {
      const double lgSRef = std::log10(showerSize);
      const double secTheta = 1 / cosTheta;
      const double cos2Theta = utl::Sqr(cosTheta);

      const double a0 = fShapeModelVector[0];
      const double a1 = fShapeModelVector[1];
      const double b0 = fShapeModelVector[2];
      const double b1 = fShapeModelVector[3];
      const double c0 = fShapeModelVector[4];
      const double c1 = fShapeModelVector[5];

      const double fo0 = fShapeModelVector[6];
      const double fo1 = fShapeModelVector[7];
      const double fa0 = fShapeModelVector[8];
      const double fa1 = fShapeModelVector[9];
      const double fp0 = fShapeModelVector[10];
      const double fp1 = fShapeModelVector[11];
      const double fs0 = fShapeModelVector[12];
      const double fs1 = fShapeModelVector[13];
      const double fb  = fShapeModelVector[14];
      const double fet = fShapeModelVector[15];
      const double fps = fShapeModelVector[16];
      const double fss = fShapeModelVector[17];

      std::vector<double> shape(4);
      double& beta = shape[0];
      double& gamma = shape[1];
      double& mu = shape[2];
      double& tau = shape[3];

      beta = a0 + a1*lgSRef + secTheta*(b0 + b1*lgSRef + secTheta*(c0 + c1*lgSRef));
      gamma =
        fo0 + fo1 * lgSRef +
        (fa0 + fa1 * lgSRef) / (exp((fs0 + fs1 * lgSRef) * (cos2Theta - (fp0 + fp1 * lgSRef))) + 1) +
        fb * pow(cos2Theta, fet) / (exp((lgSRef - fps) * fss) + 1) -
        beta; // define as modification of exponent on second term for backwards compatibility

      mu = fShapeModelVector[18];
      tau = fShapeModelVector[19];

      return shape;
    }

    double
    BetaUncertainty(const double showerSize)
      const
    {
      const double lgS = std::log10(showerSize);
      return fBetaUncertaintyModelVector[0] * exp(fBetaUncertaintyModelVector[1] * lgS);
    }

    unsigned int GetNShapeParameters() const { return 4; }

  protected:
    double
    FermiSuppression(const double r, const double rRef, const double mu, const double tau)
      const
    {
      // Fermi suppression term which conserves normalization of LDF at 1000m
      return (exp((rRef - mu)/tau) + 1) / (exp((r - mu)/tau) + 1);
    }

    double
    NKG(const double r, const double rRef, const double beta, const double gamma)
      const
    {
      const double k700 = 700*utl::meter;
      return pow(r/rRef, beta) *
        pow((k700 + r) / (rRef + k700), beta+gamma);
    }

    double
    NKGSecondDerivative(const double r, const double rRef, const double beta, const double gamma)
      const
    {
      const double k700 = 700*utl::meter;
      const double tbg = 2*beta + gamma;
      const double r2 = utl::Sqr(r);
      return NKG(r, rRef, beta, gamma) * (
        utl::Sqr(k700)*beta*(beta - 1) +
        (tbg - 1)*(2*k700*beta*r + tbg*r2)
      ) / (r2 * utl::Sqr(r+k700));
    }

  };

}


#endif