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