CFMatrixCalculator.cc

Go to the documentation of this file.

#include "CFMatrixCalculator.h"
#include "LateralLightCalculator.h"
#include "TelescopeData.h"
#include "ZetaPixel.h"

#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>
#include <fevt/EyeHeader.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeRecData.h>
#include <fevt/Pixel.h>

#include <evt/ShowerFRecData.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Pixel.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 <fwk/CentralConfig.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/UTMPoint.h>
#include <utl/AugerException.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>

#include <fwk/LocalCoordinateSystem.h>

#include <sstream>
#include <string>
#include <algorithm>

using namespace fevt;
using namespace evt;
using namespace utl;
using namespace atm;
using namespace det;
using namespace std;
using namespace FdEnergyDepositFinderKG;


#undef _CFMDEBUG_


const std::string gPrintPrefix = " ==> CFMatrixBuilder:";


void
CFMatrixCalculator::Init()
{
  delete fLateralLightCalculator;
  fLateralLightCalculator = new LateralLightCalculator();

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

  lldBranch.GetChild("multipleScattering").GetData(fDoMultipleScattering);
  ostringstream info;
  info << "\n\t multiple scattering is "
       << (fDoMultipleScattering ? "on" : "off");
  INFO(info);
}


CFMatrixCalculator::~CFMatrixCalculator()
{
  delete fLateralLightCalculator;
}


void
CFMatrixCalculator::BuildMatrix(const fevt::Eye& eye, const bool leavingAtmoIsError, const unsigned int addDense)
{
  Clear();

  if (InitCalculation(eye) && CalculateTelescopeData(eye, leavingAtmoIsError, addDense)) {

    if (fVerbosity > -1)
      cout << gPrintPrefix << " building CF-matrix, dim="
           << fNumberOfDepthBins << "..." << endl;

    CalculateDiagonalParameters();
    CalculateFluorescenceMatrix();
    CalculateDirectCherenkovMatrix();

    if (fOnlyDirect) {
      fDirectCFMatrix = fDirectCherenkovMatrix + fDirectFluorescenceMatrix;
      return;
    }

    CalculateMieAndRayScattCherenkovMatrix();

    fCFMatrix = fMieScatteredCherenkovMatrix +
      fRayScatteredCherenkovMatrix +
      fDirectCherenkovMatrix +
      fDirectFluorescenceMatrix;

    if (fVerbosity > -1)
      cout << gPrintPrefix << " done " << endl;

  }
}


bool
CFMatrixCalculator::InitCalculation(const fevt::Eye& eye)
{
  fEyeId = eye.GetId();

  // --- Assert that there's a successful aperture light reconstruction
  if (!eye.HasRecData())
    return false;
  const auto& eyeRecData = eye.GetRecData();
  if (!eyeRecData.HasLightFlux(fevt::FdConstants::eTotal))
    return false;

  // --- Xmax
  const auto& shower = eye.GetRecData().GetFRecShower();
  if (shower.HasGHParameters() && shower.GetGHParameters().GetXMax() > 0)
    fXmax = shower.GetGHParameters().GetXMax();
  else
    fXmax = 750*g/cm2;

  // --- set the Cherenkov model energy threshold for this calculation
  fElectronEnergyThreshold = shower.GetEnergyCutoff();
  const auto& atmo = det::Detector::GetInstance().GetAtmosphere();
  atmo.SetCherenkovEnergyCutoff(fElectronEnergyThreshold);

  // --- wave lengths
  fCherWaveLength = atmo.GetWavelengths(Atmosphere::eCerenkov);
  fFluoWaveLength = atmo.GetWavelengths(Atmosphere::eFluorescence);

  // remove wave lengths with no or zero efficiency
  const auto& fdDetector = Detector::GetInstance().GetFDetector();
  // note: we assume all telescope-effs are equal --> get first one...
  const auto& relEff = fdDetector.GetEye(fEyeId).TelescopesBegin()->GetMeasuredRelativeEfficiency();
  const double effMinLam = relEff.GetX(0);
  const double effMaxLam = relEff.GetX(relEff.GetNPoints() - 1);

  for (auto lIt = fFluoWaveLength.begin(); lIt != fFluoWaveLength.end(); ) {
    if (effMinLam <= *lIt && *lIt <= effMaxLam && relEff.Y(*lIt) > 0)
      ++lIt;
    else
      lIt = fFluoWaveLength.erase(lIt);
  }

  for (auto lIt = fCherWaveLength.begin(); lIt != fCherWaveLength.end(); ) {
    // in contrast to fluorescence, C-bins are integration bounds
    if (effMinLam <= *lIt && *lIt <= effMaxLam)
      ++lIt;
    else
      lIt = fCherWaveLength.erase(lIt);
  }

  return true;
}


bool
CFMatrixCalculator::CalculateTelescopeData(const fevt::Eye& eye, const bool leavingAtmoIsError, const unsigned int addDense)
{
  const auto& shower = eye.GetRecData().GetFRecShower();
  const auto& axis = shower.GetAxis();
  fShowerAxis = -shower.GetAxis(); // downward!
  const auto& core = shower.GetCorePosition();
  const auto& coreCS = fwk::LocalCoordinateSystem::Create(core);
  // should be local (in air) cos theta like in ShowerPhotonGenerator
  // fCosTheta = axis.GetCosTheta(coreCS);

  const auto& coreTime = shower.GetCoreTime();
  const auto& eyeGPSTime = eye.GetHeader().GetTimeStamp();
//#warning hardcoded trace time length!
  const auto eyeTriggerTime = eyeGPSTime - TimeInterval(1000*100*utl::ns);

  const auto& atmo = Detector::GetInstance().GetAtmosphere();
  try {
    atmo.InitSlantProfileModel(core, axis, fMethod == ePrecise ? 10*g/cm2 : 100*g/cm2);
  } catch (InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
    ERROR(e.GetMessage());
    return false;
  }

  const auto& slantDepthProfile = atmo.EvaluateSlantDepthVsDistance();
  const auto& slantDistProfile = atmo.EvaluateDistanceVsSlantDepth();
  const double atmMinDistance = slantDepthProfile.MinX() + 1*mm;
  const double atmMaxDistance = slantDepthProfile.MaxX() - 1*mm;
  const auto& verticalDepthProfile = atmo.EvaluateDepthVsHeight();
  const double atmMinHeight = verticalDepthProfile.MinX();
  const double atmMaxHeight = verticalDepthProfile.MaxX();

  for (auto telIter = eye.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
       telIter != eye.TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++telIter) {

    if (!telIter->HasRecData())
      continue;

    // Fetch light fluxes from the TelescopeRecData
    const auto& telRecData = telIter->GetRecData();
    if (!telRecData.HasLightFlux(fevt::FdConstants::eTotal) ||
        !telRecData.HasLightFlux(fevt::FdConstants::eBackground))
      continue;

    const auto& photonTrace = telRecData.GetLightFlux(fevt::FdConstants::eTotal);
    const auto& bgTrace = telRecData.GetLightFlux(fevt::FdConstants::eBackground);

    if (photonTrace.GetNPoints() != bgTrace.GetNPoints()) {
      ERROR("incompatible light fluxes --> skip tel!!");
      continue;
    }

    // Fetch the light collection efficiencies from the TelescopeRecData
    // or fake them if they're not available
    typedef unique_ptr<TabulatedFunctionErrors> TabFuncErrAutoPtr;
    TabFuncErrAutoPtr fluorLCEff(CalculateLCEff(FdConstants::eFluorDirect, telRecData));
    TabFuncErrAutoPtr directCherLCEff(CalculateLCEff(FdConstants::eCherDirect, telRecData));
    TabFuncErrAutoPtr rayCherLCEff(CalculateLCEff(FdConstants::eCherRayleighScattered, telRecData));
    TabFuncErrAutoPtr mieCherLCEff(CalculateLCEff(FdConstants::eCherMieScattered, telRecData));

    const auto& zetaPixelIds = telRecData.GetPixelsInZetaOverTime();
    if (fMethod == ePrecise && zetaPixelIds.size() != photonTrace.GetNPoints()) {
      ostringstream errMsg;
      errMsg << "inconsistent trace/zetaPixels: N(trace) = "
             << photonTrace.GetNPoints() << " N(zetaPixels)="
             << zetaPixelIds.size();
      ERROR(errMsg);
      throw NonExistentComponentException(errMsg.str());
    }

    // Currently only the largest portion of the lightflux is
    // taken into account.
    // (Which is the whole time range anyway if it was calculated
    //  by FdLightCollectionEfficiencyKG, but not necessarily if it's
    //  from FdApertureLightFinderKG.)
    bool firstTime = true;
    double maxRange = 0;
    double minTime = 0;
    double maxTime = 0;
    double maxSignal = 0;
    const auto& timeRangeList = telRecData.GetSpotFarFromBorderTimeRanges();
    for (auto timeRangeIter = timeRangeList.begin(), end = timeRangeList.end(); timeRangeIter != end; ++timeRangeIter) {

      const double thisRange = fabs(timeRangeIter->first - timeRangeIter->second);
      const double thisRangeMin = fmin(timeRangeIter->first, timeRangeIter->second);
      const double thisRangeMax = fmax(timeRangeIter->first, timeRangeIter->second);

      const auto& flux = telRecData.GetLightFlux();
      double thisSignal = 0;
      for (auto it = flux.Begin(), end = flux.End(); it != end; ++it) {
        if (thisRangeMin <= it->X() && it->X() <= thisRangeMax)
          thisSignal += it->Y() / it->YErr();
      }

      // Use the integrated signal of the time range to decide
      if (firstTime || thisSignal > maxSignal) {
        firstTime = false;
        maxSignal = thisSignal;
        maxRange = thisRange;
        minTime = thisRangeMin;
        maxTime = thisRangeMax;
      }
    } // end foreach time ranges

    if (maxRange <= 0) {
      ostringstream errMsg;
      errMsg << " zero time range for telId=" << telIter->GetId() << "??";
      ERROR(errMsg);
      continue;
    }

    const auto& fdDetector = Detector::GetInstance().GetFDetector();
    const auto& detTel = fdDetector.GetEye(eye.GetId()).GetTelescope(telIter->GetId());
    const unsigned int physicalEyeId = detTel.GetParentPhysicalEyeId();
    const auto& telEyeCS = fdDetector.GetEye(physicalEyeId).GetEyeCoordinateSystem();

    // for distance from shower core to emission point (aDist)
    // law of cosines: bDist^2 = aDist^2 + cDist^2 - 2*aDist*cDist*cos(beta)
    // cDist = distance(telescope, core)
    // bDist-aDist = (telTime-coreTime)/c = delta (see below)
    const auto coreTelescopeVec = detTel.GetPosition() - core;
    const double cDist = coreTelescopeVec.GetMag();
    const double cosBeta = CosAngle(axis, coreTelescopeVec);

    TelescopeData telData(telIter->GetId());
    //noiseTelData contains everything as telData plus bins coming from
    //light corresponding to the part of shower outside of the atmosphere - NOT FOR CFM calculation!
    TelescopeData noiseTelData(telIter->GetId());
    for (unsigned int i = 0; i < photonTrace.GetNPoints(); ++i) {

      // divide real photonTrace bins into several dense (virtual) bins for finer CFM in ProfileFitter
      for (unsigned int j = 0; j < addDense; ++j) {

        /*
        const double tMean = photonTrace.GetX(i);
        const double halfBinWidth = photonTrace.GetXErr(i);
        const double tFirst = tMean-halfBinWidth;
        const double tLast = tMean+halfBinWidth;
        */
        // definition according to dense binning
        const double halfBinWidth = photonTrace.GetXErr(i);
        const double tMean = photonTrace.GetX(i) - halfBinWidth + (2*j + 1) * halfBinWidth / addDense;
        const double tFirst = tMean-halfBinWidth / addDense;
        const double tLast = tMean+halfBinWidth / addDense;

        /*
        if (tFirst < minTime)
          continue;
        if (tLast > maxTime)
          break;
        */
        if (photonTrace.GetX(i) - halfBinWidth < minTime)
          break;
        if (photonTrace.GetX(i) + halfBinWidth > maxTime)
          break;

        try {

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

          const auto firstTimeAtTelescope = eyeTriggerTime + tFirst;
          const double deltaFirst = (firstTimeAtTelescope - coreTime).GetInterval()*kSpeedOfLight;
          const double firstDistance =
            -(cDist*cDist - deltaFirst*deltaFirst) / (2*(deltaFirst + cDist*cosBeta));
          const auto firstPointOnShower = core - firstDistance*axis;
          const double firstHeight = wgs84.PointToLatitudeLongitudeHeight(firstPointOnShower).get<2>();

          if (firstDistance < atmMinDistance || firstDistance > atmMaxDistance ||
              firstHeight < atmMinHeight || firstHeight > atmMaxHeight) {
            if (leavingAtmoIsError)
              return false;
            else {
              map<FdConstants::LightSource, double> lightEfficiencies;
              TelescopeDataBin telDataBin(0.,0.,0., Point(), Point(), Point(),
                                      tMean,
                                      halfBinWidth / addDense,
                                      photonTrace.GetY(i),
                                      photonTrace.GetYErr(i),
                                      pow(bgTrace.GetYErr(i), 2),
                                      lightEfficiencies, i);
              noiseTelData.AddTelescopeDataBin(telDataBin);
              continue;
            }
          }

          const double firstDepth = slantDepthProfile.Y(firstDistance);

          TimeStamp lastTimeAtTelescope = eyeTriggerTime + tLast;
          const double deltaLast =
            (lastTimeAtTelescope-coreTime).GetInterval()*kSpeedOfLight;
          const double lastDistance =
            -(cDist*cDist-deltaLast*deltaLast)/(2*(deltaLast+cDist*cosBeta));
          const Point lastPointOnShower = core-lastDistance*axis;
          const double lastHeight =
            wgs84.PointToLatitudeLongitudeHeight(lastPointOnShower).get<2>();

          if (lastDistance < atmMinDistance || lastDistance > atmMaxDistance ||
              lastHeight < atmMinHeight   || lastHeight > atmMaxHeight) {
            if (leavingAtmoIsError)
              return false;
            else {
              map<FdConstants::LightSource, double> lightEfficiencies;
              TelescopeDataBin telDataBin(0.,0.,0.,Point(),Point(),Point(),
                                      tMean,
                                      halfBinWidth/addDense,
                                      photonTrace.GetY(i),
                                      photonTrace.GetYErr(i),
                                      pow(bgTrace.GetYErr(i),2),
                                      lightEfficiencies, i);
              noiseTelData.AddTelescopeDataBin(telDataBin);
              continue;
            }
          }

          const double lastDepth = slantDepthProfile.Y(lastDistance);
          const double meanDepth = (firstDepth + lastDepth) / 2.;
          const double meanDistance = slantDistProfile.Y(meanDepth);
          const Point meanPointOnShower = core - meanDistance * axis;
          const double meanHeight =
            wgs84.PointToLatitudeLongitudeHeight(meanPointOnShower).get<2>();
          if (meanHeight < atmMinHeight   || meanHeight > atmMaxHeight) {
            if (leavingAtmoIsError)
              return false;
            else {
              map<FdConstants::LightSource, double> lightEfficiencies;
              TelescopeDataBin telDataBin(0.,0.,0.,Point(),Point(),Point(),
                                      tMean,
                                      halfBinWidth/addDense,
                                      photonTrace.GetY(i),
                                      photonTrace.GetYErr(i),
                                      pow(bgTrace.GetYErr(i),2),
                                      lightEfficiencies, i);
              noiseTelData.AddTelescopeDataBin(telDataBin);
              continue;
            }
          }

          if (!isfinite(firstDepth) || !isfinite(lastDepth))
            continue;

          map<FdConstants::LightSource, double> lightEfficiencies;
          if (fluorLCEff.get())
            lightEfficiencies[FdConstants::eFluorDirect] =
              fluorLCEff->GetY(i);
          if (directCherLCEff.get())
            lightEfficiencies[FdConstants::eCherDirect] =
              directCherLCEff->GetY(i);
          if (rayCherLCEff.get())
            lightEfficiencies[FdConstants::eCherRayleighScattered] =
              rayCherLCEff->GetY(i);
          if (mieCherLCEff.get())
            lightEfficiencies[FdConstants::eCherMieScattered] =
              mieCherLCEff->GetY(i);

          TelescopeDataBin telDataBin(fmin(firstDepth,lastDepth),
                                      fmax(firstDepth,lastDepth),
                                      meanHeight,
                                      firstPointOnShower,
                                      lastPointOnShower,
                                      meanPointOnShower,
                                      tMean,
                                      halfBinWidth/addDense,
                                      photonTrace.GetY(i),
                                      photonTrace.GetYErr(i),
                                      pow(bgTrace.GetYErr(i),2),
                                      lightEfficiencies, i);

          // ----  add zeta pixels
          //  pixel coordinates in CS with x-axis pointing to shower

          if (fMethod == ePrecise) {
            Vector vecToShower = meanPointOnShower-detTel.GetPosition();
            const double vecToShowerPhi = vecToShower.GetPhi(telEyeCS);
            const double vecToShowerElevation =
              kPiOnTwo-vecToShower.GetTheta(telEyeCS);
            const CoordinateSystemPtr auxCS =
              CoordinateSystem::RotationZ(vecToShowerPhi, telEyeCS);
            const CoordinateSystemPtr toShowerCS =
              CoordinateSystem::RotationY(-vecToShowerElevation, auxCS);

            const vector<unsigned int> pixelIdVec = zetaPixelIds[i];
            for (unsigned int iPix=0; iPix < pixelIdVec.size(); ++iPix) {

              const unsigned int pixId = pixelIdVec[iPix];
              const fdet::Pixel& thisPixel = detTel.GetPixel(pixId);
              const Vector& pixelDir = thisPixel.GetDirection();
              const double x = pixelDir.GetPhi(toShowerCS);
              const double y = kPiOnTwo-pixelDir.GetTheta(toShowerCS);
              // side length from solid angle
              const double omega = thisPixel.GetSolidAngle();
              const double sideLength = sqrt(2./3.*omega/kSqrt3);

              double tiltAngle = 0;

              // local tilt from slope to pixel from next column
              // (if not a single pixel telescope)
              if (detTel.GetFirstColumn() != detTel.GetLastColumn()) {
                const unsigned int thisRow = thisPixel.GetRow();
                const unsigned int thisColumn = thisPixel.GetColumn();
                const unsigned int nextColumn =
                  (thisColumn == 1 ? thisColumn + 1: thisColumn - 1);
                const fdet::Pixel& nextPixel = detTel.GetPixel(thisRow, nextColumn);
                const double nextX = nextPixel.GetDirection().GetPhi(toShowerCS);
                const double nextY =
                  kPiOnTwo-nextPixel.GetDirection().GetTheta(toShowerCS);
                tiltAngle = atan2(nextY - y, nextX -x );
                if (thisColumn == 1)
                  tiltAngle += kPi;
              }

              telDataBin.AddZetaPixel( ZetaPixel(pixId, x, y, sideLength, tiltAngle));

            }
          }

          telData.AddTelescopeDataBin(telDataBin);
          noiseTelData.AddTelescopeDataBin(telDataBin);

        } catch (UTMPoint::UTMZoneException& d) {
          ERROR(d.GetMessage());
          continue;
        }

      }

    } // end for each photon trace bin

    if (!telData.GetTelescopeDataBins().empty()) {
      telData.SortBins();
      noiseTelData.SortBins();
      SetTelescopeParameters(eye, telData);
      SetTelescopeParameters(eye, noiseTelData);
      telData.SetZeta(telRecData.GetZeta());
      noiseTelData.SetZeta(telRecData.GetZeta());
      fTelescopeData.push_back(telData);
      fNoiseTelescopeData.push_back(noiseTelData);
    }

  } // end for each telescope

  if (!fTelescopeData.empty()) {

    fFOVsize = 0;
    for (std::vector<TelescopeData>::const_iterator
           telIter = fTelescopeData.begin();
         telIter != fTelescopeData.end(); ++telIter)
      fFOVsize += telIter->GetTelescopeDataBins().size();

    AddBinsOutsideFOV(eye);

    fNumberOfDepthBins = 0;
    sort(fTelescopeData.begin(), fTelescopeData.end());
    for (std::vector<TelescopeData>::const_iterator
           telIter = fTelescopeData.begin();
         telIter != fTelescopeData.end(); ++telIter)
      fNumberOfDepthBins += telIter->GetTelescopeDataBins().size();

    sort(fNoiseTelescopeData.begin(), fNoiseTelescopeData.end());

    return true;
  }

  return false;
}


void
CFMatrixCalculator::AddBinsOutsideFOV(const fevt::Eye& eye)
{
  double minDepth = 0;
  bool first = true;

  for (std::vector<TelescopeData>::const_iterator telIter = fTelescopeData.begin();
       telIter != fTelescopeData.end(); ++telIter) {

    const double minTelDepth = telIter->GetMinDepth();
    if (first || minTelDepth < minDepth) {
      minDepth = minTelDepth;
      first = false;
    }
  }

  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  const ProfileResult& distProfile  = atmo.EvaluateDistanceVsSlantDepth();
  const double atmoMinDepth = distProfile.MinX();
  const ProfileResult& verticalDepthProfile = atmo.EvaluateDepthVsHeight();
  const double atmMinHeight = verticalDepthProfile.MinX();
  const double atmMaxHeight = verticalDepthProfile.MaxX();
  const double dX = 25.*(g/cm2); // depth step size

  const ShowerFRecData& shower = eye.GetRecData().GetFRecShower();
  const Vector& axis = shower.GetAxis();
  const Point& core = shower.GetCorePosition();

  // no C-light from below detector
  const Point& detEyePos =
    Detector::GetInstance().GetFDetector().GetEye(fEyeId).GetPosition();
  const double minHeight =
    wgs84.PointToLatitudeLongitudeHeight(detEyePos).get<2>();

  // push bins into fake telescope

  TelescopeData telData(0, TelescopeData::eOutsideFOV);
  telData.SetZeta(1.3*degree);

  double thisDepth = minDepth;
  while (thisDepth > atmoMinDepth+dX) {
    const double nextDepth = thisDepth - dX;
    const Point firstPointOnShower = core-distProfile.Y(nextDepth)*axis;
    const Point lastPointOnShower = core-distProfile.Y(thisDepth)*axis;
    const double meanDepth = (thisDepth+nextDepth)/2.;
    const double meanDistance = distProfile.Y(meanDepth);
    const Point meanPointOnShower = core-meanDistance*axis;
    const map<FdConstants::LightSource, double> dummyEff;
    const double meanHeight =
      wgs84.PointToLatitudeLongitudeHeight(meanPointOnShower).get<2>();
    if (meanHeight > minHeight && meanHeight > atmMinHeight &&
        meanHeight < atmMaxHeight && fmin(thisDepth, nextDepth) > atmoMinDepth) {
      telData.AddTelescopeDataBin(TelescopeDataBin(fmin(thisDepth,nextDepth),
                                                   fmax(thisDepth,nextDepth),
                                                   meanHeight,
                                                   firstPointOnShower,
                                                   lastPointOnShower,
                                                   meanPointOnShower,
                                                   0, 0, 0, 0, 0,
                                                   dummyEff, -1));
    }

    thisDepth = nextDepth;
  }

  if (fVerbosity > -1)
    cout << gPrintPrefix
         << " added " << telData.GetTelescopeDataBins().size()
         << " bins outside FOV (X<" << minDepth/g*cm2 << ")"
         << endl;

  if (!telData.GetTelescopeDataBins().empty()) {
    telData.SortBins();
    fTelescopeData.push_back(telData);
    fNoiseTelescopeData.push_back(telData);
  }
}


void
CFMatrixCalculator::CalculateDiagonalParameters()
{
  fCherTransmissionToTel.clear();
  fFluoTransmissionToTel.clear();

  const Atmosphere& atmo = Detector::GetInstance().GetAtmosphere();
  const fdet::Eye& detEye =
    Detector::GetInstance().GetFDetector().GetEye(fEyeId);

  // vector of C/F wave lengths
  const std::vector<double>& lambdaCVec = fCherWaveLength;
  const unsigned int nCLambda = lambdaCVec.size();
  vector<double> cherenkovTelTransmission(nCLambda);
  vector<double> cherenkovTrackTransmission(nCLambda);
  vector<double> cherenkovRayleighScattering(nCLambda);
  vector<double> cherenkovMieScattering(nCLambda);
  vector<double> cherenkovProduction(nCLambda);

  const std::vector<double>& lambdaFVec = fFluoWaveLength;
  const unsigned int nFLambda = lambdaFVec.size();
  std::vector<double> fluorescenceTransmission(nFLambda);

  const double fourPi = 4.0*kPi;
  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();

  for (vector<TelescopeData>::const_iterator telIter = fTelescopeData.begin();
       telIter != fTelescopeData.end(); ++telIter) {

    const Point& telPos = telIter->GetType() == TelescopeData::eInsideFOV?
      detEye.GetTelescope(telIter->GetId()).GetPosition():
      detEye.GetPosition();

    for (vector<TelescopeDataBin>::const_iterator binIter =
           telIter->TelDataBinsBegin();
         binIter !=  telIter->TelDataBinsEnd(); ++binIter) {

      const Point& meanPos = binIter->fMeanDepthPoint;

      // ---- geometrical factor

      const Vector showerTelVec = telPos-meanPos;
      const double distToTelSquared = showerTelVec.GetMag2();
      const double distToTel = sqrt(distToTelSquared);
      fGeometricalFactor.push_back(1./fourPi/distToTelSquared);
      fZetaDistance.push_back(tan(telIter->GetZeta())*distToTel);

      // ---- average energy deposit per electron
      const double alpha = utl::EnergyDeposit(binIter->GetMeanDepth(),
                                              fXmax,
                                              fElectronEnergyThreshold);
      fMeandEdXPerElectron.push_back(alpha);

      // ---- local cosTheta
      const double cosLocalZenith =
          (-fShowerAxis).GetCosTheta(fwk::LocalCoordinateSystem::Create(UTMPoint(meanPos, wgs84)));
      fCosTheta.push_back(cosLocalZenith);

      // ---- attenuation to telescope for cherenkov

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

      // ---- attenuation to telescope for fluorescence

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

      // --- scattering to telescope

      const double emissionAngle = Angle(fShowerAxis, showerTelVec);
      ScatteringResult rScattering =
        atmo.EvaluateRayleighScattering(binIter->fFirstPoint,
                                        binIter->fLastPoint,
                                        emissionAngle, distToTel, lambdaCVec);
      ScatteringResult mScattering =
        atmo.EvaluateMieScattering(binIter->fFirstPoint,
                                   binIter->fLastPoint,
                                   emissionAngle, distToTel, lambdaCVec);
      const utl::TabulatedFunctionErrors& rayleighScattTrackTel =
        rScattering.GetScatteringFactor();
      const utl::TabulatedFunctionErrors& mieScattTrackTel =
        mScattering.GetScatteringFactor();

      // --- transmission along track
      const AttenuationResult mAtt =
        atmo.EvaluateMieAttenuation(binIter->fFirstPoint,
                                    binIter->fLastPoint, lambdaCVec);
      const utl::TabulatedFunctionErrors& mieTrackAtt =
        mAtt.GetTransmissionFactor();
      const AttenuationResult rAtt =
        atmo.EvaluateRayleighAttenuation(binIter->fFirstPoint,
                                         binIter->fLastPoint, lambdaCVec);
      const utl::TabulatedFunctionErrors& rayTrackAtt =
        rAtt.GetTransmissionFactor();

      // --- cherenkov at track
      const double showerAge = utl::ShowerAge(binIter->GetMeanDepth(), fXmax);
      const utl::TabulatedFunction& cherenkovProd =
        atmo.EvaluateCherenkovPhotons(binIter->fFirstPoint,
                                      binIter->fLastPoint, showerAge);

      // --- finally copy to data members
      for (unsigned int j=0; j<nCLambda; ++j)  {
        cherenkovRayleighScattering[j] = rayleighScattTrackTel.GetY(j);
        cherenkovMieScattering[j] = mieScattTrackTel.GetY(j);
        cherenkovTelTransmission[j] = mieCAtt.GetY(j)*rayCAtt.GetY(j);
        cherenkovTrackTransmission[j] = rayTrackAtt.GetY(j)*mieTrackAtt.GetY(j);
        cherenkovProduction[j] = cherenkovProd.GetY(j);
      }

      for (unsigned int j=0; j<nFLambda; ++j)
        fluorescenceTransmission[j] = mieFAtt.GetY(j) * rayFAtt.GetY(j);

      fFluoTransmissionToTel.push_back(fluorescenceTransmission);
      fRayScatToTel.push_back(cherenkovRayleighScattering);
      fMieScatToTel.push_back(cherenkovMieScattering);
      fCherTransmissionToTel.push_back(cherenkovTelTransmission);
      fCherenkovAtTrack.push_back(cherenkovProduction);
      fCherTransmissionShower.push_back(cherenkovTrackTransmission);

    }
  }
}


void
CFMatrixCalculator::CalculateFluorescenceMatrix()
{
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  const double dEdX0 = atmo.GetdEdX0();

  const std::vector<double>& lambdaVec = fFluoWaveLength;
  const fdet::FDetector& fdDetector = Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = fdDetector.GetEye(fEyeId);
  fDirectFluorescenceMatrix.SetSize(fNumberOfDepthBins);

  int i = 0;
  for (vector<TelescopeData>::const_iterator telIter = fTelescopeData.begin();
       telIter != fTelescopeData.end(); ++telIter) {

    const bool realTelescope =
      (telIter->GetType() == TelescopeData::eInsideFOV);
    const unsigned int telId = realTelescope? telIter->GetId() :
      detEye.GetFirstTelescopeId();

    const TabulatedFunction& relEff =
      fdDetector.GetEye(fEyeId).GetTelescope(telId).GetMeasuredRelativeEfficiency();

    const Point& telPos = detEye.GetTelescope(telId).GetPosition();

    for (vector<TelescopeDataBin>::const_iterator binIter =
           telIter->TelDataBinsBegin();
         binIter !=  telIter->TelDataBinsEnd(); ++binIter) {

      const double height = binIter->fHeight;

      const utl::TabulatedFunction& fluoyield_tab =
        atmo.EvaluateFluorescenceYield(height);

      double epsYTsum = 0.;
      for (unsigned int j = 0; j < lambdaVec.size(); j++) {
        const double Y_j = fluoyield_tab.Y(lambdaVec[j]);
        epsYTsum += Y_j * relEff.Y(lambdaVec[j]) * fFluoTransmissionToTel[i][j];
      }

      const double deltaL =
        (binIter->fFirstPoint - binIter->fLastPoint).GetMag();
      const double lldFraction =
        fLateralLightCalculator->GetLightFraction(FdConstants::eFluorDirect,
                                                  *binIter, *binIter,
                                                  detEye.GetTelescope(telId),
                                                  fXmax,
                                                  fCosTheta[i],
                                                  telIter->GetZeta());

      const double msFactor =
        1. / (1 - MultipleScatteringFraction(*binIter,
                                             lambdaVec,
                                             fluoyield_tab,
                                             telIter->GetZeta(),
                                             telPos,
                                             relEff));
      fFluorescenceMultipleScattering.push_back(msFactor);

      fDirectFluorescenceMatrix(i,i) = lldFraction * msFactor *
        epsYTsum * fGeometricalFactor[i] / dEdX0 * deltaL;

      ++i;
    }
  }
}


void
CFMatrixCalculator::CalculateDirectCherenkovMatrix()
{
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  const std::vector<double>& lambdaVec = fCherWaveLength;
  const fdet::FDetector& fdDetector = Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = fdDetector.GetEye(fEyeId);

  fDirectCherenkovMatrix.SetSize(fNumberOfDepthBins);

  int i = 0;
  for (vector<TelescopeData>::const_iterator telIter = fTelescopeData.begin();
       telIter != fTelescopeData.end(); ++telIter) {

    const bool realTelescope =
      (telIter->GetType() == TelescopeData::eInsideFOV);
    const unsigned int telId = realTelescope? telIter->GetId() :
      detEye.GetFirstTelescopeId();

    const TabulatedFunction& relEff =
      detEye.GetTelescope(telId).GetMeasuredRelativeEfficiency();

    const Point& telPos = detEye.GetTelescope(telId).GetPosition();

    for (vector<TelescopeDataBin>::const_iterator binIter =
           telIter->TelDataBinsBegin();
         binIter !=  telIter->TelDataBinsEnd(); ++binIter) {

      const double showerAge = utl::ShowerAge(binIter->GetMeanDepth(),fXmax);
      const utl::TabulatedFunction& directCherenkov =
        atmo.EvaluateCherenkovDirect(binIter->fFirstPoint,
                                     binIter->fLastPoint, telPos, showerAge);
      double epsYTsum=0.;
      for (unsigned int j = 0; j < lambdaVec.size(); j++) {
        const double Y_j = directCherenkov.GetY(j);
        epsYTsum += Y_j * relEff.Y(lambdaVec[j]) * fCherTransmissionToTel[i][j];
      }

      const double lldFraction =
        fLateralLightCalculator->GetLightFraction(FdConstants::eCherDirect,
                                                  *binIter, *binIter,
                                                  detEye.GetTelescope(telId),
                                                  fXmax,
                                                  fCosTheta[i],
                                                  telIter->GetZeta());

      const double msFactor =
        1. / (1 - MultipleScatteringFraction(*binIter,
                                             lambdaVec,
                                             directCherenkov,
                                             telIter->GetZeta(),
                                             telPos,
                                             relEff));
      fCherenkovMultipleScattering.push_back(msFactor);

      fDirectCherenkovMatrix(i) =
        epsYTsum / fMeandEdXPerElectron[i] * msFactor * lldFraction;

      ++i;
    }
  }
}


