ShowerPhotonGenerator.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "ShowerPhotonGenerator.h"
#include "MultipleScatterer.h"

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/Math.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/Trace.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/Photon.h>
#include <utl/RandomEngine.h>
#include <utl/RandomSamplerFromPDF.h>
#include <utl/UTMPoint.h>

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

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Mirror.h>
#include <fdet/Filter.h>
#include <fdet/Corrector.h>
#include <fdet/Camera.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/LaserData.h> // Modified by D'Urso & Valore 02/03/06

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

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

#include <CLHEP/Random/Randomize.h>

#include <boost/tuple/tuple.hpp>

#include <iomanip>
#include <fstream>
#include <sstream>
#include <map>

#include <TROOT.h>
#include <TStyle.h>
#include <TH1D.h>
#include <TCanvas.h>
#include <TLegend.h>

// #define DEBUG 1
// #define DEBUG_3D 1

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

using CLHEP::RandFlat;



ShowerPhotonGenerator::ShowerPhotonGenerator() :
  fUseOnlyReferenceWavelength(false),
  fFluorescenceLDF(eOff),
  fScatteredCherenkovLDF(eOff),
  fDirectCherenkovLDF(eOff),
  fLightSourceSelector(-2),
  fMultipleScatterer(NULL)
{
}


VModule::ResultFlag
ShowerPhotonGenerator::Init()
{
  CentralConfig* const cc = CentralConfig::GetInstance();
  Branch topB = cc->GetTopBranch("ShowerPhotonGenerator");

  fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

  topB.GetChild("maxNRayTracePerBin").GetData(fMaxNRayTrace);
  topB.GetChild("minNRayTracePerBin").GetData(fMinNRayTrace);
  topB.GetChild("extraRayTraceFactor").GetData(fExtraRayTraceFactor);

  string fluorescence3D;
  topB.GetChild("FluorescenceLDF").GetData(fluorescence3D);
  if (fluorescence3D == "off")
    fFluorescenceLDF = eOff;
  else if (fluorescence3D == "NKG")
    fFluorescenceLDF = eNKG;
  else if (fluorescence3D == "Gora")
    fFluorescenceLDF = eGora;
  else {
    ostringstream err;
    err << " Unknown fluorescence LDF mode: \"" << fluorescence3D;
    ERROR(err);
    return eFailure;
  }

  string cherenkov3D;
  topB.GetChild("DirectCherenkovLDF").GetData(cherenkov3D);
  if (cherenkov3D == "off")
    fDirectCherenkovLDF = eOff;
  else if (cherenkov3D == "NKG")
    fDirectCherenkovLDF = eNKG;
  else if (cherenkov3D == "Gora")
    fDirectCherenkovLDF = eGora;
  else if (cherenkov3D == "Cherenkov")
    fDirectCherenkovLDF = eCherenkov;
  else if (cherenkov3D == "CherenkovAtAxis")
    fDirectCherenkovLDF = eCherenkovAtAxis;
  else if (cherenkov3D == "CORSIKA")
    fDirectCherenkovLDF = eCherenkovCORSIKA;
  else {
    ostringstream err;
    err << " Unknown direct cherenkov LDF option: " << cherenkov3D;
    ERROR(err);
    return eFailure;
  }
  topB.GetChild("ScatteredCherenkovLDF").GetData(cherenkov3D);
  if (cherenkov3D == "off")
    fScatteredCherenkovLDF = eOff;
  else if (cherenkov3D == "NKG")
    fScatteredCherenkovLDF = eNKG;
  else if (cherenkov3D == "Gora")
    fScatteredCherenkovLDF = eGora;
  else if (cherenkov3D == "Cherenkov")
    fScatteredCherenkovLDF = eCherenkov;
  else if (cherenkov3D == "CherenkovAtAxis")
    fScatteredCherenkovLDF = eCherenkovAtAxis;
  else {
    ostringstream err;
    err << " Unknown scattered cherenkov LDF option: " << cherenkov3D;
    ERROR(err);
    return eFailure;
  }
  topB.GetChild("MultipleScattering").GetData(fMultipleScattering);

  // ---- info output -----
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "    max raytraced photons per time bin: " << fMaxNRayTrace<< "\n"
          "    min raytraced photons per time bin: " << fMinNRayTrace<< "\n"
          "   artificial factor for extra photons: " << fExtraRayTraceFactor<< "\n"
          "                      fluorescence LDF: ";
  switch (fFluorescenceLDF) {
  case eOff: info << "off\n"; break;
  case eNKG: info << "NKG\n"; break;
  case eGora: info << "Gora et al.\n"; break;
  default: info << "UNKNOWN\n"; break;
  }
  info << "                  direct Cherenkov LDF: ";
  switch (fDirectCherenkovLDF) {
  case eOff: info << "off\n"; break;
  case eNKG: info << "NKG\n"; break;
  case eGora: info << "Gora\n"; break;
  case eCherenkov: info << "full sim. Cherenkov emission\n"; break;
  case eCherenkovAtAxis: info << "Cherenkov sim. for electrons at axis\n"; break;
  case eCherenkovCORSIKA: info << "Cherenkov input from CORSIKA\n"; break;
  }
  info << "               scattered Cherenkov LDF: ";
  switch (fScatteredCherenkovLDF) {
  case eOff: info << "off\n"; break;
  case eNKG: info << "NKG\n"; break;
  case eGora: info << "Gora\n"; break;
  case eCherenkov: info << "full sim. Cherenkov emission\n"; break;
  case eCherenkovAtAxis: info << "Cherenkov sim. for electrons at axis\n"; break;
  case eCherenkovCORSIKA: info << "eCherenkovCORSIKA?\n"; break;
  }
  info << "                   multiple scattering: " << (fMultipleScattering ? "ON" : "OFF") << "\n";
  INFO(info);

  if (fMultipleScattering)
     fMultipleScatterer = new MultipleScatterer();

  if (fDirectCherenkovLDF == eCherenkov || fDirectCherenkovLDF == eCherenkovAtAxis) {
//#warning Sim. Cherenkov LDF for direct Cherenkov emission not implemeted yet.
    ERROR("Sim. Cherenkov LDF for direct Cherenkov emission not implemeted yet.");
    return eFailure;
  }

  return eSuccess;
} // end of Init


// ************************************************************************
VModule::ResultFlag
ShowerPhotonGenerator::Run(evt::Event& event)
{
  if (!event.HasFEvent()) {
    ERROR ("Event has no FEvent.");
    return eContinueLoop;
  }
  fevt::FEvent& fEvent = event.GetFEvent();
  
  if (!event.HasSimShower() ||
      !event.GetSimShower().HasDirection() ||
      (!event.GetSimShower().HasGHParameters() && !event.GetSimShower().HasLaserData())) {
    ERROR ("Event has no sim-data.");
    return eContinueLoop;
  }
  
  evt::ShowerSimData& simShower = event.GetSimShower();
  const CoordinateSystemPtr& showerCS = simShower.GetShowerCoordinateSystem();
  
  const bool isLaser = simShower.HasLaserData();
  
  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();
  const Point showerCore(0, 0, 0, showerCS);
  const Vector showerAxis(0, 0, -1, showerCS);

  const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
  const double cosZenith = (-simShower.GetDirection()).GetCosTheta(localCS); // just for flat earth
  const double absCosZenith = std::fabs(cosZenith); // just for flat earth
  const double coreZ = wgs84.PointToLatitudeLongitudeHeight(showerCore).get<2>();
  
  Detector& detector = Detector::GetInstance();
  const FDetector& detFD = detector.GetFDetector();
  const Atmosphere& atmo = detector.GetAtmosphere();
  
  const double refWl = detFD.GetReferenceLambda();
  const vector<double>& WLengthFluo = atmo.GetWavelengths(Atmosphere::eFluorescence);
  const vector<double>& WLengthCkov = atmo.GetWavelengths(Atmosphere::eCerenkov);
  const int nWLengthCkov = WLengthCkov.size();
  
  const ProfileResult& depthProfile = atmo.EvaluateDepthVsHeight();
  const ProfileResult& temperatureProfile = atmo.EvaluateTemperatureVsHeight();
  const ProfileResult& pressureProfile = atmo.EvaluatePressureVsHeight();
  
  const double tempHeightMin = temperatureProfile.MinX();
  const double tempHeightMax = temperatureProfile.MaxX();
  const double pressHeightMin = pressureProfile.MinX();
  const double pressHeightMax = pressureProfile.MaxX();
  
  bool flatEarth = false;
  ProfileResult const * slantDepthVsDistance = 0;
  double atmDistanceMin = 0;
  double atmDistanceMax = 0;
  if (!isLaser && (fFluorescenceLDF != eOff ||
                   (fDirectCherenkovLDF != eOff) ||
                   (fScatteredCherenkovLDF != eOff))) {
    try {
      slantDepthVsDistance = &atmo.EvaluateSlantDepthVsDistance();
    } catch (InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
      flatEarth = true;
    }
    atmDistanceMin = flatEarth ? (coreZ - depthProfile.MinX()) / absCosZenith : slantDepthVsDistance->MinX();
    atmDistanceMax = flatEarth ? (coreZ - depthProfile.MaxX()) / absCosZenith : slantDepthVsDistance->MaxX();
  }

  double X1 = 0;
  double showerXmax = 0;
  double laserWavelength = 0;
  if (!isLaser) {
    X1 = simShower.GetXFirst();
    showerXmax = simShower.GetGHParameters().GetXMax();
  } else {
    laserWavelength = simShower.GetLaserData().GetLaserWavelength();
  }

  // For making nice LDF-pictures -----------
  // Add shower ages to following vector to
  // produce LDF pictures
  vector<double> ages;
  //ages.push_back(0.7);
  //ages.push_back(0.9);
  //ages.push_back(1.0);
  //ages.push_back(1.2);
  //ages.push_back(1.3);
  for (unsigned int i = 0; i < ages.size(); ++i) {
    const double agePlot = ages[i];
    PlotLDF(simShower, simShower.GetCherenkovBeamProductionPhotons(nWLengthCkov/2), agePlot);
  } // -------------------------------------------


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

    if (iEye->GetStatus() == fevt::ComponentSelector::eDeSelected)
      continue;

    unsigned int eyeId = iEye->GetId();
    fevt::Eye& eyeEvent = *iEye;

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

      const unsigned int telId = iTel->GetId();

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

      const double rDia = detTel.GetDiaphragmRadius();
      const double diaphragmArea = detTel.GetDiaphragmArea();
      
      if (!telEvent.HasSimData())
        continue;      
      fevt::TelescopeSimData& telSim = telEvent.GetSimData();
      
      if (!telSim.HasDistanceTrace())
        continue;
      const TraceD& distanceTrace = telSim.GetDistanceTrace();
      
      ostringstream info;
      info << "Photons for eye=" << eyeId
           << ", telescope=" << telId
           << ", rDia:" << setw(4) << detTel.GetDiaphragmRadius();
      unsigned int countRTPhotons = 0;
      double totalWeightPhotons = 0;
      
      
      const CoordinateSystemPtr& telCS = detTel.GetTelescopeCoordinateSystem();
      const Point telescopePos(0.0, 0.0, 0.0, telCS);
      // RU ---- nkg option ---------
      const Point p1(telescopePos);
      const Vector telescopeToShowerCore(showerCore-p1);
      // const double distTelToCore = telescopeToShowerCore.GetMag();
      // ----------------------------
      
      // definition of the light trace at the diaphragm
      const unsigned int nBins = distanceTrace.GetSize();
      const double tracebin = distanceTrace.GetBinning(); // [ns]
      
      if (!telSim.HasRayTracedPhotonTrace())
        telSim.MakeRayTracedPhotonTrace(nBins, tracebin);
      TraceI& rayTracedPhotonTrace = telSim.GetRayTracedPhotonTrace();
      
      const fdet::Camera& detCamera = detTel.GetCamera();
      const double telTraceBinWidth = detCamera.GetFADCBinSize();
      const double totalTraceDuration = tracebin * nBins;
      const unsigned int telTraceNBins = int(0.5 + totalTraceDuration / telTraceBinWidth);
      telSim.SetNumberOfPhotonBins(telTraceNBins);
      
#ifdef DEBUG
      info << " photonTraceSize=" << nBins
           << " tracebin=" << tracebin
           << " telTraceBinWidth=" << telTraceBinWidth
           << " telTraceNBins=" << telTraceNBins << endl;
#endif


      // --- main loop along the shower axis -------------------------

      if (nBins<2) 
        continue;

      double previousDistance = 0;
      for (unsigned int iTrace = 0; iTrace < nBins; ++iTrace) {
        
        unsigned int totalNRayTracedPhotonsPerBin = 0;
        
        // smear photons along distance equally
        const bool firstBin =  iTrace <= 0;
        const bool lastBin = iTrace >= nBins-1;
                
        const double currentDistance = distanceTrace[iTrace];
        const double deltaNext = (lastBin ? 0. : distanceTrace[iTrace+1]-currentDistance);
        const double deltaPrev = (firstBin ? 0. : currentDistance-previousDistance);
        const double binFraction = RandFlat::shoot(&fRandomEngine->GetEngine(), -0.5, +0.5); // where exactly in distance bin
        const double delta = (binFraction<0 ? deltaPrev : deltaNext); // distance to next/previous bin
        double distancePhoton = currentDistance + delta * binFraction;
        previousDistance = currentDistance;

        const double traceTime = tracebin * (0.5 + binFraction + iTrace);

        // some quantities in the middle of the current trace bin
        // that need not be re-calculated for each photon
        // const double averageDistancePhoton = 0.5*(distanceBin + distanceTrace[min(iTrace+1,nBins-1)]);
        const Point pointOnShower = showerCore + showerAxis * distancePhoton;
        const double height = wgs84.PointToLatitudeLongitudeHeight(pointOnShower).get<2>();
        const double temperature = temperatureProfile.Y(std::max(std::min(height, tempHeightMax),
                                                                 tempHeightMin));
        const double pressure = pressureProfile.Y(std::max(std::min(height, pressHeightMax),
                                                           pressHeightMin));
        const double cosLocalZenith =
          (-showerAxis).GetCosTheta(LocalCoordinateSystem::Create(UTMPoint(pointOnShower, wgs84)));
        const double moliereR = MoliereRadius(temperature, pressure,
                                              flatEarth ?  cosZenith : cosLocalZenith);

        // build probability distribution for wavelengths AND light sources
        map<double, double> sourcePDF;
        map<int, RandomSamplerFromPDF*> randomWavelengthSource;
        ScopeGuard wlGuard(randomWavelengthSource); // will free the RandomSamplerFromPDF objects on scope exit

        double totalSumPhotons = 0;
        for (fevt::TelescopeSimData::PhotonTraceSourceIterator iSource = telSim.PhotonTracesSourceBegin(),
             end = telSim.PhotonTracesSourceEnd();
             iSource != end; ++iSource) {

          // this can be used to select a single light source for tests:
          if (GetLightSourceSelection() >= 0 && (int)(iSource->first) != GetLightSourceSelection())
            continue;
          
          // if direct-Cherenkov photons are to be taken from CORSIKA directly, IGNORE them here: 
          if (fDirectCherenkovLDF==eCherenkovCORSIKA && iSource->first == fevt::FdConstants::eCherDirect) 
            continue; 
          
          map<double, double> sourceWavelengthPDF;
          
          double totalSumPhotonsSource = 0;
          for (fevt::TelescopeSimData::PhotonTraceIterator iBin = telSim.PhotonTracesBegin(iSource->first),
               end = telSim.PhotonTracesEnd(iSource->first);
               iBin != end; ++iBin) {
            const TraceD& photonTrace = iBin->GetTrace();
            const double nPhotonsWL = photonTrace[iTrace];
            if (nPhotonsWL > 0) {
              sourceWavelengthPDF[iBin->GetLabel()] = nPhotonsWL;
              totalSumPhotonsSource += nPhotonsWL;
            }
          }

          if (totalSumPhotonsSource <= 0)
            continue;

          if (!laserWavelength) {
            randomWavelengthSource[int(iSource->first)] = new RandomSamplerFromPDF(sourceWavelengthPDF);
          }
          sourcePDF[int(iSource->first)] = totalSumPhotonsSource;
          totalSumPhotons += totalSumPhotonsSource;
        }

        if (totalSumPhotons <= 0)    // nothing to do
          continue;

        RandomSamplerFromPDF randomSource(sourcePDF);

        // no. photons to trace, taking into account the
        // min/max per bin and the extra ray-tracing factor
        const unsigned int tmpRayTrace = (unsigned int)(fExtraRayTraceFactor*((unsigned int)totalSumPhotons + 1));
        const unsigned int nRayTrace = std::min(fMaxNRayTrace, std::max(fMinNRayTrace, tmpRayTrace));

        totalNRayTracedPhotonsPerBin += nRayTrace;
        countRTPhotons += nRayTrace;
        const double nBunchPhotons = totalSumPhotons / nRayTrace;

        map<int, RandomSamplerFromPDF*> cherenkovProductionDistance;
        ScopeGuard cherProdDistGuard(cherenkovProductionDistance);

        for (unsigned int iSample = 0; iSample < nRayTrace; ++iSample) {

          // sample light source and wavelength
          const fevt::FdConstants::LightSource lightSource =
            fevt::FdConstants::LightSource(randomSource.shoot(fRandomEngine->GetEngine()));
          const int iwl = isLaser ? 0 : int(randomWavelengthSource[lightSource]->shoot(fRandomEngine->GetEngine()));

          EPhotonSource photonSource = eUnkown;
          double wavelength = 0;
          switch (lightSource) {
          case fevt::FdConstants::eTotal:
            ERROR("DONT KNOW WHAT TO DO");
            return eFailure;
            break;
          case fevt::FdConstants::eFluorTotal:
          case fevt::FdConstants::eFluorDirect:
          case fevt::FdConstants::eFluorRayleighScattered:
          case fevt::FdConstants::eFluorMieScattered:
          case fevt::FdConstants::eFluorMultScattered:
            photonSource = eFluorescence;
            wavelength = fUseOnlyReferenceWavelength ? refWl : WLengthFluo[iwl];
            break;
          case fevt::FdConstants::eCherDirect:
            photonSource = eDirectCherenkov;
            wavelength = fUseOnlyReferenceWavelength ? refWl : WLengthCkov[iwl];
            break;
          case fevt::FdConstants::eCherRayleighScattered:
          case fevt::FdConstants::eCherMieScattered:
          case fevt::FdConstants::ePrimaryCherDirect:
          case fevt::FdConstants::ePrimaryCherRayleighScattered:
          case fevt::FdConstants::ePrimaryCherMieScattered:
            {
              photonSource = eScatteredCherenkov;
              wavelength = fUseOnlyReferenceWavelength ? refWl : WLengthCkov[iwl];
              if (!cherenkovProductionDistance.count(iwl)) {
                TabulatedFunction& cherenkovBeamProduction = simShower.GetCherenkovBeamProductionPhotons(iwl);
                TabulatedFunction cherenkovProductionBeamLDF;
                double totalBeamLDF = 0; 
                for (unsigned int iTraceBeam = 0; iTraceBeam < cherenkovBeamProduction.GetNPoints(); ++iTraceBeam) {
                  const double distanceBeam = cherenkovBeamProduction[iTraceBeam].X();
                  if (distanceBeam > distancePhoton)
                    break;
                  const double nProd = cherenkovBeamProduction[iTraceBeam].Y();
                  cherenkovProductionBeamLDF.PushBack(distanceBeam, nProd);
                  totalBeamLDF += nProd;
                }
                cherenkovProductionDistance[iwl] = 0;             
                if (totalBeamLDF>0) {
                  cherenkovProductionDistance[iwl] =
                    new RandomSamplerFromPDF(cherenkovProductionBeamLDF, RandomSamplerFromPDF::eLinear);
                }
              }
              break;
            }
          case fevt::FdConstants::eLaserRayleighScattered:
          case fevt::FdConstants::eLaserMieScattered:
            photonSource = eLaser;
            wavelength = laserWavelength;
            break;
          case fevt::FdConstants::eBackground:
            ERROR("Background light photons are not implemented");
            return eFailure;
            break;
          default:
            break;
          }

          if (wavelength == 0) {
            ERROR("No wavelength selected");
            return eFailure;
          }

          // Generating random point on the diaphragm
          // Uniform phi
          const double diaPh = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, kTwoPi);
          const double diaR = rDia * sqrt(RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 1.0));
          const double xDia = diaR * cos(diaPh);
          const double yDia = diaR * sin(diaPh);
          const Point pIn(xDia, yDia, 0.0, telCS);

          Point pointOnShower = showerCore + showerAxis * distancePhoton;
          Vector photonDir = pointOnShower - telescopePos;
          double timeShift = 0; // time shift with respect to 1D flight path

          if (!isLaser && (fFluorescenceLDF != eOff ||
                           (fDirectCherenkovLDF != eOff) ||
                           (fScatteredCherenkovLDF != eOff))) {

            const double distToShower = (pointOnShower-telescopePos).GetMag();
            const double height = wgs84.PointToLatitudeLongitudeHeight(pointOnShower).get<2>();
            const double Xvert = depthProfile.Y(std::max(depthProfile.MinX(), std::min(height, depthProfile.MaxX())));
            const double Xslant = flatEarth ? Xvert/absCosZenith : slantDepthVsDistance->Y(std::max(atmDistanceMin, std::min(distancePhoton, atmDistanceMax)));

            if (Xslant > X1) { // no LDF before the first interaction

              const double age = ShowerAge(Xslant, showerXmax);
              // const double rm = moliereR;

              switch(photonSource) {

              case eFluorescence:
                {
                  switch(fFluorescenceLDF) {
                  case eNKG: pointOnShower += LateralDistributionNKG(showerCS, age, moliereR); break;
                  case eGora: pointOnShower += LateralDistributionGora(showerCS, age, moliereR); break;
                  default: break;
                  }
                }
                break;

              case eDirectCherenkov:
                {
                  switch(fDirectCherenkovLDF) {
                  case eNKG: pointOnShower += LateralDistributionNKG(showerCS, age, moliereR); break;
                  case eGora: pointOnShower += LateralDistributionGora(showerCS, age, moliereR); break;
                  default: break;
                  }
                }
                break;

              case eScatteredCherenkov:
                {
                  switch(fScatteredCherenkovLDF) {
                  case eNKG: pointOnShower += LateralDistributionNKG(showerCS, age, moliereR); break;
                  case eGora: pointOnShower += LateralDistributionGora(showerCS, age, moliereR); break;
                  case eCherenkov:
                  case eCherenkovAtAxis:
                    {
                      if (cherenkovProductionDistance[iwl]) {
                        double emissionDistance = 0;
                        Point pointOfCkovEmission;
                        Vector directionOfCkovEmission;
                        pointOnShower +=
                          LateralDistributionScatteredCherenkov(showerCS, showerCore, showerAxis,
                                                                distancePhoton,
                                                                showerXmax,
                                                                moliereR,
                                                                emissionDistance,
                                                                pointOfCkovEmission, directionOfCkovEmission,
                                                                *cherenkovProductionDistance[iwl],
                                                                depthProfile,
                                                                slantDepthVsDistance,
                                                                absCosZenith,
                                                                coreZ,
                                                                flatEarth,
                                                                atmDistanceMin,
                                                                atmDistanceMax);
                        /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                         * RU, Mon Feb  1 13:57:32 EST 2010
                         *
                         * the time shift by the longer flight path of the direct cherenkov
                         * photons in the Cherenkov cone vs. the 1 dimensional propagation is not taken
                         * into account. It is of the order of ~<100ns, and thus similar to the
                         * effect caused by the shower front curvature, which is also neglected
                         * here.
                         *
                         */
                      }
                    }
                    break;
                  default: break;
                  }
                }
                break;
              default: break;
              }

              photonDir = pointOnShower - telescopePos; // altered photon incident direction

              const double dist3D = photonDir.GetMag();    // real length of photon flight path
              timeShift += (dist3D - distToShower) / kSpeedOfLight; // distance difference -> time difference

#ifdef DEBUG_3D
              info << " p " << Pressure/bar
                   << " rm: " << moliereR/m
                   << " age: " << age
                   << " T: " << Temperature
                   << " Height: " << height/m
                   << " X: " << Xvert/g*cm*cm
                //                 << " tel: " << telescopePos.GetZ (localCS)
                   << " dist: " << distToShower/m
                   << " anlge: " << angleToCore/deg
                   << " shift: " << timeShift/ns
                   << endl;
#endif
            } // if Xslant>X1
            else {
              // before first interaction -> no LDF (!)
            }
            
          } // 3d shower
          
          photonDir.Normalize();
          const double projectedDiaphragmArea = diaphragmArea * photonDir.GetCosTheta(telCS);
          const double weight = nBunchPhotons * projectedDiaphragmArea; // NOW IN LIGHTATDIAPHRAGMSIMULATOR
          
          /*
          // don't track photons with weight<1
          if (weight<1) {
            if (RandFlat::shoot(&fRandomEngine->GetEngine(), 0., 1.))

          }
          */      

          /*
            This can happen if photonDir.GetCosTheta(telCS) gets negative due to
            large sample LDF distance in shower disc. Sometimes photons are then
            emitted "behind" the telescope aperture ...
          */
          if (weight <= 0)
            continue;

          totalWeightPhotons += weight;
          utl::Photon photonIn(pIn, -photonDir, wavelength, weight, lightSource);
          photonIn.SetTime(TimeInterval(traceTime + timeShift)); // time at diaphragm
          telSim.AddPhoton(photonIn);
          
          if (fMultipleScattering) {
            vector<Photon> photonVector;
            fMultipleScatterer->AddPhotons(pointOnShower, photonIn,
                                           telescopePos, photonVector);
            for (unsigned int iPh = 0; iPh < photonVector.size(); ++iPh)
              telSim.AddPhoton(photonVector[iPh]);
          }

        } // loop over raytrace samples

        rayTracedPhotonTrace[iTrace] += totalNRayTracedPhotonsPerBin;
      } // loop over time bin

      info << ", generated " << setw(6) << countRTPhotons << " photons"
           << ", total weight=" << totalWeightPhotons;
      INFO(info);

    } // End loop over Telescopes

  } // End loop over Eyes

  return eSuccess;

} // end of Run


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


#include <TCanvas.h>
#include <TPolyMarker3D.h>
#include <TPolyLine3D.h>
#include <TH2D.h>

TH1D*
ShowerPhotonGenerator::ToCumu(TH1D* h, double total)
{
  ostringstream newname;
  newname << h->GetName() << "_cumulative";
  TH1D* const o = new TH1D(newname.str().c_str(), newname.str().c_str(),
                           h->GetNbinsX(), h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
  o->SetLineColor(h->GetLineColor());
  o->SetLineStyle(h->GetLineStyle());
  o->SetLineWidth(h->GetLineWidth());
  if (total <= 0)
    total = 1.0;
  double sum = 0;
  for (int i = 1; i <= h->GetNbinsX(); ++i) {
    sum += h->GetBinContent(i);
    o->SetBinContent(i, sum/total);
  }
  return o;
}


// ---------------------------------------------------------------------------------------
void
ShowerPhotonGenerator::PlotLDF(const ShowerSimData& simShower,
                               const TabulatedFunction& cherenkovProductionBeam,
                               const double age)
{
  const bool flatEarth = false;

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

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

  const ProfileResult& depthProfile = atmo.EvaluateDepthVsHeight();
  const ProfileResult& temperatureProfile = atmo.EvaluateTemperatureVsHeight();
  const double tempHeightMin = temperatureProfile.MinX();
  const double tempHeightMax = temperatureProfile.MaxX();
  const ProfileResult& distanceVsSlantDepth = atmo.EvaluateDistanceVsSlantDepth();
  const ProfileResult* const slantDepthVsDistance = &atmo.EvaluateSlantDepthVsDistance();
  const double atmDistanceMin = slantDepthVsDistance->MinX();
  const double atmDistanceMax = slantDepthVsDistance->MaxX();
  const ProfileResult& heightVsDistance = atmo.EvaluateHeightVsDistance();

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

  const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
  const double cosZenith = (-simShower.GetDirection()).GetCosTheta(localCS); // just for flat earth
  const double absCosZenith = std::fabs(cosZenith); // just for flat earth
  const double coreZ = wgs84.PointToLatitudeLongitudeHeight(showerCore).get<2>();
  const double showerXmax = simShower.GetGHParameters().GetXMax();
  const double X1 = simShower.GetXFirst();
  const double distanceX1 = atmo.EvaluateDistanceVsSlantDepth().Y(X1);
  const Point& pointX1 = showerCore + showerAxis*distanceX1;

  const double Xslant = 2. * age * showerXmax / (3 - age);
  const double distance = distanceVsSlantDepth.Y(Xslant);
  const double height = heightVsDistance.Y(distance);
  const double Xvert = depthProfile.Y(height);
  const Point& pointOnShower = showerCore + showerAxis*distance;

  const double Temperature = temperatureProfile.Y(std::max(std::min(height, tempHeightMax), tempHeightMin));
  const double Pressure = Xvert * kgEarth;
  const double rm = MoliereRadius(Temperature, Pressure,
                                  (-showerAxis).GetCosTheta(LocalCoordinateSystem::Create(UTMPoint(pointOnShower, wgs84))) );

  TabulatedFunction cherenkovProductionBeamLDF;
  for (unsigned int iTrace = 0; iTrace < cherenkovProductionBeam.GetNPoints(); ++iTrace) {
    const double distanceBin = cherenkovProductionBeam[iTrace].X();
    if (distanceBin > distance)
      break;
    cherenkovProductionBeamLDF.PushBack(cherenkovProductionBeam[iTrace].X(), cherenkovProductionBeam[iTrace].Y());
  }
  RandomSamplerFromPDF* const cherenkovProductionDistance =
    new RandomSamplerFromPDF(cherenkovProductionBeamLDF, RandomSamplerFromPDF::eLinear);

  // RU: for making LDF-pictures
  cout << " PlotLDF: age=" << age
       << " Xslant=" << Xslant/g*cm2
       << " Xvert=" << Xvert/g*cm2
       << " distance=" << distance
       << " height=" << height
       << " rm=" << rm
       << " bins=" << cherenkovProductionBeamLDF.GetNPoints()
       << endl;

  ostringstream hFluoLDFName, hCkovLDFName;
  hFluoLDFName << "hFluoLDF_" << fixed << setprecision(2) << age;
  hCkovLDFName << "hCkovLDF_" << fixed << setprecision(2) << age;
  const int nBins = 100;
  const double minR = 0;
  const double maxR = 1600;
  TH1D* const hFluoLDF = new TH1D(hFluoLDFName.str().c_str(), hFluoLDFName.str().c_str(), nBins, minR, maxR);
  TH1D* const hCkovLDF = new TH1D(hCkovLDFName.str().c_str(), hCkovLDFName.str().c_str(), nBins, minR, maxR);
  TCanvas* const c3D = new TCanvas("c3D", "c3D", 500, 1500);
  c3D->SetFillColor(kWhite);
  /*
    TH2D* frame = new TH2D("frame", "rame",
    10, -maxR, maxR, 10, distanceX1, distance, 0);
    frame->Draw();
    TH2D* projXZ = new TH2D("projXZ", "projection",
    100, -maxR, maxR, 500, distanceX1, distance, 0);
  */
  TPolyMarker3D* const p1 = new TPolyMarker3D(2);
  p1->SetPoint(0, pointOnShower.GetX(showerCS), pointOnShower.GetY(showerCS), pointOnShower.GetZ(showerCS));
  p1->SetPoint(1, pointX1.GetX(showerCS), pointX1.GetY(showerCS), pointX1.GetZ(showerCS));
  p1->SetMarkerStyle(20);
  p1->SetMarkerColor(kRed);
  p1->Draw();
  const int n_samples = 100000;
  for (int i = 0; i < n_samples; ++i) {
    const double fluoLDF = (fFluorescenceLDF == eGora ?
                            LateralDistributionGora(showerCS, age, rm).GetMag() :
                            LateralDistributionNKG(showerCS, age, rm).GetMag()) ;
    double emissionDistance = 0;
    Point pointOfCkovEmission;
    Vector directionOfCkovEmission;
    Vector ckovR = LateralDistributionScatteredCherenkov(showerCS, showerCore, showerAxis,
                                                         distance,
                                                         showerXmax,
                                                         rm,
                                                         emissionDistance,
                                                         pointOfCkovEmission,
                                                         directionOfCkovEmission,
                                                         *cherenkovProductionDistance,
                                                         depthProfile,
                                                         slantDepthVsDistance,
                                                         absCosZenith,
                                                         coreZ,
                                                         flatEarth,
                                                         atmDistanceMin,
                                                         atmDistanceMax);
    TPolyLine3D* const l3D = new TPolyLine3D(2);
    l3D->SetPoint(0, pointOfCkovEmission.GetX(showerCS), pointOfCkovEmission.GetY(showerCS), pointOfCkovEmission.GetZ(showerCS));
    l3D->SetPoint(1, (pointOnShower+ckovR).GetX(showerCS), (pointOnShower+ckovR).GetY(showerCS), (pointOnShower+ckovR).GetZ(showerCS));
    l3D->SetLineStyle(1);
    l3D->SetLineWidth(1);
    l3D->SetLineColor(kBlue);
    l3D->Draw();
    const double ckovLDF = ckovR.GetMag();
    hFluoLDF->Fill(fluoLDF/m);
    hCkovLDF->Fill(ckovLDF/m);
  } // end n-photons

  TH1D* const hFluoLDFc = ToCumu(hFluoLDF, n_samples);
  TH1D* const hCkovLDFc = ToCumu(hCkovLDF, n_samples);

  ostringstream hLDFPrint, hLDFcumu;
  hLDFPrint << "LDF_" << fixed << setfill('0') << setw(3) << setprecision(3) <<  age << ".eps";
  hLDFcumu << "LDF_" << fixed << setfill('0') << setw(3) << setprecision(3) <<  age << ".cumulative.eps";
  gROOT->SetStyle("Plain");
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  TCanvas* const c = new TCanvas("ldf", "ldf");
  TLegend* const l = new TLegend(0.5, 0.6, 0.85, 0.85);
  l->SetFillColor(0);
  l->SetBorderSize(0);
  l->SetTextFont(42);
  l->SetTextSize(0.045);
  l->AddEntry(hFluoLDF, "Fluorescence", "l");
  l->AddEntry(hCkovLDF, "Cherenkov", "l");

  double integral = hFluoLDF->Integral();
  if (integral)
    hFluoLDF->Scale(1. / integral);
  integral = hCkovLDF->Integral();
  if (integral)
    hCkovLDF->Scale(1. / integral);

  hFluoLDF->SetLineColor(2);
  hCkovLDF->SetLineColor(4);
  hFluoLDF->SetXTitle("Distance   [m]");
  hFluoLDF->SetYTitle("R dP/dR");
  hFluoLDF->GetYaxis()->SetTitleOffset(1.15);
  hFluoLDF->Draw("l");
  hCkovLDF->Draw("lsame");
  l->Draw();
  c->Print(hLDFPrint.str().c_str());

  TCanvas* const c2 = new TCanvas("ldfcumu", "ldfcumu");
  TLegend* const l2 = new TLegend(0.5, 0.2, 0.85, 0.6);
  l2->SetFillColor(0);
  l2->SetBorderSize(0);
  l2->SetTextFont(42);
  l2->SetTextSize(0.045);
  l2->AddEntry(hFluoLDF, "Fluorescence", "l");
  l2->AddEntry(hCkovLDF, "Cherenkov", "l");

  hFluoLDFc->SetLineColor(2);
  hCkovLDFc->SetLineColor(4);
  hFluoLDFc->SetXTitle("Distance   [m]");
  hFluoLDFc->SetYTitle("R dP/dR");
  hFluoLDFc->GetYaxis()->SetTitleOffset(1.15);
  hFluoLDFc->Draw("l");
  hCkovLDFc->Draw("lsame");
  l2->Draw();
  c2->Print(hLDFcumu.str().c_str());

  ostringstream h3DPrint;
  h3DPrint << "LDF_" << fixed << setfill('0') << setw(3) << setprecision(3) <<  age << ".3d.root";
  //c3D.ShoweAxis3D();
  c3D->Print(h3DPrint.str().c_str());

  delete cherenkovProductionDistance;
}


//==========================================================================
//==========================================================================
Vector
ShowerPhotonGenerator::LateralDistributionNKG(const CoordinateSystemPtr& shwCS,
                                              const double s_age,
                                              const double r_moliere)
{
  const double fMax = 0.95;
  const double r = r_moliere * NKGFunction(s_age, fMax);
  const double phi = RandFlat::shoot(&fRandomEngine->GetEngine(), 0., kTwoPi);

  return Vector(r, phi, 0., shwCS, Vector::kCylindrical);
}


//==========================================================================
//==========================================================================
Vector
ShowerPhotonGenerator::LateralDistributionGora(const CoordinateSystemPtr& shwCS,
                                               const double s_age,
                                               const double r_moliere)
{
  const double r = r_moliere * GoraFunction(s_age);
  const double phi = RandFlat::shoot(&fRandomEngine->GetEngine(), 0., kTwoPi);

  return Vector(r, phi, 0., shwCS, Vector::kCylindrical);
}


//==========================================================================
//==========================================================================
Vector
ShowerPhotonGenerator::LateralDistributionScatteredCherenkov(const utl::CoordinateSystemPtr& showerCS,
                                                             const utl::Point& core,
                                                             const utl::Vector& axis,
                                                             double distanceCoreToPointOnShower,
                                                             double Xmax,
                                                             double rMoliere,
                                                             double& emissionDistance,
                                                             utl::Point& pointOfEmission,
                                                             utl::Vector& directionOfEmission,
                                                             const utl::VRandomSampler& cherenkovProduction,
                                                             const ProfileResult& depthProfile,
                                                             const ProfileResult* slantDepthVsDistance,
                                                             const double absCosZenith,
                                                             const double coreZ,
                                                             const bool flatEarth,
                                                             const double atmDistanceMin,
                                                             const double atmDistanceMax)
{
  // --------- point of cherenkov emission at shower axis ---------
  emissionDistance = cherenkovProduction.shoot(fRandomEngine->GetEngine());

  const bool emissionIn = emissionDistance > atmDistanceMin && emissionDistance < atmDistanceMax;
  const double emissionDepthVert =
    emissionIn ? depthProfile.Y(std::max(depthProfile.MinX(), std::min(depthProfile.MaxX(), coreZ-emissionDistance*absCosZenith))) : 0;
  const double emissionDepth =
    flatEarth ? emissionDepthVert / absCosZenith :
      (emissionIn ? slantDepthVsDistance->Y(emissionDistance) : 0);
  const double emissionAge = ShowerAge(emissionDepth, Xmax);

  pointOfEmission = core + emissionDistance*axis + LateralDistributionGora(showerCS, emissionAge, rMoliere);

  // --------- Cherenkov emission angle at shower axis  ---------
  const double maxIntegrationAngle = 90.*utl::degree;
  const double deltaAngle = 0.25*utl::degree;

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

  const double rndmCDFVal = RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);
  const double cdfMax =
    atmo.AngularCherenkovCDF(maxIntegrationAngle, emissionDepthVert, emissionAge);

  // find inverse of CDF at rndmCDFVal using bisection method
  double leftAngle = 0.;
  double rightAngle = maxIntegrationAngle;
  double cdfLeft =
    atmo.AngularCherenkovCDF(leftAngle, emissionDepthVert, emissionAge)/cdfMax - rndmCDFVal;

  while (fabs(leftAngle - rightAngle) > 2*deltaAngle) {

    const double midAngle = 0.5 * (leftAngle + rightAngle);
    const double cdfMiddle =
      atmo.AngularCherenkovCDF(midAngle, emissionDepthVert, emissionAge)/cdfMax - rndmCDFVal;

    if (cdfLeft*cdfMiddle < 0)
      rightAngle = midAngle;
    else {
      cdfLeft =
        atmo.AngularCherenkovCDF(midAngle, emissionDepthVert, emissionAge)/cdfMax - rndmCDFVal;
      leftAngle = midAngle;
    }
  }

  const double cdfRight =
    atmo.AngularCherenkovCDF(rightAngle, emissionDepthVert, emissionAge)/cdfMax - rndmCDFVal;
  const double electronAngleTheta = leftAngle - (rightAngle-leftAngle)/(cdfRight-cdfLeft) * (cdfLeft);

  const double electronAnglePhi = RandFlat::shoot(&fRandomEngine->GetEngine(), 0., kTwoPi);

  const CoordinateSystemPtr electronCS =
    CoordinateSystem::RotationX(electronAngleTheta, CoordinateSystem::RotationZ(electronAnglePhi, showerCS));

  // sample from Cherenkov cone
  const double cherenkovAnglePhi = RandFlat::shoot(&fRandomEngine->GetEngine(), 0., kTwoPi);
  // Cherenkov angle for beta=1
  const double nRefrac = utl::RefractionIndex::LorentzLorentz(emissionDepthVert);
  const double cosCherenkovTheta = 1. / nRefrac;
  const double sinCherenkovTheta = sqrt(1.-cosCherenkovTheta*cosCherenkovTheta);
  /*
    //#warning turned off cherenkov cone
    const double cherenkovAnglePhi = 0;
    const double cosCherenkovTheta = 1.;
    const double sinCherenkovTheta = 0;
  */

  directionOfEmission = Vector(-sinCherenkovTheta*sin(cherenkovAnglePhi),
                               -sinCherenkovTheta*cos(cherenkovAnglePhi),
                               -cosCherenkovTheta,
                               electronCS);

  // const double cosPhotonAngle = CosAngle(directionOfEmission, axis); // not needed any more

  const double photonTravelDistance = (distanceCoreToPointOnShower-emissionDistance) / (directionOfEmission*axis); //cosPhotonAngle;

  return Vector(pointOfEmission.GetX(showerCS) + directionOfEmission.GetX(showerCS) * photonTravelDistance,
                pointOfEmission.GetY(showerCS) + directionOfEmission.GetY(showerCS) * photonTravelDistance,
                0,
                showerCS);
}


