← Back to branch summary

~maddevelopers/mg5amcnlo/2.5.3

« back to all changes in this revision

Viewing changes to tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/improve_ps.f

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/improve_ps.f
1
 
      SUBROUTINE ML5_0_IMPROVE_PS_POINT_PRECISION(P)
2
 
      IMPLICIT NONE
3
 
C     
4
 
C     CONSTANTS
5
 
C     
6
 
      INTEGER    NEXTERNAL
7
 
      PARAMETER (NEXTERNAL=5)
8
 
C     
9
 
C     ARGUMENTS 
10
 
C     
11
 
      DOUBLE PRECISION P(0:3,NEXTERNAL)
12
 
      REAL*16 QP_P(0:3,NEXTERNAL)
13
 
C     
14
 
C     LOCAL VARIABLES 
15
 
C     
16
 
      INTEGER I,J
17
 
 
18
 
C     ----------
19
 
C     BEGIN CODE
20
 
C     ----------
21
 
 
22
 
      DO I=1,NEXTERNAL
23
 
        DO J=0,3
24
 
          QP_P(J,I)=P(J,I)
25
 
        ENDDO
26
 
      ENDDO
27
 
 
28
 
      CALL ML5_0_MP_IMPROVE_PS_POINT_PRECISION(QP_P)
29
 
 
30
 
      DO I=1,NEXTERNAL
31
 
        DO J=0,3
32
 
          P(J,I)=QP_P(J,I)
33
 
        ENDDO
34
 
      ENDDO
35
 
 
36
 
      END
37
 
 
38
 
 
39
 
      SUBROUTINE ML5_0_MP_IMPROVE_PS_POINT_PRECISION(P)
40
 
      IMPLICIT NONE
41
 
C     
42
 
C     CONSTANTS
43
 
C     
44
 
      INTEGER    NEXTERNAL
45
 
      PARAMETER (NEXTERNAL=5)
46
 
C     
47
 
C     ARGUMENTS 
48
 
C     
49
 
      REAL*16 P(0:3,NEXTERNAL)
50
 
C     
51
 
C     LOCAL VARIABLES 
52
 
C     
53
 
      INTEGER I,J
54
 
      INTEGER ERRCODE,ERRCODETMP
55
 
      REAL*16 NEWP(0:3,NEXTERNAL)
56
 
C     
57
 
C     FUNCTIONS
58
 
C     
59
 
      LOGICAL ML5_0_MP_IS_PHYSICAL
60
 
C     
61
 
C     SAVED VARIABLES
62
 
C     
63
 
      INCLUDE 'MadLoopParams.inc'
64
 
C     
65
 
C     SAVED VARIABLES
66
 
C     
67
 
      INTEGER WARNED
68
 
      DATA WARNED/0/
69
 
 
70
 
      LOGICAL TOLD_SUPPRESS
71
 
      DATA TOLD_SUPPRESS/.FALSE./
72
 
C     ----------
73
 
C     BEGIN CODE
74
 
C     ----------
75
 
 
76
 
C     ERROR CODES CONVENTION
77
 
C     
78
 
C     1         ::  None physical PS point input
79
 
C     100-1000  ::  Error in the origianl method for restoring
80
 
C      precision
81
 
C     1000-9999 ::  Error when restoring precision ala PSMC
82
 
C     
83
 
      ERRCODETMP=0
84
 
      ERRCODE=0
85
 
 
86
 
      DO J=1,NEXTERNAL
87
 
        DO I=0,3
88
 
          NEWP(I,J)=P(I,J)
89
 
        ENDDO
90
 
      ENDDO
91
 
 
92
 
C     Check the sanity of the original PS point
93
 
      IF (.NOT.ML5_0_MP_IS_PHYSICAL(NEWP,WARNED)) THEN
94
 
        ERRCODE = 1
95
 
        WRITE(*,*) 'ERROR:: The input PS point is not precise enough.'
96
 
        GOTO 100
97
 
      ENDIF
98
 
 
99
 
C     Now restore the precision
100
 
      IF (IMPROVEPSPOINT.EQ.1) THEN
101
 
        CALL ML5_0_MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE
102
 
     $   ,WARNED)
103
 
      ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN
104
 
        CALL ML5_0_MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE
105
 
     $   ,WARNED)
106
 
      ENDIF
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.'
111
 
          WARNED=WARNED+1
112
 
        ENDIF
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)
119
 
        ENDIF
120
 
        IF (ERRCODETMP.NE.0) GOTO 100
121
 
      ENDIF
122
 
 
123
 
C     Report to the user or update the PS point.
124
 
 
125
 
      GOTO 101
126
 
 100  CONTINUE
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)
131
 
        WARNED = WARNED +1
132
 
      ENDIF
133
 
      GOTO 102
134
 
 101  CONTINUE
135
 
      DO J=1,NEXTERNAL
136
 
        DO I=0,3
137
 
          P(I,J)=NEWP(I,J)
138
 
        ENDDO
139
 
      ENDDO
140
 
 102  CONTINUE
141
 
 
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.'
145
 
        TOLD_SUPPRESS=.TRUE.
146
 
      ENDIF
147
 
 
148
 
      END
149
 
 
150
 
 
151
 
      FUNCTION ML5_0_MP_IS_CLOSE(P,NEWP,WARNED)
152
 
      IMPLICIT NONE
153
 
C     
154
 
C     CONSTANTS
155
 
C     
156
 
      INTEGER    NEXTERNAL
157
 
      PARAMETER (NEXTERNAL=5)
158
 
      REAL*16 ZERO
159
 
      PARAMETER (ZERO=0.0E+00_16)
160
 
      REAL*16 THRS_CLOSE
161
 
      PARAMETER (THRS_CLOSE=1.0E-02_16)
162
 
C     
163
 
C     ARGUMENTS 
164
 
C     
165
 
      REAL*16 P(0:3,NEXTERNAL), NEWP(0:3,NEXTERNAL)
166
 
      LOGICAL ML5_0_MP_IS_CLOSE
167
 
      INTEGER WARNED
168
 
C     
169
 
C     LOCAL VARIABLES 
170
 
C     
171
 
      INTEGER I,J
172
 
      REAL*16 REF,REF2
173
 
      DOUBLE PRECISION BUFFDP
174
 
 
175
 
C     NOW MAKE SURE THE SHIFTED POINT IS NOT TOO FAR FROM THE ORIGINAL
176
 
C      ONE
177
 
      ML5_0_MP_IS_CLOSE = .TRUE.
178
 
      REF  = ZERO
179
 
      REF2 = ZERO
180
 
      DO J=1,NEXTERNAL
181
 
        DO I=0,3
182
 
          REF2 = REF2 + ABS(P(I,J))
183
 
          REF = REF + ABS(P(I,J)-NEWP(I,J))
184
 
        ENDDO
185
 
      ENDDO
186
 
 
187
 
      IF ((REF/REF2).GT.THRS_CLOSE) THEN
188
 
        ML5_0_MP_IS_CLOSE = .FALSE.
189
 
        IF (WARNED.LT.20) THEN
190
 
          BUFFDP = (REF/REF2)
191
 
          WRITE(*,*) 'WARNING:: The improved PS point is too far from'
192
 
     $     //' the original one',BUFFDP
193
 
          WARNED=WARNED+1
194
 
        ENDIF
195
 
      ENDIF
196
 
 
197
 
      END
198
 
 
199
 
      FUNCTION ML5_0_MP_IS_PHYSICAL(P,WARNED)
200
 
      IMPLICIT NONE
201
 
C     
202
 
C     CONSTANTS
203
 
C     
204
 
      INTEGER    NEXTERNAL
205
 
      PARAMETER (NEXTERNAL=5)
206
 
      INTEGER    NINITIAL
207
 
      PARAMETER (NINITIAL=2)
208
 
      REAL*16 ZERO
209
 
      PARAMETER (ZERO=0.0E+00_16)
210
 
      REAL*16 MP__ZERO
211
 
      PARAMETER (MP__ZERO=ZERO)
212
 
      REAL*16 ONE
213
 
      PARAMETER (ONE=1.0E+00_16)
214
 
      REAL*16 TWO
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)
220
 
C     
221
 
C     ARGUMENTS 
222
 
C     
223
 
      REAL*16 P(0:3,NEXTERNAL)
224
 
      LOGICAL ML5_0_MP_IS_PHYSICAL
225
 
      INTEGER WARNED
226
 
C     
227
 
C     LOCAL VARIABLES 
228
 
C     
229
 
      INTEGER I,J
230
 
      REAL*16 BUFF,REF
231
 
      REAL*16 MASSES(NEXTERNAL)
232
 
      DOUBLE PRECISION BUFFDPA,BUFFDPB
233
 
C     
234
 
C     GLOBAL VARIABLES
235
 
C     
236
 
      INCLUDE 'mp_coupl.inc'
237
 
 
238
 
      MASSES(1)=MP__ZERO
239
 
      MASSES(2)=MP__ZERO
240
 
      MASSES(3)=MP__MDL_MW
241
 
      MASSES(4)=MP__MDL_MT
242
 
      MASSES(5)=MP__MDL_MB
243
 
 
244
 
C     ----------
245
 
C     BEGIN CODE
246
 
C     ----------
247
 
 
248
 
      ML5_0_MP_IS_PHYSICAL = .TRUE.
249
 
 
250
 
C     WE FIRST CHECK THAT THE INPUT PS POINT IS REASONABLY PHYSICAL
251
 
C     FOR THAT WE NEED A REFERENCE SCALE
252
 
      REF=ZERO
253
 
      DO J=1,NEXTERNAL
254
 
        REF=REF+ABS(P(0,J))
255
 
      ENDDO
256
 
      DO I=0,3
257
 
        BUFF=ZERO
258
 
        DO J=1,NINITIAL
259
 
          BUFF=BUFF-P(I,J)
260
 
        ENDDO
261
 
        DO J=NINITIAL+1,NEXTERNAL
262
 
          BUFF=BUFF+P(I,J)
263
 
        ENDDO
264
 
        IF ((BUFF/REF).GT.THRES_FOURMOM) THEN
265
 
          IF (WARNED.LT.20) THEN
266
 
            BUFFDPA = (BUFF/REF)
267
 
            WRITE(*,*) 'ERROR:: Four-momentum conservation is not'
268
 
     $       //' accurate enough, ',BUFFDPA
269
 
            CALL ML5_0_MP_WRITE_MOM(P)
270
 
            WARNED=WARNED+1
271
 
          ENDIF
272
 
          ML5_0_MP_IS_PHYSICAL = .FALSE.
273
 
        ENDIF
274
 
      ENDDO
275
 
      REF = REF / (ONE*NEXTERNAL)
276
 
      DO I=1,NEXTERNAL
277
 
        REF=ABS(P(0,I))+ABS(P(1,I))+ABS(P(2,I))+ABS(P(3,I))
278
 
        IF ((SQRT(ABS(P(0,I)**2-P(1,I)**2-P(2,I)**2-P(3,I)**2-MASSES(I)
279
 
     $   **2))/REF).GT.THRES_ONSHELL) THEN
280
 
          IF (WARNED.LT.20) THEN
281
 
            BUFFDPA=MASSES(I)
282
 
            BUFFDPB=(SQRT(ABS(P(0,I)**2-P(1,I)**2-P(2,I)**2-P(3,I)**2
283
 
     $       -MASSES(I)**2))/REF)
284
 
            WRITE(*,*) 'ERROR:: Onshellness of the momentum of'
285
 
     $       //' particle ',I,' of mass ',BUFFDPA,' is not accurate'
286
 
     $       //' enough, ',BUFFDPB
287
 
            CALL ML5_0_MP_WRITE_MOM(P)
288
 
            WARNED=WARNED+1
289
 
          ENDIF
290
 
          ML5_0_MP_IS_PHYSICAL = .FALSE.
291
 
        ENDIF
292
 
      ENDDO
293
 
 
294
 
      END
295
 
 
296
 
      SUBROUTINE ML5_0_WRITE_MOM(P)
297
 
      IMPLICIT NONE
298
 
      INTEGER    NEXTERNAL
299
 
      PARAMETER (NEXTERNAL=5)
300
 
      INTEGER    NINITIAL
301
 
      PARAMETER (NINITIAL=2)
302
 
      DOUBLE PRECISION ZERO
303
 
      PARAMETER (ZERO=0.0D0)
304
 
      DOUBLE PRECISION ML5_0_MDOT
305
 
 
306
 
      INTEGER I,J
307
 
 
308
 
C     
309
 
C     ARGUMENTS 
310
 
C     
311
 
      DOUBLE PRECISION P(0:3,NEXTERNAL),PSUM(0:3)
312
 
      DO I=0,3
313
 
        PSUM(I)=ZERO
314
 
        DO J=1,NINITIAL
315
 
          PSUM(I)=PSUM(I)+P(I,J)
316
 
        ENDDO
317
 
        DO J=NINITIAL+1,NEXTERNAL
318
 
          PSUM(I)=PSUM(I)-P(I,J)
319
 
        ENDDO
320
 
      ENDDO
321
 
      WRITE (*,*) ' Phase space point:'
322
 
      WRITE (*,*) '    ---------------------'
323
 
      WRITE (*,*) '    E | px | py | pz | m '
324
 
      DO I=1,NEXTERNAL
325
 
        WRITE (*,'(1x,5e27.17)') P(0,I),P(1,I),P(2,I),P(3,I)
326
 
     $   ,SQRT(ABS(ML5_0_MDOT(P(0,I),P(0,I))))
327
 
      ENDDO
328
 
      WRITE (*,*) '    Four-momentum conservation sum:'
329
 
      WRITE (*,'(1x,4e27.17)') PSUM(0),PSUM(1),PSUM(2),PSUM(3)
330
 
      WRITE (*,*) '   ---------------------'
331
 
      END
332
 
 
333
 
      DOUBLE PRECISION FUNCTION ML5_0_MDOT(P1,P2)
334
 
      IMPLICIT NONE
335
 
      DOUBLE PRECISION P1(0:3),P2(0:3)
336
 
      ML5_0_MDOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
337
 
      RETURN
338
 
      END
339
 
 
340
 
      SUBROUTINE ML5_0_MP_WRITE_MOM(P)
341
 
      IMPLICIT NONE
342
 
      INTEGER    NEXTERNAL
343
 
      PARAMETER (NEXTERNAL=5)
344
 
      INTEGER    NINITIAL
345
 
      PARAMETER (NINITIAL=2)
346
 
      REAL*16 ZERO
347
 
      PARAMETER (ZERO=0.0E+00_16)
348
 
      REAL*16 ML5_0_MP_MDOT
349
 
 
350
 
      INTEGER I,J
351
 
 
352
 
C     
353
 
C     ARGUMENTS 
354
 
C     
355
 
      REAL*16 P(0:3,NEXTERNAL),PSUM(0:3),DOT
356
 
      DOUBLE PRECISION DP_P(0:3,NEXTERNAL),DP_PSUM(0:3),DP_DOT
357
 
 
358
 
      DO I=0,3
359
 
        PSUM(I)=ZERO
360
 
        DO J=1,NINITIAL
361
 
          PSUM(I)=PSUM(I)+P(I,J)
362
 
        ENDDO
363
 
        DO J=NINITIAL+1,NEXTERNAL
364
 
          PSUM(I)=PSUM(I)-P(I,J)
365
 
        ENDDO
366
 
      ENDDO
367
 
 
368
 
C     The GCC4.7 compiler on SLC machines has trouble to write out
369
 
C      quadruple precision variable with the write(*,*) statement. I
370
 
C      therefore perform the cast by hand
371
 
      DO I=0,3
372
 
        DP_PSUM(I)=PSUM(I)
373
 
        DO J=1,NEXTERNAL
374
 
          DP_P(I,J)=P(I,J)
375
 
        ENDDO
376
 
      ENDDO
377
 
 
378
 
      WRITE (*,*) ' Phase space point:'
379
 
      WRITE (*,*) '    ---------------------'
380
 
      WRITE (*,*) '    E | px | py | pz | m '
381
 
      DO I=1,NEXTERNAL
382
 
        DOT=SQRT(ABS(ML5_0_MP_MDOT(P(0,I),P(0,I))))
383
 
        DP_DOT=DOT
384
 
        WRITE (*,'(1x,5e27.17)') DP_P(0,I),DP_P(1,I),DP_P(2,I),DP_P(3
385
 
     $   ,I),DP_DOT
386
 
      ENDDO
387
 
      WRITE (*,*) '    Four-momentum conservation sum:'
388
 
      WRITE (*,'(1x,4e27.17)') DP_PSUM(0),DP_PSUM(1),DP_PSUM(2)
389
 
     $ ,DP_PSUM(3)
390
 
      WRITE (*,*) '   ---------------------'
391
 
      END
392
 
 
393
 
      REAL*16 FUNCTION ML5_0_MP_MDOT(P1,P2)
394
 
      IMPLICIT NONE
395
 
      REAL*16 P1(0:3),P2(0:3)
396
 
      ML5_0_MP_MDOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
397
 
      RETURN
398
 
      END
399
 
 
400
 
C     Rotate_PS rotates the PS point PS (without modifying it)
401
 
C     stores the result in P and for the quadruple precision 
402
 
C     version , it also modifies the global variables
403
 
C     PS and MP_DONE accordingly.
404
 
 
405
 
      SUBROUTINE ML5_0_ROTATE_PS(P_IN,P,ROTATION)
406
 
      IMPLICIT NONE
407
 
C     
408
 
C     CONSTANTS 
409
 
C     
410
 
      INTEGER    NEXTERNAL
411
 
      PARAMETER (NEXTERNAL=5)
412
 
C     
413
 
C     ARGUMENTS 
414
 
C     
415
 
      DOUBLE PRECISION P_IN(0:3,NEXTERNAL),P(0:3,NEXTERNAL)
416
 
      INTEGER ROTATION
417
 
C     
418
 
C     LOCAL VARIABLES 
419
 
C     
420
 
      INTEGER I,J
421
 
 
422
 
C     ----------
423
 
C     BEGIN CODE
424
 
C     ----------
425
 
 
426
 
      DO I=1,NEXTERNAL
427
 
C       rotation=1 => (xp=z,yp=-x,zp=-y)
428
 
        IF(ROTATION.EQ.1) THEN
429
 
          P(0,I)=P_IN(0,I)
430
 
          P(1,I)=P_IN(3,I)
431
 
          P(2,I)=-P_IN(1,I)
432
 
          P(3,I)=-P_IN(2,I)
433
 
C         rotation=2 => (xp=-z,yp=y,zp=x)
434
 
        ELSEIF(ROTATION.EQ.2) THEN
435
 
          P(0,I)=P_IN(0,I)
436
 
          P(1,I)=-P_IN(3,I)
437
 
          P(2,I)=P_IN(2,I)
438
 
          P(3,I)=P_IN(1,I)
439
 
        ELSE
440
 
          P(0,I)=P_IN(0,I)
441
 
          P(1,I)=P_IN(1,I)
442
 
          P(2,I)=P_IN(2,I)
443
 
          P(3,I)=P_IN(3,I)
444
 
        ENDIF
445
 
      ENDDO
446
 
 
447
 
      END
448
 
 
449
 
 
450
 
      SUBROUTINE ML5_0_MP_ROTATE_PS(P_IN,P,ROTATION)
451
 
      IMPLICIT NONE
452
 
C     
453
 
C     CONSTANTS 
454
 
C     
455
 
      INTEGER    NEXTERNAL
456
 
      PARAMETER (NEXTERNAL=5)
457
 
C     
458
 
C     ARGUMENTS 
459
 
C     
460
 
      REAL*16 P_IN(0:3,NEXTERNAL),P(0:3,NEXTERNAL)
461
 
      INTEGER ROTATION
462
 
C     
463
 
C     LOCAL VARIABLES 
464
 
C     
465
 
      INTEGER I,J
466
 
C     
467
 
C     GLOBAL VARIABLES
468
 
C     
469
 
      LOGICAL MP_DONE
470
 
      COMMON/ML5_0_MP_DONE/MP_DONE
471
 
 
472
 
C     ----------
473
 
C     BEGIN CODE
474
 
C     ----------
475
 
 
476
 
      DO I=1,NEXTERNAL
477
 
C       rotation=1 => (xp=z,yp=-x,zp=-y)
478
 
        IF(ROTATION.EQ.1) THEN
479
 
          P(0,I)=P_IN(0,I)
480
 
          P(1,I)=P_IN(3,I)
481
 
          P(2,I)=-P_IN(1,I)
482
 
          P(3,I)=-P_IN(2,I)
483
 
C         rotation=2 => (xp=-z,yp=y,zp=x)
484
 
        ELSEIF(ROTATION.EQ.2) THEN
485
 
          P(0,I)=P_IN(0,I)
486
 
          P(1,I)=-P_IN(3,I)
487
 
          P(2,I)=P_IN(2,I)
488
 
          P(3,I)=P_IN(1,I)
489
 
        ELSE
490
 
          P(0,I)=P_IN(0,I)
491
 
          P(1,I)=P_IN(1,I)
492
 
          P(2,I)=P_IN(2,I)
493
 
          P(3,I)=P_IN(3,I)
494
 
        ENDIF
495
 
      ENDDO
496
 
 
497
 
      MP_DONE = .FALSE.
498
 
 
499
 
      END
500
 
 
501
 
C     *****************************************************************
502
 
C     Beginning of the routine for restoring precision with V.H. method
503
 
C     *****************************************************************
504
 
 
505
 
      SUBROUTINE ML5_0_MP_ORIG_IMPROVE_PS_POINT_PRECISION(P,ERRCODE
506
 
     $ ,WARNED)
507
 
      IMPLICIT NONE
508
 
C     
509
 
C     CONSTANTS 
510
 
C     
511
 
      INTEGER    NEXTERNAL
512
 
      PARAMETER (NEXTERNAL=5)
513
 
      INTEGER    NINITIAL
514
 
      PARAMETER (NINITIAL=2)
515
 
      REAL*16 ZERO
516
 
      PARAMETER (ZERO=0.0E+00_16)
517
 
      REAL*16 MP__ZERO
518
 
      PARAMETER (MP__ZERO=ZERO)
519
 
      REAL*16 ONE
520
 
      PARAMETER (ONE=1.0E+00_16)
521
 
      REAL*16 TWO
522
 
      PARAMETER (TWO=2.0E+00_16)
523
 
      REAL*16 THRS_TEST
524
 
      PARAMETER (THRS_TEST=1.0E-15_16)
525
 
C     
526
 
C     ARGUMENTS 
527
 
C     
528
 
      REAL*16 P(0:3,NEXTERNAL)
529
 
      INTEGER ERRCODE, WARNED
530
 
C     
531
 
C     FUNCTIONS
532
 
C     
533
 
      LOGICAL ML5_0_MP_IS_CLOSE
534
 
C     
535
 
C     LOCAL VARIABLES 
536
 
C     
537
 
      INTEGER I,J, P1, P2
538
 
C     PT STANDS FOR PTOT
539
 
      REAL*16 PT(0:3), NEWP(0:3,NEXTERNAL)
540
 
      REAL*16 BUFF,REF,REF2,DISCR
541
 
      REAL*16 MASSES(NEXTERNAL)
542
 
      REAL*16 SHIFTE(2),SHIFTZ(2)
543
 
C     
544
 
C     GLOBAL VARIABLES
545
 
C     
546
 
      INCLUDE 'mp_coupl.inc'
547
 
 
548
 
      MASSES(1)=MP__ZERO
549
 
      MASSES(2)=MP__ZERO
550
 
      MASSES(3)=MP__MDL_MW
551
 
      MASSES(4)=MP__MDL_MT
552
 
      MASSES(5)=MP__MDL_MB
553
 
 
554
 
C     ----------
555
 
C     BEGIN CODE
556
 
C     ----------
557
 
      ERRCODE = 0
558
 
 
559
 
C     NOW WE MAKE SURE THAT THE PS POINT CAN BE IMPROVED BY THE
560
 
C      ALGORITHM
561
 
      REF=ZERO
562
 
      DO J=1,NEXTERNAL
563
 
        REF=REF+ABS(P(0,J))
564
 
      ENDDO
565
 
 
566
 
      IF (NINITIAL.NE.2) ERRCODE = 100
567
 
 
568
 
      IF (ABS(P(1,1)/REF).GT.THRS_TEST.OR.ABS(P(2,1)/REF)
569
 
     $ .GT.THRS_TEST.OR.ABS(P(1,2)/REF).GT.THRS_TEST.OR.ABS(P(2,2)/REF)
570
 
     $ .GT.THRS_TEST) ERRCODE = 200
571
 
 
572
 
      IF (MASSES(1).NE.ZERO.OR.MASSES(2).NE.ZERO) ERRCODE = 300
573
 
 
574
 
      DO I=1,NEXTERNAL
575
 
        IF (P(0,I).LT.ZERO) ERRCODE = 400 + I
576
 
      ENDDO
577
 
 
578
 
      IF (ERRCODE.NE.0) GOTO 100
579
 
 
580
 
C     WE FIRST SHIFT ALL THE FINAL STATE PARTICLES TO MAKE THEM
581
 
C      EXACTLY ONSHELL
582
 
 
583
 
      DO I=0,3
584
 
        PT(I)=ZERO
585
 
      ENDDO
586
 
      DO I=NINITIAL+1,NEXTERNAL
587
 
        DO J=0,3
588
 
          IF (J.EQ.3) THEN
589
 
            NEWP(3,I)=SIGN(SQRT(ABS(P(0,I)**2-P(1,I)**2-P(2,I)**2
590
 
     $       -MASSES(I)**2)),P(3,I))
591
 
          ELSE
592
 
            NEWP(J,I)=P(J,I)
593
 
          ENDIF
594
 
          PT(J)=PT(J)+NEWP(J,I)
595
 
        ENDDO
596
 
      ENDDO
597
 
 
598
 
C     WE CHOOSE P1 IN THE ALGORITHM TO ALWAYS BE THE PARTICLE WITH
599
 
C      POSITIVE PZ
600
 
      IF (P(3,1).GT.ZERO) THEN
601
 
        P1=1
602
 
        P2=2
603
 
      ELSEIF (P(3,2).GT.ZERO) THEN
604
 
        P1=2
605
 
        P2=1
606
 
      ELSE
607
 
        ERRCODE = 500
608
 
        GOTO 100
609
 
      ENDIF
610
 
 
611
 
C     Now we calculate the shift to bring to P1 and P2
612
 
C     Mathematica gives
613
 
C     ptotC = {ptotE, ptotX, ptotY, ptotZ};
614
 
C     pm1C = {pm1E + sm1E, pm1X, pm1Y, pm1Z + sm1Z};
615
 
C     {pm0E + sm0E, ptotX - pm1X, ptotY - pm1Y, pm0Z + sm0Z};
616
 
C     sol = Solve[{ptotC[[1]] - pm1C[[1]] - pm0C[[1]] == 0,  
617
 
C     ptotC[[4]] - pm1C[[4]] - pm0C[[4]] == 0,
618
 
C     pm1C[[1]]^2 - pm1C[[2]]^2 - pm1C[[3]]^2 - pm1C[[4]]^2 == m1M^2,
619
 
C     pm0C[[1]]^2 - pm0C[[2]]^2 - pm0C[[3]]^2 - pm0C[[4]]^2 == m2M^2},
620
 
C     {sm1E, sm1Z, sm0E, sm0Z}] // FullSimplify;
621
 
C     (solC[[1]] /. {m1M -> 0, m2M -> 0} /. {pm1X -> 0, pm1Y -> 0})
622
 
C     END
623
 
C     
624
 
      DISCR = -PT(0)**2 + PT(1)**2 + PT(2)**2 + PT(3)**2
625
 
      IF (DISCR.LT.ZERO) DISCR = -DISCR
626
 
 
627
 
      SHIFTE(1) = (PT(0)*(-TWO*P(0,P1)*PT(0) + PT(0)**2 + PT(1)**2 +
628
 
     $  PT(2)**2) + (TWO*P(0,P1) - PT(0))*PT(3)**2 + PT(3)*DISCR)/(TWO
629
 
     $ *(PT(0) - PT(3))*(PT(0) + PT(3)))
630
 
      SHIFTE(2) = -(PT(0)*(TWO*P(0,P2)*PT(0) - PT(0)**2 + PT(1)**2 +
631
 
     $  PT(2)**2) + (-TWO*P(0,P2) + PT(0))*PT(3)**2 + PT(3)*DISCR)
632
 
     $ /(TWO*(PT(0) - PT(3))*(PT(0) + PT(3)))
633
 
      SHIFTZ(1) = (-TWO*P(3,P1)*(PT(0)**2 - PT(3)**2) + PT(3)*(PT(0)*
634
 
     $ *2 + PT(1)**2 + PT(2)**2 - PT(3)**2) + PT(0)*DISCR)/(TWO*(PT(0)
635
 
     $ **2 - PT(3)**2))
636
 
      SHIFTZ(2) = -(TWO*P(3,P2)*(PT(0)**2 - PT(3)**2) + PT(3)*(-PT(0)*
637
 
     $ *2 + PT(1)**2 + PT(2)**2 + PT(3)**2) + PT(0)*DISCR)/(TWO*(PT(0)
638
 
     $ **2 - PT(3)**2))
639
 
      NEWP(0,P1) = P(0,P1)+SHIFTE(1)
640
 
      NEWP(3,P1) = P(3,P1)+SHIFTZ(1)
641
 
      NEWP(0,P2) = P(0,P2)+SHIFTE(2)
642
 
      NEWP(3,P2) = P(3,P2)+SHIFTZ(2)
643
 
      NEWP(1,P2) = P(1,P2)
644
 
      NEWP(2,P2) = P(2,P2)
645
 
      DO J=1,2
646
 
        REF=ZERO
647
 
        DO I=NINITIAL+1,NEXTERNAL
648
 
          REF = REF + P(J,I)
649
 
        ENDDO
650
 
        REF = REF - P(J,P2)
651
 
        NEWP(J,P1) = REF
652
 
      ENDDO
653
 
 
654
 
      IF (.NOT.ML5_0_MP_IS_CLOSE(P,NEWP,WARNED)) THEN
655
 
        ERRCODE=999
656
 
        GOTO 100
657
 
      ENDIF
658
 
 
659
 
      DO J=1,NEXTERNAL
660
 
        DO I=0,3
661
 
          P(I,J)=NEWP(I,J)
662
 
        ENDDO
663
 
      ENDDO
664
 
 
665
 
 100  CONTINUE
666
 
 
667
 
      END
668
 
 
669
 
C     *****************************************************************
670
 
C     Beginning of the routine for restoring precision a la PSMC
671
 
C     *****************************************************************
672
 
 
673
 
      SUBROUTINE ML5_0_MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE
674
 
     $ ,WARNED)
675
 
      IMPLICIT NONE
676
 
C     
677
 
C     CONSTANTS 
678
 
C     
679
 
      INTEGER    NEXTERNAL
680
 
      PARAMETER (NEXTERNAL=5)
681
 
      INTEGER    NINITIAL
682
 
      PARAMETER (NINITIAL=2)
683
 
      REAL*16 ZERO
684
 
      PARAMETER (ZERO=0.0E+00_16)
685
 
      REAL*16 MP__ZERO
686
 
      PARAMETER (MP__ZERO=ZERO)
687
 
      REAL*16 ONE
688
 
      PARAMETER (ONE=1.0E+00_16)
689
 
      REAL*16 TWO
690
 
      PARAMETER (TWO=2.0E+00_16)
691
 
      REAL*16 CONSISTENCY_THRES
692
 
      PARAMETER (CONSISTENCY_THRES=1.0E-25_16)
693
 
 
694
 
      INTEGER    NAPPROXZEROS
695
 
      PARAMETER (NAPPROXZEROS=3)
696
 
 
697
 
C     
698
 
C     ARGUMENTS 
699
 
C     
700
 
      REAL*16 P(0:3,NEXTERNAL)
701
 
      INTEGER ERRCODE,ERROR,WARNED
702
 
C     
703
 
C     FUNCTIONS
704
 
C     
705
 
      LOGICAL ML5_0_MP_IS_CLOSE
706
 
C     
707
 
C     LOCAL VARIABLES 
708
 
C     
709
 
      INTEGER I,J, P1, P2
710
 
      REAL*16 NEWP(0:3,NEXTERNAL), PBUFF(0:3)
711
 
      REAL*16 BUFF, BUFF2, XSCALE, APPROX_ZEROS(NAPPROXZEROS)
712
 
      REAL*16 MASSES(NEXTERNAL)
713
 
C     
714
 
C     GLOBAL VARIABLES
715
 
C     
716
 
      INCLUDE 'mp_coupl.inc'
717
 
 
718
 
C     ----------
719
 
C     BEGIN CODE
720
 
C     ----------
721
 
 
722
 
      MASSES(1)=MP__ZERO
723
 
      MASSES(2)=MP__ZERO
724
 
      MASSES(3)=MP__MDL_MW
725
 
      MASSES(4)=MP__MDL_MT
726
 
      MASSES(5)=MP__MDL_MB
727
 
 
728
 
      ERRCODE = 0
729
 
      XSCALE = ONE
730
 
 
731
 
C     Define the seeds which should be tried
732
 
      APPROX_ZEROS(1)=1.0E+00_16
733
 
      APPROX_ZEROS(2)=1.1E+00_16
734
 
      APPROX_ZEROS(3)=0.9E+00_16
735
 
 
736
 
C     Start by copying the momenta
737
 
      DO I=1,NEXTERNAL
738
 
        DO J=0,3
739
 
          NEWP(J,I)=P(J,I)
740
 
        ENDDO
741
 
      ENDDO
742
 
 
743
 
C     First make sur that the space like momentum is exactly conserved
744
 
      DO J=0,3
745
 
        PBUFF(J)=ZERO
746
 
      ENDDO
747
 
      DO I=1,NINITIAL
748
 
        DO J=1,3
749
 
          PBUFF(J)=PBUFF(J)+NEWP(J,I)
750
 
        ENDDO
751
 
      ENDDO
752
 
      DO I=NINITIAL+1,NEXTERNAL-1
753
 
        DO J=1,3
754
 
          PBUFF(J)=PBUFF(J)-NEWP(J,I)
755
 
        ENDDO
756
 
      ENDDO
757
 
      DO J=1,3
758
 
        NEWP(J,NEXTERNAL)=PBUFF(J)
759
 
      ENDDO
760
 
 
761
 
C     Now find the 'x' rescaling factor
762
 
      DO I=1,NAPPROXZEROS
763
 
        CALL ML5_0_FINDX(NEWP,APPROX_ZEROS(I),XSCALE,ERROR)
764
 
        IF(ERROR.EQ.0) THEN
765
 
          GOTO 1001
766
 
        ELSE
767
 
          ERRCODE=ERRCODE+(10**(I-1))*ERROR
768
 
        ENDIF
769
 
      ENDDO
770
 
      IF (WARNED.LT.20) THEN
771
 
        WRITE(*,*) 'WARNING:: Could not find the proper rescaling'
772
 
     $   //' factor x. Restoring precision ala PSMC will therefore not'
773
 
     $   //' be used.'
774
 
        WARNED=WARNED+1
775
 
      ENDIF
776
 
      IF (ERRCODE.LT.1000) THEN
777
 
        ERRCODE=ERRCODE+1000
778
 
      ENDIF
779
 
      GOTO 1000
780
 
 1001 CONTINUE
781
 
      ERRCODE = 0
782
 
 
783
 
C     Apply the rescaling
784
 
      DO I=1,NEXTERNAL
785
 
        DO J=1,3
786
 
          NEWP(J,I)=NEWP(J,I)*XSCALE
787
 
        ENDDO
788
 
      ENDDO
789
 
 
790
 
C     Now restore exact onshellness of the particles.
791
 
      DO I=1,NEXTERNAL
792
 
        BUFF=MASSES(I)**2
793
 
        DO J=1,3
794
 
          BUFF=BUFF+NEWP(J,I)**2
795
 
        ENDDO
796
 
        NEWP(0,I)=SQRT(BUFF)
797
 
      ENDDO
798
 
 
799
 
C     Consistency check
800
 
      BUFF=ZERO
801
 
      BUFF2=ZERO
802
 
      DO I=1,NINITIAL
803
 
        BUFF=BUFF-NEWP(0,I)
804
 
        BUFF2=BUFF2+NEWP(0,I)
805
 
      ENDDO
806
 
      DO I=NINITIAL+1,NEXTERNAL
807
 
        BUFF=BUFF+NEWP(0,I)
808
 
        BUFF2=BUFF2+NEWP(0,I)
809
 
      ENDDO
810
 
      IF ((ABS(BUFF)/BUFF2).GT.CONSISTENCY_THRES) THEN
811
 
        IF (WARNED.LT.20) THEN
812
 
          WRITE(*,*) 'WARNING:: The consistency check in the a la PSMC'
813
 
     $     //' precision restoring algorithm failed. The result will'
814
 
     $     //' therefore not be used.'
815
 
          WARNED=WARNED+1
816
 
        ENDIF
817
 
        ERRCODE = 1000
818
 
        GOTO 1000
819
 
      ENDIF
820
 
 
821
 
      IF (.NOT.ML5_0_MP_IS_CLOSE(P,NEWP,WARNED)) THEN
822
 
        ERRCODE=999
823
 
        GOTO 1000
824
 
      ENDIF
825
 
 
826
 
      DO J=1,NEXTERNAL
827
 
        DO I=0,3
828
 
          P(I,J)=NEWP(I,J)
829
 
        ENDDO
830
 
      ENDDO
831
 
 
832
 
 1000 CONTINUE
833
 
 
834
 
      END
835
 
 
836
 
 
837
 
      SUBROUTINE ML5_0_FINDX(P,SEED,XSCALE,ERROR)
838
 
      IMPLICIT NONE
839
 
C     
840
 
C     CONSTANTS 
841
 
C     
842
 
      INTEGER    NEXTERNAL
843
 
      PARAMETER (NEXTERNAL=5)
844
 
      INTEGER    NINITIAL
845
 
      PARAMETER (NINITIAL=2)
846
 
      REAL*16 ZERO
847
 
      PARAMETER (ZERO=0.0E+00_16)
848
 
      REAL*16 MP__ZERO
849
 
      PARAMETER (MP__ZERO=ZERO)
850
 
      REAL*16 ONE
851
 
      PARAMETER (ONE=1.0E+00_16)
852
 
      REAL*16 TWO
853
 
      PARAMETER (TWO=2.0E+00_16)
854
 
      INTEGER MAXITERATIONS
855
 
      PARAMETER (MAXITERATIONS=8)
856
 
      REAL*16 CONVERGED
857
 
      PARAMETER (CONVERGED=1.0E-26_16)
858
 
C     
859
 
C     ARGUMENTS 
860
 
C     
861
 
      REAL*16 P(0:3,NEXTERNAL),SEED,XSCALE
862
 
      INTEGER ERROR
863
 
C     
864
 
C     LOCAL VARIABLES 
865
 
C     
866
 
      INTEGER I,J,ERR
867
 
      REAL*16 PVECSQ(NEXTERNAL)
868
 
      REAL*16 XN, XNP1,FVAL,DVAL
869
 
 
870
 
C     ----------
871
 
C     BEGIN CODE
872
 
C     ----------
873
 
 
874
 
      ERROR = 0
875
 
      XSCALE = SEED
876
 
      XN = SEED
877
 
      XNP1 = SEED
878
 
 
879
 
      DO I=1,NEXTERNAL
880
 
        PVECSQ(I)=P(1,I)**2+P(2,I)**2+P(3,I)**2
881
 
      ENDDO
882
 
 
883
 
      DO I=1,MAXITERATIONS
884
 
        CALL ML5_0_FUNCT(PVECSQ(1),XN,.FALSE.,ERR, FVAL)
885
 
        IF (ERR.NE.0) THEN
886
 
          ERROR=ERR
887
 
          GOTO 710
888
 
        ENDIF
889
 
        CALL ML5_0_FUNCT(PVECSQ(1),XN,.TRUE.,ERR, DVAL)
890
 
        IF (ERR.NE.0) THEN
891
 
          ERROR=ERR
892
 
          GOTO 710
893
 
        ENDIF
894
 
        XNP1=XN-(FVAL/DVAL)
895
 
        IF((ABS(((XNP1-XN)*TWO)/(XNP1+XN))).LT.CONVERGED) THEN
896
 
          XN=XNP1
897
 
          GOTO 700
898
 
        ENDIF
899
 
        XN=XNP1
900
 
      ENDDO
901
 
      ERROR=9
902
 
      GOTO 710
903
 
 
904
 
 700  CONTINUE
905
 
C     For good measure, we iterate one last time
906
 
      CALL ML5_0_FUNCT(PVECSQ(1),XN,.FALSE.,ERR, FVAL)
907
 
      IF (ERR.NE.0) THEN
908
 
        ERROR=ERR
909
 
        GOTO 710
910
 
      ENDIF
911
 
      CALL ML5_0_FUNCT(PVECSQ(1),XN,.TRUE.,ERR, DVAL)
912
 
      IF (ERR.NE.0) THEN
913
 
        ERROR=ERR
914
 
        GOTO 710
915
 
      ENDIF
916
 
 
917
 
      XSCALE=XN-(FVAL/DVAL)
918
 
 
919
 
 710  CONTINUE
920
 
 
921
 
      END
922
 
 
923
 
      SUBROUTINE ML5_0_FUNCT(PVECSQ,X,DERIVATIVE,ERROR,RES)
924
 
      IMPLICIT NONE
925
 
C     
926
 
C     CONSTANTS 
927
 
C     
928
 
      INTEGER    NEXTERNAL
929
 
      PARAMETER (NEXTERNAL=5)
930
 
      INTEGER    NINITIAL
931
 
      PARAMETER (NINITIAL=2)
932
 
      REAL*16 ZERO
933
 
      PARAMETER (ZERO=0.0E+00_16)
934
 
      REAL*16 MP__ZERO
935
 
      PARAMETER (MP__ZERO=ZERO)
936
 
      REAL*16 ONE
937
 
      PARAMETER (ONE=1.0E+00_16)
938
 
      REAL*16 TWO
939
 
      PARAMETER (TWO=2.0E+00_16)
940
 
C     
941
 
C     ARGUMENTS 
942
 
C     
943
 
      REAL*16 PVECSQ(NEXTERNAL),X,RES
944
 
      INTEGER ERROR
945
 
      LOGICAL DERIVATIVE
946
 
C     
947
 
C     LOCAL VARIABLES 
948
 
C     
949
 
      INTEGER I,J
950
 
      REAL*16 BUFF,FACTOR
951
 
      REAL*16 MASSES(NEXTERNAL)
952
 
C     
953
 
C     GLOBAL VARIABLES
954
 
C     
955
 
      INCLUDE 'mp_coupl.inc'
956
 
 
957
 
C     ----------
958
 
C     BEGIN CODE
959
 
C     ----------
960
 
 
961
 
      MASSES(1)=MP__ZERO
962
 
      MASSES(2)=MP__ZERO
963
 
      MASSES(3)=MP__MDL_MW
964
 
      MASSES(4)=MP__MDL_MT
965
 
      MASSES(5)=MP__MDL_MB
966
 
 
967
 
      ERROR=0
968
 
      RES=ZERO
969
 
      BUFF=ZERO
970
 
 
971
 
      DO I=1,NEXTERNAL
972
 
        IF (I.LE.NINITIAL) THEN
973
 
          FACTOR=-ONE
974
 
        ELSE
975
 
          FACTOR=ONE
976
 
        ENDIF
977
 
        BUFF=MASSES(I)**2+PVECSQ(I)*X**2
978
 
        IF (BUFF.LT.ZERO) THEN
979
 
          RES=ZERO
980
 
          ERROR = 1
981
 
          GOTO 800
982
 
        ENDIF
983
 
        IF (DERIVATIVE) THEN
984
 
          RES=RES + FACTOR*((X*PVECSQ(I))/SQRT(BUFF))
985
 
        ELSE
986
 
          RES=RES + FACTOR*SQRT(BUFF)
987
 
        ENDIF
988
 
      ENDDO
989
 
 
990
 
 800  CONTINUE
991
 
 
992
 
      END
993