//
//  Tree_reading.cpp
//
//
//  Created by Dilini Bulumulla on 6/3/22.
//


#include <algorithm>
#include <TRandom3.h>
#include <TF2.h>
#include <TH2D.h>
#include <TMath.h>
#include <TFile.h>
#include <TCanvas.h>
//#include "TMinuit.h"
using namespace ROOT::Math;
const double W2Min = 4.0;
const double MpSq=0.9383*0.9383;
const double kBeam=10.6;
const double QSqMin = 1.5;
const double QSqMax = 10.0;
const int nsimu = 10000;

double xBjMin(double Q2);
double xBjMax(double Q2);
//Double_t Likelihood(Double_t *val,Double_t *par);
//void myFCN(int &npar, double *gin, double &logli, double *par, int iflag);


Double_t Likelihood(Double_t *val,Double_t *par){
    double kBeam=10.6;
    const double MpSq=0.9383*0.9383;
    Float_t x = val[0];
    Float_t y = val[1];
    Double_t result=0;
    Double_t func=0;
    if ((y<QSqMin)||(y>QSqMax)) return result;
    if ((x<xBjMin(y))||(x>xBjMax(y))) return result;
    if (x>0.0) {
        func= par[4]*TMath::Power(x,par[0]) * TMath::Power((1-x),par[1])*TMath::Power((1+ (y*y)/ par[3]),-par[2]);
        Double_t y_val   = y/(x*2.0*kBeam*sqrt(MpSq));
        func*= (1-y_val+y_val*y_val/2.)/(x*y*y);
        result = func;
    }
    return result;
    
}

/*double func(const double *x){
    
    const int nxB=3, nQ2=4;
    const double QSqBinLow[nQ2+1]={1.5,2.5,4.0 ,6.0 ,10.0};
    const int nQ2Log =16;
    double Q2BinLowLog[nQ2Log+1];
    double fQ2 = pow(QSqBinLow[nQ2]/QSqBinLow[0],1./((double)nQ2Log));
    double Q2val = QSqBinLow[0];
    TF2 *f2 = new TF2("f2",Likelihood,0.0,1.0,0.0,10.0,5);
    for (int iQ2=0; iQ2<=nQ2Log; iQ2++)
    {
        Q2BinLowLog[iQ2]=Q2val;
        Q2val *= fQ2;
    }
    
    return  f2->Integral(0,1,1.,);
}*/

