FdProfileReconstructor.cc

Go to the documentation of this file.

#include "FdProfileReconstructor.h"
#include "EnergyFitter.h"
#include "CherenkovFluorescenceMatrix.h"
#include "RootCFMatrixOutput.h"
#include "OpticalHalo.h"
#include "matrixInversionExceptions.h"
#include "diagonalMatrix.h"
#include "lowerTriangularMatrix.h"
#include "upperTriangularMatrix.h"

#include "UncorrelatedUncertainties.h"

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>

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

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

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

#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/UTMPoint.h>
#include <utl/Transformation.h>

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

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

#include <boost/io/ios_state.hpp>

#include <iostream>
#include <string>
#include <sstream>
#include <sstream>
#include <vector>

using namespace FdProfileReconstructorKG;
using namespace fwk;
using namespace utl;
using namespace fevt;
using namespace fdet;
using namespace evt;
using namespace atm;
using namespace oBLAS;
using namespace std;




/******************************************************************\
 *
 *  constructor
 *
\******************************************************************/
FdProfileReconstructor::FdProfileReconstructor() :
  fInclinedModelDepthBinning(10*g/cm2),
  fYieldRefit(false),
  fAlwaysGuessXMaxFromAperture(false),
  fIsMultiEyeEvent(false),
  fXmax(0.),
  fEnergyFitter(NULL),
  fRootCFMOutput(NULL)
{
}


/******************************************************************\
 *
 *  destructor
 *
\******************************************************************/
FdProfileReconstructor::~FdProfileReconstructor()
{
}


/******************************************************************\
 *
 *  initialize and read XML cards
 *
\******************************************************************/
fwk::VModule::ResultFlag FdProfileReconstructor::Init()
{

  CentralConfig* cc = CentralConfig::GetInstance();

  Branch topB =
    cc->GetTopBranch("FdProfileReconstructor");

  Branch profCalcB = topB.GetChild("profileCalculation");

  // verbosity flag controls general output (-1: no output)
  profCalcB.GetChild("verbosity").GetData(fVerbosity);

  // energy cut off for number of electrons
  profCalcB.GetChild("energyCutoff").GetData(fEcut);

  // scaling factors for fluorescence and Cherenkov yields
  profCalcB.GetChild("fluoScaleFac").GetData(fFfac);
  profCalcB.GetChild("cherScaleFac").GetData(fCfac);

  // angle above which inclined atmosphere should be used
  double inclinedModelMinZenith;
  profCalcB.GetChild("inclinedModelMinZenith").GetData(inclinedModelMinZenith);
  fInclinedModelMaxCosZenith=cos(inclinedModelMinZenith);

  profCalcB.GetChild("inclinedModelDepthBinning").GetData(fInclinedModelDepthBinning);


  // correction of shower image width
  std::string fluoLDF, dirCherLDF, scatCherLDF, multscatLDF, haloType;

  profCalcB.GetChild("fluorescenceLDF").GetData(fluoLDF);
  profCalcB.GetChild("directCherenkovLDF").GetData(dirCherLDF);
  profCalcB.GetChild("scatteredCherenkovLDF").GetData(scatCherLDF);
  profCalcB.GetChild("multipleScatteringLDF").GetData(multscatLDF);
  profCalcB.GetChild("cherenkovCone").GetData(fCherenkovCone);
  profCalcB.GetChild("opticalHalo").GetData(haloType);

  if ( fluoLDF == "eNone")
    fFluoLDF = CherenkovFluorescenceMatrix::eNoFluoLDF;
  else if ( fluoLDF == "eGora")
    fFluoLDF = CherenkovFluorescenceMatrix::eGoraFluoLDF;
  else {
    ERROR("unknown fluorescence LDF");
    throw utl::NonExistentComponentException();
  }

  if ( dirCherLDF == "eNone")
    fDirCherLDF = CherenkovFluorescenceMatrix::eNoDirCherLDF;
  else if ( dirCherLDF == "eGora")
    fDirCherLDF = CherenkovFluorescenceMatrix::eGoraDirCherLDF;
  else {
    ERROR("unknown direct Cherenkov LDF");
    throw utl::NonExistentComponentException();
  }

  if ( scatCherLDF == "eNone")
    fScatCherLDF = CherenkovFluorescenceMatrix::eNoScatCherLDF;
  else if ( scatCherLDF == "eGora")
    fScatCherLDF = CherenkovFluorescenceMatrix::eGoraScatCherLDF;
  else if ( scatCherLDF == "eExponential")
    fScatCherLDF = CherenkovFluorescenceMatrix::eExpoScatCherLDF;
  else if ( scatCherLDF == "eGiller")
    fScatCherLDF = CherenkovFluorescenceMatrix::eGillerScatCherLDF;
  else {
    ERROR("unknown scattered Cherenkov LDF");
    throw utl::NonExistentComponentException();
  }

  if ( multscatLDF == "eNone")
    fMultScatLDF = CherenkovFluorescenceMatrix::eNoMultScatLDF;
  else if ( multscatLDF == "eRoberts")
    fMultScatLDF = CherenkovFluorescenceMatrix::eRobertsLDF;
  else if ( multscatLDF == "ePekala")
    fMultScatLDF = CherenkovFluorescenceMatrix::ePekalaLDF;
  else {
    ERROR("unknown multiple scattering model");
    throw utl::NonExistentComponentException();
  }

  if ( haloType == "eNone")
    fOpticalHalo = OpticalHalo::eNone;
  else if ( haloType == "eFlasher2005" )
    fOpticalHalo = OpticalHalo::eFlasher2005;
  else if ( haloType == "eFlasher2008" )
    fOpticalHalo = OpticalHalo::eFlasher2008;
  else if ( haloType == "eSpotGroup" )
    fOpticalHalo = OpticalHalo::eSpotGroup;
  else if ( haloType == "eSpotGroup2" )
    fOpticalHalo = OpticalHalo::eSpotGroup2;
  else {
    ERROR("unknown optical halo type");
    throw utl::NonExistentComponentException();
  }

  Branch errPropB = topB.GetChild("errorPropagation");
  // controls print level for EnergyFitter
  errPropB.GetChild("verbosity").GetData(fErrPropVerb);
  // option to propagate geometry errors
  errPropB.GetChild("propagateGeometryErrors").GetData(fPropagateGeometryErrors);
  // option to propagate atmospheric errors
  errPropB.GetChild("propagateAtmUncertainties").GetData(fPropagateAtmUncertainties);
  // option to rescale errors for propagation
  errPropB.GetChild("rescaleErrors").GetData(fRescaleErrors);
  // option to combine data points for faster error propagation
  errPropB.GetChild("maxErrPropPoints").GetData(fmaxErrPropPoints);

  Branch rootB = topB.GetChild("rootOutput");
  rootB.GetChild("saveToRoot").GetData(fSaveToRoot);

  // invisible energy options
  Branch eCalcB = topB.GetChild("invisibleEnergy");
  string compositionOption;
  eCalcB.GetChild("composition").GetData(compositionOption);

  fComposition =
    (compositionOption == "data") ?
      InvisibleEnergy::eData :
      (compositionOption == "mixed") ?
        InvisibleEnergy::eMixedComposition :
        (compositionOption == "proton") ?
          InvisibleEnergy::eProton : InvisibleEnergy::eIron;

  string interactionOption;
  eCalcB.GetChild("interactionModel").GetData(interactionOption);
  if (interactionOption == "SIBYLL21") {
    fInvisibleEnergyModel = InvisibleEnergy::eSIBYLL21;
  } else if (interactionOption == "QGSJET01") {
    fInvisibleEnergyModel = InvisibleEnergy::eQGSJET01;
  } else if (interactionOption == "DATA") {
    fInvisibleEnergyModel = InvisibleEnergy::eDATA;
  } else {
    ostringstream err;
    err << "Unknown hadronic interaction model used: " << interactionOption;
    ERROR(err);
    return eFailure;
  }

  if ((fInvisibleEnergyModel == InvisibleEnergy::eDATA) ^ (fComposition == InvisibleEnergy::eData)) {
    ERROR("inconsistent invisible energy options (data options cannot "
          "be mixed with other options)");
    return eFailure;
  }

  if (fVerbosity > -1) {
    // ----- info output
    std::string idstr("$Id$");

    std::istringstream is(idstr);
    std::string version;
    is >> version >> version >> version;

    std::ostringstream info;
    info << " Version: " << version
        << "\n Parameters:\n"
        << "\t energy cutoff:  "
        << fEcut/MeV << " MeV \n"
        << "\t yield scale factors: F="
        << fFfac << " C=" << fCfac << "\n"
        << "\t using inclined profile model above "
        << inclinedModelMinZenith/degree << " deg. (dX="
        << fInclinedModelDepthBinning/g*cm2 << "g/cm^2)\n"
        << "\t geometric error propagation is "
        << (fPropagateGeometryErrors?"on":"off") << "\n"
        << "\t atmospheric error propagation is "
        << (fPropagateAtmUncertainties?"on":"off") << "\n"
        << "\t lateral fluorescence light model: " << fluoLDF  << "\n"
        << "\t lateral direct Cherenkov light model: " << dirCherLDF  << "\n"
        << "\t lateral scattered Cherenkov light model: " << scatCherLDF
        << (fScatCherLDF == CherenkovFluorescenceMatrix::eGillerScatCherLDF?
        (fCherenkovCone?" (with Cherenkov cone)":" (without Cherenkov cone)"):
        "") << "\n"
        << "\t optical halo: " << haloType  << "\n"
        << "\t multiple scattering model: " << multscatLDF  << "\n"
        << "\t invisible energy model: "
        << interactionOption << "/"
        << compositionOption << "\n";

    INFO(info);
  }


  if (fSaveToRoot)
    fRootCFMOutput = new RootCFMatrixOutput();

  // create an EnergyFitter and initialize it

  fEnergyFitter=new EnergyFitter();
  if(fEnergyFitter->Init()) {
    return eSuccess;
  }
  else
    return eFailure;

}

/******************************************************************\
 *
 *  things to do at end of run
 *
\******************************************************************/
fwk::VModule::ResultFlag FdProfileReconstructor::Finish()
{

  DeleteCFMatrix();

  delete fEnergyFitter; fEnergyFitter=NULL;
  delete fRootCFMOutput;fRootCFMOutput=NULL;

  if(fEnergyFitter->Finish())
    return eSuccess;
  else
    return eFailure;
}


/******************************************************************\
 *
 *  steering routine for profile calculation for each eye
 *
\******************************************************************/
fwk::VModule::ResultFlag
FdProfileReconstructor::Run(evt::Event& event)
{
  if (event.HasSimShower() &&
      (fOpticalHalo != OpticalHalo::eNone || fMultScatLDF)) {
    WARNING("############################################################################\n"
            "############################################################################\n"
            "# really reconstruct simulated shower with\n"
            "# multiple scattering and/or spot halo ???\n"
            "# note: none of these effects is currently simulated,\n"
            "# i.e. unless you know what you are doing, the profile reconstruction\n"
            "# will be probably wrong...)\n"
            "############################################################################\n"
            "############################################################################");
  }

  // clean up
  if (!fYieldRefit)
    DeleteCFMatrix();

  if (!event.HasFEvent())
    return eSuccess;

  // set cout format (but remember old one)
  boost::io::ios_all_saver save(cout);
  cout.flags(std::ios::scientific);
  cout.precision(3);

  // count eyes with profile reconstruction
  int nEyes = 0;

  if (fVerbosity > -1)
    cout << "\n    --[FdProfileReconstructor]--> " << endl;

  // loop on eyes for profile reconstruction
  fevt::FEvent& fdEvent = event.GetFEvent();
  fIsMultiEyeEvent = (fdEvent.GetNEyes() > 1);
  for (fevt::FEvent::EyeIterator eyeIter = fdEvent.EyesBegin(ComponentSelector::eHasData);
       eyeIter != fdEvent.EyesEnd(ComponentSelector::eHasData);
       ++eyeIter) {

    fevt::Eye& eye = *eyeIter;
    fErrorCalcMode = false;

    // check for geometry and light at aperture
    if (!eye.HasRecData() ||
        !eye.GetRecData().HasLightFlux(fevt::FdConstants::eTotal) )
      continue;

    if (fVerbosity > -1)
      cout << "\n      calculating profile for eye " << eye.GetId() << endl;

    // set default geometry
    SetGeometry(eye, 0,0,0,0,0);
    fGeomVariance.clear();
    fAtmVariance.clear();

    // retrieve the reconstructed data
    if (GetShowerFRecData(event, eye)) {

      CalculateDetEff(eye);

      // calculate the shower geometry, depth and efficiencies
      if (!CalculateGeometryAndDepth(eye))
        continue;

      // calculate the profiles at track
      fXmax = GuessShowerMaximum(event, eye);
      if (!CalculateProfiles(eye))
        continue;

      // fit Gaisser-Hillas and calculate energy
      if (fEnergyFitter->Run(eye, fChFlPtr)) {
        ++nEyes;

        if (fSaveToRoot)
        fRootCFMOutput->Write(event, eye, *fChFlPtr);

        // calculate total errors and save them
        CalculateTotalErrors(eye);

      }

    }

  }

  fRelEff.Clear();
  fTrackPositions.clear();
  fDepth.clear();

  // check if at least one eye was sucessful
  return (nEyes > 0) ? eSuccess : eContinueLoop;
}


/******************************************************************\
 *
 *  get ShowerFRecData from data structure
 *  either from global data from event (if calculated stereo or hybrid)
 *  or from the current eye
 *
\******************************************************************/
bool
FdProfileReconstructor::GetShowerFRecData(evt::Event& /*event*/, fevt::Eye& eye)
{
  // stuff below eventually for global showers (i.e. stereo
  // currently not supported)

  /*
    fShowerFRecData=NULL;
    if (event.HasRecShower()) {
    if(event.GetRecShower().HasFRecShower()) {
    fShowerFRecData=&event.GetRecShower().GetFRecShower();
    // set the cut off value in data structure
    fShowerFRecData->SetEnergyCutoff(fEcut);
    }
    }
  */

  if (fShowerFRecData == NULL &&
      eye.GetRecData().HasFRecShower()) {
    fShowerFRecData = &eye.GetRecData().GetFRecShower();
    fShowerFRecData->SetEnergyCutoff(fEcut);
  }

  return fShowerFRecData != NULL;
}


/******************************************************************\
 *
 *  calculation of
 *
 *                - dEdX at track
 *                - N_e  at track
 *                - fluorescence photons at aperture
 *                - scattered Cherenkov photons at aperture
 *                - direct Cherenkov photons at aperture
 *
\******************************************************************/
bool
FdProfileReconstructor::CalculateProfiles(fevt::Eye& eye)
{
  namespace ublas = boost::numeric::ublas;

  if (fVerbosity > -1)
    cout << "        current shower maximum is " << fXmax/(g/cm2) << endl;

  // eye position
  const Point& eyepos =
    det::Detector::GetInstance().GetFDetector().GetEye(eye).GetPosition();

  // light flux at diaphragm
  const TabulatedFunctionErrors& lightflux =
    eye.GetRecData().GetLightFlux(fevt::FdConstants::eTotal);

  const unsigned int nF = fDepth.size();

  const unsigned int kMinPoints = 5;
  if (fDepth.size() - fOutSideFOV<kMinPoints) {
    if (fVerbosity > -1)
      cout << "        too few points for fit:  nF="
           << fDepth.size()-fOutSideFOV << endl;
    return false;
  }

  // vector of light flux and its covariance
  ublas::vector<double> n(nF);
  diagonalMatrix Vn(nF);
  for (unsigned int i = 0; i < nF ; ++i) {
    if (i < fOutSideFOV) {
      n(i) = 0;
      Vn(i, i) = 0;
    } else {
      const int j = i - fOutSideFOV;
      n(i) = lightflux.GetY(j);
      const double en = lightflux.GetYErr(j);
      Vn(i, i) = en*en;
    }
  }

  const double zeta = eye.GetRecData().GetZeta();
  if (zeta <= 1.e-12) {
    ostringstream msg;
    msg << "############################################################################\n"
           "# Eye " << eye.GetId() << "'s Zeta angle is " << zeta << ". There was likely a problem\n"
           "# in the aperture-light finding. Check your configuration! Skipping Eye...\n"
           "############################################################################";
    WARNING(msg);
    return false;
  }

  double fluoScaleFac = fFfac;
  double cherScaleFac = fCfac;
  if (fluoScaleFac < 0) {
    fluoScaleFac *= -1;
    cherScaleFac *= -1;
  }

  // solve n=CF*dEdX
  if (!fYieldRefit) {
    fChFlPtr =
      new CherenkovFluorescenceMatrix(fXmax, eyepos, fTrackPositions,
                                      fRelEff, fDepth, fEcut, zeta,
                                      cherScaleFac, fluoScaleFac,
                                      fFluoLDF, fDirCherLDF,
                                      fScatCherLDF, fMultScatLDF,
                                      fCherenkovCone, fOpticalHalo);
    fChFlPtrMap[eye.GetId()] = fChFlPtr;
  } else {
    std::map<int, CherenkovFluorescenceMatrix*>::iterator currentCFMatrix =
      fChFlPtrMap.find(eye.GetId());
    if (currentCFMatrix != fChFlPtrMap.end()) {
      fChFlPtr = currentCFMatrix->second;
      fChFlPtr->SetYieldFactors(fluoScaleFac, cherScaleFac);
    } else {
      cout << " could not find CFM for eye " << eye.GetId()
           << "\n map contains " << fChFlPtrMap.size()
           << " entries :";
      for (std::map<int, CherenkovFluorescenceMatrix*>::iterator iter = fChFlPtrMap.begin();
           iter != fChFlPtrMap.end(); ++iter)
        cout<< iter->first << " ";
      cout << endl;
      return false;
    }
  }

  CherenkovFluorescenceMatrix& ChFl = *fChFlPtr;

  if (fVerbosity > -1)
    cout << "        inverting (" << fDepth.size() << "x" <<fDepth.size()
         << ") CherenkovFluorescenceMatrix ... " << endl;

  const lowerTriangularMatrix& CF = *ChFl.GetCherenkovFluorescenceMatrix();
  const lowerTriangularMatrix* CFinvPtr = nullptr;
  try {
    CFinvPtr = CF.GetInverse();
  } catch (SingularMatrixException& s) {
    cerr << "       FdProfileReconstructor::CalculateProfiles - Error:"
         << " SingularMatrixException for CF" << endl;
    return false;
  }
  const lowerTriangularMatrix& CFinv = *CFinvPtr;

  // calculate covariance matrix

  if (fVerbosity > -1)
    cout << "        covariance matrix ... " << endl;

  // calculate diagonal elements of  CF*V*(CF^T)^(-1)

  std::vector<double> covdEdX(nF);
  for (unsigned int i = 0; i < nF; ++i) {
    double sum = 0;
    for (unsigned int k = i; k < nF; ++k)
      sum += CFinv(k, i) * CFinv(k, i) * Vn(i, i);
    covdEdX[i] = sum;
  }

  // calculate light fractions at aperture

  if (fVerbosity > -1)
    cout << "        profiles at aperture ... \n" << endl;

  // fill profiles

  if (!fErrorCalcMode) {

    InitProfiles(eye);

    fevt::EyeRecData& eyeRecData = eye.GetRecData();
    TabulatedFunctionErrors& dirCher = eyeRecData.GetLightFlux(fevt::FdConstants::eCherDirect);
    TabulatedFunctionErrors& scatMCher = eyeRecData.GetLightFlux(fevt::FdConstants::eCherMieScattered);
    TabulatedFunctionErrors& scatRCher = eyeRecData.GetLightFlux(fevt::FdConstants::eCherRayleighScattered);
    TabulatedFunctionErrors& longProf = eyeRecData.GetFRecShower().GetLongitudinalProfile();
    TabulatedFunctionErrors& dedxProf = eyeRecData.GetFRecShower().GetEnergyDeposit();
    ublas::vector<double> dEdX = ublas::prod(CFinv, n);
    ublas::vector<double> vChDir = ublas::prod(*(ChFl.GetDirectCherenkovMatrix()), dEdX);
    ublas::vector<double> vChMScat = ublas::prod(*(ChFl.GetMieScatteredCherenkovMatrix()), dEdX);
    ublas::vector<double> vChRScat = ublas::prod(*(ChFl.GetRayleighScatteredCherenkovMatrix()), dEdX);
    for (unsigned int i = 0; i < nF-fOutSideFOV; ++i) {
      const int j = i + fOutSideFOV;
      const double eDep = dEdX(j);
      const double eDepErr = sqrt(covdEdX[j]); // sqrt(covdEdX(j,j));
      const double dEdXperElectron = utl::EnergyDeposit(fDepth[j], fXmax, fEcut);
      const double Ne = eDep / dEdXperElectron;
      const double NeErr = eDepErr / dEdXperElectron;
      const double time = lightflux.GetX(i);
      dedxProf.PushBack(fDepth[j], 0, eDep, eDepErr);
      longProf.PushBack(fDepth[j], 0, Ne, NeErr);
      dirCher.PushBack(time, 0, vChDir[j], 0);
      scatMCher.PushBack(time, 0, vChMScat[j], 0);
      scatRCher.PushBack(time, 0, vChRScat[j], 0);
    }
  } else {
    fevt::EyeRecData& eyeRecData = eye.GetRecData();
    TabulatedFunctionErrors& dedxProf = eyeRecData.GetFRecShower().GetEnergyDeposit();
    ublas::vector<double> dEdX = ublas::prod(CFinv, n);
    for (unsigned int i = 0; i < nF-fOutSideFOV; ++i) {
      const int j = i + fOutSideFOV;
      const double eDep = dEdX(j);
      const double eDepErr = sqrt(covdEdX[j]); // sqrt(covdEdX(j,j));
      dedxProf.PushBack(fDepth[j], 0, eDep, eDepErr);
    }
  }

  return true;
}


/******************************************************************\
 *
 *  calculation of relative detector efficiency as
 *  function of wavelength
 *
\******************************************************************/
void FdProfileReconstructor::CalculateDetEff(const fevt::Eye& eye)
{
  // Relative detector efficiency of first telescope
  fRelEff.Clear();

  //  fevt::Eye::ConstTelescopeIterator itTel = eye.TelescopesBegin();
  // const fdet::Telescope& dettel =
  //  det::Detector::GetInstance().GetFDetector().GetTelescope(*itTel);
  //const utl::TabulatedFunction& fEfficiencyCalibration = dettel.GetModelEfficiency();

  // Note: Since the measured efficiency is really a telescope property (albeit
  //       currently the same for all telescopes), this shouldn't simply grab the
  //       efficiency of the first telescope. --Steffen 18.5.09
  fRelEff = det::Detector::GetInstance().GetFDetector().GetEye(eye).GetTelescope(1).GetMeasuredRelativeEfficiency();
}


/******************************************************************\
 *
 *  calculates vector of shower positions in local
 *  eye coordinates and fills slant depth
 *
\******************************************************************/
bool
FdProfileReconstructor::CalculateGeometryAndDepth(fevt::Eye& eye)
{
  // clear vectors
  fTrackPositions.clear();
  fDepth.clear();

  det::Detector& theDetector = det::Detector::GetInstance();

  const Atmosphere& atmo = theDetector.GetAtmosphere();
  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();

  // set the Cherenkov model energy threshold for this calculation
  atmo.SetCherenkovEnergyCutoff(fEcut);

  fevt::EyeRecData& eyeRecData = eye.GetRecData();

  const Point& eyePos = theDetector.GetFDetector().GetEye(eye).GetPosition();
  const Vector coreEyeVec = fCore - eyePos;
  const double coreEyeDistance = coreEyeVec.GetMag();
  CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(fCore);

  const double cosTheta = fAxis.GetCosTheta(coreCS);
  const bool flatEarth = (cosTheta > fInclinedModelMaxCosZenith);

  if (flatEarth && cosTheta < 0.1) {
    cerr << " FdProfileReconstructor::CalculateGeometryAndDepth() - "
        " skipping upward/horizontal track\n"
        " please use inclined model for these tracks" << endl;
    return false;
  }

  double minDistance;
  double maxDistance;
  if (!InitializeAtmosphere(flatEarth, minDistance, maxDistance))
    return false;

  const ProfileResult& verticalDepthProfile = atmo.EvaluateDepthVsHeight();
  const double minHeight = verticalDepthProfile.MinX();
  const double maxHeight = verticalDepthProfile.MaxX();

  const ProfileResult& depthProfile = (flatEarth?
                                       verticalDepthProfile :
                                       atmo.EvaluateSlantDepthVsDistance());

  // new profiles will be filled if some points are outside atmosphere
  unsigned int nSkip = 0;
  TabulatedFunctionErrors newLightFlux;

  TabulatedFunctionErrors& lightflux = eyeRecData.GetLightFlux(fevt::FdConstants::eTotal);
  const unsigned int nAperturePoints = lightflux.GetNPoints();

  vector<Point> trackBin(2);

  // real data inside field of view
  for (unsigned int i = 0; i < nAperturePoints; ++i) {

    // time at start and end of bin
    double t[2] = { lightflux.GetX(i)/ns - lightflux.GetXErr(i)/ns,
                    lightflux.GetX(i)/ns + lightflux.GetXErr(i)/ns };
    double Xdepth[2] = { 0, 0 };

    bool skipThisPoint = false;
    for (int j = 0; j < 2; ++j) {

      const double chi_i = fChi0 - 2* atan(kSpeedOfLight / fRp *(t[j] - fT0));
      double l = coreEyeDistance * sin(chi_i) / sin(fChi0 - chi_i);
      double distance = -l;

      // vector and point detector-shower at t=t[j]
      const Vector chi_i_vec  = coreEyeVec + l * fAxis;
      const Point chi_i_point = fCore      + l * fAxis;
      const double height = wgs84.PointToLatitudeLongitudeHeight(chi_i_point).get<2>();

      trackBin[j] = chi_i_point;

      if (distance<minDistance || distance>maxDistance ||
          height < minHeight   || height > maxHeight )
        skipThisPoint = true;
      else
        Xdepth[j] = (flatEarth ?
                depthProfile.Y(height)/cosTheta :
                depthProfile.Y(distance));
    }

    if (!skipThisPoint) {
      fDepth.push_back(0.5*(Xdepth[0] + Xdepth[1]));
      fTrackPositions.push_back(trackBin);
      newLightFlux.PushBack(lightflux.GetX(i), lightflux.GetXErr(i),
                            lightflux.GetY(i), lightflux.GetYErr(i));
    } else
      ++nSkip;
  }

  if (nSkip == nAperturePoints)
    return false; // zero points left ...

  // refill light at diaphragm if there were points outside atmosphere
  if (nSkip > 0) {
    if (fVerbosity > -1)
      cout << "      CalculateGeometryAndDepth: " << nSkip
           << " points outside atmosphere! -> "
           << " refilling light at aperture" << endl;

    TabulatedFunctionErrors& flightflux =
      eye.GetRecData().GetLightFlux(fevt::FdConstants::eFluorTotal);

    lightflux.Clear();
    if (!fErrorCalcMode)
      flightflux.Clear();
    for (unsigned int i = 0; i < newLightFlux.GetNPoints(); ++i) {
      const double x  = newLightFlux.GetX(i);
      const double ex = newLightFlux.GetXErr(i);
      const double y  = newLightFlux.GetY(i);
      const double ey = newLightFlux.GetYErr(i);
      lightflux.PushBack(x, ex, y, ey);
      if (!fErrorCalcMode)
        flightflux.PushBack(x, ex, y, ey);
    }
  }

  AddPointsOutsideFOV();

  return true;
}


/******************************************************************\
 *
 *  adds points outside field of view (for Cherenkov beam)
 *
\******************************************************************/
void FdProfileReconstructor::AddPointsOutsideFOV() {

  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();
  CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(fCore);
  const double coreHeight    = wgs84.PointToLatitudeLongitudeHeight(fCore).get<2>();
  const double cosTheta      = fAxis.GetCosTheta(coreCS);

  const bool flatEarth = (cosTheta>fInclinedModelMaxCosZenith);

  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  const ProfileResult& depthProfile = (flatEarth ?
                                        atmo.EvaluateDepthVsHeight() :
                                        atmo.EvaluateSlantDepthVsDistance());
  const ProfileResult& distProfile  = (flatEarth ?
                                        atmo.EvaluateHeightVsDepth() :
                                        atmo.EvaluateDistanceVsSlantDepth());

  double minDepth = distProfile.MinX();
  double maxDepth = distProfile.MaxX();

  const double dX    = 25.*(g/cm2); // depth step size
  const double xInit = 50.*(g/cm2); // start of shower

  double lastDepth = 0.;
  if ( flatEarth ) {
    const double lastHeight =
      wgs84.PointToLatitudeLongitudeHeight(fTrackPositions[0][0]).get<2>();
    lastDepth = depthProfile.Y(lastHeight)/cosTheta;
  }
  else {
    const double lastDistance = (fTrackPositions[0][0]-fCore).GetMag()
      * ((fCore-fTrackPositions[0][0])*fAxis>0 ? +1.0 : -1.0);
    lastDepth = depthProfile.Y(lastDistance);
  }

  vector<Point> trackBin(2);
  unsigned int nAperturePoints = fTrackPositions.size();
  fOutSideFOV=0;

  while ( fDepth[0] > xInit &&
          fOutSideFOV < nAperturePoints ) {  // do not add more points than seen

    trackBin[1]=fTrackPositions[0][0];  // end of this bin is start of previous bin
    const double Xdepth=lastDepth-dX;

    if (!flatEarth) {

      if (Xdepth> minDepth && Xdepth<maxDepth) {
        const double distance = distProfile.Y(Xdepth);
        const Point pDistance(fCore - distance*fAxis); // by convention
        trackBin[0] = pDistance;
        fDepth.push_front(Xdepth+dX/2.);
        fTrackPositions.push_front(trackBin);
        lastDepth = Xdepth;
        fOutSideFOV++;
      }
      else
        break;

    }
    else { // flat earth

      const double verticalDepth = Xdepth*cosTheta;
      if (verticalDepth>minDepth &&
        verticalDepth<maxDepth ) {
        const double h = distProfile.Y(verticalDepth) - coreHeight;
        const double l = h/cosTheta;
        trackBin[0] = fCore + l*fAxis;
        fDepth.push_front(Xdepth+dX/2.);
        fTrackPositions.push_front(trackBin);
        lastDepth = Xdepth;
        fOutSideFOV++;
      }
      else
        break;
    }
  }

  if(fVerbosity>-1)
    cout << "        added " << fOutSideFOV
        << " points outside FoV" << endl;

}

/******************************************************************\
 *
 *  create profile structures for current eye
 *
\******************************************************************/
void FdProfileReconstructor::InitProfiles(fevt::Eye& eye)
{
  evt::ShowerFRecData& showerFRecData=eye.GetRecData().GetFRecShower();

  // set the cut off value in data structure
  showerFRecData.SetEnergyCutoff(fEcut);

  // create longitudinal profile if necessary
  if (!showerFRecData.HasLongitudinalProfile())
    showerFRecData.MakeLongitudinalProfile();
  showerFRecData.GetLongitudinalProfile().Clear();

  // create energy deposit profile if necessary
  if (!showerFRecData.HasEnergyDeposit())
    showerFRecData.MakeEnergyDeposit();
  showerFRecData.GetEnergyDeposit().Clear();

  // create fluorescence photons at track if necessary
  if (!showerFRecData.HasFluorescencePhotons())
    showerFRecData.MakeFluorescencePhotons();
  showerFRecData.GetFluorescencePhotons().Clear();

  // direct cherenkov light at aperture
  if(!eye.GetRecData().HasLightFlux(fevt::FdConstants::eCherDirect) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eCherDirect);
  eye.GetRecData().GetLightFlux(fevt::FdConstants::eCherDirect).Clear();

  // Rayleigh scattered cherenkov light at aperture
  if(!eye.GetRecData().HasLightFlux(fevt::FdConstants::eCherRayleighScattered) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eCherRayleighScattered);
  eye.GetRecData().GetLightFlux(fevt::FdConstants::eCherRayleighScattered).Clear();

  // Mie scattered cherenkov light at aperture
  if(!eye.GetRecData().HasLightFlux(fevt::FdConstants::eCherMieScattered) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eCherMieScattered);
  eye.GetRecData().GetLightFlux(fevt::FdConstants::eCherMieScattered).Clear();

  // direct and multiple scattered fluorescence light at aperture
  if(!eye.GetRecData().HasLightFlux(fevt::FdConstants::eFluorMultScattered) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eFluorMultScattered);
  eye.GetRecData().GetLightFlux(fevt::FdConstants::eFluorMultScattered).Clear();

  if(!eye.GetRecData().HasLightFlux(fevt::FdConstants::eFluorDirect) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eFluorDirect);
  eye.GetRecData().GetLightFlux(fevt::FdConstants::eFluorDirect).Clear();

  // multiple scattered cherenkov light at aperture
  if(!eye.GetRecData().HasLightFlux(fevt::FdConstants::eCherMultScattered) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eCherMultScattered);
  eye.GetRecData().GetLightFlux(fevt::FdConstants::eCherMultScattered).Clear();

}

/******************************************************************\
 *
 * get current guess of shower maximum
 *
\******************************************************************/
double FdProfileReconstructor::GuessShowerMaximum(evt::Event& event,
                                                  const fevt::Eye& eye)
{
  double xMax=700.*(g/cm/cm);

  // first try Gaisser Hillas parameters of event ...
  if(!fAlwaysGuessXMaxFromAperture
      && event.HasRecShower()&& event.GetRecShower().HasFRecShower()
      && event.GetRecShower().GetFRecShower().HasGHParameters()) {
    xMax=event.GetRecShower().GetFRecShower().GetGHParameters().GetXMax();
  }
  // ... then Gaisser Hillas parameters of eye ...
  else if( !fAlwaysGuessXMaxFromAperture
           && eye.GetRecData().HasFRecShower()
           && eye.GetRecData().GetFRecShower().HasGHParameters()) {
    xMax=eye.GetRecData().GetFRecShower().GetGHParameters().GetXMax();
  }

  //  ... otherwise use depth where light at aperture is maximum
  else {
    int iMax=-1;
    double lightMax=0.;
    const TabulatedFunctionErrors& lightflux = eye.GetRecData().GetLightFlux(fevt::FdConstants::eTotal);
    for (unsigned int i=0; i<lightflux.GetNPoints() ; i++) {
      double flux=lightflux.GetY(i)/lightflux.GetXErr(i);
      if(flux>lightMax) {
        lightMax=flux;
        iMax=i;
      }
    }
    if(iMax>=0)
      xMax=fDepth[iMax+fOutSideFOV];

    // limit this first guess to reasonable values
    const double kxMaxMin=  400.*(g/cm/cm);
    const double kxMaxMax= 1100.*(g/cm/cm);
    if (xMax<kxMaxMin)
      xMax=kxMaxMin;
    else if (xMax>kxMaxMax)
      xMax=kxMaxMax;
  }
  return xMax;
}

