{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced McStas features:  SPLIT\n",
    "McStas uses the Monte Carlo ray-tracing technique, which allows some tricks on how the physics is sampled as long as the resulting probability distributions match the physics. This is possible as each ray has a weight, corresponding to how much intensity this ray represents. The SPLIT keyword can be used to split a ray into many equal parts, which can be useful if the remaining instrument has many different simulated and random outcomes. In this tutorial we will use the SPLIT keyword on a powder sample, as there are many powder Bragg peaks each ray could select, and splitting the ray samples this more efficiently."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up an example instrument\n",
    "First we set up an example instrument, this is taken from the basic tutorial and consists of source, guide, monochromator, sample and banana detector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mcstasscript as ms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "instrument = ms.McStas_instr(\"python_tutorial\", input_path=\"run_folder\",\n",
    "                                output_path=\"data_folder/mcstas_SPLIT\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "instrument.add_component(\"Origin\", \"Progress_bar\")\n",
    "\n",
    "src = instrument.add_component(\"source\", \"Source_div\")\n",
    "src.xwidth = 0.1\n",
    "src.yheight = 0.05\n",
    "src.focus_aw = 1.2\n",
    "src.focus_ah = 2.3\n",
    "\n",
    "wavelength = instrument.add_parameter(\"wavelength\", value=5.0, comment=\"Wavelength in [Ang]\")\n",
    "src.lambda0 = wavelength\n",
    "src.dlambda = \"0.03*wavelength\"\n",
    "\n",
    "guide = instrument.add_component(\"guide\", \"Guide_gravity\", AT=[0,0,2], RELATIVE=src)\n",
    "guide.w1 = 0.05\n",
    "guide.w2 = 0.05\n",
    "guide.h1 = 0.05\n",
    "guide.h2 = 0.05\n",
    "guide.l = 8.0\n",
    "guide.m = 3.5\n",
    "guide.G = -9.82\n",
    "\n",
    "mono_Q = instrument.add_declare_var(\"double\", \"mono_Q\", value=1.714) # Q for Ge 311\n",
    "instrument.add_declare_var(\"double\", \"wavevector\")\n",
    "instrument.append_initialize(\"wavevector = 2.0*PI/wavelength;\")\n",
    "\n",
    "mono_rotation = instrument.add_declare_var(\"double\", \"mono_rotation\")\n",
    "instrument.append_initialize(\"mono_rotation = asin(mono_Q/(2.0*wavevector))*RAD2DEG;\")\n",
    "instrument.append_initialize('printf(\"monochromator rotation = %g deg\\\\n\", mono_rotation);')\n",
    "\n",
    "mono = instrument.add_component(\"mono\", \"Monochromator_flat\")\n",
    "mono.zwidth = 0.05\n",
    "mono.yheight = 0.08\n",
    "mono.Q = mono_Q\n",
    "mono.set_AT([0, 0, 8.5], RELATIVE=guide)\n",
    "mono.set_ROTATED([0, mono_rotation, 0], RELATIVE=guide)\n",
    "\n",
    "beam_direction = instrument.add_component(\"beam_dir\", \"Arm\", AT_RELATIVE=mono)\n",
    "beam_direction.set_ROTATED([0, \"mono_rotation\", 0], RELATIVE=\"mono\")\n",
    "\n",
    "sample = instrument.add_component(\"sample\", \"PowderN\", AT=[0,0,1.1], RELATIVE=beam_direction)\n",
    "sample.radius = 0.015\n",
    "sample.yheight = 0.05\n",
    "sample.reflections = '\"Na2Ca3Al2F14.laz\"'\n",
    "\n",
    "banana = instrument.add_component(\"banana\", \"Monitor_nD\", RELATIVE=sample)\n",
    "banana.xwidth = 2.0\n",
    "banana.yheight = 0.3\n",
    "banana.restore_neutron = 1\n",
    "banana.filename = '\"banana.dat\"'\n",
    "banana.options = '\"theta limits=[5 175] bins=150, banana\"'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running the simulation\n",
    "Here we run the simulation with very few neutrons to show problematic sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "instrument.settings(ncount=1E6)\n",
    "\n",
    "instrument.set_parameters(wavelength=2.8)\n",
    "                          \n",
    "data_low = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.make_sub_plot(data_low)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding the SPLIT keyword\n",
    "Here we add the SPLIT keyword to the sample, we choose to split each ray into 30."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample.set_SPLIT(30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "# No need to set settings or parameters as these have not changed\n",
    "data_reasonable = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.make_sub_plot(data_reasonable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Caution on split value\n",
    "It is however possible to mismanage splitting, mainly by simulating a too few rays and splitting too much. Here we do this on purpose to see how such data would look. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "sample.set_SPLIT(10000)\n",
    "\n",
    "instrument.settings(ncount=1E3) # Change settings to lower ncount, but keep parameters\n",
    "                          \n",
    "data_unreasonable = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.make_sub_plot(data_unreasonable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison with high statistics run\n",
    "We here compare the different runs to a reference. The reference run is set up to have 50 times more rays than the earlier runs with 5E7 instead of 1E6 rays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "sample.set_SPLIT(1)\n",
    "instrument.settings(ncount=2E7)\n",
    "                          \n",
    "data_ref = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.make_sub_plot(data_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting data on same plot\n",
    "Here we only have one monitor in each data list, but we still use the *name_search* function to retrieve the correct data object from each. This avoids the code breaking in case additional monitors are added.\n",
    "\n",
    "Once we have the objects, we use the *xaxis*, *Intensity* and *Error* attributes to plot the data with matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "banana_low = ms.name_search(\"banana\", data_low)\n",
    "banana_reasonable = ms.name_search(\"banana\", data_reasonable)\n",
    "banana_unreasonable = ms.name_search(\"banana\", data_unreasonable)\n",
    "banana_ref = ms.name_search(\"banana\", data_ref)\n",
    "\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.errorbar(banana_low.xaxis, banana_low.Intensity, yerr=banana_low.Error, fmt=\"r\")\n",
    "plt.errorbar(banana_ref.xaxis, banana_ref.Intensity, yerr=banana_ref.Error, fmt=\"b\")\n",
    "plt.xlabel(\"2Theta [deg]\")\n",
    "plt.ylabel(\"Intensity [n/s]\")\n",
    "plt.legend([\"Low statistics\", \"High statistics reference\"])\n",
    "\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.errorbar(banana_reasonable.xaxis, banana_reasonable.Intensity,\n",
    "             yerr=banana_reasonable.Error, fmt=\"r\")\n",
    "plt.errorbar(banana_ref.xaxis, banana_ref.Intensity, yerr=banana_ref.Error, fmt=\"b\")\n",
    "plt.xlabel(\"2Theta [deg]\")\n",
    "plt.ylabel(\"Intensity [n/s]\")\n",
    "plt.legend([\"Low statistics with SPLIT\", \"High statistics reference\"])\n",
    "\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.errorbar(banana_unreasonable.xaxis, banana_unreasonable.Intensity,\n",
    "             yerr=banana_unreasonable.Error, fmt=\"r\")\n",
    "plt.errorbar(banana_ref.xaxis, banana_ref.Intensity, yerr=banana_ref.Error, fmt=\"b\")\n",
    "plt.xlabel(\"2Theta [deg]\")\n",
    "plt.ylabel(\"Intensity [n/s]\")\n",
    "l = plt.legend([\"Very low statistics with unreasonable SPLIT\", \"High statistics reference\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretation of the data\n",
    "We see that with low statistics, the data quality is so bad that noise can be mistaken for peaks. Using SPLIT improves the situation a lot, and the data is very similar to the high statistics reference which takes longer to compute. The situation with a low number of simulated rays and very high SPLIT have some erratic behavior, showing some very different peak intensities than the reference, and some peaks that shouldn't be there at all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "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.8"
  },
  "metadata": {
   "execution": {
    "timeout": 100
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}