#ifndef SDRECLPNHE_H
#define SDRECLPNHE_H

#include <unistd.h>
#include <vector>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <math.h>

#include "TMinuit.h"
#include "TF1.h"
#include "IoSd.h"
#include "Es.h"
#include "TReconstruct.h"
#include "ToolsLib.h"

#define SIGMA_T_2 (35.*35.)     ///< old time variance
#define XVGD (875.5)            ///< Xmax at ground
#define DISTREF 1000
#define DISTREFINFILL 450
#define SINTH2MAX 1.05
#define CHI2LDFMAX 50
#define DCHI2LDFMAX 50
#define DIST_BARY_CORE_MAX 500
#define SILSIGTHRESH 3.         ///< signal threshold for silent stations
#define SILSIGTHRESHNEW 1.      ///same with new triggers 
#define SIGTHRESH 0.2

namespace SdLpnhe {
enum TVarianceTs        {eR4, eT50, eT50Mod, eT70, eT90};
enum TGeometry  {eGeomNone, ePlane, eCurve, eVCurve, eGeomOptimal};
enum TLDF               {eLDFNone, eLog, eLogXmax, eNKG, eLDFOptimal};
enum TTimeRec   {eTimeNone, eRise, eFall};
enum TPrimary   {eProton, eIron, eOther};
enum TFit               {eGeom, eLDF, eGlobal};
}

using namespace SdLpnhe;

/// LPNHE reconstructor. Based on the reconstruction of Er v2*
class SdRecLpnhe: public TReconstruct {

 private:
  // Singleton stuff
  SdRecLpnhe();
  ~SdRecLpnhe();
  static SdRecLpnhe * _instance;
  
  // internal parameters
  int _errflag;                    ///< reconstruction flag
  vector <double> fWgt;  /// Weight of station [i] in time fit; default 1.0
  void SetWgt(TCalibStation *st, double w);
  // Internal Recontruction functions
  int _GeometricalRec(TGeometry geom, bool allowChanges = false);
  int _LDFRec    (TLDF ldf);
  int _CompleteRec(TGeometry geom, TLDF ldf);
  int _Reconstruct(TGeometry geom, TLDF ldf, bool global);

  // Internal auxiliary functions
  void   _EstimateLDFPar(TLDF ldf = eLog);
  void   _EstimateSRef();
  void   _EstimateXYZTCore();
  void   _CalculateXmax();
  bool   _Minimize(TMinuit*);
  bool   _FitBeta();
  int    _CountDoublets();
  int    _CompactConfig();
  int    _EstimateUV();
  int    _OptimizeStationSelection(int loop=1);
//  int    _OptimizeStationSelection2();
  int    _TimeRec(TTimeRec);
  double _StatDt(int, int, TGeometry);
  double _StatT(int, int, TGeometry);
  double _Stat2Core(int);
  double _MatPlanFit(int, int&);
  double _MatCurvedFit(int, int&);
  double _TimeResiduals(TGeometry, bool&);
  double _EstimateTCore();
  double _EstimateCurvature(double signal,double costheta,bool infillrec);
  double _GetResults(TMinuit*, TFit, TGeometry, TLDF);
  TMinuit* _NewMinuit(TFit, TGeometry, TLDF);
  TGeometry _Best(TGeometry);
  TLDF     _Best(TLDF);
  /// This function (inherited fomr TReconstruct)) is called automatically by TErEvent.Reconstruct().
  /// The reconstruction behavior is driven by the global parmaters as set with the  Set-functions above.
  bool    _Reconstruct();


public:
  
  static SdRecLpnhe* Instance();
  double Wgt(TCalibStation* t);  /// Weight of station [i] in time fit; default 1.0
  //double T50(int stat, bool doIt = false);  /// Selected stations  average T50
  // Set Functions to drive the reconstruction
  void SetRec(TGeometry rec1 = eGeomOptimal, TLDF rec2 = eLDFOptimal);
  void SetAsymmetry(bool asym = true);
  void SetPrimary(TPrimary prim = eProton);
  void SetEnergyConverter(string e = "LOG");
  void SetVerbosity(int verb = 0);
  void SetCore(double xc = 0.0, double yc = 0.0, double zc = 1400.0);
  void SetSilent(bool silent = true, double Rcut = 5000.0);
  void SetForce(bool force = false);
  void SetVarianceModel(TVarianceTs  var = eT50);
  void SetOfficial(bool official = false);
  void SetOptimize(bool optimize = true);
  void SetGlobal(bool = true);
  void SetDefault(); /// Puts back everything to original default value.

  void SetErEvent(TErEvent *ev); /// Called by ErEvent constructor to associate ErEvent to reconstructor.
  int InitReconstruction();///Called before any independent reconstruction process
  double CorrectTime(int);
  double TimeVariance(int, double, double, TGeometry);
  // Xmax stuff
  int XmaxOK();

  // User driven functions
  bool Reconstruct() {
    return fErEvent->Reconstruct(this);
  } // theLpnheRec().Reconstruct() is equivalent to TErEvent.Reconstruct(theLpnheRec());
  int  GeometricalRec(TGeometry geom = eGeomOptimal);///User-callable angular reconstruction. options ePlan, eCurve, eVCurve.
  int  LDFRec(TLDF ldf = eLDFOptimal);///User-callable LDF reconstruction. Option eLog, eXmax
  int  CompleteRec(TGeometry geom = eGeomOptimal, TLDF ldf = eLDFOptimal);///User-callable complete reconstruction; Options as above
  int  HorizontalReconstruction();

  /// dumps for the reconstruction results
  void StationDump(int,  TGeometry);
  void HeaderDump();
  void RecDump(TLDF);
  void HorizontalReconstructionDump(TGeometry = eGeomOptimal);
  void AngularDump(TGeometry, TLDF);
  void  GenResPlot(TGeometry);
  void LDFDump(TGeometry, TLDF);
  void HorizontalResultsDump();
  void CompleteDump(TGeometry, TLDF);
  
  double XmaxfromRiseTime;///< Xmax estimated from the risetime at r = 1000 m
  double XmaxfromFallTime;///< Xmax estimated from the falltime at r = 1000 m

  // max time residuals
  //  unsigned int RESMAX, RESPLANMAX, RESPLANSTATMAX, RESSTATMAX;
  
  // parameters used in fits
  double fn_sRef;
  double ftCore;
  double fdBeta,fdGamma;
  double ftau, fdtau;
  double falpha, fdalpha;

  //for timeResiduals fit 
  vector <double> fRPlotTimes;
  vector <double> fRPlotDts;
  vector <double> fRPlotDist;
  vector <double> fRPlotDtsErr;
  vector <double> fRPlotDistErr;
 vector <int> fRPlotIds;

  // chi2s...
  double angularChi2, ldfChi2, globalChi2;
  double residu;
  unsigned int LDFFitOK, GlobalFitOK, ndof;
  
  ostream & put(ostream & s, char *option);
};

SdRecLpnhe& theLpnheRec();

ostream& operator<<(ostream & s, SdRecLpnhe & r);
const char *ErSetFormat(char *option);

#endif