{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Histogram processing in Coffea\n",
"\n",
"While `coffea` provides its own [`coffea.hist`][coffeahist] package, to better\n",
"interface with the python ecosystem in general, I recommend users to start\n",
"their analysis using the [`scikit hist`][scikithist] package, as it is more\n",
"general purpose, and coffea histogram will be depreciated in newer version of\n",
"`coffea`. All the core functionalities would be the same, but `sckikit hist`\n",
"(which we will refer to as just `hist` for the remainder of this tutorial) has\n",
"more bells and whistles out of the box. A very good documentation on using\n",
"`hist` is already provided by the developers [here][scikithist], so this\n",
"documentation would mainly focus on going over the fundamental concepts of\n",
"using histogram for physics analysis. \n",
"\n",
"One thing that will not be included in this documentation would be the plotting\n",
"of histogram, as there are much too many customization options to go over. Here\n",
"we will be using the in-built `Hist.plot1d` and `Hist.plot2d` with minimal\n",
"styling settings, just to illustrate how the data in the histogram are being\n",
"stored. \n",
"\n",
"Let us start by importing the basic libraries that we will be using: `hist` and\n",
"`numpy` (`numpy` would mainly be used for pseudo-data generation and nothing\n",
"else.)\n",
"\n",
"\n",
"[coffeahist]: https://coffeateam.github.io/coffea/api/coffea.hist.Hist.html\n",
"[scikithist]: https://hist.readthedocs.io/en/latest/?badge=latest"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"import hist \n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Histogram construction in computation\n",
"\n",
"Histograms in the context of an analysis is typically used to tally the number\n",
"of entities that fall under some \"discrete\" criteria. These criteria can either\n",
"be some true discrete variable (like data sets/particle flavor), or some\n",
"continuous variable that we explicitly discrete-ize (like for particle kinematic\n",
"variables). Let use make a very simple histogram that represents the 3 particle\n",
"kinematic variables used in HEP: $p_{T}$, $\\eta$ and $\\phi$:\n",
"\n",
"The filling of the histogram can be done in parallel, supposing that you have\n",
"arrays of identical dimensions. \n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"
\n",
""
],
"text/plain": [
"Hist(\n",
" Regular(100, 20, 100, name='pt', label='pt'),\n",
" Regular(100, -3, 3, name='eta', label='eta'),\n",
" Regular(100, -3.14159, 3.14159, name='phi', label='phi'),\n",
" storage=Double()) # Sum: 5334785.0 (10000000.0 with flow)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h = hist.Hist(hist.axis.Regular(100, 20, 100, name='pt'),\n",
" hist.axis.Regular(100, -3, 3, name='eta'),\n",
" hist.axis.Regular(100, -np.pi, np.pi, name='phi'))\n",
"\n",
"N = 10_000_000\n",
"pt = np.random.exponential(50, size=N)\n",
"eta = np.random.normal(size=N)\n",
"phi = np.random.uniform(-np.pi, np.pi, size=N)\n",
"h.fill(pt=pt, eta=eta, phi=phi)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For a computer, a histogram is just an array of values, indicating the number of\n",
"entries that was filling in \"bin\". Formally, the `hist` package uses what is\n",
"known as the [\"unified histogram interface notation\" (UHI)][uhi], so the\n",
"following\n",
"\n",
"[uhi]: https://uhi.readthedocs.io/en/latest/\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"20.0"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h[0, 50, 15]\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Simply indicates the histogram value stored in the 0-th bin of the $p_T$-axis,\n",
"the 50-th bin of the $\\eta$-axis, and the 15-th bin of the $\\phi$ axis.\n",
"\n",
"As working with raw bin indices is rather cumbersome, the UHI notation allows the\n",
"use of `j` as a postfix to indicate that one want the bin that contains some\n",
"numerical value, for instance:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"26.0"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h[25j, 0.5j, -1.2j]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Would indicate the histogram value of the bin that contains $p_{T}$ value 25,\n",
"$\\eta$ value 0.5 and $\\phi$ of value -1.2.\n",
"\n",
"Here we are treating the histogram as an array, knowing that the first index\n",
"corresponds to the $p_{T}$ UHI bin index, the second corresponding to the $\\eta$\n",
"UHI bin index and so on. As long as one has declared the histogram axis with a\n",
"name, a dictionary can also be used to explicitly specify the bin index for some\n",
"named axis, rather than having to keep track of the axis ordering manually. For\n",
"example the two expressions below"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"26.0\n",
"26.0\n"
]
}
],
"source": [
"print(h[25j, 0.5j, -1.2j])\n",
"print(h[{'phi': -1.2j, 'pt': 25j, 'eta': 0.5j}])\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Means entirely the same UHI index, just one is much more human friendly with\n",
"more explicitly notations of which bins are being requested without needed to\n",
"tracking the axis ordering. \n",
"\n",
"\n",
"## Histogram slicing and truncation\n",
"\n",
"Very rarely will we need to operate directly on raw histogram values. Most of\n",
"the time, the rationale for using a high dimensional histogram like this would\n",
"be that we can consider the most general distribution among all variables of\n",
"interest, and retroactively zoom in on some region of phase space of interest.\n",
"That is where the real power of thinking of histograms as arrays comes in: the\n",
"same notation used to perform array slicing in `numpy` can be used here for\n",
"histograms as well. Supposed we are only interested in the case where $\\phi$ of\n",
"the particle lies in the same bin as some value (may that specific region of the\n",
"detector is faulty), what we can so would be something like: \n"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"