RdPlaneFit.cc

Go to the documentation of this file.

#include "RdPlaneFit.h"

#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/ErrorLogger.h>
#include <utl/AnalyticCoordinateTransformator.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerRRecDataQuantities.h>

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

#include <CLHEP/Matrix/Matrix.h>
#include <CLHEP/Matrix/Vector.h>

#include <TMinuit.h>

#include <cmath>

using namespace std;
using namespace revt;
using namespace evt;
using namespace fwk;
using namespace utl;
using CLHEP::HepMatrix;
using CLHEP::HepVector;


namespace RdPlaneFit {

  // global variables
  REvent* RdPlaneFit::fgCurrentREvent;
  Point RdPlaneFit::fgBarycenter;
  TimeStamp RdPlaneFit::fgBaryTime;
  CoordinateSystemPtr RdPlaneFit::fgLocalCS;
  Point RdPlaneFit::fgCoordinateOrigin;
  double RdPlaneFit::fgMeanEventTime;


  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();
  }


  VModule::ResultFlag
  RdPlaneFit::Init()
  {
    // Read in the configurations of the xml file
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdPlaneFit");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

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

    topBranch.GetChild("minNumberOfStations").GetData(fMinNumberOfStations);
    topBranch.GetChild("allowUnphysicalCosines").GetData(fAllowUnphysicalCosines);

    const string yesnostr = topBranch.GetChild("allowMultipleCalls").Get<string>();
    fLockParameterStorage = (yesnostr == "no");

    return eSuccess;
  }


  // analytic linear fit in local CS; approximation by neglecting
  // z-coordinate, since z ~ 0 in local bary coordinate system,
  // therefore loosing the w component of the axis. it is funny
  // though, that stations living on a plane can not distinguish
  // between upwards or downwards-going showers.

  VModule::ResultFlag
  RdPlaneFit::LinearFit(FitParameters& fit, const bool reuseFit)
    const
  {
    //const double cosTheta = (reuseFit ? fit.w : 1);
    const Vector approx_axis = (reuseFit ? Vector(fit.u, fit.v, fit.w, fgLocalCS) : Vector());
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    double s1 = 0;
    double sx = 0;
    double sy = 0;
    double st = 0;
    double sxx = 0;
    double sxy = 0;
    double syy = 0;
    double sxt = 0;
    double syt = 0;
    //double stt = 0;

    // loop over all stations
    for (auto sIt = fgCurrentREvent->SignalStationsBegin();
        sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {
      const Vector sPos = rDetector.GetStation(*sIt).GetPosition() - fgBarycenter;
      const double x = sPos.GetX(fgLocalCS);
      const double y = sPos.GetY(fgLocalCS);
      const StationRecData& sRec = sIt->GetRecData();
      //const double t = (sRec.GetSignalTime() - fgBaryTime).GetInterval();
      const double t = sRec.GetParameter(eSignalTime) - fgMeanEventTime;
      //const double invSigma2 = 1 / (sRec.GetSignalTimeError() * sRec.GetSignalTimeError());
      const double invSigma2 = 1 / (sRec.GetParameterError(eSignalTime) * sRec.GetParameterError(eSignalTime));
      s1 += invSigma2;
      const double x_s = x * invSigma2;
      sx += x_s;
      sxx += x * x_s;
      sxy += y * x_s;
      sxt += t * x_s;
      const double y_s = y * invSigma2;
      sy += y_s;
      syy += y * y_s;
      syt += t * y_s;
      st += t * invSigma2;
      //stt += t*t;
    }

    st *= kSpeedOfLight;
    sxt *= kSpeedOfLight;
    syt *= kSpeedOfLight;
    //stt *= kSpeedOfLight*kSpeedOfLight;

    HepMatrix mat(3, 3);

    mat[0][0] = sxx;
    mat[0][1] = sxy;
    mat[0][2] = -sx;
    mat[1][0] = sxy;
    mat[1][1] = syy;
    mat[1][2] = -sy;
    mat[2][0] = -sx;
    mat[2][1] = -sy;
    mat[2][2] = s1;

    int ierr = 0;
    mat.invert(ierr);

    if (ierr) {
      INFO("Matrix inversion failed.");
      return eFailure;
    }

    HepVector k(3);

    k[0] = -sxt;
    k[1] = -syt;
    k[2] = st;

    const HepVector p = mat * k;
    double w2 = 1 - Sqr(p[0]) - Sqr(p[1]);

    if (w2 < 0) {
      WARNING("Unphysical direction cosines");
      if (fAllowUnphysicalCosines) {
        WARNING("Setting z-compenent of the shower axis to 0");
        w2 = 0;
      } else
        return eFailure;
    }

    fit.u = p[0];
    fit.v = p[1];
    fit.w = sqrt(w2);
    fit.ct0 = p[2];

    fit.sigmaU2 = mat[0][0];
    fit.sigmaV2 = mat[1][1];
    fit.sigmaUV = mat[0][1];
    fit.sigmact02 = mat[2][2];

    return eSuccess;
  }


  VModule::ResultFlag
  RdPlaneFit::Run(evt::Event& event)
  {
    INFO("Plane fit");

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

    fgCurrentREvent = &event.GetREvent();
    const REvent& rEvent = *fgCurrentREvent;

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

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

    try {
      fgCoordinateOrigin = showerrrec.GetCoordinateOrigin();
      fgLocalCS = LocalCoordinateSystem::Create(fgCoordinateOrigin);
    } catch (utl::NonExistentComponentException& e) {
      ostringstream err;
      err << "Problem when getting coordinate origin \n";
      err << "Did you call the RdEventInitializer ? \n";
      err << " Caught utl::NonExistentComponentException : " << e.what();
      err << "\n Will Break loop \n";
      ERROR(err.str());
      return eBreakLoop;
    }

    const TimeStamp& eventTime = rEvent.GetHeader().GetTime();

    // weighted sums
    Vector barySum(0, 0, 0, fgLocalCS);
    TimeInterval timeSum = 0;
    double weightSum = 0;
    int nStations = 0;
    for (auto sIt = fgCurrentREvent->SignalStationsBegin();
        sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {
      const rdet::Station& dStation = rDetector.GetStation(*sIt);
      const StationRecData& sRec = sIt->GetRecData();
      const double weight = sqrt(sRec.GetParameter(eSignal));
      weightSum += weight;
      barySum += weight * (dStation.GetPosition() - fgCoordinateOrigin);
      timeSum += weight * sRec.GetParameter(eSignalTime);
      ++nStations;
    }

    ostringstream info;
    info << "Number of signal stations = " << nStations;
    INFOFinal(info);

    // not worth it
    if (nStations < fMinNumberOfStations) {
      info.str("");
      info << "eContinueLoop: Not enough stations." << endl;
      if (event.HasRecShower()) {
        if (event.GetRecShower().HasRRecShower()) {
          event.GetRecShower().GetRRecShower().SetParameter(eFitSuccess, 0, fLockParameterStorage);
          info << "Set FitSuccess to false.";
        }
      }
      INFOFinal(info);
      return eContinueLoop;
    }

    if (!weightSum) {
      info.str("");
      info << "eContinueLoop: sum of weights = 0" << endl;
      if (event.HasRecShower()) {
        if (event.GetRecShower().HasRRecShower()) {
          event.GetRecShower().GetRRecShower().SetParameter(eFitSuccess, 0, fLockParameterStorage);
          info << "Set FitSuccess to false.";
        }
      }
      INFOIntermediate(info);
      return eContinueLoop;
    }

    barySum /= weightSum;
    timeSum /= weightSum;
    fgMeanEventTime = timeSum;
    fgBarycenter = barySum + fgCoordinateOrigin;
    fgBaryTime = eventTime + timeSum;

    // site coordinate system (usually Malargue), only you used for terminal output
    const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();
    // get reference coordinate system of detector (usually PampaAmarilla)
    const CoordinateSystemPtr referenceCS = detector.GetReferenceCoordinateSystem();
    info.str("");
    info << "barycenter = " << ToString(fgBarycenter, siteCS) << " (site)" << std::endl
         << "bary time = " << timeSum / nanosecond << " [ns] (to event time)\n"
         << "local/site zenith angle diff. = "
         << acos(Vector(0, 0, 1, fgLocalCS) * Vector(0, 0, 1, siteCS)) / degree << " [deg]";
    INFOIntermediate(info);

    // this will point to the final result
    FitParameters finalFitResult;

    // first approximation
    FitParameters fitLin;
    if (LinearFit(fitLin) != eSuccess) {
      INFOIntermediate("eContinueLoop: Linear fit was not successful");
      if (event.HasRecShower()) {
        if (event.GetRecShower().HasRRecShower()) {
          event.GetRecShower().GetRRecShower().SetParameter(eFitSuccess, 0, fLockParameterStorage);
          INFOIntermediate("Set FitSuccess to false.");
        }
      }
      return eContinueLoop;
    }

    // try again with previous values as approxiamtions for sigma
    LinearFit(fitLin, true);

    // with first approximation on u and v, the nonlinear fit is attempted
    FitParameters fitNonlin = fitLin;

    bool nonLinearFitUsed = false;

    info.str("");
    if (PlaneFit3DDriver(fitNonlin) == eSuccess) {
      //fitNonlin.stage = ShowerSRecData::kPlaneFit3d;
      info << "\n\tStage " << fitLin.stage << ": linear plane fit"
           << "\n\taxis = (" << fitLin.u << ", " << fitLin.v << ", " << fitLin.w << ") (bary)"
           << "\n\ttheta = " << acos(fitLin.w) / degree
           << "\n\tphi   = " << atan2(fitLin.v, fitLin.u) / degree
           << "\n\tct0 = " << fitLin.ct0 / meter << " [m]";
      nonLinearFitUsed = true;
      finalFitResult = fitNonlin;
    } else {
      info << "  Nonlinear fit failed, using linear approximation.";
      finalFitResult = fitLin;
    }
    INFOIntermediate(info);


    // Fill the event
    const Vector axis(finalFitResult.u, finalFitResult.v, finalFitResult.w, fgLocalCS);
    showerrrec.SetParameter(eShowerAxisX, axis.GetX(fgLocalCS), fLockParameterStorage);
    showerrrec.SetParameter(eShowerAxisY, axis.GetY(fgLocalCS), fLockParameterStorage);
    showerrrec.SetParameter(eShowerAxisZ, axis.GetZ(fgLocalCS), fLockParameterStorage);

    // save Covariance matrix, copied from RdWaveFit
    double thetaError;
    double phiError;
    double thetaPhiCorrelation;
    PropagateAxisErrors(axis.GetCoordinates(fgLocalCS),
                        Vector::Triple(finalFitResult.sigmaU2, finalFitResult.sigmaUV, finalFitResult.sigmaV2),
                        thetaError, phiError, thetaPhiCorrelation);
    int r = 1;
    double theta = acos(finalFitResult.w);
    double phi = atan2(finalFitResult.v, finalFitResult.u);
    int c00 = 0;
    int c01 = 0;
    int c02 = 0; // radius not fitted
    double c11 = thetaError * thetaError;
    double c22 = phiError * phiError;
    double c12 = thetaPhiCorrelation * thetaError * phiError;

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

    showerrrec.SetParameter(eCoreTime, fgMeanEventTime + finalFitResult.ct0 / kSpeedOfLight, fLockParameterStorage);
    showerrrec.SetParameterError(eCoreTime, fgMeanEventTime + sqrt(finalFitResult.sigmact02), fLockParameterStorage);

    // set signal arrival direction for individual stations
    // For a plane wave fit the direction is of course the same for all station, but
    // for some applications the RdAntennaChannelToStationConverter has to know the individual
    // arrival directions. So the values have to be set.
    for (auto sIt = fgCurrentREvent->SignalStationsBegin();
        sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {
      Station& station = *sIt;
      station.GetRecData().SetParameter(eSignalArrivalDirectionX, axis.GetX(fgLocalCS), fLockParameterStorage);
      station.GetRecData().SetParameter(eSignalArrivalDirectionY, axis.GetY(fgLocalCS), fLockParameterStorage);
      station.GetRecData().SetParameter(eSignalArrivalDirectionZ, axis.GetZ(fgLocalCS), fLockParameterStorage);
    }

    double chi2 = 0;
    // calculates residuals and fill them in station rec
    CalculateTimeResidual(axis, finalFitResult.ct0, chi2);

    const int angleNDOF = nStations - 3;
    showerrrec.SetParameter(eWaveFitNDF, angleNDOF, fLockParameterStorage);
    showerrrec.SetParameter(eWaveFitChi2, chi2, fLockParameterStorage);
    showerrrec.SetParameter(eFitSuccess, 1, fLockParameterStorage);

    info.str("");
    info << "\n\tStage " << finalFitResult.stage << ": " << (nonLinearFitUsed ? "3d" : "linear") << " plane fit"
         << "\n\taxis  = (" << finalFitResult.u << ", " << finalFitResult.v << ", " << finalFitResult.w << ") (bary)"
         << "\n\ttheta = " << acos(finalFitResult.w) / degree << " +/- "
         << showerrrec.GetZenithError()/degree << " [deg] (bary)"
         << "\n\tphi   = " << atan2(finalFitResult.v, finalFitResult.u) / degree << " +/- "
         << showerrrec.GetAzimuthError()/degree << " [deg]"
         << "\n\tct0   = " << finalFitResult.ct0 / meter << " [m]";
    INFOFinal(info);

    return eSuccess;
  }


  VModule::ResultFlag
  RdPlaneFit::PlaneFit3DDriver(FitParameters& fit)
    const
  {
    TMinuit minuit(3);
    int errFlag = 0;
    double argList[10];
    argList[0] = -1;
    if (!fMinuitOutput) {
      minuit.mnexcm("SET PRINTOUT", argList, 1, errFlag);
      minuit.mnexcm("SET NOWARNINGS", argList, 0, errFlag);
      minuit.SetPrintLevel(-1);
    }

    if (fAllowUnphysicalCosines) {
      minuit.SetFCN(RdPlaneFit::PlaneFit3DHorizonFnc);
      argList[0] = 1; // chi2
      minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);
      minuit.mnparm(0, "theta", acos(fit.w), 0.01, 0, 0, errFlag);
      minuit.mnparm(1, "phi", atan2(fit.v, fit.u), 0.01, 0, 0, errFlag);
      minuit.mnparm(2, "ct0", fit.ct0, 10 * meter, 0, 0, errFlag);
    } else {
      minuit.SetFCN(RdPlaneFit::PlaneFit3DFnc);
      argList[0] = 1; // chi2
      minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);
      minuit.mnparm(0, "u", fit.u, 0.01, 0, 0, errFlag);
      minuit.mnparm(1, "v", fit.v, 0.01, 0, 0, errFlag);
      minuit.mnparm(2, "ct0", fit.ct0, 10 * meter, 0, 0, errFlag);
    }

    ostringstream info;

    // init fnc
    argList[0] = 1;
    minuit.mnexcm("CALI", argList, 1, errFlag);

    argList[0] = 500;
    minuit.mnexcm("MINIMIZE", argList, 1, errFlag);
    if (errFlag) {
      info.str("");
      info << "minuit minimize error: " << errFlag << '\n';
      INFOIntermediate(info);

      minuit.mnexcm("MINOS", argList, 0, errFlag);
      if (errFlag) {
        info.str("");
        info << "minuit minos error: " << errFlag << '\n';
        INFOIntermediate(info);

        return eFailure;
      }
    }

    // Get results
    double foo = 0;
    double par1 = 0;
    double par2 = 0;

    minuit.GetParameter(0, par1, foo);
    minuit.GetParameter(1, par2, foo);

    if (fAllowUnphysicalCosines) {
      fit.u = sin(par1) * cos(par2);
      fit.v = sin(par1) * sin(par2);
      fit.w = cos(par1);
    } else {
      fit.u = par1;
      fit.v = par2;
      double w2 = 1 - Sqr(fit.u) - Sqr(fit.v);
      if (w2 < 0) {
        info.str("");
        info << "eFailure: Fit failed due to unphysical angle cosines";
        INFOFinal(info);

        return eFailure;
      }
      fit.w = sqrt(w2);
    }

    minuit.GetParameter(2, fit.ct0, foo);

    // error matrix
    double emat[3][3];
    minuit.mnemat(&emat[0][0], 3);
    fit.sigmaU2 = emat[0][0];
    fit.sigmaV2 = emat[1][1];
    fit.sigmaUV = emat[0][1];
    fit.sigmact02 = emat[2][2];

    return eSuccess;
  }


  VModule::ResultFlag
  RdPlaneFit::Finish()
  {
    INFODebug("");
    return eSuccess;
  }


  void
  RdPlaneFit::PlaneFit3DFnc(int& /*nPar*/, double* const /*grad*/, double& value, double* const par, const int flag)
  {
    const double& u = par[0];
    const double& v = par[1];
    const double& ct0 = par[2];
    static int nStations = 0;
    static vector<double> sTiming;     // station timing
    static vector<double> sSignal;     // station signal
    static vector<Vector> sPosition;   // station position
    static vector<double> sTimeSigma2; // Uncurtainty on the pulse timing
    //static const sdet::STimeVariance* timeVariance = 0;

    if (flag == 1) {  // fnc init

      sTiming.clear();
      sSignal.clear();
      sPosition.clear();
      sTimeSigma2.clear();

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

      for (auto sIt = fgCurrentREvent->SignalStationsBegin();
           sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {
        const StationRecData& sRec = sIt->GetRecData();
        //sTiming.push_back(kSpeedOfLight * (sRec.GetSignalTime() - fgBaryTime).GetInterval());
        sTiming.push_back(kSpeedOfLight * (sRec.GetParameter(eSignalTime) - fgMeanEventTime));
        sSignal.push_back(sRec.GetParameter(eSignal));
        sPosition.push_back(rDetector.GetStation(*sIt).GetPosition() - fgBarycenter);
        sTimeSigma2.push_back(sRec.GetParameterError(eSignalTime) * sRec.GetParameterError(eSignalTime));
      }

      nStations = sTiming.size();

      //timeVariance = &sdet::STimeVariance::GetInstance();

    } // fnc init

    const double w2 = 1 - Sqr(u) - Sqr(v);

    if (w2 < 0) {
      value = 1e30;
      return;
    }

    const double w = sqrt(w2);

    const Vector axis(u, v, w, fgLocalCS);

    double chi2 = 0;

    for (int i = 0; i < nStations; ++i)
      chi2 += Sqr(sTiming[i] - ct0 + axis * sPosition[i]) / sTimeSigma2[i];

    value = chi2 / kSpeedOfLight2;
  }


  void
  RdPlaneFit::PlaneFit3DHorizonFnc(int& /*nPar*/, double* const /*grad*/, double& value,
                                   double* const par, const int flag)
  {
    const double& theta = par[0];
    const double& phi = par[1];
    const double& ct0 = par[2];
    static int nStations = 0;
    static vector<double> sTiming;     // station timing
    static vector<double> sSignal;     // station signal
    static vector<Vector> sPosition;   // station position
    static vector<double> sTimeSigma2; // Uncurtainty on the pulse timing
    //static const sdet::STimeVariance* timeVariance = 0;

    if (flag == 1) {  // fnc init

      sTiming.clear();
      sSignal.clear();
      sPosition.clear();
      sTimeSigma2.clear();

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

      for (auto sIt = fgCurrentREvent->SignalStationsBegin();
           sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {
        const StationRecData& sRec = sIt->GetRecData();
        sTiming.push_back(kSpeedOfLight * (sRec.GetParameter(eSignalTime) - fgMeanEventTime));
        sSignal.push_back(sRec.GetParameter(eSignal));
        sPosition.push_back(rDetector.GetStation(*sIt).GetPosition() - fgBarycenter);
        sTimeSigma2.push_back((sRec.GetParameterError(eSignalTime) * sRec.GetParameterError(eSignalTime)));
      }

      nStations = sTiming.size();

      //timeVariance = &sdet::STimeVariance::GetInstance();

    } // fnc init

    const Vector axis(cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta), fgLocalCS);

    double chi2 = 0;

    for (int i = 0; i < nStations; ++i)
      chi2 += Sqr(sTiming[i] - ct0 + axis * sPosition[i]) / sTimeSigma2[i];

    value = chi2 / kSpeedOfLight2;
  }


  void
  RdPlaneFit::CalculateTimeResidual(const utl::Vector& axis, const double ct0, double& chi2)
    const
  {
    const double t0 = ct0 / kSpeedOfLight;
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
    chi2 = 0;
    for (auto sIt = fgCurrentREvent->SignalStationsBegin();
         sIt != fgCurrentREvent->SignalStationsEnd(); ++sIt) {
      StationRecData& sRec = sIt->GetRecData();
      const Vector stLoc_c = (rDetector.GetStation(*sIt).GetPosition() - fgBarycenter) / kSpeedOfLight;
      const double t = (sRec.GetParameter(eSignalTime) - fgMeanEventTime);
      double sigma_t =sRec.GetParameterError(eSignalTime);
      double timeResidual = t - t0 + (stLoc_c * axis);
      sRec.SetParameter(eTimeResidual, timeResidual, fLockParameterStorage);
      sRec.SetParameterError(eTimeResidual, sigma_t, fLockParameterStorage);
      chi2 += pow(timeResidual / sigma_t,2);
    }
  }

}