278 lines
9.5 KiB
Fortran
278 lines
9.5 KiB
Fortran
c------------------------------------------------------------------------------
|
|
c ~/mc/mcv3k/src/read_mcv3k_out_nt.f
|
|
c------------------------------------------------------------------------------
|
|
c
|
|
program read_mcv3k_out_nt
|
|
c
|
|
c program to read the unformatted output file of the MCV3K_TEST program
|
|
c data are written into a ntuple to use PAW to display data
|
|
c
|
|
c at the moment the data "format" is
|
|
c 12 entries per record
|
|
c
|
|
c 1.entry number of particle I*4
|
|
c 2.entry number of layer I*4
|
|
c 3.-5.entry x,y,z R*4
|
|
c 6.-8.entry px,py,pz R*4
|
|
c 9.entry kinetic energy R*4
|
|
c 10.entry theta R*4
|
|
c 11.entry phi R*4
|
|
c 12.entry total time R*4
|
|
c
|
|
c at the moment particle data at the exit of seven layers are written
|
|
c to the data file (the name must be FOR0##.DAT);
|
|
c
|
|
c TP, PSI, 13-jun-1995
|
|
c
|
|
c TP, 26-jun-1995: modified number of NTUPLE entries for
|
|
c use in analysing of the piE5 data
|
|
c
|
|
c TP, 27-jun-1995: modified again number of ntuple entries
|
|
c only the data in front of the last layer are
|
|
c selected to minimize file space
|
|
c there is a little error now: the first entry
|
|
c of the NTUPLE will be ID=1 and all other variables
|
|
c equal to 0, if the first ID read from file is not
|
|
c 1
|
|
c prompt now for the input file name without extension
|
|
c ( the default extension is .DAT !). the root file
|
|
c name will be used to create the output filename
|
|
c with the extension .NT
|
|
c
|
|
c TP, 10-Aug-1995: read filename without extension from command
|
|
c line using the CLI$GET_VALUE function
|
|
c
|
|
c TP, 30-Jan-1996: read initial and final data from file, so, add data
|
|
c to NTUPLE
|
|
c
|
|
c TP, 30-Jan-1996: changed back
|
|
c TP, 13-Sep-2000: changed to run on Windows NT 4.0 and Unix
|
|
c
|
|
c------------------------------------------------------------------------------
|
|
c
|
|
implicit none ! inquire variable definition
|
|
c
|
|
character*80 fname ! string containing filename without
|
|
integer*4 length ! extension and length
|
|
c
|
|
integer*4 id, il ! particle and layer number
|
|
real*4 x,y,z,px,py,pz,ekin ! location, momentum and energy of
|
|
c ! particle
|
|
real*4 theta, phi,time ! direction angulars and total time of
|
|
c ! flight
|
|
c
|
|
c for the NTUPLE
|
|
c
|
|
c first declarations for the initialisation of the HBOOK package
|
|
c
|
|
integer*4 MAXPAGES
|
|
parameter (MAXPAGES = 1000)
|
|
integer*4 ipawc(128 * MAXPAGES)
|
|
c
|
|
common /pawc/ ipawc
|
|
c
|
|
integer*4 iquest(100)
|
|
common /quest/iquest !needed to get enough space for NTUPLE
|
|
c
|
|
c for handling the NTUPLE file
|
|
c
|
|
integer*4 lunin /20/, lunnt /21/ ,istat
|
|
character*80 infile, ntfile
|
|
c
|
|
c for booking and filling of the NTUPLE
|
|
c
|
|
integer*4 nvar, nprime
|
|
c parameter nvar = 71 !number of variables in NTUPLE
|
|
parameter (nvar = 11) !number of variables in NTUPLE
|
|
parameter (nprime = 1000) !NPRIME words are allocated
|
|
character*8 tags(nvar) !short name of each variable
|
|
real*4 xtup(nvar) !array with the current data
|
|
c !points
|
|
data tags
|
|
1 / 'ID',
|
|
2 'X0', 'Y0', 'Z0', 'PX0', 'PY0', 'PZ0', 'E0', 'TH0', 'PH0', 'T0'/
|
|
c 3 'X1', 'Y1', 'Z1', 'PX1', 'PY1', 'PZ1', 'E1', 'TH1', 'PH1', 'T1'/
|
|
c 5 'X2', 'Y2', 'Z2', 'PX2', 'PY2', 'PZ2', 'E2', 'TH2', 'PH2', 'T2'/
|
|
c 6 'X3', 'Y3', 'Z3', 'PX3', 'PY3', 'PZ3', 'E3', 'TH3', 'PH3', 'T3'/
|
|
c 7 'X4', 'Y4', 'Z4', 'PX4', 'PY4', 'PZ4', 'E4', 'TH4', 'PH4', 'T4',
|
|
c 8 'X5', 'Y5', 'Z5', 'PX5', 'PY5', 'PZ5', 'E5', 'TH5', 'PH5', 'T5',
|
|
c 9 'X6', 'Y6', 'Z6', 'PX6', 'PY6', 'PZ6', 'E6', 'TH6', 'PH6', 'T6'
|
|
c 1 /
|
|
c
|
|
c other variables
|
|
c
|
|
integer*4 temp, n, i
|
|
logical*1 cont
|
|
c
|
|
c-----------------------------------------------------------------------------
|
|
c
|
|
c read first the filename for data input from an MCV3K run and make
|
|
c the NTUPLE filename
|
|
c
|
|
write(6,'(/,$,a)')
|
|
1 'Enter the filename with the MCV3K data without extension : '
|
|
read(5,'(a)') fname
|
|
call str_trim (fname, length)
|
|
c
|
|
c write(6,'('' length of string = '',i)') length
|
|
c
|
|
infile = fname(1:length)//'.dat'
|
|
ntfile = fname(1:length)//'.nt'
|
|
c
|
|
c open input file
|
|
c
|
|
open (lunin, file=infile, status = 'OLD', form = 'UNFORMATTED',
|
|
1 err= 1)
|
|
goto 2
|
|
1 continue
|
|
write(6,'('' %OPEN_MCV3K-E- error opening input file '',a)')
|
|
1 infile
|
|
call exit
|
|
2 continue
|
|
write(6,'('' %OPEN_MCV3K-I- successfully opened input file '',a)')
|
|
1 infile
|
|
c
|
|
c initialize HBOOK
|
|
c
|
|
call hlimit ( 128 * MAXPAGES )
|
|
c
|
|
c open NTUPLE file and book NTUPLE
|
|
c
|
|
iquest(10) = 65000 ! to allocate enough space for the NTUPLE
|
|
c
|
|
call hropen(lunnt,'MCV3K',ntfile,'NQ',4096,istat)
|
|
c
|
|
if (istat.eq.0) then ! everything is all right
|
|
write(6,
|
|
1 '('' %READ_MCV3K-I- succesfully opened NTUPLE file '',a)') ntfile
|
|
else
|
|
write(6,'('' %READ_MCV3K-E- error opening NTUPLE file '',a)')
|
|
1 ntfile
|
|
call exit
|
|
endif
|
|
c
|
|
call hbookn(111,'MCV3KDAT',nvar,'//MCV3K',nprime,tags)
|
|
c
|
|
temp = 1 ! initialize the temporary storage with the
|
|
c ! the first particle number
|
|
id = 1
|
|
n = 0
|
|
c
|
|
c endless loop on reading the file and filling the ntuple
|
|
c
|
|
do while ( .true. )
|
|
c
|
|
xtup(1) = temp
|
|
cont = .true.
|
|
c
|
|
do while ( cont ) ! read until a particle is lost
|
|
c ! then, the array XTUP will be written
|
|
c ! into the NTUPLE, and then XTUP will
|
|
c ! be cleared and the data of the next
|
|
c ! particle will be read
|
|
c
|
|
read(lunin,err=11,end=10) id,il,x,y,z,px,py,pz,ekin,theta,
|
|
1 phi,time
|
|
c write(6,'(I6,2x,I3,10f10.3)') id,il,x,y,z,px,py,pz,
|
|
c 1 ekin,theta,phi,time
|
|
c
|
|
c copy data into the ntuple array xtup
|
|
c
|
|
if ( id .eq. temp) then ! still data from one particle
|
|
c
|
|
xtup( 2 + n*10 ) = x
|
|
xtup( 3 + n*10 ) = y
|
|
xtup( 4 + n*10 ) = z
|
|
xtup( 5 + n*10 ) = px
|
|
xtup( 6 + n*10 ) = py
|
|
xtup( 7 + n*10 ) = pz
|
|
xtup( 8 + n*10 ) = ekin
|
|
xtup( 9 + n*10 ) = theta
|
|
xtup( 10 + n*10 ) = phi
|
|
xtup( 11 + n*10 ) = time
|
|
n = n + 1
|
|
else
|
|
call hfn (111, xtup)
|
|
c
|
|
c reset the xtup array to be ready for the next particle
|
|
c
|
|
do i = 1, nvar
|
|
xtup (i) = 0.
|
|
enddo
|
|
c
|
|
c we just read the initial paramter of the following
|
|
c particle, so, write them to the XTUP array
|
|
c
|
|
xtup( 2 ) = x
|
|
xtup( 3 ) = y
|
|
xtup( 4 ) = z
|
|
xtup( 5 ) = px
|
|
xtup( 6 ) = py
|
|
xtup( 7 ) = pz
|
|
xtup( 8 ) = ekin
|
|
xtup( 9 ) = theta
|
|
xtup( 10) = phi
|
|
xtup( 11) = time
|
|
temp = id
|
|
cont = .false.
|
|
n = 1
|
|
endif
|
|
enddo
|
|
c
|
|
enddo
|
|
10 write(6,'('' Succesfully read data from file '',a)') infile
|
|
close (lunin)
|
|
call hfn (111, xtup) ! save also data from the last particle
|
|
call end_ntuple( lunnt, ntfile )
|
|
call exit
|
|
11 write(6,'('' Error during reading...'')')
|
|
close (lunin)
|
|
call end_ntuple( lunnt, ntfile )
|
|
call exit
|
|
c
|
|
end
|
|
c
|
|
c------------------------------------------------------------------------------
|
|
c
|
|
subroutine end_ntuple(lunnt, ntfile)
|
|
c
|
|
c routine to close the ntuple file
|
|
c
|
|
implicit none
|
|
c
|
|
integer*4 lunnt ! LUN of Ntuple file
|
|
character*(*) ntfile
|
|
integer*4 i
|
|
c
|
|
c------------------------------------------------------------------------------
|
|
c
|
|
call hrout (111, i, ' ')
|
|
call hrend ('MCV3K')
|
|
close (lunnt)
|
|
write(6,'('' %READ_MCV3K-I- closed NTUPLE file '',a,/)') ntfile
|
|
return
|
|
end
|
|
c
|
|
c------------------------------------------------------------------------------
|
|
c
|
|
subroutine str_trim(string, length)
|
|
c
|
|
c routine to get the length of string without blanks
|
|
c
|
|
implicit none
|
|
c
|
|
character*(*) string
|
|
integer*4 i, length
|
|
c
|
|
i = 1
|
|
do while (string(i:i) .ne. ' ')
|
|
i = i + 1
|
|
enddo
|
|
length = i-1
|
|
return
|
|
end
|
|
c
|
|
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
c
|
|
c EOF READ_MCV3K_OUT.FOR
|