diff --git a/docs/user/ufo-model-coverage.md b/docs/user/ufo-model-coverage.md index ba27450c7..a4a1f1f96 100644 --- a/docs/user/ufo-model-coverage.md +++ b/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. | diff --git a/src/pyamplicol/models/tensors.py b/src/pyamplicol/models/tensors.py index 05bd1d7a9..fa44a392f 100644 --- a/src/pyamplicol/models/tensors.py +++ b/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): diff --git a/tests/integration/test_ufo_chiral_scalar_coupling.py b/tests/integration/test_ufo_chiral_scalar_coupling.py new file mode 100644 index 000000000..28fe76630 --- /dev/null +++ b/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) diff --git a/tests/unit/test_model_tensors.py b/tests/unit/test_model_tensors.py index ef9db47c4..4a65db9da 100644 --- a/tests/unit/test_model_tensors.py +++ b/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