{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# McStasScript introduction\n", "This notebook shows how to use McStas and McStasScript to perform a basic simulation of a neutron diffractometer. The following software is required:\n", "- McStas (www.mcstas.org)\n", "- McStasScript (can be installed with python -m pip install McStasScript)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Anatomy of a McStas instrument\n", "\n", "In McStas a simulation is described using an instrument file. Such an instrument has five sections where code can be added to define the simulation to be perfomed.\n", "\n", "- Instrument definition\n", "- Declare section\n", "- Initialize section\n", "- Trace section\n", "- Finally section\n", "\n", "##### Instrument definition\n", "In the instrument definition it is possible to define *instrument parameters* which can be specified at run time and used in the remaining sections for either calculations or as direct input to the components.\n", "\n", "##### Declare section\n", "Here internal variables can be declared with C syntax.\n", "\n", "##### Initialize section\n", "The initialize section is used for performing calculations, typically using both instrument parameters and declared variables to calculate for example chopper phases, angles and similar. The calculations are performed using C syntax. These calculations are performed before the raytracing simulation, and thus only performed once in a given simulation.\n", "\n", "##### Trace section\n", "In the trace section McStas *components* are added, these are the building blocks of the simulation and correspond to different c codes that describe parts of neutron instruments or samples. Each component have a set of available parameters, some of which may be required. These will set the behavior of a component, a guide component may for example have parameters describing the physical shape and mirror reflectivity. Components also need to be placed in 3D space, and can be placed either in the absolute coordinate system or relative to a previously defined component.\n", "\n", "##### Finally section\n", "The finally section is very similar to the initialize section, here calculations can be performed after the raytracing has been completed, again using C syntax. This may be some brief data analysis or print of some status.\n", "\n", "### McStasScript python package and this tutorial\n", "The McStasScript python package provides an API to build and run such instruments files, but it is still necessary to have a basic understanding of the structure of the underlying instrument file and its capabilities and limitations. These tutorials will teach basic use of McStas through the McStasScript API without assuming expertise in the underlying McStas software." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Import the McStasScript package\n", "McStasScript needs to be imported into the users python environment, you can abbreviate the name for easier access." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mcstasscript as ms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### McStasScript configuration\n", "Before the first use of McStasScript it is necessary to configure the package so it can locate the McStas installation and call the binaries. One way to find the path is to open a terminal with the McStas environment and run:\n", "\n", "which mcrun\n", "\n", "This should return the path for the binary, and the mcstas path is usually just one step back." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "configurator = ms.Configurator()\n", "configurator.set_mcrun_path(\"/Applications/McStas-2.7.1.app/Contents/Resources/mcstas/2.7.1/bin/\")\n", "configurator.set_mcstas_path(\"/Applications/McStas-2.7.1.app/Contents/Resources/mcstas/2.7.1\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create an instrument object\n", "A McStas instrument is described with a McStas instrument object which is created using the *McStas_instr* method on the instr class. Creating an instrument object also reads available components, both in the work folder and from the McStas installation. By default, the work folder is the current work directory, but using the input_path keyword argument this can be change to avoid cluttering the folder containing notebooks.\n", "\n", "Here our instrument object for this tutorial is created, we give it the name python_tutorial." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument = ms.McStas_instr(\"python_tutorial\", input_path=\"run_folder\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Requesting help on source components\n", "The main building blocks used for creating a McStas simulation are the components. One can ask an instrument object which components are available, and get help for each component. Here we check what sources are available, and ask for help on the Source_div component." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.available_components()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.available_components(\"sources\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.component_help(\"Source_div\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding a component\n", "Now we are ready to add a component to our simulation which is done with the *add_component* method on our instrument. This method requires two inputs:\n", "\n", "- Nickname for the component used to refer to this component instance\n", "- Name of the component type to be used\n", "\n", "Here we want to make a component nicknamed \"source\" of type \"Source_div\".\n", "\n", "We also use the *print_components* method to confirm our component was added successfully. Running this code block multiple times result in an error, as McStas does not allow two components with the same nickname." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "src = instrument.add_component(\"source\", \"Source_div\")\n", "instrument.show_components()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Working with component objects\n", "The src object created by *add_component* can be used to modify the component. It also holds the information on the component, which can be shown by printing the object. This will tell us for example if any required parameters are yet to be set and the position of the component." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(src)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Modifying a component object\n", "The parameters of a component object can be modified as attributes. From the above print we know there are four required parameters, so we start by setting these and then print the resulting component status." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "src.xwidth = 0.1\n", "src.yheight = 0.05\n", "src.focus_aw = 1.2\n", "src.focus_ah = 2.3\n", "\n", "print(src)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Getting status of all parameters\n", "Printing a component only show the required parameters and user specified parameters, but it is also possible to see all parameters with the *show_parameters* method. This reminds us to set an energy or wavelength range for the source, as it is necessary to set one of these even though they are technically not required parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "src.show_parameters()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding an instrument parameter to control wavelength\n", "Controlling the wavelength range emitted by the source is best done with an instrument parameter, then this same parameter can be used to for example rotate a monochromator or set the range for an wavelength sensitive monitor. Adding an instrument parameter is done using the instrument method *add_parameter*, and it is possible to set a default value and comment. The method returns a parameter object that can be used to assign the parameter to a component. The current instrument parameters can be viewed with the *show_parameters* method on the isntrument object.\n", "\n", "The default type for instrument parameters is a double (floating point number), but other types can be selected if necessary by providing a type string before, here we also provide an example of an integer." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wavelength = instrument.add_parameter(\"wavelength\", value=5.0,\n", " comment=\"Wavelength in [Ang]\")\n", "\n", "order = instrument.add_parameter(\"int\", \"order\", value=1,\n", " comment=\"Monochromator order, integer\")\n", "instrument.show_parameters()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now our source component can have its parameters assigned to a instrument parameter, or even a mathematical expression using the variable. This allows us to set a reasonable wavelength range for our source component." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "src.lambda0 = wavelength\n", "src.dlambda = \"0.01*wavelength\" # When performing math use a string and the parameter name\n", "print(src)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using keyword arguments when adding a component\n", "When adding a component, several keyword arguments are available, for example for setting the position of the component.\n", "\n", "- AT set position with list of x,y,z coordinates\n", "- AT_RELATIVE set reference point for position (name of component instance or object)\n", "- ROTATED set rotation around x,y,z axis\n", "- ROTATED_RELATIVE set reference rotation (name of component instance or object)\n", "- RELATIVE set both reference position and rotation (name of component instance or object)\n", "\n", "We use this to set up a guide 2 meters after the source. The McStas coordinate system convention is such that the nominal beam direction is in the Z direction and with Y vertical against gravity. We use the component instance name as a string to refer to our source. The RELATIVE could also have been specified as src, which is our source object." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "guide = instrument.add_component(\"guide\", \"Guide_gravity\", AT=[0,0,2], RELATIVE=\"source\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we set the parameters for our guide component." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "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", "print(guide)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding calculations to an instrument file\n", "One of the advantages of McStas is the ease of adding calculations to the instrument. Here we calculate the rotation of a monochromator so that its scatters the wavelengths from our source. We need to declare variables using *add_declare_var* and append C code to initialize using *append_initialize*.\n", "\n", "For *add_declare_var* the first argument is the C type, usually double or int, the next is the variable name. A default value can be specified with the value keyword. Like when adding a parameter, a *add_declare* also returns an object that can be used to refer to this variable later.\n", "\n", "*append_initialize* just adds the given C code to the initialize section of the McStas instrument file. It is necessary to follow C syntax, for example remember semicolon at the end of statements." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "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);')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding the monochromator\n", "Here the monochromator is added, and we use the declared variables *mono_Q* and *mono_rotation* prepared above. Setting position and rotation can also be done using the *set_AT* and *set_ROTATED* methods on the component objects. Here it is also demonstrated how one can use either component objects or component names for the relative keyword.\n", "\n", "Rotation is specified around each axis, so rotation of our monochromator should be around the Y axis in order to keep the beam in the usual X-Z plane." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mono = instrument.add_component(\"mono\", \"Monochromator_flat\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "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\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(mono)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using an arm to define the beam direction\n", "As the beam changes direction at the monochromator, we wish to define the new direction to simplify adding latter components. This can be done with an Arm component, which performs no simulation but can be used as new coordinate reference. The outgoing direction correspond to one more rotation of *mono_rotation*." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "beam_direction = instrument.add_component(\"beam_dir\", \"Arm\", AT_RELATIVE=\"mono\")\n", "beam_direction.set_ROTATED([0, mono_rotation, 0], RELATIVE=\"mono\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding a sample\n", "We now add a powder sample using the PowderN component placed relative to our newly defiend beam direction. The chosen powder is Na2Ca3Al2F14 which is a standard sample due to its large number of available reflections." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sample = instrument.add_component(\"sample\", \"PowderN\",\n", " AT=[0, 0, 1.1], RELATIVE=beam_direction)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sample.radius = 0.015\n", "sample.yheight = 0.05\n", "sample.reflections = '\"Na2Ca3Al2F14.laz\"'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sample.print_long()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding a cylindrical monitor\n", "The flexible Monitor_nD component can be used to add a banana monitor (part of a cylinder). The component shape is specified using an option string. The restore_neutron parameter is set to 1 to allow other monitors to record each neutron.\n", "\n", "We have to specify a filename and option string here, and if we just use a string like \"banana.dat\" it would be interpreted as an instrument parameter called *banana.dat* and fail, so it is necessary to add single quotes around, '\"banana.dat\"'." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "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": [ "### Adding a psd monitor\n", "We also add a simple PSD (position sensitive detector) monitor to see the transmitted beam." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mon = instrument.add_component(\"monitor\", \"PSD_monitor\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mon.nx = 100\n", "mon.ny = 100\n", "mon.filename = '\"psd.dat\"'\n", "mon.xwidth = 0.05\n", "mon.yheight = 0.08\n", "mon.restore_neutron = 1\n", "mon.set_AT([0,0,0.1], RELATIVE=sample)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Print the components contained in an instrument\n", "Before performing the simulation, it is a good idea to check that the instrument contains the expected components and that they are appropriately placed in space. The *print_components* method is useful for this purpose." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.show_components()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Show instrument diagram\n", "While show components provides a text based overview of the instrument file, the instrument diagram can provide a visual overview. The diagram shows the sequence of components and color them based on their category. Arrows show the AT and ROTATED references within the instrument." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.show_diagram()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running the simulation\n", "Running the simulation is done in three steps\n", "\n", "- Setting the parameters with *set_parameters*\n", "- Setting the settings with *settings*\n", "- Running the McStas simulation with *backengine*\n", "\n", "The *set_parameters* method takes a value for each of the parameters defined in the instrument, here wavelength.\n", "\n", "Settings adjust settings for the simulations, a few examples can be seen here\n", "\n", "- ncount sets the number of rays\n", "- mpi sets the number of CPU cores used for execution (requires mpi installed)\n", "- output_path sets the name of the output folder\n", "- increment_folder_name if set to True, automatically changes the foldername if it already exists (default).\n", "\n", "The *backengine* method takes no parameters and just performs the simulation and returns the generated data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.set_parameters(wavelength=2.8)\n", "instrument.show_parameters()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.settings(ncount=5E6, output_path=\"data_folder/mcstas_basics\")\n", "instrument.show_settings()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "scroll-output" ] }, "outputs": [], "source": [ "data = instrument.backengine() # Perform simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting the data\n", "The *run_full_instrument* method returned a list of McStasData objects which can be plotted by the McStasScript plotter module. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ms.make_sub_plot(data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adjusting plots\n", "The McStasData objects contain preferences for how the data should be plotted, which can be modified using the functions module and the *name_plot_options* function. The function arguments are the name of the monitor component and a list of McStasData objects, then options are provided with the keyword arguments.\n", "\n", "The following plot options are often useful:\n", "- log [True or False] For plotting on logarithmic axis\n", "- orders_of_mag [number] When using logarithmic plotting, limits the maximum orders of magnitudes shown\n", "- left_lim [number] lower limit of plot x axis\n", "- right_lim [number] upper limit of plot x axis\n", "- bottom_lim [number] lower limit of plot y axis\n", "- top_lim [number] upper limit of plot y axis\n", "- colormap [string] name of matplotlib colormap to use" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ms.name_plot_options(\"monitor\", data, log=True)\n", "ms.name_plot_options(\"banana\", data, left_lim=90, right_lim=150)\n", "ms.make_sub_plot(data, fontsize=16)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Behind the scenes \n", "McStasScript writes the instrument file and uses mcrun to compile and run it. The file can be found in the input_path selected when the instrument object were created. We can print it here to see what was done behind the scenes." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "instrument.show_instrument_file()" ] }, { "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 }