Merge branch 'Hans.Rabus-master-patch-80115' into 'master'

Add new directory

See merge request !1

Simulation_Track_Structure/GDL_Data_Analysis_Codes/Ptra4ICSD.pro

+function get_one_file,file,eventf
+  line=strarr(file_lines(file)) & openr,u,/get,file & readf,u,line & free_lun,u
+  nhline=where(strpos(line,'GEANT4') gt 0,nw)
+  if nw gt 0 then begin
+    eventf=strsplit(line[nhline],'= ',/extract) 
+    eventf=float(eventf[10+3*indgen(7)])/float(eventf[4])
+  endif
+  nhline=max(where(strmid(line,0,1) eq '#')) & line=line[nhline+1:*] 
+  dd=float([[strmid(line,2,7)],[strmid(line,9,7)],[strmid(line,16,12)]])
+  xpos=dd[*,0] & zpos=dd[*,1]
+  qq=where(xpos eq 0 and zpos eq 133.8,nq)
+  line=strsplit(line[qq[0]],/extract)
+  return,float(line[3:*])
+end ;function
+
+function headline,tabstart=tabstart,param,what
+  if keyword_set(tabstart) then txt=string(9b) else txt=''
+  if vtype(param) eq 'FLOAT' then thick=string(param,f='(f11.7)') else thick=param
+  if n_params() eq 1 or what eq '' then begin
+    bra='' & ket=''
+  endif else begin 
+    bra='[' & ket=']'
+  endelse
+  txt += string(what+bra+thick[0]+ket,f='(a11)')
+  for j=1,n_elements(thick)-1 do txt += string(9b)+string(what+bra+thick[j]+ket,f='(a11)')
+  return,txt
+  end ; function headline
+
+pro output_ICSD,subdir,which,thick,events,dat,yplot=yplot,the_case
+  nfiles=n_elements(thick)
+  ; write results for M1(d) -------------------------------------------
+  openw,u,/get,subdir+'PTra4HIT_ICSD'+which+'.dat' 
+  
+  printf,u,headline(thick,'ev',/tab)
+  printf,u,headline(events,'',/tab)
+  printf,u,headline(thick,'P_n',/tab)
+  
+  txt=string(indgen(n_elements(dat[*,0])),f='(i3)')
+  for j=0,nfiles-1 do txt=txt+string(9b)+string(dat[*,j],f='(f11.5)')
+  printf,u,transpose(txt)
+  free_lun,u
+  print,f='(/,a,/)','# Output written to '+subdir+'PTra4HIT_ICSD'+which+'.dat '
+  
+  ; plot ICSD  -------------------------------------------------------
+  if keyword_set(yplot) then begin
+    plot,dat[*,0],/ylog,ps=2
+    for i=0,nfiles-1 do oplot,dat[*,i],ps=2,col=i+1
+  endif
+end ; ------------------------------------------------------------------
+
+pro Ptra4ICSDcase,subdir,date,bor=bor,yplot=yplot
+
+  if not keyword_set(date) then date='*'
+
+; identify which files to read --------------------------------------
+  the_case=keyword_set(bor)+2*keyword_set(out)+4*keyword_set(other)
+  case the_case of
+    0: files=file_search(subdir+'*mm','*_HIT_both_ICSD_*'+date+'*.dat')
+    1: files=file_search(subdir+'*mm','*_HIT_bothB_ICSD_*'+date+'*.dat')
+    2: files=file_search(subdir+'*mm','*_HIT_out_ICSD_*'+date+'*.dat')
+    3: files=file_search(subdir+'*mm','*_HIT_outB_ICSD_*'+date+'*.dat')
+    else: files=file_search(subdir+'*mm','*_HIT_other_ICSD_*'+date+'*.dat')
+  endcase
+
+; find the most recent ones -------------------------------------
+  files=files[uniq(strmid(files,0,max(strpos(files,'/P'))))]
+  nfiles=n_elements(files)
+;
+  
+  if nfiles eq 0 then begin & print,'ERROR: not enough files' & stop & end 
+
+; sort file names ascending in absorber thickness ---------------'
+  thick=files 
+  the_date=''
+  for i=0,n_elements(files)-1 do begin
+    thick[i]=strmid(thick[i],strpos(thick[i],'/')+1) 
+    the_date=max([the_date,strmid(files[i],strlen(files[i])-12,8)])
+    if strpos(files[0],'1/') gt 0 then begin
+      thick[i]=strmid(thick[i],0,strpos(thick[i],'mm')) 
+    endif else begin
+      case strmid(thick[i],0,1) of
+        'b': thick[i]='1'
+        'g': thick[i]='2'
+        'p': thick[i]='3'
+        'r': thick[i]='0'
+      endcase
+    endelse 
+  endfor
+  
+  so=sort(float(thick))
+  files=files[so] & thick=thick[so]
+  
+  if strpos(files[0],'2/') gt 0 then begin
+    for i=0,n_elements(files)-1 do begin
+      case thick[i] of
+        '0': thick[i]='red'
+        '1': thick[i]='blue'
+        '2': thick[i]='green'
+        '3': thick[i]='pink'
+      endcase
+    endfor
+  endif 
+
+; define the suffices of the output files
+  which=strmid(files[0],0,strlen(files[0])-17)
+  which=strmid(which,strpos(which,'HIT')+3)
+  if which eq '_both_' then which='_bothC_'
+  if which eq '_out_' then which='_outC_'
+  which=which+the_date
+ 
+; read the data into array "what" ------------------------------
+  for i=0,n_elements(files)-1 do begin
+    print,f='(20X,a)',files[i]
+    if i eq 0 then begin
+      what=get_one_file(files[i],eventf) 
+    endif else begin
+      was=get_one_file(files[i],events)
+      eventf=[[eventf],[events]]  
+      if i eq 1 then begin
+        n1=n_elements(what)
+        n2=n_elements(was)
+        for j=n1+1,n2 do what=[what,0]
+        for j=n2+1,n1 do was=[was,0]
+        what=[[what],[was]] 
+      endif else begin
+        siz=size(what) & n1=siz[1]
+        n2=n_elements(was)
+        for j=n1+1,n2 do what=transpose([[transpose(what)],[replicate(0,i)]])
+        for j=n2+1,n1 do was=[was,0]
+        what=[[what],[was]] 
+      endelse
+    endelse
+  endfor
+
+; write results for ICSD -------------------------------------------
+  output_ICSD,subdir,which,thick,eventf[the_case,*], $
+                 what,yplot=yplot,the_case
+  
+end ; ------------------------------------------------------------------
+
+pro unita,subdir,date,what
+  files1=file_search(subdir,'P*ICSD_*'+what+'C*'+date+'*')
+  files2=file_search(subdir,'P*ICSD_*'+what+'B*'+date+'*')
+  for i=0,n_elements(files1)-1 do begin
+    pos=strpos(files1[i],what)+strlen(what)
+    file=strmid(files1[i],0,pos)+strmid(files1[i],pos+1)
+    lines1=strarr(file_lines(files1[i])) & lines2=lines1
+    openr,u,/get,files1[i] & readf,u,lines1 & free_lun,u
+    ev1=float(strsplit(lines1[1],string(9b),/extra))
+    openr,u,/get,files2[i] & readf,u,lines2 & free_lun,u
+    ev2=float(strsplit(lines2[1],string(9b),/extra))
+    
+    nfiles=n_elements(ev1)
+      
+    openw,u,/get,file
+    printf,u,lines1[0]
+    txt=string(ev1[0]+ev2[0],f='(f11.7)')
+    if strpos(files1[i],'(0)') eq -1 then txt = string(9b)+txt
+    for j=1,nfiles-1 do txt += string(9b)+string(ev1[j]+ev2[j],f='(f11.7)')
+    printf,u,txt
+    printf,u,lines1[2]
+    
+    lines1=lines1[3:*] & lines2=lines2[3:*]
+    
+    x1=dblarr(nfiles) & x2=x1
+    pos=1.
+    for j=0,n_elements(lines1)-1 do begin
+      reads,lines1[j],pos,x1
+      reads,lines2[j],pos,x2
+      xx=(ev1*x1+ev2*x2)/(ev1+ev2)
+      txt=string(pos,f='(f11.5)')
+      for k=0,nfiles-1 do txt += string(9b)+string(xx[k],f='(f9.5)')
+      printf,u,txt
+    endfor
+    free_lun,u
+    print,f='(20X,a,/)',file
+  endfor
+end
+
+pro Ptra4ICSD,date,yplot=yplot
+
+  if not keyword_set(date) then date='*'
+  
+  subdirs=['largeBeamExp1','largeBeamExp2','smallBeamExp1','smallBeamExp2']+'/'
+  
+  for i=0,0 do begin
+    subdir=subdirs[i]
+    Ptra4ICSDcase,subdir,date,bor=0,yplot=yplot
+    Ptra4ICSDcase,subdir,date,bor=1,yplot=yplot
+    unita,subdir,date,'both'
+  endfor
+
+end
\ No newline at end of file

Simulation_Track_Structure/GDL_Data_Analysis_Codes/README.md

+## How to use the GDL codes
+The data analysis of the output files produced by the PTra simulations was performed using the GNU data language environment which can be generally directly installed in most Linux distributions.
+
+1. Download all *.pro files into a directory in your default GDL path (encoded in the GDL constant !PATH) or into the directory mentioned in the next point. <br />
+2. Open a terminal and navigate to the directory containing as subfolders the directories containing the data produced by the PTra simulations for the different cases (thickness of PMMA absorber and primary carbon ion energy). <br />
+3. Start GDL and in the GDL prompt execute the routines Ptra4MS and Ptra4ICSD. <br />
+   GDL> Ptra4MS  <br />
+   GDL> Ptra4ICSD  <br />
+The results are written into the current directory. (The data produced in the frame of the paper have been uploaded as *results.zip* in this subdirectory of the repository.<br />
+For producing the plots for the publication, the data were copy-pasted from these Ascii-files into **Origin** worksheets.
+
+## Contact
+
+Hans Rabus - hans.rabus@ptb.de <br />
+
+## License
+
+GNU GENERAL PUBLIC LICENSE. See `LICENSE.txt` for more information.
+
+## Citation
+
+If you use the code for your research, please consider citing:
+
+Miriam Schwarze, Gerhard Hilgers, and Hans Rabus. "Nanodosimetric investigation of the track structure of therapeutic carbon ion radiation part 2: Detailed simulation." 10.1088/2057-1976/ad9152.<br />
+
+

Simulation_Track_Structure/PTra_HIT/FilterExtrAp.f

+      PROGRAM FilterExtrAp
+      IMPLICIT REAL (A-H,O-Z)
+      IMPLICIT INTEGER (I-N)
+      INTEGER*8 NEV,IEV
+      CHARACTER*100 FILNAM,FULLNM,BUFFER
+      CHARACTER*20 STRING
+      CHARACTER*20 DATEN(8,1000),PARTIC,LINES(1000)*80
+      CHARACTER*80 FPATH,FPATT,OPATH,PPATH       ! File name pattern
+      CHARACTER*2 PREFIX(2)
+      CHARACTER ADJUSTL,TRIM
+
+      LOGICAL ACCEPT,BORON,CARBON,GOOD,THERE
+
+      DIMENSION LUNOUT(2)
+
+      DATA PREFIX / 'F_','O_' /
+      DATA LUNOUT / 12,13 /
+
+      PPATH='-'
+      
+      PRINT*,'# files per run'
+      READ(*,*) NF
+
+      PRINT*, 'Number of source directories'
+      READ(*,*) NDIRS
+
+      PRINT*, 'Replace existing output files? (1/0)'
+      READ(*,*) KILL
+      IF(KILL.EQ.0) THEN
+        PRINT*,'Existing output files will be extended.'
+      ELSE
+        PRINT*,'Existing output files will be replaced'
+      ENDIF
+
+      FPATT='ExtractionScoringPlane_t#.csv'
+      LOCF=INDEX(FPATT,'#')
+      
+      Loop_Directories: DO IDIR=1,NDIRS
+
+      READ(*,'(a)') FPATH
+      READ(*,'(a)') OPATH
+      PRINT*,'|'//TRIM(FPATH)//'|'//TRIM(OPATH)//'|'
+      
+      IEV=-1
+      EX1=0.
+      EY1=0.
+      EZ1=0.
+      E1=0.
+      N1=0
+      EX2=0.
+      EY2=0.
+      EZ2=0.
+      E2=0.
+      N2=0
+      IDUM=1
+
+      Loop_Files: DO I=1,NF
+        WRITE(STRING,'(I8)') I-1
+        FILNAM=FPATT(1:LOCF-1)//TRIM(ADJUSTL(STRING))
+     &                        //FPATT(LOCF+1:80)
+        PRINT*, TRIM(FILNAM)//'|'
+
+        FULLNM=TRIM(FPATH)//FILNAM
+        INQUIRE(FILE=FULLNM,EXIST=THERE)
+        IF(.NOT.THERE) THEN
+          PRINT*,'File '//TRIM(FULLNM)//' not found'
+          CYCLE Loop_Files
+        END IF
+
+        OPEN(11,FILE=FULLNM,STATUS='OLD')
+        DO J=1,13
+          READ(11,'(a)') LINES(J)
+        END DO
+
+        DO ICASE=1,2
+          FULLNM=TRIM(OPATH)//PREFIX(ICASE)//FILNAM
+          PRINT*,TRIM(FULLNM)//'|'
+
+          INQUIRE(FILE=FULLNM,EXIST=THERE)
+          IF(PPATH.EQ.OPATH.OR.(THERE.AND.KILL.EQ.0)) THEN
+            OPEN(LUNOUT(ICASE),FILE=FULLNM,STATUS='OLD',ACCESS='APPEND')
+          ELSE
+            OPEN(LUNOUT(ICASE),FILE=FULLNM,STATUS='UNKNOWN')
+            DO J=1,13
+              WRITE(LUNOUT(ICASE),'(a)') LINES(J)
+            END DO
+          END IF
+        END DO
+
+        J=1 
+        READ(11,'(A)',END=10) LINES(J)
+        READ(LINES(J),*) (DATEN(K,J),K=1,8),NEV
+        NEVENT=NEV
+
+        Process_this_File: DO 
+          IDUM=IDUM+1
+          IF(MOD(IDUM,100000).EQ.0) WRITE(*,'(a,$)') '.'
+
+          DO WHILE(NEV.EQ.NEVENT)
+            J=J+1
+            READ(11,'(A)',END=10) LINES(J)
+            READ(LINES(J),*) (DATEN(K,J),K=1,8),NEV
+          END DO
+
+          GOOD=CARBON(DATEN(8,1))
+C!        Make sure the first output line is carbon ion (if there)
+          IF(.NOT.GOOD) THEN 
+            DO K=2,J-1
+              GOOD=CARBON(DATEN(8,K))
+              IF(GOOD) THEN 
+                BUFFER=LINES(1)
+                LINES(1)=LINES(K)
+                LINES(K)=BUFFER
+                DO L=1,8
+                  BUFFER=DATEN(L,1)
+                  DATEN(L,1)=DATEN(L,K)
+                  DATEN(L,K)=BUFFER
+                END DO
+                EXIT
+              ENDIF 
+            END DO
+          ENDIF
+
+          IF(GOOD) THEN
+            DO K=1,J-1
+              IF(ACCEPT(DATEN(8,K))) WRITE(LUNOUT(1),'(a)') LINES(K)
+            END DO
+          ELSE
+            DO K=1,J-1
+              IF(ACCEPT(DATEN(8,K))) WRITE(LUNOUT(2),'(a)') LINES(K)
+            END DO
+          END IF
+
+          LINES(1)=LINES(J)
+          J=1
+          READ(LINES(J),*) (DATEN(K,J),K=1,8),NEV
+          NEVENT=NEV
+        END DO Process_this_file
+
+  10    CONTINUE
+        CLOSE(11)
+        CLOSE(LUNOUT(ICASE))
+        PRINT*,'Done!'
+      END DO Loop_Files
+      END DO Loop_Directories
+      END
+
+      LOGICAL FUNCTION ACCEPT(PARTIC)
+      CHARACTER*20 PARTIC
+      LOGICAL CARBON
+      ACCEPT=CARBON(PARTIC).OR.PARTIC(1:2).EQ.'e-'
+     &                     .OR.PARTIC(1:5).EQ.'alpha'
+     &                     .OR.PARTIC(1:5).EQ.'proto'
+     &                     .OR.PARTIC(1:5).EQ.'deute'
+     &                     .OR.PARTIC(1:5).EQ.'trito'
+     &                     .OR.PARTIC(1:2).EQ.'He'
+     &                     .OR.PARTIC(1:2).EQ.'Li'
+     &                     .OR.PARTIC(1:2).EQ.'Be' ! 26-Apr-2024
+     &                     .OR.PARTIC(1:1).EQ.'B'  ! 26-Apr-2024
+     &                     .OR.PARTIC(1:1).EQ.'N'  ! 26-Apr-2024
+      RETURN
+      END
+
+      LOGICAL FUNCTION BORON(PARTIC)
+      CHARACTER*20 PARTIC
+      BORON=PARTIC(1:1).EQ.'B'.AND.PARTIC(2:2).NE.'e'
+      RETURN
+      END
+
+      LOGICAL FUNCTION CARBON(PARTIC)
+      CHARACTER*20 PARTIC
+      CARBON=PARTIC(1:2).EQ.'C1'.OR.PARTIC(1:2).EQ.'C9'
+      RETURN
+      END
\ No newline at end of file

Simulation_Track_Structure/PTra_HIT/HITEIN_240513.txt

+# ### Input command file for PTra_c3h8_20170608HIT.T ##################
+# Blank lines and lines starting with a hashtag are comments (and ignored during reading)
+
+# Program name, input file version (YYMMDD.HH) and verbosity flag
+PTra_C3H8_HIT        240513.07   .F
+# The next line is a header and must not start with a hashtag!!!!
+ 1: IONIZATION CLUSTERS IN PROPANE: ENERGY/SPECTRA, EFF./MAP, OUTPUT/CS
+
+# ### Source definition ###############################################
+#+++ NMAX, NEPRIM, iEnergy +++
+ 1000000    1      0                             
+#++++++++++ EPMIN,EPMAX in MeV; ZATOM,ZPROJ,APROJ,QPROJ +++++++++++++++
+ 0.00000E+00 -1.0000E-03 -1.0000E+00 -1.0000E-00 5.50000E-04 -1.0000E+00
+#++ WEIGHT OF SPECIFIED PARTICLE ENERGY (NEPRIM VALUES: FORMAT 6E12.5) ++
+ 1.00000E+00
+#+++++++++++ BEAM Parameters depending on first input parameter (IBEAM)
+#   IBEAM=0 X0 Y0 Z0 RADAL0 (start position and radius of circular beam)
+#   IBEAM=1 XBEAM, YBEAM (full height and width of rectangular beam)
+#   IBEAM>1 ?????
+#   IBEAM=-1 File pattern and first and last number read from terminal
+      -1  
+# ### End source definition ###########################################
+
+# ### Detector definition #############################################
+# >>> GAS PRESSURE in mbar AND TEMPERATURE in °C: DRUCK,TEMP <<<<<<<<<<
+ 1.20000E+00 2.50000E+01
+# >>> INTERACTION VOLUME: XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX in mm <<<<<<<<<
+-33.75 33.75 -33.75 33.75 13.80 243.80
+# >>> LOCATION of extraction aperture: XEXTAP, YEXTAP, ZEXTAP in mm <<<<
+3.0 -16.25 133.8
+# >>> DETECTION EFFICIENCY <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+#   IEFFIC, EFFIC0, iTS, INTRSQ, EFFICE 
+#   IEFFIC: 0 -> constant efficiency, 1 -> efficiency map 
+#   EFFIC0: Efficiency value (IEFFIC=0), else efficiency scaling factor (for test purposes)
+#   ITS:    Index of the drift time window ! Position of 1st character of efficiency value (in data file with efficiency map)
+#   INTRSQ: 1 -> efficiency map is interpolated in radial direction as
+#                function of the square of the radius 
+#   EFFICE: File name of the efficiency map
+ 1  1.00000E+00  2  1
+../../nanodosimeter3_eff_HIT_1.2mbar-C3H8_drift_V4c_part-1-2_ds2v4508.dat
+# >>> MULTIPLE DETECTORS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+#   MTARGS,XTMAX,DXTARG,ZTMIN,ZTMAX,DZTARG
+#   MTARGS: 0 -> single target, 1 -> multiple targets 
+#   XTMAX:  Maximum impact parameter in mm
+#   DXTARG: Increment of impact parameter in mm     
+#   ZTMIN:  Position of first detector along beam direction in mm
+#   ZTMAX:  Position of last detector along beam direction in mm
+#   DZTARG: Increment of detector position in beam direction in mm
+1 0.70000E+01 2.50000E-01 88.8000E+00 229.800E+00 1.50000E+00 
+# ### End detector definition #########################################
+ 
+# ### SIMULATION OPTIONS ##############################################
+#  ISECEL: electron interactions are included 
+#  IAUGER: Auger relaxation on  |  ICONST: Keep ion energy constant             
+       1       1       0
+#  Output options: IAUS, IOUT -----------------------------------------
+       1       0
+#  Further output options: IONOUT  PREFIX
+#  IONOUT=1 -> ionizations are written to a file PREFIX_'datafilename'
+  0   More_Ionizations
+#
+#+++++++++++ DUMP FILE NAME (written after each processed file)
+DumpfileNew.dat
+#+++++++++++ DEBUG FLAGS for MAIN
+ .F .F .F .F .F .F .F .F
+#+++++++++++ DEBUG FLAGS for READ_PSD
+.F .F .F .F
+# ### End simulation options ##########################################

Simulation_Track_Structure/PTra_HIT/Obsolete/Collimator.f

+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!    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Simulation_Track_Structure/PTra_HIT/Obsolete/ConverSi.f

+      PROGRAM ConvertSi
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      IMPLICIT INTEGER (I-N)
+      INTEGER*8 NEV,IEV
+      CHARACTER*80 FILNAM,LINE
+      CHARACTER*20 STRING
+      CHARACTER*32 FPATT            ! File name pattern
+      CHARACTER ADJUSTL,TRIM
+      
+      FPATT='SiliconDetector_t#.csv'
+      LOC=INDEX(FPATT,'#')
+      
+      PRINT*,'# files'
+      READ(*,*) NF
+      
+      IEV=-1
+      EX1=0.
+      EY1=0.
+      EZ1=0.
+      E1=0.
+      N1=0
+      EX2=0.
+      EY2=0.
+      EZ2=0.
+      E2=0.
+      N2=0
+      
+      DO I=1,NF
+        WRITE(STRING,'(I8)') I-1
+        FILNAM=FPATT(1:LOC-1)//TRIM(ADJUSTL(STRING))
+     &                       //FPATT(LOC+1:32)
+        PRINT*, TRIM(FILNAM)//'#'
+        OPEN(11,FILE=TRIM(FILNAM),STATUS='OLD')
+        DO J=1,10
+          READ(11,*) LINE
+        END DO
+        
+        OPEN(12,FILE='H_'//FILNAM,STATUS='UNKNOWN')
+        WRITE(12,*) 'iev,x/mm,y/mm,z/mm,dX/dZ,dY/dZ'
+        
+        NEVENT=0
+        DO WHILE(1.EQ.1)
+          READ(11,*,END=10) X,Y,Z,EDEP,EKIN,NEV
+          IF(NEV.EQ.IEV) THEN
+            IF(Z.LE.100) THEN
+              N1=N1+1
+              E1=E1+EDEP
+              EX1=EX1+EDEP*X
+              EY1=EY1+EDEP*Y
+              EZ1=EZ1+EDEP*Z
+            ELSE
+              N2=N2+1
+              E2=E2+EDEP
+              EX2=EX2+EDEP*X
+              EY2=EY2+EDEP*Y
+              EZ2=EZ2+EDEP*Z
+            ENDIF
+          ELSE
+            NEVENT=NEVENT+1
+            IF(MOD(NEVENT,1000).EQ.0) WRITE(*,'(1H.,a,$)') ''
+            IF(IEV.GT.-1) THEN
+              IF(N1*N2.GT.0) THEN
+                X1=EX1/E1
+                Y1=EY1/E1
+                Z1=EZ1/E1
+                X2=EX2/E2
+                Y2=EY2/E2
+                Z2=EZ2/E2
+                DIS=(Z2-Z1)
+                W1=(Z2-133.8)/DIS
+                W2=(133.8-Z1)/DIS
+                XC= X1*W1+X2*W2
+                YC= Y1*W1+Y2*W2
+                TX=(X2-X1)/DIS
+                TY=(Y2-Y1)/DIS
+                WRITE(STRING,'(I10)') IEV
+                LINE=ADJUSTL(STRING)
+                LINE=ADJUSTL(LINE)
+                WRITE(STRING,'(F9.3)') XC
+                LINE=TRIM(LINE)//','//ADJUSTL(STRING)
+                WRITE(STRING,'(F9.3)') YC
+                LINE=TRIM(LINE)//','//ADJUSTL(STRING)
+                WRITE(STRING,'(F10.6)') TX
+                LINE=TRIM(LINE)//','//ADJUSTL(STRING)
+                WRITE(STRING,'(F10.6)') TY
+                LINE=TRIM(LINE)//','//ADJUSTL(STRING)
+                WRITE(12,'(a)') TRIM(LINE)
+              ENDIF
+            ENDIF
+            IEV=NEV
+            EX1=0.
+            EY1=0.
+            EZ1=0.
+            E1=0.
+            N1=0
+            EX2=0.
+            EY2=0.
+            EZ2=0.
+            E2=0.
+            N2=0
+          ENDIF  
+        END DO
+  10    CONTINUE
+        CLOSE(11)
+        CLOSE(12)
+        PRINT*,'Done!'
+      END DO
+      END
\ No newline at end of file

Simulation_Track_Structure/PTra_HIT/Obsolete/ConverSiNew.f

+      PROGRAM ConvertSi
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      IMPLICIT INTEGER (I-N)
+      INTEGER*8 NEV,IEV
+      CHARACTER*80 FILNAM,LINE,STRING*20
+      CHARACTER*80 FPATT,OPATT            ! File name pattern
+      
+      NF=28
+      
+      Loop_Directories: DO IDIR=1,4
+      
+        IF(IDIR.EQ.1) THEN
+          FPATT='/home/rabus01/dblue/SiliconDetector_t#.csv'
+          OPATT='/home/rabus01/hpcshare/HIT/blue/E_'//
+     &        'SiliconDetector.csv'
+        ELSE IF(IDIR.EQ.2) THEN
+          FPATT='/home/rabus01/dgreen/SiliconDetector_t#.csv'
+          OPATT='/home/rabus01/hpcshare/HIT/green/E_'//
+     &        'SiliconDetector.csv'
+        ELSE IF(IDIR.EQ.3) THEN
+          FPATT='/home/rabus01/dpink/SiliconDetector_t#.csv'
+          OPATT='/home/rabus01/hpcshare/HIT/pink/E_'//
+     &        'SiliconDetector.csv'
+        ELSE
+          FPATT='/home/rabus01/dred/SiliconDetector_t#.csv'
+          OPATT='/home/rabus01/hpcshare/HIT/red/E_'//
+     &        'SiliconDetector.csv'
+        ENDIF
+      
+        OPEN(12,FILE=OPATT,STATUS='UNKNOWN')
+        WRITE(12,*) 'iev,iHIT,x/mm,y/mm,E'
+          
+        LOCF=INDEX(FPATT,'#')
+        IEV=-1
+
+        Loop_Files: DO I=1,NF
+          WRITE(STRING,'(I8)') I-1
+          FILNAM=FPATT(1:LOCF-1)//TRIM(ADJUSTL(STRING))
+     &                         //FPATT(LOCF+1:80)
+          PRINT*, TRIM(FILNAM)
+          OPEN(11,FILE=TRIM(FILNAM),STATUS='OLD')
+          
+          LINE='#'
+          DO WHILE(LINE(1:1).EQ.'#')
+            READ(11,'(a)') LINE
+          END DO
+          
+          READ(LINE,*) X,Y,Z,EDEP,EKIN,NEV
+          IEV=NEV  
+        
+          NEVENT=1
+          Process_this_File: DO WHILE(NEV.EQ.IEV)
+            IF(MOD(NEVENT,1000).EQ.0) WRITE(*,'(1H.,a,$)') ''
+            
+            READ(11,*,END=10) X,Y,Z,EDEP,EKIN,NEV
+            X1=X
+            Y1=Y
+            E1=EKIN
+            IF(Z.LT.100.) THEN
+              IHIT=1
+              IEV=NEV
+              DO WHILE(NEV.EQ.IEV.AND.Z.LT.100.)  ! Read remaining entries of this event in 1st detector
+                READ(11,*,END=10) X,Y,Z,EDEP,EKIN,NEV
+              END DO
+              IF(NEV.EQ.IEV) IHIT=3
+            ELSE
+              IHIT=2
+            ENDIF
+            
+            DO WHILE(NEV.EQ.IEV) ! Read remaining entries of this event
+              READ(11,*,END=10) X,Y,Z,EDEP,EKIN,NEV
+            END DO
+            CALL SCHREIB(12,IEV,IHIT,X1,Y1,E1)
+
+            IEV=NEV
+            NEVENT=NEVENT+1
+          END DO Process_this_File
+  10      CONTINUE
+  
+          CALL SCHREIB(12,IEV,IHIT,X1,Y1,E1)
+          CLOSE(11)
+          PRINT*,'Done!'
+        END DO Loop_Files
+        CLOSE(12)
+      END DO Loop_Directories
+      END
+      
+      SUBROUTINE SCHREIB(LUN,IEV,IHIT,XC,YC,EC)
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      IMPLICIT INTEGER (I-N)
+      CHARACTER STRING*20,LINE*80
+      CHARACTER ADJUSTL,TRIM
+      INTEGER*8 IEV
+                
+      WRITE(STRING,'(I10)') IEV
+      LINE=ADJUSTL(STRING)
+      WRITE(STRING,'(I1)') IHIT
+      LINE=TRIM(LINE)//' '//STRING(1:1)
+      WRITE(STRING,'(F9.3)') XC
+      LINE=TRIM(LINE)//' '//ADJUSTL(STRING)
+      WRITE(STRING,'(F9.3)') YC
+      LINE=TRIM(LINE)//' '//ADJUSTL(STRING)
+      WRITE(STRING,'(F10.3)') EC
+      LINE=TRIM(LINE)//' '//ADJUSTL(STRING)
+      WRITE(LUN,'(a)') TRIM(LINE)
+      
+      END SUBROUTINE SCHREIB 
\ No newline at end of file

Simulation_Track_Structure/PTra_HIT/Obsolete/Electrode_hits_subs_PSF.f

+      
+************************************************************************
+      SUBROUTINE READPSX(IRFILE,APROJ,ZATOM)
+************************************************************************
+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=84)
+      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
+      COMMON /LUNMBS/ LUNINP(3),LUNBAD(2),LUNOUT,LUNPLT(2),LUNRES
+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/ XCION,COSTHC,SINTHC y-coordinate and direction 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.
+C!19-MAY-2024      COMMON /EVINFO/ YCION,COSTHC,SINTHC,X,Y,Z,E,U,V,W,NRCASE !13-MAY-2024
+      COMMON /EVINFO/ XCION,COSTHC,SINTHC,X,Y,Z,E,U,V,W,NRCASE !19-MAY-2024
+      COMMON /CAVEAT/ ELOWP !13-MAY-2024
+C!    /EVENTS/ IEVC,IEVI, IEVEL are the event IDs associated with the 
+C!             three groups of data.
+      COMMON /EVENTS/ IEVC,IEVI,IEVEL,IEVBAD,MXG4EV
+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,YEXTAP,ZEXTAP
+      
+      PARAMETER(NCASES=7)     ! 13-MAY-2024
+      DIMENSION NEVTS(0:NCASES) ! 13-MAY-2024
+      CHARACTER*9 CASES(0:NCASES)  ! 19-MAY-2024 
+      COMMON /HISTOR/ NG4EVT,NEVTS,NIONS,NNULL,CASES
+
+      LOGICAL NEWEV1, NEWEV2, NODATA ! Flags for new events and no data read (yet)
+      COMMON /NEWST/ NEWEV1, NEWEV2, NODATA
+      
+      LOGICAL LDEBUG(4)
+      COMMON /PSDEBG/ LDEBUG  
+
+      LOGICAL ROTATE ! 13-MAY-2024
+      COMMON /MYSCRT/ ROTATE  ! 13-MAY-2024
+
+      DATA IEVI,IEVEL / -1, -1 / 
+
+      PARAMETER(HWDETX=5.0,HWDETY=1.0) ! 13-MAY-2024 
+C!     Note: Miriam's coordinate system is such that impact parameter is along y
+
+      CHARACTER LINE*100
+C     ------------------------------------------------------------------
+C     Code starts here 
+C     ------------------------------------------------------------------
+      NODATA=.TRUE. 
+     
+C!    Read record from input file
+      READ(LINE,*,END=10,ERR=40) Y2,X2,Z2,Y1,X1,Z1,E,PARTIC,NEV ! 19-MAY-2024
+C
+C!    Identify particle
+      CALL PARTYP(PARTIC,APROJ,ZATOM)
+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
+      
+
+        IZATOM=NINT(ZATOM)
+        IF((IZATOM.EQ.6).OR.(IZATOM.EQ.5)) THEN ! Event with carbon or boron ion
+          IF(E.LE.(ELOWP*APROJ)) THEN ! 13-NOV-2023 Energy is below 2nd ionization threshold.
+            NNULL=NNULL+1
+            NRCASE=12-IZATOM
+            IEVBAD=IEVC
+            IF(LUNBAD(2).GT.0) WRITE(LUNBAD(2),'(i2,i8)') IFILE-1,IEVC !19-MAY-2024
+          ELSE
+C!          Find intersection points with second detector plane 
+            XATSD2=X1+UP/WP*(ZMAX-Z1)
+            YATSD2=Y1+VP/WP*(ZMAX-Z1)
+            IF(ABS(XATSD2).LE.HWDETX.AND.ABS(YATSD2).LE.HWDETY) THEN 
+              NRCASE=7-IZATOM ! C ion trajectory intersects trigger detector
+            ELSE
+              NRCASE=9-IZATOM ! C ion trajectory misses trigger detector
+            END IF            
+          END IF
+        ELSE
+          NRCASE=5 ! Events without C and B ion in nanodosimeters
+        END IF
+      
+        U=UP
+        V=VP
+        W=WP
+        X=X1 
+        Y=Y1
+        Z=Z1
+      
+      IF(APROJ.GT.0.5) THEN ! Move ions back to stART OF INTERACTION VOLUME
+        CALL RELECT(X,Y,Z,U,V,W,E,ZATOM,APROJ,NEV,NRCASE,IFILE) ! 13-MAY-2024 - check whether capacitor plate is hit
+      ENDIF
+      
+      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 ! SUBROUTINE READPSX
+
+************************************************************************
+      SUBROUTINE RELECT(X,Y,Z,U,V,W,E,ZATOM,APROJ,NEV,NRCASE,IFILE)
+************************************************************************
+      IMPLICIT DOUBLE PRECISION (A-H)
+      IMPLICIT INTEGER*4 (I-N)
+      IMPLICIT DOUBLE PRECISION (O-Z)
+      
+      COMMON /EXTRAP/ XEXTAP,YEXTAP,ZEXTAP
+      COMMON /LUNMBS/ LUNINP(3),LUNBAD(2),LUNOUT,LUNPLT(2),LUNRES
+      
+      IZATOM=NINT(ZATOM)
+      
+      YMIN=YEXTAP
+      YMAX=YEXTAP+50.
+      
+      IF(Y.GT.YMIN.AND.Y.LT.YMAX) THEN
+        IF(V.LT.0) THEN
+          FMIN=(YMIN-Y)/V
+          XPLT=ABS(X+FMIN*U-XEXTAP)
+          ZPLT=ABS(Z+FMIN*W-ZEXTAP)
+          IF(DMAX1(XPLT,ZPLT).LT.125.) THEN
+            RADIUS=SQRT(XPLT*XPLT+ZPLT*ZPLT)
+            IF(RADIUS.LE.125.) THEN
+              WRITE(LUNPLT(1),'(i1,2i3,f9.3,f7.2,i8,i3)') NRCASE,
+     &              NINT(ZATOM),NINT(APROJ),E/1.0e3,RADIUS,NEV,IFILE-1
+            ENDIF
+          ENDIF
+        ELSE IF(V.GT.0) THEN
+          FMAX=(YMAX-Y)/V
+          XPLT=ABS(X+FMAX*U-XEXTAP)
+          ZPLT=ABS(Z+FMAX*W-ZEXTAP)
+          IF(DMAX1(XPLT,ZPLT).LT.125.) THEN
+            RADIUS=SQRT(XPLT*XPLT+ZPLT*ZPLT)
+            IF(RADIUS.LE.125.) THEN
+              WRITE(LUNPLT(2),'(i1,2i3,f9.3,f7.2,i8,i3)') NRCASE,
+     &              NINT(ZATOM),NINT(APROJ),E,RADIUS,NEV,IFILE-1
+            ENDIF
+          ENDIF
+        ENDIF
+      ENDIF
+      
+      RETURN
+      END