SdReconstruction/LDFFinderKG/LDFFinder.cc

Go to the documentation of this file.

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSimData.h>

#include <sevt/SEvent.h>
#include <sevt/StationRecData.h>
#include <sevt/EventTrigger.h>
#include <sevt/Station.h>
#include <sevt/StationConstants.h>

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

#include <det/Detector.h>

#include <sdet/SDetector.h>
#include <sdet/SDetectorConstants.h>
#include <sdet/Station.h>
#include <sdet/STimeVariance.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>

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

#include <utl/CoordinateSystem.h>
#include <utl/PhysicalConstants.h>
#include <utl/ErrorLogger.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/TimeStamp.h>
#include <utl/AxialVector.h>
#include <utl/Math.h>
#include <utl/Reader.h>
#include <utl/Accumulator.h>
#include <utl/NumericalErrorPropagator.h>
#include <utl/String.h>
#include <utl/Is.h>

#include <Minuit2/MnMinimize.h>
#include <Minuit2/MnHesse.h>
#include <Minuit2/FCNBase.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnUserParameters.h>
#include <Minuit2/MnUserCovariance.h>
#include <Minuit2/MnPrint.h>

#include "LDFFinder.h"
#include "ROptFitter.h"

#include <boost/range/iterator_range.hpp>
#include <cmath>

using namespace LDFFinderKG;
using namespace sevt;
using namespace fevt;
using namespace evt;
using namespace fwk;
using namespace std;
using namespace utl;

using CLHEP::HepRotation;


#define OUT(_x_) if ((_x_) <= fInfoLevel) cerr


namespace LDFFinderKG {

  // local helpers

  template<typename T>
  inline
  bool
  HasNaN(const T& result)
  {
    const auto npar = result.fPar.size();
    for (size_t i = 0; i < npar; ++i) {
      if (std::isnan(result.fPar[i]))
        return true;
      for (size_t j = i; j < npar; ++j)
        if (std::isnan(result.fCov(i, j)))
          return true;
    }
    return false;
  }


  inline
  string
  GetShowerSizeLabel(const LDFFitConfig& config)
  {
    ostringstream label;
    const int refdist = round(config.fLDF.GetReferenceDistance());
    label << "S" << refdist/meter;
    if (refdist < 1000*meter)
      label << ' ';
    return label.str();
  }


  inline
  const char*
  GetShapeLabel(const int iShape)
  {
    static const char* const labels[] = { "beta", "gamma", "mu", "tau" };
    return Is(iShape).InRange(0, 3) ? labels[iShape] : "UNKNOWN COMPONENT!!!";
  }


  class CoreShiftPropagator : public NumericalErrorPropagator {
  public:
    CoreShiftPropagator(const double coreX,
                        const double coreY,
                        const vector<double>& axisParameters) :
      fCore(coreX, coreY, 0, gBaryCS),
      fAxisParameters(axisParameters)
    { }

    void
    Transform(vector<double>& output, const vector<double>& input)
      const
    {
      const double& u = fAxisParameters[0];
      const double& v = fAxisParameters[1];
      const double& rCore = fAxisParameters[2];
      const double& ctCore = fAxisParameters[3];
      const double& newCoreX = input[0];
      const double& newCoreY = input[1];

      const double w = sqrt(1 - u*u - v*v);

      const Point origin = fCore + rCore*Vector(u, v, w, gBaryCS);
      const Point newCore = Point(newCoreX, newCoreY, 0, gBaryCS);
      const Vector newAxis = origin - newCore;
      const double newRCore = newAxis.GetMag();

      if (output.size() != 4)
        output = vector<double>(4);

      output[LDFFinder::eU] = newAxis.GetX(gBaryCS)/newRCore;
      output[LDFFinder::eV] = newAxis.GetY(gBaryCS)/newRCore;
      output[LDFFinder::eRCore] = newRCore;
      output[LDFFinder::eCTCore] = ctCore + newRCore - rCore;
    }

  protected:
    const Point fCore;
    const vector<double> fAxisParameters;
  };


  inline
  double
  GetTimeVariance(const LDFFitConfig& fconf, const sdet::STimeVariance& timeVariance,
                  const double signal, const StationRecData& sRec,
                  const double cosTheta, const double rho)
  {
    double var = timeVariance.GetTimeSigma2(signal, sRec.GetT50(), cosTheta, rho);
    if (fconf.fUseAdditionalGPSTimeVariance)
      var += sRec.GetGPSTimeVariance();
    return var;
  }

}


void
LDFFinder::OutputResults(const Event& event)
  const
{
  const ShowerSRecData& recShower = event.GetRecShower().GetSRecShower();

  OUT(eFinal)
    << "***** Final reconstructed shower data *****\n"
       "  LDF stage = " << recShower.GetLDFRecStage() << "\n"
       "  " << GetShowerSizeLabel(fLDFFitConfig) << " = "
            << recShower.GetShowerSize() << " +/- " << recShower.GetShowerSizeError() << " +/- "
            << recShower.GetShowerSizeSystematics() << " (sys.) VEM\n"
       "  core = ("
          << String::Make(recShower.GetCorePosition(), det::Detector::GetInstance().GetSiteCoordinateSystem(), meter, ", ") << ") m (site)\n"
       "       = (" << String::Make(recShower.GetCorePosition(), gBaryCS, meter, ", ") << ") +/- "
                "(" << String::Make(recShower.GetCoreError(), gBaryCS, meter, ", ") << ") m (bary)\n"
       "  x-y correlation = " << recShower.GetCorrelationXY() << "\n"
       "  beta = " << recShower.GetBeta() << " +/- " << recShower.GetBetaError() << "\n"
       "  gamma = " << recShower.GetGamma() << " +/- " << recShower.GetGammaError() << "\n"
       "  energy = " << recShower.GetEnergy()/EeV << " +/- " << recShower.GetEnergyError()/EeV << " EeV\n";
  const double curv = recShower.GetCurvature();
  if (curv) {
    const Vector& axis = recShower.GetAxis();
    const double rc = 1 / curv;
    OUT(eFinal)
      << "  axis = (" << String::Make(axis, gBaryCS, 1, ", ") << ") (bary)\n"
         "  theta = " << axis.GetTheta(gBaryCS)/degree << " +/- " << recShower.GetThetaError()/degree << " deg (bary)\n"
         "  phi = " << axis.GetPhi(gBaryCS)/degree << " +/- " << recShower.GetPhiError()/degree << " deg\n"
         "  theta-phi correlation = " << recShower.GetCorrelationThetaPhi() << "\n"
         "  Rc = " << rc/km << " +- " << recShower.GetCurvatureError()*rc*rc/km << " km\n"
         "  time chi2 = " << recShower.GetAngleChi2() << " / " << recShower.GetAngleNdof() << "\n"
         "  time residual = " << recShower.GetTimeResidualMean()/nanosecond << " +/- "
                              << recShower.GetTimeResidualSpread()/nanosecond << " ns\n";
  }
  const double rOpt = recShower.GetLDFOptimumDistance();
  if (rOpt)
    OUT(eFinal)
      << "  Ropt = " << rOpt/meter << " m\n";

  ++RunController::GetInstance().GetRunData().GetNamedCounters()[
    "LDFFinder/" + (boost::format("%|.5|") % recShower.GetLDFRecStage()).str()
  ];
}


VModule::ResultFlag
LDFFinder::Init()
{
  auto topB = CentralConfig::GetInstance()->GetTopBranch("LDFFinderKG");

  topB.GetChild("infoLevel").GetData(fInfoLevel);

  {
    const auto arrayType = topB.GetChild("arrayType").Get<string>();
    if (arrayType == "SD1500")
      fArrayType = sdet::Array::e1500;
    else if (arrayType == "SD750")
      fArrayType = sdet::Array::e750;
    else if (arrayType == "SD433")
      fArrayType = sdet::Array::e433;
    else {
      ERROR("Unknown array type");
      return eFailure;
    }
  }

  {
    const auto sigVariance = topB.GetChild("signalVariance").Get<string>();
    if (sigVariance == "eGAP2012_012")
      fLDFFitConfig.fSignalVarianceModel = wcd::eGAP2012_012;
    else if (sigVariance == "eGAP2014_035")
       fLDFFitConfig.fSignalVarianceModel = wcd::eGAP2014_035;
    else {
      ERROR("Unknown signal variance");
      return eFailure;
    }

    fLDFFitConfig.fUseAdditionalGPSTimeVariance =
      topB.GetChild("useAdditionalGPSTimeVariance").Get<bool>();

    // setup LDF
    const auto ldfType = topB.GetChild("ldfType").Get<string>();

    const double ldfReferenceDistance = topB.GetChild("ldfReferenceDistance").Get<double>();

    const auto ldfParameters = topB.GetChild("ldfParameters").Get<vector<double>>();

    const auto ldfUncertaintyParameters =
      topB.GetChild("ldfUncertaintyParameters").Get<vector<double>>();

    fLDFFitConfig.fLDF = LDF(ldfType, ldfReferenceDistance, ldfParameters, ldfUncertaintyParameters);
    fLDFFitConfig.fLDFShapeTreatment.resize(fLDFFitConfig.fLDF.GetNShapeParameters(), eModeled);
  }

  {
    // setup energy calibration
    auto energyBranch = topB.GetChild("energyCalibration");

    fEnergyConversion.fCicReferenceAngle = energyBranch.GetChild("referenceAngle").Get<double>();
    fEnergyConversion.fCicReferenceS38 = energyBranch.GetChild("referenceS38").Get<double>();
    energyBranch.GetChild("s38Range").GetData(fEnergyConversion.fCicS38Range);
    const auto& r = fEnergyConversion.fCicS38Range;
    if (r.first <= 0 || r.second <= 0 || r.first > r.second) {
      ERROR("Invalid <s38Range>!");
      return eFailure;
    }
    fEnergyConversion.SetCICParameters(
      energyBranch.GetChild("attenuationPar1").Get<vector<double>>(),
      energyBranch.GetChild("attenuationPar2").Get<vector<double>>(),
      energyBranch.GetChild("attenuationPar3").Get<vector<double>>()
    );

    fEnergyConversion.fEnergyConstant = energyBranch.GetChild("energyS38Const").Get<double>();
    fEnergyConversion.fEnergySlope = energyBranch.GetChild("energyS38Slope").Get<double>();
  }

  auto silentsBranch = topB.GetChild("silentStations");
  silentsBranch.GetChild("maxLDFValue").GetData(fSilentsVars.fMaxLDFValue);
  silentsBranch.GetChild("maxDistance").GetData(fSilentsVars.fMaxDistance);
  silentsBranch.GetChild("minNumberOfStations").GetData(fSilentsVars.fMinNumberOfStations);

  {
    const auto coreType = topB.GetChild("coreType").Get<string>();
    if (coreType == "MC")
      fCoreType = eMC;
    else if (coreType == "Hybrid")
      fCoreType = eHybrid;
    else if (coreType == "Fit")
      fCoreType = eFit;
    else {
      ERROR("Unknown core type");
      return eFailure;
    }
  }

  fLDFFitConfig.fFixCore = (fCoreType != eFit);

  topB.GetChild("fixedAxis").GetData(fFixedAxis);

  topB.GetChild("globalFit").GetData(fGlobalFit);

  topB.GetChild("maxTheta").GetData(fMaxTheta);

  const auto minimizationMethod = topB.GetChild("minimizationMethod").Get<string>();
  fUseMaxLikelihood = (minimizationMethod == "MaxLike");

  topB.GetChild("maxChi2").GetData(fMaxChi2);

  topB.GetChild("useSilentStations").GetData(fLDFFitConfig.fUseSilentStations);

  topB.GetChild("useSilentStationsFromOtherGrids").GetData(fUseSilentStationsFromOtherGrids);

  topB.GetChild("silentMaxRadius").GetData(fLDFFitConfig.fSilentMaxRadius);

  topB.GetChild("silentRadiusTransition").GetData(fLDFFitConfig.fSilentRadiusTransition);

  topB.GetChild("silentSignalThreshold").GetData(fLDFFitConfig.fSilentSignalThreshold);

  topB.GetChild("useSaturationRecovery").GetData(fLDFFitConfig.fUseSaturationRecovery);

  topB.GetChild("recoveryThreshold").GetData(fLDFFitConfig.fRecoveryThreshold);

  topB.GetChild("maxBaryToCoreDistance").GetData(fMaxBaryToCoreDistance);
  if (fMaxBaryToCoreDistance <= 0) {
    ERROR("<maxBaryToCoreDistance> must be > 0");
    return eFailure;
  }

  topB.GetChild("RoptN").GetData(fRoptN);
  if (fRoptN < 2)
    fRoptN = 0;

  auto stage0B = topB.GetChild("stage0");
  stage0B.GetChild("minNumberOfStations").GetData(fStage0.fMinNumberOfStations);
  stage0B.GetChild("useCorePosition").GetData(fStage0.fUseCorePosition);

  auto stage1B = topB.GetChild("stage1");
  stage1B.GetChild("minNumberOfStations").GetData(fStage1.fMinNumberOfStations);

  auto stage2B = topB.GetChild("stage2");
  stage2B.GetChild("minNumberRelaxBeta").GetData(fStage2.fMinNumberRelaxBeta);
  stage2B.GetChild("fixBetaRange").GetData(fStage2.fFixBetaRange);
  stage2B.GetChild("fixBetaLeverArmRange").GetData(fStage2.fFixBetaLeverArmRange);
  stage2B.GetChild("minNumberRelaxGamma").GetData(fStage2.fMinNumberRelaxGamma);
  stage2B.GetChild("fixGammaRange").GetData(fStage2.fFixGammaRange);
  stage2B.GetChild("fixGammaLeverArmRange").GetData(fStage2.fFixGammaLeverArmRange);

  auto stage3B = topB.GetChild("stage3");
  stage3B.GetChild("minNumberOfStations").GetData(fStage3.fMinNumberOfStations);
  stage3B.GetChild("minNumberForFullCurvatureFit").GetData(fStage3.fMinNumberForFullCurvatureFit);
  stage3B.GetChild("maxAxisAngleDifference").GetData(fStage3.fMaxAxisAngleDifference);
  stage3B.GetChild("outliersRejection").GetData(fStage3.fOutliersRejection);

  {
    ostringstream info;
    const auto& p = fEnergyConversion.fCicParameters;
    info << "\n"
            "  LDF type: " << fLDFFitConfig.fLDF.GetType() << "\n"
            "  Shower size estimator: " << GetShowerSizeLabel(fLDFFitConfig) << "\n"
            "  Minimization method: " << (fUseMaxLikelihood ? "Maximum-Likelihood" : "Chi2") << "\n"
            "  S38 conversion: " << "S38 = " << GetShowerSizeLabel(fLDFFitConfig) << " / CIC(x), "
            "  CIC(x) = ";
    for (unsigned int ix = 0, nx = Length(p)+1; ix < nx; ++ix) {
      if (!ix)
        info << 1;
      else {
        info << " + (";
        for (unsigned int iy = 0; iy < 3; ++iy) {
          if (!iy)
            info << p[ix-1][iy];
          else {
            info << " + " << p[ix-1][iy] << "*y";
            if (iy > 1)
              info << '^' << iy;
          }
        }
        info << ") * x";
        if (ix > 1)
          info << '^' << ix;
      }
    }
    info << " where "
            "x = -(sin^2(theta) - sin^2(38deg)), "
            "y = lg(s / S38_ref), "
            "s = min(max(S38, " << fEnergyConversion.fCicS38Range.first << "), "
                                << fEnergyConversion.fCicS38Range.second << "), "
            "S38_ref = " << fEnergyConversion.fCicReferenceS38 << " VEM\n"
            "  energy conversion: E = " << fEnergyConversion.fEnergyConstant << " eV * "
            "pow(S38 / VEM, " << fEnergyConversion.fEnergySlope << ")\n"
            "  Ropt fitting: " << (fRoptN < 2 ? "off" : "on");
    INFO(info);
  }

  return eSuccess;
}


VModule::ResultFlag
LDFFinder::Run(evt::Event& event)
{
  // nothing to do if there is no SEvent
  if (!(event.HasSEvent() && event.HasRecShower() &&
        event.GetRecShower().HasSRecShower()))
    return eContinueLoop;

  const auto& sEvent = event.GetSEvent();

  auto& showerRec = event.GetRecShower().GetSRecShower();

  int nStations = 0;
  {
    int nSaturated = 0;
    int nSilent = 0;
    for (const auto& station : sEvent.StationsRange())
      if (station.IsCandidate()) {
        ++nStations;
        if (station.IsLowGainSaturation())
          ++nSaturated;
      } else if (station.IsSilent())
        ++nSilent;

    OUT(eIntermediate)
      << "# candidate stations = " << nStations
      << " (" << nSaturated << " saturated), " << nSilent << " silent\n";
  }

  // not worth it
  if (nStations < fStage0.fMinNumberOfStations) {
    OUT(eFinal)
      << "Not enough stations.\n";
    return eContinueLoop;
  }

  // barycenter is done by the SdPlaneFitOG module
  fBarycenter = showerRec.GetBarycenter();
  fBaryTime = showerRec.GetBarytime();
  gBaryCS = showerRec.GetBarycenterCoordinateSystem();

  if (!fBaryTime)
    throw utl::NonExistentComponentException("needs barycenter!");

  LDFFitConfig config = fLDFFitConfig;
  LDFFitResult best(config);
  CurvFitResult cbest;

  // setup initial values
  {
    Point core;
    TimeStamp coreTime;
    Vector showerAxis;

    if (!FixCore(core, coreTime, showerAxis, event))
      return eContinueLoop;

    if (fCoreType != eFit) {
      showerRec.SetAxis(showerAxis);
      showerRec.SetCorePosition(core);
      showerRec.SetCoreTime(coreTime, TimeInterval(0));
      showerRec.SetAngleErrors(showerAxis.GetCoordinates(gBaryCS), Vector::Triple(0,0,0));
      showerRec.SetCurvature(0, 0);
    }

    const double theta = showerAxis.GetTheta(gBaryCS);
    const double phi = showerAxis.GetPhi(gBaryCS);
    if (theta > fMaxTheta) {
      OUT(eFinal)
        << "  theta = " << theta/degree << " is > than limit "
        << fMaxTheta/degree << " [deg], skipping...\n";
      return eContinueLoop;
    }

    const auto& siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
    const auto& eventTime = sEvent.GetTrigger().GetTime();

    OUT(eIntermediate)
      << "  barycenter = (" << String::Make(fBarycenter, siteCS, meter, ", ") << ") m (site)\n"
         "  bary time = " << (fBaryTime - eventTime).GetInterval()/nanosecond << " ns (to event time)\n"
         "  axis = (" << String::Make(showerAxis, siteCS, 1, ", ") << ") (site)\n"
         "       = (" << String::Make(showerAxis, gBaryCS, 1, ", ") << ") (bary)\n"
         "  theta = " << theta/degree << " deg (bary)\n"
         "  phi = " << phi/degree << " deg (bary)\n";

    OUT(eObscure)
      << "  local/site zenith angle diff = "
      << Angle(Vector(0,0,1, gBaryCS), Vector(0,0,1, siteCS))/degree << " deg\n";

    best.fPar[eCoreX] = core.GetX(gBaryCS);
    best.fPar[eCoreY] = core.GetY(gBaryCS);

    cbest.fPar[eU] = sin(theta) * cos(phi);
    cbest.fPar[eV] = sin(theta) * sin(phi);
    cbest.fPar[eRCore]  = 0;
    cbest.fPar[eCTCore] = 0;
  }

  auto stationFitData = MakeStationFitData(sEvent);
  bool hasSaturatedStation = false;
  for (const auto& station : stationFitData) {
    if (station.fIsSaturated) {
      hasSaturatedStation = true;
      break;
    }
  }

  // --- Stage: fit shower size, LDF shape & core fixed --------------------

  best.fRecStage = ShowerSRecData::kLDFEstimate;
  OUT(eIntermediate)
    << "Stage " << best.fRecStage << ": fit " << GetShowerSizeLabel(config)
    << " - LDF shape & core fixed\n";

  FitLDFSimplified(best, config, cbest.GetAxis(), stationFitData);

  if (nStations < fStage1.fMinNumberOfStations) {
    SetRecData(event, best, cbest);
    OutputResults(event);
    return eSuccess;
  }

  // try with silents and without if first fit fails
  for (int useSilents = fLDFFitConfig.fUseSilentStations; useSilents >= 0; --useSilents) {

    config.fUseSilentStations = bool(useSilents);

    const string stationType = useSilents ? "candidate+silents" : "candidates";

    // --- Stage: fit shower size & core, LDF shape fixed --------------
    {
      LDFFitResult result = best;

      for (size_t i = 0; i < config.fLDF.GetNShapeParameters(); ++i)
        config.fLDFShapeTreatment[i] = eModeled;

      OUT(eIntermediate)
        << "Stage " << result.fRecStage << ": fit "
        << GetShowerSizeLabel(config) << ", core"
        << " - LDF shape fixed [" << stationType << "]\n";

      if (useSilents)
        fBadSilents = CleanSilentStation(result, config, cbest.GetAxis(), stationFitData);

      if (!FitLDF(result, config, cbest.GetAxis(), stationFitData))
        continue;

      result.fRecStage = useSilents ? ShowerSRecData::kLDFSilents : ShowerSRecData::kLDF;
      best = result;

    } // stage

    // --- Stage: shower size, core, shape fitting (optional) ----------
    {
      LDFFitResult result = best;

      config.fLDFShapeTreatment[0] =
        (nStations < fStage2.fMinNumberRelaxBeta ||
         FixBeta(result, cbest.GetAxis(), stationFitData)) ? eModeled : eFitted;
      config.fLDFShapeTreatment[1] =
        (nStations < fStage2.fMinNumberRelaxGamma ||
         FixGamma(result, cbest.GetAxis(), stationFitData)) ? eModeled : eFitted;

      if (config.fLDFShapeTreatment[0] == eFitted ||
          config.fLDFShapeTreatment[1] == eFitted) {
        OUT(eIntermediate)
          << "Stage " << result.fRecStage << ": fit " << stationType
          << " for shower size, core"
          << (config.fLDFShapeTreatment[0] == eFitted ? ", beta" : "")
          << (config.fLDFShapeTreatment[1] == eFitted ? ", gamma\n" : "\n");

        if (FitLDF(result, config, cbest.GetAxis(), stationFitData)) {
          const bool fixBeta = FixBeta(result, cbest.GetAxis(), stationFitData);
          const bool fixGamma = FixGamma(result, cbest.GetAxis(), stationFitData);
          if ((config.fLDFShapeTreatment[0] == eModeled) == fixBeta &&
              (config.fLDFShapeTreatment[1] == eModeled) == fixGamma) {
            result.fRecStage = (useSilents ? ShowerSRecData::kLDFSilents : ShowerSRecData::kLDF) +
                 (config.fLDFShapeTreatment[0] == eFitted) * ShowerSRecData::kLDFBeta +
                 (config.fLDFShapeTreatment[1] == eFitted) * ShowerSRecData::kLDFGamma +
                 (hasSaturatedStation && !config.fUseSaturationRecovery) * ShowerSRecData::kLDFLowerLimit;
            best = result;
            break;
          } else {
            if (fixBeta)
              config.fLDFShapeTreatment[0] = eModeled;
            if (fixGamma)
              config.fLDFShapeTreatment[1] = eModeled;

            if (FitLDF(result, config, cbest.GetAxis(), stationFitData)) {
              result.fRecStage = (useSilents ? ShowerSRecData::kLDFSilents : ShowerSRecData::kLDF) +
                (config.fLDFShapeTreatment[0] == eFitted) * ShowerSRecData::kLDFBeta +
                (config.fLDFShapeTreatment[1] == eFitted) * ShowerSRecData::kLDFGamma +
                (hasSaturatedStation && !config.fUseSaturationRecovery) * ShowerSRecData::kLDFLowerLimit;
              best = result;
              break;
            }
          }
        }
      }

    }

    // if we arrive here: there was no shape fit
    // but the initial fit with fixed shape was ok
    break;
  } // fit LDF with/without silents

  // --- Stage: shower size, core, beta (gamma), curvature fitting -----
  if (nStations < fStage3.fMinNumberOfStations) {
    SetRecData(event, best, cbest);
    OutputResults(event);
    return eSuccess;
  }

  {
    const CurvFitResult cinit = cbest;
    CurvFitResult cres = cbest;
    const double rCoreModel =
      ParameterizedRc(cres.GetAxis().GetCosTheta(gBaryCS),
                      best.fPar[eShowerSize]);

    OUT(eIntermediate)
      << "Fit for axis with fixed curvature radius = " << rCoreModel/kilometer << " km\n";

    fFixedCurvature = true;
    cres.fPar[eRCore] = rCoreModel;
    if (FitCurvature(cres, best.GetCore(), stationFitData)) {
      cbest = cres;

      if (nStations >= fStage3.fMinNumberForFullCurvatureFit) {
        OUT(eIntermediate)
          << "Fit for" << (!fFixedAxis ? " axis and ":" ") << "shower-front curvature.\n";

        fFixedCurvature = false;
        cres.fPar[eRCore] = rCoreModel;
        if (FitCurvature(cres, best.GetCore(), stationFitData))
          cbest = cres;
        else
          OUT(eIntermediate) << "  Curvature fit failed.\n";
      }
    }

    if (fStage3.fOutliersRejection &&
        nStations > fStage3.fMinNumberForFullCurvatureFit) {

      OUT(eIntermediate)
        << "  Performing outlier search and rejection\n";

      for (const auto& station1 : stationFitData) {

        if (!station1.fSignal)
          continue;

        vector<StationFitData> stationData;

        for (const auto& station2 : stationFitData) {
          if (!station2.fSignal || station1.fId == station2.fId)
            continue;
          stationData.push_back(station2);
        }

        cres.fPar[eRCore] = rCoreModel;
        if (FitCurvature(cres, best.GetCore(), stationData)) {

          const double newChi2Probability = utl::Chi2Probability(cres.fChi2, cres.fNdof);
          const double prevChi2Probability = utl::Chi2Probability(cbest.fChi2, cbest.fNdof);

          if (prevChi2Probability < 1e-100 ||
              (newChi2Probability > 1e-100 && newChi2Probability/prevChi2Probability > 2)) {

            OUT(eIntermediate)
              << "  Removed one outlier:\n"
                 "    previous chi2 probability = " << prevChi2Probability << "\n"
                 "    previous RCore = " << cbest.fPar[eRCore]/kilometer << " +/- " << cbest.fCov.Std(eRCore)/kilometer << " km\n"
                 "    previous CTCore = " << cbest.fPar[eCTCore]/meter << " +/- " << cbest.fCov.Std(eCTCore)/meter << " m\n"
                 "    new chi2 probability = " << newChi2Probability << "\n"
                 "    new RCore = " << cres.fPar[eRCore]/kilometer << " +/- " << cres.fCov.Std(eRCore)/kilometer << " km\n"
                 "    new CTCore = " << cres.fPar[eCTCore]/meter << " +/- " << cres.fCov.Std(eCTCore)/meter << " m\n";

            cbest = cres;
          }
        }
      }
    }

    const double openingAngle = Angle(cbest.GetAxis(), cinit.GetAxis());

    OUT(eIntermediate)
      << "  axis = (" << String::Make(cbest.GetAxis(), gBaryCS, meter, ", ") << ") m (bary)\n"
         "  axis diff = " << openingAngle/degree << " deg (to plain-front fit)\n"
         "  RCore = " << cbest.fPar[eRCore]/kilometer << " +/- " << cbest.fCov.Std(eRCore)/kilometer << " km\n"
         "  CTCore = " << cbest.fPar[eCTCore]/meter << " +/- " << cbest.fCov.Std(eCTCore)/meter << " m\n";

    if (openingAngle > fStage3.fMaxAxisAngleDifference || std::isnan(openingAngle)) {
      OUT(eIntermediate)
        << "  Opening-angle axis difference is above the limit.\n";
      cbest = cinit;
    }

    const double newTheta = cbest.GetTheta();
    if (newTheta > fMaxTheta) {
      OUT(eFinal)
        << "  theta = " << newTheta/degree << " deg is > than limit " << fMaxTheta/degree << " deg, skipping...\n";
      return eContinueLoop;
    }
  }

  if (cbest.fPar[eRCore] > 0) {
    // redo LDF fit after the axis change

    const CurvFitResult cinit = cbest;

    LDFFitResult lres = best;

    const bool useSilents = config.fUseSilentStations;
    const string stationType = (useSilents ? "candidate+silents" : "candidates");

    lres.fRecStage =
      (useSilents ? ShowerSRecData::kLDFSilents : ShowerSRecData::kLDF) +
      (config.fLDFShapeTreatment[0] == eFitted) * ShowerSRecData::kLDFBeta +
      (config.fLDFShapeTreatment[1] == eFitted) * ShowerSRecData::kLDFGamma +
      (hasSaturatedStation && !config.fUseSaturationRecovery) * ShowerSRecData::kLDFLowerLimit +
      0.5;

    OUT(eIntermediate)
      << "Redo Stage " << best.fRecStage << ": fit " << GetShowerSizeLabel(config) << ", core"
      << (config.fLDFShapeTreatment[0] == eFitted ? ", beta" : "")
      << (config.fLDFShapeTreatment[1] == eFitted ? ", gamma" : "")
      << " [" << stationType << "]\n";

    if (FitLDF(lres, config, cbest.GetAxis(), stationFitData)) {
      CoreShiftPropagator propagator(best.fPar[eCoreX], best.fPar[eCoreY], cbest.fPar);
      const vector<double> input({lres.fPar[eCoreX], lres.fPar[eCoreY]});

      CovarianceMatrix inputCov(2);
      inputCov[0] = lres.fCov[eCoreX];
      inputCov[1] = lres.fCov[eCoreY];
      inputCov(0, 1) = lres.fCov(eCoreX, eCoreY);

      CovarianceMatrix outputCov(4);
      propagator(cbest.fPar, outputCov, input, inputCov);

      // if some axis properties were not fitted, do not propagate uncertainty
      for (int i = 0; i < 4; ++i)
        for (int j = i; j < 4; ++j)
          if (cbest.fCov(i, j) > 0)
            cbest.fCov(i, j) += outputCov(i, j);

      const double openingAngle = Angle(cbest.GetAxis(), cinit.GetAxis());

      // AS: Check again if new axis is reasonable after redone fit
      // this fails in some rare cases (example: infill 201322403641)
      // in this case an inclined event is wrongly reconstructed as a vertical event!
      if (openingAngle > fStage3.fMaxAxisAngleDifference || std::isnan(openingAngle)) {
        OUT(eIntermediate)
          << "  Opening-angle axis difference is above the limit.\n";
        cbest = cinit;
      } else
        best = lres;

      if (fGlobalFit) {
        lres.fRecStage =
          ShowerSRecData::kLDFGlobalFit +
          (config.fLDFShapeTreatment[0] == eFitted) * ShowerSRecData::kLDFBeta +
          (config.fLDFShapeTreatment[1] == eFitted) * ShowerSRecData::kLDFGamma +
          (hasSaturatedStation && !config.fUseSaturationRecovery) * ShowerSRecData::kLDFLowerLimit;

        OUT(eIntermediate)
          << "Stage " << lres.fRecStage << ": global fit [" << stationType << "]\n";

        CurvFitResult cres = cbest;
        if (FitGlobal(cres, lres, config, stationFitData)) {
          cbest = cres;
          best = lres;
        }
      }
    } // if LDF fit fails, curvature fit is also rejected
  }

  // calculation of Ropt, see GAP-2006-045
  if (fRoptN >= 2) {
    OUT(eIntermediate)
      << "Ropt determination:\n";

    ROptFitter rOptFitter(config.fLDF, best.GetLDFShapeParameters(), fInfoLevel >= eMinuit);

    const double finalBeta = best.fPar[eShapeParameters + 0];
    // beta spread parameterization from Talianna
    const double betaFixSigma = config.fLDF.BetaUncertainty(best.fPar[eShowerSize]);

    const double startp = NormalCDF(-1);
    const double stopp = NormalCDF(1);
    const double deltap = (stopp - startp) / (fRoptN - 1);

    LDFFitConfig thisConfig = config;
    for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i)
      thisConfig.fLDFShapeTreatment[i] = eFixed;

    for (int i = 0; i < fRoptN; ++i) {
      LDFFitResult result = best;
      {
        double& beta = result.fPar[eShapeParameters + 0];
        const double p = startp + i*deltap;
        beta = finalBeta + InverseNormalCDF(p, betaFixSigma);
      }
      FitLDF(result, thisConfig, cbest.GetAxis(), stationFitData);
      const double beta = result.fPar[eShapeParameters + 0];
      rOptFitter.PushBack(beta, result.fPar[eShowerSize]);
    }

    // simple numerical derivative (finite difference)
    const auto halfN = fRoptN / 2;
    const double s1 = rOptFitter.GetShowerSizeVector()[halfN-1];
    const double s2 = rOptFitter.GetShowerSizeVector()[halfN];
    const double b1 = rOptFitter.GetBetaVector()[halfN-1];
    const double b2 = rOptFitter.GetBetaVector()[halfN];
    const double dS1000_dBeta = (s2 - s1) / (b2 - b1);

    auto& shRec = event.GetRecShower().GetSRecShower();
    shRec.SetBetaSystematics(betaFixSigma);
    shRec.SetShowerSizeSystematics(fabs(betaFixSigma * dS1000_dBeta));

    double rOpt = 0;
    if (rOptFitter.GetResult(rOpt))
      shRec.SetLDFOptimumDistance(rOpt);
  }

  SetRecData(event, best, cbest);
  OutputResults(event);

  return eSuccess;
}


bool
LDFFinder::FixCore(Point& core, TimeStamp& refTime, Vector& axis,
                   const evt::Event& event)
  const
{
  const Vector localNormal(0,0,1, gBaryCS);
  Point coreReference;

  switch (fCoreType) {
  case eFit:
    {
      core = fBarycenter;
      refTime = fBaryTime;
      const auto& previousAxis = event.GetRecShower().GetSRecShower().GetAxis();
      if (!previousAxis) {
        ERROR("No previous axis found.");
        return false;
      }
      axis = previousAxis;
    }
    break;
  case eMC:
    {
      if (!event.HasSimShower()) {
        ERROR("The event is not a simulated shower.");
        return false;
      }
      const auto& simShower = event.GetSimShower();
      axis = -simShower.GetDirection();
      coreReference = simShower.GetPosition();
      refTime = simShower.GetTimeStamp();
    }
    break;
  case eHybrid:
    {
      if (!event.HasFEvent()) {
        ERROR("The event has no FD shower.");
        return false;
      }

      const auto& fEvent = event.GetFEvent();
      // remember area of FD core error ellipse to choose best FD eye
      double bestErrorEllipse = 0;
      const evt::ShowerFRecData* showerFRec = nullptr;

      for (const auto& eye :
           boost::make_iterator_range(fEvent.EyesBegin(ComponentSelector::eHasData),
                                      fEvent.EyesEnd(ComponentSelector::eHasData))) {

        // skip eyes without reconstruction
        if (!eye.HasRecData())
          continue;

        const auto& eyeRec = eye.GetRecData();
        if (eyeRec.GetSDPFitNDof() <= 0 || eyeRec.GetTimeFitNDof() <= 0)
          continue;

        if (!eyeRec.HasFRecShower())
          continue;
        showerFRec = &eyeRec.GetFRecShower();

        const auto& dEye = det::Detector::GetInstance().GetFDetector().GetEye(eye);
        const auto& eyeCS = dEye.GetEyeCoordinateSystem();
        const auto eyePos = dEye.GetPosition();
        const auto vecEyeCore = eyePos - showerFRec->GetCorePosition();

        const double rp = eyeRec.GetRp();
        const double rpErr = eyeRec.GetRpError();
        const double t0 = eyeRec.GetTZero();
        const double chi0 = eyeRec.GetChiZero();
        const double chi0Err = eyeRec.GetChiZeroError();
        const double sdpPhiErr = eyeRec.GetSDPPhiError();
        const double rhoChi0Rp = eyeRec.GetRpChi0Correlation();
        const double distEyeCore = vecEyeCore.GetMag();

        // calculation of core error in SDP
        // rp error projected on ground+
        // chi0 error+ correlation of errors in rp and chi0
        const double sinChi0 = sin(chi0);
        const double errorInSDP = sqrt(
          Sqr(rpErr / sinChi0) +
          Sqr((chi0Err * distEyeCore) / (cos(chi0) * sinChi0)) +
          rhoChi0Rp * rpErr * chi0Err
        );

        // error due to SDP azimuth uncertainty
        const double errorPerpSDP = distEyeCore * sin(sdpPhiErr);

        // area of core error ellipse of this FD event
        const double coreErrorEllipse = errorInSDP * errorPerpSDP * kPi;

        // if error ellipse set AND error ellipse larger than previous eye(s)
        // don't update core with actual eye
        if (bestErrorEllipse && bestErrorEllipse < coreErrorEllipse)
          continue;

        bestErrorEllipse = coreErrorEllipse;
        axis = showerFRec->GetAxis();
        coreReference = showerFRec->GetCorePosition();

        const Vector vertical(0,0,1, eyeCS);
        const double alpha = 0.5*kPi - Angle(vertical, vecEyeCore);
        const double cosBeta = CosAngle(axis, vecEyeCore);
        // 1e5ns have to be subtracted to jump from the end to the
        // start of the pixel trace (don't even think of removing the
        // comment here)
        refTime = eye.GetHeader().GetTimeStamp() +
          TimeInterval(-1e5*ns + t0 + (rp*tan(alpha) + distEyeCore*cosBeta)/kSpeedOfLight);
      } // loop over eyes

      if (!showerFRec)
        return false;
    }
    break;
  }

  if (fCoreType != eFit) {
    // shift core along shower axis into local (barycenter) ground plane
    core = coreReference -
      (localNormal * (coreReference - fBarycenter)) / (localNormal*axis) * axis;

    const double d = (core - fBarycenter).GetMag();
    if (d > fMaxBaryToCoreDistance) {
      ostringstream err;
      err << "Bary-centric distance to core " << d/kilometer << " km is above the limit.";
      ERROR(err);
      return false;
    }

    OUT(eIntermediate)
      << "  using fixed core = (" << String::Make(core, gBaryCS, meter, ", ") << ") m (bary)\n";
  } else
    OUT(eIntermediate)
      << "  using barycenter as initial core\n";

  return true;
}


// From Pierre Billoir Thoughts....
// - if there are enough stations (especially at distances around 1000 m), a fit
// with a variable LDF slope is attempted (B is set as adjustable in (3)); the
// conditions are:
// * at least 5 stations, with one of the following requirements:
// * at least 2 within (400, 1600), with a difference at least 900 in between
// * at least 3 within (400, 1600), with a maximal difference at least 800
// * at least 4 within (400, 1600), with a maximal difference at least 700
//
// New version much stricter
//
// Adjustments made for the SD750 and SD433
//
bool
LDFFinder::FixBeta(const LDFFitResult& result,
                   const Vector& showerAxis,
                   const vector<StationFitData>& data)
  const
{
  size_t nStations = 0;
  for (const auto& station : data) {
    if (!station.fSignal || (station.fIsSaturated && !station.fRecoveredSignalErr))
      continue;
    ++nStations;
  }

  if (nStations < 5) {
    OUT(eIntermediate)
      << "FixBeta = 1: nStations = " << nStations << '\n';
    return true;
  }

  int nStationsInRange = 0;
  double maxDistance = 0;

  const auto core = result.GetCore();

  for (const auto& station : data) {

    if (!station.fSignal || (station.fIsSaturated && !station.fRecoveredSignalErr))
      continue;

    const double r = RPerp(showerAxis, station.fPos - core);

    if (fStage2.fFixBetaRange[0] < r && r < fStage2.fFixBetaRange[1]) {
      ++nStationsInRange;

      for (const auto& s : data) {
        const double rr = RPerp(showerAxis, s.fPos - core);

        if (fStage2.fFixBetaRange[0] < rr && rr < fStage2.fFixBetaRange[1]) {
          const double distance = fabs(r - rr);
          if (distance > maxDistance)
            maxDistance = distance;
        }
      }
    }
  }

  const bool fixBeta = !(
    (nStationsInRange > 1 && maxDistance > fStage2.fFixBetaLeverArmRange[0]) ||
    (nStationsInRange > 2 && maxDistance > fStage2.fFixBetaLeverArmRange[1]) ||
    (nStationsInRange > 3 && maxDistance > fStage2.fFixBetaLeverArmRange[2])
  );

  OUT(eIntermediate)
    << "FixBeta = " << fixBeta << ": "
       "nStations (InRange) = " << nStations << " (" << nStationsInRange << "), "
       "maxDistance = " << maxDistance/meter << " m\n";

  return fixBeta;
}


bool
LDFFinder::FixGamma(const LDFFitResult& result,
                    const Vector& showerAxis,
                    const vector<StationFitData>& data)
  const
{
  size_t nStations = 0;
  for (const auto& station : data) {
    if (!station.fSignal || (station.fIsSaturated && !station.fRecoveredSignalErr))
      continue;
    ++nStations;
  }

  if (nStations < 5) {
    OUT(eIntermediate)
      << "FixGamma = 1: nStations = " << nStations << '\n';
    return true;
  }

  int nStationsInRange = 0;
  double maxDistance = 0;
  const auto core = result.GetCore();
  for (const auto& station : data) {
    if (!station.fSignal || (station.fIsSaturated && !station.fRecoveredSignalErr))
      continue;

    const double r = RPerp(showerAxis, station.fPos - core);

    if (fStage2.fFixGammaRange[0] < r && r < fStage2.fFixGammaRange[1]) {

      ++nStationsInRange;

      for (const auto& s : data) {

        const double rr = RPerp(showerAxis, s.fPos - core);

        if (fStage2.fFixGammaRange[0] < rr && rr < fStage2.fFixGammaRange[1]) {
          const double distance = fabs(r - rr);
          if (distance > maxDistance)
            maxDistance = distance;
        }
      }
    }
  }

  const bool fixGamma = !(
    (nStationsInRange > 1 && maxDistance > fStage2.fFixGammaLeverArmRange[0]) ||
    (nStationsInRange > 2 && maxDistance > fStage2.fFixGammaLeverArmRange[1]) ||
    (nStationsInRange > 3 && maxDistance > fStage2.fFixGammaLeverArmRange[2])
  );

  OUT(eIntermediate)
    << "FixGamma = " << fixGamma << ": "
       "nStations (InRange) = " << nStations << " (" << nStationsInRange << "), "
       "maxDistance = " << maxDistance/meter << " m\n";

  return fixGamma;
}


void
LDFFinder::SetRecData(Event& event,
                      const LDFFitResult& lres,
                      const CurvFitResult& cres)
  const
{
  auto& recShower = event.GetRecShower().GetSRecShower();

  const auto& sDetector = det::Detector::GetInstance().GetSDetector();

  const auto& ldf = fLDFFitConfig.fLDF;
  const double showerSize = lres.fPar[eShowerSize];
  const auto& ldfShapePars = lres.GetLDFShapeParameters();

  const auto& core = lres.GetCore();
  const auto& coreCS = LocalCoordinateSystem::Create(core);
  const auto& axis = cres.GetAxis();
  const double cosTheta = axis.GetCosTheta(coreCS);

  auto& sEvent = event.GetSEvent();

  const double rCurv = cres.fPar[eRCore];
  const bool haveCurvature = (rCurv != 0);
  const Point showerOrigin = core + rCurv * axis;
  const double cTCore = cres.fPar[eCTCore];

  const double theta = axis.GetTheta(coreCS);
  const double phi = axis.GetPhi(coreCS);
  const CoordinateSystemPtr showerPlaneCS(
    CoordinateSystem::RotationY(
      theta,
      CoordinateSystem::RotationZ(
        phi,
        coreCS
      )
    )
  );

  double coreCovSP[2][2]; // core covariance in shower plane coordinates
  {
    // uncertainty propagation for trafo rSP = RM r (RM = Rotation Matrix)
    // see Documentation/ReferenceGuide/SdReconstruction
    HepRotation r;
    r.rotateZ(-phi);
    r.rotateY(-theta);

    for (int i : { 0, 1 }) {
      const auto& r_i = r[i];
      auto& cc_i = coreCovSP[i];
      for (int j : { 0, 1 }) {
        const auto& r_j = r[j];
        auto& cc_ij = cc_i[j];
        cc_ij = 0;
        for (int k : { 0, 1 }) {
          const auto& r_ik = r_i[k];
          for (int m : { 0, 1 })
            cc_ij += r_ik * r_j[m] * lres.fCov(1+k, 1+m);
        }
      }
    }
  }

  const auto& timeVariance = sdet::STimeVariance::GetInstance();
  Accumulator::MinMax<double> minMaxRho(1000*kilometer, 0);
  Accumulator::SampleStandardDeviation timeStat;

  for (auto& station : sEvent.StationsRange()) {
    if (!station.HasRecData())
      continue;

    const auto& sPos = sDetector.GetStation(station).GetPosition();
    const double x = sPos.GetX(showerPlaneCS);
    const double y = sPos.GetY(showerPlaneCS);
    const double rho = std::sqrt(x*x + y*y);

    if (station.IsCandidate())
      minMaxRho(rho);

    const double rhoErr =
      std::sqrt(Sqr(x) * coreCovSP[0][0] + Sqr(y) * coreCovSP[1][1] + 2*x*y * coreCovSP[0][1]) / rho;

    auto& sRec = station.GetRecData();
    const double signal = sRec.GetTotalSignal();
    const double funcval = showerSize * ldf(rho, ldfShapePars);
    const double sigma = utl::wcd::SignalUncertainty(fLDFFitConfig.fSignalVarianceModel, cosTheta, funcval);
    sRec.SetAzimuthShowerPlane(sPos.GetPhi(showerPlaneCS));
    sRec.SetTotalSignal(signal, sigma);
    sRec.SetLDFResidual((signal - funcval)/sigma);
    sRec.SetSPDistance(rho, rhoErr);

    if (haveCurvature) {
      const double measuredTime = (sRec.GetSignalStartTime() - fBaryTime).GetInterval();
      const Vector stationToOrigin = showerOrigin - sPos;
      const double rStation = stationToOrigin.GetMag();
      const double predictedTime = (rStation + cTCore - rCurv) / kSpeedOfLight;
      const double cosThetaStation = stationToOrigin.GetZ(coreCS) / rStation;
      double sigma2 = GetTimeVariance(fLDFFitConfig, timeVariance, signal, sRec, cosThetaStation, rho);
      if (fLDFFitConfig.fUseAdditionalGPSTimeVariance)
        sigma2 += sRec.GetGPSTimeVariance();
      const double residual = (measuredTime - predictedTime) / std::sqrt(sigma2);
      sRec.SetResidual(residual);
      if (station.IsCandidate())
        timeStat(residual);
    }

  } // loop over stations with RecData

  if (!recShower.HasLDF())
    recShower.MakeLDF();
  auto& tabulatedLDF = recShower.GetLDF();
  tabulatedLDF.Clear();

  {
    const double x0 = log(1);
    const double x1 = log(10000);
    const int nLDFPoints = 20;

    for (int i = 0; i < nLDFPoints; ++i) {
      const double z = double(i) / (nLDFPoints - 1);
      const double rho = exp((1 - z)*x0 + z*x1);
      const double funcval = showerSize * ldf(rho, ldfShapePars);
      const double sigma = utl::wcd::SignalUncertainty(fLDFFitConfig.fSignalVarianceModel, cosTheta, funcval);
      tabulatedLDF.PushBack(rho, 0, funcval, sigma);
    }
  }

  // set the reconstruction data
  if (haveCurvature) {
    const double rCurvErr = cres.fCov.Std(eRCore);
    const double cTOffsetErr = cres.fCov.Std(eCTCore);
    const double cTOffset = cres.fPar[eCTCore];
    recShower.SetCurvature(1/rCurv, rCurvErr/Sqr(rCurv));
    recShower.SetCurvatureTimeOffset(cTOffset, cTOffsetErr);
    recShower.SetCoreTime(fBaryTime + TimeInterval(cTCore/kSpeedOfLight),
                          TimeInterval(cTOffsetErr/kSpeedOfLight));

    const double timeMean = timeStat.GetAverage();
    const double timeSpread = timeStat.GetStandardDeviation();
    recShower.SetTimeResidualMean(timeMean);
    recShower.SetTimeResidualSpread(timeSpread);

    recShower.SetAngleChi2(cres.fChi2, cres.fNdof);
    const double u = cres.fPar[eU];
    const double v = cres.fPar[eV];
    const double w = sqrt(1 - u*u - v*v);
    const double suu = cres.fCov[eU];
    const double suv = cres.fCov(eU, eV);
    const double svv = cres.fCov[eV];

    const double suw = -(suu*u + suv*v) / w;
    const double svw = -(suv*u + svv*v) / w;
    const double sww = (u*u*suu + 2*u*v*suv + v*v*svv) / (w*w);

    recShower.SetParameter(eShowerAxisX, u);
    recShower.SetParameter(eShowerAxisY, v);
    recShower.SetParameter(eShowerAxisZ, w);
    recShower.SetParameterCovariance(eShowerAxisX, eShowerAxisX, suu);
    recShower.SetParameterCovariance(eShowerAxisY, eShowerAxisY, svv);
    recShower.SetParameterCovariance(eShowerAxisZ, eShowerAxisZ, sww);
    recShower.SetParameterCovariance(eShowerAxisX, eShowerAxisY, suv);
    recShower.SetParameterCovariance(eShowerAxisX, eShowerAxisZ, suw);
    recShower.SetParameterCovariance(eShowerAxisY, eShowerAxisZ, svw);

    recShower.SetAxis(axis);
    recShower.SetAngleErrors(Vector::Triple(u, v, w), Vector::Triple(suu, suv, svv));

  } else {
    // we have moved the core but have axis from plane fit, recalc t0
    const TimeInterval cdiff = recShower.GetAxis() * (recShower.GetCorePosition() - core);
    const TimeStamp newCoreTime = recShower.GetCoreTime() + cdiff/kSpeedOfLight;
    recShower.SetCurvature(0, 0);
    recShower.SetCoreTime(newCoreTime, recShower.GetCoreTimeError());
  }

  {
    const double showerSizeErr = lres.fCov.Std(eShowerSize);
    vector<double> ldfShapeParsErr;
    for (size_t i = 0, n = fLDFFitConfig.fLDF.GetNShapeParameters(); i < n; ++i)
      ldfShapeParsErr.push_back(lres.fCov.Std(eShapeParameters + i));
    recShower.SetNumberOfActiveStations(timeStat.GetN());

    switch (fArrayType) {
    case sdet::Array::e1500:
      recShower.SetShowerSizeType(ShowerSRecData::eS1000);
      break;
    case sdet::Array::e750:
      recShower.SetShowerSizeType(ShowerSRecData::eS450);
      break;
    case sdet::Array::e433:
      recShower.SetShowerSizeType(ShowerSRecData::eS300);
      break;
    default:
      recShower.SetShowerSizeType(ShowerSRecData::eUndefined);
    }

    recShower.SetShowerSize(showerSize, showerSizeErr);
    recShower.SetCorePosition(core);
    const double coreXErr = lres.fCov.Std(eCoreX);
    const double coreYErr = lres.fCov.Std(eCoreY);
    recShower.SetCoreError(Vector(coreXErr, coreYErr, 0, gBaryCS));
    recShower.SetCorrelationXY((coreXErr > 0 && coreYErr > 0) ?
                               lres.fCov(eCoreX,eCoreY)/(coreXErr*coreYErr) : 0.);
    recShower.SetBeta(ldfShapePars[0], ldfShapeParsErr[0]);
    if (ldfShapePars.size() > 1)
      recShower.SetGamma(ldfShapePars[1], ldfShapeParsErr[1]);
    else
      recShower.SetGamma(0, 0);
    if (fLDFFitConfig.fLDF.GetType() == "NKGFermi")
      recShower.SetNKGFermiParameters(ldfShapePars[2], ldfShapePars[3]);
    recShower.SetLDFChi2(lres.fChi2, lres.fNdof);
    recShower.SetLDFLikelihood(lres.fLogLikelihood);

    recShower.SetCicRefAngle(fEnergyConversion.fCicReferenceAngle);
    recShower.SetS38(fEnergyConversion.GetS38(showerSize, cosTheta));

    double energy = 0;
    double energyErr = 0;
    double energySys = 0;
    fEnergyConversion.GetEnergy(cosTheta, showerSize, showerSizeErr,
                                recShower.GetShowerSizeSystematics(),
                                energy, energyErr, energySys);

    if (!energy)
      OUT(eFinal)
        << "  theta = " << theta/degree << " is outside the range of the CIC. The energy is not valid.\n";

    recShower.SetEnergy(energy, energyErr);
    recShower.SetEnergyLDFSystematics(energySys);
    recShower.SetLDFRecStage(lres.fRecStage);

    for (const auto sId : fBadSilents) {
      if (sEvent.HasStation(sId))
        sEvent.GetStation(sId).SetRejected(StationConstants::eBadSilent);
    }

  }
}


vector<StationFitData>
LDFFinder::MakeStationFitData(const SEvent& sEvent)
  const
{
  vector<StationFitData> sFitData;

  const auto& sDetector = det::Detector::GetInstance().GetSDetector();

  for (const auto& station : sEvent.StationsRange()) {
    StationFitData current;
    if (station.IsCandidate()) {
      current.fId = station.GetId();
      current.fPos = sDetector.GetStation(station).GetPosition();
      current.fIsSaturated = station.IsLowGainSaturation();
      const auto& sRec = station.GetRecData();
      current.fSignal = sRec.GetTotalSignal();
      current.fRecoveredSignal = sRec.GetRecoveredSignal();
      current.fRecoveredSignalErr = sRec.GetRecoveredSignalError();
      current.fCTime = kSpeedOfLight * (sRec.GetSignalStartTime() - fBaryTime).GetInterval();
      current.fT50 = sRec.GetT50();
      current.fGPSTimeVariance = fLDFFitConfig.fUseAdditionalGPSTimeVariance ? sRec.GetGPSTimeVariance() : 0;
      sFitData.push_back(current);
    } else if (station.IsSilent() ||
               (fUseSilentStationsFromOtherGrids &&
                station.GetRejectionStatus() == StationConstants::eOffGrid &&
                !station.HasRecData() &&
                sDetector.GetStation(station).IsInGrid(sdet::SDetectorConstants::eAny))) {
      current.fId = station.GetId();
      current.fPos = sDetector.GetStation(station).GetPosition();
      current.fIsSaturated = false;
      current.fSignal = 0;
      current.fRecoveredSignalErr = 0;
      current.fCTime = 0;
      current.fT50 = 0;
      current.fGPSTimeVariance = 0;
      sFitData.push_back(current);
    }
  } // loop over stations

  return sFitData;
}


void
LDFFinder::FitLDFSimplified(LDFFitResult& result,
                            LDFFitConfig& config,
                            const Vector& showerAxis,
                            const vector<StationFitData>& data)
  const
{
  result.fCov = 0;

  const auto core = result.GetCore();

  const double cosTheta = showerAxis.GetCosTheta(gBaryCS);

  const auto shape = config.fLDF.ShapeModel(cosTheta, 10);

  LDFLikelihoodFunction likeFcn(config, showerAxis, data);

  LDFFitResult resultA = result;
  LDFFitResult resultB = result;

  {
    // estimate shower size from station closest to reference distance

    const double kRefDist = config.fLDF.GetReferenceDistance();
    double minDistance = numeric_limits<double>::max();
    double rhoClosest = kRefDist;
    double signalClosest = 1;

    // loop over normal stations, ignore silent/saturated stations
    for (const auto& station : data) {
      if (!station.fSignal)
        continue;

      const double rho = RPerp(showerAxis, station.fPos - core);
      const double distance = fabs(rho - kRefDist);

      if (distance < minDistance) {
        minDistance = distance;
        rhoClosest = rho;
        signalClosest = station.fSignal;
      }
    }

    resultA.fPar[eShowerSize] = signalClosest / config.fLDF(rhoClosest, shape);
    resultA.fLogLikelihood = likeFcn(resultA.fPar);

    OUT(eIntermediate)
      << "  " << GetShowerSizeLabel(config) << " [estimate A] = "
      << resultA.fPar[eShowerSize] << " VEM (logLike = " << resultA.fLogLikelihood << ")\n";
  }

  {
    /*
      If core and shape of the LDF are fixed, the least-squares
      estimate of shower size S0 can be calculated analytically.

      S_i                ... data
      S(r_i) = S0 f(r_i) ... model
      sigma^2(S) = S     ... simplest signal uncertainty model

      chi^2 = sum_i { (S_i - S(r_i))^2 / S(r_i) }

      With maximum condition d chi^2/d S0 = 0 one obtains
      S0^2 = sum_i { S_i^2 / f(r_i)) } / sum_i { f(r_i) }
    */

    double sumA = 0;
    double sumB = 0;
    for (const auto& station : data) {
      const double s = station.fSignal;

      if (!s)
        continue;

      const double rho = RPerp(showerAxis, station.fPos - core);
      const double f = config.fLDF(rho, shape);

      sumA += s*s / f;
      sumB += f;
    }

    resultB.fPar[eShowerSize] = sqrt(sumA / sumB);
    resultB.fLogLikelihood = likeFcn(resultB.fPar);

    OUT(eIntermediate)
      << "  " << GetShowerSizeLabel(config) << " [estimate B] = "
      << resultB.fPar[eShowerSize] << " VEM (logLike = " << resultB.fLogLikelihood << ")\n";
  }

  // use better estimate of the two
  result = resultA.fLogLikelihood < resultB.fLogLikelihood ? resultA : resultB;

  // Finally, perform full-featured fit with fixed core
  LDFFitConfig thisConfig = config;
  thisConfig.fFixCore = true;
  for (size_t i = 0, n = thisConfig.fLDF.GetNShapeParameters(); i < n; ++i)
    thisConfig.fLDFShapeTreatment[i] = eModeled;
  FitLDF(result, thisConfig, showerAxis, data);
  config.fUseSaturationRecovery = thisConfig.fUseSaturationRecovery; // remember to possibly ignore sat. rec.
}


bool
LDFFinder::FitLDF(LDFFitResult& result,
                  LDFFitConfig& config /* needs to be mutable */,
                  const Vector& showerAxis,
                  const vector<StationFitData>& data)
  const
{
  using namespace ROOT::Minuit2;

  MnUserParameters pars;

  pars.Add("showerSize", result.fPar[0], 1.0);
  pars.SetLowerLimit("showerSize", 0);
  pars.Add("coreX", result.fPar[1], 200*meter);
  pars.Add("coreY", result.fPar[2], 200*meter);

  for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i) {
    pars.Add(GetShapeLabel(i), result.fPar[3+i], 0.1);
    if (config.fLDFShapeTreatment[i] != eFitted)
      pars.Fix(GetShapeLabel(i));
  }

  if (config.fFixCore) {
    pars.Fix("coreX");
    pars.Fix("coreY");
  }

  LDFChi2Function chi2Fcn(config, showerAxis, data);
  LDFLikelihoodFunction likeFcn(config, showerAxis, data);

  MnMinimize m(fUseMaxLikelihood ?
               static_cast<const FCNBase&>(likeFcn) :
               static_cast<const FCNBase&>(chi2Fcn), pars, 1);

  FunctionMinimum fmin = m(10000);

  if (fmin.IsAboveMaxEdm() || fmin.HasReachedCallLimit() || !fmin.IsValid()) {
    OUT(eIntermediate)
      << "\n**************** Minimize FAILED ****************\n"
      << fmin
      << "**************** Minimize FAILED ****************\n\n";
    return false;
  }

  {
    ostringstream msg;
    msg << "Found minimum after " << fmin.NFcn() << " evaluations";
    INFO(msg);
  }

  if (config.fUseSaturationRecovery) {
    // was use of saturation recovery ok? fit again without if
    // a) second derivative of LDF at sat. station > threshold,
    //    but only for events with more than three candidates
    // b) signal from LDF < (unrecovered signal - 5 sigma) of sat. station,
    //    such a fit would be inconsistent
    size_t nCandidates = 0;
    for (const auto& station : data) {
      if (station.fSignal > 0)
        ++nCandidates;
    }

    const auto& params = fmin.UserParameters().Params();
    const Point core(params[eCoreX], params[eCoreY], 0, gBaryCS);
    auto shape = config.fLDF.ShapeModel(showerAxis.GetCosTheta(gBaryCS), params[eShowerSize]);
    for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i)
      if (config.fLDFShapeTreatment[i] != eModeled)
        shape[i] = params[eShapeParameters + i];

    for (const auto& station : data) {
      if (station.fIsSaturated) {
        if (station.fRecoveredSignal > 0 && station.fRecoveredSignalErr > 0) {
          const Vector coreToStation = station.fPos - core;
          const double rho = RPerp(showerAxis, coreToStation);

          const double secDerivative = config.fLDF.SecondDerivative(rho, shape);
          const double prediction = params[eShowerSize]*config.fLDF(rho, shape);
          if ((nCandidates > 3 && secDerivative > config.fRecoveryThreshold) ||
              prediction < (station.fSignal - 5*sqrt(prediction))) {
            ostringstream msg;
            msg << "Saturation recovery rejected - fit again without\n"
                   "  2nd derivative of LDF = " << secDerivative << " [Max: " << config.fRecoveryThreshold << "]\n"
                   "  LDF prediction = " << prediction << " VEM [Min: " << (station.fSignal - 5*sqrt(prediction)) << " VEM]\n"
                   "  " << GetShowerSizeLabel(config) << " from previous minimization = " << params[eShowerSize] << " VEM";
            INFO(msg);

            config.fUseSaturationRecovery = false;

            fmin = m(10000);

            msg.str("");
            msg << "Found minimum after " << fmin.NFcn() << " evaluations";
            INFO(msg);

            break;
          }
        } else
          // remember for later that no saturation recovery was available
          config.fUseSaturationRecovery = false;
      }
    }
  }

  MnHesse hesse(1);
  hesse(m.Fcnbase(), fmin, 10000);

  if (!fmin.IsValid() || fmin.HesseFailed()) {
    OUT(eIntermediate)
      << "\n**************** Hesse FAILED ****************\n"
      << fmin
      << "**************** Hesse FAILED ****************\n\n";
    return false;
  }

  if (fInfoLevel >= eMinuit)
    OUT(eIntermediate)
      << "\n**************** Minuit Output ****************\n"
      << fmin
      << "**************** Minuit Output ****************\n\n";

  // get fit results
  result.fPar = fmin.UserParameters().Params();

  const auto modeledShape =
    config.fLDF.ShapeModel(showerAxis.GetCosTheta(gBaryCS), result.fPar[0]);
  for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i)
    if (config.fLDFShapeTreatment[i] == eModeled)
      result.fPar[3+i] = modeledShape[i];

  const auto nVar = pars.VariableParameters();

  // get covariance
  result.fCov = 0;

  for (size_t i = 0; i < nVar; ++i)
    for (size_t j = 0; j < nVar; ++j)
      result.fCov(m.ExtOfInt(i), m.ExtOfInt(j)) = fmin.UserCovariance()(i,j);

  const double minValue = fmin.Fval();
  if (fUseMaxLikelihood) {
    result.fLogLikelihood = minValue;
    bool save = config.fUseSilentStations;
    config.fUseSilentStations = false;
    result.fChi2 = chi2Fcn(result.fPar);
    config.fUseSilentStations = save;
  } else {
    result.fLogLikelihood = likeFcn(result.fPar);
    result.fChi2 = minValue;
  }

  result.fNdof = 0;
  for (const auto& station : data)
    if (station.fSignal > 0)
      ++result.fNdof;
  result.fNdof -= nVar;

  const double reducedChi2 = (result.fNdof > 0) ? result.fChi2/result.fNdof : 0;

  const auto core = result.GetCore();

  const utl::Point coreErr(result.fCov.Std(eCoreX), result.fCov.Std(eCoreY), 0, gBaryCS);
  OUT(eIntermediate)
    << "  chi2 = " << result.fChi2 << " / " << result.fNdof << "\n"
       "  " << GetShowerSizeLabel(config) << " = "
            << result.fPar[eShowerSize] << " +/- " << result.fCov.Std(eShowerSize) << " VEM\n"
       "  core = (" << String::Make(core, gBaryCS, meter, ", ") << ") +/- (" << String::Make(coreErr, gBaryCS, meter, ", ") << ") m (bary)";
  for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i)
    OUT(eIntermediate) << "\n"
      "  " << GetShapeLabel(i) << " = " << result.fPar[eShapeParameters + i] << " +/- " << result.fCov.Std(eShapeParameters + i);
  OUT(eIntermediate) << '\n';

  if (reducedChi2 > fMaxChi2) {
    OUT(eIntermediate)
      << "  chi2 / Ndof = " << reducedChi2 << " exceeds limit " << fMaxChi2 << '\n';
    return false;
  }

  const double d = (core - fBarycenter).GetMag();
  if (d > fMaxBaryToCoreDistance) {
    OUT(eIntermediate)
      << "  Distance to core " << d/km << " km exceeds limit.\n";
    return false;
  }

  return !HasNaN(result);
}


double
LDFFinder::ParameterizedRc(const double cosTheta, const double showerSize)
  const
{
  // TODO: move the parametrization of the radius of curvature to the configuration files
  // TODO: implement a parametrization of the radius of curvature for the SD433

  const double sec1 = 1/cosTheta - 1;
//#warning: IM Curvature parametrisation does not work with very inclined showers

  if (fArrayType == sdet::Array::e750 || fArrayType == sdet::Array::e433) {
    const double x = Sqr(cosTheta)- Sqr(0.81932);
    const double attenuationFactor = 1 + x*(1.55 - 1.14*x);  // Alex attenuation (Feb 2012)
    const double s35 = (attenuationFactor > 0) ? showerSize/attenuationFactor : showerSize;
    const double s = log(s35);
    const double eta = 0.1963; //+-0.002
    const double mu = -0.0135; //+-0.0007
    const double zeta = 0.0002866; //+- 7.6e-05
    const double sDependence = eta + s*(mu + zeta*s);

    const double a = -1.068; //+-0.006
    const double b = 0.5059; //+-0.007
    const double thetaDependence = 1 + sec1*(a + b*sec1);

    const double curvature = sDependence * thetaDependence / kilometer;

    return curvature ? 1/curvature : 0;
  } else {
    const double x = Sqr(cosTheta) - Sqr(0.78821);
    const double attenuationFactor = 1 + x*(0.87 - 1.49*x);  // ICRC2011 attenuation
    const double s38 = (attenuationFactor > 0) ? showerSize/attenuationFactor : showerSize;
    const double s = log(s38);
    const double eta = 0.1785; //+-  0.001408
    const double mu = -0.02149; //+- 0.0009256
    const double zeta = 0.0016; //+- 0.000149
    const double sDependence = eta + s*(mu + zeta*s);

    const double a = -1.047; //+- 0.006598
    const double b =  0.473; //+- 0.007029
    const double thetaDependence = 1 + sec1*(a + b*sec1);

    const double curvature = sDependence * thetaDependence / kilometer;

    return curvature ? 1/curvature : 0;

    /* old parametrization
    const double offset = 9.03 - 3.08*s + 1.08*s2;
    const double slope = 7.28 + 4.03*s - 1.01*s2;

    const double s = log10(showerSize);
    const double s2 = Sqr(s);
    const double offset = 7.418 - 3.0613*s + 2.5262*s2;
    const double slope = 3.9779 + 3.2245*s + 4.5705*s2;

    return (offset + slope*sec1) * kilometer;*/
  }
}


double
LDFFinder::FitCurvature(CurvFitResult& result,
                        const Point& core,
                        const vector<StationFitData>& data)
  const
{
  using namespace ROOT::Minuit2;

  MnUserParameters pars;

  // MINUIT requires unbounded variables, but
  // u and v are bounded, so I define new variables
  // pu and pv here that are unbounded.
  // The mapping is one-on-one and has the property
  // pu ~ u and pv ~ v for small u and v, thus the
  // non-linear properties appear only for very inclined events.
  const double u = result.fPar[eU];
  const double v = result.fPar[eV];
  const double w = sqrt(1 - u*u - v*v);
  pars.Add("pu", u/w, 0.01);
  pars.Add("pv", v/w, 0.01);
  pars.Add("rCore", result.fPar[eRCore], 10*meter);
  pars.SetLowerLimit("rCore", 0);
  pars.Add("cTCore", result.fPar[eCTCore], 10*meter);

  if (fFixedAxis) {
    pars.Fix("pu");
    pars.Fix("pv");
  }

  if (fFixedCurvature)
    pars.Fix("rCore");

  CurvatureChi2Function chi2Fcn(core, data);

  MnMinimize m(chi2Fcn, pars, 1);

  FunctionMinimum fmin = m(10000);

  if (fmin.IsAboveMaxEdm() || fmin.HasReachedCallLimit() || !fmin.IsValid()) {
    OUT(eIntermediate)
      << "\n**************** Minimize FAILED ****************\n"
      << fmin
      << "**************** Minimize FAILED ****************\n\n";
    return false;
  }

  {
    ostringstream msg;
    msg << "Found minimum after " << fmin.NFcn() << " evaluations";
    INFO(msg);
  }

  MnHesse hesse(1);
  hesse(m.Fcnbase(), fmin, 10000);

  if (!fmin.IsValid() || fmin.HesseFailed()) {
    OUT(eIntermediate)
      << "\n**************** Hesse FAILED ****************\n"
      << fmin
      << "**************** Hesse FAILED ****************\n\n";
    return false;
  }

  if (fInfoLevel >= eMinuit)
    OUT(eIntermediate)
      << "\n**************** Minuit Output ****************\n"
      << fmin
      << "**************** Minuit Output ****************\n\n";

  // get fit results
  const auto& params = fmin.UserParameters().Params();
  const double pu = params[0];
  const double pv = params[1];
  const double pw = sqrt(1 + pu*pu + pv*pv);
  result.fPar[eU] = pu/pw;
  result.fPar[eV] = pv/pw;
  result.fPar[eRCore] = params[eRCore];
  result.fPar[eCTCore] = params[eCTCore];

  if (!fFixedCurvature && result.fPar[eRCore] > 100*kilometer) {
    OUT(eIntermediate)
      << "  curvature radius too large: " << result.fPar[eRCore]/kilometer << " km\n";
    return false;
  }

  if (!fFixedCurvature && result.fPar[eRCore] < 1*kilometer) {
    OUT(eIntermediate)
      << "  curvature radius too small: " << result.fPar[eRCore]/kilometer << " km\n";
    return false;
  }

  const auto nVar = pars.VariableParameters();

  // get covariance, propagate from (pu,pv) to (u,v)
  CovarianceMatrix cov(4);

  for (size_t i = 0; i < nVar; ++i)
    for (size_t j = 0; j < nVar; ++j)
      cov(m.ExtOfInt(i), m.ExtOfInt(j)) = fmin.UserCovariance()(i, j);

  const double pupv = pu*pv;
  const double pwi3 = 1/(pw*pw*pw);
  const double jacobian[4][4] = {
    { (1 + pv*pv)*pwi3, -pupv*pwi3, 0, 0 },
    { -pupv*pwi3, (1 + pu*pu)*pwi3, 0, 0 },
    { 0, 0, 1, 0 },
    { 0, 0, 0, 1 }
  };

  for (size_t i = 0; i < 4; ++i) {
    const auto& jac_i = jacobian[i];
    for (size_t j = i; j < 4; ++j) {
      const auto& jac_j = jacobian[j];
      auto& c_ij = result.fCov(i, j);
      c_ij = 0;
      for (size_t k = 0; k < 4; ++k) {
        const auto& jac_ik = jac_i[k];
        for (size_t m = 0; m < 4; ++m)
          c_ij += jac_ik * jac_j[m] * cov(k, m);
      }
    }
  }

  result.fChi2 = fmin.Fval();

  result.fNdof = 0;
  for (const auto& sf : data)
    if (sf.fSignal > 0)
      ++result.fNdof;
  result.fNdof -= nVar;

  return !HasNaN(result);
}


vector<int>
LDFFinder::CleanSilentStation(const LDFFitResult& ldfres,
                              const LDFFitConfig& config,
                              const utl::Vector& axis,
                              std::vector<StationFitData>& data)
  const
{
  vector<int> badStations;
  const utl::Point core(ldfres.fPar[1], ldfres.fPar[2], 0, gBaryCS);
  int nrCandidates = 0;
  for (const auto& station : data) {
    if (station.fSignal > 0)
      ++nrCandidates;
  }
  if (nrCandidates < fSilentsVars.fMinNumberOfStations)
    return badStations;

  const std::vector<double> shape(ldfres.fPar.begin()+3, ldfres.fPar.end());

  for (auto sIt = data.begin(); sIt != data.end(); ) {
    if (!sIt->fSignal) {
      const Vector coreToStation = sIt->fPos - core;
      const double rho = RPerp(axis, coreToStation);
      const double ldfValue = ldfres.fPar[0] * config.fLDF(rho, shape);
      if (rho < fSilentsVars.fMaxDistance &&
          ldfValue > fSilentsVars.fMaxLDFValue) {
        badStations.push_back(sIt->fId);
        sIt = data.erase(sIt);
        continue;
      }
    }
    ++sIt;
  }
  return badStations;
}


bool
LDFFinder::FitGlobal(CurvFitResult& cres, LDFFitResult& lres,
                     const LDFFitConfig& config,
                     const vector<StationFitData>& data)
  const
{
  using namespace ROOT::Minuit2;

  MnUserParameters pars;

  pars.Add("u", cres.fPar[0], 0.01);
  pars.Add("v", cres.fPar[1], 0.01);
  pars.Add("rCore", cres.fPar[2], 10*meter);
  pars.SetLowerLimit("rCore", 0);
  pars.Add("cTCore", cres.fPar[3], 10*meter);

  if (fFixedAxis) {
    pars.Fix("u");
    pars.Fix("v");
  }

  if (fFixedCurvature)
    pars.Fix("rCore");

  pars.Add("showerSize", lres.fPar[0], 1.0);
  pars.SetLowerLimit("showerSize", 0);
  pars.Add("coreX", lres.fPar[1], 200*meter);
  pars.Add("coreY", lres.fPar[2], 200*meter);
  for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i) {
    pars.Add(GetShapeLabel(i), lres.fPar[3+i], 0.1);
    if (config.fLDFShapeTreatment[i] != eFitted)
      pars.Fix(GetShapeLabel(i));
  }

  if (config.fFixCore) {
    pars.Fix("coreX");
    pars.Fix("coreY");
  }

  GlobalFitFunction likeFcn(config, data);

  MnMinimize m(likeFcn, pars, 2);

  FunctionMinimum fmin = m();

  if (fmin.IsAboveMaxEdm() || fmin.HasReachedCallLimit() || !fmin.IsValid()) {
    OUT(eIntermediate)
      << "\n**************** Minimize FAILED ****************\n"
      << fmin
      << "**************** Minimize FAILED ****************\n\n";
    return false;
  }

  {
    ostringstream msg;
    msg << "Found minimum after " << fmin.NFcn() << " evaluations";
    INFO(msg);
  }

  MnHesse hesse(2);
  hesse(m.Fcnbase(), fmin);

  if (!fmin.IsValid() || fmin.HesseFailed() || !fmin.HasAccurateCovar()) {
    OUT(eIntermediate)
      << "\n**************** Hesse FAILED ****************\n"
      << fmin
      << "**************** Hesse FAILED ****************\n\n";
    return false;
  }

  if (fInfoLevel >= eMinuit)
    OUT(eIntermediate)
      << "\n**************** Minuit Output ****************\n"
      << fmin
      << "**************** Minuit Output ****************\n\n";

  // get fit results
  const auto nPar = fmin.UserParameters().Params().size();
  for (size_t i = 0; i < nPar; ++i) {
    if (i < 4)
      cres.fPar[i] = fmin.UserParameters().Params()[i];
    else
      lres.fPar[i-4] = fmin.UserParameters().Params()[i];
  }

  const auto modeledShape = config.fLDF.ShapeModel(cres.GetAxis().GetCosTheta(gBaryCS), lres.fPar[0]);
  for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i)
    if (config.fLDFShapeTreatment[i] == eModeled)
      lres.fPar[3+i] = modeledShape[i];

  if (!fFixedCurvature && cres.fPar[eRCore] > 100*kilometer) {
    OUT(eIntermediate)
      << "  curvature radius too large.\n";
    return false;
  }

  const auto core = lres.GetCore();
  const double d = (core - fBarycenter).GetMag();
  if (d > fMaxBaryToCoreDistance) {
    OUT(eIntermediate)
      << "  Distance to core " << d/km << " km exceeds limit.\n";
    return false;
  }

  const auto nVar = pars.VariableParameters();

  // get covariance - not all measured covariances are stored for now
  cres.fCov = 0;
  lres.fCov = 0;

  for (size_t i = 0; i < nVar; ++i)
    for (size_t j = 0; j < nVar; ++j) {
      const auto iext = m.ExtOfInt(i);
      const auto jext = m.ExtOfInt(j);
      if (iext < 4 && jext < 4)
        cres.fCov(iext, jext) = fmin.UserCovariance()(i, j);
      if (iext >= 4 && jext >= 4)
        lres.fCov(iext-4, jext-4) = fmin.UserCovariance()(i, j);
    }

  {
    LDFLikelihoodFunction likeFcn(config, cres.GetAxis(), data);
    lres.fLogLikelihood = likeFcn(lres.fPar);

    LDFFitConfig chi2Config = config;
    chi2Config.fUseSilentStations = false;
    LDFChi2Function chi2Fcn(chi2Config, cres.GetAxis(), data);
    lres.fChi2 = chi2Fcn(lres.fPar);

    lres.fNdof = 0;
    for (const auto& station : data)
      if (station.fSignal > 0)
        ++lres.fNdof;

    for (size_t i = 4; i < nPar; ++i)
      if (!pars.Parameter(i).IsFixed())
        --lres.fNdof;

    CurvatureChi2Function curvChi2Fcn(lres.GetCore(), data);
    cres.fChi2 = curvChi2Fcn(cres.fPar);

    cres.fNdof = 0;
    for (const auto& station : data)
      if (station.fSignal > 0)
        ++cres.fNdof;

    for (size_t i = 0; i < 4; ++i)
      if (!pars.Parameter(i).IsFixed())
        --cres.fNdof;
  }

  const double reducedChi2 = (lres.fNdof > 0) ? lres.fChi2/lres.fNdof : 0;

  if (reducedChi2 > fMaxChi2) {
    OUT(eIntermediate)
      << "  ldf chi2 / Ndof = " << reducedChi2 << " exceeds limit " << fMaxChi2 << '\n';
    return false;
  }

  const utl::Point coreErr(lres.fCov.Std(eCoreX), lres.fCov.Std(eCoreY), 0, gBaryCS);
  OUT(eIntermediate)
    << "  curv Chi2 = " << cres.fChi2 << " / " << cres.fNdof << "\n"
       "  axis = (" << String::Make(cres.GetAxis(), gBaryCS, 1, ", ") << ") (bary)\n"
       "  RCore = " << cres.fPar[eRCore]/kilometer << " +/- " << cres.fCov.Std(eRCore)/kilometer << " km\n"
       "  CTCore = " << cres.fPar[eCTCore]/meter << " +/- " << cres.fCov.Std(eCTCore)/meter << " m\n"
       "  ldf chi2 = " << lres.fChi2 << " / " << lres.fNdof << "\n"
       "  " << GetShowerSizeLabel(config) << " = " << lres.fPar[eShowerSize] << " +/- " << lres.fCov.Std(eShowerSize) << " VEM\n"
       "  core = (" << String::Make(core, gBaryCS, meter, ", ") << ") +/- (" << String::Make(coreErr, gBaryCS, meter, ", ") << ") m (bary)";
  for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i)
    OUT(eIntermediate) << "\n"
      "  " << GetShapeLabel(i) << " = " << lres.fPar[eShapeParameters + i] << " +/- " << lres.fCov.Std(eShapeParameters + i);
  OUT(eIntermediate) << '\n';

  return !HasNaN(cres) && !HasNaN(lres);
}