{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Week 2 Notebook: Data Formats and Exploration\n", "===============================================================\n", "\n", "Particle physics uses a data format (and analysis ecosystem) based on the `ROOT` library: https://root.cern.ch/.\n", "`ROOT` is a framework for data processing, born at CERN, at the heart of the research on high-energy physics. \n", "A `ROOT` file is a compressed binary file in which we can save objects of any type.\n", "There are python bindings built into `ROOT`, so called `PyROOT`, but for now we won't discuss this.\n", "\n", "Recently, a different library called `uproot` has been developed allowing python users to do `ROOT` I/O directly: https://uproot.readthedocs.io/en/latest/.\n", "`uproot` is a reader and a writer of the `ROOT` file format using only `Python` and `Numpy`. \n", "Unlike the standard C++ `ROOT` implementation, uproot is only an I/O library, primarily intended to stream data into machine learning libraries in `Python`. \n", "Unlike `PyROOT` and `root_numpy`, uproot does not depend on C++ `ROOT`. \n", "Instead, it uses `Numpy` to cast blocks of data from the `ROOT` file as `Numpy` arrays.\n", "It can also make jagged or awkward arrays based on this library: https://github.com/scikit-hep/awkward-1.0 {cite:p}`Pivarski:2020qcb`.\n", "\n", "\n", "For a tutorial on `uproot`, see: https://hsf-training.github.io/hsf-training-uproot-webpage/ or here: https://github.com/jpivarski-talks/2020-07-13-pyhep2020-tutorial." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import uproot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Remote files\n", "One nice feature of `ROOT` and `uproot` is the ability to read a file remotely. \n", "This is done using the `XRootD` library: https://xrootd.slac.stanford.edu/.\n", "`XRootD` is a generic suite for fast, low latency and scalable data access, which can serve natively any kind of data, organized as a hierarchical filesystem-like namespace, based on the concept of directory. \n", "\n", "## XRootD filenames\n", "We specify that this is a remote file accessed via `XRootD` using the `root://` protocol in front of the filename. \n", "The redirector `eospublic.cern.ch` specifies that we're accessing data from the CERN open data repository. \n", "And finally the rest of the file name `/eos/opendata/cms/datascience/HiggsToBBNtupleProducerTool/HiggsToBBNTuple_HiggsToBB_QCD_RunII_13TeV_MC/train/ntuple_merged_10.root` specifies the exact file." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = uproot.open('root://eospublic.cern.ch//eos/opendata/cms/datascience/HiggsToBBNtupleProducerTool/HiggsToBBNTuple_HiggsToBB_QCD_RunII_13TeV_MC/train/ntuple_merged_10.root')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Trees\n", "One common object is a tree. \n", "Trees in `ROOT` are basically just tables of information. \n", "Trees are composed of branches, which are the columns of the table. \n", "The rows usually represent events (individual bunch crossings).\n", "However, in this case, each row represents a jet (a localized collection of particles within a single event).\n", "\n", "First we assign the tree to a variable (named `tree` here).\n", "We can see how many rows (jets) are contained in the tree, by checking its `num_entries`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tree = f['deepntuplizer/tree']\n", "print(tree.num_entries)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's list the contents (branches) of the tree." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tree.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A note on \"jaggedness\"\n", "\n", "Some branches are listed has having an interpretation of `asdtype('>f4')` while some others are listed as `asjagged(asdtype('>f4'), 10)`.\n", "The former means there is only one number per jet.\n", "The latter means there may a variable number per jet. \n", "\n", "First, let's get just look at non-jagged arrays, starting with the ground truth labels." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Returns a dictionary\n", "labels = tree.arrays(['label_QCD_b',\n", " 'label_QCD_bb', \n", " 'label_QCD_c', \n", " 'label_QCD_cc', \n", " 'label_QCD_others', \n", " 'label_H_bb',\n", " 'sample_isQCD'], \n", " entry_stop=20000,\n", " library='np')\n", "labels" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# label QCD: require the sample to be QCD and any of the QCD flavors\n", "label_QCD = labels['sample_isQCD'] * (labels['label_QCD_b'] + \\\n", " labels['label_QCD_bb'] + \\\n", " labels['label_QCD_c'] + \\\n", " labels['label_QCD_cc'] + \\\n", " labels['label_QCD_others'])\n", "\n", "# label Hbb\n", "label_Hbb = labels['label_H_bb']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What is the signal and background?\n", "\n", "Signal (Higgs) events for this task look like this: \n", "\n", "$$gg \\to X\\to HH \\to b\\bar{b}b\\bar{b}$$\n", "\n", "\n", "\n", "Background (QCD) events for this task look like this: \n", "\n", "$$gg \\to q \\bar{q} / gg / b\\bar{b} / gb\\bar{b} / ggg / \\cdots$$\n", "\n", "\n", "\n", "We will look at differentiating signal and background at the \"jet\" level." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# number of overlap jets\n", "sum(label_QCD*label_Hbb) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# fraction of jets with some truth label defined \n", "sum(label_QCD+label_Hbb)/len(label_QCD+label_Hbb)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# jet features, namely transverse momemntum (pt) and soft-drop mass (sdmass or msd)\n", "jet_features = tree.arrays(['fj_pt', \n", " 'fj_sdmass'],\n", " entry_stop=20000,\n", " library='np')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "plt.figure()\n", "\n", "plt.hist(jet_features['fj_pt'],weights=label_QCD,bins=np.linspace(0,4000,101),density=True,alpha=0.7,label='QCD')\n", "plt.hist(jet_features['fj_pt'],weights=label_Hbb,bins=np.linspace(0,4000,101),density=True,alpha=0.7,label='H(bb)')\n", "plt.xlabel(r'Jet $p_{T}$ [GeV]')\n", "plt.ylabel('Fraction of jets')\n", "plt.legend()\n", "\n", "plt.figure()\n", "\n", "plt.hist(jet_features['fj_sdmass'],weights=label_QCD,bins=np.linspace(0,300,101),density=True,alpha=0.7,label='QCD')\n", "plt.hist(jet_features['fj_sdmass'],weights=label_Hbb,bins=np.linspace(0,300,101),density=True,alpha=0.7,label='H(bb)')\n", "plt.xlabel(r'Jet $m_{SD}$ [GeV]')\n", "plt.ylabel('Fraction of jets')\n", "plt.legend()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Jagged arrys\n", "\n", "Now let's look at an jagged (or awkward) array, like those related to track features, where there can be a variable number of tracks per jet." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import awkward as ak\n", "track_features = tree.arrays(['track_pt',\n", " 'label_H_bb'],\n", " entry_stop=20000,\n", " library='ak')\n", "track_features" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note the difference between a \"flat\" array, where these a fixed number per jet (like 1 per jet in the case of jet properties) and a jagged array, where there are a variable number. \n", "\n", "We can demonstrate this by looking at the first jet in the dataset. As we'll see there are 21 tracks, each with their own pt." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jet_features['fj_pt'][0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "track_features['track_pt'][0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note, behind the scenes, jagged arrays are just like normal numpy arrays, except there's additional structure, retreivable from the `ak.num` function, which tells us the (variable) number of tracks per jet." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ak.num(track_features['track_pt'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Operations with jaggedness\n", "\n", "We can plot the full distibution of track pts. But what if we want to find the highest track pt per jet and plot only that?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure()\n", "plt.hist(ak.flatten(track_features['track_pt']),bins=np.linspace(0,100,101),density=True,alpha=0.7,label='All tracks')\n", "plt.hist(ak.max(track_features['track_pt'], axis=-1),bins=np.linspace(0,100,101),density=True,alpha=0.7,label=r'Max. $p_{T}$ track per jet')\n", "\n", "plt.xlabel(r'Track $p_{T}$')\n", "plt.ylabel('Fraction of tracks')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conversion to regular array\n", "\n", "Sometimes we want to turn a jagged array into a regular array to make it easier to accomodate into a machine learning algorithm (like a fully conneted or convolutional neural network). \n", "The simplest way to do this is to use zero-padding and truncation to \"cap\" the number of objects at some fixed value, and zero-pad if there are less objects.\n", "\n", "To do this, we can first plot how many tracks there are per jet and choose a reasonable number to cap." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure()\n", "plt.hist(ak.num(track_features['track_pt']),bins=np.linspace(0,100,101),density=True,alpha=0.7)\n", "\n", "plt.xlabel(r'Number of tracks')\n", "plt.ylabel('Fraction of jets')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, 60 seems to be a reasonable number." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "max_tracks = 60\n", "pad_value = 0\n", "a = ak.fill_none(ak.pad_none(track_features['track_pt'], max_tracks, clip=True, axis=-1), pad_value, axis=-1).to_numpy()\n", "print(a.shape)\n", "a" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To recreate the previous plot, we can do the following (note: we expect minor differences):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure()\n", "plt.hist(a[a>0],bins=np.linspace(0,100,101),density=True,alpha=0.7, label='All tracks')\n", "plt.hist(np.max(a,axis=-1),bins=np.linspace(0,100,101),density=True,alpha=0.7, label=r'Max. $p_{T}$ track per jet')\n", "\n", "plt.xlabel(r'Number of tracks')\n", "plt.ylabel('Fraction of jets')\n", "plt.legend()\n", "plt.show()\n", "print(a.nonzero())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.9.4" } }, "nbformat": 4, "nbformat_minor": 4 }