Radiation detection with pulse pile-up non-parametric decomposition

Abstract

A method of determining an energy spectrum of individual radiation quanta received in a radiation detector is disclosed. An energy spectrum sensitive statistic is computed from a time series of digital observations from the radiation detector, defining a mapping from a pulse amplitude density to the energy spectrum sensitive statistic. The density of pulse amplitudes is estimated by applying an inverse of the mapping to the energy spectrum sensitive statistic, thereby determining the energy spectrum. The statistic can be based on a first set of non-overlapping time intervals of constant length L, which is at least as long as a duration of a pulse, without taking into account an entirety of a pulse cluster, and a second set of non-overlapping time intervals of constant length L1, which is less than L, without taking into account the entirety of the pulse cluster. A method of estimating a count rate is also disclosed.

Description

Technical Field

[0001] A novel method is provided for estimating the energy distribution of radiation quanta, such as photons incident on a detector in a spectral system such as X-ray or gamma-ray spectroscopy. This method is particularly useful for count rate methods where pulse stacking is a problem. A key step in deriving the estimators in embodiments of the invention is to novelly reformulate the problem as a decomposition problem of a composite Poisson process. This method can be applied to any type of radiation detector that detects quanta or other particles of radiation, such as X-rays, gamma rays, or other photons, neutrons, atoms, molecules, or seismic pulses. Spectroscopic applications from such detectors are well known. Such applications are widely described in the prior art, including in international patent applications PCT / AU2005 / 001423, PCT / AU2009 / 000393, PCT / AU2009 / 000394, PCT / AU2009 / 000395, PCT / AU2009 / 001648, PCT / AU2012 / 000678, PCT / AU2014 / 050420, PCT / AU2015 / 050752, PCT / AU2017 / 050514 and PCT / AU2017 / 050512, each of which is incorporated herein by reference in its entirety for the purpose of describing the potential applications of the present invention and any other background material necessary for understanding the present invention. Background Technology

[0002] X-ray and gamma-ray spectroscopy underpins a wide range of scientific, industrial, and commercial processes. One goal of spectroscopy is to estimate the energy distribution of photons incident on a detector. From a signal processing perspective, the challenge lies in converting the pulse stream output by the detector into a histogram of the area under each pulse. Pulses are generated according to a Poisson distribution, the ratio of which corresponds to the intensity of the X-rays or gamma rays used to irradiate the sample. Increasing the intensity results in more pulses per second on average, and thus a shorter time before an accurate histogram is obtained. In applications such as airport baggage scanning, this translates directly to greater throughput. Pulse stacking occurs when two or more pulses overlap in the time domain. The incidence of pulse stacking increases with the count rate (the average number of pulses per second). This increases the difficulty in determining the number of pulses present and the area under each pulse. In the limiting case, the problem is ill-conditioned: if two pulses start substantially simultaneously, their superposition becomes indistinguishable from a single pulse. The response of an X-ray or gamma-ray detector to incident photons can be modeled as a pulse shape function. Superposition of convolutions with Dirac-delta pulses,

[0003] .

[0004] Arrival time The unknown, and forms a Poisson process. The arrival of each photon is modeled as time. The Dirac delta, and has a pulse shape response proportional to the photon energy that causes the detector. amplitude Amplitude It is its common probability density function Unknown identically distributed random variables Implementation of the pulse shape function. The shape of the pulse is determined by the geometry of the detector and the interaction of photons. In some systems, the variation in pulse shape is minimal and negligible, while in others (e.g., HPGe), the individual pulse shapes may differ significantly from each other [1]. It is assumed that all pulse shapes are causal, i.e., for t<0, Single peak, finite energy, and as It decays exponentially toward zero. The integral of the pulse shape function is normalized to units, i.e. This causes the area under the pulse to change from The observed signal consists of the detector output corrupted by noise, i.e.,

[0005] .

[0006] The mathematical goal of pulse stacking correction is to achieve the following given... In the case of uniform sampling and finite length version We perform estimation. We always assume the noise distribution... It has a known variance The photon arrival time follows a zero-mean Gaussian distribution. We also assume that the photon arrival time follows a homogeneous Poisson process with a known ratio. The detector response R consists of a noise process W that corrupts the signal, where... , and Represents the first of each sequence k There are elements, and among them .let S , R and W It corresponds to the uniform sampling time series of these signals.

[0007]

[0008] in .

[0009] Overview of Pulse Processing Methods: Over the decades, numerous methods have been proposed to address the problem of pulse stacking. These methods can be broadly categorized into two types: time-domain based and energy-domain based. A popular strategy is to attempt to detect when stacking occurs in the time domain and reject or compensate for the affected pulses. Early spectroscopic systems employed rejection-based methods along with matched filtering. The drawback of this approach is that as the stacking probability increases, an increasing proportion of pulses are rejected. The system quickly becomes unusable, thus setting an upper limit on the count rate [1]. With the increase in inexpensive computing power, the number of strategies for compensating or correcting stacking has grown. These include template fitting [2], baseline subtraction [3], adaptive filtering [4, 5], sparse regression [6, 7], etc. These methods all attempt to identify and compensate for stacking in the time domain and are generally best suited for systems with low pulse shape variations. The complexity of these methods increases with the pulse shape. The variability between them increases significantly. It can be shown that any method that attempts to characterize a single pulse will be affected by stacking. The best these methods can do is reduce the occurrence of stacking. Energy-based methods attempt to solve stacking based on the statistics of the pulse set rather than the individual pulse. They typically operate on a histogram of the estimated energy (the area under the pulse cluster). The earlier work of Wielopolski and Gardner [8] and more recent extensions of their ideas [9] operate primarily in the energy domain using set-based strategies. Trigano et al. [10, 11] used edge density to estimate the incident spectrum, which comes from the joint distribution of the statistical properties of the variable-length pulse clusters, where the start and end of each cluster are detected. This avoids the need to characterize individual pulses and is robust to changes in pulse shape. Ilhe et al.

[12] examined the exponential shot noise process, restricting the pulse shape to a simple exponential one for tractable results. Further work has been done

[13] to allow for a wider range of pulse shapes. In both cases, knowledge of the pulse shape is required, along with estimates of the characteristic function and its derivative. Summary of the Invention

[0010] We chose an energy-based stacking correction method to i) avoid the limitations associated with single-pulse detection

[14] and ii) be able to handle pulse shape variations without excessively increasing computational complexity. Instead of utilizing joint distribution [10, 11] or shot process

[12] methods, we recast the stacking problem as a “decomposition” problem of a composite Poisson process. A composite Poisson process is a discrete-time stochastic process in which each component consists of a sum of random numbers of independent and identically distributed random variables, where the number of random variables in each sum is a Poisson distribution

[15] . The “decomposition” of a composite Poisson process is the task of estimating the distribution of the random variables already extracted from it using random sums. Buchmann and

[16] The decomposition of compound Poisson processes is presented in the context of insurance claims and queuing theory. The decomposition of uniformly sampled compound Poisson processes has received some attention recently [16, 17, 18, 12, 19]. The context of these derivations often (reasonably) assumes that each event is detectable (i.e., there is no ambiguity about the number of events), or that the density estimator is conditioned on the fact that at least one event occurs in each observation

[20] . These assumptions have limited value in solving the problem of spectral packing.

[0011] Unconditional nonparametric decomposition of event detection has received relatively little attention in the literature. Gugushvili

[18] proposed a nonparametric, kernel-based estimator for decomposition problems in the presence of Gaussian noise. In embodiments of the present invention, the inventors have envisioned that once a method for selecting the kernel bandwidth is obtained, along with a method for transforming the observed detector output to fit a mathematical model, the estimator can be readily extended and applied to the reconstruction of spectral stacking problems.

[0012] According to a first general aspect of the invention, a method is provided for determining the energy spectrum of a single radiation quantum received in a radiation detector, the method comprising the steps of: (1) obtaining a time series of digital observations from the radiation detector, the time series including pulses corresponding to the detection of a single quantum; (2) calculating an energy spectrum sensitivity statistic based on the detector signal, the energy spectrum sensitivity statistic defining a mapping from pulse amplitude density to the energy spectrum sensitivity statistic; and (3) determining the energy spectrum by estimating the pulse amplitude density by applying the inverse of the mapping to the energy spectrum sensitivity statistic.

[0013] In an embodiment, the spectral sensitivity statistics can be based on the sum of digital observations over multiple time intervals, and the mapping can be defined using an approximate composite Poisson process, which can be enhanced by modeling noise. This mapping can be expressed as the relationship between the characteristic function of the amplitude, the spectral sensitivity statistics, and the modeling noise. The characteristic function of the spectral sensitivity statistics can be computed using a histogram of the sum of digital observations, to which an inverse Fourier transform is applied. The computation of the characteristic function of the amplitude may include the use of a low-pass filter.

[0014] In the first embodiment, the multiple time intervals are non-overlapping and have a constant length L, and each interval is chosen to contain zero or more approximately complete pulse clusters. This can be achieved by requiring the maximum value of the detector signal at the beginning and end of each time interval. In this embodiment, the composite Poisson process can be defined as the sum of the pulse amplitudes within each time interval. This mapping can be expressed as defined in equations (40) and (41), which can be enhanced by a windowing function.

[0015] In the second embodiment, the multiple intervals include a first set of non-overlapping time intervals of constant length L, selected without considering the pulse cluster as a whole, and a second set of non-overlapping time intervals of constant length L1 less than L, also selected without considering the pulse cluster as a whole. L is at least as long as the pulse duration, and preferably, L1 is less than the pulse duration. In this embodiment, the composite Poisson process can be defined as in Section 6. The mapping can be expressed as defined in Section 6. The second embodiment can utilize the processing and calculation of each set of time intervals defined for the set of time intervals as in the first embodiment.

[0016] In this embodiment, a data-driven strategy is used, which results in a near-optimal selection of the nuclear parameters that minimizes the integral squared error (ISE) of the estimated probability density function of the incident photon energy.

[0017] According to a second generalized aspect of the invention, a method is provided for estimating the count rate of a single radiation quantum received in a radiation detector, the method comprising the steps of: (1) obtaining a time series of digital observations from the radiation detector, the time series including pulses corresponding to the detection of a single quantum; (2) calculating a spectral sensitivity statistic based on the detector signal, the spectral sensitivity statistic using an interval of constant length L and constant length L1 as described above with respect to the first generalized aspect; (3) determining an estimate of the characteristic function of the composite Poisson process using formula (109); and (4) estimating the count rate based on the estimate of the characteristic function. Step (4) above can be achieved by fitting a curve using an optimization routine or some other means, estimating the DC offset of the logarithm of the estimated characteristic function, or fitting the curve to the logarithm of the estimated characteristic function.

[0018] The remainder of the application is organized as follows. Sections 3, 4, and 5 relate to a first embodiment of the first aspect of the invention. Section 3 provides preliminary background, defines notation, outlines the mathematical model, and gives the derivation of the estimator of the first embodiment, including modifications. Section 4 shows the performance of the modified estimator of the first embodiment on both simulated and experimental data, and discusses the results. Section 5 provides the conclusions of the first embodiment. Section 6 describes a second embodiment, referring to the first embodiment where relevant. Section 7 describes a second aspect of the invention, a novel method for estimating count rates.

[0019] 3. Derivation of the estimator in the first embodiment

[0020] Our general approach to solving the stacking problem is based on the following strategy: i) from Obtain statistics sensitive to the energy distribution of incident photons, and use... The methods are: 1) estimating the statistics using a finite-length sampled version of the observations; 2) obtaining a mapping from the incident photon energy density to the statistical properties of the observed statistics; and 3) estimating the incident photon energy density by inverting this mapping. Section 3.1 describes our choice of statistics. Section 3.2 assumes that these statistics (approximately) have the same distribution as the composite Poisson process. Section 3.3 introduces the decomposition technique used to recover the energy spectrum from these statistics. It is based on the decomposition algorithm in

[18] , but further developed to achieve near-optimal performance in terms of integral squared error. The general theme followed by our approach to the stacking problem is to find some of the statistical properties that are sensitive to the underlying energy spectrum. Statistics, from These statistics are estimated in a finite-time sampling version, and then the mapping describing the statistical properties of these statistics is inverted given the base energy spectrum to produce an estimate of the energy spectrum.

[0021] 3.1 Selection of Statistics

[0022] We aim to obtain an estimate of the photon energy from the observed signal given in (2). In a typical modern spectroscopic system, the detector output... Uniform sampling by ADC. Without loss of generality, we assume that the raw observations available to the algorithm are... Since identifying a single pulse can be difficult, we instead look for fixed-length pulse clusters containing zero or more pulse groups. The intervals. More precisely, we define these intervals as... ,in

[0023] .

[0024] here, It was chosen as a trade-off between the error in energy estimation and the probability of creating the interval. The value of should be small enough to ensure an acceptablely low error in estimating the total photon energy arriving in each interval, but large enough relative to the noise variance to ensure a large number of intervals are obtained. Although the probability of partitioning the observed data into intervals approaches zero as the count rate approaches infinity, this method falls into a deadlock at higher count rates compared to a stacking rejection strategy based on a single pulse, because quasi-numerous photons are stacked in each interval. Section 4.2 describes L and for real data. The choice is as follows: Each interval contains an unknown, random number of pulses, and may contain zero pulses.

[0025] We use the original observations from the sampling to estimate the interval. Total photon energy Since the area under each pulse is equal to the photon energy defined in (1). Proportionate, we let

[0026] .

[0027] Assuming each interval The number of photons arriving, the energy of each arriving photon, and the detector output noise are random and independent of other intervals. For pulse shapes with exponential decay, the energy of a small number of photons arriving in one interval can be recorded in the next interval. Leakage is related to... Proportional and for sufficiently small Negligible. Therefore, the estimate is... This can be viewed as the realization of a weakly correlated stationary process, where each estimate is identically distributed according to the random variable X. This relationship is... Figure 1 The typical pulse shape is used to illustrate the noise-free case.

[0028] 3.2 Using the approximation of the composite Poisson process

