RdWaveFit.cc

Go to the documentation of this file.

#include "RdWaveFit.h"

#include <cmath>

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

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

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/PhysicalConstants.h>
#include <utl/AxialVector.h>
#include <utl/Math.h>
#include <utl/AnalyticCoordinateTransformator.h>

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

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


#include <TMinuit.h>
#include <TRandom3.h>

#include <Math/Minimizer.h>
#include <Minuit2/Minuit2Minimizer.h>
#include "Chi2ForSphericalWaveFitVarC.h"

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

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

namespace RdWaveFit {

template<typename G>
inline string ToString(const G& g, const CoordinateSystemPtr cs)
{
  ostringstream os;
  os << '(' << g.GetX(cs) / meter << ", " << g.GetY(cs) / meter << ", " << g.GetZ(cs) / meter
     << ") [m]";
  return os.str();
}

RdWaveFit::RdWaveFit() :
    fInfoLevel(0),
    fMinuitOutput(0),
    fMinNumberOfStations(0),
    fAllowUnphysicalCosines(0),
    fFitType(0),
    fThetaMax(0),
    fRMax(0),
    fRhoMax(0),
    fNInitsTheta(0),
    fNInitsPhi(0),
    fSphericalWaveFit(0),
    fWaveFrontModel(0),
    currentWaveFrontModel(0),
    fSphericalWaveFitVarC(0),
    fComputeSigmaGamma(0),
    fSigmaGamma(0),
    fPenaltyTerm(0),
    fThetaVarC(0),
    fNumberOfStations(0),
    fEventTimeNanoSecond(0),
    fEventTimeSecond(0),
    ChiPWF(0),
    ChiSWF(0),
    ChiSWFVarC(0),
    ChiCWF(0),
    fUsedDirection(""),
    fUsedCorePosition("")
{
}

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

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

  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

  fMinuitOutput = (fInfoLevel >= eMinuit) ? true : false;

  topBranch.GetChild("minNumberOfStations").GetData(fMinNumberOfStations);
  topBranch.GetChild("fitType").GetData(fFitType);
  topBranch.GetChild("thetaMax").GetData(fThetaMax);
  topBranch.GetChild("rMax").GetData(fRMax);
  topBranch.GetChild("rhoMax").GetData(fRhoMax);
  topBranch.GetChild("nInitsTheta").GetData(fNInitsTheta);
  topBranch.GetChild("nInitsPhi").GetData(fNInitsPhi);
  topBranch.GetChild("waveFrontModel").GetData(fWaveFrontModel);
  topBranch.GetChild("computeSigmaGamma").GetData(fComputeSigmaGamma);
  topBranch.GetChild("sigmaGamma").GetData(fSigmaGamma);
  topBranch.GetChild("penaltyTerm").GetData(fPenaltyTerm);
  topBranch.GetChild("thetaVarC").GetData(fThetaVarC);
  topBranch.GetChild("UsedDirection").GetData(fUsedDirection);
  topBranch.GetChild("UsedCorePosition").GetData(fUsedCorePosition);

  if (fNInitsPhi < 1) {
    fNInitsPhi = 1;
  }
  if (fNInitsTheta < 1) {
    fNInitsTheta = 1;
  }

  ChiPWF = new Chi2ForPlaneWaveFit();
  ChiSWF = new Chi2ForSphericalWaveFit();
  ChiSWFVarC = new Chi2ForSphericalWaveFitVarC();
  ChiCWF = new Chi2ForConicalWaveFit();
  ChiSWFVarC->SetPenalty(fPenaltyTerm);
  if (!fComputeSigmaGamma) {
    ChiSWFVarC->SetSigmaGamma(fSigmaGamma);
  }

  return eSuccess;
}

VModule::ResultFlag RdWaveFit::Run(evt::Event& event)
{
  currentWaveFrontModel = fWaveFrontModel;
  if (fInfoLevel > 0) {
    INFO("RdWaveFit::Run()");
    switch (currentWaveFrontModel) {
      case ePlane:
        INFO("Plane Wave Fit");
        break;
      case eSpherical:
        INFO("Spherical Wave Fit");
        break;
      case eSphericalVarC:
        INFO("Spherical Wave Fit with variable speed of light");
        break;
      case eConical:
        INFO("Conical Wave Fit");
        break;
    }
  }

  // nothing to do if there is no REvent
  if (!event.HasREvent()) {
    OUT(eIntermediate) << "eContinueLoop: No REvent, so no reconstruction" << endl;
    return eContinueLoop;
  }

  REvent& rEvent = event.GetREvent();

  /* check if event was reconstructed before.   */
  if (!event.HasRecShower()) {
    if (fInfoLevel > eNone) {
      WARNING(
          "RdWaveFit: No RecShower found. Plead use RdPreWaveFitter before RdWaveFit... aborting");
    }
    return eBreakLoop;
  }
  if (!event.GetRecShower().HasRRecShower()) {
    if (fInfoLevel > eNone) {
      WARNING(
          "RdWaveFit: No RRecShower found. Plead use RdPreWaveFitter before RdWaveFit... aborting");
    }
    return eBreakLoop;
  }

  ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();
  fNumberOfStations = showerrrec.GetParameter(eNumberOfStationsWithPulseFound);

  if (fNumberOfStations < fMinNumberOfStations) {
    OUT(eFinal) << "eBreakLoop: Not enough stations. (Event contains only " << fNumberOfStations
                << " station(s). Min. number of stations = " << fMinNumberOfStations << ")" << endl;
    showerrrec.SetParameter(eFitSuccess, 0);
    return eBreakLoop;
  }

  // variable to indicate if this event was previously reconstructed
  bool reuseFit = false;
  // check if reuse fit?
  if (showerrrec.HasParameter(eRecStage)) {
    if (showerrrec.GetParameter(eRecStage) >= ShowerRRecData::kPlaneFit3d) {
      reuseFit = true;    
      OUT(eMinuit) << "Fit has been performed previously, setting reusefit = true" << endl;
    }
  }


  utl::Point baryCenter;
  VModule::ResultFlag resultBaryCenter = ComputeBaryCenter(rEvent, baryCenter);
  if(resultBaryCenter != eSuccess) {
    WARNING("BaryCenter could not be determined!");
    return resultBaryCenter;
  }

  // get reference coordinate system of detector (usually PampaAmarilla)
  const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
  showerrrec.SetParameter(eCoreX, baryCenter.GetX(referenceCS));
  showerrrec.SetParameter(eCoreY, baryCenter.GetY(referenceCS));
  showerrrec.SetParameter(eCoreZ, baryCenter.GetZ(referenceCS));


  // get coordinate origin depending on xml settings
  stringstream fMessage;
  // check if needed reconstructed quantities exist
  if(fUsedDirection == "SD" or fUsedCorePosition == "SD") {
    if(!event.GetRecShower().HasSRecShower()) {
      WARNING("Event has not SRecShower. Nothing will be done!!!");
      return eSuccess;
    }
  }
  if(fUsedDirection == "MC" or fUsedCorePosition == "MC") {
    if(!event.HasSimShower()) {
      WARNING("Event has no simulated shower. Nothing will be done!!!");
      return eSuccess;
    }
  }
  if(fUsedDirection == "FD" or fUsedCorePosition == "FD") {
    if(!event.GetRecShower().HasFRecShower()) {
      WARNING("Event has no FD reconstructed shower. Nothing will be done!!!");
      return eSuccess;
    }
  }

  // getting shower direction
  double zenith(0), azimuth(0);
  utl::Vector InitialShowerAxis = utl::Vector();
  if (fUsedDirection == "SD") {
    if (event.GetRecShower().HasSRecShower()) {
      InitialShowerAxis = event.GetRecShower().GetSRecShower().GetAxis();
      Point core = event.GetRecShower().GetSRecShower().GetCorePosition();
      utl::CoordinateSystemPtr coordinateSystem = LocalCoordinateSystem::Create(core);
      zenith = event.GetRecShower().GetSRecShower().GetAxis().GetTheta(coordinateSystem);
      azimuth = event.GetRecShower().GetSRecShower().GetAxis().GetPhi(coordinateSystem);
      OUT(eMinuit) << "Zenith = " << zenith << "\t azimuth = " << azimuth << endl;
    } else {
      WARNING("No SD Shower... aborting");
      return eSuccess;
    }
  }
  if (fUsedDirection == "MC") {
    cout << "using MC shower axis" << endl;
    InitialShowerAxis = -1 * event.GetSimShower().GetDirection();
  }
  if (fUsedDirection == "FD") {
    InitialShowerAxis = event.GetRecShower().GetFRecShower().GetAxis();
  }

  // getting shower core information
  if (fUsedCorePosition == "MC") {
    cout << "using MC shower core" << endl;
    fReferenceCore = event.GetSimShower().GetPosition();
  }
  if (fUsedCorePosition == "SD") {
    fReferenceCore = event.GetRecShower().GetSRecShower().GetCorePosition();
  }
  if (fUsedCorePosition == "Radio") {
    fReferenceCore = event.GetRecShower().GetRRecShower().GetCorePosition();
  }
  if (fUsedCorePosition == "FD") {
    fReferenceCore = event.GetRecShower().GetFRecShower().GetCorePosition();
  }
  if(fUsedCorePosition == "RadioCoordinateOrigin") {
    fReferenceCore = showerrrec.GetCoordinateOrigin();
  }
  fLocalCS = LocalCoordinateSystem::Create(fReferenceCore);

  // this will point to the final result
  FitParameters sFitResults;
  // set some standard values
  sFitResults.gamma = 1;
  sFitResults.gammaError = 0;
  sFitResults.r = 1;
  sFitResults.rError = 0;
  sFitResults.rho = 0;
  sFitResults.rhoError = 0;

  // same fit results from previous fit to the struct sFitResults
  if (reuseFit) {
    //   get values from RRecShower
    utl::Vector showeraxis = showerrrec.GetAxis();
    sFitResults.r = showeraxis.GetR(fLocalCS);
    sFitResults.theta = showeraxis.GetTheta(fLocalCS);
    sFitResults.phi = showeraxis.GetPhi(fLocalCS);
//       sFitResults.gamma = showerrrec.GetGammaC();
    if (showerrrec.HasParameter(eGamma)) {
      sFitResults.gamma = showerrrec.GetParameter(eGamma);
    } else {
      sFitResults.gamma = 1;
    }
  }

  // Call Wave Fit for different Models
  VModule::ResultFlag result = eFailure;
  switch (currentWaveFrontModel) {
    case ePlane:
      result = CallPlaneWaveFit(rEvent, sFitResults, reuseFit);
      break;
    case eSpherical:
      result = CallSphericalWaveFit(rEvent, sFitResults, reuseFit);
      break;
    case eSphericalVarC:
      result = CallSphericalWaveFitVarC(rEvent, sFitResults, reuseFit);
      break;
    case eConical:
      result = CallConicalWaveFit(rEvent, sFitResults, reuseFit);
      break;
  }

  // calculate and save residuals
  const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
  vector<utl::Vector> AntennaPositions;
  vector<double> AntennaTimes;
  vector<double> AntennaTimesError;
  //loop over all stations
  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    const Vector sPos = rDetector.GetStation(*sIt).GetPosition() - fReferenceCore;
    const StationRecData& sRec = sIt->GetRecData();
    AntennaPositions.push_back(sPos);
    const double t = sRec.GetParameter(revt::eSignalTime);
//        OUT(eObscure) << "SignalTime of Station " << sIt->GetId() << ": " << t << endl;
    AntennaTimes.push_back(t);
    AntennaTimesError.push_back(sRec.GetParameterError(eSignalTime));
  }
  const utl::Vector axis(sFitResults.r, sFitResults.theta, sFitResults.phi, fLocalCS, Vector::kSpherical);
  const std::vector<double> timeResiduals = GetTimeResiduals(AntennaPositions, AntennaTimes, axis,
                                                             currentWaveFrontModel);
  // delete time residuals of all stations
  for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
    StationRecData& sRec = sIt->GetRecData();
    sRec.DeleteParameter(eTimeResidual);
  }

  // set time residuals of signal stations
  int iStation = 0;
  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    StationRecData& sRec = sIt->GetRecData();
    sRec.SetParameter(eTimeResidual, timeResiduals[iStation]);
    sRec.SetParameterError(eTimeResidual, AntennaTimesError[iStation]);
    OUT(eObscure) << "TimeResidual of Station " << sIt->GetId() << ": " << timeResiduals[iStation] << " +/- "
        << AntennaTimesError[iStation] << " ns" << endl;
    ++iStation;
  }


  // break loop if the fit was not successful
  if (result != eSuccess) {
    showerrrec.SetParameter(eFitSuccess, 0);
    return result;
  }


  showerrrec.SetParameter(eShowerAxisX, axis.GetX(fLocalCS));
  showerrrec.SetParameter(eShowerAxisY, axis.GetY(fLocalCS));
  showerrrec.SetParameter(eShowerAxisZ, axis.GetZ(fLocalCS));

  // save Covariance matrix
  double tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(
      0, 0, sFitResults.r, sFitResults.theta, sFitResults.phi, sFitResults.c00, sFitResults.c11,
      sFitResults.c22, sFitResults.c01, sFitResults.c02, sFitResults.c12);
  showerrrec.SetParameterCovariance(eShowerAxisX, eShowerAxisX, tempCov);
  tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(1, 1,
                                                                                 sFitResults.r,
                                                                                 sFitResults.theta,
                                                                                 sFitResults.phi,
                                                                                 sFitResults.c00,
                                                                                 sFitResults.c11,
                                                                                 sFitResults.c22,
                                                                                 sFitResults.c01,
                                                                                 sFitResults.c02,
                                                                                 sFitResults.c12);
  showerrrec.SetParameterCovariance(eShowerAxisY, eShowerAxisY, tempCov);
  tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(2, 2,
                                                                                 sFitResults.r,
                                                                                 sFitResults.theta,
                                                                                 sFitResults.phi,
                                                                                 sFitResults.c00,
                                                                                 sFitResults.c11,
                                                                                 sFitResults.c22,
                                                                                 sFitResults.c01,
                                                                                 sFitResults.c02,
                                                                                 sFitResults.c12);
  showerrrec.SetParameterCovariance(eShowerAxisZ, eShowerAxisZ, tempCov);
  tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(0, 1,
                                                                                 sFitResults.r,
                                                                                 sFitResults.theta,
                                                                                 sFitResults.phi,
                                                                                 sFitResults.c00,
                                                                                 sFitResults.c11,
                                                                                 sFitResults.c22,
                                                                                 sFitResults.c01,
                                                                                 sFitResults.c02,
                                                                                 sFitResults.c12);
  showerrrec.SetParameterCovariance(eShowerAxisX, eShowerAxisY, tempCov);
  tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(0, 2,
                                                                                 sFitResults.r,
                                                                                 sFitResults.theta,
                                                                                 sFitResults.phi,
                                                                                 sFitResults.c00,
                                                                                 sFitResults.c11,
                                                                                 sFitResults.c22,
                                                                                 sFitResults.c01,
                                                                                 sFitResults.c02,
                                                                                 sFitResults.c12);
  showerrrec.SetParameterCovariance(eShowerAxisX, eShowerAxisZ, tempCov);
  tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(1, 2,
                                                                                 sFitResults.r,
                                                                                 sFitResults.theta,
                                                                                 sFitResults.phi,
                                                                                 sFitResults.c00,
                                                                                 sFitResults.c11,
                                                                                 sFitResults.c22,
                                                                                 sFitResults.c01,
                                                                                 sFitResults.c02,
                                                                                 sFitResults.c12);
  showerrrec.SetParameterCovariance(eShowerAxisY, eShowerAxisZ, tempCov);

  showerrrec.SetParameter(eGamma, sFitResults.gamma);
  showerrrec.SetParameterError(eGamma, sFitResults.gammaError);

  showerrrec.SetParameter(eWaveFitChi2, sFitResults.minChi);
  showerrrec.SetParameter(eWaveFitNDF, sFitResults.NDF);
  showerrrec.SetParameter(eRecStage, sFitResults.recStage);
  showerrrec.SetParameter(eWaveFrontModel,sFitResults.recStage);

  if (sFitResults.recStage >= ShowerRRecData::kConicalWaveFit) {
    showerrrec.SetParameter(eRho, sFitResults.rho);
    showerrrec.SetParameterError(eRho, sFitResults.rhoError);
  }
  
  // set signal arrival direction for individual stations

  for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
    Station& station = *sIt;
    
    // for spherical wave use direction to point source
    // therefore, emission point is transfered to local coordniate system of station
    if ( (currentWaveFrontModel == eSpherical) || (currentWaveFrontModel == eSphericalVarC) ) {
      // has to be implemented
      const CoordinateSystemPtr stationCS = rDetector.GetStation(station.GetId()).GetLocalCoordinateSystem();
      station.GetRecData().SetParameter(eSignalArrivalDirectionX, axis.GetX(stationCS));
      station.GetRecData().SetParameter(eSignalArrivalDirectionY, axis.GetY(stationCS));
      station.GetRecData().SetParameter(eSignalArrivalDirectionZ, axis.GetZ(stationCS));
    } else { 
      // temporary solution: use shower direction for all other models.
      // this is exact for plane wave, but typically leads to 1-2 degree error for conical or hyperbolical wavefront
      station.GetRecData().SetParameter(eSignalArrivalDirectionX, axis.GetX(fLocalCS));
      station.GetRecData().SetParameter(eSignalArrivalDirectionY, axis.GetY(fLocalCS));
      station.GetRecData().SetParameter(eSignalArrivalDirectionZ, axis.GetZ(fLocalCS));
    }
  }

  showerrrec.SetParameter(eFitSuccess, 1);

  OUT(eFinal) << "RecStage: " << showerrrec.GetParameter(eRecStage) << NL
      << " axis in cartesian coord.= (" << axis.GetX(fLocalCS) << ", " << axis.GetY(fLocalCS)
      << ", " << axis.GetZ(fLocalCS) << ")" << NL << " axis in spherical coord.= ("
      << axis.GetR(fLocalCS) << ", " << axis.GetTheta(fLocalCS) << ", " << axis.GetPhi(fLocalCS)
      << ")" << NL << "theta = " << sFitResults.theta / utl::degree << " +/- "
      << sFitResults.thetaError / utl::degree << "\t phi = " << sFitResults.phi / utl::degree
      << " +/- " << sFitResults.phiError / utl::degree << NL << "Radius = "
      << sFitResults.r / utl::m << "m +/- " << sFitResults.rError / utl::m << "m" << NL
      << "Gamma = " << sFitResults.gamma << " +/- " << sFitResults.gammaError << NL << endl;
  OUT(eIntermediate) << "Chi2/NDF = " << sFitResults.minChi << "/" << sFitResults.NDF << NL << endl;

  return eSuccess;
}

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

VModule::ResultFlag RdWaveFit::PlaneWaveFit(revt::REvent& rEvent, FitParameters& fit,
                                            int type) const
{
  if (fNumberOfStations < 3) {
    OUT(eBreakLoop) << "ERROR: At least 3 Stations are needed for a plane wave fit" << endl;
    return eBreakLoop;
  }

  ROOT::Minuit2::Minuit2Minimizer* minimizer = 0;
  if (type == 0) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kMigrad);
  } else if (type == 1) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kSimplex);
  } else if (type == 2) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kCombined);
  } else if (type == 3) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kScan);
  }

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

  vector<utl::Vector> AntennaPositions;
  vector<double> AntennaTimes;
  vector<double> AntennaTimesError;
  //loop over all stations
  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    const Vector sPos = rDetector.GetStation(*sIt).GetPosition() - fReferenceCore;
    const StationRecData& sRec = sIt->GetRecData();
    AntennaPositions.push_back(sPos);
    const double t = sRec.GetParameter(revt::eSignalTime);
//        OUT(eObscure) << "SignalTime of Station " << sIt->GetId() << ": " << t << endl;
    AntennaTimes.push_back(t);
    AntennaTimesError.push_back(sRec.GetParameterError(eSignalTime));
  }
  fit.nStations = AntennaTimes.size();

// Give the Antenna Positions and Times to objective function
  ChiPWF->Set(AntennaPositions, AntennaTimes, AntennaTimesError, fLocalCS);
  minimizer->SetFunction(*ChiPWF);

// Set variables
  minimizer->SetLimitedVariable(0, "theta", fit.theta, 0.1, 0., fThetaMax);
  minimizer->SetVariable(1, "phi", fit.phi, 0.1);

// THE FIT:

  minimizer->SetPrintLevel(0);
  if (!minimizer->Minimize()) {
    fit.status = minimizer->Status();
    OUT(eObscure) << "Fit failed: FitStatus  = " << minimizer->Status() << endl;
    return eBreakLoop;
  }

// Save Results
  fit.status = minimizer->Status();
  fit.phi = NormalizeAngle(minimizer->X()[1]);
  fit.theta = NormalizeAngle(minimizer->X()[0]);
  fit.phiError = minimizer->Errors()[1];
  fit.thetaError = minimizer->Errors()[0];

  fit.nCalls = minimizer->NCalls();
  fit.minChi = minimizer->MinValue();
  fit.recStage = ShowerRRecData::kPlaneFit3d;

// save covariance matrix
  fit.c00 = 0;
  fit.c11 = minimizer->CovMatrix(0, 0);
  fit.c22 = minimizer->CovMatrix(1, 1);
  fit.c01 = 0;
  fit.c02 = 0;
  fit.c12 = minimizer->CovMatrix(0, 1);

  fit.NDF = fit.nStations - 3;

  /* the plane wave fit can be performed with just 3 stations. The number of degrees
   * of freedom is 0, though. To be able to compute chiSquare/NDF the NDF is set to 1.
   */
  if (fit.nStations < 4) {
    fit.chiSquareNDF = minimizer->MinValue() / (fit.nStations - 2);
  } else {
    fit.chiSquareNDF = minimizer->MinValue() / (fit.nStations - 3);
  }

// save Time Residual of Stations
  const utl::Vector axis(1. * utl::m, fit.theta, fit.phi, fLocalCS, Vector::kSpherical);
  const std::vector<double> timeResiduals = GetTimeResiduals(AntennaPositions, AntennaTimes, axis,
                                                             ePlane);
//loop over all stations
  int iStation = 0;
  for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
    StationRecData& sRec = sIt->GetRecData();
    sRec.DeleteParameter(eTimeResidual);
  }

  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    StationRecData& sRec = sIt->GetRecData();
    sRec.SetParameter(eTimeResidual, timeResiduals[iStation]);
    sRec.SetParameterError(eTimeResidual, AntennaTimesError[iStation]);
//        OUT(eObscure) << "TimeResidual of Station " << sIt->GetId() << ": "
//            << timeResiduals[iStation] << " +/- " << AntennaTimesError[iStation] << " ns" << endl;
    ++iStation;
  }

  delete minimizer;

  return eSuccess;
}

VModule::ResultFlag RdWaveFit::SphericalWaveFit(revt::REvent& rEvent, FitParameters& fit,
                                                int type) const
{

  if (fNumberOfStations < 4) {
    OUT(eBreakLoop) << "ERROR: At least 4 Stations are needed for a spherical wave fit" << endl;
    return eBreakLoop;
  }
  // Initialize the Minuit2 Minimizer
  ROOT::Minuit2::Minuit2Minimizer* minimizer = 0;
  if (type == 0) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kMigrad);
  } else if (type == 1) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kSimplex);
  } else if (type == 2) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kCombined);
  } else if (type == 3) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kScan);
  }

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

  vector<utl::Vector> AntennaPositions;
  vector<double> AntennaTimes;
  vector<double> AntennaTimesError;

  //loop over all stations
  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    const Vector sPos = rDetector.GetStation(*sIt).GetPosition() - fReferenceCore;

    const StationRecData& sRec = sIt->GetRecData();
    AntennaPositions.push_back(sPos);
    const double t = sRec.GetParameter(revt::eSignalTime);
    OUT(eObscure) << "SignalTime of Station " << sIt->GetId() << ": " << t << endl;
    AntennaTimes.push_back(t);
    AntennaTimesError.push_back(sRec.GetParameterError(eSignalTime));
  }
  fit.nStations = AntennaTimes.size();
  if (AntennaPositions.size() < 4) {
    OUT(eBreakLoop) << "Event has less than 3 Stations, can't perform spherical wave fit. " << endl;
    return eBreakLoop;
  }

  // Give the Antenna Positions and Times to objective function
  ChiSWF->Set(AntennaPositions, AntennaTimes, AntennaTimesError, fLocalCS);
  minimizer->SetFunction(*ChiSWF);

  minimizer->SetLimitedVariable(0, "r", fit.r, 10., 0., fRMax);
  minimizer->SetLimitedVariable(1, "theta", fit.theta, 0.1, 0., fThetaMax);
  minimizer->SetVariable(2, "phi", fit.phi, 0.1);

  // THE FIT:
  minimizer->SetPrintLevel(0);
  if (!minimizer->Minimize()) {
    fit.status = minimizer->Status();
    return eBreakLoop;
  }

  // Save Results
  fit.status = minimizer->Status();
  fit.r = minimizer->X()[0];
  fit.phi = NormalizeAngle(minimizer->X()[2]);
  fit.theta = NormalizeAngle(minimizer->X()[1]);
  fit.rError = minimizer->Errors()[0];
  fit.phiError = minimizer->Errors()[2];
  fit.thetaError = minimizer->Errors()[1];

  fit.nCalls = minimizer->NCalls();
  fit.minChi = minimizer->MinValue();
  fit.recStage = ShowerRRecData::kSphericalWaveFit;

  // save covariance matrix
  fit.c00 = minimizer->CovMatrix(0, 0);
  fit.c11 = minimizer->CovMatrix(1, 1);
  fit.c22 = minimizer->CovMatrix(2, 2);
  fit.c01 = minimizer->CovMatrix(0, 1);
  fit.c02 = minimizer->CovMatrix(0, 2);
  fit.c12 = minimizer->CovMatrix(1, 2);

  fit.NDF = fit.nStations - 4;

  /* the spherical wave fit can be performed with just 4 stations. The number of degrees
   * of freedom is 0, though. To be able to compute chiSquare/NDF the NDF is set to 1.
   */
  if (fit.nStations < 5) {
    fit.chiSquareNDF = minimizer->MinValue() / (fit.nStations - 3);
  } else {
    fit.chiSquareNDF = minimizer->MinValue() / (fit.nStations - 4);
  }

  // save Time Residual of Stations
  const utl::Vector axis(fit.r, fit.theta, fit.phi, fLocalCS, Vector::kSpherical);
  const std::vector<double> timeResiduals = GetTimeResiduals(AntennaPositions, AntennaTimes, axis,
                                                             eSpherical);

  for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
    StationRecData& sRec = sIt->GetRecData();
    sRec.DeleteParameter(eTimeResidual);
  }

  //loop over all stations
  int iStation = 0;
  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    StationRecData& sRec = sIt->GetRecData();
    sRec.SetParameter(eTimeResidual, timeResiduals[iStation]);
    sRec.SetParameterError(eTimeResidual, AntennaTimesError[iStation]);
//        OUT(eObscure) << "TimeResidual of Station " << sIt->GetId() << ": " << timeResiduals[iStation] << " +/- "
//            << AntennaTimesError[iStation] << " ns" << endl;
    ++iStation;
  }

  delete minimizer;

  return eSuccess;
}

