Lesson 5, Part A: find a peak and lineup#
In this lesson, alignment to a narrow diffraction peak is using capabilities provided by the ophyd package. The simulation consists of a simulated motor and simulated noisy detector.
The noisy detector is configured to describe a narrow diffraction peak with Gaussian profile based on the value of the motor position. The peak is centered randomly somewhere between motor values -1 and +1. The width is less than 0.05 in the same units. The peak intensity is expected to be approximately 100,000 (counts/sec are typical units).
Preparation
Make sure the instrument package is in the same directory as this jupyter notebook. The instrument package included with this lesson is a brief version of the standard package used with any APS instrument. Since the notebook is for teaching, it does not connect with any mongodb database. The scans are not kept by the databroker. However, every scan is saved to a SPEC data file as described when the instrument package is loaded.
Objective#
Use Bluesky (the tools provided by the various packages of the Bluesky framework) to find the center and width of a simulated diffraction peak. Move the motor to the peak center.
Use interactive ophyd commands to find the peak and assess the width.
Use the RunEngine and bluesky plans to find the peak and assess the width.
Advanced#
Make a custom plan and run it to find the peak center.
Add that new plan to the instrument package, then restart the notebook’s kernel and try it.
Add the simulated motor and noisy detector to the instrument package, then restart the notebook’s kernel and find the peak again.
Initialize#
Load the instrument controls (which sets up the Bluesky framework for collection: RE, bec, bp, …). This defines more than we need but works as a simple start, just like regular data acquisition at a beamline.
[1]:
from instrument.collection import *
%matplotlib inline
I Thu-11:06:11 - ############################################################ startup
I Thu-11:06:11 - logging started
I Thu-11:06:11 - logging level = 10
I Thu-11:06:11 - c:\Users\Pete\Documents\projects\bluesky_training\lessons\instrument\collection.py
I Thu-11:06:11 - c:\Users\Pete\Documents\projects\bluesky_training\lessons\instrument\mpl\notebook.py
Activating auto-logging. Current session state plus future input saved.
Filename : c:\Users\Pete\Documents\projects\bluesky_training\lessons\.logs\ipython_console.log
Mode : rotate
Output logging : True
Raw input log : False
Timestamping : True
State : active
I Thu-11:06:12 - bluesky framework
I Thu-11:06:12 - c:\Users\Pete\Documents\projects\bluesky_training\lessons\instrument\framework\check_python.py
I Thu-11:06:12 - c:\Users\Pete\Documents\projects\bluesky_training\lessons\instrument\framework\check_bluesky.py
I Thu-11:06:14 - c:\Users\Pete\Documents\projects\bluesky_training\lessons\instrument\framework\initialize.py
I Thu-11:06:17 - c:\Users\Pete\Documents\projects\bluesky_training\lessons\instrument\framework\metadata.py
I Thu-11:06:17 - c:\Users\Pete\Documents\projects\bluesky_training\lessons\instrument\framework\callbacks.py
I Thu-11:06:17 - writing to SPEC file: c:\Users\Pete\Documents\projects\bluesky_training\lessons\20200521-110617.dat
I Thu-11:06:17 - >>>> Using default SPEC file name <<<<
I Thu-11:06:17 - file will be created when bluesky ends its next scan
I Thu-11:06:17 - to change SPEC file, use command: newSpecFile('title')
Numpy provides the random number generator we’ll use.
[2]:
import numpy as np
Load the ophyd simulators
[3]:
from ophyd.sim import motor, SynGauss
Make a noisy detector#
Make a new noisy, replacing the one from the simulator.
[4]:
noisy = SynGauss(
'noisy',
motor, 'motor',
# center somewhere between -1 and 1
center=2 * (np.random.random()-0.5),
# randomize these parameters
Imax=100000 + 20000 * (np.random.random()-0.5),
noise='poisson',
sigma=0.016 + 0.015 * (np.random.random()-0.5),
noise_multiplier=0.1 + 0.02 * (np.random.random()-0.5),
labels={'detectors'})
TODO: work this into the tutorial
[5]:
%ct detectors motors
[This data will not be saved. Use the RunEngine to collect data.]
noisy 0
motor 0
motor_setpoint 0
Define the reported precision of the motor and detector.
[6]:
motor.precision = 5
noisy.precision = 0
Print the values just configured
[7]:
print(f"motor: {motor.position}")
print(f"center: {noisy.center.get()}")
print(f"sigma: {noisy.sigma.get()}")
print(f"Imax: {noisy.Imax.get()}")
print(f"noise: {noisy.noise.get()}")
print(f"noise_multiplier : {noisy.noise_multiplier.get()}")
# tell the "detector" to "count" (get a new value based on the motor position)
%ct noisy
print(f"noisy_det : {noisy.val.get()}")
motor: 0
center: -0.9922233044405593
sigma: 0.02276535181863142
Imax: 96204.44084875335
noise: poisson
noise_multiplier : 0.10293011351852485
[This data will not be saved. Use the RunEngine to collect data.]
noisy_det : 0
1. Use ophyd commands#
The first thing to learn is how to move the motor. As typical of many control systems, there are several ways to do this. We will focus on only two of these since we can apply them to both simulated motors and EPICS motors.
Different from a real motor, our simulated motor moves immediately, so motor velocity and acceleration are not involved. Since the motor moves immediately, there is no discernable delay due to short or long motor motions.
tip: There are several ways to do most things but this tutorial focuses on just a few with the hope that these ways are both simple and reuseable.
motor, move and get position#
For a motor named motor, any of these commands can be used in an interactive session to move the motor from its current position to 1.0:
motor.set(1)
motor.set(1).wait()
%mov motor 1
The %mov command is simplest so that is what we will use here.
tip: The %mov command is absolute move, while %movr is a relative move.
Move the motor from where it is now, to 1, and print its position.
[8]:
%mov motor 1
print(motor.position)
1
[9]:
motor.readback.get()
[9]:
1
Move the motor to 0 (this time using a relative move).
[10]:
%movr motor -1
motor.position
[10]:
0
detector, count and get value#
Above, we created a detector named noisy that simulates a scaler (detector that records a single integer number of detection events, usually the number of X-ray photons received). Real scalers are told to measure for a fixed time interval. Our simulator does not have that feature.
The simulated noisy detector computes its value for counts based on the position of the configured motor.
Show the name of the motor configured into the noisy detector. Check that its name is motor.
[11]:
print(noisy._motor)
SynAxis(prefix='', name='motor', read_attrs=['readback', 'setpoint'], configuration_attrs=['velocity', 'acceleration'])
The value (number of counts) is kept in noisy.val. Show it’s value now.
tip: We can drop the print() wrapper if the command we use returns the value we’d print anyway. Use this convenient shortcut.
[12]:
noisy.val
[12]:
SynSignal(name='noisy', parent='noisy', value=0, timestamp=1590077189.4523768)
We need to tell the detector to acquire data. To acquire data, our simulator will re-compute its value based on the motor position (as with a real detector, the value does not update without something that compels this computation), since that may have changed since the last computation.
For interactive use with ophyd Devices, the command to call is %ct. We labeled the noisy object as detectors so it will be counted when %ct is called.
TODO: improve the labels explanation here.
[13]:
%ct
[This data will not be saved. Use the RunEngine to collect data.]
noisy 0
Find the simulated peak#
With tools to move the motor and acquire data from the detector, we can try to find the simulated peak. It may take some retries since the peak is narrow. Take big steps first (such as 0.1) to find non-zero counts, then smaller steps to find the peak.
First, move to one end of the range, then start stepping until you find non-zero counts. Then use %movr and execute the same notebook cell repeatedly. Call the detector’s .get() method to only print the number of counts and not the other information.
[14]:
%mov motor -1
%ct
print(f"motor={motor.position:.5f} noisy={noisy.val.get()}")
[This data will not be saved. Use the RunEngine to collect data.]
noisy 90261
motor=-1.00000 noisy=90261
The next cell will probably show a very small step size. Change it to 0.1 and execute the cell repeatedly with . Once you have reached a peak (or passed it), change the sign and make the step size smaller. Repeat until you are satisfied.
[15]:
%movr motor .1
%ct
print(f"{motor.position:.5f}, {noisy.val.get()}")
[This data will not be saved. Use the RunEngine to collect data.]
noisy 31
-0.90000, 31
Compare the peak center you found with the value printed after the noisy detector was configured (above). Probably, they will differ by a small amount since the simulator applies random noise to the signal.
2. Use RunEngine and bluesky plans#
Here we use the RunEngine (RE) and standard bluesky plans to locate the simulated diffraction peak.
Since we will do a series of scans, let’s make a list to collect the results from each scan. We’ll report that list later.
[16]:
results=[]
Next, some variables are defined which will make calling the scans more consistent.
The first variable, k, is used to expand (only slightly) the range of the next scan to capture the full width of the peak. We’ll also define n as the number of points in each scan.
[17]:
k = 1.5 # range expansion factor
n = 29 # number of points per scan
Locate the approximate peak position#
Scan from -2 to 2 to find the peak. Since it is a Gaussian (which decays rapidly away from the peak), we may need to increase n, the number of points in the scan. All we need is one point above the background to find it!
Use the scan plan from bluesky.plans (provided here as bp) to locate at least one point on the peak that is above the background of 0 counts.
[18]:
RE(bp.scan([noisy], motor, -2, 2, n))
Transient Scan ID: 1 Time: 2020-05-21 11:10:04
Persistent Unique Scan ID: '000f574f-cd64-4e2e-9375-3816ffef0146'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | motor | noisy |
+-----------+------------+------------+------------+
| 1 | 11:10:04.7 | -2.00000 | 0 |
| 2 | 11:10:04.8 | -1.85714 | 0 |
| 3 | 11:10:04.9 | -1.71429 | 0 |
| 4 | 11:10:05.0 | -1.57143 | 0 |
| 5 | 11:10:05.1 | -1.42857 | 0 |
| 6 | 11:10:05.1 | -1.28571 | 0 |
| 7 | 11:10:05.2 | -1.14286 | 0 |
| 8 | 11:10:05.3 | -1.00000 | 90943 |
| 9 | 11:10:05.4 | -0.85714 | 0 |
| 10 | 11:10:05.5 | -0.71429 | 0 |
| 11 | 11:10:05.6 | -0.57143 | 0 |
| 12 | 11:10:05.7 | -0.42857 | 0 |
| 13 | 11:10:05.8 | -0.28571 | 0 |
| 14 | 11:10:05.9 | -0.14286 | 0 |
| 15 | 11:10:06.0 | 0.00000 | 0 |
| 16 | 11:10:06.1 | 0.14286 | 0 |
| 17 | 11:10:06.2 | 0.28571 | 0 |
| 18 | 11:10:06.3 | 0.42857 | 0 |
| 19 | 11:10:06.3 | 0.57143 | 0 |
| 20 | 11:10:06.4 | 0.71429 | 0 |
| 21 | 11:10:06.6 | 0.85714 | 0 |
| 22 | 11:10:06.7 | 1.00000 | 0 |
| 23 | 11:10:06.7 | 1.14286 | 0 |
| 24 | 11:10:06.8 | 1.28571 | 0 |
| 25 | 11:10:06.9 | 1.42857 | 0 |
| 26 | 11:10:07.0 | 1.57143 | 0 |
| 27 | 11:10:07.1 | 1.71429 | 0 |
| 28 | 11:10:07.2 | 1.85714 | 0 |
| 29 | 11:10:07.3 | 2.00000 | 0 |
+-----------+------------+------------+------------+
generator scan ['000f574f'] (scan num: 1)
[18]:
('000f574f-cd64-4e2e-9375-3816ffef0146',)
One of the tools that works in the background of the Bluesky framework is the BestEffortCallback, known here as bec. When bec is configured (see instrument/framework/initialize.py) as part of scanning with the RunEngine, it will assess peak parameters from each scan, where applicable. The parameters are available in bec.peaks which we have, for convenience, defined as peaks. We access a couple of
those parameters here for peak center and width.
[19]:
peaks
[19]:
{
'com':
{'noisy': -1.0}
,
'cen':
{'noisy': -1.0}
,
'max':
{'noisy': (-1.0,
90943)}
,
'min':
{'noisy': (-2.0,
0)}
,
'fwhm':
{'noisy': 0.1428571428571428}
,
}
We’ll grab the values we need from this dictionary. The term sigma is a measure of the peak width apparent from the data available. Since the step size of the scan is large with respect to the width of the peak shown above, this is a low precision finding. We should repeat this scan with finer step size near the peak to make a more precise assessment.
[20]:
cen = peaks["cen"]["noisy"]
sigma = peaks["fwhm"]["noisy"]
results.append((RE.md["scan_id"], cen, sigma))
print(f"center={cen:.7f}, FWHM={sigma:.7f}")
center=-1.0000000, FWHM=0.1428571
Refine the peak position#
Refine the scan to the range of (-sigma .. +sigma) near the center of the previous scan. Repeat as often as necessary (using ) to get the peak center and width. Use the relative scan plan from bluesky.plans to find the peak.
tip: Look for the plot in the cell above. Replots will be drawn in different colors. The legend indicates the scan_id.
[23]:
%mov motor cen
RE(bp.rel_scan([noisy], motor, -k*sigma, k*sigma, n))
# TODO: plt.gcf()
cen = peaks["cen"]["noisy"]
sigma = peaks["fwhm"]["noisy"]
results.append((RE.md["scan_id"], cen, sigma))
print(f"center={cen:.7f}, FWHM={sigma:.7f}")
Transient Scan ID: 4 Time: 2020-05-21 11:12:24
Persistent Unique Scan ID: '2631ce75-ce77-4da3-acba-4f109b29e29a'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | motor | noisy |
+-----------+------------+------------+------------+
| 1 | 11:12:24.2 | -1.07724 | 94 |
| 2 | 11:12:24.3 | -1.07116 | 250 |
| 3 | 11:12:24.3 | -1.06509 | 555 |
| 4 | 11:12:24.4 | -1.05902 | 1303 |
| 5 | 11:12:24.5 | -1.05294 | 2729 |
| 6 | 11:12:24.5 | -1.04687 | 5285 |
| 7 | 11:12:24.6 | -1.04080 | 9916 |
| 8 | 11:12:24.7 | -1.03472 | 16745 |
| 9 | 11:12:24.7 | -1.02865 | 26837 |
| 10 | 11:12:24.8 | -1.02258 | 39703 |
| 11 | 11:12:24.8 | -1.01651 | 54063 |
| 12 | 11:12:24.9 | -1.01043 | 69659 |
| 13 | 11:12:25.0 | -1.00436 | 83210 |
| 14 | 11:12:25.0 | -0.99829 | 92494 |
| 15 | 11:12:25.1 | -0.99221 | 96544 |
| 16 | 11:12:25.1 | -0.98614 | 92315 |
| 17 | 11:12:25.2 | -0.98007 | 83636 |
| 18 | 11:12:25.3 | -0.97399 | 70214 |
| 19 | 11:12:25.3 | -0.96792 | 54588 |
| 20 | 11:12:25.4 | -0.96185 | 39321 |
| 21 | 11:12:25.4 | -0.95578 | 26623 |
| 22 | 11:12:25.5 | -0.94970 | 16971 |
| 23 | 11:12:25.5 | -0.94363 | 10057 |
| 24 | 11:12:25.6 | -0.93756 | 5380 |
| 25 | 11:12:25.7 | -0.93148 | 2675 |
| 26 | 11:12:25.7 | -0.92541 | 1239 |
| 27 | 11:12:25.8 | -0.91934 | 607 |
| 28 | 11:12:25.8 | -0.91327 | 249 |
| 29 | 11:12:25.9 | -0.90719 | 72 |
+-----------+------------+------------+------------+
generator rel_scan ['2631ce75'] (scan num: 4)
center=-0.9921816, FWHM=0.0535156
Report the results#
Print a nice table with the results from each of our scans.
[24]:
tbl = pyRestTable.Table()
tbl.addLabel("scan ID")
tbl.addLabel("center")
tbl.addLabel("sigma")
for sid, cen, sigma in results:
tbl.addRow((sid, cen, sigma))
print(tbl)
======= =================== ====================
scan ID center sigma
======= =================== ====================
1 -1.0 0.1428571428571428
3 -0.9922136823434071 0.05668097101727676
4 -0.9921815602844319 0.053515606485192824
======= =================== ====================