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