#ifndef __UnivCalibConstants_h_
#define __UnivCalibConstants_h_

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fstream>

#define UNDEF -100

namespace UnivCalibConstantsNS {

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

  static const double logE_Calib = 19.;
  static const double deltaXmax_pFe_Calib = 90.;

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

  static const double MeanXmax_pFe_Calib = 755., ER_pFe_Calib = 50.;
  static const double OffsetM_Mu_Calib[2] = { 0., 0. };

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

  static const int nSys = 15;
  static const double OffsetM_Mu_Auger[2][nSys] =  { { 0.11, 0.10, 0.11, 0.08,  0.13,  0.10,  0.11, 0.11,  0.11,  0.11, 0.11, 0.11, 0.11  , 0.09, 0.13 },  // Overall
    { 0.  , 0., 0., 0., 0.,  0., 0., 0., 0., 0., 0., 0., 0., 0., 0. }
  };   // Theta dependence
  // NOTE: gRecSys=13/14 corresponds to more than 2 Sigma of the fitted OffsetM_Mu value from Golden events ( in the rest of the cases is 1 sigma )

  static const double MeanXmax_pFe_Auger = 760., ER_pFe_Auger = 50.;

  static const double fNmuSys_Auger[nSys]   =           { 1., 0.85, 1.15, 1., 1.,   1., 1., 1., 1., 1. ,        1., 1., 1.  , 1., 1. };
  static const double fEnergySys_Auger[nSys] =           { 1., 1.15, 0.85, 1., 1.,   1., 1., 1., 1., 1. ,        1., 1., 1.  , 1., 1. };
  static const double fXmaxSys_Auger[nSys] =             { 0., 0., 0., +10., -10.,     0., 0., 0., 0., 0.,         0., 0., 0.  , 0., 0. };
  static const double Offset_MeanXmax_pFe_Auger[nSys] =  { 0., 0., 0., 0., 0.,         +15., -15., 0., 0., 0.,       0., 0., 0.  , 0., 0. };
  static const double Offset_ER_pFe_Auger[nSys] =        { 0., 0., 0., 0., 0.,           0., 0., +10., -10., 0.,     0., 0., 0.  , 0., 0.};

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

  static const int nRecMixtures = 6;
  static const int nModels = 4;

  static const double MeanXmax_pFe_Models[nModels] = { 745., 767., 743., 757. }, ER_pFe_Models[nModels] = { 45., 60., 55., 57. };

  static const double OffsetM_Mu_Models[2][nModels] = { {  0.038, 0.091, 0.000, 0.000  } ,         // Overall
    { -0.011, -0.041, 0.000, 0.000  }
  };        // Theta dependence

  static const std::string fmodel_c[nModels] = {"QGSJetII-03", "EPOS-1.99", "QGSJetII-04", "EPOS-LHC"};
  static const std::string fprim_c[2] = {"p", "Fe"};

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

  static const double OffsetM_Mu_Photon[2] = { 0., 0. };

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


  class UnivCalibConstants {

  private:

    int CalibOpt;       // Calibration constants  [-3] Calib Data  [-2] Data [-1] Photon [0] QGSJetII-03  [1] EPOS1.99  [2] QGSJetII-04 [3] EPOS-LHC
    //      Calib Data refers to the analysis of the Golden Data set that is used to produce the calibration constants.


    //---------------------
    // Appliable to MC
    //---------------------
    int HadronicModel;                  //  [0] QGSjetII-03 [1] EPOS1.99 [2] QGSJetII-04 [3] EPOS-LHC
    int RecMixture;                     // Mixtures assumed  (2 component p/Fe)
    //            [0] fp=1.0   [1] fp=0.0  [2] fp=0.5 [3] fp=0.1 [4] fp=0.9  [5] fp changing with energy
    double ProtonFraction[2];           // Energy dependent proton fraction

    //---------------------
    // Appliable to Data
    //---------------------
    //  [0] Standard
    //  [1] FD energy +15%
    //  [2] FD energy -15%
    //  [3] FD Xmax             +10 [g/cm2]
    //  [4] FD Xmax             -10 [g/cm2]
    //  [5] <Xmax>_{f_p=0.5}    +15 [g/cm2]
    //  [6] <Xmax>_{f_p=0.5}    -15 [g/cm2]
    //  [7] ER_{f_p=0.5}        +10 [g/cm2]
    //  [8] ER_{f_p=0.5}        -10 [g/cm2]
    //  [9]-[12]  Use only Eye 1/2/3/4

    int RecSys;

  public:

    UnivCalibConstants(int CalibOpt_i, int RecSys_i,  int RecMixture_i);
    ~UnivCalibConstants() {}


    double GetProtonFraction(double logE);

    // <Xmax>/<Nmu>
    double GetMeanXmax(double logE);
    double GetMeanNmu(double logE, double theta);
    double GetMeanNmu(double logE, double theta, double Xmax);

    // OffsetM_Mu
    double GetOffsetM_Mu(double theta);
    double GetOffsetM_Mu(double logE, double theta, double Xmax);

    // Assumed  <Xmax>(logE) for a 50% p/Fe mixture
    double GetMeanXmax_pFe(double logE);

    //Xmax_SD / Xmax_FD / XmaxMu_Golden / Nmu_FD resolution
    double GetSigmaXmax_SD(double logE, int ndev, bool IsInfill);
    double GetSigmaXmax_SD(double logE, int ndev);
    double GetSigmaXmax_FD(double logE);
    double GetSigmaNmu_FD(double logE);
    double GetTrueRMSXmax(double rms, double logE, bool IsInfill);

  };
}
#endif // __UnivCalibConstants_h_