void Fitting(){
    TRandom3 ran3;
    //   gStyle->SetPalette(kRainBow);
    //open 2D raw and simu histograms
    // TFile *hFile=TFile::Open("xBQ2binning_uniform_generator_new_bin.root");
    TFile *hFile=TFile::Open("xBQ2binning_SIDIS_simu_sumw2.root");
    TH2F *h_QSq_xBj_log_bin;
    TH2F *h_QSq_xBj_log_bin_MC;
    TH2F *h_ratio;
    hFile->GetObject("h_QSq_xBj_log_bin", h_QSq_xBj_log_bin);
    hFile->GetObject("h_QSq_xBj_log_bin_MC", h_QSq_xBj_log_bin_MC);
    
    h_ratio = (TH2F*)h_QSq_xBj_log_bin_MC->Clone();
    h_ratio->GetXaxis()->SetTitle(" ");
    h_ratio->GetYaxis()->SetTitle(" ");
    h_ratio->SetTitle("h_QSq_xBj_log_bin_MC(raw)/h_QSq_xBj_log_bin(simu)");
    h_ratio->Divide(h_QSq_xBj_log_bin);
    
    const int nxB=3, nQ2=4;
    const double QSqBinLow[nQ2+1]={1.5,2.5,4.0 ,6.0 ,10.0};
    const int nQ2Log =16;
    double Q2BinLowLog[nQ2Log+1];
    double fQ2 = pow(QSqBinLow[nQ2]/QSqBinLow[0],1./((double)nQ2Log));
    double Q2val = QSqBinLow[0];
    for (int iQ2=0; iQ2<=nQ2Log; iQ2++)
    {
        Q2BinLowLog[iQ2]=Q2val;
        Q2val *= fQ2;
    }
    
    TH2F* h_QSq_xBj_log_bin_data=new TH2F("h_QSq_xBj_log_bin_data", "QSq vs xBj log binning cross section", 50,0.0,1.0,nQ2Log,Q2BinLowLog);
    h_QSq_xBj_log_bin_data->Sumw2();
    TH2F* h_err_QSq_xBj_log_bin_data=new TH2F("h_err_QSq_xBj_log_bin_data", "QSq vs xBj log binning cross section error", 50,0.0,1.0,nQ2Log,Q2BinLowLog);
    
    TH2F *h_ratio_errors = (TH2F*)h_ratio->Clone("h_ratio_errors");
    for (int i = 0; i < h_ratio_errors->GetSize(); i++)
        //  h_ratio_errors->GetArray()[i] = h_ratio_errors->GetBinError(i);
        h_ratio_errors->GetArray()[i] = h_ratio_errors->GetBinError(i)/(h_ratio_errors->GetBinContent(i)+0.1);
    h_ratio_errors->Sumw2(kFALSE);
    h_ratio_errors->ResetStats();
    h_ratio_errors->SetTitle("Errors of the ratio histogram");
    
    TFile *input =new TFile("xBQ2binning_data.root","read");
    
    int nparMINUIT=4;
    TCanvas *c1= new TCanvas();
    
    ROOT::Math::Minimizer* minimum[nQ2][nxB];
    
    char fName[64];
    TH1F *hist=new TH1F("hist","Histogram",100, 0.0,10);
    
    double xBj,QSq;
    int ixB,jQ2;
    double weight_total,Sigma_real,weight_systamatic;
    double Full_error,error_weight,error_systamatic,Rel_error;
    double zeta_min=0.0, zeta_max=1.0;
    double zeta;
    double xBMinBin, xBMaxBin;
    double Norm[nQ2][nxB], PhSp[nQ2][nxB];
    float integrated_luminosity=54.6325; //this is in femto barns.
    
    
    TTree * xBQ2binning_tree= (TTree*) input->Get("xBQ2binning_tree");
    
    
    
    xBQ2binning_tree->SetBranchAddress("xBj", &xBj);
    xBQ2binning_tree->SetBranchAddress("QSq", &QSq);
    xBQ2binning_tree->SetBranchAddress("ixB", &ixB);
    xBQ2binning_tree->SetBranchAddress("jQ2", &jQ2);
    xBQ2binning_tree->SetBranchAddress("weight_total", &weight_total);
    xBQ2binning_tree->SetBranchAddress("Sigma_real",&Sigma_real);
    xBQ2binning_tree->SetBranchAddress("weight_systamatic", &weight_systamatic);
    
    double Q2Min[nQ2][nxB],Q2Max[nQ2][nxB],Q2;
    double xMin0,xMax0,xMin[nQ2][nxB], xMax[nQ2][nxB];
    double Q2sim, xBjsim, ysim, yEvt;
    
    TH1F *hist_1Bin= new TH1F("hist_1Bin", "xBj", 100, 0.0,1.0);
    // TH2D * h2_Q2_vs_xB = new TH2D("h2_Q2_vs_xB","Q^{2} vs x_{B} ; x_{B}  ; Q^{2} (GeV^{2})  ",
    // 60,0.1,0.6, 85,1.5,10.0);
    
    
    
    TF1 *func = new TF1("func",Likelihood,0.0,1.0,5);
    func->SetParameters(-0.5,2.0,2.0,1.0,30000.);
    TF2 * fLikeli[nQ2][nxB];
    for (int jQ2Bin=0; jQ2Bin<nQ2; jQ2Bin++) {
        for (int ixBin=0; ixBin<nxB; ixBin++) {
            sprintf(fName,"fLikeli_Q2_%02d_xB_%02d",jQ2Bin, ixBin);
            fLikeli[jQ2Bin][ixBin]= new TF2(fName,Likelihood,0.0,1.0,0.0,10.0,5);
            fLikeli[jQ2Bin][ixBin]->SetParameters(-0.5,2.0,2.0,1.0,30000.);
        }
    }
    
    TF2 *f2 = new TF2("f2",Likelihood,0.0,1.0,0.0,10.0,5);
    
    
    double LogLi[nQ2][nxB], NBin[nQ2][nxB];
    
    //    Put event loop inside bin loop for minimization purposes
    int entries= xBQ2binning_tree->GetEntries();
    
    
    
    //  if (fChain == 0) return;
    for(int ientry=0; ientry < entries;++ientry)
    {
        
        xBQ2binning_tree->GetEntry(ientry);
     /*   if(ixB==1){
            hist_1Bin->Fill(xBj);
        }*/
        
        
        //  if(ixB==1 && jQ2==2){
        if(weight_total  < std::numeric_limits<float>::min()  || weight_total  > std::numeric_limits<float>::max()) continue;
        if(weight_total>0){
            if(TMath::Finite(weight_total) && !TMath::IsNaN(weight_total)) {
                
                h_QSq_xBj_log_bin_data->Fill(xBj,QSq,weight_total);
                
            }
        }
        if(weight_systamatic > 0.0 && weight_total > 0.0){
            if (TMath::Finite(weight_systamatic)  && !TMath::IsNaN(weight_systamatic) )  {
                h_err_QSq_xBj_log_bin_data->Fill(xBj,QSq,weight_systamatic);
            }
        }
        //  }
    }  //end of event loop
    f2->SetParameter(0,-0.5);
    f2->SetParameter(1,2.0);
    f2->SetParameter(2,0.0);
    f2->FixParameter(3,1.0);
    f2->SetParameter(4,7e6);
    
    
    
    gStyle->SetPaintTextFormat("4.2f");
    TCanvas* C_QSqxBj= new TCanvas();
    C_QSqxBj->Divide(2,2);
    C_QSqxBj->cd(1);
    h_QSq_xBj_log_bin->SetMarkerSize(1.1);
    h_QSq_xBj_log_bin->Draw("COLZ");
    C_QSqxBj->cd(2);
    h_QSq_xBj_log_bin_MC->SetMarkerSize(1.1);
    h_QSq_xBj_log_bin_MC->Draw("COLZ");
    C_QSqxBj->cd(3);
    h_ratio->SetMarkerSize(1.1);
    h_ratio->SetMarkerColor(kRed);
    h_ratio->Draw("COLZ");
    C_QSqxBj->cd(4);
    h_ratio_errors->SetMarkerSize(1.1);
    h_ratio_errors->Draw(" COLZ TEXT45");
    
    
    
    
    //cross section calculation.
    //events*efficiency factor*Acceptance factor/Integrated luminosity
    TCanvas* C_Cross_section=new TCanvas();
    
    //multiply by acceptance factor
    TH2F *h2_Q2_vs_xB_efficiency_factor_acceptance_factor = (TH2F*)h_QSq_xBj_log_bin_data->Clone("h2_Q2_vs_xB_efficiency_factor_acceptance_factor");
    TH2F *h_QSq_xBj_log_bin_data_Stat_Sys_err=new TH2F("h_QSq_xBj_log_bin_data_Stat_Sys_err", "Stat+Sys SetBinContent Relative error", 50,0.0,1.0,nQ2Log,Q2BinLowLog);
    h2_Q2_vs_xB_efficiency_factor_acceptance_factor->Multiply(h_ratio); //h_ratio histogram is from SIDIS exclusive Simulation
    //divide by Iintegrated Luminosity
    TH2F *h2_Q2_vs_xB_efficiency_factor_acceptance_factor_Inte_Lumi = (TH2F*)h2_Q2_vs_xB_efficiency_factor_acceptance_factor->Clone("h2_Q2_vs_xB_efficiency_factor_acceptance_factor_Inte_Lumi");
    h2_Q2_vs_xB_efficiency_factor_acceptance_factor_Inte_Lumi->Scale(1/integrated_luminosity);
    h2_Q2_vs_xB_efficiency_factor_acceptance_factor_Inte_Lumi->Draw("colz");
    
    
    
    TH2F * h_QSq_xBj_log_bin_data_Stat_Sys= (TH2F*)h_QSq_xBj_log_bin_data->Clone("h_QSq_xBj_log_bin_data_Stat_Sys");
    
    for(int ixbin=1;ixbin<= h_QSq_xBj_log_bin_data->GetNbinsX();++ixbin){
        for(int iybin=1;iybin<= h_QSq_xBj_log_bin_data->GetNbinsY();++iybin){
            error_weight = h_QSq_xBj_log_bin_data->GetBinError(ixbin,iybin);
            error_systamatic= h_err_QSq_xBj_log_bin_data->GetBinError(ixbin,iybin);
            
            Full_error= sqrt(error_weight * error_weight  + error_systamatic* error_systamatic);
            Rel_error = h_QSq_xBj_log_bin_data_Stat_Sys->GetBinContent(ixbin,iybin);
            
            if(Rel_error >0.0){
                Rel_error = Full_error/Rel_error;
                h_QSq_xBj_log_bin_data_Stat_Sys_err->SetBinContent(ixbin,iybin,Rel_error);
            } else {
                h_QSq_xBj_log_bin_data_Stat_Sys_err->SetBinContent(ixbin,iybin,0.0);
            }
            
            
            
            if(Full_error >0.0){
                
                
                h_QSq_xBj_log_bin_data_Stat_Sys->SetBinError(ixbin,iybin,Full_error);
             //   h_QSq_xBj_log_bin_data_Stat_Sys_err->SetBinContent(ixbin,iybin,Full_error);
            }
        }
    }
    
    TCanvas *C_sys=new TCanvas();
    h_err_QSq_xBj_log_bin_data->Draw("LEGO E");
    
    TCanvas *C_stat_sys= new TCanvas();
    h_QSq_xBj_log_bin_data_Stat_Sys->Draw("colz");
    
     TCanvas *C_stat_sys_err= new TCanvas();
    h_QSq_xBj_log_bin_data_Stat_Sys_err->Draw("COLZ");
    
    
    TCanvas *C_cross_section_2= new TCanvas();
    h2_Q2_vs_xB_efficiency_factor_acceptance_factor_Inte_Lumi->Draw("LEGO");
    f2->SetLineColor(1);
    f2->DrawClone("lego same");
    f2->SetLineColor(2);
    h2_Q2_vs_xB_efficiency_factor_acceptance_factor_Inte_Lumi->Fit("f2","I","");
    
} //void Tree_reading()



double xBjMin(double Q2){
    double result = Q2*0.75/12.0;
    return result;
}
double xBjMax(double Q2){
    double result = 1.0;
    result = Q2/(W2Min+Q2-MpSq);
    return result;
}