A convolutional neural network approach for reconstructing polarization information of photoelectric X-ray polarimeters

TL;DR

This paper introduces a CNN-based data processing algorithm for extracting polarization information from photoelectron track images in X-ray polarimeters, improving sensitivity and reducing systematic errors.

Contribution

It develops a novel CNN approach with a specialized loss function to accurately predict polarization parameters and enhance event selection in X-ray polarimetry.

Findings

CNN models achieve less than 2 pixel position accuracy across energies

The method reduces systematic modulation to below 1%

Polarization sensitivity improves by 10-20% over previous methods

Abstract

This paper presents a data processing algorithm with machine learning for polarization extraction and event selection applied to photoelectron track images taken with X-ray polarimeters. The method uses a convolutional neural network (CNN) classification to predict the azimuthal angle and 2-D position of the initial photoelectron emission from a 2-D track image projected along the X-ray incident direction. Two CNN models are demonstrated with data sets generated by a Monte Carlo simulation: one has a commonly used loss function calculated by the cross entropy and the other has an additional loss term to penalize nonuniformity for an unpolarized modulation curve based on the $H$-test, which is used for periodic signal search in X-ray/$\gamma$-ray astronomy. The modulation curve calculated by the former model with unpolarized data has several irregular features, which can be canceled out by unfolding the angular response or simulating the detector rotation. On the other hand, the latter model can predict a flat modulation curve with a residual systematic modulation down to $\lesssim1$%. Both models show almost the same modulation factors and position accuracy of less than 2 pixel (or 240 $\mu$m) for all four test energies of 2.7, 4.5, 6.4, and 8.0 keV. In addition, event selection is performed based on probabilities from the CNN to maximize the polarization sensitivity considering a trade-off between the modulation factor and signal acceptance. The developed method with machine learning improves the polarization sensitivity by 10-20%, compared to that determined with the image moment method developed previously.

Equations

$$L=\sum H+{\lambda }_2\sum w^2,$$

$${\bar{\phi }}_j(p_{j,i})=atan2\left(\sum _{i=0}^{N_{\mathrm{c}}-1}{p}_{j,i}\sin {\phi }_i,\sum _{i=0}^{N_{\mathrm{c}}-1}{p}_{j,i}\cos {\phi }_i\right),$$

$$Z_m^2({\bar{\phi }}_j)=\frac{2}{N}\sum _{k=1}^m\left[{\left(\sum _{j=0}^{N-1}\cos k{\bar{\phi }}_j\right)}^2+{\left(\sum _{j=0}^{N-1}\sin k{\bar{\phi }}_j\right)}^2\right],$$

$$\frac{\partial Z_m^2}{\partial {\bar{\phi }}_j}=\frac{4}{N}\sum _{k=1}^{m}k\left[-\sin k{\bar{\phi }}_j\left(\sum _{j=1}^N\cos k{\bar{\phi }}_j\right)+\cos k{\bar{\phi }}_j\left(\sum _{j=1}^N\sin k{\bar{\phi }}_j\right)\right].$$

$$H_M=\underset{1\le m\le 20}{\max }(Z_m^2-4m+4)=Z_M^2-4M+4\ge 0,$$

$$L=\sum H+{\lambda }_2\sum w^2+{\lambda }_H{\bar{H}}_M,$$

Full text

1 Introduction

Recent developments in micropattern gas detectors enable us to track an electron with an energy down to 1 keV. These detectors can be used as X-ray polarimeters by imaging a photoelectron track induced by an X-ray and extracting the polarization from its initial direction. Since the photoelectric effect is the dominant interaction for X-rays, the micropattern gas polarimeter combined with an X-ray focusing optics is expected to achieve a much higher sensitivity than the Bragg diffraction polarimeter [1] that is the only satellite-borne polarimeter flown in space so far in the soft X-ray energy band less than 10 keV.

The differential cross section of K-shell photoelectrons in the non-relativistic region is proportional to $\cos^{2}(\phi)$, where $\phi$ is the azimuthal angle with respect to the X-ray electric vector [2]. The modulation amplitude, $a$, which is defined as the ratio of sinusoidal amplitude to unmodulated offset derived from the photoelectron angular distribution (or so-called the modulation curve), is intrinsically 100% for perfectly linearly polarized X-rays. However, the observed modulation amplitude declines due to imperfect angular reconstruction of photoelectron track images which are curved by Coulomb scattering and are blurred by electron diffusion (see Fig. 1 in Ref. 3). As accurate an angular reconstruction method as possible is preferred to achieve the modulation amplitude close to 100%. The modulation amplitude for perfectly polarized X-rays is called the modulation factor designated as $\mu$. The observed modulation amplitude is then give by $a=a_{p}\mu$, where $a_{p}$ is the polarization degree of a source.

So far, two major reconstruction methods have been demonstrated: one is with image moments (e.g. [4, 3]) and the other is based on the shortest path problem in graph theory [5]. When applying the two methods to track images from the polarimeter we have developed [6, 7], the results were found to be similar [3]. Since both methods eventually use image moments to reconstruct the initial photoelectron emission direction, it may show a limitation of the image moments method, which is representative measure but with the loss of some information. Furthermore, image pre-processing to calculate image moments loses a part of pixel information mainly due to image thresholding. Therefore, a further improvement of the polarimeter sensitivity requires a novel reconstruction algorithm which can deal with all the pixels of a track image.

