AnalyticalCherenkovModel.cc

Go to the documentation of this file.

#include <sstream>
#include <iostream>

#include <utl/TabulatedFunction.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Transformation.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/Point.h>
#include <utl/StringCompare.h>
#include <utl/TimeStamp.h>

#include <atm/ProfileResult.h>
#include <atm/AnalyticalCherenkovModel.h>

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/MagneticFieldModel.h>
#include <fwk/RunController.h>

#include <det/Detector.h>

#include <evt/Event.h>
#include <evt/Header.h>

using namespace det;
using namespace utl;
using namespace atm;
using namespace fwk;
using namespace std;


const double AnalyticalCherenkovModel::fdE = 10*MeV;
const double AnalyticalCherenkovModel::fjMax = 800*MeV;

const double AnalyticalCherenkovModel::fNerlingMinValidShowerAge = 0.1;
const double AnalyticalCherenkovModel::fNerlingMaxValidShowerAge = 2.;

// model parameters used in AngularCDF(), AngularPDF() and CherenkovIntegral()



// --------------- Hillas parameters -----------------
inline double HillasAngle(const double energyThreshold) {

  const double kAngHillas_aPar = 0.83;
  const double kAngHillas_bPar = -0.67;

   // Parametrization with et in MeV
   return kAngHillas_aPar * pow(energyThreshold/MeV, kAngHillas_bPar);

}

// --------------- Nerling parameters --------------
inline double NerlingNormA(const double s) {
  const double kAngNerling_a0 =  4.2489e-1;
  const double kAngNerling_a1 =  5.8371e-1;
  const double kAngNerling_a2 = -8.2373e-2;
  return kAngNerling_a0 + kAngNerling_a1 * s + kAngNerling_a2 * s*s;
}

inline double NerlingNormB(const double s) {
  const double kAngNerling_b0 =  5.5108e-2;
  const double kAngNerling_b1 = -9.5587e-2;
  const double kAngNerling_b2 =  5.6952e-2;
  return kAngNerling_b0 + kAngNerling_b1 * s + kAngNerling_b2 * s*s;
}

inline double NerlingAngle(const double energyThreshold) {
  const double kAngNerling_thc0 =  6.2694e-1;
  const double kAngNerling_thc1 = -6.0590e-1;
  return kAngNerling_thc0 * pow(energyThreshold/MeV, kAngNerling_thc1);
}

inline double NerlingAngleFactor(const double s) {
  const double kAngNerling_thcc0 =  1.0509e+1;
  const double kAngNerling_thcc1 = -4.9644;
  return kAngNerling_thcc0 + kAngNerling_thcc1 * s;
}

// ---------- Giller parameters --------------
enum EGillerAngularParameter {
  eTheta0 = 0,
  eA1,
  eC1,
  eC2,
  eAlpha,
  eNGillerAngular, // number of independent parameters
  eA2              // derived parameter
};

double
GetGillerAngularParameter(const EGillerAngularParameter par,
                          const double showerAge,
                          const double height)
{
  const unsigned int nValues = 4;
  const double parArray[eNGillerAngular][nValues] =
    { {6.058, -1.103e-3, -2.886e-3, -5.447},
      {2.905, -3.851e-2, 1.072e-1, 1.066e+1},
      {7.320, -3.778e-1, 7.202e-2, 7.143},
      {-2.754, -4.242e-1, 1.084e-1, 3.344},
      {2.620, 6.837e-2, -2.247e-2, 1.110} };
  const double additionalC2Par = -4.294;

  const double* p = parArray[par];
  const double s = showerAge;
  const double h = height / km;

  if (par != eC2)
    return p[0] * exp(p[1] * s + p[2] * h) + p[3];
  else {
    const double s2 = s*s;
    return p[0] * exp(p[1] * s2 + p[2] * h) +
      (p[3] * s2 + additionalC2Par * s);
  }
}


// utility function for CherenkovIntegral - int(sin(x)*exp(-x/t)/t
inline double IntegralExpSin(const double x, const double t)  {

  return -exp(-x/t)*(t*cos(x)+sin(x))/(t*t+1);

}

// utility function for CherenkovIntegral - int(sin(x)*exp(-x/t)*exp(-ax)/t 
inline double IntegralExpSinExp(const double x, const double t, const double a)  {

  return -exp(-x*(1./t+a))*(a*t*sin(x)+t*cos(x)+sin(x))/(t*t*(a*a+1)+2*t*a+1);

}

AnalyticalCherenkovModel::AnalyticalCherenkovModel() :
  fArbeletcheAngularDistribution(1e6, kPiOnTwo),
  fShowerAgePrevious(-10)
{
  SetEnergyCutoff(1.*MeV);
}


void
AnalyticalCherenkovModel::Init()
{
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("AnalyticalCherenkovModel");

  topB.GetChild("wavelength").GetData(fWavelength);

  const string param = topB.GetChild("parametrization").Get<string>();
  if (param == "eHillas")
    fParam = eHillas;
  else if (param == "eGiller")
    fParam = eGiller;
  else if (param == "eNerling")
    fParam = eNerling;
  else {
    ERROR("Parametrization not implemented");
    throw utl::NonExistentComponentException();
  }

  const string paramAngular =
    topB.GetChild("parametrization_angular").Get<string>();
  if (paramAngular == "eAngHillas")
    fParamAngular = eAngHillas;
  else if (paramAngular == "eAngNerling")
    fParamAngular = eAngNerling;
  else if (paramAngular == "eAngGiller")
    fParamAngular = eAngGiller;
  else if (paramAngular == "eAngArbeletche")
        fParamAngular = eAngArbeletche;
  else {
    ERROR("Parametrization Angular not implemented");
    throw utl::NonExistentComponentException();
  }

  const string cherenkovDistribution =
    topB.GetChild("cherenkovDistribution").Get<string>();
  if (cherenkovDistribution == "symmetric")
    fCorrectCherenkovAsymmetry = false;
  else if (cherenkovDistribution == "asymmetric")
    fCorrectCherenkovAsymmetry = true;
  else {
    ERROR("Unknown \"cherenkovDistribution\"");
    throw utl::NonExistentComponentException();
  }

  Branch asymmLimits = topB.GetChild("asymmetryCorrectionLimits");
  asymmLimits.GetChild("minAge").GetData(fCherenkovAsymmetryMinAge);
  asymmLimits.GetChild("maxAge").GetData(fCherenkovAsymmetryMaxAge);
  asymmLimits.GetChild("minAngle").GetData(fCherenkovAsymmetryMinAngle);
  asymmLimits.GetChild("maxAngle").GetData(fCherenkovAsymmetryMaxAngle);
  asymmLimits.GetChild("minA").GetData(fCherenkovAsymmetryMinA); 
  asymmLimits.GetChild("maxA").GetData(fCherenkovAsymmetryMaxA); 

  topB.GetChild("wavelengthDepRefraction").GetData(fWlRefrac);
  topB.GetChild("lowEnergyCorrectedNerling").GetData(fLowENerling);

}


utl::TabulatedFunction&
AnalyticalCherenkovModel::EvaluateCherenkovDirect(const utl::Point& xA,
                                                  const utl::Point& xB,
                                                  const utl::Point& xEye,
                                                  const double showerAge)
  const
{
  const double prob = fWlRefrac ? 1. : EvaluateDirectCherenkovProbability (xA, xB, xEye, showerAge);

  const TabulatedFunction& photons = EvaluateCherenkovPhotons(xA, xB, showerAge);

  const int numWaves = fWavelength.size();
  vector<double> directCher;
  for (int i = 0; i < numWaves; ++i) {
    const double probWl = fWlRefrac ? EvaluateDirectCherenkovProbability (xA, xB, xEye, showerAge, fWavelength[i]) : 1.;
    directCher.push_back(photons[i].Y()*prob*probWl);
  }

  fDirectCherenkovPhotons = utl::TabulatedFunction(fWavelength, directCher);
  return fDirectCherenkovPhotons;
}


double
AnalyticalCherenkovModel::EvaluateDirectCherenkovProbability(const utl::Point& xA,
                                                             const utl::Point& xB,
                                                             const utl::Point& xEye,
                                                             const double showerAge)
  const
{
  Detector& theDetector = Detector::GetInstance();
  const CoordinateSystemPtr cs = theDetector.GetSiteCoordinateSystem();
  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();

  const Point xMean((xA.GetX(cs) + xB.GetX(cs))*0.5,
                    (xA.GetY(cs) + xB.GetY(cs))*0.5,
                    (xA.GetZ(cs) + xB.GetZ(cs))*0.5, cs);
  const Vector xAxEyeDir = xEye - xA;
  const Vector xBxEyeDir = xEye - xB;

  double aBin = Angle(xAxEyeDir, xBxEyeDir);

  // Calculate impact parameter rp
//#warning this is not generally true! Better change the interface to the Cherenkov Model.
  const Vector shwDir = Normalized(xB - xA);
  const double projection = shwDir*xAxEyeDir;
  const double rp = (xA + projection*shwDir - xEye).GetMag();
  const Vector xMeanVec = xMean - xEye;

  // Calculate emission angle.
  const double emissAngle = Angle(-shwDir,xMeanVec);

  const double hA = (wgs84.PointToLatitudeLongitudeHeight(xA)).get<2>();
  const double hB = (wgs84.PointToLatitudeLongitudeHeight(xB)).get<2>();

  double lowAngle = emissAngle - aBin*0.5;
  if (lowAngle <= 0.)
    lowAngle = emissAngle*0.5;
  const double highAngle = emissAngle + aBin*0.5;
  aBin = highAngle-lowAngle;

  const ProfileResult& dvh = theDetector.GetAtmosphere().EvaluateDepthVsHeight();
  const double atmHeightMin = dvh.MinX();
  const double atmHeightMax = dvh.MaxX();
  double depthA = 0.;
  double depthB = 0.;
  if (hA < atmHeightMin)
    depthA = dvh.Y(atmHeightMin);
  else if (hA > atmHeightMax)
    depthA = dvh.Y(atmHeightMax);
  else
    depthA = dvh.Y(hA);

  if (hB < atmHeightMin)
    depthB = dvh.Y(atmHeightMin);
  else if (hB > atmHeightMax)
    depthB = dvh.Y(atmHeightMax);
  else
    depthB = dvh.Y(hB);

  const double meanDepth = (depthA + depthB)*0.5;

  // Refraction Index for Middle of the Bin
  const double nIndex = RefractionIndex::LorentzLorentz(meanDepth);

  double eTimesSine = CherenkovIntegral(lowAngle, highAngle, showerAge, nIndex);

  if (CorrectCherenkovAsymmetry(showerAge, emissAngle))
    eTimesSine *= AsymmCorrection(xA, xB, xEye);

  // Area of tilted annulus on ground that receives
  // light = dA, is area along line of sight = dA = 2*pi*Rp*dx
  // where dx = distance to the track * abin = Rp/sin(emissAngle)*abin
  // The sin(emissAngle) term is contained inside eTimesSine.

  const double dAModified = kTwoPi*rp*rp*aBin;

  return eTimesSine/dAModified;
}

double
AnalyticalCherenkovModel::EvaluateDirectCherenkovProbability(const utl::Point& xA,
                                                             const utl::Point& xB,
                                                             const utl::Point& xEye,
                                                             const double showerAge,
                                                             const double wavelength)
  const
{
  Detector& theDetector = Detector::GetInstance();
  const CoordinateSystemPtr cs = theDetector.GetSiteCoordinateSystem();
  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();

  const Point xMean((xA.GetX(cs) + xB.GetX(cs))*0.5,
                    (xA.GetY(cs) + xB.GetY(cs))*0.5,
                    (xA.GetZ(cs) + xB.GetZ(cs))*0.5, cs);
  const Vector xAxEyeDir = xEye - xA;
  const Vector xBxEyeDir = xEye - xB;

  double aBin = Angle(xAxEyeDir, xBxEyeDir);

  // Calculate impact parameter rp
//#warning this is not generally true! Better change the interface to the Cherenkov Model.
  const Vector shwDir = Normalized(xB - xA);
  const double projection = shwDir*xAxEyeDir;
  const double rp = (xA + projection*shwDir - xEye).GetMag();
  const Vector xMeanVec = xMean - xEye;

  // Calculate emission angle.
  const double emissAngle = Angle(-shwDir,xMeanVec);

  const double hA = (wgs84.PointToLatitudeLongitudeHeight(xA)).get<2>();
  const double hB = (wgs84.PointToLatitudeLongitudeHeight(xB)).get<2>();

  double lowAngle = emissAngle - aBin*0.5;
  if (lowAngle <= 0.)
    lowAngle = emissAngle*0.5;
  const double highAngle = emissAngle + aBin*0.5;
  aBin = highAngle-lowAngle;

  const ProfileResult& dvh = theDetector.GetAtmosphere().EvaluateDepthVsHeight();
  const double atmHeightMin = dvh.MinX();
  const double atmHeightMax = dvh.MaxX();
  /*
  double depthA = 0.;
  double depthB = 0.;
  if (hA < atmHeightMin)
    depthA = dvh.Y(atmHeightMin);
  else if (hA > atmHeightMax)
    depthA = dvh.Y(atmHeightMax);
  else
    depthA = dvh.Y(hA);
  
  if (hB < atmHeightMin)
    depthB = dvh.Y(atmHeightMin);
  else if (hB > atmHeightMax)
    depthB = dvh.Y(atmHeightMax);
  else
    depthB = dvh.Y(hB);

  const double meanDepth = (depthA + depthB)*0.5;
  */

  const double meanHeight = ((hA+hB)/2.) < atmHeightMin ? atmHeightMin : ( ((hA+hB)/2.) > atmHeightMax ? atmHeightMax : ((hA+hB)/2.) );
  const double pres = theDetector.GetAtmosphere().EvaluatePressureVsHeight().Y(meanHeight);
  const double vapPres = theDetector.GetAtmosphere().EvaluateVaporPressureVsHeight().Y(meanHeight);
  const double temp = theDetector.GetAtmosphere().EvaluateTemperatureVsHeight().Y(meanHeight);
  // Refraction Index for Middle of the Bin
  const double nIndex = RefractionIndex::Ciddor95(wavelength, temp, pres, vapPres);

  double eTimesSine = CherenkovIntegral(lowAngle, highAngle, showerAge, nIndex);

  if (CorrectCherenkovAsymmetry(showerAge, emissAngle))
    eTimesSine *= AsymmCorrection(xA, xB, xEye);

  // Area of tilted annulus on ground that receives
  // light = dA, is area along line of sight = dA = 2*pi*Rp*dx
  // where dx = distance to the track * abin = Rp/sin(emissAngle)*abin
  // The sin(emissAngle) term is contained inside eTimesSine.

  const double dAModified = kTwoPi*rp*rp*aBin;

  return eTimesSine/dAModified;
}

bool
AnalyticalCherenkovModel::CorrectCherenkovAsymmetry(double age, double theta) const
{
  return fCorrectCherenkovAsymmetry
    && (age > fCherenkovAsymmetryMinAge)
    && (age < fCherenkovAsymmetryMaxAge)
    && (theta > fCherenkovAsymmetryMinAngle)
    && (theta < fCherenkovAsymmetryMaxAngle); 
  // a accounted for separately in the correction routine, as it requires geometry calculation
}

utl::TabulatedFunction&
AnalyticalCherenkovModel::EvaluateCherenkovPhotons(const utl::Point& xA,
                                                   const utl::Point& xB,
                                                   const double showerAge)
  const
{

  if (xA == fxAprevious && xB == fxBprevious && showerAge == fShowerAgePrevious)
    return fCherenkovPhotons;

  fxAprevious = xA;
  fxBprevious = xB;
  fShowerAgePrevious = showerAge;

  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();

  const double distAB = (xA - xB).GetMag();

  const double hA = (wgs84.PointToLatitudeLongitudeHeight(xA)).get<2>();
  const double hB = (wgs84.PointToLatitudeLongitudeHeight(xB)).get<2>();

  Detector& theDetector = Detector::GetInstance();
  const ProfileResult& dvh = theDetector.GetAtmosphere().EvaluateDepthVsHeight();
  const double atmHeightMin = dvh.MinX();
  const double atmHeightMax = dvh.MaxX();
  double depthA = 0.;
  double depthB = 0.;
  if (hA < atmHeightMin)
    depthA = dvh.Y(atmHeightMin);
  else if (hA > atmHeightMax)
    depthA = dvh.Y(atmHeightMax);
  else
    depthA = dvh.Y(hA);

  if (hB < atmHeightMin)
    depthB = dvh.Y(atmHeightMin);
  else if (hB > atmHeightMax)
    depthB = dvh.Y(atmHeightMax);
  else
    depthB = dvh.Y(hB);

  const double meanDepth = (depthA + depthB)*0.5;

  const double meanHeight = ((hA+hB)/2.) < atmHeightMin ? atmHeightMin : ( ((hA+hB)/2.) > atmHeightMax ? atmHeightMax : ((hA+hB)/2.) );
  const double pres = theDetector.GetAtmosphere().EvaluatePressureVsHeight().Y(meanHeight);
  const double vapPres = theDetector.GetAtmosphere().EvaluateVaporPressureVsHeight().Y(meanHeight);
  const double temp = theDetector.GetAtmosphere().EvaluateTemperatureVsHeight().Y(meanHeight);
  // Refraction Index for Middle of the Bin
  // Assumes that fWavelength[0] is the lowest wavelength (means the lowest thresholdEnergy)
  const double nIndex = fWlRefrac ? RefractionIndex::Ciddor95(fWavelength[0], temp, pres, vapPres) : RefractionIndex::LorentzLorentz(meanDepth);

  // Emission Threshold Energy (assumed electron mass)
  const double thresholdEnergy = fmin(utl::CherenkovThreshold(nIndex), 2.e3*MeV);
  
  // From Threshold to Ultra-relativistic
  const int numWaves = fWavelength.size();

  vector<double> nPhotons(numWaves,0);

  for (double energy = thresholdEnergy; energy < fjMax; energy += fdE) {

    double deltaLength = 0.;
    if (fParam == eHillas) {
      const double tl1 = CherenkovTrackLength(energy, showerAge);
      const double tl2 = CherenkovTrackLength(energy+fdE, showerAge);
      deltaLength = tl1 - tl2;
    } else if (fParam == eGiller || fParam == eNerling)
      deltaLength = CherenkovTrackLength(energy+fdE*0.5, showerAge);
    else {
      ERROR("Non-existent parametrization");
      throw utl::NonExistentComponentException();
    }

    const double eelec = energy + fdE*0.5;

    // Note that wavelength bins do not have to be constant
    // RU: THIS IS A FAILURE OF THE CHERNKOV MODEL DESIGN !!!!!
    double minWaveBin;
    double maxWaveBin = 0;
    for (int i = 0; i < numWaves; ++i) {

      const double nIndexWl = RefractionIndex::Ciddor95(fWavelength[i], temp, pres, vapPres);
      const double thresholdEnergyWl = fmin(utl::CherenkovThreshold(nIndexWl), 2.e3*MeV);
      if (fWlRefrac && (energy < thresholdEnergyWl)) continue;

      if (!i)
        minWaveBin = fWavelength[i] - 0.5*(fWavelength[i+1] - fWavelength[i]);
      else
        minWaveBin = maxWaveBin;

      if (i == numWaves-1)
        maxWaveBin = fWavelength[i] + 0.5*(fWavelength[i]-fWavelength[i-1]);
      else
        maxWaveBin = 0.5*(fWavelength[i] + fWavelength[i+1]);

      nPhotons[i] += deltaLength * DeltaPhotons(fWlRefrac ? nIndexWl : nIndex, eelec, distAB,
                                                minWaveBin, maxWaveBin);
    }
  } // energy loop


  const double ehigh = 10.*GeV;

  // Adding ultra-relativistic energy contribution
  double trackLength = 0.;
  if (fParam == eHillas)
    trackLength = CherenkovTrackLength(fjMax, showerAge);
  //F.N. BUG 329
  else if (fParam == eGiller || fParam == eNerling) {
    trackLength = 0;
    for (double jj = fjMax; jj < ehigh; jj += fdE) {
      trackLength += CherenkovTrackLength(jj + fdE*0.5, showerAge);
    }
  } else {
    ERROR("Non-existent parametrization");
    throw utl::NonExistentComponentException();
  }

  double minWaveBin;
  double maxWaveBin = 0;
  for (int i = 0; i < numWaves; ++i) {

    if (!i)
      minWaveBin = fWavelength[i] - 0.5*(fWavelength[i+1] - fWavelength[i]);
    else
      minWaveBin = maxWaveBin;

    if (i == numWaves-1)
      maxWaveBin = fWavelength[i] + 0.5*(fWavelength[i] - fWavelength[i-1]);
    else
      maxWaveBin = 0.5*(fWavelength[i] + fWavelength[i+1]);

    const double nIndexWl = RefractionIndex::Ciddor95(fWavelength[i], temp, pres, vapPres);

    nPhotons[i] += trackLength * DeltaPhotons(fWlRefrac ? nIndexWl : nIndex, ehigh, distAB,
                                              minWaveBin, maxWaveBin);
  }

  fCherenkovPhotons = utl::TabulatedFunction(fWavelength, nPhotons);

  return fCherenkovPhotons;
}

