UnivFitterKG.cc

Go to the documentation of this file.

#include "UnivFitterKG.h"

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <limits>
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <fstream>
#include <vector>

#include "Math/PdfFuncMathCore.h"
#include "TF1.h"
#include "Math/WrappedTF1.h"
#include "Math/RootFinder.h"
#include <TMath.h>

#include <boost/format.hpp>

#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnHesse.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnStrategy.h"

#include <tls/UnivTimeKGLogNormal.h>
#include <tls/UnivTimeKGGamma.h>

#include <utl/PhysicalFunctions.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/TimeInterval.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/PhysicalFunctions.h>
#include <utl/AugerUnits.h>
#include <utl/CovarianceMatrix.h>
#include <utl/NumericalErrorPropagator.h>
#include <utl/GeometryUtilities.h>

//#include <evt/Event.h>
//#include <sevt/SEvent.h>

using namespace std;
using namespace utl;
using namespace ROOT::Minuit2;

using boost::format;
using FitterUtil::sc_mu;
using utl::NormalizeAngleMinusPiPi;


namespace UnivFitterKG {

  const double gpercm2 = utl::gram / utl::cm2;


  inline
  double
  OldSmoothPoissonFactor(const double sigmaFactor)
  {
    const double z = (sigmaFactor - 1.05) / 0.05;
    const double t = 1 / (1 + std::exp(z));
    return t + (1 - t) / Sqr(sigmaFactor);
  }


  double
  RelativeTrackLength(const double cosTheta)
  {
    const double kRadius = 1.8 * meter;
    const double kHeight = 1.2 * meter;
    const double sinTheta = sqrt(1 - pow(cosTheta, 2.));
    return 1 / (cosTheta + 2 * kHeight * sinTheta / (kPi * kRadius));
  }


  double
  PlaneFrontTime(const CoordinateSystemPtr showerCS, const Point& corePosition, const Point& position)
  {
    return (corePosition.GetZ(showerCS) - position.GetZ(showerCS)) / kSpeedOfLight;
  }


  // calculates energy for proton, QGSJet 2.3, this is needed to get unbiased predictions with universality
  double
  MCEnergyCalibration(const double showerSize, const double theta, bool isInfill)
  {
    if (!isInfill) {
      const double calPars[2] = { 2.51594443, 1.03997067 };  // uncs: 0.02302776  0.00182409
      const double cicPars[3] = { 1.07546532, 2.10875764, -0.37034705 };  // uncs: 0.01036417  0.01998492  0.10950651

      const double cossqref = pow(cos(theta / radian), 2) - pow(cos(38 * degree), 2);
      const double s38 = showerSize / (1 + cossqref * (cicPars[0] + cossqref * (cicPars[1] + cossqref * cicPars[2])));
      return calPars[0] * 1e17 * pow(s38, calPars[1]);
    }
    return 0;
  }


  double
  ConvertToMCEnergyScale(const double sdenergy)
  {
    // MC calibration for p, QGSJet-II.03
    const double calPars_p[2] = { 2.51594443, 1.03997067 };
    // data calibration (ICRC 2013) - needs to be changed later
    const double calPars_data[2] = { 1.89541, 1.02465 };

    const double s38 = pow(sdenergy / (calPars_data[0] * 1e17), 1 / calPars_data[1]);

    return calPars_p[0] * 1e17 * pow(s38, calPars_p[1]);
  }


  WCDFitFunction::WCDFitFunction(const Fitter& theEf) : ef(theEf) { }


  double
  WCDFitFunction::operator()(const std::vector<double>& pars)
  const
  {
    // exception from Offline units:
    // grammage and energy are rescaled such that all parameters
    // are roughly on the same order of magnitude, see setupMinuit()

    // check sanity of parameters
    for (unsigned int i = 0; i < pars.size(); ++i) {
      if (std::isnan(pars[i]))
        return -ef.fFailLikelihood;
    }

    const Point core = Point(pars[0], pars[1], pars[2], ef.fBaryCS);
    const CoordinateSystemPtr coreCS = AugerCoordinateSystem::Create(core, ReferenceEllipsoid::eWGS84, ef.fSiteCS);
    const double relCoreTime = pars[3];
    const TimeStamp coreTime = ef.initPars.coretime + relCoreTime;

    const double theta = pars[4];
    const double phi = pars[5];
    const Vector axis(1, theta, phi, ef.fBaryCS, Vector::kSpherical);
    const CoordinateSystemPtr showerPlaneCS(
      CoordinateSystem::RotationY(theta, CoordinateSystem::RotationZ(phi, coreCS))
    );

    const double xmax = pars[6] * ef.fGrammageScale;
    double Nmu = pars[7];
    double totalEnergy = pars[8] * ef.fEnergyScale;

    if (ef.fUseNmuModel)
      Nmu = FitterUtil::getMeanParametrizedNmuData(log10(totalEnergy), theta, xmax / gpercm2);

    const double mmu_offset = pars[9];

    double x0 = pars[10] * ef.fGrammageScale;

    // TODO remove this? Remove also in Offline interface?
    if (ef.fRescaleEnergy) {
      using namespace utl::InvisibleEnergy;
      const double mf = ModelFactor(eQGSJET01, eProton);
      const double* const p = FitParameters(eProton);

      TF1 func("f", "[3]*([0]-[1]*(x/1e18)**(-[2])) - x/[4]");  // formula from PhysicalFunctions.cc
      func.SetParameter(0, p[0]);
      func.SetParameter(1, p[1]);
      func.SetParameter(2, p[2]);
      func.SetParameter(3, mf);
      func.SetParameter(4, totalEnergy);

      ROOT::Math::RootFinder rf;
      ROOT::Math::WrappedTF1 wf1(func);
      rf.SetFunction(wf1, 0.5 * totalEnergy, 1.5 * totalEnergy);
      rf.Solve();
      double calEnergy = rf.Root();

      double dummycostheta=0.7071067; //since this uses montecarlo, dummycostheta is not used,
      totalEnergy = calEnergy * (1 + Nmu * (Factor(calEnergy, eQGSJET01, eProton,dummycostheta) - 1));
    }

    const double logTotalEnergy = log10(totalEnergy);

    const double X0FromXmax = FitterUtil::getX0(ef.fX0XmaxModelVersion, xmax, logTotalEnergy);
    if (ef.x0Mode == eCoupled)
      x0 = X0FromXmax;
    const double hfi = ef.fAtmosphere->SlantDepthToHeight(x0 / gpercm2, theta) * cm;

    if (ef.fVerbosityLevel > 1) {
      printf("x %.4f, y %.4f, z %.4f, relct %.3f, ", pars[0], pars[1], pars[2], pars[3]);
      printf("theta %.3f°, phi %.3f°, ", theta / degree, NormalizeAngleMinusPiPi(phi) / degree);
      printf("xmax %.2f, x0 %.2f, ", xmax / gpercm2, x0 / gpercm2);
      printf("Nmu %.3f, loge %.3f, ", Nmu, logTotalEnergy);
      printf("\n");
    }

    const double sigmaFactor = wcd::SignalUncertaintyFactor(wcd::eGAP2014_035, cos(theta));
    const double poissonFactor = wcd::PoissonFactor(sigmaFactor);

    double TotalLikelihood = 0;
    ef.fLDFLikelihood = 0;
    ef.fShapeLikelihood = 0;
    ef.fStartTimeLikelihood = 0;

    for (std::vector<StationFitData>::const_iterator sIt = ef.fitData.begin(); sIt != ef.fitData.end(); ++sIt) {

      if (!(sIt->useForLDFFit || sIt->useForShapeFit || sIt->useForStartTimeFit))
        continue;

      sIt->LDFLikelihood = 0;
      sIt->ShapeLikelihood = 0;
      sIt->StartTimeLikelihood = 0;
      sIt->planeFrontTime = PlaneFrontTime(showerPlaneCS, core, sIt->position);  // pf time referred to core position
      const double stationRho = sIt->position.GetRho(showerPlaneCS);
      const double stationPsi = sIt->position.GetPhi(showerPlaneCS);
      sIt->height = ef.fSiteCSHeight + sIt->position.GetZ(ef.fSiteCS);

      sIt->DX_diffusive = ef.fAtmosphere->Get_DX_DL_i(stationRho / cm, stationPsi, xmax / gpercm2, theta, sIt->height / cm, false, true) * gpercm2;
      sIt->DX_rectilinear = ef.fAtmosphere->Get_DX_DL_i(stationRho / cm, stationPsi, xmax / gpercm2, theta, sIt->height / cm, false, false) * gpercm2;
      sIt->distToFirstInteraction = FitterUtil::getDistToFirstInteraction(hfi, sIt->height, theta, stationRho, stationPsi);

      sIt->firstParticleTimes.clear();
      if (ef.fTimeModelVersion == 1) {
        const double fpt_global = FitterUtil::getFirstParticleArrivalTime(stationRho, sIt->distToFirstInteraction);
        for (int icomp = 0; icomp < ef.fNSigComps; ++icomp)
          sIt->firstParticleTimes.push_back(fpt_global);
      } else if (ef.fTimeModelVersion > 1) {
        for (int icomp = 0; icomp < ef.fNSigComps; ++icomp) {

          const double D_TO = ef.fUnivParamTime->GetD_TO(xmax / gpercm2, theta, logTotalEnergy, icomp);
          const double h1 = ef.fAtmosphere->SlantDepthToHeight(xmax / gpercm2, theta) + D_TO * 1e5 * cos(theta);

          sIt->firstParticleTimes.push_back(ef.fAtmosphere->GetTimeCF_h(stationRho / cm, stationPsi, h1, theta, sIt->height / cm));
        }
      }

      const Vector stationToCore = core - sIt->position;
      const Vector stationToCenter = stationToCore + sIt->distToFirstInteraction * axis;
      const double rStation = stationToCenter.GetMag();
      const double cosThetaStation = stationToCenter.GetZ(ef.fBaryCS) / rStation;  // should be stationCS to be really exact

      // need to add offsets for MC events (depends on primary, interaction model)
      double mmoff = 0;
      if (ef.TimeModelOffsetFixed)
        if (ef.isMCEvent)
          mmoff = FitterUtil::GetMeanMuonOffsetMCCalibrated(stationRho / cm, logTotalEnergy, theta);
        else
          mmoff = FitterUtil::GetMeanMuonOffsetDataCalibrated(stationRho / cm, logTotalEnergy, theta);
      else
        mmoff = FitterUtil::GetMeanMuonOffset(stationRho / cm, mmu_offset);

      for (int icomp = 0; icomp < ef.fNSigComps; ++icomp) {

        // set mean muon offset for muonic comp. and muon decay products
        if (ef.fUseTimeModelOffset && (icomp == 0 || icomp == 2))
          ef.timeModels[icomp]->setParameterOffsets(mmoff, 0);

        ef.timeModels[icomp]->setShapeParameters(sIt->DX_diffusive / gpercm2, stationRho, stationPsi, theta, logTotalEnergy);
      }

      double predTotalSignal = 0;
      sIt->compSignals.clear();
      for (int icomp = 0; icomp < ef.fNSigComps; ++icomp) {

        const double sig =
          ef.fUnivParam->GetSignal(stationRho / cm, xmax / gpercm2, logTotalEnergy,
                                   theta, stationPsi, ef.fDensity / (gram / cm3),
                                   sIt->height / cm, Nmu, icomp, ef.atmModel);
        sIt->compSignals.push_back(sig);
        predTotalSignal += sig;
      }
      predTotalSignal *= ef.fUPFudgeFactor;

      const double totalSignalUncertainty = poissonFactor * sigmaFactor * std::sqrt(predTotalSignal);
      sIt->signalUncertainty = totalSignalUncertainty;

      const double stationTime = (sIt->startTime - coreTime).GetInterval();

      // start time correction as it is also applied in UnivRec.cc
      // implementation in UnivParamTime.cc, details in GAP2013_072
      // if (stationRho < 1000 * meter)
      //   stationTime -= ef.fUnivParamTime->tStartCorrection(stationRho / cm, logTotalEnergy, theta, false) * nanosecond;

      const double stationTimePF = stationTime - sIt->planeFrontTime - ef.fStartTimeBias;

      double StationLikelihood = 0;
      const double fem = sIt->compSignals[1] / predTotalSignal;

      if (sIt->useForShapeFit && !ef.fOnlySignalLikelihood) {

        double StationShapeLikelihood = 0;

        if (ef.fQuantileShapeFit) {

          const double qs[2] = { 0.1, 0.5 };
          const double qts[2] = { sIt->t10, sIt->t50 };

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

            const double pfac = ef.fUnivParamTime->GetPoissonFactor(fem, qs[qi]);
            const double nps = predTotalSignal * pfac;
            const double ih = nps * qs[qi];

            double pdfval = 0;
            double cdfval = 0;

            for (int icomp = 0; icomp < ef.fNSigComps; ++icomp) {
              pdfval += ef.timeModels[icomp]->pdf(stationTimePF + qts[qi] - sIt->firstParticleTimes[icomp]);
              cdfval += ef.timeModels[icomp]->cdf(stationTimePF + qts[qi] - sIt->firstParticleTimes[icomp]);
            }

            pdfval /= 4;
            cdfval /= 4;

            const double norm = 1 / TMath::Beta(nps - ih, ih + 1);
            const double LLq = log(norm * pow(1 - cdfval, nps - ih - 1) * pow(cdfval, ih) * pdfval);

            StationShapeLikelihood += (std::isnan(LLq) || std::isinf(LLq)) ? 200 : LLq;

          }

        } else {

          for (unsigned int curBin = 0; curBin < sIt->trace.size(); ++curBin) {
            const double curBinContent = sIt->trace[curBin];

            if (sIt->isSaturated) {
              if (curBinContent < ef.fMinBinContentForSatShapeFit ||
                  curBinContent > ef.fMaxBinContentForSatShapeFit ||
                  curBinContent > 0.95 * sIt->Smaxbin)
                continue;
            } else {
              if (curBinContent < ef.fMinBinContentForShapeFit ||
                  curBinContent > ef.fMaxBinContentForShapeFit)
                continue;
            }

            if (ef.doPartialShapeFit &&
                (int(curBin) < sIt->firstFitBin || int(curBin) > sIt->lastFitBin))
              continue;

            const double binPFTime = stationTimePF + (curBin + 1 - sIt->startBin) * sIt->binwidth;

            if (ef.fEarlyTraceFit && (binPFTime - stationTimePF > sIt->t50 + sIt->binwidth)) {
              if (ef.fVerbosityLevel > 3)
                cout << "Bin " << curBin << " is skipped because it is later than t50 " << sIt->t50 << endl;
              continue;
            }

            double predictedBinContent = 0;
            double predictedBinComponentContent = 0;
            const double totalBinSigma = poissonFactor * sigmaFactor * std::sqrt(curBinContent);
            double StationShapeBinLikelihood = 0;

            // calculation from shape model
            for (int icomp = 0; icomp < ef.fNSigComps; ++icomp) {
              predictedBinComponentContent =
                (ef.timeModels[icomp]->cdf(binPFTime - sIt->firstParticleTimes[icomp] + sIt->binwidth / 2) -
                 ef.timeModels[icomp]->cdf(binPFTime - sIt->firstParticleTimes[icomp] - sIt->binwidth / 2)) *
                sIt->compSignals[icomp];
              predictedBinContent += predictedBinComponentContent;
            }

            StationShapeBinLikelihood += LogarithmOfNormalPDF(curBinContent, predictedBinContent, totalBinSigma);
            if (ef.fVerbosityLevel > 3) {
              cout << "id=" << sIt->stationId << " bin=" << curBin
                   << " S=" << curBinContent << " pred=" << predictedBinContent
                   << " LL=" << -StationShapeBinLikelihood << endl;
            }

            StationShapeLikelihood += StationShapeBinLikelihood;

          }
        }

        StationLikelihood += StationShapeLikelihood;
        ef.fShapeLikelihood += StationShapeLikelihood;
        sIt->ShapeLikelihood = StationShapeLikelihood;

      }

      if (sIt->useForStartTimeFit && !ef.fOnlySignalLikelihood) {
        const double muonsInStation = sIt->compSignals[sc_mu] / RelativeTrackLength(cosThetaStation);

        const double startTimeProb =
          ef.smearFirstParticlePdf ?
          ef.timeModels[sc_mu]->firstParticlePdfSmeared(stationTimePF - sIt->firstParticleTimes[sc_mu], muonsInStation) :
          ef.timeModels[sc_mu]->firstParticlePdf(stationTimePF - sIt->firstParticleTimes[sc_mu], muonsInStation);
        const double StationStartTimeLikelihood = log(startTimeProb);

        if (ef.fVerbosityLevel > 2) {
          cout
              << " id " << sIt->stationId
              << " npart " << muonsInStation
              << " fpt " << sIt->firstParticleTimes[sc_mu] << " stpf " << stationTimePF << " stprob " << startTimeProb
              << " pft " << format(" %.10e ") % sIt->planeFrontTime
              << " stl " << StationStartTimeLikelihood
              << endl;
        }
        StationLikelihood += StationStartTimeLikelihood;
        ef.fStartTimeLikelihood += StationStartTimeLikelihood;
        sIt->StartTimeLikelihood = StationStartTimeLikelihood;
      }

      if (sIt->useForLDFFit || ((sIt->useForLDFFit || sIt->useForShapeFit) && ef.fOnlySignalLikelihood && !sIt->isSaturated)) {

        double StationLDFLikelihood = 0;
        if (ef.fVerbosityLevel > 2) {
          cout
              << " id " << sIt->stationId
              << " signal " << sIt->totalSignal
              << " pred " << predTotalSignal
              << " sigma " << totalSignalUncertainty
              << endl;
        }
        StationLDFLikelihood +=
          TMath::Log(TMath::Poisson(poissonFactor * sIt->totalSignal, poissonFactor * predTotalSignal));

        if (!(std::isnan(StationLDFLikelihood) || std::isinf(StationLDFLikelihood))) {
          StationLikelihood += StationLDFLikelihood;
          ef.fLDFLikelihood += StationLDFLikelihood;
          sIt->LDFLikelihood = StationLDFLikelihood;
        }
      }

      // constraining the total predicted signal to saturation recovered measured signal
      if (ef.addRecovSignalConstraint && ef.doLDFFit && sIt->isSaturated && sIt->recoveredSignal > 0) {
        const double constrLL = LogarithmOfNormalPDF(sIt->recoveredSignal, predTotalSignal, sIt->recoveredSignalUncertainty);
        StationLikelihood += constrLL;
      }

      if (sIt->useForMuonTimeFit) {
        double StationMuonTimeLikelihood = 0;

        for (unsigned int k = 0; k < sIt->muonTimes.size(); k++) {
          const double muonTimePF = sIt->muonTimes[k] - sIt->planeFrontTime;

          const double MuonTimeLikelihood = log(ef.timeModels[sc_mu]->pdf(muonTimePF + 50 * nanosecond - sIt->firstParticleTimes[sc_mu]));
          StationMuonTimeLikelihood += MuonTimeLikelihood;
          if (ef.fVerbosityLevel > 2) {
            cout << " id " << sIt->stationId
                 << " timepf " << muonTimePF + 50 * nanosecond
                 << " fpt " << sIt->firstParticleTimes[sc_mu]
                 << " muon LL " << MuonTimeLikelihood
                 << endl;
          }
        }
        StationLikelihood += StationMuonTimeLikelihood;
      }

      if (ef.addMuonLDFConstraint && ef.doLDFFit && (sIt->useForLDFFit || sIt->useForShapeFit)) {
        const double predMuonSignal = sIt->compSignals[sc_mu];
        //const double muonConstraint = LogarithmOfNormalPDF(sIt->muonSignal, predMuonSignal, sIt->muonSignalError);
        const double muonSignal = sIt->muonTimes.size() * RelativeTrackLength(cosThetaStation);
        const double muonConstraint = LogarithmOfNormalPDF(muonSignal, predMuonSignal, sqrt(muonSignal));
        StationLikelihood += muonConstraint;
      }

      TotalLikelihood += StationLikelihood;
    }

    if (ef.x0Mode != eFixed && ef.fX0Constraint.sigma > 0)
      TotalLikelihood += LogarithmOfNormalPDF(x0, ef.initPars.x0, ef.fX0Constraint.sigma);

    if (!ef.coreFixed) {
      if (ef.fCoreConstraint.sigma > 0) {
        const double distToInit = (core - ef.initPars.core).GetMag();
        TotalLikelihood += LogarithmOfNormalPDF(distToInit, 0, ef.fCoreConstraint.sigma);
      }
      if (ef.fCoreTimeConstraint.sigma > 0)
        TotalLikelihood += LogarithmOfNormalPDF(relCoreTime, 0, ef.fCoreTimeConstraint.sigma);
    }

    if (!ef.energyFixed && ef.fEnergyConstraint.sigma > 0)
      TotalLikelihood += LogarithmOfNormalPDF(totalEnergy, ef.fEnergyConstraint.mean, ef.fEnergyConstraint.sigma);

    if (!ef.xmaxFixed && ef.fXmaxConstraint.sigma > 0)
      TotalLikelihood += LogarithmOfNormalPDF(xmax, ef.fXmaxConstraint.mean, ef.fXmaxConstraint.sigma);

    if (!ef.axisFixed) {
      if (ef.fThetaConstraint.sigma > 0)
        TotalLikelihood += LogarithmOfNormalPDF(theta, ef.fThetaConstraint.mean, ef.fThetaConstraint.sigma);

      if (ef.fPhiConstraint.sigma > 0)
        TotalLikelihood += LogarithmOfNormalPDF(phi, ef.fPhiConstraint.mean, ef.fPhiConstraint.sigma);
    }

    if (std::isnan(TotalLikelihood) || std::isinf(TotalLikelihood))
      TotalLikelihood = ef.fFailLikelihood;

    if (ef.fVerbosityLevel > 1)
      printf("ll %.3f\n", TotalLikelihood);

    return -TotalLikelihood;
  }


  double
  AMIGAFitFunction::operator()(const std::vector<double>& /*pars*/)
  const
  {
    // TODO implement
    return 0;
  }


  double
  ASCIIFitFunction::operator()(const std::vector<double>& /*pars*/)
  const
  {
    // TODO implement
    return 0;
  }


  FitFunction::FitFunction(const Fitter& theEf): ef(theEf) { }


  double
  FitFunction::Up()
  const
  {
    return 0.5; // 1.0 for chi2, 0.5 for negative log likelihood
  }


  double
  FitFunction::operator()(const std::vector<double>& pars)
  const
  {
    WCDFitFunction& wcd = *(ef.fWCDFunc);
    //AMIGAFitFunction& amiga = *(ef.fAMIGAFunc);
    //ASCIIFitFunction& ascii = *(ef.fASCIIFunc);
    return wcd(pars);
  }


  // Note the "Minuit units" here (for grammage and energy, see above)
  minScan1D
  Fitter::scanMinimum1D(const int parIndex, const double val, const double delta, const int nSteps)
  {
    const double minScanValue = val - delta;
    const double maxScanValue = val + delta;
    const double scanStep = (maxScanValue - minScanValue) / nSteps;

    vector<double> params = fFittedMinuitParams;

    minScan1D scanVals;

    // suppress output temporarily
    const int verblevel = fVerbosityLevel;
    if (fVerbosityLevel <= 2)
      fVerbosityLevel = 0;

    for (double scanValue = minScanValue; scanValue < maxScanValue; scanValue += scanStep) {
      scanVals.x.push_back(scanValue);
      params[parIndex] = scanValue;
      const double funcval = (*fFunc)(params);
      scanVals.f.push_back(funcval);
      scanVals.fLDF.push_back(-fLDFLikelihood);
      scanVals.fShape.push_back(-fShapeLikelihood);
      scanVals.fStartTime.push_back(-fStartTimeLikelihood);
    }
    if (fVerbosityLevel <= 2)
      fVerbosityLevel = verblevel;
    return scanVals;
  }


  // Note the "Minuit units" here (for grammage and energy, see above)
  minScan2D
  Fitter::scanMinimum2D(const int parIndexx, const int parIndexy, const double valx, const double valy,
                        const double deltaX, const double deltaY, const int nSteps)
  {
    const double minScanValuex = valx - deltaX;
    const double maxScanValuex = valx + deltaX;
    const double minScanValuey = valy - deltaY;
    const double maxScanValuey = valy + deltaY;
    const double scanStepx = (maxScanValuex - minScanValuex) / nSteps;
    const double scanStepy = (maxScanValuey - minScanValuey) / nSteps;

    vector<double> params = fFittedMinuitParams;

    minScan2D scanVals;

    // suppress output temporarily
    const int verblevel = fVerbosityLevel;
    if (fVerbosityLevel < 3)
      fVerbosityLevel = 0;

    for (int i = 0; i < nSteps; ++i) {
      scanVals.x.push_back(minScanValuex + i * scanStepx);
      scanVals.y.push_back(minScanValuey + i * scanStepy);
    }

    for (int iy = 0; iy < nSteps; ++iy) {
      for (int ix = 0; ix < nSteps; ++ix) {
        params[parIndexx] = scanVals.x[ix];
        params[parIndexy] = scanVals.y[iy];
        scanVals.f.push_back(-(*fFunc)(params));
      }
    }
    if (fVerbosityLevel < 3)
      fVerbosityLevel = verblevel;
    return scanVals;
  }


  void
  Fitter::writeMinScan()
  {
    const int nSteps = 1000;

    //minScan1D scan = scanMinimum1D(0, fittedPars.core.GetX(fBaryCS), 50, nSteps);
    //minScan1D scan = scanMinimum1D(1, fittedPars.core.GetY(fBaryCS), 50, nSteps);
    //minScan1D scan = scanMinimum1D(4, fittedPars.axis.GetTheta(fBaryCS), .2*degree, nSteps);
    //minScan1D scan = scanMinimum1D(5, fittedPars.axis.GetPhi(fBaryCS), .05*degree, nSteps);
    minScan1D scan = scanMinimum1D(6, fittedPars.xmax / gpercm2, 300, nSteps);

    std::ofstream out("scan_out.txt");
    for (unsigned int i = 0; i < scan.x.size(); ++i) {
      out << format("%.15e ") % scan.x[i]
          << format("%.15e ") % scan.f[i]
          << format("%.15e ") % scan.fLDF[i]
          << format("%.15e ") % scan.fShape[i]
          << format("%.15e ") % scan.fStartTime[i]
          << endl;
    }
  }


  Fitter::Fitter()
  {
    fUnivParam.reset(new UnivParamNS::UnivParam(0));  // 0 for WCD
    fUnivParamTime.reset(new UnivParamTimeNS::UnivParamTime(0));  // 0 for WCD

    fFunc.reset(new FitFunction(*this));
    fWCDFunc.reset(new WCDFitFunction(*this));
    //fAMIGAFunc.reset(new AMIGAFitFunction(*this));
    //fASCIIFunc.reset(new ASCIIFitFunction(*this));

    fSiteCSHeight = 1400 * m;

    CoordinateSystemPtr ecef = CoordinateSystem::GetRootCoordinateSystem();
    ReferenceEllipsoid wgs84 = ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

    UTMPoint malargueOrigin(6060000 * m, 440000 * m, 0, 19, 'H', wgs84);
    CoordinateSystemPtr malargueCS = AugerBaseCoordinateSystem::Create(malargueOrigin, ecef);

    UTMPoint pampaAmarillaOrigin(-35.25 * deg, -69.25 * deg, fSiteCSHeight, wgs84);
    CoordinateSystemPtr pampaAmarillaCS = AugerCoordinateSystem::Create(pampaAmarillaOrigin, malargueCS);

    fSiteCS = pampaAmarillaCS;

    fNSigComps = 4;
    fUPFudgeFactor = 1.0;  // new plots confirm correct prediction of signals when using local muon number
    fStartTimeBias = 14.5 * nanosecond; // bias in start time that needs to be corrected in order to unbias MC

    fMinBinsAboveShapeFitLimit = 5;
    fMinBinContentForShapeFit = 0.7;
    fMaxBinContentForShapeFit = std::numeric_limits<double>::max();
    fMinBinContentForSatShapeFit = 0.7;
    fMaxBinContentForSatShapeFit = std::numeric_limits<double>::max();

    fLowerRadiusCutLDF = 100 * meter;
    fUpperRadiusCutLDF = 2500 * meter;
    fLowerRadiusCutShape = 100 * meter;
    fUpperRadiusCutShape = 2200 * meter;
    fLowerRadiusCutTime = 100 * meter;
    fUpperRadiusCutTime = 2200 * meter;

    fXmaxLowerLimit = 200 * gpercm2;
    fXmaxUpperLimit = 2000 * gpercm2;
    fAxisFitLimit = 10 * degree;
    fCoreFitLimit = 0 * meter;

    fMigradTol = 0.1;
    // the likelihood function returns -1*total likelihood
    // so this value becomes a huge positive value for minuit
    fFailLikelihood = -1e8;
    fMigradStrategy = 1;
    fEnergyScale = 1e18 * eV;
    fGrammageScale = gram / cm2;
    fVerbosityLevel = 0;

    fPrintResults = true;

    // apply missing energy correction (off for MC showers)
    // Not necessary since Observer v9r3 due to new missing energy correction (GAP 2013-026)
    fRescaleEnergy = false;

    // for old time model (part of the bias stems from problems with start time fit)
    fApplyXmaxBiasCorrection = false;
    fUseMCEnergyScale = false;
    fEarlyTraceFit = false;
    fQuantileShapeFit = false;
    fIterativeFit = false;
    fCalibMode = false;
    fOnlySignalLikelihood = false;
    fUseTimeModelOffset = true;
    fUseNmuModel = false;

    isMCEvent = false;

    fX0XmaxModelVersion = 2;  // see FitterUtil.cc

    doLDFFit = true; // fit to the total signal
    doShapeFit = true; // bin-by-bin fit to VEM trace
    doPartialShapeFit = false; // fit to VEM trace between predefined quantiles (e.g. to avoid long tails)
    doStartTimeFit = false; // fit start time to pdf of first particle, leads to an additional bias with the current impl.
    doSaturatedStartTimeFit = false;
    addMuonLDFConstraint = false;
    doMuonTimeFit = false;
    smearFirstParticlePdf = true;
    addRecovSignalConstraint = false;

    coreFixed = false;
    axisFixed = false;
    xmaxFixed = false;
    NmuFixed = false;
    energyFixed = true;
    TimeModelOffsetFixed = true;

    x0Mode = eCoupled;

    reset();
  }


  Fitter::~Fitter()
  {
    clearTimeModel();
  }


  void
  Fitter::reset()
  {
    fitData.clear();

    mcPars.valid = false;
    sdPars.valid = false;
    fdPars.valid = false;

    eventHasSaturation = false;
    atmModel = 0;
    atmType = "monthly";
    atmUseGDAS = false;

    fNStationsLDF = 0;
    fNStationsShape = 0;
    fNStationsStartTime = 0;
  }


  void
  Fitter::initTimeModel(unsigned int version)
  {
    fTimeModelVersion = version;
    if (version == 1) {
      for (int icomp = 0; icomp < fNSigComps; ++icomp)
        timeModels.push_back(new UnivTimeKG::LogNormalTimeModel(icomp));
    } else if (version == 2) {
      for (int icomp = 0; icomp < fNSigComps; ++icomp)
        timeModels.push_back(new UnivTimeKG::GammaTimeModel(icomp));
    } else {
      cout << "ERROR: invalid time model version" << endl;
      throw;
    }
  }


  void
  Fitter::clearTimeModel()
  {
    for (int icomp = 0; icomp < fNSigComps; ++icomp)
      delete timeModels[icomp];
    timeModels.clear();
  }


  void
  Fitter::setRecType(unsigned int type)
  {
    // all options are set to useful values in the constructor of this class
    // they are replaced depending on the type of reconstruction

    switch (type) {

      // simulated SD: energy fixed to MC energy
      case 1:
        initMode = eMC;
        break;

      // real SD events: energy fixed to SD energy
      case 2:
        initMode = eSD;
        break;

      // real SD/FD events: energy fixed to FD energy
      case 3:
        initMode = eFD;
        break;

      // MC events: calibration mode
      case 4:
        initMode = eMC;

        coreFixed = false;
        axisFixed = true;
        xmaxFixed = true;
        NmuFixed = false;
        energyFixed = true;
        TimeModelOffsetFixed = false;

        fCalibMode = true;
        fIterativeFit = false;

        x0Mode = eCoupled;

        break;

      // real golden hybrid events: parameters from FD reconstruction, Nmu fitted, no time fit
      case 5:
        initMode = eFD;

        coreFixed = false;
        axisFixed = true;
        xmaxFixed = true;
        NmuFixed = false;
        energyFixed = true;
        TimeModelOffsetFixed = true;

        fIterativeFit = false;

        x0Mode = eCoupled;

        break;

      // real golden hybrid events: parameters from FD reconstruction, Nmu fitted, fit TimeModelOffset
      case 6:
        initMode = eFD;

        coreFixed = false;
        axisFixed = true;
        xmaxFixed = true;
        NmuFixed = false;
        energyFixed = true;
        TimeModelOffsetFixed = false;

        fCalibMode = true;
        fIterativeFit = false;

        x0Mode = eCoupled;

        break;

    }
  }


  void
  Fitter::addStationFitData(const StationFitData& sdata)
  {
    fitData.push_back(sdata);
  }


  StationFitData&
  Fitter::getStationFitData(const int stationId)
  {
    for (std::vector<StationFitData>::iterator sIt = fitData.begin(); sIt != fitData.end(); ++sIt) {
      if (sIt->stationId == stationId)
        return *sIt;
    }
    throw;
  }


  bool
  Fitter::isShapeFitCandidate(const int stationId)
  {
    StationFitData& s(getStationFitData(stationId));
    return s.useForShapeFit;
  }


  bool
  Fitter::isLDFFitCandidate(const int stationId)
  {
    StationFitData& s(getStationFitData(stationId));
    return s.useForLDFFit;
  }


  bool
  Fitter::isStartTimeFitCandidate(const int stationId)
  {
    StationFitData& s(getStationFitData(stationId));
    return s.useForStartTimeFit;
  }


  // return true if the station was selected for the fit
  bool
  Fitter::isCandidateStation(const int stationId)
  {
    StationFitData& s(getStationFitData(stationId));
    return s.useForShapeFit || s.useForLDFFit || s.useForStartTimeFit;
  }


  void
  Fitter::setupMinuit(MnUserParameters& upar)
  {
    upar.Add("coreX", initPars.core.GetX(fBaryCS), 50.*meter);
    upar.Add("coreY", initPars.core.GetY(fBaryCS), 50.*meter);
    upar.Add("coreZ", initPars.core.GetZ(fBaryCS), 50.*meter);

    upar.Add("relCoreTime", 0, 25 * nanosecond);
    upar.SetLimits("relCoreTime", -1e6 * nanosecond, 1e6 * nanosecond);

    if (fCoreFitLimit > 0) {
      upar.SetLimits("coreX", upar.Value("coreX") - fCoreFitLimit, upar.Value("coreX") + fCoreFitLimit);
      upar.SetLimits("coreY", upar.Value("coreY") - fCoreFitLimit, upar.Value("coreY") + fCoreFitLimit);
    }

    if (coreFixed) {
      upar.Fix("coreX");
      upar.Fix("coreY");
      upar.Fix("relCoreTime");
    }
    upar.Fix("coreZ");

    upar.Add("theta", initPars.axis.GetTheta(fBaryCS), 1 * degree);
    upar.Add("phi", initPars.axis.GetPhi(fBaryCS), 1 * degree);
    if (fAxisFitLimit > 0) {
      double minTheta = upar.Value("theta") - fAxisFitLimit;
      if (minTheta < 0)
        minTheta = 0;
      upar.SetLimits("theta", minTheta, upar.Value("theta") + fAxisFitLimit);
      upar.SetLimits("phi", upar.Value("phi") - fAxisFitLimit, upar.Value("phi") + fAxisFitLimit);
    }
    if (axisFixed) {
      upar.Fix("theta");
      upar.Fix("phi");
    }

    // The numeric values of eV and g/cm^2 are very large.
    // Energy and grammage are rescaled for Minuit such that
    // the values are roughly on the same order of magnitude.
    upar.Add("xmax", initPars.xmax / fGrammageScale, initPars.xmaxUnc / fGrammageScale);
    upar.SetLimits("xmax", fXmaxLowerLimit / fGrammageScale, fXmaxUpperLimit / fGrammageScale);
    if (xmaxFixed)
      upar.Fix("xmax");

    upar.Add("Nmu", initPars.Nmu, initPars.NmuUnc);
    upar.SetLimits("Nmu", 0., 4.);
    if (NmuFixed)
      upar.Fix("Nmu");

    upar.Add("energy", initPars.energy / fEnergyScale, .2 * initPars.energy / fEnergyScale);
    upar.SetLimits("energy", 0.5 * initPars.energy / fEnergyScale, 1.5 * initPars.energy / fEnergyScale);
    if (energyFixed)
      upar.Fix("energy");

    upar.Add("TimeModelOffset", initPars.TimeModelOffset, 1);
    upar.SetLimits("TimeModelOffset", -3., 3.);
    if (TimeModelOffsetFixed)
      upar.Fix("TimeModelOffset");

    upar.Add("x0", initPars.x0 / fGrammageScale, 10.);
    upar.SetLimits("x0", 0., 500);
    if (x0Mode == eFixed || x0Mode == eCoupled)
      upar.Fix("x0");
  }


  void
  Fitter::init()
  {
    fAtmosphere.reset(new AtmosphereNS::Atmosphere(atmUseGDAS, atmType));
    fAtmosphere->SetCurrentAtmosphere(atmModel);
    fUnivCalibConstants.reset(new UnivCalibConstantsNS::UnivCalibConstants(fCalibOpt, 0, fRecMixture));

    if (isMCEvent && mcPars.valid) {
      mcPars.coreCS = AugerCoordinateSystem::Create(mcPars.core, ReferenceEllipsoid::eWGS84, fSiteCS);
      mcPars.showerCS = CoordinateSystem::RotationY(mcPars.axis.GetTheta(fSiteCS),
                        CoordinateSystem::RotationZ(mcPars.axis.GetPhi(fSiteCS), mcPars.coreCS));

      const CoordinateSystemPtr rotatedMCCoreCS = CoordinateSystem::RotationZ(mcPars.axis.GetPhi(fSiteCS), mcPars.coreCS);

      for (std::vector<StationFitData>::iterator sIt = fitData.begin(); sIt != fitData.end(); ++sIt) {
        // This is needed only if UnviFitterKG is used outside of the Offline framework
        // since the position is only known in (r,psi) referred to the shower axis.
        // If used within Offline, the position (utl::Point) can be set directly.
        if (sIt->isDense) {
          sIt->position = Point(sIt->denseRho * cos(sIt->densePsi) / cos(mcPars.axis.GetTheta(fSiteCS)),
                                sIt->denseRho * sin(sIt->densePsi), 0, rotatedMCCoreCS);
          //cout << sIt->stationId << " " << sIt->denseRho << " " << sIt->densePsi/degree << " "
          //<< sIt->position.GetRho(mcPars.showerCS) << " " << NormalizeAngle(sIt->position.GetPhi(mcPars.showerCS))/degree << endl;
        }
      }
    }

    if (sdPars.valid) {
      sdPars.coreCS = AugerCoordinateSystem::Create(sdPars.core, ReferenceEllipsoid::eWGS84, fSiteCS);
      sdPars.showerCS = CoordinateSystem::RotationY(sdPars.axis.GetTheta(fSiteCS),
                        CoordinateSystem::RotationZ(sdPars.axis.GetPhi(fSiteCS), sdPars.coreCS));
      sdPars.xmax = FitterUtil::getMeanParametrizedXmax(log10(sdPars.energy));
      sdPars.xmaxUnc = 40 * gpercm2;
      sdPars.Nmu = 1.9;
    }

    if (fdPars.valid) {
      fdPars.coreCS = AugerCoordinateSystem::Create(fdPars.core, ReferenceEllipsoid::eWGS84, fSiteCS);
      fdPars.showerCS = CoordinateSystem::RotationY(fdPars.axis.GetTheta(fSiteCS),
                        CoordinateSystem::RotationZ(fdPars.axis.GetPhi(fSiteCS), fdPars.coreCS));
      fdPars.Nmu = 1.9;
    }

    if (initMode == eMC && mcPars.valid)
      initPars = mcPars;
    else if (initMode == eSD && sdPars.valid)
      initPars = sdPars;
    else if (initMode == eFD && fdPars.valid)
      initPars = fdPars;
    else {
      cout << "ERROR: invalid initMode" << endl;
      throw;
    }

    initPars.x0 = FitterUtil::getX0(fX0XmaxModelVersion, initPars.xmax, log10(initPars.energy));
    initPars.x0Unc = 20 * gpercm2;
    initPars.NmuUnc = 0.3;
    initPars.TimeModelOffset = 0.;

    const Point siteOrigin(0, 0, 0, fSiteCS);
    Vector barySum(0, 0, 0, fSiteCS);
    double weightSum = 0;

    for (std::vector<StationFitData>::const_iterator sIt = fitData.begin(); sIt != fitData.end(); ++sIt) {
      const double weight = sqrt(sIt->totalSignal);   // GAP 2007-035
      weightSum += weight;
      barySum += weight * (sIt->position - siteOrigin);
    }
    barySum /= weightSum;
    fBarycenter = siteOrigin + barySum;
    fBaryCS = AugerCoordinateSystem::Create(fBarycenter, ReferenceEllipsoid::eWGS84, fSiteCS);

    double hottestStationSignal = 0;
    Point hottestStationPosition;
    for (std::vector<StationFitData>::iterator sIt = fitData.begin(); sIt != fitData.end(); ++sIt) {
      if (sIt->totalSignal > hottestStationSignal) {
        hottestStationPosition = sIt->position;
        hottestStationSignal = sIt->totalSignal;
        if (sIt->isSaturated)
          cout << "hottest station " << sIt->stationId << " is saturated" << endl;
      }
    }

    fHottestStationCS = AugerCoordinateSystem::Create(hottestStationPosition, ReferenceEllipsoid::eWGS84, fBaryCS);

    // recalculate energy with respect to p, QGSJet2.3 reference
    if (fUseMCEnergyScale && (initMode == eSD || initMode == eFD)) {
      if (initMode == eFD && !sdPars.valid) {
        cout << "ERROR: Energy rescaling requires Sd parameters at the moment!" << endl;
        throw;
      }
      // initPars.energy = MCEnergyCalibration(sdPars.showersize, sdPars.axis.GetTheta(fBaryCS), false);
      initPars.energy = ConvertToMCEnergyScale(initPars.energy);
    }

    if (fVerbosityLevel > 0) {
      cout << "\n"
           << "station selection:\n"
           << "-------------------------\n";
    }

    for (std::vector<StationFitData>::iterator sIt = fitData.begin(); sIt != fitData.end(); ++sIt) {

      const double stationRho = sIt->position.GetRho(initPars.showerCS);

      int nBinsAboveShapeFitLimit = 0;

      double Smaxbin(0);

      if (sIt->trace.size() > 0) {
        for (unsigned int curBin = 0; curBin < sIt->trace.size(); ++curBin) {
          if (sIt->trace[curBin] >= fMinBinContentForShapeFit)
            ++nBinsAboveShapeFitLimit;
          if (sIt->trace[curBin] > Smaxbin)
            Smaxbin = sIt->trace[curBin];
        }

        for (int curBin = sIt->startBin; curBin <= sIt->stopBin; ++curBin) {
          sIt->traceSum += sIt->trace[curBin];  // different for saturated stations from total signal
        }                                       // when unsaturated traces are used

        double traceCumulative = 0;
        for (unsigned int curBin = 0; curBin < sIt->trace.size(); ++curBin) {
          traceCumulative += sIt->trace[curBin];
          if (!sIt->firstFitBin && traceCumulative >= 0.05 * sIt->traceSum)
            sIt->firstFitBin = curBin;
          if (!sIt->lastFitBin && traceCumulative >= 0.8 * sIt->traceSum)
            sIt->lastFitBin = curBin - 1;
        }
      }

      sIt->Smaxbin = Smaxbin;

      eventHasSaturation |= sIt->isSaturated;

      // sIt->useForShapeFit = doShapeFit && (!sIt->isSaturated || sIt->useUnsaturatedTrace)
      sIt->useForShapeFit =
        doShapeFit &&
        nBinsAboveShapeFitLimit >= fMinBinsAboveShapeFitLimit &&
        stationRho > fLowerRadiusCutShape &&
        stationRho < fUpperRadiusCutShape;

      sIt->useForLDFFit =
        doLDFFit &&
        !sIt->useForShapeFit &&
        !sIt->isSaturated &&
        stationRho > fLowerRadiusCutLDF &&
        stationRho < fUpperRadiusCutLDF;

      sIt->useForStartTimeFit =
        doStartTimeFit &&
        !sIt->useForShapeFit &&
        stationRho > fLowerRadiusCutTime &&
        stationRho < fUpperRadiusCutTime &&
        (!sIt->isSaturated || doSaturatedStartTimeFit);

      if (sIt->useForLDFFit)
        ++fNStationsLDF;
      if (sIt->useForShapeFit)
        ++fNStationsShape;
      if (sIt->useForStartTimeFit)
        ++fNStationsStartTime;

      if (doMuonTimeFit)
        sIt->useForMuonTimeFit = true;

      if (fVerbosityLevel > 0) {
        cout << "id " << sIt->stationId
             << "  r " << showpoint << stationRho // same length for all numbers (shows zero at the end)
             << "  S " << showpoint << sIt->totalSignal << "  "
             << (sIt->useForLDFFit ? "ldf " : "    ")
             << (sIt->useForShapeFit ? "shape " : "      ")
             << (sIt->useForStartTimeFit ? "time " : "     ")
             << (sIt->isSaturated ? "saturated " : "          ")
             << (sIt->useUnsaturatedTrace ? "(use simulated unsaturated trace)" : "                                 ");

        if (sIt->useForMuonTimeFit)
          cout << sIt->muonTimes.size() << " muons";

        if (doPartialShapeFit)
          cout << " fit bins " << sIt->firstFitBin << " " << sIt->lastFitBin;
        cout << endl;
      }
    }
  }


  bool
  Fitter::run()
  {
    init();

    if (fVerbosityLevel > 0)
      cout << "-------------------------\n";

    bool fitOk = true;
    bool hesseOk = true;
    MnUserParameters upar;
    setupMinuit(upar);

    if (upar.VariableParameters() > 0) {
      const int maxSteps = 100000;

      MnMigrad migrad(*fFunc, upar, MnStrategy(fMigradStrategy));
      FunctionMinimum* min = 0;

      // can also add and remove constraints during iterative fitting
      // just set respective sigma values of the ef.xxxConstraint values

      if (fCalibMode) {

        migrad.Fix("xmax");
        migrad.Fix("theta");
        migrad.Fix("phi");
        migrad.Fix("energy");
        migrad.Fix("TimeModelOffset");
        migrad.Fix("relCoreTime");

        migrad.Release("coreX");
        migrad.Release("coreY");
        migrad.Release("Nmu");

        fOnlySignalLikelihood = true;

        min = new FunctionMinimum(migrad(maxSteps, fMigradTol));

        migrad.Fix("Nmu");
        migrad.Fix("coreX");
        migrad.Fix("coreY");
        // migrad.Fix("relCoreTime");

        migrad.Release("relCoreTime");
        migrad.Release("TimeModelOffset");

        TimeModelOffsetFixed = false;
        fOnlySignalLikelihood = false;

        *min = migrad(maxSteps, fMigradTol);

      } else {

        // warning: fixing of certain parameters is overwritten here
        if (fIterativeFit) {

          const bool dbgprint(true);

          coreFixed = false;
          axisFixed = false;
          xmaxFixed = false;
          NmuFixed = false;
          energyFixed = false;
          TimeModelOffsetFixed = true;

          fUseTimeModelOffset = true;

          const double lgE0 = log10(migrad.Value("energy") * fEnergyScale);
          //const double theta0 = migrad.Value("theta");
          //const double phi0 = migrad.Value("phi");

          // stage 1: core + energy, energy constrained
          // only signal model in use

          // constraining
          fEnergyConstraint.mean = pow(10, lgE0);
          fEnergyConstraint.sigma = 0.2 * pow(10, lgE0);

          fOnlySignalLikelihood = true;

          // need to add new parametrizations for MC events
          // should be general enough to incorporate various interaction models and primary types

          if (!isMCEvent) {

            const double lgE = log10(migrad.Value("energy") * fEnergyScale);
            const double theta = migrad.Value("theta");

            const double xmax = FitterUtil::getMeanParametrizedXmaxData(lgE);
            const double nmu = FitterUtil::getMeanParametrizedNmuData(lgE, theta, xmax);

            migrad.SetValue("xmax", xmax);
            migrad.SetValue("Nmu", nmu);
          }

          if (dbgprint) {
            cout << "Energy before stage 1: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax before stage 1: " << migrad.Value("xmax") << endl;
            cout << "Nmu before stage 1: " << migrad.Value("Nmu") << endl;
            cout << "Zenith before stage 1: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime before stage 1: " << migrad.Value("relCoreTime") << endl;
          }

          migrad.Fix("TimeModelOffset");

          migrad.Fix("xmax");
          migrad.Fix("Nmu");
          migrad.Fix("theta");
          migrad.Fix("phi");
          migrad.Fix("relCoreTime");

          if (!eventHasSaturation)
            migrad.Release("energy");
          else
            migrad.Fix("energy");

          migrad.Release("coreX");
          migrad.Release("coreY");

          min = new FunctionMinimum(migrad(maxSteps, fMigradTol));

          if (dbgprint) {
            cout << "Energy after stage 1: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax after stage 1: " << migrad.Value("xmax") << endl;
            cout << "Nmu after stage 1: " << migrad.Value("Nmu") << endl;
            cout << "Zenith after stage 1: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime after stage 1: " << migrad.Value("relCoreTime") << endl;
          }

          // stage 2: core + energy (updated model values after fitting energy), tighter energy constraint
          // only signal model is used

          fEnergyConstraint.mean = migrad.Value("energy") * fEnergyScale;
          fEnergyConstraint.sigma = 0.2 * fEnergyConstraint.mean;

          if (!isMCEvent) {
            const double lgE = log10(migrad.Value("energy") * fEnergyScale);
            const double theta = migrad.Value("theta");

            const double xmax = FitterUtil::getMeanParametrizedXmaxData(lgE);
            const double nmu = FitterUtil::getMeanParametrizedNmuData(lgE, theta, xmax);

            migrad.SetValue("xmax", xmax);
            migrad.SetValue("Nmu", nmu);
          }

          if (dbgprint) {
            cout << "Energy before stage 2: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax before stage 2: " << migrad.Value("xmax") << endl;
            cout << "Nmu before stage 2: " << migrad.Value("Nmu") << endl;
            cout << "Zenith before stage 2: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime after stage 2: " << migrad.Value("relCoreTime") << endl;
          }

          migrad.Fix("xmax");
          migrad.Fix("Nmu");
          migrad.Fix("theta");
          migrad.Fix("phi");
          migrad.Fix("relCoreTime");

          if (!eventHasSaturation)
            migrad.Release("energy");
          else
            migrad.Fix("energy");

          migrad.Release("coreX");
          migrad.Release("coreY");

          *min = migrad(maxSteps, fMigradTol);

          if (dbgprint) {
            cout << "Energy after stage 2: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax after stage 2: " << migrad.Value("xmax") << endl;
            cout << "Nmu after stage 2: " << migrad.Value("Nmu") << endl;
            cout << "Zenith after stage 2: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime after stage 2: " << migrad.Value("relCoreTime") << endl;
          }

          fOnlySignalLikelihood = false;

          // stage 3a: relative core time

          fEnergyConstraint.sigma = 0;

          migrad.Fix("energy");
          migrad.Fix("coreX");
          migrad.Fix("coreY");
          migrad.Fix("Nmu");
          migrad.Fix("theta");
          migrad.Fix("phi");

          migrad.Release("relCoreTime");

          *min = migrad(maxSteps, fMigradTol);

          // stage3b: xmax + relative core time (+ Nmu model)

          // if (!isMCEvent)
          //   fUseNmuModel = true;

          migrad.Release("relCoreTime");
          migrad.Release("xmax");

          *min = migrad(maxSteps, fMigradTol);

          if (dbgprint) {
            cout << "Energy after stage 3: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax after stage 3: " << migrad.Value("xmax") << endl;
            cout << "Nmu after stage 3: " << migrad.Value("Nmu") << endl;
            cout << "Zenith after stage 3: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime after stage 3: " << migrad.Value("relCoreTime") << endl;
          }

          // stage 4: xmax, geometry, relative core time, geometry constrained

          fThetaConstraint.mean = migrad.Value("theta");
          fThetaConstraint.sigma = 1 * degree;

          fPhiConstraint.mean = migrad.Value("phi");
          fPhiConstraint.sigma = 1 * degree;

          // fXmaxConstraint.sigma = 0;

          /*if (!isMCEvent) {

            const double lgE = log10(migrad.Value("energy") * fEnergyScale);
            const double theta = migrad.Value("theta");
            const double xmax = migrad.Value("xmax");

            // double nmu(0);

            // if ((xmax > 700) & (xmax < 900)) // validity range of the model
            //   nmu = FitterUtil::getMeanParametrizedNmuData(lgE, theta, xmax);
            // else
            //   nmu = FitterUtil::getMeanParametrizedNmuData(lgE, theta);

            // migrad.SetValue("Nmu", nmu);
          }*/

          if (dbgprint) {
            cout << "Energy before stage 4: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax before stage 4: " << migrad.Value("xmax") << endl;
            cout << "Nmu before stage 4: " << migrad.Value("Nmu") << endl;
            cout << "Zenith before stage 4: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime before stage 4: " << migrad.Value("relCoreTime") << endl;
          }

          migrad.Fix("Nmu");
          migrad.Fix("energy");
          migrad.Fix("coreX");
          migrad.Fix("coreY");

          migrad.Release("theta");
          migrad.Release("phi");

          *min = migrad(maxSteps, fMigradTol);

          if (!isMCEvent && fUseNmuModel) {
            fUseNmuModel = false;

            const double lgE = log10(migrad.Value("energy") * fEnergyScale);
            const double theta = migrad.Value("theta");
            const double xmax = migrad.Value("xmax");

            const double nmu = FitterUtil::getMeanParametrizedNmuData(lgE, theta, xmax);

            migrad.SetValue("Nmu", nmu);
          }

          if (dbgprint) {
            cout << "Energy after stage 4: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax after stage 4: " << migrad.Value("xmax") << endl;
            cout << "Nmu after stage 4: " << migrad.Value("Nmu") << endl;
            cout << "Zenith after stage 4: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime after stage 4: " << migrad.Value("relCoreTime") << endl;
          }

          // catching misreconstructions by fixing the energy in previous stages
          // if (migrad.Value("xmax") > 1200) {

          //   cout << "Warning: Xmax is larger than 1100, trying to recover with fixed energy !!!" << endl;

          //   migrad.SetValue("energy", pow(10, lgE0) / fEnergyScale);

          //   const double xmax = FitterUtil::getMeanParametrizedXmaxData(lgE0);
          //   const double nmu = FitterUtil::getMeanParametrizedNmuData(lgE0, theta0, xmax);

          //   migrad.SetValue("xmax", xmax);
          //   migrad.SetValue("Nmu", nmu);
          //   migrad.SetValue("theta", theta0);
          //   migrad.SetValue("phi", phi0);

          //   fThetaConstraint.sigma = 0;
          //   fPhiConstraint.sigma = 0;

          //   // refitting core with fixed energy

          //   migrad.Fix("theta");
          //   migrad.Fix("phi");
          //   migrad.Fix("xmax");
          //   migrad.Fix("relCoreTime");

          //   migrad.Release("coreX");
          //   migrad.Release("coreY");

          //   *min = migrad(maxSteps, fMigradTol);

          //   // refitting xmax and timing

          //   migrad.Fix("coreX");
          //   migrad.Fix("coreY");

          //   // migrad.Release("theta");
          //   // migrad.Release("phi");
          //   migrad.Release("xmax");
          //   migrad.Release("relCoreTime");

          //   *min = migrad(maxSteps, fMigradTol);

          //   // refitting xmax, timing and geometry

          //   fThetaConstraint.mean = theta0;
          //   fThetaConstraint.sigma = 0.5 / 180. * M_PI;

          //   fPhiConstraint.mean = phi0;
          //   fPhiConstraint.sigma = 0.5 / 180. * M_PI;

          //   migrad.Release("theta");
          //   migrad.Release("phi");
          //   migrad.Release("xmax");
          //   migrad.Release("relCoreTime");

          //   *min = migrad(maxSteps, fMigradTol);

          // }

          // stage 5: re-fit of Nmu

          fThetaConstraint.sigma = 0;
          fPhiConstraint.sigma = 0;
          fEnergyConstraint.sigma = 0;

          fOnlySignalLikelihood = true;

          migrad.Fix("energy");
          migrad.Fix("coreX");
          migrad.Fix("coreY");
          migrad.Fix("theta");
          migrad.Fix("phi");
          migrad.Fix("relCoreTime");
          migrad.Fix("xmax");

          migrad.Release("Nmu");

          *min = migrad(maxSteps, fMigradTol);

          if (dbgprint) {
            cout << "Energy after stage 5: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax after stage 5: " << migrad.Value("xmax") << endl;
            cout << "Nmu after stage 5: " << migrad.Value("Nmu") << endl;
            cout << "Zenith after stage 5: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime after stage 5: " << migrad.Value("relCoreTime") << endl;
          }

          fOnlySignalLikelihood = false;

          // stage 6: xmax, geometry, geometry constrained (updated Nmu value)

          fThetaConstraint.mean = migrad.Value("theta");
          fThetaConstraint.sigma = 1 * degree;

          fPhiConstraint.mean = migrad.Value("phi");
          fPhiConstraint.sigma = 1 * degree;

          migrad.Fix("Nmu");

          migrad.Release("theta");
          migrad.Release("phi");
          migrad.Release("xmax");

          *min = migrad(maxSteps, fMigradTol);

          if (dbgprint) {
            cout << "Energy after stage 6: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax after stage 6: " << migrad.Value("xmax") << endl;
            cout << "Nmu after stage 6: " << migrad.Value("Nmu") << endl;
            cout << "Zenith after stage 6: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime after stage 6: " << migrad.Value("relCoreTime") << endl;
          }

          // stage 7: re-fit energy as previous last step (using signal and time models)

          fThetaConstraint.sigma = 0;
          fPhiConstraint.sigma = 0;

          fEnergyConstraint.mean = migrad.Value("energy") * fEnergyScale;
          fEnergyConstraint.sigma = 0.1 * fEnergyConstraint.mean;

          migrad.Fix("Nmu");
          migrad.Fix("theta");
          migrad.Fix("phi");
          migrad.Fix("xmax");
          migrad.Release("energy");

          *min = migrad(maxSteps, fMigradTol);

          if (dbgprint) {
            cout << "Energy after stage 7: " << log10(migrad.Value("energy") * fEnergyScale) << endl;
            cout << "Xmax after stage 7: " << migrad.Value("xmax") << endl;
            cout << "Nmu after stage 7: " << migrad.Value("Nmu") << endl;
            cout << "Zenith after stage 7: " << migrad.Value("theta") << endl;
            cout << "RelCoreTime after stage 7: " << migrad.Value("relCoreTime") << endl;
          }

          // stage 8: re-fit Nmu as last step (using signal model)

          // fEnergyConstraint.sigma = 0;

          // fOnlySignalLikelihood = true;

          // migrad.Fix("energy");
          // migrad.Fix("theta");
          // migrad.Fix("phi");
          // migrad.Fix("xmax");
          // migrad.Release("Nmu");

          // *min = migrad(maxSteps, fMigradTol);

          // fOnlySignalLikelihood = false;

          // if (dbgprint)
          //   cout << "Xmax after stage 8: " << migrad.Value("xmax") << endl;

        } else
          min = new FunctionMinimum(migrad(maxSteps, fMigradTol));
      }

      if (fVerbosityLevel > 1)
        cout << *min;

      // The likelihood near the minimum is slightly non-parabolic in the x/y direction in rare cases.
      // In this case, Minuit overestimates the EDM and pretends it didn't converge
      // although the fitted value is fine (checked by hand many cases).
      // After fixing the core and running again, the change in the remaining free parameters is very small.
      if (!isFitOk(*min) && !fIterativeFit) {
        cout << "\n"
             "No convergence according to Minuit. Fixing the core position and running again."
             << endl;
        migrad.Fix("coreX");
        migrad.Fix("coreY");
        *min = migrad(maxSteps, fMigradTol);
      }

      if (fVerbosityLevel > 1)
        cout << *min;
      fitOk = isFitOk(*min);

      if (fVerbosityLevel > 1)
        cout << "\n\n### calculating hesse matrix\n" << endl;
      const int verblevel = fVerbosityLevel;
      fVerbosityLevel = 0;
      MnHesse hesse(fMigradStrategy);
      hesse(migrad.Fcnbase(), *min, maxSteps);
      fVerbosityLevel = verblevel;
      hesseOk = !min->HesseFailed();

      const MnUserParameterState& us = min->UserState();
      fFittedMinuitParams = us.Params();
      if (fVerbosityLevel > 1)
        cout << min->UserCovariance() << endl;

      fittedPars.axis = Vector(1., us.Value("theta"), us.Value("phi"), fBaryCS, Vector::kSpherical);
      fittedPars.axisUncTheta = us.Error("theta");
      fittedPars.axisUncPhi = us.Error("phi");

      fittedPars.core = Point(us.Value("coreX"), us.Value("coreY"), us.Value("coreZ"), fBaryCS);
      fittedPars.coreUncX = us.Error("coreX");
      fittedPars.coreUncY = us.Error("coreY");
      fittedPars.coretime = initPars.coretime + us.Value("relCoreTime");
      fittedPars.coreCS = AugerCoordinateSystem::Create(fittedPars.core, ReferenceEllipsoid::eWGS84, fSiteCS);
      fittedPars.showerCS = CoordinateSystem::RotationY(fittedPars.axis.GetTheta(fSiteCS),
                            CoordinateSystem::RotationZ(fittedPars.axis.GetPhi(fSiteCS), fittedPars.coreCS));

      fittedPars.Nmu = us.Value("Nmu");
      fittedPars.NmuUnc = us.Error("Nmu");
      fittedPars.energy = us.Value("energy") * fEnergyScale;
      fittedPars.energyUnc = us.Error("energy") * fEnergyScale;

      fittedPars.xmax = us.Value("xmax") * fGrammageScale;
      fittedPars.xmaxUnc = us.Error("xmax") * fGrammageScale;

      fittedPars.TimeModelOffset = us.Value("TimeModelOffset");
      fittedPars.TimeModelOffsetUnc = us.Error("TimeModelOffset");

      if (fApplyXmaxBiasCorrection)
        fittedPars.xmax -= FitterUtil::getMeanXmaxBias(log10(fittedPars.energy),
                           fittedPars.axis.GetTheta(fBaryCS), eventHasSaturation);

      fittedPars.x0 = us.Value("x0") * fGrammageScale;
      fittedPars.x0Unc = us.Error("x0") * fGrammageScale;

      if (x0Mode == eCoupled) {
        fittedPars.x0 = FitterUtil::getX0(fX0XmaxModelVersion, fittedPars.xmax, log10(fittedPars.energy));
        fittedPars.x0Unc = 0;
      }

      fittedPars.hfi = fAtmosphere->SlantDepthToHeight(fittedPars.x0 / gpercm2, fittedPars.axis.GetTheta(fBaryCS)) * cm;

      // call fcn at the end to set final likelihood values (total + per station) etc...
      (*fFunc)(us.Params());

      delete min;

    } else {
      // this is done so that fittedPars can be used from outside in the same way
      // even if all parameters are fixed (e.g. to check the prediction of the MC values)
      fittedPars = initPars;
      (*fFunc)(upar.Params());
    }

    if (fPrintResults)
      printResults();

    return fitOk && hesseOk;
  }


  bool
  Fitter::isFitOk(const FunctionMinimum& min)
  {
    const MnUserParameterState& us = min.UserState();
    if (!min.IsValid() || !us.IsValid() || min.IsAboveMaxEdm() || min.HasReachedCallLimit())
      return false;

    const vector<double>& pars = us.Params();
    for (unsigned int i = 0; i < pars.size(); ++i) {
      if (std::isnan(pars[i]))
        return false;
    }

    return true;
  }


  void
  Fitter::printResults()
  {
    const double eoutsc = 1e18;

    cout << endl;

    cout << format("E/eV") << endl;
    if (isMCEvent)
      cout << format("             MC   %.2f * 10^%i") % (mcPars.energy / eoutsc) % log10(eoutsc) << endl;
    cout << format("             init %.2f * 10^%i") % (initPars.energy / eoutsc) % log10(eoutsc) << endl;
    if (!energyFixed)
      cout << format("             fit  (%.2f ± %.2f) * 10^%i") % (fittedPars.energy / eoutsc) % (fittedPars.energyUnc / eoutsc) % log10(eoutsc) << endl;
    cout << (fEnergyConstraint.sigma > 0 ? "(constrained)" : "");
    cout << endl;


    cout << "Nmu" << endl;
    if (isMCEvent)
      cout << format("             MC   %.2f") % (mcPars.Nmu) << endl;
    cout << format("             init %.2f") % (initPars.Nmu) << endl;
    if (!NmuFixed) {
      cout << format("             fit  %.2f ± %.3f") % (fittedPars.Nmu) % (fittedPars.NmuUnc);
      if (isMCEvent)
        cout << format(" (%+.2f to MC)") % (fittedPars.Nmu - mcPars.Nmu);
      cout << endl;
    }
    cout << endl;


    cout << "xmax/gcm^-2" << endl;
    if (isMCEvent)
      cout << format("             MC   %.2f") % (mcPars.xmax / gpercm2) << endl;
    cout << format("             init %.2f") % (initPars.xmax / gpercm2) << endl;
    if (!xmaxFixed) {
      cout << format("             fit  %.2f ± %.2f") % (fittedPars.xmax / gpercm2) % (fittedPars.xmaxUnc / gpercm2);
      if (isMCEvent)
        cout << format(" (%+.2f to MC)") % ((fittedPars.xmax - mcPars.xmax) / gpercm2);
      cout << endl;
    }
    cout << endl;


    if (x0Mode == eFree) {
      cout << "x0/gcm^-2" << endl;
      if (isMCEvent)
        cout << format("             MC   %.2f") % (mcPars.x0 / gpercm2) << endl;
      cout << format("             init %.2f") % (initPars.x0 / gpercm2) << endl;
      cout << format("             fit  %.2f ± %.2f") % (fittedPars.x0 / gpercm2) % (fittedPars.x0Unc / gpercm2);
      if (isMCEvent)
        cout << format(" (%+.2f to MC)") % ((fittedPars.x0 - mcPars.x0) / gpercm2);
      cout << endl;
    }


    cout << "core (site) " << endl;
    if (isMCEvent)
      cout << format("             MC   %.2f, %.2f") % mcPars.core.GetX(fSiteCS) % mcPars.core.GetY(fSiteCS) << endl;
    cout << format("             init %.2f, %.2f") % initPars.core.GetX(fSiteCS) % initPars.core.GetY(fSiteCS) << endl;
    if (!coreFixed) {
      cout << format("             fit  %.2f ± %.2f, %.2f ± %.2f") % fittedPars.core.GetX(fSiteCS) % fittedPars.coreUncX
           % fittedPars.core.GetY(fSiteCS) % fittedPars.coreUncY;
      if (isMCEvent)
        cout << format(" (%.2fm to MC)") % (fittedPars.core - mcPars.core).GetMag();
      cout << endl;
    }
    cout << endl;


    cout << "core time " << endl;
    if (isMCEvent)
      cout << format("             MC   s %i ns %10.3f") % mcPars.coretime.GetGPSSecond()
           % mcPars.coretime.GetGPSNanoSecond() << endl;
    cout << format("             init s %i ns %10.3f") % initPars.coretime.GetGPSSecond()
         % initPars.coretime.GetGPSNanoSecond() << endl;
    if (!coreFixed) {
      cout << format("             fit  s %i ns %10.3f") % fittedPars.coretime.GetGPSSecond()
           % fittedPars.coretime.GetGPSNanoSecond();
      if (isMCEvent)
        cout << format(" (%.3f ns to MC)") % (fittedPars.coretime - initPars.coretime).GetInterval();
      cout << endl;
    }
    cout << endl;


    cout << "axis (site)" << endl;
    if (isMCEvent)
      cout << format("             MC   theta %.2f°, phi %.2f°") % (mcPars.axis.GetTheta(fSiteCS) / degree)
           % (NormalizeAngleMinusPiPi(mcPars.axis.GetPhi(fSiteCS)) / degree) << endl;
    cout << format("             init theta %.2f°, phi %.2f°") % (initPars.axis.GetTheta(fSiteCS) / degree)
         % (NormalizeAngleMinusPiPi(initPars.axis.GetPhi(fSiteCS)) / degree) << endl;
    if (!axisFixed) {
      cout << format("             fit  theta %.2f ± %.2f, phi %.2f ± %.2f")
           % (fittedPars.axis.GetTheta(fSiteCS) / degree) % (fittedPars.axisUncTheta / degree)
           % (NormalizeAngleMinusPiPi(fittedPars.axis.GetPhi(fSiteCS)) / degree) % (fittedPars.axisUncPhi / degree);
      if (isMCEvent)
        cout << format(" (%.2f° to MC)") % (getAngle(fittedPars.axis, mcPars.axis) / degree);
      cout << endl;
    }

    cout << "TimeModelOffset (dM_mu)" << endl;
    if (!TimeModelOffsetFixed)
      cout << format("             fit  %.2f ± %.3f") % (fittedPars.TimeModelOffset) % (fittedPars.TimeModelOffsetUnc);
    else
      cout << "fixed to calibrated value (see code)!";
    cout << endl;

    if (fVerbosityLevel > 0) {
      cout << "number of stations: LDF " << fNStationsLDF
           << " shape " << fNStationsShape << " time " << fNStationsStartTime << endl
           << "likelihood: LDF " << fLDFLikelihood
           << " shape " << fShapeLikelihood << " time " << fStartTimeLikelihood << endl
           << endl
           << "likelihood contributions:" << endl
           <<  "id    ldf    shape  time" << endl;

      for (std::vector<StationFitData>::iterator sIt = fitData.begin(); sIt != fitData.end(); ++sIt) {
        if (!(sIt->useForLDFFit || sIt->useForShapeFit || sIt->useForStartTimeFit)) continue;
        cout << format("%i  %.3f  %.3f  %.3f") % sIt->stationId
             % sIt->LDFLikelihood % sIt->ShapeLikelihood
             % sIt->StartTimeLikelihood << endl;
      }
      cout << endl;
    }
  }

}