Proton beam with realistic geometry

2014-08-24T07:45:59Z

78.146.154.109:

== Introduction ==

This example shows the dose distribution in water along the incident proton beam. This example is very similar to the monoenergetic proton pencil beam example. The difference is that the beam is defined with realistic geometry. For the generation of the proton beam instead of particle gun, we use general particle source. For more details about '''G4GeneralParticleSource''' class look [https://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch02s07.html here].

The volume of the water cube is divided into slices perpendicular to the incident beam. The slices are created using class '''G4PVReplica'''. The energy and dose are scored using classes '''G4UserSteppingAction''' and '''G4UserRunAction'''. Alternatively, the energy and dose are scored using class '''G4ScoringManager''' by defining scoring mesh. There is an option to chose among several EM and '''QGSP_BIC_EMY''' physics lists.

== Setting up the environment ==

; Connect to HEP cluster and create folder ProtonGBFolder in your area

<pre style="color: #800000; background-color: #dcdcdc">
ssh -X username@plus1.hep.ucl.ac.uk

password: type your password here

cd /home/username/

mkdir ProtonGBFolder

cd ProtonGBFolder
</pre>

; Setup your environment

<pre style="color: #800000; background-color: #dcdcdc">
source /unix/pbt/software/dev/bin/pbt-dev.sh
</pre>

== How to get the code ==

; Copy the code to your working directory and rename it

<pre style="color: #800000; background-color: #dcdcdc">
cp -r /unix/pbt/tutorials/basic/ProtonGB .

mv ProtonGB ProtonGB_source
</pre>

== How to run the code ==

; Inside /home/username/ProtonGBFolder/ create a directory

<pre style="color: #800000; background-color: #dcdcdc">
mkdir ProtonGB_build
</pre>

; To compile the code enter this directory and run cmake and make

<pre style="color: #800000; background-color: #dcdcdc">
cd ProtonGB_build

cmake -DGeant4_DIR=/unix/pbt/software/dev /home/username/ProtonGBFolder/ProtonGB_source

make
</pre>

; Run macro proton.mac.

<pre style="color: #800000; background-color: #dcdcdc">
./protonGB proton.mac
</pre>

== How to analyze data ==

The macro produces a root file '''Proton.root''' with two histograms. The first histogram shows the energy deposition in water box along the beam line, the second histogram shows zoomed energy deposition around the peak. The macro also produces several text files.

* The data in files '''DoseFile.txt''' and '''PlotDose.txt''' was created using classes '''G4UserSteppingAction''' and '''G4UserRunAction'''. The file '''DoseFile.txt''' contains energy and dose deposition for every layer. The file '''PlotDose.txt''' contains only depth vs dose for each layer. These text files can be analyzed with MATLAB or ROOT.

* The files '''DoseLongitudinalMesh.txt''', '''EnergyLongitudinalMesh.txt''', '''DoseLateralMesh.txt''' and '''EnergyLateralMesh.txt''' contain information about the dose and energy deposition in voxels in longitudinal and lateral direction of the beam. The data was created using class '''G4ScoringManager''' and commands '''/score/''' in proton.mac. These text files can be analyzed with MATLAB or ROOT.

* The two ways to record data should give similar result.

=== Text files ===

This is an example output for '''DoseFile.txt''' with physics process '''QGSP_BIC_EMY''' and incident proton energy of '''62 MeV'''.
Use your favorite editor '''pico''', '''vi''', '''emacs''' etc. to open text files. Now, open the text file '''DoseFile.txt''':

<pre style="color: #800000; background-color: #dcdcdc">
pico DoseFile.txt
</pre>

<pre style="color: #800000; background-color: #dcdcdc">
Layers x[mm] Edep Edep/Ebeam[%] Dose Dose/MaxDose[%]
layer 1: 1.33333 8.80235 GeV 2.366 1.058e-05 Gy 21.8
layer 2: 2.667 8.898 GeV 2.392 1.069e-05 Gy 22.03
layer 3: 4 9.213 GeV 2.476 1.107e-05 Gy 22.81
layer 4: 5.333 9.374 GeV 2.52 1.126e-05 Gy 23.21
layer 5: 6.667 9.594 GeV 2.579 1.153e-05 Gy 23.76
layer 6: 8 9.91 GeV 2.664 1.191e-05 Gy 24.54
layer 7: 9.333 10.07 GeV 2.708 1.211e-05 Gy 24.95
layer 8: 10.67 10.26 GeV 2.758 1.233e-05 Gy 25.41
layer 9: 12 10.55 GeV 2.835 1.267e-05 Gy 26.11
layer 10: 13.33 10.88 GeV 2.925 1.308e-05 Gy 26.95
layer 11: 14.67 11.21 GeV 3.012 1.347e-05 Gy 27.75
layer 12: 16 11.64 GeV 3.13 1.399e-05 Gy 28.83
layer 13: 17.33 11.96 GeV 3.215 1.437e-05 Gy 29.62
layer 14: 18.67 12.5 GeV 3.361 1.502e-05 Gy 30.96
layer 15: 20 12.99 GeV 3.492 1.561e-05 Gy 32.17
layer 16: 21.33 13.71 GeV 3.686 1.648e-05 Gy 33.95
layer 17: 22.67 14.5 GeV 3.899 1.743e-05 Gy 35.91
layer 18: 24 15.37 GeV 4.131 1.847e-05 Gy 38.06
layer 19: 25.33 16.58 GeV 4.457 1.993e-05 Gy 41.06
layer 20: 26.67 18.18 GeV 4.887 2.185e-05 Gy 45.02
layer 21: 28 20.45 GeV 5.496 2.457e-05 Gy 50.63
layer 22: 29.33 24.39 GeV 6.555 2.93e-05 Gy 60.38
layer 23: 30.67 32.75 GeV 8.803 3.935e-05 Gy 81.09
layer 24: 32 40.39 GeV 10.86 4.853e-05 Gy 100
layer 25: 33.33 3.931 GeV 1.057 4.724e-06 Gy 9.734
layer 26: 34.67 1.427 MeV 0.0003836 1.715e-09 Gy 0.003534
layer 27: 36 903.1 keV 0.0002427 1.085e-09 Gy 0.002236
layer 28: 37.33 7.17 MeV 0.001927 8.616e-09 Gy 0.01775
layer 29: 38.67 972.7 keV 0.0002615 1.169e-09 Gy 0.002408
layer 30: 40 0 eV 0 0 Gy 0

The run consists of 6000 protons of 62 MeV through 4 cm of Water (density: 1 g/cm3 ) divided into 30 slices.

Edep is the deposited energy in every slice.
Total incident energy(Ebeam)= 372 GeV
Total energy deposit= 366.7 GeV
Dose is the deposited dose in every slice.
MaxDose is the highest dose value from all slices.
</pre>

The corresponding '''PlotDose.txt''' is:

<pre style="color: #800000; background-color: #dcdcdc">
1.33333 21.7953
2.66667 22.033
4 22.8122
5.33333 23.2101
6.66667 23.7561
8 24.5372
9.33333 24.9451
10.6667 25.4071
12 26.114
13.3333 26.9454
14.6667 27.7468
16 28.8325
17.3333 29.6189
18.6667 30.9594
20 32.1666
21.3333 33.9549
22.6667 35.9146
24 38.0563
25.3333 41.0607
26.6667 45.0231
28 50.6268
29.3333 60.3833
30.6667 81.0937
32 100
33.3333 9.73398
34.6667 0.00353384
36 0.00223614
37.3333 0.0177538
38.6667 0.00240845
40 0
</pre>

[http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonGB/DoseLongitudinalMesh.txt '''DoseLongitudinalMesh.txt'''] and [http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonGB/EnergyLongitudinalMesh.txt '''EnergyLongitudinalMesh.txt'''] contain information about the dose and energy deposition in 30 voxels along the beam.

[http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonGB/DoseLateralMesh.txt '''DoseLateralMesh.txt'''] and [http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonGB/EnergyLateralMesh.txt '''EnergyLateralMesh.txt'''] contain information about the dose and energy deposition in 30 voxels in direction perpendicular to the beam at its peak location along the beam.

=== Root file ===

Open the '''Proton.root''' file in the following way:

<pre style="color: #800000; background-color: #dcdcdc">
root -l Proton.root

new TBrowser

Select ROOT Files and Proton.root
</pre>

 '''This is the energy deposition along the beam in the absorber:''' 

http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonGB/BraggPeak_GB.png

 '''This is the energy deposition along the beam in the absorber, zoomed around the peak:''' 

http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonGB/BraggPeak_GBzoom.png

You can use script '''PlotSimulation.C''' from folder '''RootScripts''' to plot the dose deposition along the absorber. This script uses '''PlotDose.txt'''. Copy the script from '''/ProtonGB_source/RootScripts/''' to your build directory as it was done in the '''Monoenergetic proton pencil beam''' tutorial and run it:

<pre style="color: #800000; background-color: #dcdcdc">
root -l

.x PlotSimulation.C
</pre>

This will create root file with the following plot:

http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonGB/Simulation.png

=== Changes in proton.mac ===

You can change the physics process, incident proton energy and number of slices by
modifying the macro proton.mac. In addition, you can configure the particle source by using '''/gps/''' commands. The beam characteristics used in this macro are similar to the ones used at the Laboratori Nazionali del Sud (INFN) in Catania, Italy and the Clatterbridge Cancer Center(energy distribution is similar). The proton beam has Guassian energy distribution. [https://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch02s07.html Here] you can learn more about the different /gps/ commands.

Open the macro proton.mac:

<pre style="color: #800000; background-color: #dcdcdc">
pico proton.mac
</pre>

You will see:

<pre style="color: #800000; background-color: #dcdcdc">
# proton.mac
#
/control/verbose 2
/run/verbose 2
/tracking/verbose 0
/run/particle/verbose 1
/run/particle/dumpList
#
# set geometry
/protonGB/det/setSizeX 4 cm
/protonGB/det/setSizeYZ 4 cm
/protonGB/det/setSliceSizeYZ 1 cm
/protonGB/det/sliceNumber 30
#
# define longitudinal scoring mesh
# along the beam
/score/create/boxMesh waterMeshlongitudinal
/score/mesh/boxSize 2. 2. 2. cm
/score/mesh/nBin 30 1 1
/score/mesh/translate/xyz 0. 0. 0. cm
/score/quantity/energyDeposit energyDeposit
/score/quantity/doseDeposit doseDeposit
/score/close
#
# define lateral scoring mesh
# centered at the Bragg peak
/score/create/boxMesh waterMeshlateral
/score/mesh/boxSize 1. 2. 2. cm
/score/mesh/nBin 1 30 1
/score/mesh/translate/xyz 1.2 0. 0. cm
/score/quantity/energyDeposit energyDeposit
/score/quantity/doseDeposit doseDeposit
/score/close
#
# set physics process
/protonGB/phys/addPhysics QGSP_BIC_EMY
#/protonGB/phys/addPhysics emlivermore
#/protonGB/phys/addPhysics empenelope
#
# production tresholds (recommended range
#cut off not bigger than 10% of slice thickness)
/protonGB/phys/setCuts 0.2 mm
#/protonGB/phys/setGCut 1 um
#/protonGB/phys/setECut 1 um
#/protonGB/phys/setPCut 1 um
#
# initialize
/run/initialize
# General particle source
# proton circle source
/gps/pos/shape Circle
/gps/pos/centre -4. 0. 0. cm
/gps/pos/radius 0. mm
/gps/pos/sigma_r 2. mm
/gps/particle proton
/gps/pos/type Beam
#
# the incident surface is in the y-z plane
/gps/pos/rot1 0 1 0
/gps/pos/rot2 0 0 1
#
# the beam is travelling along the x-axis without any angular
#dispersion (angular despersion set to 0.0)
/gps/ang/rot1 0 0 1
/gps/ang/rot2 0 1 0
/gps/ang/type beam1d
/gps/ang/sigma_r 0. deg
#
# the beam energy is in gaussian profile
/gps/ene/type Gauss
/gps/ene/mono 62 MeV
/gps/ene/sigma 0.3 MeV
#
# step limit (recommended not bigger than 5% of
# slice thickness)
/protonGB/stepMax 0.1 mm
#
/protonGB/event/printModulo 50
#
# output file
/analysis/setFileName Proton
#
/analysis/h1/set 2 50 25 35 mm
# number of events
/run/beamOn 6000
#
# dump scores to a file
/score/dumpQuantityToFile waterMeshlongitudinal doseDeposit DoseLongitudinalMesh.txt
/score/dumpQuantityToFile waterMeshlongitudinal energyDeposit EnergyLongitudinalMesh.txt
/score/dumpQuantityToFile waterMeshlateral doseDeposit DoseLateralMesh.txt
/score/dumpQuantityToFile waterMeshlateral energyDeposit EnergyLateralMesh.txt
</pre>

You can modify '''proton.mac''' as it is done in the tutorial '''Monoenergetic proton pencil beam'''.

=== Visualisation ===

Under construction (problems with OpenGL installed on plus1)

== Comparison with data from The Clatterbridge Cancer Centre ==

Compare simulation with data using ROOT macros in folder '''RootScripts'''. These scripts are similar to the ones used in the tutorial '''Monoenergetic proton pencil beam'''. For example, by using '''PlotDataAndSim.C''' you can compare proton data from Clatterbridge with simulation (PlotDose.txt).

http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonGB/ClatterbridgeSimulation.png

== Files ==

[[List of proton beam with realistic geometry files with brief description]]

UCL HEP PBT Wiki - User contributions [en]

Proton beam with realistic geometry