Monoenergetic proton pencil beam

== <span style="color:#000080"> Introduction </span> ==

This example shows the dose distribution in water along the incident proton beam. The beam hits the water cube surface and deposits a dose under the surface of the water. The volume of the water cube is divided into slices perpendicular to the incident beam. At each slice the deposited dose and energy is computed. 

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 two scoring meshes in logitudinal and lateral direction of the beam. More information about scoring meshes can be found [https://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch04s08.html here]. Protons are generated using '''G4ParticleGun''' class. There is an option to chose among several '''EM''' and '''QGSP_BIC_EMY''' physics lists.

== <span style="color:#000080"> Setting up the environment </span> == 

; Connect to HEP cluster and create folder ProtonPBFolder 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 ProtonPBFolder 

cd ProtonPBFolder  
</pre>

; Setup your environment

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

== <span style="color:#000080"> How to get the code </span> ==

; Copy the code to your working directory and rename it

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

== <span style="color:#000080"> How to run the code </span> == 

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

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

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

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

cmake -DGeant4_DIR=/unix/pbt/software/dev /home/username/ProtonPBFolder/ProtonPB_source 

make  
</pre>

; Run macro proton.mac. The macro generates 6000 events.

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

== <span style="color:#000080"> How to analyze data </span> ==

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 values for each layer. These values can be imported in MATLAB or analyzed with ROOT macro.

* 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. 

* The two ways should give similar result.   

=== <span style="color:#000080"> Text files </span> ===

This is output from '''DoseFile.txt''' with physics process '''QGSP_BIC_EMY''' and incident proton energy of '''62 MeV'''.

<pre style="color: #800000; background-color: #dcdcdc">
 Layers :	 x[mm]    Edep 	    Edep/Ebeam[%]  Dose 	Dose/MaxDose[%]
 layer 1: 	1.33333	8.84429 GeV	2.3775	1.06276e-05 Gy	20.6613	
 layer 2: 	2.66667	8.91694 GeV	2.39703	1.07149e-05 Gy	20.831	
 layer 3: 	4	9.17299 GeV	2.46586	1.10226e-05 Gy	21.4291	
 layer 4: 	5.33333	9.38711 GeV	2.52342	1.12799e-05 Gy	21.9294	
 layer 5: 	6.66667	9.62373 GeV	2.58702	1.15642e-05 Gy	22.4821	
 layer 6: 	8	9.75588 GeV	2.62255	1.1723e-05 Gy	22.7908	
 layer 7: 	9.33333	9.98464 GeV	2.68404	1.19979e-05 Gy	23.3252	
 layer 8: 	10.6667	10.3128 GeV	2.77225	1.23922e-05 Gy	24.0918	
 layer 9: 	12	10.5176 GeV	2.82732	1.26383e-05 Gy	24.5704	
 layer 10: 	13.3333	10.8275 GeV	2.91061	1.30107e-05 Gy	25.2942	
 layer 11: 	14.6667	11.164 GeV	3.00107	1.3415e-05 Gy	26.0803	
 layer 12: 	16	11.6363 GeV	3.12804	1.39825e-05 Gy	27.1837	
 layer 13: 	17.3333	12.0232 GeV	3.23203	1.44474e-05 Gy	28.0874	
 layer 14: 	18.6667	12.6138 GeV	3.3908	1.51571e-05 Gy	29.4672	
 layer 15: 	20	12.9498 GeV	3.48112	1.55609e-05 Gy	30.2521	
 layer 16: 	21.3333	13.6288 GeV	3.66365	1.63768e-05 Gy	31.8383	
 layer 17: 	22.6667	14.5047 GeV	3.89912	1.74294e-05 Gy	33.8847	
 layer 18: 	24	15.4695 GeV	4.15847	1.85887e-05 Gy	36.1385	
 layer 19: 	25.3333	16.5838 GeV	4.45801	1.99276e-05 Gy	38.7416	
 layer 20: 	26.6667	18.3271 GeV	4.92663	2.20224e-05 Gy	42.8141	
 layer 21: 	28	20.6828 GeV	5.5599	2.48532e-05 Gy	48.3174	
 layer 22: 	29.3333	24.4737 GeV	6.57896	2.94084e-05 Gy	57.1734	
 layer 23: 	30.6667	32.4582 GeV	8.72532	3.90028e-05 Gy	75.826	
 layer 24: 	32	42.8062 GeV	11.507	5.14373e-05 Gy	100	
 layer 25: 	33.3333	1.99409 GeV	0.536046     2.39617e-06 Gy	4.65842	
 layer 26: 	34.6667	1.01076 MeV	0.000271708  1.21456e-09 Gy	0.00236124	
 layer 27: 	36	503.954 keV	0.000135472  6.05568e-10 Gy	0.00117729	
 layer 28: 	37.3333	33.8047 keV	9.08729e-06  4.06209e-11 Gy	7.89716e-05	
 layer 29: 	38.6667	2.03563 MeV	0.000547213  2.44608e-09 Gy	0.00475547	
 layer 30: 	40	0 eV 	        0	     0 Gy	        0	


 The run consists of 6000 proton 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= 367.244 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	20.6613	
2.66667	20.831	
4	21.4291	
5.33333	21.9294	
6.66667	22.4821	
8	22.7908	
9.33333	23.3252	
10.6667	24.0918	
12	24.5704	
13.3333	25.2942	
14.6667	26.0803	
16	27.1837	
17.3333	28.0874	
18.6667	29.4672	
20	30.2521	
21.3333	31.8383	
22.6667	33.8847	
24	36.1385	
25.3333	38.7416	
26.6667	42.8141	
28	48.3174	
29.3333	57.1734	
30.6667	75.826	
32	100	
33.3333	4.65842	
34.6667	0.00236124	
36	0.00117729	
37.3333	7.89716e-05	
38.6667	0.00475547	
40	0	
</pre> 

[http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/DoseLongitudinalMesh.txt '''DoseLongitudinalMesh.txt'''] and [http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/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/ProtonPB/DoseLateralMesh.txt '''DoseLateralMesh.txt'''] and [http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/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. 

=== <span style="color:#000080"> Root file </span> ===

Open '''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>

<span style="color:#000080"> '''This is the energy deposition along the beam in the absorber:''' </span>
 
http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/BraggPeak_PB.png 
 
<span style="color:#000080"> '''This is the energy deposition along the beam in the absorber, zoomed around the peak:''' </span>

http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/BraggPeak_PBzoom.png  

You can use macro '''PlotSimulation.C''' to plot the dose deposition along the absorber. This macro uses '''PlotDose.txt'''.

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

=== <span style="color:#000080"> Changes in proton.mac </span> ===

You can change the physics process, incident proton energy and number of slices etc. by
modifying the macro proton.mac. Use your favorite editor '''pico''', '''vi''', '''emacs''' etc. For example open the macro with editor '''pico''':

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

This is what 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 
/protonPB/det/setSizeX  4 cm
/protonPB/det/setSizeYZ 4 cm
/protonPB/det/setSliceSizeYZ 1 cm
/protonPB/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
/protonPB/phys/addPhysics QGSP_BIC_EMY
#/protonPB/phys/addPhysics emlivermore
#/protonPB/phys/addPhysics empenelope
#
# production tresholds (recommended range 
#cut off not bigger than 10% of slice thickness)
/protonPB/phys/setCuts 0.2 mm
#/protonPB/phys/setGCut 1 um
#/protonPB/phys/setECut 1 um
#/protonPB/phys/setPCut 1 um
#
# initialize
/run/initialize
#
/gun/particle proton
# particle energy used in Clatterbridge Centre  
/gun/energy 62 MeV
#
# beam size
#/photonPB/gun/rndm 3 mm
#
# step limit (recommended not bigger than 5% of 
# slice thickness)
/protonPB/stepMax 0.1 mm
#
/protonPB/event/printModulo 50
#
# output file
/analysis/setFileName Proton
#
# histogram 
/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>

'''Change dimensions of the water box'''

The default size is 4x4x4 cm. You can change the dimensions by modifying the lines

<pre style="color: #800000; background-color: #dcdcdc">
/protonPB/det/setSizeX  4 cm
/protonPB/det/setSizeYZ 4 cm 
</pre>

'''Change the physics process'''

The default physics process is '''QGSP_BIC_EMY'''. This is a physics list recommended for proton therapy. You can check what will be the dose deposition if you change the physics list. In proton.mac change

<pre style="color: #800000; background-color: #dcdcdc">
/photonPB/phys/addPhysics QGSP_BIC_EMY
</pre>

to 

<pre style="color: #800000; background-color: #dcdcdc">
/photonPB/phys/addPhysics emlivermore
</pre>

'''Change the incident proton energy'''

The default energy is 62 MeV. This is one of the energies used in radiotherapy. In proton.mac you can change the value of 62 MeV 

<pre style="color: #800000; background-color: #dcdcdc">
/gun/energy 62 MeV
</pre>

to, for example, 50 MeV

<pre style="color: #800000; background-color: #dcdcdc">
/gun/energy 50 MeV
</pre>

'''Change the diameter of the beam'''

You can set the diameter of the beam with the command:

<pre style="color: #800000; background-color: #dcdcdc">
/protonPB/gun/rndm 3 mm
</pre>

'''Change the number of slices'''

You can change the number of slices. The default number is 30. Keep in mind that 
if you want to increase the number of slices you need to modify the file DetectorConstruction.hh in
/ProtonPB_source/include/.

In DetectorConstruction.hh set MaxLayer to value bigger than the number of your slices. The default value is MaxLayer=40. For example, if you want to have 35 slices you do not need to modify MaxLayer. Then, only in proton.mac change the number of slices 

<pre style="color: #800000; background-color: #dcdcdc">
/protonPB/det/sliceNumber 35 
</pre>

'''Modify the mesh'''

You can change the size of the mesh and the number of voxels by modifying 

<pre style="color: #800000; background-color: #dcdcdc">
/score/mesh/boxSize 2. 2. 2. cm
/score/mesh/nBin 30 1 1
</pre>

'''After modifications in proton.mac'''

After modifying the macro proton.mac you can run it

<pre style="color: #800000; background-color: #dcdcdc">
./protonPB proton.mac
</pre>
 
If you modify files in directory ProtonPB_source (for example you change the value of MaxLayer) you need to comppile the code. In directory PhotonPB_build do

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

then run the macro proton.mac.

=== <span style="color:#000080"> Visualisation </span> ===

Under construction

== <span style="color:#000080"> Data from The Clatterbridge Cancer Centre </span> ==

This is [http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/Clatterbridge/ClatterbridgeBraggPeak.txt data] from the Clatterbridge Cancer Center. The first column represents the values in mm in depth, the second column represents the normalized values of deposited dose and normalized at the peak. You can use macro '''PlotData.C''' to plot the data. 

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

== <span style="color:#000080"> Comparison with data from The Clatterbridge Cancer Centre </span> ==

You can use macro '''PlotDataAndSim.C''' to compare data (from the Clatterbridge Cancer Center) and simulation ('''PlotDose.txt'''). 
 
http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/ClatterbridgeGeant4.png

You can also compare data (from the Clatterbridge Cancer Center) with simulation done with scoring mesh ('''DoseLongitudinalMesh.txt''').
This can be done with macro '''PlotDataAndSimMesh.C'''. At the moment this macro works only if before running it you substitude the commas 
in '''DoseLongitudinalMesh.txt''' with spaces. This is the result:

http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/ClatterbridgeGeant4Mesh.png

There is discrepancy between the data in '''PlotDose.txt''' and '''DoseLongitudinalMesh.txt''' (two different ways of retrieving the dose deposition). Probably wrong position of the scoring mesh. This is will be fixed after visualisation.  

== <span style="color:#000080"> Files </span> ==

[[List of monoenergetic proton pencil beam files with brief description]]