RdHASLDFFitter.cc

Go to the documentation of this file.

#include "RdHASLDFFitter.h"
#include "LDFFitFunction.h"

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

// utilities
#include <utl/ErrorLogger.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Minou.h>

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

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

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

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

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


using namespace evt;
using namespace fwk;
using namespace std;
using namespace utl;


namespace RdHASLDFFitter {

  void
  RdHASLDFFitter::RadiationDensityCorrection(
    const double densityMax, const double densityMaxErr,
    double& radiationDensityCorrection, double& radiationDensityCorrectionErr) const
  {
    /* Takes density at shower maximum in kg / m3
    Returns correction factor for E_geo_rad / sin(alpha) ** 2*/

    const double densityAvg = fAvgDensity;  // * kg/m3;
    const double p0 = fp0;
    const double p1 = fp1;  // / (kg/m3);

    radiationDensityCorrection = pow(1 - p0 + p0 * exp(p1 * (densityMax - densityAvg)), -2);

    radiationDensityCorrectionErr =
      abs(-2 * p0 * p1 * exp(p1 * (densityMax - densityAvg)) * densityMaxErr) /
      pow(1 - p0 + p0 * exp(p1 * (densityMax - densityAvg)), 3);
  }


  double
  GetSineAlphaSquareErr(const Event& event)
  {
    // assumes that magnetic field is known perfectly. uncertainty of sin2(alpha). gaussian propagation without correlation
    const ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();
    const utl::Point& refShowerCore = showerrrec.GetReferenceCorePosition(event);
    const CoordinateSystemPtr refCoreCS = LocalCoordinateSystem::Create(refShowerCore);

    double zenith = 0;
    double azimuth = 0;
    double zenithErr = 0;
    double azimuthErr = 0;

    if (showerrrec.GetReferenceAxisFlag() == evt::ShowerRRecData::eUseSDReferenceAxis) {
      const evt::ShowerSRecData& sShowerRec = event.GetRecShower().GetSRecShower();
      const Vector sShowerAxis = sShowerRec.GetAxis();

      azimuthErr = sShowerRec.GetPhiError();
      zenithErr = sShowerRec.GetThetaError();
      zenith = sShowerAxis.GetTheta(refCoreCS);
      azimuth = sShowerAxis.GetPhi(refCoreCS);
    } else if (showerrrec.GetReferenceAxisFlag() == evt::ShowerRRecData::eUseRDReferenceAxis) {
      const Vector rShowerAxis = showerrrec.GetAxis();

      azimuthErr = showerrrec.GetAzimuthError();
      zenithErr = showerrrec.GetZenithError();
      zenith = rShowerAxis.GetTheta(refCoreCS);
      azimuth = rShowerAxis.GetPhi(refCoreCS);
    } else if (showerrrec.GetReferenceAxisFlag() == evt::ShowerRRecData::eUseMCReferenceAxis) {
      return 0;
    } else {
      ERROR("Selected ReferenceAxis (in RdEventInitalizer) is not supported! Return sineAlphaSquareErr = 0.");
      return 0;
    }

    const Vector magneticField = showerrrec.GetMagneticFieldVector();
    const Vector magneticFieldNormalized = magneticField / magneticField.GetMag();
    const double mx = magneticFieldNormalized.GetX(refCoreCS);
    const double my = magneticFieldNormalized.GetY(refCoreCS);
    const double mz = magneticFieldNormalized.GetZ(refCoreCS);

    // Derivertives are based on sin2alpha = (m x a) (m x a), with m=mag. vector and a=shower axis
    const double dsineAlphaSquareDzenith =
      2 * ((mz * cos(zenith) * cos(azimuth) + mx * sin(zenith)) * (-mx * cos(zenith) + mz * cos(azimuth) * sin(zenith)) +
      cos(zenith) * sin(zenith) * Sqr(my * cos(azimuth) - mx * sin(azimuth)) -
      (my * sin(zenith) + mz * cos(zenith) * sin(azimuth)) * (my * cos(zenith) - mz * sin(zenith) * sin(azimuth)));

    const double dsineAlphaSquareDazimuth =
      2 * sin(zenith) * (-my * cos(azimuth) + mx * sin(azimuth)) *
      (mz * cos(zenith) + sin(zenith) * (mx * cos(azimuth) + my * sin(azimuth)));

    const double sineAlphaSquareErr =
      sqrt(Sqr(dsineAlphaSquareDzenith * zenithErr) + Sqr(dsineAlphaSquareDazimuth * azimuthErr));
    return sineAlphaSquareErr;
  }


  void
  RdHASLDFFitter::ElectromagneticEnergy(
    const double sRad, const double sRadErr,
    double& electromagneticEnergyRec, double& electromagneticEnergyRecErr, const double zenith) const
  {
    /* Takes corrected geomagnetic radiation energy in eV
    Returns electromagnetic shower energy in eV */

    double s19 = fS19;
    const double gamma = fGamma;

    // See explaination for this correction in XML
    s19 = s19 / (1 + fSTEPFCCorrectionToSimulatedRadioEnergyScale);

    electromagneticEnergyRec = pow(sRad / s19, 1 / gamma) * 1e19 * eV;
    electromagneticEnergyRecErr = pow(sRad / s19, 1 / gamma - 1) * 1e19 * eV / (s19 * gamma) * sRadErr;

    if (fApplyParametricCorrectionToElectromagneticEnergyUncertainty) {
      // Parameterization explained in PhD Thesis Felix Schlüter Chap. 8.4.2.1.
      // Determined with an amplitude smearing of 5% (also see XML)
      const double correctionFactorEem =
        1.98593279e+03 * pow(Sqr(sin(zenith)) - Sqr(sin(65 * utl::deg)), 4) + 1.44615405e+00;
      electromagneticEnergyRecErr *= correctionFactorEem;
    }

  }


  bool
  RdHASLDFFitter::GetXmaxEstimator(const evt::Event& event, double& xmax)
    const
  {
    if (fXmaxType == "MC") {
      if (!event.HasSimShower()) {
        ERROR("Event has no sim shower but fXmaxType == simulation was choose.");
        return false;
      }
      const evt::ShowerSimData& simShower = event.GetSimShower();
      const auto& gh = simShower.GetGHParameters();
      xmax = gh.GetXMax();
    } else if (fXmaxType == "Param") {
      if (!event.GetRecShower().HasSRecShower()) {
        ERROR("Event has no SD rec shower. Can not estimate Xmax with utl::XmaxParam::Mean.");
        return false;
      }
      const evt::ShowerSRecData& sdShower = event.GetRecShower().GetSRecShower();
      const double crEnergy = sdShower.GetEnergy();
      xmax = 0.5 * (XmaxParam::Mean(crEnergy, 1, fHadronicInteractionModel)
                  + XmaxParam::Mean(crEnergy, 56, fHadronicInteractionModel));
    } else if (fXmaxType == "Average") {
      xmax = 750 * g / cm2;  // ~ 10 EeV, 50% proton, 50% iron composition
    } else {
      ERROR("Wrong fXmaxType choosen!");
      return false;
    }

    return true;
  }


  VModule::ResultFlag
  RdHASLDFFitter::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdHASLDFFitter");

    // Read in the configurations of the xml file
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    topBranch.GetChild("MinimumNumberOfStations").GetData(fMinimumNumberOfStations);
    topBranch.GetChild("FitDistanceToXmax").GetData(fFitDistanceToXmax);
    topBranch.GetChild("FitCore").GetData(fFitCore);
    topBranch.GetChild("AddLowestSignalAsError").GetData(fAddLowestSignalAsError);
    topBranch.GetChild("AddRelativeMaxSignalError").GetData(fAddRelativeMaxSignalError);

    topBranch.GetChild("XmaxEstimator").GetData(fXmaxType);

    std::string hadronicInteractionModule =
      topBranch.GetChild("HadronicInteractionModel").Get<string>();
    if (hadronicInteractionModule == "Sibyll2.3d")
      fHadronicInteractionModel = utl::HadronicInteractionModel::eSibyll23d;
    else if (hadronicInteractionModule == "EPOS-LHC")
      fHadronicInteractionModel = utl::HadronicInteractionModel::eEPOSLHC;
    else if (hadronicInteractionModule == "QGSJETII-04")
      fHadronicInteractionModel = utl::HadronicInteractionModel::eQGSJETII04;
    else {
      throw InvalidConfigurationException(
        "You used an invalid HadronicInteractionModel in the XML config.");
      return eFailure;
    }

    topBranch.GetChild("FitDirection").GetData(fFitDirection);
    topBranch.GetChild("UseSaturatedStations").GetData(fUseSaturatedStations);
    topBranch.GetChild("UseParametrizationToDisentanglePolarisation").GetData(
      fUseParametrizationToDisentanglePolarisation);
    topBranch.GetChild("UseSoftExponent").GetData(fUseSoftExponent);

    topBranch.GetChild("S19").GetData(fS19);
    topBranch.GetChild("Gamma").GetData(fGamma);
    topBranch.GetChild("AvgDensity").GetData(fAvgDensity);
    topBranch.GetChild("p0").GetData(fp0);
    topBranch.GetChild("p1").GetData(fp1);
    topBranch.GetChild("STEPFCCorrectionToSimulatedRadioEnergyScale")
      .GetData(fSTEPFCCorrectionToSimulatedRadioEnergyScale);

    topBranch.GetChild("ApplyParametricCorrectionToElectromagneticEnergyUncertainty")
      .GetData(fApplyParametricCorrectionToElectromagneticEnergyUncertainty);

    ostringstream info;
    info << "\n"
            "\tMinimum number of required stations: " << fMinimumNumberOfStations << "\n"
            "\tUse saturated stations: " << fUseSaturatedStations << "\n"
            "\tFit distance to xmax / core: (" << fFitDistanceToXmax << "/" << fFitCore << ")\n"
            "\tXmax estimator used to determine inital fit values: " << fXmaxType << "\n";
    if (fXmaxType == "Param")
      info << "\tUsed high-energy hadronic interaction model : " << hadronicInteractionModule << "\n";
    info << "\tDirection used in fit: " << fFitDirection << "\n"
            "\tUse parametrization to determine geomagnetic energy fluence: "
            << fUseParametrizationToDisentanglePolarisation << "\n"
            "\tUse ";
    if (fUseSoftExponent)
      info << "\"soft\" (default) ";
    else
      info << "\"hard\" (not default) ";
    info << "exponent in geomagnetic ldf" << "\n"
         << "\tUse correction for uncertainty: " << fApplyParametricCorrectionToElectromagneticEnergyUncertainty << "\n"
         << "\tApply STEPFC energy scale correction factor of: " << fSTEPFCCorrectionToSimulatedRadioEnergyScale << "\n";
    INFOFinal(info);

    return eSuccess;
  }


  VModule::ResultFlag
  RdHASLDFFitter::Run(evt::Event& event)
  {
    if (!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }
    revt::REvent& rEvent = event.GetREvent();

    // initialize data structures
    if (!event.HasRecShower())
      event.MakeRecShower();

    if (!event.GetRecShower().HasRRecShower())
      event.GetRecShower().MakeRRecShower();
    ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();

    // get reference coordinate system of detector (usually PampaAmarilla)
    const CoordinateSystemPtr referenceCS =
      det::Detector::GetInstance().GetReferenceCoordinateSystem();

    // takes the reference geometry set in event initializer
    const utl::Point& refShowerCore = showerrrec.GetReferenceCorePosition(event);
    utl::Vector refShowerAxis;

    if (fFitDirection == "Reference") {
      refShowerAxis = showerrrec.GetReferenceAxis(event);
    } else if (fFitDirection == "Rd") {
      if (!showerrrec.HasAxis()) {
        INFO("Rd axis requested for reconstruction. Not available. Skip event ...");
        return eContinueLoop;
      }
      refShowerAxis = showerrrec.GetAxis();
    }

    const CoordinateSystemPtr refCoreCS = LocalCoordinateSystem::Create(refShowerCore);
    const double zenith = refShowerAxis.GetTheta(refCoreCS);
    const double azimuth = refShowerAxis.GetPhi(refCoreCS);

    ostringstream info;
    info << "Use geometry:\n"
            "\tzenith = " << zenith / utl::deg << " deg, azimuth = " << azimuth / utl::deg << " deg\n"
            "\t" << "core in ref cs (PA): (" << refShowerCore.GetX(referenceCS) / meter << ", "
                                             << refShowerCore.GetY(referenceCS) / meter << ", "
                                             << refShowerCore.GetZ(referenceCS) / meter << ')';
    INFOIntermediate(info);

    const Vector magneticField = showerrrec.GetMagneticFieldVector();
    utl::RadioGeometryUtilities transformation = RadioGeometryUtilities(refShowerAxis, refCoreCS, magneticField);
    const double sineAlpha = RadioGeometryUtilities::GetLorentzVector(refShowerAxis, magneticField).GetMag();

    // get detector instances
    const det::Detector& detector = det::Detector::GetInstance();
    const rdet::RDetector& rDetector = detector.GetRDetector();

    vector<StationFitData> vStationFitData;
    unsigned int numberOfSignalStations = 0;

    double addErrMin = 0;
    if (fAddLowestSignalAsError) {
      addErrMin = 1e25;
      // iteration over stations
      for (const auto& station : rEvent.SignalStationsRange()) {
        const auto& sRec = station.GetRecData();

        // skip saturated stations if it is configured
        if (!fUseSaturatedStations && station.IsSaturated())
          continue;

        if (addErrMin > sRec.GetParameter(revt::eSignalEnergyFluenceVxB))
            addErrMin = sRec.GetParameter(revt::eSignalEnergyFluenceVxB);
      }
    }

    double addErrMax = 0;
    if (fAddRelativeMaxSignalError) {
      if (!event.HasSimShower()) {
        throw InvalidConfigurationException(
            "fAddRelativeMaxSignalError is true, requieres SimShower.");
      }
      const evt::ShowerSimData &simShower = event.GetSimShower();
      const double dmaxMc = simShower.GetDistanceOfShowerMaximum();
      const double eemMc = simShower.GetElectromagneticEnergy();
      const double sineAlphaMC = RadioGeometryUtilities::GetLorentzVector(
        -simShower.GetDirection(), magneticField).GetMag();
      // simple parameterisation of f_geo_max
      const double fGeoMax = (6.70443081e+07 * meter / dmaxMc + 1.07044401e+02) * Sqr(eemMc / EeV * sineAlphaMC);
      addErrMax = fAddRelativeMaxSignalError * fGeoMax;
    }

    // iteration over stations - Sofar this ldf module can not make use of NoSignal stations
    for (auto& station : rEvent.SignalStationsRange()) {
      auto& sRec = station.GetRecData();

      // skip saturated stations if it is configured
      if (!fUseSaturatedStations && station.IsSaturated())
        continue;

      if (fAddLowestSignalAsError) {
        sRec.SetParameterError(revt::eSignalEnergyFluenceVxB,
                               sqrt(Sqr(sRec.GetParameterError(revt::eSignalEnergyFluenceVxB)) +
                                    Sqr(addErrMin)));
      }

      if (fAddRelativeMaxSignalError) {
        sRec.SetParameterError(revt::eSignalEnergyFluenceVxB,
                               sqrt(Sqr(sRec.GetParameterError(revt::eSignalEnergyFluenceVxB)) +
                                    Sqr(addErrMax)));
      }

      StationFitData stationFitData;
      stationFitData.fId = station.GetId();
      stationFitData.fPosition = rDetector.GetStation(station).GetPosition();
      stationFitData.fEnergyFluenceVxB = sRec.GetParameter(revt::eSignalEnergyFluenceVxB);
      stationFitData.fEnergyFluenceVxBErr = sRec.GetParameterError(revt::eSignalEnergyFluenceVxB);
      stationFitData.fEnergyFluenceVxVxB = sRec.GetParameter(revt::eSignalEnergyFluenceVxVxB);
      stationFitData.fEnergyFluenceVxVxBErr = sRec.GetParameterError(revt::eSignalEnergyFluenceVxVxB);
      stationFitData.fRejected = false;
      vStationFitData.push_back(stationFitData);

      ++numberOfSignalStations;
    }

    info.str("");
    info << "Number of Signal Stations: " << numberOfSignalStations;
    INFOIntermediate(info);

    if (numberOfSignalStations < fMinimumNumberOfStations) {
      WARNING("Number of signal stations after rejection " +
              to_string(fMinimumNumberOfStations) + " -> doing nothing!");
      showerrrec.SetParameter(revt::eHASLDFHasFit, false);
      return eSuccess;
    }

    // Estimate Xmax to estimate start value for dmax
    double showerXmaxGuess = 0;
    if (!GetXmaxEstimator(event, showerXmaxGuess) || showerXmaxGuess / (g/cm2) <= 1) {
      info.str("");
      info << "Could not get xmax estimator \"" << fXmaxType << "\": "
           << showerXmaxGuess / (g / cm2);
      WARNING(info);
      showerrrec.SetParameter(revt::eHASLDFHasFit, false);
      return eSuccess;
    }

    // initiate atm + tables (which model is used is decieded in the config)
    const atm::Atmosphere& theAtm = detector.GetAtmosphere();
    theAtm.InitSlantProfileModel(refShowerCore, refShowerAxis, 10 * g / cm2);
    //theAtm.InitSlantProfileModel(utl::Point(0., 0., 0., referenceCS), refShowerAxis, 10.*g/cm2);

    const atm::ProfileResult& distanceFromDepth = theAtm.EvaluateDistanceVsSlantDepth();
    const atm::ProfileResult& densityFromHeight = theAtm.EvaluateDensityVsHeight();
    const atm::ProfileResult& heightFromDistance = theAtm.EvaluateHeightVsDistance();  // needed within fit
    const atm::ProfileResult& refracFromHeight = theAtm.EvaluateRefractionIndexVsHeight();  // needed within fit

    // shower maximum
    const double distanceToXmaxGuess = abs(distanceFromDepth.Y(showerXmaxGuess));  // yep, the result is negative, convention thing...

    const utl::Point showerMax = refShowerCore + refShowerAxis * distanceToXmaxGuess;  // -1 * axis ???
    info.str("");
    info << "Position of shower maximum (estimation with \"" << fXmaxType << "\": "
         << showerXmaxGuess / (g/cm2) << ") (for PA): (" << showerMax.GetX(refCoreCS) / meter
         << ", " << showerMax.GetY(refCoreCS) / meter
         << ", " << showerMax.GetZ(refCoreCS) / meter << ')';
    INFODebug(info);

    // set shower related quantities use in the fit
    ShowerFitData fitShowerData;
    fitShowerData.fSineAlpha = sineAlpha;
    fitShowerData.fZenith = zenith;
    fitShowerData.fDistanceTo750 = abs(distanceFromDepth.Y(750 * g / cm2));
    fitShowerData.fRefShowerAxis = refShowerAxis;
    fitShowerData.fMagneticField = magneticField;
    fitShowerData.fRefShowerCore = refShowerCore;
    fitShowerData.fDensityFromHeight = densityFromHeight;
    fitShowerData.fHeightFromDistance = heightFromDistance;
    fitShowerData.fRefracFromHeight = refracFromHeight;

    // set config
    FitConfig fitConfig;
    fitConfig.fUseParam = fUseParametrizationToDisentanglePolarisation;
    fitConfig.fUseSoftExponent = fUseSoftExponent;

    // determine if core fit is feasible
    const bool fitDXmax = (fFitDistanceToXmax && numberOfSignalStations >= 3);
    const bool fitCore = (fFitCore && !fFitDistanceToXmax && numberOfSignalStations >= 3) ||
      (fFitCore && numberOfSignalStations >= 4);

    // get start shape parameter from parametrisation
    double crEnergy;
    if (event.GetRecShower().HasSRecShower()) {
      crEnergy = event.GetRecShower().GetSRecShower().GetEnergy();
    } else {
      INFOFinal("No SD rec shower information available. "
                "Use a cosmic ray energy 1 EeV to determine start value for radiation energy ");
      crEnergy = 1 * EeV;
    }
    const double startEgeo = 26.86e6 * eV * pow(crEnergy / (1e18 * eV), 1.989) * Sqr(sineAlpha);

    info.str("");
    info << "Attempt LDF fit: Egeo = " << startEgeo << " eV, dmax = " << distanceToXmaxGuess << " m.";
    INFOFinal(info);

    // Initalize parameters (WARNING: Take care of correct ordering)
    vector<utl::Minou::ParameterDef> parameters;

    // name, value, step, min, max, fixed
    parameters.emplace_back("E_geo", startEgeo, 10, 1e3, 1e13, false);

    // distance to xmax (fixed)
    // boundaries are strange because of table convention (see -1 in fit function)
    const double dxmaxMin = max(-heightFromDistance.MaxX(), 0.7 * distanceToXmaxGuess / meter);
    const double dxmaxMax = min(-heightFromDistance.MinX(), 1.5 * distanceToXmaxGuess / meter);
    parameters.emplace_back("distance_xmax_geometric", distanceToXmaxGuess / meter, 100, dxmaxMin, dxmaxMax, true);

    // fixed core
    parameters.emplace_back("core_x", 0, 1, -1000, 1000, true);
    parameters.emplace_back("core_y", 0, 1, -1000, 1000, true);

    // init fitfunction
    const double egeoCorrectionFactor = 1;
    LDFFitFunctionGaussSigmoid fitFunction(
      vStationFitData, fitShowerData, parameters, fitConfig, egeoCorrectionFactor
    );

    // to store fit results
    vector<double> tmpFitResult;
    double chi2 = 1e99;
    int ndf = 0;
    vector<pair<double, double>> bestParametersAndErrors;
    {
      // fit
      Minou::Minimizer<LDFFitFunction> min(fitFunction);
      INFOIntermediate("Inital fit: only Egeo is free");
      if (fInfoLevel >= eInfoIntermediate)
        min.SetVerbose(true);

      if (min.Minimize(50000)) {
        INFOFinal("Fit falied ...");
        showerrrec.SetParameter(revt::eHASLDFHasFit, false);
        ++fFailedFits;
        return eSuccess;
      }

      // store station data
      tmpFitResult = min.GetParameters();
      fitFunction.GetPrediction(tmpFitResult);
      chi2 = fitFunction.GetChi2(tmpFitResult);
      ndf = fitFunction.GetNDF();
      bestParametersAndErrors = min.GetParametersAndErrors();
    }

    if (fitCore) {
      fitFunction.SetParameterDefValues(tmpFitResult);

      // A, Dxmax coreX, coreY
      const vector<int> fixed({0, 1, 0, 0});
      fitFunction.SetParameterDefFixed(fixed);

      // fit
      Minou::Minimizer<LDFFitFunction> min(fitFunction);
      INFOIntermediate("Egeo, coreX, coreY are free");

      if (fInfoLevel >= eInfoIntermediate)
        min.SetVerbose(true);

      if (min.Minimize(50000)) {
        INFOFinal("Fit falied ...");
        showerrrec.SetParameter(revt::eHASLDFHasFit, false);
        ++fFailedFits;
        return eSuccess;
      }

      // store station data
      tmpFitResult = min.GetParameters();
      fitFunction.GetPrediction(min.GetParameters());
      chi2 = fitFunction.GetChi2(min.GetParameters());
      ndf = fitFunction.GetNDF();
      bestParametersAndErrors = min.GetParametersAndErrors();
    }

    if (fitDXmax) {
      fitFunction.SetParameterDefValues(tmpFitResult);

      // A, Dxmax coreX, coreY
      const vector<int> fixed({0, 0, 1, 1});
      fitFunction.SetParameterDefFixed(fixed);

      // fit
      Minou::Minimizer<LDFFitFunction> min(fitFunction);
      INFOIntermediate("Egeo, dmax are free");

      if (fInfoLevel >= eInfoIntermediate)
        min.SetVerbose(true);

      if (min.Minimize(50000)) {
        INFOFinal("Fit falied ...");
        showerrrec.SetParameter(revt::eHASLDFHasFit, false);
        ++fFailedFits;
        return eSuccess;
      }

      // store station data
      tmpFitResult = min.GetParameters();
      fitFunction.GetPrediction(min.GetParameters());
      chi2 = fitFunction.GetChi2(min.GetParameters());
      ndf = fitFunction.GetNDF();
      bestParametersAndErrors = min.GetParametersAndErrors();
    }

    if (fitCore) {
      fitFunction.SetParameterDefValues(tmpFitResult);

      // A, Dxmax coreX, coreY
      const vector<int> fixed({0, !fitDXmax, 0, 0});
      fitFunction.SetParameterDefFixed(fixed);

      // fit
      Minou::Minimizer<LDFFitFunction> min(fitFunction);
      if (fitDXmax) {
        INFOIntermediate("Egeo, coreX, coreY, dmax are free");
      } else {
        INFOIntermediate("Egeo, coreX, coreY are free");
      }

      if (fInfoLevel >= eInfoIntermediate)
        min.SetVerbose(true);

      if (min.Minimize(50000)) {
        INFOFinal("Fit falied");
        showerrrec.SetParameter(revt::eHASLDFHasFit, false);
        ++fFailedFits;
        return eSuccess;
      }

      // store station data
      tmpFitResult = min.GetParameters();
      fitFunction.GetPrediction(min.GetParameters());
      chi2 = fitFunction.GetChi2(min.GetParameters());
      ndf = fitFunction.GetNDF();
      bestParametersAndErrors = min.GetParametersAndErrors();
    }

    // ------ get fit results ---------
    const double egeo = bestParametersAndErrors.at(0).first;
    const double egeoErr = bestParametersAndErrors.at(0).second;

    const double distanceToXmaxFit = bestParametersAndErrors.at(1).first;
    double distanceToXmaxFitErr;
    if (fitDXmax) {
      distanceToXmaxFitErr = bestParametersAndErrors.at(1).second;
    } else {
      // error when not fitting dxmax but take the estimated from the param. Xmax (is actually not gaussian!)
      distanceToXmaxFitErr = 0.06 * distanceToXmaxFit;
    }

    const double coreX = bestParametersAndErrors.at(2).first;
    const double coreXerr = bestParametersAndErrors.at(2).second;

    const double coreY = bestParametersAndErrors.at(3).first;
    const double coreYerr = bestParametersAndErrors.at(3).second;

    double densityAtXmaxFit = 0;
    double densityAtXmaxFitErr = 0;
    try {
      const double dxmax_sysdown = distanceToXmaxFit - distanceToXmaxFitErr;
      const double dxmax_sysup = distanceToXmaxFit + distanceToXmaxFitErr;
      // here the "-" signs are due to table convention
      densityAtXmaxFit = densityFromHeight.Y(heightFromDistance.Y(-distanceToXmaxFit));
      densityAtXmaxFitErr =
        (densityFromHeight.Y(heightFromDistance.Y(-dxmax_sysdown)) -
         densityFromHeight.Y(heightFromDistance.Y(-dxmax_sysup))) / 2;
    } catch (const std::exception& e) {
      INFOFinal("No proper error propagation possible, fit is discarded."
                " distanceToXmaxFit +/- distanceToXmaxFitErr out of table range for heightFromDistance");
      showerrrec.SetParameter(revt::eHASLDFHasFit, false);
      ++fFailedFits;
      return eSuccess;
    }

    info.str("");
    info << "LDF fit succesfull:\n"
            "\tEgeo = " << egeo/MeV << " +- " <<  egeoErr / MeV << " MeV (" << egeoErr / egeo * 100 << " %)\n"
            "\tdmax = " << distanceToXmaxFit / km << " +- " <<  distanceToXmaxFitErr / km << " km "
            "(" << distanceToXmaxFitErr / distanceToXmaxFit * 100 << " %)\n"
            "\tCoreX = " << coreX / m << " +- " <<  coreXerr / m << " m, CoreY = "
            << coreY / m << " +- " <<  coreYerr / m << " m\n"
            "\tchi2 / ndf = " << chi2 << " / " << ndf << " = " << chi2 / ndf;
    INFOFinal(info);

    showerrrec.SetParameter(revt::eHASLDFHasFit, true);
    showerrrec.SetParameter(revt::eHASLDFChi2, chi2);
    showerrrec.SetParameter(revt::eHASLDFNDF, ndf);

    for (const auto& s : vStationFitData) {
      revt::Station& station = rEvent.GetStation(s.fId);
      revt::StationRecData& sRec = station.GetRecData();

      sRec.SetParameter(revt::eEarlyLateCorrectionFactor, s.fCEarlyLate);
      sRec.SetParameter(revt::eChargeExcessFraction, s.fCeFraction);
      sRec.SetParameter(revt::ePredictedEnergyFluenceVxB, s.fEnergyFluenceVxBPredict);
      sRec.SetParameter(revt::eGeomagneticEnergyFluence, s.fEnergyFluenceGeomagnetic);
      sRec.SetParameterError(revt::eGeomagneticEnergyFluence, s.fEnergyFluenceGeomagneticErr);
    }

    // store shower data in Auger units
    showerrrec.SetParameter(revt::eGeomagneticRadiationEnergy, egeo); // is already in eV
    showerrrec.SetParameterError(revt::eGeomagneticRadiationEnergy, egeoErr); // is already in eV

    showerrrec.SetParameter(revt::eHASLDFDensityAtXmax, densityAtXmaxFit / (kg/m3));
    showerrrec.SetParameterError(revt::eHASLDFDensityAtXmax, densityAtXmaxFitErr / (kg/m3));
    showerrrec.SetParameter(revt::eHASLDFDistanceToXmaxGeometric, distanceToXmaxFit / meter);  // in meter
    showerrrec.SetParameterError(revt::eHASLDFDistanceToXmaxGeometric, distanceToXmaxFitErr / meter);  // in meter

    showerrrec.SetParameter(revt::eHASLDFCoreX, coreX);
    showerrrec.SetParameter(revt::eHASLDFCoreY, coreY);
    showerrrec.SetParameterError(revt::eHASLDFCoreX, coreXerr);
    showerrrec.SetParameterError(revt::eHASLDFCoreY, coreYerr);

    // original core is origin in vB system
    utl::Point radioCore = transformation.GetVectorFromShowerPlaneVxB(coreX, coreY);
    showerrrec.SetParameter(revt::eCoreX, radioCore.GetX(referenceCS));
    showerrrec.SetParameter(revt::eCoreY, radioCore.GetY(referenceCS));
    showerrrec.SetParameter(revt::eCoreZ, radioCore.GetZ(referenceCS));

    utl::RadioGeometryUtilities transformationRadioCore =
      RadioGeometryUtilities(refShowerAxis, LocalCoordinateSystem::Create(radioCore), magneticField);

    // iteration over stations - store position in vB,vvB with final radio core
    // we use the station iterator here to save the position in the vxB frame for all stations
    for (auto& station : rEvent.StationsRange()) {
      revt::StationRecData& sRec = station.GetRecData();
      const Point sPos = rDetector.GetStation(station).GetPosition();

      // due to the choice of coordinate system the station position will be evaluated relative to the core position
      double xvxB = 0;
      double yvxB = 0;
      double zvxB = 0;
      transformationRadioCore.GetVectorInShowerPlaneVxB(xvxB, yvxB, zvxB, sPos);

      sRec.SetParameter(revt::eLDFFitStationPositionVxB, xvxB);
      sRec.SetParameter(revt::eLDFFitStationPositionVxVxB, yvxB);
      sRec.SetParameter(revt::eLDFFitStationPositionV, zvxB);
    }

    // get energy reconstruction
    double densityCorrection = 0;
    double densityCorrectionErr = 0;
    RadiationDensityCorrection(densityAtXmaxFit, densityAtXmaxFitErr, densityCorrection, densityCorrectionErr);
    showerrrec.SetParameter(revt::eDensityCorrectionFactor, densityCorrection);
    showerrrec.SetParameterError(revt::eDensityCorrectionFactor, densityCorrectionErr);

    double radioEnergyEstimator = egeo / Sqr(sineAlpha) * densityCorrection;

    const double sineAlphaSquareErr = GetSineAlphaSquareErr(event);
    const double radioEnergyEstimatorErr =
      radioEnergyEstimator *
      sqrt(Sqr(egeoErr / egeo) +
           Sqr(densityCorrectionErr / densityCorrection) +
           Sqr(sineAlphaSquareErr / Sqr(sineAlpha)));

    showerrrec.SetParameter(revt::eCorrectedGeomagneticRadiationEnergy, radioEnergyEstimator);
    showerrrec.SetParameterError(revt::eCorrectedGeomagneticRadiationEnergy, radioEnergyEstimatorErr);

    double electromagneticEnergyRec = 0;
    double electromagneticEnergyRecErr = 0;
    ElectromagneticEnergy(
      radioEnergyEstimator, radioEnergyEstimatorErr,
      electromagneticEnergyRec, electromagneticEnergyRecErr, zenith);
    showerrrec.SetParameter(revt::eReconstructedElectromagneticEnergy, electromagneticEnergyRec);
    showerrrec.SetParameterError(revt::eReconstructedElectromagneticEnergy, electromagneticEnergyRecErr);

    /* Calibration from electromagnetic to absolut shower/cosmic ray energy. See PhD Thesis from Felix Schüter Chap. 8.4.3.
       Calibration determined with QGSJetII-04 and 4 primaries: p, He, N, Fe. Results with Sibyll2.3d were very similar */
    const double cosmicRayEnergy = electromagneticEnergyRec * (1.14263 - 0.03279 * log10(electromagneticEnergyRec / (10 * utl::EeV)));
    const double cosmicRayEnergyErr = electromagneticEnergyRecErr * (cosmicRayEnergy / electromagneticEnergyRec);
    showerrrec.SetParameter(revt::eCosmicRayEnergy, cosmicRayEnergy);
    showerrrec.SetParameterError(revt::eCosmicRayEnergy, cosmicRayEnergyErr);

    // fix or parametrized
    const double r0 = fitFunction.CalculateR0();
    const double sig = fitFunction.GetSig(distanceToXmaxFit);
    const double p = fitFunction.GetP(distanceToXmaxFit);
    const double arel = fitFunction.GetArel(distanceToXmaxFit);
    const double r02 = fitFunction.GetR02(distanceToXmaxFit);

    showerrrec.SetParameter(revt::eHASLDFNorm, fitFunction.GetNormalization());

    showerrrec.SetParameter(revt::eHASLDFr0, r0);
    showerrrec.SetParameter(revt::eHASLDFsig, sig);
    showerrrec.SetParameter(revt::eHASLDFp, p);
    showerrrec.SetParameter(revt::eHASLDFarel, arel);
    showerrrec.SetParameter(revt::eHASLDFslope, 5);
    showerrrec.SetParameter(revt::eHASLDFr02, r02);

    // If it reaches this point: successful fit
    ++fSuccessfulFits;
    return eSuccess;
  }


  VModule::ResultFlag
  RdHASLDFFitter::Finish()
  {
    ostringstream info;
    const double total = fSuccessfulFits + fFailedFits;  // double to prevent integer devision
    if (total) {
      info << "Number of successful and failed fits: " << fSuccessfulFits << " / " << fFailedFits
           << " (" << fSuccessfulFits / total * 100 << "%, " << fFailedFits / total * 100 << "%)";
    } else
      info << "No radio has ldf fit was performed!";

    INFOFinal(info);
    return eSuccess;
  }

}