#include <iostream>
#include <stdlib.h>
#include <sstream>
#include <fstream>
#include <vector>
#include <string>
#include <unistd.h>

#include <RooRandom.h>
#include <RooMsgService.h>
#include <RooRealVar.h>
#include <RooArgList.h>
#include <RooBernstein.h>
#include <RooGenericPdf.h>
#include <RooGaussian.h>
#include <RooDataHist.h>
#include <RooAddPdf.h>
#include <RooPlot.h>
#include <RooNLLVar.h>
#include <RooChi2Var.h>
#include <RooMinimizer.h>
#include <RooDataHist.h>
#include <RooHistPdf.h>

#include <TApplication.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TPad.h>
#include <TTree.h>
#include <TChain.h>

#include "generateSimulation.h"

using namespace std;
using namespace RooFit;

TApplication Runner("gui",0,NULL);

string bkg_generating="B",bkg_fitting="B";
int num_simulations=0;



RooRealVar *x;

RooRealVar *a0_gen;
RooRealVar *a1_gen;
RooRealVar *a2_gen;
RooRealVar *a3_gen;
RooRealVar *a4_gen;

RooRealVar *e0_gen;
RooRealVar *e1_gen;
RooRealVar *e2_gen;


RooRealVar *nb;
RooRealVar *nZ;


RooRealVar *a0_fit;
RooRealVar *a1_fit;
RooRealVar *a2_fit;
RooRealVar *a3_fit;
RooRealVar *a4_fit;

RooRealVar *e0_fit;
RooRealVar *e1_fit;
RooRealVar *e2_fit;

RooRealVar *mean_fit;
RooRealVar *sigma_fit;


int n_bkg=250000,n_signal=800;


void inizializeParameters();


