LongitudinalXmaxScanner.cc

Go to the documentation of this file.

#include <utl/config.h>
#include "LongitudinalXmaxScanner.h"

#include <adst/RecEvent.h>

#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>
#include <utl/Transformation.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/PhysicalFunctions.h>

#include <fwk/CoordinateSystemRegistry.h>

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

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

#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/ReferenceEllipsoid.h>

#include <evt/Event.h>
#include <evt/GaisserHillas4Parameter.h>
#include <evt/ShowerFRecData.h>
#include "evt/ShowerRecData.h"

#include <TMatrixD.h>

#include <iostream>

using namespace std;
using namespace fwk;
using namespace atm;
using namespace det;
using namespace fdet;
using namespace utl;
using namespace otoa;


//--> magic constants ==> needs fix!

const double gcm2 = utl::g/utl::cm2;

// GH constraints
const double kAverageLambda = 61*gcm2;
const double kAverageX0 = -121*gcm2;
// variance of GH constraints
const double kLambdaVar = 13*13*gcm2*gcm2;
const double kX0Var = 172*172*gcm2*gcm2;
// value returned for singular matrices etc
const double kMaxVariance = 999*999*gcm2*gcm2;

//--> end magic constants


inline
double
AverageErr(const double x0, const double x1, const double x2)
{
  const double e1 = x1 - x0;
  const double e2 = x0 - x2;
  return 0.5*(e1 + e2);
}


// 68% space angle given angular track length (MC study from Jose, 09/05/28) updated (2017) taking into account the PCGf reco
inline
double
SpaceAngle(const double angularLength, const FDEvent& fdEvent)
{
  const double refAngle = 30*degree;

  const FdRecGeometry& fdRecGeometry = fdEvent.GetFdRecGeometry();

  // hybrid
  if (fdEvent.IsHybridEvent()) {
    const double par = 0.213 * degree;
    const double alpha = 2.1;
    const double offset = 0.185 * degree;
    return offset + par * pow(refAngle/angularLength, alpha);
  }

  // pcgf
  if (fdRecGeometry.HasPCGFData()) {
    const double par = 0.266 * degree;
    const double alpha = 0.515;
    const double offset = 0.265 * degree;
    return offset + par * pow(refAngle/angularLength, alpha);
  }

  return std::numeric_limits<double>::infinity();
}


void
LongitudinalXmaxScanner::EstimateXmaxErrors(const RecEvent& theRecEvent)
{
  fLongitudinalScan.clear();
  fErrorAtXmax.clear();

  for (auto iEye = theRecEvent.EyesBegin(); iEye != theRecEvent.EyesEnd(); ++iEye) {

    fCalorimetricEnergy = iEye->GetFdRecShower().GetCalorimetricEnergy() * eV;
    if (fCalorimetricEnergy > 0) {

      GetTelescopeProperties(iEye->GetEyeId());

      if (FillLightFactors(*iEye)) {
        CalculateXmaxUncertainties(*iEye);
        continue;
      }
    }

    // fill empty errors if reach here
    fLongitudinalScan.push_back(LongitudinalScan());
    fErrorAtXmax.push_back(sqrt(kMaxVariance));

  }
}


bool
LongitudinalXmaxScanner::FillLightFactors(const FDEvent& fdEvent)
{
  fDepth.clear();
  fGeomVariance.clear();
  fTime.clear();
  fViewingAngle.clear();
  fDeltaAngle.clear();
  fDeltaTime.clear();
  fPhotoElectronFactor.clear();
  fNoiseVariance.clear();
  fCloseToBoundary.clear();
  fOutsideTelescope.clear();

  det::Detector& detector = det::Detector::GetInstance();
  CoordinateSystemPtr referenceCS = detector.GetReferenceCoordinateSystem();

  const FdRecPixel& recPixel = fdEvent.GetFdRecPixel();
  const double maxRMS = 60; // protect against crazy values
  const int nAverageZetaPix = 2;
  const double averagePhotonVariance =
    pow(fmin(maxRMS, recPixel.GetMeanPixelRMS()), 2) * nAverageZetaPix;
  const double averagePhotoElectronVariance =
    averagePhotonVariance * pow(fPhotonToPhotoElectron, 2);
  fPhotoElectronBGRMS = sqrt(averagePhotoElectronVariance/nAverageZetaPix);

  const unsigned int eyeId = fdEvent.GetEyeId();
  const FdRecShower& recShower = fdEvent.GetFdRecShower();
  double zeta = fdEvent.GetFdRecApertureLight().GetZeta();
  if (zeta <= 0) {
    INFO("zeta not filled!!!");
    zeta = 1.4*degree;
  }
  const TVector3& adstCore = recShower.GetCoreSiteCS();
  const TVector3& adstAxis = recShower.GetAxisSiteCS();
  const Vector axis(adstAxis.X(), adstAxis.Y(), adstAxis.Z(), referenceCS);
  const Point core(adstCore.X()*m, adstCore.Y()*m, adstCore.Z()*m, referenceCS);
  const fdet::Eye& detEye = detector.GetFDetector().GetEye(eyeId);
  const Point& eyePos  = detEye.GetPosition();
  const Vector coreEyeVec = core - eyePos;

  const FdRecGeometry& theGeometry = fdEvent.GetFdRecGeometry();
  const double t0 = theGeometry.GetT0();
  const double chi0 = theGeometry.GetChi0();
  const double rp = theGeometry.GetRp();

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

  try {
    atmo.InitSlantProfileModel(core, axis, 30.*g/cm2);
  } catch (InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
    INFO(e.GetMessage());
    return false;
  }
  const ProfileResult& distVsSlantDepth = atmo.EvaluateDistanceVsSlantDepth();
  const double maxDepth = distVsSlantDepth.MaxX();
  const double minDepth = distVsSlantDepth.MinX();

  const double xLow = recShower.GetXFOVMin() * g/cm2;
  const double xUp = recShower.GetXFOVMax() * g/cm2;

  const int nBins = int((xUp - xLow) / fDepthBinWidth);
  const double dX = (xUp - xLow) / nBins;

  double currX = xLow;
  for (int i = 0; i < nBins+1; ++i) {

    fDepth.push_back(currX);

    const double thisX = currX - 0.5*dX;
    const double nextX = currX + 0.5*dX;

    if (thisX >= minDepth && thisX <= maxDepth &&
        nextX >= minDepth && nextX <= maxDepth) {

      const Point thisPoint = core - axis * distVsSlantDepth.Y(thisX);
      const Point nextPoint = core - axis * distVsSlantDepth.Y(nextX);
      const Point currPoint = core - axis * distVsSlantDepth.Y(currX);
      const Vector thisVec = thisPoint - eyePos;
      const Vector nextVec = nextPoint - eyePos;
      const Vector currVec = currPoint - eyePos;
      const double chi_i = Angle(currVec, coreEyeVec);
      const double t_i = t0 + rp/kSpeedOfLight * std::tan(0.5*(chi0 - chi_i));
      fTime.push_back(t_i);
      fViewingAngle.push_back(chi0 - chi_i);
      const double chi_1 = Angle(thisVec, coreEyeVec);
      const double chi_2 = Angle(nextVec, coreEyeVec);
      const double t_1 = t0 + rp/kSpeedOfLight * std::tan(0.5*(chi0 - chi_1));
      const double t_2 = t0 + rp/kSpeedOfLight * std::tan(0.5*(chi0 - chi_2));
      fDeltaAngle.push_back(fabs(chi_1 - chi_2));
      fDeltaTime.push_back(fabs(t_1 - t_2));
      
      const double lightFactor =
        CalculateLightFactor(thisPoint, nextPoint, eyeId);
      const double timeLength =
        CalculateTimeLength(thisPoint, nextPoint, eyeId);
      const auto telAndPixId = GetTelescopeAndPixelId(currVec, eyeId);
      const auto telId = telAndPixId[0];
      if (telId) {
        fCloseToBoundary.push_back(IsNearBorder(currVec, eyeId, telId, zeta));
        fOutsideTelescope.push_back(false);
        const double area = detEye.GetTelescope(telId).GetDiaphragmArea();
        const double photoElectronFactor =
          lightFactor * area * fPhotonToPhotoElectron;
        fPhotoElectronFactor.push_back(photoElectronFactor);
        if (fUseBGLoop) {
          const fdet::Telescope& detTel = detEye.GetTelescope(telId);
          const fdet::Channel& detChannel = detTel.GetChannel(telAndPixId[1]);
          const double adcVar = detChannel.GetADCVariance();
          const fdet::Pixel& detPixel = detTel.GetPixel(telAndPixId[1]);
          const double calConst =
            detPixel.GetEndToEndCalibrationAtReferenceWavelength();
          const double photonVar = adcVar * pow(calConst, 2);
          const double var = pow(fPhotonToPhotoElectron, 2) * photonVar;
          const double timeBin = detChannel.GetFADCBinSize();
          const double nTimeBins = timeLength / timeBin;
          fNoiseVariance.push_back(nTimeBins * var);
          /*
          cout << detEye.GetId() << " " << detTel.GetId() << " "
               << timeBin << " " << timeLength << " "
               << calConst << " "
               << sqrt(var / pow(timeBin, 2))
               << " " << averagePhotoElectronVariance/2 << endl;
          */
        } else {
          // wrong for HECO and HEAT, but kept for historical reasons (ICRC17)
          const double timeBin = 100*nanosecond;
          fNoiseVariance.push_back(timeLength/timeBin *
                                   averagePhotoElectronVariance);
        }
      } else {
        fCloseToBoundary.push_back(true);
        fOutsideTelescope.push_back(true);
        fPhotoElectronFactor.push_back(0);
        fNoiseVariance.push_back(0);
      }
    } else {
      fPhotoElectronFactor.push_back(-1);
      fTime.push_back(-1);
      fViewingAngle.push_back(-1);
      fDeltaAngle.push_back(-1);
      fDeltaTime.push_back(-1);
      fNoiseVariance.push_back(999);
      fCloseToBoundary.push_back(true);
      fOutsideTelescope.push_back(true);
    }
    
    currX += dX;
  }
  
  return true;
}


vector<double>
LongitudinalXmaxScanner::GetChangedDepth(const FDEvent& fdEvent,
                                         const double nSigRp,
                                         const double nSigChi0,
                                         const double nSigT0,
                                         const double nSigPhi,
                                         const double nSigTheta)
{
  const unsigned int eyeId = fdEvent.GetEyeId();
  const det::Detector& detector = det::Detector::GetInstance();
  const atm::Atmosphere& atmo = detector.GetAtmosphere();
  const fdet::Eye& detEye = detector.GetFDetector().GetEye(eyeId);

  const FdRecGeometry& theGeometry = fdEvent.GetFdRecGeometry();

  const double t0 = theGeometry.GetT0() + nSigT0 * theGeometry.GetT0Error();
  const double chi0 = theGeometry.GetChi0() + nSigChi0 * theGeometry.GetChi0Error();
  const double rp = theGeometry.GetRp() + nSigRp * theGeometry.GetRpError();

  const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();

  const double sdpTheta = theGeometry.GetSDPTheta() + nSigTheta * theGeometry.GetSDPThetaError();
  const double sdpPhi = theGeometry.GetSDPPhi() + nSigPhi * theGeometry.GetSDPPhiError();
  const Vector sdp = Vector(1, sdpTheta, sdpPhi, eyeCS, Vector::kSpherical);

  const Vector vertical(0,0,1, eyeCS);
  const Vector horizontalInSDP = Normalized(Cross(sdp, vertical));

  const Transformation rot(Transformation::Rotation(-chi0, sdp, eyeCS));
  const Vector axis(rot * horizontalInSDP);
  const Vector coreEyeVec = rp / sin(kPi - chi0) * horizontalInSDP;
  const Point core = detEye.GetPosition() + coreEyeVec;

  const ProfileResult* slantDepthProfile = nullptr;
  try {
    // less precision for error propagation, see GAP-2007-007, Fig. 5
    atmo.InitSlantProfileModel(core, axis, 100*g/cm2);
    slantDepthProfile = &atmo.EvaluateSlantDepthVsDistance();
  } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException&) {
    return vector<double>();
  }

  const double minDistance = slantDepthProfile->MinX();
  const double maxDistance = slantDepthProfile->MaxX();

  vector<double> depthVec;
  for (unsigned int i = 0; i < fTime.size(); ++i) {
    const double t_i = fTime[i];
    double depth = -1*g/cm2;
    if (t_i > 0) {
      const double chi_i = chi0 - 2*atan(kSpeedOfLight / rp * (t_i - t0));
      const double distance = -coreEyeVec.GetMag() * sin(chi_i) / sin(chi0 - chi_i);
      if (distance > minDistance && distance < maxDistance)
        depth = slantDepthProfile->Y(distance);
    }
    depthVec.push_back(depth);
  }

  return depthVec;
}


void
LongitudinalXmaxScanner::PropagateGeometryUncertainty(const FDEvent& fdEvent)
{
  fXmaxGeomVar =  kMaxVariance;

  vector<double> depthDefault = GetChangedDepth(fdEvent,0,0,0,0,0);
  vector<double> depthRp1 = GetChangedDepth(fdEvent,1,0,0,0,0);
  vector<double> depthRp2 = GetChangedDepth(fdEvent,-1,0,0,0,0);
  vector<double> depthChi1 = GetChangedDepth(fdEvent,0,1,0,0,0);
  vector<double> depthChi2 = GetChangedDepth(fdEvent,0,-1,0,0,0);
  vector<double> depthT01 = GetChangedDepth(fdEvent,0,0,1,0,0);
  vector<double> depthT02 = GetChangedDepth(fdEvent,0,0,-1,0,0);
  vector<double> depthPhi1 = GetChangedDepth(fdEvent,0,0,0,1,0);
  vector<double> depthPhi2 = GetChangedDepth(fdEvent,0,0,0,-1,0);
  vector<double> depthTheta1 = GetChangedDepth(fdEvent,0,0,0,0,1);
  vector<double> depthTheta2 = GetChangedDepth(fdEvent,0,0,0,0,-1);

  // check for InclinedModel exceptions
  const unsigned int nDepth = fDepth.size();
  bool allValid = (depthDefault.size() == nDepth &&
                   depthRp1.size() == nDepth &&
                   depthRp2.size() == nDepth &&
                   depthChi1.size() == nDepth &&
                   depthChi2.size() == nDepth &&
                   depthT01.size() == nDepth &&
                   depthT02.size() == nDepth &&
                   depthPhi1.size() == nDepth &&
                   depthPhi2.size() == nDepth &&
                   depthTheta1.size() == nDepth &&
                   depthTheta2.size() == nDepth);

  if (!allValid) {
    INFO(" GetChangedDepth failed!");
    for (unsigned int i = 0; i < nDepth; ++i)
      fGeomVariance.push_back(kMaxVariance);
    return;
  }

  // correlation coefficients
  const FdRecGeometry& theGeometry = fdEvent.GetFdRecGeometry();
  const double rhoRT = theGeometry.GetRpT0Correlation();
  const double rhoRC = theGeometry.GetChi0RpCorrelation();
  const double rhoCT = theGeometry.GetChi0T0Correlation();
  const double rhoPT = theGeometry.GetSDPThetaPhiCorrelation();

  for (unsigned int i = 0; i < nDepth; ++i) {
    if (fDepth[i] > 0) {

      const double rpErr = AverageErr(depthDefault[i],
                                      depthRp1[i],
                                      depthRp2[i]);
      const double t0Err = AverageErr(depthDefault[i],
                                      depthT01[i],
                                      depthT02[i]);

      const double chi0Err = AverageErr(depthDefault[i],
                                        depthChi1[i],
                                        depthChi2[i]);

      const double phiErr = AverageErr(depthDefault[i],
                                       depthPhi1[i],
                                       depthPhi2[i]);

      const double thetaErr = AverageErr(depthDefault[i],
                                         depthTheta1[i],
                                         depthTheta2[i]);

      const double timeFitVariance = pow(rpErr,2)
                                     + pow(t0Err,2)
                                     + pow(chi0Err,2)
                                     + 2 * rpErr * t0Err * rhoRT
                                     + 2 * rpErr * chi0Err * rhoRC
                                     + 2 * t0Err * chi0Err * rhoCT;

      const double sdpFitVariance = pow(phiErr,2) + pow(thetaErr,2)
                                    + 2 * phiErr * thetaErr * rhoPT;

      // chi2 scale factors as in FdProfileReconstructor
      const double scaleFac =
        theGeometry.GetTimeFitChi2() / theGeometry.GetTimeFitNdof();
      fGeomVariance.push_back(scaleFac*timeFitVariance + sdpFitVariance);

    } else {
      fGeomVariance.push_back(kMaxVariance);
    }
  }

  if (nDepth < 2)
    return;

  // uncertainty at Xmax
  const double xLow = fDepth[0];
  const double xUp = fDepth[nDepth-1];
  const double xMax = fXmax;
  const double dX = (xUp - xLow) / (nDepth - 1);
  const int index = int((xMax - xLow) / dX);
  if (index >= 0 && index < int(nDepth) - 1) {
    const double deltaX = xMax - xLow - index*dX;
    const double z1 = fGeomVariance[index];
    const double z2 = fGeomVariance[index+1];
    fXmaxGeomVar = z1 + deltaX*(z2 - z1)/dX;
  } else if (index == int(nDepth) - 1)
    fXmaxGeomVar = fGeomVariance[index];
}


void
LongitudinalXmaxScanner::CalculateXmaxUncertainties(const FDEvent& fdEvent)
{
  const FdRecShower& theShower = fdEvent.GetFdRecShower();
  const FdRecApertureLight& theAp = fdEvent.GetFdRecApertureLight();
  fXmax = theShower.GetXmax()*g/cm2;
  fShowerAngularLength = theAp.GetMaxAngle() - theAp.GetMinAngle();

  // calculate geometry uncertainty
  PropagateGeometryUncertainty(fdEvent);

  // uncertainty at Xmax
  fXmaxProfVar = PropagateProfileUncertainty(fdEvent,fXmax);
  fXmaxAngularLength = fAngularLength;
  fXmaxTotErr = CalculateTotalError(fXmaxGeomVar,fXmaxProfVar,fdEvent);

  if (fVerbosity > 1)
    PrintCurrentVariables(fdEvent);

  LongitudinalScan longScan;
  for (unsigned int i = 0; i < fDepth.size(); ++i) {

    const double X = fDepth[i];
    const double geomVar = fGeomVariance[i];
    const double profVar = PropagateProfileUncertainty(fdEvent, X);
    const double totErr = CalculateTotalError(geomVar, profVar, fdEvent);

    longScan.AddScanResult(X, totErr, fMinViewingAngle);

    if (fVerbosity > 0)
      cout << " X= " << X/gcm2 << ' '
           << sqrt(profVar)/gcm2 << ' '
           << sqrt(fGeomVarScaleFac*geomVar)/gcm2 << ' '
           << totErr/gcm2 << ' ' << theShower.GetXmaxError() << ' '
           << fAngularLength/degree
           << endl;
  }

  fLongitudinalScan.push_back(longScan);
  fErrorAtXmax.push_back(fXmaxTotErr);
}


double
LongitudinalXmaxScanner::PropagateProfileUncertainty(const FDEvent& fdEvent,
                                                     const double Xmax)
{
  const FdRecShower& theShower = fdEvent.GetFdRecShower();

  // chi2 scale factor as in FdProfileReconstructor
  const double scaleFac = theShower.GetGHChi2()/theShower.GetGHNdf();

  using namespace evt;
  const GaisserHillas4Parameter::ShapeParameter x0(gh::eX0, kAverageX0);
  const GaisserHillas4Parameter::ShapeParameter lambda(gh::eLambda, kAverageLambda);
  GaisserHillas4Parameter ghFunction(1, Xmax, x0, lambda);

  const double ghIntegral = ghFunction.GetIntegral();
  const double dEdXmax = fCalorimetricEnergy/ghIntegral;

  vector<double> npeVec;
  vector<double> ghFunc;
  vector<double> dEdXErr;

  for (unsigned int i = 0; i < fDepth.size(); ++i) {

    const double dEdX = dEdXmax * ghFunction.Eval(fDepth[i]);
    const double npe = dEdX * fPhotoElectronFactor[i];
    const double variance = npe*(1 + fGainVariance) + fNoiseVariance[i];
    const double dEdX_Err = sqrt(variance) / fPhotoElectronFactor[i];
    dEdXErr.push_back(dEdX_Err);
    ghFunc.push_back(dEdX);
    npeVec.push_back(npe);

  }

  vector<unsigned short> triggerVec = CalculatePixelTrigger(npeVec);

  // "deactivate" bins without trigger by setting large uncertainty
  // and fill some variables
  fTrackMin = 0;
  fTrackMax = 0;
  fMinViewingAngle = 0;
  fAngularLength = 0;
  fBins = 0;

  double angSum = 0;
  for (unsigned int i = 0; i < fDepth.size(); ++i) {
    // MU: triggerVec can not be commented out, otherwise MVA will be wrong
    if (!triggerVec[i] || fCloseToBoundary[i] || fOutsideTelescope[i])
      dEdXErr[i] = ghFunc[i]*100;
    else {
      if (!fBins || fDepth[i] < fTrackMin)
        fTrackMin = fDepth[i];
      if (!fBins || fDepth[i] > fTrackMax)
        fTrackMax = fDepth[i];
      if (!fBins || fMinViewingAngle > fViewingAngle[i])
        fMinViewingAngle = fViewingAngle[i];
      angSum += fDeltaAngle[i];
      ++fBins;
    }
  }

  fAngularLength = angSum;

  double profVar = kMaxVariance;
  if (fBins > 0 && angSum > 0) {

    profVar = EstimateXmaxVariance(fDepth, ghFunc, dEdXErr,
                                   dEdXmax, Xmax,
                                   kAverageX0, kAverageLambda);
    profVar *= scaleFac;
  }

  if (fVerbosity > 0) {
    if (Xmax == fXmax) {
      for (unsigned int i = 0; i < fDepth.size(); ++i) {
        cout << "-- X=" << fDepth[i]/g*cm2
             << " dEdX=" << ghFunc[i]/(PeV/g*cm2)
             << " s(dEdX)=" << dEdXErr[i]/(PeV/g*cm2)
             << " relErr=" << dEdXErr[i]/ghFunc[i]
             << " trig=" << triggerVec[i]
             << " telBound=" << fCloseToBoundary[i]
             << " isInTelescope=" << !fOutsideTelescope[i]
             << " dAng=" << fDeltaAngle[i]/degree
             << " dT=" << fDeltaTime[i]/(100*ns)
             << " viewAngle= " << fViewingAngle[i]/degree
             << '\n';
      }
      cout << flush;
    }
  }

  return profVar;
}


