]>
Commit | Line | Data |
---|---|---|
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 |