Upload New 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