ISAJET/code/setdky.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE SETDKY(LPRINT)
   3 C
   4 C          Read in decay table from tape ITDKY and set up /DKYTAB/.
   5 C          Then append forced decay modes and set LOOK to negative
   6 C          number pointing to LOOK2, which points to table.
   7 C          Forced decays for antiparticles are stored in conjugated
   8 C          form so that DECAY can always conjugate them.
   9 C
  10 C          Logical flag LPRINT controls printing of table.
  11 C
  12 C          Ver 7.41: Check version of decay table. Also read matrix
  13 C          element flags and save in MELEM:
  14 C          MELEM=0: Phase space
  15 C          MELEM=1: Dalitz decay
  16 C          MELEM=2: omega/phi decay
  17 C          MELEM=3: V-A
  18 C          MELEM=4: V-A plus W propagator (for top)
  19 C          MELEM=5: tau -> ell nu nu
  20 C          MELEM=6: tau -> nu pi/K
  21 C          MELEM=7: tau -> nu rho/a1
  22 C
  23 #if defined(CERNLIB_IMPNONE)
  24       IMPLICIT NONE
  25 #endif
  26 #include "isajet/itapes.inc"
  27 #include "isajet/force.inc"
  28 #include "isajet/dkytab.inc"
  29 #include "isajet/nodcay.inc"
  30 #include "isajet/ssmode.inc"
  31 #include "isajet/sstype.inc"
  32 #include "isajet/xmssm.inc"
  33 #include "isajet/keys.inc"
  34 C
  35       INTEGER IMODE(6),LOOP,IOLD,I,IRES,ITYPE,K,J,IPOINT
  36       INTEGER IFL1,IFL2,IFL3,JSPIN,INDEX,ID1,IDANTI,KTYPE,IRES2
  37       REAL    BR
  38       CHARACTER*8 LABEL,LMODE(6),LRES
  39       CHARACTER*8 IBLANK,LREAD(10),IQUIT
  40       LOGICAL LPRINT
  41       INTEGER NOUT,NTHAD
  42       PARAMETER (NOUT=33)
  43       PARAMETER (NTHAD=12)
  44       INTEGER IDOUT(NOUT),ITHAD(NTHAD),IDUMMY(5),MEOUT
  45       REAL SUMBR,SUMBR2,SUMGAM
  46       CHARACTER*40,V,VOLD,VISAJE
  47 C
  48       DATA IDOUT/
  49      $IDTP,ISGL,ISUPL,ISDNL,ISSTL,ISCHL,ISBT1,ISTP1,ISUPR,ISDNR,
  50      $ISSTR,ISCHR,ISBT2,ISTP2,ISEL,ISMUL,ISTAU1,ISNEL,ISNML,ISNTL,
  51      $ISER,ISMUR,ISTAU2,ISZ1,ISZ2,ISZ3,ISZ4,ISW1,ISW2,
  52      $ISHL,ISHH,ISHA,ISHC/
  53       DATA IQUIT/'////'/,IBLANK/' '/
  54       DATA ITHAD/-160,-260,-360,
  55      $  1160,1260,2260,2160,1360,2360,3160,3260,3360/
  56 C
  57 C          Print header for table.
  58 C
  59       IF(LPRINT) WRITE(ITLIS,10)
  60 10    FORMAT('1',30('*')/' *',28X,'*'/
  61      1' *',5X,'ISAJET DECAY TABLE',5X,'*'/
  62      2' *',28X,'*'/' ',30('*')//
  63      33X,'PART',16X,'DECAY MODE',16X,'CUM BR',10X,'IDENT',18X,
  64      4'DECAY IDENT'/)
  65 C
  66 C          Initialize. LOOP is the decay mode counter.
  67 C
  68       LOOP=0
  69       IOLD=0
  70       DO 100 I=1,MXLOOK
  71         LOOK(I)=0
  72 100   CONTINUE
  73       DO 110 I=1,MXFORC
  74         LOOK2(1,I)=0
  75         LOOK2(2,I)=0
  76 110   CONTINUE
  77 C
  78 C          Read in table, checking for valid version.
  79 C
  80       IF(NODCAY.OR.ITDKY.EQ.0) RETURN
  81       REWIND ITDKY
  82 C
  83       VOLD=VISAJE()
  84       READ(ITDKY,*) V
  85       IF(V.NE.VOLD) THEN
  86         WRITE(ITLIS,2000) V,VOLD
  87 2000    FORMAT(//
  88      $  '  ***WARNING: DECAY TABLE DOES NOT MATCH ISAJET VERSION'/
  89      $  '  ***DECAY VERSION  : ',A40/
  90      $  '  ***PROGRAM VERSION: ',A40)
  91       ENDIF
  92 C
  93 200   LOOP=LOOP+1
  94       IF(LOOP.GT.MXDKY) GO TO 9999
  95 220   DO 210 I=1,5
  96         IMODE(I)=0
  97         LMODE(I)=IBLANK
  98 210   CONTINUE
  99       READ(ITDKY,*) IRES,ITYPE,BR,IMODE
 100 C
 101       IF(IRES.NE.0) THEN
 102         IF(NOPI0.AND.IRES.EQ.110) GO TO 220
 103         IF(NOETA.AND.IRES.EQ.220) GO TO 220
 104         IF(IRES.NE.IOLD) THEN
 105           CALL FLAVOR(IRES,IFL1,IFL2,IFL3,JSPIN,INDEX)
 106           LOOK(INDEX)=LOOP
 107         ENDIF
 108         IOLD=IRES
 109         CBR(LOOP)=BR
 110         MELEM(LOOP)=ITYPE
 111         DO 240 I=1,5
 112           MODE(I,LOOP)=IMODE(I)
 113           IF(IMODE(I).NE.0) LMODE(I)=LABEL(IMODE(I))
 114 240     CONTINUE
 115         LRES=LABEL(IRES)
 116         IF(LPRINT) WRITE(ITLIS,20) LRES,(LMODE(K),K=1,5),
 117      1  BR,IRES,(IMODE(K),K=1,5)
 118 20      FORMAT(3X,A5,4X,5(A5,2X),F8.5,10X,I5,4X,5(I7,2X))
 119         GO TO 200
 120       ENDIF
 121 C
 122 C          Add HIGGS FOR WHIGGS
 123 C
 124       IF(KEYS(10).AND..NOT.GOMSSM) THEN
 125           SUMGAM=0
 126           SUMBR=0
 127           DO 244 J=1,NSSMOD
 128             IF(ISSMOD(J).EQ.81.AND.GSSMOD(J).GT.0) THEN
 129               SUMGAM=SUMGAM+GSSMOD(J)
 130             ENDIF
 131 244       CONTINUE
 132           DO 245 J=1,NSSMOD
 133             IF(ISSMOD(J).EQ.81.AND.GSSMOD(J).GT.0) THEN
 134               BSSMOD(J)=GSSMOD(J)/SUMGAM
 135             ENDIF
 136 245       CONTINUE
 137           DO 246 J=1,NSSMOD
 138             IF(ISSMOD(J).EQ.81.AND.BSSMOD(J).GT.0) THEN
 139               SUMBR=SUMBR+BSSMOD(J)
 140             ENDIF
 141 246       CONTINUE
 142 C          If modes exist, add them
 143           IF(SUMBR.LE.0) GO TO 249
 144           IRES=81
 145           LRES=LABEL(IRES)
 146           CALL FLAVOR(IRES,IFL1,IFL2,IFL3,JSPIN,INDEX)
 147           LOOK(INDEX)=LOOP+1
 148           SUMBR2=0
 149           DO 247 J=1,NSSMOD
 150             IF(ISSMOD(J).EQ.81.AND.BSSMOD(J).GT.0) THEN
 151               LOOP=LOOP+1
 152               SUMBR2=SUMBR2+BSSMOD(J)
 153               BR=SUMBR2/SUMBR
 154               CBR(LOOP)=BR
 155               MELEM(LOOP)=0
 156               DO 248 K=1,5
 157                 MODE(K,LOOP)=JSSMOD(K,J)
 158                 LMODE(K)=LABEL(MODE(K,LOOP))
 159 248           CONTINUE
 160               IF(LPRINT) WRITE(ITLIS,20) LRES,(LMODE(K),K=1,5),
 161      $        BR,IRES,(MODE(K,LOOP),K=1,5)
 162             ENDIF
 163 247       CONTINUE
 164 249     CONTINUE
 165       END IF
 166 C
 167 C          Add MSSM decay modes if applicable, OR H_SM FOR WHIGGS
 168 C
 169       IF(GOMSSM) THEN
 170         DO 250 I=1,NOUT
 171 C          Check for modes
 172           SUMBR=0
 173           DO 251 J=1,NSSMOD
 174             IF(ISSMOD(J).EQ.IDOUT(I).AND.BSSMOD(J).GT.0) THEN
 175               SUMBR=SUMBR+BSSMOD(J)
 176             ENDIF
 177 251       CONTINUE
 178 C          If modes exist, add them
 179           IF(SUMBR.LE.0) GO TO 250
 180           IRES=IDOUT(I)
 181           LRES=LABEL(IRES)
 182           CALL FLAVOR(IRES,IFL1,IFL2,IFL3,JSPIN,INDEX)
 183           LOOK(INDEX)=LOOP+1
 184           SUMBR2=0
 185           DO 252 J=1,NSSMOD
 186             IF(ISSMOD(J).EQ.IDOUT(I).AND.BSSMOD(J).GT.0) THEN
 187               LOOP=LOOP+1
 188               SUMBR2=SUMBR2+BSSMOD(J)
 189               BR=SUMBR2/SUMBR
 190               CBR(LOOP)=BR
 191               MELEM(LOOP)=MSSMOD(J)
 192               DO 253 K=1,5
 193                 MODE(K,LOOP)=JSSMOD(K,J)
 194                 LMODE(K)=LABEL(MODE(K,LOOP))
 195 253           CONTINUE
 196               IF(LPRINT) WRITE(ITLIS,20) LRES,(LMODE(K),K=1,5),
 197      $        BR,IRES,(MODE(K,LOOP),K=1,5)
 198             ENDIF
 199 252       CONTINUE
 200 250     CONTINUE
 201 C
 202 C          Top hadron decays
 203 C
 204         DO 260 I=1,NTHAD
 205 C          Check for modes
 206           SUMBR=0
 207           DO 261 J=1,NSSMOD
 208             IF(ISSMOD(J).EQ.6.AND.BSSMOD(J).GT.0) THEN
 209               SUMBR=SUMBR+BSSMOD(J)
 210             ENDIF
 211 261       CONTINUE
 212 C          If modes exist, add them -- conjugate for antimesons
 213           IF(SUMBR.LE.0) GO TO 260
 214           IRES=IABS(ITHAD(I))
 215           LRES=LABEL(IRES)
 216           CALL FLAVOR(IRES,IFL1,IFL2,IFL3,JSPIN,INDEX)
 217           LOOK(INDEX)=LOOP+1
 218           SUMBR2=0
 219           DO 262 J=1,NSSMOD
 220             IF(ISSMOD(J).EQ.6.AND.BSSMOD(J).GT.0) THEN
 221               LOOP=LOOP+1
 222               SUMBR2=SUMBR2+BSSMOD(J)
 223               BR=SUMBR2/SUMBR
 224               CBR(LOOP)=BR
 225               IF(IABS(JSSMOD(1,J)).LT.20.AND.IABS(JSSMOD(2,J)).LT.20
 226      $        .AND.IABS(JSSMOD(3,J)).LT.20.AND.IABS(JSSMOD(4,J)).LT.20
 227      $        .AND.IABS(JSSMOD(5,J)).LT.20) THEN
 228                 MELEM(LOOP)=4
 229               ELSE
 230                 MELEM(LOOP)=0
 231               ENDIF
 232               DO 263 K=1,5
 233                 IF(ITHAD(I).GT.0) THEN
 234                   MODE(K,LOOP)=JSSMOD(K,J)
 235                 ELSE
 236                   MODE(K,LOOP)=IDANTI(JSSMOD(K,J))
 237                 ENDIF
 238                 LMODE(K)=LABEL(MODE(K,LOOP))
 239 263           CONTINUE
 240               IF(LPRINT) WRITE(ITLIS,20) LRES,(LMODE(K),K=1,5),
 241      $        BR,IRES,(MODE(K,LOOP),K=1,5)
 242             ENDIF
 243 262       CONTINUE
 244 260     CONTINUE
 245       ENDIF
 246 C
 247 C          Set forced decay modes.
 248 C          LOOK(INDEX) = -IRES, where LOOK2(K,IRES) points to entries in
 249 C          decay table for IDENT>0 and IDENT<0.
 250 C          LOOKST(IRES) = standard LOOK value.
 251 C
 252       IF(NFORCE.EQ.0) GO TO 400
 253 C          Append each forced decay to table
 254       IRES=0
 255       DO 310 I=1,NFORCE
 256         IF(IFORCE(I).EQ.0) GO TO 310
 257         LOOP=LOOP+1
 258         IF(LOOP.GT.MXDKY) GO TO 9999
 259         CALL FLAVOR(IFORCE(I),IFL1,IFL2,IFL3,JSPIN,INDEX)
 260         IF(IFORCE(I).GT.0) THEN
 261           KTYPE=1
 262         ELSE
 263           KTYPE=2
 264         ENDIF
 265 C
 266         IF(LOOK(INDEX).GE.0) THEN
 267           IRES=IRES+1
 268           IF(IRES.GT.MXFORC) GO TO 9998
 269           LOOKST(IRES)=LOOK(INDEX)
 270           LOOK2(KTYPE,IRES)=LOOP
 271           LOOK2(3-KTYPE,IRES)=LOOKST(IRES)
 272           LOOK(INDEX)=-IRES
 273         ELSE
 274           IRES2=-LOOK(INDEX)
 275           IF(IRES2.GT.MXFORC) GO TO 9998
 276           LOOK2(KTYPE,IRES2)=LOOP
 277         ENDIF
 278 C          Set forced decay mode - conjugate if necessary
 279         IF(KTYPE.EQ.1) THEN
 280           DO 320 K=1,5
 281 320       MODE(K,LOOP)=MFORCE(K,I)
 282         ELSE
 283           DO 330 K=1,5
 284 330       MODE(K,LOOP)=IDANTI(MFORCE(K,I))
 285         ENDIF
 286         CBR(LOOP)=1.
 287 C          Set matrix element flag
 288         CALL ORDER(IFORCE(I),MFORCE(1,I),IDUMMY,MEOUT)
 289         MELEM(LOOP)=MEOUT
 290         MEFORC(I)=MEOUT
 291 310   CONTINUE
 292 C
 293 400   RETURN
 294 C
 295 C          Errors
 296 C
 297 9999  WRITE(ITLIS,3001) LOOP
 298 3001  FORMAT(//' ***** ERROR IN SETDKY ... DECAY COUNTER LOOP = ',
 299      $I6,' *****')
 300       STOP 99
 301 9998  WRITE(ITLIS,3002) IRES
 302 3002  FORMAT(//' ***** ERROR IN SETDKY ... FORCE COUNTER IRES = ',
 303      $I6,' *****')
 304       STOP 99
 305       END