Antiproton annihilation physics in the Monte Carlo particle transport code SHIELD-HIT12A

Nuclear Instruments and Methods in Physics Research B 347 (2015) 65–71

Contents lists available at ScienceDirect

Nuclear Instruments and Methods in Physics Research B journal homepage: www.elsevier.com/locate/nimb

Antiproton annihilation physics in the Monte Carlo particle transport code SHIELD-HIT12A Vicki Trier Taasti a, Helge Knudsen a, Michael H. Holzscheiter a,b, Nikolai Sobolevsky c,d, Bjarne Thomsen a, Niels Bassler a,⇑ a

Dept. of Physics and Astronomy, Aarhus University, Denmark Dept. of Physics and Astronomy, University of New Mexico, USA Institute for Nuclear Research of the Russian Academy of Sciences (INR), Moscow, Russia d Moscow Institute of Physics and Technology (MIPT), Dolgoprudny, Russia b c

a r t i c l e

i n f o

Article history: Received 11 November 2014 Received in revised form 30 January 2015 Accepted 1 February 2015

Keywords: Antiprotons Monte Carlo simulations SHIELD-HIT12A AD-4/ACE Annihilation physics

a b s t r a c t The Monte Carlo particle transport code SHIELD-HIT12A is designed to simulate therapeutic beams for cancer radiotherapy with fast ions. SHIELD-HIT12A allows creation of antiproton beam kernels for the treatment planning system TRiP98, but first it must be benchmarked against experimental data. An experimental depth dose curve obtained by the AD-4/ACE collaboration was compared with an earlier version of SHIELD-HIT, but since then inelastic annihilation cross sections for antiprotons have been updated and a more detailed geometric model of the AD-4/ACE experiment was applied. Furthermore, the Fermi–Teller Z-law, which is implemented by default in SHIELD-HIT12A has been shown not to be a good approximation for the capture probability of negative projectiles by nuclei. We investigate other theories which have been developed, and give a better agreement with experimental findings. The consequence of these updates is tested by comparing simulated data with the antiproton depth dose curve in water. It is found that the implementation of these new capture probabilities results in an overestimation of the depth dose curve in the Bragg peak. This can be mitigated by scaling the antiproton collision cross sections, which restores the agreement, but some small deviations still remain. Best agreement is achieved by using the most recent antiproton collision cross sections and the Fermi–Teller Z-law, even if experimental data conclude that the Z-law is inadequately describing annihilation on compounds. We conclude that more experimental cross section data are needed in the lower energy range in order to resolve this contradiction, ideally combined with more rigorous models for annihilation on compounds. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Since 2002 the AD-4/ACE research collaboration at CERN has investigated whether antiprotons provide benefits for radiation therapy compared to protons [1–4]. As described in [5] the antiproton depth dose curve is quite similar to the one for protons, but an extra dose is deposited in the Bragg peak due to the antiproton annihilation products. When an antiproton comes to rest, it is captured by a nucleus and here it annihilates on either a proton or a neutron, whereby an energy of 1.88 GeV (two times the rest mass of the proton) is released. The annihilation energy is typically converted into the rest mass of 4–5 pions – a mixture of pþ , p and p0 ⇑ Corresponding author. E-mail address: [email protected] (N. Bassler). http://dx.doi.org/10.1016/j.nimb.2015.02.002 0168-583X/Ó 2015 Elsevier B.V. All rights reserved.

to satisfy charge conservation – and most of the surplus annihilation energy is carried away by the pions as kinetic energy. A minor fraction of the annihilation energy,  30 MeV , is deposited locally and causes a doubling of the dose in the Bragg peak compared to a proton beam [6]. This local dose deposition comes from short-range ions created by nuclear fragmentation of the target atom, since on average 1–2 of the created pions will penetrate the remaining nucleus causing it to break up. Even if no pion is absorbed by the target nucleus, the instantaneous removal of a nucleon causes the target nucleus to recoil with a kinetic energy of several MeV. Besides the extra physical dose for antiprotons in the Bragg peak compared with protons, a higher RBE (Relative Biological Effectiveness) is also expected here, since the annihilation fragments are assumed to have a higher LET (Linear Energy Transfer).

66

V.T. Taasti et al. / Nuclear Instruments and Methods in Physics Research B 347 (2015) 65–71

The heaviest fragments contribute only little to the physical dose whence the dose averaged LET, and therefore the RBE, is dominated by the contribution from the lighter ones [7]. Antiprotons deposit a slightly higher dose in the plateau region compared to protons. This is caused by in-flight annihilation of the antiprotons, i.e. some antiprotons will collide with a target atom before coming to rest and then they may annihilate on this atom. However, due to the much higher dose deposition from the antiproton annihilation in the Bragg peak, the peak to plateau ratio of antiprotons is improved by a factor of 2 or more in a spread out Bragg peak [8]. Furthermore, an additional dose is deposited after the Bragg peak due to the annihilation products. These are distributed isotropically and they give a dose beyond the tumor region as well as a wider penumbra compared to a proton beam [9]. Pions, energetic protons and neutrons have low collision cross sections and consequently long ranges. They can therefore lead to dose depositions even far from the tumor, which is a potential drawback. However, the particles escaping the body can be used for imaging; the annihilation point can be determined by extrapolating the particle trajectory in a detector back to its origin [10]. SHIELD-HIT12A is a Monte Carlo particle transport code, which was specifically developed for transport of ions in biological tissues [11–13] for research purposes [14,15], and is based on earlier versions as presented in [16,17]. The aim of this work is to benchmark SHIELD-HIT12A to simulate the antiproton depth dose curve in water correctly. SHIELD-HIT12A provides an interface to the treatment planning system TRiP98 [18]. This way SHIELD-HIT12A can be used as a base for creating treatment plans for antiproton therapy, once the radiobiological effects have been assessed by the AD-4/ACE collaboration. Along with a benchmarked RBE model, such treatment plans can be used to realistically examine the clinical consequences of antiproton therapy by comparing the efficacy with other beam modalities. 2. Theoretical models and scaling laws

2.1.1. The Fermi–Teller Z-law Fermi and Teller [19] argue that the capture probability is proportional to the energy loss of the ‘‘negative mesotron’’ near the atom, which is again proportional to the number of atomic electrons. They therefore suggest

rðZÞ / Z:

ð4Þ

This is the so-called ‘‘Z-law’’. However, the Z-law has later been shown to be a poor approximation to real capture cross sections. Baijal et al. [20] measured a substantial number of capture probability ratios WðZÞ=WðZ 0 Þ for

l projectiles. Their results corresponded to rðZÞ / Z k , with values of k varying from 1 to +1.5. Further, since the Z-law does not take the chemical structure into account, it predicts rðZÞ, and thereby also AZ , to be independent of the type of bonds in the molecule. Zinov et al. [21] measured the capture probability ratio for metal oxides. They found that in different oxides containing the same atoms, but in different ratios, the capture probabilities varied significantly. This reveals a dependence of the capture probability on the chemical composition, thus proving the inadequacy of the Z-law. Plotting the capture probability ratios for oxides a variation following the periods of the periodic table was evident. Krumshtein et al. [22] measured the probability for negative pions to be captured to the hydrogen nuclei of hydrogen-containing compounds Z m Hn . Here the Z-law predicts

WðHÞ ¼

n n  ; n þ mZ mZ

ð5Þ

but the proper fit to their data showed a much larger suppression of the capture to hydrogen. Their results rather gave

WðHÞ 

n mZ 3

:

ð6Þ

For water, they find WðHÞ ¼ 3:5  0:6  103 , which is much smaller than the 0.2 expected from the Z-law. In conclusion the Fermi–Teller Z-law does not hold in general.

2.1. Capture of very slow antiprotons in compounds By default in SHIELD-HIT12A the probability for annihilation on a specific atom in a compound target is based on the Fermi–Teller Z-law, but as explained below, this law has been shown to be inconsistent with experimental findings. A literature survey has therefore been performed and a replacement for the Fermi–Teller Z-law has been identified and implemented in the code. We will now discuss the information available on the capture of , by heavy, negatively charged particles, such as l , p , K  and p the nuclei of the material in which they come to rest. We regard a compound Z 0n Z m and wish to investigate the probability W(Z) for capture to take place on the nuclei with charge Z. The rest of the particles are then captured to the other nuclei, therefore we must require

WðZÞ þ WðZ 0 Þ ¼ 1:

ð1Þ

We may define a ‘‘cross section for capture’’ integrated over the appropriate antiproton energies, rðZÞ. We then have

WðZÞ ¼

mrðZÞ ; mrðZÞ þ nrðZ 0 Þ

ð2Þ 0

where m and n is the number of Z and Z atoms in the compound. From this the ratio of the capture probabilities per atom can be defined as

AZ ðZ; Z 0 Þ 

nWðZÞ rðZÞ : ¼ mWðZ 0 Þ rðZ 0 Þ

ð3Þ

2.1.2. Capture to protons in hydrogenous compounds The strong suppression of capture to hydrogen is explained by Ponomarev [23]. His theory gives a way of calculating the capture probability to hydrogen in molecules containing hydrogen, Z m Hn . A base for Ponomarev’s theory is the experimental finding that in chemical compounds the mesons are often captured to the outer molecular orbits rather than to the individual levels of the single atoms. The theory proposes the following model for calculating WðHÞ: in order to have annihilation on hydrogen, the heavy, negative particle must be captured to a common molecular valence state by replacing an electron (the only electron it can replace around the hydrogen nucleus is a valence electron) – this happens with probability W 1 . Then the heavy, negative particle valence state must decay to a bound hydrogen state, probability W 2 , and the following cascade towards the nucleus must not involve transfer of the meson to a state around the Z-nucleus, probability W 3 . Therefore, the total probability for annihilation on H is

WðHÞ ¼ aW 1 W 2 W 3 ;

ð7Þ

where a is a coefficient which takes into account the features of the chemical bond of hydrogen in the molecules. Krumshtein et al. [22] measured WðHÞ for a number of hydrogen compounds and found a ¼ 1:3 to give a good fit to the results. They suggested that in general a can be stated as

a ¼ 2s þ 1:3ð1  sÞ;

ð8Þ

67

V.T. Taasti et al. / Nuclear Instruments and Methods in Physics Research B 347 (2015) 65–71

where s is the degree of hydrogen bond ionicity. Here the first term should be abandoned if the hydrogen electronegativity is smaller than that of the Z-atom. There are 2n valence electrons and n þ mZ electrons in total, so if it is assumed that all the molecular electrons are equivalent, W 1 is given by

W1 ¼

2n : n þ mZ

argues that this probability is proportional to Z 2 , and he gets

W 2  Z 2 :

ð10Þ

When discussing W 3 , it should be remembered that if a heavy, negative particle-hydrogen atom is created, it can interact with its own Z-atom, or with other Z-atoms. This means that W 3 depends on the target density, in case we have a gaseous target. For solids or liquids, the density is not a parameter. Ponomarev gives arguments that W 3 does not depend on Z, and ignores the effect, letting W 3 ¼ 1, which means that no transfer happens. Therefore his result is

a n : Z 2 n þ mZ

ð11Þ

This is in agreement with the empirical result in Eq. (6). For more complex molecules Z 0k Z m Hn Eq. (11) can be expanded into

WðHÞ ¼

Projectile

Pt



p

0.9 0.7 0.0

K  p

ð9Þ

A heavy, negative particle captured into a valence state will have the same binding energy as that of the replaced electron, and therefore also an orbital with the same extension. The likelihood for Auger transitions with another electron is proportional to the spatial overlap of the two wave functions, and Ponomarev

WðHÞ ¼

Table 1 Transfer probabilities, P t , for the negative projectile particles.

amZ 2 þ a0 m0 Z 02 0 : n þ mZ þ kZ

ð12Þ

2.1.3. Capture in compounds not containing hydrogen For compounds not containing hydrogen, the Ponomarev theory cannot be used. Schneuwly et al. [27] introduced a detailed model which takes into account the different chemical structures of the compounds. The same was done in another and more simple way by Daniel [28]. Since this work is aiming to obtain a fairly simple algorithm for calculation of the capture probability, focus will be on the simpler model by Daniel. He proposed an extension of the work by Fermi and Teller [19]. Here, the heavy, negative particle travels in the screened atomic potential

UðrÞ ¼ 0:354 a0

e2 Z 2=3 r2

for r < r 0 ;

ð14Þ

where r 0 is a cut-off distance taken to be the same for all atoms, and approximately equal to the distance between the atoms. The projectile loses energy to the outer (valence) electrons, which have a density like a Fermi gas. Daniel used Fermi and Teller’s calculation of the energy loss DEðTÞ and found

AZ ðZ; Z 0 Þ ¼

Z 1=3 lnð0:57ZÞ Z 01=3 lnð0:57Z 0 Þ

ð15Þ

;

where the constant 0.57 does not depend on the projectile mass. Later Daniel [29] modified his result, Eq. (15), to take into account that the cut-off distance in his model is not the same for all the target atoms. He introduced the real atomic radii R(Z)

Z 1=3 lnð0:57ZÞRðZ 0 Þ

Here m and m0 are the number of valence hydrogen bonds with the Z and Z 0 atoms, respectively.

AZ ðZ; Z 0 Þ ¼

Comparison between Ponomarev’s theory and antiproton measurements: Pawlewicz et al. [24] measured WðHÞ for antiprotons stopped in propane (C3H8) and found WðHÞ ¼

Daniel observed the result (15) to fit experimental negative muon data for 27 different metal halogens better than the Z-law. Von Egidy et al. [30] measured AZ ðZ; Z 0 Þ for negative muons stopped in 57 different oxide targets. They found that Daniel’s result (15) gives a good, simple, analytic approximation to all data for oxides. A better fit to AZ ðZ; Z 0 Þ will probably be obtained with Eq. (16), but to get a general result out of this equation, large tables of ionic radii are needed.

ð110  30Þ  103 . This corresponds to a value of a ¼ 13. Similar measurements on propane with K  gave WðHÞ ¼ ð36  4Þ  103 [25,26]. This dependence on the projectile may be explained by the fact that Ponomarev’s theory ignores the transfer probability, P t ¼ 1  W 3 . According to Pawlewicz et al. [24], this shortcoming can be remedied by replacing Eq. (11) by

WðHÞ ¼ ð1  Pt Þ

AL n ; Z 2 n þ mZ

ð13Þ

where P t is the probability for transfer of the negative particle to the Z-atom during its cascade towards the proton. Since the cascade time of fairly light particles such as l and p is considerably  it is easy to imagine that Pt ðl Þ  Pt ðp Þ. Pt ðp Þ is longer than for p therefore set to zero. Pawlewicz et al. suggests that if we use Eq. (13) for all heavy, negative projectiles, applying AL ¼ 13 overall and the values of Pt listed in Table 1, then Eq. (13) agrees with ‘‘all’’ data for Z m Hn targets. For water, we get WðHÞ ¼ 4:0%. This is far from the result of the Z-law which predicts 20%. As far as there seems to be no other measurements of WðHÞ for antiprotons, we are left with Eq. (13) combined with Pawlewicz’s suggested values for AL and P t as our best educated guess for this case. Clearly, further investigations are needed.

Z 01=3 lnð0:57Z 0 ÞRðZÞ

:

ð16Þ

2.1.4. The PPD-law Based on a multitude of experimental data for negative muon and pion capture, and a few data for kaon and antiproton capture, we conclude: 1. For compounds not containing hydrogen, the formula by Daniel (15) gives a suitable fit to the data, very likely accurate within 10–20%. 2. For hydrogen compounds, we should instead use the formula (13) by Ponomarev and expanded by Pawlewicz et al., which is accurate within 30%. 3. The ‘‘Z-law’’ is wrong in all cases but especially for capture to hydrogen, where it often is wrong by an order of magnitude or more. It should be emphasized that the conclusions for antiproton capture are based on an expansion of the Ponomarev model by Pawlewicz et al. to antiprotons, based on one measured datum (for propane). That Daniel’s formulae can be used for antiprotons

68

V.T. Taasti et al. / Nuclear Instruments and Methods in Physics Research B 347 (2015) 65–71

is inferred by the fact that it does not depend on the projectile mass. In the following the combination of the laws by Ponomarev/ Pawlewicz and by Daniel used in the implementation in SHIELDHIT12A will be referred to as the PPD-law. In this implementation the following assumptions are made: 1. Eq. (12) can be extended to apply for compounds containing more than three elements, and the Pawlewicz correction can be used here as well. 2. A constant value of AL ¼ 13 can be used for all compounds. 3. Calculating the capture probability for elements different from hydrogen in hydrogen-containing compounds the Daniel law can be applied, when a proper normalization is used.

the ratio inelastic/total is the same as for protons. For the remaining atoms in the periodic table an interpolation of the form A2=3 is used. Furthermore, the cross sections in [31] are only tabulated for energies above 10 MeV. An extrapolation to lower energies is performed in the code, but since the residual ranges of these antiprotons are negligible this does not introduce any relevant uncertainty. The accuracy of the Sychev tables is also given and are of the order of 3–10%. The cross sections for interactions of antiprotons with nuclei govern how many antiprotons will reach the Bragg peak region. Large cross sections give a high probability for in-flight annihilation and a smaller energy deposition in the Bragg peak. Thus by scaling the cross sections it can be controlled where the energy deposition takes place. This gives one way of fitting the simulated depth dose curve to the experimental data.

2.2. Implementation of the PPD-law The new calculation for the capture probabilities is implemented in the following way: Should the compound target medium not contain hydrogen, then Daniel’s simple formula, Eq. (15), is interpreted

WðZÞ ¼ Z 1=3 lnð0:57ZÞNðZÞ;

ð17Þ

where N(Z) is the partial number concentration of the Z-atom. Should the target compound contain hydrogen the capture probability to hydrogen is calculated

8 13 n > < ð1  Pt Þ Z2 nþmZ P j WðHÞ ¼ m Z 2 i¼1 i i > : ð1  Pt Þ 13 P j nþ

i¼1

mi Z i

for Hn Z m :

ðaÞ

for Hn Z 1;m1    Z j;mj : ðbÞ

ð18Þ

With the transfer probabilities, Pt , given in Table 1. In the case Eq. (18b) the target medium is a compound which contains more than one element besides hydrogen. The number of bonds between hydrogen and the other elements, mi , is estimated in a precalculation. A routine for calculation of mi is built by making an initial calculation of the maximal number of bonds for each Z i atom, based on their highest possible oxidation state. Then these numbers are P reduced until the constraint n ¼ mi (that is the sum of bonds to hydrogen equals the number of hydrogen atoms) is met. The mi values are first reduced for the least electronegative atom and later for the atom initially having the highest number of bonds to hydrogen. The reduction procedure is iterated until the constraint is fulfilled. The capture probabilities to the Z-atoms is found by

WðZÞ ¼

8 1  WðHÞ > > > < 1=3 > > > :

ð1WðHÞÞZ j X 1=3

Zi

for Hn Z m : lnð0:57ZÞNðZÞ

ðaÞ

for Hn Z 1;m1    Z j;mj : ðbÞ

ð19Þ

lnð0:57Z i ÞNðZ i Þ

i¼1

2.3. Antiproton collision cross sections In SHIELD-HIT12A the collision cross sections for different hadrons with nuclei are given in the form of a table for the inelastic and total cross sections, respectively. These table values are then interpolated for the given target nucleus and projectile energy. The cross sections for antiprotons come from tabulated data given in the book by Sychev [31]. The tabulated data for antiprotons listed in this book only include the total and the elastic cross sections for the proton–antiproton interactions. Inelastic cross sections are given for an antiproton interacting with the nucleus of 12 other elements. For these 12 elements the total cross sections have been tabulated in the source code of SHIELD-HIT12A applying a scaling of the inelastic cross sections, under the assumption that

3. The experiment The experiment is described in [32] and there all necessary experimental details can be found. A beam of 502 MeV/c antiprotons left the vacuum beam line through a 0.025 mm thick titanium window and was directed at a water tank made of PMMA. On its way the antiproton beam passed through normal air and through an ionization chamber used for normalization of the pulse fluctuations, referred to as ‘‘the monitor chamber’’. The beam entered the water tank through a PMMA window having a thickness of 3 mm and a diameter of 70 mm. In the water tank the dose was measured with a second ionization chamber (‘‘the measuring chamber’’) which recorded the depth dose curve. Both ionization chambers where of the ‘‘Advanced Roos’’ type from PTW Freiburg with graphite electrodes. The Advanced Roos chamber is similar to the commercially available plane parallel PTW Roos chamber, but has an effective diameter of 3.96 cm. Both ionization chambers were calibrated absolutely in photon beams. In this work we reanalyzed the data with the aim of absolute quality, instead of just relative as in [32]. Again we follow the IAEA TRS-398 ion dosimetry protocol as closely as possible [33]. From the collected charge the deposited dose-to-water can be calculated using the standard formalism

Dw;Q ¼ M Q ND;w;Q 0 kQ ;Q 0 :

ð20Þ

The calibration factor for dose-to-water was N D;w;Q 0 ¼ 1:32  107 Gy/C, with an uncertainty of 4% (see [34]), and the beam quality correction factor kQ ;Q 0 , which takes into account the differences between measurements using the reference 60Co beam and an antiproton beam, was set to 1.00 as explained in [32]. M Q is the reading of the dosimeter times a product of correction factors for influence quantities such as pressure, temperature, polarization voltage and collection efficiency, Pk [33]. Here a correction for the collection efficiency and the pressure is needed. The reading from the chamber is corrected for recombination losses using the Boag Two Voltage Method, and the pressure correction is given by the ratio PP0 , where P 0 is the atmospheric pressure during chamber calibration (1013 hPa) [34]. Eq. (20) gives the total dose stemming from all antiprotons in one spill, but in order to compare the measured dose with the simulated dose, the number of antiprotons exiting the beam nozzle then has to be known. Since the monitor chamber is calibrated in absolute terms, we know the deposited dose. Using a model of the experimental setup in SHIELD-HIT12A we can relate the measured dose in the monitor chamber with the actual number of antiprotons which left the antiproton decelerator.

V.T. Taasti et al. / Nuclear Instruments and Methods in Physics Research B 347 (2015) 65–71

One may dispute finding the number of antiprotons in the experiment based on the simulations, which are to be tested against the experimental data, as this could be considered circular. However, it is the simulated dose depositions in the measuring chamber which are to be compared with the experimental data, while it is the simulated dose in the monitor chamber which is used for finding the number of antiprotons in the experiment. The models to be investigated in this work, the Fermi–Teller Zlaw and the antiproton collision cross sections, only have a significant influence in the Bragg peak and very little influence on the monitor chamber. Independent checks confirmed that the dose measured in the monitor chamber is almost entirely based on the stopping power of the antiprotons alone, and therefore we apply this method given the lack of alternative absolute dose monitoring. A coarse estimate of the number of antiprotons per spill is provided by a beam transformer located a few meters upstream the exit window. Our estimate of the number of primary particles agreed within  2—4%, which is comparable to the calibration uncertainty of the ionization chambers.

69

Fig. 1. Measured and simulated absolute dose-to-water. Each simulated point is the result of a full simulation including a complete modeling of the measuring ionization chamber. The simulations use the default Z-law.

3.1. Stopping power ratios In the original simulation approach presented in [32], the dose was scored directly in the water tank without modeling the ionization chamber. The dose was scored in cylindrical bins with a thickness of 0.01 cm and with a radius of 1.98 cm, corresponding to the radius of the sensitive volume of the Advanced Roos ionization chamber. The dose simulated is therefore in terms of dose-towater. Contrary to [32], here we have modeled the experiment and ionization chamber design precisely. The Advanced Roos ionization chamber is a PMMA encapsulated air cavity with graphite electrodes, where the exact geometry was kindly provided by the manufacturer PTW Freiburg. The dose was scored inside the air gap of the sensitive volume of the chamber. The simulated dose is returned in terms of dose-to-air. Bragg–Gray cavity theory was applied in order to convert dose-to-air into dose-to-water, summing over all charged ions crossing the sensitive volume as simulated by SHIELD-HIT12A. Secondary electrons are not considered, as SHIELD-HIT12A does not provide electron transport. SHIELDHIT12A provides a feature which allows to score fluence multiplied with stopping power for any user-specified medium during runtime, which eases the calculation of stopping power ratios. SHIELD-HIT12A has previously been used for assessing stopping power ratios of various ion beams for particle therapy [17,35]. Since SHIELD-HIT12A simulates all relevant secondaries, we do not need to apply fluence corrections factors as described in [36].

4. Results 4.1. The effect of the geometrical description First a simulation of the antiproton depth dose curve was performed with the current version of SHIELD-HIT12A with the default implementation of the Z-law. Each simulated point is the result of a full simulation including a complete modeling of the measuring ionization chamber at the actual measurement position. This means, the entire simulation must be repeated for each ionization chamber position along the depth dose curve, contrary to the original approach in [32] where the entire depth dose curve could be acquired in a single simulation. It is seen in Fig. 1 that even though, as explained before, the antiproton annihilation physics is not simulated correctly due to

Fig. 2. Comparison of simulations modeling the measuring ionization chamber and scoring the dose directly in water. Both simulations apply the Z-law.

