Manipulaciones y ajustes de histogramas

Librerías de histogramas

Python 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())

[[172, 89, 29, 69, 277, 1.64e+03, 24, 0, 2, 2],
 [0.389, 17.6, 34.7, 51.9, 69.1, 86.2, 103, 121, 138, 155, 172]]
----------------------------------------------------------------
type: 2 * var * float64

Debido a la prominencia de NumPy, este 2-tupla de arreglos (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 alta energía esperan (sub/sobre flujo, 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())

(array([ 172.,   89.,   29.,   69.,  277., 1640.,   24.,    0.,    2.,
           2.]),
 array([  0.38905792,  17.56032889,  34.73159987,  51.90287084,
         69.07414181,  86.24541279, 103.41668376, 120.58795473,
        137.75922571, 154.93049668, 172.10176766]),
 <BarContainer object of 10 artists>)

Además de los mismos contenidos de bins y bordes que NumPy, Matplotlib incluye un gráfico que se puede graficar.

Boost-histogram y hist

La característica principal que falta en estas funciones (sin algún esfuerzo) es la posibilidad de rellenar nuevamente. Los físicos de altas energías generalmente quieren llenar histogramas con más datos de los que pueden caber en memoria, lo que significa establecer intervalos de bins en un contenedor vacío y llenarlo en lotes (secuencialmente o en paralelo).

Boost-histogram es una biblioteca diseñada para ese propósito. Está destinada a ser 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 gráficos, por ejemplo) es proporcionada por una biblioteca llamada “hist.”

import hist

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

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

h.plot()

[StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff605a82d50>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]

Indexación Universal de Histogramas (UHI)

Hay un intento dentro de Scikit-HEP para estandarizar lo que significan las rebanadas similares a arreglos para un histograma. (Ver documentación.)

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

h[10:110].plot()

[StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff60267de80>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]

pero a menudo deseamos seleccionar bins por valor de coordenada.

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

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

[StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff60258be60>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]

o rebinar por un factor,

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

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

[StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff60288c500>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]

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://pivarski-princeton.s3.amazonaws.com/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_vertices[:, :, sum].plot2d_full()

(ColormeshArtists(pcolormesh=<matplotlib.collections.QuadMesh object at 0x7ff5fa54e990>, cbar=None, text=[]),
 [StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff5fa54efc0>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)],
 [StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff6026a0e30>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)])

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

(ColormeshArtists(pcolormesh=<matplotlib.collections.QuadMesh object at 0x7ff5fa57e6f0>, cbar=None, text=[]),
 [StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff5fa3f3350>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)],
 [StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff602238f50>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)])

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

[StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff5fa4b8da0>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]

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

[StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff5fa368260>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]

Un objeto de histograma puede tener más dimensiones de las que se pueden visualizar razonablemente; puedes rebanar, rebin y proyectar en algo visual más tarde.

Ajuste de histogramas

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

import numpy as np
import iminuit.cost

norm = len(h.axes[0].widths) / (h.axes[0].edges[-1] - h.axes[0].edges[0]) / h.sum()


def f(x, background, mu, gamma):
    return (
        background
        + (1 - background) * gamma**2 / ((x - mu) ** 2 + gamma**2) / np.pi / gamma
    )


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

[<matplotlib.lines.Line2D at 0x7ff5f9f7eea0>]

O a través de zfit, un ajustador similar a RooFit en Python:

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)
model.to_hist().plot(density=1)

/home/runner/miniconda3/envs/skhep-tutorial/lib/python3.12/site-packages/zfit/__init__.py:60: 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(

[StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7ff5d02f0fe0>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]