MadLoop: keep an external leg exactly at rest through improve_ps

.github/workflows/acceptancetest.yml

@@ -1075,6 +1075,22 @@ jobs:
             cd $GITHUB_WORKSPACE
             ./tests/test_manager.py test_ML_gauge_invariance_epem_ttx_all_modes -pA -t0 -l INFO
 
+  acceptancetest_improve_ps_at_rest:
+    # improve_ps must not move an external leg off rest: a polarised matrix
+    # element is evaluated in a frame where the selected leg has an exactly
+    # zero three-momentum, and HELAS reads its spin axis off that zero.
+    runs-on: ubuntu-24.04
+    if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true
+    steps:
+      - uses: actions/checkout@v4
+      - uses: ./.github/actions/checkout_mg5
+      - uses: ./.github/actions/restore_all
+
+      - name: test improve_ps keeps a leg at rest
+        run: |
+            cd $GITHUB_WORKSPACE
+            ./tests/test_manager.py test_improve_ps_keeps_a_leg_at_rest -pA -t0 -l INFO
+
   acceptancetest_check_madevent_me_vs_standalone:
     # Previously uncovered acceptance test in tests/acceptance_tests/test_cmd_madloop.py
     runs-on: ubuntu-24.04

Template/loop_material/StandAlone/Cards/MadLoopParams.dat

@@ -84,6 +84,9 @@
 !     You can add other evaluation method to check for the stability in DP and QP.
 !     Below you can chose if you want to use zero, one or two rotations of the PS point 
 !     in QP.
+!     WARNING:: do not use rotations for polarised matrix elements with an
+!     external particle at rest: the spin axis is tied to the frame axes, so
+!     a rotated evaluation legitimately differs and looks 'unstable'.
 #NRotations_DP
 !0
 ! Default :: 0

madgraph/iolibs/template_files/loop/improve_ps.inc

@@ -381,6 +381,11 @@ C stores the result in P and for the quadruple precision
 c version , it also modifies the global variables
 C PS and MP_DONE accordingly.
 
+C     NOTE:: an exact axis permutation, so exact zeros survive. But
+C     it is NOT a symmetry of a polarised amplitude with a leg at
+C     rest: HELAS pins that leg's spin axis to the frame z.
+C     NRotations_DP/QP default to 0; enabling them on a polarised
+C     run gives spurious instability flags.
 	  SUBROUTINE %(proc_prefix)sROTATE_PS(P_IN,P,ROTATION)
 	  IMPLICIT NONE
 C  
@@ -426,6 +431,11 @@ C rotation=2 => (xp=-z,yp=y,zp=x)
 	  END
 
 
+C     NOTE:: an exact axis permutation, so exact zeros survive. But
+C     it is NOT a symmetry of a polarised amplitude with a leg at
+C     rest: HELAS pins that leg's spin axis to the frame z.
+C     NRotations_DP/QP default to 0; enabling them on a polarised
+C     run gives spurious instability flags.
 	  SUBROUTINE %(proc_prefix)s%(mp_prefix)sROTATE_PS(P_IN,P,ROTATION)
 	  IMPLICIT NONE
 C  
@@ -666,7 +676,7 @@ C
 C
 C LOCAL VARIABLES 
 C
-      INTEGER I,J, P1, P2 
+      INTEGER I,J, P1, P2, IABSORB
       %(real_format)s NEWP(0:3,NEXTERNAL), PBUFF(0:3)
 	  %(real_format)s BUFF, BUFF2, XSCALE, APPROX_ZEROS(NAPPROXZEROS)
 	  %(real_format)s MASSES(NEXTERNAL)
@@ -706,13 +716,30 @@ C     First make sur that the space like momentum is exactly conserved
 		  PBUFF(J)=PBUFF(J)+NEWP(J,I)
 		ENDDO
       ENDDO
-	  DO I=NINITIAL+1,NEXTERNAL-1
-	    DO J=1,3
-		  PBUFF(J)=PBUFF(J)-NEWP(J,I)
-		ENDDO
+C     Absorb the residual in the final-state leg with the largest
+C     |p|, not in leg NEXTERNAL: that leg is never at rest, so an
+C     exactly vanishing three-momentum survives. HELAS reads the
+C     spin axis of a polarised leg off 'pp.eq.0', so a 1e-14
+C     residual there builds the longitudinal vector along noise.
+C     Largest |p| is also the smallest relative distortion.
+      IABSORB=NEXTERNAL
+	  BUFF=ZERO
+      DO I=NINITIAL+1,NEXTERNAL
+	    BUFF2=NEWP(1,I)**2+NEWP(2,I)**2+NEWP(3,I)**2
+		IF (BUFF2.GT.BUFF) THEN
+		  BUFF=BUFF2
+		  IABSORB=I
+		ENDIF
+	  ENDDO
+	  DO I=NINITIAL+1,NEXTERNAL
+	    IF (I.NE.IABSORB) THEN
+		  DO J=1,3
+		    PBUFF(J)=PBUFF(J)-NEWP(J,I)
+		  ENDDO
+		ENDIF
       ENDDO
       DO J=1,3
-        NEWP(J,NEXTERNAL)=PBUFF(J)
+        NEWP(J,IABSORB)=PBUFF(J)
       ENDDO
 
 C     Now find the 'x' rescaling factor

tests/acceptance_tests/test_cmd_madloop.py

@@ -186,6 +186,76 @@ def test_ML_check_brs_gd_gd(self):
             raise
         self.setup_logFile_for_logger('madgraph.check_cmd',restore=True)
 
+    def test_improve_ps_keeps_a_leg_at_rest(self):
+        """IMPROVE_PS_POINT_PRECISION must not destroy an exactly-zero
+        three-momentum.
+
+        A polarised matrix element is evaluated in a frame where the selected
+        leg is at rest, and boost_to_frame puts its three-momentum at exactly
+        zero on purpose: HELAS reads the spin quantisation axis of a massive
+        vector off 'pp.eq.rZero'. The PSMC branch of improve_ps used to
+        restore momentum conservation by dumping the residual into leg
+        NEXTERNAL unconditionally, so whenever the at-rest leg was the last
+        one its exact zero came back as a ~1e-14 vector and the longitudinal
+        polarisation vector was built along rounding noise. On g g > z{0}
+        z{0} [noborn=QCD] that split me_frame=[3] from me_frame=[4] by 26%%.
+
+        This drives the routine directly, so the check is exact rather than
+        statistical, and it covers both entry points: the double precision
+        wrapper and the quad routine reached through SET_MP_PS on the
+        stability-escalation path.
+        """
+        out_dir = pjoin(self.tmpdir, 'ML_improve_ps')
+        self.do('import model loop_sm')
+        self.do('generate g g > z z [sqrvirt=QCD]')
+        self.do('output standalone %s -f' % out_dir)
+
+        proc_dir = pjoin(out_dir, 'SubProcesses', 'P0_gg_zz')
+        prefix = open(pjoin(proc_dir, 'proc_prefix.txt')).read().strip()
+
+        # SET_MP_PS is the second entry point; it improves the raw momenta in
+        # quad, so it has to reach the same routine the driver calls.
+        loop_matrix = open(pjoin(proc_dir, 'loop_matrix.f')).read()
+        set_mp_ps = loop_matrix.split('SUBROUTINE %sSET_MP_PS(P)' % prefix)[1]
+        set_mp_ps = set_mp_ps.split('\n      END\n')[0]
+        self.assertIn('CALL %sMP_IMPROVE_PS_POINT_PRECISION' % prefix,
+                      set_mp_ps)
+
+        driver = open(pjoin(_pickle_path,
+                            'improve_ps_at_rest_driver.f')).read()
+        driver = driver.replace('__PFX__', prefix)
+        driver_path = pjoin(proc_dir, 'improve_ps_at_rest_driver.f')
+        open(driver_path, 'w').write(driver)
+
+        exe = pjoin(proc_dir, 'improve_ps_at_rest')
+        compile = subprocess.Popen(
+            ['gfortran', '-ffixed-line-length-132',
+             '-I%s' % proc_dir, '-I%s' % pjoin(out_dir, 'SubProcesses'),
+             driver_path, pjoin(proc_dir, 'improve_ps.f'), '-o', exe],
+            stdout=subprocess.PIPE, stderr=subprocess.STDOUT, cwd=proc_dir)
+        (out, _) = compile.communicate()
+        self.assertEqual(compile.returncode, 0,
+                         'could not build the improve_ps driver:\n%s'
+                         % out.decode())
+
+        run = subprocess.Popen([exe], stdout=subprocess.PIPE,
+                               stderr=subprocess.STDOUT, cwd=proc_dir)
+        (out, _) = run.communicate()
+        out = out.decode()
+        self.assertEqual(run.returncode, 0, out)
+
+        results = [l.split() for l in out.split('\n')
+                   if l.startswith('RESULT ')]
+        # 2 at-rest legs x 2 ImprovePSPoint settings x 2 entry points
+        self.assertEqual(len(results), 8, out)
+        for _, entry, leg, mode, px, py, pz in results:
+            for label, value in (('px', px), ('py', py), ('pz', pz)):
+                self.assertEqual(
+                    float(value), 0.0,
+                    'improve_ps moved leg %s off rest (%s, entry %s, '
+                    'ImprovePSPoint=%s): %s = %s'
+                    % (leg, label, entry, mode, label, value))
+
     def test_generate_output_semicolon_preserves_loop_process(self):
         """Regression test for combined generate/output in one command line."""
 

tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%improve_ps.f

@@ -403,6 +403,11 @@ REAL*16 FUNCTION MP_MDOT(P1,P2)
 C     version , it also modifies the global variables
 C     PS and MP_DONE accordingly.
 
+C     NOTE:: an exact axis permutation, so exact zeros survive. But
+C     it is NOT a symmetry of a polarised amplitude with a leg at
+C     rest: HELAS pins that leg's spin axis to the frame z.
+C     NRotations_DP/QP default to 0; enabling them on a polarised
+C     run gives spurious instability flags.
       SUBROUTINE ROTATE_PS(P_IN,P,ROTATION)
       IMPLICIT NONE
 C     
@@ -448,6 +453,11 @@ SUBROUTINE ROTATE_PS(P_IN,P,ROTATION)
       END
 
 
+C     NOTE:: an exact axis permutation, so exact zeros survive. But
+C     it is NOT a symmetry of a polarised amplitude with a leg at
+C     rest: HELAS pins that leg's spin axis to the frame z.
+C     NRotations_DP/QP default to 0; enabling them on a polarised
+C     run gives spurious instability flags.
       SUBROUTINE MP_ROTATE_PS(P_IN,P,ROTATION)
       IMPLICIT NONE
 C     
@@ -706,7 +716,7 @@ SUBROUTINE MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE,WARNED)
 C     
 C     LOCAL VARIABLES 
 C     
-      INTEGER I,J, P1, P2
+      INTEGER I,J, P1, P2, IABSORB
       REAL*16 NEWP(0:3,NEXTERNAL), PBUFF(0:3)
       REAL*16 BUFF, BUFF2, XSCALE, APPROX_ZEROS(NAPPROXZEROS)
       REAL*16 MASSES(NEXTERNAL)
@@ -750,13 +760,30 @@ SUBROUTINE MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE,WARNED)
           PBUFF(J)=PBUFF(J)+NEWP(J,I)
         ENDDO
       ENDDO
-      DO I=NINITIAL+1,NEXTERNAL-1
-        DO J=1,3
-          PBUFF(J)=PBUFF(J)-NEWP(J,I)
-        ENDDO
+C     Absorb the residual in the final-state leg with the largest
+C     |p|, not in leg NEXTERNAL: that leg is never at rest, so an
+C     exactly vanishing three-momentum survives. HELAS reads the
+C     spin axis of a polarised leg off 'pp.eq.0', so a 1e-14
+C     residual there builds the longitudinal vector along noise.
+C     Largest |p| is also the smallest relative distortion.
+      IABSORB=NEXTERNAL
+      BUFF=ZERO
+      DO I=NINITIAL+1,NEXTERNAL
+        BUFF2=NEWP(1,I)**2+NEWP(2,I)**2+NEWP(3,I)**2
+        IF (BUFF2.GT.BUFF) THEN
+          BUFF=BUFF2
+          IABSORB=I
+        ENDIF
+      ENDDO
+      DO I=NINITIAL+1,NEXTERNAL
+        IF (I.NE.IABSORB) THEN
+          DO J=1,3
+            PBUFF(J)=PBUFF(J)-NEWP(J,I)
+          ENDDO
+        ENDIF
       ENDDO
       DO J=1,3
-        NEWP(J,NEXTERNAL)=PBUFF(J)
+        NEWP(J,IABSORB)=PBUFF(J)
       ENDDO
 
 C     Now find the 'x' rescaling factor

tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%improve_ps.f

@@ -403,6 +403,11 @@ REAL*16 FUNCTION MP_MDOT(P1,P2)
 C     version , it also modifies the global variables
 C     PS and MP_DONE accordingly.
 
+C     NOTE:: an exact axis permutation, so exact zeros survive. But
+C     it is NOT a symmetry of a polarised amplitude with a leg at
+C     rest: HELAS pins that leg's spin axis to the frame z.
+C     NRotations_DP/QP default to 0; enabling them on a polarised
+C     run gives spurious instability flags.
       SUBROUTINE ROTATE_PS(P_IN,P,ROTATION)
       IMPLICIT NONE
 C     
@@ -448,6 +453,11 @@ SUBROUTINE ROTATE_PS(P_IN,P,ROTATION)
       END
 
 
+C     NOTE:: an exact axis permutation, so exact zeros survive. But
+C     it is NOT a symmetry of a polarised amplitude with a leg at
+C     rest: HELAS pins that leg's spin axis to the frame z.
+C     NRotations_DP/QP default to 0; enabling them on a polarised
+C     run gives spurious instability flags.
       SUBROUTINE MP_ROTATE_PS(P_IN,P,ROTATION)
       IMPLICIT NONE
 C     
@@ -706,7 +716,7 @@ SUBROUTINE MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE,WARNED)
 C     
 C     LOCAL VARIABLES 
 C     
-      INTEGER I,J, P1, P2
+      INTEGER I,J, P1, P2, IABSORB
       REAL*16 NEWP(0:3,NEXTERNAL), PBUFF(0:3)
       REAL*16 BUFF, BUFF2, XSCALE, APPROX_ZEROS(NAPPROXZEROS)
       REAL*16 MASSES(NEXTERNAL)
@@ -750,13 +760,30 @@ SUBROUTINE MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE,WARNED)
           PBUFF(J)=PBUFF(J)+NEWP(J,I)
         ENDDO
       ENDDO
-      DO I=NINITIAL+1,NEXTERNAL-1
-        DO J=1,3
-          PBUFF(J)=PBUFF(J)-NEWP(J,I)
-        ENDDO
+C     Absorb the residual in the final-state leg with the largest
+C     |p|, not in leg NEXTERNAL: that leg is never at rest, so an
+C     exactly vanishing three-momentum survives. HELAS reads the
+C     spin axis of a polarised leg off 'pp.eq.0', so a 1e-14
+C     residual there builds the longitudinal vector along noise.
+C     Largest |p| is also the smallest relative distortion.
+      IABSORB=NEXTERNAL
+      BUFF=ZERO
+      DO I=NINITIAL+1,NEXTERNAL
+        BUFF2=NEWP(1,I)**2+NEWP(2,I)**2+NEWP(3,I)**2
+        IF (BUFF2.GT.BUFF) THEN
+          BUFF=BUFF2
+          IABSORB=I
+        ENDIF
+      ENDDO
+      DO I=NINITIAL+1,NEXTERNAL
+        IF (I.NE.IABSORB) THEN
+          DO J=1,3
+            PBUFF(J)=PBUFF(J)-NEWP(J,I)
+          ENDDO
+        ENDIF
       ENDDO
       DO J=1,3
-        NEWP(J,NEXTERNAL)=PBUFF(J)
+        NEWP(J,IABSORB)=PBUFF(J)
       ENDDO
 
 C     Now find the 'x' rescaling factor

tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%improve_ps.f

@@ -404,6 +404,11 @@ REAL*16 FUNCTION MP_MDOT(P1,P2)
 C     version , it also modifies the global variables
 C     PS and MP_DONE accordingly.
 
+C     NOTE:: an exact axis permutation, so exact zeros survive. But
+C     it is NOT a symmetry of a polarised amplitude with a leg at
+C     rest: HELAS pins that leg's spin axis to the frame z.
+C     NRotations_DP/QP default to 0; enabling them on a polarised
+C     run gives spurious instability flags.
       SUBROUTINE ROTATE_PS(P_IN,P,ROTATION)
       IMPLICIT NONE
 C     
@@ -449,6 +454,11 @@ SUBROUTINE ROTATE_PS(P_IN,P,ROTATION)
       END
 
 
+C     NOTE:: an exact axis permutation, so exact zeros survive. But
+C     it is NOT a symmetry of a polarised amplitude with a leg at
+C     rest: HELAS pins that leg's spin axis to the frame z.
+C     NRotations_DP/QP default to 0; enabling them on a polarised
+C     run gives spurious instability flags.
       SUBROUTINE MP_ROTATE_PS(P_IN,P,ROTATION)
       IMPLICIT NONE
 C     
@@ -708,7 +718,7 @@ SUBROUTINE MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE,WARNED)
 C     
 C     LOCAL VARIABLES 
 C     
-      INTEGER I,J, P1, P2
+      INTEGER I,J, P1, P2, IABSORB
       REAL*16 NEWP(0:3,NEXTERNAL), PBUFF(0:3)
       REAL*16 BUFF, BUFF2, XSCALE, APPROX_ZEROS(NAPPROXZEROS)
       REAL*16 MASSES(NEXTERNAL)
@@ -753,13 +763,30 @@ SUBROUTINE MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE,WARNED)
           PBUFF(J)=PBUFF(J)+NEWP(J,I)
         ENDDO
       ENDDO
-      DO I=NINITIAL+1,NEXTERNAL-1
-        DO J=1,3
-          PBUFF(J)=PBUFF(J)-NEWP(J,I)
-        ENDDO
+C     Absorb the residual in the final-state leg with the largest
+C     |p|, not in leg NEXTERNAL: that leg is never at rest, so an
+C     exactly vanishing three-momentum survives. HELAS reads the
+C     spin axis of a polarised leg off 'pp.eq.0', so a 1e-14
+C     residual there builds the longitudinal vector along noise.
+C     Largest |p| is also the smallest relative distortion.
+      IABSORB=NEXTERNAL
+      BUFF=ZERO
+      DO I=NINITIAL+1,NEXTERNAL
+        BUFF2=NEWP(1,I)**2+NEWP(2,I)**2+NEWP(3,I)**2
+        IF (BUFF2.GT.BUFF) THEN
+          BUFF=BUFF2
+          IABSORB=I
+        ENDIF
+      ENDDO
+      DO I=NINITIAL+1,NEXTERNAL
+        IF (I.NE.IABSORB) THEN
+          DO J=1,3
+            PBUFF(J)=PBUFF(J)-NEWP(J,I)
+          ENDDO
+        ENDIF
       ENDDO
       DO J=1,3
-        NEWP(J,NEXTERNAL)=PBUFF(J)
+        NEWP(J,IABSORB)=PBUFF(J)
       ENDDO
 
 C     Now find the 'x' rescaling factor