void
CFMatrixCalculator::CalculateMieAndRayScattCherenkovMatrix()
{
  fMieScatteredCherenkovMatrix.SetSize(fNumberOfDepthBins);
  fRayScatteredCherenkovMatrix.SetSize(fNumberOfDepthBins);

  const fdet::FDetector& fdDetector = Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = fdDetector.GetEye(fEyeId);

  const unsigned int nLambda = fCherWaveLength.size();
  vector<double> transmission(nLambda);

  int i = 0;
  for (vector<TelescopeData>::const_iterator telIter_i = fTelescopeData.begin();
       telIter_i != fTelescopeData.end(); ++telIter_i) {
    const bool realTelescope = telIter_i->GetType() == TelescopeData::eInsideFOV;
    const unsigned int telId =
      realTelescope? telIter_i->GetId() : detEye.GetFirstTelescopeId();

    const TabulatedFunction& relEff =
      fdDetector.GetEye(fEyeId).GetTelescope(telId).GetMeasuredRelativeEfficiency();

    for (vector<TelescopeDataBin>::const_iterator binIter_i = telIter_i->TelDataBinsBegin();
         binIter_i !=  telIter_i->TelDataBinsEnd(); ++binIter_i) {

      const double msFactor = fCherenkovMultipleScattering[i];

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

      const bool realTelescope =
        (telIter_i->GetType() == TelescopeData::eInsideFOV);
      const unsigned int telId = realTelescope? telIter_i->GetId():
        detEye.GetFirstTelescopeId();

      for (int j = i; j >= 0; --j) {

        // check for telescope overlaps that might (should!) occur
        // in lightCollectionEfficiency-mode
        if (IsOverlapBin(i,j)) {
          fMieScatteredCherenkovMatrix(i, j) = 0;
          fRayScatteredCherenkovMatrix(i, j) = 0;
          continue;
        }

        double epsYTsumMie = 0;
        double epsYTsumRay = 0;
        for (unsigned int k = 0; k < nLambda; ++k) {

          if (j != i)
            transmission[k]*=fCherTransmissionShower[j][k];

          const double detEff = relEff.Y(fCherWaveLength[k]);
          const double factor =
            transmission[k]*detEff*fCherenkovAtTrack[j][k]*
            fCherTransmissionToTel[i][k];

          epsYTsumMie += factor*fMieScatToTel[i][k];
          epsYTsumRay += factor*fRayScatToTel[i][k];

        }

        const double oneOverAlpha = 1 / fMeandEdXPerElectron[j];

        fMieScatteredCherenkovMatrix(i, j) = epsYTsumMie * oneOverAlpha;
        fRayScatteredCherenkovMatrix(i, j) = epsYTsumRay * oneOverAlpha;

        const TelescopeDataBin& emissionBin = *(GetTelescopeDataBin(j).second);
        const TelescopeDataBin& detectionBin = *binIter_i;

        // no difference between eCherMieScattered and eCherRayScattered
        FdConstants::LightSource lightType = FdConstants::eCherMieScattered;

        const double lldFraction =
          fLateralLightCalculator->GetLightFraction(lightType,
                                                    emissionBin, detectionBin,
                                                    detEye.GetTelescope(telId),
                                                    fXmax,
                                                    fCosTheta[i],
                                                    telIter_i->GetZeta());
        fRayScatteredCherenkovMatrix(i,j) *= (lldFraction*msFactor);
        fMieScatteredCherenkovMatrix(i,j) *= (lldFraction*msFactor);

      }

      ++i;
    }
  }
}


utl::TabulatedFunctionErrors*
CFMatrixCalculator::CalculateLCEff(const FdConstants::LightSource lightSource,
                                   const TelescopeRecData& telRecData)
  const
{
  utl::TabulatedFunctionErrors* lcEff = nullptr;

  // fetch lceff from TelescopeRecData if available...
  if (telRecData.HasLightCollectionEfficiency()) {
    const MultiTabulatedFunctionErrors& lcEffs =
      telRecData.GetLightCollectionEfficiency();
    if (lcEffs.HasLabel(lightSource)) {
      const utl::TabulatedFunctionErrors& tmp =
        lcEffs.GetTabulatedFunctionErrors(lightSource);
      lcEff = new utl::TabulatedFunctionErrors(tmp);
    }
  }

  return lcEff;
}


void
CFMatrixCalculator::Clear()
{
  fCFMatrix.Clear();
  fMieScatteredCherenkovMatrix.Clear();
  fRayScatteredCherenkovMatrix.Clear();
  fDirectCherenkovMatrix.Clear();
  fDirectFluorescenceMatrix.Clear();
  fTelescopeData.clear();
  fNoiseTelescopeData.clear();
  fCherWaveLength.clear();
  fFluoWaveLength.clear();
  fCherTransmissionToTel.clear();
  fFluoTransmissionToTel.clear();
  fCherTransmissionShower.clear();
  fRayScatToTel.clear();
  fMieScatToTel.clear();
  fFluorescenceMultipleScattering.clear();
  fCherenkovMultipleScattering.clear();
  fCherenkovAtTrack.clear();
  fGeometricalFactor.clear();
  fZetaDistance.clear();
  fMeandEdXPerElectron.clear();
}


pair<const TelescopeData*, const TelescopeDataBin*>
CFMatrixCalculator::GetTelescopeDataBin(const unsigned int i)
  const
{
  int j = 0;
  for (std::vector<TelescopeData>::const_iterator
         telIter = fTelescopeData.begin();
       telIter != fTelescopeData.end(); ++telIter) {
    const unsigned int thisSize =
      telIter->TelDataBinsEnd()-telIter->TelDataBinsBegin();
    if (j + thisSize > i) {
      const TelescopeData* tel = &(*telIter);
      const TelescopeDataBin* data = &(*(telIter->TelDataBinsBegin()+(i-j)));
      return pair<const TelescopeData*, const TelescopeDataBin*>(tel, data);
    } else
      j += thisSize;
  }

  // if reach this point --> error
  ostringstream errMsg;
  errMsg << " request bin " << i << " in ";
  for (std::vector<TelescopeData>::const_iterator telIter = fTelescopeData.begin();
       telIter != fTelescopeData.end(); ++telIter)
    errMsg << "(tel " << telIter->GetId() << " n="
           <<  telIter->TelDataBinsEnd()-telIter->TelDataBinsBegin()
           << ") ";

  ERROR(errMsg.str());
  throw OutOfBoundException(errMsg.str());

  // to make compiler happy
  return pair<const TelescopeData*, const TelescopeDataBin*>
    (&fTelescopeData[0], &(*fTelescopeData[0].TelDataBinsBegin()));
}


