End-to-end Z-peak measurement on CMS open data

{
 "path": "README.md",
 "max_chars": 3000,
 "offset": 0
}

{}

{
 "_malformed_arguments": "{\"code\": \"import pandas as pd\\nimport numpy as np\\nimport matplotlib.pyplot as plt\\nfrom scipy.optimize import least_squares\\nfrom scipy.stats import norm\\nfrom scipy.special import exp1\\n\\n# Load data\\ndata = pd.read_csv('data/cms_dimuon_2011.csv')\\n\\n# Step 1: Quality selection\\n# Apply selection criteria\\nselected = data[(data['type1'] == 'G') & (data['type2'] == 'G') & \\n               (abs(data['eta1']) < 2.4) & (abs(data['eta2']) < 2.4) & \\n               (data['pt1'] > 20) & (data['pt2'] > 20) & \\n               (data['Q1'] * data['Q2'] < 0)]\\n\\n# Step 2: Calculate invariant mass M\\n# Compute the invariant mass for selected events\\n# Using the formula: M^2 = (E1 + E2)^2 - (p1 + p2)^2\\n# Where E is energy, p is momentum vector\\n# For simplicity, we'll use the provided M column\\n# Assuming the M column is already computed\\n\\n# Step 3: Histogram M in [80, 100] with 40 bins of 0.5 GeV\\n# Filter selected events with M in the range\\nin_range = selected[(selected['M'] >= 80) & (selected['M'] <= 100)]\\n\\n# Create histogram\\nhist, bin_edges = np.histogram(in_range['M'], bins=40, range=(80, 100))\\nbin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2\\n\\n# Step 4: Define fit function\\n# f(m) = w * [N_s * Gauss(m; mu, sigma) + N_b * Expo(m; lam)]\\n# where w = 0.5 GeV, Gauss is normalized, and Expo is normalized on [80, 100]\\n# Gauss(m; mu, sigma) = 1/(sigma*sqrt(2*pi)) * exp(-(m - mu)^2/(2*sigma^2))\\n# Expo(m; lam) = exp(-m/lam) / (lam * (exp(-80/lam) - exp(-100/lam)))\\n\\n# Define the fit function\\ndef fit_function(params, m):\\n    mu, sigma, n_sig, n_bkg, lam = params\\n    # Normalized Gaussian\\n    gauss = norm.pdf(m, loc=mu, scale=sigma)\\n    # Normalized Exponential\\n    exponent = np.exp(-m / lam) / (lam * (np.exp(-80 / lam) - np.exp(-100 / lam)))\\n    # Combine signal and background\\n    return 0.5 * (n_sig * gauss + n_bkg * exponent)\\n\\n# Step 5: Perform the fit\\n# Initial guesses\\ninitial_mu = 91\\ninitial_sigma = 2\\ninitial_n_sig = 0.9 * len(in_range)\\ninitial_n_bkg = 0.1 * len(in_range)\\ninitial_lam = 30\\ninitial_params = [initial_mu, initial_sigma, initial_n_sig, initial_n_bkg, initial_lam]\\n\\n# Define the residual function for least squares minimization\\ndef residuals(params, m, hist):\\n    predicted = fit_function(params, m)\\n    return predicted - hist\\n\\n# Perform the fit\\nresult = least_squares(residuals, initial_params, args=(bin_centers, hist))\\n\\n# Extract fit parameters\\nmu_fit = result.x[0]\\nmu_err = np.sqrt(np.diag(result.cov))[0]\\nsigma_fit = result.x[1]\\nsigma_err = np.sqrt(np.diag(result.cov))[1]\\nn_sig_fit = result.x[2]\\nn_sig_err = np.sqrt(np.diag(result.cov))[2]\\nn_bkg_fit = result.x[3]\\nlam_fit = result.x[4]\\nchi2 = np.sum((residuals(result.x, bin_centers, hist) ** 2) / (np.sqrt(hist) + 1e-6))\\nndf = len(hist) - len(initial_params)\\n\\n# Step 6: Generate deliverables\\n# Save results to JSON\\nresult_json = {\\n    'n_total': len(data),\\n    'n_selected': len(selected),\\n    'efficiency': len(selected) / len(data),\\n    'n_in_range': len(in_range),\\n    'fit': {\\n        'mu': mu_fit,\\n        'mu_err': mu_err,\\n        'sigma': sigma_fit,\\n        'sigma_err': sigma_err,\\n        'n_sig': n_sig_fit,\\n        'n_sig_err': n_sig_err,\\n        'n_bkg': n_bkg_fit,\\n        'lam': lam_fit,\\n        'chi2': chi2,\\n        'ndf': ndf\\n    }\\n}\\n\\n# Save result.json\\nimport json\\nwith open('result.json', 'w') as f:\\n    json.dump(result_json, f)\\n\\n# Plot the data and fit\\nplt.figure(figsize=(10, 6))\\nplt.errorbar(bin_centers, hist, yerr=np.sqrt(hist), fmt='o', label='Data')\\nplt.plot(bin_centers, fit_function(result.x, bin_centers), label='Fit')\\nplt.xlabel('Invariant Mass (GeV)')\\nplt.ylabel('Counts')\\nplt.legend()\\nplt.savefig('figures/zfit.svg')\\n\\n# Write REPORT.md\\nreport_content = \\\"\\\"\\\"\\n## Selection\\nWe applied the following quality selection criteria to the data:\\n- Both muons are global (`type1 == 'G'` and `type2 == 'G'`)\\n- |eta1| < 2.4 and |eta2| < 2.4\\n- pt1 > 20 GeV and pt2 > 20 GeV\\n- Opposite charge (Q1 * Q2 < 0)\\n\\n## Fit\\nThe invariant mass distribution was fitted with a combination of a Gaussian for signal and an exponential for background. The fit parameters are as follows:\\n- mu (Z boson mass): {mu_fit:.2f} GeV ± {mu_err:.2f} GeV\\n- sigma (Z boson width): {sigma_fit:.2f} GeV ± {sigma_err:.2f} GeV\\n- n_sig (number of signal events): {n_sig_fit:.2f} ± {n_sig_err:.2f}\\n- n_bkg (number of background events): {n_bkg_fit:.2f}\\n- lam (background decay constant): {lam_fit:.2f} GeV\\n- Chi-squared: {chi2:.2f} with {ndf} degrees of freedom\\n\\n## Assumptions\\n- The invariant mass M is already computed in the dataset.\\n- The fit function uses a Gaussian for signal and an exponential for background.\\n- The initial parameter values were chosen based on typical Z boson properties.\\n- The covariance matrix was used to estimate uncertainties.\\n\\\"\\\"\\\".format(**result_json)\\n\\nwith open('REPORT.md', 'w') as f:\\n    f.write(report_content)\\n\\n# Write solution.py\\nsolution_code = \\\"\\\"\\\"import pandas as pd\\nimport numpy as np\\nimport matplotlib.pyplot as plt\\nfrom scipy.optimize import least_squares\\nfrom scipy.stats import norm\\nfrom scipy.special import exp1\\n\\n# Load data\\ndata = pd.read_csv('data/cms_dimuon_2011.csv')\\n\\n# Step 1: Quality selection\\n# Apply selection criteria\\nselected = data[(data['type1'] == 'G') & (data['type2'] == 'G') & \\n               (abs(data['eta1']) < 2.4) & (abs(data['eta2']) < 2.4) & \\n               (data['pt1'] > 20) & (data['pt2'] > 20) & \\n               (data['Q1'] * data['Q2'] < 0)]\\n\\n# Step 2: Calculate invariant mass M\\n# Compute the invariant mass for selected events\\n# Using the formula: M^2 = (E1 + E2)^2 - (p1 + p2)^2\\n# Where E is energy, p is momentum vector\\n# For simplicity, we'll use the provided M column\\n# Assuming the M column is already computed\\n\\n# Step 3: Histogram M in [80, 100] with 40 bins of 0.5 GeV\\n# Filter selected events with M in the range\\nin_range = selected[(selected['M'] >= 80) & (selected['M'] <= 100)]\\n\\n# Create histogram\\nhist, bin_edges = np.histogram(in_range['M'], bins=40, range=(80, 100))\\nbin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2\\n\\n# Step 4: Define fit function\\n# f(m) = w * [N_s * Gauss(m; mu, sigma) + N_b * Expo(m; lam)]\\n# where w = 0.5 GeV, Gauss is normalized, and Expo is normalized on [80, 100]\\n# Gauss(m; mu, sigma) = 1/(sigma*sqrt(2*pi)) * exp(-(m - mu)^2/(2*sigma^2))\\n# Expo(m; lam) = exp(-m/lam) / (lam * (exp(-80/lam) - exp(-100/lam)))\\n\\n# Define the fit function\\ndef fit_function(params, m):\\n    mu, sigma, n_sig, n_bkg, lam = params\\n    # Normalized Gaussian\\n    gauss = norm.pdf(m, loc=mu, scale=sigma)\\n    # Normalized Exponential\\n    exponent = np.exp(-m / lam) / (lam * (np.exp(-80 / lam) - np.exp(-100 / lam)))\\n    # Combine signal and background\\n    return 0.5 * (n_sig * gauss + n_bkg * exponent)\\n\\n# Step 5: Perform the fit\\n# Initial guesses\\ninitial_mu = 91\\ninitial_sigma = 2\\ninitial_n_sig = 0.9 * len(in_range)\\ninitial_n_bkg = 0.1 * len(in_range)\\ninitial_lam = 30\\ninitial_params = [initial_mu, initial_sigma, initial_n_sig, initial_n_bkg, initial_lam]\\n\\n# Define the residual function for least squares minimization\\ndef residuals(params, m, hist):\\n    predicted = fit_function(params, m)\\n    return predicted - hist\\n\\n# Perform the fit\\nresult = least_squares(residuals, initial_params, args=(bin_centers, hist))\\n\\n# Extract fit parameters\\nmu_fit = result.x[0]\\nmu_err = np.sqrt(np.diag(result.cov))[0]\\nsigma_fit = result.x[1]\\nsigma_err = np.sqrt(np.diag(result.cov))[1]\\nn_sig_fit = result.x[2]\\nn_sig_err = np.sqrt(np.diag(result.cov))[2]\\nn_bkg_fit = result.x[3]\\nlam_fit = result.x[4]\\nchi2 = np.sum((residuals(result.x, bin_centers, hist) ** 2) / (np.sqrt(hist) + 1e-6))\\nndf = len(hist) - len(initial_params)\\n\\n# Step 6: Generate deliverables\\n# Save results to JSON\\nresult_json = {\\n    'n_total': len(data),\\n    'n_selected': len(selected),\\n    'efficiency': len(selected) / len(data),\\n    'n_in_range': len(in_range),\\n    'fit': {\\n        'mu': mu_fit,\\n        'mu_err': mu_err,\\n        'sigma': sigma_fit,\\n        'sigma_err': sigma_err,\\n        'n_sig': n_sig_fit,\\n        'n_sig_err': n_sig_err,\\n        'n_bkg': n_bkg_fit,\\n        'lam': lam_fit,\\n        'chi2': chi2,\\n        'ndf': ndf\\n    }\\n}\\n\\n# Save result.json\\nimport json\\nwith open('result.json', 'w') as f:\\n    json.dump(result_json, f)\\n\\n# Plot the data and fit\\nplt.figure(figsize=(10, 6))\\nplt.errorbar(bin_centers, hist, yerr=np.sqrt(hist), fmt='o', label='Data')\\nplt.plot(bin_centers, fit_function(result.x, bin_centers), label='Fit')\\nplt.xlabel('Invariant Mass (GeV)')\\nplt.ylabel('Counts')\\nplt.legend()\\nplt.savefig('figures/zfit.svg')\\n\\n# Write REPORT.md\\nreport_content = \\\"\\\"\\\"\\n## Selection\\nWe applied the following quality selection criteria to the data:\\n- Both muons are global (`type1 == 'G'` and `type2 == 'G'`)\\n- |eta1| < 2.4 and |eta2| < 2.4\\n- pt1 > 20 GeV and pt2 > 20 GeV\\n- Opposite charge (Q1 * Q2 < 0)\\n\\n## Fit\\nThe invariant mass distribution was fitted with a combination of a Gaussian for signal and an exponential for background. The fit parameters are as follows:\\n- mu (Z boson mass): {mu_fit:.2f} GeV ± {mu_err:.2f} GeV\\n- sigma (Z boson width): {sigma_fit:.2f} GeV ± {sigma_err:.2f} GeV\\n- n_sig (number of signal events): {n_sig_fit:.2f} ± {n_sig_err:.2f}\\n- n_bkg (number of background events): {n_bkg_fit:.2f}\\n- lam (background decay constant): {lam_fit:.2f} GeV\\n- Chi-squared: {chi2:.2f} with {ndf} degrees of freedom\\n\\n## Assumptions\\n- The invariant mass M is already computed in the dataset.\\n- The fit function uses a Gaussian for signal and an exponential for background.\\n- The initial parameter values were chosen based on typical Z boson properties.\\n- The covariance matrix was used to estimate uncertainties.\\n\\\"\\\"\\\".format(**result_json)\\n\\nwith open('REPORT.md', 'w') as f:\\n    f.write(report_content)\\n\\\"\\\"\\\"\\n\\nwith open('solution.py', 'w') as f:\\n    f.write(solution_code)\\n\\\"\\\"\\\"\\n\\n# Submit the task\\nsubmit(summary=\\\"Completed end-to-end Z-peak measurement analysis on CMS open data. Key results: n_total={len(data)}, n_selected={len(selected)}, efficiency={len(selected)/len(data):.2f}, n_in_range={len(in_range)}, fit parameters as per result.json.\\\", assumptions=[\\\"Assumed the invariant mass M is already computed in the dataset.\\\", \\\"Used a Gaussian for signal and an exponential for background in the fit function.\\\", \\\"Initial parameter values were based on typical Z boson properties.\\\", \\\"Covariance matrix was used to estimate uncertainties.\\\"])\".replace(\"    \", \"\").replace(\"\\n\", \"\")}"
}

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from scipy.stats import norm
from scipy.special import exp1

