PowerLawLDF.h

Go to the documentation of this file.

#ifndef _LDFFinderKG_PowerLawLDF_h_
#define _LDFFinderKG_PowerLawLDF_h_


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


namespace LDFFinderKG {

  class PowerLawLDF: public VLDF {

  public:
    PowerLawLDF(const double refDistance) : VLDF(refDistance, 6, 2) { }

    double
    Value(const double r, const std::vector<double>& shape)
      const
    {
      const double kNearRadius = 300*utl::meter;
      const double rRef = fReferenceDistance;
      const double beta = shape[0];
      const double gamma = shape[1];

      const double nearCore = kNearRadius / rRef;
      const double relR = r / rRef;
      if (relR > nearCore)
        return std::pow(relR, beta + gamma*std::log(relR));

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

    double
    SecondDerivative(const double r, const std::vector<double>& shape)
      const
    {
      const double kNearRadius = 300*utl::meter;
      const double rRef = fReferenceDistance;
      const double beta = shape[0];
      const double gamma = shape[1];

      const double nearCore = kNearRadius / rRef;
      const double relR = r / rRef;
      if (relR > nearCore) {
        const double tgl = 2 * gamma * std::log(relR);
        return Value(r, shape) * (
          (beta-1)*beta + 2*gamma + tgl*(tgl + 2*beta - 1)
        ) / utl::Sqr(r);
      }

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

    std::vector<double>
    ShapeModel(const double cosTheta, const double showerSize)
      const
    {
      const double lgSize = std::log10(showerSize);
      const double secTheta = 1 / 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];

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

      beta = a0 + a1*lgSize + secTheta*(b0 + b1*lgSize + secTheta*(c0 + c1*lgSize));

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

      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