Compare revisions

## About The Project

This repository contains the code on which the papers 

PAPER 1<br />

PAPER 2

are based on. A two-step simulation of nanodosimetric measurements with the PTB Ion Counting Nanodosimeter of therapeutic-energy carbon ions penetrating simulated tissue is presented and discussed.

![image](setup.png)

Detailed simulations using the Geant4 toolkit were performed to investigate the radiation field in the nanodosimeter and to provide input data for track structure simulations, which were performed with a developed version of the PTra code.

If you are interested in the simulation data, please contact us.

### Structure

1. Simulation_RadiationField <br />

Condensed-history simulation using Geant4 of the experiment, obtaining information about the particle field in TVP and hits in the PSDs

2. Simulation_TrackStructure_Preparation <br />

Condensed history simulation using Geant4 of the experiment, obtaining input data for the track structure simulations

3. Track_Structure_Simulations <br />

Track structure simulation using PTra of the particle field in the target volume of the nanodosimeter

## Requirements

- Geant4 version 11.00.02 

- python version 3.9.6

- numpy version 1.21.2

- matplotlib version 3.4.3 

- pandas version  1.3.1 

- scipy version  1.7.1 

## Contact

Miriam Schwarze - miriam.schwarze@ptb.de <br />

Hans Rabus - hans.rabus@ptb.de <br />

Gerhard Hilgers - gerhard.hilgers@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:

PAPER

This repository contains the code on which the papers 

Hilgers, Gerhard, Miriam Schwarze, and Hans Rabus. "Nanodosimetric investigation of the track structure of therapeutic carbon ion radiation part 1: measurement of ionization cluster size distributions." Biomedical Physics & Engineering Express 10.6 (2024): 065030.<br />

PAPER 2

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 />

are based on. A two-step simulation of nanodosimetric measurements with the PTB Ion Counting Nanodosimeter of therapeutic-energy carbon ions penetrating simulated tissue is presented and discussed.

If you use the code for your research, please consider citing:

Hilgers, Gerhard, Miriam Schwarze, and Hans Rabus. "Nanodosimetric investigation of the track structure of therapeutic carbon ion radiation part 1: measurement of ionization cluster size distributions." Biomedical Physics & Engineering Express 10.6 (2024): 065030. <br />

PAPER2

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 />

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

## 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 />

      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

# ### 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