Manipulación y ajuste de histogramas - Tutorial de Scikit-HEP

Librerías de histogramas¶

El Python convencional ya tiene librerías para llenar histogramas.

NumPy¶

NumPy, por ejemplo, tiene una función np.histogram.

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'>)

Debido a la prominencia de NumPy, esta tupla de dos arrays (contenidos de los bins y bordes) es un formato de histograma ampliamente reconocido, aunque carece de muchas de las características que los físicos de altas energías esperan (underflow/overflow, etiquetas de ejes, incertidumbres, etc.).

Matplotlib¶

Matplotlib también tiene una función plt.hist.

import matplotlib.pyplot as plt

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

Además de los mismos contenidos de bins y bordes que NumPy, Matplotlib incluye un gráfico listo para mostrar.

Boost-histogram y hist¶

La característica principal que les falta a estas funciones (sin algo de esfuerzo) es la posibilidad de volver a llenarlos. Los físicos de altas energías generalmente quieren llenar histogramas con más datos de los que caben en memoria, lo que significa establecer los intervalos de los bins en un contenedor vacío y llenarlo por lotes (secuencialmente o en paralelo).

Boost-histogram es una librería diseñada para ese propósito. Está pensada como un componente de infraestructura. Puedes explorar su funcionalidad “de bajo nivel” al importarla:

import boost_histogram as bh

Una capa más amigable para el usuario (con graficación, por ejemplo) la proporciona una librería llamada “hist”.

import hist

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

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

h.plot();

Indexación Universal de Histogramas (UHI)¶

Dentro de Scikit-HEP hay un esfuerzo por estandarizar lo que significan las rebanadas tipo array para un histograma. (Ver documentación.)

Naturalmente, las rebanadas con enteros deberían seleccionar un rango de bins,

h[10:110].plot();

pero a menudo quieres seleccionar bins por valor de coordenada

# Versión explícita
# h[hist.loc(90) :].plot();

# Versión corta
h[90j:].plot();

o reagrupar los bins por un factor,

# Versión explícita
# h[:: hist.rebin(2)].plot();

# Versión corta
h[::2j].plot();

o sumar sobre un rango.

# Versión explícita
# h[hist.loc(80) : hist.loc(100) : sum]

# Versión corta
h[90j:100j:sum]

1102.0

Las cosas se vuelven más interesantes cuando un histograma tiene múltiples dimensiones.

import uproot
import hist
import awkward as ak

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

hist_vertices = 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"),
)

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

hist_vertices.fill(
    ak.flatten(datos_vertices["Event.mPrimaryVertexX"]),
    ak.flatten(datos_vertices["Event.mPrimaryVertexY"]),
    ak.flatten(datos_vertices["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

Este histograma tiene tres ejes, así que cada uno de los siguientes gráficos escoge una vista distinta de él. Poner sum en la posición de un eje lo integra y lo elimina del resultado: sumar sobre z deja la distribución de los puntos de colisión en el plano transversal al haz. plot2d_full dibuja esa distribución 2D junto con sus proyecciones en x y en y.

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

La misma vista, con un acercamiento a la región de interacción del haz. Igual que en los ejemplos 1D de arriba, el sufijo j selecciona por valor de coordenada en lugar de por número de bin, así que esto conserva la región de -0.25 a 0.25 en ambos ejes.

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

Sumar en cambio sobre x e y deja la distribución a lo largo de la línea del haz, que es mucho más ancha que la dispersión transversal.

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

Seleccionar y sumar se pueden combinar en una sola rebanada: -0.25j:0.25j:sum conserva solo ese rango de coordenadas y luego suma sobre él. Esta es la distribución en z de las colisiones dentro de la región de interacción del haz. Resulta casi idéntica a la del gráfico anterior, porque casi todas las colisiones ya estaban en esa región central.

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

Un objeto histograma puede tener más dimensiones de las que puedes visualizar razonablemente; puedes rebanarlo, reagrupar sus bins y proyectarlo más adelante en algo visual.

Ajuste de histogramas¶

Escribiendo directamente una función de pérdida en Minuit:

import numpy as np
import iminuit.cost

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

# reescala los contenidos de los bins a una densidad de probabilidad, para que
# se puedan comparar con un modelo normalizado
norm = len(h.axes[0].widths) / (xmax - xmin) / h.sum()


def f(x, background, mu, gamma):
    # pico de Cauchy (Breit-Wigner), normalizado en el rango del ajuste
    pico = gamma / ((x - mu) ** 2 + gamma**2) / np.pi
    pico /= (np.arctan((xmax - mu) / gamma) - np.arctan((xmin - mu) / gamma)) / np.pi
    # fondo plano, también normalizado en el rango del ajuste, de modo que
    # `background` es la fracción de eventos que hay en el fondo
    return background / (xmax - xmin) + (1 - background) * pico


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));

O a través de zfit, un ajustador pythónico al estilo de RooFit. Esto construye el mismo modelo (un pico de Cauchy más un fondo plano, con background como la fracción de fondo), así que los parámetros ajustados salen parecidos a los de arriba. No son idénticos, porque este ajuste minimiza una log-verosimilitud negativa calculada sobre los bins, en lugar de un costo de mínimos cuadrados.

import zfit

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

binning = zfit.binned.RegularBinning(120, 60, 120, name="masa")
space = zfit.Space("masa", 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)

# El modelo es una pdf normalizada, así que sus bins no tienen incertidumbre
# estadística. Lo dibujamos como una curva: pedirle a hist que lo grafique
# intentaría poner barras de error de Poisson sobre varianzas nulas, lo que
# implica una división por cero.
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(