11.4. Custom Fit Functions
Table of Contents
Defining your own model
In many practical analyses you will need fit functions that are not covered by ROOT built ins such as Gaussian or polynomial functions. ROOT lets you define completely custom models and use them with the same fitting machinery as standard functions.
There are two main ways to define a custom fit model with TF1:
- As a mathematical expression string.
- As a C++ function that you write yourself.
Both approaches use parameters $[0], [1], \dots$ inside the model. During a fit, ROOT changes these parameters to minimize the chi square (for binned fits) or the negative log likelihood (for unbinned fits) and find the best description of your data.
In a custom fit function the free parameters are always referred to by index as $[0], [1], [2], \dots$ inside expression strings or by the array par[] inside C++ functions.
Expression based custom functions
The quickest way to define a custom model is to pass a formula string to TF1. The string can contain variables and parameters. The variable is always called x and parameters are written as [0], [1] and so on.
For example, a simple exponential with an offset can be written as
$$
f(x) = [0] \, e^{-[1] x} + [2].
$$
In ROOT C++ you can create this function as:
TF1 *fexp = new TF1("fexp", "[0]*exp(-[1]*x) + [2]", 0.0, 10.0);
The fourth and fifth arguments specify the valid range in x. Outside this range the function is not evaluated for fitting or drawing.
You can combine many standard functions and operators. Some examples of allowed elements are collected in the table.
| Type | Examples |
|---|---|
| Arithmetic | +, -, *, /, ^ |
| Constants | pi, e |
| Elementary | sin(x), cos(x), tan(x) |
| Exponential | exp(x), log(x), log10(x) |
| Power, root | sqrt(x), pow(x,y) |
| Conditional | x>0 ? [0]*x : 0 |
With expression strings you do not need to write any extra C++ code or compile anything. For many simple models, for example sums of exponentials or piecewise parametrizations, this approach is sufficient.
Keep in mind that expression based functions are interpreted. For very complex models or for heavy repeated use, compiled C++ functions are usually faster and more flexible.
C++ function based custom models
For more complicated models you often want a real C++ function. This is especially helpful when the expression would become long and hard to read, when you want to branch on conditions with if, or when you want to reuse other C++ code.
A custom function used with TF1 must have a specific signature. The most common one is:
Double_t myModel(Double_t *x, Double_t *par)Here:
x[0]is the coordinate at which the function is evaluated.par[i]are the fit parameters.
A common pattern is:
Double_t myModel(Double_t *x, Double_t *par) {
Double_t xx = x[0];
// Example: Gaussian peak plus linear background
Double_t gaus = par[0] * TMath::Gaus(xx, par[1], par[2], true);
Double_t bkg = par[3] + par[4]*xx;
return gaus + bkg;
}
In a ROOT macro you would then connect this function to a TF1:
TF1 *f = new TF1("f", myModel, 0.0, 10.0, 5); // 5 parameters: par[0]..par[4]
The last argument of the TF1 constructor is the number of parameters. Inside the function body you are free to compute anything as long as you return a single Double_t value.
If your function depends on more than one variable, for example x and y for use with a 2D histogram or a graph with errors, you can still use the same interface. In this case ROOT will fill x[0], x[1], etc. For the current chapter we focus on one dimensional fits, so only x[0] is used.
Because C++ functions are compiled, they are generally faster than expression strings, especially if you enable macro compilation with ACLiC or use a compiled library. For large data sets or complex physics models, that performance difference becomes important.
Connecting histograms or graphs to custom functions
Once you have a TF1 built from either an expression or a C++ function, you use it with histograms or graphs in the same way as built in shapes.
For a histogram h:
h->Fit(f, "R"); // "R" enforces the function range
For a graph gr:
gr->Fit("f", "R");
ROOT will vary the parameters [0], [1], etc. or par[0], par[1] and so on to find the best fit.
It is important that your function is numerically stable over the chosen range. If the function returns very large values or NaN for some x inside the fit region, the minimization can fail or give meaningless results.
Every custom model used for fitting must return a finite Double_t for all x values inside the fit range, otherwise the fit can fail or converge to wrong parameter values.
Setting initial parameters
Before running a fit ROOT needs starting values for each parameter. These initial parameters have a strong practical impact on the quality and speed of the fit, especially for non linear and multi parameter models. Good starting values help the minimizer find the correct minimum and avoid unphysical solutions.
Setting parameter values and names
When you create a TF1, the parameters are initialized to default values, usually zero. You almost always want to overwrite these defaults with values that are close to what you expect from the data.
For a TF1 *f you can do:
f->SetParameter(0, 100.0); // par[0]
f->SetParameter(1, 1.0); // par[1]
f->SetParameter(2, 0.5); // par[2]You can also give each parameter a name. This affects both the display on plots and how results are printed:
f->SetParName(0, "Norm");
f->SetParName(1, "Slope");
f->SetParName(2, "Offset");
For a model defined as a C++ function, the indices of par[] in your function body and the indices you set with SetParameter must match.
To read back a parameter value, for example after a fit, you use:
double p0 = f->GetParameter(0);
double p0err = f->GetParError(0);Reading parameter values and uncertainties in more detail is treated in the dedicated chapter on fit parameters.
Choosing reasonable starting values
ROOT can perform a fit even with very rough starting values, but there are a few simple strategies to choose good initial parameters.
For a Gaussian peak model, a typical choice is:
- Estimate the mean by inspection of the histogram. Use the bin center of the visible peak or approximate it by eye.
- Estimate the width from the spread of the peak, for example take a value comparable to the full width at half maximum divided by $2.35$.
- Estimate the amplitude from the peak height.
Translated into code for a histogram h:
int binMax = h->GetMaximumBin();
double xMax = h->GetXaxis()->GetBinCenter(binMax);
double yMax = h->GetMaximum();
// Example: Gaussian + constant background
TF1 *g = new TF1("g", "[0]*exp(-0.5*((x-[1])/[2])^2) + [3]", xMax-5, xMax+5);
g->SetParameter(0, yMax); // amplitude
g->SetParameter(1, xMax); // mean
g->SetParameter(2, 1.0); // width guess
g->SetParameter(3, 0.0); // background
For linear or polynomial models, you can often estimate slopes and intercepts from a quick look at the data. For example, if your y values increase by roughly 10 when x increases by 1, an initial slope of 10 is reasonable.
If you are unsure, you can start with simple generic values such as 1.0 for scale parameters and 0.0 for offsets. However, for models with many parameters or strong correlations between parameters, poor initial values can lead to a failed fit or convergence to a local minimum that does not represent the data.
Always set initial parameter values that are physically sensible and roughly compatible with the data. Poor starting values are one of the most common reasons for unstable or failed fits with custom functions.
Setting parameter limits
Sometimes you know that a parameter must stay within a certain range. For example, a width must be positive, or a physical efficiency must lie between 0 and 1. You can tell ROOT about such constraints using parameter limits:
f->SetParLimits(0, 0.0, 1e6); // par[0] restricted to [0, 1e6]
f->SetParLimits(1, -10.0, 10.0);If you set both limits to the same value, the parameter is fixed and not varied in the fit:
f->SetParLimits(2, 0.5, 0.5); // par[2] fixed at 0.5It is often better to use sensible limits than to leave parameters completely free, especially in custom models where unphysical parameter combinations can lead to numerical problems.
Parameter limits should not be too tight, otherwise the best fit might lie outside the allowed region, and the reported uncertainties will be misleading. Use them to exclude values you know are impossible, not to force the fit to your expectation.
Using different starting values for repeated fits
In some workflows you may want to fit many histograms or many subsets of data with the same model. In those cases you should update the initial parameters for each fit according to the current data, for example by reusing the previous fit result as the new start, or by estimating peak positions and scales for each histogram.
Within a loop over histograms h[i], you can do:
// For each histogram, set a new starting mean from the maximum bin
int binMax = h[i]->GetMaximumBin();
double xMax = h[i]->GetXaxis()->GetBinCenter(binMax);
f->SetParameter(1, xMax); // update mean guess for this histogram
h[i]->Fit(f, "RQ"); // quiet fit with rangeBy adapting initial parameters to each data subset, you make it easier for ROOT to converge quickly and reliably, even when peaks move or amplitudes change from one histogram to the next.
With well chosen initial parameters and, when necessary, reasonable parameter limits, your custom fit functions become powerful and robust tools for modeling real experimental data.
Views: 13
KAHIBARO