Maxim Moraru

Mina Warnet

Julien Loiseau

Vinay Ramakrishnaiah

Nirmal Prajapati

Hyun Lim

Sumathi Lakshmiranganatha

Jamal Mohd-Yusof

Karen Tsai

Richard Berger

Patrick Mccormick

Transformations for Energy Efficient Accelerated Chain Matrix Multiplication (TEE-ACM 2 )

Matrix Chain Multiplication plays a key role in the training of deep learning models. They also appear in physics, computer graphics, image processing, etc. Matrix Multiplications often cause a bottleneck in terms of performance and energy because of the heavy costs in computations and memory operations. While the runtime performance has been studied for years, significantly less effort has been expended into optimizing energy efficiency.

Most processors weren't built with energy in mind [3], only performance. For a supercomputer, the cost of power is the most significant expense, and even with the most affordable energy, the cost of a MW ends up at $1M. Moreover, the Exascale Computing Project [4] defines a 'capable' exascale system as one that operates in a power envelope of 20-30 MW.

Thus, reducing the energy cost of these types of computations is a major challenge.

Matrix Chain Multiplication

Problem: Given a sequence of matrices {A 1, A 2, ... A n } with sizes {P 0, P 1, ... P n } , compute ς 𝑘=1 𝑛 𝐴 𝑘 • Multiplication order can significantly impact the performance of the algorithm • The Optimal Parenthesization (OP_Count) algorithm in "Cormen et al. Introduction to Algorithms." outputs an order of matrix multiplications that minimizes the total number of operations.

Energy Efficient GPU implementation

• Energy consumption of a GPU can be broken down into 2 major parts:

• Energy of the operations executing on the GPU • Energy of the GPU itself • Memory operations dominate the executed operations

Our goal : optimize the total energy consumption for Matrix Chain Multiplications

Introduction

• Less power consumption may be preferred over shorter execution times for applications such as edge computing. This work focuses on such use cases • For our fused kernel implementation, the energy consumption was primarily dictated by the execution time (static power) • Performing fewer global loads and stores significantly decreases the GPU power consumption at the cost of increased execution time • By using the OP_Count tree we optimize the computation order in a sequence of matrices.

• The hardware counters and APIs supplied by Nvidia provide accurate results with precision comparable to the measurements made using the hardware power capture analysis tool (PCAT)

• Future work: we plan to optimize the fused kernels by removing atomic operations, thus reducing computation time. Once the fused kernels are implemented, we will generate the fusion tree and compare the time and energy between TEE-ACM 2 with and without fusion

Conclusions and Future Work

Extending the approach to three matrices Optimal Algorithm for Matrix Chain Multiplication

• The OP_Count algorithm produces a tree decomposing a matrix chain, which indicates the order of matrix products to minimize operations • OP_DM_Fuse is created from OP_Count, which reduces the number of operations and then minimizes the off-chip data transfers by using fusion • Minimize over the options: left-fuse, right-fuse, and no-fusion Example for fusion decision algorithm • The output tree consists of the minimum number of computations and data transfers for multiplying the given matrix chain Fused Matrix Multiplication [START_REF] Cormen | Introduction to algorithms[END_REF] P 0 Energy efficiency

P 1 P 2 x A 1 A 2 A 3

GPU resources usage

• TEE-ACM 2 uses more than two times on-chip memory than cuBLAS

The impact of varying the tile size

• By using larger tile sizes, the number of global loads decrease significantly • cuBLAS performs 3 times more global loads on average than TEE-ACM 2

Single Matrix Multiplication Results

Single Matrix Multiplication : Blocking Strategy Single Matrix Multiplication [START_REF] Cormen | Introduction to algorithms[END_REF] Theorem: Given the matrix sizes -P 0 ,P 1 ,P 2 and shared memory capacity -M, the minimum number of data transfers is given by:

Fused-Kernel Results

*Present results were obtained on a Darwin node equipped with v100 GPUs and Cascade Lake CPUs (averaged over 10 samples of 10 runs each)

• Our implementation of fused-kernel approach uses atomic adds that drastically impact the execution time • Also, when multiplying three matrices, the GPU resources are distributed across 4 tiles compared to 3 in the case of a single matrix multiplication

R = A 1 A 2 A 3 = 10×30×5 + 10×5×60 = 4,500 multiplications (A 1 A 2 ) A 3
Sizes 10×30, 30×5, 5×60

= 27,000 multiplications • When computing with the best case, cuBLAS without the tree is 2% faster because of the overhead of building and traversing the tree • In the worst case, the creation of the matrices is longer without the tree because the parenthesizing isn't optimal, the temporary matrices are bigger, thus the longer execution time

A 1 (A 2 A 3 ) = 30×5×60 + 10×30×60
• The blue percentages indicate the energy consumption ratio between TEE-ACM 2 and cuBLAS using OP_Count. By using our optimized kernel, we end up with 5 to 10% energy savings with the optimal tile sizes of 224 • The orange percentages indicate the energy consumption ratio between cuBLAS with and without the tree. In the best case, we consume 1.5% more energy with the tree, but in the worst case, using the tree allows us to save up to 50% energy on average

•

  Input: 5 matrices with sizes 936, 1008, 552, 368, 1016, 616, 544 and M = 65,536 1. Visit each node of the OP_Count tree 2. Consider options: no fuse, left fuse, right fuse 3. Calculate the minimum number of data transfers for each node 4. Use the optimal over all possible options

Figure 11 :

 11 Figure 11 : Energy savings of TEE-ACM 2 over cuBLAS (lower is better)

Figure 8 :

 8 Figure 8 : Installation of PCAT device on GPU for hardware power measurement

Figure 1 :Figure 2 Figure 3 :Figure 4 :

 1234 Figure 1: Data transfers vs arithmetic operations (Bill Dally, Challenges of Future Computing Systems, HiPEAC 2015)

Figure 5 :

 5 Figure 5: Fused matrix multiplication approach, where three matrices are multiplied in a single GPU kernel.

Figure 6 :

 6 Figure 6: Left fuse strategy

Figure 7 :Figure 16 :

 716 Figure 7: Left fuse strategy

Figure 20 :

 20 Figure 20 : Energy consumption of computation for a sequence of 6 matrices


Julia Bauer
email: julia.bauer@med.uni-heidelberg.de

Meret Hildebrandt

Michael Baumgartl

Fine Fiedler

Charlotte Robert

Irène Buvat

Wolfgang Enghardt

Katia Parodi

Quantitative assessment of radionuclide production yields in in-beam and offline PET measurements at different proton irradiation facilities

Keywords:

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Objective: Reliable radionuclide production-yield data are a prerequisite for positron-emissiontomography (PET) based in-vivo proton treatment verification. In this context, activation data acquired at two different treatment facilities with different imaging systems were analyzed to provide experimentally determined radionuclide yields in thick targets and were compared with each other to investigate the impact of the respective imaging technique.

Approach: Homogeneous thick targets made of three different materials (PMMA, gelatine, and graphite) were irradiated with mono-energetic proton pencil beams at two distinct energies. Material activation was measured (i) in-beam during and after beam delivery with a double-head prototype PET camera and (ii) offline shortly after beam delivery with a commercial full-ring PET/CT scanner. Integral as well as depth-resolved b + -emitter yields were calculated from the PET data for the dominant positron-emitting radionuclides 11 C, 15 O, 13 N and (in-beam only) 10 C. Moreover, the in-beam data were used to investigate the qualitative impact of the time scheme of common monitoring workflows on activity depth profiles and their quantitative impact on count rates and total activity.

Introduction

Ion beam therapy, in particular proton therapy, has nowadays become a recognized treatment modality complementing conventional external photon irradiation for the treatment of various cancer diseases. Due to the characteristic inverted depth dose profile of ions stopping in the irradiated tissue, the beam range needs to be well controlled in the patient for a safe delivery of highly conformal dose distributions to complex target volumes. Range uncertainties can be minimized either by using high-energetic transmission beams to accurately determine the particle stopping power of the tissue [START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF], or by secondary radiation detection to reconstruct the range of the primary beam particles [START_REF] Knopf | In vivo proton range verification: a review[END_REF]. Based on the latter concept, two different techniques for in-vivo beam range monitoring have so far been studied clinically: (i) prompt gamma imaging (PGI) [START_REF] Richter | First clinical application of a prompt gamma based in vivo proton range verification system[END_REF][START_REF] Xie | Prompt Gamma Imaging for In Vivo Range Verification of Pencil Beam Scanning Proton Therapy[END_REF] and (ii) positron-emission-tomography (PET) imaging [START_REF] Ferrero | Online proton therapy monitoring: clinical test of a Silicon-photodetector-based in-beam PET[END_REF][START_REF] Min | Clinical application of in-room positron emission tomography for in vivo treatment monitoring in proton radiation therapy[END_REF][START_REF] Miyatake | Measurement and verification of positron emitter nuclei generated at each treatment site by target nuclear fragment reactions in proton therapy[END_REF][START_REF] Nischwitz | Clinical implementation and range evaluation of in vivo PET dosimetry for particle irradiation in patients with primary glioma[END_REF][START_REF] Nishio | The development and clinical use of a beam ON-LINE PET system mounted on a rotating gantry port in proton therapy[END_REF]Parodi et al., 2007b). While PGI operates during beam delivery, the PET-based approach can be implemented in different schemes [START_REF] Shakirin | Implementation and workflow for PET monitoring of therapeutic ion irradiation: a comparison of in-beam, in-room, and off-line techniques[END_REF]: for in-beam monitoring the PET camera is installed directly at the treatment place and detects the decay of irradiation induced PET-radionuclides during patient irradiation; in-room installations operate a PET-scanner in the treatment room but not integrated in the beam nozzle; offline installations use a nearby PETscanner located outside the treatment room. In-beam installations are superior concerning the signal rate, thus requiring shorter acquisition time. However, the high demands of geometric compatibility of the camera with the (often rotating) beam nozzle have so far only been met by the use of limitedangle double-head cameras, with considerably affected imaging quality compared to full-ring PET cameras [START_REF] Crespo | On the detector arrangement for in-beam PET for hadron therapy monitoring[END_REF]. Offline installations, on the other hand, usually operate commercial full-ring PET-scanners with state-of-the-art imaging performance. However, the signal rate is low and spatial information is blurred due to physical decay and biological wash-out of the patient activation during the time delay caused by patient transport and repositioning after beam delivery.

Besides imaging performance and signal quality aspects, a reliable range verification is furthermore challenged by uncertainties in the production cross-section data of relevant radionuclides [START_REF] Seravalli | Monte Carlo calculations of positron emitter yields in proton radiotherapy[END_REF]. These data are required to predict the expected patient activation and thus to assess the observed beam range. An experimental determination of radionuclide-specific production yields is therefore an important contribution to increase the accuracy of the verification strategy.

Treatment centers implementing a PET-based verification concept usually conduct various preclinical activity measurements in phantoms for system commissioning (Bauer et al., 2013a;[START_REF] Bisogni | INSIDE in-beam positron emission tomography system for particle range monitoring in hadrontherapy[END_REF][START_REF] Horst | Measurement of PET isotope production cross sections for protons and carbon ions on carbon and oxygen targets for applications in particle therapy range verification[END_REF][START_REF] Meißner | Towards offline PET monitoring of proton therapy at MedAustron[END_REF]Parodi et al., 2007a;[START_REF] Shao | In-beam PET imaging for on-line adaptive proton therapy: an initial phantom study[END_REF][START_REF] Zhu | Monitoring proton radiation therapy with in-room PET imaging[END_REF]). An inter-facility comparison of the obtained production yield data, however, has so far been impeded by differences in the respective irradiation setup. To overcome this issue, an experimental campaign was conducted at the Heidelberg Ion-Beam Therapy Center (HIT, Heidelberg, Germany) under irradiation conditions comparable to a previously conducted phantom activation measurement performed at GSI (Helmholtzzentrum für Schwerionenforschung, Darmstadt, Germany). At both facilities activation data were acquired in different homogeneous block-shaped thick targets with different PET imaging concepts: an in-beam PET camera installed at GSI [START_REF] Enghardt | Dose quantification from in-beam positron emission tomography[END_REF] and an offline PET/CT scanner installed at HIT (Bauer et al., 2013b). With these data sets we were able to:

(i)
independently assess radionuclide production yield data obtained at different facilities (ii) study the qualitative and quantitative impact of the respective imaging technique Both aspects are considered relevant for the commissioning of future PET-monitoring installations and workflow implementations relying on activity predictions for beam range evaluation. While activation data obtained at the two irradiation facilities HIT and GSI for ions heavier than protons have already been published separately [START_REF] Bauer | Offline imaging of positron emitters induced by therapeutic helium, carbon and oxygen ion beams with a full-ring PET/CT scanner: experiments in reference targets[END_REF][START_REF] Fiedler | Anwendung des in-beam PET Therapiemonitorings auf Präzisionsbestrahlungen mit Helium-Ionen[END_REF][START_REF] Priegnitz | An Experiment-Based Approach for Predicting Positron Emitter Distributions Produced During Therapeutic Ion Irradiation[END_REF][START_REF] Sommerer | In-beam PET monitoring of mono-energetic 16 O and 12 C beams: experiments and FLUKA simulations for homogeneous targets[END_REF], this study reports the resulting integral and depth-resolved radionuclidespecific production yields induced by therapeutic proton beams and, for the first time, an interfacility comparison of these datasets.

Materials and Methods

Phantom description

Table 1 provides an overview of relevant PET radionuclides ( 11 C, 10 C, 13 N, 15 O) produced in therapeutic human tissue irradiation, the dominant production channels, and their respective half-life. Phantoms made of gelatine (GEL) and graphite (GR) were used to study the respective production channels separately, while the used PMMA (polymethyl methacrylate, C5H8O2) phantoms are considered a good surrogate for soft human tissue. Similar rectangular shaped homogeneous phantoms were used at both centers, GSI and HIT. They were block-shaped with dimensions between 15 cm to 35 cm in length and (5-10) cm in width and height. Geometrical details of the phantoms as well as the respective photon attenuation coefficients required for the GSI image reconstruction pipeline are reported in Table 2.

Gelatine phantoms at GSI were put on an acrylic glass plate of (2-5) mm thickness, respectively, for irradiation and imaging to cope with their unstable structure. The plate was considered for activity calculation. The gelatine phantoms at HIT were prepared in a dedicated PMMA box of 1 cm wall thickness and kept in there to ensure its stability during transportation to and repositioning at the offline PET/CT scanner after irradiation.

Proton irradiation setup

At GSI, the phantoms were positioned at the experimental treatment place on the treatment table

with the long side aligned to the beam direction for irradiation. Mono-energetic pencil-like proton beams were extracted from the synchrotron at a lower energy of 125.4 MeV/u (E1) in 100 spills and a higher energy of 176.95 MeV/u (E2) in 120 spills. The resulting irradiation times (tirr) are reported in Table 3 together with the number of delivered protons (Np) and the calculated intensity of the beam (I). The time structure of the beam delivery was reasonably stable with mean spill and pause times of

𝑡 ! # = 1.
3 𝑠 and 𝑡 " # = 3.0 𝑠, respectively, with variations below 11%. At the time of experiments, a 4 mm thick protective PMMA plate, which corresponds to a range shift of 4.6 mm water-equivalent path length (WEPL), had accidentally been forgotten in the GSI beam line during all phantom irradiation experiments. Thus, the energy of the primaries reaching the phantoms was slightly reduced and the production yields were expected to be lower. On the other hand, the WEPL of the GSI beam nozzle is about 1mm smaller as compared to the HIT beam nozzle. The impact of these two factors on the expected production yield was estimated with a Monte Carlo simulation using the FLUKA-based framework of [START_REF] Bauer | Integration and evaluation of automated Monte Carlo simulations in the clinical practice of scanned proton and carbon ion beam therapy[END_REF].

Phantom irradiation at HIT was performed under comparable experimental conditions (Bauer et al., 2013a). The phantoms were irradiated in the experimental room and positioned with the long side aligned to the beam direction such that the geometrical center of the phantoms coincided with the isocenter of the experimental room. Mono-energetic pencil-like proton beams of two energies, 125.67 MeV/u (E1) and 176.75 MeV/u (E2), respectively, were delivered over a time period of about 3 minutes by the synchrotron. Since the energy difference of about 0.2 MeV/u between the GSI and the HIT experiments is considered negligible, the labels E1 and E2 are in the following used for both experimental series. Beam delivery at HIT was stable within negligible deviations: the spill length and pause times between the spills were on average 𝑡 ! # = 1.7 𝑠 and 𝑡 " # = 4.3 𝑠 with variations of below 5% and 2%, respectively (Bauer et al., 2013a).

PET data acquisition

Imaging concepts as well as workflow implementation were fundamentally different at the two facilities, as sketched in Figure 1. At GSI, a double-head PET camera system was directly installed at the treatment place. This setup allowed for PET data acquisition during irradiation (in-beam) and after beam delivery (offline) without the need for transportation or repositioning of the phantoms.

Each camera head consisted of 4 x 8 position-sensitive scintillation block detectors made of bismuth germanate (BGO). The centre of the PET camera field-of-view (FOV), where the imaging performance is best, coincided with the isocentre of the treatment room. Therefore, the phantoms were positioned such that the area of the expected distal activity fall-off, thus the activity maximum, was located at the isocenter (cf. Figure 1 left). Correct positioning was checked by means of a point sourced mounted on a reference point on the phantom prior to the irradiation [START_REF] Parodi | In-beam PET measurements of β + radioactivity induced by proton beams[END_REF]. A synchronization of the PET data acquisition with the beam control system, mapping the time course of irradiation and the acquired PET data, was used to suppress prompt gamma background radiation.

PET data acquired in-spill could thus be discarded later on for the analysis. In order to study the activity decay, the in-beam PET acquisition was prolonged for 25 to 35 minutes after beam delivery.

At HIT, the irradiated phantoms were immediately transported to a nearby installed commercial fullring PET/CT scanner (Siemens Biograph mCT ®) based on lutetium oxyorthosilicate (LSO)-crystal technology. PET measurements began with a typical time delay Dt of 3 to 4 minutes after beam delivery (cf. Table 3), needed for transport and repositioning of the phantom. At the PET/CT, the phantoms were positioned with the long side transversal to the scanner axis, centered in the FOV using the internal laser alignment system and the indexed patient support table. On the measurement day, the PET/CT scanner was calibrated according to the standard procedure of the daily quality assurance to ensure the quantification reproducibility at clinically required level. The examination protocol was set up such that the 30 min PET acquisition directly followed the initial topogram scan, while the CT scan necessary to derive the attenuation and scatter correction of the PET data was performed at the end of the examination. In order to eliminate the scanner applied correction for the decay of the radionuclide as required for diagnostic PET applications, the administered radionuclide was set to 22 Na with a half-life of about 2.6 years, considering the b + branching ratio of this radionuclide of 90% later on in the data analysis.

Time structure information relevant for the PET imaging of both experimental series is reported in Table 3. At both facilities, PET data were acquired in "list mode", imposing a time stamp every millisecond (every 10 ms) on the coincidences record acquired at HIT (GSI), which allows for a retrospective reconstruction of counts and PET images in arbitrary time frames.

Data analysis of the GSI experiment

The GSI data analysis pipeline has in general been adopted from [START_REF] Parodi | In-beam PET measurements of β + radioactivity induced by proton beams[END_REF]. In order to determine the absolute production yield, the signal counts detected by the GSI PET camera have to be reconstructed in several subsequent steps. First, the activity decay measured after beam delivery is used to extrapolate the radionuclide-specific activity present at the end of the irradiation. From this point, it is possible to calculate the radionuclide-specific production rate per second and eventually the radionuclide production yield per primary incoming proton considering the relevant irradiation parameters. Since the GSI analysis chain is based on raw PET camera data, the obtained production yields are relative numbers. Additional corrections for detector efficiency, annihilation photon attenuation and geometrical acceptance due to the solid-angle camera field-of-view need to be considered to obtain absolute production yields.

Fit of activity decay

The radionuclide-specific contribution to the total activity present at the end of irradiation, A0,i (i = 10 C, 15 O, 13 N and 11 C) can in general be determined by a nonlinear fit of the measured activity decay according to the law of radioactive decay:

𝐴(𝑡) = / A #,% • 𝑒 ('( ! ⋅*) % (1)
with li being the half-life of the respective radionuclide. A customized software coded in IDL (Interactive Data Language, Version 5.3) was used for the regression analysis.

Various consistency checks were performed to ensure the reliability of the fitting procedure: (i) regression with fixed or variable radionuclide half-lives, (ii) regression using the full activity decay or only a selected time interval optimized for reasonable radionuclide abundances giving the best results for the fitted half-lives compared to the expected half-lives of the radionuclides, and (iii)

applying the regression code used for the HIT data (cf. section 2.5) as cross-check. The intermediate results of these checks have been documented in [START_REF] Hildebrandt | Experimental investigation of beta plus emitter yields induced by proton beams in different materials for improved modeling of in-vivo PET verification[END_REF], the final activity values resulting from the investigation of Hildebrandt are summarized in Table 4 together with the deviations found for the cross-check. These A0 values were used as input for the subsequent analysis pipeline (cf. equation 2). The radionuclide-specific production rate Pi can be calculated from the activity present at the end of irradiation, A0,i (cf. Table 4), considering the recorded time structure of the pulsed synchrotron irradiation, with beam delivery during spills (ts) and pauses in between (tp). Two different approaches have been pursued to calculate relative production rates Pi, analogous to [START_REF] Parodi | In-beam PET measurements of β + radioactivity induced by proton beams[END_REF]. The first approach assumes of a constant spill structure over Ns spills, thus using mean spill (Ts) and mean pause times (Tp), and a constant beam intensity. The production rate Pi is then given by:

𝑃 % = 𝐴 #,% (1 -𝑒 ', ! ⋅-" ) ⋅ ∑ 𝑒 ', ! ⋅(.'/)⋅(-" 0-# ) 1 " .2/
(2)

A second approach considers the activity Aj measured at the end of each spill and the relative contribution fi of the i th radionuclide to the estimated amount of activity 𝒜 𝒿 produced in the j th spill, thus allowing energy and time variations during irradiation to be considered. The production rate is then calculated as follows:

𝑃 %,. = 𝒜 𝒿 ⋅ 𝑓 % 1 -𝑒 '(⋅* ",% (3)
with

𝒜 𝒿 = 𝐴 . -/ / 𝑃 %,4 ⋅ 91 -𝑒 '(⋅* ",& : .'/ 42/ % ⋅ 𝑒 '( ' ⋅∑ 6* ",()& 0* #,()&*+ 7 %*+ (,+ (4)
The factor fi is a weighting factor for the contribution of each considered radionuclide of the activity in each spill and is calculated as follows:

𝑓 % = 𝑃 % ⋅ 91 -𝑒 '( ! * " : ∑ 9𝑃 % ⋅ (1 -𝑒 '( ! * " ): % (5)

Quantitative estimation of integral production yields

For an absolute estimation of the GSI production yields, spatial dependent correction factors for the material-specific g-ray attenuation, the solid angle coverage relative to the centre of the camera as well as the measured absolute detector efficiency in the centre of the scanner have to be applied.

Voxel-wise attenuation factor maps were created by a dedicated IDL routine for each phantom. The space around the phantom was filled with air with an attenuation factor of 𝜇 8%9 = 1.03677 ⋅ 10 ': [1/𝑐𝑚]. Inside the phantom the respective attenuation factor of GR, GEL or PMMA was applied (cf. Table 2). For the border voxels containing both phantom material and air, the resulting attenuation factor had to be calculated proportionally. Based on these maps, attenuation corrections for each line of response (LOR) of the scanner were generated, using a home-made software [START_REF] Pönisch | The modelling of positron emitter production and PET imaging during carbon ion therapy[END_REF]. These attenuation corrections are inputs to the filtered backprojection image reconstruction in order to generate an attenuation corrected PET image, which can then be further corrected for the limited solid angle detection. An additional PET image was reconstructed using filtered backprojection without the attenuation and the solid angle correction. Based on these images the combined correction factor C was calculated similar to [START_REF] Parodi | In-beam PET measurements of β + radioactivity induced by proton beams[END_REF] as:

𝐶 = ∫ 𝐴 ; 𝑑𝑟 FF⃗ ∫ 𝐴 𝑑𝑟 ⃗ ⋅ α ( 6
)
where Au is the uncorrected activity distribution, A the one corrected for attenuation and solid angle geometrical acceptance, both integrated over the whole volume where non-negligible activity is present. The detector efficiency a at the centre of the scanner was measured with a 𝑁𝑎 << point source of known activity as the ratio of true measured events over the number of radioactive decays, yielding a=0.022.

The number of radionuclides produced by a single primary beam particle can then be calculated from the production rate 𝑃 % and the number of beam particles delivered per second 𝐼 ! :

𝑁 % = 𝑃 % 𝐶 • 𝐼 ! (7)
2.4.4. Determination of depth dependent production rates

In order to obtain production rates that are spatially resolved along the irradiation axis (depth profiles), the phantoms were subdivided into several slices perpendicular to the beam direction, socalled regions of interest (ROIs). The full analysis pipeline (2.4.1.-2.4.3.) was then applied to these ROIs separately. The thickness of the ROIs varied in size between 3 mm to 50 mm. The choice represents a tradeoff between limited signal statistics and the desired resolution, which should be higher in regions of particular interest like the edge before the distal fall-off of the production yield.

Therefore, smaller ROIs were chosen in this part, contrary to coarser ones in the entrance channel.

Time-dependent analysis of the reconstructed PET images

In-beam data acquired at GSI during and after irradiation were reconstructed in several time frames for hypothetical in-room and offline measurement scenarios. For the in-beam scenario only those data acquired during the irradiation pauses between the spills were considered, hence the time frame for this reconstruction began after the first and ended with the last spill. For the in-room scenario the activity was reconstructed for a time frame of 3 minutes, starting 2 minutes after the irradiation. The offline scenario started at the timepoint where the ratios of the dominating b +emitting radionuclides matched the ratios of the offline HIT experiments and comprised the full PET acquisition time (cf. Table 3). It turned out that this timepoint was 2 to 3 minutes after completion of beam delivery, whereas realistic delays for patient examinations are in the order of 6 to 12 minutes.

Therefore, an additional late-offline scenario being closer to clinical practice was defined starting 10 minutes after the end of irradiation.

Data analysis of the HIT experiment

Offline activity measurements at HIT were performed with a commercial PET/CT scanner, providing integrated state-of-the-art image reconstruction functionality. The analysis pipeline was therefore more straightforward, as reported in [START_REF] Bauer | Offline imaging of positron emitters induced by therapeutic helium, carbon and oxygen ion beams with a full-ring PET/CT scanner: experiments in reference targets[END_REF](Bauer et al., , 2013a)). For image reconstruction, the full 30 min acquisition time was split into 30 single frames of one minute each. Volumetric activity images were obtained from the scanner software, using an iterative ordered subset expectation maximization reconstruction (3 iterations, 24 subsets, voxel size: 2.036 x 2.036 x 3.0 mm³) which automatically considered random, scatter and attenuation correction for activity quantification, and applying a Gaussian image smoothing with an isotropic kernel of σ = 5 mm. These parameters were chosen according to the results of (Bauer et al., 2013a), who tested extensively the impact of several image reconstruction parameters on the measured activity signal. The CT-based attenuation factors determined by the scanner software were checked and found to agree well with the literature value for PMMA and gelatine. The graphite attenuation factor, however, mismatched and was therefore corrected by a manipulation of the CT as reported in (Bauer et al., 2013a). Resulting PET images were exported from the scanner database for further external processing.

Data regression for the determination of the radionuclide-specific activity contribution 𝐴 #,% =>-present at the beginning of each PET measurement was conducted with an in-house fitting program based on python and implemented in MeVisLab® (MeVisLab Medical Solutions AG and Frauenhofer MEVIS)

analogous to equation (1):

𝐴(𝑟 ⃗, 𝑡) = / 𝐴 #,% =>-(𝑟 ⃗) ⋅ % 𝑒 '( ! ⋅* (8)
with the spatially resolved 𝐴 #,% =>-(𝑟 ⃗) as regression parameter, thus using tabulated half-lives for the considered radionuclides 𝑖 = C // , N /? , O /@ . The activity values were already quantitative at this point, including corrections for random coincidences, photon scattering and attenuation as well as measured detection efficiency. Consideration of the experiment-specific transfer time Δ𝑡 allowed to calculate the activity contribution of radionuclide species 𝑖 at the end of beam delivery:

𝐴 #,% (𝑟 ⃗, 𝑡) = 𝐴 #,% =>-(𝑟 ⃗) ⋅ 9𝑒 '( ! ⋅A* : '/ (9)
Absolute production yield maps were determined according to the GSI analysis pipeline part 2.4.2 using the spatially resolved 𝐴 #,% images as input and voxel-level application of equation ( 7) with C=1.

Integral production yields were obtained from an integration of the full image, depth resolved distributions from lateral integration of the activity for each voxel position in depth.

Results

Comparison of integral production yields

Integral production yields obtained from the PET activation measurements at HIT and GSI together with their uncertainties are reported in Table 5. Uncertainties on the HIT total yields are conservatively estimated as 10% for all radionuclides, considering uncertainties in the calibration of the beam monitor system, thus the number of delivered primaries, the choice of reconstruction time points, the accuracy of the manual time delay measurement and fit uncertainties [START_REF] Bauer | Offline imaging of positron emitters induced by therapeutic helium, carbon and oxygen ion beams with a full-ring PET/CT scanner: experiments in reference targets[END_REF](Bauer et al., , 2013a)).

Uncertainties on the GSI yields were explicitly estimated for each radionuclide and experiment by quadratically summarizing uncertainties from different regression methods (cf. Table 4), variations in the time structure of the irradiation and the number of delivered protons (2%-40%, depending on the experiment), and accuracy of the derived correction factor C (cf. Table 2). Regarding the impact of the different beam nozzle WEPLs and the PMMA plate in the GSI beam line, the conducted MC simulations predicted a small systematic effect of ~3% higher yields expected at HIT.

Ratios were calculated as the GSI yields divided by the corresponding HIT yields, with uncertainties obtained from maximum error estimate. GSI yields were found to be similar to or higher as the HIT yields, with the strength of the deviation depending on the phantom material.

Comparison of depth-resolved production yields

Depth-resolved production yields, stratified by radionuclide type are reported in Figure 2 In general, a reasonably good quantitative agreement was observed between the GSI and HIT data within uncertainties, in particular regarding the aforementioned unconsidered scatter background.

Larger deviations for 13 N as compared to the other radionuclide yields are mainly due to the scarce abundance of that radionuclide, resulting in higher uncertainties in the production yield quantification. The observed shorter activity range of the GSI data compared to the HIT profiles corresponds to the expected range shift due to the PMMA slab of 4.6 mm water-equivalent path length, which was accidentally present in the beamline during irradiation at GSI.

Despite the quantitative compatibility, the production yield curves of the GSI samples show a slightly convex shape whereas the HIT production yields rise continuously towards the maximum. A check of the impact of the PET image reconstruction, using filtered back-projection as well as an iterative reconstruction method for the GSI data, confirmed that the observed discrepancies are not due to the different reconstruction methods. The feature stems most likely from an insufficient correction of the limited-angle imaging artefacts of the double-head PET camera used at GSI.

Comparison of PET imaging scenarios

Figure 4 shows activity depth profiles reconstructed for the low-energy GSI experiments in the different monitoring scenarios in-beam, in-room, offline and late-offline. The profiles were normalized in the proximal entrance region to the same activity level. Since these concepts start measuring at different times during or after the irradiation the radionuclide fractions contributing to the resulting activity curves are different and determine the shape of the activity depth profile in the respective material.

Activity profiles for the graphite and the PMMA phantoms shown in Figure 4 (a) and (b) are rather stable over time due to the predominance of 11 C with a half-life of about 20 min. Differences observed in graphite are due to the 10 C contribution in the in-beam scenario, which is not present in the later post-irradiation imaging frames. With a rather flat cross section and a higher production threshold energy for 10 C of 34.5 MeV (compared to 20.61 MeV for 11 C, cf. Table 1) the in-beam activity depth-profile shows a more convex plateau and a slightly shorter range compared to the other scenarios. Similar effects were observed in PMMA, though less pronounced since the measured signal is composed of several radionuclides. As expected, activity depth-profiles measured in gelatine showed the largest differences over the four imaging frames. While the in-beam signal is largely dominated by the decay of 15 O, its fraction decreases with time and the shape of the depthprofile is more and more determined by the signal from the 11 C decay. Since the curve shapes of the energy-dependent production cross-sections are similar for 15 O and 11 C, the shape of the activity profile remains stable over time. Considerable changes were observed, however, in the activity range, defined as the position of the 50% activity fall-off. The largest difference found between inbeam and late-offline was 3.5 mm in case of the GEL-GSI-E1-1 experiment (cf. Figure 4 (c)).

Differences in activity range are due to different production threshold energies, as reported in Table 1. In gelatine, 3), 30min-Dt), late offline (10 min after beam delivery, 20 min).

In addition to variations in the profile shape, the chosen imaging scenario also impacts the total number of measured coincidences, Nc, which are reported in Table 6. This number depends on the duration of the measurement as well as on the half-lives of the radionuclides. For the graphite phantoms, for example, the number of late-offline coincidences is more than doubled compared to the in-beam scenario. The reason is the long acquisition time of the late-offline scenario and the long half-life of 11 C whose decay can still be detected after several minutes. In the gelatine phantoms the signal is dominated by 15 O decays with a half-life of about two minutes. Thus, the late-offline coincidence number is smaller than the coincidence number of the earlier in-beam and in-room scenarios.

For PMMA, late-offline and in-beam coincidence numbers are almost the same because the longerliving radionuclide 11 C is still detected in the longer lasting late-offline scenario, while short-living radionuclides like 10 C and 15 O are detected in the shorter lasting in-beam scenario.

Table 46 also reports the coincidence rates Rc [1/s] that were obtained on average over the respective monitoring scenarios. Generally, the average count rate decreases with time delay to measurement start. For the graphite phantoms, with a high number of long-living 11 C radionuclides produced, the difference in coincidence rates between in-beam and late-offline is small with a ratio of about 1.5. For the gelatine phantoms, on the contrary, the number of in-beam coincidences is about a factor of 10 to 17 times higher than in the late-offline scenario. This is due to a transition in the radionuclide predominance with time.

Early scenarios are dominated by the short-lived 15 O, with a high production yield, while the longlived 11 C radionuclides, generated in 16 O target fragmentation reactions with a much scarcer production yield, dominate the late-offline scenario. For the PMMA phantoms the in-beam coincidence rate is about three times higher than the late-offline coincidence rate. This ratio reflects the mixed radionuclide yield obtained in this composed material. The ratio is therefore larger than for graphite, due to the presence of 15 O and also due to the smaller amount of the longer-living 11 C radionuclides in PMMA compared with graphite. Compared with the ratio in gelatine, on the other hand, the PMMA coincidence ratio is smaller because the abundance of 11 C radionuclides in gelatine is negligible.

Discussion

A comprehensive data analysis has been performed to calculate radionuclide-specific production yields from in-beam PET data acquired at GSI. While methodological uncertainties in the calculation of relative quantities could be well controlled, the determination of reliable quantitative numbers remained difficult due to the research character of the infrastructure. The prototype PET camera installed at GSI is a limited-angle double-head system with well-known limitations regarding a reliable activity quantification [START_REF] Crespo | On the detector arrangement for in-beam PET for hadron therapy monitoring[END_REF]. Moreover, the beam delivered by the GSI research accelerator suffered from intensity fluctuations and empty spills for some experiments. This resulted in additional uncertainties regarding the determination of the actual irradiation time and number of delivered protons and thus of the eventual production rate. Overall, a quantification of the production yield was possible, as integral number and spatially resolved along the beam path, but with non-negligible uncertainties.

In a second step, these results were compared to radionuclide-specific production yield data determined at HIT, thus providing the first inter-facility comparison of long-lived proton-induced PET radionuclide yields in terms of both integral figures and depth-resolved distributions. GSI production yields were found to be similar to or higher than HIT yields, with the strength of the deviation depending on the phantom material. This phenomenon has also been observed in a recent production yield comparison conducted for similar phantom activation experiments using heavy ion beams [START_REF] Bauer | Offline imaging of positron emitters induced by therapeutic helium, carbon and oxygen ion beams with a full-ring PET/CT scanner: experiments in reference targets[END_REF]. One presumable reason for the overestimation of production yields determined at GSI is an incorrect handling of scattered coincidences. A Monte-Carlo based in-silico study has recently been performed by Baumgartl [START_REF] Baumgartl | PET-based Hadrontherapy Monitoring: Monte Carlo Simulations for Nuclear Interaction Studies in Phantoms[END_REF] investigating the fraction of scattered coincidences (SC) in the here presented GSI proton experiments, as well as the overall signal quantification bias due to SC after full PET image reconstruction. The author determined scatter fractions of 38.6%, 45.7% and 26.3% in the GSI PMMA, graphite and gelatine phantoms, respectively, for the different GSI experiments. This SC fraction leads to an eventual scatter bias of the reconstructed signal of 7.7% (PMMA), 10.6% (graphite) and 6.6% (gelatin). More details on the study are provided in the supplementary material. These signal quantification uncertainties are largely consistent with the discrepancies observed between the HIT and GSI production yields. They have not been accounted for in the present study, since the analysis pipeline was chosen to be as consistent as possible to previous works. However, the found scatter bias should be considered in potential future quantitative analyses of in-beam activation data acquired with the GSI PET camera.

The overall agreement between GSI and HIT data within uncertainties, in particular for PMMA and gelatine, is insofar remarkable as the imaging system used for activity detection and the analysis pipeline to determine production yields were fundamentally different from each other. GSI data were acquired with a prototype double-head PET camera and yields were calculated from raw count rates. HIT data, on the other hand, were acquired with a commercial full-ring PET scanner providing clinical accuracy and yields were calculated from spatially reconstructed activity maps.

Despite the applied signal corrections for the limited-angle FOV, qualitative differences remained in the depth-profile shape between HIT and GSI data, which could not be attributed to the chosen reconstruction algorithm.

Regarding the monitoring scenario, offline systems provide highest imaging performance since standard full-ring scanners can be utilized. However, the time delay of this workflow causes a considerable signal reduction and loss of spatial tissue activation information, resulting in a low clinical verification power. As we also demonstrated by comparing different hypothetical monitoring scenarios, the signal loss due to the physical activity decay during the time delay can be compensated by considerably prolonged acquisition time frames. However, in the clinical workflow this often conflicts with patient comfort and patient throughput, and the issue of signal distortion due to activity washout from the target volume remains unsolved. Thus, the current trend in clinical PET-imaging installations for beam range verification is towards in-beam solutions [START_REF] Bisogni | INSIDE in-beam positron emission tomography system for particle range monitoring in hadrontherapy[END_REF]. These systems typically require an open-ring geometry to ensure compatibility with the rotating beam delivery and robotic patient positioning system, thus have to accept a limited-angle signal detection with a compromised imaging quality. The quantitative and qualitative discrepancies in production yields and depth-profile shape, respectively, that we observed between HIT and GSI data underline the necessity of a thorough calibration of the imaging system. In particular geometrical image corrections for prototype devices have to be carefully considered to avoid artificial qualitative signal distortions that might impact clinical beam range verification.

In general, verification workflows involving a comparison of the measured signal to a pre-calculated expectation have a high demand on the calculation accuracy. It has been demonstrated in the past that available Monte Carlo based calculation engines used for PET-based beam range verification have to be at least benchmarked against and usually further tuned with reliable experimental data (Bauer et al., 2013a;[START_REF] Parodi | In-beam PET measurements of β + radioactivity induced by proton beams[END_REF][START_REF] Pinto | A filtering approach for PET and PG predictions in a proton treatment planning system[END_REF][START_REF] Seravalli | Monte Carlo calculations of positron emitter yields in proton radiotherapy[END_REF]. The availability of reliable production yield data, not only integral but also depth-resolved, is therefore of high interest for this specific application. However, our findings also demonstrate that production yields calculated from quantitative activation data can be very sensitive to the imaging and data analysis workflow.

Thus, a determination of energy dependent production cross sections from such thick target experiments can be problematic. The conduction of activation measurements for the tuning of the in-vivo activity signal prediction therefore remains an essential component of system commissioning for the integration of PET-based beam range verification in the clinic.

Conclusion

We provide experimentally determined production yield data in thick homogenous phantoms acquired at two treatment facilities with different PET imaging installations and processed by different analysis pipelines. The observed quantitative and qualitative differences between the resulting datasets underline the importance of a thorough system commissioning prior to clinical use to ensure a reliable beam range verification. In-house phantom activation measurements can thus be considered a gold-standard for careful system response characterization preceding the integration of the PET-based range verification workflow into clinical routine.

for each experiment, which allowed to correlate the time dependent annihilation signal to the corresponding spills (i.e. beam-on times). Combining the time correlation with the information of the signal origin it was possible to investigate the impact of scatter coincidences for particular phantoms and irradiation setups.

The simulated scatter fraction determined after image reconstruction, Asf, was defined as:

𝐴 !G = 𝐴 !H 𝐴 *I* ( 10
)
where Asc is the number of scattered coincidences and Atot the sum of scattered and true coincidences (random coincidences were neglected). The found scatter fractions after image reconstruction for a PMMA (9x9x20cm 3 ), graphite (9x9x15cm 3 ) and gelatine (9x7.3x20cm 3 ) phantom were 38.6%, 45.7% and 26.3%, respectively, assuming the default energy window of the GSI PET camera (E=[250-850] keV).

In order to estimate the scatter bias on the signal quantification, AS,bias, two different PET images were reconstructed from the simulated b + emitter maps. First, an image was reconstructed according to the experimental data analysis workflow: the photons penetrated the respective phantom material and were thus corrected for attenuation, both true and scattered photons by the same factor. The total amount of activity in the resulting reconstructed PET image is denoted as AS. In a second step, the image was reconstructed in absence of the phantom, thus the annihilation photons were propagated in air neglecting the attenuation correction. The total amount of activity in this second image, containing only unscattered events, is denoted as ANoM.

The activity bias due to scattered coincidences was then determined by comparing the amount of activity in the two reconstructed images:

𝐴 J,K%8! = 𝐴 ∑ -𝐴 1IL 𝐴 1IL (11)
The scatter bias found for the aforementioned observed scatter fractions Asf were 7.7%, 10.6% and 6.6% in the PMMA, graphite and gelatine phantom, respectively.

Figure 1 :

 1 Figure 1: Schematic drawing of the irradiation place and the in-beam PET measurement at GSI (left) and of the offline PET/CT scanner at HIT (right).

  ( 11 C, 15 O) and Figure 3 ( 13 N, 10 C), with each radionuclide-specific figure containing all relevant experiments, respectively. The 15 O depth curves from different targets are normalized to the respective 16 O density in gelatine and PMMA and should thus ideally be at the same level. Such normalization could not be performed for 11 C and 13 N due to multiple production channels (from 12 C and 16 O) in the case of 11 C or large signal fluctuations in the case of 13 N. The curves for 11 C in graphite are by definition expected to be ideally at the same level.

Figure 2 :

 2 Figure 2: Production yields measured after phantom irradiation at GSI (triangles) and HIT (squares) in gelatine (GEL), graphite (GR) and polymethyl methacrylate (PMMA). Top: 11 C production yield measured in GEL, GR and PMMA after irradiation with E1 (a) and E2 (b). Bottom: 15 O production yield measured in GEL and PMMA after irradiation with E1 (c) and E2 (d). The x axis shows the phantom depth and is displayed in water equivalent path length in [mm], the y axis reports the number of radionuclides produced per proton per [mm] penetration depth. For the 15 O production yield the y axis is normalized to the 16 O density of the respective phantom material.

Figure 3 :

 3 Figure 3: Production yields measured after phantom irradiation at GSI (triangles) and HIT (squares) in gelatine (GEL), graphite (GR) and polymethyl methacrylate (PMMA). Top: 10 C production yield measured at GSI in GR and PMMA after irradiation with E1 (a) and E2 (b). Bottom: 13 N production yield measured in GEL and PMMA after irradiation with E1 (c) and E2 (d). The x axis shows the phantom depth and is displayed in water equivalent path length in [mm], the y axis reports the number of radionuclides produced per proton per [mm] penetration depth. For the 10 C production yield the y axis is normalized to the 12 C density of the respective phantom material.

  11 C and 15 O radionuclides are generated in proton-induced fragmentation reactions of the target oxygen nuclei only. The threshold energies are 27.50 MeV and 16.79 MeV, respectively, thus explaining the observed range decrease over time.

Figure 4 :

 4 Figure 4: GSI activity profiles as a function of depth measured in the graphite (a), PMMA (b) and gelatine (c) phantoms after mono-energetic pencil beam irradiation with E1=125.4 MeV/u for the different hypothetical in-room, offline and late-offline scenarios with the following time structure (acquisition start / acquisition time): in-beam: (end of first spill, duration of beam

Table 1 :

 1 Production channels for the main radionuclides produced in graphite, PMMA and gelatine, corresponding production threshold energies Ethr and half-lives T1/2[START_REF] Beebe-Wang | Dependence of the Production Yields of Positron Emitters in Proton Therapyon the Cross Section Data Variations[END_REF].

	Radionuclide	Production channel	Ethr [MeV]	T1/2 [s]
	11 C	12 C(p,pn) 11 C	20.61	1,223
	11 C	16 O(p,3p3n) 11 C	27.50	1,223
	15 O	16 O(p,pn) 15 O	16.79	122
	13 N	16 O(p,2p2n) 13 N	5.66	568
	10 C	12 C(p,p2n) 10 C	34.5	19

Table 2 :

 2 Overview of the homogeneous phantoms used at GSI (upper section) and HIT (lower section). Reported are the

	C used for the quantitative production yield calculation together with their uncertainty s . Attenuation coefficients of the HIT
	phantoms were automatically determined during image reconstruction from a CT scan for the polymethyl methacrylate
	(PMMA) and the gelatine (GEL) phantom, and manually adjusted ('+corr', cf section 2.5) for the graphite phantom (GR).
	Experiment ID	W x H x L [cm]	µ [1/cm]	C [•10³] / s
	GR-GSI-E1	9 x 9 x 15	0.156	4.6 / 4%
	GEL-GSI-E1-1	9 x 9 x 20	0.096	7.0 / 3%
	GEL-GSI-E1-2	9 x 9 x 20	0.096	6.8 / 3%
	PMMA-GSI-E1	9 x 9 x 20	0.111	6.6 / 3%
	GR-GSI-E2	9 x 9 x 31	0.142	4.1 / 5%
	GEL-GSI-E2	9 x 9 x 30	0.096	5.8 / 3%
	PMMA-GSI-E2	9 x 9 x 30	0.111	5.4 / 4%
	GR-HIT-E1/-E2	7 x 7 x 20	CT-based+corr	-
	GEL-HIT-E1/-E2	10 x 5 x 35	CT-based	-
	PMMA-HIT-E1/-E2	9 x 9 x 30	CT-based	-

geometrical size

(width W, height H and length L)

and, for the GSI phantoms, the material-specific photon attenuation coefficients µ (at a photon energy of 511 keV) used for activity reconstruction, and the experiment-specific correction factors

Table 3

 3

	: Parameters of different phantom experiments: experiment ID, proton beam energy E, total number of delivered
	particles Np, average beam intensity I, and duration of irradiation tirr, as well as the time delay Dt between the end of
	irradiation and start of PET data acquisition for the HIT experiments and (*) the time delay used to mimic an "offline"
	scenario for the GSI experiments (cf. section 2.4.5), and the PET acquisition frame time tdecay.

Exp ID E [MeV/u] Np [•10 11 ] I [10 8 /s] tirr [min:s] Dt [min:s] tdecay [min] GR

	-GSI-E1	125.4	2.50	16.3	10:00	02:40*	30
	GEL-GSI-E1-1	125.4	1.05	6.6	10:00	03:03*	35
	GEL-GSI-E1-2	125.4	0.99	4.9	08:20	03:40*	25
	PMMA-GSI-E1 125.4	1.16	8.5	10:00	01:54*	30
	GR-GSI-E2	176.95	1.30	13.2	08:22	02:40*	30
	GEL-GSI-E2	176.95	1.41	76.6	08:24	02:03*	30
	PMMA-GSI-E2 176.95	1.46	13.8	08:19	02:04*	30
	GR-HIT-E1	125.67	1.12	32.0	03:20	04:02	30
	GEL-HIT-E1	125.67	1.12	32.0	04:35	04:35	30
	PMMA-HIT-E1 125.67	1.12	32.0	03:43	03:43	30
	GR-HIT-E2	176.75	1.12	32.0	04:11	04:11	30
	GEL-HIT-E2	176.75	1.12	32.0	03:42	03:42	30
	PMMA-HIT-E2 176.75	1.12	32.0	03:48	03:48	30

Table 4 :

 4 GSI radionuclide-resolved activity values at the end of irradiation A0 [10³/s] as obtained from the time analysis of the decay data for all considered phantoms and deviations D in[%] to the results from the fit cross-check.

	A0 • [10 3 /s]

Table 5 :

 5 Production yields 𝑁 - " with 𝑖 = 𝐶 .. , 𝑂 ./ , 𝐶 .0 , 𝑁 .1 obtained from GSI and HIT data and their respective sum 𝑁 233 ) . Ratios are calculated as GSI yields divided by HIT yields. All yields are provided in units of [⋅ 10 45 𝑝𝑒𝑟 𝑝𝑟𝑖𝑚𝑎𝑟𝑦 𝑝𝑟𝑜𝑡𝑜𝑛].Uncertainties s on the yields are provided in percentage, uncertainties on the ratio in absolute numbers.

	Experiment ID	𝑵 𝑪 𝟏𝟏 W / s	𝑵 𝑶 𝟏𝟓 W / s	𝑵 𝑪 𝟏𝟎 W / s	𝑵 𝑵 𝟏𝟑 W / s	𝑵 𝒂𝒍𝒍 X / s
	GR-GSI-E1	42,057 / 39%	-1,590 / 41%	-	43,647 / 41%
	GR-HIT-E1	36,435 / 10%	-	-	-	36,435 / 10%
	Ratio	1.1 ± 0.46	-	-	-	1.2 ± 0.51
	GR-GSI-E2	94,215 / 12%	-5,248 / 13%	-	99,463 / 13%
	GR-HIT-E2	57,540 / 10%	-	-	-	57,540 / 10%
	Ratio	1.6 ± 0.16	-	-	-	1.7 ± 0.17
	GEL-GSI-E1-1	6,332 / 8%	24,894 / 8%	-	2,581 / 8%	33,807 / 8%
	GEL-GSI-E1-2	7,490 / 27%	31,374 / 27%	-3,706 / 27%	42,570 / 28%
	GEL-HIT-E1	6,860 / 10%	23,733 / 10%	-2,412 / 10%	33,005 / 10%
	Ratio-1	0.9 ± 0.11	1.0 ± 0.13	-	1.0 ± 0.13	1.0 ± 0.13
	Ratio-2	1.0 ± 0.32	1.3 ± 0.38	-	1.5 ± 0.44	1.2 ± 0.37
	GEL-GSI-E2	10,879 / 8%	39,482 / 5%	-4,963 / 14%	55,324 / 16%
	GEL-HIT-E2	10,313 / 10%	38,369 / 10%	-3,868 / 10%	52,550 / 10%
	Ratio	1.0 ± 0.11	1.0 ± 0.11	-	1.2 ± 0.13	1.0 ± 0.11
	PMMA-GSI-E1	28,196 / 10%	10,929 / 10% 1,607 / 15% 1,085 / 10%	41,817 / 15%
	PMMA-HIT-E1	22,513 / 10%	8,101 / 10%	-1,451 / 10%	32,065 / 10%
	Ratio	1.2 ± 0.17	1.3 ± 0.19	-	0.7 ± 0.1	1.3 ± 0.18
	PMMA-GSI-E2	52,637 / 12%	18,620 / 12% 2,090 / 22% 1,835 / 12%	75,182 / 23%
	PMMA-HIT-E2	36,085 / 10%	14,030 / 10%	-1,564 / 10%	51,679 / 10%
	Ratio	1.4 ± 0.22	1.3 ± 0.20	-	1.1 ± 0.18	1.4 ± 0.22

Table 6 :

 6 Total number of coincidences Nc and corresponding coincidence rates Rc [1/s] averaged over the respective acquisition time for the hypothetical in-room, offline and late-offline scenarios of all GSI experiments with the following time structure (acquisition start / acquisition time): in-beam: (with first spill, duration of beam delivery), in-room: (2 min after beam delivery, 3 min), offline: (Dt after beam delivery when radionuclide ratio matches HIT experiments (cf. Table3), 30 min-Dt), late offline (10 min after beam delivery, 20 min).

	Exp ID	in-beam	in-room	offline	late-offline
		Nc	Rc [1/s]	Nc	Rc [1/s]	Nc	Rc [1/s]	Nc	Rc [1/s]
	GR-GSI-E1	9.8	22.6	4.7	26.2	29.9	17.1	19.4	14.8
	GR-GSI-E2	6.5	17.0	2.7	15.0	21.8	12.9	14.0	11.1
	GEL-GSI-E1-1	13.5	30.7	2.8	15.7	6.5	3.2	3.0	1.8
	GEL-GSI-E1-2	11.9	32.7	3.2	17.6	4.0	3.5	2.6	3.0
	GEL-GSI-E2	21.2	57.0	5.2	29.1	13.6	7.9	4.8	3.9
	PMMA-GSI-E1	7.6	17.8	2.7	14.9	12.8	7.7	7.0	5.9
	PMMA-GSI-E2	12.4	34.2	4.9	27.0	25.1	12.2	14.9	9.4

Acknowledgements

The authors would like to thank Dieter Schardt and Peter Heeg for their assistance with the experiments at GSI, and Stephan Brons for his support of experiments at HIT. JB acknowledges funding from the German Federal Ministry for Research and Education (grant agreement numbers 01IB08002F, 01IB13001G, 13GW0436A) and by the University of Heidelberg (Olympia-Morata fellowship 2021/2022).

Conflicts of Interest

The authors have no conflicts of interest to declare.

Supplement

Determination of scatter correction bias on production yields

Sophisticated simulation studies are required to realistically estimate the scatter fraction, which were not available at the time of previously published GSI data analyses [START_REF] Fiedler | Anwendung des in-beam PET Therapiemonitorings auf Präzisionsbestrahlungen mit Helium-Ionen[END_REF][START_REF] Parodi | In-beam PET measurements of β + radioactivity induced by proton beams[END_REF][START_REF] Sommerer | In-beam PET monitoring of mono-energetic 16 O and 12 C beams: experiments and FLUKA simulations for homogeneous targets[END_REF]. Such a dedicated Monte-Carlo (MC) study has been conducted by [START_REF] Baumgartl | PET-based Hadrontherapy Monitoring: Monte Carlo Simulations for Nuclear Interaction Studies in Phantoms[END_REF]. The author reproduced the full experimental GSI set-up using the GATE MC package (V6.2, based on Geant4 V9.5) [START_REF] Jan | GATE V6: a major enhancement of the GATE simulation platform enabling modelling of CT and radiotherapy[END_REF][START_REF] Jan | GATE: a simulation toolkit for PET and SPECT[END_REF]. The simulated annihilation signal was processed by the same data analysis and image reconstruction pipeline as the measured data. The simulation framework considered also the macroscopic time structure of the pulsed beam delivery


Jorge Ruiz-Garcia

Christos Bilitos

Enrica Martini

Giovanni Toso

Stefano Maci

David Gonzalez-Ovejero

Metal only Multi-beam Fabry-Perot Antenna

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Suat Gumussoy
email: suat.gumussoy@cs.kuleuven.be

Didier Henrion
email: henrion@laas.fr

Marc Millstone
email: millstone@cims.nyu.edu

Michael L Overton
email: overton@cims.nyu.edu

Multiobjective Robust Control with HIFOO 2.0

Keywords: robust control, multiobjective control, optimization

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Multiobjective control aims at designing a feedback control law meeting potentially conflicting specifications defined on various input/output channels.

In the context of linear systems, a standard approach to multiobjective control is the Lyapunov shaping paradigm proposed in the mid 1990s by [START_REF] Scherer | Multiobjective output feedback control via LMI optimization[END_REF], as an outgrowth of the LMI (linear matrix inequality) formalism of [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF]. Within this scope, multiobjective controller design boils down to semidefinite programming (linear programming over the cone of positive semidefinite matrices) provided all the closedloop specifications are certified simultaneously by a unique quadratic Lyapunov function. Moreover, the controller is retrieved a posteriori via tedious linear algebra, and its order is equal to the order of the open-loop plant plus the order of the weighting functions, which can be quite high in practice, in contradiction with simplicity of implementability requirements of embedded control laws. Another computational approach to multiobjective control exploits the parametrization of all stabilizing controllers described e.g. by [START_REF] Vidyasagar | Control system synthesis: a factorization approach[END_REF]. Linear programming can be used in this context to design controllers, see e.g. [START_REF] Boyd | Linear control design: limits of performance[END_REF], but they are typically of very high order.

Following a decade of research efforts, these restrictions have been gradually relaxed (distinct Lyapunov functions for distinct performance channels, parameter-dependent Lyapunov functions, decoupling between Lyapunov and controller variables, lower-order controller design) at the price of an increased computational burden. The ROMULOC (robust multiobjective control) toolbox is a recent public-domain matlab implementation of these techniques, see [START_REF] Peaucelle | Robust multiobjective control toolbox[END_REF].

Particular cases of multiobjective robust control problems include strong stabilization (where a plant must be stabilized by a controller which is itself stable, see [START_REF] Vidyasagar | Control system synthesis: a factorization approach[END_REF]) or simultaneous stabilization (where a single controller must be found that stabilizes several plants, see [START_REF] Blondel | Simultaneous stabilization of linear systems[END_REF]). Most of the algorithms or heuristics available to solve these problems also typically result in very high order controllers.

In this paper, we introduce the new release 2.0 of our freely available package hifoo, which is aimed at removing the above mentioned limitations in the context of multiobjective controller design. First, the controller order is fixed at the outset, allowing for lowerorder controller design. Second, no Lyapunov or lifting variables are introduced to deal with the conflicting specifications. The resulting optimization problem is formulated on the controller coefficients only, resulting in a typically small-dimensional, nonsmooth, nonconvex optimization problem that does not require the solution of any large convex subproblems, relieving the computational burden typical of Lyapunov LMI techniques. Because finding the global minimum of this optimization problem may be hard, we use an algorithm that searches only for local minima. While no guarantee can be made about the behaviour of this algorithm, in practice it is often possible to determine an acceptable controller quite efficiently.

See also [START_REF] Apkarian | Mixed H 2 /H ∞ control via nonsmooth optimization[END_REF] for nonsmooth nonconvex optimization techniques applied to multiobjective robust control. As far as we know, no software implementation of these techniques is publicly available at present.

Problem Formulation

The ith generalized plant

P i = (A i , B i , C i , D i ) describes the state-space equations ẋi (t) = A i x i (t) + B i 1 w i (t) + B i 2 u i (t), z i (t) = C i 1 x i (t) + D i 11 w i (t) + D i 12 u i (t), y i (t) = C i 2 x i (t) + D i 21 w i (t) + D i 22 u i (t),
where

A i ∈ R n i ×n i , D i 12 ∈ R p i 1 ×m 2 , D i 21 ∈ R p 2 ×m i 1
, with other matrices having compatible dimensions. The signals (z i , w i , y i , u i ) respectively represent the regulated outputs, the exogenous inputs (including disturbance and commands), the measured (or sensor) inputs, and the control inputs. Let N be the number of plants.

The problem is to choose a single controller

K = (A K , B K , C K , D K ) with state-space equations ẋK (t) = A K x K (t) + B K y(t), u(t) = C K x K (t) + D K y(t),
where A K ∈ R n K ×n K , with B K , C K , D K having dimensions that are compatible with A K and the generalized plant matrices. The controller order n K is fixed, so it can be specified by the designer. The H ∞ norm of the ith closed loop system is the norm of the transfer function from input w i to output z i ; see [START_REF] Zhou | Robust and optimal control[END_REF] for details. The complex stability radius is a useful alternative measure when no w i and z i performance channels are specified and hence the H ∞ norm is not defined: for a stable closed loop system, this is the largest 2-norm perturbation to the closed loop system that can be tolerated while guaranteeing that the perturbed system remains stable. The spectral abscissa of a closed loop system is the largest of the real parts of its poles (eigenvalues).

Let β j , j = 1, ..., N each be a real number or ∞, and consider the following optimization problem:

min K max j=1,...,N {g j (K) : β j = ∞} (1) subject to g j (K) ≤ β j , j = 1, 2, . . . , N, (2)
where each g j is one of the following supported functions of the closed-loop system for P j , abbreviated by a single letter as follows:

• 'h': H ∞ norm (∞ if unstable);
• 'r': reciprocal of complex stability radius (∞ if unstable);

• 's': spectral abscissa.

The functions g j are all nonconvex, nonsmooth functions of the controller matrices. Thus, the optimization problem is potentially quite difficult. We focus on two scenarios.

Scenario 1. All β j = ∞, and all g j are 's'. Thus the problem is to minimize the maximum of the real parts of the closed loop poles of all plants. This approach is suitable for simultaneous stabilization, as the goal is to move the poles as far left in the complex plane as possible (in the minmax sense). If the final objective value is negative, all closed loop plants are stable.

Scenario 2. All β j = ∞, and all g j are 'h'. The problem is to stabilize all plants and minimize the maximum of the H ∞ norms of the closed loop plants. The advantage over Scenario 1 is that the goal is to not only stabilize the plants, but also optimize their H ∞ performance (in the minmax sense). The disadvantage is that evaluating the H ∞ norm repeatedly is more time-consuming than computing the spectral abssissa repeatedly, and I/O performance channels must be specified.

These two scenarios are the ones for which benchmarking is done in Section 5. However, the interface to HIFOO is sufficiently flexible that there are many other alternative ways to call it which may be of interest to users. For example:

Alternative 3. All β j = ∞, and all g j are 'r'. The problem is to stabilize all plants and maximize the minimum of the complex stability radii of the closed loop plants (minimize the maximum of their reciprocals). The advantage is that the complex stability radius is a more robust measure of stability than the spectral abscissa, so this is an appropriate alternative to Scenario 2 when I/O performance channels are not specified.

Alternative 4. This is a more specific example. Suppose

β 1 = ∞, β 2 = 100, g 1 is 'h' and g 2 is 'r'.
The problem is to stabilize both plants and minimize the H ∞ norm of the first closed loop plant subject to the complex stability radius of the second being at least 0.01.

If one wishes to impose restrictions, such as stability, on the controller (so-called strong stabilization), it is possible to do so using the multiple-plant model above by defining a plant P i so that the closed loop plant is equivalent to the system described by the controller. However, we provide a more convenient way to specify controller stability directly, as described below. Benchmarks assessing the value of a previous version of hifoo for strong stabilization appear in Gumussoy et al. (2008).

3 Optimization method hifoo 2.0 searches for local minimizers of ( 1)-( 2). The algorithm has two phases. In each phase the main workhorse is the BFGS optimization algorithm, which is surprisingly effective for nonconvex, nonsmooth optimization, see [START_REF] Lewis | Nonsmooth optimization via BFGS[END_REF]. The user can provide an initial guess for the desired controller (see below); if this is not provided, hifoo generates randomly generated initial controllers, and even when an initial guess is provided, hifoo generates some additional randomly generated initial controllers in case they provide better results.

The first phase is stabilization: BFGS is used to minimize the maximum of the spectral abscissae of the closed loop plants for which g j is either 'h' or 'r'. This process terminates as soon as a controller is found that stabilizes these plants, thus providing a starting point for which the objective function for the second phase is finite.

The second phase is optimization: BFGS is used to look for a local minimizer of the following unconstrained problem:

min K F (K) + ρ N j=1 max(0, g j (K) -β j ),
where F is the objective function defined in (1) and ρ is a positive penalty parameter multiplying the sum of the constraint violations. If BFGS is unable to find a point for which the constraint violations are zero, the penalty parameter ρ is increased and the optimization is repeated as needed (unless F is identically zero, that is all β j are finite).

Although there are no guarantees, very often this process is quite effective and reasonably fast. By default, hifoo invokes the gradient sampling method of [START_REF] Burke | A robust gradient sampling algorithm for nonsmooth, nonconvex optimization[END_REF] after BFGS terminates, but this is generally more time consuming, and can be avoided as explained in the next section.

4 User interface hifoo is written in matlab. The simple call K = hifoo(P, order) looks for a controller K solving Scenario 2 above: stabilize the plants described by P and minimize the sum of the H ∞ norms of the closed loop plants, using a controller of the specified order. Here P is a cell array of plants specified in any of several formats, typically using the ss class of the matlab Control System Toolbox. A useful abbreviation is P{j}='K', which specifies that the closed loop system for the jth plant is the controller itself. If the order argument is omitted, the default order 0 is used (static output feedback).

A more general calling sequence is [K, F, viol] = hifoo(P, order, init, fun, ... upperbnd, options)

where P and order are as above, init is an initial guess for the controller (several formats are supported), fun is a string specifying the characters defining the supported functions g j (see Section 2) or a single character if all g j are the same (default: 'h'), upperbnd is an array specifying the upper bounds β j (default: all bounds set to ∞), and options is a structure with various optional fields, some of which are described below. The order of the input arguments does not matter except that P must be first. The outputs are, in addition to the controller K, the value of the objective function F (K) in (1) and a vector of constraint violations max(0, g j (K)β j ). There is a fourth output argument loc ("local optimality certificate") that can be requested if gradient sampling is used in addition to BFGS.

Some of the more useful fields in options are:

• options.cpumax: requests hifoo to quit when the CPU time in seconds exceeds this quantity (default: ∞)

• options.fast: 1 to use a fast optimization method only (BFGS), 0 to finish optimization with a slower method (gradient sampling, which may give a better answer) (default: 0, as long as quadprog is in the path; see below)

• options.prtlevel: one of 0 (no printing), 1 (minimal printing), 2 or 3 (more verbose) (default: 1)

• options.struct: specifies sparsity structure to be imposed on the controller (see documentation for details)

• options.weightNormK: weight for adding a penalty on the size of the controller to the objective function, specifically (

A K 2 + B K 2 + C K 2 + D K 2 ) 1/2 (default: 0)
• options.augmentHinf: weight for adding the reciprocal of the complex stability radius to the H ∞ norm to avoid closed loop plants that are only marginally stabilized: applies to all plants for which g j is 'h' (default: 0) hifoo 2.0 uses the following external software:

• required: hanso 1.0, a hybrid algorithm for nonsmooth optimization, freely available from the hifoo web page;

• required: the matlab Control System Toolbox, for H ∞ norm and complex stability radius computation; this also provides user-friendly system modeling with the ss class;

• optional: quadprog from either the matlab Optimization Toolbox or mosek; needed only by the gradient sampling part of the algorithm, which is not required.

Benchmarks

We consider a number of simultaneous stabilization problems from the literature; in each case we wish to stabilize multiple plants with a single controller. We consider both Scenario 1 (optimizing the spectral abscissa, by a call such as K = hifoo(P,'s',n), where n is the order of the controller), and Scenario 2 (optimizing H ∞ performance, by a call such as K = hifoo(P,'h',n), or equivalently K = hifoo(P,n)). Notice that we do not need to explicitly set the upper bounds to ∞, since that is the default value.

The following matlab script illustrates the process, using Scenario 1. There are three first-order plants, given in [START_REF] Jia | Condition and algorithm for simultaneous stabilization of linear plants[END_REF], and we wish to stabilize them with a first-order controller: Note that a third-order controller was designed in [START_REF] Jia | Condition and algorithm for simultaneous stabilization of linear plants[END_REF].

>> P =
The output of hifoo may differ on different runs since the initialization is done randomly. However, the output from one run may be used to initialize a second run on the same problem; the result cannot be worse. For example, following the run above with >> K = hifoo(P,K,'s',1); results in ... hifoo: best order 1 controller found has spectral abscissa -0.286524 which is a slight improvement. Depending on the problem and the initial randomization, several successive calls to hifoo may be required to obtain a stabilizing controller.

We have collected various academic and application examples. We now give benchmark results for Scenarios 1 and 2.

Scenario 1: Simultaneous Stabilization

To evaluate hifoo, we consider 31 benchmark problems for simultaneous stabilization as shown in Table 1. There are 11 problems from applications and 20 academic test problems. The benchmark problems include dynamic, state-feedback and static-output feedback controllers.

The performance of hifoo is very good considering the large variety and number of benchmark examples. hifoo successfully solves the simultaneous stabilization problems with a low-order single controller compared to the existing methods in the literature.

Note that contrary to the other methods in the literature, hifoo allows the user to set the controller order a priori.

In particular, hifoo is very successful in application problems. hifoo performs better than existing methods for the application benchmarks Lane-Keeping of Automated Vehicles and Helicopter Toy. The Rationalized Helicopter Model (RHM14) shows the conservativeness of some of the methods in the literature, which produce a 26th order stabilizing controller. In contrast, hifoo stabilizes the same benchmark by a static controller.

hifoo shows similar performance for academic test problems. For almost all problems, the results are better than or equivalent to the results of existing methods. Some of these are Bredemann Ex.5.6, Chen et al., F.-Anaya et al. Ex.1 and Ex.2, and Gündeş-Kabuli, for all of which the existing methods stabilize the benchmark problem with a high-order controller and hifoo solves the same problem by a static or first-order controller.

There are two benchmark examples for which hifoo performs slightly worse than the existing methods, namely Ex.1-1 and 2-1 from Bhattacharyya et al. which are interval plants. hifoo stabilizes both benchmark problems with a controller whose order is one more than the order of the controllers in the literature.

Scenario 2: Simultaneous H ∞ Optimization

We start with two academic benchmark problems for simultaneous H ∞ optimization as shown in Table 2. In the first problem, there are two plants of order two, with optimal fullorder H ∞ performance (using different second-order controllers) equal to 1.290 and 1.245 respectively. hifoo finds a single first-order controller with the same H ∞ performance (as measured by the maximum of the norms of the two closed loop plants). Also, hifoo finds a static controller with performance 1.530. The previously best known H ∞ performance using a single controller was 1.806, using a fourth-order controller.

The second problem consists of three plants of order two with respective optimal H ∞ norms 1.290, 1.245 and 1.038 using full-order controllers. In the literature a 6th order simultaneously stabilizing controller is known, with H ∞ performance 1.833. hifoo achieves the optimal H ∞ performance with a single first-order controller. where the weighting functions are W 1 (s) = 1 s+1 and W 2 (s) = 0.2. The dimensions of the plant P i are A i ∈ R n i +1×n i +1 , D i 12 ∈ R 2p 2 ×m 2 , D i 21 ∈ R p 2 ×m 2 .

Table 3

 3 extends the benchmarks for simultaneous stabilization in Table1to H ∞ performance optimization. Performance channels are added to each plant in Table1using the augw function in the Control System Toolbox by

	G{k} = augw(P{k},tf(1,[1 1]),[],0.2);

The research of D. Henrion was partly supported by project MSM6840770038 of the Ministry of Education of the Czech Republic. The work of M. Millstone and M.L. Overton was partly supported by the U.S. National Science Foundation under grant DMS-0714321; views expressed in the paper are those of the authors and not of the NSF. The work of M.L. Overton was also partly funded by Université Paul Sabatier, Toulouse, France.

Table 1: Benchmarks on simultaneous stabilization: 11 industrial and 20 academic examples. The third and fourth columns show the lowest stabilizing controller orders in the literature and the lowest found by hifoo.

Our benchmark results are given in Table 1, showing the problem name, dimensions (with n max = max{n i }), the lowest known controller order from the literature (third column), and the controller orders obtained by hifoo (final column). For each controller order, hifoo was run 10 times; the best results are reported. Runs for order k with k > 1 were initialized with the best controller found for order k -1. We used options.fast = 0 (the default) but all runs were repeated with options.fast = 1, for which the results were almost the same. Further information on references for the benchmark problems and other methods will be provided in a technical report that will be made available at this website.

In conclusion, hifoo, by allowing the designer to specify the controller order, is a useful tool for simultaneous stabilization and simultaneous H ∞ optimization.