GaisserHillas4Parameter.cc

Go to the documentation of this file.

#include <evt/GaisserHillas4Parameter.h>
#include <utl/Math.h>
#include <utl/MathConstants.h>
#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/LambertW.h>
#include <utl/Deprecator.h>

#include <cmath>
#include <limits>

using namespace utl;
using namespace std;
using namespace evt::gh;

namespace evt {

  utl::TabulatedFunction GaisserHillas4Parameter::fRvsAsymTable;

  GaisserHillas4Parameter::GaisserHillas4Parameter(const EFunctionType functionType) :
    fFunctionType(functionType),
    fRhoShapePar1Par2(0)
  {
    for (int i = 0; i <= eLast; ++i) {
      fShapePar[i] = 0;
      fShapeParError[i] = 0;
      fRhoNMaxShapePar[i] = 0;
      fRhoXMaxShapePar[i] = 0;
    }
  }

  GaisserHillas4Parameter::GaisserHillas4Parameter(const double nMax, const double xMax,
                                                   const ShapeParameter par1,
                                                   const ShapeParameter par2,
                                                   const EFunctionType functionType) :
    VGaisserHillasParameter(nMax, xMax),
    fFunctionType(functionType),
    fRhoShapePar1Par2(0)
  {
    for (int i = 0; i <= eLast; ++i) {
      fShapeParError[i] = 0;
      fRhoNMaxShapePar[i] = 0;
      fRhoXMaxShapePar[i] = 0;
    }
    SetShapeParameter(par1.first, par1.second, 0);
    SetShapeParameter(par2.first, par2.second, 0);
  }

  double
  GaisserHillas4Parameter::Eval(const double depth)
    const
  {
    
    if (fFunctionType == eUSP) {
      // According to:  Eq. (4), APP 34 (2011) 360
      const double L = GetShapeParameter(eUSP_L);
      const double R = GetShapeParameter(eUSP_R);
      const double Y = (depth - fXMax)/L;

      double GH = 0;
      if(std::abs(R) < 5.e-7)
        GH = exp(-0.5*pow(Y, 2));

      if(R*Y > -1.){
        const double p1 = log(1.+R*Y)-R*Y;
        const double p2 = p1*pow(R,-2.);
        GH = fNMax * exp(p2);
      }

      double x0 = -(L/R) + fXMax;
      if (depth <= x0)
        return 0;
      else
        return GH;
    } 
    else {
      const double x0 = GetShapeParameter(eX0);
      const double lam = GetShapeParameter(eLambda);

      if (depth <= x0)
        return 0;
      else {
        const double tmp = fXMax - x0;
        return fNMax * pow((depth - x0) / tmp, tmp/lam) *
          exp((fXMax - depth) / lam);
      }
    }
  }

  void
  GaisserHillas4Parameter::Dump(ostream& os)
    const
  {
    os << "GH function with 4 parameters" << '\n'
       << "  Nmax  = " << fNMax << " +/- "<< fNMaxError << "\n"
       << "  Xmax  = " << fXMax/g*cm2 << " +/- "<< fXMaxError/g*cm2
       << " g/cm^2 \n";
    for (int i = 0; i < eLast; ++i)
      os << gh::GetShapeParameterName(EShapeParameter(i)) << " = "
         << fShapePar[i]/g*cm2 << " +/- "<< fShapeParError[i]/g*cm2;
    os << endl;
  }

  double
  GaisserHillas4Parameter::GetIntegral()
    const
  {
    const double w = (fXMax - GetShapeParameter(eX0)) / GetShapeParameter(eLambda);

    if (w <= 0)
      return 0;

    const double exponent = w + LogGamma(w+1);

    if (exponent < std::numeric_limits<double>::max_exponent10 * kLn10)
      return fNMax * GetShapeParameter(eLambda) * pow(w, -w) * exp(exponent);
    else
      return std::numeric_limits<double>::max();
  }


  double
  GaisserHillas4Parameter::GetIntegralError()
    const
  {
    const double w0 = (fXMax - GetShapeParameter(eX0)) / GetShapeParameter(eLambda);

    const double wError =
      sqrt(Sqr(fXMaxError) + Sqr(GetShapeParameterError(eX0))) /
      GetShapeParameter(eLambda);

    const double w1 = w0 + wError;
    const double w2 = w0 - wError;

    if (w0 <= 0 || w2 <= 0)
      return std::numeric_limits<double>::max();

    const double exponent0 = w0 * (1 - log(w0)) + LogGamma(w0 + 1);
    const double exponent1 = w1 * (1 - log(w1)) + LogGamma(w1 + 1);
    const double exponent2 = w2 * (1 - log(w2)) + LogGamma(w2 + 1);

    const double maxExp = std::numeric_limits<double>::max_exponent10 * kLn10;

    if (w0 > maxExp || w1 > maxExp || w2 > maxExp)
      return std::numeric_limits<double>::max();

    const double gamma0 = exp(exponent0);
    const double gamma1 = exp(exponent1);
    const double gamma2 = exp(exponent2);

    const double f1 = 0.5 * (gamma1 - gamma2) / gamma0;
    const double f2 = fNMaxError / fNMax;

    return sqrt(Sqr(f1) + Sqr(f2));
  }


  double
  GaisserHillas4Parameter::GetShapeParameter(const EShapeParameter par)
    const
  {

    switch (fFunctionType) {
    case eClassic:
      {
        switch (par)
          {
          case eX0:
          case eLambda:
            return fShapePar[ExternalToInternal(par)];
          case eFWHM:
            return InverseRight(0.5) - InverseLeft(0.5);
          case eAsym:
            {
              const double left = fXMax - InverseLeft(0.5);
              const double right = InverseRight(0.5) - fXMax;
              return left / (left + right);
            }
          case eUSP_L:
            {
              const double X0 = fShapePar[ExternalToInternal(eX0)];
              const double lambda = fShapePar[ExternalToInternal(eLambda)];
              const double X0prime = X0 - fXMax;
              return sqrt(abs(X0prime*lambda));
            }
          case eUSP_R:
            {
              const double X0 = fShapePar[ExternalToInternal(eX0)];
              const double lambda = fShapePar[ExternalToInternal(eLambda)];
              const double X0prime = X0 - fXMax;
              return sqrt(lambda/abs(X0prime));
            }
          }
        break;
      }
    case eWidth:
      {
        switch (par)
          {
          case eX0:
            {
              const double asym = fShapePar[ExternalToInternal(eAsym)];
              const double fwhm = fShapePar[ExternalToInternal(eFWHM)];
              const double R = CalculateR(asym);
              return fXMax - fwhm / R;
            }
          case eLambda:
            {
              const double asym = fShapePar[ExternalToInternal(eAsym)];
              const double fwhm = fShapePar[ExternalToInternal(eFWHM)];
              const double R = CalculateR(asym);
              return fwhm * (log(1 - asym * R) + asym * R) / R / log(0.5);
            }
          case eFWHM:
          case eAsym:
            return fShapePar[ExternalToInternal(par)];
          case eUSP_L:
            {
              const double asym = fShapePar[ExternalToInternal(eAsym)];
              const double fwhm = fShapePar[ExternalToInternal(eFWHM)];
              const double R = CalculateR(asym);
              const double X0 = fXMax - fwhm / R;
              const double lambda = fwhm * (log(1 - asym * R) + asym * R) / R / log(0.5);
              const double X0prime = X0 - fXMax;
              return sqrt(abs(X0prime*lambda));
            }
          case eUSP_R:
            {
              const double asym = fShapePar[ExternalToInternal(eAsym)];
              const double fwhm = fShapePar[ExternalToInternal(eFWHM)];
              const double R = CalculateR(asym);
              const double X0 = fXMax - fwhm / R;
              const double lambda = fwhm * (log(1 - asym * R) + asym * R) / R / log(0.5);
              const double X0prime = X0 - fXMax;
              return sqrt(lambda/abs(X0prime));
            }
          }
        break;
      }
    case eUSP:
      {
        switch (par)
          {
          case eX0:
            {
              const double L = fShapePar[ExternalToInternal(eUSP_L)];
              const double R = fShapePar[ExternalToInternal(eUSP_R)];
              //L = sqrt(abs(X0prime*lambda));
              //R = sqrt(lambda/abs(X0prime));
              //L2 = X0prime * lambda;
              //R2 = lambda / abs(X0prime);
              //L2 = X0prime * R2 * abs(X0prime) = R2 * X0prime2
              const double X0prime = -sqrt(L*L/(R*R));
              //const double lambda = R*R * abs(X0prime);
              return X0prime + fXMax;
            }
          case eLambda:
            {
              const double L = fShapePar[ExternalToInternal(eUSP_L)];
              const double R = fShapePar[ExternalToInternal(eUSP_R)];
              const double X0prime = -sqrt(L*L/(R*R));
              const double lambda = R*R * abs(X0prime);
              return lambda;
            }
          case eFWHM:
            return InverseRight(0.5) - InverseLeft(0.5);
          case eAsym:
            {
              const double left = fXMax - InverseLeft(0.5);
              const double right = InverseRight(0.5) - fXMax;
              return left / (left + right);
            }
          case eUSP_L:
          case eUSP_R:
            return fShapePar[ExternalToInternal(par)];
          }
        break;
      }
    }

    ostringstream errMsg;
    errMsg << "no shape parameter available for "
           << GetShapeParameterName(par)
           << ", GH type " << GetFunctionTypeName(fFunctionType);
    throw DoesNotComputeException(errMsg.str());

  }


  double
  GaisserHillas4Parameter::GetShapeParameterError(const EShapeParameter par)
    const
  {

    switch (fFunctionType) {
    case eClassic:
      {
        switch (par)
          {
          case eX0:
          case eLambda:
            return fShapeParError[ExternalToInternal(par)];
          case eFWHM: 
          case eAsym:
          case eUSP_L:
          case eUSP_R:
            return 0; // error propagation not yet implemented
          }
        break;
      }
    case eWidth:
      {
        switch (par)
          {
          case eFWHM:
          case eAsym:
            return fShapeParError[ExternalToInternal(par)];
          case eX0:
          case eLambda:
          case eUSP_L:
          case eUSP_R:
            return 0; // error propagation not yet implemented
          }
        break;
      }
    case eUSP:
      {
        switch (par)
          {
          case eUSP_L:
          case eUSP_R:
            return fShapeParError[ExternalToInternal(par)];
          case eX0:
          case eAsym:
          case eLambda:
          case eFWHM:
            return 0; // error propagation not yet implemented
          }
        break;
      }
    }

    ostringstream errMsg;
    errMsg << "no shape parameter available for "
           << GetShapeParameterName(par)
           << ", GH type " << GetFunctionTypeName(fFunctionType);
    throw DoesNotComputeException(errMsg.str());

  }

  double
  GaisserHillas4Parameter::GetCorrelationXMaxShapeParameter(const EShapeParameter par)
    const
  {
    return fRhoXMaxShapePar[ExternalToInternal(par)];
  }

  double
  GaisserHillas4Parameter::GetCorrelationNMaxShapeParameter(const EShapeParameter par)
    const
  {
    return fRhoNMaxShapePar[ExternalToInternal(par)];
  }

  double
  GaisserHillas4Parameter::Inverse(const double h, const int branch)
    const
  {
    if (h <= 0 || h > 1) {
      ostringstream errMsg;
      errMsg << "Value "<< h <<" for fraction of NMax not in (0,1]";
      throw OutOfBoundException(errMsg.str());
    }

    const double x0 = GetShapeParameter(eX0);
    const double lambda = GetShapeParameter(eLambda);
    const double xmax = (fXMax - x0) / lambda;
    const double z = - pow(h, 1/xmax) / kE;
    const double xB = -xmax * (branch < 0 ? LambertW<-1>(z) :
                               LambertW<0>(z));
    return x0 + lambda * xB ;
  }

  double
  GaisserHillas4Parameter::InverseLeft(const double h)
    const
  {
    return Inverse(h, 0);
  }

  double
  GaisserHillas4Parameter::InverseRight(const double h)
    const
  {
    return Inverse(h, -1);
  }

  inline
  double
  CalcFfromR(const double R)
  {
    // see Eq.(6) of width paper
    return (exp(-R) * (1+R) - 1) / (R * (exp(-R) - 1));
  }

