CMake: prepend_prefix() to every elm in a list
[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)*(1.+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,'(/1x,79(''=''))')
602       WRITE (6,'(1x,''I'',36x,''TPHIC'',36x,      ''I'')')
603       WRITE (6,'(1x,''I'',34x,''Process '',i1,34x,''I'')') iproc
604       WRITE (6,'(1x,''I'',34x,''---------'',  34x,''I'')')
605       WRITE (6,'(1x,''I'',77x,                    ''I'')')
606       WRITE (6,'(1x,''I'',20x,''Events thrown   : '',i8,31x,''I'')')
607      &       ntry
608       WRITE (6,'(1x,''I'',20x,''Events accepted : '',i8,31x,''I'')')
609      &       nevent
610       WRITE (6,'(1x,''I'',20x,''Cross section   : '',e10.3,'' +- '','//
611      &      'e10.3,'' nb'',12x ''I'')') xstot, xstote
612       WRITE (6,'(1x,''I'',77x,                    ''I'')')
613       WRITE (6,'(1x,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-amprt(2))*(eb+amprt(2)))
830         pc = pb
831         costhe = 2*ggrnd(0) - 1.
832         sinthe = sqrt((1.-costhe)*(1.+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)*(1.+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-1.)*(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)*(1.+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-am4)*(e4+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-amas)*(eb+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       save ifl,aky,akm,stpy,stpm,rmin,rmax,dr
1269 C
1270       IF (ifl.EQ.0) THEN
1271         IF (ymin.LT.ymn .OR. ymin.GT.ymx .OR.
1272      *      ymax.LT.ymn .OR. ymax.GT.ymx .OR.
1273      *      xmin.LT.amn .OR. xmin.GT.amx .OR.
1274      *      xmax.LT.amn .OR. xmax.GT.amx) THEN
1275           WRITE (6,*) 'Luminosity function is not defined in this range'
1276           STOP
1277         ENDIF
1278         aky = ymax-ymin
1279         akm = xmax-xmin
1280         stpy = (ymx-ymn)/nny
1281         stpm = (amx-amn)/nnms
1282         rmin = exp(cexp*(xmax))
1283         rmax = exp(cexp*(xmin))
1284         dr = rmax-rmin
1285         ifl = 1
1286       ENDIF
1287 C
1288     1 CONTINUE
1289       IF (.NOT.lres) THEN
1290         rm  = ggrnd(0)                     ! random mass
1291         rm  = rmin + rm*dr
1292         rm  = alog(rm)/cexp                ! exponetial law
1293       ELSE
1294         wsq = xmres**2 + xmres*xgres * tan(pi*(ggrnd(0)-0.5))
1295         IF (wsq .LT. xmin*xmin) GOTO 1
1296         IF (wsq .GT. xmax*xmax) GOTO 1
1297         rm  = sqrt (wsq)
1298       END IF
1299 C
1300       ry = ggrnd(0)                        ! random Y
1301 *
1302 *       Find y-limits for generated gg-mass: y_max = log(E_max/M_gg)
1303 *
1304       ym_max  = alog (egg_max/rm)
1305       ygg_max = min (ymax, ym_max)
1306       ygg_min = max (ymin,-ym_max)
1307       ry = ygg_min + ry*(ygg_max-ygg_min)
1308       iy = (ry-ymn)/stpy+1.
1309 C
1310       im = (rm-amn)/stpm
1311       r = anorm*exp(cexp*rm)*ggrnd(0)      ! normalization
1312       p = (ry-(iy-1)*stpy-ymn)/stpy        ! four point interpolation
1313       q = (rm-im*stpm-amn)/stpm
1314       k = iy+im*(nny+1)
1315       f = (1.-p)*(1.-q)*alfun(k)       + p*(1.-q)*alfun(k+1) +
1316      *         q*(1.-p)*alfun(k+nny+1) +      q*p*alfun(k+2+nny)
1317       IF (r .GT. f) GOTO 1
1318       wlum =  1000.*f                      ! differential luminosity in GeV^-1
1319 *
1320       RETURN
1321       END
1322 *
1323 *=======================================================================
1324 CDECK  ID>, GGDATA. 
1325       BLOCK DATA ggdata
1326 *-----------------------------------------------------------------------
1327 *       Some useful constants and masses of MSSM particles
1328 *       Author: Yu.Kharlov
1329 *-----------------------------------------------------------------------
1330       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1331      &               gvconst(4,10)
1332       REAL nc
1333       COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
1334      &               amin,amax,ymin,ymax,nmas,ny, kferm,
1335      &               kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
1336      &               thetamin, costhv1, kv1,kv2,gvpar(4)
1337       CHARACTER lumfil*80
1338       COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
1339      &  xstote, ssbr(10)
1340       COMMON /ggmssm/ xm1,   xm2,    xmg,    xms,    xmtl,   xmtr,
1341      &                xmll,  xmlr,   xmnl,   xtanb,  xmha,   xmu,
1342      &                xmt,   xat,    xmbr,   xab,    u11,    v11
1343       COMMON/D2LParam/lz,la,Rion,Gamma
1344       INTEGER lz,la
1345       REAL    Rion,Gamma
1346       PARAMETER (al = 1./128.)
1347       DATA pi /3.14159276/, alpha /al/, gev2nb /389379./,
1348      &     hbarc /0.197327/
1349       DATA iproc, nevent, ilumf /1, 10000, -1/
1350       DATA lumfil /'gglum.dat'/
1351       DATA Rion /0.0/
1352       DATA amprt /110.,10.,150.,-1.,0.106/
1353       DATA thetamin /1.e-03/
1354       DATA costhv1 /1.0/
1355       DATA xsbra /1./
1356       DATA ny, nmas /50, 50/
1357       DATA xm1,   xm2,    xmg,    xms,    xmtl,   xmtr,
1358      &     xmll,  xmlr,   xmnl,   xtanb,  xmha,   xmu,
1359      &     xmt,   xat,    xmbr,   xab
1360      &     /30.,   50.,   250.,   700.,   600.,   600.,
1361      &     400.,  400.,   400.,   4.,     300.,  -300.,
1362      &     174.,  300.,   300.,   300./
1363       DATA gvconst / 63.0, 0.22, 6.0, 1.0,    ! rho rho
1364      &               14.0, 0.22, 6.0, 1.0,    ! rho omega
1365      &                7.0, 0.22, 4.0, 1.0,    ! rho phi
1366      &               0.05, 0.80, 2.5, 0.0,    ! rho J/psi
1367      &                0.0, 0.00, 0.0, 0.0,    ! omega omega
1368      &                0.8, 0.22, 4.0, 1.0,    ! omega phi
1369      &              6.E-3, 0.80, 2.5, 0.0,    ! omega J/psi
1370      &                0.2, 0.22, 2.0, 1.0,    ! phi phi
1371      &              3.E-3, 0.80, 1.5, 0.0,    ! phi J/psi
1372      &              6.E-5, 2.00, 1.5, 0.0/    ! J/psi J/psi
1373       END
1374 *
1375 *=======================================================================
1376 CDECK  ID>, EVPRINT.
1377       SUBROUTINE evprint(Lun,Iev,P3,Ygg,ifl,Eb,ia,iz,amas,ptag1,ptag2)
1378 *-----------------------------------------------------------------------
1379 *       Event output for internal use by GGHIC authors
1380 *-----------------------------------------------------------------------
1381       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
1382       DOUBLE PRECISION P,V
1383       SAVE  /PYJETS/
1384       REAL ptag1(4),ptag2(4),p3(5)
1385 C
1386       WRITE(lun,100) iev,n,Eb,p3(5),ygg,ifl,-ifl
1387       WRITE(lun,101) 90000+ia,iz,ptag1,amas
1388       WRITE(lun,101) 90000+ia,iz,ptag2,amas
1389       WRITE(lun,101) 90000+ 3, 0,(p3(j),j=1,5)
1390       IF (n.NE.0) THEN
1391         DO i=1,n
1392           ichg=pychge(k(i,2))/3
1393           WRITE(lun,101) k(i,2),ichg,(p(i,j),j=1,5)
1394         ENDDO
1395       ENDIF
1396       RETURN
1397   100 FORMAT(i6,i3,3e16.8,2x,2i8)
1398   101 FORMAT(i6,i3,3(1x,e14.8),e13.8,e13.8)
1399       END
1400 *
1401 *=======================================================================
1402 C
1403 CDECK  ID>, TWOINI. 
1404       SUBROUTINE TWOINI(LOUTX,IPAR,XPAR,YPAR)                           TWOG0965
1405 C***********************************************************************TWOG0966
1406 C                                                                       TWOG0967
1407 C     Initialises the gamma-gamma generator of Ion collisions TWOGAM.   TWOG0968
1408 C                                                                       TWOG0969
1409 C   INPUTS:             LOUTX    Logical unit for print-out             TWOG0970
1410 C   ======= XPAR(1)     EBM      Beam energy                        GeV TWOG0971
1411 C           XPAR(2)     WMN      Minimum two-gamma mass             GeV TWOG0972
1412 C           XPAR(3)     WMX      Maximum two-gamma mass             GeV TWOG0973
1413 C                                                                       TWOG0974
1414 C           XPAR(8)     WGHTMX   Maximum weight                         TWOG0978
1415 C           XPAR(9)     XMRES    Mass of resonance                  GeV TWOG0979
1416 C           XPAR(10)    XGRES    Total width of resonance           GeV TWOG0980
1417 C                                                                       TWOG0981
1418 C           IPAR(1)     0        Unweighted events                      TWOG0982
1419 C                       1        Weighted events                        TWOG0983
1420 C           IPAR(2)     0        Continuum production                   TWOG0984
1421 C                       1        Resonance formation                    TWOG0985
1422 C                                                                       TWOG0986
1423 C   EXTERNALS:   none                                                   TWOG0987
1424 C   ==========                                                          TWOG0988
1425 C                                                                       TWOG0989
1426 C   AUTHOR/DATE:  S.A. Sadovsky,   17 January 1995                      TWOG0990
1427 C   ===========                                                         TWOG0991
1428 C                                                                       TWOG0992
1429 C***********************************************************************TWOG0993
1430 C     IMPLICIT NONE                                                     TWOG0994
1431       INTEGER   ICNTER,NCNTRS,NEVTS,LOUT,DEBUG                          TWOG0995
1432       PARAMETER (NCNTRS =7)                                             TWOG0996
1433       COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX,        TWOG0997
1434      +   ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS,     TWOG0998
1435      +   WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES,                       TWOG0999
1436      +   ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG                     TWOG1000
1437       SAVE   /TWGENC/                                                   TWOG1001
1438       REAL      ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT               TWOG1002
1439       REAL      WMIN,WMAX,OMMIN,OMRAT,EPS,TH1PMN,TH1PMX                 TWOG1003
1440       REAL      TH2PMN,TH2PMX,S,EB,WGHTSM,WGHTMX,WGHTGT,WTMXFD          TWOG1004
1441       REAL      XMRES,XGRES                                             TWOG1005
1442       LOGICAL   LRES,LWATE                                              TWOG1006
1443 C                                                                       TWOG1007
1444       REAL      ALPHA,PI,TWOPI,XME,XMU                                  TWOG1008
1445       PARAMETER (ALPHA=1./137.036,PI=3.14159265)                        TWOG1009
1446       PARAMETER (TWOPI=2.*PI,XME=0.000511,XMU=0.105658)                 TWOG1010
1447 C                                                                       TWOG1011
1448       REAL     XPAR(11),YPAR(6)                                         TWOG1012
1449       INTEGER  LOUTX,IPAR(10),I                                         TWOG1013
1450 C                                                                       TWOG1014
1451       EB    = XPAR(1)                                                   TWOG1015
1452       WMIN  = XPAR(2)                                                   TWOG1016
1453       WMAX  = XPAR(3)                                                   TWOG1017
1454       WGHTMX= XPAR(8)                                                   TWOG1022
1455       WTMXFD= 0                                                         TWOG1023
1456       XMRES = XPAR(9)                                                   TWOG1024
1457       XGRES = XPAR(10)                                                  TWOG1025
1458       LOUT  = LOUTX                                                     TWOG1026
1459       LWATE = IPAR(1) .EQ. 1                                            TWOG1027
1460       LRES  = IPAR(2) .EQ. 1                                            TWOG1028
1461       EPS   = 1./64.*(XME*WMIN**2/EB**3)**2/100.                        TWOG1029
1462 C+                                                                      TWOG1030
1463 C EPS is a small constant. We throw 1/2 COS(th/2)/(eps+SIN(th/2))       TWOG1031
1464 C rather than 1/2 COT(th/2) to be able to start theta from zero.        TWOG1032
1465 C-                                                                      TWOG1033
1466       S     = (2.*EB)**2                                                TWOG1054
1467       OMMIN = WMIN**2/(4.*EB)                                           TWOG1055
1468       OMRAT = EB/OMMIN                                                  TWOG1056
1469 C+                                                                      TWOG1057
1470 C Print out initial conditions.                                         TWOG1058
1471 C-                                                                      TWOG1059
1472       WRITE(LOUT,1000)                                                  TWOG1060
1473  1000 FORMAT(1H ,10X,27('-'),' TWINIT ',27('-'))                        TWOG1061
1474       WRITE(LOUT,1001)                                                  TWOG1062
1475  1001 FORMAT(1H ,10X,'I',60X,'I')                                       TWOG1063
1476 C
1477       WRITE(LOUT,1002)EB,YPAR(3),YPAR(4),YPAR(1),YPAR(2),WMIN,WMAX,
1478      +                                   YPAR(5),YPAR(6)
1479 C                                                                       TWOG1065
1480  1002 FORMAT                                                            TWOG1066
1481      +  (1H ,10X,'I    Beam energy               ',E15.5,               TWOG1067
1482      +  ' GeV',11X,'I',/,                                               TWOG1068
1483      +  1H ,10X,'I    Mass of ion                  ',F12.5,             TWOG1067
1484      +  ' GeV',11X,'I',/,                                               TWOG1068
1485      +  1H ,10X,'I    Gamma of ion                 ',F12.5,             TWOG1067
1486      +  '    ',11X,'I',/,                                               TWOG1068
1487      +  1H ,10X,'I    Charge of ion                ',F12.5,             TWOG1067
1488      +  '    ',11X,'I',/,                                               TWOG1068
1489      +  1H ,10X,'I    Atomic number of ion         ',F12.5,             TWOG1067
1490      +  '    ',11X,'I',/,                                               TWOG1068
1491      +  1H ,10X,'I    Minimum gamma-gamma mass     ',F12.5,             TWOG1069
1492      +  ' GeV',11X,'I',/,                                               TWOG1070
1493      +  1H ,10X,'I    Maximum gamma-gamma mass     ',F12.5,             TWOG1071
1494      +  ' GeV',11X,'I',/,                                               TWOG1080
1495      +  1H ,10X,'I    Minimum rapidity             ',F12.5,             TWOG1069
1496      +  '    ',11X,'I',/,                                               TWOG1070
1497      +  1H ,10X,'I    Maximum rapidity             ',F12.5,             TWOG1071
1498      +  '    ',11X,'I'  )                                               TWOG1080
1499       IF (LRES) THEN                                                    TWOG1081
1500       WRITE(LOUT,1001)                                                  TWOG1082
1501       WRITE(LOUT,1003)XMRES,XGRES                                       TWOG1083
1502       ENDIF                                                             TWOG1084
1503  1003 FORMAT                                                            TWOG1085
1504      +  (1H ,10X,'I    Resonance mass               ',F12.5,            TWOG1086
1505      +  ' GeV',11X,'I',/,                                               TWOG1087
1506      +  1H ,10X,'I    Resonance total width        ',F12.5,             TWOG1088
1507      +  ' GeV',11X,'I')                                                 TWOG1089
1508 C
1509       IF (LWATE) THEN                                                   TWOG1090
1510       WRITE(LOUT,1001)                                                  TWOG1091
1511       WRITE(LOUT,1004)                                                  TWOG1092
1512       ELSE                                                              TWOG1095
1513       WRITE(LOUT,1001)                                                  TWOG1096
1514       WRITE(LOUT,1005)WGHTMX                                            TWOG1097
1515  1004 FORMAT(1H ,10X,'I    Produce weighted events      ',27X,'I')      TWOG1094
1516  1005 FORMAT(1H ,10X,'I    Maximum weight               ',F12.5,        TWOG1099
1517      +                                            '    ',11X,'I')       TWOG1100
1518       ENDIF                                                             TWOG1101
1519 C                                                                       TWOG1102
1520       WRITE(LOUT,1001)                                                  TWOG1103
1521       WRITE(LOUT,1000)                                                  TWOG1104
1522 C                                                                       TWOG1105
1523       DO 10 I=1,NCNTRS                                                  TWOG1106
1524    10 ICNTER(I)=0                                                       TWOG1107
1525       WGHTGT = 0.                                                       TWOG1108
1526       WGHTSM = 0.                                                       TWOG1109
1527       RETURN                                                            TWOG1110
1528       END                                                               TWOG1111
1529 C=========================================================================
1530 CDECK  ID>, DLDMX.  
1531       FUNCTION DLDMX(Rm,Ny)
1532       PARAMETER (isize=100000)
1533 C                                                      S.A.Sadovsky
1534 C   Calculate the Luminosity integral over dY            5.02.1995
1535 C             in the rapidity range ymin-ymax
1536 C             using luminosity interpolation table
1537 C      Ny/2 - number of subintervals for Simpson integration
1538 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1539       COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
1540      +                           xmin,xmax,amn,amx,nnms,alfun(isize)
1541       stpy=(ymx-ymn)/nny
1542       stpm=(amx-amn)/nnms                 ! Interpolation steps
1543 C
1544       im=(rm-amn)/stpm
1545       q =(rm-im*stpm-amn)/stpm
1546 C
1547       My = Ny/2
1548       Hy =(ymax-ymin)/My                  ! Integration step over Y
1549       ry = ymin
1550       iy =(ry-ymn)/stpy+1.
1551       p  =(ry-(iy-1)*stpy-ymn)/stpy       ! Four point interpolation:  1
1552       k  = iy+im*(nny+1)
1553       F1 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1554      *        q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1555 C
1556       SIM = 0.0
1557       DO 10 i=1,My
1558 C
1559         ry = ymin + Hy*i
1560         iy =(ry-ymn)/stpy+1.
1561         p  =(ry-(iy-1)*stpy-ymn)/stpy       ! Four point interpolation:  3
1562         k  = iy+im*(nny+1)
1563         F3 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1564      *        q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1565 C
1566         ry = ry   - Hy*0.5
1567         iy =(ry-ymn)/stpy+1.
1568         p  =(ry-(iy-1)*stpy-ymn)/stpy       ! Four point interpolation:  2
1569         k  = iy+im*(nny+1)
1570         F2 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
1571      *        q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
1572 C
1573         SIM = SIM + (F1 + 4.*F2 + F3)/6.
1574    10 F1  = F3
1575 C
1576       DLDMX = 1000.*SIM*HY                ! dL/dMx  in GeV^-1
1577       RETURN
1578       END
1579 *=======================================================================
1580 CDECK  ID>, LUMFUN. 
1581       SUBROUTINE LUMFUN(ifl,file,ymin,ymax,ny,xmn,xmx,nms)
1582 C-----------------------------------------------------------------------
1583 C preparation of luminosity function for generator
1584 C
1585 C     ifl=-1 - calculate function + save in file 'file'
1586 C     ifl= 0 - calculate function
1587 C     ifl= 1 - read function from file
1588 C       ymin,ymax,ny - rapidity range & number of points
1589 C       xmn,xmx,nms  - mass range & number of points
1590 C       Attention: dimensions of alfun are alfun(ny+1,nms+1)
1591 C
1592 C     Called by GGINIT
1593 C
1594 C     Author: G.V.Khaoustov 12.12.1994
1595 C     Modified: Yu.Kharlov 18-JUN-1997
1596 C-----------------------------------------------------------------------
1597       COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
1598      &               gvconst(4,10)
1599       REAL nc
1600       COMMON/D2LParam/lz,la,Rion,Gamma
1601       INTEGER lz,la
1602       REAL    Rion,Gamma
1603       PARAMETER (isize=100000)
1604       CHARACTER*(*) file
1605       REAL work(1000)
1606       COMMON/lumfunct/anorm,cexp,umn,umx,ymn,ymx,nny,
1607      +                           wmn,wmx,amn,amx,nnms,alfun(isize)
1608 C
1609       umn = ymin                    !  parameters for luminosity generation
1610       umx = ymax
1611       wmn = xmn
1612       wmx = xmx
1613 C
1614 C Define R : radii of the ions (in fm) (if set non-zero, use it)
1615       IF (Rion .LE. 0) THEN
1616         rion = 1.2*float(la)**(1.0/3.0)
1617       ENDIF
1618 *
1619 *       Fing max energy of gg-system (with safity factor 2)
1620 *
1621       egg_max = 2.*hbarc * gamma / rion
1622 *
1623       IF(ifl.EQ.1) THEN             !  read lum. function from file
1624         OPEN(unit=10,err=1,file=file,status='old',
1625      *                                        form='formatted')
1626         READ(10, *,err=3) anorm,cexp
1627         READ(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
1628         READ(10,100,err=3) (alfun(I),i=1,(nny+1)*(nnms+1))
1629         CLOSE(unit=10)
1630         PRINT *,'The Luminosity has been read from File '
1631         RETURN
1632       ENDIF
1633 c
1634       IF (ny*nms .GT. isize) GOTO 2
1635       PRINT *,'Start the Luminosity calculation'
1636       nny = ny
1637       ymn = ymin
1638       ymx = ymax
1639       nnms= nms
1640       amn = xmn
1641       amx = xmx
1642 c
1643       stepy=(ymx-ymn)/nny
1644       stepm=(amx-amn)/nnms
1645 *
1646 *         Lum. function calculation
1647 *
1648       DO IM = 0,nnms
1649         am=amn+im*stepm
1650 *         Find Y-limits for mass AM (Yu.Kharlov)
1651         ym_max  = alog (egg_max/am)
1652         ym_max  = max(ym_max,5.*stepy)
1653         ygg_max = min (ymax, ym_max)
1654         ygg_min = max (ymin,-ym_max)
1655         DO IY = 1,nny+1
1656           Y = ymn+(iy-1)*stepy
1657           IF (Y.LT.ygg_max .AND. Y.GT.ygg_min) THEN
1658             dl = D2LDMDY(1000.*am,Y)   ! Differential Luminosity in MeV^-1
1659           ELSE
1660             dl = 0.0
1661           END IF
1662           alfun(iy+(nny+1)*im) = DL
1663         ENDDO
1664       ENDDO
1665 c
1666       s=0.
1667       DO i=1,nnms+1   ! find max values for exp. parameters calcul.
1668         am=0.
1669         DO j=1,nny+1
1670           IF(am.LT.alfun(j+(nny+1)*(i-1))) am=alfun(j+(nny+1)*(i-1))
1671         ENDDO
1672         work(i)=am
1673         s=s+am
1674       ENDDO
1675       s=(s-work(nnms+1))*stepm
1676 c
1677       smax=0.      ! calculation optimized parameters for exp. generator
1678       iflag=-1
1679       DO 101 i=1,nnms
1680         DO 102 j=i+1,nnms+1
1681           IF (work(j) .LE. 0.0) GOTO 102
1682           c=alog(work(i)/work(j))/((i-j)*stepm)
1683           am=work(i)/exp(c*(i-1)*stepm)
1684           se=am/c*(exp(c*(amx-amn))-1.)
1685           IF(s/se.lt.smax) GOTO 102
1686           DO k=1,nnms+1
1687             amp=am*exp(c*(k-1)*stepm)
1688             IF((work(k)-amp)/amp.GT.0.00001) THEN  !  Accuracy problem
1689               GOTO 102
1690             ENDIF
1691           ENDDO
1692           smax =s/se
1693           anorm=am
1694           cexp =c
1695           iflag=1
1696   102   CONTINUE
1697   101 CONTINUE
1698       IF(iflag.LT.0) THEN
1699         PRINT *,'Fault to find best exponential approximation'
1700         STOP
1701       ENDIF
1702 c
1703       IF(ifl.EQ.0) RETURN        ! save lum. function
1704       OPEN(unit=10,err=1,file=file,status='unknown',
1705      *                                          form='formatted')
1706       WRITE(10, *,err=3) anorm,cexp
1707       WRITE(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
1708       WRITE(10,100,err=3) (alfun(i),i=1,(nny+1)*(nnms+1))
1709       CLOSE(unit=10)
1710       PRINT *,'The Luminosity calculation is finished'
1711       RETURN
1712     1 PRINT *,' file open error, file:',file
1713       STOP
1714     2 PRINT *,' array "ALFUN" to small to contain luminosity function'
1715       STOP
1716     3 PRINT *,' file read/writr error, file:',file
1717       STOP
1718   100 FORMAT(5e13.5)
1719    99 FORMAT(2e13.5,i5,2e13.5,i5)
1720       END
1721 C======================================================================
1722 CDECK  ID>, D2LDMDY.
1723       REAL FUNCTION  D2LDMDY(M,Y)
1724 C-----------------------------------------------------------------------
1725 C Calculate the Luminosity-function d^2L/dMdY
1726 C where M is the invariant mass and Y the rapidity
1727 C of the photon-photon-system
1728 C
1729 C Calculations are based on the formulae of
1730 C N. Baron, G. Baur, Phys. Rev. C
1731 C C. Cahn, J. D. Jackson, Nucl. Phys.
1732 C
1733 C Kai Hencken, 29. 9. 1994
1734 C
1735 C M : invariant mass (MeV)
1736 C Y : rapidity
1737 C
1738 C D2LDMDY : diff. Lum. (MeV**-1)
1739 C
1740 C     Called by LUMFUN
1741 C-----------------------------------------------------------------------
1742       IMPLICIT NONE
1743       COMMON/D2LParam/lz,la,Rion,Gamma
1744       INTEGER lz,la
1745       REAL    Rion,Gamma
1746       REAL M,Y
1747       REAL RM,W1,W2
1748       COMMON/D2LIParam/RM,W1,W2
1749
1750       REAL NClass,DeltaL
1751       EXTERNAL NClass,DeltaL
1752
1753 C
1754 C physical constants hbar*c and 1/alpha
1755 C
1756       REAL hbarc,falphai
1757       PARAMETER (hbarc=197.327053,falphai=137.0359895)
1758
1759       REAL LClass,DL
1760 C
1761 C calculate photon frequencies from M,Y
1762 C
1763       W1 = M/2.0 * EXP ( Y)
1764       W2 = M/2.0 * EXP (-Y)
1765 C
1766 C calculate Rho in MeV^-1 instead of fm
1767 C
1768       RM = Rion / hbarc
1769 C
1770 C calculate the "classical" value first.
1771 C
1772
1773       LClass = NClass (W1,RM,Gamma) * NClass (W2,RM,Gamma)
1774
1775 C
1776 C substract from this the Delta-value
1777 C
1778       DL = DeltaL ()
1779
1780       D2LDMDY = 2.0 / M * float(LZ)**4 / falphai**2 * (LClass - DL)
1781
1782       RETURN
1783       END
1784 *=======================================================================
1785 CDECK  ID>, NCLASS. 
1786       REAL FUNCTION NClass (W, R, Gamma)
1787 C
1788 C the classical photon-number without impact parameter cutoff:
1789 C
1790       IMPLICIT NONE
1791       REAL W,R, Gamma
1792
1793       REAL Pi
1794       PARAMETER (Pi = 3.141592)
1795
1796       REAL BesselK0,BesselK1
1797       EXTERNAL BesselK0,BesselK1
1798
1799       REAL Xi
1800
1801       Xi = W*R/Gamma
1802
1803       NClass = 2.0/PI * (Xi*BesselK0 (Xi)*BesselK1 (Xi) -
1804      &          Xi**2/2.0 *(BesselK1(Xi)**2 - BesselK0(Xi)**2))
1805
1806       RETURN
1807       END
1808 *=======================================================================
1809 CDECK  ID>, DELTAL. 
1810       REAL FUNCTION DeltaL ()
1811 C
1812 C the difference to the classical value
1813 C     Called by d2LdMdY
1814 C
1815       IMPLICIT NONE
1816       COMMON/D2LParam/lz,la,Rion,Gamma
1817       INTEGER lz,la
1818       REAL    Rion,Gamma
1819       REAL RM,W1,W2
1820       COMMON/D2LIParam/RM,W1,W2
1821
1822       DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1823       COMMON /GAUSSD/XGAUSS,WGAUSS
1824
1825       INTEGER N,I
1826       REAL    b1
1827       REAL*8  t,tmin,tmax
1828       REAL*8  Sum,Int,Int2
1829
1830       REAL Itgb1
1831       EXTERNAL Itgb1
1832
1833       REAL   Pi
1834       REAL*8 Accur
1835       PARAMETER (Pi = 3.141592,Accur = 1D -2)
1836
1837 C
1838 C integrate first over b1
1839 C
1840
1841 C
1842 C Loop incrementing the boundary
1843 C
1844       tmin = 0.00
1845       tmax = 0.25
1846       Sum  = 0.00
1847
1848   100 CONTINUE
1849 C
1850 C Loop for the Gauss integration
1851 C
1852       Int=0.0
1853       DO N=1,6
1854         Int2 = Int
1855         Int=0.0
1856         DO I=2**N-1,2**(N+1)-2
1857           t   = (tmax-tmin)/2.0*XGAUSS(I) + (tmax+tmin)/2.0
1858           b1  = RM * DEXP (t)
1859           Int = Int + WGAUSS(I) * Itgb1(b1) * b1**2
1860         ENDDO
1861         Int = (tmax-tmin)/2.0*Int
1862         IF (Int .NE. 0.) THEN
1863           IF (DABS ((Int2-Int)/Int) .LT. Accur) GOTO 200
1864         END IF
1865       ENDDO
1866       IF (int .LT. 1.D-20) THEN
1867         int = 0.0D0
1868       ELSE
1869         WRITE(*,*) ' (b1) GAUSS MAY BE INACCURATE'
1870       END IF
1871
1872   200 CONTINUE
1873       Sum = Sum + Int
1874       IF (sum .NE. 0.0) THEN
1875         IF (ABS (Int2/Sum) .LT. Accur) GOTO 300
1876       ELSE
1877         IF (b1 .GT. 1.E+06)            GOTO 300
1878       END IF
1879       tmin = tmax
1880       tmax = tmax + 0.5
1881       GOTO 100
1882
1883   300 CONTINUE
1884
1885       DeltaL = 4.0*Pi * Sum
1886       RETURN
1887       END
1888
1889 CDECK  ID>, ITGB1.  
1890       REAL FUNCTION Itgb1(b)
1891 C
1892 C integration then over b2
1893 C     Called by DeltaL
1894 C
1895 C                                  Corrected by S.A.Sadovsky
1896 C                                             25.10.1994
1897       IMPLICIT NONE
1898 c  +seq,d2lpar. Not used Yu.Kh. 18-JUN-1997
1899       REAL     b
1900       REAL*8   b1,bmin,bmax
1901       REAL     RM,W1,W2
1902       COMMON/D2LIParam/RM,W1,W2
1903       DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
1904       COMMON /GAUSSD/XGAUSS,WGAUSS
1905       INTEGER  N,I
1906       REAL*8   b2,Int,Int2
1907       REAL*8   Itgb2
1908       EXTERNAL Itgb2
1909       REAL*8   Accur
1910       PARAMETER (Accur = 1.D-2)
1911
1912       b1 = b
1913       bmin = b1 - 2.0*RM
1914       IF (RM .GT. bmin) THEN
1915         bmin = RM
1916       ENDIF
1917       bmax = b1 + 2.0 * RM
1918       Int = 0.0
1919       DO N=1,6
1920         Int2 = Int
1921         Int = 0.0
1922         DO I=2**N-1,2**(N+1)-2
1923           b2 = (bmax-bmin)/2.0*XGAUSS(I) + (bmax+bmin)/2.0
1924           Int = Int + WGAUSS(I) * Itgb2(b1,b2) * b2
1925         END DO
1926         Int = (bmax-bmin)/2.0*Int
1927         IF (Int.NE.0.) THEN
1928           IF (DABS((Int2 - Int)/Int) .LT. Accur) GOTO 100
1929         ENDIF
1930       ENDDO
1931       IF (int .LT. 1.D-20) THEN
1932         int = 0.0D0
1933       ELSE
1934         WRITE(*,*) ' (b2) GAUSS MAY BE INACCURATE, Itgb1=',Int
1935       END IF
1936   100 CONTINUE
1937       Itgb1 = Int
1938       RETURN
1939       END
1940
1941 CDECK  ID>, ITGB2.  
1942       REAL*8 FUNCTION Itgb2 (b1,b2)
1943 C
1944 C       The function to be integrated over
1945 C       Called by Itgb1
1946 C
1947 C                                  Corrected by S.A.Sadovsky
1948 C                                             25.10.1994
1949       IMPLICIT NONE
1950       COMMON/D2LParam/lz,la,Rion,Gamma
1951       INTEGER lz,la
1952       REAL    Rion,Gamma
1953       REAL*8 b1,b2,arg
1954       REAL RM,W1,W2
1955       COMMON/D2LIParam/RM,W1,W2
1956
1957       REAL Nphoton
1958       EXTERNAL Nphoton
1959
1960       arg = (b1**2 + b2**2 - 4.0 * RM**2) / (2.0*b1*b2)
1961       IF (arg .GT. 1.D0) arg = 1.D0
1962       Itgb2 = Nphoton(W1,sngl(b1),Gamma) *
1963      &        Nphoton(W2,sngl(b2),Gamma) *
1964      &        DACOS (arg)
1965
1966       RETURN
1967       END
1968
1969 CDECK  ID>, NPHOTON.
1970       REAL FUNCTION Nphoton (W,Rho,Gamma)
1971 C
1972 C       The differential photonnumber for a nucleus
1973 C       (without form factor)
1974 C
1975       IMPLICIT NONE
1976       REAL W,Rho,Gamma
1977
1978       REAL BESSELK1
1979       EXTERNAL BESSELK1
1980
1981       REAL Wphib,WGamma
1982
1983       REAL PI
1984       PARAMETER (PI=3.141592)
1985
1986       WGamma = W/Gamma
1987 C
1988 C without form factor
1989 C
1990       Wphib = WGamma*BESSELK1 (WGamma*Rho)
1991
1992       Nphoton = 1.0/PI**2 * Wphib**2
1993
1994       RETURN
1995       END
1996
1997 CDECK  ID>, BESSELK0.   
1998       REAL FUNCTION BESSELK0(X)
1999 C
2000 C  Special functions I0,K0,I1,K1
2001 C  see Abramowitz&Stegun
2002 C
2003       IMPLICIT NONE
2004       REAL X,Y
2005
2006       REAL BESSELI0
2007       EXTERNAL BESSELI0
2008
2009       IF (X .LT. 2.0) THEN
2010
2011         Y = X**2/4.0
2012
2013         BESSELK0 =
2014      &    (-LOG(X/2.0)*BESSELI0(X))+(-.57721566+Y*(0.42278420
2015      &    +Y*(0.23069756+Y*(0.3488590E-1+Y*(0.262698E-2
2016      &    +Y*(0.10750E-3+Y*0.740E-5))))))
2017
2018       ELSE
2019
2020         Y = 2.0/X
2021
2022         BESSELK0 =
2023      &    (EXP(-X)/SQRT(X))*(1.25331414+Y*(-0.7832358E-1
2024      &    +Y*(0.2189568E-1+Y*(-0.1062446E-1+Y*(0.587872E-2
2025      &    +Y*(-0.251540E-2+Y*0.53208E-3))))))
2026
2027       ENDIF
2028
2029       RETURN
2030       END
2031
2032
2033 CDECK  ID>, BESSELI0.   
2034       REAL FUNCTION BESSELI0(X)
2035       IMPLICIT NONE
2036
2037       REAL X,Y,AX
2038
2039       AX = ABS(X)
2040
2041       IF (AX .LT. 3.75) THEN
2042
2043         Y = (X/3.75)**2
2044
2045         BESSELI0 =
2046      &    1.0+Y*(3.5156229+Y*(3.0899424+Y*(1.2067492
2047      &    +Y*(0.2659732+Y*(0.360768E-1+Y*0.45813E-2)))))
2048
2049       ELSE
2050
2051         Y = 3.75/AX
2052
2053         BESSELI0 =
2054      &    (EXP(AX)/SQRT(AX))*(0.39894228+Y*(0.1328592E-1
2055      &    +Y*(0.225319E-2+Y*(-0.157565E-2+Y*(0.916281E-2
2056      &    +Y*(-0.2057706E-1+Y*(0.2635537E-1+Y*(-0.1647633E-1
2057      &    +Y*0.392377E-2))))))))
2058
2059       ENDIF
2060
2061       RETURN
2062       END
2063
2064
2065 CDECK  ID>, BESSELK1.   
2066       REAL FUNCTION BESSELK1(X)
2067       IMPLICIT NONE
2068
2069       REAL X,Y
2070
2071       REAL BESSELI1
2072       EXTERNAL BESSELI1
2073
2074       IF (X .LT. 2.0) THEN
2075
2076         Y = X**2/4.0
2077         BESSELK1 =
2078      &    (LOG(X/2.0)*BESSELI1(X))+(1.0/X)*(1.0+Y*(0.15443144
2079      &    +Y*(-0.67278579+Y*(-0.18156897+Y*(-0.1919402E-1
2080      &    +Y*(-0.110404E-2+Y*(-0.4686E-4)))))))
2081
2082       ELSE
2083
2084         Y=2.0/X
2085         BESSELK1 =
2086      &    (EXP(-X)/SQRT(X))*(1.25331414+Y*(0.23498619
2087      &    +Y*(-0.3655620E-1+Y*(0.1504268E-1+Y*(-0.780353E-2
2088      &    +Y*(0.325614E-2+Y*(-0.68245E-3)))))))
2089
2090       ENDIF
2091
2092       RETURN
2093       END
2094
2095 CDECK  ID>, BESSELI1.   
2096       REAL FUNCTION BESSELI1(X)
2097       IMPLICIT NONE
2098
2099       REAL X,Y,AX
2100
2101       AX = ABS(X)
2102
2103       IF (AX .LT. 3.75) THEN
2104
2105         Y = (X/3.75)**2
2106
2107         BESSELI1 =
2108      &    AX*(0.5+Y*(0.87890594+Y*(0.51498869+Y*(0.15084934
2109      &    +Y*(0.2658733E-1+Y*(0.301532E-2+Y*0.32411E-3))))))
2110
2111
2112       ELSE
2113
2114         Y = 3.75/AX
2115
2116         BESSELI1 =
2117      &    0.2282967E-1+Y*(-0.2895312E-1+Y*(0.1787654E-1
2118      &    -Y*0.420059E-2))
2119
2120         BESSELI1 =
2121      &    0.39894228+Y*(-0.3988024E-1+Y*(-0.362018E-2
2122      &    +Y*(0.163801E-2+Y*(-0.1031555E-1+Y*BESSELI1))))
2123
2124         BESSELI1 = BESSELI1 *
2125      &    EXP(AX)/SQRT(AX)
2126
2127       ENDIF
2128
2129       IF (X .LT. 0) THEN
2130         BESSELI1 = -BESSELI1
2131       ENDIF
2132
2133       RETURN
2134       END
2135
2136 CDECK  ID>, GAUSSDAT.   
2137 C DATA for the Gauss integration, x-values and weight for a
2138 C N=2,4,8,16,32,64 point Gauss integration.
2139 C
2140 C based on program GAUSS64 of D. Trautmann
2141 C data from Abramowitz & Stegun
2142 C
2143 C KAI HENCKEN, FEBRUAR 1993
2144 C
2145       BLOCKDATA GAUSSDATA
2146       IMPLICIT NONE
2147       DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
2148       COMMON /GAUSSD/XGAUSS,WGAUSS
2149
2150       DATA XGAUSS(1)/ .57735026918962576D0/
2151       DATA XGAUSS(2)/-.57735026918962576D0/
2152       DATA WGAUSS(1)/ 1.00000000000000000D0/
2153       DATA WGAUSS(2)/ 1.00000000000000000D0/
2154
2155       DATA XGAUSS(3)/ .33998104358485627D0/
2156       DATA XGAUSS(4)/ .86113631159405258D0/
2157       DATA XGAUSS(5)/-.33998104358485627D0/
2158       DATA XGAUSS(6)/-.86113631159405258D0/
2159       DATA WGAUSS(3)/ .65214515486254613D0/
2160       DATA WGAUSS(4)/ .34785484513745385D0/
2161       DATA WGAUSS(5)/ .65214515486254613D0/
2162       DATA WGAUSS(6)/ .34785484513745385D0/
2163
2164       DATA XGAUSS(7)/ .18343464249564981D0/
2165       DATA XGAUSS(8)/ .52553240991632899D0/
2166       DATA XGAUSS(9)/ .79666647741362674D0/
2167       DATA XGAUSS(10)/ .96028985649753623D0/
2168       DATA XGAUSS(11)/-.18343464249564981D0/
2169       DATA XGAUSS(12)/-.52553240991632899D0/
2170       DATA XGAUSS(13)/-.79666647741362674D0/
2171       DATA XGAUSS(14)/-.96028985649753623D0/
2172       DATA WGAUSS(7)/ .36268378337836198D0/
2173       DATA WGAUSS(8)/ .31370664587788727D0/
2174       DATA WGAUSS(9)/ .22238103445337448D0/
2175       DATA WGAUSS(10)/ .10122853629037627D0/
2176       DATA WGAUSS(11)/ .36268378337836198D0/
2177       DATA WGAUSS(12)/ .31370664587788727D0/
2178       DATA WGAUSS(13)/ .22238103445337448D0/
2179       DATA WGAUSS(14)/ .10122853629037627D0/
2180
2181       DATA XGAUSS(15)/ .0950125098376374402D0/
2182       DATA XGAUSS(16)/ .281603550779258913D0/
2183       DATA XGAUSS(17)/ .458016777657227386D0/
2184       DATA XGAUSS(18)/ .617876244402643748D0/
2185       DATA XGAUSS(19)/ .755404408355003034D0/
2186       DATA XGAUSS(20)/ .865631202387831744D0/
2187       DATA XGAUSS(21)/ .944575023073232576D0/
2188       DATA XGAUSS(22)/ .989400934991649933D0/
2189       DATA XGAUSS(23)/-.0950125098376374402D0/
2190       DATA XGAUSS(24)/-.281603550779258913D0/
2191       DATA XGAUSS(25)/-.458016777657227386D0/
2192       DATA XGAUSS(26)/-.617876244402643748D0/
2193       DATA XGAUSS(27)/-.755404408355003034D0/
2194       DATA XGAUSS(28)/-.865631202387831744D0/
2195       DATA XGAUSS(29)/-.944575023073232576D0/
2196       DATA XGAUSS(30)/-.989400934991649933D0/
2197       DATA WGAUSS(15)/ .189450610455068496D0/
2198       DATA WGAUSS(16)/ .182603415044923589D0/
2199       DATA WGAUSS(17)/ .169156519395002538D0/
2200       DATA WGAUSS(18)/ .149595988816576732D0/
2201       DATA WGAUSS(19)/ .124628971255533872D0/
2202       DATA WGAUSS(20)/ .0951585116824927848D0/
2203       DATA WGAUSS(21)/ .0622535239386478929D0/
2204       DATA WGAUSS(22)/ .0271524594117540949D0/
2205       DATA WGAUSS(23)/ .189450610455068496D0/
2206       DATA WGAUSS(24)/ .182603415044923589D0/
2207       DATA WGAUSS(25)/ .169156519395002538D0/
2208       DATA WGAUSS(26)/ .149595988816576732D0/
2209       DATA WGAUSS(27)/ .124628971255533872D0/
2210       DATA WGAUSS(28)/ .0951585116824927848D0/
2211       DATA WGAUSS(29)/ .0622535239386478929D0/
2212       DATA WGAUSS(30)/ .0271524594117540949D0/
2213
2214       DATA XGAUSS(31)/ .0483076656877383162D0/
2215       DATA XGAUSS(32)/ .144471961582796493D0/
2216       DATA XGAUSS(33)/ .239287362252137075D0/
2217       DATA XGAUSS(34)/ .331868602282127650D0/
2218       DATA XGAUSS(35)/ .421351276130635345D0/
2219       DATA XGAUSS(36)/ .506899908932229390D0/
2220       DATA XGAUSS(37)/ .587715757240762329D0/
2221       DATA XGAUSS(38)/ .663044266930215201D0/
2222       DATA XGAUSS(39)/ .732182118740289680D0/
2223       DATA XGAUSS(40)/ .794483795967942407D0/
2224       DATA XGAUSS(41)/ .849367613732569970D0/
2225       DATA XGAUSS(42)/ .896321155766052124D0/
2226       DATA XGAUSS(43)/ .934906075937739689D0/
2227       DATA XGAUSS(44)/ .964762255587506430D0/
2228       DATA XGAUSS(45)/ .985611511545268335D0/
2229       DATA XGAUSS(46)/ .997263861849481564D0/
2230       DATA XGAUSS(47)/-.0483076656877383162D0/
2231       DATA XGAUSS(48)/-.144471961582796493D0/
2232       DATA XGAUSS(49)/-.239287362252137075D0/
2233       DATA XGAUSS(50)/-.331868602282127650D0/
2234       DATA XGAUSS(51)/-.421351276130635345D0/
2235       DATA XGAUSS(52)/-.506899908932229390D0/
2236       DATA XGAUSS(53)/-.587715757240762329D0/
2237       DATA XGAUSS(54)/-.663044266930215201D0/
2238       DATA XGAUSS(55)/-.732182118740289680D0/
2239       DATA XGAUSS(56)/-.794483795967942407D0/
2240       DATA XGAUSS(57)/-.849367613732569970D0/
2241       DATA XGAUSS(58)/-.896321155766052124D0/
2242       DATA XGAUSS(59)/-.934906075937739689D0/
2243       DATA XGAUSS(60)/-.964762255587506430D0/
2244       DATA XGAUSS(61)/-.985611511545268335D0/
2245       DATA XGAUSS(62)/-.997263861849481564D0/
2246       DATA WGAUSS(31)/ .0965400885147278006D0/
2247       DATA WGAUSS(32)/ .0956387200792748594D0/
2248       DATA WGAUSS(33)/ .0938443990808045654D0/
2249       DATA WGAUSS(34)/ .0911738786957638847D0/
2250       DATA WGAUSS(35)/ .0876520930044038111D0/
2251       DATA WGAUSS(36)/ .0833119242269467552D0/
2252       DATA WGAUSS(37)/ .0781938957870703065D0/
2253       DATA WGAUSS(38)/ .0723457941088485062D0/
2254       DATA WGAUSS(39)/ .0658222227763618468D0/
2255       DATA WGAUSS(40)/ .0586840934785355471D0/
2256       DATA WGAUSS(41)/ .0509980592623761762D0/
2257       DATA WGAUSS(42)/ .0428358980222266807D0/
2258       DATA WGAUSS(43)/ .0342738629130214331D0/
2259       DATA WGAUSS(44)/ .0253920653092620595D0/
2260       DATA WGAUSS(45)/ .0162743947309056706D0/
2261       DATA WGAUSS(46)/ .00701861000947009660D0/
2262       DATA WGAUSS(47)/ .0965400885147278006D0/
2263       DATA WGAUSS(48)/ .0956387200792748594D0/
2264       DATA WGAUSS(49)/ .0938443990808045654D0/
2265       DATA WGAUSS(50)/ .0911738786957638847D0/
2266       DATA WGAUSS(51)/ .0876520930044038111D0/
2267       DATA WGAUSS(52)/ .0833119242269467552D0/
2268       DATA WGAUSS(53)/ .0781938957870703065D0/
2269       DATA WGAUSS(54)/ .0723457941088485062D0/
2270       DATA WGAUSS(55)/ .0658222227763618468D0/
2271       DATA WGAUSS(56)/ .0586840934785355471D0/
2272       DATA WGAUSS(57)/ .0509980592623761762D0/
2273       DATA WGAUSS(58)/ .0428358980222266807D0/
2274       DATA WGAUSS(59)/ .0342738629130214331D0/
2275       DATA WGAUSS(60)/ .0253920653092620595D0/
2276       DATA WGAUSS(61)/ .0162743947309056706D0/
2277       DATA WGAUSS(62)/ .00701861000947009660D0/
2278
2279       DATA XGAUSS(63)/ .02435029266342443250D0/
2280       DATA XGAUSS(64)/ .0729931217877990394D0/
2281       DATA XGAUSS(65)/ .121462819296120554D0/
2282       DATA XGAUSS(66)/ .169644420423992818D0/
2283       DATA XGAUSS(67)/ .217423643740007084D0/
2284       DATA XGAUSS(68)/ .264687162208767416D0/
2285       DATA XGAUSS(69)/ .311322871990210956D0/
2286       DATA XGAUSS(70)/ .357220158337668116D0/
2287       DATA XGAUSS(71)/ .402270157963991604D0/
2288       DATA XGAUSS(72)/ .446366017253464088D0/
2289       DATA XGAUSS(73)/ .489403145707052957D0/
2290       DATA XGAUSS(74)/ .531279464019894546D0/
2291       DATA XGAUSS(75)/ .571895646202634034D0/
2292       DATA XGAUSS(76)/ .611155355172393250D0/
2293       DATA XGAUSS(77)/ .648965471254657340D0/
2294       DATA XGAUSS(78)/ .685236313054233243D0/
2295       DATA XGAUSS(79)/ .719881850171610827D0/
2296       DATA XGAUSS(80)/ .752819907260531897D0/
2297       DATA XGAUSS(81)/ .783972358943341408D0/
2298       DATA XGAUSS(82)/ .813265315122797560D0/
2299       DATA XGAUSS(83)/ .840629296252580363D0/
2300       DATA XGAUSS(84)/ .865999398154092820D0/
2301       DATA XGAUSS(85)/ .889315445995114106D0/
2302       DATA XGAUSS(86)/ .910522137078502806D0/
2303       DATA XGAUSS(87)/ .929569172131939576D0/
2304       DATA XGAUSS(88)/ .946411374858402816D0/
2305       DATA XGAUSS(89)/ .961008799652053719D0/
2306       DATA XGAUSS(90)/ .973326827789910964D0/
2307       DATA XGAUSS(91)/ .983336253884625957D0/
2308       DATA XGAUSS(92)/ .991013371476744321D0/
2309       DATA XGAUSS(93)/ .996340116771955279D0/
2310       DATA XGAUSS(94)/ .999305041735772139D0/
2311       DATA XGAUSS(95)/-.02435029266342443250D0/
2312       DATA XGAUSS(96)/-.0729931217877990394D0/
2313       DATA XGAUSS(97)/-.121462819296120554D0/
2314       DATA XGAUSS(98)/-.169644420423992818D0/
2315       DATA XGAUSS(99)/-.217423643740007084D0/
2316       DATA XGAUSS(100)/-.264687162208767416D0/
2317       DATA XGAUSS(101)/-.311322871990210956D0/
2318       DATA XGAUSS(102)/-.357220158337668116D0/
2319       DATA XGAUSS(103)/-.402270157963991604D0/
2320       DATA XGAUSS(104)/-.446366017253464088D0/
2321       DATA XGAUSS(105)/-.489403145707052957D0/
2322       DATA XGAUSS(106)/-.531279464019894546D0/
2323       DATA XGAUSS(107)/-.571895646202634034D0/
2324       DATA XGAUSS(108)/-.611155355172393250D0/
2325       DATA XGAUSS(109)/-.648965471254657340D0/
2326       DATA XGAUSS(110)/-.685236313054233243D0/
2327       DATA XGAUSS(111)/-.719881850171610827D0/
2328       DATA XGAUSS(112)/-.752819907260531897D0/
2329       DATA XGAUSS(113)/-.783972358943341408D0/
2330       DATA XGAUSS(114)/-.813265315122797560D0/
2331       DATA XGAUSS(115)/-.840629296252580363D0/
2332       DATA XGAUSS(116)/-.865999398154092820D0/
2333       DATA XGAUSS(117)/-.889315445995114106D0/
2334       DATA XGAUSS(118)/-.910522137078502806D0/
2335       DATA XGAUSS(119)/-.929569172131939576D0/
2336       DATA XGAUSS(120)/-.946411374858402816D0/
2337       DATA XGAUSS(121)/-.961008799652053719D0/
2338       DATA XGAUSS(122)/-.973326827789910964D0/
2339       DATA XGAUSS(123)/-.983336253884625957D0/
2340       DATA XGAUSS(124)/-.991013371476744321D0/
2341       DATA XGAUSS(125)/-.996340116771955279D0/
2342       DATA XGAUSS(126)/-.999305041735772139D0/
2343       DATA WGAUSS(63)/ .0486909570091397204D0/
2344       DATA WGAUSS(64)/ .0485754674415034269D0/
2345       DATA WGAUSS(65)/ .0483447622348029572D0/
2346       DATA WGAUSS(66)/ .0479993885964583077D0/
2347       DATA WGAUSS(67)/ .0475401657148303087D0/
2348       DATA WGAUSS(68)/ .0469681828162100173D0/
2349       DATA WGAUSS(69)/ .0462847965813144172D0/
2350       DATA WGAUSS(70)/ .0454916279274181445D0/
2351       DATA WGAUSS(71)/ .0445905581637565631D0/
2352       DATA WGAUSS(72)/ .0435837245293234534D0/
2353       DATA WGAUSS(73)/ .0424735151236535890D0/
2354       DATA WGAUSS(74)/ .0412625632426235286D0/
2355       DATA WGAUSS(75)/ .0399537411327203414D0/
2356       DATA WGAUSS(76)/ .0385501531786156291D0/
2357       DATA WGAUSS(77)/ .0370551285402400460D0/
2358       DATA WGAUSS(78)/ .0354722132568823838D0/
2359       DATA WGAUSS(79)/ .0338051618371416094D0/
2360       DATA WGAUSS(80)/ .0320579283548515535D0/
2361       DATA WGAUSS(81)/ .0302346570724024789D0/
2362       DATA WGAUSS(82)/ .0283396726142594832D0/
2363       DATA WGAUSS(83)/ .0263774697150546587D0/
2364       DATA WGAUSS(84)/ .0243527025687108733D0/
2365       DATA WGAUSS(85)/ .0222701738083832542D0/
2366       DATA WGAUSS(86)/ .0201348231535302094D0/
2367       DATA WGAUSS(87)/ .0179517157756973431D0/
2368       DATA WGAUSS(88)/ .0157260304760247193D0/
2369       DATA WGAUSS(89)/ .0134630478967186426D0/
2370       DATA WGAUSS(90)/ .0111681394601311288D0/
2371       DATA WGAUSS(91)/ .00884675982636394772D0/
2372       DATA WGAUSS(92)/ .00650445796897836286D0/
2373       DATA WGAUSS(93)/ .00414703326056246764D0/
2374       DATA WGAUSS(94)/ .00178328072169643295D0/
2375       DATA WGAUSS(95)/ .0486909570091397204D0/
2376       DATA WGAUSS(96)/ .0485754674415034269D0/
2377       DATA WGAUSS(97)/ .0483447622348029572D0/
2378       DATA WGAUSS(98)/ .0479993885964583077D0/
2379       DATA WGAUSS(99)/ .0475401657148303087D0/
2380       DATA WGAUSS(100)/ .0469681828162100173D0/
2381       DATA WGAUSS(101)/ .0462847965813144172D0/
2382       DATA WGAUSS(102)/ .0454916279274181445D0/
2383       DATA WGAUSS(103)/ .0445905581637565631D0/
2384       DATA WGAUSS(104)/ .0435837245293234534D0/
2385       DATA WGAUSS(105)/ .0424735151236535890D0/
2386       DATA WGAUSS(106)/ .0412625632426235286D0/
2387       DATA WGAUSS(107)/ .0399537411327203414D0/
2388       DATA WGAUSS(108)/ .0385501531786156291D0/
2389       DATA WGAUSS(109)/ .0370551285402400460D0/
2390       DATA WGAUSS(110)/ .0354722132568823838D0/
2391       DATA WGAUSS(111)/ .0338051618371416094D0/
2392       DATA WGAUSS(112)/ .0320579283548515535D0/
2393       DATA WGAUSS(113)/ .0302346570724024789D0/
2394       DATA WGAUSS(114)/ .0283396726142594832D0/
2395       DATA WGAUSS(115)/ .0263774697150546587D0/
2396       DATA WGAUSS(116)/ .0243527025687108733D0/
2397       DATA WGAUSS(117)/ .0222701738083832542D0/
2398       DATA WGAUSS(118)/ .0201348231535302094D0/
2399       DATA WGAUSS(119)/ .0179517157756973431D0/
2400       DATA WGAUSS(120)/ .0157260304760247193D0/
2401       DATA WGAUSS(121)/ .0134630478967186426D0/
2402       DATA WGAUSS(122)/ .0111681394601311288D0/
2403       DATA WGAUSS(123)/ .00884675982636394772D0/
2404       DATA WGAUSS(124)/ .00650445796897836286D0/
2405       DATA WGAUSS(125)/ .00414703326056246764D0/
2406       DATA WGAUSS(126)/ .00178328072169643295D0/
2407
2408       END
2409 *=======================================================================
2410       SUBROUTINE LORENB (U,PS,PI,PF)
2411 C
2412 C CERN PROGLIB# U102    LORENB          .VERSION KERNFOR  4.04  821124
2413 C ORIG. 20/08/75 L.PAPE
2414 C
2415       DOUBLE PRECISION PF4, FN
2416       DIMENSION      PS(4),PI(4),PF(4)
2417
2418       IF (PS(4).EQ.U) GO TO 17
2419       PF4  = (PI(4)*PS(4)+PI(3)*PS(3)+PI(2)*PS(2)+PI(1)*PS(1)) / U
2420       FN   = (PF4+PI(4)) / (PS(4)+U)
2421       PF(1)= PI(1) + FN*PS(1)
2422       PF(2)= PI(2) + FN*PS(2)
2423       PF(3)= PI(3) + FN*PS(3)
2424       PF(4)= PF4
2425       GO TO 18
2426 C
2427    17 PF(1)= PI(1)
2428       PF(2)= PI(2)
2429       PF(3)= PI(3)
2430       PF(4)= PI(4)
2431 C
2432    18 CONTINUE
2433 C
2434       RETURN
2435 C
2436       END