ISAJET/code/qcdini.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE QCDINI(JIN1,JIN2)
   3 C
   4 C          GENERATE INITIAL-STATE QCD CASCADE USING BACKWARDS
   5 C          EVOLUTION OF GOTTSCHALK AND OF SJOSTRAND.
   6 C
   7 C          IF QCDINI FAILS WHEN ATTEMPTING TO FORCE GL-->QK+QB FOR
   8 C          HEAVY QUARKS, THEN RETURN NJSET=-1.
   9 C
  10 C          VER. 6.40: TRAP W1LIM > 0 TO PREVENT ROUNDING ERRORS.
  11 C
  12 #include "isajet/itapes.inc"
  13 #include "isajet/idrun.inc"
  14 #include "isajet/pinits.inc"
  15 #include "isajet/jetpar.inc"
  16 #include "isajet/qcdpar.inc"
  17 #include "isajet/jetset.inc"
  18 #include "isajet/jwork.inc"
  19 #include "isajet/jwork2.inc"
  20 #include "isajet/const.inc"
  21 #include "isajet/primar.inc"
  22 #include "isajet/keys.inc"
  23 C
  24       DIMENSION BOOST1(5),BOOST2(5),B2B1(5),DBL1(5),DBL2(5)
  25       DIMENSION FXOLD(2),FXNEW(2)
  26       DIMENSION PJKEEP(5,12),JINS(2),JLIST(16),PFKEEP(5)
  27 #if defined(CERNLIB_DOUBLE)
  28       DOUBLE PRECISION DBL1,DBL2,DBLM
  29 #endif
  30 C
  31 C          CONVERT IDENT+7 TO JETTYP
  32       DATA JLIST/13,11,9,7,5,3,0,2,4,6,8,10,12,0,0,1/
  33       ALAMF(A,B,C)=SQRT((A-B-C)**2-4.*B*C)
  34 C
  35 C          INITIALIZE
  36 C
  37       JINS(1)=JIN1
  38       JINS(2)=JIN2
  39       DO 97 K=1,4
  40 97    PFKEEP(K)=PJSET(K,JIN1)+PJSET(K,JIN2)
  41 C          EXCEPT FOR HIGGS, PFKEEP**2=SHAT
  42       IF(KEYS(7).OR.KEYS(9)) THEN
  43         S1KEEP=PFKEEP(4)**2-PFKEEP(1)**2-PFKEEP(2)**2-PFKEEP(3)**2
  44         PFKEEP(5)=SQRT(S1KEEP)
  45         PPKEEP=PFKEEP(4)+PFKEEP(3)
  46         PMKEEP=PFKEEP(4)-PFKEEP(3)
  47       ELSE
  48         S1KEEP=SHAT
  49         PFKEEP(5)=SQRT(S1KEEP)
  50         IF(PFKEEP(3).GT.0.) THEN
  51           PPKEEP=PFKEEP(4)+PFKEEP(3)
  52           PMKEEP=(S1KEEP+PFKEEP(1)**2+PFKEEP(2)**2)/PPKEEP
  53         ELSE
  54           PMKEEP=PFKEEP(4)-PFKEEP(3)
  55           PPKEEP=(S1KEEP+PFKEEP(1)**2+PFKEEP(2)**2)/PMKEEP
  56         ENDIF
  57         PFKEEP(4)=.5*(PPKEEP+PMKEEP)
  58         PFKEEP(3)=.5*(PPKEEP-PMKEEP)
  59       ENDIF
  60       DO 98 I=1,NJSET
  61       DO 98 K=1,5
  62 98    PJKEEP(K,I)=PJSET(K,I)
  63       NJKEEP=NJSET
  64       NPASS=0
  65       NPASS1=0
  66 C
  67 1     CONTINUE
  68       NPASS1=NPASS1+1
  69       IF(NPASS1.GT.100) GO TO 9999
  70       NJSET=NJKEEP
  71       DO 99 I=1,NJSET
  72       DO 99 K=1,5
  73 99    PJSET(K,I)=PJKEEP(K,I)
  74 C
  75       DO 100 K=1,5
  76 100   PFINAL(K)=PFKEEP(K)
  77       S1=S1KEEP
  78       PTOTPL=PPKEEP
  79       PTOTMN=PMKEEP
  80       TCUT=CUTJET**2
  81       DO 101 I=1,2
  82       JI=JINS(I)
  83       XOLD=(PJSET(4,JI)+ABS(PJSET(3,JI)))/ECM
  84       JT=JLIST(JTYPE(JI)+7)
  85       FXOLD(I)=STRUC(XOLD,QSQ,JT,IDIN(I))
  86 101   CONTINUE
  87 C
  88 C          DO FIRST EVOLUTION
  89       DO 110 I=1,2
  90       SGN=3-2*I
  91       JET=10+I
  92       JI=JINS(I)
  93       ZMIN=(PJSET(4,JI)+ABS(PJSET(3,JI)))/ECM
  94       ZMAX=1./(1.+TCUT/S1)
  95 C          DZMAX=1.-ZMAX
  96       DZMAX=ZMAX*TCUT/S1
  97       IF(ZMIN.GE.ZMAX) ZMIN=.5*ZMAX
  98       CALL QCDINT(JI)
  99       JVIR(I)=JI
 100 110   CONTINUE
 101 C
 102 C          SOLVE INITIAL KINEMATICS
 103       AM1SQ=PJSET(5,JVIR(1))**2*SIGN(1.,PJSET(5,JVIR(1)))
 104       AM2SQ=PJSET(5,JVIR(2))**2*SIGN(1.,PJSET(5,JVIR(2)))
 105       P1PL=(S1+AM1SQ-AM2SQ+ALAMF(S1,AM1SQ,AM2SQ))/(2.*PTOTMN)
 106       P1MN=AM1SQ/P1PL
 107       P2MN=(S1+AM2SQ-AM1SQ+ALAMF(S1,AM1SQ,AM2SQ))/(2.*PTOTPL)
 108       P2PL=AM2SQ/P2MN
 109       PJSET(3,JVIR(1))=.5*(P1PL-P1MN)
 110       PJSET(4,JVIR(1))=.5*(P1PL+P1MN)
 111       PJSET(3,JVIR(2))=.5*(P2PL-P2MN)
 112       PJSET(4,JVIR(2))=.5*(P2PL+P2MN)
 113 C
 114 C          TEST WHETHER NEW MASS IS PLAUSIBLE
 115       DO 111 I=1,2
 116       JI=JINS(I)
 117       XNEW=(PJSET(4,JI)+ABS(PJSET(3,JI)))/ECM
 118       IF(XNEW.GE.1.) THEN
 119         FXNEW(I)=0.
 120       ELSE
 121         JT=JLIST(JTYPE(JI)+7)
 122         FXNEW(I)=STRUC(XNEW,QSQ,JT,IDIN(I))
 123       ENDIF
 124 111   CONTINUE
 125       DO 112 I=1,2
 126       IF(FXNEW(I).LT.FXOLD(I)*RANF()) GO TO 1
 127 112   CONTINUE
 128 C
 129 C          FIND JVIR (SPACE-LIKE PARTON) WITH LARGER (-MASS) FOR NEXT
 130 C          BRANCHING.
 131 10    IF(JDCAY(JVIR(1)).GE.0.AND.JDCAY(JVIR(2)).GE.0) RETURN
 132       NPASS=NPASS+1
 133       IF(NPASS.GT.20*NJSET) GO TO 9999
 134       IF(-PJSET(5,JVIR(1)).GE.-PJSET(5,JVIR(2))) THEN
 135         IVIR=JVIR(1)
 136         IVIR2=JVIR(2)
 137         SGN=+1.
 138         JET=11
 139       ELSE
 140         IVIR=JVIR(2)
 141         IVIR2=JVIR(1)
 142         SGN=-1.
 143         JET=12
 144       ENDIF
 145 C
 146       T1=PJSET(5,IVIR)**2
 147       ZMIN=(PJSET(4,IVIR)+SGN*PJSET(3,IVIR))/ECM
 148       ZMAX=1./(1.+T1/S1)
 149       DZMAX=ZMAX*T1/S1
 150       IF(ZMIN.GE.ZMAX) GO TO 1
 151 C
 152 C          GENERATE Z AND NEW PARTONS.
 153 C          NEWV=SPACELIKE, NEWF=TIMELIKE.
 154       NEWV=NJSET+1
 155       NEWF=NJSET+2
 156       CALL QCDINZ(IVIR)
 157 C
 158 C          IF Z FAILS (BECAUSE OF STRUCTURE FUNCTION) SET NEWV=IVIR,
 159 C          NEWF=NULL AND RE-SOLVE KINEMATICS.
 160 15    IF(.NOT.ZGOOD) THEN
 161         CALL QCDINT(IVIR)
 162 C
 163         PP1PL=PJSET(4,IVIR2)+PJSET(3,IVIR2)
 164         PP1MN=PJSET(4,IVIR2)-PJSET(3,IVIR2)
 165         AMSQ=PJSET(5,IVIR)**2*SIGN(1.,PJSET(5,IVIR))
 166         AMPSQ=PJSET(5,IVIR2)**2*SIGN(1.,PJSET(5,IVIR2))
 167         IF(SGN.GT.0) THEN
 168           P2PL=(S1-AMSQ-AMPSQ+ALAMF(S1,AMSQ,AMPSQ))/(2.*PP1MN)
 169           P2MN=AMSQ/P2PL
 170         ELSE
 171           P2MN=(S1-AMSQ-AMPSQ+ALAMF(S1,AMSQ,AMPSQ))/(2.*PP1PL)
 172           P2PL=AMSQ/P2MN
 173         ENDIF
 174         PJSET(3,IVIR)=.5*(P2PL-P2MN)
 175         PJSET(4,IVIR)=.5*(P2PL+P2MN)
 176 C
 177         NEWV=IVIR
 178         DO 120 K=1,5
 179 120     PJSET(K,NEWF)=0.
 180         GO TO 30
 181       ENDIF
 182 C
 183 C          EVOLVE NEW SPACELIKE PARTON.
 184       PJSET(5,NEWV)=PJSET(5,IVIR)
 185       S2=S1/ZZC(IVIR)
 186       ZMIN=ZMIN/ZZC(IVIR)
 187       ZMAX=1./(1.+TCUT/S2)
 188       DZMAX=ZMAX*TCUT/S2
 189       IF(ZMIN.GE.ZMAX) GO TO 1
 190       CALL QCDINT(NEWV)
 191 C
 192 C          CALCULATE APPROXIMATE MASS LIMIT AND DO TIMELIKE EVOLUTION.
 193 C          VER. 6.40: TRAP W1LIM < 0 FROM ROUNDING ERRORS.
 194       W1LIM=T1*(1./(ZZC(IVIR)*(1.+T1/S1))-1.)
 195       W1LIM=AMIN1(W1LIM,T1)
 196       PJSET(5,NEWF)=SQRT(ABS(W1LIM))
 197       JDCAY(NEWF)=-1
 198 20    CALL QCDT(NEWF)
 199 C
 200 C          SOLVE KINEMATICS USING +(PL) AND -(MN) COMPONENTS FOR
 201 C          PJSET(K,NEWV)+PJSET(K,IVIR2)-->PJSET(K,NEWF)+PFINAL
 202 C          STEP 1: SOLVE FOR P2=PJSET(K,NEWV)
 203       PP1PL=PJSET(4,IVIR2)+PJSET(3,IVIR2)
 204       PP1MN=PJSET(4,IVIR2)-PJSET(3,IVIR2)
 205       AMSQ=PJSET(5,NEWV)**2*SIGN(1.,PJSET(5,NEWV))
 206       AMPSQ=PJSET(5,IVIR2)**2*SIGN(1.,PJSET(5,IVIR2))
 207       W1=PJSET(5,NEWF)**2
 208       IF(SGN.GT.0) THEN
 209         P2PL=(S2-AMSQ-AMPSQ+ALAMF(S2,AMSQ,AMPSQ))/(2.*PP1MN)
 210         P2MN=AMSQ/P2PL
 211       ELSE
 212         P2MN=(S2-AMSQ-AMPSQ+ALAMF(S2,AMSQ,AMPSQ))/(2.*PP1PL)
 213         P2PL=AMSQ/P2MN
 214       ENDIF
 215 C
 216 C          STEP 2: SOLVE FOR Q1(K)=PJSET(K,IVIR)
 217       DEN=P2PL*PP1MN-P2MN*PP1PL
 218       Q1PL=(+P2PL*(S1+T1-AMPSQ)+PP1PL*(W1+T1-AMSQ))/DEN
 219       Q1MN=(-P2MN*(S1+T1-AMPSQ)-PP1MN*(W1+T1-AMSQ))/DEN
 220       WPL=P2PL-Q1PL
 221       WMN=P2MN-Q1MN
 222 C          CALCULATE TRANSVERSE MOMENTUM AND REJECT IF UNPHYSICAL.
 223       Q1TR2=T1+Q1PL*Q1MN
 224       IF(Q1TR2.LT.0.) THEN
 225         IF(JDCAY(NEWF).EQ.-1) GO TO 20
 226         ZGOOD=.FALSE.
 227         GO TO 15
 228       ENDIF
 229 C
 230 C          DO ONE TIMELIKE BRANCHING TO INSURE CORRECT MASS. MUST FIRST
 231 C          SHIFT NJSET TO PUT DECAY PRODUCTS IN CORRECT PLACE.
 232       IF(JDCAY(NEWF).EQ.-1) THEN
 233         NJSET=NJSET+2
 234         CALL QCDZ(NEWF)
 235         NJSET=NJSET-2
 236         Z1=ZZC(NEWF)
 237         E0=.5*(WPL+WMN)
 238         P0=SQRT(.25*(WPL-WMN)**2+Q1TR2)
 239         WM0=PJSET(5,NEWF)
 240         ZLIM=AMAX1((WM0/(E0+P0))**2,CUTJET/(E0+P0))
 241         IF(Z1.LE.ZLIM.OR.Z1.GE.1.-ZLIM) GO TO 20
 242         NEWF1=NEWF+1
 243         NEWF2=NEWF+2
 244         JDCAY(NEWF)=NEWF1*JPACK+NEWF2
 245         CALL QCDT(NEWF1)
 246         CALL QCDT(NEWF2)
 247         JORIG(NEWF1)=JPACK*JET+NEWF
 248         JORIG(NEWF2)=JORIG(NEWF1)
 249         DO 130 K=1,4
 250         PJSET(K,NEWF1)=0.
 251 130     PJSET(K,NEWF2)=0.
 252       ENDIF
 253 C
 254 C          GOOD BRANCHING!
 255       PHIQ1=2.*PI*RANF()
 256       Q1TR=SQRT(Q1TR2)
 257
 258       Q1X=Q1TR*COS(PHIQ1)
 259       Q1Y=Q1TR*SIN(PHIQ1)
 260 C
 261       PJSET(1,IVIR)=Q1X
 262       PJSET(2,IVIR)=Q1Y
 263       PJSET(3,IVIR)=.5*(Q1PL-Q1MN)
 264       PJSET(4,IVIR)=.5*(Q1PL+Q1MN)
 265       JDCAY(IVIR)=JPACK*NEWV+NEWF
 266 C
 267       PJSET(1,NEWV)=0.
 268       PJSET(2,NEWV)=0.
 269       PJSET(3,NEWV)=.5*(P2PL-P2MN)
 270       PJSET(4,NEWV)=.5*(P2PL+P2MN)
 271       JORIG(NEWV)=JPACK*JET+IVIR
 272 C
 273       PJSET(1,NEWF)=-Q1X
 274       PJSET(2,NEWF)=-Q1Y
 275       PJSET(3,NEWF)=.5*(WPL-WMN)
 276       PJSET(4,NEWF)=.5*(WPL+WMN)
 277       JORIG(NEWF)=JPACK*JET+IVIR
 278 C
 279 C          BOOST ALL FINAL VECTORS (EXCEPT NEW ONES) AND RECALCULATE
 280 C          VIRTUAL MOMENTA.  BOOST IS DETERMINED BY DIFFERENCE OF
 281 C          NEW AND OLD TOTAL FINAL MOMENTA, B2B1=BOOST2-BOOST1.
 282 C
 283 30    CONTINUE
 284       DO 201 K=1,4
 285 201   BOOST1(K)=PFINAL(K)
 286       BMASS=PFINAL(5)
 287       DO 202 K=1,4
 288 202   BOOST2(K)=PJSET(K,NEWV)+PJSET(K,IVIR2)-PJSET(K,NEWF)
 289 C
 290 C          PARAMETERS FOR COMBINED BOOSTS.
 291 #if defined(CERNLIB_SINGLE)
 292       BDOTB=BOOST1(4)*BOOST2(4)-BOOST1(1)*BOOST2(1)-BOOST1(2)*BOOST2(2)
 293      $-BOOST1(3)*BOOST2(3)
 294       DO 203 K=1,4
 295 203   B2B1(K)=BOOST2(K)-BOOST1(K)
 296 #endif
 297 #if defined(CERNLIB_DOUBLE)
 298 C          DOUBLE PRECISION FOR 32-BIT MACHINES USING 3-VECTORS AND MASS
 299 C          AS EXACT.
 300       DO 204 K=1,3
 301       DBL1(K)=BOOST1(K)
 302 204   DBL2(K)=BOOST2(K)
 303       DBLM=BMASS
 304       DBL1(4)=DSQRT(DBL1(1)**2+DBL1(2)**2+DBL1(3)**2+DBLM**2)
 305       DBL2(4)=DSQRT(DBL2(1)**2+DBL2(2)**2+DBL2(3)**2+DBLM**2)
 306       BDOTB=DBL1(4)*DBL2(4)-DBL1(1)*DBL2(1)-DBL1(2)*DBL2(2)
 307      $-DBL1(3)*DBL2(3)
 308       DO 205 K=1,4
 309 205   B2B1(K)=DBL2(K)-DBL1(K)
 310 #endif
 311       B44=BDOTB/BMASS**2
 312       BI41=1./BMASS
 313       BI42=(BDOTB-BMASS**2-B2B1(4)*BMASS)/(BMASS**2*(BOOST2(4)+BMASS))
 314       B4K1=BI41
 315       B4K2=(BMASS**2-BDOTB-B2B1(4)*BMASS)/(BMASS**2*(BOOST1(4)+BMASS))
 316       BIK1=-1./(BMASS*(BOOST1(4)+BMASS))
 317       BIK2=1./(BMASS*(BOOST2(4)+BMASS))
 318       BIK3=(BMASS**2-BDOTB)/(BMASS**2*(BOOST1(4)+BMASS)
 319      $*(BOOST2(4)+BMASS))
 320 C
 321 C          BOOST FINAL JETS
 322       DO 210 J=1,NJSET
 323       IF(J.EQ.IVIR.OR.J.EQ.IVIR2) GO TO 210
 324       IF(PJSET(5,J).LT.0.) GO TO 210
 325       IF(JDCAY(J).EQ.-1) GO TO 210
 326       BP1=0.
 327       BP21=0.
 328       DO 215 K=1,3
 329       BP1=BP1+BOOST1(K)*PJSET(K,J)
 330 215   BP21=BP21+B2B1(K)*PJSET(K,J)
 331       DO 220 K=1,3
 332 220   PJSET(K,J)=PJSET(K,J)
 333      $+(B2B1(K)*BI41+BOOST2(K)*BI42)*PJSET(4,J)
 334      $+B2B1(K)*BP1*BIK1+BOOST2(K)*BP21*BIK2+BOOST2(K)*BP1*BIK3
 335       PJSET(4,J)=B44*PJSET(4,J)+BP21*B4K1+BP1*B4K2
 336 210   CONTINUE
 337 C
 338 C          SET PFINAL TO BOOST2
 339       DO 230 K=1,4
 340 230   PFINAL(K)=BOOST2(K)
 341       PFINAL(5)=BMASS
 342 C
 343 C          RESET REMAINING VECTORS
 344       DO 240 J=NJSET,1,-1
 345       IF(J.EQ.IVIR.OR.J.EQ.IVIR2) GO TO 240
 346       IF(PJSET(5,J).GE.0.) GO TO 240
 347       JX1=JDCAY(J)/JPACK
 348       JX2=JDCAY(J)-JPACK*JX1
 349       DO 250 K=1,4
 350       PJSET(K,J)=PJSET(K,JX1)-PJSET(K,JX2)
 351 250   DBL1(K)=PJSET(K,J)
 352 #if defined(CERNLIB_SINGLE)
 353       AMJ=SQRT(ABS(DBL1(4)**2-DBL1(1)**2-DBL1(2)**2-DBL1(3)**2))
 354 #endif
 355 #if defined(CERNLIB_DOUBLE)
 356       AMJ=DSQRT(ABS(DBL1(4)**2-DBL1(1)**2-DBL1(2)**2-DBL1(3)**2))
 357 #endif
 358       PJSET(5,J)=-AMJ
 359 240   CONTINUE
 360 C
 361 C          RESET PFINAL, ETC.
 362 #if defined(CERNLIB_SINGLE)
 363       DO 300 K=1,4
 364 300   PFINAL(K)=PFINAL(K)+PJSET(K,NEWF)
 365       S1=PFINAL(4)**2-PFINAL(1)**2-PFINAL(2)**2-PFINAL(3)**2
 366       IF(S1.LT.0.) GO TO 9999
 367       PFINAL(5)=SQRT(S1)
 368       PTOTPL=PJSET(4,NEWV)+PJSET(3,NEWV)+PJSET(4,IVIR2)+PJSET(3,IVIR2)
 369       PTOTMN=PJSET(4,NEWV)-PJSET(3,NEWV)+PJSET(4,IVIR2)-PJSET(3,IVIR2)
 370 #endif
 371 #if defined(CERNLIB_DOUBLE)
 372 C          NEED DOUBLE PRECISION ON 32-BIT MACHINES
 373       CALL DBLVEC(PFINAL,DBL1)
 374       CALL DBLVEC(PJSET(1,NEWF),DBL2)
 375       DO 300 K=1,4
 376       DBL1(K)=DBL1(K)+DBL2(K)
 377 300   PFINAL(K)=DBL1(K)
 378       S1=DBL1(4)**2-DBL1(1)**2-DBL1(2)**2-DBL1(3)**2
 379       PFINAL(5)=SQRT(S1)
 380       IF(S1.LT.0.) GO TO 9999
 381       PFINAL(5)=SQRT(S1)
 382       PTOTPL=PJSET(4,NEWV)+PJSET(3,NEWV)+PJSET(4,IVIR2)+PJSET(3,IVIR2)
 383       PTOTMN=PJSET(4,NEWV)-PJSET(3,NEWV)+PJSET(4,IVIR2)-PJSET(3,IVIR2)
 384 #endif
 385 C
 386 C          SET NJSET AND POINTERS IF Z WAS GOOD
 387       IF(.NOT.ZGOOD) GO TO 10
 388       NJSET=NJSET+2
 389       IF(JDCAY(NEWF).GT.0) NJSET=NJSET+2
 390       JVIR(JET-10)=NEWV
 391       GO TO 10
 392 C          ERROR -- DISCARD EVENT.
 393 9999  CONTINUE
 394       WRITE(ITLIS,9998) IEVT
 395 9998  FORMAT(/' ***** ERROR IN QCDINI ... EVENT',I8,' DISCARDED *****')
 396       NJSET=-1
 397       RETURN
 398       END