SdHorizontalReconstruction.cc

Go to the documentation of this file.

#include "SdHorizontalReconstruction.h"
#include "FitInterface.h"
#include "ShowerFrontFunction.h"
#include "ShowerSizeFunction.h"
#include "Utilities.h"

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

#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerUnits.h>
#include <utl/String.h>
#include <utl/NumericalErrorPropagator.h>
#include <utl/Transformation.h>
#include <utl/EquationSolver.h>

#include <tls/MuonProfile.h>
#include <tls/VTankResponse.h>
#include <tls/TankResponseFactory.h>
#include <tls/EMComponent.h>

#include <det/Detector.h>

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

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

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

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>
#include <fevt/FdComponentSelector.h>

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

#include <Minuit2/MnMinimize.h>
#include <Minuit2/MnHesse.h>
#include <Minuit2/MnMinos.h>
#include <Minuit2/MnMigrad.h>
#include <Minuit2/FCNBase.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnUserParameters.h>
#include <Minuit2/MnUserCovariance.h>
#include <Minuit2/MnPrint.h>

#include <iomanip>

using namespace std;
using namespace SdHorizontalReconstructionNS;
using namespace fwk;
using namespace utl;
using namespace evt;
using namespace sevt;
using namespace fevt;
using namespace tls;

class FdGeometryPropagator : public NumericalErrorPropagator {
public:
  FdGeometryPropagator(const SdHorizontalReconstruction& parent,
                       const CoordinateSystemPtr& eyeCS)
    : fParent(parent), fEyeCS(eyeCS)
  {}

  void
  Transform(vector<double>& output, const vector<double>& input)
    const
  {
    const double sDPtheta = input[0];
    const double sDPphi   = input[1];
    const double rp       = input[2];
    const double chi0     = input[3];

    const Point eyePos(0, 0, 0, fEyeCS);

    // sdp vector: vertical on SDP plane
    const Vector sdp(1.0, sDPtheta, sDPphi, fEyeCS, Vector::kSpherical);

    const Vector vertical(0, 0, 1, fEyeCS);
    Vector horizontalInSDP = cross(sdp, vertical); // horizontal in eyeCS
    horizontalInSDP.Normalize();

    const Vector core_eye_vec = rp / sin(kPi - chi0) * horizontalInSDP;
    Point core = eyePos + core_eye_vec; // preliminary core

    Transformation rot(Transformation::Rotation(-chi0, sdp, fEyeCS));

    const Vector axis(rot * horizontalInSDP); // is normalized

    // shift core along axis to barycenter altitude
    core -= axis * core.GetZ(fParent.fBaryCS)/axis.GetZ(fParent.fBaryCS);

    output[eCoreXExt] = core.GetX(fParent.fBaryCS);
    output[eCoreYExt] = core.GetY(fParent.fBaryCS);
    output[eThetaExt] = axis.GetTheta(fParent.fBaryCS);
    output[ePhiExt] = axis.GetPhi(fParent.fBaryCS);
  }

private:
  const SdHorizontalReconstruction& fParent;
  const CoordinateSystemPtr& fEyeCS;
};


class FinalPropagator : public NumericalErrorPropagator {
public:
  FinalPropagator(const SdHorizontalReconstruction& parent,
                  const vector<Point>& posList)
    : fParent(parent), fPosList(posList)
  {}

  void
  Transform(vector<double>& output, const vector<double>& input)
    const
  {
    const double coreX = input[0];
    const double coreY = input[1];
    const double theta = input[2];
    const double phi = input[3];
    const double distance = input[4];
    const double ctime = input[5];

    Point core, origin;
    fParent.GetShowerAxis(core, origin, coreX, coreY, theta, phi, distance);
    const Vector showerAxis = origin - core;

    const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(core);

    output[0] = showerAxis.GetTheta(coreCS);
    output[1] = showerAxis.GetPhi(coreCS);
    output[2] = showerAxis.GetMag();
    output[3] = ctime + showerAxis.GetMag();

    const size_t n = fPosList.size();
    for (size_t i = 0; i < n; ++i)
      output[4 + i] = GetRho(fPosList[i], core, origin);
  }

protected:
  const SdHorizontalReconstruction& fParent;
  const vector<Point>& fPosList;
};

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

  Branch topB = cc->GetTopBranch("SdHorizontalReconstruction");

  topB.GetChild("MaxRadiusNonTriggeringStations").GetData(fMaxSilentRadius);

  topB.GetChild("MaxRelativeDistanceNonTriggeringStations").GetData(fMaxSilentRelDistance);

  topB.GetChild("SilentSignalThreshold").GetData(fSilentSignalThreshold);

  topB.GetChild("DropStationsBelowThreshold").GetData(fDropStationsBelowThreshold);

  topB.GetChild("MinTheta").GetData(fMinTheta);
  topB.GetChild("MaxTheta").GetData(fMaxTheta);

  string showerFrontModel;
  topB.GetChild("ShowerFrontModel").GetData(showerFrontModel);
  if (showerFrontModel == "MuonArrival")
    fShowerFrontModel = eMuonArrival;
  else
    fShowerFrontModel = eSphere;

  topB.GetChild("MaxDistanceCoreToBarycenter").GetData(fMaxDistanceCoreToBarycenter);

  topB.GetChild("MaxThetaDifference").GetData(fMaxThetaDifference);

  topB.GetChild("MaxTimeSigma").GetData(fMaxTimeSigma);

  topB.GetChild("MaxDistanceStationToExternalCore").GetData(fMaxDistanceStationToExternalCore);

  topB.GetChild("MinStationsForEnergyFit").GetData(fMinStationsForEnergyFit);

  topB.GetChild("MinStationsForProductionDistanceFit").GetData(fMinStationsForProductionDistanceFit);

  string muonProfile;
  topB.GetChild("MuonProfile").GetData(muonProfile);
  const Branch muonProfileBranch = cc->GetTopBranch("MuonProfile").GetChild(muonProfile.c_str());
  fMuonProfile = boost::shared_ptr<MuonProfile>(new MuonProfile(muonProfileBranch));

  string tankResponse;
  topB.GetChild("TankResponse").GetData(tankResponse);
  const Branch tankResponseBranch = cc->GetTopBranch("TankResponse").GetChild(tankResponse.c_str());
  fTankResponse = &TankResponseFactory::GetInstance().
    GetTankResponse(tankResponseBranch);

  string emComponent;
  topB.GetChild("EMComponent").GetData(emComponent);
  const Branch EMComponentBranch = cc->GetTopBranch("EMComponent").GetChild(emComponent.c_str());
  fEMComponent = boost::shared_ptr<EMComponent>(new EMComponent(EMComponentBranch));

  topB.GetChild("N19Correction").GetData(fN19Correction);

  topB.GetChild("EnergyConstant").GetData(fEnergyConstant);
  topB.GetChild("Gamma")         .GetData(fGamma);

  fExternalGeometry = eNone;
  string externalGeometry;
  topB.GetChild("ExternalGeometry").GetData(externalGeometry);
  if (externalGeometry == "MC") fExternalGeometry = eMC;
  if (externalGeometry == "FixHybrid") fExternalGeometry = eFixHybrid;
  if (externalGeometry == "SoftHybrid") fExternalGeometry = eSoftHybrid;

  topB.GetChild("PropagateAxisUncertainty").GetData(fPropagateAxisUncertainty);

  topB.GetChild("CurvatureModelParameters").GetData(fCurvModelPars);

  ostringstream msg;
  msg << "\nParameters:\n"
      << "  Accepted theta angle range               : [" << fMinTheta/degree << " , " << fMaxTheta/degree << " ] deg" << '\n'
      << "  Max. radius of non-triggering stations   : " << fMaxSilentRadius/km << " km\n"
      << "  Max. rel. distance of non-trig. stations : " << fMaxSilentRelDistance/km << " km\n"
      << "  Silent signal threshold                  : " << fSilentSignalThreshold << " VEM\n"
      << "  Drop stations below threshold            : " << (fDropStationsBelowThreshold ? "yes" : "no") << "\n"
      << "  Allowed distance core,barycenter         : " << fMaxDistanceCoreToBarycenter/km << " km"
      << (fMaxDistanceCoreToBarycenter == 0 ? " (check disabled)\n" : "\n")
      << "  Allowed theta change                     : " << fMaxThetaDifference/degree << " deg"
      << (fMaxThetaDifference == 0 ? " (check disabled)\n" : "\n")
      << "  Required stations for energy fit         : " << fMinStationsForEnergyFit << '\n'
      << "  Required stations for prod. distance fit : " << fMinStationsForProductionDistanceFit << '\n'
      << "  ShowerFrontModel                         : " << showerFrontModel << '\n'
      << "  MuonProfile                              : " << muonProfileBranch.GetName()
      << " Theta limits [" << fMuonProfile->GetThetaMin()/degree << " , "
      << fMuonProfile->GetThetaMax()/degree << " ] deg" << "\n"
      << "  TankResponse                             : " << tankResponseBranch.GetName()
      << " Theta limits [" << fTankResponse->GetThetaMin()/degree << " , "
      << fTankResponse->GetThetaMax()/degree << " ] deg" << "\n"
      << "  EMComponent                              : " << EMComponentBranch.GetName()
      << " Theta limits [" << fEMComponent->GetThetaMin()/degree << " , "
      << fEMComponent->GetThetaMax()/degree << " ] deg" << "\n"
      << "  Applying empirical correction N19-unbiased = (1 + A + B*lg(N19))*N19\n"
      << "    A : " << fN19Correction[0] << "\n"
      << "    B : " << fN19Correction[1] << "\n"
      << "  Energy calibration\n"
      << "    EnergyConstant: " << fEnergyConstant/EeV << " EeV\n"
      << "    Gamma         : " << fGamma << "\n"
      << "  External geometry                        : " << externalGeometry << "\n"
      << "  Max. time offset [ext.geom.]             : " << fMaxTimeSigma << " sigma\n"
      << "  Allowed station distance [ext.geom.]     : " << fMaxDistanceStationToExternalCore/km << " km"
      << (fMaxDistanceStationToExternalCore == 0 ? " (check disabled)" : "");
  INFO(msg);

  return eSuccess;
}

VModule::ResultFlag
SdHorizontalReconstruction::Run(evt::Event & event)
{
  {
    ostringstream msg;
    msg << "Event " << event.GetHeader().GetId();
    INFO(msg);
  }

  // nothing to do if there is no SEvent with trigger
  if (!event.HasSEvent())
    return eContinueLoop;
  if (!event.GetSEvent().HasTrigger())
    return eContinueLoop;

  SEvent& sEvent = event.GetSEvent();

  sEvent.SortStations(ByDecreasingSignal());

  SizeData sizeData;
  AxisData axisData;
  StationList list;
  SilentStationList slist;
  ExternalGeometryData extGeomData;

  size_t nCandidates = 0;

  if (!UpdateBarycenter(nCandidates, sEvent))
  {
    ERROR("Could not find barycenter!");
    return eContinueLoop;
  }

  if (nCandidates < fMinStationsForEnergyFit) {
    ostringstream msg;
    msg << "Event has less than " << fMinStationsForEnergyFit
        << " stations. Dropping it...";
    INFO(msg);
    return eContinueLoop;
  }

  if (fExternalGeometry == eNone) {
    // angular reco is required prior to this module
    if (!event.HasRecShower() && !event.GetRecShower().HasSRecShower())
      return eContinueLoop;
    const ShowerSRecData& sRecSh = event.GetRecShower().GetSRecShower();

    // initialize core to surface-shifted barycenter
    // initialize axis to whatever the previous module set
    sizeData.fPar[eCoreX] = 0;
    sizeData.fPar[eCoreY] = 0;
    const Vector& axis = sRecSh.GetAxis();
    axisData.fPar[eTheta] = axis.GetTheta(fBaryCS);
    axisData.fPar[ePhi] = axis.GetPhi(fBaryCS);
    axisData.fPar[eDistance] = ExpectedCurvatureRadius(axisData.fPar[eTheta]);
    axisData.fPar[eCTime] = -axisData.fPar[eDistance];
    axisData.fCov[eTheta] = Sqr(sRecSh.GetThetaError());
    axisData.fCov[ePhi] = Sqr(sRecSh.GetPhiError());
    axisData.fCov(eTheta,ePhi) = sRecSh.GetCorrelationThetaPhi()*sRecSh.GetThetaError()*sRecSh.GetPhiError();
    axisData.fChi2 = event.GetRecShower().GetSRecShower().GetAngleChi2();
  }
  else
  {
    if (!GetExternalGeometry(extGeomData, event))
      return eContinueLoop;

    if (!event.HasRecShower())
      event.MakeRecShower();

    if (!event.GetRecShower().HasSRecShower())
      event.GetRecShower().MakeSRecShower();

    // initialize core to external core
    // initialize axis to external axis
    sizeData.fPar[eCoreX] = extGeomData.fPar[eCoreXExt];
    sizeData.fPar[eCoreY] = extGeomData.fPar[eCoreYExt];
    axisData.fPar[eTheta] = extGeomData.fPar[eThetaExt];
    axisData.fPar[ePhi] = extGeomData.fPar[ePhiExt];
    axisData.fPar[eDistance] = ExpectedCurvatureRadius(axisData.fPar[eTheta]);
    axisData.fPar[eCTime] = -axisData.fPar[eDistance];

    if (!CleanEvent(nCandidates, sEvent, axisData, sizeData))
      return eContinueLoop;
  }

  // check whether event is outside of valid theta range
  if (axisData.fPar[eTheta] < fMinTheta || fMaxTheta < axisData.fPar[eTheta]) {
    ostringstream msg;
    msg << "  theta = " << axisData.fPar[eTheta]/degree << " deg outside range ["
        << fMinTheta/degree << ", " << fMaxTheta/degree << "] deg, skipping...";
    WARNING(msg);
    return eContinueLoop;
  }

  Point core, origin;
  GetShowerAxis(core, origin, sizeData, axisData);
  PrepareStationData(list, slist, core, origin, sEvent);

  ostringstream msg;
  msg << "\n  SD event " << sEvent.GetHeader().GetId() << "\n"
      << "  Theta " << axisData.fPar[eTheta]/deg << " deg (bary)" << '\n'
      << "  # candidate stations = " << nCandidates << '\n'
      << "  # included non-triggering stations = " << slist.size();
  INFO(msg);

  {
    const bool fixAxis = (fExternalGeometry == eMC || fExternalGeometry == eFixHybrid);

    AxisData update = axisData;
    INFO("Fit shower front with fixed production distance");
    if (FitShowerFront(core, update, extGeomData, list, eFixCurvature | (fixAxis ? eFixAxis : 0)))
      axisData = update;

    if (nCandidates >= fMinStationsForProductionDistanceFit)
    {
      AxisData update = axisData;
      INFO("Fit shower front with free production distance");
      if (FitShowerFront(core, update, extGeomData, list, fixAxis ? eFixAxis : 0))
        axisData = update;
    }
  }

  const bool haveCurv = axisData.fCov[eDistance] > 0;

  INFO("Fit N19 with core fixed to barycenter, fixed axis, approximate chi2 method");
  sizeData.fPar[eN19] = 1.0;
  if (!FitShowerSize(sizeData, axisData, extGeomData, list, slist, eFixCore | eFixAxis | eApproximate))
  {
    ERROR("N19 estimate failed, this should never happen! Skipping event...");
    return eContinueLoop;
  }

  if (fExternalGeometry == eMC || fExternalGeometry == eFixHybrid)
  {
    INFO("Fit N19 with fixed external geometry (core, axis)");
    SizeData update = sizeData;
    if (FitShowerSize(update, axisData, extGeomData, list, slist, eFixCore | eFixAxis))
    {
      update.fRecStage = haveCurv ?
        ShowerSRecData::kLDFSilentsAndCurvature : ShowerSRecData::kLDFSilents;
      sizeData = update;
    }
  }
  else
  {
    INFO("Fit N19 with free core, fixed axis, approximate chi2 method");
    if (!FitShowerSize(sizeData, axisData, extGeomData, list, slist, eFixAxis | eApproximate))
    {
      ERROR("N19 estimate failed, this should never happen! Skipping event...");
      return eContinueLoop;
    }

    SizeData update = sizeData;
    INFO("Fit N19 with free core, fixed axis");
    if (FitShowerSize(update, axisData, extGeomData, list, slist, eFixAxis))
    {
      update.fRecStage = haveCurv ?
        ShowerSRecData::kLDFSilentsAndCurvature : ShowerSRecData::kLDFSilents;
      sizeData = update;

      if (fPropagateAxisUncertainty)
      {
        INFO("Fit N19 with free core, free axis");
        if (FitShowerSize(update, axisData, extGeomData, list, slist))
        {
          update.fRecStage = ShowerSRecData::kLDFGlobalFit;
          sizeData = update;
        }
      }
    }
  }

  // apply empirical N19 correction
  const double correctionFactor =
    1.0 + fN19Correction[0] + fN19Correction[1]*log10(sizeData.fPar[eN19]);
  sizeData.fPar[eN19] *= correctionFactor;
  for (size_t i = 0; i < SizeData::size; ++i)
    sizeData.fCov(eN19,i) *= correctionFactor;
  sizeData.fCov(eN19,eN19) *= correctionFactor;

  SetReconstructedValues(event, sizeData, axisData, list, slist, extGeomData);

  return eSuccess;
}


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


bool
SdHorizontalReconstruction::UpdateBarycenter(size_t& nCandidates,
                                             const SEvent& sEvent)
{
  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
  const CoordinateSystemPtr siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();
  const Point siteOrigin(0,0,0, siteCS);
  const TimeStamp eventTime(sEvent.GetTrigger().GetTime());
  Vector barySum(0,0,0, siteCS);
  double timeSum = 0;
  double weightSum = 0;
  nCandidates = 0;
  for (SEvent::ConstCandidateStationIterator
         sIt = sEvent.CandidateStationsBegin(),
         sEnd = sEvent.CandidateStationsEnd();
       sIt != sEnd; ++sIt)
  {
    const StationRecData& sRec = sIt->GetRecData();
    const double weight = sqrt(sRec.GetTotalSignal());
    weightSum += weight;
    barySum += weight * (sDetector.GetStation(*sIt).GetPosition() - siteOrigin);
    timeSum += weight * (sRec.GetSignalStartTime() - eventTime).GetInterval();
    if (!(fDropStationsBelowThreshold && (sRec.GetTotalSignal() < fSilentSignalThreshold)))
      ++nCandidates;
  }

  if (!weightSum) {
    FATAL("sum of weights = 0!");
    return false;
  }

  fBary = siteOrigin + barySum/weightSum;
  fBaryTime = eventTime + timeSum/weightSum;

  // barycenter usually is under Earth's surface
  // shift barycenter vertically to altitude of closest station
  double bestR2 = numeric_limits<double>::max();
  Point bestPos;
  for (sdet::SDetector::StationIterator
         sIt = sDetector.AllStationsBegin(),
         sEnd = sDetector.AllStationsEnd();
       sIt != sEnd; ++sIt)
  {
    const double r2 = (sIt->GetPosition() - fBary).GetMag2();
    if (r2 < bestR2)
    {
      bestR2 = r2;
      bestPos = sIt->GetPosition();
    }
  }

  const CoordinateSystemPtr prelimCS = LocalCoordinateSystem::Create(fBary);
  fBary += Vector(0, 0, bestPos.GetZ(prelimCS), prelimCS);

  fBaryCS = LocalCoordinateSystem::Create(fBary);
  // approximately...
  fEarthCenter = Point(0, 0, -(kEarthRadius + 1.4*km), fBaryCS);

  ostringstream msg;
  msg << "\n  # candidate stations = " << nCandidates << "\n"
      << "  barycenter = ( " << fBary.GetX(siteCS)/km
      << ", " << fBary.GetY(siteCS)/km << ", " << fBary.GetZ(siteCS)/km
      << ") km (site)\n"
      << "  bary time = " << (timeSum/weightSum)/nanosecond
      << " [ns] (to event time)";
  INFO(msg);

  return true;
}


bool
SdHorizontalReconstruction::GetExternalGeometry(ExternalGeometryData& extGeomData,
                                                const Event& event)
{
  // need to reset baryCS: event may still contain accidental stations
  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

  switch (fExternalGeometry)
  {
    case eMC:
    {
      if (!event.HasSimShower()) {
        ERROR("Event contains no simulation data.");
        return false;
      }

      const ShowerSimData& simShower = event.GetSimShower();

      const Vector axis = -simShower.GetDirection();
      Point core = simShower.GetPosition();

      const CoordinateSystemPtr prelimCS = LocalCoordinateSystem::Create(core);

      // core may be floating, shift to altitude of closest station
      double bestR2 = numeric_limits<double>::max();
      double bestZ = 0.0;
      for (sdet::SDetector::StationIterator
             sIt = sDetector.AllStationsBegin(),
             sEnd = sDetector.AllStationsEnd();
           sIt != sEnd; ++sIt)
      {
        const double r2 = (sIt->GetPosition() - core).GetMag2();
        if (r2 < bestR2)
        {
          bestR2 = r2;
          bestZ = sIt->GetPosition().GetZ(prelimCS);
        }
      }

      core += axis * bestZ/axis.GetZ(prelimCS);
      fBary = core;
      fBaryCS = LocalCoordinateSystem::Create(fBary);

      extGeomData.fPar[eCoreXExt] = 0;
      extGeomData.fPar[eCoreYExt] = 0;
      extGeomData.fPar[eThetaExt] = axis.GetTheta(fBaryCS);
      extGeomData.fPar[ePhiExt] = axis.GetPhi(fBaryCS);
    }
    break;

    case eFixHybrid:
    case eSoftHybrid:
    {
      if (!event.HasFEvent()) {
        ERROR("Event contains no data from a hybrid reconstruction.");
        return false;
      }

      // compute average over available Fd geometries

      Point core(0, 0, 0, fBaryCS);
      Vector axis(0, 0, 0, fBaryCS);

      int nEyes = 0;
      for (FEvent::ConstEyeIterator
             eIt = event.GetFEvent().EyesBegin(ComponentSelector::eHasData),
             eEnd = event.GetFEvent().EyesEnd(ComponentSelector::eHasData);
           eIt != eEnd; ++eIt)
      {
        const EyeRecData& r = eIt->GetRecData();
        if (r.GetTimeFitChiSquare() >= 10*sqrt(2*r.GetTimeFitNDof())) continue;
        ++nEyes;
        const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(*eIt);
        const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
        const Point eyePos(0, 0, 0, eyeCS);

        // sdp vector: vertical on SDP plane
        const Vector sdp(1.0,
                         r.GetSDP().GetTheta(eyeCS),
                         r.GetSDP().GetPhi(eyeCS),
                         eyeCS, Vector::kSpherical);

        const Vector vertical(0, 0, 1, eyeCS);
        Vector horizontalInSDP = cross(sdp, vertical); // horizontal in eyeCS
        horizontalInSDP.Normalize();

        const Vector core_eye_vec = r.GetRp() / sin(kPi - r.GetChiZero()) * horizontalInSDP;
        core += (eyePos-fBary) + core_eye_vec;

        Transformation rot(Transformation::Rotation(-r.GetChiZero(), sdp, eyeCS));

        axis += rot * horizontalInSDP; // is normalized
      }

      if (!nEyes)
      {
        ERROR("Event contains no data from a hybrid reconstruction.");
        return false;
      }

      core /= nEyes;
      axis /= nEyes;

      const CoordinateSystemPtr prelimCS = LocalCoordinateSystem::Create(core);

      // core may be floating, shift to altitude of closest station
      const SEvent& sEvent = event.GetSEvent();
      const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
      Point bestPos = fEarthCenter;
      for (SEvent::ConstCandidateStationIterator
             sIt = sEvent.CandidateStationsBegin(),
             sEnd = sEvent.CandidateStationsEnd();
           sIt != sEnd; ++sIt)
      {
        const Point& pos = sDetector.GetStation(sIt->GetId()).GetPosition();
        if ((pos - core).GetMag2() < (bestPos - core).GetMag2())
          bestPos = pos;
      }

      core += axis * bestPos.GetZ(prelimCS)/axis.GetZ(prelimCS);
      if (fMaxDistanceStationToExternalCore > 0 &&
          (bestPos - core).GetMag() > fMaxDistanceStationToExternalCore)
      {
        ostringstream msg;
        msg << "Event rejected: station with signal " << (bestPos-core).GetMag()/km
            << " km away from hybrid core [limit: " <<  fMaxDistanceStationToExternalCore/km << " km]";
        ERROR(msg);
        return false;
      }

      fBary = core;
      fBaryCS = LocalCoordinateSystem::Create(fBary);

      for (FEvent::ConstEyeIterator
             eIt = event.GetFEvent().EyesBegin(ComponentSelector::eHasData),
             eEnd = event.GetFEvent().EyesEnd(ComponentSelector::eHasData);
           eIt != eEnd; ++eIt)
      {
        const EyeRecData& r = eIt->GetRecData();
        if (r.GetTimeFitChiSquare() >= 10*sqrt(2*r.GetTimeFitNDof())) continue;
        const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(*eIt);
        const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
        FdGeometryPropagator propagate(*this, eyeCS);

        vector<double> iPar(4), oPar(4);
        CovarianceMatrix iCov(4), oCov(4);

        enum Parameters { eTheta = 0, ePhi = 1, eRp = 2, eChi0 = 3 };
        iPar[eTheta] = r.GetSDP().GetTheta(eyeCS);
        iPar[ePhi] = r.GetSDP().GetPhi(eyeCS);
        iPar[eRp] = r.GetRp();
        iPar[eChi0] = r.GetChiZero();
        iCov[eTheta] = Sqr(r.GetSDPThetaError());
        iCov[ePhi] = Sqr(r.GetSDPPhiError());
        iCov[eRp] = Sqr(r.GetRpError());
        iCov[eChi0] = Sqr(r.GetChiZeroError());
        iCov(eTheta, ePhi) = r.GetSDPCorrThetaPhi()*sqrt(iCov[eTheta]*iCov[ePhi]);
        iCov(eRp, eChi0) = r.GetRpChi0Correlation()*sqrt(iCov[eRp]*iCov[eChi0]);

        propagate(oPar, oCov, iPar, iCov);

        for (size_t i = 0; i < 4; ++i)
        {
          extGeomData.fPar[i] += oPar[i];
          for (size_t j = 0; j < 4; ++j)
            extGeomData.fInvCov[i][j] += oCov(i,j);
        }
      }

      for (size_t i = 0; i < 4; ++i)
      {
        extGeomData.fPar[i] /= nEyes;
        for (size_t j = 0; j < 4; ++j)
          extGeomData.fInvCov[i][j] /= nEyes;
      }

      try {
        InvertMatrix(4, extGeomData.fInvCov);
      }
      catch (const DoesNotComputeException& e)
      {
        return false;
      }
    }
    break;
    default:
    return false;
  }

  return true;
}


bool
SdHorizontalReconstruction::CleanEvent(std::size_t& nCandidates,
                                          sevt::SEvent& sEvent,
                                          const AxisData& axisData,
                                          const SizeData& sizeData)
  const
{ // reject stations that are grossly inconsistent with external geometry
  Point core, origin;
  GetShowerAxis(core, origin, sizeData, axisData);
  ShowerFrontFunction frontFcn(*this, core, StationList(), ExternalGeometryData());

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

  // pick most probable station as reference
  double bestWeight = 0;
  int bestId = 0;
  double bestSigmaCT = 0;
  for (SEvent::ConstCandidateStationIterator
         sIt = sEvent.CandidateStationsBegin(),
         sEnd = sEvent.CandidateStationsEnd();
       sIt != sEnd; ++sIt)
  {
    const int sId = sIt->GetId();
    const double s = sIt->GetRecData().GetTotalSignal();
    const double rho = GetRho(sDetector.GetStation(*sIt).GetPosition(), core, origin);
    const double weight = s/rho;
    if (bestWeight < weight)
    {
      bestWeight = weight;
      bestId = sId;
    }
  }

  // use it to set originCT
  double originCT = 0;
  if (bestId > 0)
  {
    ostringstream msg;
    msg << "Station " << bestId << " serves as time reference to clean event";
    INFO(msg);
    const StationRecData& sRec = sEvent.GetStation(bestId).GetRecData();
    StationData sd;
    sd.fPos = sDetector.GetStation(bestId).GetPosition();
    sd.fSignal = sRec.GetTotalSignal();
    sd.fCTime = kSpeedOfLight *
      (sRec.GetSignalStartTime() - fBaryTime).GetInterval();
    sd.fT50       = sRec.GetT50();
    double meanCT;
    frontFcn.Predict(meanCT, bestSigmaCT, sd, origin, 0);
    originCT = sd.fCTime-meanCT;
  }
  else
  {
    ERROR("No reference station found");
    return false;
  }

  // reject stations based on shower front model
  for (SEvent::CandidateStationIterator
         sIt = sEvent.CandidateStationsBegin(),
         sEnd = sEvent.CandidateStationsEnd();
       sIt != sEnd; ++sIt)
  {
    const StationRecData& sRec = sIt->GetRecData();

    StationData sd;
    sd.fPos = sDetector.GetStation(*sIt).GetPosition();
    sd.fSignal    = sRec.GetTotalSignal();
    sd.fCTime     = kSpeedOfLight *
      (sRec.GetSignalStartTime() - fBaryTime).GetInterval();
    sd.fT50       = sRec.GetT50();

    double meanCT, sigmaCT;
    frontFcn.Predict(meanCT, sigmaCT, sd, origin, originCT);
    sigmaCT = sqrt(Sqr(sigmaCT) + Sqr(bestSigmaCT));

    if (abs(sd.fCTime - meanCT) > fMaxTimeSigma*sigmaCT)
    {
      ostringstream msg;
      msg << "Station " << sIt->GetId() << " rejected:"
             " time residual to external geometry = "
          << fixed << setprecision(2)
          << (sd.fCTime - meanCT)/sigmaCT << " sigma [limit: " << fMaxTimeSigma << "]";
      INFO(msg);
      sIt->SetRejected(sevt::StationConstants::eOutOfTime);
      if (!(fDropStationsBelowThreshold && (sd.fSignal < fSilentSignalThreshold)))
        --nCandidates;
    }
  }

  return nCandidates > 0;
}


bool
SdHorizontalReconstruction::FitShowerFront(const Point& core,
                                           AxisData& ad,
                                           const ExternalGeometryData& gd,
                                           const StationList& list,
                                           const int option)
  const
{
  using namespace ROOT::Minuit2;

  MnUserParameters pars;

  if (ad.fPar[eTheta] >= 89*degree) ad.fPar[eTheta] = 88*degree;

  pars.Add("theta", ad.fPar[eTheta], 0.01);         // 0
  pars.Add("phi", ad.fPar[ePhi], 0.01);             // 1
  pars.Add("distance", ad.fPar[eDistance]/km, 0.1); // 2
  pars.Add("ctime", ad.fPar[eCTime]/km, 0.1);       // 3

  pars.SetLimits("theta", 0, 89*degree);
  pars.SetLowerLimit("distance", 0);

  ad.fNFree = 4;
  if (option & eFixDirection)
  {
    pars.Fix("theta");
    pars.Fix("phi");
    ad.fNFree -= 2;
  }

  if (option & eFixCurvature)
  {
    pars.Fix("distance");
    ad.fNFree -= 1;
  }

  ShowerFrontFunction fcn(*this, core, list, gd);

  MnMinimize m(fcn, pars, 1);

  FunctionMinimum fmin = m(10000);

  if (fmin.IsAboveMaxEdm() || fmin.HasReachedCallLimit() || !fmin.IsValid())
  {
    ERROR("Fit failed: no convergence");
    return false;
  }

  MnHesse hesse(2);
  hesse(fcn, fmin, 10000);

  if (!fmin.IsValid())
  {
    ERROR("Fit failed: error in computation of covariance matrix");
    return false;
  }

  const double oldTheta = ad.fPar[eTheta];

  // get best parameters
  ad.fPar[eTheta] = fmin.UserParameters().Value(0);
  ad.fPar[ePhi] = fmin.UserParameters().Value(1);
  ad.fPar[eDistance] = fmin.UserParameters().Value(2)*km;
  ad.fPar[eCTime] = fmin.UserParameters().Value(3)*km;

  const size_t nVar = pars.VariableParameters();

  // get covariance
  ad.fCov = 0;
  for (size_t i = 0; i < nVar; ++i)
    for (size_t j = 0; j < nVar; ++j)
  {
    if (std::isnan(fmin.UserCovariance()(i, j)))
      return false;

    const size_t k = m.ExtOfInt(i);
    const size_t l = m.ExtOfInt(j);
    ad.fCov(k, l) = fmin.UserCovariance()(i, j);
    if (k == eDistance || k == eCTime)
      ad.fCov(k, l) *= km;
    if (l == eDistance || l == eCTime)
      ad.fCov(k, l) *= km;
  }

  {
    ostringstream msg;
    msg << fixed << setprecision(2)
        << "Found minimum after " << fmin.NFcn() << " evaluations\n"
        << "  theta    = " << ad.fPar[eTheta]/degree
                << " +/- " << ad.fCov.Std(eTheta)/degree << " deg (bary)\n"
        << "  phi      = " << ad.fPar[ePhi]/degree
                << " +/- " << ad.fCov.Std(ePhi)/degree << " deg (bary)\n"
        << "  distance = " << ad.fPar[eDistance]/km
                << " +/- " << ad.fCov.Std(eDistance)/km << " km (bary)\n"
        << "  ctime    = " << ad.fPar[eCTime]/km
                << " +/- " << ad.fCov.Std(eCTime)/km << " km (bary)";
    INFO(msg);
  }

  if (fMaxThetaDifference > 0 && (abs(oldTheta - ad.fPar[eTheta]) > fMaxThetaDifference))
  {
    ostringstream msg;
    msg << fixed << setprecision(2)
        << "Fit rejected: theta changed by " << abs(oldTheta - ad.fPar[eTheta])/degree << " deg, [limit: "
        << fMaxThetaDifference/degree << " deg]";
    ERROR(msg);
    return false;
  }

  const double newChi2 = fmin.Fval();

  // do not allow new chi2 to be larger than old chi2 + 1 sigma + 1
  if (ad.fChi2 > 1e-3 && (newChi2 > (ad.fChi2 + sqrt(2*max(ad.fChi2,newChi2)) + 1)))
  {
    ostringstream msg;
    msg << fixed << setprecision(2)
        << "Fit rejected: new chi^2 = " << newChi2 << ", old chi^2 = " << ad.fChi2;
    ERROR(msg);
    return false;
  }

  ad.fChi2 = newChi2;

  return true;
}


bool
SdHorizontalReconstruction::FitShowerSize(SizeData& sd,
                                          const AxisData& ad,
                                          const ExternalGeometryData& gd,
                                          const StationList& list,
                                          const SilentStationList& slist,
                                          const int option)
  const
{
  using namespace ROOT::Minuit2;

  MnUserParameters pars;

  pars.Add("n19", sd.fPar[eN19], 1.0);               // 0
  pars.Add("coreX", sd.fPar[eCoreX]/km, 0.1);        // 1
  pars.Add("coreY", sd.fPar[eCoreY]/km, 0.1);        // 2
  pars.Add("theta", ad.fPar[eTheta], 0.01);          // 3
  pars.Add("phi", ad.fPar[ePhi], 0.01);              // 4
  pars.Add("distance", ad.fPar[eDistance]/km, 0.1);  // 5

  pars.SetLowerLimit("n19", 0);
  pars.SetLimits("theta", 0, 89*degree);
  pars.SetLowerLimit("distance", 0);

  sd.fNFree = 3;
  if (option & eFixCore)
  {
    pars.Fix("coreX");
    pars.Fix("coreY");
    sd.fNFree -= 2;
  }

  if (option & eFixAxis || ad.fCov[eTheta] < 1e-8 || ad.fCov[ePhi] < 1e-8)
  {
    pars.Fix("theta");
    pars.Fix("phi");
  }

  if (option & eFixAxis || ad.fCov[eDistance] < 1e-8)
    pars.Fix("distance");

  ShowerSizeFunction fcn(*this, list, slist, ad, gd);
  if (option & eApproximate)
    fcn.SetType(ShowerSizeFunction::eApprox);

  MnMinimize m(fcn, pars, 1);

  FunctionMinimum fmin = m(10000);

  if (fmin.IsAboveMaxEdm() || fmin.HasReachedCallLimit() || !fmin.IsValid())
  {
    ERROR("Fit failed: no convergence");
    return false;
  }

  MnHesse hesse(2);
  hesse(fcn, fmin, 10000);

  if (!fmin.IsValid())
  {
    ERROR("Fit failed: error in computation of covariance matrix");
    return false;
  }

  sd.fNLogLike = fmin.Fval();

  // get best parameters
  sd.fPar[eN19] = fmin.UserParameters().Value(0);
  sd.fPar[eCoreX] = fmin.UserParameters().Value(1)*km;
  sd.fPar[eCoreY] = fmin.UserParameters().Value(2)*km;

  const size_t nVar = pars.VariableParameters();

  // get covariance
  sd.fCov = 0;
  for (size_t i=0;i<nVar;++i)
    for (size_t j=0;j<nVar;++j)
  {
    if (std::isnan(fmin.UserCovariance()(i, j)))
      return false;

    const size_t k = m.ExtOfInt(i);
    const size_t l = m.ExtOfInt(j);
    if (k < 3 && l < 3)
    {
      sd.fCov(k, l) = fmin.UserCovariance()(i, j);
      if (k == eCoreX || k == eCoreY)
        sd.fCov(k, l) *= km;
      if (l == eCoreX || l == eCoreY)
        sd.fCov(k, l) *= km;
    }
  }

  {
    ostringstream msg;
    msg << fixed << setprecision(2)
        << "Found minimum after " << fmin.NFcn() << " evaluations\n"
        << "  N19   = " << sd.fPar[eN19] <<
               " +/- " << sd.fCov.Std(eN19) << "\n"
        << "  coreX = " << sd.fPar[eCoreX]/km <<
               " +/- " << sd.fCov.Std(eCoreX)/km << " km (bary)\n"
        << "  coreY = " << sd.fPar[eCoreY]/km <<
               " +/- " << sd.fCov.Std(eCoreY)/km << " km (bary)";
    INFO(msg);
  }

  if (fMaxDistanceCoreToBarycenter > 0)
  {
    Point core, origin;
    GetShowerAxis(core, origin, sd, ad);
    Point bary(0, 0, 0, fBaryCS);

    // check for unreasonable core
    const double d = (core - bary).GetMag();
    if (d > fMaxDistanceCoreToBarycenter) {
      ostringstream err;
      err << "Fit rejected: Distance between core and barycenter " << d/kilometer
          << " km too large [limit: " << fMaxDistanceCoreToBarycenter/km << " km]";
      ERROR(err);
      return false;
    }
  }

  return true;
}


void
SdHorizontalReconstruction::PrepareStationData(StationList& list,
                                               SilentStationList& slist,
                                               const Point& core,
                                               const Point& origin,
                                               const SEvent& sEvent)
  const
{
  const sdet::SDetector& sDetector =
    det::Detector::GetInstance().GetSDetector();

  for (SEvent::ConstCandidateStationIterator
         sIt = sEvent.CandidateStationsBegin(),
         sEnd = sEvent.CandidateStationsEnd();
       sIt != sEnd; ++sIt)
  {
    list.push_back(StationData());
    StationData& sd = list.back();
    const StationRecData& sRec = sIt->GetRecData();

    sd.fId        = sIt->GetId();
    sd.fPos       = sDetector.GetStation(*sIt).GetPosition();
    sd.fSignal    = sRec.GetTotalSignal();
    sd.fCTime     = kSpeedOfLight *
      (sRec.GetSignalStartTime() - fBaryTime).GetInterval();
    sd.fT50       = sRec.GetT50();
    sd.fSaturated = sIt->IsLowGainSaturation();
    sd.fRecoveryErr = sRec.GetRecoveredSignalError();
    if (sd.fSaturated && sd.fRecoveryErr > 0)
      sd.fSignal = sRec.GetRecoveredSignal();
    sd.fRejected = !sIt->IsCandidate();
  }

  for (SEvent::ConstSilentStationIterator
         sIt = sEvent.SilentStationsBegin(),
         sEnd = sEvent.SilentStationsEnd();
       sIt != sEnd; ++sIt)
  {
    const Point& pos = sDetector.GetStation(*sIt).GetPosition();

    if (GetRho(pos, core, origin) > fMaxSilentRadius)
      continue;

    for (StationList::const_iterator it = list.begin();
         it != list.end(); ++it)
    {
      if ((pos - it->fPos).GetMag() < fMaxSilentRelDistance)
      {
        slist.push_back(SilentStationData());
        SilentStationData& sd = slist.back();
        sd.fId  = sIt->GetId();
        sd.fPos = pos;
        // finding one candidate station nearby is enough
        break;
      }
    }
  }
}


void
SdHorizontalReconstruction::SetReconstructedValues(Event& event,
                                                      const SizeData& sizeData,
                                                      const AxisData& axisData,
                                                      const StationList& list,
                                                      const SilentStationList& slist,
                                                      const ExternalGeometryData& extGeomData)
  const
{
  // Convert N19 to energy
  const double energy =
    fEnergyConstant * pow(sizeData.fPar[eN19],fGamma);
  const double energyErr =
    fGamma * sizeData.fCov.Std(eN19) / sizeData.fPar[eN19] * energy;

  {
    ostringstream msg;

    msg << "Result:\n"
        << fixed << setprecision(2)
        << "  N19 (corrected) : "
        << sizeData.fPar[eN19] << " +/- " << sizeData.fCov.Std(eN19) << "\n"
        << "  Energy / EeV    : "
        << energy/EeV << " +/- " << energyErr/EeV;

    if (event.HasSimShower()) {
      const double energyMC = event.GetSimShower().GetEnergy();
      const double diffToMC = energy/energyMC-1.;
      const double diffToMCSigma =
        (energy-energyMC)/energyErr;
      msg << "\n  Energy (MC-input) / EeV : "
          << energyMC/EeV
          << " ( " << diffToMC*100.
          << " % | " << diffToMCSigma << " std.dev. )";
    }

    INFO(msg);
  }

  const bool haveCurvature = axisData.fCov[eDistance] > 0;

  Point core, origin;
  GetShowerAxis(core, origin, sizeData, axisData);
  Vector axis = origin - core;
  axis.Normalize();

  vector<double> oPar;
  CovarianceMatrix oCov;
  const size_t nAxis = 4;
  { // propagate core uncertainty into shower axis parameters and station radii and
    // propagate axis parameters from bary CS to core CS
    // Beware: eCTime is time *at core* (not at origin) after propagation
    vector<double> iPar(6);
    iPar[0] = sizeData.fPar[eCoreX];
    iPar[1] = sizeData.fPar[eCoreY];
    iPar[2] = axisData.fPar[eTheta];
    iPar[3] = axisData.fPar[ePhi];
    iPar[4] = axisData.fPar[eDistance];
    iPar[5] = axisData.fPar[eCTime];
    CovarianceMatrix iCov(6);
    iCov[0]   = sizeData.fCov[eCoreX];
    iCov[1]   = sizeData.fCov[eCoreY];
    iCov(0,1) = sizeData.fCov(eCoreX, eCoreY);
    for (size_t i = 0; i < AxisData::size; ++i)
      for (size_t j = i; j < AxisData::size; ++j)
        iCov(2+i, 2+j) = axisData.fCov(i, j);

    vector<Point> posList;
    for (StationList::const_iterator
           sIt = list.begin(),
           sEnd = list.end();
         sIt != sEnd; ++sIt)
      posList.push_back(sIt->fPos);

    oPar.resize(nAxis + posList.size());
    oCov.SetExtent(nAxis + posList.size());

    FinalPropagator prop(*this, posList);
    prop(oPar, oCov, iPar, iCov);

    // if parameter was fixed to begin with
    // the error propagation may not assign a non-zero uncertainty
    for (size_t i = 0; i < 4; ++i)
      for (size_t j = i; j < 4; ++j)
        if (iCov(2+i,2+j) == 0) oCov(i,j) = 0;
  }

  const double n19 = sizeData.fPar[eN19];

  const CoordinateSystemPtr coreCS =
    LocalCoordinateSystem::Create(core);

  const double theta = axis.GetTheta(coreCS);
  const double phi = axis.GetPhi(coreCS);

  const CoordinateSystemPtr showerPlaneCS =
    GetShowerCoordinateSystem(theta, phi, coreCS);

  double timeSum = 0;
  double energyChi2 = 0;
  double axisChi2 = 0;
  int nUsed = 0;
  // calculate signal and time residuals and corresponding chi2s
  ShowerFrontFunction frontFcn(*this, core, list, extGeomData);
  ShowerSizeFunction sizeFcn(*this, list, slist, axisData, extGeomData);
  for (size_t i = 0; i < list.size(); ++i)
  {
    const StationData& sd = list[i];
    StationRecData& sRec = event.GetSEvent().GetStation(sd.fId).GetRecData();

    sRec.SetAzimuthShowerPlane(sd.fPos.GetPhi(showerPlaneCS));
    sRec.SetSPDistance(oPar[nAxis + i], oCov.Std(nAxis + i));

    double signalPred = 0, signalSigma = 0;
    sizeFcn.Predict(signalPred, signalSigma, sd.fPos, origin, n19, core, sd.fRecoveryErr);

    if (signalSigma == 0)
    {
      ostringstream msg;
      msg << fixed << setprecision(2)
          << "Station " << sd.fId << ":"
          << "signal = " << sd.fSignal << " "
          << "pred = " << signalPred << " "
          << "sigma = " << signalSigma;
      WARNING(msg);
    }
    else
    {
      const double residual = (sd.fSignal - signalPred)/signalSigma;
      // sd.fSignal may be recovered, set original signal here (like in LDFFinder)
      sRec.SetTotalSignal(sRec.GetTotalSignal(), signalSigma);
      sRec.SetLDFResidual(residual);
      if (!sd.fRejected)
        energyChi2 += Sqr(residual);
    }

    if (axisData.fNFree > 0) {
      double meanCT, sigmaCT;
      frontFcn.Predict(meanCT, sigmaCT, sd, origin, axisData.fPar[eCTime]);
      const double timeResidual = (sd.fCTime - meanCT)/sigmaCT;
      sRec.SetResidual(timeResidual);
      if (!sd.fRejected)
      {
        axisChi2 += Sqr(timeResidual);
        timeSum += timeResidual;
      }
    }

    if (!sd.fRejected)
      ++nUsed;
  }

  ShowerSRecData& recShower = event.GetRecShower().GetSRecShower();

  // make average 1D LDF for the lack of something better
  if (!recShower.HasLDF())
    recShower.MakeLDF();

  {
    TabulatedFunctionErrors& tabulatedLDF = recShower.GetLDF();
    tabulatedLDF.Clear();

    const double x0 = log(1);
    const double x1 = log(10000);

    const size_t nLDFPoints = 20;
    for (size_t i = 0; i < nLDFPoints; ++i)
    {
      const double z = float(i)/(nLDFPoints-1);
      double mean = 0, smin = 0, smax = 0;
      const size_t nPsi = 4;
      const double rho0 = exp((1.0-z)*x0 + z*x1);
      for (size_t j = 0; j < nPsi; ++j) {
        const double psi = 2*kPi*double(j)/(nPsi-1);
        Point pos(rho0, psi, 0, showerPlaneCS, Point::kCylindrical);
        pos -= pos.GetZ(coreCS)/axis.GetZ(coreCS)*axis;

        double signalPred = 0, signalSigma = 0;
        sizeFcn.Predict(signalPred, signalSigma, pos, origin, n19, core);

        mean += signalPred/nPsi;
        if (signalPred < smin || smin == 0) smin = signalPred;
        if (signalPred > smax || smax == 0) smax = signalPred;
      }
      tabulatedLDF.PushBack(rho0, 0, mean, 0.5*(smax-smin));
    }
  }

  if (axisData.fNFree > 0)
  {
    // oPar[eCTime] is time at core, not origin time!
    const TimeStamp finalCoreTime = fBaryTime + oPar[eCTime]/kSpeedOfLight;
    const double finalCoreTimeErr = oCov.Std(eCTime)/kSpeedOfLight;
    recShower.SetAxis(axis);
    recShower.SetCoreTime(finalCoreTime, finalCoreTimeErr);
    recShower.SetAngleChi2(axisChi2, max(nUsed - axisData.fNFree, 0));
    const double timeMean = timeSum / nUsed;
    const double timeSpread = sqrt(axisChi2/(nUsed - 1));
    recShower.SetTimeResidualMean(timeMean);
    recShower.SetTimeResidualSpread(timeSpread);

    recShower.SetThetaError(oCov.Std(eTheta));
    recShower.SetPhiError(oCov.Std(ePhi));
    recShower.SetCorrelationThetaPhi( oCov(eTheta,ePhi)/sqrt(oCov[eTheta]*oCov[ePhi]) );

    const double rc = oPar[eDistance];
    if (haveCurvature)
    {
      const double rcErr = oCov.Std(eDistance);
      recShower.SetCurvature(1.0/rc, rcErr/Sqr(rc));
    }
    else
      recShower.SetCurvature(1.0/rc, 0);
  }
  else
  {
    // core moved but have axis from plane fit: shift core time
    const TimeInterval cdiff =
      recShower.GetAxis() * (recShower.GetCorePosition() - core);
    const TimeStamp newCoreTime =
      recShower.GetCoreTime() + cdiff/kSpeedOfLight;
    recShower.SetCurvature(0.0, 0.0);
    recShower.SetCoreTime(newCoreTime, recShower.GetCoreTimeError());
  }

  const Vector coreErr(sizeData.fCov.Std(eCoreX),
                       sizeData.fCov.Std(eCoreY),
                       0, coreCS);

  recShower.SetCorePosition (core);
  recShower.SetCoreError    (coreErr);
  recShower.SetCorrelationXY(sizeData.fCov(eCoreX, eCoreY)/
                             sqrt(sizeData.fCov[eCoreX]*sizeData.fCov[eCoreY]));
  recShower.SetEnergy       (energy, energyErr);
  recShower.SetShowerSizeType(ShowerSRecData::eN19);
  recShower.SetShowerSize   (sizeData.fPar[eN19], sizeData.fCov.Std(eN19));
  recShower.SetLDFChi2      (energyChi2, max(nUsed - sizeData.fNFree, 0));
  recShower.SetLDFLikelihood(sizeData.fNLogLike);

  recShower.SetLDFRecStage  (sizeData.fRecStage);
}


