Skip to content
Snippets Groups Projects
Commit a4e9f1e8 authored by Hans Rabus's avatar Hans Rabus
Browse files

Upload New File

parent 8a6ccceb
No related branches found
No related tags found
1 merge request!1Add new directory
Hide whitespace changes
Inline Side-by-side
Simulation_Track_Structure/PTra_HIT/Obsolete/PTra_c3h8_20170608HIT_subs_PSF0Carbon.f 0 → 100644
+ 813
− 0
View file @ a4e9f1e8
************************************************************************
SUBROUTINE HSCORE(ICALL)
************************************************************************
IMPLICIT DOUBLE PRECISION (A-H)
IMPLICIT INTEGER*4 (I-N)
IMPLICIT DOUBLE PRECISION (O-Z)
C
PARAMETER(NETRPM=1000000)
DIMENSION NEVID(NETRPM),XION(NETRPM),
& YION(NETRPM),ZION(NETRPM),
& MAPROJ(NETRPM),NZATOM(NETRPM),EPROJ(NETRPM)
COMMON /SPUREN/ XION,YION,ZION,EPROJ,NEVID,MAPROJ,NZATOM,NETRP ! HR@19-Mar-2023
PARAMETER(MTARGY=50,MTARGZ=100)
DIMENSION YTARG(-MTARGY:MTARGY),ZTARG(MTARGZ)
COMMON /TARPOS/ DYTARG,DZTARG,YTARG,ZTARG,NTARGY,NTARGZ ! HR@05-Aug-2023
COMMON /TGSIZE/ RTARG
COMMON /DETVOL/ XTMIN,XTMAX,YTMIN,YTMAX,ZTMIN,ZTMAX ! HR@05-Aug-2023
C! /LSCALE/ DCF conversion factor from length in mm to mass per area
C! in µg/nm**2.
COMMON /LSCALE/ DCF
COMMON /INTVOL/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX ! dimensions in mm
PARAMETER(MAXNU=40)
DIMENSION PNU(-MTARGY:MTARGY,MTARGZ,0:MAXNU), ! P_nu averaged over events
& TM1(-MTARGY:MTARGY,MTARGZ) ! M1 averaged over events
COMMON /TSCORE/ PNU,TM1
COMMON /HISTOR/ MXG4EV,NEVTS,NIONS,NNULL
DIMENSION APNU(-MTARGY:MTARGY,MTARGZ,0:MAXNU), ! P_nu of current event in detectors above
& BPNU(-MTARGY:MTARGY,MTARGZ,0:MAXNU), ! P_nu of current event in detectors below
& AM1(-MTARGY:MTARGY,MTARGZ),BM1(-MTARGY:MTARGY,MTARGZ), ! M1 of current event in detectors above and below
& NUMAXA(-MTARGY:MTARGY,MTARGZ), ! Max. cluster size of current event in detectors above
& NUMAXB(-MTARGY:MTARGY,MTARGZ) ! Max. cluster size of current event in detectors below
COMMON /LSCORE/ APNU,BPNU,AM1,BM1,NUMAXA,NUMAXB,LEVENT
PARAMETER(NM1=(2*MTARGY+1)*MTARGZ,NPNU=NM1*MAXNU)
DATA AM1,BM1,TM1 / NM1*0.0d0, NM1*0.0d0, NM1*0.0d0 /
DATA APNU,BPNU / NM1*1.0d0,NPNU*0.0d0,NM1*1.0d0,NPNU*0.0d0 /
DATA PNU / NM1*0.0d0,NPNU*0.0d0 /
C! /LUNMBS/ Logical unit numbers of input and output files and their names
CHARACTER DUMPFN*40,OUTFNM*40,PREFIX*20
COMMON /LUNMBS/ LUNINP(3),LUNOUT,LUNRES,DUMPFN,OUTFNM,PREFIX
C -----------------------------------------------------------------
C WRITE(*,'(2(a,i8),a,i10,a,i8,a,i2,a)') ' HSCORE | MXG4EV=',
C & MXG4EV,' | NEVTS=',NEVTS,' | NIONS=',NIONS,
C & ' | NNULL=',NNULL, ' | ICALL=',ICALL, ' | ENTRY'
LEVENT=NEVID(1)
LPOINT=NETRP
IF(ICALL.EQ.0) THEN ! File was not yet completely read --> ignore last event which may be only partially read
DO WHILE(NEVID(LPOINT).EQ.NEVID(NETRP).AND.LPOINT.GT.0)
LPOINT=LPOINT-1
END DO
END IF
!PRINT*,LPOINT,NEVID(1),NEVID(LPOINT),NEVID(NETRP)
DO ITY=-NTARGY,NTARGY
DO ITZ=1,NTARGZ
AM1(ITY,ITZ)=0.0d0
BM1(ITY,ITZ)=0.0d0
APNU(ITY,ITZ,0)=1.0d0
BPNU(ITY,ITZ,0)=1.0d0
DO NU=1,MAXNU
APNU(ITY,ITZ,NU)=0.0d0
BPNU(ITY,ITZ,NU)=0.0d0
END DO !NU=1,MAXNU
END DO ! ITZ=1,NTARGZ
END DO ! ITY=1,NTARGY
Process_ETP_stack: DO I=1,NETRP
!Print*,'@Loop',LEVENT,I,J,NEVID(I)
C! Scoring
NewEvent: IF(NEVID(I).NE.LEVENT) THEN ! Transfer info to global counters
!Print*,'@',NEVID(I),LEVENT
DO ITY=-NTARGY,NTARGY
!WRITE(*,'(a,i2,$)') 'ITY',ITY
DO ITZ=1,NTARGZ
TM1(ITY,ITZ)=TM1(ITY,ITZ)+0.5d0*
& (AM1(ITY,ITZ)+BM1(ITY,ITZ))
AM1(ITY,ITZ)=0.0d0
BM1(ITY,ITZ)=0.0d0
NUMAX=MAX(NUMAXA(ITY,ITZ),NUMAXB(ITY,ITZ))
DO NU=0,NUMAX
PNU(ITY,ITZ,NU)=PNU(ITY,ITZ,NU)+0.5d0*
& (APNU(ITY,ITZ,NU)+BPNU(ITY,ITZ,NU))
APNU(ITY,ITZ,NU)=0.0d0
BPNU(ITY,ITZ,NU)=0.0d0
END DO !NU=0,NUMAX
APNU(ITY,ITZ,0)=1.0d0
NUMAXA(ITY,ITZ)=0
BPNU(ITY,ITZ,0)=1.0d0
NUMAXB(ITY,ITZ)=0
END DO ! ITZ=1,NTARGZ
END DO ! ITY=1,NTARGY
LEVENT=NEVID(I) ! store new event ID
ENDIF NewEvent
IF(I.GT.LPOINT) EXIT Process_ETP_stack
C* WRITE(*,'(3f10.2)') XTMIN,YTMIN,ZTMIN
C* WRITE(*,'(3f10.2)') XION(I),YION(I),ZION(I)
C* WRITE(*,'(3f10.2)') XTMAX,YTMAX,ZTMAX
C! If the energy transfer point is outside the target region, then take next one
IF(XION(I).LT.XTMIN.OR.XION(I).GT.XTMAX.OR.
& ABS(YION(I)).GT.YTMAX.OR.
& ZION(I).LT.ZTMIN.OR.ZION(I).GT.ZTMAX) THEN
!PRINT*,XION(I)
CYCLE
ELSE
NIONS=NIONS+1
ENDIF
C! Get the first and last indices of relevant targets along y and z
ITYMIN=NINT((YION(I)-RTARG)/DYTARG)-1
IF(ITYMIN.LT.-MTARGY) ITYMIN=-MTARGY
ITYMAX=NINT((YION(I)+RTARG)/DYTARG)+1
IF(ITYMAX.GT.MTARGY) ITYMAX=MTARGY
ITZMIN=INT((ZION(I)-ZTARG(1)-RTARG)/DZTARG)-1
IF(ITZMIN.LT.1) ITZMIN=1
ITZMAX=INT((ZION(I)-ZTARG(1)+RTARG)/DZTARG)+1
IF(ITZMAX.GT.MTARGZ) ITZMAX=MTARGZ
DO ITY=ITYMIN,ITYMAX ! Loop through all revelant targets
DO ITZ=ITZMIN,ITZMAX
Y=YION(I)-YTARG(ITY)
Z=ZION(I)-ZTARG(ITZ)
R2=Y*Y+Z*Z
TPROB=CEFFIC(XION(I),R2) ! Above detector
IF(TPROB.GT.0) THEN
NUMAXA(ITY,ITZ)=NUMAXA(ITY,ITZ)+1
IF(NUMAXA(ITY,ITZ).GT.MAXNU) NUMAXA(ITY,ITZ)=MAXNU
NUMAX=NUMAXA(ITY,ITZ)
AM1(ITY,ITZ)=AM1(ITY,ITZ)+TPROB
IF(NUMAX.EQ.1) THEN
APNU(ITY,ITZ,0)=1.0d0-TPROB
APNU(ITY,ITZ,1)=TPROB
ELSE
IF(NUMAX.EQ.MAXNU) THEN
APNU(ITY,ITZ,NUMAX)=APNU(ITY,ITZ,NUMAX)+
& APNU(ITY,ITZ,NUMAX-1)*TPROB
ELSE
APNU(ITY,ITZ,NUMAX)=APNU(ITY,ITZ,NUMAX-1)*TPROB
ENDIF
DO NU=NUMAX-1,1,-1 ! was NUMAX,1,-1 ! 7-NOV-2023
APNU(ITY,ITZ,NU)=APNU(ITY,ITZ,NU)*(1.0d0-TPROB)+
& APNU(ITY,ITZ,NU-1)*TPROB
END DO
APNU(ITY,ITZ,0)=APNU(ITY,ITZ,0)*(1.0d0-TPROB)
ENDIF
END IF
TPROB=CEFFIC(-XION(I),R2) ! Below detector
IF(TPROB.GT.0) THEN
NUMAXB(ITY,ITZ)=NUMAXB(ITY,ITZ)+1
IF(NUMAXB(ITY,ITZ).GT.MAXNU) NUMAXB(ITY,ITZ)=MAXNU
NUMAX=NUMAXB(ITY,ITZ)
BM1(ITY,ITZ)=BM1(ITY,ITZ)+TPROB
IF(NUMAX.EQ.1) THEN
BPNU(ITY,ITZ,0)=1.0d0-TPROB
BPNU(ITY,ITZ,1)=TPROB
ELSE
IF(NUMAX.EQ.MAXNU) THEN
BPNU(ITY,ITZ,NUMAX)=BPNU(ITY,ITZ,NUMAX)+
& BPNU(ITY,ITZ,NUMAX-1)*TPROB
ELSE
BPNU(ITY,ITZ,NUMAX)=BPNU(ITY,ITZ,NUMAX-1)*TPROB
ENDIF
DO NU=NUMAX-1,1,-1 ! was NUMAX,1,-1 ! 7-NOV-2023
BPNU(ITY,ITZ,NU)=BPNU(ITY,ITZ,NU)*(1.0d0-TPROB)+
& BPNU(ITY,ITZ,NU-1)*TPROB
END DO
BPNU(ITY,ITZ,0)=BPNU(ITY,ITZ,0)*(1.0d0-TPROB)
ENDIF
END IF
END DO
END DO
C! Save the ionization to the list file
IF(LUNOUT.GT.0) WRITE(LUNOUT,'(I9,3f9.3,2i4,f12.3)') NEVID(I),
& XION(I),YION(I),ZION(I),NZATOM(I),MAPROJ(I),EPROJ(I)
C*! & XION(I)/DCF,YION(I)/DCF,ZION(I)/DCF ! convert back to mm
C! zero the counters
NEVID(I)=-1
XION(I)=0.D+00
YION(I)=0.D+00
ZION(I)=0.D+00
EPROJ(I)=0.D+00
MAPROJ(I)=0
NZATOM(I)=0
END DO Process_ETP_stack
NETRP=NETRP-LPOINT ! remaining energy transfer points
DO I=1,NETRP
NEVID(I)=NEVID(LPOINT+I)
XION(I)=XION(LPOINT+I)
YION(I)=YION(LPOINT+I)
ZION(I)=ZION(LPOINT+I)
EPROJ(I)=EPROJ(LPOINT+I)
MAPROJ(I)=MAPROJ(LPOINT+I)
NZATOM(I)=NZATOM(LPOINT+I)
NEVID(LPOINT+I)=-1
XION(LPOINT+I)=0.D+00
YION(LPOINT+I)=0.D+00
ZION(LPOINT+I)=0.D+00
EPROJ(LPOINT+I)=0.D+00
MAPROJ(LPOINT+I)=0
NZATOM(LPOINT+I)=0
END DO
WRITE(*,'(2(a,i8),a,i10,a,i8,a,i2,a)') ' HSCORE | MXG4EV=',
& MXG4EV,' | NEVTS=',NEVTS,' | NIONS=',NIONS,
& ' | NNULL=',NNULL, ' | ICALL=',ICALL, ' | EXIT'
CALL WRTRES(ICALL)
END ! HSCORE
************************************************************************
SUBROUTINE ISTACK(IEV,X1,Y1,Z1,E,NZPRJ,MAPRJ)
************************************************************************
C! Hans Rabus, 03-JUN-2023
C! This procedure handles the scoring of ionizations into the Arrays
C! of COMMON Block .
C
IMPLICIT DOUBLE PRECISION (A-H)
IMPLICIT INTEGER*4 (I-N)
IMPLICIT DOUBLE PRECISION (O-Z)
C
PARAMETER(NETRPM=1000000)
DIMENSION NEVID(NETRPM),XION(NETRPM),
& YION(NETRPM),ZION(NETRPM),
& MAPROJ(NETRPM),NZATOM(NETRPM),EPROJ(NETRPM)
COMMON /SPUREN/ XION,YION,ZION,EPROJ,NEVID,MAPROJ,NZATOM,NETRP ! HR@19-Mar-2023
C! /LSCALE/ DCF conversion factor from length in mm to mass per area
C! in µg/nm**2.
COMMON /LSCALE/ DCF
COMMON /INTVOL/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX ! dimensions in mm
IF(IEV.EQ.0) PRINT*, 'ISTACK',X1/DCF,Y1/DCF,Z1/DCF,IEV
C! Score the ionization
IF((X1.GT.XMIN).AND.(X1.LT.XMAX).AND.(Y1.GT.YMIN).AND.(Y1.LT.YMAX)
& .AND.(Z1.GT.ZMIN).AND.(Z1.LT.ZMAX)) THEN
IF(NETRP.EQ.NETRPM) THEN
CALL HSCORE(0)
ELSE
NETRP=NETRP+1
ENDIF
NEVID(NETRP)=IEV
XION(NETRP)=X1
YION(NETRP)=Y1
ZION(NETRP)=Z1
EPROJ(NETRP)=E
MAPROJ(NETRP)=MAPRJ
NZATOM(NETRP)=NZPRJ
ENDIF
RETURN
END
************************************************************************
SUBROUTINE OPENPSF(IRFILE) ! Version 23-SEP-2023
************************************************************************
C! Hans Rabus, 23-SEP-2023
C! Modified the routine to account for the new approach of using
C! particles passing the plane perpendicular to the intiial C12-ion
C! trajectory that contains the center of the extraction aperture
C
C! Hans Rabus, 03-JUN-2023
C! This procedure handles the files containing data from Miriam's
C! simulations of carbon ion beam experiments with the Ion Counter
C! nanodosimeter at HIT.
C! The simulations were run in parallel threads, which each produced
C! three files with different information, which are flagged by the
C! variable IRFILE (from 1 to 3):
C! 1. Event ID and x- and y- position of the carbon ion when crossing
C! the plane through the center of the extraction aperture that is
C! perpendicular to the nominal beam direction.
C! 2. A list of records consisting of event ID, coordinates of an
C! energy transfer point where an ionization occurred and the
C! deposited energy for all such energy transfer points that are
C! within a plane-parellel slab which is symmetrically bi-sected
C! by the plane mentioned in point 1.
C! 3. A list of records consisting of the event ID and the phase-
C! space coordinates (position, energy, direction of motion) of
C! electrons stopped inside the gas volume of the nanodosimeter
C! when their energy drops below 10 keV.
C! For the initial value of IRFILE of 0, this routine prompts for and
C! reads the number of threads in the simulations as well as the name
C! of a file listing the files produced in the simulations.
C! Otherwise it closes the files from the previous call and opens
C! those from the present pass.
C
IMPLICIT DOUBLE PRECISION (A-H)
IMPLICIT INTEGER*4 (I-N)
IMPLICIT DOUBLE PRECISION (O-Z)
C
CHARACTER*80 HEADER
LOGICAL CLOSIT
C
PARAMETER(MBATCH=32)
CHARACTER*32 FLIST(MBATCH,3) ! List of simulation data files
C! /FILES/ File names, number of threads, current thread number
COMMON /FILES/ FLIST,IFILE,NFILES,NFIRST
C! /LUNMBS/ Logical unit numbers of input and output files and their names
CHARACTER DUMPFN*40,OUTFNM*40,PREFIX*20
COMMON /LUNMBS/ LUNINP(3),LUNOUT,LUNRES,DUMPFN,OUTFNM,PREFIX
C! /FLAGS/ Flags indicating that end of file was not yet reached
LOGICAL ISOPEN(3) ! Flags whether files are open
COMMON /FLAGS/ ISOPEN
C! /EVENTS/ IEVC,IEVI, IEVEL are the event IDs associated with the
C! three groups of data.
COMMON /EVENTS/ IEVC,IEVI,IEVEL,IEVBAD
C! Open the next three files
PRINT*, 'OPENPSF',IFILE,NFILES
IF(IFILE.LE.NFILES) THEN
!PRINT*
DO I=1,1 !3 !23-SEP-2023
OPEN(LUNINP(I),FILE=FLIST(IFILE,I),STATUS='OLD',ERR=30)
PRINT*, 'Opened '//FLIST(IFILE,I)
C! Read header to avoid crash
NHEADL=0
READ(LUNINP(I),'(a)') HEADER
DO WHILE(HEADER(1:1).EQ.'#'.OR.HEADER(1:1).EQ.'e')
NHEADL=NHEADL+1
READ(LUNINP(I),'(a)') HEADER
END DO
REWIND(LUNINP(I))
DO J=1,NHEADL
READ(LUNINP(I),*) HEADER
END DO
ISOPEN(I)=.TRUE.
END DO
IRFILE=1
IEVBAD=-1
IEVC=-1
IEVEL=-1
IEVI=-1
IF(LUNOUT.GT.0) THEN
INQUIRE(UNIT=LUNOUT,OPENED=CLOSIT)
IF(CLOSIT) CLOSE(LUNOUT)
OPEN(LUNOUT,FILE=TRIM(PREFIX)//'_'//FLIST(IFILE,1),
& STATUS='UNKNOWN',ERR=40)
PRINT*, 'Opened '//TRIM(PREFIX)//'_'//FLIST(IFILE,1)
ENDIF
ELSE
PRINT*, 'Batch number exceeds number of batches.'
STOP
ENDIF
RETURN
30 PRINT*, 'ERROR: File '//FLIST(IFILE,I)//' could not be opened'
STOP
40 PRINT*, 'ERROR: File '//TRIM(PREFIX)//'_'//FLIST(IFILE,I)//
& ' could not be opened'
STOP
END
************************************************************************
SUBROUTINE PSFLIST() ! Version 23-SEP-2023
************************************************************************
C! Hans Rabus, 23-SEP-2023
C! Modified the routine to account for the new approach of using
C! particles passing the plane perpendicular to the intiial C12-ion
C! trajectory that contains the center of the extraction aperture
C
C! Hans Rabus, 07-AUG-2023
C! This procedure initializes the COMMON block /FILES/ with a list
C! of output files from Miriam's simulations of carbon ion beam
C! experiments with the Ion Counter nanodosimeter at HIT.
C! The user is equested to input the number of files to process
C! and three name patterns (where a # represents the file number).
C
IMPLICIT DOUBLE PRECISION (A-H)
IMPLICIT INTEGER*4 (I-N)
IMPLICIT DOUBLE PRECISION (O-Z)
C
CHARACTER*80 HEADER
CHARACTER TRIM,ADJUSTL
C
PARAMETER(MBATCH=32)
CHARACTER*32 FLIST(MBATCH,3) ! List of simulation data files
CHARACTER*32 FPATT ! File name pattern
CHARACTER*8 STRING ! Receives the file number
C! /FILES/ File names, number of threads, current thread number
COMMON /FILES/ FLIST,IFILE,NFILES,NFIRST
C! Get number of files
PRINT*, 'Enter first and last file numbers:'
READ(*,*) NFIRST,NLAST
NFILES=NLAST-NFIRST+1
PRINT*, NFILES
IF(NFILES.GT.MBATCH) THEN
PRINT*, 'ERROR: Number of batches exceeds maximum of ',MBATCH
PRINT*, 'Edit code and change parameter MBATCH.'
STOP
ENDIF
PRINT*, 'Pattern of the data file names in the form AAAA#.BBB'
PRINT*, 'Particles passing the extraction aperture plane ' !23-SEP-2023
C!23-SEP-2023 PRINT*, 'Pattern 1: energy deposits in silicon detectors'
C!23-SEP-2023 PRINT*, 'Pattern 2: loci of ionizations in gas volume'
C!23-SEP-2023 PRINT*, 'Pattern 3 phase space data of electrons in gas volume'
C!23-SEP-2023 DO J=1,3
J=1 !23-SEP-2023
READ(*,'(a)') FPATT
LOC=INDEX(FPATT,'#')
DO I=1,NFILES
WRITE(STRING,'(I8)') I-1+NFIRST
FLIST(I,J)=FPATT(1:LOC-1)//TRIM(ADJUSTL(STRING))
& //FPATT(LOC+1:80)
END DO
C!23-SEP-2023 END DO
PRINT*, 'Creating list of file names was successful.'
END
************************************************************************
SUBROUTINE READPZC(IRFILE,APROJ,ZATOM)
************************************************************************
C! Hans Rabus, 23-SEP-2023
C! This procedure reads the data from Miriam's simulations of carbon
C! ion beam experiments with the Ion Counter nanodosimeter at HIT.
C! The simulations produces files with the following information:
C! X1,Y1,Z1: Last position of the particle before it passes the plane
C! through the center of the extraction aperture that is
C! perpendicular to the nominal beam direction.
C! X2,Y2,Z2: First position of the particle after it passes the plane
C! through the center of the extraction aperture that is
C! perpendicular to the nominal beam direction.
C! EP: Energy of the particle at point (X1,Y1,Z1)
C! PARTIC: String describing the particle type.
C! NEVENT: Number of the event.
C! The information read from the files is passed to the main program
C! through the common block EVINFO. The flag IRFILE indicates that
C! data was read (IRFILE=1) or the end of file was reached (IRFILE=4)
C
IMPLICIT DOUBLE PRECISION (A-H)
IMPLICIT INTEGER*4 (I-N)
IMPLICIT DOUBLE PRECISION (O-Z)
C
CHARACTER*40 PARTIC ! particle type
C
PARAMETER(MBATCH=32)
CHARACTER*32 FLIST(MBATCH,3) ! List of simulation data files
C! /FILES/ File names, number of threads, current thread number
COMMON /FILES/ FLIST,IFILE,NFILES,NFIRST
C
C! /LUNMBS/ Logical unit numbers of input and output files and their names
CHARACTER DUMPFN*40,OUTFNM*40,PREFIX*20
COMMON /LUNMBS/ LUNINP(3),LUNOUT,LUNRES,DUMPFN,OUTFNM,PREFIX
C! /FLAGS/ Flags indicating that end of file was not yet reached
LOGICAL ISOPEN(3) ! Flags whether files are open
COMMON /FLAGS/ ISOPEN
C
C! /EVINFO/ YC,COSTHC,SINTHC y-coordinate and durection cosines of
C! original C ion trajectory when passing plane through
C! extraction aperture. (Are changed to 0, 1, 0.)
C! X1P...EP : start position and start energy of particle,
C! UP,VP,WP: direction cosines of momentum,
C! APROJ,ZATOM: mass number and charge number of particle.
COMMON /EVINFO/ YC,COSTHC,SINTHC,X,Y,Z,E,U,V,W
C! /EVENTS/ IEVC,IEVI, IEVEL are the event IDs associated with the
C! three groups of data.
COMMON /EVENTS/ IEVC,IEVI,IEVEL,IEVBAD
C! /LSCALE/ DCF conversion factor from length in mm to mass per area
C! in µg/nm**2.
COMMON /LSCALE/ DCF
COMMON /INTVOL/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX ! dimensions in mm
COMMON /EXTRAP/ XEXTAP,ZEXTAP
COMMON /HISTOR/ MXG4EV,NEVTS,NIONS,NNULL
LOGICAL NEWEV1, NEWEV2, NODATA ! Flags for new events and no data read (yet)
COMMON /NEWST/ NEWEV1, NEWEV2, NODATA
LOGICAL LDEBUG(4)
COMMON /PSDEBG/ LDEBUG
PARAMETER(NETRPM=1000000)
DIMENSION NEVID(NETRPM),XION(NETRPM),YION(NETRPM),ZION(NETRPM),
& MAPROJ(NETRPM),NZATOM(NETRPM),EPROJ(NETRPM)
COMMON /SPUREN/ XION,YION,ZION,EPROJ,NEVID,MAPROJ,NZATOM,NETRP ! HR@19-Mar-2023
DATA IEVI, IEVEL, NEVTS, NIONS, NNULL / -1, -1, 0, 0, 0 /
C ------------------------------------------------------------------
C Code starts here
C ------------------------------------------------------------------
NODATA=.TRUE.
C! The following while loop is meant to make sure that events without
C! ionizations and electrons are properly handledmore
IF(LDEBUG(4).AND.IRFILE.LT.4) PRINT*, IRFILE,FLIST(IFILE, IRFILE)
C! Read record from input file
C! 13-NOV-2023 Ensuing changes are to omit data from tracks with C energy below ionization threshold
NEV=IEVBAD ! 13-NOV-2023
PARTIC='C?' ! 27-JAN-2024
C!27-JAN-2024 DO WHILE(NEV.EQ.IEVBAD) ! 13-NOV-2023
DO WHILE(NEV.EQ.IEVBAD.OR.PARTIC(1:1).EQ.'C') ! 27-JAN-2024
READ(LUNINP(1),*,END=10,ERR=40) X2,Y2,Z2,X1,Y1,Z1,E,PARTIC,NEV
END DO ! 13-NOV-2023
IF(NEV.NE.IEVC) THEN ! 27-JAN-2024 >>>>>>>>>>>>>>>>>>>>>>>
IF(NEVID(NETRP).EQ.IEVC) THEN
NEVTS=NEVTS+1
ELSE
NNULL=NNULL+1
ENDIF
IEVC=NEV ! event number
IF(MXG4EV.LT.IEVC) MXG4EV=IEVC ! Score largest event ID
ENDIF ! 27-JAN-2024 <<<<<<<<<<<<<<<<<<<<<<<
IF(LDEBUG(1)) WRITE(*,'(a,2i8,1X,a,1X,6f8.2,f8.0)') 'READPZC 615',
& NEVTS,NEV,TRIM(PARTIC),X2,Y2,Z2,X1,Y1,Z1,E
C
C! Calculate direction information
DX=X2-X1
DY=Y2-Y1
DZ=Z2-Z1
DIST=SQRT(DX*DX+DY*DY+DZ*DZ)
UP=DX/DIST
VP=DY/DIST
WP=DZ/DIST
IF(LDEBUG(2)) WRITE(*,'(a,i8,1X,a,1X,6f8.2,f8.0)') 'READPZC 625',
& NEV,TRIM(PARTIC),DX,DY,DZ,DIST,UP,VP,WP
C!27-JAN-2024 IF(PARTIC(1:2).EQ.'C1'.OR.PARTIC(1:2).EQ.'C9') THEN
C!27-JAN-2024 ZATOM=6.
C!27-JAN-2024 IF(PARTIC(2:2).EQ.'9') THEN
C!27-JAN-2024 APROJ=9.
C!27-JAN-2024 ELSE
C!27-JAN-2024 READ(PARTIC(2:3),*) APROJ
C!27-JAN-2024 ENDIF
C!27-JAN-2024 IEVC=NEV ! event number
C!27-JAN-2024 VNORM=1./SQRT(VP*VP+WP*WP)
C!27-JAN-2024 COSTHC=WP*VNORM
C!27-JAN-2024 SINTHC=VP*VNORM
C!27-JAN-2024 YC=COSTHC*Y1-SINTHC*(Z1-ZEXTAP) ! impact parameter in target plane
C!27-JAN-2024 NEVTS=NEVTS+1
C!27-JAN-2024 IF(MXG4EV.LT.IEVC) MXG4EV=IEVC ! Score largest event ID
C!27-JAN-2024 ELSE IF(PARTIC(1:5).EQ.'proto') THEN
IF(PARTIC(1:5).EQ.'proto') THEN !27-JAN-2024
ZATOM=1.
APROJ=1.
ELSE IF(PARTIC(1:5).EQ.'deute') THEN
ZATOM=1.
APROJ=2.
ELSE IF(PARTIC(1:5).EQ.'trito') THEN
ZATOM=1.
APROJ=3.
ELSE IF(PARTIC(1:5).EQ.'alpha') THEN
ZATOM=2.
APROJ=4.
ELSE IF(PARTIC(1:2).EQ.'He') THEN
ZATOM=2.
APROJ=3.
ELSE IF(PARTIC(1:2).EQ.'Li') THEN
ZATOM=3.
READ(PARTIC(3:3),*) NAPROJ
IF(NAPROJ.EQ.1) READ(PARTIC(3:4),*) NAPROJ
APROJ=REAL(NAPROJ)
ELSE IF(PARTIC(1:2).EQ.'Be') THEN
ZATOM=4.
READ(PARTIC(3:3),*) NAPROJ
IF(NAPROJ.EQ.1) READ(PARTIC(3:4),*) NAPROJ
APROJ=REAL(NAPROJ)
ELSE IF(PARTIC(1:1).EQ.'B') THEN
ZATOM=5.
READ(PARTIC(2:2),*) NAPROJ
IF(NAPROJ.EQ.1) READ(PARTIC(2:3),*) NAPROJ
APROJ=REAL(NAPROJ)
ELSE IF(PARTIC(1:2).EQ.'Na') THEN
C! ZATOM=11.
C! READ(PARTIC(3:4),*) APROJ
Print*,'Particle ignored: '//PARTIC
ELSE IF(PARTIC(1:2).EQ.'Ne') THEN
C! ZATOM=10.
C! READ(PARTIC(3:4),*) APROJ
Print*,'Particle ignored: '//PARTIC
ELSE IF(PARTIC(1:1).EQ.'N') THEN
ZATOM=7.
READ(PARTIC(2:2),*) NAPROJ
IF(NAPROJ.EQ.1) READ(PARTIC(2:3),*) NAPROJ
APROJ=REAL(NAPROJ)
ELSE IF(PARTIC(1:2).EQ.'e-') THEN
ZATOM=-1.
APROJ=5.50000E-04
ELSE
Print*,'WARNING: Unknown particle ignored: '//PARTIC
ENDIF
U=UP
C!27-JAN-2024 V=COSTHC*VP-SINTHC*WP ! rotate
C!27-JAN-2024 W=SINTHC*VP+COSTHC*WP ! rotate
V=VP ! !27-JAN-2024
W=WP ! !27-JAN-2024
X=X1 ! convert mm to µg/nm**2
C!27-JAN-2024 Y=(COSTHC*Y1-SINTHC*(Z1-ZEXTAP)-YC) ! convert to distance from C ion
C!27-JAN-2024 Z=(SINTHC*Y1+COSTHC*(Z1-ZEXTAP)+ZEXTAP) ! rotate
Y=Y1 ! 27-JAN-2024
Z=Z1 ! 27-JAN-2024
IF(APROJ.LT.0.5) THEN
E=E*1.0d6 ! convert MeV to eV for electrons
ELSE
E=E*1.0d3 ! convert MeV to keV for ions
C! start ion from begin of interaction volume
X=X+(ZMIN-Z)*U/W
Y=Y+(ZMIN-Z)*V/W
Z=ZMIN ! Bugfix 13-JAN-2024 - was missing before
ENDIF
IF(LDEBUG(3)) WRITE(*,'(a,i8,1X,a,1X,6f8.2,f8.0)') 'READPZC 681',
& NEV,TRIM(PARTIC),X,Y,Z,U,V,W,E
!!X=X! *DCF ! Omit these unneccesary conversions ! convert mm to µg/nm**2
!!Y=Y! *DCF ! Omit these unneccesary conversions ! convert to distance from C ion and from mm to µg/nm**2
!!Z=Z! *DCF ! Omit these unneccesary conversions ! rotate and convert mm to µg/nm**2
NODATA=.FALSE.
10 IF(NODATA) THEN ! End of file was reached
CLOSE(LUNINP(1))
ISOPEN(1)=.FALSE.
IRFILE=4 ! Flags all data have been read
ENDIF
RETURN
40 PRINT*,'ERROR reading from file '//FLIST(IFILE,1)
STOP
END
************************************************************************
CHARACTER*24 FUNCTION TSTAMP()
************************************************************************
CHARACTER DATE*8, TIME*10, ZONE*5, TIMEST*24
INTEGER*4 IDT(8)
CALL DATE_AND_TIME(DATE,TIME,ZONE,IDT)
C! 123456789 123456789 1234
TIMEST = 'DD-MMM-YYYY HH:MM +HH:MM'
TSTAMP=DATE(7:8)//'-'//DATE(5:6)//'-'//DATE(1:4)//' '//
& TIME(1:2)//':'//TIME(3:4)//' '//ZONE(1:3)//':'//ZONE(4:5)
RETURN
END
************************************************************************
SUBROUTINE UNDUMP(LUNRES,DUMPFN)
************************************************************************
IMPLICIT DOUBLE PRECISION (A-H)
IMPLICIT INTEGER*4 (I-N)
IMPLICIT DOUBLE PRECISION (O-Z)
C
CHARACTER DUMPFN*40
PARAMETER(MTARGY=50,MTARGZ=100)
DIMENSION YTARG(-MTARGY:MTARGY),ZTARG(MTARGZ)
COMMON /TARPOS/ DYTARG,DZTARG,YTARG,ZTARG,NTARGY,NTARGZ ! HR@05-Aug-2023
PARAMETER(MAXNU=40)
DIMENSION PNU(-MTARGY:MTARGY,MTARGZ,0:MAXNU), ! P_nu averaged over events
& TM1(-MTARGY:MTARGY,MTARGZ) ! M1 averaged over events
COMMON /TSCORE/ PNU,TM1
COMMON /HISTOR/ MXG4EV,NEVTS,NIONS,NNULL
C -----------------------------------------------------------------
C! Read dump file
OPEN(LUNRES,STATUS='OLD',ERR=10,FILE=DUMPFN)
C! Read metainfo
READ(LUNRES,*) NTY,NTZ,NUMAX,MG4,NEV,NIO,NNU
IF(NTY.EQ.NTARGY.AND.NTZ.EQ.NTARGZ) THEN
MXG4EV=MG4
NEVTS=NEV
NIONS=NIO
NNULL=NNU
C! Read the results
DO ITY=-NTARGY,NTARGY
DO ITZ=1,NTARGZ
READ(LUNRES,*) TM1(ITY,ITZ),(PNU(ITY,ITZ,N),N=0,NUMAX)
END DO
END DO
CLOSE(LUNRES,STATUS='KEEP')
ELSE
CLOSE(LUNRES,STATUS='KEEP')
PRINT*,'ERROR: Dump file is not matching present target numbers'
PRINT*,'You should delete the dump file and modify the first'//
& ' file number in the command file.'
PRINT*,'NTARGY NTARGZ ',NTARGY,NTARGZ
STOP
ENDIF
10 CONTINUE
END ! UNDUMP
************************************************************************
SUBROUTINE WRTRES(IDUMP)
************************************************************************
IMPLICIT DOUBLE PRECISION (A-H)
IMPLICIT INTEGER*4 (I-N)
IMPLICIT DOUBLE PRECISION (O-Z)
PARAMETER(MTARGY=50,MTARGZ=100)
DIMENSION YTARG(-MTARGY:MTARGY),ZTARG(MTARGZ)
COMMON /TARPOS/ DYTARG,DZTARG,YTARG,ZTARG,NTARGY,NTARGZ ! HR@05-Aug-2023
C! /LSCALE/ DCF conversion factor from length in mm to mass per area
C! in µg/nm**2.
COMMON /LSCALE/ DCF
COMMON /INTVOL/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX ! dimensions in mm
PARAMETER(MAXNU=40)
DIMENSION PNU(-MTARGY:MTARGY,MTARGZ,0:MAXNU), ! P_nu averaged over events
& TM1(-MTARGY:MTARGY,MTARGZ) ! M1 averaged over events
COMMON /TSCORE/ PNU,TM1
COMMON /HISTOR/ MXG4EV,NEVTS,NIONS,NNULL
C! /LUNMBS/ Logical unit numbers of input and output files and their names
CHARACTER DUMPFN*40,OUTFNM*40,PREFIX*20
COMMON /LUNMBS/ LUNINP(3),LUNOUT,LUNRES,DUMPFN,OUTFNM,PREFIX
CHARACTER*80 FILENM,HEADER
COMMON /OUTPT0/ FILENM,HEADER
CHARACTER*50 FORMT,FORMT2,WRLINE*24,TXTPNU(0:MAXNU)*12
CHARACTER*24 TSTAMP
DIMENSION PNULO(0:MTARGY,MTARGZ,0:MAXNU), ! P_nu averaged over events
& TM1LOC(0:MTARGY,MTARGZ) ! M1 averaged over events
C-----------------------------------------------------------------------
C! Find maximum NU with non-zero values
NUMAX=MAXNU+1
SUMME=0.0d0
DO WHILE((SUMME.EQ.0.).AND.(NUMAX.GT.0))
NUMAX=NUMAX-1
DO ITY=-NTARGY,NTARGY
DO ITZ=1,NTARGZ
SUMME=SUMME+PNU(ITY,ITZ,NUMAX)
END DO
END DO
END DO
IF(NUMAX.LT.0) NUMAX=0
C! When file was accomplished, save intermediate results to dump file
IF(IDUMP.EQ.1) THEN
OPEN(LUNRES,STATUS='UNKNOWN',FILE=DUMPFN)
C! Write metainfo
WRITE(LUNRES,*) NTARGY,NTARGZ,NUMAX,MXG4EV,NEVTS,NIONS,NNULL
C! Write the results
DO ITY=-NTARGY,NTARGY
DO ITZ=1,NTARGZ
WRITE(LUNRES,*) TM1(ITY,ITZ),(PNU(ITY,ITZ,N),N=0,NUMAX)
END DO
END DO
CLOSE(LUNRES,STATUS='KEEP')
ENDIF
C! Normalize
FNORM=0.5d0/NEVTS
PNULL=FLOAT(NNULL)*FNORM
DO ITY=0,NTARGY
DO ITZ=1,NTARGZ
TM1LOC(ITY,ITZ)=(TM1(ITY,ITZ)+TM1(-ITY,ITZ))*FNORM
DUMMY=0.
DO NU=0,NUMAX
PNULO(ITY,ITZ,NU)=(PNU(ITY,ITZ,NU)+PNU(-ITY,ITZ,NU))*FNORM
DUMMY=DUMMY+PNU(ITY,ITZ,NU)
END DO
C! IF(ITY.EQ.0) PRINT*,'CHECK ME',ITY,ITZ,NEVTS,NINT(DUMMY)
END DO
END DO
C! Column header entries
DO NU=0,NUMAX
WRITE(TXTPNU(NU),'(I3)') NU
TXTPNU(NU)='P'//ADJUSTL(TXTPNU(NU))
TXTPNU(NU)=ADJUSTR(TXTPNU(NU))
END DO
OPEN(LUNRES,STATUS='UNKNOWN',FILE=OUTFNM)
WRITE(LUNRES,'(A)') '# Results from programm '//FILENM
WRITE(LUNRES,'(A)') '# File name: '//OUTFNM//' (produced at '
& //TSTAMP()//')'
WRITE(LUNRES,'(2(a,i8),a,i10,2(a,i8),a)') '# EventID | MXG4EV=',
& MXG4EV,' | NEVTS=',NEVTS,' | NIONS=',NIONS,
& ' | NNULL=',NNULL, ' | '
C! Write column headers
FORMT='(a, a)'
WRITE(FORMT(4:6),'(I3)') NUMAX+1
WRITE(LUNRES,FORMT) '# y/mm z/mm M1',
& (TXTPNU(NU),NU=0,NUMAX)
C! Write the results
FORMT='(2F8.3, E12.5)'
WRITE(FORMT(8:10),'(I3)') NUMAX+2
DO ITY=0,NTARGY
DO ITZ=1,NTARGZ
C*! WRITE(LUNRES,FORMT) YTARG(ITY)/DCF,ZTARG(ITZ)/DCF,
WRITE(LUNRES,FORMT) YTARG(ITY),ZTARG(ITZ),
& TM1LOC(ITY,ITZ),(PNULO(ITY,ITZ,N),N=0,NUMAX)
END DO
END DO
CLOSE(LUNRES,STATUS='KEEP')
END
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment