LateralLightCalculator.cc

Go to the documentation of this file.

#include "LateralLightCalculator.h"
#include "TelescopeData.h"
//#warning MU move OpticalHalo to FDetector
#include "../FdProfileReconstructorKG/OpticalHalo.h"

#include <fwk/CentralConfig.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Point.h>
#include <utl/Integrator.h>
#include <utl/TabulatedFunction.h>
#include <det/Detector.h>
#include <fdet/Telescope.h>
#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>
#include <atm/AttenuationResult.h>
#include <atm/ScatteringResult.h>
#include <atm/InclinedAtmosphericProfile.h>

#include <string>
#include <sstream>
using namespace std;
using namespace utl;
using namespace fdet;
using namespace fevt;
using namespace atm;
using namespace FdEnergyDepositFinderKG;


class ScattCherFunctor {

public:
  ScattCherFunctor(const double distanceToTelescope,
                   const double distanceAlongShower,
                   const double showerAge,
                   const double verticalDepth,
                   const TabulatedFunction& arcLengthFunction) :
    fDistanceToTelescope(distanceToTelescope),
    fDistanceAlongShower(distanceAlongShower),
    fShowerAge(showerAge),
    fVerticalDepth(verticalDepth),
    fArcLengthFunction(arcLengthFunction),
    fSmallAngleApprox(false)
  { }

  double
  operator()(const double zeta)
    const
  {
    const double distanceRatio = fDistanceToTelescope / fDistanceAlongShower;
    const auto& atmo = det::Detector::GetInstance().GetAtmosphere();

    double pdf;
    double jacobian;

    if (fSmallAngleApprox) {

      const double theta = zeta * distanceRatio;
      jacobian = distanceRatio;
      pdf = atmo.AngularCherenkovPDF(theta, fVerticalDepth, fShowerAge);

    } else {

      const double cosZeta = cos(zeta);
      const double sinZeta = sin(zeta);
      const double tanZeta = sinZeta / cosZeta;

      const double cosZeta2 = cosZeta * cosZeta;
      const double sinZeta2 = sinZeta * sinZeta;
      const double theta = atan(tanZeta * distanceRatio);
      jacobian = 1 / (cosZeta2/distanceRatio + distanceRatio*sinZeta2);
      pdf = atmo.AngularCherenkovPDF(theta, fVerticalDepth, fShowerAge);

    }

    const double arcLength = fArcLengthFunction.Y(zeta);

    return pdf * jacobian * arcLength;
  }

private:
  double fDistanceToTelescope = 0;
  double fDistanceAlongShower = 0;
  double fShowerAge = 0;
  double fVerticalDepth = 0;
  const TabulatedFunction& fArcLengthFunction;
  bool fSmallAngleApprox = false;

};


class FluorescenceFunctor {

public:
  FluorescenceFunctor(const double showerAge,
                      const double distance,
                      const double moliereRadius,
                      const TabulatedFunction& arcLengthFunction):
    fShowerAge(showerAge),
    fDistance(distance),
    fMoliereRadius(moliereRadius),
    fArcLengthFunction(arcLengthFunction),
    fSmallAngleApprox(false)
  { }

  double
  operator()(const double zeta)
    const
  {
    const double age = fShowerAge;
    const double d = fDistance;
    const double rM = fMoliereRadius;
    const double distanceRatio = d / rM;

    double pdf;
    double jacobian;
    if (fSmallAngleApprox) {
      jacobian = distanceRatio;
      pdf = GoraPDF(zeta * distanceRatio, age);
    } else {
      const double tanZeta = tan(zeta);
      const double rStar = distanceRatio * tanZeta;
      jacobian = distanceRatio*(1 + Sqr(tanZeta));
      pdf = GoraPDF(rStar, age);
    }

    const double arcLength = fArcLengthFunction.Y(zeta);

    return pdf * jacobian * arcLength;
  }

private:
  double fShowerAge = 0;
  double fDistance = 0;
  double fMoliereRadius = 0;
  const TabulatedFunction& fArcLengthFunction;
  bool fSmallAngleApprox = false;

};


LateralLightCalculator::LateralLightCalculator() :
  fPrecision(1e-4),
  fOpticalHalo(nullptr)
{
  Init();
}


LateralLightCalculator::~LateralLightCalculator()
{
  delete fOpticalHalo;
}


void
LateralLightCalculator::Init()
{
  auto topBranch =
    fwk::CentralConfig::GetInstance()->GetTopBranch("FdEnergyDepositFinder");
  auto lldBranch = topBranch.GetChild("lateralLight");

  ostringstream info;

  string lldOption;
  lldBranch.GetChild("method").GetData(lldOption);
  fMethod =  (lldOption == "eCircle") ?
    eCircle :
      (lldOption == "eBinByBin") ?
        eBinByBin : ePixelIntegrated;
  info << "\n\t method: " << lldOption;

  bool lldIsOn;
  lldBranch.GetChild("fluorescence").GetData(lldIsOn);
  fLLDIsOn[FdConstants::eFluorDirect] = lldIsOn;
  info << "\n\t lateral fluorescence is "
       << (lldIsOn ? "on" : "off");

  lldBranch.GetChild("directCherenkov").GetData(lldIsOn);
  fLLDIsOn[FdConstants::eCherDirect] = lldIsOn;
  info << "\n\t lateral direct Cherenkov is "
       << (lldIsOn ? "on" : "off");

  lldBranch.GetChild("directCherenkovCorrectedGora").GetData(fDirCherCorrGora);
  info << "\n\t lateral direct Cherenkov correction for Gora is "
       << (fDirCherCorrGora ? "on" : "off");

  lldBranch.GetChild("scatteredCherenkov").GetData(lldIsOn);
  fLLDIsOn[FdConstants::eCherMieScattered] = lldIsOn;
  fLLDIsOn[FdConstants::eCherRayleighScattered] = lldIsOn;
  info << "\n\t lateral scattered Cherenkov is "
       << (lldIsOn ? "on" : "off");

  if (lldBranch.GetChild("lateralWidthCorrection").Get<bool>()) {
    if (fMethod == ePixelIntegrated) {
      info << "\n\t lateralWidthCorrection in PixelIntegrated mode is "
              "calculated in the same way as for eCircle option";
    } else if (fMethod != eCircle) {
      throw NonExistentComponentException("lateralWidthCorrection can currently"
                                          "only be used with eCircle option...");
    }

    using namespace FdProfileReconstructorKG;
    fOpticalHalo = new OpticalHalo(OpticalHalo::eSpotGroup2);
  }
  info << "\n\t lateral width correction is "
       << (fOpticalHalo ? "on" : "off");

  INFO(info);
}


double
LateralLightCalculator::GetLightFraction(fevt::FdConstants::LightSource source,
                                         const TelescopeDataBin& emissionBin,
                                         const TelescopeDataBin& detectionBin,
                                         const fdet::Telescope& telescope,
                                         const double xMax,
                                         const double cosTheta,
                                         const double zeta)
  const
{
  const auto sourceIter = fLLDIsOn.find(source);

  if (sourceIter != fLLDIsOn.end()) {

    if (sourceIter->second) {
      if (source == FdConstants::eFluorDirect)
        return FluoLightFraction(detectionBin, telescope, xMax, cosTheta, zeta);
      else if (source == FdConstants::eCherDirect)
        return DirCherLightFraction(detectionBin, telescope, xMax, cosTheta, zeta);
      else if (source == FdConstants::eCherMieScattered ||
               source == FdConstants::eCherRayleighScattered)
        return ScattCherLightFraction(emissionBin, detectionBin, telescope, xMax, zeta);
    } else {
      for (unsigned int i = 0; i < emissionBin.GetZetaPixels().size(); ++i)
        emissionBin.GetZetaPixels()[i].fLightFraction = emissionBin.GetZetaPixels()[i].IsInside(0, 0);
      return 1;  // option for no ldd
    }
  }

  ostringstream errMsg;
  errMsg << " source " << source << " not implemented";
  ERROR(errMsg);
  throw NonExistentComponentException(errMsg.str());
}


double
LateralLightCalculator::FluoLightFraction(const TelescopeDataBin& telData,
                                          const fdet::Telescope& detTel,
                                          const double xMax,
                                          const double cosTheta,
                                          const double zeta)
  const
{
  if (fMethod ==  eBinByBin) {
    if (telData.HasEfficiency(FdConstants::eFluorDirect))
      return telData.GetEfficiency(FdConstants::eFluorDirect);
    else {
      const string errMsg = "no bin-by-bin efficiency for fluorescence light";
      ERROR(errMsg);
      throw NonExistentComponentException(errMsg);
    }
  } else {
    const auto& atmo = det::Detector::GetInstance().GetAtmosphere();

    const auto& tempProfile = atmo.EvaluateTemperatureVsHeight();
    const auto& pressureProfile = atmo.EvaluatePressureVsHeight();

    const double height = telData.fHeight;
    const double xSlant = telData.GetMeanDepth();
    const double pressure = pressureProfile.Y(height);
    const double temperature = tempProfile.Y(height);
    const double age = ShowerAge(xSlant, xMax);
    const double rMoliere = MoliereRadius(temperature, pressure, cosTheta);

    const auto& meanPos = telData.fMeanDepthPoint;
    const auto& telescopePos = detTel.GetPosition();
    const double showerTelDist = (meanPos - telescopePos).GetMag();
    const double distanceToShowerAxis = tan(zeta) * showerTelDist;

    const double rStar = distanceToShowerAxis / rMoliere;

    if (fMethod ==  eCircle || telData.GetZetaPixels().empty()) {
      const double spotFraction = fOpticalHalo ?
        fOpticalHalo->GetHaloFraction(zeta, age, showerTelDist) : 1;
      return GoraCDF(rStar, age) * spotFraction;
    } else if (fMethod ==  ePixelIntegrated) {
      FluorescenceFunctor fluoFunc(age, showerTelDist, rMoliere,
                                   telData.GetArcLengthFunction());
      const Integrator<FluorescenceFunctor> integrator(fluoFunc, fPrecision);
      const double fraction = integrator.GetIntegral(0., telData.GetMaxZeta());

      for (unsigned int i = 0; i < telData.GetZetaPixels().size(); ++i) {
        if (telData.GetZetaPixels()[i].fLightFraction == -1) {
          FluorescenceFunctor fluoFuncPixel(age, showerTelDist, rMoliere,
                                            telData.GetZetaPixels()[i].GetArcLengthFunction());
          const Integrator<FluorescenceFunctor> integratorPixel(fluoFuncPixel, fPrecision);
          telData.GetZetaPixels()[i].fLightFraction = integratorPixel.GetIntegral(0., telData.GetMaxZeta())/fraction;
        }
      }
      const double spotFraction = fOpticalHalo ?
        fOpticalHalo->GetHaloFraction(zeta, age, showerTelDist) : 1;
      return fraction * spotFraction;
    }
  }

  return 1;
}


double
LateralLightCalculator::ScattCherLightFraction(const TelescopeDataBin& emissionBin,
                                               const TelescopeDataBin& detectionBin,
                                               const fdet::Telescope& detTel,
                                               const double xMax,
                                               const double zeta)
  const
{
  if (fMethod ==  eBinByBin) {

    int nEff = 0;
    double eff = 0;
    if (detectionBin.HasEfficiency(FdConstants::eCherMieScattered)) {
      eff += detectionBin.GetEfficiency(FdConstants::eCherMieScattered);
      ++nEff;
    }
    if (detectionBin.HasEfficiency(FdConstants::eCherRayleighScattered)) {
      eff += detectionBin.GetEfficiency(FdConstants::eCherRayleighScattered);
      ++nEff;
    }
    if (nEff == 0) {
      const string errMsg = "no bin-by-bin efficiency for fluorescence light";
      ERROR(errMsg);
      throw NonExistentComponentException(errMsg);
    } else
      return eff / nEff;

  } else {

    const double pathLength =
      (detectionBin.fMeanDepthPoint - emissionBin.fMeanDepthPoint).GetMag();

    if (pathLength < 1*m)
      return 1;

    const auto& meanPos = detectionBin.fMeanDepthPoint;
    const auto& telescopePos = detTel.GetPosition();
    const double showerTelDist = (meanPos - telescopePos).GetMag();
    const double distanceToShowerAxis = tan(zeta) * showerTelDist;

    const auto& atmo = det::Detector::GetInstance().GetAtmosphere();
    const auto& depthProfile = atmo.EvaluateDepthVsHeight();
    const double height = emissionBin.fHeight;
    const double xSlant = emissionBin.GetMeanDepth();
    const double xVert = depthProfile.Y(height);
    const double showerAge = utl::ShowerAge(xSlant, xMax);

    if (fMethod ==  eCircle || detectionBin.GetZetaPixels().empty()) {
      const double maxEmissionAngle = atan2(distanceToShowerAxis, pathLength);
      const double fraction =
        atmo.AngularCherenkovCDF(maxEmissionAngle, xVert, showerAge);
      const double spotFraction = fOpticalHalo ?
        fOpticalHalo->GetHaloFraction(zeta, showerAge, showerTelDist) : 1;
      return fraction * spotFraction;
    } else if (fMethod ==  ePixelIntegrated) {
      ScattCherFunctor cherFunc(showerTelDist, pathLength,
                                showerAge, xVert,
                                detectionBin.GetArcLengthFunction());
      const Integrator<ScattCherFunctor> integrator(cherFunc, fPrecision);
      const double fraction =
        integrator.GetIntegral(0., detectionBin.GetMaxZeta());
      return fraction;
    }
  }

  return 1;
}


