#include <memory>
#include <iostream>
#include <string>

#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TFile.h"
#include "TTree.h"
#include "TLatex.h"
#include "TMath.h"

// user include files
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooPolynomial.h"
#include "RooDataHist.h"
#include "RooAbsPdf.h"
#include "RooAddPdf.h"
#include "RooArgSet.h"
#include "RooArgList.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooNLLVar.h"
#include "RooChi2Var.h"
#include "RooMinuit.h"

using std::cout;
using std::endl;
using std::string;
using namespace RooFit;


void FitMassPeakRooFit( TH1F* h, double xlo, double xhi,  int HistoIndex, int ngaus, bool isEB, bool is_2011_);


void FitPi0( ){

    TFile* inputEpsilonFile_ = TFile::Open( "2012D_epsilonPlots.root" );
    if( !inputEpsilonFile_ ) cout<<"WARNING: File 2012D_epsilonPlots.root do not exist"<<endl;
    bool Are_pi0_ = true;

    for(int iR=0; iR < 1999; iR++){
	char line[100];
	sprintf(line,"Barrel/epsilon_EB_iR_%d",iR);
	TH1F *h1 = (TH1F*)inputEpsilonFile_->Get(line);
	FitMassPeakRooFit( h1, Are_pi0_? 0.08:0.4, Are_pi0_? 0.21:0.65, iR, 1, true, true); //0.05-0.3
    }
}


void FitMassPeakRooFit(TH1F* h, double xlo, double xhi,  int HistoIndex, int ngaus, bool isEB, bool is_2011_) 
{
    //-----------------------------------------------------------------------------------
    bool Are_pi0_ = true;
    RooRealVar x("x","#gamma#gamma invariant mass",xlo, xhi, "GeV/c^2");
    //x.setRange(xlo,xhi);

    RooDataHist dh("dh","#gamma#gamma invariant mass",RooArgList(x),h);

    RooRealVar mean("mean","#pi^{0} peak position", Are_pi0_? 0.116:0.57,  Are_pi0_? 0.105:0.5, Are_pi0_? 0.150:0.62,"GeV/c^{2}");
    RooRealVar sigma("sigma","#pi^{0} core #sigma",0.013, 0.005,0.020,"GeV/c^{2}");


    if(!isEB)  {
	  mean.setRange( Are_pi0_? 0.10:0.45, Are_pi0_? 0.140:0.62); // 0.200
	  mean.setVal(Are_pi0_? 0.120:0.55);
	  sigma.setRange(0.005, 0.060);
    }

    RooRealVar Nsig("Nsig","#pi^{0} yield",1000.,0.,1.e7);
    Nsig.setVal( h->GetSum()*0.1);

    RooGaussian gaus("gaus","Core Gaussian",x, mean,sigma);

    RooRealVar sigmaTail("sigmaTail","#pi^{0} tail #sigma",0.040, 0.020,0.065,"GeV/c^{2}");

    RooGaussian gaus2("gaus2","Tail Gaussian",x, mean,sigmaTail);

    RooRealVar fcore("fcore","f_{core}",0.9,0.,1.);
    RooAddPdf  signal("signal","signal model",RooArgList(gaus,gaus2),fcore);

    RooRealVar p0("p0","p0", 1000.,-1.e5,1.e5);
    RooRealVar p1("p1","p1", -3000.,-1.e5,1.e5);
    RooRealVar p2("p2","p2", 10000.,-1.e5,1.e5);
    RooRealVar p3("p3","p3", -10000.,-1.e5,1.e5);
    RooRealVar p4("p4","p4",-4000.,-1.e5,1.e5);
    RooRealVar p5("p5","p5", 5.,-1.e5,1.e5);
    RooRealVar p6("p6","p6", 6.,-1.e5,1.e5);

    RooRealVar cb0("cb0","cb0", 0.2, -1.,1.);
    RooRealVar cb1("cb1","cb1",-0.1, -1.,1.);
    RooRealVar cb2("cb2","cb2", 0.1, -1.,1.);
    RooRealVar cb3("cb3","cb3",-0.1, -0.5,0.5);
    RooRealVar cb4("cb4","cb4", 0.1, -1.,1.);
    RooRealVar cb5("cb5","cb5", 0.1, -1.,1.);
    RooRealVar cb6("cb6","cb6", 0.3, -1.,1.);

    RooArgList cbpars(cb0,cb1,cb2);
    cbpars.add( cb3);

    RooChebychev bkg("bkg","bkg model", x, cbpars );

    RooRealVar Nbkg("Nbkg","background yield",1.e3,0.,1.e8);
    Nbkg.setVal( h->GetSum()*0.8 );

    RooAbsPdf* model=0;

    RooAddPdf model1("model","sig+bkg",RooArgList(gaus,bkg),RooArgList(Nsig,Nbkg));
    RooAddPdf model2("model","sig+bkg",RooArgList(signal,bkg),RooArgList(Nsig,Nbkg));

    if(ngaus==1)      model = &model1;
    else if(ngaus==2) model = &model2;


    RooNLLVar nll("nll","log likelihood var",*model,dh, Extended());
    RooMinuit m(nll);
    m.setVerbose(kFALSE);
    //m.setVerbose(kTRUE);
    m.migrad();
    //m.hesse();
    RooFitResult* res = m.save() ;

    RooChi2Var chi2("chi2","chi2 var",*model,dh, true);
    int ndof = h->GetNbinsX() - res->floatParsFinal().getSize();

    // compute S/B

    //compute S/B and chi2
    x.setRange("sobRange",mean.getVal()-3.*sigma.getVal(), mean.getVal()+3.*sigma.getVal());
    RooAbsReal* integralSig = gaus.createIntegral(x,NormSet(x),Range("sobRange"));

    RooAbsReal* integralBkg = bkg.createIntegral(x,NormSet(x),Range("sobRange"));

    float normSig = integralSig->getVal();
    float normBkg = integralBkg->getVal();

    RooPlot*  xframe = x.frame(h->GetNbinsX());
    xframe->SetTitle(h->GetTitle());
    dh.plotOn(xframe);
    model->plotOn(xframe,Components(bkg),LineStyle(kDashed), LineColor(kRed));
    model->plotOn(xframe);

    xframe->Draw();

    TLatex lat;
    char line[300];
    lat.SetNDC();
    lat.SetTextSize(0.040);
    lat.SetTextColor(1);

    float xmin(0.55), yhi(0.80), ypass(0.05);
    sprintf(line,"Yield: %.0f #pm %.0f", Nsig.getVal(), Nsig.getError() );
    lat.DrawLatex(xmin,yhi, line);

    sprintf(line,"m_{#gamma#gamma}: %.2f #pm %.2f", mean.getVal()*1000., mean.getError()*1000. );
    lat.DrawLatex(xmin,yhi-ypass, line);

    sprintf(line,"#sigma: %.2f #pm %.2f (%.2f%s)", sigma.getVal()*1000., sigma.getError()*1000., sigma.getVal()*100./mean.getVal(), "%" );
    lat.DrawLatex(xmin,yhi-2.*ypass, line);

	  std::stringstream ind;
	  ind << (int) HistoIndex;
	  TString nameHistofit = "Fit_n_" + ind.str();
	  xframe->SetName(nameHistofit.Data());
	  xframe->Write();
}