[0029] In this section, we discuss... The distribution of X is described. We will then invert this in Section 3.3 to obtain the density. The estimator. Using (9), (2), (1) and It is a fact of cause and effect, we have

[0030] .

[0031] The following proof simplifies to

[0032]

[0033] in

[0034] and .

[0035] and Both are iid sequences of random variables. We use Y and Z to represent their distributions. The distribution of Z is completely derived from... The distribution is determined, assuming The distribution has zero mean and a known variance. The Gaussian distribution. Furthermore, Y is a composite Poisson process because the number of terms in the sum (the number of photon arrivals in an interval of length L) has Poisson statistics. Equations (11)-(13) are proved as follows. The first term in (10) represents leakage from the earlier intervals and is approximately zero. For Gaussian noise, by relating... Performing a Taylor expansion, this is easily shown.

[0036]

[0037] .

[0038] Therefore, the probability that some energy belonging to previous intervals will be included in the current estimate is finite but small. In fact, this contribution is sufficiently small. The noise is considerable. The third term in (10) is zero because It is causal. The second term in (10) can be written as

[0039]

[0040] Here, we assume the pulse shape Smooth enough to make It approximates all at intervals. The total energy of the photons arriving inside. Let Specify at interval The number of photons arriving. We assume... It has a ratio parameter The realization of the homogeneous Poisson process, in which Let the expected number of photons per interval of length L be used. We will then assume that (11) holds exactly, and write...

[0041] .

[0042] Finally, we will Written

[0043]

[0044] Here, we assume that Z has a known variance. In this subsection, we model the statistics from Section 3.1 using a compound Poisson process. This allows us to derive the density from observable quantities. The estimate. In the interval The number of photons arriving inside is what we specify The total energy in interval Y is a Poisson random variable. This total energy can be modeled as a composite Poisson process, i.e.,

[0045]

[0046] in, It is the index of the arrival time of the first photon in the interval, the arrival times are assumed to be ordered, and it represents the photon energy. It has a density function Independent realizations of random variable A. Formed with ratio parameters A homogeneous Poisson process. Poisson's law. Expressed in terms of the expected number of photons per interval of length L.

[0047] .

[0048] The relationship between the implementation of Y and the response of the sampling detector is in Figure 1 The diagram is shown in the image. By substituting (2) into (9), we can use... To approximate observation ,

[0049]

[0050] in It is a realization of the unobservable random variable Y, which represents the photon energy in the interval of the discrete-time detector response.

[0051]

[0052] in Z is a realization of an independent random variable, representing the sampling process and estimation. The error in Z. We assume that Z has a known variance. Under these definitions of X and Y, the number of intervals that can be found in the detector output of a finite length is a random variable N. At high count rates, the method falls into paralysis because the probability of being able to partition the observation data into intervals approaches zero. The onset of paralysis occurs at higher count rates compared to the stacking rejection-based strategy because a quasi-many photons are stacked in each interval. Assume that the time series defined in (3)-(6) has been uniformly sampled. Without loss of generality, assume that the unit sampling interval is from Beginning, that is, , Let R be a discrete-time random process, and let represent the sampling detector response of (1). It is a discrete-time random process, and its components This represents the total energy of photons arriving during a fixed time interval. The composite Poisson process can be used to model Y, i.e.,

[0053]

[0054] in, They are independent Poisson random variables, and It has a density function . Independent and identically distributed random variables. Forming with ratio parameters The homogeneous Poisson process. Process Y is not directly observable. Assume the pulse shape... It has limited support. Let Let A be an indicator function for set A. Let the pulse length... Depend on Give it. Let. It is a discrete-time random process, represented by the detector output of the observation given in (2). It consists of the detector response R corrupted by the noisy process W. Without loss of generality, we assume a unit sampling interval. Based on the observation S, we form the process X, where

[0055]

[0056] And among them, It comes from the known variance The random variables of the independent noise process. When we let the pulse shape in (1) At that time, a simple test theoretical model was obtained, and in this case, we let ,and N Only the sample length K .from S get It's more complex with real data. In that case, we'll discuss the process. S The partition is of length L Non-overlapping blocks, where Poisson's ratio L Represented using photons in each block. The start of each block is chosen. This ensures that the total energy of any pulse is completely contained within the block it arrives at.

[0057] .

[0058] Figure 1 Show We let

[0059]

[0060] in yes The estimate. Section 4.2 describes the actual data. L and The choice. In Under this definition, for a given sample length K , Y The number of components becomes a random variable. At high count rates, the method fails because the probability of creating a block approaches zero. Compared to stacking rejection-based strategies, the failure occurs at even higher count rates because quasi-numerous photons accumulate within each block. It is a discrete-time random process, and its components Given by the following formula

[0061]

[0062] in It is a constant, chosen such that h and d are small thresholds close to zero. Therefore, the random variable This represents the total photon energy arriving during a fixed time interval of length L. The value of d ensures that the signal associated with photon arrival is very small at the beginning and end of each interval. This is in... Figure 1 The diagram in the middle shows that a compound Poisson process can be used to model Y, i.e.,

[0063]

[0064]

[0065] in It has a ratio parameter The homogeneous Poisson process, and It has a density function Let be independent and identically distributed random variables. It is a discrete-time random process, represented by the sampling detector output given in (2). It consists of the detector response R, which is corrupted by the noisy process W. The process Y is not directly observable. Using (2), (25), and (32), we obtain the process Y through the process Y. Modeling observations, that is,

[0066]

[0067] Where Z is the known variance The noise process. In a given observation All random variables involved in the modeling They are all assumed to be independent. Let There are N independent and identically distributed observations. Let... become The collection. Let the corresponding feature function be... .

[0068] 3.3 Basic Forms of Estimators

[0069] We seek to invert the mapping from the distribution of photon energy A to the distribution of X. Our strategy is to first, based on... Obtain the characteristic function of X, then, assuming the count rate and noise characteristics are known, invert this mapping. Let... yes The characteristic function. As is well known

[15] for ratios The complex Poisson process Y,

[0070]

[0071] And since X = Y + Z, therefore

[0072] .

[0073] Given observations We can form an empirical estimate of the characteristic function of X. Treating this as the true characteristic function, we can invert (40) and (41) to obtain the characteristic function of A, and then perform a Fourier transform to find the amplitude spectrum. Specifically, using (40) and (41) and leveraging the assumption that Z is a Gaussian distribution, we can ensure... It will be non-zero. We let The curve described by the following formula

[0074] .

[0075] Temporary assumption After taking the difference logarithm of (43) and rearranging it, we have

[0076] .

[0077] Ideally, it can be recovered by performing a Fourier transform.

[0078] .

[0079] The basic form of our proposed estimator is given in (88) and derived from (45) through a series of steps. First, based on the data... Make an estimate (step 1). Simply replace the estimate in (42). No production The best estimate of ISE. The approximate ISE is derived from... The approximate estimate of the error distribution is obtained (step 2). Then, we determine a reasonable windowing function. (In step 3), and by the following formula Make an estimate

[0080] .

