Source code for PyMca5.PyMcaIO.MEDFile

#!/usr/bin/python
#    Copyright (c) 2010 Matthew Newville, The University of Chicago
#
#   Permission to use and redistribute the source code or binary forms of
#   this software and its documentation, with or without modification is
#   hereby granted provided that the above notice of copyright, these
#   terms of use, and the disclaimer of warranty below appear in the
#   source code and documentation, and that none of the names of The
#   University of Chicago, The University of Washington, or the authors
#   appear in advertising or endorsement of works derived from this
#   software without specific prior written permission from all parties.
#
#   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
#   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
#   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
#   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
#   CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
#   TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
#   SOFTWARE OR THE USE OR OTHER DEALINGS IN THIS SOFTWARE.
#
__license__ = "BSD"
__author__ = "M. Newville - The University of Chicago"
"""
Simple interface to M. River's Multi-Element MCA Data Format

M. Newville
"""
import numpy as np
import re

MIN_SLOPE   = 1.e-7

[docs]def str_converter(strin, delim=None, converter=None): """convert a string of a delimited array to a list""" if delim is None: arr = strin.split() else: arr = re.split(delim, strin) conv = converter if hasattr(conv, '__call__'): return [conv(elem) for elem in arr] else: return arr
[docs]def str2float(strin, delim=None): "string of floats to array of floats" return str_converter(strin, delim=delim, converter=float)
[docs]def str2int(strin, delim=None): "string of integers to array of ints" return str_converter(strin, delim=delim, converter=int)
[docs]def str2str(strin, delim=None): "string to array of strings" return str_converter(strin, delim=delim)
[docs]class ROI(object): "simple Region of Interest" def __init__(self, index=0, left=-1, right=-1, name=None, spectra=None): self.index = index self.left = left self.right = right self.name = name self.spectra = np.array(spectra) self.__counts = -1 def __repr__(self): return "<ROI('%s' chan:[%i, %i])>" % (self.name, self.left, self.right)
[docs] def counts(self): "total counts in roi" if self.spectra is None: return None return self.spectra[self.left:self.right+1].sum()
[docs]class MCA(object): """ basic MCA spectra""" def __init__(self, data=None): self.npts = 1 self.data = data self.energy = None self.realtime = 0 self.livetime = 0 self.deadtime_correction = 1.0 self.cal_offset = 0 self.cal_slope = 0.010 self.cal_quad = 0 self.cal_tth = 0 self.rois = [] def __repr__(self): return "<MCA(%i points) %s>" % (self.npts, hex(id(self)))
[docs] def chan2energy(self, i): "get energy from a channel number" if self.energy is not None: self.get_energy() return self.energy[i]
[docs] def get_calibration(self): "return calibration constants" self.cal_slope = max(MIN_SLOPE, self.cal_slope) return (self.cal_offset, self.cal_slope, self.cal_quad)
[docs] def get_energy(self): "return full energy array" if self.energy is not None: return self.energy idx = np.arange(self.npts) self.cal_slope = max(MIN_SLOPE, self.cal_slope) self.energy = self.cal_offset + idx * (self.cal_slope + idx * self.cal_quad) return self.energy
[docs]class MEDFile(object): """MultiElement XRF Data File Format """ def __init__(self, filename=None): self.default_detector = 0 # "good" detector for energy calibration self.env = [] self.mcas = [] self.filename = filename if filename is not None: self.mca_read_file(filename)
[docs] def get_calibration(self, detector=None): "get calibration constants" if detector is None: detector = self.default_detector return self.mcas[detector].get_calibration()
[docs] def chan2energy(self, i, detector=None): "get energy from a channel number" if detector is None: detector = self.default_detector return self.mcas[detector].chan2energy(i)
[docs] def get_energy(self, detector=None): "get energy array" if detector is None: detector = self.default_detector return self.mcas[detector].get_energy()
[docs] def get_data(self, detector=None, sum_all=True): """ get detector data, if sum_all == False, just the 1 array is returned if sum_all == True, the sum of all detectors is returned, aligned to the energy of the specified detector""" if detector is None: detector = self.default_detector dat = self.mcas[detector].data if sum_all: enref = self.mcas[detector].get_energy() dat = np.zeros(len(enref)) for mca in self.mcas: et = mca.get_energy() dt = mca.data[:] dat = dat + np.interp(enref, et, dt) return dat
[docs] def mca_read_file(self, fname): "read MCA data file" self.filename = fname f = open(fname) lines = f.readlines() f.close() mode = 'HEADER' nelem = 1 # tmp data for data and headers, and rois tmpdat = [] header = {} _roi_0, _roi_1, _roi_n = {}, {}, {} for line in lines: line = line.strip() if len(line) < 1: continue if mode == 'DATA': # data mode tmpdat.append(str2int(line)) else: words = [x.strip() for x in line.split(' ', 1)] if len(words) < 2: # note that 'Data:' line as 1 word. words.append('') tag, val = words[0], words[1] tag = tag.replace(':', '').lower() if tag == 'data': mode = 'DATA' elif tag == 'elements': nelem = int(val) elif tag in ('rois', 'real_time', 'live_time', 'cal_offset', 'cal_slope', 'cal_quad', 'two_theta'): header[tag] = str2float(val) elif tag == 'environment': self.env.append(val) elif tag.startswith('roi_'): x, sroi, item = tag.split('_') iroi = int(sroi) if item == "label": labels = str2str(val, delim='\&') if labels[-1] == '': labels = labels[:-1] _roi_n[iroi] = labels elif item == "left": _roi_0[iroi] = str2int(val) elif item == "right": _roi_1[iroi] = str2int(val) else: header[tag] = val # find first valid detector, identify bad detectors self.mcas = [MCA() for i in range(nelem)] tmpdat = np.transpose(np.array(tmpdat)) for imca, mca in enumerate(self.mcas): mca.npts = int(header['channels']) mca.nrois = int(header['rois'][imca]) mca.start_time = header['date'] mca.realtime = header['real_time'][imca] mca.livetime = header['live_time'][imca] mca.cal_offset = header['cal_offset'][imca] mca.cal_slope = header['cal_slope'][imca] mca.cal_quad = header['cal_quad'][imca] mca.cal_tth = header['two_theta'][imca] mca.data = 1 * tmpdat[imca, :] for iroi in _roi_n: name = _roi_n[iroi][imca].strip() ileft = _roi_0[iroi][imca] iright = _roi_1[iroi][imca] mca.rois.append(ROI(index=iroi, left=ileft, right=iright, name=name, spectra=mca.data))
[docs] def write_ascii(self, fname, elem=None, sum_all=True): """write data to ASCII column file""" out = [] out.append("# XRF data from %s\n" % (self.filename)) if len(self.env)>0: out.append("# Extra PVs:\n") for i in self.env: out.append("# %s\n" % i) out.append("#-------------------------\n") out.append("# energy counts\n") en = self.get_energy() if elem is not None: dat = self.get_data(detector=elem) elif sum_all: dat = self.get_data() for i in ("%8.4f %i\n" % (ei, di) for ei, di in zip(en, dat)): out.append("%s"%i) f = open(fname, "w+") f.writelines(out) f.close()
if __name__ == '__main__': try: import pylab HAS_PYLAB = True except ImportError: HAS_PYLAB = False xrf = MEDFile('test.xrf') energy = xrf.get_energy() d0 = xrf.get_data(detector=0, sum_all=False) d1 = xrf.get_data(detector=1, sum_all=False) d2 = xrf.get_data(detector=2, sum_all=False) d3 = xrf.get_data(detector=3, sum_all=False) dsum = xrf.get_data(detector=0, sum_all=True) xrf.write_ascii('test.dat') print(' ROIs from Element 2:') print(' ------------------') print(' Name | Sum ') for roi in xrf.mcas[1].rois: print(' %s = %d ' % (roi.name, roi.counts())) if HAS_PYLAB: pylab.plot(energy, dsum) pylab.show()