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.

results matching ""

    No results matching ""