NKGAGLDF.h

Go to the documentation of this file.

#ifndef _LDFFinderKG_NKGAGLDF_h_
#define _LDFFinderKG_NKGAGLDF_h_


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


namespace LDFFinderKG {

  class NKGAGLDF : public VLDF {

  public:
    NKGAGLDF(const double refDistance) : VLDF(refDistance, 10, 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 k700 = 700*utl::meter;
      const double k3000 = 3000*utl::meter;

      return std::pow(r/rRef, beta) *
        std::pow((k700 + r)/(rRef + k700), beta)  *
        std::pow((1 + std::pow(r/k3000, 3)) / (1 + std::pow(rRef/k3000, 3)), gamma);
    }

    double
    SecondDerivative(const double r, const std::vector<double>& shape)
      const
    {
      const double beta = shape[0];
      const double gamma = shape[1];

      const double k700 = 700*utl::meter;
      const double k3000 = 3000*utl::meter;

      const double r2 = utl::Sqr(r);
      const double r700_2 = utl::Sqr(r + k700);
      const double r3000_2 = utl::Sqr(r + k3000);

      return Value(r, shape) *
        (std::pow(beta/r + beta/(r + k700) + 3*gamma/(r + k3000), 2) - (beta/r2 + beta/r700_2 + 3*gamma/r3000_2));
    }

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

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

      const double a0_Beta = fShapeModelVector[0];
      const double a1_Beta = fShapeModelVector[1];
      const double b0_Beta = fShapeModelVector[2];
      const double c0_Beta = fShapeModelVector[3];
      const double d0_Beta = fShapeModelVector[4];

      const double a0 = fShapeModelVector[5];
      const double a1 = fShapeModelVector[6];
      const double b0 = fShapeModelVector[7];
      const double c0 = fShapeModelVector[8];
      const double d0 = fShapeModelVector[9];

      beta = a0_Beta + a1_Beta*lgSRef + secTheta*(b0_Beta + secTheta*(c0_Beta + secTheta*d0_Beta));

      gamma = a0 + a1*lgSRef + secTheta*(b0 + secTheta*(c0 + secTheta*d0));

      if (gamma > 0)
        gamma = 0;

      if (secTheta > 3)
        beta = a0_Beta + a1_Beta*lgSRef + 3*b0_Beta;

      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 2; }

  };

}


#endif