~maddevelopers/mg5amcnlo/3.2.0

C     --=========================================--
C     Main subroutine
C     --=========================================--

      SUBROUTINE ML5_0_SLOOPMATRIX(P_USER,ANS)
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: u u~ > d d~ [ QCD QED ] QED^2=4
C     
C     Modules
C     
      USE ML5_0_POLYNOMIAL_CONSTANTS
C     
      IMPLICIT NONE
C     
C     USER CUSTOMIZABLE OPTIONS
C     
C     This parameter is designed for the check timing command of MG5.
C      It skips the loop reduction.
      LOGICAL SKIPLOOPEVAL
      PARAMETER (SKIPLOOPEVAL=.FALSE.)
C     For timing checks. Stops the code after having only initialized
C      its arrays from the external data files 
      LOGICAL BOOTANDSTOP
      PARAMETER (BOOTANDSTOP=.FALSE.)
      INTEGER TIR_CACHE_SIZE
C     To change memory foot-print of MadLoop, you can change this
C      parameter to be 0,1 or 2 *and recompile*.
C     Notice that this will impact MadLoop speed performances in the
C      context of stability checks.
      INCLUDE 'tir_cache_size.inc'
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=7)
      INTEGER    NLOOPS, NLOOPGROUPS, NCTAMPS
      PARAMETER (NLOOPS=787, NLOOPGROUPS=120, NCTAMPS=235)
      INTEGER    NLOOPAMPS
      PARAMETER (NLOOPAMPS=1022)
      INTEGER    NCOLORROWS
      PARAMETER (NCOLORROWS=NLOOPAMPS)
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER NINITIAL
      PARAMETER (NINITIAL=2)
      INTEGER    NWAVEFUNCS,NLOOPWAVEFUNCS
      PARAMETER (NWAVEFUNCS=18,NLOOPWAVEFUNCS=1215)
      INTEGER    NCOMB
      PARAMETER (NCOMB=16)
      REAL*8     ZERO
      PARAMETER (ZERO=0D0)
      REAL*16     MP__ZERO
      PARAMETER (MP__ZERO=0E0_16)
      COMPLEX*16 IMAG1
      PARAMETER (IMAG1=(0D0,1D0))
C     These are constants related to the split orders
      INCLUDE 'nsqso_born.inc'
      INTEGER    NSO, NSQUAREDSO, NAMPSO
      PARAMETER (NSO=2, NSQUAREDSO=4, NAMPSO=5)
      INTEGER ANS_DIMENSION
      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
      INTEGER NSQSOXNLG
      PARAMETER (NSQSOXNLG=NSQUAREDSO*NLOOPGROUPS)
      INTEGER NSQUAREDSOP1
      PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1)
C     The total number of loop reduction libraries
C     At present, there are only
C      CutTools,PJFry++,IREGI,Golem95,Samurai, Ninja and COLLIER
      INTEGER NLOOPLIB
      PARAMETER (NLOOPLIB=7)
C     Only CutTools or possibly Ninja (if installed with qp support)
C      provide QP
      INTEGER QP_NLOOPLIB
      PARAMETER (QP_NLOOPLIB=1)
      INTEGER MAXSTABILITYLENGTH
      DATA MAXSTABILITYLENGTH/20/
      COMMON/ML5_0_STABILITY_TESTS/MAXSTABILITYLENGTH

C     
C     ARGUMENTS 
C     
      REAL*8 P_USER(0:3,NEXTERNAL)
C     
C     The zeroth component of the second dimension is the result
C      summed over all
C     contributing split orders. The zeroth component of the first one
C      is the Born.
C     Notice that the upper bound of the second integer is not number
C      of squared orders
C     combination for the loops but the maximum between this number
C      for the Born
C     contributions and the loop ones. There are some cases for which
C      the Born contrib.
C     has squared split order contributions than the loop does. For
C      example
C     
C     generate u u~ > d d~ QCD^2<=2 QED^2<=99 [virt=QCD]
C     
C     It is however somehow academical. This is why ANS_DIMENSION is
C      not just NSQSO but rather MAX(NSQSO,NSQSO_BORN)
C     
      REAL*8 ANS(0:3,0:ANS_DIMENSION)
C     
C     LOCAL VARIABLES 
C     
      INTEGER I,J,K,L,H,HEL_MULT,I_QP_LIB,DUMMY

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

      INTEGER CTMODEINIT_BU
      REAL*8 MLSTABTHRES_BU
      INTEGER NEWHELREF
      LOGICAL HEL_INCONSISTENT
      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


      REAL*8 ACC(0:NSQUAREDSO)
      REAL*8 DP_RES(3,0:NSQUAREDSO,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,0:NSQUAREDSO,MAXSTABILITYLENGTH)
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
      INTEGER NATTEMPTS
      DATA NATTEMPTS/0/
      DATA IC/NEXTERNAL*1/
      REAL*8 HELSAVED(3,NCOMB)
      INTEGER ITEMP
      LOGICAL LTEMP
      REAL*8 BORNBUFF(0:NSQSO_BORN),TMPR
      REAL*8 BUFFR(3,0:NSQUAREDSO),BUFFR_BIS(3,0:NSQUAREDSO),TEMP(0:3
     $ ,0:NSQUAREDSO),TEMP1(0:NSQUAREDSO)
      COMPLEX*16 CTEMP
      REAL*8 TEMP2
      COMPLEX*16 COEFS(MAXLWFSIZE,0:VERTEXMAXCOEFS-1,MAXLWFSIZE)
      COMPLEX*16 CFTOT
      LOGICAL FOUNDHELFILTER,FOUNDLOOPFILTER
      DATA FOUNDHELFILTER/.TRUE./
      DATA FOUNDLOOPFILTER/.TRUE./
      LOGICAL LOOPFILTERBUFF(NSQUAREDSO,NLOOPGROUPS)
      DATA ((LOOPFILTERBUFF(J,I),J=1,NSQUAREDSO),I=1,NLOOPGROUPS)
     $ /NSQSOXNLG*.FALSE./

      LOGICAL AUTOMATIC_CACHE_CLEARING
      DATA AUTOMATIC_CACHE_CLEARING/.TRUE./
      COMMON/ML5_0_RUNTIME_OPTIONS/AUTOMATIC_CACHE_CLEARING

      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
      INTEGER OLD_GOODHEL(NCOMB)
      LOGICAL OLD_GOODAMP(NSQUAREDSO,NLOOPGROUPS)
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
      COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
C     
C     FUNCTIONS
C     
      INTEGER ML5_0_TIRCACHE_INDEX
      INTEGER ML5_0_ML5SOINDEX_FOR_BORN_AMP
      INTEGER ML5_0_ML5SOINDEX_FOR_LOOP_AMP
      INTEGER ML5_0_ML5SQSOINDEX
      INTEGER ML5_0_ISSAME
      LOGICAL ML5_0_ISZERO
      LOGICAL ML5_0_IS_HEL_SELECTED
      INTEGER SET_RET_CODE_U
C     
C     GLOBAL VARIABLES
C     
      INCLUDE 'process_info.inc'
      INCLUDE 'unique_id.inc'
      INCLUDE 'coupl.inc'
      INCLUDE 'mp_coupl.inc'
      INCLUDE 'MadLoopParams.inc'


      LOGICAL CHOSEN_SO_CONFIGS(NSQUAREDSO)
      DATA CHOSEN_SO_CONFIGS/.TRUE.,.TRUE.,.TRUE.,.FALSE./
      COMMON/ML5_0_CHOSEN_LOOP_SQSO/CHOSEN_SO_CONFIGS

      INTEGER N_DP_EVAL, N_QP_EVAL
      DATA N_DP_EVAL/1/
      DATA N_QP_EVAL/1/
      COMMON/ML5_0_N_EVALS/N_DP_EVAL,N_QP_EVAL

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

      LOGICAL MP_DONE
      DATA MP_DONE/.FALSE./
      COMMON/ML5_0_MP_DONE/MP_DONE
C     A FLAG TO DENOTE WHETHER THE CORRESPONDING LOOPLIBS ARE
C      AVAILABLE OR NOT
      LOGICAL LOOPLIBS_AVAILABLE(NLOOPLIB)
      DATA LOOPLIBS_AVAILABLE/.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.
     $ ,.TRUE.,.TRUE./
      COMMON/ML5_0_LOOPLIBS_AV/ LOOPLIBS_AVAILABLE
C     A FLAG TO DENOTE WHETHER THE CORRESPONDING DIRECTION TESTS
C      AVAILABLE OR NOT IN THE LOOPLIBS
      LOGICAL LOOPLIBS_DIRECTEST(NLOOPLIB)
      DATA LOOPLIBS_DIRECTEST /.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.
     $ ,.TRUE.,.TRUE./
C     Specifying for which reduction tool quadruple precision is
C      available.
C     The index 0 is dummy and simply means that the corresponding
C      loop_library is not available
C     in which case neither is its quadruple precision version.
      LOGICAL LOOPLIBS_QPAVAILABLE(0:7)
      DATA LOOPLIBS_QPAVAILABLE /.FALSE.,.TRUE.,.FALSE.,.FALSE.
     $ ,.FALSE.,.FALSE.,.FALSE.,.FALSE./
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
      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 MP_PS_SET
      DATA MP_PS_SET/.FALSE./

C     The parameter below sets the convention for the helicity filter
C     For a given helicity, the attached integer 'i' means
C     'i' in ]-inf;-HELOFFSET[ -> Helicity is equal, up to a sign, to
C      helicity number abs(i+HELOFFSET)
C     'i' == -HELOFFSET        -> Helicity is analytically zero
C     'i' in ]-HELOFFSET,inf[  -> Helicity is contributing with weight
C      'i'. If it is zero, it is skipped.
C     Typically, the hel_offset is 10000
      INTEGER HELOFFSET
      DATA HELOFFSET/10000/
      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NSQUAREDSO,NLOOPGROUPS)
      COMMON/ML5_0_FILTERS/GOODAMP,GOODHEL,HELOFFSET

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

C     This integer can be accessed by an external user to set its
C      target squared split order.
C     If set to a value different than -1, the code will try to avoid
C      computing anything which
C     does not contribute to contributions of squared split orders
C      SQSO_TARGET and below.
      INTEGER SQSO_TARGET
      DATA SQSO_TARGET/-1/
      COMMON/ML5_0_SOCHOICE/SQSO_TARGET
C     The following logical are used to broadcast the fact that the
C      target 'required' CT and
C     loop split orders contributions have been reached already and
C      the rest can be skipped.
      LOGICAL UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE,CT_REQ_SO_DONE
     $ ,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE,MP_LOOP_REQ_SO_DONE
     $ ,CTCALL_REQ_SO_DONE,FILTER_SO
      DATA UVCT_REQ_SO_DONE/.FALSE./
      DATA MP_UVCT_REQ_SO_DONE/.FALSE./
      DATA CT_REQ_SO_DONE/.FALSE./
      DATA MP_CT_REQ_SO_DONE/.FALSE./
      DATA LOOP_REQ_SO_DONE/.FALSE./
      DATA MP_LOOP_REQ_SO_DONE/.FALSE./
      DATA CTCALL_REQ_SO_DONE/.FALSE./
      DATA FILTER_SO/.FALSE./
      COMMON/ML5_0_SO_REQS/UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE
     $ ,CT_REQ_SO_DONE,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE
     $ ,MP_LOOP_REQ_SO_DONE,CTCALL_REQ_SO_DONE,FILTER_SO

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

      INTEGER I_SO
      DATA I_SO/1/
      COMMON/ML5_0_I_SO/I_SO
      INTEGER I_LIB
      DATA I_LIB/1/
      COMMON/ML5_0_I_LIB/I_LIB
      LOGICAL QP_TOOLS_AVAILABLE
      DATA QP_TOOLS_AVAILABLE/.FALSE./
      INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB+1)
      COMMON/ML5_0_LOOP_TOOLS/QP_TOOLS_AVAILABLE,INDEX_QP_TOOLS

      COMPLEX*16 AMP(NBORNAMPS)
      COMMON/ML5_0_AMPS/AMP
      COMPLEX*16 W(20,NWAVEFUNCS)
      COMMON/ML5_0_W/W

      COMPLEX*32 MPW(20,NWAVEFUNCS)
      COMMON/ML5_0_MP_W/MPW

      COMPLEX*16 WL(MAXLWFSIZE,0:LOOPMAXCOEFS-1,MAXLWFSIZE,
     $ -1:NLOOPWAVEFUNCS)
      COMPLEX*16 PL(0:3,-1:NLOOPWAVEFUNCS)
      COMMON/ML5_0_WL/WL,PL

      COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
      COMMON/ML5_0_LCOEFS/LOOPCOEFS


      LOGICAL TIR_DONE(NLOOPGROUPS)
      COMMON/ML5_0_TIRCACHING/TIR_DONE

      COMPLEX*16 AMPL(3,NCTAMPS)
      COMMON/ML5_0_AMPL/AMPL

      COMPLEX*16 LOOPRES(3,NSQUAREDSO,NLOOPGROUPS)
      LOGICAL S(NSQUAREDSO,NLOOPGROUPS)
      COMMON/ML5_0_LOOPRES/LOOPRES,S

      INTEGER CF_D(NCOLORROWS,NBORNAMPS)
      INTEGER CF_N(NCOLORROWS,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

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

      CHARACTER(512) MLPATH
      COMMON/MLPATH/MLPATH

C     This is just so that if the user disabled the computation of
C      poles by COLLIER
C     using the MadLoop subroutine, we don't overwrite his choice when
C      reading the parameters
      LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
      DATA  FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION
     $ ,FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION/.FALSE.,.FALSE./
      COMMON/ML5_0_COLLIERPOLESFORCEDCHOICE
     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE

C     This variable controls the general initialization which is
C      *common* between all MadLoop SubProcesses.
C     For example setting the MadLoopPath or reading the ML runtime
C      parameters.
      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./

      LOGICAL WARNED_LORENTZ_STAB_TEST_OFF
      DATA WARNED_LORENTZ_STAB_TEST_OFF/.FALSE./
      INTEGER NROTATIONS_DP_BU,NROTATIONS_QP_BU

      LOGICAL FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
      DATA FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION/.FALSE.,.FALSE./
      COMMON/ML5_0_FPE_IN_REDUCTION/FPE_IN_DP_REDUCTION,
     $  FPE_IN_QP_REDUCTION

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
        ML_INIT = .FALSE.
        CALL PRINT_MADLOOP_BANNER()
        TMP = 'auto'
        CALL SETMADLOOPPATH(TMP)
        CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN)
        CALL MADLOOPPARAMREADER(PARAMFN,.TRUE.)
        IF (FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION) THEN
          COLLIERCOMPUTEUVPOLES = COLLIER_UV_POLE_COMPUTATION_CHOICE
        ENDIF
        IF (FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION) THEN
          COLLIERCOMPUTEIRPOLES = COLLIER_IR_POLE_COMPUTATION_CHOICE
        ENDIF
        IF (FORBID_HEL_DOUBLECHECK) THEN
          DOUBLECHECKHELICITYFILTER = .FALSE.
        ENDIF

C       Make sure that HELFILTERLEVEL is at most 1 if the beam is
C        polarized
        IF (POLARIZATIONS(0,0).EQ.0) THEN
          IF (HELICITYFILTERLEVEL.GT.1) THEN
            WRITE(*,*) '##INFO: When using polarized beam, the'
     $       //' helicity filter of MadLoop can be at most 1. Now'
     $       //' setting HELICITYFILTERLEVEL to 1.'
            HELICITYFILTERLEVEL = 1
          ENDIF
        ENDIF


C       Make sure that NROTATIONS_QP and NROTATIONS_DP are set to zero
C        if AUTOMATIC_CACHE_CLEARING is disabled.
        IF(.NOT.AUTOMATIC_CACHE_CLEARING) THEN
          IF(NROTATIONS_DP.NE.0.OR.NROTATIONS_QP.NE.0) THEN
            WRITE(*,*) '##INFO: AUTOMATIC_CACHE_CLEARING is disabled,'
     $       //' so MadLoop automatically resets NROTATIONS_DP and'
     $       //' NROTATIONS_QP to 0.'
            NROTATIONS_QP=0
            NROTATIONS_DP=0
          ENDIF
        ENDIF

      ENDIF

      IF (LOCAL_ML_INIT) THEN
        LOCAL_ML_INIT = .FALSE.
        QP_TOOLS_AVAILABLE=.FALSE.
        INDEX_QP_TOOLS(1:QP_NLOOPLIB+1)=0
C       SKIP THE ONES THAT NOT AVAILABLE
        J=1
        DO I=1,NLOOPLIB
          IF(MLREDUCTIONLIB(J).EQ.0)EXIT
          IF(.NOT.LOOPLIBS_AVAILABLE(MLREDUCTIONLIB(J)))THEN
            MLREDUCTIONLIB(J:NLOOPLIB-1)=MLREDUCTIONLIB(J+1:NLOOPLIB)
            MLREDUCTIONLIB(NLOOPLIB)=0
          ELSE
            J=J+1
          ENDIF
        ENDDO
        IF(MLREDUCTIONLIB(1).EQ.0)THEN
          STOP 'No available loop reduction lib is provided. Make sure'
     $     //' MLReductionLib is correct.'
        ENDIF
        J=0
        DO I=1,NLOOPLIB
          IF(LOOPLIBS_QPAVAILABLE(MLREDUCTIONLIB(I)))THEN
            J=J+1
            IF(.NOT.QP_TOOLS_AVAILABLE) THEN
              QP_TOOLS_AVAILABLE=.TRUE.
            ENDIF
            INDEX_QP_TOOLS(J)=I
          ENDIF
        ENDDO

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)

        CALL ML5_0_SET_N_EVALS(N_DP_EVAL,N_QP_EVAL)

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,*) 1
          ENDDO
 6116     CONTINUE
          CLOSE(1)
        ENDIF

        OPEN(1, FILE=COLORNUMFN, ERR=104, STATUS='OLD',          
     $    ACTION='READ')
        DO I=1,NCOLORROWS
          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,NCOLORROWS
          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)

C       SETUP OF THE COMMON STARTING EXTERNAL LOOP WAVEFUNCTION
C       IT IS ALSO PS POINT INDEPENDENT, SO IT CAN BE DONE HERE.
C       The index -1 is for the charge-conjugated fermions with
C        flipped fermion flow.
        DO I=0,3
          PL(I,-1)=DCMPLX(0.0D0,0.0D0)
          PL(I,0)=DCMPLX(0.0D0,0.0D0)
        ENDDO
        DO I=1,MAXLWFSIZE
          DO J=0,LOOPMAXCOEFS-1
            DO K=1,MAXLWFSIZE
              WL(I,J,K,-1)=(0.0D0,0.0D0)
              IF(I.EQ.K.AND.J.EQ.0) THEN
                WL(I,J,K,0)=(1.0D0,0.0D0)
              ELSE
                WL(I,J,K,0)=(0.0D0,0.0D0)
              ENDIF
            ENDDO
          ENDDO
        ENDDO
        IF(BOOTANDSTOP) THEN
          WRITE(*,*) '##Stopped by user request.'
          STOP
        ENDIF
      ENDIF

C     This is the chare conjugate version of the unit 4-currents in
C      the canonical cartesian basis.
C     This, for now, is only defined for 4-fermionic currents.
      WL(1,0,2,-1) = DCMPLX(-1.0D0,0.0D0)
      WL(2,0,1,-1) = DCMPLX(1.0D0,0.0D0)
      WL(3,0,4,-1) = DCMPLX(1.0D0,0.0D0)
      WL(4,0,3,-1) = DCMPLX(-1.0D0,0.0D0)

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 = .TRUE.
        ENDIF
        NROTATIONS_QP=0
        NROTATIONS_DP=0
        CALL ML5_0_SET_N_EVALS(N_DP_EVAL,N_QP_EVAL)
      ENDIF

      IF(NTRY.EQ.0) THEN
        HELDOUBLECHECKED=(.NOT.DOUBLECHECKHELICITYFILTER)
     $   .OR.(HELICITYFILTERLEVEL.EQ.0)
        OPEN(1, FILE=LOOPFILTERFN, ERR=100, STATUS='OLD',          
     $    ACTION='READ')
        DO J=1,NLOOPGROUPS
          READ(1,*,END=101) (GOODAMP(I,J),I=1,NSQUAREDSO)
        ENDDO
        GOTO 101
 100    CONTINUE
        FOUNDLOOPFILTER=.FALSE.
        DO J=1,NLOOPGROUPS
          DO I=1,NSQUAREDSO
            GOODAMP(I,J)=(.NOT.USELOOPFILTER)
          ENDDO
        ENDDO
 101    CONTINUE
        CLOSE(1)

        IF (.NOT.USELOOPFILTER) THEN
          DO J=1,NLOOPGROUPS
            DO I=1,NSQUAREDSO
              GOODAMP(I,J)=.TRUE.
            ENDDO
          ENDDO
        ENDIF

        IF (HELICITYFILTERLEVEL.EQ.0) THEN
          FOUNDHELFILTER=.TRUE.
          DO J=1,NCOMB
            GOODHEL(J)=1
          ENDDO
          GOTO 122
        ENDIF
        OPEN(1, FILE=HELFILTERFN, ERR=102, STATUS='OLD',          
     $    ACTION='READ')
        DO I=1,NCOMB
          READ(1,*,END=103) GOODHEL(I)
        ENDDO
        GOTO 103
 102    CONTINUE
        FOUNDHELFILTER=.FALSE.
        DO J=1,NCOMB
          GOODHEL(J)=1
        ENDDO
 103    CONTINUE
        CLOSE(1)
        IF (HELICITYFILTERLEVEL.EQ.1) THEN
C         We must make sure to remove the matching-helicity
C          optimisation, as requested by the user.
          DO J=1,NCOMB
            IF ((GOODHEL(J).GT.1).OR.(GOODHEL(J).LT.-HELOFFSET)) THEN
              GOODHEL(J)=1
            ENDIF
          ENDDO
        ENDIF
 122    CONTINUE
      ENDIF

C     First compute the borns, it will store them in ANS(0,I)
C     It is left untouched for the rest of MadLoop evaluation.
C     Notice that the squared split order index I does NOT
C     correspond to the same ordering of J for the loop ME 
C     results stored in ANS(K,J), with K in [1-3].The ordering 
C     of each can be obtained with ML5SOINDEX_FOR_SQUARED_ORDERS 
C     and SQSOINDEX_FROM_ORDERS for the loop ME and born ME 
C     respectively. For this to work, we assume that there is 
C     always more squared split orders in the loop ME than in the
C     born ME, which is practically always true. In any case, only
C     the split_order summed value I=0 is used in ML5 code.
      DO I=0,NSQSO_BORN
        BORNBUFF(I)=0.0D0
      ENDDO
      CALL ML5_0_SMATRIXHEL_SPLITORDERS(P_USER,USERHEL,BORNBUFF(0))
      DO I=0,NSQSO_BORN
        ANS(0,I)=BORNBUFF(I)
      ENDDO

C     We set here the reference to the born summed over all split
C      orders
      REF=0.0D0
      DO I=1,NSQSO_BORN
        REF=REF+ANS(0,I)
      ENDDO

      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

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

      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(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 ((HELICITYFILTERLEVEL.NE.0).AND.(.NOT. CHECKPHASE)
     $   .AND.(.NOT.FOUNDHELFILTER)) THEN
          OPEN(1, FILE=HELFILTERFN, ERR=110, STATUS='NEW'
     $     ,ACTION='WRITE')
          DO I=1,NCOMB
            WRITE(1,*) GOODHEL(I)
          ENDDO
 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,NLOOPGROUPS
            WRITE(1,*) (GOODAMP(I,J),I=1,NSQUAREDSO)
          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)=1
        ENDDO
        DO I=1,NSQUAREDSO
          DO J=1,NLOOPGROUPS
            OLD_GOODAMP(I,J)=GOODAMP(I,J)
            GOODAMP(I,J)=.TRUE.
          ENDDO
        ENDDO
      ENDIF

      IF(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
        HELPICKED=1
        CTMODE=CTMODEINIT
      ELSE
        IF (USERHEL.NE.-1) THEN
          IF(GOODHEL(USERHEL).EQ.-HELOFFSET) THEN
            DO I=0,NSQUAREDSO
              ANS(1,I)=0.0D0
              ANS(2,I)=0.0D0
              ANS(3,I)=0.0D0
            ENDDO
            GOTO 9999
          ENDIF
        ENDIF
        HELPICKED=USERHEL
        IF (CTMODERUN.NE.-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

C     Make sure we start with empty caches
      IF (AUTOMATIC_CACHE_CLEARING) THEN
        CALL ML5_0_CLEAR_CACHES()
      ENDIF

C     Now make sure to turn on the global COLLIER cache if applicable
      CALL ML5_0_SET_COLLIER_GLOBAL_CACHE(.TRUE.)

      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
        DO I=0,NSQUAREDSO
          BUFFR(K,I)=0.0D0
        ENDDO
        DO I=1,NCTAMPS
          AMPL(K,I)=(0.0D0,0.0D0)
        ENDDO
      ENDDO

C     Start by using the first available loop reduction library and qp
C      library.
      I_LIB=1
      I_QP_LIB=1

      GOTO 208
C     MadLoop jumps to this label during stability checks when it
C      recomputes a rotated PS point
 200  CONTINUE
C     For the computation of a rotated version of this PS point we
C      must reset the all MadLoop cache since this changes the
C      definition of the loop denominators.
C     We don't check for AUTOMATIC_CACHE_CLEARING here because the
C      Lorentz test should anyway be disabled if the flag is turned
C      off.
      CALL ML5_0_CLEAR_CACHES()
 208  CONTINUE

C     MadLoop jumps to this label during initialization when it goes
C      to the computation of the next helicity.
 205  CONTINUE

      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

      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))

      CTCALL_REQ_SO_DONE=.FALSE.
      FILTER_SO = (.NOT.CHECKPHASE)
     $ .AND.HELDOUBLECHECKED.AND.(SQSO_TARGET.NE.-1)

      DO I=1,NLOOPGROUPS
        DO J=0,LOOPMAXCOEFS-1
          DO K=1,NSQUAREDSO
            LOOPCOEFS(J,K,I)=(0.0D0,0.0D0)
          ENDDO
        ENDDO
      ENDDO

      DO I=1,NLOOPGROUPS
        DO J=1,3
          DO K=1,NSQUAREDSO
            LOOPRES(J,K,I)=(0.0D0,0.0D0)
          ENDDO
        ENDDO
      ENDDO

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

C     Check if we directly go to multiple precision
      IF (CTMODE.GE.4) THEN
        IF (.NOT.MP_DONE) THEN
          CALL ML5_0_MP_COMPUTE_LOOP_COEFS(MP_P,BUFFR_BIS)
C         It should be safe to directly set MP_DONE to true already
C          here. But maybe I overlooked something.
          MP_DONE=.TRUE.
        ENDIF
C       Even if MP_DONE is .TRUE. we should anyway skip the
C       double precision evaluation as it as already been
C       computed in quadruple precision.
        GOTO 300
      ENDIF

      DO H=1,NCOMB
        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1)
     $   .AND.(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H)
     $   .GT.-HELOFFSET.AND.GOODHEL(H).NE.0)))) 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

          DO I=1,NEXTERNAL
            NHEL(I)=HELC(I,H)
          ENDDO

          UVCT_REQ_SO_DONE=.FALSE.
          CT_REQ_SO_DONE=.FALSE.
          LOOP_REQ_SO_DONE=.FALSE.

          IF (.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.HELPICKED.EQ.-1)
     $      THEN
            HEL_MULT=GOODHEL(H)
          ELSE
            HEL_MULT=1
          ENDIF


C         Helas calls for the born amplitudes and counterterms
C          associated to given loops
          CALL ML5_0_HELAS_CALLS_AMPB_1(P,NHEL,H,IC)
 2000     CONTINUE
          CT_REQ_SO_DONE=.TRUE.

C         Helas calls for the counterterm of type 'UVtree' in the UFO.
C         These are generated irrespectively of the produced loops.
C         In general, only wavefunction renormalization counterterms
C         (if needed by the loop UFO model) are of this type.
C         Quite often and in principle for all loop UFO models from 
C         FeynRules, there are none of these type of counterterms.
          CALL ML5_0_HELAS_CALLS_UVCT_1(P,NHEL,H,IC)
 3000     CONTINUE
          UVCT_REQ_SO_DONE=.TRUE.

          DO J=1,NBORNAMPS
            CTEMP = 2.0D0*HEL_MULT*DCONJG(AMP(J))
            DO I=1,NCTAMPS
              CFTOT=DCMPLX(CF_N(I,J)/DBLE(ABS(CF_D(I,J))),0.0D0)
              IF(CF_D(I,J).LT.0) CFTOT=CFTOT*IMAG1
              ITEMP =
     $          ML5_0_ML5SQSOINDEX(ML5_0_ML5SOINDEX_FOR_LOOP_AMP(I)
     $         ,ML5_0_ML5SOINDEX_FOR_BORN_AMP(J))
              IF (.NOT.FILTER_SO.OR.SQSO_TARGET.EQ.ITEMP) THEN
                DO K=1,3
                  TEMP2 = DBLE(CFTOT*AMPL(K,I)*CTEMP)
                  ANS(K,ITEMP)=ANS(K,ITEMP)+TEMP2
                  ANS(K,0)=ANS(K,0)+TEMP2
                ENDDO
              ENDIF
            ENDDO
          ENDDO

          CALL ML5_0_COEF_CONSTRUCTION_1(P,NHEL,H,IC)
          CALL ML5_0_COEF_CONSTRUCTION_2(P,NHEL,H,IC)
 4000     CONTINUE
          LOOP_REQ_SO_DONE=.TRUE.



        ENDIF
      ENDDO

C     Grouping of loop diagrams now done directly when creating the
C      LOOPCOEFS.

      DO I=0,NSQUAREDSO
        DO J=1,3
          BUFFR_BIS(J,I)=ANS(J,I)
        ENDDO
      ENDDO


C     MadLoop jumps to this label during stability checks when it
C      recomputes the same PS point with a different CTMode
 300  CONTINUE

C     Make sure that the loop calls are performed since this is new
C      evaluation.
      CTCALL_REQ_SO_DONE=.FALSE.

C     Free cache when using IREGI
      IF(IREGIRECY.AND.MLREDUCTIONLIB(I_LIB).EQ.3) THEN
        CALL IREGI_FREE_PS()
      ENDIF

C     Even if the user did ask to turn off the automatic TIR cache
C      clearing, we must do it now if the CTModeIndex rolls over the
C      size of the TIR cache employed.
C     Notice that we must do that only when processing a new CT mode
C      as part of the stability test and not when computing a new
C      helicity as part of the filtering process.
C     This we check that we are not in the initialization phase.
C     If we are not in CTModeRun=-1, then we never need to clear the
C      cache since the TIR will always be used for a unique
C      computation (not stab test).
C     Also, it is clear that if we are running OPP when reaching this'
C     //' line, then we shouldn't clear the TIR cache as it might
C      still be useful later.
C     Finally, notice that the conditional statement below should
C      never be true except you have TIR library supporting quadruple
C      precision or when TIR_CACHE_SIZE<2.
      IF((.NOT.CHECKPHASE.AND.(HELDOUBLECHECKED)).AND.CTMODERUN.EQ.
     $ -1.AND.(MLREDUCTIONLIB(I_LIB).NE.1.AND.MLREDUCTIONLIB(I_LIB)
     $ .NE.5).AND.(ML5_0_TIRCACHE_INDEX(CTMODE).EQ.(TIR_CACHE_SIZE+1)))
     $  THEN
        CALL ML5_0_CLEAR_TIR_CACHE()
      ENDIF


      DO I=0,NSQUAREDSO
        DO J=1,3
          ANS(J,I)=BUFFR_BIS(J,I)
        ENDDO
      ENDDO


      IF(SKIPLOOPEVAL.OR.(.NOT.LOOP_REQ_SO_DONE.AND..NOT.MP_LOOP_REQ_SO
     $_DONE)) THEN
        GOTO 1226
      ENDIF

      DO I_SO=1,NSQUAREDSO
        DO J=1,NLOOPGROUPS
          S(I_SO,J)=.TRUE.
        ENDDO
        IF (FILTER_SO.AND.SQSO_TARGET.NE.I_SO) GOTO 5001
        CALL ML5_0_LOOP_CT_CALLS_1(P,NHEL,H,IC)
        GOTO 5001
 5000   CONTINUE
        CTCALL_REQ_SO_DONE=.TRUE.
 5001   CONTINUE
      ENDDO

      DO I=1,NLOOPGROUPS
        LTEMP=.TRUE.
        DO K=1,NSQUAREDSO
          IF (.NOT.FILTER_SO.OR.SQSO_TARGET.EQ.K) THEN
            IF (.NOT.S(K,I)) LTEMP=.FALSE.
            DO J=1,3
              ANS(J,K)=ANS(J,K)+LOOPRES(J,K,I)
              ANS(J,0)=ANS(J,0)+LOOPRES(J,K,I)
            ENDDO
          ENDIF
        ENDDO
        IF((CTMODERUN.NE.-1).AND..NOT.CHECKPHASE.AND.(.NOT.LTEMP)) THEN
          WRITE(*,*) '##W03 WARNING Contribution ',I,' is unstable.'
        ENDIF
      ENDDO

C     Make sure that no NaN is present in the result
      DO K=1,NSQUAREDSO
        DO J=1,3
          IF (.NOT.(ANS(J,K).EQ.ANS(J,K))) THEN
            IF (DOING_QP_EVALS) THEN
              FPE_IN_QP_REDUCTION = .TRUE.
            ELSE
              FPE_IN_DP_REDUCTION = .TRUE.
            ENDIF
          ENDIF
        ENDDO
      ENDDO

 1226 CONTINUE

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
        IF((USERHEL.EQ.-1).OR.(USERHEL.EQ.HELPICKED)) THEN
C         Make sure that that no polarization constraint filters out
C          this helicity
          IF (POLARIZATIONS(0,0).EQ.
     $     -1.OR.ML5_0_IS_HEL_SELECTED(HELPICKED)) THEN
C           TO KEEP TRACK OF THE FINAL ANSWER TO BE RETURNED DURING
C            CHECK PHASE
            DO I=0,NSQUAREDSO
              DO K=1,3
                BUFFR(K,I)=BUFFR(K,I)+ANS(K,I)
              ENDDO
            ENDDO
          ENDIF
        ENDIF
C       SAVE RESULT OF EACH INDEPENDENT HELICITY FOR COMPARISON DURING
C        THE HELICITY FILTER SETUP
        HELSAVED(1,HELPICKED)=0.0D0
        HELSAVED(2,HELPICKED)=0.0D0
        HELSAVED(3,HELPICKED)=0.0D0
        DO I=1,NSQUAREDSO
          IF (CHOSEN_SO_CONFIGS(I)) THEN
            HELSAVED(1,HELPICKED)=HELSAVED(1,HELPICKED)+ANS(1,I)
            HELSAVED(2,HELPICKED)=HELSAVED(2,HELPICKED)+ANS(2,I)
            HELSAVED(3,HELPICKED)=HELSAVED(3,HELPICKED)+ANS(3,I)
          ENDIF
        ENDDO

        IF (CHECKPHASE.AND.NTRY.NE.0) THEN
C         SET THE HELICITY FILTER
          IF(.NOT.FOUNDHELFILTER) THEN
            HEL_INCONSISTENT=.FALSE.
            IF(ML5_0_ISZERO(DABS(HELSAVED(1,HELPICKED))
     $       +DABS(HELSAVED(2,HELPICKED))+DABS(HELSAVED(3,HELPICKED))
     $       ,REF/DBLE(NCOMB),-1,-1)) THEN
              IF(NTRY.EQ.1) THEN
                GOODHEL(HELPICKED)=-HELOFFSET
              ELSEIF(GOODHEL(HELPICKED).NE.-HELOFFSET) THEN
                WRITE(*,*) '##W02A WARNING Inconsistent zero helicity'
     $           //' ',HELPICKED
                IF(HELINITSTARTOVER) THEN
                  WRITE(*,*) '##I01 INFO Initialization starting over'
     $             //' because of inconsistency in the helicity filter'
     $             //' setup.'
                  NTRY=0
                ELSE
                  HEL_INCONSISTENT=.TRUE.
                ENDIF
              ENDIF
            ELSEIF(HELICITYFILTERLEVEL.GT.1) THEN
              DO H=1,HELPICKED-1
                IF(GOODHEL(H).GT.-HELOFFSET) THEN
C                 Be looser for helicity check, bring a factor 100
                  DUMMY=ML5_0_ISSAME(HELSAVED(1,HELPICKED),HELSAVED(1
     $             ,H),REF,.FALSE.)
                  IF(DUMMY.NE.0) THEN
                    IF(NTRY.EQ.1) THEN
C                     Set the matching helicity to be contributing
C                      once more
                      GOODHEL(H)=GOODHEL(H)+DUMMY
C                     Use an offset to clearly show it is linked to an
C                      other one and to avoid overlap
                      GOODHEL(HELPICKED)=-H-HELOFFSET
C                     Make sure we have paired this hel config to the
C                      same one last PS point
                    ELSEIF(GOODHEL(HELPICKED).NE.(-H-HELOFFSET)) THEN
                      WRITE(*,*) '##W02B WARNING Inconsistent matching'
     $                 //' helicity ',HELPICKED
                      IF(HELINITSTARTOVER) THEN
                        WRITE(*,*) '##I01 INFO Initialization starting'
     $                   //' over because of inconsistency in the'
     $                   //' helicity filter setup.'
                        NTRY=0
                      ELSE
                        HEL_INCONSISTENT=.TRUE.
                      ENDIF
                    ENDIF
                  ENDIF
                ENDIF
              ENDDO
            ENDIF
            IF(HEL_INCONSISTENT) THEN
C             This helicity has unstable filter so we will always
C              compute it by itself.
C             We therefore also need to remove it from the
C              multiplicative factor of the corresponding helicity.
              IF(GOODHEL(HELPICKED).LT.-HELOFFSET) THEN
                GOODHEL(-GOODHEL(HELPICKED)-HELOFFSET)=GOODHEL(
     $           -GOODHEL(HELPICKED)-HELOFFSET)-1
              ENDIF
C             If several helicities were matched to that one, we need
C              to chose another one as reference and redirect the
C              others to this new one
C             Of course if it is one, then we do not need to do
C              anything (because with HELINITSTARTOVER=.FALSE. we only
C              support exactly identical Hels.)
              IF(GOODHEL(HELPICKED).GT.
     $         -HELOFFSET.AND.GOODHEL(HELPICKED).NE.1) THEN
                NEWHELREF=-1
                DO H=1,NCOMB
                  IF (GOODHEL(H).EQ.(-HELOFFSET-HELPICKED)) THEN
                    IF (NEWHELREF.EQ.-1) THEN
                      NEWHELREF=H
                      GOODHEL(H)=GOODHEL(HELPICKED)-1
                    ELSE
                      GOODHEL(H)=-NEWHELREF-HELOFFSET
                    ENDIF
                  ENDIF
                ENDDO
              ENDIF
C             In all cases, from now on this helicity will be computed
C              independantly of the others.
C             In particular, it is the only thing to do if the
C              helicity was flagged not contributing.
              GOODHEL(HELPICKED)=1
            ENDIF
          ENDIF

C         SET THE LOOP FILTER
          IF(.NOT.FOUNDLOOPFILTER.AND.USELOOPFILTER) THEN
            DO I=1,NLOOPGROUPS
              DO J=1,NSQUAREDSO
                IF(.NOT.ML5_0_ISZERO(ABS(LOOPRES(1,J,I))+ABS(LOOPRES(2
     $           ,J,I))+ABS(LOOPRES(3,J,I)),(REF*1.0D-4),I,J)) THEN
                  IF(NTRY.EQ.1) THEN
                    GOODAMP(J,I)=.TRUE.
                    LOOPFILTERBUFF(J,I)=.TRUE.
                  ELSEIF(.NOT.LOOPFILTERBUFF(J,I)) THEN
                    WRITE(*,*) '##W02 WARNING Inconsistent loop amp '
     $               ,I,'.'
                    IF(LOOPINITSTARTOVER) THEN
                      WRITE(*,*) '##I01 INFO Initialization starting'
     $                 //' over because of inconsistency in the loop'
     $                 //' filter setup.'
                      NTRY=0
                    ELSE
                      GOODAMP(J,I)=.TRUE.
                    ENDIF
                  ENDIF
                ENDIF
              ENDDO
            ENDDO
          ENDIF
        ELSEIF (.NOT.HELDOUBLECHECKED.AND.NTRY.NE.0)THEN
C         DOUBLE CHECK THE HELICITY FILTER
          IF (GOODHEL(HELPICKED).EQ.-HELOFFSET) THEN
            IF (.NOT.ML5_0_ISZERO(DABS(HELSAVED(1,HELPICKED))
     $       +DABS(HELSAVED(2,HELPICKED))+DABS(HELSAVED(2,HELPICKED))
     $       ,REF/DBLE(NCOMB),-1,-1)) THEN
              WRITE(*,*) '##W15 Helicity filter could not be'
     $         //' successfully double checked.'
              WRITE(*,*) '##One reason for this is that you might 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(29,FILE=HELFILTERFN,ERR=348)
 348          CONTINUE
              CLOSE(29,STATUS='delete')
            ENDIF
          ENDIF
          IF (GOODHEL(HELPICKED).LT.-HELOFFSET.AND.NTRY.NE.0) THEN
            IF(ML5_0_ISSAME(HELSAVED(1,HELPICKED),HELSAVED(1
     $       ,ABS(GOODHEL(HELPICKED)+HELOFFSET)),REF,.TRUE.).EQ.0) THEN
              WRITE(*,*) '##W15 Helicity filter could not be'
     $         //' successfully double checked.'
              WRITE(*,*) '##One reason for this is that you might have'
     $         //' changed sensible parameters which affected the'
     $         //' helicity dependance relations.'
              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
          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 205
        ELSE
C         Useful printout
C         do I=1,NCOMB
C         write(*,*) 'HELSAVED(1,',I,')=',HELSAVED(1,I)
C         write(*,*) 'HELSAVED(2,',I,')=',HELSAVED(2,I)
C         write(*,*) 'HELSAVED(3,',I,')=',HELSAVED(3,I)
C         write(*,*) '   GOODHEL(',I,')=',GOODHEL(I)
C         ENDDO
          DO I=0,NSQUAREDSO
            DO K=1,3
              ANS(K,I)=BUFFR(K,I)
            ENDDO
          ENDDO
          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 1
            ENDIF
          ENDIF
        ENDIF
      ENDIF

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


      IF(.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.(CTMODERUN.EQ.-1))
     $  THEN
        STAB_INDEX=STAB_INDEX+1
        IF(DOING_QP_EVALS.AND.LOOPLIBS_QPAVAILABLE(MLREDUCTIONLIB(I_LIB)
     $   )) THEN
C         Only run over the reduction algorithms which support
C          quadruple precision
          DO I=0,NSQUAREDSO
            DO K=1,3
              QP_RES(K,I,STAB_INDEX)=ANS(K,I)
            ENDDO
          ENDDO
        ELSE
          DO I=0,NSQUAREDSO
            DO K=1,3
              DP_RES(K,I,STAB_INDEX)=ANS(K,I)
            ENDDO
          ENDDO
        ENDIF

        IF(DOING_QP_EVALS.AND.LOOPLIBS_QPAVAILABLE(MLREDUCTIONLIB(I_LIB)
     $   )) 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.
          IF(LOOPLIBS_DIRECTEST(MLREDUCTIONLIB(I_LIB)))THEN
            CTMODE=BASIC_CT_MODE+1
            GOTO 300
          ELSE
C           If some TIR library would not support the loop direction
C            test (they all do for now), then we would just copy the
C            answer from mode 1 and carry on.
            STAB_INDEX=STAB_INDEX+1
            IF(DOING_QP_EVALS)THEN
              DO I=0,NSQUAREDSO
                DO K=1,3
                  QP_RES(K,I,STAB_INDEX)=ANS(K,I)
                ENDDO
              ENDDO
            ELSE
              DO I=0,NSQUAREDSO
                DO K=1,3
                  DP_RES(K,I,STAB_INDEX)=ANS(K,I)
                ENDDO
              ENDDO
            ENDIF
          ENDIF
        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.AND.LOOPLIBS_QPAVAILABLE(MLREDUCTIONLIB(I_LIB)
     $   )) THEN
          CALL ML5_0_COMPUTE_ACCURACY(QP_RES,N_QP_EVAL,ACC,ANS)
C         If a floating point exception was encountered during the
C          reduction,
C         the result cannot be trusted at all and we hardset all
C          accuracies to 1.0
          IF(FPE_IN_QP_REDUCTION) THEN
            DO I=0,NSQUAREDSO
              ACC(I)=1.0D0
            ENDDO
          ENDIF
          DO I=0,NSQUAREDSO
            ACCURACY(I)=ACC(I)
          ENDDO
          RET_CODE_H=3
          RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.TRUE.
     $     ,.TRUE.)
          IF(MAXVAL(ACC).GE.MLSTABTHRES) THEN
            I_QP_LIB=I_QP_LIB+1
            IF(I_QP_LIB.GT.QP_NLOOPLIB.OR.INDEX_QP_TOOLS(I_QP_LIB)
     $       .EQ.0)THEN
              RET_CODE_H=4
              RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.TRUE.
     $         ,.FALSE.)
              NEPS=NEPS+1
              CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,TEMP1,TEMP)
              CALL ML5_0_COMPUTE_ACCURACY(QP_RES,N_QP_EVAL,ACC,ANS)
              IF(NEPS.LE.10) THEN
                WRITE(*,*) '##W03 WARNING An unstable PS point was',  
     $                ' detected.'
                IF(FPE_IN_QP_REDUCTION) THEN
                  WRITE(*,*) '## The last QP reduction was deemed'
     $             //' unstable because a floating point exception was'
     $             //' encountered.'
                ENDIF
                IF (NSQUAREDSO.NE.1) THEN
                  WRITE(*,*) '##Accuracies for each split order,'
     $             //' starting with the summed case'
                  WRITE(*,*) '##DP accuracies (for each split order):'
     $             //' ',(TEMP1(I),I=0,NSQUAREDSO)
                  WRITE(*,*) '##QP accuracies (for each split order):'
     $             //' ',(ACC(I),I=0,NSQUAREDSO)
                ELSE
                  WRITE(*,*) '##DP accuracy: ',TEMP1(1)
                  WRITE(*,*) '##QP accuracy: ',ACC(1)
                ENDIF
                DO J=0,NSQUAREDSO
                  IF (NSQUAREDSO.NE.1.OR.J.NE.0) THEN
                    IF (J.EQ.0) THEN
                      WRITE(*,*) 'Details for all split orders summed'
     $                 //' :'
                    ELSE
                      WRITE(*,*) 'Details for split order index : ',J
                    ENDIF
                    WRITE(*,*) 'Best estimate (fin,1eps,2eps):',(ANS(I
     $               ,J),I=1,3)
                    WRITE(*,*) 'Finite double precision evaluations :'
     $               ,(DP_RES(1,J,I),I=1,N_DP_EVAL)
                    WRITE(*,*) 'Finite quad   precision evaluations :'
     $               ,(QP_RES(1,J,I),I=1,N_QP_EVAL)
                  ENDIF
                ENDDO
                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
            ELSE
C             A new reduction tool will be used. Reinitialize the FPE
C              flags.
              FPE_IN_DP_REDUCTION=.FALSE.
              FPE_IN_QP_REDUCTION=.FALSE.
              I_LIB=INDEX_QP_TOOLS(I_QP_LIB)
              EVAL_DONE(1)=.TRUE.
              DO I=2,MAXSTABILITYLENGTH
                EVAL_DONE(I)=.FALSE.
              ENDDO
              STAB_INDEX=0
              IF(NROTATIONS_QP.GE.1)THEN
                GOTO 200
              ELSE
                GOTO 300
              ENDIF
            ENDIF
          ENDIF
        ELSEIF(.NOT.DOING_QP_EVALS)THEN
          CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,ACC,ANS)
C         If a floating point exception was encountered during the
C          reduction,
C         the result cannot be trusted at all and we hardset all
C          accuracies to 1.0
          IF(FPE_IN_DP_REDUCTION) THEN
            DO I=0,NSQUAREDSO
              ACC(I)=1.0D0
            ENDDO
          ENDIF
          IF(MAXVAL(ACC).GE.MLSTABTHRES) THEN
            I_LIB=I_LIB+1
            IF((I_LIB.GT.NLOOPLIB.OR.MLREDUCTIONLIB(I_LIB).EQ.0)
     $       .AND.QP_TOOLS_AVAILABLE)THEN
              I_LIB=INDEX_QP_TOOLS(1)
C             A new reduction tool will be used. Reinitialize the FPE
C              flags.
              FPE_IN_DP_REDUCTION=.FALSE.
              FPE_IN_QP_REDUCTION=.FALSE.
              I_QP_LIB=1
              DOING_QP_EVALS=.TRUE.
              EVAL_DONE(1)=.TRUE.
              DO I=2,MAXSTABILITYLENGTH
                EVAL_DONE(I)=.FALSE.
              ENDDO
              STAB_INDEX=0
              CTMODE=4
              GOTO 200
            ELSEIF(I_LIB.LE.NLOOPLIB.AND.MLREDUCTIONLIB(I_LIB).GT.0)
     $       THEN
C             A new reduction tool will be used. Reinitialize the FPE
C              flags.
              FPE_IN_DP_REDUCTION=.FALSE.
              FPE_IN_QP_REDUCTION=.FALSE.
              EVAL_DONE(1)=.TRUE.
              DO I=2,MAXSTABILITYLENGTH
                EVAL_DONE(I)=.FALSE.
              ENDDO
              STAB_INDEX=0
              IF(NROTATIONS_DP.GE.1)THEN
                GOTO 200
              ELSE
                GOTO 300
              ENDIF
            ELSE
              DO I=0,NSQUAREDSO
                ACCURACY(I)=ACC(I)
              ENDDO
              RET_CODE_H=4
              RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE.
     $         ,.FALSE.)
              NEPS=NEPS+1
              IF(NEPS.LE.10) THEN
                WRITE(*,*) '##W03 WARNING An unstable PS point was',  
     $                ' detected.'
                WRITE(*,*) '##W03 WARNING No quadruple precision will'
     $           //' be used.'
                IF(FPE_IN_DP_REDUCTION) THEN
                  WRITE(*,*) '## The last DP reduction was deemed'
     $             //' unstable because a floating point exception was'
     $             //' encountered.'
                ENDIF
                CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,ACC,ANS)
                IF (NSQUAREDSO.NE.1) THEN
                  WRITE(*,*) 'Accuracies for each split order,'
     $             //' starting with the summed case'
                  WRITE(*,*) 'DP accuracies (for each split order): '
     $             ,(ACC(I),I=0,NSQUAREDSO)
                ELSE
                  WRITE(*,*) 'DP accuracy: ',ACC(1)
                ENDIF
                DO J=0,NSQUAREDSO
                  IF (NSQUAREDSO.NE.1.OR.J.NE.0) THEN
                    IF (J.EQ.0) THEN
                      WRITE(*,*) 'Details for all split orders summed'
     $                 //' :'
                    ELSE
                      WRITE(*,*) 'Details for split order index : ',J
                    ENDIF
                    WRITE(*,*) 'Best estimate (fin,1eps,2eps):',(ANS(I
     $               ,J),I=1,3)
                    WRITE(*,*) 'Finite double precision evaluations :'
     $               ,(DP_RES(1,J,I),I=1,N_DP_EVAL)
                  ENDIF
                ENDDO
                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
            DO I=0,NSQUAREDSO
              ACCURACY(I)=ACC(I)
            ENDDO
            RET_CODE_H=2
            RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE.
     $       ,.TRUE.)
          ENDIF
        ENDIF
      ELSE
        RET_CODE_H=1
        DO I=0,NSQUAREDSO
          ACCURACY(I)=-1.0D0
        ENDDO
        RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE.
     $   ,.FALSE.)
      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_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE.
     $   ,.FALSE.)
        RET_CODE_T=RET_CODE_T+2
        DO I=0,NSQUAREDSO
          ACCURACY(I)=-1.0D0
        ENDDO
      ENDIF

