COSMOS v7.655  COSMOSv7655
(AirShowerMC)
mknrfai.f
Go to the documentation of this file.
1 ! suppose you have .nrfai file with invalid "rec" component while
2 ! "all" and "dE/dx" part are correct.
3 ! Also you have size reduced .dat file.
4 ! This program read .dat file and count data in each (r, fai) bin
5 ! and make correct .nrfai data.
6 ! The program structure is the same as reduceSize.f but we accept
7 ! all particles in the .dat
8 ! input file: .nrfai data. give it by environmental
9 ! varialbe NRFAIFILE
10 ! .dat file main output from sim.
11 ! standard input.
12 ! output file: rfaifile with .nrfai-r extension.
13 
14 #include "../FleshHist/asinfo.f"
15 #include "../FleshHist/asdensity.f"
16 #include "../FleshHist/crecprob.f"
17 
18  implicit none
19 #include "Zmaxdef.h"
20 #include "Zobs.h"
21 #include "../FleshHist/Zprivate0.h"
22 #include "./realLimit.h"
23  integer ldep, code, subcode,
24  * charge, ridx, faiidx, codex
25  real rinmu, fai, Ek, time,
26  * wx, wy, wz
27  real*8 u
28  integer ndepth
29  real E0
30  parameter(ndepth= nsites)
31  real nrfaiRec0(nrbin, nfai, 4, ndepth)
32  real rnrfaiRec(nrbin, nfai, 4, ndepth)
33  real pnrfaiRec(nrbin, nfai, 4, ndepth)
34  real nrfaiAll0(nrbin, nfai, 4, ndepth)
35  real dErfai(nrbin, nfai, ndepth)
36  real recprob(nrbin, 4, ndepth)
37 
38  integer EvNo0
39 
40  real limit(4), Emin
41  integer maxsites, nr, newfmt
42  real rat, all, rec, prob
43  integer klena
44  real intdep(ndepth)
45  integer nrbina, nfaia, ansites0
46  integer leng, i, j, k, l
47  integer i0, j0, k0, l0
48  integer NN
49  integer icon0, iconx, icont
50  character*128 input0
51  character*100 nrfaifile, nrfaifile2
52  character*3 id
53  character*5 id5
54  character*20 field(15)
55  integer fnonrfai, fnonrfai2
56  real cosz, age, sum, Nx, depth
57  real nptcls(nrbin, 4, ndepth)
58  integer indivdep(ndepth), packeddepidx(ndepth), depidx
59 ! real rnptcls(nrbin, 4, ndepth)
60  integer icon, kgetenv2, nrec
61 
62 
63  fnonrfai=11
64  fnonrfai2=21
65  nrfaifile=" "
66  leng = kgetenv2("NRFAIFILE", nrfaifile)
67  if(leng .le. 0) then
68  write(0,*) "Env. NRFAIFILE not given"
69  stop 11111
70  endif
71  call copenfw2(fnonrfai, nrfaifile, 1, icon)
72  if(icon .ne. 1) then
73  write(0,*) ' error cannot open', nrfaifile
74  stop 0000
75  endif
76  nrfaifile2 =" "
77  nrfaifile2 = nrfaifile(1:leng)//"-r"
78  call copenfw2(fnonrfai2, nrfaifile2, 1, icon)
79  if(icon .ne. 0) then
80  write(0,*) ' error ', nrfaifile2, ' cannot be created'
81  stop 1235
82  endif
83 
84  do k = 1, ndepth
85  packeddepidx(k)=0
86  do j = 1, 4
87  do l = 1,nfai
88  do i = 1, nrbin
89  rnrfairec(i, l, j, k)=0
90  enddo
91  enddo
92  enddo
93  enddo
94 !
95  do while(.true.)
96  input0 = ' '
97  read(fnonrfai, '(a)', end=1000 ) input0
98  if(input0 .ne. " ") then
99  call ksplit(input0, 20, 15, field, nr)
100  if(nr .eq. 8) then
101  read(input0(1:klena(input0)), *)
102  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0
103  maxsites = ansites0
104  emin=500.d-6
105  newfmt = 0
106  elseif(nr .eq. 9 ) then
107  read(input0(1:klena(input0)), *)
108  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0,
109  * maxsites
110  emin=500.d-6
111  newfmt = 1
112  elseif( nr .eq. 13) then
113  read(input0(1:klena(input0)), *)
114  * evno0, e0,nn, cosz, emin, nrbina, nfaia, ansites0,
115  * maxsites, limit
116  newfmt = 2
117  endif
118  write(0,*) input0
119  limit(1) = reallimitg
120  limit(2) = reallimite
121  limit(3) = reallimitmu
122  limit(4) = reallimith
123 
124  if(nrbina .ne. nrbin .or. nfaia .ne. nfai) then
125  write(0,*)' nrbina=',nrbina, 'or nfaia=',nfaia,
126  * ' differ from the def. in this prog'
127  stop 5555
128  endif
129  if(ansites0 .gt. ndepth) then
130  write(0,*) ' too many depths'
131  stop 6666
132  endif
133 ! ********
134  do i = 1, ansites0
135  do j = 1, 4
136  do k = 1, nfai
137  read(fnonrfai, '(a, f7.1, 4i4)' )
138  * id, intdep(i), l0, i0, j0, k0
139  indivdep(i)=l0
140  packeddepidx(l0)=i ! original dep index to packed indx
141 !
142 ! when the above is written
143 ! l = indivdep(i)
144 ! write(fnonrfai, '("rec",f7.1, 4i4)' )
145 ! * ASDepthList(l)*0.1, l, i, j, k
146 
147 
148  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
149  write(0,*) ' intdep, i0,j0,k0=',
150  * intdep(i), i0, j0, k0, ' strange'
151  stop 8888
152  endif
153  if( id .ne. "rec") then
154  write(0,*) 'id=',id, ' strange'
155  stop 5678
156  endif
157 !
158  read(fnonrfai, *)
159  * ( nrfairec0(l,k,j,i), l=1,nrbin )
160  enddo
161  enddo
162  enddo
163 ! ************
164  do i = 1, ansites0
165  do j = 1, 4
166  do k = 1, nfai
167  read(fnonrfai, '(a,f7.1, 4i4)' )
168  * id, intdep(i), l0, i0, j0, k0
169  indivdep(i)=l0
170  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
171  write(0,*) ' intdep, i0,j0,k0=',
172  * intdep(i), i0, j0, k0, ' strange'
173  stop 9876
174  endif
175  if( id .ne. "all") then
176  write(0,*) 'id=',id,' strange'
177  stop 9999
178  endif
179  read(fnonrfai, *)
180  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
181  enddo
182  enddo
183  enddo
184 ! dErfai
185  do i = 1, ansites0
186  do k = 1, nfai
187  read(fnonrfai, '(a,f7.1, 3i4)' )
188  * id5, intdep(i), l0, i0, k0
189  indivdep(i)=l0
190  if(i0 .ne. i .or. k .ne. k0) then
191  write(0,*) ' intdep, i0,k0=',
192  * intdep(i), i0, k0, ' strange'
193  stop 98765
194  endif
195  if( id5 .ne. "dE/dx") then
196  write(0,*) 'id=',id5,' strange'
197  stop 9999
198  endif
199  read(fnonrfai, *)
200  * ( derfai(l,k,i), l=1,nrbin )
201  enddo
202  enddo
203  else
204 !
205  write(0,*) ' all nrfai data has been read'
206  if(newfmt .eq. 0) then
207  write(0, '(i2, 1pE11.3,i3, 0pf7.1, 1pE11.3,0p, 3i4)')
208  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0
209  elseif(newfmt .eq. 1) then
210  write(0, '(i2, 1pE11.3,i3, 0pf7.1, 1pE11.3,0p, 4i4)')
211  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0,
212  * maxsites
213  else
214  write(0,
215  7 '(i2, 1pE11.3,i3, 0pf7.1, 1pE11.3, 0p,4i4, 4f10.0)')
216  * evno0, e0,nn, cosz, emin, nrbina, nfaia, ansites0,
217  * maxsites, limit
218  endif
219  endif
220 
221  enddo
222  1000 continue
223  input0= ' '
224 !************* reset limit
225 ! limit=min(4000.0, limit)
226  if(newfmt .ne. 2) then
227  write(0,*) ' limit=',limit(1), ' is being used'
228  else
229  write(0,*) ' limit=',limit, ' are being used'
230  endif
231 !************
232  do i = 1, ansites0
233  do j = 1, 4
234  do k = 1,nfai
235  do l = 1, nrbin
236  if(nrfairec0(l, k, j, i) .gt. limit(j)) then
237 ! accept with this prob.
238  pnrfairec(l, k, j, i)=
239  * limit(j)/nrfairec0(l, k, j, i)
240  else
241  pnrfairec(l, k, j, i)=1.0
242  endif
243  enddo
244 !///////
245  write(0, '(f7.1, 4i4,a)' )
246  * intdep(i), indivdep(i), i, j, k, ' prob'
247  write(0, '(1p10E11.3)')
248  * (pnrfairec(l,k,j,i), l=1, nrbin)
249 !c/////
250  enddo
251  enddo
252  enddo
253 ! -------------
254 ! main input data; header; obsolete
255 ! read(*,'(a)') input0
256 ! write(*,'(a)') input0(1:klena(input0))
257  nrec= 0
258  do while(.true.)
259  read(*,'(a)', end=100, err=500) input0
260  if( index(input0(1:2), "i") .gt. 0 ) cycle
261 ! read(*, *, end=100, Err=500)
262  read(input0,*)
263  * ldep, code, subcode,
264  * charge, ridx, faiidx,
265  * rinmu, fai,
266  * ek, time,
267  * wx, wy, wz
268  nrec= nrec+1
269  depidx = packeddepidx(ldep)
270 
271  if(depidx .le. 0) then
272  write(0,*) ' should not happen. depidx=',depidx
273  write(0,*) ' ldep=',ldep, ' code=',code, 'nrec=',nrec
274  stop 9875
275  endif
276  codex=min(code, 4)
277 ! if(u .lt.
278 ! * pnrfaiRec(ridx, faiidx, codex, depidx) ) then
279  rnrfairec(ridx, faiidx, codex, depidx)=
280  * rnrfairec(ridx, faiidx, codex, depidx) + 1
281  enddo
282  100 continue
283 
284  if( newfmt .eq. 0) then
285  write( fnonrfai2,
286  * '(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3,0p,3i4)' )
287  * evno0, e0, nn, cosz, limit(1), nrbin, nfai, ansites0
288  elseif( newfmt .eq. 1) then
289  write( fnonrfai2,
290  * '(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3,0p,4i4)' )
291  * evno0, e0, nn, cosz, limit(1), nrbin, nfai, ansites0,
292  * maxsites
293  else
294  write( fnonrfai2,
295  * '(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3,0p, 4i4,4f10.0)' )
296  * evno0, e0, nn, cosz, emin, nrbin, nfai, ansites0,
297  * maxsites, limit
298  endif
299 
300 
301  do i = 1, ansites0
302  do j = 1, 4
303  do k = 1, nfai
304  l = indivdep(i)
305  write(fnonrfai2, '("rec",f7.1, 4i4)' )
306  * intdep(i), l, i, j, k
307  write(fnonrfai2, '(1p10E11.3)')
308  * ( rnrfairec(l,k,j,i), l=1,nrbin )
309  enddo
310  enddo
311  enddo
312  do i = 1, ansites0
313  do j = 1, 4
314  do k = 1, nfai
315  l = indivdep(i)
316  write(fnonrfai2, '("all",f7.1, 4i4)' )
317  * intdep(i), l, i, j, k
318  write(fnonrfai2, '(1p10E11.3)')
319  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
320  enddo
321  enddo
322  enddo
323 ! dErfai
324  do i = 1, ansites0
325  do k = 1, nfai
326  l = indivdep(i)
327  write(fnonrfai2, '("dE/dx",f7.1, 3i4)' )
328  * intdep(i), l, i, k
329  write(fnonrfai2, '(1p10E11.3)')
330  * ( derfai(l,k,i), l=1,nrbin )
331  enddo
332  enddo
333 
334  write(0,*) 'end of run'
335  write(fnonrfai2,*)
336  stop
337  500 continue
338  write(0,*) ' input error at record =', nrec
339  read(*,'(a)') input0
340  write(0,*) input0
341  end
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
ksplit
subroutine ksplit(a, m, n, b, nr)
Definition: ksplit.f:2
kgetenv2
integer function kgetenv2(envname, envresult)
Definition: cgetLoginN.f:77
i
nodes i
Definition: cstdatmos0.fNew2.h:93
time
subroutine time(xxx)
Definition: chook.f:5
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
leng
integer leng
Definition: interface2.h:1
nrbin
! timing nrbin
Definition: Zprivate2.h:12
j
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
copenfw2
subroutine copenfw2(io, fnin, form, icon)
Definition: copenf.f:205
code
*Zfirst p fm *Zfirst p Zfirst p code
Definition: ZavoidUnionMap.h:1
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
charge
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p charge
Definition: ZavoidUnionMap.h:1
klena
integer function klena(cha)
Definition: klena.f:20
nfai
! timing nfai
Definition: Zprivate2.h:12
subcode
*Zfirst p fm *Zfirst p Zfirst p Zfirst p subcode
Definition: ZavoidUnionMap.h:1