<jupyter_start><jupyter_text>Study of  $D^0 \to K^- \pi^+ $ decay- Mass- Lifetimehttp://lhcb-public.web.cern.ch/lhcb-public/en/LHCb-outreach/masterclasses/en/D0Lifetime.html<jupyter_code>%%bash
pip install iminuit probfit seaborn cython
import pandas as pd
import matplotlib.pylab as plt
%matplotlib inline
import numpy as np
import pandas as pd
import seaborn
import iminuit, probfit
from __future__ import print_function
%%bash
if [ ! -f masterclass_all.csv.gz ] ; then
  curl -L -o masterclass_all.csv.gz https://www.dropbox.com/s/57i1ulucgeqhjmz/masterclass_all.csv.gz?dl=0
fi
ls -l masterclass_all.csv.gz
#df = pd.DataFrame(root_numpy.root2array(datafile_2012, treename='DecayTree'))
df0 = pd.read_csv("masterclass_all.csv.gz")
print ("Size:" , df0.shape)
df0.head()

plt.hist(df0['D0_MM'], bins=100, histtype='stepfilled', alpha=0.3);
plt.xlabel("Mass, MeV");<jupyter_output><empty_output><jupyter_text>Preselection<jupyter_code>df = df0[df0['D0_MINIP'] < 2][df0['Kplus_ProbNNk']> 0.2][df0['piminus_ProbNNpi']> 0.2]
print ("Size:" , df0.shape)
plt.hist(df['D0_MM'], bins=100, histtype='stepfilled', alpha=0.3);
df.shape<jupyter_output><empty_output><jupyter_text>D0 mass measurements* Fitting the D0 mass    - First plot the D0 mass distribution    - Now fit to it. In this fit, leave the signal and mass ranges to their default values.    - Define the signal mass region as +-3 sigma around the mean value    - What is the signal significance<jupyter_code>from iminuit import Minuit
from probfit import UnbinnedLH, gaussian, linear
fit_range = (1830, 1900)
normalized_poly = probfit.Normalized(probfit.Polynomial(1), fit_range)
normalized_poly = probfit.Extended(normalized_poly, extname='NBkg')

gauss1 = probfit.Extended(probfit.rename(probfit.gaussian, ['x', 'mu1', 'sigma1']), 
                          extname='N1')

# Define an extended PDF consisting of three components
pdf = probfit.AddPdf(normalized_poly, gauss1)

print('normalized_poly: {}'.format(probfit.describe(normalized_poly)))
print('gauss1:          {}'.format(probfit.describe(gauss1)))
print('pdf:             {}'.format(probfit.describe(pdf)))
binned_likelihood = probfit.BinnedLH(pdf, df['D0_MM'], bins=200, extended=True, bound=fit_range)

# This is a quite complex fit (11 free parameters!), so we need good starting values.
# Actually we even need to set an initial parameter error
# for 'mu1' and 'mu2' to make MIGRAD converge.
# The initial parameter error is used as the initial step size in the minimization.
pars = dict(mu1=1865, sigma1=10, N1=35000,
            c_0=0.1, c_1=0.01, NBkg=20000)
minuit = iminuit.Minuit(binned_likelihood, pedantic=False, print_level=0, **pars)
# You can see that the model already roughly matches the data
binned_likelihood.draw(minuit, parts=True);
minuit.migrad();
#minuit.minos();
binned_likelihood.show(minuit, parts=True);
minuit.print_fmin()
minuit.print_matrix()
print ("D0 mass:", minuit.values['mu1'], "+/-", minuit.errors['mu1'])<jupyter_output>D0 mass: 1866.63834187 +/- 0.0680419690226<jupyter_text>D0 lifetime measurementWelcome to the LHCb masterclass exercise on measuring the lifetime of the D0 meson. The goal of this exercise is to measure the lifetime of the D0 meson, a fundamentalparticle made of a charm quark and an up anti-quark. In order to do so, you willfirst learn how to separate signal D0 mesons from backgrounds. Finally, you willcompare your results to the values found by the Particle Data Group (http://pdgLive.lbl.gov). Step-by-step instructions :3. Plot the variable distributions. You will see three further plots appearing, andin each one the blue points represent the distribution of the signal in that variablewhile the red points represent the distribution of the background. The plot is logarithmicin the Y axis, and each point represents the fraction of the total signal in that bin.Which regions of each variable contain mostly signal? Which contain mostly background?4. Fit the lifetime distribution. Save the resultsof your fit and compare them to the PDG value. Do they agree?5. Repeat step 4 but now varying the upper D0 log(IP) variable rangefrom 1.5 to -2 in steps of 0.2. Do you notice a pattern?6. Does the D0 lifetime with an log(IP) cut of-1.5 agree better or worse with the PDG than the lifetime with an log(IP) cut of 1.5?;  LifetimeExponential decay:$$\frac{dN}{dt} = -\lambda N$$where $\lambda$ is decay rate, or $$\frac{dN}{N} = -\lambda dt$$integrating:$$\mathrm{ln} N = - \lambda t + C$$$$ N(t) = e^Ce^{-\lambda t} = N_0 e^{-\lambda t} $$Mean lifetime (https://en.wikipedia.org/wiki/Exponential_decayDerivation_of_the_mean_lifetime):$$\tau = \langle t \rangle = \int_0^\infty t \cdot c \cdot N_0 e^{-\lambda t}\, dt = \int_0^\infty \lambda t e^{-\lambda t}\, dt = \frac{1}{\lambda}$$  Variables distribution<jupyter_code>plt.figure(figsize=(18, 6))
plt.subplot(1,3,1)
plt.hist(df['D0_TAU']*1000, bins=50, histtype='stepfilled', alpha=0.5)
plt.yscale('log', basey=10)
plt.xlabel('t, ps')
plt.ylabel('log10(N)')
plt.title('D0_TAU, Decay time')

plt.subplot(1,3,2)
plt.hist(np.log10(df['D0_MINIP']), bins=50, histtype='stepfilled', alpha=0.5);
plt.yscale('log')
plt.title('log10(D0_MINIP)')
plt.xlabel('log10(MINIP)')

plt.subplot(1,3,3)
plt.hist(df['D0_MM'], bins=50, histtype='stepfilled', alpha=0.5);
plt.xlabel('MeV')
plt.yscale('log')
plt.title('D0_MM')<jupyter_output><empty_output><jupyter_text>What is MINIP?__IP__ - Impact parameter, distance between two tracks, the less the more likely they come from the same vertex;__MINIP__ - minimum IP.  Background substractionidea: 1. We are interested in signal/background sample distribution (not individual events) wrt variable TAU.2. Check that variable MM is not correlated with TAU.3. For the whole sample define signal region R_s wrt mass variable (MM), i.e. interval with majority of signal. It splits whole data sample into 2 subsamples: "signal region events" and "sideband region events".4. Since TAU and MM are not correlated, distribution of a variable TAU for background to be the distribution of _sideband_ events. 5. To plot normalized distribution of signal wrt variable TAU, one have to substract histogram  of background events (4) from histogram of _signal region_ events.  Checking correlation between MM & TAU<jupyter_code>np.cov(df['D0_MM'],df['D0_TAU'])
np.cov(df['D0_MM'],df['D0_PT'])
np.cov(df['D0_MM'],df['D0_MINIP'])
seaborn.jointplot(df['D0_MM'],np.log(df['D0_TAU']*1000), kind='hex').\
  set_axis_labels("D0_MM", "log(D0_tau)");<jupyter_output><empty_output><jupyter_text>Splitting into Signal and Sideband regions<jupyter_code>signal_region = (1844, 1890)
signal_region_mask = df['D0_MM'] >= signal_region[0]
signal_region_mask &= df['D0_MM'] < signal_region[1]
signal_region_events = df[signal_region_mask]
sideband_region_mask = signal_region_mask == False
sideband_region_events = df[sideband_region_mask]
back_scaling_factor = float(signal_region_events.shape[0]) / sideband_region_events.shape[0]
back_scaling_factor
# fig, (a1, a2) = plt.subplots(1, 2, figsize=(18, 6))
a = seaborn.jointplot(signal_region_events['D0_MM'], signal_region_events['D0_TAU']*1000,
                  kind='hex');
a.ax_joint.set_title("aasd")
seaborn.jointplot(sideband_region_events['D0_MM'],sideband_region_events['D0_TAU']*1000, kind='hex', color='r');
n_bins = 30
h_sr, bar_x, _ = plt.hist(signal_region_events['D0_TAU']*1000, bins=n_bins, histtype='stepfilled', alpha=0.5)
h_sb, _, _ = plt.hist(sideband_region_events['D0_TAU']*1000, bins=n_bins, histtype='step', alpha=1., lw=2, color='r')
h_sig = h_sr - h_sb / back_scaling_factor
plt.yscale('log')
plt.xlabel('ps')
plt.title('D0_TAU, Decay time');
plt.bar(bar_x[:-1], h_sig, width=bar_x[1] - bar_x[0])
plt.yscale('log')
plt.xlabel('ps')
plt.title("Distribution Difference");
def line(x, m, c): # define it to be parabolic or whatever you like
    return m * x + c
err = np.maximum(np.nan_to_num(np.log(h_sig + np.sqrt(h_sig)) - np.log(h_sig)), 0.01)
#err = np.ones(len(h_sig))
chi2 = probfit.Chi2Regression(line, (bar_x[:-1] + bar_x[1:])/2, np.log(h_sig), err)
# plt.errorbar((bar_x[:1] + bar_x[1:])/2*1000, np.log(h_sig), yerr=err, fmt='none');
# see iminuit tutorial on how to give initial value/range/error
minuit = iminuit.Minuit(chi2, m=-1, c=10, limit_m=(-5,0), error_m=0.01, error_c=0.1) 
minuit.migrad();<jupyter_output><empty_output><jupyter_text>Result<jupyter_code>chi2.draw(minuit)
print ("D^{0} lifetime %.4f #pm %.6f (ps)" % 
  (-1./minuit.values["m"], minuit.errors["m"]/(minuit.values["m"]**2)))<jupyter_output>D^{0} lifetime 0.7932 #pm 0.003445 (ps)<jupyter_text>Tasks* MINIP dependence    - Repeat the previous steps (starting from Preselection) for different values of MINIP (from 10 to 0.01)    - Plot the trend of lifetime vs. MINIP threshold.<jupyter_code># TODO, make the plot described above<jupyter_output><empty_output><jupyter_text>Quiz* Does the Gaussian fit the mass distribution well?    - Explain the reason why in the fits where there is less background, the Gaussian undershoots the data points on the left.* Why does the original lifetime fit not agree with the PDG value?* Why does cutting on MINIP help, but cutting on D0 PT or TAU not help?* How would you estimate the systematic uncertainty on the measurement from the trend plot?<jupyter_code># Put your answers below<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13214821160848148, 'fraction_numerical': 0.03394731205825659, 'mean_word_length': 5.637526652452026, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 38, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '260355', 'license': 'Apache-2.0', 'n_tokens_mistral': 3473, 'n_tokens_neox': 3239, 'n_words': 925, 'path': 'D0/D0_example.ipynb', 'repo_name': 'yandexdataschool/mlhep-course-2017', 'size': '9338'}