Interpolation¶

The Interpolation class implements the \(\eta\)-interpolation method. This interpolation technique is based on charge sharing: for detected photon hits (e.g. clusters), it refines the estimated photon hit using information from neighboring pixels.

The method relies on the so-called \(\eta\)-functions, which describe the relationship between the energy measured in the central cluster pixel (the initially estimated photon hit) and the energies measured in its neighboring pixels. Depending on how much energy each neighboring pixel receives relative to the central pixel, the estimated photon hit is shifted toward that neighbor by a certain offset to the actual photon hit position in the pixel \((x, y)\).

The mapping between the \(\eta\) values and the corresponding spatial photon position \((x,y)\) can be viewed as an optimal transport problem.

One can readily compute the probability distribution \(P_{\eta}\) of the \(\eta\) values by forming a 2D histogram. However, the probability distribution \(P_{x,y}\) of the true photon positions is generally unknown unless the detector is illuminated uniformly (i.e. under flat-field conditions). In a flat-field, the photon positions are uniformly distributed.

With this assumption, the problem reduces to determining a transport map \(T:(\eta_x,\eta_y) \rightarrow (x,y)\), that pushes forward the distribution of \((\eta_x, \eta_y)\) to the known uniform distribution of photon positions of a flatfield.

The map \(T\) is given by:

\[\begin{split}\begin{align*} T_1: & F_{x}^{-1} F_{\eta_x|\eta_y} \\ T_2: & F_{y}^{-1} F_{\eta_y|\eta_x}, \end{align*}\end{split}\]

where \(F_{\eta_x|\eta_y}\) and \(F_{\eta_y|\eta_x}\) are the conditional cumulative distribution functions e.g. \(F_{\eta_x|\eta_y}(\eta_x', \eta_y') = P_{\eta_x, \eta_y}(\eta_x \leq \eta_x' | \eta_y = \eta_y')\). And \(F_{x}\) and \(F_{y}\) are the cumulative distribution functions of \(x\) and \(y\). Note as \(x\) and \(y\) are uniformly distributed \(F_{x}\) and \(F_{y}\) are the identity functions. The map \(T\) thus simplifies to

\[\begin{split}\begin{align*} T_1: & F_{\eta_x|\eta_y} \\ T_2: & F_{\eta_y|\eta_x}. \end{align*}\end{split}\]

Note that for the implementation \(P_{\eta}\) is not only a distribution of \(\eta_x\), \(\eta_y\) but also of the estimated photon energy \(e\). The energy level correlates slightly with the z-depth. Higher z-depth leads to more charge sharing and a different \(\eta\) distribution. Thus we create a mapping \(T\) for each energy level.

\(\eta\)-Functions:¶

template<typename T> struct Eta2¶

eta struct

Public Members

double x = {}¶: eta in x direction

double y = {}¶: eta in y direction

corner c = {0}¶: index of subcluster with highest energy value (given as corner relative to cluster center)

T sum = {}¶: photon energy (cluster sum)

Note

The corner value c is only relevant when one uses calculate_eta_2 or calculate_full_eta2. Otherwise its default value is cTopLeft.

Supported are the following \(\eta\)-functions:

\(\eta\)-Function on 2x2 Clusters:¶

\[\begin{equation*} {\color{blue}{\eta_x}} = \frac{Q_{1,1}}{Q_{1,0} + Q_{1,1}} \quad \quad {\color{green}{\eta_y}} = \frac{Q_{1,1}}{Q_{0,1} + Q_{1,1}} \end{equation*}\]

The \(\eta\) values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change. If the center pixel is in the bottom left pixel \(\eta_x\) will be close to zero. If the center pixel is in the bottom right pixel \(\eta_y\) will be close to 1.

One can apply this \(\eta\) not only on 2x2 clusters but on clusters with any size. Then the 2x2 subcluster with maximum energy is choosen and the \(\eta\) function applied on the subcluster.

template<typename ClusterType, typename = std::enable_if_t<is_cluster_v<ClusterType>>> std::vector<Eta2<typename ClusterType::value_type>> aare::calculate_eta2(const ClusterVector<ClusterType> &clusters)¶: Calculate the eta2 values for all clusters in a ClusterVector.

template<typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY, typename CoordType = uint16_t> Eta2<T> aare::calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl)¶: Calculate the eta2 values for a generic sized cluster and return them in a Eta2 struct containing etay, etax and the index (as corner) of the respective 2x2 subcluster relative to the cluster center.

Full \(\eta\)-Function on 2x2 Clusters:¶

\[\begin{equation*} {\color{blue}{\eta_x}} = \frac{Q_{0,1} + Q_{1,1}}{\sum_i^{1}\sum_j^{1}Q_{i,j}} \quad \quad {\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_i^{1}\sum_j^{1}Q_{i,j}} \end{equation*}\]

template<typename ClusterType, typename = std::enable_if_t<is_cluster_v<ClusterType>>> std::vector<Eta2<typename ClusterType::value_type>> aare::calculate_full_eta2(const ClusterVector<ClusterType> &clusters)¶: Calculate the full eta2 values for all clusters in a ClusterVector.

template<typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY, typename CoordType> Eta2<T> aare::calculate_full_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl)¶: Calculate the eta2 values for a generic sized cluster and return them in a Eta2 struct containing etay, etax and the index (as corner) of the respective 2x2 subcluster relative to the cluster center.

Full \(\eta\)-Function on 3x3 Clusters:¶

\[\begin{equation*} {\color{blue}{\eta_x}} = \frac{\sum_{i=0}^{2} Q_{i,2} - \sum_{i=0}^{2} Q_{i,0}}{\sum_{i=0}^{2}\sum_{j=0}^{2} Q_{i,j}} \quad \quad {\color{green}{\eta_y}} = \frac{\sum_{j=0}^{2} Q_{2,j} - \sum_{j=0}^{2} Q_{0,j}}{\sum_{i=0}^{2}\sum_{j=0}^{2} Q_{i,j}} \end{equation*}\]

The \(\eta\) values can range between -0.5,0.5.

template<typename ClusterType, typename = std::enable_if_t<is_cluster_v<ClusterType>>> std::vector<Eta2<typename ClusterType::value_type>> aare::calculate_eta3(const ClusterVector<ClusterType> &clusters)¶: Calculate eta3 for all 3x3 clusters in a ClusterVector.

template<typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY, typename CoordType = uint16_t> Eta2<T> aare::calculate_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl)¶

Cross \(\eta\)-Function on 3x3 Clusters:¶

\[\begin{equation*} {\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} + Q_{1,1} + Q_{1,2}} \quad \quad {\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} + Q_{1,1} + Q_{2,1}} \end{equation*}\]

The \(\eta\) values can range between -0.5,0.5.

template<typename ClusterType, typename = std::enable_if_t<is_cluster_v<ClusterType>>> std::vector<Eta2<typename ClusterType::value_type>> aare::calculate_cross_eta3(const ClusterVector<ClusterType> &clusters)¶: Calculate cross eta3 for all 3x3 clusters in a ClusterVector.

template<typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY, typename CoordType = uint16_t> Eta2<T> aare::calculate_cross_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl)¶

Interpolation class:¶

Warning

The interpolation might lead to erroneous photon positions for clusters at the borders of a frame. Make sure to filter out such cases.

Warning

Make sure to use the same \(\eta\)-function during interpolation as given by the joint \(\eta\)-distribution passed to the constructor.

Note

Make sure to use resonable energy bins, when constructing the joint distribution. If data is too sparse for a given energy the interpolation will lead to erreneous results.

class Interpolator¶

Public Functions

Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins, NDView<double, 1> ybins, NDView<double, 1> ebins)¶

Constructor for the Interpolator class.

Parameters:

etacube – joint distribution of etaX, etaY and photon energy (note first dimension is etaX, second etaY, third photon energy)
xbins – bin edges for etaX
ybins – bin edges for etaY
ebins – bin edges for photon energy

Interpolator(NDView<double, 1> xbins, NDView<double, 1> ybins, NDView<double, 1> ebins)¶

Constructor for the Interpolator class.

Parameters:

xbins – bin edges for etaX
ybins – bin edges for etaY
ebins – bin edges for photon energy

void rosenblatttransform(NDView<double, 3> etacube)¶

transforms the joint eta distribution of etaX and etaY to the two independant uniform distributions based on the Roseblatt transform for each energy level

Parameters:: etacube – joint distribution of etaX, etaY and photon energy (first dimension is etaX, second etaY, third photon energy)

inline NDArray<double, 3> get_ietax()¶

inline NDArray<double, 3> get_ietay()¶

template<auto EtaFunction = calculate_eta2, typename ClusterType, typename Eanble = std::enable_if_t<is_cluster_v<ClusterType>>> std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters) const¶

interpolates the cluster centers for all clusters to a better precision

Template Parameters:

ClusterType – Type of Clusters to interpolate
Etafunction – Function object that calculates desired eta default: calculate_eta2

Returns:

interpolated photons (photon positions are given as double but following row column format e.g. x=0, y=0 means top row and first column of frame) (An interpolated photon position of (1.5, 2.5) corresponds to an estimated photon hit at the pixel center of pixel (1,2))

template<typename T> Coordinate2D transform_eta_values(const Eta2<T> &eta) const¶

transforms the eta values to uniform coordinates based on the CDF ieta_x and ieta_y

Template Parameters:: eta – Eta to transform
Returns:: uniform coordinates {x,y}

Private Functions

template<typename T> std::pair<double, double> bilinear_interpolation(const size_t ix, const size_t iy, const size_t ie, const Eta2<T> &eta) const¶

bilinear interpolation of the transformed eta values

Parameters:

ix – index of etaX bin
iy – index of etaY bin
ie – index of energy bin

Returns:

pair of interpolated transformed eta values (ietax, ietay)

Private Members

NDArray<double, 3> m_ietax¶: marginal CDF of eta_x (if rosenblatt applied), conditional CDF of eta_x conditioned on eta_y value at (i, j, e): F(eta_x[i] | eta_y[j], energy[e])

NDArray<double, 3> m_ietay¶: conditional CDF of eta_y conditioned on eta_x value at (i,j,e): F(eta_y[j] | eta_x[i], energy[e])

NDArray<double, 1> m_etabinsx¶

NDArray<double, 1> m_etabinsy¶

NDArray<double, 1> m_energy_bins¶