Coverage for larch/io/xsp3_hdf5.py: 9%

142 statements  

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

1#!/usr/bin/python 

2""" 

3support for netcdf file output files containing MCA spectra 

4from Epics Mapping Mode with XIA xMXAP electronics 

5""" 

6import numpy as np 

7import time 

8import h5py 

9import sys 

10import os 

11 

12from .. import Group 

13 

14# Default tau values for xspress3 

15XSPRESS3_TAU = 80.e-9 

16 

17def estimate_icr(ocr, tau, niter=3): 

18 "estimate icr from ocr and tau" 

19 maxicr = 1.0/tau 

20 maxocr = 1/(tau*np.exp(1.0)) 

21 ocr[np.where(ocr>2*maxocr)[0]] = 2*maxocr 

22 icr = 1.0*ocr 

23 for c in range(niter): 

24 delta = (icr - ocr*np.exp(icr*tau))/(icr*tau - 1) 

25 delta[np.where(delta < 0)[0]] = 0.0 

26 icr = icr + delta 

27 icr[np.where(icr>5*maxicr)[0]] = 5*maxicr 

28 #endfor 

29 return icr 

30#enddef 

31 

32 

33class XSP3Data(object): 

34 def __init__(self, npix, ndet, nchan): 

35 self.firstPixel = 0 

36 self.numPixels = 0 

37 self.realTime = np.zeros((npix, ndet), dtype='f8') 

38 self.liveTime = np.zeros((npix, ndet), dtype='f8') 

39 self.outputCounts = np.zeros((npix, ndet), dtype='f8') 

40 self.inputCounts = np.zeros((npix, ndet), dtype='f8') 

41 # self.counts = np.zeros((npix, ndet, nchan), dtype='f4') 

42 

43def get_counts_carefully(h5link): 

44 """ 

45 get counts array with some error checking for corrupted files, 

46 especially those that give 

47 'OSError: Can't read data (inflate() failed)' 

48 because one data point is bad. 

49 

50 This seems to be a common enough failure mode that this function looks 

51 for such points and replaces offending points by the average of the 

52 neighboring points 

53 """ 

54 # will usually succeed, of course. 

55 try: 

56 return h5link[()] 

57 except OSError: 

58 pass 

59 

60 # find bad point 

61 npts = h5link.shape[0] 

62 ilo, ihi = 0, npts 

63 imid = (ihi-ilo)//2 

64 while (ihi-ilo) > 1: 

65 try: 

66 _tmp = h5link[ilo:imid] 

67 ilo, imid = imid, ihi 

68 except OSError: 

69 ihi, imid = imid, (imid+ilo)//2 

70 

71 # retest to make sure there is one bad point 

72 for i in range(ilo, ihi): 

73 try: 

74 _tmp = h5link[i] 

75 except OSError: 

76 ibad = i 

77 break 

78 

79 # its not unusual to have two bad points in a row: 

80 p1bad = m1bad = False 

81 try: 

82 _tmp = h5link[i+1] 

83 except OSError: 

84 p1bad = True 

85 try: 

86 _tmp = h5link[i-1] 

87 except OSError: 

88 m1bad = True 

89 

90 if m1bad: 

91 ibad = ibad - 1 

92 p1bad = True 

93 

94 counts = np.zeros(h5link.shape, dtype=h5link.dtype) 

95 counts[:ibad] = h5link[:ibad] 

96 counts[ibad+2:] = h5link[ibad+2:] 

97 if p1bad: # two in a row 

98 print("fixing 2 bad points in h5 file") 

99 counts[ibad] = h5link[ibad-1] 

100 counts[ibad+1] = h5link[ibad+2] 

101 else: 

102 print("fixing bad point in h5 file") 

103 counts[ibad+1] = h5link[ibad+1] 

104 if ibad == 0: 

105 counts[ibad] = counts[ibad+1] 

106 elif ibad == npts - 1: 

107 counts[ibad] = counts[ibad-1] 

108 else: 

109 counts[ibad] = ((counts[ibad-1]+counts[ibad+1])/2.0).astype(h5link.dtype) 

110 return counts 

111 

112 

113def read_xsp3_hdf5(fname, npixels=None, verbose=False, 

114 estimate_dtc=False, **kws): 

115 # Reads a HDF5 file created with the Xspress3 driver 

116 npixels = None 

117 

118 clockrate = 12.5e-3 # microseconds per clock tick: 80MHz clock 

119 t0 = time.time() 

120 h5file = h5py.File(fname, 'r') 

121 

122 root = h5file['entry/instrument'] 

123 counts = get_counts_carefully(root['detector/data']) 

124 

125 # support bother newer and earlier location of NDAttributes 

126 ndattr = None 

127 try: 

128 ndattr = root['NDAttributes'] 

129 except KeyError: 

130 pass 

131 

132 if 'CHAN1SCA0' not in ndattr: 

133 try: 

134 ndattr = root['detector/NDAttributes'] 

135 except KeyError: 

136 pass 

137 if 'CHAN1SCA0' not in ndattr: 

138 raise ValueError("cannot find NDAttributes for '%s'" % fname) 

139 

140 # note: sometimes counts has npix-1 pixels, while the time arrays 

141 # really have npix... So we take npix from the time array, and 

142 npix = ndattr['CHAN1SCA0'].shape[0] 

143 ndpix, ndet, nchan = counts.shape 

144 if npixels is None: 

145 npixels = npix 

146 if npixels < ndpix: 

147 ndpix = npixels 

148 

149 out = XSP3Data(npixels, ndet, nchan) 

150 out.numPixels = npixels 

151 t1 = time.time() 

152 

153 if ndpix < npix: 

154 out.counts = np.zeros((npix, ndet, nchan), dtype='f8') 

155 out.counts[:ndpix, :, :] = counts 

156 else: 

157 out.counts = counts 

158 

159 if estimate_dtc: 

160 dtc_taus = [XSPRESS3_TAU]*ndet 

161 

162 for i in range(ndet): 

163 chan = "CHAN%i" %(i+1) 

164 clock_ticks = ndattr['%sSCA0' % chan][()] 

165 reset_ticks = ndattr["%sSCA1" % chan][()] 

166 all_events = ndattr["%sSCA3" % chan][()] 

167 if "%sEventWidth" in ndattr: 

168 event_width = 1.0 + ndattr['%sEventWidth' % chan][()] 

169 else: 

170 event_width = 6.0 

171 

172 clock_ticks[np.where(clock_ticks<10)] = 10.0 

173 rtime = clockrate * clock_ticks 

174 out.realTime[:, i] = rtime 

175 out.liveTime[:, i] = rtime 

176 ocounts = out.counts[:, i, 1:-1].sum(axis=1) 

177 ocounts[np.where(ocounts<0.1)] = 0.1 

178 out.outputCounts[:, i] = ocounts 

179 

180 denom = clock_ticks - (all_events*event_width + reset_ticks) 

181 denom[np.where(denom<2.0)] = 1.0 

182 dtfactor = clock_ticks/denom 

183 out.inputCounts[:, i] = dtfactor * ocounts 

184 

185 if estimate_dtc: 

186 ocr = ocounts/(rtime*1.e-6) 

187 icr = estimate_icr(ocr, dtc_taus[i], niter=3) 

188 out.inputCounts[:, i] = icr * (rtime*1.e-6) 

189 

190 h5file.close() 

191 t2 = time.time() 

192 if verbose: 

193 print(' time to read file = %5.1f ms' % ((t1-t0)*1000)) 

194 print(' time to extract data = %5.1f ms' % ((t2-t1)*1000)) 

195 print(' read %i pixels ' % npixels) 

196 print(' data shape: ' , out.counts.shape) 

197 return out 

198 

199def test_read(fname): 

200 print( fname, os.stat(fname)) 

201 fd = read_xsp3_hdf5(fname, verbose=True) 

202 print(fd.counts.shape)