Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code

ARTICLE IN PRESS Physica Medica ■■ (2016) ■■–■■

Contents lists available at ScienceDirect

Physica Medica j o u r n a l h o m e p a g e : h t t p : / / w w w. p h y s i c a m e d i c a . c o m

Review Paper

Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code C.Q.M. Reis *, P. Nicolucci Departamento de Física, Faculdade de Filosofia Ciências e Letras de Ribeirão Preto, Universidade de São Paulo, Av. Bandeirantes 3900, 14040-901 Ribeiro Preto, SP, Brazil

A R T I C L E

I N F O

Article history: Received 7 October 2015 Received in revised form 20 January 2016 Accepted 22 January 2016 Available online Keywords: Dosimetry Correction factors Photon beams Monte Carlo simulation

A B S T R A C T

The purpose of this study was to investigate Monte Carlo-based perturbation and beam quality correction factors for ionization chambers in photon beams using a saving time strategy with PENELOPE code. Simulations for calculating absorbed doses to water using full spectra of photon beams impinging the whole water phantom and those using a phase-space file previously stored around the point of interest were performed and compared. The widely used NE2571 ionization chamber was modeled with PENELOPE using data from the literature in order to calculate absorbed doses to the air cavity of the chamber. Absorbed doses to water at reference depth were also calculated for providing the perturbation and beam quality correction factors for that chamber in high energy photon beams. Results obtained in this study show that simulations with phase-space files appropriately stored can be up to ten times shorter than using a full spectrum of photon beams in the input-file. Values of kQ and its components for the NE2571 ionization chamber showed good agreement with published values in the literature and are provided with typical statistical uncertainties of 0.2%. Comparisons to kQ values published in current dosimetry protocols such as the AAPM TG-51 and IAEA TRS-398 showed maximum percentage differences of 0.1% and 0.6% respectively. The proposed strategy presented a significant efficiency gain and can be applied for a variety of ionization chambers and clinical photon beams. © 2016 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.

water, Dw, can be related do the dose in the air cavity of the chamber, Dch , by [4–6]

Introduction According to the dosimetry protocol TRS-398 of the International Atomic Energy Agency (IAEA) [1], ionization chambers are the standard dosimeters for assessing absorbed dose to water in radiotherapy. The theoretical formalism adopted in that international document is based on the Spencer–Attix cavity theory, which allows one to know the absorbed dose to the water, Dw, from the average dose in an ideal air cavity, Dair , that fulfills the conditions of the theory. For an ionization chamber that follows those requirements the dose to the medium can be given by [2,3]

Dw = Dair ⋅ sw ,air

(1)

where sw ,air is the Spencer–Attix water to air stopping-power ratio. However, practical ion chambers in general will not fulfill the conditions of an ideal cavity. In this sense many corrections must be made when using Eq. (1) in order to make it useful for real ionization chambers. Taking into account these corrections, the dose to

* Corresponding author. Departamento de Física, Faculdade de Filosofia Ciências e Letras de Ribeirão Preto, Universidade de São Paulo, Av. Bandeirantes 3900, 14040-901 Ribeiro Preto, SP, Brazil. Tel.: +5516988211879; fax: +551633159101. E-mail address: fi[email protected] (C.Q.M. Reis).

Dw = Dch ⋅ sw ,air ⋅ pwall ⋅ pcel ⋅ pdis ⋅ pcav = Dch ⋅ sw ,air ⋅ pQ

(2)

where pQ = pwall ⋅ pcel ⋅ pdis ⋅ pcav is referred in the TRS-398 protocol as the overall perturbation factor and accounts for all different effects related to the constructive details of the chamber, which are assumed to be independent of each other [1,4,6]. The factor pwall is the correction factor which accounts for non-medium equivalence of the chamber wall and any waterproofing material. pcel is the factor that corrects the response of an ionization chamber with a central electrode material different from the cavity medium. pdis accounts for the smaller attenuation in the air cavity that replaces the water and also causes an upstream shift of the effective point of measurement. The factor pcav accounts for other changes in the electron fluence spectrum due to the fact that the cavity material (air) is not the same as the medium (water) where it is inserted. The product pdis ⋅ pcav is referred as the replacement correction factor in the TG51 protocol of the American Association of Physicists in Medicine (AAPM) [7]. On the other hand, absorbed dose determination with ionization chambers makes use of the beam quality correction factor, kQ ,Q 0 (or kQ), besides an absorbed dose to water calibration coefficient N D,w ,Q 0 in the reference quality Q0 (usually 60Co). Based on this, the

http://dx.doi.org/10.1016/j.ejmp.2016.01.482 1120-1797/© 2016 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.

Please cite this article in press as: C.Q.M. Reis, P. Nicolucci, Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code, Physica Medica (2016), doi: 10.1016/j.ejmp.2016.01.482

ARTICLE IN PRESS C.Q.M. Reis, P. Nicolucci/Physica Medica ■■ (2016) ■■–■■

2

absorbed dose to water in an arbitrary beam quality Q can be given by [1]

Dw ,Q = MQ N D,w ,Q 0 kQ ,Q 0 = MQ N D,w ,Q

(3)

where MQ is the reading of the chamber fully corrected for influence quantities other than beam quality. A theoretical approach for calculating this factor is given in dosimetry protocols according to

( sw ,air )Q [ pwall pcel pcav pdis ]Q ( sw ,air )Q pQ kQ = = ( sw ,air )Q 0[ pwall pcel pcav pdis ]Q 0 ( sw ,air )Q 0 pQ 0

(4)

where the assumption is made that the mean energy in air per ion pair, W/e, is independent of the beam quality. Comparison of Eqs. (2) and (4) shows that kQ can be given by [8,9] Q

⎛D ⎞ kQ = ⎜ w ⎟ ⎝ Dch ⎠ 60Co

(5)

which allows one to calculate kQ directly from the evaluation of the absorbed dose-to-water in the absence of the chamber, Dw, at the depth corresponding to the point of measurement, and the absorbed dose-to-air averaged over the cavity volume of the ion chamber, Dch , in the qualities Q and 60Co. Monte Carlo (MC) simulation has been considered as a “gold” standard for validating radiation dosimetry quantities as well as providing dosimetry data for the current clinical dosimetry protocols such as the AAPM TG-51 [7] or the IAEA TRS-398 [1], due to its great potential for accurately simulating the transport of radiation in the matter [10–12]. In this sense MC simulation is a useful tool for the study of correction factors used in dosimetry protocols as well as for the evaluation of the response of different ionization chambers at different radiation spectra without the costs of experimental procedures. Recent publications from the literature have shown that MC calculations have had an important role in the investigation of correction factors for ionization chambers in reference dosimetry [8,9,13–15]. However, performing simulations to obtain correction factors or comparing calculated and measured quantities requires a very long computation time. In this sense, the use of methodologies to save CPU time can improve the efficiency of simulations of ion chamber responses. Based on this, the goal of this work is to calculate beam quality and perturbation correction factors for ionization chambers in reference conditions using a strategy for reducing time simulations with PENELOPE code.

all regions except for the water phantom, where WCC = WCR = 100 keV was used. Absorption energies for electrons, positrons and photons were also set to Eabs (e − , e + ) = Eabs (γ ) = 10 keV for the regions of interest and 100 keV for the rest of the water phantom. Radiation sources Radiation sources were defined in the input-file user.in of PENELOPE by using published photon spectra found in the literature for 60Co [18] and 4, 6, 10 and 25 MV nominal energy linear accelerators [19]. In this sense, radiation sources were assumed to be point sources emitting photons with those corresponding published 20 was used as a beam quality specifier. spectra. In this study TPR10 Since the radiation sources used do not correspond to full beam Monte Carlo models of accelerators, where contaminant electrons can be taken into account, values of %dd (10)x were calculated ac20 was cording to the geometrical set-up described in TG-51. The TPR10 then calculated from the %dd (10)x by using the equation provided by Kalach and Rogers [20]: 20 TPR10 = −0.8228 + 0.0342( %dd (10)x ) − 0.0001776 ( %dd (10)x )

2

(6)

For the determination of %dd (10)x values, absorbed doses in a 30 × 30 × 30 cm3 water phantom were calculated using square scoring voxels of 0.2 cm thickness and 0.5 cm width. Phase-space file strategy for reducing time simulation A strategy based on the use of phase-space files (psf’s) is presented in order to save CPU time when calculating absorbed doses at reference depths in a water phantom in radiation therapy beams. The strategy comprises a two-step MC simulation and it is illustrated in Fig. 1. On a first step simulation, radiation beams reach the surface of a water phantom and particles are tracked until they reach an artificial volume surrounding the region of interest (ionization chamber or scoring dose voxels illustrated in Fig. 1 as a solid ellipsoid) and where their transport is immediately terminated. The

Materials and methods Simulation packages Simulations were performed using the Monte Carlo code PENELOPE 2008 version [16]. In particular, for calculating chamber correction factors, simulations were also performed by using the package clonEasy [17], which allows the use of parallel simulations in different computers. In this study clonEasy was used together with PENELOPE for performing parallel simulations on a LINUX cluster in which up to 10 CPUs were used simultaneously. Simulation parameters in PENELOPE were set in order to guarantee a reasonable compromise between speed and accuracy for all calculations. The average angular deflection parameter, C1, and the maximum average fractional energy loss parameter, C2, between two consecutive hard elastic events were set to C1 = C 2 = 0.02 in all regions of interest, i.e., dose scoring voxels regions and ion chamber geometries. For other regions of the water phantom surrounding the regions of interest, these parameters were fixed to C1 = C 2 = 0.1. The threshold energies for hard inelastic interactions, WCC , and hard bremsstrahlung emission, WCR , were set to WCC = WCR = 10 keV for

Figure 1. Illustration of the two-step Monte Carlo simulation developed in this study for absorbed dose calculations in water phantoms irradiated with photon beams. On a first step simulation (a), the radiation beam reaches the surface of the water phantom and particles are tracked until they reach an artificial volume define by a cylindrical surface (dashed line) which in turn involves the region of interest (solid ellipsoid). In a second step (b) the phase-space file (psf) stored in the artificial surface is used as source inside the phantom for all calculation within the region of interest.

Please cite this article in press as: C.Q.M. Reis, P. Nicolucci, Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code, Physica Medica (2016), doi: 10.1016/j.ejmp.2016.01.482

ARTICLE IN PRESS C.Q.M. Reis, P. Nicolucci/Physica Medica ■■ (2016) ■■–■■

surface that delimits this volume is illustrated as a cylinder in Fig. 1a (dashed line) and Fig. 1b (solid line) involving the region of interest. A phase-space file is created storing information of the particles, such as type, energy, position, direction and statistical weight, as soon as they enter the artificial volume. In a second step (Fig. 1b), the stored psf is used as a source inside the phantom for all single chamber and dose to water calculations whose geometries are involved by the artificial volume. Despite using the whole water phantom again and allowing the particles to leave and/or re-enter the artificial volume, simulation time will be reduced because the tracking of particles far from the region of interest (and that could not reach it) will not be performed again since the particles will not start their transport at the surface of the phantom. In addition, different simulation parameters can be assigned to several parts of the geometry in order to avoid using very small values in the whole phantom. The artificial volume should be chosen to accommodate different chambers and was defined by a cylinder of 5 cm diameter and 1.6 cm height centered at reference depths. Absorbed energies to water in a disk with 1 cm radius and 0.025 cm thickness centered at reference depths in a 30 × 30 × 30 cm3 water phantom were calculated using photon beams (full spectrum) with N = 3 × 109 primary histories and a 10 × 10 cm2 field at an SSD = 100 cm. Similar simulations were also performed using the stored psf in the artificial volume instead of the incident beam. Results obtained with both methods were then compared in order to test the strategy. Calculations were done at a 5 cm depth for 60Co and at a 10 cm depth for mega-voltage photon beams. Small values of simulation parameters were assigned in order to guarantee accuracy of the simulation in all regions within the cylindrical surface used to store the psf. For others regions of the water phantom (outside of the artificial volume) larger values of those parameters were used. Absorbed dose calculations Calculations of the beam quality correction factor, kQ, were accomplished by calculating the absorbed doses for each of the four terms in Eq. (5) and using the psf’s previously stored for each of the five energy photon beams and using N = 5 × 108 primary histories. Besides looking for good statistics on dose calculations, this number of primary histories was chosen in order to guarantee the accuracy of the sampled spectrum and at the same time to save memory space in the computer hard disk. All absorbed energy quantities were calculated in a 30 × 30 × 30 cm3 water phantom with the point of measurement on the central axis of the beam and at depths of 5 cm for 60Co and 10 cm for mega-voltage beams. In order to improve statistics, splitting and Russian roulette variance reduction techniques

3

were used with Nsplit = 100 and weight window WGMIN = 1 × 10−8 and WGMAX = 1 × 10−1. In this way, particles in the phase-space file which had initial weights, WGHT, less than WGMIN were subjected to Russian Roulette (with a killing probability = 1 − WGHT/ WGMIN), and those with WGHT larger than WGMAX were split into NSPLIT equivalent particles, with weights equal to WGHT/Nsplit. The values for Nsplit, WGMIN and WGMAX were optimized after several runs in order to get the best efficiency in absorbed energies. Furthermore, parallel simulations in 10 different computers were also used with the package clonEasy according to the strategy described by Badal and Sempau [17]. PENELOPE code was used to model the widely used cylindrical ionization chamber NE2571. This chamber has a volume of 0.6 cm3 and was modeled by means of quadric surfaces through the package PENGEOM. Dimensions were taken from data available in the literature [1,8,21]. Its cavity has a diameter of 0.64 cm and 2.4 cm length. The chamber has also a 2.06 cm aluminum central electrode with diameter of 1 mm. Its graphite wall has a thickness of 0.061 g/cm2. Since the NE2571 itself is not waterproof, the model also includes a 1 mm PMMA waterproofing sleeve. The chamber stem modeled here also includes portions of aluminum, graphite, and polytetrafluoroethylene (PTFE or TEFLON) [8]. Models of the NE2571 ionization chamber were built with different constructive details as illustrated in Fig. 2. A full model of the chamber was used for calculating the dose to its air cavity, Dch . The dose to the air cavity of the chamber without the central electrode, Dnoel , was calculated replacing the aluminum electrode by air. The dose to a bare air cavity of the chamber, Dair , was also calculated by considering all the other parts of the chamber made of water (without wall, sleeve and stem). Absorbed doses to water, Dw, were also calculated using a disk of water of 0.025 cm thickness and 1 cm radius as described previously. Studies in the literature have shown that no difference is observed within a statistical uncertainty ≤0.05% in the calculated dose using voxels with a thickness below 0.05 cm [8]. Perturbation factors, as defined in Eq. (2) were independently evaluated using different geometry configurations of the NE2571 ionization chamber as shown in Fig. 2 and by making the following ratios [8]:

⎛D ⎞ pQ = ⎜ w ⎟ sw ,air ⎝ Dch ⎠

(7)

⎛D ⎞ prepl = ⎜ w ⎟ sw ,air ⎝ Dair ⎠

(8)

⎛D ⎞ pcel = ⎜ noel ⎟ ⎝ Dch ⎠

(9)

Figure 2. Simulation geometries modeled with PENELOPE package for calculating correction factors. In (a) a complete geometry of the NE2571 ionization chamber (for calculating Dch ) is shown in comparison with the geometry without electrode (b) (for calculating Dnoel ) and also with the geometry corresponding only to the bare air cavity of the chamber in (c) (for calculating Dair ). A small disk of water for calculating Dw is shown in (d).

Please cite this article in press as: C.Q.M. Reis, P. Nicolucci, Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code, Physica Medica (2016), doi: 10.1016/j.ejmp.2016.01.482

ARTICLE IN PRESS C.Q.M. Reis, P. Nicolucci/Physica Medica ■■ (2016) ■■–■■

4

Table 1 20 Radiation sources and Monte Carlo calculated TPR10 quality index for the 5 energy beams used in this study. Comparisons are made between our results (third column) and those obtained by other authors that used the same energy spectra. The values between parentheses represent the relative percentage difference with our values. Our results were obtained from the calculated %dd (10)x and using Eq. (6). Nominal energy (MV)

60Co

-

Varian Clinac

4 6 10 25

Elekta SL25

⎛D ⎞ pwall = ⎜ air ⎟ ⎝ Dnoel ⎠

20 TPR10

This work

Muir and Rogers [9]

Wulff et al. [8].

Erazo and Lallena [14]

0.555 ± 1.6% 0.640 ± 1.2% 0.668 ± 1.0% 0.743 ± 0.7% 0.795 ± 0.3%

0.569 (2.5%) 0.623 (2.7%) 0.666 (0.4%) 0.734 (1.2%) 0.791 (0.5%)

0.572 (3%) 0.621 (3%) 0.662 (1%) 0.736 (0.9%) 0.797 (0.3%)

0.571 (2.8%) 0.644 (0.6%) 0.670 (0.2%) 0.739 (0.5%) 0.799 (0.5%)

(10)

Values of sw ,air needed in Eqs. (7) and (8) were also calculated in this study using the PENELOPE code. Results and discussion Beam quality index and stopping-power ratios 20 Table 1 shows the TPR10 beam quality specifier obtained in this study from calculated values of %dd (10)x and by using Eq. (6). Comparison of MC calculated values from Wulff et al. with our values shows maximum deviations of 3% for 60Co and 4 MV beams. For other energy beams deviations were less than or equal to 1%. We found a maximum deviation of 2.7% for the 4 MV beam when comparing with results from Muir and Rogers. A maximum deviation of 2.8% is found for the 60Co beam and less than 1% for the other energies when comparing our results to those of Erazo and Lallena. It is worth to point out that Eq. (6) introduces maximum deviations of 0.7% in 20 [20]. Except for 60Co where deviations were calculations of the TPR10 up to 3%, the best agreement was observed when comparing our data to those obtained by Erazo and Lallena, who used the same MC code for calculating this quality index. On the other hand, all deviations presented here contrast with the agreement found when comparing our calculated beam quality correction factors with those obtained by the same authors. Values of water-to-air ratios of the Spencer–Attix stopping power, 20 sw ,air , with Δ = 10 keV, as a function of the TPR10 are shown in Fig. 3. For comparison, data from the TRS-398 dosimetry protocol [1] were also plotted using the cubic fit provided in that report. The fit data from Kalach and Rogers [20] for TPR20,10 ≥ 0.62 is also shown in addition with a stopping-power value calculated for a 60Co beam from Rogers and Yang [22]. A cubic fit for providing sw ,air values for 20 20 0.555 ≤ TPR10 ≤ 0.795 is given to our data. For values of TPR10 near to 0.557 this fit presents differences of 0.1% and 0.03% with the values 20 of 0.8 of the TRS-398 and Rogers and Yang, respectively. For a TPR10 these differences are both approximately 0.1%. The root mean square deviation (rms) of the data from the fit is around 0.07%. Therefore, our calculated values for the water-to-air restricted stoppingpower ratios are in excellent agreement with the values provided by the TRS-398 dosimetry protocol and also with the data published by Kalach and Rogers [20] with differences less than 0.1%.

Phase-space file vs. beam simulations Table 2 compares the absorbed energy values calculated at reference depths in water using the full spectrum beam and the phasespace file as radiation sources in PENELOPE. The absorbed energies calculated using the psf’s show agreement with those calculated using the full spectrum within maximum statistical uncertainties of 0.5%. The maximum percentage relative difference was 0.3% for the 25 MV beam, which is still less than the combined uncertain-

ties of 0.4% for the absorbed energy values calculated using the two methodologies. For all energy beams investigated, the differences between the values obtained by the two methods were smaller than their respective total uncertainties. Table 2 also shows a significant difference between the simulation time using the two methods. It is important to point out that the time needed to store psf’s is on average 50% of the CPU time used to run a simulation with a full spectrum beam, except for 25 MV, where that time is approximately 80%. Therefore, for most of the photon beams used in this study 50% of the simulation time is spent out of the region of interest, which shows the importance of using phase-space files around the region of interest instead of full spectrum beams reaching the surface of the phantom. For these particular dose calculations, a direct comparison of simulation time shown in the third and fifth columns is not fair since we should take into account the time to store the psf. However, results show that once psf’s are previously stored, further calculations can be efficiently accomplished including those with several types of ionization chambers inserted into the region of interest. Chamber perturbation factors Figure 4 shows perturbation factors for the NE2571 ionization chamber in photon beams as a function of the TPR20,10 . Linear fits for all perturbation factors calculated in this study are proposed for 20 ≤ 0.795 . Fit paproviding their values within the range 0.555 ≤ TPR 10

sw,air

Beam

1.135 1.130 1.125 1.120 1.115 1.110 1.105 1.100 1.095 1.090 1.085 1.080 1.075 1.070 1.065 1.060 0.45

this work TRS-398 Kalach and Rogers fit to this work

0.50

0.55

0.60

0.65 TPR20,10

0.70

0.75

0.80

Figure 3. Spencer–Attix water to air stopping-power ratios calculated by Monte Carlo simulation. Our values calculated with PENELOPE (solid circles) are shown together with a cubic fit (straight line) given by 20 sw ,air = 1.7736 − 3.1476 (TPR10 ) + 5.2889 (TPR1020 )2 − 3.0527 (TPR1020 )3 with rms deviation of 0.07%. Data from the fit provided in TRS-398 are shown as a dashed line. A fit from 20 ≥ 0.63 plus the stopping power ratio for a 60Co beam Rogers and Kalach for TPR10 from Rogers and Yang is also shown (dotted line with triangles).

Please cite this article in press as: C.Q.M. Reis, P. Nicolucci, Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code, Physica Medica (2016), doi: 10.1016/j.ejmp.2016.01.482

ARTICLE IN PRESS C.Q.M. Reis, P. Nicolucci/Physica Medica ■■ (2016) ■■–■■

5

Table 2 Absorbed energies to water calculated in a disk with 1.0 cm radius and 0.025 cm thickness at and SSD = 100 cm using a full spectrum beam with N = 3 × 109 primary histories and a phase-space file (psf) in PENELOPE input-files. Simulations were run independently in a single computer. The last column shows the time needed to store the psf’s for each source. Beam

Absorbed energy (eV/histories) Simulations with beam

Time (days)

Simulations with psf

Time (days)

Relative difference

Time to store psf (days)

60Co

17.976 ± 0.5% 21.783 ± 0.5% 24.269 ± 0.5% 37.265 ± 0.4% 59.184 ± 0.3%

12.9 17 17 20.1 4.2

17.997 ± 0.5% 21.822 ± 0.5% 24.218 ± 0.5% 37.175 ± 0.4% 59.371 ± 0.3%

