#include "TF2.h"
#include "TH2F.h"
#include "TSpectrum2.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TMath.h"
#include "TROOT.h"

TSpectrum2 *s;
void  Owen2DFitExampleFunc()
{

TFile* fitexhist = new TFile("fitexhist.root", "RECREATE", "2DFitExampleFunc histograms");
s = new TSpectrum2();

//define canvases
c1 = new TCanvas("c1","c1",1);
c2 = new TCanvas("c2","c2",1);
c3 = new TCanvas("c3","c3",1);

std::cout<<"making 2D gaussian"<<std::endl;

  //Make A Random 2D Gaussian and set titles and other fun stuff
  TH2F h1("h1","2DGaussian",100,-3,3,100,-3,3);
  h1.SetOption("surf3z");
  h1.GetXaxis()->SetTitle("My X Title");
  h1.GetYaxis()->SetTitle("My Y Title");
  int binx1 = h1.GetXaxis()->FindBin(-.3);
  int binx2 = h1.GetXaxis()->FindBin(.3);
  int biny1 = h1.GetYaxis()->FindBin(-.3);
  int biny2 = h1.GetYaxis()->FindBin(.3);
  TRandom3 foo(12345);
  for(int i=0;i<10000;i++)
    {
      h1.Fill(foo.Gaus(0,1),foo.Gaus(0,1));
    }

std::cout<<"create 2D fit function"<<std::endl;
  //This is just defining a 2D gaussian explicity
  //TF2 *my2dfit = new TF2("my2dfit","[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])",-3,3,-3,3);
  TF2 *my2dfit = new TF2("my2dfit","(([0]/(TMath::Sqrt(TMath::TwoPi())))*(1/([2]*[4])))*(TMath::Exp(-0.5*((x-[1])/[2])^2))*(TMath::Exp(-0.5*((y-[3])/[4])^2))");
  my2dfit->SetParameter(2,1);
  my2dfit->SetParameter(4,1);
  my2dfit->SetParameter(1,0);
  my2dfit->SetParameter(3,0);

  
  TH2F *hnew = (TH2F*)h1.Clone("hnew");
 
  c1->cd();
  h1.Draw();
  h1.Fit("my2dfit");
  h1.Fit("my2dfit");
  h1.Fit("my2dfit");
  h1.Fit("my2dfit");
  h1.Write();
  
  //find peaks
  std::cout<<"finding peaks"<<std::endl;
  int nfound = s->Search(h1,2,"col");
  
  c2->cd();
  my2dfit->Draw();

  //  hnew->Divide(myfit);
int intmult = my2dfit->Integral(-100,100,-100,100);
//std::cout<<"Integral of function: "<<my2dfit->Integral(-100,100,-100,100)<<std::endl;  
std::cout<<"Integral of function: "<<intmult*TMath::Sqrt(TMath::TwoPi())*1<<std::endl;  
std::cout<<"Integral of histo: "<<h1.Integral(binx1,binx2,biny1,biny2,"")<<std::endl;
std::cout<<"maximum = "<<h1.GetMaximum()<<std::endl;
  c3->cd();
  hnew->Divide(my2dfit);
  hnew->Draw();
  my2dfit->Write();
  hnew->Write();
}