2198
Anal. Chem. 1988, 60.2198-2203
Statistical Evaluation of Dead Time Effects and Pulse Pileup in Fast Photon Counting. Introduction of the Sequential Model John A. Williamson* Department of Mathematics, University of Colorado, Boulder, Colorado 80309
Michael W. Kendall-Tobias,’ Marshall Buhl, and Michael Seibert Solar Energy Research Institute, Golden, Colorado 80401
Corrections to photocount dlstrlbutlons due to dead-the effects are evaluated by the use of a new pulse-plleup model lncorporattng a paralyzable portlon and nonparalyzable component In seqwnce. In the oknple model, the two dead thnes are regarded as nonvarylng. A more complete model is presented in the Appendlx, which Incorporates dead tlmes that are themselves random varlables.
Electron multipliers are widely used as flux transducers, whether incorporated into photomultiplier tubes (PMT) for measuring photon flux, or used “bare” for detecting ions, electrons, or electremagnetic radiation in the vacuum ultraviolet. The flux can be determined either by measuring the average direct current output of the electron multiplier (i.e., the “dc” method) or by counting individual charge pulses. The latter technique is often referred to as photon counting, and it is increasingly being used for analytical spectrophotometric measurements. There are several inherent advantages of this technique over the conventional dc method. First, the results are obtained in digital form providing direct access to digital processing of discrete spectral information. In addition, photon counting provides the ability to discriminate against the electron multiplier dark current, to improve the signal to noise ratio, to increase sensitivity at very low input levels, and to decrease effects related to drift and noise resulting from changes in voltage and temperature (I). Despite these attractive features there are some drawbacks to photon counting (2). By far the most serious is the phenomenon of pulse pileup or counting losses noticeable a t moderate to high flux levels. The applicability of this technique could be extended considerably if count losses due to dead time effects could be satisfactorily modeled and corrected. The literature deals with two main kinds of affected counting mechanisms. In the first, a pulse is registered if, and only if, no event has occurred during the preceding interval of time of length equal to the dead time after a single counting event. Each pulse extends the dead time whether or not the pulse is counted. Thus the dead time can be extended indefinitely and the mechanism involved is defined as being “paralyzable” (Figure la, top). The second type is termed “nonparalyzable”. In this mechanism, after each registration, the counter is locked for a constant dead time (Figure Ib, top). An event is registered only if the dead time has elapsed since the last registration. There have been many attempts to correct for pulse pileup using statistical methods and these usually have involved assuming one or the other of the above-mentioned counting mechanisms (3-6), although compromise models have been suggested (7, 8). Current address: Oriel Corp., 250 Longbeach Blvd., Stratford,
CT 06497-0827.
0003-2700/88/0360-2196$01.50/0
We propose here that a real counting system consists of both the paralyzable and nonparalyzable types of counters in sequence. We have analyzed the circuitry of a photon counting system (utilizing an EG&G Princeton Applied Research Model 1121 amplifier discriminator and a Model 1109 photon counter) and suggest that the counting system can be considered to be composed of four simplified parts: a detector (including the diode string), an amplifier discriminator, a pulse shaper (within the amplifier discriminator), and a counter. We propose that the detector and amplifier discriminator are both paralyzable systems and can be represented by a single rate-limiting paralyzable counting mechanism. The pulse shaper and counter have electronic components that are nonparalyzable and slower than the previous circuits. They can be represented by a single, rate-limiting nonparalyzable mechanism. The pulse-pileup effects can then be derived for the two sequential systems, each with its own dead time (see also ref 6).
BACKGROUND AND THEORY Theory of Counting Losses. A number of groups (5-7) have summarized the theory of counting losses by defining two sequences of events, Ri and R,. Ri is the mean number, or expected number, of photons in the primary sequence arriving at a detector per unit of time. R, is the number of events or counts per unit of time in the secondary sequence as recorded by the detectar. We will assume a quantum efficiency of 1 for the purposes of the discussion. The photon flux is regarded as random, and the probability of observing j events in a given time interval, t, is given by the Poisson distribution (6, 9, 10)
R, is only observed after the Poisson distribution with intensity Ri is altered upon passing through both the paralyzable and nonparalyzable systems. We make use of two random variables in obtained our estimation of Ri. The first, Y, is the period of time that the paralyzable system is sensitive to an arriving photon, i.e., the time between which the system turns on and will record an event and the time that it turns back off again. It is seen by setting j = 0 in eq 1that Y has an exponential distribution with density Rie7&,with y greater than or equal to 0. In other words
P(Y I t) = ltRie-YRidy 0
A second random variable, X,is defined as the time that the paralyzable system remains insensitive between successive periods of sensitivity. From its definition we see that X must be equal to or greater than T‘, a constant defined as the dead 0 1988 American Chemical Society
ANALYTICAL CHEMISTRY, VOL. 60, NO. 20, OCTOBER 15, 1988
-2.0 O
-8.0 0
2199
h
2
4
6
8
10
‘og1oR1
Flgure 2. Plot of log W v s log R,. (l/Rl)@lT‘= W = 7”). 7’ = 4.42 ns, T” = 16.01 ns.
+ ((l/R,) -
7’
This follows from integrating both sides of eq 2 from 7’ to infinity and using the fact that an integration by parts gives 1/ T ‘
E(X)= l m 0 [ -l F ( x ) ] dx
B
Integrating eq 2 from f to infinity gives
J,-f(x) dx = J,-[l - F(x - ~ ’ ) ] R ~ e -dx ~ i “ (4)
I”
Flgure 1. Representatlon of a paralyzable (A) or nonparalyzable (B) countlng system. The stream of photons (top), shown as R , pulses, is observed as a series of R , pulses after passing through a T’ (A) or a T” (B) filter. The plots (bottom) of R , vs R , are for a paralyzable (A) and nonparalyzable (E) system.
time of the paralyzable system. We will write F(x) = P(X 5 x ) and denote the density of X on the interval ( d ,+a) by f ( x ) . It is possible to write down expressions for F(x) and f ( x ) , but they are complicated and, at this point, not needed in our investigation. Instead we establish and use the following result. For x greater than 7’, we have f ( x ) = [I - F(x - 7’)]RieVRif
(2)
This identity could be established by using the explicit expressions for F ( x ) and f ( x ) , but a more direct argument is available for use. Using the notation of differentials, let dx be the length of a very short interval on the line x (i.e., dx is positive but small). Neglecting higher order terms we can write P(x < X < x dx) = f ( x ) dx. But for X to be between x and x + dx,the following three conditions must be met: (i) X must have been greater than x - 7’;(ii) there must have been a photon arrival between x - 7’ and x - T‘ + dx;and (iii) there must have been no photon arrivals between x - T’ + dx and x . The probabilities for these three events are 1 - F(x - #), Ridx, and respectively. These three events are independent and hence eq 2 follows. We can now use eq 2 to find E ( X ) ,the desired theoretical average of X
+
E ( X ) = (eRir’- l ) / R i
(3)
Setting x ‘= x - d on the right side of eq 4 shows that the right side equals E(X)Rie-Rif.The left-hand side of eq 4 must be examined more closely. The distribution function of X is made up of two parts. There is a positive probability that X takes on the value 7’. It is the probability that following the arrival of the initial photon, no photon arrives in the next interval of time of length 7’. This probability is P ( X = 7’) = e-Rif. The remaining probability is distributed over the interval x > T’ and has the density, f ( x ) . Therefore the value of the integral on the left side of eq 4 is not 1 but 1 - e-Rif. Replacing the integrals in eq 4 by their values and solving for E ( X ) then gives eq 3. We now estimate Ri in the following way. At t = 0 a photon arrives a t the detector causing both the paralyzable and nonparalyzablesystems to become insensitive. Let T%be the random time at which the nonparalyzable system is rendered insensitive for the (no 1)th time so that n,/Tno= R,. If we denote the dead time for the nonparalyzable system by f ’ , then Tn0- no#’ is the sum of 2n, independent random variables, and no of them have the same distribution of Y . The other no of these random variables have the same distribution as X - 7’. To see that X - 7’ is the correct expression, first notice that for our detector system 7’’ is greater than T ’ , so that when the nonparalyzable system becomes sensitive, the paralyzable system will be sensitive with probability e-&f,that is, P(X - f = 0 ) = e-Ri“. If when the nonparalyzable system becomes sensitive, the paralyzable system is insensitive due to the arrival of a photon during the previous 7’ time unit, then the additional time that the system remains insensitive will have the same distribution as X - f . The expected value of each of the no random variables having the same distribution as Y is l / R i . The expected value of each of the no random variables having the same distribution as X - 7’ is found by subtracting f from the right side of eq 3. Therefore
+
Adding T‘ to (Trio - norff)/no provides us with an estimate of (l/Ri)$ff. This approach actually yields two estimates of Ri because the function g(RJ = (l/Ri)eRiffor positive Ri is concave upward with a minimum at Ri = L/T’ (Figure 2 ) . Physical considerations, however, normally lead us to take as our estimate of Ri the smallest root of the equation 7’ (Trio - noTf’)/no= (l/Ri)eRiT’
+
2200
ANALYTICAL CHEMISTRY, VOL. 60,
NO. 20,
OCTOBER 15, 1988 6.0 I
By rearranging and noting that no/Tn0= R, we have
Equation 6 is the basic formula for the estimation of Ri in a sequential system. Note again that this equation assumes T”(nonpara1yzable)> T’(para1yzable) are constant. If 7” I T’, then eq 6 reduces to R, = Rie-Rif, This is in agreement with the result for a single paralyzable system (see Figure la, bottom). Letting T’ approach zero in eq 6 gives Ri = R,/(l - R,T”), which is in agreement with the result for a single nonparalyzable system (see Figure lb, bottom). Variance in the Estimation of Ri.We turn to the question of the amount of error in our estimate of Ri. The first identity that we require is
J m x ( l - F ( x ) ) dx = E ( X 2 ) / 2
(7)
This follows from an integration by parts. Multiplying both sides of eq 2 by x , integrating from 7’ to infinity, and applying eq 7 yields the formula for the second moment of x
E(X ’) = 2E(X )(eRiT’ - 7’Ri)/Ri - 2(7’/Ri)
(8)
More specifically eq 8 follows from writing
S,,-(x - #)(I- F(x - 7’)) dx
+
L m 7 ’ ( 1-F(x - 7’)) dx = E ( X 2 ) / 2
+ T’E(X)
and from writing
x , m x f ( x )dx = E(X) - 7’e-Rir’ If from both sides of eq 8 we now subtract (E(X))’, we obtain Var(X - 7’) = Var(X) = [ ( e 2 R i f - 1)/Ri2] - [27’eRi“/Ri] (9) The variance of Y is Var(Y) = l/RP. If we let W = T’ + (Tm0 - r”no)/no,then our estimate, W , of (l/Ri)& is seen to have a variance given by Var(W) = [(eZRir‘/Rt)- (27’eRir’/Ri)]/n, (10) The standard error for W is given by [(W- 27’W)/n0]’’2 (11) If we restrict the domain of the function g(RJ = (l/Ri)eRif to either Ri values between 0 and 1 / ~or’ to Ri values between 1 / ~and ’ +a, then we can write our estimate of Ri as g-’(W). A measure of the possible error in our estimate of Ri is then given by E-’( - ( - 2T’w)’i2/n0’/2) - g-’( w)l (12) Taking square roots in eq 10 gives the standard deviation (SD) of w
w w2
Dependence of the Standard Error on Counting Time and Rate. Equation 13 gives us an estimate of how much error there may be in our calculation of the incident photon flux. Before any practical application of the sequential model is attempted, it is important to quantify the magnitude of the variance over a range of count rates and total counting times, since this will (theoretically) determine the limits of the model’s usefulness. We can define the perceived signal to noise ratio as follows: incident photon count (ni) S/N = standard error in ni(W)
0
2
6
4
8
10
Log R ,
Flgure 3. Calculated plots of the log of the signal/noise ratio vs log of the input count rate R , for count times of 0.1, 1, 10, 100, and 1000 s (in ascending order). T‘ = 5.79 ns, 7“ = 12.15 ns.
Table I. Theoretical Relationships”
incident std error count ni obsd count no in ni (W) 10’ 102 103 104 105 106 107 108
10 100
4.62
“Count time = 1 s.
3.16 10.0 31.62
11.1
1000
9.998 X 9.984 X 9.842 X 8.612 X 3.683 X
S I N ratio
lo3 lo4 lo5 lo6 lo7 7’ =
32.7 101 317 1009 3398 19416 4.42 ns.
7”
100
316.27 1000 3162.3 10000 =
2.16 9 31 99 315 991 2943 5150
16.01 ns.
where the expression “standard error in ni(W ) ”is to be understood to mean the standard error in ni calculated by using eq 12. We can now plot the logarithm of the SIN against the logarithm of the incident count rate (Ri) for differing count times (see Figure 3). It can be seen that the log S I N ratio rises linearly with log Ri until it dips suddenly. This dip is the point at which the maximum observable count rate occurs and at which Ri = l / ~ ‘It. is also the point at which g(Ri) is at a minimum and is equal to e# (see Figure 2). It is evident from Figure 2 that the estimation of the variance will be very large at this point. After Rj = l / ~ ’log , ( S / N )rises rapidly to a new maximum and falls back again toward zero as pulse pileup becomes increasingly large. Table I shows the theoretical relationahip between ni for a 1-s counting time, the observed count no for a 1-s counting time, the standard error in ni(W),ni1I2,and the S I N ratio. The is the standard error when there is no pulse pileup, and it can be seen when the incident photon flux is kept below approximately 1/2r’ that the true error closely follows the ni1j2. Monte Carlo Simulation of the Sequential Model. The applicability of the sequential model for pulse-pileup correction in a real system is primarily dependent on the algorithm’s accuracy (eq 6) for a “perfectly”Poisson distributed light flux and a “perfect“ detector with known dead times. This arrangement was simulated on a CDC Cyber 70-76 computer. A series of random time intervals, distributed according to an exponential law, was generated to simulate the times between photons arriving at the detector. These “photons” were passed through a sequential system with two known nonvarying dead times, T’(2 ns) and 7’‘ (10 ns). We used a photon flux of (1X 106)/sand ”stopped”the simulated experiment after the detection of 50000 photons. The “experiment” was repeated a hundred times. The estimated input count rate determined from eq 6 was found to have a mean error of f0.03670 from the actual input count rate. Of the one hundred trials, 67 of the estimated rates were found to lie within the standard deviation (eq 13), 17 estimated were above, and 16 were below these bounds. The error in detection, had no correction been applied, would have been about
ANALYTICAL CHEMISTRY, VOL. 60, NO. 20, OCTOBER 15, 1988
-1%. Thus eq 6 can accurately be applied within the constraints of our theoretical model. EXPERIMENTAL SECTION Count Correction in Real Counting Systems. Some Practical Considerations. Considering the complexity of applying the sequential model for pulse-pileup correction of real data, we have limited our experimental application to the simple sequentialmodel described above. This necessarily involves many assumptions that are approximationsbut that may be acceptable due to their small effect on the accuracy and the usefulness of the simple model. Specifically, we have excluded considerations of the following realities: 1. As described previously, we are regarding all detectable pulses as being of nonvarying shape. The width of a pulse as recorded by the detector after passing through a paralyzable electronic circuit is the dead time, T’. 2. We are assuming all detectable pulses are of the same height and that the amplifer discriminatorthreshold voltage is irrelevant. In practice, we have set the threshold at an arbitrary value that would eliminate most background pulses but would let through most pulses due to photon arrivals. 3. Background counts have been ignored; since these generally are less than 0.01% of the recorded count, we feel it is justifiable. It is worth noting that at low light levels where background is an important consideration, application of the sequential model is not necessary. Whereas at the higher count rates where the model becomes useful, the background is virtually nonexistent. 4. We have not attempted to measure the dead times T’ and f’in our experiment, since this is beyond our capabilities. Thus we have no absolute proof of the applicability of our model. We have, instead, derived T’ and 7’‘ from calibration experiments assuming the sequential model is a good approximation to reality. Although this method has its obvious criticisms, we have relied on the closeness of fit of the model and on the believability of the derived dead times to press our case. Methods. The photon counting equipment consisted of a Hamamatsu R955P photomultiplier tube and a Princeton Applied Research 1109 photon counter with a Princeton Applied Research 1121A amplifier discriminator (EG&G, Princeton, NJ). The photomultiplier tube was cooled to -15 OC in order to reduce thermal counts to below 100 counts/s. It was operated at 1200 V, uniformly distributed between the cathode, anode, and the nine diodes. The power supply was a Kepco APH 2000M (Kepco, Flushing, NY). The photon counter was used in its fast counting mode, and the amplifier discriminator threshold voltage was usually set at 1.5 mV. The light source was a 250-W quartz halogen tungsten lamp rendered monochromatic (1 nm bandwidth) by passing the light through an ISA 430 monochromator (Instruments S.A., Metuchen, NJ). A Keithley 427 current to voltage amplifier (Keithley Instruments, Cleveland,OH) was used for all dc measurements of the photocurrent. Where possible, all operations were directly controlled by a MNC PDP 11/23 computer (Digital Equipment Corp.) using a program called PCCALB written for this purpose. However cable connections for reading the photocurrent or photocounts were separate and had to be changed manually in order to eliminate interferencesbetween the two measuring systems. The light flux either was computer controlled by varying the current supply to the light lamp from a Kepco ATE 25-10M power supply via a 12-bitdigital-tc-analog interface (DEC MNCAA) or was controlled manually by opening and closing the monochromator slits. In all cases, the data presented were taken under conditions most likely to resemble operating circumstances. No adaptation period was allowed between readings and runs, other than to allow for conditions to stabilize. RESULTS AND DISCUSSION Calculation of T’ and 7’’from Experimental Data. As mentioned in the Experimental Section, the dead times must be derived empirically from a set of data, obtained for the specific purpose of calibrating the instrumentation. The relevant data are the observed photocount rate and the observed photocurrent at different light fluxes. If the exact dead times, f and f’,were known for eq 6, then a linear relationship
2201
Table 11. Results of a Three-Parameter Fit for Calibration Data from Three Separate Experimentsa experiment parameter
1
slope (counts&pA-’) T’, ns estd uncertainty, ns T”, ns estd uncertainty, ns % std dev in data % corr of nearest neighbor devns
3
2
1300778 1311203 1314808 4.33 4.44 4.51 0.11 0.09 0.10 16.07 16.01 15.94 0.07 0.07 0.07 0.24 0.22 0.23 -56.9 -33.6 -33.0
“The amplifier discriminator threshold was set at 1.5 mV, and the photon counting time was 10 s. ? 1’5
8
a
r
0.5
’‘oLf&?Ro 00
20
40
60
80
100
PMT Current b A )
Figure 4. Plot of corrected photocounts vs photocurrent for the paralyzable (P), sequential (S), and nonparalyzable (NP) models shown along with a plot of the observed photocounts R,. is to be expected between Ri and the photocurrent, having a constant slope, s, with units of counts/kA. We can find the values of the three unknown parameters f,d’, and s by using an iterative process and deducing the values that give the best linear least-squares fit for Ri,or equivalently, by deducing the values that give the best least-squares fit of the curve R, = (T” - T’ e m S ~ / m s )to - l the R, data. Here the variable, m, is microamperes (11). A computer program, TAUFIT, that performs this fit can be supplied by the authors on request. Table I1 shows the results for three sets of calibration data taken under the same instrumental conditions. There is no significant drift or variation in the values of T’ or T”. Overall the standard deviation of Ri from the linear function is fairly constant, a t about 0.23%. However, this figure can be misleading since there is no indication of whether the variation is evenly distributed throughout the count range or is more pronounced at low or high count rates. In fact the variation is most significant at low count rates and seems to be due to the difficulty that was experienced in obtaining stable readings at low light levels. The instability was predominantly in the light source. Accuracy i n t h e Applied Model. Let us now consider accuracy in the count correction. A plot of photocounts versus photocurrent is shown in Figure 4 (R,). The phenomenon of pulse pileup is evident and amounts to about 60% of the input counts at lo8 counts/s. These data were applied to a threeparameter fitting technique to solve for d and T” (4.5 and 15.9 ns, respectively, in this experiment). Furthermore, Figure 4 depicts a plot of Ri against photocurrent for three different pulse-pileup corrections. The plots are the best fits for the calculations of Ri using the paralyzable model (P), the sequential model (S),and the nonparalyzable model (NP). It can be seen that only the sequential model developed in this study gives values of Ri that are close to a linear function of the photocurrent. The value determined experimentally for T’, 4.5 ns, is close to the 3.6 ns fwhm pulse width of the photomultiplier tube, whereas the value of 9.65 ns for f from the simple paralyzable model is much larger (Table 111). The value for T” for our
+
2202
ANALYTICAL CHEMISTRY, VOL. 60, NO. 20, OCTOBER 15, 1988
Table 111. Comparison of the Dead Times for the Three-Pulse Pileup Modelsn 7‘,
approx specifications sequential model paralyzable model nonparalyzable model
ns
3.6 4.5 9.65
T”,
ns
10 15.9
17.3
“Note: This table is based on the data from Table 11, experiment 1. model, of 15.9 ns, is less than, but similar to, the value of 17.3 ns for the nonparalyzable model, both of which are higher than the nominal dead time of 10 ns for the output of the photon counter which is rated at lo8 counts/s. The major source of inaccuracy comes from the uncertainty in the values determined for 7’ and T”. From Table I1 we see that the estimated uncertainty in r‘ is 0.1 ns or 2.2%,whereas that in 7” is 0.07 ns or 0.43%. This implies that in fitting T’ and 7’’ to our data only a small amount of uncertainty in 7” can be tolerated compared to the uncertainty in T’, and, therefore, that f’ plays a more significant role in our sequential model than does T’. This is borne out in Figure 4. The nonparalyzable (7”) model is closer to our sequential model than is the paralyzable (7’) model.
CONCLUSION While the simplified sequential model that is presented here overlooks certain aspects of the pulse-countingprocess (pulse shape, threshold levels, etc.), it is sufficient fundamentally to describe the count losses in present day counting systems. We have shown that the derived dgorithm (eq 6 ) is applicable to real systems with repeatable results for the calculation of the dead times required to fit the algorithm to sets of calibration data (Table 11). The standard error in the estimation of Ri, based on 100 simulated readings of R, and on eq 12, was found to be very similar to the root mean square for the 100 estimates of Ri indicating that the error introduced by the count corrections technique was not significant. This is in agreement with the theoretical calculations found in Table I. Figure 4 shows that the sequential model gives a closer fit to a linear relationship between Ri and photocurrent than do either the paralyzable or nonparalyzable models alone. Furthermore the values derived for the dead times 7’ and T” are realistic values reinforcing the validity of the proposed models. The nonparalyzable portion of the model dominates, and the uncertainty in 7” has a more marked effect on the calibration than that of 7’. When considering modern electronic counting systems, it must be recognized that count losses cannot be described by simple models and that pulse pileup is the result of several mechanisms superimposed on one another. This paper shows that combining two well-known models for pileup into a sequential model can better portry the photon counting process than either of the simple models. More involved analysis, such as that presented in the Appendix, is likely to lead to even more accurate representations of the pulse counting process.
ACKNOWLEDGMENT The authors thank Dr. J. S. Connolly for his helpful discussions about this work.
APPENDIX We turn next to the more general problem where the values T’ and 7’’ are no longer fixed but are themselves random variables. Following the arrival of a photon and in the absence of any later arriving photons, the paralyzable system will be dead for a random length of time. This random length of time
we will denote by T! After the nonparalyzable system records an event, it will be dead for a random length of time. This second random time we denote by T”. In these paragraphs we will use T’ and 7” to denote the expectations of these two random variables: E(T’) = T’ and E(T”) = 7”. As was done earlier, we denote by Y the length of a random time interval between dead times in the paralyzable system. Y then has an exponential distribution with E(Y) = l/Rb This is as before and is not altered by the fact that we are now allowing T’and T”to be random. The random variable 2 is the total time that the paralyzable system is dead beyond the time T”p1us the value of T”, An alternative way to describe 2 is as follows. A photon arrives and both systems record the event. The time until both systems are able to record another photon arrival is the random time 2. If T‘and T”are nonrandom with fixed values 7‘ and 7” and if X is as defined in the text, then 2 has the same distribution as the distribution of X + (7’’ - 7’). In order to formulate our main result, we need the following notation: H ( t ) = P(T’I t ) and H*(t) = P ( T ” It).
Result:
E ( Z ) = eR,r‘Sme-R,ia’(l-H(u))du(l 0 - H(t)H*(t))dt
(14)
Proof We define a Markov stochastic process, X ( t ) , as follows: X ( t ) = 0 if both systems are able to record a photon arrival a t time t. If at time t both systems are not able to record a photon arrival, then X ( t ) is equal to the amount of time that must pass, in the absence of any future photon arrivals, until the two systems can both record a photon arrival. A trajectory or sample path for X ( t ) might look like the following: For t < t,, X ( t ) = 0. At time tl a photon arrives. X ( t l ) is then equal to max(T’,T”) where T’and T”are the dead times in the two systems caused by the arrival of this photon. After t,, X ( t ) decreases linearly with slope -1 until the arrival of the next photon at time t2 If the dead time, T’, for this second photon is 5 X ( t 2 - ) ,then X ( t 2 )= X(t,-). If T’is greater than X(t,-), then X ( t , ) = T’. In either case, for t > t2, X ( t ) again decreases linearly with a slope of -1 until the time of the next photon arrival. A t the time of this third photon arrival X ( t ) either will jump to a random value T’ for this third photon or will continue to move continuously toward 0 with slope -1. The random time from t , to the time X ( t ) again touches 0 has the same distribution as 2. We denote the stationary distribution for this Markov stochastic process by G(x). This stationary distribution will have an atom of probability at 0 given by G(0) and a density over positive x values that we will denote by g(x). Denoting the densities of H ( t ) and H*(t) by h(t)and h*(t),and letting dt be the length of a small positive increment of time, permits us to write for positive x g ( x ) = Ri dt[G(x) - G(O)]h(x) + g(x + dt)[l - Ri d t ] + Ri d t G(O)(H(x)H*(x))’ + g(x + dt)Rj d t H ( x ) + o(dt) (15) Moving g(x + dt) to the left side in eq 15, dividing both sides by dt, and letting dt approach 0 then gives -g’(x) =
RiG(O)[(H(x)H*(x))’- h(x)l
+ Ri[G(x)(H(x) - 111’ (16)
Integrating both sides gives -&)
=
RiG(O)[H(x)H*(x) - H ( x ) ] + RiG(x)[H(x) - 11 + C
(17) In eq 17 letting x increase without bound shows that C = 0. Writing g ( x ) = G’(x) and multiplying both sides of eq 17 by
ANALYTICAL CHEMISTRY, VOL. 60, NO. 20, OCTOBER 15, 1988
the integrating factor eRiJ-ot(H(u)-l) du then gives
G(t)eRiJO’(H(U)-l) dU =
+
-RiG(0)S’(H(s)H*(s)- H(s))eRiJo’(H(u)-l)du ds C (18) 0
Letting t approach 0 from above in eq 18 gives G(0) = C. Letting t approach infinity in eq 18 gives
G(0) = @I#[
1 - RiJ
+m
(H(s)H*(s)- H(S))e-RiJo‘(l-H(U)) du ds
1-’
(19) But G(0) is the limiting probability of the event that both systems are able to record a photon arrival and therefore
G(O) = (1/Ri)/[E(Z)
+ (1/Ri)l
(20) Equating the right sides of (19) and (20),writing H(s)H*(s) - H ( s ) = H(s)H*(s) - 1 + 1 - H ( s ) and integrating where possible, and then solving for E ( 2 ) establishes (14). We use eq 14 to estimate Ri in the following way. As before we let Tnobe the time from the beginning of the first dead time interval to the beginning of the (no + 1)st dead time interval (or from the time the investigator first begins observing to the end of the n,th dead time. With this second definition, the last exponential wait is replaced by a first exponential wait, the waiting time until the first photon arrives). Tn0is then the sum of 2n, independent random variables, no of them having the same distribution as Y and the other no having the same distribution as 2. Hence E(T,Jn,,) = E(Z) E ( Y ) = E ( 2 ) + l / R i (21)
+
The right-hand side of (21) is a function of Ri and we denote it byg*(RJ. This is not the same function as the functiong(Ri) defined in eq 4-6 but we use the same notation because we use the function in the same way in order to obtain our estimate of Ri. If Tno/no can be expected to be near g*(Ri), we take as our estimate of Ri
fii = g*-’(T,,/n,)
(22) The graph of our newly defined g*(Ri) is similar to the graph of eq 11in that it is concave upward and approaches infinity as Ri nears 0 and as Ri grows large. There will be, then, in general two estimates of Ri corresponding to a given Tno/no and, as before, other considerations must be introduced in order to decide which estimate is the meaningful one. If we let the variances of T’ and T” approach zero, eq 14 should take on a familiar form. To see that this is indeed the case, we first let T’ = T” and let the variances tend to zero. Then for t greater than T’, H ( t ) = H*(t) = 1, and for t < T ’ , H ( t ) = H*(t) = 0. In this case eq 14 becomes
E ( 2 ) = eRifSfe-RItd t = (eRif- 1)/Ri 0
(23)
But the right-hand side of eq 23 is E ( X ) as expected. If we now let the variances approach zero but with T’ < r” we get
The right-hand side of eq 24 is E ( X ) + T” - T’. If we remember that 2 is the total amount of dead time while X is only the time that the paralyzable system is dead following the recording of a photon arrival, we see also that eq 24 is as expected. One would like to have an estimate of the error involved in our estimate of Ri. To obtain such an estimate we would
2203
have to have an expression for the variance of 2. Unfortunatly, the method developed to derive eq 14 cannot be used to obtain the variance of 2. Without a formula for the variance of 2, it is necessary to turn to a simulation for error estimates. It should be pointed out that there are a number of formulas in the literature that claim to give information on the amount of error that can be expected when T’ and T“ are random variables. These results deal with systems that are either paralyzable or nonparalyzable but not with systems that are a combination of both kinds of instrumentation. However, that is not our objection to these formulas. The formulas that we are referring to are formulas that begin by writing
Ri = R o / ( l - ROT’’)(nonparalyzable case)
(25)
R, = Rie-Rir (paralyzable case)
(26)
Using these equations, the authors then relate the variation in T’ (or 2’”) and the variation in the estimatesof Ri. By doing this, these authors are using an incorrect model. They are assuming that a value for T’ (or T”j is chosen according to the distribution governing T’ (or T’? and then this value is used in the system for each arriving photon. In this incorrect model, there is a variability in T’ (or T ” ) and this variability does introduce variation in the estimate of Ri. It is the variation that arises in the initial selection of T r(or T ” ) . Once the initial selection is made, there is no more variability in T’ (or T ” ) in this incorrect model. The correct model, the one that we are working with, takes into account the fact that each photon arrival produces a possibly different T’ (or T ” ) value. To get a t this more complicated, but more realistic, variation, either one needs an analytic expression for Var(Z) or one needs to do simulations. It is true that in our more complicated model where T’ (or T ” ) vary from photon to photon, if there is only a nonparalyzable system or if there is only a paralyzable system, then our formulas for the estimate of Ri reduce to eq 25 and eq 26, respectively, but with T’ and T“ replaced by T’ and 7“. In either the sequential model, or paralyzable model, or nonparalyzablemodel, because T’or T”will vary from photon to photon, the variability in our estimate of Ri will be greater than that given in eq 10 and hence the perceived signal-tonoise ratio will be reduced from the value given by the simpler model.
LITERATURE CITED (1) Franklln, M. L.; Horllck. G.; Malmstadt, H. V. Anal. Chem. 1968, 4 1 , 2-10. (2) Ingle, J. D.; Crouch, S. R. Anal. Chem. 1972, 44, 785-794. (3) Srinlvasan, S. K. J. Phys. A 1978, 1 1 , 2333-2340. (4) Hayes, J. M.; Matthews, D. E.; Schoeller, D. A. Anal. Chem. 1978, 50, 25-32. (5) Hayes, J. M.; Schoeller, D. A. Anal. Chem. 1977, 4 9 , 306-311. (6) Ingle, J. D.; Crouch, S. R. Anal. Chem. 1972. 44, 777-784. (7) Feller, W. Studies and Essays; Interscience: New York, 1948; pp 105- 115. (8) Albert, G. E.; Nelson, L. Ann. Math. Stet. 1953, 24, 9-22. (9) Mandel, L. Progress in Optics; North-Holland: Amsterdam, 1963; Vol. 2, Chapter 5. (10) Berdard, G. Roc. Phys. SOC. London 1967, 90, 131-141. (1 1) Gorman, D. S.; Connolly. J. S. International Journal of Chemlcal Kinetics; 1973; Vol. V, pp 977-989.
RECEIVED for review November 24,1987. Resubmitted June 20,1988. Accepted June 20,1988. This work was sponsored in part by the Division of Energy Biosciences, Office of Basic Energy Sciences, U.S. Department of Energy (MS). The Solar Energy Research Institute is a Division of the Midwest Research Institute and is operated for the United States Department of Energy under Contract No. DE-ACO283CH10093.
