HybridGeometryFinderWG/HybridGeometryFinder.cc

Go to the documentation of this file.

#include "HybridGeometryFinder.h"

#include <det/Detector.h>

#include <evt/Event.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerRecData.h>

#include <fdet/Channel.h>
#include <fdet/Eye.h>
#include <fdet/FDetector.h>
#include <fdet/FTimeFitModel.h>
#include <fdet/Pixel.h>
#include <fdet/Telescope.h>

#include <fevt/Eye.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeRecData.h>
#include <fevt/FdComponentSelector.h>
#include <fevt/FEvent.h>
#include <fevt/Pixel.h>
#include <fevt/PixelRecData.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeRecData.h>

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

#include <sdet/SDetector.h>

#include <sevt/SEvent.h>
#include <sevt/Station.h>

#include <utl/AxialVector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/CorrelationMatrix.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/NumericalErrorPropagation.h>
#include <utl/Point.h>
#include <utl/PhysicalConstants.h>
#include <utl/Reader.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/Transformation.h>
#include <utl/UTMPoint.h>
#include <utl/Vector.h>

#include <TMinuit.h>

#include <algorithm>
#include <cmath>
#include <iostream>
#include <map>
#include <set>
#include <vector>
#include <iostream>

using namespace std;

namespace HybridGeometryFinderWG {

  fevt::Eye* HybridGeometryFinder::fgCurEye = NULL;

  sevt::Station* HybridGeometryFinder::fgFitStation = NULL;
  vector<unsigned short int> HybridGeometryFinder::fgStationList;

  int HybridGeometryFinder::fgHybridFitMode = 0;
  unsigned int HybridGeometryFinder::fgNDofAxis = 0;

  int HybridGeometryFinder::fgVerbosity = -1;
  int HybridGeometryFinder::fgLaserShots = 0;

  double HybridGeometryFinder::fgVerticalCurvature = 0.0001111 / utl::m;
  double HybridGeometryFinder::fgSDPError = 0.35 * utl::degree;

  double HybridGeometryFinder::fgCoreHeight = 1400. * utl::m;

  fwk::VModule::ResultFlag
  HybridGeometryFinder::Init()
  {

    fwk::CentralConfig* cc = fwk::CentralConfig::GetInstance();

    utl::Branch topB = cc->GetTopBranch("HybridGeometryFinderWG");

    if (!topB) {
      ERROR("Could not find branch HybridGeometryFinderWG");
      return eFailure;
    }

    topB.GetChild("verbosityLevel").GetData(fgVerbosity);
    topB.GetChild("SDPError").GetData(fgSDPError);
    topB.GetChild("TimeFitModel").GetData(fTimeFitModel);
    topB.GetChild("AxisMinPixel").GetData(fAxisMinPixel);
    topB.GetChild("MaxSDPDist").GetData(fMaxSDPDist);
    topB.GetChild("VerticalCurvature").GetData(fgVerticalCurvature);
    topB.GetChild("MaxStationToSDPDist").GetData(fMaxStationToSDPDist);
    topB.GetChild("MaxTimeResidue").GetData(fMaxTimeResidue);
    topB.GetChild("LaserShots").GetData(fgLaserShots);
    topB.GetChild("FdOnly").GetData(fFdOnly);
    topB.GetChild("PassThrough").GetData(fPassThrough);

    // Set the right time fit model
    fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();
    switch (fTimeFitModel) {
      case 0:
        timeFitModel.SetModel(fdet::FTimeFitModel::eVacuumAtm);
        break;
      case 1:
        timeFitModel.SetModel(fdet::FTimeFitModel::eRealisticAtm);
        break;
      default:
        ERROR("Invalid input parameter for TimeFitModel");
        return eFailure;
    }

    // Configuration output
    if (fgVerbosity >= 0) {
      ostringstream info;
      info << " Version: "
        << GetVersionInfo(VModule::eRevisionNumber) << "\n";
      info << "                                Verbosity level: " << fgVerbosity << "\n";
      info << "                                      SDP Error: " << fgSDPError << "\n";
      info << "                                 Time Fit Model: " << fTimeFitModel << "\n";
      info << "          Minimum number of pixels for axis fit: " << fAxisMinPixel << "\n";
      info << "       Maximal difference in Chi_i for the axis: " << fMaxSDPDist << "\n";
      info << "                             Vertical curvature: " << fgVerticalCurvature << "\n";
      info << "                  Maximal Axis-Station Distance: " << fMaxStationToSDPDist << "\n";
      info << "                           Maximal Time Residue: " << fMaxTimeResidue << "\n";
      info << "                        Reconstruct laser shots: " << fgLaserShots << "\n";
      info << "                                  Allow Fd Only: " << fFdOnly << "\n";

      INFO(info);
    }

    return eSuccess;

  } // fwk::VModule::ResultFlag HybridGeometryFinder::Init()


  fwk::VModule::ResultFlag
  HybridGeometryFinder::Run(evt::Event& event)
  {

    INFO("Run HybridGeometryFinderWG");

    if (!event.HasFEvent()) {
      if (fgVerbosity >= 0) {
        INFO("This event has no FEvent --- Skip event!");
      }
      return eContinueLoop;
    }

    if (!event.HasSEvent() && !fFdOnly) {
      if (fgVerbosity >= 0) {
        INFO("This event has no SEvent --- Skip event!");
      }
      return eSuccess;
    }

    // Store the event for later use
    fEvent = &event;

    // Run the reconstruction for each eye
    for (fevt::FEvent::EyeIterator eyeIter = event.GetFEvent().EyesBegin(fevt::ComponentSelector::eHasData);
        eyeIter != event.GetFEvent().EyesEnd(fevt::ComponentSelector::eHasData);
        ++eyeIter) {

      if (fgVerbosity >= 0) {
        ostringstream info;
        info << "\n Start reconstructing eye " << eyeIter->GetId() << "\n";
        INFO(info);
      }

      // Store the eye for later use
      fgCurEye = &*eyeIter;
      const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(*fgCurEye);
      bool foundAxis = false;

      if (!fgCurEye->HasRecData())
        continue;
      // Skip this eye if there are not enough pixel
      if (fgCurEye->GetRecData().GetNPulsedPixels() < fAxisMinPixel) {
        if (fgVerbosity >= 0) {
          ostringstream info;
          info << "There are not enough pulsed pixels (";
          info << fgCurEye->GetRecData().GetNPulsedPixels();
          info << " < ";
          info << fAxisMinPixel;
          info << ") for an axis fit in eye ";
          info << fgCurEye->GetId();
          info << " --- jump to next eye!";
          INFO(info);
        }

        continue;
      }

      // Get the data from previous fits and take it as first guess
      const int prevEye = GetDataFromPreviousFit();
      if (prevEye == 0) {
        if (fgVerbosity >= 0) {
          INFO("Could not find results from previous fits --- jump to next eye!");
        }

        continue;
      }
      else if (!(SelectHybridStation(prevEye)
            || fFdOnly
            || detEye.IsVirtual())) {
        if (fgVerbosity >= 0) {
          INFO("No hybrid station found --- jump to next eye!");
        }

        continue;
      }

      // Recalculate the geometry for each telescope
      CalculateGeometryForTels();

      // Select the pixels for the time fit
      SelectPixels();

      // If there are not enough time fit pixels stop this eye
      if (fgCurEye->GetRecData().GetNTimeFitPixels() < fAxisMinPixel) {
        if (fgVerbosity >= 0) {
          ostringstream info;
          info << "There are not enough time fit pixels (";
          info << fgCurEye->GetRecData().GetNTimeFitPixels();
          info << " < ";
          info << fAxisMinPixel;
          info << ") for an axis fit in eye ";
          info << fgCurEye->GetId();
          info << " --- jump to next eye!";
          INFO(info);
        }

        continue;
      }

      if (fPassThrough) {
        foundAxis = true;
      }

      // First fit using the station core distance
      fgHybridFitMode = 0;
      if (FitAxis()) {
        foundAxis = true;
      }

      // Second fit using the SD/FD times
      fgHybridFitMode = 1;
      if (FitAxis()) {
        foundAxis = true;
      }

      // Recalculate the geometry for each telescope
      CalculateGeometryForTels();

      // Recalculate the time chi square
      RecalculateChiSquare();

      // Readmit pixels that might have been left out before but can now be used in the fit
      while (ReadmitPixel()) {
        fgHybridFitMode = 1;
        if (FitAxis()) {
          foundAxis = true;
        }

        CalculateGeometryForTels();

        RecalculateChiSquare();
      }

      // Remove pixels which have a too high time residual
      while (RemovePixel()) {
        fgHybridFitMode = 1;
        if (FitAxis()) {
          foundAxis = true;
        }

        CalculateGeometryForTels();

        RecalculateChiSquare();
      }

      // If no axis was found, don't store the data
      if (!foundAxis) {
        if (fgVerbosity >= 0) {
          ostringstream info;
          info << "No axis for eye " << eyeIter->GetId() << " found --- jump to next eye.";
          INFO(info);
        }

        continue;
      }

      // Store the axis data
      StoreData();

    } // for all eyes with data

    return eSuccess;

  } // fwk::VModule::ResultFlag HybridGeometryFinder::Run(evt::Event& event)


  fwk::VModule::ResultFlag
  HybridGeometryFinder::Finish()
  {

    return eSuccess;

  } // fwk::VModule::ResultFlag HybridGeometryFinder::Finish()


  int
  HybridGeometryFinder::GetDataFromPreviousFit()
  {

    const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
    fevt::FEvent& fevent = fEvent->GetFEvent();

    set<unsigned int> eyeIds;

    // if the eye is virtual, take the data from any physical eye
    if (detFD.GetEye(*fgCurEye).IsVirtual()) {
      for (fdet::FDetector::EyeIterator eyeIter = detFD.EyesBegin(fdet::FDetComponentSelector::ePhysical);
          eyeIter != detFD.EyesEnd(fdet::FDetComponentSelector::ePhysical);
          ++eyeIter) {

        if (fevent.HasEye(eyeIter->GetId())) {
          eyeIds.insert(eyeIter->GetId());
        }

      }
    }
    // if the eye is physical, take the data from this eye
    // or in case no previous data is found take it from any eye
    // UPDATE - treat physical eyes the same way as virtual eyes
    else {

      bool hasData = false;
      if (fgCurEye->HasRecData()) {
        if ((fgCurEye->GetRecData().HasFRecShower()) && (fgCurEye->GetRecData().GetRp() != 0.)) {
          hasData = true;
        }
      }

      if (hasData && fPassThrough) {
        eyeIds.insert(fgCurEye->GetId());
      }
      else {

        for (fdet::FDetector::EyeIterator eyeIter = detFD.EyesBegin(fdet::FDetComponentSelector::ePhysical);
            eyeIter != detFD.EyesEnd(fdet::FDetComponentSelector::ePhysical);
            ++eyeIter) {

          if (fevent.HasEye(eyeIter->GetId())) {
            eyeIds.insert(eyeIter->GetId());
          }

        }
      }
    }

    // Find the most time fit pixels
    int maxNFitPixels = 0;
    unsigned int maxNFitPixelsId = 0;
    for (set<unsigned int>::const_iterator IdIter = eyeIds.begin();
        IdIter != eyeIds.end();
        ++IdIter) {

      // Check if the previous modules found something
      if (!fevent.GetEye(*IdIter).HasRecData()) {
        continue;
      }
      else if (!fevent.GetEye(*IdIter).GetRecData().HasFRecShower()) {
        continue;
      }

      const fevt::EyeRecData& eyeRecData = fevent.GetEye(*IdIter).GetRecData();
      double nFitPixels = eyeRecData.GetNTimeFitPixels();

      if (nFitPixels > maxNFitPixels) {
        maxNFitPixels = nFitPixels;
        maxNFitPixelsId = *IdIter;
      }

      if (fgVerbosity >= 2) {
        cout << "Eye " << *IdIter << " has " << nFitPixels << " pixels in the time fit" << endl;
      }

    } // for eyeIds

    int takeFromEye = 0;
    if (maxNFitPixelsId != 0) {
      takeFromEye = maxNFitPixelsId;
    }
    else {
      return takeFromEye;
    }

    const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(takeFromEye);
    fevt::Eye& eventEye = fEvent->GetFEvent().GetEye(takeFromEye);

    const fevt::EyeRecData& eyeRecData = eventEye.GetRecData();

    // Get the data
    fSDP = eyeRecData.GetSDP();
    fChiZero = eyeRecData.GetChiZero();
    fRp = eyeRecData.GetRp();
    fTZero = eyeRecData.GetTZero();
    fAxis = eyeRecData.GetFRecShower().GetAxis();
    fCore = eyeRecData.GetFRecShower().GetCorePosition();

    // Readjust the TZero according to the trigger offset of the different eyes
    const utl::TimeInterval offset = eventEye.GetHeader().GetTimeStamp() - fgCurEye->GetHeader().GetTimeStamp();
    fTZero += offset.GetInterval();

    // Shift the core to a height of fgCoreHeight (1400m)
    const utl::ReferenceEllipsoid& refEllipsoid = utl::ReferenceEllipsoid::GetWGS84();
    const double coreHeight = refEllipsoid.PointToLatitudeLongitudeHeight(fCore).get<2>();
    utl::CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(fCore);
    fCore -= fAxis * (coreHeight - fgCoreHeight) / fAxis.GetZ(coreCS);
    coreCS = fwk::LocalCoordinateSystem::Create(fCore);
    double height;
    try {
      const utl::UTMPoint UTMCore(fCore, refEllipsoid);
      fNorthing = UTMCore.GetNorthing();
      fEasting = UTMCore.GetEasting();
      height = UTMCore.GetHeight();
    }
    catch (utl::UTMPoint::UTMZoneException& zone_e) {
      ERROR ("UTMZoneException: core invalid");
      return 0;
    }
    catch (utl::UTMPoint::UTMException& e) {
      ERROR ("UTMException: core invalid");
      return 0;
    }


    // Check the direction of the shower
    const bool downGoing = fAxis.GetTheta(coreCS) < utl::kPiOnTwo;

    // Find the point at which T0 is defined
    const utl::Point telPos = detEye.GetPosition();
    utl::Vector telToRp = utl::cross(fSDP, fAxis);
    telToRp.Normalize();
    telToRp *= fRp;
    const utl::Point pointRp = telPos + telToRp;

    // Shift T0 from the pointRp to the core
    double tCore = fTZero;
    if (pointRp.GetZ(coreCS) > fCore.GetZ(coreCS)) {
      if (downGoing) {
        tCore += (pointRp - fCore).GetMag() / utl::kSpeedOfLight;
      }
      else {
        tCore -= (pointRp - fCore).GetMag() / utl::kSpeedOfLight;
      }
    }
    else {
      if (downGoing) {
        tCore -= (pointRp - fCore).GetMag() / utl::kSpeedOfLight;
      }
      else {
        tCore += (pointRp - fCore).GetMag() / utl::kSpeedOfLight;
      }
    }

    // Store the core and axis for the hybrid fit
    fAxisTheta = fAxis.GetTheta(coreCS);
    fAxisPhi = fAxis.GetPhi(coreCS);
    fTCore = tCore;

    if (fgVerbosity >= 1) {
      cout << "The following parameters from the previous fits in eye "
        << takeFromEye
        << " are taken as first guess\n"
        << "           ChiZero = " << fChiZero / utl::degree << "deg,\n"
        << "             TZero = " << fTZero / utl::ns << "ns,\n"
        << "                Rp = " << fRp / utl::m << "m,\n"
        << "    Core(Northing) = " << fNorthing / utl::m << "m,\n"
        << "     Core(Easting) = " << fEasting / utl::m << "m,\n"
        << "    Core(Altitude) = " << height / utl::m << "m,\n"
        << " AxisTheta(coreCS) = " << fAxisTheta / utl::degree << "deg,\n"
        << "   AxisPhi(coreCS) = " << fAxisPhi / utl::degree << "deg,\n"
        << " SDPsTheta(coreCS) = " << fSDP.GetTheta(coreCS) / utl::degree << "deg,\n"
        << "   SDPsPhi(coreCS) = " << fSDP.GetPhi(coreCS) / utl::degree << "deg,\n"
        << "             TCore = " << fTCore / utl::ns << "ns,\n"
        << " The Shower is " << (downGoing ? "down-going" : "up-going") << endl;
    }

    return takeFromEye;

  } // int HybridGeometryFinder::GetDataFromPreviousFit()


  void
  HybridGeometryFinder::CalculateGeometryForTels()
  {

    const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(fgCurEye->GetId());
    for (fdet::Eye::TelescopeIterator telIter = detEye.TelescopesBegin();
        telIter != detEye.TelescopesEnd();
        ++telIter) {

      const int telId = telIter->GetId();
      fTelParameters[telId] = CalculateAxisForTel(fCore, fAxis, fTCore, telId);

    }

  } // void HybridGeometryFinder::CalculateGeometryForTels()


  vector<double>
  HybridGeometryFinder::CalculateAxisForTel(const utl::Point core,
                                            utl::Vector axis,
                                            const double tCore,
                                            const int telId)
  {

    const utl::Point telPos = det::Detector::GetInstance().GetFDetector().GetEye(*fgCurEye).GetTelescope(telId).GetPosition();

    // Calculate everything in the local coordinate system at the telescope position
    const utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);

    // Check the pointing of the axis
    const bool downGoing = axis.GetTheta(telCS) < utl::kPiOnTwo;

    // Shift the core along the axis such that it is at the same height as the telescope
    const utl::Point telCore = core - axis * core.GetZ(telCS) / axis.GetZ(telCS);

    // Get the SDP for this telescope
    const utl::Vector telToCore = telCore - telPos;
    utl::AxialVector sdp;
    sdp = utl::cross(axis, telToCore);
    sdp.Normalize();

    // Get the Chi0
    double chiZero = utl::Angle(axis, telToCore);

    // Get Rp
    const double rp = telToCore.GetMag() * sin(chiZero);

    // If the shower is up-going Chi0 < 0
    if (!downGoing) {
      chiZero = utl::kPi - chiZero;
      chiZero *= -1;
    }

    // Get the Point of Rp
    utl::Vector telToRp = utl::cross(sdp, axis);
    telToRp.Normalize();
    telToRp *= rp;
    const utl::Point pointRp = telPos + telToRp;

    // Shift T0 from the core to pointRp
    double tZero = tCore;
    if (pointRp.GetZ(telCS) > core.GetZ(telCS)) {
      if (downGoing) {
        tZero -= (pointRp - core).GetMag() / utl::kSpeedOfLight;
      }
      else {
        tZero += (pointRp - core).GetMag() / utl::kSpeedOfLight;
      }
    }
    else {
      if (downGoing) {
        tZero += (pointRp - core).GetMag() / utl::kSpeedOfLight;
      }
      else {
        tZero -= (pointRp - core).GetMag() / utl::kSpeedOfLight;
      }
    }

    vector<double> returnValues;

    returnValues.push_back(telCore.GetX(telCS));
    returnValues.push_back(telCore.GetY(telCS));
    returnValues.push_back(axis.GetTheta(telCS));
    returnValues.push_back(axis.GetPhi(telCS));
    returnValues.push_back(sdp.GetTheta(telCS));
    returnValues.push_back(sdp.GetPhi(telCS));
    returnValues.push_back(chiZero);
    returnValues.push_back(tZero);
    returnValues.push_back(rp);

