Double_t gauss2d( Double_t *x, Double_t *par)
{
   Double_t sqrt2pi = sqrt(2*TMath::Pi()), sqrt2 = sqrt(2.);
   Double_t cx = par[0];    // const x
   Double_t mx = par[1];    // mean x
   if (par[2] == 0) par[2] = 1;
   Double_t sx = par[2];    // sigma x
   Double_t cy = par[3];    // const y
   Double_t my = par[4];    // mean y
   if (par[5] == 0) par[5] = 1;
   Double_t sy = par[5];    // sigma 
   Double_t xij = x[0] - mx;  
   Double_t argx = xij/(sqrt2*sx);
   Double_t xval = cx * exp(-argx * argx)/(sqrt2pi*sx);
   Double_t yij = x[1] - my;  
   Double_t argy = yij/(sqrt2*sy);
   Double_t yval = cy * exp(-argy * argy)/(sqrt2pi*sy);
   return xval * yval;
};

//fit_user_function(const char *hname)
fit_user_function()
// This is a template macro to fit a user defined 2d function
// As an example a 2d gauss is provided
{
//   gROOT->Reset();
   TCanvas * c1 = new TCanvas("c1","2D Example",200,10,600,600);
   TH2* hist  = new TH2F("hpxpy1","py vs px",40,-4,4,40,-8,8);
   gRandom->SetSeed();
   Float_t px, py, pz;
   for (Int_t i = 0; i < 25000; i++) {
     gRandom->Rannor(px,py);
     hpxpy1->Fill(0.5*px, 1.5*py);
	}
	hist->SetOption("lego2");
	hist->Draw();
//   TH2* hist = (TH1*)gROOT->FindObject(hname);
//   if(!hist){
//     cout << "histogram not found" << endl;
//     return 0;
//   }
//  try to use same canvas dimensions for fitted function as used for histogram 
   TCanvas * mc  = (TCanvas *) gPad;
   Int_t wtopx = 100, wtopy = 100;
   UInt_t ww = 700, wh = 700;
   TString cname("canvas_2d");
   if (mc) {
      mc->GetCanvasPar(wtopx, wtopy, ww, wh);
      cname = mc->GetName();
      cname += "_2df";
      wtopx = TMath::Max(0,wtopx - 100);
      wtopy = TMath::Max((Int_t)wh, wtopy + 100);
//       cout << wtopx << " "  << wtopy << " "  <<  ww << " " <<  wh << endl;
   }

//  nb: fit range may be different from histo range
   Float_t x_from = -4; 
   Float_t x_to   =  4; 
   Float_t y_from = -7; 
   Float_t y_to   =  7; 
   TF2 * f2 = new TF2("gauss2d", gauss2d, x_from, x_to, y_from, y_to , 6);
   f2->SetParName(0, "const_x");
   f2->SetParName(1, "mean_x");
   f2->SetParName(2, "sigma_x");
   f2->SetParName(3, "const_y");
   f2->SetParName(4, "mean_y");
   f2->SetParName(5, "sigma_y");

   f2->SetParameters(100, 0, 1, 100, 0, 1);
 //  f2->Draw();
   hist->Fit("gauss2d","R0","same");  // dont draw in same picture 
   TCanvas *cc = new TCanvas(cname , cname, wtopx, wtopy, ww, wh);
//  draw in same range as histogram with same drawing options
   f2->SetRange(hist->GetXaxis()->GetXmin(), hist->GetYaxis()->GetXmin(),
                hist->GetXaxis()->GetXmax(), hist->GetYaxis()->GetXmax()); 
   f2->SetNpx(hist->GetNbinsX());
   f2->SetNpy(hist->GetNbinsY());
   f2->Draw(hist->GetOption()); 
}