To ease the reading, we provide a brief description of the retrigger mechanism principles, which have already been carefully addressed in [3, 6], with the help of Figure 1. In single X-ray photon counting circuitry, an event is detected and counted as soon as the rising edge of the incoming analog signal (shaperout) exceeds a predefined threshold (vth). Without a retrigger mechanism, during the time-over-threshold, the system is “dead” or insensitive to further incoming signals (pileup), leading to paralyzation at high fluxes (count paralyzed). Conversely, with the retrigger mechanism, as soon as the first analog signal exceeds the threshold, an additional internal digital dead time signal (deadgenout) is generated and periodically evaluated—with a configurable period that we call retrigger time τ _R (t _dead ). At every evaluation, if the incoming signal is still high, the counter is increased by a unit; if it has fallen below the threshold, the dead time signal is deactivated. Possible threshold crossings during the τ _R period do not give rise to counts. In this way, a non-paralyzable counting mode is achieved (count retriggered) with a limiting recorded rate of 1/τ _R in the high flux regime. The value of τ _R is customarily set slightly longer than the time-over-threshold of the input signal to prevent the system counting individual pulses multiple times.

The model is derived based on the following assumptions:

i. The events’ arrival follows Poisson statistics with an average true incoming rate n.

As for the monochromatic case described in [6], to the purpose of an analytical treatment of the problem, it is easier to focus on the rate loss n − m. This can be expressed for convenience as the sum of four separate contributions that we call l ₁, l ₂, l ₃, and l ₄ such that n − m = l 1 + l 2 + l 3 + l 4 . (1) l ₁ and l ₂ deal with the rate loss occurring in the “dead time” period after the detection of a signal or after an evaluation of the retrigger signal; l ₃ deals with the possibility of an incoming event being undetected in the first place; and l ₄ deals with the rate loss occurring during the time periods in between detected events.

In the following list, we describe the derivation of l ₁ to l ₄:

l 1 After a positive detection of an incoming event during the “live period” of the detector or after an evaluation of the retrigger signal, for a time interval of length Δτ = τ _R − τ _P , all incoming events are lost since they would not last until the next evaluation of the retrigger signal. Therefore, on average, nΔτ incoming events are lost for every recorded count m, yielding

Before making l ₂ explicit, it is necessary to introduce the following notation. Given s _E (ϵ) as the probability density function of the amplitude (which is proportional to its energy) of an individual signal, by definition the corresponding cumulative distribution function, S _E (ϵ) is S E ϵ = Pr ( E ≤ ϵ ) = ∫ 0 ϵ s E x d x . (3)

When i signals pileup, the stochastic variable E i = ∑ k = 1 i E k is distributed according to the convolution of the individual probability density function repeated i times, which is indicated as s E i ϵ = s E ϵ ∗ s E ϵ ∗ … ∗ s E ϵ ⏟ i = def s E ϵ * i , (4)

from which the corresponding cumulative distribution function can be deduced as S E i ϵ = Pr ( E i ≤ ϵ ) = ∫ 0 ϵ s E x * i d x . (5)

Let us also recall the probability of having i incoming events following the Poisson distribution in the time interval τ _P . Pr i | τ P = n τ P i e − n τ P i ! . (6)

We can now write l ₂ as l 2 = m ∑ i = 1 + ∞ i Pr i | τ P S E i E t h + ∑ i = 1 + ∞ i − 1 Pr i | τ P 1 − S E i E t h , (7)

where E_th is the threshold energy of the counting circuitry. Rearranging the terms, we can write that l 2 = m ∑ i = 1 + ∞ i Pr i | τ P S E i E t h + ∑ i = 1 + ∞ i Pr i | τ P 1 − S E i E t h − ∑ i = 1 + ∞ Pr i | τ P 1 − S E i E t h = m ∑ i = 1 + ∞ i Pr i | τ P − ∑ i = 1 + ∞ Pr i | τ P 1 − S E i E t h . (8)

Considering that ∑ i = 1 + ∞ i Pr ( i | τ P ) = ∑ i = 0 + ∞ i Pr ( i | τ P ) = n τ P and defining A = def ∑ i = 1 + ∞ Pr ( i | τ P ) S E i ( E t h ) , we can rewrite l ₂ as follows: l 2 = m n τ P + 1 − Pr 0 | τ P − A , (9)

and since Pr ( 0 | τ P ) = e − n τ P , we obtain l 2 = m n τ P − m + m e − n τ P + m A . (10)

l 3 An additional term has to be added to account for the possibility that an incoming event actually passes undetected; in the case when the time interval τ _P precedes each recorded count m, the overall signal lies below the threshold energy.

with the term 1/T arising for normalization purposes. Defining B = def ∑ i = 1 + ∞ i Pr ( i | τ P ) S E i ( E t h ) , we can rewrite l ₄ as follows: l 4 = 1 − m τ R τ P B . (13)

We can finally substitute the expressions of l ₁ to l ₄ of Eqs. 2, 10, 11, and 13, respectively, into Eq. 1, and by solving it for m, we obtain the final relation. m = n − B τ P n τ R + e − n τ P + 2 A − τ R τ P B . (14)

A few remarks were observed. Since the energy spectrum s _E (ϵ) is a free parameter of the model that usually lacks an analytical formulation, the terms A and B need to be evaluated numerically (typically with negligible computational burden). In addition, we left implicit in the notation the main functional dependencies of m = m(n, Eth). In the following, we will refer to the quantity m(n; Eth) (the semicolon specifies the fixed parameter) as the count rate curve, m(Eth; n) as the recorded integral energy spectrum, and ∂ m ( E t h ; n ) ∂ E t h as the recorded differential energy spectrum.

To validate the model, we used the numerical results of the custom-developed Monte Carlo model already used in [6, 7], consisting of the following. First, a temporal axis is created, spanning a statistically relevant length of several tens of ms and with a sufficiently small time step, typically on the order of ns. a number pulse-like events, computed according to the chosen value of the incoming rate, are distributed across the temporal axis, following a Poisson distribution. This resulting signal, which corresponds to the input of the detector system, is then convolved with the pulse shape of the analog front-end electronics, which is, in principle, arbitrary but is assumed to be a rectangle function in this context. To mimic the behavior of the counting electronics, a count is recorded for every rising edge crossing a specified threshold. If the retrigger capability is enabled, a count is also added for every time interval of length τ _R for which the signal lies above the threshold. The total recorded counts are then normalized to the overall temporal axis length to obtain the recorded rate. The entire operation is repeated for several values of the incoming rate and as a function of the threshold energy. The simulation error is estimated by subdividing the total time span into a set of smaller intervals of equal size, computing the variance over their ensemble, and scaling it to the overall interval time (see Appendix A).

The amplitudes of the incoming signals are randomly and independently extracted from a probability distribution function that corresponds to the specified energy spectrum. The energy spectrum was obtained as the convolution of the direct spectrum of an X-ray tube, computed with the model described in [18], with the response function of a pixelated counting detector, estimated based on the numerical results from FLUKA ³ —a Monte Carlo particle transport and interaction suite [19, 20]—and a custom-developed pixelization and counting routine. Specifically, the FLUKA code computes, for a large number of incoming X-ray photons, the individual energy depositions in the sensor volume at the microscopic level, including all secondary effects like fluorescence and escape. The pixelization routine translates deposited energy into electric charge, which is drift-diffused toward the collecting electrodes. Signals collected at the pixel side (comprehensive of charge sharing effects) are added to a Gaussian random signal representing the electronic noise and subsequently histogrammed to obtain the sought energy spectrum.

The predictions of the analytical model were compared with the results of Monte Carlo simulations for a realistic case study consisting of a direct beam from an X-ray tube with a W-anode, an acceleration voltage of 120 kV_p, and 6.8 mm Al filtering, detected using a CdTe sensor with a pixel size of 330 μm and a thickness of 1500 μm, and read out by a counting ASIC. The spectrum and the threshold energy were sampled with a step of 1 keV. The front-end electronics was assumed to introduce a noise of 1 keV rms. The pulse time was τ _P = 80 ns, and the retrigger time was τ _R = 100 ns. Incoming fluxes n were investigated up to 10⁸ ph/s/pix, which corresponds to roughly 10⁹ ph/s/mm² for the used pixel size, a value that can be considered the upper boundary in clinical CT applications. Figure 2 shows the energy spectrum as detected by the pixelated detector and used as input for both the model and Monte Carlo simulations. The spectrum of the original direct X-ray tube beam is also shown as a reference.

Figure 3 shows the comparison between the count rate curves predicted by the model and simulated with the Monte Carlo tool for incoming rates up to 10⁸ ph/s/pix and several values of threshold energies ranging from 5 keV to 140 keV. The rate curve corresponding to the total counts—intended as the ideal (no pileup and no electronic noise) rate obtained with a threshold energy of 0 keV—is also shown as a reference. For increasing threshold energies, the deviation of the recorded curves from one of the total counts, even at low incoming rates, becomes larger as the obvious consequence of losing an increasing fraction of low-energy events. On the other hand, all the recorded curves saturate for increasing incoming rates to the limit 1/τ _R , set by the retrigger mechanism. As a quantitative measure of the agreement between predicted and simulated results, we chose the L ² relative error norm (L²REN), which is defined as follows: L 2 R E N = ∑ i = 1 NoP x i sim − x i ana 2 ∑ i = 1 NoP x i sim 2 , (15) where NoP is the number of analyzed points. The values are reported in Table 1. An excellent match on the sub-percent level is obtained up to the last threshold energy (90 keV in this case), whose value lies below the highest energy actually contained in the input spectrum, i.e., 120 keV. For threshold energies higher than that, only pure piled-up events are detected, and the increase in L²REN seems to suggest that the accuracy of the approximations used in the derivation of the model (in particular, the term l ₄) begins, in these conditions, to fade.

