5_mdx2_merging.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cdf2cde1-161d-48a9-b1e6-ab3802a46d2a",
   "metadata": {},
   "source": [
    "# Merging diffuse scattering with *mdx2*\n",
    "\n",
    "In this tutorial, we will:\n",
    "\n",
    "- Merge unscaled intensities for a quick look at reciprocal space\n",
    "- Fit a scaling model\n",
    "- Merge scaled intensities for a better quality map\n",
    "\n",
    "Requirements:\n",
    "- `geometry.nxs` from the *mdx2* geometry tutorial\n",
    "- `corrected.nxs` from the *mdx2* integration tutorial\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>NOTE: This tutorial assumes you have already procesed the full dataset (50 degrees).</b> If you've only processed the first 6 degrees so far (started from results_1_60.expt), you can use a pre-processed corrected.nxs (attached to this release on GitHub), which has the full 50 degrees.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab842cdb-9705-4905-a8ed-c7298d96009f",
   "metadata": {},
   "source": [
    "Imports (run this before you do anything else)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b09faee4-35ae-4bcc-92e6-4350f689e124",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for loading mdx2 objects in python\n",
    "from mdx2.utils import loadobj"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85698290-5cae-4437-a20d-475d12915813",
   "metadata": {},
   "source": [
    "### Background\n",
    "\n",
    "#### Merging\n",
    "\n",
    "The goal of merging is to combine redundant observations (equivalent by e.g. Laue group symmetry). Cubic insulin as 24 Laue group operators, so even with only 50 degrees of data we have good redundancy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7dc310e3-2dcd-4b90-8060-e713cc211fb0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading hkl_table from corrected.nxs\n",
      "  importing as HKLTable from mdx2.data\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD7CAYAAAB0d9PAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAU1UlEQVR4nO3dfbBdVXnH8e9DojQWiUm4ICbUy0jUgq0oaWCqrS+xSZQqtAPT1KmkTtpMLRbry9hQO5OKjQ12FGUsVCqRgC8Q0Q4opZgG0VoRuAgSQ2QSBSEDQjQR8Y0x+PSPvW45Od677j439y3J9zOz5+yz9l7PXee+nN9Ze+9zbmQmkiQN55DJHoAkaWozKCRJVQaFJKnKoJAkVRkUkqQqg0KSVDV9sgcw1o444ojs7++f7GFI0n7l9ttv/35m9g217YALiv7+fgYGBiZ7GJK0X4mI7w63zUNPkqQqg0KSVGVQSJKqDApJUpVBIUmqMigkSVUGhSSpyqCQJFUdcG+4U3v9q64bcZ/71p46ASORNJU5o5AkVRkUkqQqg0KSVGVQSJKqDApJUpVBIUmqMigkSVUGhSSpyjfcaZ/5xj3pwOaMQpJUZVBIkqoMCklSlUEhSaoyKCRJVQaFJKmqdVBExLSIuCMiPl/uz46IjRGxrdzO6tj33IjYHhH3RMSSjvaTImJz2XZhRERpPzQirirtt0REf0ef5eVrbIuI5WPyqCVJrfUyo3gLsLXj/ipgU2bOBzaV+0TE8cAy4ARgKXBRREwrfS4GVgLzy7K0tK8AdmfmccAFwPml1mxgNXAysBBY3RlIkqTx1yooImIecCrw0Y7m04D1ZX09cHpH+5WZ+Xhm3gtsBxZGxNHA4Zl5c2YmcHlXn8FaVwOLymxjCbAxM3dl5m5gI0+GiyRpArSdUXwQeCfwy462ozLzIYBye2Rpnws80LHfjtI2t6x3t+/VJzP3AI8Ccyq19hIRKyNiICIGdu7c2fIhSZLaGDEoIuIPgUcy8/aWNWOItqy0j7bPkw2Zl2Tmgsxc0NfX13KYkqQ22swoXgK8LiLuA64EXhkRHwceLoeTKLePlP13AMd09J8HPFja5w3RvlefiJgOzAR2VWpJkibIiEGRmedm5rzM7Kc5SX1jZv4ZcC0weBXScuCasn4tsKxcyXQszUnrW8vhqcci4pRy/uGsrj6Dtc4oXyOBG4DFETGrnMReXNokSRNkXz49di2wISJWAPcDZwJk5paI2ADcDewBzs7MJ0qfNwGXATOA68sCcClwRURsp5lJLCu1dkXEe4Dbyn7nZeaufRizJKlHPQVFZt4E3FTWfwAsGma/NcCaIdoHgBcM0f5zStAMsW0dsK6XcUqSxo7vzJYkVR0U/7jIf6wjSaPnjEKSVGVQSJKqDApJUpVBIUmqMigkSVUGhSSpyqCQJFUZFJKkKoNCklRlUEiSqgwKSVKVQSFJqjooPhTwQOQHHUqaKM4oJElVBoUkqcqgkCRVGRSSpCqDQpJUZVBIkqoMCklSlUEhSaoyKCRJVQaFJKnKoJAkVRkUkqQqg0KSVGVQSJKqDApJUpX/j0JTxkj/Y8P/ryFNDmcUkqQqg0KSVGVQSJKqDApJUpVBIUmqMigkSVUGhSSpasSgiIhfi4hbI+IbEbElIt5d2mdHxMaI2FZuZ3X0OTcitkfEPRGxpKP9pIjYXLZdGBFR2g+NiKtK+y0R0d/RZ3n5GtsiYvmYPnpJ0ojazCgeB16ZmS8ETgSWRsQpwCpgU2bOBzaV+0TE8cAy4ARgKXBRREwrtS4GVgLzy7K0tK8AdmfmccAFwPml1mxgNXAysBBY3RlIkqTxN2JQZOPH5e5TypLAacD60r4eOL2snwZcmZmPZ+a9wHZgYUQcDRyemTdnZgKXd/UZrHU1sKjMNpYAGzNzV2buBjbyZLhIkiZAq3MUETEtIu4EHqF54r4FOCozHwIot0eW3ecCD3R031Ha5pb17va9+mTmHuBRYE6lVvf4VkbEQEQM7Ny5s81DkiS11CooMvOJzDwRmEczO3hBZfcYqkSlfbR9Osd3SWYuyMwFfX19laFJknrV01VPmflD4Caawz8Pl8NJlNtHym47gGM6us0DHizt84Zo36tPREwHZgK7KrUkSROkzVVPfRHxjLI+A3gV8C3gWmDwKqTlwDVl/VpgWbmS6Viak9a3lsNTj0XEKeX8w1ldfQZrnQHcWM5j3AAsjohZ5ST24tImSZogbT5m/Ghgfbly6RBgQ2Z+PiJuBjZExArgfuBMgMzcEhEbgLuBPcDZmflEqfUm4DJgBnB9WQAuBa6IiO00M4llpdauiHgPcFvZ77zM3LUvD1iS1JsRgyIz7wJeNET7D4BFw/RZA6wZon0A+JXzG5n5c0rQDLFtHbBupHFKksaH78yWJFUZFJKkKv8Vaksj/ZtO8F91SjowOaOQJFUZFJKkKoNCklRlUEiSqgwKSVKVQSFJqjIoJElVBoUkqcqgkCRVGRSSpCqDQpJU5Wc96YDiZ3JJY88ZhSSpyqCQJFUZFJKkKoNCklRlUEiSqgwKSVKVl8dOMC/flLS/cUYhSaoyKCRJVQaFJKnKoJAkVRkUkqQqg0KSVGVQSJKqDApJUpVBIUmqMigkSVUGhSSpyqCQJFUZFJKkKoNCklRlUEiSqgwKSVLViEEREcdExBcjYmtEbImIt5T22RGxMSK2ldtZHX3OjYjtEXFPRCzpaD8pIjaXbRdGRJT2QyPiqtJ+S0T0d/RZXr7GtohYPqaPXpI0ojYzij3A2zPzN4FTgLMj4nhgFbApM+cDm8p9yrZlwAnAUuCiiJhWal0MrATml2VpaV8B7M7M44ALgPNLrdnAauBkYCGwujOQJEnjb8SgyMyHMvPrZf0xYCswFzgNWF92Ww+cXtZPA67MzMcz815gO7AwIo4GDs/MmzMzgcu7+gzWuhpYVGYbS4CNmbkrM3cDG3kyXCRJE6CncxTlkNCLgFuAozLzIWjCBDiy7DYXeKCj247SNresd7fv1Scz9wCPAnMqtbrHtTIiBiJiYOfOnb08JEnSCFoHRUQcBnwG+NvM/FFt1yHastI+2j5PNmRekpkLMnNBX19fZWiSpF61CoqIeApNSHwiMz9bmh8uh5Mot4+U9h3AMR3d5wEPlvZ5Q7Tv1ScipgMzgV2VWpKkCdLmqqcALgW2ZuYHOjZdCwxehbQcuKajfVm5kulYmpPWt5bDU49FxCml5lldfQZrnQHcWM5j3AAsjohZ5ST24tImSZog01vs8xLgDcDmiLiztP09sBbYEBErgPuBMwEyc0tEbADuprli6uzMfKL0exNwGTADuL4s0ATRFRGxnWYmsazU2hUR7wFuK/udl5m7RvdQJUmjMWJQZOZXGPpcAcCiYfqsAdYM0T4AvGCI9p9TgmaIbeuAdSONU5I0PnxntiSpyqCQJFUZFJKkKoNCklRlUEiSqgwKSVKVQSFJqjIoJElVBoUkqcqgkCRVGRSSpCqDQpJUZVBIkqoMCklSlUEhSaoyKCRJVQaFJKnKoJAkVRkUkqQqg0KSVGVQSJKqDApJUpVBIUmqmj7ZA5Cmov5V1424z31rT52AkUiTzxmFJKnKoJAkVRkUkqQqg0KSVGVQSJKqDApJUpVBIUmqMigkSVUGhSSpyqCQJFUZFJKkKoNCklRlUEiSqgwKSVLViEEREesi4pGI+GZH2+yI2BgR28rtrI5t50bE9oi4JyKWdLSfFBGby7YLIyJK+6ERcVVpvyUi+jv6LC9fY1tELB+zRy1Jaq3NjOIyYGlX2ypgU2bOBzaV+0TE8cAy4ITS56KImFb6XAysBOaXZbDmCmB3Zh4HXACcX2rNBlYDJwMLgdWdgSRJmhgjBkVmfhnY1dV8GrC+rK8HTu9ovzIzH8/Me4HtwMKIOBo4PDNvzswELu/qM1jramBRmW0sATZm5q7M3A1s5FcDS5I0zkZ7juKozHwIoNweWdrnAg907LejtM0t693te/XJzD3Ao8CcSq1fERErI2IgIgZ27tw5yockSRrKWJ/MjiHastI+2j57N2ZekpkLMnNBX19fq4FKktoZbVA8XA4nUW4fKe07gGM69psHPFja5w3RvlefiJgOzKQ51DVcLUnSBBptUFwLDF6FtBy4pqN9WbmS6Viak9a3lsNTj0XEKeX8w1ldfQZrnQHcWM5j3AAsjohZ5ST24tImSZpA00faISI+BbwcOCIidtBcibQW2BARK4D7gTMBMnNLRGwA7gb2AGdn5hOl1JtorqCaAVxfFoBLgSsiYjvNTGJZqbUrIt4D3Fb2Oy8zu0+qS1NW/6rrRtznvrWnTsBIpH0zYlBk5p8Os2nRMPuvAdYM0T4AvGCI9p9TgmaIbeuAdSONUZI0fnxntiSpyqCQJFUZFJKkKoNCklRlUEiSqgwKSVKVQSFJqjIoJElVBoUkqcqgkCRVGRSSpCqDQpJUZVBIkqoMCklSlUEhSaoyKCRJVQaFJKnKoJAkVRkUkqQqg0KSVGVQSJKqDApJUtX0yR6ApLr+VdeNuM99a0+dgJHoYOWMQpJUZVBIkqoMCklSlUEhSaoyKCRJVQaFJKnKoJAkVRkUkqQq33AnHSRGeuOeb9rTcJxRSJKqDApJUpVBIUmq8hyFpNb8gMKDkzMKSVKVMwpJE86Zyf5lvwiKiFgKfAiYBnw0M9dO8pAkTQFe8jsxpnxQRMQ04F+BPwB2ALdFxLWZeffkjkzSgcDZzcimfFAAC4HtmfkdgIi4EjgNMCgkTRljEThjFVpjHX6Rma13ngwRcQawNDP/otx/A3ByZr65Y5+VwMpy93nAPSOUPQL4/hgMbyrVcSzjW8exjG+dqTSWsaqzv43l2ZnZN9SG/WFGEUO07ZVumXkJcEnrghEDmblgnwc2heo4lvGt41jGt85UGstY1TmQxrI/XB67Azim4/484MFJGoskHXT2h6C4DZgfEcdGxFOBZcC1kzwmSTpoTPlDT5m5JyLeDNxAc3nsuszcso9lWx+m2o/qOJbxreNYxrfOVBrLWNU5YMYy5U9mS5Im1/5w6EmSNIkMCklSlUEhSaoyKHoQEc+PiEURcVhX+9IeaiyMiN8p68dHxNsi4jVjMLbLx6DGS8t4FvfQ5+SIOLysz4iId0fE5yLi/IiY2UOdcyLimJH3rNZ4akScFRGvKvdfHxEfjoizI+IpPdZ6TkS8IyI+FBHvj4i/6uXxSAeSg/pkdkS8MTM/1nLfc4Czga3AicBbMvOasu3rmfniFjVWA6+mudpsI3AycBPwKuCGzFzTcizdlwcH8ArgRoDMfF3LOrdm5sKy/pc0j+8/gMXA59p8+GJEbAFeWK5OuwT4KXA1sKi0/3HLsTwK/AT4NvAp4NOZubNN344an6D53j4N+CFwGPDZMpbIzOUt65wDvBb4EvAa4E5gN/BHwF9n5k29jEsji4gjM/ORyR7HoIiYk5k/mOxxTBmZedAuwP097LsZOKys9wMDNGEBcEcPNabRPJH9CDi8tM8A7uphLF8HPg68HHhZuX2orL+shzp3dKzfBvSV9V8HNressbVzXF3b7uxlLDQz3MXApcBO4L+A5cDTW9a4q9xOBx4GppX70eP3d3NH36cBN5X132j7sy77zwTWAt8CflCWraXtGWP0O3x9y/0OB/4ZuAJ4fde2i3r4es8ELqb5oM45wD+W79cG4OiWNWZ3LXOA+4BZwOwexrK063t9KXAX8EngqB7qrAWOKOsLgO8A24Hvtv17Kn+T/wA8Zx9/nguAL5a/72NoXlA+Wv4+X9RDncOA84Atpf9O4GvAn49mXAf8oaeIuGuYZTNwVA+lpmXmjwEy8z6aJ+dXR8QHGPpjRoayJzOfyMyfAt/OzB+Vej8DftnDWBYAtwPvAh7N5hXuzzLzS5n5pR7qHBIRsyJiDs0r7p1lPD8B9rSs8c2IeGNZ/0ZELACIiOcCv+hhLJmZv8zML2TmCuBZwEXAUpo/3DYOKW/KfDrNE/zgoaJDgZ4OPfHke4wOLfXIzPt7rLOBZiby8syck5lzaGZ+u4FPty0SES8eZjmJZnbbxsdofk8/AyyLiM9ExKFl2yltxwJcRvOBnA/QPKH9DDgV+B/g31rW+D7N7+/gMgDMpXmyHehhLO/tWH8/zYul19I8qX6khzqnZubg5yD9C/AnmXkczSdWv79ljVnAM4AvRsStEfHWiHhWD2MYdBHwPuA64KvARzJzJrCqbGvrEzR/N0uAdwMXAm8AXhER7611HNK+pN/+sNC8sjwReHbX0g882EOdG4ETu9qmA5cDT7SscQvwtLJ+SEf7TLpejbesN4/mCefD9DA76uh/X/llurfcPjOffDVyZ8saM2mePL5dHt8vSq0v0Rx6ajuWOyrbZrSs8dbytb8LnANsAv6d5hXv6h7G8haaV6aX0MwG3lja+4Av91DnntFsG2LfJ8rv3xeHWH7WssadXfffBfwvzav51r977D0Lvb/2NSo13kEzW/ytjrZ7R/H7+/XhvnbbsZR9vwVML+tf69rWdmbdOZbfo3lS/175Ga0co+/vHT3U+UbX/dvK7SHAt3r+XvfaYX9baKajLx1m2yd7qDNv8Il0iG0vaVnj0GHaj+j8oxnFYzwVeO8Yfs+eBhzbY5+nAy8ETqKHaX9H/+eO0difBTyrrD8DOANYOIo6J5S+z9+HsXwBeGfn94NmFvt3wH/3UOebwPxhtj3QssZWOl6clLblNIcmvtvDWL7Rsf5PXdtaPamWfQdf5Hyg/O58ZxTf3x3A24C307xAiI5tvRxq/Jvys3olzaG0DwK/T/NK/IqWNX4lbGkOMy8FPtbDWG6mOfx6Js0LntNL+8uAgR7qfHXweY9mlnVDx7bWL1L+v0+vHVxcXNotNIcjzqd5xbqrLFtL26we6pwBPG+Ybae3rPE+4FVDtC8FtvUwlvMo5+q62o8Drh7F9+i1NMfOvzeKvqu7lsFzbM8ELu+x1suBq2jOlW0G/pPmXxdMb9n/yjH6nXkhzccVXQ88n+Y/e/6wBPrv9lDnt4FbS9+vUF6I0cyKz+l1XAf1VU/SZOnlirvxrjPZY4mIGTQngb852WMZjzpTaSyjrWNQSJMgIu7PzN+YCnUcy/jWmUpjGW2dKf/psdL+KiLuGm4TPVxxNxZ1HMv41plKYxnLOoMMCmn8HEVzeeLurvagOdk4kXUcy/jWmUpjGcs6gEEhjafP05z4vbN7Q0TcNMF1HMv41plKYxnLOk0fz1FIkmoO+HdmS5L2jUEhSaoyKCRJVQaFJKnKoJAkVf0fZ5jkiOPi7S4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# load table\n",
    "tab = loadobj('corrected.nxs','hkl_table')\n",
    "\n",
    "# compute number of times each h,k,l is observed\n",
    "c = tab.to_frame().groupby(['h','k','l'])['phi'].count().value_counts().sort_index()\n",
    "c[c>1000].plot.bar();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ceb11242-eb9a-49a3-98aa-4492caae6b4d",
   "metadata": {},
   "source": [
    "#### Scaling\n",
    "\n",
    "In order to merge a high quality map, it's usually necessary to correct for non-ideal aspects of the experiment by scaling.\n",
    "\n",
    "Scaling is important in Bragg data reduction also. Programs like *aimless*, *xscale*, and *dials.scale* are commonly used, and they have similar scaling models.\n",
    "\n",
    "One key difference between diffuse and Bragg data is that Bragg intensities have an internal background reference (the intensity near the peak), so separate background measurements are not required. In diffuse scattering, the experimental background subtraction is not always perfect. So a scaling model for diffuse scattering might also include a correction to the background.\n",
    "\n",
    "We previously implemented a scaling model for diffuse scattering in *mdx-lib* that includes four terms:\n",
    "\n",
    "- $a =$ absorption correction (a function of phi, x, and y) \n",
    "- $b =$ overall scale (a function of phi) \n",
    "- $c =$ isotropic background (a function of phi and s) \n",
    "- $d =$ detector chip gain (a function of x and y)\n",
    "\n",
    "The model is refined using regularization to ensure that the correction factors vary smoothly (e.g. as a function of phi).\n",
    "\n",
    "In *mdx2*, we have only implemented the $b$ term so far. That turns out to be good enough for the insulin dataset, so we'll use it here. The procedure for refining the full scaling model is similar, and will be available in future releases.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4f5e289-1a30-4df7-ae14-fd7bef01fe49",
   "metadata": {},
   "source": [
    "## 1) Quick slice\n",
    "\n",
    "### *mdx2.merge*\n",
    "\n",
    "This command lets us merge all redundant observations to produce a table of intensities in the asymmetric unit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "17fa6043-20c7-4742-90cf-6080f8502cb8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: mdx2.merge [-h] [--scale SCALE] [--outfile OUTFILE] hkl\n",
      "\n",
      "Merge corrected intensities using a scaling model\n",
      "\n",
      "positional arguments:\n",
      "  hkl                NeXus file with hkl_table\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help         show this help message and exit\n",
      "  --scale SCALE      NeXus file with scaling_model (default: None)\n",
      "  --outfile OUTFILE  name of the output NeXus file (default: merged.nxs)\n"
     ]
    }
   ],
   "source": [
    "!mdx2.merge --help"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cace334e-ee76-4855-90ae-761e0c51b565",
   "metadata": {},
   "source": [
    "#### Parameter choice\n",
    "\n",
    "The scaling model is optional. Since we haven't scaled the data let, let's just run it on the unscaled data and see what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8fc3af8f-fc1f-48de-b74b-78adad1c135f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading hkl_table from corrected.nxs\n",
      "  importing as HKLTable from mdx2.data\n",
      "Grouping redundant observations\n",
      "Merging\n",
      "Exporting <class 'mdx2.data.HKLTable'> to nexus object\n",
      "  writing hkl_table to merged_not_scaled.nxs\n",
      "done!\n"
     ]
    }
   ],
   "source": [
    "!mdx2.merge corrected.nxs --outfile merged_not_scaled.nxs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "640304cf-ea4c-4874-8506-5a16cebbd872",
   "metadata": {},
   "source": [
    "### *mdx2.map*\n",
    "\n",
    "The result from merging is just another table of values. To visualize this table we'll generate a slice through reciprocal space by symmetry expansion and plot it using *nexpy*.\n",
    "\n",
    "The *mdx.map* command can make slices for plotting (it can also generate a 3D volume - more on this later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "45b077a5-ffc5-4605-8216-ea5f12089a48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: mdx2.map [-h] [--symmetry TF] [--limits H1 H2 K1 K2 L1 L2]\n",
      "                [--signal SIGNAL] [--outfile OUTFILE]\n",
      "                geom hkl\n",
      "\n",
      "Create a map from data in an hkl table\n",
      "\n",
      "positional arguments:\n",
      "  geom                  NeXus file with symmetry and crystal\n",
      "  hkl                   NeXus file with hkl_table\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  --symmetry TF         apply symmetry operators? (default: True)\n",
      "  --limits H1 H2 K1 K2 L1 L2\n",
      "                        region to map (default: [0, 10, 0, 10, 0, 10])\n",
      "  --signal SIGNAL       column in hkl_table to map (default: intensity)\n",
      "  --outfile OUTFILE     name of the output NeXus file (default: map.nxs)\n"
     ]
    }
   ],
   "source": [
    "!mdx2.map --help"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d12d21e-3213-4af7-9c7e-b183a6f02476",
   "metadata": {},
   "source": [
    "#### Parameter choice\n",
    "\n",
    "We need to set the limits of the map. We'll look at a slice with $l=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f3b8de3d-e75f-4550-9f2f-32c764684099",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading hkl_table from merged_not_scaled.nxs\n",
      "  importing as HKLTable from mdx2.data\n",
      "Reading crystal from geometry.nxs\n",
      "  importing as Crystal from mdx2.geometry\n",
      "Reading symmetry from geometry.nxs\n",
      "  importing as Symmetry from mdx2.geometry\n",
      "map region:\n",
      "  h from -50.0 to 50.0 (401 grid points)\n",
      "  k from -50.0 to 50.0 (401 grid points)\n",
      "  l from 0.0 to 0.0 (1 grid points)\n",
      "generating Miller index array\n",
      "mapping Miller indices to asymmetric unit\n",
      "looking up intensity in data table\n",
      "preparing output array\n",
      "Exporting <class 'mdx2.geometry.GridData'> to nexus object\n",
      "  writing intensity to slice_not_scaled.nxs\n",
      "done!\n"
     ]
    }
   ],
   "source": [
    "!mdx2.map geometry.nxs merged_not_scaled.nxs --limits -50 50 -50 50 0 0 --outfile slice_not_scaled.nxs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11a618b2-1bd4-4f0b-b063-84f19d5c3949",
   "metadata": {},
   "source": [
    "### Visualization in *nexpy*\n",
    "\n",
    "Open `slice_not_scaled.nxs` in *nexpy*. Adjust the sliders in the \"signal\" tab to clearly see the background.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>Question:</b> Can you see any artifacts in the map?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84eab100-7824-4cf9-ab53-7182d98a8346",
   "metadata": {},
   "source": [
    "## 2) Scaling\n",
    "\n",
    "### *mdx2.scale*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "e9319ea3-1e73-44af-a89d-dc5b9191a8be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: mdx2.scale [-h] [--smoothness SMOOTHNESS] [--phi_increment DEGREES]\n",
      "                  [--iter ITER] [--outfile OUTFILE]\n",
      "                  hkl\n",
      "\n",
      "Fit scaling model to unmerged corrected intensities\n",
      "\n",
      "positional arguments:\n",
      "  hkl                   NeXus file with hkl_table\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  --smoothness SMOOTHNESS\n",
      "                        amount to multiply the regularization parameter\n",
      "                        (default: 1)\n",
      "  --phi_increment DEGREES\n",
      "                        spacing of phi control points in degrees (default: 1)\n",
      "  --iter ITER           number of iterations (default: 5)\n",
      "  --outfile OUTFILE     name of the output NeXus file (default: scales.nxs)\n"
     ]
    }
   ],
   "source": [
    "!mdx2.scale --help"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29075362-b485-4b60-a0fc-718ad47da6e5",
   "metadata": {},
   "source": [
    "#### Parameter choice\n",
    "\n",
    "The defaults are fine to get started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "af7fee5b-6162-486d-945c-6a303690450a",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading hkl_table from corrected.nxs\n",
      "  importing as HKLTable from mdx2.data\n",
      "Grouping redundant observations\n",
      "initializing scaling model with 51 samples\n",
      "iteration 1 of 5\n",
      "  re-calculating scales\n",
      "  merging\n",
      "  fitting the model\n",
      "  current x2: 1.3362474224487717\n",
      "iteration 2 of 5\n",
      "  re-calculating scales\n",
      "  merging\n",
      "  fitting the model\n",
      "  current x2: 1.0108451187022505\n",
      "iteration 3 of 5\n",
      "  re-calculating scales\n",
      "  merging\n",
      "  fitting the model\n",
      "  current x2: 0.9876241832268707\n",
      "iteration 4 of 5\n",
      "  re-calculating scales\n",
      "  merging\n",
      "  fitting the model\n",
      "  current x2: 0.9859305872556283\n",
      "iteration 5 of 5\n",
      "  re-calculating scales\n",
      "  merging\n",
      "  fitting the model\n",
      "  current x2: 0.9858042534411342\n",
      "finished refining\n",
      "Exporting <class 'mdx2.scaling.ScalingModel'> to nexus object\n",
      "  writing scaling_model to scales.nxs\n",
      "done!\n"
     ]
    }
   ],
   "source": [
    "!mdx2.scale corrected.nxs --smoothness 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b9bca26-3e31-45d7-8ebb-9d7f6a9fc0f1",
   "metadata": {},
   "source": [
    "### Results\n",
    "\n",
    "You should see $\\chi^2$ (`x2`) drop with each iteration, until it settles down to ~1 by iteration 5."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a73bb0df-3817-4c35-9b7c-a54a85578a7a",
   "metadata": {},
   "source": [
    "The scale model parameters are saved in a NeXus file. You can plot them by opening the file in *nexpy* and double-clicking. For convenience I'll plot them in the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "9ce57650-c0da-4626-b885-cd8188553577",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading scaling_model from scales.nxs\n",
      "  importing as ScalingModel from mdx2.scaling\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcVklEQVR4nO3df7wcdX3v8debEODwQ05CDkoOCeFSDMXy0wP0IVgBuQa4UiLWKvoQFSyX3mLFq1zAtoDFRxPEH20f1lKuIgUkqBfMxSIiV4p5WPzBiUHAAooiTU6QBEhAS8QQPvePmWM2y85mN2dnZ3bm/Xw89nFmZ2Z3vt89u/OZ789RRGBmZvW1XdEJMDOzYjkQmJnVnAOBmVnNORCYmdWcA4GZWc05EJiZ1ZwDgU2ZpEslXZ8uz5X0K0nTik7XtpB0rKRVHe7723z3g6S7JL03XX6HpG/069i91s1n15hvy4cDgfVURPxHROwaEZuKTksRJL1d0g15HycivhARb8j7OFYPDgRmvXUy8LWiE2HWDQeCmpJ0gaQJSb+U9LCk10uaJunDkn6arl8uaU66/99JWinp2XT9azPed56kkLR9+vwuSZdJ+rf0Pb8haVbD/mdIekzSU5L+StLPJZ2wlbRfKunLkq5P3/N+Sa+UdJGkNWk639Cw/2xJt0h6WtIjkv6kYduQpGskrZP078ARTceaLekmSWslPSrpz9ukazvgvwJfl7RTmr6nJK2XdI+kl6f7zZT0eUmr0+MuTdfPkPQv6bHWpct7Zxzr3ZK+3fA8JJ0j6Sfpa/9BktJt0yR9QtKTaR7ObfwftcnPXZI+KunutLrvq5L2kPSF9Htwj6R5Dfu/Jl33TPr3NQ3b9pX0rfT/dQcwq+lYv58eZ72kH0o6tl3arLccCGpI0nzgXOCIiNgNWAD8HPifwOkkV7UvA84Enktfdg9wKDATuAH4sqSdOjzk24H3AHsCOwAfStNxIPAZ4B3AXsDuwGiH73kKcB0wA1gB3E7yfR4F/hr4p4Z9lwCrgNnAHwF/I+n16bZLgP3SxwLgXZMvSk/sXwV+mL7v64HzJC3ISNORwM8i4sn0fXYH5gB7AOcAG9L9rgN2Bl6VfiafStdvB3we2AeYm+7/6Q4/D4A3kgSyQ4A/TvMD8CfASST/v8OBhV2859uAd5Lkfz/gO2kaZwIPknx+SJoJ3Ar8PUl+PwncKmmP9H1uAJaTBIDL2PJzHk1f+9H0fT8E3CRppIt02lREhB81ewC/A6wBTgCmN6x/GDi1w/dYBxySLl8KXJ8uzwMC2D59fhfwlw2v+x/A19Pli4ElDdt2Bn4DnLCVY18K3NHw/BTgV8C09PluaRqGSU7Em4DdGvZfBFyTLv8MOLFh29nAqnT5KOA/mo59EfD55nynzy8D/ipdPhO4Gzi46fV7AS8CMzr4jA8F1jU8vwt4b7r8buDbDdsCOKbh+ZeAC9PlO4H/3rDthMb/UZvj3wX8RcPzTwC3NX3u96bL7wS+3/T676TpnAu8AOzSsO2Ghu/MBcB1Ta+9HXhXc779yOfhEkENRcQjwHkkJ7I1km6UNJvkpPnTVq+R9EFJD6bF/vUkV7uzWu3bwi8alp8Ddk2XZwMrG9L1HPBUh+/5RMPyBuDJ2NxAPXnlvWt6jKcj4pcN+z/G5pLHFmlIt03aB5idVlesT/P9YeDlGWlqbB+4juRkdmNaBfQxSdNJPuOnI2Jd84sl7Szpn9KqsmeBZcCwOu+B1dHn3LS8Nc2fc/PzxmM0fnaw+XOeTRLQ/rNp26R9gLc0fc7HkARN6wMHgpqKiBsi4hiSH2EAl5OcIPZr3ldJe8AFJNUNMyJiGHgG0BST8Tjw2zpwSUMk1Qq9tBqYKWm3hnVzgYmGNMxp2jZpJfBoRAw3PHaLiJObDyLpFSQnrh8ARMTGiPhIRBwIvIak2uaM9D1nShpukdYPAvOBoyLiZcAfTL59Vzl+qS0+Z7bMb6+sJvkuNZr8nB8HZkjapWnbpJUkJYLGz3mXiFicQzqtBQeCGpI0X9LxknYEfk1yZbcJ+CxwmaT9lTg4rePdjaRovxbYXtLFJG0IU/V/gFPSRsYdgI8w9ZPeFiJiJUkVzaK0Afdg4CzgC+kuXwIuShtq9wbe1/Dy7wPPKmlYH0obXX9P0hYNyqmTSaq8AkDScZIOSq/mnwU2Apsi4nHgNuAz6TGnS5o84e9G8r9Yn9a5X9Kjj+FLwPsljaYB6IIevW+jrwGvVNJ9dntJbwUOBP4lIh4DxoGPSNpB0jEk1UqTrif5HixIP+OdlIznaNlQbr3nQFBPOwKLgSdJqhP2JKny+CTJSeMbJCevzwFDJFUctwE/JinS/5ruqhdaiogfkZx4byS5avwlSdvF81N97yank7RdrAa+AlwSEXek2z5CkqdHSfJ9XUP6NpGcsA5Ntz9JEix3b3GM5m6jryAJdM+SNKp+i+SEB0l9+kbgIZL8npeu/1uSz/tJ4LvA17chr638b5K83UfSsP41ksDes7EeEfEUSanngyTVe/8LeGMkDeeQdBg4CniaJMBd2/DalcCpJN/BtSTfrfPx+alvlF7AmBVO0q7AemD/iHi04OR0LO2G+Qtgv4h4puj0bI2kk4ArI6K5KsdqyhHXCiXplLSRdBfg48D9JF1ZB8lMkt5CpQwCabXWyWmVzSjJFflXik6XlYcDgRXtVJIqm9XA/sDbIiIk3ZYOYmp+fLjY5L5URKyJiH8sOh1tiKQKbB1J1dCDJF13yfiMf6WMAYNWTa4aMjOrOZcIzMxqru1cI2U0a9asmDdvXtHJMDMbKMuXL38yIlpO2zFwgWDevHmMj48XnQwzs4EiqXnk92/lVjUk6WolM0E+sJX9jpC0SdIf5ZUWMzPLlmcbwTXAie12SEddXk4yYMnMzAqQWyCIiGUkowjbeR9wE8noSjMzK0BhvYbSgS1vAq7sYN+zJY1LGl+7dm3+iTMzq5Eiu4/+LXBBdHBv24i4KiLGImJsZMT3qjAz66Uiew2NkczVDsm89idLeiEilvb6QEtXTHDF7Q+zev0GZg8Pcf6C+Sw8rNMbYZmZVVthgSAi9p1clnQNyXS1S3t9nKUrJrjo5vvZsDEpeEys38BFN98P4GBgZka+3UeXkNyqbr6kVZLOUnJz7XPyOmYrV9z+8G+DwKQNGzdxxe0P9zMZZmallVuJICJO72Lfd+eVjtXrN3S13sysbio/19Ds4aGu1puZ1U3lA8H5C+YzNH3Le38PTZ/G+QvmF5QiM7NyGbi5hro12SDsXkNmZq1VPhBAEgx84jcza63yVUNmZtaeA4GZWc05EJiZ1ZwDgZlZzTkQmJnVnAOBmVnN1aL7aDuemdTM6q7WgcAzk5qZ1bxqyDOTmpnVPBB4ZlIzs5oHAs9MamZW80DgmUnNzGreWOyZSc3Mah4IwDOTmpnVumrIzMwcCMzMas+BwMys5hwIzMxqzoHAzKzmHAjMzGrOgcDMrOYcCMzMas6BwMys5mo/sjiLb1hjZnXhQNCCb1hjZnWSW9WQpKslrZH0QMb2UyXdJ+leSeOSjskrLd3yDWvMrE7ybCO4BjixzfZvAodExKHAmcBnc0xLV3zDGjOrk9wCQUQsA55us/1XERHp012AyNq333zDGjOrk0J7DUl6k6SHgFtJSgVZ+52dVh+Nr127Nvd0+YY1ZlYnhQaCiPhKRBwALAQua7PfVRExFhFjIyMjuadr4WGjLDrtIEaHhxAwOjzEotMOckOxmVVSKXoNRcQySftJmhURTxadHvANa8ysPgorEUj6HUlKlw8HdgCeKio9ZmZ1lVuJQNIS4FhglqRVwCXAdICIuBJ4M3CGpI3ABuCtDY3HZmbWJ7kFgog4fSvbLwcuz+v4ZmbWGc81ZGZWc6VoLDYzs2x5z33mQGBmVmL9mPvMVUNmZiXWj7nPXCLokqenNrN+6sfcZy4RdGGyiDaxfgPB5iLa0hUTRSfNzCqqH3OfORB0wdNTm1m/9WPuM1cNdcHTU5tZv01WPbvXUEnMHh5iosVJ39NTm1me8p77zFVDXfD01GZWRS4RdKEfRTQzs35zIOiSp6c2s6px1ZCZWc05EJiZ1ZwDgZlZzTkQmJnVnAOBmVnNuddQj3gyOjMbVA4EPdCP+cLNzPLiqqEe8GR0ZjbIHAh6wJPRmdkgcyDogX7MF25mlhcHgh7wZHRmNsjcWNwDnozOzAaZA0GPeDI6MxtUrhoyM6s5BwIzs5pz1VDOPOLYzMrOgSBHHnFsZoMgt6ohSVdLWiPpgYzt75B0X/q4W9IheaWlKB5xbGaDIM8SwTXAp4FrM7Y/CrwuItZJOgm4Cjgqx/T0nUccm1k3iqpKzq1EEBHLgKfbbL87ItalT78L7J1XWoriEcdm1qnJquSJ9RsINlclL10xkfuxy9Jr6CzgtqyNks6WNC5pfO3atX1M1tR4xLGZdarIquTCG4slHUcSCI7J2iciriKpOmJsbCz6lLQp84hjM+tUkVXJhQYCSQcDnwVOioinikxLXjzi2Mw6MXt4iIkWJ/1+VCUXVjUkaS5wM/DOiPhxUekwMyuDIquScysRSFoCHAvMkrQKuASYDhARVwIXA3sAn5EE8EJEjOWVHjOzMiuyKlkRA1PlDiRtBOPj40Unw8xsoEhannWxXXhjcV156gkzKwsHggJ46gkzK5OyjCOoFU89YWZl4kBQAE89YWZl4kBQAE89YWZl4kBQAE89YWZl4sbiAnjqCTMrEweCgnjqCTMrC1cNmZnVnAOBmVnNORCYmdWc2whKxlNPmFm/ORCUiKeeMLMiuGqoRDz1hJkVwYGgRDz1hJkVwYGgRDz1hJkVwYGgRDz1hJkVwY3FJeKpJ8ysCA4EJeOpJ8ys3xwIzMz6rGzjhRwIzMz6qIzjhToKBJIubrU+Iv66t8kxM6u2duOFSh0IgP9sWN4JeCPwYO+TY2ZWbWUcL9RRIIiITzQ+l/Rx4JZcUmQtla1O0cy2zezhISZanPSLHC+0reMIdgb+Sy8TYtkm6xQn1m8g2FynuHTFRNFJM7MulXG8UKdtBPcDkT6dBowAbh/okzLWKZrZtinjeKFO2wje2LD8AvBERLyQQ3qshTLWKZrZtivbeKFO2wgeyzshlq2MdYpmVh25zTUk6WpJayQ9kLH9AEnfkfS8pA/llY4qKGOdoplVR56Tzl0DnNhm+9PAnwMfzzENlbDwsFEWnXYQo8NDCBgdHmLRaQeVqmhpZoMrt5HFEbFM0rw229cAayT9t7zSUCVlq1M0s+rwFBMDzuMLzGyqBiIQSDobOBtg7ty5BaemPMo4Z4mZDZ6BuDFNRFwVEWMRMTYyMlJ0ckrD9zg2s14YiEBgrXl8gZn1Qm5VQ5KWAMcCsyStAi4BpgNExJWSXgGMAy8DXpR0HnBgRDybV5qqxuMLzKwX8uw1dPpWtv8C2Duv49fB+Qvmb9FGAB5fYGbdG4jGYmutjHOWmNngcSAYcB5fYGZT5UBQUR5fYGadciCoII8vMCuHQbkgc/fRCvL4ArPiDdINpRwIKsjjC8yKN0gXZA4EFZQ1jsDjC8z6Z5AuyBwIKsj3LzAr3iBdkDkQVFC7+xcsXTHB0YvvZN8Lb+XoxXeWsr7SrAoG6YLMvYYqqtX4AvcmMuufQRrw6UBQI+0ar8r45TQbdIMy4NNVQzUySI1XZtY/DgQ1MkiNV2bWPw4ENTJIjVdm1j9uI6iRQWq8MrP+cSComazGq0GZE8XMes+BwNyt1Kzm3EZgAzUnipn1nksE5m6lZlNQhWpVlwjM3UrNttEgTTXdjgOBuVup2TaqSrWqq4bM3UrNtlFVqlUdCAwYnDlRzMpk9vAQEy1O+oNWreqqIdsqT11t1lpVqlVdIrC2PMbALFtVqlUdCKwtT11t1l4VqlVdNWRtVaUxzMyyORBYWx5jYFZ9DgTWVlUaw8ymqsqdJnILBJKulrRG0gMZ2yXp7yU9Iuk+SYfnlRbbdgsPG2XRaQcxOjyEgNHhIRaddtDA14madaMqI4iz5NlYfA3waeDajO0nAfunj6OAf0z/WslUoTHMbCqq3mkit0AQEcskzWuzy6nAtRERwHclDUvaKyIezytN1ltVmGzLrBNV7zRRZBvBKLCy4fmqdN1LSDpb0rik8bVr1/YlcdZe1YvKZo2q3mmiyECgFuui1Y4RcVVEjEXE2MjISM7Jsk5UZbIts05UvdNEkQPKVgFzGp7vDawuKC3WpaoXlc0aVWUEcZYiA8EtwLmSbiRpJH7G7QODoyqTbZk1y2r7qnKnidwCgaQlwLHALEmrgEuA6QARcSXwNeBk4BHgOeA9eaXFeu/8BfO3mIMINheV3Yhsg6quc2vl2Wvo9K1sD+DP8jq+5SurqAzU8odk1VD1bqJZPOmcbbNWReWjF99Zyx+SVUNd2748xYT1VF1/SFYNVe8mmsWBwHqqrj8kq4aqdxPN4qoh66l2jchmZdKuU0PdOjs4EFhPtfshuTeRlcXWegfV7XvpQGA91+qHVNdueVZOde0dlMVtBNYXnpLCysSdGrbkQGB94R+elYk7NWzJgcD6wj88K5O69g7K4kBgfeEfnpWJ77y3JTcWW1+4N5EVpY6TyHXLgcD6xr2JrN/8/eqMq4asUO5NZHny96szDgRWKPcmsjz5+9UZVw1Zodrd4MZtBzZVvoFSZ1wisEJl9SY67oARLrr5fibWbyDYXLe7dMVEMQm1geTeap1xicAKldWbyFMAWLc8idy2U3KjsMExNjYW4+PjRSfDcrbvhbfS6psp4FNvPdQ/bNtCc+8gSK786zw2oJmk5REx1mqbq4aslLLqcHcfmu4qI3sJ9w6aGgcCK6Wsul0J/+DtJdw7aGocCKyUsqYAWP/cxpb7+wdfb57LamrcWGyl1Wok8hW3P+zupjXX6v/sO+NNjRuLbaBkNQq++dWj3LR8omVjIbjXSFW0axQG/5/baddY7EBgA6fVFWFWSWF4aDrPv/Cie5NUxNGL72z5fx4dHuLfLjy+gBQNjnaBwFVDNnBaVRl94Iv3ttx3/YaXtil4PMJgaBXw3SicDwcCq4SsqQSyrF6/wW0KJZY1a+jwztNZ16LDgBuFp8a9hqwSsrqbzth5esv9PR6h3LLGBUTgKSNy4BKBVULWVAJAy8bFrY1HcEmhWFlVPc9s2OiR5TlwILDKaHfHqeYTR1abwmTJwDcy6Z9WVXTtZg31ncV6L9deQ5JOBP4OmAZ8NiIWN22fAVwN7Af8GjgzIh5o957uNWS9kNX7ZJrEpha/CfdKyce2dAd2ENg2hcw1JGka8A/AScCBwOmSDmza7cPAvRFxMHAGSdAwy11Wm0KrIACbqyqWrpjg6MV3su+Ft3L04jvdpjBFWW0B//rQWt9cvo/yrBo6EngkIn4GIOlG4FTg3xv2ORBYBBARD0maJ+nlEfFEjukyazv9dbuRy6426q123UFdBdQ/eQaCUWBlw/NVwFFN+/wQOA34tqQjgX2AvYEtAoGks4GzAebOnZtXeq1msk40WVMVbG2GSzdgttdtW4D1T57dR9ViXXO5ezEwQ9K9wPuAFcALL3lRxFURMRYRYyMjIz1PqNmkrMnuFh42mnn1OlkycFfUbJOlqebP6LgDRtwdtATyLBGsAuY0PN8bWN24Q0Q8C7wHQJKAR9OHWWGySgpZV6/TJJcUtmJrbQH+jIqVZyC4B9hf0r7ABPA24O2NO0gaBp6LiN8A7wWWpcHBrHSyZrhsPsFNatcVFeoVINwWUG65BYKIeEHSucDtJN1Hr46IH0k6J91+JfC7wLWSNpE0Ip+VV3rMpqrbBuasksKlt/xoi4nwqtbo7LaAwePZR82mKKsvfFZJIUsVxip4XEB5+Z7FZjnKamAe7fJqd5Bm0MwaT+FxAYPJU0yY9UA3XVF3mr5d5gyagzAjarvxFG4LGEwOBGY56XYivOMOGCld43LWTYCyekm5LWAwuY3ArACDcJe1bts+BHzqrYdm3krSpYFi+Q5lZiXTq7usQW9KCt1c+WdNzDc5M2iv0mT94xKBWUlkzYjaTvMV+tZu5N7qhA+tq6ra9XrKOq5P+OXlm9ebDYCsqpisxuWsK/OsqqSsLpzdvv9oQ4nBV/2Dw1VDZgOg28blrCv2rKqkJd9b+ZIT+4aNmzLfZ1NEyyv/yZO+T/zV4UBgViLd3GUtq3E5S9a9FrL4yr8+HAjMBkAvxil0W5XkK//68MhiswGVNaL5klNe1XJq59OPmtNy/aV/+CqP+q05lwjMBlg3VUkLDxtlbJ+ZmVU9PvHXl3sNmZnVgCedMzOzTA4EZmY150BgZlZzDgRmZjXnQGBmVnMD12tI0lrgsW18+SzgyR4mZxA4z/XgPNfDVPK8T0SMtNowcIFgKiSNZ3WfqirnuR6c53rIK8+uGjIzqzkHAjOzmqtbILiq6AQUwHmuB+e5HnLJc63aCMzM7KXqViIwM7MmDgRmZjVXqUAg6WpJayQ90LT+fZIelvQjSR9rWH+RpEfSbQv6n+Kpa5VnSV+UdG/6+Lmkexu2VTXPh0r6bprncUlHNmyrap4PkfQdSfdL+qqklzVsG+g8S5oj6V8lPZj+bt+frp8p6Q5JP0n/zmh4TVXz/Jb0+YuSxppe05s8R0RlHsAfAIcDDzSsOw74f8CO6fM9078HAj8EdgT2BX4KTCs6D73Ic9P2TwAXVz3PwDeAk9Llk4G7apDne4DXpctnApdVJc/AXsDh6fJuwI/TfH0MuDBdfyFweQ3y/LvAfOAuYKxh/57luVIlgohYBjzdtPpPgcUR8Xy6z5p0/anAjRHxfEQ8CjwCHMmAycgzAJIE/DGwJF1V5TwHMHlFvDuwOl2ucp7nA8vS5TuAN6fLA5/niHg8In6QLv8SeBAYJcnbP6e7/TOwMF2ubJ4j4sGIeLjFS3qW50oFggyvBF4r6XuSviXpiHT9KLCyYb9V6boqeS3wRET8JH1e5TyfB1whaSXwceCidH2V8/wA8Ifp8luAOelypfIsaR5wGPA94OUR8TgkJ05gz3S3Kuc5S8/yXIdAsD0wA/h94HzgS+mVslrsW7W+tKezuTQA1c7znwIfiIg5wAeAz6Xrq5znM4E/k7ScpCrhN+n6yuRZ0q7ATcB5EfFsu11brHOeO1SHQLAKuDkS3wdeJJm4aRWbr6AA9mZzdcLAk7Q9cBrwxYbVVc7zu4Cb0+Uvs7mIXNk8R8RDEfGGiHg1ScD/abqpEnmWNJ3khPiFiJj83z4haa90+17AZFVvlfOcpWd5rkMgWAocDyDplcAOJLP33QK8TdKOkvYF9ge+X1Qic3AC8FBErGpYV+U8rwZely4fD0xWh1U2z5L2TP9uB/wlcGW6aeDznJbaPwc8GBGfbNh0C0nQJ/37fxvWVzXPWXqX56Jbynvc6r4EeBzYSBItzyI58V9PUp/6A+D4hv3/guQq6mHSHieD9miV53T9NcA5LfavZJ6BY4DlJL0ovge8ugZ5fj9Jz5IfA4tJZwqoQp7T/2cA9wH3po+TgT2Ab5IE+m8CM2uQ5zel//PngSeA23udZ08xYWZWc3WoGjIzszYcCMzMas6BwMys5hwIzMxqzoHAzKzmHAjMzGrOgcDMrOYcCMymSNIRku6TtJOkXdK543+v6HSZdcoDysx6QNJHgZ2AIWBVRCwqOElmHXMgMOsBSTuQ3Cjm18BrImJTwUky65irhsx6YyawK8l00DsVnBazrrhEYNYDkm4BbiS5ZeBeEXFuwUky69j2RSfAbNBJOgN4ISJukDQNuFvS8RFxZ9FpM+uESwRmZjXnNgIzs5pzIDAzqzkHAjOzmnMgMDOrOQcCM7OacyAwM6s5BwIzs5r7//Db9XiKmUSaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Model = loadobj('scales.nxs','scaling_model')\n",
    "Model.to_nexus().plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cfc56d4-bcdd-4af0-9f80-0eb18249e036",
   "metadata": {},
   "source": [
    "## 3) Merging scaled data\n",
    "\n",
    "Now that we have a scaling model, we can merge and then generate slices and volumes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "a12833c0-ea4a-4e0c-a491-59d886edbea4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading hkl_table from corrected.nxs\n",
      "  importing as HKLTable from mdx2.data\n",
      "Grouping redundant observations\n",
      "Reading scaling_model from scales.nxs\n",
      "  importing as ScalingModel from mdx2.scaling\n",
      "Calculating scales\n",
      "Merging\n",
      "Exporting <class 'mdx2.data.HKLTable'> to nexus object\n",
      "  writing hkl_table to merged.nxs\n",
      "done!\n"
     ]
    }
   ],
   "source": [
    "!mdx2.merge corrected.nxs --scale scales.nxs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b382863b-b5c4-467c-80a9-187b358631ad",
   "metadata": {},
   "source": [
    "### 2D slice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "5726928b-42e1-40e2-810d-03135557e5e2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading hkl_table from merged.nxs\n",
      "  importing as HKLTable from mdx2.data\n",
      "Reading crystal from geometry.nxs\n",
      "  importing as Crystal from mdx2.geometry\n",
      "Reading symmetry from geometry.nxs\n",
      "  importing as Symmetry from mdx2.geometry\n",
      "map region:\n",
      "  h from -50.0 to 50.0 (401 grid points)\n",
      "  k from -50.0 to 50.0 (401 grid points)\n",
      "  l from 0.0 to 0.0 (1 grid points)\n",
      "generating Miller index array\n",
      "mapping Miller indices to asymmetric unit\n",
      "looking up intensity in data table\n",
      "preparing output array\n",
      "Exporting <class 'mdx2.geometry.GridData'> to nexus object\n",
      "  writing intensity to slice.nxs\n",
      "done!\n"
     ]
    }
   ],
   "source": [
    "!mdx2.map geometry.nxs merged.nxs --limits -50 50 -50 50 0 0 --outfile slice.nxs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98c9432c-d40f-487b-b95c-a0c3dbdb194d",
   "metadata": {},
   "source": [
    "#### Visualization in *nexpy*\n",
    "\n",
    "Open `slice.nxs` in *nexpy* and adjust the sliders in the signal tab to make the background scattering variations clear.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>Questions:</b> Has scaling reduced the artifacts in the map? What artifacts remain?\n",
    "</div>\n",
    "\n",
    "The result should look something like this:\n",
    "\n",
    "![INSULIN](resources/insulin_slice.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0322454-2c9f-41f3-b073-59e73bf655fd",
   "metadata": {},
   "source": [
    "### 3D map\n",
    "\n",
    "We can also use *mdx2.map* to generate a 3D volume rather than a slice. This can be a bit resource intensive. So for the totorial we'll reconstruct one octant of reciprocal space (which is the same as the other seven by symmetry)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "90a97226-3589-4097-94c9-c32afe476ca4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading hkl_table from merged.nxs\n",
      "  importing as HKLTable from mdx2.data\n",
      "Reading crystal from geometry.nxs\n",
      "  importing as Crystal from mdx2.geometry\n",
      "Reading symmetry from geometry.nxs\n",
      "  importing as Symmetry from mdx2.geometry\n",
      "map region:\n",
      "  h from 0.0 to 50.0 (201 grid points)\n",
      "  k from 0.0 to 50.0 (201 grid points)\n",
      "  l from 0.0 to 50.0 (201 grid points)\n",
      "generating Miller index array\n",
      "mapping Miller indices to asymmetric unit\n",
      "looking up intensity in data table\n",
      "preparing output array\n",
      "Exporting <class 'mdx2.geometry.GridData'> to nexus object\n",
      "  writing intensity to map.nxs\n",
      "done!\n"
     ]
    }
   ],
   "source": [
    "# map of positive octant (h>=0, k>=0, l>=0)\n",
    "!mdx2.map geometry.nxs merged.nxs --limits -0 50 0 50 0 50 --outfile map.nxs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f8b235b-a42b-473e-a4bb-d498a1f3132f",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Visualization in *nexpy*\n",
    "\n",
    "Open `map.nxs` in *nexpy*. Scroll through z-slices. It helps to unselect \"autoscale\" in the z tab so that intensity scale is fixed, and then set adjust it in the \"signal\" tab.\n",
    "\n",
    "Spend some time looking through the map.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>Questions:</b> What sorts of diffuse scattering features can you find? How might you test whether a feature is real or an artifact of processing? How might you further process the data?\n",
    "</div>\n",
    "\n",
    "## Next steps\n",
    "\n",
    "In the workflow we developed with *mdx-lib*, the next steps are absolute intensity scaling and 3D-$\\Delta$PDF estimation. They have not been implemented yet in *mdx2*, so this is as far as we'll go in this tutorial. I encourage you to read more about them in our 2020 Nature Communications paper, and to try *mdx-lib* if you are interested.\n",
    "\n",
    "The final tutorial contains examples of how *mdx2* can be used in a python environment for crystallographic computing. It is intended for someone familiar with python. It complements the first 5 tutorials, but is optional.\n",
    "\n",
    "Alternatively, you might consider going back through the tutorials to explore different parameter choices. For example, what happens to the scaling model if you set \"--smoothness 0\"?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05a3e684-9f29-45b6-a99b-ae9402fc6fec",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}