#ifndef PROF_H
#define PROF_H

#include <iostream>
#include <string>
#include <vector>
#include <TH1D.h>
#include <TF1.h>
#include <TRandom3.h>
#include <TMinuit.h>
#include <Math/Random.h>
#include <Math/GSLRndmEngines.h>

using namespace std;

class Prof {

  public: 

    Prof (int N, int M); 
    Prof (string inputFile);
    TMinuit* minuit() { return m_minuit; }

    void getStartPar();
    double lnL(vector<double> parVec, int M);
    double q();
    double q(int M, int N);
    int fit(int M);
    int N(){ return m_N; }
    int M(){ return m_M; }
    void setN(int N) { m_N = N; m_npar = m_M + 2*m_N; }     // num measurements
    void setM(int M) { m_M = M; m_npar = m_M + 2*m_N; }     // num free poi
    int npar(){ return m_npar; }
    void generateData(int N, int M, vector<double> parVec,
		      TRandom3* ran, 
		      ROOT::Math::Random<ROOT::Math::GSLRngMT>* rgsl);
    vector<double> mu(vector<double> theta);
    double y(int i){ return m_yVec[i]; }
    double u(int i){ return m_uVec[i]; }
    double s(int i){ return m_sVec[i]; }
    double r(int i){ return m_rVec[i]; }
    double sigma_y(int i){ return m_sigma_y_Vec[i]; }
    void set_y(double y, int i){ m_yVec[i] = y; }
    void set_u(double u, int i){ m_uVec[i] = u; }
    void set_s(double s, int i){ m_sVec[i] = s; }
    void set_r(double r, int i){ m_rVec[i] = r; }
    void set_sigma_y(double sigma_y, int i){ m_sigma_y_Vec[i] = sigma_y; }
    double startPar(int i){ return m_startPar[i]; }
    double stepSize(int i){ return m_stepSize[i]; }
    void setRanGen(TRandom3* ran, 
		   ROOT::Math::Random<ROOT::Math::GSLRngMT>* rgsl){
      m_ran = ran;
      m_rgsl = rgsl;
    }
    void setMinosErrors(double eparab, double eplus, double eminus){
      m_eparab = eparab;
      m_eplus = eplus;
      m_eminus = eminus;      // note this is negative
    }
    double getMuHat() { return m_muHat; }
    double getSigmaMuHat(){ return m_sigmaMuHat; }
    double getEparab(){ return m_eparab; }
    double getEplus(){ return m_eplus; }
    double getEminus(){ return m_eminus; }

    vector<double> parVec(){  return m_parVec; }
    vector<double> sigParVec(){ return m_sigParVec; }
    void setParVec(vector<double> parVec) { m_parVec = parVec; }
    void setSigParVec(vector<double> sigParVec) { m_sigParVec = sigParVec; }

  private:

    int                  m_N;        // number of input measurements
    int                  m_M;        // number of param. of interest
    int                  m_npar;

    double               m_muHat;
    double               m_sigmaMuHat;
    double               m_eparab;
    double               m_eplus;
    double               m_eminus;

    vector<double>       m_parVec;
    vector<double>       m_sigParVec;
    vector<double>       m_yVec;
    vector<double>       m_sigma_y_Vec;
    vector<double>       m_sVec;
    vector<double>       m_rVec;
    vector<double>       m_uVec;
    vector<double>       m_startPar;
    vector<double>       m_stepSize;

    TMinuit* m_minuit;
    TRandom3* m_ran;
    ROOT::Math::Random<ROOT::Math::GSLRngMT>* m_rgsl;

};

#endif