← Back to branch summary

~maddevelopers/mg5amcnlo/2.5.3

« 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: 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/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