← Back to branch summary

~maddevelopers/mg5amcnlo/MS_with_decayevents

« back to all changes in this revision

Viewing changes to tests/input_files/IOTestsComparison/IOTestsComparison/short_ML_SMQCD_LoopInduced/gg_hh/loop_matrix.f

  • Committer: olivier Mattelaer
  • Date: 2014-03-31 16:58:22 UTC
  • mfrom: (252.1.59 2.1.1)
  • Revision ID: olivier.mattelaer@uclouvain.be-20140331165822-4eb1040586ic1y1v
pass to 2.1.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/input_files/IOTestsComparison/IOTestsComparison/short_ML_SMQCD_LoopInduced/gg_hh/loop_matrix.f
 
1
      SUBROUTINE SLOOPMATRIXHEL(P,HEL,ANS)
 
2
      IMPLICIT NONE
 
3
C     
 
4
C     CONSTANTS
 
5
C     
 
6
      INTEGER    NEXTERNAL
 
7
      PARAMETER (NEXTERNAL=4)
 
8
C     
 
9
C     ARGUMENTS 
 
10
C     
 
11
      REAL*8 P(0:3,NEXTERNAL)
 
12
      REAL*8 ANS(3)
 
13
      INTEGER HEL, USERHEL
 
14
      COMMON/USERCHOICE/USERHEL
 
15
C     ----------
 
16
C     BEGIN CODE
 
17
C     ----------
 
18
      USERHEL=HEL
 
19
      CALL SLOOPMATRIX(P,ANS)
 
20
      END
 
21
 
 
22
      LOGICAL FUNCTION ISZERO(TOTEST, REFERENCE_VALUE, AMPLN)
 
23
      IMPLICIT NONE
 
24
C     
 
25
C     CONSTANTS
 
26
C     
 
27
      INTEGER    NLOOPAMPS
 
28
      PARAMETER (NLOOPAMPS=20)
 
29
C     
 
30
C     ARGUMENTS 
 
31
C     
 
32
      REAL*8 TOTEST, REFERENCE_VALUE
 
33
      INTEGER AMPLN
 
34
C     
 
35
C     GLOBAL 
 
36
C     
 
37
      INCLUDE 'MadLoopParams.inc'
 
38
 
 
39
      COMPLEX*16 AMPL(3,NLOOPAMPS)
 
40
      LOGICAL S(NLOOPAMPS)
 
41
      COMMON/AMPL/AMPL,S
 
42
C     ----------
 
43
C     BEGIN CODE
 
44
C     ----------
 
45
      IF(ABS(REFERENCE_VALUE).EQ.0.0D0) THEN
 
46
        ISZERO=.FALSE.
 
47
        WRITE(*,*) '##E02 ERRROR Reference value for comparison is
 
48
     $    zero.'
 
49
        STOP
 
50
      ELSE
 
51
        ISZERO=((ABS(TOTEST)/ABS(REFERENCE_VALUE)).LT.ZEROTHRES)
 
52
      ENDIF
 
53
      IF(AMPLN.NE.-1) THEN
 
54
        IF((.NOT.ISZERO).AND.(.NOT.S(AMPLN))) THEN
 
55
          WRITE(*,*) '##W01 WARNING Contribution ',AMPLN,' is detected
 
56
     $      as contributing with CR=',(ABS(TOTEST)/ABS(REFERENCE_VALUE
 
57
     $     )),' but is unstable.'
 
58
        ENDIF
 
59
      ENDIF
 
60
 
 
61
      END
 
62
 
 
63
      SUBROUTINE SLOOPMATRIX(P_USER,ANS)
 
64
C     
 
65
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
 
66
C     By the MadGraph5_aMC@NLO Development Team
 
67
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
 
68
C     
 
69
C     Returns amplitude squared summed/avg over colors
 
70
C     and helicities for the point in phase space P(0:3,NEXTERNAL)
 
71
C     and external lines W(0:6,NEXTERNAL)
 
72
C     
 
73
C     Process: g g > h h QED=2 QCD=2 [ virt = QCD ] WEIGHTED=12
 
74
C     
 
75
      IMPLICIT NONE
 
76
C     
 
77
C     CONSTANTS
 
78
C     
 
79
      CHARACTER*64 PARAMFILENAME
 
80
      PARAMETER ( PARAMFILENAME='MadLoopParams.dat')
 
81
 
 
82
      INTEGER    NLOOPAMPS, NCTAMPS
 
83
      PARAMETER (NLOOPAMPS=20, NCTAMPS=4)
 
84
      INTEGER    NEXTERNAL
 
85
      PARAMETER (NEXTERNAL=4)
 
86
      INTEGER    NWAVEFUNCS
 
87
      PARAMETER (NWAVEFUNCS=5)
 
88
      INTEGER    NCOMB
 
89
      PARAMETER (NCOMB=4)
 
90
      REAL*8     ZERO
 
91
      PARAMETER (ZERO=0D0)
 
92
      REAL*16     MP__ZERO
 
93
      PARAMETER (MP__ZERO=0E0_16)
 
94
      COMPLEX*16 IMAG1
 
95
      PARAMETER (IMAG1=(0D0,1D0))
 
96
C     This parameter is designed for the check timing command of MG5
 
97
      LOGICAL SKIPLOOPEVAL
 
98
      PARAMETER (SKIPLOOPEVAL=.FALSE.)
 
99
      LOGICAL BOOTANDSTOP
 
100
      PARAMETER (BOOTANDSTOP=.FALSE.)
 
101
      INTEGER MAXSTABILITYLENGTH
 
102
      DATA MAXSTABILITYLENGTH/20/
 
103
      COMMON/STABILITY_TESTS/MAXSTABILITYLENGTH
 
104
C     
 
105
C     ARGUMENTS 
 
106
C     
 
107
      REAL*8 P_USER(0:3,NEXTERNAL)
 
108
      REAL*8 ANS(3)
 
109
C     
 
110
C     LOCAL VARIABLES 
 
111
C     
 
112
      INTEGER I,J,K,H
 
113
      INTEGER HELPICKED_BU, CTMODEINIT_BU
 
114
      REAL*8 MLSTABTHRES_BU
 
115
C     P is the actual PS POINT used for the computation, and can be
 
116
C      rotated for the stability test purposes.
 
117
      REAL*8 P(0:3,NEXTERNAL)
 
118
C     DP_RES STORES THE DOUBLE PRECISION RESULT OBTAINED FROM
 
119
C      DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY.
 
120
C     THE STAB_STAGE COUNTER I CORRESPONDANCE GOES AS FOLLOWS
 
121
C     I=1 -> ORIGINAL PS, CTMODE=1
 
122
C     I=2 -> ORIGINAL PS, CTMODE=2, (ONLY WITH CTMODERUN=-1)
 
123
C     I=3 -> PS WITH ROTATION 1, CTMODE=1, (ONLY WITH CTMODERUN=-2)
 
124
C     I=4 -> PS WITH ROTATION 2, CTMODE=1, (ONLY WITH CTMODERUN=-3)
 
125
C     I=5 -> POSSIBLY MORE EVALUATION METHODS IN THE FUTURE, MAX IS
 
126
C      MAXSTABILITYLENGTH
 
127
C     IF UNSTABLE IT GOES TO THE SAME PATTERN BUT STAB_INDEX IS THEN
 
128
C      I+20.
 
129
      LOGICAL EVAL_DONE(MAXSTABILITYLENGTH)
 
130
      LOGICAL DOING_QP_EVALS
 
131
      INTEGER STAB_INDEX,BASIC_CT_MODE
 
132
      INTEGER N_DP_EVAL, N_QP_EVAL
 
133
      DATA N_DP_EVAL/1/
 
134
      DATA N_QP_EVAL/1/
 
135
C     This is used for loop-induced where the reference scale for
 
136
C      comparisons is infered from
 
137
C     the previous points
 
138
      REAL*8 NEXTREF
 
139
      DATA NEXTREF/ZERO/
 
140
      INTEGER NPSPOINTS
 
141
      DATA NPSPOINTS/0/
 
142
 
 
143
      REAL*8 ACC
 
144
      REAL*8 DP_RES(3,MAXSTABILITYLENGTH)
 
145
C     QP_RES STORES THE QUADRUPLE PRECISION RESULT OBTAINED FROM
 
146
C      DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY.
 
147
      REAL*8 QP_RES(3,MAXSTABILITYLENGTH)
 
148
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
 
149
      INTEGER NATTEMPTS
 
150
      DATA NATTEMPTS/0/
 
151
      DATA IC/NEXTERNAL*1/
 
152
      REAL*8 BUFFR(3),TEMP(3),TEMP1,TEMP2
 
153
      COMPLEX*16 CFTOT
 
154
      LOGICAL FOUNDHELFILTER,FOUNDLOOPFILTER
 
155
      DATA FOUNDHELFILTER/.TRUE./
 
156
      DATA FOUNDLOOPFILTER/.TRUE./
 
157
      INTEGER IDEN
 
158
      DATA IDEN/512/
 
159
      INTEGER HELAVGFACTOR
 
160
      DATA HELAVGFACTOR/4/
 
161
      LOGICAL DONEHELDOUBLECHECK
 
162
      DATA DONEHELDOUBLECHECK/.FALSE./
 
163
      INTEGER NEPS
 
164
      DATA NEPS/0/
 
165
C     Below are variables to bypass the checkphase and insure
 
166
C      stability check to take place
 
167
      LOGICAL OLD_CHECKPHASE, OLD_HELDOUBLECHECKED
 
168
      LOGICAL OLD_GOODHEL(NCOMB)
 
169
      LOGICAL OLD_GOODAMP(NLOOPAMPS,NCOMB)
 
170
 
 
171
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
 
172
      COMMON/BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
 
173
C     
 
174
C     FUNCTIONS
 
175
C     
 
176
      LOGICAL ISZERO
 
177
C     
 
178
C     GLOBAL VARIABLES
 
179
C     
 
180
      INCLUDE 'coupl.inc'
 
181
      INCLUDE 'mp_coupl.inc'
 
182
      INCLUDE 'MadLoopParams.inc'
 
183
 
 
184
      INTEGER NTRY
 
185
      DATA NTRY/0/
 
186
      LOGICAL CHECKPHASE
 
187
      DATA CHECKPHASE/.TRUE./
 
188
      LOGICAL HELDOUBLECHECKED
 
189
      DATA HELDOUBLECHECKED/.FALSE./
 
190
      REAL*8 REF
 
191
      DATA REF/0.0D0/
 
192
      COMMON/INIT/NTRY,CHECKPHASE,HELDOUBLECHECKED,REF
 
193
 
 
194
C     THE LOGICAL BELOWS ARE JUST TO KEEP TRACK OF WHETHER THE MP_PS
 
195
C      HAS BEEN SET YET OR NOT AND WHETER THE MP EXTERNAL WFS HAVE
 
196
C      BEEN COMPUTED YET.
 
197
      LOGICAL MP_DONE
 
198
      DATA MP_DONE/.FALSE./
 
199
      COMMON/MP_DONE/MP_DONE
 
200
      LOGICAL MP_PS_SET
 
201
      DATA MP_PS_SET/.FALSE./
 
202
      COMMON/MP_PS_SET/MP_PS_SET
 
203
 
 
204
C     PS CAN POSSIBILY BE PASSED THROUGH IMPROVE_PS BUT IS NOT
 
205
C      MODIFIED FOR THE PURPOSE OF THE STABILITY TEST     
 
206
C     EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT
 
207
C      USED ANYWHERE ELSE
 
208
      REAL*8 PS(0:3,NEXTERNAL)
 
209
      COMMON/PSPOINT/PS
 
210
C     AGAIN BELOW, MP_PS IS THE FIXED (POSSIBLY IMPROVED) MP PS POINT
 
211
C      AND MP_P IS THE ONE WHICH CAN BE MODIFIED (I.E. ROTATED ETC.)
 
212
C      FOR STABILITY PURPOSE
 
213
C     EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT
 
214
C      USED ANYWHERE ELSE THAN HERE AND SET_MP_PS()
 
215
      REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL)
 
216
      COMMON/MP_PSPOINT/MP_PS,MP_P
 
217
 
 
218
      REAL*8 LSCALE
 
219
      INTEGER CTMODE
 
220
      COMMON/CT/LSCALE,CTMODE
 
221
 
 
222
      LOGICAL GOODHEL(NCOMB)
 
223
      LOGICAL GOODAMP(NLOOPAMPS,NCOMB)
 
224
      COMMON/FILTERS/GOODAMP,GOODHEL
 
225
 
 
226
      INTEGER HELPICKED
 
227
      DATA HELPICKED/-1/
 
228
      COMMON/HELCHOICE/HELPICKED
 
229
      INTEGER USERHEL
 
230
      DATA USERHEL/-1/
 
231
      COMMON/USERCHOICE/USERHEL
 
232
 
 
233
 
 
234
      COMPLEX*16 W(20,NWAVEFUNCS)
 
235
      INTEGER VALIDH
 
236
      COMMON/WFCTS/W
 
237
      COMMON/VALIDH/VALIDH
 
238
 
 
239
      COMPLEX*16 AMPL(3,NLOOPAMPS)
 
240
      LOGICAL S(NLOOPAMPS)
 
241
      COMMON/AMPL/AMPL,S
 
242
 
 
243
      INTEGER CF_D(NLOOPAMPS,NLOOPAMPS)
 
244
      INTEGER CF_N(NLOOPAMPS,NLOOPAMPS)
 
245
      COMMON/CF/CF_D,CF_N
 
246
 
 
247
      INTEGER HELC(NEXTERNAL,NCOMB)
 
248
      COMMON/HELCONFIGS/HELC
 
249
 
 
250
      REAL*8 PREC,USER_STAB_PREC
 
251
      DATA USER_STAB_PREC/-1.0D0/
 
252
      COMMON/USER_STAB_PREC/USER_STAB_PREC
 
253
 
 
254
C     Return codes H,T,U correspond to the hundreds, tens and units
 
255
C     building returncode, i.e.
 
256
C     RETURNCODE=100*RET_CODE_H+10*RET_CODE_T+RET_CODE_U
 
257
 
 
258
      INTEGER RET_CODE_H,RET_CODE_T,RET_CODE_U
 
259
      REAL*8 ACCURACY
 
260
      DATA ACCURACY/1.0D0/
 
261
      DATA RET_CODE_H,RET_CODE_T,RET_CODE_U/1,1,0/
 
262
      COMMON/ACC/ACCURACY,RET_CODE_H,RET_CODE_T,RET_CODE_U
 
263
 
 
264
      LOGICAL MP_DONE_ONCE
 
265
      DATA MP_DONE_ONCE/.FALSE./
 
266
      COMMON/MP_DONE_ONCE/MP_DONE_ONCE
 
267
 
 
268
C     ----------
 
269
C     BEGIN CODE
 
270
C     ----------
 
271
 
 
272
      IF(NTRY.EQ.0) THEN
 
273
        CALL MADLOOPPARAMREADER(PARAMFILENAME,.TRUE.)
 
274
        CALL SET_N_EVALS(N_DP_EVAL,N_QP_EVAL)
 
275
        HELDOUBLECHECKED=.NOT.DOUBLECHECKHELICITYFILTER
 
276
        DO J=1,NCOMB
 
277
          DO I=1,NCTAMPS
 
278
            GOODAMP(I,J)=.TRUE.
 
279
          ENDDO
 
280
        ENDDO
 
281
        OPEN(1, FILE='LoopFilter.dat', ERR=100, STATUS='OLD',         
 
282
     $     ACTION='READ')
 
283
        DO J=1,NCOMB
 
284
          READ(1,*,END=101) (GOODAMP(I,J),I=NCTAMPS+1,NLOOPAMPS)
 
285
        ENDDO
 
286
        GOTO 101
 
287
 100    CONTINUE
 
288
        FOUNDLOOPFILTER=.FALSE.
 
289
        DO J=1,NCOMB
 
290
          DO I=NCTAMPS+1,NLOOPAMPS
 
291
            GOODAMP(I,J)=(.NOT.USELOOPFILTER)
 
292
          ENDDO
 
293
        ENDDO
 
294
 101    CONTINUE
 
295
        CLOSE(1)
 
296
        OPEN(1, FILE='HelFilter.dat', ERR=102, STATUS='OLD',          
 
297
     $    ACTION='READ')
 
298
        READ(1,*,END=103) (GOODHEL(I),I=1,NCOMB)
 
299
        GOTO 103
 
300
 102    CONTINUE
 
301
        FOUNDHELFILTER=.FALSE.
 
302
        DO J=1,NCOMB
 
303
          GOODHEL(J)=.TRUE.
 
304
        ENDDO
 
305
 103    CONTINUE
 
306
        CLOSE(1)
 
307
        OPEN(1, FILE='ColorNumFactors.dat', ERR=104, STATUS='OLD'
 
308
     $   ,           ACTION='READ')
 
309
        DO I=1,NLOOPAMPS
 
310
          READ(1,*,END=105) (CF_N(I,J),J=1,NLOOPAMPS)
 
311
        ENDDO
 
312
        GOTO 105
 
313
 104    CONTINUE
 
314
        STOP 'Color factors could not be initialized from file
 
315
     $    ColorNumFactors.dat. File not found'
 
316
 105    CONTINUE
 
317
        CLOSE(1)
 
318
        OPEN(1, FILE='ColorDenomFactors.dat', ERR=106, STATUS='OLD'
 
319
     $   ,           ACTION='READ')
 
320
        DO I=1,NLOOPAMPS
 
321
          READ(1,*,END=107) (CF_D(I,J),J=1,NLOOPAMPS)
 
322
        ENDDO
 
323
        GOTO 107
 
324
 106    CONTINUE
 
325
        STOP 'Color factors could not be initialized from file
 
326
     $    ColorDenomFactors.dat. File not found'
 
327
 107    CONTINUE
 
328
        CLOSE(1)
 
329
        OPEN(1, FILE='HelConfigs.dat', ERR=108, STATUS='OLD',         
 
330
     $            ACTION='READ')
 
331
        DO H=1,NCOMB
 
332
          READ(1,*,END=109) (HELC(I,H),I=1,NEXTERNAL)
 
333
        ENDDO
 
334
        GOTO 109
 
335
 108    CONTINUE
 
336
        STOP 'Color helictiy configurations could not be initialized
 
337
     $    from file HelConfigs.dat. File not found'
 
338
 109    CONTINUE
 
339
        CLOSE(1)
 
340
        IF(BOOTANDSTOP) THEN
 
341
          WRITE(*,*) 'Stopped by user request.'
 
342
          STOP
 
343
        ENDIF
 
344
      ENDIF
 
345
 
 
346
      MP_DONE=.FALSE.
 
347
      MP_PS_SET=.FALSE.
 
348
      STAB_INDEX=0
 
349
      DOING_QP_EVALS=.FALSE.
 
350
      EVAL_DONE(1)=.TRUE.
 
351
      DO I=2,MAXSTABILITYLENGTH
 
352
        EVAL_DONE(I)=.FALSE.
 
353
      ENDDO
 
354
 
 
355
      IF (USER_STAB_PREC.GT.0.0D0) THEN
 
356
        MLSTABTHRES_BU=MLSTABTHRES
 
357
        MLSTABTHRES=USER_STAB_PREC
 
358
C       In the initialization, I cannot perform stability test and
 
359
C        therefore guarantee any precision
 
360
        CTMODEINIT_BU=CTMODEINIT
 
361
C       So either one choses quad precision directly
 
362
C       CTMODEINIT=4
 
363
C       Or, because this is very slow, we keep the orignal value. The
 
364
C        accuracy returned is -1 and tells the MC that he should not
 
365
C        trust the evaluation for checks.
 
366
        CTMODEINIT=CTMODEINIT_BU
 
367
      ENDIF
 
368
 
 
369
      IF(.NOT.BYPASS_CHECK) THEN
 
370
        NTRY=NTRY+1
 
371
      ENDIF
 
372
 
 
373
      IF(DONEHELDOUBLECHECK.AND.(.NOT.HELDOUBLECHECKED)) THEN
 
374
        HELDOUBLECHECKED=.TRUE.
 
375
        DONEHELDOUBLECHECK=.FALSE.
 
376
      ENDIF
 
377
 
 
378
      CHECKPHASE=(NTRY.LE.CHECKCYCLE).AND.(((.NOT.FOUNDLOOPFILTER
 
379
     $ ).AND.USELOOPFILTER).OR.(.NOT.FOUNDHELFILTER))
 
380
 
 
381
      IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDHELFILTER)) THEN
 
382
        OPEN(1, FILE='HelFilter.dat', ERR=110, STATUS='NEW',      
 
383
     $    ACTION='WRITE')
 
384
        WRITE(1,*) (GOODHEL(I),I=1,NCOMB)
 
385
 110    CONTINUE
 
386
        CLOSE(1)
 
387
        FOUNDHELFILTER=.TRUE.
 
388
      ENDIF
 
389
 
 
390
      IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDLOOPFILTER).AND.USELOOPFILT
 
391
     $ ER) THEN
 
392
        OPEN(1, FILE='LoopFilter.dat', ERR=111, STATUS='NEW',      
 
393
     $    ACTION='WRITE')
 
394
        DO J=1,NCOMB
 
395
          WRITE(1,*) (GOODAMP(I,J),I=NCTAMPS+1,NLOOPAMPS)
 
396
        ENDDO
 
397
 111    CONTINUE
 
398
        CLOSE(1)
 
399
        FOUNDLOOPFILTER=.TRUE.
 
400
      ENDIF
 
401
 
 
402
      IF (BYPASS_CHECK) THEN
 
403
        OLD_CHECKPHASE = CHECKPHASE
 
404
        OLD_HELDOUBLECHECKED = HELDOUBLECHECKED
 
405
        CHECKPHASE = .FALSE.
 
406
        HELDOUBLECHECKED = .TRUE.
 
407
        DO I=1,NCOMB
 
408
          OLD_GOODHEL(I)=GOODHEL(I)
 
409
          GOODHEL(I) = .TRUE.
 
410
        ENDDO
 
411
        DO I=1,NCOMB
 
412
          DO J=1,NLOOPAMPS
 
413
            OLD_GOODAMP(J,I)=GOODAMP(J,I)
 
414
            GOODAMP(J,I) = .TRUE.
 
415
          ENDDO
 
416
        ENDDO
 
417
      ENDIF
 
418
 
 
419
      IF(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
 
420
        HELPICKED=1
 
421
        CTMODE=CTMODEINIT
 
422
      ELSE
 
423
        IF (USERHEL.NE.-1) THEN
 
424
          IF(.NOT.GOODHEL(USERHEL)) THEN
 
425
            ANS(1)=0.0D0
 
426
            ANS(2)=0.0D0
 
427
            ANS(3)=0.0D0
 
428
            GOTO 9999
 
429
          ENDIF
 
430
        ENDIF
 
431
        HELPICKED=USERHEL
 
432
        IF (CTMODERUN.GT.-1) THEN
 
433
          CTMODE=CTMODERUN
 
434
        ELSE
 
435
          CTMODE=1
 
436
        ENDIF
 
437
      ENDIF
 
438
 
 
439
      DO I=1,NEXTERNAL
 
440
        DO J=0,3
 
441
          PS(J,I)=P_USER(J,I)
 
442
        ENDDO
 
443
      ENDDO
 
444
 
 
445
      IF (IMPROVEPSPOINT.GE.0) THEN
 
446
C       Make the input PS more precise (exact onshell and energy-moment
 
447
C       um conservation)
 
448
        CALL IMPROVE_PS_POINT_PRECISION(PS)
 
449
      ENDIF
 
450
 
 
451
      DO I=1,NEXTERNAL
 
452
        DO J=0,3
 
453
          P(J,I)=PS(J,I)
 
454
        ENDDO
 
455
      ENDDO
 
456
 
 
457
      DO K=1, 3
 
458
        BUFFR(K)=0.0D0
 
459
        DO I=1,NLOOPAMPS
 
460
          AMPL(K,I)=(0.0D0,0.0D0)
 
461
        ENDDO
 
462
      ENDDO
 
463
 
 
464
      LSCALE=DSQRT(ABS((P(0,1)+P(0,2))**2-(P(1,1)+P(1,2))**2-(P(2,1)
 
465
     $ +P(2,2))**2-(P(3,1)+P(3,2))**2))
 
466
 
 
467
C     For loop-induced, the reference for comparison is set later from
 
468
C      the total contribution of the previous PS point considered.
 
469
C     But you can edit here the value to be used for the first PS
 
470
C      point.
 
471
      IF (NPSPOINTS.EQ.0) THEN
 
472
        REF=1.0D-50
 
473
      ELSE
 
474
        REF=NEXTREF/DBLE(NPSPOINTS)
 
475
      ENDIF
 
476
 
 
477
 200  CONTINUE
 
478
 
 
479
      IF (CTMODE.EQ.0.OR.CTMODE.GE.4) THEN
 
480
        CALL MP_UPDATE_AS_PARAM()
 
481
      ENDIF
 
482
 
 
483
      IF (.NOT.MP_PS_SET.AND.(CTMODE.EQ.0.OR.CTMODE.GE.4)) THEN
 
484
        CALL SET_MP_PS(P_USER)
 
485
        MP_PS_SET = .TRUE.
 
486
      ENDIF
 
487
 
 
488
      DO K=1,3
 
489
        ANS(K)=0.0D0
 
490
      ENDDO
 
491
 
 
492
      VALIDH=-1
 
493
      DO H=1,NCOMB
 
494
        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(
 
495
     $   .NOT.HELDOUBLECHECKED).OR.GOODHEL(H)))) THEN
 
496
          IF (VALIDH.EQ.-1) VALIDH=H
 
497
          DO I=1,NEXTERNAL
 
498
            NHEL(I)=HELC(I,H)
 
499
          ENDDO
 
500
C         Check if we are in multiple precision and compute wfs and
 
501
C          amps accordingly if needed   
 
502
          IF (CTMODE.GE.4) THEN
 
503
C           Force that only current helicity is used in the routine
 
504
C            below
 
505
C           This should always be done, even if MP_DONE is True
 
506
C           because the AMPL of the R2 MUST be recomputed for loop
 
507
C            induced.
 
508
C           (because they are not saved for each hel configuration)
 
509
C           (This is not optimal unlike what is done int the loop
 
510
C            optimized output)
 
511
            HELPICKED_BU = HELPICKED
 
512
            HELPICKED = H
 
513
            CALL MP_BORN_AMPS_AND_WFS(MP_P)
 
514
            HELPICKED = HELPICKED_BU
 
515
            GOTO 300
 
516
          ENDIF
 
517
          CALL VXXXXX(P(0,1),ZERO,NHEL(1),-1*IC(1),W(1,1))
 
518
          CALL VXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2))
 
519
          CALL SXXXXX(P(0,3),+1*IC(3),W(1,3))
 
520
          CALL SXXXXX(P(0,4),+1*IC(4),W(1,4))
 
521
C         Counter-term amplitude(s) for loop diagram number 1
 
522
          CALL R2_GGHH_0(W(1,1),W(1,2),W(1,4),W(1,3),R2_GGHHB,AMPL(1
 
523
     $     ,1))
 
524
          CALL SSS1_1(W(1,3),W(1,4),GC_30,MDL_MH,MDL_WH,W(1,5))
 
525
C         Counter-term amplitude(s) for loop diagram number 3
 
526
          CALL VVS1_0(W(1,1),W(1,2),W(1,5),R2_GGHB,AMPL(1,2))
 
527
C         Counter-term amplitude(s) for loop diagram number 9
 
528
          CALL R2_GGHH_0(W(1,1),W(1,2),W(1,4),W(1,3),R2_GGHHT,AMPL(1
 
529
     $     ,3))
 
530
C         Counter-term amplitude(s) for loop diagram number 11
 
531
          CALL VVS1_0(W(1,1),W(1,2),W(1,5),R2_GGHT,AMPL(1,4))
 
532
 300      CONTINUE
 
533
          HELPICKED_BU=HELPICKED
 
534
          HELPICKED=H
 
535
          MP_DONE=.FALSE.
 
536
          IF(SKIPLOOPEVAL) THEN
 
537
            GOTO 1227
 
538
          ENDIF
 
539
C         Loop amplitude for loop diagram with ID 1
 
540
          CALL LOOP_4_4(1,1,2,4,3,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB
 
541
     $     ,KIND=16),DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16)
 
542
     $     ,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16),DCMPLX(MDL_MB)
 
543
     $     ,CMPLX(MP__MDL_MB,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5
 
544
     $     ,GC_33,MP__GC_33,GC_33,MP__GC_33,4,1,5,AMPL(1,5),S(5))
 
545
C         Loop amplitude for loop diagram with ID 2
 
546
          CALL LOOP_4_4(1,1,2,3,4,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB
 
547
     $     ,KIND=16),DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16)
 
548
     $     ,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16),DCMPLX(MDL_MB)
 
549
     $     ,CMPLX(MP__MDL_MB,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5
 
550
     $     ,GC_33,MP__GC_33,GC_33,MP__GC_33,4,1,6,AMPL(1,6),S(6))
 
551
C         Loop amplitude for loop diagram with ID 3
 
552
          CALL LOOP_3_3(2,1,2,5,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB
 
553
     $     ,KIND=16),DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16)
 
554
     $     ,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16),GC_5,MP__GC_5
 
555
     $     ,GC_5,MP__GC_5,GC_33,MP__GC_33,3,1,7,AMPL(1,7),S(7))
 
556
C         Loop amplitude for loop diagram with ID 4
 
557
          CALL LOOP_3_3(3,1,2,5,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB
 
558
     $     ,KIND=16),DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16)
 
559
     $     ,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16),GC_5,MP__GC_5
 
560
     $     ,GC_5,MP__GC_5,GC_33,MP__GC_33,3,1,8,AMPL(1,8),S(8))
 
561
C         Loop amplitude for loop diagram with ID 5
 
562
          CALL LOOP_4_4(4,1,2,4,3,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB
 
563
     $     ,KIND=16),DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16)
 
564
     $     ,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16),DCMPLX(MDL_MB)
 
565
     $     ,CMPLX(MP__MDL_MB,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5
 
566
     $     ,GC_33,MP__GC_33,GC_33,MP__GC_33,4,1,9,AMPL(1,9),S(9))
 
567
C         Loop amplitude for loop diagram with ID 6
 
568
          CALL LOOP_4_4(5,1,3,2,4,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB
 
569
     $     ,KIND=16),DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16)
 
570
     $     ,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16),DCMPLX(MDL_MB)
 
571
     $     ,CMPLX(MP__MDL_MB,KIND=16),GC_5,MP__GC_5,GC_33,MP__GC_33
 
572
     $     ,GC_5,MP__GC_5,GC_33,MP__GC_33,4,1,10,AMPL(1,10),S(10))
 
573
C         Loop amplitude for loop diagram with ID 7
 
574
          CALL LOOP_4_4(4,1,2,3,4,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB
 
575
     $     ,KIND=16),DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16)
 
576
     $     ,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16),DCMPLX(MDL_MB)
 
577
     $     ,CMPLX(MP__MDL_MB,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5
 
578
     $     ,GC_33,MP__GC_33,GC_33,MP__GC_33,4,1,11,AMPL(1,11),S(11))
 
579
C         Loop amplitude for loop diagram with ID 8
 
580
          CALL LOOP_4_4(6,1,3,2,4,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB
 
581
     $     ,KIND=16),DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16)
 
582
     $     ,DCMPLX(MDL_MB),CMPLX(MP__MDL_MB,KIND=16),DCMPLX(MDL_MB)
 
583
     $     ,CMPLX(MP__MDL_MB,KIND=16),GC_5,MP__GC_5,GC_33,MP__GC_33
 
584
     $     ,GC_5,MP__GC_5,GC_33,MP__GC_33,4,1,12,AMPL(1,12),S(12))
 
585
C         Loop amplitude for loop diagram with ID 9
 
586
          CALL LOOP_4_4(1,1,2,4,3,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT
 
587
     $     ,KIND=16),DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16)
 
588
     $     ,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16),DCMPLX(MDL_MT)
 
589
     $     ,CMPLX(MP__MDL_MT,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5
 
590
     $     ,GC_37,MP__GC_37,GC_37,MP__GC_37,4,1,13,AMPL(1,13),S(13))
 
591
C         Loop amplitude for loop diagram with ID 10
 
592
          CALL LOOP_4_4(1,1,2,3,4,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT
 
593
     $     ,KIND=16),DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16)
 
594
     $     ,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16),DCMPLX(MDL_MT)
 
595
     $     ,CMPLX(MP__MDL_MT,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5
 
596
     $     ,GC_37,MP__GC_37,GC_37,MP__GC_37,4,1,14,AMPL(1,14),S(14))
 
597
C         Loop amplitude for loop diagram with ID 11
 
598
          CALL LOOP_3_3(2,1,2,5,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT
 
599
     $     ,KIND=16),DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16)
 
600
     $     ,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16),GC_5,MP__GC_5
 
601
     $     ,GC_5,MP__GC_5,GC_37,MP__GC_37,3,1,15,AMPL(1,15),S(15))
 
602
C         Loop amplitude for loop diagram with ID 12
 
603
          CALL LOOP_3_3(3,1,2,5,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT
 
604
     $     ,KIND=16),DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16)
 
605
     $     ,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16),GC_5,MP__GC_5
 
606
     $     ,GC_5,MP__GC_5,GC_37,MP__GC_37,3,1,16,AMPL(1,16),S(16))
 
607
C         Loop amplitude for loop diagram with ID 13
 
608
          CALL LOOP_4_4(4,1,2,4,3,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT
 
609
     $     ,KIND=16),DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16)
 
610
     $     ,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16),DCMPLX(MDL_MT)
 
611
     $     ,CMPLX(MP__MDL_MT,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5
 
612
     $     ,GC_37,MP__GC_37,GC_37,MP__GC_37,4,1,17,AMPL(1,17),S(17))
 
613
C         Loop amplitude for loop diagram with ID 14
 
614
          CALL LOOP_4_4(5,1,3,2,4,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT
 
615
     $     ,KIND=16),DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16)
 
616
     $     ,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16),DCMPLX(MDL_MT)
 
617
     $     ,CMPLX(MP__MDL_MT,KIND=16),GC_5,MP__GC_5,GC_37,MP__GC_37
 
618
     $     ,GC_5,MP__GC_5,GC_37,MP__GC_37,4,1,18,AMPL(1,18),S(18))
 
619
C         Loop amplitude for loop diagram with ID 15
 
620
          CALL LOOP_4_4(4,1,2,3,4,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT
 
621
     $     ,KIND=16),DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16)
 
622
     $     ,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16),DCMPLX(MDL_MT)
 
623
     $     ,CMPLX(MP__MDL_MT,KIND=16),GC_5,MP__GC_5,GC_5,MP__GC_5
 
624
     $     ,GC_37,MP__GC_37,GC_37,MP__GC_37,4,1,19,AMPL(1,19),S(19))
 
625
C         Loop amplitude for loop diagram with ID 16
 
626
          CALL LOOP_4_4(6,1,3,2,4,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT
 
627
     $     ,KIND=16),DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16)
 
628
     $     ,DCMPLX(MDL_MT),CMPLX(MP__MDL_MT,KIND=16),DCMPLX(MDL_MT)
 
629
     $     ,CMPLX(MP__MDL_MT,KIND=16),GC_5,MP__GC_5,GC_37,MP__GC_37
 
630
     $     ,GC_5,MP__GC_5,GC_37,MP__GC_37,4,1,20,AMPL(1,20),S(20))
 
631
          HELPICKED=HELPICKED_BU
 
632
          DO I=NCTAMPS+1,NLOOPAMPS
 
633
            IF((CTMODERUN.NE.-1).AND..NOT.CHECKPHASE.AND.(.NOT.S(I))
 
634
     $       ) THEN
 
635
              WRITE(*,*) '##W03 WARNING Contribution ',I
 
636
              WRITE(*,*) ' is unstable for helicity ',H
 
637
            ENDIF
 
638
C           IF(.NOT.ISZERO(ABS(AMPL(2,I))+ABS(AMPL(3,I)),REF,-1,H))
 
639
C            THEN
 
640
C           WRITE(*,*) '##W04 WARNING Contribution ',I,' for helicity
 
641
C            ',H,' has a contribution to the poles.'
 
642
C           WRITE(*,*) 'Finite contribution         = ',AMPL(1,I)
 
643
C           WRITE(*,*) 'single pole contribution    = ',AMPL(2,I)
 
644
C           WRITE(*,*) 'double pole contribution    = ',AMPL(3,I)
 
645
C           ENDIF
 
646
          ENDDO
 
647
 1227     CONTINUE
 
648
          DO I=1,NLOOPAMPS
 
649
            DO J=1,NLOOPAMPS
 
650
              CFTOT=DCMPLX(CF_N(I,J)/DBLE(ABS(CF_D(I,J))),0.0D0)
 
651
              IF(CF_D(I,J).LT.0) CFTOT=CFTOT*IMAG1
 
652
              ANS(1)=ANS(1)+DBLE(CFTOT*AMPL(1,I)*DCONJG(AMPL(1,J)))
 
653
              IF (J.EQ.1) THEN
 
654
                ANS(2)=ANS(2)+DBLE(CFTOT*AMPL(2,I))+DIMAG(CFTOT
 
655
     $           *AMPL(2,I))
 
656
                ANS(3)=ANS(3)+DBLE(CFTOT*AMPL(3,I))+DIMAG(CFTOT
 
657
     $           *AMPL(3,I))
 
658
              ENDIF
 
659
            ENDDO
 
660
          ENDDO
 
661
        ENDIF
 
662
      ENDDO
 
663
 
 
664
C     WHEN CTMODE IS >=4, then the MP computation of wfs and amps is
 
665
C      automatically done.
 
666
      IF (CTMODE.GE.4) THEN
 
667
        MP_DONE = .TRUE.
 
668
      ENDIF
 
669
 
 
670
      IF(SKIPLOOPEVAL) THEN
 
671
        GOTO 1226
 
672
      ENDIF
 
673
 
 
674
 
 
675
 
 
676
      IF(.NOT.ISZERO(ABS(ANS(2))+ABS(ANS(3)),REF*(10.0D0**-2),
 
677
     $ -1,H)) THEN
 
678
        WRITE(*,*) '##W05 WARNING Found a PS point with a contribution
 
679
     $    to the single pole.'
 
680
        WRITE(*,*) 'Finite contribution         = ',ANS(1)
 
681
        WRITE(*,*) 'single pole contribution    = ',ANS(2)
 
682
        WRITE(*,*) 'double pole contribution    = ',ANS(3)
 
683
      ENDIF
 
684
 
 
685
 1226 CONTINUE
 
686
 
 
687
      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN
 
688
C       Update of NEXTREF, will be used for loop induced only.
 
689
        NEXTREF = NEXTREF + ANS(1) + ANS(2) + ANS(3)
 
690
        IF((USERHEL.EQ.-1).OR.(USERHEL.EQ.HELPICKED)) THEN
 
691
          BUFFR(1)=BUFFR(1)+ANS(1)
 
692
          BUFFR(2)=BUFFR(2)+ANS(2)
 
693
          BUFFR(3)=BUFFR(3)+ANS(3)
 
694
        ENDIF
 
695
 
 
696
        IF (CHECKPHASE) THEN
 
697
C         SET THE HELICITY FILTER
 
698
          IF(.NOT.FOUNDHELFILTER) THEN
 
699
            IF(ISZERO(ABS(ANS(1))+ABS(ANS(2))+ABS(ANS(3)),REF
 
700
     $       /DBLE(NCOMB),-1)) THEN
 
701
              IF(NTRY.EQ.1) THEN
 
702
                GOODHEL(HELPICKED)=.FALSE.
 
703
              ELSEIF(GOODHEL(HELPICKED)) THEN
 
704
                WRITE(*,*) '##W02A WARNING Inconsistent helicity '
 
705
     $           ,HELPICKED
 
706
                IF(HELINITSTARTOVER) THEN
 
707
                  WRITE(*,*) '##I01 INFO Initialization starting over
 
708
     $              because of inconsistency in the helicity filter
 
709
     $              setup.'
 
710
                  NTRY=0
 
711
                ENDIF
 
712
              ENDIF
 
713
            ELSE
 
714
              IF(.NOT.GOODHEL(HELPICKED)) THEN
 
715
                WRITE(*,*) '##W02B WARNING Inconsistent helicity '
 
716
     $           ,HELPICKED
 
717
                IF(HELINITSTARTOVER) THEN
 
718
                  WRITE(*,*) '##I01 INFO Initialization starting over
 
719
     $              because of inconsistency in the helicity filter
 
720
     $              setup.'
 
721
                  NTRY=0
 
722
                ELSE
 
723
                  GOODHEL(HELPICKED)=.TRUE.
 
724
                ENDIF
 
725
              ENDIF
 
726
            ENDIF
 
727
          ENDIF
 
728
 
 
729
C         SET THE LOOP FILTER
 
730
          IF(.NOT.FOUNDLOOPFILTER.AND.USELOOPFILTER) THEN
 
731
            DO I=NCTAMPS+1,NLOOPAMPS
 
732
              IF(.NOT.ISZERO(ABS(AMPL(1,I))+ABS(AMPL(2,I))+ABS(AMPL(3
 
733
     $         ,I)),(REF*1.0D-4),I)) THEN
 
734
                IF(NTRY.EQ.1) THEN
 
735
                  GOODAMP(I,HELPICKED)=.TRUE.
 
736
                ELSEIF(.NOT.GOODAMP(I,HELPICKED)) THEN
 
737
                  WRITE(*,*) '##W02 WARNING Inconsistent loop amp ',I
 
738
     $             ,' for helicity ',HELPICKED,'.'
 
739
                  IF(LOOPINITSTARTOVER) THEN
 
740
                    WRITE(*,*) '##I01 INFO Initialization starting
 
741
     $                over because of inconsistency in the loop filter
 
742
     $                setup.'
 
743
                    NTRY=0
 
744
                  ELSE
 
745
                    GOODAMP(I,HELPICKED)=.TRUE.
 
746
                  ENDIF
 
747
                ENDIF
 
748
              ENDIF
 
749
            ENDDO
 
750
          ENDIF
 
751
        ELSEIF (.NOT.HELDOUBLECHECKED)THEN
 
752
          IF ((.NOT.GOODHEL(HELPICKED)).AND.(.NOT.ISZERO(ABS(ANS(1))
 
753
     $     +ABS(ANS(2))+ABS(ANS(3)),REF/DBLE(NCOMB),-1))) THEN
 
754
            WRITE(*,*) '##W15 Helicity filter could not be successfully
 
755
     $        double checked.'
 
756
            WRITE(*,*) 'One reason for this is that you have changed
 
757
     $        sensible parameters which affected what are the zero
 
758
     $        helicity configurations.'
 
759
            WRITE(*,*) 'MadLoop will try to reset the Helicity filter
 
760
     $        with the next PS points it receives.'
 
761
            NTRY=0
 
762
            OPEN(30,FILE='HelFilter.dat',ERR=349)
 
763
 349        CONTINUE
 
764
            CLOSE(30,STATUS='delete')
 
765
          ENDIF
 
766
C         SET HELDOUBLECHECKED TO .TRUE. WHEN DONE
 
767
C         even if it failed we do not want to redo the check afterwards
 
768
C          if HELINITSTARTOVER=.FALSE.
 
769
          IF (HELPICKED.EQ.NCOMB.AND.(NTRY.NE.0.OR..NOT.HELINITSTARTOVE
 
770
     $     R)) THEN
 
771
            DONEHELDOUBLECHECK=.TRUE.
 
772
          ENDIF
 
773
        ENDIF
 
774
 
 
775
C       GOTO NEXT HELICITY OR FINISH
 
776
        IF(HELPICKED.NE.NCOMB) THEN
 
777
          HELPICKED=HELPICKED+1
 
778
          MP_DONE=.FALSE.
 
779
          GOTO 200
 
780
        ELSE
 
781
          ANS(1)=BUFFR(1)
 
782
          ANS(2)=BUFFR(2)
 
783
          ANS(3)=BUFFR(3)
 
784
C         We add one here to the number of PS points used for building
 
785
C          the reference scale for comparison (used only for loop-induc
 
786
C         ed processes).
 
787
          NPSPOINTS = NPSPOINTS+1
 
788
          IF(NTRY.EQ.0) THEN
 
789
            NATTEMPTS=NATTEMPTS+1
 
790
            IF(NATTEMPTS.EQ.MAXATTEMPTS) THEN
 
791
              WRITE(*,*) '##E01 ERROR Could not initialize the filters
 
792
     $          in ',MAXATTEMPTS,' trials'
 
793
              STOP
 
794
            ENDIF
 
795
          ENDIF
 
796
        ENDIF
 
797
 
 
798
      ENDIF
 
799
 
 
800
      DO K=1,3
 
801
        ANS(K)=ANS(K)/DBLE(IDEN)
 
802
        IF (USERHEL.NE.-1) THEN
 
803
          ANS(K)=ANS(K)*HELAVGFACTOR
 
804
        ENDIF
 
805
      ENDDO
 
806
 
 
807
      IF(.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.(CTMODERUN.LE.
 
808
     $ -1)) THEN
 
809
        STAB_INDEX=STAB_INDEX+1
 
810
        IF(DOING_QP_EVALS) THEN
 
811
          QP_RES(1,STAB_INDEX)=ANS(1)
 
812
          QP_RES(2,STAB_INDEX)=ANS(2)
 
813
          QP_RES(3,STAB_INDEX)=ANS(3)
 
814
        ELSE
 
815
          DP_RES(1,STAB_INDEX)=ANS(1)
 
816
          DP_RES(2,STAB_INDEX)=ANS(2)
 
817
          DP_RES(3,STAB_INDEX)=ANS(3)
 
818
        ENDIF
 
819
 
 
820
        IF(DOING_QP_EVALS) THEN
 
821
          BASIC_CT_MODE=4
 
822
        ELSE
 
823
          BASIC_CT_MODE=1
 
824
        ENDIF
 
825
 
 
826
C       BEGINNING OF THE DEFINITIONS OF THE DIFFERENT EVALUATION
 
827
C        METHODS
 
828
 
 
829
        IF(.NOT.EVAL_DONE(2)) THEN
 
830
          EVAL_DONE(2)=.TRUE.
 
831
          CTMODE=BASIC_CT_MODE+1
 
832
          GOTO 200
 
833
        ENDIF
 
834
 
 
835
        CTMODE=BASIC_CT_MODE
 
836
 
 
837
        IF(.NOT.EVAL_DONE(3).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
 
838
     $   .1).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
 
839
          EVAL_DONE(3)=.TRUE.
 
840
          CALL ROTATE_PS(PS,P,1)
 
841
          IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,1)
 
842
          GOTO 200
 
843
        ENDIF
 
844
 
 
845
        IF(.NOT.EVAL_DONE(4).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
 
846
     $   .2).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
 
847
          EVAL_DONE(4)=.TRUE.
 
848
          CALL ROTATE_PS(PS,P,2)
 
849
          IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,2)
 
850
          GOTO 200
 
851
        ENDIF
 
852
 
 
853
        CALL ROTATE_PS(PS,P,0)
 
854
        IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,0)
 
855
 
 
856
C       END OF THE DEFINITIONS OF THE DIFFERENT EVALUATION METHODS
 
857
 
 
858
        IF(DOING_QP_EVALS) THEN
 
859
          CALL COMPUTE_ACCURACY(QP_RES,N_QP_EVAL,ACC,ANS)
 
860
          ACCURACY=ACC
 
861
          RET_CODE_H=3
 
862
          IF(ACC.GE.MLSTABTHRES) THEN
 
863
            RET_CODE_H=4
 
864
            NEPS=NEPS+1
 
865
            CALL COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,TEMP1,TEMP)
 
866
            WRITE(*,*) '##W03 WARNING An unstable PS point was'
 
867
     $       ,       ' detected.'
 
868
            WRITE(*,*) '(DP,QP) accuracies : (',TEMP1,',',ACC,')'
 
869
            WRITE(*,*) 'Best estimate (fin,1eps,2eps) :',(ANS(I),I=1,3)
 
870
            IF(NEPS.LE.10) THEN
 
871
              WRITE(*,*) 'Double precision evaluations :',(DP_RES(1,I)
 
872
     $         ,I=1,N_DP_EVAL)
 
873
              WRITE(*,*) 'Quad   precision evaluations :',(QP_RES(1,I)
 
874
     $         ,I=1,N_QP_EVAL)
 
875
              WRITE(*,*) 'PS point specification :'
 
876
              WRITE(*,*) 'Renormalization scale MU_R=',MU_R
 
877
              DO I=1,NEXTERNAL
 
878
                WRITE (*,'(i2,1x,4e27.17)') I, P(0,I),P(1,I),P(2,I)
 
879
     $           ,P(3,I)
 
880
              ENDDO
 
881
            ENDIF
 
882
            IF(NEPS.EQ.10) THEN
 
883
              WRITE(*,*) 'Further output of the details of these
 
884
     $          unstable PS points will now be suppressed.'
 
885
            ENDIF
 
886
          ENDIF
 
887
        ELSE
 
888
          CALL COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,ACC,ANS)
 
889
          IF(ACC.GE.MLSTABTHRES) THEN
 
890
            DOING_QP_EVALS=.TRUE.
 
891
            EVAL_DONE(1)=.TRUE.
 
892
            DO I=2,MAXSTABILITYLENGTH
 
893
              EVAL_DONE(I)=.FALSE.
 
894
            ENDDO
 
895
            STAB_INDEX=0
 
896
            CTMODE=4
 
897
            GOTO 200
 
898
          ELSE
 
899
            RET_CODE_H=2
 
900
            ACCURACY=ACC
 
901
          ENDIF
 
902
        ENDIF
 
903
      ELSE
 
904
        RET_CODE_H=1
 
905
        ACCURACY=-1.0D0
 
906
      ENDIF
 
907
 
 
908
 9999 CONTINUE
 
909
 
 
910
C     Finalize the return code
 
911
      IF (MP_DONE_ONCE) THEN
 
912
        RET_CODE_T=2
 
913
      ELSE
 
914
        RET_CODE_T=1
 
915
      ENDIF
 
916
      IF(CHECKPHASE.OR..NOT.HELDOUBLECHECKED) THEN
 
917
        RET_CODE_H=1
 
918
        RET_CODE_T=RET_CODE_T+2
 
919
        ACCURACY=-1.0D0
 
920
      ENDIF
 
921
      RET_CODE_U=0
 
922
 
 
923
C     Reinitialize the default threshold if it was specified by the
 
924
C      user
 
925
      IF (USER_STAB_PREC.GT.0.0D0) THEN
 
926
        MLSTABTHRES=MLSTABTHRES_BU
 
927
        CTMODEINIT=CTMODEINIT_BU
 
928
      ENDIF
 
929
 
 
930
C     Reinitialize the check phase logicals and the filters if check
 
931
C      bypassed
 
932
      IF (BYPASS_CHECK) THEN
 
933
        CHECKPHASE = OLD_CHECKPHASE
 
934
        HELDOUBLECHECKED = OLD_HELDOUBLECHECKED
 
935
        DO I=1,NCOMB
 
936
          GOODHEL(I)=OLD_GOODHEL(I)
 
937
        ENDDO
 
938
        DO I=1,NCOMB
 
939
          DO J=1,NLOOPAMPS
 
940
            GOODAMP(J,I)=OLD_GOODAMP(J,I)
 
941
          ENDDO
 
942
        ENDDO
 
943
      ENDIF
 
944
 
 
945
      END
 
946
 
 
947
      SUBROUTINE COMPUTE_ACCURACY(FULLLIST, LENGTH, ACC, ESTIMATE)
 
948
      IMPLICIT NONE
 
949
C     
 
950
C     PARAMETERS 
 
951
C     
 
952
      INTEGER MAXSTABILITYLENGTH
 
953
      COMMON/STABILITY_TESTS/MAXSTABILITYLENGTH
 
954
C     
 
955
C     ARGUMENTS 
 
956
C     
 
957
      REAL*8 FULLLIST(3,MAXSTABILITYLENGTH)
 
958
      INTEGER LENGTH
 
959
      REAL*8 ACC, ESTIMATE(3)
 
960
C     
 
961
C     LOCAL VARIABLES 
 
962
C     
 
963
      LOGICAL MASK(MAXSTABILITYLENGTH)
 
964
      LOGICAL MASK3(3)
 
965
      DATA MASK3/.TRUE.,.TRUE.,.TRUE./
 
966
      INTEGER I,J
 
967
      REAL*8 AVG
 
968
      REAL*8 DIFF
 
969
      REAL*8 ACCURACIES(3)
 
970
      REAL*8 LIST(MAXSTABILITYLENGTH)
 
971
 
 
972
C     ----------
 
973
C     BEGIN CODE
 
974
C     ----------
 
975
      DO I=1,LENGTH
 
976
        MASK(I)=.TRUE.
 
977
      ENDDO
 
978
      DO I=LENGTH+1,MAXSTABILITYLENGTH
 
979
        MASK(I)=.FALSE.
 
980
C       For some architectures, it is necessary to initialize all the
 
981
C        elements of fulllist(i,j)
 
982
C       Beware that if the length provided is incorrect, then this can
 
983
C        corrup the fulllist given in argument.
 
984
        DO J=1,3
 
985
          FULLLIST(J,I)=0.0D0
 
986
        ENDDO
 
987
      ENDDO
 
988
 
 
989
      DO I=1,3
 
990
        DO J=1,MAXSTABILITYLENGTH
 
991
          LIST(J)=FULLLIST(I,J)
 
992
        ENDDO
 
993
        DIFF=MAXVAL(LIST,1,MASK)-MINVAL(LIST,1,MASK)
 
994
        AVG=(MAXVAL(LIST,1,MASK)+MINVAL(LIST,1,MASK))/2.0D0
 
995
        ESTIMATE(I)=AVG
 
996
        IF (AVG.EQ.0.0D0) THEN
 
997
          ACCURACIES(I)=DIFF
 
998
        ELSE
 
999
          ACCURACIES(I)=DIFF/ABS(AVG)
 
1000
        ENDIF
 
1001
      ENDDO
 
1002
 
 
1003
C     The technique below is too sensitive, typically to
 
1004
C     unstablities in very small poles
 
1005
C     ACC=MAXVAL(ACCURACIES,1,MASK3)
 
1006
C     The following is used instead
 
1007
      ACC = 0.0D0
 
1008
      AVG = 0.0D0
 
1009
      DO I=1,3
 
1010
        ACC = ACC + ACCURACIES(I)*ABS(ESTIMATE(I))
 
1011
        AVG = AVG + ESTIMATE(I)
 
1012
      ENDDO
 
1013
      ACC  = ACC / ( ABS(AVG) / 3.0D0)
 
1014
 
 
1015
      END
 
1016
 
 
1017
      SUBROUTINE SET_N_EVALS(N_DP_EVALS,N_QP_EVALS)
 
1018
 
 
1019
      IMPLICIT NONE
 
1020
      INTEGER N_DP_EVALS, N_QP_EVALS
 
1021
 
 
1022
      INCLUDE 'MadLoopParams.inc'
 
1023
 
 
1024
      IF(CTMODERUN.LE.-1) THEN
 
1025
        N_DP_EVALS=2+NROTATIONS_DP
 
1026
        N_QP_EVALS=2+NROTATIONS_QP
 
1027
      ELSE
 
1028
        N_DP_EVALS=1
 
1029
        N_QP_EVALS=1
 
1030
      ENDIF
 
1031
 
 
1032
      IF(N_DP_EVALS.GT.20.OR.N_QP_EVALS.GT.20) THEN
 
1033
        WRITE(*,*) 'ERROR:: Increase hardcoded maxstabilitylength.'
 
1034
        STOP
 
1035
      ENDIF
 
1036
 
 
1037
      END
 
1038
 
 
1039
 
 
1040
C     THIS SUBROUTINE SIMPLY SET THE GLOBAL PS CONFIGURATION GLOBAL
 
1041
C      VARIABLES FROM A GIVEN VARIABLE IN DOUBLE PRECISION
 
1042
      SUBROUTINE SET_MP_PS(P)
 
1043
 
 
1044
      INTEGER    NEXTERNAL
 
1045
      PARAMETER (NEXTERNAL=4)
 
1046
      REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL)
 
1047
      COMMON/MP_PSPOINT/MP_PS,MP_P
 
1048
      REAL*8 P(0:3,NEXTERNAL)
 
1049
 
 
1050
      DO I=1,NEXTERNAL
 
1051
        DO J=0,3
 
1052
          MP_PS(J,I)=P(J,I)
 
1053
        ENDDO
 
1054
      ENDDO
 
1055
      CALL MP_IMPROVE_PS_POINT_PRECISION(MP_PS)
 
1056
      DO I=1,NEXTERNAL
 
1057
        DO J=0,3
 
1058
          MP_P(J,I)=MP_PS(J,I)
 
1059
        ENDDO
 
1060
      ENDDO
 
1061
 
 
1062
      END
 
1063
 
 
1064
      SUBROUTINE FORCE_STABILITY_CHECK(ONOFF)
 
1065
C     
 
1066
C     This function can be called by the MadLoop user so as to always
 
1067
C      have stability
 
1068
C     checked, even during initialisation, when calling the *_thres
 
1069
C      routines.
 
1070
C     
 
1071
      LOGICAL ONOFF
 
1072
 
 
1073
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
 
1074
      DATA BYPASS_CHECK, ALWAYS_TEST_STABILITY /.FALSE.,.FALSE./
 
1075
      COMMON/BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
 
1076
 
 
1077
      ALWAYS_TEST_STABILITY = ONOFF
 
1078
 
 
1079
      END
 
1080
 
 
1081
      SUBROUTINE SLOOPMATRIXHEL_THRES(P,HEL,ANS,PREC_ASKED,PREC_FOUND
 
1082
     $ ,RET_CODE)
 
1083
      IMPLICIT NONE
 
1084
C     
 
1085
C     CONSTANTS
 
1086
C     
 
1087
      INTEGER    NEXTERNAL
 
1088
      PARAMETER (NEXTERNAL=4)
 
1089
C     
 
1090
C     ARGUMENTS 
 
1091
C     
 
1092
      REAL*8 P(0:3,NEXTERNAL)
 
1093
      REAL*8 ANS(3)
 
1094
      INTEGER HEL,RET_CODE
 
1095
      REAL*8 PREC_ASKED,PREC_FOUND
 
1096
C     
 
1097
C     GLOBAL VARIABLES
 
1098
C     
 
1099
      REAL*8 USER_STAB_PREC
 
1100
      COMMON/USER_STAB_PREC/USER_STAB_PREC
 
1101
 
 
1102
      INTEGER H,T,U
 
1103
      REAL*8 ACCURACY
 
1104
      COMMON/ACC/ACCURACY,H,T,U
 
1105
 
 
1106
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
 
1107
      COMMON/BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
 
1108
 
 
1109
C     ----------
 
1110
C     BEGIN CODE
 
1111
C     ----------
 
1112
      USER_STAB_PREC = PREC_ASKED
 
1113
 
 
1114
      CALL SLOOPMATRIXHEL(P,HEL,ANS)
 
1115
      IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY.LT.0.0D0)) THEN
 
1116
        BYPASS_CHECK = .TRUE.
 
1117
        CALL SLOOPMATRIXHEL(P,HEL,ANS)
 
1118
        BYPASS_CHECK = .FALSE.
 
1119
C       Make sure we correctly return an initialization-type T code
 
1120
        IF (T.EQ.2) T=4
 
1121
        IF (T.EQ.1) T=3
 
1122
      ENDIF
 
1123
 
 
1124
      PREC_FOUND=ACCURACY
 
1125
      RET_CODE=100*H+10*T+U
 
1126
C     Reset it to default value not to affect next runs
 
1127
      USER_STAB_PREC = -1.0D0
 
1128
 
 
1129
      END
 
1130
 
 
1131
      SUBROUTINE SLOOPMATRIX_THRES(P,ANS,PREC_ASKED,PREC_FOUND
 
1132
     $ ,RET_CODE)
 
1133
C     
 
1134
C     Inputs are:
 
1135
C     P(0:3, Nexternal)  double  :: Kinematic configuration (E,px,py,pz
 
1136
C     )
 
1137
C     PEC_ASKED          double  :: Target relative accuracy, -1 for
 
1138
C      default
 
1139
C     
 
1140
C     Outputs are:
 
1141
C     ANS(3)             double  :: Result (finite, single pole,
 
1142
C      double pole) 
 
1143
C     PREC_FOUND         double  :: Relative accuracy estimated for
 
1144
C      the result
 
1145
C     Returns -1 if no stab test could be performed.
 
1146
C     RET_CODE                   integer :: Return code. See below for details
 
1147
C     
 
1148
C     Return code conventions: RET_CODE = H*100 + T*10 + U
 
1149
C     
 
1150
C     H == 1
 
1151
C     Stability unknown.
 
1152
C     H == 2
 
1153
C     Stable PS (SPS) point.
 
1154
C     No stability rescue was necessary.
 
1155
C     H == 3
 
1156
C     Unstable PS (UPS) point.
 
1157
C     Stability rescue necessary, and successful.
 
1158
C     H == 4
 
1159
C     Exceptional PS (EPS) point.
 
1160
C     Stability rescue attempted, but unsuccessful.
 
1161
C     
 
1162
C     T == 1
 
1163
C     Default computation (double prec.) was performed.
 
1164
C     T == 2
 
1165
C     Quadruple precision was used for this PS point.
 
1166
C     T == 3
 
1167
C     MadLoop in initialization phase. Only double precision used.
 
1168
C     T == 4
 
1169
C     MadLoop in initialization phase. Quadruple precision used.
 
1170
C     
 
1171
C     U is a number left for future use (always set to 0 for now).
 
1172
C     example: TIR vs OPP usage.
 
1173
C     
 
1174
      IMPLICIT NONE
 
1175
C     
 
1176
C     CONSTANTS
 
1177
C     
 
1178
      INTEGER    NEXTERNAL
 
1179
      PARAMETER (NEXTERNAL=4)
 
1180
C     
 
1181
C     ARGUMENTS 
 
1182
C     
 
1183
      REAL*8 P(0:3,NEXTERNAL)
 
1184
      REAL*8 ANS(3)
 
1185
      REAL*8 PREC_ASKED,PREC_FOUND
 
1186
      INTEGER RET_CODE
 
1187
C     
 
1188
C     GLOBAL VARIABLES
 
1189
C     
 
1190
      REAL*8 USER_STAB_PREC
 
1191
      COMMON/USER_STAB_PREC/USER_STAB_PREC
 
1192
 
 
1193
      INTEGER H,T,U
 
1194
      REAL*8 ACCURACY
 
1195
      COMMON/ACC/ACCURACY,H,T,U
 
1196
 
 
1197
      LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY
 
1198
      COMMON/BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY
 
1199
 
 
1200
C     ----------
 
1201
C     BEGIN CODE
 
1202
C     ----------
 
1203
      USER_STAB_PREC = PREC_ASKED
 
1204
 
 
1205
      CALL SLOOPMATRIX(P,ANS)
 
1206
      IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY.LT.0.0D0)) THEN
 
1207
        BYPASS_CHECK = .TRUE.
 
1208
        CALL SLOOPMATRIX(P,ANS)
 
1209
        BYPASS_CHECK = .FALSE.
 
1210
C       Make sure we correctly return an initialization-type T code
 
1211
        IF (T.EQ.2) T=4
 
1212
        IF (T.EQ.1) T=3
 
1213
      ENDIF
 
1214
 
 
1215
C     Reset it to default value not to affect next runs
 
1216
      USER_STAB_PREC = -1.0D0
 
1217
      PREC_FOUND=ACCURACY
 
1218
      RET_CODE=100*H+10*T+U
 
1219
 
 
1220
      END
 
1221