//
//  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"

//   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
//
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;
    
    double Norm[nQ2][nxB], PhSp[nQ2][nxB];
    double  NBin[nQ2][nxB];
    TRandom3 ran3;
    double xBj,QSq;
    
    double logli=0.0;
    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++) {
            /*
            tntuple.Loop(Norm[jQ2Bin][ixBin]/PhSp[jQ2Bin][ixBin], jQ2Bin, ixBin,
                         fLikeli[jQ2Bin][ixBin]);
             */
            
            gMinuit[jQ2Bin][ixBin] = new TMinuit(nparMINUIT); //4 parameters
            
            gMinuit[jQ2Bin][ixBin]->SetFCN(tntuple.Loop);
            
            double arglist[4];
            int ierflg = 0;
            
            arglist[0] = 1;
            gMinuit[jQ2Bin][ixBin]->mnexcm("SET ERR", arglist ,1,ierflg);
            
            // Set starting values and step sizes for parameters
            // Double_t alpha[nPar] = {-1.0,-0.8,-0.6,-0.4,-0.2,0.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};
            // Double_t gamma[nPar] = {1.,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.7,3.0};
            // Double_t lamda[nPar] = {1.,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.7,3.0};
            
            static double vstart[4] = {-0.5, 2.0 , 2.0 , 1.0};
            static double step[4] = {0.2 , 0.2 , 0.2 , 0.2};
            gMinuit[jQ2Bin][ixBin]->mnparm(0, "alpha", vstart[0], step[0], -1.0,0.8,ierflg);
            gMinuit[jQ2Bin][ixBin]->mnparm(1, "beta", vstart[1], step[1], 1.0,3.0,ierflg);
            gMinuit[jQ2Bin][ixBin]->mnparm(2, "gamma", vstart[2], step[2], 1.0,3.0,ierflg);
            gMinuit[jQ2Bin][ixBin]->mnparm(3, "lamda", vstart[3], step[3], 1.0,3.0,ierflg);
            
            
            
            /*     // Now ready for minimization step
             arglist[0] = 500;
             arglist[1] = 1.;
             gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
             
             // Print results
             double amin,edm,errdef;
             int nvpar,nparx,icstat;
             gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
             //gMinuit->mnprin(3,amin);*/
            
        }
        
        
    }
}





void fcn(int &npar, double *gin, double &f, double *par, int iflag)
 {
 
 }