  double
  GaisserHillas4Parameter::CalculateR(const double asym)
    const
  {

    const double asymMin = 0.01;
    const double asymMax = 0.5;

    if (asym >= asymMin && asym <= asymMax) {

      // fill table for R vs. asym
      if (!fRvsAsymTable.GetNPoints()) {
        vector<double> rVal;
        vector<double> fVal;

        // initial values at R = 0
        double currR = 0;
        double currAsym = 0.5;
        rVal.push_back(currR);
        fVal.push_back(currAsym);

        // fine binning for steep part
        const double deltaRfine = 0.05;
        while (currAsym > 0.4) {
          currR += deltaRfine;
          currAsym = CalcFfromR(currR);
          rVal.push_back(currR);
          fVal.push_back(currAsym);
        }

        // coarse binning for flat part
        const double deltaRCoarse = 1;
        while (currAsym > asymMin) {
          currR += deltaRCoarse;
          currAsym = CalcFfromR(currR);
          rVal.push_back(currR);
          fVal.push_back(currAsym);
        }
        const int nPol = 2;
        fRvsAsymTable = TabulatedFunction(fVal, rVal, nPol);
      }

      return fRvsAsymTable.Y(asym);

    }
    else {
      ostringstream errMsg;
      errMsg << "Asymmetry value " << asym <<" out of range [" << asymMin
             << "," << asymMax << "]";
      throw OutOfBoundException(errMsg.str());
    }
  }

  GaisserHillas4Parameter::EInternalShapeParameter
  GaisserHillas4Parameter::ExternalToInternal(const EShapeParameter par)
    const
  {
    switch (fFunctionType) {
    case eClassic: {
      switch (par) {
      case eX0:
        return eFirst;
      case eLambda:
        return eLast;
      default:
        break;
      }
      break;
    }
    case eWidth:
      {
        switch (par) {
        case eFWHM:
          return eFirst;
        case eAsym:
          return eLast;
        default:
          break;
        }
        break;
      }
    case eUSP:
      {
        switch (par) {
        case eUSP_L:
          return eFirst;
        case eUSP_R:
          return eLast;
        default:
          break;
        }
        break;
      }
    }
    ostringstream errMsg;
    errMsg << "parameter " << gh::GetShapeParameterName(par)
           << " is not a shape parameter of GH type "
           << GetFunctionTypeName(fFunctionType);
    throw DoesNotComputeException(errMsg.str());
  }

  EShapeParameter
  GaisserHillas4Parameter::InternalToExternal(const EInternalShapeParameter par)
    const
  {
    return InternalToExternal(fFunctionType, par);
  }

  EShapeParameter
  GaisserHillas4Parameter::InternalToExternal(const EFunctionType type,
                                              const EInternalShapeParameter par)
  {
    switch (type) {
    case eClassic:
      {
        switch (par) {
        case eFirst:
          return eX0;
        case eLast:
          return eLambda;
        }
        break;
      }
    case eWidth:
      {
        switch (par) {
        case eFirst:
          return eFWHM;
        case eLast:
          return eAsym;
        }
        break;
      }
    case eUSP:
      {
        switch (par) {
        case eFirst:
          return eUSP_L;
        case eLast:
          return eUSP_R;
        }
        break;
      }
    }

    ostringstream errMsg;
    errMsg << "no external shape parameter for par =  " << par
           << ", GH type " << GetFunctionTypeName(type);
    throw DoesNotComputeException(errMsg.str());

  }

  void
  GaisserHillas4Parameter::SetCorrelationXMaxShapeParameter(const EShapeParameter par,
                                                            const double rho)
  {
    fRhoXMaxShapePar[ExternalToInternal(par)] = rho;
  }

  void
  GaisserHillas4Parameter::SetCorrelationNMaxShapeParameter(const EShapeParameter par,
                                                            const double rho)
  {
    fRhoNMaxShapePar[ExternalToInternal(par)] = rho;
  }

  void
  GaisserHillas4Parameter::SetShapeParameter(const EShapeParameter par,
                                             const double value, const double error)
  {
    fShapePar[ExternalToInternal(par)] = value;
    fShapeParError[ExternalToInternal(par)] = error;
  }


