COSMOS v7.655  COSMOSv7655
(AirShowerMC)
getST.f
Go to the documentation of this file.
1 #include "BlockData/cblkGene.h"
2  implicit none
3 #include "Zglobalc.h"
4 #include "Zmanagerp.h"
5 #include "ZrigCut.h"
6 #include "Zptcl.h"
7 #include "Zobs.h"
8 #include "Zobsp.h"
9 #include "Zprimary.h"
10 #include "Zprimaryc.h"
11 #include "Zprimaryv.h"
12 #include "Zincidentp.h"
13 #ifdef NEXT486
14 #define IMAG_P dimag
15 #else
16 #define IMAG_P imag
17 #endif
18  include 'Zflux.h'
19 
20  integer i
21 
22  real*8 cosmin, cosmax, cylr, cylh
23  real*8 ans, sum1, sum2, sum3
24  logical deg
25  call cerrormsg('**********IMPORTANT**********',1)
26  call cerrormsg(
27  * 'You must write radius and height of a cylinder',1)
28  call cerrormsg('or 0 0 after namelist parameter',1)
29  call cerrormsg(
30  * '(1 space line may be needed before data)', 1)
31  call creadparam(5)
32  call cbeginrun
33  call cprintprim(errorout)
34  cosmax = imag_p(coszenith)
35  cosmin = real(coszenith)
36  azmmax = imag_p(azimuth) + xaxisfromsouth
37  azmmin = real(Azimuth) + XaxisFromSouth
38 
39  if(cutofffile .eq. ' ') then
40  rigc = 0.
41  deg = .false.
42  else
43  rigc = 100.
44  if(zenvalue .eq. 'deg') then
45  deg = .true.
46  else
47  deg = .false.
48  endif
49  endif
50  zen1 = cosmin
51  zen2 = cosmax
52 
53  call csetcosdeg(0, .false.)
54  write(*,*)
55  write(*,*)
56  * ' cos region=',zen1, ' to ', zen2
57  write(*,*)
58  * ' fai region=', real(Azimuth),' to ',
59  * imag_p(azimuth), ' deg'
60  write(*,*) ' X-axis From South=', xaxisfromsouth, ' deg'
61 
62  write(*,*)
63  if(rigc .eq. 0) then
64  write(*,*) ' No rigidity cut is assumed'
65  write(*,*)
66  " ' primary Int(dI/dE) sum(cumlative);',
67  * ' sphere'
68  else
69  write(*,*) ' Rigidity cut is applied'
70  write(*,*)
71  " ' primary Int(dI/dE*RigCut) sum1(cumlative);',
72  * ' sphere'
73  endif
74  sum1 = 0.
75  do i = 1, prim%no_of_comps
76  call inteflux(prim%each(i), ans)
77  sum1 = sum1 + ans
78  write(*,*)' ', prim%each(i)%symb,' ', sngl(ans),
79  * ' ', sngl(sum1)
80  enddo
81 
82 ! ---------------------
83 
84 
85  call csetcosdeg(1, deg)
86  if(rigc .eq. 0) then
87  write(*,*)
88  * ' primary Int(cos x dI/dE) sum2(cumlative)',
89  * ' sum2/sum1 flat '
90  else
91  write(*,*)
92  * ' primary Int(cos x dI/dE*RigCut) sum2(cumlative)',
93  * ' sum2/sum1 flat'
94  endif
95  sum2 = 0.
96  do i = 1, prim%no_of_comps
97  call inteflux(prim%each(i), ans)
98  sum2 = sum2 + ans
99  write(*,*)' ', prim%each(i)%symb,' ', sngl(ans),
100  * ' ', sngl(sum2), sngl(sum2/sum1)
101  enddo
102 
103 
104  call csetcosdeg(2, deg)
105  if(rigc .eq. 0) then
106  write(*,*)
107  * ' primary Int((1+cos)/2 dI/dE) sum3(cumlative);',
108  * ' sum3/sum1 hemisphere'
109  else
110  write(*,*)
111  * ' primary Int((1+cos)/2 dI/dE*Rigcut) sum3(cumlative);',
112  * ' sum3/sum1 hemisphere'
113  endif
114  sum3 = 0.
115  do i = 1, prim%no_of_comps
116  call inteflux(prim%each(i), ans)
117  sum3 = sum3 + ans
118  write(*,*)' ', prim%each(i)%symb,' ', sngl(ans),
119  * ' ', sngl(sum3), sngl(sum3/sum1)
120  enddo
121  write(*,*)
122  read(*,*) cylr, cylh
123  if(cylr .gt. 0. .and. cylh .gt. 0.) then
124  ratio =2*cylr*cylh/(3.141592*cylr**2)
125  call csetcosdeg(3, deg)
126  if(rigc .eq. 0) then
127  write(*,*)
128  * ' primary Int((cos+ratio*sin)/sqrt(1+ratio**2)dI/dE)',
129  * ' sum3(cumlative); sum3/sum1 cylinder'
130  else
131  write(*,*) ' cylinder: 2rh/pir^2=ratio=',ratio
132  write(*,*)
133  * ' primary Int((cos+ratio*sin)/sqrt(1+ratio**2)dI/dE*Rc)',
134  * ' sum3(cumlative);',' sum3/sum1 cylinder'
135  endif
136  sum3 = 0.
137  do i = 1, prim%no_of_comps
138  call inteflux(prim%each(i), ans)
139  sum3 = sum3 + ans
140  write(*,*)' ', prim%each(i)%symb,' ', sngl(ans),
141  * ' ', sngl(sum3), sngl(sum3/sum1)
142  enddo
143  endif
144  write(*,*)
145  * 'If N primaries are generated in simulation, ST= N/sum1'
146 
147  end
148 ! ******************
149  subroutine inteflux(comp, ans)
150  implicit none
151 #include "Zglobalc.h"
152 #include "Zptcl.h"
153 #include "Zprimary.h"
154  type(component):: comp
155  real*8 ans
156  include 'Zflux.h'
157 
158  real*8 eps, error, ans2, Eth, e_or_p
159  integer icon
160 
161 
162  integer imax
163  external primdn
164  real*8 primdN
165  type(ptcl)::aptcl
166  integer j
167  data eps/1.d-4/
168 
169 
170 
171  compx = comp
172  imax = comp%no_of_seg
173  if(rigc .eq. 0.) then
174 ! no rigidy cut. integrate segmented power functions
175  call inteprim2(comp, 1, imax, ans)
176  if(cosfactor .eq. 0) then
177  if(abs(ans/comp%inte_value-1.d0) .gt. 1.d-3 ) then
178  write(*,*) ' ans=',ans, ' internal integral=',
179  * comp%inte_value
180  stop 9999
181  endif
182  ans = ans * (zen2-zen1)* abs(azmmax - azmmin)*torad
183  elseif(cosfactor .eq. 1) then
184 ! take into account the horizontal area.
185  ans = comp%inte_value * (zen2**2- zen1**2)/2 *
186  * abs(azmmax - azmmin)*torad
187  elseif(cosfactor .eq. 2) then
188 !
189  ans = comp%inte_value
190  * *(zen2-zen1)/2.d0* (1.d0 + (zen1 + zen2)/2)*
191  * abs(azmmax - azmmin)*torad
192  elseif(cosfactor .eq. 3) then
193 ! cos + rat*sin
194  ans = comp%inte_value *
195  * ( (zen2-zen1)/2.0d0 + ratio* (
196  * acos(zen1)-acos(zen2) +
197  * 0.5*( zen2*sqrt(1.-zen2**2) -
198  * zen1*sqrt(1.-zen1**2) )) )
199  * /sqrt(1.+ratio**2)
200  * * abs(azmmax - azmmin)*torad
201  else
202  write(*,*) ' error of cosin'
203  stop 9999
204  endif
205  else
206  call cmkptc(comp%code, comp%subcode, comp%charge, aptcl)
207  eth =sqrt( (rigc*comp%charge)**2 + aptcl%mass**2 )
208  call cconv_prim_e2(comp, eth, e_or_p)
209  call kdwhereis(e_or_p, comp%no_of_seg+1, comp%energy, 1, j)
210  if(j .le. imax) then
211 ! energy integral;
212  call kdexpintfb(primdn, comp%energy(1), comp%energy(j+1),
213  * eps, ans, error, icon)
214  ans = ans * torad
215  endif
216 ! add E> comp.energy(j+1)
217  if(j+1 .lt. imax ) then
218  call inteprim2(comp, j+1, imax, ans2)
219  if(cosfactor .eq. 0 ) then
220  ans2 = ans2 * (zen2- zen1)*abs(azmmax - azmmin)*torad
221  elseif(cosfactor .eq. 1) then
222  ans2 = ans2 * (zen2**2- zen1**2)/2 *
223  * abs(azmmax - azmmin)*torad
224  elseif(cosfactor .eq. 2) then
225 ! fai is from 0 to 2pi;
226  ans2 =
227  * ans2 *(zen2-zen1)/2.d0* (1.d0 + (zen1 + zen2)/2.d0)
228  else
229  ans2 = ans2 *
230  * ( (zen2-zen1)/2.0d0 + ratio* (
231  * acos(zen1)-acos(zen2) +
232  * 0.5*( zen2*sqrt(1.-zen2**2) -
233  * zen1*sqrt(1.-zen1**2) )) )
234  * /sqrt(1.+ratio**2)
235  * * abs(azmmax - azmmin)*torad
236  endif
237  else
238  ans2 = 0.
239  endif
240  ans = ans + ans2
241  endif
242  end
243 
244 ! ****************************
245  real*8 function primdn(eorp)
246  implicit none
247 
248 #include "Zglobalc.h"
249 #include "Zptcl.h"
250 #include "Zprimary.h"
251 ! primary flux at E
252  real*8 eorp
253 
254  include 'Zflux.h'
255 
256  real*8 funczen, ans
257  external funczen
258 
259  e = eorp
260 ! integral on zenith
261  call k16pgaussleg(funczen, zen1, zen2, 16, ans)
262  primdn = ans
263  end
264 ! ********************************
265  real*8 function funczen(zenx)
266  implicit none
267 #include "Zptcl.h"
268 #include "Zprimary.h"
269 
270  real*8 zenx
271 
272  include 'Zflux.h'
273 
274  integer icon
275  real*8 eps, ans, error
276  real*8 funcazim
277  external funcazim
278  data eps/1.d-5/
279 
280  zen = zenx
281  call kdexpintf(funcazim, azmmin, azmmax, eps, ans, error, icon)
282  funczen = ans
283  end
284 !.......................
285  real*8 function funcazim(azm)
286  implicit none
287 #include "Zglobalc.h"
288 #include "Zptcl.h"
289 #include "Zprimary.h"
290  include 'Zflux.h'
291 
292  type(ptcl):: aptcl
293  real*8 rig, azm, prob, flux, zeny
294 
295 ! E= E_or_P must be converted to rigidity
296  call cconv_prim_e(compx, e, aptcl)
297 
298  rig =sqrt( aptcl%fm%p(4)**2 - aptcl%mass**2)/aptcl%charge
299  if(degree) then
300  zeny = zen/torad
301  else
302  zeny = zen
303  endif
304  call crigcut(azm, zeny, rig, prob)
305 !/////////////
306 ! prob = 1.
307 !/////////////
308  call cprimflux0(compx, e, flux) ! E = EorP
309  if(cosfactor .eq. 0 ) then
310  funcazim = flux * prob
311  elseif(cosfactor .eq. 1) then
312  funcazim = flux * prob * zen
313  elseif(cosfactor .eq. 2) then
314  funcazim = flux * prob *(1.0d0+ zen)/2.d0
315  else
316  funcazim = flux * prob *
317  * (zen + ratio*sqrt(1.-zen**2)) / sqrt(1.+ratio**2)
318  endif
319  end
320  subroutine chooktrace
321  end
322  subroutine chookceren
323  end
324  subroutine chookcerens
325  end
326  subroutine chookcerene
327  end
328  subroutine chookbgrun
329  end
330 ! ***************************
331  subroutine csetcosdeg(cosin, degin)
332  implicit none
333  logical cosin, degin
334 #include "Zptcl.h"
335 #include "Zprimary.h"
336 
337  include 'Zflux.h'
338 
339  cosfactor = cosin
340  degree = degin
341  end
342 
343 
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
zen2
real zen2
Definition: Zflux.h:1
primdn
real *8 function primdn(eorp)
Definition: getST.f:246
zen1
real zen1
Definition: Zflux.h:1
e
dE dx *! Nuc Int sampling table e
Definition: cblkMuInt.h:130
ans
atmos%rho(atmos%nodes) **exp(-(z-atmos%z(atmos%nodes))/Hinf) elseif(z .lt. atmos%z(1)) then ans=atmos%rho(1) **exp((atmos%z(1) -z)/atmos%H(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) a=atmos%a(i) if(a .ne. 0.d0) then ans=atmos%rho(i) **(1+a *(z-atmos%z(i))/atmos%H(i)) **(-1.0d0-1.d0/a) else ans=*atmos%rho(i) *exp(-(z-atmos%z(i))/atmos%H(i)) endif endif ! zsave=z ! endif cvh2den=ans end ! ---------------------------------- real *8 function cvh2temp(z) implicit none ! vettical height to temperatur(Kelvin) real *8 z ! input. vertical height in m ! output is temperature of the atmospher in Kelvin real *8 ans integer i if(z .gt. atmos%z(atmos%nodes)) then ans=atmos%T(atmos%nodes) elseif(z .lt. atmos%z(1)) then ans=atmos%T(1)+atmos%b(1) *(z - atmos%z(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) ans=atmos%T(i)+atmos%b(i) *(z-atmos%z(i)) endif cvh2temp=ans end !--------------------------------------------- real *8 function cthick2h(t) implicit none real *8 t ! input. air thickness in kg/m^2 real *8 logt, ans integer i real *8 dod0, fd, a logt=log(t) if(t .ge. atmos%cumd(1)) then ans=atmos%z(1) - *(logt - atmos%logcumd(1)) *atmos%H(1) elseif(t .le. atmos%cumd(atmos%nodes)) then ans=atmos%z(atmos%nodes) - *Hinf *log(t/atmos%cumd(atmos%nodes)) else call kdwhereis(t, atmos%nodes, atmos%cumd, 1, i) ! i is such that X(i) > x >=x(i+1) ans
Definition: cstdatmos0.fNew2.h:82
funcazim
real *8 function funcazim(azm)
Definition: getST.f:286
i
nodes i
Definition: cstdatmos0.fNew2.h:93
crigcut
subroutine crigcut(azmin, zen, rig, prob)
Definition: crigCut.f:6
cconv_prim_e2
subroutine cconv_prim_e2(comp, E, e_or_p)
Definition: csampPrimary.f:230
kdwhereis
subroutine kdwhereis(x, in, a, step, loc)
Definition: kdwhereis.f:27
azmmin
real * azmmin
Definition: Zflux.h:1
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
cprintprim
subroutine cprintprim(out)
Definition: cprintPrim.f:3
cconv_prim_e
subroutine cconv_prim_e(comp, e_or_p, aPtcl)
Definition: csampPrimary.f:128
funczen
real *8 function funczen(zenx)
Definition: getST.f:266
csetcosdeg
subroutine csetcosdeg(cosin, degin)
Definition: cprimAcceptance.f:104
chookceren
subroutine chookceren
Definition: det2Exyz.f:63
rigc
real rigc
Definition: Zflux.h:1
creadparam
subroutine creadparam(io)
Definition: creadParam.f:5
chooktrace
subroutine chooktrace
Definition: chook.f:275
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
kdexpintf
subroutine kdexpintf(func, a, b, eps, ans, error, icon)
Definition: kdexpIntF.f:55
k16pgaussleg
subroutine k16pgaussleg(func, a, b, n, ans)
Definition: k16pGaussLeg.f:20
kdexpintfb
subroutine kdexpintfb(func, a, b, eps, ans, error, icon)
Definition: kdexpIntFb.f:10
chookcerene
subroutine chookcerene
Definition: det2Exyz.f:67
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
coszenith
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p *Zfirst pos xyz Zfirst pos xyz *Zfirst pos *Zfirst pos Zfirst pos *Zfirst pos *Zfirst Zfirst vec w *Zfirst vec w Zfirst vec coszenith
Definition: ZavoidUnionMap.h:1
chookcerens
subroutine chookcerens(no, primary, angle)
Definition: ctemplCeren.f:19
inteprim2
subroutine inteprim2(comp, i1, i2, ans)
Definition: intePrim2.f:2
cprimflux0
subroutine cprimflux0(comp, e_or_p, flux)
Definition: cprimFlux0.f:2
inteflux
subroutine inteflux(comp, ans)
Definition: getST.f:150
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
azmmax
real azmmax
Definition: Zflux.h:1
ptcl
Definition: Zptcl.h:75
cbeginrun
subroutine cbeginrun
Definition: cbeginRun.f:7
zen
real zen
Definition: Zflux.h:1
chookbgrun
subroutine chookbgrun
Definition: chook.f:15
degree
real zen logical degree
Definition: Zflux.h:1