void FitTESTE(){
    TFile *ficheiro = new TFile("histogramTESTE.root","READ");
  //  TFile *ficheiroLight = new TFile("HistogramasFitted2Gauss.root","RECREATE");
    Double_t minX, maxX, ErrorSigma, minimumX, maximumX, min2, max2, min3, max3, media, sumSST, sumSSE, R2;
    TH1D *histogramas; //Declare an array of histogram pointers
	Int_t BinValue, dummyNBins;

    TCanvas *Cajustes = new TCanvas("CanvasAjustes");


     energias=35;
	TString hname="histoMediaPosX35";
        TKey *key = ficheiro->FindKey(hname); //Ve se encontra os histogramas com aquele nome

        if (key ==0){
            cout << "!!Histogram "<< hname <<" does not exist!!" << endl;
            throw 1;
        }

        histogramas =  (TH1D*)ficheiro->Get(hname);


        histogramas->Draw();

        minX=histogramas->GetXaxis()->GetXmin();
        cout<<"minimum "<<minX<<endl;
        maxX=histogramas->GetXaxis()->GetXmax();
        cout<<"maximum "<<maxX<<endl;

//Fit 3 single gaussians and use their output as starting point for a function composed by e gaussians
           	minX=-2000;
            maxX=2000;
            TF1 *g1    = new TF1("g1","gaus",minX,maxX);
            g1->SetLineColor(1);
            min2=-100;
            max2=100;
            TF1 *g2    = new TF1("g2","gaus",min2,max2);
            g2->SetLineColor(3);
			min3=-10;
            max3=10;

			TF1 *g3    = new TF1("g3","gaus",min3,max3);
            g3->SetLineColor(5);


            histogramas->Fit(g1,"Rq");
//            histogramas->GetFunction("g1")->SetBit(TF1::kNotDraw);
            histogramas->Fit(g2,"Rq+");
//            histogramas->GetFunction("g2")->SetBit(TF1::kNotDraw);
            histogramas->Fit(g3,"Rq+");
//            histogramas->GetFunction("g3")->SetBit(TF1::kNotDraw);
            
            
            Double_t par3[9];
            g1->GetParameters(&par3[0]);
            g2->GetParameters(&par3[3]);
			g3->GetParameters(&par3[6]);
			
            TF1 *fitGauss = new TF1("fitGauss","gaus(0)+gaus(3)+gaus(6)",minX,maxX);

            fitGauss->SetParameters(par3);
            histogramas->Fit("fitGauss","R+");


   
/************************R2 CALCULATION****************************************/
// I have tried to follow the definition present here: http://www.mathworks.com/help/curvefit/evaluating-goodness-of-fit.html
//First copy the histogram ROI
Int_t binmin=histogramas->GetXaxis()->FindBin(minX);//Find bin index in wich is the minimim value
Int_t binmax=histogramas->GetXaxis()->FindBin(maxX);
cout<<"min "<<binmin<< " max "<<binmax<<endl;

dummyNBins=(binmax-binmin);

TH1D *dummyHisto=new TH1D("dummyHisto","dummyHisto",dummyNBins,minX,maxX);
Int_t bino=0;
TVectorD Yvalues(dummyNBins);
TVectorD Yexpected(dummyNBins);
TAxis *xaxis = histogramas->GetXaxis();

for (Int_t bin=binmin;bin<binmax;bin++){

	BinValue=histogramas->GetBinContent(bin);
	dummyHisto->SetBinContent(bino,BinValue);
	Yvalues[bino]=BinValue;
	Double_t binCenter = xaxis->GetBinCenter(bin);
	Yexpected[bino]=fitGauss(binCenter);
	bino=bino+1;}

//Average calcullation
double integral = dummyHisto->Integral("width");
Int_t Nbins=dummyHisto->GetNbinsX();
cout<<"Integral: "<<integral<<endl;
cout<<"NBins: "<<Nbins<<endl;
media=integral/Nbins;
cout<<"Average: "<<media<<endl;
sumSST=0;
sumSSE=0;
for (int k=0;k<dummyNBins;k++){
sumSSE=sumSSE+((Yvalues[k]-Yexpected[k])**2);
sumSST=sumSST+((Yvalues[k]-media)**2);
}
cout<<"sumSSE "<<sumSSE<<endl;
cout<<"sumSST "<<sumSST<<endl;
R2=1-(sumSSE/sumSST);
cout<<"R2 "<<R2<<endl;

		
			
/******************************************************************************/



}//void