using namespace std;
#include <iostream>
#include <fstream>
//ROOT
#include "TApplication.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TH1.h"

Double_t poly2(Double_t *, Double_t *);
double computeFWHM(TH1D *);
//int main(int argc, char **argv) {
int testError() {
   	int iNBins = 100;

	TH1D *histo = new TH1D("h1","Current Distribution",iNBins,-2,2);
	cout << "Histogram created" << endl;
	histo->FillRandom("gaus",100000);
	cout << "Histogram filled" << endl;
	//on normalize la distribution
        histo->Scale(1.0/histo->Integral("width"));
	cout << "Histogram scaled" << endl;
	
	cout << "Largeur a mi-hauteur : " << computeFWHM(histo) << endl;
	
	//delete histo;
}

/******************************************************************************
 * computeFWHM
 * Cette fonction calcule la largeur a mi-hauteur (Full Width at Half Maximum)
 * a partir d'un histogramme. Elle assume qu'il n'y a qu'un seul maximum
 * dans la courbe.
 *
 * Elle commence par trouver le maximum de l'histogramme, puis elle trouve la
 * premiere case dont la valeur se situe en-dessous de max/2 de chaque cote du
 * maximum. Ensuite, elle fait un fit parabolic sur les 2 points de chaque
 * cote pour determiner la position reelle de la mi-hauteur.
 *
 * Entree :
 * 	Un pointeur vers l'histogramme
 * Sortie :
 * 	double : la largeur a mi-hauteur de l'histogramme.
 * ***************************************************************************/
double computeFWHM(TH1D * histo) {
	//on trouve le maximum
	Double_t max = histo->GetMaximum();
	Int_t maxbin = histo->GetMaximumBin();
	//Double_t xmax = histo->GetBinCenter(maxbin);
	
	//on trouve le point milieu du cote gauche
	histo->GetXaxis()->SetRange(0,maxbin);
	Double_t max_2_L = histo->GetMaximum(max/2.0);
	Int_t maxbin_2_L=0;
	Double_t xmax_2_L;
	for (int i=0; i<histo->GetNbinsX() && i<maxbin; i++) {
		if (histo->GetBinContent(i) == max_2_L)
			{	maxbin_2_L = i; 	break;		}
	}
	xmax_2_L = histo->GetBinCenter(maxbin_2_L);
	
	Double_t a,b,c;
	//on fait un fit quadratique du cote gauche
	TF1 * poly2L = new TF1("Poly2L",poly2,0,histo->GetBinCenter(histo->GetNbinsX()),3);
	Double_t penteL = (histo->GetBinContent(maxbin_2_L)-histo->GetBinContent(maxbin_2_L+1)) / (xmax_2_L - histo->GetBinCenter(maxbin_2_L+1));
	poly2L->SetParameters(0.0,penteL,max_2_L);
	histo->Fit("Poly2L","","",histo->GetBinCenter(maxbin_2_L-2),histo->GetBinCenter(maxbin_2_L+2));
	//on calcule le xmax_2_L reel
	a = poly2L->GetParameter(0);
	b = poly2L->GetParameter(1);
	c = poly2L->GetParameter(2) - max/2.0;
	Double_t xmax_2_LR = (-b + sqrt(b*b-4*a*c))/(2*a);
        //delete poly2L;

	//on trouve le point milieu du cote droit
	histo->GetXaxis()->SetRange(maxbin,histo->GetNbinsX());
	Double_t max_2_R = histo->GetMaximum(max/2.0);
	Int_t maxbin_2_R=0;
	Double_t xmax_2_R;
	for (int i=maxbin; i<histo->GetNbinsX(); i++) {
		if (histo->GetBinContent(i) == max_2_R)
			{	maxbin_2_R = i; 	break;		}
	}
	xmax_2_R = histo->GetBinCenter(maxbin_2_R);
	
	//on fait un fit quadratique du cote droit
	TF1 * poly2R = new TF1("Poly2R",poly2,0,histo->GetBinCenter(histo->GetNbinsX()),3);
	Double_t penteR = (histo->GetBinContent(maxbin_2_R)-histo->GetBinContent(maxbin_2_R+1)) / (xmax_2_R - histo->GetBinCenter(maxbin_2_R+1));
	poly2R->SetParameters(0.0,penteR,max_2_R);
	histo->Fit("Poly2R","+","",histo->GetBinCenter(maxbin_2_R-2),histo->GetBinCenter(maxbin_2_R+2));
	//on calcule le xmax_2_R reel
	a = poly2R->GetParameter(0);
	b = poly2R->GetParameter(1);
	c = poly2R->GetParameter(2) - max/2.0;
	Double_t xmax_2_RR = (-b - sqrt(b*b-4*a*c))/(2*a);
	//delete poly2R;
	
	histo->GetXaxis()->SetRange(0,histo->GetNbinsX());
	return xmax_2_RR - xmax_2_LR;
}

/******************************************************************************
 * poly2
 * Cette fonction calcule un polynome de degre 2 dans l'optique de faire
 * un fit.
 * Entree :
 * 	*x : Un pointeur vers un Double_t x
 * 	*p : Un tableau de 3 parametres : y = p[0]*x^2+p[1]*x+p[2]
 * Sortie :
 * 	Double_t y = p[0]*x^2+p[1]*x+p[2]
 * ***************************************************************************/
Double_t poly2(Double_t *x, Double_t *p) {
	Double_t a = p[0];
	Double_t b = p[1];
	Double_t c = p[2];

	return a*x[0]*x[0] + b*(x[0]) + c;
}