  bool
  GaisserHillas4Parameter::IsInternal(const EShapeParameter par)
    const
  {
    switch (fFunctionType) {
    case eClassic:
      {
        switch (par)
          {
          case eX0:
          case eLambda:
            return true;
          default:
            return false;
          }
        break;
      }
    case eWidth:
      {
        switch (par)
          {
          case eFWHM:
          case eAsym:
            return true;
          default:
            return false;
          }
        break;
      }
    case eUSP:
      {
        switch (par)
          {
          case eUSP_L:
          case eUSP_R:
            return true;
          default:
            return false;
          }
        break;
      }
    }

    ostringstream errMsg;
    errMsg << "no internal shape parameter defined for "
           << GetShapeParameterName(par)
           << ", GH type " << GetFunctionTypeName(fFunctionType);
    throw DoesNotComputeException(errMsg.str());
  }

  // -----------------------------------------------------------
  // below this point only deprecated functions
  double
  GaisserHillas4Parameter::GetXZero()
    const
  {
    Deprecator::GetInstance().Deprecated(
                                         "evt::GaisserHillas4Parameter::GetXZero()",
                                         "v2r9",
                                         "GaisserHillas4Parameter::GetShapeParameter(gh::eX0)");
    return GetShapeParameter(eX0);
  }

  double
  GaisserHillas4Parameter::GetXZeroError()
    const
  {
    Deprecator::GetInstance().Deprecated(
                                         "evt::GaisserHillas4Parameter::GetXZeroError()",
                                         "v2r9",
                                         "GaisserHillas4Parameter::GetShapeParameterError(gh::eX0)");
    return GetShapeParameterError(eX0);
  }

  double
  GaisserHillas4Parameter::GetLambda()
    const
  {
    Deprecator::GetInstance().Deprecated(
                                         "evt::GaisserHillas4Parameter::GetLambda()",
                                         "v2r9",
                                         "GaisserHillas4Parameter::GetShapeParameter(gh::eLambda)");
    return GetShapeParameter(eLambda);
  }

  double
  GaisserHillas4Parameter::GetLambdaError()
    const
  {
    Deprecator::GetInstance().Deprecated(
                                         "evt::GaisserHillas4Parameter::GetLambdaError()",
                                         "v2r9",
                                         "GaisserHillas4Parameter::GetShapeParameterError(gh::eLambda)");
    return GetShapeParameterError(eLambda);
  }

  double
  GaisserHillas4Parameter::GetNMaxLambdaCorrelation()
    const
  {
    Deprecator::GetInstance().Deprecated(
                                         "evt::GaisserHillas4Parameter::GetNMaxLambdaCorrelation()",
                                         "v2r9",
                                         "GaisserHillas4Parameter::GetCorrelationXMaxShapeParameter(gh::eLambda)");
    return GetCorrelationNMaxShapeParameter(eLambda);
  }

  double
  GaisserHillas4Parameter::GetNMaxX0Correlation()
    const
  {
    Deprecator::GetInstance().Deprecated(
                                         "evt::GaisserHillas4Parameter::GetNMaxX0Correlation()",
                                         "v2r9",
                                         "GaisserHillas4Parameter::GetCorrelationXMaxShapeParameter(gh::eX0)");
    return GetCorrelationNMaxShapeParameter(eX0);
  }

  double
  GaisserHillas4Parameter::GetXMaxLambdaCorrelation()
    const
  {
    Deprecator::GetInstance().Deprecated(
                                         "evt::GaisserHillas4Parameter::GetXMaxLambdaCorrelation()",
                                         "v2r9",
                                         "GaisserHillas4Parameter::GetCorrelationXMaxShapeParameter(gh::eLambda)");
    return GetCorrelationXMaxShapeParameter(eLambda);
  }

  double
  GaisserHillas4Parameter::GetXMaxX0Correlation()
    const
  {
    Deprecator::GetInstance().Deprecated(
                                         "evt::GaisserHillas4Parameter::GetXMaxX0Correlation()",
                                         "v2r9",
                                         "GaisserHillas4Parameter::GetCorrelationXMaxShapeParameter(gh::eX0)");
    return GetCorrelationXMaxShapeParameter(eX0);
  }

  double
    GaisserHillas4Parameter::GetLambdaX0Correlation()
    const
  {
    Deprecator::GetInstance().Deprecated(
    "evt::GaisserHillas4Parameter::GetLambdaX0Correlation()",
    "v2r9",
    "GaisserHillas4Parameter::GetCorrelationShapeParameters()");
    return GetCorrelationShapeParameters();
  }
}