// This program performs a gaussian plus gaussian fit to a histogram
// 
// Author: Sean Kuvin
// Date: 3-18-08
//

#include "TApplication.h"
#include "TH1.h"
#include "TF1.h"
#include "TMath.h"
#include <iostream>

using namespace std;

//Gaussian fit

double Gaussian(double *x, double *par){
  return TMath::Exp(-0.5*(x[0]-par[2])*(x[0]-par[2]) /
    TMath::Max(1e-10,(par[1]*par[1])))*par[0] /
      TMath::Max(1e-10,(par[1]*TMath::Sqrt(2*3.14)));
}


double fitFunction(double *x,double *par){
  return Gaussian(x,par)+Gaussian(x,&par[3]);
}


void doFit(){
  const int nbins = 60;
  double data[60]={  6, 1,10,12, 6,13,23,22,15,21,
		    23,26,36,25,27,35,40,44,66,81,
		    75,57,48,45,46,41,35,36,53,32,
		    40,37,38,31,36,44,42,37,32,32,
		    43,44,35,33,33,39,29,41,32,44,
	            26,39,29,35,32,21,21,15,25,15};

  // Look up TH1F at 
  //    http://root.cern.ch/root/Reference.html
  TH1F *histogram = new TH1F("histo", "title", 60, 0, 3);

  for(int i =0; i<nbins;i++) {
    histogram->SetBinContent(i+1, data[i]);
    histogram->SetBinError(i+1, TMath::Sqrt(data[i]));
  }

  // Look up TF1 at 
  //	http://root.cern.ch/root/Reference,html
  TF1* fitFcn = new TF1("fitFcn", fitFunction,0,3,6);
  fitFcn->SetParameter(4,0.2);
  fitFcn->SetParameter(5,1);
  fitFcn->SetParameter(0,0.1);
  fitFcn->SetParameter(1,0.1);
  fitFcn->SetParameter(2,.5);
  fitFcn->SetParameter(3,10);

  // try fitting w/ starting values for the parameter
  histogram->Fit("fitFcn","V+");

  TF1*gaussFcn=new TF1("gaussFcn", Gaussian, 0,3,3);
  TF1*gaussFcn2=new TF1("gaussFcn2",Gaussian,0,3,3);

  gaussFcn->SetLineColor(3);
  gaussFcn2->SetLineColor(4);

  double par[6];
  fitFcn->GetParameters(par);
  gaussFcn->SetParameters(par);
  gaussFcn->Draw("Same");
  gaussFcn2->SetParameters(&par[3]);
  gaussFcn2->Draw("Same");
}

int main(int argc,char **argv) {
  TApplication theApp("App",&argc,argv);
  doFit();
  theApp.Run();
}