ZombieGraveyard/LDFFinderOG/LDFFinder.cc

Go to the documentation of this file.

#include <cmath>

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

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

#include <fevt/FEvent.h>
#include <fevt/Header.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeRecData.h>

#include <det/Detector.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <sdet/STimeVariance.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>

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

#include <utl/CoordinateSystem.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/ErrorLogger.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/TimeStamp.h>
#include <utl/AxialVector.h>
#include <utl/Math.h>
#include <utl/Reader.h>
#include <utl/Accumulator.h>

#include <TMinuit.h>

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

#include "RoptFit.h"
#include "LDFFinder.h"

using namespace LDFFinderOG;
using namespace sevt;
using namespace fevt;
using namespace evt;
using namespace fwk;
using namespace std;
using namespace utl;

using CLHEP::HepMatrix;
using CLHEP::HepVector;
using CLHEP::HepRotation;


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


const SEvent* LDFFinder::fgCurrentSEvent;
Point LDFFinder::fgBarycenter;
TimeStamp LDFFinder::fgBaryTime;
Vector LDFFinder::fgShowerAxis;
CoordinateSystemPtr LDFFinder::fgLocalCS;
LDFFunctionType LDFFinder::fgLDFType;
bool LDFFinder::fgUseSilentStations;
double LDFFinder::fgSilentMaxRadius;
double LDFFinder::fgSilentRadiusTransition;
double LDFFinder::fgSilentSignalThreshold;
bool LDFFinder::fgFixLDF;
double LDFFinder::fgNoRecovery;
bool LDFFinder::fgLowerLimit;
double LDFFinder::fgNKGFermiMu;
double LDFFinder::fgNKGFermiTau;

namespace LDFFinderOG {

  // square of the perpendicular distance from the axis
  inline
  double
  RPerp2(const Vector& axis, const Vector& station)
  {
    const double scal = axis*station;
    return station*station - scal*scal;
  }


  // perpendicular distance from the axis
  inline
  double
  RPerp(const Vector& axis, const Vector& station)
  {
    return sqrt(RPerp2(axis, station));
  }


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

}


LDFFinder::LDFFitInterface::LDFFitInterface() :
  fS1000(0), fS1000Err(0),
  // fCore & fCoreErr ctor() ok
  fFixedBeta(false), fBeta(0), fBetaErr(0),
  fFixedGamma(false), fGamma(0), fGammaErr(0),
  fEnergy(0), fEnergyErr(0),
  fChi2(0), fNdof(0), fRecStage(0)
{ }


LDFFinder::CurvatureFitInterface::CurvatureFitInterface() :
  // fCore
  fFixedAxis(false),
  // fAxis, fAxisErr
  // fSigmaU2, fSigmaV2, fSigmaUV
  fFixedRc(false),
  fRc(0), fRcErr(0),
  fct0(0), fct0Err(0)
{ }


void
LDFFinder::OutputStations(const LDFFitInterface& fi)
  const
{
  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

  // loop on stations
  for (SEvent::ConstStationIterator sIt = fgCurrentSEvent->StationsBegin();
       sIt != fgCurrentSEvent->StationsEnd(); ++sIt)
    if (sIt->IsCandidate()) {

      const double r =
        RPerp(fgShowerAxis, sDetector.GetStation(*sIt).GetPosition() - fgBarycenter);

      cout << "  candidate " << sIt->GetId() << ": " << r << ' '
           << sIt->GetRecData().GetTotalSignal() << ' '
           << fi.fS1000 * LDFFunction(fgLDFType, r, fi.fBeta, fi.fGamma, fgNKGFermiMu, fgNKGFermiTau) << '\n';

    } else if (sIt->IsSilent()) {

      const double r =
        RPerp(fgShowerAxis, sDetector.GetStation(*sIt).GetPosition() - fgBarycenter);

      cout << "  silent " << sIt->GetId() << ": " << r << ' '
           << fi.fS1000 * LDFFunction(fgLDFType, r, fi.fBeta, fi.fGamma, fgNKGFermiMu, fgNKGFermiTau) << '\n';
    }
}


void
LDFFinder::OutputResults(const Event& event)
  const
{
  const ShowerSRecData& recShower = event.GetRecShower().GetSRecShower();

  OUT(eFinal)
    << "***** Final reconstructed shower data *****\n"
       "  LDF stage             = " << recShower.GetLDFRecStage() << "\n"
       "  S1000                 = " << recShower.GetShowerSize() << " +/- " << recShower.GetShowerSizeError() << "\n"
       "  core                  = " << ToString(recShower.GetCorePosition(),
                det::Detector::GetInstance().GetSiteCoordinateSystem()) << " (site)\n"
       "                        = " << ToString(recShower.GetCorePosition(), fgLocalCS) << " (bary) +/-\n"
       "                          " << ToString(recShower.GetCoreError(), fgLocalCS) << "\n"
       "  x-y correlation       = " << recShower.GetCorrelationXY() << "\n"
       "  beta                  = " << recShower.GetBeta() << " +/- " << recShower.GetBetaError() << "\n"
       "  gamma                 = " << recShower.GetGamma() << " +/- " << recShower.GetGammaError() << "\n"
       "  energy                = " << recShower.GetEnergy()/EeV << " +/- " << recShower.GetEnergyError()/EeV
    << " [EeV]\n";
  const double curv = recShower.GetCurvature();
  if (curv) {
    const Vector& axis = recShower.GetAxis();
    const double rc = 1/curv;
    const Vector::Triple uvw(axis.GetCoordinates(fgLocalCS));
    OUT(eFinal)
      << "axis                  = (" << uvw.get<0>() << ", " << uvw.get<1>() << ", " << uvw.get<2>()
                      << ") (bary)" << "\n"
         "  theta                 = " << axis.GetTheta(fgLocalCS)/degree << " +/- "
                     << recShower.GetThetaError()/degree << " [deg] (bary)\n"
         "  phi                   = " << axis.GetPhi(fgLocalCS)/degree << " +/- "
                     << recShower.GetPhiError()/degree << " [deg]\n"
         "  theta-phi correlation = " << recShower.GetCorrelationThetaPhi() << "\n"
         "  Rc                    = " << rc/km << " +- " << recShower.GetCurvatureError()*rc*rc/km << " [km]\n"
         "  time chi2             = " << recShower.GetAngleChi2() << " / "
                            << recShower.GetAngleNdof() << "\n"
         "  time residual         = " << recShower.GetTimeResidualMean()/nanosecond << " +/- "
                            << recShower.GetTimeResidualSpread()/nanosecond << " [ns]\n";
  }
  const double rOpt = recShower.GetLDFOptimumDistance();
  if (rOpt)
    OUT(eFinal)
      << "Ropt   = " << rOpt/meter << " [m]\n";
  cerr << flush;
}


VModule::ResultFlag
LDFFinder::Init()
{
  CentralConfig* const cc = CentralConfig::GetInstance();

  const Branch topB = cc->GetTopBranch("LDFFinder");

  topB.GetChild("infoLevel").GetData(fInfoLevel);
  if (fInfoLevel < eNone || fInfoLevel > eMinuit) {
    ERROR("<infoLevel> out of bounds");
    return eFailure;
  }
  fMinuitOutput = (fInfoLevel >= eMinuit);

  string tmp;
  topB.GetChild("coreType").GetData(tmp);
  fCoreType = (tmp == "MC" ? eMC : (tmp == "Hybrid" ? eHybrid : eFit));

  topB.GetChild("maxTheta").GetData(fMaxTheta);
  if (fMaxTheta <= 0 || fMaxTheta > 90*degree) {
    ERROR("<maxTheta> does not make sense");
    return eFailure;
  }

  //  string tmp;
  topB.GetChild("ldfType").GetData(tmp);
  fgLDFType = (tmp == "PL" ? ePL : (tmp == "NKG") ? eNKG : eNKGFermi);

  topB.GetChild("NKGFermiMu").GetData(fgNKGFermiMu);
  if (fgNKGFermiMu <= 0) {
    ERROR("<NKGFermiMu> must be > 0");
    return eFailure;
  }

  topB.GetChild("NKGFermiTau").GetData(fgNKGFermiTau);
  if (fgNKGFermiTau <= 0) {
    ERROR("<NKGFermiTau> must be > 0");
    return eFailure;
  }

  topB.GetChild("minimizationMethod").GetData(tmp);
  fUseMaxLike = (tmp == "MaxLike");

  topB.GetChild("maxChi2").GetData(fMaxChi2);

  topB.GetChild("useSilentStations").GetData(fUseSilentStations);

  topB.GetChild("silentMaxRadius").GetData(fgSilentMaxRadius);
  if (fgSilentMaxRadius <= 0) {
    ERROR("<silentMaxRadius> must be > 0");
    return eFailure;
  }

  topB.GetChild("silentRadiusTransition").GetData(fgSilentRadiusTransition);
  if (fgSilentRadiusTransition <= 0) {
    ERROR("<silentRadiusTransition> must be > 0");
    return eFailure;
  }

  topB.GetChild("silentSignalThreshold").GetData(fgSilentSignalThreshold);

  topB.GetChild("maxBaryToCoreDistance").GetData(fMaxBaryToCoreDistance);
  if (fMaxBaryToCoreDistance <= 0) {
    ERROR("<maxBaryToCoreDistance> must be > 0");
    return eFailure;
  }

  topB.GetChild("RoptN").GetData(fRoptN);

  const Branch stage0B = topB.GetChild("stage0");
  stage0B.GetChild("minNumberOfStations").GetData(fStage0.fMinNumberOfStations);
  stage0B.GetChild("useCorePosition").GetData(fStage0.fUseCorePosition);

  const Branch stage1B = topB.GetChild("stage1");
  stage1B.GetChild("minNumberOfStations").GetData(fStage1.fMinNumberOfStations);
  stage1B.GetChild("maxChi2").GetData(fStage1.fMaxChi2);

  const Branch stage2B = topB.GetChild("stage2");
  stage2B.GetChild("minNumberRelaxBeta").GetData(fStage2.fMinNumberRelaxBeta);
  stage2B.GetChild("minNumberRelaxGamma").GetData(fStage2.fMinNumberRelaxGamma);

  const Branch stage3B = topB.GetChild("stage3");
  stage3B.GetChild("minNumberOfStations").GetData(fStage3.fMinNumberOfStations);
  stage3B.GetChild("minNumberForFullCurvatureFit").GetData(fStage3.fMinNumberForFullCurvatureFit);
  stage3B.GetChild("maxAxisAngleDifference").GetData(fStage3.fMaxAxisAngleDifference);

  if (fInfoLevel != eNone) {
    ostringstream info;
    info << "LDF type: " << (fgLDFType == ePL ? "PL" : (fgLDFType == eNKG ? "NKG" : "NKGFermi"))
         << ", minimization: " << (fUseMaxLike ? "Maximum-Likelihood" : "Chi2");
    INFO(info);
  }

  const Branch energyBranch = topB.GetChild("energyCalibration");
  energyBranch.GetChild("attenuationPar1").GetData(fEnergyConversion.fAttenuationPar1);
  energyBranch.GetChild("attenuationPar2").GetData(fEnergyConversion.fAttenuationPar2);
  energyBranch.GetChild("attenuationPar3").GetData(fEnergyConversion.fAttenuationPar3);
  energyBranch.GetChild("energyS38Const").GetData(fEnergyConversion.fEnergyS38Const);
  energyBranch.GetChild("energyS38Slope").GetData(fEnergyConversion.fEnergyS38Slope);

  if (fInfoLevel != eNone) {
    ostringstream info;
    info << "S38 conversion: " << "1+(" << fEnergyConversion.fAttenuationPar1 << ")*x+("
         << fEnergyConversion.fAttenuationPar2 << ")*x^2+("
         << fEnergyConversion.fAttenuationPar3 << ")*x^3 ;"
            " energy [eV] = " << fEnergyConversion.fEnergyS38Const << "*"
            " pow(S38," << fEnergyConversion.fEnergyS38Slope << ')';
    INFO(info);
  }

  //the limit for the second derivative of the ldf
  fgNoRecovery = 1;
  return eSuccess;
}


inline
void
LDFFinder::EstimateLDF(LDFFitInterface& fi, const Point& core)
  const
{
  EstimateBetaGamma(fi.fBeta, fi.fGamma, fgLDFType, fgShowerAxis.GetCosTheta(fgLocalCS));
  EstimateS1000(fi);
  OUT(eIntermediate)
    << "  S1000  = " << fi.fS1000 << ", "
       "beta = " << fi.fBeta << ", "
       "gamma = " << fi.fGamma << '\n';
  if (fStage0.fUseCorePosition) {
    // project core to the local tangent plane
    Vector localNormal(0,0,1, fgLocalCS);
    fi.fCore = core -
      (localNormal*(core-fgBarycenter))/(localNormal*fgShowerAxis)*fgShowerAxis;
    OUT(eIntermediate)
      << "  core   = " << ToString(fi.fCore, fgLocalCS) << " (bary)\n";
  } else {
    fi.fCore = fgBarycenter;
    OUT(eIntermediate)
      << "  using barycenter as core\n";
  }
  EstimateEnergy(fi);
  OUT(eIntermediate)
    << "  energy = " << fi.fEnergy/EeV << " [EeV]\n";
}


inline
bool
LDFFinder::FitLDF(LDFFitInterface& fi, const int nStations)
  const
{
  fgUseSilentStations = fi.fUseSilents;
  fi.fLikelihood = FitLDFDriver(fi);
  
  if (fUseMaxLike && fi.fLikelihood)
    fi.fChi2 = EstimateChi2(fi);

  if (fi.fLikelihood < 0) {
    OUT(eIntermediate)
      << "  Fit failed.\n";
    return false;
  } else {

    fi.fNdof = nStations - 1 - (fCoreType == eFit ? 2 : 0) - !fi.fFixedBeta - !fi.fFixedGamma;

    const double nchi2 = fi.fNdof > 0 ? fi.fChi2/fi.fNdof : 0;
    
    if (nchi2 > fi.fMaxChi2) {
      OUT(eIntermediate)
        << "  Chi2/NDOF = " << nchi2 << " exceeds the limit "
        << fi.fMaxChi2 << '\n';
      return false;
    } else {

      const double d = (fi.fCore - fgBarycenter).GetMag();
      if (d > fMaxBaryToCoreDistance) {
        OUT(eIntermediate)
          << "  Distance to core " << d/km << " [km] is above the limit.\n";
        return false;
      } else {
        
        EstimateEnergy(fi);

        OUT(eIntermediate)
          << "  chi2   = " << fi.fChi2 << " / " << fi.fNdof << "\n  "
             "  S1000  = " << fi.fS1000 << " +/- " << fi.fS1000Err << "\n  "
             "  core   = " << ToString(fi.fCore, fgLocalCS) << " +/-\n  "
             "           " << ToString(fi.fCoreErr, fgLocalCS) << " (bary)\n";
        if (!fi.fFixedBeta)
          OUT(eIntermediate)
            << "  beta   = " << fi.fBeta << " +/- " << fi.fBetaErr << '\n';
        if (!fi.fFixedGamma)
          OUT(eIntermediate)
            << "  gamma = " << fi.fGamma << " +/- " << fi.fGammaErr << '\n';
        OUT(eIntermediate)
          << "  energy = "
          << fi.fEnergy/EeV << " +/- "
          << fi.fEnergyErr/EeV << " [EeV]\n";
      }
    }
  }

  return true;
}


VModule::ResultFlag
LDFFinder::Run(evt::Event& event)
{
  INFO("");

  // nothing to do if there is no SEvent
  if (!(event.HasSEvent() && event.HasRecShower() &&
        event.GetRecShower().HasSRecShower()))
    return eContinueLoop;

  const SEvent& sEvent = event.GetSEvent();
  fgCurrentSEvent = &sEvent;

  const ShowerSRecData& showerRec = event.GetRecShower().GetSRecShower();

  int nStations = 0;
  fgLowerLimit = false;
  {
    int nSaturated = 0;
    int nSilent = 0;
    for (SEvent::ConstStationIterator sIt = sEvent.StationsBegin();
         sIt != sEvent.StationsEnd(); ++sIt)
      if (sIt->IsCandidate()) {
        ++nStations;
        if (sIt->IsLowGainSaturation())
          ++nSaturated;
      } else if (sIt->IsSilent())
        ++nSilent;

    OUT(eIntermediate)
      << "# candidate stations = " << nStations
      << " (" << nSaturated << " saturated), " << nSilent << " silent\n";
  }

  // not worth it
  if (nStations < fStage0.fMinNumberOfStations) {
    OUT(eFinal)
      << "Not enough stations.\n";
    return eContinueLoop;
  }

  fgBarycenter = showerRec.GetBarycenter();
  fgBaryTime = showerRec.GetBarytime();
  fgLocalCS = showerRec.GetBarycenterCoordinateSystem();

  if (!fgBaryTime)
    throw utl::NonExistentComponentException("needs barycenter!");

  if (fCoreType != eFit) {

    Point core;
    TimeStamp coreTime;

    if (!FixCore(event, core, coreTime, fgShowerAxis))
      return eContinueLoop;

    ShowerSRecData& showerRec = event.GetRecShower().GetSRecShower();
    showerRec.SetAxis(fgShowerAxis);
    showerRec.SetCorePosition(core);
    showerRec.SetCoreTime(coreTime, TimeInterval(0));
    showerRec.SetAngleErrors(fgShowerAxis.GetCoordinates(fgLocalCS), Vector::Triple(0,0,0));
    showerRec.SetCurvature(0, 0);

  }

  // assume SdPlaneFit module has been run already
  fgShowerAxis = showerRec.GetAxis();
  const double theta = fgShowerAxis.GetTheta(fgLocalCS);

  if (theta > fMaxTheta) {
    OUT(eFinal)
      << "  theta = " << theta/degree << " is > than limit "
      << fMaxTheta/degree << " [deg], skipping...\n";
    return eContinueLoop;
  }

  const Point& core = showerRec.GetCorePosition();

  const CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
  const TimeStamp eventTime(sEvent.GetTrigger().GetTime());

  OUT(eIntermediate)
    << "barycenter = " << ToString(fgBarycenter, siteCS) << " (site)\n  "
       "bary time  = " << (fgBaryTime - eventTime).GetInterval()/nanosecond
    << " [ns] (to event time)\n  "
       "axis       = " << ToString(fgShowerAxis, siteCS, false) << " (site)\n  "
       "           = " << ToString(fgShowerAxis, fgLocalCS, false) << " (bary)\n  "
       "theta      = " << theta/degree << " [deg] (bary)\n  "
       "phi        = " << fgShowerAxis.GetPhi(fgLocalCS)/degree << " [deg] (bary)\n";

  OUT(eObscure)
    << "local/site zenith angle diff. = "
    << acos(Vector(0,0,1, fgLocalCS)*Vector(0,0,1, siteCS))/degree
    << " [deg]\n";

  LDFFitInterface best;
  CurvatureFitInterface cbest;
  cbest.fAxis = fgShowerAxis;

  // --- Stage: LDF estimation ------------------------------------------------

  best.fRecStage = ShowerSRecData::kLDFEstimate;
  OUT(eIntermediate)
    << "Stage " << best.fRecStage << ": estimate beta, gamma, S1000\n";
  EstimateLDF(best, core);

  if (nStations < fStage1.fMinNumberOfStations) {
    SetRecData(event, best, cbest);
    OutputResults(event);
    return eSuccess;
  }

  bool useSilentStations = (nStations > 4);
  // loop: without/with silents
  do {
    useSilentStations = !useSilentStations;
    const bool failOnError = !fUseSilentStations || useSilentStations;
    const string stationType = useSilentStations ? "candidate+silents" : "candidates";

    // --- Stage: S1000 & core fitting ----------------------------------------
    {
      LDFFitInterface fi = best;
      fi.fRecStage = (useSilentStations ? ShowerSRecData::kLDFSilents : ShowerSRecData::kLDF);

      fi.fFixedBeta = true;
      fi.fFixedGamma = true;
      fi.fUseSilents = useSilentStations;
      fi.fMaxChi2 = fStage1.fMaxChi2;
      fgFixLDF = (fi.fFixedBeta && fi.fFixedGamma);

      OUT(eIntermediate)
        << "Stage " << fi.fRecStage << ": fit " << stationType << " for S1000, core\n";
      fgLowerLimit = false;
      if (FitLDF(fi, nStations))
        best = fi;
      else if (failOnError)
        break;
      
    } // stage

    // --- Stage: S1000, core, beta (gamma) fitting ---------------------------
    {
      LDFFitInterface fi = best;

      fi.fFixedBeta = (nStations < fStage2.fMinNumberRelaxBeta || FixBeta(fi)) ? true : false;
      fi.fFixedGamma = (nStations < fStage2.fMinNumberRelaxGamma || FixGamma(fi)) ? true : false;
      fgFixLDF = (fi.fFixedBeta && fi.fFixedGamma);
      if (!fi.fFixedBeta || !fi.fFixedGamma) {

        fi.fRecStage = (fgUseSilentStations ? ShowerSRecData::kLDFSilents : ShowerSRecData::kLDF) +
          (!fi.fFixedBeta) * ShowerSRecData::kLDFBeta +
          (!fi.fFixedGamma) * ShowerSRecData::kLDFGamma;
        fi.fUseSilents = useSilentStations;
        fi.fMaxChi2 = fMaxChi2;

        OUT(eIntermediate)
          << "Stage " << fi.fRecStage << ": fit " << stationType << " for S1000, core"
          << (!fi.fFixedBeta ? ", beta" : "")
          << (!fi.fFixedGamma ? ", gamma\n" : "\n");
        fgLowerLimit = false;
        if (FitLDF(fi, nStations)) {
          const bool fixBeta = FixBeta(fi);
          const bool fixGamma = FixGamma(fi);
          if (fi.fFixedBeta == fixBeta && fi.fFixedGamma == fixGamma)
            best = fi;
          else {
            fi.fFixedBeta = fi.fFixedBeta ? fi.fFixedBeta : fixBeta;
            fi.fFixedGamma = fi.fFixedGamma ? fi.fFixedGamma : fixGamma;
            fgFixLDF = (fi.fFixedBeta && fi.fFixedGamma);
            fi.fRecStage = (fgUseSilentStations ? ShowerSRecData::kLDFSilents : ShowerSRecData::kLDF) +
                           (!fi.fFixedBeta) * ShowerSRecData::kLDFBeta +
                           (!fi.fFixedGamma) * ShowerSRecData::kLDFGamma;
            fgLowerLimit = false;
            if (FitLDF(fi, nStations))
              best = fi;
            else if (failOnError)
              break;
            
          }
        } else if (failOnError)
          break;
        
      }

    } // stage

  } while (fUseSilentStations && !useSilentStations);

  if (nStations < fStage3.fMinNumberOfStations) {
    SetRecData(event, best, cbest);
    OutputResults(event);
    return eSuccess;
  }

  {
    CurvatureFitInterface cfi;
    cfi.fCore = best.fCore;
    cfi.fAxis = cbest.fAxis;
    const double parameterizedRc = ParameterizedRc(best);
    cfi.fRc = parameterizedRc;
    cfi.fct0 = -cfi.fRc;
    double chi2 = 0;

    if (nStations < fStage3.fMinNumberForFullCurvatureFit) {
      cfi.fFixedRc = true;
      OUT(eIntermediate)
        << "Fit for axis with fixed Rc = " << cfi.fRc/kilometer << " [km]\n";
      chi2 = FitCurvatureDriver(cfi);
    } else {
      cfi.fFixedRc = false;
      OUT(eIntermediate)
        << "Fit for shower-front curvature and axis.\n";
      chi2 = FitCurvatureDriver(cfi);
      if (chi2 <= 0) {
        OUT(eIntermediate)
          << "  Fit failed, retrying with fixed Rc = " << parameterizedRc << " [km]\n";
        cfi.fFixedRc = true;
        cfi.fRc = parameterizedRc;
        cfi.fRcErr = 0;
        cfi.fct0 = -cfi.fRc;
        chi2 = FitCurvatureDriver(cfi);
      }
    }

    if (chi2 <= 0) {
      OUT(eIntermediate)
        << "  Last fit failed.\n";
      cbest.fRc = 0;
    } else {
      const double diff = acos(cfi.fAxis * fgShowerAxis);
      OUT(eIntermediate)
        << "  axis = " << ToString(cfi.fAxis, fgLocalCS) << " (bary)\n  "
           "  axis diff = " << diff/degree << " [deg] (to plain-front fit)\n  "
           "  Rc   = " << cfi.fRc/meter << " +/- " << cfi.fRcErr/meter << " [m]\n  "
           "  ct0  = " << cfi.fct0/meter << " +/- " << cfi.fct0Err/meter << " [m]\n";
      if (diff > fStage3.fMaxAxisAngleDifference) {
        OUT(eIntermediate)
          << "  Opening-angle axis difference is above the limit.\n";
        cbest.fRc = 0;
      } else
        cbest = cfi;
    }
  }

  const double newTheta = cbest.fAxis.GetTheta(fgLocalCS);
  if (newTheta > fMaxTheta) {
    OUT(eFinal)
      << "  theta = " << newTheta/degree << " is > than limit "
      << fMaxTheta/degree << " [deg], skipping...\n";
    return eContinueLoop;
  }

  // redo LDF fit after the axis change
  if (cbest.fRc) {
    LDFFitInterface fi = best;
    const Vector oldAxis = fgShowerAxis;
    fgShowerAxis = cbest.fAxis;

    if (best.fRecStage == ShowerSRecData::kLDFEstimate) {
      OUT(eIntermediate)
        << "Redo Stage " << best.fRecStage << ": estimate beta, gamma, S1000\n";
      EstimateLDF(best, core);
      best.fRecStage = ShowerSRecData::kLDFEstimateAndCurvature;
    } else {
      const string stationType = fi.fUseSilents ? "candidate+silents" : "candidates";
      OUT(eIntermediate)
        << "Redo Stage " << fi.fRecStage << ": fit " << stationType << " for S1000, core"
        << (!fi.fFixedBeta ? ", beta" : "")
        << (!fi.fFixedGamma ? ", gamma\n" : "\n");

      EstimateBetaGamma(fi.fBeta, fi.fGamma,
                        fgLDFType, fgShowerAxis.GetCosTheta(fgLocalCS), fi.fS1000);
      fgLowerLimit = false;
      if (FitLDF(fi, nStations)) {
        // showerCenter = oldCore + oldAxis * oldRc
        const Point showerCenter = best.fCore + cbest.fAxis * cbest.fRc;
        // newAxis = normalized(showerCenter - newCore)
        const Vector coreCenter = showerCenter - fi.fCore;
        // newRc = |showerCenter - newCore|
        const double newRc = coreCenter.GetMag();
        cbest.fRc = newRc;
        cbest.fAxis = coreCenter / newRc;
        best = fi;
        best.fRecStage += 0.5;
      } else {
        // return to the last successful results if the LDF fit fails
        fgShowerAxis = oldAxis;
        cbest.fRc = 0;
      }
    }
  }

  SetRecData(event, best, cbest);

  // calculation of Ropt, see GAP-2006-045
  const double finalBeta = best.fBeta;
  const double betaFixSigma = SigmaForFixedBeta(best.fS1000);
  //  const int fRoptN = 2;
  double vS1000[fRoptN];
  double vBeta[fRoptN];
  const double startp = NormalCDF(-1);
  const double stopp = NormalCDF(1);
  const double deltap = (stopp - startp) / (fRoptN - 1);
  OUT(eIntermediate)
    << "Ropt determination: \n"; 
  for (int i = 0; i < fRoptN; ++i) {
    LDFFitInterface fi = best;
    const double p = startp + i*deltap;
    fi.fBeta = vBeta[i] = finalBeta + InverseNormalCDF(p, betaFixSigma);
    fi.fFixedBeta = true;
    fgFixLDF = false;
    fi.fGamma = best.fGamma;
    fi.fFixedGamma = true;
    FitLDF(fi, nStations);
    vS1000[i] = fi.fS1000;
  }

  RoptFit roptFit;
  pair<double, double> ropt = roptFit.Fit(fgLDFType, fRoptN, vS1000, vBeta, fMinuitOutput);
  ShowerSRecData& shRec = event.GetRecShower().GetSRecShower();
  shRec.SetLDFOptimumDistance(ropt.first);
  // "derivative" (finite difference)
  const double dS1000_dBeta = (vS1000[0] - vS1000[fRoptN-1]) / (vBeta[0] - vBeta[fRoptN-1]);
  // beta spread parameterization from Talianna
  shRec.SetBetaSystematics(betaFixSigma);
  shRec.SetShowerSizeSystematics(fabs(betaFixSigma * dS1000_dBeta));

  OutputResults(event);

  return eSuccess;
}