double
AnalyticalCherenkovModel::CherenkovTrackLength(const double energy,
                                               const double showerAge)
  const
{

  if (fParam == eHillas) {
    /* A.M. Hillas parametrization.

       "Angular and energy distributions of charged
        particles in electron-photon cascades in air"

       J.Phys.G: Nucl.Phys. 8 --> 1982 <-- (pag: 1461-1473)
    */

    const double e0 = (showerAge < 0.4) ?
      26.0*MeV :
      (44.0 - 17.0*(showerAge - 1.46)*(showerAge - 1.46))*MeV;

    // The following formula is valid for energy in units of MeV
    const double trackLength =
      pow((0.89*e0/MeV - 1.2)/(e0/MeV + energy/MeV), showerAge) /
      pow((1.0 + 1.0e-4*showerAge*energy/MeV), 2);

    return trackLength;

  } else if (fParam == eGiller) {
    /* Maria Giller parametrization.

       "Energy spectra of electrons in the extensive
        air showers of ultra-high energy"

       J.Phys.G: Nucl.Part.Phys 30 --> 2004 <-- (pag:97-105)
    */

    const double a = 1.005;
    const double b = 0.06;
    const double c = 189.0;

    const double cs = 0.111*showerAge + 0.134;
    const double ds = 7.06*showerAge + 12.48;

    const double energyCr = 80.0*MeV;

    const double auxExponent = -(showerAge + b*log(energy/(c*energyCr)));

    const double dist_giller = cs*(1.0 - a*exp(-ds*energy/energyCr)) *
      pow(1.0 + energy/energyCr, auxExponent);

    const double trackLength = dist_giller *(fdE/energy);

    return trackLength;

  } else if (fParam == eNerling) {
    const double truncatedAge = fmax(fmin(fNerlingMaxValidShowerAge, showerAge),
                                     fNerlingMinValidShowerAge);
    const double a1 = 6.42522 - 1.53183 * truncatedAge;
    const double a2 = 168.168 - 42.1368 * truncatedAge;


    const double a0 = fNerlingFactor[0] *
      exp(fNerlingFactor[1]*truncatedAge + fNerlingFactor[2]*truncatedAge*truncatedAge);

    // The following formula is valid for energy in units of MeV !!!
    double dist_nerling = a0 / (energy/MeV + a1) / pow(energy/MeV + a2, truncatedAge);

    const double trackLength = dist_nerling * (fdE/MeV);

    return trackLength;

  } else {
    ERROR("Non-existent parametrization");
    throw utl::NonExistentComponentException();
  }
}


double
AnalyticalCherenkovModel::DeltaPhotons(const double nIndex,
                                       const double eelec, const double dSlant,
                                       const double wl1, const double wl2)
  const
{
  const double xme = kElectronMass;
  const double delta = nIndex - 1.;
  const double alpha = 1./137.035989;
  const double photons = (kTwoPi*alpha) * abs(1.0/wl1 - 1.0/wl2) *
    (2.*delta - xme*xme/(eelec*eelec))*dSlant;
  //F.N.    *(1.0 - 1.0/((beta*nIndex)*(beta*nIndex)))*dSlant;
  return photons;
}


double
AnalyticalCherenkovModel::CherenkovIntegral(const double lowAngle,
                                            const double highAngle,
                                            const double showerAge,
                                            const double refractiveIndex)
  const
{
        // Threshold energy for the emission of Cherenkov light
        const double energyThreshold = fmin(utl::CherenkovThreshold(refractiveIndex),  2.e3*MeV);

  if (fParamAngular == eAngHillas) {

    const double t = HillasAngle(energyThreshold);
    return IntegralExpSin(highAngle,t)-IntegralExpSin(lowAngle,t);

  }
  else if (fParamAngular == eAngNerling) {

    const double truncatedAge = fmax(fmin(fNerlingMaxValidShowerAge, showerAge),
                                     fNerlingMinValidShowerAge);
    const double A = NerlingNormA(truncatedAge);
    const double B = NerlingNormB(truncatedAge);
    const double c = NerlingAngleFactor(truncatedAge);
    const double t1  = NerlingAngle(energyThreshold);
    const double t2 = c * t1;

    //standard Nerling parametrization
    if (!fLowENerling || highAngle > 30*deg) {
      return A*(IntegralExpSin(highAngle,t1)-IntegralExpSin(lowAngle,t1))+
           B*(IntegralExpSin(highAngle,t2)-IntegralExpSin(lowAngle,t2));
    }

    //correction to Nerling model based on 3D CORSIKA simulations at low energy (10^15-10^17 eV)
    double par[5] = {6.98169e-01, -2.86420e-01, 1.76411e-01/deg, -2.36286e+00, -7.70755e-02/deg};
    return par[0]*( A*(IntegralExpSin(highAngle,t1)-IntegralExpSin(lowAngle,t1))+
           B*(IntegralExpSin(highAngle,t2)-IntegralExpSin(lowAngle,t2)) ) +
           exp(par[1])*( A*(IntegralExpSinExp(highAngle,t1,par[2])-IntegralExpSinExp(lowAngle,t1,par[2]))+
           B*(IntegralExpSinExp(highAngle,t2,par[2])-IntegralExpSinExp(lowAngle,t2,par[2])) ) +
           exp(par[3])*( A*(IntegralExpSinExp(highAngle,t1,par[4])-IntegralExpSinExp(lowAngle,t1,par[4]))+
           B*(IntegralExpSinExp(highAngle,t2,par[4])-IntegralExpSinExp(lowAngle,t2,par[4])) );

  }
  else if (fParamAngular == eAngArbeletche) {
    return fArbeletcheAngularDistribution.Integral(lowAngle, highAngle, showerAge, refractiveIndex);
  }
  else {
    ERROR("Non-existent parametrization");
    throw utl::NonExistentComponentException();
    return 0.;
  }

}

double
AnalyticalCherenkovModel::AngularCDF(const double theta,
                                     const double verticalDepth,
                                     const double showerAge)
  const
{

  if (theta < 0.)
    return 0.;
  else if (theta > kPiOnTwo)
    return 1.;

  // emission energy threshold
  const double refractiveIndex = RefractionIndex::LorentzLorentz(verticalDepth);
  const double energyThreshold =
    fmin(utl::CherenkovThreshold(refractiveIndex), 2.e3*MeV);

  if (fParamAngular == eAngHillas) {

    const double elecAngle = HillasAngle(energyThreshold);
    const double hillasCDF = 1.- exp(-theta/elecAngle);

    return hillasCDF;

  }
  else if (fParamAngular == eAngNerling) {

    /* Parametrisation of Cherenkov Light Production by Air Showers
       F. Nerling et al., Astropart.Phys.24:421-437,2006
    */

    // shower age
    const double s = fmax(fmin(fNerlingMaxValidShowerAge, showerAge),
                          fNerlingMinValidShowerAge);
    const double A = NerlingNormA(s);
    const double B = NerlingNormB(s);
    const double c = NerlingAngleFactor(s);
    const double th_c  = NerlingAngle(energyThreshold);
    const double th_cc = c * th_c;

    const double nerlingCDF = A*(1 - exp(-theta/th_c)) +
                              B*(1 - exp(-theta/th_cc));

    // normalize to 0-90 deg.
    const double nerlingNorm = A*(1 - exp(-kPiOnTwo/th_c)) +
                               B*(1 - exp(-kPiOnTwo/th_cc));

    return nerlingCDF/nerlingNorm;

  }
  else if (fParamAngular == eAngArbeletche) {
    return fArbeletcheAngularDistribution.CDF(theta, showerAge, refractiveIndex);
  }
  else {
    ERROR("Non-existent CDF parametrization");
    throw utl::NonExistentComponentException();
  }
}

/* Helper functions for asymmetric Cherenkov correction */
enum charged_particles {electrons, positrons};

inline double 
b1param (double theta, int species)
{
  static const double b11[]=
    {/*[electrons]=*/ -8.2e-5/(degree*degree), /*[positrons]=*/ -6.2e-5/(degree*degree)};
  static const double b12[]=
    {/*[electrons]=*/ 2.8e-3/degree, /*[positrons]=*/ 2.2e-3/degree};
  static const double b13[]=
    {/*[electrons]=*/ 1.24e-2, /*[positrons]=*/ 1.36e-2};
  
  return b11[species] * theta * theta  +  b12[species] * theta  +  b13[species];
}

