//Programma di analisi dei file .root prodotti dai TLG

//#include <stdlib.h>
#include <TMath.h>
#include "TCanvas.h"
#include "TH1.h"
#include "TF1.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"

Double_t Background(Double_t *x, Double_t *par);
Double_t Signal(Double_t *x, Double_t *par);
Double_t Total(Double_t *x, Double_t *par); 

main(char *name)
{
	gROOT->Reset();

	Int_t npeaks;

    TFile *file = new TFile(name);

	MyWinSpectrumPeaks = new TCanvas("MyWinSpectrumPeaks", "spectrum peaks", 200, 10, 900, 700);

 
	//TH1F *h_slope = new TH1F("h_slope","Peaks slopes",100,-0.7,0.7);

	TH1F *h_slope = (TH1F*)file->Get("h_slope;1");

	MyWinSpectrumPeaks->Divide(1,2);

	MyWinSpectrumPeaks->cd(1);
	h_slope->Draw();

	//Use TSpectrum to find the peak candidates
	TSpectrum *spc = new TSpectrum(2*npeaks);
    Int_t nfound = spc->Search(h_slope,1,"new");

    printf("Found %d candidate peaks to fit\n",nfound);

	MyWinSpectrumPeaks->Update();
    MyWinSpectrumPeaks->cd(2);

	//signal
	Double_t *par;
	par = new Double_t[3+3*nfound];

	npeaks = nfound;

	//Get peaks position and set initial parameters:
	Float_t *xpeaks = spc->GetPositionX();
	
	for(Int_t i=0; i<nfound; i++)
	{
		printf("%f \n", xpeaks[i]);
	}
	
	
    for(Int_t p=0;p<nfound;p++)
    {
	  printf("Fit picco: %d \n", p);
      
	  Float_t xp = xpeaks[p];

	  printf("Peaks position %d is: %f \n",p,xp);

      Int_t bin = h_slope->GetXaxis()->FindBin(xp);

	  printf("Global bin number corresponding to %f is: %d \n", xp,bin);

      Float_t yp = h_slope->GetBinContent(bin);

	  printf("Content of bin number %f is: %d \n", xp,yp);

	  par[3*p+3] = yp;
      par[3*p+4] = xp;
      par[3*p+5] = 3;

	  printf("Fit parameters: %f %f %f \n", par[3*p+3],par[3*p+4],par[3*p+5]);

	}

	//Fit background using a gaussian
	TF1 *fbackg = new TF1("fbackg","gaus",-0.7,0.7);

	//h_slope->Fit("fbackg","q0");
	h_slope->Fit("fbackg");
	
    par[0] = fbackg->GetParameter(0);
    par[1] = fbackg->GetParameter(1);
	par[2] = fbackg->GetParameter(2);

	printf("/////////////////////////////////////////////////// \n");
	printf("Background fit parameters: %f %f %f \n", par[0],par[1],par[2]);
	printf("/////////////////////////////////////////////////// \n");

    	
	//Subtract estimated background to signal
	//Creates a temporary histogram
	TH1F *h_temp = (TH1F*)h_slope->Clone();
	h_temp->Reset();

	TF1 *eback = h_slope->GetFunction("fbackg");

	for(Int_t bin = 1;bin <= 100; bin++)
	{
		Float_t x = h_slope->GetBinCenter(bin);

		Double_t fval = eback->Eval(x);

		Double_t diff = TMath::Abs(fval - h_slope->GetBinContent(bin));

		h_temp->Fill(x,diff);

	}


	//Fit signal with my function
	TF1 *fsignal = new TF1("fsignal",Signal,-0.7,0.7,3*npeaks);

	h_temp->GetListOfFunctions()->Add(fsignal);
	h_temp->Fit("fsignal");

	//TF1 *esig = h_temp->GetFunction("gaus");
	//TF1 *esig = h_temp->GetFunction("fsignal");

	//Fit background + signal
	printf("Now be fitting, be patient \n");

	eback->GetParameters(&par[0]);
	//esig->GetParameters(&par[3]);
	fsignal->GetParameters(&par[3]);

	for(Int_t m=0; m<3+3*nfound; m++)
	{
		printf("Check parameters: %f \n", par[m]);

	}

	//Make the final function for fit:
	TF1 *theory = new TF1("theory",Total,-0.7,0.7,3+3*nfound);
	theory->SetParameters(par);
	
	//TVirtualFitter::Fitter(h_slope,10+3*npeaks); //we may have more than the default 25 parameters
	//TVirtualFitter::Fitter(h_slope)->SetFCN("theory");
	//theory->SetParameters(par);
    //theory->SetNpx(1000);


	//gStyle->SetOptFit(1111); 

	//h_slope->GetListOfFunctions()->Add(theory); 

	h_slope->Fit("theory");


  
	
}

Double_t Background(Double_t *x, Double_t *par) 
{
   //The background function
   Double_t arg = 0;

   if(par[2]) arg = (x[0] - par[1])/par[2];

   Double_t val = par[0]*TMath::Exp(-0.5*arg*arg)*x[0]*x[0];

   return val;
}

Double_t Signal(Double_t *x, Double_t *par) 
{
   //The signal function

   Double_t result = 0;
   Double_t e = 0;
  
   Int_t npeaks = 6;

   for(Int_t p=0;p<npeaks;p++)
   {
      Double_t norm  = par[3*p+3];
      Double_t mean  = par[3*p+4];
      Double_t sigma = par[3*p+5];

	  if(sigma != 0) e = (x[0] - mean)/sigma;

	  result += norm*TMath::Exp(-0.5*e*e);
     
   }
   
   return result;
}

Double_t Total(Double_t *x, Double_t *par)
{
	//Combined background + signal

	Double_t tot = Background(x,par) + Signal(x,&par[3]);

	return tot;

}