|
|
| (20 intermediate revisions by 2 users not shown) |
| Line 1: |
Line 1: |
| 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 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 physics list used is QGSP_BIC_HP, the standard for simulating clinical proton beams.
| |
|
| |
|
| = Geometry =
| | 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. |
| === Treatment room ===
| | The simulation was built using: |
| The schematic below shows the Clatterbridge treatment room, illustrating the layout of the geometries defined in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/DetectorConstruction.cc <code>DetectorConstruction</code>]. 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.
| | * [[/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]. |
| | | * [[/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]]
| |
| | |
| === Dose monitor ===
| |
| [[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]]
| |
| | |
| = Generating primary protons =
| |
| === Using the <code>/gps/</code> commands ===
| |
| 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_sim/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 [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/gps.mac <code>gps.mac</code>] macro as follows:
| |
| <pre>
| |
| #-------------------------------------------------------------------------------
| |
| # Set particle gun settings
| |
| #-------------------------------------------------------------------------------
| |
| /gps/particle proton
| |
| /gps/number 1 # 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>
| |
| | |
| The macro is called in <code>proton.mac</code> as follows:
| |
| <pre>
| |
| #-------------------------------------------------------------------------------
| |
| # Primary generator settings
| |
| #-------------------------------------------------------------------------------
| |
| ...
| |
| ## Use macro to set properties of primaries
| |
| /control/execute gps.mac
| |
| </pre>
| |
| | |
| === Reading primaries from a phase space file ===
| |
| Primaries can be generated from an input phase space file using [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/FileReader.cc <code>FileReader</code>] in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/PrimaryGeneratorAction.cc <code>PrimaryGeneratorAction</code>]. The input must be of the same format as the output files generated by the simulation and its first line is assumed to be a header. The path to the input file should be given in <code>proton.mac</code> and the position at which the particles should be generated:
| |
| <pre>
| |
| #-------------------------------------------------------------------------------
| |
| # Primary generator settings
| |
| #-------------------------------------------------------------------------------
| |
| ...
| |
| ## Reading from phase space file
| |
| /primarygenerator/input path/to/file
| |
| | |
| # Set position in z particles are generated at
| |
| /primarygenerator/generateAt 0
| |
| #/primarygenerator/generateAt 81
| |
| #/primarygenerator/generateAt 357
| |
| #/primarygenerator/generateAt 1700
| |
| </pre>
| |
| Specifying where particles should be generated in useful when running staged simulations as it allows particles to be generated at the same position they were initially recorded at.
| |
| | |
| = Gathering data =
| |
| == Tracking ==
| |
| === Accumulating hits over a run ===
| |
| Several quantities can be monitored at different stages in the beamline using the <code>SensitiveDetector</code> defined in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/PhaseSpaceSD.cc <code>PhaseSpaceSD.cc</code>] (<code>PhaseSpaceSD</code> extends [http://www.apc.univ-paris7.fr/~franco/g4doxy/html/classG4VSensitiveDetector.html <code>G4VSensitiveDetector</code>]). When a particle hits a boundary on which a <code>PhaseSpaceSD</code> is defined, <code>PhaseSpaceSD::ProcessHits</code> is called and generates a <code>PhaseSpaceHit</code> (extension of [http://www.apc.univ-paris7.fr/~franco/g4doxy/html/classG4VHit.html <code>G4VHit</code>]) that is then stored in the hits collection of the sensitive detector. A unique ID belonging to the volume the hit was registered on is passed to the hit along with any required quantity. Hits are accumulated in the sensitive detector's hits collection over an event (an event corresponds to the full simulation of one primary).
| |
| | |
| The [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/RunAccumulator.cc <code>RunAccumulator</code>] (extension of [http://www.apc.univ-paris7.fr/~franco/g4doxy/html/classG4Run.html <code>G4Run</code>]) accumulates hits over a run (the collection of all events). It is constructed by <code>RunAction::GenerateRun</code> at the start of a run. At the end of an event, <code>RunAccumulator::RecordEvent</code> stores all hits that were registered during the event in the run hits vector and increments the event counter. Once that counter reaches the buffer, the hits are written to their corresponding output file. The buffer size can be set using the following line in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/proton.mac <code>proton.mac</code>]:
| |
| <pre>
| |
| #-------------------------------------------------------------------------------
| |
| # Detector settings
| |
| #-------------------------------------------------------------------------------
| |
| ...
| |
| # Set buffer size
| |
| /user/run/buffer 1000
| |
| </pre>
| |
| | |
| === Recording at regular intervals using a parallel world ===
| |
| To characterise the beam, the evolution of several properties (energy spectrum, emittance, spatial profile, etc.) must be monitored at regular intervals, or at least at several positions, along the beamline. To achieve this, many volumes would need to be defined in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/DetectorConstruction.cc <code>DetectorConstruction</code>] to be used as sensitive detectors. However, volumes overlapping or sharing boundaries can lead to errors in the simulation and must therefore be avoided. One approach in achieving this would be to carefully construct each detector volume such that it fits in, or around, a given region. An example illustrating why this may be challenging is given by the following; suppose detection was required inside and outside of the first aluminium box on a plane transverse to the direction of the beam. Separate detector volumes would be required on the inside and outside of the box just to detect particles incident on the same plane. This does not yet address the issue of how to handle the boundaries between the box and the detectors as particles would be likely to pass through the gaps and hence go undetected. Even for this simple example, the geometry of the detector volumes would be rather complex and hence not easily reusable.
| |
| | |
| To avoid these problems, the main geometry (called the ''mass geometry'') can be overlaid with another parallel geometry that does not interact physically with particles in the mass geometry, as described [http://geant4.web.cern.ch/geant4/workAreaUserDocKA/Backup/Docbook_UsersGuides_beta/ForApplicationDeveloper/html/ch04s07.html here]. Different volumes can then be defined at arbitrary positions and have sensitive detectors assigned to them in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/ParallelWorldConstruction.cc <code>ParallelWorldConstruction</code>]. To track particles at regular intervals along the beamline, an array of thin boxes is defined by passing the parametrisation [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/ParallelWorldParam.cc <code>ParallelWorldParam</code>] (extension of [http://www.apc.univ-paris7.fr/~franco/g4doxy/html/classG4VPVParameterisation.html <code>G4VPVParametrisation</code>]) to the parametrised volume constructor [http://www.apc.univ-paris7.fr/~franco/g4doxy/html/classG4PVParameterised.html <code>G4PVParametrised</code>]. This creates boxes of a given size at given intervals spanning a given distance from the source. The size of the box can be adjusted in <code>ParallelWorldConstruction</code>, while the spacing and total distance spanned can be set in the <code>proton.mac</code> macro using the following settings:
| |
| <pre>
| |
| #-------------------------------------------------------------------------------
| |
| # Detector settings
| |
| #-------------------------------------------------------------------------------
| |
| ...
| |
| ## If using 'all'
| |
| /parallel/detector/spacing 25 mm # default is 25 mm
| |
| /parallel/detector/distance 1900 mm # default is 1900 mm
| |
| | |
| # Pass spacing and distance to run action and set buffer size
| |
| /user/run/detector/spacing 25
| |
| /user/run/detector/distance 1900
| |
| </pre>
| |
| | |
| === Recording in a single position (for staged simulations) ===
| |
| Alternatively, particles can be recorded in a single position by choosing one of the components and uncommenting accordingly. The available components are
| |
| * scatterfoil1 at 81 mm
| |
| *tube1 at 357 mm
| |
| *nozzle at 1700 mm
| |
| *outside at 1768 m
| |
| Components can be added to the components map in <code>ParallelWorldConstruction</code>. Only recording in one position can be useful when the beam has been characterised and the simulation is know to be running as intended. Writing to output less often improves the simulation's performance significantly. This can be used to speed up the simulation further by performing simulations in a staged manner. Sections upstream, where no changes to the geometry are made (e.g. from the source to after the first collimator where the particle count reduces by about 90%), then only need to be simulated once and the output from that simulation can be used as the input to the rest of the simulation.
| |
| | |
| For example, if the simulation should only run up to the first scatter foil (situated after the first collimator), the commands would look as follows:
| |
| <pre>
| |
| #-------------------------------------------------------------------------------
| |
| # Detector settings
| |
| #-------------------------------------------------------------------------------
| |
| ...
| |
| ## Choose where particles are recorded
| |
| #/parallel/detector all
| |
| /parallel/detector scatterfoil1
| |
| #/parallel/detector tube1
| |
| #/parallel/detector nozzle
| |
| #/parallel/detector outside
| |
| ...
| |
| ## If using a component, set dump mode and detector position
| |
| /user/run/dump/single true
| |
| /user/run/detector/position 81
| |
| #/user/run/detector/position 357
| |
| #/user/run/detector/position 1700
| |
| #/user/run/detector/position 1768
| |
| | |
| ...
| |
| | |
| #-------------------------------------------------------------------------------
| |
| # Stepping action settings
| |
| #-------------------------------------------------------------------------------
| |
| ...
| |
| ## Kill particles 1 mm after being recorded
| |
| /steppingAction/kill 82
| |
| #/steppingAction/kill 358
| |
| #/steppingAction/kill 1701
| |
| #/steppingAction/kill 1769
| |
| </pre>
| |
| | |
| The command <code>/user/run/dump/single</code> sets a flag used in <code>RunAccumulator::InitialiseOutputFiles</code>.
| |
| | |
| == Scoring ==
| |
| Scorers are defined in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/score_init.mac <code>score_init.mac</code>] and the quantities they record are dumped to files as defined in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/score_dump.mac <code>score_dump.mac</code>].
| |
| | |
| === 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. It is defined as follows:
| |
| | |
| <pre>
| |
| /score/create/boxMesh detEnergyLon
| |
| /score/mesh/boxSize 20. 20. 20. mm
| |
| /score/mesh/nBin 1 1 200
| |
| /score/mesh/translate/xyz 0. 0. -2340 mm
| |
| /score/quantity/energyDeposit energyDeposit
| |
| /score/close
| |
| </pre>
| |
| | |
| On the third line in the snippet above you can see that there is only a single bin in both the x and y directions but 200 in the z-direction.
| |
| | |
| === Lateral scoring mesh ===
| |
| Similarly to the longitudinal meshes, a mesh is said to be ''lateral'' when it records data transverse to the direction of travel of the protons. For a mesh to be lateral it must have bins defined perpendicularly to its z-axis. The scorer defined on the mesh <code>braggEnergyLat</code> is an example of such a mesh. It is defined as follows:
| |
| | |
| <pre>
| |
| /score/create/boxMesh braggEnergyLat
| |
| /score/mesh/boxSize 20. 20. 0.5 mm
| |
| /score/mesh/nBin 200 1 1
| |
| /score/mesh/translate/xyz 0. 0. -2329.2 mm
| |
| /score/quantity/energyDeposit energyDeposit
| |
| /score/close
| |
| </pre>
| |
| | |
| Again, you can see that there is only a single bin in the y and z directions but 200 in the x-direction.
| |
| | |
| = Running simulations =
| |
| | |
| === Build 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/users/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 -j4
| |
| </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@pc1XX build]$ mkdir data
| |
| [username@pc1XX build]$ ./clatterbridgeSim proton.mac
| |
| </pre>
| |
| | |
| === Running simulations in parallel ===
| |
| 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 parallel simulations, some parameters need to be set in the macro [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/proton.mac <code>proton.mac</code>].
| |
| *When running simulations in parallel, each simulation is executed two directories further down from where the parallel run is initialised. Hence, the correct path to the <code>gps.mac</code> macro should be set in <code>proton.mac</code>:
| |
| <pre>
| |
| #-------------------------------------------------------------------------------
| |
| # Primary generator settings
| |
| #-------------------------------------------------------------------------------
| |
| ...
| |
| # Parallel version
| |
| /control/execute ../../gps.mac
| |
| </pre>
| |
| | |
| *Separate simulations should have different random number generator seeds so that they are independent but in order for simulations to be reproducible, the same sequence of seeds should be assigned for each run. To achieve this, the seeds in <code>proton.mac</code> are set by the job submission script [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/cb_parallel.sh <code>cb_parallel.sh</code>] using the index of the current simulation. For the <code>proton.mac</code> macro to use the seeds set by the submission script, the following line should be set:
| |
| <pre>
| |
| #-------------------------------------------------------------------------------
| |
| # Set seeds for randon number generators
| |
| #-------------------------------------------------------------------------------
| |
| ...
| |
| # Parallel version
| |
| /random/setSeeds ${seed1} ${seed2}
| |
| </pre>
| |
| | |
| * If generating primaries from a phase space file that was split using [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/split_input.sh <code>split_input.sh</code>] as described below, the number of events is set by <code>cb_parallel.sh</code> in <code>proton.mac</code> and the following line should be uncommented: | |
| <pre>
| |
| #-------------------------------------------------------------------------------
| |
| # Run simulation
| |
| #-------------------------------------------------------------------------------
| |
| ...
| |
| # Parallel version with split input file
| |
| /run/beamOn ${nevents}
| |
| </pre>
| |
| | |
| As mentioned above, simulations can be run in parallel while reading primaries from an input file. This file should be split and distributed evenly across the simulations using [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/split_input.sh <code>split_input.sh</code>] in order to make use of the parallelism. The script calculates the number of lines each simulation should receive given a number of simulations and then writes chunks of the initial file to separate files which will then be used as input files for the parallel simulation. It should be called in the simulation's base directory as follows:
| |
| <pre>
| |
| ./split_input.sh /path/to/file nsimulations
| |
| </pre>
| |
| | |
| The script [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/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_sim/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_sim/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_sim/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_sim/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]$ ./run.sh
| |
| </pre>
| |
| | |
| === Monitoring job status ===
| |
| | |
| 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_sim/source/combine_tracking.py <code>combine_tracking.py</code>] and [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/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_sim/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_sim/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_sim/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_sim/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_sim/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_sim/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]]
| |
| | |
| == Changing parameters ==
| |
| | |
| === Beam parameters ===
| |
| These can be modified in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/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 ===
| |
| Filters can be applied to scoring meshes so as to only record a given type of particle using the command <code>/score/filter/particle protonFilter proton</code> in the definition of the scoring mesh. This can be useful as, for example, one may want to monitor the flux of different particles along the beamline. This would require a large mesh positioned over the whole beamline as follows:
| |
| | |
| <pre>
| |
| /score/create/boxMesh beamlineLongitudinal
| |
| /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>
| |
| | |
| This mesh records two separate quantities, the proton flux and the neutron flux, which can then be written to separate output files.
| |
| | |
| === Beamline components ===
| |
| Components of the beamline can be added, removed or modified in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/DetectorConstruction.cc <code>DetectorConstruction.cc</code>]. For example, the position of the second dose monitor can be changed by changing the z-value of the <code>G4ThreeVector</code> on the line below. As the dose monitor is placed within the dosimetry box <code>lAlBox2</code> it's position is given relative to the centre of that box.
| |
| | |
| <pre>
| |
| pDoseMonitor2 = new G4PVPlacement(0, G4ThreeVector(0, 0, -169.5*CLHEP::mm), lDoseMonitor1, "DoseMonitor2", lAlBox2, false, 0, checkOverlaps);
| |
| </pre>
| |
| | |
| '''After any changes made to the source or header files (ending in <code>.cc</code> or <code>.hh</code>), the code will need to be recompiled. In the build directory, write:'''
| |
| | |
| <pre>
| |
| [username@pc1XX build]$ make
| |
| </pre>
| |
| | |
| === Physics list ===
| |
| A user-defined physics list can be specified in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/PhysicsList.cc <code>PhysicsList.cc</code>] and then selected in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/clatterbridgeSim.cc <code>clatterbridgeSim.cc</code>] using the line
| |
| | |
| <pre>
| |
| runManager->SetUserInitialization(new PhysicsList);
| |
| </pre>
| |
| | |
| === Tracking positions ===
| |
| As described in the section "[http://www.hep.ucl.ac.uk/pbt/wiki/Clatterbridge#Tracking Tracking] ", protons are tracked in [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/SteppingAction.cc <code>SteppingAction.cc</code>]. This class contains a vector <code>trackingPosVec</code> which contains the positions of interest where protons are to be tracked. Positions can be changed by changing <code>trackingPosFiller</code> in the header file [http://www.hep.ucl.ac.uk/pbt/wikiData/code/Clatterbridge_sim/source/src/SteppingAction.hh <code>SteppingAction.hh</code>] on the following line:
| |
| | |
| <pre>
| |
| trackingPosFiller = { 0, 45, 79, 81, 100, 150, ..., 1759, 1800, 1839, 1881 };
| |
| </pre>
| |
| | |
| '''Don't forget to recompile after any changes made.'''
| |
| | |
| = Files =
| |
| | |
| [[Clatterbridge files|List of files for Clatterbridge simulation]].
| |