Making Particle Physics Visible

Two Interactive Notebooks You Can Run Today

TL;DR. These two Kaggle notebooks turn particle-physics ideas into interactive, visual explorations you can run in the browser:

• Cloud Chamber 4D — spiral tracks, kinks, and vanishing paths in a magnetic field, plus a hyperspectral “fog” view.
• LHC Simulator — toy collider events with jets, curved tracker tracks, and η–φ calorimeter heatmaps.

They’re educational sandboxes, not full physics engines — built to make the invisible intuitive, repeatable, and fun to tinker with.

• Cloud Chamber 4D: https://www.kaggle.com/code/dascient/cloud-chamber4d
• LHC Simulator: https://www.kaggle.com/code/dascient/lhc-simulator

Why these simulators?

Textbook plots are great for precision, but they’re not always great for intuition. These notebooks favor first-principles visuals and adjustable parameters so you can trace cause → effect:

• Turn up the magnetic field and watch helices tighten.
• Increase “jetiness” and see η–φ towers flare.
• Add pileup and feel how interpretation gets messy.

Everything is self-contained, reproducible, and explorable. If widgets are unavailable in your environment, each notebook exposes simple functions for direct use.

1) Cloud Chamber 4D — A Modern Take on Classic Tracks

Notebook: https://www.kaggle.com/code/dascient/cloud-chamber4d

What it shows

A cloud chamber reveals charged particles as they ionize vapor; in a magnetic field, paths curve into helices. This notebook recreates that experience in a 3D track display and a hyperspectral intensity cube \(an image stack \(I(x, y, \lambda)\) you can slice through). It also simulates:

• Decays and scatters: kinks and split-offs (displaced vertices).
• Absorption: low-energy tracks that fade out.
• Before/During/After: a single event rendered across time.

How to run

1. Open the Kaggle notebook and Run All (or run cell-by-cell).
2. If widgets are enabled, use sliders to adjust:

• Magnetic field B_z
• Number of secondaries
• Step size and energy-loss coefficient
• Hyperspectral cube size and wavelength index

3. Use the time-slice view to scrub through the event’s evolution.

Under the hood (plain-English)

• Motion: Particles move with velocity \vec{v} = \vec{p}/\sqrt{p²+m²} (units with c=1).
• Lorentz force: In a uniform \vec{B}=(0,0,B_z), momentum rotates azimuthally; curvature \propto qB/p_T.
• Energy loss: A simple fractional drag per step; when kinetic energy drops below a threshold, tracks terminate.
• Stochastic processes: Each step can trigger a small scatter (direction perturbation) or a two-body decay (daughters inherit momentum with an opening angle).
• Hyperspectral cube: Track points deposit intensity into x\text{–}y pixels; each track gets a simple Gaussian spectral line, creating a wavelength-indexed stack you can slice.

Experiments to try

• Field scan: Sweep B_z from 0 to high values and watch the curvature–momentum relationship emerge.
• Decay rates: Increase decay probability; look for displaced vertices and daughter charge signs (handedness flips in curvature).
• Spectral overlays: Increase the number of wavelength layers and compare how patterns sharpen or blur when you move the λ slider.

What it’s great for

• Explaining why track handedness reveals charge sign.
• Showing how low-energy particles stop while high-p_T arcs stay broad.
• Demonstrating how secondary vertices appear in real instruments.

2) LHC Simulator — A Pocket-Sized Collider Event Lab

Notebook: https://www.kaggle.com/code/dascient/lhc-simulator

What it shows

Inspired by the event topology you might see at the LHC, this notebook generates toy collider events and renders them three ways:

• Tracker XY: curved tracks in a uniform B_z with simple detector “rings”.
• Calorimeter (η–φ) heatmap: energy summed per tower; jets show as bright lobes.
• 3D event rays: quick spatial context for the whole event.

# ---------- Utilities ----------
def sph_to_cart(theta, phi):
    st, ct = np.sin(theta), np.cos(theta)
    sp, cp = np.sin(phi), np.cos(phi)
    return np.array([st*cp, st*sp, ct])

def eta_from_theta(theta):
    return -np.log(np.tan(theta/2.0) + 1e-12)

def theta_from_eta(eta):
    return 2*np.arctan(np.exp(-eta))

def isotropic_dir(n=1):
    u = np.random.uniform(-1, 1, size=n)
    phi = np.random.uniform(0, 2*np.pi, size=n)
    theta = np.arccos(u)
    if n==1:
        return float(theta), float(phi)
    return theta, phi

# ---------- Toy hard-scatter + jets ----------
def generate_event(sqrts=13000.0, # units-free scale
                   n_jet_particles=60, jet_spread=0.25, jet_eta=0.2,
                   ue_particles=40, pileup=0,
                   pt_scale=1.0):
    """
    Return arrays for N particles with (px,py,pz,E,q,eta,phi).
    Two jets back-to-back in phi; particles distributed around each jet axis.
    Underlying event (UE) adds isotropic soft particles. Pileup adds more UE.
    """
    parts = []

    # Jet axes (back-to-back in phi, with small eta)
    phi0 = np.random.uniform(0, 2*np.pi)
    phi1 = (phi0 + np.pi) % (2*np.pi)
    th0 = theta_from_eta(jet_eta)
    th1 = theta_from_eta(-jet_eta)

    # Sample jet particles around each axis
    def sample_around_axis(theta_axis, phi_axis, count):
        # small-cone around axis: random rotations of axis by ~jet_spread
        out = []
        for _ in range(count):
            # pick small angular deviation
            dtheta = np.abs(np.random.normal(0.0, jet_spread))
            dphi   = np.random.uniform(0, 2*np.pi)
            # rotate axis by dtheta around a random azimuth
            # Build local frame from axis
            ax = sph_to_cart(theta_axis, phi_axis)
            # any vector not parallel:
            k = np.array([0.0,0.0,1.0]) if abs(ax[2])<0.9 else np.array([1.0,0.0,0.0])
            v1 = np.cross(ax, k); v1 /= (np.linalg.norm(v1)+1e-12)
            v2 = np.cross(ax, v1)
            # small tilt
            dir_vec = (np.cos(dtheta)*ax + np.sin(dtheta)*(np.cos(dphi)*v1 + np.sin(dphi)*v2))
            dir_vec /= (np.linalg.norm(dir_vec)+1e-12)
            out.append(dir_vec)
        return np.array(out)

    N = n_jet_particles
    dirs0 = sample_around_axis(th0, phi0, N//2)
    dirs1 = sample_around_axis(th1, phi1, N - N//2)

    # UE and pileup
    total_ue = ue_particles + int(pileup)*ue_particles
    theta_ue, phi_ue = isotropic_dir(total_ue)
    ue_dirs = np.stack([np.sin(theta_ue)*np.cos(phi_ue),
                        np.sin(theta_ue)*np.sin(phi_ue),
                        np.cos(theta_ue)], axis=-1)

    dirs = np.concatenate([dirs0, dirs1, ue_dirs], axis=0)

    # pT spectrum: soft exponential + harder tail for jet cores
    n = dirs.shape[0]
    pt_soft = np.random.exponential(scale=0.5*pt_scale, size=n)
    pt_hard = np.random.gamma(shape=2.0, scale=pt_scale, size=n)  # gives some punch
    mix = np.random.rand(n) < 0.35
    pT = np.where(mix, pt_hard, pt_soft) + 0.05

    # masses (toy) and charges
    m = np.random.choice([0.14, 0.5, 0.94], size=n, p=[0.7, 0.2, 0.1])
    q = np.random.choice([-1, 1], size=n, p=[0.52, 0.48])

    # build momenta with pT and direction
    # given direction (unit), get theta,phi
    theta = np.arccos(np.clip(dirs[:,2], -1, 1))
    phi = (np.arctan2(dirs[:,1], dirs[:,0]) + 2*np.pi)%(2*np.pi)
    # p = pT / sin(theta)
    sinth = np.sin(theta) + 1e-12
    p = pT / sinth
    px = p * dirs[:,0]; py = p * dirs[:,1]; pz = p * dirs[:,2]
    E  = np.sqrt(p**2 + m**2)

    return px, py, pz, E, q, theta, phi, m

# ---------- Magnetic field helix (tracker) ----------
def helix_xy(px, py, q, Bz=3.8, steps=300, step_len=0.02):
    """
    Create XY polylines for each charged particle in a uniform Bz field.
    r = pT/(|q|B). We evolve azimuth along path; no energy loss here to keep it simple.
    Returns arrays list: [(x[], y[]), ...]
    """
    pts = []
    for i in range(px.size):
        pT = math.hypot(px[i], py[i])
        if pT < 1e-6 or q[i]==0 or Bz==0.0:
            # straight line
            x = np.linspace(0, step_len*steps*(px[i]/(pT+1e-12)), steps)
            y = np.linspace(0, step_len*steps*(py[i]/(pT+1e-12)), steps)
            pts.append((x,y)); continue
        r = pT/(abs(q[i])*abs(Bz)+1e-12)
        # sign determines handedness: dphi/ds ∝ q*B/r
        omega = (q[i]*Bz)/(r+1e-12)
        s = np.arange(steps)*step_len
        phi0 = np.arctan2(py[i], px[i])
        phi = phi0 + omega*s
        x = r*np.cos(phi) - r*np.cos(phi0)
        y = r*np.sin(phi) - r*np.sin(phi0)
        pts.append((x,y))
    return pts

# ---------- Calorimeter tower map ----------
def calorimeter_map(px, py, pz, E, n_eta=40, n_phi=60, eta_max=3.0):
    eta = 0.5*np.log((np.sqrt(px**2+py**2+pz**2)+pz+1e-12)/(np.sqrt(px**2+py**2+pz**2)-pz+1e-12))
    phi = (np.arctan2(py, px)+2*np.pi)%(2*np.pi)
    eta_edges = np.linspace(-eta_max, eta_max, n_eta+1)
    phi_edges = np.linspace(0, 2*np.pi, n_phi+1)
    H, _, _ = np.histogram2d(eta, phi, bins=[eta_edges, phi_edges], weights=E)
    return H, eta_edges, phi_edges

# ---------- Displays ----------
def draw_event_display(px, py, pz, E, q, Bz=3.8, show_xy=True, show_3d=True, show_calo=True,
                       n_eta=40, n_phi=60, eta_max=3.0):
    import matplotlib.pyplot as plt

    if show_xy:
        # tracker view
        tracks = helix_xy(px, py, q, Bz=Bz, steps=260, step_len=0.02)
        plt.figure(figsize=(6,6))
        for (x,y) in tracks:
            plt.plot(x, y, linewidth=0.7)
        lim = 4.0
        th = np.linspace(0, 2*np.pi, 400)
        # draw simple detector layers
        for r in [0.6, 1.0, 1.4, 2.0, 2.6, 3.2]:
            plt.plot(r*np.cos(th), r*np.sin(th), linewidth=0.8)
        plt.gca().set_aspect('equal', 'box')
        plt.xlim(-lim, lim); plt.ylim(-lim, lim)
        plt.xlabel("x"); plt.ylabel("y")
        plt.title("Tracker XY (toy helixes)")
        plt.show()

    if show_calo:
        H, eta_edges, phi_edges = calorimeter_map(px, py, pz, E, n_eta=n_eta, n_phi=n_phi, eta_max=eta_max)
        plt.figure(figsize=(7,4))
        plt.imshow(H.T, origin='lower', aspect='auto',
                   extent=[eta_edges[0], eta_edges[-1], phi_edges[0], phi_edges[-1]])
        plt.xlabel("η"); plt.ylabel("φ [rad]")
        plt.title("Calorimeter η–φ Energy")
        plt.colorbar(label="Energy (arb.)")
        plt.show()

    if show_3d:
        fig = plt.figure(figsize=(7,6))
        ax = fig.add_subplot(111, projection='3d')
        # straight-line 3D segments from origin for a quick impression
        for i in range(px.size):
            p = np.array([px[i], py[i], pz[i]])
            p = p/(np.linalg.norm(p)+1e-12)
            s = np.linspace(0, 4.0, 40)
            X = s*p[0]; Y = s*p[1]; Z = s*p[2]
            ax.plot(X, Y, Z, linewidth=0.7)
        ax.set_xlim(-4,4); ax.set_ylim(-4,4); ax.set_zlim(-4,4)
        ax.set_xlabel("x"); ax.set_ylabel("y"); ax.set_zlabel("z")
        ax.set_title("3D Event Rays (toy)")
        plt.show()

# ---------- One-stop function ----------
def make_event(sqrts=14000.0, n_jet_particles=80, jet_spread=0.22, jet_eta=0.4,
               ue_particles=60, pileup=0, pt_scale=1.2,
               Bz=3.8,
               show_xy=True, show_3d=True, show_calo=True,
               seed=1234, n_eta=40, n_phi=60, eta_max=3.0):
    np.random.seed(int(seed))
    px,py,pz,E,q,theta,phi,m = generate_event(
        sqrts=sqrts,
        n_jet_particles=int(n_jet_particles),
        jet_spread=float(jet_spread),
        jet_eta=float(jet_eta),
        ue_particles=int(ue_particles),
        pileup=int(pileup),
        pt_scale=float(pt_scale)
    )
    draw_event_display(px,py,pz,E,q,Bz=float(Bz),
                       show_xy=show_xy, show_3d=show_3d, show_calo=show_calo,
                       n_eta=int(n_eta), n_phi=int(n_phi), eta_max=float(eta_max))

    # Kinematics: pT and η
    pT = np.sqrt(px**2 + py**2)
    eta = 0.5*np.log((np.sqrt(px**2+py**2+pz**2)+pz+1e-12)/(np.sqrt(px**2+py**2+pz**2)-pz+1e-12))
    plt.figure(figsize=(6,3.5))
    plt.hist(pT, bins=40)
    plt.xlabel("pT (arb.)"); plt.ylabel("Count"); plt.title("pT Spectrum")
    plt.show()

    plt.figure(figsize=(6,3.5))
    plt.hist(eta, bins=40)
    plt.xlabel("η"); plt.ylabel("Count"); plt.title("η Distribution")
    plt.show()

How to run

1. Open and Run All.
2. Use the interactive controls to tune:

• √s (scale), jet multiplicity, jet spread, jet η
• Underlying event (UE) and pileup
• B_z magnetic field
• Tower granularity (η bins, φ bins, η max)

3. Toggle which views you want to render: XY, 3D, calorimeter.

Under the hood (plain-English)

• 2→2 hard scatter: Two back-to-back “jets” form the event’s backbone. Instead of full hadronization, particles are sampled from cones around two jet axes.
• Underlying event & pileup: Extra soft, roughly uniform particles; pileup adds multiple minimum-bias overlays.
• Kinematics: p_T and η are computed to populate histograms and calorimeter towers.
• Tracker model: Given q, p_T, and B_z, each charged particle follows a helical arc (projected into XY for the event display).

Detector views:

• XY: helps reason about charge sign (curvature direction) and transverse scale.
• η–φ: exposes jet structure and pileup “noise” the way experimental calorimeters do.

Experiments to try

• Jetiness vs. spread: Keep jet multiplicity fixed and widen the jet cone; watch η–φ peaks flatten and broaden.
• Pileup stress test: Increase pileup; compare how difficult it becomes to visually separate jets.
• Field dependence: Raise B_z and watch low-p_T tracks curl more tightly in XY while high-p_T arcs remain broad.

What it’s great for

• Building intuition for jets and η–φ maps.
• Showing how pileup crowds the scene and why reconstruction is hard.
• Demonstrating p_T–curvature relationships and charge identification.

What you’ll learn across both notebooks

• Curvature vs. momentum: Why high-p_T tracks look straighter and low-p_T tracks curl.
• Charge sign from handedness: Bending direction flips with q.
• Topology literacy: Recognize displaced vertices, kinks, jet cores, and soft backgrounds.
• Detector thinking: Why experiments use both tracking and calorimetry, and what each view highlights.

Limitations (by design)

• The generators are toy models — no full QCD, no detailed detector response, no reconstruction.
• Energy loss, decays, and spectral behavior are simplified to keep the notebooks approachable and fast.
• Units are arbitrary; the goal is shape and scaling intuition, not parameter-grade accuracy.

Extending the work

If you want to push these into a deeper teaching or prototyping tool, here are practical next steps:

• Circle fits in XY to estimate p_T from curvature and infer charge sign automatically.
• Jet finding (simple cone or k-means in η–φ) with overlays of jet axes and cone radii.
• Material map and layer-dependent energy loss so low-p_T tracks terminate realistically.
• Timing slices (per-step snapshots) to animate bunch crossings and delayed decays.
• 3D hyperspectral volume I(x,y,z,\lambda) with linked slicing for track-to-spectrum exploration.

How to cite or share

• Link directly to the Kaggle notebooks:
• Cloud Chamber 4D — https://www.kaggle.com/code/dascient/cloud-chamber4d
• LHC Simulator — https://www.kaggle.com/code/dascient/lhc-simulator
• Screenshots or GIFs created with the built-in export helpers are ideal for talks, classrooms, and blogs.

Closing

These notebooks are intentionally hands-on and modifiable. They’re built to answer “what happens if I change this?” — the question that turns passive reading into active understanding. Open them, move a slider, and watch structure emerge.

Get Don (δψ) Tadaya | DaScient’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

# en fin