RdLDFFitter.cc

Go to the documentation of this file.

#include "RdLDFFitter.h"

#include <cmath>
#include <sstream>
#include <fstream>
#include <iostream>
#include <iomanip>  // std::setprecision
#include <utility>  // for pair
#include <string>
#include <vector>

#include <boost/filesystem/operations.hpp>

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

#include <utl/config.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/RectangleFilter.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/RadioGeometryUtilities.h>
#include <utl/Probability.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 <det/Detector.h>
#include <rdet/RDetector.h>

// include root stuff
#include "TGraphErrors.h"
#include "TF1.h"
#include "TH2D.h" // for plotting of Likelihood scan
#include "TCanvas.h"
#include "TLegend.h"
#include "TEllipse.h"
#include "TMath.h"
#include "TLatex.h"
#include "TLine.h"
#include "TStyle.h"
#include "TAxis.h"


#include "Minuit2/Minuit2Minimizer.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MinimumParameters.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnPlot.h"
#include "Minuit2/MnUserCovariance.h"
//#include "TFitterMinuit.h"

#include "LikelihoodFunction.h"

using namespace revt;
using namespace fwk;
//using namespace utl;

#define OUT(x) if ((x) <= fInfoLevel) std::cerr << "  "
#define INFO2(x,y) if ((x) <= fInfoLevel) INFO(y)

namespace RdLDFFitter {

RdLDFFitter::RdLDFFitter() :
    fFitConfig(),
    fInfoLevel(eNone),
    fSaveContour(false),
    fPerformScan(false),
    fPlotSigmaContours(true),
    fNStepsX(11),
    fNStepsY(11),
    fStepWidthX(10*utl::meter),
    fStepWidthY(10*utl::meter),
    fScanCenter(),
    fMCUncertaintyOnEFieldVector(0.),
    fPlotStyle(0),
    fREvent(NULL),
    fConstEvent(NULL),
    fEvent(NULL),
    fShowerAxis(),
    fMagneticFieldAxis()
{
}

RdLDFFitter::~RdLDFFitter()
{
}

VModule::ResultFlag RdLDFFitter::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("RdLDFFitter::Init()");

  // Read in the configurations of the xml file
  utl::Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdLDFFitter");

  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
  topBranch.GetChild("SaveContour").GetData(fSaveContour);
  topBranch.GetChild("performScan").GetData(fPerformScan);
  topBranch.GetChild("nStepsX").GetData(fNStepsX);
  topBranch.GetChild("nStepsY").GetData(fNStepsY);
  topBranch.GetChild("stepWidthX").GetData(fStepWidthX);
  topBranch.GetChild("stepWidthY").GetData(fStepWidthY);
  topBranch.GetChild("outputFolder").GetData(fOutputFolder);
  topBranch.GetChild("MCUncertaintyOnEFieldVector").GetData(fMCUncertaintyOnEFieldVector);

  std::string temp;
  topBranch.GetChild("scanCenter").GetData(temp);
  if(temp == "Radio") 
    { fScanCenter = eRadio; }
  if(temp == "SD") 
    { fScanCenter = eSD; }
  if(temp == "FD") 
    { fScanCenter = eFD; }
  if(temp == "MC") 
    { fScanCenter = eMC; }


//    int tmp = 0;
//    topBranch.GetChild("LDFModel").GetData(tmp);
//    fFitConfig.LDFModel = tmp;
  utl::Branch bFitConfig = topBranch.GetChild("FitConfig");
  bFitConfig.GetChild("LDFModel").GetData(fFitConfig.LDFModel);
  bFitConfig.GetChild("chargeExcessStrength").GetData(fFitConfig.chargeExcessStrength);
  bFitConfig.GetChild("fitChargeExcess").GetData(fFitConfig.fitChargeExcess);
  bFitConfig.GetChild("useChargeExcessCorrectionInLDFFit").GetData(
      fFitConfig.useChargeExcessCorrectionInLDFFit);
  bFitConfig.GetChild("usePolarisation").GetData(fFitConfig.usePolarisation);
  bFitConfig.GetChild("useSDCoreToImproveRadioCore").GetData(
      fFitConfig.useSDCoreToImproveRadioCore);
  bFitConfig.GetChild("useStationsWithoutSignal").GetData(fFitConfig.useStationsWithoutSignal);
  bFitConfig.GetChild("useScintillator").GetData(fFitConfig.useScintillator);

  return eSuccess;
}

VModule::ResultFlag RdLDFFitter::Run(evt::Event& event)
{

  // Check if there are events at all
  if (!event.HasREvent()) {
    WARNING("No radio event found!");
    return eSuccess;
  }
  std::stringstream info;
  fConstEvent = & event;   // setting pointers to event structure to access it from other functions
  fEvent = & event;   // setting pointers to event structure to access it from other functions
  fREvent = &event.GetREvent();

  revt::REvent& rEvent = event.GetREvent();
  const unsigned int runNumber = fREvent->GetHeader().GetRunNumber();
  const unsigned int eventId = fREvent->GetHeader().GetId();

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

  if(!(showerrrec.HasParameter(eRecStage) and showerrrec.GetParameter(eRecStage) >= evt::ShowerRRecData::kPlaneFit3d)) {
    WARNING("The reconstruction stage is not set or smaller that kPlaneFit3d. LDF Fit can not be performed.");
    return eSuccess;
  }

  //Get Radio detector for station loop
  const det::Detector& detector = det::Detector::GetInstance();
  const rdet::RDetector& rDetector = detector.GetRDetector();

  fShowerAxis = showerrrec.GetAxis();
  utl::Point coordinateOrigin = showerrrec.GetCoordinateOrigin();
  fLocalCS = LocalCoordinateSystem::Create(coordinateOrigin);

  fMagneticFieldAxis = event.GetRecShower().GetRRecShower().GetMagneticFieldVector();
  // use magnetic field from simulations if present
  if(event.HasSimShower()) {
    if(fInfoLevel > eFinal) {
      INFO("Using magnetic field from simulated air shower.");
    }
    const double magFieldAzimuth = event.GetSimShower().GetMagneticFieldAzimuth();
    const double magFieldZenith = event.GetSimShower().GetMagneticFieldZenith();
    info.str("");
    info << "Difference in MC and Offline magnetic field:\n"
         << fMagneticFieldAxis.GetPhi(fLocalCS)/utl::deg << " vs. " << magFieldAzimuth / utl::deg << "\n"
         << fMagneticFieldAxis.GetTheta(fLocalCS)/utl::deg << " vs. " << magFieldZenith / utl::deg;
    INFO2(eDebug, info.str());
    fMagneticFieldAxis = utl::Vector(1., magFieldZenith, magFieldAzimuth, fLocalCS,
                                     utl::Vector::kSpherical);
    fShowerAxis = -1. * event.GetSimShower().GetDirection();
  }

  // set EventFitData
  EventFitData eventFitData;
  eventFitData.fShowerAxis = fShowerAxis;
  eventFitData.fLocalCS = fLocalCS;

  if (event.GetRecShower().HasSRecShower()) {
    evt::ShowerSRecData& showersrec = event.GetRecShower().GetSRecShower();  // get reconstructed SD shower for SD core position
    eventFitData.SDCore = showersrec.GetCorePosition();
    eventFitData.SDCoreXError = showersrec.GetCoreError().GetX(fLocalCS);
    eventFitData.SDCoreYError = showersrec.GetCoreError().GetY(fLocalCS);
    eventFitData.SDCoreXYCorrelation = showersrec.GetCorrelationXY();
  }

  if (fFitConfig.useSDCoreToImproveRadioCore) {
    if (not event.GetRecShower().HasSRecShower()) {
      WARNING(
          "RdLDFFitter: No reconstructed SD shower found. SD core position can't be used to improve radio core. This option will be deactivated for this event.");
      fFitConfig.useSDCoreToImproveRadioCore = false;
    }
  }

  std::vector<StationFitData> vStationFitData;

  INFO2(eDebug, "starting station loop");

  // loop through stations and fill all relevant parameters into the StationFitData struct
  info.str("");
  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    Station& station = *sIt;
    StationFitData tmpStationFitData;
    const StationRecData& sRec = station.GetRecData();
    const utl::Point sPos = rDetector.GetStation(station).GetPosition();
    tmpStationFitData.antennaPosition = sPos;
    tmpStationFitData.signal = sRec.GetParameter(revt::eSignal);
    tmpStationFitData.signalError = sRec.GetParameterError(revt::eSignal);
    utl::Vector efield = utl::Vector(sRec.GetParameter(revt::eEFieldVectorEW),
                                     sRec.GetParameter(revt::eEFieldVectorNS),
                                     sRec.GetParameter(revt::eEFieldVectorV), fLocalCS);
    if (fMCUncertaintyOnEFieldVector) {
      efield = utl::Probability::GetInstance().GetRandomFisher(
          efield, 1. / (fMCUncertaintyOnEFieldVector * fMCUncertaintyOnEFieldVector));
      info << "smeer out efield by " << fMCUncertaintyOnEFieldVector / utl::deg << " degree, ";
    }
    tmpStationFitData.efield = efield;
    if (sRec.HasParameterError(revt::eAngleToLorentzVector)) {
//      tmpStationFitData.lorentzAngleError = fMCUncertaintyOnEFieldVector;
      tmpStationFitData.lorentzAngleError = sqrt(pow(1 * utl::degree, 2) +
                                                 pow(fMCUncertaintyOnEFieldVector, 2));
      info << "setting uncertainty to " << tmpStationFitData.lorentzAngleError/utl::deg << " degree\n";
    } else {
      WARNING(
          "RdLDFFitter: Parameter eAngleToLorentzVector not set. Please run the RdStationEFieldReconstructor Module before the LDF Fitter. This event will be skipped!");
      return eFailure;
    }

    vStationFitData.push_back(tmpStationFitData);
  }
  INFO2(eDebug, info.str());

  bool hasScintillator =
      showerrrec.HasParameter(revt::eNumberOfAERATopScintWithPulseFound)
          && (showerrrec.GetParameter(revt::eNumberOfAERATopScintWithPulseFound) >= 3) ?
          true : false;
  ParameterScintillatorLDFModel preScintillatorFit;
  std::vector<ScintillatorFitData> vScintillatorFitData;
  if(hasScintillator) {
    // fill scintillator fit data
    int nSignalScintiallators = 0;
    for (REvent::StationIterator sIt = rEvent.StationsBegin();
          sIt != rEvent.StationsEnd(); ++sIt) {
        Station& station = *sIt;
        const StationRecData& sRec = station.GetRecData();
        ScintillatorFitData tmpScintillatorData;
        if(sRec.HasParameter(revt::eAERAScintStationVEM) and sRec.GetParameter(revt::eAERAScintStationVEM) > 1.) {
          tmpScintillatorData.signal = sRec.GetParameter(revt::eAERAScintStationVEM);
          tmpScintillatorData.silent = false;
          nSignalScintiallators++;
        } else {
          tmpScintillatorData.silent = true;
        }
        const utl::Point sPos = rDetector.GetStation(station).GetPosition();
        tmpScintillatorData.scintillatorPosition = sPos;
        vScintillatorFitData.push_back(tmpScintillatorData);
    }
    if(nSignalScintiallators < 3) {
      hasScintillator = false;
      INFO2(eIntermediate, "number of scintillators is less than three, no scintillator LDF will be fitted");
    } else {
      // perform pre fit
      preScintillatorFit = PreScintillatorLDFFit(vScintillatorFitData, fShowerAxis, coordinateOrigin);
      info.str("");

      if(preScintillatorFit.status == 0 or preScintillatorFit.status >= 1000) {
        info << "Pre scintillator LDF fit was successful. N_e = " << preScintillatorFit.N_e
             << " +/- " << preScintillatorFit.N_eError;
      } else {
        info << "Pre scintillator LDF fit was not successful. Error code is " << preScintillatorFit.status;
      }
      INFO2(eIntermediate, info.str());
    }
  }

  ParameterLDFModel PreFitResult;
  if (fFitConfig.LDFModel) {
    INFO2(eDebug, "performing pre LDF fit");
    // perform a pre LDF fit with core = coordinate-origin (barycenter of the antennas)
    PreFitResult = PreLDFFit(fFitConfig.LDFModel, vStationFitData, fShowerAxis, coordinateOrigin);
  }

  // the FIT

  info.str("");
  info << "Creating LDF fit function with the following parameters: \n " << "LDF Model:\t"
       << fFitConfig.LDFModel << "\n" << "usePolarisation:\t" << fFitConfig.usePolarisation << "\n"
       << "chargeExcessStrength:\t" << fFitConfig.chargeExcessStrength << "\n"
       << "useChargeExcessCorrectionInLDFFit:\t" << fFitConfig.useChargeExcessCorrectionInLDFFit
       << "\n" << "useSDCoreToImproveRadioCore:\t" << fFitConfig.useSDCoreToImproveRadioCore
       << "\n" << "useScintillator:\t" << fFitConfig.useScintillator
       << std::endl;
  INFO2(eDebug, info.str());

  LDFLikelihoodFunction fitFunction = LDFLikelihoodFunction(fFitConfig, eventFitData,
                                                            vStationFitData,
                                                            vScintillatorFitData,
                                                            fMagneticFieldAxis);

  INFO2(eDebug, "setting fit parameters");

  ROOT::Minuit2::MnUserParameters upar;
  upar.Add("coreX", 0., 20);
  upar.Add("coreY", 0., 20);
  upar.Add("coreZ", 0., 10);

  upar.Add("a", fFitConfig.chargeExcessStrength, 0.01);
  upar.Fix("a");

  // add scintillator fit variable
  upar.Add("N_e", preScintillatorFit.N_e, 1.e5);  // number of charged particles
  upar.Add("s", 1.7, .1);    // shower age
  upar.Add("Rm", 35., 1.);  // moliere radius
  upar.Fix("s");
  upar.Fix("Rm");
  if(!hasScintillator) {
    upar.Fix("N_e");  // fix all scintialltor LDF parameters if not enough scnintillators are present
  }

  // please add parameters for different LDF models here
  if (fFitConfig.LDFModel == 1) {
    upar.Add("A", PreFitResult.A, PreFitResult.AError);
    upar.Add("R0", PreFitResult.R0, PreFitResult.R0Error);
  }

  ROOT::Minuit2::MnMigrad migrad(fitFunction, upar);
  INFO2(eDebug, "fixing parameter 2 (coreZ)");
  migrad.Fix("coreZ");

  // fix LDF parameters for first minimization step
  if (fFitConfig.LDFModel) {
    if (fFitConfig.LDFModel == 1) {
      upar.Fix("A");
      upar.Fix("R0");
      migrad();  // perform minimization
      upar.Release("A");
      upar.Release("R0");
    }
  }

  // perform minimization
  ROOT::Minuit2::FunctionMinimum minimum = migrad();  // minimize calles migrad algorithm. If migrad fails it calls simplex and then migrad again.
  info.str("");
  info << "minimum: " << minimum;
  INFO2(eIntermediate, info.str());

  ROOT::Minuit2::MnUserParameterState minUstate_tmp = minimum.UserState();
  utl::Point core_tmp(minUstate_tmp.Value("coreX"), minUstate_tmp.Value("coreY"), minUstate_tmp.Value("coreZ"), fLocalCS);
  if(event.HasSimShower()) {
    info.str("");
    double distance = (event.GetSimShower().GetPosition() - core_tmp).GetMag();
    info << "distance to MC core is " << distance / utl::m << " m";
    INFO2(eIntermediate, info.str());
  }
  if(event.GetRecShower().HasSRecShower()) {
    info.str("");
    double distance = (event.GetRecShower().GetSRecShower().GetCorePosition() - core_tmp).GetMag();
    info << "distance to SD core is " << distance / utl::m << " m";
    INFO2(eIntermediate, info.str());
  }

  // abort if fit was not successful
  if (!minimum.IsValid()) {
    WARNING("Core position fit not converged.");
    info.str("");
    info << "current rec stage is " << showerrrec.GetParameter(eRecStage) << std::endl;
    INFO2(eIntermediate, info.str());
    return eSuccess;
  }
  showerrrec.SetParameter(eRecStage, 2.0);

  double optimalChargeExcessStrength = fFitConfig.chargeExcessStrength;

  if(fFitConfig.fitChargeExcess) {
    INFO2(eIntermediate, "charge excess strength will be fitted.");
    migrad.Release("a");
    ROOT::Minuit2::FunctionMinimum minimum2 = migrad();
    if (!minimum2.IsValid()) {
        WARNING("Core position fit with free charge excess strength not converged. Saving result without variation of a");
    } else {
      INFO2(eIntermediate, "fit of charge excess strength was successful.");
      minimum = minimum2;
      info.str("");
      info << "minimum: " << minimum;
      INFO2(eIntermediate, info.str());

      showerrrec.SetParameter(revt::eChargeExcessStrength, minimum2.UserState().Value("a"));
      showerrrec.SetParameter(eRecStage, 2.1);
      optimalChargeExcessStrength = minimum2.UserState().Value("a");
    }
  }


  // get the minimal parameters
  ROOT::Minuit2::MnUserParameterState minUstate = minimum.UserState();
  ROOT::Minuit2::MnUserCovariance cov = minUstate.Covariance();

  std::cout << "covariance matrix \n" << cov << std::endl;

  // get probabilty that SD core is correct
  const double arr[] = {eventFitData.SDCore.GetX(fLocalCS),eventFitData.SDCore.GetY(fLocalCS), eventFitData.SDCore.GetZ(fLocalCS)};
  std::vector<double> params (arr, arr + sizeof(arr) / sizeof(arr[0]) );
  const double likelihoodValueAtSDCore = fitFunction(params);
  showerrrec.SetParameter(revt::eProbSDCore, TMath::Prob(likelihoodValueAtSDCore, 2));

  // save fit result
  // get reference coordinate system of detector (usually PampaAmarilla)
  const utl::CoordinateSystemPtr referenceCS = detector.GetReferenceCoordinateSystem();
  INFO2(eDebug, "save fit results");
  utl::Point core(minUstate.Value("coreX"), minUstate.Value("coreY"), minUstate.Value("coreZ"),
                  fLocalCS);

  showerrrec.SetParameter(revt::eCoreX, core.GetX(referenceCS));
  showerrrec.SetParameter(revt::eCoreY, core.GetY(referenceCS));
  showerrrec.SetParameter(revt::eCoreZ, core.GetZ(referenceCS));

  showerrrec.SetParameterCovariance(revt::eCoreX, revt::eCoreY, cov(0, 1));
  showerrrec.SetParameterCovariance(revt::eCoreX, revt::eCoreZ, 0);
  showerrrec.SetParameterCovariance(revt::eCoreY, revt::eCoreZ, 0);
  showerrrec.SetParameterCovariance(revt::eCoreX, revt::eCoreX, cov(0, 0));
  showerrrec.SetParameterCovariance(revt::eCoreY, revt::eCoreY, cov(1, 1));
  showerrrec.SetParameterCovariance(revt::eCoreZ, revt::eCoreZ, 0);

  info.str("");
  info << "core (local CS) = (" << minUstate.Value("coreX") << ","
                     << minUstate.Value("coreY") << "," << minUstate.Value("coreZ") << ")";
  INFO2(eIntermediate, info.str());
  info.str("");
  info << "core (reference CS) = (" << core.GetX(referenceCS) << ","
      << core.GetY(referenceCS) << "," << core.GetZ(referenceCS) << ")";
  INFO2(eIntermediate, info.str());
  info.str("");
  info << "optimal charge excess is " << minUstate.Value("a")*100. << "%";
  INFO2(eIntermediate, info.str());
  if(event.GetRecShower().HasSRecShower()) {
    info.str("");
    double distance = (event.GetRecShower().GetSRecShower().GetCorePosition() - core).GetMag();
    info << "distance to SD core is " << distance / utl::m << " m";
    INFO2(eIntermediate, info.str());
  }
  if(event.HasSimShower()) {
    info.str("");
    double distance = (event.GetSimShower().GetPosition() - core).GetMag();
    info << "distance to MC core is " << distance / utl::m << " m";
    INFO2(eIntermediate, info.str());
  }



  // save LDF
  if(fFitConfig.LDFModel){
    showerrrec.SetParameter(revt::eLDFParameter1, minUstate.Value("A"));
    showerrrec.SetParameter(revt::eLDFParameter2, minUstate.Value("R0"));
    showerrrec.SetParameterCovariance(revt::eLDFParameter1, revt::eLDFParameter1, cov(2, 2));
    showerrrec.SetParameterCovariance(revt::eLDFParameter2, revt::eLDFParameter2, cov(3, 3));
    showerrrec.SetParameterCovariance(revt::eLDFParameter1, revt::eLDFParameter2, cov(2, 3));
    // correlation with core position
    showerrrec.SetParameterCovariance(revt::eCoreX, revt::eLDFParameter1, cov(0, 2));
    showerrrec.SetParameterCovariance(revt::eCoreY, revt::eLDFParameter1, cov(1, 2));
    showerrrec.SetParameterCovariance(revt::eCoreX, revt::eLDFParameter2, cov(0, 3));
    showerrrec.SetParameterCovariance(revt::eCoreY, revt::eLDFParameter2, cov(1, 3));
    OUT(eIntermediate) << "LDF parameters at minimum: A = " << minUstate.Value("A") << "+/-"
                       << minUstate.Error("A") << "R0 = " << minUstate.Value("R0") << "+/-"
                       << minUstate.Error("R0") << std::endl;
  }

  if(fFitConfig.useScintillator and hasScintillator) {
    TF1 NKG = TF1("NKG","[0]*TMath::Gamma(4.5-[1])/(2.*TMath::Pi()*[2]*[2]*TMath::Gamma([1])*TMath::Gamma(4.5-(2.*[1])))*TMath::Power((x/[2]),([1]-2.))*TMath::Power((1+(x/[2])),[1]-4.5)");
    NKG.SetParameter(0, minUstate.Value("N_e"));
    NKG.SetParameter(1, minUstate.Value("s"));
    NKG.SetParameter(2, minUstate.Value("Rm"));
    showerrrec.SetScintLDF(NKG);
  }

  // scan
  utl::Point scanCenter;
  if(fScanCenter == eRadio) {
    scanCenter = showerrrec.GetCorePosition();
  }
  if(fScanCenter == eSD) {
    if (not event.GetRecShower().HasSRecShower()) {
      WARNING(
          "RdLDFFitter: No reconstructed SD shower found. SD core position can't be used as center of scan. Scan will not be performed.");
      fPerformScan = false;
    } else {
      evt::ShowerSRecData& showersrec = event.GetRecShower().GetSRecShower();  // get reconstructed SD shower for SD core position
      scanCenter = showersrec.GetCorePosition();
    }
  }
  if(fScanCenter == eFD) {
    if (not event.GetRecShower().HasFRecShower()) {
      WARNING(
          "RdLDFFitter: No reconstructed FD shower found. FD core position can't be used as center of scan. Scan will not be performed.");
      fPerformScan = false;
    } else {
      evt::ShowerFRecData& showerfrec = event.GetRecShower().GetFRecShower();  // get reconstructed SD shower for SD core position
      scanCenter = showerfrec.GetCorePosition();
    }
  }
  if(fScanCenter == eMC) {
    if (not event.HasSimShower()) {
      WARNING(
          "RdLDFFitter: No reconstructed MC shower found. MC core position can't be used as center of scan. Scan will not be performed.");
      fPerformScan = false;
    } else {
      scanCenter  = event.GetSimShower().GetPosition();
    }
  }
  utl::Point coreScan;
  if(fPerformScan) {
    upar.Fix("a");
    fNStepsX = (fNStepsX % 2) == 0 ? fNStepsX + 1 : fNStepsX; // make n-steps odd
    fNStepsY = (fNStepsY % 2) == 0 ? fNStepsY + 1 : fNStepsY;
    std::vector< std::vector<double> > scanResult(fNStepsX, std::vector<double>(fNStepsY));
    coreScan = Scan(scanResult, fitFunction, scanCenter, fNStepsX, fNStepsY,
                                fStepWidthX, fStepWidthY, optimalChargeExcessStrength);
    showerrrec.SetParameter(revt::eCoreScanX, coreScan.GetX(referenceCS));
    showerrrec.SetParameter(revt::eCoreScanY, coreScan.GetY(referenceCS));
    showerrrec.SetParameter(revt::eCoreScanZ, coreScan.GetZ(referenceCS));

    std::stringstream filename;
    filename << fOutputFolder << "/scan_" << std::setfill('0') << std::setw(6) << runNumber << "_"
             << std::setfill('0') << std::setw(7) << eventId << ".png";
    SaveScan(scanResult, filename.str(), scanCenter.GetX(referenceCS),
             scanCenter.GetY(referenceCS), fNStepsX, fNStepsY, fStepWidthX, fStepWidthY);
    filename.str("");
    filename << fOutputFolder << "/scan_" << std::setfill('0') << std::setw(6) << runNumber << "_"
             << std::setfill('0') << std::setw(7) << eventId << ".png";
    PlotScan(scanResult, filename.str(),
             scanCenter, fLocalCS,
             scanCenter, core,
             fNStepsX, fNStepsY, fStepWidthX, fStepWidthY,
             vStationFitData, eventFitData, optimalChargeExcessStrength);
    filename.str("");
    filename << std::setfill('0') << std::setw(6) << runNumber << "_" << std::setfill('0')
             << std::setw(7) << eventId;
    PlotGoodnessOfFit(fOutputFolder, filename.str(), core);
  }

  // save polarization angles for different cores
  std::vector<std::pair<double, double> > anglesCore = GetAnglesToEFieldExpectation(core, optimalChargeExcessStrength);
  std::vector<std::pair<double, double> > anglesCoreScan = GetAnglesToEFieldExpectation(coreScan, optimalChargeExcessStrength);
  std::vector<std::pair<double, double> > anglesReferenceCore = GetAnglesToEFieldExpectation(scanCenter, optimalChargeExcessStrength);
  unsigned int i = 0;
  info.str("");
  info << "betas at fit optimum are: ";
  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    Station& station = *sIt;
    StationRecData& sRec = station.GetRecData();
    sRec.SetParameter(revt::eBetaFit, anglesCore[i].first);
    info << anglesCore[i].first / utl::deg << "deg  ";
    // probability to measure the deviation beta or something greater
    sRec.SetParameter(
        revt::eProbBetaFit,1 - utl::Probability::GetFisherCDF(anglesCore[i].first,
                                             1 / (anglesCore[i].second * anglesCore[i].second)));

    if(fPerformScan) {
      sRec.SetParameter(revt::eBetaScan, anglesCoreScan[i].first);
      // probability to measure the deviation beta or something greater
      sRec.SetParameter(
          revt::eProbBetaScan, 1 - utl::Probability::GetFisherCDF(anglesCoreScan[i].first,
                  1 / (anglesCoreScan[i].second * anglesCoreScan[i].second)));
    }

    sRec.SetParameter(revt::eBetaReferenceCore, anglesReferenceCore[i].first);
    // probability to measure the deviation beta or something greater
    sRec.SetParameter(
        revt::eProbBetaReferenceCore, 1 - utl::Probability::GetFisherCDF(anglesReferenceCore[i].first,
                1 / (anglesReferenceCore[i].second * anglesReferenceCore[i].second)));
    ++i;
  }
  INFO2(eDebug, info.str());



  if(fFitConfig.LDFModel) {
    // compute cosmic ray energy
    TF1 fLDF = TF1("fLDF", "[0]*exp(-x/[1])");
    fLDF.SetParameter(0, minUstate.Value("A"));
    fLDF.SetParameter(1, minUstate.Value("R0"));
    double energy_estimator = fLDF.Eval(110);
    double energy = GetEnergy(energy_estimator);
    showerrrec.SetParameter(revt::eEnergy, energy);

    std::stringstream sstr;
    sstr << "energy = " << energy << " eV";
    INFO2(eIntermediate, sstr.str());
  }

  if (fSaveContour) {
    // get uncertainty and plot confidence contour
    ROOT::Minuit2::MnContours contours(fitFunction, minimum);
    fitFunction.SetErrorDef(2.3);   // 68% confidence
    std::vector<std::pair<double, double> > cont = contours(0, 1, 40);

    fitFunction.SetErrorDef(5.99);  // 95% confidence
    std::vector<std::pair<double, double> > cont2 = contours(0, 1, 40);

    // save contour
    VModule::ResultFlag result = SaveContours(runNumber, eventId, cont, cont2, cov,
                                              minUstate.Value("coreX"), minUstate.Value("coreY"));
    if (result != eSuccess) {
      return result;
    }

//    ROOT::Minuit2::MnPlot plot;
//    plot(minimum.UserState().Value("coreX"), minimum.UserState().Value("coreY"), cont);
  }

  return eSuccess;
}


VModule::ResultFlag RdLDFFitter::Finish()
{
  // Put any termination or cleanup code here.
  // This method is called once at the end of the run.

  INFO("RdLDFFitter::Finish()");

  return eSuccess;
}



utl::Point RdLDFFitter::Scan(std::vector< std::vector<double> >& scanResult, const LDFLikelihoodFunction L,
                                      const utl::Point core,
                                      const unsigned int nStepsX, const unsigned int nStepsY,
                                      const double widthX, const double widthY, const double a)
{
  const double coreX = core.GetX(fLocalCS);
  const double coreY = core.GetY(fLocalCS);
  const double coreZ = core.GetZ(fLocalCS);

  std::stringstream msg;
  msg.str("");
  msg << "Scan of Likelihood function around " << coreX << ", " << coreY << " will be performed\n"
      << nStepsX << " steps in x direction with width " << widthX << "\n"
      << nStepsY << " steps in y direction with width " << widthY;
  INFO2(eIntermediate, msg.str().c_str());

  double Lmin = 9999999999;
  double xmin(0), ymin(0);

  msg.str("");
  msg << "loop from " << - (int) (nStepsX / 2) << " to " << nStepsX / 2;
  INFO2(eDebug, msg.str().c_str());
  msg.str("");
  for (unsigned int iX = 0; iX < nStepsX; ++iX) {
    const double x = coreX + static_cast<int>(iX - nStepsX/2) * widthX;
    for (unsigned int iY = 0; iY < nStepsY; ++iY) {
      const double y = coreY + static_cast<int>(iY - nStepsY/2) * widthY;
      const double arr[] = {x, y, coreZ, a};
      std::vector<double> params (arr, arr + sizeof(arr) / sizeof(arr[0]) );
      double temp = L(params);
      msg << temp << " ";
      scanResult[iX][iY] = temp;
      if(temp < Lmin) {
        Lmin = temp;
        xmin = x;
        ymin = y;
      }
    }
  }
  msg << "\bLmin is " << Lmin;
  INFO2(eDebug, msg.str());

  msg.str("");
  msg << "Scan result: delta -2*logLikelihood at given core position is "
      << scanResult[nStepsX/2+1][nStepsY/2+1]<< "\n"
      << "this corresponds to a probability of " << TMath::Prob(scanResult[nStepsX/2+1][nStepsY/2+1], 2) << "\n"
      << "this corresponds to sigma = " << TMath::NormQuantile(TMath::Prob(scanResult[nStepsX/2+1][nStepsY/2+1], 2));
  INFO2(eIntermediate, msg.str());

  for(std::vector< std::vector<double> >::iterator it = scanResult.begin(); it != scanResult.end(); ++it) {
    for(std::vector<double>::iterator it2 = it->begin(); it2 != it->end(); ++it2) {
      *it2 = *it2 - Lmin;
    }
  }
  return utl::Point(xmin, ymin, coreZ, fLocalCS);
}


std::string GetFormattedNumber(const double number) {
  std::stringstream msg;
  msg.str("");
  if(number >= 0.01) {
    msg << std::setprecision(2) << std::fixed << number;
  } else {
    msg << std::setprecision(1) << std::scientific << number;
  }
  return msg.str();
}

fwk::VModule::ResultFlag RdLDFFitter::PlotScan(const std::vector<std::vector<double> >& scanResult,
                                                 const std::string filename,
                                                 const utl::Point scanCenter, const utl::CoordinateSystemPtr cs,
                                                 const utl::Point referenceCore, const utl::Point coreRd,
                                                 const unsigned int nStepsX, const unsigned int nStepsY,
                                                 const double widthX, const double widthY,
                                                 const std::vector<StationFitData>& vStationFitData,
                                                 const EventFitData& eventFitData, const double a)
{
  bool plotDebug = false;
  const det::Detector& detector = det::Detector::GetInstance();
  const rdet::RDetector& rDetector = detector.GetRDetector();
  const double coreX = scanCenter.GetX(cs);
  const double coreY = scanCenter.GetY(cs);
  // define sigma contours in 2D

  TH2D h = TH2D("h", "", nStepsX,
                coreX - static_cast<int>(nStepsX / 2) * widthX - 0.5 * widthX,
                coreX + static_cast<int>(nStepsX / 2) * widthX + 0.5 * widthX, nStepsY,
                coreY - static_cast<int>(nStepsY / 2) * widthY - 0.5 * widthY,
                coreY + static_cast<int>(nStepsY / 2) * widthY + 0.5 * widthY);
  double minimum = 9999999.;
  std::pair<double, double> minimum_bin (0,0);
  std::stringstream msg;
  msg.str("");
  for (unsigned int iY = 0; iY < nStepsY; ++iY) {
    const double y = coreY + static_cast<int>(iY - nStepsY/2) * widthY;
    for (unsigned int iX = 0; iX < nStepsX; ++iX) {
      const double x = coreX + static_cast<int>(iX - nStepsX/2) * widthX;
      msg << scanResult[iX][iY] << " ";
      h.Fill(x, y, scanResult[iX][iY]);
      if(scanResult[iX][iY] < minimum) {
        minimum = scanResult[iX][iY];
        minimum_bin.first = x;
        minimum_bin.second = y;
      }
    }
  }
  INFO2(eDebug, msg.str());

  TCanvas c1("c1", "", 1000,822);
  c1.cd();
  c1.GetPad(0)->SetTopMargin(0.02);
  c1.GetPad(0)->SetLeftMargin(0.10);
  c1.GetPad(0)->SetBottomMargin(0.08);
  c1.GetPad(0)->SetRightMargin(0.16);
  h.GetXaxis()->SetTitle("x [m]");
  h.GetYaxis()->SetTitle("y [m]");
  h.GetYaxis()->SetTitleOffset(1.2);
  h.GetZaxis()->SetTitle("-2 ln(L)");
  h.GetZaxis()->SetTitleOffset(1.3);
  if(fPlotSigmaContours) {
    const unsigned int nContours = 9;
    const double contours[nContours] = { -0.00001, 2.29574893, 6.18007431, 11.82915808, 19.33390861, 28.74370243,
          40.08724343, 53.3822385, 68.50378784};
    h.SetContour(nContours, contours);
    h.GetZaxis()->SetRangeUser(-0.00001, 1.1*68.50378784);
  }
  h.SetStats(false);
  h.Draw("COLZ");

  TLegend leg = TLegend(0.6,0.8,0.84,0.98);


  // plot radio stations
  if(plotDebug) {
    const std::vector<int>& stationList = rDetector.GetFullStationList();
    TGraph gPosAll = TGraph(stationList.size());
    for(unsigned int i = 0; i< stationList.size(); ++i) {
      const utl::Point position = rDetector.GetStation(stationList[i]).GetPosition();
      gPosAll.SetPoint(i, position.GetX(cs), position.GetY(cs));
    }
    gPosAll.SetMarkerColor(1);
    gPosAll.SetMarkerStyle(8);
    gPosAll.SetMarkerSize(2);
    gPosAll.Draw("Psame");
    leg.AddEntry(&gPosAll,"RDS", "P");
  }

  TGraph gPos = TGraph(vStationFitData.size());
  for(unsigned int i = 0; i< vStationFitData.size(); ++i) {
    gPos.SetPoint(i, vStationFitData[i].antennaPosition.GetX(cs), vStationFitData[i].antennaPosition.GetY(cs));
  }
  if(plotDebug) {
    gPos.SetMarkerColor(1);
    gPos.SetMarkerStyle(22);
    gPos.SetMarkerSize(2);
  } else {
    gPos.SetMarkerColor(1);
    gPos.SetMarkerStyle(8);
    gPos.SetMarkerSize(2);
  }
  gPos.Draw("Psame");
  leg.AddEntry(&gPos,"Rd station", "P");


  // compute distance of SD core position to likelihood minima in terms of sigma
  const double sigma = TMath::Prob(scanResult[nStepsX/2][nStepsY/2], 2);

  // plot SD core
  TGraph gSDCore(1);
  gSDCore.SetPoint(0, referenceCore.GetX(cs), referenceCore.GetY(cs));
  gSDCore.SetMarkerStyle(29);
  gSDCore.SetMarkerColor(6);
  gSDCore.SetMarkerSize(4);
  gSDCore.Draw("Psame");
  msg.str("");
  if(fConstEvent->GetRecShower().HasSRecShower()) {
    msg << "SD";
  }
  if(fConstEvent->HasSimShower()) {
    msg << "MC";
  }
  msg << " core";
  if(plotDebug) {
    msg << "(p=" << GetFormattedNumber(sigma) << ")";
  }
  leg.AddEntry(&gSDCore, msg.str().c_str(), "P");

  // draw incoming azimuth direction
  TLine line = TLine();
  line.SetLineWidth(2);
  const double azimuth = fConstEvent->GetRecShower().GetRRecShower().GetAxis().GetPhi(fLocalCS);
  const double length = nStepsX * widthX * 0.25;
  line.DrawLine(referenceCore.GetX(cs), referenceCore.GetY(cs), referenceCore.GetX(cs) + cos(azimuth) * length,
                referenceCore.GetY(cs) + sin(azimuth) * length);

  // plot improved SD core
//  utl::Point coreSdImproved;
//  if(fConstEvent->GetRecShower().HasSRecShower()) {
//    coreSdImproved = fConstEvent->GetRecShower().GetSRecShower().GetCorrectedCorePosition();
//    TGraph gSDCoreImproved(1);  // graph needs to be defined to make it variable available in other if clauses
//    if(fPlotStyle > 0) {
//      gSDCoreImproved.SetPoint(0, coreSdImproved.GetX(cs), coreSdImproved.GetY(cs));
//      gSDCoreImproved.SetMarkerStyle(29);
//      gSDCoreImproved.SetMarkerColor(2);
//      gSDCoreImproved.SetMarkerSize(4);
//      gSDCoreImproved.Draw("Psame");
//      msg.str("");
//      msg << "SD core improved (" << std::setprecision(0) << std::fixed
//          << coreSdImproved.GetX(cs) / utl::m << "m," << coreSdImproved.GetY(cs) / utl::m << "m)";
//      leg.AddEntry(&gSDCoreImproved, msg.str().c_str(), "P");
//    }
//  }

  // plot Rd core
  TGraph gRdCore(1);
  gRdCore.SetPoint(0, coreRd.GetX(cs), coreRd.GetY(cs));
  gRdCore.SetMarkerStyle(23);
  gRdCore.SetMarkerColor(6);
  gRdCore.SetMarkerSize(3);
  gRdCore.Draw("Psame");

  msg.str("");
  msg << "Rd core";
  if(plotDebug) {
    msg << "(" << std::setprecision(0) << std::fixed << coreRd.GetX(cs) / utl::m << "m," << coreRd.GetY(cs) / utl::m << "m)";
  }
  leg.AddEntry(&gRdCore, msg.str().c_str(), "P");

  if(fConstEvent->GetRecShower().HasSRecShower()) {
    double sizeFactor = 1.51;  // to get 68 % coverage for one sigma intervals in 2D
    double sx2 = eventFitData.SDCoreXError * eventFitData.SDCoreXError;
    double sy2 = eventFitData.SDCoreYError * eventFitData.SDCoreYError;
    double sxy = eventFitData.SDCoreXError * eventFitData.SDCoreYError * eventFitData.SDCoreXYCorrelation;
    double errorEllipseAngle = 0.5 * atan2(2 * sxy, sx2 - sy2);
    double sa = sin(errorEllipseAngle);
    double ca = cos(errorEllipseAngle);
    double ru = sizeFactor * sqrt(sx2 * ca * ca + 2 * sxy * sa * ca + sy2 * sa * sa);
    double rv = sizeFactor * sqrt(sx2 * sa * sa - 2 * sxy * sa * ca + sy2 * ca * ca);
    TEllipse cc = TEllipse(referenceCore.GetX(cs), referenceCore.GetY(cs), ru, rv, 0., 360.,
                  errorEllipseAngle / utl::deg);
    msg.str("");
    msg << "plotting SD error ellipse with " << referenceCore.GetX(cs) << ", "
        << referenceCore.GetY(cs) << ", " << ru << ", " <<  rv << ", " << errorEllipseAngle / utl::deg;
    INFO2(eDebug, msg.str());
    cc.SetFillColor(16);
    cc.SetLineColor(15);
    cc.SetFillStyle(3017);
    cc.Draw("F");
  }

  msg.str("");
  msg << "scan min (" << std::setprecision(0) << std::fixed << minimum_bin.first << "m,"
      << minimum_bin.second << "m)";
  if(plotDebug) {
    leg.AddEntry(&gRdCore, msg.str().c_str(), "");
  }
  leg.Draw();

  // rd core from fit
  std::vector<std::pair<double, double> > angles = GetAnglesToEFieldExpectation(coreRd, a);
  double meanAngle(0.);
  //double meanAngleSigma(0.);
  unsigned int n = angles.size();
  double p_value_combined(0.);
  double p_Likelihood = 0.;
  for(unsigned int i = 0; i < n; ++i) {
    meanAngle += angles[i].first;
    //meanAngleSigma += angles[i].first/angles[i].second;
    const double p_value = 1 - utl::Probability::GetFisherCDF(angles[i].first,
                                         1 / (angles[i].second * angles[i].second));
    p_value_combined += log(p_value);
    p_Likelihood += -2 * log(utl::Probability::GetFisherPDF(angles[i].first,
                                           1 / (angles[i].second * angles[i].second)));
  }
  p_value_combined *= -2.;
  meanAngle /= n;
  //meanAngleSigma /= n;

  TLatex text;
  text.SetNDC(true);
  msg.str("");
  msg << "#bar{#beta}_{fit} = " << std::setprecision(1) << std::fixed << TMath::RadToDeg() * meanAngle
      << "#circ, p=" << GetFormattedNumber(TMath::Prob(p_value_combined, 2 * n));
  if(fPlotStyle > 0) {
    msg << ", L=" << p_Likelihood;
  }
  text.SetTextSize(0.030);
  double yy = 0.77;
  double xx = 0.6;
  if(plotDebug) {
    yy = 0.67;
    xx = 0.68;
  }
  const double dy = 0.035;
  //  text.DrawLatex(xx, yy, msg.str().c_str());
  //  yy -= dy;
  fEvent->GetRecShower().GetRRecShower().SetParameter(revt::eProbRdCorePol, TMath::Prob(p_value_combined, 2 * n));


  // scan minima
  angles = GetAnglesToEFieldExpectation(
      utl::Point(minimum_bin.first, minimum_bin.second, coreRd.GetZ(cs), cs), a);
  meanAngle = 0.;
  //meanAngleSigma = 0.;
  n = angles.size();
  p_value_combined = 0.;
  p_Likelihood = 0.;
  for(unsigned int i = 0; i < n; ++i) {
    meanAngle += angles[i].first;
    //meanAngleSigma += angles[i].first/angles[i].second;
    const double p_value = 1 - utl::Probability::GetFisherCDF(angles[i].first,
                                         1 / (angles[i].second * angles[i].second));
    p_value_combined += log(p_value);
    p_Likelihood += -2 * log( utl::Probability::GetFisherPDF(angles[i].first,
                                           1 / (angles[i].second * angles[i].second)));
  }
  p_value_combined *= -2.;
  meanAngle /= n;
  //meanAngleSigma /= n;
  msg.str("");
  msg << "#bar{#beta}_{scan} = " << std::setprecision(1) << std::fixed << TMath::RadToDeg() * meanAngle
      << "#circ, p=" << GetFormattedNumber(TMath::Prob(p_value_combined, 2 * n));
  if(fPlotStyle > 0) {
    msg << ", L=" << p_Likelihood;
  }
  //  text.DrawLatex(xx, yy, msg.str().c_str());
  //  yy -= dy;
  fEvent->GetRecShower().GetRRecShower().SetParameter(revt::eProbRdCoreScanPol, TMath::Prob(p_value_combined, 2 * n));

  // SD core
  angles = GetAnglesToEFieldExpectation(referenceCore, a);
  meanAngle = 0.;
  //meanAngleSigma = 0.;
  n = angles.size();
  p_value_combined = 0.;
  p_Likelihood = 0.;
  for(unsigned int i = 0; i < n; ++i) {
    meanAngle += angles[i].first;
    //meanAngleSigma += angles[i].first/angles[i].second;
    const double p_value = 1
        - utl::Probability::GetFisherCDF(angles[i].first,
                                         1 / (angles[i].second * angles[i].second));
    p_value_combined += log(p_value);
    p_Likelihood += -2
        * log(
            utl::Probability::GetFisherPDF(angles[i].first,
                                           1 / (angles[i].second * angles[i].second)));
  }
  p_value_combined *= -2.;
  meanAngle /= n;
  //meanAngleSigma /= n;
  msg.str("");
  msg << "#bar{#beta}_{SD} = " << std::setprecision(1) << std::fixed << TMath::RadToDeg() * meanAngle
      << "#circ, p=" << GetFormattedNumber(TMath::Prob(p_value_combined, 2 * n));
  if(fPlotStyle > 0) {
    msg << ", L=" << p_Likelihood;
  }
//  text.DrawLatex(xx, yy, msg.str().c_str());
//  yy -= dy;
  fEvent->GetRecShower().GetRRecShower().SetParameter(revt::eProbSdCorePol, TMath::Prob(p_value_combined, 2 * n));

//  // improved SD core
//  if(fConstEvent->GetRecShower().HasSRecShower() and fPlotStyle > 0) {
//    angles = GetAnglesToEFieldExpectation(coreSdImproved);
//    meanAngle = 0.;
//    meanAngleSigma = 0.;
//    n = angles.size();
//    p_value_combined = 0.;
//    p_Likelihood = 0.;
//    for(unsigned int i = 0; i < n; ++i) {
//      meanAngle += angles[i].first;
//      meanAngleSigma += angles[i].first/angles[i].second;
//      const double p_value = 1
//          - utl::Probability::GetFisherCDF(angles[i].first,
//                                           1 / (angles[i].second * angles[i].second));
//      p_value_combined += log(p_value);
//      p_Likelihood += -2
//          * log(
//              utl::Probability::GetFisherPDF(angles[i].first,
//                                             1 / (angles[i].second * angles[i].second)));
//    }
//    p_value_combined *= -2.;
//    meanAngle /= n;
//    meanAngleSigma /= n;
//    msg.str("");
//    msg << "#bar{#beta}_{SD+} = " << std::setprecision(1) << std::fixed << TMath::RadToDeg() * meanAngle
//        << "#circ, p=" << GetFormattedNumber(TMath::Prob(p_value_combined, 2 * n));
//    if(fPlotStyle > 0) {
//      msg << ", L=" << p_Likelihood;
//    }
//    text.DrawLatex(xx, yy, msg.str().c_str());
//    yy -= dy;
//  }
//  fEvent->GetRecShower().GetRRecShower().SetParameter(revt::eProbSdCoreImprovedPol, TMath::Prob(p_value_combined, 2 * n));

  // angel to Lorentz
  angles = GetAnglesToLorentzVector();
  meanAngle = 0.;
  double meanSigma = 0.;
  double minSigma = 999.;
  n = angles.size();
  for(unsigned int i = 0; i < n; ++i) {
    meanAngle += angles[i].first;
    meanSigma += angles[i].second;
    if(angles[i].second < minSigma) {
      minSigma = angles[i].second;
    }
  }
  meanAngle /= n;
  meanSigma /= n;
  msg.str("");
  msg << "#delta = " << std::setprecision(1) << std::fixed << TMath::RadToDeg() * meanAngle
      << "#circ"
      << " #sigma_{min}=" << minSigma * TMath::RadToDeg() << "#circ #bar{#sigma}=" << meanSigma  * TMath::RadToDeg() << "#circ";
//  text.DrawLatex(xx, yy, msg.str().c_str());
//  yy -= dy;


  if(fConstEvent->GetRecShower().HasSRecShower()) {
    msg.str("");
    const double energy = fConstEvent->GetRecShower().GetSRecShower().GetEnergy();
    msg << "E = " << std::setprecision(2) << std::fixed << energy / utl::EeV << "EeV";
    text.DrawLatex(xx, yy, msg.str().c_str());
    yy -= dy;
    // plot zenith
    msg.str("");
    const double zenith = fConstEvent->GetRecShower().GetSRecShower().GetAxis().GetTheta(fLocalCS);
    msg << "#Theta = " << std::setprecision(1) << std::fixed << zenith / utl::deg << "#circ";
    text.DrawLatex(xx, yy, msg.str().c_str());
    yy -= dy;
  }
  if(fConstEvent->HasSimShower()) {
    msg.str("");
    const double energy = fConstEvent->GetSimShower().GetEnergy();
    msg << "E = " << std::setprecision(2) << std::fixed << energy / utl::EeV << "EeV";
    text.DrawLatex(xx, yy, msg.str().c_str());
    yy -= dy;
    // plot zenith
    msg.str("");
    const utl::CoordinateSystemPtr localCS = fConstEvent->GetSimShower().GetLocalCoordinateSystem();
    const double zenith = (-fConstEvent->GetSimShower().GetDirection()).GetTheta(localCS); 
    msg << "#Theta = " << std::setprecision(1) << std::fixed << zenith / utl::deg << "#circ";
    text.DrawLatex(xx, yy, msg.str().c_str());
    yy -= dy;
  }

  msg.str("");
  msg << "a = " << std::setprecision(2) << a;
  text.DrawLatex(xx, yy, msg.str().c_str());
  yy -= dy;

  double distance = (referenceCore - coreRd).GetMag();
  msg.str("");
  msg << "D_{to SD} = " << std::setprecision(1) << distance / utl::m << " m";
  text.DrawLatex(xx, yy, msg.str().c_str());
  // yy -= dy;  // not useful here



  c1.SaveAs(filename.c_str());
  c1.SaveAs((filename+".C").c_str());

  return eSuccess;
}

std::vector<std::pair<double, double> > RdLDFFitter::GetAnglesToEFieldExpectation(const utl::Point core, const double a)
{
  const det::Detector& detector = det::Detector::GetInstance();
  const rdet::RDetector& rDetector = detector.GetRDetector();
  std::vector<std::pair<double, double> > angles;
  for (REvent::ConstSignalStationIterator sIt = fREvent->SignalStationsBegin();
      sIt != fREvent->SignalStationsEnd(); ++sIt) {
    const StationRecData& sRec = sIt->GetRecData();
    utl::Vector currentEField = utl::Normalized(
        utl::Vector(sRec.GetParameter(revt::eEFieldVectorEW),
                    sRec.GetParameter(revt::eEFieldVectorNS),
                    sRec.GetParameter(revt::eEFieldVectorV), fLocalCS));
    const utl::Point stationPosition = rDetector.GetStation(*sIt).GetPosition();
    double angleToExpectedEField= utl::RadioGeometryUtilities::GetAngleToEFieldExpectation2D(currentEField,
        core, fShowerAxis, stationPosition, fMagneticFieldAxis, a, fLocalCS);
    if (angleToExpectedEField > TMath::Pi() / 2.) {
      angleToExpectedEField = TMath::Pi() - angleToExpectedEField;
    }
    angles.push_back(std::pair<double, double>(angleToExpectedEField,
                                  sRec.GetParameterError(revt::eAngleToLorentzVector)));
  }
  return angles;
}

