Important

You can run this notebook in a live session or view it on nbviewer or GitHub.

# Monochromator Optimization¶

In this notebook you will:

• Use a custom plan that is tuned to a specific, common task at an XAFS (X-ray Absorbance Fine Structure) beamline.

• This plan will take some readings, do some prompt analysis on those readings, and then immediately take an action based on the result.

## Science Background¶

Here is 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.

The first crystal monochromates the X-ray beam by Bragg diffraction. The crystal is on a high-preceision rotation stage and the angle between the incident beam and the lattice planes of the crystal is chosen so that a specific wavelength meets the Bragg condition, $$\lambda = 2d\sin(\Theta)$$, where $$d$$ is the spacing between the lattice planes of the crystal and $$\Theta$$ is the angle between the incident beam and the lattice planes of the first crystal. 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 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 needs to be 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 first and second crystals must be parallel, 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 – it is prudent to verify that this difference in angle is small. This is done by making a scan of the pitch of the second crustal, 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 onochromator is maximized. Thus this pitch scan of the second crystal will be a peaked function. We want to place the second crystal pitch at the maximum of this peak.

## Running a beamline¶

BlueSky provides all the tools we need to perform this optimization. Here’s the game plan:

1. Perform a scan of the pitch motor over a range that will include the optimal pitch position

2. Record an intensity signal from a detector downstream of the DCM

3. When the scan is finished, analyze the intensity measurements and determine the position of peak intensity

4. Move the second crystal pitch motor to its optimized position

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:

• The name of the motor to be scanned

• The name of the detector to be monitored

• The syntax of the mathematical tools used to analyze the measured data

• The names of the standard plans for scanning and moving motors

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 modestly complicated 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.

## Tutorial¶

We start by setting up our BlueSky environment:

[1]:

%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)

/home/travis/virtualenv/python3.7.1/lib/python3.7/site-packages/pims/image_reader.py:26: RuntimeWarning: PIMS image_reader.py could not find scikit-image. Falling back to matplotlib's imread(), which uses floats instead of integers. This may break your scripts.
(To ignore this warning, include the line "warnings.simplefilter("ignore", RuntimeWarning)" in your script.)
warnings.warn(RuntimeWarning(ski_preferred))


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.

[2]:

from plans import rocking_curve

help(rocking_curve)

Help on function rocking_curve in module plans:

rocking_curve(start=-0.1, stop=0.1, nsteps=101, choice='peak')
Perform a relative scan of the DCM 2nd crystal pitch around the current
position to find the peak of the crystal rocking curve.  Begin by opening
the hutch slits to 3 mm. At the end, move to the position of maximum
intensity on I0, then return to the hutch slits to their original height.
Input:
start:    (float)  starting position relative to current  [-0.1]
end:      (float)  ending position relative to current    [0.1]
nsteps:   (int)    number of steps                        [101]
choice:   (string) 'peak', fit' or 'com' (center of mass) ['peak']
If choice is fit, the fit is performed using the
SkewedGaussianModel from lmfit, which works pretty well for this
measurement at BMM.  The line shape is a bit skewed due to the
convolution with the slightly misaligned entrance slits.


[3]:

plt.figure()

[3]:

<Figure size 640x480 with 0 Axes>

[4]:

RE(rocking_curve())

rocking curve scan: pitch, I0, -0.100, 0.100, 101 -- starting at 4.000
uid = 12dd79c0-c0e7-462d-b4c4-7ab67348db89, scan_id = 1
Found and moved to peak at 4.04 via method peak

[4]:

('12dd79c0-c0e7-462d-b4c4-7ab67348db89',)


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 second crystal pitch to that position. At this point, the DCM is optimized and the beamline is ready to begin measuring 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.

[5]:

plt.gcf()  # Display a snapshot of the current state of the figure.

[5]:


By default, the rocking_curve() plan simply looks for the position of highest intesity, then moves to that position. In practice, this works quite well 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:

1. peak: find the position of maximum intensity

2. com: find the position of the center of mass of the masured peak

3. fit: 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 vairable will contain the data from the rocking_curve() performed above.