C     Finally for the summed result in ANS(1:3,0), make sure to only
C     consider the squared order asked for by the user.
C     Notice that this filtering using CHOSEN_SO_CONFIGS happens
C     here only while everywhere else one always considers the sum.
      DO J=1,3
        ANS(J,0)=0.0D0
      ENDDO
      DO I=1,NSQUAREDSO
        IF (CHOSEN_SO_CONFIGS(I)) THEN
          DO J=1,3
            ANS(J,0)=ANS(J,0)+ANS(J,I)
          ENDDO
        ENDIF
      ENDDO

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     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,NSQUAREDSO
          DO J=1,NLOOPGROUPS
            GOODAMP(I,J)=OLD_GOODAMP(I,J)
          ENDDO
        ENDDO
      ENDIF

C     Make sure that we finish by emptying caches
      IF (AUTOMATIC_CACHE_CLEARING) THEN
        CALL ML5_0_CLEAR_CACHES()
      ENDIF

C     Now make sure to turn off the global COLLIER cache if applicable
      CALL ML5_0_SET_COLLIER_GLOBAL_CACHE(.FALSE.)

      END

      SUBROUTINE ML5_0_CLEAR_CACHES()
C     
C     This routine can be called directly from the user if
C     AUTOMATIC_CACHE_CLEARING is set to False. It must then be called
C      after
C     ech event
C     
      CALL ML5_0_CLEAR_TIR_CACHE()
      CALL NINJA_CLEAR_INTEGRAL_CACHE()
      CALL ML5_0_CLEAR_COLLIER_CACHE()
      END

C     --=========================================--
C     General Helper functions and subroutine
C     for the main sloopmatrix subroutine
C     --=========================================--

      LOGICAL FUNCTION ML5_0_IS_HEL_SELECTED(HELID)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NCOMB
      PARAMETER (NCOMB=16)
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, LOOP,
     $  SOINDEX)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NLOOPGROUPS
      PARAMETER (NLOOPGROUPS=120)
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=4)
C     
C     ARGUMENTS 
C     
      REAL*8 TOTEST, REFERENCE_VALUE
      INTEGER LOOP, SOINDEX
C     
C     GLOBAL 
C     
      INCLUDE 'MadLoopParams.inc'
      COMPLEX*16 LOOPRES(3,NSQUAREDSO,NLOOPGROUPS)
      LOGICAL S(NSQUAREDSO,NLOOPGROUPS)
      COMMON/ML5_0_LOOPRES/LOOPRES,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 1
      ELSE
        ML5_0_ISZERO=((ABS(TOTEST)/ABS(REFERENCE_VALUE)).LT.ZEROTHRES)
      ENDIF

      IF(LOOP.NE.-1) THEN
        IF((.NOT.ML5_0_ISZERO).AND.(.NOT.S(SOINDEX,LOOP))) THEN
          WRITE(*,*) '##W01 WARNING Contribution ',LOOP,' of split'
     $     //' order ',SOINDEX,' is detected as contributing with CR='
     $     ,(ABS(TOTEST)/ABS(REFERENCE_VALUE)),' but is unstable.'
        ENDIF
      ENDIF

      END

      INTEGER FUNCTION ML5_0_ISSAME(RESA,RESB,REF,USEMAX)
      IMPLICIT NONE
C     This function compares the result from two different helicity
C      configuration A and B
C     It returns 0 if they are not related and (+/-wgt) if
C      A=(+/-wgt)*B.
C     For now, the only wgt implemented is the integer 1 or -1.
C     If useMax is .TRUE., it uses all implemented weights no matter
C      what is HELINITSTARTOVER
C     
C     CONSTANTS
C     
      INTEGER MAX_WGT_TO_TRY
      PARAMETER (MAX_WGT_TO_TRY=2)
C     
C     ARGUMENTS 
C     
      REAL*8 RESA(3), RESB(3)
      REAL*8 REF
      LOGICAL USEMAX
C     
C     LOCAL VARIABLES
C     
      LOGICAL ML5_0_ISZERO
      INTEGER I,J
      INTEGER N_WGT_TO_TRY
      INTEGER WGT_TO_TRY(MAX_WGT_TO_TRY)
      DATA WGT_TO_TRY/1,-1/
C     
C     INCLUDES
C     
      INCLUDE 'MadLoopParams.inc'
C     ----------
C     BEGIN CODE
C     ----------
      ML5_0_ISSAME=0

C     If the helicity can be constructed progressively while allowing
C      inconsistency, then we only allow for weight one comparisons.
      IF (.NOT.HELINITSTARTOVER.AND..NOT.USEMAX) THEN
        N_WGT_TO_TRY=1
      ELSE
        N_WGT_TO_TRY=MAX_WGT_TO_TRY
      ENDIF

      DO I=1,N_WGT_TO_TRY
        DO J=1,3
          IF (ML5_0_ISZERO(ABS(RESB(J)),REF,-1,-1)) THEN
            IF(.NOT.ML5_0_ISZERO(ABS(RESB(J))+ABS(RESA(J)),REF,-1,-1))
     $        THEN
              GOTO 1231
            ENDIF
C           Be looser for helicity comparison, so bring a factor 100
          ELSEIF(.NOT.ML5_0_ISZERO(ABS((RESA(J)/RESB(J))
     $     -DBLE(WGT_TO_TRY(I))),1.0D0,-1,-1)) THEN
            GOTO 1231
          ENDIF
        ENDDO
        ML5_0_ISSAME = WGT_TO_TRY(I)
        RETURN
 1231   CONTINUE
      ENDDO
      END

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

C     
C     GLOBAL VARIABLES
C     
      LOGICAL CHOSEN_SO_CONFIGS(NSQUAREDSO)
      COMMON/ML5_0_CHOSEN_LOOP_SQSO/CHOSEN_SO_CONFIGS
      INTEGER I_LIB
      COMMON/ML5_0_I_LIB/I_LIB
      INCLUDE 'MadLoopParams.inc'

C     ----------
C     BEGIN CODE
C     ----------
      DO I=1,LENGTH
        MASK(I)=.TRUE.
      ENDDO
      DO I=LENGTH+1,MAXSTABILITYLENGTH
        MASK(I)=.FALSE.
C       For some architectures, it is necessary to initialize all the
C        elements of fulllist(i,j)
C       Beware that if the length provided is incorrect, then this can
C        corrup the fulllist given in argument.
        DO J=0,NSQUAREDSO
          DO K=1,3
            FULLLIST(K,J,I)=0.0D0
          ENDDO
        ENDDO
      ENDDO

      DO K=0,NSQUAREDSO

        DO I=1,3
          DO J=1,MAXSTABILITYLENGTH
            LIST(J)=FULLLIST(I,K,J)
          ENDDO
          DIFF=MAXVAL(LIST,1,MASK)-MINVAL(LIST,1,MASK)
          AVG=(MAXVAL(LIST,1,MASK)+MINVAL(LIST,1,MASK))/2.0D0
          ESTIMATE(I,K)=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(k)=MAXVAL(ACCURACIES,1,MASK3)
C       The following is used instead
        ACC(K) = 0.0D0
        AVG = 0.0D0
        DO I=1,3
          ACC(K) = ACC(K) + ACCURACIES(I)*ABS(ESTIMATE(I,K))
          AVG = AVG + ESTIMATE(I,K)
        ENDDO
        IF (AVG.NE.0.0D0) THEN
          ACC(K)  = ACC(K) / ( ABS(AVG) / 3.0D0)
        ENDIF

C       When using COLLIER with the internal stability test, the first
C        evaluation is typically more reliable so we do not want to
C        use the average but rather the first evaluation.
        IF (MLREDUCTIONLIB(I_LIB)
     $   .EQ.7.AND.COLLIERUSEINTERNALSTABILITYTEST) THEN
          DO I=1,3
            ESTIMATE(I,K) = FULLLIST(I,K,1)
          ENDDO
        ENDIF

C       Make sure to hard-set to zero accuracies of coupling orders
C        not included
        IF (K.NE.0) THEN
          IF (.NOT.CHOSEN_SO_CONFIGS(K)) THEN
            ACC(K) = 0.0D0
          ENDIF
        ENDIF

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,K,J))) THEN
              ACC(K) = 1.0D99
            ENDIF
          ENDDO
        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 1
      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=4)
      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

C     --=========================================--
C     Functions for dealing with the ordering
C     and indexing of split order contributions
C     --=========================================--

      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     
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=4)

      INTEGER NSQSO

      NSQSO=NSQUAREDSO

      END

      SUBROUTINE ML5_0_GET_ANSWER_DIMENSION(ANS_DIM)
C     
C     MadLoop subroutines return an array of dimension
C      ANS(0:3,0:ANS_DIM)
C     In order for the user program to be able to correctly declare
C      this
C     array when calling MadLoop, this subroutine returns its dimension
C     
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=4)
      INTEGER ANS_DIM

      INCLUDE 'nsqso_born.inc'

      ANS_DIM=MAX(NSQSO_BORN,NSQUAREDSO)

      END

      INTEGER FUNCTION ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS(ORDERS)
C     
C     This functions returns the integer index identifying the split
C      orders list passed in argument which correspond to the values
C      of the following list of couplings (and in this order):
C     ['QCD', 'QED']
C     
C     CONSTANTS
C     
      INTEGER    NSO, NSQSO
      PARAMETER (NSO=2, NSQSO=4)
C     
C     ARGUMENTS
C     
      INTEGER ORDERS(NSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I,J
      INTEGER SQPLITORDERS(NSQSO,NSO)
      DATA (SQPLITORDERS(  1,I),I=  1,  2) /    6,    0/
      DATA (SQPLITORDERS(  2,I),I=  1,  2) /    4,    2/
      DATA (SQPLITORDERS(  3,I),I=  1,  2) /    2,    4/
      DATA (SQPLITORDERS(  4,I),I=  1,  2) /    0,    6/
      COMMON/ML5_0_ML5SQPLITORDERS/SQPLITORDERS
C     
C     BEGIN CODE
C     
      DO I=1,NSQSO
        DO J=1,NSO
          IF (ORDERS(J).NE.SQPLITORDERS(I,J)) GOTO 1009
        ENDDO
        ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS = I
        RETURN
 1009   CONTINUE
      ENDDO

      WRITE(*,*) 'ERROR:: Stopping function'
     $ //' ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS'
      WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
      STOP

      END

      INTEGER FUNCTION ML5_0_GETORDPOWFROMINDEX_ML5(IORDER, INDX)
C     
C     Return the power of the IORDER-th order appearing at position
C      INDX
C     in the split-orders output
C     
C     ['QCD', 'QED']
C     
C     CONSTANTS
C     
      INTEGER    NSO, NSQSO
      PARAMETER (NSO=2, NSQSO=4)
C     
C     ARGUMENTS
C     
      INTEGER ORDERS(NSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I,J
      INTEGER SQPLITORDERS(NSQSO,NSO)
      DATA (SQPLITORDERS(  1,I),I=  1,  2) /    6,    0/
      DATA (SQPLITORDERS(  2,I),I=  1,  2) /    4,    2/
      DATA (SQPLITORDERS(  3,I),I=  1,  2) /    2,    4/
      DATA (SQPLITORDERS(  4,I),I=  1,  2) /    0,    6/
C     
C     BEGIN CODE
C     
      IF (IORDER.GT.NSO.OR.IORDER.LT.1) THEN
        WRITE(*,*) 'INVALID IORDER ML5', IORDER
        WRITE(*,*) 'SHOULD BE BETWEEN 1 AND ', NSO
        STOP
      ENDIF

      IF (INDX.GT.NSQSO.OR.INDX.LT.1) THEN
        WRITE(*,*) 'INVALID INDX ML5', INDX
        WRITE(*,*) 'SHOULD BE BETWEEN 1 AND ', NSQSO
        STOP
      ENDIF

      ML5_0_GETORDPOWFROMINDEX_ML5=SQPLITORDERS(INDX, IORDER)

      END

      INTEGER FUNCTION ML5_0_ML5SOINDEX_FOR_BORN_AMP(AMPID)
C     
C     For a given born amplitude number, it returns the ID of the
C      split orders it has
C     
C     CONSTANTS
C     
      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=7)
C     
C     ARGUMENTS
C     
      INTEGER AMPID
C     
C     LOCAL VARIABLES
C     
      INTEGER BORNAMPORDERS(NBORNAMPS)
      DATA (BORNAMPORDERS(I),I=  1,  5) /    1,    2,    2,    2,    2/
      DATA (BORNAMPORDERS(I),I=  6,  7) /    2,    2/
C     -----------
C     BEGIN CODE
C     -----------
      IF (AMPID.GT.NBORNAMPS) THEN
        WRITE(*,*) 'ERROR:: Born amplitude ID ',AMPID,' above the'
     $   //' maximum ',NBORNAMPS
      ENDIF
      ML5_0_ML5SOINDEX_FOR_BORN_AMP = BORNAMPORDERS(AMPID)

      END

      INTEGER FUNCTION ML5_0_ML5SOINDEX_FOR_LOOP_AMP(AMPID)
C     
C     For a given loop amplitude number, it returns the ID of the
C      split orders it has
C     
C     CONSTANTS
C     
      INTEGER    NLOOPAMPS
      PARAMETER (NLOOPAMPS=1022)
C     
C     ARGUMENTS
C     
      INTEGER AMPID
C     
C     LOCAL VARIABLES
C     
      INTEGER LOOPAMPORDERS(NLOOPAMPS)
      DATA (LOOPAMPORDERS(I),I=  1,  5) /    3,    3,    3,    3,    3/
      DATA (LOOPAMPORDERS(I),I=  6, 10) /    3,    3,    3,    3,    3/
      DATA (LOOPAMPORDERS(I),I= 11, 15) /    3,    3,    3,    3,    3/
      DATA (LOOPAMPORDERS(I),I= 16, 20) /    3,    3,    3,    3,    3/
      DATA (LOOPAMPORDERS(I),I= 21, 25) /    3,    3,    3,    3,    3/
      DATA (LOOPAMPORDERS(I),I= 26, 30) /    3,    3,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I= 31, 35) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I= 36, 40) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I= 41, 45) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I= 46, 50) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I= 51, 55) /    4,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I= 56, 60) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I= 61, 65) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I= 66, 70) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I= 71, 75) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I= 76, 80) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I= 81, 85) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I= 86, 90) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I= 91, 95) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I= 96,100) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=101,105) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=106,110) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=111,115) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=116,120) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=121,125) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=126,130) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=131,135) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=136,140) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=141,145) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=146,150) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=151,155) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=156,160) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=161,165) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=166,170) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=171,175) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=176,180) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=181,185) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=186,190) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=191,195) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=196,200) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=201,205) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=206,210) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=211,215) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=216,220) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=221,225) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=226,230) /    5,    5,    5,    3,    4/
      DATA (LOOPAMPORDERS(I),I=231,235) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=236,240) /    3,    3,    3,    3,    3/
      DATA (LOOPAMPORDERS(I),I=241,245) /    3,    3,    3,    3,    3/
      DATA (LOOPAMPORDERS(I),I=246,250) /    3,    3,    3,    3,    4/
      DATA (LOOPAMPORDERS(I),I=251,255) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=256,260) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=261,265) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=266,270) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=271,275) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=276,280) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=281,285) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=286,290) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=291,295) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=296,300) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=301,305) /    4,    4,    4,    4,    4/
      DATA (LOOPAMPORDERS(I),I=306,310) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=311,315) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=316,320) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=321,325) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=326,330) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=331,335) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=336,340) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=341,345) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=346,350) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=351,355) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=356,360) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=361,365) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=366,370) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=371,375) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=376,380) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=381,385) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=386,390) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=391,395) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=396,400) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=401,405) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=406,410) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=411,415) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=416,420) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=421,425) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=426,430) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=431,435) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=436,440) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=441,445) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=446,450) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=451,455) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=456,460) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=461,465) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=466,470) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=471,475) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=476,480) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=481,485) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=486,490) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=491,495) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=496,500) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=501,505) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=506,510) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=511,515) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=516,520) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=521,525) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=526,530) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=531,535) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=536,540) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=541,545) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=546,550) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=551,555) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=556,560) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=561,565) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=566,570) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=571,575) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=576,580) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=581,585) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=586,590) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=591,595) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=596,600) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=601,605) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=606,610) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=611,615) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=616,620) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=621,625) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=626,630) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=631,635) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=636,640) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=641,645) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=646,650) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=651,655) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=656,660) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=661,665) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=666,670) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=671,675) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=676,680) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=681,685) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=686,690) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=691,695) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=696,700) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=701,705) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=706,710) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=711,715) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=716,720) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=721,725) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=726,730) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=731,735) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=736,740) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=741,745) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=746,750) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=751,755) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=756,760) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=761,765) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=766,770) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=771,775) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=776,780) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=781,785) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=786,790) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=791,795) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=796,800) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=801,805) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=806,810) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=811,815) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=816,820) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=821,825) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=826,830) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=831,835) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=836,840) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=841,845) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=846,850) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=851,855) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=856,860) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=861,865) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=866,870) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=871,875) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=876,880) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=881,885) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=886,890) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=891,895) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=896,900) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=901,905) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=906,910) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=911,915) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=916,920) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=921,925) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=926,930) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=931,935) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=936,940) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=941,945) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=946,950) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=951,955) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=956,960) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=961,965) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=966,970) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=971,975) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=976,980) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=981,985) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=986,990) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=991,995) /    5,    5,    5,    5,    5/
      DATA (LOOPAMPORDERS(I),I=996,1000) /    5,    5,    5,    5,   
     $  5/
      DATA (LOOPAMPORDERS(I),I=1001,1005) /    5,    5,    5,    5,   
     $  5/
      DATA (LOOPAMPORDERS(I),I=1006,1010) /    5,    5,    5,    5,   
     $  5/
      DATA (LOOPAMPORDERS(I),I=1011,1015) /    5,    5,    5,    5,   
     $  5/
      DATA (LOOPAMPORDERS(I),I=1016,1020) /    5,    5,    5,    5,   
     $  5/
      DATA (LOOPAMPORDERS(I),I=1021,1022) /    5,    5/
C     -----------
C     BEGIN CODE
C     -----------
      IF (AMPID.GT.NLOOPAMPS) THEN
        WRITE(*,*) 'ERROR:: Loop amplitude ID ',AMPID,' above the'
     $   //' maximum ',NLOOPAMPS
      ENDIF
      ML5_0_ML5SOINDEX_FOR_LOOP_AMP = LOOPAMPORDERS(AMPID)

      END


      INTEGER FUNCTION ML5_0_ML5SQSOINDEX(ORDERINDEXA, ORDERINDEXB)
C     
C     This functions plays the role of the interference matrix. It can
C      be hardcoded or 
C     made more elegant using hashtables if its execution speed ever
C      becomes a relevant
C     factor. From two split order indices, it return the
C      corresponding index in the squared 
C     order canonical ordering.
C     
C     CONSTANTS
C     
      INTEGER    NSO, NSQUAREDSO, NAMPSO
      PARAMETER (NSO=2, NSQUAREDSO=4, NAMPSO=5)
C     
C     ARGUMENTS
C     
      INTEGER ORDERINDEXA, ORDERINDEXB
C     
C     LOCAL VARIABLES
C     
      INTEGER I, SQORDERS(NSO)
      INTEGER AMPSPLITORDERS(NAMPSO,NSO)
      DATA (AMPSPLITORDERS(  1,I),I=  1,  2) /    2,    0/
      DATA (AMPSPLITORDERS(  2,I),I=  1,  2) /    0,    2/
      DATA (AMPSPLITORDERS(  3,I),I=  1,  2) /    4,    0/
      DATA (AMPSPLITORDERS(  4,I),I=  1,  2) /    2,    2/
      DATA (AMPSPLITORDERS(  5,I),I=  1,  2) /    0,    4/
      COMMON/ML5_0_ML5AMPSPLITORDERS/AMPSPLITORDERS
C     
C     FUNCTION
C     
      INTEGER ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS
C     
C     BEGIN CODE
C     
      DO I=1,NSO
        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)
     $   +AMPSPLITORDERS(ORDERINDEXB,I)
      ENDDO
      ML5_0_ML5SQSOINDEX=ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
      END

C     This is the inverse subroutine of ML5SOINDEX_FOR_SQUARED_ORDERS.
C      Not directly useful, but provided nonetheless.
      SUBROUTINE ML5_0_ML5GET_SQUARED_ORDERS_FOR_SOINDEX(SOINDEX
     $ ,ORDERS)
C     
C     This functions returns the orders identified by the squared
C      split order index in argument. Order values correspond to
C      following list of couplings (and in this order):
C     ['QCD', 'QED']
C     
C     CONSTANTS
C     
      INTEGER    NSO, NSQSO
      PARAMETER (NSO=2, NSQSO=4)
C     
C     ARGUMENTS
C     
      INTEGER SOINDEX, ORDERS(NSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I
      INTEGER SQPLITORDERS(NSQSO,NSO)
      COMMON/ML5_0_ML5SQPLITORDERS/SQPLITORDERS
C     
C     BEGIN CODE
C     
      IF (SOINDEX.GT.0.AND.SOINDEX.LE.NSQSO) THEN
        DO I=1,NSO
          ORDERS(I) =  SQPLITORDERS(SOINDEX,I)
        ENDDO
        RETURN
      ENDIF

      WRITE(*,*) 'ERROR:: Stopping function'
     $ //' ML5_0_ML5GET_SQUARED_ORDERS_FOR_SOINDEX'
      WRITE(*,*) 'Could not find squared orders index ',SOINDEX
      STOP

      END SUBROUTINE

C     This is the inverse subroutine of getting amplitude SO orders.
C      Not directly useful, but provided nonetheless.
      SUBROUTINE ML5_0_ML5GET_ORDERS_FOR_AMPSOINDEX(SOINDEX,ORDERS)
C     
C     This functions returns the orders identified by the split order
C      index in argument. Order values correspond to following list of
C      couplings (and in this order):
C     ['QCD', 'QED']
C     
C     CONSTANTS
C     
      INTEGER    NSO, NAMPSO
      PARAMETER (NSO=2, NAMPSO=5)
C     
C     ARGUMENTS
C     
      INTEGER SOINDEX, ORDERS(NSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I
      INTEGER AMPSPLITORDERS(NAMPSO,NSO)
      COMMON/ML5_0_ML5AMPSPLITORDERS/AMPSPLITORDERS
C     
C     BEGIN CODE
C     
      IF (SOINDEX.GT.0.AND.SOINDEX.LE.NAMPSO) THEN
        DO I=1,NSO
          ORDERS(I) =  AMPSPLITORDERS(SOINDEX,I)
        ENDDO
        RETURN
      ENDIF

      WRITE(*,*) 'ERROR:: Stopping function'
     $ //' ML5_0_ML5GET_ORDERS_FOR_AMPSOINDEX'
      WRITE(*,*) 'Could not find amplitude split orders index ',SOINDEX
      STOP

      END SUBROUTINE


C     This function is not directly useful, but included for
C      completeness
      INTEGER FUNCTION ML5_0_ML5SOINDEX_FOR_AMPORDERS(ORDERS)
C     
C     This functions returns the integer index identifying the
C      amplitude split orders passed in argument which correspond to
C      the values of the following list of couplings (and in this
C      order):
C     ['QCD', 'QED']
C     
C     CONSTANTS
C     
      INTEGER    NSO, NAMPSO
      PARAMETER (NSO=2, NAMPSO=5)
C     
C     ARGUMENTS
C     
      INTEGER ORDERS(NSO)
C     
C     LOCAL VARIABLES
C     
      INTEGER I,J
      INTEGER AMPSPLITORDERS(NAMPSO,NSO)
      COMMON/ML5_0_ML5AMPSPLITORDERS/AMPSPLITORDERS
C     
C     BEGIN CODE
C     
      DO I=1,NAMPSO
        DO J=1,NSO
          IF (ORDERS(J).NE.AMPSPLITORDERS(I,J)) GOTO 1009
        ENDDO
        ML5_0_ML5SOINDEX_FOR_AMPORDERS = I
        RETURN
 1009   CONTINUE
      ENDDO

      WRITE(*,*) 'ERROR:: Stopping function'
     $ //' ML5_0_ML5SOINDEX_FOR_AMPORDERS'
      WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
      STOP

      END

C     --=========================================--
C     Definition of additional access routines
C     --=========================================--

      SUBROUTINE ML5_0_COLLIER_COMPUTE_UV_POLES(ONOFF)
C     
C     This function can be called by the MadLoop user so as to chose
C      to have COLLIER
C     compute the UV pole or not (it costs more time).
C     
      LOGICAL ONOFF

      INCLUDE 'MadLoopParams.inc'

      LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
      COMMON/ML5_0_COLLIERPOLESFORCEDCHOICE
     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE

      COLLIERCOMPUTEUVPOLES                        = ONOFF
C     This is just so that if we read the param again, we don't
C      overwrite the choice made here
      FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION = .TRUE.
      COLLIER_UV_POLE_COMPUTATION_CHOICE           = ONOFF

      END SUBROUTINE

      SUBROUTINE ML5_0_COLLIER_COMPUTE_IR_POLES(ONOFF)
C     
C     This function can be called by the MadLoop user so as to chose
C      to have COLLIER
C     compute the IR pole or not (it costs more time).
C     
      LOGICAL ONOFF

      INCLUDE 'MadLoopParams.inc'

      LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
      COMMON/ML5_0_COLLIERPOLESFORCEDCHOICE
     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE

      COLLIERCOMPUTEIRPOLES         = ONOFF
C     This is just so that if we read the param again, we don't
C      overwrite the choice made here
      FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION = .TRUE.
      COLLIER_IR_POLE_COMPUTATION_CHOICE           = ONOFF

      END SUBROUTINE

      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

      SUBROUTINE ML5_0_SET_AUTOMATIC_CACHE_CLEARING(ONOFF)
C     
C     This function can be called by the MadLoop user so as to
C      manually chose when
C     to reset the TIR cache.
C     
      IMPLICIT NONE

      INCLUDE 'MadLoopParams.inc'

      LOGICAL ONOFF

      LOGICAL AUTOMATIC_CACHE_CLEARING
      DATA AUTOMATIC_CACHE_CLEARING/.TRUE./
      COMMON/ML5_0_RUNTIME_OPTIONS/AUTOMATIC_CACHE_CLEARING

      INTEGER N_DP_EVAL, N_QP_EVAL
      COMMON/ML5_0_N_EVALS/N_DP_EVAL,N_QP_EVAL


      AUTOMATIC_CACHE_CLEARING = ONOFF

      IF (NROTATIONS_DP.NE.0.OR.NROTATIONS_QP.NE.0) THEN
        WRITE(*,*) 'Warning: One cannot remove the TIR cache automatic'
     $   //' clearing while at the same time keeping Lorentz rotations'
     $   //' for stability tests.'
        WRITE(*,*) 'MadLoop will therefore automatically set'
     $   //' NRotations_DP and NRotations_QP to 0.'
        NROTATIONS_DP = 0
        NROTATIONS_QP = 0
        CALL ML5_0_SET_N_EVALS(N_DP_EVAL,N_QP_EVAL)
      ENDIF
      END SUBROUTINE

      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     If set to a value different than -1, the code will try to avoid
C      computing anything which
C     does not contribute to contributions of squared split orders
C      SQSO_TARGET and below.
C     This can considerably speed up the code. However, keep in mind
C      that any contribution of
C     'squared order index' larger than SQSO_TARGET cannot be trust.
C     
C     ARGUMENTS
C     
      INTEGER SOTARGET
C     
C     GLOBAL
C     
      INTEGER SQSO_TARGET
      COMMON/ML5_0_SOCHOICE/SQSO_TARGET
C     ----------
C     BEGIN CODE
C     ----------
      SQSO_TARGET = SOTARGET
      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=4)
      INTEGER NPOLENTRIES
      PARAMETER (NPOLENTRIES=(NEXTERNAL+1)*6)
      INTEGER    NCOMB
      PARAMETER (NCOMB=16)
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(P,HEL,ANS)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=4)
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

      SUBROUTINE ML5_0_SLOOPMATRIXHEL_THRES(P,HEL,ANS,PREC_ASKED
     $ ,PREC_FOUND,RET_CODE)
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=4)
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     LOCAL VARIABLES
C     
      INTEGER I
C     
C     GLOBAL VARIABLES
C     
      REAL*8 USER_STAB_PREC
      COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC

      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 == 0
C     Not stable.
C     U == 1
C     Stable with CutTools in double precision.
C     U == 2
C     Stable with PJFry++.
C     U == 3
C     Stable with IREGI.
C     U == 4
C     Stable with Golem95
C     U == 5
C     Stable with Samurai
C     U == 6
C     Stable with Ninja in double precision
C     U == 8
C     Stable with Ninja in quadruple precision
C     U == 9
C     Stable with CutTools in quadruple precision.         
C     
      IMPLICIT NONE
C     
C     CONSTANTS
C     
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NSQUAREDSO
      PARAMETER (NSQUAREDSO=4)
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     LOCAL VARIABLES
C     
      INTEGER I
C     
C     GLOBAL VARIABLES
C     
      REAL*8 USER_STAB_PREC
      COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC

      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

C     The subroutine below perform clean-up duties for MadLoop like
C      de-allocating
C     arrays
      SUBROUTINE ML5_0_EXIT_MADLOOP()
      CONTINUE
      END