COSMOS v7.655  COSMOSv7655
(AirShowerMC)
reduceEachSize.bin.f
Go to the documentation of this file.
1 #define FNODATDEF 33
2 #include "./realLimit.h"
3 
4 ! reduce the file size; exmine the number of particles (N)
5 ! in each (r,fai) bin and if it exceeds a limit,
6 ! each partilce is discarded with the probabilty of
7 ! 1-limit/N.
8 ! input file: combined .nrfai data. give it by environmental
9 ! varialbe NRFAIFILE
10 ! each .dat file (binary file) main output from sim.
11 ! give it by env. vari. DATFILE
12 ! env. vari. HEADER for the first set give it "yes"
13 ! else 'no'.
14 ! output file: size reduced file corresponding to .dat; this is put stdout (ascii)
15 ! so this must be concatinated to get combied .dat file.
16 ! nrfaifile for each .dat; -r is attached as .nrfai-r.
17 ! you must combine .nrfai-r files;
18 ! rm old each .nrfai file. and rename .nrfai-r to .nrfai
19 ! and use assemNrfai.csh in Assemble.
20 ! ********"all" data in .nrfai-r is really the "all" so this one should not be
21 ! added any more in assemNrfai
22 #include "../FleshHist/asinfo.f"
23 #include "../FleshHist/asdensity.f"
24 #include "../FleshHist/crecprob.f"
25  implicit none
26 #include "Zmaxdef.h"
27 #include "Zobs.h"
28 #include "../FleshHist/Zprivate0.h"
29 #include "../FleshHist/Zprivate1.h"
30 #include "../FleshHist/Zprivate3.h"
31 
32  integer ldep,
33  * charge, ridx, faiidx, codex
34 
35  real*8 u
36  integer ndepth
37  integer EventNo
38  integer*2 code
39  real*8 Et, wx, wy, wz
40  real E0
41  parameter(ndepth= nsites)
42  real nrfaiRec0(nrbin, nfai, 4, ndepth)
43  real rnrfaiRec(nrbin, nfai, 4, ndepth)
44  real pnrfaiRec(nrbin, nfai, 4, ndepth)
45  real nrfaiAll0(nrbin, nfai, 4, ndepth)
46  real dErfai(nrbin, nfai, ndepth)
47  real dErfai2(nrbin, nfai, ndepth)
48  real recprob(nrbin, 4, ndepth)
49  character*20 field(15)
50  integer EvNo0
51  integer l0(ndepth), recl0(ndepth)
52  real limit(4), Emin
53  real rat, all, rec
54  integer klena
55  real intdep(ndepth), recdep(ndepth)
56 
57 
58 
59  integer nrbina, nfaia, ansites0
60  integer leng, i, j, k, l, leng2, lengh
61  integer maxsites, nr
62  integer i0, j0, k0, reci0, recj0, reck0
63  integer NN
64  integer icon0, iconx, icont
65  character*128 input0
66  character*100 nrfaifile, nrfaifile2, datfile, datfile2
67  character*3 id
68  character*5 id5
69 ! character*5 header
70  integer fnonrfai, fnonrfai2, fnodat2
71  real cosz, age, sum, Nx, depth, wgt, wwgt, prob
72  real nptcls(nrbin, 4, ndepth)
73  integer indivdep(ndepth), packeddepidx(ndepth), depidx
74 ! real rnptcls(nrbin, 4, ndepth)
75  integer icon, kgetenv2, nrec
76  logical SeeLowdE, accept
77  integer newfmt
78  fnonrfai=11
79  fnonrfai2=21
80 ! fnodat = 6
81 ! fnodat2 = 8
82  nrfaifile=" "
83  leng = kgetenv2("NRFAIFILE", nrfaifile)
84  if(leng .le. 0) then
85  write(0,*) "Env. NRFAIFILE not given"
86  stop 11111
87  endif
88  call copenfw2(fnonrfai, nrfaifile, 1, icon)
89  if(icon .ne. 1) then
90  write(0,*) ' error cannot open', nrfaifile
91  stop 0000
92  endif
93  leng = kgetenv2("SeeLowdE", datfile)
94  if(leng .le. 0) then
95  write(0,*) ' SeeLowdE not given'
96  stop
97  endif
98  seelowde = datfile(1:leng) .eq. "yes"
99 
100  leng2 = kgetenv2("DATFILE", datfile)
101  if(leng2 .le. 0) then
102  write(0,*) "Env. DATFILE not given"
103  stop 11112
104  endif
105  call copenfw2(fnodat, datfile, 2, icon)
106  if(icon .ne. 1) then
107  write(0,*) ' input main data not exist'
108  write(0,*) ' file=',datfile
109  stop 1113
110  endif
111 
112  datfile2 =" "
113  datfile2 = datfile(1:leng2)//"-r"
114  nrfaifile2 =" "
115  nrfaifile2 = datfile(1:leng2-3)//"nrfai-r"
116  call copenfw2(fnonrfai2, nrfaifile2, 1, icon)
117  if(icon .ne. 0) then
118  write(0,*) ' error ', nrfaifile2, ' cannot be created'
119  stop 1235
120  endif
121 ! lengh= kgetenv2("HEADER", header)
122 
123  do k = 1, ndepth
124  packeddepidx(k)=0
125  do j = 1, 4
126  do l = 1,nfai
127  do i = 1, nrbin
128  rnrfairec(i, l, j, k)=0
129  enddo
130  enddo
131  enddo
132  enddo
133 !
134  do while(.true.)
135  input0 = ' '
136  read(fnonrfai, '(a)', end=1000 ) input0
137  if(input0 .ne. " ") then
138  call ksplit(input0, 20, 15, field, nr)
139  if(nr .eq. 8) then
140  read(input0(1:klena(input0)), *)
141  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0
142  maxsites = ansites0
143  emin = 500.d-6
144  newfmt = 0
145  elseif(nr .eq. 9 ) then
146  read(input0(1:klena(input0)), *)
147  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0,
148  * maxsites
149  emin = 500.d-6
150  newfmt = 1
151  elseif(nr .eq. 13) then
152  read(input0(1:klena(input0)), *)
153  * evno0, e0,nn, cosz, emin, nrbina, nfaia, ansites0,
154  * maxsites, limit
155  newfmt = 2
156  else
157  write(0,*) ' header fmt strage '
158  write(0,*) input0
159  stop 01234
160  endif
161 
162 ! write(0,*) input0
163 
164  if(nrbina .ne. nrbin .or. nfaia .ne. nfai) then
165  write(0,*)' nrbina=',nrbina, 'or nfaia=',nfaia,
166  * ' differ from the def. in this prog'
167  stop 5555
168  endif
169  if(ansites0 .gt. ndepth) then
170  write(0,*) ' too many depths'
171  stop 6666
172  endif
173 ! ********
174  do i = 1, ansites0
175  do j = 1, 4
176  do k = 1, nfai
177  read(fnonrfai, '(a, f7.1, 4i4)' )
178  * id, recdep(i), recl0(i), i0, j0, k0
179  indivdep(i)=recl0(i)
180  packeddepidx(recl0(i))=i ! original dep index to packed indx
181 
182 !
183 ! when the above is written
184 ! l = indivdep(i)
185 ! write(fnonrfai, '("rec",f7.1, 4i4)' )
186 ! * ASDepthList(l)*0.1, l, i, j, k
187 
188 
189  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
190  write(0,*) ' intdep, i0,j0,k0=',
191  * recdep(i), i0, j0, k0, ' strange'
192  stop 8888
193  endif
194  if( id .ne. "rec") then
195  write(0,*) 'id=',id, ' strange'
196  stop 5678
197  endif
198 !
199  read(fnonrfai, *)
200  * ( nrfairec0(l,k,j,i), l=1,nrbin )
201  enddo
202  enddo
203  enddo
204 ! ************
205  do i = 1, maxsites
206  do j = 1, 4
207  do k = 1, nfai
208  read(fnonrfai, '(a,f7.1, 4i4)' )
209  * id, intdep(i), l0(i), i0, j0, k0
210  indivdep(i)=l0(i)
211  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
212  write(0,*) ' intdep, i0,j0,k0=',
213  * intdep(i), i0, j0, k0, ' strange'
214  stop 9876
215  endif
216  if( id .ne. "all") then
217  write(0,*) 'id=',id,' strange'
218  stop 9999
219  endif
220  read(fnonrfai, *)
221  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
222  enddo
223  enddo
224  enddo
225 ! dErfai
226  do i = 1, maxsites
227  do k = 1, nfai
228  read(fnonrfai, '(a,f7.1, 3i4)' )
229  * id5, intdep(i), l0(i), i0, k0
230  indivdep(i)=l0(i)
231  if(i0 .ne. i .or. k .ne. k0) then
232  write(0,*) ' intdep, i0,k0=',
233  * intdep(i), i0, k0, ' strange'
234  stop 98765
235  endif
236  if( id5 .ne. "dE/dx") then
237  write(0,*) 'id=',id5,' strange'
238  stop 9999
239  endif
240  read(fnonrfai, *)
241  * ( derfai(l,k,i), l=1,nrbin )
242  if(seelowde) read(fnonrfai, *)
243  * ( derfai2(l,k,i), l=1,nrbin )
244  enddo
245  enddo
246  else
247 !
248  write(0,*) ' all nrfai data has been read'
249  limit(1) = reallimitg
250  limit(2) = reallimite
251  limit(3) = reallimitmu
252  limit(4) = reallimith
253  if(newfmt .eq. 0) then
254  write(0, '(i2, 1pE11.3,i3, 0pf7.1, 1pE11.3,0p, 3i4)')
255  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0
256  elseif( newfmt .eq. 1 ) then
257  write(0, '(i2, 1pE11.3,i3, 0pf7.1, 1pE11.3, 0p,4i4)')
258  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0,
259  * maxsites
260  else
261  write(0, '(i2, 1pE11.3,i3, 0pf7.1,1pE11.3,0p,4i4,4f8.9)')
262  * evno0, e0,nn, cosz, emin, nrbina, nfaia, ansites0,
263  * maxsites, limit
264  endif
265  endif
266  enddo
267  1000 continue
268  input0= ' '
269 !************* reset limit
270 ! limit=min(4000.0, limit)
271  limit(1) = reallimitg
272  limit(2) = reallimite
273  limit(3) = reallimitmu
274  limit(4) = reallimith
275 
276  if(newfmt .ne. 2) then
277  write(0,*) ' limit=',limit(1), ' is being used'
278  else
279  write(0,*) ' limit=',limit, ' ares being used'
280  endif
281 !************
282  do i = 1, ansites0
283  do j = 1, 4
284  do k = 1,nfai
285  do l = 1, nrbin
286  if(nrfairec0(l, k, j, i) .gt. limit(j)) then
287 ! accept with this prob.
288  pnrfairec(l, k, j, i)=
289  * limit(j)/nrfairec0(l, k, j, i)
290  else
291  pnrfairec(l, k, j, i)=1.0
292  endif
293  enddo
294 !///////
295  write(0, '(f7.1, 4i4,a)' )
296  * recdep(i), indivdep(i), i, j, k, ' prob'
297  write(0, '(1p10E11.3)')
298  * (pnrfairec(l,k,j,i), l=1, nrbin)
299 !c/////
300  enddo
301  enddo
302  enddo
303 ! -------------
304 ! main input data; header Obsolete
305 ! read(fnodat) EventNo, code, Et, wx, wy, wz
306 ! if( header(1:lengh) .eq. "yes") then
307 ! write(*,'("i ", i3, i4, g13.4,3f11.7)')
308 ! * EventNo, code, Et, wx, wy, wz
309 ! endif
310 
311  nrec= 0
312  do while(.true.)
313  read(fnodat,end=100) bufc, (buf(i), i=1, bufc)
314 !/////////////
315 ! write(0,*) ' bufc=',bufc
316 !///////////
317  do i = 1, bufc
318  nrec= nrec+1
319  ldep = buf(i).ldep
320  depidx = packeddepidx(ldep)
321  faiidx= buf(i).faiidx
322  ridx = buf(i).ridx
323  codex = min(buf(i).code, 4)
324  wgt = buf(i).wgt
325  prob = pnrfairec(ridx, faiidx, codex, depidx)
326  if( prob .gt. 1.) then
327  wwgt = wgt
328  accept = .true.
329  else
330  prob = prob * wgt
331  if(prob .gt. 1.) then
332  accept = .true.
333  wwgt = prob
334  else
335  call rndc(u)
336  if(u .lt. prob) then
337  accept = .true.
338  wwgt= 1.
339  else
340  accept = .false.
341  endif
342  endif
343  endif
344  if(accept) then
345  rnrfairec(ridx, faiidx, codex, depidx)=
346  * rnrfairec(ridx, faiidx, codex, depidx) + wwgt
347  write(*,
348  * '(6i3, 1pE11.3, 0p,f6.1,1p2E11.3,0p, 2f8.4,f10.6,1pE11.3)')
349  * buf(i).ldep, buf(i).code, buf(i).subcode,
350  * buf(i).charge, buf(i).ridx, buf(i).faiidx,
351  * buf(i).rinmu, buf(i).fai, buf(i).ek,
352  * buf(i).t, buf(i).wx, buf(i).wy, buf(i).wz,
353  * wwgt
354  endif
355  enddo
356  enddo
357  100 continue
358 
359  write( fnonrfai2,
360  * '(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3,0p,4i4,4f8.0)' )
361  * evno0, e0, nn, cosz, emin, nrbin, nfai, ansites0,
362  * maxsites, limit
363 
364  do i = 1, ansites0
365  do j = 1, 4
366  do k = 1, nfai
367  l = recl0(i)
368  write(fnonrfai2, '("rec",f7.1, 4i4)' )
369  * recdep(i), l, i, j, k
370  write(fnonrfai2, '(1p10E11.3)')
371  * ( rnrfairec(l,k,j,i), l=1,nrbin )
372  enddo
373  enddo
374  enddo
375  do i = 1, maxsites
376  do j = 1, 4
377  do k = 1, nfai
378 ! l = indivdep(i)
379  l = l0(i)
380  write(fnonrfai2, '("all",f7.1, 4i4)' )
381  * intdep(i), l, i, j, k
382  write(fnonrfai2, '(1p10E11.3)')
383  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
384  enddo
385  enddo
386  enddo
387 ! dErfai
388  do i = 1, maxsites
389  do k = 1, nfai
390 ! l = indivdep(i)
391  l = l0(i)
392  write(fnonrfai2, '("dE/dx",f7.1, 3i4)' )
393  * intdep(i), l, i, k
394  write(fnonrfai2, '(1p10E11.3)')
395  * ( derfai(l,k,i), l=1,nrbin )
396  if(seelowde) write(fnonrfai2, '(1p10E11.3)')
397  * ( derfai2(l,k,i), l=1,nrbin )
398  enddo
399  enddo
400 
401  write(0,*) 'end of run'
402  write(fnonrfai2, *)
403  stop
404  500 continue
405  write(0,*) ' input error at record =', nrec
406  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
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
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
t
nodes t
Definition: cstdatmos0.fNew2.h:84
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