Coverage for larch/xrf/xrf_calib.py: 9%

112 statements  

« prev     ^ index     » next       coverage.py v7.6.0, created at 2024-10-16 21:04 +0000

1#/usr/bin/env python 

2""" 

3XRF Calibration routines 

4""" 

5 

6import numpy as np 

7from lmfit.models import GaussianModel, ConstantModel 

8 

9from xraydb import xray_line 

10 

11from ..math import index_of, linregress, fit_peak 

12from .roi import split_roiname 

13from .mca import isLarchMCAGroup 

14 

15def xrf_calib_init_roi(mca, roiname): 

16 """initial calibration step for MCA: 

17 find energy locations for one ROI 

18 """ 

19 if not isLarchMCAGroup(mca): 

20 print( 'Not a valid MCA') 

21 return 

22 energy = 1.0*mca.energy 

23 chans = 1.0*np.arange(len(energy)) 

24 counts = mca.counts 

25 bgr = getattr(mca, 'bgr', None) 

26 if bgr is not None: 

27 counts = counts - bgr 

28 if not hasattr(mca, 'init_calib'): 

29 mca.init_calib = {} 

30 

31 roi = None 

32 for xroi in mca.rois: 

33 if xroi.name == roiname: 

34 roi = xroi 

35 break 

36 if roi is None: 

37 return 

38 words = roiname.split() 

39 elem = words[0].title() 

40 family = 'Ka' 

41 if len(words) > 1: 

42 family = words[1].title() 

43 if family == 'Lb': 

44 family = 'Lb1' 

45 try: 

46 eknown = xray_line(elem, family).energy/1000.0 

47 except: 

48 eknown = 0.001 

49 llim = max(0, roi.left - roi.bgr_width) 

50 hlim = min(len(chans)-1, roi.right + roi.bgr_width) 

51 segcounts = counts[llim:hlim] 

52 maxcounts = max(segcounts) 

53 ccen = llim + np.where(segcounts==maxcounts)[0][0] 

54 ecen = ccen * mca.slope + mca.offset 

55 bkgcounts = counts[llim] + counts[hlim] 

56 if maxcounts < 2*bkgcounts: 

57 mca.init_calib[roiname] = (eknown, ecen, 0.0, ccen, None) 

58 else: 

59 model = GaussianModel() + ConstantModel() 

60 params = model.make_params(amplitude=maxcounts, 

61 sigma=(chans[hlim]-chans[llim])/2.0, 

62 center=ccen-llim, c=0.00) 

63 params['center'].min = -10 

64 params['center'].max = hlim - llim + 10 

65 params['c'].min = -10 

66 out = model.fit(counts[llim:hlim], params, x=chans[llim:hlim]) 

67 ccen = llim + out.params['center'].value 

68 ecen = ccen * mca.slope + mca.offset 

69 fwhm = out.params['fwhm'].value * mca.slope 

70 mca.init_calib[roiname] = (eknown, ecen, fwhm, ccen, out) 

71 

72 

73def xrf_calib_fitrois(mca): 

74 """initial calibration step for MCA: 

75 find energy locations for all ROIs 

76 

77 """ 

78 if not isLarchMCAGroup(mca): 

79 print( 'Not a valid MCA') 

80 return 

81 

82 energy = 1.0*mca.energy 

83 chans = 1.0*np.arange(len(energy)) 

84 counts = mca.counts 

85 bgr = getattr(mca, 'bgr', None) 

86 if bgr is not None: 

87 counts = counts - bgr 

88 calib = {} 

89 for roi in mca.rois: 

90 words = roi.name.split() 

91 elem = words[0].title() 

92 family = 'ka' 

93 if len(words) > 1: 

94 family = words[1] 

95 try: 

96 eknown = xray_line(elem, family).energy/1000.0 

97 except: 

98 continue 

99 llim = max(0, roi.left - roi.bgr_width) 

100 hlim = min(len(chans)-1, roi.right + roi.bgr_width) 

101 fit = fit_peak(chans[llim:hlim], counts[llim:hlim], 

102 'Gaussian', background='constant') 

103 

104 

105 ccen = fit.params['center'].value 

106 ecen = ccen * mca.slope + mca.offset 

107 fwhm = 2.354820 * fit.params['sigma'].value * mca.slope 

108 calib[roi.name] = (eknown, ecen, fwhm, ccen, fit) 

109 mca.init_calib = calib 

110 

111def xrf_calib_compute(mca, apply=False): 

112 """compute linear energy calibration 

113 from init_calib dictionary found from xrf_calib_fitrois() 

114 

115 To exclude lines from the calibration, first run 

116 >>> xrf_calib_fitrois(mca) 

117 then remove items (by ROI name) from the mca.init_calib dictionay 

118 

119 """ 

120 if not isLarchMCAGroup(mca): 

121 print( 'Not a valid MCA') 

122 return 

123 if not hasattr(mca, 'init_calib'): 

124 xrf_calib_fitrois(mca) 

125 

126 # current calib 

127 offset, slope = mca.offset, mca.slope 

128 x = np.array([c[3] for c in mca.init_calib.values()]) 

129 y = np.array([c[0] for c in mca.init_calib.values()]) 

130 _s, _o, r, p, std = linregress(x, y) 

131 

132 mca.new_calib = (_o, _s) 

133 mca.cal_x = x 

134 mca.cal_y = y 

135 

136 if apply: 

137 xrf_calib_apply(mca, offset=_o, slope=_s) 

138 

139def xrf_calib_apply(mca, offset=None, slope=None): 

140 """apply calibration to MCA 

141 

142 either supply offset and slope arguments (in keV and keV/chan) 

143 or run xrf_calib_compute(mca) to estimate these from ROI peaks 

144 """ 

145 # print(" xrf calib apply ", mca, mca.offset, mca.slope, offset, slope) 

146 if not isLarchMCAGroup(mca): 

147 print( 'Not a valid MCA') 

148 return 

149 if (offset is None or slope is None) and not hasattr(mca, 'new_calib'): 

150 print( 'must supply offset and slope or run xrf_calib_compute()!') 

151 return 

152 

153 if (offset is None or slope is None): 

154 offset, slope = mca.new_calib 

155 mca.offset = offset 

156 mca.slope = slope 

157 npts = len(mca.energy) 

158 mca.energy = (offset + slope*np.arange(npts))