From d7269f4075d9105278b14cb1c389e7a2caa476f8 Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Wed, 9 Sep 2026 23:04:46 +0200 Subject: [PATCH 1/3] MadSpin: six fixes for the frame-projected polarised density Rescued verbatim from the worktree they were developed in, on the base they were written against (PR #103's head). Not yet ported to main -- that base is ~100 commits behind and the gap contains PR #63, the MadSpin performance rewrite, so the port is a real one and is kept separate from preserving the work. Together these are what make a polarised NLO production decayable at all; before them `p p > z{0} z{0} [QCD]` + MadSpin ran at ~450 s per event and was killed at event 40 of 4000. After, 20000 events decay in 15.6 s. 1. _production_polarization parsed the process line with MadSpin's tree-level mg5cmd, which refuses the perturbation-order bracket. MadSpin only warned, so the density was left unrestricted and every NLO polarised production got UNPOLARISED decay angles -- with the right cross section and no crash. Strip the bracket; it is irrelevant to the polarisation. 2. _frame_boost called get_momenta without merged_map, raising ValueError under flavour grouping. 3. get_density called get_pdg (which matches by exact momentum equality) on the BOOSTED momenta, so no particle was found. Read the pdgs off the lab momenta before boosting. 4. get_density's cached branch recomputed the tag with a bare get_tag_and_order(), losing get_pdir's 1->N antiparticle fallback. Cache the tag get_pdir resolved instead. 5. The offshell mass-stage weight and the joint weight divided a me_frame-boosted Tr(rho_off) by a LAB-frame |M_prod|^2. A helicity-restricted matrix element is not Lorentz invariant, so the two are different projections; for an unpolarised production the trace is invariant and they agree, which is why only polarised runs were affected. _onshell_production_norm returns calculate_matrix_ element unchanged when there is no frame boost -- so unpolarised runs are bit-for-bit identical -- and otherwise the trace of the on-shell rho in the same frame. This is the load-bearing one: reverting it alone puts the mass-stage bound back to 7.9e6 while leaving the two angular bounds bit-identical. 6. The same frame fix applied to the joint path. Verified safe there: over 5414 production events the denominator was bit-identical across every joint trial of the same event, so it cancels out of the accept/reject; Route B's published numbers move by <= 0.5 sigma. Upstream 3.8.0 carries an independent version of fix 5 (mg5amcnlo#407) and of the joint symmetry factor (#408/#409). When MadGraph7 merges 3.8.0 these need reconciling rather than both being applied -- the MadGraph7 get_density signature carries flavour-grouping arguments that 3.8.0 has no equivalent of. Co-Authored-By: Claude Opus 5 --- MadSpin/interface_madspin.py | 137 ++++++++++++++++++++++++++++++++--- 1 file changed, 128 insertions(+), 9 deletions(-) diff --git a/MadSpin/interface_madspin.py b/MadSpin/interface_madspin.py index 979858fdf2..12bc05ed7a 100755 --- a/MadSpin/interface_madspin.py +++ b/MadSpin/interface_madspin.py @@ -7501,8 +7501,16 @@ def _production_polarization(self): source = {} multi_id = set() for line in lines: + # The perturbation-order bracket of an NLO process line + # ('p p > z{0} z{0} [QCD]') is irrelevant to the polarisation + # braces, which sit on the legs -- but MadSpin's mg5cmd holds + # the TREE-level model, so extract_process refuses it with + # "Perturbation order QCD is not among the perturbation orders + # allowed for by the loop model" and every NLO polarised + # production silently lost its restriction. Strip it. + bare = re.sub(r'\[[^\]]*\]', '', line).strip() try: - procdef = self.mg5cmd.extract_process(line) + procdef = self.mg5cmd.extract_process(bare) except Exception as error: logger.warning('MadSpin could not re-read the polarisation of ' 'the production process "%s" (%s); the density ' @@ -8801,6 +8809,52 @@ def _upfront_production(self, production, order, particles, slot_to_index, parents = {slot: finals[slot_to_index[slot]] for slot in order} return rho_off, jac_reshuffle, slot_mass, parents, frame_boost + def _onshell_production_norm(self, production, prod_static): + """|M_prod|^2 on shell, in the frame the density modes quantise in. + + The denominator of the offshell mass-set weight. Without a frame boost + (the unpolarised case) this is exactly ``calculate_matrix_element`` -- + same value, same normalisation -- so nothing changes there. With one + (a polarisation brace, polarised beams, the pure-interference mode) the + matrix element is *not* frame independent and the lab-frame value is a + different projection from the one Tr(rho_off) is built from, so the + trace of the on-shell rho *in that frame* is taken instead. + + The two sides agree where they must: with the boost left out + (``frame_boost=None``) this trace reproduces ``calculate_matrix_element`` + to 1.4e-7 at worst over 600 `p p > z z [QCD]` production events -- they + are the same helicity sum through two entry points. What that measurement + also showed is that the *boosted* trace is not that number: with no + helicity restriction at all it is analytically the full (Lorentz + invariant) helicity sum, yet the boosted evaluation returns a value that + is <= the lab one, by more than 1e-3 on 5.8% of production events and by + more than 10% on 0.8%, with only a weak correlation with the size of the + boost. That deficit lives in the boosted density evaluation itself, not + in this choice of denominator -- it is already inside the *numerator* of + every ``me_frame`` run -- and taking the denominator in the same frame is + what makes it cancel out of the accept/reject weight rather than + multiply it. Its origin is not understood: GET_DENSITY applies no + good-helicity filter and sums the whole NHEL table, so the trace ought + to be invariant. Worth a separate look. + + Evaluated on a round-tripped copy, like ``_upfront_production``'s + ``prod_off``, so the two sides of the ratio see the same %.10e + truncation -- and so that the ``_ms_orig_order_for_density`` cache is + never left on the shared production event. + """ + frame_boost = self._frame_boost(production) + if frame_boost is None: + return self.calculate_matrix_element(production) + prod_on = lhe_parser.Event(str(production)) + rho_on = self.get_density(prod_on, prod_static['position'], + prod_static['allowed_hel'], + prod_static['ncomb'], + prod_static['dimension'], + frame_boost=self._frame_boost(prod_on), + hel_restriction=prod_static.get('hel_restriction'), + hel_restriction_trace=prod_static.get('hel_restriction_trace')) + return rho_on.trace().real + def _sequential_offshell(self): """Whether the sequential accept/reject runs its offshell (madspin/full) branch: the production density is evaluated at reshuffled momenta, so the @@ -9487,12 +9541,29 @@ def sequential_accept_reject(self, production, evt_decayfile, maxwgts, # the production density and by every decay contracted against it. The # offshell branch gets its own from _upfront_production, derived from # the reshuffled production rho is evaluated at. + # + # It has to be the *same* quantity as the numerator, in the *same* + # frame. calculate_matrix_element hands the matrix element the LAB + # momenta, and a polarised matrix element is not Lorentz invariant: for + # a braced production (p p > z{0} z{0}) the lab-frame projection is a + # different object from the me_frame one Tr(rho_off) is built from, and + # it can be six orders of magnitude smaller on a boosted event (measured + # on p p > z{0} z{0} [QCD], 25 probe events: Tr(rho) in the lab spans + # 1.8e7 against 1.4e4 in the me_frame, with single events at 4.4e-10 + # against 2.2e-3 in the frame). The probe's single bound is then set by + # those events -- 8.1e6 against a median mass-set weight of 29 -- and the + # mass stage of every ordinary event accepts at ~4e-6, i.e. ~1e5 redraws + # and ~450 s per decayed event. Taking the trace of the on-shell rho in + # the frame instead makes w_mass exactly the offshell/onshell ratio. + # Unpolarised runs have no frame boost (_frame_boost returns None) and + # keep the matrix-element call, so they are bit-for-bit unchanged. frame_boost = None me_prod_on = 1.0 if offshell: me_prod_on = getattr(production, 'me_wgt', None) if not me_prod_on: - me_prod_on = self.calculate_matrix_element(production) + me_prod_on = self._onshell_production_norm(production, + prod_static) production.me_wgt = me_prod_on if not me_prod_on: # a production event with no matrix element cannot be normalised @@ -10414,10 +10485,38 @@ def calculate_matrix_element_from_density(self, production, decays, decay_dict, #VALENTIN: except for the mode "full", we should not compute the matrix element here MEdenom_prod, MEdenom_decay = None, None if not density_pole_approximation: - # compute the denominator and then reshuffle the event before - # computing the numerator - MEdenom_prod = self.calculate_matrix_element(production) - MEdenom_decay = 1.0 + # compute the denominator and then reshuffle the event before + # computing the numerator + # + # The denominator has to be the same quantity as the numerator, + # in the same frame: the numerator is the contraction of the + # (possibly helicity-restricted) production density built in the + # me_frame, and a *restricted* matrix element is not Lorentz + # invariant, so the lab-frame calculate_matrix_element is a + # different projection. _onshell_production_norm returns exactly + # calculate_matrix_element when there is no frame boost, and the + # trace of the on-shell rho in the frame when there is. Same fix + # as on the sequential mass stage. + # + # This call site is NOT polarised-only: an unpolarised + # production that only asks for keep_weight_for_polarization_* + # still switches the frame axis on, and 'unweighting = auto' + # sends an unbraced production here. That case is safe, and + # measured to be so rather than assumed. The denominator is a + # per-production-event constant either way (checked directly: + # over 5414 production events of `p p > z z [QCD]` the value was + # bit-identical across every joint trial of the same event), and + # a constant factor cancels out of an accept/reject -- so the + # accepted decay distribution is unchanged in law. End to end, + # redecaying 150000 production events of one such sample here + # gave f0 = 0.18199 +- 0.00054, f00 0.05857, fTT 0.69460, Ckk + # -0.57691 against 0.18176 +- 0.00054 / 0.05830 / 0.69478 / + # -0.57706 from the pre-fix run on the SAME production events; + # and a seed-matched 45000-event pair decayed both ways cost + # 10.09 against 10.15 accept/reject trials per event. + MEdenom_prod = self._onshell_production_norm(production, + prod_static) + MEdenom_decay = 1.0 for key in decays: for dec in decays[key]: MEdenom_decay *= self.calculate_matrix_element(dec) @@ -10762,7 +10861,14 @@ def _frame_boost(self, event): if frame_id <= 0: return None _, orig_order, _, _, _ = self.get_pdir(event) - momenta = event.get_momenta(orig_order) + # merged_map: orig_order comes out of all_me, which is keyed by the + # MERGED-pdg tag when apply_flavor_grouping is on (81/-81 in place of + # u/d/...). Without the map the raw event pdgs are looked up in a + # merged block and get_mapping dies with + # "ValueError: list.index(x): x not in list" -- exactly as the matrix + # element call twenty lines below already guards against. + momenta = event.get_momenta(orig_order, + merged_map=self._revert_merged or None) selected = [n for n in range(1, len(momenta) + 1) if frame_id >> n & 1] if not selected: return None @@ -10837,8 +10943,17 @@ def get_density(self, event, position, allow_hel, ncomb, dimension, if orig_order is None: _, orig_order, _, _, tag = self.get_pdir(event) event._ms_orig_order_for_density = orig_order + # cache the tag get_pdir resolved rather than recomputing it: it is + # the MERGED-pdg tag (all_me is keyed by 81/-81, not by the raw + # event pdgs, so a bare get_tag_and_order() raises KeyError on the + # lookup below -- e.g. ((-1, 1), (23, 23)) against ((-81, 81), + # (23, 23))), and get_pdir also owns the 1 -> N antiparticle + # fallback, which no recomputation here would reproduce. + event._ms_tag_for_density = tag else: #in any case, we need tag to differentiate between production and decay - tag, _ = event.get_tag_and_order() + tag = getattr(event, '_ms_tag_for_density', None) + if tag is None: + tag, _ = event.get_tag_and_order(self._revert_merged or None) # Fast path: single-point momentum extraction without permutation @@ -10853,6 +10968,10 @@ def get_density(self, event, position, allow_hel, ncomb, dimension, all_p = event.get_all_momenta(orig_order, merged_map=self._revert_merged or None) assert len(all_p) == 1, "Error: get_density can only be called for a single phase-space point" p = all_p[0] + # get_pdg below identifies the particles by EXACT momentum equality + # against the event record, so it has to see the lab-frame momenta: + # keep them before the frame boost rewrites every component. + p_lab = p if frame_boost is not None: p = self._boost_momenta(p, frame_boost, rest_leg=frame_rest_leg) P = rwgt_interface.ReweightInterface.invert_momenta(p) @@ -10866,7 +10985,7 @@ def get_density(self, event, position, allow_hel, ncomb, dimension, need_raw_pdg = (self._revert_merged and any(abs(pid) in merged_particles for pid in pdg_template)) if need_raw_pdg: - pdgs = event.get_pdg(p) + pdgs = event.get_pdg(p_lab) else: pdgs = pdg_template n_changing = len(position) From 02cda9c636bd718396f05e3267080bc1e81d5a09 Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Wed, 9 Sep 2026 23:15:51 +0200 Subject: [PATCH 2/3] tests: let the MadSpin stubs follow the interface they double The six density fixes were written against PR #103's head; #65's head has nine commits on top, three of which added stubs that the fixes now outgrow. All three are the double falling behind the real class, not a problem with the fixes, so the stubs move rather than the production code: * _FrameStub gains _revert_merged. MadSpinInterface carries it as a class attribute ({} unless flavour grouping is on) and _frame_boost now passes it to get_momenta as merged_map; the stub does not inherit from MadSpinInterface, so it did not have it. Guarding the production side with getattr would have hidden a real attribute behind a test-only default. * _MomentaEvent.get_momenta gains merged_map=None, matching lhe_parser.Event.get_momenta, which has carried it all along. * _offshell_fixture gains _onshell_production_norm, delegating to the calculate_matrix_element it already fakes. That is exactly what the real helper does when there is no frame boost, which is the only case a stub with no boost can represent -- so the fixture's own comment, "|M_prod|^2 on shell, the denominator of the offshell mass-set weight", still describes what the test measures. Borrowing the real helper instead would have dragged _frame_boost and get_density into the stub's hand-kept borrow list. 502 tests, OK -- the same count and result as the base without the fixes, so nothing was skipped to get there. Co-Authored-By: Claude Opus 5 --- tests/unit_tests/madspin/test_madspin.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/madspin/test_madspin.py b/tests/unit_tests/madspin/test_madspin.py index e41a0355b4..1625045644 100755 --- a/tests/unit_tests/madspin/test_madspin.py +++ b/tests/unit_tests/madspin/test_madspin.py @@ -438,6 +438,10 @@ def __init__(self, frame_id, beampol, prodpol=None, vector=(), fermion=(), self.options['keep_weight_for_polarization_fermion'] = list(fermion) self.options['pure_interference'] = pure_interference self.model = _PIModelStub() + # MadSpinInterface carries this as a class attribute (it is {} unless + # flavour grouping is on); _frame_boost passes it to get_momenta as + # merged_map, so the stub has to have it too. + self._revert_merged = {} # what _production_polarization would have parsed out of the banner's # proc_card: {} for a brace-free production process self._production_polarization_cache = prodpol if prodpol else {} @@ -467,7 +471,9 @@ class _MomentaEvent(object): def __init__(self, momenta): self.momenta = momenta - def get_momenta(self, orig_order): + def get_momenta(self, orig_order, merged_map=None): + # merged_map mirrors lhe_parser.Event.get_momenta: _frame_boost has + # to pass it so the ME ordering resolves under flavour grouping. return self.momenta @@ -6703,6 +6709,12 @@ def _offshell_fixture(self): stub._slot_of = {index: slot for slot, index in enumerate(slots)} # |M_prod|^2 on shell, the denominator of the offshell mass-set weight stub.calculate_matrix_element = lambda *args, **opts: 1.0 + # The mass stage takes that denominator through + # _onshell_production_norm, which returns calculate_matrix_element + # unchanged when there is no frame boost -- which is this stub's + # case, and the only one it can represent. + stub._onshell_production_norm = \ + lambda production, prod_static: stub.calculate_matrix_element(production) def _no_pool(*args, **opts): raise TestPAUpFrontMass._NoDecayPool() From 309691166777c0a401f4128c8f254bf6c8ff24d5 Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Mon, 14 Sep 2026 23:43:00 +0200 Subject: [PATCH 3/3] MadSpin: address the #119 review - Test the perturbation-order bracket stripping in _production_polarization: the stub parser now refuses a bracket the way MadSpin's tree-level mg5cmd does, and the new test fails ({} != ...) with the stripping reverted. - Trim the long benchmark comments on the sequential mass stage, the joint path and _onshell_production_norm down to the invariant; drop the reference to a FINDINGS.md that is not in the repository. The measurements stay in the PR description. - get_tag_and_order: pass the merged map by keyword. - p_lab: say why a reference is enough (_boost_momenta returns new tuples and never mutates its input) instead of copying. madspin unit tests: 609 OK. Co-Authored-By: Claude Opus 5 --- MadSpin/interface_madspin.py | 82 +++++++----------------- tests/unit_tests/madspin/test_madspin.py | 19 ++++++ 2 files changed, 42 insertions(+), 59 deletions(-) diff --git a/MadSpin/interface_madspin.py b/MadSpin/interface_madspin.py index b952a7e6b1..2903e4876e 100755 --- a/MadSpin/interface_madspin.py +++ b/MadSpin/interface_madspin.py @@ -9925,16 +9925,9 @@ def _onshell_production_norm(self, production, prod_static): ``prod_off``, so the two sides of the ratio see the same %.10e truncation. - With the boost left out (``frame_boost=None``) this trace reproduces - ``calculate_matrix_element`` to 1.4e-7 at worst over 600 `p p > z z - [QCD]` production events -- the same helicity sum through two entry - points. The *boosted* trace is a separate, open question: with no - helicity restriction at all it is analytically that same Lorentz - invariant sum, yet it comes back <= the lab value, by more than 1e-3 - on 5.8% of production events. That deficit is in the boosted density - evaluation, not in this choice of denominator -- it is already inside - the *numerator* of every ``me_frame`` run, and taking the denominator - in the same frame is what makes it cancel. See FINDINGS.md. + Any residual inaccuracy of the boosted density evaluation is already + inside the numerator of every ``me_frame`` run; taking the denominator + in the same frame is what makes it cancel. """ frame_boost = self._frame_boost(production) if frame_boost is None: @@ -10643,33 +10636,19 @@ def sequential_accept_reject(self, production, evt_decayfile, maxwgts, # the production density and by every decay contracted against it. The # offshell branch gets its own from _upfront_production, derived from # the reshuffled production rho is evaluated at. - # - # It has to be the *same* quantity as the numerator, in the *same* - # frame. calculate_matrix_element hands the matrix element the LAB - # momenta, and a polarised matrix element is not Lorentz invariant: for - # a braced production (p p > z{0} z{0}) the lab-frame projection is a - # different object from the me_frame one Tr(rho_off) is built from, and - # it can be six orders of magnitude smaller on a boosted event (measured - # on p p > z{0} z{0} [QCD], 25 probe events: Tr(rho) in the lab spans - # 1.8e7 against 1.4e4 in the me_frame, with single events at 4.4e-10 - # against 2.2e-3 in the frame). The probe's single bound is then set by - # those events -- 8.1e6 against a median mass-set weight of 29 -- and the - # mass stage of every ordinary event accepts at ~4e-6, i.e. ~1e5 redraws - # and ~450 s per decayed event. Taking the trace of the on-shell rho in - # the frame instead makes w_mass exactly the offshell/onshell ratio. - # Unpolarised runs have no frame boost (_frame_boost returns None) and - # keep the matrix-element call, so they are bit-for-bit unchanged. frame_boost = None me_prod_on = 1.0 if offshell: me_prod_on = getattr(production, 'me_wgt', None) if not me_prod_on: # The denominator has to be the same quantity as the numerator, - # in the same frame: Tr(rho_off) is built in the me_frame while - # calculate_matrix_element hands the matrix element the LAB - # momenta, and a helicity-restricted matrix element is not - # Lorentz invariant. Unpolarised runs have no frame boost and - # keep the matrix-element call, bit for bit. + # in the same frame: Tr(rho_off) is built in the me_frame, and a + # helicity-restricted matrix element is not Lorentz invariant, + # so the lab-frame calculate_matrix_element is a different + # projection (on a boosted polarised event, by orders of + # magnitude -- which then sets the mass-stage bound). + # Unpolarised runs have no frame boost and keep the + # matrix-element call, bit for bit. me_prod_on = self._onshell_production_norm(production, prod_static) production.me_wgt = me_prod_on @@ -11652,32 +11631,15 @@ def calculate_matrix_element_from_density(self, production, decays, decay_dict, # compute the denominator and then reshuffle the event before # computing the numerator # - # The denominator has to be the same quantity as the numerator, - # in the same frame: the numerator is the contraction of the - # (possibly helicity-restricted) production density built in the - # me_frame, and a *restricted* matrix element is not Lorentz - # invariant, so the lab-frame calculate_matrix_element is a - # different projection. _onshell_production_norm returns exactly - # calculate_matrix_element when there is no frame boost, and the - # trace of the on-shell rho in the frame when there is. Same fix - # as on the sequential mass stage. + # The denominator has to be taken in the same frame as the + # numerator (see _onshell_production_norm), as on the + # sequential mass stage. # - # This call site is NOT polarised-only: an unpolarised - # production that only asks for keep_weight_for_polarization_* - # still switches the frame axis on, and 'unweighting = auto' - # sends an unbraced production here. That case is safe, and - # measured to be so rather than assumed. The denominator is a - # per-production-event constant either way (checked directly: - # over 5414 production events of `p p > z z [QCD]` the value was - # bit-identical across every joint trial of the same event), and - # a constant factor cancels out of an accept/reject -- so the - # accepted decay distribution is unchanged in law. End to end, - # redecaying 150000 production events of one such sample here - # gave f0 = 0.18199 +- 0.00054, f00 0.05857, fTT 0.69460, Ckk - # -0.57691 against 0.18176 +- 0.00054 / 0.05830 / 0.69478 / - # -0.57706 from the pre-fix run on the SAME production events; - # and a seed-matched 45000-event pair decayed both ways cost - # 10.09 against 10.15 accept/reject trials per event. + # Not polarised-only: keep_weight_for_polarization_* or + # 'unweighting = auto' can bring an unbraced production here with + # the frame axis on. That is safe: the denominator is a constant + # per production event, identical across its joint trials, so it + # cancels out of the accept/reject. MEdenom_prod = self._onshell_production_norm(production, prod_static) MEdenom_decay = 1.0 @@ -12095,7 +12057,8 @@ def get_density(self, event, position, allow_hel, ncomb, dimension, else: #in any case, we need tag to differentiate between production and decay tag = getattr(event, '_ms_tag_for_density', None) if tag is None: - tag, _ = event.get_tag_and_order(self._revert_merged or None) + tag, _ = event.get_tag_and_order( + merged_particle=self._revert_merged or None) # Fast path: single-point momentum extraction without permutation @@ -12111,8 +12074,9 @@ def get_density(self, event, position, allow_hel, ncomb, dimension, assert len(all_p) == 1, "Error: get_density can only be called for a single phase-space point" p = all_p[0] # get_pdg below identifies the particles by EXACT momentum equality - # against the event record, so it has to see the lab-frame momenta: - # keep them before the frame boost rewrites every component. + # against the event record, so it has to see the lab-frame momenta. + # _boost_momenta builds new tuples and never touches its input, so + # keeping a reference is enough -- no copy needed. p_lab = p if frame_boost is not None: p = self._boost_momenta(p, frame_boost, rest_leg=frame_rest_leg) diff --git a/tests/unit_tests/madspin/test_madspin.py b/tests/unit_tests/madspin/test_madspin.py index 72e784d193..78ba5d6886 100755 --- a/tests/unit_tests/madspin/test_madspin.py +++ b/tests/unit_tests/madspin/test_madspin.py @@ -5235,6 +5235,11 @@ def __init__(self, lines): self.mg5cmd = self def extract_process(self, line): + if '[' in line: + # like MadSpin's tree-level mg5cmd on an NLO process line + raise self.InvalidCmd('Perturbation order QCD is not among the ' + 'perturbation orders allowed for by the ' + 'loop model') legs = [] initial, final = line.split('>') for state, part in ([(False, p) for p in initial.split()] + @@ -5272,6 +5277,20 @@ def test_partially_polarised_same_pdg(self): self.assertEqual(self.polarization('generate p p > z{0} z'), {23: ((0,), None)}) + def test_nlo_perturbation_bracket_is_ignored(self): + """'p p > z{0} z{0} [QCD]': the tree-level parser refuses the bracket, + which used to leave every NLO polarised production unrestricted (with + only a warning). The bracket says nothing about the braces on the legs, + so the restriction must be the bracket-free line's.""" + for bracket in ('[QCD]', '[virt=QCD]', '[noborn=QCD]'): + self.assertEqual( + self.polarization('generate p p > z{0} z{0} %s' % bracket), + {23: ((0,),)}) + self.assertEqual( + self.polarization('generate p p > w+{0} w+{T} %s' % bracket, + 'add process p p > w+{0} w+{T} j %s' % bracket), + {24: ((0,), (-1, 1))}) + def test_broadcast_survives_extra_subprocesses(self): """The multiplicity of a broadcast pdg does not have to match between subprocesses -- this is the common 'generate X; add process X j' case."""