Monoenergetic proton pencil beam
From UCL HEP PBT Wiki
for
Monoenergetic proton pencil beam
Jump to:
navigation
,
search
== <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 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. === <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'''. Use your favorite editor '''pico''', '''vi''', '''emacs''' etc to open text files. For example, open this text file with editor '''pico''': <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: 0.8 5.17382 GeV 1.39081 6.47608e-07 Gy 17.3315 layer 2: 1.6 5.24175 GeV 1.40907 6.5611e-07 Gy 17.5591 layer 3: 2.4 5.29862 GeV 1.42436 6.63229e-07 Gy 17.7496 layer 4: 3.2 5.41481 GeV 1.45559 6.77772e-07 Gy 18.1388 layer 5: 4 5.43787 GeV 1.46179 6.80658e-07 Gy 18.216 layer 6: 4.8 5.52101 GeV 1.48414 6.91065e-07 Gy 18.4946 layer 7: 5.6 5.64633 GeV 1.51783 7.06751e-07 Gy 18.9144 layer 8: 6.4 5.63969 GeV 1.51605 7.0592e-07 Gy 18.8921 layer 9: 7.2 5.71744 GeV 1.53695 7.15652e-07 Gy 19.1526 layer 10: 8 5.78086 GeV 1.55399 7.2359e-07 Gy 19.365 layer 11: 8.8 5.94371 GeV 1.59777 7.43975e-07 Gy 19.9106 layer 12: 9.6 6.02518 GeV 1.61967 7.54172e-07 Gy 20.1835 layer 13: 10.4 6.10292 GeV 1.64057 7.63903e-07 Gy 20.4439 layer 14: 11.2 6.18071 GeV 1.66148 7.7364e-07 Gy 20.7045 layer 15: 12 6.2621 GeV 1.68336 7.83827e-07 Gy 20.9771 layer 16: 12.8 6.37762 GeV 1.71441 7.98286e-07 Gy 21.3641 layer 17: 13.6 6.52458 GeV 1.75392 8.16682e-07 Gy 21.8564 layer 18: 14.4 6.67805 GeV 1.79517 8.35892e-07 Gy 22.3705 layer 19: 15.2 6.86252 GeV 1.84476 8.58982e-07 Gy 22.9884 layer 20: 16 6.95226 GeV 1.86889 8.70215e-07 Gy 23.289 layer 21: 16.8 7.11679 GeV 1.91312 8.90809e-07 Gy 23.8402 layer 22: 17.6 7.15125 GeV 1.92238 8.95122e-07 Gy 23.9556 layer 23: 18.4 7.43274 GeV 1.99805 9.30356e-07 Gy 24.8986 layer 24: 19.2 7.58811 GeV 2.03981 9.49804e-07 Gy 25.419 layer 25: 20 7.8156 GeV 2.10097 9.78279e-07 Gy 26.1811 layer 26: 20.8 7.94754 GeV 2.13643 9.94794e-07 Gy 26.6231 layer 27: 21.6 8.35363 GeV 2.2456 1.04562e-06 Gy 27.9834 layer 28: 22.4 8.44564 GeV 2.27033 1.05714e-06 Gy 28.2917 layer 29: 23.2 8.74817 GeV 2.35166 1.09501e-06 Gy 29.3051 layer 30: 24 9.08194 GeV 2.44138 1.13679e-06 Gy 30.4232 layer 31: 24.8 9.50886 GeV 2.55615 1.19022e-06 Gy 31.8533 layer 32: 25.6 9.93302 GeV 2.67017 1.24332e-06 Gy 33.2741 layer 33: 26.4 10.5627 GeV 2.83943 1.32213e-06 Gy 35.3834 layer 34: 27.2 11.1563 GeV 2.99902 1.39644e-06 Gy 37.3721 layer 35: 28 12.0025 GeV 3.22647 1.50235e-06 Gy 40.2065 layer 36: 28.8 13.1124 GeV 3.52485 1.64128e-06 Gy 43.9247 layer 37: 29.6 14.5158 GeV 3.90209 1.81694e-06 Gy 48.6257 layer 38: 30.4 16.8932 GeV 4.54119 2.11452e-06 Gy 56.5898 layer 39: 31.2 20.9827 GeV 5.64051 2.62641e-06 Gy 70.289 layer 40: 32 29.8521 GeV 8.02475 3.73658e-06 Gy 100 layer 41: 32.8 14.922 GeV 4.0113 1.86779e-06 Gy 49.9866 layer 42: 33.6 330.005 MeV 0.088711 4.13067e-08 Gy 1.10547 layer 43: 34.4 2.92265 MeV 0.000785658 3.65828e-10 Gy 0.00979044 layer 44: 35.2 864.413 keV 0.000232369 1.08199e-10 Gy 0.00289566 layer 45: 36 673.958 keV 0.000181172 8.43594e-11 Gy 0.00225766 layer 46: 36.8 1.18851 MeV 0.000319491 1.48765e-10 Gy 0.00398132 layer 47: 37.6 3.81319 MeV 0.00102505 4.77297e-10 Gy 0.0127736 layer 48: 38.4 4.58219 MeV 0.00123177 5.73552e-10 Gy 0.0153496 layer 49: 39.2 7.29449 MeV 0.00196088 9.13052e-10 Gy 0.0244355 layer 50: 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 50 slices. Edep is the deposited energy in every slice. Total incident energy(Ebeam)= 372 GeV Total energy deposit= 367.368 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"> 0.8 17.3315 1.6 17.5591 2.4 17.7496 3.2 18.1388 4 18.216 4.8 18.4946 5.6 18.9144 6.4 18.8921 7.2 19.1526 8 19.365 8.8 19.9106 9.6 20.1835 10.4 20.4439 11.2 20.7045 12 20.9771 12.8 21.3641 13.6 21.8564 14.4 22.3705 15.2 22.9884 16 23.289 16.8 23.8402 17.6 23.9556 18.4 24.8986 19.2 25.419 20 26.1811 20.8 26.6231 21.6 27.9834 22.4 28.2917 23.2 29.3051 24 30.4232 24.8 31.8533 25.6 33.2741 26.4 35.3834 27.2 37.3721 28 40.2065 28.8 43.9247 29.6 48.6257 30.4 56.5898 31.2 70.289 32 100 32.8 49.9866 33.6 1.10547 34.4 0.00979044 35.2 0.00289566 36 0.00225766 36.8 0.00398132 37.6 0.0127736 38.4 0.0153496 39.2 0.0244355 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 50 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 50 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_PB1.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_PBzoom1.png You can close your ROOT session by typing <pre style="color: #800000; background-color: #dcdcdc"> .q </pre> Folder '''RootScripts''' contains several ROOT scripts which plot dose deposition in data and simulation. You can use script '''PlotSimulation.C''' to plot the dose deposition along the absorber. This script uses '''PlotDose.txt'''. Copy the script to your current ProtonPB_build directory: <pre style="color: #800000; background-color: #dcdcdc"> cp /home/username/ProtonPBFolder/ProtonPB_source/RootScripts/PlotSimulation.C . </pre> Then run the script in the following way: <pre style="color: #800000; background-color: #dcdcdc"> root -l .x PlotSimulation.C </pre> This will create '''Simulation.root''' file with the following plot: http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/Simulation1.png You can also plot the data in '''DoseLongitudinalMesh.txt''' and '''DoseLateralMesh.txt''' which were created using commands '''/score/''' in the macro '''proton.mac'''. The file '''DoseLongitudinalMesh.txt''' will be used later to compare with data from the Clatterbridge Cancer Center. Now, use script '''PlotLateralDoseMesh.C''' to plot the lateral dose distribution. <pre style="color: #800000; background-color: #dcdcdc"> cp /home/username/ProtonPBFolder/ProtonPB_source/RootScripts/PlotLateralDoseMesh.C . root -l .x PlotLateralDoseMesh.C </pre> This will create '''LateralDose_Mesh.root''' file with the following plot: http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/SimulationLateralMesh.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. 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 4 cm /protonPB/det/sliceNumber 50 # # define longitudinal scoring mesh # along the beam /score/create/boxMesh waterMeshlongitudinal /score/mesh/boxSize 2. 2. 2. cm /score/mesh/nBin 50 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 0.1 2. 2. cm /score/mesh/nBin 1 1 50 /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 # # visualisation #/control/execute visualisation.mac # /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 # # drawing projections #/score/drawProjection waterMeshlongitudinal doseDeposit #/score/drawProjection waterMeshlateral doseDeposit # # 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 50. 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=60. For example, if you want to have 55 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 55 </pre> '''Modify the mesh''' You can change the size of the mesh (longitudinal and lateral) and the number of voxels by modifying their corresponding lines <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> === If you want to use visualisation, in macro '''proton.mac''' uncomment the line '''/control/execute visualisation.mac'''. This will run macro '''visualisation.mac''' with a specific visualisation settup. If you uncomment the lines '''/score/drawProjection waterMeshlongitudinal doseDeposit''' and '''/score/drawProjection waterMeshlateral doseDeposit''' you will draw the dose projections. In this example, we use the '''DAWN''' event display. Before running the visualisation have a look at this [http://geant4.slac.stanford.edu/Presentations/vis/G4DAWNTutorial/G4DAWNTutorial.html DAWN tutorial]. To run the visualisation, first uncomment lines '''/control/execute visualisation.mac''' and '''/score/drawProjection waterMeshlongitudinal doseDeposit'''. Then, run the proton.mac <pre style="color: #800000; background-color: #dcdcdc"> ./protonPB proton.mac </pre> In addition to the text files it will create two .prim files, g4_00.prim and g4_01.prim. The first one with geometry and particle interactions and the second with dose projections. While running the proton.mac it will ask you to open g4_00.prim in DAWN. In the opened window press OK. This will create your first visualisation: http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/g4_00_6000e.eps The second file g4_01.prim will not open automatically. Open it in the following way: <pre style="color: #800000; background-color: #dcdcdc"> dawn g4_01.prim </pre> In the opened window press OK and this will create the image: http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/g4_02_6000e.eps You can modify the .prim files in DAWN. For example, change the polar and azimuthal angles: (polar angle, azimuthal angle) = (0,90) will create this image http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/g4_01_6000e.eps (polar angle, azimuthal angle) = (90,0) will create this image http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/g4_04_6000e.eps Now, in proton.mac macro uncomment '''/score/drawProjection waterMeshlateral doseDeposit''' and comment '''/score/drawProjection waterMeshlongitudinal doseDeposit'''. Run the proton.mac with the new settings. The file with the lateral dose projections will look like that: http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/g4_03_6000e.eps == <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 script '''PlotData.C''' from folder'''RootScripts''' to plot the data. Copy files '''ClatterbridgeData.txt''' and '''PlotData.C''' to your current ProtonPB_build directory and run the script: <pre style="color: #800000; background-color: #dcdcdc"> cp /home/username/ProtonPBFolder/ProtonPB_source/RootScripts/PlotData.C . cp /home/username/ProtonPBFolder/ProtonPB_source/RootScripts/ClatterbridgeData.txt . root -l .x PlotData.C </pre> This will create '''ClatterbridgeData.root''' file with the following plot: http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/Clatterbridge/Clatterbridge.png == <span style="color:#000080"> Comparison between data and simulation </span> == The scripts '''PlotDataAndSim.C''' and '''PlotDataAndSimMesh.C''' in folder '''RootScripts''' compare data with simulation. You can use script '''PlotDataAndSim.C''' to compare data (ClatterbridgeData.txt) and simulation ('''PlotDose.txt'''). Both text files must be in the folder where you run the script. <pre style="color: #800000; background-color: #dcdcdc"> cp /home/username/ProtonPBFolder/ProtonPB_source/RootScripts/PlotDataAndSim.C . root -l .x PlotDataAndSim.C </pre> This will create '''BraggPeakComparison.root''' file with the following plot: http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/DataSimulation.png You can also compare data (ClatterbridgeData.txt) with simulation done with scoring mesh ('''DoseLongitudinalMesh.txt'''). This can be done with script '''PlotDataAndSimMesh.C'''. This script works only if before running it you substitude the commas in '''DoseLongitudinalMesh.txt''' with spaces. Remove also the header in the text file. <pre style="color: #800000; background-color: #dcdcdc"> cp /home/username/ProtonPBFolder/ProtonPB_source/RootScripts/PlotDataAndSimMesh.C . root -l .x PlotDataAndSimMesh.C </pre> This creates '''BraggPeakComparison_Mesh.root''' file with the following plot: http://www.hep.ucl.ac.uk/pbt/RadiotherapyWorkbook/skins/common/images/ProtonPB/DataSimulationMesh.png == <span style="color:#000080"> Files </span> == [[List of monoenergetic proton pencil beam files with brief description]]
Return to
Monoenergetic proton pencil beam
.
Views
Page
Discussion
View source
History
Personal tools
18.117.255.197
Talk for this IP address
Log in
Navigation
Main page
Community portal
Current events
Recent changes
Random page
Help
Search
Toolbox
What links here
Related changes
Special pages