COSMOS v7.655  COSMOSv7655
(AirShowerMC)
reduceEachSize.binbyL.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*160 Lfile(ndepth)
68  character*80 baseDir
69  integer lengBaseDir
70  character*32 execID
71  integer lengExecID
72  character*3 id
73  character*5 id5
74 ! character*5 header
75  integer fnonrfai, fnonrfai2, fnodat2
76  real cosz, age, sum, Nx, depth, wgt, wwgt, prob
77  real nptcls(nrbin, 4, ndepth)
78  integer indivdep(ndepth), packeddepidx(ndepth), depidx
79 ! real rnptcls(nrbin, 4, ndepth)
80  integer icon, kgetenv2, nrec
81  logical SeeLowdE, accept
82  integer newfmt
83  fnonrfai=11
84  fnonrfai2=21
85 ! fnodat = 6
86 ! fnodat2 = 8
87  nrfaifile=" "
88  leng = kgetenv2("NRFAIFILE", nrfaifile)
89  if(leng .le. 0) then
90  write(0,*) "Env. NRFAIFILE not given"
91  stop 11111
92  endif
93  call copenfw2(fnonrfai, nrfaifile, 1, icon)
94  if(icon .ne. 1) then
95  write(0,*) ' error cannot open', nrfaifile
96  stop 0000
97  endif
98  leng = kgetenv2("SeeLowdE", datfile)
99  if(leng .le. 0) then
100  write(0,*) ' SeeLowdE not given'
101  stop
102  endif
103  seelowde = datfile(1:leng) .eq. "yes"
104 
105  leng2 = kgetenv2("DATFILE", datfile)
106  if(leng2 .le. 0) then
107  write(0,*) "Env. DATFILE not given"
108  stop 11112
109  endif
110 
111  lengbasedir=kgetenv2("BASEDIR", basedir)
112  if(lengbasedir .eq. 0) then
113  write(0,*) 'BASEDIR not given'
114  stop 8888
115  endif
116  lengexecid=kgetenv2("EXECID", execid)
117  if(lengexecid .eq. 0) then
118  write(0,*) 'EXECID not given'
119  stop 9999
120  endif
121 
122 
123 
124 
125  call copenfw2(fnodat, datfile, 2, icon)
126  if(icon .ne. 1) then
127  write(0,*) ' input main data not exist'
128  write(0,*) ' file=',datfile
129  stop 1113
130  endif
131 
132  datfile2 =" "
133  datfile2 = datfile(1:leng2)//"-r"
134  nrfaifile2 =" "
135  nrfaifile2 = datfile(1:leng2-3)//"nrfai-r"
136  call copenfw2(fnonrfai2, nrfaifile2, 1, icon)
137  if(icon .ne. 0) then
138  write(0,*) ' error ', nrfaifile2, ' cannot be created'
139  stop 1235
140  endif
141 ! lengh= kgetenv2("HEADER", header)
142 
143  do k = 1, ndepth
144  packeddepidx(k)=0
145  do j = 1, 4
146  do l = 1,nfai
147  do i = 1, nrbin
148  rnrfairec(i, l, j, k)=0
149  enddo
150  enddo
151  enddo
152  enddo
153 !
154  do while(.true.)
155  input0 = ' '
156  read(fnonrfai, '(a)', end=1000 ) input0
157  if(input0 .ne. " ") then
158  call ksplit(input0, 20, 15, field, nr)
159  if(nr .eq. 8) then
160  read(input0(1:klena(input0)), *)
161  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0
162  maxsites = ansites0
163  emin = 500.d-6
164  newfmt = 0
165  elseif(nr .eq. 9 ) then
166  read(input0(1:klena(input0)), *)
167  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0,
168  * maxsites
169  emin = 500.d-6
170  newfmt = 1
171  elseif(nr .eq. 13) then
172  read(input0(1:klena(input0)), *)
173  * evno0, e0,nn, cosz, emin, nrbina, nfaia, ansites0,
174  * maxsites, limit
175  newfmt = 2
176  else
177  write(0,*) ' header fmt strage '
178  write(0,*) input0
179  stop 01234
180  endif
181 
182 ! write(0,*) input0
183 
184  if(nrbina .ne. nrbin .or. nfaia .ne. nfai) then
185  write(0,*)' nrbina=',nrbina, 'or nfaia=',nfaia,
186  * ' differ from the def. in this prog'
187  stop 5555
188  endif
189  if(ansites0 .gt. ndepth) then
190  write(0,*) ' too many depths'
191  stop 6666
192  endif
193 ! ********
194  do i = 1, ansites0
195  do j = 1, 4
196  do k = 1, nfai
197  read(fnonrfai, '(a, f7.1, 4i4)' )
198  * id, recdep(i), recl0(i), i0, j0, k0
199  indivdep(i)=recl0(i)
200  packeddepidx(recl0(i))=i ! original dep index to packed indx
201 
202 !
203 ! when the above is written
204 ! l = indivdep(i)
205 ! write(fnonrfai, '("rec",f7.1, 4i4)' )
206 ! * ASDepthList(l)*0.1, l, i, j, k
207 
208 
209  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
210  write(0,*) ' intdep, i0,j0,k0=',
211  * recdep(i), i0, j0, k0, ' strange'
212  stop 8888
213  endif
214  if( id .ne. "rec") then
215  write(0,*) 'id=',id, ' strange'
216  stop 5678
217  endif
218 !
219  read(fnonrfai, *)
220  * ( nrfairec0(l,k,j,i), l=1,nrbin )
221  enddo
222  enddo
223  enddo
224 ! ************
225  do i = 1, maxsites
226  do j = 1, 4
227  do k = 1, nfai
228  read(fnonrfai, '(a,f7.1, 4i4)' )
229  * id, intdep(i), l0(i), i0, j0, k0
230  indivdep(i)=l0(i)
231  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
232  write(0,*) ' intdep, i0,j0,k0=',
233  * intdep(i), i0, j0, k0, ' strange'
234  stop 9876
235  endif
236  if( id .ne. "all") then
237  write(0,*) 'id=',id,' strange'
238  stop 9999
239  endif
240  read(fnonrfai, *)
241  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
242  enddo
243  enddo
244  enddo
245 ! dErfai
246  do i = 1, maxsites
247  do k = 1, nfai
248  read(fnonrfai, '(a,f7.1, 3i4)' )
249  * id5, intdep(i), l0(i), i0, k0
250  indivdep(i)=l0(i)
251  if(i0 .ne. i .or. k .ne. k0) then
252  write(0,*) ' intdep, i0,k0=',
253  * intdep(i), i0, k0, ' strange'
254  stop 98765
255  endif
256  if( id5 .ne. "dE/dx") then
257  write(0,*) 'id=',id5,' strange'
258  stop 9999
259  endif
260  read(fnonrfai, *)
261  * ( derfai(l,k,i), l=1,nrbin )
262  if(seelowde) read(fnonrfai, *)
263  * ( derfai2(l,k,i), l=1,nrbin )
264  enddo
265  enddo
266  else
267 !
268  write(0,*) ' all nrfai data has been read'
269  limit(1) = reallimitg
270  limit(2) = reallimite
271  limit(3) = reallimitmu
272  limit(4) = reallimith
273  if(newfmt .eq. 0) then
274  write(0, '(i2, 1pE11.3,i3, 0pf7.1, 1pE11.3,0p, 3i4)')
275  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0
276  elseif( newfmt .eq. 1 ) then
277  write(0, '(i2, 1pE11.3,i3, 0pf7.1, 1pE11.3, 0p,4i4)')
278  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0,
279  * maxsites
280  else
281  write(0, '(i2, 1pE11.3,i3, 0pf7.1,1pE11.3,0p,4i4,4f8.9)')
282  * evno0, e0,nn, cosz, emin, nrbina, nfaia, ansites0,
283  * maxsites, limit
284  endif
285  endif
286  enddo
287  1000 continue
288  input0= ' '
289 !************* reset limit
290 ! limit=min(4000.0, limit)
291  limit(1) = reallimitg
292  limit(2) = reallimite
293  limit(3) = reallimitmu
294  limit(4) = reallimith
295 
296  if(newfmt .ne. 2) then
297  write(0,*) ' limit=',limit(1), ' is being used'
298  else
299  write(0,*) ' limit=',limit, ' ares being used'
300  endif
301 !************
302  do i = 1, ansites0
303  do j = 1, 4
304  do k = 1,nfai
305  do l = 1, nrbin
306  if(nrfairec0(l, k, j, i) .gt. limit(j)) then
307 ! accept with this prob.
308  pnrfairec(l, k, j, i)=
309  * limit(j)/nrfairec0(l, k, j, i)
310  else
311  pnrfairec(l, k, j, i)=1.0
312  endif
313  enddo
314 !///////
315  write(0, '(f7.1, 4i4,a)' )
316  * recdep(i), indivdep(i), i, j, k, ' prob'
317  write(0, '(1p10E11.3)')
318  * (pnrfairec(l,k,j,i), l=1, nrbin)
319 !c/////
320  enddo
321  enddo
322  enddo
323 ! -------------
324 ! main input data; header Obsolete
325 ! read(fnodat) EventNo, code, Et, wx, wy, wz
326 ! if( header(1:lengh) .eq. "yes") then
327 ! write(*,'("i ", i3, i4, g13.4,3f11.7)')
328 ! * EventNo, code, Et, wx, wy, wz
329 ! endif
330 !
331 ! open files with -Lxx.dat (xx=recl0)
332 !
333  do i = 1, ansites0
334  if(recl0(i) .gt. 0) then
335  write(lfile(recl0(i)),'(a,"/", a, "-L",i2.2,".dat")')
336  * basedir(1:lengbasedir), execid(1:lengexecid), recl0(i)
337  open(100+recl0(i), file=lfile(recl0(i)), access="append" ,
338  * form="formatted", iostat=icon)
339  if(icon .ne. 0) then
340  write(0,*) lfile(recl0(i)), " could not be opend"
341  stop 9999
342  endif
343  endif
344  enddo
345 !
346 !
347  nrec= 0
348  do while(.true.)
349  read(fnodat,end=100) bufc, (buf(i), i=1, bufc)
350 !/////////////
351 ! write(0,*) ' bufc=',bufc
352 !///////////
353  do i = 1, bufc
354  nrec= nrec+1
355  ldep = buf(i).ldep
356  depidx = packeddepidx(ldep)
357  faiidx= buf(i).faiidx
358  ridx = buf(i).ridx
359  codex = min(buf(i).code, 4)
360  wgt = buf(i).wgt
361  prob = pnrfairec(ridx, faiidx, codex, depidx)
362  if( prob .gt. 1.) then
363  wwgt = wgt
364  accept = .true.
365  else
366  prob = prob * wgt
367  if(prob .gt. 1.) then
368  accept = .true.
369  wwgt = prob
370  else
371  call rndc(u)
372  if(u .lt. prob) then
373  accept = .true.
374  wwgt= 1.
375  else
376  accept = .false.
377  endif
378  endif
379  endif
380  if(accept) then
381  rnrfairec(ridx, faiidx, codex, depidx)=
382  * rnrfairec(ridx, faiidx, codex, depidx) + wwgt
383 ! write(*,
384  write(100+buf(i).ldep,
385  * '(6i3, 1pE11.3, 0p,f6.1,1p2E11.3,0p, 2f8.4,f10.6,1pE11.3)')
386  * buf(i).ldep, buf(i).code, buf(i).subcode,
387  * buf(i).charge, buf(i).ridx, buf(i).faiidx,
388  * buf(i).rinmu, buf(i).fai, buf(i).ek,
389  * buf(i).t, buf(i).wx, buf(i).wy, buf(i).wz,
390  * wwgt
391  endif
392  enddo
393  enddo
394  100 continue
395 
396  write( fnonrfai2,
397  * '(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3,0p,4i4,4f8.0)' )
398  * evno0, e0, nn, cosz, emin, nrbin, nfai, ansites0,
399  * maxsites, limit
400 
401  do i = 1, ansites0
402  do j = 1, 4
403  do k = 1, nfai
404  l = recl0(i)
405  write(fnonrfai2, '("rec",f7.1, 4i4)' )
406  * recdep(i), l, i, j, k
407  write(fnonrfai2, '(1p10E11.3)')
408  * ( rnrfairec(l,k,j,i), l=1,nrbin )
409  enddo
410  enddo
411  enddo
412  do i = 1, maxsites
413  do j = 1, 4
414  do k = 1, nfai
415 ! l = indivdep(i)
416  l = l0(i)
417  write(fnonrfai2, '("all",f7.1, 4i4)' )
418  * intdep(i), l, i, j, k
419  write(fnonrfai2, '(1p10E11.3)')
420  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
421  enddo
422  enddo
423  enddo
424 ! dErfai
425  do i = 1, maxsites
426  do k = 1, nfai
427 ! l = indivdep(i)
428  l = l0(i)
429  write(fnonrfai2, '("dE/dx",f7.1, 3i4)' )
430  * intdep(i), l, i, k
431  write(fnonrfai2, '(1p10E11.3)')
432  * ( derfai(l,k,i), l=1,nrbin )
433  if(seelowde) write(fnonrfai2, '(1p10E11.3)')
434  * ( derfai2(l,k,i), l=1,nrbin )
435  enddo
436  enddo
437 
438  write(0,*) 'end of run'
439  write(fnonrfai2, *)
440  stop
441  500 continue
442  write(0,*) ' input error at record =', nrec
443  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