--- /dev/null
+CDECK ID>, GGCDES.
+*=======================================================================
+CDECK ID>, GGINIT.
+ SUBROUTINE gginit
+*-----------------------------------------------------------------------
+* GGHIC initialization
+* Author: Yu.Kharlov
+*-----------------------------------------------------------------------
+ COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
+ & amin,amax,ymin,ymax,nmas,ny, kferm,
+ & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
+ & thetamin, costhv1, kv1,kv2,gvpar(4)
+ CHARACTER lumfil*80
+ COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
+ & gvconst(4,10)
+ REAL nc
+ COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
+ & xstote, ssbr(10)
+ COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
+ & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
+ COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr,
+ & xmll, xmlr, xmnl, xtanb, xmha, xmu,
+ & xmt, xat, xmbr, xab, u11, v11
+ COMMON/D2LParam/lz,la,Rion,Gamma
+ INTEGER lz,la
+ REAL Rion,Gamma
+ COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
+ DOUBLE PRECISION CKIN
+ SAVE /PYSUBS/
+ COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
+ DOUBLE PRECISION PARP,PARI
+ SAVE /PYPARS/
+ COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
+ DOUBLE PRECISION P,V
+ SAVE /PYJETS/
+ COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
+ DOUBLE PRECISION PARU,PARJ
+ SAVE /PYDAT1/
+ COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
+ DOUBLE PRECISION PMAS,PARF,VCKM
+ SAVE /PYDAT2/
+ COMMON/PYDAT4/CHAF(500,2)
+ CHARACTER CHAF*16
+ SAVE /PYDAT4/
+C MXSS = maximum number of modes
+C NSSMOD = number of modes
+C ISSMOD = initial particle
+C JSSMOD = final particles
+C GSSMOD = width
+C BSSMOD = branching ratio
+ INTEGER MXSS
+ PARAMETER (MXSS=1000)
+ COMMON/SSMODE/NSSMOD,ISSMOD(MXSS),JSSMOD(5,MXSS),GSSMOD(MXSS)
+ $,BSSMOD(MXSS)
+ INTEGER NSSMOD,ISSMOD,JSSMOD
+ REAL GSSMOD,BSSMOD
+ SAVE /SSMODE/
+C SUSY parameters
+C AMGLSS = gluino mass
+C AMQKSS = squark mass
+C AMLLSS = left-slepton mass
+C AMLRSS = right-slepton mass
+C AMNLSS = sneutrino mass
+C TWOM1 = Higgsino mass = - mu
+C RV2V1 = ratio v2/v1 of vev's
+C AMTLSS,AMTRSS = left,right stop masses
+C AMT1SS,AMT2SS = light,heavy stop masses
+C AMBLSS,AMBRSS = left,right sbottom masses
+C AMB1SS,AMB2SS = light,heavy sbottom masses
+C AMZiSS = signed mass of Zi
+C ZMIXSS = Zi mixing matrix
+C AMWiSS = signed Wi mass
+C GAMMAL,GAMMAR = Wi left, right mixing angles
+C AMHL,AMHH,AMHA = neutral Higgs h0, H0, A0 masses
+C AMHC = charged Higgs H+ mass
+C ALFAH = Higgs mixing angle
+C AAT = stop trilinear term
+C THETAT = stop mixing angle
+C AAB = sbottom trilinear term
+C THETAB = sbottom mixing angle
+ COMMON/SSPAR/AMGLSS,AMULSS,AMURSS,AMDLSS,AMDRSS,AMSLSS
+ $,AMSRSS,AMCLSS,AMCRSS,AMBLSS,AMBRSS,AMB1SS,AMB2SS
+ $,AMTLSS,AMTRSS,AMT1SS,AMT2SS,AMELSS,AMERSS,AMMLSS,AMMRSS
+ $,AMLLSS,AMLRSS,AML1SS,AML2SS,AMN1SS,AMN2SS,AMN3SS
+ $,TWOM1,RV2V1,AMZ1SS,AMZ2SS,AMZ3SS,AMZ4SS,ZMIXSS(4,4)
+ $,AMW1SS,AMW2SS
+ $,GAMMAL,GAMMAR,AMHL,AMHH,AMHA,AMHC,ALFAH,AAT,THETAT
+ $,AAB,THETAB,AAL,THETAL,AMGVSS
+ REAL AMGLSS,AMULSS,AMURSS,AMDLSS,AMDRSS,AMSLSS
+ $,AMSRSS,AMCLSS,AMCRSS,AMBLSS,AMBRSS,AMB1SS,AMB2SS
+ $,AMTLSS,AMTRSS,AMT1SS,AMT2SS,AMELSS,AMERSS,AMMLSS,AMMRSS
+ $,AMLLSS,AMLRSS,AML1SS,AML2SS,AMN1SS,AMN2SS,AMN3SS
+ $,TWOM1,RV2V1,AMZ1SS,AMZ2SS,AMZ3SS,AMZ4SS,ZMIXSS
+ $,AMW1SS,AMW2SS
+ $,GAMMAL,GAMMAR,AMHL,AMHH,AMHA,AMHC,ALFAH,AAT,THETAT
+ $,AAB,THETAB,AAL,THETAL,AMGVSS
+ REAL AMZISS(4)
+ EQUIVALENCE (AMZISS(1),AMZ1SS)
+ SAVE /SSPAR/
+ REAL pgam1(4),pgam2(4),xpar(10),ypar(6)
+ INTEGER ipar(2)
+ INTEGER pycomp
+ REAL*8 ggrnd
+ EXTERNAL pydata
+ la = ia
+ lz = iz
+* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+*
+* Check bounds on 2 photon mass
+*
+ IF (iproc .EQ. 1) THEN
+ IF (amin .LE. 2.31) THEN
+ WRITE (6,*)
+ & ' %GGINIT: minimal 2-photon mass is too low, ',
+ & 'set to 3m(pho)'
+ amin = 2.31
+ END IF
+*
+ ELSE IF (iproc .EQ. 2) THEN
+ kc = pycomp (kf_onium)
+ xmres = pmas(kc,1)
+ IF (kf_onium .EQ. 441) THEN
+ xgtres = 0.013
+ ELSE
+ xgtres = pmas(kc,2)
+ END IF
+c. xggres = ggwid(kf_onium)
+ IF (xgtres .EQ. 0.) xgtres = 0.1
+ IF (xmres.LE.amin .OR. xmres.GE.amax) THEN
+ WRITE(6,*)
+ & ' %GGINIT: Xmres is out of the Luminosity mass range'
+ STOP
+ ENDIF
+*
+ ELSE IF (iproc .EQ. 3) THEN
+ kc = pycomp (kferm)
+ amprt(1) = pmas(kc,1)
+ IF (amin .LT. 2*amprt(1)) THEN
+ WRITE (6,*)
+ & ' %GGINIT: minimal 2-photon mass is too low, ',
+ & 'set to 2m(fermion)'
+ amin = 2*amprt(1)
+ END IF
+*
+ ELSE IF (iproc .EQ. 4) THEN
+ IF (amprt(4) .LE. 0.) amprt(4) = pmas(pycomp(24),1)
+ IF (amin .LT. 2*amprt(4)) then
+ WRITE (6,*)
+ & ' %GGINIT: minimal 2-photon mass is too low, ',
+ & 'set to 2m(W+-)'
+ amin = 2*amprt(4)
+ END IF
+*
+ ELSE IF (iproc .EQ. 5) THEN
+ WRITE (6,*) 'Process 5 is not implemented yet'
+ STOP
+C CALL ssmssm (xm1, xm2, xmg, xms, xmtl, xmtr, xmll, xmlr, xmnl,
+C & xtanb, xmha, xmu, xmt, xat, xmbr, xab, iallow)
+C IF (iallow .NE. 0) THEN
+C WRITE (6,*) 'Lightest neutralino is not LSP, ',
+C & 'input other MSSM parameters'
+C STOP
+C END IF
+C amprt(1) = abs (amw1ss)
+C amprt(2) = abs (amz1ss)
+C WRITE (6,'('' %GGINIT'','//
+C & ''', M(chi+) = '',f5.1,'//
+C & ''', M(chi0) = '',f5.1)') amprt(1), amprt(2)
+C *
+C IF (amin .LT. 2*amprt(1)) THEN
+C WRITE (6,*)
+C & ' %GGINIT: minimal 2-photon mass is too low, ',
+C & 'set to 2M(chi+)'
+C amin = 2*amprt(1)
+C END IF
+*
+ ELSE IF (iproc .EQ. 6) THEN
+ kc1 = pycomp (kv1)
+ amprt(1) = pmas(kc1,1) + pmas(kc1,3)
+ kc2 = pycomp (kv2)
+ amprt(2) = pmas(kc2,1) + pmas(kc2,3)
+ IF (amin .LT. amprt(1)+amprt(2)) THEN
+ WRITE (6,*)
+ & ' %GGINIT: minimal 2-photon mass is too low, ',
+ & 'set to m(kv1)+m(kv2)'
+ amin = amprt(1)+amprt(2)
+ END IF
+*
+ END IF
+*
+ IF (amin .GT. amax) THEN
+ WRITE (6,*) ' %GGINIT: AMIN > AMAX, execution stopped'
+ STOP
+ END IF
+*
+* ----- TWOGAM Initialization -----
+*
+ gamma = ebmn*ia/amas
+ xpar(1)= ebmn*ia ! Beam energy
+ xpar(2)= amin ! Minimum gg mass
+ xpar(3)= amax ! Maximum gg mass
+ xpar(8)= 100. ! Weight, but not too long
+*
+ ipar(1)= 1 ! Unweighted events
+ ipar(2)= 0 ! Continuum generation
+*
+ IF (iproc .EQ. 2) THEN
+ xpar(9) = xmres ! Resonance mass
+ xpar(10)= xgtres ! Resonance total width
+ ipar(2) = 1 ! Resonance generation
+ END IF
+*
+ ypar(1)= iz ! Ion charge
+ ypar(2)= ia ! Ion atomic number
+ ypar(3)= amas ! Ion mass
+ ypar(4)= gamma ! Gamma factor of ion
+ ypar(5)= ymin ! Minimum gg-rapidity
+ ypar(6)= ymax ! Maximum gg-rapidity
+ eb = xpar(1)
+*
+ CALL twoini (6,ipar,xpar,ypar)
+*
+* The Luminosity function initialization
+*
+* ilumf = +1 ! Reading the Luminosity from file
+* ilumf = -1 ! New calculation of the Luminosity function
+*
+ CALL lumfun (ilumf,lumfil,ymin,ymax,ny,amin,amax,nmas)
+*
+* Integral of lum. function over the (M,y) region
+ xlumint = 0.0
+ hmas = (amax-amin)/nmas
+ f1 = dldmx(amin ,ny)
+ DO i=1,nmas
+ xm = amin + i*hmas
+ f2 = dldmx(xm-hmas/2.,ny)
+ f3 = dldmx(xm ,ny)
+ xlumint = xlumint +(f1 + 4.*f2 + f3)/6.
+ f1 = f3
+ END DO
+ xlumint = xlumint * hmas
+* ----- end of TWOGAM initialization -----
+*
+ IF (iproc .EQ. 1) THEN
+* ----- PYTHIA initialization -----
+ parp(2) = dble(amin)
+ mstp(14) = 20
+ mstp(82) = 1
+ mstp(171) = 1
+ mstp(172) = 1
+ DO ip = 1,2
+ p(1,ip) = 0.D0
+ p(2,ip) = 0.D0
+ END DO
+ p(1,3) = dble( amax/2)
+ p(2,3) = dble(-amax/2)
+ p(1,4) = dble( amax/2)
+ p(2,4) = dble( amax/2)
+ p(1,5) = 0.D0
+ p(2,5) = 0.D0
+ CALL pyinit ('5MOM','gamma','gamma',0.0D0)
+* ----- end of PYTHIA initialization -----
+*
+ ELSE IF (iproc .EQ. 2) THEN
+* --- Cross section for narrow Resonance production ---
+ totxsc = dldmx(xmres,ny) ! dL/dMx in GeV^-1
+ njres = kf_onium / 10
+ njres = kf_onium - 10 * njres
+ xnorm = 4. * pi**2 * njres * xggres * gev2nb / xmres**2
+ xstot = xnorm * totxsc
+ xstote = 0.
+*
+ ELSE IF (iproc .EQ. 3) THEN
+ qf = abs(kchg(kc,1))/3.
+ nc = abs(kchg(kc,2))*2. + 1.
+ ELSE IF (iproc .EQ. 5) THEN
+* Chargino storing in JETSET
+ kc = pycomp (41)
+ chaf(kc,1) = 'chi_1'
+ kchg(kc,1) = -3
+ kchg(kc,2) = 0
+ kchg(kc,3) = 1
+ pmas(kc,1) = amprt(1)
+ qf = 1.
+ nc = 1.
+* Neutralino storing in JETSET
+ kc = pycomp (42)
+ chaf(kc,1) = 'chi_1'
+ kchg(kc,1) = 0
+ kchg(kc,2) = 0
+ kchg(kc,3) = 1
+ pmas(kc,1) = amprt(2)
+*
+* Gaugino mixing parameters
+ u11 = -zmixss(4,1)*sin(gammar)/sqrt(2.) +
+ & zmixss(2,1)*cos(gammar)
+ v11 = zmixss(3,1)*sin(gammal)/sqrt(2.) +
+ & zmixss(2,1)*cos(gammal)
+ DO idc = 1,nssmod
+ IF (issmod(idc) .EQ. 39) THEN
+ IF (moddcy .EQ. 1) THEN
+ DO jmod = 1,5
+ IF (jssmod(jmod,idc) .EQ. -12) THEN
+ xsbra = bssmod(idc)
+ xsbra = xsbra*xsbra*4
+ GOTO 10
+ END IF
+ END DO
+ ELSE
+ DO jmod = 1,5
+ IF (jssmod(jmod,idc) .EQ. -12) THEN ! e nu_e
+ ssbr(1) = bssmod(idc)
+ GOTO 20
+ ELSE IF (jssmod(jmod,idc) .EQ. -14) THEN ! mu nu_mu
+ ssbr(2) = bssmod(idc)
+ GOTO 20
+ ELSE IF (jssmod(jmod,idc) .EQ. -16) THEN ! tau nu_tau
+ ssbr(3) = bssmod(idc)
+ GOTO 20
+ ELSE IF (jssmod(jmod,idc) .EQ. -2) THEN ! u dbar
+ ssbr(4) = bssmod(idc)
+ GOTO 20
+ ELSE IF (jssmod(jmod,idc) .EQ. 4) THEN ! c sbar
+ ssbr(5) = bssmod(idc)
+ GOTO 20
+ END IF
+ END DO
+ 20 CONTINUE
+ END IF
+ END IF
+ END DO
+*
+ 10 CONTINUE
+ ELSE IF (iproc .EQ. 6) THEN
+* Fix parameters of process gg -> V1V2
+ IF (kv1 .GT. kv2) THEN
+ kvtmp = kv1
+ kv1 = kv2
+ kv2 = kvtmp
+ END IF
+ IF (kv1.EQ.113 .AND. kv2.EQ.113) THEN
+ idx_vv = 1
+ ELSE IF (kv1.EQ.113 .AND. kv2.EQ.223) THEN
+ idx_vv = 2
+ ELSE IF (kv1.EQ.113 .AND. kv2.EQ.333) THEN
+ idx_vv = 3
+ ELSE IF (kv1.EQ.113 .AND. kv2.EQ.443) THEN
+ idx_vv = 4
+ ELSE IF (kv1.EQ.223 .AND. kv2.EQ.223) THEN
+ idx_vv = 5
+ ELSE IF (kv1.EQ.223 .AND. kv2.EQ.333) THEN
+ idx_vv = 6
+ ELSE IF (kv1.EQ.223 .AND. kv2.EQ.443) THEN
+ idx_vv = 7
+ ELSE IF (kv1.EQ.333 .AND. kv2.EQ.333) THEN
+ idx_vv = 8
+ ELSE IF (kv1.EQ.333 .AND. kv2.EQ.443) THEN
+ idx_vv = 9
+ ELSE IF (kv1.EQ.443 .AND. kv2.EQ.443) THEN
+ idx_vv = 10
+ ELSE
+ WRITE (6,*) ' %GGINIT: unknown (V1 V2) combination'
+ WRITE (6,'(3x,''KV1 = '',i3,'', KV2 = '',i3)') kv1, kv2
+ STOP
+ END IF
+ DO 101 igvpar=1,4
+ 101 gvpar(igvpar) = gvconst(igvpar,idx_vv)
+ IF (gvpar(1) .EQ. 0.0) THEN
+ WRITE (6,*) ' %GGINIT: (V1 V2) combination not implemented'
+ WRITE (6,'(3x,''KV1 = '',i3,'', KV2 = '',i3)') kv1, kv2
+ STOP
+ END IF
+ END IF
+*
+ IF (iproc.EQ.1 .OR. iproc.EQ.3 .OR.
+ & iproc.EQ.4 .OR. iproc.EQ.5 .OR. iproc.EQ.6) THEN
+* Max cross section estimation
+ xsmax0 = 0.
+ DO nest = 1,100
+ xm = amin + ggrnd(0) * (amax-amin)
+ wsq = xm*xm
+ CALL ggxsec
+ END DO
+ xssum = 0.
+ ntry = 0
+ END IF
+*
+ RETURN
+ END
+*
+*=======================================================================
+CDECK ID>, GGWID.
+ FUNCTION ggwid (kf_res)
+*-----------------------------------------------------------------------
+* Routine to define 2-gamma width of a resonance
+* with JETSET code KF_RES
+* Author: Yu.Kharlov
+*-----------------------------------------------------------------------
+ REAL wid(5)
+ DATA wid /9.E-09, 1.E-06, 5.E-06, 6.3E-06, 0.41E-06/
+ IF (kf_res .EQ. 111) THEN
+ ggwid = wid(1)
+ ELSE IF (kf_res .EQ. 221) THEN
+ ggwid = wid(2)
+ ELSE IF (kf_res .EQ. 331) THEN
+ ggwid = wid(3)
+ ELSE IF (kf_res .EQ. 441
+ & .OR. kf_res .EQ. 10441
+ & .OR. kf_res .EQ. 445) THEN
+ ggwid = wid(4)
+ ELSE IF (kf_res .EQ. 551
+ & .OR. kf_res .EQ. 10551
+ & .OR. kf_res .EQ. 555) THEN
+ ggwid = wid(5)
+ END IF
+ RETURN
+ END
+*
+*=======================================================================
+CDECK ID>, GGRUN.
+ SUBROUTINE ggrun
+*-----------------------------------------------------------------------
+* Single event generation
+* Author: Yu.Kharlov
+*-----------------------------------------------------------------------
+ COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
+ & amin,amax,ymin,ymax,nmas,ny, kferm,
+ & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
+ & thetamin, costhv1, kv1,kv2,gvpar(4)
+ CHARACTER lumfil*80
+ COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
+ & gvconst(4,10)
+ REAL nc
+ COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
+ & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
+ COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
+ & xstote, ssbr(10)
+ COMMON /ggkin/ xkin(10)
+ COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
+ DOUBLE PRECISION P,V
+ SAVE /PYJETS/
+ COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
+ DOUBLE PRECISION PARU,PARJ
+ SAVE /PYDAT1/
+ COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
+ DOUBLE PRECISION PMAS,PARF,VCKM
+ SAVE /PYDAT2/
+ REAL pgam1(4),pgam2(4),pev(4,10)
+ REAL p3(4), p4(4)
+ REAL*8 dbetaz
+ REAL*8 ggrnd
+*
+ 1 CALL gg2gam (amas,wsq,ygg,q1sq,q2sq,
+ & pgam1,pgam2,ptag1,ptag2,weight)
+*
+* Cross section calculation
+ IF (iproc.EQ.1 .OR. iproc.EQ.3 .OR.
+ & iproc.EQ.4 .OR. iproc.EQ.5 .OR. iproc.EQ.6) THEN
+ CALL ggxsec
+ xssum = xssum + xscur
+ ntry = ntry + 1
+* Select event on cross section of process gg -> XX
+ IF (xscur0/xsmax0 .LT. sngl (ggrnd(0))) GOTO 1
+* Total cross section calculation due to current statistics
+ xstot = xssum/ntry
+ xstote = xstot / sqrt(float(ntry))
+ END IF
+*
+ xmgg = sqrt (wsq)
+ xmg1 = sqrt (q1sq)
+ xmg2 = sqrt (q2sq)
+ dbetaz = dtanh(1.d0*ygg)
+*
+* calculation of the 4-vector of two-photon system
+* and filling arrays pgg(1:2,1:5), kgg(1:2)
+ DO ip=1,4
+ p2g(ip) = pgam1(ip) + pgam2(ip)
+ pgg(1,ip) = pgam1(ip)
+ pgg(2,ip) = pgam2(ip)
+ END DO
+ p2g(5) = xmgg
+ pgg(1,5) = -xmg1
+ pgg(2,5) = -xmg2
+ kgg(1) = 22
+ kgg(2) = 22
+ ngg = 2
+*
+* =======> PROCESS # 1
+ IF (iproc .EQ. 1) THEN
+ DO ip = 1,4
+ p(1,ip) = dble (pgam1(ip))
+ p(2,ip) = dble (pgam2(ip))
+ END DO
+ p(1,5) = dble (-xmg1)
+ p(2,5) = dble (-xmg2)
+ CALL pyevnt
+*
+* =======> PROCESS # 2
+ ELSE IF (iproc .EQ. 2) THEN
+ ngg = ngg + 1
+ DO ip = 1,5
+ pgg(ngg,ip) = p2g(ip)
+ END DO
+ kgg(ngg) = kf_onium
+* Fill the JETSET common block
+ CALL gg2py (ngg)
+ CALL pyexec
+*
+* =======> PROCESS # 3, # 4, # 5, # 6
+ ELSE IF (iproc.EQ.3 .OR. iproc.EQ.4 .OR.
+ & iproc.EQ.5 .OR. iproc.EQ.6) THEN
+ phi = 2*pi*ggrnd(0)
+ IF (iproc.EQ.3 .OR. iproc.EQ.4 .OR. iproc.EQ.5) THEN
+ e3 = 0.5*sqrt(wsq)
+ e4 = e3
+ beta = xkin(1)
+ costhe = xkin(2)
+ am1 = xkin(3)
+ am2 = xkin(3)
+ p3abs = e3 * beta
+ ELSE IF (iproc.EQ.6) THEN
+ am1 = xkin(1)
+ am2 = xkin(2)
+ costhe = xkin(3)
+ p3abs = xkin(4)
+ e3 = xkin(5)
+ e4 = xmgg - e3
+ END IF
+ sinthe = sqrt (1 - costhe*costhe)
+ p3(1) = p3abs * cos(phi) * sinthe
+ p3(2) = p3abs * sin(phi) * sinthe
+ p3(3) = p3abs * costhe
+ p3(4) = e3
+ p4(4) = e4
+ DO ip = 1,3
+ p4(ip) = -p3(ip)
+ END DO
+ igg1 = ngg + 1
+ igg2 = ngg + 2
+ DO ip = 1,4
+ pgg(igg1,ip) = p3(ip)
+ pgg(igg2,ip) = p4(ip)
+ END DO
+ pgg(igg1,5) = am1
+ pgg(igg2,5) = am2
+ ngg = ngg + 2
+*
+ IF (iproc .EQ. 3) THEN
+ kgg(igg1) = kferm
+ kgg(igg2) = -kferm
+ CALL gg2py (ngg)
+ ELSE IF (iproc .EQ. 4) THEN
+ kgg(igg1) = 24
+ kgg(igg2) = -24
+ CALL gg2py (ngg)
+ ELSE IF (iproc .EQ. 5) THEN
+ kgg(igg1) = 41
+ kgg(igg2) = -41
+* Decay of charginos into W+neutralino
+ CALL ggdecy (igg1)
+ CALL ggdecy (igg2)
+*
+ CALL gg2py (ngg)
+ k(igg1,1) = 11
+ k(igg2,1) = 11
+ ELSE IF (iproc.EQ.6) THEN
+ kgg(igg1) = kv1
+ kgg(igg2) = kv2
+ CALL gg2py (ngg)
+ END IF
+*
+* Boost event into initial frame
+ mstu(33) = 1
+ CALL pyrobo(3,ngg, 0.d0,0.d0, 0.d0,0.d0, dbetaz)
+ CALL pyexec
+*
+ END IF
+*
+ RETURN
+ END
+*
+*=======================================================================
+CDECK ID>, GGEXIT.
+ SUBROUTINE ggexit
+*-----------------------------------------------------------------------
+* Cross section calculation and end of run and
+* print out cross section and related variables
+* Author: Yu.Kharlov
+*-----------------------------------------------------------------------
+ COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
+ & amin,amax,ymin,ymax,nmas,ny, kferm,
+ & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
+ & thetamin, costhv1, kv1,kv2,gvpar(4)
+ CHARACTER lumfil*80
+ COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
+ & xstote, ssbr(10)
+ COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
+ & gvconst(4,10)
+ REAL nc
+*
+ WRITE (6,'(/x,79(''=''))')
+ WRITE (6,'(x,''I'',36x,''TPHIC'',36x, ''I'')')
+ WRITE (6,'(x,''I'',34x,''Process '',i1,34x,''I'')') iproc
+ WRITE (6,'(x,''I'',34x,''---------'', 34x,''I'')')
+ WRITE (6,'(x,''I'',77x, ''I'')')
+ WRITE (6,'(x,''I'',20x,''Events thrown : '',i8,31x,''I'')')
+ & ntry
+ WRITE (6,'(x,''I'',20x,''Events accepted : '',i8,31x,''I'')')
+ & nevent
+ WRITE (6,'(x,''I'',20x,''Cross section : '',e10.3,'' +- '','//
+ & 'e10.3,'' nb'',12x ''I'')') xstot, xstote
+ WRITE (6,'(x,''I'',77x, ''I'')')
+ WRITE (6,'(x,79(''=''))')
+*
+ IF (iproc .EQ. 1) THEN
+* --- Alternative cross section for minimum bias events ---
+ totxsc = 0.0
+ hmas = (amax-amin)/nmas
+ f1 = dldmx(amin ,ny) * dchad(amin)
+ DO i=1,nmas
+ xm = amin + i*hmas
+ f2 = dldmx(xm-hmas/2.,ny) * dchad(xm-hmas/2.)
+ f3 = dldmx(xm ,ny) * dchad(xm)
+ totxsc = totxsc +(f1 + 4.*f2 + f3)/6.
+ f1 = f3
+ END DO
+ sigtt = totxsc * hmas
+ sigmx = sigtt/sqrt(float(nevent))
+ WRITE (6,*) ' CrosSec =',sigtt,' +/-',sigmx,' mb'
+ END IF
+*
+ RETURN
+ END
+*
+ REAL FUNCTION dchad(xm)
+*-----------------------------------------------------------------------
+* Total cross section for gamma-gamma collisions vs. Xm in mb
+* Parametrization is given by G.A.Schuler, T.Sjostrand,
+* CERN-TH.7193/94, formula (9).
+* Coded by Yu.Kharlov 20.03.1995.
+*-----------------------------------------------------------------------
+ DATA epsilon, eta /0.0808, 0.4525/
+ dchad = 0.211e-3 * (xm**2)**epsilon + 0.297e-3/(xm**2)**eta
+ RETURN
+ END
+*
+*=======================================================================
+CDECK ID>, GGXSEC.
+ SUBROUTINE ggxsec
+*-----------------------------------------------------------------------
+* Cross section calculation in current event for iproc=1,3,4,5,6
+* Author: Yu.Kharlov
+*-----------------------------------------------------------------------
+ COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
+ & amin,amax,ymin,ymax,nmas,ny, kferm,
+ & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
+ & thetamin, costhv1, kv1,kv2,gvpar(4)
+ CHARACTER lumfil*80
+ COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
+ & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
+ COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
+ & gvconst(4,10)
+ REAL nc
+ COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
+ & xstote, ssbr(10)
+ COMMON /ggkin/ xkin(10)
+ COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
+ DOUBLE PRECISION PMAS,PARF,VCKM
+ SAVE /PYDAT2/
+ PARAMETER (epsilon = 0.0808, eta = -0.4525)
+ REAL*8 ggrnd
+ REAL*8 pymass
+*
+ IF (iproc .EQ. 1) THEN
+* Parametrization by G.A.Schuler, T.Sjostrand,
+* CERN-TH.7193/94, formula (9).
+ xscur0 = 211. * wsq**epsilon + 297. * wsq**eta
+ xscur = xscur0 * xlumint
+ ELSE IF (iproc .EQ. 3 .OR. iproc .EQ. 5) THEN
+ beta2 = 1.0 - 4.0*amprt(1)**2/wsq
+ beta2 = min (beta2, 1.0)
+ beta = sqrt (beta2)
+ IF (beta2 .LT. 0.9) THEN
+ 1 costhe = 2.*ggrnd(0) - 1.
+ costh2 = costhe*costhe
+ xnum = 1. + 2.*beta2 * (1.-beta2) * (1.-costh2) -
+ & (beta2*costh2)**2
+ xden = (1 - beta2*costh2)**2
+ ffcur = xnum / xden
+ ffmax = (1. + beta2) / (1. - beta2)
+ rff = ggrnd(0)
+ IF (ffcur .LE. rff*ffmax) GOTO 1
+ ELSE
+ small = 1.0 - cos(thetamin)
+ 2 zz = exp (ggrnd(0) * alog(2./small))
+C sgn = sign (1., 2.*ggrnd(0)-1.)
+ sgn = 2.*ggrnd(0)-1.
+ sgn = sign (1., sgn)
+ costhe = sgn * (zz-1.) / (zz+1.)
+ ffcur = (1.+costhe*costhe) / (1.-beta2*costhe*costhe)
+ ffmax = 2. / (1.- costhe*costhe)
+ rff = ggrnd(0)
+ IF (ffcur .LE. rff*ffmax) GOTO 2
+ END IF
+*
+ xkin(1) = beta
+ xkin(2) = costhe
+ xkin(3) = amprt(1)
+*
+* Yu.Kharlov, 3.06.1995
+ totfac = 4. * pi * alpha**2 * gev2nb * qf**4 * nc / wsq
+ IF (beta2 .LT. 0.99) THEN
+ xscur0 =((3.-beta2**2)/(2.*beta) * alog ((1.+beta)/(1.-beta))-
+ & 2. + beta2) * beta
+ ELSE
+ xscur0 = 2. * alog (sqrt(wsq)/amprt(1)) - 1.0
+ END IF
+ xscur0 = totfac * xscur0 * xsbra
+ xscur = xscur0 * xlumint
+ ELSE IF (iproc .EQ. 4) THEN
+ xmw = amprt(4)
+ beta2 = 1. - 4.*xmw**2/wsq
+ beta = sqrt (beta2)
+ 3 costhe = 2.*ggrnd(0) - 1.
+ costh2 = costhe*costhe
+ xnum = 19. - 6.*beta2 * (1.-beta2) +
+ & 2.*beta2 * (8. - 3.*beta2) * costh2 +
+ & 3.*(beta2*costh2)**2
+ xden = (1. - beta2*costh2)**2
+ ffcur = xnum / xden
+ ffmax = (19. + 10.*beta2 + beta2**2) / (1. - beta2)**2
+ rff = ggrnd(0)
+ IF (ffcur .LE. rff*ffmax) GOTO 3
+*
+ xkin(1) = beta
+ xkin(2) = costhe
+ xkin(3) = xmw
+*
+* Yu.Kharlov, 3.06.1995
+ totfac = pi * alpha**2 * beta * gev2nb/ wsq
+ xscur0 = 2. * (22. - 9.*beta2 + 3.*beta2**2) / (1. - beta2) -
+ & 3. * (1.-beta2**2)/beta * alog ((1.+beta)/(1.-beta))
+ xscur0 = totfac * xscur0 * xsbra
+ xscur = xscur0 * xlumint
+ ELSE IF (iproc .EQ. 6) THEN
+ xmgg = sqrt (wsq)
+ bb = gvpar(3) + gvpar(4)*alog(xmgg/50.0)
+ am1 = sngl (pymass(kv1))
+ am2 = sngl (pymass(kv2))
+ p1 = pfork (xmgg,am1,am2)
+ e1 = sqrt (p1*p1 + am1*am1)
+ tmin = am1*am1 - xmgg*(e1+p1*costhv1)
+ tmax = am1*am1 - xmgg*(e1-p1*costhv1)
+ fmin = exp (bb*tmin)
+ fmax = exp (bb*tmax)
+ IF (fmax .GT. 1.e-10) THEN
+ tvar = 1./bb * dlog (fmax - ggrnd(0)*(fmax-fmin))
+ ELSE
+ tvar = tmax + 1./bb * dlog (1. - ggrnd(0))
+ END IF
+ costhe = (e1*xmgg - am1*am1 + tvar) / xmgg / p1
+ IF (costhe .GT. 1.0) costhe = 1.0
+ IF (costhe .LT.-1.0) costhe = -1.0
+ IF (ggrnd(0) .LT. 0.5) costhe = -costhe
+*
+ xkin(1) = am1
+ xkin(2) = am2
+ xkin(3) = costhe
+ xkin(4) = p1
+ xkin(5) = e1
+*
+ IF (kv1 .EQ. 443 .AND. kv2 .EQ. 443) THEN
+ factor = 1.0 / alog(xmgg/3.097)
+ ELSE
+ factor = 1.0
+ END IF
+ xscur0 = factor * gvpar(1) * xmgg**gvpar(2) / bb * exp(bb*tmax)
+ xscur = xscur0 * xlumint
+*
+ END IF
+*
+* Find max cross section of gg-process
+*
+ IF (xscur0 .GT. xsmax0) xsmax0 = xscur0
+*
+ RETURN
+ END
+*
+*=======================================================================
+CDECK ID>, GGDECY.
+ SUBROUTINE ggdecy (igg)
+*-----------------------------------------------------------------------
+* Routine to decay a particle number I in GGEVNT common block
+* Input : IGG - particle number in common GGEVNT arrays
+* Author : Yu.Kharlov
+*-----------------------------------------------------------------------
+ COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
+ & amin,amax,ymin,ymax,nmas,ny, kferm,
+ & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
+ & thetamin, costhv1, kv1,kv2,gvpar(4)
+ CHARACTER lumfil*80
+ COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
+ & gvconst(4,10)
+ REAL nc
+ COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
+ & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
+ COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
+ & xstote, ssbr(10)
+ COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
+ DOUBLE PRECISION PMAS,PARF,VCKM
+ SAVE /PYDAT2/
+ COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr,
+ & xmll, xmlr, xmnl, xtanb, xmha, xmu,
+ & xmt, xat, xmbr, xab, u11, v11
+ REAL p1(4), pcm(4), p2(4), p3(4), p4(4), ptmp(3)
+ REAL o1(3,3), o2(3,3)
+ INTEGER pycomp
+ DATA ifl /0/
+ REAL*8 ggrnd
+*
+ DO ip = 1,4
+ p1(ip) = pgg(igg,ip)
+ END DO
+*
+ IF (abs(kgg(igg)) .EQ. 41 .AND. iproc .EQ. 3) THEN
+* Decay chi_1^+- --> chi_1^0 W^+- in CMS, uniformly
+ eb = (amprt(1)**2 + amprt(2)**2 - amprt(4)**2) / (2 * amprt(1))
+ ec = (amprt(1)**2 - amprt(2)**2 + amprt(4)**2) / (2 * amprt(1))
+ pb = sqrt (eb*eb - amprt(2)*amprt(2))
+ pc = pb
+ costhe = 2*ggrnd(0) - 1.
+ sinthe = sqrt (1 - costhe*costhe)
+ phi = 2*pi*ggrnd(0)
+ pcm(1) = pb * sinthe * cos(phi)
+ pcm(2) = pb * sinthe * sin(phi)
+ pcm(3) = pb * costhe
+ pcm(4) = eb
+ CALL lorenb (amprt(1),p1,pcm,p2)
+ DO ip = 1,3
+ pcm(ip) = -pcm(ip)
+ END DO
+ pcm(ip) = ec
+ CALL lorenb (amprt(1),p1,pcm,p3)
+ DO ip = 1,4
+ pgg(ngg+1,ip) = p2(ip)
+ pgg(ngg+2,ip) = p3(ip)
+ END DO
+ pgg(ngg+1,5) = amprt(2)
+ pgg(ngg+2,5) = amprt(4)
+ kgg(ngg+1) = 42 * sign(1.,1.*kgg(igg))
+ kgg(ngg+2) = 24 * sign(1.,1.*kgg(igg))
+ ngg = ngg + 2
+*
+ ELSE IF (abs(kgg(igg)) .EQ. 41 .AND. iproc .EQ. 5) THEN
+* Decay chi_1^+- --> e^+- nu_e chi_1^0 in CMS,
+* according to the exact matrix element
+*
+ xm23max = amprt(1) - amprt(2)
+ xmw = sngl (pmas(pycomp(24),1))
+ gwt = sngl (pmas(pycomp(24),2))
+ atmin = - atan (xmw/gwt)
+ atmax = atan ((xm23max**2 - xmw**2) / (xmw*gwt))
+ am1 = amprt(1)
+ am4 = amprt(2)
+*
+ 1 CONTINUE
+*
+* Random mass of (2+3) according to Breit-Wigner
+*
+ xm23 = xmw**2 + xmw*gwt *
+ * tan (atmin + ggrnd(0)*(atmax - atmin))
+ xm23 = max (0., xm23)
+ xm23 = sqrt (xm23)
+ IF (xm23 .GT. xm23max) GOTO 1 ! to avoid boundary effect
+ e23 = (am1**2 + xm23**2 - am4**2) / 2. / am1
+*
+* Rest frame of (2+3)
+*
+ phi = 2.*pi * ggrnd(0)
+ costhe = 2.*ggrnd(0) - 1.
+ sinthe = sqrt (1. - costhe*costhe)
+*
+ pp = xm23 / 2.
+ p2(1) = pp * sinthe * cos(phi)
+ p2(2) = pp * sinthe * sin(phi)
+ p2(3) = pp * costhe
+ p2(4) = pp
+*
+ DO ip = 1,3
+ p3(ip) = -p2(ip)
+ END DO
+ p3(4) = pp
+*
+* Transformation from rest frame of (2+3) to rest frame of (1),
+* 3-momenta p2, p3 are along z-axis
+*
+ gamma = e23 / xm23
+ gambe = sqrt (gamma*gamma - 1.)
+*
+ p2z = gamma*p2(3) + gambe*p2(4)
+ e2 = gamma*p2(4) + gambe*p2(3)
+ p3z = gamma*p3(3) + gambe*p3(4)
+ e3 = gamma*p3(4) + gambe*p3(3)
+*
+ p2(3) = p2z
+ p2(4) = e2
+ p3(3) = p3z
+ p3(4) = e3
+*
+* Arbitrary rotation of rest frame of (1)
+*
+ phi = 2.*pi * ggrnd(0)
+ costhe = 2.*ggrnd(0) - 1.
+ sinthe = sqrt (1. - costhe*costhe)
+ cosphi = cos(phi)
+ sinphi = sin(phi)
+*
+ o1(1,1) = costhe
+ o1(1,2) = 0.
+ o1(1,3) = -sinthe
+ o1(2,1) = 0.
+ o1(2,2) = 1.
+ o1(2,3) = 0.
+ o1(3,1) = sinthe
+ o1(3,2) = 0.
+ o1(3,3) = costhe
+*
+ o2(1,1) = cosphi
+ o2(1,2) = -sinphi
+ o2(1,3) = 0.
+ o2(2,1) = sinphi
+ o2(2,2) = cosphi
+ o2(2,3) = 0.
+ o2(3,1) = 0.
+ o2(3,2) = 0.
+ o2(3,3) = 1.
+*
+ DO ix = 1,3
+ ptmp(ix) = 0.
+ DO jx = 1,3
+ ptmp(ix) = ptmp(ix) + o1(ix,jx) * p2(jx)
+ END DO
+ END DO
+ DO ix = 1,3
+ p2(ix) = ptmp(ix)
+ END DO
+ DO ix = 1,3
+ ptmp(ix) = 0.
+ DO jx = 1,3
+ ptmp(ix) = ptmp(ix) + o2(ix,jx) * p2(jx)
+ END DO
+ END DO
+ DO ix = 1,3
+ p2(ix) = ptmp(ix)
+ END DO
+*
+ DO ix = 1,3
+ ptmp(ix) = 0.
+ DO jx = 1,3
+ ptmp(ix) = ptmp(ix) + o1(ix,jx) * p3(jx)
+ END DO
+ END DO
+ DO ix = 1,3
+ p3(ix) = ptmp(ix)
+ END DO
+ DO ix = 1,3
+ ptmp(ix) = 0.
+ DO jx = 1,3
+ ptmp(ix) = ptmp(ix) + o2(ix,jx) * p3(jx)
+ END DO
+ END DO
+ DO ix = 1,3
+ p3(ix) = ptmp(ix)
+ END DO
+*
+ p12 = am1 * p2(4)
+*
+ e4 = am1 - e23
+ pp =sqrt (e4*e4 - am4*am4)
+ p4(1) = pp * sinthe * cos(phi)
+ p4(2) = pp * sinthe * sin(phi)
+ p4(3) = -pp * costhe
+ p4(4) = e4
+*
+ width = (- 4.*p12**2*u11**2 - 4.*p12**2*v11**2 +
+ & 4.*p12*xm23**2*u11**2 - 2.*p12*am4**2*u11**2 -
+ & 2.*p12*am4**2*v11**2 + 2.*p12*am1**2*u11**2 +
+ & 2.*p12*am1**2*v11**2 - xm23**4*u11**2 +
+ & xm23**2*am4**2*u11**2 - 2.*xm23**2*am4*am1*u11*v11 -
+ & xm23**2*am1**2*u11**2) /
+ & ((xm23**2-xmw**2)**2 + (xmw*gwt)**2)
+ width = width * pp * xm23
+*
+ widmax = (2.*am1*xm23max *
+ & (2.*(xm23max*u11)**2 + am1**2*(u11**2 + v11**2)) +
+ & (xm23max*am4*u11)**2) * (am1**2 - am4**2) / 2./am1 *
+ & xm23 / ((xm23**2-xmw**2)**2 + (xmw*gwt)**2)
+*
+* Select current phase space configuration
+*
+ IF (width .GT. widmax) WRITE (6,*)
+ & '%GGDECY: Maximum chi+ width violated'
+ rwid = widmax * ggrnd(0)
+ IF (width .LE. rwid) GOTO 1
+*
+* Pick up final lepton flavour
+ rflav = ggrnd(0)
+ IF (moddcy .EQ. 1) THEN
+ IF (rflav .LE. 0.5) THEN
+ lflav1 = 11
+ lflav2 = -12
+ ELSE
+ lflav1 = 13
+ lflav2 = -14
+ END IF
+ ELSE
+ rr1 = ssbr(1)
+ rr2 = ssbr(1)+ssbr(2)
+ rr3 = ssbr(1)+ssbr(2)+ssbr(3)
+ rr4 = ssbr(1)+ssbr(2)+ssbr(3)+ssbr(4)
+ rr5 = ssbr(1)+ssbr(2)+ssbr(3)+ssbr(4)+ssbr(5)
+ IF (rflav .LE. rr1) THEN
+ lflav1 = 11
+ lflav2 = -12
+ ELSE IF (rflav .LE. rr2) THEN
+ lflav1 = 13
+ lflav2 = -14
+ ELSE IF (rflav .LE. rr3) THEN
+ lflav1 = 15
+ lflav2 = -16
+ ELSE IF (rflav .LE. rr4) THEN
+ lflav1 = -2
+ lflav2 = 1
+ ELSE IF (rflav .LE. rr5) THEN
+ lflav1 = -4
+ lflav2 = 3
+ ELSE
+ GOTO 1
+ END IF
+ END IF
+*
+* Transformation to CMS of 2-gamma
+ CALL lorenb (amprt(1),p1,p2,pcm)
+ DO ip = 1,4
+ pgg(ngg+1,ip) = pcm(ip)
+ END DO
+ pgg(ngg+1,5) = 0.
+ kgg(ngg+1) = lflav1 * sign(1.,1.*kgg(igg))
+ ngg = ngg + 1
+*
+ CALL lorenb (amprt(1),p1,p3,pcm)
+ DO ip = 1,4
+ pgg(ngg+1,ip) = pcm(ip)
+ END DO
+ pgg(ngg+1,5) = 0.
+ kgg(ngg+1) = lflav2 * sign(1.,1.*kgg(igg))
+ ngg = ngg + 1
+*
+ CALL lorenb (amprt(1),p1,p4,pcm)
+ DO ip = 1,4
+ pgg(ngg+1,ip) = pcm(ip)
+ END DO
+ pgg(ngg+1,5) = amprt(2)
+ kgg(ngg+1) = 42 * sign(1.,1.*kgg(igg))
+ ngg = ngg + 1
+ END IF
+*
+ RETURN
+ END
+*
+*=======================================================================
+CDECK ID>, GG2PY.
+ SUBROUTINE gg2py (nngg)
+*-----------------------------------------------------------------------
+* Routine to fill PYTHIA event record by particles from GG
+* Input: NNGG : number of particles
+* Author: Yu.Kharlov
+*-----------------------------------------------------------------------
+ COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5),
+ & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5)
+ COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
+ DOUBLE PRECISION P,V
+ SAVE /PYJETS/
+ DO igg = 1,nngg
+ IF (igg .LE. 2) THEN
+ k(igg,1) = 21
+ ELSE
+ k(igg,1)=1
+ END IF
+ k(igg,2) = kgg(igg)
+ DO ip = 1,5
+ p(igg,ip) = pgg(igg,ip)
+ END DO
+ END DO
+ n = nngg
+ RETURN
+ END
+*
+*=======================================================================
+CDECK ID>, PFORK.
+ FUNCTION pfork (am0,am1,am2)
+*-----------------------------------------------------------------------
+* Routine to find a momentum of a secondary particle in the decay
+* of a particle with mass AM0 in the rest to 2 particles with
+* masses AM1 and AM2.
+*-----------------------------------------------------------------------
+ denom = am0**4 + am1**4 + am2**4 -
+ & 2.*(am0*am1)**2 - 2.*(am0*am2)**2 - 2.*(am1*am2)**2
+ IF (denom .GT. 0.) THEN
+ pfork = sqrt (denom) / 2./am0
+ ELSE
+ pfork = -1.
+ END IF
+ RETURN
+ END
+*=======================================================================
+CDECK ID>, GG2GAM.
+ SUBROUTINE gg2gam (amas,wsq,y3,qs1,qs2,pgam1,pgam2,ptag1,ptag2,wt)
+C***********************************************************************
+C
+C Generator of the Direct Two Photon Processes in S.A.Sadovsky
+C coherent heavy ion collisions: A+B --> 1+2+3 17.01.1995
+C with double differencial gg-luminosity function Yu.Kharlov
+C 20.03.1995
+C
+C LUMFUN must be called first to initialise the generator
+C
+C***********************************************************************
+C
+C DETAILED DESCRIPTION
+C ======== ===========
+C
+C INPUTS: only through common /TWOGENC/
+C =======
+C
+C OUTPUTS: AMAS The mass of ion.
+C ======== WSQ Square mass of the gamma-gamma system
+C Y3 The rapidity of the gamma-gamma system
+C QS1 Virtual mass squared of photon 1.
+C QS2 Virtual mass squared of photon 2.
+C PGAM1 Four-vector of photon 1
+C PGAM2 Four-vector of photon 2
+C PTAG1 Four-vector of tag 1
+C PTAG2 Four-vector of tag 2
+C WT The weight of the event (dummy so far)
+C
+C The remaining variables are internal to the program.
+C
+C***********************************************************************
+ INTEGER icnter,ncntrs,nevts,lout,debug
+ PARAMETER (ncntrs =7)
+ COMMON /twgenc/s,eb,wmin,wmax,th1pmn,th1pmx,th2pmn,th2pmx,
+ + st12mn,st12mx,st12rt,st22mn,st22mx,st22rt,ommin,omrat,eps,
+ + wghtsm,wghtgt,wghtmx,wtmxfd,xmres,xgres,
+ + icnter(ncntrs),nevts,lout,lres,lwate,debug
+ SAVE /twgenc/
+ REAL st12mn,st12mx,st12rt,st22mn,st22mx,st22rt
+ REAL wmin,wmax,ommin,omrat,eps,th1pmn,th1pmx
+ REAL th2pmn,th2pmx,s,eb,wghtsm,wghtmx,wghtgt,wtmxfd
+ REAL xmres,xgres
+ LOGICAL lres,lwate
+C
+ REAL alpha,pi,wlum,weight
+ PARAMETER (alpha=1./137.036,pi=3.14159265)
+C
+ REAL ptag1(10),ptag2(10),pgam1(10),pgam2(10),wt
+ REAL amas,wsq,q1sq,q2sq
+C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ REAL*4 xmin,xmax,ymin,ymax,gmma,y3,xm,q0,qs,qs1,qs2
+ REAL*8 pa(5),pb(5),p1(5),p2(5),p3(5)
+C
+ REAL*8 ph1,ph2,ak1,ak2,qk1,qk2,q1,q2,q3,e1,e2,e3
+ REAL*8 y,y1,y2,phi,pi2,pae,aeq,beq,wps,tm3,dts,wpt
+C
+ LOGICAL lstr
+ DATA q0,pi2,lstr / 0.060, 6.28318530718 d0, .false. /
+ REAL*8 ggrnd
+C
+ 999 CONTINUE
+ icnter(1) = icnter(1) + 1
+C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ IF (lstr) go to 10
+ lstr =.true.
+ qs = q0/sqrt(2.)
+C
+ pa(1)= 0D0
+ pa(2)= 0D0
+ pa(3)= dsqrt(1D0*(eb*eb-amas*amas))
+ pa(4)= eb*1D0
+ pa(5)= amas*1D0
+C
+ pb(1)= 0D0
+ pb(2)= 0D0
+ pb(3)=-pa(3)
+ pb(4)= pa(4)
+ pb(5)= amas*1D0
+C
+ 10 xm = xmres
+* Generate 2gamma rapidity y3 and mass xm
+ CALL ranlum(y3,xm,wt)
+C
+ wsq = xm*xm
+C
+* Energy and z-momentum of 2gamma system
+ egg = xm*cosh(y3)
+ pgg = xm*sinh(y3)
+*
+* Squared masses of the photons
+ qs1 = - qs**2 * dlog(ggrnd(0)) ! is it true that squared mass
+ qs2 = - qs**2 * dlog(ggrnd(0)) ! is distributed in exponent?
+ xm1 = -qs1
+ xm2 = -qs2
+*
+* Energies and z-momenta of the photons
+ xfun12 = wsq + xm1 - xm2
+*
+ eg1 = (egg + pgg*sqrt(1 - 4*wsq*xm1/(xfun12*xfun12))) *
+ & xfun12/(2*wsq)
+ eg2 = egg - eg1
+ pg1 = sqrt (eg1**2 - xm1)
+ pg2 = -sqrt (eg2**2 - xm2)
+*
+* 4-momenta of the photons
+ DO ip=1,2
+ pgam1(ip) = 0.
+ pgam2(ip) = 0.
+ END DO
+ pgam1(3) = pg1
+ pgam1(4) = eg1
+ pgam2(3) = pg2
+ pgam2(4) = eg2
+*
+* 4-momenta of the recoil ions
+ DO ip=1,4
+ ptag1(ip) = pa(ip) - pgam1(ip)
+ ptag2(ip) = pb(ip) - pgam2(ip)
+ END DO
+*
+ RETURN
+ END
+*
+*=======================================================================
+CDECK ID>, RANLUM.
+ SUBROUTINE ranlum (ry,rm,wlum)
+C G.V.Khaoustov
+C 12.12.1994
+C Corrected: Yu.Kharlov
+C 23.05.1995
+C 24-JUN-1997
+C
+C generate two random numbers (ry and rm) with probability proportional
+C to luminosity function in ymin-ymax rapidity and xmin-xmax mass range
+C
+ COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
+ & gvconst(4,10)
+ REAL nc
+ PARAMETER (isize=100000)
+ COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
+ + xmin,xmax,amn,amx,nnms,alfun(isize)
+ PARAMETER (NCNTRS =7)
+ COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX,
+ + ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS,
+ + WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES,
+ + ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG
+ LOGICAL lres,lwate
+ REAL*4 wlum
+ DATA ifl/0/, stpy/0./, stpm/0./
+ REAL*8 ggrnd
+C
+ IF (ifl.EQ.0) THEN
+ IF (ymin.LT.ymn .OR. ymin.GT.ymx .OR.
+ * ymax.LT.ymn .OR. ymax.GT.ymx .OR.
+ * xmin.LT.amn .OR. xmin.GT.amx .OR.
+ * xmax.LT.amn .OR. xmax.GT.amx) THEN
+ WRITE (6,*) 'Luminosity function is not defined in this range'
+ STOP
+ ENDIF
+ aky = ymax-ymin
+ akm = xmax-xmin
+ stpy = (ymx-ymn)/nny
+ stpm = (amx-amn)/nnms
+ rmin = exp(cexp*(xmax))
+ rmax = exp(cexp*(xmin))
+ dr = rmax-rmin
+ ifl = 1
+ ENDIF
+C
+ 1 CONTINUE
+ IF (.NOT.lres) THEN
+ rm = ggrnd(0) ! random mass
+ rm = rmin + rm*dr
+ rm = alog(rm)/cexp ! exponetial law
+ ELSE
+ wsq = xmres**2 + xmres*xgres * tan(pi*(ggrnd(0)-0.5))
+ IF (wsq .LT. xmin*xmin) GOTO 1
+ IF (wsq .GT. xmax*xmax) GOTO 1
+ rm = sqrt (wsq)
+ END IF
+C
+ ry = ggrnd(0) ! random Y
+*
+* Find y-limits for generated gg-mass: y_max = log(E_max/M_gg)
+*
+ ym_max = alog (egg_max/rm)
+ ygg_max = min (ymax, ym_max)
+ ygg_min = max (ymin,-ym_max)
+ ry = ygg_min + ry*(ygg_max-ygg_min)
+ iy = (ry-ymn)/stpy+1.
+C
+ im = (rm-amn)/stpm
+ r = anorm*exp(cexp*rm)*ggrnd(0) ! normalization
+ p = (ry-(iy-1)*stpy-ymn)/stpy ! four point interpolation
+ q = (rm-im*stpm-amn)/stpm
+ k = iy+im*(nny+1)
+ f = (1.-p)*(1.-q)*alfun(k) + p*(1.-q)*alfun(k+1) +
+ * q*(1.-p)*alfun(k+nny+1) + q*p*alfun(k+2+nny)
+ IF (r .GT. f) GOTO 1
+ wlum = 1000.*f ! differential luminosity in GeV^-1
+*
+ RETURN
+ END
+*
+*=======================================================================
+CDECK ID>, GGDATA.
+ BLOCK DATA ggdata
+*-----------------------------------------------------------------------
+* Some useful constants and masses of MSSM particles
+* Author: Yu.Kharlov
+*-----------------------------------------------------------------------
+ COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
+ & gvconst(4,10)
+ REAL nc
+ COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas,
+ & amin,amax,ymin,ymax,nmas,ny, kferm,
+ & kf_onium,xmres,xgtres,xggres, xlumint, moddcy,
+ & thetamin, costhv1, kv1,kv2,gvpar(4)
+ CHARACTER lumfil*80
+ COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot,
+ & xstote, ssbr(10)
+ COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr,
+ & xmll, xmlr, xmnl, xtanb, xmha, xmu,
+ & xmt, xat, xmbr, xab, u11, v11
+ COMMON/D2LParam/lz,la,Rion,Gamma
+ INTEGER lz,la
+ REAL Rion,Gamma
+ PARAMETER (al = 1./128.)
+ DATA pi /3.14159276/, alpha /al/, gev2nb /389379./,
+ & hbarc /0.197327/
+ DATA iproc, nevent, ilumf /1, 10000, -1/
+ DATA lumfil /'gglum.dat'/
+ DATA Rion /0.0/
+ DATA amprt /110.,10.,150.,-1.,0.106/
+ DATA thetamin /1.e-03/
+ DATA costhv1 /1.0/
+ DATA xsbra /1./
+ DATA ny, nmas /50, 50/
+ DATA xm1, xm2, xmg, xms, xmtl, xmtr,
+ & xmll, xmlr, xmnl, xtanb, xmha, xmu,
+ & xmt, xat, xmbr, xab
+ & /30., 50., 250., 700., 600., 600.,
+ & 400., 400., 400., 4., 300., -300.,
+ & 174., 300., 300., 300./
+ DATA gvconst / 63.0, 0.22, 6.0, 1.0, ! rho rho
+ & 14.0, 0.22, 6.0, 1.0, ! rho omega
+ & 7.0, 0.22, 4.0, 1.0, ! rho phi
+ & 0.05, 0.80, 2.5, 0.0, ! rho J/psi
+ & 0.0, 0.00, 0.0, 0.0, ! omega omega
+ & 0.8, 0.22, 4.0, 1.0, ! omega phi
+ & 6.E-3, 0.80, 2.5, 0.0, ! omega J/psi
+ & 0.2, 0.22, 2.0, 1.0, ! phi phi
+ & 3.E-3, 0.80, 1.5, 0.0, ! phi J/psi
+ & 6.E-5, 2.00, 1.5, 0.0/ ! J/psi J/psi
+ END
+*
+*=======================================================================
+CDECK ID>, EVPRINT.
+ SUBROUTINE evprint(Lun,Iev,P3,Ygg,ifl,Eb,ia,iz,amas,ptag1,ptag2)
+*-----------------------------------------------------------------------
+* Event output for internal use by GGHIC authors
+*-----------------------------------------------------------------------
+ COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
+ DOUBLE PRECISION P,V
+ SAVE /PYJETS/
+ REAL ptag1(4),ptag2(4),p3(5)
+C
+ WRITE(lun,100) iev,n,Eb,p3(5),ygg,ifl,-ifl
+ WRITE(lun,101) 90000+ia,iz,ptag1,amas
+ WRITE(lun,101) 90000+ia,iz,ptag2,amas
+ WRITE(lun,101) 90000+ 3, 0,(p3(j),j=1,5)
+ IF (n.NE.0) THEN
+ DO i=1,n
+ ichg=pychge(k(i,2))/3
+ WRITE(lun,101) k(i,2),ichg,(p(i,j),j=1,5)
+ ENDDO
+ ENDIF
+ RETURN
+ 100 FORMAT(i6,i3,3e16.8,2x,2i8)
+ 101 FORMAT(i6,i3,3(1x,e14.8),e13.8,e13.8)
+ END
+*
+*=======================================================================
+C
+CDECK ID>, TWOINI.
+ SUBROUTINE TWOINI(LOUTX,IPAR,XPAR,YPAR) TWOG0965
+C***********************************************************************TWOG0966
+C TWOG0967
+C Initialises the gamma-gamma generator of Ion collisions TWOGAM. TWOG0968
+C TWOG0969
+C INPUTS: LOUTX Logical unit for print-out TWOG0970
+C ======= XPAR(1) EBM Beam energy GeV TWOG0971
+C XPAR(2) WMN Minimum two-gamma mass GeV TWOG0972
+C XPAR(3) WMX Maximum two-gamma mass GeV TWOG0973
+C TWOG0974
+C XPAR(8) WGHTMX Maximum weight TWOG0978
+C XPAR(9) XMRES Mass of resonance GeV TWOG0979
+C XPAR(10) XGRES Total width of resonance GeV TWOG0980
+C TWOG0981
+C IPAR(1) 0 Unweighted events TWOG0982
+C 1 Weighted events TWOG0983
+C IPAR(2) 0 Continuum production TWOG0984
+C 1 Resonance formation TWOG0985
+C TWOG0986
+C EXTERNALS: none TWOG0987
+C ========== TWOG0988
+C TWOG0989
+C AUTHOR/DATE: S.A. Sadovsky, 17 January 1995 TWOG0990
+C =========== TWOG0991
+C TWOG0992
+C***********************************************************************TWOG0993
+C IMPLICIT NONE TWOG0994
+ INTEGER ICNTER,NCNTRS,NEVTS,LOUT,DEBUG TWOG0995
+ PARAMETER (NCNTRS =7) TWOG0996
+ COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX, TWOG0997
+ + ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS, TWOG0998
+ + WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES, TWOG0999
+ + ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG TWOG1000
+ SAVE /TWGENC/ TWOG1001
+ REAL ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT TWOG1002
+ REAL WMIN,WMAX,OMMIN,OMRAT,EPS,TH1PMN,TH1PMX TWOG1003
+ REAL TH2PMN,TH2PMX,S,EB,WGHTSM,WGHTMX,WGHTGT,WTMXFD TWOG1004
+ REAL XMRES,XGRES TWOG1005
+ LOGICAL LRES,LWATE TWOG1006
+C TWOG1007
+ REAL ALPHA,PI,TWOPI,XME,XMU TWOG1008
+ PARAMETER (ALPHA=1./137.036,PI=3.14159265) TWOG1009
+ PARAMETER (TWOPI=2.*PI,XME=0.000511,XMU=0.105658) TWOG1010
+C TWOG1011
+ REAL XPAR(11),YPAR(6) TWOG1012
+ INTEGER LOUTX,IPAR(10),I TWOG1013
+C TWOG1014
+ EB = XPAR(1) TWOG1015
+ WMIN = XPAR(2) TWOG1016
+ WMAX = XPAR(3) TWOG1017
+ WGHTMX= XPAR(8) TWOG1022
+ WTMXFD= 0 TWOG1023
+ XMRES = XPAR(9) TWOG1024
+ XGRES = XPAR(10) TWOG1025
+ LOUT = LOUTX TWOG1026
+ LWATE = IPAR(1) .EQ. 1 TWOG1027
+ LRES = IPAR(2) .EQ. 1 TWOG1028
+ EPS = 1./64.*(XME*WMIN**2/EB**3)**2/100. TWOG1029
+C+ TWOG1030
+C EPS is a small constant. We throw 1/2 COS(th/2)/(eps+SIN(th/2)) TWOG1031
+C rather than 1/2 COT(th/2) to be able to start theta from zero. TWOG1032
+C- TWOG1033
+ S = (2.*EB)**2 TWOG1054
+ OMMIN = WMIN**2/(4.*EB) TWOG1055
+ OMRAT = EB/OMMIN TWOG1056
+C+ TWOG1057
+C Print out initial conditions. TWOG1058
+C- TWOG1059
+ WRITE(LOUT,1000) TWOG1060
+ 1000 FORMAT(1H ,10X,27('-'),' TWINIT ',27('-')) TWOG1061
+ WRITE(LOUT,1001) TWOG1062
+ 1001 FORMAT(1H ,10X,'I',60X,'I') TWOG1063
+C
+ WRITE(LOUT,1002)EB,YPAR(3),YPAR(4),YPAR(1),YPAR(2),WMIN,WMAX,
+ + YPAR(5),YPAR(6)
+C TWOG1065
+ 1002 FORMAT TWOG1066
+ + (1H ,10X,'I Beam energy ',E15.5, TWOG1067
+ + ' GeV',11X,'I',/, TWOG1068
+ + 1H ,10X,'I Mass of ion ',F12.5, TWOG1067
+ + ' GeV',11X,'I',/, TWOG1068
+ + 1H ,10X,'I Gamma of ion ',F12.5, TWOG1067
+ + ' ',11X,'I',/, TWOG1068
+ + 1H ,10X,'I Charge of ion ',F12.5, TWOG1067
+ + ' ',11X,'I',/, TWOG1068
+ + 1H ,10X,'I Atomic number of ion ',F12.5, TWOG1067
+ + ' ',11X,'I',/, TWOG1068
+ + 1H ,10X,'I Minimum gamma-gamma mass ',F12.5, TWOG1069
+ + ' GeV',11X,'I',/, TWOG1070
+ + 1H ,10X,'I Maximum gamma-gamma mass ',F12.5, TWOG1071
+ + ' GeV',11X,'I',/, TWOG1080
+ + 1H ,10X,'I Minimum rapidity ',F12.5, TWOG1069
+ + ' ',11X,'I',/, TWOG1070
+ + 1H ,10X,'I Maximum rapidity ',F12.5, TWOG1071
+ + ' ',11X,'I' ) TWOG1080
+ IF (LRES) THEN TWOG1081
+ WRITE(LOUT,1001) TWOG1082
+ WRITE(LOUT,1003)XMRES,XGRES TWOG1083
+ ENDIF TWOG1084
+ 1003 FORMAT TWOG1085
+ + (1H ,10X,'I Resonance mass ',F12.5, TWOG1086
+ + ' GeV',11X,'I',/, TWOG1087
+ + 1H ,10X,'I Resonance total width ',F12.5, TWOG1088
+ + ' GeV',11X,'I') TWOG1089
+C
+ IF (LWATE) THEN TWOG1090
+ WRITE(LOUT,1001) TWOG1091
+ WRITE(LOUT,1004) TWOG1092
+ ELSE TWOG1095
+ WRITE(LOUT,1001) TWOG1096
+ WRITE(LOUT,1005)WGHTMX TWOG1097
+ 1004 FORMAT(1H ,10X,'I Produce weighted events ',27X,'I') TWOG1094
+ 1005 FORMAT(1H ,10X,'I Maximum weight ',F12.5, TWOG1099
+ + ' ',11X,'I') TWOG1100
+ ENDIF TWOG1101
+C TWOG1102
+ WRITE(LOUT,1001) TWOG1103
+ WRITE(LOUT,1000) TWOG1104
+C TWOG1105
+ DO 10 I=1,NCNTRS TWOG1106
+ 10 ICNTER(I)=0 TWOG1107
+ WGHTGT = 0. TWOG1108
+ WGHTSM = 0. TWOG1109
+ RETURN TWOG1110
+ END TWOG1111
+C=========================================================================
+CDECK ID>, DLDMX.
+ FUNCTION DLDMX(Rm,Ny)
+ PARAMETER (isize=100000)
+C S.A.Sadovsky
+C Calculate the Luminosity integral over dY 5.02.1995
+C in the rapidity range ymin-ymax
+C using luminosity interpolation table
+C Ny/2 - number of subintervals for Simpson integration
+C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny,
+ + xmin,xmax,amn,amx,nnms,alfun(isize)
+ stpy=(ymx-ymn)/nny
+ stpm=(amx-amn)/nnms ! Interpolation steps
+C
+ im=(rm-amn)/stpm
+ q =(rm-im*stpm-amn)/stpm
+C
+ My = Ny/2
+ Hy =(ymax-ymin)/My ! Integration step over Y
+ ry = ymin
+ iy =(ry-ymn)/stpy+1.
+ p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 1
+ k = iy+im*(nny+1)
+ F1 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
+ * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
+C
+ SIM = 0.0
+ DO 10 i=1,My
+C
+ ry = ymin + Hy*i
+ iy =(ry-ymn)/stpy+1.
+ p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 3
+ k = iy+im*(nny+1)
+ F3 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
+ * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
+C
+ ry = ry - Hy*0.5
+ iy =(ry-ymn)/stpy+1.
+ p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 2
+ k = iy+im*(nny+1)
+ F2 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+
+ * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny)
+C
+ SIM = SIM + (F1 + 4.*F2 + F3)/6.
+ 10 F1 = F3
+C
+ DLDMX = 1000.*SIM*HY ! dL/dMx in GeV^-1
+ RETURN
+ END
+*=======================================================================
+CDECK ID>, LUMFUN.
+ SUBROUTINE LUMFUN(ifl,file,ymin,ymax,ny,xmn,xmx,nms)
+C-----------------------------------------------------------------------
+C preparation of luminosity function for generator
+C
+C ifl=-1 - calculate function + save in file 'file'
+C ifl= 0 - calculate function
+C ifl= 1 - read function from file
+C ymin,ymax,ny - rapidity range & number of points
+C xmn,xmx,nms - mass range & number of points
+C Attention: dimensions of alfun are alfun(ny+1,nms+1)
+C
+C Called by GGINIT
+C
+C Author: G.V.Khaoustov 12.12.1994
+C Modified: Yu.Kharlov 18-JUN-1997
+C-----------------------------------------------------------------------
+ COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max,
+ & gvconst(4,10)
+ REAL nc
+ COMMON/D2LParam/lz,la,Rion,Gamma
+ INTEGER lz,la
+ REAL Rion,Gamma
+ PARAMETER (isize=100000)
+ CHARACTER*(*) file
+ REAL work(1000)
+ COMMON/lumfunct/anorm,cexp,umn,umx,ymn,ymx,nny,
+ + wmn,wmx,amn,amx,nnms,alfun(isize)
+C
+ umn = ymin ! parameters for luminosity generation
+ umx = ymax
+ wmn = xmn
+ wmx = xmx
+C
+C Define R : radii of the ions (in fm) (if set non-zero, use it)
+ IF (Rion .LE. 0) THEN
+ rion = 1.2*float(la)**(1.0/3.0)
+ ENDIF
+*
+* Fing max energy of gg-system (with safity factor 2)
+*
+ egg_max = 2.*hbarc * gamma / rion
+*
+ IF(ifl.EQ.1) THEN ! read lum. function from file
+ OPEN(unit=10,err=1,file=file,status='old',
+ * form='formatted')
+ READ(10, *,err=3) anorm,cexp
+ READ(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
+ READ(10,100,err=3) (alfun(I),i=1,(nny+1)*(nnms+1))
+ CLOSE(unit=10)
+ PRINT *,'The Luminosity has been read from File '
+ RETURN
+ ENDIF
+c
+ IF (ny*nms .GT. isize) GOTO 2
+ PRINT *,'Start the Luminosity calculation'
+ nny = ny
+ ymn = ymin
+ ymx = ymax
+ nnms= nms
+ amn = xmn
+ amx = xmx
+c
+ stepy=(ymx-ymn)/nny
+ stepm=(amx-amn)/nnms
+*
+* Lum. function calculation
+*
+ DO IM = 0,nnms
+ am=amn+im*stepm
+* Find Y-limits for mass AM (Yu.Kharlov)
+ ym_max = alog (egg_max/am)
+ ym_max = max(ym_max,5.*stepy)
+ ygg_max = min (ymax, ym_max)
+ ygg_min = max (ymin,-ym_max)
+ DO IY = 1,nny+1
+ Y = ymn+(iy-1)*stepy
+ IF (Y.LT.ygg_max .AND. Y.GT.ygg_min) THEN
+ dl = D2LDMDY(1000.*am,Y) ! Differential Luminosity in MeV^-1
+ ELSE
+ dl = 0.0
+ END IF
+ alfun(iy+(nny+1)*im) = DL
+ ENDDO
+ ENDDO
+c
+ s=0.
+ DO i=1,nnms+1 ! find max values for exp. parameters calcul.
+ am=0.
+ DO j=1,nny+1
+ IF(am.LT.alfun(j+(nny+1)*(i-1))) am=alfun(j+(nny+1)*(i-1))
+ ENDDO
+ work(i)=am
+ s=s+am
+ ENDDO
+ s=(s-work(nnms+1))*stepm
+c
+ smax=0. ! calculation optimized parameters for exp. generator
+ iflag=-1
+ DO 101 i=1,nnms
+ DO 102 j=i+1,nnms+1
+ IF (work(j) .LE. 0.0) GOTO 102
+ c=alog(work(i)/work(j))/((i-j)*stepm)
+ am=work(i)/exp(c*(i-1)*stepm)
+ se=am/c*(exp(c*(amx-amn))-1.)
+ IF(s/se.lt.smax) GOTO 102
+ DO k=1,nnms+1
+ amp=am*exp(c*(k-1)*stepm)
+ IF((work(k)-amp)/amp.GT.0.00001) THEN ! Accuracy problem
+ GOTO 102
+ ENDIF
+ ENDDO
+ smax =s/se
+ anorm=am
+ cexp =c
+ iflag=1
+ 102 CONTINUE
+ 101 CONTINUE
+ IF(iflag.LT.0) THEN
+ PRINT *,'Fault to find best exponential approximation'
+ STOP
+ ENDIF
+c
+ IF(ifl.EQ.0) RETURN ! save lum. function
+ OPEN(unit=10,err=1,file=file,status='unknown',
+ * form='formatted')
+ WRITE(10, *,err=3) anorm,cexp
+ WRITE(10,99,err=3) ymn,ymx,nny,amn,amx,nnms
+ WRITE(10,100,err=3) (alfun(i),i=1,(nny+1)*(nnms+1))
+ CLOSE(unit=10)
+ PRINT *,'The Luminosity calculation is finished'
+ RETURN
+ 1 PRINT *,' file open error, file:',file
+ STOP
+ 2 PRINT *,' array "ALFUN" to small to contain luminosity function'
+ STOP
+ 3 PRINT *,' file read/writr error, file:',file
+ STOP
+ 100 FORMAT(5e13.5)
+ 99 FORMAT(2e13.5,i5,2e13.5,i5)
+ END
+C======================================================================
+CDECK ID>, D2LDMDY.
+ REAL FUNCTION D2LDMDY(M,Y)
+C-----------------------------------------------------------------------
+C Calculate the Luminosity-function d^2L/dMdY
+C where M is the invariant mass and Y the rapidity
+C of the photon-photon-system
+C
+C Calculations are based on the formulae of
+C N. Baron, G. Baur, Phys. Rev. C
+C C. Cahn, J. D. Jackson, Nucl. Phys.
+C
+C Kai Hencken, 29. 9. 1994
+C
+C M : invariant mass (MeV)
+C Y : rapidity
+C
+C D2LDMDY : diff. Lum. (MeV**-1)
+C
+C Called by LUMFUN
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ COMMON/D2LParam/lz,la,Rion,Gamma
+ INTEGER lz,la
+ REAL Rion,Gamma
+ REAL M,Y
+ REAL RM,W1,W2
+ COMMON/D2LIParam/RM,W1,W2
+
+ REAL NClass,DeltaL
+ EXTERNAL NClass,DeltaL
+
+C
+C physical constants hbar*c and 1/alpha
+C
+ REAL hbarc,falphai
+ PARAMETER (hbarc=197.327053,falphai=137.0359895)
+
+ REAL LClass,DL
+C
+C calculate photon frequencies from M,Y
+C
+ W1 = M/2.0 * EXP ( Y)
+ W2 = M/2.0 * EXP (-Y)
+C
+C calculate Rho in MeV^-1 instead of fm
+C
+ RM = Rion / hbarc
+C
+C calculate the "classical" value first.
+C
+
+ LClass = NClass (W1,RM,Gamma) * NClass (W2,RM,Gamma)
+
+C
+C substract from this the Delta-value
+C
+ DL = DeltaL ()
+
+ D2LDMDY = 2.0 / M * float(LZ)**4 / falphai**2 * (LClass - DL)
+
+ RETURN
+ END
+*=======================================================================
+CDECK ID>, NCLASS.
+ REAL FUNCTION NClass (W, R, Gamma)
+C
+C the classical photon-number without impact parameter cutoff:
+C
+ IMPLICIT NONE
+ REAL W,R, Gamma
+
+ REAL Pi
+ PARAMETER (Pi = 3.141592)
+
+ REAL BesselK0,BesselK1
+ EXTERNAL BesselK0,BesselK1
+
+ REAL Xi
+
+ Xi = W*R/Gamma
+
+ NClass = 2.0/PI * (Xi*BesselK0 (Xi)*BesselK1 (Xi) -
+ & Xi**2/2.0 *(BesselK1(Xi)**2 - BesselK0(Xi)**2))
+
+ RETURN
+ END
+*=======================================================================
+CDECK ID>, DELTAL.
+ REAL FUNCTION DeltaL ()
+C
+C the difference to the classical value
+C Called by d2LdMdY
+C
+ IMPLICIT NONE
+ COMMON/D2LParam/lz,la,Rion,Gamma
+ INTEGER lz,la
+ REAL Rion,Gamma
+ REAL RM,W1,W2
+ COMMON/D2LIParam/RM,W1,W2
+
+ DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
+ COMMON /GAUSSD/XGAUSS,WGAUSS
+
+ INTEGER N,I
+ REAL b1
+ REAL*8 t,tmin,tmax
+ REAL*8 Sum,Int,Int2
+
+ REAL Itgb1
+ EXTERNAL Itgb1
+
+ REAL Pi
+ REAL*8 Accur
+ PARAMETER (Pi = 3.141592,Accur = 1D -2)
+
+C
+C integrate first over b1
+C
+
+C
+C Loop incrementing the boundary
+C
+ tmin = 0.00
+ tmax = 0.25
+ Sum = 0.00
+
+ 100 CONTINUE
+C
+C Loop for the Gauss integration
+C
+ Int=0.0
+ DO N=1,6
+ Int2 = Int
+ Int=0.0
+ DO I=2**N-1,2**(N+1)-2
+ t = (tmax-tmin)/2.0*XGAUSS(I) + (tmax+tmin)/2.0
+ b1 = RM * DEXP (t)
+ Int = Int + WGAUSS(I) * Itgb1(b1) * b1**2
+ ENDDO
+ Int = (tmax-tmin)/2.0*Int
+ IF (Int .NE. 0.) THEN
+ IF (DABS ((Int2-Int)/Int) .LT. Accur) GOTO 200
+ END IF
+ ENDDO
+ IF (int .LT. 1.D-20) THEN
+ int = 0.0D0
+ ELSE
+ WRITE(*,*) ' (b1) GAUSS MAY BE INACCURATE'
+ END IF
+
+ 200 CONTINUE
+ Sum = Sum + Int
+ IF (sum .NE. 0.0) THEN
+ IF (ABS (Int2/Sum) .LT. Accur) GOTO 300
+ ELSE
+ IF (b1 .GT. 1.E+06) GOTO 300
+ END IF
+ tmin = tmax
+ tmax = tmax + 0.5
+ GOTO 100
+
+ 300 CONTINUE
+
+ DeltaL = 4.0*Pi * Sum
+ RETURN
+ END
+
+CDECK ID>, ITGB1.
+ REAL FUNCTION Itgb1(b)
+C
+C integration then over b2
+C Called by DeltaL
+C
+C Corrected by S.A.Sadovsky
+C 25.10.1994
+ IMPLICIT NONE
+c +seq,d2lpar. Not used Yu.Kh. 18-JUN-1997
+ REAL b
+ REAL*8 b1,bmin,bmax
+ REAL RM,W1,W2
+ COMMON/D2LIParam/RM,W1,W2
+ DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
+ COMMON /GAUSSD/XGAUSS,WGAUSS
+ INTEGER N,I
+ REAL*8 b2,Int,Int2
+ REAL*8 Itgb2
+ EXTERNAL Itgb2
+ REAL*8 Accur
+ PARAMETER (Accur = 1.D-2)
+
+ b1 = b
+ bmin = b1 - 2.0*RM
+ IF (RM .GT. bmin) THEN
+ bmin = RM
+ ENDIF
+ bmax = b1 + 2.0 * RM
+ Int = 0.0
+ DO N=1,6
+ Int2 = Int
+ Int = 0.0
+ DO I=2**N-1,2**(N+1)-2
+ b2 = (bmax-bmin)/2.0*XGAUSS(I) + (bmax+bmin)/2.0
+ Int = Int + WGAUSS(I) * Itgb2(b1,b2) * b2
+ END DO
+ Int = (bmax-bmin)/2.0*Int
+ IF (Int.NE.0.) THEN
+ IF (DABS((Int2 - Int)/Int) .LT. Accur) GOTO 100
+ ENDIF
+ ENDDO
+ IF (int .LT. 1.D-20) THEN
+ int = 0.0D0
+ ELSE
+ WRITE(*,*) ' (b2) GAUSS MAY BE INACCURATE, Itgb1=',Int
+ END IF
+ 100 CONTINUE
+ Itgb1 = Int
+ RETURN
+ END
+
+CDECK ID>, ITGB2.
+ REAL*8 FUNCTION Itgb2 (b1,b2)
+C
+C The function to be integrated over
+C Called by Itgb1
+C
+C Corrected by S.A.Sadovsky
+C 25.10.1994
+ IMPLICIT NONE
+ COMMON/D2LParam/lz,la,Rion,Gamma
+ INTEGER lz,la
+ REAL Rion,Gamma
+ REAL*8 b1,b2,arg
+ REAL RM,W1,W2
+ COMMON/D2LIParam/RM,W1,W2
+
+ REAL Nphoton
+ EXTERNAL Nphoton
+
+ arg = (b1**2 + b2**2 - 4.0 * RM**2) / (2.0*b1*b2)
+ IF (arg .GT. 1.D0) arg = 1.D0
+ Itgb2 = Nphoton(W1,sngl(b1),Gamma) *
+ & Nphoton(W2,sngl(b2),Gamma) *
+ & DACOS (arg)
+
+ RETURN
+ END
+
+CDECK ID>, NPHOTON.
+ REAL FUNCTION Nphoton (W,Rho,Gamma)
+C
+C The differential photonnumber for a nucleus
+C (without form factor)
+C
+ IMPLICIT NONE
+ REAL W,Rho,Gamma
+
+ REAL BESSELK1
+ EXTERNAL BESSELK1
+
+ REAL Wphib,WGamma
+
+ REAL PI
+ PARAMETER (PI=3.141592)
+
+ WGamma = W/Gamma
+C
+C without form factor
+C
+ Wphib = WGamma*BESSELK1 (WGamma*Rho)
+
+ Nphoton = 1.0/PI**2 * Wphib**2
+
+ RETURN
+ END
+
+CDECK ID>, BESSELK0.
+ REAL FUNCTION BESSELK0(X)
+C
+C Special functions I0,K0,I1,K1
+C see Abramowitz&Stegun
+C
+ IMPLICIT NONE
+ REAL X,Y
+
+ REAL BESSELI0
+ EXTERNAL BESSELI0
+
+ IF (X .LT. 2.0) THEN
+
+ Y = X**2/4.0
+
+ BESSELK0 =
+ & (-LOG(X/2.0)*BESSELI0(X))+(-.57721566+Y*(0.42278420
+ & +Y*(0.23069756+Y*(0.3488590E-1+Y*(0.262698E-2
+ & +Y*(0.10750E-3+Y*0.740E-5))))))
+
+ ELSE
+
+ Y = 2.0/X
+
+ BESSELK0 =
+ & (EXP(-X)/SQRT(X))*(1.25331414+Y*(-0.7832358E-1
+ & +Y*(0.2189568E-1+Y*(-0.1062446E-1+Y*(0.587872E-2
+ & +Y*(-0.251540E-2+Y*0.53208E-3))))))
+
+ ENDIF
+
+ RETURN
+ END
+
+
+CDECK ID>, BESSELI0.
+ REAL FUNCTION BESSELI0(X)
+ IMPLICIT NONE
+
+ REAL X,Y,AX
+
+ AX = ABS(X)
+
+ IF (AX .LT. 3.75) THEN
+
+ Y = (X/3.75)**2
+
+ BESSELI0 =
+ & 1.0+Y*(3.5156229+Y*(3.0899424+Y*(1.2067492
+ & +Y*(0.2659732+Y*(0.360768E-1+Y*0.45813E-2)))))
+
+ ELSE
+
+ Y = 3.75/AX
+
+ BESSELI0 =
+ & (EXP(AX)/SQRT(AX))*(0.39894228+Y*(0.1328592E-1
+ & +Y*(0.225319E-2+Y*(-0.157565E-2+Y*(0.916281E-2
+ & +Y*(-0.2057706E-1+Y*(0.2635537E-1+Y*(-0.1647633E-1
+ & +Y*0.392377E-2))))))))
+
+ ENDIF
+
+ RETURN
+ END
+
+
+CDECK ID>, BESSELK1.
+ REAL FUNCTION BESSELK1(X)
+ IMPLICIT NONE
+
+ REAL X,Y
+
+ REAL BESSELI1
+ EXTERNAL BESSELI1
+
+ IF (X .LT. 2.0) THEN
+
+ Y = X**2/4.0
+ BESSELK1 =
+ & (LOG(X/2.0)*BESSELI1(X))+(1.0/X)*(1.0+Y*(0.15443144
+ & +Y*(-0.67278579+Y*(-0.18156897+Y*(-0.1919402E-1
+ & +Y*(-0.110404E-2+Y*(-0.4686E-4)))))))
+
+ ELSE
+
+ Y=2.0/X
+ BESSELK1 =
+ & (EXP(-X)/SQRT(X))*(1.25331414+Y*(0.23498619
+ & +Y*(-0.3655620E-1+Y*(0.1504268E-1+Y*(-0.780353E-2
+ & +Y*(0.325614E-2+Y*(-0.68245E-3)))))))
+
+ ENDIF
+
+ RETURN
+ END
+
+CDECK ID>, BESSELI1.
+ REAL FUNCTION BESSELI1(X)
+ IMPLICIT NONE
+
+ REAL X,Y,AX
+
+ AX = ABS(X)
+
+ IF (AX .LT. 3.75) THEN
+
+ Y = (X/3.75)**2
+
+ BESSELI1 =
+ & AX*(0.5+Y*(0.87890594+Y*(0.51498869+Y*(0.15084934
+ & +Y*(0.2658733E-1+Y*(0.301532E-2+Y*0.32411E-3))))))
+
+
+ ELSE
+
+ Y = 3.75/AX
+
+ BESSELI1 =
+ & 0.2282967E-1+Y*(-0.2895312E-1+Y*(0.1787654E-1
+ & -Y*0.420059E-2))
+
+ BESSELI1 =
+ & 0.39894228+Y*(-0.3988024E-1+Y*(-0.362018E-2
+ & +Y*(0.163801E-2+Y*(-0.1031555E-1+Y*BESSELI1))))
+
+ BESSELI1 = BESSELI1 *
+ & EXP(AX)/SQRT(AX)
+
+ ENDIF
+
+ IF (X .LT. 0) THEN
+ BESSELI1 = -BESSELI1
+ ENDIF
+
+ RETURN
+ END
+
+CDECK ID>, GAUSSDAT.
+C DATA for the Gauss integration, x-values and weight for a
+C N=2,4,8,16,32,64 point Gauss integration.
+C
+C based on program GAUSS64 of D. Trautmann
+C data from Abramowitz & Stegun
+C
+C KAI HENCKEN, FEBRUAR 1993
+C
+ BLOCKDATA GAUSSDATA
+ IMPLICIT NONE
+ DOUBLE PRECISION XGAUSS(126),WGAUSS(126)
+ COMMON /GAUSSD/XGAUSS,WGAUSS
+
+ DATA XGAUSS(1)/ .57735026918962576D0/
+ DATA XGAUSS(2)/-.57735026918962576D0/
+ DATA WGAUSS(1)/ 1.00000000000000000D0/
+ DATA WGAUSS(2)/ 1.00000000000000000D0/
+
+ DATA XGAUSS(3)/ .33998104358485627D0/
+ DATA XGAUSS(4)/ .86113631159405258D0/
+ DATA XGAUSS(5)/-.33998104358485627D0/
+ DATA XGAUSS(6)/-.86113631159405258D0/
+ DATA WGAUSS(3)/ .65214515486254613D0/
+ DATA WGAUSS(4)/ .34785484513745385D0/
+ DATA WGAUSS(5)/ .65214515486254613D0/
+ DATA WGAUSS(6)/ .34785484513745385D0/
+
+ DATA XGAUSS(7)/ .18343464249564981D0/
+ DATA XGAUSS(8)/ .52553240991632899D0/
+ DATA XGAUSS(9)/ .79666647741362674D0/
+ DATA XGAUSS(10)/ .96028985649753623D0/
+ DATA XGAUSS(11)/-.18343464249564981D0/
+ DATA XGAUSS(12)/-.52553240991632899D0/
+ DATA XGAUSS(13)/-.79666647741362674D0/
+ DATA XGAUSS(14)/-.96028985649753623D0/
+ DATA WGAUSS(7)/ .36268378337836198D0/
+ DATA WGAUSS(8)/ .31370664587788727D0/
+ DATA WGAUSS(9)/ .22238103445337448D0/
+ DATA WGAUSS(10)/ .10122853629037627D0/
+ DATA WGAUSS(11)/ .36268378337836198D0/
+ DATA WGAUSS(12)/ .31370664587788727D0/
+ DATA WGAUSS(13)/ .22238103445337448D0/
+ DATA WGAUSS(14)/ .10122853629037627D0/
+
+ DATA XGAUSS(15)/ .0950125098376374402D0/
+ DATA XGAUSS(16)/ .281603550779258913D0/
+ DATA XGAUSS(17)/ .458016777657227386D0/
+ DATA XGAUSS(18)/ .617876244402643748D0/
+ DATA XGAUSS(19)/ .755404408355003034D0/
+ DATA XGAUSS(20)/ .865631202387831744D0/
+ DATA XGAUSS(21)/ .944575023073232576D0/
+ DATA XGAUSS(22)/ .989400934991649933D0/
+ DATA XGAUSS(23)/-.0950125098376374402D0/
+ DATA XGAUSS(24)/-.281603550779258913D0/
+ DATA XGAUSS(25)/-.458016777657227386D0/
+ DATA XGAUSS(26)/-.617876244402643748D0/
+ DATA XGAUSS(27)/-.755404408355003034D0/
+ DATA XGAUSS(28)/-.865631202387831744D0/
+ DATA XGAUSS(29)/-.944575023073232576D0/
+ DATA XGAUSS(30)/-.989400934991649933D0/
+ DATA WGAUSS(15)/ .189450610455068496D0/
+ DATA WGAUSS(16)/ .182603415044923589D0/
+ DATA WGAUSS(17)/ .169156519395002538D0/
+ DATA WGAUSS(18)/ .149595988816576732D0/
+ DATA WGAUSS(19)/ .124628971255533872D0/
+ DATA WGAUSS(20)/ .0951585116824927848D0/
+ DATA WGAUSS(21)/ .0622535239386478929D0/
+ DATA WGAUSS(22)/ .0271524594117540949D0/
+ DATA WGAUSS(23)/ .189450610455068496D0/
+ DATA WGAUSS(24)/ .182603415044923589D0/
+ DATA WGAUSS(25)/ .169156519395002538D0/
+ DATA WGAUSS(26)/ .149595988816576732D0/
+ DATA WGAUSS(27)/ .124628971255533872D0/
+ DATA WGAUSS(28)/ .0951585116824927848D0/
+ DATA WGAUSS(29)/ .0622535239386478929D0/
+ DATA WGAUSS(30)/ .0271524594117540949D0/
+
+ DATA XGAUSS(31)/ .0483076656877383162D0/
+ DATA XGAUSS(32)/ .144471961582796493D0/
+ DATA XGAUSS(33)/ .239287362252137075D0/
+ DATA XGAUSS(34)/ .331868602282127650D0/
+ DATA XGAUSS(35)/ .421351276130635345D0/
+ DATA XGAUSS(36)/ .506899908932229390D0/
+ DATA XGAUSS(37)/ .587715757240762329D0/
+ DATA XGAUSS(38)/ .663044266930215201D0/
+ DATA XGAUSS(39)/ .732182118740289680D0/
+ DATA XGAUSS(40)/ .794483795967942407D0/
+ DATA XGAUSS(41)/ .849367613732569970D0/
+ DATA XGAUSS(42)/ .896321155766052124D0/
+ DATA XGAUSS(43)/ .934906075937739689D0/
+ DATA XGAUSS(44)/ .964762255587506430D0/
+ DATA XGAUSS(45)/ .985611511545268335D0/
+ DATA XGAUSS(46)/ .997263861849481564D0/
+ DATA XGAUSS(47)/-.0483076656877383162D0/
+ DATA XGAUSS(48)/-.144471961582796493D0/
+ DATA XGAUSS(49)/-.239287362252137075D0/
+ DATA XGAUSS(50)/-.331868602282127650D0/
+ DATA XGAUSS(51)/-.421351276130635345D0/
+ DATA XGAUSS(52)/-.506899908932229390D0/
+ DATA XGAUSS(53)/-.587715757240762329D0/
+ DATA XGAUSS(54)/-.663044266930215201D0/
+ DATA XGAUSS(55)/-.732182118740289680D0/
+ DATA XGAUSS(56)/-.794483795967942407D0/
+ DATA XGAUSS(57)/-.849367613732569970D0/
+ DATA XGAUSS(58)/-.896321155766052124D0/
+ DATA XGAUSS(59)/-.934906075937739689D0/
+ DATA XGAUSS(60)/-.964762255587506430D0/
+ DATA XGAUSS(61)/-.985611511545268335D0/
+ DATA XGAUSS(62)/-.997263861849481564D0/
+ DATA WGAUSS(31)/ .0965400885147278006D0/
+ DATA WGAUSS(32)/ .0956387200792748594D0/
+ DATA WGAUSS(33)/ .0938443990808045654D0/
+ DATA WGAUSS(34)/ .0911738786957638847D0/
+ DATA WGAUSS(35)/ .0876520930044038111D0/
+ DATA WGAUSS(36)/ .0833119242269467552D0/
+ DATA WGAUSS(37)/ .0781938957870703065D0/
+ DATA WGAUSS(38)/ .0723457941088485062D0/
+ DATA WGAUSS(39)/ .0658222227763618468D0/
+ DATA WGAUSS(40)/ .0586840934785355471D0/
+ DATA WGAUSS(41)/ .0509980592623761762D0/
+ DATA WGAUSS(42)/ .0428358980222266807D0/
+ DATA WGAUSS(43)/ .0342738629130214331D0/
+ DATA WGAUSS(44)/ .0253920653092620595D0/
+ DATA WGAUSS(45)/ .0162743947309056706D0/
+ DATA WGAUSS(46)/ .00701861000947009660D0/
+ DATA WGAUSS(47)/ .0965400885147278006D0/
+ DATA WGAUSS(48)/ .0956387200792748594D0/
+ DATA WGAUSS(49)/ .0938443990808045654D0/
+ DATA WGAUSS(50)/ .0911738786957638847D0/
+ DATA WGAUSS(51)/ .0876520930044038111D0/
+ DATA WGAUSS(52)/ .0833119242269467552D0/
+ DATA WGAUSS(53)/ .0781938957870703065D0/
+ DATA WGAUSS(54)/ .0723457941088485062D0/
+ DATA WGAUSS(55)/ .0658222227763618468D0/
+ DATA WGAUSS(56)/ .0586840934785355471D0/
+ DATA WGAUSS(57)/ .0509980592623761762D0/
+ DATA WGAUSS(58)/ .0428358980222266807D0/
+ DATA WGAUSS(59)/ .0342738629130214331D0/
+ DATA WGAUSS(60)/ .0253920653092620595D0/
+ DATA WGAUSS(61)/ .0162743947309056706D0/
+ DATA WGAUSS(62)/ .00701861000947009660D0/
+
+ DATA XGAUSS(63)/ .02435029266342443250D0/
+ DATA XGAUSS(64)/ .0729931217877990394D0/
+ DATA XGAUSS(65)/ .121462819296120554D0/
+ DATA XGAUSS(66)/ .169644420423992818D0/
+ DATA XGAUSS(67)/ .217423643740007084D0/
+ DATA XGAUSS(68)/ .264687162208767416D0/
+ DATA XGAUSS(69)/ .311322871990210956D0/
+ DATA XGAUSS(70)/ .357220158337668116D0/
+ DATA XGAUSS(71)/ .402270157963991604D0/
+ DATA XGAUSS(72)/ .446366017253464088D0/
+ DATA XGAUSS(73)/ .489403145707052957D0/
+ DATA XGAUSS(74)/ .531279464019894546D0/
+ DATA XGAUSS(75)/ .571895646202634034D0/
+ DATA XGAUSS(76)/ .611155355172393250D0/
+ DATA XGAUSS(77)/ .648965471254657340D0/
+ DATA XGAUSS(78)/ .685236313054233243D0/
+ DATA XGAUSS(79)/ .719881850171610827D0/
+ DATA XGAUSS(80)/ .752819907260531897D0/
+ DATA XGAUSS(81)/ .783972358943341408D0/
+ DATA XGAUSS(82)/ .813265315122797560D0/
+ DATA XGAUSS(83)/ .840629296252580363D0/
+ DATA XGAUSS(84)/ .865999398154092820D0/
+ DATA XGAUSS(85)/ .889315445995114106D0/
+ DATA XGAUSS(86)/ .910522137078502806D0/
+ DATA XGAUSS(87)/ .929569172131939576D0/
+ DATA XGAUSS(88)/ .946411374858402816D0/
+ DATA XGAUSS(89)/ .961008799652053719D0/
+ DATA XGAUSS(90)/ .973326827789910964D0/
+ DATA XGAUSS(91)/ .983336253884625957D0/
+ DATA XGAUSS(92)/ .991013371476744321D0/
+ DATA XGAUSS(93)/ .996340116771955279D0/
+ DATA XGAUSS(94)/ .999305041735772139D0/
+ DATA XGAUSS(95)/-.02435029266342443250D0/
+ DATA XGAUSS(96)/-.0729931217877990394D0/
+ DATA XGAUSS(97)/-.121462819296120554D0/
+ DATA XGAUSS(98)/-.169644420423992818D0/
+ DATA XGAUSS(99)/-.217423643740007084D0/
+ DATA XGAUSS(100)/-.264687162208767416D0/
+ DATA XGAUSS(101)/-.311322871990210956D0/
+ DATA XGAUSS(102)/-.357220158337668116D0/
+ DATA XGAUSS(103)/-.402270157963991604D0/
+ DATA XGAUSS(104)/-.446366017253464088D0/
+ DATA XGAUSS(105)/-.489403145707052957D0/
+ DATA XGAUSS(106)/-.531279464019894546D0/
+ DATA XGAUSS(107)/-.571895646202634034D0/
+ DATA XGAUSS(108)/-.611155355172393250D0/
+ DATA XGAUSS(109)/-.648965471254657340D0/
+ DATA XGAUSS(110)/-.685236313054233243D0/
+ DATA XGAUSS(111)/-.719881850171610827D0/
+ DATA XGAUSS(112)/-.752819907260531897D0/
+ DATA XGAUSS(113)/-.783972358943341408D0/
+ DATA XGAUSS(114)/-.813265315122797560D0/
+ DATA XGAUSS(115)/-.840629296252580363D0/
+ DATA XGAUSS(116)/-.865999398154092820D0/
+ DATA XGAUSS(117)/-.889315445995114106D0/
+ DATA XGAUSS(118)/-.910522137078502806D0/
+ DATA XGAUSS(119)/-.929569172131939576D0/
+ DATA XGAUSS(120)/-.946411374858402816D0/
+ DATA XGAUSS(121)/-.961008799652053719D0/
+ DATA XGAUSS(122)/-.973326827789910964D0/
+ DATA XGAUSS(123)/-.983336253884625957D0/
+ DATA XGAUSS(124)/-.991013371476744321D0/
+ DATA XGAUSS(125)/-.996340116771955279D0/
+ DATA XGAUSS(126)/-.999305041735772139D0/
+ DATA WGAUSS(63)/ .0486909570091397204D0/
+ DATA WGAUSS(64)/ .0485754674415034269D0/
+ DATA WGAUSS(65)/ .0483447622348029572D0/
+ DATA WGAUSS(66)/ .0479993885964583077D0/
+ DATA WGAUSS(67)/ .0475401657148303087D0/
+ DATA WGAUSS(68)/ .0469681828162100173D0/
+ DATA WGAUSS(69)/ .0462847965813144172D0/
+ DATA WGAUSS(70)/ .0454916279274181445D0/
+ DATA WGAUSS(71)/ .0445905581637565631D0/
+ DATA WGAUSS(72)/ .0435837245293234534D0/
+ DATA WGAUSS(73)/ .0424735151236535890D0/
+ DATA WGAUSS(74)/ .0412625632426235286D0/
+ DATA WGAUSS(75)/ .0399537411327203414D0/
+ DATA WGAUSS(76)/ .0385501531786156291D0/
+ DATA WGAUSS(77)/ .0370551285402400460D0/
+ DATA WGAUSS(78)/ .0354722132568823838D0/
+ DATA WGAUSS(79)/ .0338051618371416094D0/
+ DATA WGAUSS(80)/ .0320579283548515535D0/
+ DATA WGAUSS(81)/ .0302346570724024789D0/
+ DATA WGAUSS(82)/ .0283396726142594832D0/
+ DATA WGAUSS(83)/ .0263774697150546587D0/
+ DATA WGAUSS(84)/ .0243527025687108733D0/
+ DATA WGAUSS(85)/ .0222701738083832542D0/
+ DATA WGAUSS(86)/ .0201348231535302094D0/
+ DATA WGAUSS(87)/ .0179517157756973431D0/
+ DATA WGAUSS(88)/ .0157260304760247193D0/
+ DATA WGAUSS(89)/ .0134630478967186426D0/
+ DATA WGAUSS(90)/ .0111681394601311288D0/
+ DATA WGAUSS(91)/ .00884675982636394772D0/
+ DATA WGAUSS(92)/ .00650445796897836286D0/
+ DATA WGAUSS(93)/ .00414703326056246764D0/
+ DATA WGAUSS(94)/ .00178328072169643295D0/
+ DATA WGAUSS(95)/ .0486909570091397204D0/
+ DATA WGAUSS(96)/ .0485754674415034269D0/
+ DATA WGAUSS(97)/ .0483447622348029572D0/
+ DATA WGAUSS(98)/ .0479993885964583077D0/
+ DATA WGAUSS(99)/ .0475401657148303087D0/
+ DATA WGAUSS(100)/ .0469681828162100173D0/
+ DATA WGAUSS(101)/ .0462847965813144172D0/
+ DATA WGAUSS(102)/ .0454916279274181445D0/
+ DATA WGAUSS(103)/ .0445905581637565631D0/
+ DATA WGAUSS(104)/ .0435837245293234534D0/
+ DATA WGAUSS(105)/ .0424735151236535890D0/
+ DATA WGAUSS(106)/ .0412625632426235286D0/
+ DATA WGAUSS(107)/ .0399537411327203414D0/
+ DATA WGAUSS(108)/ .0385501531786156291D0/
+ DATA WGAUSS(109)/ .0370551285402400460D0/
+ DATA WGAUSS(110)/ .0354722132568823838D0/
+ DATA WGAUSS(111)/ .0338051618371416094D0/
+ DATA WGAUSS(112)/ .0320579283548515535D0/
+ DATA WGAUSS(113)/ .0302346570724024789D0/
+ DATA WGAUSS(114)/ .0283396726142594832D0/
+ DATA WGAUSS(115)/ .0263774697150546587D0/
+ DATA WGAUSS(116)/ .0243527025687108733D0/
+ DATA WGAUSS(117)/ .0222701738083832542D0/
+ DATA WGAUSS(118)/ .0201348231535302094D0/
+ DATA WGAUSS(119)/ .0179517157756973431D0/
+ DATA WGAUSS(120)/ .0157260304760247193D0/
+ DATA WGAUSS(121)/ .0134630478967186426D0/
+ DATA WGAUSS(122)/ .0111681394601311288D0/
+ DATA WGAUSS(123)/ .00884675982636394772D0/
+ DATA WGAUSS(124)/ .00650445796897836286D0/
+ DATA WGAUSS(125)/ .00414703326056246764D0/
+ DATA WGAUSS(126)/ .00178328072169643295D0/
+
+ END
+*=======================================================================
+ SUBROUTINE LORENB (U,PS,PI,PF)
+C
+C CERN PROGLIB# U102 LORENB .VERSION KERNFOR 4.04 821124
+C ORIG. 20/08/75 L.PAPE
+C
+ DOUBLE PRECISION PF4, FN
+ DIMENSION PS(4),PI(4),PF(4)
+
+ IF (PS(4).EQ.U) GO TO 17
+ PF4 = (PI(4)*PS(4)+PI(3)*PS(3)+PI(2)*PS(2)+PI(1)*PS(1)) / U
+ FN = (PF4+PI(4)) / (PS(4)+U)
+ PF(1)= PI(1) + FN*PS(1)
+ PF(2)= PI(2) + FN*PS(2)
+ PF(3)= PI(3) + FN*PS(3)
+ PF(4)= PF4
+ GO TO 18
+C
+ 17 PF(1)= PI(1)
+ PF(2)= PI(2)
+ PF(3)= PI(3)
+ PF(4)= PI(4)
+C
+ 18 CONTINUE
+C
+ RETURN
+C
+ END