|
|
| (322 intermediate revisions by 3 users not shown) |
| Line 1: |
Line 1: |
| == Simulation of the Clatterbridge beamline ==
| |
|
| |
|
| This simulation models the monoenergetic 62.5 MeV proton beam at Clatterbridge Cancer Centre, as it traverses the components of the beamline and is then deposited into a volume of water. The beamline components are contained within a geometry modelling the Clatterbridge treatment room. | | 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 protons are generated using the G4ParticleGun class, and the physics list used is QGSP_BIC_HP, standard for simulating clinical proton beams.
| | * [[/Geant4 | Geant4]] |
| | | ** 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 energy of the beam after travelling through the beamline components is measured by tracking the energy deposition of individual protons within the water volume, using an implementation of the <code> G4VSensitiveDetector.cc </code> class. This simulation produces a post-beamline energy of 60.08 MeV and a Bragg peak at a depth of 31 mm in water.
| | * [[/TOPAS | TOPAS]] |
| | |
| '''[insert image of visualisation]'''
| |
| | |
| == Running the simulation ==
| |
| | |
| '''Run macro proton.mac'''
| |
| | |
| This will run the simulation and produce the required output files.
| |
| | |
| <pre>
| |
| [username@plus1 ProtonPB_build]$ ./protonPB proton.mac | |
| </pre>
| |
| | |
| == Output files ==
| |
| | |
| The simulation code and <code> proton.mac </code> produce several output files:
| |
| | |
| === <code> kin.out </code> ===
| |
| | |
| This text file contains the output information from <code> SteppingAction.cc </code>, printed on a step-by-step basis for each proton (event). The first column contains the <code> z </code> position of the particle, relative to the position of the source at the inner room boundary. The second column contains the energy (MeV) of the proton at this <code> z </code> position.
| |
| | |
| <pre>
| |
| 300 62.3248
| |
| 300 62.0776
| |
| 300 62.2042
| |
| 300 62.4347
| |
| 300 62.2698
| |
| 300 62.2164
| |
| </pre>
| |
| | |
| === <code> FluxLongitudinal.txt </code> ===
| |
| | |
| This file is automatically written by <code> proton.mac </code> 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 cm<sup>2</sup>:
| |
| | |
| <pre>
| |
| # 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
| |
| </pre>
| |
| | |
| In <code> simulation_analysis.C </code> the data in <code> FluxLongitudinal.txt </code> is read in and written to another text file, <code> FluxLongitudinal_Mod.txt </code> to modify the standard formatting into a suitable format for further analysis:
| |
| | |
| <pre>
| |
| 0 0 0 25.97935669321147
| |
| 0 0 1 25.90015344712939
| |
| 0 0 2 25.84461190407653
| |
| </pre>
| |
| | |
| === <code> EnergyLongitudinalMesh.txt </code> ===
| |
| | |
| This text file is produced by <code> proton.mac </code> in the same way as <code> FluxLongitudinal.txt </code> above, and contains information about the proton energy deposited within the water volume, divided into bins parallel to the beam:
| |
| | |
| <pre>
| |
| # 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
| |
| </pre>
| |
| | |
| === <code> DoseLateralMesh.txt </code> ===
| |
| | |
| This text file is produced by <code> proton.mac </code> in the same way as the other scorers above, and contains information about the proton dose in bins perpendicular to the beam:
| |
| | |
| <pre>
| |
| # mesh name: waterMeshlateral
| |
| # primitive scorer name: doseDeposit
| |
| # iX, iY, iZ, value [Gy]
| |
| ...
| |
| 95,0,0,3.307349132489631e-09
| |
| 96,0,0,3.296548479981562e-09
| |
| 97,0,0,1.01396339893476e-09
| |
| </pre>
| |
| | |
| == Data Analysis ==
| |
| | |
| | |
| '''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 ===
| |
| | |
| This 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 from the water box sensitive detector.
| |
| | |
| === Proton stopping distance in water ===
| |
| | |
| | |
| === Proton flux along beamline ===
| |
| | |
| This plot is produced from <code> FluxLongitudinal_Mod.txt </code> and shows the number of protons per cm<sup>2</sup> 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).
| |
| | |
| === Kinetic energy of beam ===
| |
| | |
| This 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
| |
| | |
| == Changing parameters ==
| |
| | |
| === Initial beam parameters ===
| |
| | |
| Initial parameters of the proton beam can be modified in <code> proton.mac </code>
| |
| | |
| ==== Beam radius ====
| |
| | |
| <pre>
| |
| /gps/pos/radius 3 mm
| |
| </pre>
| |
| | |
| ==== Beam energy ====
| |
| | |
| This simulation models the proton beam source with 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 = -420 cm </code> relative to the centre of the inner room (the mother volume), which translates as the wall surface of the inner room.
| |
| | |
| <pre>
| |
| /gps/pos/type Plane
| |
| /gps/pos/shape Circle
| |
| /gps/pos/centre 0.0 0.0 -420 cm
| |
| </pre>
| |
| | |
| === 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 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 dose distribution of the protons. 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>.
| |
| | |
| === Kinetic energy readings ===
| |
| [could go under data analysis section] | |
| | |
| 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 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));
| |
| componentMap.insert( std::pair<std::string,double>("KaptonWindow",356.05));
| |
| componentMap.insert( std::pair<std::string,double>("AntiScatterCollimator1",1140.));
| |
| componentMap.insert( std::pair<std::string,double>("MonitorChambers",1150.02));
| |
| componentMap.insert( std::pair<std::string,double>("NozzleBefore",1692.));
| |
| componentMap.insert( std::pair<std::string,double>("NozzleAfter",1766.));
| |
| componentMap.insert( std::pair<std::string,double>("WaterSurface",1840.));
| |
| </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 ==
| |