The Geant4 STCyclotron Advanced Example

Introduction

The STCyclotron example simulates the solid target of the South Australian Health and Medical Research Institute (SAHMRI), Adelaide, South Australia. The simulation aims at modelling the production of radio-isotopes and undesired by-products. An example of the use of the Geant4 STCyclotron application in research is available in F. Poignant, S. Penfold, J. Asp, P. Takhar, P. Jackson, GEANT4 simulation of cyclotron radioisotope production in a solid target, Physica Medica, Volume 32, Issue 5, 2016, Pages 728-734, ISSN 1120-1797 (https://doi.org/10.1016/j.ejmp.2016.04.006)

Authors: F. Poignant (IN2P3, France), S. Penfold (University of Adelaide, Australia), P. Takhar and J. Asp (SAHMRI, Adelaide). Contributor: S. Guatelli (University of Wollongong, Australia).

For more information on the simulated cyclotron, please visit: https://www.comecer.com/new-research-study-simulation-solid-target-performance

Pre-requisites

In order to execute the Advanced Example, Geant4 should be installed with the following options: GEANT4_BUILD_MULTITHREADED (preferable)

GEANT4_INSTALL_DATA

GEANT4_USE_OPENGL_X11

GEANT4_USE_QT

GEANT4_USE_SYSTEM_CLHEP

GEANT4_USE_EXPAT

For more information on how to install GEANT4, please visit: https://geant4.web.cern.ch/

ROOT may be used as analysis tool. To download and install the latest version, please visit: https://root.cern.ch/

Set-up of the proton-low energy database

Users should make sure that the database used for inelastic collisions of primary particles is set up correctly. This database is not installed by default. TENDL and ENDF data libraries can be downloaded at the following links:

http://geant4.web.cern.ch/geant4/support/download.shtml (TENDL)
ftp://gdo-nuclear.ucllnl.org/pub/G4LEND/ (G4 Low Energy Nuclear Data)

In your simulation environment, add the following:

export G4PROTONHPDATA=/PATHTO_TENDL____ *OR* ____ENDF_DATABASE/Proton
export G4NEUTRONHPDATA=/PATHTO_GEANT4_INSTALLATION_FOLDER/share/Geant4-vXX.XX/data/G4NDL_VERSION
export G4PHP_DO_NOT_ADJUST_FINAL_STATE=1export G4PHP_MULTIPLICITY_METHOD=Poisson

Compile the example

The compilation of this example is analogous to the other examples.

Create a directory 'STCyclotron-build'. Build the example by typing

cmake ../STCyclotron 
make

Input parameters

The simulation parameters may be changed either in the GUI session or in the Macro files in 'STCyclotron/Macro'.

Macro/init_parameters.mac

To study the optimization of the isotope production, a list of parameters can be changed in the Macro/init_parameters.mac (or using the Geant4 User Interface):

Beam parameters

Different parameters can be changed for the beam: type of particle, energy, energy distribution, shape of the beam, etc ... To design your own beam, please refer to the Geant4 User's Guide for Application Developers. In addition you may change:

• setBeamCurrent/beamCurrent 30E-6 -- current (ampere).

• setTimeOfIrradiation/time 6. -- time of irradiation (hour(s))

As Geant4 doesn't model any time scale explicitly, the number of primary particles per event represents 10E-11 seconds of a real experiment. It is calculated as follows:

NumberOfParticle = beamCurrent*timePerEvent/chargeParticle

where beamCurrent is a parameter that can be changed (ampere); timePerEvent is set to 10E-11 second; chargeParticle is the charge of the particle (Coulomb, e.g. 1.9E-19 C for proton). For 30 μA, one event will then correspond to 1 875 protons. This value was chosen so the number of protons per event is not too high. If you work on high current (over 100 microA) you might need to change the set up. This time is defined in the PrimaryGeneratorAction class. Note that, for one event, all the protons are initialized at the same primary coordinates.

The user can choose the time of irradiation by changing the parameter. Note that a simulation for a real time of irradiation would be too long. The final number of produced radioisotopes is analytically calculated using the number of radio-isotopes obtained during the simulation, which represents a fraction of the total time of the experiment.

Target parameters:

• /changeTarget/diameter 6. (mm) -- the target thickness. Must be smaller than 38.32 mm, which is the length of the tube containing the target.

• /changeTarget/thickness 0.6 (mm) -- the target diameter. Must be smaller than 15 mm, which is the diameter of the tube containing the target.

• the target material: there are two ways to change the material of the target. If the material is natural, the user can choose to use the NIST database. If the material is enriched, the user can set up their own material. The material created is made of a number of elements. The number of elements can be set up by the user. Then, one element can be created in two different ways : it can be a natural element, using the NIST database, or it can be made of a number of isotopes that the user can set up. Please, note that the order to declare parameters is important and has to be the following :
  1. Material settings
  2. Element i settings
  3. Isotopes settings for the element i
  4. Element i+1 settings
  5. Isotopes settings for the element i+1, .... and so on.

For example, to create an target of nickel enriched to 60% of Ni-64 : one will create a new material, made of two elements : natural nickel and pure Ni-64. The natural nickel will be an element created using the NIST database. The pure Ni-64 will be an element made of one isotope : Ni-64. In case you want to create a pure Ni-64 target, the materials made of one element, made itself of one isotope (Ni-64). Few examples are provided in the folder "Macro/Material/Target" that can be executed in the init_parameters.mac.

Foil parameters:

• /changeFoil/thickness 0.32 (mm) -- the foil thickness.

• the foil material. With this parameter, it is possible to change the foil material in order to study some specific aspects of the reaction. Refers to the change of target material for more detail.

Histograms:

The histograms parameters can be changed in order to fit to the expected range for a given data. For example, for a proton beam with an expected energy of 16 MeV when reaching the target, the beam energy profile range can be set up between 15. to 17. MeV.

Macro/Vis/vis.mac

To visualize all tracks you may want to increase the visualization display limit by adding the following GUI command in vis.mac: /vis/ogl/set/displayListLimit 1000000

Macro/GUI/gui.mac

This file sets up the tool bars and buttons that enables to modify the parameters using the Geant4 User Interface.

Run the simulation

To launch the simulation, simply type

./STCyclotron

If the Graphical User Interface is activated and theGeant4 environment is correctly set, a Graphical User Interface should open.

On the top, the tool bar enables the user to execute the different commands. There are few menus: one for the beam parameters, one for the target material, one for the geometry of the target, one for the foil material, and one for the foil geometry. Through this tool bar, the main parameters of the simulation set-up may be changed. Remember! Don't set up the diameter of the particle beam larger than the diameter of the target.

To generate the proton beam, type in the session part

/run/beamOn 1

One event, which contains a given number of particles (protons), will be executed. The simulation is displayed in the viewer window, as illustrated below.

Green particles are photons, blue ones are proton, red ones are electrons and yellow ones are positrons. Displaying the interactions is very time consuming, so when running more than 10 events it is strongly advised to disable the viewer.

Results

In order to have enough statistics, it is needed to run several events (typically at least of the order of 10 000). Once the simulation is finished, the results are saved in the build directory. To exit Geant4, type “exit” in the session part or use the tool bar.

In this example ROOT is used to analyse the results. The user may choose another analysis tool. If using ROOT, execute the file Plot.C by typing

root Plot.C

It will create many PDF files. To exit ROOT, type

.q

Different types of output are available. Examples below are given for a 18 MeV proton beam, a nickel target of thickness equal to 0.6 mm and diameter of 6 mm, and a foil of thickness equal to 0.32 mm. The time of irradiation was set equal to 6h and the current at 30 microA.

ROOT file (STCyclotron.root)

This ROOT file is generated as output of the simulation and contains a list of histograms representing the following data :

• 1D histograms

1. The energy distribution of primary particles (e.g. protons) when reaching the target (MeV).
2. The energy distribution of primary particles (e.g. protons) when reaching the foil (MeV).
3. The energy distribution of primary particles (e.g. protons) going out of the target (MeV).
4. The energy distribution of primary particles (e.g. protons) going out of the foil (MeV).
5. The depth of isotope production in the target (number of particles as a function of the foil thickness in mm).
6. Energy spectrum of particles produced in the target following inelastic collision of primary particles (e.g. protons) with the target material (MeV). In order: 5 = positrons; 6 = electrons; 7 = gammas; 8 = neutrons.
7. Energy spectrum of particles produced in the target following decay of isotopes produced in the target (MeV). In order: 9 = positrons; 10 = electrons; 11 = gammas; 12 = neutrons; 13 = nu; 14 = anti_nu (electron (anti)neutrinos).

• 2D histograms :

1. The beam intensity profile before hitting the target (mm x mm).
2. The beam intensity profile before hitting the foil (mm x mm).
3. The radioisotopes produced according to their Z and A number.
4. The energy of the primary particles (e.g. protons) according to depth in the target (mm x MeV).
5. The beam intensity going out from the target (mm x mm).
6. The beam intensity going out from the foil (mm x mm).

/!\ the histograms are not normalized /!\. The file 'Plot.C' renormalize the histograms and plot them into PDFs as explained below.

Nevertheless, it is possible to visualize the outputs by typing

root (which launchs the software ROOT)
new TBrowser (which opens a window to visualize the ROOT file)

Text files

• Output_General.txt

This file summarizes the parameters used during the simulation:

1. Beam parameters: primary particles (by default protons), energy of the primary particles (MeV), current of the beam (Ampere), irradiation time (hour(s)), and current factor. This last factor is a rescaling factor: in the simulation, the number of particles sent is calculated for a current obtained before the foil, while the current in the actual cyclotron the current is measured after the foil. This parameter therefore rescales the number of particles to match the current arriving at the target.
2. Simulation parameters: equivalent time per event (by default set at 10^-11 second), number of events run during the simulation, number of primaries per event (calculated according to the time per event, the beam current and the charge of the primary particle), total number of particles sent during the simulation.
3. Geometry parameters: target thickness, diameter and foil thickness.
4. It also provides the heating of the target and the foil (W/mm3).

Example of output:

//-----------------------------------//
// Parameters of the simulation: //
//-----------------------------------//
Beam parameters:
proton - Name of beam primary particles.
18 - Energy of beam primary particles (MeV).
3e-05 - Beam current (Ampere).
6 - Irradiation time in hour(s).
0.731972 - Current factor.
//-----------------------------------//
Simulation parameters:
1e-11 - Equivalent time per event (s).
1000 - Number of events
1875 - Primaries per event
1875000 - Total number of particles sent.
//-----------------------------------//
Geometry parameters:
0.6 - target thickness (mm).
6 - target diameter (mm).
0.32 - foil thickness (mm).
//-------------------------------------------------//
// Heating, total activity and process data //
//-------------------------------------------------//
Total heating in the target : 26.0042 W/mm3
The total heating during the irradiation is 4.72859e+19J/mm3
Total heating in the foil : 0.000713828 W/mm3

---

• Output_ParentIsotopes.txt

This file provides a list of radioisotopes produced during the irradiation of the target. For each isotope, it contains:

1. Name of the isotope.
2. Number of isotopes created during the simulation. Can be used to evaluate the accuracy of your predictions.
3. Decay constant (s-1).
4. Half life time (hour(s)).
5. Process that induced its creation.
6. Number of isotopes produced per second of irradiation.
7. Number of isotopes produced at the end of the beam.
8. Activity induced by the isotope at the end of the beam (mCi).

Example of output:

//-----------------------------------//
// Data for parent isotopes //
//-----------------------------------//

Co55 - name of parent isotope.
60 - number of isotopes created during the simulation.
1.09835e-05 - decay constant in s-1.
17.53 - half life time in hour(s).
protonInelastic - creation process.
4.39183e+09 - isotope per sec.
8.44502e+13 - yield EOB.
25.0441 - activity (mCi) at the EOB.
------------------------
Co57 - name of parent isotope.
223 - number of isotopes created during the simulation.
2.95228e-08 - decay constant in s-1.
6521.76 - half life time in hour(s).
protonInelastic - creation process.
1.6323e+10 - isotope per sec.
3.52464e+14 - yield EOB.
0.280955 - activity (mCi) at the EOB.
------------------------
Cu64 - name of parent isotope.
26 - number of isotopes created during the simulation.
1.51595e-05 - decay constant in s-1.
12.701 - half life time in hour(s).
protonInelastic - creation process.
1.90313e+09 - isotope per sec.
3.50555e+13 - yield EOB.
14.3484 - activity (mCi) at the EOB.
------------------------

Etc ...

To have sufficient statistics on the radio-isotope production, the number of isotopes created during the simulation should be larger than about one thousand.

---

• Output_DaughterIsotopes.txt

This file provides a list of unstable daughter radioisotopes produced due to the decay on unstable primary (parent) radiosotopes (if applicable). As for the file Output_ParentIsotopes.txt, it contains:

1. Name of the daughter isotope.
2. Name of the parent isotope.
3. Decay constant of the parent isotope (s-1).
4. Decay constant of the daughter isotope (s-1).
5. Half life time of the parent isotope (hour(s)).
6. Half life time of the daughter isotope (hour(s)).
7. Number of daughter isotopes produced per second of irradiation.
8. Number of daughter isotopes produced at the end of the beam.
9. Activity induced by the daughter isotope at the end of the beam (mCi).

---

• Output_StableIsotopes.txt

For information, this file provides a list of stable isotopes (name and number of isotopes produced during the simulation) that are produced in the target due to the decay of radioisotopes.

---

• Output_Particles.txt

For information, this file provides a list of other particles such as electrons, etc., (name and number of isotopes produced
during the simulation) that are produced in the target.

These text files are read by the Plot.C file to extract the parameters and results.

PDF files

After running the 'Plot.C' file, PDF files are created in a folder 'Results'. This code reads the different outputs from the simulation (.root file and .txt files), normalize the results and plot them in PDFs in various folders:

• Results/BeamData folder

1. BeamEnergyInFoil.pdf and BeamEnergyInTarget.pdf: beam energy distribution before entering the foil/target using histograms 1D0 and 1D1, normalized to the number of primary protons and the bin width.
2. BeamEnergyOutFoil.pdf and BeamEnergyOutTarget.pdf: beam energy distribution when exiting the foil/target, using histograms 1D2 and 1D3, normalized to the number of primary protons and the bin width.
3. BeamIntensityInFoil.pdf and BeamIntensityInTarget.pdf: beam intensity before entering the foil/target using histograms 2D0 and 2D1, normalized per primary particle and to the bins widths.
4. BeamIntensityOutTarget.pdf: beam intensity when exiting the target using histogram 2D4, normalized per primary particle and to the bins widths.
5. EnergyDepth.pdf: energy of protons as a function of the depth in the target.

Example of outputs are shown here:

• Results/IsotopesProduction

1. ActivityOfXX.pdf and YieldOfXX.pdf: Shows the production of the isotope XX (number of nuclei or activity) as a function of the time, starting from the beginning of irradiation and up to 30 hours. Note that if the time of irradiation is longer than 30 hours, you must change the maximum time to display the activity or yield by opening the file 'Plot.C' and changing tMax.
2. ActivitySaturationOfXX.pdf and YiedSaturationOfXX.pdf: Shows the saturation reached for the production of the isotope XX (number of nuclei or activity) as a function of the time, if the time of irradiation is set 'infinite'.
3. Activity.pdf/Activity.jpg and Yield.pdf/Yield.jpg: Shows the activity (or yield) of all the isotopes produced during the irradiation as a function of the time up to 30 hours on the same graph.
4. TotalActivity.pdf: Shows the sum of the activities induced by all the radioisotope up to 30 hours.
5. RadioisotopeProduction.pdf/RadioisotopeProduction.jpg: Shows the number of isotopes produced per primary particles, as a function of Z and A.
6. DepthCreation.pdf: Shows the depth at which radioisotopes were created.

An example of output is here below:

• ParticlesEnergySpectra folder

1. Subfolder: beam. Energy spectra (normalized per primary particles and bin width) of particles created following the inelastic interaction of the beam with the target (1D 5->8).
2. Subfolder: decay. Energy spectra (normalized per primary particles and bin width) of particles created following the decay of radioisotopes created in the target (1D 9->14).

Analysis of the results

The accuracy and consistency of the simulation results can be verified by checking the nuclear database you are using. Go to the following website : http://www.oecd-nea.org/janis/book/

In the web access part, click on the “protons” to access the database of protons. Click on the atomic species of the target. For example, for the production of Cu-64, Ni-64 is used, so you will click on 28-Ni. The list of isotopes of Nickel is available. Click on Ni-64 (Z=28) and select the nuclear reaction of interest. The cross sections will be displayed for different nuclear databases and experiments.

The year of the nuclear databases displayed might be different from the one you are using. To display the ones used in the simulation, click on the following link :
http://www.oecd-nea.org/janisweb/search/endf

Enter the parameters of the irradiation. Here is an example:

Then click on “Open selected”.

In “Energy-angle distribution / Ni64” then MT=5, look for the product of interest. For the case of Cu-64 production, tick P and T for production cross section. It will display the plot of the cross section according to the specific incident energy and the computed data.

The computed values can be used therefore to be compared to experimental ones using the EXFOR website. Go on this website : https://www-nds.iaea.org/exfor/exfor.htm

Enter the parameters of the irradiation:

Then click on submit. There will be different experimental data available. Tick the ones of interest. Tick “Quick plot” and then click on retrieve. A plot will be displayed with different experimental values. On the right, tick “use my data”. Do it and add the data from the JANIS Book website. Make sure the units are coherent. Rename your data. Tick “authors” and “legend”. Then click on repaint:

Plots can be saved by clicking on PostScript and selecting a PDF format. These tools can be used to check the accuracy of the simulation. Yields are calculated using Computed Data which can differ from experimental data. This may induce inaccuracy. Always remember to verify how good the cross section are compared to the experimental data if possible.

-- SusannaGuatelli - 2019-11-25

• Figure1:
  

• IMAGE3.png:
  

• IMAGE5.png:
  

• IMAGE6.jpg:
  

• IMAGE7.jpeg:
  

• IMAGE10.png:
  

• IMAGE11.png:
  

• IMAGE12.png:
  

• IMAGE13.png:
  

• IMAGE14.png:
  

• IMAGE15.png: