NLO+PS polarisation: put the MC-counterterm azimuth in the me_frame

.github/workflows/acceptancetest.yml

@@ -1410,3 +1410,24 @@ jobs:
         run: |
             cd $GITHUB_WORKSPACE
             ./tests/test_manager.py test_polarised_nlo_me_frame -pA -t0 -l INFO
+
+  acceptancetest_106:
+    # 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 NLO+PS (MC@NLO matching) in a chosen rest frame. Asserts that
+      # the total cross-section is the fixed-order one while the absolute one
+      # stays low, which is what the azimuthal phase of the MC counterterm
+      # controls.
+      - name: test one of the test test_polarised_nlo_ps_me_frame
+        run: |
+            cd $GITHUB_WORKSPACE
+            ./tests/test_manager.py test_polarised_nlo_ps_me_frame -pA -t0 -l INFO

Template/NLO/SubProcesses/boost_to_frame.f

@@ -94,6 +94,111 @@ subroutine azifact_me_frame(p_born_in, imother, azifact, boosted)
       end
 
 
+      subroutine me_frame_born_boost(p_born_in, pboost, boosted)
+c**************************************************************************
+c     The boost sborn_frame() would apply to these Born momenta, exposed so
+c     the MC counterterms can put their azimuthal ingredients in the same
+c     frame. boosted=.false. means no frame was requested.
+c**************************************************************************
+      implicit none
+      include 'nexternal.inc'
+      double precision p_born_in(0:3,nexternal-1), pboost(0:3)
+      logical boosted
+      integer ids(nexternal-1)
+      logical trivial
+
+      call get_frame_mask_born(ids)
+      call get_me_frame_boost(p_born_in, nexternal-1, ids, pboost,
+     &                        trivial)
+      boosted=.not.trivial
+
+      return
+      end
+
+
+      subroutine azifact_from_spinors(pi, pj, azifact)
+c**************************************************************************
+c     <ij>/[ij] from two massless momenta -- the expression the legacy
+c     branches build inline. Those keep their own copies on purpose:
+c     rewriting them would break the byte-identity of the no-frame path.
+c**************************************************************************
+      USE ALOHA_OBJECT
+      implicit none
+      double precision pi(0:3), pj(0:3)
+      double complex azifact
+      double precision zero
+      parameter (zero=0d0)
+      TYPE(ALOHA) W1,W2,W3,W4
+      double complex Wij_angle,Wij_recta
+      integer i
+
+      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
+         Wij_angle = Wij_angle + W1%W(i)*W2%W(i)
+         Wij_recta = Wij_recta + W3%W(i)*W4%W(i)
+      enddo
+      azifact=Wij_angle/Wij_recta
+
+      return
+      end
+
+
+      subroutine azifact_mc_frame(p_born_in, imother, p_i, p_j,
+     &                            y_ij_fks, pboost, azifact,
+     &                            cphi_mother, sphi_mother)
+c**************************************************************************
+c     Azimuthal ingredients of the MC counterterms, in the me_frame, scaled
+c     so that the caller's existing multiply comes out right:
+c         ISR   -(cphi_mother + i sphi_mother)**2 * conjg(azifact)
+c         FSR   -(cphi_mother - i sphi_mother)**2 *       azifact
+c     both = +-exp(-+2 i psi), psi the emission azimuth about the mother.
+c     The boost (pboost, from me_frame_born_boost) is the Born's own, since
+c     the phase multiplies borntilde. azifact is <ij>/[ij] only in the
+c     generic branch; the collinear one pre-cancels the caller's mother
+c     factor instead, so do not assume it if you change the caller.
+c**************************************************************************
+      implicit none
+      include 'nexternal.inc'
+      double precision p_born_in(0:3,nexternal-1)
+      integer imother
+      double precision p_i(0:3), p_j(0:3), y_ij_fks, pboost(0:3)
+      double complex azifact
+      double precision cphi_mother, sphi_mother
+      double precision pm(0:3), kp(0:3), pib(0:3), pjb(0:3)
+      double precision xij_kperp(0:3)
+      common/cxij_kperp/xij_kperp
+
+      call boostx(p_born_in(0,imother), pboost, pm)
+      call getaziangles(pm, cphi_mother, sphi_mother)
+
+      if (1d0-y_ij_fks.le.0d0) then
+c        Exactly collinear: <ij>/[ij] is 0/0 in every frame, so rebuild it
+c        from the stored emission direction, using
+c            lim_{kt->0} <ij>/[ij] = -exp(2 i psi) * exp(2 i phi_m),
+c        the identity that makes this branch agree with the generic one;
+c        azifact_from_kperp returns only the first factor.
+c
+c        Taken at kt exactly zero and nowhere else, on purpose: the spinors
+c        are exact at any kt>0, while this reconstruction carries an
+c        O((kt/E)/betaT) error, betaT the transverse velocity of the boost.
+         call boostx(xij_kperp, pboost, kp)
+         call azifact_from_kperp(pm, kp, azifact)
+         azifact=azifact*dcmplx(cphi_mother,sphi_mother)**2
+      else
+         call boostx(p_i, pboost, pib)
+         call boostx(p_j, pboost, pjb)
+         call azifact_from_spinors(pib, pjb, azifact)
+      endif
+
+      return
+      end
+
+
       subroutine azifact_from_kperp(pmother, kperp, azifact)
 c**************************************************************************
 c     Rebuild the collinear limit of <ij>/[ij] from the emission direction.

Template/NLO/SubProcesses/montecarlocounter.f

@@ -2867,6 +2867,9 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
       double complex ximag
       parameter (ximag=(0.d0,1.d0))
 
+      double precision pboost(0:3)
+      logical me_boosted
+
       double precision xi_i_fks_ev,y_ij_fks_ev,t
       double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
       common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
@@ -2908,9 +2911,13 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
 
 c
 c BORN/BORNTILDE
+      call me_frame_born_boost(p_born,pboost,me_boosted)
 C check if momenta have to be rotated
+c The R_y(pi) flip only put a j_fks=2 mother back on +z, meaningless once the
+c mother is off the beam axis; dropped under a frame, as in sborncol_isr. Do
+c not reinstate: it would flip the polarisation axis of a leg at rest.
       if ((ileg.eq.1.or.ileg.eq.2) .and.
-     &    (j_fks.eq.2 .and. nexternal-1.ne.3)) then
+     &    (j_fks.eq.2 .and. nexternal-1.ne.3) .and. .not.me_boosted) then
 c Rotation according to innerpin.m. Use rotate_invar() if a more 
 c general rotation is needed.
 c Exclude 2->1 (at the Born level) processes: matrix elements are
@@ -2969,6 +2976,13 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
 c BORN TILDE
       if(ileg.eq.1.or.ileg.eq.2)then
 c Insert <ij>/[ij] which is not included by sborn()
+       if (me_boosted) then
+c borntilde was computed in the me_frame, so the phase multiplying it has to
+c be as well: emission spinors and mother azimuth both from boosted momenta.
+          call azifact_mc_frame(p_born,j_fks,p_i_fks_ev,p(0,j_fks),
+     &                          y_ij_fks,pboost,azifact,cphi_mother,
+     &                          sphi_mother)
+       else
          if (1d0-y_ij_fks.lt.vtiny)then
             azifact=xij_aor
          else
@@ -3004,6 +3018,7 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
             cphi_mother=1.d0
             sphi_mother=0.d0
          endif
+       endif
          do iord=1, nsplitorders
            borntilde(iord) = -(cphi_mother+ximag*sphi_mother)**2 *
      #                borntilde(iord) * dconjg(azifact)
@@ -3024,6 +3039,13 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
             enddo
          elseif((m_type.eq.8.or.m_type.eq.1).and.ch_m.eq.0d0)then
 c Insert <ij>/[ij] which is not included by sborn()
+            imother_fks=min(i_fks,j_fks)
+            if (me_boosted) then
+c Same as the ISR case: both ingredients in the frame borntilde lives in.
+               call azifact_mc_frame(p_born,imother_fks,p_i_fks_ev,
+     &                               p(0,j_fks),y_ij_fks,pboost,azifact,
+     &                               cphi_mother,sphi_mother)
+            else
             if(1.d0-y_ij_fks.lt.vtiny)then
                azifact=xij_aor
             else
@@ -3044,9 +3066,9 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
                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)
+            endif
             do iord=1, nsplitorders
                borntilde(iord) = -(cphi_mother-ximag*sphi_mother)**2 *
      #                  borntilde(iord) * azifact

Template/NLO/SubProcesses/montecarlocounter_alt.f

@@ -1329,6 +1329,9 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
       double complex ximag
       parameter (ximag=(0.d0,1.d0))
 
+      double precision pboost(0:3)
+      logical me_boosted
+
       double precision xi_i_fks_ev,y_ij_fks_ev,t
       double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
       common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
@@ -1372,6 +1375,7 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
 
 c
 C BORN
+      call me_frame_born_boost(p_born,pboost,me_boosted)
       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
@@ -1391,8 +1395,10 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
 
 c BORN TILDE
 C check if momenta have to be rotated
+c The R_y(pi) flip only ever served to put a j_fks=2 mother back on +z; in
+c the me_frame the mother is off the beam axis anyway. See montecarlocounter.f.
       if (ileg.eq.1.or.ileg.eq.2) then
-        if (j_fks.eq.2 .and. nexternal-1.ne.3) then
+        if (j_fks.eq.2 .and. nexternal-1.ne.3 .and. .not.me_boosted) then
 c Rotation according to innerpin.m. Use rotate_invar() if a more 
 c general rotation is needed.
 c Exclude 2->1 (at the Born level) processes: matrix elements are
@@ -1447,6 +1453,13 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
         enddo
 
 c Insert <ij>/[ij] which is not included by sborn()
+        if (me_boosted) then
+c borntilde was computed in the me_frame, so the phase multiplying it has to
+c be as well: emission spinors and mother azimuth both from boosted momenta.
+           call azifact_mc_frame(p_born,j_fks,p_i_fks_ev,p(0,j_fks),
+     &                           y_ij_fks,pboost,azifact,cphi_mother,
+     &                           sphi_mother)
+        else
         if (1d0-y_ij_fks.lt.vtiny)then
            azifact=xij_aor
         else
@@ -1482,6 +1495,7 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
            cphi_mother=1.d0
            sphi_mother=0.d0
         endif
+        endif
         do iord=1, nsplitorders
            borntilde(iord) = -(cphi_mother+ximag*sphi_mother)**2 *
      #                borntilde(iord) * dconjg(azifact)
@@ -1502,6 +1516,13 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
             enddo
          elseif((m_type.eq.8.or.m_type.eq.1).and.ch_m.eq.0d0)then
 c Insert <ij>/[ij] which is not included by sborn()
+            imother_fks=min(i_fks,j_fks)
+            if (me_boosted) then
+c Same as the ISR case: both ingredients in the frame borntilde lives in.
+               call azifact_mc_frame(p_born,imother_fks,p_i_fks_ev,
+     &                               p(0,j_fks),y_ij_fks,pboost,azifact,
+     &                               cphi_mother,sphi_mother)
+            else
             if(1.d0-y_ij_fks.lt.vtiny)then
                azifact=xij_aor
             else
@@ -1522,9 +1543,9 @@ subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
                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)
+            endif
             do iord=1, nsplitorders
                borntilde(iord) = -(cphi_mother-ximag*sphi_mother)**2 *
      #                  borntilde(iord) * azifact