Histogram manipulations and fitting

Histogram libraries¶

Mainstream Python has libraries for filling histograms.

NumPy¶

NumPy, for instance, has an np.histogram function.

import skhep_testdata, uproot

tree = uproot.open(skhep_testdata.data_path("uproot-Zmumu.root"))["events"]

import numpy as np

np.histogram(tree["M"].array())

(<Array [172, 89, 29, 69, 277, 1640, 24, 0, 2, 2] type='10 * int64'>,
 <Array [0.389, 17.6, 34.7, 51.9, ..., 121, 138, 155, 172] type='11 * float64'>)

Because of NumPy’s prominence, this 2-tuple of arrays (bin contents and edges) is a widely recognized histogram format, though it lacks many of the features high-energy physicists expect (under/overflow, axis labels, uncertainties, etc.).

Matplotlib¶

Matplotlib also has a plt.hist function.

import matplotlib.pyplot as plt

plt.hist(tree["M"].array());

In addition to the same bin contents and edges as NumPy, Matplotlib includes a plottable graphic.

Boost-histogram and hist¶

The main feature that these functions lack (without some effort) is refillability. High-energy physicists usually want to fill histograms with more data than can fit in memory, which means setting bin intervals on an empty container and filling it in batches (sequentially or in parallel).

Boost-histogram is a library designed for that purpose. It is intended as an infrastructure component. You can explore its “low-level” functionality upon importing it:

import boost_histogram as bh

A more user-friendly layer (with plotting, for instance) is provided by a library called “hist.”

import hist

h = hist.Hist(hist.axis.Regular(120, 60, 120, name="mass"))

h.fill(tree["M"].array())

h.plot();

Universal Histogram Indexing (UHI)¶

There is an attempt within Scikit-HEP to standardize what array-like slices mean for a histogram. (See documentation.)

Naturally, integer slices should select a range of bins,

h[10:110].plot();

but often you want to select bins by coordinate value

# Explicit version
# h[hist.loc(90) :].plot();

# Short version
h[90j:].plot();

or rebin by a factor,

# Explicit version
# h[:: hist.rebin(2)].plot();

# Short version
h[::2j].plot();

or sum over a range.

# Explicit version
# h[hist.loc(80) : hist.loc(100) : sum]

# Short version
h[90j:100j:sum]

1102.0

Things get more interesting when a histogram has multiple dimensions.

import uproot
import hist
import awkward as ak

picodst = uproot.open(
    "https://zenodo.org/records/21777191/files/pythia_ppZee_run17emb.picoDst.root:PicoDst"
)

vertexhist = hist.Hist(
    hist.axis.Regular(600, -1, 1, label="x"),
    hist.axis.Regular(600, -1, 1, label="y"),
    hist.axis.Regular(40, -200, 200, label="z"),
)

vertex_data = picodst.arrays(filter_name="*mPrimaryVertex[XYZ]")

vertexhist.fill(
    ak.flatten(vertex_data["Event.mPrimaryVertexX"]),
    ak.flatten(vertex_data["Event.mPrimaryVertexY"]),
    ak.flatten(vertex_data["Event.mPrimaryVertexZ"]),
)

Hist(
  Regular(600, -1, 1, label='x'),
  Regular(600, -1, 1, label='y'),
  Regular(40, -200, 200, label='z'),
  storage=Double()) # Sum: 8004.0

This histogram has three axes, so each of the following plots picks out a different view of it. Putting sum in an axis’s slot integrates that axis away: summing over z leaves the distribution of collision points transverse to the beam. plot2d_full draws that 2D distribution together with its x and y projections.

vertexhist[:, :, sum].plot2d_full();

The same view, zoomed in on the beam spot. As with the 1D examples above, the j suffix selects by coordinate value rather than by bin number, so this keeps the region from -0.25 to 0.25 on both axes.

vertexhist[-0.25j:0.25j, -0.25j:0.25j, sum].plot2d_full();

Summing over x and y instead leaves the distribution along the beam line, which is much broader than the transverse spread.

vertexhist[sum, sum, :].plot();

Selecting and summing can be combined in one slice: -0.25j:0.25j:sum keeps only that coordinate range and then sums over it. This is the z distribution of the collisions inside the beam spot. It comes out nearly identical to the previous plot, because almost all of the collisions were already in that central region.

vertexhist[-0.25j:0.25j:sum, -0.25j:0.25j:sum, :].plot();

A histogram object can have more dimensions than you can reasonably visualize—you can slice, rebin, and project it into something visual later.

Fitting histograms¶

By directly writing a loss function in Minuit:

import numpy as np
import iminuit.cost

xmin, xmax = h.axes[0].edges[0], h.axes[0].edges[-1]

# rescale the bin counts into a probability density, so that they can be
# compared against a normalized model
norm = len(h.axes[0].widths) / (xmax - xmin) / h.sum()


def f(x, background, mu, gamma):
    # Cauchy (Breit-Wigner) peak, normalized over the fitted range
    peak = gamma / ((x - mu) ** 2 + gamma**2) / np.pi
    peak /= (np.arctan((xmax - mu) / gamma) - np.arctan((xmin - mu) / gamma)) / np.pi
    # flat background, also normalized over the fitted range, so that
    # `background` is the fraction of events in the background
    return background / (xmax - xmin) + (1 - background) * peak


loss = iminuit.cost.LeastSquares(
    h.axes[0].centers, h.values() * norm, np.sqrt(h.variances()) * norm, f
)
loss.mask = h.variances() > 0

minimizer = iminuit.Minuit(loss, background=0, mu=91, gamma=4)

minimizer.migrad()
minimizer.hesse()

(h * norm).plot()
plt.plot(loss.x, f(loss.x, *minimizer.values));

Or through zfit, a Pythonic RooFit-like fitter. This builds the same model — a Cauchy peak plus a flat background, with background as the background fraction — so the fitted parameters come out similar to the ones above. They are not identical, because this fit minimizes a binned negative log-likelihood instead of a least-squares cost.

import zfit

binned_data = zfit.data.BinnedData.from_hist(h)

binning = zfit.binned.RegularBinning(120, 60, 120, name="mass")
space = zfit.Space("mass", binning=binning)

background = zfit.Parameter("background", 0)
mu = zfit.Parameter("mu", 91)
gamma = zfit.Parameter("gamma", 4)
unbinned_model = zfit.pdf.SumPDF(
    [zfit.pdf.Uniform(60, 120, space), zfit.pdf.Cauchy(mu, gamma, space)], [background]
)

model = zfit.pdf.BinnedFromUnbinnedPDF(unbinned_model, space)
loss = zfit.loss.BinnedNLL(model, binned_data)

minimizer = zfit.minimize.Minuit()
result = minimizer.minimize(loss)

binned_data.to_hist().plot(density=1)

# The model is a normalized pdf, so its bins carry no statistical uncertainty.
# Draw it as a curve: asking hist to plot it would try to put Poisson error
# bars on zero variances, which divides by zero.
model_hist = model.to_hist()
model_axis = model_hist.axes[0]
plt.plot(
    model_axis.centers,
    model_hist.values() / (model_hist.values().sum() * model_axis.widths),
);

/home/runner/miniconda3/envs/skhep-tutorial/lib/python3.12/site-packages/zfit/__init__.py:93: UserWarning: TensorFlow warnings are by default suppressed by zfit. In order to show them, set the environment variable ZFIT_DISABLE_TF_WARNINGS=0. In order to suppress the TensorFlow warnings AND this warning, set ZFIT_DISABLE_TF_WARNINGS=1.
  warnings.warn(