In [1]:
ls
In [2]:
# the name of the file that you wish to open
specfilename = '20151111'
# the name of the x column
x = 'MCMY'
# the name of the detector (y column)
y = 'PD21'
# the name of the monitor column
monitor = 'SRcur'
# the scans that you wish to process
scans = [108, 110, 112, 114]
# the name of the output file that you are going to write
output_file_name = '-'.join([specfilename, x, y, monitor]) + '-' + '_'.join([str(scan) for scan in scans])
print('output file name is going to be: %s' % output_file_name)
Specify the kind of interpolation you want to use as a string
'linear'
'nearest'
'zero'
'slinear'
'quadratic
'cubic'
where 'slinear', 'quadratic' and 'cubic' refer to a spline interpolation of first, second or third order) or as an integer specifying the order of the spline interpolator to use.
In [3]:
interpolation_mode = 'linear'
# The number to divide the step size by
# use a value < 1 for more interpolated points
# use a value > 1 for less interpolated points
densify_interpolated_axis_factor = 1
In [4]:
# Some imports that are required for this notebook
import matplotlib.pyplot as plt
from lmfit.models import LorentzianModel, LinearModel
import numpy as np
import os
import pandas as pd
%matplotlib inline
In [5]:
class Specfile:
def __init__(self, filename):
self.filename = os.path.abspath(filename)
with open(self.filename, 'r') as f:
scan_data = f.read().split('#S')
scan_data = [section.split('\n') for section in scan_data]
self.header = scan_data.pop(0)
self.scans = {}
for scan in scan_data:
sid = int(scan[0].split()[0])
self.scans[sid] = Specscan(self, scan)
def __getitem__(self, key):
return self.scans[key]
def __len__(self):
return len(self.scans)-1
def __iter__(self):
return (self.scans[sid] for sid in sorted(self.scans.keys()))
class Specscan:
def __init__(self, specfile, raw_scan_data):
self.specfile = specfile
self.raw_scan_data = raw_scan_data
header_row = self.raw_scan_data.pop(0).split()
self.scan_id = int(header_row.pop(0))
self.scan_command = header_row.pop(0)
self.scan_args = header_row
for row in self.raw_scan_data:
if row.startswith('#L'):
self.col_names = row.split()[1:]
scan_data = [row.split() for row in self.raw_scan_data
if not row.startswith('#') if row]
self.scan_data = pd.DataFrame(data=scan_data, columns=self.col_names, dtype=float)
def __repr__(self):
return 'Specfile("%s")[%s]' % (self.specfile.filename, self.scan_id)
def __str__(self):
return str(self.scan_data)
def __len__(self):
return len(self.scan_data)
def plot(self, column_names=None, x=None):
if x is None:
x = self.scan_data.columns[0]
if column_names is None:
column_names = self.scan_data.columns
ncols = 2
nrows = int(np.ceil(len(column_names)/ncols))
try:
self.ncols
self.nrows
except AttributeError:
self.ncols = 0
self.nrows = 0
if self.ncols != ncols or self.nrows != nrows:
self.ncols, self.nrows = ncols, nrows
self.fig, self.axes = plt.subplots(nrows=nrows,
ncols=ncols,
figsize=(5*ncols, 2*nrows))
self.arts = {}
for data, ax in zip(column_names, self.axes.ravel()):
ax.cla()
self.arts[data] = ax.plot(self.scan_data[x], self.scan_data[data], label=data)
ax.legend(loc=0)
def fit(x, y, bounds=None):
"""Fit a lorentzian + linear background to `field` in `scan`
Parameters
----------
scan : Specscan object
field : The field to fit
bounds : The +/- range to fit the data to
Returns
-------
fit : lmfit.model.ModelFit
The results of fitting the data to a linear + lorentzian peak
Examples
--------
>>> fit = fit_lorentzian(scan.scan_data)
>>> fit.plot()
"""
lorentzian = LorentzianModel()
linear = LinearModel()
center = x[np.argmax(y)]
if bounds is None:
lower, upper = 0, len(x)
else:
lower = center - bounds
upper = center + bounds
if lower < 0:
lower = 0
if upper > len(x):
upper = len(x)
bounds = slice(lower, upper)
# print("Using bounds = {}".format(bounds))
y = y[bounds]
x = x[bounds]
# print("Using x = {}".format(x))
# print("Using y = {}".format(y))
lorentzian_params = lorentzian.guess(y, x=x, center=center)
linear_params = linear.guess(y, x=x)
lorentzian_params.update(linear_params)
model = lorentzian + linear
return model.fit(y, x=x, params=lorentzian_params)
def plotter(xy):
fig, ax = plt.subplots()
arts = {}
for sid, (xdata, ydata) in zip(scans, xy):
arts[sid] = ax.plot(xdata, ydata, '-o', label=sid)
ax.legend(loc=0)
In [6]:
f = Specfile(specfilename)
In [7]:
raw = [(
f[scan_id].scan_data[x].values,
f[scan_id].scan_data[y].values
) for scan_id in scans]
In [8]:
# Use the plotter helper function defined above
plotter(raw)
In [9]:
normalized = [(
x,
y / f[scan_id].scan_data[monitor].values
) for scan_id, (x, y) in zip(scans, raw)]
In [10]:
# Use the plotter helper function defined above
plotter(normalized)
In [11]:
fits = [fit(xdata, ydata) for (xdata, ydata) in normalized]
In [12]:
fig_kws = {'figsize': [15, 15]}
for sid, f in zip(scans, fits):
title_dict = {'title': 'scan_id: %s' % sid}
f.plot(numpoints=len(f.data)*10, ax_res_kws=title_dict,
ax_fit_kws=title_dict, fig_kws=fig_kws)
In [13]:
zeroed = [(np.array(f.userkws['x']-f.params['center'], dtype=float), f.data) for f in fits]
In [14]:
plotter(zeroed)
In [15]:
diff = np.average([np.average(np.diff(x)) for x, y in zeroed])
minval = np.min([np.min(x) for x, y in zeroed])
maxval = np.max([np.max(x) for x, y in zeroed])
In [16]:
diff, minval, maxval
Out[16]:
In [17]:
new_axis = np.arange(minval, maxval, diff / densify_interpolated_axis_factor)
In [18]:
from scipy.interpolate import interp1d
In [19]:
interpolaters = [interp1d(x, y, kind=interpolation_mode,
bounds_error=False,
fill_value=np.nan)
for x, y in zeroed]
In [20]:
# Create a dict of the interpolated values so it can easily be passed to pandas
interpolated = {sid: interpolator(new_axis)
for sid, interpolator in zip(scans, interpolaters)}
In [21]:
df = pd.DataFrame(interpolated, index=new_axis)
In [22]:
df.plot(style='-o')
Out[22]:
In [23]:
df.sum(axis=1).plot(style='-o')
Out[23]:
In [24]:
(df.sum(axis=1) / df.count(axis=1)).plot(style='-o')
Out[24]:
In [25]:
df.dropna().sum(axis=1).plot(style='-o')
Out[25]:
In [26]:
df.dropna().sum(axis=1).to_csv(output_file_name, encoding='ascii', sep=',')
In [27]:
!cat 20151111-MCMY-PD21-SRcur-108_110_112_114
!pwdprint(output_file_name)!pwd and open the file named output_file_name, shown above
In [28]:
!pwd
In [29]:
print(output_file_name)
In [ ]: