//Code solely for fitting the distributions, this uses FCN method
// Output is from Momentum AOD analyzer


// This is needed for calling standalone classes (not needed on RACF)
#define _VANILLA_ROOT_

// C++ headers
#include <iostream>

// ROOT headers
#include <cstdio>
#include <vector>
#include <utility>
#include "Riostream.h"
#include "TObject.h"
#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TString.h"
#include "TVector3.h"
#include "TF1.h"
#include "TH1D.h"
#include "TMinuit.h"

Double_t z[50],x[50],y[50],errorz[50];

//______________________________________________________________________________
Double_t gamma(Double_t x, Double_t *par){
        
    Double_t alpha = par[0];
    Double_t beta = par[1];
    Double_t f1 = TMath::Gamma(alpha);
    Double_t f2 = pow(x,(alpha-1))*exp((-x)/beta); 
    Double_t f3 = pow(beta,alpha);
    Double_t f4 = f2/(f1*f3);
    return f4;
}

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

    //calculate chisquare
   double chisq = 0;
   double delta;
Int_t nbins =50;

   for (Int_t i=0;i<nbins; i++) {

        delta = (z[i]- gamma(x[i],par))/errorz[i];
        chisq+= delta*delta;
        }
    f = chisq;        

}



//______________________________________________________________________________



void myPadSetUp(TPad *currentPad, float currentLeft=0.11, float currentTop=0.04, float currentRight=0.04, float currentBottom=0.15){
                currentPad->SetLeftMargin(currentLeft);
                currentPad->SetTopMargin(currentTop); 
                currentPad->SetRightMargin(currentRight); 
                currentPad->SetBottomMargin(currentBottom); 
                return;
}
//______________________________________________________________________________





void Pro_Gamma_fitting(){

TFile *file = TFile::Open("./Avghists_charge_dists.root");
//file->ls(); 

TH1F *hAvgPt;

hAvgPt = (TH1F*)file->Get("hAvgPt_0");


TCanvas *c1 = new TCanvas("c1","Avg Momentum Distribution",20,15,800,600);
c1->cd();
TPad *myPad1 = new TPad("myPad1", "The pad",0,0,1,1); myPadSetUp(myPad1,0.15,0.04,0.04,0.15);
myPad1->SetLogy();
myPad1->Draw();
myPad1->SetTickx(1);
myPad1->SetTicky(1);
myPad1->cd();

TH2F *fhist = new TH2F("fhist","Average Momentum Distributions",50,0.4,.8,200,1e0,1e4);
  fhist->GetXaxis()->SetTitle("<p_{T}> (GeV/c) ");
  fhist->GetYaxis()->SetTitle(" Counts");
  fhist->GetXaxis()->CenterTitle();
  fhist->GetYaxis()->CenterTitle();
  fhist->SetTickLength(0.02,"x");
  fhist->SetTickLength(0.02,"y");
  fhist->SetTitleOffset(1.30,"y");
  fhist->SetTitleOffset(1.5,"x");
  fhist->GetYaxis()->SetNdivisions(5);
  fhist->Draw("");

gStyle->SetOptStat(0);
c1->GetFrame()->SetFillColor(0);
c1->GetFrame()->SetBorderSize(12);
c1->SetLogy();

Int_t NBins_x = hAvgPt->GetNbinsX();
cout << NBins_x <<  endl;
TGraphErrors *grAvgPt[1];
//Double_t z[NBins_x],errorz[NBins_x],x[NBins_x];
TAxis *xaxis = hAvgPt->GetXaxis();


for(Int_t i =0;i<NBins_x;i++){
    
    Double_t binCenter = xaxis->GetBinCenter(i+1);
    z[i] = hAvgPt->GetBinContent(i+1);
    errorz[i] = hAvgPt->GetBinError(i+1);
    x[i] = binCenter;

    cout << x[i] << z[i] << errorz[i] << endl;
}



//grAvgPt[0] = new TGraphErrors(NBins_x, x,y,0,y_err);

TMinuit *gMinuit = new TMinuit(3);  //initialize TMinuit with a maximum of 5 params
   gMinuit->SetFCN(fcn);

double arglist[10];
   int ierflg = 0;
 
   arglist[0] = 1;
   gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);


// Set starting values and step sizes for parameters
static double vstart[2] = {561.451, 0.00110568};
static double step[2] = {1 ,0.001};

gMinuit->mnparm(0, "alpha", vstart[0], step[0], 0,0,ierflg);
gMinuit->mnparm(1, "beta", vstart[1], step[1], 0,0,ierflg);


// Now ready for minimization step
arglist[0] = 5000;
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);




}