#include "cdas_univ_config.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vector>
#include <fstream>

#ifdef CDAS
#include <Atmosphere.h>
#include <Er.h>
#ifdef HAVE_IOMD
#include <MdEvent.h>
#endif
#else
#include <ToyShower.h>
#endif

#include <TVector3.h>
#include <UnivParam.h>
#include <UnivParamTime.h>

#ifndef UNDEF
#define UNDEF -100
#endif

namespace UnivRecNS {

  const int debug=0;

  //--- Inherited from DetectorResponse.h ( change accordingly)

  // Auger Water Cherenkov
  const double tankh = 120.;
  const double tankr = 180.;

  // Scintillator 3 m^2
  const double scinlength= 170.;
  const double scinwidth = 170.;
  const double scinheight= 1.;


  //--- MC/Auger settings
  const double fProton[3]={ 0.5, 0.1, 0.9 };

  const double Xmax_PhotonSearch=1050.;
  const double Nmu_PhotonSearch=0.15;
  const double OffsetT1Mu_PhotonSearch=0.;

  const double edm_tolerance=0.1;

  //--- Auger specific settings
  //  [0] Standard
  //  [1] FD energy -15%
  //  [2] SD DT1Mu calib
  //  [3] C_XMaxMu=0.85
  //  [4] rho=1.06 kg/m3
  //  [5] (m.s) modulation Off
  //  [6]-[9]  Use only Eye 1/2/3/4
  //  [10] FD energy +15%
  //  [11] FD Xmax  +10 g/cm2
  //  [12] FD Xmax  -10 g/cm2

  //  [6]/[7]/[8]/[9]/[11]/[12]  ( Only has an effect when gRecDataSet=1, gRecType==0 )

  const int nSys=13;
  const double C_XmaxMu_Auger[nSys]=  { 0.925, 0.925, 0.925, 0.85, 0.925, 0.925 ,   0.925, 0.925, 0.925, 0.925 ,    0.925, 0.925, 0.925 };
  const double OffsetT1Mu_Auger[nSys]={ 17.50, 16.00, 10.5 , 20.0 ,17.5, 17.5,    17.5, 17.5, 17.5, 17.5,           18.0,  14.0 , 22.0  };
  //const double OffsetT1Mu_Auger[nSys]={ 7.50, 16.00, 10.5 , 20.0 ,17.5, 17.5,    17.5, 17.5, 17.5, 17.5,           18.0,  14.0 , 22.0  };
  const double fNmuSys_Auger[nSys]   ={ 1., 0.85,  1., 1., 1.,1.,  1.,1.,1.,1. ,   1.15 , 1., 1. };
  const double fEnergySys_Auger[nSys]={ 1., 1.15,  1., 1.,1.,1.,   1.,1.,1.,1. ,   0.85 , 1., 1. };
  const double fXmaxSys_Auger[nSys]=  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,0., +10., -10. };


  //---- MC specific settings
  const double SmearingAngle=0.25; // Smearing angle for MC (Hybrid angular resolution)

  // [im][p,Fe,Mixtures]
  const double C_XmaxMu_Models=0.925;
  const double OffsetT1Mu_Models[2][6]={  {-2.600,0.500, 0.75, 0.35,-1.00, 1.50}, { 2.50, 12.50, 6.75,12.00, 3.00, 5.50} };

  //-----
  const int nQ=3;
  const double lnP_inf=-200.;

  const int nDetectorType=3;

  const unsigned int npar=9;
  const double unit[npar]={ 20.e2, 20.e2, 0.03, 0.20, 30.,  100., 2.e-3, 2.e-3, 1.0 };
  const double low[npar]= {  -1500.e2,  -1500.e2, 17.0, 0.05,  200. , -500., -sin(25.*M_PI/180.), -sin(25.*M_PI/180.),-50.  };
  const double high[npar]={   1500.e2,   1500.e2, 22.5, 3.5 , 1200. ,  500.,  sin(25.*M_PI/180.),  sin(25.*M_PI/180.),+50.  };
  const std::string ParName[npar]={ "xgcore_relative","ygcore_relative","logE", "Nmu" , "Xmax", "T0", "l", "m", "DT1Mu" };

  const int nAoP=4;
  const double AreaOverPeak_Run[nAoP]={3.60,3.60,4.30,2.93};

  const double SaturationLevelVEM_WCD[2]={200.,2000.};
  const double SaturationLevelVEM_Scin[2]={200.,800.};

  const double InterGPSAccuracy[2]={10.,4.};
  const double ParameterizationRMSLimit=5.;
  const double StartAccuracy[2]={6.25,2.0};


  //---
  const double rQuantileCut[2]= { 0.e2,2000.e2};
  const double vemTimeCut[3]={ 30., 5., 5. };      // [0] T1 [1] T10 [2] T50

  //---
  const double rSignalCut[2]= {0.e2,2500.e2};
  const double vemSignalCut=3.,vemSignalCut_low=1.;      // WCD  ( ASCII and AMIGA are always used if WCD is flagged as useable )
  const double vemSignalCut_NmuSD=10.;                   //    gRecType=3

  const int nTimeStationsMin=4;
  const int nSignalStationsMin=3;

  const double f_md=0.6;  // The difference between Nmu >~250 MeV and Nmu>~1 GeV
  const double AreaMD= 10.e4;// [cm²]
  const double SatLevelMD=200.;//

  //----------------------------------
  //  Shower Reconstruction results
  //----------------------------------

  typedef struct RecInfo {

    int id;
    double AugerDensity;
    int iatm;

    double xgcore,xgcore_err;
    double ygcore,ygcore_err;
    double zgcore;

    double logE, logE_err;

    double Xmax, Xmax_err;

    double Nmu, Nmu_err;

    double theta;
    double azi;

    double T0,T0_err;

    double logE_ref,logE_ref_err;
    double theta_ref,azi_ref;

    double C_XmaxMu;
    double OffsetT1Mu,OffsetT1Mu_err;

    double RA,Dec;

    int FitStatus[2];
    float MinEigenVal,edm;

    double lnP[3];  //Time/Signal/Constrains

    int nStationsScin,nStationsWCD;          // non-accidental triggered stations above signal threshold ( Different threshold for Xmax/Signal rec. )
    int nTimeStations, nSignalStations;

    int seed [3]; // The ids of the WCD detectors that form the triangle with the highest signal
    TVector3 theDir_Ref;

    int FitStage;

    //  Time reconstruction information
    bool AllTimesOK;           // All start time and time quantiles positive
    double TimeRejected[2];    // Rejected station  [0] Maximum standard deviation in time likelihood [1] Id of the station
    double TimeMaxdev[2];      // [0] Maximum standard deviation in time likelihood [1] Id of the station

    // Additional information
    double rHot;
    int gDetHot;
    int nFCNcalls;
    bool IsSaturated[nDetectorType]; // WCD / Scintillator

  } RecInfo_t;

  //----------------------------------
  //  MC/FD info
  //----------------------------------
  typedef struct RefInfo{

    int im,ip,ie,ith,ish,iatm;

    double xgcore,ygcore,zgcore;
    double theta,theta_fluct;
    double azi,azi_fluct;
    double logE,logE_err,logE_fluct;
    double Xmax,Xmax_err,Xmax_fluct;
    double Nmu,T0;
    double logE_SD,logE_SD_err; //
    double theta_SD,azi_SD;     //
    double RA, dec;

    bool isGoodSD; // (Data) relax T5
    bool isGoodFD; // (Data) fiducial cuts
  } RefInfo_t;


  //----------------------------------
  //  Rec Station
  //----------------------------------
  typedef struct RecStation {
    int type; // [0] WCD [1] Scin [2] MD

    int id;

    double xg,yg,zg;
    double GPSnano;
    double AreaOverPeak;

    bool mask_quantiles;
    double quantile[nQ],tquantile[nQ];
    double tq_rec[nQ],tq_pred[nQ],tq_rms[nQ],tq_corr[nQ];
    bool UseGausPDF_quantiles[nQ];
    double lnP_quantiles;

    bool mask_starttime,  dummy;
    double starttime;
    double start_rec,start_pred,start_rms,start_corr;
    bool UseGausPDF_start;
    double lnP_start;

    bool mask_signal;
    double signalsize;
    double signalsize_pred,signalsize_rms;
    double lnP_signal;

    bool isSaturated;

    double r,psi,DX[4],T0_CF[4],T0_PF;
  } RecStation_t;



  //----------------------------------
  //  UnivRec
  //----------------------------------

  class UnivRec  {

  public:
    UnivRec();
    ~UnivRec() {}

    bool InitRec_MC( int RecType, bool RecInfill, bool RecSDEUpgrade, int RecAoP, int RecDetector, int RecMixture);
    bool InitRec_Data(  int RecType, bool RecInfill,  int RecDataSet, int RecSys, int year, std::string &FileName );
#ifndef CDAS
    bool InitEvent(ToyShowerNS::ToyShower *fevt);
#else
    bool InitEvent(TErEvent *fevt
#ifdef HAVE_IOMD
    , md::Event *mdevt=NULL
#endif
);
    bool InitRec_CDAS(int RecType, bool RecInfill, int RecSys);
#endif
    bool InitRecParameters(bool UseMC);
    bool InitRecParameters(bool UseMC, double OffsetT1Mu_Custom );

    float DistPoints(float x1, float y1, float x2, float y2);
    float NormalizeAngle(const float angle);

    void SetSPCoordinates(int idet);
    void SetAllSPCoordinates();

    bool GetRecEstimate(int itype);
    bool PlaneFit( bool CurvCorr, bool AltitudeCorr, bool RecSeed, int itype );
    bool PlaneFit(std::vector<double> fX,std::vector<double>fY,std::vector<double>fT,std::vector<double>fW);
    void GetRecSeed(int itype);

    void PrintRecInfo( bool PrintResiduals);
    void WriteTextFile(FILE *fp);
    bool RWRecinfo(FILE *fp, bool RWStations, bool readFlag);
    void GetChi2(double *Chi2, int *nChi2);
    void GetResiduals( std::vector<double> &mean, std::vector<double> &rms, std::vector<double> &rec, int optY, std::vector<double> &fx, int optX,bool OnlyNotSaturated );
    void GetResiduals( std::vector<double> &mean, std::vector<double> &rms, std::vector<double> &rec, int optY, std::vector<double> &fx, int optX );
    void GetMeanRMS_AoP(double &mean, double &rms);
    bool 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);
    bool 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);
    void 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);
    bool IsGoodRec();
    bool IsGoodRec(bool CutBadYears);
    bool CheckEnergyCut(bool PrintOut);

    void SetT0FromHot(bool UseT0_relative, double T0_relative );

    double GetTimeLikelihood();
    double GetSignalLikelihood();

    bool SetFitStage(int FitStage);
    bool StationSelection();
    bool TimeSignalOK();
    bool RejectTimeOutliers();
    void SetOptPDF();

    void SaveFCNParameters(std::vector<double> par, std::vector<int> status );
    int  InitFCNParameters(std::vector<double> &par, std::vector<double> &epar , std::vector<int> &status ,std::vector<bool> &hasLimits,bool debugXmax);

    void Unrotate(double l_p, double m_p, double &theta, double &azi);
    int GetErrors(std::vector<double> par,  std::vector<double> &epar, std::vector<int> parStatus  );
    double GetLikelihood(std::vector<double> par, std::vector<int> status);

    //---
    RecStation_t *GetRecStation(int idet) { return &fstation[idet]; }
    RecInfo_t *GetRecInfo()  { return &recinfo; }
    RefInfo_t *GetRefInfo()  { return &refinfo; }
    unsigned int GetGPSsec() { return GPSsec;   }
    unsigned int GetnRecStations() { return nRecStations;}
    bool isUpgrade() { return gRecSDE_Upgrade;}
    bool isMC()      { return gRecMC; }
    void GetXmaxRange(double &Xmax0, double &Xmax1){ Xmax0=XmaxLow; Xmax1=XmaxHigh;}
    void SetXmaxRange(double Xmax0, double Xmax1){XmaxLow=Xmax0; XmaxHigh=Xmax1;}
    bool IsTimeReco(){ return IsTimeRec; }
    void SetIsSignalCut_low(bool flag) { IsSignalCut_low=flag; }

  private:
    std::vector<RecStation_t> fstation;
    int nRecStations;
    RefInfo_t refinfo;
    RecInfo_t recinfo;

    int gRecType;
    bool gRecMC,gRecInfill;
    bool gRecSDE_Upgrade;
    int gRecAoP;
    int gRecMixture;
    bool gDetectorTypeMask[nDetectorType];
    int gRecDataSet; //[0] SD standard/infill [1] Golden standard/infill [2] Golden standard/infill (no FOV) [3] SD high energy

    bool IsTimeRec;
    bool IsSignalCut_low; //   FitStage=-1 use all triggered stations and silent in 3 rings around hot   /   Fitstage=0  vemSignalCut=1 [vem]

    int gRecSys;

    unsigned int GPSsec;
    double XmaxLow,XmaxHigh;

    UnivParamNS::UnivParam *fsignal[nDetectorType];
    UnivParamTimeNS::UnivParamTime *ftime[nDetectorType];
    AtmosphereNS::Atmosphere *fatm;
  };

}