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