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