Coverage for larch/xrmmap/xrm_mapfile.py: 6%
2253 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
1import os
2import sys
3import uuid
4import time
5import h5py
6import numpy as np
7from pathlib import Path
8import scipy.stats as stats
9import json
10import multiprocessing as mp
11from functools import partial
13import larch
14from larch.utils import (isotime, fix_varname, fix_filename,
15 bytes2str, version_ge, unixpath)
17from larch.io import (new_filename, read_xrf_netcdf,
18 read_xsp3_hdf5, read_xrd_netcdf, read_xrd_hdf5)
20from larch.xrf import MCA, ROI
21from larch.xrd import (XRD, E_from_lambda, integrate_xrd_row, q_from_twth,
22 q_from_d, lambda_from_E, read_xrd_data, read_poni)
24from larch.math.tomography import tomo_reconstruction, reshape_sinogram, trim_sinogram
26from .configfile import FastMapConfig
27from .asciifiles import (readASCII, readMasterFile, readROIFile,
28 readEnvironFile, parseEnviron)
30from .gsexrm_utils import (GSEXRM_MCADetector, GSEXRM_Area, GSEXRM_Exception,
31 GSEXRM_MapRow, GSEXRM_FileStatus)
34DEFAULT_XRAY_ENERGY = 39987.0 # probably means x-ray energy was not found in meta data
35NINIT = 64
36COMPRESSION_OPTS = 2
37COMPRESSION = 'gzip'
38#COMPRESSION = 'lzf'
39DEFAULT_ROOTNAME = 'xrmmap'
40VALID_ROOTNAMES = ('xrmmap', 'xrfmap')
41EXTRA_DETGROUPS = ('scalars', 'work', 'xrd1d', 'xrd2d')
42NOT_OWNER = "Not Owner of HDF5 file %s"
43READ_ONLY = "HDF5 file %s is open read-only"
44QSTEPS = 2048
46H5ATTRS = {'Type': 'XRM 2D Map',
47 'Version': '2.1.0',
48 # 'Version': '2.0.0',
49 # 'Version': '1.0.1',
50 'Title': 'Epics Scan Data',
51 'Beamline': 'GSECARS, 13-IDE / APS',
52 'Start_Time': '',
53 'Stop_Time': '',
54 'Map_Folder': '',
55 'Dimension': 2,
56 'Process_Machine': '',
57 'Process_ID': 0,
58 'Compression': ''}
60def h5str(obj):
61 '''strings stored in an HDF5 from Python2 may look like
62 "b'xxx'", that is containg "b". strip these out here
63 '''
64 out = str(obj)
65 if out.startswith("b'") and out.endswith("'"):
66 out = out[2:-1]
67 return out
69def get_machineid():
70 "machine id / MAC address, independent of hostname"
71 return hex(uuid.getnode())[2:]
73def strlist(alist):
74 return [a.encode('utf-8') for a in alist]
76def isGSEXRM_MapFolder(fname):
77 "return whether folder a valid Scan Folder (raw data)"
78 if (fname is None or not Path(fname).exists() or
79 not Path(fname).is_dir()):
80 return False
81 flist = os.listdir(fname)
82 for f in ('Master.dat', 'Environ.dat', 'Scan.ini'):
83 if f not in flist:
84 return False
86 has_xrmdata = False
88 header, rows = readMasterFile(Path(fname, 'Master.dat'))
89 try:
90 for f in rows[0]:
91 if f in flist:
92 has_xrmdata = True
93 except:
94 pass
95 return has_xrmdata
98def test_h5group(group, folder=None):
99 "test h5 group as a XRM Map"
100 valid = ('config' in group and 'roimap' in group)
101 for attr in ('Version', 'Map_Folder', 'Dimension', 'Start_Time'):
102 valid = valid and attr in group.attrs
103 if not valid:
104 return None, None
105 status = GSEXRM_FileStatus.hasdata
106 vers = h5str(group.attrs.get('Version',''))
107 fullpath = group.attrs.get('Map_Folder','.')
108 _folder = Path(fullpath).name
109 if folder is not None and folder != _folder:
110 status = GSEXRM_FileStatus.wrongfolder
111 return status, vers
114def create_xrmmap(h5root, root=None, dimension=2, folder='', start_time=None):
115 '''creates a skeleton '/xrmmap' group in an open HDF5 file
117 This is left as a function, not method of GSEXRM_MapFile below
118 because it may be called by the mapping collection program
119 (ie, from collector.py) when a map is started
121 This leaves a structure to be filled in by
122 GSEXRM_MapFile.init_xrmmap(),
123 '''
124 attrs = {}
125 attrs.update(H5ATTRS)
126 if start_time is None:
127 start_time = time.ctime()
129 attrs.update({'Dimension':dimension, 'Start_Time':start_time,
130 'Map_Folder': folder, 'Last_Row': -1})
131 if root in ('', None):
132 root = DEFAULT_ROOTNAME
133 xrmmap = h5root.create_group(root)
135 for key, val in attrs.items():
136 xrmmap.attrs[key] = str(val)
138 g = xrmmap.create_group('roimap')
139 g.attrs['type'] = 'roi maps'
140 g.attrs['desc'] = 'ROI data, including summed and deadtime corrected maps'
142 g = xrmmap.create_group('config')
143 g.attrs['type'] = 'scan config'
144 g.attrs['desc'] = '''scan configuration, including scan definitions,
145 ROI definitions, MCA calibration, Environment Data, etc'''
147 g = xrmmap.create_group('areas')
148 g.attrs['type'] = 'areas'
150 g = xrmmap.create_group('positions')
151 g.attrs['type'] = 'position arrays'
153 g = xrmmap.create_group('scalars')
154 g.attrs['type'] = 'scalar detectors'
156 g = xrmmap.create_group('work')
157 g.attrs['type'] = 'virtual detectors for work/analysis arrays'
159 conf = xrmmap['config']
160 for name in ('scan', 'general', 'environ', 'positioners', 'notes',
161 'motor_controller', 'rois', 'mca_settings', 'mca_calib'):
162 conf.create_group(name)
163 h5root.flush()
165def ensure_subgroup(subgroup, group, dtype='virtual detector'):
166 if subgroup not in group.keys():
167 g = group.create_group(subgroup)
168 g.attrs['type'] = dtype
169 return g
170 else:
171 g = group[subgroup]
172 if 'type' not in g.attrs:
173 g.attrs['type'] = dtype
174 return g
176def toppath(pname, n=4):
177 words = []
178 for i in range(n):
179 words.append(Path(pname).name)
180 return '/'.join(words)
183def remove_zigzag(map, zigzag=0):
184 if zigzag == 0:
185 return map
186 nrows, ncols = map.shape
187 tmp = 1.0 * map
188 even = 0
189 if zigzag < 0:
190 even = 1
191 zigzag = -zigzag
192 for i in range(nrows):
193 if (i % 2) == even:
194 tmp[i, :] = map[i, :]
195 else:
196 tmp[i, zigzag:] = map[i, :-zigzag]
197 return tmp
200class GSEXRM_MapFile(object):
201 '''
202 Access to GSECARS X-ray Microprobe Map File:
204 The GSEXRM Map file is an HDF5 file built from a folder containing
205 'raw' data from a set of sources
206 xmap: XRF spectra saved to NetCDF by the Epics MCA detector
207 struck: a multichannel scaler, saved as ASCII column data
208 xps: stage positions, saved as ASCII file from the Newport XPS
210 The object here is intended to expose an HDF5 file that:
211 a) watches the corresponding folder and auto-updates when new
212 data is available, as for on-line collection
213 b) stores locking information (Machine Name/Process ID) in the top-level
215 For extracting data from a GSEXRM Map File, use:
217 >>> from epicscollect.io import GSEXRM_MapFile
218 >>> map = GSEXRM_MapFile('MyMap.001')
219 >>> fe = map.get_roimap('Fe', det='mca2')
220 >>> as = map.get_roimap('As Ka', det='mca1', dtcorrect=True)
221 >>> rgb = map.get_rgbmap('Fe', 'Ca', 'Zn', dtcorrect=True, scale_each=False)
222 >>> en = map.get_energy(det=1)
224 All these take the following options:
226 det: which detector element to use (1, 2, 3, 4, None), [None]
227 None means to use the sum of all detectors
228 dtcorrect: whether to return dead-time corrected spectra [True]
230 '''
232 ScanFile = 'Scan.ini'
233 EnvFile = 'Environ.dat'
234 ROIFile = 'ROI.dat'
235 XRDCALFile = 'XRD.poni'
236 MasterFile = 'Master.dat'
238 def __init__(self, filename=None, folder=None, create_empty=False,
239 hotcols=False, zigzag=0, dtcorrect=True, root=None,
240 chunksize=None, xrdcal=None, xrd2dmask=None, xrd2dbkgd=None,
241 xrd1dbkgd=None, azwdgs=0, qstps=QSTEPS, flip=True,
242 bkgdscale=1., has_xrf=True, has_xrd1d=False, has_xrd2d=False,
243 compression=COMPRESSION, compression_opts=COMPRESSION_OPTS,
244 facility='APS', beamline='13-ID-E', run='', proposal='',
245 user='', scandb=None, all_mcas=False, **kws):
247 self.filename = filename
248 self.folder = folder
249 self.root = root
250 self.chunksize = chunksize
251 # whether to remove first and last columns from data
252 self.hotcols = hotcols
253 # whether to shift rows to fix zig-zag
254 self.zigzag = zigzag
255 self.dtcorrect = dtcorrect
256 self.scandb = scandb
257 self.envvar = None
258 self.write_access = True
259 self.status = GSEXRM_FileStatus.err_notfound
260 self.dimension = None
261 self.nmca = None
262 self.npts = None
263 self.start_time = None
264 self.xrmmap = None
265 self.h5root = None
266 self.last_row = -1
267 self.rowdata = []
268 self.roi_names = {}
269 self.roi_slices = None
270 self._pixeltime = None
271 self.masterfile = None
272 self.force_no_dtc = False
273 self.all_mcas = all_mcas
274 self.detector_list = None
276 self.compress_args = {'compression': compression}
277 if compression != 'lzf':
278 self.compress_args['compression_opts'] = compression_opts
280 self.incident_energy = None
281 self.has_xrf = has_xrf
282 self.has_xrd1d = has_xrd1d
283 self.has_xrd2d = has_xrd2d
284 self.pos_desc = []
285 self.pos_addr = []
286 ## used for XRD
287 self.bkgd_xrd2d = None
288 self.bkgd_xrd1d = None
289 self.mask_xrd2d = None
290 self.xrdcalfile = None
291 self.xrd2dmaskfile = None
292 self.xrd2dbkgdfile = None
293 self.xrd1dbkgdfile = None
294 self.bkgdscale = bkgdscale if bkgdscale > 0 else 1.
295 self.azwdgs = 0 if azwdgs > 36 or azwdgs < 2 else int(azwdgs)
296 self.qstps = int(qstps)
297 self.flip = flip
298 self.master_modtime = -1
300 ## used for tomography orientation
301 self.x = None
302 self.reshape = None
303 self.notes = {'facility' : facility,
304 'beamline' : beamline,
305 'run' : run,
306 'proposal' : proposal,
307 'user' : user}
309 nmaster = -1
310 # initialize from filename or folder
311 if self.filename is not None:
312 self.getFileStatus(root=root)
313 # print("Get File Status ", root, self.status)
314 if self.status == GSEXRM_FileStatus.empty:
315 ftmp = open(self.filename, 'r')
316 self.folder = ftmp.readlines()[0][:-1].strip()
317 if '/' in self.folder:
318 self.folder = self.folder.split('/')[-1]
319 ftmp.close()
320 os.unlink(self.filename)
322 if (self.status==GSEXRM_FileStatus.err_notfound and
323 self.filename is not None and self.folder is None and
324 isGSEXRM_MapFolder(self.filename)):
325 self.folder = self.filename
326 self.filename = None
328 if isGSEXRM_MapFolder(self.folder):
329 nmaster = self.read_master()
330 if self.filename is None:
331 raise GSEXRM_Exception(
332 "'%s' is not a valid GSEXRM Map folder" % self.folder)
333 self.getFileStatus(root=root)
335 # for existing file, read initial settings
336 if self.status in (GSEXRM_FileStatus.hasdata,
337 GSEXRM_FileStatus.created):
338 self.open(self.filename, root=self.root, check_status=True)
339 self.reset_flags()
340 return
342 # file exists but is not hdf5
343 if self.status == GSEXRM_FileStatus.err_nothdf5:
344 raise GSEXRM_Exception(
345 "'%s' is not a readable HDF5 file" % self.filename)
347 # file has no write permission
348 # if self.status == GSEXRM_FileStatus.err_nowrite:
349 # raise GSEXRM_Exception(
350 # "'%s' does not have write access" % self.filename)
352 # create empty HDF5 if needed
353 if self.status == GSEXRM_FileStatus.empty and Path(self.filename).exists():
354 try:
355 flines = open(self.filename, 'r').readlines()
356 if len(flines) < 3:
357 os.unlink(self.filename)
358 self.status = GSEXRM_FileStatus.err_notfound
359 except (IOError, ValueError):
360 pass
362 if (self.status in (GSEXRM_FileStatus.err_notfound,
363 GSEXRM_FileStatus.wrongfolder) and
364 self.folder is not None and isGSEXRM_MapFolder(self.folder)):
366 if nmaster < 1:
367 self.read_master()
368 if self.status == GSEXRM_FileStatus.wrongfolder:
369 self.filename = new_filename(self.filename)
370 cfile = FastMapConfig()
371 cfile.Read(Path(self.folder, self.ScanFile))
372 cfile.config['scan']['filename'] = self.filename
373 print("Create HDF5 File ")
374 self.h5root = h5py.File(self.filename, 'w')
375 self.write_access = True
376 if self.dimension is None and isGSEXRM_MapFolder(self.folder):
377 if nmaster < 1:
378 self.read_master()
379 create_xrmmap(self.h5root, root=self.root, dimension=self.dimension,
380 folder=self.folder, start_time=self.start_time)
382 # cfile = FastMapConfig()
383 # self.add_map_config(cfile.config, nmca=nmca)
385 self.notes['h5_create_time'] = isotime()
386 self.status = GSEXRM_FileStatus.created
387 self.open(self.filename, root=self.root, check_status=False)
389 for xkey, xval in zip(self.xrmmap.attrs.keys(), self.xrmmap.attrs.values()):
390 if xkey == 'Version':
391 self.version = xval
393 self.add_XRDfiles(xrdcalfile=xrdcal,
394 xrd2dmaskfile=xrd2dmask,
395 xrd2dbkgdfile=xrd2dbkgd,
396 xrd1dbkgdfile=xrd1dbkgd)
397 elif (self.filename is not None and
398 self.status == GSEXRM_FileStatus.err_notfound and create_empty):
399 print("Create HDF5 File")
400 self.h5root = h5py.File(self.filename, 'w')
401 self.write_access = True
402 create_xrmmap(self.h5root, root=None, dimension=2, start_time=self.start_time)
403 self.notes['h5_create_time'] = isotime()
404 self.xrmmap = self.h5root[DEFAULT_ROOTNAME]
405 self.take_ownership()
406 # self.status = GSEXRM_FileStatus.created
407 cfile = FastMapConfig()
408 self.add_map_config(cfile.config, nmca=1)
409 self.open(self.filename, root=self.root, check_status=False)
410 self.reset_flags()
412 else:
413 raise GSEXRM_Exception('GSEXMAP Error: could not locate map file or folder')
414 print("Initialized done ", self.status, self.version, self.root)
416 def __repr__(self):
417 fname = '' if self.filename is None else self.filename
418 return f"GSEXRM_MapFile('{Path(self.filename).name}')"
421 def getFileStatus(self, filename=None, root=None, folder=None):
422 '''return status, top-level group, and version'''
424 if filename is not None:
425 self.filename = filename
426 filename = self.filename
427 folder = self.folder
428 # print("getFileStatus 0 ", filename, folder)
429 if folder is not None:
430 folder = Path(folder).absolute().name
431 self.status = GSEXRM_FileStatus.err_notfound
432 self.root, self.version = '', ''
433 if root not in ('', None):
434 self.root = root
435 # see if file exists:
436 if not (Path(filename).exists() and Path(filename).is_file()):
437 return
438 # see if file is empty/too small(signifies "read from folder")
439 if os.stat(filename).st_size < 1024:
440 self.status = GSEXRM_FileStatus.empty
441 return
443 if not os.access(filename, os.W_OK):
444 self.write_access = False
445 # print("getFileStatus 1 ", filename, self.write_access)
446 # see if file is an H5 file
447 try:
448 fh = h5py.File(filename, 'r')
449 except IOError:
450 self.status = GSEXRM_FileStatus.err_nothdf5
451 return
453 status = GSEXRM_FileStatus.no_xrfmap
454 if root is not None and root in fh:
455 stat, vers = test_h5group(fh[root], folder=folder)
456 if stat is not None:
457 self.status = stat
458 self.root, self.version = root, vers
459 else:
460 for root, group in fh.items():
461 stat, vers = test_h5group(group, folder=folder)
462 if stat is not None:
463 self.status = stat
464 self.root, self.version = root, vers
465 break
466 # print("getFileStatus 2 ", self.status)
467 fh.close()
468 return
470 def get_det(self, index):
471 return GSEXRM_MCADetector(self.xrmmap, index=index)
473 def area_obj(self, index, det=None):
474 return GSEXRM_Area(self.xrmmap, index, det=det)
476 def get_coarse_stages(self):
477 '''return coarse stage positions for map'''
478 stages = []
479 env_addrs = [h5str(s) for s in self.xrmmap['config/environ/address']]
480 env_vals = [h5str(s) for s in self.xrmmap['config/environ/value']]
481 for addr, pname in self.xrmmap['config/positioners'].items():
482 name = h5str(pname[()])
483 addr = h5str(addr)
484 val = ''
485 if not addr.endswith('.VAL'):
486 addr = '%s.VAL' % addr
487 if addr in env_addrs:
488 val = env_vals[env_addrs.index(addr)]
490 stages.append((addr, val, name))
492 return stages
494 def open(self, filename, root=None, check_status=True):
495 '''open GSEXRM HDF5 File :
496 with check_status=False, this **must** be called
497 for an existing, valid GSEXRM HDF5 File!!
498 '''
499 if root in ('', None):
500 root = DEFAULT_ROOTNAME
501 if check_status:
502 self.getFileStatus(filename, root=root)
503 if self.status not in (GSEXRM_FileStatus.hasdata,
504 GSEXRM_FileStatus.created):
505 raise GSEXRM_Exception(
506 "'%s' is not a valid GSEXRM HDF5 file" % self.filename)
507 self.filename = filename
508 # print("OPEN ", filename, self.write_access, self.h5root)
509 if self.h5root is None:
510 mode = 'a' if self.write_access else 'r'
511 try:
512 self.h5root = h5py.File(self.filename, mode)
513 except PermissionError:
514 self.write_access = False
515 self.h5root = h5py.File(self.filename, 'r')
516 print("Warning : file opened as read only")
517 self.xrmmap = self.h5root[root]
518 if self.folder is None:
519 self.folder = bytes2str(self.xrmmap.attrs.get('Map_Folder',''))
520 self.last_row = int(self.xrmmap.attrs.get('Last_Row',0))
522 try:
523 self.dimension = self.xrmmap['config/scan/dimension'][()]
524 except:
525 pass
527 if (len(self.rowdata) < 1 or
528 (self.dimension is None and isGSEXRM_MapFolder(self.folder))):
529 self.read_master()
531 if self.nmca is None:
532 self.nmca = self.xrmmap.attrs.get('N_Detectors', 1)
534 def close(self):
535 if self.check_hostid() and self.write_access:
536 self.xrmmap.attrs['Process_Machine'] = ''
537 self.xrmmap.attrs['Process_ID'] = 0
538 self.xrmmap.attrs['Last_Row'] = self.last_row
539 try:
540 self.h5root.close()
541 except RuntimeError:
542 print("Got Runtime Error ")
543 print(sys.exc_info())
545 self.h5root = None
547 def add_XRDfiles(self, flip=None, xrdcalfile=None, xrd2dmaskfile=None,
548 xrd2dbkgdfile=None, xrd1dbkgdfile=None):
549 '''
550 adds mask file to exisiting '/xrmmap' group in an open HDF5 file
551 mkak 2018.02.01
552 '''
553 xrd1dgrp = ensure_subgroup('xrd1d', self.xrmmap)
555 if xrdcalfile is not None:
556 self.xrdcalfile = xrdcalfile
557 if self.xrdcalfile is not None:
558 pcal = Path(self.xrdcalfile).absolute()
559 if pcal.exists():
560 self.xrdcalfile = pcal.as_posix()
561 print(f'Calibration file loaded: {self.xrdcalfile}')
562 xrd1dgrp.attrs['calfile'] = self.xrdcalfile
565 self.flip = flip if flip is not None else self.flip
567 if xrd1dbkgdfile is not None:
568 self.xrd1dbkgdfile= xrd1dbkgdfile
569 if self.xrd1dbkgdfile is not None:
570 if Path(self.xrd1dbkgdfile).exists():
571 print(f'xrd1d background file loaded: {self.xrd1dbkgdfile}')
572 xrd1dgrp.attrs['1Dbkgdfile'] = f'{self.xrd1dbkgdfile}'
573 self.bkgd_xrd1d = read_xrd_data(self.xrd1dbkgdfile)*self.bkgdscale
575 if xrd2dbkgdfile is not None:
576 self.xrd2dbkgdfile= xrd2dbkgdfile
577 if os.path.exists(str(self.xrd2dbkgdfile)):
578 print('2DXRD background file loaded: %s' % self.xrd2dbkgdfile)
579 xrd1dgrp.attrs['2Dbkgdfile'] = '%s' % (self.xrd2dbkgdfile)
580 self.bkgd_xrd2d = read_xrd_data(self.xrd2dbkgdfile)*self.bkgdscale
582 if xrd2dmaskfile is not None:
583 self.xrd2dmaskfile= xrd2dmaskfile
584 if os.path.exists(str(self.xrd2dmaskfile)):
585 print('Mask file loaded: %s' % self.xrd2dmaskfile)
586 xrd1dgrp.attrs['maskfile'] = '%s' % (self.xrd2dmaskfile)
587 self.mask_xrd2d = read_xrd_data(self.xrd2dmaskfile)
589 self.h5root.flush()
591 def add_data(self, group, name, data, attrs=None, **kws):
592 ''' creata an hdf5 dataset, replacing existing dataset if needed'''
593 if not self.check_hostid():
594 raise GSEXRM_Exception(NOT_OWNER % self.filename)
595 if not self.write_access:
596 raise GSEXRM_Exception(READ_ONLY % self.filename)
598 kws.update(self.compress_args)
599 if name in group:
600 del group[name]
601 d = group.create_dataset(name, data=data, **kws)
602 if isinstance(attrs, dict):
603 for key, val in attrs.items():
604 d.attrs[key] = val
605 return d
607 def add_map_config(self, config, nmca=None):
608 '''add configuration from Map Folder to HDF5 file
609 ROI, DXP Settings, and Config data
610 '''
611 if not self.check_hostid():
612 raise GSEXRM_Exception(NOT_OWNER % self.filename)
613 if not self.write_access:
614 raise GSEXRM_Exception(READ_ONLY % self.filename)
616 group = self.xrmmap['config']
617 for name, sect in (('scan', 'scan'),
618 ('general', 'general'),
619 ('positioners', 'slow_positioners'),
620 ('motor_controller', 'xps')):
621 for key, val in config[sect].items():
622 if key in group[name]:
623 del group[name][key]
624 group[name].create_dataset(key, data=val)
626 scanfile = os.path.join(self.folder, self.ScanFile)
627 if os.path.exists(scanfile):
628 scantext = open(scanfile, 'r').read()
629 else:
630 scantext = ''
631 if 'text' in group['scan']:
632 del group['scan']['text']
633 group['scan'].create_dataset('text', data=scantext)
635 roifile = os.path.join(self.folder, self.ROIFile)
636 self.nmca = 0
637 if nmca is not None:
638 self.nmca = nmca
640 if os.path.exists(roifile):
641 roidat, calib, extra = readROIFile(roifile)
642 self.xrmmap.attrs['N_Detectors'] = self.nmca = len(calib['slope'])
643 roi_desc, roi_addr, roi_lim = [], [], []
644 roi_slices = []
646 for iroi, label, lims in roidat:
647 roi_desc.append(label)
648 try:
649 xrf_prefix = config['xrf']['prefix']
650 except KeyError: # very old map files
651 xrf_prefix = 'xrf_det:'
652 roi_addr.append("%smca%%i.R%i" % (xrf_prefix, iroi))
653 roi_lim.append([lims[i] for i in range(self.nmca)])
654 roi_slices.append([slice(lims[i][0], lims[i][1]) for i in range(self.nmca)])
655 roi_lim = np.array(roi_lim)
656 self.add_data(group['rois'], 'name', strlist(roi_desc))
657 self.add_data(group['rois'], 'address', strlist(roi_addr))
658 self.add_data(group['rois'], 'limits', roi_lim)
660 for key, val in calib.items():
661 self.add_data(group['mca_calib'], key, val)
663 for key, val in extra.items():
664 try:
665 self.add_data(group['mca_settings'], key, val)
666 except TypeError:
667 pass
669 self.roi_desc = roi_desc
670 self.roi_addr = roi_addr
671 self.roi_slices = roi_slices
672 self.calib = calib
673 else:
674 nmca = self.nmca
675 roilims = np.array([ [[0, 1]] for i in range(nmca)])
676 self.add_data(group['rois'], 'name', [b'roi1']*nmca)
677 self.add_data(group['rois'], 'address', [b'addr%i.ROI1']*nmca)
678 self.add_data(group['rois'], 'limits', roilims)
679 self.add_data(group['mca_calib'], 'offset', [0.00]*nmca)
680 self.add_data(group['mca_calib'], 'slope', [0.01]*nmca)
681 self.add_data(group['mca_calib'], 'quad', [0.00]*nmca)
683 # add env data
684 envfile = os.path.join(self.folder, self.EnvFile)
685 if os.path.exists(envfile):
686 envdat = readEnvironFile(envfile)
687 else:
688 envdat = ['Facility.Ring_Current (UnknownPV) = 0']
689 env_desc, env_addr, env_val = parseEnviron(envdat)
691 self.add_data(group['environ'], 'name', strlist(env_desc))
692 self.add_data(group['environ'], 'address', strlist(env_addr))
693 self.add_data(group['environ'], 'value', strlist(env_val))
695 cmprstr = '%s' % self.compress_args['compression']
696 if self.compress_args['compression'] != 'lzf':
697 cmprstr = '%s-%s' % (cmprstr,self.compress_args['compression_opts'])
698 self.xrmmap.attrs['Compression'] = cmprstr
700 self.h5root.flush()
702 def initialize_xrmmap(self, callback=None):
703 ''' initialize '/xrmmap' group in HDF5 file, generally
704 possible once at least 1 row of raw data is available
705 in the scan folder.
706 '''
707 self.starttime = time.time()
708 if self.status == GSEXRM_FileStatus.hasdata:
709 return
710 if self.status != GSEXRM_FileStatus.created:
711 print( 'Warning, cannot initialize xrmmap yet.')
712 return
714 if not self.check_hostid():
715 raise GSEXRM_Exception(NOT_OWNER % self.filename)
716 if not self.write_access:
717 raise GSEXRM_Exception(READ_ONLY % self.filename)
719 if (len(self.rowdata) < 1 or
720 (self.dimension is None and isGSEXRM_MapFolder(self.folder))):
721 self.read_master()
722 if len(self.rowdata) < 1:
723 return
725 self.last_row = -1
726 self.add_map_config(self.mapconf)
728 self.process_row(0, flush=True, callback=None,
729 nrows_expected=len(self.rowdata))
731 self.status = GSEXRM_FileStatus.hasdata
733 def process_row(self, irow, flush=False, complete=False, offset=None,
734 nrows_expected=None, callback=None):
735 row = self.read_rowdata(irow, offset=offset)
736 if irow == 0:
737 nmca, nchan = 0, 2048
738 if row.counts is not None:
739 nmca, xnpts, nchan = row.counts.shape
740 xrd2d_shape = None
741 if row.xrd2d is not None:
742 xrd2d_shape = rows.xrd2d.shape
743 self.build_schema(row.npts, nmca=nmca, nchan=nchan,
744 scaler_names=row.scaler_names,
745 scaler_addrs=row.scaler_addrs,
746 xrd2d_shape=xrd2d_shape, verbose=True,
747 nrows_expected=nrows_expected)
748 if row.read_ok:
749 self.add_rowdata(row, callback=callback)
751 if flush or complete:
752 # print("Flush, ", irow, self.last_row, flush, complete)
753 self.resize_arrays(self.last_row+1, force_shrink=True)
754 self.h5root.flush()
755 if self._pixeltime is None:
756 self.calc_pixeltime()
757 if callable(callback):
758 status = 'complete' if complete else 'flush'
759 callback(filename=self.filename, status=status)
762 def process(self, maxrow=None, force=False, callback=None, offset=None,
763 force_no_dtc=False, all_mcas=None):
764 "look for more data from raw folder, process if needed"
765 self.force_no_dtc = force_no_dtc
766 if all_mcas is not None:
767 self.all_mcas = all_mcas
769 if not self.check_hostid():
770 raise GSEXRM_Exception(NOT_OWNER % self.filename)
771 if not self.write_access:
772 raise GSEXRM_Exception(READ_ONLY % self.filename)
774 self.reset_flags()
775 if self.status == GSEXRM_FileStatus.created:
776 self.initialize_xrmmap(callback=callback)
777 if (force or len(self.rowdata) < 1 or
778 (self.dimension is None and isGSEXRM_MapFolder(self.folder))):
779 self.read_master()
781 nrows = len(self.rowdata)
782 if maxrow is not None:
783 nrows = min(nrows, maxrow)
785 if force or self.folder_has_newdata():
786 irow = self.last_row + 1
787 while irow < nrows:
788 flush = irow < 2 or (irow % 64 == 0)
789 complete = irow >= nrows-1
790 self.process_row(irow, flush=flush, offset=offset,
791 complete=complete, callback=callback)
792 irow = irow + 1
793 if callable(callback):
794 callback(filename=self.filename, status='complete')
797 def set_roidata(self, row_start=0, row_end=None):
798 if row_end is None:
799 row_end = self.last_row
801 # print("Process ROIs for rows %d to %d " % (row_start+1, row_end+1))
802 rows = slice(row_start, row_end+1)
803 roigrp = self.xrmmap['roimap']
804 conf = self.xrmmap['config']
805 roi_names = [h5str(s) for s in conf['rois/name']]
806 roi_limits = conf['rois/limits'][()]
807 # print("roi names ", roi_names, roi_limits)
809 for roiname, lims in zip(roi_names, roi_limits):
810 # dt = debugtimer()
811 roi_slice = lims[0]
812 # dt.add('get slice')
813 sumraw = roigrp['mcasum'][roiname]['raw'][rows,]
814 # dt.add('get sum raw')
815 sumcor = roigrp['mcasum'][roiname]['cor'][rows,]
816 # dt.add('get sum cor')
818 for detname in self.mca_dets:
819 mcaraw = self.xrmmap[detname]['counts'][rows,:,roi_slice].sum(axis=2)
820 # dt.add(" mcaraw %s " % detname)
821 mcacor = mcaraw*self.xrmmap[detname]['dtfactor'][rows,:]
822 # dt.add(" mcacor %s " % detname)
823 roigrp[detname][roiname]['raw'][rows,] = mcaraw
824 roigrp[detname][roiname]['cor'][rows,] = mcacor
825 # dt.add(" set roigrps for %s " % detname)
826 sumraw += mcaraw
827 sumcor += mcacor
828 # dt.add(" sum %s " % detname)
829 #dt.show()
831 roigrp['mcasum'][roiname]['raw'][rows,] = sumraw
832 roigrp['mcasum'][roiname]['cor'][rows,] = sumcor
833 self.h5root.flush()
835 def calc_pixeltime(self):
836 scanconf = self.xrmmap['config/scan']
837 rowtime = float(scanconf['time1'][()])
838 start = float(scanconf['start1'][()])
839 stop = float(scanconf['stop1'][()])
840 step = float(scanconf['step1'][()])
841 npts = int((abs(stop - start) + 1.1*step)/step)
842 self._pixeltime = rowtime/(npts-1)
843 return self._pixeltime
845 @property
846 def pixeltime(self):
847 """Return the pixel time"""
848 if self._pixeltime is None:
849 self.calc_pixeltime()
850 return self._pixeltime
852 def read_rowdata(self, irow, offset=None):
853 '''read a row worth of raw data from the Map Folder
854 returns arrays of data
855 '''
856 if self.dimension is None or irow > len(self.rowdata):
857 self.read_master()
859 if self.folder is None or irow >= len(self.rowdata):
860 return
862 if self.has_xrd1d and self.xrdcalfile is None:
863 self.xrdcalfile = bytes2str(self.xrmmap['xrd1d'].attrs.get('calfile',''))
864 if self.xrdcalfile in (None, ''):
865 calfile = os.path.join(unixpath(self.folder), self.XRDCALFile)
866 if os.path.exists(calfile):
867 self.xrdcalfile = calfile
869 scan_version = getattr(self, 'scan_version', 1.00)
870 print(" read row data, scan version ", scan_version, self.xrdcalfile)
871 xrdcal_dat = bytes2str(self.xrmmap['xrd1d'].attrs.get('caldata','{}'))
872 if self.xrdcalfile is not None and len(xrdcal_dat) < 10:
873 xrdcal_dat = read_poni(self.xrdcalfile)
874 self.xrmmap['xrd1d'].attrs['caldata'] = json.dumps(xrdcal_dat)
876 # if not self.has_xrf and not self.has_xrd2d and not self.has_xrd1d:
877 # raise IOError('No XRF or XRD flags provided.')
878 # return
880 if scan_version > 1.35 or self.has_xrd2d or self.has_xrd1d:
881 yval, xrff, sisf, xpsf, xrdf, etime = self.rowdata[irow]
882 if xrff.startswith('None'):
883 xrff = xrff.replace('None', 'xsp3')
884 if sisf.startswith('None'):
885 sisf = sisf.replace('None', 'struck')
886 if xpsf.startswith('None'):
887 xpsf = xpsf.replace('None', 'xps')
888 if xrdf.startswith('None'):
889 xrdf = xrdf.replace('None', 'pexrd')
890 else:
891 yval, xrff, sisf, xpsf, etime = self.rowdata[irow]
892 xrdf = '_unused_'
894 if '_unused_' in xrdf:
895 self.has_xrd1d = False
896 self.has_xrd2d = False
898 # eiger XRD maps with 1D data
899 if (xrdf.startswith('eig') and xrdf.endswith('.h5') or
900 xrdf.startswith('pexrd')):
901 self.has_xrd2d = False
902 self.has_xrd1d = True
904 if '_unused_' in xrff:
905 self.has_xrf = False
907 reverse = True
908 ioffset = 0
909 if scan_version > 1.35:
910 ioffset = 1
911 if scan_version >= 2.0:
912 reverse = False
913 if offset is not None:
914 ioffset = offset
915 self.has_xrf = self.has_xrf and xrff != '_unused_'
916 return GSEXRM_MapRow(yval, xrff, xrdf, xpsf, sisf, self.folder,
917 irow=irow, nrows_expected=self.nrows_expected,
918 ixaddr=0, dimension=self.dimension,
919 npts=self.npts,
920 reverse=reverse,
921 ioffset=ioffset,
922 force_no_dtc=self.force_no_dtc,
923 masterfile=self.masterfile, flip=self.flip,
924 xrdcal=self.xrdcalfile,
925 xrd2dmask=self.mask_xrd2d,
926 xrd2dbkgd=self.bkgd_xrd2d, wdg=self.azwdgs,
927 steps=self.qstps, has_xrf=self.has_xrf,
928 has_xrd2d=self.has_xrd2d,
929 has_xrd1d=self.has_xrd1d)
932 def add_rowdata(self, row, callback=None, flush=True):
933 '''adds a row worth of real data'''
934 # dt = debugtimer()
935 if not self.check_hostid():
936 raise GSEXRM_Exception(NOT_OWNER % self.filename)
937 if not self.write_access:
938 raise GSEXRM_Exception(READ_ONLY % self.filename)
940 thisrow = self.last_row + 1
941 if hasattr(callback, '__call__'):
942 callback(row=(thisrow+1), maxrow=len(self.rowdata), filename=self.filename)
944 pform = 'Add row %4i, yval=%s' % (thisrow+1, row.yvalue)
945 if self.has_xrf:
946 pform = '%s, xrffile=%s' % (pform, row.xrffile)
947 if self.has_xrd2d or self.has_xrd1d:
948 pform = '%s, xrdfile=%s' % (pform, row.xrdfile)
949 print(pform)
951 # dt.add(" ran callback, print, version %s" %self.version)
953 if version_ge(self.version, '2.0.0'):
955 mcasum_raw,mcasum_cor = [],[]
956 nrows = 0
957 map_items = sorted(self.xrmmap.keys())
958 # dt.add(" get %d map items" % len(map_items))
959 for gname in map_items:
960 g = self.xrmmap[gname]
961 if bytes2str(g.attrs.get('type', '')).startswith('scalar detect'):
962 first_det = list(g.keys())[0]
963 nrows, npts = g[first_det].shape
965 # dt.add(" got %d map items" % len(map_items))
966 if thisrow >= nrows:
967 self.resize_arrays(NINIT*(1+nrows/NINIT), force_shrink=False)
969 # dt.add(" resized ")
970 sclrgrp = self.xrmmap['scalars']
971 for ai, aname in enumerate(row.scaler_names):
972 sclrgrp[aname][thisrow, :npts] = row.sisdata[:npts].transpose()[ai]
973 # dt.add(" add scaler group")
974 if self.has_xrf:
976 npts = min([len(p) for p in row.posvals])
977 pos = self.xrmmap['positions/pos']
978 rowpos = np.array([p[:npts] for p in row.posvals])
980 tpos = rowpos.transpose()
981 pos[thisrow, :npts, :] = tpos[:npts, :]
982 nmca, xnpts, nchan = row.counts.shape
983 mca_dets = []
984 # dt.add(" map xrf 1")
985 for gname in map_items:
986 g = self.xrmmap[gname]
987 if bytes2str(g.attrs.get('type', '')).startswith('mca detect'):
988 mca_dets.append(gname)
989 nrows, npts, nchan = g['counts'].shape
990 # dt.add(" map xrf 2")
991 _nr, npts, nchan = self.xrmmap['mcasum']['counts'].shape
992 npts = min(npts, xnpts, self.npts)
993 # dt.add(" map xrf 3")
994 # print("ADD ROW ", self.all_mcas, mca_dets, self.nmca)
995 if self.all_mcas:
996 for idet, gname in enumerate(mca_dets):
997 grp = self.xrmmap[gname]
998 grp['counts'][thisrow, :npts, :] = row.counts[idet, :npts, :]
999 grp['dtfactor'][thisrow, :npts] = row.dtfactor[idet, :npts]
1000 grp['realtime'][thisrow, :npts] = row.realtime[idet, :npts]
1001 grp['livetime'][thisrow, :npts] = row.livetime[idet, :npts]
1002 grp['inpcounts'][thisrow, :npts] = row.inpcounts[idet, :npts]
1003 grp['outcounts'][thisrow, :npts] = row.outcounts[idet, :npts]
1005 livetime = np.zeros(npts, dtype=np.float64)
1006 realtime = np.zeros(npts, dtype=np.float64)
1007 dtfactor = np.zeros(npts, dtype=np.float32)
1008 inpcounts = np.zeros(npts, dtype=np.float32)
1009 outcounts = np.zeros(npts, dtype=np.float32)
1010 # dt.add(" map xrf 4a: alloc ")
1011 # print("ADD ", inpcounts.dtype, row.inpcounts.dtype)
1012 for idet in range(self.nmca):
1013 realtime += row.realtime[idet, :npts]
1014 livetime += row.livetime[idet, :npts]
1015 inpcounts += row.inpcounts[idet, :npts]
1016 outcounts += row.outcounts[idet, :npts]
1017 livetime /= (1.0*self.nmca)
1018 realtime /= (1.0*self.nmca)
1019 # dt.add(" map xrf 4b: time sums")
1021 sumgrp = self.xrmmap['mcasum']
1022 sumgrp['counts'][thisrow, :npts, :nchan] = row.total[:npts, :nchan]
1023 # dt.add(" map xrf 4b: set counts")
1024 # print("add realtime ", sumgrp['realtime'].shape, self.xrmmap['roimap/det_raw'].shape, thisrow)
1025 sumgrp['realtime'][thisrow, :npts] = realtime
1026 sumgrp['livetime'][thisrow, :npts] = livetime
1027 sumgrp['dtfactor'][thisrow, :npts] = row.total_dtfactor[:npts]
1028 sumgrp['inpcounts'][thisrow, :npts] = inpcounts
1029 sumgrp['outcounts'][thisrow, :npts] = outcounts
1030 # dt.add(" map xrf 4c: set time data ")
1032 if version_ge(self.version, '2.1.0'): # version 2.1
1033 det_raw = self.xrmmap['roimap/det_raw']
1034 det_cor = self.xrmmap['roimap/det_cor']
1035 sum_raw = self.xrmmap['roimap/sum_raw']
1036 sum_cor = self.xrmmap['roimap/sum_cor']
1038 detraw = list(row.sisdata[:npts].transpose())
1039 detcor = detraw[:]
1040 sumraw = detraw[:]
1041 sumcor = detraw[:]
1043 if self.roi_slices is None:
1044 lims = self.xrmmap['config/rois/limits'][()]
1045 nrois, nmca, nx = lims.shape
1047 self.roi_slices = []
1048 for iroi in range(nrois):
1049 x = [slice(lims[iroi, i, 0],
1050 lims[iroi, i, 1]) for i in range(nmca)]
1051 self.roi_slices.append(x)
1053 for slices in self.roi_slices:
1054 iraw = [row.counts[i, :npts, slices[i]].sum(axis=1)
1055 for i in range(nmca)]
1056 icor = [row.counts[i, :npts, slices[i]].sum(axis=1)*row.dtfactor[i, :npts]
1057 for i in range(nmca)]
1058 detraw.extend(iraw)
1059 detcor.extend(icor)
1060 sumraw.append(np.array(iraw).sum(axis=0))
1061 sumcor.append(np.array(icor).sum(axis=0))
1062 # dt.add(" map xrf 5a: got simple ROIS")
1063 det_raw[thisrow, :npts, :] = np.array(detraw).transpose()
1064 det_cor[thisrow, :npts, :] = np.array(detcor).transpose()
1065 sum_raw[thisrow, :npts, :] = np.array(sumraw).transpose()
1066 sum_cor[thisrow, :npts, :] = np.array(sumcor).transpose()
1069 else: # version 2.0
1070 roigrp = self.xrmmap['roimap']
1071 en = self.xrmmap['mcasum']['energy'][:]
1072 for roiname in roigrp['mcasum'].keys():
1073 en_lim = roigrp['mcasum'][roiname]['limits'][:]
1074 roi_slice = slice(np.abs(en-en_lim[0]).argmin(),
1075 np.abs(en-en_lim[1]).argmin())
1076 sumraw = roigrp['mcasum'][roiname]['raw'][thisrow,]
1077 sumcor = roigrp['mcasum'][roiname]['cor'][thisrow,]
1078 for detname in mca_dets:
1079 mcaraw = self.xrmmap[detname]['counts'][thisrow,][:,roi_slice].sum(axis=1)
1080 mcacor = mcaraw*self.xrmmap[detname]['dtfactor'][thisrow,]
1081 roigrp[detname][roiname]['raw'][thisrow,] = mcaraw
1082 roigrp[detname][roiname]['cor'][thisrow,] = mcacor
1083 sumraw += mcaraw
1084 sumcor += mcacor
1085 roigrp['mcasum'][roiname]['raw'][thisrow,] = sumraw
1086 roigrp['mcasum'][roiname]['cor'][thisrow,] = sumcor
1087 # dt.add(" map xrf 6")
1088 else: # version 1.0.1
1089 if self.has_xrf:
1090 nmca, xnpts, nchan = row.counts.shape
1091 xrm_dets = []
1093 nrows = 0
1094 map_items = sorted(self.xrmmap.keys())
1095 for gname in map_items:
1096 g = self.xrmmap[gname]
1097 if bytes2str(g.attrs.get('type', '')).startswith('mca detect'):
1098 xrm_dets.append(g)
1099 nrows, npts, nchan = g['counts'].shape
1101 if thisrow >= nrows:
1102 self.resize_arrays(NINIT*(1+nrows/NINIT), force_shrink=False)
1104 _nr, npts, nchan = xrm_dets[0]['counts'].shape
1105 npts = min(npts, xnpts, self.npts)
1106 for idet, grp in enumerate(xrm_dets):
1107 grp['dtfactor'][thisrow, :npts] = row.dtfactor[idet, :npts]
1108 grp['realtime'][thisrow, :npts] = row.realtime[idet, :npts]
1109 grp['livetime'][thisrow, :npts] = row.livetime[idet, :npts]
1110 grp['inpcounts'][thisrow, :npts] = row.inpcounts[idet, :npts]
1111 grp['outcounts'][thisrow, :npts] = row.outcounts[idet, :npts]
1112 grp['counts'][thisrow, :npts, :] = row.counts[idet, :npts, :]
1114 # here, we add the total dead-time-corrected data to detsum.
1115 self.xrmmap['detsum']['counts'][thisrow, :npts, :nchan] = row.total[:npts, :nchan]
1117 pos = self.xrmmap['positions/pos']
1118 rowpos = np.array([p[:npts] for p in row.posvals])
1120 tpos = rowpos.transpose()
1122 pos[thisrow, :npts, :] = tpos[:npts, :]
1124 # now add roi map data
1125 roimap = self.xrmmap['roimap']
1126 det_raw = roimap['det_raw']
1127 det_cor = roimap['det_cor']
1128 sum_raw = roimap['sum_raw']
1129 sum_cor = roimap['sum_cor']
1131 detraw = list(row.sisdata[:npts].transpose())
1133 detcor = detraw[:]
1134 sumraw = detraw[:]
1135 sumcor = detraw[:]
1137 if self.roi_slices is None:
1138 lims = self.xrmmap['config/rois/limits'][()]
1139 nrois, nmca, nx = lims.shape
1141 self.roi_slices = []
1142 for iroi in range(nrois):
1143 x = [slice(lims[iroi, i, 0],
1144 lims[iroi, i, 1]) for i in range(nmca)]
1145 self.roi_slices.append(x)
1147 for slices in self.roi_slices:
1148 iraw = [row.counts[i, :npts, slices[i]].sum(axis=1)
1149 for i in range(nmca)]
1150 icor = [row.counts[i, :npts, slices[i]].sum(axis=1)*row.dtfactor[i, :npts]
1151 for i in range(nmca)]
1152 detraw.extend(iraw)
1153 detcor.extend(icor)
1154 sumraw.append(np.array(iraw).sum(axis=0))
1155 sumcor.append(np.array(icor).sum(axis=0))
1157 det_raw[thisrow, :npts, :] = np.array(detraw).transpose()
1158 det_cor[thisrow, :npts, :] = np.array(detcor).transpose()
1159 sum_raw[thisrow, :npts, :] = np.array(sumraw).transpose()
1160 sum_cor[thisrow, :npts, :] = np.array(sumcor).transpose()
1162 if self.has_xrd1d and row.xrdq is not None:
1163 if thisrow < 2:
1164 if len(row.xrdq.shape) == 1:
1165 self.xrmmap['xrd1d/q'][:] = row.xrdq
1166 else:
1167 self.xrmmap['xrd1d/q'][:] = row.xrdq[0]
1169 if self.bkgd_xrd1d is not None:
1170 self.xrmmap['xrd1d/counts'][thisrow,] = row.xrd1d - self.bkgd_xrd1d
1171 else:
1172 _ni, _nc, _nq = self.xrmmap['xrd1d/counts'].shape
1173 _rc, _rq = row.xrd1d.shape
1174 _nc = min(_nc, _rc)
1175 _nq = min(_nq, _rq)
1176 self.xrmmap['xrd1d/counts'][thisrow, :_nc, :_nq] = row.xrd1d[:_nc,:_nq]
1178 if self.azwdgs > 1 and row.xrd1d_wdg is not None:
1179 for iwdg,wdggrp in enumerate(self.xrmmap['work/xrdwedge'].values()):
1180 try:
1181 wdggrp['q'] = row.xrdq_wdg[0,:,iwdg]
1182 except:
1183 pass
1185 ## does not yet subtract a background here BECAUSE q-range different
1186 ## for each wedge - won't be same size or shape array
1187 ## mkak 2018.02.26
1188 wdggrp['counts'][thisrow,] = row.xrd1d_wdg[:,:,iwdg]
1191 if self.has_xrd2d and row.xrd2d is not None:
1192 self.xrmmap['xrd2d/counts'][thisrow,] = row.xrd2d
1193 # dt.add("xrd done")
1194 self.last_row = thisrow
1195 self.xrmmap.attrs['Last_Row'] = thisrow
1196 #self.h5root.flush()
1197 # dt.add("flushed h5 file")
1198 # dt.show()
1201 def build_schema(self, npts, nmca=1, nchan=2048, scaler_names=None,
1202 scaler_addrs=None, xrd2d_shape=None, nrows_expected=None,
1203 verbose=False):
1204 '''build schema for detector and scan data'''
1205 self.t0 = time.time()
1206 if not self.check_hostid():
1207 raise GSEXRM_Exception(NOT_OWNER % self.filename)
1208 if not self.write_access:
1209 raise GSEXRM_Exception(READ_ONLY % self.filename)
1211 xrmmap = self.xrmmap
1212 for xkey, xval in xrmmap.attrs.items():
1213 if xkey == 'Version':
1214 self.version = xval
1216 self.xrmmap.attrs['N_Detectors'] = nmca
1217 if scaler_names is None:
1218 scaler_names = []
1219 if scaler_addrs is None:
1220 scaler_addrs = [''] * len(scaler_names)
1222 conf = xrmmap['config']
1223 for key in self.notes:
1224 conf['notes'].attrs[key] = self.notes[key]
1226 if self.npts is None:
1227 self.npts = npts
1229 if self.chunksize is None:
1230 self.chunksize = (1, min(2048, npts), nchan)
1232 if nrows_expected is not None:
1233 NSTART = max(NINIT, nrows_expected)
1234 NSTART = NINIT*2
1236 # positions
1237 pos = xrmmap['positions']
1238 for pname in ('mca realtime', 'mca livetime'):
1239 self.pos_desc.append(pname)
1240 self.pos_addr.append(pname)
1241 npos = len(self.pos_desc)
1242 self.add_data(pos, 'name', strlist(self.pos_desc))
1243 self.add_data(pos, 'address', strlist(self.pos_addr))
1244 pos.create_dataset('pos', (NSTART, npts, npos), np.float32,
1245 maxshape=(None, npts, npos), **self.compress_args)
1247 # cfile = FastMapConfig()
1248 # print(" build_schema -> mapconfig")
1249 # self.add_map_config(cfile.config, nmca=nmca)
1250 conf = xrmmap['config']
1252 offset = conf['mca_calib/offset'][()]
1253 slope = conf['mca_calib/slope'][()]
1254 quad = np.array([0.0]*len(offset))
1255 if 'quad' in conf['mca_calib']:
1256 quad = conf['mca_calib/quad'][()]
1258 if len(offset) != nmca:
1259 raise GSEXRM_Exception("incorrect XRF calibration data: need %d MCAs, not %d" % (nmca, len(offset)))
1262 _ex = np.arange(nchan, dtype=np.float64)
1263 enarr = []
1264 for i in range(len(offset)):
1265 enarr.append(offset[i] + slope[i]*_ex + quad[i]*_ex**2)
1267 if version_ge(self.version, '2.0.0'):
1268 sismap = xrmmap['scalars']
1269 sismap.attrs['type'] = 'scalar detector'
1270 for aname in scaler_names:
1271 sismap.create_dataset(aname, (NSTART, npts), np.float32,
1272 chunks=self.chunksize[:-1],
1273 maxshape=(None, npts), **self.compress_args)
1275 roishape = conf['rois/name'].shape
1276 if roishape[0] > 0:
1277 roi_names = [h5str(s) for s in conf['rois/name']]
1278 roi_limits = np.einsum('jik->ijk', conf['rois/limits'][()])
1279 else:
1280 roi_names = ['_']
1281 roi_limits = np.array([[[0, 2]]])
1283 if verbose:
1284 allmcas = 'with all mcas' if self.all_mcas else 'with only mca sum'
1285 msg = '--- Build XRF Schema: %i ---- MCA: (%i, %i) %s'
1286 print(msg % (npts, nmca, nchan, allmcas))
1288 ## mca1 to mcaN
1289 if self.all_mcas:
1290 for i in range(nmca):
1291 imca = "mca%d" % (i+1)
1292 for grp in (xrmmap, xrmmap['roimap']):
1293 dgrp = grp.create_group(imca)
1294 dgrp.attrs['type'] = 'mca detector'
1295 dgrp.attrs['desc'] = imca
1297 dgrp = xrmmap[imca]
1298 self.add_data(dgrp, 'energy', enarr[i],
1299 attrs={'cal_offset':offset[i],
1300 'cal_slope': slope[i],
1301 'cal_quad': quad[i]})
1302 dgrp.create_dataset('counts', (NSTART, npts, nchan), np.uint32,
1303 chunks=self.chunksize,
1304 maxshape=(None, npts, nchan), **self.compress_args)
1306 for name, dtype in (('realtime', np.int64),
1307 ('livetime', np.int64),
1308 ('dtfactor', np.float32),
1309 ('inpcounts', np.float32),
1310 ('outcounts', np.float32)):
1311 dgrp.create_dataset(name, (NSTART, npts), dtype,
1312 maxshape=(None, npts), **self.compress_args)
1314 dgrp = xrmmap['roimap'][imca]
1315 for rname, rlimit in zip(roi_names,roi_limits[i]):
1316 rgrp = dgrp.create_group(rname)
1317 for aname,dtype in (('raw', np.uint32),
1318 ('cor', np.float32)):
1319 rgrp.create_dataset(aname, (1, npts), dtype,
1320 chunks=self.chunksize[:-1],
1321 maxshape=(None, npts), **self.compress_args)
1323 rlimit = [max(0, rlimit[0]), min(len(enarr[i])-1, rlimit[1])]
1324 lmtgrp = rgrp.create_dataset('limits', data=enarr[i][rlimit])
1325 lmtgrp.attrs['type'] = 'energy'
1326 lmtgrp.attrs['units'] = 'keV'
1328 ## mcasum
1329 for grp in (xrmmap, xrmmap['roimap']):
1330 dgrp = grp.create_group('mcasum')
1331 dgrp.attrs['type'] = 'virtual mca detector'
1332 dgrp.attrs['desc'] = 'sum of detectors'
1334 dgrp = xrmmap['mcasum']
1335 if len(enarr) == 0:
1336 enarr = [np.zeros(nchan)]
1337 counts = [np.zeros(nchan)]
1338 offset = slope = quad = [0.0]
1339 self.add_data(dgrp, 'energy', enarr[0],
1340 attrs={'cal_offset':offset[0],
1341 'cal_slope': slope[0],
1342 'cal_quad': quad[0]})
1343 dgrp.create_dataset('counts', (NSTART, npts, nchan), np.float64,
1344 chunks=self.chunksize,
1345 maxshape=(None, npts, nchan), **self.compress_args)
1347 for name, dtype in (('realtime', np.int64),
1348 ('livetime', np.int64),
1349 ('dtfactor', np.float32),
1350 ('inpcounts', np.float32),
1351 ('outcounts', np.float32)):
1352 dgrp.create_dataset(name, (NSTART, npts), dtype,
1353 maxshape=(None, npts), **self.compress_args)
1356 dgrp = xrmmap['roimap']['mcasum']
1357 for rname, rlimit in zip(roi_names, roi_limits[0]):
1358 rgrp = dgrp.create_group(rname)
1359 for aname,dtype in (('raw', np.uint32),
1360 ('cor', np.float32)):
1361 rgrp.create_dataset(aname, (1, npts), dtype,
1362 chunks=self.chunksize[:-1],
1363 maxshape=(None, npts), **self.compress_args)
1364 rlimit = [max(0, rlimit[0]), min(len(enarr[0])-1, rlimit[1])]
1366 lmtgrp = rgrp.create_dataset('limits', data=enarr[0][rlimit],
1367 **self.compress_args)
1368 lmtgrp.attrs['type'] = 'energy'
1369 lmtgrp.attrs['units'] = 'keV'
1371 if version_ge(self.version, '2.1.0'): # NEW ROI RAW DAT
1372 rdat = xrmmap['roimap']
1373 det_addr = [i.strip() for i in scaler_addrs]
1374 det_desc = [i.strip() for i in scaler_names]
1375 for addr in conf['rois/address']:
1376 det_addr.extend([h5str(addr) % (i+1) for i in range(nmca)])
1378 for nam in roi_names:
1379 det_desc.extend(["%s (mca%i)" % (nam, i+1) for i in range(nmca)])
1381 sums_map = {}
1382 sums_desc = []
1383 nsum = 0
1384 for idet, addr in enumerate(det_desc):
1385 if '(mca' in addr:
1386 addr = addr.split('(mca')[0].strip()
1388 if addr not in sums_map:
1389 sums_map[addr] = []
1390 sums_desc.append(addr)
1391 sums_map[addr].append(idet)
1392 nsum = max([len(s) for s in sums_map.values()])
1393 sums_list = []
1394 for sname in sums_desc:
1395 slist = sums_map[sname]
1396 if len(slist) < nsum:
1397 slist.extend([-1]*(nsum-len(slist)))
1398 sums_list.append(slist)
1400 nsum = len(sums_list)
1401 nsca = len(det_desc)
1402 sums_list = np.array(sums_list)
1404 self.add_data(rdat, 'det_name', strlist(det_desc))
1405 self.add_data(rdat, 'det_address', strlist(det_addr))
1406 self.add_data(rdat, 'sum_name', strlist(sums_desc))
1407 self.add_data(rdat, 'sum_list', sums_list)
1409 for name, nx, dtype in (('det_raw', nsca, np.uint32),
1410 ('det_cor', nsca, np.float32),
1411 ('sum_raw', nsum, np.uint32),
1412 ('sum_cor', nsum, np.float32)
1413 ):
1414 rdat.create_dataset(name, (NSTART, npts, nx), dtype,
1415 chunks=(2, npts, nx),
1416 maxshape=(None, npts, nx), **self.compress_args)
1418 else: # version 1.*
1419 if verbose:
1420 msg = '--- Build XRF Schema: {:d} MCA: ({:d}, {:d})'
1421 print(msg.format(npts, nmca, nchan))
1423 roi_names = [h5str(s) for s in conf['rois/name']]
1424 roi_addrs = [h5str(s) for s in conf['rois/address']]
1425 roi_limits = conf['rois/limits'][()]
1426 for imca in range(nmca):
1427 dname = 'det%i' % (imca+1)
1428 dgrp = xrmmap.create_group(dname)
1429 dgrp.attrs['type'] = 'mca detector'
1430 dgrp.attrs['desc'] = 'mca%i' % (imca+1)
1431 self.add_data(dgrp, 'energy', enarr[imca],
1432 attrs={'cal_offset':offset[imca],
1433 'cal_slope': slope[imca],
1434 'cal_quad': quad[imca]})
1435 self.add_data(dgrp, 'roi_name', strlist(roi_names))
1436 self.add_data(dgrp, 'roi_address', strlist([s % (imca+1) for s in roi_addrs]))
1437 self.add_data(dgrp, 'roi_limits', roi_limits[:,imca,:])
1439 dgrp.create_dataset('counts', (NSTART, npts, nchan), np.uint32,
1440 chunks=self.chunksize,
1441 maxshape=(None, npts, nchan), **self.compress_args)
1442 for name, dtype in (('realtime', np.int64),
1443 ('livetime', np.int64),
1444 ('dtfactor', np.float32),
1445 ('inpcounts', np.float32),
1446 ('outcounts', np.float32)):
1447 dgrp.create_dataset(name, (NSTART, npts), dtype,
1448 maxshape=(None, npts), **self.compress_args)
1450 # add 'virtual detector' for corrected sum:
1451 dgrp = xrmmap.create_group('detsum')
1452 dgrp.attrs['type'] = 'virtual mca'
1453 dgrp.attrs['desc'] = 'deadtime corrected sum of detectors'
1454 self.add_data(dgrp, 'energy', enarr[0],
1455 attrs={'cal_offset':offset[0],
1456 'cal_slope': slope[0],
1457 'cal_quad': quad[0]})
1458 self.add_data(dgrp, 'roi_name', strlist(roi_names))
1459 self.add_data(dgrp, 'roi_address', strlist((s % 1) for s in roi_addrs))
1460 self.add_data(dgrp, 'roi_limits', roi_limits[: ,0, :])
1461 dgrp.create_dataset('counts', (NSTART, npts, nchan), np.uint32,
1462 chunks=self.chunksize,
1463 maxshape=(None, npts, nchan), **self.compress_args)
1464 # roi map data
1465 scan = xrmmap['roimap']
1466 det_addr = [i.strip() for i in scaler_addrs]
1467 det_desc = [i.strip() for i in scaler_names]
1468 for addr in roi_addrs:
1469 det_addr.extend([addr % (i+1) for i in range(nmca)])
1471 for desc in roi_names:
1472 det_desc.extend(["%s (mca%i)" % (desc, i+1)
1473 for i in range(nmca)])
1475 sums_map = {}
1476 sums_desc = []
1477 nsum = 0
1478 for idet, addr in enumerate(det_desc):
1479 if '(mca' in addr:
1480 addr = addr.split('(mca')[0].strip()
1482 if addr not in sums_map:
1483 sums_map[addr] = []
1484 sums_desc.append(addr)
1485 sums_map[addr].append(idet)
1486 nsum = max([len(s) for s in sums_map.values()])
1487 sums_list = []
1488 for sname in sums_desc:
1489 slist = sums_map[sname]
1490 if len(slist) < nsum:
1491 slist.extend([-1]*(nsum-len(slist)))
1492 sums_list.append(slist)
1494 nsum = len(sums_list)
1495 nsca = len(det_desc)
1497 sums_list = np.array(sums_list)
1499 self.add_data(scan, 'det_name', strlist(det_desc))
1500 self.add_data(scan, 'det_address', strlist(det_addr))
1501 self.add_data(scan, 'sum_name', strlist(sum_desc))
1502 self.add_data(scan, 'sum_list', sums_list)
1504 nxx = min(nsca, 8)
1505 for name, nx, dtype in (('det_raw', nsca, np.uint32),
1506 ('det_cor', nsca, np.float32),
1507 ('sum_raw', nsum, np.uint32),
1508 ('sum_cor', nsum, np.float32)):
1509 scan.create_dataset(name, (NSTART, npts, nx), dtype,
1510 chunks=(2, npts, nx),
1511 maxshape=(None, npts, nx), **self.compress_args)
1514 if self.has_xrd2d or self.has_xrd1d:
1515 if self.has_xrd2d:
1516 xrdpts, xpixx, xpixy = xrd2d_shape # row.xrd2d.shape
1517 if verbose:
1518 prtxt = '--- Build XRD Schema: %i ---- 2D XRD: (%i, %i)'
1519 print(prtxt % (npts, xpixx, xpixy))
1521 xrdgrp = ensure_subgroup('xrd2d', xrmmap, dtype='2DXRD')
1523 xrdgrp.attrs['type'] = 'xrd2d detector'
1524 xrdgrp.attrs['desc'] = '' #'add detector name eventually'
1526 xrdgrp.create_dataset('mask', (xpixx, xpixy), np.uint16, **self.compress_args)
1527 xrdgrp.create_dataset('background', (xpixx, xpixy), np.uint16, **self.compress_args)
1529 chunksize_2DXRD = (1, npts, xpixx, xpixy)
1530 xrdgrp.create_dataset('counts', (NSTART, npts, xpixx, xpixy), np.uint32,
1531 chunks = chunksize_2DXRD,
1532 maxshape=(None, npts, xpixx, xpixy), **self.compress_args)
1534 if self.has_xrd1d:
1535 xrdgrp = ensure_subgroup('xrd1d', xrmmap)
1536 xrdgrp.attrs['type'] = 'xrd1d detector'
1537 xrdgrp.attrs['desc'] = 'pyFAI calculation from xrd2d data'
1539 xrdgrp.create_dataset('q', (self.qstps,), np.float32, **self.compress_args)
1540 xrdgrp.create_dataset('background', (self.qstps,), np.float32, **self.compress_args)
1542 chunksize_xrd1d = (1, npts, self.qstps)
1543 xrdgrp.create_dataset('counts',
1544 (NSTART, npts, self.qstps),
1545 np.float32,
1546 chunks = chunksize_xrd1d,
1547 maxshape=(None, npts, self.qstps), **self.compress_args)
1549 if self.azwdgs > 1:
1550 xrmmap['work'].create_group('xrdwedge')
1551 for azi in range(self.azwdgs):
1552 wdggrp = xrmmap['work/xrdwedge'].create_group('wedge_%02d' % azi)
1554 wdggrp.create_dataset('q', (self.qstps,), np.float32, **self.compress_args)
1556 wdggrp.create_dataset('counts',
1557 (NSTART, npts, self.qstps),
1558 np.float32,
1559 chunks = chunksize_xrd1d,
1560 maxshape=(None, npts, self.qstps), **self.compress_args)
1562 #wdggrp.create_dataset('limits', (2,), np.float32)
1563 wdg_sz = 360./self.azwdgs
1564 wdg_lmts = np.array([azi*wdg_sz, (azi+1)*wdg_sz]) - 180
1565 wdggrp.create_dataset('limits', data=wdg_lmts)
1567 try:
1568 print('\nStart: %s' % isotime(self.starttime))
1569 except:
1570 pass
1572 self.h5root.flush()
1574 def add_xrd1d(self, qstps=None):
1575 xrd1dgrp = ensure_subgroup('xrd1d',self.xrmmap)
1576 xrdcalfile = bytes2str(xrd1dgrp.attrs.get('calfile', ''))
1577 if os.path.exists(xrdcalfile):
1578 print('Using calibration file : %s' % xrdcalfile)
1579 try:
1580 nrows, npts , xpixx, xpixy = self.xrmmap['xrd2d/counts'].shape
1581 except:
1582 return
1584 if qstps is not None: self.qstps = qstps
1586 pform ='\n--- Build XRD1D Schema (%i, %i, %i) from 2D XRD (%i, %i, %i, %i) ---'
1587 print(pform % (nrows, npts, self.qstps, nrows, npts, xpixx, xpixy))
1589 try:
1590 xrd1dgrp.attrs['type'] = 'xrd1d detector'
1591 xrd1dgrp.attrs['desc'] = 'pyFAI calculation from xrd2d data'
1592 xrd1dgrp.create_dataset('q', (self.qstps,), np.float32)
1593 xrd1dgrp.create_dataset('background', (self.qstps,), np.float32)
1595 chunksize_xrd1d = (1, npts, self.qstps)
1596 xrd1dgrp.create_dataset('counts',
1597 (nrows, npts, self.qstps),
1598 np.float32,
1599 chunks = chunksize_xrd1d)
1601 attrs = {'steps':self.qstps,'mask':self.xrd2dmaskfile,'flip':self.flip}
1602 print('\nStart: %s' % isotime())
1603 for i in np.arange(nrows):
1604 rowq, row1d = integrate_xrd_row(self.xrmmap['xrd2d/counts'][i],xrdcalfile,**attrs)
1605 if i == 0:
1606 self.xrmmap['xrd1d/q'][:] = rowq[0]
1607 self.xrmmap['xrd1d/counts'][i,] = row1d
1609 self.has_xrd1d = True
1610 # print('End: %s' % isotime())
1611 except:
1612 print('xrd1d data already in file.')
1613 return
1615 def get_slice_y(self):
1616 for name, val in zip([h5str(a) for a in self.xrmmap['config/environ/name']],
1617 [h5str(a) for a in self.xrmmap['config/environ/value']]):
1618 name = str(name).lower()
1619 if name.startswith('sample'):
1620 name = name.replace('samplestage.', '')
1621 if name.lower() == 'fine y' or name.lower() == 'finey':
1622 return float(val)
1624 def get_datapath_list(self, remove='raw'):
1625 def find_detector(group):
1626 sub_list = []
1627 if 'counts' in group.keys():
1628 sub_list += [group['counts'].name]
1629 elif 'scal' in group.name:
1630 for key, val in dict(group).items():
1631 sub_list += [group[key].name]
1632 return sub_list
1634 dlist = []
1635 for det in self.get_detector_list():
1636 for idet in find_detector(self.xrmmap[det]):
1637 if not (remove in idet):
1638 dlist.append(idet)
1640 return dlist
1642 def get_roi_list(self, det_name, force=False):
1643 """
1644 get a list of rois from detector
1645 """
1646 detname = self.get_detname(det_name)
1647 if not force and (detname not in EXTRA_DETGROUPS):
1648 roilist = self.roi_names.get(detname, None)
1649 if roilist is not None:
1650 return roilist
1652 roigrp = ensure_subgroup('roimap', self.xrmmap, dtype='roi maps')
1653 def sort_roi_limits(roidetgrp):
1654 roi_name, roi_limits = [],[]
1655 for name in roidetgrp.keys():
1656 roi_name.append(name)
1657 roi_limits.append(roidetgrp[name]['limits'][0])
1658 return [y for (x,y) in sorted(zip(roi_limits,roi_name))]
1660 def sort_orderby(det, orderby):
1661 roi_name, roi_order = [], []
1662 for name in det.keys():
1663 roi_name.append(name)
1664 roi_order.append(det[name].attrs.get(orderby, 9999))
1665 return [y for (x, y) in sorted(zip(roi_order, roi_name))]
1667 rois = []
1668 if version_ge(self.version, '2.0.0'):
1669 if detname in roigrp.keys():
1670 rois = sort_roi_limits(roigrp[detname])
1671 else:
1672 det = self.xrmmap[detname]
1673 if (detname in EXTRA_DETGROUPS or
1674 'detector' in det.attrs.get('type')):
1675 rois = list(det.keys())
1676 orderby = det.attrs.get('orderby', None)
1677 if orderby is not None:
1678 rois = sort_orderby(det, orderby)
1679 else:
1680 if detname in EXTRA_DETGROUPS:
1681 rois = list(self.xrmmap[detname].keys())
1682 elif detname in self.xrmmap.keys():
1683 rois = list(roigrp['sum_name']) + rois
1684 try:
1685 rois = sort_roi_limits(roigrp[detname]) + rois
1686 except:
1687 pass
1688 rois.append('1')
1689 self.roi_names[detname] = [h5str(a) for a in rois]
1690 return self.roi_names[detname]
1692 def get_detector_list(self, use_cache=False):
1693 """get a list of detector groups,
1694 ['mcasum', 'mca1', ..., 'scalars', 'work', 'xrd1d', ...]
1695 """
1696 workgroup = ensure_subgroup('work', self.xrmmap)
1697 if use_cache and self.detector_list is not None:
1698 return self.detector_list
1699 def build_dlist(group):
1700 detlist, sumslist = [], []
1701 for key, grp in group.items():
1702 if ('det' in bytes2str(grp.attrs.get('type', '')) or
1703 'mca' in bytes2str(grp.attrs.get('type', ''))):
1704 if 'sum' in key.lower():
1705 sumslist.append(key)
1706 else:
1707 detlist.append(key)
1708 return sumslist + detlist
1710 xrmmap = self.xrmmap
1711 det_list = []
1712 if version_ge(self.version, '2.0.0'):
1713 det_list = build_dlist(xrmmap['roimap'])
1714 else:
1715 det_list = build_dlist(xrmmap)
1716 for det in build_dlist(xrmmap['roimap']):
1717 if det not in det_list:
1718 det_list.append(det)
1720 # add any other groups with 'detector' in the `type` attribute:
1721 for det, grp in xrmmap.items():
1722 attrs = getattr(grp, 'attrs', {'type': ''})
1723 if det not in det_list and 'detector' in h5str(attrs.get('type', '')):
1724 det_list.append(det)
1725 self.detector_list = det_list
1726 if len(det_list) < 1:
1727 det_list = ['']
1728 self.detector_list = None
1729 return det_list
1731 def reset_flags(self):
1732 '''
1733 Reads hdf5 file for data and sets the flags.
1734 mkak 2016.08.30 // rewritten mkak 2017.08.03 // rewritten mkak 2017.12.05
1735 '''
1736 for det in self.get_detector_list():
1737 if det in self.xrmmap:
1739 detgrp = self.xrmmap[det]
1741 dettype = bytes2str(detgrp.attrs.get('type', '')).lower()
1742 if 'mca' in dettype:
1743 self.has_xrf = 'counts' in detgrp
1744 elif 'xrd2d' in dettype:
1745 self.has_xrd2d = 'counts' in detgrp
1746 elif 'xrd1d' in dettype:
1747 self.has_xrd1d = 'counts' in detgrp
1748 elif det == 'xrd': ## compatible with old version
1749 try:
1750 detgrp['data1d']
1751 self.has_xrd1d = True
1752 except:
1753 pass
1754 try:
1755 detgrp['data2D']
1756 self.has_xrd2d = True
1757 except:
1758 pass
1760 def print_flags(self):
1761 print(' HAS XRF, XRD1D, XRD2D: %s, %s, %s' % (self.has_xrf,
1762 self.has_xrd1d,
1763 self.has_xrd2d))
1765 def resize_arrays(self, nrow, force_shrink=True):
1766 "resize all arrays for new nrow size"
1767 if not self.check_hostid():
1768 raise GSEXRM_Exception(NOT_OWNER % self.filename)
1769 if not self.write_access:
1770 raise GSEXRM_Exception(READ_ONLY % self.filename)
1771 if version_ge(self.version, '2.0.0'):
1773 g = self.xrmmap['positions/pos']
1774 old, npts, nx = g.shape
1775 if nrow < old and not force_shrink:
1776 return
1777 g.resize((nrow, npts, nx))
1779 for g in self.xrmmap.values():
1780 type_attr = bytes2str(g.attrs.get('type', ''))
1781 if type_attr.find('det') > -1:
1782 if type_attr.startswith('scalar'):
1783 for aname in g.keys():
1784 oldnrow, npts = g[aname].shape
1785 g[aname].resize((nrow, npts))
1786 elif type_attr.startswith('mca'):
1787 oldnrow, npts, nchan = g['counts'].shape
1788 g['counts'].resize((nrow, npts, nchan))
1789 for aname in ('livetime', 'realtime',
1790 'inpcounts', 'outcounts', 'dtfactor'):
1791 g[aname].resize((nrow, npts))
1792 elif type_attr.startswith('virtual mca'):
1793 oldnrow, npts, nchan = g['counts'].shape
1794 g['counts'].resize((nrow, npts, nchan))
1795 for aname in ('livetime', 'realtime',
1796 'inpcounts', 'outcounts', 'dtfactor'):
1797 if aname in g:
1798 g[aname].resize((nrow, npts))
1800 elif type_attr.startswith('xrd2d'):
1801 oldnrow, npts, xpixx, xpixy = g['counts'].shape
1802 g['counts'].resize((nrow, npts, xpixx, xpixy))
1803 elif type_attr.startswith('xrd1d'):
1804 oldnrow, npts, qstps = g['counts'].shape
1805 g['counts'].resize((nrow, npts, qstps))
1807 if self.azwdgs > 1:
1808 for g in self.xrmmap['work/xrdwedge'].values():
1809 g['counts'].resize((nrow, npts, qstps))
1811 if version_ge(self.version, '2.1.0'):
1812 for bname in ('det_raw', 'det_cor', 'sum_raw', 'sum_cor'):
1813 g = self.xrmmap['roimap'][bname]
1814 old, npts, nx = g.shape
1815 g.resize((nrow, npts, nx))
1816 else:
1817 for g in self.xrmmap['roimap'].values(): # loop through detectors in roimap
1818 for h in g.values(): # loop through rois in roimap
1819 for aname in ('raw','cor'):
1820 oldnrow, npts = h[aname].shape
1821 h[aname].resize((nrow, npts))
1823 else: ## old file format method
1825 realmca_groups = []
1826 virtmca_groups = []
1827 for g in self.xrmmap.values():
1828 # include both real and virtual mca detectors!
1829 type_attr = bytes2str(g.attrs.get('type', ''))
1830 if type_attr.find('det') > -1 or type_attr.find('mca') > -1:
1831 if type_attr.startswith('mca'):
1832 realmca_groups.append(g)
1833 elif type_attr.startswith('virtual mca'):
1834 virtmca_groups.append(g)
1835 elif type_attr.startswith('xrd2d'):
1836 oldnrow, npts, xpixx, xpixy = g['counts'].shape
1837 g['counts'].resize((nrow, npts, xpixx, xpixy))
1838 elif type_attr.startswith('xrd1d'):
1839 oldnrow, npts, qstps = g['counts'].shape
1840 g['counts'].resize((nrow, npts, qstps))
1842 if self.azwdgs > 1:
1843 for g in self.xrmmap['work/xrdwedge'].values():
1844 g['counts'].resize((nrow, npts, qstps))
1846 oldnrow, npts, nchan = realmca_groups[0]['counts'].shape
1847 for g in realmca_groups:
1848 g['counts'].resize((nrow, npts, nchan))
1849 for aname in ('livetime', 'realtime',
1850 'inpcounts', 'outcounts', 'dtfactor'):
1851 g[aname].resize((nrow, npts))
1853 for g in virtmca_groups:
1854 g['counts'].resize((nrow, npts, nchan))
1856 g = self.xrmmap['positions/pos']
1857 old, npts, nx = g.shape
1858 g.resize((nrow, npts, nx))
1860 for bname in ('det_raw', 'det_cor', 'sum_raw', 'sum_cor'):
1861 g = self.xrmmap['roimap'][bname]
1862 old, npts, nx = g.shape
1863 g.resize((nrow, npts, nx))
1865 self.h5root.flush()
1867 def add_work_array(self, data, name, parent='work',
1868 dtype='virtual detector', **kws):
1869 '''
1870 add an array to the work group of processed arrays
1871 '''
1872 workgroup = ensure_subgroup(parent, self.xrmmap, dtype=dtype)
1873 if name is None:
1874 name = 'array_%3.3i' % (1+len(workgroup))
1875 if name in workgroup:
1876 raise ValueError("array name '%s' exists in '%s" % (name, parent))
1877 ds = workgroup.create_dataset(name, data=data)
1878 for key, val in kws.items():
1879 ds.attrs[key] = val
1880 self.h5root.flush()
1882 def add_work_arrays(self, arraydict, parent='work'):
1883 if parent in self.xrmmap:
1884 del self.xrmmap[parent]
1885 self.h5root.flush()
1887 workgroup = ensure_subgroup(parent, self.xrmmap)
1888 workgroup.attrs['orderby'] = 'order'
1889 count = 0
1890 for key, data in arraydict.items():
1891 count += 1
1892 name = fix_varname(key)
1893 if name in workgroup:
1894 del workgroup[name]
1895 ds = workgroup.create_dataset(name, data=data)
1896 ds.attrs['order'] = count
1897 self.h5root.flush()
1899 def del_work_array(self, name, parent='work'):
1900 '''
1901 delete an array to the work group of processed arrays
1902 '''
1903 workgroup = ensure_subgroup(parent, self.xrmmap)
1904 name = h5str(name)
1905 if name in workgroup:
1906 del workgroup[name]
1907 self.h5root.flush()
1909 def get_work_array(self, name, parent='work'):
1910 '''
1911 get an array from the work group of processed arrays by index or name
1912 '''
1913 workgroup = ensure_subgroup(parent, self.xrmmap)
1914 dat = None
1915 name = h5str(name)
1916 if name in workgroup:
1917 dat = workgroup[name]
1918 return dat
1920 def work_array_names(self, parent='work'):
1921 '''
1922 return list of work array descriptions
1923 '''
1924 workgroup = ensure_subgroup(parent, self.xrmmap)
1925 return [h5str(g) for g in workgroup.keys()]
1927 def add_area(self, amask, name=None, desc=None):
1928 '''add a selected area, with optional name
1929 the area is encoded as a boolean array the same size as the map
1931 '''
1932 if not self.check_hostid():
1933 raise GSEXRM_Exception(NOT_OWNER % self.filename)
1934 if not self.write_access:
1935 raise GSEXRM_Exception(READ_ONLY % self.filename)
1937 area_grp = ensure_subgroup('areas', self.xrmmap, dtype='areas')
1938 if name is None:
1939 name = 'area_001'
1940 if len(area_grp) > 0:
1941 count = len(area_grp)
1942 while name in area_grp and count < 9999:
1943 name = 'area_%3.3i' % (count)
1944 count += 1
1945 ds = area_grp.create_dataset(name, data=amask)
1946 if desc is None:
1947 desc = name
1948 ds.attrs['description'] = desc
1949 self.h5root.flush()
1950 return name
1952 def export_areas(self, filename=None):
1953 '''export areas to datafile '''
1954 if filename is None:
1955 file_str = '%s_Areas.npz'
1956 filename = file_str % self.filename
1958 areas = ensure_subgroup('areas', self.xrmmap, dtype='areas')
1959 kwargs = {key: val[:] for key, val in areas.items()}
1960 np.savez(filename, **kwargs)
1961 return filename
1963 def import_areas(self, filename, overwrite=False):
1964 '''import areas from datafile exported by export_areas()'''
1965 fname = Path(filename).name
1966 if fname.endswith('.h5_Areas.npz'):
1967 fname = fname.replace('.h5_Areas.npz', '')
1969 areas = ensure_subgroup('areas', self.xrmmap, dtype='areas')
1971 npzdat = np.load(filename)
1972 for aname in npzdat.files:
1973 desc = name = aname
1974 if name in areas and not overwrite:
1975 name = '%s_%s' % (name, fname)
1976 desc = '%s from %s' % (aname, fname)
1977 self.add_area(npzdat[aname], name=name, desc=desc)
1979 def get_area(self, name=None):
1980 '''
1981 get area group by name or description
1982 '''
1983 area_grp = ensure_subgroup('areas', self.xrmmap, dtype='areas')
1984 if name is not None and name in area_grp:
1985 return area_grp[name]
1986 else:
1987 for aname in area_grp:
1988 if name == bytes2str(area_grp[aname].attrs.get('description','')):
1989 return area_grp[aname]
1990 return None
1992 def get_area_stats(self, name=None, desc=None):
1993 '''return statistics for all raw detector counts/sec values
1995 for each raw detector returns
1996 name, length, mean, standard_deviation,
1997 median, mode, minimum, maximum,
1998 gmean, hmean, skew, kurtosis
2000 '''
2001 area = self.get_area(name=name, desc=desc)
2002 if area is None:
2003 return None
2005 if 'roistats' in area.attrs:
2006 return json.loads(area.attrs.get('roistats',''))
2008 amask = area[()]
2010 roidata = []
2011 d_addrs = [d.lower() for d in self.xrmmap['roimap/det_address']]
2012 d_names = [d for d in self.xrmmap['roimap/det_name']]
2013 # count times
2014 ctime = [1.e-6*self.xrmmap['roimap/det_raw'][:,:,0][amask]]
2015 for i in range(self.xrmmap.attrs.get('N_Detectors',0)):
2016 tname = 'det%i/realtime' % (i+1)
2017 ctime.append(1.e-6*self.xrmmap[tname][()][amask])
2019 for idet, dname in enumerate(d_names):
2020 daddr = d_addrs[idet]
2021 det = 0
2022 if 'mca' in daddr:
2023 det = 1
2024 words = daddr.split('mca')
2025 if len(words) > 1:
2026 det = int(words[1].split('.')[0])
2027 if idet == 0:
2028 d = ctime[0]
2029 else:
2030 d = self.xrmmap['roimap/det_raw'][:,:,idet][amask]/ctime[det]
2032 try:
2033 hmean, gmean = stats.gmean(d), stats.hmean(d)
2034 skew, kurtosis = stats.skew(d), stats.kurtosis(d)
2035 except ValueError:
2036 hmean, gmean, skew, kurtosis = 0, 0, 0, 0
2037 mode = stats.mode(d)
2038 roidata.append((dname, len(d), d.mean(), d.std(), np.median(d),
2039 stats.mode(d), d.min(), d.max(),
2040 gmean, hmean, skew, kurtosis))
2042 if 'roistats' not in area.attrs:
2043 area.attrs['roistats'] = json.dumps(roidata)
2044 self.h5root.flush()
2046 return roidata
2048 def get_translation_axis(self, hotcols=None):
2049 if hotcols is None:
2050 hotcols = self.hotcols
2051 posnames = [bytes2str(n.lower()) for n in self.xrmmap['positions/name']]
2052 if 'x' in posnames:
2053 x = self.get_pos('x', mean=True)
2054 elif 'fine x' in posnames:
2055 x = self.get_pos('fine x', mean=True)
2056 else:
2057 x = self.get_pos(0, mean=True)
2059 if hotcols and x is not None:
2060 if len(x) == self.xrmmap[self.get_detname()]['counts'].shape[1]:
2061 x = x[1:-1]
2063 return x
2065 def get_rotation_axis(self, axis=None, hotcols=None):
2066 if hotcols is None:
2067 hotcols = self.hotcols
2068 posnames = [bytes2str(n.lower()) for n in self.xrmmap['positions/name']]
2069 omega = None
2070 if axis is not None:
2071 if axis in posnames or type(axis) == int:
2072 omega = self.get_pos(axis, mean=True)
2073 else:
2074 omega = None
2075 for pname in posnames:
2076 for aname in ('theta', 'phi', 'omega', 'chi'):
2077 if aname in pname:
2078 omega = self.get_pos(pname, mean=True)
2079 # print(f"using positoner '{pname:s}' as rotation axis ")
2080 break
2081 if omega is not None:
2082 break
2083 if hotcols and omega is not None:
2084 if len(omega) == self.xrmmap[self.get_detname()]['counts'].shape[1]:
2085 omega = omega[1:-1]
2086 return omega
2088 def get_tomography_center(self):
2089 tomogrp = ensure_subgroup('tomo', self.xrmmap)
2090 try:
2091 return tomogrp['center'][()]
2092 except:
2093 self.set_tomography_center()
2095 return tomogrp['center'][()]
2097 def set_tomography_center(self,center=None):
2098 if center is None:
2099 center = len(self.get_translation_axis())/2.
2101 tomogrp = ensure_subgroup('tomo', self.xrmmap)
2102 try:
2103 del tomogrp['center']
2104 except:
2105 pass
2106 tomogrp.create_dataset('center', data=center)
2108 self.h5root.flush()
2110 def get_sinogram(self, roi_name, det=None, trim_sino=False,
2111 hotcols=None, dtcorrect=None, **kws):
2112 '''extract roi map for a pre-defined roi by name
2114 Parameters
2115 ---------
2116 roiname : str ROI name
2117 det : str detector name
2118 dtcorrect : None or bool [None] deadtime correction
2119 hotcols : None or bool [None] suppress hot columns
2121 Returns
2122 -------
2123 sinogram for ROI data
2124 sinogram_order (needed for knowing shape of sinogram)
2126 Notes
2127 -----
2128 if dtcorrect or hotcols is None, they are taken from
2129 self.dtcorrect and self.hotcols
2130 '''
2131 if hotcols is None:
2132 hotcols = self.hotcols
2133 if dtcorrect is None:
2134 dtcorrect = self.dtcorrect
2136 sino = self.get_roimap(roi_name, det=det, hotcols=hotcols, **kws)
2137 x = self.get_translation_axis(hotcols=hotcols)
2138 omega = self.get_rotation_axis(hotcols=hotcols)
2140 if omega is None:
2141 print('\n** Cannot compute tomography: no rotation axis specified in map. **')
2142 return
2144 if trim_sino:
2145 sino, x, omega = trim_sinogram(sino, x, omega)
2147 return reshape_sinogram(sino, x, omega)
2149 def get_tomograph(self, sino, omega=None, center=None, hotcols=None, **kws):
2150 '''
2151 returns tomo_center, tomo
2152 '''
2153 if hotcols is None:
2154 hotcols = self.hotcols
2155 if center is None:
2156 center = self.get_tomography_center()
2157 if omega is None:
2158 omega = self.get_rotation_axis(hotcols=hotcols)
2159 if omega is None:
2160 print('\n** Cannot compute tomography: no rotation axis specified in map. **')
2161 return
2163 center, recon = tomo_reconstruction(sino, omega, center=center, **kws)
2164 self.set_tomography_center(center=center)
2165 return recon
2167 def save_tomograph(self, datapath, algorithm='gridrec',
2168 filter_name='shepp', num_iter=1, dtcorrect=None,
2169 hotcols=None, **kws):
2170 '''
2171 saves group for tomograph for selected detector
2172 '''
2173 if hotcols is None:
2174 hotcols = self.hotcols
2175 if dtcorrect is None:
2176 dtcorrect = self.dtcorrect
2178 ## check to make sure the selected detector exists for reconstructions
2179 detlist = self.get_datapath_list(remove=None)
2180 if datapath not in detlist:
2181 print("Detector '%s' not found in data." % datapath)
2182 print('Known detectors: %s' % detlist)
2183 return
2184 datagroup = self.xrmmap[datapath]
2186 ## check to make sure there is data to perform tomographic reconstruction
2187 center = self.get_tomography_center()
2189 x = self.get_translation_axis(hotcols=hotcols)
2190 omega = self.get_rotation_axis(hotcols=hotcols)
2192 if omega is None:
2193 print('\n** Cannot compute tomography: no rotation axis specified in map. **')
2194 return
2196 ## define detector path
2197 detgroup = datagroup
2198 while isinstance(detgroup,h5py.Dataset):
2199 detgroup = detgroup.parent
2200 detpath = detgroup.name
2202 ## create path for saving data
2203 tpath = datapath.replace('/xrmmap','/tomo')
2204 tpath = tpath.replace('/scalars','')
2205 if tpath.endswith('raw'):
2206 tpath = tpath.replace('_raw','')
2207 dtcorrect = False
2208 elif tpath.endswith('counts'):
2209 tpath = Path(tpath).absolute().parent.as_posix()
2211 ## build path for saving data in tomo-group
2212 grp = self.xrmmap
2213 for kpath in tpath.split('/'):
2214 if len(kpath) > 0:
2215 grp = ensure_subgroup(kpath ,grp)
2216 tomogrp = grp
2218 ## define sino group from datapath
2219 if 'scalars' in datapath or 'xrd' in datapath:
2220 sino = datagroup[()]
2221 elif dtcorrect:
2222 if 'sum' in datapath:
2223 sino = np.zeros(np.shape(np.einsum('jki->ijk', datagroup[()])))
2224 for i in range(self.nmca):
2225 idatapath = datapath.replace('sum', str(i+1))
2226 idatagroup = self.xrmmap[idatapath]
2227 idetpath = detpath.replace('sum', str(i+1))
2228 idetgroup = self.xrmmap[idetpath]
2229 sino += np.einsum('jki->ijk', idatagroup[()]) * idetgroup['dtfactor'][()]
2231 else:
2232 sino = np.einsum('jki->ijk', datagroup[()]) * detgroup['dtfactor'][()]
2233 else:
2234 sino = datagroup[()]
2236 sino,order = reshape_sinogram(sino, x, omega)
2238 center, tomo = tomo_reconstruction(sino, algorithm=algorithm,
2239 filter_name=filter_name,
2240 num_iter=num_iter, omega=omega,
2241 center=center, sinogram_order=order)
2243 tomogrp.attrs['tomo_alg'] = '-'.join([str(t) for t in (algorithm, filter_name)])
2244 tomogrp.attrs['center'] = '%0.2f pixels' % (center)
2246 try:
2247 tomogrp.create_dataset('counts', data=np.swapaxes(tomo,0,2), **self.compress_args)
2248 except:
2249 del tomogrp['counts']
2250 tomogrp.create_dataset('counts', data=np.swapaxes(tomo,0,2), **self.compress_args)
2252 for data_tag in ('energy','q'):
2253 if data_tag in detgroup.keys():
2254 try:
2255 tomogrp.create_dataset(data_tag, data=detgroup[data_tag])
2256 del tomogrp[data_tag]
2257 except:
2258 del tomogrp[data_tag]
2259 tomogrp.create_dataset(data_tag, data=detgroup[data_tag])
2261 for key, val in dict(detgroup.attrs).items():
2262 tomogrp.attrs[key] = val
2264 self.h5root.flush()
2266 def take_ownership(self):
2267 "claim ownership of file"
2268 if self.xrmmap is None or not self.write_access:
2269 return
2270 self.xrmmap.attrs['Process_Machine'] = get_machineid()
2271 self.xrmmap.attrs['Process_ID'] = os.getpid()
2272 self.h5root.flush()
2274 def release_ownership(self):
2275 if not self.write_access:
2276 return
2277 self.xrmmap.attrs['Process_Machine'] = ''
2278 self.xrmmap.attrs['Process_ID'] = 0
2279 self.xrmmap.attrs['Last_Row'] = self.last_row
2281 def check_ownership(self, take_ownership=True):
2282 return self.check_hostid(take_ownership=take_ownership)
2284 def check_hostid(self, take_ownership=True):
2285 '''checks host and id of file:
2286 returns True if this process the owner of the file
2288 By default, this takes ownership if it can.
2289 '''
2290 if self.xrmmap is None:
2291 return False
2292 if not self.write_access:
2293 return True
2295 attrs = self.xrmmap.attrs
2296 self.folder = attrs['Map_Folder']
2297 file_mach = attrs['Process_Machine']
2298 file_pid = attrs['Process_ID']
2299 if len(file_mach) < 1 or file_pid < 1:
2300 if take_ownership:
2301 self.take_ownership()
2302 return True
2303 return (file_mach == get_machineid() and file_pid == os.getpid())
2305 def folder_has_newdata(self):
2306 if self.folder is not None and isGSEXRM_MapFolder(self.folder):
2307 self.read_master()
2308 return (self.last_row < len(self.rowdata)-1)
2309 return False
2311 def read_master(self):
2312 "reads master file for toplevel scan info"
2313 if self.folder is None or not isGSEXRM_MapFolder(self.folder):
2314 return
2315 self.masterfile = unixpath(Path(self.folder, self.MasterFile))
2316 header, rows, mtime = [], [], -1
2317 if self.scandb is not None:
2318 # check that this map folder is the one currently running from scandb:
2319 try:
2320 db_folder = toppath(self.scandb.get_info('map_folder'))
2321 except:
2322 db_folder = None
2323 disk_folder = toppath(Path(self.folder).absolute().as_posix())
2325 if db_folder == disk_folder: # this is the current map
2326 mastertext = self.scandb.get_slewscanstatus()
2327 mtime = time.time()
2328 header, rows = [], []
2329 for srow in mastertext:
2330 line = str(srow.text.strip())
2331 if line.startswith('#'):
2332 header.append(line)
2333 else:
2334 rows.append(line.split())
2336 if len(header) < 1 or mtime < 0: # this is *not* the map that is currently being collected:
2337 # if file the master file is not new, the current row data is OK:
2339 try:
2340 header, rows = readMasterFile(self.masterfile)
2341 except IOError:
2342 raise GSEXRM_Exception("cannot read Master file from '%s'" %
2343 self.masterfile)
2345 mtime = os.stat(self.masterfile).st_mtime
2346 if mtime < (self.master_modtime+1.0) and len(self.rowdata) > 1:
2347 # print("READ MASTER not a new masterfile ", len(self.rowdata), len(rows))
2348 return len(self.rowdata)
2350 self.master_modtime = mtime
2352 self.notes['end_time'] = isotime(mtime)
2353 self.master_header = header
2354 # carefully read rows to avoid repeated rows due to bad collection
2355 self.rowdata = []
2356 _yl, _xl, _s1 = None, None, None
2357 for row in rows:
2358 yval, xrff, sisf = row[0], row[1], row[2]
2359 il = len(self.rowdata)-1
2360 if il > -1:
2361 _yl, _xl, _s1 = (self.rowdata[il][0],
2362 self.rowdata[il][1],
2363 self.rowdata[il][2])
2364 # skip repeated rows in master file
2365 if yval != _yl and (xrff != _xl or sisf != _s1):
2366 self.rowdata.append(row)
2368 self.scan_version = 1.00
2369 self.nrows_expected = None
2370 self.start_time = time.ctime()
2371 for line in header:
2372 words = line.split('=')
2373 if 'scan.starttime' in words[0].lower():
2374 self.start_time = words[1].strip()
2375 self.notes['scan_start_time'] = self.start_time
2376 elif 'scan.version' in words[0].lower():
2377 self.scan_version = words[1].strip()
2378 elif 'scan.nrows_expected' in words[0].lower():
2379 self.nrows_expected = int(words[1].strip())
2380 self.scan_version = float(self.scan_version)
2382 # print("read_master scan version = ", self.scan_version)
2383 self.stop_time = time.ctime(self.master_modtime)
2384 try:
2385 last_file = Path(self.folder, rows[-1][2])
2386 self.stop_time = time.ctime(os.stat(last_file).st_ctime)
2387 except:
2388 pass
2389 self.notes['scan_end_time'] = self.stop_time
2390 self.notes['scan_version'] = self.scan_version
2392 if self.scan_version < 1.35 and (self.has_xrd2d or self.has_xrd1d):
2393 xrd_files = [fn for fn in os.listdir(self.folder) if fn.endswith('nc')]
2394 for i, addxrd in enumerate(xrd_files):
2395 self.rowdata[i].insert(4, addxrd)
2397 cfile = FastMapConfig()
2398 cfile.Read(Path(self.folder, self.ScanFile))
2400 mapconf = self.mapconf = cfile.config
2402 if self.filename is None:
2403 self.filename = mapconf['scan']['filename']
2404 if not self.filename.endswith('.h5'):
2405 self.filename = "%s.h5" % self.filename
2407 slow_pos = mapconf['slow_positioners']
2408 fast_pos = mapconf['fast_positioners']
2410 scanconf = mapconf['scan']
2411 self.dimension = scanconf['dimension']
2412 start = mapconf['scan']['start1']
2413 stop = mapconf['scan']['stop1']
2414 step = mapconf['scan']['step1']
2415 span = abs(stop-start)
2416 self.npts = int(abs(abs(step)*1.01 + span)/abs(step))
2418 pos1 = scanconf['pos1']
2419 self.pos_addr = [pos1]
2420 self.pos_desc = [slow_pos[pos1]]
2421 # note: XPS gathering file now saving ONLY data for the fast axis
2423 if self.dimension > 1:
2424 yaddr = scanconf['pos2']
2425 self.pos_addr.append(yaddr)
2426 self.pos_desc.append(slow_pos[yaddr])
2428 return len(rows)
2430 def get_detname(self, det=None):
2431 "return XRMMAP group for a detector"
2433 mcastr = 'mca' if version_ge(self.version, '2.0.0') else 'det'
2434 detname = '%ssum' % mcastr
2436 if isinstance(det, str):
2437 for d in self.get_detector_list():
2438 if det.lower() == d.lower():
2439 detname = d
2440 elif isinstance(det, int):
2441 if det in range(1, self.nmca+1):
2442 detname = '%s%i' % (mcastr, det)
2444 return detname
2446 def get_detgroup(self, det=None):
2447 "return XRMMAP group for a detector"
2448 return self.xrmmap[self.get_detname(det)]
2450 def get_energy(self, det=None):
2451 '''return energy array for a detector'''
2452 try:
2453 group = self.xrmmap[det]
2454 except:
2455 group = self.get_detgroup(det)
2456 return group['energy'][()]
2458 def get_shape(self):
2459 '''returns NY, NX shape of array data'''
2460 ny, nx, npos = self.xrmmap['positions/pos'].shape
2461 return ny, nx
2463 def get_envvar(self, name=None):
2464 """get environment value by name"""
2465 if self.envvar is None:
2466 self.envvar = {}
2467 env_names = [h5str(a) for a in self.xrmmap['config/environ/name']]
2468 env_vals = [h5str(a) for a in self.xrmmap['config/environ/value']]
2469 for name, val in zip(env_names, env_vals):
2470 name = h5str(name)
2471 val = h5str(val)
2472 try:
2473 fval = float(val)
2474 except:
2475 fval = val
2476 self.envvar[name] = fval
2477 if name is not None:
2478 if name in self.envvar:
2479 return self.envvar[name]
2480 elif name.lower() in self.envvar:
2481 return self.envvar[name.lower()]
2483 def get_incident_energy(self):
2484 """ special case of get_envvar"""
2485 env_names = [h5str(a) for a in self.xrmmap['config/environ/name']]
2486 env_vals = [h5str(a) for a in self.xrmmap['config/environ/value']]
2487 for name, val in zip(env_names, env_vals):
2488 name = name.lower().replace('.', ' ')
2489 if name.startswith('mono energy'):
2490 return float(val)
2491 return DEFAULT_XRAY_ENERGY
2493 def get_counts_rect(self, ymin, ymax, xmin, xmax, mapdat=None,
2494 det=None, dtcorrect=None):
2495 '''return counts for a map rectangle, optionally
2496 applying area mask and deadtime correction
2498 Parameters
2499 ---------
2500 ymin : int low y index
2501 ymax : int high y index
2502 xmin : int low x index
2503 xmax : int high x index
2504 mapdat : optional, None or map data
2505 det : optional, None or int index of detector
2506 dtcorrect : optional, bool [None] dead-time correct data
2508 Returns
2509 -------
2510 ndarray for XRF counts in rectangle
2512 Does *not* check for errors!
2514 Note: if mapdat is None, the map data is taken from the 'det' parameter
2515 '''
2516 if dtcorrect is None:
2517 dtcorrect = self.dtcorrect
2518 if mapdat is None:
2519 mapdat = self.get_detgroup(det)
2521 if 'counts' not in mapdat:
2522 mapdat = self.get_detgroup(None)
2524 nx, ny = (xmax-xmin, ymax-ymin)
2525 sx = slice(xmin, xmax)
2526 sy = slice(ymin, ymax)
2528 if len(mapdat['counts'].shape) == 4:
2529 counts = mapdat['counts'][sy, sx, :, :]
2530 else:
2531 counts = mapdat['counts'][sy, sx, :]
2533 if dtcorrect and 'dtfactor' in mapdat:
2534 counts = counts*mapdat['dtfactor'][sy, sx].reshape(ny, nx, 1)
2535 return counts
2537 def get_mca_area(self, areaname, det=None, dtcorrect=None):
2538 '''return XRF spectra as MCA() instance for
2539 spectra summed over a pre-defined area
2541 Parameters
2542 ---------
2543 areaname : str name of area
2544 dtcorrect : optional, bool [None] dead-time correct data
2546 Returns
2547 -------
2548 MCA object for XRF counts in area
2550 '''
2551 try:
2552 area = self.get_area(areaname)[()]
2553 except:
2554 raise GSEXRM_Exception("Could not find area '%s'" % areaname)
2555 if dtcorrect is None:
2556 dtcorrect = self.dtcorrect
2557 npixels = area.sum()
2558 if npixels < 1:
2559 return None
2561 dgroup = self.get_detname(det)
2563 # first get data for bounding rectangle
2564 _ay, _ax = np.where(area)
2565 ymin, ymax, xmin, xmax = _ay.min(), _ay.max()+1, _ax.min(), _ax.max()+1
2566 opts = {'dtcorrect': dtcorrect, 'det': det}
2567 counts = self.get_counts_rect(ymin, ymax, xmin, xmax, **opts)
2568 ltime, rtime = self.get_livereal_rect(ymin, ymax, xmin, xmax, **opts)
2569 ltime = ltime[area[ymin:ymax, xmin:xmax]].sum()
2570 rtime = rtime[area[ymin:ymax, xmin:xmax]].sum()
2571 counts = counts[area[ymin:ymax, xmin:xmax]]
2572 while(len(counts.shape) > 1):
2573 counts = counts.sum(axis=0)
2574 return self._getmca(dgroup, counts, areaname, npixels=npixels,
2575 real_time=rtime, live_time=ltime)
2577 def get_mca_rect(self, ymin, ymax, xmin, xmax, det=None, dtcorrect=None):
2578 '''return mca counts for a map rectangle, optionally
2580 Parameters
2581 ---------
2582 ymin : int low y index
2583 ymax : int high y index
2584 xmin : int low x index
2585 xmax : int high x index
2586 det : optional, None or int index of detector
2587 dtcorrect : optional, bool [None] dead-time correct data
2589 Returns
2590 -------
2591 MCA object for XRF counts in rectangle
2593 '''
2594 if dtcorrect is None:
2595 dtcorrect = self.dtcorrect
2596 dgroup = self.get_detname(det)
2597 mapdat = self.get_detgroup(det)
2598 if 'counts' not in mapdat:
2599 mapdat = self.get_detgroup(None)
2600 dgroup = self.get_detname(None)
2601 counts = self.get_counts_rect(ymin, ymax, xmin, xmax, mapdat=mapdat,
2602 det=det, dtcorrect=dtcorrect)
2603 name = 'rect(y=[%i:%i], x==[%i:%i])' % (ymin, ymax, xmin, xmax)
2604 npix = (ymax-ymin+1)*(xmax-xmin+1)
2605 ltime, rtime = self.get_livereal_rect(ymin, ymax, xmin, xmax, det=det,
2606 dtcorrect=dtcorrect)
2607 counts = counts.sum(axis=0).sum(axis=0)
2608 return self._getmca(dgroup, counts, name, npixels=npix,
2609 real_time=rtime.sum(), live_time=ltime.sum())
2612 def get_livereal_rect(self, ymin, ymax, xmin, xmax, det=None, **kws):
2613 '''return livetime, realtime for a map rectangle, optionally
2614 applying area mask and deadtime correction
2616 Parameters
2617 ---------
2618 ymin : int low y index
2619 ymax : int high y index
2620 xmin : int low x index
2621 xmax : int high x index
2622 det : optional, None or int index of detector
2624 Returns
2625 -------
2626 realtime, livetime in seconds
2628 '''
2629 tshape = self.get_shape()
2630 dmap = self.get_detgroup(det)
2632 if ymax < 0: ymax += tshape[0]
2633 if xmax < 0: xmax += tshape[1]
2634 nx, ny = (xmax-xmin, ymax-ymin)
2635 sx = slice(xmin, xmax)
2636 sy = slice(ymin, ymax)
2638 if 'livetime' in dmap:
2639 livetime = 1.e-6*dmap['livetime'][sy, sx]
2640 realtime = 1.e-6*dmap['realtime'][sy, sx]
2641 else:
2642 livetime = self.pixeltime * np.ones((ny, nx))
2643 realtime = self.pixeltime * np.ones((ny, nx))
2645 return livetime, realtime
2647 def _getmca(self, dgroup, counts, name, npixels=None, **kws):
2648 '''return an MCA object for a detector group
2649 (map is one of the 'mca1', ... 'mcasum')
2650 with specified counts array and a name
2653 Parameters
2654 ---------
2655 det : detector object (one of det1, det2, ..., detsum)
2656 counts : ndarray array of counts
2657 name : name for MCA
2659 Returns
2660 -------
2661 MCA object
2663 '''
2664 map = self.xrmmap[dgroup]
2665 cal = map['energy'].attrs
2666 _mca = MCA(counts=counts, offset=cal['cal_offset'],
2667 slope=cal['cal_slope'], **kws)
2668 if self.incident_energy is None:
2669 self.incident_energy = self.get_incident_energy()
2671 _mca.incident_energy = 0.001*self.incident_energy
2672 _mca.energy = map['energy'][()]
2673 env_names = [h5str(a) for a in self.xrmmap['config/environ/name']]
2674 env_addrs = [h5str(a) for a in self.xrmmap['config/environ/address']]
2675 env_vals = [h5str(a) for a in self.xrmmap['config/environ/value']]
2676 for desc, val, addr in zip(env_names, env_vals, env_addrs):
2677 _mca.add_environ(desc=desc, val=val, addr=addr)
2679 if npixels is not None:
2680 _mca.npixels=npixels
2682 if version_ge(self.version, '2.0.0'):
2683 for roi in self.xrmmap['roimap'][dgroup]:
2684 emin, emax = self.xrmmap['roimap'][dgroup][roi]['limits'][:]
2685 Eaxis = map['energy'][:]
2687 imin = (np.abs(Eaxis-emin)).argmin()
2688 imax = (np.abs(Eaxis-emax)).argmin()
2689 _mca.add_roi(roi, left=imin, right=imax)
2690 else:
2691 # a workaround for poor practice -- some '1.3.0' files
2692 # were built with 'roi_names', some with 'roi_name'
2693 roiname = 'roi_name'
2694 if roiname not in map: roiname = 'roi_names'
2696 roinames = [h5str(a) for a in map[roiname]]
2697 roilim0 = [lim[0] for lim in map['roi_limits']]
2698 roilim1 = [lim[1] for lim in map['roi_limits']]
2699 for roi, lim0, lim1 in zip(roinames, roilim0, roilim1):
2700 _mca.add_roi(roi, left=lim0, right=lim1)
2702 _mca.areaname = _mca.title = name
2703 fname = Path(self.filename).name
2704 _mca.filename = fix_filename(fname)
2705 _mca.info = f"Data from File '{self.filename}', detector '{dgroup}', area '{name}'"
2707 return _mca
2709 def get_xrd1d_area(self, areaname, callback=None, **kws):
2710 '''return 1D XRD pattern for a pre-defined area
2712 Parameters
2713 ---------
2714 areaname : str name of area
2716 Returns
2717 -------
2718 diffraction pattern for given area
2720 '''
2721 try:
2722 area = self.get_area(areaname)[()]
2723 except:
2724 raise GSEXRM_Exception("Could not find area '%s'" % areaname)
2725 return
2726 npix = area.sum()
2727 if npix < 1:
2728 return None
2730 stps, xpix, ypix, qdat = 0, 0, 0, None
2731 sy, sx = [slice(min(_a), max(_a)+1) for _a in np.where(area)]
2732 xmin, xmax, ymin, ymax = sx.start, sx.stop, sy.start, sy.stop
2733 nx, ny = (xmax-xmin), (ymax-ymin)
2735 xrdgroup = 'xrd1d'
2736 mapdat = self.xrmmap[xrdgroup]
2737 counts = self.get_counts_rect(ymin, ymax, xmin, xmax,
2738 mapdat=mapdat, dtcorrect=False)
2739 counts = counts[area[ymin:ymax, xmin:xmax]]
2741 name = '%s: %s' % (xrdgroup, areaname)
2742 kws['energy'] = energy = 0.001 * self.get_incident_energy()
2743 kws['wavelength'] = lambda_from_E(energy, E_units='keV')
2745 counts = counts.sum(axis=0)
2746 xrd = XRD(data1D=counts, steps=len(counts), name=name, **kws)
2747 if xrdgroup != 'xrd1d':
2748 xpix, ypix = counts.shape
2749 xrd = XRD(data2D=counts, xpixels=xpix, ypixels=ypix, name=name, **kws)
2751 fname = Path(self.filename).name
2752 xrd.filename = fname
2753 xrd.areaname = xrd.title = name
2754 xrd.mapname = mapdat.name
2755 xrd.info = f"Data from File '{self.filename}', detector '{mapdat.name}', area '{name}'"
2756 xrd.q = mapdat['q'][()]
2757 return xrd
2759 def get_xrd2d_rect(self, xmin, xmax, ymin, ymax, **kws):
2760 tshape = self.get_shape()
2762 if ymax < 0: ymax += tshape[0]
2763 if xmax < 0: xmax += tshape[1]
2764 nx, ny = (xmax-xmin, ymax-ymin)
2765 sx = slice(xmin, xmax)
2766 sy = slice(ymin, ymax)
2767 xrdgroup = 'xrd2d'
2768 xrdgrp = ensure_subgroup('xrd2d', self.xrmmap, dtype='2DXRD')
2770 xmin, xmax, ymin, ymax = sx.start, sx.stop, sy.start, sy.stop
2772 xrd_file = Path(self.folder, self.rowdata[0][4]).as_posix()
2773 if Path(xrd_file).exists():
2774 print("Reading XRD Patterns for rows %d to %d" %(ymin, ymax-1))
2775 data = None
2776 for yrow in range(ymin, ymax+1):
2777 xrd_file = Path(self.folder, self.rowdata[yrow][4]).as_posix()
2778 print(f"read XRD for row {yrow:d}: {xrd_file:s}")
2779 h5file = h5py.File(xrd_file, 'r')
2780 rowdat = h5file['entry/data/data'][1:,:,:]
2781 h5file.close()
2782 if (yrow % 2) == 1:
2783 rowdat = rowdat[::-1, :, :]
2784 rowdat = rowdat[sx,:,:].sum(axis=0)
2785 if data is None:
2786 data = rowdat*1.0
2787 else:
2788 data += rowdat
2790 aname = f'X[{xmin:d}:{xmax-1:d}]/Y[{ymin:d}:{ymax-1:d}]'
2791 bname = f'X{xmin:03d}Y{ymin:03d}'
2792 name = f'{xrdgroup:s} {aname:s}'
2793 energy = 0.001 * self.get_incident_energy()
2794 kws = {'energy': energy,
2795 'wavelength': lambda_from_E(energy, E_units='keV')}
2797 xrd = XRD(data2D=data, name=name, **kws)
2798 # print("made xrd ", xrd, kws)
2799 xrd.filename = self.fname
2800 xrd.areaname = xrd.title = aname
2801 xrd.info = f"Data from File '{self.filename}', XRD2D '{aname}'"
2802 xrd.ponifile = self.xrdcalfile
2803 return xrd
2805 def get_xrd2d_area(self, areaname, callback=None, **kws):
2806 '''return 2D XRD pattern for a pre-defined area
2808 Parameters
2809 ---------
2810 areaname : str name of area
2812 Returns
2813 -------
2814 diffraction pattern for given area
2816 Notes
2817 ------
2818 slow because it really reads from the raw XRD h5 files
2819 '''
2820 try:
2821 area = self.get_area(areaname)[()]
2822 except:
2823 raise GSEXRM_Exception("Could not find area '%s'" % areaname)
2824 return
2825 npix = area.sum()
2826 if npix < 1:
2827 return None
2829 xrdgroup = 'xrd2d'
2830 xrdgrp = ensure_subgroup('xrd2d', self.xrmmap, dtype='2DXRD')
2832 stps, xpix, ypix, qdat = 0, 0, 0, None
2833 sy, sx = [slice(min(_a), max(_a)+1) for _a in np.where(area)]
2834 xmin, xmax, ymin, ymax = sx.start, sx.stop, sy.start, sy.stop
2835 nx, ny = (xmax-xmin), (ymax-ymin)
2836 xrd_file = Path(self.folder, self.rowdata[0][4]).as_posix()
2837 if Path(xrd_file).exists():
2838 print("Reading XRD Patterns for rows %d to %d" %(ymin+1, ymax))
2839 data = None
2840 for yrow in range(ymin, ymax+1):
2841 xrd_file = Path(self.folder, self.rowdata[yrow][4]).as_posix()
2842 print(f"read XRD for row {yrow:d}: {xrd_file:s}")
2843 h5file = h5py.File(xrd_file, 'r')
2844 rowdat = h5file['entry/data/data'][1:,:,:]
2845 h5file.close()
2846 if (yrow % 2) == 1:
2847 rowdat = rowdat[::-1, :, :]
2848 rowdat = rowdat[np.where(area[yrow])[0], :, :].sum(axis=0)
2849 if data is None:
2850 data = rowdat * 1.0
2851 else:
2852 data += rowdat * 1.0
2854 name = '%s: %s' % (xrdgroup, areaname)
2855 kws = {}
2856 kws['energy'] = energy = 0.001 * self.get_incident_energy()
2857 kws['wavelength'] = lambda_from_E(energy, E_units='keV')
2858 # print("MAKE XRD ", data.shape, data.dtype, data.min(),
2859 # data.max(), data.mean())
2860 xrd = XRD(data2D=data, name=name, **kws)
2861 # print("made xrd ", xrd, kws)
2862 fname = Path(self.filename).name
2863 xrd.filename = fname
2864 xrd.areaname = xrd.title = areaname
2865 xrd.info = f"Data from File '{self.filename}', XRD 2d, area '{areaname}'"
2866 xrd.ponifile = self.xrdcalfile
2867 return xrd
2870 def get_pos(self, name, mean=True):
2871 '''return position by name (matching 'roimap/pos_name' if
2872 name is a string, or using name as an index if it is an integer
2874 Parameters
2875 ---------
2876 name : str ROI name
2877 mean : optional, bool [True] return mean x-value
2879 with mean=True, and a positioner in the first two position,
2880 returns a 1-d array of mean x-values
2882 with mean=False, and a positioner in the first two position,
2883 returns a 2-d array of x values for each pixel
2884 '''
2885 index = -1
2886 if isinstance(name, int):
2887 index = name
2888 else:
2889 for ix, nam in enumerate(self.xrmmap['positions/name']):
2890 if bytes2str(nam.lower()) == name.lower():
2891 index = ix
2892 break
2893 if index == -1:
2894 raise GSEXRM_Exception("Could not find position %s" % repr(name))
2895 pos = self.xrmmap['positions/pos'][:, :, index]
2896 if index in (0, 1) and mean:
2897 pos = pos.sum(axis=index)/pos.shape[index]
2898 return pos
2901 def build_xrd_roimap(self, xrd='1d'):
2902 detname = None
2903 xrdtype = 'xrd%s detector' % xrd
2905 roigroup = ensure_subgroup('roimap', self.xrmmap, dtype='roi maps')
2906 for det, grp in self.xrmmap.items():
2907 if bytes2str(grp.attrs.get('type', '')).startswith(xrdtype):
2908 detname = det
2909 ds = ensure_subgroup(det, roigroup)
2910 ds.attrs['type'] = xrdtype
2911 return roigroup, detname
2913 def add_xrd2droi(self, xyrange, roiname, unit='pixels'):
2915 if version_ge(self.version, '2.0.0'):
2916 if not self.has_xrd2d:
2917 return
2919 roigroup, detname = self.build_xrd_roimap(xrd='2d')
2920 xrmdet = self.xrmmap[detname]
2922 if roiname in roigroup[detname]:
2923 raise ValueError("Name '%s' exists in 'roimap/%s' arrays." % (roiname,detname))
2925 xyrange = [int(x) for x in xyrange]
2926 xmin,xmax,ymin,ymax = xyrange
2928 xrd2d_counts = xrmdet['counts'][:,:,slice(xmin,xmax),slice(ymin,ymax)]
2929 xrd2d_counts = xrd2d_counts.sum(axis=2).sum(axis=2)
2930 if abs(xmax-xmin) > 0 and abs(ymax-ymin) > 0:
2931 xrd2d_cor = xrd2d_counts/(abs(xmax-xmin)*abs(ymax-ymin))
2932 else:
2933 xrd2d_cor = xrd2d_counts
2935 self.save_roi(roiname,detname,xrd2d_counts,xrd2d_cor,xyrange,'area','pixels')
2936 self.get_roi_list(detname, force=True)
2938 else:
2939 print('Only compatible with newest hdf5 mapfile version.')
2942 def read_xrd1d_ROIFile(self,filename,verbose=False):
2944 roidat = readROIFile(filename,xrd=True)
2945 print('\nReading 1D-XRD ROI file: %s' % filename)
2946 for iroi, label, xunit, xrange in roidat:
2947 print(' Adding ROI: %s' % label)
2948 self.add_xrd1droi(label, xrange, unit=xunit)
2949 print('Finished.\n')
2951 def add_xrd1droi(self, roiname, xrange, unit='q', subtract_bkg=False):
2952 if not version_ge(self.version, '2.0.0'):
2953 print('Only compatible with newest hdf5 mapfile version.')
2954 return
2956 if not self.has_xrd1d:
2957 print('No 1D-XRD data in file')
2958 return
2960 if self.incident_energy is None:
2961 self.incident_energy = self.get_incident_energy()
2962 xrange = np.array(xrange)
2963 if unit.startswith('2th'): ## 2th to 1/A
2964 qrange = q_from_twth(xrange,
2965 lambda_from_E(self.incident_energy, E_units='eV'))
2966 elif unit == 'd': ## A to 1/A
2967 qrange = q_from_d(xrange)
2968 else:
2969 qrange = xrange
2971 roigroup, detname = self.build_xrd_roimap(xrd='1d')
2973 if roiname in roigroup[detname]:
2974 raise ValueError("Name '%s' exists in 'roimap/%s' arrays." % (roiname, detname))
2976 counts = self.xrmmap[detname]['counts']
2977 q = self.xrmmap[detname]['q'][:]
2979 imin = (np.abs(q-qrange[0])).argmin()
2980 imax = (np.abs(q-qrange[1])).argmin()+1
2982 xrd1d_sum = xrd1d_cor = counts[:, :, imin:imax].sum(axis=2)
2983 if subtract_bkg and imax > imin:
2984 ibkglo = max(0, imin-3)
2985 ibkghi = min(len(q), imax+3)
2986 bkglo = counts[:,:,ibkglo:imin].sum(axis=2)/(imin-ibkglo)
2987 bkghi = counts[:,:,imax::ibkghi].sum(axis=2)/(ibkghi-imax)
2988 xrd1d_cor -= (imax-imin)* (bkglo + bkghi)/2.0
2990 self.save_roi(roiname, detname, xrd1d_sum, xrd1d_cor,
2991 qrange,'q', '1/A')
2992 self.get_roi_list(detname, force=True)
2995 def del_all_xrd1droi(self):
2997 ''' delete all 1D-XRD ROI'''
2999 roigrp_xrd1d = ensure_subgroup('xrd1d', self.xrmmap['roimap'])
3001 for roiname in roigrp_xrd1d.keys():
3002 self.del_xrd1droi(roiname)
3004 def del_xrd1droi(self, roiname):
3006 ''' delete a 1D-XRD ROI'''
3008 roigrp_xrd1d = ensure_subgroup('xrd1d', self.xrmmap['roimap'])
3010 if roiname not in roigrp_xrd1d.keys():
3011 print("No ROI named '%s' found to delete" % roiname)
3012 return
3014 roiname = h5str(roiname)
3015 if roiname in roigrp_xrd1d:
3016 del roigrp_xrd1d[roiname]
3017 self.h5root.flush()
3020 def save_roi(self, roiname, det, raw, cor, drange, dtype, units):
3021 ds = ensure_subgroup(roiname, self.xrmmap['roimap'][det])
3022 ds.create_dataset('raw', data=raw )
3023 ds.create_dataset('cor', data=cor )
3024 ds.create_dataset('limits', data=drange )
3025 ds['limits'].attrs['type'] = dtype
3026 ds['limits'].attrs['units'] = units
3028 self.h5root.flush()
3030 def build_mca_roimap(self):
3031 det_list = []
3032 sumdet = None
3033 roigroup = ensure_subgroup('roimap', self.xrmmap, dtype='roi map')
3034 for det, grp in zip(self.xrmmap.keys(),self.xrmmap.values()):
3035 if bytes2str(grp.attrs.get('type', '')).startswith('mca det'):
3036 det_list += [det]
3037 ds = ensure_subgroup(det, roigroup)
3038 ds.attrs['type'] = 'mca detector'
3039 if (bytes2str(grp.attrs.get('type', '')).startswith('virtual mca')
3040 and 'sum' in bytes2str(grp.attrs.get('desc', ''))):
3041 sumdet = det
3042 ds = ensure_subgroup(det,roigroup)
3043 ds.attrs['type'] = 'virtual mca detector'
3045 return roigroup, det_list, sumdet
3047 def add_xrfroi(self, roiname, Erange, unit='keV'):
3048 if not self.has_xrf:
3049 return
3051 if unit == 'eV':
3052 Erange[:] = [x/1000. for x in Erange] ## eV to keV
3054 roigroup, det_list, sumdet = self.build_mca_roimap()
3055 if sumdet not in det_list:
3056 det_list.append(sumdet)
3058 for det in det_list:
3059 if roiname in self.xrmmap['roimap'][det]:
3060 print(f"ROI '{roiname:s}' exists for detector '{det:s}'. Delete before adding")
3061 continue
3062 mapdat = self.xrmmap[det]
3063 if unit.startswith('chan'):
3064 emin, emax = Erange
3065 else:
3066 en = mapdat['energy'][:]
3067 emin = (np.abs(en-Erange[0])).argmin()
3068 emax = (np.abs(en-Erange[1])).argmin()+1
3069 raw = mapdat['counts'][:, :, emin:emax].sum(axis=2)
3070 cor = raw * mapdat['dtfactor']
3071 self.save_roi(roiname, det, raw, cor, Erange, 'energy', unit)
3072 self.get_roi_list('mcasum', force=True)
3074 def del_xrfroi(self, roiname):
3075 roigroup, det_list, sumdet = self.build_mca_roimap()
3076 if sumdet not in det_list:
3077 det_list.append(sumdet)
3079 for det in det_list:
3080 if roiname in self.xrmmap['roimap'][det]:
3081 del self.xrmmap['roimap'][det][roiname]
3082 self.get_roi_list('mcasum', force=True)
3083 self.h5root.flush()
3086 def check_roi(self, roiname, det=None, version=None):
3087 if version is None:
3088 version = self.version
3090 if (type(det) is str and det.isdigit()) or type(det) is int:
3091 det = int(det)
3092 detname = 'mca%i' % det
3093 else:
3094 detname = det
3096 if roiname is not None: roiname = roiname.lower()
3097 # print("Check ROI ", roiname, detname, version)
3099 if version_ge(version, '2.0.0'):
3100 for d in self.get_detector_list():
3101 if detname.lower() == d.lower():
3102 for rname in self.xrmmap[d]:
3103 if roiname.lower() == rname.lower():
3104 return rname, d
3106 if detname is not None:
3107 detname = detname.replace('det', 'mca')
3109 if detname is None:
3110 detname = 'roimap/mcasum'
3111 elif not detname.startswith('roimap'):
3112 detname = 'roimap/%s' % detname
3114 try:
3115 roi_list = [r for r in self.xrmmap[detname]]
3116 if roiname is None:
3117 return roi_list, detname
3118 if roiname not in roi_list:
3119 for roi in roi_list:
3120 if roi.lower().startswith(roiname):
3121 roiname = roi
3122 except:
3123 ## provide summed output counts if fail
3124 detname = 'roimap/mcasum'
3125 roiname = 'outputcounts'
3127 else:
3128 if detname is not None:
3129 detname = detname.replace('mca','det')
3131 sum_roi = [h5str(r).lower() for r in self.xrmmap['roimap/sum_name']]
3132 det_roi = [h5str(r).lower() for r in self.xrmmap['roimap/det_name']]
3134 if roiname not in sum_roi:
3135 if roiname is None:
3136 return np.arange(len(sum_roi)),'roimap/sum_'
3137 else:
3138 for roi in sum_roi:
3139 if roi.startswith(roiname):
3140 roiname = roi
3142 if detname in ['det%d' % (i+1) for i in range(self.nmca)]:
3143 idet = int(''.join([i for i in detname if i.isdigit()]))
3144 detname = 'roimap/det_'
3146 if roiname not in det_roi:
3147 roiname = '%s (mca%i)' % (roiname,idet)
3148 roiname = det_roi.index(roiname)
3150 else:
3151 detname = 'roimap/sum_'
3152 try:
3153 roiname = sum_roi.index(roiname)
3154 except:
3155 roiname = sum_roi.index('outputcounts')
3157 return roiname, detname
3160 def get_roimap(self, roiname, det=None, hotcols=None, zigzag=None,
3161 dtcorrect=None, minval=None, maxval=None):
3162 '''extract roi map for a pre-defined roi by name
3163 Parameters
3164 ---------
3165 roiname : str ROI name
3166 det : str detector name
3167 dtcorrect : optional, bool [None] dead-time correct data
3168 hotcols : optional, bool [None] suppress hot columns
3169 minval: float, trim to minimum value
3170 maxval: float, trim to maximum value
3172 Returns
3173 -------
3174 ndarray for ROI data
3175 '''
3176 if hotcols is None:
3177 hotcols = self.hotcols
3178 if zigzag is None:
3179 zigzag = self.zigzag
3180 if dtcorrect is None:
3181 dtcorrect = self.dtcorrect
3183 nrow, ncol, npos = self.xrmmap['positions']['pos'].shape
3184 out = np.zeros((nrow, ncol))
3186 det = self.get_detname(det)
3187 dtcorrect = dtcorrect and ('mca' in det or 'det' in det)
3189 if roiname == '1' or roiname == 1:
3190 out = np.ones((nrow, ncol))
3191 if hotcols:
3192 out = out[1:-1]
3193 return out
3195 roi, detaddr = self.check_roi(roiname, det)
3196 ext = ''
3197 if detaddr.startswith('roimap'):
3198 ext = 'raw'
3199 if dtcorrect:
3200 ext = 'cor'
3203 # print("GetROIMAP roiname=%s|roi=%s|det=%s" % (roiname, roi, det))
3204 # print("detaddr=%s|ext=%s|version=%s" % (detaddr, ext, self.version))
3205 if version_ge(self.version, '2.0.0'):
3206 if detaddr.startswith('roimap'):
3207 roi_ext = '%s/' + ext
3208 else:
3209 roi_ext = '%s_' + ext if ext == 'raw' else '%s'
3210 roiaddr = roi_ext % roi
3211 # print("looking for detattr, roiaddr ", detaddr, roiaddr)
3212 try:
3213 out = self.xrmmap[detaddr][roiaddr][:]
3214 except (KeyError, OSError):
3215 _roiname, _roic = roiaddr.split('/')
3216 try:
3217 del self.xrmmap[detaddr][_roiname]
3218 except:
3219 pass
3220 rgrp = self.xrmmap[detaddr].create_group(_roiname)
3221 for aname,dtype in (('raw', np.uint32),
3222 ('cor', np.float32)):
3223 rgrp.create_dataset(aname, (1, ncol), dtype,
3224 chunks=(1, ncol),
3225 maxshape=(None, ncol), **self.compress_args)
3226 lmtgrp = rgrp.create_dataset('limits', data=[0., 0.], **self.compress_args)
3227 lmtgrp.attrs['type'] = 'energy'
3228 lmtgrp.attrs['units'] = 'keV'
3230 out = np.zeros([1, ncol])
3231 # print("found roi data ", out.shape, nrow, ncol)
3232 if version_ge(self.version, '2.1.0') and out.shape != (nrow, ncol):
3233 _roi, _detaddr = self.check_roi(roiname, det, version='1.0.0')
3234 detname = '%s%s' % (_detaddr, ext)
3235 out = self.xrmmap[detname][:, :, _roi]
3236 # print("from v1, got roi map ", _roi, detname, out.shape)
3237 if self.write_access:
3238 self.xrmmap[detaddr][roiaddr].resize((nrow, ncol))
3239 self.xrmmap[detaddr][roiaddr][:, :] = out
3241 else: # version1
3242 if det in EXTRA_DETGROUPS:
3243 detname = "%s/%s" % (det, roiname)
3244 out = self.xrmmap[detname][:,:]
3245 else:
3246 detname = '%s%s' % (detaddr, ext)
3247 out = self.xrmmap[detname][:, :, roi]
3249 if zigzag is not None and zigzag != 0:
3250 out = remove_zigzag(out, zigzag)
3251 elif hotcols:
3252 out = out[:, 1:-1]
3253 if minval is not None:
3254 out[np.where(out<minval)] = minval
3255 if maxval is not None:
3256 out[np.where(out>maxval)] = maxval
3257 return out
3260 def get_mca_erange(self, det=None, dtcorrect=None,
3261 emin=None, emax=None, by_energy=True):
3262 '''extract map for an ROI set here, by energy range:
3264 not implemented
3265 '''
3266 pass
3268 def get_rgbmap(self, rroi, groi, broi, det=None, rdet=None, gdet=None, bdet=None,
3269 hotcols=None, dtcorrect=None, scale_each=True, scales=None):
3270 '''return a (NxMx3) array for Red, Green, Blue from named
3271 ROIs (using get_roimap).
3273 Parameters
3274 ----------
3275 rroi : str name of ROI for red channel
3276 groi : str name of ROI for green channel
3277 broi : str name of ROI for blue channel
3278 det : optional, None or int [None] index for detector
3279 dtcorrect : optional, bool [None] dead-time correct data
3280 hotcols : optional, bool [None] suppress hot columns
3281 scale_each : optional, bool [True]
3282 scale each map separately to span the full color range.
3283 scales : optional, None or 3 element tuple [None]
3284 multiplicative scale for each map.
3286 By default (scales_each=True, scales=None), each map is scaled by
3287 1.0/map.max() -- that is 1 of the max value for that map.
3289 If scales_each=False, each map is scaled by the same value
3290 (1/max intensity of all maps)
3292 '''
3293 if hotcols is None:
3294 hotcols = self.hotcols
3295 if dtcorrect is None:
3296 dtcorrect = self.dtcorrect
3297 if det is not None:
3298 rdet = gdet = bdet = det
3300 kws = dict(hotcols=hotcols, dtcorrect=dtcorrect)
3301 rmap = self.get_roimap(rroi, det=rdet, **kws)
3302 gmap = self.get_roimap(groi, det=gdet, **kws)
3303 bmap = self.get_roimap(broi, det=bdet, **kws)
3305 if scales is None or len(scales) != 3:
3306 scales = (1./rmap.max(), 1./gmap.max(), 1./bmap.max())
3307 if scale_each:
3308 rmap *= scales[0]
3309 gmap *= scales[1]
3310 bmap *= scales[2]
3311 else:
3312 scale = min(scales[0], scales[1], scales[2])
3313 rmap *= scale
3314 bmap *= scale
3315 gmap *= scale
3317 return np.array([rmap, gmap, bmap]).swapaxes(0, 2).swapaxes(0, 1)
3319 def add_roi(self, name, high, low, address='', det=1,
3320 overwrite=False, **kws):
3321 '''add named ROI to an XRMMap file.
3322 These settings will be propogated through the
3323 ROI maps and all detectors.
3325 '''
3326 # data structures affected:
3327 # config/rois/address
3328 # config/rois/name
3329 # config/rois/limits
3330 # roimap/det_address
3331 # roimap/det_name
3332 # roimap/det_raw
3333 # roimap/det_cor
3334 # roimap/sum_list
3335 # roimap/sum_name
3336 # roimap/sum_raw
3337 # roimap/sum_cor
3338 # det{I}/roi_address for I = 1, N_detectors (xrmmap attribute)
3339 # det{I}/roi_name for I = 1, N_detectors (xrmmap attribute)
3340 # det{I}/roi_limits for I = 1, N_detectors (xrmmap attribute)
3341 # detsum/roi_address for I = 1, N_detectors (xrmmap attribute)
3342 # detsum/roi_name for I = 1, N_detectors (xrmmap attribute)
3343 # detsum/roi_limits for I = 1, N_detectors (xrmmap attribute)
3345 ''
3346 roi_names = [i.lower().strip() for i in self.xrmmap['config/rois/name']]
3347 if name.lower().strip() in roi_names:
3348 if overwrite:
3349 self.del_roi(name)
3350 else:
3351 print("An ROI named '%s' exists, use overwrite=True to overwrite" % name)
3352 return
3354 def del_roi(self, name):
3355 ''' delete an ROI'''
3356 roi_names = [i.lower().strip() for i in self.xrmmap['config/rois/name']]
3357 if name.lower().strip() not in roi_names:
3358 print("No ROI named '%s' found to delete" % name)
3359 return
3360 iroi = roi_names.index(name.lower().strip())
3361 roi_names = [i in self.xrmmap['config/rois/name']]
3362 roi_names.pop(iroi)
3365def read_xrmmap(filename, root=None, **kws):
3366 '''read GSE XRF FastMap data from HDF5 file or raw map folder'''
3367 key = 'filename'
3368 if Path(filename).is_dir():
3369 key = 'folder'
3370 kws.update({key: filename, 'root': root})
3372 return GSEXRM_MapFile(**kws)
3374def process_mapfolder(path, take_ownership=False, **kws):
3375 """process a single map folder
3376 with optional keywords passed to GSEXRM_MapFile
3377 """
3378 try:
3379 kws['xrdcal'] = kws.pop('poni')
3380 except:
3381 pass
3382 if Path(path).is_dir() and isGSEXRM_MapFolder(path):
3383 print( '\n build map for: %s' % path)
3384 try:
3385 g = GSEXRM_MapFile(folder=path, **kws)
3386 except:
3387 print( 'Could not create MapFile')
3388 print( sys.exc_info() )
3389 return
3390 try:
3391 if take_ownership:
3392 g.take_ownership()
3393 if g.check_ownership():
3394 g.process()
3395 else:
3396 print( 'Skipping file %s: not owner' % path)
3397 except KeyboardInterrupt:
3398 sys.exit()
3399 except:
3400 print( 'Could not convert %s' % path)
3401 print( sys.exc_info() )
3402 return
3403 finally:
3404 g.close()
3406def process_mapfolders(folders, ncpus=None, take_ownership=False, **kws):
3407 """process a list of map folders
3408 with optional keywords passed to GSEXRM_MapFile
3409 """
3410 try:
3411 kws['xrdcal'] = kws.pop('poni')
3412 except:
3413 pass
3414 if ncpus is None:
3415 ncpus = max(1, mp.cpu_count()-1)
3416 if ncpus == 0:
3417 for path in folders:
3418 process_mapfolder(path, **kws)
3419 else:
3420 pool = mp.Pool(ncpus)
3421 kws['take_ownership'] = take_ownership
3422 myfunc = partial(process_mapfolder, **kws)
3423 pool.map(myfunc, folders)