← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

Viewing changes to tests/input_files/IOTestsComparison/TestCmdMatchBox.moved/MatchBoxOutput/%TEST%SubProcesses%P1_uux_uux%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/TestCmdMatchBox.moved/MatchBoxOutput/%TEST%SubProcesses%P1_uux_uux%improve_ps.f
 
1
      SUBROUTINE MG5_1_IMPROVE_PS_POINT_PRECISION(P)
 
2
      IMPLICIT NONE
 
3
C     
 
4
C     CONSTANTS
 
5
C     
 
6
      INTEGER    NEXTERNAL
 
7
      PARAMETER (NEXTERNAL=4)
 
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 MG5_1_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 MG5_1_MP_IMPROVE_PS_POINT_PRECISION(P)
 
40
      IMPLICIT NONE
 
41
C     
 
42
C     CONSTANTS
 
43
C     
 
44
      INTEGER    NEXTERNAL
 
45
      PARAMETER (NEXTERNAL=4)
 
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 MG5_1_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.MG5_1_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 MG5_1_MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE
 
102
     $   ,WARNED)
 
103
      ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN
 
104
        CALL MG5_1_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 precisio'
 
110
     $     //'n improvement with an alternative method.'
 
111
          WARNED=WARNED+1
 
112
        ENDIF
 
113
        IF (IMPROVEPSPOINT.EQ.1) THEN
 
114
          CALL MG5_1_MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODETMP
 
115
     $     ,WARNED)
 
116
        ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN
 
117
          CALL MG5_1_MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODETMP
 
118
     $     ,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 MG5_1_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_p'
 
144
     $   //'s routine will now be supressed.'
 
145
        TOLD_SUPPRESS=.TRUE.
 
146
      ENDIF
 
147
 
 
148
      END
 
149
 
 
150
 
 
151
      FUNCTION MG5_1_MP_IS_CLOSE(P,NEWP,WARNED)
 
152
      IMPLICIT NONE
 
153
C     
 
154
C     CONSTANTS
 
155
C     
 
156
      INTEGER    NEXTERNAL
 
157
      PARAMETER (NEXTERNAL=4)
 
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 MG5_1_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
      MG5_1_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
        MG5_1_MP_IS_CLOSE = .FALSE.
 
189
        IF (WARNED.LT.20) THEN
 
190
          BUFFDP = (REF/REF2)
 
191
          WRITE(*,*) 'WARNING:: The improved PS point is too far fro'
 
192
     $     //'m the original one',BUFFDP
 
193
          WARNED=WARNED+1
 
194
        ENDIF
 
195
      ENDIF
 
196
 
 
197
      END
 
198
 
 
199
      FUNCTION MG5_1_MP_IS_PHYSICAL(P,WARNED)
 
200
      IMPLICIT NONE
 
201
C     
 
202
C     CONSTANTS
 
203
C     
 
204
      INTEGER    NEXTERNAL
 
205
      PARAMETER (NEXTERNAL=4)
 
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 MG5_1_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__ZERO
 
241
      MASSES(4)=MP__ZERO
 
242
 
 
243
C     ----------
 
244
C     BEGIN CODE
 
245
C     ----------
 
246
 
 
247
      MG5_1_MP_IS_PHYSICAL = .TRUE.
 
248
 
 
249
C     WE FIRST CHECK THAT THE INPUT PS POINT IS REASONABLY PHYSICAL
 
250
C     FOR THAT WE NEED A REFERENCE SCALE
 
251
      REF=ZERO
 
252
      DO J=1,NEXTERNAL
 
253
        REF=REF+ABS(P(0,J))
 
254
      ENDDO
 
255
      DO I=0,3
 
256
        BUFF=ZERO
 
257
        DO J=1,NINITIAL
 
258
          BUFF=BUFF-P(I,J)
 
259
        ENDDO
 
260
        DO J=NINITIAL+1,NEXTERNAL
 
261
          BUFF=BUFF+P(I,J)
 
262
        ENDDO
 
263
        IF ((BUFF/REF).GT.THRES_FOURMOM) THEN
 
264
          IF (WARNED.LT.20) THEN
 
265
            BUFFDPA = (BUFF/REF)
 
266
            WRITE(*,*) 'ERROR:: Four-momentum conservation is no'
 
267
     $       //'t accurate enough, ',BUFFDPA
 
268
            CALL MG5_1_MP_WRITE_MOM(P)
 
269
            WARNED=WARNED+1
 
270
          ENDIF
 
271
          MG5_1_MP_IS_PHYSICAL = .FALSE.
 
272
        ENDIF
 
273
      ENDDO
 
274
      REF = REF / (ONE*NEXTERNAL)
 
275
      DO I=1,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
 
278
     $   -MASSES(I)**2))/REF).GT.THRES_ONSHELL) THEN
 
279
          IF (WARNED.LT.20) THEN
 
280
            BUFFDPA=MASSES(I)
 
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 o'
 
284
     $       //'f particle ',I,' of mass ',BUFFDPA,' is not accurat'
 
285
     $       //'e enough, ',BUFFDPB
 
286
            CALL MG5_1_MP_WRITE_MOM(P)
 
287
            WARNED=WARNED+1
 
288
          ENDIF
 
289
          MG5_1_MP_IS_PHYSICAL = .FALSE.
 
290
        ENDIF
 
291
      ENDDO
 
292
 
 
293
      END
 
294
 
 
295
      SUBROUTINE MG5_1_WRITE_MOM(P)
 
296
      IMPLICIT NONE
 
297
      INTEGER    NEXTERNAL
 
298
      PARAMETER (NEXTERNAL=4)
 
299
      INTEGER    NINITIAL
 
300
      PARAMETER (NINITIAL=2)
 
301
      DOUBLE PRECISION ZERO
 
302
      PARAMETER (ZERO=0.0D0)
 
303
      DOUBLE PRECISION MG5_1_MDOT
 
304
 
 
305
      INTEGER I,J
 
306
 
 
307
C     
 
308
C     ARGUMENTS 
 
309
C     
 
310
      DOUBLE PRECISION P(0:3,NEXTERNAL),PSUM(0:3)
 
311
      DO I=0,3
 
312
        PSUM(I)=ZERO
 
313
        DO J=1,NINITIAL
 
314
          PSUM(I)=PSUM(I)+P(I,J)
 
315
        ENDDO
 
316
        DO J=NINITIAL+1,NEXTERNAL
 
317
          PSUM(I)=PSUM(I)-P(I,J)
 
318
        ENDDO
 
319
      ENDDO
 
320
      WRITE (*,*) ' Phase space point:'
 
321
      WRITE (*,*) '    ---------------------'
 
322
      WRITE (*,*) '    E | px | py | pz | m '
 
323
      DO I=1,NEXTERNAL
 
324
        WRITE (*,'(1x,5e27.17)') P(0,I),P(1,I),P(2,I),P(3,I),SQRT(ABS(M
 
325
     $   G5_1_MDOT(P(0,I),P(0,I))))
 
326
      ENDDO
 
327
      WRITE (*,*) '    Four-momentum conservation sum:'
 
328
      WRITE (*,'(1x,4e27.17)') PSUM(0),PSUM(1),PSUM(2),PSUM(3)
 
329
      WRITE (*,*) '   ---------------------'
 
330
      END
 
331
 
 
332
      DOUBLE PRECISION FUNCTION MG5_1_MDOT(P1,P2)
 
333
      IMPLICIT NONE
 
334
      DOUBLE PRECISION P1(0:3),P2(0:3)
 
335
      MG5_1_MDOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
 
336
      RETURN
 
337
      END
 
338
 
 
339
      SUBROUTINE MG5_1_MP_WRITE_MOM(P)
 
340
      IMPLICIT NONE
 
341
      INTEGER    NEXTERNAL
 
342
      PARAMETER (NEXTERNAL=4)
 
343
      INTEGER    NINITIAL
 
344
      PARAMETER (NINITIAL=2)
 
345
      REAL*16 ZERO
 
346
      PARAMETER (ZERO=0.0E+00_16)
 
347
      REAL*16 MG5_1_MP_MDOT
 
348
 
 
349
      INTEGER I,J
 
350
 
 
351
C     
 
352
C     ARGUMENTS 
 
353
C     
 
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
 
356
 
 
357
      DO I=0,3
 
358
        PSUM(I)=ZERO
 
359
        DO J=1,NINITIAL
 
360
          PSUM(I)=PSUM(I)+P(I,J)
 
361
        ENDDO
 
362
        DO J=NINITIAL+1,NEXTERNAL
 
363
          PSUM(I)=PSUM(I)-P(I,J)
 
364
        ENDDO
 
365
      ENDDO
 
366
 
 
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
 
370
      DO I=0,3
 
371
        DP_PSUM(I)=PSUM(I)
 
372
        DO J=1,NEXTERNAL
 
373
          DP_P(I,J)=P(I,J)
 
374
        ENDDO
 
375
      ENDDO
 
376
 
 
377
      WRITE (*,*) ' Phase space point:'
 
378
      WRITE (*,*) '    ---------------------'
 
379
      WRITE (*,*) '    E | px | py | pz | m '
 
380
      DO I=1,NEXTERNAL
 
381
        DOT=SQRT(ABS(MG5_1_MP_MDOT(P(0,I),P(0,I))))
 
382
        DP_DOT=DOT
 
383
        WRITE (*,'(1x,5e27.17)') DP_P(0,I),DP_P(1,I),DP_P(2,I),DP_P(3
 
384
     $   ,I),DP_DOT
 
385
      ENDDO
 
386
      WRITE (*,*) '    Four-momentum conservation sum:'
 
387
      WRITE (*,'(1x,4e27.17)') DP_PSUM(0),DP_PSUM(1),DP_PSUM(2)
 
388
     $ ,DP_PSUM(3)
 
389
      WRITE (*,*) '   ---------------------'
 
390
      END
 
391
 
 
392
      REAL*16 FUNCTION MG5_1_MP_MDOT(P1,P2)
 
393
      IMPLICIT NONE
 
394
      REAL*16 P1(0:3),P2(0:3)
 
395
      MG5_1_MP_MDOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
 
396
      RETURN
 
397
      END
 
398
 
 
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.
 
403
 
 
404
      SUBROUTINE MG5_1_ROTATE_PS(P_IN,P,ROTATION)
 
405
      IMPLICIT NONE
 
406
C     
 
407
C     CONSTANTS 
 
408
C     
 
409
      INTEGER    NEXTERNAL
 
410
      PARAMETER (NEXTERNAL=4)
 
411
C     
 
412
C     ARGUMENTS 
 
413
C     
 
414
      DOUBLE PRECISION P_IN(0:3,NEXTERNAL),P(0:3,NEXTERNAL)
 
415
      INTEGER ROTATION
 
416
C     
 
417
C     LOCAL VARIABLES 
 
418
C     
 
419
      INTEGER I,J
 
420
 
 
421
C     ----------
 
422
C     BEGIN CODE
 
423
C     ----------
 
424
 
 
425
      DO I=1,NEXTERNAL
 
426
C       rotation=1 => (xp=z,yp=-x,zp=-y)
 
427
        IF(ROTATION.EQ.1) THEN
 
428
          P(0,I)=P_IN(0,I)
 
429
          P(1,I)=P_IN(3,I)
 
430
          P(2,I)=-P_IN(1,I)
 
431
          P(3,I)=-P_IN(2,I)
 
432
C         rotation=2 => (xp=-z,yp=y,zp=x)
 
433
        ELSEIF(ROTATION.EQ.2) THEN
 
434
          P(0,I)=P_IN(0,I)
 
435
          P(1,I)=-P_IN(3,I)
 
436
          P(2,I)=P_IN(2,I)
 
437
          P(3,I)=P_IN(1,I)
 
438
        ELSE
 
439
          P(0,I)=P_IN(0,I)
 
440
          P(1,I)=P_IN(1,I)
 
441
          P(2,I)=P_IN(2,I)
 
442
          P(3,I)=P_IN(3,I)
 
443
        ENDIF
 
444
      ENDDO
 
445
 
 
446
      END
 
447
 
 
448
 
 
449
      SUBROUTINE MG5_1_MP_ROTATE_PS(P_IN,P,ROTATION)
 
450
      IMPLICIT NONE
 
451
C     
 
452
C     CONSTANTS 
 
453
C     
 
454
      INTEGER    NEXTERNAL
 
455
      PARAMETER (NEXTERNAL=4)
 
456
C     
 
457
C     ARGUMENTS 
 
458
C     
 
459
      REAL*16 P_IN(0:3,NEXTERNAL),P(0:3,NEXTERNAL)
 
460
      INTEGER ROTATION
 
461
C     
 
462
C     LOCAL VARIABLES 
 
463
C     
 
464
      INTEGER I,J
 
465
C     
 
466
C     GLOBAL VARIABLES
 
467
C     
 
468
      LOGICAL MP_DONE
 
469
      COMMON/MG5_1_MP_DONE/MP_DONE
 
470
 
 
471
C     ----------
 
472
C     BEGIN CODE
 
473
C     ----------
 
474
 
 
475
      DO I=1,NEXTERNAL
 
476
C       rotation=1 => (xp=z,yp=-x,zp=-y)
 
477
        IF(ROTATION.EQ.1) THEN
 
478
          P(0,I)=P_IN(0,I)
 
479
          P(1,I)=P_IN(3,I)
 
480
          P(2,I)=-P_IN(1,I)
 
481
          P(3,I)=-P_IN(2,I)
 
482
C         rotation=2 => (xp=-z,yp=y,zp=x)
 
483
        ELSEIF(ROTATION.EQ.2) THEN
 
484
          P(0,I)=P_IN(0,I)
 
485
          P(1,I)=-P_IN(3,I)
 
486
          P(2,I)=P_IN(2,I)
 
487
          P(3,I)=P_IN(1,I)
 
488
        ELSE
 
489
          P(0,I)=P_IN(0,I)
 
490
          P(1,I)=P_IN(1,I)
 
491
          P(2,I)=P_IN(2,I)
 
492
          P(3,I)=P_IN(3,I)
 
493
        ENDIF
 
494
      ENDDO
 
495
 
 
496
      MP_DONE = .FALSE.
 
497
 
 
498
      END
 
499
 
 
500
C     *****************************************************************
 
501
C     Beginning of the routine for restoring precision with V.H. method
 
502
C     *****************************************************************
 
503
 
 
504
      SUBROUTINE MG5_1_MP_ORIG_IMPROVE_PS_POINT_PRECISION(P,ERRCODE
 
505
     $ ,WARNED)
 
506
      IMPLICIT NONE
 
507
C     
 
508
C     CONSTANTS 
 
509
C     
 
510
      INTEGER    NEXTERNAL
 
511
      PARAMETER (NEXTERNAL=4)
 
512
      INTEGER    NINITIAL
 
513
      PARAMETER (NINITIAL=2)
 
514
      REAL*16 ZERO
 
515
      PARAMETER (ZERO=0.0E+00_16)
 
516
      REAL*16 MP__ZERO
 
517
      PARAMETER (MP__ZERO=ZERO)
 
518
      REAL*16 ONE
 
519
      PARAMETER (ONE=1.0E+00_16)
 
520
      REAL*16 TWO
 
521
      PARAMETER (TWO=2.0E+00_16)
 
522
      REAL*16 THRS_TEST
 
523
      PARAMETER (THRS_TEST=1.0E-15_16)
 
524
C     
 
525
C     ARGUMENTS 
 
526
C     
 
527
      REAL*16 P(0:3,NEXTERNAL)
 
528
      INTEGER ERRCODE, WARNED
 
529
C     
 
530
C     FUNCTIONS
 
531
C     
 
532
      LOGICAL MG5_1_MP_IS_CLOSE
 
533
C     
 
534
C     LOCAL VARIABLES 
 
535
C     
 
536
      INTEGER I,J, P1, P2
 
537
C     PT STANDS FOR PTOT
 
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)
 
542
C     
 
543
C     GLOBAL VARIABLES
 
544
C     
 
545
      INCLUDE 'mp_coupl.inc'
 
546
 
 
547
      MASSES(1)=MP__ZERO
 
548
      MASSES(2)=MP__ZERO
 
549
      MASSES(3)=MP__ZERO
 
550
      MASSES(4)=MP__ZERO
 
551
 
 
552
C     ----------
 
553
C     BEGIN CODE
 
554
C     ----------
 
555
      ERRCODE = 0
 
556
 
 
557
C     NOW WE MAKE SURE THAT THE PS POINT CAN BE IMPROVED BY THE
 
558
C      ALGORITHM
 
559
      REF=ZERO
 
560
      DO J=1,NEXTERNAL
 
561
        REF=REF+ABS(P(0,J))
 
562
      ENDDO
 
563
 
 
564
      IF (NINITIAL.NE.2) ERRCODE = 100
 
565
 
 
566
      IF (ABS(P(1,1)/REF).GT.THRS_TEST.OR.ABS(P(2,1)/REF).GT.THRS_TEST.
 
567
     $ OR.ABS(P(1,2)/REF).GT.THRS_TEST.OR.ABS(P(2,2)/REF).GT.THRS_TEST
 
568
     $ ) ERRCODE = 200
 
569
 
 
570
      IF (MASSES(1).NE.ZERO.OR.MASSES(2).NE.ZERO) ERRCODE = 300
 
571
 
 
572
      DO I=1,NEXTERNAL
 
573
        IF (P(0,I).LT.ZERO) ERRCODE = 400 + I
 
574
      ENDDO
 
575
 
 
576
      IF (ERRCODE.NE.0) GOTO 100
 
577
 
 
578
C     WE FIRST SHIFT ALL THE FINAL STATE PARTICLES TO MAKE THEM
 
579
C      EXACTLY ONSHELL
 
580
 
 
581
      DO I=0,3
 
582
        PT(I)=ZERO
 
583
      ENDDO
 
584
      DO I=NINITIAL+1,NEXTERNAL
 
585
        DO J=0,3
 
586
          IF (J.EQ.3) THEN
 
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))
 
589
          ELSE
 
590
            NEWP(J,I)=P(J,I)
 
591
          ENDIF
 
592
          PT(J)=PT(J)+NEWP(J,I)
 
593
        ENDDO
 
594
      ENDDO
 
595
 
 
596
C     WE CHOOSE P1 IN THE ALGORITHM TO ALWAYS BE THE PARTICLE WITH
 
597
C      POSITIVE PZ
 
598
      IF (P(3,1).GT.ZERO) THEN
 
599
        P1=1
 
600
        P2=2
 
601
      ELSEIF (P(3,2).GT.ZERO) THEN
 
602
        P1=2
 
603
        P2=1
 
604
      ELSE
 
605
        ERRCODE = 500
 
606
        GOTO 100
 
607
      ENDIF
 
608
 
 
609
C     Now we calculate the shift to bring to P1 and P2
 
610
C     Mathematica gives
 
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})
 
620
C     END
 
621
C     
 
622
      DISCR = -PT(0)**2 + PT(1)**2 + PT(2)**2 + PT(3)**2
 
623
      IF (DISCR.LT.ZERO) DISCR = -DISCR
 
624
 
 
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)
 
627
     $ /(TWO*(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)
 
633
     $ **2 - PT(3)**2))
 
634
      SHIFTZ(2) = -(TWO*P(3,P2)*(PT(0)**2 - PT(3)**2) + PT(3)*(
 
635
     $ -PT(0)**2 + PT(1)**2 + PT(2)**2 + PT(3)**2) + PT(0)*DISCR)/(TWO
 
636
     $ *(PT(0)**2 - PT(3)**2))
 
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)
 
641
      NEWP(1,P2) = P(1,P2)
 
642
      NEWP(2,P2) = P(2,P2)
 
643
      DO J=1,2
 
644
        REF=ZERO
 
645
        DO I=NINITIAL+1,NEXTERNAL
 
646
          REF = REF + P(J,I)
 
647
        ENDDO
 
648
        REF = REF - P(J,P2)
 
649
        NEWP(J,P1) = REF
 
650
      ENDDO
 
651
 
 
652
      IF (.NOT.MG5_1_MP_IS_CLOSE(P,NEWP,WARNED)) THEN
 
653
        ERRCODE=999
 
654
        GOTO 100
 
655
      ENDIF
 
656
 
 
657
      DO J=1,NEXTERNAL
 
658
        DO I=0,3
 
659
          P(I,J)=NEWP(I,J)
 
660
        ENDDO
 
661
      ENDDO
 
662
 
 
663
 100  CONTINUE
 
664
 
 
665
      END
 
666
 
 
667
C     *****************************************************************
 
668
C     Beginning of the routine for restoring precision a la PSMC
 
669
C     *****************************************************************
 
670
 
 
671
      SUBROUTINE MG5_1_MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE
 
672
     $ ,WARNED)
 
673
      IMPLICIT NONE
 
674
C     
 
675
C     CONSTANTS 
 
676
C     
 
677
      INTEGER    NEXTERNAL
 
678
      PARAMETER (NEXTERNAL=4)
 
679
      INTEGER    NINITIAL
 
680
      PARAMETER (NINITIAL=2)
 
681
      REAL*16 ZERO
 
682
      PARAMETER (ZERO=0.0E+00_16)
 
683
      REAL*16 MP__ZERO
 
684
      PARAMETER (MP__ZERO=ZERO)
 
685
      REAL*16 ONE
 
686
      PARAMETER (ONE=1.0E+00_16)
 
687
      REAL*16 TWO
 
688
      PARAMETER (TWO=2.0E+00_16)
 
689
      REAL*16 CONSISTENCY_THRES
 
690
      PARAMETER (CONSISTENCY_THRES=1.0E-25_16)
 
691
 
 
692
      INTEGER    NAPPROXZEROS
 
693
      PARAMETER (NAPPROXZEROS=3)
 
694
 
 
695
C     
 
696
C     ARGUMENTS 
 
697
C     
 
698
      REAL*16 P(0:3,NEXTERNAL)
 
699
      INTEGER ERRCODE,ERROR,WARNED
 
700
C     
 
701
C     FUNCTIONS
 
702
C     
 
703
      LOGICAL MG5_1_MP_IS_CLOSE
 
704
C     
 
705
C     LOCAL VARIABLES 
 
706
C     
 
707
      INTEGER I,J, P1, P2
 
708
      REAL*16 NEWP(0:3,NEXTERNAL), PBUFF(0:3)
 
709
      REAL*16 BUFF, BUFF2, XSCALE, APPROX_ZEROS(NAPPROXZEROS)
 
710
      REAL*16 MASSES(NEXTERNAL)
 
711
C     
 
712
C     GLOBAL VARIABLES
 
713
C     
 
714
      INCLUDE 'mp_coupl.inc'
 
715
 
 
716
C     ----------
 
717
C     BEGIN CODE
 
718
C     ----------
 
719
 
 
720
      MASSES(1)=MP__ZERO
 
721
      MASSES(2)=MP__ZERO
 
722
      MASSES(3)=MP__ZERO
 
723
      MASSES(4)=MP__ZERO
 
724
 
 
725
      ERRCODE = 0
 
726
      XSCALE = ONE
 
727
 
 
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
 
732
 
 
733
C     Start by copying the momenta
 
734
      DO I=1,NEXTERNAL
 
735
        DO J=0,3
 
736
          NEWP(J,I)=P(J,I)
 
737
        ENDDO
 
738
      ENDDO
 
739
 
 
740
C     First make sur that the space like momentum is exactly conserved
 
741
      DO J=0,3
 
742
        PBUFF(J)=ZERO
 
743
      ENDDO
 
744
      DO I=1,NINITIAL
 
745
        DO J=1,3
 
746
          PBUFF(J)=PBUFF(J)+NEWP(J,I)
 
747
        ENDDO
 
748
      ENDDO
 
749
      DO I=NINITIAL+1,NEXTERNAL-1
 
750
        DO J=1,3
 
751
          PBUFF(J)=PBUFF(J)-NEWP(J,I)
 
752
        ENDDO
 
753
      ENDDO
 
754
      DO J=1,3
 
755
        NEWP(J,NEXTERNAL)=PBUFF(J)
 
756
      ENDDO
 
757
 
 
758
C     Now find the 'x' rescaling factor
 
759
      DO I=1,NAPPROXZEROS
 
760
        CALL MG5_1_FINDX(NEWP,APPROX_ZEROS(I),XSCALE,ERROR)
 
761
        IF(ERROR.EQ.0) THEN
 
762
          GOTO 1001
 
763
        ELSE
 
764
          ERRCODE=ERRCODE+(10**(I-1))*ERROR
 
765
        ENDIF
 
766
      ENDDO
 
767
      IF (WARNED.LT.20) THEN
 
768
        WRITE(*,*) 'WARNING:: Could not find the proper rescalin'
 
769
     $   //'g factor x. Restoring precision ala PSMC will therefor'
 
770
     $   //'e not be used.'
 
771
        WARNED=WARNED+1
 
772
      ENDIF
 
773
      IF (ERRCODE.LT.1000) THEN
 
774
        ERRCODE=ERRCODE+1000
 
775
      ENDIF
 
776
      GOTO 1000
 
777
 1001 CONTINUE
 
778
      ERRCODE = 0
 
779
 
 
780
C     Apply the rescaling
 
781
      DO I=1,NEXTERNAL
 
782
        DO J=1,3
 
783
          NEWP(J,I)=NEWP(J,I)*XSCALE
 
784
        ENDDO
 
785
      ENDDO
 
786
 
 
787
C     Now restore exact onshellness of the particles.
 
788
      DO I=1,NEXTERNAL
 
789
        BUFF=MASSES(I)**2
 
790
        DO J=1,3
 
791
          BUFF=BUFF+NEWP(J,I)**2
 
792
        ENDDO
 
793
        NEWP(0,I)=SQRT(BUFF)
 
794
      ENDDO
 
795
 
 
796
C     Consistency check
 
797
      BUFF=ZERO
 
798
      BUFF2=ZERO
 
799
      DO I=1,NINITIAL
 
800
        BUFF=BUFF-NEWP(0,I)
 
801
        BUFF2=BUFF2+NEWP(0,I)
 
802
      ENDDO
 
803
      DO I=NINITIAL+1,NEXTERNAL
 
804
        BUFF=BUFF+NEWP(0,I)
 
805
        BUFF2=BUFF2+NEWP(0,I)
 
806
      ENDDO
 
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 PSM'
 
810
     $     //'C precision restoring algorithm failed. The result wil'
 
811
     $     //'l therefore not be used.'
 
812
          WARNED=WARNED+1
 
813
        ENDIF
 
814
        ERRCODE = 1000
 
815
        GOTO 1000
 
816
      ENDIF
 
817
 
 
818
      IF (.NOT.MG5_1_MP_IS_CLOSE(P,NEWP,WARNED)) THEN
 
819
        ERRCODE=999
 
820
        GOTO 1000
 
821
      ENDIF
 
822
 
 
823
      DO J=1,NEXTERNAL
 
824
        DO I=0,3
 
825
          P(I,J)=NEWP(I,J)
 
826
        ENDDO
 
827
      ENDDO
 
828
 
 
829
 1000 CONTINUE
 
830
 
 
831
      END
 
832
 
 
833
 
 
834
      SUBROUTINE MG5_1_FINDX(P,SEED,XSCALE,ERROR)
 
835
      IMPLICIT NONE
 
836
C     
 
837
C     CONSTANTS 
 
838
C     
 
839
      INTEGER    NEXTERNAL
 
840
      PARAMETER (NEXTERNAL=4)
 
841
      INTEGER    NINITIAL
 
842
      PARAMETER (NINITIAL=2)
 
843
      REAL*16 ZERO
 
844
      PARAMETER (ZERO=0.0E+00_16)
 
845
      REAL*16 MP__ZERO
 
846
      PARAMETER (MP__ZERO=ZERO)
 
847
      REAL*16 ONE
 
848
      PARAMETER (ONE=1.0E+00_16)
 
849
      REAL*16 TWO
 
850
      PARAMETER (TWO=2.0E+00_16)
 
851
      INTEGER MAXITERATIONS
 
852
      PARAMETER (MAXITERATIONS=8)
 
853
      REAL*16 CONVERGED
 
854
      PARAMETER (CONVERGED=1.0E-26_16)
 
855
C     
 
856
C     ARGUMENTS 
 
857
C     
 
858
      REAL*16 P(0:3,NEXTERNAL),SEED,XSCALE
 
859
      INTEGER ERROR
 
860
C     
 
861
C     LOCAL VARIABLES 
 
862
C     
 
863
      INTEGER I,J,ERR
 
864
      REAL*16 PVECSQ(NEXTERNAL)
 
865
      REAL*16 XN, XNP1,FVAL,DVAL
 
866
 
 
867
C     ----------
 
868
C     BEGIN CODE
 
869
C     ----------
 
870
 
 
871
      ERROR = 0
 
872
      XSCALE = SEED
 
873
      XN = SEED
 
874
      XNP1 = SEED
 
875
 
 
876
      DO I=1,NEXTERNAL
 
877
        PVECSQ(I)=P(1,I)**2+P(2,I)**2+P(3,I)**2
 
878
      ENDDO
 
879
 
 
880
      DO I=1,MAXITERATIONS
 
881
        CALL MG5_1_FUNCT(PVECSQ(1),XN,.FALSE.,ERR, FVAL)
 
882
        IF (ERR.NE.0) THEN
 
883
          ERROR=ERR
 
884
          GOTO 710
 
885
        ENDIF
 
886
        CALL MG5_1_FUNCT(PVECSQ(1),XN,.TRUE.,ERR, DVAL)
 
887
        IF (ERR.NE.0) THEN
 
888
          ERROR=ERR
 
889
          GOTO 710
 
890
        ENDIF
 
891
        XNP1=XN-(FVAL/DVAL)
 
892
        IF((ABS(((XNP1-XN)*TWO)/(XNP1+XN))).LT.CONVERGED) THEN
 
893
          XN=XNP1
 
894
          GOTO 700
 
895
        ENDIF
 
896
        XN=XNP1
 
897
      ENDDO
 
898
      ERROR=9
 
899
      GOTO 710
 
900
 
 
901
 700  CONTINUE
 
902
C     For good measure, we iterate one last time
 
903
      CALL MG5_1_FUNCT(PVECSQ(1),XN,.FALSE.,ERR, FVAL)
 
904
      IF (ERR.NE.0) THEN
 
905
        ERROR=ERR
 
906
        GOTO 710
 
907
      ENDIF
 
908
      CALL MG5_1_FUNCT(PVECSQ(1),XN,.TRUE.,ERR, DVAL)
 
909
      IF (ERR.NE.0) THEN
 
910
        ERROR=ERR
 
911
        GOTO 710
 
912
      ENDIF
 
913
 
 
914
      XSCALE=XN-(FVAL/DVAL)
 
915
 
 
916
 710  CONTINUE
 
917
 
 
918
      END
 
919
 
 
920
      SUBROUTINE MG5_1_FUNCT(PVECSQ,X,DERIVATIVE,ERROR,RES)
 
921
      IMPLICIT NONE
 
922
C     
 
923
C     CONSTANTS 
 
924
C     
 
925
      INTEGER    NEXTERNAL
 
926
      PARAMETER (NEXTERNAL=4)
 
927
      INTEGER    NINITIAL
 
928
      PARAMETER (NINITIAL=2)
 
929
      REAL*16 ZERO
 
930
      PARAMETER (ZERO=0.0E+00_16)
 
931
      REAL*16 MP__ZERO
 
932
      PARAMETER (MP__ZERO=ZERO)
 
933
      REAL*16 ONE
 
934
      PARAMETER (ONE=1.0E+00_16)
 
935
      REAL*16 TWO
 
936
      PARAMETER (TWO=2.0E+00_16)
 
937
C     
 
938
C     ARGUMENTS 
 
939
C     
 
940
      REAL*16 PVECSQ(NEXTERNAL),X,RES
 
941
      INTEGER ERROR
 
942
      LOGICAL DERIVATIVE
 
943
C     
 
944
C     LOCAL VARIABLES 
 
945
C     
 
946
      INTEGER I,J
 
947
      REAL*16 BUFF,FACTOR
 
948
      REAL*16 MASSES(NEXTERNAL)
 
949
C     
 
950
C     GLOBAL VARIABLES
 
951
C     
 
952
      INCLUDE 'mp_coupl.inc'
 
953
 
 
954
C     ----------
 
955
C     BEGIN CODE
 
956
C     ----------
 
957
 
 
958
      MASSES(1)=MP__ZERO
 
959
      MASSES(2)=MP__ZERO
 
960
      MASSES(3)=MP__ZERO
 
961
      MASSES(4)=MP__ZERO
 
962
 
 
963
      ERROR=0
 
964
      RES=ZERO
 
965
      BUFF=ZERO
 
966
 
 
967
      DO I=1,NEXTERNAL
 
968
        IF (I.LE.NINITIAL) THEN
 
969
          FACTOR=-ONE
 
970
        ELSE
 
971
          FACTOR=ONE
 
972
        ENDIF
 
973
        BUFF=MASSES(I)**2+PVECSQ(I)*X**2
 
974
        IF (BUFF.LT.ZERO) THEN
 
975
          RES=ZERO
 
976
          ERROR = 1
 
977
          GOTO 800
 
978
        ENDIF
 
979
        IF (DERIVATIVE) THEN
 
980
          RES=RES + FACTOR*((X*PVECSQ(I))/SQRT(BUFF))
 
981
        ELSE
 
982
          RES=RES + FACTOR*SQRT(BUFF)
 
983
        ENDIF
 
984
      ENDDO
 
985
 
 
986
 800  CONTINUE
 
987
 
 
988
      END
 
989