COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cmigdFunc.f
Go to the documentation of this file.
1 ! ***************************************
2 ! cmigdG: G(s)
3 ! cmigdPsi: Psi(s)
4 ! cmigdGzai: Gsai(s)
5 ! cmigdSforBrem: s for Bremsstrahlung
6 ! cmigdSforPair: s for pair creation
7 !
8 ! ***************************************
9 ! --------------------------------
10 ! to test cmigdSforBrem/cmigdSforPair
11 !
12 ! real*8 cmigdSforPair, s, v, e, rho, x
13 !
14 !
15 ! e = 1.e19/1.e9 ! GeV
16 !
17 ! x = 1.e19/1.e9 * 1.e-2 ! GeV kg/m3
18 ! do while (x .lt. 1.e22/1.e9 * 1.2)
19 ! v = 1.e-4
20 ! do while (v .lt. 1.)
21 ! rho = x / e !in assume kg/m3
22 ! s = cmigdSforPair(e, rho, v)
23 ! v = v*10.**0.01
24 ! enddo
25 ! write(*, *)
26 ! x = x * 10.**0.5
27 ! enddo
28 ! end
29 !
30 ! --------------------------------
31 ! to test cmigdPsi:
32 ! real*8 cmigdPsi, psimig, s
33 ! do s = 0., 1., 0.01
34 ! write(*,*) cmigdPsi(s), psimig(s, 1.d-5)
35 ! enddo
36 ! end
37 ! --------------------------------
38 ! to test cmigdG:
39 ! real*8 cmigdG, gmigdl, s
40 ! do s = 0., 1., 0.01
41 ! write(*,*) cmigdG(s), gmigdl(s, 1.d-5)
42 ! enddo
43 ! end
44 
45 
46 ! ******************************************
47 !
48 ! G(s): better than 1% accuracy.
49 !
50 ! ******************************************
51  real*8 function cmigdg(s)
52  implicit none
53  real*8 s ! input
54 !
55 ! G(s)/( 12pi*s^2/(1+12pi*s^2)) is approximated by polinomials
56 !
57  real*8 xxx1(5) ! pol coeff. for s < 0.12 G(s) error is 3 x 10^-5
58  real*8 xxx2(7) ! pol coeff. for s > 0.11 error is 1 x 10^-3
59 
60  real*8 pi, const, f
61  integer i
62  parameter(pi= 3.141592654d0, const = 12.d0 * pi)
63 
64  data ( xxx1(i), i= 1, 5)/
65  1 0.99995995 , -6.2614068 , 57.629644 ,
66  2 -242.27058 , 450.83728
67  * /
68 
69  data ( xxx2(i), i= 1, 7)/
70  1 0.78283668 , -1.3317861 , 12.312840 ,
71  2 -33.848229 , 44.733673 , -29.201775 ,
72  3 7.5561750
73  * /
74 
75 
76  if(s .lt. 0.12) then
77  f = 0.
78  do i = 5, 2, -1
79  f = (f + xxx1(i))* s
80  enddo
81  f = f + xxx1(1)
82  cmigdg = f * const*s**2/(1. + const*s**2)
83  elseif(s .lt. 1.1) then
84  f = 0.
85  do i = 7, 2, -1
86  f = (f + xxx2(i))* s
87  enddo
88  f = f + xxx2(1)
89  cmigdg = f * const*s**2/(1. + const*s**2)
90  else
91  cmigdg = 1. - 0.022/s**4
92  endif
93 
94  end
95 ! *********************************************
96 !
97 ! psi(s): since psi--> 6s for s -->0, we approximate
98 ! psi(s)/ ( 6s/(1+6s) ) by a polynomial.
99 !
100 !
101 ! *********************************************
102  real*8 function cmigdpsi(s)
103  implicit none
104  real*8 s ! input
105 
106  real*8 f
107  real*8 xxx1(4) ! coef. of pol. s<0.3 better than 0.1 % for Psi
108  real*8 xxx2(5) ! // s> 0.3 //
109  integer i
110 
111  data ( xxx1(i), i= 1, 4)/
112  1 1.00 , 2.6978063 , -9.4242869 ,
113  2 11.468298
114  * /
115 
116  data ( xxx2(i), i= 1, 5)/
117  1 1.2095058 , 0.57895055 , -1.6531094 ,
118  2 1.4846143 , -0.46392960
119  * /
120 
121  if(s .lt. 0.3) then
122  f = 0.
123  do i = 4, 2, -1
124  f = (f + xxx1(i)) * s
125  enddo
126  f = f + xxx1(1)
127  cmigdpsi = f * 6.d0 *s /(1.d0 + 6.d0*s)
128  elseif(s .lt. 1.2) then
129  f = 0.
130  do i = 5, 2, -1
131  f = (f + xxx2(i)) * s
132  enddo
133  f = f + xxx2(1)
134  cmigdpsi = f * 6.d0 *s /(1.d0 + 6.d0*s)
135  else
136  cmigdpsi = 1. - 0.012/s**4
137  endif
138  end
139 ! ***************************************
140 !
141 ! cmigdGzai; this is for air
142 !
143 ! **************************************
144 !
145  real*8 function cmigdgzai(s)
146  implicit none
147  real*8 s
148  common /cmigdc/ s1, logs1, vm, logvm, smpa1, smpa2, suma,
149  * a2prob, x0
150  real*8 vm, logvm, smpa1, smpa2, suma, a2prob, X0
151  real*8 s1 ! /1.1185e-4/ ! (Z^0.333/183)^2 ; Z=7.25
152  real*8 logs1 ! /-9.0983/ ! ln(s1)
153 
154 
155 !
156  if( s .gt. 1.) then
157  cmigdgzai = 1.
158  elseif(s .gt. s1) then
159  cmigdgzai = log(s)/logs1 + 1.
160  else
161  cmigdgzai =2.
162  endif
163  end
164 ! **************************************
165 !
166 ! cmigdSforBrem; get s for brems defined recursively
167 !
168 ! **************************************
169  real*8 function cmigdsforbrem(ee, rho, v)
170  implicit none
171 #include "Zmass.h"
172 ! #include "Zelemagp.h"
173 !
174  real*8 ee ! input. Electron energy in GeV.
175  real*8 rho ! input. air density in kg/m3
176  real*8 v ! input. Fractional energy of gamma, Eg/Ee
177 !
178  real*8 x, ss, ss2, cmigdGzai, temp
179  integer i
180 
181  common /cmigdc/ s1, logs1, vm, logvm, smpa1, smpa2, suma,
182  * a2prob, x0
183  real*8 vm, logvm, smpa1, smpa2, suma, a2prob, X0
184  real*8 s1 ! /1.1185e-4/ ! (Z^0.333/183)^2 ; Z=7.25
185  real*8 logs1 ! /-9.0983/ ! ln(s1)
186 
187  if(v .eq. 1.) then
188  cmigdsforbrem = 1.
189  else
190  x = ee * rho ! GeV kg/m3
191  ss = 1.
192  temp = masele* x0/x *100. * v/(1.-v) !100: X0/rho in cm
193 
194  do i = 1, 10
195 ! normally 3 to 4 iteration is ok
196  ss2 = 1.37e3 * sqrt(temp/cmigdgzai(ss))
197  if(abs(ss/ss2 -1.d0) .lt. 3.e-3) goto 200
198  ss = ss2
199  enddo
200  ss2 = 1.
201  200 continue
202  cmigdsforbrem = ss2
203 !
204 ! write(*, *) i, sngl(ss2), sngl(v), sngl(ee*rho)
205 !
206  endif
207  end
208 ! **************************************
209 !
210 ! cmigdSforPair; get s for pair defined recursively
211 !
212 ! **************************************
213  real*8 function cmigdsforpair(eg, rho, v)
214  implicit none
215 #include "Zmass.h"
216 ! #include "Zelemagp.h"
217 !
218  real*8 eg ! input. Gamma energy in GeV.
219  real*8 rho ! input. air density in kg/m3
220  real*8 v ! input. Fractional energy of e-/e+ Ee/Eg
221 !
222  real*8 x, ss, ss2, cmigdGzai, temp
223  integer i
224  common /cmigdc/ s1, logs1, vm, logvm, smpa1, smpa2, suma,
225  * a2prob, x0
226  real*8 vm, logvm, smpa1, smpa2, suma, a2prob, X0
227  real*8 s1 ! /1.1185e-4/ ! (Z^0.333/183)^2 ; Z=7.25
228  real*8 logs1 ! /-9.0983/ ! ln(s1)
229 
230  if(v .eq. 1.) then
231  cmigdsforpair = 1.
232  else
233  x = eg * rho ! GeV kg/m3
234  ss = 1.
235  temp = masele* x0/x *100./v/(1.-v) !100: X0/rho in cm
236  do i = 1, 10
237 ! normally 3 to 4 iteration is ok
238  ss2 = 1.37e3 * sqrt(temp/cmigdgzai(ss))
239  if(abs(ss/ss2 -1.d0) .lt. 3.e-3) goto 200
240  ss = ss2
241  enddo
242  ss2 = 1.
243  200 continue
244  cmigdsforpair = ss2
245 !
246 ! write(*, *) i, sngl(ss2), sngl(v), sngl(eg*rho)
247 !
248  endif
249  end
250 ! ******************************************
251 ! this must be called before using cbremErgLPM,
252 ! cpairErgLPM
253 
254  subroutine csetlpmcnst(s1in, logs1in, vmin, X0in)
255  implicit none
256  real*8 s1in ! Migdals s1 = (Z^0.333/183)^2
257  real*8 logs1in ! log(s1)
258  real*8 vmin ! Egmin/Ee
259  real*8 X0in ! radiation length in kg/m^2. if it is
260  ! g/cm2, multiply it by 0.1/.00001 = 10.0
261  common /cmigdc/ s1, logs1, vm, logvm, smpa1, smpa2, suma,
262  * a2prob, x0
263  real*8 vm, logvm, smpa1, smpa2, suma, a2prob, X0
264  real*8 s1 ! /1.1185e-4/ ! (Z^0.333/183)^2 ; Z=7.25
265  real*8 logs1 ! /-9.0983/ ! ln(s1)
266 
267  s1 = s1in
268  logs1 = logs1in
269  vm = vmin
270  logvm = log(vm)
271  smpa1 = 0.5
272  smpa2 = -4.0*logvm/3.0
273  suma = smpa1 + smpa2
274  a2prob = smpa2/suma
275  x0 = x0in
276  end
277 
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
cmigdsforpair
real *8 function cmigdsforpair(eg, rho, v)
Definition: cmigdFunc.f:214
csetlpmcnst
subroutine csetlpmcnst(s1in, logs1in, vmin, X0in)
Definition: cmigdFunc.f:255
e
dE dx *! Nuc Int sampling table e
Definition: cblkMuInt.h:130
cmigdsforbrem
real *8 function cmigdsforbrem(ee, rho, v)
Definition: cmigdFunc.f:170
cmigdg
real *8 function cmigdg(s)
Definition: cmigdFunc.f:52
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
cmigdgzai
real *8 function cmigdgzai(s)
Definition: cmigdFunc.f:146
e3
block data cblkElemag data *AnihiE e3
Definition: cblkElemag.h:7
masele
const double masele
Definition: Zmass.h:2
cmigdpsi
real *8 function cmigdpsi(s)
Definition: cmigdFunc.f:103