COSMOS v7.655  COSMOSv7655
(AirShowerMC)
ckaonDecay.f
Go to the documentation of this file.
1 ! ******************************************************************
2 ! * *
3 ! * ckaonDecay
4 ! * *
5 ! ******************************************************************
6 !
7  subroutine ckaondecay(pj, mupol, a, np, polari)
8  implicit none
9 !---- include '../../Zptcl.h'
10 #include "Zptcl.h"
11 !---- include '../../Zcode.h'
12 #include "Zcode.h"
13 
14  integer np !output. no. of ptcls produced
15  type(ptcl):: pj ! input. kaon
16  logical mupol ! input. if T, muon polarization is considered
17  type(ptcl):: a(*) ! output. produced ptcls
18  real*8 polari ! output. polarization of the muon. if a containes
19 ! muon. muon is put in a(np), if any.
20 !
21 !
22  if(pj%charge .ne. 0) then
23 ! k+-
24  call ckchgdcy(pj, mupol, a, np, polari)
25  elseif(abs(pj%subcode) .eq. k0s) then
26 ! k0 short
27  call ckshortdecay(pj,mupol, a, np, polari)
28  else
29 ! k0 long
30  call cklongdecay(pj, mupol, a, np, polari)
31  endif
32  end
33  subroutine ckchgdcy(pj, mupol, a, np, polari)
34 !
35 ! k+- decay
36 ! -- process --
37 ! 1) k---->mu+neu (63.5 %)
38 ! 2) ---->pic + pi0 21
39 ! 3) ---->pi0+e+neu 4.8
40 ! 4) ---->pi0+mu+neu 3.2
41 ! 5) ---->pic+pic+pic 5.6
42 ! 6) ---->pic+pi0+pi0 1.7
43  implicit none
44 !
45 !---- include '../../Zptcl.h'
46 #include "Zptcl.h"
47 
48  integer np
49  type(ptcl):: pj
50  type(ptcl):: a(*)
51  logical mupol
52  real*8 polari
53 !
54  real*8 u
55  integer icp(2)/1, 0/, icm(2)/-1, 0/
56  integer icp3(3)/1, 1, -1/
57  integer icm3(3)/-1, 1, -1/
58  integer ic3p0(3)/1, 0, 0/, ic3m0(3)/-1,0,0/
59 !
60  call rndc(u)
61 
62  if(u .lt. 0.635) then
63 !//////////
64 ! write(0,'(a, 1p,e14.4)') 'kch 1 ', pj.fm.p(4)
65 !////////////
66 ! k-->mu + neu
67  call ckmudecay(pj, mupol, a, np, polari)
68  elseif(u .lt. .845) then
69 ! k+ --> pi+ + pi0 or c.c
70  if(pj%charge .gt. 0) then
71  call ck2pidecay(pj, icp, a, np)
72  else
73  call ck2pidecay(pj, icm, a, np)
74  endif
75  elseif(u .lt. .893) then
76 ! k+ ---> pi0 + e+ + neue
77 !//////////
78 ! write(0,'(a, 1p,e14.4)') 'kch 2 ', pj.fm.p(4)
79 !////////////
80  call ckpieneudecay(pj, a, np)
81  elseif(u .lt. .925) then
82 ! k+ --->mu+ + Neumu + pi0
83 !//////////
84 ! write(0,'(a, 1p,e14.4)') 'kch 3 ', pj.fm.p(4)
85 !////////////
86  if(mupol) then
87  call ckmuneupidcy(pj, a, np, polari)
88  else
89  call ckmuneupidcy2(pj, a, np)
90  polari = 0.
91  endif
92  elseif(u .lt. .981) then
93 ! k+ ---> 3*pic
94  if(pj%charge .eq. 1) then
95  call ck3pidecay(pj, icp3, a, np)
96  else
97  call ck3pidecay(pj, icm3, a, np)
98  endif
99  else
100  if(pj%charge .eq. 1.) then
101  call ck3pidecay(pj, ic3p0, a, np)
102  else
103  call ck3pidecay(pj, ic3m0, a, np)
104  endif
105  endif
106  end
107 ! *******************************************
108 ! k0s decay
109 ! -- process -- 2014.sep
110 ! 1) ---->pi+ + pi- 68.61%-->69.20
111 ! 2) ---->pi0 + pi0 31.39 -->30.69
112 ! 3) ---->pic + e + neu_e 7.04e-4
113 ! 4) ---->pic + mu + neu_mu 4.69e-4
114  subroutine ckshortdecay(pj, mupol, a, np, polari)
115  implicit none
116 !---- include '../../Zptcl.h'
117 #include "Zptcl.h"
118 #include "Zevhnp.h"
119  integer np
120  type(ptcl):: pj, a(*)
121  logical,intent(in):: mupol
122  real(8),intent(out):: polari
123 
124  integer ic1(2)/1, -1/, ic2(2)/0, 0/
125  real*8 u
126 
127  real(8),parameter::br(4)=(/69.20d0,30.69d0,7.04d-2,4.69d-2/)
128  logical,save::first=.true.
129  real(8),save:: cbr(4)
130  integer::i
131  if( first ) then
132  cbr(:)=br(:)
133  do i = 2, 4
134  cbr(i) = cbr(i-1)+ cbr(i)
135  enddo
136  if( k0ssemild == 1) then
137 ! 3) pi+e +nue is considered
138  cbr(:) = cbr(:)/cbr(3)
139  elseif(k0ssemild == 2) then
140 ! 4) pi+mu +numu is considered
141  cbr(3) = cbr(2)
142  cbr(:) = cbr(:)/cbr(4)
143  elseif( k0ssemild == 12) then
144 ! 3),4) are considerd
145  cbr(:) = cbr(:)/cbr(4)
146  elseif( k0ssemild == 0) then
147 ! 3),4) are neglected
148  cbr(:) = cbr(:)/cbr(2)
149  else
150  write(0,*) ' K0sSemiLD =', k0ssemild, ' invalid'
151  stop
152  endif
153 !///////////////
154 ! write(0,*) ' cbr =', cbr(:)
155 !/////////////
156  first = .false.
157  endif
158 
159  call rndc(u)
160  do i = 1, 4
161  if(u .lt. cbr(i)) exit
162  enddo
163 !/////
164 ! write(0,'(a, i3, 1p, e14.4)') 'k0s ', i, pj.fm.p(4)
165 !///////////
166  polari = 0.
167  if(i == 1 ) then
168  call ck2pidecay(pj, ic1, a, np)
169  elseif( i == 2 ) then
170  call ck2pidecay(pj, ic2, a, np)
171  elseif( i == 3 ) then
172  call ckpieneudecay(pj, a, np) ! use K0l case
173  elseif( i == 4) then ! use k0l case
174  if(mupol) then
175  call ckmuneupidcy(pj, a, np, polari)
176  else
177  call ckmuneupidcy2(pj, a, np)
178  polari = 0.
179  endif
180  endif
181  end
182 ! k0l decay
183 ! -- process --
184 ! 1) ---->e pi neue 38.7
185 ! 2) ---->mu pi neum 27.1 %. (k0==>mu+, k0bar==>mu-)
186 ! 3) ---->3 pi0 21.5
187 ! 4) ---->pi+ pi- pi0 12.4
188 !
189  subroutine cklongdecay(pj, mupol, a, np, polari)
190  implicit none
191 !---- include '../../Zptcl.h'
192 
193 #include "Zptcl.h"
194  logical mupol
195  integer np
196  type(ptcl):: pj ! kaon
197  type(ptcl):: a(*) ! outputn
198  real*8 polari ! output
199 
200  real*8 u
201  integer ic(3)/0, 0, 0/, ic2(3)/1, -1, 0/
202 !
203  call rndc(u)
204 
205  if(u .lt. .387) then
206 ! e + neue + pi
207 !//////////
208 ! write(0,'(a, 1p, e14.4)') 'k0l 1 ', pj.fm.p(4)
209 !///////////
210  call ckpieneudecay(pj, a, np)
211  elseif(u .lt. .658) then
212 !//////////
213 ! write(0,'(a, 1p, e14.4)') 'k0l 2 ', pj.fm.p(4)
214 !//////////
215 ! mu + neumu + pi
216  if(mupol) then
217  call ckmuneupidcy(pj, a, np, polari)
218  else
219  call ckmuneupidcy2(pj, a, np)
220  polari = 0.
221  endif
222  elseif(u .lt. .873) then
223 ! 3 pi0
224  call ck3pidecay(pj, ic, a, np)
225  else
226 ! pi+ pi- pi0
227  call ck3pidecay(pj, ic2, a, np)
228  endif
229  end
230 ! ****************************************************
231  subroutine ckmuneupidcy(pj, a, np, polari)
232 ! ****************************************************
233 ! k-> mu + neum + pi (parent may be k charge,or k0)
234 !
235  implicit none
236 !---- include '../../Zptcl.h'
237 #include "Zptcl.h"
238 !---- include '../../Zcode.h'
239 #include "Zcode.h"
240 
241  integer np
242  type(ptcl):: pj ! kaon
243  type(ptcl):: a(*) ! outputn
244  real*8 polari ! output
245 
246  real*8 u, ecm, cosa, f, pcm
247  integer jpa, i
248  logical ok
249 ! make ptcl
250 ! k+==>pi0 mu+ nue(m) k- ==> pi0 mu- nue_b(m)
251 ! k0l==> pi- mu+ nue(m) or => pi+ mu- nue_b(m)
252 
253 
254  if(pj%charge .eq. -1.) then
255  call cmkptc(kneumu, antip, 0, a(1))
256  call cmkptc(kpion, 0, 0, a(2))
257  call cmkptc(kmuon, 0, -1, a(3))
258  jpa = -1
259  elseif(pj%charge .eq. 1) then
260  call cmkptc(kneumu, regptcl, 0, a(1))
261  call cmkptc(kpion, 0, 0, a(2))
262  call cmkptc(kmuon, 0, 1, a(3))
263  jpa = 1
264  else
265 ! k0l; anti or regptcl
266  call rndc(u)
267  if(u .lt. 0.5) then
268  call cmkptc(kneumu, antip, 0, a(1))
269  call cmkptc(kpion, 0, 1, a(2))
270  call cmkptc(kmuon, 0, -1, a(3))
271  jpa = 1
272  else
273  call cmkptc(kneumu, regptcl, 0, a(1))
274  call cmkptc(kpion, 0, -1, a(2))
275  call cmkptc(kmuon, 0, 1, a(3))
276  jpa = -1
277  endif
278  endif
279 
280  ok = .false.
281  do while (.not. ok)
282 ! sample energy of neum at rest of k
283  call csampneuekl3(f)
284 
285  ecm=f*pj%mass
286 ! angle
287  call rndc(u)
288  cosa=2*u-1.
289 ! set px,py,pz
290  call cpcos2pxyz(cosa, ecm, a(1)%fm)
291  a(1)%fm%p(4) = ecm
292 
293 ! muon ; should be put in the last place in a
294  np =3
295  call csampmuekl3(f)
296  ecm=max(f*pj%mass, a(np)%mass*1.0001d0)
297  pcm=sqrt(ecm**2- a(np)%mass**2)
298 ! angle
299  call rndc(u)
300  cosa=2*u-1.
301  call cpcos2pxyz(cosa, pcm, a(np)%fm)
302  a(np)%fm%p(4) = ecm
303  a(2)%fm%p(4) = pj%mass- ecm - a(1)%fm%p(4)
304  if( a(2)%fm%p(4) .gt. a(2)%mass) then
305 ! simply to satisfy the conservation
306  a(2)%fm%p(1) = -a(1)%fm%p(1) - a(np)%fm%p(1)
307  a(2)%fm%p(2) = -a(1)%fm%p(2) - a(np)%fm%p(2)
308  a(2)%fm%p(3) = -a(1)%fm%p(3) - a(np)%fm%p(3)
309  ok = .true.
310  endif
311  enddo
312 
313  do i = 1, np
314  call cibst1(i, pj, a(i), a(i))
315  enddo
316 ! set muon polarization
317  call cmupolatlabk(jpa, a(np), pj, polari)
318 
319 
320  end
321 ! ******************************
322  subroutine ckmuneupidcy2(pj, a, np)
323 ! ******************************
324 ! k-> mu + neummu + pi (parent may be k charge,or k0)
325 ! all is considered but not polarization
326  implicit none
327 !---- include '../../Zptcl.h'
328 #include "Zptcl.h"
329 !---- include '../../Zcode.h'
330 #include "Zcode.h"
331 
332  integer np
333  type(ptcl):: pj ! kaon
334  type(ptcl):: a(*) ! outputn
335 
336  real*8 w, u
337  integer i, icon
338 
339 ! make 3 ptcls
340 ! k+==>pi0 mu+ nue(m) k- ==> pi0 mu- nue_b(m)
341 ! k0l==> pi- mu+ nue(m) or => pi+ mu- nue_b(m)
342  if(pj%charge .eq. 1) then
343  call cmkptc(kneumu, regptcl, 0, a(1))
344  call cmkptc(kpion, 0, 0, a(2))
345  call cmkptc(kmuon, 0, 1, a(3))
346  elseif(pj%charge .eq. -1) then
347  call cmkptc(kneumu, antip, 0, a(1))
348  call cmkptc(kpion, 0, 0, a(2))
349  call cmkptc(kmuon, 0, -1, a(3))
350  else
351  call rndc(u)
352  if(u .lt. 0.5) then
353  call cmkptc(kneumu, antip, 0, a(1))
354  call cmkptc(kpion, 0, 1, a(2))
355  call cmkptc(kmuon, 0, -1, a(3))
356  else
357  call cmkptc(kneumu, regptcl, 0, a(1))
358  call cmkptc(kpion, 0, -1, a(2))
359  call cmkptc(kmuon, 0, 1, a(3))
360  endif
361  endif
362  np=3
363 ! 3 body pure phase space
364  call cnbdcy(3, pj%mass, a, 0, w, icon)
365 ! boost to lab
366  do i = 1, np
367  call cibst1(i, pj, a(i), a(i))
368  enddo
369  end
370  subroutine ckpieneudecay(pj, a, np )
371 ! e + neue + pi ;
372 !
373 ! k+==>pi0 e+ neue; k-==>pi0 e- neue_b
374 ! k0l==>pi+ e- +neue_b or ==>pi- e+ neue
375  implicit none
376 !---- include '../../Zptcl.h'
377 #include "Zptcl.h"
378 !---- include '../../Zcode.h'
379 #include "Zcode.h"
380 
381  integer np
382  type(ptcl):: pj, a(*)
383  type(ptcl)::piesys
384  real*8 u, f, ecm, cosa
385  integer i, echg
386  logical ok
387 
388 
389  ok = .false.
390  echg = pj%charge
391  if(pj%charge .eq. -1.) then
392  call cmkptc(kneue, antip, 0, a(1))
393  call cmkptc(kelec, 0, echg, a(2))
394  call cmkptc(kpion, 0, 0, a(3))
395  elseif(pj%charge .eq. 1) then
396  call cmkptc(kneue, regptcl, 0, a(1))
397  call cmkptc(kelec, 0 , echg, a(2))
398  call cmkptc(kpion, 0, 0, a(3))
399  else
400  call rndc(u)
401  if(u .lt. 0.5) then
402  call cmkptc(kneue, regptcl, 0, a(1))
403  call cmkptc(kelec, 0, 1, a(2))
404  call cmkptc(kpion, 0, -1, a(3))
405  else
406  call cmkptc(kneue, antip, 0, a(1))
407  call cmkptc(kelec, 0, -1, a(2))
408  call cmkptc(kpion, 0, 1, a(3))
409  endif
410  endif
411 
412  do while ( .not. ok)
413 ! sample energy of neue at rest of k
414  call csampneuekl3(f)
415  ecm=f*pj%mass
416 ! angle
417  call rndc(u)
418  cosa=2*u-1.
419 ! set px,py,pz
420  a(1)%fm%p(4) = ecm
421  call cpcos2pxyz(cosa, ecm, a(1)%fm)
422 
423 ! make pi + e system
424  piesys%fm%p(4) = pj%mass - ecm
425  if( piesys%fm%p(4) .gt. a(2)%mass+ a(3)%mass) then
426  piesys%fm%p(1) = -a(1)%fm%p(1)
427  piesys%fm%p(2) = -a(1)%fm%p(2)
428  piesys%fm%p(3) = -a(1)%fm%p(3)
429  piesys%mass = piesys%fm%p(4)**2 -(
430  * piesys%fm%p(1)**2+ piesys%fm%p(2)**2+piesys%fm%p(3)**2 )
431  if(piesys%mass .gt. 0.) then
432  piesys%mass = sqrt(piesys%mass)
433  call c2bdcy(piesys, a(2), a(3))
434  ok = .true.
435  endif
436  endif
437  enddo
438  np=3
439 ! boost to lab
440  do i = 1, np
441  call cibst1(i, pj, a(i), a(i))
442  enddo
443  end
444 ! ******************************************************************
445 ! * ckMuDecay: k -> mu + neumu
446 ! ******************************************************************
447 !
448 !
449 ! decay of k ---> mu + neu ( b.r=.635 ) is treated.
450 !
451  subroutine ckmudecay(pj, mupol, a, np, polari)
452  implicit none
453 !---- include '../../Zptcl.h'
454 #include "Zptcl.h"
455 !---- include '../../Zcode.h'
456 #include "Zcode.h"
457 
458  integer np
459  type(ptcl):: pj ! input. kaon
460  type(ptcl):: a(*) ! output. ptcls produced
461  logical mupol ! input. T==> muon polarization taken into acc.
462  real*8 polari ! output. muon polarizaton.
463 ! muon must be put a(np).
464  integer charge, subcode
465 !
466 !
467 
468 ! make muon neutrino : muon set last
469  subcode = -pj%charge
470  call cmkptc(kneumu, subcode, 0, a(1))
471 ! make muon
472  charge = pj%charge
473  call cmkptc(kmuon, 0, charge, a(2))
474 !
475 ! k-->mu + neu
476  call c2bdcy(pj, a(1), a(2))
477 ! set polarization of muon
478  if(mupol) then
479  call ckmupolari(pj, a(2), polari)
480  else
481  polari=0.
482  endif
483  np=2
484  end
485 ! **************************************
486  subroutine ck2pidecay(pj, ic, a, np )
487 ! **************************************
488 ! k--> 2 pi (k+-->pi+ pi0 or k- --> pi- pi0)
489 ! (k0-->pi+ pi- or k0 --> pi0 +pi0)
490  implicit none
491 !---- include '../../Zptcl.h'
492 #include "Zptcl.h"
493 !---- include '../../Zcode.h'
494 #include "Zcode.h"
495 
496  integer np, ic(2)
497  type(ptcl):: pj, a(*)
498  call cmkptc(kpion, 0, ic(1), a(1))
499  call cmkptc(kpion, 0, ic(2), a(2))
500 
501  call c2bdcy(pj, a(1), a(2))
502  np=2
503  end
504 ! **************************************
505  subroutine ck3pidecay(pj, ic, a, np )
506 ! **************************************
507 ! k--> 3 pi;
508 
509  implicit none
510 !---- include '../../Zptcl.h'
511 #include "Zptcl.h"
512 !---- include '../../Zcode.h'
513 #include "Zcode.h"
514 
515  integer np, ic(3)
516  type(ptcl):: pj, a(*)
517 
518  real*8 w
519  integer icon, i
520 !
521  call cmkptc(kpion, 0, ic(1), a(1))
522  call cmkptc(kpion, 0, ic(2), a(2))
523  call cmkptc(kpion, 0, ic(3), a(3))
524  call cnbdcy(3, pj%mass, a, 0, w, icon)
525  np = 3
526  do i=1, np
527  call cibst1(i, pj, a(i), a(i))
528  enddo
529  end
csampmuekl3
subroutine csampmuekl3(f)
Definition: cpimuPolari.f:137
ckmudecay
subroutine ckmudecay(pj, mupol, a, np, polari)
Definition: ckaonDecay.f:452
k0s
max ptcl codes in the kseethru ! subcode integer k0s
Definition: Zcode.h:2
ckmuneupidcy2
subroutine ckmuneupidcy2(pj, a, np)
Definition: ckaonDecay.f:323
cibst1
subroutine cibst1(init, p1, p2, po)
Definition: cibst1.f:29
kelec
max ptcl codes in the kelec
Definition: Zcode.h:2
ck2pidecay
subroutine ck2pidecay(pj, ic, a, np)
Definition: ckaonDecay.f:487
true
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
kneue
max ptcl codes in the kneue
Definition: Zcode.h:2
ckmupolari
subroutine ckmupolari(kaon, muon, polari)
Definition: cpimuPolari.f:58
cnbdcy
subroutine cnbdcy(n, ecm, p, jw, w, icon)
Definition: cnbdcy.f:48
rndc
subroutine rndc(u)
Definition: rnd.f:91
regptcl
max ptcl codes in the kseethru ! subcode integer regptcl
Definition: Zcode.h:2
ckpieneudecay
subroutine ckpieneudecay(pj, a, np)
Definition: ckaonDecay.f:371
ckmuneupidcy
subroutine ckmuneupidcy(pj, a, np, polari)
Definition: ckaonDecay.f:232
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
kneumu
max ptcl codes in the kneumu
Definition: Zcode.h:2
cklongdecay
subroutine cklongdecay(pj, mupol, a, np, polari)
Definition: ckaonDecay.f:190
cpcos2pxyz
subroutine cpcos2pxyz(cosa, p, pxyz)
Definition: cpCos2pxyz.f:2
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
csampneuekl3
subroutine csampneuekl3(f)
Definition: cpimuPolari.f:98
ckaondecay
subroutine ckaondecay(pj, mupol, a, np, polari)
Definition: ckaonDecay.f:8
c2bdcy
subroutine c2bdcy(p, p1, p2)
Definition: c2bdcy.f:44
false
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon ! knockon is considered Obsolete *PhotoProd false
Definition: cblkElemag.h:7
cmkptc
subroutine cmkptc(code, subcode, charge, p)
Definition: cmkptc.f:15
ckchgdcy
subroutine ckchgdcy(pj, mupol, a, np, polari)
Definition: ckaonDecay.f:34
ckshortdecay
subroutine ckshortdecay(pj, mupol, a, np, polari)
Definition: ckaonDecay.f:115
ptcl
Definition: Zptcl.h:75
antip
max ptcl codes in the kseethru ! subcode integer antip
Definition: Zcode.h:2
kpion
max ptcl codes in the kpion
Definition: Zcode.h:2
cmupolatlabk
subroutine cmupolatlabk(jpa, muon, kaon, p)
Definition: cpimuPolari.f:200
kmuon
max ptcl codes in the kmuon
Definition: Zcode.h:2
ck3pidecay
subroutine ck3pidecay(pj, ic, a, np)
Definition: ckaonDecay.f:506