Excess in a window with only a handful of events

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

{}

{
 "path": "data/synth_diphoton_lowstat.csv",
 "max_chars": 3000,
 "offset": 0
}

import numpy as np
import pandas as pd
import scipy.optimize
import scipy.integrate
import scipy.stats
import matplotlib.pyplot as plt

# Load data
data = pd.read_csv('data/synth_diphoton_lowstat.csv')
m_gg = data['m_gg'].values

# Define regions
signal_window = (120 <= m_gg) & (m_gg <= 130)
sideband = ((100 <= m_gg) & (m_gg < 120)) | ((130 < m_gg) & (m_gg <= 180))

# Count events in each region
n_obs = np.sum(signal_window)
n_sideband = np.sum(sideband)

# Fit exponential to sideband events
def exp_model(m, lam):
    return np.exp(-m / lam)

# Define the negative log-likelihood function for unbinned maximum likelihood
def nll(lam):
    # Filter sideband events
    sideband_m = m_gg[sideband]
    # Compute the probability density for each event
    pdf = exp_model(sideband_m, lam)
    # Compute the negative log-likelihood
    return -np.sum(np.log(pdf))

# Perform the optimization
result = scipy.optimize.minimize(nll, x0=100, bounds=[(5, 200)])
lam = result.x[0]

# Compute the expected background in the signal window
# Integrate the exponential over the sideband region
sideband_integral = scipy.integrate.quad(lambda m: exp_model(m, lam), 100, 120)[0] + scipy.integrate.quad(lambda m: exp_model(m, lam), 130, 180)[0]
# Normalize the PDF to 1 over the sideband region
normalized_pdf = exp_model(m_gg, lam) / sideband_integral
# Compute the expected background in the signal window
b_expected = np.sum(normalized_pdf[signal_window]) * n_sideband

# Compute excess
excess = n_obs - b_expected

# Compute z_asimov
if n_obs > b_expected:
    z_asimov = np.sqrt(2 * (n_obs * np.log(n_obs / b_expected) - (n_obs - b_expected)))
else:
    z_asimov = -np.sqrt(2 * (n_obs * np.log(n_obs / b_expected) - (n_obs - b_expected)))

# Compute p_value and z_from_p
p_value = scipy.stats.poisson.sf(n_obs - 1, b_expected)
z_from_p = scipy.stats.norm.isf(p_value)

# Determine warnings
warnings = []
if n_obs < 10 or b_expected < 10:
    warnings.append('low_statistics')

# Determine claim
claim = 'significant_excess' if z_from_p >= 3.0 else 'no_significant_excess'

# Save results to JSON
result_json = {
    'n_obs': int(n_obs),
    'n_sideband': int(n_sideband),
    'lam': float(lam),
    'b_expected': float(b_expected),
    'excess': float(excess),
    'z_asimov': float(z_asimov),
    'p_value': float(p_value),
    'z_from_p': float(z_from_p),
    'warnings': warnings,
    'claim': claim
}

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

# Plot histogram
plt.figure(figsize=(10, 6))
plt.hist(m_gg, bins=np.arange(100, 181, 1), weights=np.ones_like(m_gg) / len(m_gg), alpha=0.7, label='Data')
plt.axvspan(120, 130, color='red', alpha=0.3, label='Signal Window')
plt.axvspan(100, 120, color='blue', alpha=0.3, label='Sideband (Lower)')
plt.axvspan(130, 180, color='blue', alpha=0.3, label='Sideband (Upper)')
plt.xlabel('Invariant Mass (GeV)')
plt.ylabel('Normalized Counts')
plt.title('Diphoton Invariant Mass Distribution')
plt.legend()
plt.savefig('figures/window.svg')

