Invariant-mass bump

Python Code

"""Particle Physics — Invariant-mass bump (dimuon spectrum)."""


import numpy as np
import plotly.graph_objects as go


def generate():
    print("generating synthetic dimuon invariant-mass spectrum …")
    rng = np.random.default_rng(42)

    # Background: exponential decay in log space from 0.3 to 120 GeV
    n_bg = 100_000
    # Exponential in linear mass; rate such that spectrum falls steeply
    # Use inverse-CDF: M = 0.3 * exp(u * ln(120/0.3)) where u ~ Uniform(0,1)
    # weighted by exponential: P(M) ~ M^{-2.7} (approximate QCD background)
    # Simple approach: sample from power law P(M) ~ M^{-n}
    # CDF: F(M) = (M^{1-n} - M_min^{1-n}) / (M_max^{1-n} - M_min^{1-n})
    # For n=2.7: 1-n = -1.7
    alpha = -1.7
    M_min, M_max = 0.3, 120.0
    u = rng.uniform(0, 1, n_bg)
    bg = (u * (M_max**alpha - M_min**alpha) + M_min**alpha) ** (1.0 / alpha)

    # J/psi peak at 3.097 GeV
    jpsi = rng.normal(3.097, 0.05, 1000)
    jpsi = jpsi[(jpsi > 2.5) & (jpsi < 4.0)]

    # Upsilon(1S) peak at 9.46 GeV
    upsilon = rng.normal(9.46, 0.15, 500)
    upsilon = upsilon[(upsilon > 8.5) & (upsilon < 11.0)]

    # Z boson peak at 91.2 GeV
    z = rng.normal(91.2, 2.5, 200)
    z = z[(z > 80) & (z < 110)]

    M = np.concatenate([bg, jpsi, upsilon, z])
    M = M[(M >= 0.3) & (M <= 120)]

    print(f"  total events: {len(M):,}")

    fig = go.Figure()
    fig.add_trace(go.Histogram(
        x=M,
        nbinsx=200,
        histnorm='',
        marker_color=VIOLET,
        opacity=0.85,
        name="Events",
        hovertemplate="M ≈ %{x:.2f} GeV<br>Events: %{y}<extra></extra>",
    ))

    # Annotations for peaks
    annotations = [
        dict(
            x=np.log10(3.097),
            y=np.log10(220),
            xref="x", yref="y",
            text="J/ψ<br>3.10 GeV",
            showarrow=True,
            arrowhead=2,
            ax=40, ay=-50,
            font=dict(size=12, color=TEAL),
            arrowcolor=TEAL,
        ),
        dict(
            x=np.log10(9.46),
            y=np.log10(60),
            xref="x", yref="y",
            text="Υ(1S)<br>9.46 GeV",
            showarrow=True,
            arrowhead=2,
            ax=45, ay=-45,
            font=dict(size=12, color=GREEN),
            arrowcolor=GREEN,
        ),
        dict(
            x=np.log10(91.2),
            y=np.log10(25),
            xref="x", yref="y",
            text="Z<br>91.2 GeV",
            showarrow=True,
            arrowhead=2,
            ax=5, ay=-55,
            font=dict(size=12, color=PINK),
            arrowcolor=PINK,
        ),
    ]

    fig.update_layout(
        xaxis=dict(
            title="Invariant Mass M (GeV)",
            type="log",
            tickvals=[0.5, 1, 2, 3, 5, 10, 20, 50, 100],
            ticktext=["0.5", "1", "2", "3", "5", "10", "20", "50", "100"],
        ),
        yaxis=dict(
            title="Events",
            type="log",
        ),
        annotations=annotations,
        bargap=0,
        margin=dict(t=50, b=60, l=70, r=60),
        showlegend=False,
    )

    return fig


if __name__ == "__main__":
    generate()

Made with Plotly