[6]:

run = catalog[-1]
run

12dd79c0-c0e7-462d-b4c4-7ab67348db89:
args:
entry: !!python/object:databroker.core.Entry
args: []
cls: databroker.core.Entry
kwargs:
name: 12dd79c0-c0e7-462d-b4c4-7ab67348db89
description: {}
driver: databroker.core.BlueskyRunFromGenerator
direct_access: forbid
args:
gen_args: !!python/tuple
- /home/travis/.local/share/bluesky/12dd79c0-c0e7-462d-b4c4-7ab67348db89.msgpack
gen_func: &id003 !!python/name:databroker._drivers.msgpack.gen ''
gen_kwargs: {}
get_filler: &id004 !!python/object/apply:functools.partial
args:
- &id001 !!python/name:event_model.Filler ''
state: !!python/tuple
- *id001
- !!python/tuple []
- handler_registry: !!python/object:event_model.HandlerRegistryView
_handler_registry:
NPY_SEQ: !!python/name:ophyd.sim.NumpySeqHandler ''
newton: !!python/name:bluesky_tutorial_utils._newton.NewtonHandler ''
npy: !!python/name:bluesky_tutorial_utils._old_handlers.NpyHandler ''
npy_FRAMEWISE: !!python/name:bluesky_tutorial_utils._old_handlers.NpyFrameWise ''
inplace: false
root_map: {}
- null
transforms:
descriptor: &id002 !!python/name:databroker.core._no_op ''
resource: *id002
start: *id002
stop: *id002
cache: null
parameters: []
start:
detectors:
- I0
hints:
dimensions:
- - - pitch
- primary
motors:
- pitch
num_intervals: 100
num_points: 101
plan_args:
args:
configuration_attrs=['velocity', 'acceleration'])
- -0.1
- 0.1
detectors:
- SynSignal(name='I0', value=1.6168615132709439, timestamp=1594844987.100452)
num: 101
per_step: None
plan_name: rel_scan
plan_pattern: inner_product
plan_pattern_args:
args:
configuration_attrs=['velocity', 'acceleration'])
- -0.1
- 0.1
num: 101
plan_pattern_module: bluesky.plan_patterns
plan_type: generator
scan_id: 1
time: 1594844987.1752076
uid: 12dd79c0-c0e7-462d-b4c4-7ab67348db89
versions:
bluesky: 1.6.4
ophyd: 1.5.2
stop:
exit_status: success
num_events:
primary: 101
reason: ''
run_start: 12dd79c0-c0e7-462d-b4c4-7ab67348db89
time: 1594844993.5087304
catalog_dir: null
getenv: true
getshell: true
catalog:
cls: databroker._drivers.msgpack.BlueskyMsgpackCatalog
args:
- /home/travis/.local/share/bluesky/*.msgpack
kwargs: {}
gen_args: !!python/tuple
- /home/travis/.local/share/bluesky/12dd79c0-c0e7-462d-b4c4-7ab67348db89.msgpack
gen_func: *id003
gen_kwargs: {}
get_filler: *id004
transforms:
descriptor: *id002
resource: *id002
start: *id002
stop: *id002
description: ''
driver: databroker.core.BlueskyRunFromGenerator
catalog_dir: null
start: !!python/object/new:databroker.core.Start
dictitems:
detectors: &id005
- I0
hints: &id006
dimensions:
- - - pitch
- primary
motors: &id007
- pitch
num_intervals: 100
num_points: 101
plan_args: &id008
args:
configuration_attrs=['velocity', 'acceleration'])
- -0.1
- 0.1
detectors:
- SynSignal(name='I0', value=1.6168615132709439, timestamp=1594844987.100452)
num: 101
per_step: None
plan_name: rel_scan
plan_pattern: inner_product
plan_pattern_args: &id009
args:
configuration_attrs=['velocity', 'acceleration'])
- -0.1
- 0.1
num: 101
plan_pattern_module: bluesky.plan_patterns
plan_type: generator
scan_id: 1
time: 1594844987.1752076
uid: 12dd79c0-c0e7-462d-b4c4-7ab67348db89
versions: &id010
bluesky: 1.6.4
ophyd: 1.5.2
state:
detectors: *id005
hints: *id006
motors: *id007
num_intervals: 100
num_points: 101
plan_args: *id008
plan_name: rel_scan
plan_pattern: inner_product
plan_pattern_args: *id009
plan_pattern_module: bluesky.plan_patterns
plan_type: generator
scan_id: 1
time: 1594844987.1752076
uid: 12dd79c0-c0e7-462d-b4c4-7ab67348db89
versions: *id010
stop: !!python/object/new:databroker.core.Stop
dictitems:
exit_status: success
num_events: &id011
primary: 101
reason: ''
run_start: 12dd79c0-c0e7-462d-b4c4-7ab67348db89
time: 1594844993.5087304
state:
exit_status: success
num_events: *id011
reason: ''
run_start: 12dd79c0-c0e7-462d-b4c4-7ab67348db89
time: 1594844993.5087304



Access the saved data.

[7]:

data = run.primary.read()
data

[7]:

<xarray.Dataset>
Dimensions:                   (time: 101)
Coordinates:
* time                      (time) float64 1.595e+09 1.595e+09 ... 1.595e+09
Data variables:
pitch                     (time) float64 3.9 3.902 3.904 ... 4.096 4.098 4.1
pitch_setpoint            (time) float64 3.9 3.902 3.904 ... 4.096 4.098 4.1
I0                        (time) float64 0.001037 0.001037 ... 0.5701 0.4802
pitch:pitch_velocity      (time) int64 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1
pitch:pitch_acceleration  (time) int64 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1
I0:I0                     (time) float64 0.001037 0.001037 ... 0.001037
seq_num                   (time) int64 1 2 3 4 5 6 7 ... 96 97 98 99 100 101
uid                       (time) <U36 '54c3ae2d-c414-48bb-9a5f-87c890b3ed...

Plot I0 vs pitch.

[8]:

plt.figure()  # New figure

[8]:

<Figure size 640x480 with 0 Axes>

[9]:

data.plot.scatter("pitch", "I0")
plt.gcf()  # Display a snapshot of the current state of the figure.

[9]:


We could have gone straight to the plot in one line by chaining all of this together.

[10]:

catalog[-1].primary.read().plot.scatter("pitch", "I0")
plt.gcf()  # Display a snapshot of the current state of the figure.

[10]:


The variable data contains the result of our just completed rocking curve measurement. The type of this data is xarray Dataset. In the following, we will work instead with a pandas DataFrame. Here, we convert the xarray Dataset to a pandas DataFrame:

[11]:

df = data.to_dataframe()
df

[11]:

pitch pitch_setpoint I0 pitch:pitch_velocity pitch:pitch_acceleration I0:I0 seq_num uid
time
1.594845e+09 3.900 3.900 0.001037 1 1 0.001037 1 54c3ae2d-c414-48bb-9a5f-87c890b3ed06
1.594845e+09 3.902 3.902 0.001037 1 1 0.001037 2 3a9a76bd-4b4d-49c2-b3bc-0140029aea50
1.594845e+09 3.904 3.904 0.001037 1 1 0.001037 3 90ded646-734d-43b3-855a-2c044829e48d
1.594845e+09 3.906 3.906 0.000835 1 1 0.001037 4 a9e7760a-ffd5-4d32-9e54-6af72cd3121c
1.594845e+09 3.908 3.908 0.000925 1 1 0.001037 5 facc7b75-a3fb-401b-9da5-702fee9dc3ba
... ... ... ... ... ... ... ... ...
1.594845e+09 4.092 4.092 0.988044 1 1 0.001037 97 c5bbb236-1c0a-47e9-afdc-f368003f7ec8
1.594845e+09 4.094 4.094 0.843324 1 1 0.001037 98 6584cba9-2120-4032-8fb6-1ef0962d13fc
1.594845e+09 4.096 4.096 0.690857 1 1 0.001037 99 583767fd-88ac-40c0-bcd3-af8faf409b67
1.594845e+09 4.098 4.098 0.570050 1 1 0.001037 100 e258fa0f-b407-480f-a0a8-8fd6a1faea88
1.594845e+09 4.100 4.100 0.480227 1 1 0.001037 101 e8366fbb-3a0a-43fb-874d-735d564f4502

101 rows × 8 columns

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):

[12]:

measurement = df.loc[:, ['I0', 'pitch']]
measurement

[12]:

I0 pitch
time
1.594845e+09 0.001037 3.900
1.594845e+09 0.001037 3.902
1.594845e+09 0.001037 3.904
1.594845e+09 0.000835 3.906
1.594845e+09 0.000925 3.908
... ... ...
1.594845e+09 0.988044 4.092
1.594845e+09 0.843324 4.094
1.594845e+09 0.690857 4.096
1.594845e+09 0.570050 4.098
1.594845e+09 0.480227 4.100

101 rows × 2 columns

### Moving to the peak of the rocking curve¶

First we find the time index of the data point which has the highest I0 value, i.e. the peak of the plot above.

[13]:

i0max = measurement['I0'].idxmax()
i0max

[13]:

1594844991.682537


Now we find the value of pitch at which the I0 intensity was maximum:

[14]:

optimal_pitch = measurement.loc[i0max, 'pitch']
optimal_pitch

[14]:

4.042


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))


### Moving to the center of mass of the rocking curve¶

SciPy’s center of mass calculation is an N-dimensional generalization of the sort of gravitational center of mass calculation you might remember from your mechanics class in physics. 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 we have opted for the simpler solution here.

[15]:

from scipy.ndimage.measurements import center_of_mass
import numpy
arr = numpy.array(measurement['I0'])
val = int(round(center_of_mass(arr)[0]))
val

[15]:

71

[16]:

optimal_pitch = measurement['pitch'].iloc[val]
optimal_pitch

[16]:

4.042


Again, we want to move the pitch motor to its optimal value. So something like like:

RE(mv(pitch, optimal_pitch))


### Fitting a peak function to the rocking curve measurement¶

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.

[17]:

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

[[Model]]
Model(skewed_gaussian)
[[Fit Statistics]]
# fitting method   = leastsq
# function evals   = 64
# data points      = 101
# variables        = 4
chi-square         = 2.05451858
reduced chi-square = 0.02118060
Akaike info crit   = -385.402975
Bayesian info crit = -374.942493
[[Variables]]
amplitude:  0.75779889 +/- 0.00212435 (0.28%) (init = 1.017352)
center:     4.03788514 +/- 0.00480338 (0.12%) (init = 4.042)
sigma:      0.02048252 +/- 0.00108204 (5.28%) (init = 0.022)
gamma:      0.28674007 +/- 0.31991790 (111.57%) (init = 0)
[[Correlations]] (unreported correlations are < 0.000)
C(center, gamma)     = -1.000
C(center, sigma)     = -0.998
C(sigma, gamma)      =  0.998
C(amplitude, sigma)  =  0.074
C(amplitude, center) = -0.041
C(amplitude, gamma)  =  0.039

[17]:

4.037885139917573


Finally, we move the pitch motor to its optimal value. So something like like:

RE(mv(pitch, optimal_pitch))


We can examine the result of the fit by showing the plot:

[18]:

plt.gcf()

[18]:


The peak and center-of-mass algorithms only return the optimal position. An advantage of a fitted analysis of the measurement is that we terms for the amplitide, width, and (in this case) peak assymetry. 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.

## In conclusion¶

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:

RE(rocking_curve(choice='peak'))


Here is the center of mass calculation:

RE(rocking_curve(choice='com'))


and here is 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.

### Hierarchies of plans in BlueSky¶

In truth, 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

1. verifies that the photon delivery system is in the correct configuration for the specified energy range

2. performs a rocking curve scan just like the one in this tutorial

3. optimizes the position of the beam-definition slits

4. makes sure detectors are properly configured to make measurements at the chosen absorption edge

leaving the beamline completely optimized 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 (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 these smaller 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. In fact, at BMM, we allow and encourage our users to run the change_edge() command all by themselves!