LDFFunctions.h

Go to the documentation of this file.

#ifndef _LDFFinderOG_LDFFunctions_h_
#define _LDFFinderOG_LDFFunctions_h_


#include <cmath>
#include <utl/AugerUnits.h>
#include <utl/Math.h>

namespace LDFFinderOG {

  using namespace std;
  using namespace utl;


  enum LDFFunctionType {
    ePL,
    eNKG,
    eNKGFermi
  };


  namespace {

    inline
    double
    PowerLawFunction(const double r, const double beta, const double gamma = 0)
    {
      const double kNearRadius = 300*meter;
      const double k1000m = 1000*meter;
      const double nearCore = kNearRadius / k1000m;
      const double relR = r / k1000m;
      if (relR > nearCore)
        return pow(relR, beta + gamma*log(relR));

      const double gf = gamma*log(nearCore);
      return pow(relR, beta + 2*gf) * pow(nearCore, -gf);
    }

    inline
    double
    NKGFunction(const double r, const double beta, const double gamma = 0)
    {
      const double k700 = 700*meter;
      const double k1000 = 1000*meter;
      return pow(r/k1000, beta) *
        pow((k700 + r)/(k1000 + k700), beta+gamma);
    }

    inline
    double
    NKGFermiFunction(const double r, const double beta, const double gamma = 0,
                     const double mu = 2660.*meter, const double tau = 242.*meter)
    {
      const double k700 = 700*meter;
      const double k1000 = 1000*meter;
      return pow(r/k1000, beta) *
        pow((k700 + r)/(k1000 + k700), beta+gamma) / (exp((r - mu)/tau) + 1) *
        (exp((k1000 - mu)/tau) + 1); // normalization of LDF at 1000m
    }

  } // NS


  inline
  double
  LDFFunction(const LDFFunctionType type,
              const double r, const double beta, const double gamma = 0,
              const double mu = 2660.*meter, const double tau = 242.*meter)
  {
    if (type == ePL)
      return PowerLawFunction(r, beta, gamma);

    if (type == eNKG)
      return NKGFunction(r, beta, gamma);

    return NKGFermiFunction(r, beta, gamma, mu, tau);
  }


  namespace {

    inline
    double
    PowerLawSecondDerivative(const double r, const double beta, const double gamma = 0)
    {
      const double kNearRadius = 300*meter;
      const double k1000m = 1000*meter;
      const double nearCore = kNearRadius / k1000m;
      const double relR = r / k1000m;
      if (relR > nearCore) {
        const double tgl = 2 * gamma * log(relR);
        return PowerLawFunction(r, beta, gamma) * (
          (beta-1)*beta + 2*gamma + tgl*(tgl + 2*beta - 1)
        ) / Sqr(r);
      }

      const double bgl = beta + gamma * log(nearCore);
      return (bgl - 1)*bgl * pow(relR, bgl) / Sqr(r);
    }


    inline
    double
    NKGSecondDerivative(const double r, const double beta, const double gamma = 0)
    {
      const double k700 = 700*meter;
      const double tbg = 2*beta + gamma;
      const double r2 = Sqr(r);
      return NKGFunction(r, beta, gamma) * (
        Sqr(k700)*beta*(beta-1) +
        (tbg-1)*(2*k700*beta*r + tbg*r2)
      ) / (r2 * Sqr(r+k700));
    }

    inline
    double
    NKGFermiSecondDerivative(const double r, const double beta, const double gamma = 0,
                             const double mu = 2660.*meter, const double tau = 242.*meter)
    {
      const double k700 = 700*meter;
      //const double tbg = 2*beta + gamma;
      const double expofunc = exp((r - mu)/tau);
      const double fermi = 1 / (expofunc + 1);
      const double fermi2 = Sqr(fermi);
      const double nkgFirstDerivative = NKGFunction(r, beta, gamma) * (2.*beta*r + beta*k700 + gamma*r) / (r * (r + k700));
      return NKGSecondDerivative(r, beta, gamma) * fermi +
        nkgFirstDerivative * (expofunc/tau * fermi2) +
        NKGFunction(r, beta, gamma) * fermi2 / (tau*tau) * expofunc * (2*fermi + 1);
    }

  } // NS


  inline
  double
  LDFFunctionSecondDerivative(const LDFFunctionType type,
                              const double r, const double beta, const double gamma = 0,
                              const double mu = 2660.*meter, const double tau = 242.*meter)
  {
    if (type == ePL)
      return PowerLawSecondDerivative(r, beta, gamma);

    if (type == eNKG)
      return NKGSecondDerivative(r, beta, gamma);

    return NKGFermiSecondDerivative(r, beta, gamma, mu, tau);

  }


  inline
  double
  SigmaForFixedBeta(const double s1000)
  {
    const double lgS1000 = log10(s1000);
    return 0.428623 /* +-0.49688 */ * exp(-0.406203 /* +-0.612831*/ * lgS1000);
  }


  namespace {

    inline
    void
    PowerLawEstimateBetaGamma(double& beta, double& gamma,
                              const double cosTheta, const double s1000)
    {
      if (s1000 > 0) {
        const double log10s1000 = log10(s1000);
        const double secTheta = 1/cosTheta;
        const double sec2Theta = Sqr(secTheta);
        const double a0 = -4.73;
        const double a1 = -0.519;
        const double b0 = 1.32;
        const double b1 = 0.405;
        const double c0 = -0.105;
        const double c1 = -0.117;
        beta =
          a0 + a1*log10s1000 +
          (b0 + b1*log10s1000) * secTheta +
          (c0 + c1*log10s1000) * sec2Theta;
      } else
        beta = 0.7 * atan(6*(0.65 - cosTheta)) - 3;

      gamma = 0.05*sin(8*(cosTheta - 0.6)) - 0.5;
    }


    inline
    void
    NKGEstimateBetaGamma(double& beta, double& gamma,
                         const double cosTheta, const double s1000)
    {
      if (s1000 > 0) {
        const double log10s1000 = log10(s1000);
        const double secTheta = 1/cosTheta;
        const double sec2Theta = Sqr(secTheta);
        const double cos2Theta = Sqr(cosTheta);

        // parameters from GAP_2010_046
        const double a0 = -3.72;
        const double a1 = 0.0967;
        const double b0 = 1.74;
        const double b1 = -0.242;
        const double c0 = -0.274;
        const double c1 = 0.0349;
        /* parameters from GAP_2007_106
           const double a0 = -3.35;
           const double a1 = -0.125;
           const double b0 = 1.33;
           const double b1 = -0.0324;
           const double c0 = -0.191;
           const double c1 = -0.00573;
        */

        // parameters from GAP_2010_046
        const double fo0 = -1.87;
        const double fo1 = -0.183;
        const double fa0 =  0.490;
        const double fa1 = -0.0650;
        const double fp0 =  0.483;
        const double fp1 =  0.005;
        const double fs0 =  19.6;
        const double fs1 = -2.10;
        const double fb  = -0.272;
        const double fet =  2.32;
        const double fps =  1.95;
        const double fss =  18.01;

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

      } else {
        beta = 0.9 / cosTheta - 3.3;
        if (beta > -0.5)
          beta = -0.5;
        gamma = beta;
      }
    }

  } // NS

  inline
  void
  EstimateBetaGamma(double& beta, double& gamma,
                    const LDFFunctionType type,
                    const double cosTheta, const double s1000 = 0)
  {
    if (type == ePL)
      return PowerLawEstimateBetaGamma(beta, gamma, cosTheta, s1000);

    return NKGEstimateBetaGamma(beta, gamma, cosTheta, s1000);
  }

} // NS LDFFinderOG


#endif