double
LateralLightCalculator::DirCherLightFraction(const TelescopeDataBin& telData,
                                             const fdet::Telescope& detTel,
                                             const double xMax,
                                             const double cosTheta,
                                             const double zeta)
  const
{
  if (fMethod ==  eBinByBin) {
    if (telData.HasEfficiency(FdConstants::eCherDirect))
      return telData.GetEfficiency(FdConstants::eCherDirect);
    else {
      const string errMsg = "no bin-by-bin efficiency for direct Cherenkov light";
      ERROR(errMsg);
      throw NonExistentComponentException(errMsg);
    }
  } else {

    if (!fDirCherCorrGora) {

      // currently same lateral light distribution as fluorescence
      return FluoLightFraction(telData, detTel, xMax, cosTheta, zeta);

    } else {

      // corrections to Gora LDF distributions based on 3D simulations in CORSIKA
      // at low energies (10^15 - 10^17 eV) - parametrized by pol5 in viewing angle
      static const double p_VA[30][6] = {
        { 0.501418, 2.94216, -7.30234, 9.31207, -5.99745, 1.54418 },
        { -0.00333121, 4.5624, -9.71922, 11.7979, -7.69457, 2.05705 },
        { -0.00994375, 1.81732, 0.719025, -4.27035, 4.22282, -1.48334 },
        { -0.00107321, 0.503863, 5.99427, -13.9105, 12.9082, -4.49541 },
        { 0.000745283, 0.148745, 6.33147, -13.439, 12.1988, -4.23893 },
        { 0.000709519, 0.0885937, 5.10991, -8.54954, 6.0705, -1.71655 },
        { 0.000880354, 0.0530137, 4.34411, -5.90561, 3.24284, -0.732068 },
        { 0.000579476, 0.0499766, 3.56326, -3.08367, -0.234635, 0.708838 },
        { 0.000467793, 0.0390384, 3.02424, -1.45418, -1.85989, 1.25452 },
        { 0.000499658, 0.0371965, 2.43002, 0.258876, -3.57741, 1.85539 },
        { 0.00040399, 0.0208659, 2.08504, 1.23074, -4.54824, 2.21649 },
        { 8.86504e-05, 0.0440711, 1.35947, 3.24868, -6.701, 3.05502 },
        { -0.000288326, 0.055371, 0.907594, 4.31152, -7.62045, 3.35124 },
        { -0.000158639, 0.0429035, 0.800877, 4.33574, -7.43096, 3.2566 },
        { -0.000312911, 0.0631904, 0.162198, 6.2208, -9.62415, 4.1835 },
        { -0.000578045, 0.0752885, 0.0811873, 5.96277, -8.99537, 3.87796 },
        { -0.000428683, 0.0672003, -0.0811488, 6.14087, -9.08156, 3.9582 },
        { -0.000585311, 0.083422, -0.437892, 6.90878, -9.8617, 4.30713 },
        { -0.000603162, 0.0796437, -0.646415, 7.40142, -10.4691, 4.63569 },
        { -0.000975556, 0.116029, -0.977511, 7.81595, -10.8033, 4.84644 },
        { -0.00104516, 0.121701, -1.13216, 8.16602, -11.4698, 5.30535 },
        { -0.000692655, 0.0743631, -0.839957, 7.01192, -10.029, 4.76265 },
        { -0.000873929, 0.0974708, -1.08584, 7.59646, -11.055, 5.43051 },
        { -0.000278001, 0.0384464, -0.475365, 5.31638, -8.13554, 4.22311 },
        { -0.000373472, 0.0491568, -0.658973, 5.52961, -8.70434, 4.73939 },
        { -0.000784299, 0.0923149, -1.20043, 7.06114, -10.5462, 5.52973 },
        { -0.000272599, 0.0343134, -0.602254, 5.44353, -9.1348, 5.18153 },
        { -0.000153735, 0.0203647, -0.445935, 4.48833, -7.8483, 4.69983 },
        { -0.00365289, 0.175263, -2.03252, 9.50892, -14.4687, 7.745 },
        { -0.00219899, 0.110025, -1.4529, 7.59706, -12.3696, 6.98898 }
      };

      const auto& meanPos = telData.fMeanDepthPoint;
      const auto& telPos = detTel.GetPosition();
      const auto showerDir = Normalized(telData.fLastPoint - telData.fFirstPoint);
      const double va = Angle(telPos - meanPos, showerDir) / degree;
      const double* p = p_VA[int(round(max(0.,  min(29., floor(va)))))];
      const double x = FluoLightFraction(telData, detTel, xMax, cosTheta, zeta);
      return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*(p[4] + x*p[5]))));
    }
  }
}