bool
LDFFinder::FixCore(const evt::Event& event, Point& core,
                   TimeStamp& refTime, Vector& axis)
  const
{
  const Vector localNormal(0,0,1, fgLocalCS);
  Point coreReference;


  switch (fCoreType) {
  case eMC: {
    if (!event.HasSimShower()) {
      ERROR("The event is not a simulated shower.");
      return false;
    }
    const ShowerSimData& simShower = event.GetSimShower();
    axis = -simShower.GetDirection();
    coreReference = simShower.GetPosition();
    refTime = simShower.GetTimeStamp();
    break;
  }
  case eHybrid:
    if (!event.HasFEvent()) {
      ERROR("The event has no FD shower.");
      return false;
    }

    const fevt::FEvent& fEvent = event.GetFEvent();
    const bool stereoRec =
      event.HasRecShower() && event.GetRecShower().HasFRecShower();
    // remember area of FD core error ellipse to choose best FD eye
    double bestErrorEllipse = 0;

    const evt::ShowerFRecData* showerFRec = 0;
    if (stereoRec)
      showerFRec = &event.GetRecShower().GetFRecShower();

    for (FEvent::ConstEyeIterator eyeIt = fEvent.EyesBegin(ComponentSelector::eHasData);
         eyeIt != fEvent.EyesEnd(ComponentSelector::eHasData); ++eyeIt) {

      // skip eyes without reconstruction
      if (!eyeIt->HasRecData())
        continue;

      const fevt::EyeRecData& eyeRec = eyeIt->GetRecData();
      if (eyeRec.GetSDPFitNDof() <= 0 || eyeRec.GetTimeFitNDof() <= 0)
        continue;

      if (!eyeRec.HasFRecShower())
        continue;

      if (!stereoRec)
        showerFRec = &eyeRec.GetFRecShower();

      const fdet::Eye& dEye = det::Detector::GetInstance().GetFDetector().GetEye(*eyeIt);
      const CoordinateSystemPtr eyeCS = dEye.GetEyeCoordinateSystem();
      const Point eyePos = dEye.GetPosition();
      const Vector vecEyeCore = eyePos - showerFRec->GetCorePosition();

      const double rp = eyeRec.GetRp();
      const double rpErr = eyeRec.GetRpError();
      const double t0 = eyeRec.GetTZero();
      const double chi0 = eyeRec.GetChiZero();
      const double chi0Err = eyeRec.GetChiZeroError();
      const double sdpPhiErr = eyeRec.GetSDPPhiError();
      const double rhoChi0Rp = eyeRec.GetRpChi0Correlation();
      const double distEyeCore = vecEyeCore.GetMag();

      // calculation of core error in SDP
      // rp error projected on ground+
      // chi0 error+ correlation of errors in rp and chi0
      const double sinChi0 = sin(chi0);
      const double errorInSDP = sqrt(Sqr(rpErr/sinChi0) +
                                     Sqr((chi0Err*distEyeCore)/(cos(chi0)*sinChi0)) +
                                     rhoChi0Rp * rpErr * chi0Err);

      // error due to SDP azimuth uncertainty
      const double errorPerpSDP = distEyeCore * sin(sdpPhiErr);

      // area of core error ellipse of this FD event
      const double coreErrorEllipse = errorInSDP * errorPerpSDP * kPi;

      // if error ellipse set AND error ellipse larger than previous eye(s)
      // don't update core with actual eye
      if (bestErrorEllipse && bestErrorEllipse < coreErrorEllipse)
        continue;

      bestErrorEllipse = coreErrorEllipse;
      axis = showerFRec->GetAxis();
      coreReference = showerFRec->GetCorePosition();

      const Vector vertical(0,0,1, eyeCS);
      const double alpha = 0.5*kPi - Angle(vertical, vecEyeCore);
      const double cosBeta = CosAngle(axis, vecEyeCore);
      // 1e5ns have to be subtracted to jump from the end to the
      // start of the pixel trace (don't even think of removing the
      // comment here)
      refTime = eyeIt->GetHeader().GetTimeStamp() +
        TimeInterval(-1e5*ns + t0 + (rp*tan(alpha) + distEyeCore*cosBeta)/kSpeedOfLight);
      break;
    }

    if (!showerFRec)
      return false;
    break;
  }

  core = coreReference -
    (localNormal * (coreReference - fgBarycenter)) / (localNormal*axis) * axis;
 
  const double d = (core - fgBarycenter).GetMag();
  if (d > fMaxBaryToCoreDistance) {
    ostringstream err;
    err << "Bary-centric distance to core " << d/kilometer
        << " [km] is above the limit.";
    ERROR(err);
    return false;
  }

  return true;
}


// From Pierre Billoir Thoughts....
// - if there are enough stations (especially at distances around 1000 m), a fit
// with a variable LDF slope is attempted (B is set as adjustable in (3)); the
// conditions are:
// * at least 5 stations, with one of the following requirements:
// * at least 2 within (500, 1500), with a difference at least 500 in between
// * at least 3 within (500, 1500), with a maximal difference at least 400
// * at least 4 within (500, 1500), with a maximal difference at least 300
bool
LDFFinder::FixBeta(const LDFFitInterface& fi)
  const
{
  vector<int> candidates;
  const SEvent::ConstCandidateStationIterator sEnd = fgCurrentSEvent->CandidateStationsEnd();
  for (SEvent::ConstCandidateStationIterator sIt = fgCurrentSEvent->CandidateStationsBegin();
       sIt != sEnd; ++sIt)
    if (!sIt->IsLowGainSaturation() || sIt->GetRecData().GetRecoveredSignal())
      candidates.push_back(sIt->GetId());

  const int nStations = candidates.size();
  if (nStations < 5) {
    OUT(eIntermediate)
      << "FixBeta = " << true
      << ": nStations = " << nStations << '\n';
    return true;
  }

  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
  int nStationsInRange = 0;
  double maxDistance = 0;

  for (vector<int>::const_iterator sIt = candidates.begin();
       sIt != candidates.end(); ++sIt) {

    const double r =
      RPerp(fgShowerAxis, sDetector.GetStation(*sIt).GetPosition() - fi.fCore);

    if (500*meter < r && r < 1500*meter) {

      ++nStationsInRange;

      for (vector<int>::const_iterator sJt = sIt + 1;
           sJt != candidates.end(); ++sJt) {

        const double rr =
          RPerp(fgShowerAxis, sDetector.GetStation(*sJt).GetPosition() - fi.fCore);

        if (500*meter < rr && rr < 1500*meter) {
          const double distance = fabs(r - rr);
          if (distance > maxDistance)
            maxDistance = distance;
        }
      }
    }
  }

  const bool fixBeta = !((nStationsInRange > 1 && maxDistance > 500*meter) ||
                         (nStationsInRange > 2 && maxDistance > 400*meter) ||
                         (nStationsInRange > 3 && maxDistance > 300*meter));
  OUT(eIntermediate)
    << "FixBeta = " << fixBeta
    << ": nStations (InRange) = " << nStations << " (" << nStationsInRange
    << "), maxDistance = " << maxDistance/meter << " [m]\n";

  return fixBeta;
}


// Criterion if gamma can be fitted analogue to FixBeta.
bool
LDFFinder::FixGamma(const LDFFitInterface& fi)
  const
{
  vector<int> candidates;
  const SEvent::ConstCandidateStationIterator sEnd = fgCurrentSEvent->CandidateStationsEnd();
  for (SEvent::ConstCandidateStationIterator sIt = fgCurrentSEvent->CandidateStationsBegin();
       sIt != sEnd; ++sIt)
    if (!sIt->IsLowGainSaturation() || sIt->GetRecData().GetRecoveredSignal())
      candidates.push_back(sIt->GetId());

  const int nStations = candidates.size();
  if (nStations < 5) {
    OUT(eIntermediate)
      << "FixGamma = " << true
      << ": nStations = " << nStations << '\n';
    return true;
  }

  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
  int nStationsInRange = 0;
  double maxDistance = 0;

  for (vector<int>::const_iterator sIt = candidates.begin();
       sIt != candidates.end(); ++sIt) {

    const double r =
      RPerp(fgShowerAxis, sDetector.GetStation(*sIt).GetPosition() - fi.fCore);

    if (1000*meter < r && r < 2000*meter) {

      ++nStationsInRange;

      for (vector<int>::const_iterator sJt = sIt + 1;
           sJt != candidates.end(); ++sJt) {

        const double rr =
          RPerp(fgShowerAxis, sDetector.GetStation(*sJt).GetPosition() - fi.fCore);

        if (1000*meter < rr && rr < 2000*meter) {
          const double distance = fabs(r - rr);
          if (distance > maxDistance)
            maxDistance = distance;
        }
      }
    }
  }

  const bool fixGamma = !((nStationsInRange > 1 && maxDistance > 500*meter) ||
                         (nStationsInRange > 2 && maxDistance > 400*meter) ||
                         (nStationsInRange > 3 && maxDistance > 300*meter));
  OUT(eIntermediate)
    << "FixGamma = " << fixGamma
    << ": nStations (InRange) = " << nStations << " (" << nStationsInRange
    << "), maxDistance = " << maxDistance/meter << " [m]\n";

  return fixGamma;
}


void
LDFFinder::EstimateS1000(LDFFitInterface& fi)
  const
{
  const double k1000m = 1000*meter;
  double rMin = 1e20;
  double rClosest = k1000m;
  double signal = 1;

  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

  // loop on stations
  for (SEvent::ConstCandidateStationIterator sIt = fgCurrentSEvent->CandidateStationsBegin();
       sIt != fgCurrentSEvent->CandidateStationsEnd(); ++sIt) {

    const double r =
      RPerp(fgShowerAxis, sDetector.GetStation(*sIt).GetPosition() - fgBarycenter);

    if (abs(r - k1000m) < rMin) {
      rMin = fabs(r - k1000m);
      signal = sIt->GetRecData().GetTotalSignal();
      rClosest = r;
    }

  }

  fi.fS1000 = signal / LDFFunction(fgLDFType, rClosest, fi.fBeta, fi.fGamma, fgNKGFermiMu, fgNKGFermiTau);
}


double
LDFFinder::EstimateNStationsInFit(const LDFFitInterface& fi)
  const
{
  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

  double nStations = 0;
  // loop on stations
  for (SEvent::ConstStationIterator sIt = fgCurrentSEvent->StationsBegin();
       sIt != fgCurrentSEvent->StationsEnd(); ++sIt) {

    if (sIt->IsCandidate())
      ++nStations;
    else if (fgUseSilentStations && sIt->IsSilent()) {
      const double r =
        RPerp(fgShowerAxis, sDetector.GetStation(*sIt).GetPosition() - fi.fCore);
      nStations +=
        Fermi(r, LDFFinder::fgSilentMaxRadius, LDFFinder::fgSilentRadiusTransition);
    }

  }

  return nStations;
}


void
LDFFinder::EstimateEnergy(LDFFitInterface& fi)
  const
{
  fEnergyConversion.GetEnergy(fgShowerAxis.GetCosTheta(fgLocalCS), fi.fS1000, fi.fS1000Err,
                              fi.fEnergy, fi.fEnergyErr);
}


double
LDFFinder::ParameterizedRc(const LDFFitInterface& fi)
  const
{
  const double s = log10(fi.fS1000);
  const double s2 = Sqr(s);
  /* markus' fit
  const double offset = 9.03 - 3.08*s + 1.08*s2;
  const double slope = 7.28 + 4.03*s - 1.01*s2;
  */
  const double offset = 7.418 - 3.0613*s + 2.5262*s2;
  const double slope = 3.9779 + 3.2245*s + 4.5705*s2;
  const double sec1 = 1/fgShowerAxis.GetCosTheta(fgLocalCS) - 1;
  return (offset + slope*sec1) * kilometer;
}


void
LDFFinder::SetRecData(Event& event, const LDFFitInterface& lfi,
                      const CurvatureFitInterface& cfi)
  const
{
  ShowerSRecData& recShower =
    event.GetRecShower().GetSRecShower();

  const sdet::SDetector& sDetector =
    det::Detector::GetInstance().GetSDetector();

  const double s1000 = lfi.fS1000;
  const double beta = lfi.fBeta;
  const double gamma = lfi.fGamma;
  const double cosTheta = fgShowerAxis.GetCosTheta(fgLocalCS);
  const double sigmaFactor = wcd::SignalUncertaintyFactor(wcd::eGAP2012_012, cosTheta);

  SEvent& sEvent = event.GetSEvent();

  const bool haveCurvature = (cfi.fRc != 0);
  const Point showerOrigin(lfi.fCore + cfi.fRc * cfi.fAxis);
  const double ct0 = cfi.fct0;

  const double theta = fgShowerAxis.GetTheta(fgLocalCS);
  const double phi = fgShowerAxis.GetPhi(fgLocalCS);
  const CoordinateSystemPtr showerPlaneCS(
    CoordinateSystem::RotationY(
      theta,
      CoordinateSystem::RotationZ(
        phi,
        LocalCoordinateSystem::Create(lfi.fCore)
      )
    )
  );

  HepMatrix hepErrMatrix(3,3,0);
  if (fCoreType == eFit) {
    hepErrMatrix[0][0] = lfi.fErrMatrix[1][1];
    hepErrMatrix[0][1] = lfi.fErrMatrix[1][2];
    hepErrMatrix[1][0] = lfi.fErrMatrix[2][1];
    hepErrMatrix[1][1] = lfi.fErrMatrix[2][2];
  } else
    for (int i = 0; i < 2; ++i)
      for (int j = 0; j < 2; ++j)
        hepErrMatrix[i][j] = (i == j);

  HepRotation r;
  r.rotateZ(-phi);
  r.rotateY(-theta);

  HepMatrix hepRotationMatrix(3,3);
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      hepRotationMatrix[i][j] = r[i][j];

  HepMatrix hepRotatedErrMatrix(3,3);
  hepRotatedErrMatrix = hepRotationMatrix * hepErrMatrix *  hepRotationMatrix.T();

  const sdet::STimeVariance& timeVariance = sdet::STimeVariance::GetInstance();
  Accumulator::MinMax<double> minMaxRho(1000*kilometer, 0);
  Accumulator::SampleStandardDeviation timeStat;

  for (SEvent::StationIterator sIt = sEvent.StationsBegin();
       sIt != sEvent.StationsEnd(); ++sIt)

    if (sIt->HasRecData()) {

      const Point& sPos = sDetector.GetStation(*sIt).GetPosition();
      const double rho2 = RPerp2(cfi.fAxis, sPos - lfi.fCore);
      const double rho = sqrt(rho2);
      minMaxRho(rho);

      StationRecData& sRec = sIt->GetRecData();
      const double signal = sRec.GetTotalSignal();
      const double funcval = s1000 * LDFFunction(fgLDFType, rho, beta, gamma, fgNKGFermiMu, fgNKGFermiTau);
      const double sigma = sigmaFactor * sqrt(funcval);

      const double drho2 =
        (Sqr(sPos.GetX(showerPlaneCS)) * hepRotatedErrMatrix[0][0] +
         Sqr(sPos.GetY(showerPlaneCS)) * hepRotatedErrMatrix[1][1] +
         2 * sPos.GetX(showerPlaneCS) * sPos.GetY(showerPlaneCS) * hepRotatedErrMatrix[0][1]) / rho2;

      const double drho = sqrt(drho2);

      sRec.SetAzimuthShowerPlane(sPos.GetPhi(showerPlaneCS));
      sRec.SetTotalSignal(signal,sigma);
      sRec.SetLDFResidual((signal - funcval)/sigma);
      sRec.SetSPDistance(rho, drho);

      if (haveCurvature) {
        const double measuredTime = (sRec.GetSignalStartTime() - fgBaryTime).GetInterval();
        const Vector stationOrigin = showerOrigin - sPos;
        const double distOrigin = stationOrigin.GetMag();
        const double predictedTime = (distOrigin + ct0) / kSpeedOfLight;
        const double stationCosTheta = stationOrigin.GetZ(fgLocalCS) / distOrigin;
        const double sigma =
          sqrt(timeVariance.GetTimeSigma2(signal, sRec.GetT50(), stationCosTheta, rho));
        const double residual = (measuredTime - predictedTime)/sigma;
        sRec.SetResidual(residual);
        if (sIt->IsCandidate())
          timeStat(residual);
      }

    }

  double minRho = 0.5 * minMaxRho.GetMin();
  double maxRho = 1.2 * minMaxRho.GetMax();
  if (maxRho < 1200*meter)
    maxRho = 1200*meter;

  if (!recShower.HasLDF())
    recShower.MakeLDF();
  TabulatedFunctionErrors& ldf = recShower.GetLDF();
  ldf.Clear();

  const int nLDFPoints = 200;
  const double dRho = (maxRho - minRho)/(nLDFPoints - 1);
  for (int i = 0; i < nLDFPoints; ++i) {
    const double rho = minRho + i*dRho;
    const double funcval = s1000 * LDFFunction(fgLDFType, rho, beta, gamma, fgNKGFermiMu, fgNKGFermiTau);
    const double sigma = sigmaFactor * sqrt(funcval);
    ldf.PushBack(rho, 0, funcval, sigma);
  }

  // set the reconstruction data
  if (haveCurvature) {
    recShower.SetAxis(cfi.fAxis);
    const double rc = cfi.fRc;
    recShower.SetCurvature(1/rc, cfi.fRcErr/Sqr(rc));
    recShower.SetCoreTime(fgBaryTime + TimeInterval((cfi.fct0 + rc)/kSpeedOfLight),
                          TimeInterval(sqrt(Sqr(cfi.fct0Err) + Sqr(cfi.fRcErr))/kSpeedOfLight));
    const double timeMean = timeStat.GetAverage();
    const double timeSpread = timeStat.GetStandardDeviation();
    recShower.SetTimeResidualMean(timeMean);
    recShower.SetTimeResidualSpread(timeSpread);
    const int ndof = timeStat.GetN()  // number of candidate stations
      - (cfi.fFixedAxis ? 0 : 2)      // for u and v
      - (cfi.fFixedRc ? 0 : 1)        // for Rc
      - 1;                            // for t0
    recShower.SetAngleChi2(timeStat.GetSumOfSquares(), ndof);
    recShower.SetAngleErrors(cfi.fAxis.GetCoordinates(fgLocalCS),
                             Vector::Triple(cfi.fSigmaU2, cfi.fSigmaUV, cfi.fSigmaV2));
  } else {
    // we have moved the core but have axis from plane fit, recalc t0
    const TimeInterval cdiff =
      recShower.GetAxis() * (recShower.GetCorePosition() - lfi.fCore);
    const TimeStamp newCoreTime =
      recShower.GetCoreTime() + cdiff/kSpeedOfLight;
    recShower.SetCoreTime(newCoreTime, recShower.GetCoreTimeError());
    /*OUT(eIntermediate)
      << "Core moved, core time change: " << cdiff/meter << " [m].\n";*/
  }
  recShower.SetNumberOfActiveStations(timeStat.GetN());
  recShower.SetShowerSizeType(ShowerSRecData::eS1000);
  recShower.SetShowerSize(lfi.fS1000, lfi.fS1000Err);
  recShower.SetCorePosition(lfi.fCore);
  recShower.SetCoreError(lfi.fCoreErr);
  recShower.SetCorrelationXY(hepErrMatrix[0][1]/sqrt(hepErrMatrix[0][0]*hepErrMatrix[1][1]));
  recShower.SetBeta(lfi.fBeta, lfi.fBetaErr);
  recShower.SetGamma(lfi.fGamma, lfi.fGammaErr);
  recShower.SetNKGFermiParameters(fgNKGFermiMu, fgNKGFermiTau);
  recShower.SetLDFChi2(lfi.fChi2, lfi.fNdof);
  recShower.SetLDFLikelihood(lfi.fLikelihood);
  double energy = 0;
  double energyErr = 0;
  
  fEnergyConversion.GetEnergy(fgShowerAxis.GetCosTheta(fgLocalCS), lfi.fS1000, lfi.fS1000Err,
                              energy, energyErr);

  recShower.SetEnergy(energy, energyErr);
  recShower.SetLDFRecStage(lfi.fRecStage + fgLowerLimit*ShowerSRecData::kLDFLowerLimit);
}


VModule::ResultFlag
LDFFinder::Finish()
{
  return eSuccess;
}


double
LDFFinder::FitLDFDriver(LDFFitInterface& fi)
  const
{
  TMinuit minuit(5);
  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 (fUseMaxLike) {
    minuit.SetFCN(LDFFinder::LDFFitMaxLikeFnc);
    argList[0] = 0.5;
  } else {
    minuit.SetFCN(LDFFinder::LDFFitChi2Fnc);
    argList[0] = 1;
  }
  minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);

  minuit.mnparm(0, "s1000", fi.fS1000, 1., 0, 0, errFlag);
  minuit.mnparm(1, "x", fi.fCore.GetX(fgLocalCS), 200*meter, 0, 0, errFlag);
  minuit.mnparm(2, "y", fi.fCore.GetY(fgLocalCS), 200*meter, 0, 0, errFlag);
  minuit.mnparm(3, "beta", fi.fBeta, 0.1, 0, 0, errFlag);
  minuit.mnparm(4, "gamma", fi.fGamma, 0.02, 0, 0, errFlag);

  if (fCoreType != eFit) {
    minuit.FixParameter(1);
    minuit.FixParameter(2);
  }

  if (fi.fFixedBeta)
    minuit.FixParameter(3);
  if (fi.fFixedGamma)
    minuit.FixParameter(4);

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

  argList[0] = 500;
  minuit.mnexcm("MINIMIZE", argList, 1, errFlag);
  if (errFlag) {
    OUT(eObscure)
      << "  minuit minimize error: " << errFlag << '\n';
    minuit.mnexcm("MINOS", argList, 0, errFlag);
    if (errFlag) {
      OUT(eObscure)
        << "  minuit minos error: " << errFlag << '\n';
      return -1;
    }
  }

  // Get results
  minuit.GetParameter(0, fi.fS1000, fi.fS1000Err);
  double x;
  double y;
  double xErr;
  double yErr;
  if (fCoreType == eFit) {
    minuit.GetParameter(1, x, xErr);
    minuit.GetParameter(2, y, yErr);
    fi.fCore = Point(x, y, 0, fgLocalCS);
    fi.fCoreErr = Vector(xErr, yErr, 0, fgLocalCS);
  }
  if (!fi.fFixedBeta)
    minuit.GetParameter(3, fi.fBeta, fi.fBetaErr);
  else {
    fi.fBetaErr = 0;
    EstimateBetaGamma(fi.fBeta, fi.fGamma,
                      fgLDFType, fgShowerAxis.GetCosTheta(fgLocalCS), fi.fS1000);
  }
  if (!fi.fFixedGamma)
    minuit.GetParameter(4, fi.fGamma, fi.fGammaErr);
  else
    fi.fGammaErr = 0;

  // error matrix
  if (minuit.fNpar > 5) {
    ostringstream err;
    err << "Number of parameters used in Minuit covarince matrix larger "
           "than size of matrix fErrMatrix";
    ERROR(err);

    throw utl::OutOfBoundException(err.str());
  }

  minuit.mnemat(&fi.fErrMatrix[0][0], 5);

  return minuit.fAmin;
}


double
LDFFinder::EstimateChi2(const LDFFitInterface& fi)
  const
{
  TMinuit minuit(5);
  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);
  }

  // disable silent stations for chi2 calculation
  const bool saveSilent = LDFFinder::fgUseSilentStations;
  LDFFinder::fgUseSilentStations = false;

  minuit.SetFCN(LDFFinder::LDFFitChi2Fnc);
  argList[0] = 1;
  minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);

  minuit.mnparm(0, "s1000", fi.fS1000, 0, 0, 0, errFlag);
  minuit.mnparm(1, "x", fi.fCore.GetX(fgLocalCS), 0, 0, 0, errFlag);
  minuit.mnparm(2, "y", fi.fCore.GetY(fgLocalCS), 0, 0, 0, errFlag);
  minuit.mnparm(3, "beta", fi.fBeta, 0, 0, 0, errFlag);
  minuit.mnparm(4, "gamma", fi.fGamma, 0, 0, 0, errFlag);

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

  // restore
  LDFFinder::fgUseSilentStations = saveSilent;

  if (errFlag) {
    OUT(eObscure)
      << "  minuit get chi2 error: " << errFlag << '\n';
    return 0;
  }

  return minuit.fAmin;
}


