LightAtDiaphragmSimulator.cc

Go to the documentation of this file.

#include <iomanip>
#include <algorithm>
#include <vector>
#include <locale>

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/TransformationMatrix.h>
#include <utl/TimeStamp.h>
#include <utl/AugerCoordinateSystem.h>

#include <utl/Particle.h>

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

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/Telescope.h>

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

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

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

#include "LightAtDiaphragmSimulator.h"

#include <TH2D.h> // for debug
#include <TH1D.h> // for debug
#include <TCanvas.h>// for debug
#include <TGraph.h>// for debug
#include <TFile.h>// for debug

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


VModule::ResultFlag
LightAtDiaphragmSimulator::Init()
{
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("LightAtDiaphragmSimulatorKG");

  topB.GetChild("fluorDirect").GetData(fFluorDirect);

  topB.GetChild("cherDirect").GetData(fCherDirect);
  topB.GetChild("cherDirectCORSIKA").GetData(fCherDirectCORSIKA);
  topB.GetChild("cherMieScattered").GetData(fCherMieScattered);
  topB.GetChild("cherRayleighScattered").GetData(fCherRayleighScattered);

  topB.GetChild("laserMieScattered").GetData(fLaserMieScattered);
  topB.GetChild("laserRayleighScattered").GetData(fLaserRayleighScattered);
  topB.GetChild("wavelengthDepRefraction").GetData(fWlRefrac);

  const bool showerMode = fFluorDirect ||
    fCherDirect || fCherDirectCORSIKA || fCherRayleighScattered || fCherMieScattered;

  const bool laserMode = fLaserMieScattered || fLaserRayleighScattered;

  if (showerMode && laserMode) {
    ERROR("You cannot propagate shower and laser light simultanously. "
           "Check your xml data card!");
    return eFailure;
  }

  if (fCherDirect && fCherDirectCORSIKA) {
    ERROR("You cannot simulate direct Cherenkov from parameterization and CORSIKA at the same time. Check your xml settings!");
    return eFailure;
  }

  // -------- do some output ------------
  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n";
  if (showerMode) {
    info << "            -- shower mode -- \n"
            "         fluorescence: " << fFluorDirect << "\n"
            "     cherenkov direct: " << fCherDirect << " " << (fCherDirectCORSIKA ? "[CORSIKA]" : "[model]") << "\n"
            "        cherenkov Mie: " << fCherMieScattered << "\n"
            "   cherenkov Rayleigh: " << fCherRayleighScattered << "\n"
            "    wavelength dep. n: " << fWlRefrac << "\n";
  } else if (laserMode) {
    info << "            -- laser mode -- \n"
            "  laser Mie scattered: " << fLaserMieScattered << "\n"
            "  laser Rayleigh scat: " << fLaserRayleighScattered << "\n";
  } else {
    info << "           -- unknown mode !!! -- \n";
  }
  INFO(info);

  return eSuccess;
}


VModule::ResultFlag
LightAtDiaphragmSimulator::Run(evt::Event& event)
{
  if (!event.HasFEvent()) {
    ERROR("Missing FEvent. Check your module configurations and ModuleSequence.");
    return eFailure;
  }

  if (!event.HasSimShower()) {
    ERROR("Event has no simulated shower.");
    return eFailure;
  }

  evt::ShowerSimData& simShower = event.GetSimShower();
  fIsLaserEvent = simShower.HasLaserData();
  const TimeStamp timeAtCore = simShower.GetTimeStamp();

  if (!simShower.HasLongitudinalProfile() && !fIsLaserEvent) {
    ERROR("Event has no longitudinal profile for the simulated shower.");
    return eFailure;
  }

  const CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();

  const Point core(0, 0, 0, showerCS);
  const Vector showerAxis(0, 0, -1, showerCS);
  const double zenith = (-simShower.GetDirection()).GetTheta(localCS); // only for flat earth ....
  const double absCosZenith = std::fabs(cos(zenith));

  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();
  const double coreZ = wgs84.PointToLatitudeLongitudeHeight(core).get<2>();

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

  bool flatEarth = false;
  ProfileResult const * slantDepthVsDistance = 0;
  ProfileResult const * depthProfile = 0;
  try {
    slantDepthVsDistance = &atmo.EvaluateSlantDepthVsDistance();
  } catch (InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
    WARNING("using flat earth!");
    flatEarth = true;
    depthProfile = &atmo.EvaluateDepthVsHeight();
  }

  // get shower/laser specific information
  double Xmax = 0;
  vector <double> WLengthLaser;
  if (fIsLaserEvent) {
    double laserWavelength =  event.GetSimShower().GetLaserData().GetLaserWavelength();
    WLengthLaser.push_back (laserWavelength);
  } else {
    if (!simShower.HasGHParameters()) {
      ERROR("Shower without GH-fit: Cannot compute 3D shower structure ! SKIPPING !");
      // return eFailure;
      return eContinueLoop;
    }
    Xmax = simShower.GetGHParameters().GetXMax();
  }

  const vector<double>& WLengthFluo = atmo.GetWavelengths(Atmosphere::eFluorescence);
  const vector<double>& WLengthCkov = atmo.GetWavelengths(Atmosphere::eCerenkov);
  const unsigned int NWLengthFluo = WLengthFluo.size();
  const unsigned int NWLengthCkov = WLengthCkov.size();
  //const unsigned int nWaveLengths = (fIsLaserEvent ? 1 : NWLengthCkov+NWLengthFluo);

  /*
    ----------------------------------------------------------------------------------
    Now loop over all diaphragm time traces, and fill them
  */

  INFO(" ... computing, this may take a while ... ");

  const FDetector& detFD = detector.GetFDetector();
  fevt::FEvent& fevent = event.GetFEvent();

  for (fevt::FEvent::EyeIterator iEye = fevent.EyesBegin(fevt::ComponentSelector::eInDAQ);
       iEye != fevent.EyesEnd(fevt::ComponentSelector::eInDAQ) ; ++iEye) {

    fevt::Eye& eyeEvent = *iEye;

    if (eyeEvent.GetStatus() == fevt::ComponentSelector::eDeSelected)
      continue;

    for (fevt::Eye::TelescopeIterator iTel = eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
         iTel != eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++iTel) {

      fevt::Telescope& telEvent = *iTel;
      const fdet::Telescope& detTel = detFD.GetTelescope(telEvent);
      const Point& telPosition = detTel.GetPosition();

      if (!telEvent.HasSimData())
        continue;
      fevt::TelescopeSimData& telSim = telEvent.GetSimData();

      if (!telSim.HasDistanceTrace())
        continue;
      TraceD& distanceTrace = telSim.GetDistanceTrace();

      const double tracebin = distanceTrace.GetBinning();
      const double nBins = distanceTrace.GetSize();

      if (nBins<2)
        continue;

      const TimeStamp photonStartTimeAtDia = telSim.GetPhotonsStartTime();
      const double tMinAtDia = photonStartTimeAtDia - timeAtCore;

      const Vector telToCore = telPosition - core;
      const double distTelToCore = telToCore.GetMag();
      const double cosBeta = CosAngle(telToCore, showerAxis);
      const double Rp = cross(showerAxis, telToCore).GetMag() / showerAxis.GetMag();
      const double T0 = distTelToCore * cosBeta / kSpeedOfLight;

      // USE: dist0 = T0 * kSpeedOfLight
      // const double dist0 = std::fabs(distStartCore) + distTelToCore * cosBeta;
      // OLD CALCULATION: sqrt (pow ((showerInitialPos-telPosition).GetMag(), 2) - Rp*Rp);

      // loop over all bins of trace at diaphragm
      for (int iBin = 0; iBin < nBins; ++iBin) {

        const double timeBinStart = tMinAtDia + tracebin * iBin; // at diaphragm, relative to photonStartTime
        const double timeBin      = timeBinStart + 0.5*tracebin;
        const double timeBinEnd   = timeBinStart + tracebin;

        //const double delta = CalculatePointOnAxis(timeBin, Rp, dist0);
        const double distance = CalculateDistanceFromCore(timeBin, Rp, T0);
        const Point pave = core + showerAxis * distance;

        // cout << " timeBinStart=" << timeBinStart << " " << ((telPosition-pave).GetMag()-(pave-core).GetMag())/kSpeedOfLight << " beta=" << acos(cosBeta)/deg << endl;

        //const double deltaStart = CalculatePointOnAxis(timeBinStart, Rp, dist0);
        const double distanceStart = CalculateDistanceFromCore(timeBinStart, Rp, T0);
        const Point stepStart = core + showerAxis * distanceStart;

        //const double deltaEnd = CalculatePointOnAxis(timeBinEnd, Rp, dist0);
        const double distanceEnd = CalculateDistanceFromCore(timeBinEnd, Rp, T0);
        const Point stepEnd = core + showerAxis * distanceEnd;

        const double xStepWidth = std::fabs(distanceEnd-distanceStart);
        distanceTrace[iBin] = distance;

        if (slantDepthVsDistance &&
            (distance > slantDepthVsDistance->MaxX() ||
             distance < slantDepthVsDistance->MinX())) {
          /*
          ostringstream warn;
          warn << " skipping trace bin #" << iBin << " ["
               << timeBinStart/ns << ", " << timeBinEnd/ns << "] ns, "
               << " d = " << distance/km << " km, slant table = ["
               << slantDepthVsDistance->MinX()/km << ", "
               << slantDepthVsDistance->MaxX()/km << "] km";
          WARNING(warn);
          */
          continue;
        }

        const Vector paveDir = pave - telPosition;
        const double paveDist2 = paveDir.GetMag2();
        const double paveDist = sqrt(paveDist2);
        const double geoFactor = 1. / (4.*kPi*paveDist2);

        //const double diaphragmArea = detTel.GetDiaphragmArea();
        //const double projectedDiaphragmArea = diaphragmArea * paveDir.GetCosTheta(telCS); // unused -- silence warning

        /*
        cout << " light: eye/tel: " << eyeId << "/" << telId
             << " i " << iBin
             << " distance/m=" << distance/m
             << " r=" << paveDist/m
             << " t=" << timeBin/ns
             << " Rp=" << Rp/m
             << " T0=" << T0/ns
             << " tMinAtDia=" << tMinAtDia/ns
             << " Xmax=" << slantDepthVsDistance->MaxX()/m
             << " Xmin=" << slantDepthVsDistance->MinX()/m
             << endl;
        */

        if (!fIsLaserEvent) {

          // FLUORESCENCE PART

          // prepare fluorescence
          // --- read all needed atmospheric properties -------
          const TabulatedFunctionErrors* rAttFluo = 0;
          const TabulatedFunctionErrors* mAttFluo = 0;

          // this is needed HERE because the Atmosphere does not return references ! TODO: change !
          static AttenuationResult rAttenFluo;
          static AttenuationResult mAttenFluo;
          // -----------------------

          if (fFluorDirect) {
            rAttenFluo = atmo.EvaluateRayleighAttenuation(pave, telPosition, WLengthFluo);
            rAttFluo = &rAttenFluo.GetTransmissionFactor();
            mAttenFluo = atmo.EvaluateMieAttenuation(telPosition, pave, WLengthFluo);
            mAttFluo = &mAttenFluo.GetTransmissionFactor();
          }
          // --- properties read end -

          // All the calculations are now wavelength dependent
          for (unsigned int iwl = 0; iwl < NWLengthFluo; ++iwl) {

            if (!simShower.HasFluorescencePhotons(iwl))
              continue;

            const TabulatedFunction& fluoPhotons =
              simShower.GetFluorescencePhotons(iwl);

            // DIRECT FLUORESCENCE ----------------------------------
            if (fFluorDirect) {

              const double att =  (*rAttFluo)[iwl].Y() * (*mAttFluo)[iwl].Y();

              if (!telSim.HasPhotonTrace(fevt::FdConstants::eFluorDirect, iwl))
                telSim.MakePhotonTrace(fevt::FdConstants::eFluorDirect, iwl, nBins, tracebin);
              TraceD& fluorDirectTrace = telSim.GetPhotonTrace(fevt::FdConstants::eFluorDirect, iwl);

              double auxPh = 0;
              if (fluoPhotons.GetNPoints() &&
                  *(fluoPhotons.XBegin()) < distance &&
                  *(fluoPhotons.XEnd() - 1) > distance)
                auxPh = fluoPhotons.Y(distance);
              const double photons = auxPh * geoFactor * att * xStepWidth; // * projectedDiaphragmArea;
              fluorDirectTrace[iBin] += photons;
            } // end of DirectFluorescence

          } // end of fluorescence only wavelength loop

          //cerr << " ph_ax " << photonsAtAxis << " ph_dia " << photonsAtDia;

          // CHERENKOV PART

          if (fCherDirect || fCherMieScattered || fCherRayleighScattered) {

            // prepare cherenkov

            const double Xslant = (flatEarth ?
                                   depthProfile->Y(coreZ-distance*absCosZenith)/absCosZenith :
                                   slantDepthVsDistance->Y(distance));
            const double showerAge = utl::ShowerAge(Xslant, Xmax);

            // --- read all needed atmospheric properties -------

            // Attenuation and geometrical factors
            const TabulatedFunctionErrors* rAttCkov = 0;
            const TabulatedFunctionErrors* mAttCkov = 0;
            const TabulatedFunctionErrors* mFactor = 0;
            const TabulatedFunctionErrors* rFactor = 0;
            double directCher = 0;

            // this is needed because the Atmosphere does not return references ! TODO: change !
            static AttenuationResult rAttenCkov; // static is evil, mkay?
            static AttenuationResult mAttenCkov;
            static ScatteringResult rScat;
            static ScatteringResult mScat;

            rAttenCkov = atmo.EvaluateRayleighAttenuation(pave, telPosition, WLengthCkov);
            mAttenCkov = atmo.EvaluateMieAttenuation(telPosition, pave, WLengthCkov);
            // -----------------------
            rAttCkov = &rAttenCkov.GetTransmissionFactor();
            mAttCkov = &mAttenCkov.GetTransmissionFactor();

            const double scattAngle = acos(-showerAxis*paveDir / paveDir.GetMag());

            // mie scattering factors
            if (fCherMieScattered) {
              mScat = atmo.EvaluateMieScattering(stepStart, stepEnd,
                                                 scattAngle, paveDist,
                                                 WLengthCkov);
              mFactor = &mScat.GetScatteringFactor();
            }

            // rayleight scattering factors
            if (fCherRayleighScattered) {
              rScat = atmo.EvaluateRayleighScattering(stepStart, stepEnd,
                                                      scattAngle, paveDist,
                                                      WLengthCkov);
              rFactor = &rScat.GetScatteringFactor();
            }

            // direct cerenkov
            if (fCherDirect)
              directCher = fWlRefrac ? 1. : atmo.EvaluateDirectCherenkovProbability(stepStart, stepEnd, telPosition,
                                                                   showerAge);
            // --- properties read end -

            // All the calculations are now wavelength dependent
            for (unsigned int iwl = 0; iwl < NWLengthCkov; ++iwl) {

              if (!simShower.HasCherenkovPhotons(iwl))
                continue;

              const double att = (*rAttCkov)[iwl].Y() * (*mAttCkov)[iwl].Y();  // * projectedDiaphragmArea;

              double beamPhotons = 0;
              if (simShower.HasCherenkovBeamPhotons(iwl)) {
                const TabulatedFunction& cherPhotons = simShower.GetCherenkovBeamPhotons(iwl);
                if (cherPhotons.GetNPoints() &&
                    *cherPhotons.XBegin() < distance &&
                    *(cherPhotons.XEnd() - 1) > distance)
                  beamPhotons += cherPhotons.Y(distance);
              }

              double directPhotons = 0;
              if (simShower.HasCherenkovPhotons(iwl)) {
                const TabulatedFunction& cherPhotons = simShower.GetCherenkovPhotons(iwl);
                if (cherPhotons.GetNPoints() &&
                    *cherPhotons.XBegin() < distance &&
                    *(cherPhotons.XEnd() - 1) > distance)
                  directPhotons += cherPhotons.Y (distance);
              }

              // MIE-SCATTERED CHERENKOV LIGHT --------------------------
              if (fCherMieScattered) {
                if (!telSim.HasPhotonTrace(fevt::FdConstants::eCherMieScattered, iwl))
                  telSim.MakePhotonTrace(fevt::FdConstants::eCherMieScattered, iwl, nBins, tracebin);
                TraceD& cherMieTrace = telSim.GetPhotonTrace(fevt::FdConstants::eCherMieScattered, iwl);
                const double photons = beamPhotons * (*mFactor)[iwl].Y() * att;
                cherMieTrace[iBin] += photons;
              }

              // RAYLEIGH-SCATTERED CERENKOV LIGHT -----------------------
              if (fCherRayleighScattered) {
                if (!telSim.HasPhotonTrace(fevt::FdConstants::eCherRayleighScattered, iwl))
                  telSim.MakePhotonTrace(fevt::FdConstants::eCherRayleighScattered, iwl, nBins, tracebin);
                TraceD& cherRayleighTrace = telSim.GetPhotonTrace (fevt::FdConstants::eCherRayleighScattered, iwl);
                const double photons = beamPhotons * (*rFactor)[iwl].Y() *  att;
                cherRayleighTrace[iBin] += photons;
              }

              // DIRECT CERENKOV -----------------------------------------
              if (fCherDirect) {
                const double directCherWl = fWlRefrac ? atmo.EvaluateDirectCherenkovProbability(stepStart, stepEnd, telPosition,
                                                                   showerAge, WLengthCkov[iwl]) : 1.;
                if (!telSim.HasPhotonTrace(fevt::FdConstants::eCherDirect, iwl))
                  telSim.MakePhotonTrace(fevt::FdConstants::eCherDirect, iwl, nBins, tracebin);
                TraceD& cherDirectTrace = telSim.GetPhotonTrace (fevt::FdConstants::eCherDirect, iwl);
                const double photons = directPhotons * directCher*directCherWl * att * xStepWidth;
                cherDirectTrace[iBin] += photons;
              }

            } // end Cerenkov only loop over wavelengths

          } // end of cherenkov

        } else { // if LASER

          // -----------------------------------------------
          // LASER PART

          // prepare laser

          // --- read all needed atmospheric properties -------
          const TabulatedFunctionErrors* rAttLaser;
          const TabulatedFunctionErrors* mAttLaser;

          const double scattAngle = acos(-showerAxis*paveDir / paveDir.GetMag());

          AttenuationResult rAttenLaser = atmo.EvaluateRayleighAttenuation(pave, telPosition,
                                                                           WLengthLaser);
          rAttLaser = &rAttenLaser.GetTransmissionFactor();

          AttenuationResult mAttenLaser = atmo.EvaluateMieAttenuation(telPosition, pave, WLengthLaser);
          mAttLaser = &mAttenLaser.GetTransmissionFactor();

          ScatteringResult rScat = atmo.EvaluateRayleighScattering(stepStart, stepEnd,
                                                                   scattAngle, paveDist,
                                                                   WLengthLaser);
          const TabulatedFunctionErrors& rFactor = rScat.GetScatteringFactor();

          ScatteringResult mScat = atmo.EvaluateMieScattering(stepStart, stepEnd,
                                                              scattAngle, paveDist,
                                                              WLengthLaser);
          const TabulatedFunctionErrors& mFactor = mScat.GetScatteringFactor();

          // there is only one laser wavelength
          double att =  (*rAttLaser)[0].Y() * (*mAttLaser)[0].Y();// * projectedDiaphragmArea;
          double laserPhotons = 0;

          const TabulatedFunction& laserLight = simShower.GetFluorescencePhotons(0);
          if (laserLight.GetNPoints() &&
              *laserLight.XBegin() < distance &&
              *(laserLight.XEnd() - 1) > distance)
            laserPhotons += laserLight.Y(distance);

          if (fLaserMieScattered) {

            // Mie-scattered Laser light simulation
            if (!telSim.HasPhotonTrace(fevt::FdConstants::eLaserMieScattered, 0))
              telSim.MakePhotonTrace(fevt::FdConstants::eLaserMieScattered, 0, nBins, tracebin);
            TraceD& laserMieTrace = telSim.GetPhotonTrace(fevt::FdConstants::eLaserMieScattered, 0);
            const double photons = laserPhotons * att * mFactor[0].Y();
            laserMieTrace[iBin] += photons;
          }

          if (fLaserRayleighScattered) {

            if (!telSim.HasPhotonTrace(fevt::FdConstants::eLaserRayleighScattered, 0))
              telSim.MakePhotonTrace(fevt::FdConstants::eLaserRayleighScattered, 0, nBins, tracebin);
            TraceD& laserRayleighTrace = telSim.GetPhotonTrace(fevt::FdConstants::eLaserRayleighScattered, 0);
            const double photons = laserPhotons * att * rFactor[0].Y();
            laserRayleighTrace[iBin] += photons;
          }

        } // end laser sim

        // cerr << endl;

      } // loop bins

    } // loop tels

  } // end loop eyes

  if (fCherDirectCORSIKA) {
    EvaluateDirectCherenkovHits(event);
  }

  return eSuccess;
}



VModule::ResultFlag
LightAtDiaphragmSimulator::Finish()
{
  return eSuccess;
}


/*
inline double
LightAtDiaphragmSimulator::CalculatePointOnAxis(const double time,
                                                const double Rp,
                                                const double dist0)
{
  // corresponding distance: d = t*c
  const double dist = time * utl::kSpeedOfLight;

  const double delta =
    0.5*(dist0*dist0 + Rp*Rp - dist*dist) / (dist0 - dist);

  return delta;
}
*/


inline double
LightAtDiaphragmSimulator::CalculateDistanceFromCore(const double tDia,
                                                     const double Rp,
                                                     const double T0)
{
  // corresponding distance: d = t*c
  // const double dist = time * utl::kSpeedOfLight;
  //return (pow(T0,2) + pow(Rp/kSpeedOfLight,2) - pow(tDia,2)) / (T0 - tDia) / 2;
  return kSpeedOfLight/2 * (pow(Rp/kSpeedOfLight,2) - pow(tDia,2) + pow(T0,2)) / (T0 - tDia);
}



struct my_facet : public std::numpunct<char> {
    explicit my_facet(size_t refs = 0) : std::numpunct<char>(refs) {}
    virtual char do_thousands_sep() const { return '.'; }
    virtual std::string do_grouping() const { return "\003"; }
};

