← Back to branch summary

~alifson/chiralityflow/trunk

« back to all changes in this revision

Viewing changes to tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa.f

  • Committer: andrew.lifson at lu
  • Date: 2021-09-01 15:34:39 UTC
  • Revision ID: andrew.lifson@thep.lu.se-20210901153439-7fasjhav4cp4m88r
testing a new repository of a madgraph folder

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_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 QED=2)'
 
257
          WRITE (*,*) '---------------------------------'
 
258
          WRITE(*,*) 'All loop contributions are of split orders'
 
259
     $     //' (QCD=2 QED=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 2'
 
365
          WRITE (69,*) 'SO_Born BORN ',MATELEM(0,1)
 
366
          WRITE (69,*) 'Split_Orders_Names QCD QED'
 
367
          WRITE (69,*) 'Loop_SO_Results 2 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