Building a mapping for a Feynman diagram

The previous examples built a mapping from a flat list of masses. This example instead follows one Feynman diagram, using Diagram and Topology. The diagram is a single-quark t-channel exchange feeding a resonant W boson, u d~ > mu+ vm g. It has one t-channel propagator and one resonant s-channel propagator, so the mapping can sample the W resonance directly instead of missing it.

Building the diagram

Diagram takes the incoming and outgoing masses, a list of Propagator and the vertices connecting them. A vertex is a triple of line labels: i for an incoming line, o for an outgoing one and p for a propagator, each followed by its 0-based index:

import numpy as np
import madspace as ms

M_W, W_W = 80.379, 2.085  # W mass and width, in GeV

diagram = ms.Diagram(
    [0.0, 0.0],            # incoming u, d~
    [0.0, 0.0, 0.0],       # outgoing g, mu+, vm
    [
        ms.Propagator(mass=0.0, width=0.0, pdg_id=2),   # t-channel u-quark
        ms.Propagator(mass=M_W, width=W_W, pdg_id=24),  # s-channel W+
    ],
    [["i0", "o0", "p0"], ["p0", "i1", "p1"], ["p1", "o1", "o2"]],
)

The first vertex is the incoming u quark radiating the gluon and turning into the t-channel propagator. The second vertex absorbs the d~ and produces the W. The third vertex is the W decay.

Inspecting the topology

Topology turns a diagram into one integration channel. It can be printed directly, which shows the t-channel propagator and the s-channel decay tree it builds:

topology = ms.Topology(diagram)
print(topology)

├── incoming: index=0, mass=0
├── incoming: index=1, mass=0
└── t-channel: (order=0, mass=0, width=0)
    ├── outgoing: index=0, mass=0
    └── decay: order=0, mass=80.379, width=2.085, e_min=0, e_max=0, pdg_id=24, on_shell=true, on_shell_boundary=false
        ├── outgoing: index=1, mass=0
        └── outgoing: index=2, mass=0

Sampling the resonance

Building a PhaseSpaceMapping from the topology samples the W invariant mass around its pole, using the propagator’s mass and width:

mapping = ms.PhaseSpaceMapping(topology, 13000.0)
rng = np.random.default_rng(0)
r = rng.random((100000, mapping.random_dim()))
momenta, x1, x2, det = mapping.map_forward([r])

p_mu, p_vm = momenta[:, 3], momenta[:, 4]
p_sum = p_mu + p_vm
m_inv = np.sqrt(np.abs(p_sum[:, 0] ** 2 - np.sum(p_sum[:, 1:] ** 2, axis=1)))
print(f"fraction within 5 GeV of the W mass: {(np.abs(m_inv - M_W) < 5).mean():.3f}")

Most of the sample lands close to the pole, even though the window is narrow compared to the full available energy:

fraction within 5 GeV of the W mass: 0.877

For comparison, a flat mapping built from the same external masses, with mode="rambo", has no notion of this resonance and only finds a tiny fraction of its events there:

flat = ms.PhaseSpaceMapping([0.0, 0.0, 0.0, 0.0, 0.0], 13000.0, mode="rambo")
momenta_flat, _, _, _ = flat.map_forward([r])
p_sum_flat = momenta_flat[:, 3] + momenta_flat[:, 4]
m_inv_flat = np.sqrt(np.abs(p_sum_flat[:, 0] ** 2 - np.sum(p_sum_flat[:, 1:] ** 2, axis=1)))
print(f"fraction within 5 GeV of the W mass: {(np.abs(m_inv_flat - M_W) < 5).mean():.3f}")

fraction within 5 GeV of the W mass: 0.009