| 0795afa3 |
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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame