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

Upload New File

parent ca961ad4
No related branches found
No related tags found
1 merge request!1Add new directory
C! *****************************************************************
PROGRAM Collimator
C! -----------------------------------------------------------------
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT INTEGER (I-N)
PARAMETER(ZCOLL=-112.5,ZEXTAP=133.8)
INTEGER*8 NEV,IEVS,IEVX
LOGICAL CARBON,THERE
CHARACTER*80 FILNAM,LINE
CHARACTER*20 PARTIC,STRING,PATH*80
CHARACTER*80 FPATT,OPATT,SPATT ! File name pattern
C! -----------------------------------------------------------------
PRINT*, 'Number of directories'
READ(*,*) NDIRS
PRINT*, NDIRS
PRINT*,'# files per run'
READ(*,*) NF
PRINT*, NF
PRINT*, 'Replace existing output files? (1/0)' ! 25-APR-2024
READ(*,*) KILL ! 25-APR-2024
IF(KILL.EQ.0) THEN ! 25-APR-2024
PRINT*,'Existing output files will be extended.' ! 25-APR-2024
ELSE ! 25-APR-2024
PRINT*,'Existing output files will be replaced' ! 25-APR-2024
ENDIF ! 25-APR-2024
Loop_Directories: DO IDIR=1,NDIRS
READ(*,'(a)') PATH ! Path to data files
PRINT*,PATH
FPATT=TRIM(PATH)//'ExtractionScoringPlane_t#.csv'
SPATT=TRIM(PATH)//'SiliconDetector_t#.csv'
READ(*,'(a)') PATH ! Path for output
PRINT*,PATH
OPATT=TRIM(PATH)//'Collimator.dat'
PRINT*,OPATT
INQUIRE(FILE=TRIM(OPATT),EXIST=THERE) ! 25-APR-2024
IF(THERE.AND.KILL.EQ.0) THEN ! 25-APR-2024
OPEN(16,FILE=OPATT,STATUS='OLD',ACCESS='APPEND') ! 25-APR-2024
ELSE ! 25-APR-2024
OPEN(16,FILE=OPATT,STATUS='UNKNOWN')
END IF ! 25-APR-2024
WRITE(16,*) 'iev mC xc/mm yc/mm Ec/MeV Edep/keV '//
& 'xt/mm yt/mm Et/MeV'
IEVS=-1
Loop_Files: DO I=1,NF
C! Open data file with interaction of C ion in silicon detectors
CALL OPENPY(11,I-1,SPATT)
C! ... and read first LINE
READ(11,*) X1,Y1,Z1,EDEP,EKIN,NEV
C! Open data file with info on particles passing the target plane
CALL OPENPY(12,I-1,FPATT)
NEVENT=0
Process_this_File: DO
NEVENT=NEVENT+1
IF(MOD(NEVENT,10000).EQ.0) WRITE(*,'(a,$)') '.'
C! Get the next data point in the silicon detectors
DO WHILE(Z1.GT.ZEXTAP) ! Discard events where the C ion missed the first detector
READ(11,*,END=10) X1,Y1,Z1,EDEP,EKIN,NEV
END DO
IEVS=NEV
READ(11,*,END=10) X2,Y2,Z2,EDEP2,EKIN2,NEV
DO WHILE(Z2.EQ.Z1)
READ(11,*,END=10) X2,Y2,Z2,EDEP2,EKIN2,NEV
END DO
IF(NEV.EQ.IEVS) THEN
FACTOR=(ZCOLL-Z1)/(Z2-Z1)
XCOLL=X1+(X2-X1)*FACTOR
YCOLL=Y1+(Y2-Y1)*FACTOR
ELSE ! No second interaction in 1st silicon detector -> skip this event
X1=X2
Y1=Y2
Z1=Z2
EDEP=EDEP2
EKIN=EKIN2
CYCLE Process_this_File
ENDIF
IEVX=IEVS+1 ! The event ID in the extraction plane was counted starting with 1 instead of 0
C! Find the track containing the current event
READ(12,*) X2,Y2,Z2,X1,Y1,Z1,E1,PARTIC,NEV
DO WHILE(NEV.LT.IEVX)
READ(12,*,END=10) X2,Y2,Z2,X1,Y1,Z1,E1,PARTIC,NEV
END DO
C! Find the entry for the carbon ion
IF(NEV.EQ.IEVX) THEN
DO WHILE(.NOT.CARBON(PARTIC).AND.NEV.EQ.IEVX)
READ(12,*,END=10) X2,Y2,Z2,X1,Y1,Z1,E1,PARTIC,NEV
END DO
ENDIF
IF(NEV.EQ.IEVX) THEN ! Carbon ion was found
C! calculate coordinates in target plane
FACTOR=(ZEXTAP-Z1)/(Z2-Z1)
XT=X1+(X2-X1)*FACTOR
YT=Y1+(Y2-Y1)*FACTOR
C! Get the mass number of the C ion
IF(PARTIC(2:2).EQ.'9') THEN
M=9
ELSE
READ(PARTIC(2:3),'(i2)') M
ENDIF
C! write data to output file
CALL SCHREIB(16,IEVS,M,XCOLL,YCOLL,EKIN,EDEP,XT,YT,E1)
C! Advance silicon file to next event
NEV=IEVS
DO WHILE(NEV.EQ.IEVS)
READ(11,*,END=10) X1,Y1,Z1,EDEP,EKIN,NEV
END DO
ENDIF
END DO Process_this_File
10 CONTINUE
CLOSE(11)
CLOSE(12)
PRINT*,'Done!'
END DO Loop_Files
CLOSE(16)
END DO Loop_Directories
END
C! *****************************************************************
LOGICAL FUNCTION CARBON(PARTIC)
C! -----------------------------------------------------------------
CHARACTER*20 PARTIC
CARBON=PARTIC(1:2).EQ.'C1'.OR.PARTIC(1:2).EQ.'C9'
RETURN
END
C! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C! *****************************************************************
SUBROUTINE OPENPY(LUN,NFILE,FPATT)
C! -----------------------------------------------------------------
C! Input variables
CHARACTER*(*) FPATT
INTEGER LUN,NFILE
C! Local variables
CHARACTER*80 FILNAM,LINE,SNFILE*2
INTEGER LINES,LOCF
C! -----------------------------------------------------------------
C! Get the filename
LOCF=INDEX(FPATT,'#')
WRITE(SNFILE,'(I2)') NFILE
FILNAM=FPATT(1:LOCF-1)//TRIM(ADJUSTL(SNFILE))
& //FPATT(LOCF+1:LEN(FPATT))
PRINT*, 'OPENPY '//TRIM(FILNAM)
C! Open file and read the header lines
LINE='#'
LINES=0
OPEN(LUN,FILE=FILNAM,STATUS='OLD')
DO WHILE(LINE(1:1).EQ.'#')
READ(LUN,'(a)') LINE
LINES=LINES+1
END DO
REWIND(LUN)
DO I=1,LINES
READ(LUN,*)
END DO
END SUBROUTINE OPENPY
C! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C! *****************************************************************
SUBROUTINE SCHREIB(LUN,IEV,MPART,XC,YC,EC,ED,XT,YT,ET)
C! -----------------------------------------------------------------
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT INTEGER (I-N)
CHARACTER STRING*20,LINE*80
CHARACTER ADJUSTL,TRIM
CHARACTER SFLOAT*16,SINTEG*8
INTEGER*8 IEV
WRITE(SINTEG,'(I8)') IEV
LINE=ADJUSTL(SINTEG)
WRITE(SINTEG,'(I2)') MPART
LINE=TRIM(LINE)//' '//ADJUSTL(SINTEG)
LINE=TRIM(LINE)//' '//SFLOAT(XC,3)
LINE=TRIM(LINE)//' '//SFLOAT(YC,3)
LINE=TRIM(LINE)//' '//SFLOAT(EC,3)
LINE=TRIM(LINE)//' '//SFLOAT(ED*1000.,3)
LINE=TRIM(LINE)//' '//SFLOAT(XT,3)
LINE=TRIM(LINE)//' '//SFLOAT(YT,3)
LINE=TRIM(LINE)//' '//SFLOAT(ET,3)
WRITE(LUN,'(a)') TRIM(LINE)
END SUBROUTINE SCHREIB
C! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C! ******************************************************************
CHARACTER*16 FUNCTION SFLOAT(X,NDIG)
C! ------------------------------------------------------------------
INTEGER NDIG
REAL*8 X
CHARACTER FMZ*7
FMZ='(F16.0)'
IF(NDIG.GE.10) NDIG=9
WRITE(FMZ(6:6),'(i1)') NDIG
SFLOAT=''
WRITE(SFLOAT,FMZ) X
N=16
DO WHILE(SFLOAT(N:N).EQ.'0')
SFLOAT(N:N)=' '
N=N-1
END DO
SFLOAT=ADJUSTL(SFLOAT)
RETURN
END FUNCTION SFLOAT
C! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment