UnivV2Rec.cc

Go to the documentation of this file.

#include "UnivV2Rec.h"

#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <random>
#include <string>
#include <vector>

#include <boost/filesystem.hpp>
#include <boost/format.hpp>

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

#include <evt/Event.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/ShowerUnivRecData.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>
#include <fevt/FEvent.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <sevt/SEvent.h>

#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/UTCDateTime.h>

// local
#include <src/AtmosphereUtilities.h>
#include <src/GetSignal.h>
#include <src/PowerConfig.h>

// fd energy long term correction
#include <src/Long_term_correction_header.h>

using namespace utl;
using namespace std;

namespace un2 {

void SimpleReco(evt::Event &event) {

  // ------------------------------- PROLOG ---------------
  // ------------------------------------------------------
  //
  // the managment advises to listen to some relaxing music:
  // https://www.youtube.com/watch?v=cSkqscB9j4I
  //
  const CoordinateSystemPtr augerCS =
      det::Detector::GetInstance().GetReferenceCoordinateSystem();
  const auto &sDet = det::Detector::GetInstance().GetSDetector();
  const bool verbose = false;
  const Branch topB =
      fwk::CentralConfig::GetInstance()->GetTopBranch("UniversalityFitterV2");

  string hadrIntModel;
  topB.GetChild("HadronicInteractionModel").GetData(hadrIntModel);

  bool fixXmaxToFd;
  bool fixEnergyToFd;
  topB.GetChild("FixXmaxToFd").GetData(fixXmaxToFd);
  topB.GetChild("FixEnergyToFd").GetData(fixEnergyToFd);

  INFO(boost::format("Fix Xmax to FD: %b") % (fixXmaxToFd));
  INFO(boost::format("Fix energy to FD: %b") % (fixEnergyToFd));

  if (!event.GetRecShower().HasUnivRecShower())
    event.GetRecShower().MakeUnivRecShower();

  const double gcm = gram / cm2;

  if (!event.HasRecShower()) {
    ERROR("No RecShower found!");
    return;
  }

  // reconstruction algorithm of the universality-v2 framework
  // based on the model developed in this PhD thesis:
  // https://publikationen.bibliothek.kit.edu/1000144727
  //
  // PLEASE CONSIDER PULLING THE FULL REPOSITORY
  // https://gitlab.iap.kit.edu/shower-universality/universality-v2.git
  //
  // the reconstruction algorithm creates an object for each considered detector
  // in a station. the detectors are evaluated for saturation and are eventually
  // removed. for each detector a best-fit for the start time according to an
  // estimate of xmax is evaluated for the later fit, to tackle uncertainties
  // from the geometry that affect the timing. detectors with bad chi^2 are
  // removed. only detectors within 2000m in the shower plane are considered,
  // more conservative cuts are applied for time-trace fitting.
  //
  // trace fit is done on the integrated trace to achieve smoothing of the data,
  // which is one of the main differences in this approach. this allows you to
  // write down a signal ucnertainty model more easily than bin-wise for a trace
  //
  // energy is constrained to MC energy +/- 10% (15%) or taken directly either
  // from FD or SD reconstruction in case of data (instead of simulation data).
  // Rmu is evaluated using a fixed energy estimate from a first fit stage
  // and should never be fit simultaneously with the primary energy!
  // Xmax is esimated mostly from the trace information, using a constraint fit
  // around expected Xmax values based on Rmu (this can be removed, if needed.
  // removing this constraint does not change precision for protons, but for
  // iron!). lnA is estimated from results of Rmu and Xmax. data is stored in
  // ADST and python-like dictionary text file (open using the ast module).
  //
  // MC value for Rmu is calculated from weighted average of dense stations.
  // please pay attention that eHistory has to be available for the corsika
  // shower you are considering for Rmu to be calculated! If no dense stations
  // are available, Rmu is set to 666!

  // INTRODUCTION
  // ---------------------------------------------------------------------------------
  const bool realEvent = !event.HasSimShower();
  const auto &recShower = event.GetRecShower().GetSRecShower();
  double recEnergy =
      !realEvent ? event.GetSimShower().GetEnergy() : recShower.GetEnergy();
  double avgXmax = avgXmax0(log10(recEnergy), realEvent) - 45.;
  const Vector &recAxis = recShower.GetAxis();
  const Point &recCore = recShower.GetCorePosition();
  const double recTheta = recAxis.GetTheta(augerCS);
  const bool realAtmosphere = !event.HasSimShower();
  double realXmax = event.HasSimShower()
                        ? event.GetSimShower().GetGHParameters().GetXMax() / gcm
                        : avgXmax;
  const unsigned int pId =
      event.HasSimShower() ? event.GetSimShower().GetPrimaryParticle() : 666;
  const int month =
      UTCDateTime(det::Detector::GetInstance().GetTime()).GetMonth();
  const double recPhi = recAxis.GetPhi(augerCS);
  const double recAzimuth = (recPhi < 0) ? recPhi + 2 * M_PI : recPhi;

  vector<vector<double>> FdXmax = {};
  double FdMeanXmax = 0;
  double FdMeanEnergy = 0;
  double FdMeanLgEnergy = 0;
  double wEnergySum = 0;
  double wXmaxSum = 0;

  if (event.HasFEvent()) {

    int nRecEyes = 0;
    for (fevt::FEvent::ConstEyeIterator eye =
             event.GetFEvent().EyesBegin(fevt::ComponentSelector::eHasData);
         eye != event.GetFEvent().EyesEnd(fevt::ComponentSelector::eHasData);
         ++eye) {

      if (!(event.GetFEvent().HasEye(eye->GetId()) && eye->HasRecData()))
        continue;

      if (!eye->GetRecData().HasFRecShower())
        continue;

      const auto &eyeRec = eye->GetRecData();
      if (!eyeRec.GetFRecShower().HasGHParameters())
        continue;

      double eyeEnergy = eyeRec.GetFRecShower().GetTotalEnergy();
      const double eyeXmax = eyeRec.GetFRecShower().GetGHParameters().GetXMax();
      const double eyeXmaxErr =
          eyeRec.GetFRecShower().GetGHParameters().GetXMaxError();
      const double eyeEnergyErr = eyeRec.GetFRecShower().GetTotalEnergyError();

      if (true /* correct long term drift fd energy*/) {
        const double energyRatioSdFd = recEnergy / eyeEnergy;
        const double corrFactor = LongTermCorrectionPerEye(energyRatioSdFd, event.GetHeader().GetTime().GetGPSSecond() / pow(10, 6), eye->GetId() - 1);
        eyeEnergy = eyeEnergy * energyRatioSdFd/corrFactor;
        INFO(boost::format("eye (%i) correction factor %.2f") % eye->GetId() % energyRatioSdFd);
      }

      FdXmax.push_back({eyeXmax / gcm, eyeXmaxErr / gcm});
      nRecEyes += 1;
      FdMeanXmax += eyeXmax / gcm / pow(eyeXmaxErr / gcm, 2);
      FdMeanEnergy += eyeEnergy / pow(eyeEnergyErr, 2);
      wXmaxSum += 1 / pow(eyeXmaxErr / gcm, 2);
      wEnergySum += 1 / pow(eyeEnergyErr, 2);
      // cout << eyeXmax / gcm << endl;
    }
    FdMeanXmax /= wXmaxSum;
    FdMeanEnergy /= wEnergySum;
    FdMeanLgEnergy = log10(FdMeanEnergy);
    INFO(boost::format("FD mean Xmax   %.2f g/cm^2") % FdMeanXmax);
    INFO(boost::format("FD mean lgE/eV %.2f") % FdMeanLgEnergy);
  }

  if (fixXmaxToFd && FdMeanXmax > 0)
    avgXmax = FdMeanXmax;

  if (fixXmaxToFd && FdMeanXmax > 0 && realEvent)
    realXmax = FdMeanXmax;

  if (fixEnergyToFd && FdMeanEnergy > 0)
    recEnergy = FdMeanEnergy;

  INFO(boost::format("lg(E/eV) = %.2f   θ = %.2f°") % (log10(recEnergy)) %
       un2::Rad2Deg(recTheta));

  if (log10(recEnergy) < 18.5) {
    WARNING("Energy too low for universality reconstruction!");
    WARNING("Energy too low for universality reconstruction!");
    WARNING("Energy too low for universality reconstruction!");
    return; // speedup
  }

  auto &univRecShower = event.GetRecShower().GetUnivRecShower();

  univRecShower.SetCorePosition(recCore);
  univRecShower.SetAxis(recAxis);

  const auto &sEvent = event.GetSEvent();

  double RmuSum = 0;
  double wRmuSum = 0;

  ssd::ESignalVarianceModel stdSSDmodel = ssd::eGAP2019_000;
  wcd::ESignalVarianceModel stdWCDmodel = wcd::eGAP2014_035;

  // -----------------------------------------------------
  // PREAMBLE---------------------------------------------
  RmuXmaxFit powerFit;
  LDFFit LdfFit;
  LDFFit EnergyFit;
  SignalRatioFit RatioFit;

  double useCorrelationPenalty;
  bool normalizeCDF;
  topB.GetChild("UseCorrelation").GetData(useCorrelationPenalty);
  topB.GetChild("NormalizeCDF").GetData(normalizeCDF);

  powerFit.AddMetaData(avgXmax, log10(recEnergy), recTheta,
                       event.HasSimShower(), month,
                       normalizeCDF,           // normalized CDFs
                       useCorrelationPenalty); // penalty
  LdfFit.AddMetaData(avgXmax, log10(recEnergy), recTheta, event.HasSimShower(),
                     month);
  EnergyFit.AddMetaData(avgXmax, log10(recEnergy), recTheta,
                        event.HasSimShower(), month);

  RatioFit.AddMetaData(avgXmax, log10(recEnergy), recTheta,
                        event.HasSimShower(), month);

  unsigned short int numDet = 0;
  unsigned short int numEnDet = 0;
  unsigned short int numTimeDet = 0;
  unsigned short int numDenseDet = 0;
  unsigned short int numSatDet = 0;
  const int numDenseDetTot =
      sDet.GetDenseStationList().size() * 2; // AugerPrime

  const string basepath = ".";
  const string dumpfileId =
      event.GetHeader().GetId(); // boost::str(boost::format("%08X") %
                                 // (distribution(generator)));
  const string pythonDumpFilepath =
      boost::str(boost::format("%s/event_dump_universalityv2_%s.txt") %
                 basepath % dumpfileId);

  INFO(boost::format("Dumpfile will be at %s") % pythonDumpFilepath);
  // ofstream pythonDump(pythonDumpFilepath.c_str());
  ofstream pythonDump("dumpfile.txt");

  pythonDump << "{\n";
  pythonDump << boost::format("\"ID\": \"%s\",\n") % (dumpfileId);
  pythonDump << boost::format("\"GPSSecond\": %i,\n") %
                    event.GetHeader().GetTime().GetGPSSecond();
  pythonDump << boost::format("\"reclgE\": %.5f,\n") % (log10(recEnergy));
  pythonDump << boost::format("\"recTheta\": %.5f,\n") % recTheta;
  pythonDump << boost::format("\"recPhi\": %.5f,\n") % recPhi;
  pythonDump << boost::format("\"recAzimuth\": %.5f,\n") % recAzimuth;
  pythonDump << boost::format("\"realXmax\": %.5f,\n") % realXmax;
  pythonDump << boost::format("\"particleId\": %i,\n") % pId;
  pythonDump << boost::format("\"detectors\":\n[\n");

  vector<un2Detector> un2Detectors;

  // stations iteration
  for (const auto &station : boost::make_iterator_range(sEvent.StationsBegin(),
                                                        sEvent.StationsEnd())) {

    // avoid dereferencing NULL
    if (!station.HasRecData())
      continue;

    const unsigned int stationId = station.GetId();
    const bool isDense = sDet.GetStation(stationId).IsDense();
    const auto &statRecData = station.GetRecData();
    const double rsp = statRecData.GetSPDistance();

    // hard cut needed, universality parametrizations are
    // not valid beyond ~ 2000 m
    if (rsp > 2000 * utl::m)
      continue;

    const double az = statRecData.GetAzimuthShowerPlane();
    const double height =
        sDet.GetStation(stationId).GetPosition().GetZ(augerCS) + 1400.;

    const size_t sigStartSlot = statRecData.GetSignalStartSlot();
    const size_t sigStopSlot = statRecData.GetSignalEndSlot();

    const double coreTime = recShower.GetCoreTime().GetGPSNanoSecond();
    const double stationTime =
        statRecData.GetSignalStartTime().GetGPSNanoSecond();
    const double stationPFTime =
        coreTime - 1. / 0.2998 * (cos(az) * rsp * tan(recTheta));
    //    station.HasSimData()
    //        ? station.GetSimData().GetPlaneFrontTime().GetGPSNanoSecond()
    //        : coreTime - 1. / 0.2998 * (cos(az) * rsp * tan(recTheta));
    const double planeFrontTimeResidual = stationPFTime - stationTime;

    un2Detector un2WCD(station.GetId(), station, "WCD", rsp, az, height,
                       planeFrontTimeResidual, isDense);
    un2WCD.signalStartSlot = sigStartSlot;
    un2WCD.signalStopSlot = sigStopSlot;
    un2WCD.referenceDt40 =
        GetTimeQuantile(recTheta, rsp, az, height, realXmax, "WCD", "t40",
                        month, realAtmosphere)[0];
    un2WCD.referenceDt0 =
        GetTimeQuantile(recTheta, rsp, az, height, realXmax, "WCD", "t0", month,
                        realAtmosphere)[0];

    un2Detectors.push_back(un2WCD);

    // add a similar block of code for further detectors
    if (station.HasScintillator() && station.GetScintillator().HasRecData()) {
      un2Detector un2SSD(station.GetId(), station, "SSD", rsp, az, height,
                         planeFrontTimeResidual, isDense);
      un2SSD.signalStartSlot = sigStartSlot;
      un2SSD.signalStopSlot = sigStopSlot; // not part of constructor because
                                           // these might change in the future!
      un2SSD.referenceDt40 =
          GetTimeQuantile(recTheta, rsp, az, height, realXmax, "SSD", "t40",
                          month, realAtmosphere)[0];
      un2SSD.referenceDt0 =
          GetTimeQuantile(recTheta, rsp, az, height, realXmax, "SSD", "t0",
                          month, realAtmosphere)[0];

      un2Detectors.push_back(un2SSD);

    }
  }

  for (auto &un2det : un2Detectors) {

    INFO(boost::format("Detector %s %i") %un2det.type %un2det.Id);
    un2det.CalcUnit();

    if (0==un2det.signalUnit) {
      WARNING(boost::format("Could not recover signal unit! Skipping Detector!"));
      continue;
    }

    const auto &trace = un2det.GetTrace(); // not part of the detector object

    if ("WCD" == un2det.type)
      un2det.totalSignal = un2det.station.GetRecData().GetTotalSignal();

    if ("SSD" == un2det.type)
      un2det.totalSignal =
          un2det.station.GetScintillator().GetRecData().GetTotalSignal();

    un2det.stationNSTime =
        un2det.station.GetRecData().GetSignalStartTime().GetGPSNanoSecond();

    un2det.referenceSignal =
        GetSignal(hadrIntModel, log10(recEnergy), recTheta, un2det.spRadius,
                  un2det.spAzimuth, un2det.heightAboveSL, realXmax, 1.,
                  un2det.type, month, realAtmosphere, true)[0];
    un2det.uncertaintySignal =
        ("WCD" == un2det.type)
            ? utl::wcd::SignalUncertainty(stdWCDmodel, cos(recTheta),
                                          un2det.totalSignal)
            : utl::ssd::SignalUncertainty(stdSSDmodel, cos(recTheta),
                                          un2det.totalSignal);

    // small signals have close to zero uncertainty in SSDs, which
    // results in crazy outliers in dense ring Rmu calculation using the
    // weighted signal sum, so here I set a low threshold until SSD
    // signal uncertainty model is fixed
    if ("SSD" == un2det.type && un2det.uncertaintySignal < 1)
      un2det.uncertaintySignal = 1;

    un2det.DX =
        CalcDX(realAtmosphere, false /*useDL*/, un2det.spRadius,
               un2det.spAzimuth, un2det.heightAboveSL, realXmax /*estimator*/,
               recTheta, month, false /*silent*/);

    un2det.nsBin = trace.GetBinning(); // depends on electronics!

    if (un2det.isVirtual) {

      un2det.mcSpRadius = sDet.GetStation(un2det.Id).GetPosition().GetRho(
          event.GetSimShower().GetShowerCoordinateSystem());

      const auto &muonTrace =
          un2det.HasTrace(sevt::StationConstants::eMuon)
              ? un2det.GetTrace(sevt::StationConstants::eMuon)
              : un2det.GetTrace();
      const double muonSignal =
          accumulate(muonTrace.Begin(), muonTrace.End(), 0.);

      // SCALING PROBLEM NEEDS TO BE FIXED
      // Rmu from PURE muonic signal is not Rmu using tot signal
      const double muonPrediction =
          un2det.HasTrace(sevt::StationConstants::eMuon)
              ? GetSignal(hadrIntModel, log10(event.GetSimShower().GetEnergy()), recTheta,
                          un2det.mcSpRadius, un2det.spAzimuth,
                          un2det.heightAboveSL, realXmax, 1., un2det.type,
                          month, realAtmosphere, true)[2]
              : GetSignal(hadrIntModel, log10(event.GetSimShower().GetEnergy()), recTheta,
                          un2det.mcSpRadius, un2det.spAzimuth,
                          un2det.heightAboveSL, realXmax, 1., un2det.type,
                          month, realAtmosphere, true)[0];

      const double stationRmu = muonSignal * un2det.signalUnit / muonPrediction;

      RmuSum += stationRmu / pow(un2det.uncertaintySignal, 2);
      wRmuSum += 1 / pow(un2det.uncertaintySignal, 2);

      INFO(boost::format(
               "Station level Rmu for %i - %s: %.2f +/- %.2f at %.2f m") %
           un2det.Id % un2det.type % stationRmu %
           (un2det.uncertaintySignal / un2det.referenceSignal) %
           un2det.mcSpRadius);

      numDenseDet += 1;
    }

    if (un2det.station.IsLowGainSaturation() && !un2det.isVirtual)
      numSatDet += 1;

    vector<double> predIntegratedSignal;
    vector<double> fitIntegratedSignal;
    vector<double> stdIntegratedSignal;

    // detector selection
    if (!un2det.isVirtual && !un2det.station.IsLowGainSaturation() &&
        un2det.referenceSignal > 3 &&
        un2det.totalSignal > 1) { // this is needed because the SSD signal
      // uncertainty model is fucked up for small signals
      if (un2det.type == "SSD")
        un2det.weight = 1.;

      vector<double> timeVec = {0};

      // traces should be of the SAME LENGTH! otherwise distant stations (short
      // signals) do not get accounted correctly, use 4500ns, that is roughly
      // 540 (or 180) bins
      for (unsigned int i = un2det.signalStartSlot;
           i < (un2det.signalStartSlot + 4500 / un2det.nsBin); ++i) {

        timeVec.push_back(timeVec.back() + un2det.nsBin);

        // use all signal between start/stop slots
        if (i < trace.GetSize()-1 && i < un2det.signalStopSlot) {

          un2det.integratedSignal.push_back(un2det.integratedSignal.back() +
                                            trace[i] * un2det.signalUnit);
          un2det.integratedNormalizedSignal.push_back(
              un2det.integratedNormalizedSignal.back() +
              trace[i] * un2det.signalUnit / un2det.totalSignal);

        } else { // this ensures uniform length

          un2det.integratedSignal.push_back(un2det.integratedSignal.back());
          un2det.integratedNormalizedSignal.push_back(
              un2det.integratedNormalizedSignal.back());
        }

        if (un2det.integratedSignal.back() < 0.4 * un2det.totalSignal)
          un2det.t40RelToStart =
              timeVec.back(); // reassign t40 until final value is reached
      }

      un2det.Dt40 = un2det.t40RelToStart - un2det.pfTimeResidual;

      // FEEDING THE FIT MACHINE
      const string detT0ParName =
          "t0_" + un2det.type + "_" + to_string(un2det.Id);

      // ALL detectors here go into LDF fit, however only converging stations go
      // to power fit
      if (true) {
        INFO(boost::format("Detector %s %i added to LDFFit") %un2det.type %un2det.Id);
        LdfFit.AddDetectorData(un2det);
        numDet += 1;
      }

      if (true) {
        INFO(boost::format("Detector %s %i added to EnergyFit") %un2det.type %un2det.Id);
        EnergyFit.AddDetectorData(un2det);
        numEnDet += 1;
      }

      if (un2det.referenceSignal > 5) { // really just setting arbitrary limits
         // here because I think the SSD uncertainty model is not working!
         // really large signals needed, otherweise = crap
        INFO(boost::format("Detector %s %i added to RatioFit") %un2det.type %un2det.Id);
        RatioFit.AddDetector(un2det);
      }

      if (un2det.spRadius < 1800. * utl::m && un2det.referenceSignal > 5) {

        // single station level check
        RmuXmaxFit detFit;

        detFit.SetParameterDefs({{"Xmax", 750., 5.0, 650., 1050., false},
                                 {"Rmu", 1.10, 0.1, 0., 4., true},
                                 {"lgE", log10(recEnergy), .1, 18, 21., true}});

        detFit.AddMetaData(avgXmax, log10(recEnergy), recTheta,
                           event.HasSimShower(), month, true, 0.);
        detFit.AddDetectorData(un2det);

        const double t0start = un2det.referenceDt0[0] + un2det.pfTimeResidual;
        const double t0sig = un2det.referenceDt0[1];
        auto &parDefs = detFit.GetParameterDefs();
        parDefs.push_back({detT0ParName, t0start, 25., t0start - t0sig,
                           t0start + 0.5 * t0sig, false});

        utl::Minou::Minimizer<RmuXmaxFit> detMin(detFit);
        detMin.SetVerbose(verbose);
        detMin.Minimize(500000);

        const auto detRes = detMin.GetParametersAndErrors();

        vector<double> timeVecDet = {};
        const double t0 = detRes.at(3).first;

        for (unsigned int t = 0; t < un2det.integratedSignal.size(); ++t) {
          timeVecDet.push_back(t * un2det.nsBin - t0);
        }

        // predIntegratedSignal = GetIntegratedTotalSignalSum(
        predIntegratedSignal = GetIntegratedTotalSignalSum(
            timeVec, hadrIntModel, log10(recEnergy), recTheta, un2det.spRadius,
            un2det.spAzimuth, un2det.heightAboveSL, realXmax, 1.0, un2det.type,
            month, realAtmosphere, true, un2det.pfTimeResidual, un2det.nsBin,
            /*"cdf", 0, */ un2det.integratedSignal.back());

        fitIntegratedSignal = GetIntegratedTotalSignalSum(
            timeVecDet, hadrIntModel, log10(recEnergy), recTheta,
            un2det.spRadius, un2det.spAzimuth, un2det.heightAboveSL,
            detRes.at(0).first, detRes.at(1).first, un2det.type, month,
            realAtmosphere, true, un2det.pfTimeResidual, un2det.nsBin,
            /*"cdf", 0, */ un2det.integratedSignal.back());

        stdIntegratedSignal = GetIntegratedTotalSignal(
            timeVec, hadrIntModel, log10(recEnergy), recTheta, un2det.spRadius,
            un2det.spAzimuth, un2det.heightAboveSL, realXmax, 1.0, un2det.type,
            month, realAtmosphere, true, un2det.pfTimeResidual, un2det.nsBin,
            "unc", 0, un2det.integratedSignal.back());

        for (unsigned int i = 0; i < timeVecDet.size(); ++i) {

          un2det.detChiSq +=
              utl::Sqr((fitIntegratedSignal[i] - un2det.integratedSignal[i]) /
                       stdIntegratedSignal[i]);
        }

        un2det.detChiSq =
            (un2det.detChiSq == un2det.detChiSq) ? un2det.detChiSq : FLT_MAX;

        // some output to check that everything is within a physically
        // reasonable range of values
        ostringstream msg;
        msg << boost::format("%s %i\nr = %.1f m  ψ = %.0f°\n") % un2det.type %
                   un2det.Id % un2det.spRadius % Rad2Deg(un2det.spAzimuth);
        msg << boost::format("reconstructed Signal / pred. reference Signal: "
                             "%.2f / %.2f  ± %.2f\n") %
                   un2det.totalSignal % un2det.referenceSignal %
                   un2det.uncertaintySignal;
        msg << boost::format("reference ΔX = %.2f g/cm²\n") % un2det.DX;
        // msg << boost::format("PlaneFrontTimeResidual = %.2f ns \n") %
        //(-1. * un2det.pfTimeResidual);
        msg << boost::format("Δt0 = %.2f ns (ref.: %.2f ± %.2f ns)\n") %
                   (-1. * un2det.pfTimeResidual) % un2det.referenceDt0[0] %
                   un2det.referenceDt40[1];
        msg << boost::format("Δt40 = %.2f ns (ref.: %.2f ± %.2f ns)\n") %
                   (un2det.Dt40) % un2det.referenceDt40[0] %
                   un2det.referenceDt40[1];
        msg << boost::format("Detector level Xmax: %.2f ± %.4f g/cm²\n") %
                   detRes.at(0).first % detRes.at(0).second;
        msg << boost::format("Detector level Rmu: %.2f ± %.4f g/cm²\n") %
                   detRes.at(1).first % detRes.at(1).second;
        msg << boost::format("Shift t0 = %.2f ns\n") % detRes.at(3).first;
        msg << boost::format("CDF chi squared = %.2f ") % un2det.detChiSq;

        // final cut, remove station if the fit remained on the start value
        if (detRes.at(0).first == 780. || un2det.detChiSq > 5.e2 * un2det.referenceSignal/20.) {

          INFO("CDF fit did not converge, not adding detector to PowerFit!");

        } else {

          INFO(msg);
          INFO(boost::format("Detector %s %i added to PowerFit") %un2det.type %un2det.Id);

          // station level quantities that may not be needed in future
          // versions of this fitting algorithm
          un2det.detBestFitXmax = detRes.at(0).first;
          un2det.detBestFitRmu = detRes.at(1).first;
          /* detRes.at(2) skipped, energy is fixed */
          un2det.detBestFitT0 = detRes.at(3).first;

          powerFit.AddDetectorData(un2det);

          auto &parDefsPow = powerFit.GetParameterDefs();

          // MUST be fixed
          parDefsPow.push_back({detT0ParName, t0, 0.1, t0 - 10, t0 + 10, true});

          numTimeDet += 1;

        }
      }
    }

    // PYTHON DEBUG AND PLOTTY-PLOT OUTPUT
    pythonDump << boost::format("  {\n");
    pythonDump << boost::format("    \"type\": \"%s\", \n") % un2det.type;
    pythonDump << boost::format("    \"isVirtual\": \"%s\", \n") %
                      un2det.isVirtual;
    pythonDump << boost::format("    \"isSat\": \"%s\", \n") %
                      un2det.station.IsLowGainSaturation();
    pythonDump << boost::format("    \"id\": %i, \n") % un2det.Id;
    pythonDump << boost::format("    \"totalSignal\": %.3f, \n") %
                      un2det.totalSignal;
    pythonDump << boost::format("    \"uncertaintySignal\": %.3f, \n") %
                      un2det.uncertaintySignal;
    pythonDump << boost::format("    \"referenceSignal\": %.3f, \n") %
                      un2det.referenceSignal;
    pythonDump << boost::format("    \"referencet40\": %.3f, \n") %
                      un2det.referenceDt40[0];
    pythonDump << boost::format("    \"referencet0\": %.3f, \n") %
                      un2det.referenceDt0[0];
    pythonDump << boost::format("    \"r\": %.2f, \n") % un2det.spRadius;
    pythonDump << boost::format("    \"psi\": %.2f, \n") % un2det.spAzimuth;
    pythonDump << boost::format("    \"h\": %.2f, \n") % un2det.heightAboveSL;
    pythonDump << boost::format("    \"nsBin\": %.5f, \n") % un2det.nsBin;
    pythonDump << boost::format("    \"sigStartSlot\": %.5f, \n") %
                      un2det.signalStartSlot;
    pythonDump << boost::format("    \"sigStopSlot\": %.5f, \n") %
                      un2det.signalStopSlot;
    pythonDump << boost::format("    \"stationNSTime\": %.5f, \n") %
                      un2det.stationNSTime;
    pythonDump << boost::format("    \"stationPFTimeResidual\": %.5f, \n") %
                      un2det.pfTimeResidual;
    pythonDump << boost::format("    \"AoP\": %.5f, \n") % un2det.avgAoP;
    pythonDump << boost::format("    \"nPMTs\": %.5f, \n") %
                      un2det.station.GetNPMTs();
    pythonDump
        << boost::format("    \"commissionTimeGPSSecond\": %.2f, \n") %
               sDet.GetStation(un2det.Id).GetCommissionTime().GetGPSSecond();


    if (!un2det.HasTrace(sevt::StationConstants::eTotal)) {
        pythonDump << boost::format("  }\n,");
        break;
    }

    // this part is buggy, i don't know why the trace throws an
    // exception here from time to time (this was not always the case)
    const utl::TraceD& thisTrace = un2det.GetTrace();
    pythonDump << boost::format("    \"trace\": [");
    for (unsigned int i = 0; i < thisTrace.GetSize()-1; ++i)
      pythonDump << boost::format("%.4f, ") % thisTrace[i];
    pythonDump << boost::format(" ] ,\n");

    pythonDump << boost::format("    \"integratedSignal\": [");
    for (unsigned int i = 0; i < un2det.integratedSignal.size(); ++i)
      pythonDump << boost::format("%.4f, ") % un2det.integratedSignal[i];
    pythonDump << boost::format(" ] ,\n");

    pythonDump << boost::format("    \"bestFitIntegratedSignal\": [");
    for (unsigned int i = 0; i < fitIntegratedSignal.size(); ++i)
      pythonDump << boost::format("%.4f, ") % fitIntegratedSignal[i];
    pythonDump << boost::format(" ] ,\n");

    pythonDump << boost::format("    \"predIntegratedSignal\": [");
    for (unsigned int i = 0; i < predIntegratedSignal.size(); ++i)
      pythonDump << boost::format("%.4f, ") % predIntegratedSignal[i];
    pythonDump << boost::format(" ] ,\n");

    pythonDump << boost::format("    \"stdIntegratedSignal\": [");
    for (unsigned int i = 0; i < stdIntegratedSignal.size(); ++i)
      pythonDump << boost::format("%.4f, ") % stdIntegratedSignal[i];
    pythonDump << boost::format(" ] ,\n");

    pythonDump << boost::format("    \"t0\": %.5f, \n") % un2det.detBestFitT0;
    pythonDump << boost::format("    \"n\": %.5f, \n") % un2det.detBestFitN;
    pythonDump << boost::format("    \"signalUnit\": %.5f, \n") %
                      un2det.signalUnit;
    pythonDump << boost::format("    \"bestFitXmax\": %.5f, \n") %
                      un2det.detBestFitXmax;
    pythonDump << boost::format("    \"bestFitRmu\": %.5f \n") %
                      un2det.detBestFitRmu;
    pythonDump << boost::format("  }\n,");

  } // DETECTOR ITERATOR ENDS HERE

  pythonDump << boost::format("], \n");

  const double realRmu = isnan(RmuSum / wRmuSum) ? 666 : (RmuSum / wRmuSum);

  INFO(boost::format("dense ring Rmu: %.2f") % realRmu);
  INFO(boost::format("Found %i regular and %i/%i virtual dense detectors.") %
       numDet % numDenseDet % numDenseDetTot);
  INFO(boost::format("Found %i detectors for time fit.") % numTimeDet);

  double XmaxResult = 0.;
  double XmaxResultErr = 10000.;
  double RmuResult = 0.;
  double RmuResultErr = 10000.;
  double lgEResult = 0;
  double lgEResultErr = 10000.;

  // check whether to apply bias correctino to the later fit results
  bool biasCorrection;
  topB.GetChild("ApplyBiasCorrection").GetData(biasCorrection);

  // ENERGY FIT
  // this fit is only performed properly for simulations, for data the recEnergy
  // is used
  auto &newEnergyDefs = EnergyFit.GetParameterDefs();
  newEnergyDefs[0] = {"Xmax", avgXmax, .1, avgXmax - 50, avgXmax + 50, true};
  newEnergyDefs[1] = {"Rmu", 1.1, .01, .9, 1.3, true};


  double deltaLgEnergy;
  topB.GetChild("EnergyConstraint").GetData(deltaLgEnergy);

  if (0.05 > deltaLgEnergy && !(fixEnergyToFd)) {
      INFO("You are running in BENCHMARK MODE! Energy is tightly constrained to MC or S38 value!");
  }

  const double startLgRecEnergy = (fixEnergyToFd && FdMeanLgEnergy > 0) ? FdMeanLgEnergy : log10(recShower.GetEnergy());
  newEnergyDefs[2] = {"lgE", startLgRecEnergy, .015,
                      log10(recEnergy) - deltaLgEnergy, log10(recEnergy) + deltaLgEnergy, fixEnergyToFd};

  EnergyFit.SetParameterDefs(newEnergyDefs);

  if (0 == numEnDet) {
    WARNING("ENERGY CANNOT BE ESTIMATED!");
  }

  {

    utl::Minou::Minimizer<LDFFit> EnergyMin(EnergyFit);
    EnergyMin.SetVerbose(verbose);
    EnergyMin.Minimize(50000000);

    const auto energyResult = EnergyMin.GetParametersAndErrors();
    lgEResult = energyResult.at(2).first;
    lgEResultErr = energyResult.at(2).second;
    INFO(boost::format("LDF lgE  fit result: %.2f ± %.2f") % lgEResult %
         lgEResultErr);

  }

  // LDF FIT
  auto &newLDFDefs = LdfFit.GetParameterDefs();
  // constrained, let fluctuate downwards for iron
  // avgXmax is set to FD value if the respective flag is set
  newLDFDefs[0] = {"Xmax",       avgXmax,      3.,
                   avgXmax - 70, avgXmax + 40, true};
  newLDFDefs[1] = {"Rmu", 1.15, .05, 0., 3., false};
  newLDFDefs[2] = {"lgE", lgEResult, .005, lgEResult - 0.041, lgEResult + 0.041,
                   true};

  LdfFit.SetParameterDefs(newLDFDefs);
  double XmaxResultLDFOnly, RmuResultLDFOnly;

  {

    utl::Minou::Minimizer<LDFFit> LdfMin(LdfFit);
    LdfMin.SetVerbose(verbose);
    LdfMin.Minimize(50000000);

    const auto LdfResult = LdfMin.GetParametersAndErrors();

    INFO(boost::format("LDF Xmax fit result: %.2f g/cm² ± %.2f g/cm²") %
         LdfResult.at(0).first % LdfResult.at(0).second);
    INFO(boost::format("LDF Rmu  fit result: %.2f ± %.2f") %
         LdfResult.at(1).first % LdfResult.at(1).second);

    XmaxResult = LdfResult.at(0).first;

    RmuResult = LdfResult.at(1).first;
    RmuResultErr = LdfResult.at(1).second;

    lgEResult = LdfResult.at(2).first;
    // lgEResultErr = LdfResult.at(2).second; // nonsense, since it should be tightly
    // constrained

    XmaxResultLDFOnly = LdfResult.at(0).first;
    RmuResultLDFOnly = LdfResult.at(1).first;

  }


  // AugerPrime ratio fit
  RatioFit.GeneratePairs();
  auto &newRatioDefs = RatioFit.GetParameterDefs();
  newRatioDefs[0] = {"Xmax", XmaxResult, 2., 600, 1200, true};
  newRatioDefs[1] = {"Rmu", RmuResult, .05, 0, 3, false};
  // fixed
  newRatioDefs[2] = {"lgE", lgEResult, .001, lgEResult - 0.02, lgEResult + 0.02,
                     true};
  powerFit.SetParameterDefs(newRatioDefs);

  if (2 < RatioFit.NumPairs()) {
    utl::Minou::Minimizer<SignalRatioFit> RatioMin(RatioFit);
    RatioMin.SetVerbose(verbose);
    RatioMin.Minimize(50000000);

    const auto RatioResult = RatioMin.GetParametersAndErrors();

    INFO(boost::format("Ratio Rmu  fit result: %.2f ± %.2f") %
         RatioResult.at(1).first % RatioResult.at(1).second);

    RmuResult = RatioResult.at(1).first;
    RmuResultErr = RatioResult.at(1).second;

  }

  if (biasCorrection) {
    RmuResult +=
        BiasCorrectionRmu(utl::Sqr(sin(recTheta)), log10(recEnergy), realEvent);
  }


  // POWER FIT
  auto &newPowerDefs = powerFit.GetParameterDefs();
  newPowerDefs[0] = {"Xmax", avgXmax, 5., 600, 1200, fixXmaxToFd};
  // MUST be fixed in this stage if correlation is used
  newPowerDefs[1] = {"Rmu", RmuResult, .01, RmuResult - 0.05, RmuResult + 0.05,
                     true};
  // fixed
  newPowerDefs[2] = {"lgE", lgEResult, .001, lgEResult - 0.02, lgEResult + 0.02,
                     true};
  powerFit.SetParameterDefs(newPowerDefs);

  {

    utl::Minou::Minimizer<RmuXmaxFit> powerMin(powerFit);
    powerMin.SetVerbose(verbose);
    powerMin.Minimize(50000000);

    const auto powerResult = powerMin.GetParametersAndErrors();

    XmaxResult = powerResult.at(0).first;
    XmaxResultErr = powerResult.at(0).second;

    // RmuResult = powerResult.at(1).first;
    // RmuResultErr = powerResult.at(1).second;

    lgEResult = powerResult.at(2).first;
    // lgEResultErr = powerResult.at(2).second;

  }

  if (biasCorrection) {
    XmaxResult += BiasCorrectionXmax(utl::Sqr(sin(recTheta)), log10(recEnergy),
                                     realEvent);
  }

  const double recLnA = lnA(RmuResult, XmaxResult, log10(recEnergy), realEvent);

  pythonDump << boost::format("\"numDet\": %i,\n") % numDet;
  pythonDump << boost::format("\"numEnDet\": %i,\n") % numEnDet;
  pythonDump << boost::format("\"numTimeDet\": %i,\n") % numTimeDet;
  pythonDump << boost::format("\"numSatDet\": %i,\n") % numSatDet;
  pythonDump << boost::format("\"T5Status\": %i,\n") %
                    event.GetRecShower().GetSRecShower().GetT5Trigger();
  pythonDump << boost::format("\"realRmu\": %.5f,\n") % realRmu;
  pythonDump << boost::format("\"resultRmu\": %.5f,\n") % RmuResult;
  pythonDump << boost::format("\"resultXmax\": %.5f,\n") % XmaxResult;
  pythonDump << boost::format("\"resultlgE\": %.5f,\n") % lgEResult;
  pythonDump << boost::format("\"uncRmu\": %.5f,\n") % RmuResultErr;
  pythonDump << boost::format("\"uncXmax\": %.5f,\n") % XmaxResultErr;
  pythonDump << boost::format("\"LdfRmu\": %.5f,\n") % RmuResultLDFOnly;
  pythonDump << boost::format("\"LdfXmax\": %.5f,\n") % XmaxResultLDFOnly;
  pythonDump << boost::format("\"recLnA\": %.5f,\n") % recLnA;
  pythonDump << boost::format("\"FdMeanXmax\": %.5f,\n") % FdMeanXmax;
  pythonDump << boost::format("\"FdXmax\": [");
  for (unsigned int i = 0; i < FdXmax.size(); ++i)
    pythonDump << boost::format("(%.4f, %.4f), ") % FdXmax[i].at(0) %
                      FdXmax[i].at(1);
  pythonDump << boost::format("] ,\n");
  pythonDump << boost::format("}\n");
  pythonDump.close();

  INFO(boost::format("MC/average Xmax = %.2f g/cm²") % (realXmax));
  INFO(boost::format("Combined Xmax fit result: %.2f g/cm² ± %.2f g/cm²") %
       XmaxResult % XmaxResultErr);
  INFO(boost::format("Combined Rmu  fit result: %.2f ± %.2f") % RmuResult %
       RmuResultErr);
  INFO(boost::format("Xmax/Rmu correlation lnA: %.2f ") % recLnA);

  // EPILOG
  // -----------------------------------------------------------------------------------

  if (event.HasSimShower())
    event.GetSimShower().SetNmu(realRmu);

  event.GetRecShower().GetUnivRecShower().SetGoodRec(XmaxResultErr < 300. && RmuResultErr < 1.);
  event.GetRecShower().GetUnivRecShower().SetNmu(RmuResult, RmuResultErr);
  event.GetRecShower().GetUnivRecShower().SetXmax(XmaxResult * gcm,
                                                  XmaxResultErr * gcm);
  event.GetRecShower().GetUnivRecShower().SetEnergy(pow(10., lgEResult),
                                                    pow(10, lgEResultErr));
  event.GetRecShower().GetUnivRecShower().SetLnA(
      recLnA, 0. /*error propagation needed*/);

  INFO(boost::format("Ending un2 reconstruction."));
}

} // namespace un2