← Back to branch summary

~maddevelopers/mg5amcnlo/2.5.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%check_sa.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/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%check_sa.f
1
 
      PROGRAM DRIVER
2
 
C     *****************************************************************
3
 
C     ********
4
 
C     THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT.
5
 
C     IT USES A SIMPLE PHASE SPACE GENERATOR
6
 
C     *****************************************************************
7
 
C     ********
8
 
      IMPLICIT NONE
9
 
C     
10
 
C     CONSTANTS  
11
 
C     
12
 
      REAL*8 ZERO
13
 
      PARAMETER (ZERO=0D0)
14
 
 
15
 
      LOGICAL READPS
16
 
      PARAMETER (READPS = .FALSE.)
17
 
 
18
 
      INTEGER NPSPOINTS
19
 
      PARAMETER (NPSPOINTS = 4)
20
 
 
21
 
C     integer nexternal and number particles (incoming+outgoing) in
22
 
C      the me 
23
 
      INTEGER NEXTERNAL, NINCOMING
24
 
      PARAMETER (NEXTERNAL=3,NINCOMING=2)
25
 
 
26
 
      CHARACTER(512) MADLOOPRESOURCEPATH
27
 
 
28
 
C     
29
 
C     INCLUDE FILES
30
 
C     
31
 
C     the include file with the values of the parameters and masses   
32
 
C        
33
 
      INCLUDE 'coupl.inc'
34
 
C     particle masses
35
 
      REAL*8 PMASS(NEXTERNAL)
36
 
C     integer    n_max_cg
37
 
      INCLUDE 'ngraphs.inc'
38
 
      INCLUDE 'nsqso_born.inc'
39
 
      INCLUDE 'nsquaredSO.inc'
40
 
 
41
 
C     
42
 
C     LOCAL
43
 
C     
44
 
      INTEGER I,J,K
45
 
C     four momenta. Energy is the zeroth component.
46
 
      REAL*8 P(0:3,NEXTERNAL)
47
 
      INTEGER MATELEM_ARRAY_DIM
48
 
      REAL*8 , ALLOCATABLE :: MATELEM(:,:)
49
 
      REAL*8 SQRTS,AO2PI,TOTMASS
50
 
C     sqrt(s)= center of mass energy 
51
 
      REAL*8 PIN(0:3), POUT(0:3)
52
 
      CHARACTER*120 BUFF(NEXTERNAL)
53
 
      INTEGER RETURNCODE, UNITS, TENS, HUNDREDS
54
 
      INTEGER NSQUAREDSO_LOOP
55
 
      REAL*8 , ALLOCATABLE :: PREC_FOUND(:)
56
 
 
57
 
C     
58
 
C     GLOBAL VARIABLES
59
 
C     
60
 
C     This is from ML code for the list of split orders selected by
61
 
C     the process definition
62
 
C     
63
 
      INTEGER NLOOPCHOSEN
64
 
      CHARACTER*20 CHOSEN_LOOP_SO_INDICES(NSQUAREDSO)
65
 
      LOGICAL CHOSEN_LOOP_SO_CONFIGS(NSQUAREDSO)
66
 
      COMMON/CHOSEN_LOOP_SQSO/CHOSEN_LOOP_SO_CONFIGS
67
 
      INTEGER NBORNCHOSEN
68
 
      CHARACTER*20 CHOSEN_BORN_SO_INDICES(NSQSO_BORN)
69
 
      LOGICAL CHOSEN_BORN_SO_CONFIGS(NSQSO_BORN)
70
 
      COMMON/CHOSEN_BORN_SQSO/CHOSEN_BORN_SO_CONFIGS
71
 
 
72
 
C     
73
 
C     SAVED VARIABLES
74
 
C     
75
 
      LOGICAL INIT
76
 
      DATA INIT/.TRUE./
77
 
      COMMON/INITCHECKSA/INIT
78
 
C     
79
 
C     EXTERNAL
80
 
C     
81
 
      REAL*8 DOT
82
 
      EXTERNAL DOT
83
 
 
84
 
C     
85
 
C     BEGIN CODE
86
 
C     
87
 
C     
88
 
 
89
 
      IF (INIT) THEN
90
 
        INIT=.FALSE.
91
 
        CALL GET_ANSWER_DIMENSION(MATELEM_ARRAY_DIM)
92
 
        ALLOCATE(MATELEM(0:3,0:MATELEM_ARRAY_DIM))
93
 
        CALL GET_NSQSO_LOOP(NSQUAREDSO_LOOP)
94
 
        ALLOCATE(PREC_FOUND(0:NSQUAREDSO_LOOP))
95
 
 
96
 
C       INITIALIZATION CALLS
97
 
C       
98
 
C       Call to initialize the values of the couplings, masses and
99
 
C        widths 
100
 
C       used in the evaluation of the matrix element. The primary
101
 
C        parameters of the
102
 
C       models are read from Cards/param_card.dat. The secondary
103
 
C        parameters are calculated
104
 
C       in Source/MODEL/couplings.f. The values are stored in common
105
 
C        blocks that are listed
106
 
C       in coupl.inc .
107
 
C       first call to setup the paramaters
108
 
        CALL SETPARA('param_card.dat')
109
 
C       set up masses
110
 
        INCLUDE 'pmass.inc'
111
 
 
112
 
      ENDIF
113
 
 
114
 
 
115
 
C     Start by initializing what is the squared split orders indices
116
 
C      chosen
117
 
      NLOOPCHOSEN=0
118
 
      DO I=1,NSQUAREDSO
119
 
        IF (CHOSEN_LOOP_SO_CONFIGS(I)) THEN
120
 
          NLOOPCHOSEN=NLOOPCHOSEN+1
121
 
          WRITE(CHOSEN_LOOP_SO_INDICES(NLOOPCHOSEN),'(I3,A2)') I,'L)'
122
 
        ENDIF
123
 
      ENDDO
124
 
      NBORNCHOSEN=0
125
 
      DO I=1,NSQSO_BORN
126
 
        IF (CHOSEN_BORN_SO_CONFIGS(I)) THEN
127
 
          NBORNCHOSEN=NBORNCHOSEN+1
128
 
          WRITE(CHOSEN_BORN_SO_INDICES(NBORNCHOSEN),'(I3,A2)') I,'B)'
129
 
        ENDIF
130
 
      ENDDO
131
 
 
132
 
      AO2PI=G**2/(8.D0*(3.14159265358979323846D0**2))
133
 
 
134
 
      WRITE(*,*) 'AO2PI=',AO2PI
135
 
C     Now use a simple multipurpose PS generator (RAMBO) just to get a 
136
 
C     RANDOM set of four momenta of given masses pmass(i) to be used
137
 
C      to evaluate 
138
 
C     the madgraph matrix-element.       
139
 
C     Alternatevely, here the user can call or set the four momenta at
140
 
C      his will, see below.
141
 
C     
142
 
      IF(NINCOMING.EQ.1) THEN
143
 
        SQRTS=PMASS(1)
144
 
      ELSE
145
 
        TOTMASS = 0.0D0
146
 
        DO I=1,NEXTERNAL
147
 
          TOTMASS = TOTMASS + PMASS(I)
148
 
        ENDDO
149
 
C       CMS energy in GEV
150
 
        SQRTS=MAX(1000D0,2.0D0*TOTMASS)
151
 
      ENDIF
152
 
 
153
 
      CALL PRINTOUT()
154
 
 
155
 
 
156
 
 
157
 
      DO K=1,NPSPOINTS
158
 
 
159
 
        IF(READPS) THEN
160
 
          OPEN(967, FILE='PS.input', ERR=976, STATUS='OLD',
161
 
     $      ACTION='READ')
162
 
          DO I=1,NEXTERNAL
163
 
            READ(967,*,END=978) P(0,I),P(1,I),P(2,I),P(3,I)
164
 
          ENDDO
165
 
          GOTO 978
166
 
 976      CONTINUE
167
 
          STOP 'Could not read the PS.input phase-space point.'
168
 
 978      CONTINUE
169
 
          CLOSE(967)
170
 
        ELSE
171
 
          IF ((NINCOMING.EQ.2).AND.((NEXTERNAL - NINCOMING .EQ.1)))
172
 
     $      THEN
173
 
            IF (PMASS(3).EQ.0.0D0) THEN
174
 
              STOP 'Cannot generate 2>1 kin. config. with m3=0.0d0'
175
 
            ELSE
176
 
C             deal with the case of only one particle in the final
177
 
C              state
178
 
              P(0,1) = PMASS(3)/2D0
179
 
              P(1,1) = 0D0
180
 
              P(2,1) = 0D0
181
 
              P(3,1) = PMASS(3)/2D0
182
 
              IF (PMASS(1).GT.0D0) THEN
183
 
                P(3,1) = DSQRT(PMASS(3)**2/4D0 - PMASS(1)**2)
184
 
              ENDIF
185
 
              P(0,2) = PMASS(3)/2D0
186
 
              P(1,2) = 0D0
187
 
              P(2,2) = 0D0
188
 
              P(3,2) = -PMASS(3)/2D0
189
 
              IF (PMASS(2) > 0D0) THEN
190
 
                P(3,2) = -DSQRT(PMASS(3)**2/4D0 - PMASS(1)**2)
191
 
              ENDIF
192
 
              P(0,3) = PMASS(3)
193
 
              P(1,3) = 0D0
194
 
              P(2,3) = 0D0
195
 
              P(3,3) = 0D0
196
 
            ENDIF
197
 
          ELSE
198
 
            CALL GET_MOMENTA(SQRTS,PMASS,P)
199
 
          ENDIF
200
 
        ENDIF
201
 
 
202
 
        DO I=0,3
203
 
          PIN(I)=0.0D0
204
 
          DO J=1,NINCOMING
205
 
            PIN(I)=PIN(I)+P(I,J)
206
 
          ENDDO
207
 
        ENDDO
208
 
 
209
 
C       In standalone mode, always use sqrt_s as the renormalization
210
 
C        scale.
211
 
        SQRTS=DSQRT(DABS(DOT(PIN(0),PIN(0))))
212
 
        MU_R=SQRTS
213
 
 
214
 
C       Update the couplings with the new MU_R
215
 
        CALL UPDATE_AS_PARAM()
216
 
 
217
 
C       Optionally the user can set where to find the
218
 
C        MadLoop5_resources folder.
219
 
C       Otherwise it will look for it automatically and find it if it
220
 
C        has not
221
 
C       been moved
222
 
C       MadLoopResourcePath = '<MadLoop5_resources_path>'
223
 
C       CALL SETMADLOOPPATH(MadLoopResourcePath)
224
 
C       To force the stabiliy check to also be performed in the
225
 
C        initialization phase
226
 
C       CALL FORCE_STABILITY_CHECK(.TRUE.)
227
 
C       To chose a particular tartget split order, SOTARGET is an
228
 
C        integer labeling
229
 
C       the possible squared order couplings contributions (only in
230
 
C        optimized mode)
231
 
C       CALL SET_COUPLINGORDERS_TARGET(SOTARGET)
232
 
 
233
 
 
234
 
C       
235
 
C       Now we can call the matrix element
236
 
C       
237
 
        CALL SLOOPMATRIX_THRES(P,MATELEM,-1.0D0,PREC_FOUND,RETURNCODE)
238
 
 
239
 
C       
240
 
C       write the information on the four momenta 
241
 
C       
242
 
        IF (K.EQ.NPSPOINTS) THEN
243
 
          WRITE (*,*)
244
 
          WRITE (*,*) ' Phase space point:'
245
 
          WRITE (*,*)
246
 
          WRITE (*,*) '---------------------------------'
247
 
          WRITE (*,*)  'n  E  px py pz m'
248
 
          DO I=1,NEXTERNAL
249
 
            WRITE (*,'(i2,1x,5e15.7)') I, P(0,I),P(1,I),P(2,I),P(3,I)
250
 
     $       ,DSQRT(DABS(DOT(P(0,I),P(0,I))))
251
 
          ENDDO
252
 
          WRITE (*,*) '---------------------------------'
253
 
          WRITE (*,*) 'Detailed result for each coupling orders'
254
 
     $     //' combination.'
255
 
          WRITE(*,*) 'All Born contributions are of split orders'
256
 
     $     //' (QCD=0)'
257
 
          WRITE (*,*) '---------------------------------'
258
 
          WRITE(*,*) 'All loop contributions are of split orders'
259
 
     $     //' (QCD=2)'
260
 
          WRITE (*,*) '---------------------------------'
261
 
          UNITS=MOD(RETURNCODE,10)
262
 
          TENS=(MOD(RETURNCODE,100)-UNITS)/10
263
 
          HUNDREDS=(RETURNCODE-TENS*10-UNITS)/100
264
 
          IF (HUNDREDS.EQ.1) THEN
265
 
            IF (TENS.EQ.3.OR.TENS.EQ.4) THEN
266
 
              WRITE(*,*) 'Unknown numerical stability because MadLoop'
267
 
     $         //' is in the initialization stage.'
268
 
            ELSE
269
 
              WRITE(*,*) 'Unknown numerical stability, check CTModeRun'
270
 
     $         //' value in MadLoopParams.dat.'
271
 
            ENDIF
272
 
          ELSEIF (HUNDREDS.EQ.2) THEN
273
 
            WRITE(*,*) 'Stable kinematic configuration (SPS).'
274
 
          ELSEIF (HUNDREDS.EQ.3) THEN
275
 
            WRITE(*,*) 'Unstable kinematic configuration (UPS).'
276
 
            WRITE(*,*) 'Quadruple precision rescue successful.'
277
 
          ELSEIF (HUNDREDS.EQ.4) THEN
278
 
            WRITE(*,*) 'Exceptional kinematic configuration (EPS).'
279
 
            WRITE(*,*) 'Both double an quadruple precision'
280
 
     $       //' computations, are unstable.'
281
 
          ENDIF
282
 
          IF (TENS.EQ.2.OR.TENS.EQ.4) THEN
283
 
            WRITE(*,*) 'Quadruple precision computation used.'
284
 
          ENDIF
285
 
          IF (HUNDREDS.NE.1) THEN
286
 
            IF (PREC_FOUND(0).GT.0.0D0) THEN
287
 
              WRITE(*,'(1x,a23,1x,1e10.2)') 'Relative accuracy     ='
288
 
     $         ,PREC_FOUND(0)
289
 
            ELSEIF (PREC_FOUND(0).EQ.0.0D0) THEN
290
 
              WRITE(*,'(1x,a23,1x,1e10.2,1x,a30)') 'Relative accuracy '
291
 
     $         //'    =',PREC_FOUND(0),'(i.e. beyond double precision)'
292
 
            ELSE
293
 
              WRITE(*,*) 'Estimated accuracy could not be computed for'
294
 
     $         //' an unknown reason.'
295
 
            ENDIF
296
 
          ENDIF
297
 
          WRITE (*,'(1x,a23,3x,i3)') 'MadLoop return code   ='
298
 
     $     ,RETURNCODE
299
 
          WRITE (*,*) '---------------------------------'
300
 
          IF (NBORNCHOSEN.EQ.0) THEN
301
 
            WRITE (*,*) 'No Born contribution satisfied the squared'
302
 
     $       //' order constraints.'
303
 
          ELSE IF (NBORNCHOSEN.NE.NSQSO_BORN) THEN
304
 
            WRITE (*,*) 'Selected squared coupling orders combination'
305
 
     $       //' for the Born summed result below:'
306
 
            WRITE (*,*) (CHOSEN_BORN_SO_INDICES(I),I=1,NBORNCHOSEN)
307
 
          ENDIF
308
 
          IF (NLOOPCHOSEN.NE.NSQUAREDSO) THEN
309
 
            WRITE (*,*) 'Selected squared coupling orders combination'
310
 
     $       //' for the loop summed result below:'
311
 
            WRITE (*,*) (CHOSEN_LOOP_SO_INDICES(I),I=1,NLOOPCHOSEN)
312
 
          ENDIF
313
 
          WRITE (*,*) '---------------------------------'
314
 
          WRITE (*,*) 'Matrix element born   = ', MATELEM(0,0), 
315
 
     $     ' GeV^',-(2*NEXTERNAL-8)
316
 
          WRITE (*,*) 'Matrix element finite = ', MATELEM(1,0), 
317
 
     $     ' GeV^',-(2*NEXTERNAL-8)
318
 
          WRITE (*,*) 'Matrix element 1eps   = ', MATELEM(2,0), 
319
 
     $     ' GeV^',-(2*NEXTERNAL-8)
320
 
          WRITE (*,*) 'Matrix element 2eps   = ', MATELEM(3,0), 
321
 
     $     ' GeV^',-(2*NEXTERNAL-8)
322
 
          WRITE (*,*) '---------------------------------'
323
 
          IF (MATELEM(0,0).NE.0.0D0) THEN
324
 
            WRITE (*,*) 'finite / (born*ao2pi) = ', MATELEM(1,0)
325
 
     $       /MATELEM(0,0)/AO2PI
326
 
            WRITE (*,*) '1eps   / (born*ao2pi) = ', MATELEM(2,0)
327
 
     $       /MATELEM(0,0)/AO2PI
328
 
            WRITE (*,*) '2eps   / (born*ao2pi) = ', MATELEM(3,0)
329
 
     $       /MATELEM(0,0)/AO2PI
330
 
          ELSE
331
 
            WRITE (*,*) 'finite / ao2pi      = ', MATELEM(1,0)/AO2PI
332
 
            WRITE (*,*) '1eps   / ao2pi      = ', MATELEM(2,0)/AO2PI
333
 
            WRITE (*,*) '2eps   / ao2pi      = ', MATELEM(3,0)/AO2PI
334
 
          ENDIF
335
 
          WRITE (*,*) '---------------------------------'
336
 
 
337
 
          OPEN(69, FILE='result.dat', ERR=976, ACTION='WRITE')
338
 
          DO I=1,NEXTERNAL
339
 
            WRITE (69,'(a2,1x,5ES30.15E3)') 'PS',P(0,I),P(1,I),P(2,I)
340
 
     $       ,P(3,I)
341
 
          ENDDO
342
 
          WRITE (69,'(a3,1x,i3)') 'EXP',-(2*NEXTERNAL-8)
343
 
          WRITE (69,'(a4,1x,1ES30.15E3)') 'BORN',MATELEM(0,0)
344
 
          IF (MATELEM(0,0).NE.0.0D0) THEN
345
 
            WRITE (69,'(a3,1x,1ES30.15E3)') 'FIN',MATELEM(1,0)
346
 
     $       /MATELEM(0,0)/AO2PI
347
 
            WRITE (69,'(a4,1x,1ES30.15E3)') '1EPS',MATELEM(2,0)
348
 
     $       /MATELEM(0,0)/AO2PI
349
 
            WRITE (69,'(a4,1x,1ES30.15E3)') '2EPS',MATELEM(3,0)
350
 
     $       /MATELEM(0,0)/AO2PI
351
 
          ELSE
352
 
            WRITE (69,'(a3,1x,1ES30.15E3)') 'FIN',MATELEM(1,0)/AO2PI
353
 
            WRITE (69,'(a4,1x,1ES30.15E3)') '1EPS',MATELEM(2,0)/AO2PI
354
 
            WRITE (69,'(a4,1x,1ES30.15E3)') '2EPS',MATELEM(3,0)/AO2PI
355
 
          ENDIF
356
 
          WRITE (69,'(a6,1x,1ES30.15E3)') 'ASO2PI',AO2PI
357
 
          WRITE (69,*) 'Export_Format Default'
358
 
          WRITE (69,'(a7,1x,i3)') 'RETCODE',RETURNCODE
359
 
          WRITE (69,'(a3,1x,1e10.4)') 'ACC',PREC_FOUND(0)
360
 
          WRITE (69,*) 'Born_kept',(CHOSEN_BORN_SO_CONFIGS(I),I=1
361
 
     $     ,NSQSO_BORN)
362
 
          WRITE (69,*) 'Loop_kept',(CHOSEN_LOOP_SO_CONFIGS(I),I=1
363
 
     $     ,NSQUAREDSO)
364
 
          WRITE (69,*) 'Born_SO_Results 0'
365
 
          WRITE (69,*) 'SO_Born BORN ',MATELEM(0,1)
366
 
          WRITE (69,*) 'Split_Orders_Names QCD'
367
 
          WRITE (69,*) 'Loop_SO_Results 2'
368
 
          WRITE (69,*) 'SO_Loop ACC  ',PREC_FOUND(1)
369
 
          WRITE (69,*) 'SO_Loop FIN  ',MATELEM(1,1)
370
 
          WRITE (69,*) 'SO_Loop 1EPS ',MATELEM(2,1)
371
 
          WRITE (69,*) 'SO_Loop 2EPS ',MATELEM(3,1)
372
 
          CLOSE(69)
373
 
        ELSE
374
 
          WRITE (*,*) 'PS Point #',K,' done.'
375
 
        ENDIF
376
 
      ENDDO
377
 
 
378
 
C     C
379
 
C     C      Copy down here (or read in) the four momenta as a string. 
380
 
C     C      
381
 
C     C
382
 
C     buff(1)=" 1   0.5630480E+04  0.0000000E+00  0.0000000E+00 
383
 
C      0.5630480E+04"
384
 
C     buff(2)=" 2   0.5630480E+04  0.0000000E+00  0.0000000E+00
385
 
C      -0.5630480E+04"
386
 
C     buff(3)=" 3   0.5466073E+04  0.4443190E+03  0.2446331E+04
387
 
C      -0.4864732E+04"
388
 
C     buff(4)=" 4   0.8785819E+03 -0.2533886E+03  0.2741971E+03 
389
 
C      0.7759741E+03"
390
 
C     buff(5)=" 5   0.4916306E+04 -0.1909305E+03 -0.2720528E+04 
391
 
C      0.4088757E+04"
392
 
C     C
393
 
C     C      Here the k,E,px,py,pz are read from the string into the
394
 
C      momenta array.
395
 
C     C      k=1,2          : incoming
396
 
C     C      k=3,nexternal  : outgoing
397
 
C     C
398
 
C     do i=1,nexternal
399
 
C     read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i)
400
 
C     enddo
401
 
C     
402
 
C     C print the momenta out
403
 
C     
404
 
C     do i=1,nexternal
405
 
C     write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i), 
406
 
C     &dsqrt(dabs(DOT(p(0,i),p(0,i))))
407
 
C     enddo
408
 
C     
409
 
C     CALL SLOOPMATRIX(P,MATELEM)
410
 
C     
411
 
C     write (*,*) "-------------------------------------------------"
412
 
C     write (*,*) "Matrix element = ", MATELEM(1), "
413
 
C      GeV^",-(2*nexternal-8)      
414
 
C     write (*,*) "-------------------------------------------------"
415
 
 
416
 
      DEALLOCATE(MATELEM)
417
 
      DEALLOCATE(PREC_FOUND)
418
 
 
419
 
      END
420
 
 
421
 
 
422
 
 
423
 
 
424
 
      DOUBLE PRECISION FUNCTION DOT(P1,P2)
425
 
C     *************************************************************
426
 
C     4-Vector Dot product
427
 
C     *************************************************************
428
 
      IMPLICIT NONE
429
 
      DOUBLE PRECISION P1(0:3),P2(0:3)
430
 
      DOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
431
 
      END
432
 
 
433
 
 
434
 
      SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P)
435
 
C     auxiliary function to change convention between madgraph and
436
 
C      rambo
437
 
C     four momenta.         
438
 
      IMPLICIT NONE
439
 
      INTEGER NEXTERNAL, NINCOMING
440
 
      PARAMETER (NEXTERNAL=3,NINCOMING=2)
441
 
C     ARGUMENTS
442
 
      REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT
443
 
C     LOCAL
444
 
      INTEGER I
445
 
      REAL*8 ETOT2,MOM,M1,M2,E1,E2
446
 
 
447
 
      ETOT2=ENERGY**2
448
 
      M1=PMASS(1)
449
 
      M2=PMASS(2)
450
 
      MOM=(ETOT2**2 - 2*ETOT2*M1**2 + M1**4 - 2*ETOT2*M2**2 - 2*M1**2
451
 
     $ *M2**2 + M2**4)/(4.*ETOT2)
452
 
      MOM=DSQRT(MOM)
453
 
      E1=DSQRT(MOM**2+M1**2)
454
 
      E2=DSQRT(MOM**2+M2**2)
455
 
C     write (*,*) e1+e2,mom
456
 
 
457
 
      IF(NINCOMING.EQ.2) THEN
458
 
 
459
 
        P(0,1)=E1
460
 
        P(1,1)=0D0
461
 
        P(2,1)=0D0
462
 
        P(3,1)=MOM
463
 
 
464
 
        P(0,2)=E2
465
 
        P(1,2)=0D0
466
 
        P(2,2)=0D0
467
 
        P(3,2)=-MOM
468
 
 
469
 
        CALL RAMBO(NEXTERNAL-2,ENERGY,PMASS(3),PRAMBO,WGT)
470
 
        DO I=3, NEXTERNAL
471
 
          P(0,I)=PRAMBO(4,I-2)
472
 
          P(1,I)=PRAMBO(1,I-2)
473
 
          P(2,I)=PRAMBO(2,I-2)
474
 
          P(3,I)=PRAMBO(3,I-2)
475
 
        ENDDO
476
 
 
477
 
      ELSEIF(NINCOMING.EQ.1) THEN
478
 
 
479
 
        P(0,1)=ENERGY
480
 
        P(1,1)=0D0
481
 
        P(2,1)=0D0
482
 
        P(3,1)=0D0
483
 
 
484
 
        CALL RAMBO(NEXTERNAL-1,ENERGY,PMASS(2),PRAMBO,WGT)
485
 
        DO I=2, NEXTERNAL
486
 
          P(0,I)=PRAMBO(4,I-1)
487
 
          P(1,I)=PRAMBO(1,I-1)
488
 
          P(2,I)=PRAMBO(2,I-1)
489
 
          P(3,I)=PRAMBO(3,I-1)
490
 
        ENDDO
491
 
      ENDIF
492
 
 
493
 
      RETURN
494
 
      END
495
 
 
496
 
 
497
 
      SUBROUTINE RAMBO(N,ET,XM,P,WT)
498
 
C     *****************************************************************
499
 
C     *****
500
 
C     RAMBO                                         *
501
 
C     RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)                 
502
 
C      *
503
 
C     *
504
 
C     A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR               
505
 
C      *
506
 
C     AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING                
507
 
C      *
508
 
C     THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS                     
509
 
C      *
510
 
C     -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90)   
511
 
C      *
512
 
C     *
513
 
C     N  = NUMBER OF PARTICLES                                        
514
 
C      *
515
 
C     ET = TOTAL CENTRE-OF-MASS ENERGY                                
516
 
C      *
517
 
C     XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming )                
518
 
C      *
519
 
C     P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) )           
520
 
C      *
521
 
C     WT = WEIGHT OF THE EVENT                                        
522
 
C      *
523
 
C     *****************************************************************
524
 
C     *****
525
 
      IMPLICIT REAL*8(A-H,O-Z)
526
 
      INTEGER NEXTERNAL, NINCOMING
527
 
      PARAMETER (NEXTERNAL=3,NINCOMING=2)
528
 
      DIMENSION XM(NEXTERNAL-NINCOMING),P(4,NEXTERNAL-NINCOMING)
529
 
      DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4)
530
 
     $ ,B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING)
531
 
     $ ,E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5)
532
 
      SAVE ACC,ITMAX,IBEGIN,IWARN
533
 
      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
534
 
C     
535
 
C     INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
536
 
      IF(IBEGIN.NE.0) GOTO 103
537
 
      IBEGIN=1
538
 
      TWOPI=8.*DATAN(1.D0)
539
 
      PO2LOG=LOG(TWOPI/4.)
540
 
      Z(2)=PO2LOG
541
 
      DO 101 K=3,(NEXTERNAL-NINCOMING)
542
 
 101  Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
543
 
      DO 102 K=3,(NEXTERNAL-NINCOMING)
544
 
 102  Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
545
 
C     
546
 
C     CHECK ON THE NUMBER OF PARTICLES
547
 
 103  IF(N.GT.1.AND.N.LT.101) GOTO 104
548
 
      PRINT 1001,N
549
 
      STOP
550
 
C     
551
 
C     CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
552
 
 104  XMT=0.
553
 
      NM=0
554
 
      DO 105 I=1,N
555
 
      IF(XM(I).NE.0.D0) NM=NM+1
556
 
 105  XMT=XMT+ABS(XM(I))
557
 
      IF(XMT.LE.ET) GOTO 201
558
 
      PRINT 1002,XMT,ET
559
 
      STOP
560
 
C     
561
 
C     THE PARAMETER VALUES ARE NOW ACCEPTED
562
 
C     
563
 
C     GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
564
 
 201  DO 202 I=1,N
565
 
      R1=RN(1)
566
 
      C=2.*R1-1.
567
 
      S=SQRT(1.-C*C)
568
 
      F=TWOPI*RN(2)
569
 
      R1=RN(3)
570
 
      R2=RN(4)
571
 
      Q(4,I)=-LOG(R1*R2)
572
 
      Q(3,I)=Q(4,I)*C
573
 
      Q(2,I)=Q(4,I)*S*COS(F)
574
 
 202  Q(1,I)=Q(4,I)*S*SIN(F)
575
 
C     
576
 
C     CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
577
 
      DO 203 I=1,4
578
 
 203  R(I)=0.
579
 
      DO 204 I=1,N
580
 
      DO 204 K=1,4
581
 
 204  R(K)=R(K)+Q(K,I)
582
 
      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
583
 
      DO 205 K=1,3
584
 
 205  B(K)=-R(K)/RMAS
585
 
      G=R(4)/RMAS
586
 
      A=1./(1.+G)
587
 
      X=ET/RMAS
588
 
C     
589
 
C     TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
590
 
      DO 207 I=1,N
591
 
      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
592
 
      DO 206 K=1,3
593
 
 206  P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
594
 
 207  P(4,I)=X*(G*Q(4,I)+BQ)
595
 
C     
596
 
C     CALCULATE WEIGHT AND POSSIBLE WARNINGS
597
 
      WT=PO2LOG
598
 
      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
599
 
      IF(WT.GE.-180.D0) GOTO 208
600
 
      IF(IWARN(1).LE.5) PRINT 1004,WT
601
 
      IWARN(1)=IWARN(1)+1
602
 
 208  IF(WT.LE. 174.D0) GOTO 209
603
 
      IF(IWARN(2).LE.5) PRINT 1005,WT
604
 
      IWARN(2)=IWARN(2)+1
605
 
C     
606
 
C     RETURN FOR WEIGHTED MASSLESS MOMENTA
607
 
 209  IF(NM.NE.0) GOTO 210
608
 
C     RETURN LOG OF WEIGHT
609
 
      WT=WT
610
 
      RETURN
611
 
C     
612
 
C     MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
613
 
 210  XMAX=SQRT(1.-(XMT/ET)**2)
614
 
      DO 301 I=1,N
615
 
      XM2(I)=XM(I)**2
616
 
 301  P2(I)=P(4,I)**2
617
 
      ITER=0
618
 
      X=XMAX
619
 
      ACCU=ET*ACC
620
 
 302  F0=-ET
621
 
      G0=0.
622
 
      X2=X*X
623
 
      DO 303 I=1,N
624
 
      E(I)=SQRT(XM2(I)+X2*P2(I))
625
 
      F0=F0+E(I)
626
 
 303  G0=G0+P2(I)/E(I)
627
 
      IF(ABS(F0).LE.ACCU) GOTO 305
628
 
      ITER=ITER+1
629
 
      IF(ITER.LE.ITMAX) GOTO 304
630
 
      PRINT 1006,ITMAX
631
 
      GOTO 305
632
 
 304  X=X-F0/(X*G0)
633
 
      GOTO 302
634
 
 305  DO 307 I=1,N
635
 
      V(I)=X*P(4,I)
636
 
      DO 306 K=1,3
637
 
 306  P(K,I)=X*P(K,I)
638
 
 307  P(4,I)=E(I)
639
 
C     
640
 
C     CALCULATE THE MASS-EFFECT WEIGHT FACTOR
641
 
      WT2=1.
642
 
      WT3=0.
643
 
      DO 308 I=1,N
644
 
      WT2=WT2*V(I)/E(I)
645
 
 308  WT3=WT3+V(I)**2/E(I)
646
 
      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
647
 
C     
648
 
C     RETURN FOR  WEIGHTED MASSIVE MOMENTA
649
 
      WT=WT+WTM
650
 
      IF(WT.GE.-180.D0) GOTO 309
651
 
      IF(IWARN(3).LE.5) PRINT 1004,WT
652
 
      IWARN(3)=IWARN(3)+1
653
 
 309  IF(WT.LE. 174.D0) GOTO 310
654
 
      IF(IWARN(4).LE.5) PRINT 1005,WT
655
 
      IWARN(4)=IWARN(4)+1
656
 
C     RETURN LOG OF WEIGHT
657
 
 310  WT=WT
658
 
      RETURN
659
 
C     
660
 
 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
661
 
 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',' SMALLER'
662
 
     $ //' THAN TOTAL ENERGY =',D15.6)
663
 
 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
664
 
 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
665
 
 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
666
 
     $ ' DESIRED ACCURACY =',D15.6)
667
 
      END
668
 
 
669
 
      FUNCTION RN(IDUMMY)
670
 
      REAL*8 RN,RAN
671
 
      SAVE INIT
672
 
      DATA INIT /1/
673
 
      IF (INIT.EQ.1) THEN
674
 
        INIT=0
675
 
        CALL RMARIN(1802,9373)
676
 
        END IF
677
 
C       
678
 
 10     CALL RANMAR(RAN)
679
 
        IF (RAN.LT.1D-16) GOTO 10
680
 
        RN=RAN
681
 
C       
682
 
        END
683
 
 
684
 
 
685
 
 
686
 
        SUBROUTINE RANMAR(RVEC)
687
 
C       -----------------
688
 
C       Universal random number generator proposed by Marsaglia and
689
 
C        Zaman
690
 
C       in report FSU-SCRI-87-50
691
 
C       In this version RVEC is a double precision variable.
692
 
        IMPLICIT REAL*8(A-H,O-Z)
693
 
        COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
694
 
        COMMON/ RASET2 / IRANMR,JRANMR
695
 
        SAVE /RASET1/,/RASET2/
696
 
        UNI = RANU(IRANMR) - RANU(JRANMR)
697
 
        IF(UNI .LT. 0D0) UNI = UNI + 1D0
698
 
        RANU(IRANMR) = UNI
699
 
        IRANMR = IRANMR - 1
700
 
        JRANMR = JRANMR - 1
701
 
        IF(IRANMR .EQ. 0) IRANMR = 97
702
 
        IF(JRANMR .EQ. 0) JRANMR = 97
703
 
        RANC = RANC - RANCD
704
 
        IF(RANC .LT. 0D0) RANC = RANC + RANCM
705
 
        UNI = UNI - RANC
706
 
        IF(UNI .LT. 0D0) UNI = UNI + 1D0
707
 
        RVEC = UNI
708
 
        END
709
 
 
710
 
        SUBROUTINE RMARIN(IJ,KL)
711
 
C       -----------------
712
 
C       Initializing routine for RANMAR, must be called before
713
 
C        generating
714
 
C       any pseudorandom numbers with RANMAR. The input values should
715
 
C        be in
716
 
C       the ranges 0<=ij<=31328 ; 0<=kl<=30081
717
 
        IMPLICIT REAL*8(A-H,O-Z)
718
 
        COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
719
 
        COMMON/ RASET2 / IRANMR,JRANMR
720
 
        SAVE /RASET1/,/RASET2/
721
 
C       This shows correspondence between the simplified input seeds
722
 
C        IJ, KL
723
 
C       and the original Marsaglia-Zaman seeds I,J,K,L.
724
 
C       To get the standard values in the Marsaglia-Zaman paper
725
 
C        (i=12,j=34
726
 
C       k=56,l=78) put ij=1802, kl=9373
727
 
        I = MOD( IJ/177 , 177 ) + 2
728
 
        J = MOD( IJ     , 177 ) + 2
729
 
        K = MOD( KL/169 , 178 ) + 1
730
 
        L = MOD( KL     , 169 )
731
 
        DO 300 II = 1 , 97
732
 
        S =  0D0
733
 
        T = .5D0
734
 
        DO 200 JJ = 1 , 24
735
 
        M = MOD( MOD(I*J,179)*K , 179 )
736
 
        I = J
737
 
        J = K
738
 
        K = M
739
 
        L = MOD( 53*L+1 , 169 )
740
 
        IF(MOD(L*M,64) .GE. 32) S = S + T
741
 
        T = .5D0*T
742
 
 200    CONTINUE
743
 
        RANU(II) = S
744
 
 300    CONTINUE
745
 
        RANC  =   362436D0 / 16777216D0
746
 
        RANCD =  7654321D0 / 16777216D0
747
 
        RANCM = 16777213D0 / 16777216D0
748
 
        IRANMR = 97
749
 
        JRANMR = 33
750
 
        END
751
 
 
752
 
 
753
 
 
754
 
 
755
 
 
756
 
 
757