~maddevelopers/mg5amcnlo/2.5.3

« back to all changes in this revision

Viewing changes to tests/input_files/IOTestsComparison/long_ML_SMQCD_default/dux_mumvmxg/loop_matrix.f

  • Committer: olivier-mattelaer
  • Date: 2017-03-08 12:31:17 UTC
  • Revision ID: olivier-mattelaer-20170308123117-h0zkqjyh9sihsc61
empty version to have an effective freeze of the code

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
      SUBROUTINE ML5_0_SLOOPMATRIXHEL(P,HEL,ANS)
2
 
      IMPLICIT NONE
3
 
C     
4
 
C     CONSTANTS
5
 
C     
6
 
      INTEGER    NEXTERNAL
7
 
      PARAMETER (NEXTERNAL=5)
8
 
 
9
 
      INCLUDE 'nsquaredSO.inc'
10
 
 
11
 
C     
12
 
C     ARGUMENTS 
13
 
C     
14
 
      REAL*8 P(0:3,NEXTERNAL)
15
 
      REAL*8 ANS(0:3,0:NSQUAREDSO)
16
 
      INTEGER HEL, USERHEL
17
 
      COMMON/ML5_0_USERCHOICE/USERHEL
18
 
C     ----------
19
 
C     BEGIN CODE
20
 
C     ----------
21
 
      USERHEL=HEL
22
 
      CALL ML5_0_SLOOPMATRIX(P,ANS)
23
 
      END
24
 
 
25
 
      LOGICAL FUNCTION ML5_0_IS_HEL_SELECTED(HELID)
26
 
      IMPLICIT NONE
27
 
C     
28
 
C     CONSTANTS
29
 
C     
30
 
      INTEGER    NEXTERNAL
31
 
      PARAMETER (NEXTERNAL=5)
32
 
      INTEGER    NCOMB
33
 
      PARAMETER (NCOMB=32)
34
 
C     
35
 
C     ARGUMENTS
36
 
C     
37
 
      INTEGER HELID
38
 
C     
39
 
C     LOCALS
40
 
C     
41
 
      INTEGER I,J
42
 
      LOGICAL FOUNDIT
43
 
C     
44
 
C     GLOBALS
45
 
C     
46
 
      INTEGER HELC(NEXTERNAL,NCOMB)
47
 
      COMMON/ML5_0_HELCONFIGS/HELC
48
 
 
49
 
      INTEGER POLARIZATIONS(0:NEXTERNAL,0:5)
50
 
      COMMON/ML5_0_BEAM_POL/POLARIZATIONS
51
 
C     ----------
52
 
C     BEGIN CODE
53
 
C     ----------
54
 
 
55
 
      ML5_0_IS_HEL_SELECTED = .TRUE.
56
 
      IF (POLARIZATIONS(0,0).EQ.-1) THEN
57
 
        RETURN
58
 
      ENDIF
59
 
 
60
 
      DO I=1,NEXTERNAL
61
 
        IF (POLARIZATIONS(I,0).EQ.-1) THEN
62
 
          CYCLE
63
 
        ENDIF
64
 
        FOUNDIT = .FALSE.
65
 
        DO J=1,POLARIZATIONS(I,0)
66
 
          IF (HELC(I,HELID).EQ.POLARIZATIONS(I,J)) THEN
67
 
            FOUNDIT = .TRUE.
68
 
            EXIT
69
 
          ENDIF
70
 
        ENDDO
71
 
        IF(.NOT.FOUNDIT) THEN
72
 
          ML5_0_IS_HEL_SELECTED = .FALSE.
73
 
          RETURN
74
 
        ENDIF
75
 
      ENDDO
76
 
      RETURN
77
 
 
78
 
      END
79
 
 
80
 
      LOGICAL FUNCTION ML5_0_ISZERO(TOTEST, REFERENCE_VALUE, AMPLN)
81
 
      IMPLICIT NONE
82
 
C     
83
 
C     CONSTANTS
84
 
C     
85
 
      INTEGER    NLOOPAMPS
86
 
      PARAMETER (NLOOPAMPS=39)
87
 
C     
88
 
C     ARGUMENTS 
89
 
C     
90
 
      REAL*8 TOTEST, REFERENCE_VALUE
91
 
      INTEGER AMPLN
92
 
C     
93
 
C     GLOBAL 
94
 
C     
95
 
      INCLUDE 'MadLoopParams.inc'
96
 
 
97
 
      COMPLEX*16 AMPL(3,NLOOPAMPS)
98
 
      LOGICAL S(NLOOPAMPS)
99
 
      COMMON/ML5_0_AMPL/AMPL,S
100
 
C     ----------
101
 
C     BEGIN CODE
102
 
C     ----------
103
 
      IF(ABS(REFERENCE_VALUE).EQ.0.0D0) THEN
104
 
        ML5_0_ISZERO=.FALSE.
105
 
        WRITE(*,*) '##E02 ERRROR Reference value for comparison is'
106
 
     $   //' zero.'
107
 
        STOP
108
 
      ELSE
109
 
        ML5_0_ISZERO=((ABS(TOTEST)/ABS(REFERENCE_VALUE)).LT.ZEROTHRES)
110
 
      ENDIF
111
 
      IF(AMPLN.NE.-1) THEN
112
 
        IF((.NOT.ML5_0_ISZERO).AND.(.NOT.S(AMPLN))) THEN
113
 
          WRITE(*,*) '##W01 WARNING Contribution ',AMPLN,' is detected'
114
 
     $     //' as contributing with CR=',(ABS(TOTEST)
115
 
     $     /ABS(REFERENCE_VALUE)),' but is unstable.'
116
 
        ENDIF
117
 
      ENDIF
118
 
 
119
 
      END
120
 
 
121
 
      SUBROUTINE ML5_0_SLOOPMATRIX(P_USER,ANSRETURNED)
122
 
C     
123
 
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
124
 
C     By the MadGraph5_aMC@NLO Development Team
125
 
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
126
 
C     
127
 
C     Returns amplitude squared summed/avg over colors
128
 
C     and helicities for the point in phase space P(0:3,NEXTERNAL)
129
 
C     and external lines W(0:6,NEXTERNAL)
130
 
C     
131
 
C     Process: d u~ > m- vm~ g QCD<=1 QED<=2 [ virt = QCD ]
132
 
C     
133
 
      IMPLICIT NONE
134
 
C     
135
 
C     CONSTANTS
136
 
C     
137
 
      CHARACTER*512 PARAMFNAME,HELCONFIGFNAME,LOOPFILTERFNAME
138
 
      CHARACTER*512 COLORNUMFNAME,COLORDENOMFNAME, HELFILTERFNAME
139
 
      CHARACTER*512 PROC_PREFIX
140
 
      PARAMETER ( PARAMFNAME='MadLoopParams.dat')
141
 
      PARAMETER ( HELCONFIGFNAME='HelConfigs.dat')
142
 
      PARAMETER ( LOOPFILTERFNAME='LoopFilter.dat')
143
 
      PARAMETER ( HELFILTERFNAME='HelFilter.dat')
144
 
      PARAMETER ( COLORNUMFNAME='ColorNumFactors.dat')
145
 
      PARAMETER ( COLORDENOMFNAME='ColorDenomFactors.dat')
146
 
      PARAMETER ( PROC_PREFIX='ML5_0_')
147
 
 
148
 
      INTEGER NBORNAMPS
149
 
      PARAMETER (NBORNAMPS=2)
150
 
      INTEGER    NLOOPAMPS, NCTAMPS
151
 
      PARAMETER (NLOOPAMPS=39, NCTAMPS=28)
152
 
      INTEGER    NEXTERNAL
153
 
      PARAMETER (NEXTERNAL=5)
154
 
      INTEGER NINITIAL
155
 
      PARAMETER (NINITIAL=2)
156
 
      INTEGER    NWAVEFUNCS
157
 
      PARAMETER (NWAVEFUNCS=10)
158
 
      INTEGER    NCOMB
159
 
      PARAMETER (NCOMB=32)
160
 
      REAL*8     ZERO
161
 
      PARAMETER (ZERO=0D0)
162
 
      REAL*16     MP__ZERO
163
 
      PARAMETER (MP__ZERO=0E0_16)
164
 
      COMPLEX*16 IMAG1
165
 
      PARAMETER (IMAG1=(0D0,1D0))
166
 
C     This parameter is designed for the check timing command of MG5
167
 
      LOGICAL SKIPLOOPEVAL
168
 
      PARAMETER (SKIPLOOPEVAL=.FALSE.)
169
 
      LOGICAL BOOTANDSTOP
170
 
      PARAMETER (BOOTANDSTOP=.FALSE.)
171
 
      INCLUDE 'nsquaredSO.inc'
172
 
      INTEGER NSQUAREDSOP1
173
 
      PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1)
174
 
      INTEGER MAXSTABILITYLENGTH
175
 
      DATA MAXSTABILITYLENGTH/20/
176
 
      COMMON/ML5_0_STABILITY_TESTS/MAXSTABILITYLENGTH
177
 
C     
178
 
C     ARGUMENTS 
179
 
C     
180
 
      REAL*8 P_USER(0:3,NEXTERNAL)
181
 
      REAL*8 ANSRETURNED(0:3,0:NSQUAREDSO)
182
 
C     
183
 
C     LOCAL VARIABLES 
184
 
C     
185
 
      REAL*8 ANS(0:3)
186
 
      INTEGER I,J,K,H
187
 
 
188
 
      CHARACTER*512 PARAMFN,HELCONFIGFN,LOOPFILTERFN,COLORNUMFN
189
 
     $ ,COLORDENOMFN,HELFILTERFN
190
 
      CHARACTER*512 TMP
191
 
      SAVE PARAMFN
192
 
      SAVE HELCONFIGFN
193
 
      SAVE LOOPFILTERFN
194
 
      SAVE COLORNUMFN
195
 
      SAVE COLORDENOMFN
196
 
      SAVE HELFILTERFN
197
 
 
198
 
      INTEGER HELPICKED_BU, CTMODEINIT_BU
199
 
      REAL*8 MLSTABTHRES_BU
200
 
C     P is the actual PS POINT used for the computation, and can be
201
 
C      rotated for the stability test purposes.
202
 
      REAL*8 P(0:3,NEXTERNAL)
203
 
C     DP_RES STORES THE DOUBLE PRECISION RESULT OBTAINED FROM
204
 
C      DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY.
205
 
C     THE STAB_STAGE COUNTER I CORRESPONDANCE GOES AS FOLLOWS
206
 
C     I=1 -> ORIGINAL PS, CTMODE=1
207
 
C     I=2 -> ORIGINAL PS, CTMODE=2, (ONLY WITH CTMODERUN=-1)
208
 
