04_Scikit-HEP_histogramming-visualisation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1><b><span style=\"color:blue\">Histogramming & visualisation</span></b></h1></center>\n",
    "\n",
    "&nbsp;<br>\n",
    "This is not an under-statement - **histograms are ubiquitous in Particle Physics!** This being said, and as far as the (general) scientific Python ecosystem is concerned, histograms aren't quite seen as first-class citizens ... Would you then be surprised to know that histogramming saw a lot of developments in our field?\n",
    "\n",
    "### **Quick intro to the following packages**\n",
    "- `histoprint` - pretty print of NumPy (and other) histograms to the console.\n",
    "- `boost-histogram` - fast and flexible multi-dimensional histograms and profiles.\n",
    "- `Hist` - analyst-friendly front-end for boost-histogram.\n",
    "- `mplhep` - HEP domain-specific add-ons to matplotlib.\n",
    "\n",
    "**Important note:** much of the functionality in these packages significantly enhances the general scientific ecosystem and is at all not (HEP) domain specific!\n",
    "\n",
    "**Important note (2):** the fact that there is some overlap in some of the functionality these packages provide may be overwhelming at first - no different from the fact that very many Python libraries exist for visualisation out there, actually.\n",
    "Remember, the development of these packages are largely driven by the community (needs and feedback). See this as a motivation to contribute towards the most useful and best features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&nbsp;<br>\n",
    "\n",
    "<center><h2><b><span style=\"color:green\">histoprint - pretty print of NumPy (and other) histograms to the console</span></b></h2></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An introduction to `histoprint` very largely based on the README, see https://github.com/scikit-hep/histoprint.\n",
    "\n",
    "The package is mainly meant to be used on the console; else for quick and trivial displays of 1-dimensional histograms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from histoprint import text_hist, print_hist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`text_hist(...)`, the simplest ever thin wrapper for `numpy.histogram`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "text_hist(\n",
    "    np.random.randn(1000), bins=[-5, -3, -2, -1, -0.5, 0, 0.5, 1, 2, 3, 5], title=\"Variable bin widths\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_numpy = np.histogram(np.random.normal(0, 1, 100000), bins=20, range=(-5, 5))\n",
    "h_numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.randn(1000) - 2\n",
    "B = np.random.randn(1000)\n",
    "C = np.random.randn(1000) + 2\n",
    "D = np.random.randn(500) * 2\n",
    "\n",
    "histA = np.histogram(A, bins=15, range=(-5, 5))\n",
    "histB = np.histogram(B, bins=15, range=(-5, 5))\n",
    "histC = np.histogram(C, bins=15, range=(-5, 5))\n",
    "histD = np.histogram(D, bins=15, range=(-5, 5))\n",
    "histAll = ([histA[0], histB[0], histC[0], histD[0]], histA[1])\n",
    "\n",
    "# print_hist can be used to print multiple histograms at once\n",
    "# (or just to print a single one as returned by numpy.histogram)\n",
    "print_hist(histAll, title=\"Overlays\", labels=\"ABCDE\")\n",
    "print_hist(\n",
    "    histAll,\n",
    "    title=\"Stacks\",\n",
    "    stack=True,\n",
    "    symbols=\"      \",\n",
    "    bg_colors=\"rgbcmy\",\n",
    "    labels=\"ABCDE\",\n",
    ")\n",
    "print_hist(\n",
    "    histAll,\n",
    "    title=\"Summaries\",\n",
    "    symbols=r\"=|\\/\",\n",
    "    fg_colors=\"0\",\n",
    "    bg_colors=\"0\",\n",
    "    labels=[\"AAAAAAAAAAAAAAAA\", \"B\", \"CCCCCCCCCCCCC\", \"D\"],\n",
    "    summary=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** the last example does not use terminal colors, so it can be copied as text."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **1. Command-line interface**\n",
    "\n",
    "`Histoprint` comes with a simple command-line interface to create histograms of tabulated data.\n",
    "It can read in files or take data directly from STDIN:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!histoprint --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **2. Interoperability - support for other histogram types**\n",
    "\n",
    "`Histoprint` can directly plot other (more fancy) types of histograms if they follow the PlottableProtocol conventions,\n",
    "or offer a way of being converted to the NumPy format.\n",
    "Currently this means they have to expose a numpy() or to_numpy() method.\n",
    "Both the TH1 histograms of uproot4, as well as the histograms of boost-histogram, are supported like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import boost_histogram as bh\n",
    "\n",
    "hist = bh.Histogram(bh.axis.Regular(20, -3, 3))\n",
    "hist.fill(np.random.randn(1000))\n",
    "print_hist(hist, title=\"Boost Histogram\")\n",
    "\n",
    "import uproot\n",
    "file = uproot.open(\"http://scikit-hep.org/uproot3/examples/Event.root\")\n",
    "hist = file[\"htime\"]\n",
    "print_hist(hist, title=\"uproot TH1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&nbsp;<br><center><img src=\"images/logo_boost-histogram.png\" alt=\"boost-histogram package logo\" style=\"width: 150px;\"/></center>\n",
    "\n",
    "<center><h2><b><span style=\"color:green\">Fast and flexible multi-dimensional histograms and profiles</span></b></h2></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Python bindings for the C++14 Boost.Histogram library.\n",
    "* No external dependency.\n",
    "* Very fast and feature-rich.\n",
    "* Multi-dimensional histograms and profiles.\n",
    "* Supports weighted and unweighted data.\n",
    "* And much much more, see the [docs](https://boost-histogram.readthedocs.io/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "<b>Acknowledgements</b>\n",
    "\n",
    "This mini-tutorial is largely based on material by Henry Schreiner (Princeton University), with minor modifications.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **1. Basic 1D histogram**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import boost_histogram as bh\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's generate some data and create a histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data1 = np.random.normal(3.5, 2.5, size=1_000_000)\n",
    "\n",
    "h1 = bh.Histogram(bh.axis.Regular(40, -2, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h1.fill(data1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's explicitly check to see how many entries are in the histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h1.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happened to the missing items? They are in the underflow and overflow bins ;-):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h1.sum(flow=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As is ROOT, overflow bins are on by default. We can turn them off, but they enable some powerful things such as projections.\n",
    "\n",
    "Make sure things visually look OK:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.bar(h1.axes[0].centers, h1, width=h1.axes[0].widths);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From now on, let's be lazy with a helper function for plotting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plothist = lambda h: plt.bar(h.axes[0].centers, h, width=h.axes[0].widths);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **2. Drop-in replacement for NumPy**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins2, edges2 = bh.numpy.histogram(data1, bins=10)\n",
    "bins2, edges2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **3. \"Extra dimensions\"**\n",
    "The same API works for multiple dimensions. This is an area where the library shines particularly well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2 = bh.Histogram(\n",
    "    bh.axis.Regular(150, -1.5, 1.5),\n",
    "    bh.axis.Regular(100, -1, 1)\n",
    ")\n",
    "\n",
    "h2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_2D_data(*, mean=(0,0), widths=(1,1), size=1_000_000):\n",
    "    cov = np.asarray(widths) * np.eye(2)\n",
    "    return np.random.multivariate_normal(mean, cov, size=size).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data3x = make_2D_data(mean=[-.75, .5], widths=[.2, 0.02])\n",
    "data3y = make_2D_data(mean=[.75, .5], widths=[.2, 0.02])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From here on we will be using `.reset()` before a `.fill()`, just to make sure each cell in the notebook can be rerun."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2.reset()\n",
    "h2.fill(*data3x)\n",
    "h2.fill(*data3y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, let's make plotting a little function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plothist2d(h):\n",
    "    X, Y = h.axes.edges\n",
    "    X, Y = np.broadcast_arrays(X, Y)\n",
    "\n",
    "    return plt.pcolormesh(X, Y, h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The above code is the same as doing:\n",
    "> \n",
    "> ```python\n",
    "X, Y = np.meshgrid(*(a.edges for a in h.axes, indexing='ij')\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plothist2d(h2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try a 3D histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data3d = [np.random.normal(size=1_000_000) for _ in range(3)]\n",
    "\n",
    "h3d = bh.Histogram(\n",
    "    bh.axis.Regular(150, -5, 5),\n",
    "    bh.axis.Regular(100, -5, 5),\n",
    "    bh.axis.Regular(100, -5, 5)\n",
    ")\n",
    "\n",
    "h3d.fill(*data3d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Projection of the first two axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plothist2d(h3d.project(0,1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **4. Unified Histogram Indexing**\n",
    "\n",
    "Let's explore the boost-histogram UHI syntax. We will reuse the previous 2D histogram from part 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plothist2d(h2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I can see that I want y from 0.25 to 0.75, in data coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plothist2d(h2[:, bh.loc(.25):bh.loc(.75)]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's the contents of a bin?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2[100,87]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How about in data coordinates?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2[bh.loc(.5), bh.loc(.75)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: to get the coordinates manually:\n",
    "> \n",
    "> ```python\n",
    "h2.axes[0].index(.5) == 100\n",
    "h2.axes[1].index(.75) == 87\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at one part and rebin:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plothist2d(h2[:50:bh.rebin(2), 50::bh.rebin(2)]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the value at `(-.75, .5)`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2[bh.loc(-.75), bh.loc(.5)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **5. Fancier axis types - a circular axis example**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = bh.Histogram(bh.axis.Regular(30, 0, 2*np.pi, circular=True))\n",
    "h.fill(np.random.uniform(0, np.pi*4, size=300))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.subplot(111, polar=True)\n",
    "\n",
    "plothist = lambda h: plt.bar(h.axes[0].centers, h, width=h.axes[0].widths)\n",
    "plothist(h);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **6. Advanced usage - storage and accumulators**\n",
    "\n",
    "Explore:\n",
    "- The extensive list of histogram operations available.\n",
    "- The rich set of storage types (Double, Int64, Weight, etc.).\n",
    "- The accumulators defined (Sum, WeightedSum, etc.). Storages store accumulators.\n",
    "\n",
    "A simple example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = bh.accumulators.Mean()\n",
    "mean.fill([0.3, 0.4, 0.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = bh.Histogram(bh.axis.Regular(10, 0, 10), storage=bh.storage.Mean())\n",
    "\n",
    "h.fill([2.5] * 3, sample=[0.3, 0.4, 0.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_hist(h.to_numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&nbsp;<br><center><img src=\"images/logo_mplhep.png\" alt=\"mplhep package logo\" style=\"width:170px;\"/></center>\n",
    "\n",
    "<center><h2><b><span style=\"color:green\">HEP domain-specific add-ons to matplotlib</span></b></h2></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Package motivation \n",
    "\n",
    "- Collection of wrappers for `matplotlib` to more easily produce plots typically needed in HEP.\n",
    "- Collect styles that are compatible with current collaboration recommendations (ROOT-package-like plots for large [LHC](https://home.cern/science/accelerators/large-hadron-collider) experiments ATLAS, CMS, etc.)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import mplhep"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **1. Some basics**\n",
    "\n",
    "Plotting of 1D histograms made easy for \"standard\" and commonly used looks in HEP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use([mplhep.style.ROOT, mplhep.style.firamath])  # see style down in this notebook\n",
    "\n",
    "h, bins = np.histogram(np.random.normal(10, 3, 400), bins=10)\n",
    "\n",
    "f, axs = plt.subplots(2, 2, sharex=True, sharey=True)\n",
    "axs = axs.flatten()\n",
    "\n",
    "axs[0].set_title(\"Default\", fontsize=18)\n",
    "mplhep.histplot(h, bins, ax=axs[0])\n",
    "\n",
    "axs[1].set_title(\"Plot Edges\", fontsize=18)\n",
    "mplhep.histplot(h, bins, edges=True, ax=axs[1])\n",
    "\n",
    "axs[2].set_title(\"Plot Errorbars\", fontsize=18)\n",
    "mplhep.histplot(h, bins, yerr=np.sqrt(h), ax=axs[2])\n",
    "\n",
    "axs[3].set_title(\"Filled Histogram\", fontsize=18)\n",
    "mplhep.histplot(h, bins, histtype='fill', ax=axs[3])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is just as easy to overlay 2, or more, 1D histograms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h, bins = np.histogram(np.random.normal(10, 3, 400), bins = 10)\n",
    "\n",
    "f, axs = plt.subplots(2, 2, sharex=True, sharey=True)\n",
    "axs = axs.flatten()\n",
    "\n",
    "axs[0].set_title(\"Default Overlay\", fontsize=18)\n",
    "mplhep.histplot([h, 1.5*h], bins, ax=axs[0])\n",
    "\n",
    "axs[1].set_title(\"Default Overlay w/ Errorbars\", fontsize=18)\n",
    "mplhep.histplot([h, 1.5*h], bins, yerr=[np.sqrt(h), np.sqrt(1.5*h)], ax=axs[1])\n",
    "\n",
    "\n",
    "axs[2].set_title(\"Automatic Errorbars\", fontsize=18)\n",
    "mplhep.histplot([h, 1.5*h], bins, yerr=True, ax=axs[2])\n",
    "\n",
    "axs[3].set_title(\"With Labels\", fontsize=18)\n",
    "mplhep.histplot([h, 1.5*h], bins, yerr=True, ax=axs[3], label=[\"First\", \"Second\"])\n",
    "\n",
    "axs[3].legend(fontsize=15)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In other situations one rather needs to stack a series of 1D histograms. Here is how this can be done:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axs = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(9, 9))\n",
    "axs = axs.flatten()\n",
    "\n",
    "axs[0].set_title(\"Default\", fontsize=18)\n",
    "mplhep.histplot([h, 1.5*h], bins, stack=True, ax=axs[0])\n",
    "\n",
    "axs[1].set_title(\"Plot Edges\", fontsize=18)\n",
    "mplhep.histplot([h, 1.5*h], bins, edges=True, stack=True,ax=axs[1])\n",
    "\n",
    "axs[2].set_title(\"Plot Errorbars\", fontsize=18)\n",
    "mplhep.histplot([h, 1.5*h], bins, yerr=np.sqrt(h), stack=True, ax=axs[2])\n",
    "\n",
    "axs[3].set_title(\"Filled Histogram\", fontsize=18)\n",
    "mplhep.histplot([1.5*h, h], bins, histtype='fill', stack=True, ax=axs[3])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **2. Fancier plots**\n",
    "\n",
    "Likewise, `mplhep` provides a simple way to plot 2D histograms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "xedges = np.arange(0, 11.5, 1.5)\n",
    "yedges = [0, 2, 3, 4, 6, 7]\n",
    "x = np.random.normal(5, 1.5, 100)\n",
    "y = np.random.normal(4, 1, 100)\n",
    "\n",
    "mplhep.hist2dplot(*np.histogram2d(x, y, bins=(xedges, yedges)), labels=True);\n",
    "\n",
    "# This would work just as well:\n",
    "# H, xedges, yedges = np.histogram2d(x, y, bins=(xedges, yedges))\n",
    "# mplhep.hist2dplot(H, xedges, yedges, labels=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's no trouble if starting from a `boost-histogram` histogram - just convert it on-the-fly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "h2 = bh.Histogram(\n",
    "    bh.axis.Variable(xedges),\n",
    "    bh.axis.Variable(yedges)\n",
    ")\n",
    "h2.fill(x,y)\n",
    "\n",
    "mplhep.hist2dplot(*h2.to_numpy());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, is the conversion at all necessary? Not anymore, as Boost histograms are supported:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mplhep.hist2dplot(h2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **3. Styling**\n",
    "\n",
    "As said above, several styles are predefined. Here is a simple example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.uniform(0, 10, 240)\n",
    "y = np.random.normal(512, 112, 240)\n",
    "z = np.random.normal(0.5, 0.1, 240)\n",
    "\n",
    "plt.style.use([mplhep.style.ROOT, mplhep.style.firamath])  # styles can be chained\n",
    "\n",
    "f, ax = plt.subplots()\n",
    "ax.scatter(x,y, c=z, label='XSX')\n",
    "mplhep.cms.label(loc=0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&nbsp;<br><center><img src=\"images/logo_Hist.png\" style=\"width: 150px;\"/></center>\n",
    "\n",
    "<center><h2><b><span style=\"color:green\">Analyst friendly front-end for boost-histogram</span></b></h2></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `Hist` is a powerful, friendly and analysis-focused histogramming tool based on `boost-histogram`.\n",
    "- It provides some cool shortcuts for histogram creation, plotting tools, and new ideas.\n",
    "\n",
    "Let's explore the library a bit with some examples. See also the [GitHub repository](https://github.com/scikit-hep/hist) for full documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hist\n",
    "from hist import Hist\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **1. Cool representations in notebooks**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hist.new.Reg(50, 1, 2).Double().fill(np.random.normal(1.5, 0.3, 10_000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2 = Hist.new.Reg(50, 0, 2, name='My preferred x-axis title').Reg(50, 10, 20).Double().fill(\n",
    "    np.random.normal(1, 0.5, 10_000), np.random.normal(15, 3, 10_000)\n",
    ")\n",
    "h2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hist.new.Reg(50, 0, 2).Reg(50, 10, 20).Reg(2, 3, 4).Double()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on `boost-histogram`’s Axis, `Hist` supports six types of axis (Regular, Boolean, Variable, Integer, IntCategory and StrCategory), with additional names and labels.\n",
    "\n",
    "You can checkout `hist.axis` or use the shortcut definitions instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the axes using the shortcut method\n",
    "h = (\n",
    "    Hist.new.Reg(10, -5, 5, overflow=False, underflow=False, name=\"A\")\n",
    "    .Bool(name=\"B\")\n",
    "    .Var(range(10), name=\"C\")\n",
    "    .Int(-5, 5, overflow=False, underflow=False, name=\"D\")\n",
    "    .IntCat(range(10), name=\"E\")\n",
    "    .StrCat([\"T\", \"F\"], name=\"F\")\n",
    "    .Double()\n",
    ")\n",
    "h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **2. Handy plotting functions**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data1 = np.random.normal(3.5, 2.5, size=1_000_000)\n",
    "\n",
    "h1 = Hist(hist.axis.Regular(40, -2, 10, name=\"Variable x\"), storage=hist.storage.Double())\n",
    "\n",
    "h1.fill(data1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram \"knows how to plot itself\" - very handy for quick checks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h1.plot(color=\"darkviolet\", lw=3, histtype='fill')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "# plot2d\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "\n",
    "h2.plot2d(ax=ax, cmap=\"plasma\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This being said, `Hist` talks well with the HEP plotting library `mplhep`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = bh.Histogram(\n",
    "bh.axis.Regular(2, 0, 1, metadata=\"x\"),\n",
    "bh.axis.Regular(4, 0, 1, metadata=\"y\"))\n",
    "h.fill(\n",
    "[.2, .4, .3],\n",
    "[.3, .5, .2])\n",
    "h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mplhep.hist2dplot(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did we not say \"handy plotting functions\" and a \"powerful, friendly and analysis-focused histogramming tool\"? Check this out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = Hist(\n",
    "    hist.axis.Regular(50, -5, 5, name=\"S\", label=\"s [units]\", flow=False),\n",
    "    hist.axis.Regular(50, -5, 5, name=\"W\", label=\"w [units]\", flow=False),\n",
    ")\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "s_data = np.random.normal(size=10_000) + np.ones(10_000)\n",
    "w_data = np.random.normal(size=10_000)\n",
    "s_data = np.random.normal(1, 0.5, 100_000)\n",
    "w_data = np.random.normal(13, 3, 100_000)\n",
    "\n",
    "h = Hist.new.Reg(50, -5, 5, name=\"S\", label=\"s [units]\", flow=False).Reg(50, -20, 20, name=\"W\", label=\"w [units]\", flow=False).Double().fill(\n",
    "    s_data, w_data)\n",
    "\n",
    "# plot2d_full\n",
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "h.plot2d_full(\n",
    "    main_cmap=\"coolwarm\",\n",
    "    top_ls=\"--\",\n",
    "    top_color=\"orange\",\n",
    "    top_lw=2,\n",
    "    side_ls=\":\",\n",
    "    side_lw=2,\n",
    "    side_color=\"steelblue\",\n",
    ")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we really are mixing all libraries - checking the projections with `mplhep`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mplhep\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(12, 5))\n",
    "\n",
    "mplhep.histplot(h.project(\"S\"), ax=axs[0])\n",
    "mplhep.histplot(h.project(\"W\"), ax=axs[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pull plots are very commonly used in HEP, e.g. to make sure that a toy study produces unbiased results and returns adequate uncertainties.\n",
    "`Hist` hence provides a powerful and flexible method for pull studies, called `.plot_pull()`. The method can be passed in a Callable object to be used to fit the plot.\n",
    "\n",
    "For the sake of argument, let's define a Gaussian function (following the `Hist` documentation):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The unumpy module of the uncertainties package deals with uncertainties in arrays\n",
    "from uncertainties import unumpy as unp\n",
    "\n",
    "def pdf(x, a=1 / np.sqrt(2 * np.pi), x0=0, sigma=1, offset=0):\n",
    "    exp = unp.exp if a.dtype == np.dtype(\"O\") else np.exp  #  'O' for a Python object\n",
    "    return a * exp(-((x - x0) ** 2) / (2 * sigma ** 2)) + offset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 8))\n",
    "\n",
    "# Fill a histogram with 500 entries\n",
    "hh = hist.Hist(\n",
    "    hist.axis.Regular(\n",
    "        50, -5, 5, name=\"X\", label=\"x [units]\", underflow=False, overflow=False\n",
    "    )\n",
    ").fill(np.random.normal(size=500))\n",
    "\n",
    "# Plot the histogram (data points) together with the fitted curve, and the pull plot\n",
    "main_ax_artists, subplot_ax_artists = hh.plot_pull(\n",
    "    pdf,  # Normal distributions are so commonly used in pull distributions that the str aliases \"normal\", \"gauss\", and \"gaus\" are supported as well!\n",
    "    eb_ecolor=\"steelblue\",\n",
    "    eb_mfc=\"steelblue\",\n",
    "    eb_mec=\"steelblue\",\n",
    "    eb_fmt=\"o\",\n",
    "    eb_ms=6,\n",
    "    eb_capsize=1,\n",
    "    eb_capthick=2,\n",
    "    eb_alpha=0.8,\n",
    "    fp_c=\"hotpink\",\n",
    "    fp_ls=\"-\",\n",
    "    fp_lw=2,\n",
    "    fp_alpha=0.8,\n",
    "    bar_fc=\"royalblue\",\n",
    "    pp_num=5,\n",
    "    pp_fc=\"royalblue\",\n",
    "    pp_alpha=0.5,\n",
    "    pp_ec=None,\n",
    "    ub_alpha=0.5,\n",
    ")\n",
    "\n",
    "plt.ylim((-2,2))\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}