ShowerLightSimulator.cc

Go to the documentation of this file.

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/UTMPoint.h>
#include <utl/Particle.h>
#include <utl/GeometryUtilities.h>

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

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

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>

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

#include <CLHEP/Random/Randomize.h>

#include "ShowerLightSimulator.h"

#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_io.hpp>

using namespace ShowerLightSimulatorKG;
using namespace std;
using namespace utl;
using namespace evt;
using namespace fwk;
using namespace det;
using namespace fdet;
using namespace atm;


VModule::ResultFlag
ShowerLightSimulator::Init()
{
  const Branch topB = CentralConfig::GetInstance()->GetTopBranch("ShowerLightSimulatorKG");

  fPrimaryCherenkov = false;
  topB.GetChild("fluorescence").GetData(fFluorescence);
  topB.GetChild("cherenkov").GetData(fCherenkov);
  topB.GetChild("cherenkovFromCORSIKA").GetData(fUseCherenkovProfile);
  topB.GetChild("flatEarth").GetData(fFlatEarth);
  topB.GetChild("timeBin").GetData(fTimeBin);

  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "         fluorescence: " << (fFluorescence ? "ON" : "OFF") << "\n"
          "            cherenkov: " << (fCherenkov ? "ON" : "OFF") << "\n"
          "       cherenkov mode: " << (fUseCherenkovProfile ? "CORSIKA" : "model") << "\n"
          "  atmosphere geometry: " << (fFlatEarth ? "FLAT" : "CURVED") << '\n';

  INFO(info);

  return eSuccess;
}


VModule::ResultFlag
ShowerLightSimulator::Run(evt::Event& event)
{
  if (!event.HasSimShower()) {
    ERROR("Event has no simulated shower.");
    return eFailure;
  }

  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);
  Detector& theDetector = Detector::GetInstance();

  const evt::ShowerSimData& simShower = event.GetSimShower();

  CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  Point showerCore(0, 0, 0, showerCS);
  Vector showerDir(0, 0, -1, showerCS);

  // just needed for flat geometry
  const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
  const double zenith = (-simShower.GetDirection()).GetTheta(localCS);
  const double cosTheta = cos (zenith); // JUST FOR LIMITS
  const double absCosTheta = std::fabs(cosTheta);
  const double coreZ = wgs84.PointToLatitudeLongitudeHeight (showerCore).get<2>(); // JUST FOR FLAT-EARTH

  //double cosTheta = std::fabs(cos (zenith)); // JUST FOR LIMITS
  //double Xfirst = simShower.GetXFirst();

  const ProfileResult& heightVsDepth = theDetector.GetAtmosphere().EvaluateHeightVsDepth();

  const ProfileResult* distanceVsSlantDepth = 0;
  double atmDepthMin = 0;
  double atmDepthMax = 0;
  double atmDistanceMin = 0;
  double atmDistanceMax = 0;

  bool upward = false;
  bool skimming = false;

  if (!fFlatEarth) {

    try {

      // theDir has to point to where the shower is comming from !!!!
      // this is now done in FieldOfViewCalculator:
      // theDetector.GetAtmosphere().InitSlantProfileModel(showerCore, -showerDir, 10.*g/cm2);

      distanceVsSlantDepth = &theDetector.GetAtmosphere().EvaluateDistanceVsSlantDepth();


      // get the LIMITS of the slant atmosphere
      atmDepthMin = distanceVsSlantDepth->MinX();
      atmDepthMax = distanceVsSlantDepth->MaxX();
      atmDistanceMin = distanceVsSlantDepth->Y(atmDepthMin);
      atmDistanceMax = distanceVsSlantDepth->Y(atmDepthMax);

      const ProfileResult& slantDepthVsDistance = theDetector.GetAtmosphere().EvaluateSlantDepthVsDistance();
      const ProfileResult& heightVsDistance = theDetector.GetAtmosphere().EvaluateHeightVsDistance();
      upward = (heightVsDistance.Y(atmDistanceMin) < heightVsDistance.Y(atmDistanceMax) - 1.*m); // security margin
      skimming = ((upward ? heightVsDistance.Y(atmDistanceMin) : heightVsDistance.Y(atmDistanceMax)) >
                   heightVsDistance.Y((atmDistanceMax+atmDistanceMin)/2));
      double atmHeightMax = (upward ? heightVsDistance.Y(atmDistanceMax) : heightVsDistance.Y(atmDistanceMin));
      double atmHeightMin = (skimming ? heightVsDistance.Y((atmDistanceMax+atmDistanceMin)/2) :
                              (upward ? heightVsDistance.Y(atmDistanceMin) : heightVsDistance.Y(atmDistanceMax)));

      // check against limits of vertical atmosphere
      const ProfileResult& temperaturProfile = theDetector.GetAtmosphere().EvaluateTemperatureVsHeight();
      const ProfileResult& densityProfile = theDetector.GetAtmosphere().EvaluateDensityVsHeight();

      if (atmHeightMax>temperaturProfile.MaxX() ||
          atmHeightMax>densityProfile.MaxX()) {

        atmHeightMax = std::min(temperaturProfile.MaxX(), densityProfile.MaxX()) - 1.*mm; // security margin
        const vector<Point> intersections =
          Intersection(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84), atmHeightMax,
                       Line(showerCore, -showerDir));
        if (intersections.size() != 2) {
          ostringstream warn;
          warn << "\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
                  "Number of intersections point of air shower axis with ellipsoid of height="
               << atmHeightMax/km << " km is: " << intersections.size()
               << " ---> skip event\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
          ERROR(warn);
          return eContinueLoop;
        }

        if (upward) { // -> redefine exit point

          atmDistanceMax = (intersections[0] - showerCore).GetMag();
          atmDepthMax = slantDepthVsDistance.Y(atmDistanceMax);

        } else { // -> redefine entry point

          atmDistanceMin = (intersections[1] - showerCore).GetMag();
          atmDepthMin = slantDepthVsDistance.Y(atmDistanceMin);

          if (skimming) {
            atmDistanceMax = std::min(atmDistanceMax, (intersections[0] - showerCore).GetMag());
            atmDepthMax = slantDepthVsDistance.Y(atmDistanceMax);
            atmHeightMin = heightVsDistance.Y((atmDistanceMax + atmDistanceMin)/2);
          }
        }
      }

      if (atmHeightMin<temperaturProfile.MinX() ||
          atmHeightMin<densityProfile.MinX()) {

        atmHeightMin = std::max(temperaturProfile.MinX(), densityProfile.MinX()) + 1.*mm; // security margin
        const vector<Point> intersections =
          Intersection(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84), atmHeightMin,
                       Line(showerCore, -showerDir));

        if (intersections.size() != 2) {
          ostringstream warn;
          warn << "\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
                  "Number of intersections point of air shower axis with ellipsoid of height="
               << atmHeightMin/km << " km is: " << intersections.size()
               << " ---> skip event\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
          ERROR(warn);
          return eContinueLoop;
        }

        if (upward) { // -> redefine entry point

          atmDistanceMin = (intersections[0] - showerCore).GetMag();
          atmDepthMin = slantDepthVsDistance.Y(atmDistanceMin);

        } else { // -> redefine exit point

          atmDistanceMax = (intersections[1] - showerCore).GetMag();
          atmDepthMax = slantDepthVsDistance.Y(atmDistanceMax);

          if (skimming) {
            atmDistanceMin = std::max(atmDistanceMin, (intersections[0] - showerCore).GetMag());
            atmDepthMin = slantDepthVsDistance.Y(atmDistanceMin);
            atmHeightMax = heightVsDistance.Y((atmDistanceMax + atmDistanceMin)/2);
          }
        }
      }
    } catch (InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {

      ERROR(string("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
                   "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
                   "InclinedAtmosphericProfile EXCEPTION:\n") + e.GetMessage() +
                  //                   " ---> switch back to FLAT GEOMETRY\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
                  "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
      //fFlatEarth = true;     // -> GO ON ANYWAY
      return eContinueLoop;    // -> SKIP EVENT
    }
  }

  if (fFlatEarth) {

    upward = (cosTheta < 0);
    // skimming not suppoted for flat geometry !!

    // check against limits of flat atmosphere
    const ProfileResult& temperaturProfile = theDetector.GetAtmosphere().EvaluateTemperatureVsHeight();
    const ProfileResult& densityProfile = theDetector.GetAtmosphere().EvaluateDensityVsHeight();
    const ProfileResult& depthVsHeight = theDetector.GetAtmosphere().EvaluateDepthVsHeight();

    // get the LIMITS of the atmosphere
    atmDistanceMin = (coreZ - depthVsHeight.MinX()) / absCosTheta;
    atmDistanceMin = std::max(atmDistanceMin, (coreZ - densityProfile.MinX()) / absCosTheta);
    atmDistanceMin = std::max(atmDistanceMin, (coreZ - temperaturProfile.MinX()) / absCosTheta);

    atmDistanceMax = (coreZ - depthVsHeight.MaxX()) / absCosTheta;
    atmDistanceMax = std::min(atmDistanceMax, (coreZ - densityProfile.MaxX()) / absCosTheta);
    atmDistanceMax = std::min(atmDistanceMax, (coreZ - temperaturProfile.MaxX()) / absCosTheta);

    atmDepthMin = depthVsHeight.Y(coreZ - atmDistanceMax*absCosTheta);
    atmDepthMax = depthVsHeight.Y(coreZ - atmDistanceMin*absCosTheta);

    atmDepthMin /= absCosTheta;
    atmDepthMax /= absCosTheta;

    ostringstream info;
    info << "Flat geometry with Theta=" << zenith/deg;
    INFO(info);
  }

  // print some information for the logfile
  ostringstream info;
  info << "Atmosphere starts at " << atmDepthMin/g*cm*cm << "g/cm2, "
          "d=" << atmDistanceMin/m << "m, "
          "and ends at " << atmDepthMax/g*cm*cm << "g/cm2, "
          "d=" << atmDistanceMax/m << "m";
  INFO(info);

  // RESTRICT SHOWER TO VALID ATMOSPHERE !!!!!!!!!!!!!!!!!!!!!!! ------------
  bool truncated = false;

  if (!simShower.HasdEdX()) {
    ERROR("The shower has no energy deposit profile.");
    return eFailure;
  }
  const TabulatedFunction& energyDep = simShower.GetdEdX();
  if (!energyDep.GetNPoints()) {
    ERROR("The shower has an empty energy deposit profile.");
    return eFailure;
  }
  int iFirstValidBin_dEdX = 0;
  int iLastValidBin_dEdX = energyDep.GetNPoints() - 1;

  if (!simShower.HasLongitudinalProfile()) {
    ERROR("The shower has no charge profile.");
    return eFailure;
  }
  const TabulatedFunction& chargeProfile = simShower.GetLongitudinalProfile();
  if (!chargeProfile.GetNPoints()) {
    ERROR("The shower has an empty charge profile.");
    return eFailure;
  }
  int iFirstValidBin_Ne = 0;
  int iLastValidBin_Ne = chargeProfile.GetNPoints() - 1;

  const double profDepthMin =
    std::max(energyDep[iFirstValidBin_dEdX].X(), chargeProfile[iFirstValidBin_Ne].X());
  const double profDepthMax =
    std::min(energyDep[iLastValidBin_dEdX].X(), chargeProfile[iLastValidBin_Ne].X());

  info.str("");
  info << "Shower profile starts at " << profDepthMin/g*cm*cm << " g/cm2, "
          "and ends at " << profDepthMax/g*cm*cm << " g/cm2";
  INFO(info);

  // check dEdX profile
  while (energyDep[iFirstValidBin_dEdX].X() < atmDepthMin &&
         iFirstValidBin_dEdX<iLastValidBin_dEdX) {
    ++iFirstValidBin_dEdX;
    truncated = true;
  }

  while (energyDep[iLastValidBin_dEdX].X() > atmDepthMax &&
         iLastValidBin_dEdX>iFirstValidBin_dEdX) {
    --iLastValidBin_dEdX;
    truncated = true;
  }

  // check Ne profile
  while (chargeProfile[iFirstValidBin_Ne].X() < atmDepthMin &&
         iFirstValidBin_Ne<iLastValidBin_Ne) {
    ++iFirstValidBin_Ne;
    truncated = true;
  }

  while (chargeProfile[iLastValidBin_Ne].X() > atmDepthMax &&
         iLastValidBin_Ne>iFirstValidBin_Ne) {
    --iLastValidBin_Ne;
    truncated = true;
  }

  const double truncDepthMin =
    std::max(energyDep[iFirstValidBin_dEdX].X(), chargeProfile[iFirstValidBin_Ne].X());

  const double truncDepthMax =
    std::min(energyDep[iLastValidBin_dEdX].X(), chargeProfile[iLastValidBin_Ne].X());

  //
  double truncHeightMax;
  double truncHeightMin;
  if (fFlatEarth) {
    truncHeightMax = heightVsDepth.Y(truncDepthMin*absCosTheta);
    truncHeightMin = heightVsDepth.Y(truncDepthMax*absCosTheta);
  } else {
    const ProfileResult& heightVsSlantDepth =
      theDetector.GetAtmosphere().EvaluateHeightVsSlantDepth();
    truncHeightMax = heightVsSlantDepth.Y(truncDepthMin);
    truncHeightMin = heightVsSlantDepth.Y(truncDepthMax);
  }
  // -----------------------------------------------------------------------

  // ------- determine start and stop of observation along shower axis --------
  double distInitCore;
  double distFinalCore;
  if (fFlatEarth) {
    distInitCore = (coreZ-truncHeightMax) / cosTheta;
    distFinalCore = (coreZ-truncHeightMin) / cosTheta;
  } else {
    distInitCore = distanceVsSlantDepth->Y(truncDepthMin);
    distFinalCore = distanceVsSlantDepth->Y(truncDepthMax);
  }
  // -------------------------------------------------------------------------

  /*
    for DOWNWARD shower:
    distance from core to init is smaller than to final point of atmosphere

    for UPWARD shower:
    distance from core to init is bigger than to final point of atmosphere


    DOWNWARD GOING        UPWARD GOING
                      /                  /
       init -->     \/                  ^
                    /                  /
                   /                  /
                  /                  /
       final --> /                  /
    -------------X-----------------X----
                CORE              CORE

   */

  //const bool upwardGoing = (distInitCore>distFinalCore);
  //if (upwardGoing) swap(distFinalCore, distInitCore);

  if (truncated) {
    info.str("");
    info << "Profile was truncated to fit the Atmosphere:\n"
            "            -> start at "
            "X=" << truncDepthMin/g*cm2 << "g/cm2, "
            "h=" << truncHeightMax/m << "m, "
            "d=" << distInitCore/m << "m\n"
            "            -> end at "
            "X=" << truncDepthMax/g*cm2 << "g/cm2, "
            "h=" << truncHeightMin/m << "m, "
            "d=" << distFinalCore/m << "m";
    INFO(info);
  } else {
    info.str("");
    info << "Profile fits into the Atmosphere definition:\n"
            "            -> start at "
            "h=" << truncHeightMax/m << "m, "
            "d=" << distInitCore/m << "m\n"
            "            -> end at "
            "h=" << truncHeightMin/m << "m, "
            "d=" << distFinalCore/m << "m";
    INFO(info);
  }

  if (upward) {
    INFO ("shower is UPWARD going");
  }

  if (fFluorescence) {
    // Calculate Fluorescence light
    FluorescenceLight(event, distInitCore, distFinalCore, upward);
  }

  if (fCherenkov) {

    if (!simShower.HasGHParameters()) {
      ERROR("Shower without GH-fit: Cannot compute Cherenkov light ! SKIPPING !");
      // return eFailure;
      return eContinueLoop;
    }

    // Calculate Cherenkov light beam
    CherenkovLight(event, distInitCore, distFinalCore, upward, fUseCherenkovProfile);
  }

  if (fPrimaryCherenkov) {
    // Calculate Cherenkov light of primary
    PrimaryCherenkovLight(event);
  }

  return eSuccess;
}


VModule::ResultFlag
ShowerLightSimulator::Finish()
{
  return eSuccess;

}


void
ShowerLightSimulator::FluorescenceLight(evt::Event& event,
                                        double distInitCore,
                                        double distFinalCore,
                                        bool /*upwardGoing*/)
{
  INFO ("Calculating fluorescence light");

  Detector& theDetector = Detector::GetInstance();
  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

  evt::ShowerSimData& simShower = event.GetSimShower();
  CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  Point core(0, 0, 0, showerCS);
  Vector axis(0, 0, -1, showerCS);

  // just needed for flat geometry
  const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
  const double cosTheta = (-simShower.GetDirection()).GetCosTheta(localCS);
  const double absCosTheta = std::fabs(cosTheta);
  const double coreZ = wgs84.PointToLatitudeLongitudeHeight(core).get<2>();

  const ProfileResult& slantDepthVsDistance = // for curved geometry
    theDetector.GetAtmosphere().EvaluateSlantDepthVsDistance();
  const ProfileResult& depthVsHeight = // for flat geometry
    theDetector.GetAtmosphere().EvaluateDepthVsHeight();

  const double distInitFinal = std::fabs(distFinalCore - distInitCore);
  const double delta = kSpeedOfLight*fTimeBin;
  const int nPoints = int(distInitFinal / delta);

  const TabulatedFunction& energyDep = simShower.GetdEdX();

  const double dEdX0 = theDetector.GetAtmosphere().GetdEdX0();

  // Vector of wavelengths
  const vector<double>& wavelengths =
    theDetector.GetAtmosphere().GetWavelengths(Atmosphere::eFluorescence);
  unsigned int nWaves = wavelengths.size();

  for (unsigned int iwl = 0; iwl < nWaves; ++iwl) {
    if (!simShower.HasFluorescencePhotons(iwl)) {
      TabulatedFunction fEmpty;
      simShower.AddFluorescencePhotons(fEmpty, iwl);
    } else {
      ERROR("Flourescence light produced twice !!!!!");
    }
  }

  //const double distanceStart = ( upwardGoing ? distFinalCore : distInitCore );
  const double distanceStart = distInitCore;

  for (int i = 0; i < nPoints; ++i) {

    const double distance = distanceStart + delta * (0.5 + i); // take bin middle because of Cherenkov ...
    //const double distance = distanceStart + delta*i;
    //double distance = distFinalCore - delta * (.5+i);

    double heightBin;
    double depthBin;

    if (fFlatEarth) {
      heightBin = coreZ - distance * cosTheta;
      depthBin = depthVsHeight.Y(heightBin) / absCosTheta;
    } else {
      Point p (core + axis * distance);
      heightBin = wgs84.PointToLatitudeLongitudeHeight(p).get<2>();
      depthBin = slantDepthVsDistance.Y(distance);
    }

    const TabulatedFunction& spectrum =
      theDetector.GetAtmosphere().EvaluateFluorescenceYield(heightBin);

    // double cPh = 0;

    for (unsigned int iwl = 0; iwl < nWaves; ++iwl) {

      if (!simShower.HasFluorescencePhotons(iwl)) {
        ERROR ("Data structure SimShowerData not initialized properly!");
        continue;
      }

      const double photons = spectrum[iwl].Y()*energyDep.Y(depthBin)/dEdX0; // [1/m]

      TabulatedFunction& fluoPhotons = simShower.GetFluorescencePhotons(iwl);
      fluoPhotons.Insert(distance, photons);

      //cPh += photons;

    } // end of loop over wavelengths

  } // end of loop over distance
}