    return returnValues;

  } // vector<double> HybridGeometryFinder::CalculateAxisForTel(const utl::Point core, utl::Vector axis, const double tCore, const int telId)


  void
  HybridGeometryFinder::SelectPixels()
  {

    const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();

    // if the eye is physical, leave the time fit pixels from the previous fit
    // UPDATE - treat physical eyes the same way as virtual eyes

    if (!detFD.GetEye(*fgCurEye).IsVirtual() && fPassThrough) {
      return;
    }


    // Empty TimeFitPixels list from previous fits
    fevt::EyeRecData& eyeRecData = fgCurEye->GetRecData();
    for (fevt::EyeRecData::PixelIterator pixelIter = eyeRecData.TimeFitPixelsBegin();
        pixelIter != eyeRecData.TimeFitPixelsEnd();
        ) {

      pixelIter = eyeRecData.RemoveTimeFitPixel(pixelIter);

    }

    // Check if the pulsed pixels should be used for the time fit
    const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();
    for (fevt::EyeRecData::PixelIterator pixelIter = eyeRecData.PulsedPixelsBegin();
        pixelIter != eyeRecData.PulsedPixelsEnd();
        ++pixelIter) {

      const int telId = pixelIter->GetTelescopeId();

      // Get the parameters for this telescope
      const vector<double> parameters = fTelParameters[telId];
      const double coreX = parameters[0];
      const double coreY = parameters[1];
      const double sdpTheta = parameters[4];
      const double sdpPhi = parameters[5];
      const double chiZero = parameters[6];
      const double tZero = parameters[7];
      const double rp = parameters[8];

      // The parameters for each telescope are calculated in telCS
      const utl::Point telPos = detFD.GetEye(*fgCurEye).GetTelescope(telId).GetPosition();
      const utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);

      const utl::Point core(coreX, coreY, 0, telCS);
      const utl::AxialVector sdp(1, sdpTheta, sdpPhi, telCS, utl::AxialVector::kSpherical);

      const utl::Vector vertical(0, 0, 1, telCS, utl::Vector::kCartesian);
      const utl::AxialVector horizontalInSDP = utl::cross(sdp, vertical);

      const utl::Vector pixelDirection = detFD.GetPixel(*pixelIter).GetDirection();
      const utl::Vector pixelDirInSDP = pixelDirection - (sdp * pixelDirection) * sdp;

      // Get the Chi_i for this pixel
      const double pixelChi = utl::Angle(pixelDirInSDP, horizontalInSDP);

      // Get the pixel time
      const double pixelTime = pixelIter->GetRecData().GetT_i().GetNanoSecond();
      const double pixelTimeError = pixelIter->GetRecData().GetT_iError().GetInterval();
      const double pixelTimeExpected = timeFitModel.GetTimeAtAperture(tZero, rp, chiZero, pixelChi, sdp.GetTheta(telCS), fgCurEye->GetId(), telId);
      const double tmp = pixelTime - pixelTimeExpected;
      const double timeResidue = tmp / pixelTimeError;

      const double angleWithSDP = abs(utl::kPiOnTwo - utl::Angle(sdp, pixelDirection));

      // Add the pixels if it is within the boundaries
      if ((abs(timeResidue) < 6 * fMaxTimeResidue) && (angleWithSDP < fMaxSDPDist)) {
        eyeRecData.AddTimeFitPixel(*pixelIter);
      }

      if (fgVerbosity >= 2) {
        cout << "Pixel " << pixelIter->GetId()
          << " in telescope " << pixelIter->GetTelescopeId()
          << " has time=" << pixelTime << ","
          << " residue=" << timeResidue << ","
          << " angleWithSDP=" << angleWithSDP / utl::degree
          << " and will";
        if ((abs(timeResidue) < 6 * fMaxTimeResidue) && (angleWithSDP < fMaxSDPDist)) {
          cout << " be added to the time fit" << endl;
        }
        else {
          cout << " NOT be added to the time fit" << endl;
        }
      }

    } // for pulsed pixels

    if (fgVerbosity >= 1) {
      ostringstream info;
      info << "There are " << eyeRecData.GetNTimeFitPixels()
        << " pixels selected for the time fit";
      INFO(info);
    }

  } // void HybridGeometryFinder::SelectPixels()


  bool
  HybridGeometryFinder::SelectHybridStation(const int prevEye)
  {

    // Clear list from last eye
    fgStationList.clear();
    fgFitStation = NULL;

    // Get the list of stations from the previous fit
    fevt::EyeRecData& eyeRecData = fEvent->GetFEvent().GetEye(prevEye).GetRecData();
    if (!eyeRecData.HasFRecShower()) {
      eyeRecData.MakeFRecShower();
    }
    evt::ShowerFRecData& frecShower = eyeRecData.GetFRecShower();

    fgStationList = frecShower.GetStationIds();

    if (fgStationList.size() == 0 || fgStationList.at(0) == 9999) {
      fgStationList.clear();
      return false;
    }

    sevt::SEvent& sevent = fEvent->GetSEvent();

    // Store the hybrid station for later use
    for (sevt::SEvent::StationIterator stationIter = sevent.StationsBegin();
        stationIter != sevent.StationsEnd();
        ++stationIter) {

      if (stationIter->GetId() == fgStationList.at(0)) {
        fgFitStation = &*stationIter;
      }

    }

    if (fgVerbosity >=1) {
      ostringstream info;
      info << "Use hybrid station "
        << fgFitStation->GetId()
        << " from previous fit";
      INFO(info);
    }

    return true;

  } // bool HybridGeometryFinder::SelectHybridStation(const int prevEye)


  bool
  HybridGeometryFinder::FitAxis()
  {

    // The number of parameters for the fit
    const int nParameters = 5;

    // The Minuit object
    TMinuit theMinuit(nParameters);

    // The verbosity level
    theMinuit.SetPrintLevel(-1);

    // The fit function
    theMinuit.SetFCN(HybridGeometryFinder::MinuitFitFuncAxis);

    // The arguments passed to theMinuit
    double argList[nParameters];
    int errorFlag = 0;

    // up: this is a chi^2 fit
    argList[0] = 1;
    theMinuit.mnexcm("SET ERR", argList, 1, errorFlag);

    if (errorFlag) {
      ERROR("Minuit SET ERR failed");
    }

    // The start parameters
    double startParameters[nParameters];

    startParameters[0] = fNorthing;
    startParameters[1] = fEasting;
    startParameters[2] = fAxisTheta;
    startParameters[3] = fAxisPhi;
    startParameters[4] = fTCore;

    // The step size
    double stepSize[nParameters];

    stepSize[0] = 10;
    stepSize[1] = 10;
    stepSize[2] = 0.01;
    stepSize[3] = 0.01;
    stepSize[4] = 10;

    // Give the start parameters and step sizes to theMinuit
    theMinuit.mnparm(0, "Northing", startParameters[0], stepSize[0], 0, 0, errorFlag);
    theMinuit.mnparm(1, "Easting", startParameters[1], stepSize[1], 0, 0, errorFlag);
    theMinuit.mnparm(2, "AxisTheta", startParameters[2], stepSize[2], 0, 0, errorFlag);
    theMinuit.mnparm(3, "AxisPhi", startParameters[3], stepSize[3], 0, 0, errorFlag);
    theMinuit.mnparm(4, "TCore", startParameters[4], stepSize[4], 0, 0, errorFlag);

    // The maximal number of calls for Minuit
    argList[0] = 1500;

    // The tolerance
    argList[1] = 0.01;

    // The actual fitting via MINIMIZE
    theMinuit.mnexcm("MINIMIZE", argList, 2, errorFlag);

    // In the case fit failed, don't use the fitted data
    if (errorFlag) {
      ERROR("Minuit MINIMIZE failed");
      return false;
    }

    theMinuit.mnexcm("HESSE", argList, 2, errorFlag);

    // Check the status of the covariance matrix
    // if bad use MINOS
    double fmin, fedm, errdef; // unused dummies
    int npari, nparx; //unused dummies
    int istat = 0;
    theMinuit.mnstat(fmin, fedm, errdef, npari, nparx, istat);
    if (istat < 3) { // no accurate covariance matrix
      theMinuit.mnexcm("MINOS", argList, 2, errorFlag);
      theMinuit.mnstat(fmin, fedm, errdef, npari, nparx, istat);
      if (istat < 3)
        WARNING("Neither HESSE nor MINOS found accurate covariance matrix");
    }

    // Get the fitted parameters back from Minuit
    theMinuit.GetParameter(0, fNorthing, fNorthingError);
    theMinuit.GetParameter(1, fEasting, fEastingError);
    theMinuit.GetParameter(2, fAxisTheta, fAxisThetaError);
    theMinuit.GetParameter(3, fAxisPhi, fAxisPhiError);
    theMinuit.GetParameter(4, fTCore, fTCoreError);
    theMinuit.mnemat(&fCovMatrixAxis[0][0], nParameters);
    fChiSqAxis = theMinuit.fAmin;

    // Calculate fAxis and fCore
    const utl::CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
    const utl::ReferenceEllipsoid& refEllipsoid = utl::ReferenceEllipsoid::GetWGS84();

    // The core is defined at fgCoreHeigth (1400m)
    const utl::UTMPoint UTMCore(fNorthing, fEasting, fgCoreHeight, 19, 'H', refEllipsoid);

    // Create a coordinate system at the core point
    fCore = UTMCore.GetPoint(siteCS);
    const utl::CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(fCore);

    // The axis is defined in the coreCS
    fAxis = utl::Vector(1, fAxisTheta, fAxisPhi, coreCS, utl::Vector::kSpherical);

    if (fgVerbosity >= 1) {
      ostringstream info;
      info << "The parameters of the fitted axis for eye " << fgCurEye->GetId() << " are:" << endl
        << "        Northing = " << fNorthing / utl::m
        << " +/- " << fNorthingError / utl::m << "m," << endl
        << "         Easting = " << fEasting / utl::m
        << " +/- " << fEastingError / utl::m << "m," << endl
        << "       AxisTheta = " << fAxisTheta / utl::degree
        << " +/- " << fAxisThetaError / utl::degree << "deg," << endl
        << "         AxisPhi = " << fAxisPhi / utl::degree
        << " +/- " << fAxisPhiError / utl::degree << "deg," << endl
        << "           TCore = " << fTCore / utl::ns
        << " +/- " << fTCoreError / utl::ns << "ns," << endl
        << " ChiSquare / NDF = " << fChiSqAxis << " / " << fgNDofAxis
        << " = " << fChiSqAxis / fgNDofAxis;
      INFO(info);
    }

    return true;

  } // bool HybridGeometryFinder::FitAxis()


  void
  HybridGeometryFinder::MinuitFitFuncAxis(int& /*npar*/,
                                          double* /*gin*/,
                                          double& f,
                                          double* par,
                                          int /*iflag*/)
  {

    double northing = par[0];
    double easting = par[1];
    double axisTheta = par[2];
    double axisPhi = par[3];
    double tCore = par[4];

    // Create the coordinate system at the core
    const utl::CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
    const utl::ReferenceEllipsoid& refEllipsoid = utl::ReferenceEllipsoid::GetWGS84();

    // The core is defined at fgCoreHeight (1400m)
    const utl::UTMPoint UTMCore(northing, easting, fgCoreHeight, 19, 'H', refEllipsoid);

    // Create a coordinate system at the core point
    utl::Point core = UTMCore.GetPoint(siteCS);
    const utl::CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);

    // The axis is defined in the coreCS
    const utl::Vector axis(1, axisTheta, axisPhi, coreCS, utl::Vector::kSpherical);

    // Calculate the axis parameters for each telescope
    map<int, vector<double> > telParameter;
    for (fevt::Eye::TelescopeIterator telIter = fgCurEye->TelescopesBegin(fevt::ComponentSelector::eHasData);
        telIter != fgCurEye->TelescopesEnd(fevt::ComponentSelector::eHasData);
        ++telIter) {

      const int telId = telIter->GetId();

      telParameter[telId] = CalculateAxisForTel(core, axis, tCore, telId);

    }

    const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
    const fevt::EyeRecData& eyeRecData = fgCurEye->GetRecData();
    const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();

    const int eyeId = fgCurEye->GetId();

    // Loop over all pixels and calculate their
    // contribution to the time and SDP fit
    double chiSquareTime = 0.;
    double chiSquareSDP = 0.;
    double chiSquareSD = 0.;
    double totalCharge = 0.;

    for (fevt::EyeRecData::ConstPixelIterator pixelIter = eyeRecData.TimeFitPixelsBegin();
        pixelIter != eyeRecData.TimeFitPixelsEnd();
        ++pixelIter) {

      const fevt::PixelRecData& pixelRecData = pixelIter->GetRecData();

      // Get the axis parameters for this telescope
      const int telId = pixelIter->GetTelescopeId();
      const vector<double> parameters = telParameter[telId];
      const double coreX = parameters[0];
      const double coreY = parameters[1];
      const double sdpTheta = parameters[4];
      const double sdpPhi = parameters[5];
      const double chiZero = parameters[6];
      const double tZero = parameters[7];
      const double rp = parameters[8];

      // The parameters are defined in telCS
      const utl::Point telPos = detFD.GetEye(eyeId).GetTelescope(telId).GetPosition();
      const utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);

      // Get the pulse parameters
      const double pixelTime = pixelRecData.GetT_i().GetNanoSecond();
      const double pixelTimeError = pixelRecData.GetT_iError().GetInterval();
      const double pixelCharge = pixelRecData.GetTotalCharge();

      // Recalculate Chi_i for the SDP corresponding to this axis
      const utl::Point core(coreX, coreY, 0, telCS);
      const utl::AxialVector sdp(1, sdpTheta, sdpPhi, telCS, utl::AxialVector::kSpherical);
      const utl::Vector vertical(0, 0, 1, telCS, utl::Vector::kCartesian);
      const utl::AxialVector horizontalInSDP = utl::cross(sdp, vertical);
      const utl::Vector pixelDirection = detFD.GetPixel(*pixelIter).GetDirection();
      const utl::Vector pixelDirInSDP = pixelDirection - (sdp * pixelDirection) * sdp;
      const double pixelChi = utl::Angle(pixelDirInSDP, horizontalInSDP);

      // ----------------------------------------
      // Begin the time fit part

      // Calculate the expected pixel time
      const double pixelTimeExpected = timeFitModel.GetTimeAtAperture(tZero, rp, chiZero, pixelChi, sdp.GetTheta(telCS), eyeId, telId);

      const double tmp = pixelTime - pixelTimeExpected;

      chiSquareTime += tmp * tmp / (pixelTimeError * pixelTimeError);

      // End the time fit part
      // ----------------------------------------
      // Begin the SDP fit part

      const double angDistToSDP = utl::kPiOnTwo - utl::Angle(pixelDirection, sdp);

      chiSquareSDP += angDistToSDP * angDistToSDP * pixelCharge / (fgSDPError * fgSDPError);

      totalCharge += pixelCharge;

      // End the SDP fit part
      // ----------------------------------------

      if (fgVerbosity >=2) {
        cout << "Pixel no. " << pixelIter->GetId()
          << " in Tel " << pixelIter->GetTelescopeId()
          << " ChiSquareTime = " << tmp*tmp/(pixelTimeError*pixelTimeError) << ","
          << " ChiSquareSDP = " << angDistToSDP*angDistToSDP/(fgSDPError*fgSDPError) << ","
          << " charge = " << pixelCharge
          << endl;
      }

    } // for time fit pixels

    // Normalize the SDP part
    chiSquareSDP /= totalCharge;
    chiSquareSDP *= eyeRecData.GetNTimeFitPixels();

    // Skip SD part if no station was found
    if (fgStationList.size() > 0 && fgFitStation != NULL && fgFitStation->HasRecData()) {

      // ------------------------------------------
      // Begin the SD fit part

      const sdet::SDetector& detSD = det::Detector::GetInstance().GetSDetector();
      const sdet::Station& detStation = detSD.GetStation(*fgFitStation);
      const utl::Point stationPos = detStation.GetPosition();

      // Estimate the distance of the station to the axis
      const utl::Vector coreToStation = stationPos - core;
      const double axisCoreStationAngle = utl::Angle(axis, coreToStation);
      const double deltaStationAxis = coreToStation.GetMag() * sin(axisCoreStationAngle);

      // Estimate the distance from the core to the closest point
      // of axis and station.
      // deltaCoreRp is > 0 if Rp is above the core
      // and < 0 if it is below
      const double deltaCoreRp = coreToStation.GetMag() * cos(axisCoreStationAngle);

      // Minimize the station to axis distance
      if (fgHybridFitMode == 0) {
        chiSquareSD = coreToStation.GetMag() * coreToStation.GetMag() / (10 * 10);
      }
      // Minimize the SD/FD time offset
      else if (fgHybridFitMode == 1) {
        //-----------------------------------------
        // First get the trigger time of the station as expected
        // by the FD trigger time and corrected by the axis.

        // Get the time at the Rp point for the station
        const double timeShift = deltaCoreRp / utl::kSpeedOfLight;
        // This is right regardless of up or down going showers
        double expectedStationTime = tCore - timeShift;

        // Get the time delay due to the curvature of the shower front
        // The time delay for a spherical shower front, the time delay is
        // (radius - sqrt(radius^2 - stationAxis^2)) / c
        // ~ radius * (1 - 1 + stationAxis^2 / (2 * radius^2))
        // = stationAxis^2 / (2 * radius)
        // with curvature = 1 / radius
        const double curvature = fgVerticalCurvature * axis.GetCosTheta(coreCS);
        const double timeDelay = curvature * deltaStationAxis * deltaStationAxis / (2 * utl::kSpeedOfLight);
        expectedStationTime += timeDelay;

        // The expected trigger time is relative to the beginning of the
        // FD FADC trace. Thus adding the eye trigger time (which is the
        // end of the trace) and substract the trace length
        double triggerTimeFD = fgCurEye->GetHeader().GetTimeStamp().GetGPSNanoSecond();
        const fdet::Channel& detChannel = detFD.GetEye(*fgCurEye).GetTelescope(1).GetChannel(1);
        const int numberFADCBins = detChannel.GetFADCTraceLength();
        const double lengthFADCBin = detChannel.GetFADCBinSize();
        triggerTimeFD -= numberFADCBins * lengthFADCBin;

        expectedStationTime += triggerTimeFD;

        // If laser shots are reconstructed a correction of 180ns
        // is needed to take into account the propagation from the
        // CLF to Celeste
        if (fgLaserShots) {
          expectedStationTime += 180 * utl::ns;
        }

        // Now get the real trigger time of the station

        // Get the trigger time for the station
        // Check if the GPS second is the same for FD and SD
        double timeSD = fgFitStation->GetRecData().GetSignalStartTime().GetGPSNanoSecond();
        const int secondFD = fgCurEye->GetHeader().GetTimeStamp().GetGPSSecond();
        const int secondSD = fgFitStation->GetRecData().GetSignalStartTime().GetGPSSecond();
        if (secondSD > secondFD) {
          timeSD += 1 * utl::s;
        }
        if (secondSD < secondFD) {
          timeSD -= 1 * utl::s;
        }

        // The difference between expected and real trigger time is fitted
        const double timeDiff = timeSD - expectedStationTime;

        const double timeErr = 100. * utl::ns;

        chiSquareSD = timeDiff * timeDiff / (timeErr * timeErr);
      }

      // End the SD fit part
      // ------------------------------

    }

    // Calculate the total chi^2
    double chiSquare = chiSquareTime + chiSquareSDP + chiSquareSD;

    f = chiSquare;
    fgNDofAxis = 2 * eyeRecData.GetNTimeFitPixels() - 5;

    if (fgVerbosity >= 2) {
      cout << "Minuit Run Axis FD:"
        << " Northing = " << northing / utl::m << ","
        << " Easting = " << easting / utl::m << ","
        << " AxisTheta = " << axisTheta / utl::degree << ","
        << " AxisPhi = " << axisPhi / utl::degree << ","
        << " TCore = " << tCore / utl::ns << ","
        << " ChiSquareTime = " << chiSquareTime << ","
        << " ChiSquareSDP = " << chiSquareSDP << ","
        << " ChiSquareSD = " << chiSquareSD << ","
        << " ChiSquare = " << chiSquare << ","
        << " NDF = " << fgNDofAxis
        << endl;
    }

  } // void HybridGeometryFinder::MinuitFitFuncAxis(int& npar, double* gin, double& f, double* par, int iflag)


  void
  HybridGeometryFinder::RecalculateChiSquare()
  {

    // Create the coordinate system at the core
    const utl::CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
    const utl::ReferenceEllipsoid& refEllipsoid = utl::ReferenceEllipsoid::GetWGS84();

    // The core is defined at fgCoreHeight (1400m)
    const utl::UTMPoint UTMCore(fNorthing, fEasting, fgCoreHeight, 19, 'H', refEllipsoid);

    // Create a coordinate system at the core point
    utl::Point core = UTMCore.GetPoint(siteCS);
    const utl::CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);

    // The axis is defined in the coreCS
    const utl::Vector axis(1, fAxisTheta, fAxisPhi, coreCS, utl::Vector::kSpherical);

    const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
    const fevt::EyeRecData& eyeRecData = fgCurEye->GetRecData();
    const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();

    const int eyeId = fgCurEye->GetId();

    //-------------------------------------------
    // Loop over all pixels and calculate their
    // contribution to the time and SDP fit

    double chiSquareTime = 0.;
    double chiSquareSDP = 0.;
    double chiSquareSD = 0.;
    double totalCharge = 0.;

    map<int, double> totalChargePerTel;
    map<int, int> pixelPerTel;

    // Initialise all telescope wise quantities
    for (fdet::Eye::TelescopeIterator telIter = detFD.GetEye(*fgCurEye).TelescopesBegin();
        telIter != detFD.GetEye(*fgCurEye).TelescopesEnd();
        ++telIter) {

      int telId = telIter->GetId();

      totalChargePerTel[telId] = 0.;
      fChiSqAxisPerTel[telId] = 0.;
      fNDofAxisPerTel[telId] = 0;
      fChiSqTimePerTel[telId] = 0.;
      fNDofTimePerTel[telId] = 0;
      fChiSqSDPPerTel[telId] = 0.;
      fNDofSDPPerTel[telId] = 0;

    }

    for (fevt::EyeRecData::ConstPixelIterator pixelIter = eyeRecData.TimeFitPixelsBegin();
        pixelIter != eyeRecData.TimeFitPixelsEnd();
        ++pixelIter) {

      const fevt::PixelRecData& pixelRecData = pixelIter->GetRecData();

      // Get the axis parameters for this telescope
      const int telId = pixelIter->GetTelescopeId();
      const vector<double> parameters = fTelParameters[telId];
      const double coreX = parameters[0];
      const double coreY = parameters[1];
      const double sdpTheta = parameters[4];
      const double sdpPhi = parameters[5];
      const double chiZero = parameters[6];
      const double tZero = parameters[7];
      const double rp = parameters[8];

      // The parameters are defined in telCS
      const utl::Point telPos = detFD.GetEye(eyeId).GetTelescope(telId).GetPosition();
      const utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);

      // Get the pulse parameters
      const double pixelTime = pixelRecData.GetT_i().GetNanoSecond();
      const double pixelTimeError = pixelRecData.GetT_iError().GetInterval();
      const double pixelCharge = pixelRecData.GetTotalCharge();

      // Recalculate Chi_i for the SDP corresponding to this axis
      const utl::Point core(coreX, coreY, 0, telCS);
      const utl::AxialVector sdp(1, sdpTheta, sdpPhi, telCS, utl::AxialVector::kSpherical);
      const utl::Vector vertical(0, 0, 1, telCS, utl::Vector::kCartesian);
      const utl::AxialVector horizontalInSDP = utl::cross(sdp, vertical);
      const utl::Vector pixelDirection = detFD.GetPixel(*pixelIter).GetDirection();
      const utl::Vector pixelDirInSDP = pixelDirection - (sdp * pixelDirection) * sdp;
      const double pixelChi = utl::Angle(pixelDirInSDP, horizontalInSDP);

      // ----------------------------------------
      // Begin the time fit part

      // Calculate the expected pixel time
      const double pixelTimeExpected = timeFitModel.GetTimeAtAperture(tZero, rp, chiZero, pixelChi, sdp.GetTheta(telCS), eyeId, telId);

      const double tmp = pixelTime - pixelTimeExpected;

      const double chiSquareTimePart = tmp * tmp / (pixelTimeError * pixelTimeError);
      chiSquareTime += chiSquareTimePart;
      fChiSqTimePerTel[telId] += chiSquareTimePart;

      // End the time fit part
      // ----------------------------------------
      // Begin the SDP fit part

      const double angDistToSDP = utl::kPiOnTwo - utl::Angle(pixelDirection, sdp);

      const double chiSquareSDPPart = angDistToSDP * angDistToSDP * pixelCharge / (fgSDPError * fgSDPError);
      chiSquareSDP += chiSquareSDPPart;
      fChiSqSDPPerTel[telId] += chiSquareSDPPart;

      totalCharge += pixelCharge;
      totalChargePerTel[telId] += pixelCharge;
      ++(fNDofTimePerTel[telId]);
      ++(fNDofSDPPerTel[telId]);
      fNDofAxisPerTel[telId] += 2;

      // End the SDP fit part
      // ----------------------------------------

    } // for time fit pixels

    // Normalize the SDP part
    if (totalCharge != 0)
      chiSquareSDP /= totalCharge;
    chiSquareSDP *= eyeRecData.GetNTimeFitPixels();

    // Normalize the SDP part for each telescope
    // NOTE: Due to the weighting the sum of the different
    // telescopes will not be the same as the chiSq for the eye
    for (map<int, double>::iterator iter = fChiSqSDPPerTel.begin();
        iter != fChiSqSDPPerTel.end();
        ++iter) {

      if (totalChargePerTel[iter->first] != 0)
        iter->second /= totalChargePerTel[iter->first];
      iter->second *= fNDofTimePerTel[iter->first];

    }

    double deltaStationAxis = 0.;
    // Skip SD part for a FD only reconstruction
    if (fgStationList.size() > 0 && fgFitStation != NULL && fgFitStation->HasRecData()) {

      // ------------------------------------------
      // Begin the SD fit part

      const sdet::SDetector& detSD = det::Detector::GetInstance().GetSDetector();
      const sdet::Station& detStation = detSD.GetStation(*fgFitStation);
      const utl::Point stationPos = detStation.GetPosition();

      // Estimate the distance of the station to the axis
      const utl::Vector coreToStation = stationPos - core;
      const double axisCoreStationAngle = utl::Angle(axis, coreToStation);
      deltaStationAxis = coreToStation.GetMag() * sin(axisCoreStationAngle);

      // Estimate the distance from the core to the closest point
      // of axis and station.
      // deltaCoreRp is > 0 if Rp is above the core
      // and < 0 if it is below
      const double deltaCoreRp = coreToStation.GetMag() * cos(axisCoreStationAngle);

      //-----------------------------------------
      // First get the trigger time of the station as expected
      // by the FD trigger time and corrected by the axis.

      // Get the time at the Rp point for the station
      const double timeShift = deltaCoreRp / utl::kSpeedOfLight;
      // This is right regardless of up or down going
      double expectedStationTime = fTCore - timeShift;

      // Get the time delay due to the curvature of the shower front
      // The time delay for a spherical shower front, the time delay is
      // (radius - sqrt(radius^2 - stationAxis^2)) / c
      // ~ radius * (1 - 1 + stationAxis^2 / (2 * radius^2))
      // = stationAxis^2 / (2 * radius)
      // with curvature = 1 / radius
      const double curvature = fgVerticalCurvature * axis.GetCosTheta(coreCS);
      const double timeDelay = curvature * deltaStationAxis * deltaStationAxis / (2 * utl::kSpeedOfLight);
      expectedStationTime += timeDelay;

      // The expected trigger time is relative to the beginning of the
      // FD FADC trace. Thus adding the eye trigger time (which is the
      // end of the trace) and substract the trace length
      double triggerTimeFD = fgCurEye->GetHeader().GetTimeStamp().GetGPSNanoSecond();
      const fdet::Channel& detChannel = detFD.GetEye(*fgCurEye).GetTelescope(1).GetChannel(1);
      const int numberFADCBins = detChannel.GetFADCTraceLength();
      const double lengthFADCBin = detChannel.GetFADCBinSize();
      triggerTimeFD -= numberFADCBins * lengthFADCBin;

      expectedStationTime += triggerTimeFD;

      // If laser shots are reconstructed a correction of 180ns
      // is needed to take into account the propagation from the
      // CLF to Celeste
      if (fgLaserShots) {
        expectedStationTime += 180 * utl::ns;
      }

      //-----------------------------------------
      // Now get the real trigger time of the station

      // Get the trigger time for the station
      // Check if the GPS second is the same for FD and SD
      double timeSD = fgFitStation->GetRecData().GetSignalStartTime().GetGPSNanoSecond();
      const int secondFD = fgCurEye->GetHeader().GetTimeStamp().GetGPSSecond();
      const int secondSD = fgFitStation->GetRecData().GetSignalStartTime().GetGPSSecond();
      if (secondSD > secondFD) {
        timeSD += 1 * utl::s;
      }
      if (secondSD < secondFD) {
        timeSD -= 1 * utl::s;
      }

      //-----------------------------------------
      // The difference between expected and real trigger time is fitted
      const double timeDiff = timeSD - expectedStationTime;

      const double timeErr = 100. * utl::ns;

      chiSquareSD = timeDiff * timeDiff / (timeErr * timeErr);
      fFitStationOffset = timeDiff;

      // End the SD fit part
      // ------------------------------

    }

    // Calculate the complete chisq and the ndfs per telescope
    for (map<int, double>::iterator iter = fChiSqSDPPerTel.begin();
        iter != fChiSqSDPPerTel.end();
        ++iter) {

      fChiSqAxisPerTel[iter->first] = fChiSqTimePerTel[iter->first] + fChiSqSDPPerTel[iter->first] + chiSquareSD;

      if (fNDofAxisPerTel[iter->first] >= 5)
        fNDofAxisPerTel[iter->first] -= 5;
      else
        fNDofAxisPerTel[iter->first] = 0;

      if (fNDofTimePerTel[iter->first] >= 3)
        fNDofTimePerTel[iter->first] -= 3;
      else
        fNDofTimePerTel[iter->first] = 0;

      if (fNDofSDPPerTel[iter->first] >= 2)
        fNDofSDPPerTel[iter->first] -= 2;
      else
        fNDofSDPPerTel[iter->first] = 0;

    }

    fChiSqTime = chiSquareTime;
    fNDofTime = eyeRecData.GetNTimeFitPixels() - 3;
    fChiSqSDP = chiSquareSDP;
    fNDofSDP = eyeRecData.GetNTimeFitPixels() - 2;
    fDeltaStationAxis = deltaStationAxis;

  } // void HybridGeometryFinder::RecalculateChiSquare()


  bool
  HybridGeometryFinder::ReadmitPixel()
  {

    //-------------------------------------------
    // Test for each pulsed pixel:
    // 1. Is it already in the fit -> continue
    // 2. Should it be readmitted
    const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
    fevt::EyeRecData& eyeRecData = fgCurEye->GetRecData();
    const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();

    bool readmitted = false;
    unsigned int countReadmitted = 0;

    for (fevt::EyeRecData::PixelIterator pixelIter = eyeRecData.PulsedPixelsBegin();
        pixelIter != eyeRecData.PulsedPixelsEnd();
        ++pixelIter) {

      bool inTimeFit = false;
      for (fevt::EyeRecData::ConstPixelIterator pixelIter2 = eyeRecData.TimeFitPixelsBegin();
          pixelIter2 != eyeRecData.TimeFitPixelsEnd();
          ++pixelIter2) {

        if ((pixelIter2->GetId() == pixelIter->GetId()) && (pixelIter2->GetTelescopeId() == pixelIter->GetTelescopeId())) {
          inTimeFit = true;
          break;
        }

      }

      if (inTimeFit) {
        continue;
      }

      const int telId = pixelIter->GetTelescopeId();

      const vector<double> parameters = fTelParameters[telId];
      const double coreX = parameters[0];
      const double coreY = parameters[1];
      const double sdpTheta = parameters[4];
      const double sdpPhi = parameters[5];
      const double chiZero = parameters[6];
      const double tZero = parameters[7];
      const double rp = parameters[8];

      // The parameters for each telescope are calculated in telCS
      const utl::Point telPos = detFD.GetEye(*fgCurEye).GetTelescope(telId).GetPosition();
      const utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);

      const utl::Point core(coreX, coreY, 0, telCS);
      const utl::AxialVector sdp(1, sdpTheta, sdpPhi, telCS, utl::AxialVector::kSpherical);

      const utl::Vector vertical(0, 0, 1, telCS, utl::Vector::kCartesian);
      const utl::AxialVector horizontalInSDP = utl::cross(sdp, vertical);

      const utl::Vector pixelDirection = detFD.GetPixel(*pixelIter).GetDirection();
      const utl::Vector pixelDirInSDP = pixelDirection - (sdp * pixelDirection) * sdp;

      // Get the Chi_i for this pixel
      const double pixelChi = utl::Angle(pixelDirInSDP, horizontalInSDP);

      // Get the pixel time
      const double pixelTime = pixelIter->GetRecData().GetT_i().GetNanoSecond();
      const double pixelTimeError = pixelIter->GetRecData().GetT_iError().GetInterval();
      const double pixelTimeExpected = timeFitModel.GetTimeAtAperture(tZero, rp, chiZero, pixelChi, sdp.GetTheta(telCS), fgCurEye->GetId(), telId);
      const double tmp = pixelTime - pixelTimeExpected;
      const double timeResidue = tmp / pixelTimeError;

      const double angleWithSDP = abs(utl::kPiOnTwo - utl::Angle(sdp, pixelDirection));

      // Add the pixel if it is within the boundaries
      if ((abs(timeResidue) < 3 * fMaxTimeResidue) && (angleWithSDP < fMaxSDPDist)) {
        eyeRecData.AddTimeFitPixel(*pixelIter);
        readmitted = true;
        ++countReadmitted;
      }

      if (fgVerbosity >= 2) {
        cout << "Pixel " << pixelIter->GetId()
          << " in telescope " << pixelIter->GetTelescopeId()
          << " has time=" << pixelTime << ","
          << " residue=" << timeResidue << ","
          << " angleWithSDP=" << angleWithSDP / utl::degree
          << " and will";
        if ((abs(timeResidue) < 3 * fMaxTimeResidue) && (angleWithSDP < fMaxSDPDist)) {
          cout << " be readmitted." << endl;
        }
        else {
          cout << " NOT be readmitted." << endl;
        }
      }

    }

    if (readmitted && fgVerbosity >= 1) {
      ostringstream info;
      info << "There are now " << eyeRecData.GetNTimeFitPixels()
        << " pixels for the axis fit"
        << " --- " << countReadmitted << " are readmitted.";
      INFO(info);
    }

    return readmitted;

  } // bool HybridGeometryFinder::ReadmitPixel()


  bool
  HybridGeometryFinder::RemovePixel()
  {

    bool removed = false;

    const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
    fevt::EyeRecData& eyeRecData = fgCurEye->GetRecData();
    const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();

    const int eyeId = fgCurEye->GetId();

    // Loop over all pixels and calculate their
    // contribution to the time fit

    double maxTimeResidue = 0.;
    int maxTimeResiduePixel = -1;
    int maxTimeResidueTel = -1;

    for (fevt::EyeRecData::ConstPixelIterator pixelIter = eyeRecData.TimeFitPixelsBegin();
        pixelIter != eyeRecData.TimeFitPixelsEnd();
        ++pixelIter) {

      const fevt::PixelRecData& pixelRecData = pixelIter->GetRecData();

      // Get the axis parameters for this telescope
      const int telId = pixelIter->GetTelescopeId();
      const vector<double> parameters = fTelParameters[telId];
      const double coreX = parameters[0];
      const double coreY = parameters[1];
      const double sdpTheta = parameters[4];
      const double sdpPhi = parameters[5];
      const double chiZero = parameters[6];
      const double tZero = parameters[7];
      const double rp = parameters[8];

      // The parameters are defined in telCS
      const utl::Point telPos = detFD.GetEye(eyeId).GetTelescope(telId).GetPosition();
      const utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);

      // Get the pulse parameters
      const double pixelTime = pixelRecData.GetT_i().GetNanoSecond();
      const double pixelTimeError = pixelRecData.GetT_iError().GetInterval();

      // Recalculate Chi_i for the SDP corresponding to this axis
      const utl::Point core(coreX, coreY, 0, telCS);
      const utl::AxialVector sdp(1, sdpTheta, sdpPhi, telCS, utl::AxialVector::kSpherical);
      const utl::Vector vertical(0, 0, 1, telCS, utl::Vector::kCartesian);
      const utl::AxialVector horizontalInSDP = utl::cross(sdp, vertical);
      const utl::Vector pixelDirection = detFD.GetPixel(*pixelIter).GetDirection();
      const utl::Vector pixelDirInSDP = pixelDirection - (sdp * pixelDirection) * sdp;
      const double pixelChi = utl::Angle(pixelDirInSDP, horizontalInSDP);

      // Calculate the expected pixel time
      const double pixelTimeExpected = timeFitModel.GetTimeAtAperture(tZero, rp, chiZero, pixelChi, sdp.GetTheta(telCS), eyeId, telId);

      const double tmp = pixelTime - pixelTimeExpected;
      const double timeResidue = abs(tmp / pixelTimeError);

      // Search for the maximum
      if (timeResidue > maxTimeResidue) {
        maxTimeResidue = timeResidue;
        maxTimeResiduePixel = pixelIter->GetId();
        maxTimeResidueTel = telId;
      }

      if (fgVerbosity >= 2) {
        cout << "Pixel " << pixelIter->GetId()
          << " in Telescope " << pixelIter->GetTelescopeId()
          << " has Residue = " << timeResidue << endl;
      }

    } // for pixels

    if ((maxTimeResidue > fMaxTimeResidue) && (eyeRecData.GetNTimeFitPixels() > fAxisMinPixel)) {
      for (fevt::EyeRecData::PixelIterator pixelIter = eyeRecData.TimeFitPixelsBegin();
          pixelIter != eyeRecData.TimeFitPixelsEnd(); ) {

        if ((int(pixelIter->GetId()) == maxTimeResiduePixel)
            && (int(pixelIter->GetTelescopeId()) == maxTimeResidueTel)) {
          pixelIter = eyeRecData.RemoveTimeFitPixel(pixelIter);
          removed = true;
          break;
        }
        else {
          ++pixelIter;
        }
      } // for pixel
    } // endif

    if (removed && fgVerbosity >= 1) {
      ostringstream info;
      info << "Pixel " << maxTimeResiduePixel
        << " in Telescope " << maxTimeResidueTel
        << " is removed from the time fit due to its time residue of " << maxTimeResidue;
      INFO(info);
    }

    return removed;

  } // bool HybridGeometryFinder::RemovePixel()


  void
  HybridGeometryFinder::StoreData()
  {

    // Store the SDP pixels used for the combined fit
    fevt::EyeRecData& eyeRecData = fgCurEye->GetRecData();
    for (fevt::EyeRecData::PixelIterator pixelIter = eyeRecData.SDPPixelsBegin();
        pixelIter != eyeRecData.SDPPixelsEnd();
        ) {

      pixelIter = eyeRecData.RemoveSDPPixel(pixelIter);

    }

    for (fevt::EyeRecData::PixelIterator pixelIter = eyeRecData.TimeFitPixelsBegin();
        pixelIter != eyeRecData.TimeFitPixelsEnd();
        ++pixelIter) {

      eyeRecData.AddSDPPixel(*pixelIter);

    }

    const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();

    //-------------------------------------------
    // Collect the fit parameters and the errors
    vector<utl::Parameter> fitParameterAndErrors;
    fitParameterAndErrors.push_back(utl::Parameter(fNorthing, fNorthingError));
    fitParameterAndErrors.push_back(utl::Parameter(fEasting, fEastingError));
    fitParameterAndErrors.push_back(utl::Parameter(fAxisTheta, fAxisThetaError));
    fitParameterAndErrors.push_back(utl::Parameter(fAxisPhi, fAxisPhiError));
    fitParameterAndErrors.push_back(utl::Parameter(fTCore, fTCoreError));

    vector<double> fitParameter;
    fitParameter.push_back(fNorthing);
    fitParameter.push_back(fEasting);
    fitParameter.push_back(fAxisTheta);
    fitParameter.push_back(fAxisPhi);
    fitParameter.push_back(fTCore);

    utl::CorrelationMatrix fitCorrelations (5);
    fitCorrelations.Set(0, 1, fCovMatrixAxis[0][1]/(fNorthingError*fEastingError));
    fitCorrelations.Set(0, 2, fCovMatrixAxis[0][2]/(fNorthingError*fAxisThetaError));
    fitCorrelations.Set(0, 3, fCovMatrixAxis[0][3]/(fNorthingError*fAxisPhiError));
    fitCorrelations.Set(0, 4, fCovMatrixAxis[0][4]/(fNorthingError*fTCoreError));
    fitCorrelations.Set(1, 2, fCovMatrixAxis[1][2]/(fEastingError*fAxisThetaError));
    fitCorrelations.Set(1, 3, fCovMatrixAxis[1][3]/(fEastingError*fAxisPhiError));
    fitCorrelations.Set(1, 4, fCovMatrixAxis[1][4]/(fEastingError*fTCoreError));
    fitCorrelations.Set(2, 3, fCovMatrixAxis[2][3]/(fAxisThetaError*fAxisPhiError));
    fitCorrelations.Set(2, 4, fCovMatrixAxis[2][4]/(fAxisThetaError*fTCoreError));
    fitCorrelations.Set(3, 4, fCovMatrixAxis[3][4]/(fAxisPhiError*fTCoreError));

    // Set the EyeRecData
    // These are the values for a telescope which is at the position
    // of the eye.
    const fdet::Eye& detEye = detFD.GetEye(*fgCurEye);
    const utl::Point eyePos = detEye.GetPosition();

    for (fdet::Eye::TelescopeIterator telIter = detEye.TelescopesBegin();
        telIter != detEye.TelescopesEnd();
        ++telIter) {

      const unsigned int telId = telIter->GetId();
      const utl::Point telPos = detEye.GetTelescope(telId).GetPosition();

      if ((telPos - eyePos).GetMag2() < 1.*utl::m) {
        fCurTelId = telId;
        const utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);

        const vector<double> parameters = fTelParameters[telId];
        const double coreX = parameters[0];
        const double coreY = parameters[1];
        const double axisTheta = parameters[2];
        const double axisPhi = parameters[3];
        const double sdpTheta = parameters[4];
        const double sdpPhi = parameters[5];
        const double chiZero = parameters[6];
        const double tZero = parameters[7];
        const double rp = parameters[8];

        //-----------------------------------------
        // Define the object for the error propagation
        const utl::NumericalErrorPropagation errorPropagator(AxisParameterCalculator(this), fitParameterAndErrors, fitCorrelations);

        // Get the propagated errors and correlation matrix
        const utl::CorrelationMatrix axisParameterCorrelations = errorPropagator.GetPropagatedCorrelations();
        const vector<double> axisParameterErrors = errorPropagator.GetPropagatedErrors();

        const double chiZeroError = axisParameterErrors[0];
        const double tZeroError = axisParameterErrors[1];
        const double rpError = axisParameterErrors[2];
        const double sdpThetaError = axisParameterErrors[3];
        const double sdpPhiError = axisParameterErrors[4];

        const double rTZeroRp = axisParameterCorrelations.Get(1, 2);
        const double rChiZeroRp = axisParameterCorrelations.Get(0, 2);
        const double rChiZeroTZero = axisParameterCorrelations.Get(0, 1);
        const double rSDPThetaPhi = axisParameterCorrelations.Get(3, 4);
        // FIXME There are also correlations between Chi0, T0, rp and Theta, Phi

        // Define the core, axis and SDP
        const utl::Point core(coreX, coreY, 0, telCS);
        const utl::Vector axis(1, axisTheta, axisPhi, telCS, utl::Vector::kSpherical);
        const utl::AxialVector sdp(1, sdpTheta, sdpPhi, telCS, utl::AxialVector::kSpherical);

        const utl::Vector telToCore = core - telPos;

        eyeRecData.SetChiZero(chiZero, chiZeroError);
        eyeRecData.SetTZero(tZero, tZeroError);
        eyeRecData.SetRp(rp, rpError);
        eyeRecData.SetTimeFitChiSquare(fChiSqTime, fNDofTime);
        eyeRecData.SetTimeFitCorrelations(rTZeroRp, rChiZeroRp, rChiZeroTZero);
        eyeRecData.SetSDP(sdp);
        eyeRecData.SetSDPThetaError(sdpThetaError);
        eyeRecData.SetSDPPhiError(sdpPhiError);
        eyeRecData.SetSDPCorrThetaPhi(rSDPThetaPhi);
        eyeRecData.SetSDPFitChiSquare(fChiSqSDP, fNDofSDP);
        eyeRecData.SetAxisFitCorrelations(fitCorrelations.Get(0, 1),
                                          fitCorrelations.Get(0, 2),
                                          fitCorrelations.Get(0, 3),
                                          fitCorrelations.Get(0, 4),
                                          fitCorrelations.Get(1, 2),
                                          fitCorrelations.Get(1, 3),
                                          fitCorrelations.Get(1, 4),
                                          fitCorrelations.Get(2, 3),
                                          fitCorrelations.Get(2, 4),
                                          fitCorrelations.Get(3, 4));
        eyeRecData.SetAxisFitChiSquare(fChiSqAxis, fgNDofAxis);

        // Set the ShowerRecData
        if (!fEvent->HasRecShower()) {
          fEvent->MakeRecShower();
        }
        evt::ShowerRecData& recShower = fEvent->GetRecShower();
        recShower.SetAxis(fAxis);
        recShower.SetCorePosition(core);

        // Set the ShowerFRecData
        if (!eyeRecData.HasFRecShower()) {
          eyeRecData.MakeFRecShower();
        }
        evt::ShowerFRecData& frecShower = eyeRecData.GetFRecShower();
        frecShower.SetAxis(fAxis);
        frecShower.SetCorePosition(core);

        // Shift T0 from the pointRp to the core
        double tCore = tZero;
        tCore -= (rp / tan(chiZero)) / utl::kSpeedOfLight;
        utl::TimeStamp timeAtCore = fgCurEye->GetHeader().GetTimeStamp() -
          utl::TimeInterval(1000 * 100 * utl::ns) +
          utl::TimeInterval(tCore);
        frecShower.SetCoreTime(timeAtCore, 0.);

        if (fgStationList.size() > 0) {
          frecShower.SetSDTimeResidual(fFitStationOffset);

          // Only add the station Id if it wasn't done in HybridGeometryFinderOG
          // This is needed for virtual eyes
          if (frecShower.GetStationIds().size() == 0) {
            frecShower.AddStationId(fgFitStation->GetId());
          }
        }

        if (fgVerbosity >= 0) {
          ostringstream info;
          info << "The following data for eye " << fgCurEye->GetId()
            << " is stored:" << endl
            << " Core:     X(eyeCS) = " << coreX / utl::m << endl
            << "           Y(eyeCS) = " << coreY / utl::m << endl
            << "           Northing = " << fNorthing / utl::m << " +/- " << fNorthingError / utl::m << endl
            << "            Easting = " << fEasting / utl::m << " +/- " << fEastingError / utl::m << endl
            << " Axis: Theta(coreCS) = " << fAxisTheta / utl::degree << " +/- " << fAxisThetaError / utl::degree << endl
            << "         Phi(coreCS) = " << fAxisPhi / utl::degree << " +/- " << fAxisPhiError / utl::degree << endl
            << " SDP:  Theta(eyeCS) = " << sdpTheta / utl::degree << " +/- " << sdpThetaError / utl::degree << endl
            << "         Phi(eyeCS) = " << sdpPhi / utl::degree << " +/- " << sdpPhiError / utl::degree << endl
            << "               Chi0 = " << chiZero / utl::degree << " +/- " << chiZeroError / utl::degree << endl
            << "                 T0 = " << tZero / utl::ns << " +/- " << tZeroError / utl::ns << endl
            << "                 Rp = " << rp / utl::m << " +/- " << rpError / utl::m << endl
            << "              TCore = " << tCore / utl::ns << endl
            << " ChiSq / Ndf (Time) = " << fChiSqTime << " / " << fNDofTime << endl
            << "  ChiSq / Ndf (SDP) = " << fChiSqSDP << " / " << fNDofSDP << endl
            << " ChiSq / Ndf (Axis) = " << fChiSqAxis << " / " << fgNDofAxis << endl;
          INFO(info);
        }

        break;
      } // endif telPos = eyePos

    } // for telescopes

    //-------------------------------------------
    // Calculate the axis parameters for each telescope
    for (fevt::Eye::TelescopeIterator telIter = fgCurEye->TelescopesBegin(fevt::ComponentSelector::eHasData);
        telIter != fgCurEye->TelescopesEnd(fevt::ComponentSelector::eHasData);
        ++telIter) {

      const int telId = telIter->GetId();
      fCurTelId = telId;
      const utl::Point telPos = detFD.GetEye(*fgCurEye).GetTelescope(telId).GetPosition();
      const utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);

      const vector<double> parameters = fTelParameters[telId];
      const double coreX = parameters[0];
      const double coreY = parameters[1];
      const double axisTheta = parameters[2];
      const double axisPhi = parameters[3];
      const double sdpTheta = parameters[4];
      const double sdpPhi = parameters[5];
      const double chiZero = parameters[6];
      const double tZero = parameters[7];
      const double rp = parameters[8];

      //-----------------------------------------
      // Define the object for the error propagation
      const utl::NumericalErrorPropagation errorPropagator(AxisParameterCalculator(this), fitParameterAndErrors, fitCorrelations);

      // Get the propagated errors and correlation matrix
      const utl::CorrelationMatrix axisParameterCorrelations = errorPropagator.GetPropagatedCorrelations();
      const vector<double> axisParameterErrors = errorPropagator.GetPropagatedErrors();

      const double chiZeroError = axisParameterErrors[0];
      const double tZeroError = axisParameterErrors[1];
      const double rpError = axisParameterErrors[2];
      const double sdpThetaError = axisParameterErrors[3];
      const double sdpPhiError = axisParameterErrors[4];

      const double rTZeroRp = axisParameterCorrelations.Get(1, 2);
      const double rChiZeroRp = axisParameterCorrelations.Get(0, 2);
      const double rChiZeroTZero = axisParameterCorrelations.Get(0, 1);
      const double rSDPThetaPhi = axisParameterCorrelations.Get(3, 4);
      // FIXME There are also correlations between Chi0, T0, rp and Theta, Phi

      // Define the core, axis and SDP
      const utl::Point core(coreX, coreY, 0, telCS);
      const utl::Vector axis(1, axisTheta, axisPhi, telCS, utl::Vector::kSpherical);
      const utl::AxialVector sdp(1, sdpTheta, sdpPhi, telCS, utl::AxialVector::kSpherical);

      const utl::Vector vertical(0, 0, 1, telCS, utl::Vector::kCartesian);
      const utl::AxialVector horizontalInSDP = utl::cross(sdp, vertical);

      // Recalculate Chi_i for each pulsed pixel in this telescope
      for (fevt::EyeRecData::PixelIterator pixelIter = fgCurEye->GetRecData().PulsedPixelsBegin();
          pixelIter != fgCurEye->GetRecData().PulsedPixelsEnd();
          ++pixelIter) {

        if (int(pixelIter->GetTelescopeId()) != telId) {
          continue;
        }

        const utl::Vector pixelDirection = detFD.GetPixel(*pixelIter).GetDirection();
        const utl::Vector pixelDirInSDP = pixelDirection - (sdp * pixelDirection) * sdp;
        const double pixelChi = utl::Angle(pixelDirInSDP, horizontalInSDP);

        pixelIter->GetRecData().SetChi_i(pixelChi);

      } // for pixel

      // Set the telescope RecData
      if (!telIter->HasRecData()) {
        telIter->MakeRecData();
      }
      fevt::TelescopeRecData& telRecData = telIter->GetRecData();

      telRecData.SetChiZero(chiZero, chiZeroError);
      telRecData.SetTZero(tZero, tZeroError);
      telRecData.SetRp(rp, rpError);
      telRecData.SetTimeFitChiSquare(fChiSqTimePerTel[telId], fNDofTimePerTel[telId]);
      telRecData.SetTimeFitCorrelations(rTZeroRp, rChiZeroRp, rChiZeroTZero);
      telRecData.SetSDP(sdp);
      telRecData.SetSDPThetaError(sdpThetaError);
      telRecData.SetSDPPhiError(sdpPhiError);
      telRecData.SetSDPCorrThetaPhi(rSDPThetaPhi);
      telRecData.SetSDPFitChiSquare(fChiSqSDPPerTel[telId], fNDofSDPPerTel[telId]);
      telRecData.SetAxisFitCorrelations(fitCorrelations.Get(0, 1),
                                        fitCorrelations.Get(0, 2),
                                        fitCorrelations.Get(0, 3),
                                        fitCorrelations.Get(0, 4),
                                        fitCorrelations.Get(1, 2),
                                        fitCorrelations.Get(1, 3),
                                        fitCorrelations.Get(1, 4),
                                        fitCorrelations.Get(2, 3),
                                        fitCorrelations.Get(2, 4),
                                        fitCorrelations.Get(3, 4));
      telRecData.SetAxisFitChiSquare(fChiSqAxisPerTel[telId], fNDofAxisPerTel[telId]);

      if (fgVerbosity >= 0) {
        ostringstream info;
        info << "The following data for telescope " << telId
             << " is stored:" << endl
             << " Core:     X(telCS) = " << coreX / utl::m << endl
             << "           Y(telCS) = " << coreY / utl::m << endl
             << "           Northing = " << fNorthing / utl::m << " +/- " << fNorthingError / utl::m << endl
             << "            Easting = " << fEasting / utl::m << " +/- " << fEastingError / utl::m << endl
             << " Axis: Theta(coreCS) = " << fAxisTheta / utl::degree << " +/- " << fAxisThetaError / utl::degree << endl
             << "         Phi(coreCS) = " << fAxisPhi / utl::degree << " +/- " << fAxisPhiError / utl::degree << endl
             << " SDP:  Theta(telCS) = " << sdpTheta / utl::degree << " +/- " << sdpThetaError / utl::degree << endl
             << "         Phi(telCS) = " << sdpPhi / utl::degree << " +/- " << sdpPhiError / utl::degree << endl
             << "               Chi0 = " << chiZero / utl::degree << " +/- " << chiZeroError / utl::degree << endl
             << "                 T0 = " << tZero / utl::ns << " +/- " << tZeroError / utl::ns << endl
             << "                 Rp = " << rp / utl::m << " +/- " << rpError / utl::m << endl
             << " ChiSq / Ndf (Time) = " << fChiSqTimePerTel[telId] << " / " << fNDofTimePerTel[telId] << endl
             << "  ChiSq / Ndf (SDP) = " << fChiSqSDPPerTel[telId] << " / " << fNDofSDPPerTel[telId] << endl
             << " ChiSq / Ndf (Axis) = " << fChiSqAxisPerTel[telId] << " / " << fNDofAxisPerTel[telId] << endl;
        INFO(info);
      }

    } // for telescopes

  } // void HybridGeometryFinder::StoreData()


  vector<double>
  AxisParameterCalculator::operator() (const vector<double>& input)
    const
  {

    const double northing = input[0];
    const double easting = input[1];
    const double axisTheta = input[2];
    const double axisPhi = input[3];
    const double tCore = input[4];

    //-------------------------------------------
    // Create the coordinate system at the core
    const utl::CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
    const utl::ReferenceEllipsoid& refEllipsoid = utl::ReferenceEllipsoid::GetWGS84();

    // The core is defined at fgCoreHeight (1400m)
    const utl::UTMPoint UTMCore(northing, easting, fHybridFinder->fgCoreHeight, 19, 'H', refEllipsoid);

    // Create a coordinate system at the core point
    const utl::Point core = UTMCore.GetPoint(siteCS);
    const utl::CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);

    // The axis is defined in the coreCS
    const utl::Vector axis(1, axisTheta, axisPhi, coreCS, utl::Vector::kSpherical);
    const bool downGoing = axisTheta < utl::kPiOnTwo;

    // Define the telCS
    const fevt::Eye& eye = *(fHybridFinder->fgCurEye);
    const unsigned int telId = fHybridFinder->fCurTelId;
    const utl::Point telPos = det::Detector::GetInstance().GetFDetector().GetEye(eye).GetTelescope(telId).GetPosition();
    const utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);

    // Shift the core along the axis such that it is at the same height as the telescope
    const utl::Point telCore = core - axis * core.GetZ(telCS) / axis.GetZ(telCS);

    // Get the SDP for this telescope
    const utl::Vector telToCore = telCore - telPos;
    utl::AxialVector sdp;
    sdp = utl::cross(axis, telToCore);
    sdp.Normalize();

    // Get the Chi0
    double chiZero = utl::Angle(axis, telToCore);

    // Get Rp
    const double rp = telToCore.GetMag() * sin(chiZero);

    // For up-going showers Chi0 < 0
    if (!downGoing) {
      chiZero = utl::kPi - chiZero;
      chiZero *= -1;
    }

    // Get the Point of Rp
    utl::Vector telToRp = utl::cross(sdp, axis);
    telToRp.Normalize();
    telToRp *= rp;
    const utl::Point pointRp = telPos + telToRp;

    // Shift T0 from the core to pointRp
    double tZero = tCore;
    if (pointRp.GetZ(telCS) > core.GetZ(telCS)) {
      if (downGoing) {
        tZero -= (pointRp - core).GetMag() / utl::kSpeedOfLight;
      }
      else {
        tZero += (pointRp - core).GetMag() / utl::kSpeedOfLight;
      }
    }
    else {
      if (downGoing) {
        tZero += (pointRp - core).GetMag() / utl::kSpeedOfLight;
      }
      else {
        tZero -= (pointRp - core).GetMag() / utl::kSpeedOfLight;
      }
    }

    vector<double> output;

    output.push_back(chiZero);
    output.push_back(tZero);
    output.push_back(rp);
    output.push_back(sdp.GetTheta(telCS));
    output.push_back(sdp.GetPhi(telCS));

    return output;

  } // vector<double> AxisParameterCalculator::operator() (const vector<double>& input)

} // namespace HybridGeometryFinderWG