1.7 1.6 1.7 2.2 1.0

0.1% 0.2% 0.2% 0.2% 0.3%

6.3 9.1 8.4 10.5 3.3

4 MV 6 MV 10 MV 25 MV

rameters as well as root mean square deviations of the data from the fit are give in Table 3. In Fig. 4a the overall perturbation factor calculated according to Eq. (7) is compared to data obtained by Wulff et al. [8] and Erazo and Lallena [14], who performed calculations for this factor using EGSnrc and PENELOPE codes, respectively. As defined before, this factor brings together all the effects related to the constructive elements of an ionization chamber. We can see from the plot that pQ increases linearly with the increasing energy of the photon beam and it becomes closer to unity for higher energies. This means that the total perturbation effect of the chamber decreases with increas-

ing energy. This proximity to unity for pQ with increasing energy was expected considering that at higher energies the chamber approaches the Bragg–Gray condition of not disturbing the electron fluence, since in that case the range of the secondary particles (electrons) can be larger than the cavity size of the chamber. Therefore the ratio D w /D ch approaches the value of the water-to-air restricted stopping power ratio sw ,air , and pQ approaches one. The linear behavior was also observed by Wulff et al. and approximately by Erazo and Lallena. Fig. 4b shows the replacement correction factor ( pdis ∗ pcav ) also 20 . calculated for the NE2571 ion chamber as a function of the TPR 10

0.996 1.000

0.994 this work Wulff et al Erazo and Lallena fit to this work

0.992

0.998 pdis*pcav

pQ

0.990

this work Wulff et al Wang and Rogers fit to this work

0.999

0.988 0.986

0.997 0.996

0.984 0.995

0.982

0.994

0.980 0.978

0.55

0.60

0.65 0.70 TPR20,10

0.75

0.993

0.80

0.55

0.60

(a)

0.75

0.80

0.75

0.80

0.998

1.004

this work Wulff et al Buckley et al fit to this work

0.997

1.002 0.996

1.000 0.998

pcel

pwall

0.70 TPR20,10 (b)

1.006

0.996

0.995 0.994

this work (wall + sleeve + stem) Wulff et al (wall + sleeve + stem) Wulff et al (wall + sleeve) Buckley and Rogers (wall + sleeve) fit to this work

0.994 0.992 0.990 0.988

0.65

0.55

0.60

0.65 0.70 TPR20,10 (c)

0.75

0.80

0.993 0.992 0.991

0.55

0.60

0.65 0.70 TPR20,10 (d)

Figure 4. Monte Carlo calculated perturbation factors for the NE2571 ionization chamber in photon beams as a function of the TPR20,10 . Values calculated in this study (solid circles) for the overall (a), replacement (b), wall (c) and central electrode (d) perturbation factors are compared to data from the literature.

Please cite this article in press as: C.Q.M. Reis, P. Nicolucci, Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code, Physica Medica (2016), doi: 10.1016/j.ejmp.2016.01.482

ARTICLE IN PRESS C.Q.M. Reis, P. Nicolucci/Physica Medica ■■ (2016) ■■–■■

Table 3 Fit parameters to our calculated data of perturbation factors of the NE2571 ionization chamber as a function of the TPR20,10 beam quality index for a linear function of type p = a + b ∗ TPR20,10 . Perturbation factors

a

b

r2

rms % deviation

pQ pdis ∗ pcav pwall pcel

0.952 ± 0.003 0.990 ± 0.002 0.969 ± 0.004 0.989 ± 0.01

0.053 ± 0.004 0.010 ± 0.002 0.041 ± 0.006 0.007 ± 0.002

0.98 0.85 0.93 0.83

0.06 0.03 0.1 0.03

The plot presents our data calculated by using Eq. (8) in comparison with results obtained by Wulff et al. and Wang and Rogers [23]. Despite the uncertainties (0.2% in the worst case) values of pdis ∗ pcav calculated in this study are in fair agreement with those obtained by Wulff et al. considering the total combined uncertainties. For the 60 Co beam our calculated value of 0.996 ± 0.001 for this factor was closer to that obtained by Wulff et al. of 0.995 ± 0.001 than the value provided by the TRS-398 of 0.987 ± 0.003. Our result is also in fair agreement with the value of 0.9964 ± 0.08% obtained by Wang and Rogers [23]. It is worth pointing out that pdis ∗ pcav values calculated in this study follow a theoretical definition proposed by Wulff et al. [8] and also by Wang and Rogers [23] and illustrated in Fig. 2 and Eqs. (7) to (10). Values of pdis ∗ pcav available in current dosimetry protocols are based on experimental results from Cunningham and Sontag [24] (TG-51) and Johansson et al. [25] (TRS-398). According to Huq et al. [26] differences of up to 0.5% can be found for the replacement correction factors when different approaches are used. Wulff et al. also point out that taking into account uncertainties on those experimental data, deviations from Monte Carlo calculations are not surprising but still deserve further investigations. In addition, Wang and Rogers [23] also investigated differences between Monte Carlo calculations for the pdis ∗ pcav correction factor and experimental data available in dosimetry protocols. According to these authors, the interpretation of the measured data as the replacement correction in dosimetry protocols is not correct [23]. Fig. 4c presents our calculated values of pwall for the Farmer NE2571 chamber as a function of the beam quality. It is important to mention that although our results have been compared with those ones obtained by Wulff et al., those authors calculated the perturbation of the graphite wall of the chamber plus the waterproofing sleeve (pwall+sleeve) and the effect of the chamber stem pstem in a separate way. On the other hand, what we call pwall in this study is the combined effect of the perturbation caused by the graphite wall, the stem and the sleeve of the chamber (pwall+sleeve+stem). In this sense, in order to make a comparison we took the product of pwall+sleeve and pstem calculated by the authors. The results for pwall+sleeve calculated by Buckley and Rogers [27] and Wulff et al., using the EGSnrc user codes CSnrc and egs_chamber, respectively, are also shown in Fig. 4c for comparison. As expected, all curves predict a decrease in the wall perturbation effect with increasing energy which results in increasing values of pwall and closer to unity for higher energy beams. Despite the uncertainties (maximum of 0.2%) obtained for this correction factor, the values calculated for energy beams greater than 20 = 0.668 ) are in fair agreement with those or equal to 6 MV ( TPR10 obtained by Wulff et al. and Buckley and Rogers with a maximum relative difference of 0.2%. On the other hand, the values calculated for 60Co and 4 MV beams were slightly divergent to those obtained by the referred authors even if we take into account the 20 = 0.555 ) differences were effect of the chamber steam. For 60Co ( TPR10 0.68% and 0.74% compared to the values from Buckley and Rogers and Wulff et al., respectively. This last one falls to 0.43% if we take into account the steam effect calculated separately by Wulff et al. For the 6 MV beam, our calculated value for P(wall + sleeve+ steam) differs in 0.48% and 0.61% from the P(wall + sleeve) values calculated by Buckley and Rogers and Wulff et al., respectively. In comparison to

p(wall + sleeve+ steam) = p(wall + sleeve) × p(stem) from Wulff et al. the difference is around 0.28%. Although our calculated value for pwall in the 60Co beam (0.991 ± 0.2%) is close to those reported in current dosimetry protocols such as the TRS-398 (0.992 ± 0.5%) and TG-51, the Almond–Svensson formalism [1] used for calculating pwall in those protocols is known to overestimate the wall perturbation factor [28,29]. In this sense, a better investigation of these differences for lower energy beams is necessary. Figure 4d shows our central electrode perturbation factor values in comparison with those obtained by Wulff et al. and Buckley et al. [30] for the NE2571 ion chamber. As we have already discussed, this factor accounts for the central electrode material (aluminum) being different from the cavity (air). Our results suggest that pcel varies 20 goes from 0.555 to 0.80. The data from 0.993 to 0.995 when TPR10 are in good agreement, within uncertainties (maximum of 0.2%), with those calculated by Wulff with a maximum difference of 0.14% 20 of 0.555 (60Co beam). It is also in good agreement with for a TPR10 the value of 0.993 ± 0.2% provided in the TRS-398 dosimetry protocol for this chamber in the 60Co beam. The maximum difference found between our results and Buckley et al. values is 0.08% for a 20 TPR10 of 0.640 (4 MV beam). The beam quality correction factor, kQ Figure 5 shows the beam quality correction factor, kQ, for the 20 , calculated NE2571 ionization chamber as a function of the TPR10 in this study by using Eq. (5). Data provided in current dosimetry protocols [1,7] as well as results from the literature [8,9,14] are also plotted together for comparison. The first three sets of data points correspond to our results (solid circles) and data from the dosimetry protocols TRS-398 [1] (open squares) and TG-51 [7] (triangles up). Fitting curves provided by Muir and Rogers [9] (dashed line), Wulff et al. [8] (dotted line) and Erazo and Lallena [14] (point20 range. In addition, dashed line) are also plotted for the same TPR10 a cubic polynomial fit to experimental data compiled by Andreo [31] and provided by Wulff et al. [8] is also plotted (solid line) for comparison with our Monte Carlo calculated data. A comparison is also made to values measured by McEwen [32] (triangles down) for this correction factor.

1.000 0.995 0.990 0.985 0.980

kQ

6

this work TRS-398 TG-51 Muir and Rogers Wulff et al Erazo and Lallena fit to experimental data McEwen measurements

0.975 0.970 0.965 0.960 0.955 0.950 0.945 0.940

0.56

0.60

0.64

0.68 0.72 TPR20,10

0.76

0.80

0.84

Figure 5. The beam quality correction factor for the NE2571 ionization chamber. Data points are from our results (solid circles), TRS-398 (open squares) and TG-51 (triangles up). Fit to Monte Carlo calculated data are from Muir and Rogers (dashed line), Wulff et al. (dotted line) and Erazo and Lallena (point dashed line). A fit to experimental data compiled by Andreo and given by Wulff et al. is also plotted (solid line) for comparison with Monte Carlo results. Values measured by McEwen are also shown (triangles down).

Please cite this article in press as: C.Q.M. Reis, P. Nicolucci, Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code, Physica Medica (2016), doi: 10.1016/j.ejmp.2016.01.482

ARTICLE IN PRESS C.Q.M. Reis, P. Nicolucci/Physica Medica ■■ (2016) ■■–■■

7

Table 4 20 Fit parameters to a cubic polynomial of type kQ = a + b (TPR10 ) + c (TPR1020 )2 + d (TPR1020 )3 for calculating kQ values as a function of the TPR1020 . Since in this study kQ was evaluated 20 for only four (TPR10 ) values a cubic b-spline interpolation was used in order to provide a better fit to the data. Root mean square deviations of the data from the fit are given in the last column in percentage. Parameters

Fit to this work Fit to TRS-398 data Fit to TG-51 data

a

b

c

d

rms % deviation

4.615 ± 0.024 1.64 ± 0.09 1.9 ± 0.5

−15.72 ± 0.10 −3.1 ± 0.4 −4.3 ± 2.3

22.91 ± 0.14 5.1 ± 0.7 6.9 ± 3.3

−11.20 ± 0.06 −2.8 ± 0.3 −3.8 ± 1.6

0.005 0.096 0.1

Our calculated values are in good agreement (within uncertainty of 0.2%) with those obtained by the referred authors, except for 20 of 0.743 (10 MV) where difthe kQ value corresponding to a TPR10 ferences were up to 0.4% (Muir and Rogers), 0.5% (Wulff et al.), 0.7% (Erazo and Lallena) and 0.3% (Experimental data). This result points 20 values, out the necessity of calculating kQ for a larger number of TPR10 especially in the range of 10 MV beams. Agreement was also satisfactory when comparing the data provided in dosimetry protocols with maximum percentage difference of 0.6% to the TRS-398 value 20 of 0.795 (25 MV). Differences in the values provided in for a TPR10 TG-51 were 0.1% on average. For all energy beams the larger discrepancies were found in the values quoted by Erazo and Lallena. A detailed investigation on these differences could take into account 20 calculations used in the different simulation parameters and TPR10 both works. However this is far from the main purpose of this study. Cubic polynomial fits to data from TRS-398 and TG-51 were cal20 values culated in order to compare to our results for the same TPR10 used in this study. A cubic polynomial of type 20 kQ = a + b (TPR10 ) + c (TPR1020 ) + d (TPR1020 ) 2

3

(11)

20 ≤ 0.797 . Since kQ has was also fitted to our data for 0.637 ≤ TPR10 20 values, a cubic been calculated in this study for only four TPR10 b-spline interpolation was used in order to fit the data. The fit parameters of our data as well as those for the data of the protocols shown in Fig. 5 are given in Table 4. Figure 6 provides a very enlightening plot, as suggested by Rogers [33], for understanding the dependence of the beam quality correction factor, kQ, on its respective components, such as the water-

1.02

1.00 Q

k 60

Co Q

(sw,air) 60

Co

0.98

Q

(pQ) 60

Co

0

kQ,Q and its components

1.01

0.99

Q

0.97

(pdis*pcav) 60

Co

Q

(pwall) 60

Co

0.96

Q

(pcel) 60

Co

0.95

0.55

0.60

0.65

0.70 TPR20,10

to-air restricted stopping power ratio and the chamber perturbation factors, as defined in Eq. (4). All those components are taken as the ratio between their values in the quality beam, Q, and the reference quality, Q0 ( 60Co ). We can see from this plot that the main component in kQ is the water-to-air restricted stopping power ratio sw ,air , responding for a variation of approximately 4.8% in the 20 beam quality correction factor for TPR10 between 0.555 ( 60Co ) and 0.795 (25 MV). For the same energy range the combined effect of all chamber perturbation factors (pQ) responds for a variation of 1.2% in kQ. Furthermore, Fig. 5 also shows that pwall , with a variation of 20 range investigated, is the most important chamber 0.9% in the TPR10 perturbation factor and it is responsible for 75% of the 1.2% variation in pQ mentioned above. The replacement ( pdis ∗ pcav ) and central electrode ( pcel ) perturbation factors are almost equal to unity for the energy beam range investigated, both showing a variation of around 0.2%. The beam quality correction factor, defined as the ratio, in two beam qualities, of the products of all its components, is then responsible for a variation of approximately 3.6% on the NE2571 ion 20 . chamber response with the TPR10

0.75

0.80

Figure 6. Beam quality correction factor, kQ, and its components for the NE2571 ion20 ization chamber as a function of the TPR10 . The components are given by the ratio in two beam qualities of the water to air restricted stopping power ratio sw ,air and Q the overall perturbation factor ( pQ ) 60Co , as well as the factors ( pwall )Q60Co , ( prepl )Q60Co , and ( pcel )Q60Co assumed to be independent of each other.

Conclusions In this study, Monte Carlo based perturbation and beam quality correction factors for the widely used NE2571 ionization chamber have been calculated using a time saving strategy with PENELOPE code. Results have shown that the time spent on tracking particles outside of the region of interest can be more than 50% of the total simulation time. In this sense, the use of phase-space files collected around the region of interest instead of a full spectrum reaching the surface of water phantoms shows a significant efficiency. In addition, the use of variance reduction techniques which in turn can increase time in simulations becomes more feasible when using a psf as a source instead of a full spectrum beam. Calculations of the beam quality and perturbation correction factors for the NE2571 ionization chamber using phase-space files are in fair agreement with several data from the literature. Discrepancies found for the wall perturbation factor in low energies highlight the necessity for a better investigation on calculating this factor. In addition, deviations in single perturbation factors also support the direct calculation of the beam quality correction factor or a unique chamber- and quality-dependent factor as proposed by Sempau et al. [34]. Results presented in this study confirm the use of Monte Carlo simulations to provide perturbation and beam quality correction factors for ionization chambers in radiotherapy. The proposed strategy of using phase-space files in PENELOPE code presented a significant efficiency gain and can be applied for a variety of ionization chambers and energy beams used in clinics. Acknowledgments This work was partially supported by a scholarship from the Brazilian Federal Agency for the Support and Evaluation of Graduate

Please cite this article in press as: C.Q.M. Reis, P. Nicolucci, Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code, Physica Medica (2016), doi: 10.1016/j.ejmp.2016.01.482

ARTICLE IN PRESS C.Q.M. Reis, P. Nicolucci/Physica Medica ■■ (2016) ■■–■■

8

Education – CAPES (proc. number BEX 13850/12-1) and by a scholarship from the National Council for Scientific and Technological Development - CNPq (proc. number 149119/2010-0) held by C.Q.M. Reis. References [1] IAEA. Absorbed dose determination in external beam radiotherapy: an international code of practice for dosimetry based on standards of absorbed dose to water; vol. 398 of Technical Report Series. Vienna: IAEA; 2001. [2] Spencer LV, Attix FH. A theory of cavity ionization. Radiat Res 1955;3:239–54. [3] Attix FH. Introduction to radiological physics and radiation dosimetry. New York: Wiley; 1986. [4] Andreo P. Absorbed dose beam quality factors for the dosimetry of high-energy photon beams. Phys Med Biol 1992;37:2189–211. [5] Rogers DWO. The advantages of absorbed-dose calibration factors. Med Phys 1992;19:1227–39. [6] Huq MS, Andreo P. Advances in the determination of absorbed dose to water in clinical high-energy photon and electron beams using ionization chambers. Phys Med Biol 2004;49:R49. . [7] Almond PR, Biggs PJ, Coursey BM, Hanson WF, Huq MS, Nath R, et al. AAPM’s TG–51 protocol for clinical reference dosimetry of high-energy photon and electron beams. Med Phys 1999;26:1847–70. [8] Wulff J, Heverhagen JT, Zink K. Monte-Carlo-based perturbation and beam quality correction factors for thimble ionization chambers in high-energy photon beams. Phys Med Biol 2008;53:2823–36. . [9] Muir BR, Rogers DWO. Monte Carlo calculations of kQ, the beam quality conversion factor. Med Phys 2010;37:5939–50. [10] Baro J, Sempau J, Fernandez-Varea JM, Salvat F. PENELOPE: an algorithm for Monte Carlo simulation of the penetration and energy loss of electrons and positrons in matter. Nucl Inst Meth B 1995;100:31–46. [11] Rogers DWO. Fifty years of Monte Carlo simulations for medical physics. Phys Med Biol 2006;51:R287–301. [12] Kawrakow I, Mainegra-Hing E, Rogers DWO, Tessier F, Walters BRB. The EGSnrc Code System: Monte Carlo simulation of electron and photon transport. NRC Technical Report PIRS-701 v4-2-3-2. Ottawa, Canada: National Research Council Canada; 2011 . [13] Zink K, Wulff J. Beam quality corrections for parallel-plate ion chambers in electron reference dosimetry. Phys Med Biol 2012;57:1831–54. . [14] Erazo F, Lallena AM. Calculation of beam quality correction factors for various thimble ionization chambers using the Monte Carlo code PENELOPE. Phys Med 2013;29:163–70. [15] Muir BR, Rogers DWO. Monte Carlo calculations of electron beam quality conversion factors for several ion chamber types. Med Phys 2014;41:111701. (15pp).

[16] Salvat F, Fernández-Varea JM, Sempau J. PENELOPE-2008: a code system for Monte Carlo simulation of electron and photon transport. Tech. Rep. Issy-lesMoulineaux, France: OECD Nuclear Energy Agency; 2008. [17] Badal A, Sempau J. A package of Linux scripts for the parallelization of Monte Carlo simulations. Comput Phys Commun 2006;175(6):440–50. [18] Mora G, Maio A, Rogers DWO. Monte Carlo simulation of a typical 60Co therapy source. Med Phys 1999;26:2494–502. [19] Sheikh-Bagheri D, Rogers DWO. Monte Carlo calculation of nine megavoltage photon beam spectra using the BEAM code. Med Phys 2002;29:391–402. [20] Kalach NI, Rogers DWO. Which accelerator photon beams are “clinic-like” for reference dosimetry purposes? Med Phys 2003;30:1546–55. [21] Aird EGA, Farmer FT. The design of a thimble chamber for the Farmer dosemeter. Phys Med Biol 1972;17:169–74. [22] Rogers DWO, Yang CL. Corrected relationship between %dd(10)x and stoppingpower ratios. Med Phys 1999;26:538–40. [23] Wang LLW, Rogers DWO. Calculation of the replacement correction factors for ion chambers in megavoltage beams by Monte Carlo simulation. Med Phys 2008;35:1747–55. [24] Cunningham JR, Sontag MR. Displacement corrections used in absorbed dose determinations. Med Phys 1980;7:672–6. [25] Johansson KA, Mattsson LO, Lindborg L, Svensson H. Absorbed-dose determination with ionization chambers in electron and photon beams having energies between 1 and 50 MeV. In: IAEA Symposium Proceedings. IAEA, Vienna; 1978, p. 243–70. [26] Huq MS, Andreo P, Song H. Comparison of the IAEA TRS-398 and AAPM TG–51 absorbed dose to water protocols in the dosimetry of high-energy photon and electron beams. Phys Med Biol 2001;46:2985–3006. [27] Buckley LA, Rogers DWO. Wall correction factors, Pwall, for thimble ionization chambers. Med Phys 2006;33:455–64. [28] Buckley LA, Rogers DWO. Wall correction factors, Pwall, for parallel-plate ionization chambers. Med Phys 2006;33:1788–96. [29] Nahum AE. Cavity theory, stopping power ratios, correction factors. In: Rogers DWO, Cygler JE, editors. Clinical dosimetry measurements in radiotherapy. Madison (WI): Medical Physics Publishing; 2009. [30] Buckley LA, Kawrakow I, Rogers DWO. CSnrc: correlated sampling Monte Carlo calculations using EGSnrc. Med Phys 2004;31:3425–35. [31] Andreo P. A comparison between calculated and experimental kQ photon beam quality correction factors. Phys Med Biol 2000;45:L25–38. [32] McEwen MR. Measurement of ionization chamber absorbed dose kQ factors in megavoltage photon beams. Med Phys 2010;37:2179–93. doi:10.1118/ 1.3375895. . [33] Rogers DWO. Fundamentals of dosimetry based on absorbed-dose standards. In: Palta JR, Mackie TR, editors. Teletherapy physics, present and future. Washington (DC): AAPM; 1996. p. 319–56. [34] Sempau J, Andreo P, Aldana J, Mazurier J, Salvat F. Electron beam quality correction factors for plane-parallel ionization chambers: Monte Carlo calculations using the PENELOPE system. Phys Med Biol 2004;49:4427.

Please cite this article in press as: C.Q.M. Reis, P. Nicolucci, Assessment of ionization chamber correction factors in photon beams using a time saving strategy with PENELOPE code, Physica Medica (2016), doi: 10.1016/j.ejmp.2016.01.482