UniversalityFitter.cc

Go to the documentation of this file.

#include "UniversalityFitter.h"

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

#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>

#include <evt/ShowerMRecData.h>
#include <mdet/MDetector.h>
#include <mdet/Counter.h>
#include <mevt/MEvent.h>
#include <mevt/Counter.h>
#include <mevt/Module.h>
#include <mevt/ModuleRecData.h>

#include <atm/ProfileResult.h>

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

#include <sevt/SEvent.h>
#include <sevt/Header.h>
#include <sevt/SortCriteria.h>

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

#include <utl/UTCDateTime.h>
#include <utl/TimeInterval.h>
#include <utl/UTMPoint.h>
#include <utl/GeometryUtilities.h>
#include <utl/TabulatedFunction.h>

#include <tls/Atmosphere.h>

#include "FitterUtil.h"

#include <utl/RandomEngine.h>
#include <fwk/RandomEngineRegistry.h>
#include <CLHEP/Random/Randomize.h>
#include <TMinuit.h>
#include <TVector3.h>

#include <boost/format.hpp>

#include <iostream>
#include <algorithm>

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

using namespace std;

using boost::format;

using CLHEP::RandGauss;
using CLHEP::RandFlat;


namespace UniversalityFitter {

  const double gpercm2 = utl::gram / utl::cm2;
  const double averageAugerDensity = 1.06e-3 * utl::gram / utl::cm3;


  inline
  float
  InterpolateCDF(const float vem1, const float vem2, const float t1, const float t2, const float vemQ)
  {
    return (t2 - t1) / (vem2 - vem1) * (vemQ - vem1) + t1;
  }


  void
  GetTimeQuantile(const float* const trace, const float sum, const int nT, const float timeUnit, const float f, float& tQ, float& tQ_err)
  {
    const float sumQ = f * sum;
    float timeQuantile[] = { 0, 0 }; // First/Last bin crossing the signal fraction f

    float accsum = 0;
    for (int it = 0; it <= nT; ++it) {
      const float tmax = timeUnit * float(it);
      const float binVal = (it == 0 ? 0 : trace[it - 1]);
      accsum += binVal;

      if (accsum >= sumQ  && (accsum - binVal) <= sumQ) {
        if (!timeQuantile[0])
          timeQuantile[0] = InterpolateCDF(accsum - binVal, accsum, tmax - timeUnit, tmax, sumQ);
        timeQuantile[1] = InterpolateCDF(accsum - binVal, accsum, tmax - timeUnit, tmax, sumQ);
      }
    }

    tQ = 0.5*(timeQuantile[0] + timeQuantile[1]);
    tQ_err = 0.5*(timeQuantile[1] - timeQuantile[0]);
  }


  inline
  bool
  IsCloseRel(const double a, const double b)
  {
    return fabs(a / b - 1) < 1e-3;
  }


  // 1500 m SD energy resolution model
  // Parametrization of shower-to-shower fluctuations + sampling fluctuations
  double
  SdEnergyResolution(const double energy)
  {
    const double pars[2] = { 1.090029853988826319e-1, 4.353909248484404970e-1 };
    return pars[0] + pars[1] / sqrt(energy / 1e17);
  }


  VModule::ResultFlag
  UniversalityFitter::Init()
  {
    fFitterKG.reset(new UnivFitterKG::Fitter());
    gUnivRecNS.reset(new UnivRecNS::UnivRec());

    CentralConfig* const cc = CentralConfig::GetInstance();
    const Branch topB = cc->GetTopBranch("UniversalityFitter");

    fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

    //
    // Read global configuration
    //
    topB.GetChild("minTheta").GetData(fMinTheta);
    topB.GetChild("maxTheta").GetData(fMaxTheta);

    topB.GetChild("minEnergy").GetData(fMinEnergy);
    topB.GetChild("skipNonT5").GetData(fSkipNonT5);
    topB.GetChild("selectedEye").GetData(fSelectedEye);
    topB.GetChild("skipXmaxBelowGround").GetData(fSkipXmaxBelowGround);
    topB.GetChild("minDXground").GetData(fMinDXground);
    topB.GetChild("skipXmaxOutsideProfile").GetData(fSkipXmaxOutsideProfile);
    topB.GetChild("minBinsAroundMaximum").GetData(fMinBinsAroundMaximum);
    topB.GetChild("SilentMaxRadius").GetData(fSilentMaxRadius);

    topB.GetChild("activeMethod").GetData(fActiveMethod);
    topB.GetChild("writeToTextFile").GetData(fWriteToTextFile);

    topB.GetChild("useUnsaturatedTrace").GetData(fUseUnsaturatedTrace);

    {
      ostringstream msg;
      msg << "\n"
             "======= Global configuration  =======\n"
          << format("theta cut %.1f ... %.1f deg") % (fMinTheta / degree) % (fMaxTheta / degree) << '\n'
          << format("log E cut %.1f") % (log10(fMinEnergy)) << "\n"
             "skipping non T5 events " << (fSkipNonT5 ? "yes" : "no") << "\n"
             "selected FD eye " << fSelectedEye << "\n"
             "SilentMaxRadius " << fSilentMaxRadius << "\n"
             "use unsaturated trace (only in case of simulations) " << (fUseUnsaturatedTrace ? "yes" : "no") << "\n"
             "active method " << fActiveMethod;
      INFO(msg);
    }

    const Branch& calibBranch = topB.GetChild("Calibration");
    calibBranch.GetChild("CalibOpt").GetData(fCalibOpt);
    calibBranch.GetChild("RecMixture").GetData(fRecMixture);

    fFitterKG->fCalibOpt = fCalibOpt;
    fFitterKG->fRecMixture = fRecMixture;

    {
      ostringstream msg;
      msg << "\n"
             "======= Calibration settings  =======\n"
             "CalibOpt " << fCalibOpt << "\n"
             "RecMixture " << fRecMixture;
      INFO(msg);
    }

    //
    // Read configuration for Karlsruhe reconstruction
    //
    const Branch& karlsruheBranch = topB.GetChild("KarlsruheReconstruction");
    karlsruheBranch.GetChild("outputFileName").GetData(fKarlsruheConfig.outputFileName);
    karlsruheBranch.GetChild("verbosityLevel").GetData(fFitterKG->fVerbosityLevel);
    karlsruheBranch.GetChild("useSmearedMC").GetData(fKarlsruheConfig.useSmearedMC);
    karlsruheBranch.GetChild("useRealisticEnergySmearing").GetData(fKarlsruheConfig.useRealisticEnergySmearing);
    karlsruheBranch.GetChild("RecType").GetData(fKarlsruheConfig.RecType);
    karlsruheBranch.GetChild("TimeModelVersion").GetData(fKarlsruheConfig.TimeModelVersion);
    karlsruheBranch.GetChild("StartTimeBias").GetData(fFitterKG->fStartTimeBias);
    karlsruheBranch.GetChild("applyXmaxBiasCorrection").GetData(fFitterKG->fApplyXmaxBiasCorrection);
    karlsruheBranch.GetChild("useMCEnergyScale").GetData(fFitterKG->fUseMCEnergyScale);
    karlsruheBranch.GetChild("EarlyTraceFit").GetData(fFitterKG->fEarlyTraceFit);
    karlsruheBranch.GetChild("QuantileShapeFit").GetData(fFitterKG->fQuantileShapeFit);
    karlsruheBranch.GetChild("IterativeFit").GetData(fFitterKG->fIterativeFit);

    // Possible overrides of standard settings
    const Branch& KGoverrides = karlsruheBranch.GetChild("Overrides");

    KGoverrides.GetChild("doLDFFit").GetData(fFitterKG->doLDFFit);
    KGoverrides.GetChild("doShapeFit").GetData(fFitterKG->doShapeFit);
    KGoverrides.GetChild("doStartTimeFit").GetData(fFitterKG->doStartTimeFit);
    KGoverrides.GetChild("doSaturatedStartTimeFit").GetData(fFitterKG->doSaturatedStartTimeFit);

    KGoverrides.GetChild("energyFixed").GetData(fFitterKG->energyFixed);
    KGoverrides.GetChild("coreFixed").GetData(fFitterKG->coreFixed);
    KGoverrides.GetChild("axisFixed").GetData(fFitterKG->axisFixed);
    KGoverrides.GetChild("xmaxFixed").GetData(fFitterKG->xmaxFixed);
    KGoverrides.GetChild("NmuFixed").GetData(fFitterKG->NmuFixed);
    KGoverrides.GetChild("TimeModelOffsetFixed").GetData(fFitterKG->TimeModelOffsetFixed);

    string x0Mode;
    KGoverrides.GetChild("X0Mode").GetData(x0Mode);
    if (x0Mode == "Fixed")
      fFitterKG->x0Mode = UnivFitterKG::eFixed;
    else if (x0Mode == "Coupled")
      fFitterKG->x0Mode = UnivFitterKG::eCoupled;
    else
      fFitterKG->x0Mode = UnivFitterKG::eFree;

    KGoverrides.GetChild("X0XmaxModelVersion").GetData(fFitterKG->fX0XmaxModelVersion);

    // Setting recType
    fFitterKG->setRecType(fKarlsruheConfig.RecType);

    {
      ostringstream msg;
      msg << "\n"
             "======= Configuration of Karlsruhe reconstruction =======\n";
      if (fWriteToTextFile)
        msg << "writing reconstruction summary to " << fKarlsruheConfig.outputFileName << '\n';
      msg << "\n"
             "Xmax bias correction " << (fFitterKG->fApplyXmaxBiasCorrection ? "on" : "off") << "\n"
             "rec type " << fKarlsruheConfig.RecType << "\n"
             "time model version " << fKarlsruheConfig.TimeModelVersion << "\n"
             "fit time model offset " << (!fFitterKG->TimeModelOffsetFixed ? "yes" : "no") << "\n"
             "useSmearedMC " << (fKarlsruheConfig.useSmearedMC ? "yes" : "no");
      INFO(msg);
    }

    if (fWriteToTextFile)
      fitInfoKG.open(fKarlsruheConfig.outputFileName.c_str());

    fFitterKG->initTimeModel(fKarlsruheConfig.TimeModelVersion);

    //
    // Read configuration for Bariloche reconstruction
    //
    const Branch& barilocheBranch = topB.GetChild("BarilocheReconstruction");

    barilocheBranch.GetChild("outputFileName").GetData(fBarilocheConfig.outputFileName);
    barilocheBranch.GetChild("RecType").GetData(fBarilocheConfig.RecType);
    barilocheBranch.GetChild("RecDetector").GetData(fBarilocheConfig.RecDetector);
    barilocheBranch.GetChild("RecSDEUpgrade").GetData(fBarilocheConfig.RecSDEUpgrade);
    barilocheBranch.GetChild("RecSys").GetData(fBarilocheConfig.RecSys);
    barilocheBranch.GetChild("verbose").GetData(fBarilocheConfig.verbose);
    barilocheBranch.GetChild("debug").GetData(fBarilocheConfig.debug);

    {
      ostringstream msg;
      msg << "\n"
             "======= Configuration of Bariloche reconstruction =======\n";
      if (fWriteToTextFile)
        msg << "writing reconstruction summary to " << fBarilocheConfig.outputFileName << '\n';
      msg << "\n"
             "RecType " << fBarilocheConfig.RecType << "\n"
             "RecDetector " << fBarilocheConfig.RecDetector << "\n"
             "RecSDEUpgrade " << fBarilocheConfig.RecSDEUpgrade << "\n"
             "RecSys " << fBarilocheConfig.RecSys << "\n"
             "CalibOpt " << fCalibOpt << "\n"
             "RecMixture " << fRecMixture << "\n"
             "verbose " << fBarilocheConfig.verbose << "\n"
             "debug " << fBarilocheConfig.debug;
      INFO(msg);
    }

    if (fWriteToTextFile)
      fitInfoBG.open(fBarilocheConfig.outputFileName.c_str());

    return eSuccess;
  }


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

    mcPars.reset();
    mcParsFluct.reset();
    fittedPars.reset();

    const bool needMC =
      fActiveMethod == "Karlsruhe" &&
      fKarlsruheConfig.RecType == 1;

    const bool needFD =
      (
       fActiveMethod == "Bariloche" &&
       (fBarilocheConfig.RecType == 0 || fBarilocheConfig.RecType == 7)
      ) ||
      (fActiveMethod == "Karlsruhe" && fKarlsruheConfig.RecType == 3);

    // nothing to do if there is no SEvent
    if (!(event.HasSEvent() && event.HasRecShower() && event.GetRecShower().HasSRecShower())) {
      INFO("event has no reconstructed shower");
      return eContinueLoop;
    }

    SEvent& sEvent = event.GetSEvent();
    sEvent.SortStations(sevt::ByDecreasingSignal());
    const ShowerSRecData& showerSRec = event.GetRecShower().GetSRecShower();

    // reconstruction requires successful LDF reconstruction
    if (showerSRec.GetLDFRecStage() < ShowerSRecData::kLDF) {
      INFO("LDF reconstruction stage too small");
      return eContinueLoop;
    }

    if (fSkipNonT5 && !showerSRec.IsT5()) {
      INFO("reconstruction requires T5");
      return eContinueLoop;
    }

    const CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

    const double theta = showerSRec.GetAxis().GetTheta(siteCS);
    if (theta < fMinTheta || theta > fMaxTheta) {
      INFO("Zenith angle above range of parameterization");
      return eContinueLoop;
    }

    {
      const double energy = event.HasSimShower() ? event.GetSimShower().GetEnergy() : showerSRec.GetEnergy();
      if (energy < fMinEnergy) {
        INFO("energy below cut");
        return eContinueLoop;
      }
    }

    const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
    const string& eventId = event.GetHeader().GetId();
    const int sdEventId = sEvent.GetHeader().GetId();
    const UTCDateTime eventTimeUTC = sEvent.GetHeader().GetTime();

    {
      ostringstream msg;
      msg << format("Event id %s, SD event id %i. Date %02i/%02i.") % eventId % sdEventId
          % eventTimeUTC.GetMonth() % eventTimeUTC.GetYear();
      INFO(msg);
    }

    fStations.clear();
    fEventHasSaturation = false;

    for (SEvent::ConstStationIterator sIt = sEvent.StationsBegin(), end = sEvent.StationsEnd();
         sIt != end; ++sIt) {

      if (sIt->HasSimData() && sIt->GetRejectionStatus() == StationConstants::eMCInnerRadiusCut) {
        ostringstream msg;
        msg << "Station " << sIt->GetId() << " is inside MC inner radius cut. Skipping event.";
        INFO(msg);
        return eContinueLoop;
      }

      if (!(sIt->IsCandidate() || sIt->IsSilent()))
        continue;

      StationFitData sdata;

      sdata.stationId = sIt->GetId();
      sdata.position = sDetector.GetStation(*sIt).GetPosition();

      // Values in sdata are left at the default for silent stations
      if (sIt->IsSilent()) {

        if ((sdata.position - showerSRec.GetCorePosition()).GetMag() > fSilentMaxRadius)
          continue;
        sdata.isSilent = true;

      } else {
        sdata.isSilent = false;

        const StationRecData& sRec = sIt->GetRecData();
        sdata.totalSignal = sRec.GetTotalSignal();

        sdata.isSaturated = sIt->IsLowGainSaturation();
        fEventHasSaturation |= sdata.isSaturated;

        sdata.startTime = sRec.GetSignalStartTime();
        sdata.startBin = sRec.GetSignalStartSlot();
        sdata.stopBin = sRec.GetSignalEndSlot();

        sdata.t50 = sRec.GetT50();
        sdata.t10 = sdata.t50 - sRec.GetRiseTime();
        sdata.t90 = sRec.GetFallTime() + sdata.t50;

        vector< vector<double> > vemTraces;
        double aopsum = 0;

        for (sevt::Station::ConstPMTIterator pmtIter = sIt->PMTsBegin();
             pmtIter != sIt->PMTsEnd(); ++pmtIter) {

          if (!(pmtIter->HasCalibData() && pmtIter->HasRecData()))
            continue;

          const sevt::PMTRecData& pmtRec = pmtIter->GetRecData();

          if (!pmtRec.HasVEMTrace())
            continue;

          const double calib = pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge();
          aopsum += 1 / calib;

          const bool useUnsat = sdata.isSaturated && fUseUnsaturatedTrace && pmtRec.HasVEMTrace(StationConstants::eTotalNoSaturation);

          if (useUnsat)
            sdata.useUnsaturatedTrace = true;

          const StationConstants::SignalComponent comp = useUnsat ? StationConstants::eTotalNoSaturation : StationConstants::eTotal;
          const TraceD& tr = pmtRec.GetVEMTrace(comp);

          vector<double> trace;

          for (TraceD::ConstIterator tIt = tr.Begin(), end = tr.End();
               tIt != end; ++tIt)
            trace.push_back((*tIt) * calib);

          vemTraces.push_back(trace);
        }

        const unsigned int npmt = vemTraces.size();
        sdata.AreaOverPeak = aopsum / npmt;

        if (npmt) {
          const unsigned int nBins = vemTraces[0].size();

          for (unsigned int curTrace = 0; curTrace < npmt; ++curTrace) {
            if (vemTraces[curTrace].size() != nBins) {
              ERROR("traces need to have exactly the same length");
              return eFailure;
            }
          }

          vector<double> averageTrace(nBins);

          for (unsigned int curBin = 0; curBin < nBins; ++curBin) {
            double binSum = 0;
            for (unsigned int curTrace = 0; curTrace < npmt; ++curTrace)
              binSum += vemTraces[curTrace][curBin];
            averageTrace[curBin] = binSum / npmt;
          }

          sdata.trace = averageTrace;
        }
      }
      fStations.push_back(sdata);
    }

    fHottestStation = &fStations[0];
    for (std::vector<StationFitData>::iterator sIt = fStations.begin(); sIt != fStations.end(); ++sIt) {
      if (sIt->totalSignal > fHottestStation->totalSignal)
        fHottestStation = &(*sIt);
    }

    fHottestStationHeight = utl::UTMPoint(fHottestStation->position, ReferenceEllipsoid::GetWGS84()).GetHeight();

    const atm::Atmosphere& theAtm = det::Detector::GetInstance().GetAtmosphere();
    const atm::ProfileResult& tempProfile = theAtm.EvaluateTemperatureVsHeight();
    const atm::ProfileResult& pressureProfile = theAtm.EvaluatePressureVsHeight();
    const double temperature = tempProfile.Y(fHottestStationHeight);
    const double pressure = pressureProfile.Y(fHottestStationHeight);

    if (!pressure || !temperature)
      fDensity = averageAugerDensity;
    else
      fDensity = (pressure * kDryAirMolarMass) / (kMolarGasConstant * temperature);

    if (event.HasSimShower()) {
      ShowerSimData& simShower = event.GetSimShower();

      mcPars.xmax = simShower.GetGHParameters().GetXMax();
      mcPars.x0 = simShower.GetXFirst();
      mcPars.xmaxMu = simShower.GetXmaxMu();
      mcPars.Nmu = 1.0;
      mcPars.energy = simShower.GetEnergy();
      mcPars.axis = -simShower.GetDirection();  // GetDirection() points downwards, GetAxis() points upwards
      mcPars.core = simShower.GetPosition();
      mcPars.coretime = simShower.GetTimeStamp();

      const CoordinateSystemPtr coreCS = simShower.GetLocalCoordinateSystem();
      const double theta = mcPars.axis.GetTheta(coreCS);
      const double phi = mcPars.axis.GetPhi(coreCS);

      if (fSkipXmaxBelowGround) {
        AtmosphereNS::Atmosphere fatm;
        fatm.SetCurrentAtmosphere(eventTimeUTC.GetMonth());
        const double DXground =
          fatm.Get_DX_DL(0*cm, 90*degree, mcPars.xmax / gpercm2,
                         theta, fHottestStationHeight / cm, false, true) * gpercm2;

        if (DXground < fMinDXground) {
          ostringstream msg;
          msg << format("Xmax too close or below ground (DX = %.1f), skipping event.") % (DXground / gpercm2);
          INFO(msg);
          return eContinueLoop;
        }
      }

      if (fSkipXmaxOutsideProfile) {
        if (!simShower.HasLongitudinalProfile())
          ERROR("Can not check if Xmax is inside the longitudinal profile.");
        const utl::TabulatedFunction& longProfile = simShower.GetdEdX();
        int nlow = 0;
        int nup = 0;
        for (unsigned int i = 0; i < longProfile.GetNPoints(); ++i) {
          if (longProfile[i].X() < mcPars.xmax)
            ++nlow;
          else
            ++nup;
        }

        if (nlow < fMinBinsAroundMaximum || nup < fMinBinsAroundMaximum) {
          INFO("Xmax is not contained in the longitudinal profile, skipping event.");
          return eContinueLoop;
        }
      }

      const CoordinateSystemPtr showerPlaneCS(
        CoordinateSystem::RotationY(
          theta,
          CoordinateSystem::RotationZ(
            phi,
            coreCS
          )
        )
      );

      // calculate MC Nmu
      // this should be generalized if only dense rings at other radii are available (for example infill)

      double muonSignalSum = 0;
      int nDenseStations = 0;

      for (SEvent::ConstStationIterator sIt = sEvent.StationsBegin(), end = sEvent.StationsEnd(); sIt != end; ++sIt)  {
        if (!sDetector.GetStation(*sIt).IsDense())
          continue;
        const double rhoSP = sDetector.GetStation(*sIt).GetPosition().GetRho(showerPlaneCS);
        if (!IsCloseRel(rhoSP, 1000*meter))
          continue;

        if (sIt->HasVEMTrace(sevt::StationConstants::eMuon)) {
          // assuming all PMTs have the same peak/charge (true in MC)
          const TraceD& trace = sIt->GetVEMTrace(sevt::StationConstants::eMuon);
          const double muonSignal = accumulate(trace.Begin(), trace.End(), 0.);
          const sevt::PMTRecData& pmtRec = sIt->GetPMT(1).GetRecData();
          muonSignalSum += muonSignal * pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge();
          ++nDenseStations;
        }
      }

      if (muonSignalSum > 0) {
        const double Nmu = 1.0; // proton QGSJet II-03
        UnivParamNS::UnivParam fsignal(eWCD);
        UnivCalibConstantsNS::UnivCalibConstants fcalib(0, 0, 0);
        const double sMuRef =
          fsignal.GetSignal(1000*meter / cm, mcPars.xmax / gpercm2, log10(mcPars.energy), theta, 90*degree,
                            fDensity / (gram / cm3), fHottestStationHeight / cm, Nmu, 0, eventTimeUTC.GetMonth());
        mcPars.Nmu = (muonSignalSum / nDenseStations) / sMuRef;
        simShower.SetNmu(mcPars.Nmu);
        ostringstream msg;
        msg << "True Nmu " << mcPars.Nmu;
        INFO(msg);
      } else {
        INFO("Could not set true Nmu.");
      }

      if (fKarlsruheConfig.useRealisticEnergySmearing) {
        const double sigmaE = SdEnergyResolution(mcPars.energy);
        mcParsFluct.energy = SmearValue(mcPars.energy, sigmaE * mcPars.energy);
      } else
        mcParsFluct.energy = SmearValue(mcPars.energy, 0.1 * mcPars.energy);

      mcParsFluct.xmax = SmearValue(mcPars.xmax, 20 * gpercm2);
      mcParsFluct.Nmu = SmearValue(mcPars.Nmu, 0.1);
      mcParsFluct.x0 = SmearValue(mcPars.x0, 20 * gpercm2);

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

        TVector3 theCR;
        theCR.SetMagThetaPhi(1., theta, phi);
        TVector3 rndD = theCR.Orthogonal();
        rndD.SetMag(fabs(RandGauss::shoot(&fRandomEngine->GetEngine(), 0, UnivRecNS::SmearingAngle * degree)));
        rndD.Rotate(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 2 * kPi), theCR);
        rndD = theCR + rndD;
        theCR = rndD.Unit();

        mcParsFluct.axis = Vector(1, theCR.Theta(), theCR.Phi(), coreCS, Vector::kSpherical);
      }

      // TODO check resolution of core position and time
      mcParsFluct.core = mcPars.core + Vector(RandGauss::shoot(&fRandomEngine->GetEngine(), 0, 50*meter),
                                              RandFlat::shoot(&fRandomEngine->GetEngine(), 0, kPi),
                                              0, coreCS, Vector::kCylindrical);

      mcParsFluct.coretime = mcPars.coretime + TimeInterval(RandGauss::shoot(&fRandomEngine->GetEngine(), 0, 100*ns));

    } else if (needMC)
      ERROR("The event has no simulated shower");

    fdParameters.clear();

    if (event.HasFEvent()) {

      const fevt::FEvent& fEvent = event.GetFEvent();
      double bestErrorEllipse = 0;

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

        if (!eyeIt->HasRecData())
          continue;

        const int eyeId = eyeIt->GetId();
        if (fSelectedEye > 0 && eyeId != fSelectedEye)
          continue;

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

        if (!eyeRec.HasFRecShower())
          continue;

        const evt::ShowerFRecData& showerFRec = eyeRec.GetFRecShower();

        if (!showerFRec.HasGHParameters())
          continue;

        const fdet::Eye& dEye = det::Detector::GetInstance().GetFDetector().GetEye(*eyeIt);

        fdParameters[eyeId].eyeName = dEye.GetName();
        fdParameters[eyeId].xmax = showerFRec.GetGHParameters().GetXMax();
        fdParameters[eyeId].xmaxUnc = showerFRec.GetXmaxError();
        fdParameters[eyeId].energy = showerFRec.GetTotalEnergy();
        fdParameters[eyeId].energyUnc = showerFRec.GetTotalEnergyError();

        const Point eyePos = dEye.GetPosition();
        const Vector vecEyeCore = eyePos - showerFRec.GetCorePosition();

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

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

        const double errorPerpSDP = distEyeCore * sin(sdpPhiErr);
        const double coreErrorEllipse = errorInSDP * errorPerpSDP * kPi;

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

        bestErrorEllipse = coreErrorEllipse;

        fdParameters[0].eyeName = dEye.GetName();
        fdParameters[0].axis = showerFRec.GetAxis();
        fdParameters[0].core = showerFRec.GetCorePosition();
        fdParameters[0].coretime = showerFRec.GetCoreTime();
      }

      if (fdParameters.empty() && needFD) {
        INFO("Event has no suitable FD data");
        return eContinueLoop;
      } else {
        ostringstream msg;
        msg << "shower geometry (axis, core position + time) taken from " << fdParameters[0].eyeName;
        INFO(msg);
      }

      double wghtAvEnergy = 0;
      double wghtAvXmax = 0;

      if (!fdParameters.empty()) {
        {
          double weightSum = 0;
          for (map<int, FDParameters>::iterator fd = fdParameters.begin(); fd != fdParameters.end(); ++fd) {
            if (!fd->second.energy)
              continue;
            const double unc = fd->second.energyUnc;
            const double weight = 1 / (unc * unc);
            weightSum += weight;
            wghtAvEnergy += fd->second.energy * weight;
          }
          wghtAvEnergy /= weightSum;
        }

        {
          double weightSum = 0;
          for (map<int, FDParameters>::iterator fd = fdParameters.begin(); fd != fdParameters.end(); ++fd) {
            if (!fd->second.xmax)
              continue;
            const double unc = fd->second.xmaxUnc;
            const double weight = 1 / (unc * unc);
            weightSum += weight;
            wghtAvXmax += fd->second.xmax * weight;
          }
          wghtAvXmax /= weightSum;
        }

        ostringstream msg;
        msg << "weighted averages\n"
               "log10 E = " << format("%.2f") % log10(wghtAvEnergy) << "  "
               "Xmax = " << format("%.2f") % (wghtAvXmax / gpercm2) << "  "
               "nEyes = " << format("%d") % (fdParameters.size() - 1); // one "virtual" eye for weighted values

        INFO(msg);

        fdParameters[0].xmax = wghtAvXmax;
        fdParameters[0].energy = wghtAvEnergy;
      }
    }

    if (fActiveMethod == "Karlsruhe")
      RunKarlsruheReconstruction(event);
    else if (fActiveMethod == "Bariloche")
      RunBarilocheReconstruction(event);

    WriteRecParameters(event);
    SetRecData(event);

    return eSuccess;
  }


  double
  UniversalityFitter::SmearValue(const double mean, const double spread)
  {
    return RandGauss::shoot(&fRandomEngine->GetEngine(), mean, spread);
  }


  void
  UniversalityFitter::SwitchMethod()
  {
    string newMethod;
    if (fActiveMethod == "Karlsruhe")
      newMethod = "Bariloche";
    else
      newMethod = "Karlsruhe";

    ostringstream msg;
    msg << "Switching to " << newMethod << " reconstruction";
    INFO(msg);
    fActiveMethod = newMethod;
  }


  void
  UniversalityFitter::RunKarlsruheReconstruction(evt::Event& event)
  {
    fFitterKG->reset();
    fFitterKG->isMCEvent = event.HasSimShower();

    if (event.HasSimShower()) {
      MCParameters& pars = fKarlsruheConfig.useSmearedMC ? mcParsFluct : mcPars;

      fFitterKG->mcPars.valid = true;
      fFitterKG->mcPars.energy = pars.energy;
      fFitterKG->mcPars.Nmu = pars.Nmu;
      fFitterKG->mcPars.xmax = pars.xmax;
      fFitterKG->mcPars.xmaxMu = pars.xmaxMu;
      fFitterKG->mcPars.x0 = pars.x0;
      fFitterKG->mcPars.axis = pars.axis;
      fFitterKG->mcPars.core = pars.core;
      fFitterKG->mcPars.coretime = pars.coretime;
    }

    const SEvent& sEvent = event.GetSEvent();
    const UTCDateTime eventTimeUTC = sEvent.GetHeader().GetTime();
    fFitterKG->atmModel = eventTimeUTC.GetMonth();
    fFitterKG->fDensity = fDensity;

    ShowerSRecData& showerSRec = event.GetRecShower().GetSRecShower();
    fFitterKG->sdPars.valid = true;
    fFitterKG->sdPars.showersize = showerSRec.GetShowerSize();
    fFitterKG->sdPars.energy = showerSRec.GetEnergy();
    fFitterKG->sdPars.axis = showerSRec.GetAxis();
    fFitterKG->sdPars.core = showerSRec.GetCorePosition();
    fFitterKG->sdPars.coretime = showerSRec.GetCoreTime();

    fFitterKG->sdPars.showersizeUnc = showerSRec.GetShowerSizeError();
    fFitterKG->sdPars.energyUnc = showerSRec.GetEnergyError();
    fFitterKG->sdPars.axisUncTheta = showerSRec.GetThetaError();
    fFitterKG->sdPars.axisUncPhi = showerSRec.GetPhiError();

    if (!fdParameters.empty()) {
      // always uses weighted average of Xmax and energy and best core / axis
      fFitterKG->fdPars.valid = true;
      fFitterKG->fdPars.energy = fdParameters[0].energy;
      fFitterKG->fdPars.xmax = fdParameters[0].xmax;
      fFitterKG->fdPars.axis = fdParameters[0].axis;
      fFitterKG->fdPars.core = fdParameters[0].core;
      fFitterKG->fdPars.coretime = fdParameters[0].coretime;
    }

    for (std::vector<StationFitData>::iterator sIt = fStations.begin(); sIt != fStations.end(); ++sIt) {

      if (sIt->isSilent || sIt->stationType != eWCD)
        continue;

      UnivFitterKG::StationFitData sdata;

      sdata.stationId = sIt->stationId;
      sdata.totalSignal = sIt->totalSignal;
      sdata.isSaturated = sIt->isSaturated;
      sdata.useUnsaturatedTrace = sIt->useUnsaturatedTrace;
      sdata.trace = sIt->trace;
      sdata.startTime = sIt->startTime;
      sdata.startBin = sIt->startBin;
      sdata.stopBin = sIt->stopBin;
      sdata.position = sIt->position;
      sdata.t10 = sIt->t10;
      sdata.t50 = sIt->t50;
      sdata.t90 = sIt->t90;

      fFitterKG->addStationFitData(sdata);
    }

    const bool fitSuccessful = fFitterKG->run();

    fittedPars.energy = fFitterKG->fittedPars.energy;
    fittedPars.energyUnc = fFitterKG->fittedPars.energyUnc;

    fittedPars.xmax = fFitterKG->fittedPars.xmax;
    fittedPars.xmaxUnc = fFitterKG->fittedPars.xmaxUnc;

    fittedPars.Nmu = fFitterKG->fittedPars.Nmu;
    fittedPars.NmuUnc = fFitterKG->fittedPars.NmuUnc;

    fittedPars.TimeModelOffset = fFitterKG->fittedPars.TimeModelOffset;
    fittedPars.TimeModelOffsetUnc = fFitterKG->fittedPars.TimeModelOffsetUnc;

    fittedPars.axis = fFitterKG->fittedPars.axis;
    fittedPars.axisUncTheta = fFitterKG->fittedPars.axisUncTheta;
    fittedPars.axisUncPhi = fFitterKG->fittedPars.axisUncPhi;

    fittedPars.core = fFitterKG->fittedPars.core;
    fittedPars.coreUncX = fFitterKG->fittedPars.coreUncX;
    fittedPars.coreUncY = fFitterKG->fittedPars.coreUncY;

    fittedPars.coretime = fFitterKG->fittedPars.coretime;

    fittedPars.fNCandShape = fFitterKG->fNStationsShape;
    fittedPars.fNCandStartTime = fFitterKG->fNStationsStartTime;
    fittedPars.fNCandLDF = fFitterKG->fNStationsLDF;

    fittedPars.isGoodRec = fitSuccessful;

    if (!fitSuccessful) {
      ostringstream msg;
      msg << "\n\n"
             "#############################\n"
             "### FIT DID NOT CONVERGE! ###\n"
             "#############################";
      INFO(msg);
    }
  }


  void
  UniversalityFitter::UnivRecNS_FCN(Int_t& /*npar*/, Double_t* /*gin*/, Double_t& f, Double_t* par, Int_t /*iflag*/)
  {
    std::vector<double> par_v;
    for (unsigned int ipar = 0; ipar < UnivRecNS::npar; ++ipar)
      par_v.push_back(par[ipar]);
    f = -gUnivRecNS->GetLikelihood(par_v, parStatus);
    if (std::isnan(f) || std::isinf(f))
      f = 100;
  }


  bool
  UniversalityFitter::Rec(UnivRecNS::UnivRec& theRec)
  {
    double XmaxLow = 0;
    double XmaxHigh = 0;
    theRec.GetXmaxRange(XmaxLow, XmaxHigh);

    //int status=0;
    std::vector<double> par;
    std::vector<double> epar;
    std::vector<double> low;
    std::vector<double> high;
    std::vector<bool> hasLimits;
    theRec.InitFCNParameters(par, epar, parStatus, hasLimits);

    // TODO switch to Minuit 2 (UnivFCN), gUnivRecNS -> local

    int ierflg = 0;
    double arglist[10];
    TMinuit* const minuit = new TMinuit(UnivRecNS::npar);
    arglist[0] = -1;
    minuit->mnexcm("SET PRINTOUT", arglist, 1, ierflg);
    minuit->mnexcm("SET NOWARNINGS", arglist, 0, ierflg);
    minuit->SetFCN(UnivRecNS_FCN);
    arglist[0] = 0.5;
    minuit->mnexcm("SET ERR", arglist, 1, ierflg);

    for (unsigned int ipar = 0; ipar < UnivRecNS::npar; ++ipar) {
      double low = 0;
      double high = 0;
      if (hasLimits[ipar]) {
        low = (ipar == 4 || ipar == 8 ? XmaxLow / UnivRecNS::unit[ipar]  :  UnivRecNS::low[ipar] / UnivRecNS::unit[ipar]);
        high = (ipar == 4 || ipar == 8 ? XmaxHigh / UnivRecNS::unit[ipar] :  UnivRecNS::high[ipar] / UnivRecNS::unit[ipar]);
      }
      minuit->mnparm(ipar, UnivRecNS::ParName[ipar].c_str(), par[ipar], epar[ipar], low, high, ierflg);
    }

    for (unsigned int ipar = 0; ipar < UnivRecNS::npar; ++ipar)
      if (parStatus[ipar] > 1)
        minuit->FixParameter(ipar);

    arglist[0] = 5000;
    arglist[1] = 0.001;

    int stdo = 0;
    if (!fBarilocheConfig.verbose) { // removing minuit output
      stdo = dup(1);
      close(1);
    }
    minuit->mnexcm("SIMPLEX", arglist , 2, ierflg);
    if (!fBarilocheConfig.verbose) { // removing minuit output
      dup2(stdo, 1);
      close(stdo);
    }

    for (unsigned int ipar = 0; ipar < UnivRecNS::npar; ++ipar)
      minuit->GetParameter(ipar, par[ipar], epar[ipar]);
    theRec.SaveFCNParameters(par, parStatus);

    delete minuit;
    return true;
  }


  bool
  UniversalityFitter::InitBarilocheReconstruction(evt::Event& event)
  {
    const int gRecSys = fBarilocheConfig.RecSys;
    const int gRecType = fBarilocheConfig.RecType;
    //const bool gRecSDE_Upgrade = fBarilocheConfig.RecSDEUpgrade;
    const bool gRecMC = !fBarilocheConfig.IsData;

    // SEvent
    const SEvent& sEvent = event.GetSEvent();
    int sdEventId = sEvent.GetHeader().GetId();
    const ShowerSRecData& showerSRecData = event.GetRecShower().GetSRecShower();
    const utl::Point& core = showerSRecData.GetCorePosition();

    // Coordinate system
    const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(core);
    const CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    const CoordinateSystemPtr hottestStationCS = LocalCoordinateSystem::Create(fHottestStation->position);

    // UnivRec
    UnivRecNS::UnivRec& theRec = *gUnivRecNS;
    UnivRecNS::MCInfo_t& MCinfo = theRec.MCinfo;
    UnivRecNS::AugerRecInfo_t& augerinfo = theRec.augerinfo;
    UnivRecNS::RecInfo_t& recinfo = theRec.recinfo;

    theRec.ClearRecEvent();

    theRec.nRecStations = 0;
    theRec.XmaxLow = UnivRecNS::low[4];
    theRec.XmaxHigh = UnivRecNS::high[4];

    //-----------------------------------------------
    //  Set Id, Check T4/6T5
    //-----------------------------------------------
    recinfo.id = sdEventId;
    const bool okT4 = showerSRecData.IsT4();
    const bool ok6T5 = showerSRecData.IsT5();

    //-----------------------------------------------
    //  Load MC info
    //-----------------------------------------------

    if (gRecMC) {
      MCinfo.theta = mcPars.axis.GetTheta(hottestStationCS);
      MCinfo.azi = mcPars.axis.GetPhi(hottestStationCS);

      MCinfo.xgcore = mcPars.core.GetX(hottestStationCS) / cm;
      MCinfo.ygcore = mcPars.core.GetY(hottestStationCS) / cm;
      MCinfo.zgcore = (mcPars.core.GetZ(hottestStationCS) + fHottestStationHeight) / cm;

      MCinfo.T0 = 0.;// mcPars.coretime;

      MCinfo.logE = log10(mcPars.energy);
      MCinfo.Xmax = mcPars.xmax / gpercm2;
    }

    //-----------------------------------------------
    //Set : Atmosphere/Ground Density
    //            Warning:  1) Ground density taken from the Offline atmosphere (not weather stations)
    //                      2) We assume that the month in the simulations was simulated with the corresponding atmospheric monthly model
    //-----------------------------------------------

    const TimeStamp& eventTime = sEvent.GetHeader().GetTime();
    const UTCDateTime eventTimeUTC = eventTime;
    theRec.GPSsec = eventTime.GetGPSSecond();

    const int imonth = eventTimeUTC.GetMonth();
    theRec.fatm->SetCurrentAtmosphere(imonth);
    recinfo.iatm = imonth;
    recinfo.AugerDensity = fDensity / (gram / cm3);

    //-----------------------------------------------
    // Set true Nmu in MC ( using the WCD response and proton QGSJEtII-03 at r=1000 [m]  and psi=90 [deg] as reference )
    //-----------------------------------------------
    if (gRecMC)
      MCinfo.Nmu = mcPars.Nmu;

    //-----------------------------------------------
    // Set OffsetM_Mu
    //-----------------------------------------------
    recinfo.OffsetM_Mu = 0;
    for (int itype = 0; itype < 3; ++itype)
      theRec.ftime[itype]->SetOffsetM_Mu(recinfo.OffsetM_Mu);

    //-----------------------------------------------
    //Set Auger Info
    //-----------------------------------------------

    // SD (Observer reconstruction)    core position is not saved (Do it in Offline!)
    augerinfo.logE_SD = log10(showerSRecData.GetEnergy());
    augerinfo.logE_SD_err = log10(1 + showerSRecData.GetEnergyError() / showerSRecData.GetEnergy());
    augerinfo.theta_SD = showerSRecData.GetAxis().GetTheta(hottestStationCS);
    augerinfo.azi_SD = showerSRecData.GetAxis().GetPhi(hottestStationCS);
    augerinfo.xgcore_SD = showerSRecData.GetCorePosition().GetX(hottestStationCS) / cm;
    augerinfo.ygcore_SD = showerSRecData.GetCorePosition().GetY(hottestStationCS) / cm;
    augerinfo.zgcore_SD = (showerSRecData.GetCorePosition().GetZ(hottestStationCS) + fHottestStationHeight) / cm;
    augerinfo.logE_SD += log10(UnivCalibConstantsNS::fEnergySys_Auger[gRecSys]);
    augerinfo.isGoodSD = (okT4 && ok6T5);

    recinfo.RA = 0;
    recinfo.Dec = 0;

    augerinfo.Xmax_FD_err = 20;
    augerinfo.logE_FD_err = log10(1.1);

    if (gRecMC) {
      //--- Fluctuated values of Energy and Xmax
      augerinfo.Xmax_FD = mcParsFluct.xmax / gpercm2;
      augerinfo.logE_FD = log10(mcParsFluct.energy);
      augerinfo.theta_FD = mcParsFluct.axis.GetTheta(hottestStationCS);
      augerinfo.azi_FD = mcParsFluct.axis.GetPhi(hottestStationCS);

      augerinfo.isGoodFD = true;

    } else {
      augerinfo.isGoodFD = false;
      if (!fdParameters.empty()) {
        const int EyeSel = (gRecSys == 6 || gRecSys == 7 || gRecSys == 8 || gRecSys == 9 ? gRecSys - 5 : 0);
        augerinfo.isGoodFD = true;

        for (int i = 0; i < 5; ++i)
          cout << "  Eyes " << i << "  " << fdParameters[i].energy << "  " << fdParameters[EyeSel].axis.GetTheta(coreCS) << endl;
        augerinfo.logE_FD = log10(fdParameters[EyeSel].energy); // If Offline is using M. Tueros et al
        augerinfo.logE_FD += log10(UnivCalibConstantsNS::fEnergySys_Auger[gRecSys]);

        augerinfo.Xmax_FD = fdParameters[EyeSel].xmax / gpercm2;
        augerinfo.Xmax_FD += UnivCalibConstantsNS::fXmaxSys_Auger[gRecSys];

        augerinfo.theta_FD = fdParameters[EyeSel].axis.GetTheta(coreCS);
        augerinfo.azi_FD = fdParameters[EyeSel].axis.GetPhi(coreCS);
      }
    }

    //-----------------------------------------------
    // Selection
    //-----------------------------------------------
    if (!augerinfo.isGoodSD) {
      printf(" Event is not T4/6T5  %d \n", recinfo.id);
      return false;
    }

    if (!augerinfo.isGoodFD && gRecType == 0) {
      printf(" Nmu_FD reconstruction but no FD (Xmax,Energy) %d \n", recinfo.id);
      return false;
    }

    //int nSignalStations = 0;

    //-----------------------------------------------
    // Set station vector  ( WCD )
    //-----------------------------------------------
    for (std::vector<StationFitData>::iterator sIt = fStations.begin(); sIt != fStations.end(); ++sIt) {
      const int itype = eWCD;
      if (theRec.gDetectorTypeMask[itype])
        continue;

      //----Stations rejected
      if (sIt->stationId == 94 && recinfo.id == 21354417)
        continue;

      //----Station/Event trigger and selection given by WCD
      UnivRecNS::RecStation_t station;
      theRec.InitRecStation(&station);

      station.type = itype;
      station.id = sIt->stationId;

      station.xg = sIt->position.GetX(hottestStationCS) / cm;
      station.yg = sIt->position.GetY(hottestStationCS) / cm;
      station.zg = (sIt->position.GetZ(hottestStationCS) + fHottestStationHeight) / cm;

      station.GPSnano = (sIt->startTime - fHottestStation->startTime).GetInterval();
      station.GPSnano -= sIt->startBin * 25.;

      if (sIt->isSilent) {
        theRec.fstation.push_back(station);
        continue;
      }

      if (sIt->trace.empty())
        continue;

      //----------------
      // T1/T2 stations
      //----------------

      //---- Area over peak
      station.AreaOverPeak = sIt->AreaOverPeak;

      //---- Traces
      const float timeUnit = 25.;
      const vector<float> vemtrace(sIt->trace.begin(), sIt->trace.end());

      //---- Saturation/Signal
      station.isSaturated = (sIt->isSaturated && !sIt->useUnsaturatedTrace);
      const int start = sIt->startBin;
      const int end = sIt->stopBin;
//#warning by offline trace convention here should be stop+1
      station.signalsize = accumulate(&vemtrace.front() + sIt->startBin, &vemtrace.front() + sIt->stopBin, 0.f);

      if (station.isSaturated)
        recinfo.IsSaturated[itype] = true;
      if (station.signalsize > UnivRecNS::vemSignalCut)
        station.mask_signal = false;

      //----Start time
      station.starttime = (start + 0.5) * timeUnit;

      //----Quantiles
      station.quantile[0] = 0.01 * int(100. / station.signalsize);
      if (station.quantile[0] < 0.01)
        station.quantile[0] = 0.01;
      station.quantile[1] = 0.100;
      station.quantile[2] = 0.500;
      station.quantile[3] = 0.900;

      for (int iq = 0; iq < UnivRecNS::nQ; ++iq) {
        float tq = 0;
        float tq_err = 0;
        const float signal = station.signalsize;
        GetTimeQuantile(&vemtrace[start], signal, end - start, timeUnit,  station.quantile[iq], tq, tq_err);
        station.tquantile[0][iq] = tq + start * timeUnit;
      }
      for (int iq = 0; iq < UnivRecNS::nQ; ++iq)
        if (station.signalsize < UnivRecNS::vemTimeCut[iq] || station.isSaturated) {
          station.quantile[iq] = UNDEF;
          station.tquantile[0][iq] = UNDEF;
          station.tquantile[1][iq] = UNDEF;
        }

      station.quantile[3] = UNDEF;
      station.tquantile[0][3] = UNDEF; // Disable t_90 until new Time model is produced

      if (station.isSaturated)
        station.mask_starttime = false;
      if (station.quantile[0] != UNDEF || station.quantile[1] != UNDEF ||
          station.quantile[2] != UNDEF || station.quantile[3] != UNDEF)
        station.mask_quantiles = false;

      const double vemCut = !theRec.IsTimeRec ? UnivRecNS::vemSignalCut : UnivRecNS::vemTimeCut[2];
      if (station.signalsize > vemCut)
        ++recinfo.nStationsWCD;

      theRec.fstation.push_back(station);
    } // Stations


    //-----------------------------------------------
    // Set station vector  ( MD )
    //-----------------------------------------------
    if (event.HasMEvent()) {
      // Select counters with data
      mevt::MEvent& me = event.GetMEvent();
                  const mdet::MDetector& mDetector = det::Detector::GetInstance().GetMDetector();

      for (mevt::MEvent::CounterIterator ic = me.CountersBegin(), ec = me.CountersEnd(); ic != ec; ++ic) {

        if (!ic->HasRecData() || !sEvent.HasStation(ic->GetId()) || ic->IsRejected() || !ic->IsCandidate())
          continue;
        //Reject dense array stations
       /* if (ic->GetId() > 1080000)
            continue;*/
        const int itype = eAMIGA;
        if (theRec.gDetectorTypeMask[itype])
          continue;

        //----Station info
        UnivRecNS::RecStation_t station;
        theRec.InitRecStation(&station);

        station.type = itype;
        station.id = ic->GetId();

                                utl::Point position = mDetector.GetCounter(station.id).GetPosition();
        station.xg = position.GetX(hottestStationCS) / cm;
        station.yg = position.GetY(hottestStationCS) / cm;
        station.zg = (position.GetZ(hottestStationCS) + fHottestStationHeight) / cm;

        station.GPSnano = 0;
        //Detector Area should be in cm^2
        station.DetectorArea = ic->GetActiveArea()*10000;

        //Saturation and signal
        if (ic->IsSaturated())
          station.isSaturated = true;
        else
          station.isSaturated = false;
        if (station.isSaturated)
          recinfo.IsSaturated[itype] = true;

        station.signalsize = ic->GetNumberOfEstimatedMuons();  // Info: GAP GAP2014_037.
        station.mask_signal = false;

        const double vemCut = (!theRec.IsTimeRec ? UnivRecNS::vemSignalCut : UnivRecNS::vemTimeCut[2]);
        if (station.signalsize > vemCut)
          recinfo.nStationsMD++;

        theRec.fstation.push_back(station);

        printf("id=%d %lf %lf %lf Area=%lf Nmu=%lf \n",station.id,station.xg/1.e2,station.yg/1.e2,station.zg/1.e2,station.DetectorArea/1.e4,station.signalsize);
      } // end counter loop

    } // has MEvent

    //-----------------------------------------------
    // Set station vector  ( ASCII )
    //-----------------------------------------------
    /*loop
    {
    itype=1;
    double vemCut=( !theRec.IsTimeRec ? UnivRecNS::vemSignalCut : UnivRecNS::vemTimeCut[2] );
    if (  station.signalsize >  vemCut ) recinfo.nStationsScin++;
    */

    //-----------------------------------------------
    // Selection according to the number of stations
    //-----------------------------------------------

    theRec.nRecStations = theRec.fstation.size();

    printf("Number of stations loaded: %d (%d/%d/%d) (id=%d)  | log(E_SD)=%5.2lf\n",
           theRec.nRecStations, recinfo.nStationsWCD, recinfo.nStationsScin, recinfo.nStationsMD, recinfo.id, augerinfo.logE_SD);

    const int nMin = !theRec.IsTimeRec ? UnivRecNS::nSignalStationsMin : UnivRecNS::nTimeStationsMin;
    const bool okStations = (recinfo.nStationsWCD >= nMin);
    if (!okStations) {
      printf("Not processing: Number of stations loaded %d | WCD above threshold %d  \n", theRec.nRecStations, recinfo.nStationsWCD);
      return false;
    }

    //-----------------------------------------------
    // Set AMIGA/ASCII mask_signal according to WCD
    //-----------------------------------------------
    for (int idet_i = 0; idet_i < theRec.nRecStations; ++idet_i) {
      UnivRecNS::RecStation_t* st = &theRec.fstation[idet_i];
      if (st->type == eWCD)
        continue;
      st->mask_signal = false;

      int idet_wcd = UNDEF;
      for (int idet_j = 0; idet_j < theRec.nRecStations; ++idet_j)
        if (theRec.fstation[idet_j].type == eWCD && theRec.fstation[idet_j].id == st->id)
          idet_wcd = idet_j;

      if (idet_wcd == UNDEF)
        continue;
      st->mask_signal = theRec.fstation[idet_wcd].mask_signal;
    }

    //-----------------------------------------------
    // Hot station, Height of the core position (WCD)
    //-----------------------------------------------

    double vemMax = -1e10;
    for (int idet = 0; idet < theRec.nRecStations; ++idet)
      if (theRec.fstation[idet].type == eWCD && theRec.fstation[idet].signalsize > vemMax) {
        vemMax = theRec.fstation[idet].signalsize;
        recinfo.gDetHot = idet;
      }
    recinfo.zgcore = theRec.fstation[recinfo.gDetHot].zg;

    return true;
  }


  void
  UniversalityFitter::RunBarilocheReconstruction(evt::Event& event)
  {
    const int verbose = fBarilocheConfig.verbose;
    const int debug = fBarilocheConfig.debug;
    fBarilocheConfig.IsData = (!event.HasSimShower());

    UnivRecNS::UnivRec& theRec = *gUnivRecNS;
    if (!theRec.InitRec(fBarilocheConfig.IsData, fBarilocheConfig.RecType,
                        fCalibOpt, fBarilocheConfig.RecDetector,
                        fBarilocheConfig.RecSys, fRecMixture))
      return;

    if (!InitBarilocheReconstruction(event))
      return;

    if (debug)
      cerr << "InitEvent ok" << endl;
    if (verbose > 1)
      cerr << "Number of stations loaded: "
           << theRec.GetnRecStations()
           << " (" << theRec.GetRecInfo()->nStationsWCD
           << ")" << endl;

    if (!theRec.InitRecParameters())
      return;

    if (debug)
      cerr << "InitRecParameters ok" << endl;

    //------- Fit stage -1/0
    if (theRec.SetFitStage(-1))
      Rec(theRec);

    if (!theRec.StationSelection())
      return;

    if (debug)
      cerr << "StationSelection ok" << endl;

    if (theRec.SetFitStage(0)) {
      Rec(theRec);
      if (theRec.GetRecInfo()->FitStatus[0] != 3)
        Rec(theRec);
    }

    if (debug)
      theRec.PrintRecInfo(true);

    //--------- Fit stage 1  :   Time Likelihood, Gaus PDF, (Xmax,T0)  estimation
    if (theRec.SetFitStage(1)) {
      double lnP_max = -1e10;
      double Xmax_init[4] = { 650, 750, 850, 950 };
      double Xmax0 = UNDEF;
      for (int iscan = 0; iscan < 4; ++iscan) {
        theRec.GetRecInfo()->Xmax = Xmax_init[iscan];
        Rec(theRec);
        if (debug)
          cerr << "SetFitStage 1 Search Xmax Rec ok" << endl;
        if (theRec.GetRecInfo()->lnP[0] > lnP_max) {
          lnP_max = theRec.GetRecInfo()->lnP[0];
          Xmax0 = Xmax_init[iscan];
        }
      }
      theRec.GetRecInfo()->Xmax = Xmax0;
      Rec(theRec);
      if (debug)
        theRec.PrintRecInfo(true);
    }

    //--------- Fit stage 2  ( Signal+Time+Constrains, Binomial PDF  ) ,
    if (theRec.SetFitStage(2))  {   // Set GausPDF
      Rec(theRec);
      if (!theRec.IsGoodRec())
        Rec(theRec);
      if (theRec.RejectTimeOutliers())
        Rec(theRec);
      if (!theRec.IsGoodRec())
        Rec(theRec);
      theRec.PrintRecInfo(true);
    }

    //bool WriteStations=true;
    //theRec.RWRecinfo(fdataOut,WriteStations,false);

    fittedPars.energy = pow(10., theRec.GetRecInfo()->logE);
    fittedPars.energyUnc = pow(10., theRec.GetRecInfo()->logE_err);

    fittedPars.xmax = theRec.GetRecInfo()->Xmax * gpercm2;
    fittedPars.xmaxUnc = theRec.GetRecInfo()->Xmax_err * gpercm2;

    fittedPars.Nmu = theRec.GetRecInfo()->Nmu;
    fittedPars.NmuUnc = theRec.GetRecInfo()->Nmu_err;

    CoordinateSystemPtr hottestStationCS = LocalCoordinateSystem::Create(fHottestStation->position);
    fittedPars.axis = Vector(1, theRec.GetRecInfo()->theta, theRec.GetRecInfo()->azi, hottestStationCS, Vector::kSpherical);
    fittedPars.axisUncTheta = 0.; // TODO not implemented
    fittedPars.axisUncPhi = 0.; // TODO not implemented

    fittedPars.core = Point(theRec.GetRecInfo()->xgcore * cm, theRec.GetRecInfo()->ygcore * cm, theRec.GetRecInfo()->zgcore * cm - fHottestStationHeight, hottestStationCS);
    fittedPars.coreUncX = theRec.GetRecInfo()->xgcore_err * cm;
    fittedPars.coreUncY = theRec.GetRecInfo()->ygcore_err * cm;

    fittedPars.coretime = fHottestStation->startTime + TimeInterval(theRec.GetRecInfo()->T0);

    fittedPars.isGoodRec = theRec.IsGoodRec();
  }


  void
  UniversalityFitter::SetRecData(evt::Event& event)
  {
    if (!fittedPars.isGoodRec) {
      ostringstream msg;
      msg << "No reconstruction, not writing data to event.";
      INFO(msg);
      return;
    }

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

    ShowerUnivRecData& recShower = event.GetRecShower().GetUnivRecShower();

    CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(fittedPars.core);
    CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

    recShower.SetGoodRec(fittedPars.isGoodRec);

    recShower.SetEnergy(fittedPars.energy, fittedPars.energyUnc);
    recShower.SetNmu(fittedPars.Nmu, fittedPars.NmuUnc);
    recShower.SetXmax(fittedPars.xmax, fittedPars.xmaxUnc);

    recShower.SetTimeModelOffset(fittedPars.TimeModelOffset, fittedPars.TimeModelOffsetUnc);

    recShower.SetCorePosition(fittedPars.core);
    recShower.SetCoreError(Vector(fittedPars.coreUncX, fittedPars.coreUncY, 0, coreCS));
    recShower.SetCoreTime(fittedPars.coretime, 0);  // TODO coretime error not implemented yet

    recShower.SetAxis(fittedPars.axis);
    recShower.SetThetaError(fittedPars.axisUncTheta);
    recShower.SetPhiError(fittedPars.axisUncPhi);

    // TODO: add filling of these numbers when Bariloche reco is used
    recShower.SetNumberOfShapeCandidates(fittedPars.fNCandShape);
    recShower.SetNumberOfStartTimeCandidates(fittedPars.fNCandStartTime);
    recShower.SetNumberOfLDFCandidates(fittedPars.fNCandLDF);

    ostringstream msg;
    msg << "\nReconstructed shower parameters\n"
           "\n"
           "    energy       = " << recShower.GetEnergy() / EeV << " +/- " << recShower.GetEnergyError() / EeV << " [EeV]\n"
           "    Nmu          = " << recShower.GetNmu() << " +/- " << recShower.GetNmuError() << "\n"
           "    Xmax         = " << recShower.GetXmax() / gpercm2 << " +/- " << recShower.GetXmaxError() / gpercm2 << " [g/cm^2]\n"
           "\n"
           "    core         = (" << recShower.GetCorePosition().GetX(siteCS) / meter << ", "
        << recShower.GetCorePosition().GetY(siteCS) / meter << ", "
        << recShower.GetCorePosition().GetZ(siteCS) / meter << ") [m] (site) +/-\n"
           "                   (" << recShower.GetCoreError().GetX(siteCS) / meter << ", "
        << recShower.GetCoreError().GetY(siteCS) / meter << ", "
        << recShower.GetCoreError().GetZ(siteCS) / meter << ") [m]\n"
           "\n"
           "    theta        = " << recShower.GetAxis().GetTheta(coreCS) / degree << " +/- " << recShower.GetThetaError() / degree << " [deg] (core)\n"
           "    phi          = " << NormalizeAngleMinusPiPi(recShower.GetAxis().GetPhi(coreCS)) / degree << " +/- " << recShower.GetPhiError() / degree << " [deg] (core)\n";
    INFO(msg);
  }


  void
  UniversalityFitter::WriteRecParameters(evt::Event& event)
  {
    const SEvent& sEvent = event.GetSEvent();
    const int sdEventId = sEvent.GetHeader().GetId();

    const CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

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

    ofstream* textOut = 0;
    if (fActiveMethod == "Karlsruhe")
      textOut = &fitInfoKG;
    else if (fActiveMethod == "Bariloche")
      textOut = &fitInfoBG;
    else {
      cerr << "unknown method" << endl;
      exit(1);
    }

    *textOut
      << format("%i ") % (fittedPars.isGoodRec)
      << format("%i ") % (sdEventId)
      << format("%i ") % (fEventHasSaturation)

      << format("%.5f ") % (fittedPars.energy / EeV)
      << format("%.5f ") % (fittedPars.xmax / gpercm2)
      << format("%.5f ") % (fittedPars.Nmu)

      << format("%.5f ") % (fittedPars.core.GetX(showerPlaneCS) / m)
      << format("%.5f ") % (fittedPars.core.GetY(showerPlaneCS) / m)

      << format("%.5f ") % (fittedPars.coretime.GetGPSSecond())
      << format("%.5f ") % (fittedPars.coretime.GetGPSNanoSecond())

      << format("%.5f ") % (fittedPars.axis.GetTheta(siteCS) / degree)
      << format("%.5f ") % (NormalizeAngleMinusPiPi(fittedPars.axis.GetPhi(siteCS)) / degree)
      << endl;
  }


  VModule::ResultFlag
  UniversalityFitter::Finish()
  {
    if (fWriteToTextFile) {
      fitInfoKG.close();
      fitInfoBG.close();
    }

    return eSuccess;
  }

}