{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`reciprocalspaceship` provides methods for reading and writing MTZ files, and can be easily used to join reflection data by Miller indices. We will demonstrate these uses by loading diffraction data of tetragonal hen egg-white lysozyme (HEWL)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9.9\n"
     ]
    }
   ],
   "source": [
    "import reciprocalspaceship as rs\n",
    "print(rs.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This diffraction data was collected at the Sector 24-ID-C beamline at [NE-CAT](https://lilith.nec.aps.anl.gov/) at APS. Diffraction images were collected at ambient room temperature (295K), and low energy (6550 eV) in order to collect native sulfur anomalous diffraction for experimental phasing. The diffraction images were processed in [DIALS](https://dials.github.io/) for indexing, geometry refinement, and spot integration, and scaling and merging was done in [AIMLESS](http://www.ccp4.ac.uk/html/aimless.html). This data reduction yielded an MTZ file that is included in the `data/` subdirectory. Here, we will load the MTZ file and inspect its contents."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Loading reflection data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reflection tables can be loaded using the top-level function, `rs.read_mtz()`. This returns a `DataSet` object, that is analogous to a `pandas.DataFrame`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'DataSet'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "refltable = rs.read_mtz(\"data/HEWL_SSAD_24IDC.mtz\")\n",
    "type(refltable).__name__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This reflection table was produced directly from `AIMLESS`, and contains several different data columns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>FreeR_flag</th>\n",
       "      <th>IMEAN</th>\n",
       "      <th>SIGIMEAN</th>\n",
       "      <th>I(+)</th>\n",
       "      <th>SIGI(+)</th>\n",
       "      <th>I(-)</th>\n",
       "      <th>SIGI(-)</th>\n",
       "      <th>N(+)</th>\n",
       "      <th>N(-)</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>H</th>\n",
       "      <th>K</th>\n",
       "      <th>L</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">0</th>\n",
       "      <th rowspan=\"4\" valign=\"top\">0</th>\n",
       "      <th>4</th>\n",
       "      <td>14</td>\n",
       "      <td>661.29987</td>\n",
       "      <td>21.953098</td>\n",
       "      <td>661.29987</td>\n",
       "      <td>21.953098</td>\n",
       "      <td>661.29987</td>\n",
       "      <td>21.953098</td>\n",
       "      <td>16</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>4</td>\n",
       "      <td>3229.649</td>\n",
       "      <td>105.980934</td>\n",
       "      <td>3229.649</td>\n",
       "      <td>105.980934</td>\n",
       "      <td>3229.649</td>\n",
       "      <td>105.980934</td>\n",
       "      <td>16</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>6</td>\n",
       "      <td>1361.8672</td>\n",
       "      <td>43.06085</td>\n",
       "      <td>1361.8672</td>\n",
       "      <td>43.06085</td>\n",
       "      <td>1361.8672</td>\n",
       "      <td>43.06085</td>\n",
       "      <td>16</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>19</td>\n",
       "      <td>4124.393</td>\n",
       "      <td>196.89108</td>\n",
       "      <td>4124.393</td>\n",
       "      <td>196.89108</td>\n",
       "      <td>4124.393</td>\n",
       "      <td>196.89108</td>\n",
       "      <td>8</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <td>16</td>\n",
       "      <td>559.33685</td>\n",
       "      <td>8.6263</td>\n",
       "      <td>559.33685</td>\n",
       "      <td>8.6263</td>\n",
       "      <td>559.33685</td>\n",
       "      <td>8.6263</td>\n",
       "      <td>64</td>\n",
       "      <td>64</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        FreeR_flag     IMEAN   SIGIMEAN      I(+)    SIGI(+)      I(-)  \\\n",
       "H K L                                                                    \n",
       "0 0 4           14 661.29987  21.953098 661.29987  21.953098 661.29987   \n",
       "    8            4  3229.649 105.980934  3229.649 105.980934  3229.649   \n",
       "    12           6 1361.8672   43.06085 1361.8672   43.06085 1361.8672   \n",
       "    16          19  4124.393  196.89108  4124.393  196.89108  4124.393   \n",
       "1 0 1           16 559.33685     8.6263 559.33685     8.6263 559.33685   \n",
       "\n",
       "          SIGI(-)  N(+)  N(-)  \n",
       "H K L                          \n",
       "0 0 4   21.953098    16    16  \n",
       "    8  105.980934    16    16  \n",
       "    12   43.06085    16    16  \n",
       "    16  196.89108     8     8  \n",
       "1 0 1      8.6263    64    64  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "refltable.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of reflections: 12542\n"
     ]
    }
   ],
   "source": [
    "print(f\"Number of reflections: {len(refltable)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Internally, each of these data columns is stored using a custom `dtype` that was added to the conventional `pandas` and `numpy` datatypes. This enables `DataSet` reflection tables to be written back to MTZ files. There is a `dtype` for each of the possible datatypes listed in the [MTZ file specification](http://www.ccp4.ac.uk/html/f2mtz.html#CTYPOUT). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FreeR_flag              MTZInt\n",
       "IMEAN                Intensity\n",
       "SIGIMEAN                Stddev\n",
       "I(+)          FriedelIntensity\n",
       "SIGI(+)         StddevFriedelI\n",
       "I(-)          FriedelIntensity\n",
       "SIGI(-)         StddevFriedelI\n",
       "N(+)                    MTZInt\n",
       "N(-)                    MTZInt\n",
       "dtype: object"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "refltable.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additional crystallographic metadata is read from the MTZ file and can be stored as attributes of the `DataSet`. These include the crystallographic spacegroup and unit cell parameters, which are stored as `gemmi.SpaceGroup` and `gemmi.UnitCell` objects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<gemmi.SpaceGroup(\"P 43 21 2\")>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "refltable.spacegroup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<gemmi.UnitCell(79.3439, 79.3439, 37.8099, 90, 90, 90)>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "refltable.cell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Plotting reflection data\n",
    "\n",
    "For illustrative purposes, let's plot the $I(+)$ data against the $I(-)$ data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAFzCAYAAADPISX/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABFkUlEQVR4nO3deXCc52Hn+e/TaHTjRqMBEGgSAAmQ4AEetkiK1jmW7ciWLMfyrhOPxjWxE3ujqsTZzVR21keyk8xOjSuZrUlqxslYZTnjWI7tcWRlFGs1lhxJGUVR6SRpkSIJEhdxN87G9QJgv2jg2T/Q3W6CuEihATTw+1Sh0P30+3a/ry32D89trLWIiIgsx7PRFyAiIpufwkJERFaksBARkRUpLEREZEUKCxERWZHCQkREVuTd6AtIl7KyMrtnz56NvgwRkYxy5syZIWtt+cLyLRsWe/bs4fTp0xt9GSIiGcUY07FYuZqhRERkRQoLERFZkcJCRERWpLAQEZEVKSxERGRFCgsREVmRwkJERFaksBARkRUpLEREZEUKCxERWZHCQkRki3BdF8dxcF13zd97y64NJSKynbiuS0dHB9ZajDHs3r0bn8+3Zu+vmoWIyBbgui7WWgoKCrDWrnntQmEhIrIF+Hw+jDE4joMxZk1rFaBmKBGRjOW6Lq7r4vP58Pl87N69+7rna0lhISKSgZbqo1jrkEhIazOUMSZgjHnKGHPZGNNojLnTGBM0xrxgjGmO/y5JOf5rxpgWY8wVY8zHUspPGGPejb/2DWOMSed1i4hsdunuo1go3X0W/xl43lp7EHgf0Ah8FXjJWlsPvBR/jjGmAXgEOAw8AHzTGJMVf5/HgEeB+vjPA2m+bhGRTS3dfRQLpS0sjDFFwD8D/iuAtda11o4CDwNPxA97AvhU/PHDwI+stVFr7VWgBThljAkBRdba1621FvheyjkiIttSoo9i586daz5MdjHprFnUAYPAXxljfm6M+UtjTD5QYa0NA8R/74gfvwvoSjm/O162K/54YfkNjDGPGmNOG2NODw4Oru3diIhsMj6fj4KCgrQHBaQ3LLzAceAxa+1twCTxJqclLNYPYZcpv7HQ2settSettSfLy8tv9npFRGQJ6QyLbqDbWvtm/PlTzIdHf7xpifjvgZTjq1POrwJ64+VVi5SLiMg6SVtYWGv7gC5jzIF40UeAS8AzwOfjZZ8HfhJ//AzwiDHGb4ypZb4j+614U9WEMeaO+Cioz6WcIyIi6yDd8yz+d+AHxhgf0Ab8BvMB9aQx5otAJ/CrANbai8aYJ5kPlBjwJWvtbPx9fgv4LpALPBf/ERGRFAsn6a0lMz/AaOs5efKkPX369EZfhojIulirhQSNMWestScXlmttKBGRLUALCYqIbAHp3GsCtJCgiEjGS/deE0DaFxJUzUJEJM1WaiJaq1pHOifpqWYhIpJmyzURrUetYy0oLERE0my5JqLUWkeidqGwEBHZppbqR1jv1WNvlcJCRGQDpbtjeq2og1tEZAOlc9b1WlLNQkRknSwMhkzp3AaFhYjIunAch5aWFrxeL36/P9n0lAmd26BmKBGRtHNdl5aWFvr7+xkdHSUajSaDIRM6t0E1CxGRtHNdF6/XS35+PpOTk8mJc5nSuQ0KCxGRtPP5fPj9fgKBAAUFBezbty8ZDJs9JBIUFiIiaZZJNYilKCxERNbYYsNhMzUkEhQWIiJrKJOGw94MjYYSEVlD6d6EaKMoLERE1lAmDYe9GWqGEhFZQ1uhM3sxCgsR2dbSsTbTVgqJBIWFiGxbW7UzOh3UZyEi29ZW7YxOB4WFiGxbW7UzOh3UDCUi29atdEZnyv4Ta01hISLb2s186S/s4wiFQjf9HplKYSEiskqpfRyRSISWlhby8vK2Ree4+ixERFYptY8jFovh9Xq3Tee4ahYisq3dTB9Eah8HQDgc3jad4woLEdm2VjvPYmGgJI4JhUI4jpPczGgrU1iIyLa1mj2wlwoU13UJh8NYa5mYmFCfhYjIVrWaeRZLTdzbbhP6VLMQkW1rNfMslgqU7TahT2EhItvaSh3bSwXKVl1ddikKCxGRFSxX69jqIZGgPgsR2fZc1012cMviVLMQkW1jsTkVWqZ8dRQWIrItLDcEdqXhs6JmKBHZJpYa6rrdRjXdqrSGhTGm3RjzrjHmHWPM6XhZ0BjzgjGmOf67JOX4rxljWowxV4wxH0spPxF/nxZjzDeMMSad1y0iW89yQ2B3797Nzp071QS1jPWoWXzIWvt+a+3J+POvAi9Za+uBl+LPMcY0AI8Ah4EHgG8aY7Li5zwGPArUx38eWIfrFpEtZLlQ8Pl822LJjvdiI5qhHgaeiD9+AvhUSvmPrLVRa+1VoAU4ZYwJAUXW2tettRb4Xso5IiKrplC4dekOCwv8vTHmjDHm0XhZhbU2DBD/vSNevgvoSjm3O162K/54YfkNjDGPGmNOG2NODw4OruFtiIhsb+keDXW3tbbXGLMDeMEYc3mZYxfrh7DLlN9YaO3jwOMAJ0+eXPQYEcl8N7Os+HbdBnWtpTUsrLW98d8DxpingVNAvzEmZK0Nx5uYBuKHdwPVKadXAb3x8qpFykVkG7qZeRGaQ7F20tYMZYzJN8YUJh4DHwUuAM8An48f9nngJ/HHzwCPGGP8xpha5juy34o3VU0YY+6Ij4L6XMo5IrLN3Mxqr9ttZdh0SmfNogJ4Oj7K1Qv80Fr7vDHmbeBJY8wXgU7gVwGstReNMU8Cl4AY8CVr7Wz8vX4L+C6QCzwX/xGRbehm5kVoDsXaMfMDjLaekydP2tOnT2/0ZYhIGqjPIn2MMWdSpjokabkPEck4N/PFr5BYG1ruQ0QErTy7EtUsRGTb06iplalmISLrbi3/il+L99KoqZWpZiEi62ot/4pfq/fSqKmVKSxEZF2t5f4Ra/Ve220/7VuhsBCRdbXYX/E3O7w1cTywZjUChcTyFBYikjaLhcDCv+KBm2pK6uzs5Ny5cwSDQXbs2EEoFEq+r77s00dhISJpsVx/QuoXu+M4q25K6uvr4zvf+Q6jo6Pk5ubyiU98glAoREFBwbrd13al0VAikharHWF0M53LAwMDGGOorq5menqaSCSi2sQ6Uc1CRNJitSFwM53LO3bswOfz4TgOgUCA973vfQqLdaKwEJG0uJkQWOr1SCTCyMgIJSUlBINBKisr+Y3f+A0GBgbYsWMHlZWV6bwFSaGwEJG0eS+dzpFIhJ/85CfMzs6SlZXFww8/nAwMhcT6U5+FiGxK/f39TE5OUllZyezsLCMjIxt9SduawkJENpTrukQiESKRSLITPLGEx8TEBOfOncNaS0lJyQZf6famZigR2TCu69Lc3ExPTw8Au3btor6+Htd1KS4u5pd/+Zfp6upi9+7dGh67wVSzEJEN47ouk5OTWGvxer1MTk4SiUSA+ZnZAF6vl7m5OTo6OrTA3wZSzUJENtTQ0BDd3d1Eo1FCoRB5eXlMTEwQCoVwHAeAYDD4nteRkvdGNQsR2VChUIiKigpycnKYmZlJTuKD+ZDw+/1aDXYTUM1CRDaMz+cjFosxNDSE3+9nZGSE9vZ2ysvLk8FQWloKQEFBgcJiAyksRGTD+Hw+qqur6enpIS8vj46OjuRrrusSDoeTa0upg3tjqRlKRNZd6u52wWCQ+vp68vPz2blzJ3v27CE7OxvHcYhGoxhjiEaj6tzeYKpZiMi6Wmw12vr6ehzHIRwO47pusn8iHA4zNzeHx+Nhz549G33p25rCQkTSYrG9LBzHob+/H8dxqKioSNYuCgoKCAaDFBQUJM9xXZdQKITf7ycajW7w3YjCQkTW3GK1h0gkwj/+4z+SlZXFyMgIDQ0NBAKB6zqtF64l5ff7sdbi9/vVub3BFBYisuYW7o0diUR4+eWXaWlpoaioCGstU1NTBAKBJd9D+2JvLurgFpE1t3AvC9d18fv9BINBBgcHmZubY/fu3WRnZy/bce3z+TRkdpNQzUJE1tzCWkGiXyIUClFUVERdXR2AJtplEIWFiNyyxTqxE1LLfD4fd955J47jJGsKqecl3mexc2VzUFiIyC1ZrBN74Rd8apgUFBRcN7EucWzifaLRKOFwODkCarH3k42jsBCRW+K6LtFoNDm0NXWRv8Sku3A4THZ29pJhkjg2MeJpbm4uOQJKiwZuLgoLEbllCyfNpYaE67oMDg6yf//+62oYCyU6w6PRKB6PJxlACorNRWEhIjcl8cW/cNJcJBKhu7sbay2Tk5PU1tYyODjIyMjIsiOaUjvDE7O01Wex+SgsRGTVUvspZmZmiMViyd9XrlxhZGQEv99PYWEhjuOwa9cuQqHQisNfFQ6bn8JCRBa12Ail1H4K13WJxWJ4vV6mpqbweDyUlZUxNDRETU1NcivU1BBYbvSUbG4KCxG5wVIjlEpLS5P9FNFoNLnXhOu6yc7p8vJyjh49SjAYXPQ9lxs9JZuXwkJEbrDcCKVEcPT09PD6668zPT2N1+vlrrvuora2Nrkg4EKO4+A4DiUlJct2eMvmpOU+ROQGi41QSmxAlNjmtKOjg7m5OfLz89mxYwder3fJoEhsZDQ4OEhTUxMzMzMKigyjmoWI3GC5EUqlpaW0tLQwOzvL7Ows09PT5OTkkJ+fv2QAuK5LdnY2+/fvZ2RkhFAopLDIMGmvWRhjsowxPzfGPBt/HjTGvGCMaY7/Lkk59mvGmBZjzBVjzMdSyk8YY96Nv/YNY4xJ93WLbHeps65Tl+hobW2lq6uLkZER8vLyOH78OJ/4xCdoaGhYdnhsYkHBhTO5JTOsRzPU7wKNKc+/Crxkra0HXoo/xxjTADwCHAYeAL5pjMmKn/MY8ChQH/95YB2uW0RSuK5LJBJJfuGXlZVRUVHB7bffTk1NzYpDY3fv3s3OnTvVsZ2h0hoWxpgq4CHgL1OKHwaeiD9+AvhUSvmPrLVRa+1VoAU4ZYwJAUXW2tettRb4Xso5IrIOHMfh0qVLhMNhxsbGyMvLIxAIcPDgwRtGPS1Fy41ntnT3Wfwn4MtAYUpZhbU2DGCtDRtjdsTLdwFvpBzXHS+biT9eWC4it2ix+Q5LzYFwXZeWlhb6+/vJz88nFApRWlp6XfOUbH1pCwtjzCeAAWvtGWPMfas5ZZEyu0z5Yp/5KPPNVdTU1KzuQkW2mcXmOwBLzoFwXRev10t+fj6jo6N4PJ4lRz0t/BxNwNs60tkMdTfwSWNMO/Aj4MPGmO8D/fGmJeK/B+LHdwPVKedXAb3x8qpFym9grX3cWnvSWnuyvLx8Le9FZMtI3fI0MXdiYVliTkRiBncsFsPn8zE7O0tubm5yocDlPqOjo4Pe3l46OjqWPVYyQ9rCwlr7NWttlbV2D/Md1/9grf2XwDPA5+OHfR74SfzxM8Ajxhi/MaaW+Y7st+JNVhPGmDvio6A+l3KOiNykxMikSCTC1NTUdWWO4zAzM0M4HKa3t5fm5mY6Ojrwer3EYjFqamooKSlJhslSFgskyWwbMSnvT4D7jTHNwP3x51hrLwJPApeA54EvWWtn4+f8FvOd5C1AK/Dcel+0yFaRmCsxPj6OtZZwOAyQHK0UCoXIzs6moKAgWesIBoMUFRURi8VoampicHBw2drFwj241QyV+dZlUp619mXg5fjjYeAjSxz3deDri5SfBo6k7wpFtg/Xdens7GRiYoK5uTmys7OTw2ETr8/MzOA4TvJL3nEc/H4/ZWVleL3eFZfsWLgHt8Ii82kGt8g2k9phPTk5mRzR5DgOLS0teL1ejDGUlZVdFyCJL/yJiQlc1yUxNzYRKsvtwS2ZT2Ehso0kagPGGAKBAAUFBezbtw/XdTl//jz9/f3JCXepX/apX/qJGgPM75SnVWS3B4WFyDaROmQWSAaC67pcvnyZjo4O2tvbKSwsZHp6mvr6+kXfI1HLSO3EToycUlhsXQoLkW0gsVRHNBolGAzS399PW1sbubm5dHV1ce3aNdrb25mamiIQCFBZWbnoe6TOxQiFQurE3kYUFiJbmOu6OI6TbC4Kh8PMzMzQ1dWFx+PBGMPk5CQ+n4+ioqJkUCzW37CwJgGoE3sbUViIbAFLLd/R0dGB4zgMDg6yf/9+QqFQcmKd67r09PQQi8UoKSmhqKiIuro68vPz2bdv36Id1gtrEgqJ7UNhIZLhltquNLFfts/nY2ZmhpGREYwxdHR0MDQ0xNTUFOXl5Rw8eJD29nbKy8uTQbHYUh4aDru9KSxEMtxyHc2J/bITfQyu69Lf38+RI0fo6OigoqICay15eXns2rVrxZnWContS2EhkuGWmy2d2C97ZGSEkZERsrOzmZ2dZXBwEL/fT21tLUByYyJ1VMtSFBYiGc7n8xEKhXAcJ9l8lNikKBaLMTY2xvnz58nKysLr9bJ3714AcnNzuXLlSnJ5j8QkPIWFLEZhIZLhXNclHA4TjUbp6uoCYGRkhJmZGa5du8bU1BQDAwPU1taSlZWF67rJYBgZGcHv92OtvakmJi0/vv0oLEQyXKIje3R0lMHBQaLRKDt27GB6epqmpiZGRkZoa2ujq6uL48ePU1tbSywWw3EcPB4P0WgUv99/U0Gx1N4XsnUpLEQy0MK/7GOxGJOTkwQCAcbGxujv76enp4fe3l56enoIBALk5uZy5MgRiouLsdYSi8U4ceLETQ+B1czt7UlhIZJBUifZZWdnJ/+y37dvHwBer5eCggLC4TAVFRXJvSkKCwsJBAIEAoHk8uOJBQBX2vFuIS0/vj0pLEQyxGKT7BI1jIKCAhoaGnAch9bWVmZnZ4nFYlRUVFBbW0swGCQQCLBjxw4mJibe0xe95ltsTwoLkQyRaP4pKSm5rnkpdWVYn8+H1+vl2rVruK6L3+/noYceIhAIJJuOgsHge/6iV0hsPwoLkQyxsPknOzsb+EX/ReJxLBajsLCQvLy85Jd6Xl6eluiQ90RhIbKJLezIDoVC9Pf3EwqFqKioIBKJ0NLSgrWWtrY2gsEgHo8Hv99Pe3s7Ho8H13Wpq6sjGAwqJOSWKSxENqnFlgRPzKcYGhoiOzub6elpXNdlaGiIS5cu4fV6KS0tpbKykkAgQE1NDaOjowA3BIXmSsjN8Gz0BYjI4lKHqEajUfr7+6/bC7uoqAiArq4uWltbMcbg8Xjwer0UFRXh8/kYHR3F4/Ekt01Nba7q6Oigt7eXjo6OFdeEElHNQmSTSvRRRCIRwuEwxcXFXLx4kfLycvx+P6WlpQwMDJCbm4sxhqqqKqampggGg5SXl3Ps2LFkzWF4ePi6SXQ3O1dCtRBRWIhsUql9FInlOcrLy5MrxTqOQ39/P36/n4KCAk6dOpXsl0jMnUh8yS8WDKudK6EZ2wIKC5FNY+Ff76lrPg0NDQHg8XgYGBggOzubWCxGcXExXq+X/Px8KisrCQaDyfdKfMHPzMwA3LBp0WrnSmjGtoDCQmRTWOyv98SXdCIAysrKkp3cJSUljIyMADA3N3fD2k4Lv+DLyspuGDa72iYlzdgWUFiIbAorNRVNT0/T399Pfn5+stZhjKGuri7ZAZ5q4Rf8e1l6XDO2BRQWIpvCwi93mA+QUChEJBLhpZdeIhKJ4PF4+OAHP0hdXR1FRUXJjuuFNYu1/oJXSIjCQmSDJfoqQqFQsiwcDmOtZWpqivHxcSYnJ6msrCQSiTA6OkpfXx9er5fp6Wmqq6sXnXCnL3hZSwoLkXW0WCf2wol3juMkm5Vee+01cnJy6O7uZmxsjLy8PPLz87HWMjg4mJxwl+jXEEkXhYXIOlmuE7ugoID+/n7Onz9Pbm4ura2thMNhhoaGOHLkCB/4wAfYtWsXe/bsIRgM0tjYSFtbG7m5uQwMDOA4znUjodS/IGtNYSGyTpbrxI5EInR2dpKVlYUxhsuXLxONRpmamqK/v5/q6mqOHz+enD9RXV1Nb28vhYWFxGKx5D4X8IsmLM2JkLWksBBZJ4sNQU10REciEaampmhra+Ptt99mZGSEnTt3kpuby86dO7n33nuv26QoGAyyb9++5Kio4eFhxsfHmZqawuv1EgwGNSdC1pTCQmSdpI5QApI1gcSXeVtbGy+//DIDAwPEYjF27NhBYWHhorUDn89HfX19sslpaGiIgoKC5BLlmhMha01hIZJmqftNJJ63trYyMDDAzMwMs7OzVFZW4jgOpaWlFBYWcvXqVWZmZrh27RqdnZ1EIhHuvPPO62oXqZ3kw8PDyUUG9+zZc93rImtBYSGSRolO7Wg0mlwM8OrVq8zNzTE9Pc3o6CiTk5NMTExgrWVsbIxYLEZZWRm33347AwMDBAKBZJ/EYvtla9KcrAeFhUgaJTq1/X4/0WiUzs5Ouru7mZubY2BggP7+fgYGBggGgxw4cID777+f/Pz85FpQk5OTNDc3U1VVtWhQJCgkJN1WHRbGmB3A3cBOYBq4AJy21s6l6dpEMtLCbU5nZmaw1nLt2jV6enro6elJlg8MDBAOh5meniYQCODxeCgtLSUvLw+v18uRI0cYGxtj7969y4aFSLqtGBbGmA8BXwWCwM+BASAH+BSw1xjzFPCn1trxNF6nSEZIbXbq7OxMLuBXU1NDUVERL730EjDfud3Y2MjY2BhTU1MYYxgfH6enp4eCggIKCgrIy8sjOzubHTt2aNKdbLjV1Cw+DvymtbZz4QvGGC/wCeB+4G/X+NpEMk6i2ckYQzgcxuud/ydWWVlJQUEBlZWVNDU1MT09zfT0NF6vl8LCQvLy8qipqWF2dpaZmRlGRkaorq5OLgC4cEVZ9U/IelsxLKy1/9cyr8WAv1vLCxLJZIm5FImO6+npaay1dHV1UVxcTF5eHoWFhczNzZGVlUUsFsPv93P48GFqa2ux1uL1eonFYslNjBKd24lg0EZEshFuqYPbGPOstfYTa30xIpkusbvd4OAgXq+X7u5uysrKkiOdzpw5g8fjwRjD+973PjweD2VlZTz00EPk5uYSi8XIzc29LihSw6G0tFQbEcmGuNXRULtWOsAYkwO8Avjjn/OUtfaPjDFB4G+APUA78Blr7Uj8nK8BXwRmgf/DWvuzePkJ4LtALvBT4HettfYWr10kbRK1AGstO3fu5Nq1awCcPXuWoaEhGhsbkx3ZHo+HQ4cOEQgEKCoqoqCgILnybKKJKfFeiXAAtBGRbIhbDYufr+KYKPBha61jjMkGXjXGPAf8r8BL1to/McZ8lfnO868YYxqAR4DDzI+4etEYs99aOws8BjwKvMF8WDwAPHeL1y6SFolagOM4nD17lp6eHgYHB4lGo5SVleHxeHBdl6tXr5KVlUV9fT0HDx7k9ttvv6FvInWtp4WbGCVqHOqzkPV0S2Fhrf3CKo6xgBN/mh3/scDDwH3x8ieAl4GvxMt/ZK2NAleNMS3AKWNMO1BkrX0dwBjzPeZHYiksZMOlDpNNLC3u8/mScyuGhobIyclhYmIiObopLy8Pv99PSUkJwWAwGQCp77lw2fLUGeAKCdkIqxk6+/8BjwPPW2tnFrxWB/w60G6t/c4i52YBZ4B9wH+x1r5pjKmw1oYBrLXh+PwNmG/aeiPl9O542Uz88cLyxa71UeZrINTU1Kx0ayI3Zam9KBKzs8vKygiHwxhjiEajyT6Fubk5cnJyCIVCHDt2DMdxaGpqIhaL0dfXx6FDh274nIX9Eokd8YaHh9WpLRtiNTWL3wR+D/hPxpgIMMj8PIs9QCvwF9banyx2YrwJ6f3GmADwtDHmyDKfYxZ7i2XKF/u8x5kPNk6ePKk+DVkzi+1F4ThOcpRSNBolFoslh8Qmlg2vrKxkbGyMoqIiampqOHDgAI7jUF5eTmFhIVlZWXR2dl63R/ZiW6yqU1s22mqGzvYBXwa+bIzZA4SYn8HdZK2dWs2HWGtHjTEvM9/X0G+MCcVrFSHmJ/nBfI2hOuW0KqA3Xl61SLnIuln4135fXx9NTU0MDQ3h9XoZGRnh2rVrtLe3k5+fz+zsLHl5eVRUVDA2NsZ9992H1+uloqIiORs7EomQn5+P1+u9LgAWrvUEJBcKVKe2bJSb6rOw1rYzP4JpRcaYcmAmHhS5wC8B/wF4Bvg88Cfx34layTPAD40xf8Z8B3c98Ja1dtYYM2GMuQN4E/gc8Oc3c90i71XqX/tTU1OcPXuW1tZWPB4PJSUl7N27F5/Px7PPPsvY2Biu63L48GEqKiooKirCWovH40n2Ufh8PlpaWrDWEovFFv281FDQQoGy0VbTZzHB4s0+hvl+7KIlTg0BT8T7LTzAk9baZ40xrwNPGmO+CHQCv8r8G100xjwJXAJiwJfizVgAv8Uvhs4+hzq3ZZ2l/rXf19fH9PQ0Ho+H4eFhvF4vAwMD9Pb2Mj4+TmFhIaOjo/h8PoqKivjYxz6WnDeRGA7r8/nYt28fLS0teL1ewuHwsn0RCgnZaKtphiq8lTe21p4HblukfBj4yBLnfB34+iLlp4Hl+jtE0mZhx/bIyAjd3d10d3cTi8U4cuRIMiyi0SgTExNkZ2eza9cuamtrk7WJxSbY5eXlqS9CMoKWKBdZhuu6NDc3J7/IQ6EQXq+X48ePk5eXx8jICK+++ipjY2NYa5mbm8MYg9frZWpqiqGhoeSciUS/RyJwCgsLNcFOMobCQmQZjuPQ09NDbm4u4+Pj+P1+jDHMzMzg9XqZnJzk0qVL+P1+YrEYgUCAo0ePkp+fz65duyguLqazs5O8vDxmZmaYmZmhvb0d+EXTVuKxwkI2M4WFyCrMzMzQ399PIBDA5/OxY8cOXn/9dRobGxkdHaW4uBiAoqIiPB4P1tpkc5PX6002NRUWFmKtpaSkJDnRTvtUSCZQWIgsw+fzUVJSwuTkJMFgENd16ezs5Ny5c/T19TE7O0tZWRmxWIxgMMjx48cpLi4mGAxSUlJCdnb2dU1NwWCQaDSK67pqepKMorAQWSB1CY/EnhTGGEpKSohEIrzzzju8++67DA8PMzs7SzAY5ODBg+zbt4+ysjJmZ2eTNYuenh6OHz9OMBhMNjVpGKxkIoWFSIrUEUtTU1NYa5mcnGRiYoKsrCwGBwdpbW3FcRxmZ2epr69n//79NDQ04Pf7ASgpKWFkZCQ5xLa/v5/KysrrJt0pJCTTKCxEUriuSzQaxe/3Y61lenqa0dFRBgcHGRoaoqmpiXA4TFZWFllZWRQVFeH1evF4POzYsYO6ujoKCgpobGykp6eHUChEbm6uhsVKxlNYiCwQDoeZm5tjdnaW2tpa+vr6aGtrY2hoiMnJSUpKSvD5fAQCAe68805CoRB5eXnk5uYmV5A9dOgQ2dnZeL1e/H6/gkIynsJCZIHi4mLGx8cJh8NEIhGuXLnC6dOnmZubY2BggFAoRGVlJadOnWL//v04jsP09DTBYBAgOUO7oaFBfROyZSgsRFI4jsObb75JV1cX4XCYQCBALBZjZmaGoqIiQqEQR48e5fDhwxw6dIhDhw4lO8N9Ph/hcPi6lWk1LFa2CoWFSJzrurS1tTExMUFHR0dy+Y7EPhXT09OUlZVx4MAB9u3bx6FDh64Lg4VboKqfQrYShYVsKwvXeVr4WmJRwMuXLzMzM4PH42F0dJSSkhIqKyuprq5ObmK0sNawcB8KBYVsJQoL2TYW28Ao9Qvd5/MxPT3N5OQk2dnZzMzMYIxhbm6OkpISysvLqaio4PDhw4s2L2kOhWxlCgvZUlaqOSSaiSKRCJFIJDlZDqCvr4+LFy8mm5v6+/uprq4mNzeXvLw88vPzufPOO6msrFzy8xUSslUpLGTLWE3NwRhDJBIhHA4D8zvQhUIhIpEI3/72t+nt7aW3t5eioiJ27drFzp072bt3L2VlZRQWFlJbW6swkG1JYSFbxsKaQ19fX3LeQ+pSG319fUxNTRGLxbhw4QItLS1cvHiR1tZWCgoKMMZQXFyc3ArV5/PR39+P67pcuHAhuT+FyHaisJAtI7Xm0NnZmVysb/fu3TQ0NAAQiURoamqivb2dS5cuMTY2hjGGyclJRkdHmZ2dpaKigoaGBnJycsjLy8Pn8zE0NERpaSlzc3M4jqOwkG1HYSFbRmJzov7+fvLy8hgdHcVaS3t7O5WVlfT19XH69GnOnj1LTk4OV65cITs7G9d1KSws5NixY2RnZ1NbW0t5eTkzMzMEAgG8Xi95eXnJJisFhWxHCgvZMlzXJRwOE41G6evrY2hoiEAgAJAcDtvU1MSVK1fweDx0dXVRXl5OdnY21dXVNDQ0sGfPHnJzc5mbmwPA4/Ek38d1XU6dOqWwkG1JYSFbRmqfRU5OznXLb3R3d9Pc3Mzly5dpbGzEcRy8Xi+FhYXs2LGDAwcO0NDQwMGDBwmHw/T09ADzK8jW1NRQUFBANBpVUMi2pbCQjJe6/4QxhpGREfLy8ti3bx/vvPMOAwMDRCIRmpubaWtrY3JykqysLPLz85mbm6O0tJQDBw5QVFSEz+ejvr6eUCgEXL+EhxYElO1MYSEZbeFw2VAolJxn4bpucjvUgYEBzpw5Q3t7O9PT02RlZZGbm0swGOTee+/F6/USi8WS75kYQQVoop0ICgvJcKlNT47jACSHtnZ2dtLZ2UlzczNXr16lra2Na9euAZCVlcXJkyf55Cc/SXV1Nfn5+dTU1NywEGAiIBQSst15NvoCRN4Ln8/HzMwM3d3djI2NAfMB4jgOzc3NjI+PMzc3l5wnkeD3+7nrrrvYuXMneXl55OXlASSDx1p73fEi251qFpKRUvsppqamaGtrIycnh+zsbLKzsxkdHeXMmTP09vZy5coVsrKyKCgoYGpqCr/fz6FDhwiFQoRCISoqKpK1Ei0EKLI4hYVknNR+irGxMVpbWxkfH2diYoL8/HyKi4sZGhri4sWLnDlzhvHxcbKzs9m1axfXrl1j//793H777ezevfu6cEjM9lb/hMiNFBaScRL7ZMdiMTo7O7l69SpjY2NEo9HkjOvGxkbeffddHMchOzubiooKKisrOXDgALOzsxQWFhIIBNi3bx9w/QKACgmRGyksJCO1trZy/vx5JiYmmJ2dpaamhry8PAoKCnj33Xf5+7//e4aHh5mZmWF2djbZ6T05OUllZSU7d+6ksrJy1fMmllvNVmQ7UFhIRsrLy8Pj8VBRUcHQ0FDyC/x//I//QXNzM729vcB8R7YxhoaGBqqqqpK1ipycnFV/6a+0mq3IdqDRUJJxEgsGjo2NMTg4SE5ODiUlJXR2dvL222/T29vLxMQExhhKSko4dOgQ1dXVjI+PMzU1xcjISPJ9ViN1eK5GScl2pZqFZBTXdYlEIgwPDzM7O0tbWxsVFRWMjIzQ1NREJBIhFotRXFxMVVUVt912G7FYjNLS0uTop8LCQsrLy1f9mdouVURhIRnEcRxaWloYHR2lu7ubsrIyRkZGmJiY4NKlSwwPD2OtJScnh7179/LQQw9RWFiY3OjI4/EwMTFBYWEh+fn5q/7S13apIgoL2aRSO5RhfsvTd955h8nJSXw+H+3t7YyMjHDhwgVmZ2dxXRePx0Nubi4HDx7kM5/5DLfffjt9fX2UlJQwOjpKIBAgPz+furo6Kisrb+pLXyEh253CQjaN1Il2HR0dyc2LZmZmeOONN2hpacHr9SYX9JuZmSESiTA7O8vs7Czl5eUEAgHuvfdejh49SkVFBdeuXWN0dJRoNApAeXn5TQeFiCgsZJNINDF5vV6mp6cZGBigqKiISCQCzC/DkZuby+DgIH19fclVZGdnZ5Nbn9bU1HDo0CE+/OEPs3v3bmB+EcDUxQVTFwgUkdVTWMiGc12XlpYW+vv7yc/PT9YmZmZmmJ6eZnR0lI6ODrq6uujr66Ovrw/HcRgdHU3uMxEKhbjzzjs5cOAAFRUViy4IKCK3TkNnZcO5rovX6yU/Pz/ZJ7Fr1y7Gx8eJxWI0NTUxPj7O8PAww8PDTE5OEo1GMcbg8XgIBAI8/PDDHD58mD179uDz+TTUVWSNqWYhG87n8+H3+wkEAvj9fqqrqwHwer309fVx5coVIpEILS0tyS9/r9dLSUkJpaWl1NTUUFdXR1VVFYcOHcLn8zE8PKyhriJrSGEhGy4xNDUSidDV1UU4HCYWi5GdnZ2sSXR0dBCJRPB4PGRlZZGdnU1ZWRkHDhzg3nvvpba2lr179yaX79BQV5G1lbZmKGNMtTHmfxpjGo0xF40xvxsvDxpjXjDGNMd/l6Sc8zVjTIsx5oox5mMp5SeMMe/GX/uGMcak67plY7iuS1tbGy0tLbzxxhtcvXqVoaEhGhsb6evrIxwOMzc3x+zsLDk5OdTW1nLkyBF27dpFKBSivLz8unWe1JktsrbSWbOIAf+ntfasMaYQOGOMeQH4deAla+2fGGO+CnwV+IoxpgF4BDgM7AReNMbst9bOAo8BjwJvAD8FHgCeS+O1yzpKdHC3t7fz6quv4rpucp+J8+fPMzU1lTzWWktWVhZFRUUcOnSI+vp6Tp06RXV1tYJBJI3SFhbW2jAQjj+eMMY0AruAh4H74oc9AbwMfCVe/iNrbRS4aoxpAU4ZY9qBImvt6wDGmO8Bn0JhsWUk5le0tbXR29uLx+Ph6tWrGGO4cOECk5OTwPwM7KKiIm6//XaOHDnC8ePH2bdvn4JCZB2sS5+FMWYPcBvwJlARDxKstWFjzI74YbuYrzkkdMfLZuKPF5Yv9jmPMl8DoaamZg3vQNZK6sQ7+MVifk1NTVy4cIH+/n5GR0cJh8NkZ2czOTmZ3PI0sS9FQUEBwWCQQCCgYbEi6yTtYWGMKQD+FvhX1trxZbobFnvBLlN+Y6G1jwOPA5w8eXLRY2TjJJb6jkajhMNhQqEQxhhycnIYHBxkfHycaDRKd3d3cmjs7OwseXl51NTUUFJSQkVFBcXFxZSUlFBUVLTRtySybaQ1LIwx2cwHxQ+stf89XtxvjAnFaxUhYCBe3g1Up5xeBfTGy6sWKZcMklgtNhECjuMwNTWVrEFcvnyZoaEh+vr6mJqaYm5uDiC5Her9999PUVERXq+Xubk5ysvLNSxWZB2lLSziI5b+K9Borf2zlJeeAT4P/En8909Syn9ojPkz5ju464G3rLWzxpgJY8wdzDdjfQ7483Rdt6y9RI0isYy41+slEolw4cIFAEpKSmhubqavry+51wTMN1GVlJRQX1/Pnj17OHDgADC/NEhFRQX79u1TWIisk3TWLO4Gfg141xjzTrzs95kPiSeNMV8EOoFfBbDWXjTGPAlcYn4k1ZfiI6EAfgv4LpDLfMe2OrcziOu6OI5Da2srQ0NDTE9PU1JSQkdHBzMzM4yMjHDlyhWi0SjXrl0jJycHgFAoRDAYZO/evWRlZZGfn09RURGxWIx9+/atektUEXnv0jka6lUW728A+MgS53wd+Poi5aeBI2t3dbLehoeHiUQiZGVl0dTUxMzMDP39/QQCAXp6epicnEwu31FYWEhRUREf+tCH2L17N0eOHCE7O5u6urrk3AnVKETWl2ZwS1q5rks4HCY3N5fOzk6amppobW1ldHQUx3GYnZ0lFosxMzMDzA+PLSgoYP/+/fzSL/0SdXV1yVVlg8GgQkJkgygsJG1SO7UDgQCFhYXMzs6SlZWF4zi4rsu1a9eSHdWJPbMPHjzIkSNHsNZSXV2t2oTIJqCwkDWTOociUaOw1tLZ2UkkEuHs2bPJFWSvXbtGLBbDGIO1lrm5OfLy8igtLSU/P5/CwkJycnKSy3aIyMZSWMiaWDiHori4mNHRUXbu3Mn4+DhvvfUWra2tOI7DtWvXkjWJrKwscnJyKCsrY9euXRw7diy57tPC9Z5EZOMoLGRNuK6LtRa/38/c3Bw+n49IJMKZM2dobGzkjTfeoKurK9lHkViSvLi4mOrqahoaGsjLy+PAgQM88MADanoS2WQUFrImEjWFaDSK4zj09vYyPDxMW1tbcq/sa9eu3XBeYWEhd999N/v376e4uJidO3dSUFCgGoXIJqOwkDWR2JMiMbGuqamJgYEB2tvb6enpYXh4+Lrjs7Oz8fv9hEIhjh8/TigUIjc3F7/fr9qEyCaksJD3LLVj++rVq3R3d9PR0UF/fz8dHR309fVdd7zf76e8vJyqqioOHjxIaWkpe/fu3YhLF5FVUljIexKJRLh8+TK5ubmMj4/z7rvv0tjYyNWrV+nt7WVwcPC6VWb9fj933HEHhYWFlJWVkZubmxw6mxg9NTw8rNVkRTYZhYXcMsdxePXVV+ntnV/Xsb+/n5aWFrq7u5meniYSiVwXFAUFBTQ0NPDQQw8xPj5OQ0NDcmY2zG9sVFBQkJyDobAQ2TwUFnJTUpuc+vv7k6vHXrlyhba2NsbGxmhububatWvXdWjn5ORw4sQJTp06xT333IPX671uZjaQXI1Wq8mKbD4KC1m1hXMp8vLyaGpqYnp6moGBAYaGhhgYGMBxHGKxGABer5esrCwCgQDZ2dl4vV6stRw6dAjguuGxu3fvTtYoFBYim4tnoy9AMsdicyn27t1LMBgkKyuL0dFRent7k0EB87WF8vJydu/eTSgUorq6GmstjuPcEAqJ2doKCpHNRzULWdbCJTxmZmZwXZeRkRE6Ozvp7Ozk3Xff5bXXXmNwcPC6c7OysigtLaW2tpaqqirKy8vxer0MDQ1RXFzMxMSEOrJFMoTCQpaUaHZyHIfOzk5qamowxjA9Pc21a9d47bXX6Ozs5OLFizcEBXDd8NiDBw9SXFxMTU0NrutSUFCQDCKFhcjmp7CQJbmuSzQaZXBwMDlXIi8vj+npacbHx7lw4QKXL1/GcZwbzvV6vZSWlnL8+PFkU1UgEGD37t20t7czMjKiJieRDKKwkCX5fD5isRiTk5OMj48nNzAqLi7mn/7pn2hsbGRycvKG88rKyigqKqKysvK6WdqJTY527dpFKBRSWIhkEIWFLCrRRFRTU8PQ0BDWWgYGBhgcHKSzs5N33nlnyaAoLy+npKSEO+64g7q6OmpqasjOzlZIiGQwhYXcIHWI7Pj4OFNTU0xMTNDT00NnZyfNzc2MjIzccF5NTQ233347x44do6ioiN27d1NQUMDRo0eTI58UEiKZSWEhN0jtq7h06RLj4+O0tLRw8eJFrl69uujqsYlVZwOBAHfddRdHjsxvma4VZEW2BoWFLGp8fJzu7u7k2k8///nPaWpqwlp7w7E7d+7k6NGjFBUVceTIEaLRKMPDw9TX16smIbJFKCwkyXXd5DDZ0dFRXnvtNc6fP09HRwc9PT2LnhMIBCgrKyMnJ4eqqir27dtHUVGRhsWKbDEKCwHmg+LSpUtcvHiRt99+G9d1+fnPf87Vq1cZGBhY9Jz6+np27NhBcXExhw4d4sCBA8zOzjI9PU0wGFRQiGwhCgsBoK+vj+eff54XX3yR5uZmJicncV2XiYmJRY8vLy9PDo0NBAKUl5ezf/9+6urqtGyHyBaksBBc1+Xq1avJXe3Gx8cZHR1d8viKigruuece9u7dS3Z2Nvv27WP//v0cO3ZMndkiW5TCQnAcJ7k6bFdX16LzJxKqq6vZu3cvH/zgBzlx4gRVVVXJEU+qSYhsXQqLbc5xHM6ePcsbb7zB66+/vmxQFBcXc+zYMU6cOMHu3bs5ePBgci8KEdnaFBbbVGLk06uvvsrf/M3f8O6773L58uUlj0/MzM7Pz6eyshKPx0NnZ6dqFCLbhMJiG3Ich8bGRrq6unjqqac4d+4cra2tzMzMLHp8fn4+OTk5lJWV4fV6GR0d5eDBg3i9Xg2PFdkmFBbbjOu6nDlzhtdee41z585x9uxZWltbmZubW/T4wsJCTp48yezsLIcPH6a6uppQKER+fj5+v19BIbJNKCy2Edd16ezs5JVXXuHcuXOcO3eOgYGBJYMiPz+furo6jh49Sn5+PoWFhezbt4/6+npqamrUBCWyjSgstonE4oCXLl3izTffpKenh46OjiWbnrKzs9mxYwd33XUXt912G+973/uYmZmhsrKSyspKhYTINqOw2CYcx6G9vZ3z58/T2dlJV1fXkkHh8Xiora3l6NGj/PIv/zJZWVkAlJSUKChEtimFxRbnui6RSIS33nqLv/iLv+Cdd95J7k+xmPLycgoLCzl69CiHDh2itLSU/Px87UMhss0pLLYwx3E4f/48586d44UXXuC1115bdh7Fnj17qKmpobKykhMnTnDbbbdRU1OjdZ5ERGGxVbmuS2NjIy+++CKvvvoqZ8+eXTYoqqqqeP/738/evXupq6tLrvmkoBARUFhsSY7j0NrayoULF7hw4QKtra0MDw8veXwwGORTn/oU1dXVHDp0iIqKCowxHDp0SEEhIoDCYstxHIfnn38+WZvo6elZci8KgF27drFnzx6qqqooLy/n9ttvT/ZNKChEJEFhscX09fXx5ptv8vrrr9Pd3c3g4OCSo54CgQAf/ehHqa2t5a677uLw4cNa60lEFuVJ1xsbY75jjBkwxlxIKQsaY14wxjTHf5ekvPY1Y0yLMeaKMeZjKeUnjDHvxl/7hjHGpOuaM10kEuHNN9/ktddeo7Gxkd7e3iWDwufzEQwGMcZQWlpKfX29gkJElpS2sAC+CzywoOyrwEvW2nrgpfhzjDENwCPA4fg53zTGZMXPeQx4FKiP/yx8z20tMTT29ddf54//+I/59re/TWNj45KbFsH8PIqTJ09y4sQJDh48yPvf/34FhYgsK23NUNbaV4wxexYUPwzcF3/8BPAy8JV4+Y+stVHgqjGmBThljGkHiqy1rwMYY74HfAp4Ll3XnUlc16W5uZm33nqLv/qrv6KxsZGhoaEVz7v99ts5efIkVVVV3H333Rw7dkz9EyKyrPXus6iw1oYBrLVhY8yOePku4I2U47rjZTPxxwvLhfmwmJyc5Pz583R3d68qKA4cOMCnP/1pjhw5wuHDhzUjW0RWZbN0cC/WD2GXKV/8TYx5lPkmK2pqatbmyjYx13U5ffo0P/vZz+jo6Fjx+L179/Lwww8naxPaAlVEVmu9w6LfGBOK1ypCwEC8vBuoTjmuCuiNl1ctUr4oa+3jwOMAJ0+eXDJUMpXrusn9IxzH4Yc//CHf/e53aWlpWXLlWJjf4e79738/J0+e5AMf+AAHDx5UUIjITVnvsHgG+DzwJ/HfP0kp/6Ex5s+Ancx3ZL9lrZ01xkwYY+4A3gQ+B/z5Ol/zppBYNTYajTI8PMybb77JE088wdWrV5cc8QSQlZVFRUUFhw8fprS0lOrqagWFiNy0tIWFMea/Md+ZXWaM6Qb+iPmQeNIY80WgE/hVAGvtRWPMk8AlIAZ8yVo7G3+r32J+ZFUu8x3b27Jz23VdotEoly9f5tlnn+Uf//EfaW9vX/G82tpa7rzzTh588EGOHTumPgoRuSXpHA31L5Z46SNLHP914OuLlJ8GjqzhpWWstrY2fvzjH/PWW2+tGBQ5OTnk5+dz8OBB9uzZw7Fjx7ZFP46IpMdm6eCWZSSaoJqamnj77bdXDIqSkhIOHjzI4cOH+chHPsLevXuprKxcn4sVkS1JYbHJRSIR2traaG9v5+LFi1y9enXZ47OysggEAtx333089NBD7N27VyvHish7prDYxCKRCD/4wQ+4evUq//RP/8TFixeXPT47O5vc3Fzq6uoIBALk5+crKERkTSgsNrGuri7Onj3LhQsXOH369LLH+nw+Tpw4gcfj4SMf+QgFBQXa2U5E1ozCYhNyXRfHcbhy5QqvvPIKbW1tK56TGPU0ODhIMBgkOztbQ2RFZM0oLDYZ13U5d+4cL7/8Ms8+++yKndler5e6ujp+53d+hzvuuIPJyUm8Xm+yCUpEZC0oLDaZlpYW/vIv/5LnnnuOrq6uZY8NhULcc889nDhxgnvuuYeGhgaA5CxvNUGJyFpRWGwSjuPQ1dXF97//fZ5++mkGBweXPb60tJTPfvazfPKTn2T//v3XdWQrJERkrSksNgHHcXjllVd4++23VxUUWVlZ3HbbbezcuZOioiKNeBKRtFNYbDDXdWltbeWFF17gySefpLd3yXUSk+rr6/ngBz/IRz/6UXw+X7LZSUQkXRQWG8hxHM6cOcPzzz/PE088wcjIyIrnlJaWcvLkSe644w58Ph/GGAWFiKSdwmKDuK7LG2+8wRNPPMELL7ywYlAUFhayf/9+7r//fo4fP86RI0eS8ygUFiKSbgqLDRKJRHj++ef5u7/7OxzHWfH4j3/843zmM5+hqqpKM7NFZN0pLDaA4zi8/fbbPPvss6sKinvvvZeHHnqIu+66S7UJEdkQCot1FolE+N73vsdjjz1GU1PTssfm5eXx4IMP8sgjj3Do0CHVJkRkwygs1pHjOHz/+9/nj/7ojxgfH1/22EAgwK/92q/x6U9/mqNHj2qdJxHZUAqLdeK6Lj/4wQ/4N//m36wYFCUlJfzhH/4hd999N4cOHdIaTyKy4RQWC7iuu2bLZSQWBAT46U9/yu/93u8xNTW17DmVlZU8+uijnDp1SkEhIpuGwiJFYkc6ay3GGHbv3n3LgeG6Ls3NzbS3t9PY2Mi/+3f/bsWgOHz4ML/yK7/CBz/4QYqKim7pc0VE0kFhkcJ1Xay1FBQU4DjOe5oZ7bouvb29vPjii/z1X/81ExMTyx6fm5vLF77wBe644w6Kiorw+/3qoxCRTUNhkSIxI9pxnPc8MzoSifCtb32Lv/3bv13xWGMMn/nMZ/jsZz9LMBjUqrEisukoLFL4fD527979nr+sHcfhz//8z1cVFPn5+fzKr/wKX/7yl6msrExeh4jIZqKwWOBWQyLRMe44Dt/4xjf4j//xP654TlFREb/+67/OvffeS01Nza1crojIulBYrIFEx3hPTw/f/OY3+fGPf7ziOXV1dfz7f//vOXbsGNbadbhKEZFbp7BYA67rMjIywje+8Q2efvrpFY9vaGjgxz/+MdnZ2cmRV2p6EpHNTGGxBnw+Hy+++OKqgqKsrIw//dM/paGhYU3ndIiIpJPCYg08//zz/MEf/MGKxwUCAf71v/7X3HPPPcCt94+IiKw3hcV79Nxzz/Hwww+veFwoFOLxxx/nvvvu06xsEck4CotbkGg++od/+IdVBUVxcTGPPfZYchtUEZFMo7C4SYllPF5//XV+8zd/c8XjKyoq+Pa3v83HPvYxBYWIZCzPRl9ApnEch1deeWVVQXHw4EG+9a1vKShEJOMpLBZITKxzXXfR15577jl++7d/e8X3+fjHP86PfvQjHnzwQQWFiGQ8NUOlWG7VWdd1+d73vreqGsXDDz/Md77zHYLBYLovWURkXahmkSIxuW5oaIiRkZHrahfPP//8qoLin//zf66gEJEtRzWLFI7j8NRTTzE2NkZxcTG/8zu/A8CZM2dWNerpy1/+Ml/5ylcUFCKy5SgsUjQ1NfHaa69hjGFubo5jx45RWFjIpz71qRXP/drXvsbv//7vaw6FiGxJCosU3d3ddHR04PP5cF2XS5cu8cd//McrnvfYY4/xhS98QR3ZIrJlKSxSzM3NMTk5mdw3ezVB8dRTT/HpT3863ZcmIrKhFBYpxsbGGB8fZ3Z2dlXH//SnP+XBBx9M81WJiGy8jBkNZYx5wBhzxRjTYoz5ajo+4+LFi6sOiubmZgWFiGwbGVGzMMZkAf8FuB/oBt42xjxjrb20lp/z5ptvruq4iYkJdWSLyLaSKTWLU0CLtbbNWusCPwJWHst6k86ePbviMZcuXcJxHPr6+pJ9GyIiW11G1CyAXUBXyvNu4APrfREf+tCHePrpp5mbm+Pw4cMUFBRw5513qpYhIlteptQszCJlN2xcbYx51Bhz2hhzenBwcE0vIBgMcu3aNVzXJRqNkp+fz9zcnGoXIrItZEpYdAPVKc+rgN6FB1lrH7fWnrTWniwvL1/TC/B4POTm5uLz+fD7/UxOTuLxeFSrEJFtIVOaod4G6o0xtUAP8Ajw2bX+kMQCggl+v597772XkydPcvz4cerr6wkGg8nJdwUFBQoLEdkWMiIsrLUxY8zvAD8DsoDvWGsvpumz0vG2IiIZLSPCAsBa+1Pgpxt9HSIi21Gm9FmIiMgGUliIiMiKFBYiIrIihYWIiKxIYSEiIitSWIiIyIoUFiIisiKFhYiIrEhhISIiK1JYiIjIisxWXQvJGDMIdNzi6WXA0BpezkbZCvexFe4BtsZ9bIV7AN3HSnZba29YtnvLhsV7YYw5ba09udHX8V5thfvYCvcAW+M+tsI9gO7jVqkZSkREVqSwEBGRFSksFvf4Rl/AGtkK97EV7gG2xn1shXsA3cctUZ+FiIisSDULERFZkcIihTHmAWPMFWNMizHmqxt9PQsZY75jjBkwxlxIKQsaY14wxjTHf5ekvPa1+L1cMcZ8LKX8hDHm3fhr3zCpG4+n/x6qjTH/0xjTaIy5aIz53Qy9jxxjzFvGmHPx+/h/MvE+4p+fZYz5uTHm2Qy+h/b4579jjDmdwfcRMMY8ZYy5HP83cuemuQ9rrX7mm+KygFagDvAB54CGjb6uBdf4z4DjwIWUsv8X+Gr88VeB/xB/3BC/Bz9QG7+3rPhrbwF3AgZ4DnhwHe8hBByPPy4EmuLXmmn3YYCC+ONs4E3gjky7j/jn/x7wQ+DZTPxvKv757UDZgrJMvI8ngP8t/tgHBDbLfazb/wib/Sf+P+zPUp5/DfjaRl/XIte5h+vD4goQij8OAVcWu37gZ/F7DAGXU8r/BfCtDbyfnwD3Z/J9AHnAWeADmXYfQBXwEvBhfhEWGXUP8c9s58awyKj7AIqAq8T7kjfbfagZ6hd2AV0pz7vjZZtdhbU2DBD/vSNevtT97Io/Xli+7owxe4DbmP+rPOPuI9588w4wALxgrc3E+/hPwJeBuZSyTLsHAAv8vTHmjDHm0XhZpt1HHTAI/FW8WfAvjTH5bJL7UFj8wmJtepk8VGyp+9kU92mMKQD+FvhX1trx5Q5dpGxT3Ie1dtZa+37m/zo/ZYw5sszhm+4+jDGfAAastWdWe8oiZZvi/wvgbmvtceBB4EvGmH+2zLGb9T68zDczP2atvQ2YZL7ZaSnreh8Ki1/oBqpTnlcBvRt0LTej3xgTAoj/HoiXL3U/3fHHC8vXjTEmm/mg+IG19r/HizPuPhKstaPAy8ADZNZ93A180hjTDvwI+LAx5vtk1j0AYK3tjf8eAJ4GTpF599ENdMdrqABPMR8em+I+FBa/8DZQb4ypNcb4gEeAZzb4mlbjGeDz8cefZ74PIFH+iDHGb4ypBeqBt+LV2AljzB3xERKfSzkn7eKf+V+BRmvtn6W8lGn3UW6MCcQf5wK/BFzOpPuw1n7NWltlrd3D/H/v/2Ct/ZeZdA8Axph8Y0xh4jHwUeBCpt2HtbYP6DLGHIgXfQS4tGnuYz07oTb7D/Bx5kfntAJ/sNHXs8j1/TcgDMww/9fDF4FS5jsom+O/gynH/0H8Xq6QMhoCOMn8P6ZW4C9Y0KGW5nu4h/kq8XngnfjPxzPwPo4BP4/fxwXgD+PlGXUfKddwH7/o4M6oe2C+rf9c/Odi4t9upt1H/PPfD5yO/3f1d0DJZrkPzeAWEZEVqRlKRERWpLAQEZEVKSxERGRFCgsREVmRwkJERFaksBARkRUpLETWmDHGSXkcSiz9vcpzX0xdglpks1BYiKTX7wHfXlhojPm3xphfX+T4vwZ+O90XJXKzFBYi6fVp4PmbOP4Z5peUFtlUFBYiaRJfr2fEWhtd7TnW2hHAb4wpTd+Vidw870ZfgMgWFmJ+fwIAjDFHmW9mAqgEXGPMv4o//4i1djj+eADYCSSei2w4hYVI+kwDOYkn1tp3mV8oDmPMvwXarbXfXeS8nPi5IpuGmqFE0qeJ+W1wVy2+pHQl89uEimwaCguRNLHWTgKtxph9N3HaCeANa20sTZclcku0RLlIGhlj/hfghLX2/17l8f8ZeMZa+1J6r0zk5qjPQiSNrLVP3+TIpgsKCtmMVLMQEZEVqc9CRERWpLAQEZEVKSxERGRFCgsREVmRwkJERFb0/wP6vqhbFMR1KwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6, 6))\n",
    "plt.plot(refltable['I(+)'].to_numpy(), refltable['I(-)'].to_numpy(), 'k.', alpha=0.1)\n",
    "plt.xlabel(\"I(+)\")\n",
    "plt.ylabel(\"I(-)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [next example](2_mergingstats.ipynb), we will investigate this anomalous signal in more detail.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Writing Reflection Data\n",
    "\n",
    "It is also possible to write out MTZ files using `DataSet.write_mtz()`. This functionality depends on the correct setting of each column's `dtype`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "refltable.write_mtz(\"data/HEWL_SSAD_24IDC.mtz\")"
   ]
  }
 ],
 "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": 2
}