Clatterbridge: Difference between revisions

From PBTWiki
Jump to navigation Jump to search
No edit summary
 
(99 intermediate revisions by 2 users not shown)
Line 1: Line 1:
== Simulation of the Clatterbridge beamline ==


This simulation is a model of the monoenergetic 62.5 MeV proton beam at the [http://www.clatterbridgecc.nhs.uk/patients/treatment-and-support/proton-therapy Clatterbridge Cancer Centre] as it traverses the components of the beamline and finally hits a volume of water. The simulation was built using the existing example "[http://www.hep.ucl.ac.uk/pbt/wiki/Software/Geant4/Tutorials/Basic/Monoenergetic_Proton_Pencil_Beam monoenergetic proton pencil beam] ". The physics list used is QGSP_BIC_EMY, the standard for simulating clinical proton beams.
This simulation is a model of the monoenergetic 62.5 MeV proton beam at the [http://www.clatterbridgecc.nhs.uk/patients/treatment-and-support/proton-therapy Clatterbridge Cancer Centre] as it traverses the components of the beamline and finally hits a volume of water.  
 
The simulation was built using:
=== Geometry ===
* [[/Geant4 | Geant4]]
==== Treatment room ====
** The simulation was built on the example in <code>examples/extended/electromagnetic/TestEm7</code> supplied with the Geant4 package and detailed [http://www.hep.ucl.ac.uk/pbt/wiki/Software/Geant4/Tutorials/Basic/Monoenergetic_Proton_Pencil_Beam here].
The schematic below shows the Clatterbridge treatment room, illustrating the layout of the geometries defined in the [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/src/DetectorConstruction.cc  <code>DetectorConstruction</code>] class. The whole beamline is contained within the "inner room" volume, which is placed in the "outer room" volume. As the outer volume is a solid concrete box and the inner volume is a box filled with air, the small overlap models the concrete walls of the treatment room. The proton source is placed against the wall in the inner room (-4200 mm from the origin), and all beamline components are placed relative to this reference point.
* [[/TOPAS | TOPAS]]
 
 
[[File:treatmentRoomSchematic.png|centre|700px]]
 
==== Beamline ====
[[Media:Beamline_doc.pdf|The positions of the components relative to the source are shown here (annotated PDF)]]. The beamline is shown below. The top figure is a perspective view. The second figure is a top-down view which also shows the borated plastic shielding (f) that was previously hidden to allow the full beamline to be visible.
 
[[File:beamline_perspective.png|border|right|600px]]
[[File:beamline_top_down.png|border|right|600px]]
 
The beamline consists of the following components:
 
*an aluminium tube (a) containing from left to right;
**a brass collimator
**a tungsten scattering foil
**a brass beam stopper
**another tungsten scattering foil
**a mylar window
 
*an empty aluminium box (b)
*an iron block (c)
*a second aluminium tube (d)
*a second aluminium box (e) containing from left to right;
**a brass collimator
**two dose monitors
**cross wires
 
*a brass nozzle (g)
 
The water phantom (h) is also shown for completeness.
 
==== Dual scattering tube ====
[[Media:Tube_doc.pdf|The positions of the components in the tube are shown here (annotated PDF)]]. The proton beam at Clatterbridge is delivered using passive spreading. That is, the beam is spread out using dual scattering where the proton beam is scattered through a first scattering foil followed by a beam stopper and a second scattering foil. The beam stopper is used to reduce the intensity of the beam at its centre and a high-Z material such as tungsten is chosen for the foils as it is highly scattering.
 
A visualisation of the first tube with 500 primary protons represented by blue tracks is shown below. The protons first travel through the collimator, followed by the first scattering foil which spreads out the beam. The brass stopper then blocks out approximately half of the remaining protons before the beam is again spread out through the second scattering foil. The intention is to produce a wide, homogeneous beam. The particles exit the first tube through a Kapton window which is used to keep the tube under vacuum to reduce random scattering off air molecules.
 
[[File:First_tube_500_sim.png|centre|700px]]
 
==== Dosimetry box ====
[[Media:Dosimetry_box.pdf|The positions of the components in the dosimetry box are shown here (annotated PDF)]]. The aluminium box shown below contains the dosimetry equipment. The beam is first collimated using a brass collimator to remove any stray protons. The beam then travels through the dose monitors, the cross wires and finally exits the box through the brass nozzle.
 
[[File:Second_box.png|centre|700px]]
 
[[Media:Dose_monitor_doc.pdf|The positions of the components in the dose monitors are shown here (annotated PDF)]]. The dose monitors are drift chambers used to measure the dose deposited by the proton beam. An exploded view of a dose monitor as used in the dosimetry box is shown below. It consists of a set of aluminised mylar foils wedged between layers of perspex to hold them in place. The guard ring is used to create a sealed volume of air between the foils. The aluminium layers face towards the centre of the dose monitor such that the assembly acts as a drift chamber when a potential difference is applied.
 
 
[[File:Dose_monitor_exploded.png|centre|700px]]
 
==== Visualisation ====
The geometry and particle tracks in Geant4 [http://geant4.slac.stanford.edu/Presentations/vis/G4DAWNTutorial/G4DAWNTutorial.html can be visualised using the DAWN visualiser]. Other, perhaps more suitable, visualisers are available (see [http://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch08.html]) but DAWN is useful in that it produces very high-resolution visualisations. The main drawback of DAWN is that it does not allow for click-and-drag functionality. Instead, parameters must be set before each visualisation is drawn which tends to require a fair amount of trial-and-error to obtain good images. All images of the geometry on this page were created using DAWN and an example visualisation of the beamline with all particle tracks drawn can be seen below. The first tube and second box are visualised using the wireframe setting so that their inner components are visible. This example contains 500 primary protons. Protons are shown in blue, electrons are red, positrons are cyan, gamma rays are green, and neutrons are yellow.
 
[[File:Beamline_500_sim_perspective.png|border|centre|800px]]
 
 
=== Generation of primary protons ===
 
The protons are generated using the [http://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch02s07.html <code>G4GeneralParticleSource</code>] (GPS) class which allows for a range of properties of the primary protons to be set in the [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/proton.mac <code>proton.mac</code>] file (primaries are the initial particles generated by the GPS). First, the particle source is positioned against the wall on one side of the room which is at -4200 mm from the centre. In an attempt to achieve a more realistic beam, the primary protons are distributed normally in the x and y directions centred at 0 with standard deviations of 4.0 mm and 4.5 mm respectively. This gives the beam width and results in a nearly circular profile. The primaries are generated with initial energies following a Gaussian distribution with mean 62.5 MeV and standard deviation 0.082 MeV, and Gaussian radial distributions with respect to the z-axis with standard deviations of 2.3 mrad and 1.2 mrad in the x and y directions respectively.
 
The parameters are set in the following section of the [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/proton.mac <code>proton.mac</code>] macro:
 
<pre>
### particle gun settings
/gps/particle proton
/gps/number 1 # 1 proton per event
 
# energy
/gps/ene/type Gauss
/gps/ene/mono 62.5 MeV
/gps/ene/sigma 0.082 MeV
 
# source position
/gps/position 0.0 0.0 -4200 mm
 
# beam width
/gps/pos/type Beam
/gps/pos/sigma_x 4.0 mm
/gps/pos/sigma_y 4.5 mm
 
# angular distribution of primaries for realistic emittance
/gps/ang/type beam2d
/gps/ang/sigma_x 2.3 mrad
/gps/ang/sigma_y 1.2 mrad
/gps/ang/rot1 -1 0 0 # aligns gps with positive z-axis
</pre>
 
=== Gathering data ===
==== Tracking ====
The energy of the beam is monitored at several stages in the beamline through the [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/src/SteppingAction.cc <code>SteppingAction.cc</code>] class. If the protons traverse a position in z defined in the vector <code>trackingDistFiller</code> in the [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/include/SteppingAction.hh <code>SteppingAction.hh</code>] header, their parent ID, xyz-coordinates, orientation with respect to the z-axis and energy are recorded at the given position in z. The "sensitive area" is restricted to a square of 100 mm in each direction in the xy-plane.
 
==== Scoring ====
===== Longitudinal scoring mesh =====
A scoring mesh is "longitudinal" in the sense that it records data along the direction of the beamline, e.g. the energy deposited at different positions in z. Hence, for a mesh to be longitudinal it must have bins defined along its z-axis. The scorer defined on the mesh <code>detEnergyLon</code> is an example of such a mesh. The mesh utilises a filter to detect the flux of protons per cm<sup>2</sup> and writes the data to the text file <code> FluxLongitudinal.txt </code>. The location of the mesh centre can be changed in <code> proton.mac </code>, in addition to the dimensions of the mesh and the number of bins.
 
<pre>
/score/create/boxMesh waterMeshlongitudinal
/score/mesh/boxSize 10. 10. 10. cm
/score/mesh/nBin 1 1 400
/score/mesh/translate/xyz 0. 0. -226 cm
</pre>
 
The filter can also be changed to observe the flux of particles other than protons:
 
<pre>
/score/quantity/cellFlux protonFlux
/score/filter/particle protonFilter proton
</pre>
 
==== Lateral scoring mesh ====
 
A lateral scoring mesh is positioned at the end of the nozzle to record the energy distribution of the protons perpendicular to the direction of the beam. The position, size and bin number of this mesh can be modified in the same way as the longitudinal mesh example above.
 
== Running simulations ==
 
=== Make the simulation ===
To access a Linux desktop PC running Scientific Linux 6, [http://www.hep.ucl.ac.uk/pbt/wiki/Software/Geant4/UCL_HEP_Cluster follow the instructions here]. Once access has been established, follow the instructions below to copy the simulation files to your directory:
 
<pre>
[username@pc1XX ~]$ mkdir sim
[username@pc1XX ~]$ cd sim
[username@pc1XX sim]$ cp -r /unix/pbt/mhentz/Clatterbridge_sim/source .
[username@pc1XX sim]$ mkdir build
[username@pc1XX sim]$ cd build
[username@pc1XX build]$ /unix/pbt/software/scripts/pbt.sh
[username@pc1XX build]$ cmake -DGeant4_DIR=/unix/pbt/software/dev ../source
[username@pc1XX build]$ make
</pre>
 
=== Running the simulation in batch mode with macro ===
 
This will run the simulation and produce the required output files in the <code>data</code> subdirectory created.
 
<pre>
[username@pcXXX build]$ mkdir data
[username@pcXXX build]$ ./clatterbridgeSim proton.mac
</pre>
 
=== Running large simulations ===
A simulation comprising of a large number of particles can take a long time to run as events run sequentially (about 20 hours for 1e6 events). This can be reduced considerably by running a number of independent simulations with fewer events "in parallel" and then merging the output files once they have completed.  The main limitation lies in the maximum number of simulations that can be run simultaneously–which is 120. Splitting 1e6 events into 100 simulations of 10,000 events reduces the runtime to a much more bearable 20 minutes.
 
For split simulations, the parameters are set in the macro [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/proton.mac_split <code>proton.mac_split</code>]. The main difference to [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/proton.mac <code>proton.mac</code>] is the lower number of events and the variable seeds
<pre>
/random/setSeeds ${seed1} ${seed2}
</pre>
which are set in the PBS job submission script <code>clatterbridge_X.sh</code> described below (X is an integer). This ensures that all simulations are independent (since they have different seeds) but as the same set of seeds is used each time, the simulations remain repeatable.
 
The script [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/run.sh <code>run.sh</code>] (located in the <code>source</code> directory but should automatically be copied to <code>build</code> by cmake) creates a timestamped directory along with a few useful subdirectories and submits a specified number of simulations to be executed there. It is currently set to submit 100 jobs to the medium queue through the [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/submit.sh <code>submit.sh</code>] script. This can be changed by supplying different arguments on line 22 (other queues available are <code>short</code> and <code>long</code>):
<pre>
../submit.sh 100 medium
</pre>
 
The script [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/submit.sh <code>submit.sh</code>] assigns values to the variables in the template [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/clatterbridge_split.sh <code>clatterbridge_split.sh</code>] and writes the result to a job submission script <code>clatterbridge_X.sh</code> for all 100 simulations. This is done on line 31:
<pre>
cat ../clatterbridge_split.sh | sed -e "s/\${i}/$i/" >> clatterbridge_$i.sh
</pre>
 
Besides setting vital parameters and navigating to the correct folder, the submission script sets the seeds in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/proton.mac_split <code>proton.mac_split</code>] and writes the result to a new macro <code>proton.mac_simX</code>. It then runs the simulation with the corresponding macro.
 
Follow the instructions in "[http://www.hep.ucl.ac.uk/pbt/pbtWiki/index.php?title=Clatterbridge&action=submit#Make_the_simulation Make the simulation] " unless you have already done so. Split simulations are run as follows:
<pre>
[username@pc1XX build]$ cp ../source/clatterbridge_split.sh .
[username@pc1XX build]$ ./run.sh
</pre>
 
The status of the simulations can be tracked with the following command (can be used in conjunction with the <code>watch</code> command):
<pre>
[username@pc1XX build]$ qstat -u <username>
</pre>
 
If something goes wrong and you wish to stop all the jobs simultaneously you can use the command
<pre>
[username@pc1XX build]$ qselect -u <username> | xargs qdel
</pre>
 
The resulting data can be combined by running the scripts [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/combine_tracking.py <code>combine_tracking.py</code>] and [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/combine_scorers.py <code>combine_scorers.py</code>] as follows:
<pre>
[username@pc1XX build]$ cd timestamp-dir
[username@pc1XX timestamp-dir]$ cp ../../source/combine_*.py
[username@pc1XX timestamp-dir]$ ./combine_tracking.py
[username@pc1XX timestamp-dir]$ ./combine_scorers.py
</pre>
 
== Output files ==
=== Output used for beam emittance and energy spectra ===
The output files <code>output_zX.txt</code> are written to the subdirectory <code>data/</code> in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/src/SteppingAction.cc <code>SteppingAction.cc</code>] (<code>X</code> is the position in z where the protons were recorded). Other particles can be recorded by changing [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/src/SteppingAction.cc <code>SteppingAction.cc</code>] accordingly. The output files contain the parent ID which is 0 for a primary particle and non-zero for all other particles, the xyz-coordinates in mm, the angle subtended by the projection of the momentum vector onto the x-axis and the z-axis in mrad (used to plot the beam emittance) and the kinetic energy in MeV. The z-coordinates are given relative to the position of the particle source and the z-axis lies along the beamline with its positive end pointing towards the origin. An example is shown below (<code>output_z1759.txt</code>):
 
<pre>
# zposition: 1759
# parentID, x, y, z [mm], theta [mrad], ke [MeV]
0 -9.78358 4.39357 1759.0 -6.01461 60.32370
0 -7.12684 -10.49640 1758.9 -4.95641 60.12080
0 8.83137 -2.33413 1759.0 10.60610 60.07960
0 5.92333 -3.29109 1758.9 7.89279 60.20930
0 -2.21463 6.24370 1759.0 -6.14550 60.04930
...
</pre>
 
=== Scorer outputs ===
The files described below are all written using the <code>/score/dumpQuantityToFile</code> commands in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/score_dump.mac <code>score_dump.mac</code>]. The data are recorded by the scorers defined in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/score_init.mac <code>score_init.mac</code>].
 
==== Lateral energy deposition at the Bragg peak in the detector ====
The file <code>BraggEnergyDepLat.txt</code> contains data recorded by the <code>energyDeposit</code> scorer defined on the <code>braggEnergyLat</code> mesh in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/score_init.mac <code>score_init.mac</code>]. It records the energy deposition at the Bragg peak in the plane perpendicular to the direction of incidence of the beam. The columns represent the bin number in the x, y and z directions and the energy deposition in MeV (as the mesh is only binned in x, it consists of strips parallel to the y-axis):
<pre>
# mesh name: braggEnergyLat
# primitive scorer name: energyDeposit
# iX, iY, iZ, value [MeV]
0,0,0,29.931796279684590
1,0,0,38.256791245912446
2,0,0,115.051837974658625
3,0,0,185.831154261095321
...
</pre>
This file can be used to plot the lateral energy deposition profile at the Bragg peak.
 
==== Longitudinal energy deposition throughout the detector ====
The file <code>DetEnergyDepLon.txt</code> contains data recorded by the <code>energyDeposit</code> scorer defined on the <code>detEnergyLon</code> mesh in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/score_init.mac <code>score_init.mac</code>]. It records the energy deposition in the detector in bins along the beam. The columns represent the bin number in the x, y and z directions and the energy deposition in MeV:
<pre>
# mesh name: detEnergyLon
# primitive scorer name: energyDeposit
# iX, iY, iZ, value [MeV]
0,0,0,59390.147050145809771
0,0,1,59016.142052385039278
0,0,2,59209.373600437116693
0,0,3,59578.920521325460868
...
</pre>
This file can be used to plot the energy deposition profile through the detector. The Bragg peak can be identified from such a plot.
 
== Data Analysis ==
=== Beam profile, emittance and energy spectrum ===
The script <code>analysis.py</code> produces several different plots from the <code>output_zX.txt</code> files. It plots the beam profile, the projection of the beam profile onto the x-axis, the beam emittance and the energy spectrum at the position in z where the protons were recorded. This script requires the <code>SciPy</code> module which is not available on the Linux SL6 PCs. Hence, it must currently be run locally. It is run for a given position along the beamline as follows (400 mm as an example, the script must be run in the directory containing the <code>data</code> subdirectory):
 
<pre>
[user@hostname ~]$ ./analysis.py 400
</pre>
 
This runs the script using the file <code>output_z400.txt</code> and produces the corresponding plots shown below. At this stage, the beam has undergone dual scattering and is starting to spread out again to produce a uniform beam.
 
 
[[File:Tiles_400.jpg|centre|750px]]
 
 
=== Longitudinal energy deposition profile and Bragg peak ===
The script <code>bragg.py</code> plots the longitudinal energy deposition profile through the detector from <code>DetEnergyDepLon.txt</code>. It is executed using Python 2.7 in the same directory as the data file as follows:
 
<pre>
[user@hostname ~]$ ./bragg.py
</pre>
 
This produces the following plot;
[[File:Lon_energy_deposition_bragg.png|centre|600px]]
 
=== Lateral energy deposition profile at the Bragg peak ===
The script <code>lateral.py</code> plots the lateral energy deposition at the Bragg peak from <code>BraggEnergyDepLat.txt</code> described above. It is executed on Python 2.7 as follows:
 
<pre>
[user@hostname ~]$ ./lateral.py BraggEnergyDepLat.txt
</pre>
 
This produces the following plot;
[[File:Lat_energy_deposition_bragg.png|centre|600px]]
 
=== Using ROOT ===
These instructions and the script [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/ProtonPB_source/simulation_analysis.C <code>simulation_analysis.C</code>] were written by the predecessor working on this simulation. '''If you happen to perform your data analysis in ROOT, please share your wisdom below.''' The script requires output files for certain plots that are currently not being generated to save computation time. For instance, to plot the flux along the beamline, another scoring mesh needs to be defined in <code>score_init.mac</code>:
 
<pre>
/score/create/boxMesh beamlineMeshLongitudinal
/score/mesh/boxSize 50 50 1000 mm
/score/mesh/nBin 1 1 400
/score/mesh/translate/xyz 0 0 -3200 mm
/score/quantity/flatSurfaceFlux protonFlux
/score/filter/particle protonFilter proton
/score/quantity/flatSurfaceFlux neutronFlux
/score/filter/particle neutronFilter neutron
/score/close
</pre>
 
==== Open ROOT and run analysis file ====
 
The simulation analysis file reads the data in the output files and produces the associated plots.
 
<pre>
[username@plus1 ProtonPB_build]$ root -l
 
root [0] .x simulation_analysis.C
</pre>
 
==== Proton energy deposition in water ====
 
On the left is a histogram showing the frequency of incident protons to the water box at different energy values. It is produced from data in <code> stopper.txt </code> which contains the hits data from the sensitive detectors in the simulation, including the water volume sensitive detector which is used here. The peak in proton energy deposited in the water volume is found to be at 60.08 MeV.
The plot on the right shows the cumulative proton energy deposited within the water volume, plotting data collected in the longitudinal water volume scoring mesh, showing the Bragg peak at 31 mm.
 
[[File:WaterEnergyDeposition.png|400px]]
[[File:Cumulative energy deposition in water.png|400px]]
 
==== Proton stopping distance in water ====
This plot shows the stopping distance of protons in the water volume, with the stopping distance defined as the point at which the velocity of the protons is zero.
[[File:StoppingDistance.png|centre|400px]]
 
==== Proton flux along beamline ====
This plot is produced from <code> FluxLongitudinal.txt</code> and shows the number of protons per cm<sup>2</sup> at given points along the beamline (in the z direction). The decrease in proton flux is clearly shown at particular points along the beamline, for example at the location of the brass stopper (30cm).
[[File:ProtonFluxBeamlinePlot.png|centre|700px]]
 
==== Proton flux along water volume ====
Using a similar setup as the previous beamline flux plot, proton flux data from the scoring mesh within the water volume was plotted.
 
[[File:ProtonFluxWaterVolume.png|centre|500px]]
 
==== Kinetic energy of the proton beam ====
 
This left-hand plot is produced from <code> kin.out </code> and shows the kinetic energy distribution of the proton beam at a specified point along the beamline. This location is specified in <code> SteppingAction.cc </code> in which a particular position in z can be input (see Section 5.6 Kinetic Energy Readings).
 
The xy coordinates of each proton are also recorded, and are used to produce a plot which shows the spatial coordinates and the kinetic energy of each proton at the specified position in z, along the beamline.
 
The z position along the beamline chosen to produce the following plots was just after the nozzle.
[[File:ProtonKineticEnergy.png|centre|700px]]
 
==== Neutron flux along beamline and kinetic energy ====
 
This was plotted from data in the file <code>NeutronFluxBeamline.txt</code>, and shows the increase in flux of neutrons at certain points (locations of beamline components) along the beamline. Data on the energy of neutrons produced within a particular volume along the beamline was collected, with the volume specified in <code>SteppingAction.cc</code> - the plot below shows the energy of neutrons produced within the water volume.
 
[[File:NeutronFluxBeamline.png|left|400px]]
[[File:NeutronKineticEnergy.png|centre|400px]]
 
== Changing parameters ==
 
=== Beam parameters ===
These can be modified in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge/source/proton.mac <code>proton.mac</code>] under the header "particle gun settings".
 
==== Energy ====
The GPS gives the initial protons a range of energies following a Gaussian distribution.
 
<pre>
/gps/ene/type Gauss
/gps/ene/mono 62.5 MeV
/gps/ene/sigma 0.082 MeV
</pre>
 
==== Source position ====
The proton source is positioned at <code>z = -4200 mm</code> relative to the centre of the inner room, which corresponds to a position against the wall of the inner room.
 
<pre>
/gps/position 0.0 0.0 -4200 mm
</pre>
 
=== Scoring meshes ===
==== Longitudinal scoring mesh ====
A scoring mesh is "longitudinal" in the sense that it records data along the direction of the beamline, e.g. the energy deposited at different positions in z. Hence, for a mesh to be longitudinal it must have bins defined along its z-axis. The scorer defined on the mesh <code>detEnergyLon</code> is an example of such a mesh. The mesh utilises a filter to detect the flux of protons per cm<sup>2</sup> and writes the data to the text file <code> FluxLongitudinal.txt </code>. The location of the mesh centre can be changed in <code> proton.mac </code>, in addition to the dimensions of the mesh and the number of bins.
 
<pre>
/score/create/boxMesh waterMeshlongitudinal
/score/mesh/boxSize 10. 10. 10. cm
/score/mesh/nBin 1 1 400
/score/mesh/translate/xyz 0. 0. -226 cm
</pre>
 
The filter can also be changed to observe the flux of particles other than protons:
 
<pre>
/score/quantity/cellFlux protonFlux
/score/filter/particle protonFilter proton
</pre>
 
==== Lateral scoring mesh ====
 
A lateral scoring mesh is positioned at the end of the nozzle to record the energy distribution of the protons perpendicular to the direction of the beam. The position, size and bin number of this mesh can be modified in the same way as the longitudinal mesh example above.
 
=== Beamline components ===
 
Components of the beamline can be added/removed in <code> DetectorConstruction.cc </code>. If the dimensions of the water box are modified, the following lines in <code> simulation_analysis.C </code> will also need to be modified:
 
<pre>
Float_t lengthBox = 200, widthBox = 200;
</pre>
 
The <code> depthFix </code> variable to calculate the stopping distance of the protons within the analysis will also need to be adjusted if the water box dimensions or location are modified. <code> depthFix </code> is calculated by taking the z position of the centre of the water box (relative to the centre of the inner room), and subtracting the half length of the water box (calculated in mm):
 
<pre>
//depth fix - water box centred at -2260 mm, half length = 100 mm.  
Double_t depthFix = 2360;
</pre>
 
=== Sensitive detectors ===
 
In this simulation, the water volume is assigned as a sensitive detector in <code> DetectorConstruction.cc </code>:
 
<pre>
G4SDManager* SDman = G4SDManager::GetSDMpointer();
G4String name="SD";
DetectorSD = new SensitiveDetector(name);
SDman->AddNewDetector(DetectorSD);
logicWater->SetSensitiveDetector(DetectorSD);
</pre>
 
Another beamline component may be used by modifying the following line and setting its logical volume as a sensitive detector:
 
<pre>
logicWater->SetSensitiveDetector(DetectorSD);
</pre>
 
The component should be chosen such that a significant proportion of the proton beam deposits energy, such as in the brass stopper, in order to produce enough data for plots.
 
<code> SensitiveDetector.cc </code> is derived from the <code> G4VSensitiveDetector.cc </code> base class. On a step-by-step basis, the energy deposited by the proton is recorded as a "hit" and added to a <code> HitsCollection </code> object. Other parameters may be retrieved at each step in the method <code> ProcessHits </code>.
 
<pre>
G4bool SensitiveDetector::ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist)
{
  G4double edep = aStep->GetTotalEnergyDeposit();
  if(aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton"){
    ::Hit* newHit = new ::Hit();
    newHit->SetEdep(edep);
    HitID = detectorCollection->insert(newHit);
    return true;
  }
</pre>
 
=== Physics list ===
 
This simulation uses one of the standard physics lists: QGSP_BIC_HP. A user-defined physics list can be specified in the class <code> PhysicsList.cc </code>, and then selected in <code> protonPB.cc</code>
 
<pre>
runManager->SetUserInitialization(new  PhysicsList);
</pre>
 
=== Kinetic energy readings ===
 
The kinetic energy of the beam after passing through a particular component of the beamline is found in the class <code> SteppingAction.cc </code>. This class contains a map named <code> componentMap </code> with the names of the components as keys and the z coordinate of each component (taken as the furthest edge of the component from the source):
 
<pre>
componentMap.insert( std::pair<std::string,double>("ScatteringFoil1",80.025));
componentMap.insert( std::pair<std::string,double>("BrassStopper",306.6));
componentMap.insert( std::pair<std::string,double>("ScatteringFoil2",306.625));
</pre>
 
To find the kinetic energy distribution of the beam at a particular point, insert the component name from the map above into this line:
 
<pre>
double cPos = (double)componentMap.find("ScatteringFoil2")->second;
</pre>
 
The code then compares the z position of a proton with the position of the chosen component on a step-by-step basis, and if the proton is at the chosen position its kinetic energy is recorded in <code> kin.out </code>:
 
<pre>
if(fabs(zpos-cPos) < 1.0e-8){
  std::ofstream hitsfile(filename, std::ios::app);
  if(hitsfile.is_open()){
      G4double ke =  step->GetPreStepPoint()->GetKineticEnergy();
      hitsfile << zpos << "\t"
      << ke <<G4endl;
      hitsfile.close();
  }
}
</pre>
 
 
'''After any modifications to the simulation files, the code will need to be compiled. In the build directory, write:'''
 
<pre>
[username@plus1 ProtonPB_build]$ make
</pre>
 
After this <code> proton.mac </code> can be run:
 
<pre>
[username@plus1 ProtonPB_build]$ ./protonPB proton.mac
</pre>
 
== Modifying Analysis Methods ==
 
The number of bins to store the data for each histogram can be changed, here the number of bins is 800:
 
<pre>
TH1D* waterEnergy = new TH1D ("waterDep", "Energy deposition in water box", 800, 0., 65.);
</pre>
 
The range of the axes for each plot can also be specified in order to zoom into the energy peak:
 
<pre>
waterEnergy->GetXaxis()->SetRangeUser(0, 64);
</pre>
 
 
 
== Files ==
 
[[Clatterbridge files|List of files for Clatterbridge simulation]].

Latest revision as of 07:46, 14 August 2020

This simulation is a model of the monoenergetic 62.5 MeV proton beam at the Clatterbridge Cancer Centre as it traverses the components of the beamline and finally hits a volume of water. The simulation was built using:

  • Geant4
    • The simulation was built on the example in examples/extended/electromagnetic/TestEm7 supplied with the Geant4 package and detailed here.
  • TOPAS
Retrieved from "https://www.hep.ucl.ac.uk/pbt/pbtWiki/index.php?title=Clatterbridge&oldid=3224"