PhysicalFunctions.cc

Go to the documentation of this file.

#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>

using namespace std;
using namespace utl;


// For some reason, CLHEP SystemOfUnits.h #defines pascal to hep_pascal,
// putting it in the global namespace and screwing up our own definition of
// pascal from AugerUnits.h.  It appears the CLHEP units are being pulled in
// via the geometry package.  The following line undoes the CLHEP sabotage.
#undef pascal


namespace utl {

  inline
  double
  NormalizedGaisserHillas(const double x, const double xMax)
  {
    return pow(x / xMax, xMax) * exp(xMax - x);
  }


  double
  GaisserHillas(const double x, const double x0,
                const double xMax, const double nMax,
                const double lambda)
  {
    if (x < x0)
      return 0;

    return nMax *
      NormalizedGaisserHillas((x - x0) / lambda,
                              (xMax - x0) / lambda);
  }

  namespace RefractionIndex {
    
    double
    LorentzLorentz(const double verticalDepth)
    {
      const double n0 = kRefractiveIndexSeaLevel;
      const double x0 = kOverburdenSeaLevel;

      // Where the Lorentz-Lorenz relation is assumed:
      // (n^2 - 1) / (n^2 + 2) = constant*density

      const double n02 = Sqr(n0);
      const double dx = verticalDepth * (n02 - 1) / (n02 + 2) / x0;

      return sqrt((1 + 2 * dx) / (1 - dx));
    }

    double
    GladstoneDale(const double density, 
                  const double densityAtSeaLevel/* = kAirDensitySeaLevel*/,
                  const double refractiveIndexAtSeaLevel/* = kRefractiveIndexSeaLevel*/)
    {
      return 1 + (refractiveIndexAtSeaLevel - 1) * density / densityAtSeaLevel;
    }

    double
    Ciddor95(const double wl, const double temperature,
             const double pressure, const double vaporPressure)
    {
      // Refraction index: P. Ciddor, Appl Opt 35 (1996) 1566-1573
      const double invWl2 = pow(wl / micrometer, -2);
      const double tempC = temperature - kH2OFreezingPoint;
      const double tempK = temperature / kelvin;
      const double presPa = pressure / pascal;
      constexpr double R = kMolarGasConstant / (joule / (mole * kelvin));

      // Calculate the enhancement factor and mole fraction (xV) of water vapor
      const double fpt = 1.00062 + 3.14e-8 * presPa + 5.6e-7 * tempC * tempC;
      const double xV = fpt * vaporPressure / pressure;
      constexpr double Mv = kH2OMolarMass / (kg / mole);

      // Molar mass of dry air [kg/mol] with a CO2 concentration of xCO2 [ppm]
      constexpr double xCO2 = kCO2AirFraction / perMillion;
      constexpr double Ma = kDryAirMolarMass / (kg/mole) + 1.2011e-8 * (xCO2 - 400);

      // Calculate the component indices of refraction as r = 1e-8 * (n-1).
      // Standard dry air index (15 C, 101.325 kPa, 450 ppm CO2 or xCO2 ppm)
      const double rAS = 1e-8 * (5792105 / (238.0185 - invWl2) +
                                167917 / (57.362 - invWl2));
      const double rAXS = rAS * (1 + 5.34e-7 * (xCO2 - 450));

      // Standard index and density of pure water vapor (xV = 1) at 20 C, 1333 Pa
      const double rVS = 1.022e-8 * (
        295.235 + invWl2 * (2.6422 + invWl2 * (-0.03238 + invWl2 * 0.004028)));
      constexpr double rhoVS = 0.00985938;

      // Compressibility and density of dry air
      constexpr double pa = 101325;
      constexpr double Ta = 288.15;
      constexpr double Za = 0.9995922115;
      constexpr double rhoAXS = pa * Ma / (Za * R * Ta);

      // Compressibility of moist air
      const double pOnT = presPa / tempK;
      const double Zm =
        1 - pOnT * (1.58123e-6 + tempC * (-2.9311e-8 + 1.1043e-10 * tempC) +
            xV * ((5.707e-6 - 2.051e-8*tempC) + (1.9898e-4 - 2.376e-6*tempC)*xV))
          + pOnT * pOnT * (1.93e-11 - 0.765e-8 * xV * xV);

      // Density of water vapor and moist air
      const double rhoV = xV * presPa * Mv / (Zm * R * tempK);
      const double rhoA = (1 - xV) * presPa * Ma / (Zm * R * tempK);

      return 1 + (rhoA / rhoAXS) * rAXS + (rhoV / rhoVS) * rVS;
    }
  }

  double
  SaturationVaporPressure(const double temperature)
  {
    // Use the expression due to P.H. Huang, Proc. Third Int. Symp. Humidity
    // and Moisture 1 (1998) 69-76 (see also http://emtoolbox.nist.gov).
    double vaporPressure = 0;
    const double T = temperature / kelvin;

    if (T >= 273.15) {
      constexpr double k1 =  1.16705214528e+03;
      constexpr double k2 = -7.24213167032e+05;
      constexpr double k3 = -1.70738469401e+01;
      constexpr double k4 =  1.20208247025e+04;
      constexpr double k5 = -3.23255503223e+06;
      constexpr double k6 =  1.49151086135e+01;
      constexpr double k7 = -4.82326573616e+03;
      constexpr double k8 =  4.05113405421e+05;
      constexpr double k9 = -2.38555575678e-01;
      constexpr double k10 = 6.50175348448e+02;

      const double O = T + k9 / (T - k10);
      const double O2 = O*O;
      const double A = O2 + k1 * O + k2;
      const double B = k3 * O2 + k4 * O + k5;
      const double C = k6 * O2 + k7 * O + k8;
      const double X = -B + sqrt(B*B - 4*A*C);

      vaporPressure = 1e6 * pow(2*C/X, 4) * pascal;
    } else {
      const double A = -13.928169;
      const double B =  34.7078238;
      const double O = T / 273.16;
      const double Y =
        A * (1 - pow(O, -1.5)) + B * (1 - pow(O, -1.25));

      vaporPressure = 611.657 * exp(Y) * pascal;
    }

    return vaporPressure;
  }


  double
  MoliereRadius(double temperature, double pressure, const double cosTheta)
  {
//#warning MoliereRadius() does not work for inclined events
//#warning MoliereRadius() should use density from Atmosphere
    const double exponent = 1 / 5.25588;

    double correctedPressure = pressure - kgEarth * 2 * kRadiationLength * cosTheta;

    temperature /= kelvin;
    pressure /= millibar;
    correctedPressure /= millibar;

    // for security. High in the atmosphere, low pressure
    // linear asymptotic behaviour to zero
    const double p_crit = 10;  // [millibar]
    if (correctedPressure < p_crit) {
      const double dp = pressure - correctedPressure;
      correctedPressure = pressure - dp * pressure / (p_crit + dp); // should never get zero
    }

    const double rm = 272.5 * temperature *
      pow(correctedPressure/pressure, exponent) / correctedPressure;

    return rm * meter;
  }


  double
  ElectronsAboveCut(double enCut)
  {
    constexpr double maxMeVCut = 3*MeV; // maximum extrapolation

    if (enCut > maxMeVCut) {
      ostringstream warn;
      warn << "enCut " << enCut/MeV << " is larger than " << maxMeVCut/MeV << " MeV";
      WARNING(warn);
      enCut = maxMeVCut;
    }

    return 1 - 0.045*enCut/MeV;
  }


  double
  EnergyDeposit(const double age, const double enCut)
  {
    // parameters at 1 MeV
    constexpr double p[] = { 3.90883, 1.05301, 9.91717, 2.41715, 0.13180 };

    // limit to reasonable ranges
    const double scaleFactor = ElectronsAboveCut(enCut) / ElectronsAboveCut(1*MeV);

    const double eDep = p[0]/pow(p[1] + age, p[2]) + p[3] + p[4]*age;

    return eDep / scaleFactor * MeV/(g/cm2);
  }


  /****************************************************
   *
   * inversion of 3-by-3 matrix A
   *
   * (returns false if matrix singular)
   *
   ****************************************************/
  bool
  Invert3x3(double a[3][3])
  {
    const double determinant =
        + a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) 
        - a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) 
        + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);

    const double small = 1e-30;
    if (fabs(determinant) < small)
    {
      ERROR("Invert3x3: Error-matrix singular");
      return false;
    }

    double b[3][3];

    b[0][0] = a[1][1] * a[2][2] - a[1][2] * a[2][1];
    b[1][0] = a[1][2] * a[2][0] - a[2][2] * a[1][0];
    b[2][0] = a[1][0] * a[2][1] - a[1][1] * a[2][0];

    b[0][1] = a[0][2] * a[2][1] - a[2][2] * a[0][1];
    b[1][1] = a[0][0] * a[2][2] - a[2][0] * a[0][2];
    b[2][1] = a[0][1] * a[2][0] - a[0][0] * a[2][1];

    b[0][2] = a[0][1] * a[1][2] - a[1][1] * a[0][2];
    b[1][2] = a[0][2] * a[1][0] - a[1][2] * a[0][0];
    b[2][2] = a[0][0] * a[1][1] - a[0][1] * a[1][0];

    const double iDet = 1 / determinant;
    for (int i = 0; i < 3; ++i)
      for (int j = 0; j < 3; ++j)
        a[i][j] = iDet * b[i][j];

    return true;
  }


  /****************************************************
   *
   * solves linear system
   *
   *   / y[0] \   / A[0][0] A[0][1] A[0][2] \   / x[0] \
   *   | y[1] | = | A[1][0] A[1][1] A[1][2] | * | x[1] |
   *   \ y[2] /   \ A[2][0] A[2][1] A[2][2] /   \ x[2] /
   *
   *
   * Input:  y[3] and A[3][3]
   * Output: returns true when succeded (i.e. A is not singular)
   *         A is overwritten with its inverse
   *
   * M. Unger 12/1/05
   *
   *****************************************************/
  bool
  Solve3x3(const double y[3], double a[3][3], double x[3])
  {
    if (Invert3x3(a))
    {
      for (int i = 0; i < 3; ++i)
        x[i] = a[i][0] * y[0] + a[i][1] * y[1] + a[i][2] * y[2];
      return true;
    }
    return false;
  }


  bool
  QuadraticMaximumInterpolation(const std::vector<double>& x, const std::vector<double>& y,
                                double& xMax, double& yMax)
  {
    yMax = 0;
    xMax = 0;

    const unsigned int nX = x.size();

    if (nX != y.size() || nX < 5)
      return false;

    unsigned int iMax = 0;
    double tempYMax = y[iMax];

    for (unsigned int i = 1; i < nX; ++i) {
      if (y[i] > tempYMax) {
        iMax = i;
        tempYMax = y[i];
      }
    }

    if (iMax < 2 || iMax > nX - 2) {
      ERROR("QuadraticMaximumInterpolation - Error: Maximum not confined!");
      return false;
    }

    const double yin[3] = {y[iMax - 1], y[iMax], y[iMax + 1]};

    const double meanAroundMax = (yin[0] + yin[1] + yin[2]) / 3;
    const double threshold = 2;

    if (yin[1]/meanAroundMax > threshold || yin[0] < 0 || yin[2] < 0) {
      ERROR("QuadraticMaximumInterpolation - Error: Interpolation unreliable due to spike in data!");
      return false;
    }

    // quadratic "fit" around maximum to get yMax and Xmax
    const double xin[3] = {x[iMax - 1], x[iMax], x[iMax + 1]};
    double m[3][3];
    m[0][0] = Sqr(xin[0]);
    m[0][1] = xin[0];
    m[0][2] = 1;
    m[1][0] = Sqr(xin[1]);
    m[1][1] = xin[1];
    m[1][2] = 1;
    m[2][0] = Sqr(xin[2]);
    m[2][1] = xin[2];
    m[2][2] = 1;

    double a[3];

    Solve3x3(yin, m, a);

    if (a[0] <= 0) {
      xMax = -a[1] / (2 * a[0]);
      yMax = (a[0] * xMax + a[1]) * xMax + a[2];
    } else {
      ERROR("QuadraticMaximumInterpolation -  Error: no maximum found ...");
      return false;
    }

    return true;
  }


  namespace wcd {

    double
    SignalUncertainty(const ESignalVarianceModel model,
                      const double cosTheta, const double signal)
    {
      switch (model) {
      case eGAP2012_012:
        return SignalUncertaintyFactor(model, cosTheta) * sqrt(signal);
      case eGAP2014_035:
        {
          const double shFluctuations = SignalUncertaintyFactor(model, cosTheta) * sqrt(signal);
          const double detResolution = 0.023*signal;
          return sqrt(Sqr(shFluctuations) + Sqr(detResolution));
        }
      }
      ERROR("You should specify a SignalVariance model!");
      return 0;
    }


    double
    SignalUncertaintyFactor(const ESignalVarianceModel model, const double cosTheta)
    {
      switch (model) {
      case eGAP2012_012:
        return 0.34 + 0.46 / cosTheta;
      case eGAP2014_035:
        {
          constexpr double a = 0.865;
          constexpr double b = 0.593;
          const double factor = a*(1 + b*(1/cosTheta - 1.22));
          return factor;
        }
      }
      ERROR("You should specify a SignalVariance model!");
      return 0;
    }


    inline
    double
    ModeGauss(const double x, const double width)
    {
      return exp(-Sqr(x/width));
    }


    double
    TriggerProbability(const bool totdMoPSEnabled,
                       const double lgExpectedSignal, const double sin2Theta)
    {
      if (totdMoPSEnabled) {
        // suitable for the new triggers, calculated using doublet stations in infill array
        // see GAP-2018-017
        const double n = EvalPoly({0.091, -0.252}, sin2Theta);
        constexpr double sigma = 0.097;
        const double mu = EvalPoly({0.261, 0.174}, sin2Theta);
        const double x0 = EvalPoly({0.173, 0.170, -0.246}, sin2Theta);
        constexpr double delta = 0.225;
        const double gauss = n * ModeGauss(lgExpectedSignal - mu, sigma);
        const double sigmoid = Fermi(x0, lgExpectedSignal, delta);
        return gauss + sigmoid;
      }

      // based on a study using the new triggers in the regular array
      const double n = EvalPoly({0.163, -0.404}, sin2Theta);
      constexpr double sigma = 0.124;
      const double mu = EvalPoly({0.148, 0.294}, sin2Theta);
      const double x0 = EvalPoly({-0.012, 0.252, -0.148}, sin2Theta);
      constexpr double delta = 0.214;
      const double gauss = n * ModeGauss(lgExpectedSignal - mu, sigma);
      const double sigmoid = Fermi(x0, lgExpectedSignal, delta);
      return gauss + sigmoid;
    }

  }


  namespace ssd {

    double
    SignalUncertainty(const ESignalVarianceModel model,
                      const double cosTheta, const double signal)
    {
      switch (model) {
      case eGAP2018_025:
        return SignalUncertaintyFactor(model, cosTheta) * sqrt(signal);
      case eGAP2019_000:
        {
          const double shFluctuations = SignalUncertaintyFactor(model, cosTheta) * sqrt(signal);
          constexpr double baseFluctuations = 0.04;
          return sqrt(Sqr(shFluctuations) + Sqr(baseFluctuations));
        }
      }
      ERROR("You should specify a SignalVariance model!");
      return 0;
    }


    double
    SignalUncertaintyFactor(const ESignalVarianceModel model, const double cosTheta)
    {
      switch (model) {
      case eGAP2018_025:
        {
          constexpr double a = 1.398;
          constexpr double b = 0.210;
          const double factor = a*(1 + b*(1/cosTheta - 1.22));
          return factor;
        }
      case eGAP2019_000:
        {
          constexpr double a = 1.449;
          constexpr double b = 0.175;
          const double factor = a*(1 + b*(1/cosTheta - 1.22));
          return factor;
        }
      }
      ERROR("You should specify a SignalVariance model!");
      return 0;
    }

  }


  namespace InvisibleEnergy {

    constexpr double kQGSJetProton[] = { 0.958, 0.048, 0.162 };
    constexpr double kQGSJetIron[] = { 0.977, 0.109, 0.130 };
    constexpr double kQGSJetMixed[] = { 0.967, 0.078, 0.140 };

    constexpr double kQ2Sib[] = { 1.016, 1.012 };

    //constexpr double kData[] = { 1, 1.9065e-01, 8.9624e-01 };

    // ICRC 2013
    //constexpr double kData[] = { 1, 1.7425e-01, 9.1445e-01 };

    // ICRC 2017
    //constexpr double kData[] = {1, 0.1633, 0.9463, 0.8475, 0.957, 1.67};

    // PRD Paper
    //                           a      b      EcalA K  bextr  f1      f2     f3
    constexpr double kData[] = { 0.179, 0.947, 1.95, 4, 0.846, -0.265, 0.489, -0.441 };


    double
    EnergyLimits(const double energy)
    {
      constexpr double min = 0.01*EeV;
      constexpr double max = 1000*EeV;
      return std::min(std::max(min, energy), max);
    }


    double
    ModelFactor(const EInteractionModel interaction,
                const ECompositionModel composition)
    {
      switch (interaction) {
      case eQGSJET01:
      case eDATA:
        return 1;
      case eSIBYLL21:
        switch (composition) {
        case eProton:           return kQ2Sib[0];
        case eIron:             return kQ2Sib[1];
        case eData:             return 1;
        case eMixedComposition: return 0.5*(kQ2Sib[0] + kQ2Sib[1]);
        }
      }
      ERROR("You should specify a valid interaction and composition!");
      return 0;
    }


    const double*
    FitParameters(const ECompositionModel composition)
    {
      switch (composition) {
      case eData:             return kData;
      case eProton:           return kQGSJetProton;
      case eIron:             return kQGSJetIron;
      case eMixedComposition: return kQGSJetMixed;
      }
      ERROR("You should specify a valid composition!");
      return nullptr;
    }


    double
    Factor(const double energyEM,
           const EInteractionModel intMod,
           const ECompositionModel compMod,
           const double cosTheta)
    {
      const double energyEeV = InvisibleEnergy::EnergyLimits(energyEM) / EeV;

      const double mf = InvisibleEnergy::ModelFactor(intMod, compMod);

      const double* const p = InvisibleEnergy::FitParameters(compMod);

      if (intMod == eDATA) {
        //const double factor = mf * (p[0] + p[1] * pow(energyEeV, p[2]-1));

        // ICRC 2017
        //if (energyEeV > p[5])
        //  return mf * (p[0] + p[1] * pow(energyEeV, p[2]-1) * p[4]);
        //else (becouse the previous returns!)
        //return mf * (p[0] + p[1] * pow(p[5], p[2]) * pow(energyEeV/p[5], p[3]-1) /p[5] * p[4]);

        constexpr double th0 = 66*degree;
        const double cth = cosTheta - cos(th0);
        const double eCal = energyEeV;
        const double logECal = log10(energyEeV);
        const double a = p[0];
        const double b = p[1];
        const double eCalA = p[2];
        const double logECalA = log10(p[2]);
        const double k = p[3];
        const double bExt = p[4];
        const double fTheta = 1 + cth*(p[5] + cth*(p[6] + cth*p[7]));
        const double eInvH = a * pow(eCal, b);
        const double eInvL = a * pow(eCalA, b) * pow(eCal / eCalA, bExt);
        const double fTanh = 0.5 * (1 + tanh(k*(logECal - logECalA)));
        const double factor = 1 + (fTheta / eCal) * (eInvL + fTanh * (eInvH - eInvL));

        return factor;
      }

      return 1 / (mf * (p[0] - p[1] * pow(energyEeV, -p[2])));
    }


    double
    FactorDerivative(const double energyEM,
                     const EInteractionModel intMod,
                     const ECompositionModel compMod,
                     const double cosTheta)
    {
      const double energyLim = InvisibleEnergy::EnergyLimits(energyEM);
      const double mf = InvisibleEnergy::ModelFactor(intMod, compMod);
      const double* const p = InvisibleEnergy::FitParameters(compMod);
      const double ep2 = pow(energyLim/EeV, -p[2]);
      const double energyEeV = energyLim/EeV;

      if (intMod == eDATA) {
        //const double dfdE = (p[2] - 1) * p[1] * pow(energyLim/EeV, p[2]-1)/energyLim;

        // ICRC 2017
        //if (energyLim/EeV > p[5])
        //  return (p[2] - 1) * p[1] * pow(energyLim/EeV, p[2]-1) * p[4] / energyLim;

        constexpr double th0 = 66*degree;
        const double cth = cosTheta - cos(th0);
        const double eCal = energyEeV;
        const double logECal = log10(energyEeV);
        const double a = p[0];
        const double b = p[1];
        const double eCalA = p[2];
        const double logECalA = log10(p[2]);
        const double k = p[3];
        const double bExt = p[4];
        const double fTheta = 1 + cth*(p[5] + cth*(p[6] + cth*p[7]));
        const double eInvH = a * pow(eCal, b);
        const double eInvL = a * pow(eCalA, b) * pow(eCal / eCalA, bExt);
        const double fTanh = 0.5 * (1 + tanh(k*(logECal - logECalA)));
        const double fCosh = Sqr(cosh(k*(logECal - logECalA)));
        const double dFactor =
          fTheta / (eCal * energyLim) *
          ((bExt - 1)*eInvL +  fTanh*(((b - 1)*eInvH) - ((bExt - 1)*eInvL)) + (0.5/log(10))*k/fCosh * (eInvH - eInvL));

        return dFactor;
      }

      return -ep2 * p[1] * p[2] / (energyLim * mf * Sqr(p[0] - p[1]*ep2));
    }


    // ICRC 2013
    /*double
    FactorVariance(const double eTot)
    {
      const double p0 =  1.66299e+06;
      const double p1 = -6.35739e+00;
      const double uncert = eTot * p0 * pow(log10(eTot), p1);
      return Sqr(uncert);
    }*/


    // ICRC 2019
    double
    FactorVariance(const double eCal, const double eTot)
    {
      const double lgETot = log10(eTot);
      // Matias Tueros QGSJET0III p/Fe (50%/50%) composition
      // shower-to-shower fluctuations
      constexpr double pp0 = 1.20836e6;
      constexpr double pp1 = -6.42338;
      const double uncert1 = eTot * pp0 * pow(lgETot, pp1); // dEinv_1

      // get the uncertainty due to the mixed composition
      // Einv ~ A^(1-beta)
      // dEinv = Einv (1-beta) sqrt( Var(ln(A) )

      // beta from Conex EPOSLHC (from Einv Fe/p)
      constexpr double c0 = 1.558;
      constexpr double c1 = -0.0797843;
      constexpr double c2 = 0.00237286;
      const double fbeta = max(0., c0 + lgETot*(c1 + c2*lgETot));

      // lnA variance - EPOSLHC ICRC 2017 v2 - function from Michael U.
      constexpr double p0 = 1.52404;
      constexpr double p1 = 2.06276e2;
      constexpr double p2 = -2.19605e2;
      // not extrapolate below 0.9*1e17
      const double lgFeTot = log10(max(0.9*1e17, eTot)) / 18;
      const double var_lnA = p0 * pow(lgFeTot, p1 + p2*lgFeTot);
      const double sqrt_var_lnA = sqrt(var_lnA);

      const double uncert2 = (eTot - eCal) * (1 - fbeta) * sqrt_var_lnA; // dEinv_2

      const double uncertsquare = (Sqr(uncert1/eTot) + Sqr(uncert2/eTot)) * Sqr(eTot);

      return uncertsquare;
    }

  }


  namespace XmaxParam {
    
    double 
    Mean(const double energy, const double massNumber,
         const HadronicInteractionModel hadModel) 
    {
      const double lgE = log10(energy / utl::eV) - 19;
      const double lnA = log(massNumber);
      
      double x0 = 0;
      double d = 0;
      double xi = 0;
      double delta = 0;

      if (hadModel == HadronicInteractionModel::eSibyll23d) {
        x0 = 815.87 * g / cm2;
        d = 57.873 * g / cm2;
        xi = -0.3035 * g / cm2;
        delta = 0.7963 * g / cm2;
      } else if (hadModel == HadronicInteractionModel::eEPOSLHC) {
        x0 = 806.04 * g / cm2;
        d = 56.295 * g / cm2;
        xi = 0.3463 * g / cm2;
        delta = 1.0442 * g / cm2;
      } else if (hadModel == HadronicInteractionModel::eQGSJETII04) {
        x0 = 790.15 * g / cm2;
        d = 54.191 * g / cm2;
        xi = -0.4170 * g / cm2;
        delta = 0.6927 * g / cm2;
      } else {
        throw InvalidConfigurationException(
            "Invalid hadronic interaction model requested. " 
            "Implemented are Sibyll2.3d, EPOS-LHC and QGSJETII-04");
      }

      const double meanXmaxP = x0 + d * lgE;
      const double factorEnergy = xi - d / log(10) + delta * lgE   ;    
      return meanXmaxP + factorEnergy * lnA;
    }

    double
    StandardDeviation(const double energy, const double massNumber, 
                      const HadronicInteractionModel hadModel)
    {
      const double lgE = log10(energy / utl::eV) - 19;
      const double lnA = log(massNumber);
      double p0 = 0;
      double p1 = 0;
      double p2 = 0;
      double a0 = 0;
      double a1 = 0;
      double b = 0;

      if (hadModel == HadronicInteractionModel::eSibyll23d) {
        p0 = 3.727e3 * g * g / cm2 / cm2;
        p1 = -4.838e2 * g * g / cm2 / cm2;
        p2 = 1.325e2 * g * g / cm2 / cm2;
        a0 = -4.055e-1;
        a1 = -2.536e-4;
        b = 4.749e-2;
      } else if (hadModel == HadronicInteractionModel::eEPOSLHC) {
        p0 = 3.269e3 * g * g / cm2 / cm2;
        p1 = -3.053e2 * g * g / cm2 / cm2;
        p2 = 1.239e2 * g * g / cm2 / cm2;
        a0 = -4.605e-1;
        a1 = -4.751e-4;
        b = 5.838e-2;
      } else if (hadModel == HadronicInteractionModel::eQGSJETII04) {
        p0 = 3.688e3 * g * g / cm2 / cm2;
        p1 = -4.279e2 * g * g / cm2 / cm2;
        p2 = 5.096e1 * g * g / cm2 / cm2;
        a0 = -3.939e-1;
        a1 = 1.314e-3;
        b = 4.512e-2;
      } else {
        throw InvalidConfigurationException(
            "Invalid hadronic interaction model requested. " 
            "Implemented are Sibyll2.3d, EPOS-LHC and QGSJETII-04");
      }  

      double sigmaP2 = p0 + p1 * lgE + p2 * pow(lgE, 2);
      double a = a0 + a1 * lgE;

      return sqrt(sigmaP2 * (1 + a * lnA + b * pow(lnA, 2))); 
    }
  
  }


  namespace GumbelXmax {
    
    double
    Mu(const double energy, const double massNumber, 
         const HadronicInteractionModel hadModel)
    {
      const double lgE = log10(energy / utl::eV) - 19;
      const double lnA = log(massNumber);
      double p0 = 0;
      double p1 = 0;
      double p2 = 0;

      
      if (hadModel == HadronicInteractionModel::eSibyll23d) {
        p0 = 785.852 + lnA * (-15.5994 - 1.06906 * lnA);
        p1 = 60.5929 + lnA * (-0.786014 + 0.200728 * lnA);
        p2 = -0.689462 + lnA * (-0.294794 - 0.0399432 * lnA);
      } else if (hadModel == HadronicInteractionModel::eEPOSLHC) {
        p0 = 775.457 + lnA * (-10.3991 - 1.75261 * lnA);
        p1 = 58.5306 + lnA * (-0.827668 + 0.231144 * lnA);
        p2 = -1.40781 + lnA * (0.225624 - 0.10008 * lnA);
      } else if (hadModel == HadronicInteractionModel::eQGSJETII04) {
        p0 = 758.65 + lnA * (-12.3571 - 1.24539 * lnA);
        p1 =  56.5943 + lnA * (-1.01244 + 0.228689 * lnA);
        p2 = -0.534683 + lnA * (-0.17284 - 0.019159 * lnA);
      } else {
        throw InvalidConfigurationException(
            "Invalid hadronic interaction model requested. " 
            "Implemented are Sibyll2.3d, EPOS-LHC and QGSJETII-04");
      }

      const double mu = p0 + lgE * (p1 + p2 * lgE);
      return mu * g/cm2;
    }

    double
    Sigma(const double energy, const double massNumber, 
                      const HadronicInteractionModel hadModel)
    {
      const double lgE = log10(energy / utl::eV) - 19;
      const double lnA = log(massNumber);
      double p0 = 0;
      double p1 = 0;

      if (hadModel == HadronicInteractionModel::eSibyll23d) {
        p0 = 41.0345 + lnA * (-2.17329 - 0.306202 * lnA);
        p1 = -0.309466 + lnA * (-1.1649 + 0.225445 * lnA);
      } else if (hadModel == HadronicInteractionModel::eEPOSLHC) {
        p0 = 32.2632 + lnA * (3.94252 - 0.864421 * lnA);
        p1 = 1.27601 + lnA * (-1.81337 + 0.231914 * lnA);
      } else if (hadModel == HadronicInteractionModel::eQGSJETII04) {
        p0 = 35.4234 + lnA * (6.75921 - 1.46182 * lnA);
        p1 = -0.796042 + lnA * (0.201762 - 0.0142452 * lnA);
      } else {
        throw InvalidConfigurationException(
            "Invalid hadronic interaction model requested. " 
            "Implemented are Sibyll2.3d, EPOS-LHC and QGSJETII-04");
      }

      const double sigma = p0 + lgE * p1;
      return sigma * g/cm2;
    }


    double
    Lambda(const double energy, const double massNumber,
                      const HadronicInteractionModel hadModel)
    {
      const double lgE = log10(energy / utl::eV) - 19;
      const double lnA = log(massNumber);
      double p0 = 0;
      double p1 = 0;

      if (hadModel == HadronicInteractionModel::eSibyll23d) {
        p0 = 0.799493 + lnA * (0.235235 + 0.00856884 * lnA);
        p1 = 0.0632135 + lnA * (-0.0012847 + 0.000330525 * lnA);
      } else if (hadModel == HadronicInteractionModel::eEPOSLHC) {
        p0 = 0.641093 + lnA * (0.219762 + 0.171124 * lnA);
        p1 = 0.0726131 + lnA * (0.0353188 - 0.0131158 * lnA);
      } else if (hadModel == HadronicInteractionModel::eQGSJETII04) {
        p0 = 0.671545 + lnA * (0.373902 + 0.075325 * lnA);
        p1 = 0.0304335 + lnA * (0.0473985 - 0.000564531 * lnA);
      } else {
        throw InvalidConfigurationException(
            "Invalid hadronic interaction model requested. " 
            "Implemented are Sibyll2.3d, EPOS-LHC and QGSJETII-04");
      }

      const double lambda = p0 + lgE * p1;
      return lambda * g/cm2;
    }

  }


  inline
  double
  GoraAParameter(const double age)
  {
    return (((5.151*age - 28.925)*age + 60.056)*age - 56.718)*age + 22.331;
  }


  inline
  double
  GoraBParameter(const double age)
  {
    return (-1.039*age + 2.251)*age + 0.676;
  }


  double
  GoraCDF(const double rStar, const double age)
  {
    const double a = GoraAParameter(age);
    const double b = GoraBParameter(age);
    return 1- pow(1 + a*rStar, -b);
  }


  double
  InverseGoraCDF(const double fraction, const double age)
  {
    const double a = GoraAParameter(age);
    const double b = GoraBParameter(age);
    return (pow(1 - fraction, -1/b) - 1) / a;
  }


  double
  GoraPDF(const double rStar, const double age)
  {
    const double a = GoraAParameter(age);
    const double b = GoraBParameter(age);
    return b * a * pow(1 + a*rStar, -b - 1);
  }

}