← Back to branch summary

~maddevelopers/mg5amcnlo/2.7.1.3

« back to all changes in this revision

Viewing changes to tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%CT_interface.f

  • Committer: olivier Mattelaer
  • Date: 2016-05-12 11:00:18 UTC
  • mfrom: (262.1.150 2.3.4)
  • Revision ID: olivier.mattelaer@uclouvain.be-20160512110018-sevb79f0wm4g8mpp
pass to 2.4.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%CT_interface.f
 
1
C     ===========================================
 
2
C     ===== Beginning of CutTools Interface =====
 
3
C     ===========================================
 
4
      SUBROUTINE CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
 
5
C     
 
6
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
 
7
C     By the MadGraph5_aMC@NLO Development Team
 
8
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
 
9
C     
 
10
C     Interface between MG5 and CutTools.
 
11
C     
 
12
C     Process: u d~ > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
 
13
C     Process: c s~ > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
 
14
C     
 
15
C     
 
16
C     CONSTANTS 
 
17
C     
 
18
      INTEGER    NEXTERNAL
 
19
      PARAMETER (NEXTERNAL=3)
 
20
      LOGICAL CHECKPCONSERVATION
 
21
      PARAMETER (CHECKPCONSERVATION=.TRUE.)
 
22
      REAL*8 NORMALIZATION
 
23
      PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
 
24
     $ *2))
 
25
C     
 
26
C     ARGUMENTS 
 
27
C     
 
28
      INTEGER NLOOPLINE, RANK
 
29
      REAL*8 PL(0:3,NLOOPLINE)
 
30
      REAL*8 PCT(0:3,0:NLOOPLINE-1),ABSPCT(0:3)
 
31
      REAL*8 REF_P
 
32
      COMPLEX*16 M2L(NLOOPLINE)
 
33
      COMPLEX*16 M2LCT(0:NLOOPLINE-1)
 
34
      COMPLEX*16 RES(3)
 
35
      LOGICAL STABLE
 
36
C     
 
37
C     LOCAL VARIABLES 
 
38
C     
 
39
      COMPLEX*16 R1, ACC
 
40
      INTEGER I, J, K
 
41
      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
 
42
      COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT
 
43
     $ ,NINJAINIT
 
44
C     
 
45
C     EXTERNAL FUNCTIONS
 
46
C     
 
47
      EXTERNAL LOOPNUM
 
48
      EXTERNAL MPLOOPNUM
 
49
C     
 
50
C     GLOBAL VARIABLES
 
51
C     
 
52
      INCLUDE 'coupl.inc'
 
53
      INTEGER CTMODE
 
54
      REAL*8 LSCALE
 
55
      COMMON/CT/LSCALE,CTMODE
 
56
 
 
57
      INTEGER ID,SQSOINDEX,R
 
58
      COMMON/LOOP/ID,SQSOINDEX,R
 
59
 
 
60
C     ----------
 
61
C     BEGIN CODE
 
62
C     ----------
 
63
 
 
64
C     INITIALIZE CUTTOOLS IF NEEDED
 
65
      IF (CTINIT) THEN
 
66
        CTINIT=.FALSE.
 
67
        CALL INITCT()
 
68
      ENDIF
 
69
 
 
70
C     YOU CAN FIND THE DETAILS ABOUT THE DIFFERENT CTMODE AT THE
 
71
C      BEGINNING OF THE FILE CTS_CUTS.F90 IN THE CUTTOOLS DISTRIBUTION
 
72
 
 
73
C     CONVERT THE MASSES TO BE COMPLEX
 
74
      DO I=1,NLOOPLINE
 
75
        M2LCT(I-1)=M2L(I)
 
76
      ENDDO
 
77
 
 
78
C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
 
79
      DO I=0,3
 
80
        ABSPCT(I)=0.D0
 
81
        DO J=0,(NLOOPLINE-1)
 
82
          PCT(I,J)=0.D0
 
83
        ENDDO
 
84
      ENDDO
 
85
      DO I=0,3
 
86
        DO J=1,NLOOPLINE
 
87
          PCT(I,0)=PCT(I,0)+PL(I,J)
 
88
          ABSPCT(I)=ABSPCT(I)+ABS(PL(I,J))
 
89
        ENDDO
 
90
      ENDDO
 
91
      REF_P = MAX(ABSPCT(0), ABSPCT(1),ABSPCT(2),ABSPCT(3))
 
92
      DO I=0,3
 
93
        ABSPCT(I) = MAX(REF_P*1E-6, ABSPCT(I))
 
94
      ENDDO
 
95
      IF (CHECKPCONSERVATION.AND.REF_P.GT.1D-8) THEN
 
96
        IF ((PCT(0,0)/ABSPCT(0)).GT.1.D-6) THEN
 
97
          WRITE(*,*) 'energy is not conserved (flag:CT95)',PCT(0,0)
 
98
          STOP 'energy is not conserved (flag:CT95)'
 
99
        ELSEIF ((PCT(1,0)/ABSPCT(1)).GT.1.D-6) THEN
 
100
          WRITE(*,*) 'px is not conserved (flag:CT95)',PCT(1,0)
 
101
          STOP 'px is not conserved (flag:CT95)'
 
102
        ELSEIF ((PCT(2,0)/ABSPCT(2)).GT.1.D-6) THEN
 
103
          WRITE(*,*) 'py is not conserved (flag:CT95)',PCT(2,0)
 
104
          STOP 'py is not conserved (flag:CT95)'
 
105
        ELSEIF ((PCT(3,0)/ABSPCT(3)).GT.1.D-6) THEN
 
106
          WRITE(*,*) 'pz is not conserved (flag:CT95)',PCT(3,0)
 
107
          STOP 'pz is not conserved (flag:CT95)'
 
108
        ENDIF
 
109
      ENDIF
 
110
      DO I=0,3
 
111
        DO J=1,(NLOOPLINE-1)
 
112
          DO K=1,J
 
113
            PCT(I,J)=PCT(I,J)+PL(I,K)
 
114
          ENDDO
 
115
        ENDDO
 
116
      ENDDO
 
117
 
 
118
      CALL CTSXCUT(CTMODE,LSCALE,MU_R,NLOOPLINE,LOOPNUM,MPLOOPNUM,RANK
 
119
     $ ,PCT,M2LCT,RES,ACC,R1,STABLE)
 
120
      RES(1)=NORMALIZATION*2.0D0*DBLE(RES(1))
 
121
      RES(2)=NORMALIZATION*2.0D0*DBLE(RES(2))
 
122
      RES(3)=NORMALIZATION*2.0D0*DBLE(RES(3))
 
123
C     WRITE(*,*) 'CutTools: Loop ID',ID,' =',RES(1),RES(2),RES(3)
 
124
      END
 
125
 
 
126
      SUBROUTINE INITCT()
 
127
C     
 
128
C     INITIALISATION OF CUTTOOLS
 
129
C     
 
130
C     LOCAL VARIABLES 
 
131
C     
 
132
      REAL*8 THRS
 
133
      LOGICAL EXT_NUM_FOR_R1
 
134
C     
 
135
C     GLOBAL VARIABLES 
 
136
C     
 
137
      INCLUDE 'MadLoopParams.inc'
 
138
C     ----------
 
139
C     BEGIN CODE
 
140
C     ----------
 
141
 
 
142
C     DEFAULT PARAMETERS FOR CUTTOOLS
 
143
C     -------------------------------  
 
144
C     THRS1 IS THE PRECISION LIMIT BELOW WHICH THE MP ROUTINES
 
145
C      ACTIVATES
 
146
      THRS=CTSTABTHRES
 
147
C     LOOPLIB SET WHAT LIBRARY CT USES
 
148
C     1 -> LOOPTOOLS
 
149
C     2 -> AVH
 
150
C     3 -> QCDLOOP
 
151
      LOOPLIB=CTLOOPLIBRARY
 
152
C     MADLOOP'S NUMERATOR IN THE OPEN LOOP IS MUCH FASTER THAN THE
 
153
C      RECONSTRUCTED ONE IN CT. SO WE BETTER USE MADLOOP ONE IN THIS
 
154
C      CASE.
 
155
      EXT_NUM_FOR_R1=.TRUE.
 
156
C     -------------------------------     
 
157
 
 
158
C     The initialization below is for CT v1.8.+
 
159
      CALL CTSINIT(THRS,LOOPLIB,EXT_NUM_FOR_R1)
 
160
C     The initialization below is for the older stable CT v1.7, still
 
161
C      used for now in the beta release.
 
162
C     CALL CTSINIT(THRS,LOOPLIB)
 
163
 
 
164
      END
 
165
 
 
166
      SUBROUTINE BUILD_KINEMATIC_MATRIX(NLOOPLINE,P_LOOP,M2L,S_MAT)
 
167
C     
 
168
C     Helper function that compute the loop kinematic matrix with
 
169
C      proper thresholds
 
170
C     NLOOPLINE : Number of loop lines
 
171
C     P_LOOP    : List of external momenta running in the loop, i.e.
 
172
C      q_i in the denominator (l_i+q_i)**2-m_i**2
 
173
C     M2L       : List of complex-valued masses running in the loop.
 
174
C     S_MAT(N,N): Kinematic matrix output.
 
175
C     
 
176
C     ARGUMENTS
 
177
C     
 
178
      INTEGER NLOOPLINE
 
179
      REAL*8 P_LOOP(NLOOPLINE,0:3)
 
180
      COMPLEX*16 M2L(NLOOPLINE)
 
181
      COMPLEX*16 S_MAT(NLOOPLINE,NLOOPLINE)
 
182
C     
 
183
C     GLOBAL VARIABLES
 
184
C     
 
185
      INCLUDE 'MadLoopParams.inc'
 
186
C     
 
187
C     LOCAL VARIABLES
 
188
C     
 
189
      INTEGER I,J,K
 
190
      COMPLEX*16 DIFFSQ
 
191
      REAL*8 REF_NORMALIZATION
 
192
 
 
193
C     ----------
 
194
C     BEGIN CODE
 
195
C     ----------
 
196
 
 
197
      DO I=1,NLOOPLINE
 
198
        DO J=1,NLOOPLINE
 
199
 
 
200
          IF(I.EQ.J)THEN
 
201
            S_MAT(I,J)=-(M2L(I)+M2L(J))
 
202
          ELSE
 
203
            DIFFSQ = (DCMPLX(P_LOOP(I,0),0.0D0)-DCMPLX(P_LOOP(J,0)
 
204
     $       ,0.0D0))**2
 
205
            DO K=1,3
 
206
              DIFFSQ = DIFFSQ - (DCMPLX(P_LOOP(I,K),0.0D0)
 
207
     $         -DCMPLX(P_LOOP(J,K),0.0D0))**2
 
208
            ENDDO
 
209
C           Default value of the kinematic matrix
 
210
            S_MAT(I,J)=DIFFSQ-M2L(I)-M2L(J)
 
211
C           And we now test various thresholds. Normaly, at most one
 
212
C            applies.
 
213
            IF(ABS(M2L(I)).NE.0.0D0)THEN
 
214
              IF(ABS((DIFFSQ-M2L(I))/M2L(I)).LT.OSTHRES)THEN
 
215
                S_MAT(I,J)=-M2L(J)
 
216
              ENDIF
 
217
            ENDIF
 
218
            IF(ABS(M2L(J)).NE.0.0D0)THEN
 
219
              IF(ABS((DIFFSQ-M2L(J))/M2L(J)).LT.OSTHRES)THEN
 
220
                S_MAT(I,J)=-M2L(I)
 
221
              ENDIF
 
222
            ENDIF
 
223
C           Choose what seems the most appropriate way to compare
 
224
C           massless onshellness.
 
225
            REF_NORMALIZATION=0.0D0
 
226
C           Here, we chose to base the threshold only on the energy
 
227
C            component
 
228
            DO K=0,0
 
229
              REF_NORMALIZATION = REF_NORMALIZATION + ABS(P_LOOP(I,K))
 
230
     $          + ABS(P_LOOP(J,K))
 
231
            ENDDO
 
232
            REF_NORMALIZATION = (REF_NORMALIZATION/2.0D0)**2
 
233
            IF(REF_NORMALIZATION.NE.0.0D0)THEN
 
234
              IF(ABS(DIFFSQ/REF_NORMALIZATION).LT.OSTHRES)THEN
 
235
                S_MAT(I,J)=-(M2L(I)+M2L(J))
 
236
              ENDIF
 
237
            ENDIF
 
238
          ENDIF
 
239
 
 
240
        ENDDO
 
241
      ENDDO
 
242
 
 
243
      END
 
244
 
 
245
      SUBROUTINE MP_BUILD_KINEMATIC_MATRIX(NLOOPLINE,P_LOOP,M2L,S_MAT)
 
246
C     
 
247
C     Helper function that compute the loop kinematic matrix with
 
248
C      proper thresholds
 
249
C     NLOOPLINE : Number of loop lines
 
250
C     P_LOOP    : List of external momenta running in the loop, i.e.
 
251
C      q_i in the denominator (l_i+q_i)**2-m_i**2
 
252
C     M2L       : List of complex-valued masses running in the loop.
 
253
C     S_MAT(N,N): Kinematic matrix output.
 
254
C     
 
255
C     ARGUMENTS
 
256
C     
 
257
      INTEGER NLOOPLINE
 
258
      REAL*16 P_LOOP(NLOOPLINE,0:3)
 
259
      COMPLEX*32 M2L(NLOOPLINE)
 
260
      COMPLEX*32 S_MAT(NLOOPLINE,NLOOPLINE)
 
261
C     
 
262
C     GLOBAL VARIABLES
 
263
C     
 
264
      INCLUDE 'MadLoopParams.inc'
 
265
C     
 
266
C     LOCAL VARIABLES
 
267
C     
 
268
      INTEGER I,J,K
 
269
      COMPLEX*32 DIFFSQ
 
270
      REAL*16 REF_NORMALIZATION
 
271
 
 
272
C     ----------
 
273
C     BEGIN CODE
 
274
C     ----------
 
275
 
 
276
      DO I=1,NLOOPLINE
 
277
        DO J=1,NLOOPLINE
 
278
 
 
279
          IF(I.EQ.J)THEN
 
280
            S_MAT(I,J)=-(M2L(I)+M2L(J))
 
281
          ELSE
 
282
            DIFFSQ = (CMPLX(P_LOOP(I,0),0.0E0_16,KIND=16)
 
283
     $       -CMPLX(P_LOOP(J,0),0.0E0_16,KIND=16))**2
 
284
            DO K=1,3
 
285
              DIFFSQ = DIFFSQ - (CMPLX(P_LOOP(I,K),0.0E0_16,KIND=16)
 
286
     $         -CMPLX(P_LOOP(J,K),0.0E0_16,KIND=16))**2
 
287
            ENDDO
 
288
C           Default value of the kinematic matrix
 
289
            S_MAT(I,J)=DIFFSQ-M2L(I)-M2L(J)
 
290
C           And we now test various thresholds. Normaly, at most one
 
291
C            applies.
 
292
            IF(ABS(M2L(I)).NE.0.0E0_16)THEN
 
293
              IF(ABS((DIFFSQ-M2L(I))/M2L(I)).LT.OSTHRES)THEN
 
294
                S_MAT(I,J)=-M2L(J)
 
295
              ENDIF
 
296
            ENDIF
 
297
            IF(ABS(M2L(J)).NE.0.0E0_16)THEN
 
298
              IF(ABS((DIFFSQ-M2L(J))/M2L(J)).LT.OSTHRES)THEN
 
299
                S_MAT(I,J)=-M2L(I)
 
300
              ENDIF
 
301
            ENDIF
 
302
C           Choose what seems the most appropriate way to compare
 
303
C           massless onshellness.
 
304
            REF_NORMALIZATION=0.0E0_16
 
305
C           Here, we chose to base the threshold only on the energy
 
306
C            component
 
307
            DO K=0,0
 
308
              REF_NORMALIZATION = REF_NORMALIZATION + ABS(P_LOOP(I,K))
 
309
     $          + ABS(P_LOOP(J,K))
 
310
            ENDDO
 
311
            REF_NORMALIZATION = (REF_NORMALIZATION/2.0E0_16)**2
 
312
            IF(REF_NORMALIZATION.NE.0.0E0_16)THEN
 
313
              IF(ABS(DIFFSQ/REF_NORMALIZATION).LT.OSTHRES)THEN
 
314
                S_MAT(I,J)=-(M2L(I)+M2L(J))
 
315
              ENDIF
 
316
            ENDIF
 
317
          ENDIF
 
318
 
 
319
        ENDDO
 
320
      ENDDO
 
321
 
 
322
      END
 
323
 
 
324
 
 
325
 
 
326
C     ===========================================
 
327
C     ===== Beginning of Ninja interface  =====
 
328
C     ===========================================
 
329
 
 
330
      SUBROUTINE NINJA_LOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
 
331
C     
 
332
C     Module used
 
333
C     
 
334
      USE MNINJA
 
335
C     
 
336
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
 
337
C     By the MadGraph5_aMC@NLO Development Team
 
338
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
 
339
C     
 
340
C     Interface between MG5 and Ninja.
 
341
C     
 
342
C     Process: u d~ > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
 
343
C     Process: c s~ > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
 
344
C     
 
345
C     
 
346
C     CONSTANTS 
 
347
C     
 
348
      INTEGER    NEXTERNAL
 
349
      PARAMETER (NEXTERNAL=3)
 
350
      LOGICAL CHECKPCONSERVATION
 
351
      PARAMETER (CHECKPCONSERVATION=.TRUE.)
 
352
      REAL*8 NORMALIZATION
 
353
      PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
 
354
     $ *2))
 
355
      INTEGER NLOOPGROUPS
 
356
      PARAMETER (NLOOPGROUPS=1)
 
357
C     These are constants related to the split orders
 
358
      INTEGER NSQUAREDSO
 
359
      PARAMETER (NSQUAREDSO=1)
 
360
      INCLUDE 'loop_max_coefs.inc'
 
361
C     
 
362
C     ARGUMENTS 
 
363
C     
 
364
      INTEGER NLOOPLINE, RANK
 
365
      REAL*8 PL(0:3,NLOOPLINE)
 
366
      COMPLEX*16 M2L(NLOOPLINE)
 
367
      COMPLEX*16 RES(3)
 
368
      LOGICAL STABLE
 
369
C     
 
370
C     LOCAL VARIABLES 
 
371
C     
 
372
      REAL*8 P_TMP(0:3,0:NLOOPLINE-1),ABSP_TMP(0:3)
 
373
      REAL*8 REF_P
 
374
      REAL*8 P_NINJA(0:3,NLOOPLINE)
 
375
      REAL*8 P_S_MAT(NLOOPLINE,0:3)
 
376
      COMPLEX*16 M2L_NINJA(NLOOPLINE)
 
377
      COMPLEX*16 NINJA_RES(0:2)
 
378
      COMPLEX*16 R1
 
379
      INTEGER NINJA_STATUS
 
380
      INTEGER I, J, K
 
381
      REAL*8 PDEN_DUMMY(0:3,NLOOPLINE-1)
 
382
 
 
383
      COMPLEX*16 S_MAT(NLOOPLINE,NLOOPLINE)
 
384
      REAL*8 REAL_S_MAT(NLOOPLINE,NLOOPLINE)
 
385
 
 
386
      INTEGER CURR_MAXCOEF
 
387
      COMPLEX*16, ALLOCATABLE :: TENSORCOEFS(:)
 
388
 
 
389
C     
 
390
C     GLOBAL VARIABLES
 
391
C     
 
392
      INCLUDE 'coupl.inc'
 
393
 
 
394
      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
 
395
      COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT
 
396
     $ ,NINJAINIT
 
397
 
 
398
      REAL*8 LSCALE
 
399
      INTEGER CTMODE
 
400
      COMMON/CT/LSCALE,CTMODE
 
401
 
 
402
      INTEGER ID,SQSOINDEX,R
 
403
      COMMON/LOOP/ID,SQSOINDEX,R
 
404
      COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
 
405
      COMMON/LCOEFS/LOOPCOEFS
 
406
 
 
407
      LOGICAL FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
 
408
      COMMON/FPE_IN_REDUCTION/FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
 
409
 
 
410
C     ----------
 
411
C     BEGIN CODE
 
412
C     ----------
 
413
 
 
414
C     For the direction test, we must switch the direction in which
 
415
C      the loop is read for CTMode equal to 2 or 4.
 
416
      CALL SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN_DUMMY,M2L)
 
417
 
 
418
C     The CT initialization is also performed here if not done already
 
419
C      because it calls MPINIT of OneLOop which is necessary on some
 
420
C      system
 
421
      IF (CTINIT) THEN
 
422
        CTINIT=.FALSE.
 
423
        CALL INITCT()
 
424
      ENDIF
 
425
 
 
426
C     INITIALIZE NINJA IF NEEDED
 
427
      IF (NINJAINIT) THEN
 
428
        NINJAINIT=.FALSE.
 
429
        CALL INITNINJA()
 
430
      ENDIF
 
431
 
 
432
C     CONVERT THE MASSES TO BE COMPLEX
 
433
      DO I=1,NLOOPLINE
 
434
        M2L_NINJA(I)=M2L(I)
 
435
      ENDDO
 
436
 
 
437
C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO NINJA
 
438
C      CONVENTIONS
 
439
      DO I=0,3
 
440
        ABSP_TMP = 0.D0
 
441
        DO J=0,(NLOOPLINE-1)
 
442
          P_TMP(I,J)=0.D0
 
443
        ENDDO
 
444
      ENDDO
 
445
      DO I=0,3
 
446
        DO J=1,NLOOPLINE
 
447
          P_TMP(I,0)=P_TMP(I,0)+PL(I,J)
 
448
          ABSP_TMP(I) = ABSP_TMP(I)+ABS(PL(I,J))
 
449
        ENDDO
 
450
      ENDDO
 
451
      REF_P = MAX(ABSP_TMP(0), ABSP_TMP(1),ABSP_TMP(2),ABSP_TMP(3))
 
452
      DO I=0,3
 
453
        ABSP_TMP(I) = MAX(REF_P*1E-6, ABSP_TMP(I))
 
454
      ENDDO
 
455
 
 
456
      IF (CHECKPCONSERVATION.AND.REF_P.GT.1D-8) THEN
 
457
        IF ((P_TMP(0,0)/ABSP_TMP(0)).GT.1.D-6) THEN
 
458
          WRITE(*,*) 'energy is not conserved (flag:CT692)',P_TMP(0,0)
 
459
          STOP 'energy is not conserved (flag:CT692)'
 
460
        ELSEIF ((P_TMP(1,0)/ABSP_TMP(1)).GT.1.D-6) THEN
 
461
          WRITE(*,*) 'px is not conserved (flag:CT692)',P_TMP(1,0)
 
462
          STOP 'px is not conserved (flag:CT692)'
 
463
        ELSEIF ((P_TMP(2,0)/ABSP_TMP(2)).GT.1.D-6) THEN
 
464
          WRITE(*,*) 'py is not conserved (flag:CT692)',P_TMP(2,0)
 
465
          STOP 'py is not conserved (flag:CT692)'
 
466
        ELSEIF ((P_TMP(3,0)/ABSP_TMP(3)).GT.1.D-6) THEN
 
467
          WRITE(*,*) 'pz is not conserved (flag:CT692)',P_TMP(3,0)
 
468
          STOP 'pz is not conserved (flag:CT692)'
 
469
        ENDIF
 
470
      ENDIF
 
471
      DO I=0,3
 
472
        DO J=1,(NLOOPLINE-1)
 
473
          DO K=1,J
 
474
            P_TMP(I,J)=P_TMP(I,J)+PL(I,K)
 
475
          ENDDO
 
476
        ENDDO
 
477
      ENDDO
 
478
C     In Ninja, the loop line index starts at 1
 
479
      DO I=0,NLOOPLINE-1
 
480
        P_NINJA(0,I+1) = P_TMP(0,I)
 
481
        P_NINJA(1,I+1) = P_TMP(1,I)
 
482
        P_NINJA(2,I+1) = P_TMP(2,I)
 
483
        P_NINJA(3,I+1) = P_TMP(3,I)
 
484
      ENDDO
 
485
 
 
486
C     Number of coefficients for the current rank
 
487
      CURR_MAXCOEF = 0
 
488
      DO I=0,RANK
 
489
        CURR_MAXCOEF=CURR_MAXCOEF+(3+I)*(2+I)*(1+I)/6
 
490
      ENDDO
 
491
C     Now write the tensor coefficients for Ninja
 
492
C     It should never be allocated at this stage
 
493
      IF (.NOT. ALLOCATED(TENSORCOEFS)) THEN
 
494
        ALLOCATE(TENSORCOEFS(0:CURR_MAXCOEF-1))
 
495
      ENDIF
 
496
      DO I=0,CURR_MAXCOEF-1
 
497
        TENSORCOEFS(I) = LOOPCOEFS(I,SQSOINDEX,ID)
 
498
      ENDDO
 
499
C     The loop momentum is in fact q_loop -> -q_loop, so that the
 
500
C     coefficients must be changed accordingly
 
501
      CALL INVERT_MOMENTA_IN_POLYNOMIAL(CURR_MAXCOEF,TENSORCOEFS)
 
502
 
 
503
C     Compute the kinematic matrix
 
504
      DO J=1,NLOOPLINE
 
505
        DO I=0,3
 
506
          P_S_MAT(J,I)=P_NINJA(I,J)
 
507
        ENDDO
 
508
      ENDDO
 
509
      CALL BUILD_KINEMATIC_MATRIX(NLOOPLINE,P_S_MAT,M2L,S_MAT)
 
510
 
 
511
      DO I=1,NLOOPLINE
 
512
        DO J=1,NLOOPLINE
 
513
          REAL_S_MAT(I,J) = DBLE(S_MAT(I,J)+M2L(I)+M2L(J))
 
514
        ENDDO
 
515
      ENDDO
 
516
 
 
517
C     Below is the call specifying the kinematic matrix
 
518
      CALL NINJA_TENSOR_EVALUATE(TENSORCOEFS,NLOOPLINE,RANK,REAL_S_MAT
 
519
     $ ,P_NINJA,M2L,MU_R**2,NINJA_RES,R1,NINJA_STATUS)
 
520
 
 
521
C     Below is the call without specification of the kinematic matrix
 
522
C     call ninja_tensor_evaluate(TENSORCOEFS,NLOOPLINE,RANK,P_NINJA,M2L
 
523
C     ,MU_R**2,NINJA_RES,R1,NINJA_STATUS)
 
524
 
 
525
C     If a floating point exception was found in Ninja (e.g. exactly
 
526
C      zero gram. det.)
 
527
C     Then warn loop_matrix.f so that it will flag this kinematic
 
528
C      point as unstable no matter what.
 
529
      IF (NINJA_STATUS.EQ.NINJA_UNSTABLE_KINEMATICS) THEN
 
530
        FPE_IN_DP_REDUCTION = .TRUE.
 
531
      ENDIF
 
532
 
 
533
C     Make sure to deallocate the tensor of coefficients
 
534
      IF (ALLOCATED(TENSORCOEFS)) THEN
 
535
        DEALLOCATE(TENSORCOEFS)
 
536
      ENDIF
 
537
 
 
538
      RES(1)=NORMALIZATION*2.0D0*DBLE(NINJA_RES(0))
 
539
      RES(2)=NORMALIZATION*2.0D0*DBLE(NINJA_RES(1))
 
540
      RES(3)=NORMALIZATION*2.0D0*DBLE(NINJA_RES(2))
 
541
C     WRITE(*,*) 'Ninja:  Loop ID',ID,' =',RES(1),RES(2),RES(3)
 
542
      END
 
543
 
 
544
C     
 
545
C     Quadruple precision version of loop_ninja
 
546
C     
 
547
      SUBROUTINE MP_NINJA_LOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
 
548
C     
 
549
C     Module used
 
550
C     
 
551
      USE MNINJA
 
552
C     
 
553
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
 
554
C     By the MadGraph5_aMC@NLO Development Team
 
555
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
 
556
C     
 
557
C     Interface between MG5 and Ninja.
 
558
C     
 
559
C     Process: u d~ > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
 
560
C     Process: c s~ > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
 
561
C     
 
562
C     
 
563
C     CONSTANTS 
 
564
C     
 
565
      INTEGER    NEXTERNAL
 
566
      PARAMETER (NEXTERNAL=3)
 
567
      LOGICAL CHECKPCONSERVATION
 
568
      PARAMETER (CHECKPCONSERVATION=.TRUE.)
 
569
      REAL*8 NORMALIZATION
 
570
      PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
 
571
     $ *2))
 
572
      INTEGER NLOOPGROUPS
 
573
      PARAMETER (NLOOPGROUPS=1)
 
574
C     These are constants related to the split orders
 
575
      INTEGER NSQUAREDSO
 
576
      PARAMETER (NSQUAREDSO=1)
 
577
      INCLUDE 'loop_max_coefs.inc'
 
578
C     
 
579
C     ARGUMENTS 
 
580
C     
 
581
      INTEGER NLOOPLINE, RANK
 
582
      REAL*16 PL(0:3,NLOOPLINE)
 
583
      COMPLEX*16 M2L(NLOOPLINE)
 
584
      COMPLEX*16 RES(3)
 
585
      LOGICAL STABLE
 
586
C     
 
587
C     LOCAL VARIABLES 
 
588
C     
 
589
      REAL*16 NINJA_SCALE
 
590
      REAL*16 P_TMP(0:3,0:NLOOPLINE-1), ABSP_TMP(0:3)
 
591
      REAL*8 REF_P
 
592
      REAL(KI_QNIN) MP_P_NINJA(0:3,NLOOPLINE)
 
593
      REAL*16 MP_P(0:3,NLOOPLINE)
 
594
      REAL*16 P_S_MAT(NLOOPLINE,0:3)
 
595
      COMPLEX*32 MP_M2L(NLOOPLINE)
 
596
      COMPLEX(KI_QNIN) MP_M2L_NINJA(NLOOPLINE)
 
597
      COMPLEX(KI_QNIN) NINJA_RES(0:2)
 
598
      COMPLEX(KI_QNIN) NINJA_R1
 
599
      COMPLEX*16 R1
 
600
      COMPLEX*16 DP_RES(0:2)
 
601
      INTEGER NINJA_STATUS
 
602
      INTEGER I, J, K
 
603
      REAL*16 PDEN_DUMMY(0:3,NLOOPLINE-1)
 
604
 
 
605
      COMPLEX*32 MP_S_MAT(NLOOPLINE,NLOOPLINE)
 
606
      REAL*16 MP_REAL_S_MAT(NLOOPLINE,NLOOPLINE)
 
607
      REAL(KI_QNIN) MP_REAL_S_MAT_NINJA(NLOOPLINE,NLOOPLINE)
 
608
 
 
609
      INTEGER CURR_MAXCOEF
 
610
      COMPLEX*32, ALLOCATABLE :: MP_TENSORCOEFS(:)
 
611
      COMPLEX(KI_QNIN), ALLOCATABLE :: MP_NINJA_TENSORCOEFS(:)
 
612
 
 
613
C     
 
614
C     GLOBAL VARIABLES
 
615
C     
 
616
      INCLUDE 'coupl.inc'
 
617
 
 
618
      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
 
619
      COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT
 
620
     $ ,NINJAINIT
 
621
 
 
622
      REAL*8 LSCALE
 
623
      INTEGER CTMODE
 
624
      COMMON/CT/LSCALE,CTMODE
 
625
 
 
626
      INTEGER ID,SQSOINDEX,R
 
627
      COMMON/LOOP/ID,SQSOINDEX,R
 
628
      COMPLEX*32 MP_LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
 
629
      COMMON/MP_LCOEFS/MP_LOOPCOEFS
 
630
 
 
631
      LOGICAL FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
 
632
      COMMON/FPE_IN_REDUCTION/FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
 
633
 
 
634
C     ----------
 
635
C     BEGIN CODE
 
636
C     ----------
 
637
 
 
638
C     Cast the masses in complex quadruple precision
 
639
      DO I=1,NLOOPLINE
 
640
        MP_M2L(I) = CMPLX(M2L(I),KIND=16)
 
641
      ENDDO
 
642
 
 
643
C     For the direction test, we must switch the direction in which
 
644
C      the loop is read for CTMode equal to 2 or 4.
 
645
      CALL MP_SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN_DUMMY,MP_M2L)
 
646
 
 
647
C     The CT initialization is also performed here if not done already
 
648
C      because it calls MPINIT of OneLOop which is necessary on some
 
649
C      system
 
650
      IF (CTINIT) THEN
 
651
        CTINIT=.FALSE.
 
652
        CALL INITCT()
 
653
      ENDIF
 
654
 
 
655
C     INITIALIZE NINJA IF NEEDED
 
656
      IF (NINJAINIT) THEN
 
657
        NINJAINIT=.FALSE.
 
658
        CALL INITNINJA()
 
659
      ENDIF
 
660
 
 
661
C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO NINJA
 
662
C      CONVENTIONS
 
663
      DO I=0,3
 
664
        ABSP_TMP(I)=0.E0+0_16
 
665
        DO J=0,(NLOOPLINE-1)
 
666
          P_TMP(I,J)=0.E0+0_16
 
667
        ENDDO
 
668
      ENDDO
 
669
      DO I=0,3
 
670
        DO J=1,NLOOPLINE
 
671
          P_TMP(I,0)=P_TMP(I,0)+PL(I,J)
 
672
          ABSP_TMP(I)=ABSP_TMP(I)+ABS(PL(I,J))
 
673
        ENDDO
 
674
      ENDDO
 
675
      REF_P = MAX(ABSP_TMP(0), ABSP_TMP(1),ABSP_TMP(2),ABSP_TMP(3))
 
676
      DO I=0,3
 
677
        ABSP_TMP(I) = MAX(REF_P*1E-6, ABSP_TMP(I))
 
678
      ENDDO
 
679
      IF (CHECKPCONSERVATION.AND.REF_P.GT.1.E-8_16) THEN
 
680
        IF ((P_TMP(0,0)/ABSP_TMP(0)).GT.1.E-6_16) THEN
 
681
          WRITE(*,*) 'energy is not conserved (flag:CT968)'
 
682
     $     ,DBLE(P_TMP(0,0))
 
683
          STOP 'energy is not conserved (flag:CT968)'
 
684
        ELSEIF ((P_TMP(1,0)/ABSP_TMP(1)).GT.1.E-6_16) THEN
 
685
          WRITE(*,*) 'px is not conserved (flag:CT968)',DBLE(P_TMP(1,0)
 
686
     $     )
 
687
          STOP 'px is not conserved (flag:CT968)'
 
688
        ELSEIF ((P_TMP(2,0)/ABSP_TMP(2)).GT.1.E-6_16) THEN
 
689
          WRITE(*,*) 'py is not conserved (flag:CT968)',DBLE(P_TMP(2,0)
 
690
     $     )
 
691
          STOP 'py is not conserved (flag:CT968)'
 
692
        ELSEIF ((P_TMP(3,0)/ABSP_TMP(3)).GT.1.E-6_16) THEN
 
693
          WRITE(*,*) 'pz is not conserved (flag:CT968)',DBLE(P_TMP(3,0)
 
694
     $     )
 
695
          STOP 'pz is not conserved (flag:CT968)'
 
696
        ENDIF
 
697
      ENDIF
 
698
      DO I=0,3
 
699
        DO J=1,(NLOOPLINE-1)
 
700
          DO K=1,J
 
701
            P_TMP(I,J)=P_TMP(I,J)+PL(I,K)
 
702
          ENDDO
 
703
        ENDDO
 
704
      ENDDO
 
705
C     In Ninja, the loop line index starts at 1
 
706
      DO I=0,NLOOPLINE-1
 
707
        MP_P(0,I+1) = P_TMP(0,I)
 
708
        MP_P(1,I+1) = P_TMP(1,I)
 
709
        MP_P(2,I+1) = P_TMP(2,I)
 
710
        MP_P(3,I+1) = P_TMP(3,I)
 
711
      ENDDO
 
712
 
 
713
C     Number of coefficients for the current rank
 
714
      CURR_MAXCOEF = 0
 
715
      DO I=0,RANK
 
716
        CURR_MAXCOEF=CURR_MAXCOEF+(3+I)*(2+I)*(1+I)/6
 
717
      ENDDO
 
718
C     Now write the tensor coefficients for Ninja
 
719
C     It should never be allocated at this stage
 
720
      IF (.NOT. ALLOCATED(MP_TENSORCOEFS)) THEN
 
721
        ALLOCATE(MP_TENSORCOEFS(0:CURR_MAXCOEF-1))
 
722
      ENDIF
 
723
      IF (.NOT. ALLOCATED(MP_NINJA_TENSORCOEFS)) THEN
 
724
        ALLOCATE(MP_NINJA_TENSORCOEFS(0:CURR_MAXCOEF-1))
 
725
      ENDIF
 
726
      DO I=0,CURR_MAXCOEF-1
 
727
        MP_TENSORCOEFS(I) = MP_LOOPCOEFS(I,SQSOINDEX,ID)
 
728
      ENDDO
 
729
C     The loop momentum is in fact q_loop -> -q_loop, so that the
 
730
C     coefficients must be changed accordingly
 
731
      CALL MP_INVERT_MOMENTA_IN_POLYNOMIAL(CURR_MAXCOEF,MP_TENSORCOEFS)
 
732
 
 
733
C     Compute the kinematic matrix
 
734
      DO J=1,NLOOPLINE
 
735
        DO I=0,3
 
736
          P_S_MAT(J,I)=MP_P(I,J)
 
737
        ENDDO
 
738
      ENDDO
 
739
      CALL MP_BUILD_KINEMATIC_MATRIX(NLOOPLINE,P_S_MAT,MP_M2L,MP_S_MAT)
 
740
 
 
741
      DO I=1,NLOOPLINE
 
742
        DO J=1,NLOOPLINE
 
743
          MP_REAL_S_MAT(I,J) = REAL(MP_S_MAT(I,J)+MP_M2L(I)+MP_M2L(J)
 
744
     $     ,KIND=16)
 
745
        ENDDO
 
746
      ENDDO
 
747
 
 
748
C     Now typecast to Ninja's quadruple precision format
 
749
      DO I=0,CURR_MAXCOEF-1
 
750
        MP_NINJA_TENSORCOEFS(I)=CMPLX(MP_TENSORCOEFS(I),KIND=KI_QNIN)
 
751
      ENDDO
 
752
      DO I=1,NLOOPLINE
 
753
        DO J=1,NLOOPLINE
 
754
          MP_REAL_S_MAT_NINJA(I,J) = REAL(MP_REAL_S_MAT(I,J)
 
755
     $     ,KIND=KI_QNIN)
 
756
        ENDDO
 
757
      ENDDO
 
758
      DO I=1,NLOOPLINE
 
759
        MP_M2L_NINJA(I)=CMPLX(MP_M2L(I),KIND=KI_QNIN)
 
760
      ENDDO
 
761
      DO I=1,NLOOPLINE
 
762
        MP_P_NINJA(0,I) = REAL(MP_P(0,I),KIND=KI_QNIN)
 
763
        MP_P_NINJA(1,I) = REAL(MP_P(1,I),KIND=KI_QNIN)
 
764
        MP_P_NINJA(2,I) = REAL(MP_P(2,I),KIND=KI_QNIN)
 
765
        MP_P_NINJA(3,I) = REAL(MP_P(3,I),KIND=KI_QNIN)
 
766
      ENDDO
 
767
      NINJA_SCALE = REAL(MU_R**2,KIND=KI_QNIN)
 
768
 
 
769
 
 
770
C     Below is the call specifying the kinematic matrix
 
771
      CALL NINJA_TENSOR_EVALUATE(MP_NINJA_TENSORCOEFS,NLOOPLINE,RANK
 
772
     $ ,MP_REAL_S_MAT_NINJA,MP_P_NINJA,MP_M2L_NINJA,NINJA_SCALE
 
773
     $ ,NINJA_RES,NINJA_R1,NINJA_STATUS)
 
774
C     Below is the call without specification of the kinematic matrix
 
775
C     call ninja_tensor_evaluate(MP_NINJA_TENSORCOEFS,NLOOPLINE,RANK,MP
 
776
C     _P_NINJA,MP_M2L_NINJA,NINJA_SCALE,NINJA_RES,NINJA_R1,NINJA_STATUS)
 
777
C     
 
778
 
 
779
C     If a floating point exception was found in Ninja (e.g. exactly
 
780
C      zero gram. det.)
 
781
C     Then warn loop_matrix.f so that it will flag this kinematic
 
782
C      point as unstable no matter what.
 
783
      IF (NINJA_STATUS.EQ.NINJA_UNSTABLE_KINEMATICS) THEN
 
784
        FPE_IN_QP_REDUCTION = .TRUE.
 
785
      ENDIF
 
786
 
 
787
C     Typecast the result back
 
788
      R1 = DCMPLX(R1)
 
789
      DO I=0,2
 
790
        DP_RES(I)=DCMPLX(NINJA_RES(I))
 
791
      ENDDO
 
792
 
 
793
C     Make sure to deallocate the tensor of coefficients
 
794
      IF (ALLOCATED(MP_TENSORCOEFS)) THEN
 
795
        DEALLOCATE(MP_TENSORCOEFS)
 
796
      ENDIF
 
797
      IF (ALLOCATED(MP_NINJA_TENSORCOEFS)) THEN
 
798
        DEALLOCATE(MP_NINJA_TENSORCOEFS)
 
799
      ENDIF
 
800
 
 
801
      RES(1)=NORMALIZATION*2.0D0*DBLE(DP_RES(0))
 
802
      RES(2)=NORMALIZATION*2.0D0*DBLE(DP_RES(1))
 
803
      RES(3)=NORMALIZATION*2.0D0*DBLE(DP_RES(2))
 
804
C     WRITE(*,*) 'QP Ninja:  Loop ID',ID,' =',RES(1),RES(2),RES(3)
 
805
      END
 
806
 
 
807
      SUBROUTINE INITNINJA()
 
808
C     
 
809
C     Module used
 
810
C     
 
811
      USE MNINJA
 
812
C     
 
813
C     Initialization of Ninja 
 
814
C     
 
815
C     LOCAL VARIABLES 
 
816
C     
 
817
      INTEGER LOOPLIB
 
818
C     
 
819
C     GLOBAL VARIABLES 
 
820
C     
 
821
      INCLUDE 'MadLoopParams.inc'
 
822
C     ----------
 
823
C     BEGIN CODE
 
824
C     ----------
 
825
 
 
826
C     LOOPLIB SET WHAT LIBRARY NINJA USES
 
827
C     1 -> LOOPTOOLS
 
828
C     2 -> AVH
 
829
C     3 -> QCDLOOP
 
830
      IF (CTLOOPLIBRARY.EQ.1) THEN
 
831
        WRITE(*,*) 'Warning in Ninja initialization. LoopTools is not'
 
832
     $   //' supported by the Ninja interface. It will use OneLOop'
 
833
     $   //' instead.'
 
834
        LOOPLIB = 1
 
835
      ELSEIF (CTLOOPLIBRARY.EQ.3) THEN
 
836
        WRITE(*,*) 'Warning in Ninja initialization. LoopTools is not'
 
837
     $   //' supported by the Ninja interface. It will use OneLOop'
 
838
     $   //' instead.'
 
839
        LOOPLIB = 1
 
840
      ELSEIF (CTLOOPLIBRARY.EQ.2) THEN
 
841
        LOOPLIB = 1
 
842
      ELSE
 
843
        WRITE(*,*) 'Error in Ninja initialization. Loop library ID='
 
844
     $   ,CTLOOPLIBRARY,' is not supported. Change variable'
 
845
     $   //' CTLoopLibrary in MadLoopParams.dat.'
 
846
        STOP 1
 
847
      ENDIF
 
848
      CALL NINJA_SET_INTEGRAL_LIBRARY(LOOPLIB)
 
849
 
 
850
      END
 
851
 
 
852
      SUBROUTINE LOOP_3(W1, W2, W3, M1, M2, M3,  RANK, SQUAREDSOINDEX,
 
853
     $  LOOPNUM)
 
854
      INTEGER    NEXTERNAL
 
855
      PARAMETER (NEXTERNAL=3)
 
856
      INTEGER    NLOOPLINE
 
857
      PARAMETER (NLOOPLINE=3)
 
858
      INTEGER    NWAVEFUNCS
 
859
      PARAMETER (NWAVEFUNCS=3)
 
860
      INTEGER    NLOOPGROUPS
 
861
      PARAMETER (NLOOPGROUPS=1)
 
862
      INTEGER    NCOMB
 
863
      PARAMETER (NCOMB=12)
 
864
C     These are constants related to the split orders
 
865
      INTEGER    NSQUAREDSO
 
866
      PARAMETER (NSQUAREDSO=1)
 
867
C     
 
868
C     ARGUMENTS 
 
869
C     
 
870
      INTEGER W1, W2, W3
 
871
      COMPLEX*16 M1, M2, M3
 
872
 
 
873
      INTEGER RANK, LSYMFACT
 
874
      INTEGER LOOPNUM, SQUAREDSOINDEX
 
875
C     
 
876
C     LOCAL VARIABLES 
 
877
C     
 
878
      REAL*8 PL(0:3,NLOOPLINE)
 
879
      REAL*16 MP_PL(0:3,NLOOPLINE)
 
880
      COMPLEX*16 M2L(NLOOPLINE)
 
881
      INTEGER PAIRING(NLOOPLINE),WE(3)
 
882
      INTEGER I, J, K, TEMP,I_LIB
 
883
      LOGICAL COMPLEX_MASS,DOING_QP
 
884
C     
 
885
C     GLOBAL VARIABLES
 
886
C     
 
887
      INCLUDE 'MadLoopParams.inc'
 
888
      INTEGER ID,SQSOINDEX,R
 
889
      COMMON/LOOP/ID,SQSOINDEX,R
 
890
 
 
891
      LOGICAL CHECKPHASE, HELDOUBLECHECKED
 
892
      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED
 
893
 
 
894
      INTEGER HELOFFSET
 
895
      INTEGER GOODHEL(NCOMB)
 
896
      LOGICAL GOODAMP(NSQUAREDSO,NLOOPGROUPS)
 
897
      COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET
 
898
 
 
899
      COMPLEX*16 LOOPRES(3,NSQUAREDSO,NLOOPGROUPS)
 
900
      LOGICAL S(NSQUAREDSO,NLOOPGROUPS)
 
901
      COMMON/LOOPRES/LOOPRES,S
 
902
 
 
903
 
 
904
      COMPLEX*16 W(20,NWAVEFUNCS)
 
905
      COMMON/W/W
 
906
      COMPLEX*32 MP_W(20,NWAVEFUNCS)
 
907
      COMMON/MP_W/MP_W
 
908
 
 
909
      REAL*8 LSCALE
 
910
      INTEGER CTMODE
 
911
      COMMON/CT/LSCALE,CTMODE
 
912
      INTEGER LIBINDEX
 
913
      COMMON/I_LIB/LIBINDEX
 
914
 
 
915
C     ----------
 
916
C     BEGIN CODE
 
917
C     ----------
 
918
 
 
919
C     Determine it uses qp or not
 
920
      DOING_QP = (CTMODE.GE.4)
 
921
 
 
922
      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(SQUAREDSOIND
 
923
     $EX,LOOPNUM)) THEN
 
924
        WE(1)=W1
 
925
        WE(2)=W2
 
926
        WE(3)=W3
 
927
        M2L(1)=M3**2
 
928
        M2L(2)=M1**2
 
929
        M2L(3)=M2**2
 
930
        DO I=1,NLOOPLINE
 
931
          PAIRING(I)=1
 
932
        ENDDO
 
933
 
 
934
        R=RANK
 
935
        ID=LOOPNUM
 
936
        SQSOINDEX=SQUAREDSOINDEX
 
937
        DO I=0,3
 
938
          TEMP=1
 
939
          DO J=1,NLOOPLINE
 
940
            PL(I,J)=0.D0
 
941
            IF (DOING_QP) THEN
 
942
              MP_PL(I,J)=0.0E+0_16
 
943
            ENDIF
 
944
            DO K=TEMP,(TEMP+PAIRING(J)-1)
 
945
              PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
 
946
              IF (DOING_QP) THEN
 
947
                MP_PL(I,J)=MP_PL(I,J)-REAL(MP_W(1+I,WE(K)),KIND=16)
 
948
              ENDIF
 
949
            ENDDO
 
950
            TEMP=TEMP+PAIRING(J)
 
951
          ENDDO
 
952
        ENDDO
 
953
C       Determine whether the integral is with complex masses or not
 
954
C       since some reduction libraries, e.g.PJFry++ and IREGI, are
 
955
C        still
 
956
C       not able to deal with complex masses
 
957
        COMPLEX_MASS=.FALSE.
 
958
        DO I=1,NLOOPLINE
 
959
          IF(DIMAG(M2L(I)).EQ.0D0)CYCLE
 
960
          IF(ABS(DIMAG(M2L(I)))/MAX(ABS(M2L(I)),1D-2).GT.1D-15)THEN
 
961
            COMPLEX_MASS=.TRUE.
 
962
            EXIT
 
963
          ENDIF
 
964
        ENDDO
 
965
C       Choose the correct loop library
 
966
        CALL CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS,ID
 
967
     $   ,DOING_QP,I_LIB)
 
968
        IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
 
969
C         CutTools is used
 
970
          CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX
 
971
     $     ,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
 
972
        ELSEIF (MLREDUCTIONLIB(I_LIB).EQ.6) THEN
 
973
C         Ninja is used
 
974
          IF (.NOT.DOING_QP) THEN
 
975
            CALL NINJA_LOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1
 
976
     $       ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
 
977
          ELSE
 
978
            CALL MP_NINJA_LOOP(NLOOPLINE,MP_PL,M2L,RANK,LOOPRES(1
 
979
     $       ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
 
980
          ENDIF
 
981
        ELSE
 
982
C         Tensor Integral Reduction is used 
 
983
          CALL TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL,M2L
 
984
     $     ,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
 
985
     $     ,LOOPNUM))
 
986
        ENDIF
 
987
      ELSE
 
988
        LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
 
989
        LOOPRES(2,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
 
990
        LOOPRES(3,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
 
991
        S(SQUAREDSOINDEX,LOOPNUM)=.TRUE.
 
992
      ENDIF
 
993
      END
 
994