void
LightAtDiaphragmSimulator::EvaluateDirectCherenkovHits(evt::Event &event)
{
  ShowerSimData& simShower = event.GetSimShower();
  const TimeStamp timeAtCore = simShower.GetTimeStamp();
  Detector& detector = Detector::GetInstance();
  const Atmosphere& atmo = detector.GetAtmosphere();
  const FDetector& detFD = detector.GetFDetector();
  fevt::FEvent& fevent = event.GetFEvent();

  // Histograms for debugging !!!!!!!
  TH2D* hist2D = new TH2D("hist2D", "hist2D", 600, -0e3, 10e3, 600, 50e3, 60e3); // photon density [km]
  TH2D* histObs2D = new TH2D("histObs2D", "histObs2D", 600, -0e3, 10e3, 600, 50e3, 60e3); // photon density [km]
  TH2D* histCore2D = new TH2D("histCore2D", "histCore2D", 300, -40e3, 40e3, 300, -40e3, 40e3); // photon density [km]
  TH2D* hist2Demit = new TH2D("hist2Demit_site", "hist2Demit_site", 300, -5e3, 75e3, 300, 0e3, 85e3); // photon density [km]
  TH2D* histObs2Demit = new TH2D("histObs2Demit_site", "histObs2Demit_site", 300, -5e3, 75e3, 300, 0e3, 85e3); // photon density [km]
  TH2D* histCore2Demit = new TH2D("histCore2Demit_site", "histCore2Demit_site", 300, -40e3, 40e3, 300, -40e3, 40e3); // photon density [km]
  const double deltaZoom = 5.e3; // meters
  TH2D* hist2DzoomLL = new TH2D("hist2DzoomLL", "hist2DzoomLL", 300, -deltaZoom, deltaZoom, 300, -deltaZoom, deltaZoom); // photon density [km]
  TH2D* hist2DzoomLA = new TH2D("hist2DzoomLA", "hist2DzoomLA", 300, -deltaZoom, deltaZoom, 300, -deltaZoom, deltaZoom); // photon density [km]
  TH2D* hist2DzoomLM = new TH2D("hist2DzoomLM", "hist2DzoomLM", 300, -deltaZoom, deltaZoom, 300, -deltaZoom, deltaZoom); // photon density [km]
  TH2D* hist2DzoomCO = new TH2D("hist2DzoomCO", "hist2DzoomCO", 300, -deltaZoom, deltaZoom, 300, -deltaZoom, deltaZoom); // photon density [km]
  TH2D* angle2D = new TH2D("angle2D", "angle2D", 100, -1, 1, 100, -1, 1);
  TH2D* onTel2D = new TH2D("onTel2D", "onTel2D(telCS)", 500, -100, 100, 500, -100, 100); // photon density [m]
  TH1D* histCos = new TH1D("histCos", "histCos", 100,-1,1);
  TH1D* histCosEmission = new TH1D("histCosEmission", "histCosEmission", 100,-1,1);
  TH1D* histTimeCher = new TH1D("histTimeCher", "histTimeCher", 100,1,0);
  TH1D* histTimeArrival = new TH1D("histTimeArrival", "histTimeArrival", 100,1,0);
  vector<double> eyesX, eyesY;

  map<int, fevt::TelescopeSimData*> tel_sim;
  map<int, map<int, double>> tel_time_bin_trace;
  map<int, double> tel_time_bin_width;
  map<int, double> tel_time_nbins;

  /*
     keep this stupid snipplet of code, since it produces the data needed for CORSIKA steering cards:
  */
  /*
  cout << "*            \tx-pos [cm]   \ty-pos [cm]   \tz-pos[cm] \tradius [cm]    ID       comment" << endl;
  for (fevt::FEvent::EyeIterator iEye = fevent.EyesBegin(fevt::ComponentSelector::eExists);
       iEye != fevent.EyesEnd(fevt::ComponentSelector::eExists) ; ++iEye) {

    const CoordinateSystemPtr refCS = detector.GetReferenceCoordinateSystem();
    const UTMPoint theUTMPoint = UTMPoint(Point(0,0,0,refCS), ReferenceEllipsoid::eWGS84);
    const double origin_height = theUTMPoint.GetHeight();

    const fdet::Eye& detEye = detFD.GetEye(*iEye);
    const utl::Point eyePosition = detEye.GetPosition();
    //eyesX.push_back(eyePosition.GetX(refCS)/meter);
    //eyesY.push_back(eyePosition.GetY(refCS)/meter);
    cout << std::fixed << std::setprecision(2);
    cout << "TELESCOPE \t"
         << setw(10) << eyePosition.GetX(refCS)/centimeter << " \t"
         << setw(10) << eyePosition.GetY(refCS)/centimeter << " \t"
         << setw(10) << (eyePosition.GetZ(refCS)+origin_height)/centimeter << " \t"
         << setw(7) << 2*meter/centimeter << " \t" << 0 << "    \t" <<  detEye.GetName() <<endl;

  }
  */
  /* -------  */


  map<int, int> countPhot;
  map<int, double> weightPhot;
  map<int, double> weightSqrPhot;
  unsigned int countAll = 0;

  /*
    The most economic way to look photons and telescopes is to first
    make a list of all telescope eye combinations, then later start to
    read photons on telescope level first, then on eye level and last
    on global level. Telescope that have been already processed are
    skipped.
   */

  map<int, map<int, bool> > doneThis;

  vector<utl::AttributeMap> allTelEyeCombinations;
  for (fevt::FEvent::EyeIterator iEye = fevent.EyesBegin(fevt::ComponentSelector::eInDAQ);
       iEye != fevent.EyesEnd(fevt::ComponentSelector::eInDAQ) ; ++iEye) {

    fevt::Eye& eyeEvent = *iEye;

    if (eyeEvent.GetStatus() == fevt::ComponentSelector::eDeSelected)
      continue;

    unsigned int eyeId = eyeEvent.GetId();

    for (fevt::Eye::TelescopeIterator iTel = eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
         iTel != eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++iTel) {

      unsigned int telId = iTel->GetId();
      fevt::Telescope& telEvent = *iTel;

      if (!telEvent.HasSimData())
        continue;
      fevt::TelescopeSimData& telSim = telEvent.GetSimData();
      if (!telSim.HasDistanceTrace())
        continue;
      utl::AttributeMap thisDet;
      thisDet["eye"] = std::to_string(eyeId);
      thisDet["tel"] = std::to_string(telId);
      allTelEyeCombinations.push_back(thisDet);
    } // loop tels
    utl::AttributeMap thisDet;
    thisDet["eye"] = std::to_string(eyeId);
    allTelEyeCombinations.push_back(thisDet); // now catch all the remaining telescopes of this eye
  } // loop eyes
  utl::AttributeMap thisDet;
  allTelEyeCombinations.push_back(thisDet); // now catch everything left over



  for (const auto& itCombi:allTelEyeCombinations) {

    utl::AttributeMap am = itCombi;

    const bool hasTelAttribute = am.count("tel");
    const unsigned int telAttribute = hasTelAttribute ? std::stoi(am["tel"]) : 0;
    const bool hasEyeAttribute = am.count("eye");
    const unsigned int eyeAttribute = hasEyeAttribute ? std::stoi(am["eye"]) : 0;

    if (!simShower.HasGroundCherenkov(am))
      continue;

    /*
    cerr << "found iterator for " << (hasEyeAttribute ? am["eye"] : "-") << " " << (hasTelAttribute ? am["tel"] : "-") << " " << eyeAttribute << " " << telAttribute << endl;

    for (auto flagEye : doneThis)
      for (auto flagTel : flagEye.second)
        cout << " doneThis[" << flagEye.first << "][" << flagTel.first <<"] = true" << endl;
    */

    map<int, map<int, bool> > flagThis;
    unsigned int countThis = 0;

    for (utl::ShowerParticleIterator iCher = simShower.GroundCherenkovBegin(am);
         iCher != simShower.GroundCherenkovEnd(am);
         ++iCher) {

      if (iCher->GetType() != utl::Particle::ePhoton) {
        ostringstream err;
        err << "Cherenkov input particle not of type \'photon\': " << iCher->GetName();
        ERROR(err);
        continue;
      }

      ++countAll;
      ++countThis;

      const Vector& nIn = iCher->GetDirection();
      const Point& pIn = iCher->GetPosition();
      double weight = iCher->GetWeight();
      const double energy = iCher->GetTotalEnergy();

      if (true) { // debugging only
        // wow, my first C++ lambda function ! cool stuff.
        std::function<void (CoordinateSystemPtr,TH2D*,TH2D*)> func = [pIn,nIn,weight](CoordinateSystemPtr csDet, TH2D* hist, TH2D* hist2) {
          const Point site(0,0,0,csDet);
          const CoordinateSystemPtr refCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
          Vector shift = site - Point(0,0,0,refCS);
          // now move refCS origin to core
          CoordinateSystemPtr cs = CoordinateSystem::Translation(shift, refCS);

          const Point site_pos = site;
          const Vector site_normal(0,0,1,cs);
          const Vector emission_location = site_pos - pIn;
          // const double emission_distance = emission_location.GetMag(); // positive // unused -- silence warning
          const double distance_perp = emission_location * site_normal; // negative value
          const double cosTheta = site_normal*nIn;
          const double distance = distance_perp / cosTheta; // positive value
          const Point pGround = pIn + distance*nIn;
          if (hist2) {
            hist2->Fill(pIn.GetX(cs)/meter, pIn.GetY(cs)/meter, weight);
          }
          hist->Fill(pGround.GetX(cs)/meter, pGround.GetY(cs)/meter, weight);
        };

        const CoordinateSystemPtr llCS = detFD.GetEye("Los Leones").GetEyeCoordinateSystem();
        const CoordinateSystemPtr laCS = detFD.GetEye("Loma Amarilla").GetEyeCoordinateSystem();
        const CoordinateSystemPtr lmCS = detFD.GetEye("Los Morados").GetEyeCoordinateSystem();
        //const CoordinateSystemPtr coCS = detFD.GetEye("Coihueco").GetEyeCoordinateSystem();
        const CoordinateSystemPtr coCS = detFD.GetEye(4).GetEyeCoordinateSystem();

        CoordinateSystemPtr coreCS = simShower.GetLocalCoordinateSystem();
        const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();
        CoordinateSystemPtr obsCS = CoordinateSystem::Translation(Vector(0,0,1300.*meter,siteCS),siteCS);

        func(coreCS, histCore2D, histCore2Demit);
        func(siteCS, hist2D, hist2Demit);
        func(obsCS, histObs2D, histObs2Demit);
        func(llCS, hist2DzoomLL,0);
        func(laCS, hist2DzoomLA,0);
        func(lmCS, hist2DzoomLM,0);
        func(coCS, hist2DzoomCO,0);
      } // end debugging


      for (fevt::FEvent::EyeIterator iEye = fevent.EyesBegin(fevt::ComponentSelector::eInDAQ);
           iEye != fevent.EyesEnd(fevt::ComponentSelector::eInDAQ) ; ++iEye) {

        fevt::Eye& eyeEvent = *iEye;

        if (eyeEvent.GetStatus() == fevt::ComponentSelector::eDeSelected)
          continue;

        const unsigned int eyeId = eyeEvent.GetId();
        if (hasEyeAttribute && eyeAttribute != eyeId)
          continue;

        for (fevt::Eye::TelescopeIterator iTel = eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
             iTel != eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++iTel) {

          unsigned int telId = iTel->GetId();
          if (hasTelAttribute && telAttribute != telId)
            continue;

          if (doneThis.count(eyeId) && doneThis[eyeId].count(telId) && doneThis[eyeId][telId])
            continue; // already done this

          flagThis[eyeId][telId] = true; // remember this... only process each telescope exactly one time

          fevt::Telescope& telEvent = *iTel;

          if (!telEvent.HasSimData())
            continue;
          fevt::TelescopeSimData& telSim = telEvent.GetSimData();

          if (!telSim.HasDistanceTrace())
            continue;
          TraceD& distanceTrace = telSim.GetDistanceTrace();

          const int tel_unique = eyeId * 100000 + telId;

          const fdet::Telescope& detTel = detFD.GetTelescope(telEvent);
          const double rDia = detTel.GetDiaphragmRadius();
          const Point telPosition = detTel.GetPosition();
          const Vector telNormal = detTel.GetAxis();
          const CoordinateSystemPtr telCS = detTel.GetTelescopeCoordinateSystem();

          const double cosTheta = telNormal*nIn;

          const Vector emission_location = telPosition - pIn;
          const double emission_distance = emission_location.GetMag();       // positive value
          const double distance_perp = emission_location * telNormal;        // negative value

          const double cosThetaEmission = distance_perp / emission_distance; // negative

          histCos->Fill(cosTheta);
          histCosEmission->Fill(cosThetaEmission);
          angle2D->Fill(cosThetaEmission,cosTheta);
          const double cosFOV = cos(detTel.GetCamera().GetFieldOfView());
          if (-cosThetaEmission<cosFOV) { // cannot hit, FOV-Auger is ~16 deg (FModelConfig: 23.07deg (?) )
            continue;
          }
          // check photon travel direction. Is it in FOV of telescope?
          if (-cosTheta<cosFOV) { // cannot hit, FOV-Auger is ~16 deg (FModelConfig: 23.07deg (?) )
            continue;
          }

          const double distance = distance_perp / cosTheta; // positive value
          const Point onDia = pIn + distance*nIn;
          const double radial = sqrt(pow(onDia.GetX(telCS),2) + pow(onDia.GetY(telCS),2));

          onTel2D->Fill(onDia.GetX(telCS)/meter, onDia.GetY(telCS)/meter, weight);

          if (radial>rDia) { // did not hit
            continue;
          }

          // attenuation
          const double wavelength = utl::kPlanck / energy;
          const double rAttCkov = atmo.EvaluateRayleighAttenuation(pIn, onDia, wavelength);
          const double mAttCkov = atmo.EvaluateMieAttenuation(pIn, onDia, wavelength);
          weight *= rAttCkov * mAttCkov;

          countPhot[tel_unique] ++;
          weightPhot[tel_unique] += weight;
          weightSqrPhot[tel_unique] += weight*weight;

          tel_sim[tel_unique] = &telSim;

          const TimeStamp traceAtTelTime = telSim.GetPhotonsStartTime();
          const TimeStamp timePhoton = timeAtCore + iCher->GetTime() + utl::TimeInterval(distance/utl::kSpeedOfLight);
          const TimeInterval timeDia = timePhoton - traceAtTelTime;

          histTimeCher->Fill(iCher->GetTime()/ns);
          histTimeArrival->Fill(timeDia/ns);

          /*
            static int deli = 0;
            if (deli++<100) {
            ostringstream inf;
            inf << "phot-time prod_tel=" << distance/utl::kSpeedOfLight/nanosecond << "ns  coreToStart=" << coreToStart/ns << "ns, tMin=" << iTelFOV->Tmin/nanosecond << "ns ";
            //inf << "phot-time prod_tel=" << distance_perp/utl::kSpeedOfLight/nanosecond << "ns  coreToStart=" << coreToStart/m << "m, tMin=" << iTelFOV->Tmin/nanosecond << "ns ";
            //    inf << " t_abs=" << timeDia/ns << "ns ";
            inf << ", t_cher_prod=" << iCher->GetTime().GetInterval()/ns << "ns ";
            inf << ", t_tel_start=" << traceAtTelTime.GetInterval()/ns << "ns ";
            inf << ", t_cher_abs=" << timePhoton.GetInterval()/ns << "ns ";
            inf << ", t_cher_rel=" << timeDia.GetInterval()/ns << "ns ";
            inf << ", cosTheta=" << cosTheta << " cosThetaEmission=" << cosThetaEmission;
            inf << ", z_prod=" << pIn.GetZ(siteCS)/m;
            inf << ", distance_perp=" << distance_perp/m << " distance=" << distance/m << " radial=" << radial/m;
            inf << ", nIn=" << nIn.GetX(siteCS) << "," << nIn.GetY(siteCS) << "," << nIn.GetZ(siteCS);
            INFO(inf);
            }
          */

          const fevt::FdConstants::LightSource lightSource = fevt::FdConstants::eCherDirect;
          utl::Photon photonIn(onDia, nIn, wavelength, weight, lightSource);
          photonIn.SetTime(timeDia);
          telSim.AddPhoton(photonIn);


          // the code in the following code is NOT essential. It
          // produces (if working) only the sim-light traces in the
          // EventBrowser. The simualtion does not depend on it.
          {
            //do not recalculate weight to keep it for other telescopes
            //weight /= detTel.GetDiaphragmArea() * (-cosTheta);
            double weight_loc = weight / (detTel.GetDiaphragmArea() * (-cosTheta));

            const vector<double>& WLengthCkov = atmo.GetWavelengths(atm::Atmosphere::eCerenkov);
            const int CkovBin = WLengthCkov.size()/2;
            static const double telEffWL = detTel.GetModelRelativeEfficiency(WLengthCkov[CkovBin]);
            try {
              weight_loc *= detTel.GetModelRelativeEfficiency(wavelength) / telEffWL;
            } catch (utl::OutOfBoundException&) {
              weight_loc = 0;
            }

            TimeInterval timeDiaAbs = timeDia;

            if (tel_time_bin_trace.count(tel_unique) == 0) {
              tel_time_bin_width[tel_unique] = distanceTrace.GetBinning();
              tel_time_nbins[tel_unique] = distanceTrace.GetSize();
              tel_time_bin_trace[tel_unique] = map<int, double>();
            }
            const int tel_time_bin = timeDiaAbs.GetInterval() / tel_time_bin_width[tel_unique];
            tel_time_bin_trace[tel_unique][tel_time_bin] += weight_loc;
          } // end light binningv-------------------

          //restore weight to original value to process photon in other telescopes
          weight /= rAttCkov * mAttCkov;

          // done with this photon -> next photon - update: do not skip other telescopes
          //goto next_photon;

        } // loop tels
      } // loop eyes
    //next_photon: ;
    } // end loop cherenkov photons

    // remember for next cherenkov input file
    for (const auto& flagEye : flagThis)
      for (const auto& flagTel : flagEye.second)
        doneThis[flagEye.first][flagTel.first] = true; // remember this... only process each telescope exactly one time

    ostringstream info;
    info << "Read " << countThis << " photons with attributes: eye=\""
         << (hasEyeAttribute ? am["eye"] : "-") << "\", tel=\"" << (hasTelAttribute ? am["tel"] : "-") << "\" ";
    INFO(info.str());

  } // end loop cherenkov files

  // final output
  bool lowQualityProblem = false;

  for (map<int,double>::const_iterator iDT = tel_time_bin_width.begin();
       iDT != tel_time_bin_width.end(); ++iDT) {
    const int tel_unique = iDT->first;  //  eyeId * 100000 + telId;
    unsigned int eyeId = tel_unique  /100000;
    unsigned int telId = tel_unique % 100000;
    ostringstream info;
    if (countPhot.count(tel_unique)) {
      const double weight = weightPhot[tel_unique] / countPhot[tel_unique];
      const double weightRMS = sqrt(weightSqrPhot[tel_unique] / countPhot[tel_unique] - weight*weight);
      info << "Eye=" << eyeId << ", Tel=" << telId << ": ";
      info << "detected " << countPhot[tel_unique] << " photons on diaphragm with total weight= " << weightPhot[tel_unique];
      info << ": <weight>=" << weight << ", RMS(weight)=" << weightRMS;
      if (weight>10000 || weightRMS>1000) {
        lowQualityProblem = true;
        info << " [WEIGHTS TOO LARGE! CANNOT HANDLE!]";
      }
      INFO(info);
    }
  }

  if (lowQualityProblem) {
    ostringstream err;
    err << "Please consider to run CORSIKA with less thinning. Cherenkov simulation makes no sense with huge weights";
    ERROR(err);
    throw AugerException(err.str());
  }

  {
    ostringstream inf;

    std::locale global;
    std::locale withgroupings(global, new my_facet);
    inf.imbue(withgroupings);

    inf << "Found a total of " << countAll << " Cherenkov photons in input file" << endl;
    INFO(inf);
  }

  // RUDEBUG ------
  // some plots
  static int iSave = 0;
  iSave++;
  ostringstream savename;
  savename << "save_" << iSave << ".root";

  TDirectory* saveDir = gDirectory;
  TFile* file = TFile::Open(savename.str().c_str(), "recreate");
  file->cd();

  TCanvas* can2D = new TCanvas("can2D");
  can2D->Divide(4,2);

  can2D->cd(1);
  gPad->SetLogz(1);
  hist2Demit->Draw("colz");

  can2D->cd(2);
  gPad->SetLogz(1);
  hist2D->Draw("colz");

  can2D->cd(3);
  gPad->SetLogz(1);
  histCore2D->Draw("colz");

  can2D->cd(4);
  gPad->SetLogz(1);
  histObs2D->Draw("colz");

  can2D->cd(5);
  gPad->SetLogz(1);
  hist2DzoomLL->Draw("colz");

  can2D->cd(6);
  gPad->SetLogz(1);
  hist2DzoomLA->Draw("colz");

  can2D->cd(7);
  gPad->SetLogz(1);
  hist2DzoomLM->Draw("colz");

  can2D->cd(8);
  gPad->SetLogz(1);
  hist2DzoomCO->Draw("colz");

  TGraph* geyes = new TGraph(eyesX.size(), &eyesX.front(), &eyesY.front());
  geyes->SetMarkerStyle(20);
  geyes->SetMarkerSize(1.5);

  can2D->cd(1);
  geyes->Draw("p");
  can2D->cd(2);
  geyes->Draw("p");
  can2D->cd(3);
  geyes->Draw("p");
  can2D->cd(4);
  geyes->Draw("p");
  /*
    can2D->cd(3);
    geyes->Draw("p");
    can2D->cd(4);
    geyes->Draw("p");
    can2D->cd(5);
    geyes->Draw("p");
    can2D->cd(6);
    geyes->Draw("p");
  */

  can2D->Write();

  hist2D->Write();
  histObs2D->Write();
  histCore2D->Write();
  hist2Demit->Write();
  histObs2Demit->Write();
  histCore2Demit->Write();
  hist2DzoomLL->Write();
  hist2DzoomLA->Write();
  hist2DzoomLM->Write();
  hist2DzoomCO->Write();
  onTel2D->Write();
  histCosEmission->Write();
  histCos->Write();
  angle2D->Write();
  histTimeArrival->Write();
  histTimeCher->Write();

  file->Write();
  file->Close();
  saveDir->cd();
  // ------------- end plotting ----------

  // save light traces at diaphragm. This is just for EventBrowser
  map<int,double>::const_iterator iDT = tel_time_bin_width.begin();
  map<int,map<int,double>>::const_iterator iPh = tel_time_bin_trace.begin();
  map<int,fevt::TelescopeSimData*>::const_iterator iTS = tel_sim.begin();

  for (; iDT!=tel_time_bin_width.end(); ++iDT, ++iPh, ++iTS) {

    const vector<double>& WLengthCkov = atmo.GetWavelengths(atm::Atmosphere::eCerenkov);
    const int CkovBin = WLengthCkov.size()/2; // all in one central wavelength bin

    fevt::TelescopeSimData& telSim = *(iTS->second);

    if (telSim.HasPhotonTrace(fevt::FdConstants::eCherDirect, CkovBin)) {
      ERROR("ERROR too! RUDEBUG ");
      continue;
    }

    const int firstBin = (iPh->second).begin()->first;
    const int nBins = (iPh->second).rbegin()->first - firstBin + 1;

    telSim.MakePhotonTrace(fevt::FdConstants::eCherDirect, CkovBin, nBins, iDT->second);
    TraceD& cherDirectTrace = telSim.GetPhotonTrace (fevt::FdConstants::eCherDirect, CkovBin);
    for (map<int,double>::const_iterator iBin=iPh->second.begin(); iBin!=iPh->second.end(); ++iBin) {
      //cherDirectTrace[iBin->first - firstBin] = iBin->second;
      if (iBin->first >= 0 && iBin->first < nBins-1)
        cherDirectTrace[iBin->first] = iBin->second;
    }
  }
}