diff --git a/.github/workflows/acceptancetest.yml b/.github/workflows/acceptancetest.yml index 3c4fa29029..90fe310c63 100644 --- a/.github/workflows/acceptancetest.yml +++ b/.github/workflows/acceptancetest.yml @@ -1389,3 +1389,24 @@ jobs: cd $GITHUB_WORKSPACE ./tests/test_manager.py test_density_mode_vs_standalone_LI1 -pA -t0 -l INFO + + + acceptancetest_105: + # The type of runner that the job will run on + runs-on: ubuntu-24.04 + if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true + + # Steps represent a sequence of tasks that will be executed as part of the job + steps: + # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it + - uses: actions/checkout@v4 + - uses: ./.github/actions/checkout_mg5 + - uses: ./.github/actions/restore_all + + # Polarised fixed-order NLO evaluated in a chosen rest frame. Asserts on + # check_poles, which is the only check sensitive to the frame handling of + # the Born against the virtual. + - name: test one of the test test_polarised_nlo_me_frame + run: | + cd $GITHUB_WORKSPACE + ./tests/test_manager.py test_polarised_nlo_me_frame -pA -t0 -l INFO diff --git a/Template/LO/SubProcesses/genps.f b/Template/LO/SubProcesses/genps.f index 174aae4378..6e7ac3d06d 100644 --- a/Template/LO/SubProcesses/genps.f +++ b/Template/LO/SubProcesses/genps.f @@ -1771,7 +1771,8 @@ subroutine boost_to_frame(P1, frame_id, P2) integer ids(nexternal) integer i,j - logical trivial_boost + integer nsel, isel + logical trivial_boost c uncompress call mapid(frame_id, ids) @@ -1807,12 +1808,41 @@ subroutine boost_to_frame(P1, frame_id, P2) enddo endif enddo - do j=1,3 + do j=1,3 Pboost(j) = -1 * Pboost(j) - enddo + enddo do i=1, nexternal call boostx(p1(0,i), pboost, p2(0,i)) - enddo + enddo + +c A frame built from a single leg puts that leg at rest, and there the +c boost has to be exactly right rather than right to rounding. HELAS +c changes convention at exactly zero: vxxxxx builds the polarisation +c vectors of a massive vector along the z axis when pp.eq.0d0 and along +c the momentum direction otherwise. boostx only reaches p=0 up to the +c rounding of lf (it forms p(i)+q(i)*lf with lf=1 up to the rounding of +c (q(0)-m)+p(0)), so the residual is a few 1d-14 with a noise direction, +c and whether it rounds to zero varies event by event. Every event that +c misses zero gets its longitudinal polarisation vector pointed along +c rounding noise instead of along z. +c This was found through the NLO port, where it breaks the FKS +c subtraction outright: see docs/nlo_polarisation_boost_plan.md, M2, and +c the same fix in Template/NLO/SubProcesses/boost_to_frame.f. +c Only a one-leg selection needs this: with two or more selected legs it +c is their sum that is at rest, no single leg sits on the branch point. + nsel = 0 + isel = 0 + do i=1, nexternal + if (ids(i).eq.1) then + nsel = nsel + 1 + isel = i + endif + enddo + if (nsel.eq.1) then + p2(1,isel) = 0d0 + p2(2,isel) = 0d0 + p2(3,isel) = 0d0 + endif return end diff --git a/Template/NLO/Cards/run_card.dat b/Template/NLO/Cards/run_card.dat index 2b62e62847..2b0c47fe3f 100644 --- a/Template/NLO/Cards/run_card.dat +++ b/Template/NLO/Cards/run_card.dat @@ -66,6 +66,7 @@ %(lpp2)s = lpp2 ! beam 2 type (0 = no PDF) %(ebeam1)s = ebeam1 ! beam 1 energy in GeV %(ebeam2)s = ebeam2 ! beam 2 energy in GeV +$frame #*********************************************************************** # PDF choice: this automatically fixes also alpha_s(MZ) and its evol. * # pdlabel: lhapdf = LHAPDF (installation needed) [1412.7420] * diff --git a/Template/NLO/Source/run.inc b/Template/NLO/Source/run.inc index 0da1f1008d..eceed9d1ec 100644 --- a/Template/NLO/Source/run.inc +++ b/Template/NLO/Source/run.inc @@ -50,8 +50,16 @@ c double precision ebeam(2), xbk(2),q2fact(2) common/to_collider/ ebeam , xbk ,q2fact, lpp c +c Rest frame in which the matrix-elements are evaluated (polarization). +c Bit n of frame_id is set for each leg n listed in me_frame; see +c mapid() and boost_to_me_frame() in boost_to_frame.f. The default +c me_frame=[1,2] gives frame_id=6, which means "no boost". +c + integer frame_id + common/to_frame_me/frame_id +c c BW treatment -c +c double precision bwcutoff common/to_bwcutoff/ bwcutoff c diff --git a/Template/NLO/SubProcesses/add_write_info.f b/Template/NLO/SubProcesses/add_write_info.f index 7ced1cd6c9..2dd29baad2 100644 --- a/Template/NLO/SubProcesses/add_write_info.f +++ b/Template/NLO/SubProcesses/add_write_info.f @@ -275,7 +275,7 @@ subroutine add_write_info(p_born,pp,ybst_til_tolab,iconfig,Hevents if (colour_connections(1,1).lt.0) then ! colour not yet set: Get color flow that is consistent with ! iconfig from Born - call sborn(p_born,wgt1) + call sborn_frame(p_born,wgt1) sumborn=0.d0 do i=1,max_bcol if (icolamp(i,iBornGraph,1)) then diff --git a/Template/NLO/SubProcesses/boost_to_frame.f b/Template/NLO/SubProcesses/boost_to_frame.f new file mode 100644 index 0000000000..915da7426b --- /dev/null +++ b/Template/NLO/SubProcesses/boost_to_frame.f @@ -0,0 +1,645 @@ +c************************************************************************** +c Rest frame in which the matrix-elements are evaluated. +c +c Polarisation of a massive particle is not boost invariant, so for a +c polarised process the matrix-elements must be evaluated in a definite +c frame, chosen by the user through me_frame in the run_card. This file +c provides the kinematics; the selection of which legs define the frame +c is resolved by the caller and passed in as the mask ids(). +c +c IMPORTANT (see docs/nlo_polarisation_boost_plan.md, D3): the frame must +c be recomputed from the momenta of *each* kinematic configuration (Born, +c real, and every counter-event separately). Computing it once from the +c real and reusing it for the reduced Born breaks the cancellation of the +c collinear pole. +c************************************************************************** + + subroutine mapid_frame(id, npart, ids) +c************************************************************************** +c Uncompress frame_id into a 0/1 mask over npart legs. +c Bit n of id is set for each leg n selected by me_frame. This is the +c same encoding as mapid() in the LO cluster.f -- note it is 2**n, not +c 2**(n-1), whatever the stale comment above the LO boost_to_frame says. +c +c input: id compressed frame id (run_card frame_id) +c npart number of legs of this configuration +c output: ids ids(i)=1 if leg i takes part in the frame definition +c************************************************************************** + implicit none + integer id, npart, ids(npart) + integer i + + do i=1,npart + ids(i)=0 + if (btest(id,i)) then + ids(i)=1 + endif + enddo + + return + end + + + subroutine azifact_me_frame(p_born_in, imother, azifact, boosted) +c************************************************************************** +c Collinear azimuthal factor -exp(2 i psi), evaluated in the me_frame. +c +c The frame is rebuilt from the very momenta the matching Born call is +c given, so the spin-correlated Born and the phase multiplying it are +c guaranteed to live in the same frame. That is the whole point: the +c Born is boosted inside sborn_frame, and if the phase were left in the +c partonic c.m. the Q term would be wrong. +c +c The mother and the stored emission direction are boosted together. +c xij_kperp is orthogonal to the mother (k.p = 0 in the frame it was +c built in), and orthogonality is preserved by the boost, so projecting +c the boosted vector onto the helicity basis of the boosted mother +c recovers the azimuth correctly. +c +c boosted=.false. means no frame was requested (or the selection is one +c of the partonic-c.m. spellings). azifact is then left untouched and +c the caller must keep its legacy value, so that runs without a frame +c stay bit-identical. +c +c input: p_born_in(0:3,nexternal-1) Born momenta of this configuration +c imother position of the mother in them +c output: azifact, boosted +c************************************************************************** + implicit none + include 'nexternal.inc' + double precision p_born_in(0:3,nexternal-1) + integer imother + double complex azifact + logical boosted + double precision pboost(0:3), pm(0:3), kp(0:3) + integer ids(nexternal-1) + logical trivial + double precision xij_kperp(0:3) + common/cxij_kperp/xij_kperp + + call get_frame_mask_born(ids) + call get_me_frame_boost(p_born_in, nexternal-1, ids, pboost, + & trivial) + if (trivial) then + boosted=.false. + return + endif + + call boostx(p_born_in(0,imother), pboost, pm) + call boostx(xij_kperp, pboost, kp) + call azifact_from_kperp(pm, kp, azifact) + boosted=.true. + + return + end + + + subroutine azifact_from_kperp(pmother, kperp, azifact) +c************************************************************************** +c Rebuild the collinear limit of /[ij] from the emission direction. +c +c The counterterms need -exp(2 i psi), with psi the azimuth of the +c emission about the mother, measured in the HELAS transverse basis of +c the mother direction. The generator stores that limit as xij_aor, but +c xij_aor is a 0/0 limit of /[ij]: a boost maps exactly-parallel +c null vectors onto exactly-parallel null vectors, so the degeneracy +c survives in every frame and xij_aor can never be recomputed after a +c boost. The azimuthal information has to be carried covariantly +c instead, which is what xij_kperp is for. +c +c Basis, for a mother direction n = (sin(th)cos(ph), sin(th)sin(ph), +c cos(th)): +c e1 = ( cos(th)cos(ph), cos(th)sin(ph), -sin(th) ) +c e2 = ( -sin(ph), cos(ph), 0 ) +c which is the standard helicity basis; for n = +z it is (x,y) and for +c n = -z it is (-x,y). +c +c This reproduces both incoming legs with no idir bookkeeping: +c j_fks=1, n=+z : psi = phi_i -> -exp( 2 i phi_i) +c j_fks=2, n=-z : psi = pi - phi_i -> -exp(-2 i phi_i) +c the doubled angle turning the basis flip into a harmless 2 pi. The +c R_y(pi) rotation the old ISR code applied for j_fks=2 is therefore +c not needed here -- and it must not be reinstated, since after a boost +c it no longer maps the mother onto +z. +c +c input: pmother(0:3) mother 4-momentum (only its direction is used) +c kperp(0:3) emission direction transverse to the mother +c output: azifact -exp(2 i psi) +c************************************************************************** + implicit none + double precision pmother(0:3), kperp(0:3) + double complex azifact + double precision n(3), e1(3), e2(3) + double precision rn, cth, sth, cph, sph, a, b, norm2 + double complex ximag + parameter (ximag=(0d0,1d0)) + integer i + + rn=sqrt(pmother(1)**2+pmother(2)**2+pmother(3)**2) + if (rn.eq.0d0) then + write (*,*) 'ERROR in azifact_from_kperp: mother at rest' + stop 1 + endif + do i=1,3 + n(i)=pmother(i)/rn + enddo + + cth=n(3) + sth=sqrt(max(1d0-cth**2,0d0)) + if (sth.ne.0d0) then + cph=n(1)/sth + sph=n(2)/sth + else +c mother on the beam axis: its azimuth is undefined, take phi=0. +c Same convention as getaziangles(). + cph=1d0 + sph=0d0 + endif + + e1(1)= cth*cph + e1(2)= cth*sph + e1(3)=-sth + e2(1)=-sph + e2(2)= cph + e2(3)= 0d0 + + a=0d0 + b=0d0 + do i=1,3 + a=a+kperp(i)*e1(i) + b=b+kperp(i)*e2(i) + enddo + +c -exp(2 i psi) with psi=atan2(b,a), written without trigonometry so +c that only the component of kperp transverse to the mother enters. + norm2=a**2+b**2 + if (norm2.eq.0d0) then + write (*,*) 'ERROR in azifact_from_kperp: the emission' + write (*,*) 'direction is parallel to the mother, so its' + write (*,*) 'azimuth is undefined.' + stop 1 + endif + azifact=-(dcmplx(a,b)**2)/norm2 + + return + end + + + subroutine sborn_frame(p_in, ans_summed) +c************************************************************************** +c Evaluate the Born in the frame selected by me_frame. +c +c The frame is rebuilt from the momenta actually passed in, so each +c kinematic configuration (the Born of the n-body contribution, and in +c later milestones each counter-event separately) gets its own boost. +c Reusing one frame across configurations breaks the cancellation of the +c collinear pole -- see docs/nlo_polarisation_boost_plan.md, D3. +c +c The boost is applied to a local copy: /pborn/ is read by the event +c record, the cuts and the scales, none of which want the ME frame. +c +c Every caller of SBORN within one event must agree on whether it passes +c boosted or unboosted momenta: SBORN caches its amplitudes against +c (E, p_z) and that cache is shared with sborn_sf, born_hel and +c extra_cnt. Mixing the two silently returns amplitudes from the wrong +c frame. +c************************************************************************** + implicit none + include 'nexternal.inc' + double precision p_in(0:3,nexternal-1) + double precision ans_summed + double precision p_f(0:3,nexternal-1) + integer ids(nexternal-1) + + call get_frame_mask_born(ids) + call boost_to_me_frame(p_in, nexternal-1, ids, p_f) + call sborn(p_f, ans_summed) + + return + end + + + subroutine binothlha_frame(p_in, born_wgt, virt_wgt) +c************************************************************************** +c Virtual in the me_frame rest frame. +c +c BinothLHA hands its argument straight to sloopmatrix_thres, and takes +c born_wgt from the caller rather than recomputing it, so boosting the +c momenta here is enough to put the virtual in the same frame as the +c Born that bornsoftvirtual already boosted. +c +c Two MadLoop settings interact with an imposed frame, both via the +c HELAS branch at exactly zero momentum that a one-leg me_frame sits on +c (see boost_to_me_frame): +c +c - NRotations_DP / NRotations_QP re-evaluate the loop at a rotated +c phase-space point and expect |M|^2 back unchanged. That holds for +c fixed helicities in general, but NOT for a particle at rest, +c whose polarisation axis is the frame z axis rather than its own +c momentum, so a rotation of the momenta does not carry the axis +c with it. Both default to 0, so the rotation test is off unless a +c user turns it on; if they do, a polarised run will see spurious +c instability. Keep them at 0. +c - ImprovePSPoint (default 2) shifts momenta to restore exact +c onshellness, which can move the deliberately-zeroed leg off zero. +c The shift is along z, so the momentum direction stays on the axis +c and eps_L comes back along +-z, differing at most by a sign that +c cancels in |M|^2. Checked to be harmless by check_poles rather +c than assumed. +c************************************************************************** + implicit none + include 'nexternal.inc' + double precision p_in(0:3,nexternal-1) + double precision born_wgt, virt_wgt + double precision p_f(0:3,nexternal-1) + integer ids(nexternal-1) + + call get_frame_mask_born(ids) + call boost_to_me_frame(p_in, nexternal-1, ids, p_f) + call BinothLHA(p_f, born_wgt, virt_wgt) + + return + end + + + subroutine sborn_sf_frame(p_in, m, n, wgt) +c************************************************************************** +c Colour-linked Born in the me_frame rest frame. +c +c SBORN_SF reads the amplitudes SBORN left in the shared cache, so it +c must be handed exactly the momenta the matching sborn_frame call was +c handed, or it silently returns amplitudes from the other frame. +c************************************************************************** + implicit none + include 'nexternal.inc' + double precision p_in(0:3,nexternal-1) + integer m, n + double precision wgt + double precision p_f(0:3,nexternal-1) + integer ids(nexternal-1) + + call get_frame_mask_born(ids) + call boost_to_me_frame(p_in, nexternal-1, ids, p_f) + call sborn_sf(p_f, m, n, wgt) + + return + end + + + subroutine extra_cnt_frame(p_in, icnt, cnts) +c************************************************************************** +c Extra g/a -> q qbar counterterm Born in the me_frame rest frame. +c Same cache as SBORN, so the same rule applies. +c************************************************************************** + implicit none + include 'nexternal.inc' + include 'orders.inc' + double precision p_in(0:3,nexternal-1) + integer icnt + double complex cnts(2,nsplitorders) + double precision p_f(0:3,nexternal-1) + integer ids(nexternal-1) + + call get_frame_mask_born(ids) + call boost_to_me_frame(p_in, nexternal-1, ids, p_f) + call extra_cnt(p_f, icnt, cnts) + + return + end + + + subroutine smatrix_real_frame(p_in, wgt) +c************************************************************************** +c Real-emission matrix element in the me_frame rest frame. +c +c nexternal legs, so the mask is the real one; it is derived from the +c Born mask and i_fks of the current FKS configuration. +c************************************************************************** + implicit none + include 'nexternal.inc' + double precision p_in(0:3,nexternal) + double precision wgt + double precision p_f(0:3,nexternal) + integer ids(nexternal) + integer nFKSprocess + common/c_nFKSprocess/nFKSprocess + + call get_frame_mask_real(nFKSprocess, ids) + call boost_to_me_frame(p_in, nexternal, ids, p_f) + call smatrix_real(p_f, wgt) + + return + end + + + subroutine get_frame_mask_born(ids) +c************************************************************************** +c Mask over the nexternal-1 Born legs selected by me_frame. +c +c me_frame is given in the numbering of the process as the user wrote +c it, which the FKS sort does not preserve; frame_map_born (written at +c generation time into frame_info.inc) maps Born position -> user number. +c +c output: ids(nexternal-1) +c************************************************************************** + implicit none + include 'nexternal.inc' + include 'frame_info.inc' + include 'run.inc' + integer ids(nexternal-1) + integer i + + do i=1,nexternal-1 + ids(i)=0 + if (btest(frame_id, frame_map_born(i))) then + ids(i)=1 + endif + enddo + + return + end + + + subroutine get_frame_mask_real(iFKS, ids) +c************************************************************************** +c Mask over the nexternal legs of the real emission selected by +c me_frame, for FKS configuration iFKS. +c +c The FKS convention fixes the real -> underlying-Born correspondence in +c terms of i_fks alone (see set_pdg in chooser_functions.f): +c real r -> Born r for r < i_fks +c real r -> Born r-1 for r > i_fks +c real i_fks has no Born counterpart (it is the extra parton) +c so no separate table is needed for the reals. +c +c Note on r = j_fks: positionally the Born leg is the *mother* of i_fks +c and j_fks, not the same particle as the real leg j_fks. Selecting a +c splitting parton in me_frame is meaningless for polarisation anyway -- +c the particles whose polarisation is being fixed are spectators, present +c unchanged in both the real and the Born. +c +c output: ids(nexternal) +c************************************************************************** + implicit none + include 'nexternal.inc' + include 'fks_info.inc' + integer iFKS + integer ids(nexternal) + integer ids_born(nexternal-1) + integer r, b + + call get_frame_mask_born(ids_born) + + do r=1,nexternal + ids(r)=0 + if (r.eq.fks_i_D(iFKS)) then +c the extra parton: no Born counterpart, never part of the frame + cycle + elseif (r.lt.fks_i_D(iFKS)) then + b=r + else + b=r-1 + endif + if (b.ge.1 .and. b.le.nexternal-1) then + ids(r)=ids_born(b) + endif + enddo + + return + end + + + subroutine get_me_frame_boost(p, npart, ids, pboost, trivial) +c************************************************************************** +c Build the boost 4-vector that takes the selected system to rest. +c +c pboost is (E, -Pvec) of the sum of the selected momenta, i.e. exactly +c the argument boostx() wants in order to express a momentum given in +c the current frame in the rest frame of that system. +c +c trivial is returned .true. when the boost should be skipped rather +c than applied. Applying an identity boost is not free: it is a +c multiply-and-add through boostx that perturbs the momenta in the last +c bits, which is enough to move an adaptive integration. The LO code +c avoids exactly this, in two places -- the call site skips frame_id=6 +c (auto_dsig_v4.inc) and boost_to_frame() skips the all-final-state +c selection (genps.f) -- and this routine folds both in. +c +c Four cases: +c - nothing selected; +c - the selection is exactly the initial state; +c - the selection is exactly the whole final state; +c - the selected system already has exactly zero 3-momentum. +c +c The middle two are both "the partonic c.m.". Note they are *not* +c identity boosts here the way they are at LO: MadFKS works in a frame +c boosted along z (the two initial momenta carry different energies), +c so honouring them literally would apply a real longitudinal boost. It +c would also not be infrared safe -- a frame defined from the initial +c state jumps between the real emission and the reduced Born, because +c their momentum fractions differ by a finite amount even in the +c singular limit. Skipping is both the safe and the meaningful choice. +c +c input: p(0:3,npart) momenta of this configuration +c ids(npart) 0/1 mask, see mapid_frame +c output: pboost(0:3) boost 4-vector +c trivial .true. if the boost may be skipped +c************************************************************************** + implicit none + integer npart + double precision p(0:3,npart), pboost(0:3) + integer ids(npart) + logical trivial + integer i, j + integer nsel, nini, nfin + double precision m2, pvec2 + double precision betamin + parameter (betamin=1d-4) + include 'nexternal.inc' + + pboost(0:3)=0d0 + nsel=0 + nini=0 + nfin=0 + do i=1,npart + if (ids(i).eq.1) then + nsel=nsel+1 + if (i.le.nincoming) then + nini=nini+1 + else + nfin=nfin+1 + endif + do j=0,3 + pboost(j)=pboost(j)+p(j,i) + enddo + endif + enddo + +c Nothing selected: nothing to do. Do not treat this as an error, it is +c how frame_id=0 (or a mask that misses this configuration) shows up. + if (nsel.eq.0) then + trivial=.true. + return + endif + +c Exactly the initial state, or exactly the whole final state: both name +c the partonic c.m.. Skip rather than boost -- see the header. +c +c Note the second of these is multiplicity-dependent and therefore does +c not mean the same thing for the Born and for the real: me_frame=[3,4] +c on p p > z z is the whole Born final state, but in the real the same +c two legs are only two of three. It was checked that this asymmetry is +c *not* what makes that configuration fail -- removing the skip leaves +c the failure fractions bit-identical, because p_born is handed to us in +c the Born partonic c.m. and so the zero-3-momentum test below catches +c the same case anyway. The skip is kept for parity with LO. + if (nfin.eq.0 .and. nini.eq.nincoming) then + trivial=.true. + return + endif + if (nini.eq.0 .and. nfin.eq.npart-nincoming) then + trivial=.true. + return + endif + + pvec2=pboost(1)**2+pboost(2)**2+pboost(3)**2 + +c Already at rest: skip, so the default costs nothing and changes nothing. + if (pvec2.eq.0d0) then + trivial=.true. + return + endif + +c At rest to well below the precision the matrix elements can use: +c skip too, because below betamin the boost costs more accuracy than +c it delivers. +c +c What it delivers is O(beta): that is how much a boost of velocity +c beta moves a polarised |M|^2. What it costs is a loss of precision +c in the incoming legs. Before the boost they are exactly on the beam +c axis, so p(0)+p(3) of the backward one is exactly zero and the HELAS +c massless spinors -- built from sqrt(p(0)+p(3)) and pt/sqrt(p(0)+p(3)) +c -- take their exact on-axis branch. The boost tilts them by O(beta), +c making p(0)+p(3) = E*beta^2/2, and that is computed as a difference +c of two numbers of size E, so it carries an absolute error ulp(E) and +c a relative error 2*eps/beta^2. Small beta is the dangerous end. +c +c Measured on p p > z{0} z{0} [real=QCD] with me_frame=[3,4], whose +c boost degenerates to the identity as xi->0. Down the soft scan +c p(0)+p(3) of the boosted incoming leg runs +c 3.3d-9, 3.3d-11, 3.4d-13, 0, 0, 0 +c and the deviation of the soft ratio from 1 runs +c 9d-6, 6.5d-4, 2.1d-2, 2d-6, 4d-7, 2d-6, +c i.e. it grows like 1/beta^2 exactly while p(0)+p(3) is a few ulp, +c and goes clean again once it underflows to exactly zero. +c +c Balancing cost against benefit, 30*2*eps/beta^2 = beta (the factor 30 +c is the measured prefactor above), gives beta* = 2.4d-5. betamin is set +c one safety decade above that. It was checked that 1d-5 is not enough +c -- the split-order soft test still fails 0.16 -- and that 1d-4 and +c 1d-3 both give 0.01. +c +c This costs no physics: it fires only where the selected system is at +c rest to one part in 10^4, which in a real integration never happens. +c It is reached only in the artificial deep-soft/collinear scan of +c test_soft_col_limits, and there only for a frame that degenerates to +c the partonic c.m. in the limit being scanned. + if (pvec2 .lt. (betamin*pboost(0))**2) then + trivial=.true. + return + endif + +c The selected system must be timelike for its rest frame to exist. A +c lightlike or spacelike sum means the run_card selection is nonsense +c (e.g. a single massless leg); the LO boost_to_frame has no such guard. + m2=pboost(0)**2-pvec2 + if (m2.le.0d0) then + write (*,*) 'ERROR in get_me_frame_boost: the system' + write (*,*) 'selected by me_frame is not timelike, so its' + write (*,*) 'rest frame does not exist. m^2 = ', m2 + write (*,*) 'Selected legs:', (ids(i),i=1,npart) + stop 1 + endif + + do j=1,3 + pboost(j)=-pboost(j) + enddo + trivial=.false. + + return + end + + + subroutine boost_to_me_frame(p, npart, ids, p_out) +c************************************************************************** +c Express all npart momenta in the rest frame of the selected system. +c Pure boost, no rotation: for fixed external helicities |M|^2 is +c rotation invariant (HELAS builds the polarisation vectors in each +c particle's own momentum-direction basis, so a rotation only produces +c per-particle phases which cancel in the modulus). Orientation does +c matter for the collinear counterterm azimuthal phases, but that is +c handled by the caller, not here. +c +c p and p_out may be the same array. +c************************************************************************** + implicit none + integer npart + double precision p(0:3,npart), p_out(0:3,npart) + integer ids(npart) + double precision pboost(0:3) + logical trivial + integer i, nsel, isel + + call get_me_frame_boost(p, npart, ids, pboost, trivial) + + if (trivial) then + do i=1,npart + p_out(0:3,i)=p(0:3,i) + enddo + return + endif + + do i=1,npart + call boostx(p(0,i), pboost, p_out(0,i)) + enddo + +c A selection made of a single leg puts that leg at rest, and there +c the boost has to be *exactly* right, not right to rounding. +c +c HELAS switches convention at exactly zero: vxxxxx builds the +c polarisation vectors of a massive vector along the z axis when +c pp.eq.0d0, and along the momentum direction otherwise. boostx +c reaches p=0 only up to a relative rounding of order 1d-16 -- it +c forms p(i)+q(i)*lf with lf=1 up to the rounding of (q(0)-m)+p(0), +c so the residual is a few 1d-14 in absolute value and its direction +c is pure noise. When that residual happens not to round to zero, +c eps_L points along the noise instead of along z and |M|^2 changes +c by orders of magnitude. +c +c Whether it rounds to zero is decided independently for the real +c and for the reduced Born, so the two are then evaluated with +c different quantisation axes and the subtraction stops cancelling. +c Measured on p p > z{0} z{0} [real=QCD] with me_frame=[3]: the +c soft/collinear ratio is 1 to 1d-7 whenever the residual is exactly +c zero and 2.2d-3 whenever it is 3d-14, with no intermediate values. +c +c Imposing the defining property of the frame removes the choice. +c Only nsel=1 needs it: with two or more selected legs it is their +c *sum* that is at rest, no individual leg sits at the branch point, +c and the residual of the sum never reaches HELAS. + nsel=0 + isel=0 + do i=1,npart + if (ids(i).eq.1) then + nsel=nsel+1 + isel=i + endif + enddo + if (nsel.eq.1) then + p_out(1,isel)=0d0 + p_out(2,isel)=0d0 + p_out(3,isel)=0d0 + endif + + return + end diff --git a/Template/NLO/SubProcesses/check_poles.f b/Template/NLO/SubProcesses/check_poles.f index c3a2e82ad5..dc8e7296c6 100644 --- a/Template/NLO/SubProcesses/check_poles.f +++ b/Template/NLO/SubProcesses/check_poles.f @@ -221,11 +221,16 @@ Program DRIVER enddo CALL UPDATE_AS_PARAM() - call sborn(p_born, born) + call sborn_frame(p_born, born) ! extra initialisation calls: skip the first point ! as well as any other points which is used for initialization ! (according to the return code) - call BinothLHA(p_born, born, virt_wgt) +c The virtual must be evaluated in the same frame as the Born above: the +c poles are proportional to the Born, so a frame mismatch here shows up as +c a per-point constant ratio between the MadFKS and the OLP poles and the +c check fails for every point. Go through binothlha_frame, exactly as the +c integration does at fks_singular.f (bornsoftvirtual). + call binothlha_frame(p_born, born, virt_wgt) if (npointsChecked.eq.0) then if (mod(ret_code_ml,100)/10.eq.3 .or. & mod(ret_code_ml,100)/10.eq.4) then diff --git a/Template/NLO/SubProcesses/fks_singular.f b/Template/NLO/SubProcesses/fks_singular.f index 4bdb6ab524..060eeaf30e 100644 --- a/Template/NLO/SubProcesses/fks_singular.f +++ b/Template/NLO/SubProcesses/fks_singular.f @@ -34,7 +34,9 @@ subroutine compute_born call cpu_time(tBefore) if (f_b.eq.0d0) return if (xi_i_hat_ev*xiimax_cnt(0) .gt. xiBSVcut_used) return - call sborn(p_born,wgt_c) +c Born of the n-body contribution, evaluated in the me_frame rest frame +c (identity unless the run_card asks for a frame; see boost_to_frame.f). + call sborn_frame(p_born,wgt_c) do iamp=1, amp_split_size if (amp_split(iamp).eq.0d0) cycle call amp_split_pos_to_orders(iamp, orders) @@ -127,7 +129,7 @@ subroutine compute_6to5flav_cnt() endif ! compute the born - call sborn(p_born,wgtborn) + call sborn_frame(p_born,wgtborn) alphas = g**2/4d0/pi do iamp = 1, amp_split_size if (amp_split(iamp).eq.0d0) cycle @@ -216,7 +218,7 @@ subroutine compute_ewsudakov ! sud_mod = 0 do sud_mod = 0,1 - call sborn(p_born,wgt_c) + call sborn_frame(p_born,wgt_c) call sudakov_wrapper(p_born) do iamp=1, amp_split_size if (amp_split_ewsud_lsc(iamp).eq.0d0.and. @@ -350,7 +352,7 @@ subroutine compute_alpha_cnt() ! processes automatically generated by MG5aMC jsign = sign(1,alphascheme) ! compute the born - call sborn(p_born,wgtborn) + call sborn_frame(p_born,wgtborn) ! assumes alpha(MZ) for model, MSbar for PDFs ! the number of flavours depends on mur. ! here we will treat all leptons as massless, @@ -1500,7 +1502,10 @@ subroutine include_multichannel_enhance(imode) c Compute the multi-channel enhancement factor 'enhance'. enhance=1.d0 if (p_born(0,1).gt.0d0) then - call sborn(p_born,wgt_c) +c Must use the same frame as every other SBORN call in this event: the +c amplitudes are cached against (E,p_z) and shared. The multi-channel +c weights stay a partition of unity in any frame, so this is free. + call sborn_frame(p_born,wgt_c) elseif(p_born(0,1).lt.0d0)then enhance=0d0 endif @@ -1543,7 +1548,9 @@ subroutine include_multichannel_enhance(imode) pas(0:3,nexternal)=0d0 pas(0:3,1:nexternal-1)=p_born_used(0:3,1:nexternal-1) call set_alphas(pas) - call sborn(p_born_used,wgt_c) +c Own frame, rebuilt from p_born_used: this is a different kinematic +c configuration. Cache-isolated by the calculatedBorn resets around it. + call sborn_frame(p_born_used,wgt_c) call set_alphas(p_ev) calculatedBorn=.false. elseif(p_born_used(0,1).lt.0d0)then @@ -4698,7 +4705,7 @@ subroutine sreal(pp,xi_i_fks,y_ij_fks,wgt) amp_split(1:amp_split_size) = 0d0 endif else - call smatrix_real(pp,wgt) + call smatrix_real_frame(pp,wgt) wgt=wgt*xi_i_fks**2*(1d0-y_ij_fks) amp_split(1:amp_split_size) = amp_split(1:amp_split_size)*xi_i_fks**2*(1d0-y_ij_fks) endif @@ -4745,6 +4752,7 @@ subroutine sborncol_fsr(p,xi_i_fks,y_ij_fks,wgt) TYPE(ALOHA) W1,W2,W3,W4 double complex Wij_angle,Wij_recta double complex azifact + logical me_boosted c Particle types (=color/charges) of i_fks, j_fks and fks_mother integer i_type,j_type,m_type @@ -4784,9 +4792,9 @@ subroutine sborncol_fsr(p,xi_i_fks,y_ij_fks,wgt) E_i_fks = p(0,i_fks) z = 1d0 - E_i_fks/(E_i_fks+E_j_fks) t = z * shat/4d0 - call sborn(p_born,wgt_born) + call sborn_frame(p_born,wgt_born) if (iextra_cnt.gt.0) - 1 call extra_cnt(p_born, iextra_cnt, ans_extra_cnt) + 1 call extra_cnt_frame(p_born, iextra_cnt, ans_extra_cnt) call AP_reduced(j_type,i_type,ch_j,ch_i,t,z,ap) call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,t,z,Q) wgt=0d0 @@ -4796,12 +4804,12 @@ subroutine sborncol_fsr(p,xi_i_fks,y_ij_fks,wgt) C check if any extra_cnt is needed if (iextra_cnt.gt.0) then if (iord.eq.isplitorder_born) then - call sborn(p_born,wgt_born) + call sborn_frame(p_born,wgt_born) wgt1(1) = ans_cnt(1,iord) wgt1(2) = ans_cnt(2,iord) elseif (iord.eq.isplitorder_cnt) then ! this is the contribution from the extra cnt - call extra_cnt(p_born, iextra_cnt, ans_extra_cnt) + call extra_cnt_frame(p_born, iextra_cnt, ans_extra_cnt) wgt1(1) = ans_extra_cnt(1,iord) wgt1(2) = ans_extra_cnt(2,iord) else @@ -4809,7 +4817,7 @@ subroutine sborncol_fsr(p,xi_i_fks,y_ij_fks,wgt) stop endif else - call sborn(p_born,wgt_born) + call sborn_frame(p_born,wgt_born) wgt1(1) = ans_cnt(1,iord) wgt1(2) = ans_cnt(2,iord) endif @@ -4820,6 +4828,17 @@ subroutine sborncol_fsr(p,xi_i_fks,y_ij_fks,wgt) wgt1(2)=0d0 elseif (m_type.eq.8.or.ch_m.eq.0d0) then c Insert /[ij] which is not included by sborn() + imother_fks=min(i_fks,j_fks) + call azifact_me_frame(p_born,imother_fks,azifact,me_boosted) + if (me_boosted) then +c Same as the ISR case, up to the conjugation that distinguishes a timelike +c from a spacelike splitting: psi is measured in the helicity basis of the +c boosted mother, so the getaziangles factor of the legacy branch below is +c already accounted for and must not be applied again. + wgt1(2) = -azifact * wgt1(2) + amp_split_cnt(1:amp_split_size,2,iord) = -azifact + $ *amp_split_cnt(1:amp_split_size,2,iord) + else if (1d0-y_ij_fks.lt.vtiny)then azifact=xij_aor else @@ -4840,7 +4859,6 @@ subroutine sborncol_fsr(p,xi_i_fks,y_ij_fks,wgt) azifact=Wij_angle/Wij_recta endif c Insert the extra factor due to Madgraph convention for polarization vectors - imother_fks=min(i_fks,j_fks) call getaziangles(p_born(0,imother_fks), # cphi_mother,sphi_mother) wgt1(2) = -(cphi_mother-ximag*sphi_mother)**2 * @@ -4848,6 +4866,7 @@ subroutine sborncol_fsr(p,xi_i_fks,y_ij_fks,wgt) amp_split_cnt(1:amp_split_size,2,iord) = -(cphi_mother-ximag $ *sphi_mother)**2 *amp_split_cnt(1:amp_split_size,2 $ ,iord) * azifact + endif else write(*,*) 'FATAL ERROR in sborncol_fsr',i_type,j_type,i_fks $ ,j_fks @@ -4920,6 +4939,7 @@ subroutine sborncol_isr(p,xi_i_fks,y_ij_fks,wgt) TYPE(ALOHA) W1,W2,W3,W4 double complex Wij_angle,Wij_recta double complex azifact + logical me_boosted double precision zero,vtiny parameter (zero=0d0) @@ -4978,18 +4998,18 @@ subroutine sborncol_isr(p,xi_i_fks,y_ij_fks,wgt) if (iextra_cnt.gt.0) then if (iord.eq.isplitorder_born) then ! this is the contribution from the born ME - call sborn(p_born_used,wgt_born) + call sborn_frame(p_born_used,wgt_born) wgt1(1:2) = ans_cnt(1:2,iord) else if (iord.eq.isplitorder_cnt) then ! this is the contribution from the extra cnt - call extra_cnt(p_born_used, iextra_cnt, ans_extra_cnt) + call extra_cnt_frame(p_born_used, iextra_cnt, ans_extra_cnt) wgt1(1:2) = ans_extra_cnt(1:2,iord) else write(*,*) 'ERROR in sborncol_isr', iord stop endif else - call sborn(p_born_used,wgt_born) + call sborn_frame(p_born_used,wgt_born) wgt1(1:2) = ans_cnt(1:2,iord) endif amp_split_cnt_local(1:amp_split_size,1,iord)= @@ -5003,6 +5023,19 @@ subroutine sborncol_isr(p,xi_i_fks,y_ij_fks,wgt) amp_split_cnt_local(1:amp_split_size,2,iord)=dcmplx(0d0,0d0) else c Insert /[ij] which is not included by sborn() + call azifact_me_frame(p_born_used,j_fks,azifact,me_boosted) + if (me_boosted) then +c In the me_frame the mother is no longer on the beam axis, so neither the +c cphi_mother=1 shortcut nor the R_y(pi) flip of the legacy branch below is +c valid any more -- both only ever compensated for the mother lying along +c +-z. azifact_me_frame returns -exp(2 i psi) with psi measured in the +c helicity basis of the boosted mother, which carries both effects by +c construction. Do not reinstate either here. + wgt1(2) = -dconjg(azifact) * wgt1(2) + amp_split_cnt_local(1:amp_split_size,2,iord) = + $ -dconjg(azifact) + $ *amp_split_cnt_local(1:amp_split_size,2,iord) + else if (1d0-y_ij_fks.lt.vtiny)then azifact=xij_aor else @@ -5011,17 +5044,17 @@ subroutine sborncol_isr(p,xi_i_fks,y_ij_fks,wgt) pj(i)=p(i,j_fks) enddo if(j_fks.eq.2 .and. nincoming.eq.2)then -c Rotation according to innerpin.m. Use rotate_invar() if a more +c Rotation according to innerpin.m. Use rotate_invar() if a more c general rotation is needed pi(1)=-pi(1) pi(3)=-pi(3) pj(1)=-pj(1) pj(3)=-pj(3) endif - CALL IXXXXX(pi ,ZERO ,+1,+1,1,W1) - CALL OXXXXX(pj ,ZERO ,-1,+1,1,W2) - CALL IXXXXX(pi ,ZERO ,-1,+1,1,W3) - CALL OXXXXX(pj ,ZERO ,+1,+1,1,W4) + CALL IXXXXX(pi ,ZERO ,+1,+1,1,W1) + CALL OXXXXX(pj ,ZERO ,-1,+1,1,W2) + CALL IXXXXX(pi ,ZERO ,-1,+1,1,W3) + CALL OXXXXX(pj ,ZERO ,+1,+1,1,W4) Wij_angle=(0d0,0d0) Wij_recta=(0d0,0d0) do i=1,4 @@ -5039,6 +5072,7 @@ subroutine sborncol_isr(p,xi_i_fks,y_ij_fks,wgt) $ +ximag*sphi_mother)**2 $ *amp_split_cnt_local(1:amp_split_size,2,iord) * $ dconjg(azifact) + endif endif if (iord.eq.qcd_pos) then wgt=wgt+dble(wgt1(1)*ap(1)+wgt1(2)*Q(1)) @@ -5726,7 +5760,7 @@ subroutine sbornsoft(pp,xi_i_fks,y_ij_fks,wgt) c should always be done before calling the color-correlated Borns, c because of the caching of the diagrams. c - call sborn(p_born(0,1),wgt1) + call sborn_frame(p_born(0,1),wgt1) c C Reset the amp_split array amp_split(1:amp_split_size) = 0d0 @@ -5739,7 +5773,7 @@ subroutine sbornsoft(pp,xi_i_fks,y_ij_fks,wgt) if ((m.ne.n .or. (m.eq.n .and. pmass(m).ne.ZERO)) .and. & n.ne.i_fks.and.m.ne.i_fks) then C wgt includes the gs/w^2 - call sborn_sf(p_born,m,n,wgt) + call sborn_sf_frame(p_born,m,n,wgt) if (wgt.ne.0d0) then call eikonal_reduced(pp,m,n,i_fks,j_fks, # xi_i_fks,y_ij_fks,eik) @@ -6023,12 +6057,12 @@ subroutine sreal_deg(p,xi_i_fks,y_ij_fks, if (iextra_cnt.gt.0) then if (iord.eq.isplitorder_born) then ! this is the contribution from the born ME - call sborn(p_born_used,wgt_born) + call sborn_frame(p_born_used,wgt_born) wgt1(1) = ans_cnt(1,iord) wgt1(2) = ans_cnt(2,iord) else if (iord.eq.isplitorder_cnt) then ! this is the contribution from the extra cnt - call extra_cnt(p_born_used,iextra_cnt,ans_extra_cnt) + call extra_cnt_frame(p_born_used,iextra_cnt,ans_extra_cnt) wgt1(1) = ans_extra_cnt(1,iord) wgt1(2) = ans_extra_cnt(2,iord) else @@ -6036,7 +6070,7 @@ subroutine sreal_deg(p,xi_i_fks,y_ij_fks, stop endif else - call sborn(p_born_used,wgt_born) + call sborn_frame(p_born_used,wgt_born) wgt1(1) = ans_cnt(1,iord) wgt1(2) = ans_cnt(2,iord) endif @@ -6891,7 +6925,9 @@ subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt) stop endif - call sborn(p_born,wgt1) +c Born of the n-body contribution, evaluated in the me_frame rest frame +c (identity unless the run_card asks for a frame; see boost_to_frame.f). + call sborn_frame(p_born,wgt1) c Born contribution: bsv_wgt=wgt1 @@ -7017,7 +7053,7 @@ subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt) C setup the fks i/j info call fks_inc_chooser() C the following call to born is to setup the goodhel(nfksprocess) - call sborn(p_born,wgt1) + call sborn_frame(p_born,wgt1) contr=0d0 do i=1,fks_j_from_i(i_fks,0) do j=1,i @@ -7029,7 +7065,7 @@ subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt) c *always* a call to sborn(p_born,wgt) just before. This is okay, c because there is a call above in this subroutine C wgt includes the gs/w^2 - call sborn_sf(p_born,m,n,wgt) + call sborn_sf_frame(p_born,m,n,wgt) if (wgt.ne.0d0) then call eikonal_Ireg(p,m,n,xicut_used,eikIreg) contr=contr+wgt*eikIreg @@ -7060,7 +7096,7 @@ subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt) c convert to Binoth Les Houches Accord standards virt_wgt=0d0 - call sborn(p_born, wgt1) + call sborn_frame(p_born, wgt1) ! use the amp_split_cnt as the born to approximate the virtual ! check which one of the two (QCD, QED) is !=0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC @@ -7084,7 +7120,7 @@ subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt) if ((ran2().le.virtual_fraction(ichan) .and. $ abrv(1:3).ne.'nov').or.abrv(1:4).eq.'virt') then call cpu_time(tBefore) - Call BinothLHA(p_born,born_wgt,virt_wgt) + Call binothlha_frame(p_born,born_wgt,virt_wgt) do iamp=1,amp_split_size amp_split_virt(iamp)=amp_split_finite_ML(iamp) enddo @@ -7144,7 +7180,7 @@ subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt) c eq.(MadFKS.C.13) if(abrv.ne.'virt')then ! this is to update the amp_split array - call sborn(p_born,wgt1) + call sborn_frame(p_born,wgt1) bsv_wgt_mufoqes=0d0 do iamp=1,amp_split_size if (dble(amp_split_cnt(iamp,1,qcd_pos)).eq.0d0) cycle @@ -7187,7 +7223,7 @@ subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt) amp_split_wgtwnstmpmur(1:amp_split_size)=0d0 if(abrv.ne.'born' .and. abrv.ne.'grid')then - call sborn(p_born,wgt1) + call sborn_frame(p_born,wgt1) if(abrv(1:2).eq.'vi')then wgtwnstmpmur=0.d0 else @@ -7270,7 +7306,7 @@ subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt) endif if (ComputePoles) then - call sborn(p_born,wgt1) + call sborn_frame(p_born,wgt1) print*," " write(*,123)((p(i,j),i=0,3),j=1,nexternal) @@ -7619,7 +7655,7 @@ subroutine getpoles(p,xmu2,double,single,fksprefact) enddo aso2pi=g**2/(8d0*pi**2) aeo2pi=dble(gal(1))**2/(8d0*pi**2) - call sborn(p_born,wgt1) + call sborn_frame(p_born,wgt1) c QCD Born terms contr1 = 0d0 contr2 = 0d0 @@ -7697,7 +7733,7 @@ subroutine getpoles(p,xmu2,double,single,fksprefact) C setup the fks i/j info call fks_inc_chooser() C the following call to born is to setup the goodhel(nfksprocess) - call sborn(p_born,wgt1) + call sborn_frame(p_born,wgt1) contr1=0d0 do i=1,fks_j_from_i(i_fks,0) @@ -7706,7 +7742,7 @@ subroutine getpoles(p,xmu2,double,single,fksprefact) n=fks_j_from_i(i_fks,j) if( m.ne.n .and. n.ne.i_fks .and. m.ne.i_fks )then C wgt includes the gs/w^2 factor - call sborn_sf(p_born,m,n,wgt) + call sborn_sf_frame(p_born,m,n,wgt) c The factor -2 compensate for that missing in sborn_sf wgt=-2d0*wgt if(wgt.ne.0.d0)then diff --git a/Template/NLO/SubProcesses/genps_fks.f b/Template/NLO/SubProcesses/genps_fks.f index d612df42cd..72c97fd3a6 100644 --- a/Template/NLO/SubProcesses/genps_fks.f +++ b/Template/NLO/SubProcesses/genps_fks.f @@ -1100,6 +1100,8 @@ subroutine generate_FKS_kinematics(x,ndim,xjac0,xpswgt0, double complex xij_aor common/cxij_aor/xij_aor + double precision xij_kperp(0:3) + common/cxij_kperp/xij_kperp integer i_fks,j_fks common/fks_indices/i_fks,j_fks @@ -1175,6 +1177,7 @@ subroutine generate_FKS_kinematics(x,ndim,xjac0,xpswgt0, c if collinear counterevent will not be generated, the following c quantity will stay zero if (.not.only_event_phsp) xij_aor=(0.d0,0.d0) + if (.not.only_event_phsp) xij_kperp(0:3)=0d0 c c These will correspond to the vegas x's for the FKS variables xi_i, c y_ij and phi_i (changing this also requires changing folding parameters) @@ -1389,6 +1392,8 @@ subroutine generate_noevpr_kinematics(x,ndim,xjac0,xpswgt0, double complex xij_aor common/cxij_aor/xij_aor + double precision xij_kperp(0:3) + common/cxij_kperp/xij_kperp integer i_fks,j_fks common/fks_indices/i_fks,j_fks @@ -1474,6 +1479,7 @@ subroutine generate_noevpr_kinematics(x,ndim,xjac0,xpswgt0, c if collinear counterevent will not be generated, the following c quantity will stay zero if (.not.only_event_phsp) xij_aor=(0.d0,0.d0) + if (.not.only_event_phsp) xij_kperp(0:3)=0d0 ! if we do not do event projection, we first generate y/xi FKS ! and p_i_fks, then the other momenta @@ -1908,6 +1914,8 @@ subroutine generate_momenta_initial_noevpr(icountevts,i_fks,j_fks, common/cgenps_fks/veckn_ev,veckbarn_ev,xp0jfks double complex xij_aor common/cxij_aor/xij_aor + double precision xij_kperp(0:3) + common/cxij_kperp/xij_kperp logical softtest,colltest common/sctests/softtest,colltest double precision xi_i_fks_fix,y_ij_fks_fix @@ -2082,6 +2090,13 @@ subroutine generate_momenta_initial_noevpr(icountevts,i_fks,j_fks, & (icountevts.eq.1.and.xij_aor.eq.0) )then resAoR0=-exp( 2*idir*ximag*phi_i_fks ) xij_aor=resAoR0 +c Emission direction transverse to the beam axis. xij_aor above is a 0/0 +c limit of /[ij], so it cannot be recomputed in another frame; this +c vector carries the same azimuthal information covariantly and can be. + xij_kperp(0)=0d0 + xij_kperp(1)=cosphi_i_fks + xij_kperp(2)=sinphi_i_fks + xij_kperp(3)=0d0 endif c c Phase-space factor for (xii,yij,phii) * (tau,ycm) @@ -2119,6 +2134,8 @@ subroutine generate_momenta_massless_final(icountevts,i_fks,j_fks common/cgenps_fks/veckn_ev,veckbarn_ev,xp0jfks double complex xij_aor common/cxij_aor/xij_aor + double precision xij_kperp(0:3) + common/cxij_kperp/xij_kperp logical softtest,colltest common/sctests/softtest,colltest double precision xi_i_fks_fix,y_ij_fks_fix @@ -2349,6 +2366,19 @@ subroutine generate_momenta_massless_final(icountevts,i_fks,j_fks c$$$ & exp( 2*ximag*(phi_mother_fks+phi_i_fks) ) c$$$ xij_aor=resAoR0+resAoR5*sqrt(1-y_ij_fks) xij_aor=resAoR0 +c Emission direction transverse to the mother, carried covariantly so the +c azimuth can be rebuilt in another frame (xij_aor cannot be: it is a 0/0 +c limit of /[ij]). The emission is generated about a +z mother, so its +c transverse direction is (cos(phi_i),sin(phi_i),0); rotate_invar maps that +c onto the mother's actual orientation, and since it takes x,y onto the +c helicity basis vectors e1,e2 the azimuth of the result is again phi_i. + xij_kperp(0)=0d0 + xij_kperp(1)=cosphi_i_fks + xij_kperp(2)=sinphi_i_fks + xij_kperp(3)=0d0 + call rotate_invar(xij_kperp,xij_kperp, + & costh_mother_fks,sinth_mother_fks, + & cosphi_mother_fks,sinphi_mother_fks) endif c c Phase-space factor for (xii,yij,phii) @@ -2714,6 +2744,8 @@ subroutine generate_momenta_initial(icountevts,i_fks,j_fks, common/cgenps_fks/veckn_ev,veckbarn_ev,xp0jfks double complex xij_aor common/cxij_aor/xij_aor + double precision xij_kperp(0:3) + common/cxij_kperp/xij_kperp logical softtest,colltest common/sctests/softtest,colltest double precision xi_i_fks_fix,y_ij_fks_fix @@ -3098,6 +3130,13 @@ subroutine generate_momenta_initial(icountevts,i_fks,j_fks, & (icountevts.eq.1.and.xij_aor.eq.0) )then resAoR0=-exp( 2*idir*ximag*phi_i_fks ) xij_aor=resAoR0 +c Emission direction transverse to the beam axis. xij_aor above is a 0/0 +c limit of /[ij], so it cannot be recomputed in another frame; this +c vector carries the same azimuthal information covariantly and can be. + xij_kperp(0)=0d0 + xij_kperp(1)=cosphi_i_fks + xij_kperp(2)=sinphi_i_fks + xij_kperp(3)=0d0 endif c c Phase-space factor for (xii,yij,phii) * (tau,ycm) diff --git a/Template/NLO/SubProcesses/makefile_fks_dir b/Template/NLO/SubProcesses/makefile_fks_dir index 8325db8584..e3328dca98 100644 --- a/Template/NLO/SubProcesses/makefile_fks_dir +++ b/Template/NLO/SubProcesses/makefile_fks_dir @@ -51,6 +51,7 @@ FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) \ rescale_alpha_tagged.o \ fks_Sij.o $(fastjetfortran_madfks) fks_singular.o \ montecarlocounter.o reweight_xsec.o boostwdir2.o \ + boost_to_frame.o \ cluster.o splitorders_stuff.o \ iproc_map.o sudakov.o \ MC_integer.o $(reweight_xsec_events_pdf_dummy) \ diff --git a/Template/NLO/SubProcesses/montecarlocounter.f b/Template/NLO/SubProcesses/montecarlocounter.f index cda984b8fe..e6c3b9900d 100644 --- a/Template/NLO/SubProcesses/montecarlocounter.f +++ b/Template/NLO/SubProcesses/montecarlocounter.f @@ -2645,7 +2645,7 @@ subroutine assign_emsca_and_flow_statistical(xmcxsec,xmcxsec2 endif else ! use the born-bars - call sborn(p_born,dummy) + call sborn_frame(p_born,dummy) wgt1=0.d0 do i=1,max_bcol wgt1=wgt1+jamp2(i) @@ -2923,12 +2923,12 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde) p_born_rot(3,i)=-p_born(3,i) enddo calculatedBorn=.false. - call sborn(p_born_rot,wgt_born) - if (iextra_cnt.gt.0) call extra_cnt(p_born_rot, iextra_cnt, ans_extra_cnt) + call sborn_frame(p_born_rot,wgt_born) + if (iextra_cnt.gt.0) call extra_cnt_frame(p_born_rot, iextra_cnt, ans_extra_cnt) calculatedBorn=.false. else - call sborn(p_born,wgt_born) - if (iextra_cnt.gt.0) call extra_cnt(p_born, iextra_cnt, ans_extra_cnt) + call sborn_frame(p_born,wgt_born) + if (iextra_cnt.gt.0) call extra_cnt_frame(p_born, iextra_cnt, ans_extra_cnt) endif do iord = 1, nsplitorders diff --git a/Template/NLO/SubProcesses/montecarlocounter_alt.f b/Template/NLO/SubProcesses/montecarlocounter_alt.f index 5f198d2096..fe4e8fb552 100644 --- a/Template/NLO/SubProcesses/montecarlocounter_alt.f +++ b/Template/NLO/SubProcesses/montecarlocounter_alt.f @@ -1372,7 +1372,7 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde) c C BORN - call sborn(p_born,wgt_born) + call sborn_frame(p_born,wgt_born) do iord = 1, nsplitorders if (.not.split_type(iord).or.(iord.ne.qed_pos.and.iord.ne.qcd_pos)) cycle born(iord)=dble(ans_cnt(1,iord)) @@ -1405,12 +1405,12 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde) p_born_rot(3,i)=-p_born(3,i) enddo calculatedBorn=.false. - call sborn(p_born_rot,wgt_born) - if (iextra_cnt.gt.0) call extra_cnt(p_born_rot, iextra_cnt, ans_extra_cnt) + call sborn_frame(p_born_rot,wgt_born) + if (iextra_cnt.gt.0) call extra_cnt_frame(p_born_rot, iextra_cnt, ans_extra_cnt) calculatedBorn=.false. else - call sborn(p_born,wgt_born) - if (iextra_cnt.gt.0) call extra_cnt(p_born, iextra_cnt, ans_extra_cnt) + call sborn_frame(p_born,wgt_born) + if (iextra_cnt.gt.0) call extra_cnt_frame(p_born, iextra_cnt, ans_extra_cnt) endif do iord = 1, nsplitorders diff --git a/Template/NLO/SubProcesses/test_soft_col_limits.f b/Template/NLO/SubProcesses/test_soft_col_limits.f index 31f03c22ed..bb1b110701 100644 --- a/Template/NLO/SubProcesses/test_soft_col_limits.f +++ b/Template/NLO/SubProcesses/test_soft_col_limits.f @@ -275,7 +275,7 @@ program test_soft_col_limits $ /' Cannot perform ME tests properly for config',iconfig cycle endif - call sborn(p_born,wgt1) + call sborn_frame(p_born,wgt1) write (*,*) '' write (*,*) '' diff --git a/aloha/aloha_writers.py b/aloha/aloha_writers.py index c9ec6da6c6..839d93bbf5 100755 --- a/aloha/aloha_writers.py +++ b/aloha/aloha_writers.py @@ -1533,8 +1533,15 @@ def myHash(target_string): addon = '' else: name = short_name - if unknown_tag: + if unknown_tag and outgoing: addon += '%(propa)s' + elif unknown_tag: + # for an amplitude (outgoing == 0) the caller fills 'propa' with '' and + # puts the FLV_Coupling flag ('M') into 'tags' instead -- see + # HelasAmplitude.get_helas_call_dict. Same convention as the FFV1_2 + # scheme above, otherwise the call site emits FFV2_FFS1_0 while ALOHA + # writes FFV2_FFS1M_0. + addon += '%(tags)s' # if outgoing is not None: # return '_'.join((name,) + tuple(other_names)) + addon + '_%s' % outgoing diff --git a/docs/nlo_polarisation_boost_plan.md b/docs/nlo_polarisation_boost_plan.md new file mode 100644 index 0000000000..8c4f6c3496 --- /dev/null +++ b/docs/nlo_polarisation_boost_plan.md @@ -0,0 +1,1260 @@ +# Polarised cross-sections at NLO: frame-boost implementation plan + +This document records the assessment of the existing (dead) polarised-NLO code +and the plan to make `p p > z{0} z{0} j [QCD]` work. + +**It works.** `check_poles` cancels 20 of 20 points in all 12 P dirs and +`calculate_xsect NLO` gives 4.779e-01 +- 3.3e-03 pb (`me_frame=[3]`, nn23lo1, +fixed scales 91.188, ptj 30, etaj 4.0). See M3 for the gate results. + +Status: + +| milestone | state | +|---|---| +| M0 plumbing | **done** — run_card option, frame bookkeeping, boost routine, all inert | +| M1 `[LOonly=QCD]` | **done** — Born boosted, validated against LO madevent | +| M2 step 0 | **done** — ISR and FSR emission azimuth carried covariantly | +| M2 rest | **done** — reals + counterterms boosted, azimuthal wiring in, `test_soft_col_limits` passes, `[real=QCD]` enabled | +| M3 `[QCD]` | **done** — `check_poles` cancels 20/20 in every P dir, `calculate_xsect NLO` runs end to end, `[QCD]` enabled | +| M4 | **done** — guard open, `help_polarization` rewritten, acceptance test written, run and wired into CI | + +`[LOonly=QCD]`, `[real=QCD]`, `[QCD]` and `[virt=…]` all accept a polarised +massive particle today. The first three get their frame from `me_frame` in the +run_card; `[virt=…]` is standalone MadLoop, where the user supplies the +phase-space point and therefore picks the frame themselves (see M3). + +Fixed along the way, and **live at LO on its own**: the `me_frame` boost left +the selected leg at rest only to rounding, and HELAS switches quantisation axis +at exactly zero. See M2 step 2. It shifts the shipped LO polarised +cross-section by 1.8% (5 sigma). + +## 1. Assessment of the existing implementation + +This section records the state **before** this work started; the guard it +quotes has since been rewritten (M3/M4). + +### The guard + +`madgraph/interface/madgraph_interface.py:1240` + +```python +if '[' in process and '{' in process: + if 'noborn' in process or 'sqrvirt' in process: valid = True + else: raise InvalidCmd('Polarization restriction can not be used for NLO processes') + # below are the check when [QCD] will be valid for computation + # if order.strip().lower() != 'qcd': + # raise InvalidCmd('...generic NLO computations') +``` + +The commented-out block at `:1251` is the only authored trace of an intended +`[QCD]` path. The surviving `check(p)` at `:1254` rejects colour-charged **and +massive** particles, so even `p p > z{0} z{0} [noborn=QCD]` is refused today. +The live NLO polarisation path is massless-colourless loop-induced only, which +is exported through the **LO** madevent template and therefore inherits +`me_frame` for free. That is why an NLO boost was never needed. + +### Generation is alive and correct + +`FKSLeg` inherits `Leg`, so `polarization` survives; `fks_common.py:747` +(`to_leg`) copies it explicitly. With the guard bypassed in-memory, +`p p > z{0} z{0} j [QCD]` generates cleanly — 8 born processes, with `{0}` on +the Z legs of every born **and** every real: + +``` +born: d~ g > z{0} z{0} d~ legs: [(-1,[]), (21,[]), (23,[0]), (23,[0]), (-1,[])] +real: d~ g > z{0} z{0} d~ g legs: [(-1,[]), (21,[]), (23,[0]), (23,[0]), (-1,[]), (21,[])] +``` + +The front end is not the blocker. Everything downstream is. + +### The frame machinery is LO-only + +| piece | location | +|---|---| +| `boost_to_frame` (pure `boostx`, no rotation) | `Template/LO/SubProcesses/genps.f:1761` | +| `common /to_frame_me/frame_id` | `Template/LO/Source/run.inc:27` | +| call sites | `madgraph/iolibs/template_files/auto_dsig_v4.inc:183`, `:452`, `:488` | +| `mapid` (`ids(i)=btest(id,i)`) | `Template/LO/SubProcesses/cluster.f:128` | +| near-duplicate | `MadSpin/src/driver.f:1924` | + +`Template/NLO/` has **zero** occurrences of `frame_id`, `me_frame` or +`boost_to_frame`. + +### The six checks + +| | status | +|---|---| +| Born boost | **No.** `call sborn(p_born,…)` x20 in `fks_singular.f`, always raw partonic-CM momenta | +| Virtual boost | **No.** `Call BinothLHA(p_born,…)` — `fks_singular.f:7087` | +| Real boost | **No.** `call smatrix_real(pp,wgt)` — `fks_singular.f:4701` | +| Counterterm boost | **No**, for any of them — soft `sbornsoft`, collinear `sborncol_isr`/`sborncol_fsr`, soft-collinear, degenerate `sreal_deg`, `bornsoftvirtual`, `extra_cnt`. Each has its own reduced kinematics (`p1_cnt(0,1,0/1/2)`, `p_born`, `p_born_used`) and so needs its own frame | +| Boost before eikonal/splitting kernel | **No — and this is the hard part** (see M2 step 4) | +| run_card option | **No.** `me_frame`/`frame_id` are added in `RunCardLO` (`banner.py:4438-4440`), resolved at `:4869`, and the `frame` display block is appended only from the LO process-dependent setup (`:5191`). `RunCardNLO` (`:5758`) has none of it | + +### On rotations + +A rotation is **not** needed for the squared MEs. HELAS builds polarisation +vectors in each particle's own momentum-direction helicity basis, so a rotation +`R` gives `eps_lambda(Rp) = R eps_lambda(p) * exp(-i*lambda*dphi)` — a pure +per-particle phase that cancels in `|M|^2` when every external helicity is +fixed. Born, real and virtual are all rotation-invariant. + +Where orientation *does* matter is the collinear counterterm azimuthal phases +(`getaziangles`, and the ISR `cphi_mother=1` assumption). That is where the +"the two Z should be on the z axis" concern actually lands. + +## 2. Design decisions + +### D1 — Where to apply the boost + +Options: inside the generated ME wrappers (`SBORN`, `SMATRIX_REAL`, …), or +caller-side in `fks_singular.f`. + +**Recommend caller-side.** `SBORN` caches amplitudes against +`savemom(j,1)=p1(0,j), savemom(j,2)=p1(3,j)` +(`madgraph/iolibs/template_files/born_fks.inc:85-97`) and that cache is shared +with `sborn_sf`, `born_hel` and `extra_cnt` via `calculatedBorn`/`saveamp`. If +one entry point boosts and another does not, the cache silently returns +amplitudes from the wrong frame — a bug with no symptom. Boosting once per +kinematic configuration in the caller makes coherence structural, and gives the +azimuthal code access to the same boost vector. + +Concretely: one helper `boost_to_me_frame(p, npart, pboost_out, p_out)` in the +NLO template, plus a common block holding the boost vector for the current +configuration. + +### D2 — How `me_frame` indexes legs + +At LO, run_card positions == fortran positions, so `mapid` is enough. At NLO +this breaks twice: + +- `FKSLegList.sort()` (`fks_common.py:793`) reorders born legs relative to the + user's process, so user leg *n* != born position *n*. +- Born has `nexternal-1` legs, each real has `nexternal`, with the map given by + `set_pdg` (`Template/NLO/SubProcesses/chooser_functions.f:454-540`): real *k* + -> born *k* for `k < i_fks`, born *k-1* for `k > i_fks`, and `i_fks` has no + born counterpart. + +**Decided:** keep `me_frame` in **user process numbering**. + +Note the premise above was worse than stated: `sort_proc` (`fks_common.py:731`) +*renumbers* the legs (`leg['number'] = n + 1`) after permuting them, so the +user's numbering is not merely reordered, it is destroyed. Verified: the user +writes `p p > j z{0}` and the born comes out `g d > z{0} d` with the Z carrying +`number=3`, not 4. It therefore cannot be recovered downstream and has to be +captured before `sort_proc` runs. + +**Implemented in M0** (see the M0 section for status): + +- `fks_common.get_user_leg_order()` computes the FKS ordering on a copy and + returns the pre-sort numbers; `FKSProcess.__init__` calls it *before* + `sort_proc` and stores `self.user_leg_order`, which `FKSHelasProcess` carries + to the exporter. +- `frame_info.inc` per P dir holds only `frame_map_born(nexternal-1)`. +- **No `frame_map_real` table is needed.** The FKS convention already fixes the + real -> underlying-Born correspondence in terms of `i_fks` alone (`set_pdg`, + `chooser_functions.f:454-540`): real *r* -> Born *r* for `r < i_fks`, Born + *r-1* for `r > i_fks`, and `i_fks` itself has no Born counterpart. So + `get_frame_mask_real` derives the real mask at runtime from the Born one. + This drops the per-FKS-config table the plan originally called for. + +### D3 — Smoothness and equivariance + +The frame must be a continuous function of the momenta converging in the +singular limits, or the subtraction stops cancelling. "Rest frame of the sum of +the selected legs, recomputed from *each* configuration's own momenta" +satisfies this: as `xi -> 0` / `y -> 1` the Born's ZZ sum -> the real's ZZ sum. + +**Do not** compute the frame once from the real and reuse it for the Born — it +looks simpler and is wrong away from the limit. + +**The subtler requirement (equivariance).** A pure boost to the rest frame of a +system is *not* equivariant under longitudinal boosts: `B(L_z P) . L_z != B(P)` +whenever `P` has transverse momentum, the difference being a finite rotation. +So if the reduced Born lived in its own partonic CM rather than in the same +"tilde" frame as the real event, the polarisation axes of the two would differ +by an `O(xi)` rotation, the `1/(1-y)` pole would not cancel, and no azimuthal +fix could repair it. + +FKS bookkeeping gives us this for free in the ISR case, and it has been +checked: at `y=1` the ISR mapping applies to the final-state spectators only a +transverse boost whose rapidity vanishes identically +(`genps_fks.f:3055-3068`), and `xp` was initialised from the Born momenta +(`genps_fks.f:1206-1222`). In the tilde frame the real event's spectators are +therefore bit-identical to the Born's -> same sum of selected momenta -> same +`me_frame` boost for real and counter-event -> no residual Wigner rotation of +the Z polarisation axes. + +**Action:** assert this at runtime (or at minimum comment it) — check that the +boost 4-vector derived from `p_born_used` and from the real `p` agree at the +counter-event. It is cheap and it protects a property the whole subtraction +rests on. + +**Not yet confirmed for FSR:** `generate_momenta_massless_final` applies a +non-transverse `boostwdir2(chybst,…)` to spectators +(`genps_fks.f:2324-2339`). Confirm `shybst -> 0` in the collinear limit before +assuming FSR inherits the same property. + +## 3. Milestones + +### M0 — Plumbing, no physics — **DONE** + +Delivered: + +| what | where | +|---|---| +| `me_frame` + `frame_id` in the NLO run card | `banner.py` `RunCardNLO.default_setup`, `frame_id` resolved in `update_system_parameter_for_include` | +| `frame_block` exposed only for polarised massive legs | `RunCardNLO.blocks` + `create_default_for_process` | +| `$frame` placeholder | `Template/NLO/Cards/run_card.dat`, after the beam energies (mirrors the LO card) | +| `common/to_frame_me/frame_id` | `Template/NLO/Source/run.inc` | +| `mapid_frame`, `get_frame_mask_born`, `get_frame_mask_real`, `get_me_frame_boost`, `boost_to_me_frame` | `Template/NLO/SubProcesses/boost_to_frame.f` (new) | +| `boost_to_frame.o` in every executable | `makefile_fks_dir` `FILES` | +| `frame_info.inc` per P dir | `export_fks.write_frame_info_file` | +| user leg order captured pre-sort | `fks_common.get_user_leg_order`, `FKSProcess`, `FKSHelasProcess` | + +Gate met. The change is **runtime-inert by construction**, which is a stronger +statement than a bit-identical rerun: `frame_id` is assigned in `run_card.inc` +but read nowhere, and nothing calls the boost routines yet. Verified by grep +over a generated process dir. `p p > z j [QCD]` builds, links and runs clean. + +`get_me_frame_boost` returns `trivial=.true.` both when nothing is selected and +when the selected system already has exactly zero 3-momentum, so the default +`me_frame=[1,2]` costs nothing and cannot perturb existing results even once +call sites are added in M1. + +Unit-tested numerically (rest frame reached, invariant mass preserved, +4-momentum conservation preserved, in-place aliasing safe). + +Also fixed on the way: adding a block to `RunCardNLO` without a `$frame` +placeholder in the card template made `banner.py:3178` append the block at the +end of the file, which broke the `test_default` round-trip in +`tests/unit_tests/various/test_banner.py`. The placeholder is the fix. + +Original task list, for reference: + +1. `me_frame` + `frame_id` into `RunCardNLO` (`banner.py:5758`), mirroring + `banner.py:4438-4440`; replicate the `frame_id = sum(2**n)` resolution from + `:4869` and the `frame` display-block trigger from `:5191` into the NLO + process-dependent setup. +2. `integer frame_id / common/to_frame_me/frame_id` into + `Template/NLO/Source/run.inc`. +3. `boost_to_me_frame` + `mapid` into the NLO SubProcesses (port of + `genps.f:1761`, pure `boostx`, keep the trivial-boost short circuit). +4. Exporter writes `frame_info.inc` per P dir (D2). Add alongside the other + per-dir includes at `export_fks.py:502`. + +**Gate:** `frame_id=6` (default) reproduces current NLO results bit-for-bit on +an unpolarised `p p > z j [QCD]`. Nothing else changes yet. + +### M1 — `[LOonly=QCD]`: Born boost only — **DONE** + +Gate met. `p p > z{0} j`, matched PDF (nn23lo1), fixed scales (91.188), ptj 30, +etaj 4.0: + +| | LO madevent | NLO LOonly | | +|---|---|---|---| +| unpolarised (control) | 6966 +- 6 | 6964 +- 18 | 0.03% | +| `z{0}`, Z rest frame `me_frame=[3]` | 1370 +- 1.7 | 1366 +- 3.1 | 1.1 sigma | +| `z{0}`, partonic c.m. `me_frame=[1,2]` | 974.8 +- 1.0 | 977.3 +- 2.3 | 1.0 sigma | + +The unpolarised control matters as much as the polarised rows: it proves the +PDF/scale/cut/channel matching is right, so a polarised disagreement cannot be +blamed on setup. The frame dependence is large (40%) and both codes reproduce +it. + +Null test, `p p > z{0} z{0} [LOonly=QCD]`: `me_frame=[3,4]` gives +0.5392 +- 0.0016 pb and `me_frame=[1,2]` gives 0.5401 +- 0.0017 pb (0.4 sigma). +At Born level the ZZ system *is* the partonic c.m., so the boost must be the +identity and the answer must not depend on `me_frame` -- it does not. + +Incidentally this also confirms `nhel` is a variance knob, not a correctness +one, for polarised LO: `nhel=1` gives 1369 +- 1.7 against `nhel=0`'s +1370 +- 1.7. `help_polarization` currently implies `nhel=1` is required. + +Converted call sites (all three reachable with `abrv='born'`): +`compute_born` (`fks_singular.f:37`), `include_multichannel_enhance` +(`:1505` and the cache-isolated `:1550`), and `bornsoftvirtual` (`:6896`, +not reachable in LOonly but converted for M3). + +Guard relaxed for `LOonly` only (`madgraph_interface.py:1240`), including the +massive-particle rejection; `[QCD]` and the other NLO modes are still refused. + +Two things learned the hard way; record them before touching M2. + +**The fixed-order Born does not come from `bornsoftvirtual`.** It comes from +`compute_born` (`fks_singular.f:1`, the first routine in the file), called from +`driver_mintFO.f:445`, which does its own `call sborn(p_born,wgt_c)` and then +`add_wgt(2,...)`. `bornsoftvirtual` (`:494`) jumps straight to label 549 when +`abrv='born'` and leaves `amp_split_wgtnstmp` at zero, so its `add_wgt(3,...)` +contributes nothing. Converting only `bornsoftvirtual` therefore changed the +answer without ever boosting the Born that is actually integrated. Both sites +are now converted. + +**B5 is not hypothetical.** With only `bornsoftvirtual` converted, the run did +not simply stay unboosted: it moved by ~3% (974.8 -> 946.7). Two callers of +`SBORN` disagreed about the frame within one event, and they share both the +`amp_split` common and the `calculatedBorn`/`savemom` cache, so the second +caller either reused amplitudes from the other frame or thrashed the cache. +This is exactly the failure mode D1 predicts, and it is silent -- no crash, no +warning, just a wrong number. Every `SBORN` caller reachable in a given mode +must be converted together, or none. + +Original task list: + +`[LOonly=QCD]` generates a fake FKS config with no reals and no virtuals +(`export_fks.py:3676`), so this milestone touches exactly one ME call and zero +counterterms. Cheapest possible validation of D1 + D2 + the run_card. + +- Boost `p_born` before `call sborn` at the LOonly-reachable sites. +- Relax the guard at `madgraph_interface.py:1240` for `LOonly` only, and drop + the massive-particle rejection at `:1257` on that branch (today it blocks + `z{0}` even in the allowed `noborn` mode). + +**Gate — acceptance test 1 (`p p > z{0} j`):** `[LOonly=QCD]` +`calculate_xsect LO` vs. LO madevent `generate p p > z{0} j`, same `me_frame`, +same cuts/PDF/scales, `group_subprocesses False`. Must agree inside MC error. +This validates the NLO Born boost against an independently-validated +implementation. Pin `nhel` identically on both sides. + +`p p > z{0} z{0}` is a deliberate **null** test here: at Born level the ZZ +system *is* the partonic CM, so the boost must be a no-op and the polarised +LOonly result must be `me_frame`-independent. + +### M2 — `[real=QCD]`: reals + counterterms (the hard milestone) + +**Step 0 (ISR half) DONE.** The azimuthal factor can now be rebuilt covariantly. + +`genps_fks.f` stores `xij_kperp` (common `/cxij_kperp/`) next to `xij_aor` at +both ISR generation sites (`generate_momenta_initial`, +`generate_momenta_initial_noevpr`), zeroed wherever `xij_aor` is zeroed. +`azifact_from_kperp` in `boost_to_frame.f` rebuilds `-exp(2 i psi)` from the +mother direction and that vector. + +The reconstruction turned out cleaner than the plan assumed. Using the standard +helicity basis of the mother direction, + + e1 = ( cos(th)cos(ph), cos(th)sin(ph), -sin(th) ) + e2 = ( -sin(ph), cos(ph), 0 ) + +reproduces **both** incoming legs with no `idir` bookkeeping at all: + + j_fks=1, n=+z : psi = phi_i -> -exp( 2 i phi_i) + j_fks=2, n=-z : psi = pi - phi_i -> -exp(-2 i phi_i) + +because the doubled angle turns the basis flip into a harmless 2 pi. This is +*why* B3 is right that the `R_y(pi)` flip must be deleted rather than boosted: +it was only ever compensating for a basis convention that the psi form handles +by construction. + +Validated at Lambda=1 against `-exp(2 idir i phi_i)` over 37 azimuths on both +beams: worst deviation 2.4e-16, i.e. one ulp. Exact bit-identity is not +reachable through this route (`(a+ib)^2` vs `exp(2 i phi)` differ in the last +ulp), so the shipped code keeps using `xij_aor` whenever the boost is trivial +and only takes the covariant route when a frame is actually requested -- +unpolarised runs stay bit-identical by construction rather than by measurement. + +A longitudinal boost leaves the reconstruction *exactly* unchanged (0.0, not +just small), confirming the plan's prediction that a longitudinal-only +`me_frame` is a genuine null and therefore a useful intermediate test. + +**Step 0 (FSR half) DONE, and B4 is resolved -- both parts favourably.** + +*Which branch runs.* The same one as ISR. `sreal` dispatches to +`sborncol_fsr`/`sborncol_isr` on `pmass(j_fks)` and `j_fks<=nincoming`, but the +collinear and soft-collinear counter-events pass the literal `one` for +`y_ij_fks` (`fks_singular.f:793`, `:904`) in both cases. So FSR also always +takes `azifact = xij_aor`, and its `IXXXXX`/`OXXXXX` branch is live only in the +sliver `vtiny=1e-8 <= 1-y < tiny=1e-6` -- never during `colltest`, where +`tiny=1d-12`. + +*Equivariance.* Holds for FSR too. The spectator boost has +`shybst = -(shat-sumrec^2)/(2 sumrec sqrt(shat))` with +`recoil = p_total - p_mother`, which in the partonic c.m. vanishes exactly when +`E_m = |p_m|`, i.e. for a massless mother. Since +`m_mother^2 = 2 E_i E_j (1-y) -> 0`, `shybst = O(1-y) -> 0` in the collinear +limit. So the reduced Born and the real event share a frame in the singular +region and there is no residual rotation of the polarisation axes. + +*The FSR reconstruction needs no new machinery.* The emission is generated +about a `+z` mother as `(cos(phi_i),sin(phi_i),0)` and then passed through +`rotate_invar`, which is `R_z(phi) R_y(theta)` and therefore maps `x,y` onto +exactly the `e1,e2` of the helicity basis above. So the rotated vector has +azimuth `phi_i` again, and the *same* `azifact_from_kperp` serves both cases. +The only difference is a conjugation, which is just spacelike vs timelike +splitting: + + ISR net factor = -dconjg( azifact_from_kperp(mother, kperp) ) + FSR net factor = - azifact_from_kperp(mother, kperp) + +Checked at Lambda=1 against the shipped `-(cphi_m - i sphi_m)^2 * xij_aor` over +672 configurations of `(theta_m, phi_m, phi_i)`: worst deviation 1.4e-15. + +This also retires the worry that FSR would need its own analysis: the two +cases differ by a conjugation, not by structure. + +`[real=QCD]` gives the full FKS subtraction without virtuals — everything hard, +nothing MadLoop. + +1. Boost the real momenta before `call smatrix_real(pp,wgt)` + (`fks_singular.f:4701`). +2. Boost each reduced configuration with its *own* frame: `p1_cnt(0,1,0/1/2)`, + `p_born`, `p_born_used`, and the `extra_cnt` argument. Sites: + `fks_singular.f:494` (`bornsoftvirtual`), `:705`/`:793`/`:904` (the + `sreal`/`sreal_deg` counterterms), plus the ~20 `call sborn` sites. +3. Soft is easy: the eikonal `p_m.p_n/(p_m.k)(p_n.k)` is invariant and + `sborn_sf` is colour-linked only (no spin correlation), so `sbornsoft` needs + nothing beyond a boosted `p_born`. +4. **The azimuthal phases — the real work.** `SBORN` returns + `ANS(2) = BORNTILDE`, the +/- gluon-helicity interference built from + `JAMPH(1,.)`/`JAMPH(2,.)` (`born_fks.inc:277`). The correction multiplying + it collapses to `-exp(-2*i*psi)`, with `psi` the emission azimuth about the + mother in the HELAS `(e1,e2)` basis. + + **No extra Wigner phase is needed on `BORNTILDE`.** The mother's + little-group phase `exp(-+2*i*theta_W)` cancels against the `exp(+2*i*theta_W)` + picked up by a *boosted* `/[ij]` — same `theta_W`, since i, j and the + mother are collinear. The controlling identity, exact in any single frame + and for any mother direction, is `arg(azifact)/2 = phi_m + psi`. There is no + `theta_m` dependence (`R_y(theta)` has a real SU(2) matrix), which is + exactly why `getaziangles` returning only `phi` suffices. The generator's own + analytic limits confirm it independently: `genps_fks.f:2344` gives + `-exp(2i(phi_mother+phi_i))` (FSR), `:3099` gives `-exp(2*idir*i*phi_i)` (ISR). + + The massive-leg (Z) Wigner transformation is *not* a phase but a genuine + helicity mixing — that is the polarised observable itself, not an error to + correct. It is why `SBORN` must be handed boosted momenta. + + **ISR — BLOCKER B1, and the branch that actually runs.** + `fks_singular.f:5034-5035` hardcodes `cphi_mother=1, sphi_mother=0`; that is + numerically identical to `getaziangles` on a beam-axis momentum + (its `sth.ne.0` guard at `:4111` returns exactly `(1,0)` for `cth=+-1`) and + it dies under a boost. But repairing it is not the main problem: + + The ISR counter-events are called with `y_ij_fks = one` **exactly** + (`fks_singular.f:793`, `:904`; `y_ij_fks_matrix(1)=y_ij_fks_matrix(2)=1.d0` + at `genps_fks.f:2742`), so `1d0-y_ij_fks .lt. vtiny` always fires and + `azifact = xij_aor` is taken every time (`:5006-5007`). The + `IXXXXX`/`OXXXXX` spinor branch below it is live only in the sliver + `vtiny=1e-8 <= 1-y < tiny`, where `tiny=1d-6` in production and `1d-12` + under `colltest` (`:4658`) — i.e. it is **never exercised by + `test_soft_col_limits`**. Both branches still have to be made consistent or + that sliver becomes a discontinuity nothing catches. + + `xij_aor = -exp(2*idir*i*phi_i_fks)` is a precomputed partonic-CM constant + (`genps_fks.f:3096-3101`) and a `0/0` limit of `/[ij]`. A boost maps + exactly-parallel null vectors onto exactly-parallel null vectors, so **the + degeneracy survives in every frame — it can never be recomputed after + boosting.** The azimuthal direction must be carried through as a stored + 4-vector. + + Proposed fix: store `k_perp_dir = (0, cos(phi_i_fks), sin(phi_i_fks), 0)` + alongside `xij_aor`. One vector covers both beams (for `j_fks=2` the basis + is `(x, -y)`, so it yields `psi = -phi_i` — precisely the `idir` sign + already in `xij_aor`, and the `idir` special case disappears). Then boost + it, project out the mother component, take + `psi' = atan2(k'_perp . e2, k'_perp . e1)`, apply `-exp(-2*i*psi')`. At + `Lambda=1` this reproduces current numbers bit-for-bit, so it lands as a + safe refactor **before** any boost is wired in. + + `azifact` is invariant under independent positive rescalings, so + `p_i_fks_ev` being the rescaled null vector `p_i/xi` + (`genps_fks.f:3085-3094`) is harmless: `Lambda(lambda*p) = lambda*(Lambda p)`. + + **ISR is a package of three changes** — do any subset and you apply part of + the correction, which is worse than the current code: + (a) carry `k_perp_dir` through the boost as above; + (b) `getaziangles` on the boosted mother instead of the hardcode; + (c) **delete** the `j_fks.eq.2` `R_y(pi)` flip (`:5013-5020`) — after a boost + it no longer maps the mother onto +z. Note `montecarlocounter.f:3000-3006` + double-encodes the same thing (both the rotation and `cphi_mother=-1.d0`); + numerically identical today since only the square enters, but do not carry + both forward. + + **FSR:** `getaziangles(p_born(0,imother_fks), …)` (`:4844`) and `azifact` + built by `IXXXXX/OXXXXX` on `p_i_fks_ev`, `p(.,j_fks)` (`:4830`) must both be + evaluated on momenta carrying the *same* boost as the Born. The FSR + equivalent of B1 has not been checked — establish whether the FSR + counter-events also collapse onto the `vtiny` branch. + + **Scope narrowing:** `Qterms_reduced_spacelike` returns a real, z-only + number (`:5538-5584`), non-zero only for a **vector mother** (`col1=8`); the + `abs(m_type).eq.3 .or. ch_m.ne.0d0` branch that zeroes `Q` and `wgt1(2)` is + the same statement one level up. So **only gluon-mother ISR configurations + are affected at all** — a process whose ISR mothers are all quarks passes + the collinear test even with a completely wrong azimuthal phase. Choose the + validation channel deliberately. + + This is the frame/orientation consistency issue: it lands on the + counterterms, not on `|M|^2`. + +5. **Same construct elsewhere**, needing the same treatment if polarised + aMC@NLO (not just fixed-order) is in scope: `montecarlocounter.f:2969-3012` + and `montecarlocounter_alt.f:1445-1533`. `sreal_deg` (`fks_singular.f:5844`) + has no azimuthal term — it needs the boost for the ordinary Born only. + +**Step 1 (the B5 sweep) DONE.** Every ME entry point now goes through a frame +wrapper, so the boost is a pure function of the momenta passed in. + +That is what makes B5 go away *structurally* rather than by inspection: two +callers that pass the same momenta necessarily get the same boost, and two +callers that pass different momenta were already required to reset the cache +between them. No reachability argument is needed, which matters because the +reachability argument is exactly what I got wrong twice in M1. + +Wrappers in `boost_to_frame.f`: `sborn_frame`, `sborn_sf_frame`, +`extra_cnt_frame`, `smatrix_real_frame` (the last uses the *real* mask, derived +from the Born one and `i_fks` of the current FKS configuration). + +43 call sites converted: + +| file | sites | +|---|---| +| `fks_singular.f` | 30 | +| `montecarlocounter.f` | 5 | +| `montecarlocounter_alt.f` | 5 | +| `add_write_info.f` | 1 (this retires **B7**) | +| `check_poles.f` | 1 | +| `test_soft_col_limits.f` | 1 | + +Deliberately **not** converted, and why: the EW Sudakov paths +(`check_sudakov*.f`, `ewsudakov_functions_dummy.f`, +`sa_ewsudakov_dummyfcts.f`) are a separate feature, and `symmetry_fks_v3.f` is +`gensym`, a separate executable whose Born calls only seed the integration +grid. Neither shares a process with the integration, so neither can trip the +cache. + +**The sweep also exposed a wrong default, worth recording.** With all 43 sites +converted, unpolarised `p p > z j [QCD]` moved from 2.854e+04 to 2.834e+04 ++- 1.9e+02 -- only 0.74 sigma, but the two previous runs had been *bit* +identical, so the pipeline is deterministic and a change meant a real change. + +Cause: MadFKS does not hand the matrix elements momenta in their own partonic +c.m. `shy_lbst = -xi_i_fks*yijdir/bstfact` is non-zero for any real emission, +so `xp(0,1) != xp(0,2)` (`genps_fks.f:3074-3082`) and the real event lives in a +frame boosted along z. Honouring `me_frame=[1,2]` literally therefore applies a +genuine longitudinal boost to every configuration -- an identity for `|M|^2`, +but not bit for bit, and enough to send the VEGAS grids down a different path. + +Confirmed by control: the same build with `me_frame=[0]` (`frame_id=1`, every +mask empty, boost skipped) reproduces 2.854e+04 +- 1.9e+02 exactly. So the +sweep itself is sound and the whole shift came from the boost. + +The semantics were right -- boosting to the initial-parton rest frame is what +`[1,2]` means at LO, and the M1 cross-check agreed with LO at 1.0 sigma -- but +the *default* was wrong: an unpolarised run must not pay for machinery it never +asked for. `frame_id` is now 0, meaning "skip", unless `me_frame` actually +appears in the run_card: + +| case | frame_id | | +|---|---|---| +| unpolarised (no `me_frame` in the card) | 0 | boost skipped, bit-identical | +| polarised, default | 6 | partonic c.m., as at LO | +| polarised, `me_frame=[3]` | 8 | Z rest frame | + +The LO path is untouched: there the momenta really do arrive in the lab frame, +so `[1,2]` has always been a meaningful boost and still is. + +**Skip the identity boost structurally, as LO does.** Applying a boost that is +mathematically the identity is not free -- it goes through `boostx` and +perturbs the momenta in the last bits, which is enough to move an adaptive +integration. LO avoids this in *two* places, and the NLO port now folds in +both: + +- the call site skips `frame_id=6` outright (`auto_dsig_v4.inc:183`); +- `boost_to_frame()` skips the all-final-state selection (`genps.f:1782`), + whose comment says of the initial-state case: *"1 1 0 0 0 .... should not go + within this function"*. + +`get_me_frame_boost` now returns `trivial` for four cases: nothing selected, +exactly the initial state, exactly the whole final state, and an already-at- +rest system. Relying on the last one alone (as the first version did) is not +enough at NLO, because the tilde frame means it never fires. + +**`me_frame` must not be built from the initial state at NLO -- it is not +infrared safe.** This is a stronger statement than the numerical one above and +it constrains what the feature may *accept*, not just what it skips. The real +emission and the reduced Born carry momentum fractions that differ by a finite +amount even in the singular limit, so a frame defined from the initial state is +discontinuous across that limit and the subtraction stops cancelling. No +azimuthal fix can repair it -- this is the same failure mode as the +equivariance discussion in D3, but caused by the frame *definition* rather than +by the mapping. + +Consequently `RunCardNLO.check_validity` rejects an `me_frame` that mixes +initial-state legs with final-state ones. `[1,2]` and the full final state stay +accepted, since they merely name the partonic c.m. and are skipped. Frames +should be defined from final-state particles only, which is what the physics +wants anyway: the polarised system (the Z, the ZZ pair) is a set of final-state +spectators, present unchanged in both the real and the Born. + +Re-validated after this change -- unpolarised 2.854e+04 +- 1.9e+02 (baseline), +`z{0} j` Z-frame 1366 +- 3.1 and partonic c.m. 977.3 +- 2.3, both identical to +the M1 numbers. The skip reproduces the boost exactly for `[LOonly=QCD]`, as it +must: with no emission `shy_lbst=0` and the Born really is in the partonic c.m. + +**Step 2 (azimuthal wiring) DONE. The M2 gate PASSES and `[real=QCD]` is +enabled.** + +`azifact_me_frame` rebuilds the frame from the same Born momenta the matching +`sborn_frame` call gets, boosts the mother and `xij_kperp` together, and +returns `-exp(2 i psi)`. Both `sborncol_isr` and `sborncol_fsr` use it in their +boosted branch and, there, apply neither the `cphi_mother=1` shortcut, the +`R_y(pi)` flip, nor `getaziangles` -- the psi form carries all of them (**B3**). +The legacy branch is untouched, so unpolarised runs stay bit-identical. + +It first failed the gate: `test_ME` on `p p > z{0} j [real=QCD]` with +`me_frame=[3]` gave soft 0.48, collinear 0.37. The section below is the +diagnosis. It was *not* a frame-consistency bug and the azimuthal work was +never implicated. + +### The failure: the quantisation axis was decided by rounding + +**Root cause.** For a one-leg `me_frame` selection -- `me_frame=[3]`, "the Z +rest frame", which is the standard way to ask for Z polarisation -- the boost +puts that leg at rest, and *there the boost has to be exactly right rather than +right to rounding*. + +HELAS changes convention at exactly zero. `vxxxxx` (`aloha_functions.f`) builds +the polarisation vectors of a massive vector from + + pp = min(|p(0)|, sqrt(p(1)^2+p(2)^2+p(3)^2)) + +and branches on `pp.eq.rZero`: at exactly zero it uses the **z axis** of the +current frame as the quantisation axis, and otherwise it uses the **momentum +direction**, giving `eps_L = (pp/m, E*p/(m*pp))`. + +`boostx` does not reach `p=0` exactly. With `q=(E,-P)` and `p=P` it forms +`p(i)+q(i)*lf` where `lf` is 1 only up to the rounding of `(q(0)-m)+p(0)`, so +the residual is a few `1d-14` pointing in a direction that is pure noise. +Whether it rounds to zero is a coin flip, decided independently for the real +and for the reduced Born. Every configuration that misses zero gets `eps_L` +pointed along rounding noise instead of along z, and the subtraction is then +comparing two different observables. + +**Evidence.** Instrumented `boost_to_me_frame` to print `|p|` of the selected +leg after the boost, alongside the soft-limit ratio table, on +`p p > z{0} z{0} [real=QCD]` with `me_frame=[3]`. One-to-one, no intermediate +values: + +| residual `|p_Z|` after the boost | soft ratio | +|---|---| +| 0 (exact) | 0.9033, 0.9901, 0.9990, 0.9999, 1.0000, 1.0000, 1.0000 | +| 2.95e-14 | 0.0021716, 0.0021716, 0.0021716 | + +The Born landed on exact zero for every point sampled; the real missed it 30% +of the time. A factor 460 in `|M|^2`, from a 3e-14 momentum. + +**Fix** (`boost_to_me_frame`, `Template/NLO/SubProcesses/boost_to_frame.f`): +when the selection is a single leg, zero that leg's 3-momentum exactly after +boosting. This imposes the defining property of the frame instead of hoping +for it, and removes the choice. Only `nsel=1` needs it -- with two or more +selected legs it is their *sum* that is at rest, no individual leg sits on the +branch point, and the residual of the sum never reaches HELAS. + +### `p p > z z` was the right thing to run first + +It was decisive, but not in either of the two ways this document predicted. +`me_frame=[3]` on `p p > z{0} z{0} [real=QCD]` reproduces the `z{0} j` symptom +exactly -- soft 0.55/0.68, collinear 0.54/0.57, uniform -- and there every ISR +mother is a quark, so `Q` vanishes identically and the azimuthal phase is +provably not involved. That reduced the search to the boost itself on a process +that builds in a quarter of the time. + +The `me_frame=[3,4]` spelling suggested in the previous revision would have been +a much weaker probe: it fails only a split-order soft test (0.43/0.30) while the +summed test passes, and for a *different*, purely numerical reason -- see below. + +**What the D3 instrumentation actually showed.** The boost 4-vectors were +printed for the real and for the reduced Born at the same counter-event, as D3 +asks. They converge cleanly: on `me_frame=[3,4]` the real's boost runs +`(553.24, -13.00, 25.52, ...)` down to `(552.50, -1.2e-8, 2.4e-8, ...)` against +the Born's constant `(552.50, 0, 0, ...)`, i.e. the transverse part vanishes +like xi exactly as D3 requires. **The working hypothesis of a finite frame +mismatch was wrong**, and so was the reasoning that led to it. There is no need +for a permanent D3 assertion: the property holds, and the thing that broke was +downstream of it. + +### A second, independent finding: the same bug is live at LO + +`boost_to_frame` (`Template/LO/SubProcesses/genps.f:1759`) has the identical +defect, and `me_frame` is a *shipped* LO feature. Measured on LO madevent, +`p p > z{0} j`, `me_frame=[3]`, nn23lo1, fixed scales 91.188, ptj 20, etaj 5: + +| | cross-section | +|---|---| +| before the fix | 1399 +- 3.4 pb | +| after the fix | 1424 +- 3.8 pb | + +1.8%, 5 sigma. So a fraction of LO events had been getting a rounding-noise +polarisation axis all along. The same fix is applied there. + +This also matters for M1: that milestone's gate is NLO `[LOonly=QCD]` against +LO madevent, so fixing only the NLO side would have made the two disagree. +Re-validated with both fixed, same cuts and PDF: LO madevent 1424 +- 3.8 vs +NLO `[LOonly=QCD]` 1427 +- 3.0, **0.6 sigma**. + +### A precision floor, and the betamin guard + +Separately, `me_frame=[3,4]` on `p p > z{0} z{0}` failed the *split-order* soft +test at 0.43/0.30 while the summed test passed. That one is not physics. + +That selection is the whole Born final state, so its boost degenerates to the +identity as xi -> 0. Applying a near-identity boost is actively harmful: it +tilts the incoming partons off the beam axis by O(beta), and the HELAS massless +spinors are built from `sqrt(p(0)+p(3))`. On the axis that sum is exactly zero +and the exact branch is taken; tilted by beta it is `E*beta^2/2`, computed as a +difference of two numbers of size E, hence carrying absolute error `ulp(E)` and +relative error `2*eps/beta^2`. Measured down the soft scan, `p(0)+p(3)` of the +boosted incoming leg and the deviation of the ratio from 1: + +| `p(0)+p(3)` | 3.3e-9 | 3.3e-11 | 3.4e-13 | 0 | 0 | 0 | +|---|---|---|---|---|---|---| +| ratio - 1 | 9e-6 | 6.5e-4 | 2.1e-2 | 2e-6 | 4e-7 | 2e-6 | + +Exactly `1/beta^2` while `p(0)+p(3)` is a few ulp, clean again once it +underflows to exactly zero. + +Balancing what the boost delivers (O(beta)) against what it costs +(`30*2*eps/beta^2`, the 30 measured above) gives `beta* = 2.4e-5`. +`get_me_frame_boost` now skips the boost below `betamin=1d-4`, one safety +decade above that. Scanned: `1d-5` is not enough (0.16 still failing), `1d-4` +and `1d-3` both give 0.01. This costs no physics -- it fires only where the +selected system is at rest to one part in 10^4, which never happens in an +integration, only in the artificial deep-soft scan of `test_soft_col_limits`. +Confirmed neutral: the `[real=QCD]` cross-section is 2.604e+03 +- 2.8e+01 pb +without the guard and 2.600e+03 +- 2.2e+01 pb with it, 0.1 sigma. + +### Gate results + +`p p > z{0} j [real=QCD]`, `me_frame=[3]`, nn23lo1, fixed scales -- the +process this document names as the **discriminating** one, and the one whose +gluon-initiated channels exercise the azimuthal `Q` term: + +| | test_ME | +|---|---| +| before | 268 FAILED of 512 test lines, all 12 P dirs | +| after | **0 FAILED of 512**, worst failure fraction 0.08 against a 0.30 threshold | + +`calculate_xsect NLO` then runs end to end: **2.600e+03 +- 2.2e+01 pb**. + +`p p > z{0} z{0} [real=QCD]`, all frames, `P0_ddx_z0z0`: + +| `me_frame` | frame_id | before | after | +|---|---|---|---| +| none | 0 | passes (0.01) | passes (0.01), unchanged | +| `[3]` | 8 | FAILS 0.55/0.68 | passes, worst 0.08 | +| `[3,4]` | 24 | FAILS 0.43/0.30 | passes, worst 0.01 | + +Unpolarised runs are untouched **by construction**, not by measurement: +`frame_id=0` makes every mask empty, `get_me_frame_boost` returns `trivial` on +`nsel=0` before reaching either new branch, and `boost_to_me_frame` returns +before the `nsel=1` block. + +The guard at `madgraph_interface.py` now accepts `real` alongside `loonly`. + +### One thing that was checked and left alone + +The "exactly the whole final state" skip in `get_me_frame_boost` is +multiplicity-dependent, so it does not mean the same thing for the Born and for +the real: `me_frame=[3,4]` on `p p > z z` is the whole Born final state but only +two of three real legs. That asymmetry looked like the bug and is not: removing +the skip leaves the failure fractions **bit-identical**, because `p_born` is +handed to us in the Born partonic c.m., so the zero-3-momentum test catches the +same case anyway. Reverted, and the reason recorded in the source. + +**Recommended order within M2:** + +0. Refactor the ISR azimuthal factor to the `psi`-form at `Lambda=1` (the + `k_perp_dir` construction above, boost = identity). Demand **bit-identical** + `test_ME`. This de-risks the hardest change before any frame logic exists. +1. Wire in a **longitudinal-only** `me_frame` boost. Must also be unchanged — + the mother stays on the beam axis and the z-rotation Wigner phase of a + z-directed massless particle is trivial, so the existing hardcode remains + exactly correct. A useful null. +2. Then the general boost. + +**Gate — `test_soft_col_limits`.** Already runs automatically as `test_ME` on +every `calculate_xsect` (`amcatnlo_run_interface.py:5494`, +`Template/NLO/SubProcesses/test_soft_col_limits.f`; the collinear scan at +`:617-663` drives `1-y = 1e-2 … 1e-10`). The `Q(z)` azimuthal term is the +*same order* in the collinear limit as the `AP(z)` term, so a wrong phase does +not cancel — the ratio plateaus off 1. Run it before any cross-section +comparison. + +**Gate — acceptance test 2 (`p p > z{0} z{0}`):** the Born boost is trivial but +the real has `z z j`, so in the *finite* region the ZZ frame != partonic CM. +This exercises the real boost with the Born held fixed. + +**But it does not discriminate the azimuthal phase.** In the singular region — +which is precisely where `test_soft_col_limits` probes — the reduced Born is +`p p > z z` with zero pT(ZZ), so the `me_frame` boost degenerates to +longitudinal and the phase issue does not bite. `p p > z{0} z{0}` will pass +either way. The discriminating processes are those whose `me_frame` system +recoils already at Born level: **`p p > z{0} j`** and **`p p > z{0} z{0} j`**, +restricted to gluon-initiated channels (see the scope narrowing above). + +### M3 — `[QCD]`: virtual — **DONE, GATE PASSES** + +`binothlha_frame` boosts the momenta handed to `BinothLHA`, which passes them +straight to `sloopmatrix_thres` and takes `born_wgt` from the caller rather +than recomputing it, so one call site covers the virtual. + +**Resolved.** The integration was already right; the *gate* was the thing +comparing two frames. `check_poles.f` had been converted by the M2 sweep on +its Born line only — `call sborn_frame(p_born, born)` — while the line +immediately below it still called `BinothLHA(p_born, …)` directly, on the +unboosted `p_born`. So the driver handed MadFKS a Born in the Z rest frame +and MadLoop a phase-space point in the partonic c.m., and then compared the +poles the two produced. Routing it through `binothlha_frame`, exactly as +`bornsoftvirtual` already does at `fks_singular.f:7123`, makes the check +pass. The rest of this section is the original failure report followed by the +diagnosis, the evidence and the gate results. + +#### What the failure was + +`check_poles` on `p p > z{0} j [QCD]` with `me_frame=[3]`: **poles do not +cancel**, 20 miscancellations. `test_ME` was clean on the same run (0 FAILED, +40 PASSED), which correctly said M2 was unaffected and the problem was on the +virtual side. + +The pole coefficients were read as saying it is not a normalisation: + + COEFFICIENT DOUBLE POLE: MadFKS -3.23e-3 OLP -4.88e-5 ratio 66 + COEFFICIENT DOUBLE POLE: MadFKS -5.13e-3 OLP -9.43e-5 ratio 54 + COEFFICIENT SINGLE POLE: MadFKS -4.06e-3 OLP -7.56e-5 ratio 54 + +"the ratio is O(50-60) and *not* constant, so MadLoop is evaluating a +different polarisation state rather than the same one with a wrong factor." + +**That reading was wrong, and it is what sent the search into MadLoop.** The +three lines above come from three *different* phase-space points (points 4, 5 +and 6 of the scan), one line each — so the varying ratio is a comparison +across points, not within one. Grouped by point, the ratio is *exactly* +constant: point 1 gives 58.03 for the double pole and 58.03 for the single, +point 2 gives 81.01 and 81.01, and only the per-point value moves. Both poles +are proportional to the Born, so a per-point constant ratio says the two +sides disagree about the Born by an overall factor, i.e. about the frame. +Group diagnostics by phase-space point before drawing conclusions from them. + +Two MadLoop settings were suspected and both were correctly **excluded by +test**: + +- `NRotations_DP`/`NRotations_QP` re-evaluate the loop at a rotated PS point and + expect `|M|^2` unchanged -- which fails for a particle at rest, whose axis is + the frame z axis and does not rotate with the momenta. But both **default to + 0**, so the rotation test never runs. Latent hazard, not this bug. Keep them + at 0 for polarised runs. +- `ImprovePSPoint` (default 2) deforms the PS point to restore exact + onshellness and could move the deliberately-zeroed leg off zero. Setting it + to `-1` (no deformation in double precision) **does not fix it**: still 20 + miscancellations, ratio still ~58. + +#### The diagnosis + +**MadLoop's polarisation restriction was never the problem, and that was +checked before anything else was touched.** MadLoop reads its helicity +configurations from `MadLoop5_resources/HelConfigs.dat`, and for +`P0_gd_z0d` that file holds exactly the eight rows of `born.f`'s `NHEL` +table, with the Z entry pinned to `0`: + + -1 1 0 -1 -1 -1 0 -1 1 1 0 -1 1 -1 0 -1 + -1 1 0 1 -1 -1 0 1 1 1 0 1 1 -1 0 1 + +`NCOMB=8` in `loop_matrix.f` against `MAX_BHEL=8` in `born_nhel.inc`. The +Born and the virtual sum over the same helicity set. `hel_avg_factor` is +restored as `loop_exporters.py:1552` intends. Nothing to fix there. + +**Instrumentation closed it.** A `diag_frame` routine in `check_poles.f` +prints, for each point, the selected leg's `|p|` before and after the boost +and the Born evaluated both ways (`sborn` on the raw `p_born`, and +`sborn_frame`), with `calculatedborn` reset between the two so the shared +amplitude cache does not fight. The Z goes from `|p| = 318.6` to exactly `0`, +as `boost_to_me_frame` guarantees for `nsel=1`, and: + +| point | Born(me_frame)/Born(partonic c.m.) | MadFKS/OLP pole ratio | rel. diff | +|---|---|---|---| +| 1 | 58.0338377763 | 58.0338377763 | 1.7e-14 | +| 2 | 81.0065550136 | 81.0065550136 | 6.9e-14 | +| 3 | 64.0169796788 | 64.0169796788 | 4.4e-15 | +| … | | | | +| 20 | 57.8773835460 | 57.8773835460 | 1.3e-14 | + +All 20 points agree to 1e-13. The pole mismatch *is* the ratio of the two +Borns, with nothing else in it. + +The chain that produces it: `getpoles` (`fks_singular.f:7574`) ignores its +`p` argument for the Born and reads `common/pborn/p_born`, then calls +`sborn_frame` on it — so the MadFKS poles are always in the me_frame. +MadLoop's poles come from whatever momenta `BinothLHA` was handed. In the +integration that is `binothlha_frame`'s boosted copy and the two agree; in +`check_poles.f` it was the raw `p_born` and they did not. + +#### The fix + +One line, `Template/NLO/SubProcesses/check_poles.f:228`: + +```fortran +- call BinothLHA(p_born, born, virt_wgt) ++ call binothlha_frame(p_born, born, virt_wgt) +``` + +`check_poles.f` is compiled only into the `check_poles` executable (`POLES` +in the P-dir makefile; the integration uses `driver_mintFO.o`), so no +cross-section can move as a result of this change — only the diagnostic that +was misreporting. + +Unpolarised runs are unaffected, by measurement as well as by construction: +with `frame_id=0` every mask is empty, `get_me_frame_boost` returns `trivial` +and `boost_to_me_frame` copies the momenta unchanged. Rebuilt +`P0_gd_zd/check_poles` of unpolarised `p p > z j [QCD]` both ways — the two +logs are **byte-identical**. + +The two settings excluded by the previous session stay excluded, and the +passing run now confirms it positively rather than by elimination: it ran +with the shipped defaults `ImprovePSPoint = 2`, `NRotations_DP = 0`, +`NRotations_QP = 0`. + +#### Gate results + +`p p > z{0} j [QCD]`, `me_frame=[3]`, nn23lo1, fixed scales 91.188, ptj 30, +etaj 4.0, `req_acc_fo 0.05`: + +| | before | after | +|---|---|---| +| `check_poles`, `P0_gd_z0d` | 0 of 20 pass | **20 of 20** | +| `check_poles`, all 12 P dirs | — | **20 of 20 in every one** | +| `test_ME`, all 12 P dirs | passed | passed, unchanged | +| `calculate_xsect NLO` | — | **2.193e+03 +- 1.2e+01 pb** | + +Agreement is now at the level of the loop-library accuracy, not the +tolerance: e.g. MadFKS `-5.1257915456849783e-03` against OLP +`-5.1257915456856262e-03`. + +Unpolarised control, `p p > z j [QCD]`, same cuts and scales: `check_poles` +20 of 20 in all 12 P dirs, `calculate_xsect NLO` 1.029e+04 +- 5.5e+01 pb. + +That control also gives a physics sanity check the pole cancellation cannot, +since the M1 table has the Born of both at the same cuts and scales: + +| | Born (`[LOonly=QCD]`) | NLO (`[QCD]`) | K | +|---|---|---|---| +| unpolarised | 6964 +- 18 | 10290 +- 55 | 1.48 | +| `z{0}`, `me_frame=[3]` | 1366 +- 3.1 | 2193 +- 12 | 1.61 | + +Both K factors are the O(1.5) expected for Z+jet at `muR=muF=mZ` with +`ptj>30`, and the polarised one is close to but not equal to the unpolarised +one — which is what a correct polarised NLO should look like. A frame bug of +the kind M3 hit would not have landed anywhere near here. + +Two further processes, same settings: + +| process | `me_frame` | P dirs | `check_poles` | `test_ME` | xsec | +|---|---|---|---|---|---| +| `p p > z{0} z{0} [QCD]` | `[3]` | 4 | 20/20 each | passed | 1.217e+00 +- 6.6e-03 pb | +| `p p > z{0} z{0} j [QCD]` | `[3]` | 12 | 20/20 each | passed | 4.779e-01 +- 3.3e-03 pb | + +The second is the process this document names as the M3 acceptance test. + +#### What the guard now accepts, and what it still refuses + +A polarised massive particle needs a frame. There are **two** ways to have +one, and the guard now recognises both. + +*Frame from the run_card.* Modes `loonly`, `real` and `all` — the last being +what the parser calls `[QCD]` with no mode keyword (`LoopOption='all'`, +`madgraph_interface.py:5104`) — **for a purely QCD perturbation**. These are +the modes whose matrix elements this work boosts. + +*Frame from the user.* Mode `virt` outputs standalone MadLoop: the caller +supplies the phase-space point, so the frame is theirs to choose and there is +no `me_frame` and nothing for the code to get wrong. Accepted for **any** +perturbation orders, since the frame argument does not depend on them. This +was initially refused on the grounds that "it has no run_card, so it has no +frame" — which had the implication backwards. Confirmed by running it: +`p p > z{0} j [virt=QCD]` does reject `order=NLO`, `fixed_order=ON` and +`set me_frame [3]` at the launch prompt, because it is a standalone check and +not an integrable process directory. That is the point, not a problem. + +**Unblocking `virt` exposed a separate bug, which had to be fixed for the +unblock to mean anything.** `p p > z{0} j [virt=QCD]` was generating +`P0_gu_zu` — no `z0` — with 24 helicity configurations covering all three Z +polarisations. The user would have asked for a polarised loop and silently +received the unpolarised one. + +The polarisation is dropped by the **multiparticle expansion**, and only on +the MadLoop path. `u u~ > z{0} g [virt=QCD]` (explicit flavours) is correct — +`P0_uux_z0g`, 8 configurations, Z pinned to `0`. Tree level is correct with +multiparticles too. The culprit is `ProcessDefinition.__iter__` +(`base_objects.py:4670`), whose own docstring says *"not used by MG which +used some smarter version (use by ML)"*: it rebuilds each expanded leg as +`Leg({'id':id, 'state':…})` and never carries `polarization` across. +`loop_interface.py:885` iterates it and re-issues each combination as a text +command, so by the time `nice_string` runs the `{0}` is already gone. +`MultiProcess.generate_multi_amplitudes` — the path MG5 itself uses — copies +`polarization` properly, which is why only the loop side was affected. + +Fixed by carrying the multileg's `polarization` onto each expanded leg. +`p p > z{0} j [virt=QCD]` now gives `P0_gu_z0u` with 8 configurations in all +seven subprocess MEs; unpolarised `p p > z j [virt=QCD]` still gives 24, so +nothing else moves. + +Validated by running the standalone check, which has its own version of the +check_poles argument: the IR pole coefficient is a colour factor times the +Born, so `1eps/(born*ao2pi)` must be the same whether or not the helicities +are restricted, while the Born itself must differ. + +| `g u > z u` | Born | `1eps/(born*ao2pi)` | +|---|---|---| +| `z{0}` | 3.1963946124086368e-04 | -4.1666666666641756 | +| unpolarised | 0.55343246189280271 | -4.1666666666665897 | + +Both give -25/6. Had the loop summed all helicities while the Born stayed +restricted, this ratio would have been wrong — the same signature that gave +M3 away. + +**The sibling defect, and why it has to land with this change.** +`ProcessDefinition.get_process` (`base_objects.py:4849`) hardcodes +`'polarization':[]` the same way. Verified on this branch: + + multileg pol : [[], [], [0], []] + get_process : [(2, []), (-2, []), (23, []), (21, [])] shell: 0_uux_zg + +Its callers are the `check` command paths. That matters here because +`check_check` runs `check_process_format` (`madgraph_interface.py:1156`), so +**opening `virt` in the guard also newly admits +`check timing|stability|profile [virt=QCD]`**, +which was refused before. Without the `get_process` fix that check silently +runs the *unpolarised* loop. Tree-level `check` is unaffected by the guard — +it has no brackets, so the polarisation branch never fires — and its version +of the bug is purely pre-existing. + +Fixed separately on branch `claude/quizzical-ptolemy-16f0c8`, commit +`2b151185d`, which also found a **third** instance: +`process_checks.py:2032` `run_multiprocs_no_crossings` built its legs with no +polarization at all, and that is the function every tree-level `check` goes +through — so fixing `get_process` alone would have left `check +gauge/lorentz/permutation` unpolarised. Their evidence, `check permutation +p p > z{…} j --seed=1`: + +| | value | +|---|---| +| `g Q > z{0} Q` | 3.1783208585e-04 | +| `g Q > z{T} Q` | 1.0612625837e-01 | +| `g Q > z Q` (unpolarised) | 1.0644409045e-01 | + +the two polarised pieces summing exactly to the unpolarised one. + +**So `2b151185d` should land together with this change, not after it.** The +guard here opens the door; that commit is what makes what is behind it +correct. Its loop-side caller (`process_checks.py:2111`, +`generate_loop_matrix_element` — the `check timing/stability/profile` path) +is a pure data carry but is **unverified at runtime** on either branch, and +it is exactly the combination this guard newly admits. Worth one polarised +`check timing` run before relying on it. + +Two further things to know before using `check` on these processes, neither +of them caused by this work: + +- `run_multiprocs_no_crossings` builds its legs with no `flavor` key, so + `_flavor_group_pos` returns 1 for every leg and **only flavour index 1 is + actually gauge-checked** once flavour grouping is on (the default). A + passing `check` covers less than the summary line suggests; + `set apply_flavor_grouping False` restores the per-flavour table. +- FD/axial gauge on a process with a cross-base FLV routine (`FFVx_FFSy`, + classically massive-lepton Goldstone-Yukawa) was broken in every backend + until `6a8530f84` on `claude/bold-mestorf-777ea7` — the call site emitted + `FFV2_FFS3_0` against ALOHA's `FFV2_FFS3M_0`. It does **not** affect the + processes in this document: grepping the generated `p p > z{0} j [QCD]` + tree finds no cross-base combined routines at all, so `tags` is empty and + the fix is a no-op here. Checked rather than assumed. + +Still refused, each for a reason rather than by omission: + +- `[QED QCD]` and `[QED]` in the boosted modes. The frame wrappers are + order-agnostic and the QED counterterms (`particle_charge`, the charge + links, `sreal_deg`) go through the same `sborn_frame`/`sborn_sf_frame`, so + this is likely to work — but nothing in the QED sector has been run, and + the plan does not cover it. `[virt=QED]` is fine, being standalone. +- `[noborn=QCD]` with a *massive* polarised particle. Pre-existing and + unchanged: loop-induced is exported through the LO madevent template and so + already has `me_frame`, but `check()` has always rejected massive there. + Worth revisiting; not touched here because nothing in this work exercises + it. See section 1. + +`tests/unit_tests/interface/test_cmd.py::test_check_generate` encoded the old +restriction (`'u u~ > w+{L} [QCD]'` in the *invalid* list) and is updated: the +massive-colourless QCD forms and both `[virt=…]` spellings move to the +accepted list, `[QED QCD]` and `[QED]` stay in the rejected one. It was that +test that caught the first version of this change quietly opening mixed-order +NLO. + +#### What this milestone cost, and why + +The instrumentation-over-reasoning lesson from M2 held again, and more +sharply. Every hypothesis on the list — the helicity restriction, +`ImprovePSPoint`, `NRotations`, quadruple-precision promotion — was internal +to MadLoop, because the one misread number said "MadLoop is computing +something else". The bug was one line above the call, in the driver, in code +this project had already touched. Both times the search went to the most +intricate component available, and both times the answer was in the plumbing. + +Original plan for this milestone: + + +Boost `p_born` before `Call BinothLHA(p_born,…)` (`fks_singular.f:7087`). +MadLoop already handles polarised MEs (`loop_exporters.py:1552` restores the +helicity-averaging factor). Note MadLoop has its own stability rescue that +re-evaluates in rotated/boosted frames — check it does not fight the imposed +frame. + +**Gate:** `check_poles` must still pass (the poles are proportional to the +Born, so a Born/virtual frame mismatch shows up as a pole mismatch). Then +`p p > z{0} z{0} j [QCD]` end to end. + +### M4 — Unblock and document — **DONE** + +`help_polarization` now covers fixed-order NLO: which modes accept a polarised +massive particle and where each gets its frame. Two things it previously got +wrong are corrected -- `nhel=1` was presented as a requirement when it is only +a variance choice (measured: 1369 +- 1.7 against 1370 +- 1.7), and the +QCD-only restriction on the run_card modes was not stated at all, so a reader +would have assumed `[QED]` worked. It also warns that a frame must be built +from final-state particles, an initial-state one not being infrared safe. + +The acceptance test `test_polarised_nlo_me_frame` covers +`p p > z{0} z{0} j [QCD]` with `me_frame=[3,4]` and asserts on `check_poles` +rather than on a cross-section, since the poles are proportional to the Born +and are what a frame mismatch actually breaks. Run and verified to do real +work, not to pass vacuously: 12 P dirs, 12 `check_poles.log`, zero +miscancellations, `generate p p > z{0} z{0} j [QCD]` and `[3, 4] = me_frame` +read back from the generated cards, 3.084e-01 +- 3.0e-03 pb. + +Wired into CI as `acceptancetest_105` in +`.github/workflows/acceptancetest.yml`, modelled on the existing +`test_calculate_xsect_nlo` job. Without that the test would never have run: +the workflow invokes tests by explicit name. + + +Done in M3: the guard now accepts `loonly`, `real` and `all`, deriving the +mode the same way the parser does (no keyword in the brackets means `all`), +and `check()` no longer rejects a massive particle in those modes. The stale +commented-out `order != 'qcd'` guard is deleted, along with the dead `order` +it parsed. + +Remaining: + +- `help_polarization` (`:803`) does not say `me_frame` is LO-only, but it + does not say it works at fixed-order NLO either — add that, and say which + modes accept it; +- the acceptance test in section 4 is written up but not yet added to + `tests/acceptance_tests/test_cmd_amcatnlo.py`; +- decide whether `[noborn=QCD]` should accept a massive polarised particle. + It is loop-induced, exported through the LO madevent template, so it has + `me_frame` already; `check()` has always refused massive there and nothing + in this work exercises it. + +## 4. Acceptance test + +One new test in `tests/acceptance_tests/test_cmd_amcatnlo.py`, modelled on +`MECmdShell.generate` + `self.do('calculate_xsect LO -f')` (`:722`), with the +reference LO numbers produced the way `tests/acceptance_tests/test_cmd_madevent.py:2870` +does it. + +| stage | process | mode | assertion | +|---|---|---|---| +| M0 | `p p > z j [QCD]` unpolarised | NLO | xsec unchanged vs. today (regression) | +| M1 | `p p > z{0} j` | `[LOonly=QCD]` vs LO madevent | agree within MC error | +| M1 | `p p > z{0} z{0}` | `[LOonly=QCD]` | `me_frame`-independent (null test) | +| M2 | `p p > z{0} z{0}` | `[real=QCD]` | real boost: `test_ME` passes; xsec != partonic-CM value. **Weak** — see below | +| M2 | `p p > z{0} j` | `[real=QCD]` | **discriminating**: `test_ME` on gluon-initiated channels | +| M3 | `p p > z{0} j` | `[QCD]` | **run**: `check_poles` 20/20 in all 12 P dirs, 2.193e+03 +- 1.2e+01 pb | +| M3 | `p p > z{0} z{0}` | `[QCD]` | **run**: `check_poles` 20/20 in all 4 P dirs, 1.217e+00 +- 6.6e-03 pb | +| M3 | `p p > z{0} z{0} j` | `[QCD]` | **run**: `check_poles` 20/20 in all 12 P dirs, 4.779e-01 +- 3.3e-03 pb | + +**Which process tests what.** The two are complementary, but not in the obvious +way: + +- `p p > z{0} j` — the Z recoils against the jet already at Born level, so the + `me_frame` boost is non-trivial and non-longitudinal in *both* the Born and + the counter-event, including in the singular region. This is the process that + discriminates the azimuthal phase, and it doubles as the M1 cross-check + against LO madevent. +- `p p > z{0} z{0}` — the Born boost is the identity (the ZZ system *is* the + partonic CM), so it is a clean null test at M1 and exercises the real boost in + the finite region at M2. In the singular region the boost degenerates to + longitudinal, so it will **not** catch a wrong azimuthal phase. + +Keep M1/M2 in the always-run set (fast, no MadLoop) and gate M3 behind the slow +marker. Run via `./tests/test_manager.py`, not pytest. + +## 5. Open blockers + +- **B1 — `xij_aor` cannot be recomputed after a boost.** See M2 step 4. Drives + the `k_perp_dir` design. *Confirmed.* +- **B2 — `me_frame` leg remapping** between Born and real multiplicities (D2, + plus the `skip` logic at `genps_fks.f:1206-1222`, `:1542-1553`). + Sub-claim resolved: the `2**n` vs `2**(n-1)` convention is **not** a bug — + `mapid` is `ids(i)=btest(id,i)` (`cluster.f:141`), which matches + `sum(2**n)` at `banner.py:4869`. The comment at `genps.f:1758` saying + `sum 2**(N-1)` is stale; fix the comment. +- **B3 — the `R_y(pi)` flip must be deleted, not boosted**, and + `montecarlocounter.f:3000-3006` must not keep double-encoding it. +- **B4 — RESOLVED, both parts.** FSR counter-events do collapse onto the + `vtiny` branch, exactly like ISR (they are handed the literal `one` for + `y_ij_fks`). And `shybst = O(1-y) -> 0` in the FSR collinear limit, because + it vanishes iff the mother is massless and `m_mother^2 = 2 E_i E_j (1-y)`. + So FSR inherits the D3 equivariance property. See the M2 step 0 section. +- **B5 — RESOLVED structurally.** Rather than argue about which callers are + reachable, every ME entry point was routed through a frame wrapper, making + the boost a pure function of the momenta passed. Callers that share momenta + then agree by construction. Worth stressing that this was not a theoretical + hazard: in M1 a partial conversion produced a silent 3% shift once and the + `momenta not the same in Born` stop once. See M2 step 1. +- **B6 — no guard against a lightlike/null `me_frame` sum** in `boost_to_frame` + (`Template/LO/SubProcesses/genps.f:1761`). Pre-existing at LO; not replicated + into the NLO port: `get_me_frame_boost` stops with a diagnostic instead. +- **B7 — RESOLVED.** `add_write_info.f:278` now goes through `sborn_frame` + like everything else, so polarised event generation cannot write weights + computed in a different frame from the cross section. +- **B8 — RESOLVED, and it was the M2 blocker.** A one-leg `me_frame` left that + leg at rest only to rounding, and HELAS `vxxxxx` picks the quantisation axis + from the momentum direction unless the 3-momentum is *exactly* zero. Fixed in + both `boost_to_me_frame` (NLO) and `boost_to_frame` (LO). See M2 step 2. +- **B9 — mitigated, not eliminated.** A boost whose velocity degenerates to + zero tilts the beams just enough to wreck `sqrt(p(0)+p(3))` in the HELAS + massless spinors, without reaching the exact on-axis branch. Guarded by + `betamin=1d-4` in `get_me_frame_boost`. The underlying fragility lives in + HELAS and is not fixed; anything else that boosts an on-axis beam by a tiny + amount will hit it. See M2 step 2. +- **B12 — RESOLVED here, and it has a twin that must land with it.** + `ProcessDefinition.__iter__` dropped `polarization` when expanding + multiparticles, so any polarised MadLoop-standalone process written with + `p`/`j` came back silently unpolarised. Pre-existing, exposed by unblocking + `[virt=…]`. Fixed in `base_objects.py` on this branch. The same defect sits + in `ProcessDefinition.get_process` and in + `process_checks.run_multiprocs_no_crossings`, both feeding the `check` + command — fixed on `claude/quizzical-ptolemy-16f0c8` (`2b151185d`), which + **must be merged alongside this change**, because opening `virt` in the + guard is what makes `check … [virt=QCD]` reachable for a polarised massive + particle in the first place. See M3. +- **B10 — RESOLVED, and it was the M3 blocker.** `check_poles.f` had been + half-converted: its Born line went through `sborn_frame` while the + `BinothLHA` call one line below still got the unboosted `p_born`, so the + gate compared MadFKS poles in the me_frame against MadLoop poles in the + partonic c.m.. This is B5 all over again — two consumers of one kinematic + configuration disagreeing about the frame — but in a *driver* rather than + in the integration, which is why the M2 sweep's "route every ME entry + point through a wrapper" rule did not catch it: `BinothLHA` is not an ME + entry point, it is a wrapper of one. Fixed at `check_poles.f:228`. See M3. +- **B11 — latent, documented, not fixed.** `NRotations_DP`/`NRotations_QP` + re-evaluate the loop at a rotated phase-space point and expect `|M|^2` + unchanged. That fails for a particle put at rest by a one-leg `me_frame`, + whose quantisation axis is the frame z axis and does not rotate with the + momenta, so MadLoop would report spurious instability. Both default to 0, + so nothing is wrong today; a user who raises them on a polarised run will + see it. Same root as B8/B9: the HELAS branch at exactly zero momentum. + +## 6. Risk + +M2 step 4 is the only genuinely hard item; everything else is mechanical. D2 is +the decision most likely to need revisiting, and M1 proves it out cheaply. diff --git a/madgraph/core/base_objects.py b/madgraph/core/base_objects.py index a43e58fad7..c1f77d5b32 100755 --- a/madgraph/core/base_objects.py +++ b/madgraph/core/base_objects.py @@ -4675,22 +4675,36 @@ def __iter__(self): if leg['state'] == False] fsids = [leg['ids'] for leg in self['legs'] \ if leg['state'] == True] + # The polarization restriction lives on the multileg and has to be + # carried onto each expanded leg, or a polarized multiparticle + # process silently comes back unpolarized. MultiProcess. + # generate_multi_amplitudes (diagram_generation.py) does the same for + # the path MG5 itself uses. + ispols = [leg['polarization'] for leg in self['legs'] \ + if leg['state'] == False] + fspols = [leg['polarization'] for leg in self['legs'] \ + if leg['state'] == True] red_isidlist = [] # Generate all combinations for the initial state for prod in itertools.product(*isids): - islegs = [Leg({'id':id, 'state': False}) for id in prod] + # list() so that every yielded process owns its polarization + # rather than sharing the definition's list object + islegs = [Leg({'id':id, 'state': False, 'polarization': list(pol)}) + for id, pol in zip(prod, ispols)] if tuple(sorted(prod)) in red_isidlist: - continue - red_isidlist.append(tuple(sorted(prod))) + continue + red_isidlist.append(tuple(sorted(prod))) red_fsidlist = [] for prod in itertools.product(*fsids): # Remove double counting between final states if tuple(sorted(prod)) in red_fsidlist: - continue + continue red_fsidlist.append(tuple(sorted(prod))) leg_list = [copy.copy(leg) for leg in islegs] - leg_list.extend([Leg({'id':id, 'state': True}) for id in prod]) + leg_list.extend([Leg({'id':id, 'state': True, + 'polarization': list(pol)}) + for id, pol in zip(prod, fspols)]) legs = LegList(leg_list) process = self.get_process_with_legs(legs) yield process @@ -4834,21 +4848,27 @@ def get_process(self, initial_state_ids, final_state_ids): """ Return a Process object which has the same properties of this ProcessDefinition but with the specified given leg ids. """ + # The polarization is not passed in argument, it has to be read from the + # multi-legs of this process definition, which are ordered in the same + # way as the ids given in argument (initial states first). + my_islegs = [leg for leg in self.get('legs') if not leg.get('state')] + my_fslegs = [leg for leg in self.get('legs') if leg.get('state')] + # First make sure that the desired particle ids belong to those defined # in this process definition. if __debug__: - my_isids = [leg.get('ids') for leg in self.get('legs') \ - if not leg.get('state')] - my_fsids = [leg.get('ids') for leg in self.get('legs') \ - if leg.get('state')] for i, is_id in enumerate(initial_state_ids): - assert is_id in my_isids[i] + assert is_id in my_islegs[i].get('ids') for i, fs_id in enumerate(final_state_ids): - assert fs_id in my_fsids[i] - + assert fs_id in my_fslegs[i].get('ids') + return self.get_process_with_legs(LegList(\ - [Leg({'id': id, 'state':False, 'polarization':[]}) for id in initial_state_ids] + \ - [Leg({'id': id, 'state':True, 'polarization':[]}) for id in final_state_ids])) + [Leg({'id': id, 'state':False, + 'polarization': list(leg.get('polarization'))}) \ + for id, leg in zip(initial_state_ids, my_islegs)] + \ + [Leg({'id': id, 'state':True, + 'polarization': list(leg.get('polarization'))}) \ + for id, leg in zip(final_state_ids, my_fslegs)])) def __eq__(self, other): """Overloading the equality operator, so that only comparison diff --git a/madgraph/fks/fks_base.py b/madgraph/fks/fks_base.py index 3b6f941d39..f9faaab001 100755 --- a/madgraph/fks/fks_base.py +++ b/madgraph/fks/fks_base.py @@ -599,6 +599,8 @@ def __init__(self, start_proc = None, remove_reals = True, ncores_for_proc_gen=0 self.extra_cnt_amp_list = diagram_generation.AmplitudeList() self.ncores_for_proc_gen = ncores_for_proc_gen self.sudakov_amps = [] + # user number of each FKS-sorted born leg; see fks_common.get_user_leg_order + self.user_leg_order = [] if not remove_reals in [True, False]: raise fks_common.FKSProcessError(\ @@ -610,12 +612,18 @@ def __init__(self, start_proc = None, remove_reals = True, ncores_for_proc_gen=0 pertur = start_proc['perturbation_couplings'] if pertur: self.perturbation = sorted(pertur)[0] + # record the user's leg numbering before sort_proc discards it + self.user_leg_order = fks_common.get_user_leg_order(\ + start_proc, pert = self.perturbation) self.born_amp = diagram_generation.Amplitude(\ copy.copy(fks_common.sort_proc(\ start_proc, pert = self.perturbation))) #initialize with an amplitude elif isinstance(start_proc, diagram_generation.Amplitude): pertur = start_proc.get('process')['perturbation_couplings'] + self.user_leg_order = fks_common.get_user_leg_order(\ + start_proc['process'], + pert = self.perturbation) self.born_amp = diagram_generation.Amplitude(\ copy.copy(fks_common.sort_proc(\ start_proc['process'], diff --git a/madgraph/fks/fks_common.py b/madgraph/fks/fks_common.py index 7d055c6a0b..0ed83af936 100755 --- a/madgraph/fks/fks_common.py +++ b/madgraph/fks/fks_common.py @@ -721,6 +721,26 @@ def legs_to_color_link_string(leg1, leg2, pert = 'QCD'): #test written, all case return dict +def get_user_leg_order(process, pert = 'QCD'): + """Return, for each position of the FKS-sorted leg list, the number the + leg carried in the process as the user wrote it. + + sort_proc() both permutes and renumbers the legs, so after it has run the + user's numbering is gone for good. me_frame in the run_card is expressed + in the user's numbering, so this mapping is what lets the fortran get back + to it; it is written out as frame_map_born in frame_info.inc. + + The returned list is 0-indexed by FKS position: entry i is the user number + of the leg sitting at FKS position i+1. + + Call this *before* sort_proc. to_fks_legs() copies the legs, so this does + not disturb the process it is given. + """ + leglist = to_fks_legs(process.get('legs'), process.get('model')) + leglist.sort(pert = pert) + return [leg['number'] for leg in leglist] + + def sort_proc(process,pert = 'QCD'): """Given a process, this function returns the same process but with sorted FKSLegs. diff --git a/madgraph/fks/fks_helas_objects.py b/madgraph/fks/fks_helas_objects.py index 89ae3112a6..389fea409a 100755 --- a/madgraph/fks/fks_helas_objects.py +++ b/madgraph/fks/fks_helas_objects.py @@ -693,7 +693,9 @@ def __init__(self, fksproc=None, real_me_list =[], real_amp_list=[], self.real_processes = [] self.extra_cnt_me_list = [] self.perturbation = fksproc.perturbation - self.charges_born = fksproc.get_charges() + self.charges_born = fksproc.get_charges() + # user numbering of the born legs, for me_frame (frame_info.inc) + self.user_leg_order = fksproc.user_leg_order real_amps_new = [] for extra_cnt in fksproc.extra_cnt_amp_list: diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py index d7e1ef7fd0..ef22fbd400 100755 --- a/madgraph/interface/madgraph_interface.py +++ b/madgraph/interface/madgraph_interface.py @@ -807,8 +807,25 @@ def help_polarization(self): logger.info(" > Example: generate t{L} > w+{T} b{R}, w+ > ta+ vt",'$MG:color:GREEN') logger.info(" > Example: generate p p > z{T} z{A}, z > e+ e-",'$MG:color:GREEN') logger.info(" > Example: generate p p > z{0} z{T}, z > e+ e-, z > mu+ mu-",'$MG:color:GREEN') - logger.info(" > Users need to set 'group_subprocesses False', 'nhel=1' (run_card), and 'me_frame' (run_card)") + logger.info(" > At LO, users need to set 'group_subprocesses False' and 'me_frame' (run_card).") + logger.info(" 'nhel=1' is only a variance choice, not a requirement: it selects Monte-Carlo") + logger.info(" over helicities, which is an unbiased estimator of the same cross-section.") logger.info(" > For the proces 'p p > w+ z j j, w+ > l+ vl, z > l+ l-', the WZ rest frame is given by me_frame = [3,4,5,6]") + logger.info("Polarization at fixed-order NLO:",'$MG:BOLD') + logger.info(" > A polarized massive particle needs a frame, so it is allowed for the NLO modes") + logger.info(" that have one:") + logger.info(" - [QCD], [real=QCD] and [LOonly=QCD] take 'me_frame' from the run_card.") + logger.info(" QCD is the only perturbation order supported here so far: the boost") + logger.info(" machinery is order-agnostic, but nothing in the QED sector has been") + logger.info(" validated, so [QED] and mixed orders are still refused.") + logger.info(" - [virt=QCD] and standalone MadLoop take the frame from the momenta the") + logger.info(" caller supplies, so any perturbation order is allowed.") + logger.info(" A polarized massless particle needs no frame and was always allowed.") + logger.info(" > Example: generate p p > z{0} z{0} j [QCD] with me_frame = [3,4]",'$MG:color:GREEN') + logger.info(" > Define the frame from final-state particles only. A frame built from the") + logger.info(" initial state is not infrared safe at NLO -- the real emission and the reduced") + logger.info(" Born carry different momentum fractions even in the collinear limit -- and is") + logger.info(" refused. me_frame = [1,2] names the partonic c.m. and is simply skipped.") logger.info(" > For further details, see appendices of [arXiv:1912.01725] and [arXiv:2512.10015],") logger.info(" and for possibilities with loop-induced processes, see [2401.17365].") @@ -1238,24 +1255,52 @@ def check_process_format(self, process): raise self.InvalidCmd('Polarization restriction can not be used in forbidding particles') if '[' in process and '{' in process: - valid = False - if 'noborn' in process or 'sqrvirt' in process: - valid = True + # Which NLO mode was asked for, taken from inside the brackets so + # that a stray 'real' elsewhere in the process cannot match. + bracket = process.split('[')[1].split(']')[0] + if '=' in bracket: + nlo_mode, pert_orders = bracket.split('=', 1) + nlo_mode = nlo_mode.strip().lower() + else: + # no keyword, e.g. '[QCD]': the parser calls that mode 'all' + nlo_mode, pert_orders = 'all', bracket + + # A polarized massive particle needs a frame to be defined in. + # There are two ways to have one. + # + # 'virt' outputs standalone MadLoop, where the user supplies the + # phase-space point themselves and so chooses the frame. There is + # no run_card, no me_frame and nothing for the code to get wrong, + # whatever the perturbation orders. + standalone_olp = nlo_mode == 'virt' + + # Otherwise the frame comes from me_frame in the run_card, and + # every piece of the computation is boosted into it: the Born + # (M1), the real emission and the full set of FKS counterterms + # (M2), and the virtual, which goes through binothlha_frame (M3). + # check_poles cancels 20/20 in every P dir of p p > z{0} j [QCD] + # with me_frame=[3]. See docs/nlo_polarisation_boost_plan.md. + # + # Restricted to QCD corrections: the frame wrappers are + # order-agnostic and the QED counterterms go through them too, so + # this is likely to work, but nothing in the QED sector has been + # validated, so mixed or pure-QED perturbation stays refused. + frame_supported = standalone_olp or \ + (nlo_mode in ('loonly', 'real', 'all') + and pert_orders.strip().lower() == 'qcd') + if 'noborn' in process or 'sqrvirt' in process \ + or frame_supported: + pass else: raise self.InvalidCmd('Polarization restriction can not be used for NLO processes') - # below are the check when [QCD] will be valid for computation - order = process.split('[')[1].split(']')[0] - if '=' in order: - order = order.split('=')[1] -# if order.strip().lower() != 'qcd': -# raise self.InvalidCmd('Polarization restriction can not be used for generic NLO computations') - def check(p): if p.get('color') != 1: raise self.InvalidCmd('Polarization restriction can not be used for color charged particles') - elif p.get('mass') != 'ZERO': - raise self.InvalidCmd('Polarization restriction can not be used for massive particles') + elif p.get('mass') != 'ZERO' and not frame_supported: + # massive polarization needs a frame; only the modes whose + # matrix elements are boosted have one + raise self.InvalidCmd('Polarization restriction can not be used for massive particles') diff --git a/madgraph/iolibs/export_fks.py b/madgraph/iolibs/export_fks.py index eb45248cbe..0e5281ca51 100755 --- a/madgraph/iolibs/export_fks.py +++ b/madgraph/iolibs/export_fks.py @@ -608,6 +608,10 @@ def generate_directories_fks(self, matrix_element, fortran_model, me_number, matrix_element, fortran_model) + filename = 'frame_info.inc' + self.write_frame_info_file(writers.FortranWriter(filename), + matrix_element) + filename = 'maxconfigs.inc' self.write_maxconfigs_file(writers.FortranWriter(filename), max(nconfigs,matrix_element.born_me.get_number_of_amplitudes())) @@ -711,6 +715,7 @@ def generate_directories_fks(self, matrix_element, fortran_model, me_number, 'weight_lines.f', 'genps_fks.f', 'boostwdir2.f', + 'boost_to_frame.f', 'madfks_mcatnlo.inc', 'open_output_files.f', 'open_output_files_dummy.f', @@ -1748,6 +1753,36 @@ def write_maxparticles_file(self, writer, maxparticles): writer.writelines(lines) + def write_frame_info_file(self, writer, matrix_element): + """Write frame_info.inc, the bridge between the me_frame entries of the + run_card and the leg positions the fortran actually uses. + + me_frame is given in the numbering of the process as the user wrote it, + but sort_proc() reorders and renumbers the born legs, so the two do not + coincide. frame_map_born(i) is the user number of the born leg sitting + at position i. + + No table is needed for the real emissions: the FKS convention fixes the + real->born correspondence in terms of i_fks alone (see set_pdg in + chooser_functions.f), so get_frame_mask_real derives it at runtime. + """ + born_legs = matrix_element.born_me.get('processes')[0].get('legs') + nexternal_born = len(born_legs) + + user_order = getattr(matrix_element, 'user_leg_order', []) + if len(user_order) != nexternal_born: + # No recorded ordering (e.g. a matrix element built by hand in the + # tests). Fall back to the identity: correct whenever the FKS sort + # left the user's ordering alone, which is the common case. + user_order = list(range(1, nexternal_born + 1)) + + lines = [] + lines.append("integer frame_map_born(nexternal-1)") + lines.append("data frame_map_born /%s/" % \ + ','.join('%d' % n for n in user_order)) + writer.writelines(lines) + + def write_maxconfigs_file(self, writer, maxconfigs): """Write the maxconfigs.inc file for MadEvent""" lines = "integer lmaxconfigs\n" diff --git a/madgraph/various/banner.py b/madgraph/various/banner.py index 5d5b0c63dd..1296516a87 100755 --- a/madgraph/various/banner.py +++ b/madgraph/various/banner.py @@ -5760,7 +5760,7 @@ class RunCardNLO(RunCard): LO = False - blocks = [heavy_ion_block, running_block_nlo] + blocks = [heavy_ion_block, frame_block, running_block_nlo] dummy_fct_file = {"dummy_cuts": pjoin("SubProcesses","dummy_fct.f"), "user_dynamical_scale": pjoin("SubProcesses","dummy_fct.f"), @@ -5866,6 +5866,10 @@ def default_setup(self): self.add_param('systematics_program', 'none', include=False, hidden=True, comment='Choose which program to use for systematics computation: none, systematics') self.add_param('systematics_arguments', [''], include=False, hidden=True, comment='Choose the argment to pass to the systematics command. like --mur=0.25,1,4. Look at the help of the systematics function for more details.') + #frame in which to evaluate the matrix-element (polarization) + self.add_param("me_frame", [1,2], hidden=True, include=False, comment="choose lorentz frame where to evaluate the matrix-element [for non lorentz invariant matrix-element/polarization]:\n the entries are the leg numbers of the process as written by the user; the rest-frame of their momentum sum is used.\n [1,2] (the initial state) and the full final state both mean the partonic center of mass, and are skipped rather than applied.\n Define the frame from final-state particles only: a frame built from the initial state is not infrared safe at NLO.") + self.add_param('frame_id', 6, system=True) + #technical self.add_param('folding', [1,1,1], include=False) @@ -5924,9 +5928,28 @@ def default_setup(self): def check_validity(self): """check the validity of the various input""" - + super(RunCardNLO, self).check_validity() + # me_frame built out of initial-state legs is not infrared safe at NLO: + # the real emission and the reduced Born carry different momentum + # fractions, by a finite amount even in the singular limit, so the + # frame jumps across that limit and the subtraction stops cancelling. + # Selecting exactly the initial state (or, equivalently, exactly the + # whole final state) is the partonic c.m. and is simply skipped; any + # other use of an initial-state leg is a genuine mistake. + if 'me_frame' in self.user_set: + initial = [n for n in self['me_frame'] if n in (1, 2)] + if initial and len(self['me_frame']) > len(initial): + raise InvalidRunCard( + 'me_frame %s mixes initial-state legs with final-state ' + 'ones. A frame defined using the initial state is not ' + 'infrared safe at NLO: the real emission and the reduced ' + 'Born have different momentum fractions even in the ' + 'collinear limit, so the frame is discontinuous there. ' + 'Define the frame from final-state particles only.' + % self['me_frame']) + # if heavy ion mode use for one beam, forbid lpp!=1 if self['lpp1'] not in [1,2]: if self['nb_proton1'] !=1 or self['nb_neutron1'] !=0: @@ -6145,6 +6168,26 @@ def add_citation(self, cite_fn): def update_system_parameter_for_include(self): + # polarization: rest-frame in which to evaluate the matrix-element. + # Same encoding as at LO (see mapid in cluster.f): bit n of frame_id is + # set for each leg n listed in me_frame. + # + # frame_id=0 selects no leg at all, which the fortran reads as "skip + # the boost". That is the right default here, and it is not the same + # as the LO default: at LO the momenta reach the matrix element in the + # lab frame, so me_frame=[1,2] is a real boost to the partonic c.m., + # whereas MadFKS already works in a frame close to it. Close, but not + # equal -- the real emission lives in a frame boosted along z, since + # its two initial momenta carry different energies -- so honouring + # [1,2] literally would apply a longitudinal boost to every unpolarised + # run. That is an identity for |M|^2 but not bit for bit, and it is + # enough to send the adaptive grids down a different path. So the boost + # runs only when a frame was actually asked for. + if 'me_frame' in self.user_set: + self['frame_id'] = sum(2**(n) for n in self['me_frame']) + else: + self['frame_id'] = 0 + # set the pdg_for_cut fortran parameter pdg_to_cut = set(list(self['pt_min_pdg'].keys()) +list(self['pt_max_pdg'].keys())+ list(self['mxx_min_pdg'].keys())+ list(self['mxx_only_part_antipart'].keys())) @@ -6250,7 +6293,21 @@ def create_default_for_process(self, proc_characteristic, history, proc_def): # If model has running functionality add the additional parameter model = proc_def[0].get('model') if model['running_elements']: - self.display_block.append('RUNNING') + self.display_block.append('RUNNING') + + # if polarization is used, expose the choice of the frame in the run_card. + # Only needed for massive particles: for massless ones the helicity is + # boost invariant along the momentum, so the frame does not matter. + for proc in proc_def: + for l in proc.get('legs'): + if l.get('polarization'): + particle = proc.get('model').get_particle(l.get('id')) + if particle.get('mass').lower() != 'zero': + self.display_block.append('frame') + break + else: + continue + break # Check if need matching min_particle = 99 diff --git a/madgraph/various/process_checks.py b/madgraph/various/process_checks.py index 423ec5d603..f2d37d9d79 100755 --- a/madgraph/various/process_checks.py +++ b/madgraph/various/process_checks.py @@ -2030,11 +2030,8 @@ def run_multiprocs_no_crossings(function, multiprocess, stored_quantities, multiprocess, model, id_anti_id_dict): continue # Generate process based on the selected ids - process = multiprocess.get_process_with_legs(base_objects.LegList(\ - [base_objects.Leg({'id': id, 'state':False}) for \ - id in is_prod] + \ - [base_objects.Leg({'id': id, 'state':True}) for \ - id in fs_prod])) + # (get_process carries over the polarization of the multi-legs) + process = multiprocess.get_process(is_prod, fs_prod) if opt is not None: if isinstance(opt, dict): diff --git a/tests/acceptance_tests/test_cmd.py b/tests/acceptance_tests/test_cmd.py index 0f9a64b0b5..2a94943f9f 100755 --- a/tests/acceptance_tests/test_cmd.py +++ b/tests/acceptance_tests/test_cmd.py @@ -3864,7 +3864,14 @@ def test_check_gauge_epem_vevex_wpwm(self): self.assertIn('Summary: 1/1 passed, 0/1 failed', log) def test_check_pp_wpwm(self): - """Test `check p p > w+ w-` runs and gauge check succeeds.""" + """Test `check p p > w+ w-` runs and gauge check succeeds. + + With apply_flavor_grouping on (the default), the four light-quark + subprocesses are carried by the single merged matrix element + Q Qx > w+ w-, so the gauge block checks one process, not four. The + per-flavor coverage lives in the flavor-grouping block, which compares + the merged matrix element against the unmerged one for every flavor. + """ self.do('import model sm') with self.assertLogs('madgraph.check_cmd', level='DEBUG') as cm: @@ -3872,10 +3879,17 @@ def test_check_pp_wpwm(self): log = '\n'.join(cm.output) self.assertIn('Gauge results (switching between Unitary/Feynman/Axial/FD gauge):', log) - self.assertIn('Summary: 4/4 passed, 0/4 failed', log) + self.assertIn('Q Qx > w+ w-', log) + self.assertIn('Summary: 1/1 passed, 0/1 failed', log) + # the four flavors (both orderings) are still checked, here: + self.assertIn('Flavor grouping check results:', log) + self.assertIn('Summary: 8/8 passed, 0/8 failed', log) def test_check_gauge_pp_wpwm(self): - """Test `check gauge p p > w+ w-` includes axial and succeeds.""" + """Test `check gauge p p > w+ w-` includes axial and succeeds. + + See test_check_pp_wpwm for why a single merged process is checked. + """ self.do('import model sm') with self.assertLogs('madgraph.check_cmd', level='INFO') as cm: @@ -3883,7 +3897,8 @@ def test_check_gauge_pp_wpwm(self): log = '\n'.join(cm.output) self.assertIn('Gauge results (switching between Unitary/Feynman/Axial/FD gauge):', log) - self.assertIn('Summary: 4/4 passed, 0/4 failed', log) + self.assertIn('Q Qx > w+ w-', log) + self.assertIn('Summary: 1/1 passed, 0/1 failed', log) def test_check_gauge_epem_aa_includes_axial(self): """Test `check gauge e+ e- > a a` includes axial gauge and succeeds.""" diff --git a/tests/acceptance_tests/test_cmd_amcatnlo.py b/tests/acceptance_tests/test_cmd_amcatnlo.py index 54c8344120..67702a1126 100755 --- a/tests/acceptance_tests/test_cmd_amcatnlo.py +++ b/tests/acceptance_tests/test_cmd_amcatnlo.py @@ -756,6 +756,63 @@ def test_calculate_xsect_lo(self): self.assertTrue(os.path.exists('%s/Events/run_01_LO/alllogs_1.html' % self.path)) + @set_global() + def test_polarised_nlo_me_frame(self): + """Polarised fixed-order NLO in a chosen rest frame. + + p p > z{0} z{0} j [QCD] with me_frame = [3,4], i.e. the matrix + elements evaluated in the ZZ rest frame. The Born, the real, every FKS + counterterm and the virtual all have to be boosted to the *same* + frame, rebuilt from each configuration's own momenta. + + What is actually being tested is check_poles and test_soft_col_limits, + which the launch runs on its own and which fail the run if they do not + pass. They are the only checks sensitive to this: + + - the infrared poles are proportional to the Born, so any frame + mismatch between the Born and the virtual shows up as a per-point + constant ratio between the MadFKS and the OLP poles; + - the collinear Q term is the same order as the AP term, so a wrong + azimuthal phase does not cancel and the soft/collinear ratios + plateau off 1. + + A cross-section comparison would not catch either: both were wrong at + some point during development while the total stayed plausible. + """ + self.generate('p p > z{0} z{0} j [QCD]', 'loop_sm') + + # Ask for the ZZ rest frame. nn23lo1 rather than lhapdf, since the + # python lhapdf bindings are broken under some interpreters and would + # fail this test for an unrelated reason; scales fixed so the run is + # reproducible; req_acc_fo loose because the assertion is on the + # checks, not on the precision of the cross-section. + card_path = pjoin(self.path, 'Cards', 'run_card.dat') + run_card = banner.RunCardNLO(card_path) + run_card.set('me_frame', [3, 4], user=True) + run_card.set('pdlabel', 'nn23lo1', user=True) + run_card.set('fixed_ren_scale', True, user=True) + run_card.set('fixed_fac_scale', True, user=True) + run_card.set('mur_ref_fixed', 91.188, user=True) + run_card.set('muf_ref_fixed', 91.188, user=True) + run_card.set('ptj', 30.0, user=True) + run_card.set('etaj', 4.0, user=True) + run_card.set('req_acc_fo', 0.05, user=True) + run_card.write(card_path) + + self.do('calculate_xsect NLO -f') + + self.assertTrue(os.path.exists('%s/Events/run_01/summary.txt' % self.path)) + + # check_poles writes one log per P directory; none may report a + # miscancellation. This is the assertion that the frame handling is + # consistent between the Born and the virtual. + pole_logs = misc.glob(pjoin(self.path, 'SubProcesses', 'P*', + 'check_poles.log')) + self.assertTrue(pole_logs, 'check_poles did not run') + for log in pole_logs: + self.assertNotIn('MISCANCELLATION', open(log).read(), + 'poles do not cancel in %s' % log) + def test_amcatnlo_from_file(self): """ """ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_bbx_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_bbx_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_bbx_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_bxb_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_bxb_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_bxb_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_ccx_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_ccx_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_ccx_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_cxc_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_cxc_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_cxc_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_ddx_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_ddx_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_ddx_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_dxd_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_dxd_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_dxd_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_ssx_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_ssx_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_ssx_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_sxs_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_sxs_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_sxs_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_uux_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_uux_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_uux_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_uxu_zz%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_uxu_zz%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportEWSudTest/test_ppzz_ewsudakov/%SubProcesses%P0_uxu_zz%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ag_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ag_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ag_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ga_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ga_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ga_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uux_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uux_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uux_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uxu_ttx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uxu_ttx%frame_info.inc new file mode 100644 index 0000000000..2cd17818e4 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uxu_ttx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3,4/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%frame_info.inc new file mode 100644 index 0000000000..563a88365f --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%frame_info.inc new file mode 100644 index 0000000000..563a88365f --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,3/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%frame_info.inc new file mode 100644 index 0000000000..d917193a94 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,4,2,3/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%frame_info.inc new file mode 100644 index 0000000000..0c8bfad904 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,4,3/ diff --git a/tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%frame_info.inc b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%frame_info.inc new file mode 100644 index 0000000000..0c8bfad904 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%frame_info.inc @@ -0,0 +1,2 @@ + INTEGER FRAME_MAP_BORN(NEXTERNAL-1) + DATA FRAME_MAP_BORN /1,2,4,3/ diff --git a/tests/unit_tests/core/test_base_objects.py b/tests/unit_tests/core/test_base_objects.py index 65d99dce7c..8b883c7ac7 100755 --- a/tests/unit_tests/core/test_base_objects.py +++ b/tests/unit_tests/core/test_base_objects.py @@ -2166,6 +2166,29 @@ def test_get_process_with_legs(self): else: self.assertEqual(myleglist, testproc[k]) + def test_get_process_keeps_polarization(self): + """test that get_process carries the polarization of the multi-legs + over to the legs of the returned process.""" + + my_new_process_definition = copy.deepcopy(self.my_process_definition) + # one initial state and one final state leg are polarized + my_new_process_definition['legs'][1].set('polarization', [-1]) + my_new_process_definition['legs'][3].set('polarization', [0]) + + testproc = my_new_process_definition.get_process([3, 3], [4, 5, 3]) + + self.assertEqual([l.get('id') for l in testproc.get('legs')], + [3, 3, 4, 5, 3]) + self.assertEqual([l.get('state') for l in testproc.get('legs')], + [False, False, True, True, True]) + self.assertEqual([l.get('polarization') for l in testproc.get('legs')], + [[], [-1], [], [0], []]) + + # the process must not share the polarization list of the multi-leg + testproc.get('legs')[1].get('polarization').append(1) + self.assertEqual(my_new_process_definition['legs'][1].get('polarization'), + [-1]) + def test_values_for_prop(self): """Test filters for process properties""" diff --git a/tests/unit_tests/interface/test_cmd.py b/tests/unit_tests/interface/test_cmd.py index 7ef96da2fd..741a992258 100755 --- a/tests/unit_tests/interface/test_cmd.py +++ b/tests/unit_tests/interface/test_cmd.py @@ -257,7 +257,14 @@ def test_check_generate(self): cmd.check_process_format('g g > Z Z [ noborn=QCD] @1') cmd.check_process_format('u u~ > 2w+ 2j') cmd.check_process_format('u u~ > 2w+{0} 2j') - cmd.check_process_format,'u u~ > e+{L} vl [QED]' + cmd.check_process_format('u u~ > w+{L} [QCD]') + cmd.check_process_format('u u~ > z{0} g [QCD]') + cmd.check_process_format('u u~ > z{0} g [real=QCD]') + cmd.check_process_format('u u~ > z{0} g [LOonly=QCD]') + # standalone MadLoop: the user supplies the momenta, so the frame is + # theirs and the perturbation orders do not matter + cmd.check_process_format('u u~ > z{0} g [virt=QCD]') + cmd.check_process_format('u u~ > z{0} g [virt=QED QCD]') # unvalid syntax self.wrong(cmd.check_process_format, ' e+ e-') @@ -274,8 +281,9 @@ def test_check_generate(self): self.wrong(cmd.check_process_format, 'e+ e- > Z > mu+ mu- / W+{L}') self.wrong(cmd.check_process_format, 'e+ e- > Z > mu+ mu- $ W+{L}') self.wrong(cmd.check_process_format, 'u u~ > t{L} t~ [QCD]') + # massive colourless polarization at NLO QCD is supported since the + # me_frame boost reaches the virtual; mixed and pure QED are not self.wrong(cmd.check_process_format, 'u u~ > W+{L} vl [ QED QCD]') - self.wrong(cmd.check_process_format,'u u~ > w+{L} [QCD]') self.wrong(cmd.check_process_format,'u u~ > e+{L} vl [QED]') @test_aloha.set_global()