std::vector<std::pair<double, double> > RdLDFFitter::GetAnglesToLorentzVector()
{
  const det::Detector& detector = det::Detector::GetInstance();
  const rdet::RDetector& rDetector = detector.GetRDetector();
  std::vector<std::pair<double, double> > angles;
  for (REvent::ConstSignalStationIterator sIt = fREvent->SignalStationsBegin();
      sIt != fREvent->SignalStationsEnd(); ++sIt) {
    const StationRecData& sRec = sIt->GetRecData();
    utl::Vector currentEField = utl::Normalized(
        utl::Vector(sRec.GetParameter(revt::eEFieldVectorEW),
                    sRec.GetParameter(revt::eEFieldVectorNS),
                    sRec.GetParameter(revt::eEFieldVectorV), fLocalCS));
    const utl::Point stationPosition = rDetector.GetStation(*sIt).GetPosition();
    utl::Vector expectedEField = utl::RadioGeometryUtilities::GetLorentzVector(fShowerAxis, fMagneticFieldAxis);
    double angleToExpectedEField = utl::Angle(expectedEField, currentEField);
    if (angleToExpectedEField > TMath::Pi() / 2.) {
      angleToExpectedEField = TMath::Pi() - angleToExpectedEField;
    }
    if(fabs(angleToExpectedEField -  sRec.GetParameter(revt::eAngleToLorentzVector)) < 0.0001) {
      std::cout << "angle to Lorentz angle unequal " << std::endl;
    }
    angles.push_back(
        std::pair<double, double>(angleToExpectedEField,
                                  sRec.GetParameterError(revt::eAngleToLorentzVector)));
  }
  return angles;
}

