source: trunk/MagicSoft/Simulation/Corsika/Mmcs/kdecay.f

Last change on this file was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 12.5 KB
Line 
1 SUBROUTINE KDECAY( IGO )
2
3C-----------------------------------------------------------------------
4C K(AON) DECAY
5C
6C KAON DECAYS WITH FULL KINEMATIC, ENERGY AND MOMENTA CONSERVED
7C ALL SECONDARY PARTICLES ARE WRITTEN TO STACK
8C THIS SUBROUTINE IS CALLED FROM NUCINT
9C ARGUMENT: (TO CHARACTERIZE THE DECAYING KAON)
10C IGO = 1 K+
11C = 2 K-
12C = 3 K0S
13C = 4 K0L
14C
15C DESIGN : D. HECK IK3 FZK KARLSRUHE
16C-----------------------------------------------------------------------
17
18 IMPLICIT NONE
19*KEEP,CONST.
20 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
21 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
22*KEEP,DECAY.
23 COMMON /DECAY/ GAM345,COS345,PHI345
24 DOUBLE PRECISION GAM345(3),COS345(3),PHI345(3)
25*KEEP,IRET.
26 COMMON /IRET/ IRET1,IRET2
27 INTEGER IRET1,IRET2
28*KEEP,KAONS.
29 COMMON /KAONS/ CKA
30 DOUBLE PRECISION CKA(80)
31*KEEP,PAM.
32 COMMON /PAM/ PAMA,SIGNUM
33 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
34*KEEP,PARPAR.
35 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
36 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
37 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
38 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
39 INTEGER ITYPE,LEVL
40*KEEP,PARPAE.
41 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
42 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
43 * (CURPAR(4), PHI ), (CURPAR(5), H ),
44 * (CURPAR(6), T ), (CURPAR(7), X ),
45 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
46 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
47 * (CURPAR(12),ECM )
48*KEEP,POLAR.
49 COMMON /POLAR/ POLART,POLARF
50 DOUBLE PRECISION POLART,POLARF
51*KEEP,RANDPA.
52 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
53 DOUBLE PRECISION FAC,U1,U2
54 REAL RD(3000)
55 INTEGER ISEED(103,10),NSEQ
56 LOGICAL KNOR
57*KEEP,RUNPAR.
58 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
59 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
60 * MONIOU,MDEBUG,NUCNUC,
61 * CETAPE,
62 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
63 * N1STTR,MDBASE,
64 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
65 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
66 * ,GHEISH,GHESIG
67 COMMON /RUNPAC/ DSN,HOST,USER
68 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
69 REAL STEPFC
70 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
71 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
72 * N1STTR,MDBASE
73 INTEGER CETAPE
74 CHARACTER*79 DSN
75 CHARACTER*20 HOST,USER
76
77 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
78 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
79 * ,GHEISH,GHESIG
80*KEND.
81
82 DOUBLE PRECISION BETA3,COSTCM,COSTH3,GAMMA3,PHI3,RA,WORK1,WORK2
83 INTEGER I,ICHARG,IGO,IPI,J,M3
84C-----------------------------------------------------------------------
85
86 IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
87 444 FORMAT(' KDECAY: CURPAR=',1P,9E10.3)
88
89C COPY COORDINATES INTO SECPAR
90 DO 20 J = 5,8
91 SECPAR(J) = CURPAR(J)
92 20 CONTINUE
93
94C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95C DECAY OF K(+,-) (6 MODES)
96
97 IF ( IGO .LE. 2 ) THEN
98 21 CALL RMMAR( RD,1,1 )
99 RA = RD(1)
100
101C DECAY K(+,-) ----> MU(+,-) + NEUTRINO
102 IF ( RA .LE. CKA(23) ) THEN
103C NEUTRINO IS DROPPED
104 WORK1 = CKA(28) * GAMMA
105 WORK2 = CKA(29) * BETA * WORK1
106 CALL RMMAR( RD,2,1 )
107 COSTCM = RD(1) * 2.D0 - 1.D0
108C MU(+,-)
109 GAMMA3 = WORK1 + COSTCM * WORK2
110 BETA3 = SQRT( 1.D0 - 1.D0 / GAMMA3**2 )
111 COSTH3 = MIN( 1.D0, (GAMMA * GAMMA3 - CKA(28))
112 * / (BETA * GAMMA * BETA3 * GAMMA3) )
113 PHI3 = RD(2) * PI2
114 CALL ADDANG( COSTHE,PHI, COSTH3,PHI3, SECPAR(3),SECPAR(4) )
115 IF ( SECPAR(3) .GT. C(29) ) THEN
116 SECPAR(1) = 4 + IGO
117 SECPAR(2) = GAMMA3
118C DIRECTION OF PION IN THE MUON CM SYSTEM (= DIRECTION OF POLARIZATION)
119C SEE: G. BARR ET AL., PHYS. REV. D39 (1989) 3532, EQ. 5
120C POLART IS COS OF ANGLE BETWEEN KAON AND LABORATORY IN THE MU CM
121C POLARF IS ANGLE PHI AROUND THE LAB DIRECTION IN THE MU CM
122C POLART, POLARF WITH RESPECT TO THE MU DIRECTION IN THE LAB SYSTEM
123 POLART = ( 2.D0*PAMA(11)*GAMMA*C(6) / (PAMA(5)*GAMMA3)
124 * - C(6) - 1.D0 ) / ( BETA3 * (1.D0-C(6)) )
125 POLARF = PHI3 - PI
126C PION DIRECTION IS DIRECTION OF POLARIZATION FOR K+, OPPOSITE FOR K-
127 IF ( ITYPE .EQ. 12 ) THEN
128 POLART = -POLART
129 POLARF = POLARF + PI
130 ENDIF
131C GET THE POLARIZATION DIRECTION IN THE MU CM RELATIVE TO THE CORSIKA
132C COORDINATE SYSTEM
133 CALL ADDANG( SECPAR(3),SECPAR(4), POLART,POLARF,
134 * POLART,POLARF )
135 SECPAR(11) = POLART
136 SECPAR(12) = POLARF
137 CALL TSTACK
138 SECPAR(11) = 0.D0
139 SECPAR(12) = 0.D0
140 ENDIF
141
142C DECAY K(+,-) ----> PI(+,-) + PI(0)
143 ELSEIF ( RA .LE. CKA(47) ) THEN
144 M3 = ITYPE - 3
145 CALL DECAY1( ITYPE, M3, 7 )
146
147C DECAY K(+,-) ----> PI(+,-) + PI(+,-) + PI(-,+)
148 ELSEIF ( RA. LE. CKA(48) ) THEN
149 CALL DECAY6( PAMA(11), PAMA(8),PAMA(8),PAMA(8),
150 * CKA(51),CKA(52),CKA(53), CKA(54), 1 )
151C PI(+,-) AND PI(+,-) AND THIRD (ODD) PI(-,+)
152 DO 230 I = 1,3
153 CALL ADDANG( COSTHE,PHI, COS345(I),PHI345(I),
154 * SECPAR(3),SECPAR(4) )
155 IF ( SECPAR(3) .GT. C(29) ) THEN
156 IF ( I .LE. 2 ) THEN
157 SECPAR(1) = 7 + IGO
158 ELSE
159 SECPAR(1) = 10 - IGO
160 ENDIF
161 SECPAR(2) = GAM345(I)
162 CALL TSTACK
163 ENDIF
164 230 CONTINUE
165
166C DECAY K(+,-) ----> PI(0) + E(+,-) + NEUTRINO
167 ELSEIF ( RA. LE. CKA(49) ) THEN
168 CALL DECAY6( PAMA(11), PAMA(7),PAMA(2),0.D0,
169 * CKA(65),CKA(66),0.D0, CKA(67), 4 )
170C PI(0) AND E(+,-) / NEUTRINO IS DROPPED
171 DO 250 I = 1,2
172 CALL ADDANG( COSTHE,PHI, COS345(I),PHI345(I),
173 * SECPAR(3),SECPAR(4) )
174 IF ( SECPAR(3) .GT. C(29) ) THEN
175 IF ( I .EQ. 1 ) THEN
176 SECPAR(1) = 7.D0
177 ELSE
178 SECPAR(1) = 1 + IGO
179 ENDIF
180 SECPAR(2) = GAM345(I)
181 CALL TSTACK
182 ENDIF
183 250 CONTINUE
184
185C DECAY K(+,-) ----> PI(0) + MU(+,-) + NEUTRINO
186 ELSEIF ( RA. LE. CKA(50) ) THEN
187 CALL DECAY6( PAMA(11), PAMA(7),PAMA(5),0.D0,
188 * CKA(68),CKA(69),0.D0, CKA(70), 3 )
189C PI(0) AND MU(+,-) / NEUTRINO IS DROPPED
190 DO 260 I = 1,2
191 CALL ADDANG( COSTHE,PHI, COS345(I),PHI345(I),
192 * SECPAR(3),SECPAR(4) )
193 IF ( SECPAR(3) .GT. C(29) ) THEN
194 SECPAR(2) = GAM345(I)
195 IF ( I .EQ. 1 ) THEN
196 SECPAR(1) = 7.D0
197 ELSE
198 SECPAR(1) = 4 + IGO
199 IF ( SECPAR(1) .EQ. 6.D0 ) THEN
200C INVERT POLARIZATION DIRECTION FOR MU(-)
201 POLART = -POLART
202 POLARF = POLARF + PI
203 ENDIF
204C GET THE POLARIZATION DIRECTION IN THE MU CM RELATIVE TO THE CORSIKA
205C COORDINATE SYSTEM
206 CALL ADDANG( SECPAR(3),SECPAR(4), POLART, POLARF,
207 * POLART,POLARF )
208 SECPAR(11) = POLART
209 SECPAR(12) = POLARF
210 ENDIF
211 CALL TSTACK
212 ENDIF
213 SECPAR(11) = 0.D0
214 SECPAR(12) = 0.D0
215 260 CONTINUE
216
217C DECAY K(+,-) ----> PI(0) + PI(0) + PI(+,-)
218 ELSE
219 CALL DECAY6( PAMA(11), PAMA(7),PAMA(7),PAMA(8),
220 * CKA(55),CKA(56),CKA(57), CKA(58), 1 )
221C PI(0)'S AND PI(+,-)
222 DO 270 I = 1,3
223 CALL ADDANG( COSTHE,PHI, COS345(I),PHI345(I),
224 * SECPAR(3),SECPAR(4) )
225 IF ( SECPAR(3) .GT. C(29) ) THEN
226 IF ( I .LE. 2 ) THEN
227 SECPAR(1) = 7.D0
228 ELSE
229 SECPAR(1) = 7 + IGO
230 ENDIF
231 SECPAR(2) = GAM345(I)
232 CALL TSTACK
233 ENDIF
234 270 CONTINUE
235
236 ENDIF
237
238C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239C DECAY OF K0S (2 MODES)
240 ELSEIF ( IGO .EQ. 3 ) THEN
241
242 CALL RMMAR( RD,1,1 )
243C DECAY K0S ----> PI(+) + PI(-)
244 IF ( RD(1) .LE. CKA(24) ) THEN
245 CALL DECAY1( ITYPE, 8, 9 )
246
247C DECAY K0S ----> PI(0) + PI(0)
248 ELSE
249 CALL DECAY1( ITYPE, 7, 7 )
250
251 ENDIF
252
253C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254C DECAY OF K0L (4 MODES)
255 ELSEIF ( IGO .EQ. 4 ) THEN
256 CALL RMMAR( RD,1,1 )
257 RA = RD(1)
258
259C DECAY K0L ----> PI(+,-) + E(-,+) + NEUTRINO
260 IF ( RA .LE. CKA(27) ) THEN
261 CALL DECAY6( PAMA(10), PAMA(8),PAMA(2),0.D0,
262 * CKA(71),CKA(72),0.D0, CKA(73), 4 )
263 CALL RMMAR( RD,1,1 )
264C CHARGE ASYMMETRY PREFERS FORMATION OF PI(-)
265 ICHARG = INT(1.5016 + RD(1))
266C PI(+,-) AND E(-,+) / NEUTRINO IS DROPPED
267 DO 420 I = 1,2
268 CALL ADDANG( COSTHE,PHI, COS345(I),PHI345(I),
269 * SECPAR(3),SECPAR(4) )
270 IF ( SECPAR(3) .GT. C(29) ) THEN
271 SECPAR(1) = 10 - 3*I - (2*I-3)*ICHARG
272 SECPAR(2) = GAM345(I)
273 CALL TSTACK
274 ENDIF
275 420 CONTINUE
276
277C DECAY K0L ----> PI(+,-) + MU(-,+) + NEUTRINO
278 ELSEIF ( RA .LE. CKA(26) ) THEN
279 CALL DECAY6( PAMA(10), PAMA(8),PAMA(5),0.D0,
280 * CKA(74),CKA(75),0.D0, CKA(76), 3 )
281 CALL RMMAR( RD,1,1 )
282C CHARGE ASYMMETRY PREFERS FORMATION OF PI(-)
283 ICHARG = INT(1.5016 + RD(1))
284C PI(+,-) AND MU(-,+) / NEUTRINO IS DROPPED
285 DO 430 I = 1,2
286 CALL ADDANG( COSTHE,PHI, COS345(I),PHI345(I),
287 * SECPAR(3),SECPAR(4) )
288 IF ( I .EQ. 1 ) THEN
289 SECPAR(1) = 7 + ICHARG
290 IPI = SECPAR(1)
291 ENDIF
292 IF ( SECPAR(3) .GT. C(29) ) THEN
293 SECPAR(2) = GAM345(I)
294 IF ( I .EQ. 2 ) THEN
295 SECPAR(1) = 7 - ICHARG
296 IF ( SECPAR(1) .EQ. 6.D0 ) THEN
297C INVERT POLARIZATION DIRECTION FOR MU(-)
298 POLART = -POLART
299 POLARF = POLARF + PI
300 ENDIF
301C GET THE POLARIZATION DIRECTION IN THE MU CM RELATIVE TO THE CORSIKA
302C COORDINATE SYSTEM
303 CALL ADDANG( SECPAR(3),SECPAR(4), POLART,POLARF,
304 * POLART,POLARF )
305 SECPAR(11) = POLART
306 SECPAR(12) = POLARF
307 ENDIF
308 CALL TSTACK
309 ENDIF
310 SECPAR(11) = 0.D0
311 SECPAR(12) = 0.D0
312 430 CONTINUE
313
314C DECAY K0L ----> PI(0) + PI(0) + PI(0)
315 ELSEIF ( RA .LE. CKA(25) ) THEN
316C SEE: S.V. SOMALWAR ET AL., PHYS.REV.LET. 68(1992)2580
317 CALL DECAY6( PAMA(10), PAMA(7),PAMA(7),PAMA(7),
318 * CKA(59),-.00033D0,CKA(59), CKA(60), 1 )
319C PI(0)'S
320 SECPAR(1) = 7.D0
321 DO 440 I = 1,3
322 CALL ADDANG( COSTHE,PHI, COS345(I),PHI345(I),
323 * SECPAR(3),SECPAR(4) )
324 IF ( SECPAR(3) .GT. C(29) ) THEN
325 SECPAR(2) = GAM345(I)
326 CALL TSTACK
327 ENDIF
328 440 CONTINUE
329
330C DECAY K0L ----> PI(+) + PI(-) + PI(0)
331 ELSE
332 CALL DECAY6( PAMA(10), PAMA(8),PAMA(8),PAMA(7),
333 * CKA(61),CKA(62),CKA(63), CKA(64), 1 )
334C PI(+) AND PI(-) AND PI(0)
335 DO 450 I = 1,3
336 CALL ADDANG( COSTHE,PHI, COS345(I),PHI345(I),
337 * SECPAR(3),SECPAR(4) )
338 IF ( SECPAR(3) .GT. C(29) ) THEN
339 IF ( I .LE. 2 ) THEN
340 SECPAR(1) = 7 + I
341 ELSE
342 SECPAR(1) = 7.D0
343 ENDIF
344 SECPAR(2) = GAM345(I)
345 CALL TSTACK
346 ENDIF
347 450 CONTINUE
348
349 ENDIF
350 ENDIF
351
352C KILL CURRENT PARTICLE
353 IRET1 = 1
354
355 RETURN
356 END
Note: See TracBrowser for help on using the repository browser.