void
LongitudinalXmaxScanner::PrintCurrentVariables(const FDEvent& fdEvent)
  const
{
  const FdRecShower& theShower = fdEvent.GetFdRecShower();
  const FdRecApertureLight& theAp = fdEvent.GetFdRecApertureLight();
  cout << " check "
       << fdEvent.GetEyeId() << ' '
       << fdEvent.GetRunId() << ' '
       << fdEvent.GetEventId() << ' '
       << fBins << ' '
       << log10(theShower.GetEnergy()) << ' '
       << fXmax/g*cm2 << ' '
       << theShower.GetXmaxError() << ' '
       << sqrt(fXmaxGeomVar+fXmaxProfVar)/g*cm2 << ' '
       << sqrt(fXmaxProfVar)/g*cm2 << ' '
       << sqrt(fXmaxGeomVar)/g*cm2 << ' '
       << fTrackMin/g*cm2 << ' '
       << fTrackMax/g*cm2 << ' '
       << theShower.GetXTrackMin() << ' '
       << theShower.GetXTrackMax() << ' '
       << fAngularLength/degree << ' '
       << fMinViewingAngle/degree << ' '
       << theAp.GetMinAngle()/degree << ' '
       << fGeomVarScaleFac << endl;
}


// returns photons at aperture per dE/dX per area
double
LongitudinalXmaxScanner::CalculateLightFactor(const Point& pos1,
                                              const Point& pos2,
                                              const unsigned int eyeId)
  const
{
  det::Detector& detector = det::Detector::GetInstance();
  const atm::Atmosphere& atmo = detector.GetAtmosphere();
  const fdet::Eye& detEye = detector.GetFDetector().GetEye(eyeId);
  const Point& eyePos = detEye.GetPosition();

  const Point meanPos = pos1 - 0.5*(pos1 - pos2);
  const Vector Mpos = 0.5*(pos1 - pos2);
  const double lambdaRef = detector.GetFDetector().GetReferenceLambda();
  // FIXME tel 1
  const utl::TabulatedFunction& eff =
    detEye.GetTelescope(1).GetMeasuredRelativeEfficiency();
  const double lambdaEffMin = eff.GetX(0);
  const double lambdaEffMax = eff.GetX(eff.GetNPoints()-1);
  const double epsilonRef = eff.Y(lambdaRef);

  const vector<double>& fluoLambda =
    atmo.GetWavelengths(Atmosphere::eFluorescence);

  // Mie attenuation
  const AttenuationResult mAtt =
    atmo.EvaluateMieAttenuation(eyePos,meanPos,fluoLambda);
  const utl::TabulatedFunctionErrors& mieAtt = mAtt.GetTransmissionFactor();
  // Rayleigh attenuation
  const AttenuationResult rAtt =
    atmo.EvaluateRayleighAttenuation(eyePos,meanPos,fluoLambda);
  const utl::TabulatedFunctionErrors& rayAtt = rAtt.GetTransmissionFactor();

  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
  const double height = wgs84.PointToLatitudeLongitudeHeight(meanPos).get<2>();

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

  double effYTsum_F = 0;
  for (unsigned int i = 0; i < fluoLambda.size(); ++i) {
    if (fluoLambda[i] >= lambdaEffMin &&
        fluoLambda[i] <= lambdaEffMax) {
      const double lambda = fluoLambda[i];
      const double T_i = mieAtt.GetY(i) * rayAtt.GetY(i);
      const double Y_i = fluoyield_tab.Y(lambda);
      effYTsum_F += eff.Y(lambda) / epsilonRef * Y_i * T_i;
    }
  }

  // calculate Cherenkov part

  const vector<double>& cherLambda =
    atmo.GetWavelengths(Atmosphere::eCerenkov);

  // ---- attenuation to telescope for cherenkov

  // Mie attenuation
  const AttenuationResult mCAtt =
    atmo.EvaluateMieAttenuation(eyePos, meanPos, cherLambda);
  // Rayleigh attenuation
  const AttenuationResult rCAtt =
    atmo.EvaluateRayleighAttenuation(eyePos, meanPos, cherLambda);
  //const utl::TabulatedFunctionErrors& rayCAtt = rAtt.GetTransmissionFactor();// unused. LN.

  /*
    MU: FOV calculation should *not* depend on the Xmax of the event.
        In that sense a hardcoded value is better, but a proper fix
        should re-calculate the Cherenkov yield for each Xmax position
        again
  */

  //const double Xmax = recShower.GetGHParameters().GetXMax();;

  //const double Xmax = 750*g/cm2; // used for shower age for Cherenkov // unused. LN.
   
  //const double slantDepth = Mpos.GetMag();// unused. LN.
  //const double showerAge = utl::ShowerAge(slantDepth, Xmax);// unused. LN.
  //const Point& telPos = detEye.GetTelescope(1).GetPosition();// unused. LN.

  // Cherenkov yield per electron above electronEnergyThreshold
  const double electronEnergyThreshold = 1*MeV;
  atmo.SetCherenkovEnergyCutoff(electronEnergyThreshold);
  //const utl::TabulatedFunction& directCherenkov = // unused. LN.
  //  atmo.EvaluateCherenkovDirect(pos1, pos2, telPos, showerAge);

  // average energy deposit per electron above electronEnergyThreshold
  //const double alpha = utl::EnergyDeposit(showerAge, electronEnergyThreshold);// unused. LN.

  for (unsigned int i = 0; i < cherLambda.size(); ++i) {
    if (cherLambda[i] >= lambdaEffMin &&
        cherLambda[i] <= lambdaEffMax) {
      //const double lambda = cherLambda[i];// unused. LN.
      //const double T_i = mieCAtt.GetY(i) * rayCAtt.GetY(i); // unused. LN.
      //const double Y_i = directCherenkov.GetY(i); // unused. LN.
      //const double epsilon_i = eff.Y(lambda) / epsilonRef; // unused. LN.

      //effYTsum_C += epsilon_i * T_i * Y_i / alpha;
    }
  }

  const double deltaL = (pos1 - pos2).GetMag();
  const double fourPi = 4*kPi;
  const utl::Vector r = meanPos - eyePos;
  const double rSquared = r.GetMag2();

  /*
    MU: I think this should be
  
    lightFac = (effYTsum_F * deltaL / (fourPi * rSquared * dEdX0)) + effYTsum_C;
    
    because afaik EvaluateCherenkovDirect() gives already light flux at 
    telescope

    --> please check, direct C-light contribution commented for now
  */
  
  const double dEdX0 = atmo.GetdEdX0();
  const double totalsum = effYTsum_F;  // MU: FIXME  + effYTsum_C;
  const double lightFac = totalsum * deltaL / (fourPi * rSquared * dEdX0);

  return lightFac;
}


