#include <TROOT.h>
#include <TH1F.h>
#include <TF1.h>
#include <TMath.h>
#include <TStyle.h>
#include <TRandom.h>
#include <TVirtualFitter.h>
#include <TFitResult.h>

#include <iostream>
#include <fstream>
#include <vector>
#include <set>
#include <map>

float BIN_SIZE = 0;

inline Double_t my_gauss(Double_t * x, Double_t * par)
{
  // par[0]: # of events, par[1]: mean, par[2]: sigma
  Double_t arg = 0;
  if (par[2]<0) par[2]=-par[2]; // make sure sigma > 0
  if (par[2] != 0) arg = (x[0] - par[1])/par[2]; 
  return par[0]*BIN_SIZE*TMath::Exp(-0.5*arg*arg)/
    (TMath::Sqrt(2*TMath::Pi())*par[2]);
}

using namespace std;

int testIntegralError(void)
{
  gStyle->Reset();
  // mean of gaussian
  const float mean = 10;
  // upper, lower limits of distribution
  const float LOW_LIMIT = mean -4;
  const float HI_LIMIT = mean + 4;
  const float fit_limit = mean;

  // this histogram is the "data"
  TH1F * my_hist = new TH1F("my_hist", "Gauss distribution", 200, 
			    LOW_LIMIT, HI_LIMIT);
  BIN_SIZE = my_hist->GetBinWidth(0);


  // This function will be used for "sampling" values into a histogram
  TF1 * fgen = new TF1("gauss1", my_gauss, LOW_LIMIT, HI_LIMIT, 3);
  fgen->SetParameters(1, mean, 1); // # of events, mean, sigma

  const int NENTRIES = 15000;

  // this function will be used to fit the histogram
  TF1 * fhist = new TF1("gauss2", my_gauss, LOW_LIMIT, fit_limit, 3);
  fhist->SetParNames("Evt count", "Mean", "Sigma");
  // initial "guess" for the parameters
  fhist->SetParameters(NENTRIES, mean, 0.8); // # of events, mean, sigma
  fhist->SetLineColor(kRed);

  gStyle->SetOptFit(1111);
  gStyle->SetStatH(0.0);
  gStyle->SetStatW(0.0);

#if 1
  fhist->SetParLimits(1, mean, mean);
#endif

  for(int i = 0; i != NENTRIES; ++i)
    my_hist->Fill(fgen->GetRandom());
  
  TFitResultPtr  r = my_hist->Fit(fhist, "ISLR");
  
  TVirtualFitter * fitter = TVirtualFitter::Fitter(0,3);
  Double_t par[3] = {fitter->GetParameter(0), fitter->GetParameter(1), 
		     fitter->GetParameter(2)};

  cout << " Histogram integral: " << my_hist->Integral() << endl;
  cout << " Fit-function integral between " << LOW_LIMIT << " and " 
       << fit_limit << ": " << fhist->Integral(LOW_LIMIT, fit_limit)
    /BIN_SIZE
       << " +- " << fhist->IntegralError(LOW_LIMIT, fit_limit)/BIN_SIZE
       << endl;
  cout << " Fit-function integral between " << fit_limit << " and " 
       << HI_LIMIT << ": " 
       << fhist->Integral(fit_limit, HI_LIMIT)/BIN_SIZE
       << " +- " 
       << fhist->IntegralError(fit_limit, HI_LIMIT, par,
			       r->GetCovarianceMatrix().GetMatrixArray())/BIN_SIZE 
       << endl;

  return 0;
}