[0081] Windowing function Designed to be minimized Based on (44), (45) and (46), our... The approximation between the estimates of ISE, but in which (44) Replaced by (46). A similar idea was used to replace (44). Estimation: (In step 4) Find the weighting function This allows us to replace (45) with the following formula.

[0082]

[0083] The result is a comparison with the use of unweighted estimates. Better f A Estimate. Finally, the weighted function. The modification (in step 5) to account for the integral in (45) must actually be replaced by a finite sum. The following subsections elaborate on these five steps.

[0084] 3.4 Estimation

[0085] need The estimate is to estimate In this section, we define the histogram model and describe our approach based on... histogram of values The estimation. Assume that N intervals (and corresponding intervals) are obtained from a finite-length data sample. Value). Although empirical characteristic function

[0086]

[0087] It provides a consistent, asymptotically normal estimator of the characteristic function

[21] , but it has the disadvantage that the estimator becomes less reliable as the number of data points N and the number of evaluation points required increases. The computational burden increases rapidly with this approach. Instead, we use a histogram-based estimator, which has a lower computational burden. Assume the histogram of observed X values ​​is formed by... vector The count in the m-th bin is given by the following formula.

[0088] .

[0089] All bins in the histogram have the same width. The bin width is determined by... The value is selected based on the quantity. Since choosing different warehouse widths has the effect of simply scaling... Therefore, without loss of generality, we assume that the warehouse width is uniform. Warehouses are evenly distributed between non-negative and negative data values. The number of histogram warehouses is 2. M This affects the estimator in various ways, as discussed in later subsections. Currently, we assume 2... M Large enough to ensure the histogram includes all The value is sufficient. We achieve this by forming a scaling factor. histogram of values We perform an estimation and obtain the inverse discrete Fourier transform, i.e.,

[0090] .

