ScintillatorLDFFinder.cc

Go to the documentation of this file.

#include <evt/Event.h>

#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerScintillatorRecData.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 <sevt/StationConstants.h>
#include <sevt/Scintillator.h>
#include <sevt/ScintillatorRecData.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/SDetectorConstants.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 <fwk/RunController.h>

#include <utl/CoordinateSystem.h>
#include <utl/PhysicalConstants.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 <utl/NumericalErrorPropagator.h>

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

#include "ScintillatorLDFFinder.h"

#include <boost/range/iterator_range.hpp>
#include <cmath>

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

using CLHEP::HepRotation;

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


namespace ScintillatorLDFFinderKG {

  // local helpers

  template<typename T>
  inline
  bool
  HasNaN(const T& result)
  {
    const size_t npar = result.fPar.size();
    for (size_t i = 0; i < npar; ++i) {
      if (std::isnan(result.fPar[i]))
        return true;
      for (size_t j = i; j < npar; ++j)
        if (std::isnan(result.fCov(i, j)))
          return true;
    }
    return false;
  }


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


  inline
  string
  GetShowerSizeLabel(const LDFFitConfig& config)
  {
    ostringstream label;
    const int refdist = static_cast<int>(config.fLDF.GetReferenceDistance());
    label << "S" << refdist/meter;
    if (refdist < 1000*meter)
      label << ' ';
    return label.str();
  }


  inline
  const char*
  GetShapeLabel(const int iShape)
  {
    static const char* const labels[] = { "beta ", "gamma", "mu   ", "tau  " };
    return labels[iShape];
  }


    /* 
    This is the calculation of the core uncertainties from WCD measurement, 
    parametrized by Quentin L. 
    It probably shouldn't stay here, but moved to 
    a more appropriate place. 
  */ 
  inline
  double
  DeltaR(const double wcdShowerSize, const double cosTheta)
  { 
    const double f0 = 112.71;
    const double f1 = 30.85;
    const double f2 = -58.08;
    const double f3 = -12.43;
    const double f4 = 10.28;
    const double f5 = 11.96;

    const double sinTheta2 = 1 - cosTheta * cosTheta;

    const double a = f0 + f1 * sinTheta2;
    const double b = f2 + f3 * sinTheta2;
    const double c = f4 + f5 * sinTheta2;

    const double deltaR = a + b * std::log10(wcdShowerSize) 
      + c * std::log10(wcdShowerSize)* std::log10(wcdShowerSize);

    return deltaR;
  }

}




void
ScintillatorLDFFinder::OutputResults(const Event& event)
  const
{
  const ShowerSRecData& recShower = event.GetRecShower().GetSRecShower();
  const ShowerScintillatorRecData& scinRecShower = recShower.GetScintillatorRecShower();
  OUT(eFinal)
    << "***** Final reconstructed shower data *****\n"
       "  LDF stage             = " << " --- " << "\n"
       "  " << GetShowerSizeLabel(fLDFFitConfig) << "                 = "
            << scinRecShower.GetShowerSize()
            << " +/- " << scinRecShower.GetShowerSizeError() << "\n"
       "  beta                  = " << scinRecShower.GetBeta() << " +/- " << scinRecShower.GetBetaError() << "\n"
       "  gamma                 = " << scinRecShower.GetGamma() << " +/- " << scinRecShower.GetGammaError() << "\n"
       "  LDF chi2 / ndof       = " << scinRecShower.GetLDFChi2() << " / " << scinRecShower.GetLDFNdof() << "\n"
    << endl;
}


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

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

  topB.GetChild("infoLevel").GetData(fInfoLevel); // for verbosity

  {
    const string sigVariance = topB.GetChild("signalVariance").Get<string>();
    if (sigVariance == "eGAP2018_025")
      fLDFFitConfig.fSignalVarianceModel = ssd::eGAP2018_025;
    else if (sigVariance == "eGAP2019_000")
      fLDFFitConfig.fSignalVarianceModel = ssd::eGAP2019_000;
    else {
      ERROR("Unknown signal variance");
      return eFailure;
    }

    // setup LDF
    string ldfType;
    topB.GetChild("ldfType").GetData(ldfType);

    double ldfReferenceDistance;
    topB.GetChild("ldfReferenceDistance").GetData(ldfReferenceDistance);

    vector<double> ldfParameters;
    topB.GetChild("ldfParameters").GetData(ldfParameters);

    vector<double> ldfUncertaintyParameters;
    topB.GetChild("ldfUncertaintyParameters").GetData(ldfUncertaintyParameters);

    fLDFFitConfig.fLDF = LDF(ldfType, ldfReferenceDistance,
                             ldfParameters, ldfUncertaintyParameters);
    fLDFFitConfig.fLDFShapeTreatment.resize(fLDFFitConfig.fLDF.GetNShapeParameters(), eFitted);
  }

  topB.GetChild("useReconstructedCore").GetData(fFixedCore);
  fLDFFitConfig.fFixCore = fFixedCore;

  topB.GetChild("propagateCoreUncertainty").GetData(fCorePropagation); 
  fLDFFitConfig.fCorePropagationEnabled = fCorePropagation;

  string minimizationMethod;
  topB.GetChild("minimizationMethod").GetData(minimizationMethod);
  fUseMaxLikelihood = (minimizationMethod == "MaxLike");

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

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

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

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

  // Keep in mind to implement stage where beta and gamma are relaxed
  const Branch stage2B = topB.GetChild("stage2");
  stage2B.GetChild("minNumberRelaxBeta").GetData(fStage2.fMinNumberRelaxBeta);
  stage2B.GetChild("minNumberRelaxGamma").GetData(fStage2.fMinNumberRelaxGamma);

  {
    ostringstream info;
    info << "\n"
         << "             LDF type: " << fLDFFitConfig.fLDF.GetType() << "\n"
         << "            Core type: " << (fFixedCore ? "Reconstructed" : "Fitted") << "\n"
         << "   Core uncertainties: " << (fCorePropagation ? "Propagated" : "Not included") << "\n"
         << "Shower size estimator: " << GetShowerSizeLabel(fLDFFitConfig) << "\n"
         << "  Minimization method: " << (fUseMaxLikelihood ? "Maximum-Likelihood" : "Chi2") << "\n";
    INFO(info);
  }

  return eSuccess;
}


VModule::ResultFlag
ScintillatorLDFFinder::Run(evt::Event& event)
{
  // nothing to do if there is no SEvent
  if (!(event.HasSEvent() && event.HasRecShower() &&
        event.GetRecShower().HasSRecShower()))
    return eContinueLoop;

  const SEvent& sEvent = event.GetSEvent();

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

  if (!showerRec.HasLDF()) {
    OUT(eFinal)
      << "No LDF found from SD reconstruction!\n";
    return eContinueLoop;
  }

  // barycenter is done by the SdPlaneFitOG module
  fBarycenter = showerRec.GetBarycenter();
  fBaryTime = showerRec.GetBarytime();
  gBaryCS = showerRec.GetBarycenterCoordinateSystem();

  // load variables from SD reconstruction
  Point core = showerRec.GetCorePosition();
  Vector showerAxis = showerRec.GetAxis();
  const double cosTheta = showerAxis.GetCosTheta(gBaryCS);
  const double wcdShowerSize = showerRec.GetShowerSize();
  
  if (!fBaryTime)
    throw utl::NonExistentComponentException("needs barycenter!");

  OUT(eIntermediate)
  << "Distance cut for SSDs = " << DistanceCut(wcdShowerSize, cosTheta) << " m\n";

  int nStations = 0;
  {
    int nSaturated = 0;
    int nSilent = 0;
    int nDistant = 0;
    for (const auto& station : boost::make_iterator_range(sEvent.StationsBegin(), sEvent.StationsEnd()))
      if (station.IsCandidate() && station.HasScintillator() && station.GetScintillator().HasRecData()) {
        ++nStations;
        
        const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
        const utl::Point pos = sDetector.GetStation(station).GetPosition();
        const double rho = RPerp(showerAxis, pos - core);
        if (rho > DistanceCut(wcdShowerSize, cosTheta))
          ++nDistant;

        if (station.GetScintillator().IsLowGainSaturation())
          ++nSaturated;

      } else if (station.IsSilent())
        ++nSilent;

    // right now saturated SSDs are not handled in the likelihood
    nStations -= (nDistant + nSaturated);

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

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

  fLDFFitConfig.deltaR = DeltaR(wcdShowerSize,cosTheta);  // Calculate DeltaR, should be moved
  LDFFitConfig config = fLDFFitConfig; // from XML
  LDFFitResult best(config);

  best.fPar[eCoreX] = core.GetX(gBaryCS);
  best.fPar[eCoreY] = core.GetY(gBaryCS);

  /* 
     showersize, shower axis and core position 
     from SD reconstruction are needed to apply distance cut
     on SSD signals.

  */ 
  vector<StationFitData> stationFitData = MakeStationFitData(sEvent, core, wcdShowerSize, showerAxis);

  // --- Stage: fit shower size, LDF shape & core fixed --------------------

  best.fRecStage = ShowerSRecData::kLDFEstimate;
  OUT(eIntermediate)
    << "Stage " << best.fRecStage << ": fit " << GetShowerSizeLabel(config)
    << " - LDF shape fixed\n";

  FitLDFSimplified(best, config, showerAxis, stationFitData);

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

  // --- Stage: fit shower size, shape fixed --------------
  {
    LDFFitResult result = best;

    config.fLDFShapeTreatment[0] = eModeled;
    config.fLDFShapeTreatment[1] = eModeled;

    OUT(eIntermediate)
      << "Stage " << result.fRecStage << ": fit "
      << GetShowerSizeLabel(config)
      << " - core and shape fixed \n";

    if (!FitLDF(result, config, showerAxis, stationFitData)){
      OUT(eFinal)
        << "Fit of the LDF FAILED! \n";
      return eContinueLoop; // check this
    }

    result.fRecStage = ShowerSRecData::kLDF;
    best = result;

  } // stage

  // --- Stage: fit shower size, shape fitting (optional)  --------------
  {
    LDFFitResult result = best;

    config.fLDFShapeTreatment[0] =
      (nStations < fStage2.fMinNumberRelaxBeta ||
       FixBeta(result, showerAxis, stationFitData)) ? eModeled : eFitted;
    config.fLDFShapeTreatment[1] =
      (nStations < fStage2.fMinNumberRelaxGamma ||
       FixGamma(result, showerAxis, stationFitData)) ? eModeled : eFitted;

    if (config.fLDFShapeTreatment[0] == eFitted ||
        config.fLDFShapeTreatment[1] == eFitted) {
      OUT(eIntermediate)
        << "Stage " << result.fRecStage << ": fit "
        << GetShowerSizeLabel(config)
        << " - core fixed"
        << (config.fLDFShapeTreatment[0] == eFitted ? ", beta" : "")
        << (config.fLDFShapeTreatment[1] == eFitted ? ", gamma\n" : "\n");      

      if (FitLDF(result, config, showerAxis, stationFitData)) {
        const bool fixBeta = FixBeta(result, showerAxis, stationFitData);
        const bool fixGamma = FixGamma(result, showerAxis, stationFitData);
        if ((config.fLDFShapeTreatment[0] == eModeled) == fixBeta &&
            (config.fLDFShapeTreatment[1] == eModeled) == fixGamma) {
          result.fRecStage = ShowerSRecData::kLDF +
            (config.fLDFShapeTreatment[0] == eFitted) * ShowerSRecData::kLDFBeta +
            (config.fLDFShapeTreatment[0] == eFitted) * ShowerSRecData::kLDFGamma;
          best = result;
        } else {
          if (fixBeta)
            config.fLDFShapeTreatment[0] = eModeled;
          if (fixGamma)
            config.fLDFShapeTreatment[1] = eModeled;

          if (FitLDF(result, config, showerAxis, stationFitData)) {
            result.fRecStage = ShowerSRecData::kLDF +
              (config.fLDFShapeTreatment[0] == eFitted) * ShowerSRecData::kLDFBeta +
              (config.fLDFShapeTreatment[1] == eFitted) * ShowerSRecData::kLDFGamma;
            best = result;
          }
        }
      }
    }
  } // stage

  SetRecData(event, best);
  OutputResults(event);

  return eSuccess;
}

// 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
//
// - Adjustments for 750 m infill:
// - Interval: (250, 750); distances: 250, 200, 150
//
// New version much stricter
bool
ScintillatorLDFFinder::FixBeta(const LDFFitResult& result,
                   const Vector& showerAxis,
                   const vector<StationFitData>& data)
  const
{
  size_t nStations = 0;
  for (const auto& station : data) {
    if (!station.fSignal)
      continue;
    ++nStations;
  }

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

  int nStationsInRange = 0;
  double maxDistance = 0;

  double rRange[2];
  double leverArms[3];
  double refDist(fLDFFitConfig.fLDF.GetReferenceDistance());

  if (refDist == 1000*meter) {
    rRange[0] = 400*meter;
    rRange[1] = 1600*meter;
    leverArms[0] = 900*meter;
    leverArms[1] = 800*meter;
    leverArms[2] = 700*meter;
  } else if (refDist == 450*meter) {
    rRange[0] = 180*meter;
    rRange[1] = 720*meter;
    leverArms[0] = 405*meter;
    leverArms[1] = 360*meter;
    leverArms[2] = 315*meter;
  } else if (refDist == 250*meter) {
    rRange[0] = 100*meter;
    rRange[1] = 400*meter;
    leverArms[0] = 225*meter;
    leverArms[1] = 200*meter;
    leverArms[2] = 175*meter;
  } else // unknown case
    return false;

  const Point core = result.GetCore();

  for (const auto& station : data) {

    if (!station.fSignal)
      continue;

    const double r = RPerp(showerAxis, station.fPos - core);

    if (rRange[0] < r && r < rRange[1]) {
      ++nStationsInRange;

      for (const auto& s : data) {
        const double rr = RPerp(showerAxis, s.fPos - core);

        if (rRange[0] < rr && rr < rRange[1]) {
          const double distance = fabs(r - rr);
          if (distance > maxDistance)
            maxDistance = distance;
        }
      }
    }
  }

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

  return fixBeta;
}


bool
ScintillatorLDFFinder::FixGamma(const LDFFitResult& result,
                    const Vector& showerAxis,
                    const vector<StationFitData>& data)
  const
{
  size_t nStations = 0;
  for (const auto& station : data) {
    if (!station.fSignal)
      continue;
    ++nStations;
  }

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

  int nStationsInRange = 0;
  double maxDistance = 0;

  double rRange[2];
  double leverArms[3];
  double refDist(fLDFFitConfig.fLDF.GetReferenceDistance());

  if (refDist == 1000*meter) {
    rRange[0] = 1000*meter;
    rRange[1] = 2000*meter;
    leverArms[0] = 500*meter;
    leverArms[1] = 400*meter;
    leverArms[2] = 300*meter;
  } else if (refDist == 450*meter) {
    rRange[0] = 450*meter;
    rRange[1] = 1500*meter;
    leverArms[0] = 250*meter;
    leverArms[1] = 200*meter;
    leverArms[2] = 150*meter;
  } else { // unknown case
    return false;
  }

  const Point core = result.GetCore();

  for (const auto& station : data) {
    if (!station.fSignal)
      continue;

    const double r = RPerp(showerAxis, station.fPos - core);

    if (rRange[0] < r && r < rRange[1]) {

      ++nStationsInRange;

      for (const auto& s : data) {

        const double rr = RPerp(showerAxis, s.fPos - core);

        if (rRange[0] < rr && rr < rRange[1]) {
          const double distance = fabs(r - rr);
          if (distance > maxDistance)
            maxDistance = distance;
        }
      }
    }
  }

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

  return fixGamma;
}

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


void
ScintillatorLDFFinder::SetRecData(Event& event, const LDFFitResult& lres)
  const
{
  ShowerSRecData& recShower = event.GetRecShower().GetSRecShower();

  if (!recShower.HasScintillatorRecShower())
    recShower.MakeScintillatorRecShower();

  ShowerScintillatorRecData& scinRecShower = recShower.GetScintillatorRecShower();

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

  const LDF& ldf = fLDFFitConfig.fLDF;
  const double showerSize = lres.fPar[eShowerSize];
  const vector<double>& ldfShapePars = lres.GetLDFShapeParameters();
  const Point& core = lres.GetCore();
  const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(core);
  const Vector& axis = recShower.GetAxis();
  const double cosTheta = axis.GetCosTheta(coreCS);

  SEvent& sEvent = event.GetSEvent();

  const double theta = axis.GetTheta(coreCS);
  const double phi = axis.GetPhi(coreCS);
  const CoordinateSystemPtr showerPlaneCS(
    CoordinateSystem::RotationY(
      theta,
      CoordinateSystem::RotationZ(
        phi,
        coreCS
      )
    )
  );

  Accumulator::MinMax<double> minMaxRho(1000*kilometer, 0);
  //Accumulator::SampleStandardDeviation timeStat;// unused. LN.

  for (auto& station : boost::make_iterator_range(sEvent.StationsBegin(), sEvent.StationsEnd())) {
    if (!station.HasRecData())
      continue;

    if (!station.HasScintillator())
      continue;

    Scintillator& scintillator = station.GetScintillator();

    if (!scintillator.HasRecData())
      continue;

    const Point& sPos = sDetector.GetStation(station).GetPosition();
    const double x = sPos.GetX(showerPlaneCS);
    const double y = sPos.GetY(showerPlaneCS);
    const double rho = sqrt(x*x + y*y);

    if (station.IsCandidate())
      minMaxRho(rho);

    ScintillatorRecData& scinRecData = scintillator.GetRecData();
    const double signal = scinRecData.GetTotalSignal();
    const double funcval = showerSize * ldf(rho, ldfShapePars);
    double sigma = utl::ssd::SignalUncertainty(fLDFFitConfig.fSignalVarianceModel, cosTheta, funcval);
    scinRecData.SetTotalSignal(signal, sigma);
    scinRecData.SetLDFResidual((signal - funcval)/sigma);

  } // loop over stations with RecData and ScintillatorRecData

  if (!scinRecShower.HasLDF())
    scinRecShower.MakeLDF();
  TabulatedFunctionErrors& tabulatedLDF = scinRecShower.GetLDF();
  tabulatedLDF.Clear();

  {
    const double x0 = log(1);
    const double x1 = log(10000);
    const int nLDFPoints = 20;

    for (int i = 0; i < nLDFPoints; ++i) {
      const double z = double(i) / (nLDFPoints-1);
      const double rho = exp((1-z)*x0 + z*x1);
      const double funcval = showerSize * ldf(rho, ldfShapePars);
      const double sigma = utl::ssd::SignalUncertainty(fLDFFitConfig.fSignalVarianceModel, cosTheta, funcval);
      tabulatedLDF.PushBack(rho, 0, funcval, sigma);
    }
  }

  {
    const double showerSizeErr = lres.fCov.Std(eShowerSize);
    vector<double> ldfShapeParsErr;
    for (size_t i = 0, n = fLDFFitConfig.fLDF.GetNShapeParameters(); i < n; ++i)
      ldfShapeParsErr.push_back(lres.fCov.Std(eShapeParameters + i));

    switch (size_t(round(ldf.GetReferenceDistance()/meter))) {
    case 1000:
      scinRecShower.SetShowerSizeType(ShowerScintillatorRecData::eS1000);
      break;
    default:
      scinRecShower.SetShowerSizeType(ShowerScintillatorRecData::eUndefined);
      break;
    }

    scinRecShower.SetShowerSize(showerSize, showerSizeErr);
    scinRecShower.SetBeta(ldfShapePars[0], ldfShapeParsErr[0]);

    if (ldfShapePars.size() > 1)
      scinRecShower.SetGamma(ldfShapePars[1], ldfShapeParsErr[1]);
    else
      scinRecShower.SetGamma(0, 0);

    scinRecShower.SetLDFChi2(lres.fChi2, lres.fNdof);
    scinRecShower.SetLDFLikelihood(lres.fLogLikelihood);
    scinRecShower.SetLDFRecStage(lres.fRecStage);
  }
}


vector<StationFitData>
ScintillatorLDFFinder::MakeStationFitData(const SEvent& sEvent, const Point core, const double wcdShowerSize, const Vector& showerAxis)
  const
{
  vector<StationFitData> sFitData;

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

  for (const auto& station : boost::make_iterator_range(sEvent.StationsBegin(), sEvent.StationsEnd())) {
    /*
     * For now, only consider scintillators which were WCD
     * candidates. Silent WCD stations do not guarantee that
     * you would have had a zero signal in the scintillator
     * as well. We need to deal with this properly in the likelihood,
     * but here we do not store information from the scintillators for
     * now.
     */ 

    if (!station.IsCandidate())
        continue;

    if (!station.HasScintillator())
        continue;

    if (!station.GetScintillator().HasRecData())
        continue;
    
    const double cosTheta = showerAxis.GetCosTheta(gBaryCS);
    const utl::Point pos = sDetector.GetStation(station).GetPosition();
    const double rho = RPerp(showerAxis, pos - core);
    /*
      Apply distance cut to remove scintillator with low signals:
      For a given shower size and zenith angle, the cut is determined 
      in such a way that >95% of the SSDs below that distance
      have at least 1 MIP of signal.
    */

    if (rho > DistanceCut(wcdShowerSize, cosTheta))
      continue;
          
    StationFitData current;

    const ScintillatorRecData& scinRec = station.GetScintillator().GetRecData();

    current.fId = station.GetId();
    current.fIsSaturated = station.GetScintillator().IsLowGainSaturation();
    current.fPos = pos;
    current.fSignal = scinRec.GetTotalSignal();
    sFitData.push_back(current);

  } // loop over stations

  return sFitData;
}


void
ScintillatorLDFFinder::FitLDFSimplified(LDFFitResult& result,
                                        LDFFitConfig& config,
                                        const Vector& showerAxis,
                                        const vector<StationFitData>& data)
  const
{
  result.fCov = 0;

  const utl::Point core = result.GetCore();

  const double cosTheta = showerAxis.GetCosTheta(gBaryCS);

  const vector<double> shape = config.fLDF.ShapeModel(cosTheta, 10);

  LDFLikelihoodFunctionTN likeFcn(config, showerAxis, data);

  LDFFitResult resultA = result;
  LDFFitResult resultB = result;

  {
    // estimate shower size from station closest to reference distance

    const double kRefDist = config.fLDF.GetReferenceDistance();
    double minDistance = numeric_limits<double>::max();
    double rhoClosest = kRefDist;
    double signalClosest = 1;

    // loop over normal stations, ignore silent/saturated stations
    for (const auto& station : data) {

      if (!station.fSignal)
        continue;

      const double rho = RPerp(showerAxis, station.fPos - core);
      const double distance = fabs(rho - kRefDist);

      if (distance < minDistance) {
        minDistance = distance;
        rhoClosest = rho;
        signalClosest = station.fSignal;
      }
    }

    resultA.fPar[eShowerSize] = signalClosest / config.fLDF(rhoClosest, shape);
    resultA.fLogLikelihood = likeFcn(resultA.fPar);

    OUT(eIntermediate)
      << "  " << GetShowerSizeLabel(config) << " [estimate A]  = "
      << resultA.fPar[eShowerSize] << " ( logLike = " << resultA.fLogLikelihood << " )" << endl;
  }

  {
    /*
      If core and shape of the LDF are fixed, the least-squares
      estimate of shower size S0 can be calculated analytically.

      S_i                ... data
      S(r_i) = S0 f(r_i) ... model
      sigma^2(S) = S     ... simplest signal uncertainty model

      chi^2 = sum_i (S_i - S(r_i))^2/S(r_i)

      With maximum condition d chi^2/d S0 = 0 one obtains
        S0^2 = sum_i S_i^2 f^(-1)(r_i) / sum_i f(r_i)
    */

    double sumA = 0;
    double sumB = 0;
    for (const auto& station : data) {
      const double s = station.fSignal;

      if (!s)
        continue;

      const double rho = RPerp(showerAxis, station.fPos - core);
      const double f = config.fLDF(rho, shape);

      sumA += s*s / f;
      sumB += f;
    }

    resultB.fPar[eShowerSize] = sqrt(sumA / sumB);
    resultB.fLogLikelihood = likeFcn(resultB.fPar);

    OUT(eIntermediate)
      << "  " << GetShowerSizeLabel(config) << " [estimate B]  = "
      << resultB.fPar[eShowerSize] << " ( logLike = " << resultB.fLogLikelihood << " )" << endl;
  }

  // use better estimate of the two
  result = resultA.fLogLikelihood < resultB.fLogLikelihood ? resultA : resultB;

  // Finally, perform full-featured fit with fixed core
  // gamma set to WCD parameterisation
  // fit beta
  LDFFitConfig thisConfig = config;
  thisConfig.fFixCore = true;
  thisConfig.fLDFShapeTreatment[0] = eModeled;
  thisConfig.fLDFShapeTreatment[1] = eModeled;
  FitLDF(result, thisConfig, showerAxis, data);
}


bool
ScintillatorLDFFinder::FitLDF(LDFFitResult& result,
                              LDFFitConfig& config /* needs to be mutable */,
                              const Vector& showerAxis,
                              const vector<StationFitData>& data)
  const
{
  using namespace ROOT::Minuit2;

  MnUserParameters pars;

  pars.Add("showerSize", result.fPar[0], 1.0);
  pars.SetLowerLimit("showerSize", 0);
  pars.Add("coreX", result.fPar[1], 200*meter);
  pars.Add("coreY", result.fPar[2], 200*meter);

  for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i) {
    pars.Add(GetShapeLabel(i), result.fPar[3+i], 0.1);
    if (config.fLDFShapeTreatment[i] != eFitted)
      pars.Fix(GetShapeLabel(i));
  }

  if (config.fFixCore) {
    pars.Fix("coreX");
    pars.Fix("coreY");
  }

  LDFChi2Function chi2Fcn(config, showerAxis, data);
  LDFLikelihoodFunctionTN likeFcn(config, showerAxis, data);

  MnMinimize m(fUseMaxLikelihood ?
               static_cast<const FCNBase&>(likeFcn) :
               static_cast<const FCNBase&>(chi2Fcn), pars, 1);

  FunctionMinimum fmin = m(10000);

  if (fmin.IsAboveMaxEdm() || fmin.HasReachedCallLimit() || !fmin.IsValid()) {
    OUT(eIntermediate)
      << "\n**************** Minimize FAILED ****************\n"
      << fmin
      << "**************** Minimize FAILED ****************\n\n";
    return false;
  }

  {
    ostringstream msg;
    msg << "Found minimum after " << fmin.NFcn() << " evaluations";
    INFO(msg);
  }

  MnHesse hesse(1);
  hesse(m.Fcnbase(), fmin, 10000);

  if (!fmin.IsValid() || fmin.HesseFailed()) {
    OUT(eIntermediate)
      << "\n**************** Hesse FAILED ****************\n"
      << fmin
      << "**************** Hesse FAILED ****************\n\n";
    return false;
  }

  if (fInfoLevel >= eMinuit)
    OUT(eIntermediate)
      << "\n**************** Minuit Output ****************\n"
      << fmin
      << "**************** Minuit Output ****************\n\n";

  // get fit results
  result.fPar = fmin.UserParameters().Params();

  const vector<double> modeledShape =
    config.fLDF.ShapeModel(showerAxis.GetCosTheta(gBaryCS), result.fPar[0]);
  for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i)
    if (config.fLDFShapeTreatment[i] == eModeled)
      result.fPar[3+i] = modeledShape[i];

  const size_t nVar = pars.VariableParameters();

  // get covariance
  result.fCov = 0;

  for (size_t i = 0; i < nVar; ++i)
    for (size_t j = 0; j < nVar; ++j)
      result.fCov(m.ExtOfInt(i), m.ExtOfInt(j)) = fmin.UserCovariance()(i,j);

  const double minValue = fmin.Fval();
  if (fUseMaxLikelihood) {
    result.fLogLikelihood = minValue;
    result.fChi2 = chi2Fcn(result.fPar);
  } else {
    result.fLogLikelihood = likeFcn(result.fPar);
    result.fChi2 = minValue;
  }

  result.fNdof = 0;
  for (const auto& station : data)
    if (station.fSignal > 0)
      ++result.fNdof;
  result.fNdof -= nVar;

  const double reducedChi2 = (result.fNdof > 0) ? result.fChi2/result.fNdof : 0;

  const utl::Point core = result.GetCore();

  const utl::Point coreErr(result.fCov.Std(eCoreX), result.fCov.Std(eCoreY), 0, gBaryCS);
  OUT(eIntermediate)
    << "  chi2   = " << result.fChi2 << " / " << result.fNdof << "\n"
       "  " << GetShowerSizeLabel(config) << "  = "
            << result.fPar[eShowerSize] << " +/- " << result.fCov.Std(eShowerSize) << "\n"
       "  core   = " << ToString(core, gBaryCS) << " +/-\n"
       "           " << ToString(coreErr, gBaryCS) << " (bary)";
  for (size_t i = 0, n = config.fLDF.GetNShapeParameters(); i < n; ++i)
    OUT(eIntermediate) << "\n  " << GetShapeLabel(i)
      << "  = " << result.fPar[eShapeParameters + i] << " +/- " << result.fCov.Std(eShapeParameters + i);
  OUT(eIntermediate) << endl;

  if (reducedChi2 > fMaxChi2) {
    OUT(eIntermediate)
      << "  Chi2/Ndof = " << reducedChi2 << " exceeds limit "
      << fMaxChi2 << '\n';
    return false;
  }

  const double d = (core - fBarycenter).GetMag();
  if (d > fMaxBaryToCoreDistance) {
    OUT(eIntermediate)
      << "  Distance to core " << d/km << " [km] exceeds limit.\n";
    return false;
  }

  return !HasNaN(result);
}