#!/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()