Coverage for larch/xrmmap/gsexrm_utils.py: 12%

343 statements  

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

1import os 

2import numpy as np 

3import time 

4 

5import larch 

6 

7from larch.io import (read_xsp3_hdf5, read_xrf_netcdf, 

8 read_xrd_netcdf, read_xrd_hdf5) 

9from larch.utils.strutils import fix_varname 

10from .asciifiles import (readASCII, readMasterFile, readROIFile, 

11 readEnvironFile, read1DXRDFile, parseEnviron) 

12 

13from ..xrd import integrate_xrd_row 

14 

15def fix_xrd1d_filename(xrd_file): 

16 """check for 1D XRD file from Eiger or other detector -- 

17 avoids having to read hdf5 file at all 

18 """ 

19 # first append '.npy' 

20 xrd1d_file = xrd_file + '.npy' 

21 if os.path.exists(xrd1d_file): 

22 return xrd1d_file 

23 

24 # second, eiger .h5 -> .npy 

25 xrd1d_file = xrd_file.replace('.h5', '.npy').replace('_master', '') 

26 if os.path.exists(xrd1d_file): 

27 return xrd1d_file 

28 

29 return None 

30 

31def parse_sisnames(text): 

32 return [fix_varname(s.strip()) for s in text.replace('#', '').split('|')] 

33 

34class GSEXRM_FileStatus: 

35 no_xrfmap = 'hdf5 does not have top-level XRF map' 

36 no_xrdmap = 'hdf5 does not have top-level XRD map' 

37 created = 'hdf5 has empty schema' # xrm map exists, no data 

38 hasdata = 'hdf5 has map data' # array sizes known 

39 wrongfolder = 'hdf5 exists, but does not match folder name' 

40 err_notfound = 'file not found' 

41 empty = 'file is empty (read from folder)' 

42 err_nothdf5 = 'file is not hdf5 (or cannot be read)' 

43 err_nowrite = 'user does not have write permission' 

44 

45 

46class GSEXRM_Exception(Exception): 

47 '''GSEXRM Exception: General Errors''' 

48 pass 

49 

50 

51class GSEXRM_MCADetector(object): 

52 '''Detector class, representing 1 detector element (real or virtual) 

53 has the following properties (many of these as runtime-calculated properties) 

54 

55 rois list of ROI objects 

56 rois[i].name names 

57 rois[i].address address 

58 rois[i].left index of lower limit 

59 rois[i].right index of upper limit 

60 energy array of energy values 

61 counts array of count values 

62 dtfactor array of deadtime factor 

63 realtime array of real time 

64 livetime array of live time 

65 inputcounts array of input counts 

66 outputcount array of output count 

67 

68 ''' 

69 def __init__(self, xrmmap, prefix='mca', index=None): 

70 self.xrmmap = xrmmap 

71 self.prefix = prefix 

72 self.__ndet = xrmmap.attrs.get('N_Detectors', 0) 

73 self.det = None 

74 self.rois = [] 

75 detname = '%s1' % prefix 

76 if index is not None: 

77 detname = '%s%i' % (prefix, index) 

78 self.det = self.xrmmap[detname] 

79 

80 self.shape = self.xrmmap['%s/livetime' % detname].shape 

81 

82 # energy 

83 self.energy = self.xrmmap['%s/energy' % detname][()] 

84 

85 # set up rois 

86 rnames = self.xrmmap['%s/roi_names' % detname][()] 

87 raddrs = self.xrmmap['%s/roi_addrs' % detname][()] 

88 rlims = self.xrmmap['%s/roi_limits' % detname][()] 

89 for name, addr, lims in zip(rnames, raddrs, rlims): 

90 self.rois.append(ROI(name=name, address=addr, 

91 left=lims[0], right=lims[1])) 

92 

93 def __getval(self, param): 

94 if self.det is None: 

95 out = self.xrmmap['%s1/%s' % (self.prefix, param)][()] 

96 for i in range(2, self.__ndet): 

97 out += self.xrmmap['%s%i/%s' % (self.prefix, i, param)][()] 

98 return out 

99 return self.det[param][()] 

100 

101 @property 

102 def counts(self): 

103 "detector counts array" 

104 return self.__getval('counts') 

105 

106 @property 

107 def dtfactor(self): 

108 '''deadtime factor''' 

109 return self.__getval('dtfactor') 

110 

111 @property 

112 def realtime(self): 

113 '''real time''' 

114 return self.__getval('realtime') 

115 

116 @property 

117 def livetime(self): 

118 '''live time''' 

119 return self.__getval('livetime') 

120 

121 @property 

122 def inputcounts(self): 

123 '''inputcounts''' 

124 return self.__getval('inputcounts') 

125 

126 @property 

127 def outputcount(self): 

128 '''output counts''' 

129 return self.__getval('outputcounts') 

130 

131 

132class GSEXRM_Area(object): 

133 '''Map Area class, representing a map area for a detector 

134 ''' 

135 def __init__(self, xrmmap, index, det=None): 

136 self.xrmmap = xrmmap 

137 self.det = GSEXRM_Detector(xrmmap, index=det) 

138 if isinstance(index, int): 

139 index = 'area_%3.3i' % index 

140 self._area = self.xrmmap['areas/%s' % index] 

141 self.npts = self._area[()].sum() 

142 

143 sy, sx = [slice(min(_a), max(_a)+1) for _a in np.where(self._area)] 

144 self.yslice, self.xslice = sy, sx 

145 

146 def roicounts(self, roiname): 

147 iroi = -1 

148 for ir, roi in enumerate(self.det.rois): 

149 if roiname.lower() == roi.name.lower(): 

150 iroi = ir 

151 break 

152 if iroi < 0: 

153 raise ValueError('ROI name %s not found' % roiname) 

154 elo, ehi = self.det.rois[iroi].left, self.det.rois[iroi].right 

155 counts = self.det.counts[self.yslice, self.xslice, elo:ehi] 

156 

157 

158class GSEXRM_MapRow: 

159 ''' 

160 read one row worth of data: 

161 ''' 

162 def __init__(self, yvalue, xrffile, xrdfile, xpsfile, sisfile, folder, 

163 reverse=False, ixaddr=0, dimension=2, ioffset=0, 

164 npts=None, irow=None, dtime=None, nrows_expected=None, 

165 masterfile=None, xrftype=None, xrdtype=None, 

166 xrdcal=None, xrd2dmask=None, xrd2dbkgd=None, 

167 wdg=0, steps=4096, flip=True, force_no_dtc=False, 

168 has_xrf=True, has_xrd2d=False, has_xrd1d=False): 

169 

170 self.read_ok = False 

171 self.nrows_expected = nrows_expected 

172 

173 offslice = slice(None, None, None) 

174 if ioffset is None: 

175 ioffset = 0 

176 if ioffset > 0: 

177 offslice = slice(ioffset, None, None) 

178 elif ioffset < 0: 

179 offslice = slice(None, ioffset, None) 

180 self.npts = npts 

181 self.irow = irow 

182 if self.irow is None: 

183 self.irow = 0 

184 

185 self.yvalue = yvalue 

186 self.xrffile = xrffile 

187 self.xpsfile = xpsfile 

188 self.sisfile = sisfile 

189 self.xrdfile = xrdfile 

190 self.counts = None 

191 

192 self.xrd2d = None 

193 self.xrdq = None 

194 self.xrd1d = None 

195 self.xrdq_wdg = None 

196 self.xrd1d_wdg = None 

197 if masterfile is not None: 

198 header, rows = readMasterFile(masterfile) 

199 for row in header: 

200 if row.startswith('#XRF.filetype'): xrftype = row.split()[-1] 

201 if row.startswith('#XRD.filetype'): xrdtype = row.split()[-1] 

202 if has_xrf: 

203 if xrftype is None: 

204 xrftype = 'netcdf' 

205 if xrffile.startswith('xsp3'): 

206 xrftype = 'hdf5' 

207 if xrftype == 'netcdf': 

208 xrf_reader = read_xrf_netcdf 

209 else: 

210 xrf_reader = read_xsp3_hdf5 

211 

212 if has_xrd2d or has_xrd1d: 

213 if xrdtype == 'hdf5' or xrdfile.endswith('.h5'): 

214 xrd_reader = read_xrd_hdf5 

215 elif xrdtype == 'netcdf' or xrdfile.endswith('nc'): 

216 xrd_reader = read_xrd_netcdf 

217 else: 

218 xrd_reader = read_xrd_netcdf 

219 

220 # print( "xrd_reader ", xrd_reader, xrdtype, xrdfile, ' cal :%s: ' % xrdcal) 

221 # reading can fail with IOError, generally meaning the file isn't 

222 # ready for read. Try again for up to 5 seconds 

223 t0 = time.time() 

224 sis_ok, xps_ok = False, False 

