{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Advanced geometry using the Union components\n", "The Union components allow the user to construct advanced geometry from simple shapes. Each available shape has their own component, here are the currently available geometry components.\n", "- Union_box\n", "- Union_sphere\n", "- Union_cylinder\n", "- Union_cone\n", "\n", "They differ in their parameters describing the geometry, but are otherwise identical. In this notebook we will show how to construct hollow geometries with several layers, and that multiple scattering between these quickly result in complex behavior." ] }, { "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": "markdown", "metadata": {}, "source": [ "### Setting up some standard materials\n", "Before setting up the geometry, we need some material definition, here we set up aluminium and a sample powder." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Al_inc = instrument.add_component(\"Al_inc\", \"Incoherent_process\")\n", "Al_inc.sigma = 4*0.0082 # 4 atoms per unit cell\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*4*0.231/66.4 # barns [m^2 E-28]*Å^3 [m^3 E-30]=[m E-2]\n", "\n", "Sample_inc = instrument.add_component(\"Sample_inc\", \"Incoherent_process\")\n", "Sample_inc.sigma = 4*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*4*2.9464/1079.1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Set up source\n", "We will also need a source, and allow the wavelength to be tuned with a instrument parameter." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "src = instrument.add_component(\"source\", \"Source_div\")\n", "\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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Describing the geometry of a simple cryostat\n", "A cryostat is a complex geometry with several layers to consider. The way geometry is described in the Union components aims to make it easy to describe such systems. This is aciheved by allowing the simple geometries to overlap, and having a value called the priority to determine which is active in a given volume. If two geometries overlap, the overlapping region gets the physics from the geometry with the highest priority. In that way a cryostat model can be created by having a high priority for the sample in the center, and decreasing the priority as we move out.\n", "\n", "The ray tracing algorithm can however not handle if two geometries overlap perfectly, even with a single side. This could be two boxes sharing a side.\n", "\n", "Let us look at the parameters for a Union geometry component." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.component_help(\"Union_cylinder\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The geometry components have many parameters due to their flexibility, but only a few are needed for basic use.\n", "- material_string : string for selecting an available material\n", "- priority : number, in case of overlap the geometry with highest priority decides the material properties\n", "- geometrical parameters : Here radius and yheight\n", "\n", "In addition there is a focusing system where scattering of physical processes that support this can be forced to a certain direction, this is controlled with these parameters, but are rarely used:\n", "- target_index : relative component index of target\n", "- target_x : if target_index not set, relative x coordinate of target\n", "- target_y : if target_index not set, relative y coordinate of target\n", "- target_z : if target_index not set, relative z coordinate of target\n", "- focus_aw : angular width of focusing cone (either specify angular, box or circular)\n", "- focus_ah : angular height of focusing cone \n", "- focus_xw : spatial width of focusing cone (box type focusing)\n", "- focus_xh : spatial height of focusing\n", "- focus_r : spatial radius of focusing cone (circular)\n", "\n", "Finally there is p_interact, which is used for controlling Monte Carlo sampling frequency of the geometry, as it controls the probability for scattering occurring for any path before or after scattering.\n", "\n", "The remaining parameters including masks and number_of_activations are for advanced rules which will be described in a later tutorial." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### First geometry, a sample in a container\n", "\n", "We have defined the following materials that are available to us:\n", "- Al\n", "- Sample\n", "\n", "Lets start by building a simple powder container with a lid." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Set up loggers to check what is going on\n", "In order to view what geometry we have set up, we set up three loggers that view the scattering projected onto three different planes. These record the spatail distribution of scattering events." ] }, { "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.02\n", "logger_zx.D1_max = 0.02\n", "logger_zx.n1 = 300\n", "logger_zx.D_direction_2 = '\"x\"'\n", "logger_zx.D2_min = -0.02\n", "logger_zx.D2_max = 0.02\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.02\n", "logger_zy.D1_max = 0.02\n", "logger_zy.n1 = 300\n", "logger_zy.D_direction_2 = '\"y\"'\n", "logger_zy.D2_min = -0.02\n", "logger_zy.D2_max = 0.02\n", "logger_zy.n2 = 300\n", "logger_zy.filename = '\"logger_zy.dat\"'\n", "\n", "logger_xy = instrument.add_component(\"logger_space_xy\", \"Union_logger_2D_space\")\n", "logger_xy.set_RELATIVE(sample_geometry)\n", "logger_xy.D_direction_1 = '\"x\"'\n", "logger_xy.D1_min = -0.02\n", "logger_xy.D1_max = 0.02\n", "logger_xy.n1 = 300\n", "logger_xy.D_direction_2 = '\"y\"'\n", "logger_xy.D2_min = -0.02\n", "logger_xy.D2_max = 0.02\n", "logger_xy.n2 = 300\n", "logger_xy.filename = '\"logger_xy.dat\"'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add master component\n", "We need to remember to add a master component to actually perform the simulation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "master = instrument.add_component(\"master\", \"Union_master\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add banana monitor\n", "We are also interested in viewing some scattering data, here we add a banana monitor using the Monitor_nD component." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "banana = instrument.add_component(\"banana\", \"Monitor_nD\", RELATIVE=sample_geometry)\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\"'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run simulation\n", "Now we need to run the simulation to view the geometry we have built." ] }, { "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_geometry\")\n", "instrument.show_settings()\n", "\n", "data = instrument.backengine()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the data\n", "Due to the large differences between the scattered intensity from parts in the direct beam and outside, we use a logarithmic axis to display scattered intensity. We limit it to 4 orders of magnitude below the maximum intensity, otherwise a single very low intensity event can draw the intensity axis out to a large interval making it difficult to see the important nuances." ] }, { "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.name_plot_options(\"logger_space_xy\", data, log=True, orders_of_mag=4)\n", "ms.make_sub_plot(data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpretation of results\n", "The beam is narrower than the sample, but taller than the can, so some parts of the sample powder are not directly illuminated, and can thus be seen as a intensity area especially on the zx logger image. The aluminum scatters less, and so lower intensity still." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding a cryostat around the sample can\n", "We can add a crude model of a cryostat around our sample can by adding more Union geometry components. They have to be before the Union_master in the McStas instrument file, so we use the keyword argument *before* in the *add_component* method to specify this when adding the components.\n", "\n", "We also need to designate areas as empty, this is done using the default material Vacuum which has no absorption or scattering processes. In this way we can create several layers by decreasing the priority when going out." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inner_wall = instrument.add_component(\"cryostat_wall\", \"Union_cylinder\",\n", " before=\"master\")\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", " before=\"master\")\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", " before=\"master\")\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", " before=\"master\")\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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adjusting the logger view to see the larger cryostat area\n", "The loggers were only viewing a small area around the sample can, but this can be expanded as we still have access to the component objects." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "logger_zx.set_parameters(D1_min=-0.12, D1_max=0.12, D2_min=-0.12, D2_max=0.12)\n", "logger_zy.set_parameters(D1_min=-0.12, D1_max=0.12, D2_min=-0.12, D2_max=0.12)\n", "logger_xy.set_parameters(D1_min=-0.12, D1_max=0.12, D2_min=-0.12, D2_max=0.12)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Viewing the instrument\n", "The instrument can be viewed with the *show_instrument* method. The mock cryostat and detector can be seen in a 3D view." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.show_instrument()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run the updated instrument file\n", "Run the simulation with the added cryostat, since no parameters or settings are changed it is enough to just call the backengine function and grab the new data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "scroll-output" ] }, "outputs": [], "source": [ "data_cryo = instrument.backengine()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the data from the new simulation\n", "Here we increase the orders of magnitude of intensity plotted on the log plots. Try to play with these values to see how it changes the plots." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ms.make_sub_plot(data_cryo, log=[True, True, True, False], orders_of_mag=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpreting the data\n", "The different layers of the cryostat both result in scattering from the aluminium the beam has to move through, but also some increase intensity where it illuminated by scattering from the sample.\n", "\n", "### Comparing situation with and without cryostat\n", "It could be interesting to see what difference adding the cryostat did to the measured signal in the banana monitor, here we extract the numpy arrays and plot them manually with matplotlib for at direct comparison. Ensure you run the two simulations with the same wavelength in order for a comparison to be meaningful." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "banana_can = ms.name_search(\"banana\", data)\n", "banana_cryo = ms.name_search(\"banana\", data_cryo)\n", "\n", "import copy\n", "import numpy as np\n", "banana_diff = copy.deepcopy(banana_cryo)\n", "banana_diff.Intensity = banana_cryo.Intensity - banana_can.Intensity\n", "banana_diff.Error = np.sqrt(banana_cryo.Error**2 + banana_can.Error**2)\n", "\n", "import matplotlib.pyplot as plt\n", "plt.figure(figsize=(14,6))\n", "plt.plot(banana_can.xaxis, banana_can.Intensity, \"r\",\n", " banana_cryo.xaxis, banana_cryo.Intensity, \"b\",\n", " banana_diff.xaxis, banana_diff.Intensity-10.0, \"k\")\n", "plt.xlabel(\"2Theta [deg]\")\n", "plt.ylabel(\"Intensity [n/s]\")\n", "l = plt.legend([\"Sample in can\", \"Sample in can in cryostat\", \"Difference displaced to -10\"])" ] } ], "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 }