CherenkovFluorescenceMatrix.cc

Go to the documentation of this file.

#include "CherenkovFluorescenceMatrix.h"
#include "lowerTriangularMatrix.h"
#include "diagonalMatrix.h"

#include <iomanip>
#include <limits>

#include <boost/io/ios_state.hpp>

#include <utl/Vector.h>
#include <utl/UTMPoint.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/TabulatedFunctionErrors.h>

#include <det/Detector.h>

#include <atm/ProfileResult.h>
#include <atm/Atmosphere.h>
#include <atm/ScatteringResult.h>
#include <atm/AttenuationResult.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/RandomEngineRegistry.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <CLHEP/Random/RandFlat.h>
#include <utl/RandomEngine.h>

using CLHEP::RandFlat;

using namespace fwk;
using namespace det;

using std::cerr;
using std::cout;
using std::endl;

using namespace oBLAS;
using namespace atm;
using namespace det;
using namespace FdProfileReconstructorKG;


#undef _CFM_DEBUG


// constructor with geometry, relative detector efficiency and shower maximum as input

CherenkovFluorescenceMatrix::CherenkovFluorescenceMatrix(
  const double xMax, const utl::Point& eyePosition,
  const std::deque<std::vector<utl::Point>>& showerGeometry,
  const utl::TabulatedFunction& relEff,
  const std::deque<double>& depth,
  const double eCut, const double zeta, const double cFac, const double fFac,
  const eFluorescenceLDF fLDF,
  const eDirectCherenkovLDF dcLDF,
  const eScatteredCherenkovLDF scLDF,
  const eMultipleScatteringLDF msLDF,
  const bool cherenkovCone,
  const OpticalHalo::EHaloType spotHalo
) :
  fXmax(xMax),
  fEcut(eCut),
  fEyePosition(eyePosition),
  fShowerGeometry(showerGeometry),
  fRelEff(relEff),
  fDepth(depth),
  fZeta(zeta),
  fCfac(cFac),
  fFfac(fFac),
  fFluoLDF(fLDF),
  fDirCherLDF(dcLDF),
  fScatCherLDF(scLDF),
  fMultScatLDF(msLDF),
  fCherenkovCone(cherenkovCone),
  fOpticalHalo(spotHalo)
{
#ifdef _CFM_DEBUG
  cout << "        -->CherenkovFluorescenceMatrix()" << endl;
#endif

  // set number of depth points
  const unsigned int nDepth = fShowerGeometry.size();
  for (unsigned int i = 0; i < nDepth; ++i)
    fShowerPoints.push_back(fShowerGeometry[i][0] -
                            0.5*(fShowerGeometry[i][0] - fShowerGeometry[i][1]));

  // create matrices
  fMieScatteredCherenkovMatrix = new lowerTriangularMatrix(nDepth);
  fRayleighScatteredCherenkovMatrix = new lowerTriangularMatrix(nDepth);
  fDirectCherenkovMatrix = new diagonalMatrix(nDepth);
  fFluorescenceMatrix = new diagonalMatrix(nDepth);
  fCherenkovFluorescenceMatrix = new lowerTriangularMatrix(nDepth);

  fTanZeta = tan(fZeta);

  // calculate necessary ingredients ...

  // a) get wavelengths to be used for Cherenkov and Fluorescence
  CalculateWaveLengths();

  // b) geometry dependend stuff
  CalculateAttenuationToEye();
  CalculateAttenuationAlongTrack();
  CalculateGeometricalFactor();
  CalculateScatteringToEye();

  // c) age dependend stuff
  CalculateMeanEnergyDeposit();
  CalculateCherenkovAtTrack();

  // ... and finally the matrices
  CalculateFluorescenceMatrix();
  CalculateDirectCherenkovMatrix();
  CalculateMieAndRayScattCherenkovMatrix();
}


// destructor deletes allocated space

CherenkovFluorescenceMatrix::~CherenkovFluorescenceMatrix()
{
#ifdef _CFM_DEBUG
  cout << "        -->~CherenkovFluorescenceMatrix()" << endl;
#endif
  delete fMieScatteredCherenkovMatrix;
  delete fRayleighScatteredCherenkovMatrix;
  delete fDirectCherenkovMatrix;
  delete fFluorescenceMatrix;
  delete fCherenkovFluorescenceMatrix;
}


void
CherenkovFluorescenceMatrix::CalculateFluorescenceMatrix()
{
#ifdef _CFM_DEBUG
  cout << "        -->CalculateFluorescenceMatrix()" << endl;
#endif

  fFluorescenceMultipleScattering.clear();

  // current atmosphere
  const auto& atmo = det::Detector::GetInstance().GetAtmosphere();

  const double dEdX0 = atmo.GetdEdX0();

  // vector of wave lengths
  const auto& lambdaVec = fFWaveLength;

  diagonalMatrix& flMat = *fFluorescenceMatrix;

  for (unsigned int i = 0, n = fShowerPoints.size(); i < n; ++i) {
    const utl::Point& meanPos = fShowerPoints[i];

    // get vector from point to eye
    const double showerAge = utl::ShowerAge(fDepth[i], fXmax);
    const double disteye = (meanPos - fEyePosition).GetMag();
    const double spotFraction = fOpticalHalo.GetHaloFraction(fZeta, showerAge, disteye);

    double height = 0;

    try {
      utl::UTMPoint UTMp(meanPos, utl::ReferenceEllipsoid::GetWGS84());
      height = UTMp.GetHeight();
    } catch (utl::UTMPoint::UTMZoneException& d) {
      cout << d.GetExceptionName() << endl;
      const ProfileResult& densityProfile = atmo.EvaluateDensityVsHeight();
      height = densityProfile.MaxX();
    }

    const auto& fluoyield_tab = atmo.EvaluateFluorescenceYield(height);

    double epsYTsum = 0;
    for (unsigned int j = 0, n = lambdaVec.size(); j < n; ++j) {
      const double y_j = fluoyield_tab.Y(lambdaVec[j]) * fFfac;
      epsYTsum += y_j * fRelEff.Y(lambdaVec[j]) * fFT2eye[i][j];
    }

    const double msFactor =  1 / (1 - MultipleScatteringFraction(i, lambdaVec, fluoyield_tab));
    fFluorescenceMultipleScattering.push_back(msFactor);

    const double deltaL = (fShowerGeometry[i][0] - fShowerGeometry[i][1]).GetMag();
    const double zetaFraction = FluorescenceLDFFraction(i);
    flMat(i, i) = msFactor * spotFraction * zetaFraction * epsYTsum * fGfac[i] / dEdX0 * deltaL;
  }
}


// gets wavelengths to be used for Cherenkov and Fluorescence

void
CherenkovFluorescenceMatrix::CalculateWaveLengths()
{
#ifdef _CFM_DEBUG
  cout << "        -->GetWaveLengths()" << endl;
#endif

  // current atmosphere
  const auto& atmo = det::Detector::GetInstance().GetAtmosphere();

  fCWaveLength = atmo.GetWavelengths(Atmosphere::eCerenkov);
  fFWaveLength = atmo.GetWavelengths(Atmosphere::eFluorescence);

  // remove wave lengths with no or zero efficiency
  const double effMinLam = fRelEff[0].X();
  const double effMaxLam = fRelEff[fRelEff.GetNPoints()-1].X();

  for (auto it = fFWaveLength.begin(); it != fFWaveLength.end(); ) {
    const auto& lam = *it;
    if (effMinLam <= lam && lam <= effMaxLam && fRelEff.Y(lam) > 0)
      ++it;
    else
      fFWaveLength.erase(it);
  }

  for (auto it = fCWaveLength.begin(); it != fCWaveLength.end(); ) {
    const auto& lam = *it;
    if (effMinLam <= lam && lam <= effMaxLam)
      //&& fRelEff.Y(*it) > 0) // !!!C-bins are integration bounds!!!
      ++it;
    else
      fCWaveLength.erase(it);
  }
}


void
CherenkovFluorescenceMatrix::CalculateDirectCherenkovMatrix()
{
#ifdef _CFM_DEBUG
  cout << "        -->CalculateDirectCherenkovMatrix()" << endl;
#endif

  fCherenkovMultipleScattering.clear();

  // current atmosphere
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  // vector of wave lengths
  const std::vector<double>& lambdaVec = fCWaveLength;

  diagonalMatrix& cDirMat = *fDirectCherenkovMatrix;

  for (unsigned int i = 0; i < fShowerPoints.size(); ++i) {
    const double showerAge = utl::ShowerAge(fDepth[i], fXmax);

    // get vector from point to eye
    const utl::Point& meanPos = fShowerPoints[i];
    const double disteye = (meanPos - fEyePosition).GetMag();
    const double spotFraction =
      fOpticalHalo.GetHaloFraction(fZeta, showerAge, disteye);

    const utl::TabulatedFunction& directCherenkov =
      atmo.EvaluateCherenkovDirect(fShowerGeometry[i][0],
                                   fShowerGeometry[i][1],
                                   fEyePosition,showerAge);
    double epsYTsum = 0;
    for (unsigned int j = 0; j < lambdaVec.size(); ++j) {
      const double y_j = directCherenkov.GetY(j) * fCfac;
      epsYTsum += y_j * fRelEff.Y(lambdaVec[j]) * fCT2eye[i][j];
    }

    const double msFactor =
      1 / (1 - MultipleScatteringFraction(i, lambdaVec, directCherenkov));
    fCherenkovMultipleScattering.push_back(msFactor);
    const double zetaFraction = DirectCherenkovLDFFraction(i);
    cDirMat(i, i) = msFactor * spotFraction * zetaFraction * epsYTsum / fMeandEdX[i];
  }
}


void
CherenkovFluorescenceMatrix::CalculateMieAndRayScattCherenkovMatrix()
{
#ifdef _CFM_DEBUG
  cout << "        -->CalculateMieAndRayScattCherenkovMatrix()" << endl;
#endif

  // make sure CalculateDirectCherenkovMatrix() was
  // called before ....
  if (fCherenkovMultipleScattering.size() != fShowerPoints.size()) {
    cerr << " CalculateMieAndRayScattCherenkovMatrix() "
         << " Error - no multiple scattering factors " << endl;
    throw utl::OutOfBoundException();
  }

  // vector of wave lengths
  const std::vector<double>& lambdaVec = fCWaveLength;
  const unsigned int nLambda = lambdaVec.size();
  double t[nLambda];

  lowerTriangularMatrix& cMieMat = *fMieScatteredCherenkovMatrix;
  lowerTriangularMatrix& cRayMat = *fRayleighScatteredCherenkovMatrix;

  for (unsigned int i = 0; i < fShowerPoints.size(); ++i) {

    const double showerAge = utl::ShowerAge(fDepth[i],fXmax);

    // get vector from point to eye
    const utl::Point& meanPos = fShowerPoints[i];
    const double disteye = (meanPos - fEyePosition).GetMag();
    const double spotFraction =
      fOpticalHalo.GetHaloFraction(fZeta, showerAge, disteye);

    const double msFactor = fCherenkovMultipleScattering[i];

    for (unsigned int k = 0; k < nLambda; ++k)
      t[k] = 1;

    for (int j = i; j >= 0; --j) {
      double epsYTsumMie = 0;
      double epsYTsumRay = 0;
      for (unsigned int k = 0; k < nLambda; ++k) {
        if (j != int(i))
          t[k] *= fTShower[j][k];

        const double detEff = fRelEff.Y(lambdaVec[k]);
        const double factor = t[k] * detEff * fCherAtTrack[j][k] * fCT2eye[i][k];
        epsYTsumMie += factor * fMieScat2eye[i][k];
        epsYTsumRay += factor * fRayScat2eye[i][k];
      }
      const double zetaFraction = ScatteredCherenkovLDFFraction(i, j);
      const double zetaPerAlpha = zetaFraction / fMeandEdX[j];
      cMieMat(i, j) = msFactor * spotFraction * zetaPerAlpha * epsYTsumMie;
      cRayMat(i, j) = msFactor * spotFraction * zetaPerAlpha * epsYTsumRay;
    }

  }
}


// Mie- and Rayleigh-attenuation from shower to eye
void
CherenkovFluorescenceMatrix::CalculateAttenuationToEye()
{
#ifdef _CFM_DEBUG
  cout << "        -->CalculateAttenuationToEye()" << endl;
#endif

  fCT2eye.clear();
  fFT2eye.clear();

  // current atmosphere
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  // vector of C/F wave lengths
  const std::vector<double>& lambdaCVec = fCWaveLength;
  const unsigned int nCLambda = lambdaCVec.size();
  std::vector<double> tC(nCLambda);
  const std::vector<double>& lambdaFVec = fFWaveLength;
  const unsigned int nFLambda = lambdaFVec.size();
  std::vector<double> tF(nFLambda);

  for (unsigned int i = 0; i < fShowerPoints.size(); ++i) {

    const utl::Point& meanPos = fShowerPoints[i];

    // ---- for Ch. wave langth

    // Mie attenuation
    const AttenuationResult mCAtt =
      atmo.EvaluateMieAttenuation(fEyePosition, meanPos, lambdaCVec);
    const utl::TabulatedFunctionErrors& mieCAtt = mCAtt.GetTransmissionFactor();
    // Rayleigh attenuation
    const AttenuationResult rCAtt =
      atmo.EvaluateRayleighAttenuation(fEyePosition, meanPos, lambdaCVec);
    const utl::TabulatedFunctionErrors& rayCAtt = rCAtt.GetTransmissionFactor();

    for (unsigned int j = 0; j < nCLambda; ++j) {
      tC[j] = mieCAtt.GetY(j) * rayCAtt.GetY(j);
    }
    fCT2eye.push_back(tC);

    // ---- for Fl. wave langth

    // Mie attenuation
    const AttenuationResult mFAtt =
      atmo.EvaluateMieAttenuation(fEyePosition, meanPos, lambdaFVec);
    const utl::TabulatedFunctionErrors& mieFAtt = mFAtt.GetTransmissionFactor();
    // Rayleigh attenuation
    const AttenuationResult rFAtt =
      atmo.EvaluateRayleighAttenuation(fEyePosition, meanPos, lambdaFVec);
    const utl::TabulatedFunctionErrors& rayFAtt = rFAtt.GetTransmissionFactor();

    for (unsigned int j = 0; j < nFLambda; ++j) {
      tF[j] = mieFAtt.GetY(j) * rayFAtt.GetY(j);
    }
    fFT2eye.push_back(tF);

  }
}


// Mie- and Rayleigh-attenuation along the track

void
CherenkovFluorescenceMatrix::CalculateAttenuationAlongTrack()
{
#ifdef _CFM_DEBUG
  cout << "        -->CalculateAttenuationAlongTrack()" << endl;
#endif

  fTShower.clear();

  // current atmosphere
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  // vector of wave lengths
  const std::vector<double>& lambdaVec = fCWaveLength;
  const unsigned int nLambda = lambdaVec.size();
  std::vector<double> transmission(nLambda);

  for (unsigned int i = 0; i < fShowerPoints.size()-1; ++i) {

    const utl::Point& meanPos_i = fShowerPoints[i];
    const int j = i + 1;
    const utl::Point& meanPos_j = fShowerPoints[j];

    // Mie attenuation
    const AttenuationResult mAtt =
      atmo.EvaluateMieAttenuation(meanPos_i, meanPos_j, lambdaVec);
    const utl::TabulatedFunctionErrors& mieAtt = mAtt.GetTransmissionFactor();

    // Rayleigh attenuation
    const AttenuationResult rAtt =
      atmo.EvaluateRayleighAttenuation(meanPos_i, meanPos_j, lambdaVec);
    const utl::TabulatedFunctionErrors& rayAtt = rAtt.GetTransmissionFactor();

    for (unsigned int k = 0; k < nLambda; ++k) {
      transmission[k] = mieAtt.GetY(k) * rayAtt.GetY(k);
    }
    fTShower.push_back(transmission);

  }
}


// Mie- and Rayleigh scattering from shower to eye

void
CherenkovFluorescenceMatrix::CalculateScatteringToEye()
{
#ifdef _CFM_DEBUG
  cout << "        -->CalculateScatteringToEye()" << endl;
#endif

  fMieScat2eye.clear();
  fRayScat2eye.clear();

  // current atmosphere
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  // vector of wave lengths
  const std::vector<double>& lambdaVec = fCWaveLength;
  const unsigned int nLambda = lambdaVec.size();
  std::vector<double> tR(nLambda);
  std::vector<double> tM(nLambda);

  // unit vector along shower
  const utl::Vector shwDir = fShowerGeometry[0][1] - fShowerGeometry[0][0];

  for (unsigned int i = 0; i < fShowerPoints.size(); ++i) {
    const utl::Point& meanPos = fShowerPoints[i];

    // average vector from shower to eye
    utl::Vector meanVec = meanPos - fEyePosition;
    const double disteye = meanVec.GetMag();

    // emission angle.
    const double angle = Angle(-shwDir,meanVec);

    const ScatteringResult& rScattering =
      atmo.EvaluateRayleighScattering(fShowerGeometry[i][0],
                                      fShowerGeometry[i][1],
                                      angle, disteye, lambdaVec);
    const ScatteringResult& mScattering =
      atmo.EvaluateMieScattering(fShowerGeometry[i][0],
                                 fShowerGeometry[i][1],
                                 angle, disteye, lambdaVec);
    const utl::TabulatedFunctionErrors& rayleighScattTrackEye =
      rScattering.GetScatteringFactor();
    const utl::TabulatedFunctionErrors& mieScattTrackEye =
      mScattering.GetScatteringFactor();
    for (unsigned int j = 0; j < nLambda; ++j) {
      tR[j] = rayleighScattTrackEye.GetY(j);
      tM[j] = mieScattTrackEye.GetY(j);
    }
    fRayScat2eye.push_back(tR);
    fMieScat2eye.push_back(tM);

  }
}


// Cherenkov light per electron produced at track

void
CherenkovFluorescenceMatrix::CalculateCherenkovAtTrack()
{
#ifdef _CFM_DEBUG
  cout << "        -->CalculateCherenkovAtTrack()" << endl;
#endif

  // current atmosphere
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  fCherAtTrack.clear();

  // vector of wave lengths
  const std::vector<double>& lambdaVec = fCWaveLength;
  const unsigned int nLambda = lambdaVec.size();
  std::vector<double> t(nLambda);

  for (unsigned int i = 0; i < fShowerPoints.size(); ++i) {
    const double showerAge = utl::ShowerAge(fDepth[i], fXmax);
    const utl::TabulatedFunction& CherenkovProd =
      atmo.EvaluateCherenkovPhotons(fShowerGeometry[i][0], fShowerGeometry[i][1], showerAge);
    for (unsigned int j = 0; j < nLambda; ++j) {
      t[j] = CherenkovProd.GetY(j) * fCfac;
    }
    fCherAtTrack.push_back(t);
  }
}


// Calculates fraction of light 1/(4*pi*r^2) received by telescope

void
CherenkovFluorescenceMatrix::CalculateGeometricalFactor()
{
#ifdef _CFM_DEBUG
  cout << "        -->CalculateGeometricalFactor()" << endl;
#endif

  fGfac.clear();
  const double fourPi = 4*utl::kPi;
  for (unsigned int i = 0; i < fShowerPoints.size(); ++i) {
    const utl::Vector r = fShowerPoints[i] - fEyePosition;
    const double rSquared = r.GetMag2();
    fGfac.push_back(1/(fourPi*rSquared));
  }
}


// Calculates the mean energy deposit per electron

void
CherenkovFluorescenceMatrix::CalculateMeanEnergyDeposit()
{
  fMeandEdX.clear();
  for (unsigned int i = 0; i < fShowerPoints.size(); ++i) {
    const double dEdXperElectron = utl::EnergyDeposit(fDepth[i], fXmax, fEcut);
    fMeandEdX.push_back(dEdXperElectron);
  }
}


// Calculates sum of Cherenkov and fluorescence matrixes and returns pointer to total CherenkovFluorescenceMatrix

const lowerTriangularMatrix*
CherenkovFluorescenceMatrix::GetCherenkovFluorescenceMatrix()
  const
{
#ifdef _CFM_DEBUG
  cout << "        -->GetCherenkovFluorescenceMatrix()" << endl;
#endif
  if (!fDirectCherenkovMatrix ||
      !fMieScatteredCherenkovMatrix ||
      !fRayleighScatteredCherenkovMatrix ||
      !fFluorescenceMatrix ||
      !fCherenkovFluorescenceMatrix)
    return nullptr;

  if (!fUpToDate) {
    const diagonalMatrix& fl = *fFluorescenceMatrix;
    const diagonalMatrix& dc = *fDirectCherenkovMatrix;
    const lowerTriangularMatrix& rc = *fRayleighScatteredCherenkovMatrix;
    const lowerTriangularMatrix& mc = *fMieScatteredCherenkovMatrix;
    lowerTriangularMatrix& tot = *fCherenkovFluorescenceMatrix;

    for (unsigned int i = 0, n = fShowerPoints.size(); i < n; ++i) {
      tot(i, i) = fl(i, i) + dc(i, i) + rc(i, i) + mc(i, i);
      for (unsigned int j = 0; j < i; ++j)
        tot(i, j) = rc(i, j) + mc(i, j);
    }
    fUpToDate = true;
  }
  return fCherenkovFluorescenceMatrix;
}


// recalculate matrices depending on shower age

void
CherenkovFluorescenceMatrix::UpdateXmax(const double xMax)
{
#ifdef _CFM_DEBUG
  cout << "        -->updateXmax()" << endl;
#endif
  fXmax = xMax;

  CalculateMeanEnergyDeposit();
  CalculateCherenkovAtTrack();
  CalculateFluorescenceMatrix();
  CalculateDirectCherenkovMatrix();
  CalculateMieAndRayScattCherenkovMatrix();

  // flag GetCherenkovFluorescenceMatrix()
  // that it needs to recalculate
  fUpToDate = false;
}


// fraction of direct Cherenkov light shower image

double
CherenkovFluorescenceMatrix::DirectCherenkovLDFFraction(const int iPosition)
  const
{
  switch (fDirCherLDF) {
  case eNoDirCherLDF:
    return 1;
  case eGoraDirCherLDF:
    return GoraFraction(iPosition);
  }
  return 1;
}


// fraction of scattered Cherenkov light shower image

double
CherenkovFluorescenceMatrix::ScatteredCherenkovLDFFraction(const int iObserved,
                                                           const int jOrigin)
  const
{
  switch (fScatCherLDF) {
  case eNoScatCherLDF:
    return 1;
  case eGoraScatCherLDF:
    return GoraFraction(iObserved);
  case eExpoScatCherLDF:
    return ExponentialFraction(iObserved);
  case eGillerScatCherLDF:
    return GillerFraction(iObserved,jOrigin);
  }
  return 1;
}


// fraction of lateral fluorescence light shower image

double
CherenkovFluorescenceMatrix::FluorescenceLDFFraction(const int iPosition)
  const
{
  switch (fFluoLDF) {
  case eNoFluoLDF:
    return 1;
  case eGoraFluoLDF:
    return GoraFraction(iPosition);
  }
  return 1;
}


/* fraction of scattered Cherenkov light
   B. Dawson, M. Giller, G. Wieczorek, ICRC07 Merida, #0651 */

double
CherenkovFluorescenceMatrix::GillerFraction(const int iObserved,
                                            const int jOrigin)
  const
{
  const double pathLength = (fShowerGeometry[iObserved][0] -
                             fShowerGeometry[jOrigin][0]).GetMag();
  const double zetaDistance = GetZetaDistance(iObserved);

  const double maxEmissionAngle = atan2(zetaDistance, pathLength);

  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  const ProfileResult& depthProfile = atmo.EvaluateDepthVsHeight();

  double height = 0;
  try {
    const utl::UTMPoint UTMp(fShowerPoints[jOrigin],
                       utl::ReferenceEllipsoid::GetWGS84());
    height = UTMp.GetHeight();
  } catch (utl::UTMPoint::UTMZoneException& d) {
    cout << d.GetExceptionName() << endl;
    height = depthProfile.MaxX();
  }

  const double xVert = depthProfile.Y(height);
  const double xSlant = fDepth[jOrigin];
  const double showerAge = utl::ShowerAge(xSlant, fXmax);

  const double minValidCherenkovAngle = 5*utl::degree;

  const double fraction =
    (!fCherenkovCone || maxEmissionAngle>minValidCherenkovAngle) ?
      atmo.AngularCherenkovCDF(maxEmissionAngle,xVert,showerAge) :
      AngularCherenkovCDFWithCone(maxEmissionAngle, xVert, showerAge);

  return fraction;
}


// modified angular Cherenkov CDF taking into account the Cherenkov cone

double
CherenkovFluorescenceMatrix::AngularCherenkovCDFWithCone(const double maxEmissionAngle,
                                                         const double xVert, const double showerAge)
  const
{
  // number of MC evaluations
  const int nMC = 1000;

  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  // Cherenkov angle for beta=1
  const double nRefrac = utl::RefractionIndex::LorentzLorentz(xVert);
  const double cosThetaCher = 1/nRefrac;
  const double sinThetaCher = sqrt(1 - cosThetaCher*cosThetaCher);

  // array of angular CDF for random number generation
  // (yes, Hillas&Nerling are exponentials, but are we supposed to know
  //  this outside of utl::Atmosphere??)

  const double angularBinning = 0.25*utl::degree;
  std::vector<double> cdf;

  const double maxIntegrationAngle = 90*utl::degree;
  double thisAngle = 0;

  while (thisAngle <= maxIntegrationAngle) {
    const double cumulative = atmo.AngularCherenkovCDF(thisAngle, xVert, showerAge);
    cdf.push_back(cumulative);
    thisAngle += angularBinning;
  }

  utl::RandomEngine& theRandom =
    RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::ePhysics);

  std::vector<double> randomNumbers1(nMC);
  std::vector<double> randomNumbers2(nMC);

  RandFlat::shootArray(&theRandom.GetEngine(), nMC, &randomNumbers1.front());
  RandFlat::shootArray(&theRandom.GetEngine(), nMC, &randomNumbers2.front());

  // calculate fraction of photons within [maxEmissionAngle]

  const double cosThetaEmission = cos(maxEmissionAngle);

  int nInsideZeta = 0;

  for (int i = 0; i < nMC; ++i) {
    const double x = randomNumbers1[i];
    double alpha = maxIntegrationAngle;
    for (unsigned int k = 0, n = cdf.size()-1; k < n; ++k) {
      if (x > cdf[k] && x < cdf[k+1]) {
        alpha = k*angularBinning + angularBinning/(cdf[k] - cdf[k+1]) * (cdf[k] - x);
        break;
      }
    }

    const double phi = randomNumbers2[i]*utl::kTwoPi;
    const double cosTheta = sinThetaCher*cos(phi)*sin(alpha) + cos(alpha)*cosThetaCher;

    if (cosTheta > cosThetaEmission)
      ++nInsideZeta;
  }

  const double fractionWithinEmissionAngle = double(nInsideZeta) / nMC;

  //  cout << fractionWithinEmissionAngle
  //       << " " << atmo.AngularCherenkovCDF(maxEmissionAngle,xVert,showerAge)
  //       << " " << maxEmissionAngle/utl::degree
  //       << " " << acos(cosThetaCher)/utl::degree  << endl;

  return fractionWithinEmissionAngle;
}


/* fraction of multiple scattered light wrt. total signal from
   M. Roberts, J. Phys. G 31 (2005) 1291 */

double
CherenkovFluorescenceMatrix::RobertsFraction(const int iPosition,
                                             const std::vector<double>& waveLengths,
                                             const utl::TabulatedFunction& yield)
  const
{
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  // optical depth and scattering coefficient

  //         attenuation shower-->eye
  const utl::Point& meanPos = fShowerPoints[iPosition];
  const AttenuationResult mieResultAttShowerToEye =
    atmo.EvaluateMieAttenuation(fEyePosition, meanPos, waveLengths);
  const utl::TabulatedFunctionErrors& mieAttShowerToEye =
    mieResultAttShowerToEye.GetTransmissionFactor();

  const AttenuationResult rayAttResultShowerToEye =
    atmo.EvaluateRayleighAttenuation(fEyePosition, meanPos, waveLengths);
  const utl::TabulatedFunctionErrors& rayAttShowerToEye =
    rayAttResultShowerToEye.GetTransmissionFactor();

  //         attenuation within bin
  const AttenuationResult mieAttResultInBin =
    atmo.EvaluateMieAttenuation(fShowerGeometry[iPosition][0],
                                fShowerGeometry[iPosition][1],
                                waveLengths);
  const utl::TabulatedFunctionErrors& mieAttInBin =
    mieAttResultInBin.GetTransmissionFactor();

  const AttenuationResult rayAttResultInBin =
    atmo.EvaluateRayleighAttenuation(fShowerGeometry[iPosition][0],
                                     fShowerGeometry[iPosition][1],
                                     waveLengths);
  const utl::TabulatedFunctionErrors& rayAttInBin =
    rayAttResultInBin.GetTransmissionFactor();

  double yieldEffSum       = 0;
  double yieldEffAttEyeSum = 0;
  double yieldEffAttBinSum = 0;

  for (unsigned int j = 0, n = waveLengths.size(); j < n; ++j) {
    const double y_j = yield.Y(waveLengths[j]);
    const double eff = fRelEff.Y(waveLengths[j]);
    const double attEye = rayAttShowerToEye.GetY(j) * mieAttShowerToEye.GetY(j);
    const double attBin = rayAttInBin.GetY(j) * mieAttInBin.GetY(j);

    const double yEff  = y_j * eff;
    yieldEffSum += yEff;
    yieldEffAttEyeSum += yEff * attEye;
    yieldEffAttBinSum += yEff * attBin;
  }

  const double distanceToEye = (fShowerPoints[iPosition] - fEyePosition).GetMag();
  const double binLength = (fShowerGeometry[iPosition][0] -
                            fShowerGeometry[iPosition][1]).GetMag();

  const double attenuationToEye = yieldEffAttEyeSum / yieldEffSum;
  const double attenuationInBin = yieldEffAttBinSum / yieldEffSum;

  const double opticalDepth = attenuationToEye > 0 ?
    -log(attenuationToEye) :
    std::numeric_limits<double>::max();
  const double alphaTot = attenuationInBin > 0 ?
    -log(attenuationInBin)/binLength :
    std::numeric_limits<double>::max();

  // see Eq. (3) in Mike's paper
  const double argument = opticalDepth * alphaTot/utl::m *
                          sqrt(distanceToEye/utl::m) *
                          pow(fZeta/utl::degree, 1.1);
  const double fraction = 0.774 * pow(argument, 0.68);

  return fraction;
}