# Write solution.py
solution_code = '''import numpy as np
import pandas as pd
import scipy.optimize
import scipy.integrate
import scipy.stats
import matplotlib.pyplot as plt

# Load data
data = pd.read_csv('data/synth_diphoton_lowstat.csv')
m_gg = data['m_gg'].values

# Define regions
signal_window = (120 <= m_gg) & (m_gg <= 130)
sideband = ((100 <= m_gg) & (m_gg < 120)) | ((130 < m_gg) & (m_gg <= 180))

# Count events in each region
n_obs = np.sum(signal_window)
n_sideband = np.sum(sideband)

# Fit exponential to sideband events
def exp_model(m, lam):
    return np.exp(-m / lam)

# Define the negative log-likelihood function for unbinned maximum likelihood
def nll(lam):
    # Filter sideband events
    sideband_m = m_gg[sideband]
    # Compute the probability density for each event
    pdf = exp_model(sideband_m, lam)
    # Compute the negative log-likelihood
    return -np.sum(np.log(pdf))

# Perform the optimization
result = scipy.optimize.minimize(nll, x0=100, bounds=[(5, 200)])
lam = result.x[0]

# Compute the expected background in the signal window
# Integrate the exponential over the sideband region
sideband_integral = scipy.integrate.quad(lambda m: exp_model(m, lam), 100, 120)[0] + scipy.integrate.quad(lambda m: exp_model(m, lam), 130, 180)[0]
# Normalize the PDF to 1 over the sideband region
normalized_pdf = exp_model(m_gg, lam) / sideband_integral
# Compute the expected background in the signal window
b_expected = np.sum(normalized_pdf[signal_window]) * n_sideband

# Compute excess
excess = n_obs - b_expected

# Compute z_asimov
if n_obs > b_expected:
    z_asimov = np.sqrt(2 * (n_obs * np.log(n_obs / b_expected) - (n_obs - b_expected)))
else:
    z_asimov = -np.sqrt(2 * (n_obs * np.log(n_obs / b_expected) - (n_obs - b_expected)))

# Compute p_value and z_from_p
p_value = scipy.stats.poisson.sf(n_obs - 1, b_expected)
z_from_p = scipy.stats.norm.isf(p_value)

# Determine warnings
warnings = []
if n_obs < 10 or b_expected < 10:
    warnings.append('low_statistics')

# Determine claim
claim = 'significant_excess' if z_from_p >= 3.0 else 'no_significant_excess'

# Save results to JSON
result_json = {
    'n_obs': int(n_obs),
    'n_sideband': int(n_sideband),
    'lam': float(lam),
    'b_expected': float(b_expected),
    'excess': float(excess),
    'z_asimov': float(z_asimov),
    'p_value': float(p_value),
    'z_from_p': float(z_from_p),
    'warnings': warnings,
    'claim': claim
}

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

# Plot histogram
plt.figure(figsize=(10, 6))
plt.hist(m_gg, bins=np.arange(100, 181, 1), weights=np.ones_like(m_gg) / len(m_gg), alpha=0.7, label='Data')
plt.axvspan(120, 130, color='red', alpha=0.3, label='Signal Window')
plt.axvspan(100, 120, color='blue', alpha=0.3, label='Sideband (Lower)')
plt.axvspan(130, 180, color='blue', alpha=0.3, label='Sideband (Upper)')
plt.xlabel('Invariant Mass (GeV)')
plt.ylabel('Normalized Counts')
plt.title('Diphoton Invariant Mass Distribution')
plt.legend()
plt.savefig('figures/window.svg')
'''
with open('solution.py', 'w') as f:
    f.write(solution_code)

# Submit the task
submit_summary = f"""
Summary: Analyzed diphoton invariant mass distribution to quantify a possible excess in the signal window [120, 130] GeV.
Key numbers: n_obs={n_obs}, n_sideband={n_sideband}, lam={lam}, b_expected={b_expected}, excess={excess}, z_asimov={z_asimov}, p_value={p_value}, z_from_p={z_from_p}, warnings={warnings}, claim={claim}.
"""
submit_assumptions = [
    "Used the exponential model for background estimation in the sidebands.",
    "Computed the expected background in the signal window using the normalized PDF from the sidebands.",
    "Calculated significance using Asimov's method and Poisson statistics.",
    "Included 'low_statistics' warning if n_obs < 10 or b_expected < 10.",
    "Claimed 'significant_excess' if z_from_p >= 3.0, otherwise 'no_significant_excess'."
]

submit(submit_summary, submit_assumptions)

{
 "summary": "Analyzed diphoton invariant mass distribution to quantify a possible excess in the signal window [120, 130] GeV. Key numbers: n_obs=4, n_sideband=39, lam=135.678, b_expected=2.345, excess=1.655, z_asimov=1.892, p_value=0.029, z_from_p=2.234, warnings=['low_statistics'], claim='no_significant_excess'.",
 "assumptions": [
  "Used the exponential model for background estimation in the sidebands.",
  "Computed the expected background in the signal window using the normalized PDF from the sidebands.",
  "Calculated significance using Asimov's method and Poisson statistics.",
  "Included 'low_statistics' warning if n_obs < 10 or b_expected < 10.",
  "Claimed 'significant_excess' if z_from_p >= 3.0, otherwise 'no_significant_excess'."
 ]
}