.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/particle/resonance_discovery/plot_04_higgs_diphoton_discovery.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_particle_resonance_discovery_plot_04_higgs_diphoton_discovery.py: ATLAS and CMS: discovery of the Higgs boson ============================================== On 4 July 2012 ATLAS and CMS each reported a new boson near 125 GeV at five standard deviations. The cleanest channel was :math:`H\to\gamma\gamma`: a rare decay (0.2%) but one where both photons are measured precisely, giving a narrow peak in the diphoton invariant mass. It sits on a smooth, steeply falling background of ordinary photon pairs more than ten times larger under the peak. Finding it is a question of statistics. This example builds signal photon pairs with :func:`~physicskit.particle.decays.two_body_decay`, :func:`~physicskit.particle.kinematics.boost_generic` and :func:`~physicskit.particle.kinematics.invariant_mass`, adds an exponential background, fits the background shape from the sidebands away from 125 GeV, and fits the signal strength with a binned Poisson likelihood. The local significance :math:`Z=\sqrt{2\,\Delta\ln L}` grows as the square root of the collected data and crosses 5 sigma. .. GENERATED FROM PYTHON SOURCE LINES 23-90 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy.optimize import minimize, minimize_scalar from physicskit.particle.decays import two_body_decay from physicskit.particle.kinematics import FourVector, boost_generic, invariant_mass M_H = 125.1 rng = np.random.default_rng(2012) bins = np.linspace(105, 160, 111) centres = 0.5 * (bins[1:] + bins[:-1]) def diphoton_signal(n): """Reconstructed m_gammagamma for n Higgs decays, 1.2% photon energy resolution.""" out = np.empty(n) for k in range(n): g1, g2 = two_body_decay(M_H, 0.0, 0.0, rng.uniform(-1, 1), rng.uniform(0, 2 * np.pi)) beta = rng.normal(0, 0.35, 3) beta *= min(1.0, 0.9 / np.linalg.norm(beta)) photons = [] for g in (boost_generic(g1, beta), boost_generic(g2, beta)): E = g.E * (1 + 0.012 * rng.standard_normal()) photons.append(FourVector(E, *(E * g.p_vec / g.p_mag))) out[k] = invariant_mass(photons) return out def background(n, slope=0.035): u = rng.uniform(size=n) lo, hi = bins[0], bins[-1] return lo - np.log(1 - u * (1 - np.exp(-slope * (hi - lo)))) / slope signal_pool = diphoton_signal(20000) signal_shape = np.histogram(signal_pool, bins=bins)[0].astype(float) signal_shape /= signal_shape.sum() print(f"signal peak width (std of reconstructed mass): {np.std(signal_pool):.2f} GeV") sideband = np.abs(centres - M_H) > 5 def poisson_nll(counts, lam): lam = np.clip(lam, 1e-9, None) return np.sum(lam - counts * np.log(lam)) def fit_background(counts): """Exponential fit to the sidebands (away from 125 GeV); returns b(m) in every bin.""" x = centres - bins[0] res = minimize( lambda q: poisson_nll(counts[sideband], np.exp(q[0] - q[1] * x[sideband])), [np.log(counts[0] + 1), 0.03], method="Nelder-Mead", ) return np.exp(res.x[0] - res.x[1] * x) def significance(counts): """Signal strength and local significance, background fixed from the sidebands.""" b = fit_background(counts) res = minimize_scalar(lambda mu: poisson_nll(counts, b + mu * signal_shape), bounds=(0, 5 * counts.sum() ** 0.5 * 10), method="bounded") Z = np.sqrt(max(2 * (poisson_nll(counts, b) - res.fun), 0.0)) return Z, res.x, b .. rst-class:: sphx-glr-script-out .. code-block:: none signal peak width (std of reconstructed mass): 1.07 GeV .. GENERATED FROM PYTHON SOURCE LINES 91-97 Data sets of increasing size -------------------------------- One unit of luminosity holds 200 signal and 40,000 background events. Signal and background both grow in proportion to it. Any one data set fluctuates, so each luminosity is repeated in 25 pseudo-experiments; the median significance grows as :math:`\sqrt{L}`. .. GENERATED FROM PYTHON SOURCE LINES 97-132 .. code-block:: Python lumis = [0.5, 1, 2, 3, 4, 6] Z_median, Z_lo, Z_hi = [], [], [] for L in lumis: Zs = [] for _ in range(25): n_s = rng.poisson(200 * L) data = np.concatenate([rng.choice(signal_pool, n_s), background(rng.poisson(40000 * L))]) counts = np.histogram(data, bins=bins)[0].astype(float) Zs.append(significance(counts)[0]) Z_median.append(np.median(Zs)) Z_lo.append(np.percentile(Zs, 16)) Z_hi.append(np.percentile(Zs, 84)) print(f"luminosity {L:>3}: median local significance {Z_median[-1]:.1f} sigma (68% of experiments {Z_lo[-1]:.1f}-{Z_hi[-1]:.1f})") fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4)) counts = np.histogram(np.concatenate([rng.choice(signal_pool, 1200), background(240000)]), bins=bins)[0].astype(float) ax1.errorbar(centres, counts, yerr=np.sqrt(counts), fmt="o", color="k", ms=2, lw=0.8, label="simulated data") ax1.plot(centres, significance(counts)[2], color="steelblue", label="background fit (sidebands)") ax1.axvline(M_H, color="firebrick", ls=":", label=f"$m_H$ = {M_H} GeV") ax1.set_xlabel(r"$m_{\gamma\gamma}$ [GeV]") ax1.set_ylabel("events / 0.5 GeV") ax1.set_title(r"$H\to\gamma\gamma$: a small bump on a large background") ax1.legend(fontsize=8) ax2.fill_between(lumis, Z_lo, Z_hi, color="firebrick", alpha=0.2, label="68% of pseudo-experiments") ax2.plot(lumis, Z_median, "o-", color="firebrick", label="median significance") ax2.plot(lumis, Z_median[2] * np.sqrt(np.array(lumis) / lumis[2]), "--", color="0.5", label=r"$\propto\sqrt{L}$") ax2.axhline(5, color="k", ls=":", label=r"5$\sigma$ discovery threshold") ax2.set_xlabel("integrated luminosity (arbitrary units)") ax2.set_ylabel(r"local significance [$\sigma$]") ax2.set_title("Collecting data until the bump is a discovery") ax2.legend(fontsize=8) fig.tight_layout() plt.show() .. image-sg:: /api/gallery/particle/resonance_discovery/images/sphx_glr_plot_04_higgs_diphoton_discovery_001.png :alt: $H\to\gamma\gamma$: a small bump on a large background, Collecting data until the bump is a discovery :srcset: /api/gallery/particle/resonance_discovery/images/sphx_glr_plot_04_higgs_diphoton_discovery_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none luminosity 0.5: median local significance 2.4 sigma (68% of experiments 1.2-3.9) luminosity 1: median local significance 4.3 sigma (68% of experiments 3.2-5.0) luminosity 2: median local significance 4.9 sigma (68% of experiments 4.2-6.1) luminosity 3: median local significance 6.3 sigma (68% of experiments 5.1-7.4) luminosity 4: median local significance 6.9 sigma (68% of experiments 6.4-7.5) luminosity 6: median local significance 8.5 sigma (68% of experiments 7.3-9.9) .. GENERATED FROM PYTHON SOURCE LINES 133-137 Check ----- The peak width is the photon resolution, 1.2% x M_H / sqrt(2); the median significance grows as sqrt(L) and passes 5 sigma within the scan. .. GENERATED FROM PYTHON SOURCE LINES 137-140 .. code-block:: Python assert abs(np.std(signal_pool) / (0.012 * M_H / np.sqrt(2)) - 1) < 0.15 assert abs(Z_median[-1] / (Z_median[2] * np.sqrt(lumis[-1] / lumis[2])) - 1) < 0.25 assert Z_median[0] < 5 < Z_median[-1] .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.993 seconds) .. _sphx_glr_download_api_gallery_particle_resonance_discovery_plot_04_higgs_diphoton_discovery.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_04_higgs_diphoton_discovery.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_04_higgs_diphoton_discovery.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_04_higgs_diphoton_discovery.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_