Figure 4 (top) shows the comparison between the recorded integral energy spectra predicted by the model and simulated with the Monte Carlo tool as a function of the threshold energy and for increasing values of incoming rates ranging from n = 10⁵ ph/s/pix to n = 10⁸ ph/s/pix. Figure 4 (bottom) shows the corresponding deviations expressed in percentage. Looking at the integral spectra, it is possible to observe how, for increasing incoming rates, events of energy higher than the highest one contained in the input spectrum, i.e., 120 keV, and, therefore due to pure pileup effects, are detected with increasing probability. At the same time, the integral spectra tend to the limit set by the retrigger mechanism 1/τ _R , in a monotonic (i.e., non-paralyzable) way at all values of threshold energy. The relative deviation is in the order of a few percent up to the notable value of 120 keV (between −8% and +4%), a result that we can consider excellent considering the several orders of magnitudes spanned by the spectra across this threshold energy range. For threshold energies higher than this limit and high input rates (above 2⋅10⁷ ph/s/pix), the deviation slightly increases. As previously mentioned, this could be due to the approximations used in the derivation of the model starting to lose their validity. It is nevertheless worth noting that, in a realistic X-ray imaging/CT scenario, the threshold energy is not set above a certain fraction of the X-ray tube acceleration voltage. Therefore, a discrepancy in this extreme energy range is not of practical concern. For this reason, we limited the computation of the L²REN figure of merit to the threshold energy interval of 1 keV–120 keV. The values, reported in Table 2, show the agreement that lies at the (sub-)percent level.

Figure 5 shows the comparison between the recorded differential energy spectra obtained by numerical derivation of the integral spectra of Figure 4 (top) with respect to the threshold energy and the corresponding relative deviations. The L²REN metrics, computed in the threshold energy interval of 1–120 keV, are reported in Table 3. Although higher values of deviation as compared to the recorded integral spectra have to be expected due to the “amplification” caused by the derivative operation, they nevertheless lie below 10% up to incoming rates of 2⋅10⁷ ph/s/pix and at most 20% for higher incoming rates, which are anyway rather extreme as they basically lead to an almost complete count rate saturation.

Based on these results, we can conclude that the predictions of the presented analytical model faithfully represent the expectations obtained by numerical tools.

Real-world analog signals have shapes that significantly differ from the ideal rectangle. Apart from few other special cases (e.g., pulses with a triangular shape in [14] or with a decreasing exponential shape in [16]), a general analytical solution for arbitrary shapes remains, at the moment, more complex, if not out of sight.

To evaluate the impact of this non-ideality on the recorded spectra, we simulated a case study intentionally violating assumption ii. In particular, we used a more realistic pulse shape consisting of an asymmetric Gaussian-like function, as shown in Figure 6. The curve was created by using, as a rising edge, a Gaussian function with a standard deviation of 36 ns and, as a falling edge, a Gaussian function with a standard deviation of 134 ns. The overall FWHM is 80 ns, matching one of the ideal rectangles. The retrigger time was increased to τ R ′ = 150 ns in order to reduce spurious multiple counts at low threshold energies.

Figure 7 and Figure 8 show, respectively, the resulting integral and differential energy spectra (top) and corresponding relative deviations (bottom). Analytical data were recomputed with the new value of retrigger time. What emerges from the comparison is that while the overall trends are respected, significant deviations arise for values of incoming flux n ≥ 10⁷ ph/s/pix and threshold energies E_th < 10 keV. This is due to the fact that the model does not take into account the possibility that individual signals can have a time-over-threshold longer than the chosen retrigger time due to the longer pulse tail. The values of L²REN computed in the energy range of 10–120 keV for the recorded integral spectra vary in the range of 0.6%–6.4% (see Table 4), while for the recorded differential spectra, they monotonically increase from 8% to 106% with increasing incoming fluxes (see Table 5).

It is also legit to question whether the general discrepancy could be “compensated” by the use of an effective pulse time τ _Peff , as described in [6]. Unfortunately, as shown in the reference, τ _Peff would be a function of the threshold energy (to be empirically found, moreover) and, therefore, not compatible with the model presented in this work.