fit gaussiano bidimensionale

ciao a tutti,

oggi durante l'esame mi e' stato chiesto di implmentare un fit gaussiano di un istogramma bidimensionale. Ecco la mia implementazione.

#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TStyle.h"
#include "TProfile.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TF2.h"

#include <iostream>
#include <iomanip>


double randFlat (double min, double max) 
{
  return min + (max-min) * rand () / static_cast<double> (RAND_MAX) ;  
}


// ------------------------------------------------


//! gaussiana 2D senza termini misti
double func2 (double * x, double * par)
{
/*
  par 0 = ampiezza
  par 1 = mean x
  par 2 = sigma x
  par 3 = man y
  par 4 = sigma y
*/
   double argx = 0;   
   if (par[2] != 0) argx = (x[0] - par[1]) / par[2] ;
   double argy = 0;   
   if (par[4] != 0) argy = (x[1] - par[3]) / par[4] ;

   double fitval = par[0] * TMath::Exp (-1 * argx * argx - 1 * argy * argy) ;
   return fitval;
} 


// ------------------------------------------------


int main () 
{
  gROOT->SetStyle ("Plain") ; // sfondo bianco
  gStyle->SetPalette (1) ; // scala di colori "termica"
  TFile salvo ("casuali2.root","RECREATE") ;
  TH2F campana ("campana","campana",50,-5,5,50,-5,5) ;
  TF2 funzione ("gaussiana",func2,-5,5,-5,5,5) ; 
  double parameters[5] = {1,0,1,0,2} ;
  funzione.SetParameters (parameters) ;

  
  for (int i=0 ; i<100000 ; ++i)
    {
      double x,y,z ;
      do
        {
          x = randFlat (-5, 5) ;
          y = randFlat (-5, 5) ;
          z = randFlat (0, 1.5) ;
        }
      while (funzione.Eval (x,y) < z) ;
      campana.Fill (x,y) ;
    }
  salvo.Write () ;  

  double fitParameters[5] = {1,0,1,0,1} ;

  funzione.SetParameters (fitParameters) ;
  campana.Fit ("gaussiana") ;
  funzione.GetParameters (fitParameters) ;
 
  TCanvas c1 ;
  campana.Draw ("LEGO2") ;
  funzione.Draw ("surfsame") ;
  c1.Print ("fit2d.eps","eps") ;
   
  std::cout << std::setprecision (2) ;
  std::cout << "-----------------------------\n" ;
  std::cout << "param number  | gen  |  his  \n" ;
  std::cout << "-----------------------------\n" ;
  for (int i=0 ; i<5 ; ++i)
    {
      std::cout << "parametro " << i
                << "\t" << parameters[i]
                << "\t" << fitParameters[i]
                << "\n" ;
    }
  std::cout << "-----------------------------\n" ;

  return 0 ;
}