double
ShowerPhotonGenerator::NKGFunction(double s, const double fMax)
{
  const double closeToZero = 1.e-20;
  const double sMax = 1.999;

  static double lastFmax = -1;
  static double xMaxArr[200];

  // initialize xMax array on startup or if fMax changed
  if (fMax != lastFmax) {
    lastFmax = fMax;
    double dX = 0.001;
    for (int j = 0; j < 200; ++j) {
      const double a = (0.5+j)/100;
      const double b = 4.5 - 2*a;
      double frac = 0;
      double xxMax = closeToZero;
      while (frac < fMax) {
        xxMax += dX;
        const double arg = 1/(1+1/xxMax);
        frac = IncompleteBeta(a, b, arg);
      }
      xMaxArr[j] = xxMax;
      // incrase dX
      dX = xxMax/100;
    }
  }

  // random generation

  if (s >= sMax)
    s = sMax;

  const double xMin = closeToZero;
  const int ageBin = int(s*100);
  const double xMax = xMaxArr[ageBin];

  double rdm = 1;
  double Y = 0;
  double x;
  while (rdm > Y) {

    if (s > 1) {

      // get random number ~ (1.+x)^(s-4.5) between Xmin and Xmax

      const double g = s - 4.5;
      const double g1 = g + 1;

      rdm = RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);

      const double Fa = pow(xMin+1, g1);
      const double Fb = pow(xMax+1, g1);
      const double y = Fa + (Fb-Fa)*rdm;
      x = pow(y, 1/g1);

      // hit or miss for additional x*(x)^(s-2)

      x -= 1;
      Y = pow(x, s-1);
      const double Ymax = pow(xMax, s-1);
      rdm = Ymax * RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);

    } else {

      // get random number ~ x*(x)^(s-2.) between Xmin and Xmax

      const double g = s - 1;
      const double g1 = g + 1;

      rdm = RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);

      if (fabs(g1) > closeToZero) {
        const double Fa = pow(xMin, g1);
        const double Fb = pow(xMax, g1);
        const double y = Fa + (Fb - Fa)*rdm;
        x = pow(y, 1/g1);
      } else {
        const double lMin = log10(xMin);
        const double lMax = log10(xMax);
        x = pow(10, lMin + rdm*(lMax - lMin));
      }

      // hit or miss for additional (1+x)^(s-4.5)

      Y = pow(1+x, s - 4.5);
      rdm = RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);
    }

  }

  return x;
}


//==========================================================================
//==========================================================================
inline double
ShowerPhotonGenerator::GoraFunction(const double s)
{

  double rdm = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 1.0);

  // avoid div by zero
  if (rdm < std::numeric_limits<double>::min())
    rdm = std::numeric_limits<double>::min();

  // r = F(rdm)^(-1)
  return utl::InverseGoraCDF(rdm, s);
}


//==========================================================================
//==========================================================================
ScopeGuard::ScopeGuard(std::map<int, utl::RandomSamplerFromPDF*>& theMap)
  : fMap(theMap)
{
}

//==========================================================================
ScopeGuard::~ScopeGuard()
{
  for (map<int, RandomSamplerFromPDF*>::iterator iMap = fMap.begin(),
       end = fMap.end();
       iMap != end; ++iMap) {
    delete iMap->second;
  }
}

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: