RdGlobalFit.cc

Go to the documentation of this file.

#include "RdGlobalFit.h"
#include <utl/ErrorLogger.h>

#include <cmath>

#include <fwk/CentralConfig.h>
#include <revt/StationRecData.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/BasicVector.h>
#include <utl/GeometryUtilities.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/StringCompare.h>

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

#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>
#include <revt/StationRRecDataQuantities.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>

// include root stuff
#include "TGraph2DErrors.h"
#include "TF2.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TFitResultPtr.h"
#include "TRandom3.h"

#include <cstddef>
#include <string>
#include <cstdlib>
#include <iostream>

#include <sstream>
#include <fstream>
#include <iomanip>
#include <stdlib.h>
#include <sys/stat.h>
#include <math.h>
#include <algorithm>
#include <vector>

using namespace revt;
using namespace evt;
using namespace fwk;
using namespace std;
using namespace utl;
using namespace det;

//#include "DetectorGeometry.h"
#include <fwk/CentralConfig.h>
#include <det/VManager.h>
#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <utl/Branch.h>
#include <utl/RadioGeometryUtilities.h>

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

#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/StationConstants.h>
#include <revt/StationConstants.h>

#include <utl/TraceAlgorithm.h>

#include <sstream>
#include <algorithm>
#include <TMinuit.h>
#include <TRandom.h>

// atmosphere
#include <atm/AerosolDB.h>
#include <atm/AerosolZone.h>
#include <atm/ProfileResult.h>
#include <atm/MonthlyAvgDBProfileModel.h>

using namespace std;
using namespace fwk;
using namespace det;
using namespace utl;
using namespace evt;
using namespace revt;

#define OUT(x) if ((x) <= fInfoLevel) cerr << "  "
#define NL "\n  "

namespace RdGlobalFit {

RdGlobalFit::RdGlobalFit() :
    fInfoLevel(0),
    fGlobalFitData()
{
}

RdGlobalFit::~RdGlobalFit()
{
}

VModule::ResultFlag RdGlobalFit::Init()

{
  // Initialize your module here. This method
  // is called once at the beginning of the run.
  // The eSuccess flag indicates the method ended
  // successfully.  For other possible return types,
  // see the VModule documentation.

  INFO("RdGlobalFit::Init()");

  // Read in the configurations of the xml file
  Branch topBranch;
  topBranch = CentralConfig::GetInstance()->GetTopBranch("RdGlobalFit");
  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
  //B.Z.: TODO If fuseOwnSettings
  // topBranch.GetChild("UsedDirection").GetData(fUsedDirection);
  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

  topBranch.GetChild("A_scale").GetData(fLDFConstsTable.A_scale);
  topBranch.GetChild("C3").GetData(fLDFConstsTable.C3);
  topBranch.GetChild("C4").GetData(fLDFConstsTable.C4);

  topBranch.GetChild("C1theta_0_10").GetData(fLDFConstsTable.C1theta_0_10);
  topBranch.GetChild("C2theta_0_10").GetData(fLDFConstsTable.C2theta_0_10);

  topBranch.GetChild("C1theta_10_20").GetData(fLDFConstsTable.C1theta_10_20);
  topBranch.GetChild("C2theta_10_20").GetData(fLDFConstsTable.C2theta_10_20);

  topBranch.GetChild("C1theta_20_30").GetData(fLDFConstsTable.C1theta_20_30);
  topBranch.GetChild("C2theta_20_30").GetData(fLDFConstsTable.C2theta_20_30);

  topBranch.GetChild("C1theta_30_40").GetData(fLDFConstsTable.C1theta_30_40);
  topBranch.GetChild("C2theta_30_40").GetData(fLDFConstsTable.C2theta_30_40);

  topBranch.GetChild("C1theta_40_50").GetData(fLDFConstsTable.C1theta_40_50);
  topBranch.GetChild("C2theta_40_50").GetData(fLDFConstsTable.C2theta_40_50);

  topBranch.GetChild("C1theta_50_55").GetData(fLDFConstsTable.C1theta_50_55);
  topBranch.GetChild("C2theta_50_55").GetData(fLDFConstsTable.C2theta_50_55);
  topBranch.GetChild("A_scale_50_55").GetData(fLDFConstsTable.A_scale_50_55);

  topBranch.GetChild("C1theta_55").GetData(fLDFConstsTable.C1theta_55);
  topBranch.GetChild("C2theta_55").GetData(fLDFConstsTable.C2theta_55);
  topBranch.GetChild("A_scale_55").GetData(fLDFConstsTable.A_scale_55);

  utl::Branch bFitConfig = topBranch.GetChild("FitConfig");
  bFitConfig.GetChild("fitArrivalTime").GetData(fFitConfig.fitArrivalTime);
  bFitConfig.GetChild("fitTwoDLDF").GetData(fFitConfig.fitTwoDLDF);
  bFitConfig.GetChild("fixCore").GetData(fFitConfig.fixCore);
  bFitConfig.GetChild("fitGammaAndSigmaPlusIndependently").GetData(
      fFitConfig.fitGammaAndSigmaPlusIndependently);
  bFitConfig.GetChild("fitShowerMaxAsRMax").GetData(fFitConfig.fitShowerMaxAsRMax);
  bFitConfig.GetChild("fitShowerMaxAsXMAX").GetData(fFitConfig.fitShowerMaxAsXMAX);
  bFitConfig.GetChild("fitShowerMaxAsDXMAX").GetData(fFitConfig.fitShowerMaxAsDXMax);
  bFitConfig.GetChild("SaveCoreContourData").GetData(fFitConfig.SaveCoreContourData);
  bFitConfig.GetChild("printMinuitScans").GetData(fFitConfig.printMinuitScans);
  bFitConfig.GetChild("useCS").GetData(fFitConfig.useCS);

  fOrigFitConfig = fFitConfig;
  utl::Branch bGlobalFitData = topBranch.GetChild("GlobalFitData");
  utl::Branch bAplus = bGlobalFitData.GetChild("Aplus");
  bAplus.GetChild("isFix").GetData(fGlobalFitData.Aplus.isFix);
  bAplus.GetChild("hasLimits").GetData(fGlobalFitData.Aplus.hasLimits);
  bAplus.GetChild("startValueFrom").GetData(fGlobalFitData.Aplus.startValueFrom);
  bAplus.GetChild("startValue").GetData(fGlobalFitData.Aplus.startValue);
  bAplus.GetChild("startError").GetData(fGlobalFitData.Aplus.startError);
  bAplus.GetChild("lowLimit").GetData(fGlobalFitData.Aplus.lowLimit);
  bAplus.GetChild("highLimit").GetData(fGlobalFitData.Aplus.highLimit);

  utl::Branch bcore_x = bGlobalFitData.GetChild("core_x");
  bcore_x.GetChild("isFix").GetData(fGlobalFitData.core_x.isFix);
  bcore_x.GetChild("hasLimits").GetData(fGlobalFitData.core_x.hasLimits);
  bcore_x.GetChild("startValueFrom").GetData(fGlobalFitData.core_x.startValueFrom);
  bcore_x.GetChild("startValue").GetData(fGlobalFitData.core_x.startValue);
  bcore_x.GetChild("startError").GetData(fGlobalFitData.core_x.startError);
  bcore_x.GetChild("lowLimit").GetData(fGlobalFitData.core_x.lowLimit);
  bcore_x.GetChild("highLimit").GetData(fGlobalFitData.core_x.highLimit);

  utl::Branch bcore_y = bGlobalFitData.GetChild("core_y");
  bcore_y.GetChild("isFix").GetData(fGlobalFitData.core_y.isFix);
  bcore_y.GetChild("hasLimits").GetData(fGlobalFitData.core_y.hasLimits);
  bcore_y.GetChild("startValueFrom").GetData(fGlobalFitData.core_y.startValueFrom);
  bcore_y.GetChild("startValue").GetData(fGlobalFitData.core_y.startValue);
  bcore_y.GetChild("startError").GetData(fGlobalFitData.core_y.startError);
  bcore_y.GetChild("lowLimit").GetData(fGlobalFitData.core_y.lowLimit);
  bcore_y.GetChild("highLimit").GetData(fGlobalFitData.core_y.highLimit);

  utl::Branch bcore_z = bGlobalFitData.GetChild("core_z");
  bcore_z.GetChild("isFix").GetData(fGlobalFitData.core_z.isFix);
  bcore_z.GetChild("hasLimits").GetData(fGlobalFitData.core_z.hasLimits);
  bcore_z.GetChild("startValueFrom").GetData(fGlobalFitData.core_z.startValueFrom);
  bcore_z.GetChild("startValue").GetData(fGlobalFitData.core_z.startValue);
  bcore_z.GetChild("startError").GetData(fGlobalFitData.core_z.startError);
  bcore_z.GetChild("lowLimit").GetData(fGlobalFitData.core_z.lowLimit);
  bcore_z.GetChild("highLimit").GetData(fGlobalFitData.core_z.highLimit);

  utl::Branch bsigmaPlus = bGlobalFitData.GetChild("sigmaPlus");
  bsigmaPlus.GetChild("isFix").GetData(fGlobalFitData.sigmaPlus.isFix);
  bsigmaPlus.GetChild("hasLimits").GetData(fGlobalFitData.sigmaPlus.hasLimits);
  bsigmaPlus.GetChild("startValueFrom").GetData(fGlobalFitData.sigmaPlus.startValueFrom);
  bsigmaPlus.GetChild("startValue").GetData(fGlobalFitData.sigmaPlus.startValue);
  bsigmaPlus.GetChild("startError").GetData(fGlobalFitData.sigmaPlus.startError);
  bsigmaPlus.GetChild("lowLimit").GetData(fGlobalFitData.sigmaPlus.lowLimit);
  bsigmaPlus.GetChild("highLimit").GetData(fGlobalFitData.sigmaPlus.highLimit);

  utl::Branch bC2Theta = bGlobalFitData.GetChild("C2Theta");
  bC2Theta.GetChild("isFix").GetData(fGlobalFitData.C2Theta.isFix);
  bC2Theta.GetChild("hasLimits").GetData(fGlobalFitData.C2Theta.hasLimits);
  bC2Theta.GetChild("startValueFrom").GetData(fGlobalFitData.C2Theta.startValueFrom);
  bC2Theta.GetChild("startValue").GetData(fGlobalFitData.C2Theta.startValue);
  bC2Theta.GetChild("startError").GetData(fGlobalFitData.C2Theta.startError);
  bC2Theta.GetChild("lowLimit").GetData(fGlobalFitData.C2Theta.lowLimit);
  bC2Theta.GetChild("highLimit").GetData(fGlobalFitData.C2Theta.highLimit);

  utl::Branch bC1Theta = bGlobalFitData.GetChild("C1Theta");
  bC1Theta.GetChild("isFix").GetData(fGlobalFitData.C1Theta.isFix);
  bC1Theta.GetChild("hasLimits").GetData(fGlobalFitData.C1Theta.hasLimits);
  bC1Theta.GetChild("startValueFrom").GetData(fGlobalFitData.C1Theta.startValueFrom);
  bC1Theta.GetChild("startValue").GetData(fGlobalFitData.C1Theta.startValue);
  bC1Theta.GetChild("startError").GetData(fGlobalFitData.C1Theta.startError);
  bC1Theta.GetChild("lowLimit").GetData(fGlobalFitData.C1Theta.lowLimit);
  bC1Theta.GetChild("highLimit").GetData(fGlobalFitData.C1Theta.highLimit);

  utl::Branch bCTheta = bGlobalFitData.GetChild("CTheta");
  bCTheta.GetChild("isFix").GetData(fGlobalFitData.CTheta.isFix);
  bCTheta.GetChild("hasLimits").GetData(fGlobalFitData.CTheta.hasLimits);
  bCTheta.GetChild("startValueFrom").GetData(fGlobalFitData.CTheta.startValueFrom);
  bCTheta.GetChild("startValue").GetData(fGlobalFitData.CTheta.startValue);
  bCTheta.GetChild("startError").GetData(fGlobalFitData.CTheta.startError);
  bCTheta.GetChild("lowLimit").GetData(fGlobalFitData.CTheta.lowLimit);
  bCTheta.GetChild("highLimit").GetData(fGlobalFitData.CTheta.highLimit);

  utl::Branch bC1 = bGlobalFitData.GetChild("C1");
  bC1.GetChild("isFix").GetData(fGlobalFitData.C1.isFix);
  bC1.GetChild("hasLimits").GetData(fGlobalFitData.C1.hasLimits);
  bC1.GetChild("startValueFrom").GetData(fGlobalFitData.C1.startValueFrom);
  bC1.GetChild("startValue").GetData(fGlobalFitData.C1.startValue);
  bC1.GetChild("startError").GetData(fGlobalFitData.C1.startError);
  bC1.GetChild("lowLimit").GetData(fGlobalFitData.C1.lowLimit);
  bC1.GetChild("highLimit").GetData(fGlobalFitData.C1.highLimit);

  utl::Branch bC2 = bGlobalFitData.GetChild("C2");
  bC2.GetChild("isFix").GetData(fGlobalFitData.C2.isFix);
  bC2.GetChild("hasLimits").GetData(fGlobalFitData.C2.hasLimits);
  bC2.GetChild("startValueFrom").GetData(fGlobalFitData.C2.startValueFrom);
  bC2.GetChild("startValue").GetData(fGlobalFitData.C2.startValue);
  bC2.GetChild("startError").GetData(fGlobalFitData.C2.startError);
  bC2.GetChild("lowLimit").GetData(fGlobalFitData.C2.lowLimit);
  bC2.GetChild("highLimit").GetData(fGlobalFitData.C2.highLimit);

  utl::Branch bgamma = bGlobalFitData.GetChild("gamma");
  bgamma.GetChild("isFix").GetData(fGlobalFitData.gamma.isFix);
  bgamma.GetChild("hasLimits").GetData(fGlobalFitData.gamma.hasLimits);
  bgamma.GetChild("startValueFrom").GetData(fGlobalFitData.gamma.startValueFrom);
  bgamma.GetChild("startValue").GetData(fGlobalFitData.gamma.startValue);
  bgamma.GetChild("startError").GetData(fGlobalFitData.gamma.startError);
  bgamma.GetChild("lowLimit").GetData(fGlobalFitData.gamma.lowLimit);
  bgamma.GetChild("highLimit").GetData(fGlobalFitData.gamma.highLimit);

  utl::Branch bb = bGlobalFitData.GetChild("b");
  bb.GetChild("isFix").GetData(fGlobalFitData.b.isFix);
  bb.GetChild("hasLimits").GetData(fGlobalFitData.b.hasLimits);
  bb.GetChild("startValueFrom").GetData(fGlobalFitData.b.startValueFrom);
  bb.GetChild("startValue").GetData(fGlobalFitData.b.startValue);
  bb.GetChild("startError").GetData(fGlobalFitData.b.startError);
  bb.GetChild("lowLimit").GetData(fGlobalFitData.b.lowLimit);
  bb.GetChild("highLimit").GetData(fGlobalFitData.b.highLimit);

  utl::Branch bt0 = bGlobalFitData.GetChild("t0");
  bt0.GetChild("isFix").GetData(fGlobalFitData.t0.isFix);
  bt0.GetChild("hasLimits").GetData(fGlobalFitData.t0.hasLimits);
  bt0.GetChild("startValueFrom").GetData(fGlobalFitData.t0.startValueFrom);
  bt0.GetChild("startValue").GetData(fGlobalFitData.t0.startValue);
  bt0.GetChild("startError").GetData(fGlobalFitData.t0.startError);
  bt0.GetChild("lowLimit").GetData(fGlobalFitData.t0.lowLimit);
  bt0.GetChild("highLimit").GetData(fGlobalFitData.t0.highLimit);

  utl::Branch barrivalDirection_phi = bGlobalFitData.GetChild("arrivalDirection_phi");
  barrivalDirection_phi.GetChild("isFix").GetData(fGlobalFitData.arrivalDirection_phi.isFix);
  barrivalDirection_phi.GetChild("hasLimits").GetData(
      fGlobalFitData.arrivalDirection_phi.hasLimits);
  barrivalDirection_phi.GetChild("startValueFrom").GetData(
      fGlobalFitData.arrivalDirection_phi.startValueFrom);
  barrivalDirection_phi.GetChild("startValue").GetData(
      fGlobalFitData.arrivalDirection_phi.startValue);
  barrivalDirection_phi.GetChild("startError").GetData(
      fGlobalFitData.arrivalDirection_phi.startError);
  barrivalDirection_phi.GetChild("lowLimit").GetData(fGlobalFitData.arrivalDirection_phi.lowLimit);
  barrivalDirection_phi.GetChild("highLimit").GetData(
      fGlobalFitData.arrivalDirection_phi.highLimit);

  utl::Branch barrivalDirection_theta = bGlobalFitData.GetChild("arrivalDirection_theta");
  barrivalDirection_theta.GetChild("isFix").GetData(fGlobalFitData.arrivalDirection_theta.isFix);
  barrivalDirection_theta.GetChild("hasLimits").GetData(
      fGlobalFitData.arrivalDirection_theta.hasLimits);
  barrivalDirection_theta.GetChild("startValueFrom").GetData(
      fGlobalFitData.arrivalDirection_theta.startValueFrom);
  barrivalDirection_theta.GetChild("startValue").GetData(
      fGlobalFitData.arrivalDirection_theta.startValue);
  barrivalDirection_theta.GetChild("startError").GetData(
      fGlobalFitData.arrivalDirection_theta.startError);
  barrivalDirection_theta.GetChild("lowLimit").GetData(
      fGlobalFitData.arrivalDirection_theta.lowLimit);
  barrivalDirection_theta.GetChild("highLimit").GetData(
      fGlobalFitData.arrivalDirection_theta.highLimit);

  utl::Branch bshowerMax = bGlobalFitData.GetChild("showerMax");
  bshowerMax.GetChild("isFix").GetData(fGlobalFitData.showerMax.isFix);
  bshowerMax.GetChild("hasLimits").GetData(fGlobalFitData.showerMax.hasLimits);
  bshowerMax.GetChild("startValueFrom").GetData(fGlobalFitData.showerMax.startValueFrom);
  bshowerMax.GetChild("startValue").GetData(fGlobalFitData.showerMax.startValue);
  bshowerMax.GetChild("startError").GetData(fGlobalFitData.showerMax.startError);
  bshowerMax.GetChild("lowLimit").GetData(fGlobalFitData.showerMax.lowLimit);
  bshowerMax.GetChild("highLimit").GetData(fGlobalFitData.showerMax.highLimit);

  return eSuccess;
}

bool Sorter(std::vector<double> x, std::vector<double> y)
{
  return (x[7] < y[7]);
}

VModule::ResultFlag RdGlobalFit::Run(evt::Event& event)
{
  // Put your simulation, reconstruction or analysis
  // algorithm here.  This method is normally called
  // once per event, though the user can customize
  // when it is called using the sequencing XML file.

  // For debugging of the SDCore:

  INFO("RdGlobalFit::Run()");
  // Check if there are events at all
  if (!event.HasREvent()) {
    WARNING("No radio event found!");
    return eSuccess;
  }

  // initialize data structures
  REvent& rEvent = event.GetREvent();
  if (!event.HasRecShower()) {
    event.MakeRecShower();
  }
  if (!event.GetRecShower().HasRRecShower()) {
    event.GetRecShower().MakeRRecShower();
  }

  ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();

  const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

  // Initialising fit geometry.
  const Point coreRF = event.GetRecShower().GetRRecShower().GetCoordinateOrigin();
  utl::CoordinateSystemPtr RFCS = LocalCoordinateSystem::Create(coreRF);
  utl::Vector showerAxis = showerrrec.GetReferenceAxis(event);

  // B.Z. Here the used Coordinatesystem is defined! Needs OPtions like ReFCS, MCCS, SDCS
  // Later this is propagated to event.flocalCS
  fLocalCS = getlocalCSPtr(event);
  if (fLocalCS == utl::CoordinateSystemPtr()) {
    ERROR("fetching the localCSPtr according to Infile lead to NULLPtr  -> eBreakLoop");
    return (eBreakLoop);
  }

  const double zenith = showerAxis.GetTheta(fLocalCS);

  int Nant = showerrrec.GetParameter(revt::eNumberOfStationsWithPulseFound);
  if (Nant < 3)
    return eContinueLoop;

  //B.Z. Übergang in ShowerCoordinates!! localCS muss überall verwendet werden!
  Vector MagneticField = showerrrec.GetMagneticFieldVector();
  //utl::RadioGeometryUtilities transformation = utl::RadioGeometryUtilities(showerAxis, localCS,
  //                                                                         MagneticField);

  unsigned int signal_stations = 0;
  unsigned int silent_stations = 0;

  for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
      sIt != rEvent.CandidateStationsEnd(); ++sIt) {

    Station& station = *sIt;

    if (station.IsSaturated()) {   // skip saturated stations
      continue;
    }

    if (!station.HasSignal()) {   // skip silent station
      silent_stations++;
    } else {
      signal_stations++;
    }

  }

  // Only events with at least 3 stations are considered further.
  if (signal_stations < 3) {
    WARNING("Number of signalstations after rejection of saturated stations < 3 -> CONTINUE!");
    return eSuccess;
  }

  //Fetching the Data to be fitted by the Global Fit
  EventFitData eventData;
  eventData.fLocalCS = fLocalCS;
  std::vector<StationFitData> vStationFitData;
  unsigned int nrDataPoints = 0;
  double meanSignalTime = 0.;  // used as startvalue of t0 in timefit.
  double meanSignalTimeError = 0.;
  for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
      sIt != rEvent.CandidateStationsEnd(); ++sIt) {
    Station& station = *sIt;
    StationRecData& sRec = station.GetRecData();
    if (station.IsSaturated()) {   // skip saturated stations
      continue;
    }
    if (!station.HasSignal()) {   // skip silent stations
      continue;
    }
    // Fetching Position independent from coordinate system
    const Point sPos = rDetector.GetStation(station).GetPosition();

    //Fetching Energydensity for 2dLDFFit
    double integral = sRec.GetParameter(revt::eSignalEnergyFluenceMag);
    double uncertainty = sRec.GetParameterError(revt::eSignalEnergyFluenceMag);

    // the noise level is added to reflect the uncertainty of the model
    uncertainty = uncertainty * sqrt(1.79);

    //Calculating and Fetching time Parameter for ArrivalTimeFit
    double signalTime = sIt->GetRecData().GetParameter(revt::eSignalTime);
    double signalTimeError = sIt->GetRecData().GetParameterError(eSignalTime) * sqrt(1.65);

    if (station.HasSignal()) {
      ofstream errorDatei("ExperimentalErrors.txt", ios::out | ios::app);
      errorDatei << event.GetHeader().GetId() << "\t"
          << event.GetREvent().GetHeader().GetRunNumber() << "\t OfflineTimeerror: "
          << signalTimeError << "\t IntAmplError: " << uncertainty << endl;
      errorDatei.close();

      //Providing the global Fit Datastructure with data.
      StationFitData temp_StationFitData;
      temp_StationFitData.stationId = sIt->GetId();
      temp_StationFitData.stationPosition = sPos;
      temp_StationFitData.energyFluence = integral;
      temp_StationFitData.energyFluenceError = uncertainty;
      temp_StationFitData.signalTime = signalTime;
      temp_StationFitData.signalTimeError = signalTimeError;
      vStationFitData.push_back(temp_StationFitData);
      ++nrDataPoints;
      meanSignalTime = meanSignalTime + signalTime;
      meanSignalTimeError = meanSignalTimeError + signalTimeError;
    }
  }

  meanSignalTime /= nrDataPoints;
  meanSignalTimeError /= nrDataPoints;

  RdGlobalFitMinimizationCriterion fitFunction = RdGlobalFitMinimizationCriterion(fFitConfig,
                                                                                  eventData,
                                                                                  vStationFitData,
                                                                                  MagneticField,
                                                                                  fLDFConstsTable,
                                                                                  fcalcLDFConsts);

  //Needs to be uncomented when deleting the original ldffit above

  double A_scale = fLDFConstsTable.A_scale;

  if ((zenith) < 10. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_0_10;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_0_10;
  } else if ((zenith) < 20. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_10_20;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_10_20;
  } else if ((zenith) < 30. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_20_30;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_20_30;
  } else if ((zenith) < 40. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_30_40;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_30_40;
  } else if ((zenith) < 50. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_40_50;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_40_50;
  } else if ((zenith) < 55. * utl::degree) {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_50_55;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_50_55;
    A_scale = fLDFConstsTable.A_scale_50_55;
    fcalcLDFConsts.C0 = 0.46;
  } else {
    fcalcLDFConsts.C1theta = fLDFConstsTable.C1theta_55;
    fcalcLDFConsts.C2theta = fLDFConstsTable.C2theta_55;
    A_scale = fLDFConstsTable.A_scale_55;
    fcalcLDFConsts.C0 = 0.71;
  }

  // Depending on the infile settings fGlobalFitData:Values are overwritten for core, direction.
  setGlFitDataStartValues(event);

  ROOT::Minuit2::MnUserParameters upar;

  fGlobalFitData.C2Theta = FitData(false, false, fcalcLDFConsts.C2theta,
                                   0.1 * fcalcLDFConsts.C2theta, fcalcLDFConsts.C2theta,
                                   1.1 * fcalcLDFConsts.C2theta);
  fGlobalFitData.C1Theta = FitData(false, false, fcalcLDFConsts.C1theta,
                                   0.1 * fcalcLDFConsts.C1theta, fcalcLDFConsts.C1theta,
                                   1.1 * fcalcLDFConsts.C1theta);
  fGlobalFitData.CTheta = FitData(false, false, fLDFConstsTable.A_scale,
                                  0.1 * fLDFConstsTable.A_scale, A_scale,
                                  1.1 * fLDFConstsTable.A_scale);
  fGlobalFitData.C1 = FitData(false, false, fLDFConstsTable.C3, 0.1 * fLDFConstsTable.C3,
                              fLDFConstsTable.C3, 1.1 * fLDFConstsTable.C3);
  fGlobalFitData.C2 = FitData(false, false, fLDFConstsTable.C4, 0.1 * fLDFConstsTable.C4,
                              fLDFConstsTable.C4, 1.1 * fLDFConstsTable.C4);

  // Filling minuit Parameters
  // Even though, t0 startvalue is set here also, is always overwritten by the meanSignalTime later.
  upar.Add("Aplus", fGlobalFitData.Aplus.startValue, fGlobalFitData.Aplus.startError);
  upar.Add("core_x", fGlobalFitData.core_x.startValue, fGlobalFitData.core_x.startError);
  upar.Add("core_y", fGlobalFitData.core_y.startValue, fGlobalFitData.core_y.startError);
  upar.Add("core_z", fGlobalFitData.core_z.startValue, fGlobalFitData.core_z.startError);
  upar.Add("sigmaPlus", fGlobalFitData.sigmaPlus.startValue, fGlobalFitData.sigmaPlus.startError);
  upar.Add("C2Theta", fGlobalFitData.C2Theta.startValue, fGlobalFitData.C2Theta.startError);
  upar.Add("C1Theta", fGlobalFitData.C1Theta.startValue, fGlobalFitData.C1Theta.startError);
  upar.Add("CTheta", fGlobalFitData.CTheta.startValue, fGlobalFitData.CTheta.startError);
  upar.Add("C1", fGlobalFitData.C1.startValue, fGlobalFitData.C1.startError);
  upar.Add("C2", fGlobalFitData.C2.startValue, fGlobalFitData.C2.startError);
  upar.Add("gamma", fGlobalFitData.gamma.startValue, fGlobalFitData.gamma.startError);
  upar.Add("showerMax", fGlobalFitData.showerMax.startValue, fGlobalFitData.showerMax.startError);
  upar.Add("b", fGlobalFitData.b.startValue, fGlobalFitData.b.startError);
  upar.Add("t0", fGlobalFitData.t0.startValue, fGlobalFitData.t0.startError);
  upar.Add("arrivalDirection_phi", fGlobalFitData.arrivalDirection_phi.startValue,
           fGlobalFitData.arrivalDirection_phi.startError);
  upar.Add("arrivalDirection_theta", fGlobalFitData.arrivalDirection_theta.startValue,
           fGlobalFitData.arrivalDirection_theta.startError);

  //Fix Parameters which are given by Johannes 2dLDF studies. Compare ...Paper...
  upar.Fix("CTheta");
  upar.Fix("C2Theta");
  upar.Fix("C1Theta");
  upar.Fix("C1");
  upar.Fix("C2");

  if (fFitConfig.fixCore == true) {
    upar.Fix("core_x");
    upar.Fix("core_y");
    upar.Fix("core_z");
  }

  if (fFitConfig.fitTwoDLDF == false) {
    upar.Fix("sigmaPlus");
    upar.Fix("Aplus");
  }

  if (fFitConfig.fitArrivalTime == false) {
    upar.Fix("gamma");
    upar.Fix("b");
    upar.Fix("t0");
  }

  // Filling minuit Parameters
  if (fGlobalFitData.Aplus.isFix)
    upar.Fix("Aplus");
  if (fGlobalFitData.core_x.isFix)
    upar.Fix("core_x");
  if (fGlobalFitData.core_y.isFix)
    upar.Fix("core_y");
  if (fGlobalFitData.core_z.isFix)
    upar.Fix("core_z");
  if (fGlobalFitData.sigmaPlus.isFix)
    upar.Fix("sigmaPlus");
  if (fGlobalFitData.C2Theta.isFix)
    upar.Fix("C2Theta");
  if (fGlobalFitData.C1Theta.isFix)
    upar.Fix("C1Theta");
  if (fGlobalFitData.CTheta.isFix)
    upar.Fix("CTheta");
  if (fGlobalFitData.C1.isFix)
    upar.Fix("C1");
  if (fGlobalFitData.C2.isFix)
    upar.Fix("C2");
  if (fGlobalFitData.gamma.isFix)
    upar.Fix("gamma");
  if (fGlobalFitData.b.isFix)
    upar.Fix("b");
  if (fGlobalFitData.t0.isFix)
    upar.Fix("t0");
  if (fGlobalFitData.arrivalDirection_phi.isFix)
    upar.Fix("arrivalDirection_phi");
  if (fGlobalFitData.arrivalDirection_theta.isFix)
    upar.Fix("arrivalDirection_theta");
  if (fGlobalFitData.showerMax.isFix)
    upar.Fix("showerMax");

  //Setting user Parameter Limits
  if (fGlobalFitData.Aplus.hasLimits)
    upar.SetLimits("Aplus", fGlobalFitData.Aplus.lowLimit, fGlobalFitData.Aplus.highLimit);
  if (fGlobalFitData.core_x.hasLimits)
    upar.SetLimits("core_x", fGlobalFitData.core_x.lowLimit, fGlobalFitData.core_x.highLimit);
  if (fGlobalFitData.core_y.hasLimits)
    upar.SetLimits("core_y", fGlobalFitData.core_y.lowLimit, fGlobalFitData.core_y.highLimit);
  if (fGlobalFitData.core_z.hasLimits)
    upar.SetLimits("core_z", fGlobalFitData.core_z.lowLimit, fGlobalFitData.Aplus.highLimit);
  if (fGlobalFitData.sigmaPlus.hasLimits)
    upar.SetLimits("sigmaPlus", fGlobalFitData.sigmaPlus.lowLimit,
                   fGlobalFitData.sigmaPlus.highLimit);
  if (fGlobalFitData.C2Theta.hasLimits)
    upar.SetLimits("C2Theta", fGlobalFitData.C2Theta.lowLimit, fGlobalFitData.Aplus.highLimit);
  if (fGlobalFitData.C1Theta.hasLimits)
    upar.SetLimits("C1Theta", fGlobalFitData.C1Theta.lowLimit, fGlobalFitData.C1Theta.highLimit);
  if (fGlobalFitData.CTheta.hasLimits)
    upar.SetLimits("CTheta", fGlobalFitData.CTheta.lowLimit, fGlobalFitData.CTheta.highLimit);
  if (fGlobalFitData.C1.hasLimits)
    upar.SetLimits("C1", fGlobalFitData.C1.lowLimit, fGlobalFitData.C1.highLimit);
  if (fGlobalFitData.C2.hasLimits)
    upar.SetLimits("C2", fGlobalFitData.C2.lowLimit, fGlobalFitData.C2.highLimit);
  if (fGlobalFitData.gamma.hasLimits)
    upar.SetLimits("gamma", fGlobalFitData.gamma.lowLimit, fGlobalFitData.gamma.highLimit);
  if (fGlobalFitData.b.hasLimits)
    upar.SetLimits("b", fGlobalFitData.b.lowLimit, fGlobalFitData.b.highLimit);
  if (fGlobalFitData.t0.hasLimits)
    upar.SetLimits("t0", fGlobalFitData.t0.lowLimit, fGlobalFitData.t0.highLimit);
  if (fGlobalFitData.arrivalDirection_phi.hasLimits)
    upar.SetLimits("arrivalDirection_phi", fGlobalFitData.arrivalDirection_phi.lowLimit,
                   fGlobalFitData.arrivalDirection_phi.highLimit);
  if (fGlobalFitData.arrivalDirection_theta.hasLimits)
    upar.SetLimits("arrivalDirection_theta", fGlobalFitData.arrivalDirection_theta.lowLimit,
                   fGlobalFitData.arrivalDirection_theta.highLimit);
  if (fGlobalFitData.showerMax.hasLimits)
    upar.SetLimits("showerMax", fGlobalFitData.showerMax.lowLimit,
                   fGlobalFitData.showerMax.highLimit);

  // Here, t0 startvalue is always overwritten by the meansignaltime.
  upar.SetValue("t0", meanSignalTime);
  upar.SetLimits("t0", meanSignalTime + 30, meanSignalTime - 30);

  // Start the actual FIT in 5 stages
  const ROOT::Minuit2::MnStrategy strat(2);
  //Stage 1: A and sigma are free
  int NDF = 0;
  ROOT::Minuit2::MnMigrad LDFPrefit(fitFunction, upar, strat);
  if (fOrigFitConfig.fitTwoDLDF == true) {
    fFitConfig.fitTwoDLDF = true;
    fFitConfig.fitArrivalTime = false;
    upar.Fix("core_x");
    upar.Fix("core_y");
    ROOT::Minuit2::MnMigrad LDFPrefit(fitFunction, upar, strat);
    LDFPrefit.Fix("arrivalDirection_phi");
    LDFPrefit.Fix("arrivalDirection_theta");
    LDFPrefit.Fix("gamma");
    LDFPrefit.Fix("t0");
    NDF = calcNDF(LDFPrefit.Parameters(), nrDataPoints);
    if (NDF < 1) {
      INFO("NDF less than 1 -> eContinueLoop");
      return eContinueLoop;
    }
    ROOT::Minuit2::FunctionMinimum LDFPrefitMin = LDFPrefit();
    cout << "Fit result Stage 1:" << endl << LDFPrefitMin << endl;
  }

  //Stage 2: Full energy fluence fit (arrival direction stays fix)
  ROOT::Minuit2::MnMigrad LDFPrefit2 = LDFPrefit;
  if (fOrigFitConfig.fitTwoDLDF == true) {
    fFitConfig.fitTwoDLDF = true;
    fFitConfig.fitArrivalTime = false;
    LDFPrefit2.Fix("arrivalDirection_phi");
    LDFPrefit2.Fix("arrivalDirection_theta");
    LDFPrefit2.Fix("gamma");
    LDFPrefit2.Fix("t0");
    if (!fGlobalFitData.core_x.isFix)
      LDFPrefit2.Release("core_x");
    if (!fGlobalFitData.core_y.isFix)
      LDFPrefit2.Release("core_y");
    NDF = calcNDF(LDFPrefit2.Parameters(), nrDataPoints);
    if (NDF < 1) {
      INFO("NDF less than 1 -> eContinueLoop");
      return eContinueLoop;
    }
    ROOT::Minuit2::FunctionMinimum LDFPrefitMin2 = LDFPrefit2();
    cout << "Fit result Stage 2:" << endl << LDFPrefitMin2 << endl;
  }

  //Stage 3: Only impact parameter t0 is free
  ROOT::Minuit2::MnMigrad PrefitT0 = LDFPrefit2;
  if (fOrigFitConfig.fitArrivalTime == true) {
    //still use LDF chi2 to bind the fit improvement to 2dLDF already.
    if (fOrigFitConfig.fitTwoDLDF == true) {
      fFitConfig.fitTwoDLDF = true;
      PrefitT0.Fix("Aplus");
      PrefitT0.Fix("sigmaPlus");
    } else {
      fFitConfig.fitTwoDLDF = false;
    }
    fFitConfig.fitArrivalTime = true;
    PrefitT0.Fix("arrivalDirection_phi");
    PrefitT0.Fix("arrivalDirection_theta");
    PrefitT0.Fix("gamma");
    PrefitT0.Fix("core_x");
    PrefitT0.Fix("core_y");
    if (!fGlobalFitData.t0.isFix)
      PrefitT0.Release("t0");
    NDF = calcNDF(PrefitT0.Parameters(), nrDataPoints);
    if (NDF < 1) {
      INFO("NDF less than 1 -> eContinueLoop");
      return eContinueLoop;
    }
    ROOT::Minuit2::FunctionMinimum PrefitT0Min = PrefitT0();
    cout << "Fit result Stage 3:" << endl << PrefitT0Min << endl;
  }

  //Stage 4: Arrivaltime fit (core taken from energy fluence fit)
  ROOT::Minuit2::MnMigrad PrefitT = PrefitT0;
  if (fOrigFitConfig.fitArrivalTime == true) {
    //still use LDF chi2 to base the fit improvement on 2dLDF chi2 already.
    if (fOrigFitConfig.fitTwoDLDF == true) {
      fFitConfig.fitTwoDLDF = true;
      PrefitT.Fix("Aplus");
      PrefitT.Fix("sigmaPlus");
    } else {
      fFitConfig.fitTwoDLDF = false;
    }
    fFitConfig.fitArrivalTime = true;
    if (!fGlobalFitData.arrivalDirection_phi.isFix)
      PrefitT.Release("arrivalDirection_phi");
    if (!fGlobalFitData.arrivalDirection_theta.isFix)
      PrefitT.Release("arrivalDirection_theta");
    if (!fGlobalFitData.gamma.isFix)
      PrefitT.Release("gamma");
    if (!fGlobalFitData.t0.isFix)
      PrefitT.Release("t0");

    PrefitT.Fix("core_x");
    PrefitT.Fix("core_y");
    NDF = calcNDF(PrefitT.Parameters(), nrDataPoints);
    if (NDF < 1) {
      INFO("NDF less than 1 -> eContinueLoop");
      return eContinueLoop;
    }
    ROOT::Minuit2::FunctionMinimum PrefitTMin = PrefitT();
    cout << "Fit result Stage 4:" << endl << PrefitTMin << endl;
  }

  //Stage 5: Finally, all Parameters are fitted according the xml specifications.
  ROOT::Minuit2::MnMigrad Fit = PrefitT;
  if (fOrigFitConfig.fitTwoDLDF == true) {
    fFitConfig.fitTwoDLDF = true;
    if (!fGlobalFitData.Aplus.isFix)
      Fit.Release("Aplus");
    if (!fGlobalFitData.sigmaPlus.isFix)
      Fit.Release("sigmaPlus");
  } else {
    fFitConfig.fitTwoDLDF = false;
  }
  if (fOrigFitConfig.fitArrivalTime == true) {
    fFitConfig.fitArrivalTime = true;
    if (!fGlobalFitData.arrivalDirection_phi.isFix)
      Fit.Release("arrivalDirection_phi");
    if (!fGlobalFitData.arrivalDirection_theta.isFix)
      Fit.Release("arrivalDirection_theta");
    if (!fGlobalFitData.gamma.isFix)
      Fit.Release("gamma");
    if (!fGlobalFitData.t0.isFix)
      Fit.Release("t0");
  } else {
    fFitConfig.fitArrivalTime = false;
  }
  Fit.Fix("core_x");
  Fit.Fix("core_y");
  if (!fGlobalFitData.core_x.isFix)
    Fit.Release("core_x");
  if (!fGlobalFitData.core_y.isFix)
    Fit.Release("core_y");
  NDF = calcNDF(Fit.Parameters(), nrDataPoints);
  if (NDF < 1) {
    INFO("NDF less than 1 -> eContinueLoop");
    return eContinueLoop;
  }
  if (fOrigFitConfig.fitArrivalTime == true) {
    if (!fGlobalFitData.arrivalDirection_phi.isFix)
      Fit.Release("arrivalDirection_phi");
    if (!fGlobalFitData.arrivalDirection_theta.isFix)
      Fit.Release("arrivalDirection_theta");
  }
  ROOT::Minuit2::FunctionMinimum minimum = Fit(10000, 0.0001);
  cout << "Final fit result (Stage 5):" << endl << event.GetHeader().GetId() << "\t"
      << event.GetREvent().GetHeader().GetRunNumber() << "\t" << endl << "Nr. St = " << nrDataPoints
      << endl << meanSignalTime << meanSignalTime << endl << minimum << endl;
  std::cout << "Nr fct calls: " << minimum.NFcn() << std::endl;

  // get the Likelihoods of the Minimum.
  fitFunction(minimum.UserParameters().Params());
  double minTwoDLDFLikelihood = fitFunction.getTwoDLDFLikelyhood();
  double minArrivalTimeLikelihood = fitFunction.getArrivalTimeLikelihood();

  // Save core contour as text file.
  // 1.) range from -300 m to 300 m
  if (fFitConfig.SaveCoreContourData) {
    std::vector<double> contourParams = minimum.UserParameters().Params();
    std::vector<std::vector<double> > contourVec;
    ofstream Zieldatei("CoreContour.dat", ios::out | ios::app);
    for (int x = -300; x < 301; x = x + 10) {
      for (int y = -300; y < 301; y = y + 10) {
        std::vector<double> contourtriple;
        contourParams[1] = x;
        contourParams[2] = y;
        contourtriple.push_back(x);
        contourtriple.push_back(y);
        contourtriple.push_back(fitFunction(contourParams));
        contourVec.push_back(contourtriple);
        Zieldatei << contourtriple[0] << "\t" << contourtriple[1] << "\t" << contourtriple[2]
            << "\t" << event.GetHeader().GetId() << std::endl;
      }
    }
    Zieldatei.close();
  }
  // 2.) Zoom into the range from -5 m to 5 m
  if (fFitConfig.SaveCoreContourData) {
    std::vector<double> contourParams = minimum.UserParameters().Params();
    double min_x = minimum.UserParameters().Value("core_x");
    double min_y = minimum.UserParameters().Value("core_y");
    std::vector<std::vector<double> > contourVec;
    ofstream Zieldatei("CoreContourZoom.dat", ios::out | ios::app);
    for (double x = min_x - 5; x < min_x + 5; x = x + 0.25) {
      for (double y = min_y - 5; y < min_y + 5; y = y + 0.25) {
        std::vector<double> contourtriple;
        contourParams[1] = x;
        contourParams[2] = y;
        contourtriple.push_back(x);
        contourtriple.push_back(y);
        contourtriple.push_back(fitFunction(contourParams));
        contourVec.push_back(contourtriple);
        Zieldatei << contourtriple[0] << "\t" << contourtriple[1] << "\t" << contourtriple[2]
            << "\t" << event.GetHeader().GetId() << std::endl;
      }
    }
    Zieldatei.close();
  }

  //Scan & Plot the parameters to the console (partly via std::cout)
  if (fFitConfig.printMinuitScans) {
    std::vector<std::string> parameterNames;
    parameterNames.push_back("showerMax");
    parameterNames.push_back("core_x");
    parameterNames.push_back("core_y");
    parameterNames.push_back("core_z");
    parameterNames.push_back("t0");
    parameterNames.push_back("gamma");
    parameterNames.push_back("arrivalDirection_theta");
    parameterNames.push_back("arrivalDirection_phi");
    parameterNames.push_back("Aplus");
    parameterNames.push_back("sigmaPlus");
    std::cout << "Defined scan parameters" << std::endl << flush;
    std::vector<std::pair<std::vector<std::pair<double, double> >, std::string> > scanPlots(
        scanParameters(fitFunction, minimum, parameterNames));
    printPlotsVec(scanPlots);
  }

  ROOT::Minuit2::MnUserParameterState minUstate = minimum.UserState();
  ROOT::Minuit2::MnUserCovariance cov = minUstate.Covariance();

  utl::Vector fittedShoweraxis(
      sin(minUstate.Value("arrivalDirection_theta")) * cos(minUstate.Value("arrivalDirection_phi")),
      sin(minUstate.Value("arrivalDirection_theta")) * sin(minUstate.Value("arrivalDirection_phi")),
      cos(minUstate.Value("arrivalDirection_theta")), eventData.fLocalCS);
  double sineAlpha = RadioGeometryUtilities::GetLorentzVector(fittedShoweraxis,
                                                              showerrrec.GetMagneticFieldVector())
      .GetMag();
  double sigma = minUstate.Value("sigmaPlus");
  double sigma_error = minUstate.Error("sigmaPlus");
  double A = minUstate.Value("Aplus");
  double A_error = minUstate.Error("Aplus");
  double sigmaA_error = 0;  //The global Fit does not now this error

  double gl_epsilon =
      A * TMath::Pi()
          * (pow(sigma, 2)
              - fcalcLDFConsts.C0
                  * TMath::Exp(2 * fLDFConstsTable.C3 + 2 * fLDFConstsTable.C4 * sigma));
  double gl_epsilonError = pow(
      pow(TMath::Pi(), 2)
          * pow(
              fcalcLDFConsts.C0
                  * TMath::Exp(2 * fLDFConstsTable.C3 + 2 * fLDFConstsTable.C4 * sigma)
                  - pow(sigma, 2),
              2) * pow(A_error, 2)
          + pow(A, 2) * pow(TMath::Pi(), 2)
              * pow(
                  (2 * fcalcLDFConsts.C0 * fLDFConstsTable.C4
                      * TMath::Exp(2 * fLDFConstsTable.C4 * sigma + 2 * fLDFConstsTable.C3)
                      - 2 * sigma),
                  2) * pow(sigma_error, 2)
          - TMath::Pi()
              * (2 * fcalcLDFConsts.C0 * fLDFConstsTable.C4
                  * TMath::Exp(2 * fLDFConstsTable.C4 * sigma + 2 * fLDFConstsTable.C3) - 2 * sigma)
              * sigmaA_error,
      0.5);

  const double d = 2.3137 * 1e14;  // energy calibration parameters from Aab et al., PRD 93, 122005 (2016)
  const double c = 0.505;
  double gl_energy = d * std::pow(gl_epsilon / (sineAlpha * sineAlpha), c);
  double gl_energyError = c * d * std::pow(1 / (sineAlpha * sineAlpha), c)
      * std::pow(gl_epsilon, c - 1) * gl_epsilonError;

  Point coreFit(minUstate.Value("core_x"), minUstate.Value("core_y"), minUstate.Value("core_z"),
                eventData.fLocalCS);
  Point coreFitError(minUstate.Error("core_x"), minUstate.Error("core_y"),
                     minUstate.Error("core_z"), fLocalCS);

  // should be GetSiteCoordinateSystem(); but RdFiller.cc::AddRadioDetector() sets
  // the positions against the comments in ADST Detector to referenceCS.
  utl::CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
  Point coreInSiteCS = getCoreInCSPlane(coreFit, fittedShoweraxis, siteCS);

  // Get the referenceCS (which is the same as siteCS, see above)
  utl::CoordinateSystemPtr referenceCS =
      det::Detector::GetInstance().GetReferenceCoordinateSystem();
  Point coreInReferenceCS = getCoreInCSPlane(coreFit, fittedShoweraxis, referenceCS);

  // At one point, one can think of moving this to FitModels.cc
  double RmaxFrSigma = 65 * minUstate.Value("sigmaPlus") - 3129;  // in m
  double RmaxFrGamma = 200 * minUstate.Value("gamma") - 17700;  // in m
  double zenitInLocalCS = fittedShoweraxis.GetTheta(getlocalCSPtr(event));
  double XmaxFromSigma = XmaxFromRmax(event, RmaxFrSigma, zenitInLocalCS);
  double XmaxFromGamma = XmaxFromRmax(event, RmaxFrGamma, zenitInLocalCS);

  ShowerRRecData& showerReco = event.GetRecShower().GetRRecShower();
  showerReco.SetParameter(eGlobalFitA, minUstate.Value("Aplus"));
  showerReco.SetParameterError(eGlobalFitA, minUstate.Error("Aplus"));
  showerReco.SetParameter(eGlobalFitB, minUstate.Value("b"));
  showerReco.SetParameterError(eGlobalFitB, minUstate.Error("b"));
  showerReco.SetParameter(eGlobalFitChi2, minUstate.Fval());
  showerReco.SetParameter(eGlobalFitChi2A, minTwoDLDFLikelihood);
  showerReco.SetParameter(eGlobalFitChi2T, minArrivalTimeLikelihood);
  showerReco.SetParameter(eGlobalFitE, gl_energy);
  showerReco.SetParameterError(eGlobalFitE, gl_energyError);
  showerReco.SetParameter(eGlobalFitGamma, minUstate.Value("gamma"));
  showerReco.SetParameterError(eGlobalFitGamma, minUstate.Error("gamma"));
  showerReco.SetParameter(eGlobalFitPhi, fittedShoweraxis.GetPhi(RFCS));
  showerReco.SetParameter(eGlobalFitTheta, fittedShoweraxis.GetTheta(RFCS));
  showerReco.SetParameter(eGlobalFitRadE, gl_epsilon);
  showerReco.SetParameterError(eGlobalFitRadE, gl_epsilonError);
  showerReco.SetParameter(eGlobalFitSigma, minUstate.Value("sigmaPlus"));
  showerReco.SetParameterError(eGlobalFitSigma, minUstate.Error("sigmaPlus"));
  showerReco.SetParameter(eGlobalFitT0, minUstate.Value("t0"));
  showerReco.SetParameterError(eGlobalFitT0, minUstate.Error("t0"));
  showerReco.SetParameter(eGlobalFitX, coreInReferenceCS.GetX(referenceCS));
  showerReco.SetParameterError(eGlobalFitX, minUstate.Error("core_x"));
  showerReco.SetParameter(eGlobalFitY, coreInReferenceCS.GetY(referenceCS));
  showerReco.SetParameterError(eGlobalFitY, minUstate.Error("core_y"));
  showerReco.SetParameter(eGlobalFitZ, coreInReferenceCS.GetY(referenceCS));
  showerReco.SetParameterError(eGlobalFitZ, minUstate.Error("core_z"));
  showerReco.SetParameter(eGlobalFitValMin, minimum.IsValid());
  showerReco.SetParameter(eGlobalFitRmaxFrSigma, RmaxFrSigma);
  showerReco.SetParameter(eGlobalFitXmaxFrSigma, XmaxFromSigma);
  showerReco.SetParameter(eGlobalFitRmaxFrGamma, RmaxFrGamma);
  showerReco.SetParameter(eGlobalFitXmaxFrGamma, XmaxFromGamma);
  showerReco.SetParameter(eGlobalFitXsitCS, coreInSiteCS.GetX(siteCS));
  showerReco.SetParameter(eGlobalFitYsitCS, coreInSiteCS.GetY(siteCS));
  showerReco.SetParameter(eGlobalFitZsitCS, coreInSiteCS.GetZ(siteCS));
  showerReco.SetParameter(eGlobalFitNDF, NDF);
  showerReco.SetParameter(eGlobalFitNDFA, 0.);
  showerReco.SetParameter(eGlobalFitNDFT, 0.);
  return (eSuccess);
}

std::vector<std::pair<std::vector<std::pair<double, double> >, std::string> > RdGlobalFit::scanParameters(
    RdGlobalFitMinimizationCriterion fitFct,
    ROOT::Minuit2::FunctionMinimum minimum/*,  ROOT::Minuit2::MnStrategy strat*/,
    std::vector<std::string> parameterToScan)
{
  ROOT::Minuit2::MnStrategy strat(2);
  std::vector<std::pair<std::vector<std::pair<double, double> >, std::string> > scanResults;
  ROOT::Minuit2::MnScan scan(fitFct, minimum.UserState(), strat);
  for (std::vector<std::string>::iterator it = parameterToScan.begin(); it != parameterToScan.end();
      ++it) {
    std::cout << "Actual scaning of parameters." << std::endl << flush;
    int parameterIndex = minimum.UserParameters().Index(*it);
    double low = minimum.UserParameters().Value(*it) - 10 * minimum.UserParameters().Error(*it);
    double up = minimum.UserParameters().Value(*it) + 10 * minimum.UserParameters().Error(*it);
    std::vector<std::pair<double, double> > points = scan.Scan(parameterIndex, 50, low, up);
    std::pair<std::vector<std::pair<double, double> >, std::string> parameterScan(points, *it);
    scanResults.push_back(parameterScan);
  }
  return scanResults;
}

void RdGlobalFit::printPlotsVec(
    std::vector<std::pair<std::vector<std::pair<double, double> >, std::string> > plotsVec)
{
  for (std::vector<std::pair<std::vector<std::pair<double, double> >, std::string> >::iterator it =
      plotsVec.begin(); it != plotsVec.end(); ++it) {
    std::cout << "Plotting core scans" << std::endl << flush;
    std::vector<std::pair<double, double> > plotPoints = it->first;
    std::string plotComment = it->second;
    ROOT::Minuit2::MnPlot plot;
    cout << "Comment: " << endl << plotComment << endl;
    plot(plotPoints);
  }
}

utl::CoordinateSystemPtr RdGlobalFit::getlocalCSPtr(evt::Event& event)
{
  if (fFitConfig.useCS == "RLocal") {
    const Point coreRLocal = event.GetRecShower().GetRRecShower().GetCoordinateOrigin();
    utl::CoordinateSystemPtr RLocalCS = LocalCoordinateSystem::Create(coreRLocal);
    return RLocalCS;
  }

  if (fFitConfig.useCS == "SD") {
    const Point coreSD = event.GetRecShower().GetSRecShower().GetCorePosition();  // use SD-Core as reference for fit
    utl::CoordinateSystemPtr SDCS = LocalCoordinateSystem::Create(coreSD);  // use SD-Core as reference for fit
    return SDCS;
  }

  if (fFitConfig.useCS == "MC") {
    const Point coreMC = event.GetSimShower().GetPosition();
    utl::CoordinateSystemPtr MCCS = LocalCoordinateSystem::Create(coreMC);
    return MCCS;
  }

  if (fFitConfig.useCS == "REF") {
    const Point coreRF = event.GetRecShower().GetRRecShower().GetCoordinateOrigin();
    utl::CoordinateSystemPtr RFCS = LocalCoordinateSystem::Create(coreRF);
    return RFCS;
  }

  if (fFitConfig.useCS == "RD") {
    const Point coreRD = event.GetRecShower().GetRRecShower().GetCorePosition();
    utl::CoordinateSystemPtr RDCS = LocalCoordinateSystem::Create(coreRD);
    return RDCS;
  }

  return (utl::CoordinateSystemPtr());
}

// Is not using the right atmospheric model yet. Its results are not validated yet.
double RdGlobalFit::XmaxFromRmax(evt::Event& event, double rmax, double zenithInLocalCS)
{
  det::Detector& detector = det::Detector::GetInstance();
  detector.Update(event.GetHeader().GetTime().GetGPSSecond());
  const atm::Atmosphere& theAtm = det::Detector::GetInstance().GetAtmosphere();
  double hauger = 1564;  //m
  atm::ProfileResult depthVsheight = theAtm.EvaluateDepthVsHeight();
  double Xmax = ((depthVsheight.Y(hauger + rmax * m * cos(zenithInLocalCS))) / (g / cm / cm))
      / cos(zenithInLocalCS);  //Xmax
  return Xmax;
}
void RdGlobalFit::setGlFitDataStartValues(evt::Event& event)
{
  setGlFitDataCore(event);
  setGlFitDataDirection(event);
}

utl::Point RdGlobalFit::adaptSDCoreTofLocalCS(evt::Event& event)
{
  const Point originFLocalCS(0., 0., 0., fLocalCS);
  const utl::Vector zAxisFLocalCS(0., 0., 1., fLocalCS);
  const utl::Plane XYPlaneFLocalCS(originFLocalCS, zAxisFLocalCS);

  const Point coreSD = event.GetRecShower().GetSRecShower().GetCorePosition();
  const utl::Vector& axisSD = event.GetRecShower().GetSRecShower().GetAxis();
  const utl::Line SDCoreAxisLine(coreSD, axisSD);

  utl::Point CorrCoreSD(utl::Intersection(SDCoreAxisLine, XYPlaneFLocalCS));
  return CorrCoreSD;
}

utl::Point RdGlobalFit::getCoreInCSPlane(utl::Point core, utl::Vector& axis,
                                         utl::CoordinateSystemPtr targetCS)
{
  const Point originTargetCS(0., 0., 0., targetCS);
  const utl::Vector zAxisTargetCS(0., 0., 1., targetCS);
  const utl::Plane XYPlaneTargetCS(originTargetCS, zAxisTargetCS);

  const utl::Line CoreAxisLine(core, axis);

  utl::Point CoreInTargetCSPlane(utl::Intersection(CoreAxisLine, XYPlaneTargetCS));
  return CoreInTargetCSPlane;
}

void RdGlobalFit::setGlFitDataCore(evt::Event& event)
{
  /* This is dangerous. Offline may crash when asked e.g. for the RD core, even if it is not
     used at all later on. Setting start values to MC missing.
     Initiating the x, y, and z component with different SD/RD/MC reconstructions seems odd.
  */

  const Point coreRD = event.GetRecShower().GetRRecShower().GetCorePosition();
  const Point coreSD = adaptSDCoreTofLocalCS(event);

  if (fGlobalFitData.core_x.startValueFrom == "RD") {
    fGlobalFitData.core_x.startValue = coreRD.GetX(fLocalCS);
  }
  if (fGlobalFitData.core_x.startValueFrom == "SD") {
    fGlobalFitData.core_x.startValue = coreSD.GetX(fLocalCS);
  }

  if (fGlobalFitData.core_y.startValueFrom == "RD") {
    fGlobalFitData.core_y.startValue = coreRD.GetY(fLocalCS);
  }
  if (fGlobalFitData.core_y.startValueFrom == "SD") {
    fGlobalFitData.core_y.startValue = coreSD.GetY(fLocalCS);
  }

  if (fGlobalFitData.core_z.startValueFrom == "RD") {
    fGlobalFitData.core_z.startValue = coreRD.GetZ(fLocalCS);
  }
  if (fGlobalFitData.core_z.startValueFrom == "SD") {
    fGlobalFitData.core_z.startValue = coreSD.GetZ(fLocalCS);
  }
}

void RdGlobalFit::setGlFitDataDirection(evt::Event& event)
{

  if (fGlobalFitData.arrivalDirection_phi.startValueFrom == "RD") {
    fGlobalFitData.arrivalDirection_phi.startValue = event.GetRecShower().GetRRecShower().GetAxis()
        .GetPhi(fLocalCS);
  }
  if (fGlobalFitData.arrivalDirection_phi.startValueFrom == "SD") {
    fGlobalFitData.arrivalDirection_phi.startValue = event.GetRecShower().GetSRecShower().GetAxis()
        .GetPhi(fLocalCS);
  }
  if (fGlobalFitData.arrivalDirection_phi.startValueFrom == "MC") {
    //This compensates the fact, that MC-axis points in oposit direction with respect to RD, SD.
    utl::Vector sAxis = -event.GetSimShower().GetDirection();
    fGlobalFitData.arrivalDirection_phi.startValue = sAxis.GetPhi(fLocalCS);
  }

  if (fGlobalFitData.arrivalDirection_theta.startValueFrom == "RD") {
    fGlobalFitData.arrivalDirection_theta.startValue =
        event.GetRecShower().GetRRecShower().GetAxis().GetTheta(fLocalCS);
  }
  if (fGlobalFitData.arrivalDirection_theta.startValueFrom == "SD") {
    fGlobalFitData.arrivalDirection_theta.startValue =
        event.GetRecShower().GetSRecShower().GetAxis().GetTheta(fLocalCS);
  }
  if (fGlobalFitData.arrivalDirection_theta.startValueFrom == "MC") {
    //This compensates the fact, that MC-axis points in oposit direction with respect to RD, SD.
    utl::Vector sAxis = -event.GetSimShower().GetDirection();
    fGlobalFitData.arrivalDirection_theta.startValue = sAxis.GetTheta(fLocalCS);
  }
}
int RdGlobalFit::calcNDF(ROOT::Minuit2::MnUserParameters upar, int nrSigSt)
{
  const std::vector<ROOT::Minuit2::MinuitParameter>& parameter = upar.Parameters();
  int dataPerStation = 0;
  if (fOrigFitConfig.fitTwoDLDF == true)
    ++dataPerStation;
  if (fOrigFitConfig.fitArrivalTime == true)
    ++dataPerStation;
  int NDF = dataPerStation * nrSigSt;  //Assuming Time and Amplitude are independent.
  for (std::vector<ROOT::Minuit2::MinuitParameter>::const_iterator parit = parameter.begin();
      parit != parameter.end(); ++parit) {

    if ((parit->GetName() == "core_x") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "core_y") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "core_z") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "arrivalDirection_theta") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "arrivalDirection_phi") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "Aplus") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "C2Theta") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "C1Theta") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "CTheta") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "C1") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "C2") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "sigmaPlus") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "gamma") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "b") && !parit->IsFixed())
      --NDF;
    if ((parit->GetName() == "t0") && !parit->IsFixed())
      --NDF;
  }
  std::cout << "calculated NDF = " << NDF << std::endl;
  return NDF;
}

VModule::ResultFlag RdGlobalFit::Finish()
{
  // Put any termination or cleanup code here.
  // This method is called once at the end of the run.
  INFO("RdGlobalFit::Finish()");
  return eSuccess;
}

}