inline double
bparam (double a, double theta, int species)
{
  static const double b2 [] = {/*[electrons]=*/ -8.4e-3, /*[positrons]=*/ -9.3e-3};
  return b1param(theta, species) * a + b2[species];
}

double
AnalyticalCherenkovModel::AsymmCorrection(double a, double theta, double sinphi) const
{
  const double bsum = bparam (a, theta, electrons) + bparam (a, theta, positrons);
  const double bdiff = bparam (a, theta, electrons) - bparam (a, theta, positrons);
  const double Fs = 0.5 * bsum;
  const double Fa = bdiff / bsum;

  const double binning = 10./360.;
  const double A = binning - 3./8. * Fs;

  const double sinphi4=(sinphi*sinphi) * (sinphi*sinphi);

  return (1./binning) * (A  +  Fs * (1 + Fa*sinphi) * sinphi4);
}

inline utl::Vector
perpendicular (const utl::Vector& v, const utl::Vector& axis)
{
  return v - axis*v*axis;
}

inline double
vectorsine (const utl::Vector& u, const utl::Vector& v) 
{
  return Cross(u,v).GetMag() / (u.GetMag() * v.GetMag());
}

inline double
vectorsine_oriented (const utl::Vector& u, const utl::Vector& v, const utl::Vector& up) 
{
  const utl::Vector uXv = Cross(u,v);
  const double sin = uXv.GetMag() / (u.GetMag() * v.GetMag());
  const int neg = signbit (uXv * up); //negative if uxv not "upward"
  return neg?-sin:sin;
}

//inline double
//vectorsine_inplane (const utl::Vector& u, const utl::Vector& v, const utl::Vector& up)
//{
//  const Vector u_perp = perpendicular (u, up);
//  const Vector v_perp = perpendicular (v, up);
//  return vectorsine_oriented (u_perp, v_perp, up);
//}

inline utl::Vector
bFieldAt (utl::Point x)
{
  Detector& det = Detector::GetInstance();
  return det.GetGeoMagneticFieldAt (x);
}

double
AnalyticalCherenkovModel::AsymmCorrection (const utl::Point& emissionPoint,
                                           const utl::Vector& v_toeye, 
                                           const utl::Vector& axis,
                                           const double density) const
{
  Vector b = bFieldAt (emissionPoint);

  Vector v_perp = perpendicular (v_toeye, axis);
  Vector b_perp = perpendicular (b, axis);

  double sinphi = vectorsine_oriented (b_perp, v_perp, -axis);
  double theta = Angle(axis, v_toeye);

  double a = b_perp.GetMag()*(1./gauss) / density * (kg/m3);
  if (a < fCherenkovAsymmetryMinA)
    return 1.0;
  if (a > fCherenkovAsymmetryMaxA)
    a = fCherenkovAsymmetryMaxA;
  
  return AsymmCorrection (a, theta, sinphi);
}

double
AnalyticalCherenkovModel::AsymmCorrection (const utl::Point& xA,
                                           const utl::Point& xB,
                                           const utl::Point& xEye) const
{
  utl::Vector axis=Normalized(xB - xA);
  
  const CoordinateSystemPtr cs = det::Detector::GetInstance().GetSiteCoordinateSystem();
  const Point xMean((xA.GetX(cs) + xB.GetX(cs))*0.5,
                    (xA.GetY(cs) + xB.GetY(cs))*0.5,
                    (xA.GetZ(cs) + xB.GetZ(cs))*0.5, cs);
  Vector xMeanVec = xMean-xEye;
  double altitude=xMean.GetZ(cs);

  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  double density = atmo.EvaluateDensityVsHeight().Y(altitude);

  return AsymmCorrection(xMean, -xMeanVec, axis, density);
}


// Old Cherenkov asymmetry parametrisation:

double
fparam (double a, double theta, int species)
{
  static const double b11[]={/*[electrons]=*/ 0.00206, /*[positrons]=*/0.00156};
  static const double b12[]={/*[electrons]=*/ 0.0072, /*[positrons]=*/ 0.0065};
  
  const double b1=b11[species] * theta/degree + b12[species];
  return b1 * a*a;
}

double
AnalyticalCherenkovModel::AsymmCorrectionOld (const double bm,
                                           const double bp,
                                           const double sinphi) const
{
  const double binning=10./360.;
  const double sinphi4= (sinphi*sinphi)*(sinphi*sinphi);
  return (1./binning)  *  ( (binning - bp*3./8.)  +  bp * (1+bm*sinphi) * sinphi4 );
}

double
AnalyticalCherenkovModel::AsymmCorrectionOld (const utl::Vector& v_toeye, 
                                           const utl::Vector& axis,
                                           const double density) const
{
  CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
  Vector b(24.6*microtesla, 55*degree, 87*degree, siteCS, Vector::kSpherical);

  Vector v_perp = perpendicular (v_toeye, axis);
  Vector b_perp = perpendicular (b, axis);

  double sinphi = vectorsine_oriented (b_perp, v_perp, -axis);
  double theta = Angle(axis, v_toeye);

  double a = b_perp.GetMag()*(1./gauss) / density * (kg/m3);
  double fm = fparam (a, theta, electrons);
  double fp = fparam (a, theta, positrons);
  double bm = (fm - fp) / (fm + fp);
  double bp = 0.5* (fm + fp);

  return AsymmCorrectionOld (bm, bp, sinphi);
}


