ISAJET/code/hevolv.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE HEVOLV
   3 C
   4 C          CARRY OUT BACKWARDS EVOLUTION QK --> QK + W FOR LONGITUDINAL
   5 C          W-W FUSION, GENERATING Z AND KT**2 FROM RELATION OF W AND
   6 C          QUARK STRUCTURE FUNCTIONS.
   7 C
   8 #include "isajet/itapes.inc"
   9 #include "isajet/qcdpar.inc"
  10 #include "isajet/jetpar.inc"
  11 #include "isajet/pjets.inc"
  12 #include "isajet/jetset.inc"
  13 #include "isajet/primar.inc"
  14 #include "isajet/wcon.inc"
  15 #include "isajet/const.inc"
  16 #include "isajet/idrun.inc"
  17 #include "isajet/hcon.inc"
  18 C
  19       DIMENSION X(2)
  20       EQUIVALENCE (X1,X(1))
  21       DIMENSION FZIQ(13),IWPICK(2),PFINAL(5),BST1(5),BST2(5),B2B1(5)
  22       DIMENSION PSAVE(5,2)
  23 C          LAMBDA FUNCTION
  24       ALAMF(A,B,C)=SQRT((A-B-C)**2-4.*B*C)
  25 C
  26       NJSAVE=NJSET
  27       NREJ2=-1
  28 C
  29 C          INITIALIZE
  30       DO 10 I=1,2
  31       DO 10 K=1,5
  32 10    PSAVE(K,I)=PJSET(K,I)
  33 20    CONTINUE
  34       DO 30 I=1,2
  35       DO 30 K=1,5
  36 30    PJSET(K,I)=PSAVE(K,I)
  37       DO 40 K=1,5
  38 40    PFINAL(K)=QWJET(K)
  39       NJSET=NJSAVE
  40 C
  41 C          CHOOSE A W AND DO BACKWARDS EVOLUTION FOR QK -> QK + W.
  42 C
  43       IF(RANF().LT..5) THEN
  44         IWPICK(1)=1
  45         IWPICK(2)=2
  46         SGN=+1.
  47       ELSE
  48         IWPICK(1)=2
  49         IWPICK(2)=1
  50         SGN=-1.
  51       ENDIF
  52       DO 100 JJ=1,2
  53 C
  54 C          OTHER PARTICLE IS W FOR JJ=1, QUARK FOR JJ=2:
  55       IF(JJ.EQ.1) THEN
  56         J1=IWPICK(1)
  57         J2=IWPICK(2)
  58       ELSE
  59         J1=IWPICK(2)
  60         J2=NJSAVE+1
  61         SGN=-SGN
  62       ENDIF
  63       JTLV1=JTYPE(J1)
  64       IF(JTLV1.EQ.10) THEN
  65         IW=1
  66       ELSEIF(JTLV1.EQ.80) THEN
  67         IW=2
  68       ELSEIF(JTLV1.EQ.-80) THEN
  69         IW=3
  70       ELSEIF(JTLV1.EQ.90) THEN
  71         IW=4
  72       ENDIF
  73       XV=(PJSET(4,J1)+ABS(PJSET(3,J1)))/ECM
  74       AMV=AMASS(JTLV1)
  75 C
  76 C          GENERATE VARIABLES FOR BRANCHING
  77 C          FIND MAXIMUM OF INTEGRAND USING 20 POINTS IN LOG(Z)
  78       FMAX=0.
  79       ZMULT=(1./XV)**.05
  80       ZIZ=XV
  81       DO 110 IZ=1,19
  82       ZIZ=ZIZ*ZMULT
  83       FSUM=0.
  84       DO 115 IQ=2,13
  85       IF(MATCH(IQ,IW).NE.0) THEN
  86         IFL=IQ/2
  87         CIQ=AQ(IFL,IW)**2+BQ(IFL,IW)**2
  88         FSUM=FSUM+CIQ*(1.-ZIZ)/ZIZ*STRUC(XV/ZIZ,AMV**2,IQ,IDIN(J1))
  89       ENDIF
  90 115   CONTINUE
  91       FMAX=AMAX1(FMAX,FSUM)
  92 110   CONTINUE
  93 C          GENERATE Z UNIFORMLY IN (XV,1) AND TEST
  94       NREJ1=-1
  95 120   ZV=XV+(1.-XV)*RANF()
  96       FZ=0.
  97       DO 130 IQ=2,13
  98       IF(MATCH(IQ,IW).NE.0) THEN
  99         IFL=IQ/2
 100         CIQ=AQ(IFL,IW)**2+BQ(IFL,IW)**2
 101         FZIQ(IQ)=CIQ*(1.-ZV)/ZV*STRUC(XV/ZV,AMV**2,IQ,IDIN(J1))
 102       ELSE
 103         FZIQ(IQ)=0.
 104       ENDIF
 105 130   FZ=FZ+FZIQ(IQ)
 106       IF(FZ.LT.FMAX*RANF()) THEN
 107         NREJ1=NREJ1+1
 108         IF(NREJ1.GT.NTRIES) GO TO 9999
 109         GO TO 120
 110       ENDIF
 111 C          DETERMINE QUARK TYPE
 112       TRY=RANF()
 113       SUM=0.
 114       DO 140 IQ=2,13
 115       IQ1=IQ
 116       SUM=SUM+FZIQ(IQ)/FZ
 117 140   IF(SUM.GT.TRY) GO TO 150
 118 150   IQ3=MATCH(IQ1,IW)
 119       IQ3=MATCH(IQ3,4)
 120 C          GENERATE T=-K**2 AND UNIFORM PHI
 121       T=AMV**2*(1./RANF()-1.)
 122       PHIK=2.*PI*RANF()
 123 C
 124 C          SOLVE KINEMATICS FOR THIS SIDE
 125       S=(PJSET(4,J1)+PJSET(4,J2))**2-(PJSET(1,J1)+PJSET(1,J2))**2
 126      $-(PJSET(2,J1)+PJSET(2,J2))**2-(PJSET(3,J1)+PJSET(3,J2))**2
 127       SP=S/ZV
 128       IFL1=IQ1/2
 129       IFL2=JTYPE(J2)
 130       IFL3=IQ3/2
 131       AM1=AMASS(IFL1)
 132       AM2=PJSET(5,J2)
 133       AM3=AMASS(IFL3)
 134       AM1SQ=AM1**2
 135       AM2SQ=AM2**2
 136       AM3SQ=AM3**2
 137       IF(SGN.LT.0.) THEN
 138         P2PL=PJSET(4,J2)+PJSET(3,J2)
 139         P2MN=AM2SQ/P2PL
 140       ELSE
 141         P2MN=PJSET(4,J2)-PJSET(3,J2)
 142         P2PL=AM2SQ/P2MN
 143       ENDIF
 144 C          STEP 1: SOLVE FOR PP1=PJSET(K,NEWV)
 145       IF(SGN.GT.0.) THEN
 146         PP1PL=(SP-AM1SQ-AM2SQ+ALAMF(SP,AM1SQ,AM2SQ))/(2.*P2MN)
 147         PP1MN=AM1SQ/PP1PL
 148       ELSE
 149         PP1MN=(SP-AM1SQ-AM2SQ+ALAMF(SP,AM1SQ,AM2SQ))/(2.*P2PL)
 150         PP1PL=AM1SQ/PP1MN
 151       ENDIF
 152 C          STEP 2: SOLVE FOR K = VIRTUAL W MOMENTUM
 153       DEN=PP1PL*P2MN-PP1MN*P2PL
 154       AKPL=(+PP1PL*(S+T-AM2SQ)+P2PL*(T+AM3SQ-AM1SQ))/DEN
 155       AKMN=(-PP1MN*(S+T-AM2SQ)-P2MN*(T+AM3SQ-AM1SQ))/DEN
 156       WPL=PP1PL-AKPL
 157       WMN=PP1MN-AKMN
 158       AKT2=T+AKPL*AKMN
 159 C          STEP 3: START OVER IF AKT2 UNPHYSICAL
 160       IF(AKT2.LE.0..OR.PP1PL.GE.ECM.OR.PP1MN.GE.ECM.OR.
 161      $P2PL.GE.ECM.OR.P2MN.GE.ECM) THEN
 162         NREJ2=NREJ2+1
 163         IF(NREJ2.GT.NTRIES) GO TO 9999
 164         GO TO 20
 165       ENDIF
 166 C
 167 C          SAVE NEW VECTORS
 168       NJ1=NJSET+1
 169       NJ2=NJSET+2
 170       AKT=SQRT(AKT2)
 171       AKX=AKT*COS(PHIK)
 172       AKY=AKT*SIN(PHIK)
 173       PJSET(1,J1)=AKX
 174       PJSET(2,J1)=AKY
 175       PJSET(3,J1)=.5*(AKPL-AKMN)
 176       PJSET(4,J1)=.5*(AKPL+AKMN)
 177       PJSET(5,J1)=-SQRT(T)
 178       JDCAY(J1)=JPACK*NJ1+NJ2
 179       JET=IABS(JORIG(J1))/JPACK
 180 C
 181       PJSET(1,NJ1)=0.
 182       PJSET(2,NJ1)=0.
 183       PJSET(3,NJ1)=.5*(PP1PL-PP1MN)
 184       PJSET(4,NJ1)=.5*(PP1PL+PP1MN)
 185       PJSET(5,NJ1)=AM1
 186       JORIG(NJ1)=JPACK*JET+J1
 187       JTYPE(NJ1)=IFL1
 188       JDCAY(NJ1)=0
 189 C
 190       PJSET(1,NJ2)=-AKX
 191       PJSET(2,NJ2)=-AKY
 192       PJSET(3,NJ2)=.5*(WPL-WMN)
 193       PJSET(4,NJ2)=.5*(WPL+WMN)
 194       PJSET(5,NJ2)=AM3
 195       JORIG(NJ2)=JPACK*JET+J1
 196       JTYPE(NJ2)=IFL3
 197       JDCAY(NJ2)=0
 198 C
 199 C          BOOST OTHER VECTORS TO NEW FRAME GIVEN BY DIFFERENCE OF
 200 C          OLD AND NEW FINAL MOMENTA.
 201       DO 200 K=1,4
 202       BST1(K)=PFINAL(K)
 203 200   BST2(K)=PJSET(K,J1)+PJSET(K,J2)
 204       BMASS=PFINAL(5)
 205       BST1(5)=BMASS
 206       BST2(5)=BMASS
 207 C
 208 C          PARAMETERS FOR COMBINED BOOSTS.
 209       BDOTB=BST1(4)*BST2(4)-BST1(1)*BST2(1)-BST1(2)*BST2(2)
 210      $-BST1(3)*BST2(3)
 211       DO 210 K=1,4
 212 210   B2B1(K)=BST2(K)-BST1(K)
 213 C
 214       B44=BDOTB/BMASS**2
 215       BI41=1./BMASS
 216       BI42=(BDOTB-BMASS**2-B2B1(4)*BMASS)/(BMASS**2*(BST2(4)+BMASS))
 217       B4K1=BI41
 218       B4K2=(BMASS**2-BDOTB-B2B1(4)*BMASS)/(BMASS**2*(BST1(4)+BMASS))
 219       BIK1=-1./(BMASS*(BST1(4)+BMASS))
 220       BIK2=1./(BMASS*(BST2(4)+BMASS))
 221       BIK3=(BMASS**2-BDOTB)/(BMASS**2*(BST1(4)+BMASS)
 222      $*(BST2(4)+BMASS))
 223 C
 224 C          BOOST FINAL JETS
 225       DO 220 J=1,NJSET
 226       IF(J.EQ.J1.OR.J.EQ.J2) GO TO 220
 227       IF(PJSET(5,J).LT.0.) GO TO 220
 228       BP1=0.
 229       BP21=0.
 230       DO 221 K=1,3
 231       BP1=BP1+BST1(K)*PJSET(K,J)
 232 221   BP21=BP21+B2B1(K)*PJSET(K,J)
 233       DO 222 K=1,3
 234 222   PJSET(K,J)=PJSET(K,J)
 235      $+(B2B1(K)*BI41+BST2(K)*BI42)*PJSET(4,J)
 236      $+B2B1(K)*BP1*BIK1+BST2(K)*BP21*BIK2+BST2(K)*BP1*BIK3
 237       PJSET(4,J)=B44*PJSET(4,J)+BP21*B4K1+BP1*B4K2
 238 220   CONTINUE
 239 C
 240 C          RESET VIRTUAL MOMENTA
 241       DO 230 J=1,NJSET
 242       IF(J.EQ.J1.OR.J.EQ.J2) GO TO 230
 243       IF(PJSET(5,J).GE.0.) GO TO 230
 244       JX1=JDCAY(J)/JPACK
 245       JX2=JDCAY(J)-JPACK*JX1
 246       DO 231 K=1,4
 247 231   PJSET(K,J)=PJSET(K,JX1)-PJSET(K,JX2)
 248       AMJ=PJSET(4,J)**2-PJSET(1,J)**2-PJSET(2,J)**2-PJSET(3,J)**2
 249       PJSET(5,J)=-SQRT(ABS(AMJ))
 250 230   CONTINUE
 251 C
 252 C          RESET PFINAL AND NJSET
 253       DO 240 K=1,4
 254 240   PFINAL(K)=PJSET(K,J2)+PJSET(K,NJ1)
 255       PFINAL(5)=SQRT(SP)
 256       NJSET=NJSET+2
 257 100   CONTINUE
 258       RETURN
 259 C
 260 9999  CONTINUE
 261       WRITE(ITLIS,9998) IEVT
 262 9998  FORMAT(/' ***** ERROR IN HEVOLV ... EVENT',I8,' DISCARDED *****')
 263       NJSET=-1
 264       RETURN
 265       END