Applied Radiation and Isotopes Validating production of PET radionuclides in solid and liquid targets: Comparing Geant4 predictions with FLUKA and measurements

(cid:129) Geant4 useful tool for simulation of PET isotope production. (cid:129) Geant4 and FLUKA results consistent. (cid:129) Geant4 using TENDL cross sections with QGSP-AllHP model best compromise. (cid:129) Model QGSP-BERT-HP and QGSP-BIC-HP do not produce all isotopes. The Monte Carlo toolkit Geant4 is used to simulate the production of a number of positron emitting radio- nuclides: 13 N, 18 F, 44 Sc, 52 Mn, 55 Co 61 Cu, 68 Ga, 86 Y, 89 Zr and 94 Tc, which have been produced using a 13 MeV medical cyclotron. The results are compared to previous simulations with the Monte Carlo code FLUKA and experimental measurements. The comparison shows variable degrees of agreement for di ﬀ erent isotopes. The mean absolute deviation of Monte Carlo results from experiments was ± 1.4 1.6 for FLUKA and ± 0.7 0.5 for Geant4 using TENDL cross sections with QGSP-BIC-AllHP physics. Both agree well within the large error, which is due to the uncertainties present in both experimentally determined and theoretical reaction cross sections. Overall, Geant4 has been con ﬁ rmed as a tool to simulate radionuclide production at low proton energy.


Introduction
Radioisotopes play a crucial role in the diagnosis and treatment of cancer. Numerous isotope-producing nuclear reactors are due to end their operation within a few years. As a result, proton-induced reactions have attracted significant interest from the scientific community after cyclotrons proved to be a feasible alternative to reactor produced radioisotopes Schaffer et al., 2015). Currently cyclotrons can be used to produce radioisotopes for imaging techniques such as positron emission tomography (PET) and single photon emission computed tomography (SPECT). The irradiated target can be in solid, liquid or gaseous form and may be required to satisfy strict design constraints. For example, a target may have material composition restrictions to achieve a desired specific activity, proton energy constraints to avoid unwanted isotope production, or may need to survive several hours of proton irradiation without any thermal issues. As a result, cyclotron targets and materials can be very expensive. Monte Carlo (MC) simulations can be used to assess the expected yield and to optimize target design and materials to maximize yield of the isotope of interest without increasing the production of contaminants (Infantino et al., 2011;Remetti et al., 2011;Sadeghi et al., 2013;Fassbender et al., 2007). The success in using MC for yield assessment depends strongly on the cross section data used for the simulation. Despite a large number of experiments carried out with proton activation, the data available are often inconsistent and at times data from different experiments conflict with each other.
In this work, the MC package Geant4 has been used to simulate the yields of the following PET isotopes: 13 N, 18 F, 44 Sc, 52 Mn, 55 Co 61 Cu, 68 Ga, 86 Y, 89 Zr and 94 Tc. The results have been compared to our pre-

Experiments
The experimental details have been described and, where appropriate, referenced by Infantino et al. (2016). Some details are repeated here for the convenience of the reader.
The TR13 cyclotron is located at TRIUMF, Vancouver, Canada and used for routine production of medical isotopes. It is self shielded and accelerates negative hydrogen ions to 13 MeV energy with currents of routinely up to 25 μA. Extraction occurs with the use of a carbon foil which strips off the two electrons thus reversing the charge and bending trajectory of the ion in the magnetic field. The cyclotron has two extraction ports with a target selector, which can move the target into the proton beam. Further details of the cyclotron are provided by Laxdal et al. (1994); Buckley et al. (2000).
The selector has four positions, allowing eight different targets to be installed at a time. Two target assemblies were simulated in Geant4. Figs. 1 and 2 illustrate the liquid and solid target assemblies respectively with each component labeled numerically. The proton beam enters the assembly though the baffle (1) and collimator rings (2). The beam is then collimated further with a four quadrant conical collimator (3) contained within an insulator flange (4). Each quadrant of the collimator is capable of measuring beam current separately and the four readings can be used to deduce the position of the proton beam. The beam then enters the target assembly through a 25 µm thick aluminium foil (5), which separates the cyclotron vacuum from the target assembly. Due to the power deposition in the foil, helium cooling (6) is applied to the foil in the helium window (7).
The liquid target (9) is a closed volume of 0.9 ml capacity, with 8 mm depth and 12 mm diameter. The liquid target is separated from the helium cooling (6) by a HAVAR foil (8). HAVAR is a cobalt based metal alloy with high tensile strength. It is composed of 42.5% cobalt, 13.0% nickel, 20.0% chromium, 2.0% molybdenum, 0.2% carbon, 0.04% beryllium, 1.6% manganese, 2.8% tungsten and remainder iron, see Hamilton Precision Metals (2017). The target body (10) is composed of standard niobium. Target loading and unloading is performed using an automated loading system, see Hoehr et al. (2014).
In the solid target assembly, the foil target (11) is in the place of the HAVAR with helium jets for cooling on both sides (6) and (12). Due to the use of thin foils, the proton beam traverses through the target and is finally stopped by the water cooled aluminium block (13) which acts as the beam dump. The geometries were modelled as accurately as possible by using dimensions from technical drawings.
The nuclear and chemical properties of the liquid and solid target materials are listed in Table 1. After irradiation, isotopic yield measurements were performed using gamma-ray spectrometry analysis or ionization chamber measurements. All measured yields were decaycorrected to the end of bombardment (EOB). When multiple irradiations took place for the same isotope, the yield was normalized to the beam current prior to calculating the average saturation yield. The error in the yield is dominated by the standard deviation of the different irradiations. For more details see Infantino et al. (2016) and references therein.

Monte Carlo simulations
Geant4 is an all particle Monte Carlo toolkit designed for simulating particle interactions from 100 TeV down to a few eV. Geant4 is implemented in C++ and has great flexibility and expandability and thus is used in various applications such as space research, Large Hadron Collider (LHC) experiments, medical physics or microdosimetry applications (Agostinelli et al., 2003;Allison et al., 2006Allison et al., , 2016.

Saturation yields
In Geant4, the calculation of induced activity relies on the cross section library used for the inelastic nuclear reactions. These cross sections are included in the TENDL1.3 package and can directly calculate the number of isotopes produced. The production rate for each isotope is simulated taking into account primary proton impact, secondary interactions and decay of other isotopes produced in these interactions. Geant4 calculates the isotope production from the primary particle induced production and also the full Bateman solution considering the breeding of radioactive decay products. The number of isotopes at any given time t during the irradiation is given by: where n is the number of isotopes produced per unit mass and unit time (a function of the proton flux, the target density, and the nuclear cross section), I is the proton beam current from the accelerator, n p is the number of incident protons, e is the proton electric charge, and λ is the decay constant of the isotope, see Bungau et al. (2014). This equation reaches a saturation level for long irradiation, When calculating yield ratios the experimental and MC uncertainties have been added in quadrature.

Target geometry and material definition
The solid and liquid targets have been represented in Geant4 using two geometries as shown in Figs. 1 and 2. Geant4 provides a wide range of simple solid geometries that can be used. More complex geometries such as the conical collimator can be generated by combining existing shapes with Boolean operators such as G4UnionSolid and G4Sub-tractionSolid.
The target materials have been divided into two categories: liquid target, containing water solution of salts, and solid targets. While it is possible to use the natural isotopic composition of elements from the NIST database, user defined isotopic compositions were used in order to

T. Amin et al.
Applied Radiation and Isotopes 133 (2018) 61-67 match material definitions in the FLUKA model in Infantino et al. (2016). For solid and liquid targets, the mass fractions were calculated for each element and used in the definition of materials.

Physics models
Geant4 provides multiple (data-driven, parametrized and theorydriven) physics models, each applicable for different particle interactions at different energy levels. In this study, in order to model proton (and neutron) inelastic hadronic interactions in the relevant energy range, three physics lists were considered: Bertini Intranuclear Cascade High Precision (QGSP-BERT-HP) model, Binary Intranuclear Cascade High Precision (QGSP-BIC-HP) model and Binary Intranuclear Cascade All High Precision (QGSP-BIC-AllHP) model. In Geant4 QGSP-BERT-HP and QGSP-BIC-HP are well established physics list for low energy application but were not developed for predicting radionuclide production. From the three physics lists investigated, QGSP-BIC-AllHP proved to be the best approximator for our investigation. The QGSP-BERT-HP list failed to calculate any yield for 13 N and 61 Cu. QGSP-BIC-HP did not calculate any 13 N yield. Due to these limitations, mainly results from the QGSP-BIC-AllHP physics list are being discussed in this paper.
QGSP-BIC-AllHP is a new data-driven all particle, high precision physics model that uses TALYS-based Evaluated Nuclear Data Library (TENDL). TENDL is based on experimental and calculated results of TALYS nuclear model code to produce a nuclear data library for Alpha, Deuteron, 3 He, Proton, Neutron and Triton for energies below 200 MeV. The proton sub-library contains cross sections of about 2800 isotopes. This model has been validated against experimental data (Koning and Rochman, 2012). In this work TENDL 2015 cross sections were used with Geant4 10.1 for energies below 200 MeV.
To describe electromagnetic interactions, electromagnetic options 1, 2 and 3 were tested. Electromagnetic option 1 proved to produce comparable results with the benefit of reduced computation time. The production thresholds were set at 1 mm for all particles inside the target volume.

Proton beam and scoring
In this work, isotopic yields have been normalized to beam current on the target. Since collimated protons do not contribute to yields, but consume simulation time, a idealized pencil beam was used in simulations. The parameters scored in the simulation were secondary nucleons produced due to inelastic protons interactions, N and number of protons incident on the solid or liquid target, n p . For isotopes in excited states, the yield is presented as the sum of the metastable and ground states. Inside the target volume, 100 μm and 1 μm binning was used for the liquid and solid target respectively. To achieve statistical uncertainties in the yield of less than 1%, for the isotopic yield of interest, numbers of primaries simulated were between − 10 10 9 1 0 .

External cross sections
As the saturation yield is a function of the nuclear reaction cross section in the energy range between the beam entering the target to the beam exiting the target or being stopped, see Infantino et al. (2016), comparing the area under the cross section in this energy range between experimental and TENDL cross sections is a good measure of the expected yield difference. Experimental Nuclear Reaction Data cross sections were taken from EXFOR (IAEA, 2017). For reactions with multiple available sources, selections were performed taking into account error margins and the number of data points available for the energy range concerned. The source(s) of experimental reaction cross sections for every isotope investigated are listed in Table 2 under the reference column. After selecting appropriate cross sections, a curve was fitted through the cross sections, and the area under the curve was calculated for both the EXFOR and the TENDL cross sections. Comparisons between TENDL and EXFOR cross section areas are shown in Table 2.

FLUKA
FLUKA is a Fortran-based general purpose Monte Carlo code used to investigate particle transport and interaction with matter (Ferrari et al., 2005;Boehlen et al., 2014). It is applicable at energies from low (keV) energy to TeV energy levels such as shielding, target design, calorimetry, hadron therapy, neutrino physics, cosmic rays, etc. FLUKA is jointly developed by the European Organization for Nuclear Research (CERN) and the Italian Institute for Nuclear Physics (INFN). The FLUKA MC package version 2011.2b.6 was used for the isotope production at the medical cyclotron. Isotope production in FLUKA is handled inside the software package. Materials were defined by the user to match with experimental details. For more details, see Infantino et al. (2016).

Experimental reference data
Experimental yields including the reference source, and Monte Carlo values used for comparison are listed in Table 2. For each isotope, the saturation yield Y exp and the number of irradiations conducted are reported. The produced activity in each irradiation is normalized to the beam current before calculating the average saturation yield of the given isotope. The yield has been decay-corrected to EOB and normalized for a beam current of 1 μA incident on the target for a 1 h irradiation.

Direct assessment from Geant4
The isotopic yields calculated in Geant4 are compared with measurements in Table 2. Y F refers to yield from FLUKA, Y G refers to Geant4 yields using TENDL libraries. Table 2 compares the simulated and experimental yields. The table also lists cross section ratios obtained from Fig Fig 3, where the TENDL and EXFOR cross sections (X T , X E ) have been expressed as ratios of each other. This has been referred to as cross section ratio later in the paper. The cross section ratio can be used to get an approximation of the theoretical yields that can be expected when using the TENDL library. In Fig. 3 only the most probable nuclear reaction cross section was taken into account. However, in MC code incident protons and all secondaries are taken into account. For the reaction nat Sr(p,x) 86 Y, the cross section for only 86 Sr(p,n) 86 Y reaction was taken into account as it was the majority contributor to 86 Y yield. Contributions from other reactions are assumed to be insignificant to the overall yield and hence not taken into account. During experiments or routine isotope production, there are losses in the transfer system and in vials prior to measurement for liquid targets and dissolved solid targets respectively. Also MC codes do not take into account complex thermal and fluid dynamics of the liquid target. Thus a factor of 2 seems to be an acceptable limit for the ratio of saturation yield. The comparison between Geant4 and experiment for the isotopes 18 F, 44 Sc, 52 Mn, 55 Co, 61 Cu, 68 Ga, 89 Zr, and 94 Tc fulfills this criteria. Only for 13 N and 86 Y is the ratio between Geant4 and experiment larger than 2, and none is smaller than 0.5. Overall in this section, Geant4 is less than a factor of two away from the experimental yield. It is also closer to the experimental yield than FLUKA for five isotopes ( 18 F, 52 Mn, 55 Co, 61 Cu, 68 Ga), while FLUKA is closer for three isotopes ( 68 Ga, 89 Zr, 94 Tc). In general the comparison of the EXFOR cross sections with the TENDL cross sections used in Geant4 are within 25% except for 61 Cu (0.55).
While the yield of 18 F is under-calculated by a factor of 0.53 using Geant4, FLUKA over-estimates it by a factor of 1.66. The database takes into account multiple sources to provide a single unified table of cross sections that has been used. The 18 O(p,n) 18 F reaction has multiple resonances between 2 and 10 MeV with each experiment reporting slightly different peaks. This does not take into account the efficiency of the liquid target system. For 52 Mn, 55 Co and 61 Cu Geant4 performed better than FLUKA with ratios of 1.1, 0.7 and 0.6 against FLUKA's 4.62, 0.3 and 3.13 respectively. For these solid targets FLUKA appears to be less reliable than Geant4, with all yield ratios outside acceptable limits. For these three isotopes, the yield ratios of Geant4 to experimental values correlate very well to the cross section ratios between TENDL and EXFOR. For 52 Mn the yield ratio is 1.1 while the cross section ratio is 0.93, for 55 Co yield and cross section ratios are 0.7 and 0.88 respectively. 61 Cu has a yield ratio and a cross section ratio of 0.6 and 0.55 respectively. The excellent level of agreement between the yield ratio and cross section ratio indicates that while the Geant4 yield might be different from experiments, it is a consequence of mismatching TENDL and EXFOR cross sections.
For 68 Ga FLUKA performed better with a ratio of 1.03 against Geant4's 0.84. The cross section of the 68 Zn(p,n) 68 Ga reaction currently has significant discrepancies, hence multiple sources were taken and a spline fit was used to make comparisons. The TENDL library underestimates yields over the concerned energy range compared to fitted EXFOR with a ratio of 0.75. Geant4 over-calculates the yield of 86 Y by a factor of 2.5 whereas the FLUKA yield ratio was 0.9.
The yield for 89 Zr was calculated more accurately using FLUKA than Geant4, the respective yield ratios are 0.87 and 0.69 respectively. Both MC codes under-estimate the yield, with Geant4's performance disagreeing with theoretical expectations. The TENDL cross section is higher than most EXFOR tabulated cross sections. This indicates that Geant4 should calculate a yield higher than experiments, however, the yield from both MC codes is lower than that of experiments. At this moment no explanation has been found why the MC results challenge the cross sections available.
Due to Geant4 and FLUKA's inability to calculate metastable isotopes, 94 m Tc is presented as the sum of metastable and ground state. For this isotope, Geant4 calculates the yield with a factor of 1.7 whereas FLUKA ratio is 1.53 and the cross section ratio is 1.16. Both MC are able to calculate accurately 94 Tc yield, with FLUKA performing slightly better.  Zhuravlev et al. (1994) 3.2.2. 13 N, 44 Sc, and 86 Y For these isotopes the deviation from the experiment is larger than a factor of two. For 13 N and 44 Sc the deviation in the Geant4 simulation is smaller than for FLUKA, while only for 86 Y is the deviation for Geant4 larger than for FLUKA. The EXFOR cross section area is the same as the TENDL cross section area, except for 13 N which has a very large ratio of 2.34. 13 N yield was overestimated by a factor of 2.72 compared to a factor of 5.9 from FLUKA. The yield between Geant4 and experiments of 13 N results are not comparable at TR13 energy levels as TENDL does not account for the resonance at 7.9 MeV for the nat O(p,x) 13 N reaction. For energies above 8.5 MeV, TENDL cross section vastly over-estimate the yield and has large disagreements with EXFOR. This is illustrated by Fig. 3a. This phenomenon was also observed when 13 N was created inside a PMMA target under proton therapy conditions in Amin et al. (2017). T. Amin et al. Applied Radiation and Isotopes 133 (2018) 61-67 For 44 Sc FLUKA performed worse with a ratio of 2.35 against Geant4's 2.1. When comparing with experimental cross sections, the ratio for 44 Sc is 1.05. While the agreement between ratios is acceptable, EXFOR lacks sufficient good quality cross sections for the reaction of interest at these low energy levels. The contribution of 44 m Sc to the total production of 44 Sc in MC calculations is negligible for TR13 energy ranges.
Geant4 overestimates the yield of 86 Y by a factor of 2.5 whereas the FLUKA yield ratio is a very good 0.9. The yield of 86 Y has been represented here as the sum of metastable and ground states of 86 Y. As a result, a minor overestimation from Geant4 is expected when comparing simulated yields with experimental yields. The fitted tabulated cross sections for 86 Sr(p,x) 86 Y reaction were provided by the IAEA. The fit was performed using data from Roesch et al. (1993) and Levkovsky et al. (1990a), where the former had significant error bars contributing to a slightly inaccurate smoothing of the fit. Compared to EXFOR, the TENDL data had a marginally larger overall yield in the energy ranges relevant to this work. Despite the discrepancy between the MC codes, the ratio of Geant4 to experimental yields agrees well with the ratio of cross sections for 86 Y, as shown is Table 2.

Conclusions
A Geant4 simulation to model the liquid and solid target assembly for the TR13 medical cyclotron at TRIUMF has been developed. The agreement between Monte Carlo simulation and experimental yield measurements varies depending on the isotope considered. The physics list QGSP-BIC-AllHP in Geant4 was investigated in this study, a new list recently released with version 10.1. Its performance depends almost entirely on the accuracy of the TENDL cross sections utilized by the user. In our work, TENDL has proven to provide accurate cross sections for certain reactions, whereas for example the 16 O(p,x) 13 N cross section currently is incorrect. For certain isotopes such as 68 Ga, 86 Y, 89 Zr and 94 Tc FLUKA was better able to calculate yield. For 13 N, 18 F, 44 Sc, 52 Mn, 55 Co and 61 Cu Geant4 performed better, despite the MC models not accounting for thermal effects or density changes in the liquid target or loss in the transfer system. Overall, in our situation using Geant4 10.1 with the physics list QGSP-BIC-AllHP the mean absolute deviation for all targets is ± 0.7 0.5 compared to ± 1.4 1.6 for FLUKA. In addition, QGSP-BERT-HP and QGSP-BIC-HP in Geant4 were also investigated. The QGSP-BERT-HP list managed to produce a slightly lower mean absolute deviation of ± 0.6 0.4, but failed to calculate any yield for 13 N and 61 Cu. QGSP-BIC-HP had a mean absolute deviation of ± 1.1 1.2, failing to calculate any 13 N yield. Due to these limitations, neither were further considered. For some isotopes differences are still large. A wider range of isotopes needs to be examined for a better assessment.