NKGScintillatorLDF.h

Go to the documentation of this file.

#ifndef _ScintillatorLDFFinderKG_NKGScintillatorLDF_h_
#define _ScintillatorLDFFinderKG_NKGScintillatorLDF_h_


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


namespace ScintillatorLDFFinderKG {

  class NKGScintillatorLDF : public VLDF {

  public:
    NKGScintillatorLDF(const double refDistance) : VLDF(refDistance, 12, 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;

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

    double
    FirstDerivative(const double r, const std::vector<double>& shape)
      const
    {
      const double beta = shape[0];
      const double gamma = shape[1];
      const double rRef  = fReferenceDistance;
      const double k700 = 700*utl::meter;
      return  (beta*std::pow(r/rRef,beta-1)
        * std::pow((r+k700)/(k700+rRef),beta+gamma)/rRef) 
        + (std::pow(r/rRef,beta)*std::pow((r+k700)/(k700+rRef),beta+gamma-1) 
        * (beta + gamma))/(k700+rRef);
    }

    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 tbg = 2*beta + gamma;
      const double r2 = utl::Sqr(r);
      return Value(r, shape) * (
        utl::Sqr(k700)*beta*(beta - 1) +
        (tbg - 1)*(2*k700*beta*r + tbg*r2)
      ) / (r2 * utl::Sqr(r + k700));
    }

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

      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 d0 = fShapeModelVector[6];
      const double d1 = fShapeModelVector[7];
      const double e0 = fShapeModelVector[8];
      const double e1 = fShapeModelVector[9];
      const double f0 = fShapeModelVector[10];
      const double f1 = fShapeModelVector[11];

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

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

      gamma =
        d0 + d1*lgSRef +
        (e0 + e1*lgSRef) * secTheta +
        (f0 + f1*lgSRef) * sec2Theta - 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