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