fwk::VModule::ResultFlag RdLDFFitter::PlotGoodnessOfFit(const std::string folder,
                                                        const std::string eventIdentifier,
                                                        const utl::Point core)
{
  const det::Detector& detector = det::Detector::GetInstance();
  const rdet::RDetector& rDetector = detector.GetRDetector();
  TH1D hAngleToExpectedEField = TH1D("hAngleToExpectedEField", "", 18,0,90);
  hAngleToExpectedEField.GetXaxis()->SetTitle("#angle(#vec{E}_{exp}, #vec{E}_{measured})");
  hAngleToExpectedEField.GetYaxis()->SetTitle("entries");

  TH1D hAngleToExpectedEFieldSigma = TH1D("hAngleToExpectedEFieldSigma", "", 20,0,10);
  hAngleToExpectedEFieldSigma.GetXaxis()->SetTitle("#angle(#vec{E}_{exp}, #vec{E}_{measured})/#sigma_{#vec{E}_{measured}}");
  hAngleToExpectedEFieldSigma.GetYaxis()->SetTitle("entries");

  for (REvent::ConstSignalStationIterator sIt = fREvent->SignalStationsBegin();
        sIt != fREvent->SignalStationsEnd(); ++sIt) {
//      Station& station = *sIt;
      const StationRecData& sRec = sIt->GetRecData();
      utl::Vector currentEField = utl::Normalized(
        utl::Vector(sRec.GetParameter(revt::eEFieldVectorEW),
                    sRec.GetParameter(revt::eEFieldVectorNS),
                    sRec.GetParameter(revt::eEFieldVectorV), fLocalCS));
      const utl::Point stationPosition = rDetector.GetStation(*sIt).GetPosition();
      utl::Vector expectedEField = utl::RadioGeometryUtilities::GetExpectedEFieldVector(
          core, fShowerAxis, stationPosition, fMagneticFieldAxis, fFitConfig.chargeExcessStrength);
      double angleToExpectedEField = utl::Angle(expectedEField, currentEField);
      if (angleToExpectedEField > TMath::Pi() / 2) {
        angleToExpectedEField = TMath::Pi() - angleToExpectedEField;
      }
      hAngleToExpectedEField.Fill(angleToExpectedEField*TMath::RadToDeg());
      hAngleToExpectedEFieldSigma.Fill(angleToExpectedEField/sRec.GetParameterError(revt::eAngleToLorentzVector));
  }
  TCanvas c1("c1", "", 800,800);
  c1.cd();
  hAngleToExpectedEField.SetLineWidth(2);
  hAngleToExpectedEField.Draw();
  c1.SaveAs((folder+"/hAngleToEFieldExectation_"+eventIdentifier+".png").c_str());

  TCanvas c2("c2", "", 800,800);
  c2.cd();
  hAngleToExpectedEFieldSigma.SetLineWidth(2);
  hAngleToExpectedEFieldSigma.Draw();
  c2.SaveAs((folder+"/hAngleToEFieldExectation_"+eventIdentifier+"_sigma.png").c_str());
  return eSuccess;
}

