COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cmuBrsmp.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine cmubrsmpp (Emu, prob, path)
 
subroutine cmubrsmpe (Emu, Eg)
 

Function/Subroutine Documentation

◆ cmubrsmpe()

subroutine cmubrsmpe ( real*8  Emu,
real*8  Eg 
)

Definition at line 32 of file cmuBrsmp.f.

References d0, masmu, mulogf0, rndc(), and true.

Referenced by cmuinte().

32  implicit none
33 #include "Zcmuint.h"
34 #include "Zmass.h"
35 
36  real*8 emu ! input. muon total energy in GeV
37  real*8 eg ! output. sampled energy loss of muon
38 !
39 ! brems v=Eg/Emu distribution:
40 ! dv cmuBremLogf*(4(1 - v)/3 + v)/v
41 !
42  real*8 vc, vmx, cmuvmax2, term1, term2, u1, u2, x
43  real*8 delta, logf, func, u, v
44 
45  vc = mubrvmin
46  vmx = cmuvmax2(emu)
47  term1 =4.d0/3.d0 * (log(vmx/vc) - (vmx-vc))
48  term2 = (vmx-vc)*(vmx+vc)/2.d0
49 ! loop for rejection
50  100 continue
51  call rndc(u)
52  if(u .le. term1/(term1+term2)) then
53 ! (1/x -1)dx
54  do while (.true.)
55 ! average number of trials is 1.0xx; xx depends on vc
56  call rndc(u)
57  x = vc**u
58  call rndc(u)
59  if( u .lt. (1.0-x)) then
60  if(x .lt. vmx) goto 10
61  endif
62  enddo
63  10 continue
64  else
65 ! x dx in [vc:vmx]
66  do while (.true.)
67  call rndc(u1)
68  call rndc(u2)
69  x = max(u1,u2)
70  if(x .gt. vc .and. x .lt. vmx) goto 20
71  enddo
72  20 continue
73  endif
74  v = x
75 ! final rejection by cmuBremLogf.
76 ! we use specially made cmuBremLogf here
77 
78 ! mim. momentum transfer in unit of Mmu
79  delta = (masmu/emu) * v /(1.d0-v)/2
80  logf = muakm / (1.d0 + muakm2*delta)
81  if(zeff .gt. 10.) then
82  logf = logf * 2./3.d0 /zeff3
83  endif
84  func = log(logf) ! = cmuBremLogf
85  call rndc(u)
86  if(u .gt. func/mulogf0) goto 100
87  eg = v * emu
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
subroutine rndc(u)
Definition: rnd.f:91
masmu
Definition: Zmass.h:5
real *8 function cmuvmax2(E)
Definition: cmuvmax2.f:2
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
! muon interaction related variables ! real Zeff3 ! Zeff **real muPrEmax1 integer muPrTXT ! real muBrEmax1 integer muBrTXT ! real muNEmax1 integer muNTXT ! real muNpwdEdxt real real mulogf0 real muBrLEmin common muintc mulogf0
Definition: Zcmuint.h:26
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cmubrsmpp()

subroutine cmubrsmpp ( real*8  Emu,
real*8  prob,
real*8  path 
)

Definition at line 2 of file cmuBrsmp.f.

References kintp3(), and rndc().

Referenced by csampmuint().

2  implicit none
3 #include "Zcmuint.h"
4 
5  real*8 emu ! input. muon total energy in GeV
6  real*8 prob ! output. muon brems creation prob. /X0
7  real*8 path ! output. sampled path in r%l
8  ! if Emu < muBrEmin, prob=0
9  ! and path becomes big
10 
11  real*8 u, ale
12 
13  if(emu .le. mubremin) then
14  prob = 0.
15  elseif(emu .le. mubremax1) then
16  ale = log10(emu)
17  call kintp3(mubrtx, 1, mubrtxt, mubrlemin,
18  * mubrdetx, ale, prob)
19  else
20 ! prob is const for v > vmin; Emu > muBrEmax1
21  prob = mubrtx(mubrtxt)
22  endif
23  if(prob .gt. 0.) then
24  call rndc(u)
25  path =- log(u)/prob
26  else
27  path = 1.d30
28  endif
subroutine rndc(u)
Definition: rnd.f:91
subroutine kintp3(f, intv, n, x1, h, x, ans)
Definition: kintp3.f:19
Here is the call graph for this function:
Here is the caller graph for this function: