#define parameter_Fit_cxx
#include "parameter_Fit.h"
#include <TStyle.h>
#include <TCanvas.h>
#include <TChain.h>
#include <iostream>
#include <fstream>
#include <cstdio>
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TSystem.h"
#include "TH1.h"
#include "THnSparse.h"
#include "TCanvas.h"
#include "TStopwatch.h"
#include "TF1.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TVirtualFitter.h"
#include "TFitter.h"

void stats(THnSparseD* fSparse, Double_t &f, Double_t *par);

double get_bin_integral(Double_t* bin_low_edge, Double_t* bin_up_edge, Int_t* bin_n, Double_t *par);
double get_prediction(Double_t* angles, Double_t *par);

void parameter_Fit::Loop()
{

  Int_t bin_run, bin_i;
  Int_t bin_num;
  Int_t actual_bin_num;
  Int_t Dim_n;
  Int_t k;
  
  Int_t bins[5];
  Double_t COSTHET_PI_edges[4];
  Double_t COSTHET_PSI_edges[4];
  Double_t COSTHET_X_edges[4];
  Double_t DELTA_PHI_edges[5];
  Double_t SIGMA_PHI_edges[5];
  Double_t arrBinEdges[5][5];
  Double_t Data[kMaxevtarray][5];
  
  bins[0] = 4;
  bins[1] = 4;
  bins[2] = 4;
  bins[3] = 5;
  bins[4] = 5;
  
  bin_num = bins[0]*bins[1]*bins[2]*bins[3]*bins[4];
  Int_t nbins[5] = {10, 10, 10, 10, 10};
  Double_t xmin[5] = {-1,-1,-1, 0, 0};
  Double_t xmax[5] = {1, 1, 1, TMath::TwoPi(), TMath::TwoPi()};

//   In a ROOT session, you can do:
//      Root > .L parameter_Fit.C
//      Root > parameter_Fit 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
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<1;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;

      THnSparseD* fSparse = new THnSparseD("experiment", "experiment", 5, bins,xmin,xmax );
      Dim_n = fSparse->GetNdimensions();
  
      for(k=0;k<kMaxevtarray;k++){

        Data[k][0] = evtarray_COSTHET_PI[k]; 
        Data[k][1] = evtarray_COSTHET_PSI[k];
        Data[k][2] = evtarray_COSTHET_X[k];
        Data[k][3] = evtarray_DELTA_PHI[k];
        Data[k][4] = evtarray_SIGMA_PHI[k];
        
      }

        fSparse->Fill(Data[k],1.0);
           
     }
      
     TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 5);
     minuit->SetFCN(stats);
     minuit->SetParameter(0, "constant", 1,   100,  0, 1000);
     
     Double_t arglist[10];
     arglist[0] = 1;
     minuit->ExecuteCommand("SET PRINT", arglist, 1);
     minuit->ExecuteCommand("MIGRAD", arglist, 0);
      
      
      
}

void
LetsStart() {

  parameter_Fit *Ntuple = new parameter_Fit();

  if (Ntuple != NULL) {
    Ntuple->Loop();
  }
  else {
    printf("Error with ntuple\n");
  }
}

void stats(THnSparseD* fSparse, Double_t &f, Double_t *par)
{

   Double_t r;
   Int_t counter, counter1;
   Double_t c, eventn, bin_integral;
   
   double *low = new double[5];
   double *up = new double[5];
   int *idx = new int[5];
   
   Double_t COSTHET_PI_bin_low_edge[100], COSTHET_PSI_bin_low_edge[100], 
   COSTHET_X_bin_low_edge[100], DELTA_PHI_bin_low_edge[100], SIGMA_PHI_bin_low_edge[100];
   
   Double_t COSTHET_PI_bin_up_edge[100], COSTHET_PSI_bin_up_edge[100], 
   COSTHET_X_bin_up_edge[100], DELTA_PHI_bin_up_edge[100], SIGMA_PHI_bin_up_edge[100];
   
   Double_t COSTHET_PI_bin_center[100], COSTHET_PSI_bin_center[100], 
   COSTHET_X_bin_center[100], DELTA_PHI_bin_center[100], SIGMA_PHI_bin_center[100];
   
   Double_t COSTHET_PI_bin_w[100], COSTHET_PSI_bin_w[100], 
   COSTHET_X_bin_w[100], DELTA_PHI_bin_w[100], SIGMA_PHI_bin_w[100];
   
   f = 0;
   Double_t df_lnlike = 0;
   Double_t df_chi_squared = 0;
   Double_t func, norm;
   norm = 0;
   eventn = 0;
   counter = 0;
   counter1 = 0;
   
   Int_t entries_n = (int) fSparse->GetEntries();
  
   
//____________________assigning axis to histo_____________________________________________

   TAxis *COSTHET_PI = fSparse->GetAxis(0);
   TAxis *COSTHET_PSI = fSparse->GetAxis(1);
   TAxis *COSTHET_X = fSparse->GetAxis(2);
   TAxis *DELTA_PHI = fSparse->GetAxis(3);
   TAxis *SIGMA_PHI = fSparse->GetAxis(4);
   
//_____________________getting bin number+_________________________________________________
   
   Int_t n_COSTHET_PI  = COSTHET_PI->GetNbins();
   Int_t n_COSTHET_PSI = COSTHET_PSI->GetNbins();
   Int_t n_COSTHET_X   = COSTHET_X->GetNbins();
   Int_t n_DELTA_PHI   = DELTA_PHI->GetNbins();
   Int_t n_SIGMA_PHI   = SIGMA_PHI->GetNbins();
   
   Int_t bin_n[5] = {n_COSTHET_PI,n_COSTHET_PSI,n_COSTHET_X,n_DELTA_PHI,n_SIGMA_PHI};
   
   //cout<<n_DELTA_PHI<<"here";
   
//______________________getting bin edges and width______________________________________
   for (Int_t bin=1;bin<=n_COSTHET_PI;bin++) {
         COSTHET_PI_bin_low_edge[bin] = COSTHET_PI->GetBinLowEdge(bin);
         COSTHET_PI_bin_up_edge[bin] = COSTHET_PI->GetBinUpEdge(bin);
         COSTHET_PI_bin_center[bin] = COSTHET_PI->GetBinCenter(bin);
         COSTHET_PI_bin_w[bin] = COSTHET_PI->GetBinWidth(bin);
         //cout<< COSTHET_PI_bin_w[bin]<< " ";
         }
     
   for (Int_t bin=1;bin<=n_COSTHET_PSI;bin++) {  
         COSTHET_PSI_bin_low_edge[bin] = COSTHET_PSI->GetBinLowEdge(bin);
         COSTHET_PSI_bin_up_edge[bin] = COSTHET_PSI->GetBinUpEdge(bin);
         COSTHET_PSI_bin_center[bin] = COSTHET_PSI->GetBinCenter(bin);
         COSTHET_PSI_bin_w[bin] = COSTHET_PSI->GetBinWidth(bin);
         //cout<< COSTHET_PSI_bin_w[bin]<< " ";
        }
   
   for (Int_t bin=1;bin<=n_COSTHET_X;bin++) { 
         COSTHET_X_bin_low_edge[bin] = COSTHET_X->GetBinLowEdge(bin);
         COSTHET_X_bin_up_edge[bin] = COSTHET_X->GetBinUpEdge(bin);
         COSTHET_X_bin_center[bin] = COSTHET_X->GetBinCenter(bin);
         COSTHET_X_bin_w[bin] = COSTHET_X->GetBinWidth(bin);
         //cout<< COSTHET_X_bin_w[bin]<< " ";
         }
   
   for (Int_t bin=1;bin<=n_DELTA_PHI;bin++) { 
         DELTA_PHI_bin_low_edge[bin] = DELTA_PHI->GetBinLowEdge(bin);
         DELTA_PHI_bin_up_edge[bin] = DELTA_PHI->GetBinUpEdge(bin);
         DELTA_PHI_bin_center[bin] = DELTA_PHI->GetBinCenter(bin);
         DELTA_PHI_bin_w[bin] = DELTA_PHI->GetBinWidth(bin);
         //cout<< DELTA_PHI_bin_w[bin]<< " ";
         }
   
   for (Int_t bin=1;bin<=n_SIGMA_PHI;bin++) { 
         SIGMA_PHI_bin_low_edge[bin] = SIGMA_PHI->GetBinLowEdge(bin);
         SIGMA_PHI_bin_up_edge[bin]  = SIGMA_PHI->GetBinUpEdge(bin);
         SIGMA_PHI_bin_center[bin]  = SIGMA_PHI->GetBinCenter(bin);
         SIGMA_PHI_bin_w[bin]        = SIGMA_PHI->GetBinWidth(bin);
         //cout<<  SIGMA_PHI_bin_w[bin]<< " ";
         }
         

//_______________________calculating bin integral____________________

   
   for (Int_t i=1;i<=n_COSTHET_PI;i++) {
      low[0] = COSTHET_PI_bin_low_edge[i];
      up[0] = COSTHET_PI_bin_up_edge[i];
      idx[0] = i;
      
      for (Int_t j=1;j<=n_COSTHET_PSI;j++) {
         low[1] = COSTHET_PSI_bin_low_edge[j];
         up[1] = COSTHET_PSI_bin_up_edge[j];
         idx[1] = j;
         
         for (Int_t k=1;k<=n_COSTHET_X;k++) {
            low[2] = COSTHET_X_bin_low_edge[k];
            up[2] = COSTHET_X_bin_up_edge[k];
            idx[2] = k;       
            
            for (Int_t l=1;l<=n_DELTA_PHI;l++) {
               low[3] = DELTA_PHI_bin_low_edge[l];
               up[3] = DELTA_PHI_bin_up_edge[l];
               idx[3] = l;      
               
               
               for (Int_t m=1; m<=n_SIGMA_PHI; m++) {
                  low[4] = SIGMA_PHI_bin_low_edge[m];
                  up[4] = SIGMA_PHI_bin_up_edge[m];
                  idx[4] = m;      
                  
                  bin_integral = get_bin_integral(low, up, bin_n, par);
                                                     
                  func = 100 * bin_integral;
                 
                  norm += bin_integral;
                  counter ++;
                  
                  //cout << norm <<endl;
                  
                  
                  
                  if (func < 0) cout<<"theoretical value is negative!!!";
                  
                  c = fSparse->GetBinContent(idx);
                  
//_______________________calculating stats - chi^2 and log likelihood____________________                  
                  
                  
                  df_chi_squared = pow((c-func),2)/func;
                  
                  r = c/func;
                  
                  eventn += c; 
                  
                  f += df_chi_squared;
                  
               }
            }
         }
      }
   }
   
   delete [] low;
   delete [] up;
   delete [] idx;
   
}

