Warning, /python-analysis-bootcamp/Accessing Truth Information.ipynb is written in an unsupported language. File is not indexed.
0001 {"cells":[{"cell_type":"markdown","metadata":{},"source":["# Accessing Truth Information"]},{"attachments":{},"cell_type":"markdown","metadata":{},"source":["An important aspect at this stage in the design of EIC experiments is the comparison of the so-called 'truth' information with simulated variables. In this notebook we will go over how to access the truth information in an ePIC event."]},{"cell_type":"markdown","metadata":{},"source":["## Importing uproot"]},{"cell_type":"markdown","metadata":{},"source":["Depending on the versions of uproot and XRootD that you have installed, you may encouter a warning from uproot below. Nevertheless, because of the simple data format of the ATHENA ROOT files, we are able to ignore this warning."]},{"cell_type":"code","execution_count":null,"id":"crazy-lambda","metadata":{"trusted":true},"outputs":[],"source":["import uproot as ur\n","print('Uproot version: ' + ur.__version__)"]},{"cell_type":"markdown","metadata":{},"source":["## Opening a file with uproot"]},{"cell_type":"markdown","metadata":{},"source":["To test uproot, we will open a sample file (a single-particle simulation of interest to those who wish to study detector performance):"]},{"cell_type":"code","execution_count":null,"id":"d940cf1b","metadata":{},"outputs":[],"source":["server = 'root://dtn-eic.jlab.org//work/eic2/'\n","dir = 'EPIC/RECO/24.02.0/epic_craterlake/DIS/NC/18x275/minQ2=10/'\n","file = 'pythia8NCDIS_18x275_minQ2=10_beamEffects_xAngle=-0.025_hiDiv_1.0000.eicrecon.tree.edm4eic.root'"]},{"cell_type":"code","execution_count":null,"id":"wicked-amsterdam","metadata":{"trusted":true},"outputs":[],"source":["events = ur.open(server + dir + file + ':events')"]},{"attachments":{},"cell_type":"markdown","metadata":{},"source":["## Truth information in the `MCParticles` branch"]},{"attachments":{},"cell_type":"markdown","metadata":{},"source":["The truth information is stored in the `MCParticles` branch. This includes the true *generated particles* in the simulation, subject to certain conditions. For the purpose of end-user analysis of centrally produced simulation files, the conditions are essentially that only primary particles are included."]},{"cell_type":"markdown","metadata":{},"source":["Several fields are available for the truth information:"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["events['MCParticles'].keys()"]},{"attachments":{},"cell_type":"markdown","metadata":{},"source":["Besides the particle data group code `PDG`, the indices of the parent and daughter tracks, and the generator status `generatorStatus`, you will also see the creation vertex position (`vertex`) and momentum (`momentum`), as well as the endpoint position (`endpoint`) and momentum (`momentumAtEndpoint`) of the particle. Thus, `momentum.x` corresponds to the `x` component of the momentum at the creation vertex. Let's retrieve these momenta, as well as the `PDG` code and `generatorStatus`."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["PDG = events['MCParticles.PDG'].array()\n","generatorStatus = events['MCParticles.generatorStatus'].array()\n","psx,psy,psz = events['MCParticles.momentum.x'].array(), events['MCParticles.momentum.y'].array(), events['MCParticles.momentum.z'].array()"]},{"attachments":{},"cell_type":"markdown","metadata":{},"source":["As expected, for this file the `PDG` code corresponds to many different particles, e.g. for one event we find:"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["PDG[100]"]},{"attachments":{},"cell_type":"markdown","id":"05556486","metadata":{},"source":["Let's compare this with the generator status, and keep in mind the following values in the HepMC3 standard:\n","- 0: undefined and should not occur for well-formed input,\n","- 1: macroscopically stable particles that are thrown into the detector simulation,\n","- 2: unstable particles which have been decayed by the event generator,\n","- 4: incoming beam particles,\n","- > 10: reserved for event generator use, e.g. virtual photons in hard scattering."]},{"cell_type":"code","execution_count":null,"id":"95976f22","metadata":{},"outputs":[],"source":["generatorStatus[100]"]},{"attachments":{},"cell_type":"markdown","metadata":{},"source":["We can also look at the total momentum. For this, we are importing the `numpy` library to use the `sqrt` function on arrays."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["import numpy as np\n","p = np.sqrt(psx**2 + psy**2 + psz**2)"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["p[100]"]},{"cell_type":"markdown","metadata":{},"source":["## Making a simple plot"]},{"cell_type":"markdown","metadata":{},"source":["We can now create a simple plot of the angular (theta) distribution of the generated particles."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["import awkward as ak\n","theta = np.arctan2(np.sqrt(psx**2 + psy**2), psz)\n","theta[generatorStatus == 1]"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["import matplotlib.pyplot as plt\n","plt.hist(-ak.flatten(theta[generatorStatus == 1]), range = (-np.pi, 0), bins = 50)\n","plt.xlabel('Initial Scattering Angle, $-\\\\theta$ [rad]')\n","plt.ylabel('Number of events')\n","plt.show()"]},{"attachments":{},"cell_type":"markdown","id":"bf649e57","metadata":{},"source":["# Exercise"]},{"attachments":{},"cell_type":"markdown","id":"371352c4","metadata":{},"source":["Create a stacked histogram that separates the theta distributions by particle type. Think (or determine) which particle types will be most relevant here."]},{"cell_type":"code","execution_count":null,"id":"96b03433","metadata":{},"outputs":[],"source":["# your code here\n","\n"]},{"attachments":{},"cell_type":"markdown","id":"8b44fe7f","metadata":{},"source":["### Possible solution"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["stacked_theta = [\n"," -ak.flatten(theta[np.logical_and(PDG == 11, generatorStatus == 1)]),\n"," -ak.flatten(theta[np.logical_and(PDG == 22, generatorStatus == 1)]),\n"," -ak.flatten(theta[np.logical_and(np.logical_or(PDG == 211, PDG == -211), generatorStatus == 1)]),\n","]\n","plt.hist(stacked_theta, range = (-np.pi, 0), bins = 50, stacked = True, label = [\"e\", \"$\\\\gamma$\", \"$\\\\pi^\\\\pm$\"])\n","plt.xlabel('Initial Scattering Angle, $-\\\\theta$ [rad]')\n","plt.ylabel('Number of events')\n","plt.legend()\n","plt.show()"]},{"cell_type":"code","execution_count":null,"id":"0dae8be1","metadata":{},"outputs":[],"source":[]}],"metadata":{"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.8"}},"nbformat":4,"nbformat_minor":5}
