{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing what happens in Union master\n",
    "One disadvantage to collecting all the simulation in the Union_master component, is that it is not possible to insert monitors between the parts to check on the beam. This issue is addressed by adding logger components that can record scattering and absorption events that occurs during the simulation. This notebook will show examples on the usage of loggers and their features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up materials and geometry to investigate\n",
    "First we set up the same mock cryostat we created in the advanced geometry tutorial to have an interesting system to investigate using the loggers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mcstasscript as ms\n",
    "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",
    "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_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",
    "outer_wall = instrument.add_component(\"outer_cryostat_wall\", \"Union_cylinder\")\n",
    "outer_wall.set_AT([0,0,0], RELATIVE=sample_geometry)\n",
    "outer_wall.yheight = 0.15\n",
    "outer_wall.radius = 0.1\n",
    "outer_wall.material_string='\"Al\"' \n",
    "outer_wall.priority = 60\n",
    "\n",
    "outer_wall_vac = instrument.add_component(\"outer_cryostat_wall_vacuum\", \"Union_cylinder\")\n",
    "outer_wall_vac.set_AT([0,0,0], RELATIVE=sample_geometry)\n",
    "outer_wall_vac.yheight = 0.15 - 0.01\n",
    "outer_wall_vac.radius = 0.1 - 0.003\n",
    "outer_wall_vac.material_string='\"Vacuum\"' \n",
    "outer_wall_vac.priority = 61\n",
    "\n",
    "instrument.show_components()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "instrument.available_components(\"union\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding Union logger components\n",
    "Union logger components need to be added before the *Union_master* component, as the master need to record the necessary information when the simulation is being performed. There are two different kind of Union logger components, the *loggers* that record scattering and the *abs_loggers* that record absorption. They have similar parameters and user interface. Here is a list of the currently available loggers:\n",
    "\n",
    "- Union_logger_1D\n",
    "- Union_logger_2D_space\n",
    "- Union_logger_2D_space_time\n",
    "- Union_logger_3D_space\n",
    "- Union_logger_2D_kf\n",
    "- Union_logger_2D_kf_time\n",
    "- Union_logger_2DQ\n",
    "\n",
    "- Union_abs_logger_1D_space\n",
    "- Union_abs_logger_1D_space_tof\n",
    "- Union_abs_logger_2D_space\n",
    "\n",
    "The most commonly used logger is probably the *Union_logger_2D_space*, this component records spatial distribution of scattering, here are the available parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "instrument.component_help(\"Union_logger_2D_space\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up a 2D_space logger\n",
    "One can select which two axis to record using *D_direction_1* and *D_direction_2*, and the range with for example *D1_min* and *D1_max*. When spatial information is recorded it is also important to place the logger at an appropriate position, here we center it on the sample position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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.12\n",
    "logger_zx.D1_max = 0.12\n",
    "logger_zx.n1 = 300\n",
    "logger_zx.D_direction_2 = '\"x\"'\n",
    "logger_zx.D2_min = -0.12\n",
    "logger_zx.D2_max = 0.12\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.12\n",
    "logger_zy.D1_max = 0.12\n",
    "logger_zy.n1 = 300\n",
    "logger_zy.D_direction_2 = '\"y\"'\n",
    "logger_zy.D2_min = -0.12\n",
    "logger_zy.D2_max = 0.12\n",
    "logger_zy.n2 = 300\n",
    "logger_zy.filename = '\"logger_zy.dat\"'\n",
    "\n",
    "master = instrument.add_component(\"master\", \"Union_master\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running the simulation\n",
    "If mpi is installed, one can add mpi=N where N is the number of cores available to speed up the simulation."
   ]
  },
  {
   "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_loggers\")\n",
    "\n",
    "data = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.make_sub_plot(data, log=True, orders_of_mag=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting the data\n",
    "The zx logger views the cryostat from the top, while the zy loggers shows it from the side. These are histograms of scattered intensity, and it is clear the majority of the scattering happens in the direct beam. There are however scattering events in all parts of our mock cryostat, as neutrons that scattered in either the sample or cryostat walls could go in any direction due to the incoherent scattering. The aluminium and sample also have powder scattering, so some patterns can be seen from the debye scherrer cones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logger targets\n",
    "It is possible to attach a logger to a certain geometry, or even a list of geometries using the *target_geometry* parameter. In that way one can for example view the scattering in the sample environment, while ignoring the sample. It is also possible to select a number of specific scattering processes to investigate with the *target_process* parameter. This is especially useful when working with a single crystal process, that only scatters when the Bragg condition is met.\n",
    "\n",
    "Let us modify our existing loggers to view certain parts of the simulated system, and then rerun the simulation.  If mpi is installed, one can add mpi=N where N is the number of cores available to speed up the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "logger_zx.target_geometry = '\"outer_cryostat_wall,cryostat_wall\"'\n",
    "logger_zy.target_geometry = '\"sample_geometry\"'\n",
    "\n",
    "data = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.make_sub_plot(data, log=[True, False], orders_of_mag=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scattering order\n",
    "All loggers also have the option to only record given scattering orders. For example only record the second scattering.\n",
    "- order_total : Match given number of scattering events, counting all scattering events in the system\n",
    "- order_volume : Match given number of scattering events, only counting events in the current volume\n",
    "- order_volume_process : Match given number of scattering events, only counting events in current volume with current process\n",
    "\n",
    "We can modify our previous loggers to test out these features. The zx logger viewing from above will keep the target, but we remove the sample target on the zy logger, which is done by setting the *taget_geometry* to NULL. We choose to look at the second scattering event."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "logger_zx.order_total = 2\n",
    "\n",
    "logger_zy.target_geometry = '\"NULL\"'\n",
    "logger_zy.order_total = 2\n",
    "\n",
    "data = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.make_sub_plot(data, log=True, orders_of_mag=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration of additional logger components\n",
    "Here we add a few more loggers to showcase what kind of information that can be displayed.\n",
    "- 1D logger that logs scattered intensity as function of time\n",
    "- 2D abs_logger that logs absorption projected onto the scattering plane\n",
    "- 2DQ logger that logs scattering vector projected onto the scattering plane\n",
    "- 2D kf logger that logs final wavevector projected onto the scattering plane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logger_1D = instrument.add_component(\"logger_1D\", \"Union_logger_1D\", before=\"master\")\n",
    "logger_1D.variable = '\"time\"'\n",
    "logger_1D.min_value = 0.0006\n",
    "logger_1D.max_value = 0.0012\n",
    "logger_1D.n1 = 300\n",
    "logger_1D.filename = '\"logger_1D_time.dat\"'\n",
    "\n",
    "abs_logger_zx = instrument.add_component(\"abs_logger_space_zx\", \"Union_abs_logger_2D_space\",\n",
    "                                         before=\"master\")\n",
    "abs_logger_zx.set_AT([0,0,0], RELATIVE=sample_geometry)\n",
    "abs_logger_zx.D_direction_1 = '\"z\"'\n",
    "abs_logger_zx.D1_min = -0.12\n",
    "abs_logger_zx.D1_max = 0.12\n",
    "abs_logger_zx.n1 = 300\n",
    "abs_logger_zx.D_direction_2 = '\"x\"'\n",
    "abs_logger_zx.D2_min = -0.12\n",
    "abs_logger_zx.D2_max = 0.12\n",
    "abs_logger_zx.n2 = 300\n",
    "abs_logger_zx.filename = '\"abs_logger_zx.dat\"'\n",
    "\n",
    "logger_2DQ = instrument.add_component(\"logger_2DQ\", \"Union_logger_2DQ\", before=\"master\")\n",
    "logger_2DQ.Q_direction_1 = '\"z\"'\n",
    "logger_2DQ.Q1_min = -5.0\n",
    "logger_2DQ.Q1_max = 5.0\n",
    "logger_2DQ.n1 = 200\n",
    "logger_2DQ.Q_direction_2 = '\"x\"'\n",
    "logger_2DQ.Q2_min = -5.0\n",
    "logger_2DQ.Q2_max = 5.0\n",
    "logger_2DQ.n2 = 200\n",
    "logger_2DQ.filename = '\"logger_2DQ.dat\"'\n",
    "\n",
    "logger_2D_kf = instrument.add_component(\"logger_2D_kf\", \"Union_logger_2D_kf\",\n",
    "                                        before=\"master\")\n",
    "logger_2D_kf.Q_direction_1 = '\"z\"'\n",
    "logger_2D_kf.Q1_min = -2.5\n",
    "logger_2D_kf.Q1_max = 2.5\n",
    "logger_2D_kf.n1 = 200\n",
    "logger_2D_kf.Q_direction_2 = '\"x\"'\n",
    "logger_2D_kf.Q2_min = -2.5\n",
    "logger_2D_kf.Q2_max = 2.5\n",
    "logger_2D_kf.n2 = 200\n",
    "logger_2D_kf.filename = '\"logger_2D_kf.dat\"'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Runnig the simulation\n",
    "We now rerun the simulation with the new loggers. If mpi is installed, one can add mpi=N where N is the number of cores available to speed up the simulation."
   ]
  },
  {
   "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=3)\n",
    "ms.name_plot_options(\"logger_space_zy\", data, log=True, orders_of_mag=3)\n",
    "ms.name_plot_options(\"abs_logger_space_zx\", data, log=True, orders_of_mag=3)\n",
    "ms.name_plot_options(\"logger_1D\", data, log=True, orders_of_mag=3)\n",
    "ms.name_plot_options(\"logger_2DQ\", data, log=True, orders_of_mag=3)\n",
    "ms.name_plot_options(\"logger_2D_kf\", data, log=True, orders_of_mag=3)\n",
    "\n",
    "ms.make_sub_plot(data[0:2])\n",
    "ms.make_sub_plot(data[2:4])\n",
    "ms.make_sub_plot(data[4:6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting the data\n",
    "We see the scattered intensity as a function of time, here the peaks correspond to the direct beam intersecting the sides of the cryostat and sample. The source used release all neutrons at time 0, so it is a perfect pulse.\n",
    "\n",
    "The absorption monitor shows an image very similar to the scattered intensity, but this could be very different, for example when using materials meant as shielding.\n",
    "\n",
    "The 2D scattering vector is interesting, it shows a small sphere made of vertical lines, these are powder Bragg peaks. Since the wavevector is almost identical for all incoming neutrons, the first scattering can only access this smaller region of the space. The larger circle is incoherent scattering from second and later scattering events, where the incoming wavevector could be any direction since a scattering already happened.\n",
    "\n",
    "The 2D final wavevector plot shows mainly the powder Bragg peaks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animations and time series\n",
    "Several of the Union loggers sets up more than one McStas monitor, these include:\n",
    "- Union_logger_3D_space\n",
    "- Union_logger_2D_space_time\n",
    "- Union_logger_2D_kf_time\n",
    "\n",
    "The Union_logger_2D_space_time for example sets up a number of Union_logger_2D_space monitors that are limited to specific time intervals. This can be used to make an animation of the monitor, which we will demonstrate here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_2D_st = instrument.add_component(\"logger_2D_space_time\", \"Union_logger_2D_space_time\",\n",
    "                                     before=\"master\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_2D_st.time_bins = 36\n",
    "log_2D_st.time_min = 0.0007\n",
    "log_2D_st.time_max = 0.0011\n",
    "\n",
    "log_2D_st.set_AT([0,0,0], RELATIVE=sample_geometry)\n",
    "log_2D_st.D_direction_1 = '\"z\"'\n",
    "log_2D_st.D1_min = -0.12\n",
    "log_2D_st.D1_max = 0.12\n",
    "log_2D_st.n1 = 300\n",
    "log_2D_st.D_direction_2 = '\"x\"'\n",
    "log_2D_st.D2_min = -0.12\n",
    "log_2D_st.D2_max = 0.12\n",
    "log_2D_st.n2 = 300\n",
    "log_2D_st.filename = '\"logger_2D_space_time.dat\"'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "data = instrument.backengine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating an animation\n",
    "The plotter in McStasScript can create animations when supplied with many McStasData objects. We use name_search to find all the data from the relevant logger, and then a for loop to set plot options for each of them. Then the plotter can make an animation, which is saved as a gif."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ani_data = ms.name_search(\"logger_2D_space_time\", data)\n",
    "for frame in ani_data:\n",
    "    frame.set_plot_options(log=True, colormap=\"jet\")\n",
    "    \n",
    "ms.make_animation(ani_data, filename=\"animation_demo\", fps=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viewing the animation\n",
    "Some problem in the jupyter notebook prevents playing the gif directly, but it can be played from markdown. One has to refresh this cell when a new animation is written. It should be visible that the beam enters the cryostat from the left, scatters of the sample and illuminates the entire cryostat. Running this simulation with a larger ncount and more time_bins in the monitor will reveal more details in what happens. This is available below, but commented out as the simulation can take some time.\n",
    "\n",
    "![SegmentLocal](animation_demo.gif \"Animation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional: Increasing resolution\n",
    "Running this simulation with a larger ncount and more time_bins in the monitor will reveal more details in what happens. This is available below, but commented out as the simulation can take some time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "scroll-output"
    ]
   },
   "outputs": [],
   "source": [
    "log_2D_st.time_bins = 128\n",
    "instrument.settings(ncount=2E8, mpi=4)\n",
    "#data = instrument.backengine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ani_data = ms.name_search(\"logger_2D_space_time\", data)\n",
    "#for frame in ani_data:\n",
    "#    ms.set_plot_options(log=True, colormap=\"jet\", orders_of_mag=6)\n",
    "#    \n",
    "#ms.make_animation(ani_data, filename=\"animation_demo_long\", fps=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viewing the animation\n",
    "In the longer animation it is more evident that scattering from the aluminium hits the top and bottom of the outer cylinder.\n",
    "\n",
    "![SegmentLocal](animation_demo_long.gif \"Animation\")"
   ]
  }
 ],
 "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
}