Arbeletche2021CherenkovAngularModel.cc

Go to the documentation of this file.

#include <cmath>
#include <algorithm>
#include <string>

#include <atm/Arbeletche2021CherenkovAngularModel.h>

#include <utl/AugerException.h>
#include <utl/MathConstants.h>
#include <utl/Integrator.h>

// make the math constants visible
using utl::kPi;
using utl::kPiOnTwo;

// ctor
atm::Arbeletche2021CherenkovAngularModel::Arbeletche2021CherenkovAngularModel(
  const double showerEnergyTeV,
  const double maxTheta
) :
  fMaxTheta(maxTheta)
{
  fOmegaEnergyDep = 4.513 * std::pow(showerEnergyTeV, -0.008843);
  fEpsilon = 0.009528 + 0.022552 * std::pow(showerEnergyTeV, -0.4207);
}



void
atm::Arbeletche2021CherenkovAngularModel::UpdateParameters(
  const double showerAge,
  const double refractiveIndex
) const
{
  // Truncate the values of shower age and refractive index within the following bounds:
  // 0.5 <= showerAge <= 1.5
  // refractiveIndex >= 1.0 + 3e-5
        const double nMinusOne = std::max(refractiveIndex, 1.0 + 3e-5) - 1.0;
        const double truncatedAge = std::max(0.5, std::min(showerAge, 1.5));
        const double logShowerAge = std::log(truncatedAge);
        
        // Compute the parametrized parameters
        fNu = 0.21155 * std::pow(nMinusOne, -0.16639) + 1.21803 * logShowerAge;
        fOmega1 = 1.0 / (fOmegaEnergyDep * std::pow(nMinusOne, 0.45092) -       0.058687 * logShowerAge);
        fOmega2 = fOmega1 * (1.0 -  1.0 / (0.90725 + 0.41722 * truncatedAge));
        
        // Compute thetaEm (maximum Cherenkov emission angle) and its sines and cosines
        const double v = 1. - 1./refractiveIndex;
        fThetaEm = v < 0.001 ? std::sqrt(2.*v)*(1.+(v/12.)*(1.+9.*v/40.)) : std::acos(1./refractiveIndex);
        fCosThetaEm = 1./refractiveIndex;
        fSinThetaEm = std::sqrt((1.-fCosThetaEm)*(1.+fCosThetaEm));
        
        // Compute the normalization
        fNormalization = UnnormalizedIntegral(0.0, fMaxTheta);
}



double
atm::Arbeletche2021CherenkovAngularModel::FunctionK (
  const double theta
) const
{
  return std::pow(theta,fNu-1.0) *
                (std::exp(-theta*fOmega1) + fEpsilon*std::exp(-theta*fOmega2));
}



double
atm::Arbeletche2021CherenkovAngularModel::UnnormalizedIntegral(
  const double lowAngle,
  const double highAngle
) const
{
  double integral = 0;
  
  // If lowAngle is to the left of the peak, integrate between lowAngle and min(thetaEm, highAngle)
  if (lowAngle < fThetaEm) {
    integral += UnnormalizedIntegralLeft(lowAngle, std::min(highAngle, fThetaEm));
        }
        
        // If highAngle is to the right of the peak, integrate between max(thetaEm,lowAngle) and highAngle
        if (highAngle > fThetaEm) {
          integral += UnnormalizedIntegralRight(std::max(lowAngle, fThetaEm), highAngle);
        }
        
  return integral;
}



double
atm::Arbeletche2021CherenkovAngularModel::UnnormalizedIntegralLeft(
  const double lowAngle,
  const double highAngle
) const
{
  // In this region of theta < thetaEm, we provide the integral as a series in powers of 
  // (theta - thetaEm) plus some terms resulting from integration by parts. Since
  // thetaEm < 0.025, no more than 7 terms are required for the computation of the series
  
  // Helper quantities
  const double deltaLow = fThetaEm-lowAngle;
  const double deltaHigh = fThetaEm-highAngle;
  const double deltaCosLow = std::cos(lowAngle) - fCosThetaEm;
  const double deltaCosHigh = std::cos(highAngle) - fCosThetaEm;
  const double coeffs[4] = {fCosThetaEm, fSinThetaEm, -fCosThetaEm, -fSinThetaEm};

  // The power series
  double integral = 0;
  double powLow = 1;
  double powHigh = 1;
  double term;
  bool hasConverged = false;
  for (int k = 1; k < 1000; ++k) {
    powLow *= deltaLow / k;
    powHigh *= deltaHigh / k;
    term = (powLow - powHigh) / k;
    integral += coeffs[k%4] * term;
    if (std::fabs(term) <= 1e-10*std::fabs(integral)) {
      hasConverged = true;
      break;
    }
  }

  // check if the series has actually converged
  if (!hasConverged) {
    const std::string message = "Failed to integrate the angular distribution of Cherenkov photons";
    ERROR(message);
    throw utl::DoesNotComputeException(message);
  }

  // Terms from integration by parts
  integral += kPi * (deltaCosLow - deltaCosHigh);
  integral += deltaHigh < 1e-10 ? 
    std::log(std::pow(deltaHigh/fThetaEm+1e-200,deltaCosHigh)) :
    deltaCosHigh * std::log(deltaHigh/fThetaEm); 
  integral -= deltaLow < 1e-10 ? 
    std::log(std::pow(deltaLow/fThetaEm+1e-200,deltaCosLow)) :
    deltaCosLow * std::log(deltaLow/fThetaEm);
  
  // A remaining factor
  integral *= FunctionK(fThetaEm) / fSinThetaEm;
    
  return integral;
}



