Soft x-ray tomograms are consistent with the magneto-hydrodynamic equilibrium in the Wendelstein 7-X stellarator

Soft x-ray tomograms are inferred from experimental data obtained during the recent operational phases of the superconducting, optimized stellarator Wendelstein 7-X. It is shown that the reconstructed soft x-ray emission profiles of the plasma are consistent with the numerically calculated magneto-hydrodynamic equilibrium of Wendelstein 7-X. In order to obtain reliable tomograms, the full chain of electrical and geometrical influences on the x-ray observation has to be taken into account. This has been achieved by formulation and application of an extended forward model. The forward model has been verified using phantom data derived from surrogate tomograms.


Introduction
In order to assess the plasma performance in the superconducting stellarator Wendelstein 7-X (W7-X) sophisticated diagnostics are necessary. The international stellarator scaling ISS04 [1] predicts energy confinement scaling as τ E ∝ n 0.54 , which means that the triple product [2] nT i τ E scales favorably with the density as n 1.54 . Hence, high density operation in a stellarator is desirable and possible since the stellarator is 1 For a list of members of the W7-X team refer to T Klinger et al 2019 Nucl.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. not limited by the Greenwald density [3]. A high-temperature plasma emits radiation in a broad range of the electromagnetic spectrum. Contributions include bremsstrahlung, Coulomb collisions (electron-electron, electron-ion) and recombination radiation, where high-Z impurities (e.g. Fe, W) significantly contribute with line radiation in the soft x-ray (SXR) range [2]. Detection of the SXR emission of the plasma can therefore lead to insights on the shape of the plasma and on the profiles of plasma temperature and density as well as impurity density. A typical realization of a SXR tomography diagnostic utilizes a set of pinhole cameras with beryllium filters and photodiodes as detectors [4]. Such a camera array allows for the tomographic reconstruction of the SXR emission profile of the plasma.
In the recent operational phases of W7-X, first measurements using the SXR multi camera tomography system (XMCTS) have been performed [5]. The XMCTS consists of 20 pinhole cameras poloidally arranged at the triangular cross-section of W7-X. Initial data evaluation showed unexpected asymmetries in the measured raw signals and tomographic reconstruction methods (that performed well on phantom data) resulted in tomogram contours which strongly deviate from the numerically calculated magnetic flux surfaces of the plasma equilibrium. It turned out that the main issue are the mechanical shutters of the diagnostic, which partially blocked the lines-of-sight. Since the original data analysis approach did not take into account the effect of the shutter blocking, large systematic errors caused inconsistencies in the tomographic inversion process. A suitable forward model based on a Bayesian formulation to include the shutter positions was developed. The application of the Bayesian principle of Occam's razor [6] resolves the undesired effects caused by errors in the forward model (systematic errors) and errors in the measurement (statistical errors). Also, the small changes in the camera positions, caused by deformations of the plasma vessel, have been included in the forward model. Tomographic reconstruction typically is an ill-posed problem, since the continuous emission profile is sampled only along a finite number of lines-of-sight. The number of unknowns (in this case local values of the emission) is much larger than the number of measurements. Several types of approaches have been developed in the past to deal with this ill-posedness. The first type expands the emission profile into orthogonal basis functions and fits their parameters to the experimental data. When the number of free parameters of the basis functions is less or equal to the number of data points, the problem becomes wellposed again. The Cormack method uses Zernike polynomials for the radial and Fourier series for the poloidal direction [7]. It has been developed further to use Bessel functions for the radial direction [8]. For the methods to follow the plasma is discretized into a finite number of pixels, small enough that constant emission over the size of a pixel can be assumed. There are usually many more pixels than lines-of-sight, leading to an underdetermined problem which requires some kind of regularization to be well-posed. This discretized type of tomographic reconstruction methods is based on some kind of regularization to constrain an underdetermined problem. Examples are the Maximum-Entropy method or the Tikhonov regularization with the Minimum-Fisher information regularization operator [9][10][11]. In the present paper, for the tomographic reconstruction of the SXR emission, the method of Gaussian process tomography (GPT) [12] is used for the tomographic reconstruction of the SXR emission. It can be considered as a third step in the evolution of tomography models with a fully probabilistic treatment while retaining analyticity for the nonlinear components and linearity for the numerical components of its formulation. This makes it a computationally efficient method, which seamlessly integrates into the Bayesian modeling approach. The effective dimensionality of the emission profile is explicitly controlled by use of an appropriate kernel function to construct the covariance matrix of a multivariate normal probability distribution for the emission profile. The GPT method has already been applied to SXR tomography [13][14][15][16][17]. This article is organized as follows: The XMCTS diagnostic is described in section 2. The tomographic reconstruction method is discussed in section 3. First successful modeling results from experimental data of the XMCTS diagnostic on W7-X are presented in section 4. The findings are discussed in section 5. The details of the forward model for pinhole-type cameras with obstacles in the way can be found in appendix A.

Experimental setup
The SXR tomography diagnostic XMCTS is installed inside the plasma vessel of W7-X [18]. A schematic diagram of its mechanical design and location in the device is shown in figure 1. The XMCTS features 20 cameras, organized in four segments (1-4) with five cameras each (A-E). Each of these cameras utilizes an AXUV-22EL SXR silicon photodiode array, manufactured by IRD Inc. [19]. Each array comprises of 22 photodiodes in a linear arrangement with an active area of 1 mm × 4 mm per diode. 18 of the 22 photodiodes per detector array are designated to receive radiation from the plasma, resulting in a total of 360 photodiode signals [5]. Background levels of the photodiode signals are measured before plasma ignition and are subtracted from the measured signals that are acquired during the plasma discharge in order to remove offsets in the measured signals. Between the aperture slit and the photodiode array a filter made of a beryllium foil is placed that blocks the low-energy electromagnetic spectrum and is transparent for the high-energy x-ray radiation from the plasma (details in [5]). In order to protect the beryllium foils in the cameras against deposition during plasma vessel bakeout, pneumatically operated shutters have been installed to block the aperture opening. The shutters were opened at the beginning of an experimental day and closed at the end of it. Throughout the day, a gate valve with extremely low leakage rate shuts off the Bourdon spring volume from the pressurized air supply, so that pressure fluctuations in the air supply cannot influence the shutters. It is therefore reasonable to expect that the air pressure inside the Bourdon spring is constant. Other sources of mechanical forces on the shutters are not considered relevant. The shutters are also required for calibration purposes during long pulses (>100 s) using the integrated illumination system [5]. In figure 2, one of the XMCTS cameras is shown. The shutter (in black) is clearly visible. The front side of one of the XMCTS cameras. The black sheet metal piece (color from coating to avoid sticking) is the shutter arm. The blue arrow indicates the direction of shutter movement. Inside its opening, the aperture of the pinhole-type camera is visible. In the background, the Bourdon spring operating the shutter can be spotted (details in [20]).

Methods
The GPT method is an application of the Bayesian Occam's razor principle to the tomographic reconstruction problem [12]. The emission profile is modeled as a n-dimensional GP. Typically, n = 1 for an emission profile which is constant on flux surfaces, n = 2 for a free-form solution in the poloidal plane, n = 3 takes into account toroidal variations, and n = 4 is applied to model time-dependent three-dimensional emission profiles. The kernel functions of the GPs typically feature parameters (usually called hyperparameters to denote their indirect influence on the result) which have to be adjusted in order to prevent under-or overfitting of the experimental data. The power of this method lies in the fact that, when using GPT, an objective criterion is available to guide the algorithm towards a choice of hyperparameter values that explains the data as good as possible, without requiring a too complex emission profile. The so-called evidence term in Bayes' formula (see equation (18) below) acts as a penalty for small length-scales of the inferred emission profile to limit the space of solutions to those just as complex as allowed for by the uncertainties of the observations (i.e. the measured photodiode signals).

Forward model
The forward model allows one to calculate the predictions (here the photodiode currents in the XMCTS cameras) from the free parameters (here the emission profile of the plasma, see figure 3). The measured photodiode currents are assumed to be proportional to the SXR emission of the plasma. The linear XMCTS forward model is parameterized to account for diagnostic effects, e.g. the shutters partially blocking the linesof-sight. The linear forward model is represented as a (N P × N G ) matrix M(H d ) parameterized by the diagnostic-specific hyperparameters H d : Equation (1) where s c is the scaling factor of the c-th camera. The forward model matrix M can be written as: where the matrix C (contribution matrix [9] or geometric matrix [21]) contains the geometric information about the XMCTS diagnostic 2 . The (N G × N P ) contribution matrix C is constructed using the algorithm described in appendix A.2, which is implemented as a computational node in the Minerva framework [22]. Its (k, l) component represents the contribution of a unit emission in the k-th pixel of the emission profile to the l-th photodiode current. Several other forward modeling tools have been developed in previous publications. The basis for the analytical part used in VolRayTomo has been described in [23]. An analytical two-dimensional forward model for a pinhole-type camera for an emission profile specified in projection space can be found in [24]. The assumption of an axisymmetric plasma allows to take computational advantage of this symmetry in raytracing approaches to build the forward model [25,26]. The hyperparameters of the emission profile are denoted H ϵ (see below) and the full set of hyperparameters in this model is H = H d ∪ H ϵ . All hyperparameters are assumed to have uniform prior probability distributions. This choice allows one to specify reasonable upper and lower bounds but otherwise leaves the model indifferent with respect to their values.
A schematic graphical representation of the forward model is depicted in figure 3. One of the 20 cameras (namely 4B) was not working properly and has been excluded from the analysis. The forward model is formulated including this missing camera for clarity and applicability for future experiments. The code part providing the interface between measured data (from the experiment archive) is called 'datasource'. The conversion from digitized voltages to photocurrents, including calibration for transimpedance gains, is done inside the XMCTS datasource node. The datasource node also provides the diagnostic geometry, e.g. the toroidal angle and the locations and extents of the photodiode detectors and the apertures.

Representation of the SXR emission profile
For the forward model, the plasma is approximated as a slab for this first assessment of the shutter influence. This is a reasonable simplification because of the large aspect ratio of W7-X and the flip-symmetry of the triangular poloidal cross section. The shutters only influence the viewing cones of the detectors in poloidal direction. Therefore, a numerical implementation of the forward model is only necessary in the R-Z-plane. The toroidal contribution of the viewing cones is handled analytically based on equivalent model problems in two and three dimensions without shutters (see appendix A.1.3).
The grid of the emission profile is located in the R-Z-plane at the toroidal location of the XMCTS. The pixel centres are located in regular intervals on a rectangular grid: where r min and r max are the minimum and maximum R position of the emission grid pixels, z min and z max are the minimum and maximum Z position of the emission grid pixels and N R and N Z denote the number of emission grid pixels in R and Z direction. Only emission grid pixels inside the poloidal cross section of the plasma vessel at the toroidal location of the XMCTS (G for grid) are included in the tomographic reconstruction: The pixel locations in G get a one-dimensional index number k to simplify the formulation of the model as: where the symbols i(k) and j(k) denote the mapping between the index k and the indices i and j in the R and Z directions, respectively, and N G is the total number of emission grid pixels inside the poloidal cross section. The emission profile is then described by: Two different representations of the SXR emission profile have been employed in this work. For the shutter configuration inference procedure a one-dimensional artificial emission profile is used (the SXR radiation is assumed constant on flux surfaces); a full two-dimensional emission profile is chosen to represent the actual tomographic reconstruction. They mainly differ in their complexity, i.e. their number of free parameters and the diversity of emission profile shapes. This allows one to impose quite strong constraints on the shape of the emission profile as long as the shutter configuration is only coarsely known. Once the shutter configuration can be estimated sufficiently well, more complex emission profile shapes can be considered.
In the following the procedure for inferring the shutter configuration is described. Since the plasma parameters on magnetic flux surfaces are in zeroth-order approximation constant, the SXR emission is assumed constant on flux surfaces. The flux surface geometry is obtained from magneto-hydrodynamic equilibria computed using the VMEC code [27]. The corresponding VMEC coordinates (s, u, v) for each of the pixel centre locations (r i(k) , z j(k) ) are computed using an iterative coordinate inversion routine [28]. The radial location s (not to be confused with the vector of scaling factors per camera s) is defined as the normalized toroidal magnetic flux Ψ enclosed in the corresponding flux surface, i.e. s = Ψ/Ψ LCFS , where Ψ LCFS is the total toroidal magnetic flux enclosed by the last closed flux surface (LCFS). A onedimensional radial profile of the SXR emission is then mapped onto the flux surfaces, as illustrated in figure 4. For the model profile function a function is created that is simple while representing typical features of plasma parameter profiles, e.g. a maximum in the center and a decay to zero at the edge. The simplified emission profile is described as: where ε 0 is the emission at s = 0, a is a steepness parameter and b is the radial location of the inflection point of the profile function. This emission profile parameterization has no hyperparameters. The one-dimensional radial profile of the SXR emission is mapped onto the three-dimensional space by employing the VMEC flux surface geometry. The inferred radial coordinate s(r i(k) , z j(k) ) is then used as the parameter for the emission profile function (10) to define the (uniform and isotropic) SXR emission at the given grid pixel: Alternatively, the emission profile can be represented as a two-dimensional Gaussian process in the R-Z-plane. The emission profile ϵ is then modeled by a N G -dimensional multivariate normal distribution N , which acts as the Bayesian prior where µ ϵ is the mean vector and Σ ϵ is the covariance matrix. The mean vector µ ϵ generally can be nonzero, but is commonly (and also here) taken to be zero. The two-dimensional squared-exponential kernel function is chosen to define the prior covariance matrix. The value of k SE (r, z, r ′ , z ′ ) describes the covariance between two points (r, z) and (r ′ , z ′ ) in the two-dimensional emission profile. This choice of the prior covariance matrix thereby makes the prior probability distribution being a Gaussian process: The hyperparameters σ R and σ Z correspond to length-scales in the R and Z directions, respectively, over which the emission profile may vary by values up to σ f .

Model for the XMCTS diagnostic
The measured photocurrents are denoted by: The covariance matrix of the photodiode currents is with the individual uncertainties of the (uncorrelated) signals σ l . The signal uncertainties are estimated based on the measured signals, where the photocurrents are multiplied by a relative error estimate and added to an estimated absolute error. The relative and the absolute errors are estimated by the standard deviation of the noise on the measured signals during the time interval over which the mean of the measured signals is computed. The Bayesian likelihood of a set of measured photocurrents I meas given the emission profile ϵ is modeled by a N P -dimensional multivariate normal distribution:

Inference of the shutter configuration
The partial shadowing of the viewing cones of the XMCTS diagnostic by the shutters is sketched in figure 5. The parallel distance d || to the aperture center and the perpendicular distance d ⊥ to the aperture plane parameterizes the shutter position of each camera. Tilting of the shutter is not considered here, since the tilting angle of the shutter is on the order of <2 • . Note that d || and d ⊥ can vary between experiments and must be considered in the analysis of the respective data. In this sketch neither the detector nor the viewing cone opening angle are discretized for multiple origins of viewing cones on the detector surface and multiple partial viewing cones, as is done in the actual analysis to enhance the accuracy of the model. The shutter positions and distances and the scaling factors of each camera introduce a total of 60 hyperparameters, which can only be determined based on the consistency requirement of the inferred emission profile. This requires orders of magnitude more iterations in any reconstruction method. A forward model is therefore favorable which can be integrated into the inversion algorithm.
The algorithm developed to infer the shutter configuration along with the parameters of the emission profile is as follows. Two preferably different emission profiles at two different time instants in the same plasma discharge are modeled to provide two 'views' on the unknown shutter configuration to the optimizer. The shutter configuration is initialized as centered positions where the shutters do not block (parts of) the viewing cones. The intensity scaling factors s in equation (2) are initialized to unity. Then, an interactive procedure is started, which alternates between two main steps: (a) A conjugate gradient optimizer is used to infer the values of the parameters ε 0 , a and b of the emission profiles and the scaling factors s, which maximize the likelihood of the given (measured or simulated) photodiode currents. The optimizer performs steps in the parameter space until the absolute improvement in the likelihood per step is reduced to numerical accuracy (here 2.2 × 10 −16 ). (b) The shutter configuration of each camera is scanned over a 0.1 mm grid of parallel and perpendicular shutter positions. This is necessary in the first step to avoid getting stuck in local maxima of the posterior probability distribution. Additionally, this allows to cache a database of contribution matrices (one for each combination of parallel and perpendicular shutter positions) to save time during the inference process. For each individual camera, the shutter configuration leading to the highest likelihood for the corresponding photodiode currents is retained.
The shutter configuration of the XMCTS model is thereby iteratively adjusted to be compatible with the one-dimensional radial emission profile (10) up to the 0.1 mm resolution of the grid parallel and perpendicular positions of the shutters.

Tomographic reconstruction method
In the framework of GPT [12], the posterior probability for an emission profile ϵ given the measured signals I meas can be directly computed by application of Bayes' formula: where the likelihood p(I meas | ϵ) is given by (16), the prior is given by (12) and the evidence p(I meas ) is obtained by marginalization: which can be evaluated analytically for the given form of the likelihood and prior terms. The posterior probability distribution of the linear model with fixed hyperparameters is again a multivariate normal distribution with posterior meanε and posterior covariance Σε written as: The posterior mean is also the maximum-posterior (MAP) estimate of the emission profile given the fixed hyperparameters, which implicitly enter through the forward model matrix M [6]. The tomographic reconstruction is then performed by evaluating equations (19) and (20) [12,29]. Optimization of the hyperparameter values (in the sense of Occam's razor) towards the maximum of the posterior probability is then performed using the Hooke-and-Jeeves optimization algorithm [30] to infer the most comprehensive emission profile that is still justified by the match between the measurements and the predicted signals in lieu with the modeling uncertainties.

Results
The newly developed volumetric raytracing algorithm is applied to SXR emission profiles from experimental data measured by the XMCTS diagnostic on W7-X. A comparison is made between the method where the shutter influence and the scaling factors of each camera are purposefully ignored and the full forward model as described in the previous sections. Prior to the application to experimental data, the algorithm is tested with pre-defined phantom data.

Reconstruction of surrogate tomograms
In the following the procedure of the shutter configuration inference is demonstrated using artificial modeled data for the shutter configuration as well as for the emission profile. The steps are in short: • The parameters of the two modeled emission profiles, the scaling factors of each camera and the shutter configuration parameters were sampled from their respective prior probability distributions and set to the sampled values. • The photodiode currents were then predicted using the simplified forward model described in section 3.1. Artificial Gaussian noise of the same level as observed in experiments [5] was added to the signals to mimic experimental conditions and to prevent overfitting. • The initial shutter configuration was chosen as d || = 0 mm and d ⊥ = 0.15 mm for all cameras and the initial scaling factor of each camera was set to 1. Tomographic reconstruction with simultaneous inference of the shutter configuration was performed using the method outlined in section 3.2.
The result, namely the shutter configuration, the scaling factors of each camera and the inferred emission profiles, is shown in figure 6. There is a good match between the targeted emission profiles and the final state of the reconstruction. Also the parallel positions of the shutters d || are well inferred by our algorithm. The perpendicular distances of the shutters to the camera's apertures d ⊥ are not as closely reconstructed, but the overall match is acceptable. For shutters blocking parts of the viewing cones, the model is more sensitive to the parallel positions of the shutters than to the perpendicular distances. This is reflected in the more accurate reconstruction of the parallel positions than the perpendicular distances and also in the scaling factors of each camera, which can not be reconstructed to an exact match, since differences remain in the shutter  configuration. The uncertainty of the inferred shutter configuration is of the order of the spacing between the grid points of the shutter configuration (here 0.1 mm). Note that this can be further refined if necessary and the sole reason for scanning the shutter configuration over a pre-determined grid was to be able to cache the contribution matrices for speeding up this computation.

Tomograms for available emission models
The algorithm has been benchmarked on various W7-X experimental data. Here, exemplarily one tomogram is inferred for a particular W7-X plasma discharge (experiment program 20180918.45 at 4.0 s ± 20 ms). The photodiode currents have been averaged over the specified time interval. The shutter configuration was inferred using the simplified emission profile model and the algorithm outlined in section 3.2. In order to initially assess the influence of the shutters on the reconstructed emission profile, the Gaussian process model for the emission profile (as presented in section 3.1.1) was used with the shutter influence purposefully ignored. Note that using this advanced model for the emission profile it is no longer assumed that the emission is constant on flux surfaces. The result is shown in figure 7. The reconstructed emission pattern is not compatible with the assumption of constant emission on flux surfaces. Radially aligned features are present and a significant fraction of the emission profile shows negative emission. There is however a reasonable match between the predicted and the measured photodiode currents.
The inference method is applied to obtain the shutter configuration and the scaling factors for each of the cameras under the assumption of constant emission on flux surfaces. These values for the hyperparameters are subsequently used in the Gaussian process forward model for the emission profile. These results are shown in figure 8. The inferred profile clearly shows almost constant emission on the flux surfaces. The photodiode currents represent the measurements much better than in figure 7. The length-scale hyperparameters of the Gaussian process are larger by a factor ≈3, which indicates a smoother profile. Tiny artifacts of negative emission remain. This will be subject to further optimizations in the future and is likely to disappear once the shutter configuration has been inferred to a higher degree of precision.
No other diagnostics data currently supports the strong deviations from emission being constant on a flux surface as seen in figure 7. Also note that the magnetic configuration of W7-X has been optimized for a small Shafranov shift at high values of volume averaged normalized plasma pressure of up to ⟨β⟩ ≲ 5% [31]. The outward shift of the emission maximum in figure 7 is much larger than could be expected from the MHD equilibrium (as indicated by the flux surfaces shown).

Discussion and conclusion
Incomplete knowledge on the exact geometrical properties of a physical system is a common challenge in various areas of physics. In the present case of tomographic reconstruction of a plasma equilibrium in a stellarator fusion device, the proposed method of volumetric raytracing is a novel solution to deal in particular with not well-defined obstacles in the viewing cones of a pinhole camera array. The gain in speed outweights the restricting assumptions about constant emission profiles along the coordinate axis perpendicular to the reconstruction plane. Once the values of the hyperparameters that describe the unknown geometry have been inferred, other approaches, e.g. full three-dimensional Monte-Carlo raytracing methods, can be used to refine the result to even higher degrees of precision. Volumetric raytracing fills the gap between modeling tools that are too simple to include obstacles in the light path on the one hand and on the other hand methods that are too slow to reasonably allow for standard tomographic reconstruction including the reconstruction of unknown obstacles. Our method to infer the shutter configuration and the scaling factors of each camera utilizes a conjugate gradient optimizer. We note here that initial attempts to perform the (heavily nonlinear) hyperparameter optimization using the Hooke-and-Jeeves optimization algorithm failed. However, its robustness and the fact that, in contrast to conjugate-gradient, it works on the function values of the posterior probability only makes the Hooke-and-Jeeves optimizer the preferred algorithm for the optimization of the Gaussian process hyperparameters once the shutter configuration has been fixed. The squared-exponential kernel for the covariance matrix of the Gaussian process was selected since it is one of the most simple kernel functions and since it is desirable to introduce as few hyperparameters as possible in the Gaussian process.
Using our method, proper tomographic reconstructions of the plasma equilibrium become possible despite the unknown positions of the shutters leading to large systematic errors without correction. The shutter configuration is inferred based on data from two different timestamps in the same experimental program assuming that the shutters do not move throughout an experiment. We are able to consistently infer the shutter configuration with this assumption and thus conclude that the shutters do not move over the time span of an experiment in W7-X. The assumption of constant emission on flux surfaces needed to be included only for inference of the shutter configuration and reconstruction of the actual emission profile was performed using a free-form two-dimensional Gaussian process. The flux-surface-mapped arctan emission profile model represents a plasma which has its maximum emission near the magnetic axis and the emission decays in some way towards the plasma edge. This approximation is also valid when magnetic islands are present in the core. The presence of magnetic islands in the actual experimental plasma would thus only lead to a worse initial guess compared to the case without islands. Since the shutter configuration is optimized together with the emission profile in the GPT routine, the final outcome should not be affected by presence of internal islands. Constant emission on flux surfaces implies that the inferred emission profile also constrains the location and shape of the flux surfaces and thus the actual magnetic topology of the plasma equilibrium. It is left for further work to exploit this to constrain the flux surface shape in a reconstruction of the MHD equilibrium itself. It should be noted that the assumption of constant emission on flux surfaces is only true to a certain degree (especially in a stellarator), as symmetry-breaking error fields and accumulation of intrinsic impurities can easily modify the existence of flux surfaces and the poloidal impurity distribution, respectively.

Acknowledgments
This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix A. Modeling the influence of the shutters of XMCTS
The main focus of this section is the presentation of the newly developed VolRayTomo forward modeling code for pinhole cameras with an arbitrary number of obstacles, which might partially or completely block the viewing cone(s) of the detector(s). The VolRayTomo code and the data and figures presented in this article are available from the corresponding author upon reasonable request.

A.1. Analytical model problems
Model problems in two and three dimensions are presented and solved analytically. These serve both to verify the numerical approach presented in section A.2 and to re-scale the numerical solution in two dimensions into a solution applicable to the three-dimensional geometry of XMCTS. The model problems are formulated for a situation where a detector receives emission through a single aperture from an emission region which is infinitely extended in the parallel direction with respect to the aperture. This setup allows for several simplifications in the mathematical formulation of the forward model and makes it analytically traceable.
A.1.1. Two-dimensional forward model The basic geometry of the problem is shown in figure A1. Inside some emission volume, a constant isotropic emission ε is present. Part of the emission can arrive at a detector through an aperture. The purpose of the forward model is to compute the integral power present at the detector surface given the geometric parameters of the setup. The discretization used for the following derivation is shown in figure A2. An infinitesimal small detector region dx sees a part of the emission volume with a size of dL dξ through an infinitesimal small part of the aperture du. The size of the emission region for a set of (du, dL) is given by: The radiation intensity emitted from this region is then given by: Here, ε(L) is the emission with ε = 1 inside the emission region and ε = 0 otherwise. Isotropic emission is assumed and therefore, a local radiation intensity of: is present at the location of the detector. The factor (2πL) −1 accounts for the decay of the local radiation intensity with increased distance to an isotropically emitting source. The effective area of the detector, which is given by projecting the detector into a plane perpendicular to the line-of-sight, is given by cos(α) dx, where α is the angle of incidence. The contribution to the integral power at the detector surface from this set of (dx, du, dL) is then given by: The integral power P 2d at the detector is given by integration over the detector surface, the aperture width and the distance along the line-of-sight: (25) where the detector extent is [x 0 , x 1 ] and the aperture extent is [u 0 , u 1 ]. The width of the infinitesimal emission region dξ is given by the intercept theorem: Inserting this into equation (25) leads to the following expression: The line integral along the line-of-sight can be carried out directly and results in: where e w is the thickness of the emission layer along the y coordinate. Therefore, the total integral power at the detector surface is given by: An analytical solution to this integral exists: The integral power at the detector is given by: where the detector extent is [x 0 , x 1 ] and the aperture extent is [u 0 , u 1 ].

A.1.2. Three-dimensional forward model
The basic geometry of the three-dimensional tomography model problem is shown in figure A3. Note that the angle α depends now on two directions, i.e. du and dw. The size of the infinitesimal part of the emissive volume dL dξ dζ seen through an infinitesimal part of the aperture cos(α) du dw (projected into the plane perpendicular to the central line of sight), is given by: where α is the angle of incidence. The expressions for dξ/du and dζ/dw can be found by application of the intercept theorem: The projected area of an infinitesimal part of the detector dx dz is given by: The total integral power at the detector surface is found analogously to equation (32): where the integral along L was carried out as already shown in equation (31). The corresponding antiderivative was found using Mathematica [32]: The total integral power at the detector surface can now be computed by evaluating this analytical solution for the given geometry similar to equation (35): A. 1

.3. Scaling a two-dimensional solution to three dimensions
A numerical solution in three dimensions including obstacles by e.g. raytracing would require too much computational effort for inferring the shutter geometry of the XMCTS. The numerical solution to a tomography problem including obstacles in two dimensions can be applied to a three-dimensional system by re-scaling it according to the ratio of the analytical solutions to the corresponding model problems in two and three dimensions, in which obstacles are neglected. It is assumed that the emission is constant along the third coordinate which is analytically handled by this re-scaling procedure. This is the case for XMCTS due to the large aspect ratio of the W7-X plasma. The scaling factor for each of the individual detector signals is given by the analytical solutions from equations (35) and (42) for the corresponding detector and aperture geometries:

A.2. Numerical solution in two dimensions: VolRayTomo
The presence of obstacles blocking parts of the lines-of-sight of the XMCTS made it necessary to use a numerical method to obtain the forward model. Here, the volumetric raytracing approach is implemented. Figure A4. Obstacle coordinate system as a rotated two-dimensional Cartesian coordinate system. The red and blue dots indicate the edges k and l of the obstacle, which are specified by the tuple (r, φ, k, l).

A.2.1. Computation of the contribution matrix in VolRayTomo
The volumetric raytracing algorithm is assembled from the individual algorithms described in the sections A.2.2-A.2.6. All specified obstacles are merged and transformed into the detector coordinate systems according to equation (49). The widest possible viewing cones compatible with all obstacles that are specified for a given detector are determined. The extent of the detectors is specified using coordinates as introduced in figure A4 for the special case of r = 0 and φ = 0; the active detector area is then [k d , l d ] with subscripts 'd' for 'detector'. The overlap between the active detector area and the widest possible viewing cone through all associated obstacles is determined.
where the superscript 'a' means 'active'. The active detector area is evenly subdivided into n d subdivisions in order to perform a midpoint integration over evaluation points l a i : Each of the midpoints along the active part of the detector surface represents the origin for a viewing cone. For each point on each detector, the widest possible viewing cone compatible with all obstacles specified for the corresponding detector is determined. Each of the widest possible viewing cones is evenly subdivided into n a partial viewing cones (with subscript 'a' for 'aperture') to account for non-uniformity of the respective intersection polygons with the emission grid pixels. For each of the partial viewing cones, it is checked by angle comparisons for each emission grid pixel if any of the corners of the emission grid pixel is located within the viewing cone. If this is the case, the intersection polygon between the partial viewing cone and the emission grid pixel is determined. The area A of the intersection polygon determines the amount of emission from the respective emission grid pixel contributing to the respective integral power at the detector surface. The contribution from the emission of the k-th emission grid pixel to the signal of the l-th photodiode is given in the entry (C) l,k of the contribution matrix C: where c l,k i,j is the contribution from the jth partial viewing cone originating at the ith midpoint on the active detector surface to the (l, k)th contribution matrix element. The elements c l,k i,j are computed as follows: • The intersection area A between the (i, j)th viewing cone and the k-th emission grid pixel and (if there is an intersection at all) the centroid of the intersection polygon (c x , c y ) are computed according to the algorithm described in appendix A.2.5. The following quantities are only computed if there is an intersection between the viewing cone and the emission grid pixel. • The angle of incidence φ i of the emission originating at the center of the intersection polygon with respect to the orientation of the active detector surface element with width ∆l a is computed according to equation (62). • The distanceL between the centroid of the intersection polygon (c x , c y ) and the midpoint on the active detector surface (x, y) (obtained from the corresponding l a i using coordinate transforms) is computed from: (46) • The scaling factor to re-scale the numerical tomography results for a two-dimensional system into a result applicable to a three-dimensional device as XMCTS is computed from equation (43) for the l-th photodiode and corresponding geometry of the camera.
These intermediate quantities are combined to obtain c l,k i,j : • camera group (one segment of five cameras of XMCTS; four in total), • camera (five per XMCTS segment; 20 in total), • detector group (one AXUV photodiode array per camera) and Figure A5. Two local coordinate systems are nested within each other. The origin of the nested system (x, y) is given by (x ′ 0 , y ′ 0 ) and it is rotated by an angle φ ′ 0 with respect to the outer system (x ′ , y ′ ). Note that for the geometry shown k > 0, l < 0, k ′ < 0 and l ′ < 0.
Each coordinate system is specified by its origin (R 0 , Z 0 ) and a rotation angle φ around this origin in the poloidal plane. This definition allows to easily account for possible misalignments of e.g. the camera segments as a whole around their point of fixture by defining the origin of the camera group to be the point of fixture of the segment and rotating it around that point using the rotation angle. The coordinate system of the obstacles is a rotated two-dimensional Cartesian coordinate system, which is illustrated in figure A4. An obstacle is defines as a straight line with two significant points located on it (named k and l after the respective distances to the zero of the coordinate direction along the obstacle line). In between these two points the obstacle is transparent and outside of them, it is opaque.

A.2.3. Coordinate transforms and merging of obstacles
An obstacle specified in two nested coordinate systems is shown in figure A5. A coordinate system as defined in figure A4 is described within another coordinate system (x ′ , y ′ ). An obstacle (r, φ, k, l) can be expressed using the transformed coordinates (r ′ , φ ′ , k ′ , l ′ ) in the outer coordinate system: The inverse transform from an obstacle given in outer coordinates (r ′ , φ ′ , k ′ , l ′ ) into the nested coordinate system (r, φ, k, l) is given as follows: At any level of the four nested coordinate systems, obstacles can be specified in the local coordinates. This facilitates easy specification of obstacles, which are relevant for more than one detector, since any obstacle specified in a given coordinate system is taken into account for all detectors in the coordinate systems nested within the given one. One detector group coordinate system (for the one AXUV22EL diode array per camera) and inside it, 18 individual detector coordinate systems for the respective 18 photodiodes are nested within each camera coordinate system. In case of the XMCTS, the shutters are specified as obstacles in the respective local camera coordinate systems, since there is one shutter installed per camera. Each shutter is therefore taken into account for all detectors nested in the respective camera coordinate system. The volumetric raytracing is carried out locally in each individual detector coordinate system. Therefore, the coordinates of the obstacles specified in the global coordinate system have to be transformed into the respective camera group coordinate systems. The appropriately transformed global obstacles are then concatenated to the list of obstacles directly specified in the respective camera group coordinates. This step is referred to in the following as 'merging' the obstacles. The merged camera group obstacles are then transformed into each of the nested camera coordinate systems and merged with the obstacles directly specified in the respective camera coordinates. Analogously, obstacles directly specified in detector group coordinates and in detector coordinates are merged and transformed until finally, all  obstacles which are specified in any of the parent coordinate systems of each detector have been transformed into the respective detector coordinates. This process is sketched out in figure A6.
A.2.4. Volumetric raytracing through obstacles A sketch of the obstacle geometry assumed in this algorithm is shown in figure A7. A viewing cone is defined by its starting point at the active detector surface and the direction angles of its two edge rays, which are traced through the obstacles. If there is only one obstacle specified (e.g. the aperture of the XMCTS camera), the computation of the viewing cones can be carried out directly using the method described in section A.1.1.
Often times, this is not the case and one has to find the possible viewing cones (i.e. viewing cones undisturbed up to the emission region), which are compatible with the given detector geometry and all obstacles relevant to this detector. The goal is to compute the ray crossing through the k edge of one obstacle and the l edge of another obstacle, which passes through all other obstacles and the ray crossing through the l edge of one obstacle and the k edge of another obstacle, which also passes through all other obstacles. There are three nested loops, which iterate over the list of obstacles associated with the given individual detector. The first loop unconditionally loops over all these obstacles. The second loop loops over all obstacles which follow in the list after the current one of the first loop. The third loop loops over all obstacles which are not the first Figure A8. Basic geometry for plane intersection formula. In this context, a plane can be an obstacle or a detector surface. A ray (red arrow) from (r 1 , φ 1 , l 1 ) through (r 2 , φ 2 , l 2 ) hits a plane given by (r 3 , φ 3 ) at l 3 (blue arrow Here, a, b and c are the loop counters and nObs is the number of obstacles associated with the detector. Using the plane intersection formula (53), the intersection of a ray originating at k of the first loop's current obstacle and the l edge of the current obstacle in the second loop with the plane of the current obstacle in the third loop is computed. If this intersection point on the obstacle plane of the third loop is located between its k and l, the ray is considered as passing and the loops are exited. Otherwise, the loops continue to search for a passing ray. Using the plane intersection formula, the intersection of a ray originating at l of the current obstacle in the first loop and the k of the current obstacle in the second loop with the plane of the current obstacle in the third loop is computed next. If this intersection point on the obstacle plane of the third loop lies between its k and l, the ray is taken as passing and the loops are exited. Otherwise, the loops continue to search for a passing ray. Finding two passing rays in the two tests is is equivalent to having found a finite-width viewing cone passing through the set of obstacles.
The intersection of each of the rays with the detector plane is computed if two passing rays are found. This intersection region is coerced to the available detector surface area and subdivided to carry out a midpoint integration over the remaining (possibly full) detector surface area extent. Now that the detector range possibly receiving emission has been determined, the same algorithm is used to find the possible viewing cones for each of the midpoint integration points along the detector extent possibly receiving emission. Two rays are traced from each midpoint on the active detector surface. The first one continues through the k edge of all obstacles one after another and the intersection is checked with all following obstacles. The second ray continues through the l edge of all obstacles one after another and the intersection is checked with all following obstacles. The plane intersection basic geometry is shown in figure A8. The obstacle and detector planes are specified using a coordinate representation which is similar to the obstacle geometry specification illustrated in figure A4. It is overdetermined for this use case but allows to omit additional coordinate transformations. A ray r(t) for t > 0 that originates at (r 1 , φ 1 , l 1 ) and proceeds through (r 2 , φ 2 , l 2 ) is defined by the following parameterization: r(t) = r 1 cos(φ 1 ) − l 1 sin(φ 1 ) r 1 sin(φ 1 ) + l 1 cos(φ 1 ) + t r 2 cos(φ 2 ) − l 2 sin(φ 2 ) − r 1 cos(φ 1 ) + l 1 sin(φ 1 ) r 2 sin(φ 2 ) + l 2 cos(φ 2 ) − r 1 sin(φ 1 ) − l 1 cos(φ 1 ) .
The point in plane 3 s(l 3 ) is given by the following expression: The goal now is to solve r(t) = s(l 3 ) for l 3 . This is a system of two coupled equations of the form: with the coefficients: a 1 = r 1 cos(φ 1 ) − l 1 sin(φ 1 ) − r 3 cos(φ 3 ) b 1 = r 1 sin(φ 1 ) + l 1 cos(φ 1 ) − r 3 sin(φ 3 ) Figure A10. Basic geometry for the rectangle-ray-intersection algorithm in two dimensions. The numbers 1, 3, 5 and 7 denote the corners of the emission grid pixel. The black numbers 0, 2, 4 and 6 denote the sides of the emission grid pixel. The colored numbers inside the circles are denoted (1), (2), (3) and (4) in the text and refer to the intersection points between the rays 1 (red) and 2 (green) and the emission grid pixel. The centroid c of the intersection polygon (gray) between the rectangle and the viewing cone defined by the two rays is at a distance r from the point on the detector (x, y). The arrows along the rays indicate the order in which intersections with the rectangle are searched for.
the rectangle if it is located on one of the edges of the rectangle, which has to be kept in mind when interpreting the results. In order to find all corners in the opening angle of the viewing cone, based on intersections with rectangle faces, an algorithm which can find all corners between any two given faces of a rectangle was formulated and implemented. If all corners of the rectangle are within the angle range of the viewing cone, then the rectangle is fully embedded in the viewing cone and contributes completely to the detector signal. If no intersection between any of the rays and the rectangle faces or corners occurs, the rectangle is considered to lie outside of the viewing cone and does not contribute at all to the detector signal. For the partially covered area, we can use the area of a (not necessarily convex for the formula to be applicable) polygon (x 0 , y 0 ), (x 1 , y 1 ), ..., (x n−1 , y n−1 ), which is given by: where x n = x 0 and y n = y 0 is assumed, i.e. the polygon is closed (cf theorem 1.3.3 in [33]). The centroid of the polygon (c x , c y ) can be calculated with similar formulas: (y i + y i+1 )(x i y i+1 − x i+1 y i ).
(61) Figure A11. Chosen model geometry to illustrate the convergence of the two-dimensional numerical solution using volumetric raytracing towards the analytical prediction for the two-dimensional model problem. A detector receives radiation from an infinitely-extended emissive layer through an aperture. The dotted part of the aperture corresponds to the pinhole.
It should be noted that the term (x i y i+1 − x i+1 y i ) appears in both expressions for the area and the centroid of the polygon and it is therefore beneficial to compute the polygon area and the centroid in the same loop. These formulas can be motivated by decomposing a given simple polygon into triangles and summing up the individual centers weighted by the individual triangle areas [33]. The resulting formulas are found in [34]. The area of the rectangle (if fully within a viewing cone) or the intersection polygon between a viewing cone and the rectangle determines the amount of radiation originating from the rectangle which can reach the detector associated with the viewing cone. The centroid of the intersection polygon is used to determine the distanceL between the detector and the center of emission from the intersection polygon to account for the radial decay of isotropic emission from the respective emissive rectangle.
A.2.6. Angle of incidence for a viewing cone The angle of incidence of the central line of a viewing cone onto a detector is given by where φ viewingcone is the angle of the connection between the intersection polygon centroid (c x , c y ) and the viewing cone origin on the detector (x, y) in global coordinates (see figure A10). The inclination angle of the detector in global coordinates is given by φ det . The detector signal is scaled by a factor cos(φ i ) to account for the projection of the detector area into a plane perpendicular to central line-of-sight of the viewing cone. This is equivalent to the step from equations (23)-(24).

A.3. Verification against model problems
The verification is assessed by comparing the two-dimensional numerical results from the volumetric raytracing method to results from the two-dimensional analytical model problem (see section A.1.1). Also, the results from volumetric raytracing have been scaled to a three-dimensional model problem and compared to an existing forward modeling code htforward, which uses a three-dimensional raytracing approach. A convergence study for the numerical solution to the two-dimensional model problem consists e.g. of gradually refining the discretization of the active surface of the detector into n d elements and the discretization of each of the corresponding viewing cones into n a parts. The intersection region between the emissive layer and the viewing cones are finitewidth polygons extending from the front side to the back side of the emissive layer. In the analytical model, the intersection region is also subdivided along the central line-of-sight. The width of the emissive layer is successively decreased with respect to the distance d between the detector and the aperture; this corresponds to using a finer and finer grid in the tomography algorithm. It is expected that the deviations between the analytical model and the volumetric raytracing method successively decrease with a decreasing width of the emissive layer. The model geometry is depicted in figure A11 and the result of this convergence study is shown in figure A12. For a given width of the emissive layer, a certain number of subdivisions for the detector area and the viewing cone opening angle is necessary to reach the final accuracy. When the width of the emissive layer is decreased, the achievable accuracy improves as expected, but at the cost of a higher number of subdivisions required to reach convergence. The pixel size chosen in this article for the analysis of the measurements is 4 cm and therefore it was decided to use n d = n a = 10 subdivisions for both the detector areas and the viewing cone opening angles. A relative accuracy of 1.0 × 10 −6 is to be expected for this choice, which is well below the expected uncertainties of the resulting photodiode signals.