Warning, /python-analysis-bootcamp/Determining x and Q^2 from the scattered electron.ipynb is written in an unsupported language. File is not indexed.
0001 {"cells":[{"cell_type":"markdown","id":"subjective-bleeding","metadata":{},"source":["# Determining x and Q<sup>2</sup> from the scattered electron"]},{"attachments":{},"cell_type":"markdown","id":"dense-tsunami","metadata":{},"source":["In this notebook we determine DIS kinematics parameters from the scattered lepton side. For more details, see the page on the [BNL wiki](https://wiki.bnl.gov/eic/index.php/DIS_Kinematics). The following image gives the summary:\n","\n","<img src=\"https://wiki.bnl.gov/eic/upload/Dis.variables.png\" width=\"50%\"/>"]},{"cell_type":"markdown","id":"described-kernel","metadata":{},"source":["## Importing packages"]},{"attachments":{},"cell_type":"markdown","id":"indoor-desperate","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 ePIC ROOT files, we are able to ignore this warning."]},{"cell_type":"code","execution_count":7,"id":"crazy-lambda","metadata":{"trusted":true},"outputs":[],"source":["import numpy as np\n","import uproot as ur\n","import awkward as ak\n","import matplotlib.pyplot as plt"]},{"cell_type":"markdown","id":"differential-pierce","metadata":{},"source":["## Opening a file with uproot"]},{"cell_type":"markdown","id":"thick-mozambique","metadata":{},"source":["To test uproot, we will open a sample file (a DIS simulation sample):"]},{"cell_type":"code","execution_count":17,"id":"0b1c8500","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":18,"id":"wicked-amsterdam","metadata":{"trusted":true},"outputs":[],"source":["events = ur.open(server + dir + file + ':events')"]},{"cell_type":"markdown","id":"bright-ferry","metadata":{},"source":["## Accessing the reconstructed particle momentum"]},{"attachments":{},"cell_type":"markdown","id":"cognitive-standing","metadata":{},"source":["We will access the particle momentum in the `ReconstructedChargedParticles` branch. This contains the reconstructed momentum from the tracking system. As a reminder, here are the available fields:"]},{"cell_type":"code","execution_count":null,"id":"double-gardening","metadata":{"trusted":true},"outputs":[],"source":["events['ReconstructedChargedParticles'].keys()"]},{"attachments":{},"cell_type":"markdown","id":"living-architecture","metadata":{},"source":["For this analysis we will only use the three-momentum `p` and the particle identication code `PDG`. We will select only electrons (`PDG == 11`) and combine them with their initial momentum $\\vec{p}_0$ which, in the ePIC coordinate system, is in the negative $z$ direction by definition."]},{"cell_type":"code","execution_count":null,"id":"respected-ivory","metadata":{"trusted":true},"outputs":[],"source":["reconstructed_charged_particles = events['ReconstructedChargedParticles'].arrays()"]},{"cell_type":"code","execution_count":null,"id":"paperback-arrow","metadata":{"trusted":true},"outputs":[],"source":["kp1, kp2, kp3 = reconstructed_charged_particles['ReconstructedChargedParticles.momentum.x'], reconstructed_charged_particles['ReconstructedChargedParticles.momentum.y'], reconstructed_charged_particles['ReconstructedChargedParticles.momentum.z']\n","PDG = reconstructed_charged_particles['ReconstructedChargedParticles.PDG']\n","m = reconstructed_charged_particles['ReconstructedChargedParticles.mass']"]},{"attachments":{},"cell_type":"markdown","id":"spanish-wonder","metadata":{},"source":["The mass of the particle is stored, but completely defined by the `PDG` code. For the electron it is indeed always 0.511 MeV."]},{"cell_type":"code","execution_count":null,"id":"sporting-ending","metadata":{"trusted":true},"outputs":[],"source":["m[PDG==11]"]},{"cell_type":"markdown","metadata":{},"source":["## Determining the momentum transfer $Q^2$"]},{"cell_type":"markdown","id":"small-queensland","metadata":{},"source":["For all particles we can calculate the energy, which we will consider the zeroth component of the four-momentum $p$."]},{"cell_type":"code","execution_count":null,"id":"collect-secondary","metadata":{"trusted":true},"outputs":[],"source":["kp0 = np.sqrt(m**2+(kp1**2+kp2**2+kp3**2))"]},{"cell_type":"markdown","metadata":{},"source":["The four-momentum of the incoming electron beam has only a $p_z$ and $E$ component."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["k3 = -18\n","m0 = 0.000511\n","k0 = np.sqrt(m0**2 + k3**2)"]},{"cell_type":"markdown","metadata":{},"source":["We can now calculate the components of the four-momentum transfer $q_\\mu = (k_\\mu - k'_\\mu)$:"]},{"cell_type":"code","execution_count":null,"id":"boring-testing","metadata":{"trusted":true},"outputs":[],"source":["q0 = k0 - kp0\n","q1 = - kp1\n","q2 = - kp2\n","q3 = k3 - kp3"]},{"cell_type":"markdown","metadata":{},"source":["With the four components we can form the squared four-momentum transfer, a scalar quantity, which is $Q^2 = -q^2$:"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["Q2 = -(q0**2 - q1**2 - q2**2 - q3**2)"]},{"cell_type":"code","execution_count":null,"id":"boolean-alabama","metadata":{"trusted":true},"outputs":[],"source":["plt.hist(ak.flatten(Q2[PDG==11]), bins=100)\n","plt.yscale('log')\n","plt.xlabel('$Q^2$ [GeV]')\n","plt.ylabel('Number of events')\n","plt.show()"]},{"attachments":{},"cell_type":"markdown","metadata":{"trusted":false},"source":["As we can see, the $Q^2$ falls rapidly in this sample: large $Q^2$ events are less well represented than small $Q^2$ samples. Let's verify that the lower limit is indeed 10 GeV<sup>2</sup>."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["plt.hist(ak.flatten(Q2[PDG==11]), range=[0,100], bins=100)\n","plt.yscale('log')\n","plt.xlabel('$Q^2$ [GeV]')\n","plt.ylabel('Number of events')\n","plt.show()"]},{"attachments":{},"cell_type":"markdown","metadata":{},"source":["We can see that some entries fall below 10 GeV<sup>2</sup>. This is because we have not made any attempt at ensuring that only 1 electron track per event is taken into account. An extension of this exercise could involve selecting only the 'right' electron in those cases where there are multiple electrons in an event. How would you define the 'right' electron? That is indeed not trivial!"]},{"attachments":{},"cell_type":"markdown","id":"a534f9db","metadata":{},"source":["Question: Are there other reasons why we might get values of $Q^2$ below the 10 GeV$^2$ cut off in the event generator?"]},{"cell_type":"markdown","metadata":{},"source":["## Determining the momentum fraction $x$"]},{"attachments":{},"cell_type":"markdown","metadata":{},"source":["In order to determine $x$ we also need the incoming proton momentum $\\vec{p}$. While it might be appealing to think that the proton momentum must be exactly along the $z$ axis as well, this is not the case in the interaction points of the EIC. At interaction point 6 (IP6), the crossing angle is -25 mrad in the $xz$ plane. Thus, the proton four-momentum is:"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["alpha = -0.025\n","p1 = 275 * np.sin(alpha)\n","p2 = 0\n","p3 = 275 * np.cos(alpha)\n","p0 = np.sqrt(0.938**2 + p1**2 + p2**2 + p3**2)"]},{"cell_type":"markdown","metadata":{},"source":["With this proton four-momentum we can now calculate the product $p \\cdot q$, another scalar quantity:"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["pq = p0 * q0 - p1 * q1 - p2 * q2 - p3 * q3"]},{"cell_type":"markdown","metadata":{},"source":["and finally also $x = \\frac{Q^2}{2 pq}$:"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["x = 0.5 * Q2 / pq"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["plt.hist(ak.flatten(x[PDG==11]), range=[-0.5, 1.5], bins=100)\n","plt.yscale('log')\n","plt.xlabel('$x$')\n","plt.ylabel('Number of events')\n","plt.show()"]},{"cell_type":"markdown","metadata":{},"source":["Without a limit on the allowable values of $x$, between 0 and 1 for electron-proton scattering, we are once again easily fooled by outlier events that are caused by non-primary electrons in the event. We can limit our range and choose a logarithmic $x$-axis:"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["plt.hist(ak.to_numpy(x[Q2==ak.max(Q2[PDG==11],1)]), range=[0,1], bins=100)\n","plt.xscale('log')\n","plt.yscale('log')\n","plt.xlabel('$x$')\n","plt.ylabel('Number of events')\n","plt.show()"]},{"cell_type":"markdown","metadata":{},"source":["We are still not entirely happy witht his figure, so we change to logarithmic binning as well:"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["x_bins = np.logspace(-4, 0, 30) # specify exponents of 10\n","plt.hist(ak.flatten(x[PDG==11]), bins=x_bins)\n","plt.xscale('log')\n","plt.yscale('log')\n","plt.xlabel('$x$')\n","plt.ylabel('Number of events')\n","plt.show()"]},{"cell_type":"markdown","metadata":{},"source":["## Kinematic coverage plot for DIS events"]},{"cell_type":"markdown","metadata":{},"source":["Finally, let's put our work on $Q^2$ and $x$ together to calculate a characteristic plot in DIS: the coverage in the $Q^2$ vs $x$ plane."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["x_bins = np.logspace(-4, 0, 30) # specify exponents of 10\n","Q2_bins = np.logspace(-1, 3, 40) # specify exponents of 10\n","plt.hist2d(ak.to_numpy(ak.flatten(x[PDG==11])), ak.to_numpy(ak.flatten(Q2[PDG==11])), bins = [x_bins, Q2_bins])\n","plt.xscale('log')\n","plt.yscale('log')\n","plt.xlabel('$x$')\n","plt.ylabel('$Q^2$ [Gev]')\n","plt.show()"]},{"cell_type":"markdown","metadata":{},"source":["## Next steps"]},{"attachments":{},"cell_type":"markdown","metadata":{},"source":["There are several next steps you can take from here.\n","- You could compare the kinematic coverage for several different beam energy configurations (this requires more than just a change in the file you load).\n","- You could attempt to load multiple files, to increase the statistical precision of your sample.\n","- You could compare your values for $x$ and $Q^2$ with the values produced by the reconstruction, and stored in the `InclusiveKinematicsElectron` collection."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"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.12.1"}},"nbformat":4,"nbformat_minor":5}
