Coverage for larch/xrf/deadtime.py: 14%
86 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1"""
2Deadtime calculations
4Authors/Modifications:
5----------------------
6* T. Trainor (tptrainor@alaska.edu), 6-10-2006
8Description:
9------------
10The objective is to calculate a factor 'cor' that will provide
11a deadtime correction factor, and give corrected counts via:
13 counts_corrected = counts * cor
15Background:
16------------
17A correction factor can be defined as:
19 cor = (icr_t/ocr_s)*(rt/lt)
21Here icr_t = true input count rate (actual counts hitting the detector
22per detector live time - ie icr is a rate)
24ocr_s is the output rate from detector slow channel: ocr_s = TOC_s/lt
25TOC_s is the total counts output from the slow filter (~ the total counts output
26by the detector).
28rt and lt are the detector real time and live time respectively. real time is
29the requested counting time (=time elapsed from detector start to detector stop).
30live time is the amount of time the detector is active. e.g. lt/rt*100% is the
31percent of the counting time that the detector is actually live.
33icr_t is an uknown quantity. It can be determined or approximated in a few
34different ways.
36A) icr_t may be determined by inverting the following:
38 ocr_f = icr_t * exp( -icr_t * t_f)
40Here ocr_f is the fast count rate. ocr_f = TOC_f/lt, TOC_f are the total counts
41output from the fast filter.
43t_f is the fast filter deadtime. If the detector reports TOC_f and t_f is known
44(from fit to a deatime curve) then icr_t may be calculated from the above.
46Note: a detector/file may report InputCounts or ICR.
47This is likely = ocr_f rather than icr_t.
49Note: The fast filter deatime is much less substantial than the slow filter
50deadtime. Therefore for low total count rates the approximation icr_t ~ ocr_f
51may be reasonable.
53B) alternatively icr_t may be determined by inverting the following:
55 ocr_s = icr_t * exp( -icr_t * t_s)
57Here ocr_s is the slow count rate. ocr_s =TOC_s/lt, TOC_s are the total counts
58output from the slow filter.
60t_s is the slow filter deadtime. If the detector reports TOC_s and t_s is known
61(from fit to a deatime curve) then icr_t may be calculated from the above.
62If the detector does not report TOC_s ocr_s may be approximated from:
64 ocr_s = ( Sum_i ( cnts_i) ) / lt ~ TOC_s/lt
66Note the above discussion referes primarily to digital couting electronics,
67using slow and fast filters. Analog electronics using an SCA may be corrected
68using a similiar procedure. Here the analog amplifier may
69report total input counts. In this case icr_t ~ icr_a (icr amplifier).
72Measuring t's:
73--------------
75To perform a deadtime correction the characteristic deadtimes must be determined.
76In each case (t_f or t_s) a deadtime scan is run, where the input counts are
77varied (eg scanning a slit) and output counts monitored.
79If a measure of icr_t is reported (or approximated ie = ocr_f or icr_a), then
80t_f may be deterimined directly from a fit of
81 * x = icr_t
82 * y = ocr_s
84More likely a detector not suffering from deadtime will be used as a proxy for
85icr_t (eg mon -> ion chamber)
87In this case plot:
88 * x = mon/rt
89 * y = ocr_f -> TOC_f/lt reported in file (may be called icr or total counts)
90 * icr_t = a*mon/rt + off
91 * fit varying t_f and a and off (slope and offset btwn mon and icr_t)
92 ocr_f = icr_t*exp(-icr_t*t_f)
94If the detector does not report ocr_f, then use slow counts to correct
95(or just sum all the actual output counts):
96 * x = mon/rt
97 * y = ocr_s -> TOC_s/lt ~ ( Sum_i ( cnts_i) ) / lt
98 * icr_t = a*mon/rt + off
99 * fit varying t_s and a and off (slope and offset btwn mon and icr_t)
100 ocr_s = icr_t*exp(-icr_t*t_s)
102Summary of correction methods:
103------------------------------
104To apply a deadtime correction then depends on what data can be accessed from
105the detector.
1071) detector reports actual icr_t (or reports ocr_f and total counts are low so
108 ocr_f ~icr_t). Then the correction may be calculated directly from
109 detector numbers
1112) The detector reports ocr_f and t_f has been determined. icr_t is calculated
112 from the saturation equation to use in the correction.
1143) The detector does not report ocr_f, rather ocr_s is reported (or approximated
115 from total counts) and t_f has been determined. icr_t is again calculated
116 from the saturation equation and used for the correction.
117 This is probably the most straightforward method since icr_s can be
118 approximated directly from the sum of the counts in the detector (norm by lt),
119 and this approach should work for analog and digital electronics.
1214) icr_t is uknown or icr_t ~ ocr_s then just assume that icr_t/ocr_s = 1
122 (ie just correct for livetime effects)
124Example:
125--------
126# given the arrays mon (ion chamber), rt and lt (each len = npts) and
127# the (numpy) multidimensional array for counts (e.g. nptsx2048)
128>>ocr = counts.sum(1)/lt
129>>x = mon/lt
130>>(params,msg) = fit(x,ocr)
131>>tau = params[0]
132>>a = params[1]
133>>print('a_fit= ',a,' tau_fit=', tau)
134# corrected counts
135>>counts_cor = num.ones(counts.shape)
136>>ocr_cor = num.ones(mon.shape)
137>>for j in range(len(counts)):
138...>icr = calc_icr(ocr[j],tau)
139...>cor = correction_factor(rt[j],lt[j],icr[j],ocr[j])
140...>counts_cor[j] = counts[j]*cor
141...>ocr_cor[j] = counts_cor[j].sum()/lt
142>>pyplot.plot(x,ocr)
143>>pyplot.plot(x,ocr_cor)
144"""
146##############################################################################
148import numpy as np
149import scipy
150from scipy.optimize import leastsq
151from scipy.stats import linregress
153E_INV = np.exp(-1)
155##############################################################################
156def correction_factor(rt, lt, icr=None, ocr=None):
157 """
158 Calculate the deadtime correction factor.
160 Parameters:
161 -----------
162 * rt = real time, time the detector was requested to count for
163 * lt = live time, actual time the detector was active and
164 processing counts
165 * icr = true input count rate (TOC_t/lt, where TOC_t = true total counts
166 impinging the detector)
167 * ocr = output count rate (TOC_s/lt, where TOC_s = total processed
168 {slow filter for dxp} counts output by the detector)
170 If icr and/or ocr are None then only lt correction is applied
172 Outputs:
173 -------
174 * cor = (icr/ocr)*(rt/lt)
175 the correction factor. the correction is applied as:
176 corrected_counts = counts * cor
177 """
178 if icr is not None and ocr is not None:
179 return (icr/ocr)*(rt/lt)
180 return (rt/lt)
182def correct_data(data,rt,lt,icr=None,ocr=None):
183 """
184 Apply deatime correction to data
185 """
186 cor = correction_factor(rt, lt, icr, ocr)
187 return data * cor
189def calc_icr(ocr, tau):
190 """
191 Calculate the true icr from a given ocr and corresponding deadtime factor
192 tau using a Newton-Raphson algorithm to solve the following expression.
194 ocr = icr * exp(-icr*tau)
196 Returns None if the loop cannot converge
198 Note below could be improved!
199 """
200 # error checks
201 if ocr is None or tau is None or ocr <= 0:
202 return None
204 # here assume if tau = 0, icr=ocr ie icr/ocr =1
205 if tau <=0:
206 return ocr
208 # max_icr is icr val at top of deadtime curve
209 # max_ocr is the corresponding ocr value
210 # we cannot correct the data if ocr > ocr_max
211 max_icr = 1/tau
212 max_ocr = max_icr*E_INV
213 if ocr > max_ocr:
214 print( 'ocr exceeds maximum correctible value of %g cps' % max_ocr)
215 return None
217 # Newton loop
218 icr0 = ocr
219 cnt = 0
220 while cnt < 100:
221 cnt += 1
222 delta = (ocr*np.exp(icr0*tau) - icr0) / (icr0*tau - 1)
223 if abs(delta) < 0.01:
224 icr = icr0
225 break
226 else:
227 icr0 = icr0 - delta
228 if icr0 > max_icr:
229 # went over the top, we assume that
230 # the icr is less than 1/tau
231 icr0 = 1.1 * ocr
233 if cnt >= 100:
234 print( 'Warning: icr calculation failed to converge')
235 icr = None
236 return icr
238##############################################################################
239def fit_deadtime(mon, ocr, offset=True):
240 """
241 Fit a deatime curve and return optimized value of tau
243 This fits the data to the following:
244 x = icr_t = a*mon + off -> linear relation btwn icr and mon
245 y = ocr -> TOC/lt
246 fit varying 'tau', and 'a' and 'off' (slope and offset btwn mon and icr_t)
247 ocr = icr_t*exp(-icr_t*tau)
249 This function is appropriate for fitting either slow or fast ocr's.
250 If mon is icr_t the the optimized 'a' should be ~1.
252 Parameters:
253 -----------
254 * mon is an array from a linear detector. This should be in counts/sec
255 (ie mon_cnts/scaler_count_time)
257 * ocr is an array corresponding to the output count rate (either slow or fast).
258 This should be in counts per second where
259 ocr = TOC/lt,
260 TOC is the total output counts and lt is the detector live time
262 * If offset = False, off = 0.0 and only tau and a are returned
264 Example:
265 --------
266 >>params = fit_deadtime(mon/time, ocr)
267 >>tau = params[0]
268 >>a = params[1]
269 >>off = params[2]
271 """
272 mon = np.array(mon,dtype=float)
273 ocr = np.array(ocr,dtype=float)
275 npts = len(mon)
276 if len(ocr) != npts:
277 return None
279 # make a guess at tau, assume max(ocr) is the top of the deadtime curve
280 tau = E_INV/max(ocr)
282 # make a guess at linear params, ie if the detector is linear
283 # ocr = icr = a*mon + off
284 # assume that mon and ocr are arranged in ascending order, and the first 10%
285 # have little deadtime effects (or maybe a bit so scale up the average 10-20%).
287 idx = max(int(0.10*npts), 3)
288 try:
289 xx = linregress(mon[0:idx],ocr[0:idx])
290 a = xx[0]
291 off = xx[1]
292 except:
293 mon_avg = mon[:idx].sum()/idx
294 ocr_avg = ocr[:idx].sum()/idx
295 a = 1.2*ocr_avg/mon_avg
296 off = 0.0
297 if offset:
298 params = (tau, a, off)
299 else:
300 params = (tau,a)
301 return leastsq(deadtime_residual, params, args = (mon, ocr, offset))
303def calc_ocr(params, mon, offset):
304 """ compute ocr from params"""
305 tau, a = params[0], params[1]
306 off = 0
307 if offset:
308 off = params[2]
309 return (a*mon+off) * np.exp(-tau*(a*mon+off))
311def deadtime_residual(params, mon, ocr, offset):
312 """ compute residual """
313 return ocr - calc_ocr(params, mon, offset)
315##############################################################################
316if __name__ == '__main__':
317 # test fit
318 mon = 10000. * np.arange(500.0)
319 a = 0.1
320 tau = 0.00001
321 print( 'a= ', a, ' tau= ', tau)
322 ocr = a*mon*np.exp(-a*mon*tau)
323 ocr_meas = ocr + 2*np.random.normal(size=len(ocr), scale=30.0)
324 (params,msg) = fit_deadtime(mon,ocr_meas)
325 tau = params[0]
326 a = params[1]
327 #print(msg)
328 print( 'a_fit= ',a,' tau_fit=', tau)
330 ocr = 0.3 * 1/tau
331 icr = calc_icr(ocr,tau)
332 print( 'max icr = ', 1/tau)
333 print( 'max ocr = ', np.exp(-1)/tau)
334 print( 'ocr= ', ocr, ' icr_calc= ',icr)
336 rt = 1.
337 lt = 1.
338 cor = correction_factor(rt,lt,icr,ocr)
339 print( 'cor= ', cor)