In this notebook we demonstrate:
The picture below shows a schematic of a double crystal monochromator (DCM) like the one at NSLS-II's BMM (Beamline for Materials Measurement) and many other beamlines. The broadband radiation from the 3-pole wiggler source is incident upon the first crystal of the DCM. At BMM, we most often use a Si(111) DCM and have the option to use Si(311) crystals. The schematic below applies equally well to either crystal type (and to any other crystal monochrmoator, as well).
Mmonochromation of the X-ray beam happens by Bragg diffraction. The crystal is on a high-resolution rotation stage. The angle between the incident beam and the lattice planes of the crystal is chosen so that a specific wavelength $\lambda$ meets the Bragg condition, $\lambda = 2d\sin(\Theta)$
By changing the angle, we change the wavelength of the beam diffracting from the first crystal. Because there is a simple relationship between wavelength and energy, we select X-ray energy for the experiment by changing the angle of the first crystal.
The second crystal of the DCM is used to direct the beam downstream, towards the experimental hutch. The second crystal must at the same angle as the first crystal in order to meet the Bragg condition for the wavelength selected by the first crystal. In order to pass the X-rays with high efficiency through the DCM, the lattice planes of the first and second crystals must be parallel with accuracy within microradians.
Whenever we make a large change in energy -- for example, when moving between elements with absorption edge energies that are very far apart -- the parallelism of the crystals may not be maintained. So, it is prudent to minimize this difference in angle after making that large energy change. This is done by making a scan of the pitch of the second crystal, monitoring the intensity of the X-ray beam in the experimental hutch. When the second crystal is perfectly parallel to the first crystal, the intensity of the X-rays passing the the monochromator is maximized. Thus this pitch scan of the second crystal will produce a peaked lineshape. We want to place the second crystal pitch at the maximum of this peak.
Bluesky provides all the tools we need to perform this optimization. Here's the game plan:
Bluesky has you covered. It comes with scan plans that perfom steps 1 and 2. Step 3 is readily implemented using tools from Numerical and Scientific Python. Step 4 is handled by one of Bluesky's standard motor motion commands. None of this is difficult ... except that you need to know some things:
For the beamline staff, those things should be common knowledge. For a general user or a new post-doc, those things are esoteric mysteries.
Training a user to remember the names of everything and the sequence of commands to run in order to perform this scan is certainly possible, but certainly challanging. At an XAFS beamline, this optimization might be needed dozens of times each day -- and it must be done correctly *Every. Single. Time.* Even a relatively simple procedure like this optimization becomes a serious friction point in the use of the beamline.
A solution to this friction is make a bespoke measurement plan, tailored to the beamline and constructed from the basic plans provided by Bluesky. That is, we make a plan that performs all the steps of this optimization chore. Thus, instead of training the general user or new post-doc the steps of the optimization, you simply have to train them to run this one plan when it is needed.
In the XAFS community, we usually call this monochromator optimization a "rocking curve scan". (We are "rocking" the second crystal through its optimal position and measuring the resulting peak-shaped curve.) The rocking_curve()
plan explained below is used reliably everyday by users and staff at BMM.
For the sake of completeness, let's talk about the motor and detector being used in this example:
I0
, which is a good, recognizable, semantic name for that signal.We start by setting up our Bluesky environment:
%matplotlib widget
import matplotlib.pyplot as plt
from bluesky import RunEngine
from bluesky_tutorial_utils import setup_data_saving
RE = RunEngine()
catalog = setup_data_saving(RE)
--------------------------------------------------------------------------- ModuleNotFoundError Traceback (most recent call last) <ipython-input-1-4f0398f13e49> in <module> ----> 1 get_ipython().run_line_magic('matplotlib', 'widget') 2 import matplotlib.pyplot as plt 3 from bluesky import RunEngine 4 from bluesky_tutorial_utils import setup_data_saving 5 /usr/lib/python3/dist-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line, _stack_depth) 2315 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals 2316 with self.builtin_trap: -> 2317 result = fn(*args, **kwargs) 2318 return result 2319 <decorator-gen-108> in matplotlib(self, line) /usr/lib/python3/dist-packages/IPython/core/magic.py in <lambda>(f, *a, **k) 185 # but it's overkill for just that one bit of state. 186 def magic_deco(arg): --> 187 call = lambda f, *a, **k: f(*a, **k) 188 189 if callable(arg): /usr/lib/python3/dist-packages/IPython/core/magics/pylab.py in matplotlib(self, line) 97 print("Available matplotlib backends: %s" % backends_list) 98 else: ---> 99 gui, backend = self.shell.enable_matplotlib(args.gui.lower() if isinstance(args.gui, str) else args.gui) 100 self._show_matplotlib_backend(args.gui, backend) 101 /usr/lib/python3/dist-packages/IPython/core/interactiveshell.py in enable_matplotlib(self, gui) 3417 gui, backend = pt.find_gui_and_backend(self.pylab_gui_select) 3418 -> 3419 pt.activate_matplotlib(backend) 3420 pt.configure_inline_support(self, backend) 3421 /usr/lib/python3/dist-packages/IPython/core/pylabtools.py in activate_matplotlib(backend) 318 # when this function runs. 319 # So avoid needing matplotlib attribute-lookup to access pyplot. --> 320 from matplotlib import pyplot as plt 321 322 plt.switch_backend(backend) /usr/lib/python3/dist-packages/matplotlib/pyplot.py in <module> 2347 dict.__setitem__(rcParams, "backend", rcsetup._auto_backend_sentinel) 2348 # Set up the backend. -> 2349 switch_backend(rcParams["backend"]) 2350 2351 # Just to be safe. Interactive mode can be turned on without /usr/lib/python3/dist-packages/matplotlib/pyplot.py in switch_backend(newbackend) 219 else "matplotlib.backends.backend_{}".format(newbackend.lower())) 220 --> 221 backend_mod = importlib.import_module(backend_name) 222 Backend = type( 223 "Backend", (matplotlib.backends._Backend,), vars(backend_mod)) /usr/lib/python3.8/importlib/__init__.py in import_module(name, package) 125 break 126 level += 1 --> 127 return _bootstrap._gcd_import(name[level:], package, level) 128 129 ModuleNotFoundError: No module named 'ipympl'
We will import a Bluesky plan from a script in the current directory, plans.py. The plan operates on simulated hardware defined in another script, simulated_hardware.py. For the purposes of this tutorial we do not need to interact with the hardware directly; it's all done through the plan. You are encouraged to examine plans.py to understand how the rocking curve scan is implemented.
from plans import rocking_curve
help(rocking_curve)
plt.figure()
RE(rocking_curve())
Note that the rocking_curve()
plan performs the scan of the DCM second crystal while monitoring the intensity signal, then analyzes the result to find the position of maximum intensity, and finally moves the pitch of the second crystal to that position. At this point, the DCM is optimized and the beamline is ready to measure XAFS in the range of the current energy of the DCM.
Now let's look at a plot of the rocking curve data using matplotlib's gcf
(i.e. get current figure) method. At the beamline, this is typically displayed as a live plot during the rocking curve scan.
plt.gcf() # Display a snapshot of the current state of the figure.
By default, the rocking_curve()
plan simply looks for the position of highest intesity, then moves to that position. In practice, this is what we do at BMM.
There are other ways to examine and interpret a peak-like function. In fact, the rocking_curve()
plan offers three algorithms for determining the peak position:
peak
: find the position of maximum intensitycom
: find the position of the center of mass of the masured peakfit
: fit an analytic function to the measured data and find the centroid of that function. In practice at BMM, we have found that a skewed Gaussian function works well.Let's see how these work.
First, we need to get the data from the measurement in a form that is convenient to work with. One of the things we set up back in the first step of this tutorial was a "catalog", i.e. a way of accessing the live data from a measurement. In the step that follows, the run
variable will contain the data from the rocking_curve()
performed above.
run = catalog[-1]
run
Access the saved data:
data = run.primary.read()
data
Plot I0 vs pitch:
plt.figure() # New figure
data.plot.scatter("pitch", "I0")
plt.gcf() # Display a snapshot of the current state of the figure.
We could have gone straight to the plot in one line by chaining all of this together:
catalog[-1].primary.read().plot.scatter("pitch", "I0")
plt.gcf() # Display a snapshot of the current state of the figure.
The variable data
contains the result of our just completed rocking curve measurement. The type of this data is an xarray Dataset. In the following, we will work instead with a pandas DataFrame. Here, we convert the xarray Dataset to a pandas DataFrame:
df = data.to_dataframe()
df
Let's take a slice out of that DataFrame so we can focus on the most relevant parts of this data record, i.e. the values of pitch throughout the scan and the beam intensity at each value of pitch (taken from a detector named I0
):
measurement = df.loc[:, ['I0', 'pitch']]
measurement
First we find the time index of the data point which has the highest I0
value, i.e. the peak of the plot above.
i0max = measurement['I0'].idxmax()
i0max
Now we find the value of pitch at which the I0
intensity was maximum:
optimal_pitch = measurement.loc[i0max, 'pitch']
optimal_pitch
Finally, we want to move the pitch motor to its optimal value. In the rocking_curve()
plan actually used at BMM, we have a line like:
RE(mv(pitch, optimal_pitch))
SciPy's center of mass calculation is an N-dimensional generalization of the sort of gravitational center of mass calculation you might remember from a mechanics class taken so long ago.... In this case, the I0
values at each point of the measurement are used as the "mass" of each position in the array. This way of optimizing the DCM second crystal position might be useful if the measured peak is highly assymetric.
The center of mass calculation will return a real number, not an integer. That is, the center of mass will be between two of the actual measured points. The line below with int(round( ... ))
returns the index closest to the center of mass, which we then move to. Alternately, you could interpolate to the position of the actual center of mass, but the simpler solution is shown here.
from scipy.ndimage.measurements import center_of_mass
import numpy
arr = numpy.array(measurement['I0'])
val = int(round(center_of_mass(arr)[0]))
val
optimal_pitch = measurement['pitch'].iloc[val]
optimal_pitch
Again, we want to move the pitch motor to its optimal value. So something like like:
RE(mv(pitch, optimal_pitch))
Sometimes, when interpreting a peak-like function, it is preferable to bring some more prior knwoeldge to the interpretation. If you know a line shape that provides a good physical description of the measurement, then fitting that line shape might provide you with more information.
Here is how we use lmfit to fit a simple skewed Gaussian model to our rocking curve measurement.
First we convert the I0 and pitch columns of the DataFrame to NumPy arrays, which we then feed to lmfit's skewed Gaussian model. From these, we guess initial parameters. We run the fit, then print the results and prepare a plot showing those results. Finally we set the optimal_pitch
parameter for use in the same manner as for the peak and center-of-mass interpretations.
from lmfit.models import SkewedGaussianModel
signal = numpy.array(measurement['I0'])
position = numpy.array(measurement['pitch'])
mod = SkewedGaussianModel()
pars = mod.guess(signal, x=position)
out = mod.fit(signal, pars, x=position)
print(out.fit_report(min_correl=0))
out.plot()
optimal_pitch = out.params["center"].value
optimal_pitch
Finally, we move the pitch motor to its optimal value. Again:
RE(mv(pitch, optimal_pitch))
We can examine the result of the fit by showing the plot:
plt.gcf()
The peak and center-of-mass algorithms only return the optimal position. An advantage of a fitted analysis of the measurement is that we obtain terms for the amplitide, width, and (in this case) peak assymetry, as well as uncertainties. While we do not need those terms for the purpose of optimizing the rocking curve, in other situations that kind of information might be actionable in a plan like this one.
At BMM, we usually run this rocking curve scan using the peak interpretation. That is the default, so the plan is typically run like so:
RE(rocking_curve())
You can specify the kind of intepretation by using the choice
argument. This is equivalent to the default, which is to find the peak of the meaured signal. Here is the explicit call to move to the peak:
RE(rocking_curve(choice='peak'))
Here is how to select the center of mass calculation:
RE(rocking_curve(choice='com'))
and here is how to select the fitted interpretation:
RE(rocking_curve(choice='fit'))
You could go way back to the third step of this tutorial and give each a try. You will find that the three give slightly different answers for the optimized position of the DCM second crystal.
In practice, the rocking curve plan is rarely called directly either by staff or by users at BMM. It is, nonetheless, used many times each day. At BMM, we have a plan called change_edge()
that is used to automate the reconfiguration of the beamline for measurements in different energy ranges. This plan is how the staff and users typically perform the rocking curve chore.
The change_edge()
plan requires an argument specifying the element of the next experiment. So, if we want to begin measuring the iron K edge of iron-bearing samples, we do:
RE(change_edge('Fe'))
This plan
leaving the beamline completely optimized and ready for measurement at the selected absorption edge.
The reason such a magical plan is possible is because, in Bluesky, plans are composed of plans.
Just as the rocking curve plan is built from basic plans supplied by Bluesky (specifically the rel_scan()
and mv()
plans), complex user-defined plans are composed of simpler user-defined plans. At BMM, we have plans for each item change_edge()
chore list. Thus, the change_edge()
plan is composed of several smaller, user-defined plans such as rocking_curve()
.
By building complicated plans out of well-tested, single-purpose plans, we are able to perform automation at the beamline in ways that make the beamline easy to operate for staff and users alike. The details of the steps outlined above are hidden in normal operation, but without hiding their purpose. This facilitates the correct use of the beamline.
In fact, at BMM, we allow and encourage our users to run the change_edge()
command all by themselves!