int main(int narg,char* argv[]){
    
    RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
    RooRandom::randomGenerator()->SetSeed(0);
 
    if(narg!=3){
        cout<<"USE: ./Executable --nsim=zzz --sim=XY"<<endl;
        cout<<"where    zzz is the number of toys you want to generate"<<endl;
        cout<<"         X(= E or B) is the bkg-model for the TOYMC generation"<<endl;
        cout<<"         X(= E or B) is the bkg-model for the TOYMC interpolation"<<endl;
        cout<<"             E=Exponential"<<endl;
        cout<<"             B=Bernstein"<<endl<<endl;
        exit(0);
    }
    
    for(int i=1;i<narg;i++){
        string tmp(argv[i]);
        
        if (tmp.find("--nsim=")==string::npos & tmp.find("--sim=")==string::npos){cout<<"No valid commands specified"<<endl<<endl;
            cout<<"USE: ./Executable --nsim=zzz --sim=XY"<<endl;
            cout<<"where    zzz is the number of toys you want to generate"<<endl;
            cout<<"         X(= E or B) is the bkg-model for the TOYMC generation"<<endl;
            cout<<"         X(= E or B) is the bkg-model for the TOYMC interpolation"<<endl;
            cout<<"             E=Exponential"<<endl;
            cout<<"             B=Bernstein"<<endl<<endl;
            exit(0);}
    }
    
    for(int i=1;i<narg;i++){
        string tmp(argv[i]);
        
        if (tmp.find("--nsim=")!=string::npos){
            istringstream is(tmp.substr(7,tmp.length()-7));
            is>>num_simulations;
        }
        
        if (tmp.find("--sim=")!=string::npos){
            istringstream is_gen(tmp.substr(6,1));
            is_gen>>bkg_generating;
            
            istringstream is_fit(tmp.substr(7,tmp.length()-7));
            is_fit>>bkg_fitting;
        }
    }
    
    
    
    string strPath1="./Ztry/",strPath0="./Ztry";
    strPath1+=bkg_generating;
    strPath1+=bkg_fitting;
    
    if(access( strPath0.c_str(), 0 )!=0){
        ostringstream tmp_com0;
        tmp_com0<<"mkdir ./Ztry";
        system(tmp_com0.str().c_str());
    }
    
    if(access( strPath1.c_str(), 0 )!=0){
        ostringstream tmp_com1,tmp_com2,tmp_com3;
        
        tmp_com1<<"mkdir ./Ztry/"<<bkg_generating<<bkg_fitting;
        tmp_com2<<"mkdir ./Ztry/"<<bkg_generating<<bkg_fitting<<"/plots";
        tmp_com3<<"mkdir ./Ztry/"<<bkg_generating<<bkg_fitting<<"/datas";
        
        system(tmp_com1.str().c_str());
        system(tmp_com2.str().c_str());
        system(tmp_com3.str().c_str());
    }

    
    
    
    x=new RooRealVar("x","x",60,240);
    
    // Functions
    
    a0_gen=new RooRealVar("a0_gen","a0_gen",27.1794,0.,500.);
    a1_gen=new RooRealVar("a1_gen","a1_gen",10.9069,0.,500.);
    a2_gen=new RooRealVar("a2_gen","a2_gen",6.83022,0.,500.);
    a3_gen=new RooRealVar("a3_gen","a3_gen",1.96044,0.,500.);
    a4_gen=new RooRealVar("a4_gen","a4_gen",1.34874,0.,500.);
    
    RooBernstein Bern_gen("Bern_gen","Bern_gen",*x,RooArgList(*a0_gen,*a1_gen,*a2_gen,*a3_gen,*a4_gen));
    
    e0_gen=new RooRealVar("e0_gen","e0_gen",1.12188,0.,500.);
    e1_gen=new RooRealVar("e1_gen","e1_gen",112.228,0.,500.);
    e2_gen=new RooRealVar("e2_gen","e2_gen",117.536,0.,500.);
    
    RooGenericPdf Expo_gen("Expo_gen","e0_gen*exp(-x/e1_gen)-exp(-x/e2_gen)",RooArgList(*x,*e0_gen,*e1_gen,*e2_gen));
    
    
    a0_fit=new RooRealVar("a0_fit","a0_fit",27.1794,0.,500.);
    a1_fit=new RooRealVar("a1_fit","a1_fit",10.9069,0.,500.);
    a2_fit=new RooRealVar("a2_fit","a2_fit",6.83022,0.,500.);
    a3_fit=new RooRealVar("a3_fit","a3_fit",1.96044,0.,500.);
    a4_fit=new RooRealVar("a4_fit","a4_fit",1.34874,0.,500.);
    
    RooBernstein Bern_fit("Bern_fit","Bern_fit",*x,RooArgList(*a0_fit,*a1_fit,*a2_fit,*a3_fit,*a4_fit));
    
    e0_fit=new RooRealVar("e0_fit","e0_fit",1.12188,0.,500.);
    e1_fit=new RooRealVar("e1_fit","e1_fit",112.228,0.,500.);
    e2_fit=new RooRealVar("e2_fit","e2_fit",117.536,0.,500.);
    
    RooGenericPdf Expo_fit("Expo_fit","e0_fit*exp(-x/e1_fit)-exp(-x/e2_fit)",RooArgList(*x,*e0_fit,*e1_fit,*e2_fit));
    
    
    
    
    TFile *f = new TFile("Ztemplate.root");
    TH1D *myz = (TH1D*)f->Get("Zmass5");
    
    RooDataHist h_Z("h_Z","h_Z",*x,myz);
    RooHistPdf pdf_Z("pdf_Z","pdf_Z",RooArgList(*x),h_Z,2);
    
    
    nb=new RooRealVar("nb","nb",n_bkg,0.,5*n_bkg);
    nZ=new RooRealVar("nZ","nZ",n_signal,-100*n_signal,100*n_signal);
    
    
    
    RooAddPdf *fit_pdf;
    if(bkg_fitting=="B")fit_pdf=new RooAddPdf("fit_pdf","fit_pdf",RooArgList(Bern_fit,pdf_Z),RooArgList(*nb,*nZ));
    if(bkg_fitting=="E")fit_pdf=new RooAddPdf("fit_pdf","fit_pdf",RooArgList(Expo_fit,pdf_Z),RooArgList(*nb,*nZ));
    
    vector<double> results_chi2,results_nll;
    vector<double> value_chi2,err_chi2,value_nll,err_nll;
    
    for(int i=0;i<num_simulations;i++){
        
        cout <<"\x1b[A";
        printf("%c[2K", 27);
        cout<<"\rSimulation number: "<<i+1<<"/"<<num_simulations<<"     ["<<(double)(i+1)*100/num_simulations <<" %]"<<endl;
        
    simulating_part:
    
        TH1D *simulation;
        if(bkg_generating=="B")simulation=(TH1D*)generateSimulation(x,&Bern_gen,&pdf_Z,n_bkg,n_signal);
        if(bkg_generating=="E")simulation=(TH1D*)generateSimulation(x,&Expo_gen,&pdf_Z,n_bkg,n_signal);
    
        RooDataHist simulatedData("simulatedData","simulatedData",*x,simulation);
        
        inizializeParameters();
        
        RooChi2Var *Chi2=new RooChi2Var("Chi2","Chi2",*fit_pdf,simulatedData,Extended(true));
        RooMinimizer *minuit_chi2=new RooMinimizer(*Chi2);
        minuit_chi2->setPrintLevel(-1);
        int status_migrad=minuit_chi2->migrad();
        int status_hesse=minuit_chi2->hesse();
        int status_minos=minuit_chi2->minos();
        
        results_chi2.push_back((nZ->getValV()-n_signal)/nZ->getError());
        value_chi2.push_back(nZ->getValV());
        err_chi2.push_back(nZ->getError());
        
        inizializeParameters();
        
        RooNLLVar *NLL=new RooNLLVar("NLL","NLL",*fit_pdf,simulatedData,Extended(true));
        RooMinimizer *minuit_nll=new RooMinimizer(*NLL);
        minuit_nll->setPrintLevel(-1);
        status_migrad=minuit_nll->migrad();
        status_hesse=minuit_nll->hesse();
        status_minos=minuit_nll->minos();
        
        results_nll.push_back((nZ->getValV()-n_signal)/nZ->getError());
        value_nll.push_back(nZ->getValV());
        err_nll.push_back(nZ->getError());
        
        simulation->Reset("ICES");
    }
    
    
    ofstream file_chi2,file_nll;
    
    string name_chi2="./Ztry/",name_nll="./Ztry/";
    
    name_chi2+=bkg_generating;
    name_chi2+=bkg_fitting;
    name_chi2+="/datas/datas_";
    name_chi2+=bkg_generating;
    name_chi2+=bkg_fitting;
    name_chi2+="_Chi2.txt";

    name_nll+=bkg_generating;
    name_nll+=bkg_fitting;
    name_nll+="/datas/datas_";
    name_nll+=bkg_generating;
    name_nll+=bkg_fitting;
    name_nll+="_NLL.txt";
    
    file_chi2.open(name_chi2.c_str());
    file_nll.open(name_nll.c_str());
    
    for(int i=0;i<value_chi2.size();i++){
        file_chi2<<value_chi2[i]<<" "<<err_chi2[i]<<" "<<results_chi2[i]<<endl;
        file_nll<<value_nll[i]<<" "<<err_nll[i]<<" "<<results_nll[i]<<endl;
    }
    
    file_chi2.close();
    file_nll.close();
    
    
    
    TH1D *histo_chi2=new TH1D("histo_chi2","histo_chi2",60,-6,6);
    TH1D *histo_nll=(TH1D*)histo_chi2->Clone();
    
    for(int i=0;i<results_chi2.size();i++){
        histo_chi2->Fill(results_chi2[i]);
        histo_nll->Fill(results_nll[i]);
    }
    
    RooRealVar *y=new RooRealVar("y","y",-6,6);
    
    RooDataHist Histo_chi2("Histo_chi2","Histo_chi2",*y,histo_chi2);
    RooDataHist Histo_nll("Histo_nll","Histo_nll",*y,histo_nll);
    
    RooRealVar mean_chi2("mean_chi2","mean_chi2",0,-6,6);
    RooRealVar sigma_chi2("sigma_chi2","sigma_chi2",1,0.,6.);
    RooGaussian gauss_chi2("gauss_chi2","gauss_chi2",*y,mean_chi2,sigma_chi2);
    
    RooRealVar mean_nll("mean_nll","mean_nll",0,-6,6);
    RooRealVar sigma_nll("sigma_nll","sigma_nll",1,0.,6.);
    RooGaussian gauss_nll("gauss_nll","gauss_nll",*y,mean_nll,sigma_nll);
    
    gauss_chi2.fitTo(Histo_chi2,Hesse(true),Minos(true),Extended(true));
    gauss_nll.fitTo(Histo_nll,Hesse(true),Minos(true),Extended(true));
    
    
    TCanvas *c0=new TCanvas;
    
    RooPlot *final_frame=y->frame();
    Histo_chi2.plotOn(final_frame,LineColor(2),MarkerColor(2));
    gauss_chi2.plotOn(final_frame,LineColor(2),LineStyle(2));
    
    
    Histo_nll.plotOn(final_frame,LineColor(4),MarkerColor(4));
    gauss_nll.plotOn(final_frame,LineStyle(2));
    final_frame->Draw();
    
    
    gPad->SetGridx();
    gPad->SetGridy();
    
    cout<<endl;
    cout<<"Medie: "<<mean_chi2.getValV()<<" "<<mean_nll.getValV()<<endl<<endl;
    
    c0->Draw();
    
    string name_plot="./Ztry/";
    name_plot+=bkg_generating;
    name_plot+=bkg_fitting;
    name_plot+="/plots/plot_";
    name_plot+=bkg_generating;
    name_plot+=bkg_fitting;
    name_plot+=".pdf";
    
    c0->SaveAs(name_plot.c_str());
    
    
    
    Runner.Run(true);
    
    
}








void inizializeParameters(){
    
    a0_fit->setVal(a0_gen->getValV());
    a1_fit->setVal(a1_gen->getValV());
    a2_fit->setVal(a2_gen->getValV());
    a3_fit->setVal(a3_gen->getValV());
    a4_fit->setVal(a4_gen->getValV());
    
    e0_fit->setVal(e0_gen->getValV());
    e1_fit->setVal(e1_gen->getValV());
    e2_fit->setVal(e2_gen->getValV());
    
    nb->setVal(n_bkg);
    nZ->setVal(n_signal);
    
}