<jupyter_start><jupyter_text>Measuring the Z boson massLet's look at a sample of $Z$ boson candidates recorded by CMS in 2011 and published at CERN opendata portal. It comes from DoubleMuon dataset with the following selection applied:- Both muons are "global" muons- invariant mass sits in range: 60 GeV < $ M_{\mu\mu}$ < 120 GeV- |$\eta$| < 2.1 for both muons- $p_{t}$ > 20 GeVThe following columns presented in the CSV file:- `Run`, Event are the run and event numbers, respectively- `pt` is the transverse momentum $p_{t}$ of the muon- `eta` is the pseudorapidity of the muon: $\eta$- `phi` is the $\phi$ angle of the muon direction- `Q` is the charge of the muon- `dxy` is the impact parameter in the transverse plane: $d_{xy}$ - how distant is the track from the collision point- `iso` is the track isolation: $I_{track}$ - how many other tracks are there aroung given track<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from skopt import gp_minimize
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Read dataset<jupyter_code>df = pd.read_csv('./Zmumu.csv')
df.shape<jupyter_output><empty_output><jupyter_text>Let's calculate the invariant mass $M$ of the two muons using the formula$M = \sqrt{2p_{t}^{1}p_{t}^{2}(\cosh(\eta_{1}-\eta_{2}) - \cos(\phi_{1}-\phi_{2}))}$<jupyter_code>df['M'] = np.sqrt(2 * df['pt1'] * df['pt2'] * (np.cosh(df['eta1'] - df['eta2']) - np.cos(df['phi1'] - df['phi2'])))
df.head(2)<jupyter_output><empty_output><jupyter_text>ModelThe distribution of the Z boson mass has the form of a normal distribution, in addition there is a noise. The distribution of the noise has an exponential form. Thus, the resulting model, it is a result of a superposition of two distributions - normal and polinomial Let's plot the distribution of Z boson mass<jupyter_code>def plot_mass(mass, bins_count=100):
    y, x = np.histogram(mass, bins=bins_count, density=False)
    err = np.sqrt(y)

    fig = plt.figure(figsize=(15,7))
    plt.title('Z mass', fontsize=20)
    plt.xlabel("$m_{\mu\mu}$ [GeV]", fontsize=20)
    plt.ylabel("Number of events", fontsize=20)
    plt.errorbar(x[:-1], y, yerr=err, fmt='o', color='red', ecolor='grey', capthick=0.5, zorder=1, label="data")
    return y, x
plot_mass(df.M);<jupyter_output><empty_output><jupyter_text>Exercise 1. clean up dataset a bit- demand that charge of muons should be opposite- $I_{track}$ < 3 and $d_{xy}$ < 0.2 cm<jupyter_code>df.describe()
### YOUR CODE GOES HERE ###
df_sign = df[((df.Q1 > 0) & (df.Q2 < 0)) | ((df.Q1 < 0) & (df.Q2 > 0))]
df_isolation = df_sign[(df_sign.iso1 < 3) & (df_sign.iso2 < 3) & (df_sign.dxy1.abs() < 0.2) & (df_sign.dxy2.abs() < 0.2)]
df_isolation.shape
df = df_isolation
plot_mass(df.M);<jupyter_output><empty_output><jupyter_text>Let's define parametrised modelit should represent mixture of 1) Gaussian signal and 2) background that for the simplicity we consider to be flat over mass. So it gives the following set of parameters:- m0 - center of the Gaussian- sigma - standard deviation of the Gaussian- ampl - height of the peak- bck - height of the background finding those parameters is called _fitting_ model into the data. It will be the goal for the rest of the exercise. For simplicity sake we'll stick with old good binned fit.<jupyter_code>def model_predict(params, X):
    m0, sigma, ampl, bck = params
    return bck + ampl / (sigma * np.sqrt(2 * np.pi)) * np.exp((-1) * (X - m0)**2 / (2 * sigma**2))
def model_loss(params, X, y):
#     y, x = np.histogram(mass, bins=bins_count, density=False)
#     residuals = model_predict(params, (x[1:] + x[:-1])/2) - y 
    residuals = y - model_predict(params, X)
    return np.sum(residuals**2) / len(residuals)
def plot_mass_with_model(params, mass, bins_count=100):
    y, X = plot_mass(mass, bins_count=bins_count)
    X = (X[1:] + X[:-1]) / 2
    error = model_loss(params, X, y)
    plt.plot(X, model_predict(params, X), color='blue', linewidth=3.0, zorder=2, label="fit, loss=%.2f" % error)
    plt.legend(fontsize='x-large')<jupyter_output><empty_output><jupyter_text>Here you can fit model parameters by hand<jupyter_code>plot_mass_with_model((91, 2.1, 4000, 18), df.M)<jupyter_output><empty_output><jupyter_text>... but you can do it automatically of course Setting up a scikit optimizer<jupyter_code>from tqdm import tqdm
from skopt import Optimizer

search_space = [(90.0, 92.0), # m0 range
    (1.5, 2.5), # sigma range
    (3500, 4500), # amplitude range
    (10, 20) # bck range
    ]
y, X = np.histogram(df.M, bins=120, density=False)
X = (X[1:] + X[:-1]) / 2
opt = Optimizer(search_space, base_estimator="GP", acq_func="EI", acq_optimizer="lbfgs")<jupyter_output><empty_output><jupyter_text>Running it for a while. You can re-run this cell several times<jupyter_code>from skopt.utils import create_result
for i in tqdm(range(50)):
    next_x = opt.ask()
    f_val = model_loss(next_x, X, y)
    opt.tell(next_x, f_val)
    
res = create_result(Xi=opt.Xi, yi=opt.yi, space=opt.space,
                         rng=opt.rng, models=opt.models)<jupyter_output>100%|██████████████████████████████████████████████████████████████████████████████████| 50/50 [00:46<00:00,  1.54s/it]<jupyter_text>A bit of search history<jupyter_code>import skopt.plots
skopt.plots.plot_convergence(res)
print (list(zip(["m0", "sigma", "ampl", "bck"], res.x)))<jupyter_output>[('m0', 90.73004571687005), ('sigma', 2.117595996516291), ('ampl', 3548), ('bck', 17)]<jupyter_text>Let's see how well the prediction fits the data<jupyter_code>plot_mass_with_model(res.x, df.M, bins_count=120)<jupyter_output><empty_output><jupyter_text>Grader part, do not change, please<jupyter_code>import grading
grader = grading.Grader(assignment_key="<KEY>", 
                      all_parts=["VI3xu", "VuE8x", "KzmMV", "TwZBF"])

ans_part1 = round(res.x[0])
grader.set_answer("VI3xu", ans_part1)

ans_part2 = round(res.x[1], 2)
grader.set_answer("VuE8x", ans_part2)

ans_part3 = round(res.x[3])
grader.set_answer("KzmMV", ans_part3)

ans_part4 = round(res.x[2])
grader.set_answer("TwZBF", ans_part4)
# token expires every 30 min
COURSERA_TOKEN = "<KEY>"### YOUR TOKEN HERE
COURSERA_EMAIL = "<EMAIL>"### YOUR EMAIL HERE
grader.submit(COURSERA_EMAIL, COURSERA_TOKEN)<jupyter_output>Submitted to Coursera platform. See results on assignment page!
{'fraction_non_alphanumeric': 0.1554993678887484, 'fraction_numerical': 0.033343868520859674, 'mean_word_length': 4.9877010406811735, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 52, 'https://': 0, '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': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '473336', 'license': 'MIT', 'n_tokens_mistral': 2536, 'n_tokens_neox': 2325, 'n_words': 644, 'path': 'Addressing Large Hadron Collider Challenges by Machine Learning/Week1/index.ipynb', 'repo_name': 'Mohitkr95/Advanced-ML', 'size': '6369'}