/* fraction of multiple scattered light wrt. total signal from
   Jan Pekala et al., GAP-2008-087 */

double
CherenkovFluorescenceMatrix::PekalaFraction(const int iPosition,
                                            const std::vector<double>& waveLengths,
                                            const utl::TabulatedFunction& yield)
  const
{
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  // optical depth

  //         attenuation shower-->eye
  const utl::Point& meanPos = fShowerPoints[iPosition];

  const AttenuationResult mieResultAttShowerToEye =
    atmo.EvaluateMieAttenuation(fEyePosition, meanPos, waveLengths);
  const utl::TabulatedFunctionErrors& mieAttShowerToEye =
    mieResultAttShowerToEye.GetTransmissionFactor();

  const AttenuationResult rayAttResultShowerToEye =
    atmo.EvaluateRayleighAttenuation(fEyePosition, meanPos, waveLengths);
  const utl::TabulatedFunctionErrors& rayAttShowerToEye =
    rayAttResultShowerToEye.GetTransmissionFactor();

  double yieldEffSum       = 0;
  double yieldEffAttEyeSum = 0;

  for (unsigned int j = 0, n = waveLengths.size(); j < n; ++j) {
    const double y_j = yield.Y(waveLengths[j]);
    const double eff = fRelEff.Y(waveLengths[j]);
    const double attEye = rayAttShowerToEye.GetY(j) * mieAttShowerToEye.GetY(j);

    const double yEff  = y_j * eff;
    yieldEffSum += yEff;
    yieldEffAttEyeSum += yEff * attEye;
  }

  const double attenuationToEye = yieldEffAttEyeSum / yieldEffSum;

  const double opticalDepth = attenuationToEye > 0 ?
    -log(attenuationToEye) :
    std::numeric_limits<double>::max();

  double height;
  try {
    utl::UTMPoint UTMp(fShowerPoints[iPosition],utl::ReferenceEllipsoid::GetWGS84());
    height = UTMp.GetHeight();
  } catch (utl::UTMPoint::UTMZoneException& d) {
    cout << d.GetExceptionName() << endl;
    height = atmo.EvaluateDensityVsHeight().MaxX();
  }

  const double g = 5.426*utl::km;
  //  const double f = 3.318e-2;
  // height above array --> height a.s.l. exp(1.4/G) = 1.2943
  const double f = 4.2947e-2;

  // fraction wrt. direct light
  const double fraction1 = f * opticalDepth * (fZeta/utl::degree) * exp(-height/g);

  // fraction wrt. total light
  const double fraction2 = fraction1 / (fraction1 + 1);

  return fraction2;
}


/* fraction of lateral shower image according to phenomenological
   LDF proposed by Spot group
   (see talk by M. Tueros, Malargue meeting Nov 2007) */

double
CherenkovFluorescenceMatrix::ExponentialFraction(const int iPosition)
  const
{
  const double r = GetZetaDistance(iPosition);
  const double alpha = 1.65;
  const double beta = 325.*utl::m;
  const double fraction = 1 - exp(-pow(r/beta, alpha));

  return fraction;
}


/* fraction of lateral shower image according to
   Gora et. al Astropart. Phys. 24(2006) 484. */

