10.5. Function Integrals

Table of Contents

Numerical integration

In ROOT, function integrals are usually computed through the TF1 class. Whenever you integrate a function numerically, you are approximating the area under a curve by evaluating the function at many points and combining these values according to some numerical algorithm. ROOT hides most of this complexity behind simple methods, but it is important to understand the basic usage pattern.

For a TF1 named f, the most common way to perform a numerical integration is to use the Integral method. You provide a lower and an upper bound, and ROOT returns a double that approximates the integral
$$
\int_{x_{\text{min}}}^{x_{\text{max}}} f(x)\,\mathrm{d}x.
$$

For example, if you have defined a Gaussian function:

cpp

TF1 *gaus = new TF1("gaus","gaus", -5, 5);  // built-in Gaussian
gaus->SetParameters(1.0, 0.0, 1.0);        // amplitude, mean, sigma

you can compute its integral between -1 and 1 as follows:

cpp

double area = gaus->Integral(-1.0, 1.0);

The result is a numerical approximation. ROOT chooses the integration step sizes internally, and adapts them to the shape of the function. This is usually sufficient for smooth functions, which are the typical case for physics models.

There is an important distinction between integrating a pure mathematical expression and integrating a function that internally calls more complex C++ code. TF1 can represent both. When you define a function using a formula string, such as "sin(x)/x", ROOT uses a fast internal evaluator. When you define a function from a user C++ function, ROOT calls your function at many x-values during the integration. Either way, the interface through Integral is the same.

If your function depends on parameters, such as a Gaussian with amplitude, mean, and sigma, the integral automatically uses the current parameter values stored in the TF1. You can modify parameters with SetParameter or SetParameters and then call Integral again to see how the integral changes.

The integration limits that you pass to Integral do not have to match the range that you used when creating the TF1. You can integrate over any subrange or a larger range, for example

cpp

double full = gaus->Integral(-5.0, 5.0);   // more than the definition range

ROOT will still evaluate the function at the requested x values.

For some functions, especially those with sharp peaks or non-smooth behavior, default numerical settings might not be precise enough. TF1 offers variants of Integral that accept additional arguments, where you can control the maximum allowed relative and absolute error, and the maximum number of function evaluations. These are more advanced and typically become relevant only in precision studies or complicated models.

Whenever you rely on a numerical integral for a physics result, you should develop the habit of checking its stability. You can compare integrals computed with slightly different accuracy settings, or compare with an approximate analytical result when it is available. ROOT does not warn you automatically if the integration was difficult, so this manual cross check is good practice.

For a TF1 named f, the basic numerical integral over a range is obtained with
double I = f->Integral(xmin, xmax);
which approximates
$\displaystyle I \approx \int_{x_{\min}}^{x_{\max}} f(x)\,\mathrm{d}x.$

Calculating areas

In many analyses, an integral is interpreted as an area: for example, the area under a peak in a distribution, or the total probability in a given interval. TF1 is often used to model such peaks, and integration is then used to extract physical quantities like yields or probabilities.

When you use a normalized function, the total area over the full range is 1. For example, the standard Gaussian probability density function
$$
g(x) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)
$$
is normalized so that
$$
\int_{-\infty}^{\infty} g(x)\,\mathrm{d}x = 1.
$$
If your TF1 represents a properly normalized probability density, then the integral between two values $a$ and $b$ gives the probability that the variable lies in that interval:
$$
P(a < X < b) = \int_a^b g(x)\,\mathrm{d}x.
$$

In ROOT, this is simply:

cpp

double prob = g->Integral(a, b);

When your function has a free normalization parameter, the integral represents an expected count or yield instead of a probability. For instance, suppose you fit a Gaussian to a histogram of events and the Gaussian amplitude parameter is proportional to the number of events in the peak. The integral of the fitted function over an interval around the peak is then the estimated number of events in that region. This is a standard way to extract signal yields from background models.

Often, you will compare the area in a subrange to the total area. This gives a fraction, for example the fraction of events in a signal window relative to the total modeled events:

cpp

double total = f->Integral(xmin_total, xmax_total);
double window = f->Integral(xmin_win, xmax_win);
double fraction = window / total;

This is especially useful for efficiency calculations, where the function represents an ideal distribution and you want to know what fraction lies inside a detector acceptance or selection cut. The key point is that you interpret the integral as an area, and the ratio of areas as an efficiency or probability.

A frequent practical task is to convert between a function area and a histogram count. If you have a histogram with bin width $\Delta x$ and bin content $N_i$ in bin $i$, then the total number of entries in a range is approximately the sum of bin contents in that range. The integral of a TF1 fitted to that histogram should match this sum if the function parameters have been chosen appropriately. Integration helps you translate the smooth fitted curve into a count that is comparable to your discrete histogram data.

Another common area calculation is normalization. Suppose you have defined a function that is not normalized. You can compute its integral over some range, and then divide the function by that integral so that the new function has unit area in that range. In code, you can retrieve the integral and then scale the function using SetParameter or by creating a new TF1 that multiplies the original expression by a constant factor. This is the basic idea behind building probability density functions from arbitrary shapes.

Finally, note that TF1 provides methods closely related to integrals, such as IntegralError, which propagate parameter uncertainties into an uncertainty on the integral. This allows you to attach an error estimate to your area, which is essential when the area is used as a physics result, for example a measured cross section or branching ratio.

Interpreting TF1 integrals as areas:

If $f(x)$ is a probability density normalized to 1, then
$\displaystyle \int_a^b f(x)\,\mathrm{d}x$ gives the probability for $x$ in $[a,b]$.
If $f(x)$ is scaled to match event yields, then
$\displaystyle \int_a^b f(x)\,\mathrm{d}x$ gives the expected number of events in $[a,b]$.
Fractions and efficiencies are computed as ratios of integrals,
for example $f_{\text{frac}} = \dfrac{\int_a^b f(x)\,\mathrm{d}x}{\int_{x_{\min}}^{x_{\max}} f(x)\,\mathrm{d}x}$.

10.5. Function Integrals

Numerical integration

Calculating areas

Comments

Where to Move