double
atm::Arbeletche2021CherenkovAngularModel::UnnormalizedIntegralRight(
  const double lowAngle,
  const double highAngle
) const
{
  //
  // Helper functions (which are used below)
  //
  
  // Integrand (after the change t = sqrt(1-thetaEm/theta)) used for numerical integration
  auto integrand = [=](double t) {
    const double theta = fThetaEm / ((1.0-t)*(1.0+t));
    return std::pow(theta, fNu + 1) *
            (t <= 1e-10 ? t*kPiOnTwo + 1.0 - std::pow(t, t) : t * (kPiOnTwo - std::log(t))) *
            (std::exp(-fOmega1*theta) + fEpsilon*std::exp(-fOmega2*theta));
  };
  
  // Indefinite integral given as an expansion in powers of (1 - thetaEm/theta)
  auto uSeries = [=](const double theta, const double omega)->double{
    const double u = 1.0 - fThetaEm/theta;
    if (u < 1e-18) {
      return 0;
    }
  
    const double pimlog = u > 1e-5 ?
      u*(kPi - std::log(u)) :
      u*kPi - std::log( std::pow(u+1e-200, u) );
    
    double bm2 = 0;
    double bm1 = std::exp(-omega*fThetaEm);
    double sum = (pimlog+u)*bm1;
    double a1 = fNu-1-omega*fThetaEm;
    double a2 = -1;
    double upn = u;
    double b, term;
    for (int n = 2; n < 1000; ++n) {
      a1 += 2;
      ++a2;
      b = a1*bm1 - a2*(fNu+a2)*bm2*upn;
      upn = u/n;
      b *= upn;
      term = b * (pimlog + upn);
      sum += term;
      if (std::fabs(term) <= 1e-10*std::fabs(sum)) {
        return sum * std::pow(fThetaEm, fNu);
      }
      bm2 = bm1;
      bm1 = b;
    }

    // if we arrived here, the series has not converged
    throw;
  };
  
  // Indefinite integral given as an expansion in powers of thetaEm/theta
  auto rSeries = [=](const double theta, const double omega)->double{
    const double wt = omega*theta;
    const double r = fThetaEm / theta;
    const double wtem = fThetaEm * omega;
    
    double h = 0;
    
    // Compute IncGammaUpper(nu, x)/(t^nu e^(-x)) according
    // to Press et al., Numerical Recipes, 3ed, 2007
    // The quantity computed here corresponds to the first term of
    // the power series computed in the following loop (uses recursion)
    bool hasConverged = false;
    if (wt > fNu + 1) {
      // Use continued fraction expansion
      const double fpMin = 1e-30;
      double b = wt + 1 - fNu;
      double c = 1.0 / fpMin;
      double d = 1.0 / b;
      h = -d;
      for (int n = 1; n < 1000; ++n) {
        double an = -n * (n - fNu);
        b += 2.0;
        d = an * d + b;
        if (std::fabs(d) < fpMin) {
          d = fpMin;
        }
        c = b + an/c;
        if (std::fabs(c) < fpMin) {
          c = fpMin;
        }
        d = 1.0/d;
        double del = d*c;
        h *= del;
        if (std::fabs(del - 1) < 1e-10) {
          hasConverged = true;
          break;
        }
      }
    } else {
      // Use power series
      double term = 1.0/fNu;
      h = term;
      for (int n = 1; n < 1000; ++n) {
        term *= wt / (n + fNu);
        h += term;
        if (std::fabs(term) < 1e-10*std::fabs(h)) {
          hasConverged = true;
          break;
        }
      }
      h -= std::exp(wt - fNu*std::log(wt) + std::lgamma(fNu));
    }

    // only proceed if the calculation above has converged
    if (hasConverged) {
      double rPower = 1;
      double sum = kPi * h;
      for (int k = 1; k < 1000; ++k) {
        rPower *= r;
        h = (rPower + wtem*h) / (fNu - k);
        sum += h / double(k);
        if (std::fabs(h / double(k)) <= 1e-10*std::fabs(sum)) {
          return sum * std::exp(-wt + fNu * std::log(theta));
        }
      }
    }
    
    // if we arrived here, the series has not converged
    throw;
  };
  
  //
  // Computation of the integral
  //
  
  // Define a threshold used to select a power series expansion either in terms of
  // 1-thetaEm/theta (below threshold) and thetaEm/theta (above threshold)
  const double thresholdAngle = 2*fThetaEm;

  // Start the integral with a 0 value
  double integral = 0;
  
  // Region below the threshold angle: always use the expansion in powers of
  // 1-thetaEm/theta, which converges for every combination of parameters
  if (lowAngle < thresholdAngle) {
    const double actualHighAngle = std::min(highAngle, thresholdAngle);
    try {
      integral += uSeries(actualHighAngle, fOmega1) - uSeries(lowAngle, fOmega1) +
        fEpsilon * (uSeries(actualHighAngle, fOmega2) - uSeries(lowAngle, fOmega2));
    } catch (...) {
      const std::string message = "Failed to integrate the angular distribution of Cherenkov photons";
      ERROR(message);
      throw utl::DoesNotComputeException(message);
    }
  }

  // Region above the threshold angle: use the expansion in powers of thetaEm/theta
  // only if parameter Nu is not too close to an integer value, in which the series
  // may diverge due to a factor of 1/(nu-k). If nu is close to an integer value,
  // use numerical integration instead.
  if (highAngle > thresholdAngle) {
    const double actualLowAngle = std::max(lowAngle, thresholdAngle);
    if (std::fabs(fNu - std::round(fNu)) < 1e-4) {
      // Nu is close to integer value, use numerical integration
      const double tMin = std::sqrt(1.0 - fThetaEm/actualLowAngle);
      const double tMax = std::sqrt(1.0 - fThetaEm/highAngle);
      utl::Integrator<decltype(integrand)> integrator(integrand, 1e-6);
      integral += (4.0 / fThetaEm) * integrator.GetRombergIntegral(tMin, tMax);
    } else {
      // Use power series
      try {
        integral += rSeries(highAngle, fOmega1) - rSeries(actualLowAngle, fOmega1) +
          fEpsilon * (rSeries(highAngle, fOmega2) - rSeries(actualLowAngle, fOmega2));
      } catch (...) {
        const std::string message = "Failed to integrate the angular distribution of Cherenkov photons";
        ERROR(message);
        throw utl::DoesNotComputeException(message);
      }
    }
  }
  
        return integral;
}



