ISAJET/code/qcdz.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE QCDZ(J)
   3 C
   4 C          Auxiliary routine for QCDJET.  Generate Z for parton J and
   5 C          store in ZZC(J). Add possible new partons to /JETSET/.
   6 C
   7 C          Include GM, W+, W-, and Z0 radiation.
   8 C
   9 C          Ver 7.20: Anomalous dimensions were coded incorrectly!
  10 C
  11 #if defined(CERNLIB_IMPNONE)
  12       IMPLICIT NONE
  13 #endif
  14 #include "isajet/itapes.inc"
  15 #include "isajet/jetset.inc"
  16 #include "isajet/jwork.inc"
  17 #include "isajet/qcdpar.inc"
  18 #include "isajet/wcon.inc"
  19 #include "isajet/const.inc"
  20 #include "isajet/q1q2.inc"
  21 C
  22       REAL PQQ,PGQ,PQG,PGG,Z,PGSGS,PQSQS,ALFAS,QQ,AM0,ZC,AMASS
  23       REAL GAMQQ,GAMGG,PROBG,PROBQ,RND,RANF,ZGEN,GZ
  24       REAL GAMZC,GAMSUM,AM1W,AM2W,T1W,T2W,ZCW,T0,GAMZCW,TERM,SUM
  25       REAL SUMBR,BRMODE,TRY,HELPL,HELMN,HEL,PZ
  26       INTEGER NF,J,JTABS,IQ,IFL,IW,JT0,JT1,IFL1,IFL2
  27       INTEGER IWTYPE,JET,JW,IQ1,IQ2,JPAR,IFLPAR,NJ1,NJ2,IDABS1,IDABS2
  28       REAL GAMSAV(5),ZCSAV(5),BRANCH(25)
  29       INTEGER JSAV(5),LISTW(5),LISTJ(25)
  30       DATA LISTW/9,10,80,-80,90/
  31       DATA LISTJ/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
  32      $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16/
  33 C
  34 C          Altarelli-Parisi functions.
  35       PQQ(Z)=4*(1+Z**2)/(3*(1-Z))
  36       PGQ(Z)=.5*(Z**2+(1-Z)**2)
  37       PGG(Z)=6*(1-Z*(1-Z))**2/(Z*(1-Z))
  38       PGSGS(Z)=3.*(1.+Z**2)/(1.-Z)
  39       PQSQS(Z)=8./3.*Z/(1.-Z)
  40       ALFAS(QQ)=12.*PI/((33.-2.*NF)*ALOG(QQ/ALAM2))
  41 C
  42 C          Initialize.
  43 C
  44       AM0=PJSET(5,J)
  45       ZC=ZZC(J)
  46       JTABS=IABS(JTYPE(J))
  47       NF=3
  48       DO 110 IQ=4,6
  49         IF(AM0.LT.2*AMASS(IQ)) GO TO 120
  50         NF=NF+1
  51 110   CONTINUE
  52 120   CONTINUE
  53       NJ1=NJSET+1
  54       NJ2=NJSET+2
  55 C
  56 C          Initial gluon
  57 C
  58       IF (JTABS.EQ.9) THEN
  59         GAMQQ=(1-2*ZC)*(1-ZC*(1-ZC))/3.
  60         GAMGG=12*ALOG((1-ZC)/ZC)-9*(1-2*ZC)-6*GAMQQ
  61         PROBG=GAMGG/(GAMGG+2*NF*GAMQQ)
  62         PROBQ=GAMQQ/(GAMGG+2*NF*GAMQQ)
  63         RND=RANF()
  64 C          GL--->GL+GL
  65         IF(PROBG.GT.RND) THEN
  66 130       ZGEN=(ZC/(1-ZC))**(1-2*RANF())
  67           Z=ZGEN/(1.+ZGEN)
  68           GZ=6./(Z*(1.-Z))
  69           IF(PGG(Z).LT.GZ*RANF()) GO TO 130
  70           JTYPE(NJ1)=9
  71           JTYPE(NJ2)=9
  72           ZZC(J)=Z
  73 C          GL--->QK+QB
  74         ELSE
  75 140       Z=RANF()
  76           IF(PGQ(Z).LT.0.5*RANF()) GO TO 140
  77           IFL=(RND-PROBG)/PROBQ+1.
  78           IF(IFL.GT.NF) IFL=NF-IFL
  79           JTYPE(NJ1)=IFL
  80           JTYPE(NJ2)=-IFL
  81           ZZC(J)=Z
  82         ENDIF
  83 C
  84 C          Initial quark - may radiate GL, GM, W+-, Z0
  85 C
  86       ELSEIF(JTABS.LT.9) THEN
  87 C          Gluon
  88         GAMZC=2.*ALOG((1-ZC)/ZC)-1.5*(1.-2.*ZC)
  89         GAMSAV(1)=4./3.*ALFAS(AM0**2)*GAMZC
  90         ZCSAV(1)=ZC
  91         JSAV(1)=JTYPE(J)
  92 C          Photon
  93         GAMSAV(2)=ALFA*AQ(JTABS,1)**2*GAMZC
  94         ZCSAV(2)=ZC
  95         GAMSUM=GAMSAV(1)+GAMSAV(2)
  96         JSAV(2)=JTYPE(J)
  97 C          W+- and Z0
  98         IF(AM0.GT.WMASS(4)) THEN
  99           DO 200 IW=2,4
 100             GAMSAV(IW+1)=0.
 101             ZCSAV(IW+1)=.5
 102             JSAV(IW+1)=0
 103             JT0=2*IABS(JTYPE(J))
 104             IF(JTYPE(J).LT.0) JT0=JT0+1
 105             JT1=MATCH(JT0,IW)
 106             IF(JT1.EQ.0) GO TO 200
 107             JT1=MATCH(JT1,4)
 108             IFL1=JT1/2
 109             AM1W=AMASS(IFL1)
 110             AM2W=AMASS(LISTW(IW+1))
 111             IF(AM1W+AM2W.GE.AM0) GO TO 200
 112             T0=AM0**2
 113             T1W=AM1W**2
 114             T2W=AM2W**2
 115             ZCW=(T0-T1W+T2W-SQRT((T0-T1W-T2W)**2-4*T1W*T2W))/(2*T0)
 116             GAMZCW=2.*ALOG((1-ZCW)/ZCW)-2.*(1.-2.*ZCW)
 117             TERM=(AQ(JTABS,IW)**2+BQ(JTABS,IW)**2)*ALFA*GAMZCW
 118             GAMSAV(IW+1)=TERM
 119             ZCSAV(IW+1)=ZCW
 120             JSAV(IW+1)=IFL1*ISIGN(1,JTYPE(J))
 121             GAMSUM=GAMSUM+TERM
 122 200       CONTINUE
 123         ELSE
 124           DO 210 IW=2,4
 125             GAMSAV(IW+1)=0.
 126             ZCSAV(IW+1)=.5
 127             JSAV(IW+1)=0
 128 210       CONTINUE
 129         ENDIF
 130 C          Select decay mode
 131         RND=RANF()
 132         SUM=0.
 133         DO 220 IW=1,5
 134           IWTYPE=IW
 135           SUM=SUM+GAMSAV(IW)/GAMSUM
 136           IF(RND.LT.SUM) GO TO 230
 137 220     CONTINUE
 138 C          Generate Z
 139 230     CONTINUE
 140         Z=1-(ZC/(1-ZC))**RANF()*(1-ZC)
 141         GZ=8./(3.*(1-Z))
 142         IF(PQQ(Z).LT.GZ*RANF()) GO TO 230
 143         IF(Z.LT.ZCSAV(IWTYPE).OR.Z.GT.1.-ZCSAV(IWTYPE)) GO TO 230
 144         JTYPE(NJ1)=JSAV(IWTYPE)
 145         JTYPE(NJ2)=LISTW(IWTYPE)
 146         ZZC(J)=Z
 147 C
 148 C          Initial diquark
 149 C
 150       ELSEIF(MOD(JTABS,100).EQ.0) THEN
 151 300     CONTINUE
 152         Z=1-(ZC/(1-ZC))**RANF()*(1-ZC)
 153         GZ=8./(3.*(1-Z))
 154         IF(PQQ(Z).LT.GZ*RANF()) GO TO 300
 155         JTYPE(NJ1)=JTYPE(J)
 156         JTYPE(NJ2)=9
 157         ZZC(J)=Z
 158 C
 159 C          Initial gluino
 160 C
 161        ELSEIF (JTABS.EQ.29) THEN
 162 400      Z=1.-(ZC/(1.-ZC))**RANF()*(1.-ZC)
 163          GZ=6./(1.-Z)
 164          IF(PGSGS(Z) .LT. GZ*RANF()) GOTO 400
 165          JTYPE(NJ1)=JTYPE(J)
 166          JTYPE(NJ2)=9
 167          ZZC(J)=Z
 168 C
 169 C          Initial squark
 170 C
 171       ELSEIF(JTABS.GT.20.AND.JTABS.LT.29) THEN
 172 500     CONTINUE
 173         Z=1-(ZC/(1-ZC))**RANF()*(1-ZC)
 174         GZ=8./(3.*(1-Z))
 175         IF(PQSQS(Z).LT.GZ*RANF()) GO TO 500
 176         JTYPE(NJ1)=JTYPE(J)
 177         JTYPE(NJ2)=9
 178         ZZC(J)=Z
 179 C
 180 C          Initial W+, W-, or Z0
 181 C
 182       ELSEIF(JTABS.EQ.80.OR.JTABS.EQ.90) THEN
 183 C          Select decay mode
 184         IF(JTYPE(J).EQ.+80) JW=2
 185         IF(JTYPE(J).EQ.-80) JW=3
 186         IF(JTYPE(J).EQ.+90) JW=4
 187         TRY=RANF()
 188         DO 610 IQ=2,25
 189         IF(TRY.LT.CUMWBR(IQ,JW-1)) THEN
 190           IQ1=IQ
 191           IQ2=MATCH(IQ,JW)
 192           GO TO 620
 193         ENDIF
 194 610     CONTINUE
 195 620     JTYPE(NJ1)=LISTJ(IQ1)
 196         JTYPE(NJ2)=LISTJ(IQ2)
 197 C          Select W helicity
 198         JPAR=MOD(JORIG(J),JPACK)
 199         IFLPAR=IABS(JTYPE(JPAR))
 200         HELPL=(AQ(IFLPAR,JW)-BQ(IFLPAR,JW))**2
 201         HELMN=(AQ(IFLPAR,JW)+BQ(IFLPAR,JW))**2
 202         IF(RANF()*(HELPL+HELMN).LT.HELMN) THEN
 203           HEL=-ISIGN(1,JTYPE(NJ1))
 204         ELSE
 205           HEL=+ISIGN(1,JTYPE(NJ1))
 206         ENDIF
 207 630     Z=RANF()
 208         PZ=(1.+HEL*(2.*Z-1.))**2
 209         IF(PZ.LT.4.*RANF()) GO TO 630
 210         ZZC(J)=Z
 211       ENDIF
 212 C
 213 C          Set masses and flags.
 214 C
 215       JET=IABS(JORIG(J))/JPACK
 216       JORIG(NJ1)=JPACK*JET+J
 217       JORIG(NJ2)=JPACK*JET+J
 218       IDABS1=IABS(JTYPE(NJ1))
 219       IDABS2=IABS(JTYPE(NJ2))
 220       JMATCH(NJ1)=NJ2
 221       JMATCH(NJ2)=NJ1
 222 C          JDCAY=-1 implies further decay
 223       IF(IDABS1.LE.9.OR.(IDABS1.GT.20.AND.IDABS1.LT.30.).OR.
 224      $MOD(IDABS1,100).EQ.0) THEN
 225         PJSET(5,NJ1)=Z*AM0
 226         JDCAY(NJ1)=-1
 227       ELSEIF(IDABS1.GE.80.OR.IDABS1.LE.90) THEN
 228         PJSET(5,NJ1)=AMASS(IDABS1)
 229         JDCAY(NJ1)=-1
 230       ELSE
 231         PJSET(5,NJ1)=AMASS(IDABS1)
 232         JDCAY(NJ1)=0
 233       ENDIF
 234       IF(IDABS2.LE.9.OR.(IDABS2.GT.20.AND.IDABS2.LT.30.).OR.
 235      $MOD(IDABS2,100).EQ.0) THEN
 236         PJSET(5,NJ2)=(1.-Z)*AM0
 237         JDCAY(NJ2)=-1
 238       ELSEIF(IDABS2.EQ.80.OR.IDABS2.EQ.90) THEN
 239         PJSET(5,NJ2)=AMASS(IDABS2)
 240         JDCAY(NJ2)=-1
 241       ELSE
 242         PJSET(5,NJ2)=AMASS(IDABS2)
 243         JDCAY(NJ2)=0
 244       ENDIF
 245       RETURN
 246       END