#include<iostream>
#include "TH1D.h"
#include "TF1.h"
#include "TFile.h"
#include "TCanvas.h"

double total(double *x, double *par)
{
	double xx = x[0];
	double func = par[0] * TMath::Gaus(xx, par[1], par[2]) + par[3] + par[4]*xx + par[5]*xx*xx + par[6]*xx*xx*xx; 
	return func;
}	
	
double peak(double *x, double *par)
{
	double xx = x[0];
	double func = par[0] * TMath::Gaus(xx, par[1], par[2]);
	return func;
}

double background(double *x, double *par)
{
	double xx = x[0];
	double func = par[0] + par[1]*xx + par[2]*xx*xx + par[3]*xx*xx*xx;
	return func;
}	


void testFit()
{
	TFile *f = new TFile("test.root");
	TH1D *h1 = (TH1D*)f->Get("h");

	TF1 *total = new TF1("total", total, 1.12, 1.27, 7);
	total->SetParameter(0, 500);
	total->FixParameter(1, 1.189);
	total->SetParameter(2, 0.006);
	total->SetParameter(3, 1.0);
	total->SetParameter(4, 1.0);
	total->SetParameter(5, 1.0);
	total->SetParameter(6, 500);
	total->SetLineColor(2);
	h1->Fit(total, "RMEQB");
	
	TF1 *peak = new TF1("peak", peak, 1.169, 1.209, 3);
	peak->SetParameter(0, total->GetParameter(0));
	peak->SetParameter(1, total->GetParameter(1));
	peak->SetParameter(2, total->GetParameter(2));
	peak->SetParError(0, total->GetParError(0));
	peak->SetParError(1, total->GetParError(1));
	peak->SetParError(2, total->GetParError(2));
	peak->SetLineColor(4);
	
	TF1 *background = new TF1("background", background, 1.12, 1.27, 4);
	background->SetParameter(0, total->GetParameter(3));
	background->SetParameter(1, total->GetParameter(4));
	background->SetParameter(2, total->GetParameter(5));
	background->SetParameter(3, total->GetParameter(6));
	background->SetParError(0, total->GetParameter(3));
	background->SetParError(1, total->GetParameter(4));
	background->SetParError(2, total->GetParameter(5));
	background->SetParError(3, total->GetParameter(6));
	background->SetLineColor(6);
	 
	TCanvas *c = new TCanvas("c", "c", 800, 650);
	c->cd();
	h1->Draw();
	peak->Draw("same");
	background->Draw("same");
		
	Int_t np = 10000;
        double *x = new double[np];
        double *w = new double[np];
        peak->CalcGaussLegendreSamplingPoints(np, x, w, 1e-15);

	double fBinSize = h1->GetBinWidth(2);
	double xImin = 1.174;
	double xImax = 1.204;

	double numPart = peak->IntegralFast(np, x, w, xImin, xImax) / fBinSize;
	double error = total->IntegralError(xImin, xImax, 1e-6);	

	std::cout << "Integrating between " << xImin << " - " << xImax << std::endl;
        std::cout << numPart << "    " << error << std::endl;
	delete [] x;
	delete [] w;
	
//	f->Close();
}