bool
CFMatrixCalculator::IsOverlapBin(const int i, const int j)
  const
{
  // check causality
  if (j > i)
    return true;

  typedef std::pair<const TelescopeData*, const TelescopeDataBin*> TelDataPair;
  TelDataPair data_j = GetTelescopeDataBin(j);
  TelDataPair data_i = GetTelescopeDataBin(i);

  // same telescope is always OK
  if (data_i.first->GetId() == data_j.first->GetId())
    return false;

  // different telescope j, but depth_j in telescope i is not OK
  if (data_i.first->DepthInRange(data_j.second->fMinDepth) ||
      data_i.first->DepthInRange(data_j.second->fMaxDepth)) {
#ifdef _CFMDEBUG_
    cout << " overlap type1, i=" << i << " j=" << j << " depth(j)="
         << data_j.second->GetMeanDepth()/g*cm2 << " tel_j="
         << data_j.first->GetId()
         << " tel_i " << data_i.first->GetId() << " X=["
         << data_i.first->GetMinDepth()/g*cm2 << ","
         << data_i.first->GetMaxDepth()/g*cm2 << "]" << endl;
#endif
    return true;
  }

  // depth_j occurs the first time in tel_j is OK
  for (vector<TelescopeData>::const_iterator telIter = fTelescopeData.begin();
       telIter != fTelescopeData.end(); ++telIter) {
    if (telIter->DepthInRange(data_j.second->fMinDepth) ||
        telIter->DepthInRange(data_j.second->fMaxDepth)) {
      if (telIter->GetId() != data_j.first->GetId()) {
#ifdef _CFMDEBUG_
        cout << " overlap type2, i=" << i << " j=" << j << " depth(j)="
             << data_j.second->GetMeanDepth()/g*cm2 << " tel_j="
             << data_j.first->GetId()
             << " tel_k " << telIter->GetId() << " X=["
             << telIter->GetMinDepth()/g*cm2 << ","
             << telIter->GetMaxDepth()/g*cm2 << "]" << endl;
#endif
        return true;
      }
      else
        return false;
    }
  }

  // this part should never be reached
  ostringstream errMsg;
  errMsg << " logical error for i=" << i << " j=" << j;
  ERROR(errMsg.str());
  throw OutOfBoundException(errMsg.str());

  return true;
}


void
CFMatrixCalculator::SetTelescopeParameters(const fevt::Eye& eye,
                                           TelescopeData& telData)
  const
{
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const double referenceLambda = detFD.GetReferenceLambda();
  const fevt::Telescope& evtTel = eye.GetTelescope(telData.GetId());

  double avgPEfactor = 0;
  int countPix = 0;
  for (fevt::Telescope::ConstPixelIterator pixiter =
         evtTel.PixelsBegin(ComponentSelector::eHasData);
       pixiter != evtTel.PixelsEnd(ComponentSelector::eHasData);
       ++pixiter) {
    const fdet::Pixel& detPixel = detFD.GetPixel(*pixiter);
    avgPEfactor += detPixel.GetDiaPhoton2PEFactor(referenceLambda);
    ++countPix;
  }

  if (countPix)
    avgPEfactor /= countPix;

  const fdet::Telescope& detTel = detFD.GetTelescope(evtTel);

  telData.SetTelescopeParameters(avgPEfactor,
                                 detTel.GetCamera().GetGainVariance(),
                                 detTel.GetDiaphragmArea());
}


/******************************************************************\
 *
 *  MultipleScatteringFraction()
 *
 *  fraction of multiple scattered light wrt. total signal from
 *
 *   M. Roberts, J.Phys.G31:1291,2005
 *
 *
\******************************************************************/
double
CFMatrixCalculator::MultipleScatteringFraction(const TelescopeDataBin& telDataBin,
                                               const std::vector<double>& waveLengths,
                                               const utl::TabulatedFunction& yield,
                                               const double zeta,
                                               const Point& telescopePosition,
                                               const TabulatedFunction& efficiency)
  const
{
  if (!fDoMultipleScattering)
    return 0;

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

  // optical depth and scattering coefficient

  //         attenuation shower-->eye
  const utl::Point& meanPos = telDataBin.fMeanDepthPoint;
  const AttenuationResult mieResultAttShowerToEye =
    atmo.EvaluateMieAttenuation(telescopePosition,
                                meanPos,
                                waveLengths);
  const utl::TabulatedFunctionErrors& mieAttShowerToEye =
    mieResultAttShowerToEye.GetTransmissionFactor();

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

  //         attenuation within bin
  const AttenuationResult mieAttResultInBin =
    atmo.EvaluateMieAttenuation(telDataBin.fFirstPoint,
                                telDataBin.fLastPoint,
                                waveLengths);
  const utl::TabulatedFunctionErrors& mieAttInBin =
    mieAttResultInBin.GetTransmissionFactor();

  const AttenuationResult rayAttResultInBin =
    atmo.EvaluateRayleighAttenuation(telDataBin.fFirstPoint,
                                     telDataBin.fLastPoint,
                                     waveLengths);
  const utl::TabulatedFunctionErrors& rayAttInBin =
    rayAttResultInBin.GetTransmissionFactor();

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

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

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

  const double distanceToEye = (meanPos-telescopePosition).GetMag();
  const double binLength = (telDataBin.fFirstPoint-
                            telDataBin.fLastPoint).GetMag();

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

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

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

  return fraction;
}