/******************************************************************\
 *
 * calculation of total errors on energy, Xmax etc.
 *
\******************************************************************/
void FdProfileReconstructor::CalculateTotalErrors(fevt::Eye& eye)
{

  ShowerFRecData& recShower=eye.GetRecData().GetFRecShower();

  // propagate atmopsheric and geometric uncertainties if required
  if(fPropagateAtmUncertainties || fPropagateGeometryErrors)
    PropagateUncertainties(eye);

  // calculate statistical invisible energy uncertainty

  const double calEnergy = recShower.GetEmEnergy();
  const double calEnergy2 = pow(calEnergy, 2);

  double scaleFac=1.;
  if(fRescaleErrors) {
    double chisq = recShower.GetGHParameters().GetChiSquare();
    double ndof  = recShower.GetGHParameters().GetNdof();
    scaleFac = sqrt(chisq/ndof);
  }
  const double calEnergyErr  =
    recShower.GetEmEnergyError(ShowerFRecData::eStatistical) * scaleFac;
  const double calEnergyStatVar = calEnergyErr * calEnergyErr;

  double calEnergyGeomVar=0.;
  if (fPropagateGeometryErrors) {
    if (fGeomVariance.size()==eNFitParameters)
      calEnergyGeomVar = fGeomVariance[eEem];
    else
      calEnergyGeomVar = calEnergy2; // refitting failed --> set error to 100%
  }


  double calEnergyAtmVar=0.;
  if (fPropagateAtmUncertainties) {
    if (fAtmVariance.size()==eNFitParameters)
      calEnergyAtmVar = fAtmVariance[eEem];
    else
      calEnergyAtmVar = calEnergy2; // refitting failed --> set error to 100%
  }

  // ICRC 2013
  /*
  const double additionalAtmVar =
    calEnergy2 * (pow(gHorizontalUniformityError, 2) +
                  pow(gAtmVariability, 2) +
                  pow(gNightlyRelativeCalibration, 2)
                  );
  */
   // ICRC 2019
  const double dEr = gHorizontalUniformityError_par[0];
  const double logE = log10(calEnergy);
  const double logEs = gHorizontalUniformityError_par[1];
  const double fp0 = gHorizontalUniformityError_par[2];
  const double fp1 = gHorizontalUniformityError_par[3];
  const double fp2 = gHorizontalUniformityError_par[4];
  const double fE = max(fp0+fp1*(logE -18.)+fp2*(logE -18.)*(logE -18.),0.);
  const double fEs =    fp0+fp1*(logEs-18.)+fp2*(logEs-18.)*(logEs-18.);
  const double gHorizontalUniformityError_E = dEr*fE/fEs;

  const double additionalAtmVar =
    calEnergy2 * (pow(gHorizontalUniformityError_E, 2) +
                  pow(gAtmVariability, 2) +
                  pow(gNightlyRelativeCalibration, 2) +
                  pow(gEnergyTimeDrif, 2) +
                  pow(gEnergyEyeTel, 2)
                  );

  calEnergyAtmVar += additionalAtmVar;

  // ICRC 2019
  const double calorEnergyConexMissingVar =
    calEnergy2 * pow(gEnergyResolutionConexMissingFactor,2);


  // invibisible energy factors

//added by Matias Tueros
  utl::Vector fAxis=recShower.GetAxis();

  utl::Point fCore = recShower.GetCorePosition();

  CoordinateSystemPtr coreCS=fwk::LocalCoordinateSystem::Create(fCore);

  double cosZenithAngle = fAxis.GetCosTheta(coreCS);

//end addition

  const double invEnergyFac =
    utl::InvisibleEnergy::Factor(calEnergy,
                                 fInvisibleEnergyModel,
                                 fComposition, cosZenithAngle);
  const double invEnergyFacDeriv =
    utl::InvisibleEnergy::FactorDerivative(calEnergy,
                                           fInvisibleEnergyModel,
                                           fComposition, cosZenithAngle);
  // total energy
  double totEnergy = invEnergyFac * calEnergy;

  // ICRC 2013
  // shower to shower fluctuations of invisible energy
  //const double invEVar = utl::InvisibleEnergy::FactorVariance(totEnergy);

  // ICRC 2019
  const double invEVar  = utl::InvisibleEnergy::FactorVariance(calEnergy,totEnergy);

  // geometry and statistical error of CalEnergy
  const double deriv   = invEnergyFacDeriv * calEnergy + invEnergyFac;
  const double deriv2  = deriv * deriv;
  const double statVar = deriv2 * calEnergyStatVar;
  const double geomVar = deriv2 * calEnergyGeomVar;
  const double atmVar  = deriv2 * calEnergyAtmVar;
  // ICRC 2013
  //const double totEnergyStatErr = sqrt(invEVar + geomVar + statVar);
  // ICRC 2019
  const double ResolutionConexMissingFactorVar  = deriv2 * calorEnergyConexMissingVar;
  const double totEnergyStatErr = sqrt(invEVar + geomVar + statVar + ResolutionConexMissingFactorVar);

  // store total errors
  // ICRC 2013
  //recShower.SetEmEnergy(calEnergy, sqrt(calEnergyGeomVar+calEnergyStatVar),
  //                      ShowerFRecData::eStatistical);
  // ICRC 2019
  recShower.SetEmEnergy(calEnergy, sqrt(calEnergyGeomVar+calEnergyStatVar+calorEnergyConexMissingVar),
                        ShowerFRecData::eStatistical);
  recShower.SetEmEnergyError(sqrt(calEnergyAtmVar), ShowerFRecData::eAtmospheric);
  recShower.SetTotalEnergy(totEnergy, totEnergyStatErr,ShowerFRecData::eStatistical);
  recShower.SetTotalEnergyError(sqrt(atmVar),ShowerFRecData::eAtmospheric);

  GaisserHillas4Parameter* const gh4 =
    dynamic_cast<GaisserHillas4Parameter*>(&recShower.GetGHParameters());

  // this should never, never happen ...
  if (!gh4) {
    ERROR("No 4par GH structure????");
    throw MissingEventDataException();
  }

  const double dEdXmax     = gh4->GetNMax();
  const double dEdXmax_err = scaleFac * gh4->GetNMaxError();
  const double xmax        = gh4->GetXMax();
  const double xmax_err    = scaleFac * gh4->GetXMaxError();
  const double x0          = gh4->GetShapeParameter(gh::eX0);
  const double x0_err      = scaleFac * gh4->GetShapeParameterError(gh::eX0);
  const double lambda      = gh4->GetShapeParameter(gh::eLambda);
  const double lambda_err  = scaleFac * gh4->GetShapeParameterError(gh::eLambda);

  // store total GH parameter errors
  if (fPropagateGeometryErrors) {
    double x0_tot_err, xmax_tot_err, nmax_tot_err, lambda_tot_err;

    if (fGeomVariance.size() == eNFitParameters) {
      x0_tot_err    = sqrt(x0_err * x0_err + fGeomVariance[eX0]);
      xmax_tot_err  = sqrt(xmax_err * xmax_err + fGeomVariance[eXmax]);
      nmax_tot_err  = sqrt(dEdXmax_err * dEdXmax_err + fGeomVariance[edEdXmax]);
      lambda_tot_err= sqrt(lambda_err * lambda_err + fGeomVariance[eLambda]);
      recShower.SetXmaxError(sqrt(fAtmVariance[eXmax]), ShowerFRecData::eAtmospheric);
      recShower.SetXmaxError(xmax_err, ShowerFRecData::eStatistical);
    } else {
      x0_tot_err = x0;
      xmax_tot_err = xmax;
      nmax_tot_err = dEdXmax;
      lambda_tot_err = lambda;
      recShower.SetXmaxError(xmax, ShowerFRecData::eStatistical);
    }

    gh4->SetXMax(xmax, xmax_tot_err);
    gh4->SetNMax(dEdXmax, nmax_tot_err, true);
    gh4->SetShapeParameter(gh::eX0, x0, x0_tot_err);
    gh4->SetShapeParameter(gh::eLambda, lambda, lambda_tot_err);
  } else {
    gh4->SetXMax(xmax, xmax_err);
    gh4->SetNMax(dEdXmax, dEdXmax_err, true);
    gh4->SetShapeParameter(gh::eX0, x0, x0_err);
    gh4->SetShapeParameter(gh::eLambda, lambda ,lambda_err);
  }

  // do some printout ...
  if (fVerbosity > -1) {
    cout << "\n    E [EeV]    = " << totEnergy/EeV
         << " +/- " << sqrt(statVar)/EeV << " (stat.) "
         << " +/- " << sqrt(geomVar)/EeV << " (geom.) "
         << " +/- " << sqrt(invEVar)/EeV << " (inv.) "
         << " +/- " << sqrt(atmVar)/EeV << " (atm.) "
         << endl;
    cout << " Xmax [g/cm^2] = " << xmax/(g/cm/cm)
         << " +/- " << xmax_err/(g/cm/cm) << " (stat.) ";
    if (fGeomVariance.size() == eNFitParameters)
      cout << " +/- " << sqrt(fGeomVariance[eXmax])/(g/cm2) << " (geom.) "
           << " +/- " << sqrt(fAtmVariance[eXmax])/(g/cm2) << " (atm.) ";
    cout << endl;
  }
}


