#include "TF2.h" #include "TH2.h" #include "TMath.h" // Fitting a 2-D histogram // This tutorial illustrates : // - how to create a 2-d function // - fill a 2-d histogram randomly from this function // - fit the histogram // - display the fitted function on top of the histogram // // This example can be executed via the interpreter or ACLIC // root > .x fit2.C // root > .x fit2.C++ //Author: Rene Brun Double_t g2(Double_t *x, Double_t *par) { Double_t r1 = Double_t((x[0]-par[1])/par[2]); Double_t r2 = Double_t((x[1]-par[3])/par[4]); return par[0]*TMath::Exp(-0.5*(r1*r1+r2*r2)); } //This is a gaussian of the form f(x,y) = A^(-([(x-x[0])^2/2*sigmax^2]+[(y-y[0])^2/2*sigmsy^2]). You can see this format in the equation given by return above. //where A = amplitude (x[0],y[0]) is the center, and sigmax and sigmay are the spreads of the "blob" //In general a 2D elliptical Gaussian has the form f(x,y) = A^(-(a(x-x[0])^2+2b(x-x[0])(y-y[0])+c(y=y[0])^2)) where matrix [a b b c] is positive definite Double_t fun2(Double_t *x, Double_t *par) { Double_t *p1 = &par[0]; Double_t *p2 = &par[5]; Double_t *p3 = &par[10]; Double_t result = g2(x,p1) + g2(x,p2) + g2(x,p3); //this adds the three parts of the function together to create three Gaussians return result; } void fit2() { const Int_t npar = 15; Double_t f2params[npar] = {100,-3,3,3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7}; //I believe these are three sets of parameters for the three humps on the graph. //format: amplitude, x[0], sigmax, y[0], sigmay //first set: big hump at -3,-3 second set: sharp peak at 0,0 third set: small hill at 4,4 TF2 *f2 = new TF2("f2",fun2,-10,10,-10,10, npar); f2->SetParameters(f2params); //Create an histogram and fill it randomly with f2 TH2F *h2 = new TH2F("h2","from f2",40,-10,10,40,-10,10); Int_t nentries = 100000; h2->FillRandom("f2",nentries); //Fit h2 with original function f2 //This part is redefining parameters [0], [5], and [10]; in other words, the amplitudes. It appears to be scaling them by the number of entries. It has the most effect on the small hill. //to be entirely honest, I see no effect at all when I change it. However, removing it entirely screws everything up. Float_t ratio = 4*nentries/100000; std::cout<<"ratio = "<SetParameters(f2params); h2->GetXaxis()->SetTitle("X Axis"); h2->GetYaxis()->SetTitle("Y Axis"); h2->Fit("f2"); h2->Draw("SURF2"); f2->Draw("SURF same"); cout<<"Inegral of f2: "<Integral(-10,10,-10,10)<