[0091] This is a close approximation of the empirical characteristic function, but in which... The item is rounded to the nearest histogram center (and u Reduced to 1 / (times). Item It simply counts the number of rounded items with the same value. Clearly, this can be achieved using the Fast Fourier Transform (FFT) at certain discrete points. Efficiently evaluate the function.

[0092] 3.5 Error distribution

[0093] In (46) and (47) the filters and The design depends on Error statistics between the true and false characteristic functions. In this subsection, we define and describe the characteristics of these errors. We assume the density function... Smooth enough (i.e., Furthermore, the width of the histogram bin is small enough (relative to the standard deviation of additive noise Z) to allow for the transfer of data from the histogram bin to the histogram bin. The error introduced by rounding the value to the center of each histogram bin is approximately uniformly distributed across each bin, has zero mean, and exhibits minimal peak diffusion relative to Z. In other words, by The error sources caused by the loading of the value are considered negligible. This is because both the statistical properties of Poisson counts and the expected count in each bin are non-integers. Therefore, there is a difference between the number of counts observed in any given histogram bin and the expected number of counts in that bin. In our model, we combine these two error sources and call it “histogram noise”. We emphasize that this noise is different from the additive noise Z modeled in (11), which causes peak diffusion in the histogram. Let the probability that the realization of X falls into the m-th bin be given by...

[0094]

[0095] Let the normalized histogram error in the m-th warehouse be... It is the count observed in the m-th warehouse. and expected count The difference between them is relative to the total count N in the histogram, i.e.

[0096] .

[0097] Using (50), (51), and (52), we have

[0098] .

[0099] If the histogram is modeled as a Poisson vector, then it shows

[0100] .

[0101] Since the characteristics of histogram noise can be described by the total number of observation intervals N, the impact of using observation data of finite length can be mitigated by incorporating this information into... and This should be taken into account in the design.

[0102] 3.6 Estimation

[0103] get After that, the next task is to estimate. It is not used in (42). replace Instead, we use (46) as the estimator, which requires us to choose a window function. In this section, we attempt to find a function that approximates the optimal one. When considering When considering the error distribution in (46), the windowing function of the lowest ISE estimator in the form given in (46) is applied. for

[0104]

[0105] in express The real part. We cannot calculate ,because Since it's unknown, we instead try to find approximate values. We let...

[0106] .

[0107] This is achieved by considering the function. and The relative error between them is used to prove this, where

[0108] .

[0109] The magnitude of the relative error is given by the following formula.

[0110] .

[0111] because We see that the right side of (63) is The relative error is maximum at -1. Therefore, the relative error is bounded by the following formula.

[0112]

[0113] when When I was very young, or when At that time, it proved this approximation. Furthermore, we note that the above bound is quite conservative. The distribution of photon energy in a spectral system can generally be modeled as the sum of K Gaussian peaks, where the k-th peak has a position... and scale ,Right now,

[0114]

[0115] in .

[0116] Therefore, the characteristic function will have the following form

[0117]

[0118] That is, for some ,like Oscillations within the decaying envelope. (64) gives a fairly conservative upper bound because for most... value, At most evaluation points across the energy spectrum, the approximation error will be significantly smaller. (Selection) Then we can use (46) to form a pair The estimation. The windowing function reduces the impact of histogram noise caused by a finite number of data samples. For large... The effect of windowing is negligible, and the estimator is essentially the same as using (42) directly. However, in the region where the following equation holds...

[0119]

[0120] Adding windows becomes important and works to improve our lives. The estimation bound is based on the fact that the noise Z is Gaussian distributed (therefore...). And therefore ), and because We see

[0121] .

[0122] This ensures that the variables of the distinguishing logarithm in (47) remain finite, even though .

[0123] 3.7 Estimation

[0124] Once you get We will continue to use (47) pairs To perform the estimation, another windowing function is needed. In this subsection, we find a near-ISE optimal method for estimating... function We start by defining functions. To begin with easy labeling

[0125] .

[0126] when When ISE is minimized, the optimal filter is found. Given by the following formula

[0127] .

[0128] Secondly, due to , and Since the optimal filter is unknown, we cannot compute it using (73)-(74). Instead, we make the following observations to obtain an approximation of the optimal filter for the ISE.

[0129] 3.7.1 Initial Observation

[0130] Only estimates Keep close The optimal filter keeps the true value close to unit 1. For u The smaller the value, the more likely the situation will remain the same, because

[0131]

[0132] For small u, (76).

[0133] Furthermore, equation (73) shows that if ,So ,therefore For larger u values, when The value becomes equal to or less than At that time, the estimated quantity Dominated by noise and no longer providing useful information. Estimate. In extreme cases, ,therefore And therefore

[0134] .

[0135] window These regions should be excluded from the estimate because the bias introduced would be less than the variance of the unfiltered noise. Unfortunately, before reaching that boundary condition, The estimate may be severely degraded, so (77) is not particularly helpful. A more useful method for detecting when noise begins to dominate is as follows.

[0136] 3.7.2 Filter Design Function

[0137] Further manipulation of (67) shows that, for a typical spectral system, The value will have the following form

[0138]

[0139] That is: based on peak width The average component of attenuation, and based on the position of the spectral peaks. The oscillating component changes more rapidly and decays faster. (In the design window) At that time, we are interested in attenuation. The area in, where That is, where the signal power is less than the histogram noise, which is in The estimated period is achieved by removing And thus enhanced. In order to obtain The estimate will be based on the low-pass Gaussian filter. and Perform convolution to decay the division Beyond the slowly changing large-scale features. We represent that...

[0140] .

[0141] We see With scale parameters The Rayleigh distribution. Therefore

[0142] .

[0143] As is well known, Rayleigh distributed random variables The cumulative distribution function is given by the following equation.

[0144]

[0145] Therefore, in order to help calculate windows We will use functions

[0146]

[0147] To control The shape. This function. We can provide estimates The approximation in (84) is a confident indication of how much more signal energy is contained than noise energy. This approximation stems from the fact that... Also affected by noise Random variables with slight effects. Sometimes—especially for larger ones. Values—Histogram noise can lead to sufficiently large values. This leads to a false sense of confidence and potentially allows noise-driven results to undermine the outcome. The estimation. To overcome this problem, the function was modified to be unimodal in u.

[0148] .

[0149] This modification is proven on the assumption that Gaussian noise causes exist The decrease in [the value / value] is expected. Therefore, we anticipate [the decrease / decrease]. In Add to the middle. If we ignore the fact that... The peak position in the middle caused by Local oscillations in the smoothed material are then resolved. The approximate envelope will not be Add to. Equation (74) indicates the optimal window has Therefore, the overall window shape will be in the form of... The value decreases. Therefore, if at a certain point... The estimated characteristic function in the region (where the signal-to-noise ratio is high) has been determined, and the window value should be... So, refusing to be in the area (In areas where the signal-to-noise ratio will be even worse) The suggestion is reasonable. Use the following knowledge: For small It should be close to one unit, for large Approaching 0, and should "roll off" as the signal-to-noise ratio decreases—we consider two potential windowing functions as countermeasures. Approximate to .

[0150] 3.7.3 Rectangular Window

[0151] The indicator function provides a very simple windowing function.

[0152] .

[0153] threshold The cutoff point has been determined, and it can be manually selected as desired (e.g.) Once the threshold is selected, the estimator exhibits similar ISE performance regardless of the peak position in the incident spectrum. It does not require the user to depend on the window width selected from the incident spectrum (Gugushvili

[18] proposed a nonparametric estimator for general decomposition problems using a rectangular windowing scheme. It requires manual selection of the window width, which varies with...). (Change), but rather by data via. Automatic window width selection. While simplicity is a major advantage of rectangular windows, abrupt changes in the region provide a poor model for the roll-off region of the optimal filter. A second filter shape attempts to improve this.

[0154] 3.7.4 Logistic Window

[0155] A window based on the logistic function attempts to model a smoother roll-off. It is given by the following equation:

[0156]

[0157] in Again, the received signal energy is greater than estimated. The assumed threshold for noise energy. The filter roll-off ratio near the threshold region is determined by... Control. This provides a smoother transition area than a rectangular window, reducing... The final estimate of the Gibbs oscillation. Once again, although the parameters They are manually selected, but they are for It has a much smaller dependency on [the specific value] and can be used to provide near-optimal filtering for a wide range of incident spectra. A typical value used is [value]. The performance of the rectangular and logistic window functions is compared in Section 4.

[0158] 3.8 Estimation

[0159] Window functions were designed Therefore, an estimator was designed. The final task is to estimate using the inverse Fourier transform. This section describes several issues that arise in the digital implementation. First, regarding the entire solid line... and Numerical evaluation is not feasible. Instead, we estimate it at discrete points over finite intervals. The finite intervals are chosen to be large enough that the error due to excluding signal values ​​outside the intervals is quite small. This is for... It is proven to be a Gaussian mixture because for some , and The value will be as follows Attenuation. The Fast Fourier Transform (FFT) is used to evaluate at discrete points. And thus the assessment was determined. and Similarly, FFT is used to evaluate at discrete points. The final estimate. For FFT to be used, the signal outside the interval should be small enough to reduce aliasing. The evaluation points also need to be dense enough to avoid evaluation errors. Any "phase wrapping" blurring during time. Both of these objectives can be achieved by increasing the number of bins in the histogram by 2M (zero-fill) until a sufficient number of bins is obtained. As M increases, The increased sampling density allows for the detection and management of phase winding. A larger M also allows for aliasing (caused by...). The Gaussian tails (caused by the curve) are negligible. Typically, the value of M is chosen as the least power of two, large enough that the non-zero values ​​of the histogram are restricted to the "lower half" indices, i.e., Secondly, the distinguishing logarithm in (47) is undefined, if If so. Based on the data... When making the estimate, the probability that the estimate will be zero exists but is not zero. In this case, the difference logarithm in (47) is undefined, and the technique fails. With Increase, Reduced and possibly close .when and When they have similar values, (and therefore) The probability of a value close to zero may become significant. (Filter) It should roll down faster than near This is to reduce the potential impact on the estimation. Ideally, noise should be minimized. In the region close to zero It should be zero. Gugushvili has shown

[18] that for rectangular windows, as the dataset length increases... The probability of failing to invert the value is close to zero.

[0160] 3.9 Discrete Symbols

[0161] We digress for a moment to introduce additional notation. Throughout the rest of the article, bold text will be used to indicate the 2M×1 vector corresponding to the discrete-sampled version of the named function, for example, This represents a 2m×1 vector, whose value is determined by the point... Characteristic function of evaluation Given. The square bracket symbol [k] is used to index a specific element in the vector, for example, have We also use negative indices to access vector elements, similar to the Python programming language. Negative indices should be interpreted relative to the length of the vector, i.e., ... Refers to the last element in the vector (which is equivalent to) ).

[0162] 3.10 Overview of Estimators

[0163] The estimation process we used can be summarized in the following steps.

[0164] 1. Use (8) to partition the sampled time series into intervals.

[0165] 2. Calculate the value of each interval according to (9). value.

[0166] 3. According to Value generation histogram .

[0167] 4. Calculate using inverse FFT To efficiently evaluate at each sampling point (50).

[0168] 5. Calculate at appropriate points. and .

[0169] 6. Used via (46) and calculate .

[0170] 7. Calculate , The low-pass filter version.

[0171] 8. Calculated via (83) and (85) .

[0172] 9. Use Calculate with (86) or (87) .

[0173] 10. Used via (47) and calculate .if Any element of is zero, and If the corresponding element is not zero, the estimation fails because the distinguishing logarithm is undefined.

[0174] 11. Use according to the following formula FFT calculation

[0175] .

[0176] 3.11 Performance Metrics

[0177] The performance of the estimator is measured using the integral squared of the error (ISE). ISE measures the global fit of the estimated density.

[0178]

[0179] The discrete ISE metric is given by the following formula.

[0180]

[0181] in It is a 2M×1 vector whose elements contain the probability mass of each histogram region, i.e.

[0182] .

[0183] The vector This represents the corresponding estimated probability mass vector.

[0184] 4. Numerical Results of the First Embodiment

[0185] Experiments were conducted using both simulated and real data.

[0186] 4.1 Simulation

[0187] The ideal density used by Trigano et al.

[11] for these simulations. It consists of a mixture of six Gaussian distributions and one gamma distribution to simulate the Compton background. The mixture density is given by the following equation.

[0188]

