Week6/ATLAStestbeam/ATLAStestbeam.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ATLAS Test Beam Data\n",
    "\n",
    "Python script for analysis of ATLAS test beam data. \n",
    "\n",
    "The program reads an ASCII (i.e. text) data file, containting a large number of events, where a charged particle (electron or pion) passed through a slice of the ATLAS detector. Each passage is recorded by different detectors (essentially three independent ones!), boiling down to eleven numbers (some more relevant than others). The exercise is to separate electron and pion events based on these numbers, and in turn use this information to measure the interaction (i.e. detector performance) of pions and electrons seperately.\n",
    "\n",
    "NOTE 1: In this exercise, don't consider (at least at first) the details, but just rely on the advice from an old particle physicist, that in these detectors, electrons tend to give a higher signal (i.e. values) than pions do, and that none of the particle types are very rare.\n",
    "\n",
    "NOTE 2: Though the data is from particle physics, it could in principle have been from ANY other source, and the eleven numbers could for example have been indicators of cancer, key numbers for investors, or index numbers for identifying potential costumors. Thus the advice from above could be from and old doctor, economist, or marketing consultant.\n",
    "\n",
    "NOTE 3: This is real data, and while it has been prepared somewhat, it might contain surprises and does not necessarily follow known PDFs.\n",
    "\n",
    "For more information on ATLAS test beam: http://www.nbi.dk/~petersen/Teaching/Stat2023/TestBeam/TestBeamDataAnalysis.html.\n",
    "\n",
    "***\n",
    "\n",
    "### Authors: \n",
    "- Troels C. Petersen (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "- 29-12-2023 (latest update)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set parameters of the program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write extensive output\n",
    "verbose = True\n",
    "N_verbose = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open file, `DataSet_AtlasPid_ElectronPion_2GeV.txt`, and read in all the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt('DataSet_AtlasPid_ElectronPion_2GeV.txt', skiprows=1, unpack=True)\n",
    "Cher, nLT, nHT, EM0, EM1, EM2, EM3, Had0, Had1, Had2, Muon = data\n",
    "Cher.shape\n",
    "\n",
    "# It is always very \"healthy\" to inspect the data visually/manually, just to get a quick feel for it.\n",
    "if (verbose) :\n",
    "    for i in range (N_verbose) :\n",
    "        print(f\"Cher: {Cher[i]:6.1f}    nLT, nHT: {int(nLT[i]):2d}, {int(nHT[i]):2d}    EM: {EM0[i]:5.2f} {EM1[i]:5.2f} {EM2[i]:5.2f} {EM3[i]:5.2f}   Had: {Had0[i]:5.2f} {Had1[i]:5.2f} {Had2[i]:5.2f}    Muon: {Muon[i]:5.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your analysis:\n",
    "\n",
    "Through plotting and making selections in one/several variable(s) seeing the impact (e.g. how many events remain) in others, work out what this sample consists of. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "# This is where your analysis should go:\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example plot could be to see if there are any clear types difference between the entries in the data, and if this gives rise to a way of selecting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(16, 10))\n",
    "\n",
    "ax[0,0].hist(Cher, bins=200, range=(400, 1400), histtype='step', label='Cherenkov')\n",
    "ax[0,0].set(xlabel=\"Cherenkov\", ylabel=\"Frequency\")\n",
    "\n",
    "ax[1,1].hist(EM1,  bins=200, range=(-0.5, 2.0), histtype='step', label='EM1')\n",
    "ax[1,1].set(xlabel=\"EM1\", ylabel=\"Frequency\")\n",
    "\n",
    "h = ax[0,1].hist2d(EM1,  Cher, bins=(100, 100), range=((-0.5, 2.0), (400, 1400)), norm=mpl.colors.LogNorm(), cmap=\"Reds\")\n",
    "plt.colorbar(h[3], ax=ax[0, 1])                    # z-scale on the right of the figure\n",
    "ax[0,1].set(xlabel=\"EM1\", ylabel=\"Cherenkov\")\n",
    "\n",
    "# You can add your own fourth plot:\n",
    "fig.delaxes(ax[1,0])\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: The simple plotting function gives a warning for EM1 (not Cher). Ask yourself, if you did a check of the original data file input values in any way? That is always \"healthy\"!\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions to be answered:\n",
    "\n",
    "Generally, this analysis is about separating electrons and pions (and determining how well this can be done), followed by a few questions characterizing the detector response to each type of particle. Thus, you should imagine, that your new detector/equipment/questionaire gave you this output, and now it is up to you to find out, what this tells you about your experiment, and how to extract information from it in the best possible way. Typically, it will have taken you months (if not years) to get to this point.\n",
    "\n",
    "Note that this data is NOT meant for “fine tuned analysis”, but rather “crude inspection”. You should try to get simple approximate answers out - I’m sure that you will afterwards be able to fine tune them.\n",
    "\n",
    "Below are questions guiding you, some/most of which your analysis should cover, but you do **not** have to follow them blindly (I've put \"Optional\" on those that are not essential). Start by considering the data, and get a feel for the typical range of each variable. Plot the variables, both in 1D and also 2D! From considering these plots, guess/estimate an approximate knowledge of how electrons and pions distribute themselves in the variables above, and how to make a selection of these.\n",
    "\n",
    "As described on the webpage introducing the data, there are three (relevant) detectors:\n",
    "- Cherenkov (1 value),\n",
    "- TRT (Transition Radiation Tracker, 2 values) and\n",
    "- Calorimeters (4 values)\n",
    "\n",
    "They are each capable of separating electrons and pions. As they are (largely) _INDEPENDENT_ (three separate detectors), they may be used to cross check each other, and this is what you should use, in fact the essential part of this (and many other) analysis!\n",
    "\n",
    "\n",
    "Questions:\n",
    "----------\n",
    "1. Find for each of these three detector systems one variable, which seem to separate electrons and pions best. For example, start with the Cherenkov, which is only a single number, and assume/guess that the large peak at low values is mainly from pions, while the upper broad peak is from electrons (this you would know, as you designed the experiment). Now plot the TRT and Calorimeter distributions when the Cherenkov selects a pion and afterwards an electron. This should give you a good idea about how to separate pions and electrons using the TRT and Calorimeters.\n",
    "\n",
    "    Hint: Sometimes variables from a single detector are more powerful, when they are combined, e.g. taken ratios of (or used in a Fisher or ML algorithm). For the TRT this may be somewhat doable, but for the EMcalo, it is not as simple. Here, one variable caries most of the separation power, but involving other layers may enhance the separation power. However, to begin with, just consider a single number from each detector.\n",
    "\n",
    "\n",
    "2. Next you should try to see, if you can make a selection, which gives you a fairly large and clean electron and pion sample, respectively. The question is, how can you know how clean your sample is and how efficient your selection is?  This can actually be measured in the data itself, using the fact that there are three independent detectors. For example, start by making an electron and a pion selection using two of the three variables, and plot the third variable for each of these selections. Now you can directly see, how electrons and pions will distribute themselves in this third variable. Are you worried, that there are pions in your electron sample, and vice versa? Well, there will probably be, but so few, that it won't matter too much, at least not to begin with. Why? Well, let us assume that for each detector, 80% of electrons pass your requirement, but also 10% of pions do. Assuming an even number of electrons and pions (which is not really the case), then with two detector cuts, you should get a sample, which is: $\\frac{0.8\\cdot0.8} {0.8\\cdot0.8 + 0.1\\cdot0.1} = 98.5\\%$ pure.\n",
    "\n",
    "    Now with this sample based on cuts on the two other detectors, ask what fraction of electrons and pions passes your electron selection in the remaining detector. The fraction of electrons, that are not selected as electrons will be your TYPE I errors (False Negative Rate), denoted alpha, while the fraction of pions, that do pass the cut will be your TYPE II errors (False Positive Rate), denoted beta. Measure these for each of the two cuts in the three detector types, and ask yourself if they are \"reasonable\", i.e. something like in the example above. If not, then you should perhaps reconsider adjusting your cuts.\n",
    "\n",
    "By now, you should for each detector have 6 numbers (first considering electrons and then pions as \"Positive\"):\n",
    "    - The electron cut value above which you accept an electron.\n",
    "    - The efficiency (i.e. True Positive Rate) for electrons of this cut.\n",
    "    - The fake rate (i.e. False Positive Rate) for pions of this cut.\n",
    "    - The pion cut value below which you accept a pion (may be same value as above for electrons!).\n",
    "    - The efficiency (i.e. True Positive Rate) for pions of this cut.\n",
    "    - The fake rate (i.e. False Positive Rate) for electrons of this cut.\n",
    "\n",
    "\n",
    "3. Given the efficiencies and fake rates of each cut, try to combine these (again assuming that they are independent) into knowledge of your sample purities and also the total number of electrons and pions in the whole sample. Do the sum of estimated electrons and pions added actually match the number of particles in total? This is a good cross check!\n",
    "\n",
    "### More advanced questions:\n",
    "\n",
    "4. If the number of pions was suddenly 1000 times that of elections, would you still be able to get a sample of fairly pure (say 90% pure) electrons? And if so, what would the efficiency for these electrons be? That is roughly equivalent of asking, if you can get a 99.9% pure electron sample from the data given (where pions do not dominate in numbers).\n",
    "\n",
    "\n",
    "5. Expanding on problem 2), try now to calculate ROC curves for each of the three detectors. These are obtained by making a clean selection using the two other detectors for electrons and pions seperately, and then integrating over these two distributions, using the running (normalised) integral of each as x and y coordinate in a (ROC) curve. If you do not manage on your own, perhaps consider the ROC calculator example, which is posted along with this exercise.\n",
    "\n",
    "\n",
    "6. One of the purposes of the testbeam was to measure the response of the TRT detector to exactly electrons and pions. Consider for example only events that has 33 TRT hits (i.e. `nLT` $= 33$). As the High-Threshold probability (i.e. probability of passing the High-Threshold, given that the Low-Threshold was passed), is assumed to be constant in the TRT detector (but quite different for electrons and pions), what distribution should the number of High-Threshold hits (`nHT`) follow? And is that really the case, both for electrons and pions?\n",
    "\n",
    "\n",
    "7. Still considering `nLT` $=33$, and given that there are both electrons and pions in the sample (each with a different HT probability), `nHT` should in the unselected data be a combination of two distributions. Try to fit the number of HT hits with two distributions combined. Do they fit the data? And can this fit be used to estimate the fraction of pions and electrons in the sample? And does that estimate match you previous estimate? Perhaps retry with other values for the number of TRT hits.\n",
    "\n",
    "\n",
    "8. Try to select pions using three different (mutually exclusive) techniques:\n",
    "    1. Passing only a hadronic calorimeter requirement (e.g. that the sum of the three HCal values is above some minimum energy).\n",
    "    2. Passing only Cherenkov AND EMcalo requirements.\n",
    "    3. Passing both A) and B).\n",
    "\n",
    "Try to measure the HT probability (i.e. fraction of High-Threshold hits) for each of these three pion samples. Do they agree with each other?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning points:\n",
    "\n",
    "This exercise might be the closest thing to reality regarding real data that is not easily modelled. Note that the original raw data is a lot more messy, and that a lot of cleaning has taken place, before you got it!\n",
    "\n",
    "From this exercise you should realise:\n",
    "1. that real data is not as clean-cut as you might think.\n",
    "2. that the use of independent variables is very powerful in data analysis.\n",
    "3. that you need to be sharp regarding TPR, FPR, TNR, and FNR... :-)\n",
    "4. that the sign of a good analysis is, that cross checks have been thought into it (e.g. end of Q.3)"
   ]
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "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.10.13"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}