Fitting functions
The next exercise uses fitting functions. There are several ways to make functions in ROOT1; this will describe one of the most general ways to do so. Others are listed in the ROOT docs.
The fitting function
For a general function, ROOT will call it with two pointers to arrays; the first is an array of the independent variables, and the second is an array of the parameters. You will have to make sure your array lengths are correct when you are setting up your ROOT code; these are simple C arrays and as such do not tell you the number of elements they contain.
Here is an example function for a linear background plus a Gaussian signal. It evaluates a Gaussian signal on top of a linear background. The signal has the form
Arguments 0, 1, and 2 are
total
, , and . The function returns the
estimated function value per unit bin width. This particular
fit function assumes the bin width is 1 MeV, but there is a
place-hold to multiply fitval
by the bin-width to get the
normalization right. To understand this, make a histogram with, for
example, 2 MeV wide bins. Think about how the height of each bin
changes for the same data set, how the total relates to the height
of the central bin. Arguments 3 and 4 describe a linear background.
This file must be called fit1MeV_Gaussian.C
.
Double_t fit1MeV_Gaussian(Double_t *v, Double_t *par) {
// Give nice names to values in the arrays
// Using & to avoid a value copy
Double_t &x = v[0];
Double_t &par0 = par[0];
Double_t &mean = par[1];
Double_t &width = par[2];
Double_t &lin_b = par[3];
Double_t &lin_m = par[4];
// Arg is 0 if width is 0, otherwise calculate it
Double_t arg = width == 0 ? 0 : (x - mean)/width;
// Create a Gaussian
Double_t fitval = par0*TMath::Exp(-0.5*arg*arg);
// Re-normalize to make par0 the integral of the Gaussian, assuming 1 MeV bins
fitval = (fitval*1.00)/(TMath::Sqrt(TMath::TwoPi())*width);
// Add a linear background
fitval = fitval + lin_b + x*lin_m;
return fitval;
}
Preparing the TF1 object
Your goal is to make a TF1
1-dimensional
function class object for ROOT to fit with. Using
the previous function, that would look like:
TF1 *myXiFit = new TF1("myXiFit", fit1MeV_Gaussian, 1303., 1340., 5);
myXiFit->SetParameter(0, 100000.);
myXiFit->SetParameter(1, 1322.);
myXiFit->SetParameter(2, 3.);
myXiFit->SetParLimits(2, 0., 10.);
myXiFit->SetParameter(3, 500);
myXiFit->SetParameter(4, 0.);
The first line of this code creates a pointer to a TF1
object.
- The first argument of
new TF1(...)
is the "name"myXiFit
which will be used later. - The second argument,
fit1MeV_Gaussian
, is the name of an external function used to describe the data in the histogram. - The third and fourth arguments are the lower and upper boundaries of the fitting range.
- The fifth argument is the number of parameters used by the fitting function.
The next lines provide initial values for the function parameters and limits for the Gaussian width (parameter 2). By specifying that this parameter must be positive, the code avoids the mathematical ambiguity associated with using only the square of the width in the fitting function.
Bonus: Functional
If you want to make a general function that takes any width, you can use a C++ pattern called a functional. This file must be called fit_Gaussian.C
struct fit_Gaussian {
Double_t bin_width;
/// Constructor initializes bin_width
fit_Gaussian(Double_t bin_width_init) : bin_width(bin_width_init) {}
/// This is called when you call the object
Double_t operator() (Double_t *v, Double_t *par) {
Double_t arg = width == 0 ? 0 : (v[0] - par[1])/par[2];
Double_t fitval = par[0]*TMath::Exp(-0.5*arg*arg);
fitval = (fitval*bin_width)/(TMath::Sqrt(TMath::TwoPi())*par[2]);
fitval += par[3] + v[0]*par[4];
return fitval;
}
};
You would use it like the other one after you create an instance that has your bin width in it:
.L fit_Gaussian.C
auto fit1MeV_Gaussian = fit_Gaussian(1.0);
1. There are several ways to do anything in ROOT. ↩