{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced geometry component concepts: Exit geometry and number of activations\n",
    "This notebook explains the concept of exit geometry and the activation counter both of which are tied to the geometry components and how they are treated by the *Union_master*.\n",
    "\n",
    "An exit geometry is created by setting the *material_string* of a geometry to \"Exit\", and if a ray enters such a geometry, it is immediately released from the master component. Normally this only happens when the ray does not intersect any geometries. There are several uses for this, for example inserting a monitor within a Union geometry ensemble."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mcstasscript as ms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up an example with empty sample container\n",
    "First we set up an example with an empty sample container."
   ]
  },
  {
   "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",
    "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_volume = instrument.add_component(\"sample_volume\", \"Union_cylinder\")\n",
    "sample_volume.yheight = 0.03\n",
    "sample_volume.radius = 0.0075\n",
    "sample_volume.material_string='\"Vacuum\"' \n",
    "sample_volume.priority = 100\n",
    "sample_volume.set_AT([0,0,1], RELATIVE=src)\n",
    "\n",
    "container = instrument.add_component(\"sample_container\", \"Union_cylinder\")\n",
    "container.set_RELATIVE(sample_volume)\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",
    "container_lid = instrument.add_component(\"sample_container_lid\", \"Union_cylinder\")\n",
    "container_lid.set_AT([0, 0.0155, 0], RELATIVE=container)\n",
    "container_lid.yheight = 0.004\n",
    "container_lid.radius = 0.013\n",
    "container_lid.material_string='\"Al\"' \n",
    "container_lid.priority = 98\n",
    "\n",
    "inner_wall = instrument.add_component(\"cryostat_wall\", \"Union_cylinder\")\n",
    "inner_wall.set_AT([0,0,0], RELATIVE=sample_volume)\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_volume)\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_volume)\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_volume)\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\"'\n",
    "\n",
    "master = instrument.add_component(\"master\", \"Union_master\")\n",
    "\n",
    "banana = instrument.add_component(\"banana\", \"Monitor_nD\", RELATIVE=sample_volume)\n",
    "banana.xwidth = 1.5\n",
    "banana.yheight = 0.4\n",
    "banana.restore_neutron = 1\n",
    "banana.options = '\"theta limits=[5 175] bins=250, banana\"'\n",
    "\n",
    "instrument.show_diagram()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see there is one material defined, *Al*, and 3 of the 5 geometries use this as their material. The 5 geometries are all picked up by the master component."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "instrument.set_parameters(wavelength=3.0)\n",
    "instrument.settings(ncount=3E6, output_path=\"data_folder/union_external\")\n",
    "\n",
    "data_empty = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.name_plot_options(\"logger_space_zx\", data_empty, log=True, orders_of_mag=4)\n",
    "ms.name_plot_options(\"logger_space_zy\", data_empty, log=True, orders_of_mag=4)\n",
    "ms.make_sub_plot(data_empty[0:2])\n",
    "ms.make_sub_plot(data_empty[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The loggers show our empty cryostat and sample container illuminated by a beam, and the banana monitor contains some powder scattering from the aluminium. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding an exit volume\n",
    "Now we switch the sample_volume material from Vacuum to exit, ejecting rays from the simulation when they encounter it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_volume.material_string='\"Exit\"'\n",
    "instrument.show_diagram()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diagram does not change between *Vacuum* and *Exit* being used for the *sample_volume*. Let us run the simulation and see the difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "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[0:2])\n",
    "ms.make_sub_plot(data[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the logger data one can see the back part of the sample container and exit of the cryostat is no longer illuminated, as the rays are removed from the Union simulation as soon as the touch the volume inside the sample container. The removed rays do not disappear completely, they are delivered to the next McStas component in the instrument."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding an external component in the gap\n",
    "We can now see any part of the beam that intersected the exit volume is basically removed from the Union simulation. It is now however possible to insert another component within that exit volume, for example a sample not available as a Union process. Here we just use a PowderN sample in order to demonstrate. We select the same dimensions as the exit volume, but subtract 10 micrometer to avoid a perfect overlap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = instrument.add_component(\"sample\", \"PowderN\", after=\"master\")\n",
    "sample.set_AT([0,0,0], RELATIVE=sample_volume)\n",
    "sample.radius = sample_volume.radius - 1E-5\n",
    "sample.yheight = sample_volume.yheight - 2E-5\n",
    "sample.reflections = '\"Na2Ca3Al2F14.laz\"'\n",
    "\n",
    "instrument.show_diagram()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running the simulation again\n",
    "We run the simulation again, but know that the scattering within the sample wont be directly visible in the loggers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "data_wrong = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.name_plot_options(\"logger_space_zx\", data_wrong, log=True, orders_of_mag=4)\n",
    "ms.name_plot_options(\"logger_space_zy\", data_wrong, log=True, orders_of_mag=4)\n",
    "ms.make_sub_plot(data_wrong[0:2])\n",
    "ms.make_sub_plot(data_wrong[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretation of the data\n",
    "Now we have added a sample inside the Union geometry, but when the neutron reaches that sample, it is ignored by the sample environment leading to unphysical behavior. Here the beam does not illuminate the sample environment on the way out, and all rays scattered by the PowderN sample are not attenuated by the walls of the cryostat when leaving."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Allowing the rays to return to the Union_master\n",
    "Now we could recreate the entire sample environment with new geometries and insert an additional unit master to grab the neutrons after the external sample, yet this would be error prone as all geometries would need to be exactly the same. Instead it is possible to tell Union geometries that they should be simulated in several of the next *Union_master* components using the *number_of_activation* parameter on each Union geometry, which is 1 per default.\n",
    "\n",
    "Setting it to 2, we tell the geometries that they should be simulated in the two next *Union_master* components. We do not update the sample_volume which is an exit volume, as this would allow the ray to escape once more. Instead we will replace it with Vacuum, but one could also have placed something closer to the actual sample.\n",
    "\n",
    "One last necessary detail is to set the *allow_inside_start* parameter on the second *Union_master* component. This disables an error message that would occur if a neutron starts inside a Union geometry, as this is most likely an error. Here we want to do this on purpose, and we need to let the *Union_master* component know this is allowed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "container.number_of_activations = 2\n",
    "container_lid.number_of_activations = 2\n",
    "inner_wall.number_of_activations = 2\n",
    "inner_wall_vac.number_of_activations = 2\n",
    "\n",
    "sample_replacement = instrument.add_component(\"sample_volume_replace\", \"Union_cylinder\", after=sample)\n",
    "sample_replacement.yheight = sample_volume.yheight\n",
    "sample_replacement.radius = sample_volume.radius\n",
    "sample_replacement.material_string='\"Vacuum\"' \n",
    "sample_replacement.priority = 101\n",
    "sample_replacement.set_AT([0,0,0], RELATIVE=sample_volume)\n",
    "\n",
    "master_2 = instrument.add_component(\"master_after_sample\", \"Union_master\", after=sample_replacement)\n",
    "master_2.allow_inside_start=1\n",
    "\n",
    "instrument.show_diagram()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that *sample_volume* is only simulated in the first master, there is no arrow to the second master as with all the other geometries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "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[0:2])\n",
    "ms.make_sub_plot(data[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretation of the data\n",
    "Now we see evidence of the beam leaving the sample environment after interacting with the sample, and also elevated scattering in comparison to the empty can. This is now a reasonable simulation containing an external component inside a Union geometry ensemble, but there is still one problem, if the ray leaves the external component and reenters later, it will find a Vacuum instead of that sample. This can be fixed to some extend by adding a second copy of the external component and a third *Union_master* component, while incrementing the *number_of_activations* on all geometries, then the ray would be able to leave and enter the external component once before the external component effectively disappears. Even with this assumption, it is still a reasonable approximation and a flexible approach to add for example a mirror inside a sample environment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of the three datasets\n",
    "Here we compare the three datasets, the empty sample environment, the wrong simulation where rays scattered in the sample could not interact with the sample environment, and the full simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "banana_empty = ms.name_search(\"banana\", data_empty)\n",
    "banana_wrong = ms.name_search(\"banana\", data_wrong)\n",
    "banana_sample = ms.name_search(\"banana\", data)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.plot(banana_empty.xaxis, banana_empty.Intensity, \"r\",\n",
    "         banana_wrong.xaxis, banana_wrong.Intensity, \"b\",\n",
    "         banana_sample.xaxis, banana_sample.Intensity, \"k\")\n",
    "plt.xlabel(\"2Theta [deg]\")\n",
    "plt.ylabel(\"Intensity [n/s]\")\n",
    "plt.legend([\"No sample\", \"Wrong simulation, no exit\", \"Full simulation\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretation of the data\n",
    "We see that the wrong simulation have slightly lower background, and more peak intensity. We also see the peak shape is different near the aluminium Bragg peaks. Since the Union components contain a powder process, one can also recreate this example without using an external PowderN component to check the accuracy of the approach. It is however not a fair comparison, as the Union powder process will perform multiple scattering where PowderN will not. "
   ]
  },
  {
   "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
}