models: UFO vertex spinor chains keep their order and chirality (lone ProjM/ProjP were swapped)

docs/user/ufo-model-coverage.md

@@ -143,7 +143,7 @@ carry a `ValueError` message instead of a code; they are marked below.
 | Other tensor heads | rejected | `Epsilon`, `EpsilonBar`, `C`, `IdentityL`, three-argument `PSlash` and any accepted head with a wrong arity fail model compilation with `ValueError: ... contains unsupported UFO tensors`. |
 | Form factors | rejected | `form-factors`. |
 | Sigma normalisation | experimental | expanded as (i/2)(gamma^mu gamma^nu - gamma^nu gamma^mu); vertices using `Sigma` have not been validated against a reference implementation. |
-| Chiral structures | supported | `ProjM`, `ProjP` and `Gamma5` are validated through the packaged Standard Model. |
+| Chiral structures | supported | `ProjM`, `ProjP`, `Gamma5` and their products with `Gamma` and `Sigma` are checked against their UFO definitions (`tests/unit/test_model_tensors.py`); lone-projector scalar–fermion couplings are validated against MadGraph5_aMC@NLO (`tests/integration/test_ufo_chiral_scalar_coupling.py`). |
 | Index conventions | conditional | 1-based leg indices, the `1000*component + leg` notation for spin-2 fields and negative dummy indices are supported; inconsistent indices are a `ValueError`. |
 | Several structures per vertex | supported | any number of colour and Lorentz structures through the coupling matrix. |
 | Three-point vertices | supported | any combination of scalars, Dirac fermions, vectors and spin-2 fields, momentum-dependent or not, with colour 1, delta, T^a or f^abc. |

src/pyamplicol/models/tensors.py

@@ -181,64 +181,89 @@ def __init__(
         self.projm_tensor = TensorExpression.projm(4)
         self.projp_tensor = TensorExpression.projp(4)
         self._fresh_index = 0
-        # Standard UFO vertices index a fermion pair in the opposite order to
-        # the internal spinor slots, e.g. Gamma(3,2,1) for the SM FFV1.
-        dirac_legs = (
-            [leg for leg, spin in enumerate(spins, 1) if spin == 2] if vertex else []
-        )
-        self._spinor_partner = {
-            leg: partner
-            for first, second in zip(dirac_legs[::2], dirac_legs[1::2], strict=False)
-            for leg, partner in ((first, second), (second, first))
-        }
+        # A UFO vertex structure Gamma(mu,i,j)... is a matrix product with i the
+        # row (the outgoing-fermion index, e.g. particle 2 of FFV1's
+        # Gamma(3,2,1)) and j the column. spenso's bispinor tensors take
+        # (row, column), and each UFO spinor index stays on its own leg's
+        # slot, so the network evaluates the chain in its source order: the
+        # current of the row leg is the barred spinor, that of the column leg
+        # the unbarred one. A lone projector is therefore UFO's own chirality,
+        # Gamma(mu,2,-1)*ProjM(-1,1) is gamma^mu P_L, and products of several
+        # gamma matrices (Sigma, gamma^mu gamma^5) keep their order.
+        # Propagators keep their own (column, row) order.
+        self._vertex = vertex
 
     def metric(self, match: Mapping[Expression, Expression]) -> Expression:
         left = self._minkowski_slot(match[_WILDCARDS["a_"]])
         right = self._minkowski_slot(match[_WILDCARDS["b_"]])
         return TensorExpression.g(self.minkowski)(left, right).to_expression()
 
+    def _spinor_pair(self, row: Expression, column: Expression) -> tuple[str, str]:
+        """spenso's (first, second) bispinor labels of a UFO (row, column)
+        index pair: (row, column) at vertices, (column, row) in propagators."""
+
+        first = self._bispinor_label(row)
+        second = self._bispinor_label(column)
+        return (first, second) if self._vertex else (second, first)
+
     def gamma(self, match: Mapping[Expression, Expression]) -> Expression:
         lorentz = self._minkowski_label(match[_WILDCARDS["a_"]])
-        output = self._bispinor_label(match[_WILDCARDS["b_"]])
-        input_ = self._bispinor_label(match[_WILDCARDS["c_"]])
-        return _as_expression(self.gamma_tensor(input_, output, lorentz))
+        first, second = self._spinor_pair(
+            match[_WILDCARDS["b_"]], match[_WILDCARDS["c_"]]
+        )
+        return _as_expression(self.gamma_tensor(first, second, lorentz))
 
     def gamma5(self, match: Mapping[Expression, Expression]) -> Expression:
-        output = self._bispinor_label(match[_WILDCARDS["a_"]])
-        input_ = self._bispinor_label(match[_WILDCARDS["b_"]])
-        return _as_expression(self.gamma5_tensor(input_, output))
+        first, second = self._spinor_pair(
+            match[_WILDCARDS["a_"]], match[_WILDCARDS["b_"]]
+        )
+        return _as_expression(self.gamma5_tensor(first, second))
 
     def projector_minus(
         self,
         match: Mapping[Expression, Expression],
     ) -> Expression:
-        output = self._bispinor_label(match[_WILDCARDS["a_"]])
-        input_ = self._bispinor_label(match[_WILDCARDS["b_"]])
-        return _as_expression(self.projp_tensor(input_, output))
+        """UFO ``ProjM`` = (1 - gamma^5)/2, spenso's ``projm`` at vertices.
+        Propagators keep their historical exchanged head (no supported
+        propagator uses a projector)."""
+
+        first, second = self._spinor_pair(
+            match[_WILDCARDS["a_"]], match[_WILDCARDS["b_"]]
+        )
+        tensor = self.projm_tensor if self._vertex else self.projp_tensor
+        return _as_expression(tensor(first, second))
 
     def projector_plus(
         self,
         match: Mapping[Expression, Expression],
     ) -> Expression:
-        output = self._bispinor_label(match[_WILDCARDS["a_"]])
-        input_ = self._bispinor_label(match[_WILDCARDS["b_"]])
-        return _as_expression(self.projm_tensor(input_, output))
+        """UFO ``ProjP`` = (1 + gamma^5)/2 (``projector_minus``)."""
+
+        first, second = self._spinor_pair(
+            match[_WILDCARDS["a_"]], match[_WILDCARDS["b_"]]
+        )
+        tensor = self.projp_tensor if self._vertex else self.projm_tensor
+        return _as_expression(tensor(first, second))
 
     def identity(self, match: Mapping[Expression, Expression]) -> Expression:
-        left = self.bispinor(self._bispinor_label(match[_WILDCARDS["a_"]]))
-        right = self.bispinor(self._bispinor_label(match[_WILDCARDS["b_"]]))
+        first, second = self._spinor_pair(
+            match[_WILDCARDS["b_"]], match[_WILDCARDS["a_"]]
+        )
+        left = self.bispinor(first)
+        right = self.bispinor(second)
         return TensorExpression.g(self.bispinor)(left, right).to_expression()
 
     def sigma(self, match: Mapping[Expression, Expression]) -> Expression:
         mu = self._minkowski_label(match[_WILDCARDS["a_"]])
         nu = self._minkowski_label(match[_WILDCARDS["b_"]])
-        output = self._bispinor_label(match[_WILDCARDS["c_"]])
-        input_ = self._bispinor_label(match[_WILDCARDS["d_"]])
+        # (gamma^mu gamma^nu)[row][column] = sum_x gamma^mu[row][x] gamma^nu[x][column]
+        # in spenso's (first, second) order (``_spinor_pair``).
+        start, end = self._spinor_pair(match[_WILDCARDS["c_"]], match[_WILDCARDS["d_"]])
         dummy = self._fresh_label("sigma")
-        first = _as_expression(self.gamma_tensor(input_, dummy, mu))
-        second = _as_expression(self.gamma_tensor(dummy, output, nu))
-        reverse_first = _as_expression(self.gamma_tensor(input_, dummy, nu))
-        reverse_second = _as_expression(self.gamma_tensor(dummy, output, mu))
+        first = _as_expression(self.gamma_tensor(start, dummy, mu))
+        second = _as_expression(self.gamma_tensor(dummy, end, nu))
+        reverse_first = _as_expression(self.gamma_tensor(start, dummy, nu))
+        reverse_second = _as_expression(self.gamma_tensor(dummy, end, mu))
         return E("1𝑖/2") * (  # noqa: RUF001
             first * second - reverse_first * reverse_second
         )
@@ -272,12 +297,13 @@ def propagator_momentum(
         return self._momentum_tensor(leg, self._minkowski_slot(index))
 
     def p_slash(self, match: Mapping[Expression, Expression]) -> Expression:
-        output = self._bispinor_label(match[_WILDCARDS["a_"]])
-        input_ = self._bispinor_label(match[_WILDCARDS["b_"]])
+        first, second = self._spinor_pair(
+            match[_WILDCARDS["a_"]], match[_WILDCARDS["b_"]]
+        )
         lorentz_label = self._fresh_label("pslash")
         lorentz = self.minkowski(lorentz_label)
         return _as_expression(
-            self.gamma_tensor(input_, output, lorentz_label)
+            self.gamma_tensor(first, second, lorentz_label)
         ) * self._momentum_tensor(
             1,
             lorentz,
@@ -303,7 +329,7 @@ def _bispinor_label(self, index: Expression) -> str:
             )
         if component != 1:
             raise ValueError(f"Dirac particle leg {leg} has no spin index {component}")
-        return f"ufo_s_{component}_{self._spinor_partner.get(leg, leg)}"
+        return f"ufo_s_{component}_{leg}"
 
     def _momentum_tensor(self, leg: int, slot) -> Expression:
         if leg < 1 or leg > len(self.spins):

tests/integration/test_ufo_chiral_scalar_coupling.py

@@ -0,0 +1,129 @@
+# SPDX-License-Identifier: 0BSD
+"""A charged scalar with chiral UFO lone-projector couplings against MadGraph.
+
+The packaged sm UFO plus H+ with ``H+ bbar (yHL P_L + yHR P_R) t + h.c.``,
+written with the sm's own FFS1 = ProjM(2,1) and FFS3 = ProjP(2,1). In
+``b b~ > h+ w-`` (one diagram, a t-channel top between the chiral W vertex and
+the H+ vertex) |M|^2 = yHL^2 A + yHR^2 B with A != B, so the helicity sum
+tells the two chiralities apart.
+"""
+
+from __future__ import annotations
+
+import importlib.util
+import math
+import os
+import shutil
+from pathlib import Path
+from typing import NoReturn
+
+import pytest
+
+from pyamplicol import Generator, ModelSource, Runtime
+from pyamplicol.config import (
+    ColorConfig,
+    EvaluatorConfig,
+    GenerationConfig,
+    GenerationValidationConfig,
+    JITConfig,
+    RunConfig,
+)
+
+_SM_UFO = (
+    Path(__file__).resolve().parents[2]
+    / "src"
+    / "pyamplicol"
+    / "assets"
+    / "models"
+    / "ufo"
+    / "sm"
+)
+_ADDITIONS = {
+    "parameters.py": """
+MHp = Parameter(name='MHp', nature='external', type='real', value=300.,
+                texname='MHp', lhablock='MASS', lhacode=[37])
+WHp = Parameter(name='WHp', nature='external', type='real', value=0.,
+                texname='WHp', lhablock='DECAY', lhacode=[37])
+yHL = Parameter(name='yHL', nature='external', type='real', value=1.0,
+                texname='yHL', lhablock='HPLUS', lhacode=[1])
+yHR = Parameter(name='yHR', nature='external', type='real', value=0.2,
+                texname='yHR', lhablock='HPLUS', lhacode=[2])
+""",
+    "particles.py": """
+H__plus__ = Particle(pdg_code=37, name='H+', antiname='H-', spin=1, color=1,
+                     mass=Param.MHp, width=Param.WHp, texname='H+',
+                     antitexname='H-', charge=1, GhostNumber=0,
+                     LeptonNumber=0, Y=0)
+H__minus__ = H__plus__.anti()
+""",
+    "couplings.py": """
+GC_HPL = Coupling(name='GC_HPL', value='complex(0,1)*yHL', order={'QED':1})
+GC_HPR = Coupling(name='GC_HPR', value='complex(0,1)*yHR', order={'QED':1})
+""",
+    "vertices.py": """
+V_HP1 = Vertex(name='V_HP1', particles=[P.t__tilde__, P.b, P.H__plus__],
+               color=['Identity(1,2)'], lorentz=[L.FFS1, L.FFS3],
+               couplings={(0,0):C.GC_HPL,(0,1):C.GC_HPR})
+V_HP2 = Vertex(name='V_HP2', particles=[P.b__tilde__, P.t, P.H__minus__],
+               color=['Identity(1,2)'], lorentz=[L.FFS1, L.FFS3],
+               couplings={(0,0):C.GC_HPR,(0,1):C.GC_HPL})
+""",
+}
+# sqrt(s) = 1 TeV, on shell for MB = 4.7, MHp = 300 and the model's MW.
+_POINT = (
+    (500.0, 0.0, 0.0, 499.97790951201034),
+    (500.0, 0.0, 0.0, -499.97790951201034),
+    (541.7663922180795, -193.63353947054688, 216.71729887270064, 345.03693787538566),
+    (458.23360778192045, 193.63353947054688, -216.71729887270064, -345.03693787538566),
+)
+# MadGraph5_aMC@NLO 3.7.0 standalone smatrix at _POINT for the same UFO (the
+# sm defaults, MT = 173, at full precision) and (yHL, yHR).
+_MADGRAPH = {
+    (1.0, 0.2): 2.5578120350182,
+    (0.2, 1.0): 0.11322998418504657,
+}
+
+
+def _unavailable(reason: str) -> NoReturn:
+    if os.environ.get("PYAMPLICOL_REQUIRE_NATIVE_TESTS") == "1":
+        pytest.fail(reason)
+    pytest.skip(reason)
+
+
+def test_lone_projector_couplings_keep_their_chirality(tmp_path: Path) -> None:
+    if importlib.util.find_spec("pyamplicol._rusticol") is None:
+        _unavailable("the Rusticol extension has not been built")
+    model = tmp_path / "sm_hplus"
+    shutil.copytree(_SM_UFO, model, ignore=shutil.ignore_patterns("__pycache__"))
+    for card in model.glob("restrict_*.dat"):
+        card.unlink()
+    for name, text in _ADDITIONS.items():
+        with (model / name).open("a", encoding="utf-8") as stream:
+            stream.write(text)
+    compiled = ModelSource.from_path(model, restriction="none").compile(
+        use_cache=False, require_supported=True
+    )
+    artifact = tmp_path / "bb_hpw"
+    Generator(
+        RunConfig(
+            action="generate",
+            color=ColorConfig(accuracy="full"),
+            generation=GenerationConfig(
+                workers=1,
+                emit_api_bundle=False,
+                validation=GenerationValidationConfig(
+                    enabled=False, post_build_validation=False
+                ),
+            ),
+            evaluator=EvaluatorConfig(
+                execution_mode="compiled", jit=JITConfig(optimization_level=1)
+            ),
+        )
+    ).generate("b b~ > h+ w-", artifact, model=compiled)
+    runtime = Runtime.load(artifact)
+
+    for (left, right), expected in _MADGRAPH.items():
+        runtime.set_model_parameters({"yHL": left, "yHR": right})
+        value = complex(runtime.evaluate((_POINT,))[0])
+        assert math.isclose(value.real, expected, rel_tol=1e-10), (left, right)
+        assert abs(value.imag) <= 1e-12 * abs(value.real)

tests/unit/test_model_tensors.py

@@ -76,13 +76,13 @@ def test_model_tensors_keep_supported_singlet_projection() -> None:
         (
             "UFO::ProjM(UFO::idx(1,1),UFO::idx(1,2))",
             (2, 2),
-            "projp",
+            "projm",
             ("ufo_s_1_1", "ufo_s_1_2"),
         ),
         (
             "UFO::ProjP(UFO::idx(1,1),UFO::idx(1,2))",
             (2, 2),
-            "projm",
+            "projp",
             ("ufo_s_1_1", "ufo_s_1_2"),
         ),
         (
@@ -133,6 +133,77 @@ def gamma(first: int, second: int, *, vertex: bool = True) -> object:
     assert gamma(2, 1) != gamma(2, 1, vertex=False)
 
 
+@pytest.mark.parametrize(
+    ("left", "right"),
+    (
+        # A lone projector is UFO's own chirality (FFS1/FFS3 couplings).
+        ("ProjM(2,1)", "1/2*(Identity(2,1)-Gamma5(2,1))"),
+        ("ProjP(2,1)", "1/2*(Identity(2,1)+Gamma5(2,1))"),
+        # A pseudoscalar written with projectors or with gamma^5 (UFO FFS2).
+        ("ProjM(2,1)-ProjP(2,1)", "-Gamma5(2,1)"),
+        ("ProjM(2,1)+ProjP(2,1)", "Identity(2,1)"),
+        # gamma^mu gamma^5 (the axial current) both ways: the chain keeps its
+        # order next to a projector difference as next to Gamma5.
+        (
+            "Gamma(3,2,-1)*ProjP(-1,1)-Gamma(3,2,-1)*ProjM(-1,1)",
+            "Gamma(3,2,-1)*Gamma5(-1,1)",
+        ),
+        ("Gamma(3,2,-1)*ProjM(-1,1)", "ProjP(2,-1)*Gamma(3,-1,1)"),
+        # Sigma^{mu nu} as its gamma-matrix commutator.
+        (
+            "Sigma(3,4,2,1)",
+            "1/2*complex(0,1)*(Gamma(3,2,-1)*Gamma(4,-1,1)-Gamma(4,2,-1)*Gamma(3,-1,1))",
+        ),
+        # A chiral dipole coupling sigma^{mu nu} P_L.
+        (
+            "Sigma(3,4,2,-1)*ProjM(-1,1)",
+            "1/2*(Sigma(3,4,2,1)-Sigma(3,4,2,-1)*Gamma5(-1,1))",
+        ),
+    ),
+)
+def test_vertex_spinor_chains_keep_their_order_and_chirality(
+    left: str, right: str
+) -> None:
+    """UFO identities between Lorentz structures hold after lowering: a lone
+    projector keeps its chirality and a product of gamma matrices keeps its
+    order, so scalar couplings written with projectors (FFS1-3) equal their
+    gamma^5 form, and chiral currents their gamma^5 form."""
+
+    from symbolica import E
+    from symbolica.community.tensor import TensorLibrary, TensorNetwork, as_tensor
+    from ufo_model_loader.symbolica_processing import (
+        parse_python_expression_safe,
+        wrap_indices,
+    )
+
+    from pyamplicol.models.compiler_kernels import _ordered_dense_tensor_components
+
+    spins = (2, 2, 3, 3) if "4" in left else (2, 2, 3) if "(3" in left else (2, 2, 1)
+    axes = (
+        "ufo_s_1_1",
+        "ufo_s_1_2",
+        *(f"ufo_l_1_{leg}" for leg, spin in enumerate(spins, 1) if spin == 3),
+    )
+
+    def dense(structure: str) -> tuple[object, ...]:
+        wrapped = wrap_indices(parse_python_expression_safe(structure))
+        source = wrapped.to_canonical_string()
+        expression = E(normalize_lorentz_expression(source, spins).expression)
+        library = TensorLibrary.hep_lib_atom()
+        network = TensorNetwork(as_tensor(expression), library)
+        network.execute(library=library)
+        result = network.result_tensor(library)
+        return tuple(
+            E(str(value)).expand()
+            for value in _ordered_dense_tensor_components(result, axes).values
+        )
+
+    assert all(
+        (mine - theirs).expand() == E("0")
+        for mine, theirs in zip(dense(left), dense(right), strict=True)
+    )
+
+
 def test_lorentz_clifford_contraction_preserves_4d_normalization() -> None:
     from symbolica import E