An alternative approach to determine the initial direction and position of a photoelectron track is to learn the image features automatically with a machine. Machine learning is a rapidly evolving field of computer science and can provide many different approaches to make data-driven prediction by building a model from sample inputs. Among them, convolutional neural networks (CNNs) have generally shown an excellent performance in image recognition. For example, CNNs achieve a low misclassification error rate for the Modified National Institute of Standards and Technology (MNIST) database [8] on handwritten digits (see its official site111http://yann.lecun.com/exdb/mnist/ for recent results).

This paper describes polarization extraction and track selection based on a CNN. Section 2 briefly explains the CNN model developed for this work. We describe results derived from a commonly used loss function combined with the CNN model in Section 3.1 and a modified one for polarimetry in Section 3.2. In addition, track selection with CNN-derived products is demonstrated in Section 3.3. Lastly, we conclude the study in Section 4. Throughout this paper, all errors are given at the $1\sigma$ confidence level unless otherwise stated.

2 Convolutional neural network

2.1 Overview

The main goal in this paper is to extract polarization information from 2-D images of photoelectron tracks with a CNN. An artificial neural network is made up of several layers, each of which contains a number of nodes. A value of each node is computed as a bias plus a linear combination of (all or a part of) nodes in the previous layer, followed by a nonlinear transformation with a so-called activation function. The weights of the linear combination and biases are free parameters and are optimized through training with an example data set of input-output pairs by fitting outputs predicted from the network to corresponding inputs; that is supervised learning.

Supervised learning is mainly categorized into two methods: regression computing a continuous value and classification predicting a discrete label to which an input belongs. In classification, the output layer contains votes (or probabilities which are normalized votes) of each label and selects the largest vote (or highest probability) to predict the label. In this work, multi-class classification is employed because the highest probability is efficient in selecting good tracks for polarization measurements, as described in detail later. In addition, classification is generally more robust to outliers than standard regression with the least squares approach.

A CNN is a type of feed-forward artificial neural network where data pass only in the forward direction from an input layer, through multiple hidden layers, and finally to an output layer, with no loops or cycles. The hidden layers of a CNN contain a combination of three types of layers: convolutional, pooling, and fully connected layers. The convolutional layer is a key component of a CNN to automatically extract image features by performing a convolution operation on an input with a weight matrix called a filter or kernel. The operation is repeated by shifting the filter over an input image with a step size called a stride. The pooling layer performs nonlinear downsampling to reduce the size of the convolved images while maintaining their features. The fully connected layer connects every node in one layer to every node in another layer, or in mathematical words, multiplies the input nodes by a weight matrix and then adds a bias vector. Following a single- or multi-set of the convolutional and pooling layers, the fully connected layers take image features from the last pooling layer and classify them into various classes (or regress from them to a numerical value).

2.2 Network model

Figure 1 illustrates a schematic view of the CNN classification model that we designed for X-ray polarimetry. The number of network parameters is 0.74 million. The input layer consists of 900 nodes corresponding to the $30\times 30$ pixels of a photoelectron track image. Since polarization measurements requires at least the initial direction of photoelectron emission, the output layer should have angular predictions. In addition, for imaging capability, the output layer is preferred to include predictions for the initial photoelectron position or the X-ray interaction point. Therefore, the output layer is designed to simultaneously predict the three values: the initial azimuthal angle and 2-D position of photoelectron emission. Each prediction is represented by a probability histogram with 36 equal-width bins; for example, the probability distribution of the predicted angle ranging from -90 to 90 deg222Since the differential cross section of photoelectron emission is represented with a bimodal sinusoid, a twofold angular distribution is sufficient enough to determine polarization. is divided by 36 bins each having a width of 5 deg. Thus, the output layer has $36~{}\mathrm{bins}\times 3~{}\mathrm{predictions}=108$ nodes in total. Although the number of probability bins is arbitrary, its optimization is out of scope for this work.

The hidden layers are constructed based on the Visual Geometry Group (VGG) model [9], and are comprised of the convolution layers with the filer size $3\times 3$ and stride 1 and the maximum pooling layers with the filter size $2\times 2$ and stride 2, and the fully-connected layers. The VGG developers demonstrated that multiple continuous convolution layers with a size of $3\times 3$ is more effective than a single convolution layer with a wider kernel. All the layers except for the pooling layers and last hidden layer utilize the rectified linear unit (ReLU [10]) function as the activation function. The last fully-connected layer uses the softmax function to compute each of the three probability distributions.

For each event, the actual initial direction and 2-D position of the photoelectrons in the training data are separately converted into a one-hot vector where an element to which the answer belongs is 1 and all the other 35 elements are 0 to correspond with the bin size of the outputted probability histogram. In order to evaluate a difference between each of the three probability distributions and the one-hot answer vector, the cross entropy, $H$, between them is computed. The sum of the three cross entropies is used as the loss function to evaluate the fitting (or supervise the machine learning) and is minimized with an optimization algorithm. In addition, L2 regularization, which is the sum of the square of the weights, is added to $H$ to suppress overfitting. The loss function, $L$, of the network model is given by:

[TABLE]

where $w$ is a weight and $\lambda_{2}$ is a hyperparameter, which is fixed at a constant value set before the training process, for L2 regularization to control trade-off between overfitting and training efficiency.

In order to train the CNN model (or optimize its weights and biases), we employ the Adam optimizer [11], a commonly used type of stochastic gradient descent. In addition, mini-batch training with a size of 360 images is performed to efficiently find the minimum of $L$.

In general, a deeper feed-forward network can lead to better performance (e.g. [12]). Although another set of convolution and pooling layers, or fully-connected layers were added to the network in Fig. 1, the result did not significantly change. In addition, wider networks were found to show no improvement.

2.3 Input data: photoelectron track images

Supervised machine learning requires a training data set containing input-output pairs. In this case, the training set is a set of photoelectron track images with initial direction and 2-D position. The photoelectric polarimeter that we have developed can make a photoelectron track with $30\times 30$ pixels of a size of $121\times 121~{}\mu$m [7]. Although, in actual experiment, we can easily regulate the X-ray polarization direction, we cannot control the photoelectron direction which is randomized according to the differential cross section, $\cos^{2}(\phi)$. Therefore, a Monte Carlo simulation was performed to generate the training data.

The simulator designed and developed for the time projection chamber (TPC) polarimeter [13] can reproduce the passage of X-rays and electrons through dimethyl ether gas at a pressure of 190 Torr, and generate photoelectron track images with the same dimensions as physical ones. It accounts for detector effects such as image blurring due to diffusion of drifting electrons to a readout electrode, electron amplification by a gas electron multiplier (GEM), and image elongation due to signal shaping with a time constant of 50 ns. The readout scheme, which is triggered by the GEM cathode signal, is also simulated to take a 2-D track image with the TPC technique, where signal charges are read out by 1-D strip electrodes with continuous analog-to-digital converter (ADC) sampling at a frequency of 20 MHz. The readout image size is $30\times 30$ pixels, each having an ADC value.

The simulated images can be offline processed in the same manner as actual images as described in § 3 of Ref. 14. The offline image processing mainly includes deconvolution with the signal shaping response and then image thresholding to perform the conventional angular reconstruction methods. Example raw and corresponding processed images are shown in Fig. 2. The modulation factors derived from the processed images with the image moment method can be reproduced with an absolute error of $\pm 3$% (or a relative error of $\sim 5$%) in the 2.7–8.0 keV energy range, where the polarimeter performance was evaluated in detail with a synchrotron X-ray facility [14].

As input data to the network, we used the raw images elongated on the vertical axis in Fig. 2 by the signal shaping process rather than the offline processed images which were utilized for the angular reconstruction methods so far reported. This is because the image processing loses some of image information mainly by image thresholding which sets a charge amount to zero if it is below a set threshold. Furthermore, deconvolution of the signal shaping response makes the signal-to-noise ratio worse. Since a CNN is known as a robust model that needs relatively little image pre-processing, only the linear normalization is applied to our data sets so that the minimum and maximum values among all the pixels in all the images are 0 and 1, respectively.

The actual initial direction and 2-D position of photoelectron emission to supervise the network are individually discretized and grouped into 36 classes corresponding to the outputs of the CNN model. For example, the class $i$ ($i=0,1,\cdots,35$) for the photoelectron direction constitutes angles ranging from $5i-90$ to $5i-85$ deg. The emission position is originally calculated as a relative coordinate value in the detector coordinate system and is converted into a float value in the image pixel coordinate with a $30\times 30$ pixel size. Then, the emission angle ranging from $30i/36$ to $30(i+1)/36$ pixel is categorized as the class $i$. When the emission position is located outside of the image frame, the position class is floored/ceiled to 0/35.

Classification with imbalanced training data may result in poor performance (e.g. [15]). In order to avoid the imbalanced learning problem, the training data containing balanced angular classes were generated with the simulator by irradiation of unpolarized rather than polarized X-rays. In addition, a wide variety of raw track images were created by setting the primary X-ray energy ranging from 0.5 to 12 keV and irradiation distance from (or electron drift length to) the readout electrode ranging from 0.45 to 1.55 cm, the increase of each parameter makes a photoelectron track longer and more blurred, respectively. The X-ray momentum direction was fixed parallel to the optical axis of the polarimeter, while the X-ray polarization direction was set to be randomly distributed perpendicular to its momentum direction. The X-ray energy is randomly sampled between 0.5 and 12 keV with weighting by inverse of detection efficiency to make the resulting spectrum flatten. In this way, we generated 3,600,000 and 360,000 images for unpolarized X-rays as training and validation data sets, respectively.

In addition to the training and validation data sets, unpolarized and perfectly polarized data sets each containing monochromatic (2.7, 4.5, 6.4, and 8.0 keV, separately) X-rays were generated with the same simulator to test the network model after training. The irradiation position for the test is limited to between 0.6 and 1.0 cm unlike the training and verification data sets. In addition, for the polarized data, the X-ray polarization direction is inclined at $\phi=-45$ deg in the image coordinate defined in Fig. 2. These are because the polarized test data set is a simulation to reproduce those actually observed with an X-ray synchrotron radiation facility [14]. Each test data set at monochromatic energy contains 360,000 track images. Only for comparison purposes with the machine learning (ML) approach, the test raw images are processed in the conventional way separately with the image moment (IM) and graph-based (GR) methods to calculate the modulation factor.

3 Experiments

3.1 Basic model

We programmed the network model shown in Fig. 1 with Google TensorFlow [16] version 1.8.0 and trained it with a NVIDIA GeForce GTX 1080 graphics processing unit (GPU). At the beginning of training, all biases were set to zero, while all weights were initialized to uniform random values dedicated to the ReLU activation function according to Ref. 17. The hyperparameter $\lambda_{2}$ in Eq. (1) should be low to increase training efficiency without overfitting. It was optimized with the grid search method ranging from 0.0001 to 0.0010 with step 0.0001 and was set to 0.0007 because overfitting appears if $\lambda_{2}\leq 0.0006$.

We used the 3,600,000 images for training and performed stochastic mini-batch training with a size of 360 images to efficiently find the minimum of the loss function of Eq. (1) and suppress overfitting. After every epoch, which means an iteration over all images in the training data set, we evaluated the cross entropy (or goodness of fitting) with the 360,000 validation images. A learning curve shown in Fig. 3 shows that training is successful without overfitting – the loss of both training and validation data sets decreases and approaches to a constant value as training progresses. The total computation time for 600 epochs with the GPU is approximately 45 hours.

During training, we evaluated the network model after every epoch by inputting the validation data set and stored the one with the minimum of the loss. Figure 4 (a) shows the modulation curve created by predicted angles from the best (or the lowest loss) model, which were saved at the epoch of 460, by inputting the simulated test data sets of unpolarized 8.0 keV X-rays. For comparison, the modulation curve calculated with the IM method is also shown in Fig 4 (a). Although a modulation curve for unpolarized X-rays is preferred to be directionally uniform to simplify the polarimetry analysis (e.g. [18]), the unpolarized modulation curve generated by the ML model shows an obvious residual modulation with several irregular features. On the other hand, the IM modulation curve is flat with a residual modulation amplitude, $a$, of $0.9\pm 0.2$%, which is derived by fitting it with a sinusoidal model of $C\{a\cos[2(\phi-\phi_{0})]+1\}$, where $C$ is a constant factor or unmodulated offset and $\phi_{0}$ is the polarization angle. In a similar way, the polarized modulation curve shown in Fig 4 (b) is not sinusoidal, while that obtained with IM shows a unimodal sinusoidal curve with $a=\mu=58.1\pm 0.2$% and $\phi_{0}=-44.7\pm 0.1$ deg.

Two methods are suggested to cancel out the systematic irregular modulation and extract polarization information from the modulation curves. The first method is unfolding the polarized modulation curve with the unpolarized one. The simplest method of an unfolding is dividing the polarized curve by the unpolarized one333Similar method to the so-called Crab ratio to unfold a spectral response function [19], which divides a source spectrum by the Crab spectrum observed with the same detector.. The divided modulation curve in Fig. 5 shows a unimodal sinusoidal shape with $\mu=69.2\pm 0.2$% and $\phi_{0}=-45.3\pm 0.1$ deg. A $\chi^{2}$ value of the fit is 279.8 with 33 degrees of freedom (dof), indicating that there still remain small residuals from the best-fit sinusoidal model because the division method is only an approximate unfolding. To be more accurate, an unfolding with a 2-d response function of the X-ray polarization angle vs. the predicted one is needed.

The second method needs uniform rotation of polarimeter or satellite with respect to the optical axis. In fact, the PRAXyS X-ray polarimetry mission had a plan to rotate its satellite with rotations per minute of 10 to reduce a residual modulation [20]. The rotation is simulated with the same unpolarized data that the first method is applied to, by shifting predicted angles so that the X-ray polarization angle, which is initially set to a random value, goes to -45 deg in the modulation curve coordinate. The resulting modulation curve is similar to Fig. 5 (b), but has smaller residuals from the best-fit sinusoidal model with $\chi^{2}/\mathrm{dof}=33.2/33$. The modulation factor and polarization degree are respectively $69.0\pm 0.2$% and $-45.1\pm 0.1$ deg, which are consistent with those determined with the first method.

Table 3.1 summarizes the fitting results for all the monoenergetic test data sets with the different angular reconstruction methods including the conventional IM and GR methods and more sophisticated ones described later. It shows that machine learning can improve the modulation factors by $\sim 10$% and $\sim 20$% in the lower and higher energy ranges, respectively.

Best-fit parameters of the sinusoidal curve model to modulation curves generated with the different angular reconstruction methods.

Energy Method

Unpolarized

Polarizeda

(keV)

$\mu$ (%) $\chi^{2}$ b

$\mu$ (%) $\phi_{0}$ (deg) $\chi^{2}$ b

2.7 IMc

$0.4\pm 0.2$ 23.1

$23.4\pm 0.2$ $-45.1\pm 0.3$ 41.4

GRd

$2.4\pm 0.2$ 49.2

$22.6\pm 0.2$ $-47.7\pm 0.3$ 57.5

BML_dive

$26.0\pm 0.3$ $-44.3\pm 0.4$ 363

BML_rotf

$25.7\pm 0.2$ $-45.0\pm 0.3$ 41.9

MMLg

$0.8\pm 0.2$ 147

$24.9\pm 0.2$ $-45.7\pm 0.3$ 162

MML_divh

$24.8\pm 0.3$ $-44.9\pm 0.4$ 59.2

MML_roti

$24.0\pm 0.2$ $-44.5\pm 0.3$ 30.1

4.5 IM

$0.6\pm 0.2$ 145j

$42.0\pm 0.2$ $-44.6\pm 0.2$ 173j

GR

$1.4\pm 0.2$ 221j

$41.0\pm 0.2$ $-43.8\pm 0.2$ 283j

BML_div

$46.8\pm 0.3$ $-45.0\pm 0.2$ 347

BML_rot

$47.3\pm 0.2$ $-45.2\pm 0.2$ 44.2

MML

$0.8\pm 0.2$ 106

$47.6\pm 0.2$ $-44.9\pm 0.1$ 145

MML_div

$47.0\pm 0.3$ $-44.6\pm 0.2$ 49.3

MML_rot

$47.5\pm 0.2$ $-44.8\pm 0.1$ 35.6

6.4 IM

$0.8\pm 0.2$ 45.8

$54.3\pm 0.2$ $-44.6\pm 0.1$ 34.1

GR

$1.3\pm 0.2$ 114

$52.6\pm 0.2$ $-44.2\pm 0.1$ 119

BML_div

$62.6\pm 0.3$ $-45.2\pm 0.1$ 403

BML_rot

$62.3\pm 0.2$ $-45.0\pm 0.1$ 22.5

MML

$1.1\pm 0.2$ 57.5

$62.5\pm 0.2$ $-45.3\pm 0.1$ 96.3

MML_div

$62.2\pm 0.3$ $-44.9\pm 0.1$ 29.6

MML_rot

$62.3\pm 0.2$ $-45.0\pm 0.1$ 18.8

8.0 IM

$0.9\pm 0.2$ 56.8

$58.1\pm 0.2$ $-44.7\pm 0.1$ 39.7

GR

$0.3\pm 0.2$ 81.1

$58.9\pm 0.2$ $-45.0\pm 0.1$ 101

BML_div

$69.2\pm 0.2$ $-45.3\pm 0.1$ 278

BML_rot

$69.0\pm 0.2$ $-45.1\pm 0.1$ 33.2

MML

$0.6\pm 0.2$ 53.8

$68.9\pm 0.2$ $-45.2\pm 0.1$ 101

MML_div

$68.9\pm 0.2$ $-45.0\pm 0.1$ 43.7

MML_rot

$69.1\pm 0.2$ $-45.0\pm 0.1$ 31.3

• a

Perfectly polarized data containing X-rays with a polarization angle of -45 deg.

• b

All degrees of freedom for the $\chi^{2}$ test are 33.

• c

Image moment (IM) method with processing parameters adjusted by actual data [3].

• d

Graph-based (GR) method with parameters adjusted by actual data [3].

• e

Basic machine learning (BML) method with the division approach.

• f

BML method with the polarimeter rotation technique.

• g

Modified machine learning (MML) method

• h

MML method with the division approach for the response unfolding.

• i

MML method with the polarimeter rotation technique.

• j

Relatively poor fit due to difficulty in accurately simulating photoelectron tracks with an energy of $\sim 4.5$ keV, where the track range is comparable to the electron diffusion size.

In addition to the initial direction of the photoelectron, the developed network is designed to simultaneously predict the initial position. A cumulative distribution of the distance between predicted and true positions shown in Fig. 6 displays that half power diameters (HPDs) at 2.7, 4.5, 6.4, and 8.0 keV are 1.9, 1.6, 1.5, and 1.6 pixel (or 240, 190, 180, and 200 $\mu$m), respectively. The HPDs decrease with increasing energy except for 8.0 keV. The exception is because some of longer tracks have an interaction position that protrudes outside the captured image frame. The machine-predicted HPDs are improved compared with graph-based ones which are 2.3, 3.4, 3.2, and 3.7 pixel at 2.7, 4.5, 6.4, and 8.0 keV, respectively.

In order to evaluate result reproducibility provided by machine learning, the network model is retrained 10 times in the same way but with different randomized initial values. Although the irregular residual modulation pattern varies with every training, the 8.0 keV modulation factors derived with the uniform rotation method are scattered with the standard deviation of 0.1% and the average value of 68.9%. The standard deviation for the other energies are almost the same level. In addition, position accuracy is found to be unchanged. Therefore, the developed network can reproduce the results within 0.1% for all the energies.

3.2 Modified model to uniform predictions

Although machine learning can improve the modulation factors by 10–20% as described in § 3.1, the unpolarized modulation curve is not flat. The non-flat angular response can potentially cause the dependence of the modulation factor on the polarization angle. According to Ref. 18, if the polarization-angular response depends on the difference of the polarization angle and the observed angle and thus becomes flat for unpolarized X-rays, the angular response function can be separated into two independent functions. One is the 2-D energy response matrix used widely for spectral fitting. The other is the modulation function which is a 1-d vector as a function of the intrinsic energy and can be processed in the similar way as an effective area (or a so-called ancillary response) function in the XSPEC spectral fitting package [21]. Therefore, in order to simplify the polarimetry analysis and use the XSPEC heritage, the modulation curve for unpolarized X-rays is preferable to be flat.

In order to obtain a uniform (flat) angular response for unpolarized X-rays, a nonuniform penalty term should be added to the loss function. An additional loss term must have a gradient (or be differentiable) with respect to individual predictions, which are inputs to the operation to calculate a nonuniformity, to perform backpropagation. However, the histogram process to create a modulation curve and evaluate an angular uniformity is not continuously differentiable. Therefore, an unbinned method to evaluate a uniformity is required.

So far, the most probable angle is employed as the angular class with the highest predicted probability. However, this operation is also not continuously differentiable. Instead, the circular mean of the angular probability distribution is used as the predicted angle by calculating as follows:

[TABLE]

where $\bar{\varphi}_{j}$ is the circular mean angle of the $j$-th event in a mini-batch, $N_{\mathrm{c}}$ is the number of angular classes, $p_{j,i}$ is a predicted probability in the $i$-th class for the $j$-th event, and $\varphi_{i}$ is the center angle of the $i$-th class ($\varphi_{i}=\pi(2i+1)/N_{\mathrm{c}}$). It should be noted that $\varphi$ ranges from -180 to 180 deg, while the predicted angle $\phi$ ranges from -90 to 90 deg. This is because the circular mean is valid for angles ranging from -180 to 180 deg. It is also worth noting that Eq. (2) is differentiable with respect to the input $p_{j,i}$ because it consists of a combination of trigonometric functions. The angle $\bar{\varphi}_{j}$ is a continuous value unlike a discrete angle predicted in the previous section.

As an unbinned method to evaluate a directional uniformity, the $H$-test has been widely used for periodic signal search in X-ray/$\gamma$-ray astronomy [22]. The $H$-test is based on the test statistics $Z^{2}_{m}$ that is the sum of the Fourier powers of the first $m$ harmonics as given by:

[TABLE]

where $N$ is the number of events considered. The null distribution of $Z^{2}_{m}$ approaches a $\chi^{2}$ distribution with $m\times 2$ degrees of freedom as $N$ increases. The gradient of Eq. (3) with respect to $\bar{\varphi}_{j}$ can be easily described as:

[TABLE]

By using $Z^{2}_{m}$, the $H$ statistics, $H_{M}$ described to avoid confusion with the cross entropy $H$, is defined as:

[TABLE]

where $M$ is $m$ with the maximum of $Z^{2}_{m}-4m+4$. The $H$ statistics is non-negative because $H_{M}\geq H_{1}=Z^{2}_{1}\geq 0$. Furthermore, the deviation of $H_{M}$ is $\partial H_{M}/\partial\bar{\varphi}_{j}=\partial Z^{2}_{M}/\partial\bar{\varphi}_{j}$. Therefore, it can be used as the loss function penalizing angular nonuniformity – learning with backpropagation proceeds by minimizing it.

In the practical observation, spectral polarimetry with energy-resolved modulation curves will be conducted. Therefore, energy-resolved $H_{M}$ values are calculated so that each modulation curve becomes flat. For that reason, the mini-batch data containing 0.5–12 keV X-rays are separated by the incident X-ray energy into 40 ranges, each having $\sim 0.3$ keV width. Then, the mean value, $\bar{H}_{M}$, of non-zero elements is calculated along the 40 energy ranges.

A modified loss function by adding the $H$ statistics to the loss fuction in Eq. (1) is given by

[TABLE]

where the first and third terms respectively represent the classification performance and the angular uniformity. $\lambda_{H}$ is a hyperparameter for the $H$ statistics to control the trade-off between the two qualities and set to be 0.01. The hyperparameter $\lambda_{2}$ is set to be the same value as in the previous learning, 0.007, because the number of weights making up the network model is unchanged.

Machine learning is performed with the same manner and same training and verification data sets as described in § 3.1, but with consideration of the uniform angular predictions by replacing the loss function with Eq. (6) and with a larger mini-batch size. The mini-batch size should be large enough to reduce the statistical error of uniformity and is therefore set to be 18,000 limited by the GPU memory size of 8 GB. Additionally, the training images with the total size of 3,600,000 are sorted by irradiation position (or electron drift length) and then are grouped into 200 mini-batches. This is because the predicted modulation curve should be flat along the irradiation position as well as the incident X-ray energy. Although energy- and position-sorted images are inputted, learning efficiency becomes low due to biased track shapes in a mini-batch. Instead, uniformity along energies is evaluated with energy-resolved $H_{M}$ in Eq. (6).

In the same way as described in § 3.1, the network model is learned, evaluated after finishing every epoch, saved at the lowest loss, and tested by inputting the unpolarized and polarized monoenergetic test data sets. The machine learning lasts for 1500 epochs or $\sim 110$ hours because learning efficiency is lower than that in the previous model. The 8.0 keV modulation curves in Fig. 7, which are derived from predicted angles $\bar{\phi}=\bar{\varphi}/2$ with the circular mean in Eq. (2), clearly illustrate that the residual modulation in the unpolarized data is much suppressed. In addition, the polarized modulation is close to a unimodal sinusoid. However, the residual modulation still remains indicated by a reduced $\chi^{2}$ value of $\sim 3$, which is derived by fitting the sinusoidal curve model to the modulation curve. The fitting results including those for the other energies are listed in Table 3.1.

The residual systematic modulation can be canceled out by the division or rotation method as already described in § 3.1. The best-fit parameters with the two methods are also listed in Table 3.1. It shows that both models can improve $\chi^{2}$. In particular, the divided modulation curves can be reproduced by the sinusoidal curve model with reasonable $\chi^{2}$ because it is empirically known that the division method becomes better approximation for response unfolding as the denominator histogram is more featureless. The modulation factors are degraded only by $<2$% as an absolute difference, compared to those predicted by the basic network model. The best-fit polarization degrees at the four energies are scattered around the set polarization angle of -45 deg within $<1$ deg. In addition, it is also confirmed that the determination accuracy of the initial photoelectron position is unchanged.

3.3 Track selection to improve sensitivity

Reduction of tracks, from which it is difficult to extract the initial photoelectron angle, can improve the polarization sensitivity. So far, track eccentricity and circularity were demonstrated to quantify track roundness and discard round tracks (e.g. Refs. 23, 24, 14, 3). However, it is complicated to optimize a cut region in the two (or more) parameter plane (or space). Instead, the maximum value of the angular probability distribution, $p_{\mathrm{max}}$, for each track outputted from the network model is suggested.

The test data sets, which were processed with the modified network model at the minimum loss and with uniform rotation along the optical axis in § 3.1, are used to demonstrate track selection with $p_{\mathrm{max}}$. As shown in Fig. 8, the modulation factors derived from modulation curves sorted by $p_{\mathrm{max}}$ monotonically increase with $p_{\mathrm{max}}$. Therefore, we set a threshold for $p_{\mathrm{max}}$ to separate valid from invalid tracks and discard tracks with lower $p_{\mathrm{max}}$.

In order to evaluate and optimize track selection with $p_{\mathrm{max}}$, the figure of merit, $F$, is defined in the same way as Ref. 3 as given by $F\equiv\mu\sqrt{\varepsilon}$, where $\varepsilon$ is signal acceptance decreased by track selection and $\mu$ is the modulation factor of selected events. As shown in Fig. 9, $\varepsilon$ monotonically decreases with increasing the threshold for $p_{\mathrm{max}}$, while $\mu$ monotonically increases. We search for the best selection condition where $F$ for 8.0 keV X-rays is maximized by discarding tracks in order of increasing $p_{\mathrm{max}}$, and find that $F$ with $p_{\mathrm{max}}>0.0593$, where $\varepsilon=91.3$% and $\mu=73.8$%, is improved by 2.0% compared to $F$ with no track selection. On the other hand, conventional track selection with eccentricity derived from the corresponding processed image can also improve $F$ by 0.03%. The results for the other test energies are obtained in the same manner and listed in Table 3.3. It shows that $p_{\mathrm{max}}$ is a better selector than eccentricity for all the energies.

Table 3.3 also lists reduction rates of the observation time to achieve a given polarimeter sensitivity, which increases with the square root of the observation time if the background rate is negligible. The machine-learning-based method for angular reconstruction and track selection can reduce the observation time by $\sim$20–30%, compared to the conventional IM method.

Comparison of polarimeter sensitivity derived from the various angular reconstruction methods combined with track selection.

Energy Method $\mu$ $\phi_{0}$ $\chi^{2}$ a $t$b $\varepsilon$ $F$c $R_{T}$d

(keV)

(%) (deg)

2.7 IMe $23.4\pm 0.2$ $-45.1\pm 0.3$ 41.4 0 1 0.234 1

IM_ecutf $27.7\pm 0.3$ $-45.2\pm 0.3$ 33.0 0.490 0.766 0.243 0.927

MML_rotg $24.0\pm 0.2$ $-44.5\pm 0.3$ 30.1 0 1 0.240 0.951

MML_rot_ecuth $26.2\pm 0.2$ $-44.6\pm 0.3$ 29.2 0.431 0.860 0.242 0.935

MML_rot_pcuti $30.5\pm 0.3$ $-44.8\pm 0.3$ 31.0 0.0388 0.708 0.257 0.829

4.5 IM $42.0\pm 0.2$ $-44.6\pm 0.2$ 173 0 1 0.420 1

IM_ecut $45.1\pm 0.2$ $-44.6\pm 0.2$ 181 0.550 0.896 0.427 0.967

MML_rot $47.5\pm 0.2$ $-44.8\pm 0.1$ 35.6 0 1 0.475 0.782

MML_rot_ecut $50.2\pm 0.2$ $-44.7\pm 0.1$ 36.6 0.538 0.908 0.478 0.772

MML_rot_pcut $53.4\pm 0.2$ $-44.8\pm 0.1$ 38.7 0.0488 0.834 0.488 0.741

6.4 IM $54.3\pm 0.2$ $-44.6\pm 0.1$ 34.1 0 1 0.543 1

IM_ecut $57.7\pm 0.2$ $-44.5\pm 0.1$ 30.7 0.705 0.908 0.550 0.975

MML_rot $62.3\pm 0.2$ $-45.0\pm 0.1$ 18.8 0 1 0.623 0.760

MML_rot_ecut $62.6\pm 0.2$ $-45.0\pm 0.1$ 18.8 0.398 0.993 0.623 0.760

MML_rot_pcut $68.6\pm 0.2$ $-45.0\pm 0.1$ 21.3 0.0563 0.867 0.639 0.722

8.0 IM $58.1\pm 0.2$ $-44.7\pm 0.1$ 39.7 0 1 0.581 1

IM_ecut $60.2\pm 0.2$ $-44.7\pm 0.1$ 43.7 0.710 0.946 0.585 0.986

MML_rot $69.1\pm 0.2$ $-45.0\pm 0.1$ 31.3 0 1 0.691 0.707

MML_rot_ecut $69.2\pm 0.2$ $-45.0\pm 0.1$ 31.6 0.325 0.998 0.692 0.705

MML_rot_pcut $73.8\pm 0.2$ $-45.0\pm 0.1$ 36.4 0.0593 0.913 0.705 0.679

• a

All degrees of freedom for the $\chi^{2}$ test are 33.

• b

Threshold $t$ of eccentricity for the IM method and the maximum probability $p_{\mathrm{max}}$ for the MML method. Tracks below the threshold are discarded.

• c

Figure of merit, $F\equiv\mu\sqrt{\varepsilon}$ (higher is better), representing the polarimeter sensitivity.

• d

Reduction rate of the observation time to achieve a given polarimetry sensitivity, compared to the IM method.

• e

Same as the IM method in Table 3.1

• f

IM method with track selection by using eccentricity.

• g

Same as the MML_rot method in Table 3.1.

• h

MML_rot method with track selection by using eccentricity.

• i

MML_rot method with track selection by using $p_{\mathrm{max}}$.

4 Summary and Discussion

We have developed angular and positional reconstruction of a photoelectron track image for the TPC X-ray polarimeter by using supervised machine learning. Data sets for training, validation, and test are generated by the Monte Carlo simulation we have developed. We test two CNN classifications; the VGG-based model and the modified one with a penalty term for nonuniformity of angular predictions. Although the former network outputs large residual modulation for unpolarized X-rays (see Fig. 4), it can be canceled out by unfolding the response or rotating the polarimeter. The resulting modulation factor is improved by $\sim$10–20% for all the energies, compared to the conventional reconstruction methods. The latter network considers uniformity of angular predictions, and therefore residual modulation for unpolarized X-rays is drastically reduced with decreasing the modulation factor within 2% (see Fig. 7 and Table 3.1). In addition, we suggest track selection with the maximum probability of predictions to further improve the polarimeter sensitivity.

The modified network with the nonuniform penalty still shows the residual modulation, as indicated by the relatively large $\chi^{2}$ (for the MML method in Table 3.1). The residual modulation may be reduced by increasing the mini-batch size, though it is limited by the GPU memory size. Therefore, a GPU with a larger memory size is preferable. Alternatively, a more recent network with smaller memory consumption but deeper layers based on e.g. Inception [25] and/or ResNet [26] architecture could increase the mini-batch size as well as improve the modulation factor (see Ref. 27 for information on network architecture comparisons).

The developed method would be applicable to not only the other photoelectric polarimeters, such as the Gas Pixel Detector (GPD) polarimeter [28, 29], but also electron-tracking Compton polarimeters which are capable of determining the initial emission direction of recoil electrons to improve Compton kinematics reconstruction (e.g. [30, 31]). However, the GPD polarimeter takes track images in a hexagonal grid to reduce residual modulation, and therefore requires a hexagonal convolution operation. A simple way demonstrated in Ref. 32, where the CNN is applied to H.E.S.S. hexagonal images, is that the hexagonal convolution can be performed with a combination of available operations for the square convolution.

5 Future work

Observed track images can be input to the network model trained by the simulated images to extract polarization information. However, the residual modulation is most likely to appear due to systematic effects induced by a real polarimeter and an imperfect simulation. Although the residual modulation can be canceled out with the division method or the rotation method, the modulation factor is probably degraded compared to that derived from the simulated data. Training with observed rather than simulated data should be better for practical polarimetry observation.

The developed method needs the initial emission angle of each photoelectron track, and therefore demands as accurate a simulator as possible to generate a realistic image with the known angle. However, the network model can be potentially supervised with the modulation factor instead of the angles. Our method calculates directional uniformity, which is a statistic of samples, and uses it for training of machine learning. The modulation factor is also a sample statistic and can be derived from the Stokes parameters on an event-by-event basis [33] without using a histogram to create a modulation curve. The modulation factor, $\mu$, can be alternatively represented by $\mu=2\sqrt{\left(\sum^{N-1}_{i=0}\cos\bar{\varphi_{i}}\right)^{2}+\left(~{}\sum^{N-1}_{i=0}\sin\bar{\varphi_{i}}\right)^{2}}/N=\sqrt{2Z^{2}_{1}/N}$, which ranges from 0 to 1. The polarization angle, $\varphi_{0}$, can also be given by $\varphi_{0}=\operatorname{atan2}\left(\sum^{N-1}_{i=0}\sin\bar{\varphi_{i}},~{}\sum^{N-1}_{i=0}\cos\bar{\varphi_{i}}\right)/2$. These operations are differentiable and therefore are available for backpropagation to supervise machine learning. The cross entropy $H$ in Eq. (6) can be replaced with the two statistical values to represent the learning score, meaning that the revised method can input an observed data set with known Stokes parameters (or the polarization degree and angle) of incident X-rays – then a simulator is no longer needed because the initial angle of each image is not needed for training. The above revised method will be demonstrated and discussed elsewhere.

Acknowledgments

The authors would like to acknowledge helpful discussion with Keith Jahoda on polarimeter development. This work was partially supported by JSPS KAKENHI Grant Numbers JP17K18776 and JP16H02198.