19 #include "Zincidentv.h" 21 #include "Zprivate1.h" 22 #include "Zprivate2.h" 23 #include "Zprivate3.h" 25 type(
track):: incident
26 type(
coord):: AngleAtObs
34 save rspec, lossrspec, arspec, respec
35 save rzspec, zfspec, rtspec1, rtspec2, retspec1, retspec2
36 save rezspec, rzfspec, rfspec, efspec, refspec
52 integer i, j, k, icon, ir, l
54 integer iij, code, codex
57 real*8 r, Eloss, rinmu, cosang
58 real*8 dedt, rho, dist, disto, dedxmu
63 real*8 wx, wy, wz, temp
67 integer ldep, ridx, faiidx
75 data ptcl/
"Photons",
"Electrons",
"Muons",
"hadrons"/
77 data ptcl2/
"Electrons",
"Muons",
"All"/
80 real E0, cosz, limit, limite
85 character*96 evid(nsites)
86 real*8 cog, cog2, sumne, obstimes, Savederg(5)
95 #include "interface2.f" 107 read(input(1:lengenv),*) ncpu
109 read(input(1:lengenv),*) mcpu
111 read(input(1:lengenv),*) margin
113 enhance = enhance/mcpu
114 limite = limit * enhance
115 write(0,*)
' *** ncpu=',ncpu,
116 *
' mcpu=',mcpu,
' enhance factor=',enhance
117 write(0,*)
' margin=',margin
129 if( indivdep(i) .eq. 0)
exit 131 w2il( indivdep(i) ) = i
144 rminm = rmin/10.**(
bin/2.0
d0)
167 *
"Lateral Dist. of "//
ptcl(j),
168 *
"lat",
"ptcls", .
true., power(j),
"r",
"m.u")
183 *
"dE/dx lateral dist. of "//ptcl2(j),
184 *
"dEdxLat",
"GeV/(g/cm^2)", .
true., power2(j),
194 call kwhisti2(arspec(j, i),
195 * 0., 30.0, 12,
b'00010',
197 call kwhistai2(arspec(j,i),
198 *
"Lateral dits. of "//
ptcl(j)//
" with given Fai bin",
199 *
"ar",
"ptcls", .
true., power(j),
200 *
"azimuth",
"deg",
"r",
"m.u")
210 call kwhisti2( respec(j, i),
211 * 0.01, 0.2, 20,
b'00011',
212 * 500.
e-6, 0.1, 50,
b'00001')
213 call kwhiststep2(respec(j, i), 2)
214 call kwhistai2( respec(j, i),
215 *
"Energy Spec. of "//
ptcl(j)//
" at diff. r",
216 *
"re",
"ptcls", .
true., 1.,
217 *
"r",
"m.u",
"E",
"GeV")
220 call kwhisti2( respec(3, i),
221 * 0.01, 0.2, 20 ,
b'00011',
222 * 0.031627, 0.1, 38,
b'00011' )
223 call kwhiststep2(respec(3, i), 2)
224 call kwhistai2( respec(3, i),
225 *
"Energy Spec. of mu at diff. r",
226 *
"re",
"ptcls", .
true., 0.,
227 *
"r",
"m.u",
"E",
"GeV")
235 call kwhisti2( rzspec(j, i),
236 * 0.01, 0.2, 20,
b'00011',
237 * 0., 1.0, 20,
b'10000')
238 call kwhiststep2(rzspec(j, i), 2)
239 call kwhistai2( rzspec(j, i),
240 *
"Zenith angle dist. of "//
ptcl(j)//
" at diff. r",
241 *
"rz",
"ptcls", .
true., 0.,
242 *
"r",
"m.u",
"cosz",
" ")
244 call kwhisti2(rzspec(3, i),
245 * 0.01, 0.2, 20,
b'00011',
246 * 0., 1.0, 20,
b'10000' )
247 call kwhiststep2(rzspec(3, i), 2)
248 call kwhistai2( rzspec(3, i),
249 *
"Zenith angle dist. of m with diff. r",
250 *
"rz",
"ptcls", .
true., 0.,
251 *
"r",
"m.u",
"cosz",
" ")
259 call kwhisti2( zfspec(j, i),
260 * 0.0, 1.0, 10,
b'10000',
261 * -1.0, 1.0, 50,
b'10000')
262 call kwhiststep2(zfspec(j, i), 2)
263 call kwhistai2( zfspec(j, i),
264 *
"f=(wx,wy)*(x,y) spectrum of "//
ptcl(j)//
265 *
" with diff. cosz",
266 *
"zf",
"ptcls", .
true., 0.,
267 *
"cosz",
" ",
"f",
" ")
276 call kwhisti2( rfspec(j, i),
277 * 0.01, 0.2, 20,
b'00011',
278 * -1.0, 1.0, 50,
b'10000')
279 call kwhiststep2(rfspec(j, i), 2)
280 call kwhistai2( rfspec(j, i),
281 *
"f spectrum of "//
ptcl(j)//
" with diff. r",
282 *
"rf",
"ptcls", .
true., 0.,
283 *
"r",
"m.u",
"f",
" ")
292 call kwhisti2( efspec(j, i),
293 * 500.
e-6, 0.2, 20,
b'00001',
294 * -1.0, 1.0, 50,
b'10000')
295 call kwhiststep2(efspec(j, i), 2)
296 call kwhistai2( efspec(j, i),
297 *
"f spectrum of "//
ptcl(j)//
" with diff. E",
298 *
"ef",
"ptcls",.
true., 0.,
299 *
"E",
"GeV",
"f",
" ")
307 call kwhisti2( rtspec1(j, i),
308 * 1.0, 0.2, 10,
b'00011',
309 * 10., 0.1, 55,
b'00011' )
310 call kwhiststep2(rtspec1(j, i), 2)
311 call kwhistai2( rtspec1(j, i),
312 *
"Arrival time dist. of "//
ptcl(j)//
" at diff. r",
313 *
"rt",
"ptcls", .
true., 0.,
314 *
"r",
"m.u",
"time",
"ns")
316 call kwhisti2( rtspec2(j, i),
317 * 0., 1.0, 10,
b'10000',
318 * 0., 0.25, 300,
b'00000' )
319 call kwhiststep2(rtspec2(j, i), 2)
320 call kwhistai2( rtspec2(j, i),
321 *
"Arrival time dist. of "//
ptcl(j)//
" at diff. r",
322 *
"rt",
"ptcls", .
true., 0.,
323 *
"r",
"m.u",
"time",
"ns")
332 call kwhisti3(retspec1(j, i),
333 * 1.0, 0.2, 10,
b'00011',
334 * 500.
e-6, 0.25, 15,
b'01001',
335 * 10.0, 0.1, 55,
b'00011' )
336 call kwhiststep3(retspec1(j, i), 2, 2)
337 call kwhistai3(retspec1(j, i),
338 *
"Arrival time dist. of "//
ptcl(j)//
" with diff. r&E",
339 *
"ret",
"ptcls", .
true., 0.,
340 *
"r",
"m.u",
"E",
"GeV",
"Time",
"ns")
342 call kwhisti3(retspec2(j, i),
343 * 0.0, 1.0, 10,
b'10000',
344 * 500.
e-6, 0.25, 12,
b'01001',
345 * 0.0, 0.25, 200,
b'00000' )
346 call kwhiststep3(retspec2(j, i), 2,2)
347 call kwhistai3(retspec2(j, i),
348 *
"Arrival time dist. of "//
ptcl(j)//
" with diff. r&E",
349 *
"ret2",
"ptcls", .
true., 0.,
350 *
"r",
"m.u",
"E",
"GeV",
"Time",
"ns")
353 call kwhisti3(retspec1(3, i),
354 * 1.0, 0.2, 10,
b'00011',
355 * 100.
e-3, 0.25, 10,
b'01011',
356 * 10.0, 0.1, 72,
b'00011' )
357 call kwhiststep3(retspec1(3, i), 2, 2)
358 call kwhistai3(retspec1(3, i),
359 *
"Arrival time dist. of "//
ptcl(3)//
" with diff. r&E",
360 *
"ret",
"ptcls", .
true., 0.,
361 *
"r",
"m.u",
"E",
"GeV",
"Time",
"ns")
364 call kwhisti3(retspec2(3, i),
365 * 0.0, 1.0, 10,
b'10000',
366 * 100.
e-3, 0.25, 10,
b'01011',
367 * 0.0, 0.25, 200,
b'00000' )
368 call kwhiststep3(retspec2(3, i), 2,2)
369 call kwhistai3(retspec2(3, i),
370 *
"Arrival time dist. of "//
ptcl(3)//
" with diff. r&E",
371 *
"ret2",
"ptcls", .
true., 0.,
372 *
"r",
"m.u",
"E",
"GeV",
"Time",
"ns")
380 call kwhisti3( rezspec(j, i),
381 * 0.1, 0.2, 15,
b'00011',
382 * 500.
e-6, 0.25, 14,
b'01001',
383 * 0.0, 1.0, 20,
b'10000')
384 call kwhiststep3(rezspec(j, i), 3,2)
385 call kwhistai3(rezspec(j, i),
386 *
"cos zenith dist. of "//
ptcl(j)//
" with diff. r&E",
387 *
"rez",
"ptcls", .
true., 0.,
388 *
"r",
"m.u",
"E",
"GeV",
"cosz",
" ")
390 call kwhisti3(rezspec(3, i),
391 * 0.1, 0.2, 15,
b'00011',
392 * 100.
e-3, 0.25, 10,
b'01011',
393 * 0.0, 1.0, 20,
b'10000')
394 call kwhiststep3(rezspec(3, i), 3,2)
395 call kwhistai3(rezspec(3, i),
396 *
"cos zenith dist. of "//
ptcl(3)//
" with diff. r&E",
397 *
"rez",
"ptcls", .
true., 0.,
398 *
"r",
"m.u",
"E",
"GeV",
"cosz",
" ")
407 call kwhisti3(rzfspec(j, i),
408 * 0.1, 0.2, 15,
b'00011',
409 * 0.0, 1.0, 10,
b'10000',
410 * -1.0, 1.0, 20,
b'10000')
411 call kwhiststep3(rzfspec(j, i), 3,2)
412 call kwhistai3(rzfspec(j, i),
413 *
"f spectrum of "//
ptcl(j)//
" with diff r&cosz",
414 *
"rzf",
"ptcls", .
true., 0.,
415 *
"r",
"m.u",
"cosz",
" ",
"f",
" ")
423 call kwhisti3( refspec(j, i),
424 * 0.1, 0.2, 15,
b'00011',
425 * 500.
e-6, 0.25, 16,
b'01001',
426 * -1.0, 1.0, 20,
b'10000')
427 call kwhiststep3(refspec(j, i), 3,3)
428 call kwhistai3(refspec(j, i),
429 *
"f spectrum of "//
ptcl(j)//
" with diff. r&E",
430 *
"ref",
"ptcls", .
true., 0.,
431 *
"r",
"m.u",
"E",
"GeV",
"f",
" ")
433 call kwhisti3(refspec(3, i),
434 * 0.1, 0.2, 15,
b'00011',
435 * 100.
e-3, 0.25, 10,
b'01001',
436 * -1.0, 1.0, 20,
b'10000')
437 call kwhiststep3(refspec(3, i), 3,2)
438 call kwhistai3(refspec(3, i),
439 *
"f spectrum of "//
ptcl(3)//
" with diff. r&E",
440 *
"ref",
"ptcls", .
true., 0.,
441 *
"r",
"m.u",
"E",
"GeV",
"f",
" ")
463 if(inci.
p.code .eq. 9)
then 465 elseif(inci.
p.code .eq. 1)
then 472 firstz= zfirst.pos.
depth*0.1
473 write(0,
'(a,1pE11.3,a,E11.3,a,E11.3,a)')
474 *
' E0=',e0,
' cosz=',cosz,
' firstz=',
477 call howmuch(limite, e0, nn, cosz)
480 do i = 1, noofassites
508 call kwhistc2(arspec(j,i))
516 call kwhistc2(respec(j,i))
525 call kwhistc2(rzspec(j,i))
533 call kwhistc2(rfspec(j,i))
543 call kwhistc2(efspec(j,i))
554 call kwhistc2(rtspec1(j,i))
555 call kwhistc2(rtspec2(j,i))
564 call kwhistc3(retspec1(j,i))
565 call kwhistc3(retspec2(j,i))
575 call kwhistc3(rezspec(j,i))
585 call kwhistc3(rzfspec(j,i))
594 call kwhistc3(refspec(j,i))
603 nrfairec(l, k, j, i) = 0.
604 nrfaiall(l, k, j, i) = 0.
615 entry xobs(atrack, id)
621 obstimes = obstimes + 1.
d0 622 if(mod(obstimes, 500000.
d0) .eq. 0. )
then 624 do i = 1, min(4,stack_pos)
625 if(stack(i).
p.fm.
p(4) .ne. savederg(i))
then 626 savederg(i)=stack(i).
p.fm.
p(4)
633 write(0, *)
' obstimes=', obstimes,
' ptclE=',atrack.
p.fm.
p(4)
634 do i = 1, min(4,stack_pos)
635 write(0,*)
' stack tops=', stack(i).
p.fm.
p(4)
640 if(id .eq. 2 .and. code .le. 6)
then 642 wz = atrack.vec.w.r(3)
648 ng(ldep) =
ng(ldep) + enhance
649 elseif(code .eq.
kelec)
then 650 ne(ldep) =
ne(ldep) + enhance
651 elseif(code .eq.
kmuon)
then 652 nmu(ldep) =
nmu(ldep) + enhance
654 nhad(ldep)=nhad(ldep) + enhance
657 r = sqrt( atrack.pos.
xyz.
x**2 +
658 * atrack.pos.
xyz.
y**2 )
660 ek = atrack.
p.fm.
p(4) -atrack.
p.
mass 663 if(atrack.
p.
charge .ne. 0 )
then 664 rho = cvh2den(atrack.pos.
height)
671 if(abs(wz) .gt. 1.
d-2)
then 681 if(atrack.
p.code .eq.
kmuon )
then 683 call cmudedx(muni, mubr, mupr, atrack.
p.fm.
p(4),
690 call cqelossrate(dedt)
693 eloss = dedt*dist*enhance
695 sumeloss(ldep)=sumeloss(ldep) + eloss
701 if(j .eq. 0)
goto 100
703 molu =obssites(ldep).mu
706 if(rinmu .gt. rminm)
then 707 ridx= log10(rinmu/rminm)/
bin + 1
712 fai=atan2(atrack.pos.
xyz.
y, atrack.pos.
xyz.
x)*
713 * 57.29577951308230
d0 714 fai= mod(fai + 360.
d0, 360.
d0)
715 if(fai .gt. (360.
d0-dfai/2.0
d0)) fai= fai-360.
d0 716 faiidx=(fai+dfai/2.0
d0) /dfai + 1
718 if(ridx .gt. 0 .and. ridx .le.
nrbin )
then 720 accept = u .le. recprob(ridx, codex, j)
722 nrfairec(ridx, faiidx, codex, j) =
723 * nrfairec(ridx, faiidx, codex, j) + 1.0
725 nrfaiall(ridx, faiidx, codex, j) =
726 * nrfaiall(ridx, faiidx, codex, j) + 1.0
727 derfai(ridx, faiidx, j) =
728 * derfai(ridx, faiidx, j) + eloss
734 if(bufc .lt. bufsize)
then 741 buf(bufc).faiidx= faiidx
742 buf(bufc).rinmu = rinmu
745 buf(bufc).
t = atrack.
t 746 buf(bufc).wx=-atrack.vec.w.r(1)
747 buf(bufc).wy=-atrack.vec.w.r(2)
750 write(fnodat) bufc, buf
754 if(wz .lt. 0.999)
then 756 *
'(6i3, 1pE11.3, 0p,f6.1,1p2E11.3,0p, 2f8.4,f10.6)')
758 * atrack.
p.
charge, ridx, faiidx,
761 * -atrack.vec.w.r(1), -atrack.vec.w.r(2), wz
763 write(*,
'(6i3,1pE11.3, 0p, f6.1, 1p2E11.3, a)')
765 * atrack.
p.
charge, ridx, faiidx,
767 * ek, atrack.
t,
" 0 0 1" 776 if(id .eq. 2 .and. code .le. 6 .and. i .gt. 0)
then 779 call kwhist(rspec(codex, i), sr, enhance )
782 * .and. atrack.
p.
charge .ne. 0 )
then 785 call kwhist(lossrspec(3,i), sr, de)
786 if( code .eq.
kelec )
then 788 call kwhist(lossrspec(1,i), sr, de)
791 call kwhist(lossrspec(2,i), sr, de)
799 call kwhist2( arspec(codex, i), sngl(fai), sr, enhance)
802 if(id .eq. 2 .and. code .le. 3 .and. i .gt. 0)
then 805 call kwhist2( respec(code, i), sr, ek, enhance)
809 call kwhist2( rzspec(code, i), sr,
810 * sngl(wz), enhance )
813 if( wz .lt. 1.0
d0 )
then 815 if(temp .gt. 0.)
then 818 * atrack.pos.
xyz.
x* atrack.vec.w.r(1) +
819 * atrack.pos.
xyz.
y* atrack.vec.w.r(2) ) /r
827 if( tkzfspec .and. f .lt. 1.0 )
then 828 call kwhist2( zfspec(code, i),
829 * sngl(wz), f, enhance )
833 if( tkefspec .and. f .lt. 1.0 )
then 834 call kwhist2( efspec(code, i),
838 if( tkrfspec .and. f .lt. 1.0 )
then 839 call kwhist2( rfspec(code, i),
844 call kwhist2( rtspec1(code, i), sr,
845 * sngl( atrack.
t ), enhance )
846 call kwhist2( rtspec2(code, i), sr,
847 * sngl( atrack.
t ), enhance )
851 call kwhist3( retspec1(code, i), sr, ek,
852 * sngl( atrack.
t ), enhance)
853 call kwhist3( retspec2(code, i), sr, ek,
854 * sngl( atrack.
t ), enhance)
858 call kwhist3( rezspec(code,i), sr, ek,
859 * sngl(wz), enhance )
863 call kwhist3( rzfspec(code, i), sr,
864 * sngl(wz), f, enhance )
868 call kwhist3( refspec(code, i), sr,
879 write(fnodat) bufc, buf
887 firstz= zfirst.pos.
depth*0.1
893 do i = 1, noofassites
895 asobssites(i).esize = asobssites(i).esize* enhance
897 if(i .gt. 1 .and. i .lt. noofassites )
then 898 dd =(asdepthlist(i+1) - asdepthlist(i-1))/2.0
899 elseif(i .eq. 1)
then 900 dd =(asdepthlist(2) - asdepthlist(1))
902 dd =(asdepthlist(noofassites) -
903 * asdepthlist(noofassites-1))
905 cog = cog + asobssites(i).esize*dd*asdepthlist(i)
906 sumne= sumne +asobssites(i).esize*dd
913 do i = 1, noofassites
914 if( asobssites(i).age .gt.
915 * (2.0-asobssites(noofassites).age))
then 916 if(i .gt. 1 .and. i .lt. noofassites )
then 917 dd =( asdepthlist(i+1) - asdepthlist(i-1))/2.0
918 elseif(i .eq. 1)
then 919 dd =(asdepthlist(2) - asdepthlist(1))
921 dd =(asdepthlist(noofassites) -
922 * asdepthlist(noofassites-1))
925 cog2 = cog2 + asobssites(i).esize*asdepthlist(i)*dd
926 sumne= sumne +asobssites(i).esize*dd
929 if(sumne .gt. 0.)
then 930 cog2 = cog2*0.1/sumne
933 cog2 = asdepthlist(noofassites)*0.1
937 if(fnob .ge. 0 )
then 939 *
'("h ", i4, 3i3, 1pE11.3, 0p 3f11.7, f7.2, 2f7.0)')
940 * eventno, inci.
p.code,
942 * inci.
p.fm.
e, -angle.r(1), -angle.r(2), -angle.r(3),
946 *
'("h ", i4, 3i3, 1pE11.3, 0p 3f11.7, f7.2, 2f7.0)')
947 * eventno, inci.
p.code,
949 * inci.
p.fm.
e, -angle.r(1), -angle.r(2), -angle.r(3),
953 do i = 1, noofassites
955 write(fnob,
'("t ", i3, 2f7.1, 2f6.3, 958 * asdepthlist(i)*0.1, asobssites(i).mu,
959 * asobssites(i).age, asdepthlist(i)*0.1/cog2,
960 *
ng(i),
ne(i),
nmu(i), nhad(i),
961 * asobssites(i).esize, sumeloss(i)
963 write(*,
'("t ", i3, 2f7.1, 2f6.3, 966 * asdepthlist(i)*0.1, asobssites(i).mu,
967 * asobssites(i).age, asdepthlist(i)*0.1/cog2,
968 *
ng(i),
ne(i),
nmu(i), nhad(i),
969 * asobssites(i).esize, sumeloss(i)
972 if(fnob .gt. 0 )
then 982 write(fnonrfai,
'(i2,1pE11.3, 0p,i3, f8.4, 1pE11.3,3i4)')
983 * eventno, e0, nn, cosz, limit,
nrbin,
nfai, ansites
988 write(fnonrfai,
'("rec",f7.1, 4i4)' )
989 * asdepthlist(l)*0.1, l, i, j, k
990 write(fnonrfai,
'(1p10E11.3)')
991 * ( nrfairec(l,k,j,i), l=1,
nrbin )
1000 write(fnonrfai,
'("all",f7.1, 4i4)' )
1001 * asdepthlist(l)*0.1, l, i, j, k
1002 write(fnonrfai,
'(1p10E11.3)')
1003 * ( nrfaiall(l,k,j,i)*enhance, l=1,
nrbin )
1011 write(fnonrfai,
'("dE/dx",f7.1, 3i4)' )
1012 * asdepthlist(l)*0.1, l, i, k
1013 write(fnonrfai,
'(1p10E11.3)')
1014 * ( derfai(l,k,i)*enhance, l=1,
nrbin )
1024 *
'(i3, f7.1, f6.3, f6.3, 1026 *
histdep(i), asdepthlist(j)*0.1,
1027 * asobssites(j).age, asdepthlist(j)*0.1/cog,
1028 * int(asobssites(j).mu), int(cog)
1036 call kwhists( rspec(j,i), 0. )
1040 call kwhistp( rspec(j,i), fno)
1048 call kwhists( lossrspec(j,i), 0. )
1049 call kwhistev(lossrspec(j,i), eventno)
1050 call kwhistid( lossrspec(j, i), evid(i))
1051 call kwhistdir( lossrspec(j, i), ptcl2(j)//
"/")
1052 call kwhistp( lossrspec(j, i), fno)
1063 call kwhists2( arspec(j, i), 0. )
1064 call kwhistev2( arspec(j,i), eventno)
1065 call kwhistid2( arspec(j, i), evid(i))
1066 write(dirstr,
'(a,"/d",i4, "/")')
1067 *
ptcl(j), int( asdepthlist(k)*0.1 )
1068 call kseblk(dirstr,
"|", nstr)
1069 call kwhistdir2(arspec(j, i), dirstr)
1070 call kwhistp2( arspec(j, i), fno)
1078 call kwhists2( respec(j, i), 0. )
1079 call kwhistev2( respec(j,i), eventno)
1080 call kwhistid2( respec(j, i), evid(i))
1082 write(dirstr,
'(a,"/d",i4, "/")')
1083 *
ptcl(j), int( asdepthlist(k)*0.1 )
1084 call kseblk(dirstr,
"|", nstr)
1085 call kwhistdir2(respec(j, i), dirstr)
1086 call kwhistp2( respec(j, i), fno)
1095 call kwhists2( rzspec(j, i), 0.)
1096 call kwhistid2( rzspec(j, i), evid(i))
1098 write(dirstr,
'(a,"/d",i4, "/")')
1099 *
ptcl(j), int( asdepthlist(k)*0.1 )
1100 call kseblk(dirstr,
"|", nstr)
1101 call kwhistdir2( rzspec(j, i), dirstr)
1102 call kwhistp2( rzspec(j, i), fno)
1111 call kwhists2( zfspec(j, i), 0.)
1112 call kwhistid2( zfspec(j, i), evid(i) )
1114 write(dirstr,
'(a,"/d",i4, "/")')
1115 *
ptcl(j), int( asdepthlist(k)*0.1 )
1116 call kseblk(dirstr,
"|", nstr)
1117 call kwhistdir2( zfspec(j, i), dirstr)
1118 call kwhistp2( zfspec(j, i), fno)
1128 call kwhists2( rfspec(j, i), 0.)
1129 call kwhistid2( rfspec(j, i), evid(i))
1131 write(dirstr,
'(a,"/d",i4, "/")')
1132 *
ptcl(j), int( asdepthlist(k)*0.1 )
1133 call kseblk(dirstr,
"|", nstr)
1134 call kwhistdir2(rfspec(j, i), dirstr)
1135 call kwhistp2( rfspec(j, i), fno)
1144 call kwhists2( efspec(j, i), 0.)
1145 call kwhistev2(efspec(j,i), eventno)
1146 call kwhistid2( efspec(j, i), evid(i))
1148 write(dirstr,
'(a,"/d",i4, "/")')
1149 *
ptcl(j), int( asdepthlist(k)*0.1 )
1150 call kseblk(dirstr,
"|", nstr)
1151 call kwhistdir2(efspec(j, i), dirstr)
1152 call kwhistp2( efspec(j, i), fno)
1161 call kwhists2( rtspec1(j,i), 0.)
1162 call kwhistev2( rtspec1(j,i), eventno)
1163 call kwhistid2( rtspec1(j,i), evid(i))
1165 write(dirstr,
'(a,"/d",i4, "/")')
1166 *
ptcl(j), int( asdepthlist(k)*0.1 )
1167 call kseblk(dirstr,
"|", nstr)
1168 call kwhistdir2(rtspec1(j,i), dirstr)
1169 call kwhistp2( rtspec1(j,i), fno)
1174 call kwhists2( rtspec2(j,i), 0.)
1175 call kwhistev2(rtspec2(j,i), eventno)
1176 call kwhistid2( rtspec2(j,i), evid(i))
1178 write(dirstr,
'(a,"/d",i4, "/")')
1179 *
ptcl(j), int( asdepthlist(k)*0.1 )
1180 call kseblk(dirstr,
"|", nstr)
1181 call kwhistdir2(rtspec2(j,i), dirstr)
1182 call kwhistp2( rtspec2(j,i), fno)
1191 call kwhists3( retspec1(j, i), 0.)
1192 call kwhistev3(retspec1(j,i), eventno)
1193 call kwhistid3( retspec1(j, i), evid(i))
1195 write(dirstr,
'(a,"/d",i4,"/")' )
1196 *
ptcl(j), int( asdepthlist(k)*0.1 )
1197 call kseblk(dirstr,
"|", nstr)
1198 call kwhistdir3( retspec1(j, i), dirstr)
1199 call kwhistp3( retspec1(j, i), fno)
1204 call kwhists3( retspec2(j, i), 0.)
1205 call kwhistev3(retspec2(j,i), eventno)
1206 call kwhistid3( retspec2(j, i), evid(i))
1208 write(dirstr,
'(a,"/d",i4, "/")' )
1209 *
ptcl(j), int( asdepthlist(k)*0.1 )
1210 call kseblk(dirstr,
"|", nstr)
1211 call kwhistdir3(retspec2(j, i), dirstr)
1212 call kwhistp3( retspec2(j, i), fno)
1218 if( tkrezspec )
then 1221 call kwhists3( rezspec(j, i), 0.)
1222 call kwhistev3(rezspec(j,i), eventno)
1223 call kwhistid3( rezspec(j, i), evid(i))
1225 write(dirstr,
'(a,"/d",i4,"/")' )
1226 *
ptcl(j), int( asdepthlist(k)*0.1 )
1227 call kseblk(dirstr,
"|", nstr)
1228 call kwhistdir3(rezspec(j, i), dirstr)
1229 call kwhistp3( rezspec(j, i), fno)
1235 if( tkrzfspec )
then 1238 call kwhists3( rzfspec(j, i), 0.)
1239 call kwhistev3(rzfspec(j,i), eventno)
1240 call kwhistid3( rzfspec(j, i), evid(i))
1242 write(dirstr,
'(a,"/d",i4,"/")' )
1243 *
ptcl(j), int( asdepthlist(k)*0.1 )
1244 call kseblk(dirstr,
"|", nstr)
1245 call kwhistdir3(rzfspec(j, i), dirstr)
1246 call kwhistp3( rzfspec(j, i), fno)
1252 if( tkrefspec )
then 1255 call kwhists3( refspec(j, i), 0.)
1256 call kwhistev3(refspec(j,i), eventno)
1257 call kwhistid3( refspec(j, i), evid(i))
1259 write(dirstr,
'(a,"/d",i4, "/")' )
1260 *
ptcl(j), int( asdepthlist(k)*0.1 )
1261 call kseblk(dirstr,
"|", nstr)
1262 call kwhistdir3( refspec(j, i), dirstr)
1263 call kwhistp3( refspec(j, i), fno)
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p *Zfirst pos xyz Zfirst pos xyz *Zfirst pos *Zfirst pos depth
void kwhistid(struct histogram1 *h, char *id)
subroutine howmuch(limit, E0, NN, cosz)
subroutine cqincident(incident, AngleAtObs)
integer lengid integer lengdir character *dir integer kgetenv2 character *numb character *execid character *msg logical takehist save do nsites histdep(i)=0indivdep(i)=0enddo leng
integer nsites ! max real bin
integer function kgetenv2(envname, envresult)
dE dx *! Nuc Int sampling table e
latitude latitude this system is used *****************************************************************! type coord sequence union map real z z in m endmap xyz map real ! latitude in deg is to the north ! longitude in deg is to the east *h ! height in m endmap llh map real ! polar angle ! azimuthal angle *radius ! radial distance endmap sph endunion character *sys ! which system xyz
void kwhistp(struct histogram1 *h, FILE *fno)
subroutine cmudedx(flagN, flagBr, flagPr, Emu, dEdx, dEdxF)
void kwhisti(struct histogram1 *h, float ixmin, float ibinORxmax, int inbin, int itklg)
max ptcl codes in the kelec
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
void kwhistev(struct histogram1 *h, int evno)
logical KeepWeight ! see setupenv sh logical tkarspec ! get lateral histo in a web fai bin logical tkrtspec ! get time histo in each web bin logical SeeLowdE common Zprivatec2 tkrtspec
void kwhist(struct histogram1 *h, float x, float w)
void kwhistdir(struct histogram1 *h, char *dir)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
subroutine cdedxinair(aPtcl, rhoin, dedt, dedtF)
logical KeepWeight ! see setupenv sh logical tkarspec ! get lateral histo in a web fai bin logical tkrtspec ! get time histo in each web bin logical SeeLowdE common Zprivatec2 tkarspec
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p *Zfirst pos xyz Zfirst pos xyz *Zfirst pos *Zfirst pos Zfirst pos height
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
dE dx *! Nuc Int sampling table d
void kwhistc(struct histogram1 *h)
dE dx *! Nuc Int sampling table b
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p charge
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
void kwhistai(struct histogram1 *h, char *title, char *categ, char *dNunit, int logv, float pw, char *label, char *unit)
void kwhists(struct histogram1 *h, float inorm)
*Zfirst p fm *Zfirst p mass
max ptcl codes in the kmuon
subroutine kseblk(text, c, lc)
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
*Zfirst p fm *Zfirst p Zfirst p Zfirst p subcode