the implementation of the Z-law, the agreement between the simulations and the experimental data is very good. Next a simulation study has been performed to find the consequences of simply scoring the dose directly in the water tank without simulating the real setup of the measuring ionization chamber. In Fig. 2 the results of the two simulation approaches are compared. The effect of using the exact setup of the ionization chamber is clearly visible. In the Bragg peak of the water curve the dose is far larger than in the simulations where the ionization chamber is modeled correctly. The higher dose in the peak of the water curve might be due to the higher atomic weight of oxygen compared to carbon which will result in a different fragmentation spectrum. It can be concluded that water equivalence of the ionization chamber is not a good approximation, when simulating the dose. Rather, the actual composition of the ionization chamber is of importance. 4.2. The PPD-law The Fermi–Teller Z-law is found to give an insufficient description of the capture process, see Section 2.1.1. So to improve the physics implementation in SHIELD-HIT12A the set of equations found to replace the Fermi–Teller Z-law was implemented in the code and the simulations represented in Fig. 1 were redone with the changed source code. The result is presented in Fig. 3. It can be seen from the figure that changing the routine for calculation of the capture probability has an effect on the simulated depth dose curve, but unfortunately the agreement with

70

V.T. Taasti et al. / Nuclear Instruments and Methods in Physics Research B 347 (2015) 65–71

Fig. 3. Comparison of simulation results for the region around the Bragg peak with the unchanged code (solid curve), where the capture probabilities are calculated according to the Fermi–Teller Z-law, and the changed code (dashed curve), where the PPD-law is used.

Fig. 4. Simulation results around the Bragg peak area using the PPD-law and the scaling of the antiproton cross sections by a factor of 1.08.

the measured curve is worsened. The simulations implementing the PPD-law overestimate the measured curve in the top of the Bragg peak. This change is as expected since the PPD-law (in particular the Ponomarev/Pawlewicz law, see Eq. (13)) gives a smaller probability for capture to hydrogen, and therefore leads to an increase of the capture probability to oxygen, when annihilating in water, and to oxygen and carbon, when annihilating in PMMA. This increase gives a higher probability for creation of fragments, since nuclear fragmentation can only happen when the annihilation takes place on a nucleus heavier than hydrogen. In this simulation no annihilations were found to happen in air – annihilations only took place in water or in PMMA. Water and PMMA both contain hydrogen, so the Ponomarev/Pawlewicz law is used for the capture probability calculation. However, PMMA contains two elements besides hydrogen (oxygen and carbon) so the Daniel law is used for calculating the capture probabilities for those elements, see Eq. (19b).

another possibility to improve the fit to our experimental data would be to replace the expression AL ð1  Pt Þ with AL ¼ 13 in Eqs. (13) and (18) by a fitting parameter, a. Doing this, we find a fit to the experimental depth dose curve just as good as the one found by varying the in-flight cross section values when we set a ¼ 35  5. But this approach seems unsound, since if anything, the expression a ¼ 13 ð1  Pt Þ with Pt ¼ 0 suggested by Pawlewicz et al. [24] is an overestimate. This is because Pt is a probability, and hence has values between 0 and 1, and the numerical value 13 is a result of a fit to data of a parameter which originally was estimated by Ponomarev et al. [23] to be ‘‘of the order of unity’’. Ultimately the question also needs to be raised if the experimental data available are of sufficient quality and completeness to reach a final decision. Especially the description of the beam available to us contains a number of assumptions (spot size, emittance, alignment) which could not be verified experimentally, which also have an influence on the calculated depth dose curve. Clearly this calls for more experimental work and while difficult to do at CERN, where the concentration is mostly on ultra-low energy antiprotons and antihydrogen formation, a new antiproton facility, delivering high intensity beam in the energy range of interest to this work, will become available at FAIR at the Gesellschaft für Schwerionenforschung (GSI) in Darmstadt, Germany [37]. Finally, compared to the SHIELD-HIT version used in the publication of [32], a few programing issues where found in the antiproton annihilation model implemented into SHIELD-HIT12A. The final version of SHIELD-HIT12A is tagged SVN-694.

4.3. Scaling of the antiproton collision cross sections In order to retain the implementation of the PPD-law for the capture probabilities, which gives a better description of the physics, it can be tested if a scaling of the antiproton collision cross sections, as described in Section 2.3, will lead to a recovery of the agreement in the Bragg peak. By increasing the antiproton cross sections more antiprotons will annihilate in-flight, not making it to the Bragg peak and this way the dose will be lowered here. The scaling parameter for the cross sections is found by visual inspection of the agreement of the measured and simulated dose. Here the criterion is chosen to be a good agreement in the topmost part of the Bragg peak. It is found that scaling the antiproton cross sections by a factor of 1.08 makes the simulation points fit the measured curve quite well in the top point. This 8% correction is well within the bounds of the stated accuracy from the Sychev cross sections. From Fig. 4 it can be seen that the scaling has the desired effect of restoring the agreement between the measured and the simulated data. However the small deviations in the onset of the Bragg peak and in the tail remains. No conclusions about the correctness of the tabulated antiproton cross sections can therefore be made. It can only be deduced that if the physically more correct PPD-law is to be implemented a scaling of the antiproton collision cross sections should be done as well. Since clearly the PPD antiproton capture cross section model used here is based on rather vague arguments, it would seem that

