//
//  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 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);

void Tree_reading(){
    TRandom3 ran3;
    
    TFile *input =new TFile("xBQ2binning.root","read");
    
    int nparMINUIT=4;
    TCanvas *c1= new TCanvas();
    const int nxB=3, nQ2=4;
    const double QSqBinLow[nQ2+1]={1.5,2.5,4.0 ,6.0 ,10.0};
    TF2 * fLikeli[nQ2][nxB];
    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 zeta_min=0.0, zeta_max=1.0;
    double zeta;
    double xBMinBin, xBMaxBin;
    double Norm[nQ2][nxB], PhSp[nQ2][nxB];
    
    
    
    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);
    
    
    
    double Q2Min[nQ2][nxB],Q2Max[nQ2][nxB],Q2;
    double xMin0,xMax0,xMin[nQ2][nxB], xMax[nQ2][nxB];
    double Q2sim, xBjsim, ysim, yEvt;
    
    
    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[nQ2][nxB];
    
    //    Put event loop inside bin loop for minimization purposes
    int entries= xBQ2binning_tree->GetEntries();
    printf ("entries = %d \n", entries);
    int nbin;
    for (int jQ2Bin=0; jQ2Bin<4; jQ2Bin++) {
        for (int ixBin=0; ixBin<3; ixBin++) {
            double logli=0.0, model;
            // Try moving this outside jQ2Bin and ixBin loops
            auto myFCN = [&](const double *par){
                TRandom ran3;
                int jQ2Bin = par[4];
                int ixBin  = par[5];
                nbin=0;
                
                /*
                 for (int ipar=0; ipar<nparMINUIT+2; ipar++){
                 printf (" %7.3f, ", par[ipar]);
                 }
                 printf ("\n");
                 */
                // Recalculate normalization factor with each iteration
                Norm[jQ2Bin][ixBin] = 0.0;
                PhSp[jQ2Bin][ixBin]= 0.0;
                for (int iSimu=0; iSimu<nsimu; iSimu++){
                    Q2sim = QSqBinLow[jQ2Bin] +
                    (QSqBinLow[jQ2Bin+1]-QSqBinLow[jQ2Bin])*ran3.Uniform();
                    xBMinBin = (xBjMin(Q2sim)*(nxB-ixBin)+ xBjMax(Q2sim)*ixBin )/(double)nxB;
                    xBMaxBin = (xBjMin(Q2sim)*(nxB-ixBin-1)+ xBjMax(Q2sim)*(ixBin+1) )/(double)nxB;
                    xBjsim = xBMinBin+(xBMaxBin-xBMinBin)*ran3.Uniform();
                    if (xBjsim<=0.0){
                        printf("xBjsim = %10.3f \n", xBjsim);
                        break;
                    }
                    ysim   = Q2sim/(xBjsim*2.0*kBeam*sqrt(MpSq));
                    model  = TMath::Power(xBjsim,par[0]) * TMath::Power((1-xBjsim),par[1]) * TMath::Power((1+ (Q2sim*Q2sim)/ par[3]),-par[2]);
                    model *= (1-ysim+ysim*ysim/2.)/(xBjsim*Q2sim*Q2sim);
                    if(model <=0){
                        printf("model = %10.3f \n", model);
                        break;
                    }
                    Norm[jQ2Bin][ixBin] += (model)*(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;
                /*
                 printf("jQ2, ixB, model, xBMaxBin, xBMinBin = %02d, %02d, %10.5f, %10.5f, %10.5f \n",
                 jQ2Bin, ixBin, model, xBMaxBin, xBMinBin );
                 */
                //  printf("jQ2, ixB, model, Norm, PhSp = %02d, %02d, %10.5f, %10.5f, %10.5f \n",
                //   jQ2Bin, ixBin, model, Norm[jQ2Bin][ixBin], PhSp[jQ2Bin][ixBin] );
                
                
                
                //  if (fChain == 0) return;
                for(int ientry=0; ientry < entries;++ientry)
                {
                    
                    xBQ2binning_tree->GetEntry(ientry);
                    // Select events from bin (ixBin, jQ2Bin)
                    if ((ixB==ixBin)&&(jQ2==jQ2Bin)) {
                        // replace with explicit function
                        if ((QSq==0.0)||(xBj==0)){
                            printf ("event %d, QSq, xBj = %7.3f, %7.3f", ientry, QSq, xBj);
                            break;
                        }
                        yEvt = QSq/(xBj*2.0*kBeam*sqrt(MpSq));
                        logli += log(TMath::Power(xBj,par[0]) * TMath::Power((1-xBj),par[1]) * TMath::Power((1.0+ (QSq*QSq)/ par[3]),-par[2]));
                        logli += log(1.0 - yEvt + yEvt*yEvt/2.0)-log(xBj*QSq*QSq);
                        nbin++;
                    }
                }// event loop
                logli *= -1.0;
                logli += log(Norm[jQ2Bin][ixBin]/PhSp[jQ2Bin][ixBin])*nbin;
                
                return logli;
            }; // auto myFCN
            ROOT::Math::Functor fcn(myFCN,6);
            minimum[jQ2Bin][ixBin] = ROOT::Math::Factory::CreateMinimizer("Minuit", "Migrad");
            minimum[jQ2Bin][ixBin]->SetTolerance(0.01);
            minimum[jQ2Bin][ixBin]->SetPrintLevel(3);
            minimum[jQ2Bin][ixBin]->SetFunction(fcn);
            
            
            
            double step[6] = {0.2 , 0.2 , 0.2 , 0.2, 0.01, 0.01};
            double variable[6] = { -0.5, 2.0 , 0.0 , 1.0, (double)jQ2Bin, (double)ixBin }; // starting point
            minimum[jQ2Bin][ixBin]->SetVariable(0,"alpha",variable[0], step[0]);
            minimum[jQ2Bin][ixBin]->SetVariable(1,"beta",variable[1], step[1]);
            minimum[jQ2Bin][ixBin]->SetVariable(2,"gamma",variable[2], step[2]);
            //minimum[jQ2Bin][ixBin]->SetVariable(3,"lambda",variable[3], step[3]);
            //minimum[jQ2Bin][ixBin]->SetVariable(4,"jQ2Bin",variable[4], step[4]);
            //minimum[jQ2Bin][ixBin]->SetVariable(5,"ixBin",variable[5], step[5]);
            minimum[jQ2Bin][ixBin]->SetFixedVariable(3,"lambda",variable[3]);
            minimum[jQ2Bin][ixBin]->SetFixedVariable(4,"jQ2Bin",variable[4]);
            minimum[jQ2Bin][ixBin]->SetFixedVariable(5,"ixBin",variable[5]);
            minimum[jQ2Bin][ixBin]->SetVariableLimits(0,-2.0,2.0);
            minimum[jQ2Bin][ixBin]->SetVariableLimits(1,0.1,5.0);
            minimum[jQ2Bin][ixBin]->SetVariableLimits(2,-1.0,1.0);
            
            // do the minimization
            minimum[jQ2Bin][ixBin]->Minimize();
            printf("-LogLi[jQ2=%02d,ixB=%02d]= %10.4f from %d events\n", jQ2Bin, ixBin, minimum[jQ2Bin][ixBin]->MinValue(), nbin);
            
            
            
            
            
            
        } // for loop on ixBin
    } // for loop on jQ2Bin
    
} //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;
}