[0189] in It has a mean and variance The density of the normal distribution. The density of the gamma distribution is determined by... Given, the density is sampled at 8192 equally spaced integer points to produce a discrete vector of probability mass. Perform an FFT to obtain , The sample vector of the value.

[0190] A specific count rate was selected for the experiment. , corresponding to the expected number of events per observation interval. The expected stacking density is obtained via (40), i.e., the discrete vector. Scaling Take the power of 1, then scale it. And finally apply FFT.

[0191] .

[0192] Equation (93) is convolved with a Gaussian distribution to simulate the effect of erasing the noise Z in the observed spectrum.

[0193] .

[0194] This represents the expected density of the observed spectrum, including packing and additive noise. Histograms of the observations were created using random variables distributed according to (94). Experiments were conducted by... And parameterization, in which and For each parameter pair Create one thousand observation histograms. Probability mass vector. The estimate was made using (88), where both (86) and (87) were used. For both window shapes, use The threshold, and for logistic shapes, use The threshold. Record each estimate. and the true vector The discrete ISE measure of the error between the two sides. To compare with the asymptotic bandwidth results, a rectangular window is used for estimation, the bandwidth of which is based on condition 1.3 specified by Gugushvili in

[18] – i.e. ,in —Selected. We emphasize Gugushvili's filters. Do not use (87) Confusion. The asymptotic bandwidth standard is achieved by using the following formula.

[0195]

[0196] in

[0197] Experimented with Gugushvilli The three values, namely .

[0198] Estimation was also performed using a rectangular filter (95) with various values. The fixed bandwidth. Finally, time series data is created according to (1), with an ideal rectangular pulse shape and 10 7 Each pulse has an energy distribution according to (92). The pulse length and count rate are chosen to give the Poisson's law. The algorithm described by Trigano et al.

