UnivRec.cc

Go to the documentation of this file.

#include "UnivRec.h"

#include <TMath.h>
#include <TMatrixTSym.h>
#include <TMatrixDSym.h>
#include <TMatrixDSymEigen.h>
#include <TGraphErrors.h>
#include <TCanvas.h>

#include <TRandom3.h>

using namespace UnivRecNS;


//---------------------------------------------------------
//Electronics/GPS  settings (WCD/Scintillators)
//---------------------------------------------------------
//           (Upgrade)  Sampling=8  [ns], 12-bit, InterGPS Accuracy 4  [ns]
//           (Current)  Sampling=25 [ns], 10-bit, InterGPS Accuracy 10 [ns]
//
//---------------------------------------------------------
//      Saturation (WCD)
//---------------------------------------------------------
//              Anode non-linearity fo the current nominal gain starts at 100 [mA], 5 [V] at 50 Ohm , 1 nC in 10 [ns]
//                  Gain=2 x 10^5  <pe>=94  -> 1 [nC]= 6.25 10^9 electrons -> 3.12e4 photo-electrons -> ~300 [vem]
//
//              Electronic saturation :    5 anode [ADC] = 1 [vem]
//                       (Upgrade SDE) ->  4096/5= 800 [vem]
//                       (Current SDE) ->  1024/5= 200 [vem]
//
//               But, the small PMT is now part of the baseline SDE design. Saturation should happen at total signal ~20000 [vem] ( ~2000 [vem]/8 [ns] at peak )
//
//               Saturation (Total signal) :  1000/20000 [vem]   Current/Upgraded SDE
//
//---------------------------------------------------------
//         Saturation (ASCII)
//---------------------------------------------------------
//             Set to electronic saturation assuming 5 anode [ADC] = 1 [vem]
//
//                      (Upgrade SDE) ->  4096/5 =  800 [vem]
//                      (Current SDE) ->  1024/5 =  200 [vem]
//
//             For ASCII there is a shaper so the time response is matched to the one for the WCD.
//
//             Saturation (Total signal) :  1500/10000 [vem]   Current/Upgraded SDE


//---------------------------------------------------------
//         Saturation (Muon detectors)
//---------------------------------------------------------
//              (0-64 x 3 muons)

//--------------------------------------
//         Rec types
//--------------------------------------
//
//       Station selection for Signal/Time likelihood:     after Fitstage -1
//       Gaus PDF flags                              :     after Fitstage  1
//
//
//
//       Hybrid reconstructions   ( direction fixed to hybrid reconstruction )
//       -----------------------
//                Fitstage -1 :  (logE,Xmax)=(log_FD,Xmax_FD)                 ->   Nmu_FD             (core search for station selection)
//                Fitstage  0 :  (logE,Xmax) constrained to (logE_FD,Xmax_FD) ->   Nmu_FD
//
//         [0]        Nmu_FD
//         [7]        Nmu_FD and OffsetM_Mu (FitStage 1 and 2 also )
//
//
//       SD  reconstructions
//       -----------------------
//                 Fitstage -1 : (Nmu,Xmax) fixed to <Nmu>(E,theta),<Xmax>(E)  (core search for station selection)
//                 Fitstage  0 : (Nmu,Xmax) fixed to <Nmu>(E,theta),<Xmax>(E)
//                 Fitstage  1 : (core,logE,Nmu,direction) fixed from FitStage 0                      ( Xmax/T0 guess )
//                 Fitstage  2 : Nmu=<Nmu>(E,theta,Xmax)
//
//         [1]    OffsetM_Mu   fit            Nmu=<Nmu>(logE,theta,Xmax)      ( Xmax/OffsetM_mu ) free parameter
//         [2]    Rec. Stops at FitStage 0
//         [3]    Rec. Stops at FitStage 0,    Nmu free in FitStage 0,        (  Nmu_SD with <Xmax> )   only WCD+TopScin/WCD+MD
//         [4]    Xmax_SD        fit           Nmu=<Nmu>(logE,theta,Xmax)
//         [5]    Xmax_SD        fit           log(E)=log(E_ref)  Nmu free
//         [8]    Xmax_SD        fit           Nmu=<Nmu>(logE,theta,-100)
//
//       Photon search
//       -----------------
//         [6]  Same as [4] but  Nmu=Nmu_PhotonSearch , Xmax=Xmax_PhotonSearch
//                 (FitStage=-1/0 Signal Likelihood, FitStage=1/2 Only Time likelihood )
//
//
//--------------------------------------
//         Data sets     ( always 6T5 )
//--------------------------------------
//          DataSet=0    SD            ( 1500 [m] log(E_ADST)>18.6 [eV] , 750 [m]  log(E_ADST)>17.8 [eV] )
//          DataSet=1    Golden        ( log(E_FD)>17.5 [eV] )
//          DataSet=2    SD low energy ( 6T5 750 [m] events that are also 6T5 1500 [m] events)
//          DataSet=3    SD high energy( log(E_ADST)>19.5 [eV] )

//--------------------------------------
//    Likelihood setting
//--------------------------------------
//             Parameters (logE,core,direction,Xmax,T0,Nmu)
//             WCD:          : Signal/Starttime/T1/T10/T50 likelihood used
//             Scintillators : Signal/Starttime/T1/T10/T50 likelihood used
//             Muon detectors: Signal/Starttime/T1/T10/T50 likelihood used
//             At least 4 stations with r<%4.0lf [m] after core search are required for the time likelihood \n",rQuantileCut[1]/1.e2);

//-------------------------------------------------------------

//-------------------------------------------------------------------
//-------------------------------------------------------------------
// Set flags
// Init Univ. Parameterizations
// Create atmosphere
//-------------------------------------------------------------------
//-------------------------------------------------------------------
UnivRec::UnivRec()
{
  fatm = new AtmosphereNS::Atmosphere();
  for (int itype = 0; itype < 3; itype++) {
    ftime[itype] = new UnivParamTimeNS::UnivParamTime(itype);
    fsignal[itype] = new UnivParamNS::UnivParam(itype);
  }
  fcalib = NULL;
}

bool UnivRec::InitRec(bool IsData,  int RecType, int CalibOpt, int RecDetector /*all*/, int RecSys /*Only CalibOpt=-2*/ , int RecMixture /*Only CalibOpt=0,1*/)
{
  //General
  gRecMC = !IsData;
  gRecType = RecType;
  if (gRecType < 0 || gRecType > 8) {
    printf("ERROR: Wrong [RecType] \n");
    return false;
  }

  gCalibOpt = CalibOpt;
  if (gCalibOpt < -3 || gCalibOpt > 3) {
    printf("ERROR: Wrong [CalibOpt] %d \n", gCalibOpt);
    return false;
  }

  if (RecDetector < 0 || RecDetector > 3) {
    printf("Wrong [RecDetector] \n");
    return false;
  }
  if (RecDetector > 0 && (gRecType == 1 || gRecType == 4 || gRecType == 5 || gRecType == 6 || gRecType == 7 || gRecType == 8)) {
    printf("ERROR: Time model for Scintillator above below ground is not tested, you can delete this break at your own risk \n");
    return false;
  }
  if (gRecType == 3  && RecDetector == 0) {
    printf("ERROR: Nmu_SD only possible with upgrades \n");
    return false;
  }

  gDetectorTypeMask[0] = false;
  gDetectorTypeMask[1] = (RecDetector == 1 || RecDetector == 3                   ? false : true);
  gDetectorTypeMask[2] = (RecDetector == 2 || RecDetector == 3                   ? false : true);

  //Only CalibOpt=-2
  gRecSys = RecSys;
  if (CalibOpt > -2 &&   gRecSys != 0) {
    printf("ERROR: Systematics only implemented for CalibOpt=-2 \n");
    return false;
  }
  if (CalibOpt <= -2 && (gRecSys < 0 || gRecSys >= UnivCalibConstantsNS::nSys)) {
    printf("ERROR: Wrong [RecSys] \n");
    return false;
  }

  //Only CalibOpt=0,1
  gRecMixture = RecMixture;
  if (CalibOpt >= 0 && (gRecMixture < 0 || gRecMixture > UnivCalibConstantsNS::nRecMixtures)) {
    printf("ERROR: Wrong [gRecMixture] \n");
    return false;
  }

  if (gRecType == 6 && CalibOpt != -1)
    printf("WARNING: Photon search without photon calib parameters \n");
  //--
  IsTimeRec = (gRecType == 1 || gRecType == 4 || gRecType == 5 || gRecType == 6 || gRecType == 8 || gRecType == 7 ? true : false);

  if (fcalib != NULL) delete fcalib;
  fcalib = new UnivCalibConstantsNS::UnivCalibConstants(gCalibOpt, gRecSys,  gRecMixture);

  //--Set Golden rec. flag ( So FD rec. parameters are used )
  IsGoldenRec = (gRecType == 0 || gRecType == 7 ? true : false);
  if (IsGoldenRec && CalibOpt != -3)  {
    printf("ERROR: Wrong CalibOpt  \n");
    return false;
  }

  //-- Only required for Bariloche framework (Initialized in a specialized routine)
  gRecInfill = false;
  gRecSDE_Upgrade = false;
  gRecAoP = UNDEF;
  gRecDataSet = (IsGoldenRec ? 1 : UNDEF);

  return true;
}


bool
#ifndef ExternalUse
UnivRec::InitRec_InternalFramework(bool RecInfill /*all*/, int RecDataSet, int year, std::string& FileName /*Only data*/, bool RecSDEUpgrade, int RecAoP /*Only MC */)
#else
UnivRec::InitRec_InternalFramework(bool RecInfill /*all*/, int RecDataSet, int /*year*/, std::string& /*FileName*/ /*Only data*/, bool RecSDEUpgrade, int RecAoP /*Only MC */)
#endif
{
  gRecInfill = RecInfill;
  if (gRecMC) {
    gRecSDE_Upgrade = RecSDEUpgrade;
    gRecAoP = RecAoP;
    if (!(gRecAoP == 0 || gRecAoP == 2 || gRecAoP == 3)) {
      printf("Wrong [RecAoP]\n");
      return false;
    }
  } else {
    gRecDataSet = RecDataSet;
    if (!(gRecDataSet == 0 || gRecDataSet == 1 || gRecDataSet == 2 || gRecDataSet == 3)) {
      printf("Wrong [RecDataSet]\n");
      return false;
    }

    if ((gRecDataSet == 0 || gRecDataSet == 2 || gRecDataSet == 3) && IsGoldenRec)
      printf("Hybrid reconstruction on SD data not possible\n");
    if ((gRecDataSet == 2 || gRecDataSet == 3) && gRecInfill)
      printf("Wrong data set for Infill rec.\n");

#ifndef ExternalUse
    if (gRecDataSet == 0 && !gRecInfill && (year < 2004 || year > 2012)) {
      printf("Wrong SD data set\n");
      return false;
    }
    if (gRecDataSet == 0)
      FileName = ToyUtilsNS::GetDataFileName(gRecInfill ? 1 : 0 , year);
    else if (gRecDataSet == 1)
      FileName = ToyUtilsNS::GetDataFileName(gRecInfill ? 3 : 2,  year);
    else if (gRecDataSet == 2)
      FileName = ToyUtilsNS::GetDataFileName(4,  year);
    else if (gRecDataSet == 3)
      FileName = ToyUtilsNS::GetDataFileName(5,  year);
#endif
  }

  return true;
}



//---------------------------------------------------
//---------------------------------------------------
// Init/Clear event/station
//---------------------------------------------------
//---------------------------------------------------
void UnivRec::InitRecStation(RecStation* station)
{
  station->type = UNDEF;
  station->id = UNDEF;

  station->AreaOverPeak = -1.;
  station->DetectorArea = UNDEF;
  station->isSaturated = false;

  station->xg = 0.;
  station->yg = 0.;
  station->zg = 0.;

  station->GPSnano = 0.;

  station->r = 0.;
  station->psi = 0.;
  station->T0_PF = 0.;
  for (int icomp = 0; icomp < 4; icomp++) {
    station->DX[icomp] = 0.;
    station->DL[icomp] = 0.;
    station->T0_CF[icomp] = 0.;
  }

  station->mask_quantiles = true;
  for (int iq = 0; iq < UnivRecNS::nQ; iq++) {
    station->UseGausPDF_quantiles[iq] = true;
    station->quantile[iq] = UNDEF;
    station->tquantile[0][iq] = UNDEF;
    station->tquantile[1][iq] = UNDEF;
    station->tq_rec[iq] = 0.;
    station->tq_pred[iq] = 0.;
    station->tq_rms[iq] = 0.;
    station->tq_corr[iq] = 0.;
    station->lnP_quantiles = 0.;
  }

  station->UseGausPDF_start = true;
  station->mask_starttime = true;
  station->starttime = UNDEF;
  station->start_rec = 0.;
  station->start_pred = 0.;
  station->start_rms = 0.;
  station->start_corr = 0.;
  station->lnP_start = 0.;

  station->density = 0.;
  station->signalsize = 0.;
  station->mask_signal = true;
}

void UnivRec::ClearRecEvent()
{

  //RecInfo
  recinfo.id = UNDEF;
  recinfo.AugerDensity = UNDEF;
  recinfo.iatm = UNDEF;

  recinfo.xgcore = UNDEF;
  recinfo.ygcore = UNDEF;
  recinfo.xgcore_err = UNDEF;
  recinfo.ygcore_err = UNDEF;
  recinfo.zgcore = UNDEF;

  recinfo.logE = UNDEF;
  recinfo.logE_err = 0.;
  recinfo.Xmax = UNDEF;
  recinfo.Xmax_err = 0.;
  recinfo.Nmu = UNDEF;
  recinfo.Nmu_err = 0.;
  recinfo.theta = UNDEF;
  recinfo.azi = UNDEF;
  recinfo.T0 = UNDEF;
  recinfo.T0_err = 0.;
  recinfo.logE_ref = UNDEF;
  recinfo.logE_ref_err = 0.;
  recinfo.theta_ref = UNDEF;
  recinfo.azi_ref = UNDEF;

  recinfo.OffsetM_Mu = UNDEF;

  recinfo.RA = UNDEF;
  recinfo.Dec = UNDEF;

  recinfo.nStationsMD = 0;
  recinfo.nStationsScin = 0;
  recinfo.nStationsWCD = 0;
  recinfo.nTimeStations = 0;
  recinfo.nSignalStations = 0;

  for (int iloop = 0; iloop < 2; iloop++)
    recinfo.FitStatus[iloop] = -1;

  for (int iloop = 0; iloop < 3; iloop++) {
    recinfo.lnP[iloop] = 0;
    recinfo.seed[iloop] = 0;
  }
  recinfo.MinEigenVal = 0.;
  recinfo.edm = 0.;

  recinfo.FitStage = -2;

  recinfo.AllTimesOK = false;
  for (int i = 0; i < 2; i++) {
    recinfo.TimeRejected[i] = 0.;
    recinfo.TimeMaxdev[i] = 0.;
  }

  recinfo.rHot = UNDEF;
  recinfo.gDetHot = UNDEF;
  recinfo.nFCNcalls = 0;
  recinfo.IsSaturated[0] = false;
  recinfo.IsSaturated[1] = false;
  recinfo.IsSaturated[2] = false;

  // MC Info
  MCinfo.im = UNDEF;
  MCinfo.ip = UNDEF;
  MCinfo.ie = UNDEF;
  MCinfo.ith = UNDEF;
  MCinfo.ish = UNDEF;
  MCinfo.iatm = UNDEF;
  MCinfo.xgcore = UNDEF;
  MCinfo.ygcore = UNDEF;
  MCinfo.zgcore = UNDEF;
  MCinfo.theta = UNDEF;
  MCinfo.azi = UNDEF;
  MCinfo.logE = UNDEF;
  MCinfo.Xmax = UNDEF;
  MCinfo.Nmu = UNDEF;
  MCinfo.T0 = UNDEF;

  // Auger Info
  augerinfo.xgcore_SD = UNDEF;
  augerinfo.ygcore_SD = UNDEF;
  augerinfo.zgcore_SD = UNDEF;
  augerinfo.theta_SD = UNDEF;
  augerinfo.azi_SD = UNDEF;
  augerinfo.logE_SD = UNDEF;
  augerinfo.logE_SD_err = 0.;

  augerinfo.xgcore_FD = UNDEF;
  augerinfo.ygcore_FD = UNDEF;
  augerinfo.zgcore_FD = UNDEF;
  augerinfo.theta_FD = UNDEF;
  augerinfo.azi_FD = UNDEF;
  augerinfo.logE_FD = UNDEF;
  augerinfo.logE_FD_err = 0.;
  augerinfo.Xmax_FD = UNDEF;
  augerinfo.Xmax_FD_err = 0.;

  augerinfo.isGoodSD = false;
  augerinfo.isGoodFD = false;

  augerinfo.RA = UNDEF;
  augerinfo.dec = UNDEF;

  // Stations
  fstation.clear();


}

#ifndef ExternalUse

bool UnivRec::InitEvent(ToyShowerNS::ToyShower* fevt)
{
  ClearRecEvent();

  nRecStations = 0;
  XmaxLow = (gRecType == 6 ? 200. : low[4]);
  XmaxHigh = high[4];

  bool okT4 = true, okEnergyCut = true, okAtm = true, okXmax = true;

  //-----------------------------------------------
  //  Set Id, Check T4,  Init Random Number
  //-----------------------------------------------

  recinfo.id = (gRecMC ?  int(fevt->GetCorsikaInfo()->runnr) : fevt->GetAugerInfo()->sdid);
  okT4 = fevt->CheckT4(gRecInfill);
  TRandom3* rnd = new TRandom3(recinfo.id + gRecType);

  //-----------------------------------------------
  //  Load MC info
  //-----------------------------------------------
  if (gRecMC) {
    ToyUtilsNS::GetShowerIndexes(recinfo.id, MCinfo.im, MCinfo.ip, MCinfo.ie, MCinfo.ith, MCinfo.ish, MCinfo.iatm);

    MCinfo.theta = fevt->GetCorsikaInfo()->theta;
    MCinfo.azi = fevt->GetCorsikaInfo()->azi;

    float x_ref, y_ref, z_ref;
    fevt->GetCoreMC(x_ref, y_ref, z_ref);
    MCinfo.xgcore = x_ref;
    MCinfo.ygcore = y_ref;
    MCinfo.zgcore = z_ref;

    MCinfo.T0 = fevt->GetCoreT0_MC();

    MCinfo.logE = log10(fevt->GetCorsikaInfo()->E0) + 9.;
    MCinfo.Xmax = fevt->GetCorsikaInfo()->Xmax + ToyUtilsNS::XmaxCorsikaOffset;

    int ir0 = 4 , ipsi0 = 2;
    ToyShowerNS::SamplingArea_t* saRef = fevt->GetDenseSA(ir0, ipsi0);
    double density = fevt->GetAtmosphere()->GetDensity(saRef->zg);
    double sMu = saRef->Muon.vem;

    double sMuRef = fsignal[0]->GetSignal(ToyUtilsNS::rRing[ir0], MCinfo.Xmax, MCinfo.logE, MCinfo.theta, ToyUtilsNS::psiS[ipsi0],
                                          density, ToyShowerNS::hgroundRef, 1.0 , 0 , MCinfo.iatm);
    MCinfo.Nmu = sMu / sMuRef;

  }

  //-----------------------------------------------
  //Set : Atmosphere/Ground Density
  //-----------------------------------------------
  if (gRecMC) {
    GPSsec = 0;

    int imonth = fevt->GetAtmosphere()->GetCurrentAtmosphereMonth();
    fatm->SetCurrentAtmosphere(imonth);
    if (imonth != MCinfo.iatm) okAtm = false;
    recinfo.iatm = imonth;
  } else {
    ToyShowerNS::AugerInfo_t* ainfo = fevt->GetAugerInfo();
    if (ainfo->p == 0. || ainfo->T == 0.)
      recinfo.AugerDensity = 1.06 * 1.e-3; // Average density
    else {
      double R = 286.9; //  [atm][l]/[mol][K]
      recinfo.AugerDensity = ainfo->p * 100000. / 286.9 / ainfo->T / 1.e3; // [g]/[cm3]
    }
    GPSsec = ainfo->GPSsec;

    unsigned int nsec = 0;
    utl::TimeStamp time_stamp(ainfo->GPSsec, nsec);
    int imonth = time_stamp.GetMonth();
    fatm->SetCurrentAtmosphere(imonth);
    recinfo.iatm = imonth;
  }


  //-----------------------------------------------
  //Set Auger Info
  //-----------------------------------------------
  augerinfo.Xmax_FD_err = 20.;
  augerinfo.logE_FD_err = log10(1.1);

  if (gRecMC) {
    //--- Cut events with Xmax_ref below ground ( so Mean/RMS in ToyUtils is consistent )
    if (std::isnan(MCinfo.Xmax) || std::isinf(MCinfo.Xmax) || MCinfo.Xmax <= 0.)                     okXmax = false;

    //--- Fluctuated values of Energy and Xmax
    augerinfo.Xmax_FD = rnd->Gaus(MCinfo.Xmax, augerinfo.Xmax_FD_err);
    augerinfo.logE_FD = MCinfo.logE + rnd->Gaus(0., augerinfo.logE_FD_err);

    //---Fluctuated values of (theta,azi)
    TVector3 theCR;
    theCR.SetMagThetaPhi(1., MCinfo.theta, MCinfo.azi);
    TVector3 rndD = theCR.Orthogonal();
    rndD.SetMag(fabs(rnd->Gaus(0, SmearingAngle * TMath::DegToRad())));
    rndD.Rotate(rnd->Rndm()*TMath::TwoPi(), theCR);
    rndD = theCR + rndD;
    theCR = rndD.Unit();
    augerinfo.theta_FD = theCR.Theta();
    augerinfo.azi_FD = NormalizeAngle(theCR.Phi());
  } else {
    ToyShowerNS::AugerInfo_t* ainfo = fevt->GetAugerInfo();
    // SD (Observer reconstruction)    core position is not saved (Do it in Offline!)
    augerinfo.logE_SD = log10(ainfo->E_SD) + 18.;
    augerinfo.logE_SD_err = log10(1. + ainfo->E_SD_err / ainfo->E_SD);
    if (gCalibOpt <= -2) augerinfo.logE_SD += log10(UnivCalibConstantsNS::fEnergySys_Auger[gRecSys]);
    augerinfo.theta_SD = ainfo->thetaObserver;
    augerinfo.azi_SD = ainfo->aziObserver;
    augerinfo.isGoodSD = fevt->isGoodSD(true);

    recinfo.RA = ainfo->RA;
    recinfo.Dec = ainfo->Dec;

    // FD
    int EyeSel = (gRecSys > 8  ? gRecSys - 8 : UNDEF);

    float Xmax_FD, Xmax_FD_err, E_FD, E_FD_err, theta_FD, theta_FD_err, azi_FD, azi_FD_err;
    augerinfo.isGoodFD = fevt->isGoodFD(Xmax_FD, Xmax_FD_err, E_FD, E_FD_err, theta_FD, theta_FD_err, azi_FD, azi_FD_err , true, EyeSel);

    if (augerinfo.isGoodFD) {
      augerinfo.logE_FD = log10(E_FD) + 18.; // If Offline is using M. Tueros et al, this is correct (otherwise I have to do the parameterization of missing energy)
      augerinfo.Xmax_FD = Xmax_FD;

      augerinfo.theta_FD = theta_FD;
      augerinfo.azi_FD = azi_FD;
    }
    if (gRecDataSet == 1 || gRecType == 0) okXmax = augerinfo.isGoodFD;

    // Energy cuts for different rec
    okEnergyCut = CheckEnergyCut(true);
  }
  if (gCalibOpt <= -2) {
    augerinfo.logE_FD += log10(UnivCalibConstantsNS::fEnergySys_Auger[gRecSys]);
    augerinfo.Xmax_FD += UnivCalibConstantsNS::fXmaxSys_Auger[gRecSys];
  }

  //-----------------------------------------------
  // Selection
  //-----------------------------------------------
  if (!okT4)         {
    printf(" Event is not T4  %d \n", recinfo.id);
    return false;
  }
  if (!okEnergyCut)  {
    printf(" Energy below analysis threshold  %d \n", recinfo.id);
    return false;
  }
  if (!okAtm)        {
    printf("Error: inconsistency in the atmospheres %d \n", recinfo.id);
    return false;
  }
  if (!okXmax)      {
    printf(" Xmax not well defined %d \n", recinfo.id);
    return false;
  }
  if (!gRecMC) {
    if (!augerinfo.isGoodSD) {
      printf("Warning: SD event not 6T5 \n");
      return false;
    }
    if (!augerinfo.isGoodFD && (IsGoldenRec || gRecDataSet == 1)) {
      printf("Error: Golden event without FD rec parameters %d \n", recinfo.id);
      return false;
    }
  }
  if (recinfo.id == 4067414)  {
    printf("Lightening event %d \n", recinfo.id);
    return false;
  }

  //-----------------------------------------------
  // Set station vector
  //-----------------------------------------------

  int nSignalStations = 0;
  for (int idet = 0; idet < ToyUtilsNS::nDetectors; idet++)
    for (int itype = 0; itype < 3; itype++) {
      if (gDetectorTypeMask[itype]) continue;

      //----Station/Event trigger and selection given by WCD
      ToyShowerNS::SamplingArea_t* sa = fevt->GetArraySA(idet);
      if (!gRecInfill && !sa->isStandard) continue;
      if (!sa->isAlive      && sa->triggerFlag == 0) continue; // Accept all stations with trigger (even if they are thought to be dead)
      if (sa->isAccidental && sa->triggerFlag > 0) continue;

      //---- 2012 (No T2Life)   Events picked up when comparing energies with the Observer ( Warning: there might be more in 2012 )
      if (recinfo.id == 13296659 && sa->id == 203) continue;
      if (recinfo.id == 13769235 && sa->id == 131) continue;
      if (recinfo.id == 8270184 &&  sa->id == 294) continue;

      //---- Photon search (event with an accidental muon in one station)
      if (recinfo.id == 7728630 && sa->id == 752) continue;

      //---- Xmax SD outliers ( Xmax<500. )
      if (recinfo.id == 648728    && sa->id == 157) continue;
      if (recinfo.id == 1711230   && sa->id == 216) continue;
      if (recinfo.id == 10961377  && sa->id == 438) continue;

      RecStation_t station;
      InitRecStation(&station);

      station.type = itype;
      station.id = (gRecMC ? idet : sa->id);

      station.DetectorArea = (itype == 2 ? nMD * MDlength* MDwidth : UNDEF);

      station.xg = sa->xg;
      station.yg = sa->yg;
      station.zg = sa->zg;

      station.GPSnano = sa->GPSnano;

      if (sa->triggerFlag == 0)  {
        fstation.push_back(station);
        continue;
      }

      bool hasTrace = true;

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

      //---- Area over peak
      if (itype < 2) {
        if (gRecMC) { // MC
          if (itype == 0) { //WCD
            if (gRecAoP == 2 || gRecAoP == 3)      station.AreaOverPeak = AreaOverPeak_Run[gRecAoP];
            else                                 station.AreaOverPeak = -1.;
          } else //Scin
            station.AreaOverPeak = -1.;
        } else   station.AreaOverPeak = sa->AreaOverPeak; //Auger data
      }
      if (station.AreaOverPeak > 5.) continue;

      //---- Traces
      float TimeUnit = (!gRecSDE_Upgrade || itype == 2  ? 25. : 8.);
      int nTimeBins = (!gRecSDE_Upgrade || itype == 2  ? ToyShowerNS::nFADC : ToyShowerNS::nTimeBins / 4);
      int nWindow = (!gRecSDE_Upgrade || itype == 2  ? 1  : 3);
      float* vemtrace;
      double vemSat;
      int start, end;

      if (itype == 0)      {
        vemtrace = (gRecSDE_Upgrade ? sa->vemtrace_8ns : sa->vemtrace);
        vemSat = SaturationLevelVEM_WCD[gRecSDE_Upgrade];
      } else if (itype == 1) {
        vemtrace = (gRecSDE_Upgrade ? sa->vemtraceScin_8ns : sa->vemtraceScin);
        vemSat = SaturationLevelVEM_Scin[gRecSDE_Upgrade];
      } else                 {
        vemtrace = sa->vemtraceMD;
        vemSat = SaturationLevelMD;
      }

      //----Avg quantiles ( only MC )
      if (gRecMC) {
        station.tquantile[1][0] = UNDEF;
        if (itype == 0)      {
          station.tquantile[1][1] = sa->t10T;
          station.tquantile[1][2] = sa->t50T;
          station.tquantile[1][3] = sa->t90T;
        } else if (itype == 1) {
          station.tquantile[1][1] = sa->t10T_Scin;
          station.tquantile[1][2] = sa->t50T_Scin;
          station.tquantile[1][3] = sa->t90T_Scin;
        } else if (itype == 2) {
          station.tquantile[1][1] = sa->t10T_MD;
          station.tquantile[1][2] = sa->t50T_MD;
          station.tquantile[1][3] = sa->t90T_MD;
        }
      }

      //----Saturation/Signal
      //Data
      if (!gRecMC) {
        start = sa->StartBin;
        end = sa->EndBin;
        station.isSaturated = sa->isSaturated;
        station.signalsize = 0.;
        for (int itime = start; itime < end; itime++)
          station.signalsize += vemtrace[itime];
      }
      //MC
      else {
        start = 0;
        end = (itype == 2 ? nTimeBins - 3 : nTimeBins) ;
        station.isSaturated = false;
        station.signalsize = 0.;

        int itime = start;
        while (true) {
          //AMIGA
          if (itype == 2) {
            if (vemtrace[itime] > vemSat) {
              station.isSaturated = true;
              vemtrace[itime] = vemSat;
              station.signalsize += vemtrace[itime];
            } else station.signalsize += vemtrace[itime];
          }
          //WCD/ASCII
          else {
            if (vemtrace[itime] > vemSat) {
              station.isSaturated = true;
              vemtrace[itime] = vemSat;
              station.signalsize += vemSat;
            } else                              station.signalsize += vemtrace[itime];
          }

          itime++;
          if (itime >= end) break;
        }
      }

      //---- Data events with no traces
      if (!gRecMC && station.signalsize == 0. && sa->vem[0] > 0.) {
        hasTrace = false;
        station.signalsize = sa->vem[0];
        station.isSaturated = sa->isSaturated;
      }

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

      //----Start time
      //Data
      if (!gRecMC) {
        int StartTimeBin = sa->StartTimeBin;
        station.starttime = (double(StartTimeBin) + 0.5) * TimeUnit;
      }
      //MC
      else {
        int StartTimeBin;
        // AMIGA ( First muon )
        if (itype == 2) {
          for (int itime = 0; itime < nTimeBins; itime++)
            if (vemtrace[itime] > 0.) {
              StartTimeBin = itime;
              break;
            }
        }
        // WCD/ASCII ( Signal over a window )
        else {
          for (int itime = 0; itime < nTimeBins - nWindow; itime++) {
            if (vemtrace[itime] <= (0.1 / double(nWindow))) continue;
            double vemWindow = 0.;
            for (int it = itime; it < itime + nWindow; it++) vemWindow += vemtrace[it];
            if (vemWindow > 0.1) {
              StartTimeBin = itime;
              break;
            }
          }
        }
        station.starttime = (double(StartTimeBin) + 0.5) * TimeUnit;
      }

      //----Quantiles
      if (hasTrace) {
        station.quantile[0] = 0.01 * double(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 < nQ; iq++) {
          float tq, tq_err, signal = station.signalsize;
          ToyUtilsNS::GetTQ(&vemtrace[start], signal, end - start, TimeUnit,  station.quantile[iq], tq, tq_err);
          station.tquantile[0][iq] = tq + double(start) * TimeUnit;
        }
        for (int iq = 0; iq < nQ; iq++)
          if (station.signalsize < vemTimeCut[iq] || station.isSaturated) {
            station.quantile[iq] = UNDEF;
            station.tquantile[0][iq] = UNDEF;
            station.tquantile[1][iq] = UNDEF;
          }

      } else {
        station.quantile[0] = UNDEF;
        station.tquantile[0][0] = UNDEF;
        station.tquantile[1][0] = UNDEF;
        station.quantile[1] = UNDEF;
        station.tquantile[0][1] = UNDEF;
        station.tquantile[1][1] = UNDEF;
        station.quantile[2] = UNDEF;
        station.tquantile[0][2] = UNDEF;
        station.tquantile[1][2] = UNDEF;
        station.quantile[3] = UNDEF;
        station.tquantile[0][3] = UNDEF;
        station.tquantile[1][3] = UNDEF;
      }
      station.quantile[3] = UNDEF;
      station.tquantile[0][3] = UNDEF; // Disable T90

      //----

      if (UseStartTimeOnlyWhenSaturated && station.isSaturated)              station.mask_starttime = false;
      if (!UseStartTimeOnlyWhenSaturated && station.signalsize > vemStartCut) 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;


      //----GPS Jitter
      if (gRecMC) {
        double GPSjitter = rnd->Gaus(0., InterGPSAccuracy[gRecSDE_Upgrade]);
        station.GPSnano += GPSjitter;
      }
      //----

      if (IsTimeRec && !hasTrace) continue;
      double vemCut = (!IsTimeRec ? vemSignalCut : vemTimeCut[2]);
      if (itype == 0 && station.signalsize >  vemCut) recinfo.nStationsWCD++;
      if (itype == 1 && station.signalsize >  vemCut) recinfo.nStationsScin++;
      if (itype == 2 && station.signalsize >  vemCut) recinfo.nStationsMD++;

      fstation.push_back(station);
    } // Stations
  nRecStations = fstation.size();

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

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

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

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

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

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

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


  return true;
}

#endif

//-----------------------------------------------
//-----------------------------------------------
// Set initial values used in the reconstruction
//-----------------------------------------------
//-----------------------------------------------

bool UnivRec::InitRecParameters()
{
  //-----------------------------------------------
  // Set OffsetM_Mu
  //-----------------------------------------------
  recinfo.OffsetM_Mu = 0.;
  for (int itype = 0; itype < 3; itype++)
    ftime[itype]->SetOffsetM_Mu(recinfo.OffsetM_Mu);

  //-----------------------------------------------
  //--(Xmax,Nmu)
  //-----------------------------------------------
  recinfo.Xmax = 750.;
  recinfo.Nmu = 2.;

  //-----------------------------------------------
  // (theta,azi,core)
  //-----------------------------------------------

  //-- SD Auger standard reconstructions
#ifdef ExternalUse
  recinfo.theta = augerinfo.theta_SD;
  recinfo.azi = augerinfo.azi_SD;
  recinfo.xgcore = augerinfo.xgcore_SD;
  recinfo.ygcore = augerinfo.ygcore_SD;
  recinfo.zgcore = augerinfo.zgcore_SD;
#else
  //-- Own estimation ( needed for Bariloche MC )
  if (!GetRecEstimate(0)) return false;
  if (std::isnan(recinfo.theta)) return false;
#endif

  //-- Hybrid direction is use in all reconstruction stages
  if (IsGoldenRec) {
    recinfo.theta = augerinfo.theta_FD;
    recinfo.azi = augerinfo.azi_FD;
  }

  //-----------------------------------------------
  //--Xmax ranges
  //-----------------------------------------------
  XmaxLow = (gRecType == 6 ? 200. : 400.);
  XmaxHigh = 875. / cos(recinfo.theta) + 200.;

  //-----------------------------------------------
  //--Energy
  //-----------------------------------------------
  double E_sum = 0., w = 0.;
  int n_sum = 0;
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* station = &fstation[idet];
    if (station->type != 0) continue;
    if (station->isSaturated) continue;
    if (station->signalsize < vemSignalCut) continue;

    SetSPCoordinates(idet);
    if (station->r < 600.e2 || station->r > 2500.e2) continue;

    int iatm = fatm->GetCurrentAtmosphereMonth();
    double logEi = fsignal[0]->GetLogE(station->signalsize, station->r, recinfo.Xmax,  recinfo.theta, station->psi, station->density, station->zg, recinfo.Nmu, iatm);
    if (logEi < 18. || logEi > 20.5) continue;

    double wi = 1. / station->signalsize;
    E_sum += (pow(10., logEi) * wi);
    w += wi;
    n_sum++;
  }
  if (n_sum > 0)  recinfo.logE = log10(E_sum / w);
  else          recinfo.logE = 18.5;


  recinfo.theDir_Ref.SetMagThetaPhi(1., recinfo.theta, recinfo.azi);

  //-----------------------------------------------
  //T0
  //-----------------------------------------------
  SetT0FromHot(false, UNDEF);

  return true;
}

//---------------------------------------------------
//---------------------------------------------------
//  Shower plane coordinates
//---------------------------------------------------
//---------------------------------------------------

float UnivRec::DistPoints(float x1, float y1, float x2, float y2)
{
  return sqrtf((x2 - x1) * (x2 - x1)  + (y2 - y1) * (y2 - y1));
}


float UnivRec::NormalizeAngle(const float angle)
{
  if (angle > 2 * M_PI)
    return angle - int(angle / 2 / M_PI) * 2 * M_PI;
  if (angle < 0.)
    return 2 * M_PI + angle + int(-angle / 2 / M_PI) * 2 * M_PI;
  return angle;
}

void UnivRec::SetSPCoordinates(int idet)
{
  RecStation_t* station = &fstation[idet];

  //----Shower plane coordinates

  double xg0 = station->xg - recinfo.xgcore;
  double yg0 = station->yg - recinfo.ygcore;
  double zg0 = station->zg - recinfo.zgcore;

  double costheta = cos(recinfo.theta);
  double sintheta = sin(recinfo.theta);
  double cosazi = cos(recinfo.azi);
  double sinazi = sin(recinfo.azi);

  double xr = xg0 * cosazi + yg0 * sinazi;
  double yr = -xg0 * sinazi + yg0 * cosazi;
  double zr = zg0;

  station->r = sqrt((-xr * costheta + zr * sintheta) * (-xr * costheta + zr * sintheta)  + yr * yr);
  station->psi = NormalizeAngle(atan2(yr, xr * costheta - zr * sintheta));

  //----(DX.DL,T0_CF,T0_PF)
  station->T0_PF = fatm->GetTimePF_h(station->r, station->psi, recinfo.zgcore, recinfo.theta, station->zg);
  for (int icomp = 0; icomp < 4; icomp++) {

    double Xsignal = recinfo.Xmax + UnivParamNS::XmaxShift[icomp];
    double Xtime  = recinfo.Xmax + UnivParamTimeNS::XmaxShift_time[icomp];

    station->DX[icomp] = fatm->Get_DX_DL(station->r,  station->psi,  Xsignal   , recinfo.theta,  station->zg,
                                         UnivParamNS::UseDL[icomp], UnivParamNS::UseDiffusive[icomp]);
    station->DL[icomp] = fatm->Get_DX_DL(station->r,  station->psi,  Xtime     , recinfo.theta,  station->zg,
                                         UnivParamTimeNS::UseDL_time[icomp], UnivParamTimeNS::UseDiffusive_time[icomp]);


    double D_TO = ftime[station->type]->GetD_TO(recinfo.Xmax, recinfo.theta, recinfo.logE, icomp);
    double h1 = fatm->SlantDepthToHeight(recinfo.Xmax, recinfo.theta) + D_TO * 1.e5 * cos(recinfo.theta);
    station->T0_CF[icomp] = fatm->GetTimeCF_h(station->r,  station->psi,  h1,  recinfo.theta, station->zg) ;
  }
  //--- Density at height of station
  station->density = (gRecMC ? fatm->GetDensity(station->zg) : recinfo.AugerDensity);


}

void UnivRec::SetAllSPCoordinates()
{
  for (int idet = 0; idet < nRecStations; idet++)
    SetSPCoordinates(idet);
}

//---------------------------------------------------
//---------------------------------------------------
//  Rec estimate (Only used if Auger Rec standard reconstruction parameters are not used as seed )
//---------------------------------------------------
//---------------------------------------------------

bool UnivRec::GetRecEstimate(int itype)
{
  GetRecSeed(itype);

  //Set Rec core as the barycenter
  recinfo.xgcore = 0.;
  recinfo.ygcore = 0.;

  float vemsqrt = 0.;
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* station = &fstation[idet];
    if (station->type != itype) continue;

    recinfo.xgcore += (station->xg * sqrt(station->signalsize));
    recinfo.ygcore += (station->yg * sqrt(station->signalsize));
    vemsqrt += sqrt(station->signalsize);
  }

  if (vemsqrt == 0.) return false;
  recinfo.xgcore /= vemsqrt;
  recinfo.ygcore /= vemsqrt;


  //Set Rec direction using Plane fit
  if (!PlaneFit(false, false, true, itype)) return false;

  float theta_prev = recinfo.theta;
  int nIter = 0;
  while (true) {
    SetAllSPCoordinates();
    if (!PlaneFit(true, true, true, itype)) {
      PlaneFit(false, false, true, itype);
      SetAllSPCoordinates();
      return true;
    }
    if (fabs(theta_prev - recinfo.theta) < 0.1 * M_PI / 180.) break;
    theta_prev = recinfo.theta;
    nIter++;
    if (nIter > 10) break;
  }
  return true;
}


bool UnivRec::PlaneFit(bool CurvCorr, bool AltitudeCorr, bool RecSeed, int itype)
{
  std::vector<double> fX, fY, fT, fW;
  float speedC = AtmosphereNS::c_cmperns;

  //----
  int nSel = 0;
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* station = &fstation[idet];
    if (station->type != itype) continue;
    if (RecSeed && !(idet == recinfo.seed[0] || idet == recinfo.seed[1] || idet == recinfo.seed[2])) continue;
    if (station->signalsize > vemSignalCut) nSel++;
  }
  if (nSel < 3) return false;
  //-----
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* station = &fstation[idet];
    if (station->type != itype) continue;
    if (station->signalsize <= vemSignalCut) continue;
    if (RecSeed && !(idet == recinfo.seed[0] || idet == recinfo.seed[1] || idet == recinfo.seed[2]))
      continue;

    float tstation = station->GPSnano + station->starttime;
    float t50 = station->tquantile[0][2] - station->starttime;
    float theta = (recinfo.theta == 0. ? 38.*M_PI / 180. : recinfo.theta); //If iterating
    float nMuons = station->signalsize / (cos(theta) + 2.*tankh * sin(theta) / M_PI / tankr);
    if (nMuons < 2) continue;
    float TimeVariance2 = 134 + 2.4 * pow((t50 + 10.) / (nMuons + 1), 2) * nMuons / (nMuons + 2);

    float tcorr = 0.;
    if (AltitudeCorr) //Above ground-> Positive corr
      tcorr = (station->zg - recinfo.zgcore) * cos(recinfo.theta) / speedC;
    if (CurvCorr) {   //Negative correction
      float R = (1. - cos(recinfo.theta)) * 34 + 8.; //km  ( theta=0. 8 km, theta=60 25 km)
      R *= 1.e5; //cm
      float tcurv = station->r * station->r / speedC / 2 / R;
      tcorr -= tcurv;
    }
    tstation += tcorr;

    fW.push_back(1. / TimeVariance2);
    fX.push_back(station->xg);
    fY.push_back(station->yg);
    fT.push_back(tstation);
  }
  if (fX.size() < 3) return false;
  if (!PlaneFit(fX, fY, fT, fW)) return false;
  return true;
}

bool UnivRec::PlaneFit(std::vector<double> fX, std::vector<double>fY, std::vector<double>fT, std::vector<double>fW)
{
  if (fX.size() == 0.) return false;

  double speedC = AtmosphereNS::c_cmperns;

  //---Use the weighted time center as the origin of coordinates, so T0=0
  double fX0 = 0., fY0 = 0., fT0 = 0.;
  double WSum = 0.;
  for (unsigned int i = 0; i < fX.size(); i++) {
    WSum += (fW[i]);
    fX0 += (fX[i] * fW[i]);
    fY0 += (fY[i] * fW[i]);
    fT0 += (fT[i] * fW[i]);
  }
  //printf("\n");
  for (unsigned int i = 0; i < fX.size(); i++) {
    fX[i] = (fX[i] - fX0 / WSum);
    fY[i] = (fY[i] - fY0 / WSum);
    fT[i] = (fT[i] - fT0 / WSum);
  }
  //------------------

  float Det, Det1, Det2;
  float b[2][2], c[2];
  float l, m, lerr, merr, lmerr;

  //----Calculates the matrix of the linear system--------------------------------
  b[0][0] = 0.;
  b[0][1] = 0.;
  b[1][0] = 0.;
  b[1][1] = 0.;
  c[0] = 0.;
  c[1] = 0.;

  for (unsigned int i = 0; i < fX.size(); i++) {
    b[0][0] += fX[i] * fX[i]    * fW[i];
    b[0][1] += fX[i] * fY[i]    * fW[i];
    b[1][1] += fY[i] * fY[i]    * fW[i];

    c[0] -= (fT[i] * speedC * fX[i] * fW[i]);
    c[1] -= (fT[i] * speedC * fY[i] * fW[i]);
  }
  b[1][0] = b[0][1];

  //---Calculates the solution of the linear system using Cramer determinants-----
  Det = b[1][1] * b[0][0] - b[0][1] * b[1][0];
  Det1 = c[0] * b[1][1] - b[1][0] * c[1];
  Det2 = c[1] * b[0][0] - b[1][0] * c[0];
  l = Det1 / Det;
  m = Det2 / Det;

  float lm2;
  lm2 = l * l + m * m;

  if (lm2 > 1. || l == 0.) return false;

  //---Calculates the inverse of the matrix, i.e. the variance matrix-------------
  float deter = 0;
  deter = b[0][0] * b[1][1] * WSum - b[0][1] * b[1][0] * WSum;

  lerr = (speedC * speedC) * WSum * b[1][1] / deter;
  merr = (speedC * speedC) * WSum * b[0][0] / deter;
  lmerr = (speedC * speedC) * WSum * b[0][1] / deter;

  //---Calculates the angles---------------------------------------------------------
  float theta = asin(sqrt(l * l + m * m));
  float azi = atan2(m, l);
  float theta_err, azi_err;

  theta_err = cos(azi) * cos(azi) * lerr + sin(azi) * sin(azi) * merr + 2 * cos(azi) * sin(azi) * lmerr;
  if (theta_err > 0) theta_err = sqrt(theta_err) / cos(theta);
  azi_err = cos(azi) * cos(azi) * merr + sin(azi) * sin(azi) * lerr - 2 * cos(azi) * sin(azi) * lmerr ;
  if (azi_err > 0) azi_err = sqrt(azi_err) / sin(theta);

  //---Set reconstruction
  recinfo.theta = theta;
  recinfo.azi = NormalizeAngle(azi);

  return true;
};

void UnivRec::GetRecSeed(int itype)
{
  // The triangle with the larges sum of the signal

  float minDistance = (!gRecInfill ? 1200.e2 : 700.e2);
  float maxDistance = (!gRecInfill ? 1700.e2 : 800.e2);
  float SizeMax = 0.;

  for (int i = 0; i < nRecStations; i++) {
    if (fstation[i].type != itype) continue;
    for (int j = (i + 1); j < nRecStations; j++) {
      if (fstation[j].type != itype) continue;

      float dij = DistPoints(fstation[i].xg,  fstation[i].yg,  fstation[j].xg,  fstation[j].yg);
      if (dij > maxDistance || dij < minDistance) continue;

      for (int k = 0; k < nRecStations; k++) {
        if (k == i || k == j) continue;
        if (fstation[k].type != itype) continue;

        float dik = DistPoints(fstation[i].xg,  fstation[i].yg,  fstation[k].xg,  fstation[k].yg);
        if (dik > maxDistance || dik < minDistance) continue;
        float djk = DistPoints(fstation[j].xg,  fstation[j].yg,  fstation[k].xg,  fstation[k].yg);
        if (djk > maxDistance || djk < minDistance) continue;

        float Size = fstation[i].signalsize +  fstation[j].signalsize +  fstation[k].signalsize;
        if (Size > SizeMax) {
          SizeMax = Size;
          recinfo.seed[0] = i;
          recinfo.seed[1] = j;
          recinfo.seed[2] = k;
        }
      }//k
    }// j
  }// i

}

//---------------------------------------------------
//---------------------------------------------------
//  Utilities
//---------------------------------------------------
//---------------------------------------------------


void UnivRec::PrintRecInfo(bool PrintResiduals)
{
  // Print debug info
  printf("\n");
  printf("---------------------------------------\n");
  printf(" Fit results   ( FitStage = %d ) \n", recinfo.FitStage);
  printf("---------------------------------------\n");

  printf("Id=%d \n", recinfo.id);
  printf("theta=%5.2lf (%5.2lf) azi=%5.2lf (%5.2lf)  (x,y,z)=(%5.2lf,%5.2lf,%5.2lf)  (%5.2lf,%5.2lf,%5.2lf)  sigma_x,sigma_y=(%5.2lf,%5.2lf) \n",
         recinfo.theta * 180 / M_PI, (gRecMC ? MCinfo.theta : augerinfo.theta_FD) * 180. / M_PI  ,
         recinfo.azi * 180 / M_PI, (gRecMC ? MCinfo.azi : augerinfo.azi_FD) * 180. / M_PI  ,
         recinfo.xgcore / 1.e2, recinfo.ygcore / 1.e2, recinfo.zgcore / 1.e2,
         (gRecMC ? MCinfo.xgcore : augerinfo.xgcore_SD) / 1.e2 , (gRecMC ? MCinfo.ygcore : augerinfo.ygcore_SD) / 1.e2 , (gRecMC ? MCinfo.zgcore : augerinfo.zgcore_SD) / 1.e2 ,
         recinfo.xgcore_err / 1.e2, recinfo.ygcore_err / 1.e2);

  printf("logE=%5.2lf +-%5.2lf  (%5.2lf), Nmu=%5.2lf+-%5.2lf (%5.2lf) Xmax=%5.2lf +- %5.2lf   (%5.2lf)   \n",
         recinfo.logE, recinfo.logE_err, MCinfo.logE,
         recinfo.Nmu , recinfo.Nmu_err , MCinfo.Nmu,
         recinfo.Xmax, recinfo.Xmax_err, MCinfo.Xmax);
  printf(" OffsetM_Mu=%5.2lf %5.2lf \n", recinfo.OffsetM_Mu, recinfo.OffsetM_Mu_err);
  printf(" T0=%5.2lf+-%5.2lf (%5.2lf) \n" , recinfo.T0, recinfo.T0_err, MCinfo.T0);

  printf("\n");
  printf(" logE_ref=%5.2lf +-%5.2lf \n", recinfo.logE_ref, recinfo.logE_ref_err);
  printf("\n");
  printf(" logE_SD=%5.2lf theta_SD=%5.2lf azi_SD=%5.2lf   \n", augerinfo.logE_SD, augerinfo.theta_SD * 180. / M_PI, augerinfo.azi_SD * 180. / M_PI);
  printf(" logE_FD=%5.2lf Xmax_FD=%5.2lf  theta_FD=%5.2lf  azi_FD=%5.2lf \n", augerinfo.logE_FD, augerinfo.Xmax_FD, augerinfo.theta_FD * 180. / M_PI, augerinfo.azi_FD * 180. / M_PI);

  printf("\n");

  printf(" lnP=%5.2lf/%5.2lf/%5.2lf (%5.2lf) \n", recinfo.lnP[0], recinfo.lnP[1], recinfo.lnP[2],
         recinfo.lnP[0] + recinfo.lnP[1] + recinfo.lnP[2]);
  printf(" FitStatus=%d/%d  nFCNcalls=%d    Min EigenVal=%5.2lf   EDM=%5.2lf  IsTimeRec=%d \n", recinfo.FitStatus[0], recinfo.FitStatus[1], recinfo.nFCNcalls,
         recinfo.MinEigenVal, recinfo.edm, IsTimeRec);
  if (recinfo.FitStage > 0) printf("   TimeRejected=%d (%lf) | TimeMaxDev=%d (%lf) \n", int(recinfo.TimeRejected[1]), recinfo.TimeRejected[0], int(recinfo.TimeMaxdev[1]), recinfo.TimeMaxdev[0]);
  else                       printf("\n");
  printf(" AllTimesOK= %d \n", recinfo.AllTimesOK);
  printf(" gDetHot=%d r=%5.2lf  [m]\n", recinfo.gDetHot, recinfo.rHot / 1.e2);
  if (!gRecMC) printf(" Density %lf [g]/[cm3] \n", recinfo.AugerDensity);
  printf(" month=%d \n", recinfo.iatm);

  printf("nSignalStations=%d  nTimeStations=%d \n", recinfo.nSignalStations, recinfo.nTimeStations);

  // Print residuals
  if (PrintResiduals) {
    for (int iloop = 0; iloop < 2; iloop++) {
      if (iloop == 0) printf("Signal likelihood \n");
      else             printf("Time likelihood \n");
      printf("--------------------------\n");

      for (int idet = 0; idet < nRecStations; idet++) {
        RecStation* st = &fstation[idet];
        if (iloop == 0  && st->mask_signal) continue;
        if (iloop == 1  && recinfo.FitStage < 1) continue;
        if (iloop == 1  && st->mask_quantiles && st->mask_starttime) continue;
        printf("idet=%d (%d) %d | %4.0lf %4.0lf %4.0lf %4.0lf %d %d ", idet, (st->id), st->type, st->r / 1.e2, st->psi * 180. / M_PI,
               st->DX[0], st->DX[1], st->isSaturated, st->mask_signal);

        if (iloop == 0)
          printf(" | Rec=%5.0lf Pred=%5.0lf RMS=%5.0lf lnP=%5.2lf \n", st->signalsize, st->signalsize_pred, st->signalsize_rms, st->lnP_signal);
        else {
          printf(" %4.0lf %4.2lf| t0_cf=%3.0lf | %d %d St= %4.0lf (%4.0lf %4.0lf)| %d %d%d%d  T1= %4.0lf (%4.0lf %4.0lf) T10=%4.0lf (%4.0lf %4.0lf) T50=%4.0lf (%4.0lf %4.0lf)  | lnP=%3.1lf %3.1lf  \n",
                 st->signalsize, st->AreaOverPeak, st->T0_CF[0],
                 st->mask_starttime, st->UseGausPDF_start,
                 st->start_rec + st->start_corr, st->start_pred + st->start_corr, st->start_rms,
                 st->mask_quantiles, st->UseGausPDF_quantiles[0], st->UseGausPDF_quantiles[1], st->UseGausPDF_quantiles[2],
                 st->tq_rec[0] + st->tq_corr[0], st->tq_pred[0] + st->tq_corr[0], st->tq_rms[0],
                 st->tq_rec[1] + st->tq_corr[1], st->tq_pred[1] + st->tq_corr[1], st->tq_rms[1],
                 st->tq_rec[2] + st->tq_corr[2], st->tq_pred[2] + st->tq_corr[2], st->tq_rms[2],
                 st->lnP_start, st->lnP_quantiles);
        }
      }
    }// Signal/Time
  }
}

#ifndef ExternalUse
void UnivRec::WriteTextFile(FILE* fp)
{
  // Print short summary to a file

  TVector3 theDir0, theDir;
  theDir0.SetMagThetaPhi(1., (gRecMC ? MCinfo.theta : augerinfo.theta_FD), (gRecMC ? MCinfo.azi : augerinfo.azi_FD));
  theDir.SetMagThetaPhi(1., recinfo.theta, recinfo.azi);

  bool okRec = IsGoodRec();

  double Chi2[2] = {0., 0.};
  int nChi2[2] = {0, 0};
  if (okRec) GetChi2(Chi2, nChi2);


  unsigned int GPSnano = 0;
  utl::TimeStamp time_stamp(GPSsec, GPSnano);
  unsigned int month = (!gRecMC ? time_stamp.GetMonth() : recinfo.iatm);
  unsigned int year = (!gRecMC ? time_stamp.GetYear() : UNDEF);


  fprintf(fp, "%d %d %d %d   %d %d   %lf %lf %lf %lf   %lf %lf  %lf %lf   %lf %lf %lf %lf  %lf %d %d    %lf %lf %d   %lf %lf %lf   %lf %lf   %lf %lf   %lf   %lf %lf %lf  %lf %d %lf %d  %lf %lf   %d %d   %d %d %lf %lf %d %d %d %lf %lf  \n",
          recinfo.id, gRecType, gRecAoP, gRecSDE_Upgrade, MCinfo.im, MCinfo.ip,
          recinfo.logE, recinfo.logE_err, MCinfo.logE, augerinfo.logE_FD,
          recinfo.xgcore, recinfo.ygcore,
          MCinfo.xgcore, MCinfo.ygcore,
          recinfo.Xmax, recinfo.Xmax_err, MCinfo.Xmax, augerinfo.Xmax_FD,
          recinfo.OffsetM_Mu, okRec, month ,
          augerinfo.logE_SD, recinfo.logE_ref, year,
          recinfo.Nmu, recinfo.Nmu_err, MCinfo.Nmu,
          recinfo.theta, MCinfo.theta,
          recinfo.azi, MCinfo.azi,
          theDir0.Angle(theDir),
          recinfo.T0, recinfo.T0_err, MCinfo.T0,
          recinfo.lnP[0], recinfo.FitStatus[0],
          recinfo.lnP[1], recinfo.FitStatus[1],
          Chi2[0], Chi2[1],
          nChi2[0], nChi2[1],
          recinfo.IsSaturated[0], recinfo.IsSaturated[1],
          recinfo.rHot, augerinfo.Xmax_FD_err,
          recinfo.nTimeStations, recinfo.FitStage, recinfo.AllTimesOK,
          recinfo.TimeRejected[0], recinfo.TimeMaxdev[0]);
}
#endif
bool UnivRec::RWRecinfo(FILE* fp, bool WriteStations, bool readFlag)
{
  if (fp == NULL) return false;
  if (readFlag && feof(fp)) return false;

  if (!readFlag) {
    fwrite(&gRecType, 1, sizeof(int), fp);
    fwrite(&gRecMC, 1, sizeof(bool), fp);
    fwrite(&gRecInfill, 1, sizeof(bool), fp);
    fwrite(&gRecSDE_Upgrade, 1, sizeof(bool), fp);
    fwrite(&gRecAoP, 1, sizeof(int), fp);
    int DataSet_Sys = (gRecMC ? UNDEF : gRecDataSet + gRecSys * 10);
    fwrite(&DataSet_Sys, 1, sizeof(int), fp);
    fwrite(&gRecMixture, 1, sizeof(int), fp);
    fwrite(&gCalibOpt, 1, sizeof(int), fp);
    for (int itype = 0; itype < 3; itype++)
      fwrite(&gDetectorTypeMask[itype], 1, sizeof(bool), fp);

    fwrite(&recinfo, 1, sizeof(RecInfo_t), fp);
    fwrite(&MCinfo, 1, sizeof(MCInfo_t), fp);
    fwrite(&augerinfo, 1, sizeof(AugerRecInfo_t), fp);
    if (!gRecMC) fwrite(&GPSsec, 1, sizeof(unsigned int), fp);

    int nDet = 0;
    if (!WriteStations)
      fwrite(&nDet, 1, sizeof(int), fp);
    else {
      for (int idet = 0; idet < nRecStations; idet++)
        if (fstation[idet].signalsize > 0.) nDet++;
      fwrite(&nDet, 1, sizeof(int), fp);
      for (int idet = 0; idet < nRecStations; idet++)
        if (fstation[idet].signalsize > 0.)
          fwrite(&fstation[idet], 1, sizeof(RecStation_t), fp);
    }
  } else {
    [[maybe_unused]]int nread = 0; // unused. LN.
    nread += fread(&gRecType, 1, sizeof(int), fp);
    IsGoldenRec = (gRecType == 0 || gRecType == 7);
    if (feof(fp)) return false;
    nread += fread(&gRecMC, 1, sizeof(bool), fp);
    nread += fread(&gRecInfill, 1, sizeof(bool), fp);
    nread += fread(&gRecSDE_Upgrade, 1, sizeof(bool), fp);
    nread += fread(&gRecAoP, 1, sizeof(int), fp);
    int DataSet_Sys;
    nread += fread(&DataSet_Sys, 1, sizeof(int), fp);
    gRecDataSet = (gRecMC ? UNDEF : DataSet_Sys % 10);
    gRecSys = (gRecMC ? 0 : DataSet_Sys / 10);
    if (!gRecMC && (gRecDataSet < 0 || gRecDataSet > 3 || gRecSys < 0 || gRecSys > UnivCalibConstantsNS::nSys))
      return false;

    nread += fread(&gRecMixture, 1, sizeof(int), fp);
    nread += fread(&gCalibOpt, 1, sizeof(int), fp);

    for (int itype = 0; itype < 3; itype++)
      nread += fread(&gDetectorTypeMask[itype], 1, sizeof(bool), fp);

    nread += fread(&recinfo, sizeof(RecInfo_t), 1, fp);
    nread += fread(&MCinfo, sizeof(MCInfo_t), 1, fp);
    nread += fread(&augerinfo, sizeof(AugerRecInfo_t), 1, fp);
    if (!gRecMC)  nread += fread(&GPSsec, sizeof(unsigned int), 1, fp);

    nRecStations = 0;
    nread += fread(&nRecStations, sizeof(int), 1, fp);
    if (feof(fp)) return false;
    fstation.clear();

    if (nRecStations > 0) {
      for (int idet = 0; idet < nRecStations; idet++) {
        RecStation_t station;
        nread += fread(&station, 1, sizeof(RecStation_t), fp);
        fstation.push_back(station);
      }
    }

    //---
    IsTimeRec = (gRecType == 1 || gRecType == 4 || gRecType == 5 || gRecType == 6 || gRecType == 7 || gRecType == 8 ? true : false);

    if (fcalib != NULL) delete fcalib;
    fcalib = new UnivCalibConstantsNS::UnivCalibConstants(gCalibOpt, gRecSys , gRecMixture);

    fatm = new AtmosphereNS::Atmosphere();
    fatm->SetCurrentAtmosphere(recinfo.iatm);

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

    //----
  }
  return true;
}

void UnivRec::GetChi2(double* Chi2, int* nChi2)
{
  for (int i = 0; i < 2; i++) {
    Chi2[i] = 0.;
    nChi2[i] = 0;
  }
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* st = &fstation[idet];
    if (!st->mask_signal && !st->isSaturated) {
      Chi2[0] += pow((st->signalsize_pred - st->signalsize) / st->signalsize_rms, 2.);
      nChi2[0] += 1;
    }

    if (!st->mask_starttime && st->start_pred != 0. && st->start_rec != 0. &&  st->start_rms != 0.) {
      Chi2[1] += pow((st->start_pred - st->start_rec) / st->start_rms, 2.);
      nChi2[1] += 1;
    }


    if (!st->mask_quantiles) {
      for (int iq = 0; iq < nQ; iq++) {
        if (st->quantile[iq] == UNDEF) continue;
        if (st->tq_rms[iq] == 0. || st->tq_pred[iq] == 0. || st->tq_rec[iq] == 0.) continue;
        Chi2[1] += pow((st->tq_pred[iq] - st->tq_rec[iq]) / st->tq_rms[iq], 2.);
        nChi2[1] += 1;
      }
    }
  }
}


void UnivRec::GetResiduals(std::vector<double>& mean, std::vector<double>& rms, std::vector<double>& rec, int optY, std::vector<double>& fx, int optX)
{
  GetResiduals(mean, rms, rec, optY, fx, optX, false);
}

void UnivRec::GetResiduals(std::vector<double>& mean, std::vector<double>& rms, std::vector<double>& rec, int optY, std::vector<double>& fx, int optX , bool OnlyNotSaturated)
{
  // opt=0,1,2,3,4  Signal/Start time/T1/T10/T50
  if (optY < 0 || optY > 4) return;
  if (optX < 0 || optX > 2) return;
  mean.clear();
  rms.clear();
  rec.clear();
  fx.clear();

  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* st = &fstation[idet];


    if (OnlyNotSaturated && st->isSaturated) continue;
    if (optY == 0 && (st->mask_signal || st->isSaturated)) continue;
    if (optY == 1 && st->mask_starttime) continue;
    if (optY > 1  && (st->mask_quantiles || st->quantile[optY - 2] == UNDEF)) continue;

    double x_i;
    if (optX == 0) x_i = st->r;
    else if (optX == 1) x_i = st->AreaOverPeak;


    if (optY == 0) {
      mean.push_back(st->signalsize_pred);
      rms.push_back(st->signalsize_rms);
      rec.push_back(st->signalsize);
      fx.push_back(x_i);
    } else if (optY == 1) {
      mean.push_back(st->start_pred);
      rms.push_back(st->start_rms);
      rec.push_back(st->start_rec);
      fx.push_back(x_i);
    } else {
      mean.push_back(st->tq_pred[optY - 2]);
      rms.push_back(st->tq_rms[optY - 2]);
      rec.push_back(st->tq_rec[optY - 2]);
      fx.push_back(x_i);
    }
  }
}

void UnivRec::GetMeanRMS_AoP(double& mean, double& rms)
{
  mean = 0.;
  rms = 0.;
  double nStations = 0.;
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* st = &fstation[idet];
    if (st->mask_quantiles) continue;

    mean += (st->AreaOverPeak);
    rms += pow(st->AreaOverPeak, 2.);
    nStations += 1.;
  }
  if (nStations < 2) {
    mean = 0.;
    rms = 0.;
    return;
  }
  mean /= nStations;
  rms = sqrt(nStations / (nStations - 1) * (rms / nStations - mean * mean));
}

bool UnivRec::ScanRange(int iparX, int iparY, std::vector<double> par, std::vector<int> parStatus,  double* rangeX, double* rangeY, double* d, double RotAngle, double errDef, double& edm)
{
  double InitStep = 2.0;
  return ScanRange(iparX, iparY, par, parStatus, rangeX, rangeY, d, RotAngle, errDef, InitStep, edm);
}

bool UnivRec::ScanRange(int iparX, int iparY, std::vector<double> par, std::vector<int> parStatus,  double* rangeX, double* rangeY, double* d, double RotAngle, double errDef, double InitStep, double& edm)
{
  // Note: iside refers to positive increment after the coordinates are rotated.

  std::vector<bool> isFixed;
  for (unsigned int ipar = 0; ipar < par.size(); ipar++)
    isFixed.push_back(parStatus[ipar] < 2 ? false : true);

  edm = 0.;
  bool okX = (iparX >= 0 && iparX < int(par.size()));
  bool okY = (iparY >= 0 && iparY < int(par.size()));
  if (okX) okX = !isFixed[iparX];
  if (okY) okY = !isFixed[iparY];
  if (!okX && !okY) return false;

  if (okX && okY && iparX == iparY)  okY = false;
  if (!okX) RotAngle = M_PI / 2.;
  if (!okY) RotAngle = 0.;

  const int niter_max = (InitStep < 1. ? 500 : 100);

  double lnP_max = GetLikelihood(par, parStatus);
  double val_max[2] = { (okX ? par[iparX] : UNDEF) , (okY ? par[iparY] : UNDEF) };
  bool ok_range = true;

  double lowX = UNDEF, highX = UNDEF;
  if (okX) {
    lowX = low[iparX] / unit[iparX];
    highX = high[iparX] / unit[iparX];
    if (iparX == 4) {
      lowX = XmaxLow / unit[iparX];
      highX = XmaxHigh / unit[iparX];
    }
    if (iparX == 5) {
      lowX = (val_max[0] + lowX);
      highX = (val_max[0] + highX);
    }
  }
  double lowY = UNDEF, highY = UNDEF;
  if (okY) {
    lowY = low[iparY] / unit[iparY];
    highY = high[iparY] / unit[iparY];
    if (iparY == 4) {
      lowY = XmaxLow / unit[iparY];
      highY = XmaxHigh / unit[iparY];
    }
    if (iparY == 5) {
      lowY = (val_max[1] + lowY);
      highY = (val_max[1] + highY);
    }
  }

  for (int iside = 0; iside < 2; iside++) {
    double sign = (iside == 0 ? -1. : 1.);
    double t_iter = 0., lnP_iter = 0., step = InitStep;
    int niter = 0;

    while (true) {
      t_iter += (sign * step);
      if (okX)   par[iparX] = val_max[0] + cos(RotAngle) * t_iter;
      if (okY)   par[iparY] = val_max[1] + sin(RotAngle) * t_iter;
      lnP_iter = GetLikelihood(par, parStatus);

      if (UnivRecNS::debug > 1) {
        printf(" ipar=(%d,%d) iside=%d iter=%d t=%5.4lf step=%5.4lf lnP=%lf (%lf %lf)  ", iparX, iparY, iside, niter, t_iter, step, lnP_iter, lnP_max, errDef);
        if (okX) printf(" parX=%lf (%lf) ", par[iparX]*unit[iparX], val_max[0]*unit[iparX]);
        if (okY) printf(" parY=%lf (%lf) ", par[iparY]*unit[iparY], val_max[0]*unit[iparY]);
        printf("\n");
      }

      //----------
      if (lnP_iter > lnP_max) edm = lnP_iter - lnP_max;

      //----------
      if ((okX && (par[iparX] < lowX || par[iparX] > highX)) ||
          (okY && (par[iparY] < lowY || par[iparY] > highY)) ||  niter > niter_max) {
        ok_range = false;
        break;
      }

      //----------
      if (lnP_iter < (lnP_max - errDef)) {
        if ((lnP_max - errDef) - lnP_iter < 0.01) break;
        else {
          t_iter -= (sign * step);
          step /= 2.;
        }
      }
      //----------
      niter++;
    }

    if (!ok_range) break;

    //--Interpolate

    double t_i[2] = { t_iter - sign * step , t_iter };
    double lnP_i[2] = { 0., lnP_iter};

    if (okX)   par[iparX] = val_max[0] + cos(RotAngle) * t_i[0];
    if (okY)   par[iparY] = val_max[1] + sin(RotAngle) * t_i[0];

    lnP_i[0] = GetLikelihood(par, parStatus);

    double t = (t_i[1] - t_i[0]) * (lnP_max - errDef - lnP_i[0]) / (lnP_i[1] - lnP_i[0]) +  t_i[0];
    if (okX)   rangeX[iside] = val_max[0] + cos(RotAngle) * t;
    if (okY)   rangeY[iside] = val_max[1] + sin(RotAngle) * t;
    d[iside] = t;

    //--
    if (UnivRecNS::debug > 1) {

      printf(" ipar=(%d,%d) iside=%d niter=%d lnP=%lf %lf (%lf)  RotAngle=%4.2lf ", iparX, iparY, iside, niter, lnP_i[0], lnP_i[1], lnP_max, RotAngle * 180. / M_PI);
      if (okX) printf(" parX=%lf (%lf) ", rangeX[iside]*unit[iparX], val_max[0]*unit[iparX]);
      if (okY) printf(" parY=%lf (%lf) ", rangeY[iside]*unit[iparY], val_max[1]*unit[iparY]);

      if (okX) par[iparX] = rangeX[iside];
      if (okY) par[iparY] = rangeY[iside];
      printf(" | d=%5.2lf  %lf \n", d[iside], lnP_max - GetLikelihood(par, parStatus));
    }
  }

  // Set ranges to low/high if failed
  if (!ok_range) {
    if (okX) {
      rangeX[0] = lowX;
      rangeX[1] = highX;
    }
    if (okY) {
      rangeY[0] = lowY;
      rangeY[1] = highY;
    }
  }

  // Restore
  if (okX) par[iparX] = val_max[0];
  if (okY) par[iparY] = val_max[1];
  GetLikelihood(par, parStatus);

  bool ok_minima = (edm < edm_tolerance);

  return ok_range && ok_minima;
}


void UnivRec::ScanLikelihood(int iparX, int iparY, int N, std::vector<std::pair<double, double> >& parXY,  std::vector<double>& lnP,
                             std::pair<double, double>& parXY_min, double& lnP_min)
{
  parXY.clear();
  lnP.clear();

  std::vector<double> par, epar;
  std::vector<int> parStatus;
  std::vector<bool> isFixed, hasLimits;
  InitFCNParameters(par, epar, parStatus, hasLimits);
  for (unsigned int ipar = 0; ipar < par.size(); ipar++)
    isFixed.push_back(parStatus[ipar] < 2 ? false : true);

  bool okX = (iparX >= 0 && iparX < int(par.size()));
  if (okX) okX = !isFixed[iparX];
  bool okY = (iparY >= 0 && iparY < int(par.size()));
  if (okY) okY = !isFixed[iparY];
  if (!okX && !okY) return;

  double range[2][2] = { {1.e10, -1.e10}, {1.e10, -1.e10} }; // [ix/iy][min/max]
  double dum[2], edm;
  double step[2] = {0., 0.};
  unsigned int NX = 1, NY = 1;

  if (okX && okY) {
    for (int irot = 0; irot < 8; irot++) {
      double RotAngle = M_PI / 8.*double(irot);
      double rangeX[2], rangeY[2];
      ScanRange(iparX, iparY, par, parStatus, rangeX, rangeY, dum, RotAngle, 6., edm);
      for (int iside = 0; iside < 2; iside++) {
        if (okX && rangeX[iside] < range[0][0]) range[0][0] = rangeX[iside];
        if (okX && rangeX[iside] > range[0][1]) range[0][1] = rangeX[iside];
        if (okY && rangeY[iside] < range[1][0]) range[1][0] = rangeY[iside];
        if (okY && rangeY[iside] > range[1][1]) range[1][1] = rangeY[iside];
      }
    }
  } else ScanRange(iparX, iparY, par, parStatus, range[0], range[1], dum, (okX ? 0. : M_PI / 2.) , 6., edm);


  if (UnivRecNS::debug > 1)
    printf("rangeX=(%5.2lf %5.2lf )   rangeY=(%5.2lf,%5.2lf ) \n", range[0][0]*unit[iparX], range[0][1]*unit[iparX], range[1][0]*unit[iparY], range[1][1]*unit[iparY]);

  if (okX) {
    step[0] = (range[0][1] - range[0][0]) / double(N);
    NX = N;
  }
  if (okY) {
    step[1] = (range[1][1] - range[1][0]) / double(N);
    NY = N;
  }

  for (unsigned int ix = 0; ix < NX; ix++)
    for (unsigned int iy = 0; iy < NY; iy++) {
      std::vector<double> par_i;
      for (unsigned int ipar = 0; ipar < par.size(); ipar++) par_i.push_back(par[ipar]);
      std::pair<double, double> parXY_i;
      if (okX) {
        parXY_i.first = range[0][0] + double(ix) * step[0];
        par_i[iparX] = parXY_i.first;
        parXY_i.first *= unit[iparX];
      } else         parXY_i.first = 0.;

      if (okY) {
        parXY_i.second = range[1][0] + double(iy) * step[1];
        par_i[iparY] = parXY_i.second;
        parXY_i.second *= unit[iparY];
      } else         parXY_i.second = 0.;
      double lnP_i = GetLikelihood(par_i, parStatus);

      parXY.push_back(parXY_i);
      lnP.push_back(lnP_i);
    }

  lnP_min = GetLikelihood(par, parStatus); // Restore the internal values for the minimum and recalculate the likelihood

  if (okX)
    parXY_min.first = par[iparX] * unit[iparX];
  if (okY)
    parXY_min.second = par[iparY] * unit[iparY];
}


bool
UnivRec::IsGoodRec()
{
  return IsGoodRec(true);
}


bool
#ifndef ExternalUse
UnivRec::IsGoodRec(bool CutBadYears)
#else
UnivRec::IsGoodRec(bool /*CutBadYears*/)
#endif
{
#ifndef ExternalUse
  if (!gRecMC && CutBadYears && gRecDataSet == 1) {
    unsigned int GPSnano = 0;
    utl::TimeStamp time(GPSsec, GPSnano);
    if (time.GetYear() < 2008) return false;
  }
#endif
  //-----
  bool okRec;
  if (!IsTimeRec)
    okRec = (recinfo.FitStatus[0] == 3 && recinfo.FitStage == 0);
  else {
    bool okQuality = (recinfo.Xmax_err < 200. && (gRecType == 6 ? true : recinfo.TimeMaxdev[0] < 4.0));
    bool okMinimum = (recinfo.FitStatus[1] >= 2 && recinfo.MinEigenVal > -0.5);
    okRec = (recinfo.FitStatus[0] == 3 &&  recinfo.FitStage == 2 && recinfo.AllTimesOK && okMinimum && okQuality);
  }
  //---

  if (!okRec && UnivRecNS::debug > 2) {

    printf("sdid=%d  ok=%d   Status=%d/%d FitStage=%d EDM=%5.3lf MinEigenVal=%5.3lf AlltimesOK=%d theta=%5.2f (%5.2lf %5.2lf) logE=%5.2lf (%5.2lf %5.2lf) Xmax=%5.2lf+-%5.2lf (%5.2lf %5.2lf ) Nmu=%5.2lf  nStations=%d/%d\n",
           recinfo.id, okRec, recinfo.FitStatus[0], recinfo.FitStatus[1], recinfo.FitStage,
           recinfo.edm, recinfo.MinEigenVal,
           recinfo.AllTimesOK,
           recinfo.theta * 180. / M_PI, MCinfo.theta * 180. / M_PI, augerinfo.theta_FD * 180. / M_PI,
           recinfo.logE, MCinfo.logE, augerinfo.logE_FD,
           recinfo.Xmax, recinfo.Xmax_err, MCinfo.Xmax, augerinfo.Xmax_FD,
           recinfo.Nmu, recinfo.nSignalStations, recinfo.nTimeStations);

  }
  return okRec;
}



bool UnivRec::CheckEnergyCut(bool PrintOut)
{
  if (gRecMC) return true;

  double logE = (gRecDataSet == 1 ? augerinfo.logE_FD : augerinfo.logE_SD);
  double logE_Thresh = 18.60;
  if (!IsTimeRec || gRecInfill) logE_Thresh = 17.00;
  else if (gRecType == 1) logE_Thresh = 19.10;
  else if (gRecSys   > 0) logE_Thresh = 19.00;

  if (logE < logE_Thresh) {
    if (PrintOut) printf("%d Event not selected : logE=%5.2lf \n", recinfo.id, logE);
    return false;
  }
  return true;
}


//---------------------------------------------------
//---------------------------------------------------
// Time Likelihood
//---------------------------------------------------
//---------------------------------------------------
double UnivRec::GetTimeLikelihood()
{
  return GetTimeLikelihood(false);
}

double UnivRec::GetTimeLikelihood(bool UseAvgQuantiles)
{
  recinfo.TimeMaxdev[0] = -1.e10;
  recinfo.TimeMaxdev[1] = UNDEF;
  recinfo.AllTimesOK = true;
  recinfo.nTimeStations = 0;
  recinfo.lnP[0] = 0.;


  for (int idet = 0; idet < nRecStations; idet++) {

    RecStation_t* st = &fstation[idet];
    for (int iq = 0; iq < nQ; iq++) {
      st->tq_rec[iq] = UNDEF;
      st->tq_pred[iq] = UNDEF;
      st->tq_rms[iq] = UNDEF;
      st->tq_corr[iq] = 0.;
      st->lnP_quantiles = 0.;

      st->start_rec = UNDEF;
      st->start_pred = UNDEF;
      st->start_rms = UNDEF;
      st->start_corr = 0.;
      st->lnP_start = 0.;
    }
    int type = st->type;
    if (st->mask_starttime      &&  st->mask_quantiles) continue;
    recinfo.nTimeStations++;

    //------------------------------
    //
    //------------------------------

    if (!Check_DX_DL(st)) {
      st->lnP_quantiles = lnP_inf;
      st->lnP_start = lnP_inf;
      recinfo.AllTimesOK = false;
      recinfo.lnP[0] += (st->lnP_start + st->lnP_quantiles);
      continue;
    }

    double vemTot_p = 0., fcomp[4];
    for (int icomp = 0; icomp < 4; icomp++) {
      fcomp[icomp] = fsignal[type]->GetSignal(st->r, st->DX, recinfo.logE, recinfo.theta, st->psi,
                                              st->density, st->zg, recinfo.Nmu, icomp);
      vemTot_p += fcomp[icomp];
    }
    for (int icomp = 0; icomp < 4; icomp++) fcomp[icomp] /= vemTot_p;



    //---------------
    // Start time
    //---------------

    st->start_corr = ftime[type]->tStartCorrection(st->r, recinfo.logE, recinfo.theta, gRecSDE_Upgrade);
    st->start_rec = (st->GPSnano + st->starttime - st->start_corr) - (recinfo.T0 + st->T0_PF);

    double UnitPerMuon = 0;
    if (st->type == 0)
      UnitPerMuon = 1. / (cos(recinfo.theta) + 2.*tankh * sin(recinfo.theta) / M_PI / tankr);
    else if (st->type == 1)
      UnitPerMuon = 1. / cos(recinfo.theta);
    else if (st->type == 2)
      UnitPerMuon = 1.;
    double nmuDetector = fcomp[0] * st->signalsize / UnitPerMuon;
    if (nmuDetector <= 1.)
      nmuDetector = 1.;
    double P_extreme = ftime[type]->tFirstPDF(st->DL[0], st->r, recinfo.logE, st->psi, recinfo.theta, st->T0_CF[0], nmuDetector, st->start_rec
                       , false, st->start_pred, st->start_rms);
    double Accuracy = (st->type == 2 ? StartAccuracyAMIGA : StartAccuracy[gRecSDE_Upgrade]);
    st->start_rms = sqrt(pow(st->start_rms, 2.) + pow(Accuracy, 2.) + pow(InterGPSAccuracy[gRecSDE_Upgrade], 2.));

    if (!st->mask_starttime) {
      double dev = (st->start_rec - st->start_pred) / st->start_rms;
      if (st->UseGausPDF_start)
        st->lnP_start += (-pow(dev, 2.) / 2. + log(1. / sqrt(2.*M_PI) / st->start_rms));
      else if (P_extreme < exp(lnP_inf))
        st->lnP_start += lnP_inf;
      else
        st->lnP_start += log(P_extreme);

      if (fabs(dev) > recinfo.TimeMaxdev[0]) {
        recinfo.TimeMaxdev[0] = fabs(dev);
        recinfo.TimeMaxdev[1] = idet;
      }
      if (!st->UseGausPDF_start && (st->start_rec < st->T0_CF[0] || P_extreme < exp(lnP_inf)))
        recinfo.AllTimesOK = false;
    }

    //---------------
    // T1/T10/T50/T90
    //---------------
    for (int iq = 0; iq < nQ; iq++) {
      if (st->quantile[iq] == UNDEF)
        continue;
      double tq = (UseAvgQuantiles ? st->tquantile[1][iq] : st->tquantile[0][iq]);
      if (tq == UNDEF)
        continue;
      st->tq_rec[iq] = (st->GPSnano + tq) - (recinfo.T0 + st->T0_PF);
    }

    // Corrections are applied to the reconstructed quantiles since it is not possible to change the binomial probability accordingly
    // The predicted quantiles for large signals are used.

    double RiseTime = ftime[type]->GetTime(st->DL, st->r, recinfo.logE, st->psi, recinfo.theta, fcomp, st->T0_CF, 0.50) -
                      ftime[type]->GetTime(st->DL, st->r, recinfo.logE, st->psi, recinfo.theta, fcomp, st->T0_CF, 0.10);

    for (int iq = 0; iq < nQ; iq++) {
      if (st->quantile[iq] == UNDEF)
        continue;

      //Area Over Peak
      double corrAoP = ftime[type]->tQuantileCorrection_AoP(RiseTime,   iq, st->AreaOverPeak);
      st->tq_rec[iq] -= corrAoP;

      //Quantile corrections ( Pulse effects+Small signal shifts)
      double NanoSecPerVEM = RiseTime / st->signalsize / 0.4;
      double corr = ftime[type]->tQuantileCorrection(NanoSecPerVEM, iq);
      if (!UseAvgQuantiles)
        st->tq_rec[iq] -= corr;

      st->tq_corr[iq] = corrAoP + corr;
    }
    //------------------------------------
    //------------------------------------

    for (int iq = 0; iq < nQ; iq++) {
      if (st->quantile[iq] == UNDEF)        continue;

      double signalsize = (UseAvgQuantiles ? 1.e4 : st->signalsize);
      double P_binom = ftime[type]->tQuantilePDF(st->DL, st->r , recinfo.logE, st->psi, recinfo.theta, fcomp, st->T0_CF, signalsize,
                       st->tq_rec[iq], st->quantile[iq] , false, st->tq_pred[iq], st->tq_rms[iq]);
      st->tq_rms[iq] = sqrt(pow(st->tq_rms[iq], 2.) + pow(ParameterizationRMSLimit, 2.) + pow(InterGPSAccuracy[gRecSDE_Upgrade], 2.));

      if (!st->mask_quantiles) {
        double dev = (st->tq_rec[iq] - st->tq_pred[iq]) / st->tq_rms[iq];
        if (st->UseGausPDF_quantiles[iq])    st->lnP_quantiles += (-pow(dev, 2.) / 2. + log(1. / sqrt(2.*M_PI) / st->tq_rms[iq]));
        else if (P_binom < exp(lnP_inf))    st->lnP_quantiles += lnP_inf;
        else                                   st->lnP_quantiles += log(P_binom);

        if (fabs(dev) > recinfo.TimeMaxdev[0]) {
          recinfo.TimeMaxdev[0] = fabs(dev);
          recinfo.TimeMaxdev[1] = idet;
        }
        if (!st->UseGausPDF_quantiles[iq] && (st->tq_rec[iq] < st->T0_CF[0] || P_binom < exp(lnP_inf))) recinfo.AllTimesOK = false;
      }
    }//iq
    recinfo.lnP[0] += (st->lnP_start + st->lnP_quantiles);
  }
  return recinfo.lnP[0];
}


//---------------------------------------------------
//---------------------------------------------------
// Signal Likelihood
//---------------------------------------------------
//---------------------------------------------------

double UnivRec::GetSignalLikelihood()
{
  recinfo.nSignalStations = 0;
  recinfo.lnP[1] = 0.;

  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* st = &fstation[idet];
    int type = st->type;
    if (st->mask_signal) continue;
    recinfo.nSignalStations++;

    if (!Check_DX_DL(st)) {
      st->lnP_signal = lnP_inf;
      recinfo.lnP[1] += (st->lnP_signal);
      continue;
    }

    double Rec = st->signalsize;
    double Pred = fsignal[type]->GetSignal(st->r, st->DX, recinfo.logE, recinfo.theta, st->psi, st->density, st->zg, recinfo.Nmu, -1);
    if (st->type == 2) Pred *= (st->DetectorArea / nMD / MDlength / MDwidth);

    double RMS_sh = 0.02 * Pred;
    double rms = sqrt(Pred + RMS_sh * RMS_sh);

    st->signalsize_pred = Pred;
    st->signalsize_rms = rms;

    //-------------------------
    double lnP = 0.;

    //stations saturated
    if (st->isSaturated) {
      double dev = (Rec - Pred) / 2. / rms;
      double Pi = TMath::Erfc(dev) / 2.;
      if (Pi > exp(lnP_inf)) lnP = log(Pi);
      else                   lnP = lnP_inf;

    }
    // Stations with signal: Approx. gaussian
    else if (Rec > 20) {
      if (rms == 0.) lnP = lnP_inf;
      else           lnP = -pow((Rec - Pred) / rms, 2.) / 2. + log(1. / sqrt(2.*M_PI) / rms);
    }
    // Stations with signal: Poisson distrib.
    else if (Rec > 0) {
      double part = double(int(Rec + 0.5));
      double Pi = TMath::Poisson(part, Pred);
      if (Pi < exp(lnP_inf) || std::isnan(Pi) || std::isinf(Pi)) lnP = lnP_inf;
      else                                                      lnP = log(Pi);
    }
    // Silent stations
    else {
      double Pi = 0.;
      for (int i = 0; i < 3; i++) Pi += TMath::Poisson(double(i), Pred);
      if (Pi < exp(lnP_inf)) lnP = lnP_inf;
      else                 lnP = log(Pi);
    }

    if (!st->mask_signal) {
      st->lnP_signal = lnP;
      recinfo.lnP[1] += lnP;
    } else
      st->lnP_signal = 0.;
  }
  return recinfo.lnP[1];
}


//---------------------------------------------------
//---------------------------------------------------
//
//---------------------------------------------------
//---------------------------------------------------

bool UnivRec::Check_DX_DL(RecStation* st)
{
  for (int icomp = 0; icomp < 4; icomp++) {
    if (!UnivParamNS::UseDiffusive[icomp]     && st->DX[icomp] <= 0.)          return false;
    if (!UnivParamTimeNS::UseDiffusive_time[icomp] && st->DL[icomp] <= 0.)          return false;
    double DX0 = fsignal[st->type]->RFunc(st->r, icomp, 5);
    if (st->DX[icomp] <= DX0)                              return false;
  }
  return true;
}

bool UnivRec::SetFitStage(int FitStage)
{
  bool ok = true;

  //  Good Signal reconstruction
  if (FitStage == 1)
    ok = (recinfo.FitStage == 0 && recinfo.FitStatus[0] == 3 && IsTimeRec);
  //  Set PDF flags
  else if (FitStage == 2) {
    ok = (recinfo.FitStage == 1 && IsTimeRec);
    if (ok) SetOptPDF();
  }
  if (ok) recinfo.FitStage = FitStage;

  return ok;
}

bool UnivRec::StationSelection()
{
  //-----------------
  //------Signal  ( WCD )
  //-----------------
  const double vemCut = vemSignalCut;

  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* st = &fstation[idet];
    int type = st->type;
    if (type != 0) continue;

    if (!Check_DX_DL(st)) {
      st->mask_signal = true;
      st->mask_quantiles = true;
      st->mask_starttime = true;
      continue;
    }

    if (st->r > rSignalCut[1] || st->r < rSignalCut[0]) {
      st->mask_signal = true;
      continue;
    }

    double Pred = fsignal[type]->GetSignal(st->r, st->DX , recinfo.logE, recinfo.theta, st->psi, st->density, st->zg, recinfo.Nmu, -1);
    if (Pred < vemCut) {
      st->mask_signal = true;
      continue;
    }

    st->mask_signal = false;
  }

  //-----------------
  //------Signal  ( ASCII/AMIGA )
  //-----------------
  for (int idet_i = 0; idet_i < nRecStations; idet_i++) {
    RecStation_t* st = &fstation[idet_i];
    if (st->type == 0) continue;
    st->mask_signal = false;

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

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

  //-----------------
  //------Time
  //-----------------
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* st = &fstation[idet];

    if (!st->mask_quantiles && (st->r > rQuantileCut[1] || st->r < rQuantileCut[0]))
      st->mask_quantiles = true;

    if (!st->mask_starttime && (st->r > rStartCut[1] || st->r < rStartCut[0]))
      st->mask_starttime = true;
  }

  return TimeSignalOK();
}

//----------------------------------------------------------
// -Chi2 on WCD signals
// -Minimum number of stations positions in TimeLikelihood
//----------------------------------------------------------
bool UnivRec::TimeSignalOK()
{
  //Chi2  (WCD only)
  //---------
  double Chi2 = 0., nChi2 = 0.;
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* st = &fstation[idet];
    if (st->mask_signal) continue;
    if (st->type != 0) continue;

    int type = st->type;

    if (!Check_DX_DL(st)) {
      st->mask_signal = true;
      st->mask_quantiles = true;
      st->mask_starttime = true;
      continue;
    }

    double Pred = fsignal[type]->GetSignal(st->r, st->DX, recinfo.logE, recinfo.theta, st->psi, st->density, st->zg, recinfo.Nmu, -1);

    if (!st->isSaturated) {
      Chi2 += (pow(Pred - st->signalsize, 2.) / Pred);
      nChi2 += 1.;
    }
  }
  if (Chi2 / double(nChi2) > 100.) return false;

  // Minimum number of WCD in time/signal likelihood
  //---------
  int nTimeStations = 0, nSignalStations = 0;
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* st = &fstation[idet];
    if (st->type != 0) continue;
    if (!st->mask_starttime || !st->mask_quantiles) nTimeStations++;
    if (!st->mask_signal && st->signalsize > 0.) nSignalStations++;
  }

  if (IsTimeRec && nTimeStations < nTimeStationsMin) return false;
  if (nSignalStations < nSignalStationsMin) return false;

  return true;
}


bool UnivRec::RejectTimeOutliers()
{

  if (recinfo.TimeMaxdev[0] > 4. && recinfo.nTimeStations > nTimeStationsMin) {
    recinfo.TimeRejected[0] = recinfo.TimeMaxdev[0];
    recinfo.TimeRejected[1] = recinfo.TimeMaxdev[1];

    int idet_max = int(recinfo.TimeMaxdev[1]);
    fstation[idet_max].mask_quantiles = true;
    fstation[idet_max].mask_starttime = true;
    fstation[idet_max].mask_signal = true;

    if (debug == 1) printf("Rejecting station : %d %d \n", fstation[idet_max].id, idet_max);
    return true;
  }
  return false;
}

void UnivRec::SetOptPDF()
{
  for (int idet = 0; idet < nRecStations; idet++) {
    RecStation_t* st = &fstation[idet];
    if (!st->mask_starttime) {
      double Accuracy = (st->type == 2 ? StartAccuracyAMIGA : StartAccuracy[gRecSDE_Upgrade]);
      double rms_model = sqrt(pow(st->start_rms, 2.) - pow(Accuracy, 2.) - pow(InterGPSAccuracy[gRecSDE_Upgrade], 2.));
      st->UseGausPDF_start = (rms_model < sqrt(pow(Accuracy, 2.) + pow(InterGPSAccuracy[gRecSDE_Upgrade], 2.)) ? true : false);
      if (st->start_rec < -3 * st->start_rms) st->UseGausPDF_start = true;
    }
    if (st->mask_quantiles) continue;

    for (int iq = 0; iq < nQ; iq++) {
      if (st->quantile[iq] == UNDEF) continue;

      double rms_model = sqrt(pow(st->tq_rms[iq], 2.) - pow(ParameterizationRMSLimit, 2.) - pow(InterGPSAccuracy[gRecSDE_Upgrade], 2.));
      st->UseGausPDF_quantiles[iq] = (rms_model < sqrt(pow(ParameterizationRMSLimit, 2.) + pow(InterGPSAccuracy[gRecSDE_Upgrade], 2.)) ? true : false);
      if (st->tq_rec[iq] < 0.) st->UseGausPDF_quantiles[iq] = true;
    }
  }
}


//---------------------------------------------------
//---------------------------------------------------
//
//---------------------------------------------------
//---------------------------------------------------


// Load the current guess of the reconstructed parameters into the par array for a new minimization
//   (l,m) are initialized to 0, since after every minimization the direction is centered around the reconstructed one.

//   status=0   Free
//          1   Constrained
//          2   Fixed to current value in recinfo
//          3   Fixed to a value that is set in GetLikelihood routine (it can depend on other parameters)

int UnivRec::InitFCNParameters(std::vector<double>& par , std::vector<double>& epar, std::vector<int>& status , std::vector<bool>& hasLimits)
{
  par.clear();
  epar.clear();
  status.clear();
  hasLimits.clear();


  // Fill par status
  //-----------------
  //----
  if (recinfo.FitStage == -1) {
    //                                               logE                    Nmu              Xmax      T0 l  m     OffsetM_Mu
    const int status_golden[npar] =   {   0,  0,        3,                      0,              3,       3, 2, 2,       3  };
    const int status_sd[npar] =       {   0,  0,        0, (gRecType == 3 ? 0 : 3),              3,       3, 2, 2,       3  };
    for (unsigned int ipar = 0; ipar < npar; ipar++) {
      if (IsGoldenRec) status.push_back(status_golden[ipar]);
      else               status.push_back(status_sd[ipar]);
    }
  }
  //----
  if (recinfo.FitStage ==  0) {
    //                                            logE                  Nmu                        Xmax             T0 l  m   OffsetM_Mu
    const int status_golden[npar] = {   0,  0,        1,                   0,                          1,             3, 2, 2,    3     };
    const int status_sd[npar]    = {   0,  0,        0, (gRecType == 3 ? 0 : 3),                  3,             3, 2, 2,    3     };

    for (unsigned int ipar = 0; ipar < npar; ipar++) {
      if (IsGoldenRec) status.push_back(status_golden[ipar]);
      else               status.push_back(status_sd[ipar]);
    }
  }
  //----
  if (recinfo.FitStage ==  1) {
    //                                       logE  Nmu        Xmax                    T0  l  m      OffsetM_Mu
    const int status_xmax[npar] = {   2,  2,    2,   2, (gRecType == 7 ? 2 : 0),            0, 2, 2, (gRecType == 7 ? 0 : 3) };
    for (unsigned int ipar = 0; ipar < npar; ipar++) status.push_back(status_xmax[ipar]);
  }
  //----
  if (recinfo.FitStage ==  2) {
    //                                                        logE                    Nmu                             Xmax         T0  l  m          OffsetM_Mu
    const int status_xmax[npar] =              {   0,  0, (gRecType == 5 ? 2 : 0), (gRecType == 5 ? 0 : (gRecType == 3 ? 0 : 3)),     0,          0, 0 , 0, (gRecType == 1 ? 0 : 3) };
    const int status_OffsetM_mu_FD[npar] =     {   2,  2,        2,                     2,                                2,          0, 2 , 2,             0         };
    const int status_photon[npar] =            {   2,  2,        2,                     3,                                0,          0, 0 , 0,             3         };
    for (unsigned int ipar = 0; ipar < npar; ipar++) {
      if (gRecType == 6)      status.push_back(status_photon[ipar]);
      else if (gRecType == 7) status.push_back(status_OffsetM_mu_FD[ipar]);
      else                    status.push_back(status_xmax[ipar]);
    }
  }

  // Set Init values
  //--------------------------

  double T0_relative;
  SetSPCoordinates(recinfo.gDetHot);
  RecStation_t* st = &fstation[recinfo.gDetHot];
  T0_relative = (st->GPSnano + st->starttime) - (st->T0_PF + recinfo.T0);
  double val_i[npar] = {  recinfo.xgcore - fstation[recinfo.gDetHot].xg,  recinfo.ygcore - fstation[recinfo.gDetHot].yg,
                          recinfo.logE, recinfo.Nmu, recinfo.Xmax, T0_relative , 0., 0., recinfo.OffsetM_Mu
                       };



  // Fill Limits,Init values
  //--------------------------
  for (unsigned int ipar = 0; ipar < npar; ipar++) {
    par.push_back((status[ipar] == 3 ? 0. : val_i[ipar] / unit[ipar]));
    if (status[ipar] == 0 || status[ipar] == 1) {
      epar.push_back(1.);
      hasLimits.push_back((ipar == 3 || ipar == 4 || ipar == 8 ? true : false)); // Xmax/Nmu/OffsetM_Mu
    } else {
      epar.push_back(0.);
      hasLimits.push_back(false);
    }
  }
  return par.size();
}


//Save Minimization results
void UnivRec::SaveFCNParameters(std::vector<double> par, std::vector<int> status)
{
  //-- Recalculate
  GetLikelihood(par, status);

  //-- Hottest station distance / T0 from hot
  recinfo.rHot = fstation[recinfo.gDetHot].r;

  //-- Calculate parameter errors
  std::vector<double> epar;
  for (unsigned int i = 0; i < par.size(); i++) epar.push_back(0.);
  if (recinfo.FitStage == 0)                   recinfo.FitStatus[0] = GetErrors(par, epar, status);
  if (recinfo.FitStage == 2)                   recinfo.FitStatus[1] = GetErrors(par, epar, status);

  //-- Save errors in internal coordinates
  if (status[0] < 2) recinfo.xgcore_err = epar[0] * unit[0];
  if (status[1] < 2) recinfo.ygcore_err = epar[1] * unit[1];
  if (status[2] < 2) recinfo.logE_err = epar[2] * unit[2];
  if (status[3] < 2) recinfo.Nmu_err = epar[3] * unit[3];
  if (status[4] < 2) recinfo.Xmax_err = epar[4] * unit[4];
  if (status[5] < 2) recinfo.T0_err = epar[5] * unit[5];
  if (status[8] < 2) recinfo.OffsetM_Mu_err = epar[8] * unit[8];

  //-- Save energy estimated in FitStage 0
  if (recinfo.FitStage == 0) {
    recinfo.logE_ref = recinfo.logE;
    recinfo.logE_ref_err = recinfo.logE_err;
    recinfo.theta_ref = recinfo.theta;
    recinfo.azi_ref = recinfo.azi;
  }

  //-- Update direction
  recinfo.theDir_Ref.SetMagThetaPhi(1., recinfo.theta, recinfo.azi);

}


//---------------------------------------------------
//---------------------------------------------------
//
//---------------------------------------------------
//---------------------------------------------------


void UnivRec::Unrotate(double l_p, double m_p, double& theta, double& azi)
{

  double theta_p = asin(sqrt(l_p * l_p + m_p * m_p));
  double azi_p = NormalizeAngle(atan2(m_p, l_p));

  TVector3 theDir_p;
  theDir_p.SetMagThetaPhi(1., theta_p, azi_p);
  theDir_p.RotateY(recinfo.theDir_Ref.Theta());
  theDir_p.RotateZ(recinfo.theDir_Ref.Phi());

  theta = theDir_p.Theta();
  azi = NormalizeAngle(theDir_p.Phi());
}


void
UnivRec::SetT0FromHot(bool UseT0_relative, double T0_relative)
{
  SetSPCoordinates(recinfo.gDetHot);
  RecStation_t* const st = &fstation[recinfo.gDetHot];

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

  if (UseT0_relative)
    recinfo.T0 = (st->GPSnano + st->starttime) - (st->T0_PF + T0_relative);
  else {
    //---Expected signal
    double vemTot_p = 0., fcomp[4];
    for (int icomp = 0; icomp < 4; icomp++) {
      fcomp[icomp] = fsignal[st->type]->GetSignal(st->r, st->DX, recinfo.logE, recinfo.theta, st->psi, st->density, st->zg, recinfo.Nmu, icomp);
      vemTot_p += fcomp[icomp];
    }
    for (int icomp = 0; icomp < 4; icomp++)
      fcomp[icomp] /= vemTot_p;

    // Start time
    double UnitPerMuon = 0;
    if (st->type == 0)
      UnitPerMuon = 1. / (cos(recinfo.theta) + 2.*tankh * sin(recinfo.theta) / M_PI / tankr);
    else if (st->type == 1)
      UnitPerMuon = 1. / cos(recinfo.theta);
    else if (st->type == 2)
      UnitPerMuon = 1.;
    double nmuDetector = fcomp[0] * st->signalsize / UnitPerMuon;
    if (nmuDetector <= 1.)
      nmuDetector = 1.;
    double pred = 0, rms = 0;
    double corr = ftime[st->type]->tStartCorrection(st->r, recinfo.logE, recinfo.theta, gRecSDE_Upgrade);
    /*double P_extreme =*/ ftime[st->type]->tFirstPDF(st->DL[0], st->r, recinfo.logE, st->psi, recinfo.theta, st->T0_CF[0], nmuDetector, -1., false, pred, rms);
    recinfo.T0 = (st->GPSnano + st->starttime) - (st->T0_PF + pred + corr);
  }
}


// Those parameters that are fixed use the values stored in recinfo except the exceptions found below.
double UnivRec::GetLikelihood(std::vector<double> par, std::vector<int> status)
{
  if (par.size() != npar) {
    printf("ERROR (GetLikelihood): par size %ld\n", par.size());
    exit(1);
  }
  for (unsigned int ipar = 0; ipar < npar; ++ipar)
    if (std::isnan(par[ipar]) || std::isinf(par[ipar])) {
      printf("ERROR (GetLikelihood): par %d  ->%lf \n", ipar, par[ipar]);
      return lnP_inf * 4 * nRecStations;
    }
  recinfo.nFCNcalls++;

  //-- Transform parameters
  if (status[6] < 2 && status[7] < 2)
    Unrotate(par[6]*unit[6], par[7]*unit[7], recinfo.theta, recinfo.azi);

  if (status[0] < 2) recinfo.xgcore = par[0] * unit[0] + fstation[recinfo.gDetHot].xg;
  if (status[1] < 2) recinfo.ygcore = par[1] * unit[1] + fstation[recinfo.gDetHot].yg;

  if (status[2] < 2) recinfo.logE  = par[2] * unit[2];
  else if (status[2] == 3) {
    if (gRecType == 5 && recinfo.FitStage > 0) recinfo.logE = recinfo.logE_ref;
    else if (IsGoldenRec) recinfo.logE = augerinfo.logE_FD;
    else {
      printf("ERROR (GetLikelihood): setting energy \n");
      exit(1);
    }
  }

  if (status[4] < 2)         recinfo.Xmax  = par[4] * unit[4];
  else if (status[4] == 3) {
    if (gRecType == 0 || gRecType == 7)      recinfo.Xmax = augerinfo.Xmax_FD;
    else                                         recinfo.Xmax = fcalib->GetMeanXmax(recinfo.logE);
  }

  if (status[3] < 2)         recinfo.Nmu = par[3] * unit[3];
  else if (status[3] == 3) {
    double Xmax = (gRecType == 4 || gRecType == 1 ? recinfo.Xmax : -100.);
    recinfo.Nmu = fcalib->GetMeanNmu(recinfo.logE, recinfo.theta, Xmax);
  }

  if (status[8] < 2)        recinfo.OffsetM_Mu = par[8] * unit[8];
  else if (status[8] == 3)  recinfo.OffsetM_Mu = fcalib->GetOffsetM_Mu(recinfo.logE, recinfo.theta, recinfo.Xmax);

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

  if (status[5] < 2)        SetT0FromHot(true, par[5]*unit[5]);
  else if (status[5] == 3)  SetT0FromHot(false, UNDEF);

  //-- Set coordinates,DX,T0_CF,T0_PF
  SetAllSPCoordinates();

  recinfo.lnP[2] = 0.;

  //  Golden hybrid constrains
  //--------------------------

  //--- Energy constrain
  if (status[2] == 1) {
    if (!IsGoldenRec) {
      printf("ERROR (GetLikelihood): constraining energy \n");
      exit(1);
    }
    double rec = recinfo.logE;
    double mean = augerinfo.logE_FD, rms = augerinfo.logE_FD_err;
    double lnP_Energy = (-pow((mean - rec) / rms, 2.) / 2. + log(1. / sqrt(2.*M_PI) / rms));
    recinfo.lnP[2] += lnP_Energy;
  }

  //--- Xmax constrain
  if (status[4] == 1) {
    if (!IsGoldenRec) {
      printf("ERROR (GetLikelihood): constraining Xmax \n");
      exit(1);
    }
    double rec = recinfo.Xmax;
    double mean = augerinfo.Xmax_FD, rms = augerinfo.Xmax_FD_err;
    double lnP_Xmax = (-pow((mean - rec) / rms, 2.) / 2. + log(1. / sqrt(2.*M_PI) / rms));
    recinfo.lnP[2] += lnP_Xmax;
  }

  //--- Nmu constrain
  if (status[3] == 1) {
    double rec = recinfo.Nmu / fcalib->GetMeanNmu(recinfo.logE, recinfo.theta, -100.);
    double mean = 1., rms = 0.2;
    double lnP_Nmu = (-pow((mean - rec) / rms, 2.) / 2. + log(1. / sqrt(2.*M_PI) / rms));
    recinfo.lnP[2] += lnP_Nmu;
  }

  //--- Signal/Time Likelihood
  double lnP = 0.;
  lnP = recinfo.lnP[2];

  if ((gRecType == 6 || gRecType == 7) && recinfo.FitStage > 0)          recinfo.lnP[1] = 0.;
  else                                                            lnP += GetSignalLikelihood();

  if (recinfo.FitStage > 0)                                       lnP += GetTimeLikelihood();

  /*
  printf(" %5.2lf %5.2lf %5.2lf %5.2lf %5.2lf %5.2lf %5.2lf %5.2lf  | lnP=%5.2lf | %d  \n",recinfo.xgcore/1.e2,recinfo.ygcore/1.e2,recinfo.logE,recinfo.Nmu,
   recinfo.Xmax, recinfo.T0, recinfo.theta*180./M_PI, recinfo.azi*180./M_PI
   ,lnP,status[6]);
  */

  if (std::isnan(lnP) || std::isinf(lnP)) {
    printf("ERROR (GetLikelihood): lnP=%lf %lf %lf | %lf  \n", recinfo.lnP[0], recinfo.lnP[1], recinfo.lnP[2], lnP);
    return lnP_inf * 4 * nRecStations;
  }

  return lnP;
}

// Covariance Matrix
int UnivRec::GetErrors(std::vector<double> par,  std::vector<double>& epar , std::vector<int> parStatus)
{
  // If F=-log(L) is quadratic F=F_0+0.5 (t/d)^2, the second derivative is 1/d^2
  //

  int nFCNcalls_start = recinfo.nFCNcalls;

  bool ok_range = true;

  std::vector<bool> isFixed;
  for (unsigned int ipar = 0; ipar < par.size(); ipar++)
    isFixed.push_back(parStatus[ipar] < 2 ? false : true);

  int ndim = 0;
  for (unsigned int ipar = 0; ipar < par.size(); ipar++)
    if (!isFixed[ipar]) ndim++;

  // Save Error guess 1D likelihood scan ( Step is decreased to guarantee at least a good minimum in all 1D scan )
  //----------------------------------------------
  if (debug > 0) printf("1D scan of parameter errors \n\n");
  double errDef_i[npar] = { 0.5, 1.20, 1.35, 2.40, 3.00, 3.60, 4.20 , 4.75, 5.30 };
  for (unsigned int i = 0; i < par.size(); i++) {
    if (isFixed[i]) continue;
    double rangeX[2], rangeY[2];
    double d[2], edm;
    bool ok = ScanRange(i, UNDEF, par, parStatus, rangeX, rangeY, d , 0. , errDef_i[ndim - 1], edm);
    if (!ok) ok_range = false;
    epar[i] = (fabs(d[0]) + fabs(d[1])) / 2.;
    if (debug > 0) printf(" ipar=%d %lf (%lf,%lf)  edm=%5.2lf ok=%d \n", i, par[i]*unit[i], (par[i] + d[0])*unit[i], (par[i] + d[1])*unit[i], edm, ok);
  }
  if (!ok_range) return 1;

  // Calculate d_side[i][j]
  //----------------------------------------------
  double errDef = 1.0;
  if (debug > 0) printf("\n\n2D scan of parameter errors \n\n");
  recinfo.edm = 0.;
  double d_side[par.size()][par.size()][2];
  for (unsigned int i = 0; i < par.size(); i++)
    for (unsigned int j = 0; j <= i; j++) {
      double range_i[2], range_j[2];
      if (isFixed[i] || isFixed[j]) {
        d_side[i][j][0] = 0.;
        d_side[i][j][1] = 0.;
        continue;
      }
      double edm;
      bool ok = ScanRange(i, j, par, parStatus, range_i, range_j, d_side[i][j] , (i == j ? 0. : M_PI / 4.) , errDef, edm);
      if (!ok) ok_range = false;
      if (edm > 0. && edm > recinfo.edm) recinfo.edm = edm;
      if (UnivRecNS::debug > 0) {
        printf("(i,j)=(%d,%d) ok=%d d=(%lf,%lf) ", i, j, ok, d_side[i][j][0], d_side[i][j][1]);
        if (i == j)  printf(" (%lf,%lf)  \n", (par[i] + d_side[i][i][0])*unit[i], (par[i] + d_side[i][i][1])*unit[i]);
        else         printf("\n");
        if (edm > 0.) printf("Warning %lf \n", edm);
      }
    }
  if (!ok_range)  return 1;

  // Approximate d[i][j] as the left/right average
  //----------------------------------------------
  double d[par.size()][par.size()];
  for (unsigned int i = 0; i < par.size(); i++)
    for (unsigned int j = 0; j < par.size(); j++) {
      if (j <= i) d[i][j] = (fabs(d_side[i][j][0]) + fabs(d_side[i][j][1])) / 2.;
      else        d[i][j] = (fabs(d_side[j][i][0]) + fabs(d_side[j][i][1])) / 2.;
      d[i][j] /= sqrt(errDef * 2.);
      if (isFixed[i] || isFixed[j]) continue;
    }

  // Covariance matrix
  //---------------------

  TMatrixDSym G(ndim);

  int ir = 0;
  for (unsigned int i = 0; i < par.size(); i++) {
    if (isFixed[i]) continue;
    int ic = 0;
    for (unsigned int j = 0; j < par.size(); j++) {
      if (isFixed[j]) continue;
      if (i == j)  TMatrixDRow(G, ir)[ic] = 1. / d[i][j] / d[i][j];
      else          TMatrixDRow(G, ir)[ic] = 1. / d[i][j] / d[i][j] - 1. / d[i][i] / d[i][i] / 2. - 1. / d[j][j] / d[j][j] / 2. ;
      ic++;
    }
    ir++;
  }

  // Eigen value analysis
  //---------------------
  const TMatrixDSymEigen eigen(G);

  TVectorD eigenVal = eigen.GetEigenValues();
  double MaxEigenVal = -1.e10, MinEigenVal = 1.e10;
  for (int i = 0; i < ndim; i++) {
    if (eigenVal[i] < MinEigenVal) MinEigenVal = eigenVal[i];
    if (eigenVal[i] > MaxEigenVal) MaxEigenVal = eigenVal[i];
  }
  recinfo.MinEigenVal = MinEigenVal;
  if (UnivRecNS::debug > 0)
    printf(" Max/Min eigen value  %5.2lf %5.2lf \n", MaxEigenVal, MinEigenVal);

  if (MinEigenVal <= 0.) return 2;

  // Error calculation
  //-------------------------------

  TMatrixD ErrMatrix = G;
  ErrMatrix.Invert();

  bool ok_errors = true;
  for (int i = 0; i < ndim; i++)
    if (TMatrixDRow(ErrMatrix, i)[i] < 0.) ok_errors = false;
  if (!ok_errors) return 2;

  if (UnivRecNS::debug > 0)
    printf(" \n\n Errors : ");
  int ii = 0;
  for (unsigned int ipar = 0; ipar < par.size(); ipar++) {
    if (isFixed[ipar])   epar[ipar] = 0.;
    else                 {
      epar[ipar] = sqrt(TMatrixDRow(ErrMatrix, ii)[ii]);
      ii++;
    }
    if (UnivRecNS::debug > 0)
      printf(" %5.2lf %5.2lf %5.2lf (%5.2lf) | ", epar[ipar]*unit[ipar], d_side[ipar][ipar][0]*unit[ipar], d_side[ipar][ipar][1]*unit[ipar], epar[ipar]);
  }
  if (UnivRecNS::debug > 0) printf("\n");


  // Print
  //--------------------------
  if (UnivRecNS::debug > 0) {
    std::string MatrixName[2] = { " Covariance Matrix", "Error Matrix"};
    for (int iloop = 0; iloop < 2; iloop++) {
      printf("%s \n", MatrixName[iloop].c_str());
      printf("-----------------------------\n");
      for (int ir = 0; ir < ndim; ir++) {
        for (int ic = 0; ic < ndim; ic++) {
          double val = 0.;
          if (iloop == 0) val = TMatrixDRow(G, ir)[ic];
          if (iloop == 1) val = TMatrixDRow(ErrMatrix, ir)[ic];
          printf(" %5.2lf ", val);
        }
        printf("\n");
      }
      printf("\n");
    }
    printf(" FCN calls calculating errors %d Ok=%d \n", ok_range, recinfo.nFCNcalls - nFCNcalls_start);
  }

  return 3;

}