80
86
PCT(I,0)=PCT(I,0)+PL(I,J)
87
ABSPCT(I)=ABSPCT(I)+ABS(PL(I,J))
83
IF (CHECKPCONSERVATION) THEN
84
IF (PCT(0,0).GT.1.D-6) THEN
85
WRITE(*,*) 'energy is not conserved ',PCT(0,0)
86
STOP 'energy is not conserved'
87
ELSEIF (PCT(1,0).GT.1.D-6) THEN
88
WRITE(*,*) 'px is not conserved ',PCT(1,0)
89
STOP 'px is not conserved'
90
ELSEIF (PCT(2,0).GT.1.D-6) THEN
91
WRITE(*,*) 'py is not conserved ',PCT(2,0)
92
STOP 'py is not conserved'
93
ELSEIF (PCT(3,0).GT.1.D-6) THEN
94
WRITE(*,*) 'pz is not conserved ',PCT(3,0)
95
STOP 'pz is not conserved'
90
REF_P = MAX(ABSPCT(0), ABSPCT(1),ABSPCT(2),ABSPCT(3))
92
ABSPCT(I) = MAX(REF_P*1E-6, ABSPCT(I))
94
IF (CHECKPCONSERVATION.AND.REF_P.GT.1D-8) THEN
95
IF ((PCT(0,0)/ABSPCT(0)).GT.1.D-6) THEN
96
WRITE(*,*) 'energy is not conserved (flag:CT95)',PCT(0,0)
97
STOP 'energy is not conserved (flag:CT95)'
98
ELSEIF ((PCT(1,0)/ABSPCT(1)).GT.1.D-6) THEN
99
WRITE(*,*) 'px is not conserved (flag:CT95)',PCT(1,0)
100
STOP 'px is not conserved (flag:CT95)'
101
ELSEIF ((PCT(2,0)/ABSPCT(2)).GT.1.D-6) THEN
102
WRITE(*,*) 'py is not conserved (flag:CT95)',PCT(2,0)
103
STOP 'py is not conserved (flag:CT95)'
104
ELSEIF ((PCT(3,0)/ABSPCT(3)).GT.1.D-6) THEN
105
WRITE(*,*) 'pz is not conserved (flag:CT95)',PCT(3,0)
106
STOP 'pz is not conserved (flag:CT95)'
154
SUBROUTINE ML5_0_LOOP_5(W1, W2, W3, W4, W5, M1, M2, M3, M4, M5
155
$ , RANK, SQUAREDSOINDEX, LOOPNUM)
165
SUBROUTINE ML5_0_BUILD_KINEMATIC_MATRIX(NLOOPLINE,P_LOOP,M2L
168
C Helper function that compute the loop kinematic matrix with
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.
179
REAL*8 P_LOOP(NLOOPLINE,0:3)
180
COMPLEX*16 M2L(NLOOPLINE)
181
COMPLEX*16 S_MAT(NLOOPLINE,NLOOPLINE)
185
INCLUDE 'MadLoopParams.inc'
191
REAL*8 REF_NORMALIZATION
201
S_MAT(I,J)=-(M2L(I)+M2L(J))
203
DIFFSQ = (DCMPLX(P_LOOP(I,0),0.0D0)-DCMPLX(P_LOOP(J,0)
206
DIFFSQ = DIFFSQ - (DCMPLX(P_LOOP(I,K),0.0D0)
207
$ -DCMPLX(P_LOOP(J,K),0.0D0))**2
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
213
IF(ABS(M2L(I)).NE.0.0D0)THEN
214
IF(ABS((DIFFSQ-M2L(I))/M2L(I)).LT.OSTHRES)THEN
218
IF(ABS(M2L(J)).NE.0.0D0)THEN
219
IF(ABS((DIFFSQ-M2L(J))/M2L(J)).LT.OSTHRES)THEN
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
229
REF_NORMALIZATION = REF_NORMALIZATION + ABS(P_LOOP(I,K))
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))
245
SUBROUTINE ML5_0_MP_BUILD_KINEMATIC_MATRIX(NLOOPLINE,P_LOOP,M2L
248
C Helper function that compute the loop kinematic matrix with
250
C NLOOPLINE : Number of loop lines
251
C P_LOOP : List of external momenta running in the loop, i.e.
252
C q_i in the denominator (l_i+q_i)**2-m_i**2
253
C M2L : List of complex-valued masses running in the loop.
254
C S_MAT(N,N): Kinematic matrix output.
259
REAL*16 P_LOOP(NLOOPLINE,0:3)
260
COMPLEX*32 M2L(NLOOPLINE)
261
COMPLEX*32 S_MAT(NLOOPLINE,NLOOPLINE)
265
INCLUDE 'MadLoopParams.inc'
271
REAL*16 REF_NORMALIZATION
281
S_MAT(I,J)=-(M2L(I)+M2L(J))
283
DIFFSQ = (CMPLX(P_LOOP(I,0),0.0E0_16,KIND=16)
284
$ -CMPLX(P_LOOP(J,0),0.0E0_16,KIND=16))**2
286
DIFFSQ = DIFFSQ - (CMPLX(P_LOOP(I,K),0.0E0_16,KIND=16)
287
$ -CMPLX(P_LOOP(J,K),0.0E0_16,KIND=16))**2
289
C Default value of the kinematic matrix
290
S_MAT(I,J)=DIFFSQ-M2L(I)-M2L(J)
291
C And we now test various thresholds. Normaly, at most one
293
IF(ABS(M2L(I)).NE.0.0E0_16)THEN
294
IF(ABS((DIFFSQ-M2L(I))/M2L(I)).LT.OSTHRES)THEN
298
IF(ABS(M2L(J)).NE.0.0E0_16)THEN
299
IF(ABS((DIFFSQ-M2L(J))/M2L(J)).LT.OSTHRES)THEN
303
C Choose what seems the most appropriate way to compare
304
C massless onshellness.
305
REF_NORMALIZATION=0.0E0_16
306
C Here, we chose to base the threshold only on the energy
309
REF_NORMALIZATION = REF_NORMALIZATION + ABS(P_LOOP(I,K))
312
REF_NORMALIZATION = (REF_NORMALIZATION/2.0E0_16)**2
313
IF(REF_NORMALIZATION.NE.0.0E0_16)THEN
314
IF(ABS(DIFFSQ/REF_NORMALIZATION).LT.OSTHRES)THEN
315
S_MAT(I,J)=-(M2L(I)+M2L(J))
327
SUBROUTINE ML5_0_LOOP_5(W1, W2, W3, W4, W5, M1, M2, M3, M4, M5,
328
$ RANK, SQUAREDSOINDEX, LOOPNUM)
156
329
INTEGER NEXTERNAL
157
330
PARAMETER (NEXTERNAL=5)
158
331
INTEGER NLOOPLINE
258
C Determine it uses qp or not
260
IF(CTMODE.GE.4)DOING_QP=.TRUE.
261
444
C Choose the correct loop library
262
445
CALL ML5_0_CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
446
$ ,ID,DOING_QP,I_LIB)
264
447
IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
266
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOIN
267
$ DEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
449
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1
450
$ ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
452
C Tensor Integral Reduction is used
270
453
CALL ML5_0_TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL
271
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
454
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)
455
$ ,S(SQUAREDSOINDEX,LOOPNUM))
275
458
LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
385
C Determine it uses qp or not
387
IF(CTMODE.GE.4)DOING_QP=.TRUE.
388
581
C Choose the correct loop library
389
582
CALL ML5_0_CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
583
$ ,ID,DOING_QP,I_LIB)
391
584
IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
393
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOIN
394
$ DEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
586
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1
587
$ ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
589
C Tensor Integral Reduction is used
397
590
CALL ML5_0_TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL
398
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
591
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)
592
$ ,S(SQUAREDSOINDEX,LOOPNUM))
402
595
LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
511
C Determine it uses qp or not
513
IF(CTMODE.GE.4)DOING_QP=.TRUE.
514
717
C Choose the correct loop library
515
718
CALL ML5_0_CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
719
$ ,ID,DOING_QP,I_LIB)
517
720
IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
519
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOIN
520
$ DEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
722
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1
723
$ ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
725
C Tensor Integral Reduction is used
523
726
CALL ML5_0_TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL
524
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
727
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)
728
$ ,S(SQUAREDSOINDEX,LOOPNUM))
528
731
LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
632
C Determine it uses qp or not
634
IF(CTMODE.GE.4)DOING_QP=.TRUE.
635
848
C Choose the correct loop library
636
849
CALL ML5_0_CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
850
$ ,ID,DOING_QP,I_LIB)
638
851
IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
640
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOIN
641
$ DEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
853
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1
854
$ ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
856
C Tensor Integral Reduction is used
644
857
CALL ML5_0_TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL
645
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
858
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)
859
$ ,S(SQUAREDSOINDEX,LOOPNUM))
649
862
LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
754
C Determine it uses qp or not
756
IF(CTMODE.GE.4)DOING_QP=.TRUE.
757
980
C Choose the correct loop library
758
981
CALL ML5_0_CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
982
$ ,ID,DOING_QP,I_LIB)
760
983
IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
762
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOIN
763
$ DEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
985
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1
986
$ ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
988
C Tensor Integral Reduction is used
766
989
CALL ML5_0_TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL
767
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
990
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)
991
$ ,S(SQUAREDSOINDEX,LOOPNUM))
771
994
LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
878
C Determine it uses qp or not
880
IF(CTMODE.GE.4)DOING_QP=.TRUE.
881
1114
C Choose the correct loop library
882
1115
CALL ML5_0_CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
1116
$ ,ID,DOING_QP,I_LIB)
884
1117
IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
886
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOIN
887
$ DEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
1119
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1
1120
$ ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
1122
C Tensor Integral Reduction is used
890
1123
CALL ML5_0_TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL
891
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
1124
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)
1125
$ ,S(SQUAREDSOINDEX,LOOPNUM))
895
1128
LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
1002
C Determine it uses qp or not
1004
IF(CTMODE.GE.4)DOING_QP=.TRUE.
1005
1248
C Choose the correct loop library
1006
1249
CALL ML5_0_CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
1250
$ ,ID,DOING_QP,I_LIB)
1008
1251
IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
1010
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOIN
1011
$ DEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
1253
CALL ML5_0_CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1
1254
$ ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
1256
C Tensor Integral Reduction is used
1014
1257
CALL ML5_0_TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL
1015
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
1258
$ ,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)
1259
$ ,S(SQUAREDSOINDEX,LOOPNUM))
1019
1262
LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)