void
SdHorizontalReconstruction::GetShowerAxis(Point& core,
                                          Point& origin,
                                          const double coreX,
                                          const double coreY,
                                          const double theta,
                                          const double phi,
                                          const double distance)
  const
{
  core = Point(coreX, coreY, 0, fBaryCS);
  origin = Point(distance, theta, phi, fBaryCS, Point::kSpherical);
}


void
SdHorizontalReconstruction::GetShowerAxis(Point& core,
                                          Point& origin,
                                          const SizeData& sd,
                                          const AxisData& ad)
  const
{
  GetShowerAxis(core, origin, sd.fPar[eCoreX], sd.fPar[eCoreY],
                ad.fPar[eTheta], ad.fPar[ePhi], ad.fPar[eDistance]);
}


void
SdHorizontalReconstruction::LocalCoordinates(double& xProjected,
                                             double& yProjected,
                                             double& rho, double& theta,
                                             const utl::Point& pos,
                                             const utl::Point& origin,
                                             const utl::Point& core)
  const
{
  using namespace utl;

  const Vector posToOrigin = origin - pos;
  const double d = pos.GetZ(fBaryCS)/posToOrigin.GetZ(fBaryCS);
  const Vector corePlanePos = pos - core - d*posToOrigin;
  // good approximation: pretend that baryCS has same orientation as coreCS
  xProjected = corePlanePos.GetX(fBaryCS);
  yProjected = corePlanePos.GetY(fBaryCS);
  rho = GetRho(pos, core, origin);
  theta = Angle(posToOrigin, (pos - fEarthCenter));
}


double
SdHorizontalReconstruction::ExpectedCurvatureRadius(const double theta)
  const
{
  /* Empirical parameterization obtained by fitting real events,
     valid from 50 deg to 90 deg */
  return exp(fCurvModelPars[0] + fCurvModelPars[1] * pow(theta, fCurvModelPars[2]))*utl::km;
}