Clatterbridge
Simulation of the Clatterbridge beamline
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 beamline components are contained within a geometry modelling the Clatterbridge treatment room. The simulation was built using the existing example "monoenergetic proton pencil beam" [1]. The physics list used is QGSP_BIC_EMY, the standard for simulating clinical proton beams.
Geometry
The schematic below shows the Clatterbridge treatment room, illustrating the layout of the geometries defined in the DetectorConstruction 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.
The beamline shown on the right consists of the following components:
- an aluminium tube (a) under vacuum 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 is also shown (h) for completeness.
Generation of primary protons
The protons are generated using the G4GeneralParticleSource (GPS) class which allows for a range of properties of the primary protons to be set in the proton.mac 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 proton.mac macro:
### 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
The energy of the beam is monitored at several stages in the beamline through the SteppingAction class. If the protons traverse a position in z defined in the vector trackingDistFiller in the SteppingAction header, their parent ID, xyz-coordinates, orientation with respect to the z-axis and energy are recorded at the given position in z.
The figure to the right shows the visualisation output of the simulation, displaying the particle tracks. Protons are shown in blue, electrons in red, gammas in green and neutrons in yellow.
Running the simulation
Follow the instructions below to copy the simulation files to your directory:
[username@plus1 ~]$ mkdir Clatterbridge_sim [username@plus1 ~]$ cd Clatterbridge_sim [username@plus1 Clatterbridge_sim]$ /unix/pbt/software/scripts/pbt.sh [username@plus1 Clatterbridge_sim]$ cp -r /unix/pbt/rstephens/ProtonPBFolder/ProtonPB_source . [username@plus1 Clatterbridge_sim]$ mkdir ProtonPB_build [username@plus1 Clatterbridge_sim]$ cd ProtonPB_build [username@plus1 ProtonPB_build]$ make -DGeant4_DIR=/unix/pbt/software/dev /home/username/Clatterbridge_sim/ProtonPB_source [username@plus1 ProtonPB_build]$ make
Run macro proton.mac
This will run the simulation and produce the required output files.
[username@plus1 ProtonPB_build]$ ./protonPB proton.mac
Output files
The simulation code and proton.mac produce several output files:
protonKE.out
This text file contains output information from SteppingAction.cc , printed on a step-by-step basis for each proton (event). The first two columns contain the x and y coordinates (mm) of the proton. The third column contains the z position of the proton (mm), relative to the position of the source at the inner room boundary. The last column contains the energy (MeV) of the proton at this z position.
6.66882 -11.4767 1770 59.8738 -7.60842 -11.4391 1770 60.186 -9.37841 -14.363 1770 43.1861 -13.2143 1.55943 1770 60.4776 1.81337 11.558 1770 60.2989
secondariesKE.out
This text file contains information recorded in SteppingAction.cc about neutrons produced in the simulation. The values shown in the file are the kinetic energies [MeV] of the neutrons produced within the specified volume in the code.
5.12058 2.71641
sensitiveDetectors.out
This file contains the ID of the sensitive detector (first column) and the energy deposited in that detector (second column). The ID number "0" corresponds to the water volume sensitive detector, and ID number "1" corresponds to the user-defined sensitive detector.
1 62.2594 0 43.1061 1 0 0 0 1 62.2851
ProtonFluxBeamline.txt
This file is automatically written by proton.mac and contains the proton flux data from the longitudinal scorer. The columns represent the bin number in the x, y and z directions, and the number of protons per cm2:
# mesh name: waterMeshlongitudinal # primitive scorer name: protonFlux # iX, iY, iZ, value [percm2] 0,0,0,25.97935669321147 0,0,1,25.90015344712939 0,0,2,25.84461190407653
In simulation_analysis.C the data in ProtonFluxBeamline.txt is read in and written to another text file, ProtonFluxBeamline_Mod.txt to modify the standard formatting into a suitable format for further analysis:
0 0 0 25.97935669321147 0 0 1 25.90015344712939 0 0 2 25.84461190407653
A similar output file, NeutronFluxBeamline.txt is produced in the simulation, and contains information about the flux of neutrons along the entire length of the beamline.
An additional two text files are produced from another longitudinal scoring mesh, located specifically along the length of the water volume. ProtonFluxWaterLongitudinal.txt contains information about the proton flux within the water volume, and ProtonEnergyWaterLongitudinal.txt contains information on the cumulative energy deposition of the proton beam within the water volume.
ProtonEnergyWaterLongitudinal.txt
This text file is produced by proton.mac in the same way as ProtonFluxBeamline.txt above, and contains information about the proton energy deposited within the water volume, divided into bins parallel to the beam direction:
# mesh name: waterMeshlongitudinal # primitive scorer name: energyDeposit # iX, iY, iZ, value [MeV] 0,0,0,14778.66421571247 0,0,1,14787.50686991539 0,0,2,14943.0041959171
EnergyLateralMesh.txt
This text file is produced by proton.mac in the same way as the other scorers above, and contains information about the proton energy deposition in bins perpendicular to the beam:
# mesh name: waterMeshlateral # primitive scorer name: energyDeposit # iX, iY, iZ, value [MeV] 0,0,0,0 1,0,0,0 ... 16,0,0,0.002151585502782383 17,0,0,0.01565800625140965 18,0,0,9.09434626226223 19,0,0,0.01811169052807143
Data Analysis
Open ROOT and run analysis file
The simulation analysis file reads the data in the output files and produces the associated plots.
[username@plus1 ProtonPB_build]$ root -l root [0] .x simulation_analysis.C
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 stopper.txt 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.
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.
Proton flux along beamline
This plot is produced from FluxLongitudinal_Mod.txt and shows the number of protons per cm2 passing through distinct 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).
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.
Kinetic energy of the proton beam
This left-hand plot is produced from kin.out and shows the kinetic energy distribution of the proton beam at a specified point along the beamline. This location is specified in SteppingAction.cc 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.
Neutron production
Neutron flux along beamline
This was plotted from data contained in the file NeutronFluxBeamline.txt, and shows the increase in flux of neutrons at certain points (locations of beamline components) along the beamline.
Neutron kinetic energy
Data on the energy of neutrons produced within a particular volume along the beamline was collected, with the volume specified in SteppingAction.cc - the plot below shows the energy of neutrons produced within the water volume.
Changing parameters
Initial beam parameters
Initial parameters of the proton beam can be modified in proton.mac
Beam radius
/gps/pos/radius 3 mm
Beam energy
This simulation models the proton beam source with a Gaussian distribution.
/gps/ene/type Gauss /gps/ene/mono 62.5 MeV /gps/ene/sigma 0.082 MeV
Source position
The proton source is positioned at z = -420 cm relative to the centre of the inner room (the mother volume), which translates as the wall surface of the inner room.
/gps/pos/type Plane /gps/pos/shape Circle /gps/pos/centre 0.0 0.0 -420 cm
Scoring mesh
Longitudinal scoring mesh
A longitudinal scoring mesh extends along the length of the beamline from the source to the water volume. The mesh utilises a filter to detect the flux of protons per cm2 and writes the data to the text file FluxLongitudinal.txt . The location of the mesh centre can be changed in proton.mac , in addition to the dimensions of the mesh and the number of bins.
/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
The filter can also be changed to observe the flux of particles other than protons:
/score/quantity/cellFlux protonFlux /score/filter/particle protonFilter proton
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 DetectorConstruction.cc . If the dimensions of the water box are modified, the following lines in simulation_analysis.C will also need to be modified:
Float_t lengthBox = 200, widthBox = 200;
The depthFix 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. depthFix 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):
//depth fix - water box centred at -2260 mm, half length = 100 mm. Double_t depthFix = 2360;
Sensitive detectors
In this simulation, the water volume is assigned as a sensitive detector in DetectorConstruction.cc :
G4SDManager* SDman = G4SDManager::GetSDMpointer(); G4String name="SD"; DetectorSD = new SensitiveDetector(name); SDman->AddNewDetector(DetectorSD); logicWater->SetSensitiveDetector(DetectorSD);
Another beamline component may be used by modifying the following line and setting its logical volume as a sensitive detector:
logicWater->SetSensitiveDetector(DetectorSD);
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.
SensitiveDetector.cc is derived from the G4VSensitiveDetector.cc base class. On a step-by-step basis, the energy deposited by the proton is recorded as a "hit" and added to a HitsCollection object. Other parameters may be retrieved at each step in the method ProcessHits .
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;
}
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 PhysicsList.cc , and then selected in protonPB.cc
runManager->SetUserInitialization(new PhysicsList);
Kinetic energy readings
The kinetic energy of the beam after passing through a particular component of the beamline is found in the class SteppingAction.cc . This class contains a map named componentMap 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):
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));
To find the kinetic energy distribution of the beam at a particular point, insert the component name from the map above into this line:
double cPos = (double)componentMap.find("ScatteringFoil2")->second;
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 kin.out :
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();
}
}
After any modifications to the simulation files, the code will need to be compiled. In the build directory, write:
[username@plus1 ProtonPB_build]$ make
After this proton.mac can be run:
[username@plus1 ProtonPB_build]$ ./protonPB proton.mac
Modifying Analysis Methods
The number of bins to store the data for each histogram can be changed, here the number of bins is 800:
TH1D* waterEnergy = new TH1D ("waterDep", "Energy deposition in water box", 800, 0., 65.);
The range of the axes for each plot can also be specified in order to zoom into the energy peak:
waterEnergy->GetXaxis()->SetRangeUser(0, 64);