double
atm::Arbeletche2021CherenkovAngularModel::UnnormalizedIntegralSineLeft(
  const double lowAngle,
  const double highAngle
) const
{
  // In this region of theta < thetaEm, we provide the integral as a series in powers of 
  // (theta - thetaEm) plus some terms resulting from integration by parts. Since
  // thetaEm < 0.025, no more than 7 terms are required for the computation of the series
  
  // Helper quantities
  const double deltaLow = 2.0*(fThetaEm - lowAngle);
  const double deltaHigh = 2.0*(fThetaEm - highAngle);
        const double twoThetaEm = 2.*fThetaEm;
        const double cosTwoThetaEm = 2.0*fCosThetaEm*fCosThetaEm - 1.0;
        const double sinTwoThetaEm = 2.0*fSinThetaEm*fCosThetaEm;
  const double factorLow = sinTwoThetaEm - std::sin(2.0*lowAngle) - deltaLow;
  const double factorHigh = sinTwoThetaEm - std::sin(2.0*highAngle) - deltaHigh;
  const double coeffs[4] = {sinTwoThetaEm, -cosTwoThetaEm, -sinTwoThetaEm, cosTwoThetaEm};
  
  // The integral value
  double integral = 0;
  
  // Compute the series part
  double powLow = 1.0;
  double powHigh = 1.0;
  double term;
  bool hasConverged = false;
  for (int k = 1; k < 1000; ++k) {
    powLow *= deltaLow / k;
    powHigh *= deltaHigh / k;
    term = (powLow - powHigh) / k;
    integral += coeffs[k%4] * term;
    if (std::fabs(term) <= 1e-10*std::fabs(integral)) {
      hasConverged = true;
      break;
    }
  }

  // Check convergence of the power series part
  if (!hasConverged) {
    const std::string message = "Failed to integrate the angular distribution of Cherenkov photons";
    ERROR(message);
    throw utl::DoesNotComputeException(message);
  }
  
  // Terms obtained from integration by parts
  integral += deltaLow - deltaHigh;
  integral += kPi * (factorHigh-factorLow);
  integral -= deltaHigh < 1e-10 ?
    std::log(std::pow(deltaHigh/twoThetaEm+1e-200,factorHigh)) :
    factorHigh * std::log(deltaHigh/twoThetaEm);
  integral += deltaLow < 1e-10 ?
    std::log(std::pow(deltaLow/twoThetaEm+1e-200,factorLow)) :
    factorLow * std::log(deltaLow/twoThetaEm);
  
  // A last factor depending on thetaEm alone
  integral *= 0.25 * FunctionK(fThetaEm) / fSinThetaEm;
  
  return integral;
}



double
atm::Arbeletche2021CherenkovAngularModel::UnnormalizedIntegralSineRight(
  const double lowAngle,
  const double highAngle
) const
{
  //
  // Helper functions (which are used below)
  //
  
  // Integrand (after the change t = sqrt(1-thetaEm/theta)) used for numerical integration
  auto integrand = [=](double t) {
          const double theta = fThetaEm / ((1.0-t)*(1.0+t));
          return std::sin(theta) * std::pow(theta, fNu + 1) *
                  (t <= 1e-10 ? t*kPiOnTwo + 1.0 - std::pow(t, t) : t * (kPiOnTwo - std::log(t))) *
                  (std::exp(-fOmega1*theta) + fEpsilon*std::exp(-fOmega2*theta));
  };
  
  // Indefinite integral given as an expansion in powers of (1 - thetaEm/theta)
  auto uSeries = [=](const double theta, const double omega)->double{
    const double u = 1.0 - fThetaEm/theta;
    if (u < 1e-18) {
      return 0;
    }
  
    const double pimlog = u > 1e-5 ?
      u*(kPi - std::log(u)) :
      u*kPi - std::log(std::pow(u+1e-200,u));
    
    double bm2 = 0;
    double am2 = 0;
    double bm1 = std::exp(-omega*fThetaEm);
    double am1 = bm1 * fCosThetaEm;
    bm1 *= fSinThetaEm;
    double sum = (pimlog + u) * bm1;
    double c1 = fNu-1-omega*fThetaEm;
    double c2 = -1;
    double upn = u;
    double a, b, term;
    for (int n = 2; n < 1000; ++n) {
      c1 += 2;
      ++c2;
      a = c1*am1 - c2*(fNu+c2)*am2*upn - fThetaEm*bm1;
      b = c1*bm1 - c2*(fNu+c2)*bm2*upn + fThetaEm*am1;
      upn = u/n;
      a *= upn;
      b *= upn;
      term = (pimlog + upn) * b;
      sum += term;
      if (std::fabs(term) <= 1e-10*std::fabs(sum)) {
        return sum * std::pow(fThetaEm, fNu);
      }
      bm2 = bm1;
      bm1 = b;
      am2 = am1;
      am1 = a;
    }

    // if the series has not converged
    throw;
  };
  
  // Indefinite integral given as an expansion in powers of thetaEm/theta
  auto rSeries = [=](const double theta, const double omega) {
    const double r = fThetaEm/theta;
    const double sint = std::sin(theta);
    const double cost = std::cos(theta);
    const double invt = 1.0/theta;
    
    // here we compute the first coefficient of the following power series,
    // which is also computed as a power series. the following terms are
    // obtained through recursion
    double kap = 1.0/fNu;
    double sig = 0;
    double old;
    double c = kap;
    double s = sig;
    double den = fNu * invt;
    bool hasConverged = false;
    for (int n = 1; n < 1000; ++n) {
      den += invt;
      old = kap;
      kap = (omega*kap - sig) / den;
      sig = (omega*sig + old) / den;
      c += kap;
      s += sig;
      if (std::fabs(kap) <= 1e-10*std::fabs(c) && std::fabs(sig) <= 1e-10*std::fabs(s)) {
        hasConverged = true;
        break;
      }
    }

    // ensure the computation above has converged
    if (hasConverged) {
      double b = sint*c - cost*s;
      double rp = 1;
      double series = kPi * b;
      den = fNu;
      for (int k = 1; k < 1000; ++k) {
        den--;
        rp *= r;
        old = s;
        s = (omega*s + c) * fThetaEm / den;
        c = (rp + fThetaEm*(omega*c - old)) / den;
        b = (sint*c - cost*s) / k;
        series += b;
        if (std::fabs(b) < 1e-10*std::fabs(series)) {
          return series * std::exp(-omega*theta + fNu*std::log(theta));
        }
      }
    }
    
    // if the series has not converged, throw
    throw;
  };

  //
  // Computation of the integral
  //
  
  // Define a threshold used to select a power series expansion either in terms of
  // 1-thetaEm/theta (below threshold) and thetaEm/theta (above threshold)
  const double thresholdAngle = 2*fThetaEm;

  // Start the integral with a 0 value
  double integral = 0;
  
  // Region below the threshold angle: always use the expansion in powers of
  // 1-thetaEm/theta, which converges for every combination of parameters
  if (lowAngle < thresholdAngle) {
    const double actualHighAngle = std::min(highAngle, thresholdAngle);
    try {
      integral += uSeries(actualHighAngle, fOmega1) - uSeries(lowAngle, fOmega1) +
        fEpsilon * (uSeries(actualHighAngle, fOmega2) - uSeries(lowAngle, fOmega2));
    } catch (...) {
      const std::string message = "Failed to integrate the angular distribution of Cherenkov photons";
      ERROR(message);
      throw utl::DoesNotComputeException(message);
    }
  }

  // Region above the threshold angle: use the expansion in powers of thetaEm/theta
  // only if parameter Nu is not too close to an integer value, in which the series
  // may diverge due to a factor of 1/(nu-k). If nu is close to an integer value,
  // use numerical integration instead.
  if (highAngle > thresholdAngle) {
    const double actualLowAngle = std::max(lowAngle, thresholdAngle);
    if (std::fabs(fNu - std::round(fNu)) < 1e-4) {
      // Nu is close to integer value, use numerical integration
      const double tMin = std::sqrt(1.0 - fThetaEm/actualLowAngle);
      const double tMax = std::sqrt(1.0 - fThetaEm/highAngle);
      utl::Integrator<decltype(integrand)> integrator(integrand, 1e-6);
      integral += (4.0 / fThetaEm) * integrator.GetRombergIntegral(tMin, tMax);
    } else {
      // Use power series
      try {
        integral += rSeries(highAngle, fOmega1) - rSeries(actualLowAngle, fOmega1) +
          fEpsilon * (rSeries(highAngle, fOmega2) - rSeries(actualLowAngle, fOmega2));
      } catch (...) {
        const std::string message = "Failed to integrate the angular distribution of Cherenkov photons";
        ERROR(message);
        throw utl::DoesNotComputeException(message);
      }
    }
  }

        return integral;
}



double
atm::Arbeletche2021CherenkovAngularModel::PDF(
        const double theta, 
        const double showerAge,
        const double refractiveIndex
) const
{
  static const double peakWidth = 1e-5;

        // Check the value of theta
        if (theta <= 0) {
                return 0;
        } else if (theta > fMaxTheta) {
                return 0;
  }

        // Update the parameters that depend on refractive index and shower age
        UpdateParameters(showerAge, refractiveIndex);
        
        if (theta <= fThetaEm - peakWidth) {
                // Region to the left of the peak: compute the distribution function as defined
                return  FunctionK(fThetaEm) * (kPi - std::log(1.0-theta/fThetaEm)) *
                        (std::sin(theta)/fSinThetaEm) / fNormalization;
        }
        else if (theta < fThetaEm + peakWidth) {
                // Very close to the peak: compute the average in an interval containing the peak
                return UnnormalizedIntegral(fThetaEm - peakWidth, fThetaEm + peakWidth) /
                        (2.0*peakWidth*fNormalization);
        }
        else {
                // Right of the peak: compute the distribution function as defined
                return FunctionK(theta) * (kPi - std::log(1.0 - fThetaEm/theta)) / fNormalization;
        }
}



double
atm::Arbeletche2021CherenkovAngularModel::CDF(
        const double theta, 
        const double showerAge,
        const double refractiveIndex
) const
{
        // Check value of theta
        if (theta <= 0) {
                return 0;
        } else if (theta >= fMaxTheta) {
                return 1;
        }
        
        // Update the parameters that depend on refractive index and shower age
        UpdateParameters(showerAge, refractiveIndex);
        
        // Return the (normalized) integral of the distribution between 0 and theta (given)
        return UnnormalizedIntegral(0, theta) / fNormalization;
}



double
atm::Arbeletche2021CherenkovAngularModel::Integral(
        const double lowAngle, 
        const double highAngle, 
        const double showerAge,
        const double refractiveIndex
) const
{
        // Reorder angles, if necessary
        if (lowAngle > highAngle) {
                return -Integral(highAngle, lowAngle, showerAge, refractiveIndex);
        }
        
        // Update the parameters that depend on refractive index and shower age
        UpdateParameters(showerAge, refractiveIndex);
        
        double integral = 0;
        
        if (lowAngle < fThetaEm) {
          integral += UnnormalizedIntegralSineLeft(std::max(lowAngle,0.), std::min(highAngle,fThetaEm));
        }
        
        if (highAngle > fThetaEm) {
          integral += UnnormalizedIntegralSineRight(std::max(lowAngle,fThetaEm), std::min(highAngle,fMaxTheta));
        }
        
        // return the normalized integral
        return integral / fNormalization;
}