//
//  xBQ2binning.cpp
//
//
//
//

#include <stdio.h>
#define xBQ2binning_tree_cxx
#include "xBQ2binning_tree.h"
#include "xBQ2binning_tree.C"

#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"

//   In a ROOT session, you can do:
//      root> .L xBQ2binning_tree.C
//      root> xBQ2binning_tree t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//double fcn(double *par);
//void myFCN(int &npar, double *gin, double &logli, double *par, int iflag);
void xBQ2binning(){
    // contained in .h file for this Class
    //TFile *input =new TFile("xBQ2binning_simu.root","read");
    //  Put event sum of logli function here for just one jQ2bin, ixBin value
    int jQ2,ixB;
    int nbin=0;
    double Norm[nQ2][nxB], PhSp[nQ2][nxB];
    double  NBin[nQ2][nxB];
    TRandom3 ran3;
    double xBj,QSq;
    TTree          *fChain;
    
    ROOT::Math::Minimizer* minimum[nQ2][nxB];
    //xBQ2binning_tree tntuple;
    TMinuit *gMinuit[nQ2][nxB];
    int nparMINUIT=4;
    TNtuple *xBQ2binning_tree = new TNtuple("xBQ2binning_tree","Branches","xBj:QSq:ixB:jQ2");
    char fName[64];
    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);
            // revise to different parameter sets for each Q2 x xB bin
            fLikeli[jQ2Bin][ixBin]->SetParameters(-0.5,2.0,2.0,1.0);
        }
    }
    
    
    double xBMinBin, xBMaxBin;
    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;
            
        }
    }
    //    Put event loop inside bin loop for minimization purposes
    for (int jQ2Bin=0; jQ2Bin<nQ2; jQ2Bin++) {
        for (int ixBin=0; ixBin<nxB; ixBin++) {
            
            
            
            
            
            auto myFCN = [&](const double *par){
                
                double logli=0.0;
                int jQ2Bin = par[4];
                int ixBin  = par[5];
              //  if (fChain == 0) return;
                Long64_t nentries = fChain->GetEntries();
             //   Long64_t nentries = fChain->GetEntriesFast();
                fChain->GetEntriesFast();
                // int nentries= xBQ2binning_tree->GetEntries();
                
                Long64_t nbytes = 0, nb = 0;
                for (Long64_t jentry=0; jentry<nentries;jentry++) {
                    Long64_t ientry = fChain->LoadTree(jentry);  //Long64_t ientry = LoadTree(jentry);
                    
                    if (ientry < 0) break;
                    nb = fChain->GetEntry(jentry);   nbytes += nb;
                    
                    // Select events from bin (ixBin, jQ2Bin)
                    if ((ixB==ixBin)&&(jQ2==jQ2Bin)) {
                        // replace with explicit function
                        logli += log(TMath::Power(xBj,par[0]) * TMath::Power((1-xBj),par[1])*TMath::Power((1+ (QSq*QSq)/ par[3]),-par[2]))-log(Norm[jQ2Bin][ixBin]); //not sure about the normalization is correct here
                        nbin++;
                    }
                }
                logli *= -1.0/nbin;
                //       printf("LogLi[jQ2=%02d,ixB=%02d]= %10.3f \n", jQ2Bin, ixBin, logli);
                
                return logli;
                
            };
            ROOT::Math::Functor f(myFCN,2);
            minimum[jQ2Bin][ixBin] = ROOT::Math::Factory::CreateMinimizer("Minuit", "Migrad");
            minimum[jQ2Bin][ixBin]->SetTolerance(0.01);
            minimum[jQ2Bin][ixBin]->SetPrintLevel(3);
            minimum[jQ2Bin][ixBin]->SetFunction(f);
            
            double step[6] = {0.2 , 0.2 , 0.2 , 0.2, 0.01, 0.01};
            double variable[6] = { -0.5, 2.0 , 2.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);
            //  minimum[jQ2Bin][ixBin]->SetFixedVariable(4);
             //minimum[jQ2Bin][ixBin]->SetFixedVariable(5);
            
            
            // do the minimization
            minimum[jQ2Bin][ixBin]->Minimize();
            //   minimum[jQ2Bin][ixBin]->SetLowerLimitedVariable(0,"k1",variable[0], step[0],0.); //Lower value is energy=0
            // minimum[jQ2Bin][ixBin]->SetLowerLimitedVariable(1,"k2",variable[1], step[1],0.);
            //minimum[jQ2Bin][ixBin]->Minimize(); // do the minimization
            
        }
        
        
    }
}