ISAJET/code/decay.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE DECAY(IP)
   3 C
   4 C          Decay particle IP from /PARTCL/ using /DKYTAB/ branching
   5 C          ratios and add decay products to /PARTCL/ with IORIG=IP.
   6 C          Forced decay modes are flagged by LOOK<0.
   7 C
   8 C          Auxiliary routines:
   9 C          DECPS1: generate masses for phase space
  10 C          DECPS2: generate 2-body decays and boosts for phase space
  11 C          DECVA:  V-A matrix elements
  12 C          DECTAU: tau decay matrix elements with polarization
  13 C          DECSS3: 3-body SUSY matrix element using /DKYSS3/
  14 C          DECJET: Hadronize partons from decay.
  15 C
  16 C          Matrix element for Dalitz decays and W mass for TP -> W BT
  17 C          are generated explicitly. W width is included.
  18 C
  19 C          Requirements for decay modes:
  20 C          (1) For Dalitz decays, particle 1 must be GM.
  21 C          (2) For V-A quark or lepton decays, particles 1 and 2 must
  22 C              be from (virtual) W.
  23 C          (3) For any decay into quarks, they must appear last.
  24 C
  25 #if defined(CERNLIB_IMPNONE)
  26       IMPLICIT NONE
  27 #endif
  28 #include "isajet/itapes.inc"
  29 #include "isajet/wcon.inc"
  30 #include "isajet/partcl.inc"
  31 #include "isajet/dkytab.inc"
  32 #include "isajet/jetset.inc"
  33 #include "isajet/jwork.inc"
  34 #include "isajet/const.inc"
  35 #include "isajet/primar.inc"
  36 #include "isajet/idrun.inc"
  37 #include "isajet/force.inc"
  38 #include "isajet/sstype.inc"
  39 #include "isajet/dkyss3.inc"
  40 C
  41       REAL PGEN(5,5),BETA(3),REDUCE(5),WPROP,Z,TRY,RANF,AMASS,TWOME
  42       REAL PSUM(5),SUM,PREST(4,6),DOT,PCM
  43       REAL AMEE,REE,WTEE,SWAP,WT,A,B,C,GAMMA
  44       REAL SMAX,SMIN,SVAL,TANMAX,TANMIN,TANVAL
  45       LOGICAL WDECAY,DECVA,DECTAU,DECJET
  46       INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IPOINT,ID1,I1,I2
  47       INTEGER NADD,NSTART,NEW,NADD1,J,IP,I,IDABS(5)
  48       INTEGER K,JETIP,IDANTI,NPASS,MEIP,MEA
  49       REAL DBLPCM,DECSS3,VAL
  50 C
  51       DATA REDUCE/1.,1.,2.,5.,15./
  52       DATA PSUM/5*0./
  53       DATA TWOME/1.022006E-3/
  54       DATA PREST/24*0./
  55 C
  56 C          Function definitions.
  57 C          Use double precision for PCM on 32-bit machines
  58 C
  59 #if defined(CERNLIB_SINGLE)
  60       PCM(A,B,C)=SQRT((A**2-B**2-C**2)**2-(2.*B*C)**2)/(2.*A)
  61 #endif
  62 #if defined(CERNLIB_DOUBLE)
  63       PCM(A,B,C)=DBLPCM(A,B,C)
  64 #endif
  65       DOT(I1,I2)=PREST(4,I1)*PREST(4,I2)-PREST(1,I1)*PREST(1,I2)
  66      $-PREST(2,I1)*PREST(2,I2)-PREST(3,I1)*PREST(3,I2)
  67 C          Charged W propagator.
  68       WPROP(Z)=(Z-WMASS(2)**2)**2+(WMASS(2)*WGAM(2))**2
  69 C----------------------------------------------------------------------
  70 C          Select decay mode. Note IDENT(NPTCL+1)...IDENT(NPTCL+5)
  71 C          are always defined even if zero.
  72 C----------------------------------------------------------------------
  73       IF(IDCAY(IP).NE.0) RETURN
  74       IDLV1=IDENT(IP)
  75       CALL FLAVOR(IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX)
  76 C          FLAVOR returns 0 for quark, but want IFL3=6 for top
  77       IF(IABS(IDLV1).LT.10) IFL3=IDLV1
  78       NPASS=0
  79 1     CONTINUE
  80       NPASS=NPASS+1
  81       WDECAY=.FALSE.
  82       IF(NPASS.GT.NTRIES) GO TO 9998
  83       IPOINT=LOOK(INDEX)
  84       IF(IPOINT.EQ.0) RETURN
  85 C          IPOINT<0 flags a forced decay.
  86       IF(IPOINT.LT.0) THEN
  87         I=1
  88         IF(IDENT(IP).LT.0) I=2
  89         IPOINT=LOOK2(I,IABS(IPOINT))
  90       ENDIF
  91 C
  92 C          Select mode.
  93 C
  94       TRY=RANF()
  95       IPOINT=IPOINT-1
  96 100   IPOINT=IPOINT+1
  97       IF(TRY.GT.CBR(IPOINT)) GO TO 100
  98       NADD=0
  99       SUM=0.
 100       NSTART=NPTCL+1
 101       IF(NPTCL+5.GT.MXPTCL) GO TO 9999
 102 C
 103 C          Set up masses and IDENT codes.
 104 C
 105       MEIP=MELEM(IPOINT)
 106       DO 110 I=1,5
 107         NEW=NPTCL+I
 108         IDENT(NEW)=MODE(I,IPOINT)
 109         IDABS(I)=IABS(IDENT(NEW))
 110         IF(MODE(I,IPOINT).EQ.0) GO TO 110
 111         NADD=NADD+1
 112         IDLV1=IDENT(NEW)
 113         PPTCL(5,NEW)=AMASS(IDLV1)
 114         SUM=SUM+PPTCL(5,NEW)
 115 110   CONTINUE
 116       NADD1=NADD-1
 117       DO 120 J=1,5
 118         PGEN(J,1)=PPTCL(J,IP)
 119 120   CONTINUE
 120       PGEN(5,NADD)=PPTCL(5,NPTCL+NADD)
 121 C----------------------------------------------------------------------
 122 C          Carry out appropriate decay
 123 C----------------------------------------------------------------------
 124 C
 125 C          1-body decays.
 126 C
 127       IF(NADD.EQ.1) THEN
 128         DO 200 J=1,5
 129           PPTCL(J,NPTCL+1)=PPTCL(J,IP)
 130 200     CONTINUE
 131         GO TO 300
 132       ENDIF
 133 C
 134 C          2-body phase space decays
 135 C
 136       IF(NADD.EQ.2.AND.MEIP.EQ.0) THEN
 137         CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
 138         GO TO 300
 139       ENDIF
 140 C
 141 C          N-body phase space decays
 142 C
 143       IF(NADD.GT.2.AND.MEIP.EQ.0) THEN
 144         CALL DECPS1(IP,NADD,PGEN)
 145         CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
 146         GO TO 300
 147       ENDIF
 148 C
 149 C          Dalitz decays
 150 C
 151       IF(NADD.EQ.3.AND.MEIP.EQ.1) THEN
 152 210     AMEE=TWOME*(PPTCL(5,IP)/TWOME)**RANF()
 153         REE=(TWOME/AMEE)**2
 154         WTEE=(1.-(AMEE/PPTCL(5,IP))**2)**3*SQRT(1.-REE)*(1.+.5*REE)
 155         IF(WTEE.LT.RANF()) GO TO 210
 156         PGEN(5,2)=AMEE
 157         CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
 158         GO TO 300
 159       ENDIF
 160 C
 161 C          omega/phi decays (for reasons lost in history...)
 162 C
 163       IF(NADD.EQ.3.AND.MEIP.EQ.2) THEN
 164 220     CALL DECPS1(IP,NADD,PGEN)
 165         CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
 166         WT=(PPTCL(5,NPTCL+1)*PPTCL(5,NPTCL+2)*PPTCL(5,NPTCL+3))**2
 167      $  -(PPTCL(5,NPTCL+1)*DOT(2,3))**2
 168      $  -(PPTCL(5,NPTCL+2)*DOT(1,3))**2
 169      $  -(PPTCL(5,NPTCL+3)*DOT(1,2))**2
 170      $  +2.*DOT(1,2)*DOT(2,3)*DOT(1,3)
 171         IF(WT.LT.RANF()*PPTCL(5,IP)**6/108.) GO TO 220
 172         GO TO 300
 173       ENDIF
 174 C
 175 C          V-A decays
 176 C
 177       IF(NADD.EQ.3.AND.MEIP.EQ.3) THEN
 178 230     CALL DECPS1(IP,NADD,PGEN)
 179         CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
 180         IF(.NOT.DECVA(IP,NADD,IDABS,PREST)) GO TO 230
 181         GO TO 300
 182       ENDIF
 183 C
 184 C          Top decays
 185 C          Generate mass for TP -> W BT with Breit-Wigner.
 186 C          W couples to 1+2 so swap 1<->3. Then m2+m3 < m < m0-m1.
 187 C
 188       IF(NADD.EQ.3.AND.MEIP.EQ.4) THEN
 189         WDECAY=.TRUE.
 190         SWAP=PPTCL(5,NPTCL+1)
 191         PPTCL(5,NPTCL+1)=PPTCL(5,NPTCL+3)
 192         PPTCL(5,NPTCL+3)=SWAP
 193         SMAX=(PPTCL(5,IP)-PPTCL(5,NPTCL+1))**2
 194         SMIN=(PPTCL(5,NPTCL+2)+PPTCL(5,NPTCL+3))**2
 195         TANMAX=ATAN((SMAX-WMASS(2)**2)/(WMASS(2)*WGAM(2)))
 196         TANMIN=ATAN((SMIN-WMASS(2)**2)/(WMASS(2)*WGAM(2)))
 197 240     TANVAL=RANF()*(TANMAX-TANMIN)+TANMIN
 198         SVAL=WMASS(2)**2+WMASS(2)*WGAM(2)*TAN(TANVAL)
 199         IF(SVAL.LT.SMIN.OR.SVAL.GT.SMAX) GO TO 240
 200         PGEN(5,2)=SQRT(SVAL)
 201         PGEN(5,3)=PPTCL(5,NPTCL+3)
 202         CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
 203         IF(.NOT.DECVA(IP,NADD,IDABS,PREST)) GO TO 240
 204         DO 241 K=1,5
 205           SWAP=PPTCL(K,NPTCL+1)
 206           PPTCL(K,NPTCL+1)=PPTCL(K,NPTCL+3)
 207           PPTCL(K,NPTCL+3)=SWAP
 208 241     CONTINUE
 209         PGEN(5,3)=PPTCL(5,NPTCL+3)
 210         DO 242 K=1,4
 211           SWAP=PREST(K,1)
 212           PREST(K,1)=PREST(K,3)
 213           PREST(K,3)=SWAP
 214 242     CONTINUE
 215         GO TO 300
 216       ENDIF
 217 C
 218 C          TAU decays. These are special because they take polarization
 219 C          into account.
 220 C
 221       IF(MEIP.EQ.5.OR.MEIP.EQ.6.OR.MEIP.EQ.7) THEN
 222 250     CALL DECPS1(IP,NADD,PGEN)
 223         CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
 224         IF(.NOT.DECTAU(IP,NADD,MEIP,IDABS,PREST)) GO TO 250
 225         GO TO 300
 226       ENDIF
 227 C
 228 C          3-body SUSY decays
 229 C
 230       IF(MEIP.LT.0.AND.NADD.EQ.3) THEN
 231         MEA=IABS(MEIP)
 232         IF(WTSS3(MEA).LE.0) THEN
 233           DO 260 I=1,1000
 234             CALL DECPS1(IP,NADD,PGEN)
 235             CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
 236             VAL=DECSS3(IP,MEA)
 237             WTSS3(MEA)=MAX(WTSS3(MEA),VAL)
 238 260       CONTINUE
 239           IF(WTSS3(MEA).LE.0) GO TO 9998
 240         ENDIF
 241 261     CALL DECPS1(IP,NADD,PGEN)
 242         CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
 243         VAL=DECSS3(IP,MEA)
 244         WTSS3(MEA)=MAX(WTSS3(MEA),VAL)
 245         IF(VAL.LT.WTSS3(MEA)*RANF()) GO TO 261
 246         GO TO 300
 247       ENDIF
 248 C
 249 C          Should never fall through
 250 C
 251       GO TO 9998
 252 C----------------------------------------------------------------------
 253 C          Swap particles and antiparticles if IDENT(IP)<0
 254 C          Note forced modes for antiparticles are conjugated in table.
 255 C----------------------------------------------------------------------
 256 300   CONTINUE
 257       IF(IDENT(IP).LT.0.AND.IDENT(IP).NE.-20) THEN
 258         DO 310 I=1,NADD
 259           ID1=IDENT(NPTCL+I)
 260           IDENT(NPTCL+I)=IDANTI(ID1)
 261 310     CONTINUE
 262       ENDIF
 263 C
 264 C          Set IORIG and IDCAY.
 265 C
 266       NPTCL=NPTCL+NADD
 267       IDCAY(IP)=IPACK*NSTART+NPTCL
 268       JETIP=IABS(IORIG(IP))/IPACK
 269       DO 320 I=NSTART,NPTCL
 270         IORIG(I)=IP
 271         IDCAY(I)=0
 272 320   CONTINUE
 273 C
 274 C          Evolve and hadronize partons. If it fails, start over.
 275 C
 276       IF(IDABS(NADD).LT.10.OR.MOD(IDENT(NPTCL),100).EQ.0) THEN
 277         IF(.NOT.DECJET(IP,NADD,IDABS,PREST,WDECAY,BETA,GAMMA))
 278      $  GO TO 1
 279       ENDIF
 280 C
 281       RETURN
 282 C----------------------------------------------------------------------
 283 C          Error messages.
 284 C----------------------------------------------------------------------
 285 9999  CALL PRTEVT(0)
 286       WRITE(ITLIS,99990) NPTCL
 287 99990 FORMAT(//5X,'ERROR IN DECAY...NPTCL > ',I6)
 288       RETURN
 289 9998  CALL PRTEVT(0)
 290       WRITE(ITLIS,99980) IP
 291 99980 FORMAT(//5X,'ERROR IN DECAY...NO DECAY FOUND FOR PARTICLE',I6)
 292       RETURN
 293       END