<jupyter_start><jupyter_text>2: First look at data In this lesson we will look at a toy dataset simulating $J/\psi \rightarrow \mu^+ \mu^-$ events. We will discuss ways of loading the data in python,data formats and plotting with ```matplotlib``` and [mplhep](https://mplhep.readthedocs.io/en/latest/index.html). Two plotting libraries?While ```matplotlib``` is a general Python library, the latter is a [Scikit-HEP](https://scikit-hep.org/) package. Scikit-HEP is an organisation which collects, improves and maintains Python packages that are still general, yet cover (often rather special) use cases in High Energy Physics. ```mplhep``` as an example provides plotting functionality  and styles to mimic typical plots that you see in HEP. However, it's built completely on top of matplotlib and the same can be achieved with pure matplotlib (but more cumbersome). Another example is uproot that we will use to read ROOT files.  Recap: Importing modulesIt's generally seen as good practice to put imports at the top of your file:<jupyter_code>import mplhep
import numpy as np
import uproot
from matplotlib import pyplot as plt<jupyter_output><empty_output><jupyter_text>5. The toy dataset We're going to look at some fake $J/\psi \rightarrow \mu^+ \mu^-$ data, the variables available are:- `Jpsi_M` `Jpsi_P` `Jpsi_PT` `Jpsi_PE` `Jpsi_PX` `Jpsi_PY` `Jpsi_PZ`- `mum_M` `mum_PT` `mum_eta` `mum_PE` `mum_PX` `mum_PY` `mum_PZ` `mum_IP` `mum_ProbNNmu` `mum_ProbNNpi`- `mup_M` `mup_PT` `mup_eta` `mup_PE` `mup_PX` `mup_PY` `mup_PZ` `mup_IP` `mup_ProbNNmu` `mup_ProbNNpi`- `nTracks`The meaning of the suffixes are as follows:- `_M`: Invariant mass of the particle (fixed to the PDG value for muons)- `_P`: Absolute value of the particle's three momentum- `_PT`: Absolute value of the particle's momentum in the `x`-`y` plane- `_PE`, `_PX`, `_PY`, `_PZ`: Four momentum component of the particle- `_IP`: Impact parameter, i.e. the distance of closest approach between the reconstructed particle and the primary vertex- `ProbNNmu`, `ProbNNpi`: Particle identification variables which correspond to how likely it is that the particle is a muon or a pion- `nTracks`: The total number of tracks in the event  Loading data - `uproot` is a way of reading ROOT files without having ROOT installed, see the github reposity [here](https://github.com/scikit-hep/uproot) - We can look at the objects that are available in the file and access objects using dictionary style syntax - The tree class contains converters to a variety of common Python libraries, such as numpy - We will also use `pandas DataFrames` to load data in a table like format (- `root_numpy` and `root_pandas` are an alternative, yet outdated, way of reading+writing ROOT files)First let's load the data using `uproot`.Often it is convenient to access data stored on the *grid* at CERN, so you don't have to keep it locally. This can be done using the *XRootD* protocol:```pythonmy_file = uproot.open('root://eosuser.cern.ch//eos/user/l/lhcbsk/advanced-python/data/real_data.root')```or even better```pythonwith uproot.open('root://eosuser.cern.ch//eos/user/l/lhcbsk/advanced-python/data/real_data.root') as my_file:   ...   do something here with my_file```Accessing data this way requires you to have valid CERN credentials to access it. If authentication fails you will see an error message like:```OSError: [ERROR] Server responded with an error: [3010] Unable to give access - user access restricted - unauthorized identity used ; Permission denied```Credentials can be obtained by typing `kinit <EMAIL>` in your terminal and entering your CERN password.For this tutorial we will use a publicly accessible file instead, using HTTPS to access it remotely. This is significantly slower than using XRootD.<jupyter_code>my_file = uproot.open('https://cern.ch/starterkit/data/advanced-python-2018/real_data.root',
                      httpsource={'chunkbytes': 1024*1024, 'limitbytes': 33554432, 'parallel': 64})

my_file.keys()<jupyter_output><empty_output><jupyter_text>In order to get a single variable, we can use the indexing. The `array` convert<jupyter_code>tree = my_file['DecayTree']
# Get a numpy array containing the J/Ψ mass
tree['Jpsi_M'].array(library='np')
# Load data as a pandas DataFrame
data_df = tree.arrays(library='pd')
my_file.close()  # usually, it's better to open the file with a "with" statement -> needs no closing

# Show the first 5 lines of the DataFrame
data_df.head()
data_df.columns<jupyter_output><empty_output><jupyter_text>6. Plotting a simple histogram<jupyter_code># Start with a basic histogram
plt.hist(data_df['Jpsi_M'])
plt.xlabel('Jpsi mass')<jupyter_output><empty_output><jupyter_text>That's okay but we could improve it.Matplotlib - https://matplotlib.org/api/_as_gen/matplotlib.pyplot.hist.html - Internally it calls `np.histogram` and  returns an array of counts, an array of bins (as `np.histogram` does) and an array of patches. - We can use `histtype="step"` so we can plot multiple datasets on the same axis easilyHowever, this leaves us with a few things: - what if we have our data already as a histogram? - We usually want to have uncertainties on the histogram - It should match rather the usual style of plotsThat's where ```mplhep``` comes into play!<jupyter_code># One of the features of mplhep is that we can set an experiment specific style

plt.style.use(mplhep.style.LHCb2)  # or ATLAS/CMS
# plotting again
plt.hist(data_df['Jpsi_M'], bins=40)
plt.xlabel('Jpsi mass')
# And similar with mplhep
mplhep.histplot(*np.histogram(data_df['Jpsi_M'], bins=30))
plt.xlabel('Jpsi mass')
# ... except that it can do a lot more!

# we need to histogram only once!
h, bins = np.histogram(data_df['Jpsi_M'], bins=50)  # TRY OUT: change the binning
plt.subplots(1,2,figsize=(20, 6))
plt.subplot(1, 2, 1)
mplhep.histplot(h, bins, yerr=True)  # error can also be array
plt.subplot(1, 2, 2)
half_binwidths = (bins[1] - bins[0]) / 2
mplhep.histplot(h, bins, histtype='errorbar', yerr=True, xerr=half_binwidths)
def plot_mass(df):
    h, bins = np.histogram(df['Jpsi_M'], bins=100, range=[2.75, 3.5])
    mplhep.histplot(h, bins, yerr=True)  # feel free to adjust
    # You can also use LaTeX in the axis label
    plt.xlabel('$J/\\psi$ mass [GeV]')
    plt.xlim(bins[0], bins[-1])

plot_mass(data_df)<jupyter_output><empty_output><jupyter_text>Adding variables<jupyter_code># When making the ROOT file we forgot to add some variables, no bother lets add them now!
data_df.eval('Jpsi_eta = arctanh(Jpsi_PZ/Jpsi_P)', inplace=True)
data_df.head()['Jpsi_eta']<jupyter_output><empty_output><jupyter_text>**Exercise:** Add `mu_P` and `mum_P` columns to the DataFrame.<jupyter_code>data_df.eval('mup_P = sqrt(mup_PX**2 + mup_PY**2 + mup_PZ**2)', inplace=True)
data_df.eval('mum_P = sqrt(mum_PX**2 + mum_PY**2 + mum_PZ**2)', inplace=True)
# We can also get multiple columns at the same time
data_df.head()[['mum_P', 'mup_P']]<jupyter_output><empty_output><jupyter_text>Using rectangular cuts * We want to increase the 'signal significance' of our sample - this means more signal events with respect to background* To do this we can cut on certain discriminating variables* Here we will make cuts on the `Jpsi_PT` and **PID** (Particle Identification) variables<jupyter_code>plot_mass(data_df)
data_with_cuts_df = data_df.query('Jpsi_PT > 4')
plot_mass(data_with_cuts_df)
plot_mass(data_df)
data_with_cuts_df = data_df.query('Jpsi_PT > 4')
plot_mass(data_with_cuts_df)
# Lets add some PID cuts as well
data_with_cuts_df = data_df.query('(Jpsi_PT > 4) & ((mum_ProbNNmu > 0.9) & (mup_ProbNNmu > 0.9))')
plot_mass(data_with_cuts_df)<jupyter_output><empty_output><jupyter_text>Let's go back and add a label argument to our plot function. This makes it easier to identify each line.We can also use the `density` argument in `matplotlib.hist` to plot all the histograms as the same scale.<jupyter_code>def plot_mass(df, **kwargs):
    h, bins = np.histogram(df['Jpsi_M'], bins=100, range=[2.75, 3.5])
    mplhep.histplot(h, bins, yerr=True, **kwargs)  # feel free to adjust
    # You can also use LaTeX in the axis label
    plt.xlabel('$J/\\psi$ mass [GeV]')
    plt.xlim(bins[0], bins[-1])

plot_mass(data_df, label='No cuts', density=1)
data_with_cuts_df = data_df.query('Jpsi_PT > 4')
plot_mass(data_with_cuts_df, label='$J/\\psi$ p$_T$ only', density=1)
data_with_cuts_df = data_df.query('(Jpsi_PT > 4) & ((mum_ProbNNmu > 0.9) & (mup_ProbNNmu > 0.9))')
plot_mass(data_with_cuts_df, label='$J/\\psi$ p$_T$ and muon PID', density=1)
plt.legend(loc='best')<jupyter_output><empty_output><jupyter_text>For this tutorial we have a special function for testing the significance of the signal in our dataset. There are many ways to do this with real data, though we will not cover them here.<jupyter_code>data_df.columns
from python_lesson import check_truth

print('Originally the significance is')
check_truth(data_df)

print('\nCutting on pT gives us')
check_truth(data_df.query('Jpsi_PT > 4'))

print('\nCutting on pT and ProbNNmu gives us')
check_truth(data_df.query('(Jpsi_PT > 4) & ((mum_ProbNNmu > 0.9) & (mup_ProbNNmu > 0.9))'))<jupyter_output><empty_output><jupyter_text>Comparing distributionsBefore we just used the cuts that were told to you but how do we pick them?One way is to get a sample of simulated data, we have a file in `data/simulated_data.root`.**Exercise:** Load it into a pandas `DataFrame` called `mc_df`. Don't forget to add the `Jpsi_eta`, `mup_P` and `mum_P` columns!<jupyter_code>with uproot.open('https://starterkit.web.cern.ch/starterkit/data/advanced-python-2018/simulated_data.root') as mc_file:
    mc_df = mc_file['DecayTree'].arrays(library='pd')

# mc_file = uproot.open('https://starterkit.web.cern.ch/starterkit/data/advanced-python-2018/simulated_data.root')
# mc_df = mc_file['DecayTree'].arrays(library='pd')
mc_df.eval('Jpsi_eta = arctanh(Jpsi_PZ/Jpsi_P)', inplace=True)
mc_df.eval('mup_P = sqrt(mum_PX**2 + mum_PY**2 + mum_PZ**2)', inplace=True)
mc_df.eval('mum_P = sqrt(mum_PX**2 + mum_PY**2 + mum_PZ**2)', inplace=True)<jupyter_output><empty_output><jupyter_text>**QUESTION:** What can we to get a background sample?Sidebands, we know the peak is only present in $3.0~\text{GeV} < M(J/\psi) < 3.2~\text{GeV}$. If we select events outside the region we know it's a pure background sample.**Exercise:** Make a new `DataFrame` called `bkg_df` containing only events outside $3.0~\text{GeV} < M(J/\psi) < 3.2~\text{GeV}$.<jupyter_code>bkg_df = data_df.query('~(3.0 < Jpsi_M < 3.2)')
plot_mass(bkg_df)<jupyter_output><empty_output><jupyter_text>**QUESTION:** Why is there a step at 3 GeV on this plot?It's a binning effect, we've appled a cut at 3.0 GeV but the nearest bin is $[2.9975, 3.005]$ so it is only partially filled. Now let's plot the variables in MC and background to see what they look like?<jupyter_code>var = 'Jpsi_PT'
# let's first create the histograms
hsig, bins = np.histogram(mc_df[var], bins=60)
hbkg, bins = np.histogram(bkg_df[var], bins=bins)  # using the same bins here
# then plot them
mplhep.histplot((hsig, bins), label='MC Signal')
mplhep.histplot(hbkg, bins=bins, label='Data Bkg')
plt.xlabel(var)
plt.xlim(bins[0], bins[-1])
plt.legend(loc='best')
# Those are hard to compare!
# We can add the density keyword argument to normalise the distributions

mplhep.histplot(hsig, bins=bins, label='MC Signal', density=1)
mplhep.histplot(hbkg, bins=bins, label='Data Bkg', density=1)
plt.xlabel(var)
plt.xlim(bins[0], bins[-1])
plt.legend(loc='best')<jupyter_output><empty_output><jupyter_text>**Exercise:** Make a function which plots both variables with the signature `plot_comparision(var, mc_df, bkg_df)`.<jupyter_code>def plot_comparision(var, mc_df, bkg_df):
    # create histograms
    hsig, bins = np.histogram(mc_df[var], bins=60, density=1)
    hbkg, bins = np.histogram(bkg_df[var], bins=bins, density=1)

    mplhep.histplot((hsig, bins), label='MC Signal', )
    mplhep.histplot(hbkg, bins=bins, label='Data Bkg')
    plt.xlabel(var)
    plt.xlim(bins[0], bins[-1])
    plt.legend(loc='best')<jupyter_output><empty_output><jupyter_text>We can now use this function to plot all of the variables available in the data using `data_df.columns`:<jupyter_code>for var in data_df.columns:
    plt.figure()  # creates a new figure
    plot_comparision(var, mc_df, bkg_df)<jupyter_output><empty_output><jupyter_text>Things to note: - This doesn't work for the $J/\psi$ mass variable: We can only rely on this method if the variable is independent of the mass. Fortunately we often want to do this as if a variable is heavily dependent on mass it can shape our distributions and can make it hard to know what is signal and what is background. - Muon mass is a fixed value for all muons: In this sample we have assumed the PDG value of the muon mass to allow us to calculate the energy component using only the information from the tracking detectors. This is often more precise than using calorimeters to measure $P_E$. - We got a warning about `More than 20 figures have been opened.`: Opening plots uses memory so if you open too many at the same time your scripts can be become slow or even crash. In this case we can ignore it as we only produce 30 plots but be careful if you ever make thousands of plots. - Pseudorapidity (eta) only goes between about 1 and 6: This dataset is supposed to look vaguely like LHCb data where the detector only covers that region. **Exercise:** Look at the variables above and try to get a clean $J/\psi$ mass peak. Use the significance function to see how well you do. **Aside:** We will want to use some of the variables and dataframes in the next lesson. In order to do this we will *store* them in this session and reload them in the next lesson.<jupyter_code>%store bkg_df
%store mc_df
%store data_df<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11287171695246365, 'fraction_numerical': 0.011359525359959481, 'mean_word_length': 5.322964318389753, 'pattern_counts': {'":': 0, '<': 72, '<?xml version=': 0, '>': 80, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '639946', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4927, 'n_tokens_neox': 4570, 'n_words': 1614, 'path': 'advanced-python/20DataAndPlotting.ipynb', 'repo_name': 'vlukashenko/analysis-essentials', 'size': '13830'}