void
LongitudinalXmaxScanner::GetTelescopeProperties(unsigned int eyeId)
{
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = detFD.GetEye(eyeId);
  for (fdet::Eye::TelescopeIterator teliter = detEye.TelescopesBegin();
       teliter != detEye.TelescopesEnd(); ++teliter) {
    const fdet::Telescope& detTel = *teliter;
    fGainVariance = detTel.GetCamera().GetGainVariance();
    // FIXME hardcoded
    fPhotonToPhotoElectron = 0.128;
    break;
  }

  // FIXME hardcoded
  fPixelSize = 1.5*degree;
}


// returns time difference at diaphragm
double
LongitudinalXmaxScanner::CalculateTimeLength(const utl::Point& pos1,
                                             const utl::Point& pos2,
                                             const unsigned int eyeId)
  const
{
  det::Detector& detector = det::Detector::GetInstance();
  const fdet::Eye& detEye = detector.GetFDetector().GetEye(eyeId);
  const Point& eyePos = detEye.GetPosition();

  const Vector pos1Eye = pos1 - eyePos;
  const Vector pos2Eye = pos2 - eyePos;
  const Vector pos1pos2 = pos1 - pos2;

  const double tof1 = pos1Eye.GetMag() / kSpeedOfLight;
  const double tof2 = (pos1pos2.GetMag() + pos2Eye.GetMag()) / kSpeedOfLight;

  return fabs(tof1 - tof2);
}


vector<unsigned short>
LongitudinalXmaxScanner::CalculatePixelTrigger(const vector<double>& npeVec)
  const
{
  vector<unsigned short> triggerVec;

  double angSum = 0;
  double timeSum = 0;
  double npeSum = 0;

  for (unsigned int i = 0; i < npeVec.size(); ++i) {

    if (fDeltaAngle[i] < 0) {
      triggerVec.push_back(0);
      continue;
    }

    timeSum += fDeltaTime[i];
    angSum += fDeltaAngle[i];
    npeSum += npeVec[i];

    if (angSum > fPixelSize) {

      const bool triggered = IsTriggered(timeSum, angSum, npeSum);
      for (unsigned int j = triggerVec.size(); j <= i; ++j)
        triggerVec.push_back(triggered);

      angSum = 0;
      timeSum = 0;
      npeSum = 0;

    }
  }

  if (timeSum != 0) {
    const bool triggered = IsTriggered(timeSum, angSum, npeSum);
    for (unsigned int j = triggerVec.size(); j < npeVec.size(); ++j)
      triggerVec.push_back(triggered);
  }

  return triggerVec;
}


bool
LongitudinalXmaxScanner::IsTriggered(const double dT, const double dAng, const double npe)
  const
{
  const double thresholdSlope = 15.;
  const double threshold = thresholdSlope*fPhotoElectronBGRMS;
  const int nBoxCar = 10;
  const double timeBinLength = 100*ns;

  const double timeInPixel = fPixelSize / dAng * dT;
  const double signalInPixel = fPixelSize / dAng * npe;
  const double signalBoxTime = fmin(timeInPixel, nBoxCar*timeBinLength);

  return signalInPixel/timeInPixel*signalBoxTime > threshold;
}