225 

226 gdata, sdata = [], [] 

227 while not (sis_ok and xps_ok): 

228 try: 

229 ghead, gdata = readASCII(os.path.join(folder, xpsfile)) 

230 xps_ok = len(gdata) > 1 

231 except IOError: 

232 if (time.time() - t0) > 5.0: 

233 break 

234 time.sleep(0.25) 

235 try: 

236 shead, sdata = readASCII(os.path.join(folder, sisfile)) 

237 sdata = sdata[offslice] 

238 sis_ok = len(sdata) > 1 

239 except IOError: 

240 if (time.time() - t0) > 5.0: 

241 break 

242 time.sleep(0.25) 

243 

244 if not(sis_ok and xps_ok): 

245 print('Failed to read ASCII data for SIS: %s (%i), XPS: %s (%i)' % 

246 (sisfile, len(sdata), xpsfile, len(gdata)) ) 

247 return 

248 

249 # extrapolate gathering data by in case final end-point trigger was missed 

250 gather_extra = (2*gdata[-1] - gdata[-2]).reshape((1, gdata.shape[1])) 

251 gdata = np.concatenate((gdata, gather_extra)) 

252 gnpts, ngather = gdata.shape 

253 self.sishead = shead 

254 self.scaler_names = parse_sisnames(self.sishead[-1]) 

255 self.scaler_addrs = ['']*len(self.scaler_names) 

256 for sline in shead: 

257 if sline.startswith('# Column'): 

258 l, val = sline.split(': ') 

259 val = val.strip() 

260 label, addr, formula = [x.strip() for x in val.split('|')] 

261 if label in self.scaler_names: 

262 i = self.scaler_names.index(label) 

263 self.scaler_addrs[i] = addr 

264 

265 if dtime is not None: 

266 dtime.add('maprow: read ascii files') 

267 t0 = time.time() 

268 atime = -1 

269 xrf_dat, xrf_file = None, os.path.join(folder, xrffile) 

270 xrd_dat, xrd_file = None, os.path.join(folder, xrdfile) 

271 xrd1d_file = fix_xrd1d_filename(xrd_file) 

272 has_xrd1d = xrd1d_file is not None 

273 while atime < 0 and time.time()-t0 < 10: 

274 try: 

275 atime = os.stat(os.path.join(folder, sisfile)).st_ctime 

276 if has_xrf: 

277 xrf_dat = xrf_reader(xrf_file, npixels=self.nrows_expected, verbose=False) 

278 if xrf_dat is None: 

279 print( 'Failed to read XRF data from %s' % self.xrffile) 

280 if has_xrd2d or (has_xrd1d and xrd1d_file is None): 

281 xrd_dat = xrd_reader(xrd_file, verbose=False) 

282 if xrd_dat is None: 

283 print( 'Failed to read XRD data from %s' % self.xrdfile) 

284 

285 except (IOError, IndexError): 

286 time.sleep(0.025) 

287 

288 if atime < 0: 

289 print( 'Failed to read data.') 

290 return 

291 if dtime is not None: 

292 dtime.add('maprow: read XRM files') 

293 

294 ## SPECIFIC TO XRF data 

295 if has_xrf: 

296 self.counts = xrf_dat.counts[offslice] 

297 self.inpcounts = xrf_dat.inputCounts[offslice] 

298 self.outcounts = xrf_dat.outputCounts[offslice] 

299 

300 if self.inpcounts.max() < 1: 

301 self.inpcounts = self.counts.sum(axis=2) 

302 if self.outcounts.max() < 1: 

303 self.outcounts = self.inpcounts*1.0 

304 

305 self.livetime = xrf_dat.liveTime[offslice] 

306 self.realtime = xrf_dat.realTime[offslice] 

307 if self.livetime.max() < 0.01: 

308 self.livetime = 0.100 * np.ones(self.livetime.shape) 

309 if self.realtime.max() < 0.01: 

310 self.realtime = 0.100 * np.ones(self.realtime.shape) 

311 

312 dt_denom = self.outcounts*self.livetime 

313 dt_denom[np.where(dt_denom < 1)] = 1.0 

314 self.dtfactor = self.inpcounts*self.realtime/dt_denom 

315 self.dtfactor[np.where(np.isnan(self.dtfactor))] = 1.0 

316 self.dtfactor[np.where(self.dtfactor < 0.95)] = 0.95 

317 if force_no_dtc: # in case deadtime info is unreliable (some v old data) 

318 self.outcounts = self.inpcounts*1.0 

319 self.livetime = self.realtime*1.0 

320 self.dtfactor = np.ones(self.dtfactor.shape) 

321 

322 ## SPECIFIC TO XRD data 

323 if has_xrd2d or has_xrd1d: 

324 if has_xrd2d: 

325 if self.npts == xrd_dat.shape[0]: 

326 self.xrd2d = xrd_dat 

327 elif self.npts > xrd_dat.shape[0]: 

328 self.xrd2d = np.zeros((self.npts,xrd_dat.shape[1],xrd_dat.shape[2])) 

329 self.xrd2d[0:xrd_dat.shape[0]] = xrd_dat 

330 else: 

331 self.xrd2d = xrd_dat[0:self.npts] 

332 

333 ############################################################################ 

334 ## subtracts background and applies mask, row by row 

335 ## mkak 2018.02.01 

336 ## major speed up if no background or mask specified 

337 ## updated mkak 2018.03.30 

338 

339 if xrd2dmask is not None: 

340 dir = -1 if flip else 1 

341 mask2d = np.ones(self.xrd2d[0].shape) 

342 mask2d = mask2d - xrd2dmask[::dir] 

343 if xrd2dbkgd is not None: 

344 self.xrd2d = mask2d*(self.xrd2d-xrd2dbkgd) 

345 else: 

346 self.xrd2d = mask2d*(self.xrd2d) 

347 elif xrd2dbkgd is not None: 

348 self.xrd2d = self.xrd2d-xrd2dbkgd 

349 

350 ## limits all values to positive 

351 self.xrd2d[self.xrd2d < 0] = 0 

352 ############################################################################ 

353 

354 # print("read XRD1D file? ", irow, has_xrd1d, xrdcal, xrdcal is not None, xrd1d_file) 

355 if has_xrd1d and xrdcal is not None: 

356 attrs = dict(steps=steps, flip=flip) 

357 if xrd1d_file is not None: 

358 xdat = np.load(xrd1d_file) 

359 self.xrdq = xdat[0, :] 

360 self.xrd1d = xdat[1:, :] 

361 # print(">>>READ XRDQ 1D ", xrd1d_file, self.xrdq.shape, self.xrd1d.shape) 

362 if self.xrdq is None and self.xrd2d is not None: # integrate data if needed. 

363 print("will try to integrate 2DXRD data ", self.xrd2d.shape) 

364 # attrs['flip'] = True 

365 # self.xrd2d = self.xrd2d[:, 1:-1, 3:-3] 

366 # maxval = 2**32 - 2**14 

367 self.xrd2d[np.where(self.xrd2d>maxval)] = 0 

368 self.xrdq, self.xrd1d = integrate_xrd_row(self.xrd2d, xrdcal, 

369 **attrs) 

370 # print("Integrated to ", self.xrdq.shape) 

371 if wdg > 1: 

372 self.xrdq_wdg, self.xrd1d_wdg = [], [] 

373 wdg_sz = 360./int(wdg) 

374 for iwdg in range(wdg): 

375 wdg_lmts = np.array([iwdg*wdg_sz, (iwdg+1)*wdg_sz]) - 180 

376 attrs.update({'wedge_limits':wdg_lmts}) 

377 q, counts = integrate_xrd_row(self.xrd2d, xrdcal, **attrs) 

378 self.xrdq_wdg += [q] 

379 self.xrd1d_wdg += [counts] 

380 

381 self.xrdq_wdg = np.einsum('kij->ijk', self.xrdq_wdg) 

382 self.xrd1d_wdg = np.einsum('kij->ijk', self.xrd1d_wdg) 

383 

384 

385 xnpts, nmca = gnpts, 1 

386 if has_xrf: 

387 xnpts, nmca, nchan = self.counts.shape 

388 

389 snpts, nscalers = sdata.shape 

390 

391 # print("Row npts=%s, gather=%d, sis=%d, xrf=%d" % 

392 # (repr(self.npts), gnpts, snpts, xnpts)) 

393 

394 if self.npts is None: 

395 self.npts = min(gnpts, xnpts) 

396 

397 if snpts < self.npts: # extend struck data if needed 

398 print(' extending SIS data from %i to %i !' % (snpts, self.npts)) 

399 sdata = list(sdata) 

400 for i in range(self.npts+1-snpts): 

401 sdata.append(sdata[snpts-1]) 

402 sdata = np.array(sdata) 

403 snpts = self.npts 

404 self.sisdata = sdata[:self.npts] 

405 

406 if xnpts > self.npts: 

407 if has_xrf: 

408 self.counts = self.counts[:self.npts] 

409 self.realtime = self.realtime[:self.npts] 

410 self.livetime = self.livetime[:self.npts] 

411 self.dtfactor = self.dtfactor[:self.npts] 

412 self.inpcounts = self.inpcounts[:self.npts] 

413 self.outcounts = self.outcounts[:self.npts] 

414 if has_xrd2d: 

415 self.xrd2d = self.xrd2d[:self.npts] 

416 if has_xrd1d: 

417 # self.xrdq, self.xrd1d = self.xrdq[:self.npts], self.xrd1d[:self.npts] 

418 # print(" -- has xrd1d -> ", xnpts, self.npts, self.xrdq.shape) 

419 if self.xrdq_wdg is not None: 

420 self.xrdq_wdg = self.xrdq_wdg[:self.npts] 

421 self.xrd1d_wdg = self.xrd1d_wdg[:self.npts] 

422 

423 points = list(range(1, self.npts+1)) 

424 # auto-reverse: counter-intuitively (because stage is upside-down and so 

425 # backwards wrt optical view), left-to-right scans from high to low value 

426 # so reverse those that go from low to high value 

427 # first, find fast axis: 

428 ixaddr, gvarmax = 0, -1 

429 for ig in range(ngather): 

430 gvar = gdata[:, ig].std() 

431 if gvar > gvarmax: 

432 gvarmax = gvar 

433 ixaddr = ig 

434 

435 if reverse: 

436 do_reverse = gdata[0, ixaddr] < gdata[-1, ixaddr] 

437 else: 

438 do_reverse = gdata[0, ixaddr] > gdata[-1, ixaddr] 

439 

440 do_reverse = (irow % 2) == 1 

441 self.reversed = do_reverse 

442 if do_reverse: 

443 points.reverse() 

444 self.sisdata = self.sisdata[::-1] 

445 if has_xrf: 

446 self.counts = self.counts[::-1] 

447 self.realtime = self.realtime[::-1] 

448 self.livetime = self.livetime[::-1] 

449 self.dtfactor = self.dtfactor[::-1] 

450 self.inpcounts= self.inpcounts[::-1] 

451 self.outcounts= self.outcounts[::-1] 

452 if has_xrd2d and self.xrd2d is not None: 

453 self.xrd2d = self.xrd2d[::-1] 

454 if has_xrd1d and self.xrd1d is not None: 

455 self.xrd1d = self.xrd1d[::-1] 

456 if self.xrdq_wdg is not None: 

457 self.xrdq_wdg = self.xrdq_wdg[::-1] 

458 self.xrd1d_wdg = self.xrd1d_wdg[::-1] 

459 

460 if has_xrf: 

461 xvals = [(gdata[i, ixaddr] + gdata[i-1, ixaddr])/2.0 for i in points] 

462 self.posvals = [np.array(xvals)] 

463 if dimension == 2: 

464 self.posvals.append(np.array([float(yvalue) for i in points])) 

465 self.posvals.append(self.realtime.sum(axis=1).astype('float32') / nmca) 

466 self.posvals.append(self.livetime.sum(axis=1).astype('float32') / nmca) 

467 

468 self.dtfactor = self.dtfactor.astype('float32').swapaxes(0, 1) 

469 self.inpcounts= self.inpcounts.swapaxes(0, 1) 

470 self.outcounts= self.outcounts.swapaxes(0, 1) 

471 self.livetime = self.livetime.swapaxes(0, 1) 

472 self.realtime = self.realtime.swapaxes(0, 1) 

473 self.counts = self.counts.swapaxes(0, 1) 

474 iy, ix = self.dtfactor.shape 

475 

476 self.total = self.counts.sum(axis=0) 

477 # dtfactor for total 

478 total_dtc = (self.counts * self.dtfactor.reshape(iy, ix, 1)).sum(axis=0).sum(axis=1) 

479 dt_denom = self.total.sum(axis=1) 

480 dt_denom[np.where(dt_denom < 1)] = 1.0 

481 dtfact = total_dtc / dt_denom 

482 dtfact[np.where(np.isnan(dtfact))] = 1.0 

483 dtfact[np.where(dtfact < 0.95)] = 0.95 

484 dtfact[np.where(dtfact > 50.0)] = 50.0 

485 self.total_dtfactor = dtfact 

486 self.read_ok = True