{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monochromator Optimization\n",
    "\n",
    "In this notebook we demonstrate:\n",
    "\n",
    "* Using a custom plan that is tuned to a specific, common task at an XAFS (X-ray Absorption Fine Structure) beamline.\n",
    "* How this plan is used to take readings, do some prompt analysis on those readings, and immediately take an action based on the result of that analysis.\n",
    "\n",
    "## Science Background\n",
    "\n",
    "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). \n",
    "\n",
    "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)$\n",
    "\n",
    "* $d$ is the spacing between the lattice planes of the crystal\n",
    "* $\\Theta$ is the angle between the incident beam and the lattice planes of the first crystal\n",
    "\n",
    "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.\n",
    "\n",
    "![DCM](./static/doubleb.gif)  [(image source)](http://pd.chem.ucl.ac.uk/pdnn/inst2/condit.htm)\n",
    "\n",
    "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. \n",
    "\n",
    "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.\n",
    " \n",
    "## Running a beamline\n",
    "\n",
    "Bluesky provides all the tools we need to perform this optimization.  Here's the game plan:\n",
    "\n",
    "1. Perform a scan of the pitch motor over a range that will include the optimal pitch position\n",
    "2. Record an intensity signal from a detector downstream of the DCM\n",
    "3. When the scan is finished, analyze the intensity measurements and determine the position of peak intensity\n",
    "4. Move the second crystal pitch motor to its optimized position\n",
    "\n",
    "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:\n",
    "\n",
    "* The name of the motor to be scanned\n",
    "* The name of the detector to be monitored\n",
    "* The syntax of the mathematical tools used to analyze the measured data\n",
    "* The names of the standard plans for scanning and moving motors\n",
    "\n",
    "For the beamline staff, those things should be common knowledge.  For a general user or a new post-doc, those things are esoteric mysteries.\n",
    "\n",
    "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.\n",
    "\n",
    "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.\n",
    "\n",
    "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.\n",
    "\n",
    "For the sake of completeness, let's talk about the motor and detector being used in this example:\n",
    "\n",
    "* **motor**: The motor is a linear actuator pushing against a flexure controlling the pitch of the second monochrmoator crystal. It is controlled using a Delta-Tau controller provided by the instrument vendor. We talk to it through Ophyd's [EpicsMotor](https://blueskyproject.io/ophyd/generated/ophyd.epics_motor.html) interface.\n",
    "* **detector**: The detector in use is an ionization chamber.  At BMM, our ion chambers are almost always filled with nitrogen.  We read this signal using an NSLS-II-produced four-channel electrometer and Ophyd's [QuadEM](https://blueskyproject.io/ophyd/generated/ophyd.quadem.html) interface. In BMM's startup scripts, we rename the relevant signal to `I0`, which is a good, recognizable, semantic name for that signal.\n",
    "\n",
    "\n",
    "--------\n",
    "\n",
    "## Tutorial\n",
    "\n",
    "We start by setting up our Bluesky environment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'ipympl'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-4f0398f13e49>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'matplotlib'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'widget'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mbluesky\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mRunEngine\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mbluesky_tutorial_utils\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0msetup_data_saving\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_line_magic\u001b[0;34m(self, magic_name, line, _stack_depth)\u001b[0m\n\u001b[1;32m   2315\u001b[0m                 \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'local_ns'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getframe\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf_locals\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2316\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2317\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2318\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2319\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-108>\u001b[0m in \u001b[0;36mmatplotlib\u001b[0;34m(self, line)\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    185\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    186\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    189\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/core/magics/pylab.py\u001b[0m in \u001b[0;36mmatplotlib\u001b[0;34m(self, line)\u001b[0m\n\u001b[1;32m     97\u001b[0m             \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Available matplotlib backends: %s\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mbackends_list\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     98\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 99\u001b[0;31m             \u001b[0mgui\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbackend\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshell\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0menable_matplotlib\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgui\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlower\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgui\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgui\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    100\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_show_matplotlib_backend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgui\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbackend\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    101\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36menable_matplotlib\u001b[0;34m(self, gui)\u001b[0m\n\u001b[1;32m   3417\u001b[0m                 \u001b[0mgui\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbackend\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfind_gui_and_backend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpylab_gui_select\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3418\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3419\u001b[0;31m         \u001b[0mpt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mactivate_matplotlib\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbackend\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3420\u001b[0m         \u001b[0mpt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconfigure_inline_support\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbackend\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3421\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/core/pylabtools.py\u001b[0m in \u001b[0;36mactivate_matplotlib\u001b[0;34m(backend)\u001b[0m\n\u001b[1;32m    318\u001b[0m     \u001b[0;31m# when this function runs.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    319\u001b[0m     \u001b[0;31m# So avoid needing matplotlib attribute-lookup to access pyplot.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 320\u001b[0;31m     \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    321\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    322\u001b[0m     \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mswitch_backend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbackend\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/matplotlib/pyplot.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m   2347\u001b[0m     \u001b[0mdict\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__setitem__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrcParams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"backend\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrcsetup\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_auto_backend_sentinel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2348\u001b[0m \u001b[0;31m# Set up the backend.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2349\u001b[0;31m \u001b[0mswitch_backend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrcParams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"backend\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2350\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2351\u001b[0m \u001b[0;31m# Just to be safe.  Interactive mode can be turned on without\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/matplotlib/pyplot.py\u001b[0m in \u001b[0;36mswitch_backend\u001b[0;34m(newbackend)\u001b[0m\n\u001b[1;32m    219\u001b[0m         else \"matplotlib.backends.backend_{}\".format(newbackend.lower()))\n\u001b[1;32m    220\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 221\u001b[0;31m     \u001b[0mbackend_mod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mimportlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mimport_module\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbackend_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    222\u001b[0m     Backend = type(\n\u001b[1;32m    223\u001b[0m         \"Backend\", (matplotlib.backends._Backend,), vars(backend_mod))\n",
      "\u001b[0;32m/usr/lib/python3.8/importlib/__init__.py\u001b[0m in \u001b[0;36mimport_module\u001b[0;34m(name, package)\u001b[0m\n\u001b[1;32m    125\u001b[0m                 \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    126\u001b[0m             \u001b[0mlevel\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 127\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_bootstrap\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_gcd_import\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mlevel\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpackage\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlevel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    128\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    129\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'ipympl'"
     ]
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "import matplotlib.pyplot as plt\n",
    "from bluesky import RunEngine\n",
    "from bluesky_tutorial_utils import setup_data_saving\n",
    "\n",
    "\n",
    "RE = RunEngine()\n",
    "catalog = setup_data_saving(RE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will import a Bluesky *plan* from a script in the current directory, [plans.py](./plans.py). The plan operates on simulated hardware defined in another script, [simulated_hardware.py](./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](./plans.py) to understand how the rocking curve scan is implemented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from plans import rocking_curve\n",
    "\n",
    "help(rocking_curve)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RE(rocking_curve())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "Now let's look at a plot of the rocking curve data using matplotlib's `gcf` (i.e. [get current figure](https://matplotlib.org/3.2.1/api/_as_gen/matplotlib.pyplot.gcf.html)) method. At the beamline, this is typically displayed as a live plot during the rocking curve scan."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.gcf()  # Display a snapshot of the current state of the figure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "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:\n",
    "\n",
    "1. `peak`: find the position of maximum intensity\n",
    "2. `com`: find the position of the [center of mass](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.center_of_mass.html) of the masured peak\n",
    "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](https://lmfit.github.io/lmfit-py/builtin_models.html#skewedgaussianmodel) works well.\n",
    "\n",
    "Let's see how these work.\n",
    "\n",
    "-------\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run = catalog[-1]\n",
    "run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Access the saved data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = run.primary.read()\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot I0 vs pitch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()  # New figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.plot.scatter(\"pitch\", \"I0\")\n",
    "plt.gcf()  # Display a snapshot of the current state of the figure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could have gone straight to the plot in one line by chaining all of this together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog[-1].primary.read().plot.scatter(\"pitch\", \"I0\")\n",
    "plt.gcf()  # Display a snapshot of the current state of the figure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----\n",
    "\n",
    "The variable `data` contains the result of our just completed rocking curve measurement. The type of this data is an [xarray Dataset](http://xarray.pydata.org/en/stable/generated/xarray.Dataset.html).  In the following, we will work instead with a [pandas DataFrame](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html).  Here, we convert the xarray Dataset to a pandas DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = data.to_dataframe()\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "measurement = df.loc[:, ['I0', 'pitch']]\n",
    "measurement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Moving to the peak of the rocking curve\n",
    "\n",
    "First we find the time index of the data point which has the highest `I0` value, i.e. the peak of the plot above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i0max = measurement['I0'].idxmax()\n",
    "i0max"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we find the value of pitch at which the `I0` intensity was maximum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimal_pitch = measurement.loc[i0max, 'pitch']\n",
    "optimal_pitch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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:\n",
    "```\n",
    "RE(mv(pitch, optimal_pitch))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Moving to the center of mass of the rocking curve\n",
    "\n",
    "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.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.ndimage.measurements import center_of_mass\n",
    "import numpy\n",
    "arr = numpy.array(measurement['I0'])\n",
    "val = int(round(center_of_mass(arr)[0]))\n",
    "val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimal_pitch = measurement['pitch'].iloc[val]\n",
    "optimal_pitch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we want to move the pitch motor to its optimal value. So something like like:\n",
    "```\n",
    "RE(mv(pitch, optimal_pitch))\n",
    "```\n",
    "\n",
    "### Fitting a peak function to the rocking curve measurement\n",
    "\n",
    "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.\n",
    "\n",
    "Here is how we use [lmfit](https://lmfit.github.io/lmfit-py/) to fit a simple skewed Gaussian model to our rocking curve measurement.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lmfit.models import SkewedGaussianModel\n",
    "signal = numpy.array(measurement['I0'])\n",
    "position = numpy.array(measurement['pitch'])\n",
    "mod = SkewedGaussianModel()\n",
    "pars = mod.guess(signal, x=position)\n",
    "out = mod.fit(signal, pars, x=position)\n",
    "print(out.fit_report(min_correl=0))\n",
    "out.plot()\n",
    "optimal_pitch = out.params[\"center\"].value\n",
    "optimal_pitch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we move the pitch motor to its optimal value. Again:\n",
    "```\n",
    "RE(mv(pitch, optimal_pitch))\n",
    "```  \n",
    "We can examine the result of the fit by showing the plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.gcf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "-----\n",
    "\n",
    "## In conclusion\n",
    "\n",
    "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:\n",
    "```\n",
    "RE(rocking_curve())\n",
    "```\n",
    "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:\n",
    "```\n",
    "RE(rocking_curve(choice='peak'))\n",
    "```\n",
    "Here is how to select the center of mass calculation:\n",
    "```\n",
    "RE(rocking_curve(choice='com'))\n",
    "```\n",
    "and here is how to select the fitted interpretation:\n",
    "```\n",
    "RE(rocking_curve(choice='fit'))\n",
    "```\n",
    "\n",
    "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.\n",
    "\n",
    "### Hierarchies of plans in Bluesky\n",
    "\n",
    "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.\n",
    "\n",
    "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:\n",
    "```\n",
    "RE(change_edge('Fe'))\n",
    "```\n",
    "This plan \n",
    "\n",
    "1. verifies that the 14 individual axes that comprise the photon delivery system are in the correct positions for the specified energy range\n",
    "2. performs a rocking curve scan just like the one in this tutorial\n",
    "3. optimizes the position of the beam-definition slits\n",
    "4. makes sure detectors are properly configured to make measurements at the chosen absorption edge\n",
    "5. moves a rotation stage to select the correct reference foil from a collection of reference foilss\n",
    "\n",
    "leaving the beamline completely optimized and ready for measurement at the selected absorption edge.\n",
    "\n",
    "The reason such a magical plan is possible is because, in Bluesky, **plans are composed of plans**.  \n",
    "\n",
    "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()`.  \n",
    "\n",
    "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. \n",
    "\n",
    "In fact, at BMM, we allow *and* encourage our users to run the `change_edge()` command all by themselves!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}