#include <iostream>
#include <iomanip>
#include <fstream>
#include "TRandom.h"
#include "TH1F.h"
#include "TFile.h"
#include "TF1.h"
#include "TMinuit.h"
using namespace std;
#include "fw_binned.cc"


void fitfunction(Int_t &npar, Double_t *gin, Double_t &f,Double_t *par,Int_t iflag);
int binned(TH1F *h1, TH1F *h2, TH1F *DATA, double &f1, double &ef1, double &f2, double &ef2);
TH1F *hist1;
TH1F *hist2;
TH1F *hPD;

//int fw_binned(int *a1, int *a2, int *ad, double &f1, double &ef1, double &f2, double &ef2){


int binned_fit(){
  //-------------
  // Get Template
  //-------------
  TFile *fin = new TFile("template.root", "read");
  TH1F *h1;
  TH1F *h2;
  h1 = dynamic_cast<TH1F*>(fin->Get("hist1"));
  h2 = dynamic_cast<TH1F*>(fin->Get("hist2"));

  //-------------------------------------
  // Make Pseudo Data (set Nev and Frac1)
  //-------------------------------------
  int Nev     = 4000;
  float Frac1 = 0.2;
  int Nev1 = (int)(Nev*Frac1);
  int Nev2 = Nev - Nev1;
  TH1F *hist_d = new TH1F("hist_d", "hist_d", 100, 0, 100);
  hist_d->SetMarkerStyle(21);
  hist_d->SetLineWidth(2);

  for(int iev=0; iev<Nev1; iev++) hist_d->Fill( h1->GetRandom() );
  for(int iev=0; iev<Nev2; iev++) hist_d->Fill( h2->GetRandom() );
  double f1=0, ef1=0, f2=0, ef2=0;

  //HH binned minuit fit HH//
  binned(h1, h2, hist_d, f1, ef1, f2, ef2);
  cout<<" Minuit "<<endl;
  cout<<"par1 = "<<f1<<"+-"<<ef1<<endl;
  cout<<"par2 = "<<f2<<"+-"<<ef2<<endl;



  int a1[100], a2[100], ad[100];

  cout<<h1->GetEntries()<<"   "<<hist_d->GetEntries()<<endl;
  for(int i=0; i<100; i++){
    a1[i]=h1->GetBinContent(i+1);
    a2[i]=h2->GetBinContent(i+1);
    ad[i]=hist_d->GetBinContent(i+1);
  }
  fw_binned(a1, a2, ad, f1, ef1, f2, ef2);
  cout<<" My FW fit "<<endl;
  cout<<"par1 = "<<f1<<"+-"<<ef1<<endl;
  cout<<"par2 = "<<f2<<"+-"<<ef2<<endl;



  return 1;
}

 int binned(TH1F *h1, TH1F *h2, TH1F *DATA, double& f1, double& ef1, double& f2, double& ef2){
   hist1 = h1;
   hist2 = h2;
   hist1->Sumw2();
   hist2->Sumw2();
   hPD = DATA;

   //------
   // Rebin
   //------
   int Nrebin = 1;
   hist1->Rebin(Nrebin);
   hist2->Rebin(Nrebin);
   hPD  ->Rebin(Nrebin);

   //-------------
   // Minuit setup
   //-------------
   //f = -2log
   const Int_t npar = 2;
   Int_t ier=0;
   TMinuit *gMinuit = new TMinuit(npar);
   gMinuit->SetPrintLevel(-1);
   gMinuit->SetFCN(fitfunction);

   Double_t vstart = 0.5;
   Double_t step = 0.1;
   Double_t min = 0.00;
   Double_t max = 1.00;
   Double_t fpar;
   Double_t ferr;

   // l minmum? //
   gMinuit->mnparm(0, "test1", vstart, step, min, max, ier);
   gMinuit->mnparm(1, "test2", vstart, step, min, max, ier);
   Double_t arg[] = {2000,1.0};
   gMinuit->mnexcm("MIGRAD", arg, 2, ier);
   cout<< "Error flag = " << ier <<endl;

   // output //
   Int_t aa = 0;
   TString bb;
   Double_t cc;
   Double_t dd;
   Double_t ee;
   Double_t ff; 
   Int_t gg;
   gMinuit->mnpout(0,bb,cc,dd,ee,ff,gg);
   //   cout << bb << " " << cc << " " << dd << " " << ee << " " << ff << " " << gg << endl;
   f1=cc; ef1=dd;
   gMinuit->mnpout(1,bb,cc,dd,ee,ff,gg);
   //   cout << bb << " " << cc << " " << dd << " " << ee << " " << ff << " " << gg << endl;
   f2=cc; ef2=dd;
   return 0;
 }

 double BinFrac(int ih, int ibin){
   double Frac = 0;

   if(ih==0){ // hist1
     double TotNev  = hist1->Integral();
     double Content = hist1->GetBinContent(ibin);
     Frac = Content/TotNev;

   }
   if(ih==1){ // hist2
     double TotNev  = hist2->Integral();
     double Content = hist2->GetBinContent(ibin);
     Frac = Content/TotNev;
   }
   return Frac;
 }

 void fitfunction(Int_t &npar, Double_t *gin, Double_t &f,Double_t *par,Int_t iflag){
   Double_t f1=par[0];
   Double_t f2=par[1];
   Double_t logL=0;

   int Nbin = hPD->GetNbinsX();
   for(Int_t ibin=1; ibin<=Nbin ; ibin++){
     double mu = hPD->Integral() * ( f1*BinFrac(0,ibin) + f2*BinFrac(1,ibin) );
     double N  = hPD->GetBinContent(ibin);
     double logP = 0;
     if(mu>0){
       logP = -mu + N*log(mu) - log(TMath::Factorial((int)N));
       //       cout<<"logP = "<<logP<<" = "<<-mu<<" + "<<N*log(mu)<<" - "<<log(TMath::Factorial((int)N))<<endl;
     }
     logL+=logP;
   }
   f=-2*logL;
   //   cout <<"f and f1 =========== "<< f <<"     "<<f1<<endl;
 }



void QuenchedBraggFitFunc(int &npar, doulbe *gin, double &f, double *par int iflag){
  double R0      = par[0];
  double sigma   = par[1];
  double phi0    = par[2];
  double xOffset = 0;
  double birks   = par[3];

  double highL=35767;
  if (R0<=0.) return highL;
  if (sigma<=0.) return highL;
  if (phi0<=0.) return highL;
  if (xOffset<0.) return highL;
  if (birks<=0.) return highL;
  if (R0 != R0) return highL; //Check if R0 is nan
  if (z[0] != z[0]) return highL; //Check if z is nan
  z[0] = z[0] + xOffset; //add xOffset
  if (z[0] < 0) return highL;



  double zDef = 5.*sigma; //This constant determines in which area the integrator of fInt1 and fInt2 gonna be defined and integrated, see definition of a and b
  double a = min(z[0]-zDef,(1.-precision-1e-3)*R0); // integration start. Subtract constant to make sure a is smaller than b
  double b = min(z[0]+zDef,(1.-precision)*R0); // integration end
  // HH Bo integral 1 part HH// ---------------------------------------------------------------------------------------------------------
  double width=(b-a)/my_precision;
  double xx=a+width/2.;
  double integral1=0;
  for(int i=0; i<my_precision; i++){
    double preFactor = (1.+beta*(R0-xx))/(1+beta*R0);
    double foldingTerm = 1./(2.507*sigma)*exp(-(z[0]-xx)*(z[0]-xx)/(2*sigma*sigma));    
    double tdLdz=1./(0.2117*pow(R0-xx,0.425)+birks);
    double val = preFactor * tdLdz * foldingTerm;
    integral1+=val*width;
    xx+=width;
  }
  // HH Bo integral 2 part HH// ---------------------------------------------------------------------------------------------------------
  double LofZ = 0;
  double myLofZ=0;
  if(z[0]<(1-precision)*R0){
    double width=((1-precision)*R0-z[0])/my_precision;
    double xx=z[0]+width/2.;
    for(int i=0; i<my_precision; i++){
      double val  = 1./(0.2117*pow(R0-xx,0.425)+birks);
      myLofZ+=val*width;
      xx+=width;
    }
  }
  double pF = gamma_H2O*beta/(1+beta*R0);
  double integral2 = pF * myLofZ;

  // HH Eo integral 2 part HH// ---------------------------------------------------------------------------------------------------------
  double preFactor = S*phi0/rho;           // Put all parts together and add scaling factor
  double integral = integral1 + integral2; // Put all parts together and add scaling factor
  // Combine prefactor and integral of all parts
  double val = preFactor*integral;
  f=val;

}