double
CherenkovFluorescenceMatrix::GoraFraction(const int iPosition)
  const
{
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  const ProfileResult& tempProfile  = atmo.EvaluateTemperatureVsHeight();
  const ProfileResult& depthProfile = atmo.EvaluateDepthVsHeight();
  const ProfileResult& pressureProfile = atmo.EvaluatePressureVsHeight();

  // cos(zenith) for moliere radius
  const utl::CoordinateSystemPtr cs = det::Detector::GetInstance().GetSiteCoordinateSystem();
  const utl::Vector shwDir = fShowerGeometry[0][0] - fShowerGeometry[0][1];
  const double cosTheta = shwDir.GetCosTheta(cs);

  const utl::Point& meanPos = fShowerPoints[iPosition];

  double height;
  try {
    utl::UTMPoint UTMp(meanPos,utl::ReferenceEllipsoid::GetWGS84());
    height = UTMp.GetHeight();
  } catch (utl::UTMPoint::UTMZoneException& d) {
    cout << d.GetExceptionName() << endl;
    height = depthProfile.MaxX();
  }

  const double xSlant = fDepth[iPosition];
  const double pressure = pressureProfile.Y(height);
  const double temperature = tempProfile.Y(height);
  const double age = utl::ShowerAge(xSlant, fXmax);

  const double rm = utl::MoliereRadius(temperature, pressure, cosTheta);

  const double rStar = GetZetaDistance(iPosition) / rm;

  return utl::GoraCDF(rStar, age);
}


// distance to shower axis corresponding to zeta

double
CherenkovFluorescenceMatrix::GetZetaDistance(const int iPosition)
  const
{
  const utl::Point& meanPos = fShowerPoints[iPosition];
  const double r = (meanPos - fEyePosition).GetMag() * fTanZeta;
  return r;
}


double
CherenkovFluorescenceMatrix::MultipleScatteringFraction(const int iPosition,
                                                        const std::vector<double>& waveLengths,
                                                        const utl::TabulatedFunction& yield)
  const
{
  switch (fMultScatLDF) {
  case eNoMultScatLDF:
    return 0;
  case eRobertsLDF:
    return RobertsFraction(iPosition, waveLengths, yield);
  case ePekalaLDF:
    return PekalaFraction(iPosition, waveLengths, yield);
  }
  return 0;
}


// set a new zeta value and recalculate width dependend stuff

void
CherenkovFluorescenceMatrix::SetZeta(const double zeta)
{
  fUpToDate = false;
  fTanZeta = tan(zeta);
  fZeta = zeta;

  CalculateFluorescenceMatrix();
  CalculateDirectCherenkovMatrix();
  CalculateMieAndRayScattCherenkovMatrix();
}


// rescale matrix elements for new yield factors

void
CherenkovFluorescenceMatrix::SetYieldFactors(const double fluoFac,
                                             const double chkovFac)
{
  diagonalMatrix& fl = *fFluorescenceMatrix;
  diagonalMatrix& dc = *fDirectCherenkovMatrix;
  lowerTriangularMatrix& rc = *fRayleighScatteredCherenkovMatrix;
  lowerTriangularMatrix& mc = *fMieScatteredCherenkovMatrix;
  lowerTriangularMatrix& tot = *fCherenkovFluorescenceMatrix;

  const double fluoRescale = fluoFac / fFfac;
  const double chkovRescale = chkovFac / fCfac;

  for (unsigned int i = 0, n = fShowerPoints.size(); i < n; ++i) {
    fl(i, i) *= fluoRescale;
    dc(i, i) *= chkovRescale;
    rc(i, i) *= chkovRescale;
    mc(i, i) *= chkovRescale;
    tot(i, i) = fl(i, i) + dc(i, i) + rc(i, i) + mc(i, i);
    for (unsigned int j = 0; j < i; ++j) {
      rc(i, j) *= chkovRescale;
      mc(i, j) *= chkovRescale;
      tot(i, j) = rc(i, j) + mc(i, j);
    }
  }

  fCfac = chkovFac;
  fFfac = fluoFac;
}


/* print extensive informations about current data members (for debugging)
   iwhat=0  -> print transmission to eye
        =1  -> print transmission along track */

void
CherenkovFluorescenceMatrix::Print(const int iwhat)
  const
{
  boost::io::ios_all_saver save(cout);
  cout.flags(std::ios::scientific);
  cout.precision(2);

  // print transmission to eye
  const std::vector<double>& lambdaVec = fCWaveLength;
  const unsigned int nLambda = lambdaVec.size();

  if (!iwhat) {
    cout << "\n  wave lengths:\n";
    for (unsigned int j = 0; j < nLambda; ++j)
      cout << std::setw(9) << lambdaVec[j];
    cout << endl;

    cout << "\n  transmission to eye:\n";
    for (unsigned int i = 0; i < fShowerPoints.size(); ++i) {
      for (unsigned int j = 0; j < nLambda; ++j)
        cout << std::setw(9) << fCT2eye[i][j];
      cout << endl;
    }
  }  else if (iwhat == 1) {
    const int kLam = 9;
    for (unsigned int i = 0; i < fShowerPoints.size(); ++i) {
      cout << std::setw(9) << i
           << std::setw(9) << lambdaVec[kLam]
           << std::setw(9) << fCherAtTrack[i][kLam]
           << std::setw(9) << fCT2eye[i][kLam]
           << std::setw(9) << fMieScat2eye[i][kLam]
           << endl;
    }
  }
}