#include "TF2.h"
#include "TH2.h"
#include "TMath.h"

// Fitting a 2-D histogram
// This tutorial illustrates :
//  - how to create a 2-d function
//  - fill a 2-d histogram randomly from this function
//  - fit the histogram
//  - display the fitted function on top of the histogram 
//
// This example can be executed via the interpreter or ACLIC
//   root > .x fit2.C
//   root > .x fit2.C++
//Author: Rene Brun
         
Double_t g2(Double_t *x, Double_t *par) {
   Double_t r1 = Double_t((x[0]-par[1])/par[2]);  
   Double_t r2 = Double_t((x[1]-par[3])/par[4]);
   return par[0]*TMath::Exp(-0.5*(r1*r1+r2*r2));
}   
//This is a gaussian of the form f(x,y) = A^(-([(x-x[0])^2/2*sigmax^2]+[(y-y[0])^2/2*sigmsy^2]). You can see this format in the equation given by return above.
//where A = amplitude (x[0],y[0]) is the center, and sigmax and sigmay are the spreads of the "blob"
//In general a 2D elliptical Gaussian has the form f(x,y) = A^(-(a(x-x[0])^2+2b(x-x[0])(y-y[0])+c(y=y[0])^2)) where matrix [a b b c] is positive definite

Double_t fun2(Double_t *x, Double_t *par) {
   Double_t *p1 = &par[0];
   Double_t *p2 = &par[5];
   Double_t *p3 = &par[10];
   Double_t result = g2(x,p1) + g2(x,p2) + g2(x,p3); //this adds the three parts of the function together to create three Gaussians
   return result;
}

void fit2() {
   const Int_t npar = 15;
   Double_t f2params[npar] = 
      {100,-3,3,3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7}; //I believe these are three sets of parameters for the three humps on the graph.
	  //format: amplitude, x[0], sigmax, y[0], sigmay
	  //first set: big hump at -3,-3	second set: sharp peak at 0,0	third set: small hill at 4,4
   TF2 *f2 = new TF2("f2",fun2,-10,10,-10,10, npar);
   f2->SetParameters(f2params);

   //Create an histogram and fill it randomly with f2
   TH2F *h2 = new TH2F("h2","from f2",40,-10,10,40,-10,10);
   Int_t nentries = 100000;
   h2->FillRandom("f2",nentries);
   //Fit h2 with original function f2
   
   //This part is redefining parameters [0], [5], and [10]; in other words, the amplitudes. It appears to be scaling them by the number of entries. It has the most effect on the small hill.
   //to be entirely honest, I see no effect at all when I change it. However, removing it entirely screws everything up.
   Float_t ratio = 4*nentries/100000;
   std::cout<<"ratio = "<<ratio<<std::endl; //upon examination, the ratio is just 4 because nentries = 100000.
   f2params[ 0] *= ratio;
   f2params[ 5] *= ratio;
   f2params[10] *= ratio;
   f2->SetParameters(f2params);
   
   
  h2->GetXaxis()->SetTitle("X Axis");
  h2->GetYaxis()->SetTitle("Y Axis");
   h2->Fit("f2");
   h2->Draw("SURF2");
   f2->Draw("SURF same");
   cout<<"Inegral of f2: "<<f2->Integral(-10,10,-10,10)<<endl;
   
   
   
   
}