COSMOS v7.655  COSMOSv7655
(AirShowerMC)
add2nrfai.f
Go to the documentation of this file.
1  implicit none
2 ! Suppose we have n nrfai data
3 ! h1, h2, h3, ....hn
4 ! 1) cp h1 h0
5 ! 2) h0 + h2 --> h; mv h h0
6 ! 3) h0 + h3 --> h; mv h h0
7 ! ..
8 ! n) h0 + hn --> h; mv h h0
9 !
10 ! This program add two nrfai data; h0 + hx--> h
11 ! environmental variable
12 ! file h0: NRFAIFILE0
13 ! file hx: NRFAIFILEX
14 ! file h: NRFAIFILE
15 ! "all" is should not touched: GETALL
16 !
17 !
18 #include "Zobs.h"
19 #include "../FleshHist/Zprivate0.h"
20 
21  integer ndepth
22  parameter(ndepth= nsites)
23  real nrfaiRec0(nrbin, nfai, 4, ndepth)
24  real nrfaiAll0(nrbin, nfai, 4, ndepth)
25 
26  real nrfaiRecx(nrbin, nfai, 4, ndepth)
27  real nrfaiAllx(nrbin, nfai, 4, ndepth)
28 
29  real nrfaiRect(nrbin, nfai, 4, ndepth)
30  real nrfaiAllt(nrbin, nfai, 4, ndepth)
31 
32  real dErfai0(nrbin, nfai, ndepth)
33  real dErfai02(nrbin, nfai, ndepth)
34  real dErfaix(nrbin, nfai, ndepth)
35  real dErfaix2(nrbin, nfai, ndepth)
36  real dErfait(nrbin, nfai, ndepth)
37  real dErfait2(nrbin, nfai, ndepth)
38 
39  integer EvNo0, EvNox
40  integer fn0, fnx, fnt
41  integer kgetenv2
42  integer klena
43  real intdep(ndepth), recdep(ndepth)
44  real E0, cosz, limit(4)
45  integer NN
46  integer nrbina, nfaia, ansites0, ansitesx
47  integer maxsites
48  integer newfmt
49  integer leng, i, j, k, l, nr
50  integer i0, j0, k0, reci0, recj0, reck0
51  integer l0(ndepth), recl0(ndepth)
52  integer icon0, iconx, icont
53  character*128 nrfai0, nrfaix, nrfait, flesher
54  character*128 input0, inputx, inputt
55  character*20 field(15)
56  character*5 getall
57  real*4 Emin
58  logical SeeLowdE
59 
60  fn0 = 2
61  fnx = 3
62  fnt = 4
63  leng = kgetenv2("GETALL", getall)
64  getall=getall(1:leng)
65 
66  write(0,*) ' getall =', getall
67 
68  leng = kgetenv2("NRFAIFILE0", nrfai0)
69  call copenfw2(fn0, nrfai0, 1, icon0)
70  if(icon0 .ne. 1) then
71  write(0,*) nrfai0(1:leng)
72  if( icon0 .eq. 0) then
73  write(0,*) 'not exists'
74  else
75  write(0,*) ' cannot be opened '
76  endif
77  write(0,*) ' icon=',icon0
78  stop 9999
79  else
80  write(0,*) nrfai0(1:leng), ' opened'
81  endif
82  leng = kgetenv2("NRFAIFILEX", nrfaix)
83  call copenfw2(fnx, nrfaix, 1, iconx)
84  if(iconx .ne. 1) then
85  write(0,*) nrfaix(1:leng)
86  if( iconx .eq. 0) then
87  write(0,*) 'not exists'
88  else
89  write(0,*) ' cannot be opened '
90  endif
91  write(0,*) ' icon=',iconx
92  stop 9999
93  else
94  write(0,*) nrfaix(1:leng), ' opened'
95  endif
96 
97 
98  leng = kgetenv2("NRFAIFILET", nrfait)
99  call copenfw2(fnt, nrfait, 1, icont)
100  if(icont .ne. 0) then
101  write(0,*) nrfait(1:leng)
102  write(0,*) ' cannot be opened '
103  write(0,*) ' icon=',icont
104  stop 9999
105  else
106  write(0,*) nrfait(1:leng), ' opened'
107  endif
108  leng = kgetenv2("FLESHDIR", flesher)
109 
110  leng = kgetenv2("SeeLowdE", input0)
111  if(leng .le. 0) then
112  write(0,*) "SeeLowdE not given"
113  stop
114  endif
115  seelowde = input0(1:leng) .eq. "yes"
116 
117 
118 
119 
120 
121 
122  do while(.true.)
123  input0 = ' '
124  read( fn0, '(a)', end=1000 ) input0
125  inputx = ' '
126  read( fnx, '(a)' ) inputx
127 
128 
129  if(input0 .ne. " ") then
130  call ksplit(input0, 20, 15, field, nr)
131 !//////////////
132 ! write(0,*) ' input0=', input0
133 ! write(0,*) ' nr=',nr
134 !///////////
135  if(nr .eq. 8) then
136  read(input0(1:klena(input0)), *)
137  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0
138  maxsites = ansites0
139  emin=500.d-6
140  newfmt = 0
141  elseif(nr .eq. 9 ) then
142  read(input0(1:klena(input0)), *)
143  * evno0, e0,nn, cosz, limit(1), nrbina, nfaia, ansites0,
144  * maxsites
145  emin=500.d-6
146  newfmt = 1
147  elseif( nr .eq. 13) then
148  read(input0(1:klena(input0)), *)
149  * evno0, e0,nn, cosz, emin, nrbina, nfaia, ansites0,
150  * maxsites, limit
151  newfmt = 2
152  else
153  write(0,*) ' header fmt strage '
154  write(0,*) input0
155  stop 01234
156  endif
157  if(nrbina .ne. nrbin .or. nfaia .ne. nfai) then
158  write(0,*)' nrbina=',nrbina, 'or nfaia=',nfaia,
159  * ' differ from the def. in this prog'
160  stop 5555
161  endif
162  if(newfmt .eq. 0) then
163  read(inputx(1:klena(inputx)), *)
164  * evnox, e0,nn, cosz, limit(1),
165  * nrbina, nfaia, ansitesx
166  elseif(newfmt .eq. 1) then
167  read(inputx(1:klena(inputx)), *)
168  * evnox, e0,nn, cosz, limit(1), nrbina, nfaia,
169  * ansitesx, maxsites
170  elseif(newfmt .eq. 2) then
171  read(inputx(1:klena(inputx)), *)
172  * evnox, e0,nn, cosz, emin, nrbina, nfaia,
173  * ansitesx, maxsites, limit
174  endif
175 
176  if(nrbina .ne. nrbin .or. nfaia .ne. nfai) then
177  write(0,*)' nrbina=',nrbina, 'or nfaia=',nfaia,
178  * ' differ from the def. in this prog'
179  stop 6666
180  endif
181  if(ansites0 .ne. ansitesx) then
182  write(0,*) ' ansites0=', ansites0,
183  * ' ansitesx=', ansitesx, ' diff '
184  stop 7777
185  endif
186 
187 ! ********
188  do i = 1, ansites0
189  do j = 1, 4
190  do k= 1, nfai
191  read(fn0, '(3x, f7.1, 4i4)' )
192  * recdep(i), recl0(i), reci0, recj0, reck0
193 
194  if(reci0 .ne. i .or. j .ne. recj0 .or.
195  * k .ne. reck0) then
196  write(0,*) ' recdep, reci0,recj0,reck0=',
197  * recdep(i), reci0, recj0, reck0, ' strange'
198  stop 8888
199  endif
200  read(fn0, *)
201  * ( nrfairec0(l,k,j,i), l=1,nrbin )
202 
203  enddo
204  enddo
205  enddo
206 
207 ! ************
208  do i = 1, maxsites
209  do j = 1, 4
210  do k = 1, nfai
211  read(fn0, '(3x,f7.1, 4i4)' )
212  * intdep(i), l0(i), i0, j0, k0
213  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
214  write(0,*) ' intdep, i0,j0,k0=',
215  * intdep(i), i0, j0, k0, ' strange'
216  stop 9876
217  endif
218  read(fn0, *)
219  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
220  enddo
221  enddo
222  enddo
223 
224 
225  do i = 1, maxsites
226  do k = 1, nfai
227  read(fn0, '(5x,f7.1, 3i4)' )
228  * intdep(i), l0(i), i0, k0
229  if(i0 .ne. i .or. k .ne. k0) then
230  write(0,*) ' intdep, i0,k0=',
231  * intdep(i), i0, k0, ' strange'
232  stop 9875
233  endif
234  read(fn0, *)
235  * ( derfai0(l,k,i), l=1,nrbin )
236  if(seelowde) read(fn0, *)
237  * ( derfai02(l,k,i), l=1,nrbin )
238  enddo
239  enddo
240 
241 
242 ! ************
243 ! ************
244  do i = 1, ansites0
245  do j = 1, 4
246  do k= 1, nfai
247  read(fnx, '(3x, f7.1, 4i4)' )
248  * recdep(i), recl0(i), reci0, recj0, reck0
249  if(reci0 .ne. i .or. j .ne. recj0 .or.
250  * k .ne. reck0) then
251  write(0,*) ' recdep, reci0,recj0,reck0=',
252  * intdep(i), reci0, recj0, reck0, ' strange'
253  stop 9999
254  endif
255  read(fnx, *)
256  * ( nrfairecx(l,k,j,i), l=1,nrbin )
257  enddo
258  enddo
259  enddo
260 
261 ! ************
262  do i = 1, maxsites
263  do j = 1, 4
264  do k = 1, nfaia
265  read(fnx, '(3x,f7.1, 4i4)' )
266  * intdep(i), l0(i), i0, j0, k0
267  if(i0 .ne. i .or. j .ne. j0 .or. k .ne. k0) then
268  write(0,*) ' intdep, i0,j0,k0=',
269  * intdep(i), i0, j0, k0, ' strange'
270  stop 9877
271  endif
272  read(fnx, *)
273  * ( nrfaiallx(l,k,j,i), l=1,nrbin )
274  enddo
275  enddo
276  enddo
277 
278  do i = 1, maxsites
279  do k = 1, nfaia
280  read(fnx, '(5x,f7.1, 3i4)' )
281  * intdep(i), l0(i), i0, k0
282  if(i0 .ne. i .or. k .ne. k0) then
283  write(0,*) ' intdep, i0,k0=',
284  * intdep(i), i0, k0, ' strange'
285  stop 98778
286  endif
287  read(fnx, *)
288  * ( derfaix(l,k,i), l=1,nrbin )
289  if(seelowde) read(fnx, *)
290  * ( derfaix2(l,k,i), l=1,nrbin )
291  enddo
292  enddo
293 
294 ! ************
295  else
296 
297 ! 1 event read
298  if(inputx .ne. " ") then
299  write(0,*) ' event differ', inputx
300  stop 1111
301  endif
302  do i = 1, ansites0
303  do j = 1, 4
304  do k= 1, nfai
305  do l=1, nrbin
306  nrfairect(l,k,j,i) =
307  * nrfairec0(l,k,j,i) +
308  * nrfairecx(l,k,j,i)
309  enddo
310  enddo
311  enddo
312  enddo
313  if(getall .eq. 'yes') then
314  do i = 1, maxsites
315  do j = 1, 4
316  do k= 1, nfai
317  do l=1, nrbin
318  nrfaiallt(l,k,j,i) =
319  * nrfaiall0(l,k,j,i) +
320  * nrfaiallx(l,k,j,i)
321  enddo
322  enddo
323  enddo
324  enddo
325  endif
326  if(getall .eq. 'yes') then
327  do i = 1, maxsites
328  do k= 1, nfai
329  do l=1, nrbin
330  derfait(l,k,i) =
331  * derfai0(l,k,i) +
332  * derfaix(l,k,i)
333  if(seelowde) derfait2(l,k,i) =
334  * derfai02(l,k,i) +
335  * derfaix2(l,k,i)
336 
337  enddo
338  enddo
339  enddo
340  endif
341  if( newfmt .eq. 0 ) then
342  write(fnt,
343  * '(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3, 0p, 3i4)')
344  * evnox, e0, nn, cosz, limit(1),
345  * nrbina, nfaia, ansites0
346  elseif( newfmt .eq. 1) then
347  write(fnt,
348  * '(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3,0p, 4i4)')
349  * evnox, e0, nn, cosz, limit(1),
350  * nrbina, nfaia, ansites0, maxsites
351  elseif( newfmt .eq. 2) then
352  write(fnt,
353  * '(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3,0p, 4i4,4f10.0)')
354  * evnox, e0, nn, cosz, emin,
355  * nrbina, nfaia, ansites0, maxsites, limit
356  endif
357 
358  do i = 1, ansites0
359  do j = 1, 4
360  do k = 1, nfaia
361  write(fnt, '("rec",f7.1, 4i4)' )
362  * recdep(i), recl0(i), i, j, k
363  write(fnt, '(1p10E11.3)')
364  * ( nrfairect(l,k,j,i), l=1,nrbin )
365  enddo
366  enddo
367  enddo
368 
369  do i = 1, maxsites
370  do j = 1, 4
371  do k = 1, nfaia
372  write(fnt, '("all",f7.1, 4i4)' )
373  * intdep(i), l0(i), i, j, k
374  if(getall .eq. "yes") then
375  write(fnt, '(1p10E11.3)')
376  * ( nrfaiallt(l,k,j,i), l=1,nrbin )
377  else
378  write(fnt, '(1p10E11.3)')
379  * ( nrfaiall0(l,k,j,i), l=1,nrbin )
380  endif
381  enddo
382  enddo
383  enddo
384 
385  do i = 1, maxsites
386  do k = 1, nfaia
387  write(fnt, '("dE/dx",f7.1, 3i4)' )
388  * intdep(i), l0(i), i, k
389  if(getall .eq. "yes") then
390  write(fnt, '(1p10E11.3)')
391  * ( derfait(l,k,i), l=1,nrbin )
392  if(seelowde) write(fnt, '(1p10E11.3)')
393  * ( derfait2(l,k,i), l=1,nrbin )
394  else
395  write(fnt, '(1p10E11.3)')
396  * ( derfai0(l,k,i), l=1,nrbin )
397  if(seelowde) write(fnt, '(1p10E11.3)')
398  * ( derfai02(l,k,i), l=1,nrbin )
399  endif
400  enddo
401  enddo
402 
403  write(fnt, *)
404  endif
405  enddo
406  1000 continue
407 
408  write(0,*) 'all data in the event processed '
409  end
410 
411 
412 
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
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
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
klena
integer function klena(cha)
Definition: klena.f:20
nfai
! timing nfai
Definition: Zprivate2.h:12