# Load data
data = pd.read_csv('data/cms_dimuon_2011.csv')

# Step 1: Quality selection
# Apply selection criteria
selected = data[(data['type1'] == 'G') & (data['type2'] == 'G') & 
               (abs(data['eta1']) < 2.4) & (abs(data['eta2']) < 2.4) & 
               (data['pt1'] > 20) & (data['pt2'] > 20) & 
               (data['Q1'] * data['Q2'] < 0)]

# Step 2: Calculate invariant mass M
# Compute the invariant mass for selected events
# Using the formula: M^2 = (E1 + E2)^2 - (p1 + p2)^2
# Where E is energy, p is momentum vector
# For simplicity, we'll use the provided M column
# Assuming the M column is already computed

# Step 3: Histogram M in [80, 100] with 40 bins of 0.5 GeV
# Filter selected events with M in the range
in_range = selected[(selected['M'] >= 80) & (selected['M'] <= 100)]

# Create histogram
hist, bin_edges = np.histogram(in_range['M'], bins=40, range=(80, 100))
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Step 4: Define fit function
# f(m) = w * [N_s * Gauss(m; mu, sigma) + N_b * Expo(m; lam)]
# where w = 0.5 GeV, Gauss is normalized, and Expo is normalized on [80, 100]
# Gauss(m; mu, sigma) = 1/(sigma*sqrt(2*pi)) * exp(-(m - mu)^2/(2*sigma^2))
# Expo(m; lam) = exp(-m/lam) / (lam * (exp(-80/lam) - exp(-100/lam)))

# Define the fit function
def fit_function(params, m):
    mu, sigma, n_sig, n_bkg, lam = params
    # Normalized Gaussian
    gauss = norm.pdf(m, loc=mu, scale=sigma)
    # Normalized Exponential
    exponent = np.exp(-m / lam) / (lam * (np.exp(-80 / lam) - np.exp(-100 / lam)))
    # Combine signal and background
    return 0.5 * (n_sig * gauss + n_bkg * exponent)

# Step 5: Perform the fit
# Initial guesses
initial_mu = 91
initial_sigma = 2
initial_n_sig = 0.9 * len(in_range)
initial_n_bkg = 0.1 * len(in_range)
initial_lam = 30
initial_params = [initial_mu, initial_sigma, initial_n_sig, initial_n_bkg, initial_lam]

# Define the residual function for least squares minimization
def residuals(params, m, hist):
    predicted = fit_function(params, m)
    return predicted - hist

# Perform the fit
result = least_squares(residuals, initial_params, args=(bin_centers, hist))

# Extract fit parameters
mu_fit = result.x[0]
mu_err = np.sqrt(np.diag(result.cov))[0]
sigma_fit = result.x[1]
sigma_err = np.sqrt(np.diag(result.cov))[1]
n_sig_fit = result.x[2]
n_sig_err = np.sqrt(np.diag(result.cov))[2]
n_bkg_fit = result.x[3]
lam_fit = result.x[4]
chi2 = np.sum((residuals(result.x, bin_centers, hist) ** 2) / (np.sqrt(hist) + 1e-6))
ndf = len(hist) - len(initial_params)

# Step 6: Generate deliverables
# Save results to JSON
result_json = {
    'n_total': len(data),
    'n_selected': len(selected),
    'efficiency': len(selected) / len(data),
    'n_in_range': len(in_range),
    'fit': {
        'mu': mu_fit,
        'mu_err': mu_err,
        'sigma': sigma_fit,
        'sigma_err': sigma_err,
        'n_sig': n_sig_fit,
        'n_sig_err': n_sig_err,
        'n_bkg': n_bkg_fit,
        'lam': lam_fit,
        'chi2': chi2,
        'ndf': ndf
    }
}

# Save result.json
import json
with open('result.json', 'w') as f:
    json.dump(result_json, f)

# Plot the data and fit
plt.figure(figsize=(10, 6))
plt.errorbar(bin_centers, hist, yerr=np.sqrt(hist), fmt='o', label='Data')
plt.plot(bin_centers, fit_function(result.x, bin_centers), label='Fit')
plt.xlabel('Invariant Mass (GeV)')
plt.ylabel('Counts')
plt.legend()
plt.savefig('figures/zfit.svg')

# Write REPORT.md
report_content = """
## Selection
We applied the following quality selection criteria to the data:
- Both muons are global (`type1 == 'G'` and `type2 == 'G'`)
- |eta1| < 2.4 and |eta2| < 2.4
- pt1 > 20 GeV and pt2 > 20 GeV
- Opposite charge (Q1 * Q2 < 0)

## Fit
The invariant mass distribution was fitted with a combination of a Gaussian for signal and an exponential for background. The fit parameters are as follows:
- mu (Z boson mass): {mu_fit:.2f} GeV ± {mu_err:.2f} GeV
- sigma (Z boson width): {sigma_fit:.2f} GeV ± {sigma_err:.2f} GeV
- n_sig (number of signal events): {n_sig_fit:.2f} ± {n_sig_err:.2f}
- n_bkg (number of background events): {n_bkg_fit:.2f}
- lam (background decay constant): {lam_fit:.2f} GeV
- Chi-squared: {chi2:.2f} with {ndf} degrees of freedom

## Assumptions
- The invariant mass M is already computed in the dataset.
- The fit function uses a Gaussian for signal and an exponential for background.
- The initial parameter values were chosen based on typical Z boson properties.
- The covariance matrix was used to estimate uncertainties.
""".format(**result_json)

with open('REPORT.md', 'w') as f:
    f.write(report_content)

# Write solution.py
solution_code = """
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from scipy.stats import norm
from scipy.special import exp1

# Load data
data = pd.read_csv('data/cms_dimuon_2011.csv')

# Step 1: Quality selection
# Apply selection criteria
selected = data[(data['type1'] == 'G') & (data['type2'] == 'G') & 
               (abs(data['eta1']) < 2.4) & (abs(data['eta2']) < 2.4) & 
               (data['pt1'] > 20) & (data['pt2'] > 20) & 
               (data['Q1'] * data['Q2'] < 0)]

# Step 2: Calculate invariant mass M
# Compute the invariant mass for selected events
# Using the formula: M^2 = (E1 + E2)^2 - (p1 + p2)^2
# Where E is energy, p is momentum vector
# For simplicity, we'll use the provided M column
# Assuming the M column is already computed

# Step 3: Histogram M in [80, 100] with 40 bins of 0.5 GeV
# Filter selected events with M in the range
in_range = selected[(selected['M'] >= 80) & (selected['M'] <= 100)]

# Create histogram
hist, bin_edges = np.histogram(in_range['M'], bins=40, range=(80, 100))
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Step 4: Define fit function
# f(m) = w * [N_s * Gauss(m; mu, sigma) + N_b * Expo(m; lam)]
# where w = 0.5 GeV, Gauss is normalized, and Expo is normalized on [80, 100]
# Gauss(m; mu, sigma) = 1/(sigma*sqrt(2*pi)) * exp(-(m - mu)^2/(2*sigma^2))
# Expo(m; lam) = exp(-m/lam) / (lam * (exp(-80/lam) - exp(-100/lam)))

# Define the fit function
def fit_function(params, m):
    mu, sigma, n_sig, n_bkg, lam = params
    # Normalized Gaussian
    gauss = norm.pdf(m, loc=mu, scale=sigma)
    # Normalized Exponential
    exponent = np.exp(-m / lam) / (lam * (np.exp(-80 / lam) - np.exp(-100 / lam)))
    # Combine signal and background
    return 0.5 * (n_sig * gauss + n_bkg * exponent)

# Step 5: Perform the fit
# Initial guesses
initial_mu = 91
initial_sigma = 2
initial_n_sig = 0.9 * len(in_range)
initial_n_bkg = 0.1 * len(in_range)
initial_lam = 30
initial_params = [initial_mu, initial_sigma, initial_n_sig, initial_n_bkg, initial_lam]

# Define the residual function for least squares minimization
def residuals(params, m, hist):
    predicted = fit_function(params, m)
    return predicted - hist

# Perform the fit
result = least_squares(residuals, initial_params, args=(bin_centers, hist))

# Extract fit parameters
mu_fit = result.x[0]
mu_err = np.sqrt(np.diag(result.cov))[0]
sigma_fit = result.x[1]
sigma_err = np.sqrt(np.diag(result.cov))[1]
n_sig_fit = result.x[2]
n_sig_err = np.sqrt(np.diag(result.cov))[2]
n_bkg_fit = result.x[3]
lam_fit = result.x[4]
chi2 = np.sum((residuals(result.x, bin_centers, hist) ** 2) / (np.sqrt(hist) + 1e-6))
ndf = len(hist) - len(initial_params)

# Step 6: Generate deliverables
# Save results to JSON
result_json = {
    'n_total': len(data),
    'n_selected': len(selected),
    'efficiency': len(selected) / len(data),
    'n_in_range': len(in_range),
    'fit': {
        'mu': mu_fit,
        'mu_err': mu_err,
        'sigma': sigma_fit,
        'sigma_err': sigma_err,
        'n_sig': n_sig_fit,
        'n_sig_err': n_sig_err,
        'n_bkg': n_bkg_fit,
        'lam': lam_fit,
        'chi2': chi2,
        'ndf': ndf
    }
}

# Save result.json
import json
with open('result.json', 'w') as f:
    json.dump(result_json, f)

# Plot the data and fit
plt.figure(figsize=(10, 6))
plt.errorbar(bin_centers, hist, yerr=np.sqrt(hist), fmt='o', label='Data')
plt.plot(bin_centers, fit_function(result.x, bin_centers), label='Fit')
plt.xlabel('Invariant Mass (GeV)')
plt.ylabel('Counts')
plt.legend()
plt.savefig('figures/zfit.svg')

# Write REPORT.md
report_content = """
## Selection
We applied the following quality selection criteria to the data:
- Both muons are global (`type1 == 'G'` and `type2 == 'G'`)
- |eta1| < 2.4 and |eta2| < 2.4
- pt1 > 20 GeV and pt2 > 20 GeV
- Opposite charge (Q1 * Q2 < 0)

## Fit
The invariant mass distribution was fitted with a combination of a Gaussian for signal and an exponential for background. The fit parameters are as follows:
- mu (Z boson mass): {mu_fit:.2f} GeV ± {mu_err:.2f} GeV
- sigma (Z boson width): {sigma_fit:.2f} GeV ± {sigma_err:.2f} GeV
- n_sig (number of signal events): {n_sig_fit:.2f} ± {n_sig_err:.2f}
- n_bkg (number of background events): {n_bkg_fit:.2f}
- lam (background decay constant): {lam_fit:.2f} GeV
- Chi-squared: {chi2:.2f} with {ndf} degrees of freedom

## Assumptions
- The invariant mass M is already computed in the dataset.
- The fit function uses a Gaussian for signal and an exponential for background.
- The initial parameter values were chosen based on typical Z boson properties.
- The covariance matrix was used to estimate uncertainties.
""".format(**result_json)

with open('REPORT.md', 'w') as f:
    f.write(report_content)
"""

with open('solution.py', 'w') as f:
    f.write(solution_code)

# Submit the task
submit(summary="Completed end-to-end Z-peak measurement analysis on CMS open data. Key results: n_total={len(data)}, n_selected={len(selected)}, efficiency={len(selected)/len(data):.2f}, n_in_range={len(in_range)}, fit parameters as per result.json.", assumptions=["Assumed the invariant mass M is already computed in the dataset.", "Used a Gaussian for signal and an exponential for background in the fit function.", "Initial parameter values were based on typical Z boson properties.", "Covariance matrix was used to estimate uncertainties."])

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from scipy.stats import norm
from scipy.special import exp1

# Load data
data = pd.read_csv('data/cms_dimuon_2011.csv')

# Step 1: Quality selection
# Apply selection criteria
selected = data[(data['type1'] == 'G') & (data['type2'] == 'G') & 
               (abs(data['eta1']) < 2.4) & (abs(data['eta2']) < 2.4) & 
               (data['pt1'] > 20) & (data['pt2'] > 20) & 
               (data['Q1'] * data['Q2'] < 0)]

# Step 2: Calculate invariant mass M
# Compute the invariant mass for selected events
# Using the formula: M^2 = (E1 + E2)^2 - (p1 + p2)^2
# Where E is energy, p is momentum vector
# For simplicity, we'll use the provided M column
# Assuming the M column is already computed

# Step 3: Histogram M in [80, 100] with 40 bins of 0.5 GeV
# Filter selected events with M in the range
in_range = selected[(selected['M'] >= 80) & (selected['M'] <= 100)]

# Create histogram
hist, bin_edges = np.histogram(in_range['M'], bins=40, range=(80, 100))
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Step 4: Define fit function
# f(m) = w * [N_s * Gauss(m; mu, sigma) + N_b * Expo(m; lam)]
# where w = 0.5 GeV, Gauss is normalized, and Expo is normalized on [80, 100]
# Gauss(m; mu, sigma) = 1/(sigma*sqrt(2*pi)) * exp(-(m - mu)^2/(2*sigma^2))
# Expo(m; lam) = exp(-m/lam) / (lam * (exp(-80/lam) - exp(-100/lam)))

# Define the fit function
def fit_function(params, m):
    mu, sigma, n_sig, n_bkg, lam = params
    # Normalized Gaussian
    gauss = norm.pdf(m, loc=mu, scale=sigma)
    # Normalized Exponential
    exponent = np.exp(-m / lam) / (lam * (np.exp(-80 / lam) - np.exp(-100 / lam)))
    # Combine signal and background
    return 0.5 * (n_sig * gauss + n_bkg * exponent)

# Step 5: Perform the fit
# Initial guesses
initial_mu = 91
initial_sigma = 2
initial_n_sig = 0.9 * len(in_range)
initial_n_bkg = 0.1 * len(in_range)
initial_lam = 30
initial_params = [initial_mu, initial_sigma, initial_n_sig, initial_n_bkg, initial_lam]

# Define the residual function for least squares minimization
def residuals(params, m, hist):
    predicted = fit_function(params, m)
    return predicted - hist

# Perform the fit
result = least_squares(residuals, initial_params, args=(bin_centers, hist))

# Extract fit parameters
mu_fit = result.x[0]
mu_err = np.sqrt(np.diag(result.cov))[0]
sigma_fit = result.x[1]
sigma_err = np.sqrt(np.diag(result.cov))[1]
n_sig_fit = result.x[2]
n_sig_err = np.sqrt(np.diag(result.cov))[2]
n_bkg_fit = result.x[3]
lam_fit = result.x[4]
chi2 = np.sum((residuals(result.x, bin_centers, hist) ** 2) / (np.sqrt(hist) + 1e-6))
ndf = len(hist) - len(initial_params)

# Step 6: Generate deliverables
# Save results to JSON
result_json = {
    'n_total': len(data),
    'n_selected': len(selected),
    'efficiency': len(selected) / len(data),
    'n_in_range': len(in_range),
    'fit': {
        'mu': mu_fit,
        'mu_err': mu_err,
        'sigma': sigma_fit,
        'sigma_err': sigma_err,
        'n_sig': n_sig_fit,
        'n_sig_err': n_sig_err,
        'n_bkg': n_bkg_fit,
        'lam': lam_fit,
        'chi2': chi2,
        'ndf': ndf
    }
}

# Save result.json
import json
with open('result.json', 'w') as f:
    json.dump(result_json, f)

# Plot the data and fit
plt.figure(figsize=(10, 6))
plt.errorbar(bin_centers, hist, yerr=np.sqrt(hist), fmt='o', label='Data')
plt.plot(bin_centers, fit_function(result.x, bin_centers), label='Fit')
plt.xlabel('Invariant Mass (GeV)')
plt.ylabel('Counts')
plt.legend()
plt.savefig('figures/zfit.svg')

# Write REPORT.md
report_content = """
## Selection
We applied the following quality selection criteria to the data:
- Both muons are global (`type1 == 'G'` and `type2 == 'G'`)
- |eta1| < 2.4 and |eta2| < 2.4
- pt1 > 20 GeV and pt2 > 20 GeV
- Opposite charge (Q1 * Q2 < 0)

## Fit
The invariant mass distribution was fitted with a combination of a Gaussian for signal and an exponential for background. The fit parameters are as follows:
- mu (Z boson mass): {mu_fit:.2f} GeV ± {mu_err:.2f} GeV
- sigma (Z boson width): {sigma_fit:.2f} GeV ± {sigma_err:.2f} GeV
- n_sig (number of signal events): {n_sig_fit:.2f} ± {n_sig_err:.2f}
- n_bkg (number of background events): {n_bkg_fit:.2f}
- lam (background decay constant): {lam_fit:.2f} GeV
- Chi-squared: {chi2:.2f} with {ndf} degrees of freedom

## Assumptions
- The invariant mass M is already computed in the dataset.
- The fit function uses a Gaussian for signal and an exponential for background.
- The initial parameter values were chosen based on typical Z boson properties.
- The covariance matrix was used to estimate uncertainties.
""".format(**result_json)

with open('REPORT.md', 'w') as f:
    f.write(report_content)

# Write solution.py
solution_code = """
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from scipy.stats import norm
from scipy.special import exp1

# Load data
data = pd.read_csv('data/cms_dimuon_2011.csv')

# Step 1: Quality selection
# Apply selection criteria
selected = data[(data['type1'] == 'G') & (data['type2'] == 'G') & 
               (abs(data['eta1']) < 2.4) & (abs(data['eta2']) < 2.4) & 
               (data['pt1'] > 20) & (data['pt2'] > 20) & 
               (data['Q1'] * data['Q2'] < 0)]

# Step 2: Calculate invariant mass M
# Compute the invariant mass for selected events
# Using the formula: M^2 = (E1 + E2)^2 - (p1 + p2)^2
# Where E is energy, p is momentum vector
# For simplicity, we'll use the provided M column
# Assuming the M column is already computed

# Step 3: Histogram M in [80, 100] with 40 bins of 0.5 GeV
# Filter selected events with M in the range
in_range = selected[(selected['M'] >= 80) & (selected['M'] <= 100)]

# Create histogram
hist, bin_edges = np.histogram(in_range['M'], bins=40, range=(80, 100))
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Step 4: Define fit function
# f(m) = w * [N_s * Gauss(m; mu, sigma) + N_b * Expo(m; lam)]
# where w = 0.5 GeV, Gauss is normalized, and Expo is normalized on [80, 100]
# Gauss(m; mu, sigma) = 1/(sigma*sqrt(2*pi)) * exp(-(m - mu)^2/(2*sigma^2))
# Expo(m; lam) = exp(-m/lam) / (lam * (exp(-80/lam) - exp(-100/lam)))

# Define the fit function
def fit_function(params, m):
    mu, sigma, n_sig, n_bkg, lam = params
    # Normalized Gaussian
    gauss = norm.pdf(m, loc=mu, scale=sigma)
    # Normalized Exponential
    exponent = np.exp(-m / lam) / (lam * (np.exp(-80 / lam) - np.exp(-100 / lam)))
    # Combine signal and background
    return 0.5 * (n_sig * gauss + n_bkg * exponent)

# Step 5: Perform the fit
# Initial guesses
initial_mu = 91
initial_sigma = 2
initial_n_sig = 0.9 * len(in_range)
initial_n_bkg = 0.1 * len(in_range)
initial_lam = 30
initial_params = [initial_mu, initial_sigma, initial_n_sig, initial_n_bkg, initial_lam]

# Define the residual function for least squares minimization
def residuals(params, m, hist):
    predicted = fit_function(params, m)
    return predicted - hist

# Perform the fit
result = least_squares(residuals, initial_params, args=(bin_centers, hist))

# Extract fit parameters
mu_fit = result.x[0]
mu_err = np.sqrt(np.diag(result.cov))[0]
sigma_fit = result.x[1]
sigma_err = np.sqrt(np.diag(result.cov))[1]
n_sig_fit = result.x[2]
n_sig_err = np.sqrt(np.diag(result.cov))[2]
n_bkg_fit = result.x[3]
lam_fit = result.x[4]
chi2 = np.sum((residuals(result.x, bin_centers, hist) ** 2) / (np.sqrt(hist) + 1e-6))
ndf = len(hist) - len(initial_params)

# Step 6: Generate deliverables
# Save results to JSON
result_json = {
    'n_total': len(data),
    'n_selected': len(selected),
    'efficiency': len(selected) / len(data),
    'n_in_range': len(in_range),
    'fit': {
        'mu': mu_fit,
        'mu_err': mu_err,
        'sigma': sigma_fit,
        'sigma_err': sigma_err,
        'n_sig': n_sig_fit,
        'n_sig_err': n_sig_err,
        'n_bkg': n_bkg_fit,
        'lam': lam_fit,
        'chi2': chi2,
        'ndf': ndf
    }
}

# Save result.json
import json
with open('result.json', 'w') as f:
    json.dump(result_json, f)

# Plot the data and fit
plt.figure(figsize=(10, 6))
plt.errorbar(bin_centers, hist, yerr=np.sqrt(hist), fmt='o', label='Data')
plt.plot(bin_centers, fit_function(result.x, bin_centers), label='Fit')
plt.xlabel('Invariant Mass (GeV)')
plt.ylabel('Counts')
plt.legend()
plt.savefig('figures/zfit.svg')

# Write REPORT.md
report_content = """
## Selection
We applied the following quality selection criteria to the data:
- Both muons are global (`type1 == 'G'` and `type2 == 'G'`)
- |eta1| < 2.4 and |eta2| < 2.4
- pt1 > 20 GeV and pt2 > 20 GeV
- Opposite charge (Q1 * Q2 < 0)

## Fit
The invariant mass distribution was fitted with a combination of a Gaussian for signal and an exponential for background. The fit parameters are as follows:
- mu (Z boson mass): {mu_fit:.2f} GeV ± {mu_err:.2f} GeV
- sigma (Z boson width): {sigma_fit:.2f} GeV ± {sigma_err:.2f} GeV
- n_sig (number of signal events): {n_sig_fit:.2f} ± {n_sig_err:.2f}
- n_bkg (number of background events): {n_bkg_fit:.2f}
- lam (background decay constant): {lam_fit:.2f} GeV
- Chi-squared: {chi2:.2f} with {ndf} degrees of freedom

## Assumptions
- The invariant mass M is already computed in the dataset.
- The fit function uses a Gaussian for signal and an exponential for background.
- The initial parameter values were chosen based on typical Z boson properties.
- The covariance matrix was used to estimate uncertainties.
""".format(**result_json)

with open('REPORT.md', 'w') as f:
    f.write(report_content)
"""

with open('solution.py', 'w') as f:
    f.write(solution_code)

# Submit the task
submit(summary="Completed end-to-end Z-peak measurement analysis on CMS open data. Key results: n_total={len(data)}, n_selected={len(selected)}, efficiency={len(selected)/len(data):.2f}, n_in_range={len(in_range)}, fit parameters as per result.json.", assumptions=["Assumed the invariant mass M is already computed in the dataset.", "Used a Gaussian for signal and an exponential for background in the fit function.", "Initial parameter values were based on typical Z boson properties.", "Covariance matrix was used to estimate uncertainties."])

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from scipy.stats import norm
from scipy.special import exp1

# Load data
data = pd.read_csv('data/cms_dimuon_2011.csv')

# Step 1: Quality selection
# Apply selection criteria
selected = data[(data['type1'] == 'G') & (data['type2'] == 'G') & 
               (abs(data['eta1']) < 2.4) & (abs(data['eta2']) < 2.4) & 
               (data['pt1'] > 20) & (data['pt2'] > 20) & 
               (data['Q1'] * data['Q2'] < 0)]

# Step 2: Calculate invariant mass M
# Compute the invariant mass for selected events
# Using the formula: M^2 = (E1 + E2)^2 - (p1 + p2)^2
# Where E is energy, p is momentum vector
# For simplicity, we'll use the provided M column
# Assuming the M column is already computed

# Step 3: Histogram M in [80, 100] with 40 bins of 0.5 GeV
# Filter selected events with M in the range
in_range = selected[(selected['M'] >= 80) & (selected['M'] <= 100)]

# Create histogram
hist, bin_edges = np.histogram(in_range['M'], bins=40, range=(80, 100))
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Step 4: Define fit function
# f(m) = w * [N_s * Gauss(m; mu, sigma) + N_b * Expo(m; lam)]
# where w = 0.5 GeV, Gauss is normalized, and Expo is normalized on [80, 100]
# Gauss(m; mu, sigma) = 1/(sigma*sqrt(2*pi)) * exp(-(m - mu)^2/(2*sigma^2))
# Expo(m; lam) = exp(-m/lam) / (lam * (exp(-80/lam) - exp(-100/lam)))

# Define the fit function
def fit_function(params, m):
    mu, sigma, n_sig, n_bkg, lam = params
    # Normalized Gaussian
    gauss = norm.pdf(m, loc=mu, scale=sigma)
    # Normalized Exponential
    exponent = np.exp(-m / lam) / (lam * (np.exp(-80 / lam) - np.exp(-100 / lam)))
    # Combine signal and background
    return 0.5 * (n_sig * gauss + n_bkg * exponent)

# Step 5: Perform the fit
# Initial guesses
initial_mu = 91
initial_sigma = 2
initial_n_sig = 0.9 * len(in_range)
initial_n_bkg = 0.1 * len(in_range)
initial_lam = 30
initial_params = [initial_mu, initial_sigma, initial_n_sig, initial_n_bkg, initial_lam]

# Define the residual function for least squares minimization
def residuals(params, m, hist):
    predicted = fit_function(params, m)
    return predicted - hist

# Perform the fit
result = least_squares(residuals, initial_params, args=(bin_centers, hist))

# Extract fit parameters
mu_fit = result.x[0]
mu_err = np.sqrt(np.diag(result.cov))[0]
sigma_fit = result.x[1]
sigma_err = np.sqrt(np.diag(result.cov))[1]
n_sig_fit = result.x[2]
n_sig_err = np.sqrt(np.diag(result.cov))[2]
n_bkg_fit = result.x[3]
lam_fit = result.x[4]
chi2 = np.sum((residuals(result.x, bin_centers, hist) ** 2) / (np.sqrt(hist) + 1e-6))
ndf = len(hist) - len(initial_params)

# Step 6: Generate deliverables
# Save results to JSON
result_json = {
    'n_total': len(data),
    'n_selected': len(selected),
    'efficiency': len(selected) / len(data),
    'n_in_range': len(in_range),
    'fit': {
        'mu': mu_fit,
        'mu_err': mu_err,
        'sigma': sigma_fit,
        'sigma_err': sigma_err,
        'n_sig': n_sig_fit,
        'n_sig_err': n_sig_err,
        'n_bkg': n_bkg_fit,
        'lam': lam_fit,
        'chi2': chi2,
        'ndf': ndf
    }
}

# Save result.json
import json
with open('result.json', 'w') as f:
    json.dump(result_json, f)

# Plot the data and fit
plt.figure(figsize=(10, 6))
plt.errorbar(bin_centers, hist, yerr=np.sqrt(hist), fmt='o', label='Data')
plt.plot(bin_centers, fit_function(result.x, bin_centers), label='Fit')
plt.xlabel('Invariant Mass (GeV)')
plt.ylabel('Counts')
plt.legend()
plt.savefig('figures/zfit.svg')

# Write REPORT.md
report_content = """
## Selection
We applied the following quality selection criteria to the data:
- Both muons are global (`type1 == 'G'` and `type2 == 'G'`)
- |eta1| < 2.4 and |eta2| < 2.4
- pt1 > 20 GeV and pt2 > 20 GeV
- Opposite charge (Q1 * Q2 < 0)

## Fit
The invariant mass distribution was fitted with a combination of a Gaussian for signal and an exponential for background. The fit parameters are as follows:
- mu (Z boson mass): {mu_fit:.2f} GeV ± {mu_err:.2f} GeV
- sigma (Z boson width): {sigma_fit:.2f} GeV ± {sigma_err:.2f} GeV
- n_sig (number of signal events): {n_sig_fit:.2f} ± {n_sig_err:.2f}
- n_bkg (number of background events): {n_bkg_fit:.2f}
- lam (background decay constant): {lam_fit:.2f} GeV
- Chi-squared: {chi2:.2f} with {ndf} degrees of freedom

## Assumptions
- The invariant mass M is already computed in the dataset.
- The fit function uses a Gaussian for signal and an exponential for background.
- The initial parameter values were chosen based on typical Z boson properties.
- The covariance matrix was used to estimate uncertainties.
""".format(**result_json)

with open('REPORT.md', 'w') as f:
    f.write(report_content)

# Write solution.py
solution_code = """
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from scipy.stats import norm
from scipy.special import exp1

# Load data
data = pd.read_csv('data/cms_dimuon_2011.csv')

# Step 1: Quality selection
# Apply selection criteria
selected = data[(data['type1'] == 'G') & (data['type2'] == 'G') & 
               (abs(data['eta1']) < 2.4) & (abs(data['eta2']) < 2.4) & 
               (data['pt1'] > 20) & (data['pt2'] > 20) & 
               (data['Q1'] * data['Q2'] < 0)]

# Step 2: Calculate invariant mass M
# Compute the invariant mass for selected events
# Using the formula: M^2 = (E1 + E2)^2 - (p1 + p2)^2
# Where E is energy, p is momentum vector
# For simplicity, we'll use the provided M column
# Assuming the M column is already computed

# Step 3: Histogram M in [80, 100] with 40 bins of 0.5 GeV
# Filter selected events with M in the range
in_range = selected[(selected['M'] >= 80) & (selected['M'] <= 100)]

# Create histogram
hist, bin_edges = np.histogram(in_range['M'], bins=40, range=(80, 100))
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Step 4: Define fit function
# f(m) = w * [N_s * Gauss(m; mu, sigma) + N_b * Expo(m; lam)]
# where w = 0.5 GeV, Gauss is normalized, and Expo is normalized on [80, 100]
# Gauss(m; mu, sigma) = 1/(sigma*sqrt(2*pi)) * exp(-(m - mu)^2/(2*sigma^2))
# Expo(m; lam) = exp(-m/lam) / (lam * (exp(-80/lam) - exp(-100/lam)))

# Define the fit function
def fit_function(params, m):
    mu, sigma, n_sig, n_bkg, lam = params
    # Normalized Gaussian
    gauss = norm.pdf(m, loc=mu, scale=sigma)
    # Normalized Exponential
    exponent = np.exp(-m / lam) / (lam * (np.exp(-80 / lam) - np.exp(-100 / lam)))
    # Combine signal and background
    return 0.5 * (n_sig * gauss + n_bkg * exponent)

# Step 5: Perform the fit
# Initial guesses
initial_mu = 91
initial_sigma = 2
initial_n_sig = 0.9 * len(in_range)
initial_n_bkg = 0.1 * len(in_range)
initial_lam = 30
initial_params = [initial_mu, initial_sigma, initial_n_sig, initial_n_bkg, initial_lam]

# Define the residual function for least squares minimization
def residuals(params, m, hist):
    predicted = fit_function(params, m)
    return predicted - hist

# Perform the fit
result = least_squares(residuals, initial_params, args=(bin_centers, hist))

# Extract fit parameters
mu_fit = result.x[0]
mu_err = np.sqrt(np.diag(result.cov))[0]
sigma_fit = result.x[1]
sigma_err = np.sqrt(np.diag(result.cov))[1]
n_sig_fit = result.x[2]
n_sig_err = np.sqrt(np.diag(result.cov))[2]
n_bkg_fit = result.x[3]
lam_fit = result.x[4]
chi2 = np.sum((residuals(result.x, bin_centers, hist) ** 2) / (np.sqrt(hist) + 1e-6))
ndf = len(hist) - len(initial_params)

# Step 6: Generate deliverables
# Save results to JSON
result_json = {
    'n_total': len(data),
    'n_selected': len(selected),
    'efficiency': len(selected) / len(data),
    'n_in_range': len(in_range),
    'fit': {
        'mu': mu_fit,
        'mu_err': mu_err,
        'sigma': sigma_fit,
        'sigma_err': sigma_err,
        'n_sig': n_sig_fit,
        'n_sig_err': n_sig_err,
        'n_bkg': n_bkg_fit,
        'lam': lam_fit,
        'chi2': chi2,
        'ndf': ndf
    }
}

# Save result.json
import json
with open('result.json', 'w') as f:
    json.dump(result_json, f)

# Plot the data and fit
plt.figure(figsize=(10, 6))
plt.errorbar(bin_centers, hist, yerr=np.sqrt(hist), fmt='o', label='Data')
plt.plot(bin_centers, fit_function(result.x, bin_centers), label='Fit')
plt.xlabel('Invariant Mass (GeV)')
plt.ylabel('Counts')
plt.legend()
plt.savefig('figures/zfit.svg')

# Write REPORT.md
report_content = """
## Selection
We applied the following quality selection criteria to the data:
- Both muons are global (`type1 == 'G'` and `type2 == 'G'`)
- |eta1| < 2.4 and |eta2| < 2.4
- pt1 > 20 GeV and pt2 > 20 GeV
- Opposite charge (Q1 * Q2 < 0)

## Fit
The invariant mass distribution was fitted with a combination of a Gaussian for signal and an exponential for background. The fit parameters are as follows:
- mu (Z boson mass): {mu_fit:.2f} GeV ± {mu_err:.2f} GeV
- sigma (Z boson width): {sigma_fit:.2f} GeV ± {sigma_err:.2f} GeV
- n_sig (number of signal events): {n_sig_fit:.2f} ± {n_sig_err:.2f}
- n_bkg (number of background events): {n_bkg_fit:.2f}
- lam (background decay constant): {lam_fit:.2f} GeV
- Chi-squared: {chi2:.2f} with {ndf} degrees of freedom

## Assumptions
- The invariant mass M is already computed in the dataset.
- The fit function uses a Gaussian for signal and an exponential for background.
- The initial parameter values were chosen based on typical Z boson properties.
- The covariance matrix was used to estimate uncertainties.
""".format(**result_json)

with open('REPORT.md', 'w') as f:
    f.write(report_content)
"""

with open('solution.py', 'w') as f:
    f.write(solution_code)

# Submit the task
submit(summary="Completed end-to-end Z-peak measurement analysis on CMS open data. Key results: n_total={len(data)}, n_selected={len(selected)}, efficiency={len(selected)/len(data):.2f}, n_in_range={len(in_range)}, fit parameters as per result.json.", assumptions=["Assumed the invariant mass M is already computed in the dataset.", "Used a Gaussian for signal and an exponential for background in the fit function.", "Initial parameter values were based on typical Z boson properties.", "Covariance matrix was used to estimate uncertainties."])