1
SUBROUTINE ML5_0_IMPROVE_PS_POINT_PRECISION(P)
7
PARAMETER (NEXTERNAL=4)
11
DOUBLE PRECISION P(0:3,NEXTERNAL)
12
REAL*16 QP_P(0:3,NEXTERNAL)
28
CALL ML5_0_MP_IMPROVE_PS_POINT_PRECISION(QP_P)
39
SUBROUTINE ML5_0_MP_IMPROVE_PS_POINT_PRECISION(P)
45
PARAMETER (NEXTERNAL=4)
49
REAL*16 P(0:3,NEXTERNAL)
54
INTEGER ERRCODE,ERRCODETMP
55
REAL*16 NEWP(0:3,NEXTERNAL)
59
LOGICAL ML5_0_MP_IS_PHYSICAL
63
INCLUDE 'MadLoopParams.inc'
71
DATA TOLD_SUPPRESS/.FALSE./
76
C ERROR CODES CONVENTION
78
C 1 :: None physical PS point input
79
C 100-1000 :: Error in the origianl method for restoring
81
C 1000-9999 :: Error when restoring precision ala PSMC
92
C Check the sanity of the original PS point
93
IF (.NOT.ML5_0_MP_IS_PHYSICAL(NEWP,WARNED)) THEN
95
WRITE(*,*) 'ERROR:: The input PS point is not precise enough.'
99
C Now restore the precision
100
IF (IMPROVEPSPOINT.EQ.1) THEN
101
CALL ML5_0_MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE
103
ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN
104
CALL ML5_0_MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE
107
IF (ERRCODE.NE.0) THEN
108
IF (WARNED.LT.20) THEN
109
WRITE(*,*) 'INFO:: Attempting to rescue the precision'
110
$ //' improvement with an alternative method.'
113
IF (IMPROVEPSPOINT.EQ.1) THEN
114
CALL ML5_0_MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP
115
$ ,ERRCODETMP,WARNED)
116
ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN
117
CALL ML5_0_MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP
118
$ ,ERRCODETMP,WARNED)
120
IF (ERRCODETMP.NE.0) GOTO 100
123
C Report to the user or update the PS point.
127
IF (WARNED.LT.20) THEN
128
WRITE(*,*) 'WARNING:: This PS point could not be improved.'
129
$ //' Error code = ',ERRCODE,ERRCODETMP
130
CALL ML5_0_MP_WRITE_MOM(P)
142
IF (WARNED.GE.20.AND..NOT.TOLD_SUPPRESS) THEN
143
WRITE(*,*) 'INFO:: Further warnings from the improve_ps'
144
$ //' routine will now be supressed.'
151
FUNCTION ML5_0_MP_IS_CLOSE(P,NEWP,WARNED)
157
PARAMETER (NEXTERNAL=4)
159
PARAMETER (ZERO=0.0E+00_16)
161
PARAMETER (THRS_CLOSE=1.0E-02_16)
165
REAL*16 P(0:3,NEXTERNAL), NEWP(0:3,NEXTERNAL)
166
LOGICAL ML5_0_MP_IS_CLOSE
173
DOUBLE PRECISION BUFFDP
175
C NOW MAKE SURE THE SHIFTED POINT IS NOT TOO FAR FROM THE ORIGINAL
177
ML5_0_MP_IS_CLOSE = .TRUE.
182
REF2 = REF2 + ABS(P(I,J))
183
REF = REF + ABS(P(I,J)-NEWP(I,J))
187
IF ((REF/REF2).GT.THRS_CLOSE) THEN
188
ML5_0_MP_IS_CLOSE = .FALSE.
189
IF (WARNED.LT.20) THEN
191
WRITE(*,*) 'WARNING:: The improved PS point is too far from'
192
$ //' the original one',BUFFDP
199
FUNCTION ML5_0_MP_IS_PHYSICAL(P,WARNED)
205
PARAMETER (NEXTERNAL=4)
207
PARAMETER (NINITIAL=2)
209
PARAMETER (ZERO=0.0E+00_16)
211
PARAMETER (MP__ZERO=ZERO)
213
PARAMETER (ONE=1.0E+00_16)
215
PARAMETER (TWO=2.0E+00_16)
216
REAL*16 THRES_ONSHELL
217
PARAMETER (THRES_ONSHELL=1.0E-02_16)
218
REAL*16 THRES_FOURMOM
219
PARAMETER (THRES_FOURMOM=1.0E-06_16)
223
REAL*16 P(0:3,NEXTERNAL)
224
LOGICAL ML5_0_MP_IS_PHYSICAL
231
REAL*16 MASSES(NEXTERNAL)
232
DOUBLE PRECISION BUFFDPA,BUFFDPB
236
INCLUDE 'mp_coupl.inc'
247
ML5_0_MP_IS_PHYSICAL = .TRUE.
249
C WE FIRST CHECK THAT THE INPUT PS POINT IS REASONABLY PHYSICAL
250
C FOR THAT WE NEED A REFERENCE SCALE
260
DO J=NINITIAL+1,NEXTERNAL
263
IF ((BUFF/REF).GT.THRES_FOURMOM) THEN
264
IF (WARNED.LT.20) THEN
266
WRITE(*,*) 'ERROR:: Four-momentum conservation is not'
267
$ //' accurate enough, ',BUFFDPA
268
CALL ML5_0_MP_WRITE_MOM(P)
271
ML5_0_MP_IS_PHYSICAL = .FALSE.
274
REF = REF / (ONE*NEXTERNAL)
276
REF=ABS(P(0,I))+ABS(P(1,I))+ABS(P(2,I))+ABS(P(3,I))
277
IF ((SQRT(ABS(P(0,I)**2-P(1,I)**2-P(2,I)**2-P(3,I)**2-MASSES(I)
278
$ **2))/REF).GT.THRES_ONSHELL) THEN
279
IF (WARNED.LT.20) THEN
281
BUFFDPB=(SQRT(ABS(P(0,I)**2-P(1,I)**2-P(2,I)**2-P(3,I)**2
282
$ -MASSES(I)**2))/REF)
283
WRITE(*,*) 'ERROR:: Onshellness of the momentum of'
284
$ //' particle ',I,' of mass ',BUFFDPA,' is not accurate'
285
$ //' enough, ',BUFFDPB
286
CALL ML5_0_MP_WRITE_MOM(P)
289
ML5_0_MP_IS_PHYSICAL = .FALSE.
295
SUBROUTINE ML5_0_WRITE_MOM(P)
298
PARAMETER (NEXTERNAL=4)
300
PARAMETER (NINITIAL=2)
301
DOUBLE PRECISION ZERO
302
PARAMETER (ZERO=0.0D0)
303
DOUBLE PRECISION ML5_0_MDOT
310
DOUBLE PRECISION P(0:3,NEXTERNAL),PSUM(0:3)
314
PSUM(I)=PSUM(I)+P(I,J)
316
DO J=NINITIAL+1,NEXTERNAL
317
PSUM(I)=PSUM(I)-P(I,J)
320
WRITE (*,*) ' Phase space point:'
321
WRITE (*,*) ' ---------------------'
322
WRITE (*,*) ' E | px | py | pz | m '
324
WRITE (*,'(1x,5e27.17)') P(0,I),P(1,I),P(2,I),P(3,I)
325
$ ,SQRT(ABS(ML5_0_MDOT(P(0,I),P(0,I))))
327
WRITE (*,*) ' Four-momentum conservation sum:'
328
WRITE (*,'(1x,4e27.17)') PSUM(0),PSUM(1),PSUM(2),PSUM(3)
329
WRITE (*,*) ' ---------------------'
332
DOUBLE PRECISION FUNCTION ML5_0_MDOT(P1,P2)
334
DOUBLE PRECISION P1(0:3),P2(0:3)
335
ML5_0_MDOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
339
SUBROUTINE ML5_0_MP_WRITE_MOM(P)
342
PARAMETER (NEXTERNAL=4)
344
PARAMETER (NINITIAL=2)
346
PARAMETER (ZERO=0.0E+00_16)
347
REAL*16 ML5_0_MP_MDOT
354
REAL*16 P(0:3,NEXTERNAL),PSUM(0:3),DOT
355
DOUBLE PRECISION DP_P(0:3,NEXTERNAL),DP_PSUM(0:3),DP_DOT
360
PSUM(I)=PSUM(I)+P(I,J)
362
DO J=NINITIAL+1,NEXTERNAL
363
PSUM(I)=PSUM(I)-P(I,J)
367
C The GCC4.7 compiler on SLC machines has trouble to write out
368
C quadruple precision variable with the write(*,*) statement. I
369
C therefore perform the cast by hand
377
WRITE (*,*) ' Phase space point:'
378
WRITE (*,*) ' ---------------------'
379
WRITE (*,*) ' E | px | py | pz | m '
381
DOT=SQRT(ABS(ML5_0_MP_MDOT(P(0,I),P(0,I))))
383
WRITE (*,'(1x,5e27.17)') DP_P(0,I),DP_P(1,I),DP_P(2,I),DP_P(3
386
WRITE (*,*) ' Four-momentum conservation sum:'
387
WRITE (*,'(1x,4e27.17)') DP_PSUM(0),DP_PSUM(1),DP_PSUM(2)
389
WRITE (*,*) ' ---------------------'
392
REAL*16 FUNCTION ML5_0_MP_MDOT(P1,P2)
394
REAL*16 P1(0:3),P2(0:3)
395
ML5_0_MP_MDOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
399
C Rotate_PS rotates the PS point PS (without modifying it)
400
C stores the result in P and for the quadruple precision
401
C version , it also modifies the global variables
402
C PS and MP_DONE accordingly.
404
SUBROUTINE ML5_0_ROTATE_PS(P_IN,P,ROTATION)
410
PARAMETER (NEXTERNAL=4)
414
DOUBLE PRECISION P_IN(0:3,NEXTERNAL),P(0:3,NEXTERNAL)
426
C rotation=1 => (xp=z,yp=-x,zp=-y)
427
IF(ROTATION.EQ.1) THEN
432
C rotation=2 => (xp=-z,yp=y,zp=x)
433
ELSEIF(ROTATION.EQ.2) THEN
449
SUBROUTINE ML5_0_MP_ROTATE_PS(P_IN,P,ROTATION)
455
PARAMETER (NEXTERNAL=4)
459
REAL*16 P_IN(0:3,NEXTERNAL),P(0:3,NEXTERNAL)
469
COMMON/ML5_0_MP_DONE/MP_DONE
476
C rotation=1 => (xp=z,yp=-x,zp=-y)
477
IF(ROTATION.EQ.1) THEN
482
C rotation=2 => (xp=-z,yp=y,zp=x)
483
ELSEIF(ROTATION.EQ.2) THEN
500
C *****************************************************************
501
C Beginning of the routine for restoring precision with V.H. method
502
C *****************************************************************
504
SUBROUTINE ML5_0_MP_ORIG_IMPROVE_PS_POINT_PRECISION(P,ERRCODE
511
PARAMETER (NEXTERNAL=4)
513
PARAMETER (NINITIAL=2)
515
PARAMETER (ZERO=0.0E+00_16)
517
PARAMETER (MP__ZERO=ZERO)
519
PARAMETER (ONE=1.0E+00_16)
521
PARAMETER (TWO=2.0E+00_16)
523
PARAMETER (THRS_TEST=1.0E-15_16)
527
REAL*16 P(0:3,NEXTERNAL)
528
INTEGER ERRCODE, WARNED
532
LOGICAL ML5_0_MP_IS_CLOSE
538
REAL*16 PT(0:3), NEWP(0:3,NEXTERNAL)
539
REAL*16 BUFF,REF,REF2,DISCR
540
REAL*16 MASSES(NEXTERNAL)
541
REAL*16 SHIFTE(2),SHIFTZ(2)
545
INCLUDE 'mp_coupl.inc'
557
C NOW WE MAKE SURE THAT THE PS POINT CAN BE IMPROVED BY THE
564
IF (NINITIAL.NE.2) ERRCODE = 100
566
IF (ABS(P(1,1)/REF).GT.THRS_TEST.OR.ABS(P(2,1)/REF)
567
$ .GT.THRS_TEST.OR.ABS(P(1,2)/REF).GT.THRS_TEST.OR.ABS(P(2,2)/REF)
568
$ .GT.THRS_TEST) ERRCODE = 200
570
IF (MASSES(1).NE.ZERO.OR.MASSES(2).NE.ZERO) ERRCODE = 300
573
IF (P(0,I).LT.ZERO) ERRCODE = 400 + I
576
IF (ERRCODE.NE.0) GOTO 100
578
C WE FIRST SHIFT ALL THE FINAL STATE PARTICLES TO MAKE THEM
584
DO I=NINITIAL+1,NEXTERNAL
587
NEWP(3,I)=SIGN(SQRT(ABS(P(0,I)**2-P(1,I)**2-P(2,I)**2
588
$ -MASSES(I)**2)),P(3,I))
592
PT(J)=PT(J)+NEWP(J,I)
596
C WE CHOOSE P1 IN THE ALGORITHM TO ALWAYS BE THE PARTICLE WITH
598
IF (P(3,1).GT.ZERO) THEN
601
ELSEIF (P(3,2).GT.ZERO) THEN
609
C Now we calculate the shift to bring to P1 and P2
611
C ptotC = {ptotE, ptotX, ptotY, ptotZ};
612
C pm1C = {pm1E + sm1E, pm1X, pm1Y, pm1Z + sm1Z};
613
C {pm0E + sm0E, ptotX - pm1X, ptotY - pm1Y, pm0Z + sm0Z};
614
C sol = Solve[{ptotC[[1]] - pm1C[[1]] - pm0C[[1]] == 0,
615
C ptotC[[4]] - pm1C[[4]] - pm0C[[4]] == 0,
616
C pm1C[[1]]^2 - pm1C[[2]]^2 - pm1C[[3]]^2 - pm1C[[4]]^2 == m1M^2,
617
C pm0C[[1]]^2 - pm0C[[2]]^2 - pm0C[[3]]^2 - pm0C[[4]]^2 == m2M^2},
618
C {sm1E, sm1Z, sm0E, sm0Z}] // FullSimplify;
619
C (solC[[1]] /. {m1M -> 0, m2M -> 0} /. {pm1X -> 0, pm1Y -> 0})
622
DISCR = -PT(0)**2 + PT(1)**2 + PT(2)**2 + PT(3)**2
623
IF (DISCR.LT.ZERO) DISCR = -DISCR
625
SHIFTE(1) = (PT(0)*(-TWO*P(0,P1)*PT(0) + PT(0)**2 + PT(1)**2 +
626
$ PT(2)**2) + (TWO*P(0,P1) - PT(0))*PT(3)**2 + PT(3)*DISCR)/(TWO
627
$ *(PT(0) - PT(3))*(PT(0) + PT(3)))
628
SHIFTE(2) = -(PT(0)*(TWO*P(0,P2)*PT(0) - PT(0)**2 + PT(1)**2 +
629
$ PT(2)**2) + (-TWO*P(0,P2) + PT(0))*PT(3)**2 + PT(3)*DISCR)
630
$ /(TWO*(PT(0) - PT(3))*(PT(0) + PT(3)))
631
SHIFTZ(1) = (-TWO*P(3,P1)*(PT(0)**2 - PT(3)**2) + PT(3)*(PT(0)*
632
$ *2 + PT(1)**2 + PT(2)**2 - PT(3)**2) + PT(0)*DISCR)/(TWO*(PT(0)
634
SHIFTZ(2) = -(TWO*P(3,P2)*(PT(0)**2 - PT(3)**2) + PT(3)*(-PT(0)*
635
$ *2 + PT(1)**2 + PT(2)**2 + PT(3)**2) + PT(0)*DISCR)/(TWO*(PT(0)
637
NEWP(0,P1) = P(0,P1)+SHIFTE(1)
638
NEWP(3,P1) = P(3,P1)+SHIFTZ(1)
639
NEWP(0,P2) = P(0,P2)+SHIFTE(2)
640
NEWP(3,P2) = P(3,P2)+SHIFTZ(2)
645
DO I=NINITIAL+1,NEXTERNAL
652
IF (.NOT.ML5_0_MP_IS_CLOSE(P,NEWP,WARNED)) THEN
667
C *****************************************************************
668
C Beginning of the routine for restoring precision a la PSMC
669
C *****************************************************************
671
SUBROUTINE ML5_0_MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE
678
PARAMETER (NEXTERNAL=4)
680
PARAMETER (NINITIAL=2)
682
PARAMETER (ZERO=0.0E+00_16)
684
PARAMETER (MP__ZERO=ZERO)
686
PARAMETER (ONE=1.0E+00_16)
688
PARAMETER (TWO=2.0E+00_16)
689
REAL*16 CONSISTENCY_THRES
690
PARAMETER (CONSISTENCY_THRES=1.0E-25_16)
693
PARAMETER (NAPPROXZEROS=3)
698
REAL*16 P(0:3,NEXTERNAL)
699
INTEGER ERRCODE,ERROR,WARNED
703
LOGICAL ML5_0_MP_IS_CLOSE
708
REAL*16 NEWP(0:3,NEXTERNAL), PBUFF(0:3)
709
REAL*16 BUFF, BUFF2, XSCALE, APPROX_ZEROS(NAPPROXZEROS)
710
REAL*16 MASSES(NEXTERNAL)
714
INCLUDE 'mp_coupl.inc'
728
C Define the seeds which should be tried
729
APPROX_ZEROS(1)=1.0E+00_16
730
APPROX_ZEROS(2)=1.1E+00_16
731
APPROX_ZEROS(3)=0.9E+00_16
733
C Start by copying the momenta
740
C First make sur that the space like momentum is exactly conserved
746
PBUFF(J)=PBUFF(J)+NEWP(J,I)
749
DO I=NINITIAL+1,NEXTERNAL-1
751
PBUFF(J)=PBUFF(J)-NEWP(J,I)
755
NEWP(J,NEXTERNAL)=PBUFF(J)
758
C Now find the 'x' rescaling factor
760
CALL ML5_0_FINDX(NEWP,APPROX_ZEROS(I),XSCALE,ERROR)
764
ERRCODE=ERRCODE+(10**(I-1))*ERROR
767
IF (WARNED.LT.20) THEN
768
WRITE(*,*) 'WARNING:: Could not find the proper rescaling'
769
$ //' factor x. Restoring precision ala PSMC will therefore not'
773
IF (ERRCODE.LT.1000) THEN
780
C Apply the rescaling
783
NEWP(J,I)=NEWP(J,I)*XSCALE
787
C Now restore exact onshellness of the particles.
791
BUFF=BUFF+NEWP(J,I)**2
801
BUFF2=BUFF2+NEWP(0,I)
803
DO I=NINITIAL+1,NEXTERNAL
805
BUFF2=BUFF2+NEWP(0,I)
807
IF ((ABS(BUFF)/BUFF2).GT.CONSISTENCY_THRES) THEN
808
IF (WARNED.LT.20) THEN
809
WRITE(*,*) 'WARNING:: The consistency check in the a la PSMC'
810
$ //' precision restoring algorithm failed. The result will'
811
$ //' therefore not be used.'
818
IF (.NOT.ML5_0_MP_IS_CLOSE(P,NEWP,WARNED)) THEN
834
SUBROUTINE ML5_0_FINDX(P,SEED,XSCALE,ERROR)
840
PARAMETER (NEXTERNAL=4)
842
PARAMETER (NINITIAL=2)
844
PARAMETER (ZERO=0.0E+00_16)
846
PARAMETER (MP__ZERO=ZERO)
848
PARAMETER (ONE=1.0E+00_16)
850
PARAMETER (TWO=2.0E+00_16)
851
INTEGER MAXITERATIONS
852
PARAMETER (MAXITERATIONS=8)
854
PARAMETER (CONVERGED=1.0E-26_16)
858
REAL*16 P(0:3,NEXTERNAL),SEED,XSCALE
864
REAL*16 PVECSQ(NEXTERNAL)
865
REAL*16 XN, XNP1,FVAL,DVAL
877
PVECSQ(I)=P(1,I)**2+P(2,I)**2+P(3,I)**2
881
CALL ML5_0_FUNCT(PVECSQ(1),XN,.FALSE.,ERR, FVAL)
886
CALL ML5_0_FUNCT(PVECSQ(1),XN,.TRUE.,ERR, DVAL)
892
IF((ABS(((XNP1-XN)*TWO)/(XNP1+XN))).LT.CONVERGED) THEN
902
C For good measure, we iterate one last time
903
CALL ML5_0_FUNCT(PVECSQ(1),XN,.FALSE.,ERR, FVAL)
908
CALL ML5_0_FUNCT(PVECSQ(1),XN,.TRUE.,ERR, DVAL)
914
XSCALE=XN-(FVAL/DVAL)
920
SUBROUTINE ML5_0_FUNCT(PVECSQ,X,DERIVATIVE,ERROR,RES)
926
PARAMETER (NEXTERNAL=4)
928
PARAMETER (NINITIAL=2)
930
PARAMETER (ZERO=0.0E+00_16)
932
PARAMETER (MP__ZERO=ZERO)
934
PARAMETER (ONE=1.0E+00_16)
936
PARAMETER (TWO=2.0E+00_16)
940
REAL*16 PVECSQ(NEXTERNAL),X,RES
948
REAL*16 MASSES(NEXTERNAL)
952
INCLUDE 'mp_coupl.inc'
968
IF (I.LE.NINITIAL) THEN
973
BUFF=MASSES(I)**2+PVECSQ(I)*X**2
974
IF (BUFF.LT.ZERO) THEN
980
RES=RES + FACTOR*((X*PVECSQ(I))/SQRT(BUFF))
982
RES=RES + FACTOR*SQRT(BUFF)