#ifndef EOE_H #define EOE_H // Header file for declaration of Eoe class #include #include #include #include #include #include #include

#include using namespace std;

class Eoe {

  public: 

  Eoe (string inputFileName, TRandom3* ran, 
	 ROOT::Math::Random * rgsl);

    void setStartPar();
    double lnL(vector parVec);
    double lnL(double mu, double b, int i);
    double lnLPrime(vector muVec, vector thetaVec, 
		    vector rVec, vector uVec, 
		    vector vVec);

    int fit();
    int L(){ return m_L; }
    int M(){ return m_M; }
    int N(){ return m_N; }
    int npar(){ return m_npar; }
    void generateData(vector parVec);
    vector mu(vector theta);
    double y(int i){ return m_yVec[i]; }
    double v(int i){ return m_vVec[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; m_vVec[i] = s*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; }
    void set_sigma_u(double sigma_u, int i){ m_sigma_u_Vec[i] = sigma_u; }
    void set_sigma_u_Vec(vector sigma_u_Vec){
      for (int i=0; i parVec);     // for Bartlett corr. and pVal
    double pVal() { return m_p; }    // value set by setEqMC

    void setMinosErrors(double eparab, double eplus, double eminus){
      m_eparab = eparab;
      m_eplus = eplus;
      m_eminus = eminus;      // note this is negative
    }
    double getEparab(){ return m_eparab; }
    double getEplus(){ return m_eplus; }
    double getEminus(){ return m_eminus; }

    void setUp(double up){ m_up = up; }
    double up(){ return m_up; }
    double muFit(vector sigma_u_Vec);
    double phi(int i, double mu);           // special case for average
    double phi(int i, vector muVec);
    vector thetaHatHatLinVec(double mu);
    vector thetaHatHatLinVec(vector muVec);

    vector parVec(){  return m_parVec; }
    vector sigParVec(){ return m_sigParVec; }
    void setParVec(vector parVec) { m_parVec = parVec; }
    void setSigParVec(vector sigParVec) { m_sigParVec = sigParVec; }
    vector sigma_u_HatHatVec(vector thetaVec);
    double R(int i, int j){ return m_R[i][j]; }

  private:

    bool                 m_EqMCSet;
    bool                 m_pValSet;

    int                  m_L;        // number of input measurements
    int                  m_M;        // number of parameters of interest
    int                  m_N;        // number of control measurements
    int                  m_npar;

    double               m_up;
    double               m_eparab;
    double               m_eplus;
    double               m_eminus;
    double               m_EqMC;
    double               m_p;

    vector m_parVec;
    vector m_sigParVec;
    vector m_yVec;
    vector m_sigma_y_Vec;
    vector m_sigma_u_Vec;
    vector m_sVec;
    vector m_vVec;
    vector m_rVec;
    vector m_uVec;
    vector m_startPar;
    vector m_stepSize;
    vector m_sigma_u_hathat_Vec;
    vector m_bHatVec;
    vector >  m_cov_bHat;
    vector >  m_R;

    TRandom3* m_ran;
    ROOT::Math::Random * m_rgsl;

};

#endif