ISAJET/code/iswdky.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE ISWDKY
   3 C----------------------------------------------------------------------
   4 C-
   5 C-   Purpose and Methods :
   6 C-       decay W's and Z's as done in ISAJET
   7 C-
   8 C-   Created   6-MAY-1991   Serban D. Protopopescu
   9 C-
  10 C----------------------------------------------------------------------
  11 #if defined(CERNLIB_IMPNONE)
  12       IMPLICIT NONE
  13 #endif
  14 #include "isajet/const.inc"
  15 #include "isajet/frame.inc"
  16 #include "isajet/jetpar.inc"
  17 #include "isajet/jetset.inc"
  18 #include "isajet/jwork.inc"
  19 #include "isajet/pjets.inc"
  20 #include "isajet/partcl.inc"
  21 #include "isajet/primar.inc"
  22 #include "isajet/wcon.inc"
  23       REAL X(2)
  24       EQUIVALENCE (X(1),X1)
  25       REAL PREST(5),PL(5),EL(3),EML(3),EMSQL(3)
  26       REAL WTFAC(3)
  27       REAL BRANCH(29)
  28       INTEGER LISTJ(29),LISTW(4)
  29       REAL RANF,SUM,PTDEN,QDEN,ETA,
  30      $S12,SUMBR,BRMODE,AMASS,BRINV,TRY,PL12,
  31      $COSTHL,THL,PHL,PTL,SGN,BP,PLPL,PLMN,AMINI,AMFIN,PINI,PFIN,
  32      $ QPL,QMN,AM1SQ,AM2SQ,ROOT,P1PL,P1MN,P2PL,P2MN
  33       INTEGER NADD,K,IQ1,IQ2,IFL1,IFL2,IQ,IFL,I
  34       REAL EY
  35       REAL QWPL,QWMN
  36 C
  37       DATA LISTJ/
  38      $9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
  39      $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16,
  40      $10,80,-80,90/
  41       DATA LISTW/10,80,-80,90/
  42 C----------------------------------------------------------------------
  43 C
  44 C          Entry
  45 C
  46       NPTCL=0
  47 C
  48 C          Kinematics. Note that YW is the true rapidity and QW is
  49 C          the true 3-momentum. See DRLLYN.
  50 C
  51       QMW=QWJET(5)
  52       QTW=SQRT(QWJET(1)**2+QWJET(2)**2)
  53       QW=SQRT(QWJET(1)**2+QWJET(2)**2+QWJET(3)**2)
  54       IF(QTW.NE.0) THEN
  55         PHIW=ATAN2(QWJET(2),QWJET(1))
  56         IF(PHIW.LT.0) PHIW=PHIW+2*PI
  57       ELSE
  58         PHIW=0
  59       ENDIF
  60       QWPL=QWJET(4)+QWJET(3)
  61       QWMN=QWJET(4)-QWJET(3)
  62       IF(QWPL.GT.0..AND.QWMN.GT.0.) THEN
  63         YW=0.5*ALOG(QWPL/QWMN)
  64       ELSE
  65         YW=999.*SIGN(1.,QWJET(3))
  66       ENDIF
  67       IF(QW.NE.0.) THEN
  68         THW=ACOS(QWJET(3)/QW)
  69       ELSE
  70         THW=0.
  71       ENDIF
  72 C
  73 C          Select W decay mode
  74 C          QMW dependence neglected in branching ratios
  75 C          BRANCH is cum. br. with heavy modes subtracted.
  76 C
  77       S12=QMW**2
  78       BRANCH(1)=0.
  79       SUMBR=0.
  80       DO 105 IQ1=2,25
  81         IQ2=MATCH(IQ1,JWTYP)
  82         IF(IQ2.EQ.0) THEN
  83           BRMODE=0.
  84         ELSE
  85           BRMODE=WCBR(IQ1,JWTYP)-WCBR(IQ1-1,JWTYP)
  86           IFL1=LISTJ(IQ1)
  87           IFL2=LISTJ(IQ2)
  88           IF(S12.LE.(AMASS(IFL1)+AMASS(IFL2))**2) BRMODE=0.
  89         ENDIF
  90         BRANCH(IQ1)=BRANCH(IQ1-1)+BRMODE
  91         SUMBR=SUMBR+BRMODE
  92 105   CONTINUE
  93       BRINV=1./SUMBR
  94 C
  95       TRY=RANF()
  96       DO 110 IQ=1,25
  97         IF(TRY.LT.BRANCH(IQ)*BRINV.AND.MATCH(IQ,JWTYP).NE.0) THEN
  98           JETTYP(1)=IQ
  99           JETTYP(2)=MATCH(IQ,JWTYP)
 100           GO TO 120
 101         ENDIF
 102 110   CONTINUE
 103 C
 104 120   IFL1=LISTJ(JETTYP(1))
 105       IFL2=LISTJ(JETTYP(2))
 106 C
 107 C          Select masses of decay products.
 108 C
 109       EML(1)=AMASS(IFL1)
 110       EML(2)=AMASS(IFL2)
 111 C
 112 C          Generate W decay in its rest frame
 113 C          First set up momenta of decay products:
 114 C
 115       EMSQL(1)=EML(1)**2
 116       EMSQL(2)=EML(2)**2
 117       EL(1)=(S12+EMSQL(1)-EMSQL(2))/(2.*QMW)
 118       EL(2)=(S12+EMSQL(2)-EMSQL(1))/(2.*QMW)
 119       PL12=SQRT((S12-(EML(1)+EML(2))**2)*(S12-(EML(1)-EML(2))**2))
 120      $/(2.*QMW)
 121 C          W momentum
 122       DO 140 K=1,5
 123 140   PREST(K)=QWJET(K)
 124 C          Generate next W decay
 125 20    CONTINUE
 126       COSTHL=2.*RANF()-1.
 127       THL=ACOS(COSTHL)
 128       PHL=2.*PI*RANF()
 129       PTL=PL12*SIN(THL)
 130 C
 131       DO 300 I=1,2
 132         SGN=3-2*I
 133         PL(1)=SGN*PTL*COS(PHL)
 134         PL(2)=SGN*PTL*SIN(PHL)
 135         PL(3)=SGN*PL12*COSTHL
 136         PL(4)=EL(I)
 137         PL(5)=EML(I)
 138 C          Boost with W momentum
 139         BP=0.
 140         DO 310 K=1,3
 141 310     BP=BP+PL(K)*PREST(K)
 142         BP=BP/PREST(5)
 143         DO 320 K=1,3
 144 320     PL(K)=PL(K)+PREST(K)*PL(4)/PREST(5)
 145      $  +PREST(K)*BP/(PREST(4)+PREST(5))
 146         PL(4)=PL(4)*PREST(4)/PREST(5)+BP
 147 C          Fill common blocks
 148         PT(I)=SQRT(PL(1)**2+PL(2)**2)
 149         P(I)=SQRT(PT(I)**2+PL(3)**2)
 150         IF(PT(I).GT.0.) THEN
 151           PHI(I)=ATAN2(PL(2),PL(1))
 152         ELSE
 153           PHI(I)=(I-1)*PI
 154         ENDIF
 155         IF(PHI(I).LT.0.) PHI(I)=PHI(I)+2.*PI
 156         CTH(I)=PL(3)/P(I)
 157         STH(I)=PT(I)/P(I)
 158         TH(I)=ACOS(CTH(I))
 159         XJ(I)=PL(3)/HALFE
 160         IF(CTH(I).GT.0.) THEN
 161           PLPL=PL(4)+PL(3)
 162           PLMN=(PT(I)**2+EMSQL(I))/PLPL
 163         ELSE
 164           PLMN=PL(4)-PL(3)
 165           PLPL=(PT(I)**2+EMSQL(I))/PLMN
 166         ENDIF
 167         YJ(I)=.5*ALOG(PLPL/PLMN)
 168 300   CONTINUE
 169 C
 170 C          Set PJETS
 171 C
 172       DO 501 I=1,2
 173         PJETS(3,I)=P(I)*CTH(I)
 174         PJETS(1,I)=PT(I)*COS(PHI(I))
 175         PJETS(2,I)=PT(I)*SIN(PHI(I))
 176         PJETS(4,I)=SQRT(P(I)**2+EMSQL(I))
 177         PJETS(5,I)=SQRT(EMSQL(I))
 178         IDJETS(I)=LISTJ(JETTYP(I))
 179 501   CONTINUE
 180   999 RETURN
 181       END