--- /dev/null
+#include "isajet/pilot.h"
+ SUBROUTINE DECAY(IP)
+C
+C Decay particle IP from /PARTCL/ using /DKYTAB/ branching
+C ratios and add decay products to /PARTCL/ with IORIG=IP.
+C Forced decay modes are flagged by LOOK<0.
+C
+C Auxiliary routines:
+C DECPS1: generate masses for phase space
+C DECPS2: generate 2-body decays and boosts for phase space
+C DECVA: V-A matrix elements
+C DECTAU: tau decay matrix elements with polarization
+C DECSS3: 3-body SUSY matrix element using /DKYSS3/
+C DECJET: Hadronize partons from decay.
+C
+C Matrix element for Dalitz decays and W mass for TP -> W BT
+C are generated explicitly. W width is included.
+C
+C Requirements for decay modes:
+C (1) For Dalitz decays, particle 1 must be GM.
+C (2) For V-A quark or lepton decays, particles 1 and 2 must
+C be from (virtual) W.
+C (3) For any decay into quarks, they must appear last.
+C
+#if defined(CERNLIB_IMPNONE)
+ IMPLICIT NONE
+#endif
+#include "isajet/itapes.inc"
+#include "isajet/wcon.inc"
+#include "isajet/partcl.inc"
+#include "isajet/dkytab.inc"
+#include "isajet/jetset.inc"
+#include "isajet/jwork.inc"
+#include "isajet/const.inc"
+#include "isajet/primar.inc"
+#include "isajet/idrun.inc"
+#include "isajet/force.inc"
+#include "isajet/sstype.inc"
+#include "isajet/dkyss3.inc"
+C
+ REAL PGEN(5,5),BETA(3),REDUCE(5),WPROP,Z,TRY,RANF,AMASS,TWOME
+ REAL PSUM(5),SUM,PREST(4,6),DOT,PCM
+ REAL AMEE,REE,WTEE,SWAP,WT,A,B,C,GAMMA
+ REAL SMAX,SMIN,SVAL,TANMAX,TANMIN,TANVAL
+ LOGICAL WDECAY,DECVA,DECTAU,DECJET
+ INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IPOINT,ID1,I1,I2
+ INTEGER NADD,NSTART,NEW,NADD1,J,IP,I,IDABS(5)
+ INTEGER K,JETIP,IDANTI,NPASS,MEIP,MEA
+ REAL DBLPCM,DECSS3,VAL
+C
+ DATA REDUCE/1.,1.,2.,5.,15./
+ DATA PSUM/5*0./
+ DATA TWOME/1.022006E-3/
+ DATA PREST/24*0./
+C
+C Function definitions.
+C Use double precision for PCM on 32-bit machines
+C
+#if defined(CERNLIB_SINGLE)
+ PCM(A,B,C)=SQRT((A**2-B**2-C**2)**2-(2.*B*C)**2)/(2.*A)
+#endif
+#if defined(CERNLIB_DOUBLE)
+ PCM(A,B,C)=DBLPCM(A,B,C)
+#endif
+ DOT(I1,I2)=PREST(4,I1)*PREST(4,I2)-PREST(1,I1)*PREST(1,I2)
+ $-PREST(2,I1)*PREST(2,I2)-PREST(3,I1)*PREST(3,I2)
+C Charged W propagator.
+ WPROP(Z)=(Z-WMASS(2)**2)**2+(WMASS(2)*WGAM(2))**2
+C----------------------------------------------------------------------
+C Select decay mode. Note IDENT(NPTCL+1)...IDENT(NPTCL+5)
+C are always defined even if zero.
+C----------------------------------------------------------------------
+ IF(IDCAY(IP).NE.0) RETURN
+ IDLV1=IDENT(IP)
+ CALL FLAVOR(IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX)
+C FLAVOR returns 0 for quark, but want IFL3=6 for top
+ IF(IABS(IDLV1).LT.10) IFL3=IDLV1
+ NPASS=0
+1 CONTINUE
+ NPASS=NPASS+1
+ WDECAY=.FALSE.
+ IF(NPASS.GT.NTRIES) GO TO 9998
+ IPOINT=LOOK(INDEX)
+ IF(IPOINT.EQ.0) RETURN
+C IPOINT<0 flags a forced decay.
+ IF(IPOINT.LT.0) THEN
+ I=1
+ IF(IDENT(IP).LT.0) I=2
+ IPOINT=LOOK2(I,IABS(IPOINT))
+ ENDIF
+C
+C Select mode.
+C
+ TRY=RANF()
+ IPOINT=IPOINT-1
+100 IPOINT=IPOINT+1
+ IF(TRY.GT.CBR(IPOINT)) GO TO 100
+ NADD=0
+ SUM=0.
+ NSTART=NPTCL+1
+ IF(NPTCL+5.GT.MXPTCL) GO TO 9999
+C
+C Set up masses and IDENT codes.
+C
+ MEIP=MELEM(IPOINT)
+ DO 110 I=1,5
+ NEW=NPTCL+I
+ IDENT(NEW)=MODE(I,IPOINT)
+ IDABS(I)=IABS(IDENT(NEW))
+ IF(MODE(I,IPOINT).EQ.0) GO TO 110
+ NADD=NADD+1
+ IDLV1=IDENT(NEW)
+ PPTCL(5,NEW)=AMASS(IDLV1)
+ SUM=SUM+PPTCL(5,NEW)
+110 CONTINUE
+ NADD1=NADD-1
+ DO 120 J=1,5
+ PGEN(J,1)=PPTCL(J,IP)
+120 CONTINUE
+ PGEN(5,NADD)=PPTCL(5,NPTCL+NADD)
+C----------------------------------------------------------------------
+C Carry out appropriate decay
+C----------------------------------------------------------------------
+C
+C 1-body decays.
+C
+ IF(NADD.EQ.1) THEN
+ DO 200 J=1,5
+ PPTCL(J,NPTCL+1)=PPTCL(J,IP)
+200 CONTINUE
+ GO TO 300
+ ENDIF
+C
+C 2-body phase space decays
+C
+ IF(NADD.EQ.2.AND.MEIP.EQ.0) THEN
+ CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
+ GO TO 300
+ ENDIF
+C
+C N-body phase space decays
+C
+ IF(NADD.GT.2.AND.MEIP.EQ.0) THEN
+ CALL DECPS1(IP,NADD,PGEN)
+ CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
+ GO TO 300
+ ENDIF
+C
+C Dalitz decays
+C
+ IF(NADD.EQ.3.AND.MEIP.EQ.1) THEN
+210 AMEE=TWOME*(PPTCL(5,IP)/TWOME)**RANF()
+ REE=(TWOME/AMEE)**2
+ WTEE=(1.-(AMEE/PPTCL(5,IP))**2)**3*SQRT(1.-REE)*(1.+.5*REE)
+ IF(WTEE.LT.RANF()) GO TO 210
+ PGEN(5,2)=AMEE
+ CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
+ GO TO 300
+ ENDIF
+C
+C omega/phi decays (for reasons lost in history...)
+C
+ IF(NADD.EQ.3.AND.MEIP.EQ.2) THEN
+220 CALL DECPS1(IP,NADD,PGEN)
+ CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
+ WT=(PPTCL(5,NPTCL+1)*PPTCL(5,NPTCL+2)*PPTCL(5,NPTCL+3))**2
+ $ -(PPTCL(5,NPTCL+1)*DOT(2,3))**2
+ $ -(PPTCL(5,NPTCL+2)*DOT(1,3))**2
+ $ -(PPTCL(5,NPTCL+3)*DOT(1,2))**2
+ $ +2.*DOT(1,2)*DOT(2,3)*DOT(1,3)
+ IF(WT.LT.RANF()*PPTCL(5,IP)**6/108.) GO TO 220
+ GO TO 300
+ ENDIF
+C
+C V-A decays
+C
+ IF(NADD.EQ.3.AND.MEIP.EQ.3) THEN
+230 CALL DECPS1(IP,NADD,PGEN)
+ CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
+ IF(.NOT.DECVA(IP,NADD,IDABS,PREST)) GO TO 230
+ GO TO 300
+ ENDIF
+C
+C Top decays
+C Generate mass for TP -> W BT with Breit-Wigner.
+C W couples to 1+2 so swap 1<->3. Then m2+m3 < m < m0-m1.
+C
+ IF(NADD.EQ.3.AND.MEIP.EQ.4) THEN
+ WDECAY=.TRUE.
+ SWAP=PPTCL(5,NPTCL+1)
+ PPTCL(5,NPTCL+1)=PPTCL(5,NPTCL+3)
+ PPTCL(5,NPTCL+3)=SWAP
+ SMAX=(PPTCL(5,IP)-PPTCL(5,NPTCL+1))**2
+ SMIN=(PPTCL(5,NPTCL+2)+PPTCL(5,NPTCL+3))**2
+ TANMAX=ATAN((SMAX-WMASS(2)**2)/(WMASS(2)*WGAM(2)))
+ TANMIN=ATAN((SMIN-WMASS(2)**2)/(WMASS(2)*WGAM(2)))
+240 TANVAL=RANF()*(TANMAX-TANMIN)+TANMIN
+ SVAL=WMASS(2)**2+WMASS(2)*WGAM(2)*TAN(TANVAL)
+ IF(SVAL.LT.SMIN.OR.SVAL.GT.SMAX) GO TO 240
+ PGEN(5,2)=SQRT(SVAL)
+ PGEN(5,3)=PPTCL(5,NPTCL+3)
+ CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
+ IF(.NOT.DECVA(IP,NADD,IDABS,PREST)) GO TO 240
+ DO 241 K=1,5
+ SWAP=PPTCL(K,NPTCL+1)
+ PPTCL(K,NPTCL+1)=PPTCL(K,NPTCL+3)
+ PPTCL(K,NPTCL+3)=SWAP
+241 CONTINUE
+ PGEN(5,3)=PPTCL(5,NPTCL+3)
+ DO 242 K=1,4
+ SWAP=PREST(K,1)
+ PREST(K,1)=PREST(K,3)
+ PREST(K,3)=SWAP
+242 CONTINUE
+ GO TO 300
+ ENDIF
+C
+C TAU decays. These are special because they take polarization
+C into account.
+C
+ IF(MEIP.EQ.5.OR.MEIP.EQ.6.OR.MEIP.EQ.7) THEN
+250 CALL DECPS1(IP,NADD,PGEN)
+ CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
+ IF(.NOT.DECTAU(IP,NADD,MEIP,IDABS,PREST)) GO TO 250
+ GO TO 300
+ ENDIF
+C
+C 3-body SUSY decays
+C
+ IF(MEIP.LT.0.AND.NADD.EQ.3) THEN
+ MEA=IABS(MEIP)
+ IF(WTSS3(MEA).LE.0) THEN
+ DO 260 I=1,1000
+ CALL DECPS1(IP,NADD,PGEN)
+ CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
+ VAL=DECSS3(IP,MEA)
+ WTSS3(MEA)=MAX(WTSS3(MEA),VAL)
+260 CONTINUE
+ IF(WTSS3(MEA).LE.0) GO TO 9998
+ ENDIF
+261 CALL DECPS1(IP,NADD,PGEN)
+ CALL DECPS2(IP,NADD,PGEN,PREST,BETA,GAMMA)
+ VAL=DECSS3(IP,MEA)
+ WTSS3(MEA)=MAX(WTSS3(MEA),VAL)
+ IF(VAL.LT.WTSS3(MEA)*RANF()) GO TO 261
+ GO TO 300
+ ENDIF
+C
+C Should never fall through
+C
+ GO TO 9998
+C----------------------------------------------------------------------
+C Swap particles and antiparticles if IDENT(IP)<0
+C Note forced modes for antiparticles are conjugated in table.
+C----------------------------------------------------------------------
+300 CONTINUE
+ IF(IDENT(IP).LT.0.AND.IDENT(IP).NE.-20) THEN
+ DO 310 I=1,NADD
+ ID1=IDENT(NPTCL+I)
+ IDENT(NPTCL+I)=IDANTI(ID1)
+310 CONTINUE
+ ENDIF
+C
+C Set IORIG and IDCAY.
+C
+ NPTCL=NPTCL+NADD
+ IDCAY(IP)=IPACK*NSTART+NPTCL
+ JETIP=IABS(IORIG(IP))/IPACK
+ DO 320 I=NSTART,NPTCL
+ IORIG(I)=IP
+ IDCAY(I)=0
+320 CONTINUE
+C
+C Evolve and hadronize partons. If it fails, start over.
+C
+ IF(IDABS(NADD).LT.10.OR.MOD(IDENT(NPTCL),100).EQ.0) THEN
+ IF(.NOT.DECJET(IP,NADD,IDABS,PREST,WDECAY,BETA,GAMMA))
+ $ GO TO 1
+ ENDIF
+C
+ RETURN
+C----------------------------------------------------------------------
+C Error messages.
+C----------------------------------------------------------------------
+9999 CALL PRTEVT(0)
+ WRITE(ITLIS,99990) NPTCL
+99990 FORMAT(//5X,'ERROR IN DECAY...NPTCL > ',I6)
+ RETURN
+9998 CALL PRTEVT(0)
+ WRITE(ITLIS,99980) IP
+99980 FORMAT(//5X,'ERROR IN DECAY...NO DECAY FOUND FOR PARTICLE',I6)
+ RETURN
+ END
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff