FdEnergyDepositFinder.cc

Go to the documentation of this file.

#include "FdEnergyDepositFinder.h"

#include "../FdProfileReconstructorKG/UncorrelatedUncertainties.h"

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

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

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

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeHeader.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeRecData.h>
#include <evt/ShowerFRecData.h>
#include <fevt/EyeRecData.h>

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

#include <atm/Atmosphere.h>

#include <boost/io/ios_state.hpp>

#include <iostream>
#include <sstream>
#include <iomanip>

using namespace evt;
using namespace fevt;
using namespace utl;
using namespace fwk;
using namespace std;
using namespace atm;


namespace FdEnergyDepositFinderKG {

  const string gHRule = "    --------------------------------------------------------------";


  VModule::ResultFlag
  FdEnergyDepositFinder::Init()
  {
    auto& cc = *CentralConfig::GetInstance();

    auto topB = cc.GetTopBranch("FdEnergyDepositFinder");

    // invisible energy options
    auto eCalcB = topB.GetChild("invisibleEnergy");
    const string compositionOption = eCalcB.GetChild("composition").Get<string>();

    fComposition =
      (compositionOption == "eData") ?
        InvisibleEnergy::eData :
        (compositionOption == "eMixed") ?
          InvisibleEnergy::eMixedComposition :
          (compositionOption == "eProton") ?
            InvisibleEnergy::eProton : InvisibleEnergy::eIron;

    const auto interactionOption = eCalcB.GetChild("interactionModel").Get<string>();
    if (interactionOption == "eSIBYLL21")
      fInvisibleEnergyModel = InvisibleEnergy::eSIBYLL21;
    else if (interactionOption == "eQGSJET01")
      fInvisibleEnergyModel = InvisibleEnergy::eQGSJET01;
    else if (interactionOption == "eDATA")
      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;
    }

    fGHShapePars = GHShapeParameters(topB);

    auto errPropB = topB.GetChild("errorPropagation");
    errPropB.GetChild("verbosity").GetData(fErrPropVerbosity);
    errPropB.GetChild("rescaleErrors").GetData(fRescaleErrors);
    errPropB.GetChild("maxErrPropPoints").GetData(fMaxErrPropPoints);

    auto evSelB = topB.GetChild("eventSelection");
    evSelB.GetChild("requireNonMono").GetData(fOnlyNonMono);
    evSelB.GetChild("maxZenith").GetData(fMaxZenith);
    topB.GetChild("denseMatrixDim").GetData(fDenseMatrixDim);
    topB.GetChild("antiAliasingFilterCorrection").GetData(fAafCorrection);
    topB.GetChild("useNoiseBins").GetData(fUseNoiseBins);

    ostringstream info;
    info << "Parameters:\n"
            "\tinvisible energy model: "
         << interactionOption << '/'
         << compositionOption << "\n"
            "\t" << (fOnlyNonMono ? "" : "not ") << "skipping mono events\n"
            "\tskipping upward showers with theta > " << fMaxZenith/deg;
    INFO(info);

    fCFMatrixCalculator.Init();
    fCFMatrixCalculatorDense.Init();
    fProfileFitter.Init();

    return eSuccess;
  }


  VModule::ResultFlag
  FdEnergyDepositFinder::Run(evt::Event& event)
  {
    if (!event.HasFEvent())
      return eSuccess;

    auto& fdEvent = event.GetFEvent();

    for (auto eyeIt = fdEvent.EyesBegin(ComponentSelector::eHasData),
         end = fdEvent.EyesEnd(ComponentSelector::eHasData);
         eyeIt != end; ++eyeIt) {

      auto& eye = *eyeIt;

      if (IsEventCandidate(eye)) {

        ostringstream info;
        info << " calculating profile for eye " << eye.GetId();
        INFO(info);

        fRecResultsPlus.clear();
        fRecResultsMinus.clear();

        for (unsigned int i = 0; i < eNErrorPropagation; ++i) {

          fErrPropStage = EErrorPropStage(i);
          bool fullCalculation = false;

          for (unsigned int j = 0; j < 2; ++j) {

            EErrorPropType errPropType = (j == 0) ? eMinus : ePlus;
            fullCalculation = false;

            // create a copy of fevt::Eye via Event
            Event eventCopy;
            if (fErrPropStage != eDefaultReco) {
              const auto id = eye.GetId();
              eventCopy.MakeFEvent();
              eventCopy.GetFEvent().MakeEye(id);
              eventCopy.GetFEvent().GetEye(id) = eye;
              if (!PrepareEyeCopy(eventCopy.GetFEvent().GetEye(id), errPropType))
                break;
              fCFMatrixCalculator.SetMethod(CFMatrixCalculator::eFast);
              fCFMatrixCalculator.SetVerbosity(-1);
              fCFMatrixCalculatorDense.SetMethod(CFMatrixCalculator::eFast);
              fCFMatrixCalculatorDense.SetVerbosity(-1);
              fProfileFitter.SetAafCorrection(false, false);
              fProfileFitter.SetVerbosity(-1);
            } else {
              fCFMatrixCalculator.SetMethod(CFMatrixCalculator::ePrecise);
              fCFMatrixCalculator.SetVerbosity(0);
              fCFMatrixCalculatorDense.SetMethod(CFMatrixCalculator::ePrecise);
              fCFMatrixCalculatorDense.SetVerbosity(0);
              fProfileFitter.SetAafCorrection(fAafCorrection, fUseNoiseBins);
              fProfileFitter.SetVerbosity(0);
            }

            auto& eyeEvt =
              (fErrPropStage == eDefaultReco) ? eye : eventCopy.GetFEvent().GetEye(eye.GetId());

            auto& shower = eyeEvt.GetRecData().GetFRecShower();
            shower.SetEnergyCutoff(fElectronEnergyThreshold);

            // first version of matrix
            fCFMatrixCalculator.BuildMatrix(eyeEvt);
            fCFMatrixCalculator.SetVerbosity(-1);
            fCFMatrixCalculatorDense.BuildMatrix(
              eyeEvt,
              false,
              (fErrPropStage == eDefaultReco &&
               fCFMatrixCalculator.GetFOVSize() > 0 &&
               fDenseMatrixDim/fCFMatrixCalculator.GetFOVSize() > 0) ?
              fDenseMatrixDim/fCFMatrixCalculator.GetFOVSize() : 1
            );
            fCFMatrixCalculatorDense.SetVerbosity(-1);

            if (fCFMatrixCalculator.HasMatrix() && fCFMatrixCalculatorDense.HasMatrix()) {

              // first version of profile
              if (!fProfileCalculator.CalculateProfile(eyeEvt, fCFMatrixCalculator))
                break;

              if (FitProfile(eyeEvt, eGH2dEdXChi2) &&    // prefit
                  FitProfile(eyeEvt, eGH4dEdXChi2) &&    // 4-GH fit
                  FitProfile(eyeEvt, eGH4LightLogLike))  // final fit
                fullCalculation = true;

            }

            if (fErrPropStage == eDefaultReco && !fullCalculation) {
              // reset intermediate results of not successful
              shower.SetTotalEnergy(0, 1e20*eV);
              shower.SetEmEnergy(0, 1e20*eV);
              if (shower.HasGHParameters()) {
                shower.GetGHParameters().SetXMax(0, 1000*g/cm2);
                shower.GetGHParameters().SetNMax(0, 0);
                shower.SetXmaxError(1000*g/cm2, ShowerFRecData::eStatistical);
              }
            } else if (fErrPropStage != eDefaultReco) {
              const GaisserHillas4Parameter& gh4Pars =
                shower.GetGHParameters<GaisserHillas4Parameter>();
              auto& recoResult =
                (errPropType == ePlus) ? fRecResultsPlus : fRecResultsMinus;
              recoResult[fErrPropStage] = gh4Pars;
            }

            // reset mie uncertainty (may have changed in PrepareEyeCopy())
            const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
            atmo.SetUncertaintyBound(Atmosphere::eMie, 0);

            if (!fullCalculation || fErrPropStage == eDefaultReco)
              break;

          } // plus/minus loop

          if (!fullCalculation)
            break;

        } // error prop loop

        cout << gHRule << endl;

        PropagateUncertainties(eye);

      }
    }

    return eSuccess;
  }


  bool
  FdEnergyDepositFinder::FitProfile(fevt::Eye& eye,
                                    const EFitFunctionType type)
  {
    auto& shower = eye.GetRecData().GetFRecShower();
    shower.SetEnergyCutoff(fElectronEnergyThreshold);
    fProfileFitter.ResetStartParameters();

    if (!shower.HasGHParameters()) {
      shower.MakeGHParameters(GaisserHillas4Parameter(fGHShapePars.Type()));
      if (!GuessGHParameters(shower))
        return false;
      DumpCurrentParameters(shower, false, eUndefined);
    }

    fProfileFitter.SetFitMode(ProfileFitter::eNMax);
    {
      GaisserHillas4Parameter& gh4Pars =
        shower.GetGHParameters<GaisserHillas4Parameter>();
      fProfileFitter.SetStartParameters(gh4Pars);
    }
    fProfileFitter.SetLightFluxData(fCFMatrixCalculator, fCFMatrixCalculatorDense);

    fProfileFitter.SetProfileData(shower.GetEnergyDeposit());
    fProfileFitter.SetFunctionType(type);

    if (fProfileFitter.Fit()) {

      GaisserHillas4Parameter& gh4Pars = shower.GetGHParameters<GaisserHillas4Parameter>();
      fProfileFitter.FillGHParameters(gh4Pars);

      double calorEnergy = gh4Pars.GetIntegral();
      double calorEnergyErr = calorEnergy * gh4Pars.GetIntegralError();

      // refit at final step using energy as parameter instead of Xmax
      if (type == eFinalFitStep) {
        fProfileFitter.SetFitMode(ProfileFitter::eIntegral);
        fProfileFitter.SetStartParameters(gh4Pars);
        if (fProfileFitter.Fit())
          fProfileFitter.GetEnergy(calorEnergy, calorEnergyErr);
        else
          return false;
      }
      shower.SetEmEnergy(calorEnergy, calorEnergyErr, ShowerFRecData::eStatistical);
      shower.SetEmEnergyError(0, ShowerFRecData::eAtmospheric);

      // rebuilt matrix if Xmax changed too much and not in final round
      const double currXmax = gh4Pars.GetXMax();
      const double maxDeltaXmax = 10*g/cm2;
      const bool updateCFM =
        (fabs(currXmax-fCFMatrixCalculatorDense.GetXmax()) > maxDeltaXmax &&
         type != eFinalFitStep);

      if (updateCFM) {
        fCFMatrixCalculator.BuildMatrix(eye);
        fCFMatrixCalculatorDense.BuildMatrix(
          eye,
          false,
          (fErrPropStage == eDefaultReco &&
           fCFMatrixCalculator.GetFOVSize() > 0 &&
           fDenseMatrixDim/fCFMatrixCalculator.GetFOVSize() > 0) ?
          fDenseMatrixDim/fCFMatrixCalculator.GetFOVSize() : 1
        );
      }

      if (fCFMatrixCalculator.HasMatrix() && fCFMatrixCalculatorDense.HasMatrix()) {

        if (!fProfileCalculator.CalculateProfile(eye, fCFMatrixCalculator))
          return false;
        DumpCurrentParameters(shower, updateCFM, type);

        return true;
      }

    }

    return false;
  }


  void
  FdEnergyDepositFinder::PropagateUncertainties(fevt::Eye& eye)
    const
  {
    auto& eyeRecData = eye.GetRecData();
    auto& recShower = eyeRecData.GetFRecShower();

    if (!recShower.HasGHParameters())
      return;

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

    const double dEdXmax = gh4.GetNMax();
    const double xmax = gh4.GetXMax();
    const double shapePar1 = gh4.GetShapeParameter(fGHShapePars.Type(eOne));
    const double shapePar2 = gh4.GetShapeParameter(fGHShapePars.Type(eTwo));
    const double calorEnergy = recShower.GetEmEnergy();

    // added by Matias Tueros
    const auto& axis = recShower.GetAxis();
    const auto& core = recShower.GetCorePosition();
    const auto coreCS = fwk::LocalCoordinateSystem::Create(core);

    const double cosZenithAngle = axis.GetCosTheta(coreCS);
    // end addition

    const double invEnergyFac =
      utl::InvisibleEnergy::Factor(calorEnergy,
                                   fInvisibleEnergyModel,
                                   fComposition, cosZenithAngle);
    const double totalEnergy = invEnergyFac * calorEnergy;

    // if one or more propagation steps failed, set all errors to 100%
    if (fRecResultsPlus.size() != eNErrorPropagation - 1 ||
        fRecResultsMinus.size() != eNErrorPropagation - 1) {

      if (fErrPropVerbosity > -1) {
        int failStage = 0;
        for (int j = eFirstErrProp; j < eNErrorPropagation; ++j)  {
          const EErrorPropStage tmp = EErrorPropStage(j);
          if (fRecResultsPlus.find(tmp) == fRecResultsPlus.end() ||
              fRecResultsMinus.find(tmp) == fRecResultsMinus.end()) {
            failStage = j;
            break;
          }
        }
        cout << " ==[ErrProp]==> numerical error propagation failed at stage "
             << failStage<< "\n"
                "                setting uncertainties to 100%!" << endl;
      }

      gh4.SetXMax(xmax , xmax);
      gh4.SetNMax(dEdXmax, dEdXmax, true);
      gh4.SetShapeParameter(fGHShapePars.Type(eOne), shapePar1, fabs(shapePar1));
      gh4.SetShapeParameter(fGHShapePars.Type(eTwo), shapePar2, fabs(shapePar2));

      recShower.SetXmaxError(xmax, ShowerFRecData::eStatistical);

      recShower.SetEmEnergyError(calorEnergy, ShowerFRecData::eStatistical);
      recShower.SetEmEnergyError(0, ShowerFRecData::eAtmospheric);
      recShower.SetTotalEnergy(totalEnergy, totalEnergy, ShowerFRecData::eStatistical);
      recShower.SetTotalEnergyError(0, ShowerFRecData::eAtmospheric);
      return;
    }

    // --- propagate uncertainties for all shower parameters
    enum EShowerParameter {
      eShapePar1 = 0,
      eXmax,
      edEdXmax,
      eShapePar2,
      eCalorEnergy,
      eLastFitParameter = eCalorEnergy
    };

    double geomVariance[eLastFitParameter+1];
    double atmVariance[eLastFitParameter+1];

    for (int i = 0; i <= eLastFitParameter; ++i) {

      EShowerParameter currentParameter = EShowerParameter(i);

      double uncertainty[eNErrorPropagation];
      for (int j = eFirstErrProp; j < eNErrorPropagation; ++j) {
        const evt::GaisserHillas4Parameter& ghPlus =
          fRecResultsPlus.find(EErrorPropStage(j))->second;
        const evt::GaisserHillas4Parameter& ghMinus =
          fRecResultsMinus.find(EErrorPropStage(j))->second;
        switch (currentParameter) {
        case eShapePar1:
          uncertainty[j] = 0.5*(ghPlus.GetShapeParameter(fGHShapePars.Type(eOne)) -
                                ghMinus.GetShapeParameter(fGHShapePars.Type(eOne)));
          break;
        case eXmax:
          uncertainty[j] = 0.5*(ghPlus.GetXMax()-ghMinus.GetXMax());
          break;
        case edEdXmax:
          uncertainty[j] = 0.5*(ghPlus.GetNMax()-ghMinus.GetNMax());
          break;
        case eShapePar2:
          uncertainty[j] = 0.5*(ghPlus.GetShapeParameter(fGHShapePars.Type(eTwo)) -
                                ghMinus.GetShapeParameter(fGHShapePars.Type(eTwo)));
          break;
        case eCalorEnergy:
          uncertainty[j] = 0.5*(ghPlus.GetIntegral()-ghMinus.GetIntegral());
          break;
        }
      }

      const double rhoRT = eyeRecData.GetRpTZeroCorrelation();
      const double rhoRC = eyeRecData.GetRpChi0Correlation();
      const double rhoCT = eyeRecData.GetChi0TZeroCorrelation();

      const double timeFitVariance =
        Sqr(uncertainty[eRp]) +
        Sqr(uncertainty[eT0]) +
        Sqr(uncertainty[eChi0]) +
        2 * uncertainty[eRp] * uncertainty[eT0] * rhoRT +
        2 * uncertainty[eRp] * uncertainty[eChi0] * rhoRC +
        2 * uncertainty[eChi0] * uncertainty[eT0] * rhoCT;

      const double rhoPT = eyeRecData.GetSDPCorrThetaPhi();

      const double sdpFitVariance =
        Sqr(uncertainty[eSDPTheta]) +
        Sqr(uncertainty[eSDPPhi]) +
        2 * uncertainty[eSDPPhi] * uncertainty[eSDPTheta] * rhoPT;

      geomVariance[i] = timeFitVariance + sdpFitVariance;
      atmVariance[i] = Sqr(uncertainty[eMie]);

    }

    // ICRC 2013, Energy scale Table: 2, Uncorrelated uncertainties
    // that go into the FD energy.
    //const double atmVar1 = pow(gHorizontalUniformityError * calorEnergy, 2);
    //const double atmVar2 = pow(gAtmVariability * calorEnergy, 2);
    //const double atmVar3 = pow(gNightlyRelativeCalibration * calorEnergy, 2);

    // ICRC 2019 Uncorrelated uncertainties
    // that go into the FD energy.
    const double dEr = gHorizontalUniformityError_par[0];
    const double logE = log10(calorEnergy);
    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 le = logE - 18;
    const double fE = max(fp0 + le*(fp1 + fp2*le), 0.);
    const double les = logEs - 18;
    const double fEs = fp0 + les*(fp1 + fp2*les);
    const double atmVar1 = Sqr(dEr*fE/fEs * calorEnergy);

    const double atmVar2 = Sqr(gAtmVariability             * calorEnergy);
    const double atmVar3 = Sqr(gNightlyRelativeCalibration * calorEnergy);
    const double atmVar4 = Sqr(gEnergyTimeDrif             * calorEnergy);
    const double atmVar5 = Sqr(gEnergyEyeTel               * calorEnergy);

    // optionally rescale stat. uncert. of GH fit using chi2
    const double scaleFac = fRescaleErrors ?
      sqrt(recShower.GetGHParameters().GetChiSquare() /
           recShower.GetGHParameters().GetNdof()) :
      1;
    const double calorEnergyStatVar =
      Sqr(recShower.GetEmEnergyError(ShowerFRecData::eStatistical)*scaleFac);
    const double calorEnergyGeomVar = geomVariance[eCalorEnergy];

    // ICRC 2019
    const double calorEnergyConexMissingVar = Sqr(gEnergyResolutionConexMissingFactor * calorEnergy);

    // ICRC 2013
    //const double calorEnergyAtmVar =
    //atmVariance[eCalorEnergy]  + atmVar1 + atmVar2 + atmVar3;
    // ICRC 2019
    const double calorEnergyAtmVar =
      atmVariance[eCalorEnergy] + atmVar1 + atmVar2 + atmVar3 + atmVar4 + atmVar5;

    // ICRC 2013
    //recShower.SetEmEnergy(calorEnergy,sqrt(calorEnergyGeomVar+calorEnergyStatVar),
    //                      ShowerFRecData::eStatistical);
    // ICRC 2019
    recShower.SetEmEnergy(
      calorEnergy,
      sqrt(calorEnergyGeomVar + calorEnergyStatVar + calorEnergyConexMissingVar),
      ShowerFRecData::eStatistical
    );

    recShower.SetEmEnergyError(sqrt(calorEnergyAtmVar), ShowerFRecData::eAtmospheric);

    // total energy uncertainty including shower-to-shower fluctuations of invisible energy
    const double invEnergyFacDeriv =
      utl::InvisibleEnergy::FactorDerivative(calorEnergy,
                                             fInvisibleEnergyModel,
                                             fComposition, cosZenithAngle);
    // ICRC 2013
    //const double invEVar = utl::InvisibleEnergy::FactorVariance(totalEnergy);
    // ICRC 2019
    const double invEVar = utl::InvisibleEnergy::FactorVariance(calorEnergy, totalEnergy);

    const double deriv = invEnergyFacDeriv * calorEnergy + invEnergyFac;
    // invEnergyFac = (Ecal+Einv)/Ecal = 1 + Einv/Ecal
    // invEnergyFacDeriv = d(Einv/Ecal)/d(Ecal) =   d(Einv)/d(Ecal) / Ecal - Einv/Ecal/Ecal
    // deriv                                    = ( d(Einv)/d(Ecal) / Ecal - Einv/Ecal/Ecal ) * Ecal + 1 + Einv/Ecal =  d(Einv)/d(Ecal) + 1
    // Etot = Ecal + Einv
    // dEtot/dEcal = 1 + d(Einv)/d(Ecal)
    // dEtot = ( 1 + d(Einv)/d(Ecal) ) * dEcal
    // deriv =   1 + d(Einv)/d(Ecal)
    const double deriv2 = Sqr(deriv);
    const double statVar = deriv2 * calorEnergyStatVar;
    const double geomVar = deriv2 * calorEnergyGeomVar;
    // ICRC 2019
    const double ConexMissingVar = deriv2 * calorEnergyConexMissingVar;

    const double totalEnergyAtmErr = sqrt(deriv2 * calorEnergyAtmVar);

    // ICRC 2013
    //const double totalEnergyStatErr = sqrt(invEVar + geomVar+statVar);
    // ICRC 2019
    const double totalEnergyStatErr = sqrt(invEVar + geomVar + statVar + ConexMissingVar);

    recShower.SetTotalEnergy(totalEnergy, totalEnergyStatErr, ShowerFRecData::eStatistical);
    recShower.SetTotalEnergyError(totalEnergyAtmErr, ShowerFRecData::eAtmospheric);

    // finally set total GH uncertainties, note: atm. uncert. not included here
    const double nmaxErr = scaleFac * gh4.GetNMaxError();
    const double xmaxErr = scaleFac * gh4.GetXMaxError();
    const double shapePar1Err = scaleFac * gh4.GetShapeParameterError(fGHShapePars.Type(eOne));
    const double shapePar2Err = scaleFac * gh4.GetShapeParameterError(fGHShapePars.Type(eTwo));
    const double shapePar1TotErr = sqrt(Sqr(shapePar1Err) + geomVariance[eShapePar1]);
    const double xmaxTotErr = sqrt(Sqr(xmaxErr) + geomVariance[eXmax]);
    const double nmaxTotErr = sqrt(Sqr(nmaxErr) + geomVariance[edEdXmax]);
    const double shapePar2TotErr = sqrt(Sqr(shapePar2Err) + geomVariance[eShapePar2]);

    gh4.SetShapeParameter(fGHShapePars.Type(eOne), shapePar1, shapePar1TotErr);
    gh4.SetShapeParameter(fGHShapePars.Type(eTwo), shapePar2, shapePar2TotErr);
    gh4.SetXMax(xmax, xmaxTotErr);
    gh4.SetNMax(dEdXmax, nmaxTotErr, true);

    recShower.SetXmaxError(sqrt(atmVariance[eXmax]), ShowerFRecData::eAtmospheric);
    recShower.SetXmaxError(xmaxTotErr, ShowerFRecData::eStatistical);

    // do some printout ...
    if (fErrPropVerbosity >- 1) {
      boost::io::ios_all_saver save(cout);
      cout.flags(std::ios::scientific);
      cout << "\n E/EeV = " << totalEnergy/EeV;
      cout.precision(1);
      cout << " +/- " << sqrt(statVar)/EeV << " (stat.)"
              " +/- " << sqrt(geomVar)/EeV << " (geom.)"
              " +/- " << sqrt(invEVar)/EeV << " (inv.)"
              " +/- " << totalEnergyAtmErr/EeV << " (atm.)\n"
              " Xmax/(g/cm^2) = ";
      cout.precision(2);
      cout << xmax/(g/cm2);
      cout.precision(1);
      cout << " +/- " << xmaxErr/(g/cm2) << " (stat.)"
              " +/- " << sqrt(geomVariance[eXmax])/(g/cm2) << " (geom.)"
              " +/- " << sqrt(atmVariance[eXmax])/(g/cm2) << " (atm.)"
           << endl;
    }
  }


  bool
  FdEnergyDepositFinder::PrepareEyeCopy(fevt::Eye& eye,
                                        EErrorPropType errPropType)
    const
  {
    auto& eyeRecData = eye.GetRecData();
    for (auto telIt = eye.TelescopesBegin(fevt::ComponentSelector::eInDAQ),
         end = eye.TelescopesEnd(fevt::ComponentSelector::eInDAQ);
         telIt != end; ++telIt) {

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

      auto& telRecData = telIt->GetRecData();
      if (!telRecData.HasLightFlux(fevt::FdConstants::eTotal) ||
          !telRecData.HasLightFlux(fevt::FdConstants::eBackground))
        continue;

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

      // combine data points to speed up re-fits
      vector<double> x;
      vector<double> y;
      vector<double> ex;
      vector<double> ey;
      vector<double> vBkg;
      unsigned int nInRange = 0;
      for (unsigned int i = 0; i < photonTrace.GetNPoints(); ++i) {
        x.push_back(photonTrace.GetX(i));
        y.push_back(photonTrace.GetY(i));
        ex.push_back(photonTrace.GetXErr(i));
        ey.push_back(photonTrace.GetYErr(i));
        vBkg.push_back(Sqr(bgTrace.GetYErr(i)));
        if (telRecData.IsSpotFarFromBorder(photonTrace.GetX(i)))
          ++nInRange;
      }

      if (int(nInRange) > fMaxErrPropPoints) {

        photonTrace.Clear();
        bgTrace.Clear();

        const double nCombi = max(1, int(double(nInRange) / fMaxErrPropPoints + 0.5));

        int combined = 0;
        double t1 = x[0] - ex[0];
        double sum = 0;
        double sum2 = 0;
        double bgVarianceSum = 0;

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

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

          if (combined >= nCombi || i == x.size()-1) {
            const double dt = x[i] + ex[i] - t1;
            photonTrace.PushBack(t1 + 0.5*dt, 0.5*dt, sum, sqrt(sum2));
            bgTrace.PushBack(t1 + 0.5*dt, 0.5*dt, 0, sqrt(bgVarianceSum));
            sum = 0;
            sum2 = 0;
            t1 = 0;
            combined = 0;
            bgVarianceSum = 0;
          }
        }
      }

    }

    // clear energy deposit
    auto& recShower = eye.GetRecData().GetFRecShower();
    auto& dedxProf=recShower.GetEnergyDeposit();
    dedxProf.Clear();

    // change geometry
    if (fErrPropStage >= eFirstGeometryProp &&
        fErrPropStage <= eLastGeometryProp) {

      const auto& detEye = det::Detector::GetInstance().GetFDetector().GetEye(eye);
      const auto eyeCS = detEye.GetEyeCoordinateSystem();
      const auto& eyepos = detEye.GetPosition();

      const double timeFitFactor = fRescaleErrors ?
        (eyeRecData.GetTimeFitNDof() ?
           sqrt(eyeRecData.GetTimeFitChiSquare()/eyeRecData.GetTimeFitNDof()) : 1) :
        1;
      const double sdpFitFactor = fRescaleErrors ?
        (eyeRecData.GetSDPFitNDof() ?
           sqrt(eyeRecData.GetSDPFitChiSquare()/eyeRecData.GetSDPFitNDof()) : 1) :
        1;

      double t0 = eyeRecData.GetTZero();
      double chi0 = eyeRecData.GetChiZero();
      double rp = eyeRecData.GetRp();
      double sdpTheta = eyeRecData.GetSDP().GetTheta(eyeCS);
      double sdpPhi = eyeRecData.GetSDP().GetPhi(eyeCS);
      if (fErrPropStage == eRp)
        rp += eyeRecData.GetRpError()*errPropType*timeFitFactor;
      else if (fErrPropStage == eChi0)
        chi0 += eyeRecData.GetChiZeroError()*errPropType*timeFitFactor;
      else if (fErrPropStage == eT0)
        t0 += eyeRecData.GetTZeroError()*errPropType*timeFitFactor;
      else if (fErrPropStage == eSDPTheta)
        sdpTheta += eyeRecData.GetSDPThetaError()*errPropType*sdpFitFactor;
      else if (fErrPropStage == eSDPPhi)
        sdpPhi += eyeRecData.GetSDPPhiError()*errPropType*sdpFitFactor;

      const Vector sdp(1, sdpTheta, sdpPhi, eyeCS, Vector::kSpherical);
      const Vector vertical(0,0,1, eyeCS);
      const auto horizontalInSDP = Normalized(Cross(sdp, vertical));  // horizontal

      const auto rot = Transformation::Rotation(-chi0, sdp, eyeCS);
      const auto axis = Normalized(rot*horizontalInSDP);
      const double coreDistance = rp / sin(chi0);
      const auto coreEyeVec = coreDistance * horizontalInSDP;
      const auto core = eyepos + coreEyeVec;
      recShower.SetAxis(axis);
      recShower.SetCorePosition(core);
      try {
        recShower.SetCoreTime(eye.GetHeader().GetTimeStamp() -
                              TimeInterval(1000*100*ns),
                              rp, chi0, t0);
      } catch (OutOfBoundException& e) {
        ERROR(e.GetMessage());
        return false;
      }
    } else if (fErrPropStage == eMie) {
      const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
      atmo.SetUncertaintyBound(Atmosphere::eMie,errPropType);
    }

    return true;
  }


  bool
  FdEnergyDepositFinder::GuessGHParameters(ShowerFRecData& shower)
    const
  {
    if (!shower.HasGHParameters())
      shower.MakeGHParameters(GaisserHillas4Parameter());

    GaisserHillas4Parameter& gh4Pars = shower.GetGHParameters<GaisserHillas4Parameter>();

    // running average

    const auto& energyDeposit = shower.GetEnergyDeposit();
    const unsigned int nPoints = energyDeposit.GetNPoints();
    const double depthBin = 50*g/cm2;

    bool firstTime = true;
    double xMax = 0;
    double dEdXmax = 0;
    for (unsigned int i = 0; i < nPoints; ++i) {
      const double firstDepth = energyDeposit.GetX(i);
      unsigned int nAverage = 0;
      double wSum = 0;
      double ywSum = 0;
      double xSum = 0;
      for (unsigned int j = i; j < nPoints; ++j) {
        const double thisDepth = energyDeposit.GetX(j);
        ++nAverage;
        xSum += thisDepth;
        const double w = 1 / Sqr(energyDeposit.GetYErr(j));
        wSum += w;
        ywSum += energyDeposit.GetY(j)*w;

        if (thisDepth - firstDepth > depthBin) {
          const double yAverage = ywSum / wSum;
          if (firstTime || yAverage > dEdXmax) {
            dEdXmax = yAverage;
            xMax = xSum / nAverage;
            firstTime = false;
          }
          break;
        }
      }
    }

    gh4Pars.SetShapeParameter(
      fGHShapePars.Type(eOne),
      fGHShapePars.Mean(eOne, 1e18*eV),
      fGHShapePars.Sigma(eOne, 1e18*eV)
    );
    gh4Pars.SetShapeParameter(
      fGHShapePars.Type(eTwo),
      fGHShapePars.Mean(eTwo, 1e18*eV),
      fGHShapePars.Sigma(eTwo, 1e18*eV)
    );
    if (firstTime) {
      WARNING("can not estimate start parameters");
      gh4Pars.SetXMax(700*g/cm2, 0);
      gh4Pars.SetNMax(energyDeposit.GetY(0), 0, true);
      return false;
    } else {
      gh4Pars.SetXMax(xMax, 0);
      gh4Pars.SetNMax(dEdXmax, 0, true);
    }
    return true;
  }


  void
  FdEnergyDepositFinder::DumpCurrentParameters(const ShowerFRecData& shower,
                                               const bool cfmUpdate,
                                               const EFitFunctionType funcType)
    const
  {
    if (fErrPropVerbosity < 1 && fErrPropStage != eDefaultReco)
      return;
    else if (fErrPropStage > 0 && funcType == eGH2dEdXChi2)
      cout << " ==[ErrProp]==> stage " << fErrPropStage << endl;

    switch (funcType) {
    case eUndefined:
      cout << "             "
           << setw(9) << "Xmax"
           << setw(13) << "dEdXmax"
           << setw(12) << fGHShapePars.Name(eOne)
           << setw(12) << fGHShapePars.Name(eTwo)
           << "  chi2/ndf\n"
              "             "
           << setw(9) << "[g/cm^2]"
           << setw(13) << "[PeV cm^2/g]"
           << setw(12-2-fGHShapePars.UnitName(eOne).length()) << " " << "[" << fGHShapePars.UnitName(eOne) << "]"
           << setw(12-2-fGHShapePars.UnitName(eTwo).length()) << " " << "[" << fGHShapePars.UnitName(eTwo) << "]\n"
           << gHRule << "\n"
              " ==> init   :";
      break;
    case eGH2dEdXChi2:
      cout << " ==> prefit1:";
      break;
    case eGH4dEdXChi2:
      cout << " ==> prefit2:";
      break;
    case eGH4LightLogLike:
      cout << " ==> logLike:";
      break;
    default:
      return;
    }

    if (!shower.HasGHParameters()) {
      cout << " --- " << endl;
      return;
    }

    const GaisserHillas4Parameter& gh4Pars =
      shower.GetGHParameters<GaisserHillas4Parameter>();

    boost::io::ios_all_saver save(cout);
    cout.flags(std::ios::scientific);
    cout.precision(2);

    cout << setw(9) << fixed << gh4Pars.GetXMax()/g*cm2
         << setw(13) << scientific << gh4Pars.GetNMax()/(PeV/(g/cm2))
         << setw(12)
         << gh4Pars.GetShapeParameter(fGHShapePars.Type(eOne))/fGHShapePars.Unit(eOne)
         << setw(12)
         << gh4Pars.GetShapeParameter(fGHShapePars.Type(eTwo))/fGHShapePars.Unit(eTwo)
         << setw(5) << int(gh4Pars.GetChiSquare()) << "/"
         << std::left << gh4Pars.GetNdof() << std::right
         << (cfmUpdate ? "  *" : "") << endl;
  }


  bool
  FdEnergyDepositFinder::IsEventCandidate(const fevt::Eye& eye)
    const
  {
    ostringstream info;

    if (!eye.HasRecData()) {
      info << " skipping eye " << eye.GetId() << " without RecData";
      INFO(info);
      return false;
    }

    if (!eye.GetRecData().HasFRecShower()) {
      info << " skipping eye " << eye.GetId() << " without FRecShower";
      INFO(info);
      return false;
    }

    if (eye.GetRecData().GetTimeFitNDof() <= 0) {
      info << " skipping eye " << eye.GetId() << " with time-fit Ndf="
           << eye.GetRecData().GetTimeFitNDof();
      INFO(info);
      return false;
    }

    const auto& shower = eye.GetRecData().GetFRecShower();
    if (fOnlyNonMono) {
      const bool isMCgeo =
        eye.GetRecData().GetTimeFitNDof() == eye.GetRecData().GetTimeFitChiSquare() &&
        eye.GetRecData().GetRpError()/eye.GetRecData().GetRp() < 1e-4;
      if (shower.GetStationIds().empty() &&          // not hybrid
          eye.GetRecData().GetPCGF().empty() &&      // not PCGF
          !isMCgeo) {                                // not MC axis
        info << " skipping eye " << eye.GetId()
             << " with mono geometry";
        INFO(info);
        return false;
      }
    }

    if (fMaxZenith < 180*degree) {
      const auto& core = shower.GetCorePosition();
      const auto& axis = shower.GetAxis();
      const auto coreLocalCS = fwk::LocalCoordinateSystem::Create(core);
      if (axis.GetTheta(coreLocalCS) > fMaxZenith) {
        info << " skipping upward event in eye " << eye.GetId()
             << " theta=" << axis.GetTheta(coreLocalCS)/degree << " deg.";
        INFO(info);
        return false;
      }
    }

    return true;
  }

}