#include "TCanvas.h"

#include "TFile.h"

#include "TObject.h"

#include "TString.h"

#include "TF1.h"

#include "TMath.h"

#include "TStyle.h"

#include "TChain.h"

#include "TF2.h"

#include "TF12.h"

#include "TH1D.h"

#include "TList.h"

#include "TLegend.h"

#include "TLine.h"

#include <cmath>

#include "TGraphAsymmErrors.h"

#include "Math/Integrator.h"

#include <memory>


struct FitFunction {

   std::shared_ptr<TF2> main;
   std::shared_ptr<TF1> part;
   std::shared_ptr<TF12> f12A;
   std::shared_ptr<ROOT::Math::IntegratorOneDim> ig; 

  FitFunction() {
    //main = std::make_shared<TF2>("main","[0]*(pow(y,x)*exp(-y)/TMath::Gamma(x+1))/exp(([1]*pow(y,1))+([2]*pow(y,2))+ ([3]*pow(y,3)))");//don't work at high x, use Stirling's formula
    main = std::make_shared<TF2>("main","[0]*pow(y*exp(1)/x,x)/sqrt(2*TMath::Pi()*x)/(1+1/12./x)/exp(y+([1]*pow(y,1))+([2]*pow(y,2))+ ([3]*pow(y,3)))");
    f12A = std::make_shared<TF12>("f12A",main.get(),0,"y");
    part = std::make_shared<TF1>("part", "1.0/exp(([0]*pow(x,1))+([1]*pow(x,2))+ ([2]*pow(x,3)))");
    ig   = std::make_shared<ROOT::Math::IntegratorOneDim>(ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR);
  }

  double operator()(double * x, double * p) {
      main -> SetParameters(p[0], p[1], p[2], p[3]);
      f12A -> SetXY(x[0]); 
      ig -> SetFunction( * f12A);
      //double mainsum = ig -> IntegralUp(0.0); // don't work
      double b = 400; // use any b>=400
      double mainsum = ig->Integral(0.,b);

      part -> SetParameters(p[1],p[2],p[3]);
      ig -> SetFunction( * part);
      //double deno = ig -> IntegralUp(0.0);
      double deno = ig->Integral(0.,b);

      return mainsum / deno;
    }
};

void TestCodeMod() { // fit in the range n = 250-350

  gStyle -> SetOptStat(1);
  gStyle -> SetOptFit(1);

  TCanvas * canvas = new TCanvas("graph", "Graph");
  canvas -> SetLogy(); 
  //Data
  TFile * File = TFile::Open("DataFile.root");
  if (!File) printf("File Not Found"); 
  TGraphAsymmErrors * gr = dynamic_cast < TGraphAsymmErrors * > (File -> Get("MyData"));
  gr -> Draw("ap"); 

  FitFunction* ff = new FitFunction() ;
  //My Fit
  //TF1 * MyFit = new TF1("MyFit",ff, 0,70, 4);
  TF1 * MyFit = new TF1("MyFit",ff, 250,350, 4);
  MyFit ->FixParameter(0, 1.0);
  //MyFit ->FixParameter(1, 1.620*pow(10,-1));
  //MyFit ->FixParameter(2, -1.35*pow(10,-3));
  //MyFit ->FixParameter(3, 2.589*pow(10,-6));
  MyFit ->SetParameters(1, 0.2577, -1.836e-3, 3.40e-6);

  MyFit -> SetParNames("N", "#alpha", "#beta", "#gamma");
  MyFit -> SetLineColor(1);
  //MyFit->Draw("l");
  MyFit->Draw("same");
  gr -> Fit(MyFit, "RME+");

}