double
LongitudinalXmaxScanner::CalculateTotalError(const double geomVar,
                                             const double profVar,
                                             const FDEvent& fdEvent)
{
  if (geomVar == kMaxVariance ||
      profVar == kMaxVariance ||
      fAngularLength < 1*degree)
    return sqrt(kMaxVariance);

  const double xmaxAngleErr = SpaceAngle(fXmaxAngularLength, fdEvent);
  const double currAngleErr = SpaceAngle(fAngularLength, fdEvent);
  // uncomment next two lines to re-scale to ADST track length
  //const double adstAngleErr = SpaceAngle(fShowerAngularLength);
  const double adstFactor = 1; //adstAngleErr/xmaxAngleErr;
  const double currentFactor = currAngleErr / xmaxAngleErr;

  fGeomVarScaleFac = pow(adstFactor*currentFactor, 2);

  const double totErr = sqrt(fGeomVarScaleFac*geomVar + profVar);

  return totErr;
}


bool
LongitudinalXmaxScanner::IsNearBorder(const utl::Vector& direction,
                                      const unsigned int eyeId,
                                      const unsigned int telId,
                                      const double zeta)
  const
{
  det::Detector& detector = det::Detector::GetInstance();
  const fdet::Eye& detEye = detector.GetFDetector().GetEye(eyeId);
  const fdet::Telescope& detTel = detEye.GetTelescope(telId);

  for (fdet::Telescope::OutOfBorderPixelsIterator iter = detTel.OutOfBorderPixelsBegin();
       iter!= detTel.OutOfBorderPixelsEnd(); ++iter) {
    const double angle = Angle(direction, *iter);
    if (angle < zeta)
      return true;
  }

  return false;
}


std::array<unsigned int, 2>
LongitudinalXmaxScanner::GetTelescopeAndPixelId(const utl::Vector& direction,
                                                const unsigned int eyeId)
  const
{
  det::Detector& detector = det::Detector::GetInstance();
  const fdet::Eye& detEye = detector.GetFDetector().GetEye(eyeId);

  unsigned int iPix = 0;
  unsigned int iTel = 0;
  double min = std::numeric_limits<double>::infinity();

  for (fdet::Eye::TelescopeIterator tIt = detEye.TelescopesBegin();
       tIt != detEye.TelescopesEnd(); ++tIt) {
    for (unsigned int i = tIt->GetFirstPixelId(), n = tIt->GetLastPixelId();
         i <= n; ++i) {
      const double diff = Angle(direction, tIt->GetPixel(i).GetDirection());
      if (diff < min) {
        iTel = tIt->GetId();
        iPix = i;
        min = diff;
      }
    }
  }

  return { iTel, iPix };
}


double
LongitudinalXmaxScanner::EstimateXmaxVariance(const vector<double>& x,
                                              const vector<double>& ghFunc,
                                              const vector<double>& eY,
                                              const double nMax,
                                              const double xMax,
                                              const double x0,
                                              const double lambda)
  const
{
  if (xMax < x0)
    return xMax;

  TMatrixD a(4, 4);

  for (unsigned int i = 0; i < x.size(); ++i) {

    // skip point outside atmosphere (flagged as -1.)
    if (eY[i] < 0)
      continue;

    // inverse of squared relative GH errors
    const double tmp = eY[i] / ghFunc[i];
    const double invSqRelErr = 1 / pow(tmp, 2);

    const double logTerm = log((x0 - x[i]) / (x0 - xMax));
    // dGH/dNmax
    const double ghRelDeriv1 = 1 / nMax;
    // dGH/dXmax
    const double ghRelDeriv2 = logTerm / lambda;
    const double nominatorX0 = (x0 - x[i]) * logTerm + (x[i] - xMax);
    const double nominatorLambda = (x0 - xMax) * logTerm + (x[i] - xMax);
    // dGH/dX0
    const double ghRelDeriv3 = nominatorX0 / lambda / (x[i] - x0);
    // dGH/dlambda
    const double ghRelDeriv4 = nominatorLambda / (lambda * lambda);

    const double deriv[4] = { ghRelDeriv1, ghRelDeriv2, ghRelDeriv3, ghRelDeriv4 };

    for (int i = 0; i < 4; ++i)
      for (int j = 0; j <= i; ++j)
        a[i][j] += invSqRelErr * deriv[i] * deriv[j];

  }

  a[2][2] += 1 / kX0Var;
  a[3][3] += 1 / kLambdaVar;

  a[0][1] = a[1][0];
  a[0][2] = a[2][0];
  a[0][3] = a[3][0];
  a[1][2] = a[2][1];
  a[1][3] = a[3][1];
  a[2][3] = a[3][2];

  double determinant;
  a.InvertFast(&determinant);

  if (determinant <= 0) {
    ostringstream info;
    info << "error matrix is singular, det="
         << determinant
         << ", XmaxError = " << sqrt(a[1][1])/gcm2 << endl;
    INFO(info.str());
    return kMaxVariance;
  }

  return a[1][1];
}