8S50x_PROJ2_runcells.ipynb

{"cells": [{"cell_type": "markdown", "id": "4c879cea", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<hr style=\"height: 1px;\">\n", "<i>This notebook was authored by the 8.S50x Course Team, Copyright 2022 MIT All Rights Reserved.</i>\n", "<hr style=\"height: 1px;\">\n", "<br>\n", "\n", "<h1>Project 2 - Part I: Measuring Properties of the W Boson from LHC Data</h1>\n"]}, {"cell_type": "markdown", "id": "88e7c942", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_2_0'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">PROJ2.0 Overview</h2>\n"]}, {"cell_type": "markdown", "id": "038120ad", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<h3>Navigation</h3>\n", "\n", "<table style=\"width:100%\">\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_2_1\">PROJ2.1 Introduction and Data Exploration</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_2_1\">PROJ2.1 Checkpoints</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_2_2\">PROJ2.2 Event Selection and Background Mitigation</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_2_2\">PROJ2.2 Checkpoints</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_2_3\">PROJ2.3 Beginning to Look for the W Signal in the Data</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_2_3\">PROJ2.3 Checkpoints</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_2_4\">PROJ2.4 Refining our Selection to Look for the W Signal in the Data</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_2_4\">PROJ2.4 Checkpoints</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_2_5\">PROJ2.5 Fit for W Peak</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_2_5\">PROJ2.5 Checkpoints</a></td>\n", "    </tr>\n", "</table>"]}, {"cell_type": "code", "execution_count": null, "id": "9e278ce5", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.0-slides\n", "\n", "from IPython.display import IFrame\n", "IFrame(src='https://mitx-8s50.github.io/slides/L09/slides_L09_09.html', width=970, height=550)"]}, {"cell_type": "markdown", "id": "0e495862", "metadata": {"tags": ["learner", "md", "catsoop_00", "learner_chopped"]}, "source": ["<h3>Data</h3>\n", "\n", ">description: Boosted Single Jet dataset at 8TeV<br>\n", ">source: https://zenodo.org/record/8035318 <br>\n", ">attribution: Philip Harris (CMS Collaboration), DOI:10.5281/zenodo.8035318 "]}, {"cell_type": "code", "execution_count": null, "id": "2bc7f81e", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.0-runcell00\n", "\n", "# NOTE: these files are too large to include in the original repository,\n", "# so you must download them using the options below\n", "#\n", "# Ways to download:\n", "#     1. Copy/paste the link (replace =0 with =1 to download automatically)\n", "#     2. Use the wget commands below (works in Colab, but you may need to install wget if using locally)\n", "#\n", "# Location of files:\n", "#     Move the files to the directory 'data'\n", "#\n", "# Using wget: (works in Colab)\n", "#     Upon downloading, the code below will move them to the appropriate directory\n", "\n", "#3GB Data Set: data1\n", "!wget https://www.dropbox.com/s/bcyab2lljie72aj/data.tgz?dl=0\n", "!mv data.tgz?dl=0 data.tgz #rename\n", "!tar -xvf data.tgz #extract the data\n", "!rm data.tgz #clean the downloaded file\n", "\n", "#130MB Data Set: data2\n", "!wget https://www.dropbox.com/s/p756oa4mfw17lfw/data.zip?dl=0\n", "!mv data.zip?dl=0 data.zip #rename\n", "!unzip data.zip #extract the data\n", "!rm data.zip #clean the downloaded file"]}, {"cell_type": "code", "execution_count": null, "id": "9f8e7998", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.0-runcell01\n", "\n", "# pre-requisites: install now if you have not already done so\n", "# uproot High energy physics python file format: https://masonproffitt.github.io/uproot-tutorial/aio.html\n", "!pip install uproot\n", "!pip install lmfit\n", "!pip install mplhep"]}, {"cell_type": "code", "execution_count": null, "id": "2fa4f7bc", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.0-runcell02\n", "\n", "import uproot\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import curve_fit\n", "import os,sys\n", "\n", "#!pip install lmfit #install lmfit if you have not done this already\n", "import lmfit as lm\n", "\n", "#!pip install mplhep #install mplhep if you have not done this already\n", "# plotting style for High Energy physics \n", "import mplhep as hep\n", "plt.style.use(hep.style.CMS)"]}, {"cell_type": "code", "execution_count": null, "id": "27cba66a", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.0-runcell03\n", "\n", "#set plot resolution\n", "%config InlineBackend.figure_format = 'retina'\n", "\n", "#set default figure parameters\n", "plt.rcParams['figure.figsize'] = (6,6)\n", "\n", "medium_size = 12\n", "large_size = 15\n", "\n", "plt.rc('font', size=medium_size)          # default text sizes\n", "plt.rc('xtick', labelsize=medium_size)    # xtick labels\n", "plt.rc('ytick', labelsize=medium_size)    # ytick labels\n", "plt.rc('legend', fontsize=medium_size)    # legend\n", "plt.rc('axes', titlesize=large_size)      # axes title\n", "plt.rc('axes', labelsize=large_size)      # x and y labels\n", "plt.rc('figure', titlesize=large_size)    # figure title\n"]}, {"cell_type": "markdown", "id": "5c0aa877", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_2_1'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">PROJ2.1 Introduction and Data Exploration</h2>    \n", "\n", "| [Top](#section_2_0) | [Previous Section](#section_2_0) | [Checkpoints](#problems_2_1) | [Next Section](#section_2_2) |\n"]}, {"cell_type": "markdown", "id": "a2e0e045", "metadata": {"tags": ["learner", "catsoop_01", "md", "learner_chopped"]}, "source": ["<h3> Loading data & Auxiliary functions</h3>\n", "\n", "Before we start, let's define collider coordinates centered around the collision point. We tend to write our momentum 4-vector as $\\vec{p}=(p_{T},\\eta,\\phi,m)$ in place of $\\vec{p}=(p,\\theta,\\phi,m)$. You can read more in this short link <a href=\"https://www.lhc-closer.es/taking_a_closer_look_at_lhc/0.momentum\" target=\"_blank\">here</a>.\n", "\n", "Now is a good time to look at the data. Let's take a look at the different samples we have. You should have run the first cell in Section 0 to download the data, but you can also download from these links:\n", "- <a href=\"https://www.dropbox.com/s/bcyab2lljie72aj/data.tgz?dl=0\" target=\"_blank\">3 GB file</a>\n", "- <a href=\"https://www.dropbox.com/s/p756oa4mfw17lfw/data.zip?dl=0\" target=\"_blank\">130 MB file</a>\n", "\n", "These files should be in a directory called `data`.\n"]}, {"cell_type": "code", "execution_count": null, "id": "697018b9", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.1-runcell01\n", "\n", "# Now let's open the data. \n", "\n", "# DATA\n", "#-------------------------------------------------------------------------------\n", "# Our data sample is the JetHT dataset. \n", "# What that means is the data passed triggers that have a jet in one of the triggers (discussed below).\n", "data   = uproot.open(\"data/JetHT_s.root\")[\"Tree\"]\n", "\n", "\n", "# SIMULATION\n", "#-------------------------------------------------------------------------------\n", "# In addition to above we have Monte Carlo Simulation of many processes\n", "# Some of these process are well modelled in simulation and some of them are not\n", "\n", "#the process qq=>W=>qq and qq=>Z=>qq processes at 8TeV collision energy\n", "wqq    = uproot.open(\"data/WQQ_s.root\")[\"Tree\"] \n", "zqq    = uproot.open(\"data/ZQQ_s.root\")[\"Tree\"]\n", "\n", "# To train NNs and make nice plots we will use larger samples produced at a different collision energy\n", "# qq=>W=>qq and qq=>Z=>qq at 13TeV collision energy\n", "wqq13  = uproot.open(\"data/skimh/WQQ_sh.root\")[\"Tree\"]\n", "zqq13  = uproot.open(\"data/skimh/ZQQ_sh.root\")[\"Tree\"]\n", "\n", "# qq=>W=>qq and qq=>Z=>qq at 8TeV collision energy\n", "wqq_n  = uproot.open(\"data/WQQ_8TeV_Jan11_r.root\")[\"Tree\"]\n", "zqq_n  = uproot.open(\"data/ZQQ_8TeV_Jan11_r.root\")[\"Tree\"]\n", "\n", "# Now we have our worst modeled background this is also our main background. \n", "# This is is our di-jet quark and gluon background. \n", "# We just call these backgrounds QCD because they are produced with Quantum Chromo Dynamics.\n", "qcd    = uproot.open(\"data/QCD_s.root\")[\"Tree\"]\n", "\n", "# Now we have the Higgs boson sample (we might need this in the future)\n", "ggh    = uproot.open(\"data/ggH.root\")[\"Tree\"]\n", "\n", "# And top-quark pair production background. \n", "tt     = uproot.open(\"data/TT.root\")[\"Tree\"]\n", "\n", "# Finally we have the rarer double W, W+Z and Z+Z diboson samples where we have two bosons instead of one\n", "ww     = uproot.open(\"data/WW.root\")[\"Tree\"]\n", "wz     = uproot.open(\"data/WZ.root\")[\"Tree\"]\n", "zz     = uproot.open(\"data/ZZ.root\")[\"Tree\"]\n", "\n", "dataDict = {'qcd': qcd,\n", "            'tt': tt,\n", "            'data': data,\n", "            'wqq': wqq,\n", "            'zqq': zqq,\n", "            'wqq13': wqq13,\n", "            'zqq13': zqq13,\n", "            'wqq_n': wqq_n,\n", "            'zqq_n': zqq_n,\n", "            'ww': ww,\n", "            'zz': zz,\n", "            'wz': wz,\n", "            'ggh': ggh\n", "            }\n", "from collections import OrderedDict \n", "\n", "order_of_keys = ['data','qcd','tt','ww','zz','wz','wqq','wqq13','wqq_n','zqq','zqq13','zqq_n','ggh']\n", "list_of_tuples = [(key, dataDict[key]) for key in order_of_keys]\n", "OrdDataDict = OrderedDict(list_of_tuples)"]}, {"cell_type": "code", "execution_count": null, "id": "6dedeedc", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.1-runcell02\n", "\n", "#these are the datasets that we are working with\n", "#keys: ['data','qcd','tt','ww','zz','wz','wqq','wqq13','wqq_n','zqq','zqq13','zqq_n','ggh']\n", "\n", "# You can view all of these variables within each dataset using the `keys` option\n", "print('wqq keys')\n", "print(wqq.keys())\n", "print()\n", "print('data keys')\n", "print(data.keys())\n", "#note the bdt varaibles at the end  of the `data` can be ignored,\n", "#this is an old deep learning training that we will not use for this study\n"]}, {"cell_type": "code", "execution_count": null, "id": "92e21ff3", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.1-runcell03\n", "\n", "# these are the standard weights\n", "weights=[1000*18300,\"puweight\",\"scale1fb\"]\n", "\n", "def get_weights(weights,mask,key):\n", "    # the first element of the list is the scaling weight\n", "    weight = weights[0]\n", "    # this needs to be divided by 1000 if the sample is wqq_n or zqq_n\n", "    if key=='wqq_n' or key=='zqq_n': \n", "        print('divide weight by 1000.')\n", "        weight /= 1000.\n", "    if key=='ggh': weight /= 1000. #maybe ggh too?\n", "    # now let's loop over the following weights\n", "    for i in range(1,len(weights)):\n", "        weight *= OrdDataDict[key].arrays(weights[i], library=\"np\")[weights[i]][mask]\n", "    return weight\n", "\n", "# For our samples with different collision energy (13 TeV) we need to perform a little hack on the cross section weight\n", "# so we normalize them to the number of events of the 8 TeV collision energy samples after a simple mask\n", "\n", "#This computes the integral of weighted events assuming a basic mask (see below details of this basic selection)\n", "def integral(iData,iWeights,iKey):\n", "    def selection(iData):\n", "        trigger = (iData.arrays('trigger', library=\"np\")[\"trigger\"].flatten() > 0) # trigger selection\n", "        jetpt   = (iData.arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten() > 400) # require jet pT above certain threshold\n", "        allcuts = np.logical_and.reduce([trigger,jetpt]) # apply both masks at the same time\n", "        return allcuts\n", "    mask_sel=selection(iData)\n", "    # get weights and take the integral and return it\n", "    weight = get_weights(iWeights,mask_sel,iKey)\n", "    return np.sum(weight)\n", "\n", "def scale(iData8TeV,iData13TeV,iWeights,iKey8TeV,iKey13TeV):\n", "    int_8TeV  = integral(iData8TeV,iWeights,iKey8TeV)\n", "    int_13TeV = integral(iData13TeV,iWeights,iKey13TeV)\n", "    print(\"Scale %s:\"%iKey13TeV,'ratio: ',int_8TeV/int_13TeV,' 8 TeV integral: ',int_8TeV,' 13 TeV integral: ',int_13TeV)\n", "    return int_8TeV/int_13TeV\n", "\n", "# we define this extra scaling number as:\n", "wscale=scale(wqq,wqq13,weights,'wqq','wqq13')\n", "zscale=scale(zqq,zqq13,weights,'zqq','zqq13')\n", "\n", "#w_nscale=scale(wqq,wqq_n,[18300,\"puweight\",\"scale1fb\"],'wqq','wqq_n')\n", "#z_nscale=scale(zqq,zqq_n,[18300,\"puweight\",\"scale1fb\"],'zqq','zqq_n')\n", "\n", "# Note: you could apply this weight function such as\n", "# qcd: get_weights(weights,qcd_mask,'qcd')\n", "# wqq_13: get_weights(weights,w_mask,'wqq_13')*wscale"]}, {"cell_type": "code", "execution_count": null, "id": "24be8cb5", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.1-runcell04\n", "\n", "# define some labels and colors\n", "labels = {'qcd': 'QCD',\n", "          'wqq': 'W',\n", "          'zqq': 'Z',\n", "          'wqq13': 'W (13 to 8 TeV)',\n", "          'zqq13': 'Z (13 to 8 TeV)',\n", "          'wqq_n': 'W new',\n", "          'zqq_n': 'Z new',\n", "          'tt': 'tt',\n", "          'ggh': 'H',\n", "          'zz': 'ZZ',\n", "          'ww': 'WW',\n", "          'wz': 'WZ',\n", "          'data': 'Data',\n", "         }\n", "colors = {'qcd': 'orange',\n", "          'wqq': 'royalblue',\n", "          'zqq': 'r',\n", "          'wqq13': 'cornflowerblue',\n", "          'zqq13': 'salmon',\n", "          'wqq_n': 'lightsteelblue',\n", "          'zqq_n': 'lightcoral',\n", "          'tt': 'green',\n", "          'ggh': 'cyan',\n", "          'zz': 'purple',\n", "          'ww': 'brown',\n", "          'wz': 'crimson',\n", "          'data': 'black',\n", "         }\n", "\n", "# build a plot to compare/stack histograms\n", "def histErr(iVar,iLabel,iBins,iMin,iMax,iSims,iMasks,iData=None,iMaskData=None,\n", "            iLabels=None,iColors=None,\n", "            iDensity=True,iStack=False,iWeights=None):\n", "    \n", "    fig, ax = plt.subplots(1,1,figsize=(6,6),dpi=80)\n", "\n", "    # first plot the simulated data - build arrays\n", "    if isinstance(iSims,dict): # if iSims is a dict\n", "        simhists = [x.arrays(iVar, library=\"np\")[iVar][iMasks[key]] for key,x in iSims.items()] \n", "    else: # if it's a list\n", "        simhists = [iSims[i].arrays(iVar, library=\"np\")[iVar][iMasks[i]] for i in range(0,len(iSims))]\n", "        \n", "    # define labels\n", "    plot_labels = iLabels\n", "    if iLabels is None:\n", "        plot_labels = [labels[lk] for lk in list(iSims.keys())] #labels\n", "    plot_colors = iColors\n", "    if iColors is None:\n", "        plot_colors = [colors[lk] for lk in list(iSims.keys())] # colors\n", "    \n", "    # build the histogram weights\n", "    hist_weights = None\n", "    if iWeights:\n", "        hist_weights = [get_weights(weights,iMasks[key],key) for key in iSims.keys()]\n", "        if 'wqq13' in key:\n", "            hist_weights *= wscale\n", "        if 'zqq13' in key:\n", "            hist_weights *= zscale\n", "        \n", "    htype = 'bar'\n", "    if not iStack: htype='step'\n", "        \n", "    _,bins,_ = plt.hist(simhists,\n", "                        color=plot_colors, label=plot_labels, weights=hist_weights,\n", "                        range=(iMin,iMax), bins=iBins, alpha=.6, histtype=htype, \n", "                        density=iDensity,stacked=iStack)\n", "    \n", "    # now include the data points (if any)\n", "    if iData:\n", "        data = iData.arrays(iVar, library=\"np\")[iVar][iMaskData]\n", "        counts, binEdges = np.histogram(data,bins=iBins,range=(iMin,iMax),density=iDensity)\n", "        yerr = np.sqrt(counts) # let's apply Poisson uncertainties\n", "        if iDensity: yerr /= np.sqrt(sum(iMaskData)*(binEdges[1]-binEdges[0]))\n", "        binCenters = (binEdges[1:]+binEdges[:-1])*.5\n", "        plt.errorbar(binCenters, counts, yerr=yerr,fmt=\"o\",c=\"k\",label=\"Data\", ms=3)\n", "    \n", "    #if iDensity:\n", "    #   plt.ylim(0,0.015)\n", "    \n", "    #plt.legend(prop={'size': 10})\n", "    plt.legend(loc=1)\n", "    plt.xlabel(iLabel)\n", "    if iDensity: plt.ylabel(\"Normalized Counts\") \n", "    else: plt.ylabel(\"Counts\")\n", "    plt.show()"]}, {"cell_type": "markdown", "id": "5dab0be6", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_2_1'></a>     \n", "\n", "| [Top](#section_2_0) | [Restart Section](#section_2_1) | [Next Section](#section_2_2) |\n"]}, {"cell_type": "markdown", "id": "e97fae6c", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.1.1</span>\n", "\n", "Let's consider the objectives of this project. What are we doing and why are we looking for jets? Select all options below that are relevant to understanding the objectives of this project:\n", "\n", "A) The goal of this project is to find signatures of the Higgs boson.\\\n", "B) The goal of this project is to find W and/or Z bosons that decay into quarks.\\\n", "C) Quarks leave showers of particles that we reconstruct as jets.\\\n", "D) When the momentum of a W or Z boson is high enough, quarks will decay into a single jet cone.\\\n", "E) Studying jets allows us to probe the strong interaction and investigate the properties of quarks and gluons.\\\n", "F) The study of jets helps us to search for new physics phenomena, such as the production of exotic particles or particles beyond the Standard Model."]}, {"cell_type": "markdown", "id": "f93918a6", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.1.2</span>\n", "\n", "Do all the data sets have the same keys? For instance, is there a difference between the keys in the `data` file compared to the simulation files? Explore the data (optionally complete the code below to return the difference between lists of keys).\n", "\n", "Choose the correct option:\n", "\n", "A) The data sets have the same keys.\\\n", "B) The `data` files contain the same keys as the simulation files, but also some additional information.\\\n", "C) The simulation files contain the same keys as the `data` file, but also some additional information.\\\n", "D) The data sets all contain different keys and, therefore, different types of information."]}, {"cell_type": "code", "execution_count": null, "id": "0d286dd4", "metadata": {"tags": ["py", "learner_chopped", "draft"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.1.2\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "#these are the datasets that we are working with\n", "#keys: ['data','qcd','tt','ww','zz','wz','wqq','wqq13','wqq_n','zqq','zqq13','zqq_n','ggh']\n", "\n", "#find difference between lists:\n", "def diff_lists(list1, list2):\n", "    return #YOUR CODE HERE\n", "\n", "print(diff_lists(data.keys(),zqq.keys()))"]}, {"cell_type": "markdown", "id": "eb467d08", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.1.3</span>\n", "\n", "Let's familiarize ourselves with the data a little more and practice extracting some features. Use `np.mean()` to find the average number of jets (`njets`) in the `wqq` dataset. Also, find the average number of b-tags (`nbtags`) detected in the `ggh` dataset.\n", "\n", "Report your answer as a list of two numbers with precision 1e-2: `[avg njets in wqq, avg nbtags in ggh]`\n"]}, {"cell_type": "code", "execution_count": null, "id": "bab2247a", "metadata": {"tags": ["py", "draft", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.1.3\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "print('Avg number of jets in wqq:', #YOUR CODE HERE)\n", "print('Avg number of b-tags in ggh:', #YOUR CODE HERE)"]}, {"cell_type": "markdown", "id": "2806cdc7", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.1.4</span>\n", "\n", "Let's look more closely at what the `get_weights` and `integral` functions are doing.\n", "\n", "It is important to understand that weights are how we estimate the expected number of events. For the weights that we use in this sample, we have 3 numbers:\n", "\n", "- the total luminosity of the data `18300*1000` in units of $\\mathrm{fb}^{-1}$\n", "- the weight to adjust for the beam intensity known as the pileup weight, `puweight`\n", "- the production cross section of the sample in units of fb, `scale1fb`\n", "\n", "Note that the cross section for events in samples can be different due to the way samples were produced.  To select events, we apply a mask. Effectively this is just a cut requiring a certain element of the dataset to behave a certain way. \n", "\n", "Given the information above, show that by using the `get_weights` command, we can get the same value as the `integral` function once we have applied the right mask. Complete the code below, then enter your answer as a list of two numbers with precision 1e-2: `[sum of weights, integral of weighted events]`\n"]}, {"cell_type": "code", "execution_count": null, "id": "14cb4235", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.1.4\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "#build a mask to select everything for a sample\n", "sample='qcd'\n", "test_mask      = (dataDict[sample].arrays('trigger', library=\"np\")[\"trigger\"].flatten() >= 0)\n", "test_jet       = (dataDict[sample].arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten() > 400) # require jet pT above certain threshold\n", "test_comb      = np.logical_and.reduce([test_mask,test_jet]) # apply both masks at the same time\n", "\n", "#The function get_weights() returns the events with the proper weights, after the masks are applied\n", "#print(get_weights(weights,test_comb,sample))\n", "\n", "print('Sum of weights:', #YOUR CODE HERE)\n", "print('Integrated result:', #YOUR CODE HERE) #hint: use the integral() function"]}, {"cell_type": "markdown", "id": "bb7b8b02", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_2_2'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">PROJ2.2 Event Selection and Background Mitigation</h2>    \n", "\n", "| [Top](#section_2_0) | [Previous Section](#section_2_1) | [Checkpoints](#problems_2_2) | [Next Section](#section_2_3) |\n"]}, {"cell_type": "code", "execution_count": null, "id": "7c513840", "metadata": {"tags": ["learner", "py", "catsoop_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.2-runcell01\n", "\n", "# Let's select some data (note that trigger can only be > 0)\n", "\n", "# First let's build masks on our data - these will be boolean arrays\n", "alldata      = (dataDict['data'].arrays('trigger', library=\"np\")[\"trigger\"].flatten() >= -1000000)\n", "triggerdata1 = (dataDict['data'].arrays('trigger', library=\"np\")[\"trigger\"].flatten() % 2 > 0) #let's require the lowest trigger jet pT > 320\n", "triggerdata2 = (dataDict['data'].arrays('trigger', library=\"np\")[\"trigger\"].flatten() % 4 > 1) #let's require one of our standard triggers (jet pT > 370 )\n", "\n", "# Now let's make a plot of the fat jet pt  \n", "# normalized\n", "histErr('vjet0_pt','Fat jet $p_T$ [GeV]',50,300,1e3,\n", "        [dataDict['data'],dataDict['data'],dataDict['data']],\n", "        [alldata,triggerdata1,triggerdata2],\n", "        iLabels=['all','$p_T$>320','$p_T$>370'],\n", "        iColors=['black','red','blue'],\n", "        iDensity=True,iStack=False,iWeights=None)\n", "\n", "# and without density\n", "histErr('vjet0_pt','Fat jet $p_T$ [GeV]',50,300,1e3,\n", "        [dataDict['data'],dataDict['data'],dataDict['data']],\n", "        [alldata,triggerdata1,triggerdata2],\n", "        iLabels=['all','$p_T$>320','$p_T$>370'],\n", "        iColors=['black','red','blue'],\n", "        iDensity=False,iStack=False,iWeights=None)\n", "\n", "#So you can see as you cut tighter, you get much less jets, but the data will be cleaner (I suggest triggerdata1)"]}, {"cell_type": "code", "execution_count": null, "id": "41c7b2ea", "metadata": {"tags": ["learner", "py", "catsoop_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.2-runcell02\n", "\n", "# First let's define a quick selection (a simple pT cut of 400 GeV and a 320 GeV trigger)\n", "def selection(iData):\n", "    #lets apply a trigger selection\n", "    trigger = (iData.arrays('trigger', library=\"np\")[\"trigger\"].flatten() > 0)\n", "    #Now lets require the jet pt to be above a threshold\n", "    jetpt   = (iData.arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten() > 400)\n", "    standard_trig = (iData.arrays('trigger', library=\"np\")[\"trigger\"].flatten() % 2 > 0) #lets require one of our standard triggers (jet pT > 320 )\n", "    # standard_trig = (iData.arrays('trigger', library=\"np\")[\"trigger\"].flatten() % 4 > 1) #lets require one of our standard triggers (jet pT > 370 )\n", "    allcuts = np.logical_and.reduce([trigger,jetpt])\n", "    return allcuts\n", "\n", "#print(wqq.arrays())\n", "# Let's look at all the data files (except the 8 TeV W and Z samples - let's work with the 13 TeV ones)\n", "myDataDict = OrdDataDict.copy()\n", "del myDataDict['wqq_n']\n", "del myDataDict['zqq_n']\n", "del myDataDict['data']\n", "\n", "# Get masks for the selection defined above (both for simulated datasets and data)\n", "masks = {}\n", "for key in myDataDict: masks[key] = selection(myDataDict[key])\n", "maskData = selection(dataDict['data'])\n", "\n", "# Now let's plot the mass and the groomed mass (msd0) for the QCD background\n", "fig, ax = plt.subplots(1,1,figsize=(6,6),dpi=80)\n", "plt.title(\"QCD Background\")\n", "plt.hist(qcd.arrays('vjet0_mass', library=\"np\")[\"vjet0_mass\"][masks['qcd']],weights=get_weights(weights,masks['qcd'],'qcd'),\n", "         bins=50,range=(0,300), color='salmon',label=\"groomed mass\", alpha=.6)\n", "plt.hist(qcd.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][masks['qcd']], weights=get_weights(weights,masks['qcd'],'qcd'),\n", "         bins=50,range=(0,300), color='red',label=\"mass\", alpha=.6)\n", "plt.legend()\n", "plt.xlabel(\"QCD Jet mass [GeV]\")\n", "plt.ylabel(\"Counts\")\n", "plt.show()\n", "\n", "# Let's look at the W/Z samples now (8 TeV collision energy)\n", "fig, ax = plt.subplots(1,1,figsize=(6,6),dpi=80)\n", "plt.title(\"8 TeV Collision Energy\")\n", "plt.hist(myDataDict['wqq'].arrays('vjet0_mass', library=\"np\")[\"vjet0_mass\"][masks['wqq']],weights=get_weights(weights,masks['wqq'],'wqq'),\n", "         bins=50,range=(0,300), color='salmon',label=\"W mass\", alpha=.6)\n", "plt.hist(myDataDict['wqq'].arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][masks['wqq']], weights=get_weights(weights,masks['wqq'],'wqq'),\n", "         bins=50,range=(0,300), color='red',label=\"W groomed mass\", alpha=.6)\n", "plt.hist(myDataDict['zqq'].arrays('vjet0_mass', library=\"np\")[\"vjet0_mass\"][masks['zqq']],weights=get_weights(weights,masks['zqq'],'zqq'),\n", "         bins=50,range=(0,300), color='pink',label=\"Z mass\", alpha=.6)\n", "plt.hist(myDataDict['zqq'].arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][masks['zqq']], weights=get_weights(weights,masks['zqq'],'zqq'),\n", "         bins=50,range=(0,300), color='hotpink',label=\"Z groomed mass\", alpha=.6)\n", "plt.legend()\n", "plt.xlabel(\"Signal Jet mass [GeV]\")\n", "plt.ylabel(\"Counts\")\n", "plt.show()\n", "\n", "# Let's look at the W/Z samples now (13 TeV collision energy)\n", "# Note that in the weights we need to multiply by wscale\n", "fig, ax = plt.subplots(1,1,figsize=(6,6),dpi=80)\n", "plt.title(\"13 TeV Collision Energy\")\n", "plt.hist(myDataDict['wqq13'].arrays('vjet0_mass', library=\"np\")[\"vjet0_mass\"][masks['wqq13']],weights=get_weights(weights,masks['wqq13'],'wqq13')*wscale,\n", "         bins=50,range=(0,300), color='salmon',label=\"W  mass\", alpha=.6)\n", "plt.hist(myDataDict['wqq13'].arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][masks['wqq13']], weights=get_weights(weights,masks['wqq13'],'wqq13')*wscale,\n", "         bins=50,range=(0,300), color='red',label=\"W groomed mass\", alpha=.6)\n", "plt.hist(myDataDict['zqq13'].arrays('vjet0_mass', library=\"np\")[\"vjet0_mass\"][masks['zqq13']],weights=get_weights(weights,masks['zqq13'],'zqq13')*wscale,\n", "         bins=50,range=(0,300), color='pink',label=\"Z mass\", alpha=.6)\n", "plt.hist(myDataDict['zqq13'].arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][masks['zqq13']], weights=get_weights(weights,masks['zqq13'],'zqq13')*wscale,\n", "         bins=50,range=(0,300), color='hotpink',label=\"Z groomed mass\", alpha=.6)\n", "plt.legend()\n", "plt.xlabel(\"Signal Jet mass [GeV]\")\n", "plt.ylabel(\"Counts\")\n", "plt.show()"]}, {"cell_type": "code", "execution_count": null, "id": "e506f4b7", "metadata": {"tags": ["learner", "py", "catsoop_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.2-runcell03\n", "\n", "# Compute the t21 ratio\n", "# let's use the same selection set above\n", "# note that here we are going to use our 13 TeV signal samples\n", "#print(len(masks[\"qcd\"]))\n", "\n", "fig, ax = plt.subplots(1,1,figsize=(6,6),dpi=80)\n", "qcdt21 = (qcd.arrays('vjet0_t2', library=\"np\")[\"vjet0_t2\"][masks['qcd']]/\n", "          qcd.arrays('vjet0_t1', library=\"np\")[\"vjet0_t1\"][masks['qcd']])\n", "wt21 = (wqq13.arrays('vjet0_t2', library=\"np\")[\"vjet0_t2\"][masks['wqq13']]/\n", "          wqq13.arrays('vjet0_t1', library=\"np\")[\"vjet0_t1\"][masks['wqq13']])\n", "\n", "plt.hist(qcdt21, weights=get_weights(weights,masks['qcd'],'qcd'),\n", "         bins=50, color='red',label=\"QCD\", alpha=.6, density=True)\n", "plt.hist(wt21, weights=get_weights(weights,masks['wqq13'],'wqq13')*wscale,\n", "         bins=50, color='black',label=\"W\", alpha=.6, density=True)\n", "plt.legend()\n", "plt.xlabel(r\"$\\tau_{21}$\")\n", "plt.ylabel(\"Normalized Counts\")\n", "plt.show()"]}, {"cell_type": "markdown", "id": "53e2374c", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_2_2'></a>     \n", "\n", "| [Top](#section_2_0) | [Restart Section](#section_2_2) | [Next Section](#section_2_3) |\n"]}, {"cell_type": "markdown", "id": "a64cd39a", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.2.1</span>\n", "\n", "Let's understand what the trigger is doing in code cell `PROJ2.2-runcell01`. We have defined two critical triggers that we care about. The first is whether an event has a transverse momentum pT > 320 GeV, and the second is whether an event has a transverse momentum of pT > 370 GeV.\n", "\n", "To characterize the trigger, the first bit is 1 if pT > 320 GeV and 0 if it's not. The second bit is 1 if pT > 370 GeV and 0 otherwise. We can write the bit value as: `trigger = 2*(pT > 370) + (pT > 320)`. Consider the following possible scenarios:\n", "\n", "- if we have an event with pT < 320, the value of triggger is 0\n", "- if we have an event with pT > 320 but less than 370, the value of the trigger is 1\n", "- if we have an event with pT > 370, the value of trigger is 3\n", "\n", "These are the only possible values of trigger. So, we can define the criteria for selecting events with pT > 320 GeV as `trigger % 2 > 0` (i.e., trigger mod 2 = 1).\n", "\n", "What is the criteria for selecting events with pT > 370? Complete the code below."]}, {"cell_type": "code", "execution_count": null, "id": "25983ffa", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.2.1\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "#this function is defined for you\n", "def pass_320(trigger):\n", "    #return 1 for trigger value 1 or 3\n", "    if trigger % 2 > 0:\n", "        return 1\n", "    else:\n", "        return 0\n", "    \n", "#this function you must complete\n", "def pass_370(trigger):\n", "    #return 1 for trigger value 3 only\n", "    if #YOUR CODE HERE:\n", "        return 1\n", "    else:\n", "        return 0\n", "    \n", "    \n", "#Run the following to check expected output\n", "#------------------------------------------\n", "\n", "def trigger_val(pass_320,pass_370):\n", "    #return trigger value based on given formula:\n", "    #trigger = 2*(pT > 370) + (pT > 320)\n", "    return 2*pass_370 + pass_320\n", "    \n", "\n", "#list the values of pass_320 and pass_370\n", "#consider the following combinations\n", "\n", "events1 = [0,0] #particle is below 320, trigger should return 0\n", "events2 = [1,0] #particle is between 320 and 370, trigger should return 1\n", "events3 = [1,1] #particle is between 320 and 370, trigger should return 3\n", "    \n", "def test_events(ievents):\n", "    pass_320_val, pass_370_val = ievents\n", "    trigger = trigger_val(pass_320_val, pass_370_val)\n", "    print('event', ievents)\n", "    print('trigger', trigger_val(pass_320_val, pass_370_val))\n", "    print('pass_320, computed',pass_320(trigger))\n", "    print('pass_320, expected',pass_320_val)\n", "    print('pass_370, computed',pass_370(trigger))\n", "    print('pass_370, expected',pass_370_val)\n", "    print()\n", "    return\n", "\n", "test_events(events1)\n", "test_events(events2)\n", "test_events(events3)"]}, {"cell_type": "markdown", "id": "622efd46", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.2.2</span>\n", "\n", "In this project, we are using simulated data for both 8 TeV and 13 TeV energies. It turns out that our data at 8 TeV provides the most accurate predictions, however, the 13 TeV distributions appear much smoother in the plots because they have more events. Smooth shapes, particularly of invariant quanities like mass, make it it easier to plot and interpolate.\n", "\n", "How can we effectively use both simulation data sets in our analysis?\n", "\n", "A) We can use the shape of the 13 TeV distributions, but scale the normalization to 8TeV distributions. This is an approximation, but it  gets the best features of both.\\\n", "B) We can't use 13 TeV at all, just 8 TeV for 8TeV data\\\n", "C) We can separately analyze the 8 TeV and 13 TeV datasets and compare the obtained results.\\\n"]}, {"cell_type": "markdown", "id": "763e9f37", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.2.3</span>\n", "\n", "Below what value of `t2/t1` are W events dominant? Enter your answer as number with precision 1e-1."]}, {"cell_type": "markdown", "id": "2c8023d4", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.2.4</span>\n", "\n", "Which of the following statements describes what will happen if we change our cut by increasing the `t2/t1` threshold? Select all that apply.\n", "\n", "A) We will get more W events compared to background.\\\n", "B) We will be able to better distinguish W signal from background.\\\n", "C) We will no longer be able to distinguish W signal from background.\\\n", "D) Nothing because the signal is independent of `t2/t1`.\n"]}, {"cell_type": "markdown", "id": "3c4d118d", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_2_3'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">PROJ2.3 Beginning to Look for the W Signal in the Data</h2>    \n", "\n", "| [Top](#section_2_0) | [Previous Section](#section_2_2) | [Checkpoints](#problems_2_3) | [Next Section](#section_2_4) |\n"]}, {"cell_type": "code", "execution_count": null, "id": "230af238", "metadata": {"tags": ["learner", "py", "catsoop_03", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.3-runcell01\n", "\n", "try:\n", "    del myDataDict['wqq'] #let's omit 8 TeV samples from here\n", "    del myDataDict['zqq']\n", "except:\n", "    print('samples already deleted')\n", "    \n", "# let's compare shapes\n", "histErr('vjet0_msd0','Fat jet $m_{SD}$ [GeV]',50,40,200,\n", "        myDataDict,masks,\n", "        dataDict['data'],maskData,\n", "        iDensity=True,iStack=False,iWeights=True)\n", "\n", "# Let's do a stacked plot  of all simulation and data\n", "histErr('vjet0_msd0','Fat jet $m_{SD}$ [GeV]',50,40,200,\n", "        myDataDict,masks,\n", "        dataDict['data'],maskData,\n", "        iDensity=False,iStack=True,iWeights=True)"]}, {"cell_type": "markdown", "id": "1219215a", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.3.1</span>\n", "\n", "When we look at the above distributions, we see that the W, Z and other channels yield resonant bumps at various masses. However, in the bottom plot, we don't see these bumps in the data or MC simulation. Why do we not see them? \n", "\n", "A) The bumps are not visible in the bottom plot because the detector effects smear out the resonant structures.\\\n", "B) The bumps are not observed in the data or MC simulation in the bottom plot due to limitations in the modeling of certain physical processes.\\\n", "C) The bumps seen in the W, Z, and other channels might be due to statistical fluctuations or specific experimental conditions, which are not replicated in the bottom plot.\\\n", "D) The bumps are there in the plot on the bottom, but the QCD background is just so much larger than the W, Z samples and others that we just can't see them."]}, {"cell_type": "code", "execution_count": null, "id": "615a6b58", "metadata": {"tags": ["learner", "py", "learner_chopped", "catsoop_03"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.3-runcell02\n", "\n", "#Load the data, if you have not done so in Section 1\n", "\n", "wqq    = uproot.open(\"data/WQQ_s.root\")[\"Tree\"]\n", "zqq    = uproot.open(\"data/ZQQ_s.root\")[\"Tree\"]\n", "wqq13  = uproot.open(\"data/skimh/WQQ_sh.root\")[\"Tree\"]\n", "zqq13  = uproot.open(\"data/skimh/ZQQ_sh.root\")[\"Tree\"]\n", "wqq_n  = uproot.open(\"data/WQQ_8TeV_Jan11_r.root\")[\"Tree\"]\n", "zqq_n  = uproot.open(\"data/ZQQ_8TeV_Jan11_r.root\")[\"Tree\"]\n", "qcd    = uproot.open(\"data/QCD_s.root\")[\"Tree\"]\n", "tt     = uproot.open(\"data/TT.root\")[\"Tree\"]\n", "ww     = uproot.open(\"data/WW.root\")[\"Tree\"]\n", "wz     = uproot.open(\"data/WZ.root\")[\"Tree\"]\n", "zz     = uproot.open(\"data/ZZ.root\")[\"Tree\"]\n", "ggh    = uproot.open(\"data/ggH.root\")[\"Tree\"]\n", "data   = uproot.open(\"data/JetHT_s.root\")[\"Tree\"]"]}, {"cell_type": "code", "execution_count": null, "id": "c2319e9a", "metadata": {"tags": ["learner", "py", "learner_chopped", "catsoop_03"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.3-runcell03\n", "\n", "def selection(iData):\n", "    '''\n", "    Standard pre-selection\n", "    '''\n", "    #lets apply a trigger selection\n", "    trigger = (iData.arrays('trigger', library=\"np\")[\"trigger\"].flatten() > 0)\n", "\n", "    #Now lets require the jet pt to be above a threshold (400 TODO: ASK about units)\n", "    jetpt   = (iData.arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten() > 400)\n", "\n", "    #Lets apply both jetpt and trigger at the same time\n", "    #standard_trig = (iData.arrays('trigger', library=\"np\")[\"trigger\"].flatten() % 4 > 1) #lets require one of our standard triggers (jet pT > 370 )\n", "    allcuts = np.logical_and.reduce([trigger,jetpt])\n", "\n", "    return allcuts\n", "    \n", "def get_weights(iData,weights,sel):\n", "    \n", "    weight = weights[0]\n", "    \n", "    for i in range(1,len(weights)):\n", "        weight *= iData.arrays(weights[i],library=\"np\")[weights[i]][sel]\n", "        \n", "    return weight\n", "\n", "def integral(iData,iWeights):\n", "    '''\n", "    This computs the integral of weighted events \n", "    assuming a selection given by the function selection (see below)\n", "    '''\n", "    \n", "    #perform a selection on the data (\n", "    mask_sel=selection(iData)\n", "    \n", "    #now iterate over the weights not the weights are in the format of [number,variable name 1, variable name 2,...]\n", "    weight  =iWeights[0]\n", "    \n", "    for i0 in range(1,len(iWeights)):\n", "        weightarr = iData.arrays(iWeights[i0], library=\"np\")[iWeights[i0]][mask_sel].flatten()\n", "        \n", "        #multiply the weights\n", "        weight    = weight*weightarr\n", "    \n", "    #now take the integral and return it\n", "    return np.sum(weight)\n", "\n", "\n", "def scale(iData8TeV,iData13TeV,iWeights):\n", "    '''\n", "    This computes the integral of two selections for two datasets labelled 8TeV and 13TeV,\n", "    but really can be 1 and 2. Then it returns the ratio of the integrals\n", "    '''\n", "    \n", "    int_8TeV  = integral(iData8TeV,iWeights)\n", "    int_13TeV = integral(iData13TeV,iWeights)\n", "    \n", "    return int_8TeV/int_13TeV"]}, {"cell_type": "code", "execution_count": null, "id": "751587d0", "metadata": {"tags": ["learner", "py", "learner_chopped", "catsoop_03"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.3-runcell04\n", "\n", "def plotDataSim(iVar, iSelection, iVarName, iRange):\n", "    \n", "    #Lets Look at the mass\n", "    weights = [1000*18300, \"puweight\", \"scale1fb\"]\n", "    mrange = iRange #range for mass histogram [GeV]\n", "    bins=40            #bins for mass histogram\n", "    density = False     #to plot the histograms as a density (integral=1)\n", "\n", "    qcdsel      = iSelection(qcd)\n", "    wsel        = iSelection(wqq13)\n", "    zsel        = iSelection(zqq13)\n", "    datasel     = iSelection(data)\n", "    ttsel       = iSelection(tt)\n", "    wwsel       = iSelection(ww)\n", "    wzsel       = iSelection(wz)\n", "    zzsel       = iSelection(zz)\n", "    gghsel      = iSelection(ggh)\n", "\n", "    wscale=scale(wqq,wqq13,weights)\n", "    zscale=scale(zqq,zqq13,weights)\n", "\n", "    # Getting the masses of selected events\n", "    dataW = data.arrays(iVar, library=\"np\") [iVar][datasel]\n", "    qcdW  = qcd.arrays(iVar, library=\"np\")  [iVar][qcdsel]\n", "    wW    = wqq13.arrays(iVar, library=\"np\")[iVar][wsel]\n", "    zW    = zqq13.arrays(iVar, library=\"np\")[iVar][zsel]\n", "    zzW   = zz   .arrays(iVar, library=\"np\")[iVar][zzsel]\n", "    wzW   = wz   .arrays(iVar, library=\"np\")[iVar][wzsel]\n", "    wwW   = ww   .arrays(iVar, library=\"np\")[iVar][wwsel]\n", "    ttW   = tt   .arrays(iVar, library=\"np\")[iVar][ttsel]\n", "    gghW  = ggh  .arrays(iVar, library=\"np\")[iVar][gghsel]\n", "\n", "    #Define the weights for the histograms\n", "    hist_weights = [get_weights(qcd,weights,qcdsel),\n", "                    get_weights(wqq13,weights,wsel)*wscale,\n", "                    get_weights(zqq13,weights,zsel)*zscale,\n", "                    get_weights(zz,weights,zzsel),\n", "                    get_weights(wz,weights,wzsel),\n", "                    get_weights(ww,weights,wwsel),\n", "                    get_weights(tt,weights,ttsel),\n", "                   ]\n", "\n", "    #Hint: Provide a list of selected data\n", "    plt.hist([qcdW,wW, zW, zzW, wzW, wwW, ttW],\n", "             color=[\"royalblue\",\"r\", \"orange\",\"g\", \"b\", \"purple\", \"cyan\",], \n", "             label=[\"QCD\", \"W\", \"Z\", \"ZZ\", \"WZ\", \"WW\", \"tt\",], \n", "             weights=hist_weights,\n", "             range=mrange, bins=50, alpha=.6, density=density,stacked=True)\n", "\n", "    #Other configurations for the histogram\n", "    counts, bins = np.histogram(dataW, bins=bins, range=mrange, density=density)\n", "    yerr = np.sqrt(counts) / np.sqrt(len(dataW)*np.diff(bins))\n", "    binCenters = (bins[1:]+bins[:-1])*.5\n", "    plt.errorbar(binCenters, counts, yerr=yerr,fmt=\"o\",c=\"k\",label=\"data\")\n", "    plt.legend()\n", "    plt.xlabel(iVarName)\n", "    plt.ylabel(\"Counts\")\n", "    plt.show()\n", "\n", "plotDataSim(\"vjet0_msd0\", selection, \"Jet Mass\", [40,200])\n", "plotDataSim(\"vjet0_t2\", selection, r\"$\\tau_2$\", [0,0.5]) \n", "#Add some code here to compare variables"]}, {"cell_type": "code", "execution_count": null, "id": "a2756213", "metadata": {"tags": ["learner", "py", "learner_chopped", "catsoop_03"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.3-runcell06\n", "\n", "#Define tau2/tau1\n", "def t21_func(itau1,itau2):\n", "    return itau2/itau1\n", "\n", "\n", "def selectionW_firstcut(iData):\n", "    '''\n", "    This is the specific selection for selecting out events with W signal for our analysis\n", "    '''\n", "    \n", "    #Pre-selection citeria\n", "    trigger = (iData.arrays('trigger', library=\"np\")[\"trigger\"].flatten() >= 0)\n", "    jetpt   = (iData.arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten() >= 400)\n", "    \n", "    #Select the jets to compute tau2/tau1\n", "    jett2   = (iData.arrays('vjet0_t2', library=\"np\")[\"vjet0_t2\"].flatten())\n", "    jett1   = (iData.arrays('vjet0_t1', library=\"np\")[\"vjet0_t1\"].flatten())\n", "        \n", "    t21 = t21_func(jett1,jett2)\n", "                                \n", "    #And then perform the cut\n", "    #Hint: You could determine the threshold of the cut by plotting the distribution of \n", "    #t21ddt scores for W and background and then determine a ball park threshold\n", "    #where you think the W signal would be best selected\n", "    #Or more simply you could look at the given plot and determine the appropriate threshold.\n", "    \n", "    t21cut   = t21 < #YOUR CODE HERE (enter t21 threshold)\n", "    \n", "    allcuts = np.logical_and.reduce([trigger, jetpt, t21cut])\n", "    \n", "    return allcuts\n", "\n", "plotDataSim(\"vjet0_msd0\", selectionW_firstcut, \"Jet Mass\",[40,200])"]}, {"cell_type": "markdown", "id": "a7f2fe5a", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_2_3'></a>   \n", "\n", "| [Top](#section_2_0) | [Restart Section](#section_2_3) |\n"]}, {"cell_type": "markdown", "id": "a409dd23", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.3.2</span>\n", "\n", "In the above plot, there is the simulation prediction (colored histograms), and then the data (black points). You can see the data has a shape that looks like there are two bumps on each other. However, the bumps are merged and peak roughly at the same spot. Why can we not just use the simulation to extract the W and Z bumps? Select all that apply:\n", "\n", "A) We can! The best way to analyze the properties of a signal is to see where is matches simulation exactly.\\\n", "B) We cannot because our simulations use assumptions that would bias our measurement. This is effectively a kind of circular analysis.\\\n", "C) While the simulation provides a useful reference, it is essential to account for potential discrepancies between the simulation and data due to uncertainties in the theoretical models, calibration of the detectors, or unknown physics phenomena. \n"]}, {"cell_type": "markdown", "id": "4a1cbf4c", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_2_4'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">PROJ2.4 Refining our Selection to Look for the W Signal in the Data</h2>   \n", "\n", "| [Top](#section_2_0) | [Previous Section](#section_2_3) | [Checkpoints](#problems_2_4) | [Next Section](#section_2_5) |\n"]}, {"cell_type": "markdown", "id": "a504fb6f", "metadata": {"tags": ["learner", "md", "learner_chopped", "catsoop_04"]}, "source": ["<h3>Objective 1: Develop a procedure to select data and make a W boson mass plot</h3>\n", "\n", "Since finding W peak is hard, we need to use another parameter, $\\rho$, which is a scaling variable for QCD jets. This parameter adds another channels of mass and $p_T$ to our selection, helping us to refine our W peak. The parameter $\\rho$ is defined in <a href=\"https://arxiv.org/pdf/1603.00027.pdf\" target=\"_blank\">this paper.</a>\n", "\n", "Your first goal is to figure out how $\\rho$ is defined by quoting the paper, and then figure out the best selections based on a combination of $\\rho$ and $\\tau_2/\\tau_1$. The final cut is based on a parameter defined as *DT* (Deccorelated Taggers) score:\n", "\n", "$$(\\tau_2/\\tau_1)_{dt} = \\tau_2/\\tau_1 - (\\text{your correlation})*\\rho$$\n", "\n", "Where the correlation `(your correlation)` is the correlation coefficient between $\\tau_2/\\tau_1$ and $\\rho$. \n", "\n", "To figure out the correlation, let's plot $\\tau_2/\\tau_1$ and $\\rho$ in the data first!"]}, {"cell_type": "markdown", "id": "72df7b03", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.4.1</span>\n", "\n", "Complete the code below to plot $\\tau_2/\\tau_1$ vs. $\\rho$. Specifically, we will check the function `rho_func` in the answer-checker, then you should use your result within the funciton `plot_taus_and_rho` to create a plot."]}, {"cell_type": "code", "execution_count": null, "id": "844184eb", "metadata": {"tags": ["py", "learner_chopped", "draft"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.4.1\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "def rho_func(imass,ipt,mu=1.):\n", "    return #YOUR CODE HERE\n", "    \n", "\n", "def plot_taus_and_rho(iData):\n", "    \n", "    jetptnocut = (iData.arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten())\n", "    jetmass = (iData.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"].flatten())\n", "    jett2   = (iData.arrays('vjet0_t2', library=\"np\")[\"vjet0_t2\"].flatten())\n", "    jett1   = (iData.arrays('vjet0_t1', library=\"np\")[\"vjet0_t1\"].flatten())\n", "    \n", "    #Define rho according to the paper (eqn. 3.2)\n", "    #Really this is rho_prime as defined in the paper, with mu=1 in these units\n", "    rho = #YOUR CODE HERE\n", "    \n", "    #Define tau2/tau1\n", "    t21 = t21_func(jett1,jett2)\n", "    \n", "    plt.hist2d(rho, t21, bins = 40)\n", "    \n", "    plt.xlabel(r\"$\\rho$\")\n", "    plt.ylabel(r\"$\\tau_2/\\tau_1$\")\n", "    plt.show()\n", "    \n", "plot_taus_and_rho(qcd)"]}, {"cell_type": "code", "execution_count": null, "id": "1807e3dd", "metadata": {"tags": ["learner", "py", "learner_chopped", "catsoop_04"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.4-runcell02\n", "\n", "def fit_correlation(iData):\n", "    \n", "    jetptnocut = (iData.arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten())\n", "    jetmass = (iData.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"].flatten())\n", "    jett2   = (iData.arrays('vjet0_t2', library=\"np\")[\"vjet0_t2\"].flatten())\n", "    jett1   = (iData.arrays('vjet0_t1', library=\"np\")[\"vjet0_t1\"].flatten())\n", "    \n", "    #Define rho according to the paper\n", "    rho = rho_func(jetmass,jetptnocut)\n", "    \n", "    #Define tau2/tau1\n", "    t21 = t21_func(jett1,jett2)\n", "    \n", "    plt.hist2d(rho, t21, bins = 40)\n", "    plt.xlabel(r\"$\\rho$\")\n", "    plt.ylabel(r\"$\\tau_2/\\tau_1$\")\n", "    \n", "    #Fit the line\n", "    #Produce 2D histogram\n", "    H,xedges,yedges = np.histogram2d(rho,t21, bins=40,density = True)\n", "    \n", "    bin_centers_x = (xedges[:-1]+xedges[1:])/2.0\n", "    bin_centers_y = (yedges[:-1]+yedges[1:])/2.0\n", "    \n", "    #Find the non-zero indicies\n", "    non_zero_idx = np.argwhere(H > 0.6) #You can play around with this!\n", "    x_idx = non_zero_idx[:,0]\n", "    y_idx = non_zero_idx[:,1]\n", "    \n", "    x_coord = [bin_centers_x[x_idx[i]] for i in range(0,len(x_idx))]\n", "    y_coord = [bin_centers_y[y_idx[i]] for i in range(0,len(y_idx))]\n", "    \n", "    #Fit a linear model on the points plotted\n", "    def func(x, a, b):\n", "        return a * x + b\n", "    plt.scatter(x_coord, y_coord)\n", "    \n", "    popt, pcov = curve_fit(func, x_coord, y_coord)\n", "    plt.plot(bin_centers_x, func(bin_centers_x, *popt), 'r-',\n", "             label='fit: a=%5.3f, b=%5.3f' % tuple(popt))\n", "    \n", "    #Show the fit result\n", "    legend = plt.legend()\n", "    plt.setp(legend.get_texts(), color='w')\n", "    plt.show()\n", "\n", "from scipy.optimize import curve_fit\n", "fit_correlation(qcd)"]}, {"cell_type": "markdown", "id": "b1ae5470", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.4.2</span>\n", "\n", "Now use the correlation from the plot to define $(\\tau_2/\\tau_1)_{dt}$ and plot it with $\\rho$ to verify that we have successfully decorrelated the tagger. If you do it correctly, you can see that the decorrelated scores are now independent of $\\rho$ (you should see a straight-line distribution in the histogram)!\n", "\n", "Specifically, complete the function `t21ddt_func`, which should decorrelate the `t21` value as a function of `rho`. Set the default value of `iMcorr` based on your fit above."]}, {"cell_type": "code", "execution_count": null, "id": "743daa54", "metadata": {"tags": ["py", "learner_chopped", "draft"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.4.2\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "def t21ddt_func(it21,irho,iMcorr=#YOUR CODE HERE):\n", "    #iMcorr is the correlation coefficient\n", "    return #YOUR CODE HERE\n", "\n", "\n", "def plot_tausdt_and_rho(iData):\n", "    \n", "    jetptnocut = (iData.arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten())\n", "    jetmass = (iData.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"].flatten())\n", "    jett2   = (iData.arrays('vjet0_t2', library=\"np\")[\"vjet0_t2\"].flatten())\n", "    jett1   = (iData.arrays('vjet0_t1', library=\"np\")[\"vjet0_t1\"].flatten())\n", "    \n", "    #Define rho according to the paper\n", "    rho = rho_func(jetmass,jetptnocut)\n", "    \n", "    #Define tau2/tau1\n", "    t21 = t21_func(jett1,jett2)\n", "    \n", "    #decorrelated tagger score\n", "    t21ddt = t21ddt_func(t21,rho)\n", "    \n", "    plt.hist2d(rho, t21ddt, bins = 40)\n", "    \n", "    plt.xlabel(r\"$\\rho$\")\n", "    plt.ylabel(r\"$\\tau_2/\\tau_1$_dt\")\n", "                \n", "    #Fit the line\n", "    #Produce 2D histogram\n", "    H,xedges,yedges = np.histogram2d(rho,t21ddt, bins=40,density = True)\n", "    \n", "    bin_centers_x = (xedges[:-1]+xedges[1:])/2.0\n", "    bin_centers_y = (yedges[:-1]+yedges[1:])/2.0\n", "    \n", "    #Find the non-zero indicies\n", "    non_zero_idx = np.argwhere(H > 0.6) #You can play around with this!\n", "    x_idx = non_zero_idx[:,0]\n", "    y_idx = non_zero_idx[:,1]\n", "    x_coord = [bin_centers_x[x_idx[i]] for i in range(0,len(x_idx))]\n", "    y_coord = [bin_centers_y[y_idx[i]] for i in range(0,len(y_idx))]\n", "    \n", "    #Fit a linear model on the points plotted\n", "    def func(x, a, b):\n", "        return a * x + b\n", "    plt.scatter(x_coord, y_coord)\n", "    \n", "    popt, pcov = curve_fit(func, x_coord, y_coord)\n", "    plt.plot(bin_centers_x, func(bin_centers_x, *popt), 'r-',\n", "             label='fit: a=%5.3f, b=%5.3f' % tuple(popt))\n", "    \n", "    #Show the fit result\n", "    legend = plt.legend()\n", "    plt.setp(legend.get_texts(), color='w')\n", "    plt.show()\n", "    \n", "plot_tausdt_and_rho(qcd)"]}, {"cell_type": "markdown", "id": "b18695e4", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.4.3</span>\n", "\n", "Since you figured out your decorrelation, determine the best cut for the decorrelated taggers score (Refer to the plot from code cell `PROJ2.2-runcell03` and the results from `Checkpoint 2.2.3`) and use it in your selection function!\n", "\n", "Complete the function `get_t21cut_W`, which should return values of `t21ddt` below the threshold that you define."]}, {"cell_type": "code", "execution_count": null, "id": "2bf017cf", "metadata": {"tags": ["py", "learner_chopped", "draft"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.4.3\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "def get_t21cut_W(t21ddt,t21_thresh=#YOUR CODE HERE):\n", "    #return values of t21ddt that occur below the threshold\n", "    #based on previous analysis of t21 vs. rho\n", "    ### YOUR CODE HERE ### \n", "    return\n", "\n", "    \n", "def selectionW(iData):\n", "    '''\n", "    This is the specific selection for selecting out events with W signal for our analysis\n", "    '''\n", "    \n", "    #Pre-selection citeria\n", "    trigger = (iData.arrays('trigger', library=\"np\")[\"trigger\"].flatten() >= 0)\n", "    jetpt   = (iData.arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten() >= 400)\n", "    jetptnocut = (iData.arrays('vjet0_pt', library=\"np\")[\"vjet0_pt\"].flatten())\n", "    jetmass = (iData.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"].flatten())\n", "    jett2   = (iData.arrays('vjet0_t2', library=\"np\")[\"vjet0_t2\"].flatten())\n", "    jett1   = (iData.arrays('vjet0_t1', library=\"np\")[\"vjet0_t1\"].flatten())\n", "    \n", "    \n", "    #Define the parameters rho, tau2/tau1\n", "    rho = rho_func(jetmass,jetptnocut)\n", "    t21 = t21_func(jett1,jett2)\n", "    \n", "    #Define the decorrelated tagger\n", "    Mcorr = #YOUR CODE HERE\n", "    t21ddt = t21ddt_func(t21,rho,Mcorr) \n", "                            \n", "    #And then perform the cut\n", "    #Hint: You could determine the threshold of the cut by plotting the distribution of \n", "    #t21ddt scores for W and background and then determine a ball park threshold\n", "    #where you think the W signal would be best selected\n", "    t21cut   = ### YOUR CODE HERE ### \n", "    \n", "    allcuts = np.logical_and.reduce([trigger, jetpt, t21cut])\n", "    \n", "    return allcuts"]}, {"cell_type": "code", "execution_count": null, "id": "9b8e37c1", "metadata": {"tags": ["learner", "py", "learner_chopped", "catsoop_04"]}, "outputs": [], "source": ["#>>>RUN: PROJ2.4-runcell05\n", "\n", "#Lets Look at the mass\n", "weights = [1000*18300, \"puweight\", \"scale1fb\"]\n", "mrange = (45,200)  #range for mass histogram [GeV]\n", "bins=40            #bins for mass histogram\n", "density = True     #to plot the histograms as a density (integral=1)\n", "\n", "qcdsel      = selectionW(qcd)\n", "wsel        = selectionW(wqq13)\n", "zsel        = selectionW(zqq13)\n", "datasel     = selectionW(data)\n", "ttsel       = selectionW(tt)\n", "wwsel       = selectionW(ww)\n", "wzsel       = selectionW(wz)\n", "zzsel       = selectionW(zz)\n", "gghsel      = selectionW(ggh)\n", "wscale=scale(wqq,wqq13,weights)\n", "zscale=scale(zqq,zqq13,weights)\n", "\n", "dataW = data.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][datasel]\n", "qcdW = qcd.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][qcdsel]\n", "wW = wqq13.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][wsel]\n", "zW = zqq13.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][zsel]\n", "zzW = zz.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][zzsel]\n", "wzW = wz.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][wzsel]\n", "wwW = ww.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][wwsel]\n", "ttW = tt.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][ttsel]\n", "gghW = ggh.arrays('vjet0_msd0', library=\"np\")[\"vjet0_msd0\"][gghsel]\n", "\n", "hist_weights = [get_weights(qcd,weights,qcdsel),\n", "                get_weights(wqq13,weights,wsel)*wscale,\n", "                get_weights(zqq13,weights,zsel)*zscale,\n", "                get_weights(zz,weights,zzsel),\n", "                get_weights(wz,weights,wzsel),\n", "                get_weights(ww,weights,wwsel),\n", "                get_weights(tt,weights,ttsel),\n", "               ]\n", "\n", "plt.hist([qcdW,wW, zW, zzW, wzW, wwW, ttW], \n", "         color=[\"royalblue\",\"r\", \"orange\",\"g\", \"b\", \"purple\", \"cyan\",], \n", "         label=[\"QCD\", \"W\", \"Z\", \"ZZ\", \"WZ\", \"WW\", \"tt\",], \n", "         weights=hist_weights,\n", "         range=mrange, bins=50, alpha=.6, density=density,stacked=True)\n", "\n", "counts, bins = np.histogram(dataW, bins=bins, range=mrange, density=density)\n", "yerr = np.sqrt(counts) / np.sqrt(len(dataW)*np.diff(bins))\n", "binCenters = (bins[1:]+bins[:-1])*.5\n", "plt.errorbar(binCenters, counts, yerr=yerr,fmt=\"o\",c=\"k\",label=\"data\")\n", "plt.legend()\n", "plt.xlabel(r\"Mass [GeV]\")\n", "plt.ylabel(\"Normalized Counts\")\n", "plt.show()"]}, {"cell_type": "markdown", "id": "4163ab64", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_2_5'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">PROJ2.5 Fit for W Peak</h2>   \n", "\n", "| [Top](#section_2_0) | [Previous Section](#section_2_4) | [Checkpoints](#problems_2_5) | [Next Section](#section_2_6) |\n"]}, {"cell_type": "markdown", "id": "3765cbff", "metadata": {"tags": ["learner", "md", "catsoop_05", "learner_chopped"]}, "source": ["<h3>Objective 2: Fit the W mass peak with the appropriate function</h3>\n", "\n", "Now you will perform the fit on W signal. It involves a few steps:\n", "    \n", "<h4>1. Defining a model</h4>\n", "\n", "First you need to define a fit model of your own. In this case we would use some functions (gaussian, exponential) in conjuntion with a 6th order polynomial. You could see more on how the order of the polynomials are determined here: https://en.wikipedia.org/wiki/Chow_test. The concepts were also covered in previous Lessons, if you want to review them. Adding a chow test will likely allow you to improve the measurement by lowering the polynomial. Whilte its not needed here, we strongly encourage this investigation. \n", "\n", "In the extended projec you will have the chance to determine the order of the polynomial for the Z fit. You might see that we might not necessarily need a 6th order polynomial for the Z fit. The main reason for this is that we have much more data in W sample than the Z sample."]}, {"cell_type": "markdown", "id": "c05dcf99", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.5.1</span>\n", "\n", "Define a fit function `fitW()` that combines a Gaussian of the form $a\\exp{-(x-\\mu)^2/(2\\sigma^2)}$ with a 5th order polynomial. The parameters `a`, `mu`, and `sigma`, the parameters of the polynomial are left as fit parameters. "]}, {"cell_type": "code", "execution_count": null, "id": "6e3d38b3", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.5.1\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "def fitW(x, p0, p1, p2, p3, p4, p5, a, mu, sigma):\n", "    #Our model is a gaussian on top of 5th order polynomial.\n", "    \n", "    #Define the polynomial\n", "    poly  = #YOUR CODE HERE\n", "    \n", "    #Define the gaussian\n", "    gauss = #YOUR CODE HERE\n", "    \n", "    #Stick them together\n", "    y =  poly + a*gauss\n", "    \n", "    return y"]}, {"cell_type": "markdown", "id": "8f8cbdb7", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.5.2</span>\n", "\n", "Edit the fit model `p` to set the initial conditions for the fit parameters. Run the fit and see how it looks!\n", "\n", "What is the reduced chi-squared value? Is it good? Report your answer as a number with precision 1e-1."]}, {"cell_type": "code", "execution_count": null, "id": "9f72b700", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.5.2\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "# Now we get the data histogram so we can fit it\n", "bins = 50\n", "mrange=[40,140]\n", "counts, bins = np.histogram(dataW,bins=bins,range=mrange,density=False)\n", "\n", "w = (1/ #poisson unc) #Poisson uncertainty here\n", "binCenters = (bins[1:]+bins[:-1])*.5\n", "x,y = binCenters.astype(\"float32\"), counts.astype(\"float32\")\n", "\n", "#Perform the fit \n", "model = lm.Model(fitW)\n", "     \n", "#Set initial conditions for your fit\n", "#You could experiment with zeros or your intuition first.\n", "#For better fit I suggest adding restrictions to the fit.\n", "p = model.make_params(#YOUR CODE HERE) \n", "\n", "\n", "result_W = model.fit(data=y,\n", "                   params=p,\n", "                   x=x,\n", "                   weights=w)\n", "\n", "#Plot the result\n", "plt.figure()\n", "result_W.plot()\n", "plt.xlabel(\"mass[GeV]\",position=(0.92,0.1))\n", "plt.ylabel(\"Entries/bin\",position=(0.1,0.84))\n", "\n", "#Print the fit summary\n", "print(result_W.fit_report())\n", "result_W.chisqr"]}, {"cell_type": "markdown", "id": "2e5a183d", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Checkpoint 2.5.3</span>\n", "\n", "Finally, get the mass and standard error from your fit results (corresponding to the specific choice of fit function we defined in `Checkpoint 2.5.2`). This should be the parameter `mu`, corresponding to the Gaussian fit. Report the mass with precision 1e-1."]}, {"cell_type": "code", "execution_count": null, "id": "680b6bf9", "metadata": {"tags": ["py", "learner_chopped", "learner", "draft"]}, "outputs": [], "source": ["#>>>PROBLEM: PROJ2.5.3\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "mW = #YOUR CODE HERE\n", "mWerr = #YOUR CODE HERE\n", "\n", "print(mW, \"+/-\", mWerr)"]}], "metadata": {"celltoolbar": "Tags", "kernelspec": {"display_name": "Python 3 (ipykernel)", "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.9"}}, "nbformat": 4, "nbformat_minor": 5}