/******************************************************************\
 *
 * brute force propagation of axis fit, SDP and atmospheric
 * uncertainties
 *
 *  (well, actually not so brute if fmaxErrPropPoints is
 *   set properly ...)
 *
\******************************************************************/
bool FdProfileReconstructor::PropagateUncertainties(fevt::Eye& eye)
{

  // flag that we are now in error mode
  fErrorCalcMode=true;

  fevt::EyeRecData& eyeRecData  = eye.GetRecData();

  // use current Xmax for ChFl-Matrix
  fXmax = eyeRecData.GetFRecShower().GetGHParameters().GetXMax();

  // set verbosity for err prop.
  int v1Orig=fVerbosity, v2Orig=fEnergyFitter->GetVerbosity();
  fVerbosity=fErrPropVerb-1;
  fEnergyFitter->SetVerbosity(fErrPropVerb-1);

  // calculate undisturbed, but rebinned result

  std::vector<std::vector<double> > fitResults;

  const bool haveDefault = ReFitProfile(eye, 0., 0., 0., 0., 0., fitResults);

  if ( haveDefault && fPropagateGeometryErrors ) {

    const unsigned int kNGeomVars = 5;

    if(fErrPropVerb>-1)
      cout << "\n    propagating geometry uncertainties ";
    if(fErrPropVerb>0)
      cout << "...\n" << endl;

    while ( true ) {

      // a) Rp

      if(!ReFitProfile(eye, 1., 0., 0., 0., 0., fitResults)) break;
      if(!ReFitProfile(eye,-1., 0., 0., 0., 0., fitResults)) break;

      // b) Chi0

      if(!ReFitProfile(eye, 0., 1., 0., 0., 0., fitResults)) break;
      if(!ReFitProfile(eye, 0.,-1., 0., 0., 0., fitResults)) break;

      // c) T0

      if(!ReFitProfile(eye, 0., 0., 1., 0., 0., fitResults)) break;
      if(!ReFitProfile(eye, 0., 0.,-1., 0., 0., fitResults)) break;

      // d) SDP phi

      if(!ReFitProfile(eye, 0., 0., 0., 1., 0., fitResults)) break;
      if(!ReFitProfile(eye, 0., 0., 0.,-1., 0., fitResults)) break;

      // e) SDP theta

      if(!ReFitProfile(eye, 0., 0., 0., 0., 1., fitResults)) break;
      if(!ReFitProfile(eye, 0., 0., 0., 0.,-1., fitResults)) break;

      break;

    }


    const bool geomSuccess = ( fitResults.size() == (kNGeomVars*2 +1) );

    if(fErrPropVerb==0) cout << endl;


    if ( geomSuccess ) {

      // rescale errors if required
      //
      //  on average, the axis chi2/ndf is much larger than 1
      //  or in other word the centroid errors are presumably
      //  underestimated by a factor of sqrt(chi2/ndf)
      //  (see for instance introduction in PDG long writeup)
      //  This propagates directly to RpErr etc. and correspondingly
      //  to energy, Xmax errors
      //
      double scaleFac=1.;
      if(fRescaleErrors>0) {
        double chi2 = eyeRecData.GetTimeFitChiSquare();
        double ndf  = eyeRecData.GetTimeFitNDof();
        if(ndf>0)
        scaleFac=chi2/ndf;
      }

      // correlation coefficients
      const double rhoRT = eyeRecData.GetRpTZeroCorrelation();
      const double rhoRC = eyeRecData.GetRpChi0Correlation();
      const double rhoCT = eyeRecData.GetChi0TZeroCorrelation();
      const double rhoPT = eyeRecData.GetSDPCorrThetaPhi();

      for ( unsigned int i=0; i<eNFitParameters; i++ ) {
        double x0,x1,x2;

        x0=fitResults[0][i];

        double err[kNGeomVars];  // Rp,Chi0,T0,Phi and Theta errors

        for ( unsigned int j=0;j<kNGeomVars;j++) {
        x1=fitResults[2*j+1][i]; x2=fitResults[2*j+2][i];

        double e1=(x1-x0);
        double e2=(x0-x2);
        err[j]=0.5*(e1+e2);
        }

        // time fit variance
        fGeomVariance.push_back(err[0]*err[0] + err[1]*err[1] + err[2]*err[2]
                                + 2. * err[0] * err[2] * rhoRT
                                + 2. * err[0] * err[1] * rhoRC
                                + 2. * err[1] * err[2] * rhoCT);

        fGeomVariance[i]*=scaleFac;

        // add SDP fit variance
        fGeomVariance[i]+= err[3]*err[3] + err[4]*err[4]
        + 2. * err[3] * err[4] * rhoPT;

      }
    }
  }

  if ( haveDefault && fPropagateAtmUncertainties ) {

    if(fErrPropVerb>-1) cout << "    propagating atmospheric uncertainties ";
    if(fErrPropVerb>0)  cout << "...\n" << endl;

    const unsigned int kNAtmVars = 1;

    fitResults.resize(1);

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

    // change VAODs
    atmo.SetUncertaintyBound(Atmosphere::eMie,1.);
    if ( ReFitProfile(eye, 0., 0., 0., 0., 0., fitResults) ) {
      atmo.SetUncertaintyBound(Atmosphere::eMie,-1.);
      ReFitProfile(eye, 0., 0., 0., 0., 0., fitResults);
    }


    // reset original atmosphere
    atmo.SetUncertaintyBound(Atmosphere::eMie,0.);

    const bool atmSuccess = ( fitResults.size() == (kNAtmVars*2 +1) );

    if ( atmSuccess ) {

      for ( unsigned int i=0; i<eNFitParameters; i++ ) {
        double x0,x1,x2;
        x0=fitResults[0][i];

        double v=0.;
        for ( unsigned int j=0;j<kNAtmVars;j++) {
        x1=fitResults[2*j+1][i]; x2=fitResults[2*j+2][i];
        const double e1=(x1-x0);
        const double e2=(x0-x2);
        const double err=0.5*(e1+e2);
        v+=(err*err);
        }

        fAtmVariance.push_back(v);

      }
    }
  }

  if(fErrPropVerb==0) cout << endl;

  // restore verbosities
  fVerbosity=v1Orig;
  fEnergyFitter->SetVerbosity(v2Orig);

  return true;

}