C     I=3 -> PS WITH ROTATION 1, CTMODE=1, (ONLY WITH CTMODERUN=-2)
209
 
C     I=4 -> PS WITH ROTATION 2, CTMODE=1, (ONLY WITH CTMODERUN=-3)
210
 
C     I=5 -> POSSIBLY MORE EVALUATION METHODS IN THE FUTURE, MAX IS
211
 
C      MAXSTABILITYLENGTH
212
 
C     IF UNSTABLE IT GOES TO THE SAME PATTERN BUT STAB_INDEX IS THEN
213
 
C      I+20.
214
 
      LOGICAL EVAL_DONE(MAXSTABILITYLENGTH)
215
 
      LOGICAL DOING_QP_EVALS
216
 
      INTEGER STAB_INDEX,BASIC_CT_MODE
217
 
      INTEGER N_DP_EVAL, N_QP_EVAL
218
 
      DATA N_DP_EVAL/1/
219
 
      DATA N_QP_EVAL/1/
220
 
C     This is used for loop-induced where the reference scale for
221
 
C      comparisons is infered from
222
 
C     the previous points
223
 
      REAL*8 NEXTREF
224
 
      DATA NEXTREF/ZERO/
225
 
      INTEGER NPSPOINTS
226
 
      DATA NPSPOINTS/0/
227
 
      LOGICAL FOUND_VALID_REDUCTION_METHOD
228
 
      DATA FOUND_VALID_REDUCTION_METHOD/.FALSE./
229
 
 
230
 
      REAL*8 ACC
231
 
      REAL*8 DP_RES(3,MAXSTABILITYLENGTH)
232
 
C     QP_RES STORES THE QUADRUPLE PRECISION RESULT OBTAINED FROM
233
 
C      DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY.
234
 
      REAL*8 QP_RES(3,MAXSTABILITYLENGTH)
235
 
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
236
 
      INTEGER NATTEMPTS
237
 
      DATA NATTEMPTS/0/
238
 
      DATA IC/NEXTERNAL*1/
239
 
      REAL*8 BUFFR(3),TEMP(3),TEMP1,TEMP2
240
 
      COMPLEX*16 CFTOT
241
 
      LOGICAL FOUNDHELFILTER,FOUNDLOOPFILTER
242
 
      DATA FOUNDHELFILTER/.TRUE./
243
 
      DATA FOUNDLOOPFILTER/.TRUE./
244
 
      INTEGER IDEN
245
 
      DATA IDEN/36/
246
 
      INTEGER HELAVGFACTOR
247
 
      DATA HELAVGFACTOR/4/
248
 
C     For a 1>N process, them BEAMTWO_HELAVGFACTOR would be set to 1.
249
 
      INTEGER BEAMS_HELAVGFACTOR(2)
250
 
      DATA (BEAMS_HELAVGFACTOR(I),I=1,2)/2,2/
251
 
      LOGICAL DONEHELDOUBLECHECK
252
 
      DATA DONEHELDOUBLECHECK/.FALSE./
253
 
      INTEGER NEPS
254
 
      DATA NEPS/0/
255
 
C     Below are variables to bypass the checkphase and insure
256
 
C      stability check to take place
257
 
      LOGICAL OLD_CHECKPHASE, OLD_HELDOUBLECHECKED
258
 
      LOGICAL OLD_GOODHEL(NCOMB)
259
 
      LOGICAL OLD_GOODAMP(NLOOPAMPS,NCOMB)
260
 
 
261
 
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
262
 
      COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
263
 
C     
264
 
C     FUNCTIONS
265
 
C     
266
 
      LOGICAL ML5_0_ISZERO
267
 
      LOGICAL ML5_0_IS_HEL_SELECTED
268
 
C     
269
 
C     GLOBAL VARIABLES
270
 
C     
271
 
      INCLUDE 'process_info.inc'
272
 
      INCLUDE 'coupl.inc'
273
 
      INCLUDE 'mp_coupl.inc'
274
 
      INCLUDE 'MadLoopParams.inc'
275
 
 
276
 
      INTEGER NTRY
277
 
      DATA NTRY/0/
278
 
      LOGICAL CHECKPHASE
279
 
      DATA CHECKPHASE/.TRUE./
280
 
      LOGICAL HELDOUBLECHECKED
281
 
      DATA HELDOUBLECHECKED/.FALSE./
282
 
      REAL*8 REF
283
 
      DATA REF/0.0D0/
284
 
      COMMON/ML5_0_INIT/NTRY,CHECKPHASE,HELDOUBLECHECKED,REF
285
 
 
286
 
C     THE LOGICAL BELOWS ARE JUST TO KEEP TRACK OF WHETHER THE MP_PS
287
 
C      HAS BEEN SET YET OR NOT AND WHETER THE MP EXTERNAL WFS HAVE
288
 
C      BEEN COMPUTED YET.
289
 
      LOGICAL MP_DONE
290
 
      DATA MP_DONE/.FALSE./
291
 
      COMMON/ML5_0_MP_DONE/MP_DONE
292
 
      LOGICAL MP_PS_SET
293
 
      DATA MP_PS_SET/.FALSE./
294
 
      COMMON/ML5_0_MP_PS_SET/MP_PS_SET
295
 
 
296
 
C     PS CAN POSSIBILY BE PASSED THROUGH IMPROVE_PS BUT IS NOT
297
 
C      MODIFIED FOR THE PURPOSE OF THE STABILITY TEST     
298
 
C     EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT
299
 
C      USED ANYWHERE ELSE
300
 
      REAL*8 PS(0:3,NEXTERNAL)
301
 
      COMMON/ML5_0_PSPOINT/PS
302
 
C     AGAIN BELOW, MP_PS IS THE FIXED (POSSIBLY IMPROVED) MP PS POINT
303
 
C      AND MP_P IS THE ONE WHICH CAN BE MODIFIED (I.E. ROTATED ETC.)
304
 
C      FOR STABILITY PURPOSE
305
 
C     EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT
306
 
C      USED ANYWHERE ELSE THAN HERE AND SET_MP_PS()
307
 
      REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL)
308
 
      COMMON/ML5_0_MP_PSPOINT/MP_PS,MP_P
309
 
 
310
 
      REAL*8 LSCALE
311
 
      INTEGER CTMODE
312
 
      COMMON/ML5_0_CT/LSCALE,CTMODE
313
 
 
314
 
      LOGICAL GOODHEL(NCOMB)
315
 
      LOGICAL GOODAMP(NLOOPAMPS,NCOMB)
316
 
      COMMON/ML5_0_FILTERS/GOODAMP,GOODHEL
317
 
 
318
 
      INTEGER HELPICKED
319
 
      DATA HELPICKED/-1/
320
 
      COMMON/ML5_0_HELCHOICE/HELPICKED
321
 
      INTEGER USERHEL
322
 
      DATA USERHEL/-1/
323
 
      COMMON/ML5_0_USERCHOICE/USERHEL
324
 
 
325
 
      COMPLEX*16 AMP(NBORNAMPS,NCOMB)
326
 
      COMMON/ML5_0_AMPS/AMP
327
 
      COMPLEX*16 W(20,NWAVEFUNCS,NCOMB)
328
 
      INTEGER VALIDH
329
 
      COMMON/ML5_0_WFCTS/W
330
 
      COMMON/ML5_0_VALIDH/VALIDH
331
 
 
332
 
      COMPLEX*16 AMPL(3,NLOOPAMPS)
333
 
      LOGICAL S(NLOOPAMPS)
334
 
      COMMON/ML5_0_AMPL/AMPL,S
335
 
 
336
 
      INTEGER CF_D(NLOOPAMPS,NBORNAMPS)
337
 
      INTEGER CF_N(NLOOPAMPS,NBORNAMPS)
338
 
      COMMON/ML5_0_CF/CF_D,CF_N
339
 
 
340
 
      INTEGER HELC(NEXTERNAL,NCOMB)
341
 
      COMMON/ML5_0_HELCONFIGS/HELC
342
 
 
343
 
      REAL*8 PREC,USER_STAB_PREC
344
 
      DATA USER_STAB_PREC/-1.0D0/
345
 
      COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC
346
 
 
347
 
C     Return codes H,T,U correspond to the hundreds, tens and units
348
 
C     building returncode, i.e.
349
 
C     RETURNCODE=100*RET_CODE_H+10*RET_CODE_T+RET_CODE_U
350
 
 
351
 
      INTEGER RET_CODE_H,RET_CODE_T,RET_CODE_U
352
 
      REAL*8 ACCURACY(0:NSQUAREDSO)
353
 
      DATA (ACCURACY(I),I=0,NSQUAREDSO)/NSQUAREDSOP1*1.0D0/
354
 
      DATA RET_CODE_H,RET_CODE_T,RET_CODE_U/1,1,0/
355
 
      COMMON/ML5_0_ACC/ACCURACY,RET_CODE_H,RET_CODE_T,RET_CODE_U
356
 
 
357
 
C     Allows to forbid the zero helicity double check, no matter the
358
 
C      value in MadLoopParams.dat
359
 
C     This can be accessed with the SET_FORBID_HEL_DOUBLECHECK
360
 
C      subroutine of MadLoopCommons.dat
361
 
      LOGICAL FORBID_HEL_DOUBLECHECK
362
 
      COMMON/FORBID_HEL_DOUBLECHECK/FORBID_HEL_DOUBLECHECK
363
 
 
364
 
      LOGICAL MP_DONE_ONCE
365
 
      DATA MP_DONE_ONCE/.FALSE./
366
 
      COMMON/ML5_0_MP_DONE_ONCE/MP_DONE_ONCE
367
 
 
368
 
      CHARACTER(512) MLPATH
369
 
      COMMON/MLPATH/MLPATH
370
 
 
371
 
      LOGICAL ML_INIT
372
 
      COMMON/ML_INIT/ML_INIT
373
 
 
374
 
C     This variable controls the *local* initialization of this
375
 
C      particular SubProcess.
376
 
C     For example, the reading of the filters must be done
377
 
C      independently by each SubProcess.
378
 
      LOGICAL LOCAL_ML_INIT
379
 
      DATA LOCAL_ML_INIT/.TRUE./
380
 
 
381
 
C     Variables related to turning off the Lorentz rotation test when
382
 
C      spin-2 particles are external
383
 
      LOGICAL WARNED_LORENTZ_STAB_TEST_OFF
384
 
      DATA WARNED_LORENTZ_STAB_TEST_OFF/.FALSE./
385
 
      INTEGER NROTATIONS_DP_BU,NROTATIONS_QP_BU
386
 
 
387
 
C     This array specify potential special requirements on the
388
 
C      helicities to
389
 
C     consider. POLARIZATIONS(0,0) is -1 if there is not such
390
 
C      requirement.
391
 
      INTEGER POLARIZATIONS(0:NEXTERNAL,0:5)
392
 
      COMMON/ML5_0_BEAM_POL/POLARIZATIONS
393
 
 
394
 
C     ----------
395
 
C     BEGIN CODE
396
 
C     ----------
397
 
 
398
 
      IF(ML_INIT) THEN
399
 
        CALL PRINT_MADLOOP_BANNER()
400
 
        TMP = 'auto'
401
 
        CALL SETMADLOOPPATH(TMP)
402
 
        CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN)
403
 
        CALL MADLOOPPARAMREADER(PARAMFN,.TRUE.)
404
 
        IF (FORBID_HEL_DOUBLECHECK) THEN
405
 
          DOUBLECHECKHELICITYFILTER = .FALSE.
406
 
        ENDIF
407
 
        ML_INIT = .FALSE.
408
 
C       For now only CutTools is interfaced in the default mode.
409
 
C        Samurai could follow.
410
 
        DO I=1,SIZE(MLREDUCTIONLIB)
411
 
          IF (MLREDUCTIONLIB(I).EQ.1) THEN
412
 
            FOUND_VALID_REDUCTION_METHOD = .TRUE.
413
 
          ENDIF
414
 
        ENDDO
415
 
        IF (.NOT.FOUND_VALID_REDUCTION_METHOD) THEN
416
 
          WRITE(*,*) 'ERROR:: For now, only CutTools is interfaced to'
417
 
     $     //' MadLoop in the non-optimized output.'
418
 
          WRITE(*,*) 'ERROR:: Make sure to include 1 in the parameter'
419
 
     $     //' MLReductionLib of the card MadLoopParams.dat'
420
 
          STOP 1
421
 
        ENDIF
422
 
      ENDIF
423
 
      IF (LOCAL_ML_INIT) THEN
424
 
C       Setup the file paths
425
 
        CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN)
426
 
        CALL JOINPATH(MLPATH,PROC_PREFIX,TMP)
427
 
        CALL JOINPATH(TMP,HELCONFIGFNAME,HELCONFIGFN)
428
 
        CALL JOINPATH(TMP,LOOPFILTERFNAME,LOOPFILTERFN)
429
 
        CALL JOINPATH(TMP,COLORNUMFNAME,COLORNUMFN)
430
 
        CALL JOINPATH(TMP,COLORDENOMFNAME,COLORDENOMFN)
431
 
        CALL JOINPATH(TMP,HELFILTERFNAME,HELFILTERFN)
432
 
 
433
 
C       Make sure that the loop filter is disabled when there is
434
 
C        spin-2 particles for 2>1 or 1>2 processes
435
 
        IF(MAX_SPIN_EXTERNAL_PARTICLE.GT.3.AND.(NEXTERNAL.LE.3.AND.HELI
436
 
     $CITYFILTERLEVEL.NE.0)) THEN
437
 
          WRITE(*,*) '##INFO: Helicity filter deactivated for 2>1'
438
 
     $     //' processes involving spin 2 particles.'
439
 
          HELICITYFILTERLEVEL = 0
440
 
C         We write a dummy filter for structural reasons here
441
 
          OPEN(1, FILE=HELFILTERFN, ERR=6116, STATUS='NEW'
442
 
     $     ,ACTION='WRITE')
443
 
          DO I=1,NCOMB
444
 
            WRITE(1,*) 'T'
445
 
          ENDDO
446
 
 6116     CONTINUE
447
 
          CLOSE(1)
448
 
        ENDIF
449
 
 
450
 
        OPEN(1, FILE=COLORNUMFN, ERR=104, STATUS='OLD',          
451
 
     $    ACTION='READ')
452
 
        DO I=1,NLOOPAMPS
453
 
          READ(1,*,END=105) (CF_N(I,J),J=1,NBORNAMPS)
454
 
        ENDDO
455
 
        GOTO 105
456
 
 104    CONTINUE
457
 
        STOP 'Color factors could not be initialized from file'
458
 
     $   //' ML5_0_ColorNumFactors.dat. File not found'
459
 
 105    CONTINUE
460
 
        CLOSE(1)
461
 
        OPEN(1, FILE=COLORDENOMFN, ERR=106, STATUS='OLD',          
462
 
     $    ACTION='READ')
463
 
        DO I=1,NLOOPAMPS
464
 
          READ(1,*,END=107) (CF_D(I,J),J=1,NBORNAMPS)
465
 
        ENDDO
466
 
        GOTO 107
467
 
 106    CONTINUE
468
 
        STOP 'Color factors could not be initialized from file'
469
 
     $   //' ML5_0_ColorDenomFactors.dat. File not found'
470
 
 107    CONTINUE
471
 
        CLOSE(1)
472
 
        OPEN(1, FILE=HELCONFIGFN, ERR=108, STATUS='OLD',              
473
 
     $       ACTION='READ')
474
 
        DO H=1,NCOMB
475
 
          READ(1,*,END=109) (HELC(I,H),I=1,NEXTERNAL)
476
 
        ENDDO
477
 
        GOTO 109
478
 
 108    CONTINUE
479
 
        STOP 'Color helictiy configurations could not be initialized'
480
 
     $   //' from file ML5_0_HelConfigs.dat. File not found'
481
 
 109    CONTINUE
482
 
        CLOSE(1)
483
 
        IF(BOOTANDSTOP) THEN
484
 
          WRITE(*,*) '##Stopped by user request.'
485
 
          STOP
486
 
        ENDIF
487
 
        LOCAL_ML_INIT = .FALSE.
488
 
      ENDIF
489
 
 
490
 
C     Make sure that lorentz rotation tests are not used if there is
491
 
C      external loop wavefunction of spin 2 and that one specific
492
 
C      helicity is asked
493
 
      NROTATIONS_DP_BU = NROTATIONS_DP
494
 
      NROTATIONS_QP_BU = NROTATIONS_QP
495
 
      IF(MAX_SPIN_EXTERNAL_PARTICLE.GT.3.AND.USERHEL.NE.-1) THEN
496
 
        IF(.NOT.WARNED_LORENTZ_STAB_TEST_OFF) THEN
497
 
          WRITE(*,*) '##WARNING: Evaluation of a specific helicity was'
498
 
     $     //' asked for this PS point, and there is a spin-2 (or'
499
 
     $     //' higher) particle in the external states.'
500
 
          WRITE(*,*) '##WARNING: As a result, MadLoop disabled the'
501
 
     $     //' Lorentz rotation test for this phase-space point only.'
502
 
          WRITE(*,*) '##WARNING: Further warning of that type'
503
 
     $     //' suppressed.'
504
 
          WARNED_LORENTZ_STAB_TEST_OFF = .FALSE.
505
 
        ENDIF
506
 
        NROTATIONS_QP=0
507
 
        NROTATIONS_DP=0
508
 
      ENDIF
509
 
 
510
 
      IF(NTRY.EQ.0) THEN
511
 
        CALL ML5_0_SET_N_EVALS(N_DP_EVAL,N_QP_EVAL)
512
 
        HELDOUBLECHECKED=(.NOT.DOUBLECHECKHELICITYFILTER)
513
 
     $   .OR.(HELICITYFILTERLEVEL.EQ.0)
514
 
        DO J=1,NCOMB
515
 
          DO I=1,NCTAMPS
516
 
            GOODAMP(I,J)=.TRUE.
517
 
          ENDDO
518
 
        ENDDO
519
 
        OPEN(1, FILE=LOOPFILTERFN, ERR=100, STATUS='OLD',          
520
 
     $    ACTION='READ')
521
 
        DO J=1,NCOMB
522
 
          READ(1,*,END=101) (GOODAMP(I,J),I=NCTAMPS+1,NLOOPAMPS)
523
 
        ENDDO
524
 
        GOTO 101
525
 
 100    CONTINUE
526
 
        FOUNDLOOPFILTER=.FALSE.
527
 
        DO J=1,NCOMB
528
 
          DO I=NCTAMPS+1,NLOOPAMPS
529
 
            GOODAMP(I,J)=(.NOT.USELOOPFILTER)
530
 
          ENDDO
531
 
        ENDDO
532
 
 101    CONTINUE
533
 
        CLOSE(1)
534
 
        IF (HELICITYFILTERLEVEL.EQ.0) THEN
535
 
          FOUNDHELFILTER=.TRUE.
536
 
          DO J=1,NCOMB
537
 
            GOODHEL(J)=.TRUE.
538
 
          ENDDO
539
 
          GOTO 122
540
 
        ENDIF
541
 
        OPEN(1, FILE=HELFILTERFN, ERR=102, STATUS='OLD',          
542
 
     $    ACTION='READ')
543
 
        READ(1,*,END=103) (GOODHEL(I),I=1,NCOMB)
544
 
        GOTO 103
545
 
 102    CONTINUE
546
 
        FOUNDHELFILTER=.FALSE.
547
 
        DO J=1,NCOMB
548
 
          GOODHEL(J)=.TRUE.
549
 
        ENDDO
550
 
 103    CONTINUE
551
 
        CLOSE(1)
552
 
 122    CONTINUE
553
 
      ENDIF
554
 
 
555
 
      MP_DONE=.FALSE.
556
 
      MP_DONE_ONCE=.FALSE.
557
 
      MP_PS_SET=.FALSE.
558
 
      STAB_INDEX=0
559
 
      DOING_QP_EVALS=.FALSE.
560
 
      EVAL_DONE(1)=.TRUE.
561
 
      DO I=2,MAXSTABILITYLENGTH
562
 
        EVAL_DONE(I)=.FALSE.
563
 
      ENDDO
564
 
 
565
 
C     Compute the born, for a specific helicity if asked so.
566
 
      CALL ML5_0_SMATRIXHEL(P_USER,USERHEL,ANS(0))
567
 
 
568
 
 
569
 
      IF (USER_STAB_PREC.GT.0.0D0) THEN
570
 
        MLSTABTHRES_BU=MLSTABTHRES
571
 
        MLSTABTHRES=USER_STAB_PREC
572
 
C       In the initialization, I cannot perform stability test and
573
 
C        therefore guarantee any precision
574
 
        CTMODEINIT_BU=CTMODEINIT
575
 
C       So either one choses quad precision directly
576
 
C       CTMODEINIT=4
577
 
C       Or, because this is very slow, we keep the orignal value. The
578
 
C        accuracy returned is -1 and tells the MC that he should not
579
 
C        trust the evaluation for checks.
580
 
        CTMODEINIT=CTMODEINIT_BU
581
 
      ENDIF
582
 
 
583
 
      IF(.NOT.BYPASS_CHECK) THEN
584
 
        NTRY=NTRY+1
585
 
      ENDIF
586
 
 
587
 
      IF(DONEHELDOUBLECHECK.AND.(.NOT.HELDOUBLECHECKED)) THEN
588
 
        HELDOUBLECHECKED=.TRUE.
589
 
        DONEHELDOUBLECHECK=.FALSE.
590
 
      ENDIF
591
 
 
592
 
      CHECKPHASE=(NTRY.LE.CHECKCYCLE).AND.(((.NOT.FOUNDLOOPFILTER)
593
 
     $ .AND.USELOOPFILTER).OR.(.NOT.FOUNDHELFILTER))
594
 
 
595
 
      IF (WRITEOUTFILTERS) THEN
596
 
        IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDHELFILTER)) THEN
597
 
          OPEN(1, FILE=HELFILTERFN, ERR=110, STATUS='NEW'
598
 
     $     ,ACTION='WRITE')
599
 
          WRITE(1,*) (GOODHEL(I),I=1,NCOMB)
600
 
 110      CONTINUE
601
 
          CLOSE(1)
602
 
          FOUNDHELFILTER=.TRUE.
603
 
        ENDIF
604
 
 
605
 
        IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDLOOPFILTER)
606
 
     $   .AND.USELOOPFILTER) THEN
607
 
          OPEN(1, FILE=LOOPFILTERFN, ERR=111, STATUS='NEW'
608
 
     $     ,ACTION='WRITE')
609
 
          DO J=1,NCOMB
610
 
            WRITE(1,*) (GOODAMP(I,J),I=NCTAMPS+1,NLOOPAMPS)
611
 
          ENDDO
612
 
 111      CONTINUE
613
 
          CLOSE(1)
614
 
          FOUNDLOOPFILTER=.TRUE.
615
 
        ENDIF
616
 
      ENDIF
617
 
 
618
 
      IF (BYPASS_CHECK) THEN
619
 
        OLD_CHECKPHASE = CHECKPHASE
620
 
        OLD_HELDOUBLECHECKED = HELDOUBLECHECKED
621
 
        CHECKPHASE = .FALSE.
622
 
        HELDOUBLECHECKED = .TRUE.
623
 
        DO I=1,NCOMB
624
 
          OLD_GOODHEL(I)=GOODHEL(I)
625
 
          GOODHEL(I) = .TRUE.
626
 
        ENDDO
627
 
        DO I=1,NCOMB
628
 
          DO J=1,NLOOPAMPS
629
 
            OLD_GOODAMP(J,I)=GOODAMP(J,I)
630
 
            GOODAMP(J,I) = .TRUE.
631
 
          ENDDO
632
 
        ENDDO
633
 
      ENDIF
634
 
 
635
 
      IF(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
636
 
        HELPICKED=1
637
 
        CTMODE=CTMODEINIT
638
 
      ELSE
639
 
        IF (USERHEL.NE.-1) THEN
640
 
          IF(.NOT.GOODHEL(USERHEL)) THEN
641
 
            ANS(1)=0.0D0
642
 
            ANS(2)=0.0D0
643
 
            ANS(3)=0.0D0
644
 
            GOTO 9999
645
 
          ENDIF
646
 
        ENDIF
647
 
        HELPICKED=USERHEL
648
 
        IF (CTMODERUN.GT.-1) THEN
649
 
          CTMODE=CTMODERUN
650
 
        ELSE
651
 
          CTMODE=1
652
 
        ENDIF
653
 
      ENDIF
654
 
 
655
 
      DO I=1,NEXTERNAL
656
 
        DO J=0,3
657
 
          PS(J,I)=P_USER(J,I)
658
 
        ENDDO
659
 
      ENDDO
660
 
 
661
 
      IF (IMPROVEPSPOINT.GE.0) THEN
662
 
C       Make the input PS more precise (exact onshell and
663
 
C        energy-momentum conservation)
664
 
        CALL ML5_0_IMPROVE_PS_POINT_PRECISION(PS)
665
 
      ENDIF
666
 
 
667
 
      DO I=1,NEXTERNAL
668
 
        DO J=0,3
669
 
          P(J,I)=PS(J,I)
670
 
        ENDDO
671
 
      ENDDO
672
 
 
673
 
      DO K=1, 3
674
 
        BUFFR(K)=0.0D0
675
 
        DO I=1,NLOOPAMPS
676
 
          AMPL(K,I)=(0.0D0,0.0D0)
677
 
        ENDDO
678
 
      ENDDO
679
 
 
680
 
      LSCALE=DSQRT(ABS((P(0,1)+P(0,2))**2-(P(1,1)+P(1,2))**2-(P(2,1)
681
 
     $ +P(2,2))**2-(P(3,1)+P(3,2))**2))
682
 
 
683
 
C     We chose to use the born evaluation for the reference
684
 
      CALL ML5_0_SMATRIX(P,REF)
685
 
 
686
 
 200  CONTINUE
687
 
 
688
 
      IF (CTMODE.EQ.0.OR.CTMODE.GE.4) THEN
689
 
        CALL MP_UPDATE_AS_PARAM()
690
 
      ENDIF
691
 
 
692
 
      IF (.NOT.MP_PS_SET.AND.(CTMODE.EQ.0.OR.CTMODE.GE.4)) THEN
693
 
        CALL ML5_0_SET_MP_PS(P_USER)
694
 
        MP_PS_SET = .TRUE.
695
 
      ENDIF
696
 
 
697
 
      DO K=1,3
698
 
        ANS(K)=0.0D0
699
 
      ENDDO
700
 
 
701
 
      VALIDH=-1
702
 
      DO H=1,NCOMB
703
 
        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(
704
 
     $.NOT.HELDOUBLECHECKED).OR.GOODHEL(H)))) THEN
705
 
 
706
 
C         Handle the possible requirement of specific polarizations
707
 
          IF ((.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.POLARIZATIONS(
708
 
     $0,0).EQ.0.AND.(.NOT.ML5_0_IS_HEL_SELECTED(H))) THEN
709
 
            CYCLE
710
 
          ENDIF
711
 
 
712
 
          IF (VALIDH.EQ.-1) VALIDH=H
713
 
          DO I=1,NEXTERNAL
714
 
            NHEL(I)=HELC(I,H)
715
 
          ENDDO
716
 
C         Check if we are in multiple precision and compute wfs and
717
 
C          amps accordingly if needed   
718
 
          IF (CTMODE.GE.4) THEN
719
 
C           Force that only current helicity is used in the routine
720
 
C            below
721
 
C           This should always be done, even if MP_DONE is True
722
 
C           because the AMPL of the R2 MUST be recomputed for loop
723
 
C            induced.
724
 
C           (because they are not saved for each hel configuration)
725
 
C           (This is not optimal unlike what is done int the loop
726
 
C            optimized output)
727
 
            HELPICKED_BU = HELPICKED
728
 
            HELPICKED = H
729
 
            CALL ML5_0_MP_BORN_AMPS_AND_WFS(MP_P)
730
 
            HELPICKED = HELPICKED_BU
731
 
            GOTO 300
732
 
          ENDIF
733
 
          CALL IXXXXX(P(0,1),ZERO,NHEL(1),+1*IC(1),W(1,1,H))
734
 
          CALL OXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2,H))
735
 
          CALL OXXXXX(P(0,3),ZERO,NHEL(3),+1*IC(3),W(1,3,H))
736
 
          CALL IXXXXX(P(0,4),ZERO,NHEL(4),-1*IC(4),W(1,4,H))
737
 
          CALL VXXXXX(P(0,5),ZERO,NHEL(5),+1*IC(5),W(1,5,H))
738
 
          CALL FFV1_2(W(1,1,H),W(1,5,H),GC_5,ZERO,ZERO,W(1,6,H))
739
 
          CALL FFV2_3(W(1,4,H),W(1,3,H),GC_47,MDL_MW,MDL_WW,W(1,7,H))
740
 
C         Amplitude(s) for born diagram with ID 1
741
 
          CALL FFV2_0(W(1,6,H),W(1,2,H),W(1,7,H),GC_47,AMP(1,H))
742
 
          CALL FFV1_1(W(1,2,H),W(1,5,H),GC_5,ZERO,ZERO,W(1,8,H))
743
 
C         Amplitude(s) for born diagram with ID 2
744
 
          CALL FFV2_0(W(1,1,H),W(1,8,H),W(1,7,H),GC_47,AMP(2,H))
745
 
          CALL FFV2_1(W(1,2,H),W(1,7,H),GC_47,ZERO,ZERO,W(1,9,H))
746
 
C         Counter-term amplitude(s) for loop diagram number 3
747
 
          CALL R2_QQ_1_0(W(1,6,H),W(1,9,H),R2_QQQ,AMPL(1,1))
748
 
C         Counter-term amplitude(s) for loop diagram number 4
749
 
          CALL FFV2_0(W(1,6,H),W(1,2,H),W(1,7,H),R2_SXCW,AMPL(1,2))
750
 
C         Counter-term amplitude(s) for loop diagram number 5
751
 
          CALL FFV2_0(W(1,1,H),W(1,8,H),W(1,7,H),R2_SXCW,AMPL(1,3))
752
 
C         Counter-term amplitude(s) for loop diagram number 7
753
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,4)
754
 
     $     )
755
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,5)
756
 
     $     )
757
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,6)
758
 
     $     )
759
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,7)
760
 
     $     )
761
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQB,AMPL(1,8))
762
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2,9)
763
 
     $     )
764
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQT,AMPL(1,10))
765
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
766
 
     $     ,11))
767
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),UV_GQQG_1EPS,AMPL(2
768
 
     $     ,12))
769
 
          CALL FFV1_0(W(1,1,H),W(1,9,H),W(1,5,H),R2_GQQ,AMPL(1,13))
770
 
          CALL FFV2_2(W(1,1,H),W(1,7,H),GC_47,ZERO,ZERO,W(1,10,H))
771
 
C         Counter-term amplitude(s) for loop diagram number 11
772
 
          CALL R2_QQ_1_0(W(1,10,H),W(1,8,H),R2_QQQ,AMPL(1,14))
773
 
C         Counter-term amplitude(s) for loop diagram number 12
774
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
775
 
     $     ,15))
776
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
777
 
     $     ,16))
778
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
779
 
     $     ,17))
780
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
781
 
     $     ,18))
782
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQB,AMPL(1,19))
783
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
784
 
     $     ,20))
785
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQT,AMPL(1,21))
786
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQQ_1EPS,AMPL(2
787
 
     $     ,22))
788
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),UV_GQQG_1EPS,AMPL(2
789
 
     $     ,23))
790
 
          CALL FFV1_0(W(1,10,H),W(1,2,H),W(1,5,H),R2_GQQ,AMPL(1,24))
791
 
C         Amplitude(s) for UVCT diagram with ID 14
792
 
          CALL FFV2_0(W(1,6,H),W(1,2,H),W(1,7,H),GC_47,AMPL(1,25))
793
 
          AMPL(1,25)=AMPL(1,25)*(1.0D0*UVWFCT_G_2+1.0D0*UVWFCT_G_1)
794
 
C         Amplitude(s) for UVCT diagram with ID 15
795
 
          CALL FFV2_0(W(1,6,H),W(1,2,H),W(1,7,H),GC_47,AMPL(2,26))
796
 
          AMPL(2,26)=AMPL(2,26)*(2.0D0*UVWFCT_G_2_1EPS)
797
 
C         Amplitude(s) for UVCT diagram with ID 16
798
 
          CALL FFV2_0(W(1,1,H),W(1,8,H),W(1,7,H),GC_47,AMPL(1,27))
799
 
          AMPL(1,27)=AMPL(1,27)*(1.0D0*UVWFCT_G_2+1.0D0*UVWFCT_G_1)
800
 
C         Amplitude(s) for UVCT diagram with ID 17
801
 
          CALL FFV2_0(W(1,1,H),W(1,8,H),W(1,7,H),GC_47,AMPL(2,28))
802
 
          AMPL(2,28)=AMPL(2,28)*(2.0D0*UVWFCT_G_2_1EPS)
803
 
 300      CONTINUE
804
 
 
805
 
 
806
 
 
807
 
          DO I=1,NCTAMPS
808
 
            DO J=1,NBORNAMPS
809
 
              CFTOT=DCMPLX(CF_N(I,J)/DBLE(ABS(CF_D(I,J))),0.0D0)
810
 
              IF(CF_D(I,J).LT.0) CFTOT=CFTOT*IMAG1
811
 
              DO K=1,3
812
 
                ANS(K)=ANS(K)+2.0D0*DBLE(CFTOT*AMPL(K,I)*DCONJG(AMP(J
813
 
     $           ,H)))
814
 
              ENDDO
815
 
            ENDDO
816
 
          ENDDO
817
 
        ENDIF
818
 
      ENDDO
819
 
 
820
 
C     WHEN CTMODE IS >=4, then the MP computation of wfs and amps is
821
 
C      automatically done.
822
 
      IF (CTMODE.GE.4) THEN
823
 
        MP_DONE = .TRUE.
824
 
      ENDIF
825
 
 
826
 
      IF(SKIPLOOPEVAL) THEN
827
 
        GOTO 1226
828
 
      ENDIF
829
 
 
830
 
C     Loop amplitude for loop diagram with ID 3
831
 
      CALL ML5_0_LOOP_2_2(1,6,9,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
832
 
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_5
833
 
     $ ,MP__GC_5,1,1,1,29,AMPL(1,29),S(29))
834
 
C     Loop amplitude for loop diagram with ID 4
835
 
      CALL ML5_0_LOOP_3_3(2,2,7,6,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
836
 
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
837
 
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_47,MP__GC_47,GC_5
838
 
     $ ,MP__GC_5,2,1,1,30,AMPL(1,30),S(30))
839
 
C     Loop amplitude for loop diagram with ID 5
840
 
      CALL ML5_0_LOOP_3_3(3,1,7,8,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
841
 
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
842
 
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_47,MP__GC_47,GC_5
843
 
     $ ,MP__GC_5,2,1,1,31,AMPL(1,31),S(31))
844
 
C     Loop amplitude for loop diagram with ID 6
845
 
      CALL ML5_0_LOOP_4_4(4,1,2,5,7,DCMPLX(ZERO),CMPLX(MP__ZERO
846
 
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
847
 
     $ ,CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
848
 
     $ ,GC_5,MP__GC_5,GC_5,MP__GC_5,GC_5,MP__GC_5,GC_47,MP__GC_47,3,1
849
 
     $ ,1,32,AMPL(1,32),S(32))
850
 
C     Loop amplitude for loop diagram with ID 7
851
 
      CALL ML5_0_LOOP_3_3(5,1,5,9,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
852
 
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
853
 
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_4,MP__GC_4,GC_5
854
 
     $ ,MP__GC_5,2,1,1,33,AMPL(1,33),S(33))
855
 
C     Loop amplitude for loop diagram with ID 8
856
 
      CALL ML5_0_LOOP_4_4(6,1,5,2,7,DCMPLX(ZERO),CMPLX(MP__ZERO
857
 
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
858
 
     $ ,CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
859
 
     $ ,GC_5,MP__GC_5,GC_4,MP__GC_4,GC_5,MP__GC_5,GC_47,MP__GC_47,3,1
860
 
     $ ,1,34,AMPL(1,34),S(34))
861
 
C     Loop amplitude for loop diagram with ID 9
862
 
      CALL ML5_0_LOOP_4_4(7,1,2,7,5,DCMPLX(ZERO),CMPLX(MP__ZERO
863
 
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
864
 
     $ ,CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
865
 
     $ ,GC_5,MP__GC_5,GC_5,MP__GC_5,GC_47,MP__GC_47,GC_5,MP__GC_5,3,1
866
 
     $ ,1,35,AMPL(1,35),S(35))
867
 
C     Loop amplitude for loop diagram with ID 10
868
 
      CALL ML5_0_LOOP_3_3(8,1,5,9,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
869
 
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
870
 
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5,GC_5
871
 
     $ ,MP__GC_5,2,1,1,36,AMPL(1,36),S(36))
872
 
C     Loop amplitude for loop diagram with ID 11
873
 
      CALL ML5_0_LOOP_2_2(9,8,10,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16)
874
 
     $ ,DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_5
875
 
     $ ,MP__GC_5,1,1,1,37,AMPL(1,37),S(37))
876
 
C     Loop amplitude for loop diagram with ID 12
877
 
      CALL ML5_0_LOOP_3_3(10,2,5,10,DCMPLX(ZERO),CMPLX(MP__ZERO
878
 
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
879
 
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_4,MP__GC_4,GC_5
880
 
     $ ,MP__GC_5,2,1,1,38,AMPL(1,38),S(38))
881
 
C     Loop amplitude for loop diagram with ID 13
882
 
      CALL ML5_0_LOOP_3_3(11,2,5,10,DCMPLX(ZERO),CMPLX(MP__ZERO
883
 
     $ ,KIND=16),DCMPLX(ZERO),CMPLX(MP__ZERO,KIND=16),DCMPLX(ZERO)
884
 
     $ ,CMPLX(MP__ZERO,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5,GC_5
885
 
     $ ,MP__GC_5,2,1,1,39,AMPL(1,39),S(39))
886
 
 
887
 
      DO I=NCTAMPS+1,NLOOPAMPS
888
 
        ANS(1)=ANS(1)+AMPL(1,I)
889
 
        ANS(2)=ANS(2)+AMPL(2,I)
890
 
        ANS(3)=ANS(3)+AMPL(3,I)
891
 
        IF((CTMODERUN.NE.-1).AND..NOT.CHECKPHASE.AND.(.NOT.S(I))) THEN
892
 
          WRITE(*,*) '##W03 WARNING Contribution ',I,' is unstable.'
893
 
        ENDIF
894
 
      ENDDO
895
 
 
896
 
 1226 CONTINUE
897
 
 
898
 
      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
899
 
C       Update of NEXTREF, will be used for loop induced only.
900
 
        NEXTREF = NEXTREF + ANS(1) + ANS(2) + ANS(3)
901
 
        IF((USERHEL.EQ.-1).OR.(USERHEL.EQ.HELPICKED)) THEN
902
 
          BUFFR(1)=BUFFR(1)+ANS(1)
903
 
          BUFFR(2)=BUFFR(2)+ANS(2)
904
 
          BUFFR(3)=BUFFR(3)+ANS(3)
905
 
        ENDIF
906
 
 
907
 
        IF (CHECKPHASE) THEN
908
 
C         SET THE HELICITY FILTER
909
 
          IF(.NOT.FOUNDHELFILTER) THEN
910
 
            IF(ML5_0_ISZERO(ABS(ANS(1))+ABS(ANS(2))+ABS(ANS(3)),REF
911
 
     $       /DBLE(NCOMB),-1)) THEN
912
 
              IF(NTRY.EQ.1) THEN
913
 
                GOODHEL(HELPICKED)=.FALSE.
914
 
              ELSEIF(GOODHEL(HELPICKED)) THEN
915
 
                WRITE(*,*) '##W02A WARNING Inconsistent helicity '
916
 
     $           ,HELPICKED
917
 
                IF(HELINITSTARTOVER) THEN
918
 
                  WRITE(*,*) '##I01 INFO Initialization starting over'
919
 
     $             //' because of inconsistency in the helicity filter'
920
 
     $             //' setup.'
921
 
                  NTRY=0
922
 
                ENDIF
923
 
              ENDIF
924
 
            ELSE
925
 
              IF(.NOT.GOODHEL(HELPICKED)) THEN
926
 
                WRITE(*,*) '##W02B WARNING Inconsistent helicity '
927
 
     $           ,HELPICKED
928
 
                IF(HELINITSTARTOVER) THEN
929
 
                  WRITE(*,*) '##I01 INFO Initialization starting over'
930
 
     $             //' because of inconsistency in the helicity filter'
931
 
     $             //' setup.'
932
 
                  NTRY=0
933
 
                ELSE
934
 
                  GOODHEL(HELPICKED)=.TRUE.
935
 
                ENDIF
936
 
              ENDIF
937
 
            ENDIF
938
 
          ENDIF
939
 
 
940
 
C         SET THE LOOP FILTER
941
 
          IF(.NOT.FOUNDLOOPFILTER.AND.USELOOPFILTER) THEN
942
 
            DO I=NCTAMPS+1,NLOOPAMPS
943
 
              IF(.NOT.ML5_0_ISZERO(ABS(AMPL(1,I))+ABS(AMPL(2,I))
944
 
     $         +ABS(AMPL(3,I)),(REF*1.0D-4),I)) THEN
945
 
                IF(NTRY.EQ.1) THEN
946
 
                  GOODAMP(I,HELPICKED)=.TRUE.
947
 
                ELSEIF(.NOT.GOODAMP(I,HELPICKED)) THEN
948
 
                  WRITE(*,*) '##W02 WARNING Inconsistent loop amp ',I
949
 
     $             ,' for helicity ',HELPICKED,'.'
950
 
                  IF(LOOPINITSTARTOVER) THEN
951
 
                    WRITE(*,*) '##I01 INFO Initialization starting'
952
 
     $               //' over because of inconsistency in the loop'
953
 
     $               //' filter setup.'
954
 
                    NTRY=0
955
 
                  ELSE
956
 
                    GOODAMP(I,HELPICKED)=.TRUE.
957
 
                  ENDIF
958
 
                ENDIF
959
 
              ENDIF
960
 
            ENDDO
961
 
          ENDIF
962
 
        ELSEIF (.NOT.HELDOUBLECHECKED)THEN
963
 
          IF ((.NOT.GOODHEL(HELPICKED)).AND.(.NOT.ML5_0_ISZERO(ABS(ANS(
964
 
     $1))+ABS(ANS(2))+ABS(ANS(3)),REF/DBLE(NCOMB),-1))) THEN
965
 
            WRITE(*,*) '##W15 Helicity filter could not be'
966
 
     $       //' successfully double checked.'
967
 
            WRITE(*,*) '##One reason for this is that you have changed'
968
 
     $       //' sensible parameters which affected what are the zero'
969
 
     $       //' helicity configurations.'
970
 
            WRITE(*,*) '##MadLoop will try to reset the Helicity'
971
 
     $       //' filter with the next PS points it receives.'
972
 
            NTRY=0
973
 
            OPEN(30,FILE=HELFILTERFN,ERR=349)
974
 
 349        CONTINUE
975
 
            CLOSE(30,STATUS='delete')
976
 
          ENDIF
977
 
C         SET HELDOUBLECHECKED TO .TRUE. WHEN DONE
978
 
C         even if it failed we do not want to redo the check
979
 
C          afterwards if HELINITSTARTOVER=.FALSE.
980
 
          IF (HELPICKED.EQ.NCOMB.AND.(NTRY.NE.0.OR..NOT.HELINITSTARTOVE
981
 
     $R)) THEN
982
 
            DONEHELDOUBLECHECK=.TRUE.
983
 
          ENDIF
984
 
        ENDIF
985
 
 
986
 
C       GOTO NEXT HELICITY OR FINISH
987
 
        IF(HELPICKED.NE.NCOMB) THEN
988
 
          HELPICKED=HELPICKED+1
989
 
          MP_DONE=.FALSE.
990
 
          GOTO 200
991
 
        ELSE
992
 
          ANS(1)=BUFFR(1)
993
 
          ANS(2)=BUFFR(2)
994
 
          ANS(3)=BUFFR(3)
995
 
C         We add one here to the number of PS points used for building
996
 
C          the reference scale for comparison (used only for
997
 
C          loop-induced processes).
998
 
          NPSPOINTS = NPSPOINTS+1
999
 
          IF(NTRY.EQ.0) THEN
1000
 
            NATTEMPTS=NATTEMPTS+1
1001
 
            IF(NATTEMPTS.EQ.MAXATTEMPTS) THEN
1002
 
              WRITE(*,*) '##E01 ERROR Could not initialize the filters'
1003
 
     $         //' in ',MAXATTEMPTS,' trials'
1004
 
              STOP
1005
 
            ENDIF
1006
 
          ENDIF
1007
 
        ENDIF
1008
 
 
1009
 
      ENDIF
1010
 
 
1011
 
      DO K=1,3
1012
 
        ANS(K)=ANS(K)/DBLE(IDEN)
1013
 
        IF (USERHEL.NE.-1) THEN
1014
 
          ANS(K)=ANS(K)*HELAVGFACTOR
1015
 
        ELSE
1016
 
          DO J=1,NINITIAL
1017
 
            IF (POLARIZATIONS(J,0).NE.-1) THEN
1018
 
              ANS(K)=ANS(K)*BEAMS_HELAVGFACTOR(J)
1019
 
              ANS(K)=ANS(K)/POLARIZATIONS(J,0)
1020
 
            ENDIF
1021
 
          ENDDO
1022
 
        ENDIF
1023
 
      ENDDO
1024
 
 
1025
 
      IF(.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.(CTMODERUN.LE.-1))
1026
 
     $  THEN
1027
 
        STAB_INDEX=STAB_INDEX+1
1028
 
        IF(DOING_QP_EVALS) THEN
1029
 
          QP_RES(1,STAB_INDEX)=ANS(1)
1030
 
          QP_RES(2,STAB_INDEX)=ANS(2)
1031
 
          QP_RES(3,STAB_INDEX)=ANS(3)
1032
 
        ELSE
1033
 
          DP_RES(1,STAB_INDEX)=ANS(1)
1034
 
          DP_RES(2,STAB_INDEX)=ANS(2)
1035
 
          DP_RES(3,STAB_INDEX)=ANS(3)
1036
 
        ENDIF
1037
 
 
1038
 
        IF(DOING_QP_EVALS) THEN
1039
 
          BASIC_CT_MODE=4
1040
 
        ELSE
1041
 
          BASIC_CT_MODE=1
1042
 
        ENDIF
1043
 
 
1044
 
C       BEGINNING OF THE DEFINITIONS OF THE DIFFERENT EVALUATION
1045
 
C        METHODS
1046
 
 
1047
 
        IF(.NOT.EVAL_DONE(2)) THEN
1048
 
          EVAL_DONE(2)=.TRUE.
1049
 
          CTMODE=BASIC_CT_MODE+1
1050
 
          GOTO 200
1051
 
        ENDIF
1052
 
 
1053
 
        CTMODE=BASIC_CT_MODE
1054
 
 
1055
 
        IF(.NOT.EVAL_DONE(3).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
1056
 
     $.1).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
1057
 
          EVAL_DONE(3)=.TRUE.
1058
 
          CALL ML5_0_ROTATE_PS(PS,P,1)
1059
 
          IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,1)
1060
 
          GOTO 200
1061
 
        ENDIF
1062
 
 
1063
 
        IF(.NOT.EVAL_DONE(4).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
1064
 
     $.2).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
1065
 
          EVAL_DONE(4)=.TRUE.
1066
 
          CALL ML5_0_ROTATE_PS(PS,P,2)
1067
 
          IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,2)
1068
 
          GOTO 200
1069
 
        ENDIF
1070
 
 
1071
 
        CALL ML5_0_ROTATE_PS(PS,P,0)
1072
 
        IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,0)
1073
 
 
1074
 
C       END OF THE DEFINITIONS OF THE DIFFERENT EVALUATION METHODS
1075
 
 
1076
 
        IF(DOING_QP_EVALS) THEN
1077
 
          CALL ML5_0_COMPUTE_ACCURACY(QP_RES,N_QP_EVAL,ACC,ANS(1))
1078
 
          ACCURACY(0)=ACC
1079
 
          RET_CODE_H=3
1080
 
          IF(ACC.GE.MLSTABTHRES) THEN
1081
 
            RET_CODE_H=4
1082
 
            NEPS=NEPS+1
1083
 
            CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,TEMP1,TEMP)
1084
 
            WRITE(*,*) '##W03 WARNING An unstable PS point was',      
1085
 
     $        ' detected.'
1086
 
            WRITE(*,*) '##(DP,QP) accuracies : (',TEMP1,',',ACC,')'
1087
 
            WRITE(*,*) '##Best estimate (fin,1eps,2eps) :',(ANS(I),I=1
1088
 
     $       ,3)
1089
 
            IF(NEPS.LE.10) THEN
1090
 
              WRITE(*,*) '##Double precision evaluations :',(DP_RES(1
1091
 
     $         ,I),I=1,N_DP_EVAL)
1092
 
              WRITE(*,*) '##Quad   precision evaluations :',(QP_RES(1
1093
 
     $         ,I),I=1,N_QP_EVAL)
1094
 
              WRITE(*,*) '##PS point specification :'
1095
 
              WRITE(*,*) '##Renormalization scale MU_R=',MU_R
1096
 
              DO I=1,NEXTERNAL
1097
 
                WRITE (*,'(i2,1x,4e27.17)') I, P(0,I),P(1,I),P(2,I)
1098
 
     $           ,P(3,I)
1099
 
              ENDDO
1100
 
            ENDIF
1101
 
            IF(NEPS.EQ.10) THEN
1102
 
              WRITE(*,*) '##Further output of the details of these'
1103
 
     $         //' unstable PS points will now be suppressed.'
1104
 
            ENDIF
1105
 
          ENDIF
1106
 
        ELSE
1107
 
          CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,ACC,ANS(1))
1108
 
          IF(ACC.GE.MLSTABTHRES) THEN
1109
 
            DOING_QP_EVALS=.TRUE.
1110
 
            EVAL_DONE(1)=.TRUE.
1111
 
            DO I=2,MAXSTABILITYLENGTH
1112
 
              EVAL_DONE(I)=.FALSE.
1113
 
            ENDDO
1114
 
            STAB_INDEX=0
1115
 
            CTMODE=4
1116
 
            GOTO 200
1117
 
          ELSE
1118
 
            ACCURACY(0)=ACC
1119
 
            RET_CODE_H=2
1120
 
          ENDIF
1121
 
        ENDIF
1122
 
      ELSE
1123
 
        RET_CODE_H=1
1124
 
        ACCURACY=-1.0D0
1125
 
      ENDIF
1126
 
 
1127
 
 9999 CONTINUE
1128
 
 
1129
 
C     Finalize the return code
1130
 
      IF (MP_DONE_ONCE) THEN
1131
 
        RET_CODE_T=2
1132
 
      ELSE
1133
 
        RET_CODE_T=1
1134
 
      ENDIF
1135
 
      IF(CHECKPHASE.OR..NOT.HELDOUBLECHECKED) THEN
1136
 
        RET_CODE_H=1
1137
 
        RET_CODE_T=RET_CODE_T+2
1138
 
        ACCURACY=-1.0D0
1139
 
      ENDIF
1140
 
      IF (RET_CODE_H.EQ.4) THEN
1141
 
        RET_CODE_U=0
1142
 
      ELSE
1143
 
        RET_CODE_U=1
1144
 
      ENDIF
1145
 
 
1146
 
C     Reinitialize the default threshold if it was specified by the
1147
 
C      user
1148
 
      IF (USER_STAB_PREC.GT.0.0D0) THEN
1149
 
        MLSTABTHRES=MLSTABTHRES_BU
1150
 
        CTMODEINIT=CTMODEINIT_BU
1151
 
      ENDIF
1152
 
 
1153
 
C     Reinitialize the Lorentz test if it had been disabled because
1154
 
C      spin-2 particles are in the external states.
1155
 
      NROTATIONS_DP = NROTATIONS_DP_BU
1156
 
      NROTATIONS_QP = NROTATIONS_QP_BU
1157
 
 
1158
 
C     Conform to the returned synthax of split orders even though the
1159
 
C      default output does not support it (this then done only for
1160
 
C      compatibility purpose).
1161
 
      ANSRETURNED(0,0)=ANS(0)
1162
 
      ANSRETURNED(1,0)=ANS(1)
1163
 
      ANSRETURNED(2,0)=ANS(2)
1164
 
      ANSRETURNED(3,0)=ANS(3)
1165
 
 
1166
 
C     Reinitialize the check phase logicals and the filters if check
1167
 
C      bypassed
1168
 
      IF (BYPASS_CHECK) THEN
1169
 
        CHECKPHASE = OLD_CHECKPHASE
1170
 
        HELDOUBLECHECKED = OLD_HELDOUBLECHECKED
1171
 
        DO I=1,NCOMB
1172
 
          GOODHEL(I)=OLD_GOODHEL(I)
1173
 
        ENDDO
1174
 
        DO I=1,NCOMB
1175
 
          DO J=1,NLOOPAMPS
1176
 
            GOODAMP(J,I)=OLD_GOODAMP(J,I)
1177
 
          ENDDO
1178
 
        ENDDO
1179
 
      ENDIF
1180
 
 
1181
 
      END
1182
 
 
1183
 
      SUBROUTINE ML5_0_COMPUTE_ACCURACY(FULLLIST, LENGTH, ACC,
1184
 
     $  ESTIMATE)
1185
 
      IMPLICIT NONE
1186
 
C     
1187
 
C     PARAMETERS 
1188
 
C     
1189
 
      INTEGER MAXSTABILITYLENGTH
1190
 
      COMMON/ML5_0_STABILITY_TESTS/MAXSTABILITYLENGTH
1191
 
C     
1192
 
C     ARGUMENTS 
1193
 
C     
1194
 
      REAL*8 FULLLIST(3,MAXSTABILITYLENGTH)
1195
 
      INTEGER LENGTH
1196
 
      REAL*8 ACC, ESTIMATE(3)
1197
 
C     
1198
 
C     LOCAL VARIABLES 
1199
 
C     
1200
 
      LOGICAL MASK(MAXSTABILITYLENGTH)
1201
 
      LOGICAL MASK3(3)
1202
 
      DATA MASK3/.TRUE.,.TRUE.,.TRUE./
1203
 
      INTEGER I,J
1204
 
      REAL*8 AVG
1205
 
      REAL*8 DIFF
1206
 
      REAL*8 ACCURACIES(3)
1207
 
      REAL*8 LIST(MAXSTABILITYLENGTH)
1208
 
 
1209
 
C     ----------
1210
 
C     BEGIN CODE
1211
 
C     ----------
1212
 
      DO I=1,LENGTH
1213
 
        MASK(I)=.TRUE.
1214
 
      ENDDO
1215
 
      DO I=LENGTH+1,MAXSTABILITYLENGTH
1216
 
        MASK(I)=.FALSE.
1217
 
      ENDDO
1218
 
 
1219
 
      DO I=1,3
1220
 
        DO J=1,MAXSTABILITYLENGTH
1221
 
          LIST(J)=FULLLIST(I,J)
1222
 
        ENDDO
1223
 
        DIFF=MAXVAL(LIST,1,MASK)-MINVAL(LIST,1,MASK)
1224
 
        AVG=(MAXVAL(LIST,1,MASK)+MINVAL(LIST,1,MASK))/2.0D0
1225
 
        ESTIMATE(I)=AVG
1226
 
        IF (AVG.EQ.0.0D0) THEN
1227
 
          ACCURACIES(I)=DIFF
1228
 
        ELSE
1229
 
          ACCURACIES(I)=DIFF/ABS(AVG)
1230
 
        ENDIF
1231
 
      ENDDO
1232
 
 
1233
 
C     The technique below is too sensitive, typically to
1234
 
C     unstablities in very small poles
1235
 
C     ACC=MAXVAL(ACCURACIES,1,MASK3)
1236
 
C     The following is used instead
1237
 
      ACC = 0.0D0
1238
 
      AVG = 0.0D0
1239
 
      DO I=1,3
1240
 
        ACC = ACC + ACCURACIES(I)*ABS(ESTIMATE(I))
1241
 
        AVG = AVG + ESTIMATE(I)
1242
 
      ENDDO
1243
 
      ACC  = ACC / ( ABS(AVG) / 3.0D0)
1244
 
 
1245
 
      END
1246
 
 
1247
 
      SUBROUTINE ML5_0_SET_N_EVALS(N_DP_EVALS,N_QP_EVALS)
1248
 
 
1249
 
      IMPLICIT NONE
1250
 
      INTEGER N_DP_EVALS, N_QP_EVALS
1251
 
 
1252
 
      INCLUDE 'MadLoopParams.inc'
1253
 
 
1254
 
      IF(CTMODERUN.LE.-1) THEN
1255
 
        N_DP_EVALS=2+NROTATIONS_DP
1256
 
        N_QP_EVALS=2+NROTATIONS_QP
1257
 
      ELSE
1258
 
        N_DP_EVALS=1
1259
 
        N_QP_EVALS=1
1260
 
      ENDIF
1261
 
 
1262
 
      IF(N_DP_EVALS.GT.20.OR.N_QP_EVALS.GT.20) THEN
1263
 
        WRITE(*,*) '##ERROR:: Increase hardcoded maxstabilitylength.'
1264
 
        STOP
1265
 
      ENDIF
1266
 
 
1267
 
      END
1268
 
 
1269
 
 
1270
 
C     THIS SUBROUTINE SIMPLY SET THE GLOBAL PS CONFIGURATION GLOBAL
1271
 
C      VARIABLES FROM A GIVEN VARIABLE IN DOUBLE PRECISION
1272
 
      SUBROUTINE ML5_0_SET_MP_PS(P)
1273
 
 
1274
 
      INTEGER    NEXTERNAL
1275
 
      PARAMETER (NEXTERNAL=5)
1276
 
      REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL)
1277
 
      COMMON/ML5_0_MP_PSPOINT/MP_PS,MP_P
1278
 
      REAL*8 P(0:3,NEXTERNAL)
1279
 
 
1280
 
      DO I=1,NEXTERNAL
1281
 
        DO J=0,3
1282
 
          MP_PS(J,I)=P(J,I)
1283
 
        ENDDO
1284
 
      ENDDO
1285
 
      CALL ML5_0_MP_IMPROVE_PS_POINT_PRECISION(MP_PS)
1286
 
      DO I=1,NEXTERNAL
1287
 
        DO J=0,3
1288
 
          MP_P(J,I)=MP_PS(J,I)
1289
 
        ENDDO
1290
 
      ENDDO
1291
 
 
1292
 
      END
1293
 
 
1294
 
      SUBROUTINE ML5_0_SET_COUPLINGORDERS_TARGET(SOTARGET)
1295
 
      IMPLICIT NONE
1296
 
C     
1297
 
C     This routine can be accessed by an external user to set the
1298
 
C      squared split order target.
1299
 
C     This functionality is only available in the optimized mode, but
1300
 
C      for compatibility 
1301
 
C     purposes, a dummy version is also put in this default output.
1302
 
C     
1303
 
C     
1304
 
C     ARGUMENTS
1305
 
C     
1306
 
      INTEGER SOTARGET
1307
 
C     ----------
1308
 
C     BEGIN CODE
1309
 
C     ----------
1310
 
      WRITE(*,*) '##WARNING:: Ignored, the possibility of selecting'
1311
 
     $ //' specific squared order contributions is not available in'
1312
 
     $ //' the default mode.'
1313
 
 
1314
 
      END
1315
 
 
1316
 
      SUBROUTINE ML5_0_FORCE_STABILITY_CHECK(ONOFF)
1317
 
C     
1318
 
C     This function can be called by the MadLoop user so as to always
1319
 
C      have stability
1320
 
C     checked, even during initialisation, when calling the *_thres
1321
 
C      routines.
1322
 
C     
1323
 
      LOGICAL ONOFF
1324
 
 
1325
 
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
1326
 
      DATA BYPASS_CHECK, ALWAYS_TEST_STABILITY /.FALSE.,.FALSE./
1327
 
      COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
1328
 
 
1329
 
      ALWAYS_TEST_STABILITY = ONOFF
1330
 
 
1331
 
      END
1332
 
 
1333
 
      SUBROUTINE ML5_0_GET_ANSWER_DIMENSION(ANSDIM)
1334
 
C     
1335
 
C     Simple subroutine which returns the upper bound of the second
1336
 
C      dimension of the
1337
 
C     quantity ANS(0:3,0:ANSDIM) returned by MadLoop. As long as the
1338
 
C      default output
1339
 
C     cannot handle split orders, this ANSDIM will always be 0.
1340
 
C     
1341
 
      INCLUDE 'nsquaredSO.inc'
1342
 
 
1343
 
      INTEGER ANSDIM
1344
 
 
1345
 
      ANSDIM=NSQUAREDSO
1346
 
 
1347
 
      END
1348
 
 
1349
 
      SUBROUTINE ML5_0_GET_NSQSO_LOOP(NSQSO)
1350
 
C     
1351
 
C     Simple subroutine returning the number of squared split order
1352
 
C     contributions returned in ANS when calling sloopmatrix 
1353
 
C     
1354
 
      INCLUDE 'nsquaredSO.inc'
1355
 
 
1356
 
      INTEGER NSQSO
1357
 
 
1358
 
      NSQSO=NSQUAREDSO
1359
 
 
1360
 
      END
1361
 
 
1362
 
      SUBROUTINE ML5_0_SET_LEG_POLARIZATION(LEG_ID, LEG_POLARIZATION)
1363
 
      IMPLICIT NONE
1364
 
C     
1365
 
C     ARGUMENTS
1366
 
C     
1367
 
      INTEGER LEG_ID
1368
 
      INTEGER LEG_POLARIZATION
1369
 
C     
1370
 
C     LOCALS
1371
 
C     
1372
 
      INTEGER I
1373
 
      INTEGER LEG_POLARIZATIONS(0:5)
1374
 
C     ----------
1375
 
C     BEGIN CODE
1376
 
C     ----------
1377
 
 
1378
 
      IF (LEG_POLARIZATION.EQ.-10000) THEN
1379
 
        LEG_POLARIZATIONS(0)=-1
1380
 
        DO I=1,5
1381
 
          LEG_POLARIZATIONS(I)=-10000
1382
 
        ENDDO
1383
 
      ELSE
1384
 
        LEG_POLARIZATIONS(0)=1
1385
 
        LEG_POLARIZATIONS(1)=LEG_POLARIZATION
1386
 
        DO I=2,5
1387
 
          LEG_POLARIZATIONS(I)=-10000
1388
 
        ENDDO
1389
 
      ENDIF
1390
 
      CALL ML5_0_SET_LEG_POLARIZATIONS(LEG_ID,LEG_POLARIZATIONS)
1391
 
 
1392
 
      END
1393
 
 
1394
 
      SUBROUTINE ML5_0_SET_LEG_POLARIZATIONS(LEG_ID, LEG_POLARIZATIONS)
1395
 
      IMPLICIT NONE
1396
 
C     
1397
 
C     CONSTANTS
1398
 
C     
1399
 
      INTEGER    NEXTERNAL
1400
 
      PARAMETER (NEXTERNAL=5)
1401
 
      INTEGER NPOLENTRIES
1402
 
      PARAMETER (NPOLENTRIES=(NEXTERNAL+1)*6)
1403
 
      INTEGER    NCOMB
1404
 
      PARAMETER (NCOMB=32)
1405
 
C     
1406
 
C     ARGUMENTS
1407
 
C     
1408
 
      INTEGER LEG_ID
1409
 
      INTEGER LEG_POLARIZATIONS(0:5)
1410
 
C     
1411
 
C     LOCALS
1412
 
C     
1413
 
      INTEGER I,J
1414
 
      LOGICAL ALL_SUMMED_OVER
1415
 
C     
1416
 
C     GLOBALS
1417
 
C     
1418
 
C     Entry 0 of the first dimension is all -1 if there is no
1419
 
C      polarization requirement.
1420
 
C     Then for each leg with ID legID, it is either summed over if
1421
 
C     POLARIZATIONS(legID,0) is -1, or the list of helicity considered
1422
 
C      for that
1423
 
C     leg is POLARIZATIONS(legID,1: POLARIZATIONS(legID,0)   ).
1424
 
      INTEGER POLARIZATIONS(0:NEXTERNAL,0:5)
1425
 
      DATA ((POLARIZATIONS(I,J),I=0,NEXTERNAL),J=0,5)/NPOLENTRIES*-1/
1426
 
      COMMON/ML5_0_BEAM_POL/POLARIZATIONS
1427
 
 
1428
 
      INTEGER BORN_POLARIZATIONS(0:NEXTERNAL,0:5)
1429
 
      COMMON/ML5_0_BORN_BEAM_POL/BORN_POLARIZATIONS
1430
 
 
1431
 
C     ----------
1432
 
C     BEGIN CODE
1433
 
C     ----------
1434
 
 
1435
 
      IF (LEG_POLARIZATIONS(0).EQ.-1) THEN
1436
 
        DO I=0,5
1437
 
          POLARIZATIONS(LEG_ID,I)=-1
1438
 
        ENDDO
1439
 
      ELSE
1440
 
        DO I=0,LEG_POLARIZATIONS(0)
1441
 
          POLARIZATIONS(LEG_ID,I)=LEG_POLARIZATIONS(I)
1442
 
        ENDDO
1443
 
        DO I=LEG_POLARIZATIONS(0)+1,5
1444
 
          POLARIZATIONS(LEG_ID,I)=-10000
1445
 
        ENDDO
1446
 
      ENDIF
1447
 
 
1448
 
      ALL_SUMMED_OVER = .TRUE.
1449
 
      DO I=1,NEXTERNAL
1450
 
        IF (POLARIZATIONS(I,0).NE.-1) THEN
1451
 
          ALL_SUMMED_OVER = .FALSE.
1452
 
          EXIT
1453
 
        ENDIF
1454
 
      ENDDO
1455
 
      IF (ALL_SUMMED_OVER) THEN
1456
 
        DO I=0,5
1457
 
          POLARIZATIONS(0,I)=-1
1458
 
        ENDDO
1459
 
      ELSE
1460
 
        DO I=0,5
1461
 
          POLARIZATIONS(0,I)=0
1462
 
        ENDDO
1463
 
      ENDIF
1464
 
 
1465
 
      DO I=0,NEXTERNAL
1466
 
        DO J=0,5
1467
 
          BORN_POLARIZATIONS(I,J) = POLARIZATIONS(I,J)
1468
 
        ENDDO
1469
 
      ENDDO
1470
 
 
1471
 
 
1472
 
      RETURN
1473
 
 
1474
 
      END
1475
 
 
1476
 
      SUBROUTINE ML5_0_SLOOPMATRIXHEL_THRES(P,HEL,ANS,PREC_ASKED
1477
 
     $ ,PREC_FOUND,RET_CODE)
1478
 
      IMPLICIT NONE
1479
 
C     
1480
 
C     CONSTANTS
1481
 
C     
1482
 
      INTEGER    NEXTERNAL
1483
 
      PARAMETER (NEXTERNAL=5)
1484
 
      INCLUDE 'nsquaredSO.inc'
1485
 
C     
1486
 
C     ARGUMENTS 
1487
 
C     
1488
 
      REAL*8 P(0:3,NEXTERNAL)
1489
 
      REAL*8 ANS(0:3,0:NSQUAREDSO)
1490
 
      INTEGER HEL,RET_CODE
1491
 
      REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
1492
 
C     
1493
 
C     GLOBAL VARIABLES
1494
 
C     
1495
 
      REAL*8 USER_STAB_PREC
1496
 
      COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC
1497
 
 
1498
 
      INTEGER I
1499
 
 
1500
 
      INTEGER H,T,U
1501
 
      REAL*8 ACCURACY(0:NSQUAREDSO)
1502
 
      COMMON/ML5_0_ACC/ACCURACY,H,T,U
1503
 
 
1504
 
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
1505
 
      COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
1506
 
 
1507
 
C     ----------
1508
 
C     BEGIN CODE
1509
 
C     ----------
1510
 
      USER_STAB_PREC = PREC_ASKED
1511
 
      CALL ML5_0_SLOOPMATRIXHEL(P,HEL,ANS)
1512
 
      IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY(0).LT.0.0D0))
1513
 
     $  THEN
1514
 
        BYPASS_CHECK = .TRUE.
1515
 
        CALL ML5_0_SLOOPMATRIXHEL(P,HEL,ANS)
1516
 
        BYPASS_CHECK = .FALSE.
1517
 
C       Make sure we correctly return an initialization-type T code
1518
 
        IF (T.EQ.2) T=4
1519
 
        IF (T.EQ.1) T=3
1520
 
      ENDIF
1521
 
 
1522
 
C     Reset it to default value not to affect next runs
1523
 
      USER_STAB_PREC = -1.0D0
1524
 
      DO I=0,NSQUAREDSO
1525
 
        PREC_FOUND(I)=ACCURACY(I)
1526
 
      ENDDO
1527
 
      RET_CODE=100*H+10*T+U
1528
 
 
1529
 
      END
1530
 
 
1531
 
      SUBROUTINE ML5_0_SLOOPMATRIX_THRES(P,ANS,PREC_ASKED,PREC_FOUND
1532
 
     $ ,RET_CODE)
1533
 
C     
1534
 
C     Inputs are:
1535
 
C     P(0:3, Nexternal)  double  :: Kinematic configuration
1536
 
C      (E,px,py,pz)
1537
 
C     PEC_ASKED          double  :: Target relative accuracy, -1 for
1538
 
C      default
1539
 
C     
1540
 
C     Outputs are:
1541
 
C     ANS(3)             double  :: Result (finite, single pole,
1542
 
C      double pole) 
1543
 
C     PREC_FOUND         double  :: Relative accuracy estimated for
1544
 
C      the result
1545
 
C     Returns -1 if no stab test could be performed.
1546
 
C     RET_CODE                   integer :: Return code. See below for details
1547
 
C     
1548
 
C     Return code conventions: RET_CODE = H*100 + T*10 + U
1549
 
C     
1550
 
C     H == 1
1551
 
C     Stability unknown.
1552
 
C     H == 2
1553
 
C     Stable PS (SPS) point.
1554
 
C     No stability rescue was necessary.
1555
 
C     H == 3
1556
 
C     Unstable PS (UPS) point.
1557
 
C     Stability rescue necessary, and successful.
1558
 
C     H == 4
1559
 
C     Exceptional PS (EPS) point.
1560
 
C     Stability rescue attempted, but unsuccessful.
1561
 
C     
1562
 
C     T == 1
1563
 
C     Default computation (double prec.) was performed.
1564
 
C     T == 2
1565
 
C     Quadruple precision was used for this PS point.
1566
 
C     T == 3
1567
 
C     MadLoop in initialization phase. Only double precision used.
1568
 
C     T == 4
1569
 
C     MadLoop in initialization phase. Quadruple precision used.
1570
 
C     
1571
 
C     U is a number left for future use (always set to 0 for now).
1572
 
C     example: TIR vs OPP usage.
1573
 
C     
1574
 
      IMPLICIT NONE
1575
 
C     
1576
 
C     CONSTANTS
1577
 
C     
1578
 
      INTEGER    NEXTERNAL
1579
 
      PARAMETER (NEXTERNAL=5)
1580
 
      INCLUDE 'nsquaredSO.inc'
1581
 
C     
1582
 
C     ARGUMENTS 
1583
 
C     
1584
 
      REAL*8 P(0:3,NEXTERNAL)
1585
 
      REAL*8 ANS(0:3,0:NSQUAREDSO)
1586
 
      REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
1587
 
      INTEGER RET_CODE
1588
 
C     
1589
 
C     GLOBAL VARIABLES
1590
 
C     
1591
 
      REAL*8 USER_STAB_PREC
1592
 
      COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC
1593
 
 
1594
 
      INTEGER I
1595
 
 
1596
 
      INTEGER H,T,U
1597
 
      REAL*8 ACCURACY(0:NSQUAREDSO)
1598
 
      COMMON/ML5_0_ACC/ACCURACY,H,T,U
1599
 
 
1600
 
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
1601
 
      COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
1602
 
 
1603
 
C     ----------
1604
 
C     BEGIN CODE
1605
 
C     ----------
1606
 
      USER_STAB_PREC = PREC_ASKED
1607
 
      CALL ML5_0_SLOOPMATRIX(P,ANS)
1608
 
      IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY(0).LT.0.0D0))
1609
 
     $  THEN
1610
 
        BYPASS_CHECK = .TRUE.
1611
 
        CALL ML5_0_SLOOPMATRIX(P,ANS)
1612
 
        BYPASS_CHECK = .FALSE.
1613
 
C       Make sure we correctly return an initialization-type T code
1614
 
        IF (T.EQ.2) T=4
1615
 
        IF (T.EQ.1) T=3
1616
 
      ENDIF
1617
 
 
1618
 
C     Reset it to default value not to affect next runs
1619
 
      USER_STAB_PREC = -1.0D0
1620
 
      DO I=0,NSQUAREDSO
1621
 
        PREC_FOUND(I)=ACCURACY(I)
1622
 
      ENDDO
1623
 
      RET_CODE=100*H+10*T+U
1624
 
 
1625
 
      END
1626