void
LDFFinder::LDFFitMaxLikeFnc(int& /*nPar*/, double* const /*grad*/, double& value,
                            double* const par, const int iFlag)
{
  // One comment on the following conversion factor called PoissonFactor:
  // The factor reflects the electron-gamma to muon ratio and gives rise to
  // criticism of the algorithm. But as long as we would like to include
  // "active" Zero stations, we have to assume a conversion factor. Otherwise we
  // cannot employ a maximum likelihood method! By definion a \chi^2 method gives
  // biased results (more strongly stated: in a mathematical sense they are
  // wrong) !
  // This factor has to be studied in detail by means of simulated AND real
  // data. See Haverah Park analysis for more details (Lapikens? or England?)

  const double& x = par[1];
  const double& y = par[2];
  const double s1000 = par[0];
  double beta = 0;
  double gamma = 0;

  if (s1000 <= 0) {
    value = numeric_limits<double>::max();
    return;
  }

  if (LDFFinder::fgFixLDF)
    EstimateBetaGamma(beta, gamma, fgLDFType, fgShowerAxis.GetCosTheta(fgLocalCS), s1000);
  else  {
    beta = par[3];
    gamma = par[4];
  }

  static vector<bool> sIsCandidate;
  static vector<Point> sPosition;
  static vector<double> sSignal;
  static vector<double> sSignalError;
  static vector<bool> sIsSaturation;

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

    sIsCandidate.clear();
    sPosition.clear();
    sSignal.clear();
    sSignalError.clear();
    sIsSaturation.clear();

    const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

    for (SEvent::ConstStationIterator sIt = fgCurrentSEvent->StationsBegin();
         sIt != fgCurrentSEvent->StationsEnd(); ++sIt)

      if (sIt->IsCandidate()) {

        sIsCandidate.push_back(true);
        sPosition.push_back(sDetector.GetStation(*sIt).GetPosition());
        const bool saturated = sIt->IsLowGainSaturation();
        if (saturated && sIt->GetRecData().GetRecoveredSignal()) {
          sSignal.push_back(sIt->GetRecData().GetRecoveredSignal());
          sSignalError.push_back(sIt->GetRecData().GetRecoveredSignalError());
        } else {
          sSignal.push_back(sIt->GetRecData().GetTotalSignal());
          sSignalError.push_back(0);
        }
        sIsSaturation.push_back(saturated);

      } else if (LDFFinder::fgUseSilentStations && sIt->IsSilent()) {

        sIsCandidate.push_back(false);
        sPosition.push_back(sDetector.GetStation(*sIt).GetPosition());
        sSignal.push_back(0);
        sSignalError.push_back(0);
        sIsSaturation.push_back(false);

      }

  } // fnc init

  const int nStations = sIsCandidate.size();

  const double cosTheta = fgShowerAxis.GetCosTheta(fgLocalCS);
  const double sigmaFactor = wcd::SignalUncertaintyFactor(wcd::eGAP2012_012, cosTheta);
  const double poissonFactor = max(1., 1/Sqr(sigmaFactor));

  const Point core(x, y, 0, fgLocalCS);

  double logLikelihood = 0;

  for (int i = 0; i < nStations; ++i) {

    if (sIsCandidate[i]) {

      const double r = RPerp(fgShowerAxis, sPosition[i] - core);

      const double secDerivative = LDFFunctionSecondDerivative(fgLDFType, r, beta, gamma, fgNKGFermiMu, fgNKGFermiTau);
      if (!fgLowerLimit && secDerivative > fgNoRecovery)
        fgLowerLimit = true;
      const bool hasRecovery = (!fgLowerLimit && sSignalError[i]);
      
      const double particlesInStation = (sIsSaturation[i] && !hasRecovery) ?
        poissonFactor * (sSignal[i] - sSignalError[i]) :
        poissonFactor * sSignal[i];

      const double signalPrediction = s1000 * LDFFunction(fgLDFType, r, beta, gamma, fgNKGFermiMu, fgNKGFermiTau);
      const double particlePrediction = poissonFactor * signalPrediction;

      // if recovery used
      const double sigma = (sIsSaturation[i] && hasRecovery) ?
        poissonFactor * sqrt(Sqr(sSignalError[i]) + Sqr(sigmaFactor) * signalPrediction) :
        poissonFactor * sigmaFactor * sqrt(signalPrediction);

      // if saturation and no recovery use as lower limit
      if (sIsSaturation[i] && !hasRecovery && nStations > 3)

        logLikelihood += LogarithmOfNormalComplementCDF(particlesInStation, particlePrediction, sigma);

      else if (particlesInStation > 30)

        // Stations with large signals: approximately gaussian probability
        logLikelihood += LogarithmOfNormalPDF(particlesInStation, particlePrediction, sigma);

      else {

        // Stations with low signal: Poisson probability
        logLikelihood += particlesInStation * log(particlePrediction) - particlePrediction;

        double logarithmicFactor = 0;
        const int stop = int(particlesInStation + 0.5);
        for (int j = 1; j <= stop; ++j)
          logarithmicFactor += log(double(j));

        logLikelihood -= logarithmicFactor;

      }

    } else if (s1000 > 0) {

      const double r = RPerp(fgShowerAxis, sPosition[i] - core);
      const double signalPrediction = s1000 * LDFFunction(fgLDFType, r, beta, gamma, fgNKGFermiMu, fgNKGFermiTau);
      const double particlePrediction = poissonFactor * signalPrediction;

      logLikelihood += particlePrediction > 0.03 ?
        -particlePrediction +
        log(1 +                            // 0
            particlePrediction * (1 +      // 1
            0.5*particlePrediction * (1 +  // 2
            (1./3)*particlePrediction      // 3
           )))
        :
        -(1./24)*Sqr(Sqr(particlePrediction));

    }

  }

  // maximize log-likelihood = minimize -(log-likelihood)
  value = -logLikelihood;
}


void
LDFFinder::LDFFitChi2Fnc(int& /*nPar*/, double* const /*grad*/, double& value,
                         double* const par, const int iFlag)
{
  const double& x = par[1];
  const double& y = par[2];
  const double s1000 = par[0];
  double beta = 0;
  double gamma = 0;
  static int nStations = 0;
  static vector<bool> sIsCandidate;
  static vector<Point> sPosition;
  static vector<double> sSignal;
  static vector<double> sSignalError;
  static vector<double> sIsSaturation;

  const double cosTheta = fgShowerAxis.GetCosTheta(fgLocalCS);
  const double sigmaFactor = wcd::SignalUncertaintyFactor(wcd::eGAP2012_012, cosTheta);

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

    sIsCandidate.clear();
    sPosition.clear();
    sSignal.clear();
    sSignalError.clear();
    sIsSaturation.clear();

    const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

    for (SEvent::ConstStationIterator sIt = fgCurrentSEvent->StationsBegin();
         sIt != fgCurrentSEvent->StationsEnd(); ++sIt)

      if (sIt->IsCandidate()) {

        sIsCandidate.push_back(true);
        sPosition.push_back(sDetector.GetStation(*sIt).GetPosition());
        const bool saturation = sIt->IsLowGainSaturation();
        if (saturation && sIt->GetRecData().GetRecoveredSignal()) {
          sSignal.push_back(sIt->GetRecData().GetRecoveredSignal());
          sSignalError.push_back(sIt->GetRecData().GetRecoveredSignalError());
        } else {
          sSignal.push_back(sIt->GetRecData().GetTotalSignal());
          sSignalError.push_back(0.);
        }
        sIsSaturation.push_back(saturation);

      } else if (LDFFinder::fgUseSilentStations && sIt->IsSilent()) {

        sIsCandidate.push_back(false);
        sPosition.push_back(sDetector.GetStation(*sIt).GetPosition());
        sSignal.push_back(0.);

      }

    nStations = sIsCandidate.size();

  } // fnc init

  const Point core(x, y, 0, fgLocalCS);

  double chi2 = 0;

  if (LDFFinder::fgFixLDF)
    EstimateBetaGamma(beta, gamma, fgLDFType, fgShowerAxis.GetCosTheta(fgLocalCS), s1000);
  else {
    beta = par[3];
    gamma = par[4];
  }

  for (int i = 0; i < nStations; ++i)

    if (sIsCandidate[i]) {

      const double rho = RPerp(fgShowerAxis, sPosition[i] - core);
      const double predictedSignal = s1000 * LDFFunction(fgLDFType, rho, beta, gamma, fgNKGFermiMu, fgNKGFermiTau);
      const double secDerivative =
        LDFFunctionSecondDerivative(fgLDFType, rho, beta, gamma, fgNKGFermiMu, fgNKGFermiTau);
      const double sigma =
        (sIsSaturation[i] && sSignalError[i] && secDerivative < fgNoRecovery) ?
          sqrt(Sqr(sSignalError[i]) + Sqr(sigmaFactor) * predictedSignal) :
          sigmaFactor * sqrt(predictedSignal);

      // we do NOT add any contribution to the chi2, as it it would be a bias.
      if (!sIsSaturation[i] || (sSignalError[i] && secDerivative < 1)) {
        const double dev =
          (predictedSignal > numeric_limits<double>::min() ?
            (sSignal[i] - predictedSignal) / sigma :
            pow(numeric_limits<double>::max(), 0.25));
        const double dev2 = Sqr(dev);
        chi2 += dev2;
      }

    } else {

      const double rho = RPerp(fgShowerAxis, sPosition[i] - core);

      chi2 +=
        Fermi(rho, LDFFinder::fgSilentMaxRadius, LDFFinder::fgSilentRadiusTransition) *
          s1000 * LDFFunction(fgLDFType, rho, beta, gamma, fgNKGFermiMu, fgNKGFermiTau) / fgSilentSignalThreshold;

    }

  value = chi2;
}


bool
LDFFinder::EstimateCurvature(CurvatureFitInterface& cfi)
  const
{
  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
  const Point& core = cfi.fCore;

  double sx = 0;
  double sy = 0;
  double st = 0;
  double sm2 = 0;
  double dxx = 0;
  double dxy = 0;
  double dxt = 0;
  double dxm = 0;
  double dyy = 0;
  double dyt = 0;
  double dym = 0;
  double dtt = 0;
  double dtm = 0;
  int nStations = 0;

  for (SEvent::ConstCandidateStationIterator sIt = fgCurrentSEvent->CandidateStationsBegin();
       sIt != fgCurrentSEvent->CandidateStationsEnd(); ++sIt) {

    const Vector& xi = sDetector.GetStation(*sIt).GetPosition() - core;

    sx += xi.GetX(fgLocalCS);
    sy += xi.GetY(fgLocalCS);

    const double cti =
      kSpeedOfLight*(sIt->GetRecData().GetSignalStartTime() - fgBaryTime).GetInterval();

    st += cti;

    const double mi2 = 0.5*(xi.GetMag2() - Sqr(cti));

    sm2 += mi2;

    ++nStations;

    for (SEvent::ConstCandidateStationIterator sJt = fgCurrentSEvent->CandidateStationsBegin();
         sJt != fgCurrentSEvent->CandidateStationsEnd(); ++sJt)

      if (sIt->GetId() < sJt->GetId()) {

        const Vector& xj = sDetector.GetStation(*sJt).GetPosition() - core;

        const Vector xij(xi - xj);

        const double dxij = xij.GetX(fgLocalCS);
        const double dyij = xij.GetY(fgLocalCS);

        const double ctj =
          kSpeedOfLight*(sJt->GetRecData().GetSignalStartTime() - fgBaryTime).GetInterval();

        const double cdtij = cti - ctj;

        const double mj2 = 0.5*(xj.GetMag2() - Sqr(ctj));

        const double dmij2 = mi2 - mj2;

        dxx += Sqr(dxij);
        dxy += dxij*dyij;
        dxt += dxij*cdtij;
        dxm += dxij*dmij2;
        dyy += Sqr(dyij);
        dyt += dyij*cdtij;
        dym += dyij*dmij2;
        dtt += Sqr(cdtij);
        dtm += cdtij*dmij2;

    }

  }

  sx /= nStations;
  sy /= nStations;
  st /= nStations;
  sm2 /= nStations;

  HepMatrix delta(3,3);
  delta[0][0] = dxx;  delta[0][1] = dxy;  delta[0][2] = -dxt;
  delta[1][0] = dxy;  delta[1][1] = dyy;  delta[1][2] = -dyt;
  delta[2][0] = -dxt; delta[2][1] = -dyt; delta[2][2] = dtt;

  int ierr;
  delta.invert(ierr);

  HepVector k(3);
  k[0] = dxm;
  k[1] = dym;
  k[2] = -dtm;

  HepVector p(delta*k);

  const double& rcu = p[0];
  const double& rcv = p[1];
  const double& ct0 = p[2];

  const double rc2 = 2*(rcu*sx + rcv*sy - ct0*st - sm2) + Sqr(ct0);

  if (rc2 < 0) {
    cfi.fRc = 0;
    return false;
  } else {
    const double rc = sqrt(rc2);
    const double rcw2 = rc2 - Sqr(rcu) - Sqr(rcv);
    if (rcw2 < 0) {
      cfi.fRc = 0;
      return false;
    } else {
      const double rcw = sqrt(rcw2);
      Vector curvAxis(rcu,rcv,rcw, fgLocalCS);
      curvAxis.Normalize();
      cfi.fAxis = curvAxis;
      cfi.fRc = rc;
      cfi.fRcErr = 0;
      cfi.fct0 = ct0;
      cfi.fct0Err = 0;
      cfi.fSigmaU2 = 0;
      cfi.fSigmaUV = 0;
      cfi.fSigmaV2 = 0;
    }
  }

  return true;
}


namespace LDFFinderOG {

  Point gCore;

  int complexWCounter = 0;
  const char* complexWMessage = "complex w";

}


double
LDFFinder::FitCurvatureDriver(CurvatureFitInterface& cfi)
  const
{
  gCore = cfi.fCore;

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

  complexWCounter = 0;
  minuit.SetFCN(LDFFinder::CurvatureFitFnc);
  argList[0] = 1;
  minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);

  minuit.mnparm(0, "u", cfi.fAxis.GetX(fgLocalCS), 0.01, 0, 0, errFlag);
  minuit.mnparm(1, "v", cfi.fAxis.GetY(fgLocalCS), 0.01, 0, 0, errFlag);
  minuit.mnparm(2, "Rc", cfi.fRc, 10*meter, 0, 0, errFlag);
  minuit.mnparm(3, "ct0", cfi.fct0+ cfi.fRc, 10*meter, 0, 0, errFlag);

  if (cfi.fFixedAxis) {
    minuit.FixParameter(0);
    minuit.FixParameter(1);
  }

  if (cfi.fFixedRc)
    minuit.FixParameter(2);

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

  bool fatalError = false;

  argList[0] = 500;
  minuit.mnexcm("MINIMIZE", argList, 1, errFlag);
  if (errFlag) {
    OUT(eObscure)
      << "  minuit minimize error: " << errFlag << '\n';
    minuit.mnexcm("MINOS", argList, 0, errFlag);
    if (errFlag) {
      OUT(eObscure)
        << "  minuit minos error: " << errFlag << '\n';
      fatalError = true;
    }
  }

  if (complexWCounter > 1) {
    ostringstream warn;
    warn << "Message \'" << complexWMessage << "\' repeated "
         << complexWCounter << " times.";
    WARNING(warn);
  }

  if (fatalError)
    return 0;

  // Get results
  double u;
  double uErr;
  minuit.GetParameter(0, u, uErr);
  double v;
  double vErr;
  minuit.GetParameter(1, v, vErr);
  const double w2 = 1 - Sqr(u) - Sqr(v);
  if (w2 < 0) {
    OUT(eIntermediate)
      << "  Fit failed due to unphysical angle cosines.\n";
    return 0;
  }
  cfi.fAxis = Vector(u,v,sqrt(w2), fgLocalCS);
  cfi.fAxisErr = Vector(uErr,vErr,0, fgLocalCS);

  minuit.GetParameter(2, cfi.fRc, cfi.fRcErr);
  if (cfi.fRc <= 0) {
    OUT(eIntermediate)
      << "  Fit failed due to Rc <= 0.\n";
    return 0;
  }
  if (!cfi.fFixedRc && cfi.fRc > 100*kilometer) {
    OUT(eIntermediate)
      << "  Rc too large, reverting to plane fit.\n";
    return 0;
  }

  double ct0;
  double ct0Err;
  minuit.GetParameter(3, ct0, ct0Err);
  cfi.fct0 = ct0 - cfi.fRc;
  cfi.fct0Err = sqrt(Sqr(ct0Err) + Sqr(cfi.fRcErr));

  // error matrix
  double emat[4][4];
  minuit.mnemat(&emat[0][0], 4);
  cfi.fSigmaU2 = emat[0][0];
  cfi.fSigmaV2 = emat[1][1];
  cfi.fSigmaUV = emat[0][1];

  return minuit.fAmin / Sqr(kSpeedOfLight);
}


namespace LDFFinderOG {

  struct CurvatureFitStationInfo {
    double fCTiming;
    double fT50;
    double fSignal;
    Vector fPosition;
  };

}


void
LDFFinder::CurvatureFitFnc(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& rc = par[2];
  const double& ct0 = par[3];

  static vector<CurvatureFitStationInfo> sStationInfo;
  static const sdet::STimeVariance* timeVariance = 0;

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

    sStationInfo.clear();
    const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

    for (SEvent::ConstCandidateStationIterator sIt = fgCurrentSEvent->CandidateStationsBegin();
         sIt != fgCurrentSEvent->CandidateStationsEnd(); ++sIt) {

      const StationRecData& sRec = sIt->GetRecData();
      const CurvatureFitStationInfo si = {
        kSpeedOfLight * ((sRec.GetSignalStartTime() - fgBaryTime).GetInterval()),
        sRec.GetT50(),
        sRec.GetTotalSignal(),
        sDetector.GetStation(*sIt).GetPosition() - gCore
      };
      sStationInfo.push_back(si);

    }

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

  } // fnc init

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

  if (w2 < 0) {
    if (!complexWCounter)
      WARNING(complexWMessage);
    ++complexWCounter;
    value = 1e30;
    return;
  }

  const double w = sqrt(w2);

  const Vector axis = Vector(u,v,w, fgLocalCS);
  const Vector rcAxis = rc * axis;

  double chi2 = 0;

  for (vector<CurvatureFitStationInfo>::const_iterator siIt = sStationInfo.begin();
       siIt != sStationInfo.end(); ++siIt) {

    const Vector stationCenter = rcAxis - siIt->fPosition;
    const double drc = stationCenter.GetMag();

    const double stationCosTheta = stationCenter.GetZ(fgLocalCS) / drc;

    const double sigma2 =
      timeVariance->GetTimeSigma2(siIt->fSignal, siIt->fT50, stationCosTheta, sqrt(RPerp2(axis, siIt->fPosition)));

    const double cdt = siIt->fCTiming - ct0 - drc + rc;

    chi2 += Sqr(cdt) / sigma2;

  }

  value = chi2 / Sqr(kSpeedOfLight);
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: