      SUBROUTINE ML5_0_SLOOPMATRIXHEL(P,HEL,ANS)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)

      INCLUDE 'nsquaredSO.inc'

C     
C     ARGUMENTS 
C     
      REAL*8 P(0:3,NEXTERNAL)
      REAL*8 ANS(0:3,0:NSQUAREDSO)
      INTEGER HEL, USERHEL
      COMMON/ML5_0_USERCHOICE/USERHEL
C     ----------
C     BEGIN CODE
C     ----------
      USERHEL=HEL
      CALL ML5_0_SLOOPMATRIX(P,ANS)
      END

      LOGICAL FUNCTION ML5_0_IS_HEL_SELECTED(HELID)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NCOMB
      PARAMETER (NCOMB=32)
C     
C     ARGUMENTS
C     
      INTEGER HELID
C     
C     LOCALS
C     
      INTEGER I,J
      LOGICAL FOUNDIT
C     
C     GLOBALS
C     
      INTEGER HELC(NEXTERNAL,NCOMB)
      COMMON/ML5_0_HELCONFIGS/HELC

      INTEGER POLARIZATIONS(0:NEXTERNAL,0:5)
      COMMON/ML5_0_BEAM_POL/POLARIZATIONS
C     ----------
C     BEGIN CODE
C     ----------

      ML5_0_IS_HEL_SELECTED = .TRUE.
      IF (POLARIZATIONS(0,0).EQ.-1) THEN
        RETURN
      ENDIF

      DO I=1,NEXTERNAL
        IF (POLARIZATIONS(I,0).EQ.-1) THEN
          CYCLE
        ENDIF
        FOUNDIT = .FALSE.
        DO J=1,POLARIZATIONS(I,0)
          IF (HELC(I,HELID).EQ.POLARIZATIONS(I,J)) THEN
            FOUNDIT = .TRUE.
            EXIT
          ENDIF
        ENDDO
        IF(.NOT.FOUNDIT) THEN
          ML5_0_IS_HEL_SELECTED = .FALSE.
          RETURN
        ENDIF
      ENDDO
      RETURN

      END

      LOGICAL FUNCTION ML5_0_ISZERO(TOTEST, REFERENCE_VALUE, AMPLN)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NLOOPAMPS
      PARAMETER (NLOOPAMPS=39)
C     
C     ARGUMENTS 
C     
      REAL*8 TOTEST, REFERENCE_VALUE
      INTEGER AMPLN
C     
C     GLOBAL 
C     
      INCLUDE 'MadLoopParams.inc'

      COMPLEX*16 AMPL(3,NLOOPAMPS)
      LOGICAL S(NLOOPAMPS)
      COMMON/ML5_0_AMPL/AMPL,S
C     ----------
C     BEGIN CODE
C     ----------
      IF(ABS(REFERENCE_VALUE).EQ.0.0D0) THEN
        ML5_0_ISZERO=.FALSE.
        WRITE(*,*) '##E02 ERRROR Reference value for comparison is'
     $   //' zero.'
        STOP
      ELSE
        ML5_0_ISZERO=((ABS(TOTEST)/ABS(REFERENCE_VALUE)).LT.ZEROTHRES)
      ENDIF
      IF(AMPLN.NE.-1) THEN
        IF((.NOT.ML5_0_ISZERO).AND.(.NOT.S(AMPLN))) THEN
          WRITE(*,*) '##W01 WARNING Contribution ',AMPLN,' is detected'
     $     //' as contributing with CR=',(ABS(TOTEST)
     $     /ABS(REFERENCE_VALUE)),' but is unstable.'
        ENDIF
      ENDIF

      END

      SUBROUTINE ML5_0_SLOOPMATRIX(P_USER,ANSRETURNED)
C     
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C     
C     Returns amplitude squared summed/avg over colors
C     and helicities for the point in phase space P(0:3,NEXTERNAL)
C     and external lines W(0:6,NEXTERNAL)
C     
C     Process: d u~ > m- vm~ g QCD<=1 QED<=2 [ virt = QCD ]
C     
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      CHARACTER*512 PARAMFNAME,HELCONFIGFNAME,LOOPFILTERFNAME
      CHARACTER*512 COLORNUMFNAME,COLORDENOMFNAME, HELFILTERFNAME
      CHARACTER*512 PROC_PREFIX
      PARAMETER ( PARAMFNAME='MadLoopParams.dat')
      PARAMETER ( HELCONFIGFNAME='HelConfigs.dat')
      PARAMETER ( LOOPFILTERFNAME='LoopFilter.dat')
      PARAMETER ( HELFILTERFNAME='HelFilter.dat')
      PARAMETER ( COLORNUMFNAME='ColorNumFactors.dat')
      PARAMETER ( COLORDENOMFNAME='ColorDenomFactors.dat')
      PARAMETER ( PROC_PREFIX='ML5_0_')

      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=2)
      INTEGER    NLOOPAMPS, NCTAMPS
      PARAMETER (NLOOPAMPS=39, NCTAMPS=28)
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER NINITIAL
      PARAMETER (NINITIAL=2)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=10)
      INTEGER    NCOMB
      PARAMETER (NCOMB=32)
      REAL*8     ZERO
      PARAMETER (ZERO=0D0)
      REAL*16     MP__ZERO
      PARAMETER (MP__ZERO=0E0_16)
      COMPLEX*16 IMAG1
      PARAMETER (IMAG1=(0D0,1D0))
C     This parameter is designed for the check timing command of MG5
      LOGICAL SKIPLOOPEVAL
      PARAMETER (SKIPLOOPEVAL=.FALSE.)
      LOGICAL BOOTANDSTOP
      PARAMETER (BOOTANDSTOP=.FALSE.)
      INCLUDE 'nsquaredSO.inc'
      INTEGER NSQUAREDSOP1
      PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1)
      INTEGER MAXSTABILITYLENGTH
      DATA MAXSTABILITYLENGTH/20/
      COMMON/ML5_0_STABILITY_TESTS/MAXSTABILITYLENGTH
C     
C     ARGUMENTS 
C     
      REAL*8 P_USER(0:3,NEXTERNAL)
      REAL*8 ANSRETURNED(0:3,0:NSQUAREDSO)
C     
C     LOCAL VARIABLES 
C     
      REAL*8 ANS(0:3)
      INTEGER I,J,K,H

      CHARACTER*512 PARAMFN,HELCONFIGFN,LOOPFILTERFN,COLORNUMFN
     $ ,COLORDENOMFN,HELFILTERFN
      CHARACTER*512 TMP
      SAVE PARAMFN
      SAVE HELCONFIGFN
      SAVE LOOPFILTERFN
      SAVE COLORNUMFN
      SAVE COLORDENOMFN
      SAVE HELFILTERFN

      INTEGER HELPICKED_BU, CTMODEINIT_BU
      REAL*8 MLSTABTHRES_BU
C     P is the actual PS POINT used for the computation, and can be
C      rotated for the stability test purposes.
      REAL*8 P(0:3,NEXTERNAL)
C     DP_RES STORES THE DOUBLE PRECISION RESULT OBTAINED FROM
C      DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY.
C     THE STAB_STAGE COUNTER I CORRESPONDANCE GOES AS FOLLOWS
C     I=1 -> ORIGINAL PS, CTMODE=1
C     I=2 -> ORIGINAL PS, CTMODE=2, (ONLY WITH CTMODERUN=-1)
C     I=3 -> PS WITH ROTATION 1, CTMODE=1, (ONLY WITH CTMODERUN=-2)
C     I=4 -> PS WITH ROTATION 2, CTMODE=1, (ONLY WITH CTMODERUN=-3)
C     I=5 -> POSSIBLY MORE EVALUATION METHODS IN THE FUTURE, MAX IS
C      MAXSTABILITYLENGTH
C     IF UNSTABLE IT GOES TO THE SAME PATTERN BUT STAB_INDEX IS THEN
C      I+20.
      LOGICAL EVAL_DONE(MAXSTABILITYLENGTH)
      LOGICAL DOING_QP_EVALS
      INTEGER STAB_INDEX,BASIC_CT_MODE
      INTEGER N_DP_EVAL, N_QP_EVAL
      DATA N_DP_EVAL/1/
      DATA N_QP_EVAL/1/
C     This is used for loop-induced where the reference scale for
C      comparisons is infered from
C     the previous points
      REAL*8 NEXTREF
      DATA NEXTREF/ZERO/
      INTEGER NPSPOINTS
      DATA NPSPOINTS/0/
      LOGICAL FOUND_VALID_REDUCTION_METHOD
      DATA FOUND_VALID_REDUCTION_METHOD/.FALSE./

      REAL*8 ACC
      REAL*8 DP_RES(3,MAXSTABILITYLENGTH)
C     QP_RES STORES THE QUADRUPLE PRECISION RESULT OBTAINED FROM
C      DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY.
      REAL*8 QP_RES(3,MAXSTABILITYLENGTH)
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
      INTEGER NATTEMPTS
      DATA NATTEMPTS/0/
      DATA IC/NEXTERNAL*1/
      REAL*8 BUFFR(3),TEMP(3),TEMP1,TEMP2
      COMPLEX*16 CFTOT
      LOGICAL FOUNDHELFILTER,FOUNDLOOPFILTER
      DATA FOUNDHELFILTER/.TRUE./
      DATA FOUNDLOOPFILTER/.TRUE./
      INTEGER IDEN
      DATA IDEN/36/
      INTEGER HELAVGFACTOR
      DATA HELAVGFACTOR/4/
C     For a 1>N process, them BEAMTWO_HELAVGFACTOR would be set to 1.
      INTEGER BEAMS_HELAVGFACTOR(2)
      DATA (BEAMS_HELAVGFACTOR(I),I=1,2)/2,2/
      LOGICAL DONEHELDOUBLECHECK
      DATA DONEHELDOUBLECHECK/.FALSE./
      INTEGER NEPS
      DATA NEPS/0/
C     Below are variables to bypass the checkphase and insure
C      stability check to take place
      LOGICAL OLD_CHECKPHASE, OLD_HELDOUBLECHECKED
      LOGICAL OLD_GOODHEL(NCOMB)
      LOGICAL OLD_GOODAMP(NLOOPAMPS,NCOMB)

      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
      COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
C     
C     FUNCTIONS
C     
      LOGICAL ML5_0_ISZERO
      LOGICAL ML5_0_IS_HEL_SELECTED
C     
C     GLOBAL VARIABLES
C     
      INCLUDE 'process_info.inc'
      INCLUDE 'coupl.inc'
      INCLUDE 'mp_coupl.inc'
      INCLUDE 'MadLoopParams.inc'

      INTEGER NTRY
      DATA NTRY/0/
      LOGICAL CHECKPHASE
      DATA CHECKPHASE/.TRUE./
      LOGICAL HELDOUBLECHECKED
      DATA HELDOUBLECHECKED/.FALSE./
      REAL*8 REF
      DATA REF/0.0D0/
      COMMON/ML5_0_INIT/NTRY,CHECKPHASE,HELDOUBLECHECKED,REF

C     THE LOGICAL BELOWS ARE JUST TO KEEP TRACK OF WHETHER THE MP_PS
C      HAS BEEN SET YET OR NOT AND WHETER THE MP EXTERNAL WFS HAVE
C      BEEN COMPUTED YET.
      LOGICAL MP_DONE
      DATA MP_DONE/.FALSE./
      COMMON/ML5_0_MP_DONE/MP_DONE
      LOGICAL MP_PS_SET
      DATA MP_PS_SET/.FALSE./
      COMMON/ML5_0_MP_PS_SET/MP_PS_SET

C     PS CAN POSSIBILY BE PASSED THROUGH IMPROVE_PS BUT IS NOT
C      MODIFIED FOR THE PURPOSE OF THE STABILITY TEST	  
C     EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT
C      USED ANYWHERE ELSE
      REAL*8 PS(0:3,NEXTERNAL)
      COMMON/ML5_0_PSPOINT/PS
C     AGAIN BELOW, MP_PS IS THE FIXED (POSSIBLY IMPROVED) MP PS POINT
C      AND MP_P IS THE ONE WHICH CAN BE MODIFIED (I.E. ROTATED ETC.)
C      FOR STABILITY PURPOSE
C     EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT
C      USED ANYWHERE ELSE THAN HERE AND SET_MP_PS()
      REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL)
      COMMON/ML5_0_MP_PSPOINT/MP_PS,MP_P

      REAL*8 LSCALE
      INTEGER CTMODE
      COMMON/ML5_0_CT/LSCALE,CTMODE

      LOGICAL GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPAMPS,NCOMB)
      COMMON/ML5_0_FILTERS/GOODAMP,GOODHEL

      INTEGER HELPICKED
      DATA HELPICKED/-1/
      COMMON/ML5_0_HELCHOICE/HELPICKED
      INTEGER USERHEL
      DATA USERHEL/-1/
      COMMON/ML5_0_USERCHOICE/USERHEL

      COMPLEX*16 AMP(NBORNAMPS,NCOMB)
      COMMON/ML5_0_AMPS/AMP
      COMPLEX*16 W(20,NWAVEFUNCS,NCOMB)
      INTEGER VALIDH
      COMMON/ML5_0_WFCTS/W
      COMMON/ML5_0_VALIDH/VALIDH

      COMPLEX*16 AMPL(3,NLOOPAMPS)
      LOGICAL S(NLOOPAMPS)
      COMMON/ML5_0_AMPL/AMPL,S

      INTEGER CF_D(NLOOPAMPS,NBORNAMPS)
      INTEGER CF_N(NLOOPAMPS,NBORNAMPS)
      COMMON/ML5_0_CF/CF_D,CF_N

      INTEGER HELC(NEXTERNAL,NCOMB)
      COMMON/ML5_0_HELCONFIGS/HELC

      REAL*8 PREC,USER_STAB_PREC
      DATA USER_STAB_PREC/-1.0D0/
      COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC

C     Return codes H,T,U correspond to the hundreds, tens and units
C     building returncode, i.e.
C     RETURNCODE=100*RET_CODE_H+10*RET_CODE_T+RET_CODE_U

      INTEGER RET_CODE_H,RET_CODE_T,RET_CODE_U
      REAL*8 ACCURACY(0:NSQUAREDSO)
      DATA (ACCURACY(I),I=0,NSQUAREDSO)/NSQUAREDSOP1*1.0D0/
      DATA RET_CODE_H,RET_CODE_T,RET_CODE_U/1,1,0/
      COMMON/ML5_0_ACC/ACCURACY,RET_CODE_H,RET_CODE_T,RET_CODE_U

C     Allows to forbid the zero helicity double check, no matter the
C      value in MadLoopParams.dat
C     This can be accessed with the SET_FORBID_HEL_DOUBLECHECK
C      subroutine of MadLoopCommons.dat
      LOGICAL FORBID_HEL_DOUBLECHECK
      COMMON/FORBID_HEL_DOUBLECHECK/FORBID_HEL_DOUBLECHECK

      LOGICAL MP_DONE_ONCE
      DATA MP_DONE_ONCE/.FALSE./
      COMMON/ML5_0_MP_DONE_ONCE/MP_DONE_ONCE

      CHARACTER(512) MLPATH
      COMMON/MLPATH/MLPATH

      LOGICAL ML_INIT
      COMMON/ML_INIT/ML_INIT

C     This variable controls the *local* initialization of this
C      particular SubProcess.
C     For example, the reading of the filters must be done
C      independently by each SubProcess.
      LOGICAL LOCAL_ML_INIT
      DATA LOCAL_ML_INIT/.TRUE./

C     Variables related to turning off the Lorentz rotation test when
C      spin-2 particles are external
      LOGICAL WARNED_LORENTZ_STAB_TEST_OFF
      DATA WARNED_LORENTZ_STAB_TEST_OFF/.FALSE./
      INTEGER NROTATIONS_DP_BU,NROTATIONS_QP_BU

C     This array specify potential special requirements on the
C      helicities to
C     consider. POLARIZATIONS(0,0) is -1 if there is not such
C      requirement.
      INTEGER POLARIZATIONS(0:NEXTERNAL,0:5)
      COMMON/ML5_0_BEAM_POL/POLARIZATIONS

C     ----------
C     BEGIN CODE
C     ----------

      IF(ML_INIT) THEN
        CALL PRINT_MADLOOP_BANNER()
        TMP = 'auto'
        CALL SETMADLOOPPATH(TMP)
        CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN)
        CALL MADLOOPPARAMREADER(PARAMFN,.TRUE.)
        IF (FORBID_HEL_DOUBLECHECK) THEN
          DOUBLECHECKHELICITYFILTER = .FALSE.
        ENDIF
        ML_INIT = .FALSE.
C       For now only CutTools is interfaced in the default mode.
C        Samurai could follow.
        DO I=1,SIZE(MLREDUCTIONLIB)
          IF (MLREDUCTIONLIB(I).EQ.1) THEN
            FOUND_VALID_REDUCTION_METHOD = .TRUE.
          ENDIF
        ENDDO
        IF (.NOT.FOUND_VALID_REDUCTION_METHOD) THEN
          WRITE(*,*) 'ERROR:: For now, only CutTools is interfaced to'
     $     //' MadLoop in the non-optimized output.'
          WRITE(*,*) 'ERROR:: Make sure to include 1 in the parameter'
     $     //' MLReductionLib of the card MadLoopParams.dat'
          STOP 1
        ENDIF
      ENDIF
      IF (LOCAL_ML_INIT) THEN
C       Setup the file paths
        CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN)
        CALL JOINPATH(MLPATH,PROC_PREFIX,TMP)
        CALL JOINPATH(TMP,HELCONFIGFNAME,HELCONFIGFN)
        CALL JOINPATH(TMP,LOOPFILTERFNAME,LOOPFILTERFN)
        CALL JOINPATH(TMP,COLORNUMFNAME,COLORNUMFN)
        CALL JOINPATH(TMP,COLORDENOMFNAME,COLORDENOMFN)
        CALL JOINPATH(TMP,HELFILTERFNAME,HELFILTERFN)

C       Make sure that the loop filter is disabled when there is
C        spin-2 particles for 2>1 or 1>2 processes
        IF(MAX_SPIN_EXTERNAL_PARTICLE.GT.3.AND.(NEXTERNAL.LE.3.AND.HELI
     $CITYFILTERLEVEL.NE.0)) THEN
          WRITE(*,*) '##INFO: Helicity filter deactivated for 2>1'
     $     //' processes involving spin 2 particles.'
          HELICITYFILTERLEVEL = 0
C         We write a dummy filter for structural reasons here
          OPEN(1, FILE=HELFILTERFN, ERR=6116, STATUS='NEW'
     $     ,ACTION='WRITE')
          DO I=1,NCOMB
            WRITE(1,*) 'T'
          ENDDO
 6116     CONTINUE
          CLOSE(1)
        ENDIF

        OPEN(1, FILE=COLORNUMFN, ERR=104, STATUS='OLD',          
     $    ACTION='READ')
        DO I=1,NLOOPAMPS
          READ(1,*,END=105) (CF_N(I,J),J=1,NBORNAMPS)
        ENDDO
        GOTO 105
 104    CONTINUE
        STOP 'Color factors could not be initialized from file'
     $   //' ML5_0_ColorNumFactors.dat. File not found'
 105    CONTINUE
        CLOSE(1)
        OPEN(1, FILE=COLORDENOMFN, ERR=106, STATUS='OLD',          
     $    ACTION='READ')
        DO I=1,NLOOPAMPS
          READ(1,*,END=107) (CF_D(I,J),J=1,NBORNAMPS)
        ENDDO
        GOTO 107
 106    CONTINUE
        STOP 'Color factors could not be initialized from file'
     $   //' ML5_0_ColorDenomFactors.dat. File not found'
 107    CONTINUE
        CLOSE(1)
        OPEN(1, FILE=HELCONFIGFN, ERR=108, STATUS='OLD',              
     $       ACTION='READ')
        DO H=1,NCOMB
          READ(1,*,END=109) (HELC(I,H),I=1,NEXTERNAL)
        ENDDO
        GOTO 109
 108    CONTINUE
        STOP 'Color helictiy configurations could not be initialized'
     $   //' from file ML5_0_HelConfigs.dat. File not found'
 109    CONTINUE
        CLOSE(1)
        IF(BOOTANDSTOP) THEN
          WRITE(*,*) '##Stopped by user request.'
          STOP
        ENDIF
        LOCAL_ML_INIT = .FALSE.
      ENDIF

C     Make sure that lorentz rotation tests are not used if there is
C      external loop wavefunction of spin 2 and that one specific
C      helicity is asked
      NROTATIONS_DP_BU = NROTATIONS_DP
      NROTATIONS_QP_BU = NROTATIONS_QP
      IF(MAX_SPIN_EXTERNAL_PARTICLE.GT.3.AND.USERHEL.NE.-1) THEN
        IF(.NOT.WARNED_LORENTZ_STAB_TEST_OFF) THEN
          WRITE(*,*) '##WARNING: Evaluation of a specific helicity was'
     $     //' asked for this PS point, and there is a spin-2 (or'
     $     //' higher) particle in the external states.'
          WRITE(*,*) '##WARNING: As a result, MadLoop disabled the'
     $     //' Lorentz rotation test for this phase-space point only.'
          WRITE(*,*) '##WARNING: Further warning of that type'
     $     //' suppressed.'
          WARNED_LORENTZ_STAB_TEST_OFF = .FALSE.
        ENDIF
        NROTATIONS_QP=0
        NROTATIONS_DP=0
      ENDIF

      IF(NTRY.EQ.0) THEN
        CALL ML5_0_SET_N_EVALS(N_DP_EVAL,N_QP_EVAL)
        HELDOUBLECHECKED=(.NOT.DOUBLECHECKHELICITYFILTER)
     $   .OR.(HELICITYFILTERLEVEL.EQ.0)
        DO J=1,NCOMB
          DO I=1,NCTAMPS
            GOODAMP(I,J)=.TRUE.
          ENDDO
        ENDDO
        OPEN(1, FILE=LOOPFILTERFN, ERR=100, STATUS='OLD',          
     $    ACTION='READ')
        DO J=1,NCOMB
          READ(1,*,END=101) (GOODAMP(I,J),I=NCTAMPS+1,NLOOPAMPS)
        ENDDO
        GOTO 101
 100    CONTINUE
        FOUNDLOOPFILTER=.FALSE.
        DO J=1,NCOMB
          DO I=NCTAMPS+1,NLOOPAMPS
            GOODAMP(I,J)=(.NOT.USELOOPFILTER)
          ENDDO
        ENDDO
 101    CONTINUE
        CLOSE(1)
        IF (HELICITYFILTERLEVEL.EQ.0) THEN
          FOUNDHELFILTER=.TRUE.
          DO J=1,NCOMB
            GOODHEL(J)=.TRUE.
          ENDDO
          GOTO 122
        ENDIF
        OPEN(1, FILE=HELFILTERFN, ERR=102, STATUS='OLD',          
     $    ACTION='READ')
        READ(1,*,END=103) (GOODHEL(I),I=1,NCOMB)
        GOTO 103
 102    CONTINUE
        FOUNDHELFILTER=.FALSE.
        DO J=1,NCOMB
          GOODHEL(J)=.TRUE.
        ENDDO
 103    CONTINUE
        CLOSE(1)
 122    CONTINUE
      ENDIF

      MP_DONE=.FALSE.
      MP_DONE_ONCE=.FALSE.
      MP_PS_SET=.FALSE.
      STAB_INDEX=0
      DOING_QP_EVALS=.FALSE.
      EVAL_DONE(1)=.TRUE.
      DO I=2,MAXSTABILITYLENGTH
        EVAL_DONE(I)=.FALSE.
      ENDDO

C     Compute the born, for a specific helicity if asked so.
      CALL ML5_0_SMATRIXHEL(P_USER,USERHEL,ANS(0))


      IF (USER_STAB_PREC.GT.0.0D0) THEN
        MLSTABTHRES_BU=MLSTABTHRES
        MLSTABTHRES=USER_STAB_PREC
C       In the initialization, I cannot perform stability test and
C        therefore guarantee any precision
        CTMODEINIT_BU=CTMODEINIT
C       So either one choses quad precision directly
C       CTMODEINIT=4
C       Or, because this is very slow, we keep the orignal value. The
C        accuracy returned is -1 and tells the MC that he should not
C        trust the evaluation for checks.
        CTMODEINIT=CTMODEINIT_BU
      ENDIF

      IF(.NOT.BYPASS_CHECK) THEN
        NTRY=NTRY+1
      ENDIF

      IF(DONEHELDOUBLECHECK.AND.(.NOT.HELDOUBLECHECKED)) THEN
        HELDOUBLECHECKED=.TRUE.
        DONEHELDOUBLECHECK=.FALSE.
      ENDIF

      CHECKPHASE=(NTRY.LE.CHECKCYCLE).AND.(((.NOT.FOUNDLOOPFILTER)
     $ .AND.USELOOPFILTER).OR.(.NOT.FOUNDHELFILTER))

      IF (WRITEOUTFILTERS) THEN
        IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDHELFILTER)) THEN
          OPEN(1, FILE=HELFILTERFN, ERR=110, STATUS='NEW'
     $     ,ACTION='WRITE')
          WRITE(1,*) (GOODHEL(I),I=1,NCOMB)
 110      CONTINUE
          CLOSE(1)
          FOUNDHELFILTER=.TRUE.
        ENDIF

        IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDLOOPFILTER)
     $   .AND.USELOOPFILTER) THEN
          OPEN(1, FILE=LOOPFILTERFN, ERR=111, STATUS='NEW'
     $     ,ACTION='WRITE')
          DO J=1,NCOMB
            WRITE(1,*) (GOODAMP(I,J),I=NCTAMPS+1,NLOOPAMPS)
          ENDDO
 111      CONTINUE
          CLOSE(1)
          FOUNDLOOPFILTER=.TRUE.
        ENDIF
      ENDIF

      IF (BYPASS_CHECK) THEN
        OLD_CHECKPHASE = CHECKPHASE
        OLD_HELDOUBLECHECKED = HELDOUBLECHECKED
        CHECKPHASE = .FALSE.
        HELDOUBLECHECKED = .TRUE.
        DO I=1,NCOMB
          OLD_GOODHEL(I)=GOODHEL(I)
          GOODHEL(I) = .TRUE.
        ENDDO
        DO I=1,NCOMB
          DO J=1,NLOOPAMPS
            OLD_GOODAMP(J,I)=GOODAMP(J,I)
            GOODAMP(J,I) = .TRUE.
          ENDDO
        ENDDO
      ENDIF

      IF(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
        HELPICKED=1
        CTMODE=CTMODEINIT
      ELSE
        IF (USERHEL.NE.-1) THEN
          IF(.NOT.GOODHEL(USERHEL)) THEN
            ANS(1)=0.0D0
            ANS(2)=0.0D0
            ANS(3)=0.0D0
            GOTO 9999
          ENDIF
        ENDIF
        HELPICKED=USERHEL
        IF (CTMODERUN.GT.-1) THEN
          CTMODE=CTMODERUN
        ELSE
          CTMODE=1
        ENDIF
      ENDIF

      DO I=1,NEXTERNAL
        DO J=0,3
          PS(J,I)=P_USER(J,I)
        ENDDO
      ENDDO

      IF (IMPROVEPSPOINT.GE.0) THEN
C       Make the input PS more precise (exact onshell and
C        energy-momentum conservation)
        CALL ML5_0_IMPROVE_PS_POINT_PRECISION(PS)
      ENDIF

      DO I=1,NEXTERNAL
        DO J=0,3
          P(J,I)=PS(J,I)
        ENDDO
      ENDDO

      DO K=1, 3
        BUFFR(K)=0.0D0
        DO I=1,NLOOPAMPS
          AMPL(K,I)=(0.0D0,0.0D0)
        ENDDO
      ENDDO

      LSCALE=DSQRT(ABS((P(0,1)+P(0,2))**2-(P(1,1)+P(1,2))**2-(P(2,1)
     $ +P(2,2))**2-(P(3,1)+P(3,2))**2))

C     We chose to use the born evaluation for the reference
      CALL ML5_0_SMATRIX(P,REF)

 200  CONTINUE

      IF (CTMODE.EQ.0.OR.CTMODE.GE.4) THEN
        CALL MP_UPDATE_AS_PARAM()
      ENDIF

      IF (.NOT.MP_PS_SET.AND.(CTMODE.EQ.0.OR.CTMODE.GE.4)) THEN
        CALL ML5_0_SET_MP_PS(P_USER)
        MP_PS_SET = .TRUE.
      ENDIF

      DO K=1,3
        ANS(K)=0.0D0
      ENDDO

      VALIDH=-1
      DO H=1,NCOMB
        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(
     $.NOT.HELDOUBLECHECKED).OR.GOODHEL(H)))) THEN

C         Handle the possible requirement of specific polarizations
          IF ((.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.POLARIZATIONS(
     $0,0).EQ.0.AND.(.NOT.ML5_0_IS_HEL_SELECTED(H))) THEN
            CYCLE
          ENDIF

          IF (VALIDH.EQ.-1) VALIDH=H
          DO I=1,NEXTERNAL
            NHEL(I)=HELC(I,H)
          ENDDO
C         Check if we are in multiple precision and compute wfs and
C          amps accordingly if needed   
          IF (CTMODE.GE.4) THEN
C           Force that only current helicity is used in the routine
C            below
C           This should always be done, even if MP_DONE is True
C           because the AMPL of the R2 MUST be recomputed for loop
C            induced.
C           (because they are not saved for each hel configuration)
C           (This is not optimal unlike what is done int the loop
C            optimized output)
            HELPICKED_BU = HELPICKED
            HELPICKED = H
            CALL ML5_0_MP_BORN_AMPS_AND_WFS(MP_P)
            HELPICKED = HELPICKED_BU
            GOTO 300
          ENDIF
          CALL IXXXXX(P(0,1),ZERO,NHEL(1),+1*IC(1),W(1,1,H))
          CALL OXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2,H))
          CALL OXXXXX(P(0,3),ZERO,NHEL(3),+1*IC(3),W(1,3,H))
          CALL IXXXXX(P(0,4),ZERO,NHEL(4),-1*IC(4),W(1,4,H))
          CALL VXXXXX(P(0,5),ZERO,NHEL(5),+1*IC(5),W(1,5,H))
          CALL FFV1_2(W(1,1,H),W(1,5,H),GC_5,ZERO,ZERO,W(1,6,H))
          CALL FFV2_3(W(1,4,H),W(1,3,H),GC_47,MDL_MW,MDL_WW,W(1,7,H))
C         Amplitude(s) for born diagram with ID 1
          CALL FFV2_0(W(1,6,H),W(1,2,H),W(1,7,H),GC_47,AMP(1,H))
          CALL FFV1_1(W(1,2,H),W(1,5,H),GC_5,ZERO,ZERO,W(1,8,H))
C         Amplitude(s) for born diagram with ID 2
          CALL FFV2_0(W(1,1,H),W(1,8,H),W(1,7,H),GC_47,AMP(2,H))
          CALL FFV2_1(W(1,2,H),W(1,7,H),GC_47,ZERO,ZERO,W(1,9,H))
C         Counter-term amplitude(s) for loop diagram number 3
          CALL R2_QQ_1_0(W(1,6,H),W(1,9,H),R2_QQQ,AMPL(1,1))
C         Counter-term amplitude(s) for loop diagram number 4
          CALL FFV2_0(W(1,6,H),W(1,2,H),W(1,7,H),R2_SXCW,AMPL(1,2))
C         Counter-term amplitude(s) for loop diagram number 5
          CALL FFV2_0(W(1,1,H),W(1,8,H),W(1,7,H),R2_SXCW,AMPL(1,3))
C         Counter-term amplitude(s) for loop diagram number 7
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,4)
     $     )
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,5)
     $     )
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,6)
     $     )
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,7)
     $     )
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQB,AMPL(1,8))
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,9)
     $     )
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQT,AMPL(1,10))
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
     $     ,11))
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQG_1EPS,AMPL(2
     $     ,12))
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),R2_GQQ,AMPL(1,13))
          CALL FFV2_2(W(1,1,H),W(1,7,H),GC_47,ZERO,ZERO,W(1,10,H))
C         Counter-term amplitude(s) for loop diagram number 11
          CALL R2_QQ_1_0(W(1,10,H),W(1,8,H),R2_QQQ,AMPL(1,14))
C         Counter-term amplitude(s) for loop diagram number 12
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
     $     ,15))
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
     $     ,16))
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
     $     ,17))
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
     $     ,18))
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQB,AMPL(1,19))
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
     $     ,20))
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQT,AMPL(1,21))
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
     $     ,22))
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQG_1EPS,AMPL(2
     $     ,23))
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),R2_GQQ,AMPL(1,24))
C         Amplitude(s) for UVCT diagram with ID 14
          CALL FFV2_0(W(1,6,H),W(1,2,H),W(1,7,H),GC_47,AMPL(1,25))
          AMPL(1,25)=AMPL(1,25)*(1.0D0*UVWFCT_G_2+1.0D0*UVWFCT_G_1)
C         Amplitude(s) for UVCT diagram with ID 15
          CALL FFV2_0(W(1,6,H),W(1,2,H),W(1,7,H),GC_47,AMPL(2,26))
          AMPL(2,26)=AMPL(2,26)*(2.0D0*UVWFCT_G_2_1EPS)
C         Amplitude(s) for UVCT diagram with ID 16
          CALL FFV2_0(W(1,1,H),W(1,8,H),W(1,7,H),GC_47,AMPL(1,27))
          AMPL(1,27)=AMPL(1,27)*(1.0D0*UVWFCT_G_2+1.0D0*UVWFCT_G_1)
C         Amplitude(s) for UVCT diagram with ID 17
          CALL FFV2_0(W(1,1,H),W(1,8,H),W(1,7,H),GC_47,AMPL(2,28))
          AMPL(2,28)=AMPL(2,28)*(2.0D0*UVWFCT_G_2_1EPS)
 300      CONTINUE



          DO I=1,NCTAMPS
            DO J=1,NBORNAMPS
              CFTOT=DCMPLX(CF_N(I,J)/DBLE(ABS(CF_D(I,J))),0.0D0)
              IF(CF_D(I,J).LT.0) CFTOT=CFTOT*IMAG1
              DO K=1,3
                ANS(K)=ANS(K)+2.0D0*DBLE(CFTOT*AMPL(K,I)*DCONJG(AMP(J
     $           ,H)))
              ENDDO
            ENDDO
          ENDDO
        ENDIF
      ENDDO

C     WHEN CTMODE IS >=4, then the MP computation of wfs and amps is
C      automatically done.
      IF (CTMODE.GE.4) THEN
        MP_DONE = .TRUE.
      ENDIF

      IF(SKIPLOOPEVAL) THEN
        GOTO 1226
      ENDIF

C     Loop amplitude for loop diagram with ID 3
      CALL ML5_0_LOOP_2_2(1,6,9,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_5
     $ ,MP__GC_5,1,1,1,29,AMPL(1,29),S(29))
C     Loop amplitude for loop diagram with ID 4
      CALL ML5_0_LOOP_3_3(2,2,7,6,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_47,MP__GC_47,GC_5
     $ ,MP__GC_5,2,1,1,30,AMPL(1,30),S(30))
C     Loop amplitude for loop diagram with ID 5
      CALL ML5_0_LOOP_3_3(3,1,7,8,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_47,MP__GC_47,GC_5
     $ ,MP__GC_5,2,1,1,31,AMPL(1,31),S(31))
C     Loop amplitude for loop diagram with ID 6
      CALL ML5_0_LOOP_4_4(4,1,2,5,7,DCMPLX(ZERO),CMPLX(MP__ZERO
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
     $ ,CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
     $ ,GC_5,MP__GC_5,GC_5,MP__GC_5,GC_5,MP__GC_5,GC_47,MP__GC_47,3,1
     $ ,1,32,AMPL(1,32),S(32))
C     Loop amplitude for loop diagram with ID 7
      CALL ML5_0_LOOP_3_3(5,1,5,9,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_4,MP__GC_4,GC_5
     $ ,MP__GC_5,2,1,1,33,AMPL(1,33),S(33))
C     Loop amplitude for loop diagram with ID 8
      CALL ML5_0_LOOP_4_4(6,1,5,2,7,DCMPLX(ZERO),CMPLX(MP__ZERO
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
     $ ,CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
     $ ,GC_5,MP__GC_5,GC_4,MP__GC_4,GC_5,MP__GC_5,GC_47,MP__GC_47,3,1
     $ ,1,34,AMPL(1,34),S(34))
C     Loop amplitude for loop diagram with ID 9
      CALL ML5_0_LOOP_4_4(7,1,2,7,5,DCMPLX(ZERO),CMPLX(MP__ZERO
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
     $ ,CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
     $ ,GC_5,MP__GC_5,GC_5,MP__GC_5,GC_47,MP__GC_47,GC_5,MP__GC_5,3,1
     $ ,1,35,AMPL(1,35),S(35))
C     Loop amplitude for loop diagram with ID 10
      CALL ML5_0_LOOP_3_3(8,1,5,9,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5,GC_5
     $ ,MP__GC_5,2,1,1,36,AMPL(1,36),S(36))
C     Loop amplitude for loop diagram with ID 11
      CALL ML5_0_LOOP_2_2(9,8,10,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_5
     $ ,MP__GC_5,1,1,1,37,AMPL(1,37),S(37))
C     Loop amplitude for loop diagram with ID 12
      CALL ML5_0_LOOP_3_3(10,2,5,10,DCMPLX(ZERO),CMPLX(MP__ZERO
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_4,MP__GC_4,GC_5
     $ ,MP__GC_5,2,1,1,38,AMPL(1,38),S(38))
C     Loop amplitude for loop diagram with ID 13
      CALL ML5_0_LOOP_3_3(11,2,5,10,DCMPLX(ZERO),CMPLX(MP__ZERO
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5,GC_5
     $ ,MP__GC_5,2,1,1,39,AMPL(1,39),S(39))

      DO I=NCTAMPS+1,NLOOPAMPS
        ANS(1)=ANS(1)+AMPL(1,I)
        ANS(2)=ANS(2)+AMPL(2,I)
        ANS(3)=ANS(3)+AMPL(3,I)
        IF((CTMODERUN.NE.-1).AND..NOT.CHECKPHASE.AND.(.NOT.S(I))) THEN
          WRITE(*,*) '##W03 WARNING Contribution ',I,' is unstable.'
        ENDIF
      ENDDO

 1226 CONTINUE

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
C       Update of NEXTREF, will be used for loop induced only.
        NEXTREF = NEXTREF + ANS(1) + ANS(2) + ANS(3)
        IF((USERHEL.EQ.-1).OR.(USERHEL.EQ.HELPICKED)) THEN
          BUFFR(1)=BUFFR(1)+ANS(1)
          BUFFR(2)=BUFFR(2)+ANS(2)
          BUFFR(3)=BUFFR(3)+ANS(3)
        ENDIF

        IF (CHECKPHASE) THEN
C         SET THE HELICITY FILTER
          IF(.NOT.FOUNDHELFILTER) THEN
            IF(ML5_0_ISZERO(ABS(ANS(1))+ABS(ANS(2))+ABS(ANS(3)),REF
     $       /DBLE(NCOMB),-1)) THEN
              IF(NTRY.EQ.1) THEN
                GOODHEL(HELPICKED)=.FALSE.
              ELSEIF(GOODHEL(HELPICKED)) THEN
                WRITE(*,*) '##W02A WARNING Inconsistent helicity '
     $           ,HELPICKED
                IF(HELINITSTARTOVER) THEN
                  WRITE(*,*) '##I01 INFO Initialization starting over'
     $             //' because of inconsistency in the helicity filter'
     $             //' setup.'
                  NTRY=0
                ENDIF
              ENDIF
            ELSE
              IF(.NOT.GOODHEL(HELPICKED)) THEN
                WRITE(*,*) '##W02B WARNING Inconsistent helicity '
     $           ,HELPICKED
                IF(HELINITSTARTOVER) THEN
                  WRITE(*,*) '##I01 INFO Initialization starting over'
     $             //' because of inconsistency in the helicity filter'
     $             //' setup.'
                  NTRY=0
                ELSE
                  GOODHEL(HELPICKED)=.TRUE.
                ENDIF
              ENDIF
            ENDIF
          ENDIF

C         SET THE LOOP FILTER
          IF(.NOT.FOUNDLOOPFILTER.AND.USELOOPFILTER) THEN
            DO I=NCTAMPS+1,NLOOPAMPS
              IF(.NOT.ML5_0_ISZERO(ABS(AMPL(1,I))+ABS(AMPL(2,I))
     $         +ABS(AMPL(3,I)),(REF*1.0D-4),I)) THEN
                IF(NTRY.EQ.1) THEN
                  GOODAMP(I,HELPICKED)=.TRUE.
                ELSEIF(.NOT.GOODAMP(I,HELPICKED)) THEN
                  WRITE(*,*) '##W02 WARNING Inconsistent loop amp ',I
     $             ,' for helicity ',HELPICKED,'.'
                  IF(LOOPINITSTARTOVER) THEN
                    WRITE(*,*) '##I01 INFO Initialization starting'
     $               //' over because of inconsistency in the loop'
     $               //' filter setup.'
                    NTRY=0
                  ELSE
                    GOODAMP(I,HELPICKED)=.TRUE.
                  ENDIF
                ENDIF
              ENDIF
            ENDDO
          ENDIF
        ELSEIF (.NOT.HELDOUBLECHECKED)THEN
          IF ((.NOT.GOODHEL(HELPICKED)).AND.(.NOT.ML5_0_ISZERO(ABS(ANS(
     $1))+ABS(ANS(2))+ABS(ANS(3)),REF/DBLE(NCOMB),-1))) THEN
            WRITE(*,*) '##W15 Helicity filter could not be'
     $       //' successfully double checked.'
            WRITE(*,*) '##One reason for this is that you have changed'
     $       //' sensible parameters which affected what are the zero'
     $       //' helicity configurations.'
            WRITE(*,*) '##MadLoop will try to reset the Helicity'
     $       //' filter with the next PS points it receives.'
            NTRY=0
            OPEN(30,FILE=HELFILTERFN,ERR=349)
 349        CONTINUE
            CLOSE(30,STATUS='delete')
          ENDIF
C         SET HELDOUBLECHECKED TO .TRUE. WHEN DONE
C         even if it failed we do not want to redo the check
C          afterwards if HELINITSTARTOVER=.FALSE.
          IF (HELPICKED.EQ.NCOMB.AND.(NTRY.NE.0.OR..NOT.HELINITSTARTOVE
     $R)) THEN
            DONEHELDOUBLECHECK=.TRUE.
          ENDIF
        ENDIF

C       GOTO NEXT HELICITY OR FINISH
        IF(HELPICKED.NE.NCOMB) THEN
          HELPICKED=HELPICKED+1
          MP_DONE=.FALSE.
          GOTO 200
        ELSE
          ANS(1)=BUFFR(1)
          ANS(2)=BUFFR(2)
          ANS(3)=BUFFR(3)
C         We add one here to the number of PS points used for building
C          the reference scale for comparison (used only for
C          loop-induced processes).
          NPSPOINTS = NPSPOINTS+1
          IF(NTRY.EQ.0) THEN
            NATTEMPTS=NATTEMPTS+1
            IF(NATTEMPTS.EQ.MAXATTEMPTS) THEN
              WRITE(*,*) '##E01 ERROR Could not initialize the filters'
     $         //' in ',MAXATTEMPTS,' trials'
              STOP
            ENDIF
          ENDIF
        ENDIF

      ENDIF

      DO K=1,3
        ANS(K)=ANS(K)/DBLE(IDEN)
        IF (USERHEL.NE.-1) THEN
          ANS(K)=ANS(K)*HELAVGFACTOR
        ELSE
          DO J=1,NINITIAL
            IF (POLARIZATIONS(J,0).NE.-1) THEN
              ANS(K)=ANS(K)*BEAMS_HELAVGFACTOR(J)
              ANS(K)=ANS(K)/POLARIZATIONS(J,0)
            ENDIF
          ENDDO
        ENDIF
      ENDDO

      IF(.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.(CTMODERUN.LE.-1))
     $  THEN
        STAB_INDEX=STAB_INDEX+1
        IF(DOING_QP_EVALS) THEN
          QP_RES(1,STAB_INDEX)=ANS(1)
          QP_RES(2,STAB_INDEX)=ANS(2)
          QP_RES(3,STAB_INDEX)=ANS(3)
        ELSE
          DP_RES(1,STAB_INDEX)=ANS(1)
          DP_RES(2,STAB_INDEX)=ANS(2)
          DP_RES(3,STAB_INDEX)=ANS(3)
        ENDIF

        IF(DOING_QP_EVALS) THEN
          BASIC_CT_MODE=4
        ELSE
          BASIC_CT_MODE=1
        ENDIF

C       BEGINNING OF THE DEFINITIONS OF THE DIFFERENT EVALUATION
C        METHODS

        IF(.NOT.EVAL_DONE(2)) THEN
          EVAL_DONE(2)=.TRUE.
          CTMODE=BASIC_CT_MODE+1
          GOTO 200
        ENDIF

        CTMODE=BASIC_CT_MODE

        IF(.NOT.EVAL_DONE(3).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
     $.1).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
          EVAL_DONE(3)=.TRUE.
          CALL ML5_0_ROTATE_PS(PS,P,1)
          IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,1)
          GOTO 200
        ENDIF

        IF(.NOT.EVAL_DONE(4).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
     $.2).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
          EVAL_DONE(4)=.TRUE.
          CALL ML5_0_ROTATE_PS(PS,P,2)
          IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,2)
          GOTO 200
        ENDIF

        CALL ML5_0_ROTATE_PS(PS,P,0)
        IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,0)

C       END OF THE DEFINITIONS OF THE DIFFERENT EVALUATION METHODS

        IF(DOING_QP_EVALS) THEN
          CALL ML5_0_COMPUTE_ACCURACY(QP_RES,N_QP_EVAL,ACC,ANS(1))
          ACCURACY(0)=ACC
          RET_CODE_H=3
          IF(ACC.GE.MLSTABTHRES) THEN
            RET_CODE_H=4
            NEPS=NEPS+1
            CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,TEMP1,TEMP)
            WRITE(*,*) '##W03 WARNING An unstable PS point was',      
     $        ' detected.'
            WRITE(*,*) '##(DP,QP) accuracies : (',TEMP1,',',ACC,')'
            WRITE(*,*) '##Best estimate (fin,1eps,2eps) :',(ANS(I),I=1
     $       ,3)
            IF(NEPS.LE.10) THEN
              WRITE(*,*) '##Double precision evaluations :',(DP_RES(1
     $         ,I),I=1,N_DP_EVAL)
              WRITE(*,*) '##Quad   precision evaluations :',(QP_RES(1
     $         ,I),I=1,N_QP_EVAL)
              WRITE(*,*) '##PS point specification :'
              WRITE(*,*) '##Renormalization scale MU_R=',MU_R
              DO I=1,NEXTERNAL
                WRITE (*,'(i2,1x,4e27.17)') I, P(0,I),P(1,I),P(2,I)
     $           ,P(3,I)
              ENDDO
            ENDIF
            IF(NEPS.EQ.10) THEN
              WRITE(*,*) '##Further output of the details of these'
     $         //' unstable PS points will now be suppressed.'
            ENDIF
          ENDIF
        ELSE
          CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,ACC,ANS(1))
          IF(ACC.GE.MLSTABTHRES) THEN
            DOING_QP_EVALS=.TRUE.
            EVAL_DONE(1)=.TRUE.
            DO I=2,MAXSTABILITYLENGTH
              EVAL_DONE(I)=.FALSE.
            ENDDO
            STAB_INDEX=0
            CTMODE=4
            GOTO 200
          ELSE
            ACCURACY(0)=ACC
            RET_CODE_H=2
          ENDIF
        ENDIF
      ELSE
        RET_CODE_H=1
        ACCURACY=-1.0D0
      ENDIF

 9999 CONTINUE

C     Finalize the return code
      IF (MP_DONE_ONCE) THEN
        RET_CODE_T=2
      ELSE
        RET_CODE_T=1
      ENDIF
      IF(CHECKPHASE.OR..NOT.HELDOUBLECHECKED) THEN
        RET_CODE_H=1
        RET_CODE_T=RET_CODE_T+2
        ACCURACY=-1.0D0
      ENDIF
      IF (RET_CODE_H.EQ.4) THEN
        RET_CODE_U=0
      ELSE
        RET_CODE_U=1
      ENDIF

C     Reinitialize the default threshold if it was specified by the
C      user
      IF (USER_STAB_PREC.GT.0.0D0) THEN
        MLSTABTHRES=MLSTABTHRES_BU
        CTMODEINIT=CTMODEINIT_BU
      ENDIF

C     Reinitialize the Lorentz test if it had been disabled because
C      spin-2 particles are in the external states.
      NROTATIONS_DP = NROTATIONS_DP_BU
      NROTATIONS_QP = NROTATIONS_QP_BU

C     Conform to the returned synthax of split orders even though the
C      default output does not support it (this then done only for
C      compatibility purpose).
      ANSRETURNED(0,0)=ANS(0)
      ANSRETURNED(1,0)=ANS(1)
      ANSRETURNED(2,0)=ANS(2)
      ANSRETURNED(3,0)=ANS(3)

C     Reinitialize the check phase logicals and the filters if check
C      bypassed
      IF (BYPASS_CHECK) THEN
        CHECKPHASE = OLD_CHECKPHASE
        HELDOUBLECHECKED = OLD_HELDOUBLECHECKED
        DO I=1,NCOMB
          GOODHEL(I)=OLD_GOODHEL(I)
        ENDDO
        DO I=1,NCOMB
          DO J=1,NLOOPAMPS
            GOODAMP(J,I)=OLD_GOODAMP(J,I)
          ENDDO
        ENDDO
      ENDIF

      END

      SUBROUTINE ML5_0_COMPUTE_ACCURACY(FULLLIST, LENGTH, ACC,
     $  ESTIMATE)
      IMPLICIT NONE
C     
C     PARAMETERS 
C     
      INTEGER MAXSTABILITYLENGTH
      COMMON/ML5_0_STABILITY_TESTS/MAXSTABILITYLENGTH
C     
C     ARGUMENTS 
C     
      REAL*8 FULLLIST(3,MAXSTABILITYLENGTH)
      INTEGER LENGTH
      REAL*8 ACC, ESTIMATE(3)
C     
C     LOCAL VARIABLES 
C     
      LOGICAL MASK(MAXSTABILITYLENGTH)
      LOGICAL MASK3(3)
      DATA MASK3/.TRUE.,.TRUE.,.TRUE./
      INTEGER I,J
      REAL*8 AVG
      REAL*8 DIFF
      REAL*8 ACCURACIES(3)
      REAL*8 LIST(MAXSTABILITYLENGTH)

C     ----------
C     BEGIN CODE
C     ----------
      DO I=1,LENGTH
        MASK(I)=.TRUE.
      ENDDO
      DO I=LENGTH+1,MAXSTABILITYLENGTH
        MASK(I)=.FALSE.
      ENDDO

      DO I=1,3
        DO J=1,MAXSTABILITYLENGTH
          LIST(J)=FULLLIST(I,J)
        ENDDO
        DIFF=MAXVAL(LIST,1,MASK)-MINVAL(LIST,1,MASK)
        AVG=(MAXVAL(LIST,1,MASK)+MINVAL(LIST,1,MASK))/2.0D0
        ESTIMATE(I)=AVG
        IF (AVG.EQ.0.0D0) THEN
          ACCURACIES(I)=DIFF
        ELSE
          ACCURACIES(I)=DIFF/ABS(AVG)
        ENDIF
      ENDDO

C     The technique below is too sensitive, typically to
C     unstablities in very small poles
C     ACC=MAXVAL(ACCURACIES,1,MASK3)
C     The following is used instead
      ACC = 0.0D0
      AVG = 0.0D0
      DO I=1,3
        ACC = ACC + ACCURACIES(I)*ABS(ESTIMATE(I))
        AVG = AVG + ESTIMATE(I)
      ENDDO
      ACC  = ACC / ( ABS(AVG) / 3.0D0)

C     If NaN are present in the evaluation, automatically set the
C      accuracy to 1.0d99.
      DO I=1,3
        DO J=1,MAXSTABILITYLENGTH
          IF (ISNAN(FULLLIST(I,J))) THEN
            ACC = 1.0D99
          ENDIF
        ENDDO
      ENDDO

      END

      SUBROUTINE ML5_0_SET_N_EVALS(N_DP_EVALS,N_QP_EVALS)

      IMPLICIT NONE
      INTEGER N_DP_EVALS, N_QP_EVALS

      INCLUDE 'MadLoopParams.inc'

      IF(CTMODERUN.LE.-1) THEN
        N_DP_EVALS=2+NROTATIONS_DP
        N_QP_EVALS=2+NROTATIONS_QP
      ELSE
        N_DP_EVALS=1
        N_QP_EVALS=1
      ENDIF

      IF(N_DP_EVALS.GT.20.OR.N_QP_EVALS.GT.20) THEN
        WRITE(*,*) '##ERROR:: Increase hardcoded maxstabilitylength.'
        STOP
      ENDIF

      END


C     THIS SUBROUTINE SIMPLY SET THE GLOBAL PS CONFIGURATION GLOBAL
C      VARIABLES FROM A GIVEN VARIABLE IN DOUBLE PRECISION
      SUBROUTINE ML5_0_SET_MP_PS(P)

      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL)
      COMMON/ML5_0_MP_PSPOINT/MP_PS,MP_P
      REAL*8 P(0:3,NEXTERNAL)

      DO I=1,NEXTERNAL
        DO J=0,3
          MP_PS(J,I)=P(J,I)
        ENDDO
      ENDDO
      CALL ML5_0_MP_IMPROVE_PS_POINT_PRECISION(MP_PS)
      DO I=1,NEXTERNAL
        DO J=0,3
          MP_P(J,I)=MP_PS(J,I)
        ENDDO
      ENDDO

      END

      SUBROUTINE ML5_0_SET_COUPLINGORDERS_TARGET(SOTARGET)
      IMPLICIT NONE
C     
C     This routine can be accessed by an external user to set the
C      squared split order target.
C     This functionality is only available in the optimized mode, but
C      for compatibility 
C     purposes, a dummy version is also put in this default output.
C     
C     
C     ARGUMENTS
C     
      INTEGER SOTARGET
C     ----------
C     BEGIN CODE
C     ----------
      WRITE(*,*) '##WARNING:: Ignored, the possibility of selecting'
     $ //' specific squared order contributions is not available in'
     $ //' the default mode.'

      END

      SUBROUTINE ML5_0_FORCE_STABILITY_CHECK(ONOFF)
C     
C     This function can be called by the MadLoop user so as to always
C      have stability
C     checked, even during initialisation, when calling the *_thres
C      routines.
C     
      LOGICAL ONOFF

      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
      DATA BYPASS_CHECK, ALWAYS_TEST_STABILITY /.FALSE.,.FALSE./
      COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY

      ALWAYS_TEST_STABILITY = ONOFF

      END

      SUBROUTINE ML5_0_GET_ANSWER_DIMENSION(ANSDIM)
C     
C     Simple subroutine which returns the upper bound of the second
C      dimension of the
C     quantity ANS(0:3,0:ANSDIM) returned by MadLoop. As long as the
C      default output
C     cannot handle split orders, this ANSDIM will always be 0.
C     
      INCLUDE 'nsquaredSO.inc'

      INTEGER ANSDIM

      ANSDIM=NSQUAREDSO

      END

      SUBROUTINE ML5_0_GET_NSQSO_LOOP(NSQSO)
C     
C     Simple subroutine returning the number of squared split order
C     contributions returned in ANS when calling sloopmatrix 
C     
      INCLUDE 'nsquaredSO.inc'

      INTEGER NSQSO

      NSQSO=NSQUAREDSO

      END

      SUBROUTINE ML5_0_SET_LEG_POLARIZATION(LEG_ID, LEG_POLARIZATION)
      IMPLICIT NONE
C     
C     ARGUMENTS
C     
      INTEGER LEG_ID
      INTEGER LEG_POLARIZATION
C     
C     LOCALS
C     
      INTEGER I
      INTEGER LEG_POLARIZATIONS(0:5)
C     ----------
C     BEGIN CODE
C     ----------

      IF (LEG_POLARIZATION.EQ.-10000) THEN
        LEG_POLARIZATIONS(0)=-1
        DO I=1,5
          LEG_POLARIZATIONS(I)=-10000
        ENDDO
      ELSE
        LEG_POLARIZATIONS(0)=1
        LEG_POLARIZATIONS(1)=LEG_POLARIZATION
        DO I=2,5
          LEG_POLARIZATIONS(I)=-10000
        ENDDO
      ENDIF
      CALL ML5_0_SET_LEG_POLARIZATIONS(LEG_ID,LEG_POLARIZATIONS)

      END

      SUBROUTINE ML5_0_SET_LEG_POLARIZATIONS(LEG_ID, LEG_POLARIZATIONS)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER NPOLENTRIES
      PARAMETER (NPOLENTRIES=(NEXTERNAL+1)*6)
      INTEGER    NCOMB
      PARAMETER (NCOMB=32)
C     
C     ARGUMENTS
C     
      INTEGER LEG_ID
      INTEGER LEG_POLARIZATIONS(0:5)
C     
C     LOCALS
C     
      INTEGER I,J
      LOGICAL ALL_SUMMED_OVER
C     
C     GLOBALS
C     
C     Entry 0 of the first dimension is all -1 if there is no
C      polarization requirement.
C     Then for each leg with ID legID, it is either summed over if
C     POLARIZATIONS(legID,0) is -1, or the list of helicity considered
C      for that
C     leg is POLARIZATIONS(legID,1: POLARIZATIONS(legID,0)   ).
      INTEGER POLARIZATIONS(0:NEXTERNAL,0:5)
      DATA ((POLARIZATIONS(I,J),I=0,NEXTERNAL),J=0,5)/NPOLENTRIES*-1/
      COMMON/ML5_0_BEAM_POL/POLARIZATIONS

      INTEGER BORN_POLARIZATIONS(0:NEXTERNAL,0:5)
      COMMON/ML5_0_BORN_BEAM_POL/BORN_POLARIZATIONS

C     ----------
C     BEGIN CODE
C     ----------

      IF (LEG_POLARIZATIONS(0).EQ.-1) THEN
        DO I=0,5
          POLARIZATIONS(LEG_ID,I)=-1
        ENDDO
      ELSE
        DO I=0,LEG_POLARIZATIONS(0)
          POLARIZATIONS(LEG_ID,I)=LEG_POLARIZATIONS(I)
        ENDDO
        DO I=LEG_POLARIZATIONS(0)+1,5
          POLARIZATIONS(LEG_ID,I)=-10000
        ENDDO
      ENDIF

      ALL_SUMMED_OVER = .TRUE.
      DO I=1,NEXTERNAL
        IF (POLARIZATIONS(I,0).NE.-1) THEN
          ALL_SUMMED_OVER = .FALSE.
          EXIT
        ENDIF
      ENDDO
      IF (ALL_SUMMED_OVER) THEN
        DO I=0,5
          POLARIZATIONS(0,I)=-1
        ENDDO
      ELSE
        DO I=0,5
          POLARIZATIONS(0,I)=0
        ENDDO
      ENDIF

      DO I=0,NEXTERNAL
        DO J=0,5
          BORN_POLARIZATIONS(I,J) = POLARIZATIONS(I,J)
        ENDDO
      ENDDO


      RETURN

      END

      SUBROUTINE ML5_0_SLOOPMATRIXHEL_THRES(P,HEL,ANS,PREC_ASKED
     $ ,PREC_FOUND,RET_CODE)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INCLUDE 'nsquaredSO.inc'
C     
C     ARGUMENTS 
C     
      REAL*8 P(0:3,NEXTERNAL)
      REAL*8 ANS(0:3,0:NSQUAREDSO)
      INTEGER HEL,RET_CODE
      REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
C     
C     GLOBAL VARIABLES
C     
      REAL*8 USER_STAB_PREC
      COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC

      INTEGER I

      INTEGER H,T,U
      REAL*8 ACCURACY(0:NSQUAREDSO)
      COMMON/ML5_0_ACC/ACCURACY,H,T,U

      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
      COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY

C     ----------
C     BEGIN CODE
C     ----------
      USER_STAB_PREC = PREC_ASKED
      CALL ML5_0_SLOOPMATRIXHEL(P,HEL,ANS)
      IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY(0).LT.0.0D0))
     $  THEN
        BYPASS_CHECK = .TRUE.
        CALL ML5_0_SLOOPMATRIXHEL(P,HEL,ANS)
        BYPASS_CHECK = .FALSE.
C       Make sure we correctly return an initialization-type T code
        IF (T.EQ.2) T=4
        IF (T.EQ.1) T=3
      ENDIF

C     Reset it to default value not to affect next runs
      USER_STAB_PREC = -1.0D0
      DO I=0,NSQUAREDSO
        PREC_FOUND(I)=ACCURACY(I)
      ENDDO
      RET_CODE=100*H+10*T+U

      END

      SUBROUTINE ML5_0_SLOOPMATRIX_THRES(P,ANS,PREC_ASKED,PREC_FOUND
     $ ,RET_CODE)
C     
C     Inputs are:
C     P(0:3, Nexternal)  double  :: Kinematic configuration
C      (E,px,py,pz)
C     PEC_ASKED          double  :: Target relative accuracy, -1 for
C      default
C     
C     Outputs are:
C     ANS(3)             double  :: Result (finite, single pole,
C      double pole) 
C     PREC_FOUND         double  :: Relative accuracy estimated for
C      the result
C     Returns -1 if no stab test could be performed.
C     RET_CODE			 integer :: Return code. See below for details
C     
C     Return code conventions: RET_CODE = H*100 + T*10 + U
C     
C     H == 1
C     Stability unknown.
C     H == 2
C     Stable PS (SPS) point.
C     No stability rescue was necessary.
C     H == 3
C     Unstable PS (UPS) point.
C     Stability rescue necessary, and successful.
C     H == 4
C     Exceptional PS (EPS) point.
C     Stability rescue attempted, but unsuccessful.
C     
C     T == 1
C     Default computation (double prec.) was performed.
C     T == 2
C     Quadruple precision was used for this PS point.
C     T == 3
C     MadLoop in initialization phase. Only double precision used.
C     T == 4
C     MadLoop in initialization phase. Quadruple precision used.
C     
C     U is a number left for future use (always set to 0 for now).
C     example: TIR vs OPP usage.
C     
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INCLUDE 'nsquaredSO.inc'
C     
C     ARGUMENTS 
C     
      REAL*8 P(0:3,NEXTERNAL)
      REAL*8 ANS(0:3,0:NSQUAREDSO)
      REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
      INTEGER RET_CODE
C     
C     GLOBAL VARIABLES
C     
      REAL*8 USER_STAB_PREC
      COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC

      INTEGER I

      INTEGER H,T,U
      REAL*8 ACCURACY(0:NSQUAREDSO)
      COMMON/ML5_0_ACC/ACCURACY,H,T,U

      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
      COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY

C     ----------
C     BEGIN CODE
C     ----------
      USER_STAB_PREC = PREC_ASKED
      CALL ML5_0_SLOOPMATRIX(P,ANS)
      IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY(0).LT.0.0D0))
     $  THEN
        BYPASS_CHECK = .TRUE.
        CALL ML5_0_SLOOPMATRIX(P,ANS)
        BYPASS_CHECK = .FALSE.
C       Make sure we correctly return an initialization-type T code
        IF (T.EQ.2) T=4
        IF (T.EQ.1) T=3
      ENDIF

C     Reset it to default value not to affect next runs
      USER_STAB_PREC = -1.0D0
      DO I=0,NSQUAREDSO
        PREC_FOUND(I)=ACCURACY(I)
      ENDDO
      RET_CODE=100*H+10*T+U

      END