VModule::ResultFlag RdWaveFit::SphericalWaveFitVarC(revt::REvent& rEvent, FitParameters& fit,
                                                    int type) const
{
  if (fNumberOfStations < 5) {
    OUT(eBreakLoop)
        << "ERROR: At least 5 Stations are needed for a spherical wave fit with variable speed of light"
        << endl;
    return eBreakLoop;
  }
  // Initialize the Minuit2 Minimizer
  ROOT::Minuit2::Minuit2Minimizer* minimizer = 0;
  if (type == 0) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kMigrad);
  } else if (type == 1) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kSimplex);
  } else if (type == 2) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kCombined);
  } else if (type == 3) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kScan);
  }

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

  vector<utl::Vector> AntennaPositions;
  vector<double> AntennaTimes;
  vector<double> AntennaTimesError;

  //loop over all stations
  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    const Vector sPos = rDetector.GetStation(*sIt).GetPosition() - fReferenceCore;

    const StationRecData& sRec = sIt->GetRecData();
    AntennaPositions.push_back(sPos);
    const double t = sRec.GetParameter(revt::eSignalTime);
    OUT(eObscure) << "SignalTime of Station " << sIt->GetId() << ": " << t << endl;
    AntennaTimes.push_back(t);
    AntennaTimesError.push_back(sRec.GetParameterError(eSignalTime));
  }
  fit.nStations = AntennaTimes.size();

  // Give the Antenna Positions and Times to objective function
  ChiSWFVarC->Set(AntennaPositions, AntennaTimes, AntennaTimesError, fLocalCS);
  minimizer->SetFunction(*ChiSWFVarC);

  minimizer->SetLimitedVariable(0, "r", fit.r, 10., 0., fRMax);
  minimizer->SetLimitedVariable(1, "theta", fit.theta, 0.1, 0., fThetaMax);
  minimizer->SetVariable(2, "phi", fit.phi, 0.1);
  minimizer->SetLowerLimitedVariable(3, "gamma", fit.gamma, 0.01, 0.);

  // THE FIT:
  minimizer->SetPrintLevel(0);
  if (!minimizer->Minimize()) {
    fit.status = minimizer->Status();
    return eBreakLoop;
  }

  // Save Results
  fit.status = minimizer->Status();
  fit.r = minimizer->X()[0];
  fit.phi = NormalizeAngle(minimizer->X()[2]);
  fit.theta = NormalizeAngle(minimizer->X()[1]);
  fit.rError = minimizer->Errors()[0];
  fit.phiError = minimizer->Errors()[2];
  fit.thetaError = minimizer->Errors()[1];
  fit.gamma = minimizer->X()[3];
  fit.gammaError = minimizer->Errors()[3];

  fit.nCalls = minimizer->NCalls();
  fit.minChi = minimizer->MinValue();
  fit.recStage = ShowerRRecData::kSphericalWaveFitVarC;

  // save covariance matrix
  fit.c00 = minimizer->CovMatrix(0, 0);
  fit.c11 = minimizer->CovMatrix(1, 1);
  fit.c22 = minimizer->CovMatrix(2, 2);
  fit.c01 = minimizer->CovMatrix(0, 1);
  fit.c02 = minimizer->CovMatrix(0, 2);
  fit.c12 = minimizer->CovMatrix(1, 2);

  fit.NDF = fit.nStations - 5;

  /* the spherical wave fit VarC can be performed with just 5 stations. The number of degrees
   * of freedom is 0, though. To be able to compute chiSquare/NDF the NDF is set to 1.
   */
  if (fit.nStations < 6) {
    fit.chiSquareNDF = minimizer->MinValue() / (fit.nStations - 4);
  } else {
    fit.chiSquareNDF = minimizer->MinValue() / (fit.nStations - 5);
  }
  // save Time Residual of Stations
  const utl::Vector axis(fit.r, fit.theta, fit.phi, fLocalCS, Vector::kSpherical);
  const std::vector<double> timeResiduals = GetTimeResiduals(AntennaPositions, AntennaTimes, axis,
                                                             eSphericalVarC, fit.gamma);
  //loop over all stations
  int iStation = 0;

  for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
    StationRecData& sRec = sIt->GetRecData();
    sRec.DeleteParameter(eTimeResidual);
  }

  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {

    StationRecData& sRec = sIt->GetRecData();
    sRec.SetParameter(eTimeResidual, timeResiduals[iStation]);
    sRec.SetParameterError(eTimeResidual, AntennaTimesError[iStation]);
//        OUT(eObscure) << "TimeResidual of Station " << sIt->GetId() << ": " << timeResiduals[iStation] << " +/- "
//            << AntennaTimesError[iStation] << " ns" << endl;
    ++iStation;
  }

  delete minimizer;

  return eSuccess;
}

VModule::ResultFlag RdWaveFit::ConicalWaveFit(revt::REvent& rEvent, FitParameters& fit,
                                              int type) const
{
  if (fNumberOfStations < 4) {
    OUT(eBreakLoop) << "ERROR: At least 4 stations are needed for a conical wave fit" << endl;
    return eBreakLoop;
  }
  // Initiate the Minuit2 Minimizer
  ROOT::Minuit2::Minuit2Minimizer* minimizer = 0;
  if (type == 0) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kMigrad);
  } else if (type == 1) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kSimplex);
  } else if (type == 2) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kCombined);
  } else if (type == 3) {
    minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kScan);
  }

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

  vector<utl::Vector> AntennaPositions;
  vector<double> AntennaTimes;
  vector<double> AntennaTimesError;

  // loop over all stations
  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {
    const Vector sPos = rDetector.GetStation(*sIt).GetPosition() - fReferenceCore;
    const StationRecData& sRec = sIt->GetRecData();
    AntennaPositions.push_back(sPos);
    const double t = sRec.GetParameter(revt::eSignalTime);
//                        OUT(eObscure) << "SignalTime of Station " << sIt->GetId() << ": " << t << endl;
    AntennaTimes.push_back(t);
    AntennaTimesError.push_back(sRec.GetParameterError(eSignalTime));
  }
  fit.nStations = AntennaTimes.size();
  // Give the Antenna Position and Times to objective function
  ChiCWF->Set(AntennaPositions, AntennaTimes, AntennaTimesError, fLocalCS);
  minimizer->SetFunction(*ChiCWF);
  minimizer->SetVariable(0, "rho", fit.rho, 0.01);
  //minimizer->SetLimitedVariable(0, "rho", fit.rho, 0.01, 0., fRhoMax); // for non-negative rho angles
  /*minimizer->SetLimitedVariable(1, "theta", fit.theta, 0.1, 0., fThetaMax);
  minimizer->SetVariable(2, "phi", fit.phi, 0.1);*/
  // Set theta and phi fixed on values from spherical wave fit varC
  minimizer->SetFixedVariable(1, "theta", fit.theta);
  minimizer->SetFixedVariable(2, "phi", fit.phi);


  // THE FIT:
  OUT(eMinuit) << "vor minimization" << endl;
  minimizer->SetPrintLevel(1);
  if (!minimizer->Minimize()) {
    fit.status = minimizer->Status();
    return eBreakLoop;
  }

  //Save Results
  fit.status = minimizer->Status();
  fit.rho = NormalizeAngle(minimizer->X()[0]);
  fit.phi = NormalizeAngle(minimizer->X()[2]);
  fit.theta = NormalizeAngle(minimizer->X()[1]);
  fit.rhoError = minimizer->Errors()[0];
  fit.phiError = minimizer->Errors()[2];
  fit.thetaError = minimizer->Errors()[1];

  OUT(eMinuit) << "reconstructed rho = " << fit.rho / utl::degree << " +/- "
               << fit.rhoError / utl::degree << " deg" << endl;

  fit.nCalls = minimizer->NCalls();
  fit.minChi = minimizer->MinValue();
  fit.recStage = ShowerRRecData::kConicalWaveFit;

//        OUT(eMinuit) << "nCalls: " << fit.nCalls << endl;

  // save covariance matrix
  fit.c00 = 0;  // set to zero because it's the covariance between radius and theta/phi. And radius doesn't exist in conical wave fit
  fit.c11 = minimizer->CovMatrix(1, 1);
  fit.c22 = minimizer->CovMatrix(2, 2);
  fit.c01 = 0;  // set to zero because it's the covariance between radius and theta/phi. And radius doesn't exist in conical wave fit
  fit.c02 = 0;  // set to zero because it's the covariance between radius and theta/phi. And radius doesn't exist in conical wave fit
  fit.c12 = minimizer->CovMatrix(1, 2);

  fit.NDF = fit.nStations - 4;

  /* the conical wave fit can be performed with just 4 stations. The number of degrees
   * of freedom is 0, though. To be able to compute chiSquare/NDF the NDF is set to 1.
   */
  if (fit.nStations < 5) {
    fit.chiSquareNDF = minimizer->MinValue() / (fit.nStations - 3);
  } else {
    fit.chiSquareNDF = minimizer->MinValue() / (fit.nStations - 4);
  }
  // save Time Residual of Stations
  const utl::Vector axis(1, fit.theta, fit.phi, fLocalCS, Vector::kSpherical);
  const std::vector<double> timeResiduals = GetTimeResiduals(AntennaPositions, AntennaTimes, axis,
                                                             eConical, 1, fit.rho);
  //loop over all stations
  int iStation = 0;

  for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
    StationRecData& sRec = sIt->GetRecData();
    sRec.DeleteParameter(eTimeResidual);
  }

  for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
      sIt != rEvent.SignalStationsEnd(); ++sIt) {

    StationRecData& sRec = sIt->GetRecData();
    sRec.SetParameter(eTimeResidual, timeResiduals[iStation]);
    sRec.SetParameterError(eTimeResidual, AntennaTimesError[iStation]);
//                OUT(eObscure) << "TimeResidual of Station " << sIt->GetId() << ": " << timeResiduals[iStation] << " +/- "
//                    << AntennaTimesError[iStation] << " ns" << endl;
    iStation++;
  }

  delete minimizer;

  return eSuccess;
}

VModule::ResultFlag RdWaveFit::ScanPWF(revt::REvent& rEvent, FitParameters& sFitResults) const
{
  // this will store init & intermediate Parameters
  FitParameters sFitParameters;

  // necessary to save standard values Radius and Gamma
  sFitParameters = sFitResults;

  // necessary to abort if fit fails for all initial values
  sFitResults.status = 1;

  // set initial fit parameters
  double stepTheta = 90. / (fNInitsTheta + 1) * utl::degree;
  double stepPhi = 360. / (fNInitsPhi) * utl::degree;
  double theta(0.);
  double phi = -180. * utl::degree;

  OUT(eObscure) << "Stepwidth for initial values: stepTheta = " << stepTheta / utl::deg
                << " stepPhi = " << stepPhi / utl::deg << endl;

  // call PlaneWaveFit for different init values
  bool first_fit = true;
  for (int iTheta = 0; iTheta < fNInitsTheta; ++iTheta) {
    theta = stepTheta * (iTheta + 1);
    for (int iPhi = 0; iPhi < fNInitsPhi; ++iPhi) {
      phi = -180. * utl::degree + stepPhi * iPhi;
      OUT(eObscure) << " theta = " << theta / utl::deg << " deg\tphi = " << phi / utl::deg << " deg"
                    << endl;

      sFitParameters.theta = theta;
      sFitParameters.phi = phi;

      if (PlaneWaveFit(rEvent, sFitParameters, 1) != eSuccess)                          // call Fit in kSimplex mode
          {
        OUT(eObscure) << "PlaneWaveFit for theta = " << theta / utl::deg << " phi = "
                      << phi / utl::deg
                      << " in kSimplex was not successful: skip to next initial value" << endl;
        continue;
      }
      if (PlaneWaveFit(rEvent, sFitParameters, 0) != eSuccess) {
        OUT(eObscure) << "PlaneWaveFit for theta = " << theta / utl::deg << " phi = "
                      << phi / utl::deg
                      << " in kMigrad was not successful: skip to next initial value" << endl;
        continue;
      }

      OUT(eObscure) << "sFitParameters.status = " << sFitParameters.status << endl;

      // save Fit Results if better as previous fit or if its the first successful fit
      if (sFitParameters.minChi < sFitResults.minChi || first_fit) {
        first_fit = false;
        sFitResults = sFitParameters;
      }

    }
  }

  // break Function if reconstruction is not successful
  if (sFitResults.status != 0) {
    OUT(eBreakLoop) << "ERROR: reconstruction was not successful" << endl;
    return eBreakLoop;
  }

  sFitResults.corrected = false;
  if (sFitResults.theta > TMath::Pi() / 2.)  // flip minimum about horizon, check for better minimum
      {
    sFitParameters = sFitResults;  // buffer best result, so that this result won't be lost due to the new fit

    sFitParameters.theta = TMath::Pi() - sFitParameters.theta;

    if (PlaneWaveFit(rEvent, sFitParameters, 1) != eSuccess)                            // call Fit in kSimplex mode
        {
      OUT(eObscure) << "Horizon Correction: PlaneWaveFit for theta = " << sFitParameters.theta
                    << " in kSimplex was not successful: calling again in kMigrad mode" << endl;
    }

    if (PlaneWaveFit(rEvent, sFitParameters, 0) == eSuccess)            // call Fit in kMigrad mode
        {
      if (sFitParameters.minChi < sFitResults.minChi) {
        OUT(eObscure) << "Horizon corrected! " << endl;
        sFitParameters.corrected = true;
        sFitResults = sFitParameters;
      }
    } else {
      OUT(eObscure) << "Horizon Correction: PlaneWaveFit for theta "
                    << sFitParameters.theta / utl::deg << " in kMigrad was not successful" << endl;
    }
  }

  return eSuccess;
}

VModule::ResultFlag RdWaveFit::ScanSWF(revt::REvent& rEvent, FitParameters& sFitResults) const
{
  // this will store init & intermediate Parameters
  FitParameters sFitParameters;

  //save fit results from PlaneWaveFit as initial parameters for SphericalWaveFit
  sFitParameters = sFitResults;

  // initial values for the radius
  double rInits[] = { 0., 300., 1000., 3000., 10000. };
  const int nRInits = sizeof(rInits) / sizeof(rInits[0]);

  // set initial fit status to failure
  sFitResults.status = 1;

  // flag to identify the first successful fit
  bool first_fit = true;
  for (int i = 0; i < nRInits; ++i) {
    sFitParameters.r = rInits[i];

    if (SphericalWaveFit(rEvent, sFitParameters, 1) != eSuccess)                                // call Fit in kSimplex mode
        {
      OUT(eObscure) << "SphericalWaveFit for r = " << rInits[i]
                    << " in kSimplex was not successful: skip to next initial value" << endl;
      continue;
    }
    if (SphericalWaveFit(rEvent, sFitParameters, 0) != eSuccess) {
      OUT(eObscure) << "SphericalWaveFit for r = " << rInits[i]
                    << " in kMigrad was not successful: skip to next initial value" << endl;

      continue;
    }

    // save Fit Results if better as previous fit or if its the first fit

    if (sFitParameters.minChi < sFitResults.minChi || first_fit) {
      first_fit = false;
      sFitResults = sFitParameters;
    }
  }

  // break Function if reconstruction is not successful
  if (sFitResults.status > 0) {
    OUT(eBreakLoop) << "ERROR: reconstruction was not successful" << endl;
    return eBreakLoop;
  }

  sFitResults.corrected = false;
  if (sFitResults.theta > TMath::Pi() / 2.)  // flip minimum about horizon, check for better minimum
      {
    sFitParameters = sFitResults;  // buffer best result, so that this result won't be lost due to the new fit

    sFitParameters.theta = TMath::Pi() - sFitParameters.theta;

    if (SphericalWaveFit(rEvent, sFitParameters, 1) != eSuccess)                                // call Fit in kSimplex mode
        {
      OUT(eObscure) << "Horizon Correction: SphericalWaveFit for theta = " << sFitParameters.theta
                    << " in kSimplex was not successful: calling again in kMigrad mode" << endl;
    }

    if (SphericalWaveFit(rEvent, sFitParameters, 0) == eSuccess)                // call Fit in kMigrad mode
        {
      if (sFitParameters.minChi < sFitResults.minChi) {
        OUT(eObscure) << "Horizon corrected! " << endl;
        sFitParameters.corrected = true;
        sFitResults = sFitParameters;
      }
    } else {
      OUT(eObscure) << "Horizon Correction: SphericalWaveFit for theta " << sFitParameters.theta
                    << " in kMigrad was not successful" << endl;
    }
  }

  return eSuccess;
}

VModule::ResultFlag RdWaveFit::ScanSWFVarC(revt::REvent& rEvent, FitParameters& sFitResults) const
{
  // this will store init & intermediate Parameters
  FitParameters sFitParameters;

  //save fit results from SphericalWaveFit as initial parameters for SphericalWaveFitVarC
  sFitParameters = sFitResults;

  sFitResults.status = eFailure;

  if (SphericalWaveFitVarC(rEvent, sFitParameters, 1) != eSuccess)                      // call Fit in kSimplex mode
      {
    OUT(eBreakLoop) << "SphericalWaveFitVarC in kSimplexMode was not successful" << endl;
    return eBreakLoop;
  }
  if (SphericalWaveFitVarC(rEvent, sFitParameters, 0) != eSuccess) {
    OUT(eBreakLoop) << "SphericalWaveFitVarC in kMigradMode was not successful" << endl;
    return eBreakLoop;
  }

  sFitResults = sFitParameters;

  sFitResults.corrected = false;
  if (sFitResults.theta > TMath::Pi() / 2.)  // flip minimum about horizon, check for better minimum
      {
    sFitParameters = sFitResults;  // buffer best result, so that this result won't be lost due to the new fit

    sFitParameters.theta = TMath::Pi() - sFitParameters.theta;

    SphericalWaveFitVarC(rEvent, sFitParameters, 1);            // call Fit in kSimplex mode

    if (SphericalWaveFitVarC(rEvent, sFitParameters, 0) == eSuccess)            // call Fit in kMigrad mode
        {
      if (sFitParameters.minChi < sFitResults.minChi) {
        sFitParameters.corrected = true;
        sFitResults = sFitParameters;
      }
    }
  }

  return eSuccess;
}

VModule::ResultFlag RdWaveFit::ScanCWF(revt::REvent& rEvent, FitParameters& sFitResults) const
{
  OUT(eMinuit) << "ScanCWF" << endl;
  //this will store init & intermediate Parameters
  FitParameters sFitParameters;

  //save fit results from ConicalWaveFit as initial parameters for ConicalWaveFit
  sFitParameters = sFitResults;

  OUT(eMinuit) << "start values for conical wave fit : zenith = "
               << sFitParameters.theta / utl::degree << "deg\t azimuth = "
               << sFitParameters.phi / utl::degree << " deg" << endl;

  // set initial fit status to failure
  sFitResults.status = eFailure;

  double rhoInits[] = { 0., 0.5*utl::degree, 1.0*utl::degree, 1.5*utl::degree, 2.0*utl::degree, 3.0*utl::degree };
  const int nRhoInits = sizeof(rhoInits) / sizeof(rhoInits[0]);

  bool first_fit = true;
  for (int i = 0; i < nRhoInits; i++) {
    OUT(eMinuit) << "conical wave fit for rho = " << rhoInits[i] / utl::degree << " deg" << endl;

    sFitParameters.rho = rhoInits[i];

    if (ConicalWaveFit(rEvent, sFitParameters, 1) != eSuccess) {
      OUT(eBreakLoop) << "ConicalWaveFit in kSimplexMode was not successful" << endl;
      continue;
    }
    if (ConicalWaveFit(rEvent, sFitParameters, 0) != eSuccess) {
      OUT(eObscure) << "ConicalWaveFit for rho = " << sFitParameters.rho
                    << " in kMigrad was not successful: skip to next initial value" << endl;

      continue;
    }

    // save Fit Results if better as previous fit or if its the first fit

    if (sFitParameters.minChi < sFitResults.minChi || first_fit) {
      first_fit = false;
      sFitResults = sFitParameters;
    }
  }

  // break Function if reconstruction is not successful
  if (sFitResults.status != 0) {
    OUT(eBreakLoop) << "ERROR: reconstruction was not successful" << endl;
    return eBreakLoop;
  }

  sFitResults.corrected = false;
  if (sFitResults.theta > TMath::Pi() / 2.)  // flip minimum about horizon, check for better minimum
      {
    sFitParameters = sFitResults;  // buffer best result, so that this result won't be lost due to the new fit

    sFitParameters.theta = TMath::Pi() - sFitParameters.theta;

    if (ConicalWaveFit(rEvent, sFitParameters, 1) != eSuccess)                          // call Fit in kSimplex mode
        {
      OUT(eObscure) << "Horizon Correction: ConicalWaveFit for theta = " << sFitParameters.theta
                    << " in kSimplex was not successful: calling again in kMigrad mode" << endl;
    }

    if (ConicalWaveFit(rEvent, sFitParameters, 0) == eSuccess)          // call Fit in kMigrad mode
        {
      if (sFitParameters.minChi < sFitResults.minChi) {
        OUT(eObscure) << "Horizon corrected! " << endl;
        sFitParameters.corrected = true;
        sFitResults = sFitParameters;
      }
    } else {
      OUT(eObscure) << "Horizon Correction: ConicalWaveFit for theta " << sFitParameters.theta / utl::deg
                    << " in kMigrad was not successful" << endl;
    }
  }

  return eSuccess;
}

VModule::ResultFlag RdWaveFit::CallPlaneWaveFit(revt::REvent& rEvent, FitParameters& sFitResults,
                                                bool reuseFit) const
{

  VModule::ResultFlag result;

  if (fNumberOfStations < 3) {
    OUT(eBreakLoop) << "Not enough Stations" << endl;
    return eBreakLoop;
  }
  // for the first time
  if (!reuseFit) {
    return ScanPWF(rEvent, sFitResults);
  } else  // for a reuse fit
  {
    // call first in kSimplex
    result = PlaneWaveFit(rEvent, sFitResults, 1);
    if (result != eSuccess) {
      WARNING("reuse plane wave fit in kSimplex mode was not successful");
      return result;
    }

    // call in kMigrad
    result = PlaneWaveFit(rEvent, sFitResults, 0);
    return result;
  }
}

VModule::ResultFlag RdWaveFit::CallSphericalWaveFit(revt::REvent& rEvent,
                                                    FitParameters& sFitResults, bool reuseFit)
{

  VModule::ResultFlag result;
  // if event contains only 4 Stations perform just a plane wave fit
  if (fNumberOfStations < 4) {
    if (fInfoLevel > 0) {
      INFO(
          "Event contains less than 4 Stations: Only the PlaneWaveFit will be performed. Setting WaveFrontModel to ePlane");
    }
    currentWaveFrontModel = ePlane;
    return CallPlaneWaveFit(rEvent, sFitResults, reuseFit);
  }

  if (!reuseFit) {
    // perform a plane wave fit first
    result = CallPlaneWaveFit(rEvent, sFitResults, reuseFit);
    // exit function in case of a failure
    if (result != eSuccess) {
      return result;
    }

    return ScanSWF(rEvent, sFitResults);
  } else  // for a reuse fit
  {
    return SphericalWaveFit(rEvent, sFitResults, 0);
  }
}

VModule::ResultFlag RdWaveFit::CallSphericalWaveFitVarC(revt::REvent& rEvent,
                                                        FitParameters& sFitResults, bool reuseFit)
{

  VModule::ResultFlag result;

  if (fNumberOfStations < 4) {
    if (fInfoLevel > 0) {
      INFO(
          "Event contains less than 4 Stations: Only the PlaneWaveFit will be performed. Setting WaveFrontModel to ePlane");
    }
    currentWaveFrontModel = ePlane;
    return CallPlaneWaveFit(rEvent, sFitResults, reuseFit);
  }
  if (fNumberOfStations == 4) {
    if (fInfoLevel > 0) {
      INFO(
          "Event contains 4 Stations: Only the SphericalWaveFit will be performed. Setting WaveFrontModel to eSpherical");
    }
    currentWaveFrontModel = eSpherical;
    return CallSphericalWaveFit(rEvent, sFitResults, reuseFit);
  }

  if (!reuseFit) {
    // perform a plane wave fit first
    result = CallPlaneWaveFit(rEvent, sFitResults, reuseFit);
    // exit function in case of a failure
    if (result != eSuccess) {
      return result;
    }

    if (sFitResults.theta > fThetaVarC) {
      // perform a spherical wave fit first
      result = ScanSWF(rEvent, sFitResults);
      // exit function in case of a failure
      if (result != eSuccess) {
        return result;
      }

      return ScanSWFVarC(rEvent, sFitResults);
    } else {
      return ScanSWF(rEvent, sFitResults);
    }
  } else  // for a reuse fit
  {
    if (sFitResults.theta > fThetaVarC) {
      return SphericalWaveFitVarC(rEvent, sFitResults, 0);
    } else {
      return SphericalWaveFit(rEvent, sFitResults, 0);
    }
  }
}

VModule::ResultFlag RdWaveFit::CallConicalWaveFit(revt::REvent& rEvent, FitParameters& sFitResults,
                                                  bool reuseFit)
{

  VModule::ResultFlag result;

  if (fNumberOfStations < 4) {
    if (fInfoLevel > 0) {
      INFO(
          "Event contains less than 4 Stations: Only the PlaneWaveFit will be performed. Setting WaveFrontModel to ePlane");
      currentWaveFrontModel = ePlane;
      return CallPlaneWaveFit(rEvent, sFitResults, reuseFit);
    }
  }

  if (!reuseFit) {
    // perform a plane wave fit first
    /*OUT(eMinuit) << "Call plane fit before cone fit" << endl;
    result = CallPlaneWaveFit(rEvent, sFitResults, reuseFit);
    */

    // perform a spherical wave fit with variable speed of light first
    // to then use fixed shower axis from spherical wave fit for cone angle reconstruction
    OUT(eMinuit) << "Call spherical wave fit with variable speed of light" << endl;
    result = CallSphericalWaveFitVarC(rEvent, sFitResults, reuseFit);


    // exit function in case of a failure
    if (result != eSuccess) {
      return result;
    }
    OUT(eMinuit) << "call scan cone fit" << endl;

    return ScanCWF(rEvent, sFitResults);
  } else  // for a reuse fit
  {
    return ConicalWaveFit(rEvent, sFitResults, 0);
  }

}

const std::vector<double> RdWaveFit::GetTimeResiduals(const std::vector<Vector>& AntennaPositions,
                                                      const std::vector<double>& AntennaTimes,
                                                      const utl::Vector& ShowerAxis,
                                                      const unsigned int waveFrontModel,
                                                      const double gamma, const double rho) const
{
  double c = utl::kSpeedOfLight;         // in AugerUnits
  std::vector<double> ExpectedTimes;  // Expected time for each antenna if the test position is the source position
  std::vector<double> timeResiduals;

  double tau0(0);
  tau0 = TMath::Mean(AntennaTimes.begin(), AntennaTimes.end());
  /* an absolute event time for both the test and measured times. The absolute event time is choosen to be the mean of all measured/test times
   */

  // Calculate the expected times
  int n = AntennaPositions.size();
  double exp_time;

  //  calculate expected times for the different wave models
  for (int i = 0; i < n; ++i) {
    switch (waveFrontModel) {
      case ePlane:
        exp_time = (-(ShowerAxis * AntennaPositions[i])) / c;  // calculate (negative) projection
//           OUT(eObscure) << "GetTimeResiduals(): Plane Fit: Radius of ShowerAxis = " << ShowerAxis.GetMag() << std::endl;
        break;
      case eSpherical:
        exp_time = ((AntennaPositions[i] - ShowerAxis).GetMag()) / c;
//           OUT(eObscure) << "GetTimeResiduals(): SphericalWaveFit: Radius of ShowerAxis = " << ShowerAxis.GetMag() << std::endl;
        break;
      case eSphericalVarC:
        exp_time = ((AntennaPositions[i] - ShowerAxis).GetMag()) / (c * gamma);
//           OUT(eObscure) << "GetTimeResiduals(): SphericalWaveFitVarC: Radius of ShowerAxis = " << ShowerAxis.GetMag() << std::endl;
        break;
      case eConical:
        exp_time = (-ShowerAxis*AntennaPositions[i] + tan(rho) * fabs((fabs(ShowerAxis*AntennaPositions[i])*ShowerAxis-AntennaPositions[i]).GetMag())) / c;
//                       OUT(eObscure) << "GetTimeResiduals(): ConicalWaveFit: Radius of ShowerAxis = " << ShowerAxis.GetMag() << std::endl;
        break;
      default:
        WARNING(
            "GetTimeResiduals(): No method definded to calculate the expected time for that wave model");
        return timeResiduals;
    }
    ExpectedTimes.push_back(exp_time);
  }

  // get the minimum of ExpectedTimes
  double t0(0);
  t0 = TMath::Mean(ExpectedTimes.begin(), ExpectedTimes.end());

  // Calculate chi square
  for (int i = 0; i < n; ++i) {
//       OUT(eObscure) << "GetTimeResiduals(): TimeResidual of Station "<< i << " = " << (AntennaTimes[i]-tau0)-(ExpectedTimes[i]-t0) <<  std::endl;
    timeResiduals.push_back((AntennaTimes[i] - tau0) - (ExpectedTimes[i] - t0));
  }
  return timeResiduals;
}


  VModule::ResultFlag RdWaveFit::ComputeBaryCenter(const revt::REvent& rEvent,
      utl::Point& baryCenter)
  {

    const det::Detector& detector = det::Detector::GetInstance();
    const rdet::RDetector& rDetector = detector.GetRDetector();
    const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();

    // the absolute points in time and space
    const Point siteOrigin(0, 0, 0, siteCS);

    Vector barySum(0, 0, 0, siteCS);
    double weightSum = 0;

    for(REvent::ConstSignalStationIterator sIt = rEvent.SignalStationsBegin();
        sIt != rEvent.SignalStationsEnd(); ++sIt) {
      const rdet::Station& dStation = rDetector.GetStation(*sIt);
      const StationRecData& sRec = sIt->GetRecData();
      weightSum += sqrt(sRec.GetParameter(eSignal));
      barySum += (dStation.GetPosition() - siteOrigin) * sqrt(sRec.GetParameter(eSignal));
    }  //station loop

    barySum /= weightSum;
    baryCenter = siteOrigin + barySum;

    return eSuccess;
  }


}