/*


  I M P O R T A N T 
Before to execute this macro, you need to load the Class
MultiGauss:

To do it, at the ROOT prompt:
root [1] .L MultiGauss.C
root [2] .x semfit.C

To do the fit you generally need to initialise the parameter you need to get.
This is done as follows (let's suppose that gau is your function of multi gauss)

gau->SetParameter(parametro,valore);

In order to  make the parameter  vary within limits:


gau->SetParLimits(parametro,liminf,limsup);


To fix the parameter, the variability range should be null:
liminf=limsup.

 */




double multigauss(Double_t *x, Double_t *par)
  /*
    Funzione somma di piu' gaussiane.
   */
{

  TVector *Amp=MultiGauss::Ampli();
  TVector *Media=MultiGauss::Mean();
  TVector *Erro=MultiGauss::Sigma();

  int ncomp = Amp->GetNrows();
  int i=0;
  double funcval=0;
  while(i<ncomp)
    {

      //	double ampli=Amp(i);// uncomment this if you don't want to fity, but simply plot the multigauss
      	double ampli=par[i];
	double mean=Media(i);
	double sigma=Erro(i);


	funcval += ampli
	  *TMath::Exp(
		      -TMath::pow(x[0]-mean,2)
		      /(2*TMath::pow(sigma,2))
		      );
	i++;
    }

  return funcval;

  // end multigauss
}



mgaufit()
{
  /// Let's assume we want to fit the gaussians' amplitude, knowing their mean and sigma

  TFile *fi = new TFile("Spoli.root"); // I read the file with the histogram
 
  TH1F *histo = fhh1;
  histo->SetMaximum(700);
  histo->Draw();
  

  int nog =12;
  //              0   1   2   3   4   5   6   7   8   9  10  11

  //  Range of variation for the amplitude
  double iamp[]={ 0, 10,160, 40,170,220, 60, 60, 10, 10, 10, 10};
  double amp[] ={ 0, 60,170,160,175,250,166,155,140, 60, 45,  5};
  double samp[]={ 0,150,200,250,400,350,200,200,150,100, 60, 10};

  TVector *A = new TVector(nog);
  TVector *M = new TVector(nog);
  TVector *S = new TVector(nog);

  TF1 *gau= new TF1("multigauss",multigauss,0,1000,12); // 12 is the number of parameters to fit (amplitudes of 12 gaussians )

  for(int ui=0;ui<nog;ui++)
    {
      // this is my method to get the mean and the sigma... you shoulkd implement yours
      double pval = Deco::p(5+ui);
      double sigpval = Deco::sigp(pval);
      M(ui)=pval; S(ui)=sigpval;
    
      // I set the amplitude
      gau->SetParameter(ui,amp[ui]);
      gau->SetParLimits(ui,iamp[ui],samp[ui]);

    }


  // If we want to fix some amplitudes, letting the others to vary:
  gau->FixParameter(0,0);//5
  gau->FixParameter(1,0);//6
//  gau->FixParameter(2,160);//7
//  gau->FixParameter(3,40);//8
//  gau->FixParameter(4,220);//9
//  gau->FixParameter(5,120);//10
//  gau->FixParameter(6,100);//11
//  gau->FixParameter(7,134);//12
//  gau->FixParameter(8,20);//13
//  gau->FixParameter(9,47);//14
//  gau->FixParameter(10,60);//15
//  gau->FixParameter(11,5);//16



    /* 
       I set in memory the values of the sigmas and means

       memorizzo i valori dei TVectors invocando un metodo di una classe senza definire
       alcun oggetto di essa, quindi per cosi' dire, agendo dall'esterno della stessa.
     */

  //  MultiGauss::SetAmpli(A);
  MultiGauss::SetMean(M);
  MultiGauss::SetSigma(S);



  // BEWARE: it is fundamental the  name of the TF1 ("multigauss ") to be the same of the related function (multigauss(){...}).

    double cchi2=0;

    double supbinval=histo->GetBinCenter(runmax+1);

    histo->Fit(multigauss,"","",1,supbinval);
      
    
    
}