/*********************************************************************\
 *
 *  steering routine for profile refit
 *
\*********************************************************************/
bool FdProfileReconstructor::ReFitProfile(fevt::Eye& eye,
                                          double nSigRp,
                                          double nSigChi0,
                                          double nSigT0,
                                          double nSigPhi,
                                          double nSigTheta,
                                          std::vector<std::vector<double> >& fitResults)
{
  DeleteCFMatrix();

  if(fErrPropVerb == 0 && !fitResults.empty())
    cout << "." << std::flush;

  // save a copy of the original lightflux and dEdX as
  // they might be combined or modified

  fevt::EyeRecData& eyeRecData = eye.GetRecData();
  TabulatedFunctionErrors& lightFlux =
    eyeRecData.GetLightFlux(fevt::FdConstants::eTotal);
  ShowerFRecData& recShower=eye.GetRecData().GetFRecShower();
  TabulatedFunctionErrors& dedxProf=recShower.GetEnergyDeposit();

  std::vector<double> x,y,ex,ey;
  std::vector<double> dEx,dEy,dEex,dEey;
  for (unsigned int i = 0; i < lightFlux.GetNPoints(); i++) {
    x.push_back(lightFlux.GetX(i));
    y.push_back(lightFlux.GetY(i));
    ex.push_back(lightFlux.GetXErr(i));
    ey.push_back(lightFlux.GetYErr(i));
    dEx.push_back(dedxProf.GetX(i));
    dEy.push_back(dedxProf.GetY(i));
    dEex.push_back(dedxProf.GetXErr(i));
    dEey.push_back(dedxProf.GetYErr(i));
  }

  // combine data points to speed up re-fits
  if( (int) x.size() > fmaxErrPropPoints ) {

    lightFlux.Clear();

    const double nCombi=(double) x.size()/ (double) fmaxErrPropPoints;
    double combined=0.;
    double t1=x[0]-ex[0];
    double sum=0.;
    double sum2=0.;

    for (unsigned int i=0;i<x.size();i++) {
      combined+=1.;
      sum += y[i];
      sum2+= ey[i]*ey[i];

      if(t1==0.) t1=x[i]-ex[i];

      // check for time discontinuity (mirror crossing)
      bool timeGap=false;
      if ( i<x.size()-1 && x[i]+ex[i] != x[i+1]-ex[i+1] )
        timeGap = true;

      if(combined>=nCombi||i==x.size()-1||timeGap) {
        double dt    = x[i]-t1+ex[i];
        lightFlux.PushBack    (t1+dt/2.,dt/2.,sum,sqrt(sum2));
        sum=0.; sum2=0.; t1=0.; combined=0.;
      }
    }
  }

  // start refit
  SetGeometry(eye,nSigRp,nSigChi0, nSigT0, nSigPhi, nSigTheta);
  bool success=false;
  dedxProf.Clear();

  if(CalculateGeometryAndDepth(eye)) {
    if(CalculateProfiles(eye)) {
      if(fEnergyFitter->ReFitProfile(eye, fChFlPtr)) {

        std::vector<double> result(eNFitParameters);
        result[eX0]=fEnergyFitter->GetX0();
        result[eXmax]=fEnergyFitter->GetXmax();
        result[edEdXmax]=fEnergyFitter->GetdEdXmax();
        result[eLambda]=fEnergyFitter->GetLambda();
        result[eEem]=fEnergyFitter->GetEmEnergy();
        fitResults.push_back(result);

        success=true;

      }
    }
  }

  // refill original profiles
  lightFlux.Clear();
  dedxProf.Clear();
  for ( unsigned int i=0;i<x.size(); i++ ) {
    lightFlux.PushBack(x[i],ex[i],y[i],ey[i]);
    dedxProf.PushBack(dEx[i],dEex[i],dEy[i],dEey[i]);
  }
  return success;

}

/******************************************************************\
 *
 * sets geometry data members altering Rp, Chi0 and/or T0
 * by nSig of their statistical errors
 *
\******************************************************************/
void FdProfileReconstructor::SetGeometry(fevt::Eye& eye,
                                         double nSigRp,
                                         double nSigChi0,
                                         double nSigT0,
                                        double nSigPhi,
                                        double nSigTheta)
 {

  fevt::EyeRecData& eyeRecData = eye.GetRecData();
  fT0   = eyeRecData.GetTZero()+nSigT0*eyeRecData.GetTZeroError();
  fChi0 = eyeRecData.GetChiZero()+nSigChi0*eyeRecData.GetChiZeroError();
  fRp   = eyeRecData.GetRp()+nSigRp*eyeRecData.GetRpError();

  const fdet::Eye& detEye  = det::Detector::GetInstance().GetFDetector().GetEye(eye);
  CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
  const Point& eyepos       = detEye.GetPosition();



  const double sdpTheta = eyeRecData.GetSDP().GetTheta (eyeCS)+
                          nSigTheta*eyeRecData.GetSDPThetaError();
  const double sdpPhi   = eyeRecData.GetSDP().GetPhi   (eyeCS)+
                          nSigPhi*eyeRecData.GetSDPPhiError();
  Vector sdp = Vector(1,sdpTheta,sdpPhi,eyeCS, Vector::kSpherical);

  Vector vertical(0,0,1,eyeCS);
  Vector horizontalInSDP = cross(sdp,vertical); // horizontal
  horizontalInSDP.Normalize();

  //  if ( horizontalInSDP.GetPhi(eyeCS) < 0 ||
  //     horizontalInSDP.GetPhi(eyeCS) > kPi) {
  //  if(fVerbosity>-1)
  //    cout << "        swapping horizontal ... " << endl;
  //  horizontalInSDP *= -1;
  //}

  Transformation rot (Transformation::Rotation(-fChi0,sdp, eyeCS));
  Vector axis (rot* horizontalInSDP); // apply rotation
  axis.Normalize();
  fAxis = axis;

  const double core_distance = fRp / sin( kPi - fChi0);
  const Vector coreEyeVec    = core_distance * horizontalInSDP;
  fCore                      = eyepos + coreEyeVec;

}


inline void FdProfileReconstructor::DeleteCFMatrix() {

  std::map<int, CherenkovFluorescenceMatrix *>::iterator iter;
  for (iter=fChFlPtrMap.begin(); iter!=fChFlPtrMap.end(); ++iter)
    delete iter->second;
  fChFlPtrMap.clear();
  fChFlPtr = NULL;

}

void FdProfileReconstructor::SetYieldRefit(bool what) {

  fYieldRefit=what;

  if ( fEnergyFitter != NULL )
    fEnergyFitter->SetYieldRefit(what);
  else
    std::cerr << " FdProfileReconstructor::SetYieldRefit() "
        << " - EnergyFitter can not be set! " << std::endl;


}


/******************************************************************\
 *
 * Sets the energy fitter verbosity during run-time.
 *
\******************************************************************/

void FdProfileReconstructor::SetEnergyFitterVerbosity(int verbosity) {
  if (fEnergyFitter != NULL)
    fEnergyFitter->SetVerbosity(verbosity);
}


/******************************************************************\
 *
 * Gets the energy fitter verbosity during run-time.
 *
\******************************************************************/

int FdProfileReconstructor::GetEnergyFitterVerbosity() const {
  if (fEnergyFitter != NULL)
    return fEnergyFitter->GetVerbosity();
  else {
    std::cerr << " FdProfileReconstructor::GetEnergyFitterVerbosity() "
        << " - EnergyFitter not initialized! " << std::endl;
    return 0;
  }
}

/******************************************************************\
 *
 * Initializes slant atmosphere (flatEarth=false)
 *
\******************************************************************/

bool
FdProfileReconstructor::InitializeAtmosphere(bool flatEarth,
                                             double& minDistance,
                                             double& maxDistance)
{
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  if (flatEarth) {

    const atm::ProfileResult& depthProfile  = atmo.EvaluateDepthVsHeight();

    CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(fCore);
    const double cosTheta = fAxis.GetCosTheta(coreCS);

    const double maxHeight = depthProfile.MaxX();
    const double minHeight = depthProfile.MinX();
    const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();
    const double coreHeight =
      wgs84.PointToLatitudeLongitudeHeight(fCore).get<2>();
    minDistance = (coreHeight-maxHeight)/cosTheta + 1*mm;
    maxDistance = (coreHeight-minHeight)/cosTheta - 1*mm;

  } else {

    try {
      if (fErrorCalcMode) {
        // less precision for error propagation, see GAP-2007-007, Fig. 5
        atmo.InitSlantProfileModel(fCore, fAxis, 100.*g/cm2);
      } else {
        if (!fYieldRefit || fIsMultiEyeEvent)
          atmo.InitSlantProfileModel(fCore, fAxis, fInclinedModelDepthBinning);
      }
    } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
      if (fVerbosity > -1) {
        cout << "      CalculateGeometryAndDepth - "
                "error initializing InclinedAtmosphericProfile " << endl;
      }
      return false;
    }

    const atm::ProfileResult slantDepthProfile = atmo.EvaluateSlantDepthVsDistance();
    minDistance = slantDepthProfile.MinX() + 1*mm;
    maxDistance = slantDepthProfile.MaxX() - 1*mm;

  }

  return true;
}