{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Union tagging system\n",
    "The *Union_master* is capable of recording histories for each neutron in a tree like fashion and add the total intensities for each unique history together. At the end of the simulation these are sorted by intensity and written to file, and the top 20 are shown in the terminal. This system does not work with MPI, if MPI is used only part of the data is written to disk. The system can take up a large amount of memory when used, so it is disabled per default.\n",
    "\n",
    "First we set up a simple instrument with sample, container and a layer of cryostat."
   ]
  },
  {
   "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\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Al_inc = instrument.add_component(\"Al_inc\", \"Incoherent_process\")\n",
    "Al_inc.sigma = 0.0082\n",
    "Al_inc.unit_cell_volume = 66.4\n",
    "\n",
    "Al_pow = instrument.add_component(\"Al_pow\", \"Powder_process\")\n",
    "Al_pow.reflections = '\"Al.laz\"'\n",
    "\n",
    "Al = instrument.add_component(\"Al\", \"Union_make_material\")\n",
    "Al.process_string = '\"Al_inc,Al_pow\"'\n",
    "Al.my_absorption = 100*0.231/66.4 # barns [m^2 E-28]*Å^3 [m^3 E-30]=[m E-2], factor 100\n",
    "\n",
    "Sample_inc = instrument.add_component(\"Sample_inc\", \"Incoherent_process\")\n",
    "Sample_inc.sigma = 3.4176\n",
    "Sample_inc.unit_cell_volume = 1079.1\n",
    "\n",
    "Sample_pow = instrument.add_component(\"Sample_pow\", \"Powder_process\")\n",
    "Sample_pow.reflections = '\"Na2Ca3Al2F14.laz\"'\n",
    "\n",
    "Sample = instrument.add_component(\"Sample\", \"Union_make_material\")\n",
    "Sample.process_string = '\"Sample_inc,Sample_pow\"'\n",
    "Sample.my_absorption = 100*2.9464/1079.1\n",
    "\n",
    "src = instrument.add_component(\"source\", \"Source_div\")\n",
    "src.xwidth = 0.01\n",
    "src.yheight = 0.035\n",
    "src.focus_aw = 0.01\n",
    "src.focus_ah = 0.01\n",
    "src.lambda0 = instrument.add_parameter(\"wavelength\", value=5.0,\n",
    "                                       comment=\"Wavelength in [Ang]\")\n",
    "src.dlambda = \"0.01*wavelength\"\n",
    "src.flux = 1E13\n",
    "\n",
    "sample_geometry = instrument.add_component(\"sample_geometry\", \"Union_cylinder\")\n",
    "sample_geometry.yheight = 0.03\n",
    "sample_geometry.radius = 0.0075\n",
    "sample_geometry.material_string='\"Sample\"' \n",
    "sample_geometry.priority = 100\n",
    "sample_geometry.set_AT([0,0,1], RELATIVE=src)\n",
    "\n",
    "container = instrument.add_component(\"sample_container\", \"Union_cylinder\")\n",
    "container.set_RELATIVE(sample_geometry)\n",
    "container.yheight = 0.03+0.003 # 1.5 mm top and button\n",
    "container.radius = 0.0075 + 0.0015 # 1.5 mm sides of container\n",
    "container.material_string='\"Al\"' \n",
    "container.priority = 99\n",
    "\n",
    "inner_wall = instrument.add_component(\"cryostat_wall\", \"Union_cylinder\")\n",
    "inner_wall.set_AT([0,0,0], RELATIVE=sample_geometry)\n",
    "inner_wall.yheight = 0.12\n",
    "inner_wall.radius = 0.03\n",
    "inner_wall.material_string='\"Al\"' \n",
    "inner_wall.priority = 80\n",
    "\n",
    "inner_wall_vac = instrument.add_component(\"cryostat_wall_vacuum\", \"Union_cylinder\")\n",
    "inner_wall_vac.set_AT([0,0,0], RELATIVE=sample_geometry)\n",
    "inner_wall_vac.yheight = 0.12 - 0.008\n",
    "inner_wall_vac.radius = 0.03 - 0.002\n",
    "inner_wall_vac.material_string='\"Vacuum\"' \n",
    "inner_wall_vac.priority = 81\n",
    "\n",
    "logger_zx = instrument.add_component(\"logger_space_zx\", \"Union_logger_2D_space\")\n",
    "logger_zx.set_RELATIVE(sample_geometry)\n",
    "logger_zx.D_direction_1 = '\"z\"'\n",
    "logger_zx.D1_min = -0.04\n",
    "logger_zx.D1_max = 0.04\n",
    "logger_zx.n1 = 300\n",
    "logger_zx.D_direction_2 = '\"x\"'\n",
    "logger_zx.D2_min = -0.04\n",
    "logger_zx.D2_max = 0.04\n",
    "logger_zx.n2 = 300\n",
    "logger_zx.filename = '\"logger_zx.dat\"'\n",
    "\n",
    "logger_zy = instrument.add_component(\"logger_space_zy\", \"Union_logger_2D_space\")\n",
    "logger_zy.set_RELATIVE(sample_geometry)\n",
    "logger_zy.D_direction_1 = '\"z\"'\n",
    "logger_zy.D1_min = -0.04\n",
    "logger_zy.D1_max = 0.04\n",
    "logger_zy.n1 = 300\n",
    "logger_zy.D_direction_2 = '\"y\"'\n",
    "logger_zy.D2_min = -0.06\n",
    "logger_zy.D2_max = 0.06\n",
    "logger_zy.n2 = 300\n",
    "logger_zy.filename = '\"logger_zy.dat\"'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding the Union_master with tagging\n",
    "There are two important parameters to consider when setting a *Union_master* up for tagging:\n",
    "- enable_tagging [default 0] 0 for disable, 1 for enable\n",
    "- history_limit [default 300000] Limit of unique histories recorded\n",
    "\n",
    "As the *Union_master* component records each ray in succession, their unique history is added to the history tree. If a ray takes the same path in the tree, the intensity gets added to that unique history. When the history limit is reached, the tree is not built out further, but if already recorded histories occur, their intensity is still added to the existing tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "master = instrument.add_component(\"master\", \"Union_master\")\n",
    "master.enable_tagging = 1\n",
    "master.history_limit = 300000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "instrument.set_parameters(wavelength=3.0)\n",
    "instrument.settings(ncount=3E5, output_path=\"data_folder/union_tagging\")\n",
    "\n",
    "data = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.name_plot_options(\"logger_space_zx\", data, log=True, orders_of_mag=4)\n",
    "ms.name_plot_options(\"logger_space_zy\", data, log=True, orders_of_mag=4)\n",
    "ms.make_sub_plot(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding the history file\n",
    "A file called Union_history.dat is written in the run folder with all the unique histories. In some cases a bug happens that causes the file to be unreadable, for now the best fix is to rerun the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "with open(\"run_folder/union_history.dat\", \"r\", errors=\"ignore\") as file:\n",
    "    instrument_string = file.read()\n",
    "    print(instrument_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting the histories\n",
    "The history with highest intensity in my run is: V0 -> V3 -> V4 -> V2 -> V1 -> V2 -> V4 -> V3 -> V0, which means:\n",
    "- Neutron entered Volume0 (Surrounding vacuum)\n",
    "- Neutron entered Volume3 (Cryostat wall)\n",
    "- Neutron entered Volume4 (Cryostat vacuum)\n",
    "- Neutron entered Volume2 (Container)\n",
    "- Neutron entered Volume1 (Sample)\n",
    "- Neutron entered Volume2 (Container)\n",
    "- Neutron entered Volume4 (Cryostat vacuum)\n",
    "- Neutron entered Volume3 (Cryostat wall)\n",
    "- Neutron entered Volume0 (Surrounding vacuum)\n",
    "\n",
    "So the most likely occurrence is that the ray is propagated through all geometries. The next most likely history is:\n",
    "\n",
    "V0 -> V3 -> V4 -> V2 -> V1 -> P1 -> V2 -> V4 -> V3 -> V0\n",
    "\n",
    "- Neutron entered Volume0 (Surrounding vacuum)\n",
    "- Neutron entered Volume3 (Cryostat wall)\n",
    "- Neutron entered Volume4 (Cryostat vacuum)\n",
    "- Neutron entered Volume2 (Container)\n",
    "- Neutron entered Volume1 (Sample)\n",
    "- Neutron scattered on Process1 (Since we are in Volume1, that would be Sample_pow)\n",
    "- Neutron entered Volume2 (Container)\n",
    "- Neutron entered Volume4 (Cryostat vacuum)\n",
    "- Neutron entered Volume3 (Cryostat wall)\n",
    "- Neutron entered Volume0 (Surrounding vacuum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use conditional component to filter tagging\n",
    "Just like a logger can be modified by a conditional component to only record events that satisfy that condition, the tagging system of the *Union_master* component can also be modified by a conditional component. In that case, the tagging system will only record events that satisfy the condition imposed by the conditional component.\n",
    "\n",
    "This can be useful to explain an unexpected feature, as the conditional components can filter for energy, time, direction or any combination of these. Here we demonstrate this feature by adding a *Union_conditional_PSD* outside of the direct beam and enabling the *master_tagging* control parameter, which causes the condition to be applied to the tagging system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up an arm pointing to the relevant spot\n",
    "spot_dir = instrument.add_component(\"spot_dir\", \"Arm\", RELATIVE=sample_geometry,\n",
    "                                    before=\"master\")\n",
    "spot_dir.set_ROTATED([0, 60, 0], RELATIVE=sample_geometry)\n",
    "\n",
    "# Set up a conditional component targeting all our loggers\n",
    "PSD_conditional = instrument.add_component(\"space_all_PSD_conditional\", \"Union_conditional_PSD\",\n",
    "                                           before=\"master\")\n",
    "PSD_conditional.xwidth = 0.2\n",
    "PSD_conditional.yheight = 0.2\n",
    "PSD_conditional.master_tagging = 1\n",
    "PSD_conditional.set_AT([0, 0, 0.5], RELATIVE=spot_dir) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "data = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "with open(\"run_folder/union_history.dat\", \"r\", errors=\"ignore\") as file:\n",
    "    instrument_string = file.read()\n",
    "    print(instrument_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting the data\n",
    "Now all the histories contain a scattering process as this is necessary to reach the rectangle placed by the Union_conditional_PSD component set at 2theta = 60 deg. We also observe the intensities are much lower, simply because only a fraction of the simulated rays satisfy this condition."
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}