ISAJET/code/dectau.F

   1 #include "isajet/pilot.h"
   2       LOGICAL FUNCTION DECTAU(IP,NADD,MEIP,IDABS,PREST)
   3 C
   4 C          Compute matrix elements for polarized tau decay.
   5 C          Polarization determined by tau parent.
   6 C          Auxiliary routine for DECAY.
   7 C
   8 #if defined(CERNLIB_IMPNONE)
   9       IMPLICIT NONE
  10 #endif
  11 #include "isajet/itapes.inc"
  12 #include "isajet/wcon.inc"
  13 #include "isajet/partcl.inc"
  14 #include "isajet/dkytab.inc"
  15 #include "isajet/const.inc"
  16 #include "isajet/pjets.inc"
  17 #include "isajet/keys.inc"
  18 #include "isajet/xmssm.inc"
  19 #include "isajet/sspols.inc"
  20 #include "isajet/primar.inc"
  21 C
  22       REAL PREST(4,6),WT,TAUHEL,S12,S12MAX,PIP,CTHNU,PSUM(4),AMV2,WT1
  23       REAL DOT,DOT3,RANF,Z
  24       INTEGER IP,NADD,IDABS(5),IPAR,IDPAR,JET,INU,I,K,I1,I2,IDSIB
  25       INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IDIP
  26       INTEGER MEIP,IPX,IP1,IP2
  27 C
  28       DOT(I1,I2)=PREST(4,I1)*PREST(4,I2)-PREST(1,I1)*PREST(1,I2)
  29      $-PREST(2,I1)*PREST(2,I2)-PREST(3,I1)*PREST(3,I2)
  30       DOT3(I1,I2)=PREST(1,I1)*PREST(1,I2)+PREST(2,I1)*PREST(2,I2)
  31      $+PREST(3,I1)*PREST(3,I2)
  32 C
  33       IDIP=IDENT(IP)
  34       DECTAU=.TRUE.
  35       IF(IABS(IDIP).NE.16) GO TO 999
  36 C
  37 C          Use PREST(K,6) for spin vector
  38 C
  39       PIP=SQRT(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(3,IP)**2)
  40       DO 100 K=1,3
  41         PREST(K,6)=PPTCL(K,IP)/PIP
  42 100   CONTINUE
  43       PREST(4,6)=0.
  44 C
  45 C          Take helicity TAUHEL=0 unless TAU parent is TP, W+-, H+-,
  46 C          or some SUSY particles.
  47 C          Take account of 1-particle decays!
  48 C
  49
  50       IPX=IP
  51       TAUHEL=0.
  52       IPAR=0
  53       IDPAR=0
  54 110   IF(IORIG(IPX).GT.0) THEN
  55         IPAR=MOD(IORIG(IPX),IPACK)
  56         IDPAR=IDENT(IPAR)
  57         IF(IDPAR.EQ.IDIP) THEN
  58           IP1=IDCAY(IPAR)/IPACK
  59           IP2=MOD(IDCAY(IPAR),IPACK)
  60           IF(IP1.EQ.IP2) THEN
  61             IPX=IPAR
  62             GO TO 110
  63           ENDIF
  64         ENDIF
  65         IDPAR=IABS(IDPAR)
  66         IDSIB=0
  67 C          W/top parent
  68         IF(IDPAR.GT.100.AND.MOD(IDPAR/10,10).GE.6) THEN
  69           TAUHEL=-1.
  70         ELSEIF(IDPAR.EQ.80) THEN
  71           TAUHEL=-1.
  72 C          Charged Higgs parent
  73         ELSEIF(IDPAR.EQ.86) THEN
  74           TAUHEL=+1.
  75 C          SUSY parent - polarization also depends on sibling IDSIB
  76         ELSEIF(GOMSSM.AND.IDPAR.GT.20.AND.IDPAR.LT.80) THEN
  77           I1=IDCAY(IPAR)/IPACK
  78           I2=MOD(IDCAY(IPAR),IPACK)
  79           DO 120 I=I1,I2
  80             IF(IABS(IDENT(I)).GT.20.AND.IABS(IDENT(I)).LT.80)
  81      $      IDSIB=IABS(IDENT(I))
  82 120       CONTINUE
  83           IF (IDPAR.EQ.35) THEN
  84             TAUHEL=-1.
  85           ELSEIF (IDPAR.EQ.36) THEN
  86             IF (IDSIB.EQ.30) TAUHEL=PTAU1(1)
  87             IF (IDSIB.EQ.40) TAUHEL=PTAU1(2)
  88             IF (IDSIB.EQ.50) TAUHEL=PTAU1(3)
  89             IF (IDSIB.EQ.60) TAUHEL=PTAU1(4)
  90           ELSEIF (IDPAR.EQ.56) THEN
  91             IF (IDSIB.EQ.30) TAUHEL=PTAU2(1)
  92             IF (IDSIB.EQ.40) TAUHEL=PTAU2(2)
  93             IF (IDSIB.EQ.50) TAUHEL=PTAU2(3)
  94             IF (IDSIB.EQ.60) TAUHEL=PTAU2(4)
  95           ELSEIF (IDPAR.EQ.39) THEN
  96             IF(IDSIB.EQ.35) TAUHEL=-1.
  97             IF(IDSIB.EQ.30) TAUHEL=PTAUWZ
  98           ELSEIF (IDPAR.EQ.49.AND.IDSIB.EQ.35) THEN
  99             TAUHEL=-1.
 100           ELSEIF (IDPAR.EQ.40) THEN
 101             IF(IDSIB.EQ.36) TAUHEL=PTAUZ2(1)
 102             IF(IDSIB.EQ.56) TAUHEL=PTAUZ2(2)
 103             IF(IDSIB.EQ.30) TAUHEL=PTAUZZ
 104           ELSEIF (IDPAR.EQ.50) THEN
 105             IF(IDSIB.EQ.36) TAUHEL=PTAUZ3(1)
 106             IF(IDSIB.EQ.56) TAUHEL=PTAUZ3(2)
 107           ELSEIF (IDPAR.EQ.60) THEN
 108             IF(IDSIB.EQ.36) TAUHEL=PTAUZ4(1)
 109             IF(IDSIB.EQ.56) TAUHEL=PTAUZ4(2)
 110           ENDIF
 111         END IF
 112       ELSE
 113         IF(KEYS(3)) THEN
 114           IF(IABS(IDENTW).EQ.80) TAUHEL=-1.
 115         ELSE
 116           JET=IABS(IORIG(IP))/IPACK
 117           IF(JET.GT.0.AND.JET.LE.NJET) THEN
 118             IF(IDJETS(JET).EQ.80) TAUHEL=-1.
 119           ENDIF
 120         ENDIF
 121       ENDIF
 122 C
 123 C          Leptonic decays. DECTAU is always called for TAU- decay
 124 C          products, so selection is independent of IDENT(IP).
 125 C
 126       IF(MEIP.EQ.5) THEN
 127         IF(IDENT(NPTCL+1).LT.0) THEN
 128           WT=PPTCL(5,IP)*(PREST(4,1)-TAUHEL*DOT(1,6))*DOT(2,3)
 129         ELSEIF(IDENT(NPTCL+2).LT.0) THEN
 130           WT=PPTCL(5,IP)*(PREST(4,2)-TAUHEL*DOT(2,6))*DOT(1,3)
 131         ELSE
 132           WT=PPTCL(5,IP)*(PREST(4,3)-TAUHEL*DOT(3,6))*DOT(1,2)
 133         ENDIF
 134         IF(WT.LT.RANF()*PPTCL(5,IP)**4/8.) THEN
 135           DECTAU=.FALSE.
 136         ELSE
 137           DECTAU=.TRUE.
 138         ENDIF
 139         RETURN
 140 C
 141 C          Decay to PI + NUT, K + NUT
 142 C
 143       ELSEIF(MEIP.EQ.6) THEN
 144         INU=1
 145         IF(IDABS(2).EQ.15) INU=2
 146         CTHNU=DOT3(INU,6)/SQRT(DOT3(INU,INU))
 147         WT=1.-TAUHEL*CTHNU
 148         IF(WT.LT.RANF()*2.) THEN
 149           DECTAU=.FALSE.
 150         ELSE
 151           DECTAU=.TRUE.
 152         ENDIF
 153         RETURN
 154 C
 155 C          Decay to RHO + NUT, A1 + NUT, K* + NUT
 156 C
 157       ELSEIF(MEIP.EQ.7) THEN
 158         DO 210 I=1,NADD
 159 210     IF(IDABS(I).EQ.15) INU=I
 160         DO 220 K=1,4
 161           PSUM(K)=0.
 162           DO 221 I=1,NADD
 163             IF(I.EQ.INU) GO TO 221
 164             PSUM(K)=PSUM(K)+PREST(K,I)
 165 221       CONTINUE
 166 220     CONTINUE
 167         AMV2=PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2
 168         WT1=2.*AMV2/(2.*AMV2+PPTCL(5,IP)**2)
 169         CTHNU=DOT3(INU,6)/SQRT(DOT3(INU,INU))
 170         WT=WT1*(1.+TAUHEL*CTHNU)+(1.-WT1)*(1-TAUHEL*CTHNU)
 171         IF(WT.LT.RANF()*2.) THEN
 172           DECTAU=.FALSE.
 173         ELSE
 174           DECTAU=.TRUE.
 175         ENDIF
 176         RETURN
 177 C
 178 C          Ignore matrix element for all other decays
 179 C
 180       ELSE
 181         DECTAU=.TRUE.
 182         RETURN
 183       ENDIF
 184 C
 185 C          Error
 186 C
 187 999   CALL PRTEVT(0)
 188       WRITE(ITLIS,99999) IP
 189 99999 FORMAT(//5X,'ERROR IN DECTAU FOR PARTICLE',I6)
 190       END