Édouard Guégain
email: edouard.guegain@univ-lille.fr

Clément Quinton
email: clement.quinton@univ-lille.fr

Romain Rouvoy
email: romain.rouvoy@univ-lille.fr

On Reducing the Energy Consumption of Software Product Lines

Keywords: Software Product Lines, Energy, Consumption, Measurement, Mitigation

Along the last decade, several studies considered green software design as a key development concern to improve the energy efficiency of software. Yet, few techniques address this concern for Software Product Lines (SPL). In this paper, we therefore introduce two approaches to measure and reduce the energy consumption of a SPL by analyzing a limited set of products sampled from this SPL. While the first approach relies on the analysis of individual feature consumptions, the second one takes feature interactions into account to better mitigate energy consumption of resulting products.

Our experimental results on a real-world SPL indicate that both approaches succeed to produce significant energy improvements on a large number of products, while consumption data was modeled from a small set of sampled products. Furthermore, we show that taking feature interactions into account leads to more products improved with higher energy savings per product.

CCS CONCEPTS

• Software and its engineering → Software product lines.

INTRODUCTION

Energy consumption of software systems is becoming a major concern for the environment and the society, especially with the emergence of highly-configurable and large-scale distributed systems, such as Fog, IoT, or cyber-physical systems. Along the last decade, several studies showed that software has a significant impact on the energy consumed and considered green software design as a key development concern to improve the energy efficiency of software systems at large [START_REF] Islam | Measuring energy footprint of software features[END_REF][START_REF] Jagroep | Software Energy Profiling: Comparing Releases of a Software Product[END_REF][START_REF] Pereira | SPELLing out energy leaks: Aiding developers locate energy inefficient code[END_REF]. While tools and approaches have thus been published to understand, measure, and anticipate the energy consumption of software [START_REF] Noureddine | Monitoring energy hotspots in software[END_REF], few techniques address this concern for Software Product Lines (SPL) [START_REF] Couto | Products Go Green: Worst-Case Energy Consumption in Software Product Lines[END_REF]. Measuring the energy consumption of a SPL is a challenging task due to (i) the large number of features and products that must be considered, and (ii) the time and cost required to efficiently handle such a variability. To the best of our knowledge, only Couto et al. proposed an approach to measure energy consumption in SPL [START_REF] Couto | Products Go Green: Worst-Case Energy Consumption in Software Product Lines[END_REF]. Similarly to their work, we acknowledge that measuring the energy consumption of each product of the SPL is not conceivable, in particular when considering large configuration spaces of highly-configurable software systems. But, rather than simulating and predicting the energy consumption of a single product as in [START_REF] Couto | Products Go Green: Worst-Case Energy Consumption in Software Product Lines[END_REF], we propose to measure the energy consumption of a set of products.

In this paper, we thus present a method to measure and reduce the energy consumption of multiple products at once by sampling and analyzing a minimal set of products. In particular, our method distinguishes two approaches: one that considers the energy consumption of individual features and one that takes pairs of features into account. The latter thus takes feature interactions into account when measuring the energy consumed by a product to highlight pairs of features that may cooperate or obstruct each other at a behavioral level, while altering the energy spent to complete a task. This approach also suggests candidate features whose interaction with user-required features exhibits lower energy footprint than the one produced by the initial interaction.

Our method thus provides means to estimate the energy consumption of individual product features, to highlight how feature interactions impact the energy consumed by products and to propose products with lower energy consumption while still including user-required features. We implemented our method and demonstrated it on a Java-based SPL to assess and compare the approach taking energy consumption of feature interactions into consideration with the approach only focusing on energy consumption of individual features.

In the remainder of the paper, Section 2 explains fundamentals on software energy consumption and measurement. Section 3 describes our approaches to measure and reduce the energy consumption of a product from the SPL. Section 4 presents the design and results of our experiments. Section 6 analyzes related work. Section 5 discusses the outcomes of our contributions, while Section 7 concludes the paper.

MEASURING ENERGY CONSUMPTION OF SOFTWARE

While the energy consumption of hardware components has been widely studied, software consumption only recently gained interest. By driving and managing such hardware, software is now considered as a central concern when aiming at reducing energy consumption [START_REF] Noureddine | Monitoring energy hotspots in software[END_REF]. When dealing with green concerns of software systems, the impact of such systems is often measured as power or energy consumption. While power (𝑃) measures the instantaneous consumption in Watts, energy (𝐸) reports on a accumulated consumption over a given period in Joules [START_REF] Pang | What Do Programmers Know about Software Energy Consumption[END_REF]. Throughout this paper, we will propose an approach to estimate and reduce energy consumption of SPL, expressed in Joules. Thus, all our measurements represent the total energy consumed by products from this SPL, independently of their execution time, which can vary depending on products. Energy measurement tools usually estimate the energy consumption of the CPU rather than the one of a specific software [START_REF] Noureddine | A Review of Energy Measurement Approaches[END_REF]. Thus, identifying the share of the software under study among the total energy consumption is not straightforward. To address this issue, we first sample the CPU consumption for one second before the program starts. We define this measurement as the idle energy consumption 𝑃 𝑖𝑑𝑙𝑒 , which refers to the average power consumption at rest. The general idea is to measure the energy consumed by the running software, 𝐸 𝑟𝑎𝑤 , and then substract 𝐸 𝑖𝑑𝑙𝑒 = 𝑃 𝑖𝑑𝑙𝑒 ×𝑇 𝑚𝑒𝑎𝑠𝑢𝑟𝑒 from 𝐸 𝑟𝑎𝑤 to get rid of the environment consumption.

𝐸 𝑛𝑒𝑡 = 𝐸 𝑟𝑎𝑤 -(𝑃 𝑖𝑑𝑙𝑒 × 𝑇 𝑚𝑒𝑎𝑠𝑢𝑟𝑒 ) (1)
The resulting energy consumption, 𝐸 𝑛𝑒𝑡 , can then be associated to the software under study, as depicted in Figure 1 and presented in Equation (1).

ESTIMATING AND REDUCING ENERGY CONSUMPTION FOR SPL

Unlike measuring the energy consumption of a software, measuring the energy consumption of a SPL is a non-trivial task, as multiple related software-i.e., the products of the SPL-must be measured.

These products exhibit different properties, including energy consumption, while sharing several features that perform differently in different contexts. The context of a feature can either be external to the product containing this feature-i.e., the environment hosting the product-or internal to the product. That is, a feature can exhibit different performances when combined with different sets of features. Inferring the energy consumption of a single feature by measuring it while running in one given product is therefore irrelevant and does not reflect the energy consumption of that feature in the SPL. On the other hand, measuring the energy consumption of a feature in each product individually is not feasible as some products may be complex to measure, while measuring the consumption of each product from a large SPL is not an option. To tackle these issues, we thus propose two approaches that estimate the energy consumption of features by measuring products sampled from the configuration space of the SPL, and then exploiting such sampled measures to reduce the energy consumption of any product from the SPL. Both approaches improve energy consumption of products by removing features or substituting them with other ones. However, some features are included in a product to ensure its validity with regard to the feature model, e.g., in or and xor relationships. We thus define 𝐹 𝑓 ⊂ 𝐹 as the set of features that are valid substitute features for a given feature 𝑓 . These substitute features are either sibling features of 𝑓 in or and xor relationships, or features involved in or and xor used in cross-tree constraints. On the other hand, some products may contain features due to functional constraints (e.g., stakeholder's requirements). Such features cannot be removed or substituted and are hereafter referred to as required features. Thus, from the stakeholder standpoint, all products are functionally equivalent if they contains the features required by this stakeholder.

Feature-wise Energy Analysis

Energy impact of individual features. To estimate the energy consumption of each feature from the SPL, we first measure the energy consumption of every product from the sample, using the method presented in Section 2. The energy consumption of each product is then reported in a matrix 𝑛 × 𝑚 with n the features and m the products, by copying the energy consumption of the product in the columns of each included feature. For instance, Matrix (2) defines 𝑓 1 to 𝑓 𝑛 as the available features, 𝑝 1 to 𝑝 𝑚 the sampled products, and 𝐸 𝑥 𝑦 represents the energy consumption of 𝑝 𝑥 if it includes 𝑓 𝑦 , or is left empty otherwise.

𝑓 1 𝑓 2 • • • 𝑓 𝑛 𝑝 1 𝑝 2 . . . 𝑝 𝑚           
           Ẽ (𝑓 1 ) Ẽ (𝑓 2 ) • • • Ẽ (𝑓 𝑛 ) (2)
By computing the median value of each column of the matrix, the relative energy consumption Ẽ (𝑓 ) of the feature 𝑓 represented by this column can be estimated. The expected behavior is that extreme energy consumption will cancel out and all features will have similar median energy consumption. However, if the median consumption of a feature is higher or lower than the other medians, then the presence of this feature tends to impact the performance of the products that contain it. Although such a measure does not provide a very accurate reading, it can nevertheless be used to compare the energy consumption of different features and perform preliminary optimizations, e.g., by selecting the less consuming feature among the substitutes of each feature 𝐹 𝑓 . In the remainder of the paper, we will refer to this method as the feature-wise analysis.

Feature-wise mitigation. Getting the most energy-efficient product including the required features follows a two-steps process. First, all optional features of the products are removed, thus only including the required features and the features requiring a substitution. Then, by leveraging the feature-wise analysis, the energy consumption of the remaining non-required features can be compared with their respective substitutes, to identify the one with the lowest consumption among them. Each of the non-required features are replaced by the most energy efficient substitute. This approach always converges toward an efficient product composed of no optional feature, including only the features with the lowest energy consumption within each feature substitution set.

The product resulting from this mitigation strategy is the one with the lowest energy consumption that can be obtained given an initial configuration.

Pairwise Energy Analysis

Energy impact of pairwise interactions. Although measuring the energy consumption of each feature in isolation gives a general trend, it cannot be used to compute the energy consumption of a combination of features (e.g., as the mean or the median of several individual consumption) due to the feature interactions phenomenon [START_REF] Siegmund | Predicting Performance via Automated Feature-Interaction Detection[END_REF]. Indeed, numerous work have shown that features interact with each other, hence impacting performances of products [1-4, 27, 28]. Therefore, restricting the energy consumption analysis to individual features does not provide a comprehensive landscape of a feature consumption, and additional analysis that consider feature interactions must be performed to obtain additional details about the energy consumption. By analyzing how the consumption of a feature evolves when this feature is combined with different features, it is thus possible to highlight feature interactions leading to positive or negative impact on the consumption of the product.

𝑐 1 𝑐 2 • • • 𝑐 𝑛 𝑝 1 𝑝 2 . . . 𝑝 𝑚           
           Ẽ (𝑐 1 ) Ẽ (𝑐 2 ) • • • Ẽ (𝑐 𝑛 ) (3)
A possible means to take pairwise feature interactions into account is by creating a new matrix with as many columns as there are valid pairs of features in the SPL. The consumption of each pair of features can be quantified by reporting on the energy consumption of each product in the columns of the pairs of features that this product contains, as illustrated in Matrix (3). 𝑐 1 to 𝑐 𝑛 are all the valid pairs of features, 𝑝 1 to 𝑝 𝑚 the sampled products, and 𝐸 𝑥 𝑦 is the energy consumption of 𝑝 𝑥 if it contains 𝑐 𝑦 , or is left empty otherwise. To ensure a proper interaction coverage-i.e., that all valid pairs of features are measured-the sampling of products must be performed by an algorithm ensuring such coverage.

Following the same methodology as in the feature-wise analysis, the consumption of pairwise feature interactions can be inferred by computing the median energy consumption Ẽ (𝑐) of each pair of features 𝑐. In the remainder of the paper, we will refer to this method as the pairwise analysis. It is worth noting that this method is not only valid for pairwise interactions, but can also be used to deal with larger 𝑇 -wise interactions of features.

Pairwise mitigation. Instead of replacing each feature by the substitute feature with the lowest energy consumption, this second approach iteratively picks the alternative features whose interactions with other features of the product results in a more energyefficient product. At each iteration, the approach identifies a feature to remove from the product and, if required, replaces this feature.

To identify the feature 𝑓 to be removed from a given product 𝑃, our approach relies on a scoring system: the interaction score I. The interaction score of a feature 𝑓 is computed by considering all pairs of feature containing 𝑓 in the product 𝑃, and by summing the observed median energy consumption Ẽ (𝑐) of these pairs, as described in Formula (4).1

I (𝑓 , 𝑃) = 𝑔 ∈𝑃 | 𝑔≠𝑓 Ẽ (𝑔𝑓 ) (4)
The iterations of this approach are realized as described by Algorithm 1. This algorithm iterates over the set of features until no more improvement can be performed. That is, each iteration removes or changes one feature in the product. At each iteration, the feature with the highest interaction score in the product must be removed in priority. The algorithm starts by sorting features by decreasing interaction score (line 8) and considers the first feature of the list (line 9) as a removal candidate. If this removal feature is a user-required feature and cannot be removed, it is skipped (line 13). If the removal candidate is not a required feature and must be replaced (line 17), the replacement feature is identified among all possible substitutes-i.e., 𝐹 𝑓 -by computing the interaction score of alternative features with regards to all remaining features of the product-i.e., all but the removal candidate (lines 18 to 23). The selected replacement feature is the one with the lowest interaction score among all alternative features. As a result of the iteration, a new product is created by including the replacement feature (line 26).

However, if the removal candidate has the lowest interaction score, the replacement is discarded and the algorithm skips the feature, which is kept in the product. If a feature is skipped-i.e., it was either a requirement or already the best option, a new removal candidate is defined as the next feature in the ordered list (line 34 and 11). Other features of the product will be changed over the next iterations to accommodate this skipped feature. Once a modification has been applied, the algorithm proceeds to the next iteration, unless a stop criteria is met: if a same product appears twice over different iterations, or if all features were tested during an iteration and no optimization was found (line 36).

Once a stop criteria is met, the energy consumption of the product resulting from each iteration is measured in order to monitor the energy gain. As the different mutations of the product are based on empirical data, which may be subject to imprecision and noise, it is possible that a specific iteration worsens the performance of the product. For this reason, the last step of this algorithm measures the energy consumption of the products resulting from each iteration. The product finally returned by this algorithm is the one with the lowest energy consumption, which may be the initial product in the worst case scenario (lines 38).

EMPIRICAL VALIDATION

In the previous section, we introduced two approaches to reduce the energy consumption of a given product. In this section, we experimentally assess each of these approaches. In particular, we aim to answer the following research questions:

RQ 1: Do our different analysis detect feature interactions impacting energy consumption? By applying the two approaches on the same set of products, it should be possible to determine whether feature interactions have been detected as the two analysis methods should provide different results.

RQ 2: How effective are our approaches to reduce the energy consumption of a product? The two proposed approaches rely on different analysis methods to mitigate energy consumption of products. We propose two experiments to ensure both of them improve the consumption of the products given a set of required features and evaluate how they differ.

Methodology

To assess the effectiveness of our solution when measuring energy consumption of a software product line, we performed our experiments on RobocodeSPL, a software product line designed to yield robots for Robocode [START_REF] Martinez | Software Product Line Extraction from Variability-Rich Systems: The Robocode Case Study[END_REF]. Robocode is an environment in which community-developed robots fight against each other in battles. A battle is composed of several rounds, and rounds have a time granularity of turns. During a turn, each robot taking part in the battle computes its next action and sends it to the Robocode engine which executes them all and proceeds to the next turn. A round ends when only one robot survives, and the winner of a match is the robot which caused the most damages to its opponents through the different rounds.

The RobocodeSPL proposes several implementations for the 5 mandatory features a robot requires to run properly-i.e., radar, targeting, movement, enemy selection and gun. For instance, a movement can follow linear or circular patterns, follow the walls, or ram the opponent, among others. There are also 3 optional features related to resource management (e.g., not spending more in-game energy than the robots have), for a total of 92 features and 72 leaf features. The number of valid products is 1.3 × 10 6 . Figure 2 depicts an excerpt of the feature model of RobocodeSPL.

To evaluate our approach, we launched multiple robot matches and ran our mitigation techniques to minimize the energy consumption of such matches. In particular, we launched matches opposing a sampled robot and a reference robot, the sample.Wall robot, considered as the strongest robot provided by Robocode 2 . As the goal was to minimize the energy consumption, we were not interested in which robot wins or loses the match, but in the overall energy consumption of such a match. In order to fill the pairwise analysis matrix, the sample must contain several occurrences of each valid 2 According to the Robocode Wiki: https://robowiki.net/ pair of features from the feature model. To ensure such a coverage, we relied on the T-wise algorithm ICPL [START_REF] Fagereng Johansen | An Algorithm for Generating T-Wise Covering Arrays from Large Feature Models[END_REF] with 𝑇 = 2 to sample the configuration space of RobocodeSPL. Another sampling technique may provide better uniform random samples [START_REF] Kaltenecker | The interplay of sampling and machine learning for software performance prediction[END_REF][START_REF] Munoz | Uniform random sampling product configurations of feature models that have numerical features[END_REF], but such techniques do not meet our coverage requirements. This algorithm generated 602 robots, hereafter referred to as the training sample. For each couple (sampled robot, reference robot), we ran 10 matches to consolidate the performance data, resulting in a total of 6, 020 matches of 1 round.

We used JJoules,3 a Java tool using the RAPL device of Intel CPU, to measure the energy consumption of the matches. JJoules is also able to monitor the energy consumption of the DRAM, while other tools can monitor other components, such as Hard Disk Drives. As Robocode is mainly CPU-intensive, we decided to focus on the energy consumption of the CPU. The energy consumption was monitored from the start to the end of each match, thus including the energy consumed by both robots, but excluding the energy consumption of the startup and shutdown of Robocode. All measurements were obtained from a machine running the Manjaro Linux distribution with an Intel i5 CPU at 2.9GHz and 8GB of RAM. Results, input data and instructions to reproduce these experiments are available online. 4

Results

Detecting interacting features. The pairwise analysis relies on feature interactions to mitigate energy consumption of products. The goal of this first experiment is therefore to ensure that the pairwise analysis is able to detect at least one occurrence of feature interaction. The first experiment thus compares the energy consumption of the Movement and Targeting features in different contexts-i.e., in the presence of different sets of other features. Figure 3 depicts how the energy consumption of different features evolves depending on the analysis method. 5 Figures 3a and3b report on the energy consumption (measured with the feature-wise analysis) of the products from our training sample containing respectively each targeting and movement feature. Among the targeting features, T4, T6 and T17 induce an higher energy consumption than the others, but most features show similar energy consumption, around 3 Joules. Among the movement features, M1, M20 and M26 impose the highest energy consumption, around 6 Joules, while M6, M7 and M8 report on the lowest one, slightly above 2 Joules.

These differences in energy consumption can partially be explained by the functional behavior of these features. For instance, T6 (NoTargeting) performs no particular operation and always makes the robot shoot forward-i.e., in the direction it is aiming at. This is not a smart behavior and the energy consumed by matches involving this feature depends on how fast the opponent is able to destroy this robot. By contrast, T13 (TargetAdvancingVelocitySegmentation) tries to anticipate the position of the opponent based on its speed and direction to ensure that the bullet and the opponent collide. Thus, a robot configured with T13 is able to win quickly, reducing the energy consumption despite the additional computations required to anticipate the position of the opponent.

Figure 3c presents the energy consumption of each targeting feature when the feature LinearRammingMovement is selected-i.e., M7, the most energy efficient movement feature. This figure is obtained by selecting all measurements of M7 in Figure 3b, and breaking them down per targeting feature. M7 being the best movement feature, the consumption of products containing each targeting feature is either improved or unchanged when M7 is selected. However, when targeting features are sorted by median energy consumption, their rank change depending on the context. For instance, T21 is ranked 3rd by the feature-wise analysis, but becomes 16th when M7 is selected. T13 is ranked 19th out of 23 by the feature-wise analysis, but 1st in the pairwise analysis when M7 is selected. Furthermore, the median energy consumption of the couple of M7 and T13 is 1.8J, which is lower than the medians of both of these features alone, respectively 2.2 Joules and 3.6 Joules. Therefore, despite M7 being the best movement feature, its performance can still be improved by selecting a relevant targeting feature.

As shown by the feature-wise analysis in Figure 3b, products including M9 and M12 have similar median energy consumption-i.e., respectively 3.9 Joules and 3.7 Joules. However, when paired with NoTargeting (T6), one of the worst targeting features, their consumption evolve differently, as depicted in Figure 3d. The energy consumption of products including M9 is reduced from 3.9 Joules to 3.3 Joules, while the one for products including M12 dramatically increases from 3.7 Joules to 10.6 Joules. Thus, despite being considered a sub-optimal choice by the feature-wise analysis, T6 becomes an efficient choice when paired with M9. The pair composed of T6 and M8 is another occurrence of pairwise interaction outperforming both of its members (2 Joules instead of 4.5 Joules and 2.2 Joules, respectively). This result can be explained by the behavior of the features: M8 is a ramming movement feature, meaning that it is always moving toward the opponent. In this context, the behavior of NoTargeting-i.e., always shooting forward-is very efficient, as it always hits the opponent. Such changes in the resulting energy consumption with couple of features outperforming both of their members alone show that the energy consumption of targeting and movement features changes depending on how they are paired. Therefore, it highlights feature interactions between the targeting and movement features in RobocodeSPL. This experiment thus unveiled occurrences of feature interactions allowing us to answer RQ 1 positively: the pairwise analysis is able to detect interactions significantly impacting the energy consumption of products, and such interactions were not detected by the feature-wise analysis.

Behavior without required feature. The purpose of the second experiment is (i) to ensure the two mitigation approaches lead to a product different from the initial one, and (ii) to evaluate the energy consumption improvement resulting from these approaches. The first experiment showed that the feature-wise and pairwise analysis provide different results, due to their different granularity levels. It is yet to determine if the products resulting from their respective mitigation exhibit different energy consumption.

As explained in Section 3.2, the feature-wise analysis converges toward a specific product composed of no optional feature, and the features with the lowest energy consumption in each substitution set. In RobocodeSPL, considering our optimization goal, i.e., reducing the energy consumption against sample.Wall, and without any required feature, this product is composed of the features TurnMultiplierLock, DistanceSegmentation, LinearRammingMovement, StrongestSelectEnemy, and NoFireGun. Whatever the initial product considered for improvement, the feature-wise analysis will always return this product, hereafter referred to as the Best Theoretic product, 𝐵𝑇 0 .

By applying the pairwise analysis on this product, we can determine how the pairwise analysis compares to the feature-wise analysis in the absence of required features. The result of this experiment is depicted in Figure 4, where 𝐵𝑇 0 consumes 2.8 Joules. The pairwise analysis performed three iterations before reaching a stop criteria and returning the best product of the different iterations (𝐵𝑇 2 ), whose energy consumption is 1.9 Joules, i.e., 29% lower than 𝐵𝑇 0 .

This experiment provides a partial answer to our second research question RQ 2: the pairwise analysis outperforms the best product of the feature-wise analysis by 30%, when there is no required feature.

Behavior with required features. To complete this partial answer, the third experiment is a variant of the previous experiment that takes required features into account. The purpose of this experiment is (i) to ensure the changes our approaches perform on a product containing required features effectively reduce the energy consumption of such a product, and (ii) to evaluate these reductions. The products resulting from both approaches depend on the initial product, and on which features are required in this product. Therefore, by contrast to the previous experiment, it is not possible to assess our approaches with only one initial product. Thus, we used the FeatureIDE Product Generator to produce a sample of 520 random products-i.e., 1 product tested for 2, 500 products of the SPL-hereafter referred to as the validation sample. To mimic a real use case, we defined a random feature (based on the java.util.Random class) as a requirement in each of these products.

Relying on the consumption data measured on the training sample presented in Section 4.1, we applied our two energy mitigation approaches on each product from the validation sample. We evaluated how the products resulting from both approaches perform compared to their respective initial product. The feature-wise analysis generated 520 products, and the pairwise analysis generated 2,687 products-i.e., a mean of 5 iterations per initial product. By design, the result of the first iteration of the pairwise analysis is the initial product, thus all initial products are included in these 3,207 products. We computed the performance of a product as its median energy consumption over 10 matches, for a total of 32,070 matches. Figure 5 presents the energy consumption of the products resulting from each analysis (on the vertical axis) depending on the energy consumption of the initial product (on the horizontal axis). Products that are on the improvement threshold line (𝑥 = 𝑦, identity line) performed the same as the initial product, meaning that the corresponding approach failed to reduce its consumption and returned the initial product. Products that are strictly below the improvement threshold line performed better than their respective initial products.

The feature-wise analysis improved the performance of 375 products from the validation sample (72%), while the pairwise analysis improved the performance of 501 products (96%). For 127 products (24%), the pairwise analysis found improvement when the featurewise analysis failed. For 1 product (0.2%), the feature-wise analysis found improvement while the pairwise did not. Additional analysis on this specific product tend to exclude noise or measurement error as a cause for this exception.

To get a better view on the efficiency differences between the two approaches, Figure 6 depicts their respective relative gainsi.e., by how much they reduced the energy consumption of the initial products. In the feature-wise analysis, the end of the first quartile is still at 0%, as it improved 72% of the products, whereas with the pairwise analysis the end of the first quartile is already near a 24% gain. The median gain of the feature-wise analysis is 20%. Regarding the pairwise analysis, such a gain is reached before the second quartile. Therefore, only half of the products resulting from the feature-wise analysis obtained gains higher than 20%, while the pairwise analysis improved more than three quarters of products by such a gain. Similarly, half of products improved by the pairwise analysis were improved by 40% or more, while the feature-wise analysis had such a gain for only a quarter of products. The maximum gain is similar for both approaches: 82% and 86%, respectively.

Figure 7 presents how both approaches performed on two products: those with the worst and best initial energy consumption in the validation sample. The initial product is designated with the subscript 0 (e.g., 𝑊 𝑃 0 in Figure 7a). The different iterations of the pairwise analysis on the two products are designated with their respective index (𝑊 𝑃 1 to 𝑊 𝑃 4 and 𝐵𝑃 1 to 𝐵𝑃 3 ), while the result of the feature-wise analysis used as comparison is designated with the subscript FW (𝑊 𝑃𝐹𝑊 and 𝐵𝑃𝐹𝑊 , respectively). Figure 7a depicts the product 𝑊 𝑃 0 with the worst initial energy consumption, 12 Joules. The pairwise analysis performed 4 iterations before meeting a stop criteria. Most of the gains are obtained after the first iteration, with 𝑊 𝑃 1 reducing the energy consumption by 78%. 𝑊 𝑃 2 and 𝑊 𝑃 3 brought additional gains of 32% and 3% on their preceding iteration, respectively. However, 𝑊 𝑃 4 increased the energy consumption by 8%, resulting in 𝑊 𝑃 3 being returned by the pairwise analysis, with an energy consumption 86% lower than 𝑊 𝑃 0 . The feature-wise analysis returned a product 𝑊 𝑃𝐹𝑊 with an energy consumption 82% lower than for 𝑊 𝑃 0 . The energy consumption of the product 𝑊 𝑃 3 resulting from the pairwise analysis is 21% lower than the product 𝑊 𝑃𝐹𝑊 resulting from the feature-wise analysis. Figure 7b depicts the product 𝐵𝑃 0 with the best initial energy consumption, 1.6 Joules. This product is more challenging for both of our approaches, as none of them found any optimization. The pairwise analysis performed 3 iterations with energy consumption 14%, 2% and 18% higher than 𝐵𝑃 0 , respectively. The energy consumption of the product 𝐵𝑃𝐹𝑊 resulting from the feature-wise analysis is 83% higher than 𝐵𝑃 0 . As both approach fail to find optimization, they return the initial product 𝐵𝑃 0 .

These results complete the partial answer to our second research question RQ 2: Both approaches are able to improve products, with and without required features, and the pairwise analysis outperforms the features-wise approach.

Overall, both of our approaches succeed in improving products from RobocodeSPL, with or without constraints. The feature-wise and pairwise analysis thus provided useful input data about energy consumption of features and couple of features, that could then be used to improve products through the feature-wise and pairwise mitigation processes. The pairwise analysis improved more products than the feature-wise analysis, and led to higher gains. However, although less efficient than the pairwise analysis, the feature-wise analysis is more straightforward to setup, and can be used as a first intent to reduce energy consumption. It is especially relevant in the absence of feature interactions, or in systems where pairs of features are too numerous to be exhaustively measured.

DISCUSSION

Threats to Validity. To assess our approach, we ran our experiments on a specific SPL (RobocodeSPL) to measure and reduce the energy consumption of real-world products derived from this SPL. Results, such as the success rate or the relative gains, are thus only related to this single system, and cannot be generalized. Nonetheless, our contribution can be easily applied to any SPL. The improvements resulting from applying our approaches to other SPL will depend on the initial energy consumption of products and the impact of feature interactions on these products. We leave the evaluation of our approach across a larger set of domains to a future study. A second threat to validity lies in the training sample considered. To avoid measuring all products of the SPL, we sampled the configuration space and measured 602 products, which is only 0.05% of all valid products. Such a small sample may prevent the detection of some feature interactions and therefore, energy optimization hotspots. Still, it is worth noting that despite the low number of analyzed products, significant gains were obtained on the vast majority of products using our approaches.

Limitations. During the pairwise mitigation process, products are changed over several iterations. However, it might be possible that the optimal change in a given iteration prevents further improvements in the next iterations, e.g., a sub-optimal change in that iteration might lead to further and greater improvements in the long term. Furthermore, this algorithm removes all non-required optional features, without taking into account their hypothetical positive interactions in the product. It does not either consider the possibility to add an optional feature to improve the energy consumption of the product.

In the SPL community, extensions to feature models have been developed to convey information about features [START_REF] Benavides | Automated reasoning on feature models[END_REF]. Extended feature models can be used to assign consumption data on features, in order to automatically apply optimizations. However, the adoption of such extensions may raise some challenges when dealing with consumption metrics associated to pairs of features.

In the green computing domain, a commonly-used means to reduce energy consumption of software is by refactoring inefficient The energy consumed to obtain the measurements for our experiments (i.e., the training sample of 602 products) amounts to 24, 328 Joules. In comparison, the highest energy saving among the 520 products of the validation sample is 10 Joules per match. Hence, we can consider that our approach is profitable after 2, 433 matches in the best case scenario where energy savings are high. This might seem a significant number at first, but this result must be considered keeping in mind the 1.3 × 10 6 products of the SPL that can benefit from these measurements. In addition, it cannot be generalized to others SPL since this profitability threshold tightly depends on the number of features and products of the analyzed SPL.

Finally, the pairwise analysis method relies on a sample of products containing all pairs of features. As the number of features in the SPL grows, the number of pairs had a quadratic growth. For larger feature models, the use of heuristics to identify interactions between pairs of feature may proves necessary.

RELATED WORK

The approaches presented in this paper lie at the intersection of software product lines and so-called green computing. This section discusses related work belonging to both of these research fields.

Green software. Green computing approaches applied to software are mainly focused on the evaluation of energy consumption [START_REF] Noureddine | A Review of Energy Measurement Approaches[END_REF]. Rather than evaluating software artefacts, such as functions or classes, Islam et al. [START_REF] Islam | Measuring energy footprint of software features[END_REF] evaluates the energy consumption of functionalities by slicing the source code to isolate such functionalities. Then, they measure the energy consumption of these functionalities by executing and measuring the consumption of the related sliced code. While this approach works well to measure energy consumption of features in isolation, it does not take feature interaction into consideration, as proposed by our contribution.

Other studies focus on the identification of inefficient code ("energy hotspots" or anti-patterns) [START_REF] Noureddine | Monitoring energy hotspots in software[END_REF][START_REF] Pereira | SPELLing out energy leaks: Aiding developers locate energy inefficient code[END_REF]. For instance, Pereira et al. [START_REF] Pereira | SPELLing out energy leaks: Aiding developers locate energy inefficient code[END_REF] report on an approach to measure the energy consumption of a software facing different input workloads. They propose to model the consumption of methods by analyzing the consumption of the software and the number of calls for each method of this software. Once such hotspots are identified, energy consumption can be mitigated through refactoring operations, such as changes in the implemented data structures or algorithms [START_REF] Cruz | Using Automatic Refactoring to Improve Energy Efficiency of Android Apps[END_REF][START_REF] Cruz | Leafactor: Improving Energy Efficiency of Android Apps via Automatic Refactoring[END_REF][START_REF] Ournani | Tales from the Code #1: The Effective Impact of Code Refactorings on Software Energy Consumption[END_REF]. These approaches have not been designed to take variability (especially code shared among features) into account, but they can be considered as complementary to our approaches. They provide a fine-grained analysis of the energy consumption of the software at code level, while our measures rely on the behavior (energy-wise) of the features.

Green SPL. Recently, several approaches have been proposed to deal with the energy consumption of highly configurable software systems. In particular, this concern has been addressed relying on dynamic SPL to reconfigure the system depending on context changes and ensure it continues meeting its green requirements [START_REF] Dubslaff | Probabilistic Model Checking for Energy Analysis in Software Product Lines[END_REF][START_REF] Horcas | Context-aware energy-efficient applications for cyber-physical systems[END_REF][START_REF] Munoz | Energy-aware environments for the development of green applications for cyber-physical systems[END_REF]. These approaches also take feature interactions into account, but they rely on an exhaustive detection of such interactions. This detection can be done by tools implementing dedicated heuristics, or by domain experts. Thus, this detection may be error prone. In this paper, the pairwise mitigation process assumes that each feature interacts with all other features, and is thus unaffected by detection errors.

Couto et al. also proposed different techniques to evaluate energy consumption in software product lines using static analysis [START_REF] Couto | Products Go Green: Worst-Case Energy Consumption in Software Product Lines[END_REF][START_REF] Couto | Statically Analyzing the Energy Efficiency of Software Product Lines[END_REF]. These techniques estimate the energy consumption of features in the worst case scenario by analyzing the source code of features to deduce energy consumption of products, whereas our approaches measure the median energy consumption of running products. Contrarily to our approaches, they did not aim at suggesting improvements to products based on their estimations, nor took feature interactions into account.

Performances. Energy consumption can be generalized as a performance indicator. Considering the general problem of optimizing product performances, numerous work has been done to model and predict such performances. Siegmund et al. provided various contributions related to performance models and performance predictions in software product lines [START_REF] Siegmund | Performance-Influence Models for Highly Configurable Systems[END_REF][START_REF] Siegmund | Scalable prediction of non-functional properties in software product lines: Footprint and memory consumption[END_REF]. These approaches take feature interactions into account via a systematic identification. Statistical analysis of a sample of products has also already been used by Guo et al. to predict the performance of a product based on its configuration [START_REF] Guo | Variability-Aware Performance Prediction: A Statistical Learning Approach[END_REF]. In this prediction approach, feature interactions are detected using the systematic approach presented in [START_REF] Siegmund | SPL Conqueror: Toward optimization of non-functional properties in software product lines[END_REF]. Such works are designed to predict performances, but do not suggest optimizations for poorly-performing products.

Different authors provide multi-objective optimization frameworks for configurations [START_REF] Hierons | SIP: Optimal Product Selection from Feature Models Using Many-Objective Evolutionary Optimization[END_REF][START_REF] Olaechea | Modelling and Multi-Objective Optimization of Quality Attributes in Variability-Rich Software[END_REF][START_REF] Siegmund | Predicting Performance via Automated Feature-Interaction Detection[END_REF]. Such frameworks are designed to optimize multiple performance indicators, and energy consumption can be one of them. However, they also rely on a systematic identification of feature interactions. Soltani et al. [START_REF] Soltani | Automated Planning for Feature Model Configuration Based on Functional and Non-Functional Requirements[END_REF] propose an approach relying on artificial intelligence to configure products meeting stakeholders' functional requirements, preferences and performance goals. This approach is complementary to ours, as we do not take stakeholders preferences into account, and this approach does not take energy consumption or feature interactions into account. By contrast to these works, our approaches do not require a systematic identification of feature interactions. We assume that each feature interacts with all other features.

CONCLUSION

In this paper, we proposed a method to measure and reduce energy consumption in software product lines. Our contribution is twofold. First, we showed that it is possible to estimate the energy consumption of a single feature by measuring the consumption of a small set of products containing this feature. Second, we provided a means to identify energy consumption of couples of features to take feature interactions into account without detecting them in a systematic way. We applied our approach on RobocodeSPL and improved the energy consumption of 96% of randomly sampled products. In particular, half of these products have seen their energy consumption reduced by at least 40%.

As future work, we plan to improve the pairwise mitigation to tackle its current limitations. We also plan to use our approaches on widely used variable systems, such as the Linux kernel or Apache Web Server.

Figure 1 :

 1 Figure 1: Raw vs. net energy consumption.

Figure 2 :

 2 Figure 2: Excerpt of the Feature Model of RobocodeSPL.

  (a) Energy consumption per targeting feature. (b) Energy consumption per movement feature.(c) Energy consumption per targeting features when feature LinearRam-mingMovement (M7) is selected.(d) Energy consumption per movement features when feature NoTargeting (T6) is selected.

Figure 3 :

 3 Figure 3: Energy variations between the movement and targeting features and their potential interactions.

Figure 4 :

 4 Figure 4: Improving the product resulting from the featurewise analysis with the pairwise one.

Figure 5 :

 5 Figure 5: Energy consumption of the products resulting from both analysis.

Figure 6 :

 6 Figure 6: Relative gains of the pairwise and feature-wise analysis.

  (a) Worst initial product. (b) Best initial product.

Figure 7 :

 7 Figure 7: Focus on the best and worst initial products from the validation sample.

Algorithm 1 :

 1 Interaction mitigation1 Input initialProduct : a product to improve 2 Output bestProduct : the best product from all iterations 3 iterations.addProd(initialProduct)

	4 𝑖𝑚𝑝𝑟𝑜𝑣𝑎𝑏𝑙𝑒 ← 𝑡𝑟𝑢𝑒
	5 while improvable do
	6	𝑐𝑢𝑟𝑟𝑒𝑛𝑡𝑃𝑟𝑜𝑑 ← 𝑖𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛𝑠.𝑙𝑎𝑠𝑡𝐼𝑡𝑒𝑚()
	7	currentProd.removeNonRequiredOptionalFeatures()
	8	𝑠𝑜𝑟𝑡𝑒𝑑𝐹𝑒𝑎𝑡 ← 𝑠𝑜𝑟𝑡𝐵𝑦𝐼𝑛𝑡𝑒𝑟𝑎𝑐𝑡𝑖𝑜𝑛𝑆𝑐𝑜𝑟𝑒 (𝑐𝑢𝑟𝑟𝑒𝑛𝑡𝑃𝑟𝑜𝑑)
	9	𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝐼𝑛𝑑𝑒𝑥 ← 0
	10	𝑛𝑜𝐶ℎ𝑎𝑛𝑔𝑒𝐹𝑜𝑢𝑛𝑑 ← 𝑡𝑟𝑢𝑒
	11	while (𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝐼𝑛𝑑𝑒𝑥 <
		𝑐𝑢𝑟𝑟𝑒𝑛𝑡𝑃𝑟𝑜𝑑.𝑠𝑖𝑧𝑒) ∧ 𝑛𝑜𝐶ℎ𝑎𝑛𝑔𝑒𝐹𝑜𝑢𝑛𝑑 do
	12	𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒 ← 𝑠𝑜𝑟𝑡𝑒𝑑𝐹𝑒𝑎𝑡 .𝑔𝑒𝑡 (𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝐼𝑛𝑑𝑒𝑥)
	13	if ¬𝑖𝑠𝑅𝑒𝑞𝑢𝑖𝑟𝑒𝑚𝑒𝑛𝑡 (𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒) then
	22	end
	23	end
	24	if 𝑐𝑢𝑟𝑟𝑒𝑛𝑡𝐵𝑒𝑠𝑡 ≠ 𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒 then
	25	prodCandidate.addFeat(currentBest)
	26	iterations.addProd(prodCandidate)
	27	𝑛𝑜𝐶ℎ𝑎𝑛𝑔𝑒𝐹𝑜𝑢𝑛𝑑 ← 𝑓 𝑎𝑙𝑠𝑒
	28	end
	29	else
	30	iterations.addProd(prodCandidate)
	31	𝑛𝑜𝐶ℎ𝑎𝑛𝑔𝑒𝐹𝑜𝑢𝑛𝑑 ← 𝑓 𝑎𝑙𝑠𝑒
	32	end
	33	end
	34	remCandIndex++
	35	end
	36	𝑖𝑚𝑝𝑟𝑜𝑣𝑎𝑏𝑙𝑒 ← (𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝐼𝑛𝑑𝑒𝑥 <
		𝑐𝑢𝑟𝑟𝑒𝑛𝑡𝑃𝑟𝑜𝑑.𝑠𝑖𝑧𝑒) ∧ 𝑎𝑙𝑙𝐷𝑖 𝑓 𝑓 𝑒𝑟𝑒𝑛𝑡 (𝑖𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛𝑠)
	37 end
	38 return lowestEc(iterations)

14 𝑝𝑟𝑜𝑑𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒 ← 𝑐𝑜𝑝𝑦 (𝑐𝑢𝑟𝑟𝑒𝑛𝑡𝑃𝑟𝑜𝑑) 15 𝑝𝑟𝑜𝑑𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒.𝑟𝑒𝑚𝑜𝑣𝑒 (𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒) 16 𝑠𝑢𝑏𝑂𝑝𝑡𝑖𝑜𝑛𝑠 ← 𝑎𝑙𝑙𝐹 𝑠𝑢𝑏 (𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒) 17 if subOptions then 18 𝑐𝑢𝑟𝑟𝑒𝑛𝑡𝐵𝑒𝑠𝑡 ← 𝑟𝑒𝑚𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒 19 for 𝑠𝑢𝑏𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒 ∈ 𝑠𝑢𝑏𝑂𝑝𝑡𝑖𝑜𝑛𝑠 do 20 if I(subCandidate, prodCandidate) < I(currentBest, prodCandidate) then 21 𝑐𝑢𝑟𝑟𝑒𝑛𝑡𝐵𝑒𝑠𝑡 ← 𝑠𝑢𝑏𝐶𝑎𝑛𝑑𝑖𝑑𝑎𝑡𝑒

ACM acknowledges that this contribution was authored or co-authored by an employee, contractor or affiliate of a national government. As such, the Government retains a nonexclusive, royalty-free right to publish or reproduce this article, or to allow others to do so, for Government purposes only. SPLC'21, September 6-11, 2021, Leicester, United Kingdom © 2021 Association for Computing Machinery. ACM ISBN 978-1-4503-8469-8/21/09. . . $15.00 https://doi.org/10.1145/3461001.3471142

The interaction score can also be used during the configuration process, e.g., to assist the user when selecting the most energy efficient features when dealing with a partial configuration.

https://github.com/powerapi-ng/j-joules

https://doi.org/10.5281/zenodo.5048316

Mapping to real feature names available in the open data

ACKNOWLEDGMENTS

This work was partially funded by the ANR-19-CE25-0003-01 KOALA project. We thank Baptiste Lewandoski for his help during the development of the tools supporting our experiments.


Camille Adrien

Cindy Le Loirec
email: cindy.leloirec@cea.fr

Serge Dreuil

Jean-Marc Bordy

A new Monte Carlo tool for organ dose estimation in Computed Tomography

Keywords: Monte Carlo Simulation, Computed Tomography, Organ dose

The constant increase of computed tomography (CT) exams and their major contribution to the collective dose led to international concerns regarding patient dose in CT imaging. Efforts were made to manage radiation dose in CT, mostly with the use of the CT dose index (CTDI). However CTDI does not give access to organ dose information, while Monte Carlo (MC) simulation can provide it if detailed information of the patient anatomy and the source are available.

In this work the X-ray source and the geometry of the GE VCT Lightspeed 64 were modelled, based both on the manufacturer technical note and some experimental data. Simulated dose values were compared with measurements performed in homogeneous conditions with a pencil chamber and then in CIRS ATOM anthropomorphic phantom using both optically stimulated luminescence dosimeters (OSLD) for point doses and XR-QA Gafchromic® films for relative dose maps. Organ doses were ultimately estimated in the ICRP 110 numerical female phantom and compared to data reported in the literature.

Comparison of measured and simulated values show that our tool can be used for a patient specific and organ dose oriented radiation protection tool in CT medical imaging.

Introduction

Due to the significant rise of computed tomography (CT) exams in the past few years and the resulting increase of the collective dose [START_REF] Kalender | Dose in x-ray computed tomography[END_REF], patient dose in CT imaging has become a subject of interest in public health, especially for children [START_REF] Akhlaghi | Dose estimation in reference and non-reference pediatric patients undergoing computed tomography examinations: a Monte Carlo study[END_REF], Journy et al 2017a, Journy et al 2017b[START_REF] Habib Geryes | Patient dose evaluation in computed tomography: A French national study based on clinical indications[END_REF]. Considerable efforts have thus been made these past few years to manage radiation dose in CT [START_REF] Coakley | CT Radiation Dose: What can you do right now in your practice?[END_REF][START_REF] Amis | CT Radiation Dose; Trending in the Right Direction[END_REF][START_REF] Sodickson | Strategies for reducing radiation exposure in multi-detector row CT[END_REF] and Mayo-Smith et al 2014). However CT protocols are still traditionally optimized using the CT dose index (CTDI), which is not representative of the patient dose [START_REF] Kalender | Dose in x-ray computed tomography[END_REF][START_REF] Mccollough | CT Dose Index and Patient Dose: They are not the same thing[END_REF]. Patient-specific dose reports, including absorbed dose to organs, should thus be suitable for individualized protocol optimization. Because absorbed dose to organ cannot be directly measured, some research groups attempt to adapt the treatment planning system (TPS) used in radiation therapy for kV x-ray beams dosimetry. For example, Alaei et al 2000 investigated the accuracy of a convolution/superposition TPS for predicting kV beam dosimetry, but they reported dose discrepancies up to 145% in the region surrounding bone heterogeneities. Others groups developed specific software based on precomputed Monte Carlo (MC) data, such as CT imaging1 [START_REF] Kalender | A PC program for estimating organ dose and effective dose values in computed tomography[END_REF] and ImPact2 . Axial or helical protocols are included in these software, but large discrepancies can occur for helical acquisitions, since doses are approximated from contiguous axial scan. Instead of approximating dose from precomputed MC data, other research groups developed their own MC software to directly estimate organ dose using computational patient models or patient DICOM images as inputs.

Most of the time, code benchmarking was performed using either a CTDI phantom (Jarry et al 2003[START_REF] Marco | A Monte Carlo based method to estimate radiation dose from multidetector CT (MDCT): Cylindrical and anthropomorphic phantoms[END_REF][START_REF] Deak | Validation of a Monte Carlo tool for patient specific dose simulations in multi-slice computed tomography[END_REF] or a custom-designed cylindrical 3 phantom which enables dose measurements at seven different radial distances from its central axis [START_REF] Li | Patient-specific radiation dose and cancer risk estimation in CT: Part I. Development and validation of a Monte Carlo program[END_REF]. That approach has some limitations because standard and custom-made phantoms are made of PMMA and only permit dose measurements at selected distances from the central axis. Impact of heterogeneities and 2D dose gradients cannot be estimated with such phantoms. Due to the limitations of validation in homogeneous conditions, benchmark in anthropomorphic conditions is mandatory. For that purpose, DeMarco et al 2005 put MOSFET detectors on the surface of a male anthropomorphic phantom, while [START_REF] Deak | Validation of a Monte Carlo tool for patient specific dose simulations in multi-slice computed tomography[END_REF] In this work we present the dosimetric validation of a MC tool based on PENELOPE [START_REF] Sempau | Experimental benchmarks of the Monte Carlo code PENELOPE[END_REF] by comparing simulated and experimental dose estimations. Without any detailed information provided by the manufacturer, the CT scanner was first modelled only using information provided by the manufacturer technical note and the method proposed by [START_REF] Turner | A method to generate equivalent energy spectra and filtration models based on measurement for multidetector CT Monte Carlo dosimetry simulations[END_REF] In a first step the model has been validated in homogeneous conditions by comparing experimental and simulated integrated dose obtained with a pencil chamber. Simulated and measured CTDI values have also been compared. Then point dose measurements in an anthropomorphic phantom, using optically stimulated luminescence detectors (OSLD), were compared with simulated dose values. Dose distributions in the phantom were also measured with Gafchromic XR-QA2 films and compared with the simulated dose distributions to validate the calculated dose gradient in anthropomorphic conditions.

Although anthropomorphic validations are mandatory, they require a detector with weak energy dependence and an accurate calibration for the energy spectra involved in CT. Such conditions are not met nowadays with ionization chamber used in CT. Calibrations are usually performed for one normalized spectrum with uncertainties higher than 3-5% and without any correction factor for other beam qualities. Therefore, given the high number of parameters influencing the measurements and their high uncertainties (about 10% for device parameters such as mAs), 20%differences between simulated and measured results are generally considered as a good approximation of the real situation and were chosen in the present study as the success criterion.

Finally, as an application case of the tool, we estimated organ doses with a numerical anthropomorphic phantom.

2 Monte Carlo simulation

CT scanner

The VCT Lightspeed 64 CT scanner (GE Healthcare, Waukesha, WI) is a third-generation scanner. As specified in the technical note, focal spot to isocentre distance is 541 mm. Target is a tungsten-rhenium alloy with a 7 degree tilt angle regarding to reference axis. Dual focal spots are used depending on current and voltage values. According to IEC 60336, small focal spot is 0.7 mm x 0.6 mm and large focal spot is 0.9 mm x 0.9 mm. The beam full width at half maximum (FWHM) is adapted as a function of the focal spot size and the chosen aperture. For a 40 mm beam aperture, the FWHM are 42.6 and 42.9 mm for the small and the large focal spots, respectively. Use of bowtie filter is related to a maximum Scan Field of View (SFOV). The maximum SFOV is 32 cm for the small bowtie filter reported as "Ped Body" and 50 cm for the large bowtie filter reported as "Large Body". Two tube potentials (100 and 120 kVp), two bowtie filters ("Ped Body" and "Large Body") and a 40 mm beam collimation are modelled and presented in this work.

CT model

The 2006 release of the MC code PENELOPE [START_REF] Sempau | Experimental benchmarks of the Monte Carlo code PENELOPE[END_REF] is used to model the VCT Lightspeed 64 CT scanner (GE Healthcare, Waukesha, WI).

For MC simulations, all elements of the X-ray tube need to be modelled. Without any detailed information provided by the manufacturer, filtration elements are modelled by adapting the experimental method described by [START_REF] Turner | A method to generate equivalent energy spectra and filtration models based on measurement for multidetector CT Monte Carlo dosimetry simulations[END_REF] Based on an initial soft spectrum and experimental data, this method provides information on equivalent inherent filtration and equivalent bowtie filter shapes, which reproduce the attenuation of the real filtration elements.

These filtration elements are reported in table 1. For the bowtie filter only the central filtration is indicated but the complete equivalent shape is taken into account in the simulation. One of the specificities of CT scanner is the tube motion of the gantry during the acquisition. Both scanning modes, axial and helical, associated with a specific tube path, has to be carefully implemented into the MC tool. For that purpose a specific source is defined. Instead of using a direct sampling of particles along the tube path and a limited particle splitting [START_REF] Li | Patient-specific radiation dose and cancer risk estimation in CT: Part I. Development and validation of a Monte Carlo program[END_REF], it was taken advantage of the symmetries in geometry to implement a pipe shaped source with two variance reduction techniques: a circular splitting and a translational one.

The different elements of the X-ray tube are modeled with PENELOPE according to the manufacturer technical note 3 and to Turner method [START_REF] Turner | A method to generate equivalent energy spectra and filtration models based on measurement for multidetector CT Monte Carlo dosimetry simulations[END_REF] for the filtration. A PSF is created below the bowtie filter at 15 cm from the focal spot (figure 1). This PSF contains all the 3 LightSpeed TM VCT -Technical Reference Manual 6 relevant information required to perform the simulation: position, direction, energy, statistical weight, particle type. Each stored particle is read, split and released for simulation. For the circular splitting, the initial particle is splitted in NC particles. Each splitted particle is sampled on an arc of circle by a random angle (φ) in the interval [2𝜋 × 𝑖 𝑁 𝐶 ⁄ : 2𝜋 × (𝑖 + 1) 𝑁 𝐶 ⁄ ]

with i between 0 and (NC-1). User can define NC knowing that a large NC improves the statistics but increases the computation time and might introduce bias in the simulation results. For both scanning modes X and Y particle coordinates are modified, with X and Y axes the transverse and vertical directions, respectively. In axial mode simulation, the Z coordinate remains unchanged.

On the contrary, for helical acquisition, the Z coordinate is modified according to equation 1:

𝑧 ′ = 𝑧 + 𝜑 2𝜋 ⁄ × 𝑝𝑖𝑡𝑐ℎ × 𝑐𝑜𝑙𝑙𝑖𝑚𝑎𝑡𝑖𝑜𝑛 (𝑒𝑞𝑢. 1)
Once the circular splitting is realized, each particle is splitted a second time in NR particles for the translational splitting. The NR value is determined according to the acquisition parameters as:

𝑁 𝑅 = 𝑡𝑜𝑡𝑎𝑙 𝑒𝑥𝑝𝑙𝑜𝑟𝑒𝑑 𝑙𝑒𝑛𝑔𝑡ℎ 𝑝𝑖𝑡𝑐ℎ × 𝑐𝑜𝑙𝑙𝑖𝑚𝑎𝑡𝑖𝑜𝑛 (𝑒𝑞𝑢. 2)
so that it allows a complete covering of the scanning range along the Z axis. The new Z coordinate is then determined according to equation 3, placing each NR particle on a gantry rotation.

𝑧 ′′ = 𝑧 ′ + 𝑖 × 𝑝𝑖𝑡𝑐ℎ × 𝑐𝑜𝑙𝑙𝑖𝑚𝑎𝑡𝑖𝑜𝑛 𝑤𝑖𝑡ℎ 𝑖 ∈ [0, 𝑁 𝑅 ] (𝑒𝑞𝑢. 3)

Due to the lack of information regarding the tube starting position, which is known to largely influence dose results [START_REF] Zhang | Reducing radiation dose to selected organs by selecting the tube start angle in MDCT helical scans: A Monte Carlo based study[END_REF], a unique tube path cannot be determined. Choosing only one tube path can undeniably lead to a dose underestimation/overestimation for some peripheral organs if the real tube motion is largely different of the simulated path. Because a dose underestimation cannot be considered as this tool is designed for radiation protection purposes, choice is made to simulate every possible path. In such a way dose is systematically overestimated because all peripheral organs will "see" the tube along its path. Tube starting angle is thus randomly chosen for each initial particle.

Analysis of MC data

The tube loading information is used to normalize MC results. The tube loading can be related to a number of emitted electrons (Ne) according to equation 4, where I represents the CT scanning current expressed in Ampere, t the acquisition time expressed in second, and e the electron elementary charge expressed in Coulomb.

𝑁 𝑒 = 𝐼(𝐴) × 𝑡(𝑠) 𝑒 (𝑒𝑞𝑢. 4)
Emitted primary electrons are actually related to the simulated primary showers in the MC simulation.

Deposited energy in a volume

To validate the MC model in homogeneous conditions a RadCal 10X6-3CT pencil chamber (RTI electronics) is used and modeled as a 10 cm length cylinder. The cylinder radius is the one reported in the chamber documentation. The parameters used for the simulations are reported in table 2. Elastic-scattering parameters C1 and C2 are set to 0.05 to have simulations with the highest precision. Values of the cut-off energies WCC and WCR are fixed at 100 eV. Absorption energies of 100 keV for electrons and positrons is compatible with the volume of detection as the range of a 100 keV electron is about 0.15 mm in water. To simulate acquisitions realized with this pencil chamber, integrated dose in the pencil beam is estimated as :

𝐷𝐼 = 𝐸 × 𝐼 × 𝑡 × 𝐿 𝜌 𝑎𝑖𝑟 × 𝑉 × 10 -3 (𝑒𝑞𝑢. 5)
with DI the integrated dose (in Gy.cm), E the mean deposited energy (in eV/shower) obtained with the MC simulation in the volume represented the pencil chamber, I the scanning current (in A), t the acquisition time (in s), L the chamber length (in cm), 𝜌 𝑎𝑖𝑟 the air density (in g/cm 3 ) and V the chamber volume (in cm 3 ).

Uncertainty on the estimated DI is evaluated combining the contributions to the uncertainty budget using the so called "sandwich law" described in the guide to the expression of uncertainty in measurement (ISO/IEC Guide 98-3:2008):

𝑢 𝑐 2 (𝐷𝐼) = ( 𝐼 × 𝑡 × 𝐿 𝜌 𝑎𝑖𝑟 × 𝑉 × 10 -3 ) 2 × 𝑢 2 (𝐸) + ( 𝐸 × 𝑡 × 𝐿 𝜌 𝑎𝑖𝑟 × 𝑉 × 10 -3 ) 2 × 𝑢 2 (𝐼) + ( 𝐸 × 𝐼 × 𝐿 𝜌 𝑎𝑖𝑟 × 𝑉 × 10 -3 ) 2 × 𝑢 2 (𝑡) + ( 𝐸 × 𝐼 × 𝑡 𝜌 𝑎𝑖𝑟 × 𝑉 × 10 -3 ) 2 × 𝑢 2 (𝐿) + ( 𝐸 × 𝐼 × 𝑡 × 𝐿 𝜌 𝑎𝑖𝑟 × 𝑉 2 × 10 -3 ) 2 × 𝑢 2 (𝑉) (𝑒𝑞𝑢. 6)
The values of the uncertainties u(E), u(I), u(t), u(L) and u(V) are detailled in paragraph 2.3.3.

Mean dose deposited in a volume

For each CT acquisition with the anthropomorphic phantom, the DICOM images of the CIRS phantom are converted by associating to each Hounsfield Number a MC material. This conversion is done by using a calibration function determined with the CIRS Electron density phantom. The calibration function is reported in figure 2. Only four materials are considered here: air, bone, soft tissue and lung. For point doses comparison, the values obtained in the voxels corresponding to inserts containing OSL are compared to measurements.

For dose maps comparison, the values obtained in the pixels of the slices corresponding to the film position are compared to measurements. To do that, dose maps from film read-outs and from MC simulations are centered. As pixel size is not the same between experimental dose map (0.51 mm) and MC simulations (2.9 mm) normalized dose profiles are plotted to compare the results.

To improve the MC computation time for X-ray dosimetric purposes, secondary electrons are not tracked if their range is smaller than the considered voxel size [START_REF] Deak | Validation of a Monte Carlo tool for patient specific dose simulations in multi-slice computed tomography[END_REF][START_REF] Li | Patient-specific radiation dose and cancer risk estimation in CT: Part I. Development and validation of a Monte Carlo program[END_REF].

The effects of this approximation have already been investigated in detail by Chao et al 2001, who showed negligible differences between incorporating and omitting secondary electrons transport for diagnostic energy beams. The parameters used for the simulation respect this assumption.

They are reported in table 3 for each material. To validate the MC model in the anthropomorphic phantom, the mean dose in a voxel has to be estimated as:

𝐷 = 𝐷 𝑀𝐶 × 𝐼 × 𝑡 × 1000 (𝑒𝑞𝑢. 7)
with D the dose in the voxel (in Gy), DMC the estimated MC dose value (in eV/g/shower), I the scanning current (in A) and t the acquisition time (in s).

Uncertainty on the estimated D is also evaluated combining the contributions to the uncertainty budget using the so called "sandwich law" described in the guide to the expression of uncertainty in measurement:

𝑢 𝑐 2 (𝐷) = 𝑢 2 (𝐷 𝑀𝐶 ) × (I × t × 1000) 2 + 𝑢 2 (I) × (𝐷 𝑀𝐶 × t × 1000) 2 +𝑢 2 (t) × (𝐷 𝑀𝐶 × I × 1000) 2 (𝑒𝑞𝑢. 8)
The values of the uncertainties u(DMC), u(I) and u(t) are detailled in paragraph 2.3.3.

Uncertainty budget

The stochastic uncertainties u(E) and u(DMC) due to MC simulation vary for each simulation. It is given with a coverage factor k =3. Uncertainty of the current (u(I)) and the acquisition time (u(t))

are given in the technical note of the CT :

𝑢(𝐼) = ∓(10% + 0.5 𝑚𝐴) (𝑒𝑞𝑢. 9) 𝑢(𝑡) = ∓(5% + 10𝑚𝑠)(𝑒𝑞𝑢. 10)
However, the confidence interval associated with these uncertainties is not specified. According to the recommendations of the guide for the expression of measurement uncertainty, we decided to consider that the uncertainties were expressed for a confidence interval of k = 1 and that the variable follows a uniform probability law. No information about u(V) and u(L) are given. According to the accuracy needed to build such ionization chamber, we assume that these uncertainties can be considered as negligeable.

Uncertainties are combined with the same confidence interval. In the following they are presented at k=2, ie with a confidence interval of ≈ 95%.

It is also considered that all the variables are independent that is to say that the covariance are not taken into account, this can lead to overestimate the uncertainties. It is assume that this overestimate is not large.

3 Model validation in homogeneous conditions

Isocenter validation

Integral of the air kerma over 100 mm is measured with a 10X6-3CT Radcal pencil ionization chamber (figure 3a). The chamber is introduced into a home-made PMMA tube of 3 cm exterior diameter. Acquisitions have been performed at 100 kVp and 120 kVp, with both bowtie filters, a 40 mm nominal collimation, a 50 cm beam length, 300 mAs and the X-ray tube positioned at the top of the gantry (figure 3b). 3.2 "Air kerma" index validation "Air kerma" index is a dose metric of the CT which represents the amount of radiations emitted per rotation of the RX tube. To determine "air kerma" index, two cylindrical phantoms (figure 4) composed of one central inserts and 4 peripheral inserts are used:

 The "head phantom" (figure 4a) has a diameter of 16 cm. It is used to calculate the air kerma index for head and children CT acquisitions.

 The "body phantom" (figure 4b) has a diameter of 32 cm. It is used to determine the air kerma index for body CT acquisitions

The air kerma index CTDIW (for Weighting Computed Tomographic Dose Index) is defined as:

𝐶𝑇𝐷𝐼 𝑊 = 1 3 (𝐶 𝑃𝑀𝑀𝐴,100,𝑐 + 2 × 𝐶 𝑃𝑀𝑀𝐴,100,𝑝 ) (𝑒𝑞𝑢. 11)
with 𝐶 𝑃𝑀𝑀𝐴,100,𝑐 being the value of the air kerma index obtained when the pencil chamber is localized in the central location and 𝐶 𝑃𝑀𝑀𝐴,100,𝑝 being the mean value of the air kerma index when the pencil chamber is inserted in the 4 peripheral locations.

Acquisitions have been performed at 100 kVp and 120 kVp and for a 40 mm nominal collimation.

Uncertainty budget

Measurements uncertainties are evaluated from information reported in the AIEA report 457 about the dosimetric practice in diagnostic imaging (IAEA-TRS 457). Scenario 1 has been chosen to estimate the uncertainty budget. In this scenario, the device is used in conformity with the CEI 61674 norm. According to the data reported in the AIEA report and to the fact that the pencil chamber is inserted in a PPMA phantom, experimental results are presented with a relative uncertainty of 13.3%.

Model validation in clinical conditions

Several acquisitions (table 4) are performed combining different parameters (high voltage, pitch, SFOV). The ATOM adult female phantom from CIRS is used to perform measurements with Optically Stimulated Luminescence dosimeters (OSLDs). This phantom has slots in different localizations and tissues to insert OSLDs (figure 5). In table 5 are reported the positions of the inserts in which

OSLDs have been inserted for both acquisitions (head and thorax).

OSL dose assessment

In their review, Yukihara and McKeever 2008, show the possibility of using OSLD for CT dosimetry. The great advantages of these detectors is their uniformity in sensitivity because the Al2O3:C powder used in the production process is a homogenized mixture of different crystal growth runs. The NanoDot TM (Landauer Inc.) contain single circular OSLD (5.0 mm in diameter)

placed in an adapter. An effective depth of 0.1 g/cm 3 is assumed as the point of measurement.

Yukihara and McKeever mention the OSLD energy dependence, showing variations between 20

and 30% for potentials ranging from 80 to 140 kVp. We have thus developed a specific user guide to take into account this energy dependence.

Detector calibration

The read-out is performed with the semi-automatic reader MicroStar TM NanoDot TM system.

Depending on the dose level, two light intensities are possible: for low doses (<200 mGy) all 38

LEDs are used and for high doses only 6 LEDs are used. For our application, the reader is always used in the low dose regime.

A calibration of our own OSLD has been carried out on a range from 0 to 150 mGy. with UK the uncertainty associated to the air kerma and defined as :

𝑈 𝐾 = ( 1 𝑎 𝑄 ) 2 × 𝑈 𝑆 2 + ( -𝑆 𝑎 𝑄 2 ) 2 × 𝑈 𝑎𝑄 2 (𝑒𝑞𝑢. 14)
where UaQ is the uncertainty associated to aQ the slope of the calibration function corresponding to the energy spectrum observed at the point of measurement and defined as :

𝑈 𝑎𝑞 2 = ( 1 ∑ 𝑁 𝑖 × 𝐸 𝑖 𝑖 ) 2 ∑(𝑁 𝑖 × 𝐸 𝑖 × 𝑈 𝑎(𝐸𝑖) )
2 (𝑒𝑞𝑢. 15) 𝑖 and US is the uncertainty associated to the OSLD read-out :

𝑈 𝑆 2 = 𝑈 𝑅 ̅ 𝐼 2 +𝑈 𝑅 ̅ 𝑁𝐼 2 (𝑒𝑞𝑢. 16)
with 𝑈 𝑅 ̅ 𝐼 being the uncertainty associated to the irradiated detectors and 𝑈 𝑅 ̅ 𝑁𝐼 being the uncertainty associated to the non-irradiated detectors ; both of them are calculated thanks to the following equation :

𝑈 𝑅 ̅ 2 = ( 1 𝑠𝑒 ) 2 × 𝑈 + ( 𝑈 𝑟̅ 𝑠𝑒 × 𝑠𝑒 ) 2 × 𝑈 𝑠𝑒 2 (𝑒𝑞𝑢. 17)
The uncertainty on the mass interaction coefficients is taken equal to 1%, according to the NIST (https://www.nist.gov/).

Dose maps assessment with XR-QA films

When exposed to radiation, the organic based dye of radiochromic films changes color due to polymerization: the color of XR-QA films turns from orange to brownish-black depending on the level of exposure (see figure 6). Several features of these detectors have attracted the attention of the medical physcis community: insensitivity to visible light, self-developing characteristics, dose-rate independence. showed that the use of a single calibration function leads to a relative uncertainty of 14% on the dose values, if the calibration function is obtained for a beam quality taken in the middle of the investigated energy range. In the following we have considered this value for the uncertainty associated to the film analysis.

Film calibration

The films calibration has been carried out in the LM2S (Laboratoire Modélisation, Systèmes et Simulation) laboratory for a 120 kVp X-ray beam with a HVL of 7.14 mm aluminum. Films and a Farmer 30013 PTW ionization chamber have been irradiated at the same time in order to determine the air kerma associated to the film read-out. The Farmer chamber has been previously calibrated at the French national laboratory of metrology (LNHB) in terms of air kerma. Films are read before and one week after irradiation and saved as TIF files. Unirradiated film reading is also 21 necessary to obtain net optical density. It is also recommended to control time between irradiation and read-out (at least 24 h). By using only the red channel, the mean pixel values before (𝑃𝑉 𝑁𝐼 )

and after (𝑃𝑉 𝐼 ) irradiation are calculated in a mean 1 mm 2 region of interest. The net optical density is then defined in equation 17:

𝑛𝑒𝑡𝑂𝐷 = 𝑙𝑜𝑔 10 ( 𝑃𝑉 𝑁𝐼 𝑃𝑉 𝐼 ) (𝑒𝑞𝑢. 18)
The calibration function linking the air kerma and the netOD has been adjusted according the following polynomial function (equation 18):

𝐾 𝑎𝑖𝑟,𝑄0 = 𝑎 × 𝑛𝑒𝑡𝑂𝐷 + 𝑏 × 𝑛𝑒𝑡𝑂𝐷 4 (𝑒𝑞𝑢. 19)
As for OSLD the dose in the tissue is obtained by multiplying the air kerma by the ratio of the mass interaction coefficients.

Film analysis

Gafchromic XR-QA2 films are irradiated to study the dose gradient. Films are cut to fit the anatomical shapes of the female anthropomorphic phantom and placed between two phantom slices. For thorax acquisitions, films cannot be inserted into the breasts of the phantom because breast are made from a single piece without any insert or slice. Stencil of the films contours are used to ensure the reproducibility of the film positioning during the reading steps before and after irradiation. For all acquisitions, the tube speed is fixed at 0.7 s/rot and the films are read four times before and one week after the irradiation. Films are read several times to ensure that film storage and handling have been performed in good conditions (dry and dark environment, no dust or fingerprints ..). Optical density values are then converted into air kerma according to the calibration function.

A thoracic and a head acquisition (table 7) have been performed.

Isocenter experimental validation

Table 8 shows the measured and simulated integral of the air kerma over 100 mm. Results show a good agreement between the simulations and the measurements with a deviation less than 10%

and an overlap larger than 87% for the 4 cases considered here. 9 shows the measured and simulated integral of the air kerma index obtained with both phantoms. Results show a good agreement between the simulations and the measurements with a deviation less than 4% and an overlap larger than 89% for the 4 cases considered here.

Dose maps comparison

Experimental and simulated dose distributions for a head scan are reported in figure 8. Large differences can be observed between the two maps and will be discussed later in the discussion section. 9 and obtained from figures 8a, 8b, 8c and 8d. In figures 8 and 9 the MC uncertainty is about 3% and the experimental uncertainty about 15%. Deviations between simulations and experiments are less than 20% in bone heterogeneities. In soft tissue regions the deviation can be higher, especially for the head acquisition. This deviation is mainly due to the path of the x-ray tube during the experimental scan which is not simulated accurately in the MC simulation, the initial position of the tube being experimentally unavailable. this also validates the use of the tube load information to convert simulated results into Gray. By comparing integral of the air kerma in table 8 and their associated uncertainties, we note that the simulation uncertainties budget is actually deteriorated by the conversion factor contribution.

According to the manufacturer technical note, current and acquisition time have 10% and 5%

uncertainty on the displayed value, respectively. Despite these values, relative uncertainties are below 15%. Such uncertainties are compatible with medical imaging applications.

For directly irradiated areas there are sometimes some differences in the vicinity of bone structures. They are mainly due to the different pixel sizes between the film and the dose matrices.

Due to the small number of particles, the size of the voxels in the dose matrix cannot match the size of that of the film. Indeed, the smaller the size of the voxels, the more one has to increase the number of particles to converge the simulation.

In addition to the difficulties related to the difference in resolution, the dose maps from the simulation highlight problems related to voxelization of the phantom. Since the voxels are larger than those used for the phantom, a voxel in the dose grid can be composed of several tissues (lung, bone and soft tissue in our case). If a voxel is composed of several tissues, it can be considered as being composed of an hybrid tissue associated with an intermediate density according to the densities of the materials initially present and their density. The dose deposit is then affected and the separation between the tissues is less marked. However, the dose profiles show that the gradients are still well respected.

Besides, in their article Long et al 2013 showed from MC simulations that the starting angle could lead to organ dose differences between -20% and 34% compared to the average value. We also found such discrepancies when comparing simulations and measurements obtained with OSLD.

By combining the information of Long's article and the non-homogeneous dose deposit visible on the films due to the random draw of the starting position, the important differences found in the comparison between measurement and simulation for the OSLDs can be explained.

For one studied case, organ dose estimations with our software are in agreement with those published by Zhang et al 2012, attesting the reliability of the developed software. Organ dose estimation in the ICRP 110 phantoms can be thus performed in a short notice, less than one hour using 24 CPU. In the future, improvements would be considered to reduce the simulation time.

  well as Li et al 2011, choose to insert TLD (ThermoLuminescent Dosimeter) chips into Alderson-Rando phantom.

  Figure 1: Information about the CT geometry (a), the PSF location (a), the circular splitting (b) and the

Figure 2 :

 2 Figure 2: Calibration function for the DICOM images of the CIRS phantom.

Figure 3 :

 3 Figure 3: Longitudinal section of the PMMA phantom (a) and experimental set-up diagram (b).

  The bowtie filter corresponding to the phantom has been used, thus the Ped Body filter for the head phantom and the Large Body filter for the body phantom. All acquisitions have been performed with the Radcal pencil chamber, 600 mA and an acquisition time of 1 s. The phantom is positioned at the CT isocenter. The pencil chamber is positioned in one of the insert while the other inserts are filled with PMMA cylinders. Five acquisitions are needed to get the air kerma index in the 5 inserts of each phantom.

Figure 4 :

 4 Figure 4: Layout of the head (a) and body (b) phantoms used for measuring CTDIW.

Figure 5 :

 5 Figure 5: Photo of an OSLD on the left and photo of an OSLD inserted in a CIRS phantom slot on the

Figure 6 :

 6 Figure 6: Film read-out for the head acquisition: the color of XR-QA films turns from orange to

Figure 7 :

 7 Figure 7: Experimental dose values versus simulated dose values and uncertainties given at k=2 for the head (in green) and thorax (in blue) acquisitions with the female ATOM phantom. The relation between both dose values are fitted (in red) as well as the uncertainty gap (in black).

Figure 8 :

 8 Figure 8: Simulated and measured dose distributions in mGy for a head (c and d) and a thorax (a and b) scan. Dashed lines indicate the profiles used for figure 9. Validation of the dose gradient is performed by comparing simulated and experimental profiles reported in figure 9 and obtained from figures 8a, 8b, 8c and 8d. In figures 8 and 9 the MC

Figure 9 :

 9 Figure 9: Simulated and measured relative dose profiles for the head (c and d) and the thorax (a and b) acquisitions. The profiles are obtained for the dashed lines reported in figure 8. The deviation between simulation and measurement is reported in blue.7.5 Organ dose estimation

Table 1 : Filtration of the studied beams.

 1

	Potential (kVp) SFOV (cm) Total filtration (inherent + central bowtie filtration).
	100	5.20 mm Al
	50
	120	4.65 mm Al
	100	3.44 mm Al
	32
	120	2.89 mm Al

Table 2 : MC simulation parameters used for the homogeneous acquisitions.

 2

	Material Eabs,e-	Eabs,photons (eV) Eabs,e+	C1	C2	WCC (eV) WCR (eV)
		(eV)		(eV)
	Air	10 3	10 2	10 4	0.05	0.05	10 2	10 2
	PMMA	10 3	10 2	10 4	0.05	0.05	10 2	10 2

Table 3 : Parameters used for PENELOPE simulations and for each biological material used in this study.

 3

	Material	Eabs,e-(eV) Eabs,photons (eV) Eabs,e+ (eV) C1	C2	WCC (eV) WCR (eV)

Table 4 : Acquisitions performed experimentally and with the MC simulation to estimate point dose values.

 4

	kVp	Bowtie filter	Pitch	mAs
			0.984	570
		Large
	100		1.375	750
		Ped	0.516	95
	120	Ped	0.984	150
	4.1 CIRS anthropomorphic phantom

Table 5 : OSLD locations in the anthropomorphic phantom for thorax and head acquisitions.

 5

		Location	Tissue	Section	Position
		Stomach	Soft tissue	21	137
		Right Lung	Lung	20	107
	Thorax	Spine	Bone	19	102
		Heart	Soft tissue	16	83
		Rib	Bone	15	69
		Brain	Soft tissue	3	6
		Skull	Bone	3	10
	Head
		Mandible	Bone	7	18
		Thyroid	Soft tissue	10	28

Table 6 : Calibration function (OSLD reading as a function of air kerma) obtained for the different beam qualities.

 6 reading value for irradiated (𝑅 𝐼 ̅̅̅ ) and non-irradiated (𝑅 𝑁𝐼 ̅̅̅̅̅ ) OSLD, respectively. Finally, the corrected signal S, used to determine the dose absorbed by the OSLD is calculated by subtracting raw reading values for irradiated (𝑅 𝐼 ̅̅̅ ) and non-irradiated (𝑅 𝑁𝐼 ̅̅̅̅̅ ) OSLD.To take into account the energy dependence of the OSL detectors, protocol detailed in Bordy et al 2013 is adapted. The energy spectrum corresponding to the OSL position is determined by MC simulation. Calibration factors for each energy bins are convolved using the energy spectrum as weight to adapt the calibration factor to the spectrum at the point of measurement. Assuming the electronic equilibrium condition are fulfilled, the dose in the tissue is obtained by multiplying the air kerma by the ratio of the interaction coefficients:

	Five radiation

used. Five detectors have been irradiated for each air kerma value. Average OSLD readings are used to calculate the calibration factor for each beam quality, assuming that OSLD response is linear in terms of air kerma. Results are reported in table

6

.

Mean reading values for irradiated (𝑟 𝐼 ̅) and non-irradiated OSLD (𝑟 𝑁𝐼 ̅̅̅̅) are then computed. These mean raw reading values are then corrected individually from the OSLD sensitivity (se) to obtain the real

Table 7 : Acquisitions performed experimentally and with the MC simulation to estimate dose maps. kVp Bowtie filter Pitch Position in the phantom

 7

	100	Large	0.516	Between section 16 and 17
	120	Ped	1.375	Between section 7 and 8
	5 Organ dose estimation
	Organ doses for a thoracic localisation have then been estimated into the female phantom
	provided in ICRP Publication 110. Simulations are performed for the Large Body bowtie filter with
	a 50 cm SFOV, 100 and 120 kVp, a 40 mm collimation adapted to the used focal spot (42.6 and
	42.9 mm), three mAs (100, 200 and 300), a 0.7 s/rot tube speed and three pitch values (0.531 -
	0.969 -1.375). Dose values to each voxel across all the voxels belonging to each organ are
	averaged. According to the organ wT factors reported in the ICRP 103, we decid to report the dose
	absorbed by the more sensitive organs thus the left breast glandular tissue, the stomach wall, the
	left pulmonary tissue, the esophagus and the spinal cord.
	6 Comparison of experimental and simulated results
	To compare experimental results (rexp, σexp) and MC estimation (rMC, σMC) we have used 2 index :
	 the deviation defined as :
			𝑑𝑒𝑣 =	𝑟 𝑀𝐶 -𝑟 𝑒𝑥𝑝 𝑟 𝑒𝑥𝑝	× 100 (𝑒𝑞𝑢. 21)
	 and the overlap defined as :

Table 8 : Measured and simulated values obtained for the integral of the air kerma in the pencil chamber.

 8

				MC estimation	Experimental value	Comparison
				DI ± UDI (mGy.cm)	dev (%)	ovlp (%)
	Ped	Body	100 kV	289.9 ± 35.1	278.6 ± 18.4	4.0	96.0
	Bowie filter
			120 kV	400.6 ± 48.4	408.2 ± 26.9	-1.9	99.1
	Large Body	100 kV	206.7 ± 25.5	216.1 ± 14.3	-4.4	94.9
	Bowie filter
			120 kV	309.1 ± 37.8	331.4 ± 21.9	-6.7	87.8
	7.2 Validation of the air kerma index estimation
	Table

Table 9 : Measured and simulated air kerma index obtained with both phantoms.

 9

				MC estimation	Experimental value	Comparison
				CTDIw ± UCW (mGy)	dev (%)	ovlp (%)
	Ped	Body	100 kV	76.3 ± 4.4	74.2 ± 4.9	2.8	95.0
	Bowie filter
			120 kV	112.9 ± 6.4	117.6 ± 7.8	-4.0	89.4
	Large Body	100 kV	37.7 ± 2.6	36.5 ± 2.4	3.2	94.7
	Bowie filter
			120 kV	58.7 ± 4.0	59.0 ± 3.9	-0.5	99.8
	7.3 Point dose comparison
	Experimental and simulated dose values are reported in figure 7 as well as the relation between
	them:
				𝐷 𝑒𝑥𝑝 = 0.867 × 𝐷 𝑠𝑖𝑚𝑢𝑙 (𝑒𝑞𝑢. 23)
	Simulated values are on average higher than the experimental ones. The uncertainty bars plotted
	in figure 7 is obtained by fitting the MC uncertainty associated to each simulated values. This one
	is about 23.2% (k=2) for all the simulations performed here. The experimental uncertainty is about
	7.4%.

  table 10 and compared with Zhang et al 2012 results, who also estimated organ dose for the female phantom described in the ICRP 110. They simulated in details the GE VCT LightSpeed 64 thanks to accurate data provided by the manufacturer. They used a modified version of PENELOPE reported in Li et al 2011 to estimate organ doses in the ICRP 110 phantoms. Our results, obtained in less than one hour in 24 CPU for all cases, are a little bit larger than the

dose values reported by Zhang et al 2012. Deviations are less than 5.8% for the four organs considered here.

Table 10 : Comparison of organ doses for a thoracic helical acquisition obtained with our MC tool and reported by Zhang et al 2012.

 10 Despite a lack of information about the scanner geometry, the GE VCT LightSpeed 64 has been modelled by adapting the method developed by Turner et al 2009.Results obtained in homogeneous conditions validate the use of the MC model for dosimetric estimation. Measured and simulated integrals of the air kerma over 100 mm are in agreement;

		Simulated organ doses (mGy)
	Organs			Deviations(%)
		This study Zhang et al 2012
	Esophagus 9.8 (0.003)	9.3	5.4
	Lung	10.9 (0.003)	10.3	5.8
	Breast	10.2 (0.002)	10.1	1.0
	Stomach	11.1 (0.001)	10.9	1.8
	8 Discussion and Conclusion

http://ct-imaging.de/en/ct-software-e/impactdose-e.html

http://www.impactscan.org/

Acknowledgments

The authors thank warmly Helena Chesneau for the calibration of the Gafchromic films and Fabien Moignau, Marc Denoziere, Nelly Lecerf, for their help in the calibration of the OSLD.

Measured and simulated point dose obtained in anthropomorphic conditions show deviations up to 15%. However confidence intervals are overlapped allowing us to conclude that results are in agreement. The uncertainty budget for simulated doses is mainly by the conversion factor uncertainties. The relative uncertainties for the tube current and the acquisition time are respectively 10% and 5% at k=1, as reported in the technical note. Simulated dose uncertainties might seem quite large (about 22% at k=2), but such uncertainties are compatible with dosimetric purposes in medical imaging. The benefit of the conversion factor is therefore maintained.

Large differences in the simulated and the experimental dose distributions can be noted. All the experimental dose maps show an important effect of the initial tube position, as already reported by several authors [START_REF] Li | Patient-specific radiation dose and cancer risk estimation in CT: Part I. Development and validation of a Monte Carlo program[END_REF][START_REF] Zhang | Reducing radiation dose to selected organs by selecting the tube start angle in MDCT helical scans: A Monte Carlo based study[END_REF]. The surface dose distribution resulting from a helical acquisition is periodic [START_REF] Dixon | Experimental validation of a versatile system of CT dosimetry using a conventional ion chamber: Beyond CTDI [sub 100[END_REF] and the tube start angle determines the location of the high and low dose regions. It has been reported by [START_REF] Zhang | Reducing radiation dose to selected organs by selecting the tube start angle in MDCT helical scans: A Monte Carlo based study[END_REF] that the magnitude of organ dose reduction resulting from varying tube start angle varies from 10 to 30% depending on the location and size of the organs. In the experimental dose distribution, shown in figures 8b and 8d, the tube position relative to the phantom slice containing the film can be easily determined, since a higher dose is delivered to the top and the back of the phantom for the thorax and the head acquisition, respectively. However, the tube path relative to the patient cannot be fully worked out because the tube starting angle information is not provided on the GE VCT Lightspeed 64 scanner. Taking into account this lack of information it has been decided to randomly sample the tube starting angle for each simulated particle. By making this choice, all possible tube paths are simulated, leading to a more homogenous dose distribution (figures 8a and 8c) and an overestimation for some location of the real delivered dose. Nevertheless, instead of underestimating the dose for radiosensitive organs the MC simulation considers the worst case and provides the maximum dose which could be delivered.