madspace: colour-ordered mapping undersamples two-particle sides (vs Fortran AmpliCol) #292

The colour-ordered mapping of madspace (ColorOrderedMapping, madspace/src/phasespace/color_ordered_mapping.cpp) produces large weights for colour orders where a side between the two beams has exactly two final-state particles. Fortran AmpliCol (PhaseSpace/phase_space_gen23.f03) handles these orders better. This issue summarises where the two mappings differ. Found with the [amplicol] backend of #291.

Example: g g > g g g g [amplicol], colour order (0, 2, 3, 1, 4, 5)

Setup: 13 TeV, NNPDF40_lo_as_01180, jet cuts pT > 20 GeV, |eta| < 5, ΔR > 0.4, invariant_power = 0.7. The process has five channels (colour-ordered subprocesses); subprocess 2 is the one with two gluons on each side of the beams:

No numerical problems: I evaluated the subprocess's mapping and its leading-colour matrix element on 2M points without VEGAS. There are no NaNs, infinities or negative weights.

• Weight distribution: the top 0.01% of points carry 99% of the variance; the largest weight is 1.3e5 times the mean (relative spread 154).
• Where the large weights sit: at low √ŝ (170–700 GeV), with jets near the pT cut. A same-side pair is nearly collinear (s₂₃ or s₄₅ ≈ 0.002–0.005 ŝ, ΔR ≈ 0.42–0.55), or a gluon is close to the opposite beam.

Share of accepted points vs share of σ, per colour-adjacent invariant (s_ij/ŝ):

Side 1 (beam 0 – 2 – 3 – beam 1) almost never has small invariants. Side 2 (beam 1 – 4 – 5 – beam 0) is sampled very differently, although the channel is symmetric between the two sides.

Differences with Fortran AmpliCol

Both mappings have the same overall structure:

• Sets: the final state is split into set 1 (between beam 0 and beam 1 in the colour order) and set 2 (after beam 1, wrapping round to beam 0).
• Central block: a 2→2 between the two composite systems if both sets have ≥ 2 particles, a double-t if one set has 1 particle, and a single chain if one set is empty.
• Sides with ≥ 3 particles: the first particle is peeled off the side's own beam with a t-channel 2→2. Each further particle comes from a 2→3 step with the opposite beam and the previously peeled particle as references.

They differ in the following points.

1. A side with exactly two particles, while the other side is non-empty (the main difference).
   • Fortran uses a full 2→3 step, gen23_one_step(set(i), inext, i, im1 = 3-i), with the opposite beam as the previous emission. It samples the t of the first particle with its own beam (s₀₂) and the invariant of the second particle with the opposite beam (s₃₁). With the composite mass s₂₃ from the central block, all three colour-adjacent invariants of the side are sampled directly.
   • madspace only does the first peel (_lab_scattering, a 2→2 off its own beam: t and a uniform φ). The second particle takes the remainder, so s₃₁ and s₅₀ are only reached through the φ integration.
2. Composite masses of the central 2→2 and intermediate rest masses.
   • Fortran samples them ∝ s^-0.5 (generate_mass). Each rest mass is drawn inside the step that uses it, with bounds from shatminmax given the momenta already generated and the ETmin cuts.
   • madspace samples them flat in s (_uniform_invariant is a default Invariant, power 0). All rest masses of a side are drawn up front with simple monotonic bounds. Set 1's mass is drawn first over the full range, and set 2's in the remainder (√ŝ − m₁)². This strongly biases set 1 to large and set 2 to small masses, and the coupling between the two ranges can't be undone by VEGAS.
3. Densities of the other invariants.
   • Fortran: the t and s of each step ∝ 1/x (power −1); ŝ ∝ ŝ^-1.2 with flat rapidity.
   • madspace: both 1/x^p with invariant_power (0.8 by default, 0.7 in the example).
4. Two-solution choice of the 2→3 steps.
   • Fortran computes both solutions, keeps the one that passes the pT cut (with weight ½), and picks at random only if both pass.
   • madspace exposes the choice as a discrete input, sampled uniformly by the Integrand, so the half that fails the cuts is wasted.

Subprocesses 0, 1, 3 and 4 use the single-chain or double-t branch, and their momentum transfers are sampled with a power law. That is why only the 0 2 3 1 4 5 order is this bad. The same issue applies to every colour order with a two-particle side and a non-empty opposite side, e.g. 0 2 3 1 4 in g g > g g g.

Experiment (no VEGAS, subprocess 2, 2M points)

The densities alone help, but 98% of the variance stays in the top 0.01% of points. This is consistent with point 1, which needs a structural change.

Proposed changes

1. Build a two-particle side with the existing 2→3 block (_two_to_three), using the opposite beam as the previous emission, as gen23_one_step(set(i), inext, i, 3-i) does.
2. Sample composite and rest masses ∝ s^-0.5. Ideally draw each rest mass inside its step with physical bounds, as Fortran does.
3. Use powers 1 for the t/s invariants and 1.2 for ŝ in colour-ordered mode.
4. Prefer the 2→3 solution that passes the cuts.

Also noticed: points rejected by the mapping (Jacobian −0.0) can carry unphysical momenta, e.g. a negative energy or a zero momentum. The Integrand doesn't evaluate matrix elements there, but code that evaluates every mapped point gets an error from UMAMI libraries that validate their input.