c=======================================================================
c
c s h o w a n a l y h c 3 . f
c ---------------------------
c create first a text file containing printed histograms (see
c following tabular) and isecond the bin contents of all histograms.
c compilation:
c gfortran -fbounds-check showanalyhc3.f -o showanalyhc3
c f77 -fbounds-check -m32 showanalyhc3.f -o showanalyhc3
c ifort -C -check bounds showanalyhc3.f -o showanalyhc3
c execution:
c # job_submit -p1 -cp -t2000 -m1000 showanalyhc3.sh000053
c ls -1 DAT00* | grep t -v | n -v > showanalyhc3.i000053
c ./showanalyhc3 < showanalyhc3.i000053 > showanalyhc3.out000053
c mv fort.9 showanalyhc3.fort000053
c get number of histograms:
c grep " -99.8877" showanalyhc3.fort000053 | wc
c-----------------------------------------------------------------------
c input example:
c csk000053/DAT000053-882208370-000000564
c csk000053/DAT000053-882426264-000000803
c csk000053/DAT000053-882850130-000000136
c-----------------------------------------------------------------------
c runh=211285.281 evth=217433.078
c long=52815.2969 evte=3397.39185 rune=3301.33252
c-----------------------------------------------------------------------
c juergen.oehlschlaeger@kit.edu
c-----------------------------------------------------------------------
c
c 121 (1) rec density of gam vs log10(r) "always weight 1"
c 122 (1) rec density of e+ vs log10(r)
c 123 (1) rec density of e- vs log10(r)
c 124 (1) rec density of mu+ vs log10(r)
c 125 (1) rec density of mu- vs log10(r)
c 126 (1) rec density of pi+ vs log10(r)
c 127 (1) rec density of pi- vs log10(r)
c 128 (1) rec density of p vs log10(r)
c 129 (1) rec density of n vs log10(r)
c 130 (1) rec density of cerph vs log10(r)
c 131 (1) rec density of other vs log10(r)
c
c 181 (1) N rec gam vs log10(E) "always weight 1"
c 182 (1) N rec e+ vs log10(E)
c 183 (1) N rec e- vs log10(E)
c 184 (1) N rec mu+ vs log10(E)
c 185 (1) N rec mu- vs log10(E)
c 186 (1) N rec pi+ vs log10(E)
c 187 (1) N rec pi- vs log10(E)
c 188 (1) N rec p vs log10(E)
c 189 (1) N rec n vs log10(E)
c 190 (1) N rec cerph vs log10(E)
c 191 (1) N rec other vs log10(E)
c
c 221 (1) N rec gam vs log10(t) "always weight 1"
c 222 (1) N rec e+ vs log10(t)
c 223 (1) N rec e- vs log10(t)
c 224 (1) N rec mu+ vs log10(t)
c 225 (1) N rec mu- vs log10(t)
c 226 (1) N rec pi+ vs log10(t)
c 227 (1) N rec pi- vs log10(t)
c 228 (1) N rec p vs log10(t)
c 229 (1) N rec n vs log10(t)
c 230 (1) N rec cerph vs log10(t)
c 231 (1) N rec other vs log10(t)
c
c 299 (2) particle codes vs log10(r)
c
c 341 (2) N gam vs log10(E) and log10(r)
c 342 (2) N e+ vs log10(E) and log10(r)
c 343 (2) N e- vs log10(E) and log10(r)
c 344 (2) N mu+ vs log10(E) and log10(r)
c 345 (2) N mu- vs log10(E) and log10(r)
c 346 (2) N pi+ vs log10(E) and log10(r)
c 347 (2) N pi- vs log10(E) and log10(r)
c 348 (2) N p vs log10(E) and log10(r)
c 349 (2) N n vs log10(E) and log10(r)
c 350 (2) N cerph vs log10(E) and log10(r)
c 351 (2) N other vs log10(E) and log10(r)
c
c 361 (2) N gam vs log10(t) and log10(r)
c 362 (2) N e+ vs log10(t) and log10(r)
c 363 (2) N e- vs log10(t) and log10(r)
c 364 (2) N mu+ vs log10(t) and log10(r)
c 365 (2) N mu- vs log10(t) and log10(r)
c 366 (2) N pi+ vs log10(t) and log10(r)
c 367 (2) N pi- vs log10(t) and log10(r)
c 368 (2) N p vs log10(t) and log10(r)
c 369 (2) N n vs log10(t) and log10(r)
c 370 (2) N cerph vs log10(t) and log10(r)
c 371 (2) N other vs log10(t) and log10(r)
c
c 381 (2) N gam vs log10(E) and log10(t)
c 382 (2) N e+ vs log10(E) and log10(t)
c 383 (2) N e- vs log10(E) and log10(t)
c 384 (2) N mu+ vs log10(E) and log10(t)
c 385 (2) N mu- vs log10(E) and log10(t)
c 386 (2) N pi+ vs log10(E) and log10(t)
c 387 (2) N pi- vs log10(E) and log10(t)
c 388 (2) N p vs log10(E) and log10(t)
c 389 (2) N n vs log10(E) and log10(t)
c 390 (2) N cerph vs log10(E) and log10(t)
c 391 (2) N other vs log10(E) and log10(t)
c
c-----------------------------------------------------------------------
program showanalyhc3
c-----------------------------------------------------------------------
c hist vector header:
c qhistos(1,ih) = ident number of the histogram
c qhistos(2,ih) = number of bins (1st dim), nbin
c qhistos(3,ih) = lower bound of 1st dim
c qhistos(4,ih) = upper bound of 1st dim
c qhistos(5,ih) = marker for all logarithmic axes
c qhistos(6,ih) = number of bins (2nd dim), nbi2
c qhistos(7,ih) = lower bound of 2nd dim
c qhistos(8,ih) = upper bound of 2nd dim
c qhistos(9,ih) = -99.8877
c qhistos(10,ih) = not used and not written to unit 9
c qhistos(11,ih) = first element of histogram contents
c hist vector trailer 1-dim: not written to unit 9
c qhistos(10+nbin,ih) = last element of histogram contents
c qhistos(10+nbin+3,ih) = min index of 1st dim with entries
c qhistos(10+nbin+4,ih) = max index of 1st dim with entries
c qhistos(10+nbin+5,ih) = minimum of number of entries
c qhistos(10+nbin+6,ih) = maximum of number of entries
c hist vector trailer 2-dim: not written to unit 9
c qhistos(10+nbin*nbi2,ih) = last element of histogram contents
c qhistos(10+nbin*nbi2+3,ih) = min index of 1st dim with entries
c qhistos(10+nbin*nbi2+4,ih) = max index of 1st dim with entries
c qhistos(10+nbin*nbi2+5,ih) = minimum of number of entries
c qhistos(10+nbin*nbi2+6,ih) = maximum of number of entries
c qhistos(10+nbin*nbi2+7,ih) = min index of 2nd dim with entries
c qhistos(10+nbin*nbi2+8,ih) = max index of 2nd dim with entries
c-----------------------------------------------------------------------
implicit double precision (a-h,o-z), integer (i-n)
parameter (nblklen=312,nhigh=20)
character chtitle(480)*40,chaziff(0:9)*1
double precision qhistos(10020,480)
common /histch/ chtitle
common /histos/ qhistos,rpromin,rpromax
+ ,eparmin,eparmax
+ ,xmin348,xmax348,ymin348,ymax348
common /buffer/ outbuf(nblklen,21)
common /buffch/ crunh,crune,cevth,cevte,clong
character chline(0:101),chaxis(0:101),chlin2(100),cheight(50)
character*4 crunh,crune,cevth,cevte,clong
real ddata(6552),outbuf,chablk(5),curpar
equivalence(chablk(1),crunh)
data crunh/'RUNH'/, crune/'RUNE'/, clong/'LONG'/
data cevth/'EVTH'/, cevte/'EVTE'/
data chtitle/480*' '/
data chaziff/' ','1','2','3','4','5','6','7','8','9'/
data cheight/'a','b','c','d','e','f','g','h','i','j',
+ 'k','l','m','n','o','p','q','r','s','t',
+ 'u','v','w','x','y',
+ 'A','B','C','D','E','F','G','H','I','J',
+ 'K','L','M','N','O','P','Q','R','S','T',
+ 'U','V','W','X','Y'/
c-----------------------------------------------------------------------
c physical constants and run parameters used for simulation
common / phys / obslev(10),enspec(3),ffegs,ffnkg,cut(4),
* consts(230),aatm(5),batm(5),catm(5),flag(4)
c-----------------------------------------------------------------------
c additional run parameters in analysis job
c ishow/ishowu: no of showers read/used
c ishowf: no of showers read from current file
c iblk : number of blocks read
c ifile : number of current data file
c iret1/2: flags for error or runstop in subroutines particles/header
c xoff, yoff coordinate correction for inclined showers
common /runpar/ heighp,thetap,phip,xoff,yoff,tcen,dintmod,dintcrs,
* nobslv,ishow,ishowu,ishowf,iblk,ifile,iret1,iret2,
* rmax,lfnkg,lfegs
double precision heighp,thetap,phip,tcen,dintmod,dintcrs,rmax
logical*1 lfnkg,lfegs
double precision pama,pi
common /pconst/pama(6000),pi
data pi/3.1415926535979d0/
c-----------------------------------------------------------------------
c curpar = current particle :
c curpar (1) = ident * 1000 + igen*10 + iobslv
c ident = identification according to geant
c igen = generation
c iobslv = observation level
c curpar (2,3,4) = px,py,pz in GeV/c
c curpar (5,6,7) = x,y,t in cm, ns
c modified for cerenkov bunches
c curpar(1) = 99*10**5 + nphot*10 + 1
c curpar(2,3) = x and y coordinate in cm
c curpar(4,5) = direction cosini to x and y axis
c curpar(6) = t in ns
c curpar(7) = height of emission in cm
c curpar(8) = weight of particle
c epart,ppart,r = energy, momentum, distance to shower axis
c ityp = particle group
c tf = arrival time of shower front at point(x,y)
c costhe, phi = cosine of polar angle, azimuth angle
common /partic/ curpar(8),ident,igen,iobslv,ityp,
* epart,ppart,r,tf,costhe,phi
c-----------------------------------------------------------------------
c shower parameters
c factors needed to project particles into the plane perp. to shower axis
common /shower/ profak1,profak2,profak3,
* ccla
double precision profak1,profak2,profak3
character*5 ccla(11)
data ccla / 'gam ','e+ ','e- ','mu+ ','mu- ',
* 'pi+ ','pi- ','p ','n ','cerph','other' /
parameter (nmax=50000)
common /steer/ lunmon,lundeb,lunin,lunhis,nfiles
common /steerc/ chfile(nmax),chhist,chlong
character*80 cdat,chfile,chhist,chlong
logical lexist,first/.true./
cheight(nhigh) = '@'
iblklen = 273
eparmin = 1.
eparmax = 1.
rpromin = 1.
rpromax = 1.
xmin348 = 1.
xmax348 = 1.
ymin348 = 1.
ymax348 = 1.
c-----------------------------------------------------------------------
c set defaults
lunmon = 6
lundeb = 6
lunin = 3
lunhis = 6
nobslv = 1
nfiles = 0
do ihist=1,480
do i=1,10020
qhistos(i,ihist) = 0.d0
enddo
enddo
do i=1,nmax
chfile(i) = ' '
enddo
c initialize particle masses
call pamaf
c--read run parameters including file names-----------------------------
do ifile=1,nmax
read(*,'(a)',end=427,err=427) cdat
chfile(ifile) = cdat
enddo
427 continue
nfiles = ifile - 1
nfimod = 10
if ( nfiles .gt. 10000 ) then
nfimod = 200
elseif ( nfiles .gt. 1000 ) then
nfimod = 50
elseif ( nfiles .gt. 300 ) then
nfimod = 20
endif
c--read first particle data file to test record length------------------
open(lunin,file=chfile(1),status='old',
* form='unformatted',access='sequential',err=997)
read(lunin) (ddata(i),i=1,5432)
if ( 217433.0 .lt. ddata(273+1) .and.
+ ddata(273+1) .lt. 217433.2 ) then
write(*,100)
100 format(/,' ===========================================',/,
* ' standard CORSIKA particle data analysis',/,
* ' ===========================================')
iblklen = 273
elseif ( 217433.0 .lt. ddata(312+1) .and.
+ ddata(312+1) .lt. 217433.2 ) then
write(*,101)
101 format(/,' ===========================================',/,
* ' thinning CORSIKA particle data analysis',/,
* ' ===========================================')
iblklen = 312
endif
close(lunin)
c--read particle data files---------------------------------------------
do ifile=1,nfiles
if ( ifile .lt. 10 )
+ write(*,'(i6,''. file: '',a)') ifile,chfile(ifile)
if ( ifile .eq. 11 )
+ write(*,'('' ....................'')')
inquire(file=chfile(ifile),exist=lexist)
if ( .not.lexist ) then
write(*,'(10x,''nr.'',i5,4x,a44,'' NOT found.'')')
+ ifile,chfile(ifile)
chfile(ifile) = 'dummy'
goto 431
endif
open(lunin,file=chfile(ifile),status='old',
* form='unformatted',access='sequential',err=997)
read(lunin,end=429,err=430) (ddata(i),i=1,iblklen*21)
if ( ddata(1) .eq. 0. .and. ddata(iblklen+1) .eq. 0. ) then
write(*,'(10x,''nr.'',i5,4x,a44,'' only ZEROES.'')')
+ ifile,chfile(ifile)
chfile(ifile) = 'dummy'
endif
goto 431
429 continue
write(*,'(10x,''nr.'',i5,4x,a44,'' EOF detected.'')')
+ ifile,chfile(ifile)
chfile(ifile) = 'dummy'
goto 431
430 continue
write(*,'(10x,''nr.'',i5,4x,a44,'' ERR detected.'')')
+ ifile,chfile(ifile)
chfile(ifile) = 'dummy'
431 continue
close(lunin)
enddo
if ( nfiles .eq. 10 .or. ifile .gt. 10 )
+ write(*,'(i6,''. file: '',a)') ifile-1,chfile(ifile-1)
c-----------------------------------------------------------------------
c initialize some variables, histograms
ifile = 0
ishow = 0
ishowu = 0
iblk = 0
call hisini
c-----------------------------------------------------------------------
c begin of loop over all files
10 continue
ifile = ifile + 1
ishowf = 0
c check on correctly read file and quantities
if ( chfile(ifile)(1:5) .eq. 'dummy' ) then
write(*,'(10x,''nr.'',i5,'' renamed to `dummy`.'')') ifile
goto 10
endif
c open current file
open(lunin, file=chfile(ifile), status='old',
* form='unformatted', access='sequential')
if ( mod(ifile,nfimod) .eq. 0 )
+ write(*,'(10x,''nr.'',i5,'' in progress: '',a)')
+ ifile,chfile(ifile)
c read run header (i.e. first data record)
iret2 = 0
read(lunin,end=11,err=22) (ddata(i),i=1,iblklen*21)
iblk = iblk + 1
do isb=1,21
do ib=1,iblklen
outbuf(ib,isb) = ddata(ib+iblklen*(isb-1))
enddo
enddo
c not data tape start marker read
if ( .not. ( 211285.2 .lt. outbuf(1,1) .and.
+ outbuf(1,1) .lt. 211285.4 ) ) then
write(*,*) ' wrong header line of new run'
iret2 = 1
goto 88
endif
c define and print primary particle
if ( ifile .eq. 1 ) then
write(*,219) int(outbuf(2,1)),outbuf(3,2),outbuf(4,2),
+ outbuf(7,2)/100.,outbuf(48,2)/100.,2.e7+outbuf(3,1),
+ outbuf(4,1)
219 format(/,' simulation run number ',i15.6,'.',/,
* ' primary particle ',f16.0,/,
* ' energy fixed [GeV]',1p,e22.4,0p,/,
* ' height of 1st int. [m] ',f16.0,/,
* ' observation level [m] ',f16.0,/,
* ' date of simulation ',f16.0,/,
* ' version ',f16.6)
write(*,224) (outbuf(i,1),i=21,24)
224 format(' cuts [GeV] hadron, muon, electr, photon: ',4f9.5)
obslev(1) = outbuf(48,2)
lfnkg = ( ffnkg .eq. 1. )
lfegs = ( ffegs .eq. 1. )
nflran = flag(1)
nfldif = flag(2)
nflpi0 = mod( int(flag(3)),100 )
nflpif = int(flag(3)) / 100
nflche = mod( int(flag(4)),100 )
nfragm = int(flag(4)) / 100
endif
goto 88
11 continue
write(*,*) ' end of file at reading run header '
iret2 = 1
goto 88
22 continue
write(*,*) ' error at reading run header '
iret2 = 1
c = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
88 continue
if ( iret2 .ne. 0 ) then
write(*,*) ' error or eof in reading header '
goto 99
endif
c-----------------------------------------------------------------------
c read data blocks from current file and call analysis routines
iret1 = 0
call particles(iblklen)
if ( iret2 .ne. 0 ) then
write(*,*) ' error or eof in header '
goto 99
endif
c close current file
close(lunin)
if ( iret1 .ne. 0 ) then
if ( iret1 .eq. 4 ) write(*,*) ' unexpected end-of-file at ',
+ ifile,'th file: ',chfile(ifile)
c write(*,*) 'eof detected when reading lunin '
c iret1 = 5
c return
c write(*,*) ' error at reading lunin '
c iret1 = 3
c return
write(*,1001) ifile,nfiles
1001 format(' error or end in particles at ',
* i3,'. file of',i5,' input files ')
endif
c does another input file exist?
if ( nfiles .gt. ifile ) goto 10
c-----------------------------------------------------------------------
c end of event loop
c-----------------------------------------------------------------------
99 continue
write(*,'(/,5x,''total number of records'',i10,/)') iblk
c-----------------------------------------------------------------------
c printer plots of histograms at the end of the run
c histext
c nhist, nbin, xlow, xupr, vlog, [nbi2, ylow, yupr, dummy]
c plots / table of contents
c histogram quantities will be written to fort.9 (without statistics)
c-----------------------------------------------------------------------
do ihist=121,389
if ( chtitle(ihist)(1:2) .ne. ' ' ) then
c - - - - - - print title and optionally contents of the histogram:
nbin = qhistos(2,ihist)
nbi2 = max(1,int(qhistos(6,ihist)))
write(*,'(55('' =''))')
write(*,'(1x,a40)') chtitle(ihist)
write(9,'(1x,a40)') chtitle(ihist)
if ( qhistos(6,ihist) .eq. 0. ) then ! 1-dim.
write(*,'(2(0p,2f8.0,1p,2g12.4))')
+ (qhistos(i,ihist),i=1,5)
write(9,'(2(0p,2f8.0,1p,2g12.4),0p,f14.4)')
+ (qhistos(i,ihist),i=1,5), 0.,-99.,-99., -99.8877
write(9,'(1p,10e12.4)') (qhistos(i,ihist),i=11,10+nbin)
else ! 2-dim.
write(*,'(2(0p,2f8.0,1p,2g12.4))')
+ (qhistos(i,ihist),i=1,8)
write(9,'(2(0p,2f8.0,1p,2g12.4),0p,f14.4)')
+ (qhistos(i,ihist),i=1,8), -99.8877
write(9,'(1p,10e12.4)') (qhistos(i,ihist),i=11,10+nbin*nbi2)
endif
c - - - - - - calculate minimum and maximum of the histogram:
qhistos(10+nbin*nbi2+3,ihist) = -1.
qhistos(10+nbin*nbi2+4,ihist) = -1.
qhmin = qhistos(11,ihist)
qhmax = qhistos(11,ihist)
do ib=11,10+nbin*nbi2
if ( qhistos(ib,ihist) .gt. qhmax )
+ qhmax = qhistos(ib,ihist)
if ( qhistos(ib,ihist) .lt. qhmin )
+ qhmin = qhistos(ib,ihist)
if ( qhistos(10+nbin*nbi2+3,ihist) .eq. -1. .and.
+ qhistos(ib,ihist) .gt. 0. ) then
qhistos(10+nbin*nbi2+3,ihist) = -10 + ib
endif
if ( qhistos(10+nbin*nbi2+4,ihist) .eq. -1. .and.
+ qhistos(21+nbin*nbi2-ib,ihist) .gt. 0. ) then
qhistos(10+nbin*nbi2+4,ihist) = 11. + nbin*nbi2 - ib
endif
enddo
qhistos(10+nbin*nbi2+5,ihist) = qhmin
qhistos(10+nbin*nbi2+6,ihist) = qhmax
if ( qhmin .eq. qhmax ) goto 109
c-----------------------------------------------------------------------
if ( ihist .lt. 299 ) then
c - - - - - - - calculate stretching factor (5 or 2):
qhlog = log10(qhmax)
if ( qhlog .gt. 0. ) then
ihlog = int(qhlog)
qhlog = qhlog - ihlog
qscal = 1. - ihlog
else
ihlog = 1. + int(-qhlog)
qhlog = qhlog + ihlog
qscal = 1.*ihlog
endif
if ( qhlog .gt. 0.69897d0 ) then
qfact = 1.
elseif ( qhlog .gt. 0.30103d0 ) then
qfact = 2.
else
qfact = 5.
endif
qscal = 10.d0**qscal
write(*,'(42x,''hmin='',1p,g12.5,'' hmax='',g12.5,
+ 9x,''(content *'',g9.2,'')'')') qhmin,qhmax,1.d0/qscal
c - - - - - - - plot title axis and scale quantities:
do i=1,100
chaxis(i) = '-'
enddo
if ( qfact .eq. 5. ) then
do i=10,90,10
chaxis(i) = 'o'
enddo
chaxis(25) = '+'
chaxis(50) = '|'
chaxis(75) = '+'
chaxis(100) = '|'
elseif ( qfact .eq. 2. ) then
do i=10,90,20
chaxis(i) = '+'
chaxis(i+10) = '|'
enddo
else
do i=10,100,10
chaxis(i) = '|'
enddo
endif
if ( qfact .eq. 5. ) then
write(*,'(i8,4i25)') 0,(i,i=5,20,5)
elseif ( qfact .eq. 2. ) then
write(*,'(i8,5i20)') 0,(10*i,i=1,5)
else
write(*,'(i8,10i10)') 0,(10*i,i=1,10)
endif
write(*,'(6x,'' !'',100a1)') (chaxis(i),i=1,100)
c - - - - - - - plot histogram contents (to the right):
do ib=11,10+nbin
do is=1,100
chline(is) = ' '
enddo
if ( qfact .eq. 5. ) then
do is=25,100,25
chline(is) = '|'
enddo
elseif ( qfact .eq. 2. ) then
do is=20,100,20
chline(is) = '|'
enddo
else
do is=10,100,10
chline(is) = '|'
enddo
endif
is = qscal * qfact * qhistos(ib,ihist)
if ( is .eq. 0 .and. qhistos(ib,ihist) .gt. 0. ) is = 1
chline(is) = '*'
write(*,'(f6.2,'' |'',100a1)')
+ qhistos(3,ihist)+(qhistos(4,ihist)-qhistos(3,ihist))/
+ qhistos(2,ihist)*(ib-11),(chline(i),i=1,is)
enddo
c - - - - - - - plot closing axis and scale quantities:
write(*,'(6x,'' !'',100a1)') (chaxis(i),i=1,100)
if ( qfact .eq. 5. ) then
write(*,'(i8,4i25)') 0,(i,i=5,20,5)
elseif ( qfact .eq. 2. ) then
write(*,'(i8,5i20)') 0,(10*i,i=1,5)
else
write(*,'(i8,10i10)') 0,(10*i,i=1,10)
endif
c-----------------------------------------------------------------------
elseif ( ihist .eq. 299 ) then
c - - - - - - - hist.299: plot title axis and scale quantities:
write(*,'(42x,''hmin='',1p,g12.5,'' hmax='',g12.5)')
+ qhmin,qhmax
! - - - - - linear scale for contents:
qfact = 0.99999d0 * nhigh / qhmax
!write(*,'(59x,''hmax/'',i2,''='',1p,g12.5)')nhigh,qhmax/nhigh
! - - - - - logarithmic scale for contents:
qfact = (qhmax+1.d-2)**(1.d0/nhigh)
write(*,'(62x,''hmax^(1/'',i2,'')='',1p,g12.5)')
+ nhigh,qfact
if ( qhmin .eq. qhmax ) goto 105
! - - - - - print only for qhmax > qhmin:
do i=1,100
chaxis(i) = '-'
enddo
do i=20,100,20
chaxis(i) = '|'
chaxis(i-10) = '+'
enddo
write(*,'(f9.1,10f10.1)') (qhistos(7,ihist)+
+ (qhistos(8,ihist)-qhistos(7,ihist))/10.*i,i=0,10)
write(*,'(6x,'' !'',100a1)') (chaxis(i),i=1,100)
c - - - - - - - plot 2-dim-histogram contents using cheight steps:
do ib=0,nbin-1
do iy=1,nbi2
chlin2(iy) = ' '
if ( ib .gt. 0 ) then
if ( qhistos(10+ib+(iy-1)*nbin,ihist) .gt. 0. ) then
! - - - - - - - - - - linear steps:
! ih = 1. + qhistos(10+ib+(iy-1)*nbin,ihist) * qfact
! if ( ih .eq. 0 ) ih = 1
! chlin2(iy) = cheight(ih)
! if ( qhistos(10+ib+(iy-1)*nbin,ihist).le.9. ) then
! ih = qhistos(10+ib+(iy-1)*nbin,ihist)
! chlin2(iy) = chaziff(ih)
! endif
! - - - - - - - - - - logarithmic steps:
qsqua = (qhmax+1.d-2)
ih = 21
104 continue
qsqua = qsqua / qfact
ih = ih - 1
if ( qhistos(10+ib+(iy-1)*nbin,ihist) .le. qsqua )
+ goto 104
if ( ih .eq. 0 ) ih = 1
chlin2(iy) = cheight(ih)
endif
endif
enddo
write(*,'(i6,'' |'',100a1,''|'')')
+ ib,(chlin2(iy),iy=1,nbi2)
enddo
c - - - - - - - plot closing axis and scale quantities:
write(*,'(6x,'' !'',100a1)') (chaxis(i),i=1,100)
write(*,'(f9.1,10f10.1)') (qhistos(7,ihist)+
+ (qhistos(8,ihist)-qhistos(7,ihist))/10.*i,i=0,10)
105 continue
c-----------------------------------------------------------------------
else ! ihist > 299
if ( ihist .eq. 348 ) then
write(*,'(42x,''xmin='',1p,g12.5,'' xmax='',g12.5)')
+ xmin348,xmax348
write(*,'(42x,''ymin='',1p,g12.5,'' ymax='',g12.5)')
+ ymin348,ymax348
endif
c - - - - - - - print histogram 348 (2) N p vs log10(E) and log10(r)
c - - - - - - - - - - - - - - - 368 (2) N p vs log10(t) and log10(r)
c - - - - - - - - - - - - - - - 388 (2) N p vs log10(E) and log10(t)
write(*,'(42x,''hmin='',1p,g12.5,'' hmax='',g12.5)')
+ qhmin,qhmax
write(*,'(59x,''hmax/'',i2,''='',1p,g12.5)')
+ nhigh,qhmax/nhigh
qfact = 0.99999d0 * nhigh / qhmax
do i=1,100
chaxis(i) = '-'
enddo
do i=20,100,20
chaxis(i) = '|'
chaxis(i-10) = '+'
enddo
write(*,'(f9.1,10f10.1)') (qhistos(7,ihist)+
+ (qhistos(8,ihist)-qhistos(7,ihist))/10.*i,i=0,10)
write(*,'(6x,'' !'',100a1)') (chaxis(i),i=1,100)
do ib=1,nbin
do iy=1,nbi2
chlin2(iy) = ' '
if ( qhistos(10+ib+(iy-1)*nbin,ihist) .gt. 0. ) then
ih = 1. + qhistos(10+ib+(iy-1)*nbin,ihist) * qfact
chlin2(iy) = cheight(ih)
endif
enddo
write(*,'(f6.1,'' |'',100a1,''|'')')
+ qhistos(3,ihist)+(qhistos(4,ihist)-qhistos(3,ihist))/
+ qhistos(2,ihist)*(ib-1),(chlin2(iy),iy=1,nbi2)
enddo
write(*,'(6x,'' !'',100a1)') (chaxis(i),i=1,100)
write(*,'(f9.1,10f10.1)') (qhistos(7,ihist)+
+ (qhistos(8,ihist)-qhistos(7,ihist))/10.*i,i=0,10)
endif ! end-of ihist > 299
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c - - - - - end-of individual histograms.
109 continue
endif
enddo
goto 999
c-----------------------------------------------------------------------
997 write(*,'(a)') ' problems with open data file ==> run is stopped'
write(*,'(a,4x,a)') ' chfile =', chfile(ifile)
999 continue
stop
end
c=======================================================================
subroutine hisini
c-----------------------------------------------------------------------
c book histograms
c-----------------------------------------------------------------------
c definitions
implicit double precision (a-h,o-z), integer (i-n)
parameter (nblklen=312)
character chtitle(480)*40
double precision qhistos(10020,480)
common /histch/ chtitle
common /histos/ qhistos,rpromin,rpromax
+ ,eparmin,eparmax
+ ,xmin348,xmax348,ymin348,ymax348
common /buffer/ outbuf(nblklen,21)
common /buffch/ crunh,crune,cevth,cevte,clong
character*4 crunh,crune,cevth,cevte,clong
real outbuf,chablk(5),curpar
equivalence(chablk(1),crunh)
double precision pama,pi
common /pconst/pama(6000),pi
common /partic/ curpar(8),ident,igen,iobslv,ityp,
* epart,ppart,r,tf,costhe,phi
common / phys / obslev(10),enspec(3),ffegs,ffnkg,cut(4),
* consts(230),aatm(5),batm(5),catm(5),flag(4)
common /runpar/ heighp,thetap,phip,xoff,yoff,tcen,dintmod,dintcrs,
* nobslv,ishow,ishowu,ishowf,iblk,ifile,iret1,iret2,
* rmax,lfnkg,lfegs
double precision heighp,thetap,phip,tcen,dintmod,dintcrs,rmax
logical*1 lfnkg,lfegs
common /shower/ profak1,profak2,profak3,
* ccla
double precision profak1,profak2,profak3
character*5 ccla(11)
parameter (nmax=50000)
common /steer/ lunmon,lundeb,lunin,lunhis,nfiles
common /steerc/ chfile(nmax),chhist,chlong
character*80 chfile,chhist,chlong,htit
double precision tl, ul, tl2, ul2, ulog, area
integer nb,nb2
c 1-dim histograms
c histograms as function of log10(r)
nb = 100
tl = 0.
ul = 5.
c record density
htit = 'rec density of vs log10(r)'
nn = 120
ulog = 1.
do i=1,09
htit(17:21) = ccla(i)
call hbook1(nn+i,htit,nb,tl,ul,ulog)
enddo
c histograms as function of log10(E)
nb = 100
tl = -4.
ul = 6.
c energy distribution for recs
htit = 'N rec vs log10(E)'
nn = 180
ulog = 1.
do i=1,09
htit(7:11) = ccla(i)
call hbook1(nn+i,htit,nb,tl,ul,ulog)
enddo
c histograms as function of log10(t)
nb = 100
tl = -5.
ul = 5.
c time distribution for recs
htit = 'N rec vs log10(t)'
nn = 220
ulog = 1.
do i=1,09
htit(7:11) = ccla(i)
call hbook1(nn+i,htit,nb,tl,ul,ulog)
enddo
c 2-dim histograms
c particle numbers as function of particle code and r
nb = 30
tl = 0.
ul = 30.
nb2 = 100
tl2 = -4.
ul2 = 6.
ulog = 2.
htit = 'particle codes vs log10(r)'
nn = 299
call hbook2(nn,htit,nb,tl,ul,nb2,tl2,ul2,ulog)
c particle numbers as function of energy and r
nb = 80 ! 100
tl = -4. ! -4.
ul = 4. ! 6.
nb2 = 100
tl2 = -4.
ul2 = 6.
nn = 340
ulog = 3.
htit = 'N vs log10(E) and log10(r)'
do i=1,09
htit(3:7) = ccla(i)
call hbook2(nn+i,htit,nb,tl,ul,nb2,tl2,ul2,ulog)
enddo
c particle numbers as function of time and r
nb = 100
tl = -5.
ul = 5.
nb2 = 100
tl2 = -4.
ul2 = 6.
htit = 'N vs log10(t) and log10(r)'
nn = 360
ulog = 3.
do i=1,09
htit(3:7) = ccla(i)
call hbook2(nn+i,htit,nb,tl,ul,nb2,tl2,ul2,ulog)
enddo
c particle numbers as function of energy and t
nb = 80 ! 100
tl = -4. ! -4.
ul = 4. ! 6.
nb2 = 100
tl2 = -5.
ul2 = 5.
htit = 'N vs log10(E) and log10(t)'
nn = 380
ulog = 3.
do i=1,09
htit(3:7) = ccla(i)
call hbook2(nn+i,htit,nb,tl,ul,nb2,tl2,ul2,ulog)
enddo
return
end
c=======================================================================
subroutine particles(iblklen)
c-----------------------------------------------------------------------
c reads all data blocks from unit 'lunin' for current input file
c and calls analyses routines for showers to be used
c-----------------------------------------------------------------------
implicit double precision (a-h,o-z), integer (i-n)
parameter (nblklen=312)
character chtitle(480)*40
double precision qhistos(10020,480)
common /histch/ chtitle
common /histos/ qhistos,rpromin,rpromax
+ ,eparmin,eparmax
+ ,xmin348,xmax348,ymin348,ymax348
common /buffer/ outbuf(nblklen,21)
common /buffch/ crunh,crune,cevth,cevte,clong
character*4 crunh,crune,cevth,cevte,clong
real ddata(6552),outbuf,chablk(5),curpar
equivalence(chablk(1),crunh)
c-----------------------------------------------------------------------
double precision pama,pi
common /pconst/pama(6000),pi
c-----------------------------------------------------------------------
c variables are the same as used in corsika :
c czx1 = lateral distribution in x - direction for 1. level
c ...
c czxy2 = lateral distribution in xy- direction for 2. level etc.
c ...
c nenkg1(i) = Number of electrons in steps of 100 g/cm**2
c age1 (i) = age of shower in steps of 100 g/cm**2
c dist (i) = distance bins for NKG lateral dist
c lage1 (i) = local age level 1
c nenkg2(i) = nenkg1(i); dummy
c age2 (i) = age1(i); dummy
c tlev (i) = levels for NKG output in g/cm^2
c cmlev (i) = levels for NKG output in cm
c dummy = unused words
parameter (npd=10)
common /nkg/czx1 (-npd:npd), czy1 (-npd:npd),
* czxy1(-npd:npd), czyx1(-npd:npd),
* czx2 (-npd:npd), czy2 (-npd:npd),
* czxy2(-npd:npd), czyx2(-npd:npd),
* nenkg1(npd), age1(npd), dist(npd), lage1(npd),
* tlev(npd),cmlev(npd),
* dummy(2*npd)
real*4 nkgpar(248),nenkg1,age1,dist,lage1,tlev,cmlev
equivalence (nkgpar(1),czx1(-npd))
common /partic/ curpar(8),ident,igen,iobslv,ityp,
* epart,ppart,r,tf,costhe,phi
common / phys / obslev(10),enspec(3),ffegs,ffnkg,cut(4),
* consts(230),aatm(5),batm(5),catm(5),flag(4)
common /runpar/ heighp,thetap,phip,xoff,yoff,tcen,dintmod,dintcrs,
* nobslv,ishow,ishowu,ishowf,iblk,ifile,iret1,iret2,
* rmax,lfnkg,lfegs
double precision heighp,thetap,phip,tcen,dintmod,dintcrs,rmax
logical*1 lfnkg,lfegs
common /shower/ profak1,profak2,profak3,
* ccla
double precision profak1,profak2,profak3
character*5 ccla(11)
parameter (nmax=50000)
common /steer/ lunmon,lundeb,lunin,lunhis,nfiles
common /steerc/ chfile(nmax),chhist,chlong
character*80 chfile,chhist,chlong
character*8 use(2) /'used ','ignored '/
logical*1 first /.true./
c-----------------------------------------------------------------------
c beginning of new shower
c-----------------------------------------------------------------------
isb = 2
1 continue
ishow = ishow + 1
ishowf = ishowf + 1
c decide whether shower is used or ignored
if ( ishowf .lt. 0 ) then
iu = 2
else
iu = 1
ishowu=ishowu + 1
endif
if ( iu .eq. 1 ) then
c calculate offsets of shower fronts for the observation level
heighp = outbuf( 7,isb)
thetap = outbuf(11,isb)
phip = outbuf(12,isb)
xoff = -(heighp-obslev(1)) * tan(thetap) * cos(phip)
yoff = -(heighp-obslev(1)) * tan(thetap) * sin(phip)
c get maximal radius for radial thinning (in m)
rmax = 1.d-2 * outbuf(152,isb)
c arrival time of shower core at observation level
tcen = sqrt((heighp-obslev(1))**2+xoff**2+yoff**2)/29.9792425
c factors for projection to the plane perpendicular to the shower plane
profak1 = cos(phip)**2*cos(thetap)+sin(phip)**2
profak2 = cos(phip)*sin(phip)*(cos(thetap)-1.)
profak3 = sin(phip)**2*cos(thetap)+cos(phip)**2
c get interaction model and cross sections
c venus = 1, qgsjet = 2, sibyll = 3, dpmjet = 4, hdpm = 5
if ( outbuf(76,isb) .eq. 1. ) then
dintmod = 1.
elseif ( outbuf(139,isb) .eq. 1. ) then
dintmod = 3.
elseif ( outbuf(141,isb) .eq. 1. ) then
dintmod = 2.
elseif ( outbuf(143,isb) .eq. 1. ) then
dintmod = 4.
else
dintmod = 5.
endif
if ( outbuf(145,isb) .eq. 1. ) then
dintcrs = 1.
elseif ( outbuf(140,isb) .eq. 1. ) then
dintcrs = 3.
elseif ( outbuf(142,isb) .eq. 1. ) then
dintcrs = 2.
elseif ( outbuf(144,isb) .eq. 1. ) then
dintcrs = 4.
else
dintcrs = 5.
endif
endif
c-----------------------------------------------------------------------
c get next subblock
2 continue
isb = isb + 1
if ( isb .gt. 21 ) then
read(lunin,end=112,err=114) (ddata(i),i=1,iblklen*21)
iblk = iblk + 1
do isb=1,21
do ib=1,iblklen
outbuf(ib,isb) = ddata(ib+iblklen*(isb-1))
enddo
enddo
isb = 1
endif
c shower finished ?
if ( 3397.3 .lt. outbuf(1,isb) .and.
+ outbuf(1,isb) .lt. 3397.5 ) goto 5
lh = 0
if ( iu .eq. 2 ) lh = 38
c-----------------------------------------------------------------------
c next particle. transfer particle of datablock to current particle
3 continue
iblpart = 7
if ( iblklen .eq. 312 ) iblpart = 8
do i=1,iblpart
curpar(i) = outbuf(iblpart*lh+i,isb)
enddo
c shower finished if 0 in the word where particle code should be
if ( curpar(1) .eq. 0. ) then
goto 4
endif
c-----------------------------------------------------------------------
c call standard analysis routines and
c optional user analysis routines for each particles
c if the shower is to be used for the analysis
if ( iu .eq. 1 ) then
call analys(iblklen)
endif
if ( lh .eq. 38 ) goto 2
lh = lh + 1
goto 3
c-----------------------------------------------------------------------
c shower finished
c-----------------------------------------------------------------------
4 continue
isb = isb + 1
if ( isb .gt. 21 ) then
read(lunin,end=112,err=114) (ddata(i),i=1,iblklen*21)
iblk = iblk + 1
do isb=1,21
do ib=1,iblklen
outbuf(ib,isb) = ddata(ib+iblklen*(isb-1))
enddo
enddo
isb = 1
endif
c next subblock a LONG block?
if ( 52815.2 .lt. outbuf(1,isb) .and.
+ outbuf(1,isb) .lt. 52815.4 ) then
write(*,*) 'LONG block detected, ignore for the moment'
goto 4
endif
if ( .not. ( 3397.3 .lt. outbuf(1,isb) .and.
+ outbuf(1,isb) .lt. 3397.5 ) ) then
iret1 = 7
write(*,*) 'no EVTE = evt end after particle code 0 '
return
endif
5 continue
c print statistics for shower (from end of event block)
if ( ifile .eq. 1 ) write(*,'('' '')')
if ( ifile .eq. 1 .and. nfiles .eq. 1 ) then
write(*,'(i6,''. gam,el,had,mu,sum:'',1p,5e14.6)')
+ ifile,(outbuf(i,isb),i=3,7)
else
if ( ifile .lt. 10 ) then
if ( nfiles .gt. 1 ) then
if ( outbuf(7,isb) .gt. 0. )
+ write(*,'(i6,''. gam,el,had,mu,sum:'',1p,5e14.6)')
+ ifile,(outbuf(i,isb),i=3,7)
endif
elseif ( nfiles .ge. 10 .and. ifile .eq. nfiles ) then
write(*,'(i6,''. gam,el,had,mu,sum:'',1p,5e14.6)')
+ ifile,(outbuf(i,isb),i=3,7)
endif
endif
c-----------------------------------------------------------------------
c define additional shower parameter from nkg-routine
if ( lfnkg ) then
do i=1,248
nkgpar(i) = outbuf(i+7,isb)
enddo
endif
c next shower to be read ?
if ( 0 .eq. ishowf ) then
iret1 = 6
return
endif
c get next subblock
isb = isb + 1
if ( isb .gt. 21 ) then
read(lunin,end=113,err=114) (ddata(i),i=1,iblklen*21)
iblk = iblk + 1
do isb=1,21
do ib=1,iblklen
outbuf(ib,isb) = ddata(ib+iblklen*(isb-1))
enddo
enddo
isb = 1
endif
if ( 217433.0 .lt. outbuf(1,isb) .and.
+ outbuf(1,isb) .lt. 217433.2 ) then
c new event header
goto 1
elseif ( 3301.3 .lt. outbuf(1,isb) .and.
+ outbuf(1,isb) .lt. 3301.5 ) then
c run end found
iret1 = 0
return
endif
c-----------------------------------------------------------------------
c read marks: error end
111 continue
! write(*,*) 'eof detected when reading '
iret1 = 5
return
112 continue
! write(*,*) ' unexpected end-of-file '
ishowu = ishowu - 1
iret1 = 4
return
113 continue
! write(*,*) ' eof at reading lunin after "evte" '
return
114 continue
! write(*,*) ' error at reading '
ishowu = ishowu - 1
iret1 = 3
return
end
c=======================================================================
subroutine analys(iblklen)
c-----------------------------------------------------------------------
c analysis routine: called for each particle from subroutine next
c if shower is used
c-----------------------------------------------------------------------
implicit double precision (a-h,o-z), integer (i-n)
parameter (nblklen=312)
character chtitle(480)*40
double precision qhistos(10020,480)
common /histch/ chtitle
common /histos/ qhistos,rpromin,rpromax
+ ,eparmin,eparmax
+ ,xmin348,xmax348,ymin348,ymax348
common /buffer/ outbuf(nblklen,21)
common /buffch/ crunh,crune,cevth,cevte,clong
character*4 crunh,crune,cevth,cevte,clong
real outbuf,chablk(5),curpar
equivalence(chablk(1),crunh)
c-----------------------------------------------------------------------
double precision pama,pi
common /pconst/pama(6000),pi
common /partic/ curpar(8),ident,igen,iobslv,ityp,
* epart,ppart,r,tf,costhe,phi
common / phys / obslev(10),enspec(3),
* ffegs,ffnkg,cut(4),
* consts(230),aatm(5),batm(5),catm(5),flag(4)
common /runpar/ heighp,thetap,phip,xoff,yoff,tcen,dintmod,dintcrs,
* nobslv,ishow,ishowu,ishowf,iblk,ifile,iret1,iret2,
* rmax,lfnkg,lfegs
double precision heighp,thetap,phip,tcen,dintmod,dintcrs,rmax
logical*1 lfnkg,lfegs
common /shower/ profak1,profak2,profak3,
* ccla
double precision profak1,profak2,profak3
character*5 ccla(11)
parameter (nmax=50000)
common /steer/ lunmon,lundeb,lunin,lunhis,nfiles
common /steerc/ chfile(nmax),chhist,chlong
character*80 chfile,chhist,chlong
double precision dident,rrplog,timlog,timlogw
integer icount /0/
c-----------------------------------------------------------------------
c do not analyze the muon additional information
if ( int(curpar(1)/1000) .eq. 75 .or.
* int(curpar(1)/1000) .eq. 76 ) return
if ( int(curpar(1)/1000) .eq. 85 .or.
* int(curpar(1)/1000) .eq. 86 ) return
if ( int(curpar(1)/1000) .eq. 95 .or.
* int(curpar(1)/1000) .eq. 96 ) return
if ( curpar(1) .lt. 1001. ) return
c is particle from the correct observation level ?
iobslv= 1 ! mod( int(curpar(1)) , 10 )
c unpack information for normal particles
if ( curpar(1) .lt. 1.e6 ) then
igen = mod( int(curpar(1)), 1000 ) / 10
ident = curpar(1) / 1000
px = curpar(2)
py = curpar(3)
pz = curpar(4)
x = curpar(5)
y = curpar(6)
t = curpar(7)
wt = 1.
if ( iblklen .eq. 312 ) wt = curpar(8)
c calculate momentum and energy in GeV
psq = px**2 + py**2 + pz**2
ppart = sqrt( psq )
epart = sqrt( pama(ident)**2 + psq )
eplog = log10(epart)
c calculate cosine of theta and phi
costhe = pz / ppart
if ( px .eq. 0. .and. py .eq. 0. ) then
phi = 0.
else
phi = atan2( py , px )
endif
c unpack information for cherenkov bunches
c no momentum, just direction available
elseif ( curpar(1) .ge. 1.e6 ) then
igen = 0
ident = int(curpar(1)/1.e3)
xnpho = int((curpar(1)-ident*1.e3)/10.)
www = curpar(4)**2 + curpar(5)**2
costhe = -sqrt(1. - www)
theta = acos(costhe)
if ( www .le. 0. ) then
phi = 0.
else
phi = -atan2(curpar(5),curpar(4))
endif
x = curpar(2)
y = curpar(3)
cx = curpar(4) ! direction cosine to x-axis.
cy = curpar(5) ! direction cosine to y-axis.
t = curpar(6)
wt = 1.
if ( iblklen .eq. 312 ) wt = (curpar(1)-99.e5)/10.
curpar(8) = wt
endif
c calculate arrival time of spherical shower front at point(x,y)
tf = sqrt( ( heighp-obslev(1) ) **2 +
* ( x-xoff ) **2 + ( y-yoff ) **2 ) / 29.9792425
c get the arrival time of particle with respect to the spherical shower front
c this a good quantity to get a local time distribution at one detector
trel = t - tf
c however, with trel values one cannot reconstruct a shower front
c since tf is subtracted from the times and tf is dependent on the
c position in x and y.
c To get the correct arrival times at different detectors subtract from
c the simple arrival time t the time tcen when the shower axis hits the ground.
c (see evaluation of tcen in routine next)
c this difference is negative for about half the particles in inclined showers
tabs = t - tcen
if (trel .gt. 0 ) then
timlog = log10(trel)
else
timlog = -1.e30
endif
c calculate distance from shower axis in meters
xx = x * 1.e-2
yy = y * 1.e-2
rr = sqrt(xx**2 + yy**2)
c calculate the radius projected on to the plane perpendicular to the shower axis
xxproj = xx * profak1 + yy * profak2
yyproj = xx * profak2 + yy * profak3
rrproj = sqrt(xxproj**2 + yyproj**2)
c use for plotting the radius in the plane perpendicular to the shower axis
rrplog = log10(rrproj)
c now filling of histograms starts with overall 2dim-histogram
if ( 1 .le. ident .and. ident .lt. 30 ) then
dident = ident
call hfill2(299,dident,rrplog,wt)
endif
c define particle type and count particles
c ...... EGS photons
idi = 11
if ( ident .eq. 1 ) then
idi = 1
c ...... positrons
elseif ( ident .eq. 2 ) then
idi = 2
c ...... electrons
elseif ( ident .eq. 3 ) then
idi = 3
c ...... muons+
elseif ( ident .eq. 5 ) then
idi = 4
c ...... muons-
elseif ( ident .eq. 6 ) then
idi = 5
c ...... pions+
elseif ( ident .eq. 8 ) then
idi = 6
c ...... pions-
elseif ( ident .eq. 9 ) then
idi = 7
c ...... proton/antiproton
elseif ( ident .eq. 14 .or. ident .eq. 15 ) then
idi = 8
if ( rrplog .lt. ymin348 ) ymin348 = rrplog
if ( rrplog .gt. ymax348 ) ymax348 = rrplog
if ( epart .lt. eparmin ) eparmin = epart
if ( epart .gt. eparmax ) eparmax = epart
if ( eplog .lt. xmin348 ) xmin348 = eplog
if ( eplog .gt. xmax348 ) xmax348 = eplog
c ...... neutron/antineutron
elseif ( ident .eq. 13 .or. ident .eq. 25 ) then
idi = 9
c ...... cerenkov photons
elseif ( ident .ge. 9900 ) then
c idi = 10
else
c ...... other particles
c idi = 11
endif
c fill histograms
c (attention: there is a slight inconsistency in that corsika
c applies radial thinning for particles within rmax from the
c shower core in the observation level.
c here however we plot radial distributions where the radius
c is calculated in the plane perpendicular to the shower axis.)
c make the radius cut only for those histograms that do not show
c dependence on radius
call hfill1(120+idi,rrplog,0.d0,wt)
call hfill2(340+idi,eplog,rrplog,wt)
call hfill2(360+idi,timlog,rrplog,wt)
*** if ( rrproj .gt. rmax ) then
call hfill1(180+idi,eplog,0.d0,wt)
call hfill1(220+idi,timlog,0.d0,wt)
call hfill2(380+idi,eplog,timlog,wt)
*** endif
return
end
c=======================================================================
subroutine pamaf
c-----------------------------------------------------------------------
c pa(rticle) ma(ss) f(illing)
c
c fills particle mass for particle ip in array pama
c resonances and strange baryons included
c particle masses according to geant table,
c taken from the periodic table
c or calculated with the mass formula of weizsaecker
c this subroutine is called from start
c-----------------------------------------------------------------------
implicit double precision (a-h,o-z), integer (i-n)
double precision pama,pi
common /pconst/pama(6000),pi
double precision masses(75),charge(75),amus(59,14)
data masses /
* 0.0 , 0.000511 , 0.000511 , 0.0 , 0.105658 ,
* 0.105658 , 0.134973 , 0.139568 , 0.139568 , 0.497671 ,
* 0.493646 , 0.493646 , 0.939566 , 0.938272 , 0.938272 ,
* 0.497671 , 0.5488 , 1.11563 , 1.18937 , 1.19255 ,
* 1.19743 , 1.3149 , 1.32132 , 1.67243 , 0.939566 ,
* 1.11563 , 1.18937 , 1.19255 , 1.19743 , 1.3149 ,
* 1.32132 , 1.67243 , 1.7841 , 1.7841 , 1.8693 ,
* 1.8693 , 1.8645 , 1.8645 , 1.9693 , 1.9693 ,
* 2.2852 , 80.6 , 80.6 , 91.161 , 1.877 ,
* 2.817 , 3.755 , 0.0 , 0.0 , 0.0 ,
* 0.7669 , 0.7681 , 0.7681 , 1.2309 , 1.2323 ,
* 1.2336 , 1.2349 , 1.2309 , 1.2323 , 1.2336 ,
* 1.2349 , 0.89624 , 0.89209 , 0.89209 , 0.89624 ,
* 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ,
* 0.5488 , 0.5488 , 0.5488 , 0.5488 , 0.0 /
c isotope masses calculated from: atomic data and nucl.data tables 39
c (1988) 289, (wapstra's values, corrected for electron masses)
data ((amus(i,l),i=1,59),l=1,7) /
* 1.8756 , 2.8089 , 57*0. ,
* 2.8083 , 3.7273 , 4.6678 , 5.6054 , 6.5454 , 54*0. ,
* 2*0. , 5.6014 , 6.5337 , 7.4712 , 8.4067 ,
* 9.3471 , 10.2856 , 51*0. ,
* 2*0. , 6.5341 , 7.4547 , 8.3926 , 9.3253 ,
* 10.2644 , 11.2008 , 51*0. ,
* 2*0. , 7.4722 , 8.3932 , 9.3243 , 10.2524 ,
* 11.1886 , 12.1232 , 13.0618 , 13.9986 , 49*0. ,
* 2*0. , 8.4091 , 9.3274 , 10.2538 , 11.1747 , 12.1093 ,
* 13.0406 , 13.9790 , 14.9143 , 15.8531 , 48*0. ,
* 4*0. , 11.1915 , 12.1110 , 13.0400 , 13.9687 , 14.9057 ,
* 15.8394 , 16.7761 , 17.7104 , 47*0. /
data ((amus(i,l),i=1,59),l=8,14) /
* 4*0. , 12.1282 , 13.0446 , 13.9709 , 14.8948 , 15.8302 ,
* 16.7617 , 17.6973 , 18.6293 , 19.5650 , 46*0. ,
* 7*0. , 15.8325 , 16.7629 , 17.6920 , 18.6429 , 19.5564 ,
* 20.4907 , 21.4227 , 22.3587 , 44*0. ,
* 6*0. , 15.8464 , 16.7668 , 17.6947 , 18.6174 , 19.5502 ,
* 20.4794 , 21.4137 , 22.3444 , 23.2839 , 24.2138 , 43*0. ,
* 8*0. , 18.6308 , 19.5532 , 20.4817 , 21.4088 , 22.3414 ,
* 23.2720 , 24.2059 , 25.1387 , 26.0746 , 27.0099 ,
* 27.9469 , 28.8820 , 29.8173 , 30.7546 , 31.6913 , 36*0. ,
* 7*0. , 18.6410 , 19.5658 , 20.4860 , 21.4124 , 22.3354 ,
* 23.2676 , 24.1961 , 25.1292 , 26.0602 , 26.9961 ,
* 27.9291 , 28.8660 , 29.7994 , 30.7376 , 38*0. ,
* 9*0. , 21.4241 , 22.3488 , 23.2714 , 24.1996 , 25.1261 ,
* 26.0579 , 26.9880 , 27.9218 , 28.8541 , 29.7894 ,
* 30.7233 , 31.6599 , 32.5944 , 33.5316 , 36*0. ,
* 9*0. , 22.3591 , 23.2836 , 24.2041 , 25.1304 , 26.0527 ,
* 26.9838 , 27.9128 , 28.8457 , 29.7761 , 30.7111 ,
* 31.6431 , 32.5803 , 33.5128 , 34.4505 , 35.3837 , 35*0. /
c-----------------------------------------------------------------------
c geant particles including rho, k*, and delta
do ip = 1,75
pama(ip) = masses(ip)
enddo
do ip = 76,6000
pama(ip) = 0.
enddo
do ia = 1,59
do ic = 1,ia
in = ia - ic
ip = ia * 100 + ic
if ( ic .le. 14 ) then
c masses from mass table for isotopes
if ( in .eq. 0 ) then
pama(ip) = ic * pama(14)
else
pama(ip) = amus(in,ic)
endif
c simple sum of proton and neutron masses
if ( pama(ip) .eq. 0. )
* pama(ip) = pama(14) * ic + in * pama(13)
else
c weizsaeckers mass formula gives binding energy in mev
b1 = 14.1 * ia
b2 = -13. * ia**(2./3.)
b3 = -0.595 * ic**2 / ia**(1./3.)
b4 = -19. * (ic-in)**2 / ia
b5 = 33.5 / ia**0.75
if ( mod(ic,2) .eq. 0 .and. mod(in,2) .eq. 0 ) then
ss = 1.
elseif ( mod(ic,2) .eq. 1 .and. mod(in,2) .eq. 1 ) then
ss = -1.
else
ss = 0.
endif
bind = (b1 + b2 + b3 + b4 + ss*b5) * 1.e-3
pama(ip) = masses(13) * in + masses(14) * ic - bind
endif
enddo
enddo
c masses of multineutron clusters
do in = 1,59
ip = 100 * in
pama(ip) = in * pama(13)
enddo
return
end
c=======================================================================
subroutine hbook1(ihist,htit,nbin,xlow,xupr,ulog)
implicit double precision (a-h,o-z), integer (i-n)
character chtitle(480)*40,htit*(*)
double precision qhistos(10020,480)
common /histch/ chtitle
common /histos/ qhistos,rpromin,rpromax
+ ,eparmin,eparmax
+ ,xmin348,xmax348,ymin348,ymax348
chtitle(ihist) = ' '
chtitle(ihist) = htit
qhistos(1,ihist) = 1.d0*ihist
qhistos(2,ihist) = 1.d0*nbin
qhistos(3,ihist) = xlow
qhistos(4,ihist) = xupr
qhistos(5,ihist) = ulog
qhistos(6,ihist) = 0.d0
qhistos(7,ihist) = 0.d0
qhistos(8,ihist) = 0.d0
qhistos(9,ihist) = -99.8877d0
do i=10,nbin+10
qhistos(i,ihist) = 0.
enddo
return
end
c=======================================================================
subroutine hbook2(ihist,htit,nbin,xlow,xupr,nbi2,xlo2,xup2,ulog)
implicit double precision (a-h,o-z), integer (i-n)
character chtitle(480)*40,htit*(*)
double precision qhistos(10020,480)
common /histch/ chtitle
common /histos/ qhistos,rpromin,rpromax
+ ,eparmin,eparmax
+ ,xmin348,xmax348,ymin348,ymax348
chtitle(ihist) = ' '
chtitle(ihist) = htit
qhistos(1,ihist) = 1.d0*ihist
qhistos(2,ihist) = 1.d0*nbin
qhistos(3,ihist) = xlow
qhistos(4,ihist) = xupr
qhistos(5,ihist) = ulog
qhistos(6,ihist) = 1.d0*nbi2
qhistos(7,ihist) = xlo2
qhistos(8,ihist) = xup2
qhistos(9,ihist) = -99.8877d0
do i=10,nbin*nbi2+10
qhistos(i,ihist) = 0.
enddo
return
end
c=======================================================================
subroutine hfill1(ihist,xval,yval,wval)
implicit double precision (a-h,o-z), integer (i-n)
character chtitle(480)*40
double precision qhistos(10020,480)
common /histch/ chtitle
common /histos/ qhistos,rpromin,rpromax
+ ,eparmin,eparmax
+ ,xmin348,xmax348,ymin348,ymax348
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
nbin = int(qhistos(2,ihist))
if ( nbin .gt. 0 ) then
if ( 121 .le. ihist .and. ihist .le. 231 ) then
c 128 (1) rec density of p vs log10(r)
c 188 (1) N rec p vs log10(energy)
c 228 (1) N rec p vs log10(time)
if ( xval .le. qhistos(4,ihist) ) then
if ( xval .gt. qhistos(3,ihist) ) then ! regular interval.
ix = (xval - qhistos(3,ihist)) * qhistos(2,ihist) /
/ (qhistos(4,ihist) - qhistos(3,ihist))
if ( ix .ge. 0 ) then
qhistos(11+ix,ihist) = qhistos(11+ix,ihist) + wval
else
qhistos(10,ihist) = qhistos(10,ihist) + wval
endif
else ! underflow.
qhistos(10,ihist) = qhistos(10,ihist) + wval
endif
else ! overflow.
write(*,*) ' ih=',ihist,xval,wval,' overflow'
qhistos(10+nbin+1,ihist) = qhistos(10+nbin+1,ihist) + wval
endif
endif
endif
return
end
c=======================================================================
subroutine hfill2(ihist,xval,yval,wval)
implicit double precision (a-h,o-z), integer (i-n)
character chtitle(480)*40
double precision qhistos(10020,480)
common /histch/ chtitle
common /histos/ qhistos,rpromin,rpromax
+ ,eparmin,eparmax
+ ,xmin348,xmax348,ymin348,ymax348
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if ( ihist .eq. 299 ) then
c particle codes (0,...,29) vs log10(r)
c 299. 30. 0.0000 30.00 2. 100. -4.000 6.000
if ( yval .lt. qhistos(8,ihist) ) then
if ( yval .gt. qhistos(7,ihist) ) then
ix = int(xval)
iy = 1. + (yval - qhistos(7,ihist)) * qhistos(6,ihist) /
/ (qhistos(8,ihist) - qhistos(7,ihist))
iq = ix + (iy-1) * int(qhistos(2,ihist))
qhistos(10+iq,ihist) = qhistos(10+iq,ihist) + wval
else ! y-underflow
iq = 10
qhistos(iq,ihist) = qhistos(iq,ihist) + wval
endif
else ! y-overflow
iq = 10+int(qhistos(2,ihist))*int(qhistos(6,ihist))+1
qhistos(iq,ihist) = qhistos(iq,ihist) + wval
endif
endif
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if ( ( 341 .le. ihist .and. ihist .le. 349 ) .or.
+ ( 361 .le. ihist .and. ihist .le. 369 ) .or.
+ ( 381 .le. ihist .and. ihist .le. 389 ) ) then
c N p vs log10(E) and log10(r)
c 348. 80. -4.000 4.000 3. 100. -4.000 6.000
c N p vs log10(t) and log10(r)
c 368. 100. -4.000 6.000 3. 100. 0.000 5.000
c N p vs log10(E) and log10(t)
c 388. 100. -4.000 6.000 3. 100. 0.000 5.000
if ( yval .lt. qhistos(8,ihist) ) then
if ( yval .gt. qhistos(7,ihist) ) then
if ( xval .lt. qhistos(4,ihist) ) then
if ( xval .ge. qhistos(3,ihist) ) then
ix = 1. + (xval-qhistos(3,ihist)) * qhistos(2,ihist) /
/ (qhistos(4,ihist)-qhistos(3,ihist))
iy = 1. + (yval-qhistos(7,ihist)) * qhistos(6,ihist) /
/ (qhistos(8,ihist)-qhistos(7,ihist))
iq = 10 + ix + (iy-1) * int(qhistos(2,ihist))
qhistos(iq,ihist) = qhistos(iq,ihist) + wval
else ! x-underflow
iq = 10
qhistos(iq,ihist) = qhistos(iq,ihist) + wval
endif
else ! x-overflow
iq = 10+int(qhistos(2,ihist))*int(qhistos(6,ihist))+1
qhistos(iq,ihist) = qhistos(iq,ihist) + wval
endif
else ! y-underflow
iq = 10
qhistos(iq,ihist) = qhistos(iq,ihist) + wval
endif
else ! y-overflow
iq = 10+int(qhistos(2,ihist))*int(qhistos(6,ihist))+1
qhistos(iq,ihist) = qhistos(iq,ihist) + wval
endif
endif
return
end
c=======================================================================
subroutine tabhistos
c-----------------------------------------------------------------------
c Full tabular of histograms including ident number and dimension
c-----------------------------------------------------------------------
c
c 1 (1) longi dist of gammas
c 2 (1) longi dist of e+
c 3 (1) longi dist of e-
c 4 (1) longi dist of mu+
c 5 (1) longi dist of mu-
c 6 (1) longi dist of hadrons
c 7 (1) longi dist of charged
c 8 (1) longi dist of nuclei
c 9 (1) longi dist of cerenk
c
c 11 (1) longi edep of gammas
c 12 (1) longi edep of em ioniz
c 13 (1) longi edep of em cut
c 14 (1) longi edep of mu ioniz
c 15 (1) longi edep of mu cut
c 16 (1) longi edep of had ioniz
c 17 (1) longi edep of had cut
c 18 (1) longi edep of neutrino
c 19 (1) longi edep of sum
c
c 21 (1) longi dist of nucleon
c 22 (1) longi dist of p
c 23 (1) longi dist of n
c 24 (1) longi dist of pi+/-
c 25 (1) longi dist of K+/-
c 26 (1) longi dist of Kl
c 27 (1) longi dist of Ks
c 28 (1) longi dist of Kl/s
c
c 99 (2) particle codes vs log10(r) (here on 299)
c
c 101 (1) part density of gam vs log10(r)
c 102 (1) part density of e+ vs log10(r)
c 103 (1) part density of e- vs log10(r)
c 104 (1) part density of mu+ vs log10(r)
c 105 (1) part density of mu- vs log10(r)
c 106 (1) part density of pi+ vs log10(r)
c 107 (1) part density of pi- vs log10(r)
c 108 (1) part density of p vs log10(r)
c 109 (1) part density of n vs log10(r)
c 110 (1) part density of cerph vs log10(r)
c 111 (1) part density of other vs log10(r)
c
c 121 (1) rec density of gam vs log10(r) "always weight 1"
c 122 (1) rec density of e+ vs log10(r)
c 123 (1) rec density of e- vs log10(r)
c 124 (1) rec density of mu+ vs log10(r)
c 125 (1) rec density of mu- vs log10(r)
c 126 (1) rec density of pi+ vs log10(r)
c 127 (1) rec density of pi- vs log10(r)
c 128 (1) rec density of p vs log10(r)
c 129 (1) rec density of n vs log10(r)
c 130 (1) rec density of cerph vs log10(r)
c 131 (1) rec density of other vs log10(r)
c
c 141 (1) ener density of gam vs log10(r)
c 142 (1) ener density of e+ vs log10(r)
c 143 (1) ener density of e- vs log10(r)
c 144 (1) ener density of mu+ vs log10(r)
c 145 (1) ener density of mu- vs log10(r)
c 146 (1) ener density of pi+ vs log10(r)
c 147 (1) ener density of pi- vs log10(r)
c 148 (1) ener density of p vs log10(r)
c 149 (1) ener density of n vs log10(r)
c 150 (1) ener density of cerph vs log10(r)
c 151 (1) ener density of other vs log10(r)
c
c 161 (1) N gam vs log10(energy)
c 162 (1) N e+ vs log10(energy)
c 163 (1) N e- vs log10(energy)
c 164 (1) N mu+ vs log10(energy)
c 165 (1) N mu- vs log10(energy)
c 166 (1) N pi+ vs log10(energy)
c 167 (1) N pi- vs log10(energy)
c 168 (1) N p vs log10(energy)
c 169 (1) N n vs log10(energy)
c 170 (1) N cerph vs log10(energy)
c 171 (1) N other vs log10(energy)
c
c 181 (1) N rec gam vs log10(energy) "always weight 1"
c 182 (1) N rec e+ vs log10(energy)
c 183 (1) N rec e- vs log10(energy)
c 184 (1) N rec mu+ vs log10(energy)
c 185 (1) N rec mu- vs log10(energy)
c 186 (1) N rec pi+ vs log10(energy)
c 187 (1) N rec pi- vs log10(energy)
c 188 (1) N rec p vs log10(energy)
c 189 (1) N rec n vs log10(energy)
c 190 (1) N rec cerph vs log10(energy)
c 191 (1) N rec other vs log10(energy)
c
c 201 (1) N gam vs log10(time)
c 202 (1) N e+ vs log10(time)
c 203 (1) N e- vs log10(time)
c 204 (1) N mu+ vs log10(time)
c 205 (1) N mu- vs log10(time)
c 206 (1) N pi+ vs log10(time)
c 207 (1) N pi- vs log10(time)
c 208 (1) N p vs log10(time)
c 209 (1) N n vs log10(time)
c 210 (1) N cerph vs log10(time)
c 211 (1) N other vs log10(time)
c
c 221 (1) N rec gam vs log10(time) "always weight 1"
c 222 (1) N rec e+ vs log10(time)
c 223 (1) N rec e- vs log10(time)
c 224 (1) N rec mu+ vs log10(time)
c 225 (1) N rec mu- vs log10(time)
c 226 (1) N rec pi+ vs log10(time)
c 227 (1) N rec pi- vs log10(time)
c 228 (1) N rec p vs log10(time)
c 229 (1) N rec n vs log10(time)
c 230 (1) N rec cerph vs log10(time)
c 231 (1) N rec other vs log10(time)
c
c 241 (1) N gam vs theta
c 242 (1) N e+ vs theta
c 243 (1) N e- vs theta
c 244 (1) N mu+ vs theta
c 245 (1) N mu- vs theta
c 246 (1) N pi+ vs theta
c 247 (1) N pi- vs theta
c 248 (1) N p vs theta
c 249 (1) N n vs theta
c 250 (1) N cerph vs theta
c 251 (1) N other vs theta
c
c 261 (1) N rec gam vs theta "always weight 1"
c 262 (1) N rec e+ vs theta
c 263 (1) N rec e- vs theta
c 264 (1) N rec mu+ vs theta
c 265 (1) N rec mu- vs theta
c 266 (1) N rec pi+ vs theta
c 267 (1) N rec pi- vs theta
c 268 (1) N rec p vs theta
c 269 (1) N rec n vs theta
c 270 (1) N rec cerph vs theta
c 271 (1) N rec other vs theta
c
c 281 (1) N gam vs phi
c 282 (1) N e+ vs phi
c 283 (1) N e- vs phi
c 284 (1) N mu+ vs phi
c 285 (1) N mu- vs phi
c 286 (1) N pi+ vs phi
c 287 (1) N pi- vs phi
c 288 (1) N p vs phi
c 289 (1) N n vs phi
c 290 (1) N cerph vs phi
c 291 (1) N other vs phi
c
c ( 295 (1) ring areas for particle densities )
c
c 299 (2) particle codes vs log10(r)
c
c 301 (1) N rec gam vs phi "always weight 1"
c 302 (1) N rec e+ vs phi
c 303 (1) N rec e- vs phi
c 304 (1) N rec mu+ vs phi
c 305 (1) N rec mu- vs phi
c 306 (1) N rec pi+ vs phi
c 307 (1) N rec pi- vs phi
c 308 (1) N rec p vs phi
c 309 (1) N rec n vs phi
c 310 (1) N rec cerph vs phi
c 311 (1) N rec other vs phi
c
c 321 (1) N gam vs log10(weight)
c 322 (1) N e+ vs log10(weight)
c 323 (1) N e- vs log10(weight)
c 324 (1) N mu+ vs log10(weight)
c 325 (1) N mu- vs log10(weight)
c 326 (1) N pi+ vs log10(weight)
c 327 (1) N pi- vs log10(weight)
c 328 (1) N p vs log10(weight)
c 329 (1) N n vs log10(weight)
c 330 (1) N cerph vs log10(weight)
c 331 (1) N other vs log10(weight)
c
c 341 (2) N gam vs log10(E) and log10(r)
c 342 (2) N e+ vs log10(E) and log10(r)
c 343 (2) N e- vs log10(E) and log10(r)
c 344 (2) N mu+ vs log10(E) and log10(r)
c 345 (2) N mu- vs log10(E) and log10(r)
c 346 (2) N pi+ vs log10(E) and log10(r)
c 347 (2) N pi- vs log10(E) and log10(r)
c 348 (2) N p vs log10(E) and log10(r)
c 349 (2) N n vs log10(E) and log10(r)
c 350 (2) N cerph vs log10(E) and log10(r)
c 351 (2) N other vs log10(E) and log10(r)
c
c 361 (2) N gam vs log10(t) and log10(r)
c 362 (2) N e+ vs log10(t) and log10(r)
c 363 (2) N e- vs log10(t) and log10(r)
c 364 (2) N mu+ vs log10(t) and log10(r)
c 365 (2) N mu- vs log10(t) and log10(r)
c 366 (2) N pi+ vs log10(t) and log10(r)
c 367 (2) N pi- vs log10(t) and log10(r)
c 368 (2) N p vs log10(t) and log10(r)
c 369 (2) N n vs log10(t) and log10(r)
c 370 (2) N cerph vs log10(t) and log10(r)
c 371 (2) N other vs log10(t) and log10(r)
c
c 381 (2) N gam vs log10(E) and log10(t)
c 382 (2) N e+ vs log10(E) and log10(t)
c 383 (2) N e- vs log10(E) and log10(t)
c 384 (2) N mu+ vs log10(E) and log10(t)
c 385 (2) N mu- vs log10(E) and log10(t)
c 386 (2) N pi+ vs log10(E) and log10(t)
c 387 (2) N pi- vs log10(E) and log10(t)
c 388 (2) N p vs log10(E) and log10(t)
c 389 (2) N n vs log10(E) and log10(t)
c 390 (2) N cerph vs log10(E) and log10(t)
c 391 (2) N other vs log10(E) and log10(t)
c
c-----------------------------------------------------------------------
end