5. Conclusion A method for recalculating the measured dose-to-water into absolute values is found, and tests have been performed to justify its validity. It is concluded that it is very important to model the precise geometry of the ionization chamber, because if the depth dose curve is scored directly in water, neglecting the influence of the composition of the ionization chamber, the simulation result in the Bragg peak will be significantly different from the result obtained with the ionization chamber in place. This shows the importance of the materials; water equivalence cannot be taken as a good approximation, when calculating the dose. Incorporating the exact geometry of the ionization chamber yields a good agreement between the measured values and the simulations. The consistency is best in the plateau and in the top of the Bragg peak, whereas the onset of the Bragg peak and the tail show minor deviations. A literature study shows that the Fermi–Teller Z-law for the capture probabilities for negatively charged particles in com-

V.T. Taasti et al. / Nuclear Instruments and Methods in Physics Research B 347 (2015) 65–71

pounds does not give a proper description of the capture process. Since the Z-law is implemented into SHIELD-HIT12A, the program does not simulate the capture process correctly. When implementing a physically more correct description of the capture process found by Ponomarev, Pawlevicz and Daniel, a variation in the depth dose curve is seen, as expected from the theory, but it worsens the simulation agreement with experimental data in the top of the Bragg peak, and does not improve the small deviations in the onset of the Bragg peak and in the tail. This can be mitigated by using a scaling factor for the antiproton collision cross section, which restores the agreement seen before implementing the replacement of the Z-law, the PPD-law. However the situation is quite unsatisfactory and calls for more experimental evidence of low energy antiproton annihilation on compounds, and for more experimental data on the antiproton collision cross sections, since both the capture probabilities and the collision cross sections might be wrong. Nonetheless, a drastic improvement of the agreement between the measured and the simulated result has been demonstrated, compared to what was published in [32], stemming from an update of the antiproton cross sections, carefully implementing the ionization chamber geometry, and correcting issues in the antiproton annihilation model. We find that the lack of experimental evidence makes it impossible to enhance the antiproton simulation results further. Acknowledgments The authors thank Prof. Konstantin Gudima for consulting on programing issues of the antiproton annihilation model. PTW Freiburg is thanked for providing design details of their Advanced Roos ionization chamber. Dr. Alfredo Ferrari from the FLUKA team is thanked for fruitful discussions and collaboration. References [1] C. Maggiore, N. Agazaryan, J. DeMarco, M. Doser, T. Giorgio, C. Gruhn, et al., Relative Biological Effectiveness and Peripheral Damage of Antiproton Annihilation, CERN-SPSC-2002-028, 2002, SPSC-I-225. [2] C. Maggiore, N. Agazarayan, N. Bassler, E. Blackmore, G. Beyer, J.J. DeMarco, et al., Biological effectiveness of antiproton annihilation, Nucl. Instr. Meth. B 214 (2004) 181–185. [3] M.H. Holzscheiter, N. Agazarayan, N. Bassler, G. Beyer, J.J. DeMarco, M. Doser, et al., Biological effectiveness of antiproton annihilation, Nucl. Instr. Meth. B 221 (2004) 210–214. [4] H.V. Knudsen, M.H. Holzscheiter, N. Bassler, J. Alsner, G. Beyer, J.J. DeMarco, et al., Antiproton therapy, Nucl. Instr. Meth. B 266 (3) (2008) 530–534. [5] N. Bassler, J. Alsner, G. Beyer, J.J. DeMarco, M. Doser, D. Hajdukovic, et al., Antiproton radiotherapy, Radiother. Oncol. 86 (2008) 14–19. [6] A.H. Sullivan, A measurement of the local energy deposition by antiprotons coming to rest in tissue-like material, Phys. Med. Biol. 30 (12) (1985) 1297– 1303. [7] N. Bassler, M.H. Holzscheiter, Calculated LET spectrum from antiproton beams stopping in water, Acta Oncol. 48 (2009) 223–226. [8] M.H. Holzscheiter, N. Bassler, N. Agazaryan, G. Beyer, E. Blackmore, J.J. DeMarco, et al., The biological effectiveness of antiproton irradiation, Radiother. Oncol. 81 (2006) 233–242. [9] H. Paganetti, M. Goitein, K. Parodi, Spread-out antiproton beams deliver poor physical dose distributions for radiation therapy, Radiother. Oncol. 95 (2010) 79–86. [10] I. Kantemiris, A. Angelopoulos, N. Bassler, N. Giokaris, M.H. Holzscheiter, P. Karaiskos, et al., Real-time imaging for dose evaluation during antiproton irradiation, Phys. Med. Biol. 55 (5) (2010) N123.

71

[11] D.C. Hansen, A. Lühr, N. Sobolevsky, N. Bassler, Optimizing SHIELD-HIT for carbon ion treatment, Phys. Med. Biol. 57 (8) (2012) 2393–2409. [12] D.C. Hansen, A. Lühr, R. Herrmann, N. Sobolevsky, N. Bassler, Recent improvements in the SHIELD-HIT code, Int. J. Radiat. Biol. 88 (1–2) (2012) 195–199. [13] N. Bassler, D.C. Hansen, A. Lühr, B. Thomsen, J.B. Petersen, N. Sobolevsky, SHIELD-HIT12A – a Monte Carlo particle transport program for ion therapy research, J. Phys. Conf. Ser. 489 (1) (2014) 012004. [14] A. Lühr, D.C. Hansen, R. Teiwes, N. Sobolevsky, O. Jäkel, N. Bassler, The impact of modeling nuclear fragmentation on delivered dose and radiobiology in ion therapy, Phys. Med. Biol. 57 (16) (2012) 5169. [15] A. Lühr, M. Priegnitz, F. Fiedler, N. Sobolevsky, N. Bassler, Dependence of simulated positron emitter yields in ion beam cancer therapy on modeling nuclear fragmentation, Appl. Radiat. Isot. 83 (2014) 165–170. [16] I. Gudowska, N. Sobolevsky, P. Andreo, D. Belkic´, A. Brahme, Ion beam transport in tissue-like media using the Monte Carlo code SHIELD-HIT, Phys. Med. Biol. 49 (2004) 1933–1958. [17] K. Henkner, N. Bassler, N. Sobolevsky, O. Jäkel, Monte Carlo simulations on the water-to-air stopping power ratio for carbon ion dosimetry, Med. Phys. 36 (2009) 1230–1235. [18] M. Krämer, O. Jäkel, T. Haberer, G. Kraft, D. Schardt, U. Weber, Treatment planning for heavy-ion radiotherapy: physical beam model and dose optimization, Phys. Med. Biol. 45 (2000) 3299–3317. [19] E. Fermi, E. Teller, The capture of negative mesotrons in matter, Phys. Rev. 72 (1947) 399–408. [20] J.S. Baijal, J.A. Diaz, S.N. Kaplan, R.V. Pyle, Atomic capture of M-mesons in chemical compounds, Il Nuovo Cimento 30 (3) (1963) 711–726. [21] V.G. Zinov, A.D. Konin, A.I. Mukhin, Atomic capture of negative muons in chemical compounds, Sov. J. Nucl. Phys. 2 (5) (1966) 613–618. [22] Z.V. Krumshtein, V.I. Petrukhin, L.I. Ponomarev, Y.D. Prokoshkin, Investigation of 1T-meson capture by hydrogen in hydrogenous substances, Sov. Phys. JETP 27 (6) (1968) 906–909. [23] L.I. Ponomarev, Molecular structure effects on atomic and nuclear capture of mesons, Annu. Rev. Nucl. Sci. 23 (1973) 395–430. [24] W.T. Pawlewicz, C.T. Murphy, J.G. Fetkovich, T. Dombeck, M. Derrick, T. Wangler, Fraction of stopped antiprotons which annihilate on free hydrogen in propane, Phys. Rev. D 2 (1970) 2538–2544. [25] C.T. Murphy, G. Keyes, M. Saha, M. Tanaka, Fraction of stopped K-mesons which interact with free hydrogen in propane, Phys. Rev. D 9 (5) (1974) 1248. [26] C. Vander Velde-Wilquet, G.S. Keyes, J. Sacton, J.H. Wickens, D.H. Davis, D.N. Tovee, The probability of absorption of K-mesons at rest on free protons in propane, Lett. Nuovo Cimento (1971–1985) 5 (17) (1972) 1099–1103. [27] H. Schneuwly, V.I. Pokrovsky, L.I. Ponomarev, On coulomb capture ratios of negative mesons in chemical compounds, Nucl. Phys. A 312 (3) (1978) 419– 426. [28] H. Daniel, Formation of mesonic atoms in condensed matter, Phys. Rev. Lett. 35 (1975) 1649–1651. [29] H. Daniel, Coulomb capture of muons and atomic radius, Z. Physik A 291 (1) (1979) 29–31. [30] T. von Egidy, W. Denk, R. Bergmann, H. Daniel, F.J. Hartmann, J.J. Reidy, et al., Muonic Coulomb capture ratios and x-ray cascades in oxides, Phys. Rev. A 23 (1981) 427–440. [31] B.S. Sychev, Cross Sections of Interaction of High-Energy Hadrons with Atomic Nuclei, MRTI RAS, Moscow, 1999 (in Russian). [32] N. Bassler, M.H. Holzscheiter, O. Jäkel, S. Kovacevic, H.V. Knudsen, The AD-4/ ACE Collaboration, The antiproton depth-dose curve in water, Phys. Med. Biol. 53 (2008) 793–805. [33] IAEA TRS-398, Absorbed Dose Determination in External Beam Radiotherapy An International Code of Practice for Dosimetry Based on Standards of Absorbed Dose to Water, International Atomic Energy Agency, 2001. [34] O.M. Kartal, Cross-calibration of ionization chambers for heavy ion beams (Bachelor’s thesis), Institut für Medizinische Physik und Strahlenschutz, FH Giessen, Germany, 2007. [35] A. Lühr, J. Toftegaard, I. Kantemiris, D.C. Hansen, N. Bassler, Stopping power for particle therapy: the generic library libdEdx and clinically relevant stoppingpower ratios for light ions, Int. J. Radiat. Biol. 88 (1–2) (2012) 209–212. [36] A. Lühr, D.C. Hansen, N. Sobolevsky, H. Palmans, S. Rossomme, N. Bassler, Fluence correction factors and stopping power ratios for clinical ion beams, Acta Oncol. 50 (6) (2011) 797–805. [37] E. Widmann, Perspectives for low energy antiproton physics at FAIR, Hyperfine Interact. 229 (1–3) (2014) 123–128.