Removing warnings (Sun)
[u/mrichter/AliRoot.git] / TPHIC / tphic17.f
1 CDECK  ID>, GGCDES. 
2 *=======================================================================
3 CDECK  ID>, GGINIT. 
4       SUBROUTINE gginit
5 *-----------------------------------------------------------------------
6 *       GGHIC initialization
7 *       Author: Yu.Kharlov
8 *-----------------------------------------------------------------------
9       COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
10      &               amin,amax,ymin,ymax,nmas,ny, kferm,
11      &               kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
12      &               thetamin, costhv1, kv1,kv2,gvpar(4)
13       CHARACTER lumfil*80
14       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
15      &               gvconst(4,10)
16       REAL nc
17       COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
18      &  xstote, ssbr(10)
19       COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
20      &                ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
21       COMMON /ggmssm/ xm1,   xm2,    xmg,    xms,    xmtl,   xmtr,
22      &                xmll,  xmlr,   xmnl,   xtanb,  xmha,   xmu,
23      &                xmt,   xat,    xmbr,   xab,    u11,    v11
24       COMMON/D2LParam/lz,la,Rion,Gamma
25       INTEGER lz,la
26       REAL    Rion,Gamma
27       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
28       DOUBLE PRECISION CKIN
29       SAVE  /PYSUBS/
30       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
31       DOUBLE PRECISION PARP,PARI
32       SAVE /PYPARS/
33       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
34       DOUBLE PRECISION P,V
35       SAVE  /PYJETS/
36       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
37       DOUBLE PRECISION PARU,PARJ
38       SAVE  /PYDAT1/
39       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
40       DOUBLE PRECISION PMAS,PARF,VCKM
41       SAVE  /PYDAT2/
42       COMMON/PYDAT4/CHAF(500,2)
43       CHARACTER CHAF*16
44       SAVE  /PYDAT4/
45 C          MXSS                 = maximum number of modes
46 C          NSSMOD               = number of modes
47 C          ISSMOD               = initial particle
48 C          JSSMOD               = final particles
49 C          GSSMOD               = width
50 C          BSSMOD               = branching ratio
51       INTEGER MXSS
52       PARAMETER (MXSS=1000)
53       COMMON/SSMODE/NSSMOD,ISSMOD(MXSS),JSSMOD(5,MXSS),GSSMOD(MXSS)
54      $,BSSMOD(MXSS)
55       INTEGER NSSMOD,ISSMOD,JSSMOD
56       REAL GSSMOD,BSSMOD
57       SAVE /SSMODE/
58 C          SUSY parameters
59 C          AMGLSS               = gluino mass
60 C          AMQKSS               = squark mass
61 C          AMLLSS               = left-slepton mass
62 C          AMLRSS               = right-slepton mass
63 C          AMNLSS               = sneutrino mass
64 C          TWOM1                = Higgsino mass = - mu
65 C          RV2V1                = ratio v2/v1 of vev's
66 C          AMTLSS,AMTRSS        = left,right stop masses
67 C          AMT1SS,AMT2SS        = light,heavy stop masses
68 C          AMBLSS,AMBRSS        = left,right sbottom masses
69 C          AMB1SS,AMB2SS        = light,heavy sbottom masses
70 C          AMZiSS               = signed mass of Zi
71 C          ZMIXSS               = Zi mixing matrix
72 C          AMWiSS               = signed Wi mass
73 C          GAMMAL,GAMMAR        = Wi left, right mixing angles
74 C          AMHL,AMHH,AMHA       = neutral Higgs h0, H0, A0 masses
75 C          AMHC                 = charged Higgs H+ mass
76 C          ALFAH                = Higgs mixing angle
77 C          AAT                  = stop trilinear term
78 C          THETAT               = stop mixing angle
79 C          AAB                  = sbottom trilinear term
80 C          THETAB               = sbottom mixing angle
81       COMMON/SSPAR/AMGLSS,AMULSS,AMURSS,AMDLSS,AMDRSS,AMSLSS
82      $,AMSRSS,AMCLSS,AMCRSS,AMBLSS,AMBRSS,AMB1SS,AMB2SS
83      $,AMTLSS,AMTRSS,AMT1SS,AMT2SS,AMELSS,AMERSS,AMMLSS,AMMRSS
84      $,AMLLSS,AMLRSS,AML1SS,AML2SS,AMN1SS,AMN2SS,AMN3SS
85      $,TWOM1,RV2V1,AMZ1SS,AMZ2SS,AMZ3SS,AMZ4SS,ZMIXSS(4,4)
86      $,AMW1SS,AMW2SS
87      $,GAMMAL,GAMMAR,AMHL,AMHH,AMHA,AMHC,ALFAH,AAT,THETAT
88      $,AAB,THETAB,AAL,THETAL,AMGVSS
89       REAL AMGLSS,AMULSS,AMURSS,AMDLSS,AMDRSS,AMSLSS
90      $,AMSRSS,AMCLSS,AMCRSS,AMBLSS,AMBRSS,AMB1SS,AMB2SS
91      $,AMTLSS,AMTRSS,AMT1SS,AMT2SS,AMELSS,AMERSS,AMMLSS,AMMRSS
92      $,AMLLSS,AMLRSS,AML1SS,AML2SS,AMN1SS,AMN2SS,AMN3SS
93      $,TWOM1,RV2V1,AMZ1SS,AMZ2SS,AMZ3SS,AMZ4SS,ZMIXSS
94      $,AMW1SS,AMW2SS
95      $,GAMMAL,GAMMAR,AMHL,AMHH,AMHA,AMHC,ALFAH,AAT,THETAT
96      $,AAB,THETAB,AAL,THETAL,AMGVSS
97       REAL AMZISS(4)
98       EQUIVALENCE (AMZISS(1),AMZ1SS)
99       SAVE /SSPAR/
100       REAL xpar(10),ypar(6)
101       INTEGER ipar(2)
102       INTEGER pycomp
103       REAL*8 ggrnd
104       EXTERNAL pydata
105       la = ia
106       lz = iz
107 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108 *
109 *       Check bounds on 2 photon mass
110 *
111       IF (iproc .EQ. 1) THEN
112         IF (amin .LE. 2.31) THEN
113           WRITE (6,*)
114      &      '  %GGINIT: minimal 2-photon mass is too low, ',
115      &      'set to 3m(pho)'
116           amin = 2.31
117         END IF
118 *
119       ELSE IF (iproc .EQ. 2) THEN
120         kc = pycomp (kf_onium)
121         xmres = pmas(kc,1)
122         IF (kf_onium .EQ. 441) THEN
123           xgtres = 0.013
124         ELSE
125           xgtres = pmas(kc,2)
126         END IF
127 c.        xggres = ggwid(kf_onium)
128         IF (xgtres .EQ. 0.) xgtres = 0.1
129         IF (xmres.LE.amin .OR. xmres.GE.amax) THEN
130           WRITE(6,*)
131      &      '  %GGINIT: Xmres is out of the Luminosity mass range'
132           STOP
133         ENDIF
134 *
135       ELSE IF (iproc .EQ. 3) THEN
136         kc = pycomp (kferm)
137         amprt(1) = pmas(kc,1)
138         IF (amin .LT. 2*amprt(1)) THEN
139           WRITE (6,*)
140      &      '  %GGINIT: minimal 2-photon mass is too low, ',
141      &      'set to 2m(fermion)'
142           amin = 2*amprt(1)
143         END IF
144 *
145       ELSE IF (iproc .EQ. 4) THEN
146         IF (amprt(4) .LE. 0.) amprt(4) = pmas(pycomp(24),1)
147         IF (amin .LT. 2*amprt(4)) then
148           WRITE (6,*)
149      &      '  %GGINIT: minimal 2-photon mass is too low, ',
150      &      'set to 2m(W+-)'
151           amin = 2*amprt(4)
152         END IF
153 *
154       ELSE IF (iproc .EQ. 5) THEN
155          WRITE (6,*) 'Process 5 is not implemented yet'
156          STOP
157 C         CALL ssmssm (xm1, xm2, xmg, xms, xmtl, xmtr, xmll, xmlr, xmnl,
158 C      &               xtanb, xmha, xmu, xmt, xat, xmbr, xab, iallow)
159 C         IF (iallow .NE. 0) THEN
160 C           WRITE (6,*) 'Lightest neutralino is not LSP, ',
161 C      &                'input other MSSM parameters'
162 C           STOP
163 C         END IF
164 C         amprt(1)  = abs (amw1ss)
165 C         amprt(2)  = abs (amz1ss)
166 C         WRITE (6,'(''   %GGINIT'','//
167 C      &            ''', M(chi+) = '',f5.1,'//
168 C      &            ''', M(chi0) = '',f5.1)') amprt(1), amprt(2)
169 C *
170 C         IF (amin .LT. 2*amprt(1)) THEN
171 C           WRITE (6,*)
172 C      &      '  %GGINIT: minimal 2-photon mass is too low, ',
173 C      &      'set to 2M(chi+)'
174 C           amin = 2*amprt(1)
175 C         END IF
176 *
177       ELSE IF (iproc .EQ. 6) THEN
178         kc1 = pycomp (kv1)
179         amprt(1) = pmas(kc1,1) + pmas(kc1,3)
180         kc2 = pycomp (kv2)
181         amprt(2) = pmas(kc2,1) + pmas(kc2,3)
182         IF (amin .LT. amprt(1)+amprt(2)) THEN
183           WRITE (6,*)
184      &      '  %GGINIT: minimal 2-photon mass is too low, ',
185      &      'set to m(kv1)+m(kv2)'
186           amin = amprt(1)+amprt(2)
187         END IF
188 *
189       END IF
190 *
191       IF (amin .GT. amax) THEN
192         WRITE (6,*) '  %GGINIT: AMIN > AMAX, execution stopped'
193         STOP
194       END IF
195 *
196 *  -----  TWOGAM Initialization  -----
197 *
198       gamma = ebmn*ia/amas
199       xpar(1)= ebmn*ia  !  Beam energy
200       xpar(2)= amin     !  Minimum gg mass
201       xpar(3)= amax     !  Maximum gg mass
202       xpar(8)= 100.     !  Weight, but not too long
203 *
204       ipar(1)= 1        !  Unweighted events
205       ipar(2)= 0        !  Continuum generation
206 *
207       IF (iproc .EQ. 2) THEN
208         xpar(9) = xmres     !  Resonance mass
209         xpar(10)= xgtres    !  Resonance total width
210         ipar(2) = 1         !  Resonance generation
211       END IF
212 *
213       ypar(1)= iz       !  Ion charge
214       ypar(2)= ia       !  Ion atomic number
215       ypar(3)= amas     !  Ion mass
216       ypar(4)= gamma    !  Gamma factor of ion
217       ypar(5)= ymin     !  Minimum gg-rapidity
218       ypar(6)= ymax     !  Maximum gg-rapidity
219       eb = xpar(1)
220 *
221       CALL twoini (6,ipar,xpar,ypar)
222 *
223 *       The Luminosity function initialization
224 *
225 *     ilumf = +1         !   Reading the Luminosity from file
226 *     ilumf = -1         !   New calculation of the Luminosity function
227 *
228       CALL lumfun (ilumf,lumfil,ymin,ymax,ny,amin,amax,nmas)
229 *
230 *       Integral of lum. function over the (M,y) region
231       xlumint = 0.0
232       hmas = (amax-amin)/nmas
233       f1   = dldmx(amin   ,ny)
234       DO i=1,nmas
235         xm  = amin + i*hmas
236         f2  = dldmx(xm-hmas/2.,ny)
237         f3  = dldmx(xm        ,ny)
238         xlumint = xlumint +(f1 + 4.*f2 + f3)/6.
239         f1 = f3
240       END DO
241       xlumint = xlumint * hmas
242 *  ----- end of TWOGAM initialization -----
243 *
244       IF (iproc .EQ. 1) THEN
245 *  ----- PYTHIA initialization -----
246         parp(2) = dble(amin)
247         mstp(14) = 20
248         mstp(82) = 1
249         mstp(171) = 1
250         mstp(172) = 1
251         DO ip = 1,2
252           p(1,ip) = 0.D0
253           p(2,ip) = 0.D0
254         END DO
255         p(1,3) = dble( amax/2)
256         p(2,3) = dble(-amax/2)
257         p(1,4) = dble( amax/2)
258         p(2,4) = dble( amax/2)
259         p(1,5) =  0.D0
260         p(2,5) =  0.D0
261         CALL pyinit ('5MOM','gamma','gamma',0.0D0)
262 *  ----- end of PYTHIA initialization -----
263 *
264       ELSE IF (iproc .EQ. 2) THEN
265 *   ---  Cross section for narrow Resonance production  ---
266         totxsc = dldmx(xmres,ny)                    !  dL/dMx in GeV^-1
267         njres = kf_onium / 10
268         njres = kf_onium - 10 * njres
269         xnorm = 4. * pi**2 * njres * xggres * gev2nb / xmres**2
270         xstot = xnorm * totxsc
271         xstote = 0.
272 *
273       ELSE IF      (iproc .EQ. 3) THEN
274         qf = abs(kchg(kc,1))/3.
275         nc = abs(kchg(kc,2))*2. + 1.
276       ELSE IF (iproc .EQ. 5) THEN
277 *           Chargino storing in JETSET
278         kc = pycomp (41)
279         chaf(kc,1) = 'chi_1'
280         kchg(kc,1) = -3
281         kchg(kc,2) =  0
282         kchg(kc,3) =  1
283         pmas(kc,1) = amprt(1)
284         qf = 1.
285         nc = 1.
286 *           Neutralino storing in JETSET
287         kc = pycomp (42)
288         chaf(kc,1) = 'chi_1'
289         kchg(kc,1) =  0
290         kchg(kc,2) =  0
291         kchg(kc,3) =  1
292         pmas(kc,1) = amprt(2)
293 *
294 *           Gaugino mixing parameters
295         u11 = -zmixss(4,1)*sin(gammar)/sqrt(2.) +
296      &         zmixss(2,1)*cos(gammar)
297         v11 =  zmixss(3,1)*sin(gammal)/sqrt(2.) +
298      &         zmixss(2,1)*cos(gammal)
299         DO idc = 1,nssmod
300           IF (issmod(idc) .EQ. 39) THEN
301             IF (moddcy .EQ. 1) THEN
302               DO jmod = 1,5
303                 IF (jssmod(jmod,idc) .EQ. -12) THEN
304                   xsbra = bssmod(idc)
305                   xsbra = xsbra*xsbra*4
306                   GOTO 10
307                 END IF
308               END DO
309             ELSE
310               DO jmod = 1,5
311                 IF (jssmod(jmod,idc) .EQ. -12) THEN       ! e nu_e
312                   ssbr(1) = bssmod(idc)
313                   GOTO 20
314                 ELSE IF (jssmod(jmod,idc) .EQ. -14) THEN  ! mu nu_mu
315                   ssbr(2) = bssmod(idc)
316                   GOTO 20
317                 ELSE IF (jssmod(jmod,idc) .EQ. -16) THEN  ! tau nu_tau
318                   ssbr(3) = bssmod(idc)
319                   GOTO 20
320                 ELSE IF (jssmod(jmod,idc) .EQ. -2) THEN   ! u dbar
321                   ssbr(4) = bssmod(idc)
322                   GOTO 20
323                 ELSE IF (jssmod(jmod,idc) .EQ.  4) THEN   ! c sbar
324                   ssbr(5) = bssmod(idc)
325                   GOTO 20
326                 END IF
327               END DO
328    20         CONTINUE
329             END IF
330           END IF
331         END DO
332 *
333    10   CONTINUE
334       ELSE IF (iproc .EQ. 6) THEN
335 *         Fix parameters of process gg -> V1V2
336         IF (kv1 .GT. kv2) THEN
337           kvtmp = kv1
338           kv1   = kv2
339           kv2   = kvtmp
340         END IF
341         IF      (kv1.EQ.113 .AND. kv2.EQ.113) THEN
342           idx_vv = 1
343         ELSE IF (kv1.EQ.113 .AND. kv2.EQ.223) THEN
344           idx_vv = 2
345         ELSE IF (kv1.EQ.113 .AND. kv2.EQ.333) THEN
346           idx_vv = 3
347         ELSE IF (kv1.EQ.113 .AND. kv2.EQ.443) THEN
348           idx_vv = 4
349         ELSE IF (kv1.EQ.223 .AND. kv2.EQ.223) THEN
350           idx_vv = 5
351         ELSE IF (kv1.EQ.223 .AND. kv2.EQ.333) THEN
352           idx_vv = 6
353         ELSE IF (kv1.EQ.223 .AND. kv2.EQ.443) THEN
354           idx_vv = 7
355         ELSE IF (kv1.EQ.333 .AND. kv2.EQ.333) THEN
356           idx_vv = 8
357         ELSE IF (kv1.EQ.333 .AND. kv2.EQ.443) THEN
358           idx_vv = 9
359         ELSE IF (kv1.EQ.443 .AND. kv2.EQ.443) THEN
360           idx_vv = 10
361         ELSE
362           WRITE (6,*) '  %GGINIT: unknown (V1 V2) combination'
363           WRITE (6,'(3x,''KV1 = '',i3,'', KV2 = '',i3)') kv1, kv2
364           STOP
365         END IF
366         DO 101 igvpar=1,4
367   101   gvpar(igvpar) = gvconst(igvpar,idx_vv)
368         IF (gvpar(1) .EQ. 0.0) THEN
369           WRITE (6,*) '  %GGINIT: (V1 V2) combination not implemented'
370           WRITE (6,'(3x,''KV1 = '',i3,'', KV2 = '',i3)') kv1, kv2
371           STOP
372         END IF
373       END IF
374 *
375       IF (iproc.EQ.1 .OR. iproc.EQ.3 .OR.
376      &    iproc.EQ.4 .OR. iproc.EQ.5 .OR. iproc.EQ.6) THEN
377 *         Max cross section estimation
378         xsmax0 = 0.
379         DO nest = 1,100
380           xm  = amin + ggrnd(0) * (amax-amin)
381           wsq = xm*xm
382           CALL ggxsec
383         END DO
384         xssum = 0.
385         ntry  = 0
386       END IF
387 *
388       RETURN
389       END
390 *
391 *=======================================================================
392 CDECK  ID>, GGWID.  
393       FUNCTION ggwid (kf_res)
394 *-----------------------------------------------------------------------
395 *       Routine to define 2-gamma width of a resonance
396 *       with JETSET code KF_RES
397 *       Author: Yu.Kharlov
398 *-----------------------------------------------------------------------
399       REAL wid(5)
400       DATA wid /9.E-09, 1.E-06, 5.E-06, 6.3E-06, 0.41E-06/
401       IF      (kf_res .EQ. 111) THEN
402         ggwid = wid(1)
403       ELSE IF (kf_res .EQ. 221) THEN
404         ggwid = wid(2)
405       ELSE IF (kf_res .EQ. 331) THEN
406         ggwid = wid(3)
407       ELSE IF (kf_res .EQ. 441
408      &    .OR. kf_res .EQ. 10441
409      &    .OR. kf_res .EQ. 445) THEN
410         ggwid = wid(4)
411       ELSE IF (kf_res .EQ. 551
412      &    .OR. kf_res .EQ. 10551
413      &    .OR. kf_res .EQ. 555) THEN
414         ggwid = wid(5)
415       END IF
416       RETURN
417       END
418 *
419 *=======================================================================
420 CDECK  ID>, GGRUN.  
421       SUBROUTINE ggrun
422 *-----------------------------------------------------------------------
423 *       Single event generation
424 *       Author: Yu.Kharlov
425 *-----------------------------------------------------------------------
426       COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
427      &               amin,amax,ymin,ymax,nmas,ny, kferm,
428      &               kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
429      &               thetamin, costhv1, kv1,kv2,gvpar(4)
430       CHARACTER lumfil*80
431       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
432      &               gvconst(4,10)
433       REAL nc
434       COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
435      &                ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
436       COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
437      &  xstote, ssbr(10)
438       COMMON /ggkin/ xkin(10)
439       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
440       DOUBLE PRECISION P,V
441       SAVE  /PYJETS/
442       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
443       DOUBLE PRECISION PARU,PARJ
444       SAVE  /PYDAT1/
445       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
446       DOUBLE PRECISION PMAS,PARF,VCKM
447       SAVE  /PYDAT2/
448       REAL pgam1(4),pgam2(4)
449       REAL p3(4), p4(4)
450       REAL*8 dbetaz
451       REAL*8 ggrnd
452 *
453     1 CALL gg2gam (amas,wsq,ygg,q1sq,q2sq,
454      &             pgam1,pgam2,ptag1,ptag2,weight)
455 *
456 *         Cross section calculation
457       IF (iproc.EQ.1 .OR. iproc.EQ.3 .OR.
458      &    iproc.EQ.4 .OR. iproc.EQ.5 .OR. iproc.EQ.6) THEN
459         CALL ggxsec
460         xssum = xssum + xscur
461         ntry  = ntry  + 1
462 *         Select event on cross section of process gg -> XX
463         IF (xscur0/xsmax0 .LT. sngl (ggrnd(0))) GOTO 1
464 *         Total cross section calculation due to current statistics
465         xstot  = xssum/ntry
466         xstote = xstot / sqrt(float(ntry))
467       END IF
468 *
469       xmgg   = sqrt (wsq)
470       xmg1   = sqrt (q1sq)
471       xmg2   = sqrt (q2sq)
472       dbetaz = dtanh(1.d0*ygg)
473 *
474 *       calculation of the 4-vector of two-photon system
475 *       and filling arrays pgg(1:2,1:5), kgg(1:2)
476       DO ip=1,4
477         p2g(ip)   = pgam1(ip) + pgam2(ip)
478         pgg(1,ip) = pgam1(ip)
479         pgg(2,ip) = pgam2(ip)
480       END DO
481       p2g(5)   =  xmgg
482       pgg(1,5) = -xmg1
483       pgg(2,5) = -xmg2
484       kgg(1)   = 22
485       kgg(2)   = 22
486       ngg      = 2
487 *
488 * =======> PROCESS # 1
489       IF (iproc .EQ. 1) THEN
490         DO ip = 1,4
491           p(1,ip) = dble (pgam1(ip))
492           p(2,ip) = dble (pgam2(ip))
493         END DO
494         p(1,5) = dble (-xmg1)
495         p(2,5) = dble (-xmg2)
496         CALL pyevnt
497 *
498 * =======> PROCESS # 2
499       ELSE IF (iproc .EQ. 2) THEN
500         ngg = ngg + 1
501         DO ip = 1,5
502           pgg(ngg,ip) = p2g(ip)
503         END DO
504         kgg(ngg) =  kf_onium
505 *         Fill the JETSET common block
506         CALL gg2py (ngg)
507         CALL pyexec
508 *
509 * =======> PROCESS # 3, # 4, # 5, # 6
510       ELSE IF (iproc.EQ.3 .OR. iproc.EQ.4 .OR.
511      &         iproc.EQ.5 .OR. iproc.EQ.6) THEN
512         phi    = 2*pi*ggrnd(0)
513         IF (iproc.EQ.3 .OR. iproc.EQ.4 .OR. iproc.EQ.5) THEN
514           e3     = 0.5*sqrt(wsq)
515           e4     = e3
516           beta   = xkin(1)
517           costhe = xkin(2)
518           am1    = xkin(3)
519           am2    = xkin(3)
520           p3abs  = e3 * beta
521         ELSE IF (iproc.EQ.6) THEN
522           am1    = xkin(1)
523           am2    = xkin(2)
524           costhe = xkin(3)
525           p3abs  = xkin(4)
526           e3     = xkin(5)
527           e4     = xmgg - e3
528         END IF
529         sinthe = sqrt (1 - costhe*costhe)
530         p3(1) = p3abs * cos(phi) * sinthe
531         p3(2) = p3abs * sin(phi) * sinthe
532         p3(3) = p3abs * costhe
533         p3(4) = e3
534         p4(4) = e4
535         DO ip = 1,3
536           p4(ip) = -p3(ip)
537         END DO
538         igg1 = ngg + 1
539         igg2 = ngg + 2
540         DO ip = 1,4
541           pgg(igg1,ip) = p3(ip)
542           pgg(igg2,ip) = p4(ip)
543         END DO
544         pgg(igg1,5) = am1
545         pgg(igg2,5) = am2
546         ngg = ngg + 2
547 *
548         IF      (iproc .EQ. 3) THEN
549           kgg(igg1) =  kferm
550           kgg(igg2) = -kferm
551           CALL gg2py (ngg)
552         ELSE IF (iproc .EQ. 4) THEN
553           kgg(igg1) =  24
554           kgg(igg2) = -24
555           CALL gg2py (ngg)
556         ELSE IF (iproc .EQ. 5) THEN
557           kgg(igg1) =  41
558           kgg(igg2) = -41
559 *           Decay of charginos into W+neutralino
560           CALL ggdecy (igg1)
561           CALL ggdecy (igg2)
562 *
563           CALL gg2py (ngg)
564           k(igg1,1) = 11
565           k(igg2,1) = 11
566         ELSE IF (iproc.EQ.6) THEN
567           kgg(igg1) = kv1
568           kgg(igg2) = kv2
569           CALL gg2py (ngg)
570         END IF
571 *
572 *         Boost event into initial frame
573         mstu(33) = 1
574         CALL pyrobo(3,ngg, 0.d0,0.d0, 0.d0,0.d0, dbetaz)
575         CALL pyexec
576 *
577       END IF
578 *
579       RETURN
580       END
581 *
582 *=======================================================================
583 CDECK  ID>, GGEXIT. 
584       SUBROUTINE ggexit
585 *-----------------------------------------------------------------------
586 *       Cross section calculation and end of run and
587 *       print out cross section and related variables
588 *       Author: Yu.Kharlov
589 *-----------------------------------------------------------------------
590       COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
591      &               amin,amax,ymin,ymax,nmas,ny, kferm,
592      &               kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
593      &               thetamin, costhv1, kv1,kv2,gvpar(4)
594       CHARACTER lumfil*80
595       COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
596      &  xstote, ssbr(10)
597       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
598      &               gvconst(4,10)
599       REAL nc
600 *
601       WRITE (6,'(/x,79(''=''))')
602       WRITE (6,'(x,''I'',36x,''TPHIC'',36x,      ''I'')')
603       WRITE (6,'(x,''I'',34x,''Process '',i1,34x,''I'')') iproc
604       WRITE (6,'(x,''I'',34x,''---------'',  34x,''I'')')
605       WRITE (6,'(x,''I'',77x,                    ''I'')')
606       WRITE (6,'(x,''I'',20x,''Events thrown   : '',i8,31x,''I'')')
607      &       ntry
608       WRITE (6,'(x,''I'',20x,''Events accepted : '',i8,31x,''I'')')
609      &       nevent
610       WRITE (6,'(x,''I'',20x,''Cross section   : '',e10.3,'' +- '','//
611      &      'e10.3,'' nb'',12x ''I'')') xstot, xstote
612       WRITE (6,'(x,''I'',77x,                    ''I'')')
613       WRITE (6,'(x,79(''=''))')
614 *
615       IF (iproc .EQ. 1) THEN
616 *   ---  Alternative cross section for minimum bias events ---
617         totxsc = 0.0
618         hmas   = (amax-amin)/nmas
619         f1     = dldmx(amin   ,ny) * dchad(amin)
620         DO i=1,nmas
621           xm  = amin + i*hmas
622           f2  = dldmx(xm-hmas/2.,ny) * dchad(xm-hmas/2.)
623           f3  = dldmx(xm        ,ny) * dchad(xm)
624           totxsc = totxsc +(f1 + 4.*f2 + f3)/6.
625           f1 = f3
626         END DO
627         sigtt = totxsc * hmas
628         sigmx = sigtt/sqrt(float(nevent))
629         WRITE (6,*) ' CrosSec =',sigtt,' +/-',sigmx,'  mb'
630       END IF
631 *
632       RETURN
633       END
634 *
635       REAL FUNCTION dchad(xm)
636 *-----------------------------------------------------------------------
637 *       Total cross section for gamma-gamma collisions vs. Xm in mb
638 *       Parametrization is given by G.A.Schuler, T.Sjostrand,
639 *       CERN-TH.7193/94, formula (9).
640 *       Coded by Yu.Kharlov 20.03.1995.
641 *-----------------------------------------------------------------------
642       DATA epsilon, eta /0.0808, 0.4525/
643       dchad = 0.211e-3 * (xm**2)**epsilon + 0.297e-3/(xm**2)**eta
644       RETURN
645       END
646 *
647 *=======================================================================
648 CDECK  ID>, GGXSEC. 
649       SUBROUTINE ggxsec
650 *-----------------------------------------------------------------------
651 *       Cross section calculation in current event for iproc=1,3,4,5,6
652 *       Author: Yu.Kharlov
653 *-----------------------------------------------------------------------
654       COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
655      &               amin,amax,ymin,ymax,nmas,ny, kferm,
656      &               kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
657      &               thetamin, costhv1, kv1,kv2,gvpar(4)
658       CHARACTER lumfil*80
659       COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
660      &                ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
661       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
662      &               gvconst(4,10)
663       REAL nc
664       COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
665      &  xstote, ssbr(10)
666       COMMON /ggkin/ xkin(10)
667       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
668       DOUBLE PRECISION PMAS,PARF,VCKM
669       SAVE  /PYDAT2/
670       PARAMETER (epsilon = 0.0808, eta = -0.4525)
671       REAL*8 ggrnd
672       REAL*8 pymass
673 *
674       IF (iproc .EQ. 1) THEN
675 *         Parametrization by G.A.Schuler, T.Sjostrand,
676 *         CERN-TH.7193/94, formula (9).
677         xscur0 = 211. * wsq**epsilon + 297. * wsq**eta
678         xscur  = xscur0 * xlumint
679       ELSE IF (iproc .EQ. 3 .OR. iproc .EQ. 5) THEN
680         beta2  = 1.0 - 4.0*amprt(1)**2/wsq
681         beta2  = min  (beta2, 1.0)
682         beta   = sqrt (beta2)
683         IF (beta2 .LT. 0.9) THEN
684     1     costhe = 2.*ggrnd(0) - 1.
685           costh2 = costhe*costhe
686           xnum   = 1. + 2.*beta2 * (1.-beta2) * (1.-costh2) -
687      &             (beta2*costh2)**2
688           xden   = (1 - beta2*costh2)**2
689           ffcur  = xnum / xden
690           ffmax  = (1. + beta2) / (1. - beta2)
691           rff    = ggrnd(0)
692           IF (ffcur .LE. rff*ffmax) GOTO 1
693         ELSE
694           small  = 1.0 - cos(thetamin)
695     2     zz     = exp (ggrnd(0) * alog(2./small))
696 C          sgn    = sign (1., 2.*ggrnd(0)-1.)
697           sgn = 2.*ggrnd(0)-1.
698           sgn    = sign (1., sgn)
699           costhe = sgn * (zz-1.) / (zz+1.)
700           ffcur  = (1.+costhe*costhe) / (1.-beta2*costhe*costhe)
701           ffmax  =                 2. / (1.-      costhe*costhe)
702           rff    = ggrnd(0)
703           IF (ffcur .LE. rff*ffmax) GOTO 2
704         END IF
705 *
706         xkin(1) = beta
707         xkin(2) = costhe
708         xkin(3) = amprt(1)
709 *
710 *           Yu.Kharlov, 3.06.1995
711         totfac = 4. * pi * alpha**2 * gev2nb * qf**4 * nc / wsq
712         IF (beta2 .LT. 0.99) THEN
713           xscur0 =((3.-beta2**2)/(2.*beta) * alog ((1.+beta)/(1.-beta))-
714      &             2. + beta2) * beta
715         ELSE
716           xscur0 = 2. * alog (sqrt(wsq)/amprt(1)) - 1.0
717         END IF
718         xscur0 = totfac * xscur0 * xsbra
719         xscur = xscur0 * xlumint
720       ELSE IF (iproc .EQ. 4) THEN
721         xmw = amprt(4)
722         beta2 = 1. - 4.*xmw**2/wsq
723         beta = sqrt (beta2)
724     3   costhe = 2.*ggrnd(0) - 1.
725         costh2 = costhe*costhe
726         xnum = 19. - 6.*beta2 * (1.-beta2) +
727      &         2.*beta2 * (8. - 3.*beta2) * costh2 +
728      &         3.*(beta2*costh2)**2
729         xden = (1. - beta2*costh2)**2
730         ffcur = xnum / xden
731         ffmax = (19. + 10.*beta2 + beta2**2) / (1. - beta2)**2
732         rff = ggrnd(0)
733         IF (ffcur .LE. rff*ffmax) GOTO 3
734 *
735         xkin(1) = beta
736         xkin(2) = costhe
737         xkin(3) = xmw
738 *
739 *           Yu.Kharlov, 3.06.1995
740         totfac = pi * alpha**2 * beta  * gev2nb/ wsq
741         xscur0 = 2. * (22. - 9.*beta2 + 3.*beta2**2) / (1. - beta2) -
742      &           3. * (1.-beta2**2)/beta * alog ((1.+beta)/(1.-beta))
743         xscur0 = totfac * xscur0 * xsbra
744         xscur  = xscur0 * xlumint
745       ELSE IF (iproc .EQ. 6) THEN
746         xmgg = sqrt (wsq)
747         bb   = gvpar(3) + gvpar(4)*alog(xmgg/50.0)
748         am1  = sngl (pymass(kv1))
749         am2  = sngl (pymass(kv2))
750         p1   = pfork (xmgg,am1,am2)
751         e1   = sqrt (p1*p1 + am1*am1)
752         tmin = am1*am1 - xmgg*(e1+p1*costhv1)
753         tmax = am1*am1 - xmgg*(e1-p1*costhv1)
754         fmin = exp (bb*tmin)
755         fmax = exp (bb*tmax)
756         IF (fmax .GT. 1.e-10) THEN
757           tvar = 1./bb * dlog (fmax - ggrnd(0)*(fmax-fmin))
758         ELSE
759           tvar = tmax + 1./bb * dlog (1. - ggrnd(0))
760         END IF
761         costhe = (e1*xmgg - am1*am1 + tvar) / xmgg / p1
762         IF (costhe .GT. 1.0) costhe =  1.0
763         IF (costhe .LT.-1.0) costhe = -1.0
764         IF (ggrnd(0) .LT. 0.5) costhe = -costhe
765 *
766         xkin(1) = am1
767         xkin(2) = am2
768         xkin(3) = costhe
769         xkin(4) = p1
770         xkin(5) = e1
771 *
772         IF (kv1 .EQ. 443 .AND. kv2 .EQ. 443) THEN
773           factor = 1.0 / alog(xmgg/3.097)
774         ELSE
775           factor = 1.0
776         END IF
777         xscur0 = factor * gvpar(1) * xmgg**gvpar(2) / bb * exp(bb*tmax)
778         xscur  = xscur0 * xlumint
779 *
780       END IF
781 *
782 *       Find max cross section of gg-process
783 *
784       IF (xscur0 .GT. xsmax0) xsmax0 = xscur0
785 *
786       RETURN
787       END
788 *
789 *=======================================================================
790 CDECK  ID>, GGDECY. 
791       SUBROUTINE ggdecy (igg)
792 *-----------------------------------------------------------------------
793 *       Routine to decay a particle number I in GGEVNT common block
794 *       Input  : IGG - particle number in common GGEVNT arrays
795 *       Author : Yu.Kharlov
796 *-----------------------------------------------------------------------
797       COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
798      &               amin,amax,ymin,ymax,nmas,ny, kferm,
799      &               kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
800      &               thetamin, costhv1, kv1,kv2,gvpar(4)
801       CHARACTER lumfil*80
802       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
803      &               gvconst(4,10)
804       REAL nc
805       COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
806      &                ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
807       COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
808      &  xstote, ssbr(10)
809       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
810       DOUBLE PRECISION PMAS,PARF,VCKM
811       SAVE  /PYDAT2/
812       COMMON /ggmssm/ xm1,   xm2,    xmg,    xms,    xmtl,   xmtr,
813      &                xmll,  xmlr,   xmnl,   xtanb,  xmha,   xmu,
814      &                xmt,   xat,    xmbr,   xab,    u11,    v11
815       REAL p1(4), pcm(4), p2(4), p3(4), p4(4), ptmp(3)
816       REAL o1(3,3), o2(3,3)
817       INTEGER pycomp
818       DATA ifl /0/
819       REAL*8 ggrnd
820 *
821       DO ip = 1,4
822         p1(ip) = pgg(igg,ip)
823       END DO
824 *
825       IF (abs(kgg(igg)) .EQ. 41 .AND. iproc .EQ. 3) THEN
826 *         Decay chi_1^+- --> chi_1^0 W^+- in CMS, uniformly
827         eb = (amprt(1)**2 + amprt(2)**2 - amprt(4)**2) / (2 * amprt(1))
828         ec = (amprt(1)**2 - amprt(2)**2 + amprt(4)**2) / (2 * amprt(1))
829         pb = sqrt (eb*eb - amprt(2)*amprt(2))
830         pc = pb
831         costhe = 2*ggrnd(0) - 1.
832         sinthe = sqrt (1 - costhe*costhe)
833         phi = 2*pi*ggrnd(0)
834         pcm(1) = pb * sinthe * cos(phi)
835         pcm(2) = pb * sinthe * sin(phi)
836         pcm(3) = pb * costhe
837         pcm(4) = eb
838         CALL lorenb (amprt(1),p1,pcm,p2)
839         DO ip = 1,3
840           pcm(ip) = -pcm(ip)
841         END DO
842         pcm(ip) = ec
843         CALL lorenb (amprt(1),p1,pcm,p3)
844         DO ip = 1,4
845           pgg(ngg+1,ip) = p2(ip)
846           pgg(ngg+2,ip) = p3(ip)
847         END DO
848         pgg(ngg+1,5) = amprt(2)
849         pgg(ngg+2,5) = amprt(4)
850         kgg(ngg+1) = 42 * sign(1.,1.*kgg(igg))
851         kgg(ngg+2) = 24 * sign(1.,1.*kgg(igg))
852         ngg = ngg + 2
853 *
854       ELSE IF (abs(kgg(igg)) .EQ. 41 .AND. iproc .EQ. 5) THEN
855 *         Decay chi_1^+- --> e^+- nu_e chi_1^0 in CMS,
856 *         according to the exact matrix element
857 *
858         xm23max = amprt(1) - amprt(2)
859         xmw     = sngl (pmas(pycomp(24),1))
860         gwt     = sngl (pmas(pycomp(24),2))
861         atmin   = - atan (xmw/gwt)
862         atmax   =   atan ((xm23max**2 - xmw**2) / (xmw*gwt))
863         am1     = amprt(1)
864         am4     = amprt(2)
865 *
866     1   CONTINUE
867 *
868 *       Random mass of (2+3) according to Breit-Wigner
869 *
870         xm23 = xmw**2 + xmw*gwt *
871      *       tan (atmin + ggrnd(0)*(atmax - atmin))
872         xm23 = max (0., xm23)
873         xm23 = sqrt (xm23)
874         IF (xm23 .GT. xm23max) GOTO 1       ! to avoid boundary effect
875         e23  = (am1**2 + xm23**2 - am4**2) / 2. / am1
876 *
877 *       Rest frame of (2+3)
878 *
879         phi    = 2.*pi * ggrnd(0)
880         costhe = 2.*ggrnd(0) - 1.
881         sinthe = sqrt (1. - costhe*costhe)
882 *
883         pp = xm23 / 2.
884         p2(1) = pp * sinthe * cos(phi)
885         p2(2) = pp * sinthe * sin(phi)
886         p2(3) = pp * costhe
887         p2(4) = pp
888 *
889         DO ip = 1,3
890           p3(ip) = -p2(ip)
891         END DO
892         p3(4) = pp
893 *
894 *       Transformation from rest frame of (2+3) to rest frame of (1),
895 *       3-momenta p2, p3 are along z-axis
896 *
897         gamma = e23 / xm23
898         gambe = sqrt (gamma*gamma - 1.)
899 *
900         p2z = gamma*p2(3) + gambe*p2(4)
901         e2  = gamma*p2(4) + gambe*p2(3)
902         p3z = gamma*p3(3) + gambe*p3(4)
903         e3  = gamma*p3(4) + gambe*p3(3)
904 *
905         p2(3) = p2z
906         p2(4) = e2
907         p3(3) = p3z
908         p3(4) = e3
909 *
910 *       Arbitrary rotation of rest frame of (1)
911 *
912         phi    = 2.*pi * ggrnd(0)
913         costhe = 2.*ggrnd(0) - 1.
914         sinthe = sqrt (1. - costhe*costhe)
915         cosphi = cos(phi)
916         sinphi = sin(phi)
917 *
918         o1(1,1) =  costhe
919         o1(1,2) =  0.
920         o1(1,3) = -sinthe
921         o1(2,1) =  0.
922         o1(2,2) =  1.
923         o1(2,3) =  0.
924         o1(3,1) =  sinthe
925         o1(3,2) =  0.
926         o1(3,3) =  costhe
927 *
928         o2(1,1) =  cosphi
929         o2(1,2) = -sinphi
930         o2(1,3) =  0.
931         o2(2,1) =  sinphi
932         o2(2,2) =  cosphi
933         o2(2,3) =  0.
934         o2(3,1) =  0.
935         o2(3,2) =  0.
936         o2(3,3) =  1.
937 *
938         DO ix = 1,3
939           ptmp(ix) = 0.
940           DO jx = 1,3
941             ptmp(ix) = ptmp(ix) + o1(ix,jx) * p2(jx)
942           END DO
943         END DO
944         DO ix = 1,3
945           p2(ix) = ptmp(ix)
946         END DO
947         DO ix = 1,3
948           ptmp(ix) = 0.
949           DO jx = 1,3
950             ptmp(ix) = ptmp(ix) + o2(ix,jx) * p2(jx)
951           END DO
952         END DO
953         DO ix = 1,3
954           p2(ix) = ptmp(ix)
955         END DO
956 *
957         DO ix = 1,3
958           ptmp(ix) = 0.
959           DO jx = 1,3
960             ptmp(ix) = ptmp(ix) + o1(ix,jx) * p3(jx)
961           END DO
962         END DO
963         DO ix = 1,3
964           p3(ix) = ptmp(ix)
965         END DO
966         DO ix = 1,3
967           ptmp(ix) = 0.
968           DO jx = 1,3
969             ptmp(ix) = ptmp(ix) + o2(ix,jx) * p3(jx)
970           END DO
971         END DO
972         DO ix = 1,3
973           p3(ix) = ptmp(ix)
974         END DO
975 *
976         p12 = am1 * p2(4)
977 *
978         e4 = am1 - e23
979         pp =sqrt (e4*e4 - am4*am4)
980         p4(1) =  pp * sinthe * cos(phi)
981         p4(2) =  pp * sinthe * sin(phi)
982         p4(3) = -pp * costhe
983         p4(4) =  e4
984 *
985         width = (- 4.*p12**2*u11**2      - 4.*p12**2*v11**2 +
986      &           4.*p12*xm23**2*u11**2 - 2.*p12*am4**2*u11**2 -
987      &           2.*p12*am4**2*v11**2  + 2.*p12*am1**2*u11**2 +
988      &           2.*p12*am1**2*v11**2  - xm23**4*u11**2 +
989      &           xm23**2*am4**2*u11**2 - 2.*xm23**2*am4*am1*u11*v11 -
990      &           xm23**2*am1**2*u11**2) /
991      &           ((xm23**2-xmw**2)**2 + (xmw*gwt)**2)
992         width = width * pp * xm23
993 *
994         widmax = (2.*am1*xm23max *
995      &    (2.*(xm23max*u11)**2 + am1**2*(u11**2 + v11**2)) +
996      &    (xm23max*am4*u11)**2) * (am1**2 - am4**2) / 2./am1 *
997      &    xm23 / ((xm23**2-xmw**2)**2 + (xmw*gwt)**2)
998 *
999 *       Select current phase space configuration
1000 *
1001         IF (width .GT. widmax) WRITE (6,*)
1002      &    '%GGDECY: Maximum chi+ width violated'
1003         rwid = widmax * ggrnd(0)
1004         IF (width .LE. rwid) GOTO 1
1005 *
1006 *         Pick up final lepton flavour
1007         rflav = ggrnd(0)
1008         IF (moddcy .EQ. 1) THEN
1009           IF (rflav .LE. 0.5) THEN
1010             lflav1 =  11
1011             lflav2 = -12
1012           ELSE
1013             lflav1 =  13
1014             lflav2 = -14
1015           END IF
1016         ELSE
1017           rr1 = ssbr(1)
1018           rr2 = ssbr(1)+ssbr(2)
1019           rr3 = ssbr(1)+ssbr(2)+ssbr(3)
1020           rr4 = ssbr(1)+ssbr(2)+ssbr(3)+ssbr(4)
1021           rr5 = ssbr(1)+ssbr(2)+ssbr(3)+ssbr(4)+ssbr(5)
1022           IF (rflav .LE. rr1) THEN
1023             lflav1 =  11
1024             lflav2 = -12
1025           ELSE IF (rflav .LE. rr2) THEN
1026             lflav1 =  13
1027             lflav2 = -14
1028           ELSE IF (rflav .LE. rr3) THEN
1029             lflav1 =  15
1030             lflav2 = -16
1031           ELSE IF (rflav .LE. rr4) THEN
1032             lflav1 = -2
1033             lflav2 =  1
1034           ELSE IF (rflav .LE. rr5) THEN
1035             lflav1 = -4
1036             lflav2 =  3
1037           ELSE
1038             GOTO 1
1039           END IF
1040         END IF
1041 *
1042 *         Transformation to CMS of 2-gamma
1043         CALL lorenb (amprt(1),p1,p2,pcm)
1044         DO ip = 1,4
1045           pgg(ngg+1,ip) = pcm(ip)
1046         END DO
1047         pgg(ngg+1,5) = 0.
1048         kgg(ngg+1) = lflav1 * sign(1.,1.*kgg(igg))
1049         ngg = ngg + 1
1050 *
1051         CALL lorenb (amprt(1),p1,p3,pcm)
1052         DO ip = 1,4
1053           pgg(ngg+1,ip) = pcm(ip)
1054         END DO
1055         pgg(ngg+1,5) = 0.
1056         kgg(ngg+1) = lflav2 * sign(1.,1.*kgg(igg))
1057         ngg = ngg + 1
1058 *
1059         CALL lorenb (amprt(1),p1,p4,pcm)
1060         DO ip = 1,4
1061           pgg(ngg+1,ip) = pcm(ip)
1062         END DO
1063         pgg(ngg+1,5) = amprt(2)
1064         kgg(ngg+1) = 42 * sign(1.,1.*kgg(igg))
1065         ngg = ngg + 1
1066       END IF
1067 *
1068       RETURN
1069       END
1070 *
1071 *=======================================================================
1072 CDECK  ID>, GG2PY.  
1073       SUBROUTINE gg2py (nngg)
1074 *-----------------------------------------------------------------------
1075 *       Routine to fill PYTHIA event record by particles from GG
1076 *       Input: NNGG : number of particles
1077 *       Author: Yu.Kharlov
1078 *-----------------------------------------------------------------------
1079       COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
1080      &                ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
1081       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
1082       DOUBLE PRECISION P,V
1083       SAVE  /PYJETS/
1084       DO igg = 1,nngg
1085         IF (igg .LE. 2) THEN
1086           k(igg,1) = 21
1087         ELSE
1088           k(igg,1)=1
1089         END IF
1090         k(igg,2) = kgg(igg)
1091         DO ip = 1,5
1092           p(igg,ip) = pgg(igg,ip)
1093         END DO
1094       END DO
1095       n = nngg
1096       RETURN
1097       END
1098 *
1099 *=======================================================================
1100 CDECK  ID>, PFORK.  
1101       FUNCTION pfork (am0,am1,am2)
1102 *-----------------------------------------------------------------------
1103 *       Routine to find a momentum of a secondary particle in the decay
1104 *       of a particle with mass AM0 in the rest to 2 particles with
1105 *       masses AM1 and AM2.
1106 *-----------------------------------------------------------------------
1107       denom = am0**4 + am1**4 + am2**4 -
1108      &        2.*(am0*am1)**2 - 2.*(am0*am2)**2 - 2.*(am1*am2)**2
1109       IF (denom .GT. 0.) THEN
1110         pfork = sqrt (denom) / 2./am0
1111       ELSE
1112         pfork = -1.
1113       END IF
1114       RETURN
1115       END
1116 *=======================================================================
1117 CDECK  ID>, GG2GAM. 
1118       SUBROUTINE gg2gam (amas,wsq,y3,qs1,qs2,pgam1,pgam2,ptag1,ptag2,wt)
1119 C***********************************************************************
1120 C
1121 C     Generator of the Direct Two Photon Processes in      S.A.Sadovsky
1122 C     coherent heavy ion collisions:    A+B --> 1+2+3       17.01.1995
1123 C     with double differencial gg-luminosity function      Yu.Kharlov
1124 C                                                           20.03.1995
1125 C
1126 C      LUMFUN must be called first to initialise the generator
1127 C
1128 C***********************************************************************
1129 C
1130 C   DETAILED DESCRIPTION
1131 C   ======== ===========
1132 C
1133 C   INPUTS:  only through common /TWOGENC/
1134 C   =======
1135 C
1136 C   OUTPUTS:     AMAS   The mass of ion.
1137 C   ========     WSQ    Square mass of the gamma-gamma system
1138 C                Y3     The rapidity of the gamma-gamma system
1139 C                QS1    Virtual mass squared of photon 1.
1140 C                QS2    Virtual mass squared of photon 2.
1141 C                PGAM1  Four-vector of photon 1
1142 C                PGAM2  Four-vector of photon 2
1143 C                PTAG1  Four-vector of tag 1
1144 C                PTAG2  Four-vector of tag 2
1145 C                WT     The weight of the event (dummy so far)
1146 C
1147 C   The remaining variables are internal to the program.
1148 C
1149 C***********************************************************************
1150       INTEGER   icnter,ncntrs,nevts,lout,debug
1151       PARAMETER (ncntrs =7)
1152       COMMON /twgenc/s,eb,wmin,wmax,th1pmn,th1pmx,th2pmn,th2pmx,
1153      +   st12mn,st12mx,st12rt,st22mn,st22mx,st22rt,ommin,omrat,eps,
1154      +   wghtsm,wghtgt,wghtmx,wtmxfd,xmres,xgres,
1155      +   icnter(ncntrs),nevts,lout,lres,lwate,debug
1156       SAVE   /twgenc/
1157       REAL      st12mn,st12mx,st12rt,st22mn,st22mx,st22rt
1158       REAL      wmin,wmax,ommin,omrat,eps,th1pmn,th1pmx
1159       REAL      th2pmn,th2pmx,s,eb,wghtsm,wghtmx,wghtgt,wtmxfd
1160       REAL      xmres,xgres
1161       LOGICAL   lres,lwate
1162 C
1163       REAL       alpha,pi
1164       PARAMETER (alpha=1./137.036,pi=3.14159265)
1165 C
1166       REAL ptag1(10),ptag2(10),pgam1(10),pgam2(10),wt
1167       REAL amas,wsq
1168 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1169       REAL*4  y3,xm,q0,qs,qs1,qs2
1170       REAL*8  pa(5),pb(5)
1171 C
1172       REAL*8  pi2
1173 C
1174       LOGICAL     lstr
1175       DATA q0,pi2,lstr /  0.060, 6.28318530718 d0, .false. /
1176       REAL*8 ggrnd
1177 C
1178   999 CONTINUE
1179       icnter(1) = icnter(1) + 1
1180 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1181       IF    (lstr) go to 10
1182       lstr =.true.
1183       qs   = q0/sqrt(2.)
1184 C
1185       pa(1)= 0D0
1186       pa(2)= 0D0
1187       pa(3)= dsqrt(1D0*(eb*eb-amas*amas))
1188       pa(4)= eb*1D0
1189       pa(5)= amas*1D0
1190 C
1191       pb(1)= 0D0
1192       pb(2)= 0D0
1193       pb(3)=-pa(3)
1194       pb(4)= pa(4)
1195       pb(5)= amas*1D0
1196 C
1197    10 xm   = xmres
1198 *       Generate 2gamma rapidity y3 and mass xm
1199       CALL ranlum(y3,xm,wt)
1200 C
1201       wsq  = xm*xm
1202 C
1203 *       Energy and z-momentum of 2gamma system
1204       egg = xm*cosh(y3)
1205       pgg = xm*sinh(y3)
1206 *
1207 *       Squared masses of the photons
1208       qs1 = - qs**2 * dlog(ggrnd(0))    ! is it true that squared mass
1209       qs2 = - qs**2 * dlog(ggrnd(0))    ! is distributed in exponent?
1210       xm1 = -qs1
1211       xm2 = -qs2
1212 *
1213 *       Energies and z-momenta of the photons
1214       xfun12 = wsq + xm1 - xm2
1215 *
1216       eg1 = (egg + pgg*sqrt(1 - 4*wsq*xm1/(xfun12*xfun12))) *
1217      &      xfun12/(2*wsq)
1218       eg2 = egg - eg1
1219       pg1 =  sqrt (eg1**2 - xm1)
1220       pg2 = -sqrt (eg2**2 - xm2)
1221 *
1222 *       4-momenta of the photons
1223       DO ip=1,2
1224         pgam1(ip) = 0.
1225         pgam2(ip) = 0.
1226       END DO
1227       pgam1(3) = pg1
1228       pgam1(4) = eg1
1229       pgam2(3) = pg2
1230       pgam2(4) = eg2
1231 *
1232 *       4-momenta of the recoil ions
1233       DO ip=1,4
1234         ptag1(ip) = pa(ip) - pgam1(ip)
1235         ptag2(ip) = pb(ip) - pgam2(ip)
1236       END DO
1237 *
1238       RETURN
1239       END
1240 *
1241 *=======================================================================
1242 CDECK  ID>, RANLUM. 
1243       SUBROUTINE ranlum (ry,rm,wlum)
1244 C                                                     G.V.Khaoustov
1245 C                                                       12.12.1994
1246 C                                             Corrected: Yu.Kharlov
1247 C                                                       23.05.1995
1248 C                                                       24-JUN-1997
1249 C
1250 C generate two random numbers (ry and rm) with probability proportional
1251 C to luminosity function in ymin-ymax rapidity and xmin-xmax mass range
1252 C
1253       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1254      &               gvconst(4,10)
1255       REAL nc
1256       PARAMETER (isize=100000)
1257       COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
1258      +                           xmin,xmax,amn,amx,nnms,alfun(isize)
1259       PARAMETER (NCNTRS =7)
1260       COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX,
1261      +   ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS,
1262      +   WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES,
1263      +   ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG
1264       LOGICAL lres,lwate
1265       REAL*4  wlum
1266       DATA    ifl/0/, stpy/0./, stpm/0./
1267       REAL*8 ggrnd
1268 C
1269       IF (ifl.EQ.0) THEN
1270         IF (ymin.LT.ymn .OR. ymin.GT.ymx .OR.
1271      *      ymax.LT.ymn .OR. ymax.GT.ymx .OR.
1272      *      xmin.LT.amn .OR. xmin.GT.amx .OR.
1273      *      xmax.LT.amn .OR. xmax.GT.amx) THEN
1274           WRITE (6,*) 'Luminosity function is not defined in this range'
1275           STOP
1276         ENDIF
1277         aky = ymax-ymin
1278         akm = xmax-xmin
1279         stpy = (ymx-ymn)/nny
1280         stpm = (amx-amn)/nnms
1281         rmin = exp(cexp*(xmax))
1282         rmax = exp(cexp*(xmin))
1283         dr = rmax-rmin
1284         ifl = 1
1285       ENDIF
1286 C
1287     1 CONTINUE
1288       IF (.NOT.lres) THEN
1289         rm  = ggrnd(0)                     ! random mass
1290         rm  = rmin + rm*dr
1291         rm  = alog(rm)/cexp                ! exponetial law
1292       ELSE
1293         wsq = xmres**2 + xmres*xgres * tan(pi*(ggrnd(0)-0.5))
1294         IF (wsq .LT. xmin*xmin) GOTO 1
1295         IF (wsq .GT. xmax*xmax) GOTO 1
1296         rm  = sqrt (wsq)
1297       END IF
1298 C
1299       ry = ggrnd(0)                        ! random Y
1300 *
1301 *       Find y-limits for generated gg-mass: y_max = log(E_max/M_gg)
1302 *
1303       ym_max  = alog (egg_max/rm)
1304       ygg_max = min (ymax, ym_max)
1305       ygg_min = max (ymin,-ym_max)
1306       ry = ygg_min + ry*(ygg_max-ygg_min)
1307       iy = (ry-ymn)/stpy+1.
1308 C
1309       im = (rm-amn)/stpm
1310       r = anorm*exp(cexp*rm)*ggrnd(0)      ! normalization
1311       p = (ry-(iy-1)*stpy-ymn)/stpy        ! four point interpolation
1312       q = (rm-im*stpm-amn)/stpm
1313       k = iy+im*(nny+1)
1314       f = (1.-p)*(1.-q)*alfun(k)       + p*(1.-q)*alfun(k+1) +
1315      *         q*(1.-p)*alfun(k+nny+1) +      q*p*alfun(k+2+nny)
1316       IF (r .GT. f) GOTO 1
1317       wlum =  1000.*f                      ! differential luminosity in GeV^-1
1318 *
1319       RETURN
1320       END
1321 *
1322 *=======================================================================
1323 CDECK  ID>, GGDATA. 
1324       BLOCK DATA ggdata
1325 *-----------------------------------------------------------------------
1326 *       Some useful constants and masses of MSSM particles
1327 *       Author: Yu.Kharlov
1328 *-----------------------------------------------------------------------
1329       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1330      &               gvconst(4,10)
1331       REAL nc
1332       COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
1333      &               amin,amax,ymin,ymax,nmas,ny, kferm,
1334      &               kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
1335      &               thetamin, costhv1, kv1,kv2,gvpar(4)
1336       CHARACTER lumfil*80
1337       COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
1338      &  xstote, ssbr(10)
1339       COMMON /ggmssm/ xm1,   xm2,    xmg,    xms,    xmtl,   xmtr,
1340      &                xmll,  xmlr,   xmnl,   xtanb,  xmha,   xmu,
1341      &                xmt,   xat,    xmbr,   xab,    u11,    v11
1342       COMMON/D2LParam/lz,la,Rion,Gamma
1343       INTEGER lz,la
1344       REAL    Rion,Gamma
1345       PARAMETER (al = 1./128.)
1346       DATA pi /3.14159276/, alpha /al/, gev2nb /389379./,
1347      &     hbarc /0.197327/
1348       DATA iproc, nevent, ilumf /1, 10000, -1/
1349       DATA lumfil /'gglum.dat'/
1350       DATA Rion /0.0/
1351       DATA amprt /110.,10.,150.,-1.,0.106/
1352       DATA thetamin /1.e-03/
1353       DATA costhv1 /1.0/
1354       DATA xsbra /1./
1355       DATA ny, nmas /50, 50/
1356       DATA xm1,   xm2,    xmg,    xms,    xmtl,   xmtr,
1357      &     xmll,  xmlr,   xmnl,   xtanb,  xmha,   xmu,
1358      &     xmt,   xat,    xmbr,   xab
1359      &     /30.,   50.,   250.,   700.,   600.,   600.,
1360      &     400.,  400.,   400.,   4.,     300.,  -300.,
1361      &     174.,  300.,   300.,   300./
1362       DATA gvconst / 63.0, 0.22, 6.0, 1.0,    ! rho rho
1363      &               14.0, 0.22, 6.0, 1.0,    ! rho omega
1364      &                7.0, 0.22, 4.0, 1.0,    ! rho phi
1365      &               0.05, 0.80, 2.5, 0.0,    ! rho J/psi
1366      &                0.0, 0.00, 0.0, 0.0,    ! omega omega
1367      &                0.8, 0.22, 4.0, 1.0,    ! omega phi
1368      &              6.E-3, 0.80, 2.5, 0.0,    ! omega J/psi
1369      &                0.2, 0.22, 2.0, 1.0,    ! phi phi
1370      &              3.E-3, 0.80, 1.5, 0.0,    ! phi J/psi
1371      &              6.E-5, 2.00, 1.5, 0.0/    ! J/psi J/psi
1372       END
1373 *
1374 *=======================================================================
1375 CDECK  ID>, EVPRINT.
1376       SUBROUTINE evprint(Lun,Iev,P3,Ygg,ifl,Eb,ia,iz,amas,ptag1,ptag2)
1377 *-----------------------------------------------------------------------
1378 *       Event output for internal use by GGHIC authors
1379 *-----------------------------------------------------------------------
1380       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
1381       DOUBLE PRECISION P,V
1382       SAVE  /PYJETS/
1383       REAL ptag1(4),ptag2(4),p3(5)
1384 C
1385       WRITE(lun,100) iev,n,Eb,p3(5),ygg,ifl,-ifl
1386       WRITE(lun,101) 90000+ia,iz,ptag1,amas
1387       WRITE(lun,101) 90000+ia,iz,ptag2,amas
1388       WRITE(lun,101) 90000+ 3, 0,(p3(j),j=1,5)
1389       IF (n.NE.0) THEN
1390         DO i=1,n
1391           ichg=pychge(k(i,2))/3
1392           WRITE(lun,101) k(i,2),ichg,(p(i,j),j=1,5)
1393         ENDDO
1394       ENDIF
1395       RETURN
1396   100 FORMAT(i6,i3,3e16.8,2x,2i8)
1397   101 FORMAT(i6,i3,3(1x,e14.8),e13.8,e13.8)
1398       END
1399 *
1400 *=======================================================================
1401 C
1402 CDECK  ID>, TWOINI. 
1403       SUBROUTINE TWOINI(LOUTX,IPAR,XPAR,YPAR)                           TWOG0965
1404 C***********************************************************************TWOG0966
1405 C                                                                       TWOG0967
1406 C     Initialises the gamma-gamma generator of Ion collisions TWOGAM.   TWOG0968
1407 C                                                                       TWOG0969
1408 C   INPUTS:             LOUTX    Logical unit for print-out             TWOG0970
1409 C   ======= XPAR(1)     EBM      Beam energy                        GeV TWOG0971
1410 C           XPAR(2)     WMN      Minimum two-gamma mass             GeV TWOG0972
1411 C           XPAR(3)     WMX      Maximum two-gamma mass             GeV TWOG0973
1412 C                                                                       TWOG0974
1413 C           XPAR(8)     WGHTMX   Maximum weight                         TWOG0978
1414 C           XPAR(9)     XMRES    Mass of resonance                  GeV TWOG0979
1415 C           XPAR(10)    XGRES    Total width of resonance           GeV TWOG0980
1416 C                                                                       TWOG0981
1417 C           IPAR(1)     0        Unweighted events                      TWOG0982
1418 C                       1        Weighted events                        TWOG0983
1419 C           IPAR(2)     0        Continuum production                   TWOG0984
1420 C                       1        Resonance formation                    TWOG0985
1421 C                                                                       TWOG0986
1422 C   EXTERNALS:   none                                                   TWOG0987
1423 C   ==========                                                          TWOG0988
1424 C                                                                       TWOG0989
1425 C   AUTHOR/DATE:  S.A. Sadovsky,   17 January 1995                      TWOG0990
1426 C   ===========                                                         TWOG0991
1427 C                                                                       TWOG0992
1428 C***********************************************************************TWOG0993
1429 C     IMPLICIT NONE                                                     TWOG0994
1430       INTEGER   ICNTER,NCNTRS,NEVTS,LOUT,DEBUG                          TWOG0995
1431       PARAMETER (NCNTRS =7)                                             TWOG0996
1432       COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX,        TWOG0997
1433      +   ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS,     TWOG0998
1434      +   WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES,                       TWOG0999
1435      +   ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG                     TWOG1000
1436       SAVE   /TWGENC/                                                   TWOG1001
1437       REAL      ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT               TWOG1002
1438       REAL      WMIN,WMAX,OMMIN,OMRAT,EPS,TH1PMN,TH1PMX                 TWOG1003
1439       REAL      TH2PMN,TH2PMX,S,EB,WGHTSM,WGHTMX,WGHTGT,WTMXFD          TWOG1004
1440       REAL      XMRES,XGRES                                             TWOG1005
1441       LOGICAL   LRES,LWATE                                              TWOG1006
1442 C                                                                       TWOG1007
1443       REAL      ALPHA,PI,TWOPI,XME,XMU                                  TWOG1008
1444       PARAMETER (ALPHA=1./137.036,PI=3.14159265)                        TWOG1009
1445       PARAMETER (TWOPI=2.*PI,XME=0.000511,XMU=0.105658)                 TWOG1010
1446 C                                                                       TWOG1011
1447       REAL     XPAR(11),YPAR(6)                                         TWOG1012
1448       INTEGER  LOUTX,IPAR(10),I                                         TWOG1013
1449 C                                                                       TWOG1014
1450       EB    = XPAR(1)                                                   TWOG1015
1451       WMIN  = XPAR(2)                                                   TWOG1016
1452       WMAX  = XPAR(3)                                                   TWOG1017
1453       WGHTMX= XPAR(8)                                                   TWOG1022
1454       WTMXFD= 0                                                         TWOG1023
1455       XMRES = XPAR(9)                                                   TWOG1024
1456       XGRES = XPAR(10)                                                  TWOG1025
1457       LOUT  = LOUTX                                                     TWOG1026
1458       LWATE = IPAR(1) .EQ. 1                                            TWOG1027
1459       LRES  = IPAR(2) .EQ. 1                                            TWOG1028
1460       EPS   = 1./64.*(XME*WMIN**2/EB**3)**2/100.                        TWOG1029
1461 C+                                                                      TWOG1030
1462 C EPS is a small constant. We throw 1/2 COS(th/2)/(eps+SIN(th/2))       TWOG1031
1463 C rather than 1/2 COT(th/2) to be able to start theta from zero.        TWOG1032
1464 C-                                                                      TWOG1033
1465       S     = (2.*EB)**2                                                TWOG1054
1466       OMMIN = WMIN**2/(4.*EB)                                           TWOG1055
1467       OMRAT = EB/OMMIN                                                  TWOG1056
1468 C+                                                                      TWOG1057
1469 C Print out initial conditions.                                         TWOG1058
1470 C-                                                                      TWOG1059
1471       WRITE(LOUT,1000)                                                  TWOG1060
1472  1000 FORMAT(1H ,10X,27('-'),' TWINIT ',27('-'))                        TWOG1061
1473       WRITE(LOUT,1001)                                                  TWOG1062
1474  1001 FORMAT(1H ,10X,'I',60X,'I')                                       TWOG1063
1475 C
1476       WRITE(LOUT,1002)EB,YPAR(3),YPAR(4),YPAR(1),YPAR(2),WMIN,WMAX,
1477      +                                   YPAR(5),YPAR(6)
1478 C                                                                       TWOG1065
1479  1002 FORMAT                                                            TWOG1066
1480      +  (1H ,10X,'I    Beam energy               ',E15.5,               TWOG1067
1481      +  ' GeV',11X,'I',/,                                               TWOG1068
1482      +  1H ,10X,'I    Mass of ion                  ',F12.5,             TWOG1067
1483      +  ' GeV',11X,'I',/,                                               TWOG1068
1484      +  1H ,10X,'I    Gamma of ion                 ',F12.5,             TWOG1067
1485      +  '    ',11X,'I',/,                                               TWOG1068
1486      +  1H ,10X,'I    Charge of ion                ',F12.5,             TWOG1067
1487      +  '    ',11X,'I',/,                                               TWOG1068
1488      +  1H ,10X,'I    Atomic number of ion         ',F12.5,             TWOG1067
1489      +  '    ',11X,'I',/,                                               TWOG1068
1490      +  1H ,10X,'I    Minimum gamma-gamma mass     ',F12.5,             TWOG1069
1491      +  ' GeV',11X,'I',/,                                               TWOG1070
1492      +  1H ,10X,'I    Maximum gamma-gamma mass     ',F12.5,             TWOG1071
1493      +  ' GeV',11X,'I',/,                                               TWOG1080
1494      +  1H ,10X,'I    Minimum rapidity             ',F12.5,             TWOG1069
1495      +  '    ',11X,'I',/,                                               TWOG1070
1496      +  1H ,10X,'I    Maximum rapidity             ',F12.5,             TWOG1071
1497      +  '    ',11X,'I'  )                                               TWOG1080
1498       IF (LRES) THEN                                                    TWOG1081
1499       WRITE(LOUT,1001)                                                  TWOG1082
1500       WRITE(LOUT,1003)XMRES,XGRES                                       TWOG1083
1501       ENDIF                                                             TWOG1084
1502  1003 FORMAT                                                            TWOG1085
1503      +  (1H ,10X,'I    Resonance mass               ',F12.5,            TWOG1086
1504      +  ' GeV',11X,'I',/,                                               TWOG1087
1505      +  1H ,10X,'I    Resonance total width        ',F12.5,             TWOG1088
1506      +  ' GeV',11X,'I')                                                 TWOG1089
1507 C
1508       IF (LWATE) THEN                                                   TWOG1090
1509       WRITE(LOUT,1001)                                                  TWOG1091
1510       WRITE(LOUT,1004)                                                  TWOG1092
1511       ELSE                                                              TWOG1095
1512       WRITE(LOUT,1001)                                                  TWOG1096
1513       WRITE(LOUT,1005)WGHTMX                                            TWOG1097
1514  1004 FORMAT(1H ,10X,'I    Produce weighted events      ',27X,'I')      TWOG1094
1515  1005 FORMAT(1H ,10X,'I    Maximum weight               ',F12.5,        TWOG1099
1516      +                                            '    ',11X,'I')       TWOG1100
1517       ENDIF                                                             TWOG1101
1518 C                                                                       TWOG1102
1519       WRITE(LOUT,1001)                                                  TWOG1103
1520       WRITE(LOUT,1000)                                                  TWOG1104
1521 C                                                                       TWOG1105
1522       DO 10 I=1,NCNTRS                                                  TWOG1106
1523    10 ICNTER(I)=0                                                       TWOG1107
1524       WGHTGT = 0.                                                       TWOG1108
1525       WGHTSM = 0.                                                       TWOG1109
1526       RETURN                                                            TWOG1110
1527       END                                                               TWOG1111
1528 C=========================================================================
1529 CDECK  ID>, DLDMX.  
1530       FUNCTION DLDMX(Rm,Ny)
1531       PARAMETER (isize=100000)
1532 C                                                      S.A.Sadovsky
1533 C   Calculate the Luminosity integral over dY            5.02.1995
1534 C             in the rapidity range ymin-ymax
1535 C             using luminosity interpolation table
1536 C      Ny/2 - number of subintervals for Simpson integration
1537 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1538       COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
1539      +                           xmin,xmax,amn,amx,nnms,alfun(isize)
1540       stpy=(ymx-ymn)/nny
1541       stpm=(amx-amn)/nnms                 ! Interpolation steps
1542 C
1543       im=(rm-amn)/stpm
1544       q =(rm-im*stpm-amn)/stpm
1545 C
1546       My = Ny/2
1547       Hy =(ymax-ymin)/My                  ! Integration step over Y
1548       ry = ymin
1549       iy =(ry-ymn)/stpy+1.
1550       p  =(ry-(iy-1)*stpy-ymn)/stpy       ! Four point interpolation:  1
1551       k  = iy+im*(nny+1)
1552       F1 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1553      *        q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1554 C
1555       SIM = 0.0
1556       DO 10 i=1,My
1557 C
1558         ry = ymin + Hy*i
1559         iy =(ry-ymn)/stpy+1.
1560         p  =(ry-(iy-1)*stpy-ymn)/stpy       ! Four point interpolation:  3
1561         k  = iy+im*(nny+1)
1562         F3 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1563      *        q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1564 C
1565         ry = ry   - Hy*0.5
1566         iy =(ry-ymn)/stpy+1.
1567         p  =(ry-(iy-1)*stpy-ymn)/stpy       ! Four point interpolation:  2
1568         k  = iy+im*(nny+1)
1569         F2 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1570      *        q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1571 C
1572         SIM = SIM + (F1 + 4.*F2 + F3)/6.
1573    10 F1  = F3
1574 C
1575       DLDMX = 1000.*SIM*HY                ! dL/dMx  in GeV^-1
1576       RETURN
1577       END
1578 *=======================================================================
1579 CDECK  ID>, LUMFUN. 
1580       SUBROUTINE LUMFUN(ifl,file,ymin,ymax,ny,xmn,xmx,nms)
1581 C-----------------------------------------------------------------------
1582 C preparation of luminosity function for generator
1583 C
1584 C     ifl=-1 - calculate function + save in file 'file'
1585 C     ifl= 0 - calculate function
1586 C     ifl= 1 - read function from file
1587 C       ymin,ymax,ny - rapidity range & number of points
1588 C       xmn,xmx,nms  - mass range & number of points
1589 C       Attention: dimensions of alfun are alfun(ny+1,nms+1)
1590 C
1591 C     Called by GGINIT
1592 C
1593 C     Author: G.V.Khaoustov 12.12.1994
1594 C     Modified: Yu.Kharlov 18-JUN-1997
1595 C-----------------------------------------------------------------------
1596       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1597      &               gvconst(4,10)
1598       REAL nc
1599       COMMON/D2LParam/lz,la,Rion,Gamma
1600       INTEGER lz,la
1601       REAL    Rion,Gamma
1602       PARAMETER (isize=100000)
1603       CHARACTER*(*) file
1604       REAL work(1000)
1605       COMMON/lumfunct/anorm,cexp,umn,umx,ymn,ymx,nny,
1606      +                           wmn,wmx,amn,amx,nnms,alfun(isize)
1607 C
1608       umn = ymin                    !  parameters for luminosity generation
1609       umx = ymax
1610       wmn = xmn
1611       wmx = xmx
1612 C
1613 C Define R : radii of the ions (in fm) (if set non-zero, use it)
1614       IF (Rion .LE. 0) THEN
1615         rion = 1.2*float(la)**(1.0/3.0)
1616       ENDIF
1617 *
1618 *       Fing max energy of gg-system (with safity factor 2)
1619 *
1620       egg_max = 2.*hbarc * gamma / rion
1621 *
1622       IF(ifl.EQ.1) THEN             !  read lum. function from file
1623         OPEN(unit=10,err=1,file=file,status='old',
1624      *                                        form='formatted')
1625         READ(10, *,err=3) anorm,cexp
1626         READ(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
1627         READ(10,100,err=3) (alfun(I),i=1,(nny+1)*(nnms+1))
1628         CLOSE(unit=10)
1629         PRINT *,'The Luminosity has been read from File '
1630         RETURN
1631       ENDIF
1632 c
1633       IF (ny*nms .GT. isize) GOTO 2
1634       PRINT *,'Start the Luminosity calculation'
1635       nny = ny
1636       ymn = ymin
1637       ymx = ymax
1638       nnms= nms
1639       amn = xmn
1640       amx = xmx
1641 c
1642       stepy=(ymx-ymn)/nny
1643       stepm=(amx-amn)/nnms
1644 *
1645 *         Lum. function calculation
1646 *
1647       DO IM = 0,nnms
1648         am=amn+im*stepm
1649 *         Find Y-limits for mass AM (Yu.Kharlov)
1650         ym_max  = alog (egg_max/am)
1651         ym_max  = max(ym_max,5.*stepy)
1652         ygg_max = min (ymax, ym_max)
1653         ygg_min = max (ymin,-ym_max)
1654         DO IY = 1,nny+1
1655           Y = ymn+(iy-1)*stepy
1656           IF (Y.LT.ygg_max .AND. Y.GT.ygg_min) THEN
1657             dl = D2LDMDY(1000.*am,Y)   ! Differential Luminosity in MeV^-1
1658           ELSE
1659             dl = 0.0
1660           END IF
1661           alfun(iy+(nny+1)*im) = DL
1662         ENDDO
1663       ENDDO
1664 c
1665       s=0.
1666       DO i=1,nnms+1   ! find max values for exp. parameters calcul.
1667         am=0.
1668         DO j=1,nny+1
1669           IF(am.LT.alfun(j+(nny+1)*(i-1))) am=alfun(j+(nny+1)*(i-1))
1670         ENDDO
1671         work(i)=am
1672         s=s+am
1673       ENDDO
1674       s=(s-work(nnms+1))*stepm
1675 c
1676       smax=0.      ! calculation optimized parameters for exp. generator
1677       iflag=-1
1678       DO 101 i=1,nnms
1679         DO 102 j=i+1,nnms+1
1680           IF (work(j) .LE. 0.0) GOTO 102
1681           c=alog(work(i)/work(j))/((i-j)*stepm)
1682           am=work(i)/exp(c*(i-1)*stepm)
1683           se=am/c*(exp(c*(amx-amn))-1.)
1684           IF(s/se.lt.smax) GOTO 102
1685           DO k=1,nnms+1
1686             amp=am*exp(c*(k-1)*stepm)
1687             IF((work(k)-amp)/amp.GT.0.00001) THEN  !  Accuracy problem
1688               GOTO 102
1689             ENDIF
1690           ENDDO
1691           smax =s/se
1692           anorm=am
1693           cexp =c
1694           iflag=1
1695   102   CONTINUE
1696   101 CONTINUE
1697       IF(iflag.LT.0) THEN
1698         PRINT *,'Fault to find best exponential approximation'
1699         STOP
1700       ENDIF
1701 c
1702       IF(ifl.EQ.0) RETURN        ! save lum. function
1703       OPEN(unit=10,err=1,file=file,status='unknown',
1704      *                                          form='formatted')
1705       WRITE(10, *,err=3) anorm,cexp
1706       WRITE(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
1707       WRITE(10,100,err=3) (alfun(i),i=1,(nny+1)*(nnms+1))
1708       CLOSE(unit=10)
1709       PRINT *,'The Luminosity calculation is finished'
1710       RETURN
1711     1 PRINT *,' file open error, file:',file
1712       STOP
1713     2 PRINT *,' array "ALFUN" to small to contain luminosity function'
1714       STOP
1715     3 PRINT *,' file read/writr error, file:',file
1716       STOP
1717   100 FORMAT(5e13.5)
1718    99 FORMAT(2e13.5,i5,2e13.5,i5)
1719       END
1720 C======================================================================
1721 CDECK  ID>, D2LDMDY.
1722       REAL FUNCTION  D2LDMDY(M,Y)
1723 C-----------------------------------------------------------------------
1724 C Calculate the Luminosity-function d^2L/dMdY
1725 C where M is the invariant mass and Y the rapidity
1726 C of the photon-photon-system
1727 C
1728 C Calculations are based on the formulae of
1729 C N. Baron, G. Baur, Phys. Rev. C
1730 C C. Cahn, J. D. Jackson, Nucl. Phys.
1731 C
1732 C Kai Hencken, 29. 9. 1994
1733 C
1734 C M : invariant mass (MeV)
1735 C Y : rapidity
1736 C
1737 C D2LDMDY : diff. Lum. (MeV**-1)
1738 C
1739 C     Called by LUMFUN
1740 C-----------------------------------------------------------------------
1741       IMPLICIT NONE
1742       COMMON/D2LParam/lz,la,Rion,Gamma
1743       INTEGER lz,la
1744       REAL    Rion,Gamma
1745       REAL M,Y
1746       REAL RM,W1,W2
1747       COMMON/D2LIParam/RM,W1,W2
1748
1749       REAL NClass,DeltaL
1750       EXTERNAL NClass,DeltaL
1751
1752 C
1753 C physical constants hbar*c and 1/alpha
1754 C
1755       REAL hbarc,falphai
1756       PARAMETER (hbarc=197.327053,falphai=137.0359895)
1757
1758       REAL LClass,DL
1759 C
1760 C calculate photon frequencies from M,Y
1761 C
1762       W1 = M/2.0 * EXP ( Y)
1763       W2 = M/2.0 * EXP (-Y)
1764 C
1765 C calculate Rho in MeV^-1 instead of fm
1766 C
1767       RM = Rion / hbarc
1768 C
1769 C calculate the "classical" value first.
1770 C
1771
1772       LClass = NClass (W1,RM,Gamma) * NClass (W2,RM,Gamma)
1773
1774 C
1775 C substract from this the Delta-value
1776 C
1777       DL = DeltaL ()
1778
1779       D2LDMDY = 2.0 / M * float(LZ)**4 / falphai**2 * (LClass - DL)
1780
1781       RETURN
1782       END
1783 *=======================================================================
1784 CDECK  ID>, NCLASS. 
1785       REAL FUNCTION NClass (W, R, Gamma)
1786 C
1787 C the classical photon-number without impact parameter cutoff:
1788 C
1789       IMPLICIT NONE
1790       REAL W,R, Gamma
1791
1792       REAL Pi
1793       PARAMETER (Pi = 3.141592)
1794
1795       REAL BesselK0,BesselK1
1796       EXTERNAL BesselK0,BesselK1
1797
1798       REAL Xi
1799
1800       Xi = W*R/Gamma
1801
1802       NClass = 2.0/PI * (Xi*BesselK0 (Xi)*BesselK1 (Xi) -
1803      &          Xi**2/2.0 *(BesselK1(Xi)**2 - BesselK0(Xi)**2))
1804
1805       RETURN
1806       END
1807 *=======================================================================
1808 CDECK  ID>, DELTAL. 
1809       REAL FUNCTION DeltaL ()
1810 C
1811 C the difference to the classical value
1812 C     Called by d2LdMdY
1813 C
1814       IMPLICIT NONE
1815       COMMON/D2LParam/lz,la,Rion,Gamma
1816       INTEGER lz,la
1817       REAL    Rion,Gamma
1818       REAL RM,W1,W2
1819       COMMON/D2LIParam/RM,W1,W2
1820
1821       DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1822       COMMON /GAUSSD/XGAUSS,WGAUSS
1823
1824       INTEGER N,I
1825       REAL    b1
1826       REAL*8  t,tmin,tmax
1827       REAL*8  Sum,Int,Int2
1828
1829       REAL Itgb1
1830       EXTERNAL Itgb1
1831
1832       REAL   Pi
1833       REAL*8 Accur
1834       PARAMETER (Pi = 3.141592,Accur = 1D -2)
1835
1836 C
1837 C integrate first over b1
1838 C
1839
1840 C
1841 C Loop incrementing the boundary
1842 C
1843       tmin = 0.00
1844       tmax = 0.25
1845       Sum  = 0.00
1846
1847   100 CONTINUE
1848 C
1849 C Loop for the Gauss integration
1850 C
1851       Int=0.0
1852       DO N=1,6
1853         Int2 = Int
1854         Int=0.0
1855         DO I=2**N-1,2**(N+1)-2
1856           t   = (tmax-tmin)/2.0*XGAUSS(I) + (tmax+tmin)/2.0
1857           b1  = RM * DEXP (t)
1858           Int = Int + WGAUSS(I) * Itgb1(b1) * b1**2
1859         ENDDO
1860         Int = (tmax-tmin)/2.0*Int
1861         IF (Int .NE. 0.) THEN
1862           IF (DABS ((Int2-Int)/Int) .LT. Accur) GOTO 200
1863         END IF
1864       ENDDO
1865       IF (int .LT. 1.D-20) THEN
1866         int = 0.0D0
1867       ELSE
1868         WRITE(*,*) ' (b1) GAUSS MAY BE INACCURATE'
1869       END IF
1870
1871   200 CONTINUE
1872       Sum = Sum + Int
1873       IF (sum .NE. 0.0) THEN
1874         IF (ABS (Int2/Sum) .LT. Accur) GOTO 300
1875       ELSE
1876         IF (b1 .GT. 1.E+06)            GOTO 300
1877       END IF
1878       tmin = tmax
1879       tmax = tmax + 0.5
1880       GOTO 100
1881
1882   300 CONTINUE
1883
1884       DeltaL = 4.0*Pi * Sum
1885       RETURN
1886       END
1887
1888 CDECK  ID>, ITGB1.  
1889       REAL FUNCTION Itgb1(b)
1890 C
1891 C integration then over b2
1892 C     Called by DeltaL
1893 C
1894 C                                  Corrected by S.A.Sadovsky
1895 C                                             25.10.1994
1896       IMPLICIT NONE
1897 c  +seq,d2lpar. Not used Yu.Kh. 18-JUN-1997
1898       REAL     b
1899       REAL*8   b1,bmin,bmax
1900       REAL     RM,W1,W2
1901       COMMON/D2LIParam/RM,W1,W2
1902       DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1903       COMMON /GAUSSD/XGAUSS,WGAUSS
1904       INTEGER  N,I
1905       REAL*8   b2,Int,Int2
1906       REAL*8   Itgb2
1907       EXTERNAL Itgb2
1908       REAL*8   Accur
1909       PARAMETER (Accur = 1.D-2)
1910
1911       b1 = b
1912       bmin = b1 - 2.0*RM
1913       IF (RM .GT. bmin) THEN
1914         bmin = RM
1915       ENDIF
1916       bmax = b1 + 2.0 * RM
1917       Int = 0.0
1918       DO N=1,6
1919         Int2 = Int
1920         Int = 0.0
1921         DO I=2**N-1,2**(N+1)-2
1922           b2 = (bmax-bmin)/2.0*XGAUSS(I) + (bmax+bmin)/2.0
1923           Int = Int + WGAUSS(I) * Itgb2(b1,b2) * b2
1924         END DO
1925         Int = (bmax-bmin)/2.0*Int
1926         IF (Int.NE.0.) THEN
1927           IF (DABS((Int2 - Int)/Int) .LT. Accur) GOTO 100
1928         ENDIF
1929       ENDDO
1930       IF (int .LT. 1.D-20) THEN
1931         int = 0.0D0
1932       ELSE
1933         WRITE(*,*) ' (b2) GAUSS MAY BE INACCURATE, Itgb1=',Int
1934       END IF
1935   100 CONTINUE
1936       Itgb1 = Int
1937       RETURN
1938       END
1939
1940 CDECK  ID>, ITGB2.  
1941       REAL*8 FUNCTION Itgb2 (b1,b2)
1942 C
1943 C       The function to be integrated over
1944 C       Called by Itgb1
1945 C
1946 C                                  Corrected by S.A.Sadovsky
1947 C                                             25.10.1994
1948       IMPLICIT NONE
1949       COMMON/D2LParam/lz,la,Rion,Gamma
1950       INTEGER lz,la
1951       REAL    Rion,Gamma
1952       REAL*8 b1,b2,arg
1953       REAL RM,W1,W2
1954       COMMON/D2LIParam/RM,W1,W2
1955
1956       REAL Nphoton
1957       EXTERNAL Nphoton
1958
1959       arg = (b1**2 + b2**2 - 4.0 * RM**2) / (2.0*b1*b2)
1960       IF (arg .GT. 1.D0) arg = 1.D0
1961       Itgb2 = Nphoton(W1,sngl(b1),Gamma) *
1962      &        Nphoton(W2,sngl(b2),Gamma) *
1963      &        DACOS (arg)
1964
1965       RETURN
1966       END
1967
1968 CDECK  ID>, NPHOTON.
1969       REAL FUNCTION Nphoton (W,Rho,Gamma)
1970 C
1971 C       The differential photonnumber for a nucleus
1972 C       (without form factor)
1973 C
1974       IMPLICIT NONE
1975       REAL W,Rho,Gamma
1976
1977       REAL BESSELK1
1978       EXTERNAL BESSELK1
1979
1980       REAL Wphib,WGamma
1981
1982       REAL PI
1983       PARAMETER (PI=3.141592)
1984
1985       WGamma = W/Gamma
1986 C
1987 C without form factor
1988 C
1989       Wphib = WGamma*BESSELK1 (WGamma*Rho)
1990
1991       Nphoton = 1.0/PI**2 * Wphib**2
1992
1993       RETURN
1994       END
1995
1996 CDECK  ID>, BESSELK0.   
1997       REAL FUNCTION BESSELK0(X)
1998 C
1999 C  Special functions I0,K0,I1,K1
2000 C  see Abramowitz&Stegun
2001 C
2002       IMPLICIT NONE
2003       REAL X,Y
2004
2005       REAL BESSELI0
2006       EXTERNAL BESSELI0
2007
2008       IF (X .LT. 2.0) THEN
2009
2010         Y = X**2/4.0
2011
2012         BESSELK0 =
2013      &    (-LOG(X/2.0)*BESSELI0(X))+(-.57721566+Y*(0.42278420
2014      &    +Y*(0.23069756+Y*(0.3488590E-1+Y*(0.262698E-2
2015      &    +Y*(0.10750E-3+Y*0.740E-5))))))
2016
2017       ELSE
2018
2019         Y = 2.0/X
2020
2021         BESSELK0 =
2022      &    (EXP(-X)/SQRT(X))*(1.25331414+Y*(-0.7832358E-1
2023      &    +Y*(0.2189568E-1+Y*(-0.1062446E-1+Y*(0.587872E-2
2024      &    +Y*(-0.251540E-2+Y*0.53208E-3))))))
2025
2026       ENDIF
2027
2028       RETURN
2029       END
2030
2031
2032 CDECK  ID>, BESSELI0.   
2033       REAL FUNCTION BESSELI0(X)
2034       IMPLICIT NONE
2035
2036       REAL X,Y,AX
2037
2038       AX = ABS(X)
2039
2040       IF (AX .LT. 3.75) THEN
2041
2042         Y = (X/3.75)**2
2043
2044         BESSELI0 =
2045      &    1.0+Y*(3.5156229+Y*(3.0899424+Y*(1.2067492
2046      &    +Y*(0.2659732+Y*(0.360768E-1+Y*0.45813E-2)))))
2047
2048       ELSE
2049
2050         Y = 3.75/AX
2051
2052         BESSELI0 =
2053      &    (EXP(AX)/SQRT(AX))*(0.39894228+Y*(0.1328592E-1
2054      &    +Y*(0.225319E-2+Y*(-0.157565E-2+Y*(0.916281E-2
2055      &    +Y*(-0.2057706E-1+Y*(0.2635537E-1+Y*(-0.1647633E-1
2056      &    +Y*0.392377E-2))))))))
2057
2058       ENDIF
2059
2060       RETURN
2061       END
2062
2063
2064 CDECK  ID>, BESSELK1.   
2065       REAL FUNCTION BESSELK1(X)
2066       IMPLICIT NONE
2067
2068       REAL X,Y
2069
2070       REAL BESSELI1
2071       EXTERNAL BESSELI1
2072
2073       IF (X .LT. 2.0) THEN
2074
2075         Y = X**2/4.0
2076         BESSELK1 =
2077      &    (LOG(X/2.0)*BESSELI1(X))+(1.0/X)*(1.0+Y*(0.15443144
2078      &    +Y*(-0.67278579+Y*(-0.18156897+Y*(-0.1919402E-1
2079      &    +Y*(-0.110404E-2+Y*(-0.4686E-4)))))))
2080
2081       ELSE
2082
2083         Y=2.0/X
2084         BESSELK1 =
2085      &    (EXP(-X)/SQRT(X))*(1.25331414+Y*(0.23498619
2086      &    +Y*(-0.3655620E-1+Y*(0.1504268E-1+Y*(-0.780353E-2
2087      &    +Y*(0.325614E-2+Y*(-0.68245E-3)))))))
2088
2089       ENDIF
2090
2091       RETURN
2092       END
2093
2094 CDECK  ID>, BESSELI1.   
2095       REAL FUNCTION BESSELI1(X)
2096       IMPLICIT NONE
2097
2098       REAL X,Y,AX
2099
2100       AX = ABS(X)
2101
2102       IF (AX .LT. 3.75) THEN
2103
2104         Y = (X/3.75)**2
2105
2106         BESSELI1 =
2107      &    AX*(0.5+Y*(0.87890594+Y*(0.51498869+Y*(0.15084934
2108      &    +Y*(0.2658733E-1+Y*(0.301532E-2+Y*0.32411E-3))))))
2109
2110
2111       ELSE
2112
2113         Y = 3.75/AX
2114
2115         BESSELI1 =
2116      &    0.2282967E-1+Y*(-0.2895312E-1+Y*(0.1787654E-1
2117      &    -Y*0.420059E-2))
2118
2119         BESSELI1 =
2120      &    0.39894228+Y*(-0.3988024E-1+Y*(-0.362018E-2
2121      &    +Y*(0.163801E-2+Y*(-0.1031555E-1+Y*BESSELI1))))
2122
2123         BESSELI1 = BESSELI1 *
2124      &    EXP(AX)/SQRT(AX)
2125
2126       ENDIF
2127
2128       IF (X .LT. 0) THEN
2129         BESSELI1 = -BESSELI1
2130       ENDIF
2131
2132       RETURN
2133       END
2134
2135 CDECK  ID>, GAUSSDAT.   
2136 C DATA for the Gauss integration, x-values and weight for a
2137 C N=2,4,8,16,32,64 point Gauss integration.
2138 C
2139 C based on program GAUSS64 of D. Trautmann
2140 C data from Abramowitz & Stegun
2141 C
2142 C KAI HENCKEN, FEBRUAR 1993
2143 C
2144       BLOCKDATA GAUSSDATA
2145       IMPLICIT NONE
2146       DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
2147       COMMON /GAUSSD/XGAUSS,WGAUSS
2148
2149       DATA XGAUSS(1)/ .57735026918962576D0/
2150       DATA XGAUSS(2)/-.57735026918962576D0/
2151       DATA WGAUSS(1)/ 1.00000000000000000D0/
2152       DATA WGAUSS(2)/ 1.00000000000000000D0/
2153
2154       DATA XGAUSS(3)/ .33998104358485627D0/
2155       DATA XGAUSS(4)/ .86113631159405258D0/
2156       DATA XGAUSS(5)/-.33998104358485627D0/
2157       DATA XGAUSS(6)/-.86113631159405258D0/
2158       DATA WGAUSS(3)/ .65214515486254613D0/
2159       DATA WGAUSS(4)/ .34785484513745385D0/
2160       DATA WGAUSS(5)/ .65214515486254613D0/
2161       DATA WGAUSS(6)/ .34785484513745385D0/
2162
2163       DATA XGAUSS(7)/ .18343464249564981D0/
2164       DATA XGAUSS(8)/ .52553240991632899D0/
2165       DATA XGAUSS(9)/ .79666647741362674D0/
2166       DATA XGAUSS(10)/ .96028985649753623D0/
2167       DATA XGAUSS(11)/-.18343464249564981D0/
2168       DATA XGAUSS(12)/-.52553240991632899D0/
2169       DATA XGAUSS(13)/-.79666647741362674D0/
2170       DATA XGAUSS(14)/-.96028985649753623D0/
2171       DATA WGAUSS(7)/ .36268378337836198D0/
2172       DATA WGAUSS(8)/ .31370664587788727D0/
2173       DATA WGAUSS(9)/ .22238103445337448D0/
2174       DATA WGAUSS(10)/ .10122853629037627D0/
2175       DATA WGAUSS(11)/ .36268378337836198D0/
2176       DATA WGAUSS(12)/ .31370664587788727D0/
2177       DATA WGAUSS(13)/ .22238103445337448D0/
2178       DATA WGAUSS(14)/ .10122853629037627D0/
2179
2180       DATA XGAUSS(15)/ .0950125098376374402D0/
2181       DATA XGAUSS(16)/ .281603550779258913D0/
2182       DATA XGAUSS(17)/ .458016777657227386D0/
2183       DATA XGAUSS(18)/ .617876244402643748D0/
2184       DATA XGAUSS(19)/ .755404408355003034D0/
2185       DATA XGAUSS(20)/ .865631202387831744D0/
2186       DATA XGAUSS(21)/ .944575023073232576D0/
2187       DATA XGAUSS(22)/ .989400934991649933D0/
2188       DATA XGAUSS(23)/-.0950125098376374402D0/
2189       DATA XGAUSS(24)/-.281603550779258913D0/
2190       DATA XGAUSS(25)/-.458016777657227386D0/
2191       DATA XGAUSS(26)/-.617876244402643748D0/
2192       DATA XGAUSS(27)/-.755404408355003034D0/
2193       DATA XGAUSS(28)/-.865631202387831744D0/
2194       DATA XGAUSS(29)/-.944575023073232576D0/
2195       DATA XGAUSS(30)/-.989400934991649933D0/
2196       DATA WGAUSS(15)/ .189450610455068496D0/
2197       DATA WGAUSS(16)/ .182603415044923589D0/
2198       DATA WGAUSS(17)/ .169156519395002538D0/
2199       DATA WGAUSS(18)/ .149595988816576732D0/
2200       DATA WGAUSS(19)/ .124628971255533872D0/
2201       DATA WGAUSS(20)/ .0951585116824927848D0/
2202       DATA WGAUSS(21)/ .0622535239386478929D0/
2203       DATA WGAUSS(22)/ .0271524594117540949D0/
2204       DATA WGAUSS(23)/ .189450610455068496D0/
2205       DATA WGAUSS(24)/ .182603415044923589D0/
2206       DATA WGAUSS(25)/ .169156519395002538D0/
2207       DATA WGAUSS(26)/ .149595988816576732D0/
2208       DATA WGAUSS(27)/ .124628971255533872D0/
2209       DATA WGAUSS(28)/ .0951585116824927848D0/
2210       DATA WGAUSS(29)/ .0622535239386478929D0/
2211       DATA WGAUSS(30)/ .0271524594117540949D0/
2212
2213       DATA XGAUSS(31)/ .0483076656877383162D0/
2214       DATA XGAUSS(32)/ .144471961582796493D0/
2215       DATA XGAUSS(33)/ .239287362252137075D0/
2216       DATA XGAUSS(34)/ .331868602282127650D0/
2217       DATA XGAUSS(35)/ .421351276130635345D0/
2218       DATA XGAUSS(36)/ .506899908932229390D0/
2219       DATA XGAUSS(37)/ .587715757240762329D0/
2220       DATA XGAUSS(38)/ .663044266930215201D0/
2221       DATA XGAUSS(39)/ .732182118740289680D0/
2222       DATA XGAUSS(40)/ .794483795967942407D0/
2223       DATA XGAUSS(41)/ .849367613732569970D0/
2224       DATA XGAUSS(42)/ .896321155766052124D0/
2225       DATA XGAUSS(43)/ .934906075937739689D0/
2226       DATA XGAUSS(44)/ .964762255587506430D0/
2227       DATA XGAUSS(45)/ .985611511545268335D0/
2228       DATA XGAUSS(46)/ .997263861849481564D0/
2229       DATA XGAUSS(47)/-.0483076656877383162D0/
2230       DATA XGAUSS(48)/-.144471961582796493D0/
2231       DATA XGAUSS(49)/-.239287362252137075D0/
2232       DATA XGAUSS(50)/-.331868602282127650D0/
2233       DATA XGAUSS(51)/-.421351276130635345D0/
2234       DATA XGAUSS(52)/-.506899908932229390D0/
2235       DATA XGAUSS(53)/-.587715757240762329D0/
2236       DATA XGAUSS(54)/-.663044266930215201D0/
2237       DATA XGAUSS(55)/-.732182118740289680D0/
2238       DATA XGAUSS(56)/-.794483795967942407D0/
2239       DATA XGAUSS(57)/-.849367613732569970D0/
2240       DATA XGAUSS(58)/-.896321155766052124D0/
2241       DATA XGAUSS(59)/-.934906075937739689D0/
2242       DATA XGAUSS(60)/-.964762255587506430D0/
2243       DATA XGAUSS(61)/-.985611511545268335D0/
2244       DATA XGAUSS(62)/-.997263861849481564D0/
2245       DATA WGAUSS(31)/ .0965400885147278006D0/
2246       DATA WGAUSS(32)/ .0956387200792748594D0/
2247       DATA WGAUSS(33)/ .0938443990808045654D0/
2248       DATA WGAUSS(34)/ .0911738786957638847D0/
2249       DATA WGAUSS(35)/ .0876520930044038111D0/
2250       DATA WGAUSS(36)/ .0833119242269467552D0/
2251       DATA WGAUSS(37)/ .0781938957870703065D0/
2252       DATA WGAUSS(38)/ .0723457941088485062D0/
2253       DATA WGAUSS(39)/ .0658222227763618468D0/
2254       DATA WGAUSS(40)/ .0586840934785355471D0/
2255       DATA WGAUSS(41)/ .0509980592623761762D0/
2256       DATA WGAUSS(42)/ .0428358980222266807D0/
2257       DATA WGAUSS(43)/ .0342738629130214331D0/
2258       DATA WGAUSS(44)/ .0253920653092620595D0/
2259       DATA WGAUSS(45)/ .0162743947309056706D0/
2260       DATA WGAUSS(46)/ .00701861000947009660D0/
2261       DATA WGAUSS(47)/ .0965400885147278006D0/
2262       DATA WGAUSS(48)/ .0956387200792748594D0/
2263       DATA WGAUSS(49)/ .0938443990808045654D0/
2264       DATA WGAUSS(50)/ .0911738786957638847D0/
2265       DATA WGAUSS(51)/ .0876520930044038111D0/
2266       DATA WGAUSS(52)/ .0833119242269467552D0/
2267       DATA WGAUSS(53)/ .0781938957870703065D0/
2268       DATA WGAUSS(54)/ .0723457941088485062D0/
2269       DATA WGAUSS(55)/ .0658222227763618468D0/
2270       DATA WGAUSS(56)/ .0586840934785355471D0/
2271       DATA WGAUSS(57)/ .0509980592623761762D0/
2272       DATA WGAUSS(58)/ .0428358980222266807D0/
2273       DATA WGAUSS(59)/ .0342738629130214331D0/
2274       DATA WGAUSS(60)/ .0253920653092620595D0/
2275       DATA WGAUSS(61)/ .0162743947309056706D0/
2276       DATA WGAUSS(62)/ .00701861000947009660D0/
2277
2278       DATA XGAUSS(63)/ .02435029266342443250D0/
2279       DATA XGAUSS(64)/ .0729931217877990394D0/
2280       DATA XGAUSS(65)/ .121462819296120554D0/
2281       DATA XGAUSS(66)/ .169644420423992818D0/
2282       DATA XGAUSS(67)/ .217423643740007084D0/
2283       DATA XGAUSS(68)/ .264687162208767416D0/
2284       DATA XGAUSS(69)/ .311322871990210956D0/
2285       DATA XGAUSS(70)/ .357220158337668116D0/
2286       DATA XGAUSS(71)/ .402270157963991604D0/
2287       DATA XGAUSS(72)/ .446366017253464088D0/
2288       DATA XGAUSS(73)/ .489403145707052957D0/
2289       DATA XGAUSS(74)/ .531279464019894546D0/
2290       DATA XGAUSS(75)/ .571895646202634034D0/
2291       DATA XGAUSS(76)/ .611155355172393250D0/
2292       DATA XGAUSS(77)/ .648965471254657340D0/
2293       DATA XGAUSS(78)/ .685236313054233243D0/
2294       DATA XGAUSS(79)/ .719881850171610827D0/
2295       DATA XGAUSS(80)/ .752819907260531897D0/
2296       DATA XGAUSS(81)/ .783972358943341408D0/
2297       DATA XGAUSS(82)/ .813265315122797560D0/
2298       DATA XGAUSS(83)/ .840629296252580363D0/
2299       DATA XGAUSS(84)/ .865999398154092820D0/
2300       DATA XGAUSS(85)/ .889315445995114106D0/
2301       DATA XGAUSS(86)/ .910522137078502806D0/
2302       DATA XGAUSS(87)/ .929569172131939576D0/
2303       DATA XGAUSS(88)/ .946411374858402816D0/
2304       DATA XGAUSS(89)/ .961008799652053719D0/
2305       DATA XGAUSS(90)/ .973326827789910964D0/
2306       DATA XGAUSS(91)/ .983336253884625957D0/
2307       DATA XGAUSS(92)/ .991013371476744321D0/
2308       DATA XGAUSS(93)/ .996340116771955279D0/
2309       DATA XGAUSS(94)/ .999305041735772139D0/
2310       DATA XGAUSS(95)/-.02435029266342443250D0/
2311       DATA XGAUSS(96)/-.0729931217877990394D0/
2312       DATA XGAUSS(97)/-.121462819296120554D0/
2313       DATA XGAUSS(98)/-.169644420423992818D0/
2314       DATA XGAUSS(99)/-.217423643740007084D0/
2315       DATA XGAUSS(100)/-.264687162208767416D0/
2316       DATA XGAUSS(101)/-.311322871990210956D0/
2317       DATA XGAUSS(102)/-.357220158337668116D0/
2318       DATA XGAUSS(103)/-.402270157963991604D0/
2319       DATA XGAUSS(104)/-.446366017253464088D0/
2320       DATA XGAUSS(105)/-.489403145707052957D0/
2321       DATA XGAUSS(106)/-.531279464019894546D0/
2322       DATA XGAUSS(107)/-.571895646202634034D0/
2323       DATA XGAUSS(108)/-.611155355172393250D0/
2324       DATA XGAUSS(109)/-.648965471254657340D0/
2325       DATA XGAUSS(110)/-.685236313054233243D0/
2326       DATA XGAUSS(111)/-.719881850171610827D0/
2327       DATA XGAUSS(112)/-.752819907260531897D0/
2328       DATA XGAUSS(113)/-.783972358943341408D0/
2329       DATA XGAUSS(114)/-.813265315122797560D0/
2330       DATA XGAUSS(115)/-.840629296252580363D0/
2331       DATA XGAUSS(116)/-.865999398154092820D0/
2332       DATA XGAUSS(117)/-.889315445995114106D0/
2333       DATA XGAUSS(118)/-.910522137078502806D0/
2334       DATA XGAUSS(119)/-.929569172131939576D0/
2335       DATA XGAUSS(120)/-.946411374858402816D0/
2336       DATA XGAUSS(121)/-.961008799652053719D0/
2337       DATA XGAUSS(122)/-.973326827789910964D0/
2338       DATA XGAUSS(123)/-.983336253884625957D0/
2339       DATA XGAUSS(124)/-.991013371476744321D0/
2340       DATA XGAUSS(125)/-.996340116771955279D0/
2341       DATA XGAUSS(126)/-.999305041735772139D0/
2342       DATA WGAUSS(63)/ .0486909570091397204D0/
2343       DATA WGAUSS(64)/ .0485754674415034269D0/
2344       DATA WGAUSS(65)/ .0483447622348029572D0/
2345       DATA WGAUSS(66)/ .0479993885964583077D0/
2346       DATA WGAUSS(67)/ .0475401657148303087D0/
2347       DATA WGAUSS(68)/ .0469681828162100173D0/
2348       DATA WGAUSS(69)/ .0462847965813144172D0/
2349       DATA WGAUSS(70)/ .0454916279274181445D0/
2350       DATA WGAUSS(71)/ .0445905581637565631D0/
2351       DATA WGAUSS(72)/ .0435837245293234534D0/
2352       DATA WGAUSS(73)/ .0424735151236535890D0/
2353       DATA WGAUSS(74)/ .0412625632426235286D0/
2354       DATA WGAUSS(75)/ .0399537411327203414D0/
2355       DATA WGAUSS(76)/ .0385501531786156291D0/
2356       DATA WGAUSS(77)/ .0370551285402400460D0/
2357       DATA WGAUSS(78)/ .0354722132568823838D0/
2358       DATA WGAUSS(79)/ .0338051618371416094D0/
2359       DATA WGAUSS(80)/ .0320579283548515535D0/
2360       DATA WGAUSS(81)/ .0302346570724024789D0/
2361       DATA WGAUSS(82)/ .0283396726142594832D0/
2362       DATA WGAUSS(83)/ .0263774697150546587D0/
2363       DATA WGAUSS(84)/ .0243527025687108733D0/
2364       DATA WGAUSS(85)/ .0222701738083832542D0/
2365       DATA WGAUSS(86)/ .0201348231535302094D0/
2366       DATA WGAUSS(87)/ .0179517157756973431D0/
2367       DATA WGAUSS(88)/ .0157260304760247193D0/
2368       DATA WGAUSS(89)/ .0134630478967186426D0/
2369       DATA WGAUSS(90)/ .0111681394601311288D0/
2370       DATA WGAUSS(91)/ .00884675982636394772D0/
2371       DATA WGAUSS(92)/ .00650445796897836286D0/
2372       DATA WGAUSS(93)/ .00414703326056246764D0/
2373       DATA WGAUSS(94)/ .00178328072169643295D0/
2374       DATA WGAUSS(95)/ .0486909570091397204D0/
2375       DATA WGAUSS(96)/ .0485754674415034269D0/
2376       DATA WGAUSS(97)/ .0483447622348029572D0/
2377       DATA WGAUSS(98)/ .0479993885964583077D0/
2378       DATA WGAUSS(99)/ .0475401657148303087D0/
2379       DATA WGAUSS(100)/ .0469681828162100173D0/
2380       DATA WGAUSS(101)/ .0462847965813144172D0/
2381       DATA WGAUSS(102)/ .0454916279274181445D0/
2382       DATA WGAUSS(103)/ .0445905581637565631D0/
2383       DATA WGAUSS(104)/ .0435837245293234534D0/
2384       DATA WGAUSS(105)/ .0424735151236535890D0/
2385       DATA WGAUSS(106)/ .0412625632426235286D0/
2386       DATA WGAUSS(107)/ .0399537411327203414D0/
2387       DATA WGAUSS(108)/ .0385501531786156291D0/
2388       DATA WGAUSS(109)/ .0370551285402400460D0/
2389       DATA WGAUSS(110)/ .0354722132568823838D0/
2390       DATA WGAUSS(111)/ .0338051618371416094D0/
2391       DATA WGAUSS(112)/ .0320579283548515535D0/
2392       DATA WGAUSS(113)/ .0302346570724024789D0/
2393       DATA WGAUSS(114)/ .0283396726142594832D0/
2394       DATA WGAUSS(115)/ .0263774697150546587D0/
2395       DATA WGAUSS(116)/ .0243527025687108733D0/
2396       DATA WGAUSS(117)/ .0222701738083832542D0/
2397       DATA WGAUSS(118)/ .0201348231535302094D0/
2398       DATA WGAUSS(119)/ .0179517157756973431D0/
2399       DATA WGAUSS(120)/ .0157260304760247193D0/
2400       DATA WGAUSS(121)/ .0134630478967186426D0/
2401       DATA WGAUSS(122)/ .0111681394601311288D0/
2402       DATA WGAUSS(123)/ .00884675982636394772D0/
2403       DATA WGAUSS(124)/ .00650445796897836286D0/
2404       DATA WGAUSS(125)/ .00414703326056246764D0/
2405       DATA WGAUSS(126)/ .00178328072169643295D0/
2406
2407       END
2408 *=======================================================================
2409       SUBROUTINE LORENB (U,PS,PI,PF)
2410 C
2411 C CERN PROGLIB# U102    LORENB          .VERSION KERNFOR  4.04  821124
2412 C ORIG. 20/08/75 L.PAPE
2413 C
2414       DOUBLE PRECISION PF4, FN
2415       DIMENSION      PS(4),PI(4),PF(4)
2416
2417       IF (PS(4).EQ.U) GO TO 17
2418       PF4  = (PI(4)*PS(4)+PI(3)*PS(3)+PI(2)*PS(2)+PI(1)*PS(1)) / U
2419       FN   = (PF4+PI(4)) / (PS(4)+U)
2420       PF(1)= PI(1) + FN*PS(1)
2421       PF(2)= PI(2) + FN*PS(2)
2422       PF(3)= PI(3) + FN*PS(3)
2423       PF(4)= PF4
2424       GO TO 18
2425 C
2426    17 PF(1)= PI(1)
2427       PF(2)= PI(2)
2428       PF(3)= PI(3)
2429       PF(4)= PI(4)
2430 C
2431    18 CONTINUE
2432 C
2433       RETURN
2434 C
2435       END