[11] is used to estimate the fundamental amplitude density from a two-dimensional histogram containing 32 × 1024 (duration × energy) bins—the choice of bins reportedly yields optimal accuracy and reasonable execution time. The performance and processing time of the core algorithm are recorded for comparison with our proposed algorithm. Figure 2 The data-driven logistic filter pair is described with parameter pairs ( Typical estimates made from experiments The actual vectors were also drawn. (Thin solid line) and the observed histogram (The lower curve contains some noise). The packing peaks are clearly visible in the observed histogram. Although the estimated density is affected by ringing (due to the Gibbs phenomenon), it estimates the true density in another way and corrects for the packing present in the observed histogram. Figure 3 Depicting in relation to Figure 2 Typical estimations are performed at the same operation point, but the estimator has a rectangular filter, where (96) and Select bandwidth. This corresponds to Figure 6 The operating region in which the fixed bandwidth filter ( The performance of the ) is close to that of the data-driven filter. Clearly, while also correcting for stacking, the resulting estimate contains more noise. Figure 4 The distribution density of the ISE metric as a function of sample count is shown when using a rectangular filter and various fixed bandwidths. A line is plotted between the distribution mean (MISE) to aid visualization. The results for the data-driven rectangular filter (86) are also plotted and connected with a thicker curve. This clearly illustrates the weakness of fixed-bandwidth filtering. For any fixed bandwidth, the ISE decreases with increasing sample count, eventually asymptotically approaching the point where bias becomes the dominant source of error. At this point (which is noise- and bandwidth-dependent), the ISE remains largely unchanged despite increasing sample count. Fixed bandwidth excludes the use of some estimates in the final calculation. Even when they have a high signal-to-noise ratio (SNR). Figure 4 The results for our proposed data-driven bandwidth selection rectangular filter are also shown. The curve lies near the inflection point of each fixed bandwidth curve. This indicates that across the sample count range, the bandwidth selected for the data-driven rectangular filter is close to the optimal bandwidth value (for a rectangular filter). Figures 5-7 It shows the three count rates Below, the distribution density of the ISE metric is used as a function of the estimated total number N in each histogram. The MISE curves for logistic and rectangular filters are lower than those obtained using the bandwidth given in (96) for most regions of interest in the application. Various regions exist where non-data-driven bandwidth ( The proposed algorithm delivers similar performance to the data-driven bandwidth; however, this performance is not maintained across the entire sample count range. The logistic filter shape exhibits slightly better performance than the rectangular filter shape, although the difference between the two filters is relatively small for the ISE metric. Table 1 compares the results between the proposed algorithm and the algorithm recently described in

[11] . The ISE of the two methods in the test ( The operation points are similar to those in the previous algorithm, however, our proposed algorithm requires significantly less computation.

[0199] Table 1: Comparison with the algorithm described in

[11]

[0200]

[0201] 4.2 Real Data

[0202] The estimator was applied to real data to assess its usefulness in practical applications. The threshold found in (8) Selected as additive noise The standard deviation is approximately half of the value. This ensures a reasonably high probability of creating the interval while also ensuring low error in the interval energy estimation. The value of the interval length L is chosen to be approximately four times the typical pulse "length," i.e., the interval length... Four times that of manganese. An energy histogram was obtained from the manganese sample, where the photon flux rate was nominally 102 per second. 5 One event. The observed histogram shows a slight negative skew in the shape of the main peak, indicating that complex noise sources are affecting the system. This is in Figure 8 The noise is almost invisible. It is modeled as a bimodal Gaussian mixture rather than a single Gaussian peak. A simple least-squares optimization routine is used to fit the bimodal Gaussian parameters. The noise peak is located around index zero. (Manual selection) A suitable value is found. A logistic filter with data-driven bandwidth is used to estimate the true density. Figure 8 The observed and estimated probability mass vectors are plotted. The main peak (450-600) is enhanced, while the packing, though not completely eliminated, has decayed. The first-order packing peak has been reduced. The peak-packing ratio (the ratio of the main peak height to the first packing peak height) has increased from approximately 6 to approximately 120. These improvements are comparable to other state-of-the-art systems (e.g.,

[11] ). The estimator fails to fully account for several possible reasons for the presence of packing. The accuracy of the estimator depends on the correct modeling of the Gaussian noise peaks.

[0203] The bimodal Gaussian mixture model modeled the noise peaks, ensuring the maximum error was less than 1% of the noise density peak. Given that residual packing peaks in the estimated spectrum are below 1% of the main peak, the sensitivity of the estimator to errors in noise modeling may contribute to this in some parts. A second reason for unresolved packing may be uncertainty in the observed spectrum estimation. Several residual packing peaks are relatively close to the bottom of the observed histogram. These residual peaks may simply be artifacts caused by noise in the estimator. Finally, the mathematical model may be an overly simplistic approximation of the observed spectrum. The detection process includes many second-order effects not included in the model (e.g., ballistic deficit, supply charge depletion, correlated noise, nonlinearity, etc.). These minor effects may limit the accuracy of the packing correction estimator.

[0204] 5. Overview of the First Embodiment

[0205] We adopted the estimator proposed by Gugushvili

[18] for decomposition under Gaussian noise and adapted it for pulse stacking correction in X-ray spectra. We proposed an easily implemented data-driven bandwidth selection mechanism and applied it to a wide range of sample counts of interest across the spectrum (10). 4 ~10 9 The count range provides a significant reduction in ISE / MISE. Data-driven rectangular bandwidth selection is near optimal (for rectangular filters) and outperforms bandwidth selection based on asymptotic results or fixed bandwidth within the range of interest.

[0206] While the initial results appear promising, further work is needed to improve the performance of the practical implementation. The estimate still contains the "ringing" artifact associated with the Gibbs phenomenon. Additional filter shapes attempt to reduce this, and other shapes that are closer to the MSE optimum exist.

[0207] 6 Second Embodiment

[0208] This section provides an overview of the energy spectrum estimator of the second embodiment. The second embodiment addresses the requirement of the first embodiment to approximately include the entire cluster in each interval. In the second embodiment, the entire data sequence can be used if desired, and overlap is compensated for by introducing two different interval lengths, L and L1.

[0209] We need to include some additional terms that were not mentioned in the first embodiment. In particular The energy spectrum estimate is based on

[0210] .

[0211] filter The introduction of this allows us to address several implementation problems that arise. The estimation process we use can be summarized in the following steps.

[0212] 1. Divide the sampled time series into fixed-length intervals. .

[0213] 2. According to Calculate each interval value.

[0214] 3. From Value generation histogram .

[0215] 4. Use Inverse FFT calculation .

[0216] 5. Partition the sampled time series into different interval sets with length L1, and then perform similar calculations to obtain... .

[0217] 6. Calculation and .

[0218] 7. Use and calculate .

[0219]

[0220] 8. Calculate , The low-pass filter version.

[0221] 9. Calculate .

[0222] 10. Use calculate .

[0223] 11. Use and calculate .if Any element of is zero, and If the corresponding element is not zero, the estimation fails because the distinguishing logarithm is undefined.

[0224] 12. Use The FFT is calculated according to the following formula.

[0225] .

[0226] 6.1 Algorithm Details

[0227] Partition the detector output stream into lengths L The non-overlapping interval set, i.e. .let It is the sum of the detector output samples in the j-th interval, that is,

[0228] .

[0229] Assumption L If the pulse length is greater than the j-th interval, the j-th interval may contain both "complete" pulses and pulses truncated at the end of the interval. This can be illustrated as follows: By us The energy of a "complete" pulse, which we use The truncated pulse and noise are represented It is composed of the superposition of energy.

[0230] Partition the detector output stream into a second non-overlapping interval set. ,in .let Given by the following formula

[0231] .

[0232] If L1 is chosen to be slightly smaller than the pulse length, then The item will not contain "complete" pulses, but only truncated pulses. and noise The energy is composed of a superposition. The number of truncated pulses within any interval follows a Poisson distribution. We have

[0233] .

[0234] We can divide the interval The total energy in the pulse is decomposed into energy contributions from the truncated pulse. and from completely contained in the interval Energy contribution of the pulse ,Right now,

[0235]

[0236] in, This indicates that the pulse is completely contained within the interval (length is...). Noise in the region of ) and This indicates that the pulse is truncated (length is...) Noise in the region of ). Therefore,

[0237] .

[0238] By combining (103) and (105), we have

[0239] .

[0240] Rearranged

[0241] .

[0242] We can use something similar to our estimate Estimate by way of or some other method For example, via empirical characteristic functions or through the analysis of... Perform an FFT on the normalized histogram of values.

[0243] When performing a decomposition operation, the reduced interval length Poisson's ratio Used to account for the complex Poisson process that occurs thereon. The sub-spaces.

[0244] 6.2 Visualization of Internal Quantities

[0245] To help readers understand, Figure 9 The various quantities obtained during the estimation process are plotted. The blue curve above (with a value of approximately 0.3 at zero) plots... The estimated characteristic function of the observed spectrum is shown. The brown curve is used to illustrate this. The true value, which is clearly visible as a lower curve with periodic null values ​​in the region [6000, 10000]. Quantity It appears as transparent red and represents "noise," with an average density of... The area around the warehouse reached its peak. The expected value is shown by a black dashed line. This is obtained as follows: using (75), the known Values, and assume they have known values. To obtain Gaussian noise. Quantity It is shown as a transparent blue curve. This is almost invisible because it intersects with the curve at intervals [0, 4000], [12000, 16000]. Closely overlap, and in intervals [5000, 11000] with They overlap closely. Note that... The color appears to change from red to purple within an interval of [5000, 11000] because the two transparent drawings overlap. Solid black lines indicate... It is A low-pass filter version. As mentioned in the paragraph about smoothing at the beginning of this section, the low-pass filter removes the distortion caused by peak position. Any local oscillations within. (Item) Used as The estimate. It can be seen that, in In the region, Provided Reasonable and good estimates. The quality of the estimate deteriorates as the two quantities approach each other, until it eventually becomes dominated by noise. Filter Should include To obtain a good estimate, while excluding poor estimates. In order to find a good estimate... To obtain a well-estimated region, we solved the following problem: given In local areas The calculated value is mainly based on the probability of noise occurring.

[0246] 7. Count Rate Estimation

[0247] Previous estimator assumptions It is known. Without prior knowledge, the following can be obtained regarding... The estimate.

[0248] 1. Use from the previous section ,calculate

[0249] .

[0250] 2. Use Estimate the count rate. This can be done in several ways.

[0251] 3. One approach is to use optimized routines or other methods to fit the curve to... The fitted parameters can be used to obtain an estimate of the count rate.

[0252] 4. Another approach involves estimation. The DC offset. This can be achieved by adjusting the appropriate number of... This is accomplished by averaging the points. This is done using the method described in the previous section. The points obtained by filtering are usually appropriate, although fewer points can also produce a sufficient estimate.

[0253] 5. Another approach involves using an optimization engine or other methods to fit the curve to... Fitting The appropriate parameterized curve is given by the following equation.

[0254]

[0255] in

[0256] And among them Selected to allow the curve to fit sufficiently accurately. Parameters It provides an estimate of the count rate. The optimization engine does not need to... Each point in the matrix is ​​given an equal weight.

[0257] Description of each of the 8 figures

[0258] The following diagrams will help to understand the process. Figure 1 A possible scheme for partitioning the detector output is shown. The illustration depicts the response of the sampling detector to three incident photons. Noise has been removed to make the diagram clearer. The output response is partitioned into several regions of equal length (L). The number of pulses arriving in each region is unknown to the processing system. One pulse arrives in the first interval. Two pulses arrive in the second interval. No pulses arrive in the third interval. The total photon energy arriving in each interval is calculated as the statistic of interest, which is the sum of all sample values ​​within each interval. The intervals are not time-aligned with the pulse arrival times. Figure 2 The output of the estimation process is illustrated. The true probability density of the incident photon energy is plotted as a solid black line. The photon arrival rate is such that, at any given interval... On average, three photons arrive within the detector. The detector outputs a signal. The standard deviation of additive noise is equal to the width of a histogram bin. One million intervals were collected. A histogram of the total energy in each interval was plotted. This is shown in blue. The effect of packing is clear and obvious, especially around bins 75, 150, and 225. The red trajectory plots the estimate of the true incident energy spectrum after the system processed the data. Although some noise appeared in the estimation, the effect of packing has been removed. It is expected that the estimate will correctly recover the true incident spectrum on average. The results were obtained using an internal filter whose bandwidth was automatically determined based on the data. Figure 3 The diagram illustrates the relationship between the two systems under the same operating conditions. Figure 2 The same number of parameters are used; however, in this case, the bandwidth of the internal filter is determined using asymptotic results from the literature. Although the estimated incident energy probability density has been recovered, it is different from... Figure 2 In comparison, the variance is significantly larger. Figure 8 The diagram illustrates the system's operation on real data. The blue trajectory plots the probability density of the observed energy values, while the red trajectory plots the estimated true probability density of the incident photon energy. There is no black trajectory because the true probability density is unknown. In this experiment, X-ray fluorescence from a manganese sample was used as the photon source. The photon arrival rate was approximately 10-1 / s. 5 5.9 photons. The interval length was chosen such that the average time between photons corresponds to the length of two intervals. Sufficient data was collected and partitioned to form 5.9. x l0 6 The intervals are [specified]. The standard deviation of additive noise corresponds to 4.7 histogram bins. The estimation process clearly reduces the accumulated peaks and enhances the true peaks. Figure 9The illustration shows various quantities obtained during simulation of the system described in the second embodiment. These are described in Section 5.1, “Visualization of Internal Quantities.” Figure 9-12 The second embodiment is involved. Figure 10 The diagram illustrates the probability density of the observed and actual input photon energies in this experiment, from which the following is derived: Figure 9-13 The black trajectory depicts the true probability density. The red trajectory depicts the expected observed density when an average of three photons arrive within a given interval. The blue trajectory depicts the actual observed density. Packing up to orders of magnitude can be seen in the observed density. Figure 10 Includes several plots generated from typical spectral systems. The actual incident photon density (“ideal density”) is plotted with a solid black line. Observational histograms obtained from partitioned time-series data are shown in dark blue. Spectral distortion due to pulse stacking is evident. Figure 13 Various internal quantities were plotted using a logarithmic vertical axis. The dark blue curve sloping towards the center of the plot is... The amount of green that crosses the drawing horizontally is... The cyan curve sloping towards the center of the drawing is... . Figure 11 The curve is illustrated. The trajectory in the complex plane. Figure 12 The illustration is similar to Figure 9 The internal quantity is [not specified], but there are also some additional signals. The level of noise is mainly represented by the red trajectory and the corresponding black dashed line. The value of the characteristic function of the histogram noise. The transparent green plot, forming a "noise peak" at the center of the graph, is an estimate. This amount is in Figure 9 It is drawn in blue and is almost invisible because it is... Figure 12 Not shown It's obscured. The horizontal trajectory with an average value of -3 is... The drawing. The cyan trajectory is... The magnitude of the characteristic function of additive Gaussian noise, which starts from zero at zero bin and slopes to a minimum around 8000 bins. Figure 13 This involves a second, more generalized aspect of calculating the count rate. It illustrates the calculation... The internal quantity used at that time. The cyan trajectory is... The magnitude of the characteristic function of additive Gaussian noise starts from zero at cell zero and slopes to a minimum around cell zero. The dark blue trajectory sloping to the minimum at the center of the graph is... It is an estimate of the characteristic function of the observed data. The yellow / green horizontal trajectory with a mean of -3 is... The estimate.

[0259] References

[0260]

[0261]

[0262]

[0263]

Claims

1. A method for determining the energy spectrum of a single radiation quantum received in a radiation detector, the method comprising the steps of: (1) Obtaining a time series of digital observations from a radiation detector, the time series comprising pulses corresponding to a single quantum detection; (2) Calculate the energy spectrum sensitivity statistics based on the detector signal. The energy spectrum sensitivity statistics use an approximate composite Poisson process to define the mapping from pulse amplitude density to energy spectrum sensitivity statistics. and (3) The density of pulse amplitude is estimated by applying the inverse of the mapping to the energy spectrum sensitive statistics, thereby determining the energy spectrum.

2. The method of claim 1, further comprising enabling spectral sensitivity statistics based on the sum of digital observations over multiple time intervals.

3. The method of claim 2, further comprising enhancing the approximate composite Poisson process by modeling noise.

4. The method of claim 3, further comprising expressing the mapping as a relationship between a characteristic function of amplitude, energy spectrum sensitivity statistics, and modeling noise.

5. The method of claim 4, further comprising calculating the characteristic function of the energy spectrum sensitivity statistics by applying an inverse Fourier transform to the histogram of the sum of digital observations.

6. The method of claim 4 or claim 5, further comprising calculating a characteristic function of the amplitude using a low-pass filter.

7. The method according to any one of claims 2 to 5, further comprising selecting each of the plurality of time intervals to contain zero or more approximately complete pulse clusters, and defining the plurality of time intervals as non-overlapping and having a constant length L.

8. The method of claim 7, further comprising requiring the detector signal to reach its maximum value at the beginning and end of each time interval.

9. The method of claim 7, further comprising defining the approximate composite Poisson process as the sum of the pulse amplitudes in each time interval.

10. The method according to any one of claims 2 to 5, further comprising selecting the plurality of intervals to include: The first set of non-overlapping time intervals with a constant length L is considered without regard to the overall pulse cluster. And a second set of non-overlapping time intervals of constant length L1 less than L, without considering the pulse cluster as a whole; where L is at least as long as the duration of the pulse.

11. The method of claim 10, further comprising selecting L1 as a duration less than the pulse duration.

12. The method according to any one of claims 1 to 5, further comprising using a data-driven strategy selected to result in a near-optimal choice of nuclear parameters that minimizes the integral square of the error of the estimated probability density function of the energy of a single quantum of radiation.

13. A method for estimating the count rate of a single radiation quantum received in a radiation detector, the method comprising the steps of: (1) Obtaining a time series of digital observations from a radiation detector, the time series comprising pulses corresponding to a single quantum detection; (2) Based on the sum of digital observations over multiple time intervals, calculate the energy spectrum sensitivity statistics from the detector signal, the energy spectrum sensitivity statistics being defined using an approximate composite Poisson process to map from pulse amplitude density to energy spectrum sensitivity statistics, the multiple time intervals including: a first set of non-overlapping time intervals of constant length L, which is selected without considering the overall pulse cluster; And a second set of non-overlapping time intervals of constant length L1 less than L, which are also chosen without considering the pulse cluster as a whole; where L is at least as long as the pulse duration. (3) Use the following formula to determine the estimate of the characteristic function of the approximate composite Poisson process. Where G is the windowing function. It is the characteristic function of the sum of numerical observations over each non-overlapping time interval in the first group. It is a characteristic function for modeling a noise process, and It is an estimate of the characteristic function of the sum of numerical observations over each non-overlapping time interval in the second group; (4) Estimate the count rate based on the estimation of the characteristic function.

14. The method of claim 13, further comprising estimating the count rate by: fitting a curve using an optimization routine or other means, estimating the DC offset of the estimated logarithm of the characteristic function, or fitting the curve to the estimated logarithm of the characteristic function.

Images