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

#include <algorithm>
#include <TRandom3.h>
#include <TF2.h>
#include <TH2D.h>

const double W2Min = 4.0;
const double MpSq=0.9383*0.9383;
const int nxB=3, nQ2=4;
const double QSqBinLow[nQ2+1]={1.5, 2.5, 4.0, 6.0, 10.0};
const int nsimu = 10000;
char fName[32];

double xBjMin(double Q2);
double xBjMax(double Q2);
Double_t Likelihood(Double_t *val,Double_t *par);

int main(){
    TRandom3 ran3;
    TCanvas *c1= new TCanvas();
    TFile *input =new TFile("xBQ2binning_simu.root","read");
   // TF2 * fLikeli = new TF2("fLikeli",Likelihood,0.0,1.0,0.0,10.0,4);
   // fLikeli->SetParameters(-0.5,2.0,2.0,1.0);
    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,4);
            fLikeli[jQ2Bin][ixBin]->SetParameters(-0.5,2.0,2.0,1.0);
        }
    }
    TH1F *hist=new TH1F("hist","Histogram",100, 0.0,10);
    
    double xBj,QSq;
    int ixB,jQ2;
    
    double zeta_min=0.0, zeta_max=1.0;
    double zeta;
    double xBMinBin, xBMaxBin;
    double Norm[nQ2][nxB], PhSp[nQ2][nxB];
    
    for (int jQ2Bin=0; jQ2Bin<nQ2; jQ2Bin++) {
        for (int ixBin=0; ixBin<nxB; ixBin++) {
            Norm[jQ2Bin][ixBin] = 0.0;
            PhSp[jQ2Bin][ixBin]= 0.0;
            for (int iSimu=0; iSimu<nsimu; iSimu++){
                QSq = QSqBinLow[jQ2Bin]+ (QSqBinLow[jQ2Bin+1]-QSqBinLow[jQ2Bin])*ran3.Uniform();
                xBMinBin = (xBjMin(QSq)*(nxB-ixBin)+ xBjMax(QSq)*ixBin )/(double)nxB;
                xBMaxBin = (xBjMin(QSq)*(nxB-ixBin-1)+ xBjMax(QSq)*(ixBin+1) )/(double)nxB;
                xBj = xBMinBin+(xBMaxBin-xBMinBin)*ran3.Uniform();
                Norm[jQ2Bin][ixBin] += (fLikeli[jQ2Bin][ixBin]->Eval(xBj,QSq))*(xBMaxBin-xBMinBin);
                PhSp[jQ2Bin][ixBin] += (xBMaxBin-xBMinBin);
            }
            Norm[jQ2Bin][ixBin] *= (QSqBinLow[jQ2Bin+1]-QSqBinLow[jQ2Bin])/nsimu;
            PhSp[jQ2Bin][ixBin] *= (QSqBinLow[jQ2Bin+1]-QSqBinLow[jQ2Bin]) / nsimu;
        }
    }
    
    TTree * xBQ2binning_tree= (TTree*) input->Get("xBQ2binning_tree;1");
    
    xBQ2binning_tree->SetBranchAddress("xBj", &xBj);
    xBQ2binning_tree->SetBranchAddress("QSq", &QSq);
    xBQ2binning_tree->SetBranchAddress("ixB", &ixB);
    xBQ2binning_tree->SetBranchAddress("jQ2", &jQ2);
    
    zeta=(zeta_max-zeta_min)*(ran3.Uniform())+zeta_min;
    
    
    
    double Q2Min[nQ2][nxB],Q2Max[nQ2][nxB],Q2;
    double xMin0,xMax0,xMin[nQ2][nxB], xMax[nQ2][nxB];
    double xBrj;
    int nevents=0;
    
    //  const int nPar=10;
    //Double_t alpha[nPar] = {-1.0,-0.8,-0.6,-0.4,-0.2,0., 0.2,0.4,0.6,0.8};
    //Double_t beta[nPar] = {1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.7,3.0};
    
    TH2D * h2_Q2_vs_xB = new TH2D("h2_Q2_vs_xB","Q^{2} vs x_{B} ; x_{B}  ; Q^{2} (GeV^{2})  ",
                                  50,0.0,1.0, 50,0.0,10.0);
    double LogLi[nQ2][nxB], NBin[jQ2][ixB];
    
    int entries= xBQ2binning_tree->GetEntries();
    for(int ientry=0; ientry < entries;++ientry)
    {
        
        xBQ2binning_tree->GetEntry(ientry);
        nevents++;
        xMin0=(QSq*0.75)/12.;
        xMax0= QSq/(QSq+W2Min-MpSq);
        // printf("XMin0=%f \n", xMin0);
        // printf("jQ2=%d \n", jQ2);
        
        h2_Q2_vs_xB->Fill(xBj,QSq);
        if ( ixB>=0&&ixB<nxB&&jQ2>=0&&jQ2<nQ2){
            LogLi[jQ2][ixB] += log(fLikeli[jQ2][ixB]->Eval(xBj,QSq))-log(Norm[jQ2][ixB]/PhSp[jQ2][ixB]);
            NBin[jQ2][ixB] += 1.0;
        }
        
        
    }
    
    // input->Close();
    printf("LogLikelood values \n ix = ");
    for (int ixBin=0; ixBin<nxB; ixBin++) {
        printf("    %d   ",ixBin);
    }
    for (int jQ2Bin=0; jQ2Bin<nQ2; jQ2Bin++) {
        printf("\n jQ2=%2d ", jQ2Bin);
        for (int ixBin=0; ixBin<nxB; ixBin++) {
            LogLi[jQ2Bin][ixBin] *= 1.0/NBin[jQ2Bin][ixBin];
            printf("%10.3f", LogLi[jQ2Bin][ixBin] );
        }
    }
    printf("\n");
    
    h2_Q2_vs_xB->Draw("cont1z");
    
    
}

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;
}

Double_t Likelihood(Double_t *val,Double_t *par){
    
    Float_t x = val[0];
    Float_t y = val[1];
    Double_t result=0;
    Double_t func=0;
    if (x>0.0) {
        func= TMath::Power(x,par[0]) * TMath::Power((1-x),par[1])*TMath::Power((1+ (y*y)/ par[3]),-par[2]);
        result+=func;
    }
    return result;
    
}