COSMOS v7.655  COSMOSv7655
(AirShowerMC)
reducebinSize.f
Go to the documentation of this file.
1 #define FNODATDEF 33
2 ! reduce the file size; exmine the number of particles (N)
3 ! in each (r,fai) bin and if it exceeds a limit,
4 ! each partilce is discarded with the probabilty of
5 ! 1-limit/N.
6 ! input file: .nrfai data (assembled).
7 ! give it by environmental varialbe NRFAIFILE
8 ! .dat file ( main output from sim. )
9 ! give it by environmental variable DATFILE
10 ! output file: size reduced ascii output file: stdout.
11 ! and .nrfai file for that; name will be
12 ! such that: let .dat file name is
13 ! xxx.dat, then xxx.nrfai.
14 
15 #include "../FleshHist/asinfo.f"
16 #include "../FleshHist/asdensity.f"
17 #include "../FleshHist/crecprob.f"
18  implicit none
19 #include "../FleshHist/Zprivate0.h"
20 #include "../FleshHist/Zprivate1.h"
21 #include "../FleshHist/Zprivate3.h"
22 
23  type(buffer):: abuf
24 
25  integer ldep, code, subcode,
26  * charge, ridx, faiidx, codex
27  real rinmu, fai, Ek, time,
28  * wx, wy, wz
29  real*8 u
30  integer ndepth
31  real E0
32  parameter(ndepth= nsites)
33  real nrfaiRec0(nrbin, nfai, 4, ndepth)
34  real rnrfaiRec(nrbin, nfai, 4, ndepth)
35  real pnrfaiRec(nrbin, nfai, 4, ndepth)
36  real nrfaiAll0(nrbin, nfai, 4, ndepth)
37  real recprob(nrbin, 4, ndepth)
38 
39  integer EvNo0
40 
41  real limit
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, datfile
52  character*3 id
53  integer fnonrfai, fnonrfai2
54  real cosz, age, sum, Nx, depth
55  real nptcls(nrbin, 4, ndepth)
56  integer indivdep(ndepth), packeddepidx(ndepth), depidx
57 ! real rnptcls(nrbin, 4, ndepth)
58  integer icon, kgetenv2, nrec
59 
60 
61  fnonrfai=11
62  fnonrfai2=21
63  leng = kgetenv2("NRFAIFILE", nrfaifile)
64  if(leng .le. 0) then
65  write(0,*) "Env. NRFAIFILE not given"
66  stop 11111
67  endif
68  call copenfw2(fnonrfai, nrfaifile, 1, icon)
69  if(icon .ne. 1) then
70  write(0,*) ' error cannot open', nrfaifile
71  stop 0000
72  endif
73  datfile=' '
74  leng = kgetenv2("DATFILE", datfile)
75  if(leng .le. 0) then
76  write(0,*) "Env. DATFILE not given"
77  stop 11112
78  endif
79  call copenfw2(fnodat, datfile, 2, icon)
80  if(icon .ne. 1) then
81  write(0,*) ' error cannot open', datfile
82  stop 0002
83  endif
84  nrfaifile2 = datfile(1:leng-3)//"rnrfai"
85  call copenfw2(fnonrfai2, nrfaifile2, 1, icon)
86  write(0,*) ' reduced .nrfai is', trim(nrfaifile2)
87  write(0,*) ' open icon=', icon
88  do k = 1, ndepth
89  packeddepidx(k)=0
90  do j = 1, 4
91  do l = 1,nfai
92  do i = 1, nrbin
93  rnrfairec(i, l, j, k)=0
94  enddo
95  enddo
96  enddo
97  enddo
98 !
99  do while(.true.)
100  input0 = ' '
101  read(fnonrfai, '(a)', end=1000 ) input0
102  if(input0 .ne. " ") then
103  read(input0(1:klena(input0)), *)
104  * evno0, e0,nn, cosz, limit, nrbina, nfaia, ansites0
105 
106  write(0,*) input0
107 
108  if(nrbina .ne. nrbin .or. nfaia .ne. nfai) then
109  write(0,*)' nrbina=',nrbina, 'or nfaia=',nfaia,
110  * ' differ from the def. in this prog'
111  stop 5555
112  endif
113  if(ansites0 .gt. ndepth) then
114  write(0,*) ' too many depths'
115  stop 6666
116  endif
117 ! ********
118  do i = 1, ansites0
119  do j = 1, 4
120  do k = 1, nfai
121  read(fnonrfai, '(a, f7.1, 4i4)' )
122  * id, intdep(i), l0, i0, j0, k0
123  indivdep(i)=l0
124  packeddepidx(l0)=i ! original dep index to packed indx
125 
126 !
127 ! when the above is written
128 ! l = indivdep(i)
129 ! write(fnonrfai, '("rec",f7.1, 4i4)' )
130 ! * ASDepthList(l)*0.1, l, i, j, k
131 
132 
133  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
134  write(0,*) ' intdep, i0,j0,k0=',
135  * intdep(i), i0, j0, k0, ' strange'
136  stop 8888
137  endif
138  if( id .ne. "rec") then
139  write(0,*) 'id=',id, ' strange'
140  stop 5678
141  endif
142 !
143  read(fnonrfai, *)
144  * ( nrfairec0(l,k,j,i), l=1,nrbin )
145  enddo
146  enddo
147  enddo
148 ! ************
149  do i = 1, ansites0
150  do j = 1, 4
151  do k = 1, nfai
152  read(fnonrfai, '(a,f7.1, 4i4)' )
153  * id, intdep(i), l0, i0, j0, k0
154  indivdep(i)=l0
155  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
156  write(0,*) ' intdep, i0,j0,k0=',
157  * intdep(i), i0, j0, k0, ' strange'
158  stop 9876
159  endif
160  if( id .ne. "all") then
161  write(0,*) 'id=',id,' strange'
162  stop 9999
163  endif
164  read(fnonrfai, *)
165  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
166  enddo
167  enddo
168  enddo
169  else
170 !
171  write(0,*) ' all nrfai data has been read'
172  write(0, '(i2, 1pE11.3,i3, 0pf7.1, 1pE11.3, 3i4)')
173  * evno0, e0,nn, cosz, limit, nrbina, nfaia, ansites0
174  endif
175  enddo
176  1000 continue
177  input0= ' '
178 !************* reset limit
179 ! limit=min(4000.0, limit)
180  write(0,*) ' limit=',limit, ' is being used'
181 !************
182  do i = 1, ansites0
183  do j = 1, 4
184  do k = 1,nfai
185  do l = 1, nrbin
186  if(nrfairec0(l, k, j, i) .gt. limit) then
187 ! accept with this prob.
188  pnrfairec(l, k, j, i)=
189  * limit/nrfairec0(l, k, j, i)
190  else
191  pnrfairec(l, k, j, i)=1.0
192  endif
193  enddo
194 !///////
195 ! write(0, '(f7.1, 4i4,a)' )
196 ! * intdep(i), indivdep(i), i, j, k, ' prob'
197 ! write(0, '(1p10E11.3)')
198 ! * (pnrfaiRec(l,k,j,i), l=1, nrbin)
199 !c/////
200  enddo
201  enddo
202  enddo
203 ! -------------
204 ! main input data; header
205  call binreadhead
206  nrec= 0
207 
208  do while(.true.)
209  call binread1(abuf, icon)
210  if(icon .ne. 0) exit
211  nrec= nrec+1
212  depidx = packeddepidx(abuf.ldep)
213 
214  if(depidx .le. 0) then
215  write(0,*) ' should not happen. depidx=',depidx
216  write(0,*) 'nrec=',nrec, 'icon=',icon, ' ldep=',
217  * abuf.ldep, ' code=',abuf.code
218  stop 9875
219  endif
220 
221  call rndc(u)
222  codex=min(abuf.code, 4)
223  if(u .lt.
224  * pnrfairec(abuf.ridx, abuf.faiidx, codex, depidx) ) then
225  rnrfairec(abuf.ridx, abuf.faiidx, codex, depidx)=
226  * rnrfairec(abuf.ridx, abuf.faiidx, codex, depidx) + 1
227 
228  if(abuf.wz .lt. 0.999) then
229  write(*,
230  * '(6i3, 1pE11.3, 0p,f6.1,1p2E11.3,0p, 2f8.4,f10.6)')
231  * abuf.ldep, abuf.code, abuf.subcode,
232  * abuf.charge, abuf.ridx, abuf.faiidx,
233  * abuf.rinmu, abuf.fai,
234  * abuf.ek, abuf.t,
235  * abuf.wx, abuf.wy, abuf.wz
236  else
237  write(*,'(6i3,1pE11.3, 0p, f6.1, 1p2E11.3, a)')
238  * abuf.ldep, abuf.code, abuf.subcode,
239  * abuf.charge, abuf.ridx, abuf.faiidx,
240  * abuf.rinmu, abuf.fai,
241  * abuf.ek, abuf.t,
242  * ' 0 0 1'
243  endif
244  endif
245  enddo
246  100 continue
247  write( fnonrfai2,
248  * '(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3,3i4)' )
249  * evno0, e0, nn, cosz, limit, nrbin, nfai, ansites0
250 
251  do i = 1, ansites0
252  do j = 1, 4
253  do k = 1, nfai
254  l = indivdep(i)
255  write(fnonrfai2, '("rec",f7.1, 4i4)' )
256  * intdep(i), l, i, j, k
257  write(fnonrfai2, '(1p10E11.3)')
258  * ( rnrfairec(l,k,j,i), l=1,nrbin )
259  enddo
260  enddo
261  enddo
262  do i = 1, ansites0
263  do j = 1, 4
264  do k = 1, nfai
265  l = indivdep(i)
266  write(fnonrfai2, '("all",f7.1, 4i4)' )
267  * intdep(i), l, i, j, k
268  write(fnonrfai2, '(1p10E11.3)')
269  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
270  enddo
271  enddo
272  enddo
273  write(0,*) 'end of run'
274  end
275 ! ************************
276  subroutine binreadhead
277 !
278 #include "../FleshHist/Zprivate1.h"
279 #include "../FleshHist/Zprivate3.h"
280  real*8 Et, wx, wy, wz
281  integer EventNo
282  integer*2 code
283  integer i, icon
284 
285  type(buffer):: abuf
286  save
287 
288  read(fnodat)
289  * eventno, code, et, wx, wy, wz
290  write(*,'("i ", i3, i4, g13.4,3f11.7)')
291  * eventno, code, et, wx, wy, wz
292 
293  bufc = 0
294  return
295 
296 ! ******************
297  entry binread1(abuf, icon)
298 ! ***************
299  if(bufc .eq. 0) then
300  read(fnodat,end=100, err=500) bufc, (buf(i), i=1, bufc)
301 ! ****************************
302 ! bufc usage was wrong for the first 10^19 eV event made by
303 ! olixxx. it is 0~99999; while correct one should be
304 ! 1~100000. we neglect data at 0 and 100000.
305  if(bufc .eq. 100000) bufc=bufc-1
306  endif
307  if(bufc .gt. 0) then
308  abuf= buf(bufc)
309  bufc = bufc -1
310  icon = 0
311  else
312  write(0, *) ' strange bufc'
313  stop 99999
314  endif
315  return
316  100 continue
317  icon = 1
318  return
319  500 continue
320  write(0,*) ' input data read err'
321  stop
322  end
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
kgetenv2
integer function kgetenv2(envname, envresult)
Definition: cgetLoginN.f:77
i
nodes i
Definition: cstdatmos0.fNew2.h:93
binreadhead
subroutine binreadhead
Definition: reducebinSize.f:277
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
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
t
nodes t
Definition: cstdatmos0.fNew2.h:84
nfai
! timing nfai
Definition: Zprivate2.h:12
subcode
*Zfirst p fm *Zfirst p Zfirst p Zfirst p subcode
Definition: ZavoidUnionMap.h:1