VModule::ResultFlag RdLDFFitter::SaveScan(const std::vector< std::vector<double> >& scanResult,
                                            const std::string filename,
                                            const double coreX, const double coreY,
                                            const unsigned int nStepsX, const unsigned int nStepsY,
                                            const double widthX, const double widthY) {
  std::ofstream myfile (filename.c_str());
  if (myfile.is_open()) {
    for (unsigned int iX = 0; iX < nStepsX; ++iX) {
        const double x = coreX + static_cast<int>(iX - nStepsX/2) * widthX;
        myfile << "\t" << x;
    }
    myfile << "\n";

    for (unsigned int iY = 0; iY < nStepsY; ++iY) {
      const double y = coreY + static_cast<int>(iY - nStepsY/2) * widthY;
      myfile << y;
      for (unsigned int iX = 0; iX < nStepsX; ++iX) {
          myfile << "\t" << scanResult[iX][iY];
      }
      myfile << "\n";
    }
    myfile.close();
    return eSuccess;
  }
  else WARNING("Unable to open file");
  return eFailure;
}

double RdLDFFitter::GetEnergy(const double energy_estimator)
{
  return 1.35e15 * pow(energy_estimator / utl::micro / utl::volt * utl::meter, 0.98) * utl::eV;
}

VModule::ResultFlag RdLDFFitter::SaveContours(
    unsigned int runNumber, const int eventId,
    const std::vector<std::pair<double, double> > contour,
    const std::vector<std::pair<double, double> > contour2, ROOT::Minuit2::MnUserCovariance cov,
    const double coreX, const double coreY)
{
  if (!boost::filesystem::exists("output")) {
    WARNING(
        "folder \"output\" does not exist. Can't save confidence contours of radio core position fit. Please create the output folder first.");
    return eFailure;
  }
  std::stringstream sstr("");

  sstr << fOutputFolder << "/contour_" << std::setfill('0') << std::setw(6) << runNumber << "_"
       << std::setfill('0') << std::setw(7) << eventId << ".dat";
  std::fstream fout;
  fout.open(sstr.str().c_str(), std::ios::out);
  // save fit result
  fout << "# fir result: coreX  coreY \n" << coreX << "\t" << coreY << "\n";
  // save covariance matrix
  fout << "# covariance matrix \n" << cov(0, 0) << "\t" << cov(0, 1) << "\n" << cov(1, 0) << "\t"
       << cov(1, 1) << "\n";
  fout << "# 68% confidence region: x [m]\t y[m]  in local coordinates \n";
  for (std::vector<std::pair<double, double> >::const_iterator it = contour.begin();
      it != contour.end(); ++it) {
    fout << it->first << "\t" << it->second << "\n";
  }
  fout << "# 95% confidence region: x [m]\t y[m]  in local coordinates \n";
  for (std::vector<std::pair<double, double> >::const_iterator it = contour2.begin();
      it != contour2.end(); ++it) {
    fout << it->first << "\t" << it->second << "\n";
  }
  fout.close();
  return eSuccess;
}

ParameterLDFModel RdLDFFitter::PreLDFFit(int LDFModel, const std::vector<StationFitData>& vStationFitData,
                                         utl::Vector ShowerAxis, utl::Point CorePosition) {
  switch (LDFModel) {
    case 1: {
      TGraphErrors gLDF = TGraphErrors();
      for (unsigned int i = 0; i < vStationFitData.size(); ++i) {
        StationFitData data = vStationFitData[i];
        double distance = utl::RadioGeometryUtilities::GetDistanceToAxis(
            ShowerAxis, CorePosition, data.antennaPosition);
        gLDF.SetPoint(i, distance, data.signal);
        gLDF.SetPointError(i, 0., data.signalError);
      }
      TF1 fLDF = TF1("fLDF", "[0]*exp(-x/[1])");
      fLDF.SetParNames("A_0", "R_0");
      fLDF.SetParameters(0.1, 100);
      gLDF.Fit(&fLDF);

      ParameterLDFModel par;
      par.A = fLDF.GetParameter(0);
      par.R0 = fLDF.GetParameter(1);
      par.AError = fLDF.GetParError(0);
      par.R0Error = fLDF.GetParError(1);
      return par;
    }
    default: {
      WARNING("Using default LDF model, as chosen model is not implemented\n");
      TGraphErrors gLDF = TGraphErrors();
      for (unsigned int i = 0; i < vStationFitData.size(); ++i) {
        StationFitData data = vStationFitData[i];
        double distance = utl::RadioGeometryUtilities::GetDistanceToAxis(
            ShowerAxis, CorePosition, data.antennaPosition);
        gLDF.SetPoint(i, distance, data.signal);
        gLDF.SetPointError(i, 0., data.signalError);
      }
      TF1 fLDF = TF1("fLDF", "[0]*exp(-x/[1])");
      fLDF.SetParNames("A_0", "R_0");
      fLDF.SetParameters(0.1, 100);
      gLDF.Fit(&fLDF);

      ParameterLDFModel par;
      par.A = fLDF.GetParameter(0);
      par.R0 = fLDF.GetParameter(1);
      par.AError = fLDF.GetParError(0);
      par.R0Error = fLDF.GetParError(1);
      return par;
    }
  }
}

ParameterScintillatorLDFModel RdLDFFitter::PreScintillatorLDFFit(
    const std::vector<ScintillatorFitData>& vScintillatorData, utl::Vector ShowerAxis,
    utl::Point CorePosition)
{
  ParameterScintillatorLDFModel tmpFitData;
  TGraphErrors g = TGraphErrors();
  for (std::vector<ScintillatorFitData>::const_iterator sIt = vScintillatorData.begin(), end =
      vScintillatorData.end(); sIt != end; ++sIt) {
    double distanceToShowerAxis = utl::RadioGeometryUtilities::GetDistanceToAxis(
        ShowerAxis, CorePosition, sIt->scintillatorPosition);
    g.SetPoint(g.GetN(), distanceToShowerAxis, sIt->signal);
    g.SetPointError(g.GetN() - 1, 0, sqrt(sIt->signal));
  }
  TF1 *NKG_first_stage = new
      TF1("NKG_first_stage",
          "[0]*TMath::Gamma(4.5-[1])/(2.*TMath::Pi()*[2]*[2]*TMath::Gamma([1])*TMath::Gamma(4.5-(2.*[1])))*TMath::Power((x/[2]),([1]-2.))*TMath::Power((1+(x/[2])),[1]-4.5)",
          0, 1000);

  NKG_first_stage->SetParameter(0, 1. * 1e6);
  NKG_first_stage->SetParLimits(0, 100000, 100000000000);

  NKG_first_stage->FixParameter(1, 1.7);
  NKG_first_stage->SetParameter(2, 35);
  NKG_first_stage->SetParLimits(2, 5, 300);
  NKG_first_stage->FixParameter(2, 35);
  Int_t fitStatus = g.Fit(NKG_first_stage, "NM");
  tmpFitData.status = fitStatus;
  tmpFitData.N_e = NKG_first_stage->GetParameter(0);
  tmpFitData.N_eError = NKG_first_stage->GetParError(0);
  return tmpFitData;
}
} // namespace RdLDFFitter