double
AnalyticalCherenkovModel::AngularPDF(const double theta,
                                     const double verticalDepth,
                                     const double showerAge)
  const
{

  if (theta < 0.)
    return 0.;
  else if (theta > kPiOnTwo)
    return 0.;

  // emission energy threshold
  const double refractiveIndex = RefractionIndex::LorentzLorentz(verticalDepth);
  const double energyThreshold =
    fmin(utl::CherenkovThreshold(refractiveIndex), 2.e3*MeV);

  if (fParamAngular == eAngHillas) {

    const double elecAngle = HillasAngle(energyThreshold);
    // normalize to 0-90 deg.
    const double hillasNorm = 1.- exp(-kPiOnTwo/elecAngle);
    return exp(-theta/elecAngle)/hillasNorm;

  }
  else if (fParamAngular == eAngNerling) {

    // shower age
    const double s = fmax(fmin(fNerlingMaxValidShowerAge, showerAge),
                          fNerlingMinValidShowerAge);
    const double A = NerlingNormA(s);
    const double B = NerlingNormB(s);
    const double c = NerlingAngleFactor(s);
    const double th_c  = NerlingAngle(energyThreshold);
    const double th_cc = c * th_c;

    // normalize to 0-90 deg.
    const double nerlingNorm = A*(1 - exp(-kPiOnTwo/th_c)) +
                               B*(1 - exp(-kPiOnTwo/th_cc));

    return (A*exp(-theta/th_c)/th_c+B*exp(-theta/th_cc)/th_cc)/nerlingNorm;

  }
  else if (fParamAngular == eAngGiller) {

    const ProfileResult& heightVsDepth =
      Detector::GetInstance().GetAtmosphere().EvaluateHeightVsDepth();
    const double height =
      heightVsDepth.Y(fmin(fmax(verticalDepth, heightVsDepth.MinX()),
                           heightVsDepth.MaxX()));

    const double theta0 = GetGillerAngularParameter(eTheta0, showerAge, height);
    if (theta < theta0) {
      const double a1 = GetGillerAngularParameter(eA1, showerAge, height);
      const double c1 = GetGillerAngularParameter(eC1, showerAge, height);
      const double c2 = GetGillerAngularParameter(eC2, showerAge, height);
      return a1 * exp(-(c1 * theta + c2 * theta*theta));
    }
    else {
      const double alpha = GetGillerAngularParameter(eAlpha, showerAge, height);
      const double a2 = GetGillerAngularParameter(eA2, showerAge, height);
      return a2 * pow(theta, -alpha);
    }
  }
  else if (fParamAngular == eAngArbeletche) {
    return fArbeletcheAngularDistribution.PDF(theta, showerAge, refractiveIndex);
  }
  else {
    ERROR("Non-existent PDF parametrization");
    throw utl::NonExistentComponentException();
  }
}


void
AnalyticalCherenkovModel::SetEnergyCutoff(const double ecut)
  const
{
  if (!fNerlingFactor.empty() && fEcut == ecut)
    return;

  fEcut = ecut;

  /* F. Nerling et al., Astropart.Phys.24:421-437,2006 */

  const int nEnergies = 6;

  const double Ecut[nEnergies] = { 2.*MeV, 1.*MeV, 0.5*MeV,
                                   0.25*MeV, 0.1*MeV, 0.05*MeV };

  const double p0[nEnergies] = { 1.48071e-01, 1.45098e-01,
                                 1.43458e-01, 1.42589e-01,
                                 1.42049e-01, 1.41866e-01 };

  const double p1[nEnergies] = { 6.22334, 6.20114,
                                 6.18979, 6.18413,
                                 6.18075, 6.17963 };

  const double p2[nEnergies] = { -5.89710e-01, -5.96851e-01,
                                 -6.01298e-01, -6.03838e-01,
                                 -6.05484e-01, -6.06055e-01 };

  int i1, i2;
  if (fEcut > Ecut[0]) {
    i1 = 0;
    i2 = 1;
  } else {
    i1 = nEnergies - 2;
    i2 = nEnergies - 1;

    for (int i = 1; i < nEnergies; ++i)
      if (fEcut > Ecut[i]) {
        i1 = i - 1;
        i2 = i;
        break;
      }
  }

  fNerlingFactor.resize(3);
  fNerlingFactor[0] =
    p0[i1] + (p0[i2] - p0[i1])/(Ecut[i2] - Ecut[i1]) * (fEcut - Ecut[i1]);
  fNerlingFactor[1] =
    p1[i1] + (p1[i2] - p1[i1])/(Ecut[i2] - Ecut[i1]) * (fEcut - Ecut[i1]);
  fNerlingFactor[2] =
    p2[i1] + (p2[i2] - p2[i1])/(Ecut[i2] - Ecut[i1]) * (fEcut - Ecut[i1]);
}

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. -k"
// End: