ISAJET/code/decjet.F

   1 #include "isajet/pilot.h"
   2       LOGICAL FUNCTION DECJET(IP,NADD,IDABS,PREST,WDECAY,BETA,GAMMA)
   3 C
   4 C          Auxiliary routine for DECAY. Evolve and hadronize partons.
   5 C          Check conservation laws. Return TRUE if OK, FALSE otherwise.
   6 C
   7 C          IP = particle to be decayed.
   8 C          NADD = number of products (NPTCL+1, ..., NPTCL+NADD).
   9 C          IDABS = absolute values of decay IDENT's.
  10 C          PREST = 4-momenta in rest frame.
  11 C          WDECAY = logical flag for real W decay.
  12 C          BETA,GAMMA = boost parameters.
  13 C
  14 #if defined(CERNLIB_IMPNONE)
  15       IMPLICIT NONE
  16 #endif
  17 #include "isajet/itapes.inc"
  18 #include "isajet/wcon.inc"
  19 #include "isajet/partcl.inc"
  20 #include "isajet/dkytab.inc"
  21 #include "isajet/jetset.inc"
  22 #include "isajet/jwork.inc"
  23 #include "isajet/const.inc"
  24 C
  25       REAL PGEN(5,5),RND(5),U(3),BETA(3),IDQK(3),ROT(3,3),PSAVE(3)
  26      1,REDUCE(5),WPROP,Z,TRY,RANF,AMASS,TWOME,CHARGE
  27       REAL PSUM(5),POLD(4),PNEW(4),SUM,WTMAX,SUM1,SUM2
  28       REAL PREST(4,6),PWREST(5),BETAW(3),DOT,PCM
  29       REAL AMEE,REE,WTEE,SWAP,RNEW,WT,QCM,PHI,S12,S12MAX,GAMMAW,BP
  30       REAL PJET,CTHQK,STHQK,CPHIQK,SPHIQK,SUMQ,A,B,C,GAMMA
  31       REAL CHARGW
  32       LOGICAL WDECAY
  33       INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IPOINT,ID1,I1,I2,I3
  34       INTEGER NADD,NSTART,NEW,NADD1,J,IP,I,IDABS(5),NEXT
  35       INTEGER JJ1,II,K1,K,NJSAVE,NJSAV1,NJSAV2,NJ1,NPRTN,KK,NHDRN1
  36       INTEGER IFAIL,JSAVE,JETIP,JET,NJADD,NPTLV1,IDANTI,NPJET(5)
  37       INTEGER NHDRN,NPJET3,NPTCLW,NPBEG(5)
  38 C
  39 C          Copy decay products into /JETSET/ and do QCD evolution.
  40 C
  41       IF(NJSET+NADD.GT.MXJSET) GO TO 9998
  42       NJSAVE=NJSET
  43       NSTART=NPTCL-NADD+1
  44       NPTCL=NSTART-1
  45       DO 100 I=1,NADD
  46         NJSET=NJSET+1
  47         DO 110 K=1,4
  48 110     PJSET(K,NJSET)=PREST(K,I)
  49         PJSET(5,NJSET)=PPTCL(5,NPTCL+I)
  50         JORIG(NJSET)=JPACK*I
  51         JTYPE(NJSET)=IDENT(NPTCL+I)
  52         JDCAY(NJSET)=0
  53         JMATCH(NJSET)=JPACK*(NJSAVE+1)+NJSAVE+NADD
  54 100   CONTINUE
  55 C
  56 C          For heavy quarks match 1+2 and 3+(1+2). Boost 1+2 to rest.
  57 C
  58       IF(WDECAY) THEN
  59         JMATCH(NJSAVE+1)=NJSAVE+2
  60         JMATCH(NJSAVE+2)=NJSAVE+1
  61         NJSET=NJSET+1
  62         DO 120 K=1,4
  63           PWREST(K)=PJSET(K,NJSAVE+1)+PJSET(K,NJSAVE+2)
  64           PJSET(K,NJSET)=PWREST(K)
  65 120     CONTINUE
  66         PWREST(5)=SQRT(PWREST(4)**2-PWREST(1)**2-PWREST(2)**2
  67      $  -PWREST(3)**2)
  68         PJSET(5,NJSET)=PWREST(5)
  69         JMATCH(NJSAVE+3)=NJSAVE+4
  70         JMATCH(NJSAVE+4)=NJSAVE+3
  71         JORIG(NJSAVE+4)=-1
  72         IDLV1=JTYPE(NJSAVE+1)
  73         CHARGW=CHARGE(IDLV1)
  74         IDLV1=JTYPE(NJSAVE+2)
  75         CHARGW=CHARGW+CHARGE(IDLV1)
  76         JTYPE(NJSAVE+4)=80*SIGN(1.,CHARGW)
  77         JDCAY(NJSAVE+4)=0
  78 C          Boost W vectors to rest.
  79         DO 130 K=1,3
  80 130     BETAW(K)=PWREST(K)/PWREST(4)
  81         GAMMAW=PWREST(4)/PWREST(5)
  82         NJSAV1=NJSAVE+1
  83         NJSAV2=NJSAVE+2
  84         DO 140 J=NJSAV1,NJSAV2
  85           BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J)+BETAW(3)*PJSET(3,J)
  86           DO 141 K=1,3
  87 141       PJSET(K,J)=PJSET(K,J)-GAMMAW*BETAW(K)*(PJSET(4,J)
  88      $    -BP*GAMMAW/(GAMMAW+1.))
  89           PJSET(4,J)=GAMMAW*(PJSET(4,J)-BP)
  90 140     CONTINUE
  91       ENDIF
  92 C
  93 C          Do evolution and save new W momentum. Start from parent
  94 C          mass or NADD*energy.
  95       NJSAV1=NJSAVE+1
  96       DO 150 J=NJSAV1,NJSET
  97         IF(IABS(JTYPE(J)).LT.10.OR.MOD(JTYPE(J),100).EQ.0) THEN
  98           JDCAY(J)=-1
  99           PJSET(5,J)=AMIN1(PPTCL(5,IP),NADD*PJSET(4,J))
 100         ENDIF
 101 150   CONTINUE
 102 C
 103       CALL QCDJET(NJSAVE+1)
 104 C
 105       IF(WDECAY) THEN
 106         PWREST(4)=PJSET(4,NJSAVE+4)
 107         GAMMAW=PWREST(4)/PWREST(5)
 108         DO 200 K=1,3
 109           PWREST(K)=PJSET(K,NJSAVE+4)
 110           BETAW(K)=PWREST(K)/PWREST(4)
 111 200     CONTINUE
 112       ENDIF
 113 C
 114 C          Put final partons in particle table - temporary IORIG.
 115 C          Also include virtual or real W momentum for quark decays.
 116 C
 117       NJ1=NJSAVE+1
 118       IF(WDECAY) THEN
 119 C          Real or virtual W.
 120         NPTCL=NPTCL+1
 121         NPTCLW=NPTCL
 122         DO 210 K=1,5
 123 210     PPTCL(K,NPTCL)=PJSET(K,NJSAVE+4)
 124         IORIG(NPTCL)=IP
 125         IDENT(NPTCL)=JTYPE(NJSAVE+4)
 126         IDCAY(NPTCL)=0
 127 C          Jet 3
 128         NPBEG(3)=NPTCL+1
 129         DO 220 J=NJ1,NJSET
 130           IF(JDCAY(J).NE.0) GO TO 220
 131           IF(JORIG(J)/JPACK.NE.3) GO TO 220
 132           NPTCL=NPTCL+1
 133           DO 221 K=1,5
 134 221       PPTCL(K,NPTCL)=PJSET(K,J)
 135           IORIG(NPTCL)=3*IPACK+IP
 136           IDENT(NPTCL)=JTYPE(J)
 137           IDCAY(NPTCL)=0
 138 220     CONTINUE
 139 C          Jets 1 and 2
 140         NPJET3=NPTCL
 141         DO 230 JET=1,2
 142           NPBEG(JET)=NPTCL+1
 143           DO 240 J=NJ1,NJSET
 144             IF(JDCAY(J).NE.0) GO TO 240
 145             IF(JORIG(J)/JPACK.NE.JET) GO TO 240
 146             NPTCL=NPTCL+1
 147             BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J)
 148      $      +BETAW(3)*PJSET(3,J)
 149             DO 241 K=1,3
 150 241         PPTCL(K,NPTCL)=PJSET(K,J)+GAMMAW*BETAW(K)*(PJSET(4,J)
 151      $      +BP*GAMMAW/(GAMMAW+1.))
 152             PPTCL(4,NPTCL)=GAMMAW*(PJSET(4,J)+BP)
 153             PPTCL(5,NPTCL)=PJSET(5,J)
 154             IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+NPTCLW
 155             IDENT(NPTCL)=JTYPE(J)
 156             IDCAY(NPTCL)=0
 157 240       CONTINUE
 158 230     CONTINUE
 159 C          Quark decays to W plus jet 3; then W decays.
 160         IDCAY(IP)=IPACK*NPTCLW+NPJET3
 161         IDCAY(NPTCLW)=IPACK*(NPJET3+1)+NPTCL
 162       ELSE
 163 C          Not quark decay, so just copy partons.
 164         DO 250 JET=1,NADD
 165           NPBEG(JET)=NPTCL+1
 166           DO 260 J=NJ1,NJSET
 167             IF(JDCAY(J).NE.0) GO TO 260
 168             IF(JORIG(J)/JPACK.NE.JET) GO TO 260
 169             NPTCL=NPTCL+1
 170             DO 261 K=1,5
 171 261         PPTCL(K,NPTCL)=PJSET(K,J)
 172             IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+IP
 173             IDENT(NPTCL)=JTYPE(J)
 174             IDCAY(NPTCL)=0
 175 260       CONTINUE
 176 250     CONTINUE
 177         IDCAY(IP)=NSTART*IPACK+NPTCL
 178       ENDIF
 179       NHDRN=NPTCL
 180 C
 181 C          Hadronize quarks and rotate to proper angles.
 182 C
 183       DO 300 JET=1,NADD
 184         NPRTN=NPBEG(JET)-1
 185         DO 310 I=NJ1,NJSET
 186           IF(JDCAY(I).NE.0) GO TO 310
 187           IF(JORIG(I)/JPACK.NE.JET) GO TO 310
 188           NPRTN=NPRTN+1
 189           IF(IABS(JTYPE(I)).GE.10.AND.MOD(JTYPE(I),100).NE.0)
 190      $    GO TO 330
 191 C
 192 C          Fragment parton:
 193           NEXT=NPTCL+1
 194           PJET=SQRT(PJSET(1,I)**2+PJSET(2,I)**2+PJSET(3,I)**2)
 195           CTHQK=PJSET(3,I)/PJET
 196           STHQK=1.-CTHQK**2
 197           IF(STHQK.LT.1) THEN
 198             STHQK=SQRT(STHQK)
 199             CPHIQK=PJSET(1,I)/(PJET*STHQK)
 200             SPHIQK=PJSET(2,I)/(PJET*STHQK)
 201           ELSE
 202             STHQK=0
 203             CPHIQK=1
 204             SPHIQK=0
 205           ENDIF
 206           CALL JETGEN(I)
 207           IF(NEXT.GT.NPTCL) GO TO 310
 208           ROT(1,1)=CPHIQK*CTHQK
 209           ROT(2,1)=SPHIQK*CTHQK
 210           ROT(3,1)=-STHQK
 211           ROT(1,2)=-SPHIQK
 212           ROT(2,2)=CPHIQK
 213           ROT(3,2)=0.
 214           ROT(1,3)=CPHIQK*STHQK
 215           ROT(2,3)=SPHIQK*STHQK
 216           ROT(3,3)=CTHQK
 217 C
 218           DO 320 II=NEXT,NPTCL
 219             DO 321 K=1,3
 220               PSAVE(K)=PPTCL(K,II)
 221               PPTCL(K,II)=0.
 222 321         CONTINUE
 223             DO 322 K=1,3
 224             DO 322 KK=1,3
 225 322         PPTCL(K,II)=PPTCL(K,II)+ROT(K,KK)*PSAVE(KK)
 226             IORIG(II)=IPACK*JET+NPRTN
 227             IDCAY(II)=0
 228 320       CONTINUE
 229           IDCAY(NPRTN)=NEXT*IPACK+NPTCL
 230           GO TO 310
 231 C
 232 C          or add lepton:
 233 330       NPTCL=NPTCL+1
 234           DO 331 K=1,5
 235 331       PPTCL(K,NPTCL)=PJSET(K,I)
 236           IORIG(NPTCL)=IPACK*JET+NPRTN
 237           IDENT(NPTCL)=JTYPE(I)
 238           IDCAY(NPTCL)=0
 239           IDCAY(NPRTN)=NPTCL*IPACK+NPTCL
 240 310     CONTINUE
 241         NPJET(JET)=NPTCL
 242 300   CONTINUE
 243 C
 244 C          Reset NJSET so decay jets do not appear in /JETSET/
 245       NJADD=NJSET
 246       NJSET=NJSAVE
 247 C
 248 C          Check for at least two particles
 249       IF(NPTCL.LT.NHDRN+2) THEN
 250         NPTCL=NSTART-1
 251         DECJET=.FALSE.
 252         RETURN
 253       ENDIF
 254 C
 255 C          Conserve charge
 256 C
 257       SUMQ=0.
 258       NHDRN1=NHDRN+1
 259       DO 400 I=NHDRN1,NPTCL
 260         IDLV1=IDENT(I)
 261         SUMQ=SUMQ+CHARGE(IDLV1)
 262 400   CONTINUE
 263       IDLV1=IDENT(IP)
 264       SUMQ=SUMQ-CHARGE(IDLV1)
 265 C
 266       IF(SUMQ.EQ.0.) GO TO 500
 267 C
 268 C          Charge wrong--fix it by swapping UP and DN quarks.
 269       DO 410 I=NHDRN1,NPTCL
 270         ID1=IDENT(I)
 271         IF(IABS(ID1).GT.1000) GO TO 410
 272         I1=MOD(IABS(ID1)/100,10)
 273         I2=MOD(IABS(ID1)/10,10)
 274         I3=MOD(IABS(ID1),10)
 275         IF(I1.EQ.1.AND.I2.GT.2.AND.SUMQ*ID1.GT.0.) THEN
 276           IDENT(I)=ISIGN(200+10*I2+I3,ID1)
 277         ELSEIF(I1.EQ.2.AND.I2.GT.2.AND.SUMQ*ID1.LT.0.) THEN
 278           IDENT(I)=ISIGN(100+10*I2+I3,ID1)
 279         ELSEIF(I1.EQ.1.AND.I2.EQ.2.AND.SUMQ*ID1.GT.0.) THEN
 280           IDENT(I)=110+I3
 281         ELSEIF(I1.EQ.1.AND.I2.EQ.1) THEN
 282           IDENT(I)=(120+I3)*(-SIGN(1.,SUMQ))
 283         ELSE
 284           GO TO 410
 285         ENDIF
 286         SUMQ=SIGN(ABS(SUMQ)-1.,SUMQ)
 287         IDLV1=IDENT(I)
 288         PPTCL(5,I)=AMASS(IDLV1)
 289         PPTCL(4,I)=SQRT(PPTCL(1,I)**2+PPTCL(2,I)**2+PPTCL(3,I)**2
 290      $  +PPTCL(5,I)**2)
 291 C          Sum cannot vanish for fractionally charged initial particle.
 292         IF(ABS(SUMQ).LT.1.) GO TO 500
 293 410   CONTINUE
 294 C          Failed to conserve charge.
 295       NPTCL=NSTART-1
 296       DECJET=.FALSE.
 297       RETURN
 298 C
 299 C          Rescale momenta for correct mass
 300 C
 301 500   CONTINUE
 302       IF(WDECAY) THEN
 303 C          Quark decay. First rescale jet3 + W
 304         DO 510 K=1,5
 305 510     PPTCL(K,NPTCL+1)=PPTCL(K,NPTCLW)
 306         NPTLV1=NPTCL+1
 307         DO 520 K=1,3
 308 520     PSUM(K)=0.
 309         PSUM(4)=PPTCL(5,IP)
 310         PSUM(5)=PSUM(4)
 311         CALL RESCAL(NPJET(2)+1,NPTLV1,PSUM,IFAIL)
 312         IF(IFAIL.NE.0) THEN
 313           NPTCL=NSTART-1
 314           DECJET=.FALSE.
 315           RETURN
 316         ENDIF
 317         DO 530 K=1,3
 318 530     BETAW(K)=PPTCL(K,NPTCL+1)/PPTCL(4,NPTCL+1)
 319         GAMMAW=PPTCL(4,NPTCL+1)/PPTCL(5,NPTCL+1)
 320 C          Then rescale W
 321         PSUM(4)=PPTCL(5,NPTCLW)
 322         PSUM(5)=PSUM(4)
 323         CALL RESCAL(NHDRN1,NPJET(2),PSUM,IFAIL)
 324         IF(IFAIL.NE.0) THEN
 325           NPTCL=NSTART-1
 326           DECJET=.FALSE.
 327           RETURN
 328         ENDIF
 329       ELSE
 330 C          General decay with no W.
 331         DO 550 K=1,3
 332 550     PSUM(K)=0.
 333         PSUM(4)=PPTCL(5,IP)
 334         PSUM(5)=PSUM(4)
 335         NPTLV1=NPTCL
 336         CALL RESCAL(NHDRN1,NPTLV1,PSUM,IFAIL)
 337         IF(IFAIL.NE.0) THEN
 338           NPTCL=NSTART-1
 339           DECJET=.FALSE.
 340           RETURN
 341         ENDIF
 342       ENDIF
 343 C
 344 C          Boost back to lab frame. Reset IORIG.
 345 C
 346       IF(WDECAY) THEN
 347         DO 600 I=NHDRN1,NPTCL
 348           JET=IORIG(I)/IPACK
 349           IF(JET.NE.1.AND.JET.NE.2) GO TO 600
 350           BP=BETAW(1)*PPTCL(1,I)+BETAW(2)*PPTCL(2,I)+BETAW(3)*PPTCL(3,I)
 351           DO 610 J=1,3
 352 610       PPTCL(J,I)=PPTCL(J,I)+GAMMAW*BETAW(J)*(PPTCL(4,I)
 353      $    +BP*GAMMAW/(GAMMAW+1.))
 354           PPTCL(4,I)=GAMMAW*(PPTCL(4,I)+BP)
 355 600     CONTINUE
 356       ENDIF
 357 C
 358       DO 620 I=NSTART,NPTCL
 359         IORIG(I)=MOD(IORIG(I),IPACK)
 360         BP=BETA(1)*PPTCL(1,I)+BETA(2)*PPTCL(2,I)+BETA(3)*PPTCL(3,I)
 361         DO 621 J=1,3
 362           PPTCL(J,I)=PPTCL(J,I)+GAMMA*BETA(J)*(PPTCL(4,I)
 363      $    +BP*GAMMA/(GAMMA+1.))
 364 621     CONTINUE
 365         PPTCL(4,I)=GAMMA*(PPTCL(4,I)+BP)
 366 620   CONTINUE
 367 C
 368 C          Normal exit
 369 C
 370       DECJET=.TRUE.
 371       RETURN
 372 C
 373 C          Error messages.
 374 C
 375 9998  DECJET=.FALSE.
 376       CALL PRTEVT(0)
 377       WRITE(ITLIS,99980) NJSET
 378 99980 FORMAT(//5X,'ERROR IN DECJET...NJSET > ',I5)
 379       RETURN
 380       END