double get_bin_integral(Double_t* bin_low_edge, Double_t* bin_up_edge, Int_t* bin_n, Double_t *par)
{

  Double_t COSTHET_PI_bin_low_edge, COSTHET_PI_bin_up_edge;
  Double_t COSTHET_PSI_bin_low_edge, COSTHET_PSI_bin_up_edge;    
  Double_t COSTHET_X_bin_low_edge, COSTHET_X_bin_up_edge;
  Double_t DELTA_PHI_bin_low_edge, DELTA_PHI_bin_up_edge;
  Double_t SIGMA_PHI_bin_low_edge, SIGMA_PHI_bin_up_edge;
  
  Int_t COSTHET_PI_bin_n;
  Int_t COSTHET_PSI_bin_n;
  Int_t COSTHET_X_bin_n;
  Int_t DELTA_PHI_bin_n;
  Int_t SIGMA_PHI_bin_n;
  
  Double_t Integral;
  
  Integral = 0.0;
  
  COSTHET_PI_bin_low_edge =  bin_low_edge[0];
  COSTHET_PI_bin_up_edge = bin_up_edge[0];
  COSTHET_PI_bin_n = bin_n[0];

  COSTHET_PSI_bin_low_edge = bin_low_edge[1];
  COSTHET_PSI_bin_up_edge = bin_up_edge[1];
  COSTHET_PSI_bin_n = bin_n[1];
  
  COSTHET_X_bin_low_edge = bin_low_edge[2];
  COSTHET_X_bin_up_edge = bin_up_edge[2];
  COSTHET_X_bin_n = bin_n[2];

  DELTA_PHI_bin_low_edge = bin_low_edge[3];
  DELTA_PHI_bin_up_edge = bin_up_edge[3];
  DELTA_PHI_bin_n = bin_n[3];
  
  SIGMA_PHI_bin_low_edge = bin_low_edge[4];
  SIGMA_PHI_bin_up_edge = bin_up_edge[4];
  SIGMA_PHI_bin_n = bin_n[4];
    
  Integral = 9.0/((24.0+8.0*par)*TMath::Pi())*(
  
        ((COSTHET_PI_bin_up_edge - 1.0/3.0*pow(COSTHET_PI_bin_up_edge,3)) - 
        (COSTHET_PI_bin_low_edge - 1.0/3.0*pow(COSTHET_PI_bin_low_edge,3)))*
  
        1.0/2.0*(
        (DELTA_PHI_bin_up_edge - DELTA_PHI_bin_low_edge)*
  
        ((COSTHET_PSI_bin_up_edge + 1.0/3.0*pow(COSTHET_PSI_bin_up_edge,3)) - 
        (COSTHET_PSI_bin_low_edge + 1.0/3.0*pow(COSTHET_PSI_bin_low_edge,3))) +
  
        par*(TMath::Sin(DELTA_PHI_bin_up_edge)*TMath::Cos(DELTA_PHI_bin_up_edge) - 
        TMath::Sin(DELTA_PHI_bin_low_edge)*TMath::Cos(DELTA_PHI_bin_low_edge))*
  
        ((COSTHET_PSI_bin_up_edge - 1.0/3.0*pow(COSTHET_PSI_bin_up_edge,3)) - 
        (COSTHET_PSI_bin_low_edge - 1.0/3.0*pow (COSTHET_PSI_bin_low_edge,3)))
        ));
        
  Integral = Integral*((SIGMA_PHI_bin_up_edge - SIGMA_PHI_bin_low_edge)/TMath::TwoPi())*
             ((COSTHET_X_bin_up_edge - COSTHET_X_bin_low_edge)/2.0);

  return  Integral;

}