void
ShowerLightSimulator::CherenkovLight(evt::Event& event,
                                     double distInitCore,
                                     double distFinalCore,
                                     bool /*upwardGoing*/,
                                     bool useCherenkovProfile)
{
  INFO("Calculating Cherenkov light and beam");

  Detector& theDetector = Detector::GetInstance();
  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

  evt::ShowerSimData& simShower = event.GetSimShower();
  CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  Point core(0, 0, 0, showerCS);
  Vector axis(0, 0, -1, showerCS);

  // just needed for flat geometry
  const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
  const double cosTheta = (-simShower.GetDirection()).GetCosTheta(localCS);
  const double absCosTheta = std::fabs(cosTheta);
  const double coreZ = wgs84.PointToLatitudeLongitudeHeight (core).get<2>();

  const ProfileResult& slantDepthVsDistance = // for curved geometry
    theDetector.GetAtmosphere().EvaluateSlantDepthVsDistance();
  const ProfileResult& depthVsHeight = // for flat geometry
    theDetector.GetAtmosphere().EvaluateDepthVsHeight();

  const double distInitFinal = std::fabs(distFinalCore - distInitCore);
  const double delta = kSpeedOfLight * fTimeBin;
  const int nPoints = int(distInitFinal / delta);

  /* set the Cherenkov energy cutoff for electrons */

  const TabulatedFunction* chargeProfile = nullptr;
  const TabulatedFunction* cherenkovProfile = nullptr;
  double Xmax = 0;
  double wl_spect_norm_prod = 0;

  if (!useCherenkovProfile) {
    const double eCut = simShower.GetEnergyCutoff(evt::ShowerSimData::eElectron);
    theDetector.GetAtmosphere().SetCherenkovEnergyCutoff(eCut);
    chargeProfile = &simShower.GetLongitudinalProfile();
    Xmax = simShower.GetGHParameters().GetXMax();
  } else {
    cherenkovProfile = &simShower.GetLongitudinalProfile(ShowerSimData::eCherenkov);
    const double lambdaMin = simShower.GetMinCherenkovWavelength();
    const double lambdaMax = simShower.GetMaxCherenkovWavelength();
    wl_spect_norm_prod = 1./lambdaMin - 1./lambdaMax;
  }

  // Vector of wavelengths
  const vector<double>& wavelengths =
    theDetector.GetAtmosphere().GetWavelengths(Atmosphere::eCerenkov);
  const unsigned int nWaves = wavelengths.size();

  // initialize data structures
  for (unsigned int iwl = 0; iwl < nWaves; ++iwl) {

    if (!simShower.HasCherenkovPhotons(iwl)) {
      TabulatedFunction fTmp;
      simShower.AddCherenkovPhotons(fTmp, iwl);
    } else {
      ERROR("Cherenkov light produced twice !!!!!");
    }

    if (!simShower.HasCherenkovBeamPhotons(iwl)) {
      TabulatedFunction fTmp;
      simShower.AddCherenkovBeamPhotons(fTmp, iwl);
    } else {
      ERROR("Cherenkov beam produced twice !!!!!");
    }

    if (!simShower.HasCherenkovBeamProductionPhotons(iwl)) {
      TabulatedFunction fTmp;
      simShower.AddCherenkovBeamProductionPhotons(fTmp, iwl);
    } else {
      ERROR("Cherenkov beam produced twice !!!!!");
    }
  }

  //const double distanceStart = ( upwardGoing ? distFinalCore : distInitCore );
  const double distanceStart = distInitCore;

  vector<double> beam(nWaves, 0); // cherenkov beam photon number
  for (int i = 0; i < nPoints; ++i) {

    const double distance = distanceStart + delta * (0.5 + i); // use middle of bin
    const double distanceStart = distance - 0.5*delta;
    const double distanceEnd = distance + 0.5*delta;
    const Point xStepStart(core + axis * distanceStart);
    const Point xStepEnd(core + axis * distanceEnd);

    double heightBin = 0;
    double depthBin = 0;
    if (fFlatEarth) {
      heightBin = coreZ - distance * cosTheta;
      depthBin = depthVsHeight.Y(heightBin) / absCosTheta;
    } else {
      const Point p(core + axis * distance);
      heightBin = wgs84.PointToLatitudeLongitudeHeight(p).get<2>();
      depthBin = slantDepthVsDistance.Y(distance);
    }

    AttenuationResult mieAtt =
      theDetector.GetAtmosphere().EvaluateMieAttenuation(xStepStart, xStepEnd, wavelengths);
    AttenuationResult rayleightAtt =
      theDetector.GetAtmosphere().EvaluateRayleighAttenuation(xStepStart, xStepEnd, wavelengths);

    double showerAge = 0;
    const TabulatedFunction* cherPh = nullptr;
    double cherenkovProduction = 0;
    if (!useCherenkovProfile) {
      showerAge = utl::ShowerAge(depthBin,Xmax);
      cherPh = &theDetector.GetAtmosphere().EvaluateCherenkovPhotons (xStepStart, xStepEnd, showerAge);
    } else {
      cherenkovProduction = cherenkovProfile->Y(depthBin) / wl_spect_norm_prod;
      // This is cherenkov production PER _tracklength_  !
    }

    // We calculate three quantities here, all as a (binned) function
    // along the shower axis and separately for each wavelength:
    // the Cherenkov photons
    // => Simply the number of Cherenkov photons that are generated
    //    in the given distance bin per distance.
    // the CherenkovBeamPhotons
    // => These are the accumulated photons in the beam, attenuated
    //    along the track:
    //    beam[i] = (beam[i-1] + generated[i]) * attenuation[i]
    // the CherenkovBeamProductionPhotons
    // => The tricky bit. For a total of n bins:
    //    beamP[i] = generated[i] * Product[j=i..n](attenuation[j])
    //    That means the last bin simply contains generated[n]*attenuation[n],
    //    the second-to-last contains
    //      generated[n-1]*attenuation[n-1]*attenuation[n]
    //    and so on.
    //    This makes sense because in the ShowerPhotonGenerator, we walk
    //    along the axis again. When generating photons there, we draw their
    //    source distance from a PDF that is taken to be the above beamP[i]
    //    for all bins up to and including the current one.
    //    When looking at this PDF for any bin k except the last, it will
    //    have a constant factor Product[i=k+1..n](attenuation[i]), but
    //    since the PDF is automatically normalized, this factor (that is the
    //    same for all bins) drops out.
    //    Using this strategy, we can calculate these PDFs for all bins in
    //    one go.

    double minWaveBin = 0;
    double maxWaveBin = 0;
    for (unsigned int iwl = 0; iwl < nWaves; ++iwl) {

      const double mieAttenuation = mieAtt.GetTransmissionFactor()[iwl].Y();
      const double rayleightAttenuation = rayleightAtt.GetTransmissionFactor()[iwl].Y();
      const double attenuation = mieAttenuation * rayleightAttenuation;

      // -------------------------
      // WARNING:
      // THIS is very bad code, see also AnalyticalCherenkovModel:
      //      the "binning" is not defined. The user must know how this works...
      minWaveBin = (iwl == 0 ?
                    wavelengths[iwl] - 0.5*(wavelengths[iwl+1] - wavelengths[iwl]) :
                    maxWaveBin);
      maxWaveBin = (iwl == nWaves-1 ?
                    wavelengths[iwl] + 0.5*(wavelengths[iwl] - wavelengths[iwl-1]) :
                    0.5*(wavelengths[iwl] + wavelengths[iwl+1]));
      // -------------------------

      // Apply attenuation to entire beam, not just the
      // most recent contribution.

      double generatedPhotons = 0;
      if (!useCherenkovProfile) {
        const double Nch = chargeProfile->Y(depthBin);
        generatedPhotons = cherPh->Y(wavelengths[iwl]) * Nch;
        beam[iwl] += generatedPhotons;

      } else {

        // wavelength spectrum, dN/dlambda = 1/lambda2
        const double deltaDepth = fFlatEarth ?
          (depthVsHeight.Y(coreZ - distanceEnd*cosTheta) - depthVsHeight.Y(coreZ - distanceStart*cosTheta)) / absCosTheta :
          slantDepthVsDistance.Y(distanceEnd) - slantDepthVsDistance.Y(distanceStart);
        generatedPhotons = cherenkovProduction * deltaDepth * (1./minWaveBin - 1./maxWaveBin);
        beam[iwl] += generatedPhotons;
      }

      beam[iwl] *= attenuation;

      TabulatedFunction& ckovPhotons = simShower.GetCherenkovPhotons(iwl);
      TabulatedFunction& ckovBeam = simShower.GetCherenkovBeamPhotons(iwl);
      TabulatedFunction& ckovBeamProduction = simShower.GetCherenkovBeamProductionPhotons(iwl);

      ckovPhotons.Insert(distance, generatedPhotons/delta);
      ckovBeam.Insert(distance, beam[iwl]);
      ckovBeamProduction.Insert(distance, generatedPhotons/delta);

      for (unsigned int jTrace = 0; jTrace < ckovBeamProduction.GetNPoints(); ++jTrace) {
        ckovBeamProduction.GetY(jTrace) *= attenuation;
        if (ckovBeamProduction.GetY(jTrace) < 0) {
          WARNING("Cherenkov beam production negative. Forcing it to zero.");
          ckovBeamProduction.GetY(jTrace) = 0.;
        }
      }

    } // loop over wavelengths

  } // loop along shower axis
}


void
ShowerLightSimulator::PrimaryCherenkovLight(evt::Event& event)
{
  INFO("Calculating primary Cherenkov light and beam");

  Detector& theDetector = Detector::GetInstance();
  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

  evt::ShowerSimData& simShower = event.GetSimShower();
  CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  Point core(0, 0, 0, showerCS);
  Vector axis(0, 0, -1, showerCS);

  //const double cosTheta = cos (simShower.GetZenith()); // just needed for flat geometry

  const ProfileResult& depthVsHeight = // for flat geometry
    theDetector.GetAtmosphere().EvaluateDepthVsHeight();
  const ProfileResult& heightVsSlantDepth =
    theDetector.GetAtmosphere().EvaluateHeightVsSlantDepth();
  const ProfileResult& distanceVsSlantDepth =
    theDetector.GetAtmosphere().EvaluateDistanceVsSlantDepth();

  //const double distInitFinal = std::fabs(distFinalCore-distInitCore);
  //const double delta = kSpeedOfLight*(fTimeBin);
  //const int nPoints = int(distInitFinal / delta);

  // get the LIMITS of the atmosphere
  double atmDepthMin = heightVsSlantDepth.MinX();
  double atmDepthMax = heightVsSlantDepth.MaxX();
  double atmHeightMin = heightVsSlantDepth.Y (atmDepthMax);

  double Z = 0;
  if (fOverRidePrimaryCharge >= 0) {
    Z = fOverRidePrimaryCharge;
    cout << "  Override primary charge ----------> Z=" << Z << endl;
  } else {
    const int particleID = simShower.GetPrimaryParticle();
    switch (particleID) {
    case utl::Particle::ePhoton:
      Z = 0;
      break;
    case utl::Particle::eNeutron:
      Z = 0;
      break;
    case utl::Particle::eProton:
      Z = 1;
      break;
    case utl::Particle::eIron:
      Z = 26;
      break;
    default:
      ostringstream err;
      err << "Unkown primary particle type " << particleID;
      ERROR(err);
      throw AugerException(err.str());
      break;
    }
  }

  // just needed for flat geometry
  //double cosTheta = cos (simShower.GetZenith());
  //double coreZ = wgs84.PointToLatitudeLongitudeHeight (core).get<2>();

  const double Xfirst = simShower.GetXFirst();

  if (Xfirst < atmDepthMin) {
    ostringstream info;
    info << "Xfirst is out of the bounds of the atmosphere. X1=" << Xfirst/g*cm2
         << ", Xmin=" << atmDepthMin/g*cm2 ;
    INFO(info);
    return;
  }

  double height1 = heightVsSlantDepth.Y(atmDepthMin);
  double height2 = heightVsSlantDepth.Y(Xfirst);
  double dist1 = distanceVsSlantDepth.Y(atmDepthMin);
  double dist2 = distanceVsSlantDepth.Y(Xfirst);
  double distInitFinal = std::abs(dist2 - dist1);
  double delta = kSpeedOfLight * fTimeBin;
  int nPoints = std::min(10000, int(distInitFinal / delta));
  dist1 = dist2 - nPoints*delta;                  // correction to right start point

  ostringstream info;
  info << "track primary from depth="
       << atmDepthMin/g*cm2 << " g/cm2, dist=" << dist1/m << " m, height=" << height1/m << " m; "
          "point of first interaction at "
       << Xfirst/(g/cm2) << " " << dist2/m << " " << height2/m;
  INFO(info);

  // use cherenkov wl bins for now !!!! TODO check
  vector<double> wavelength = theDetector.GetAtmosphere().GetWavelengths(Atmosphere::eCerenkov);
  const unsigned int nWaves = wavelength.size();

  const ProfileResult& refractiveIndex = theDetector.GetAtmosphere().EvaluateRefractionIndexVsHeight();
  const double minHeight = refractiveIndex.MinX();
  const double maxHeight = refractiveIndex.MaxX();
  const double minRI = refractiveIndex.Y(maxHeight);
  //const double maxRI = refractiveIndex.Y(minHeight);

  const double minRIdepth = depthVsHeight.Y(maxHeight);

  double maxWaveBin = wavelength[0] - (wavelength[1] - wavelength[0])/2;
  for (unsigned int wl = 0; wl < nWaves; ++wl) {

    // wavelength bins (centered around mean value)
    double minWaveBin = maxWaveBin;

    if (wl == nWaves-1)
      maxWaveBin = wavelength[wl] + (wavelength[wl] - wavelength[wl-1])/2.;
    else
      maxWaveBin = (wavelength[wl] + wavelength[wl+1])/2.;

    double lambda1 = minWaveBin;
    double lambda2 = maxWaveBin;

    vector<double> vPhotons;
    vector<double> vPhotonsBeam;
    vector<double> vDist;

    double beam = 0; // beam photons
    for (int i = 0; i < nPoints*1.5; ++i) {

      double distBin1 = dist1 + delta * i;
      double distBin2 = dist1 + delta * (i + 1);
      double dist = (distBin2 + distBin1)/2.;

      Point p1 = core + distBin1*axis;
      Point p2 = core + distBin2*axis;
      Point p = core + dist*axis;
      double height = wgs84.PointToLatitudeLongitudeHeight (p).get<2>();

      if (height < atmHeightMin)
        break; // left atmosphere

      vDist.push_back(dist);

      if (i <= 0) {
        vPhotons.push_back(0);
        vPhotonsBeam.push_back(0);
        continue;
      }

      AttenuationResult mieAtt =
        theDetector.GetAtmosphere().EvaluateMieAttenuation(p1, p2, wavelength);
      AttenuationResult rayleightAtt =
        theDetector.GetAtmosphere().EvaluateRayleighAttenuation(p1, p2, wavelength);

      double mie = mieAtt.GetTransmissionFactor()[wl].Y();
      double rayleight = rayleightAtt.GetTransmissionFactor()[wl].Y();
      double attenuation = mie * rayleight;

      // refractive index scales with air density (approximation taken from F. Nerlings thesis)
      double vertDepth = depthVsHeight.Y(height);
      double n_height = (minRI - 0.000283) + 0.000283 * (vertDepth/minRIdepth);
      if (height > minHeight && height < maxHeight) {
        n_height = refractiveIndex.Y(height);
      }
      // wavelength dependence (approximation taken from F. Nerlings thesis)
      double lambda = (lambda1 + lambda2)/2 / (1e-10*m); // wavelength in [angstrom]
      double n = n_height + 1e-7 * (2726.42 + 12.288/(lambda*lambda*1e-8) + 0.3555/(lambda*lambda*lambda*lambda+1.e-16) );

      double beta = 1.;  // CHECK, if this is good enough

      double dNdl = Z*Z * 2.*kPi * 1./137. * (1.-1./(pow(n*beta, 2))) * (-1) * (1./lambda2 - 1./lambda1);

      if (i >= nPoints-1) { // after first interaction
        dNdl = 0;
      }

      // dN/dl = Z^2 2.*kPi*alpha*(1-1/(n*n * beta*beta)) * dlambda 1/lambda2;
      // cos (delta) = 1/(beta n)
      // beta_thresh = 1/(n)

      //double dN = dNdl * delta;
      //double dN_BinBefore = vPhotons.back() * delta;
      beam += dNdl * delta;
      beam *= attenuation;

      vPhotons.push_back(dNdl);
      vPhotonsBeam.push_back(beam);

    } // loop along axis

    TabulatedFunction cher(vDist, vPhotons);
    TabulatedFunction cherBeam(vDist, vPhotonsBeam);

#ifdef PRIM_CHER
    simShower.AddPrimaryCherenkovPhotons(cher, wl);
    simShower.AddPrimaryCherenkovBeamPhotons(cherBeam, wl);
#else
    ERROR("You cannot us primary particle cherenkov light, without compiling with -DPRIM_CHER !!! ");
    throw AugerException();
#endif

  } // loop wavelengths

}