<jupyter_start><jupyter_text>McStasScript introductionThis notebook shows how to use McStas and McStasScript to perform a basic simulation of a neutron diffractometer. The following software is required:- McStas (www.mcstas.org)- McStasScript (can be installed with python -m pip install McStasScript)  Anatomy of a McStas instrumentIn 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.- Instrument definition- Declare section- Initialize section- Trace section- Finally section Instrument definitionIn 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. Declare sectionHere internal variables can be declared with C syntax. Initialize sectionThe 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. Trace sectionIn 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. Finally sectionThe 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. McStasScript python package and this tutorialThe 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.  Import the McStasScript packageThe McStasScript modules intended for normal use is located in the interface submodule, and one usually imports the necessary modules from there.<jupyter_code>from mcstasscript.interface import instr, functions, plotter<jupyter_output><empty_output><jupyter_text>McStasScript configurationBefore 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:which mcrunThis should return the path for the binary, and the mcstas path is usually just one step back.<jupyter_code>configurator = functions.Configurator()
configurator.set_mcrun_path("/Applications/McStas-2.5.app/Contents/Resources/mcstas/2.5/bin/")
configurator.set_mcstas_path("/Applications/McStas-2.5.app/Contents/Resources/mcstas/2.5")<jupyter_output><empty_output><jupyter_text>Create an instrument objectA 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.Here our instrument object for this tutorial is created, we give it the name python_tutorial.<jupyter_code>instrument = instr.McStas_instr("python_tutorial", input_path="run_folder")<jupyter_output><empty_output><jupyter_text>Requesting help on source componentsThe 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.<jupyter_code>instrument.show_components()
instrument.show_components("sources")
instrument.component_help("Source_div")<jupyter_output><empty_output><jupyter_text>Adding a componentNow 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:- Nickname for the component used to refer to this component instance- Name of the component type to be usedHere we want to make a component nicknamed "source" of type "Source_div".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.<jupyter_code>src = instrument.add_component("source", "Source_div")
instrument.print_components()<jupyter_output><empty_output><jupyter_text>Working with component objectsThe 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 with the *print_long* method. This will tell us for example if any required parameters are yet to be set and the position of the component.<jupyter_code>src.print_long()<jupyter_output><empty_output><jupyter_text>Modifying a component objectThe 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.<jupyter_code>src.xwidth = 0.1
src.yheight = 0.05
src.focus_aw = 1.2
src.focus_ah = 2.3

src.print_long()<jupyter_output><empty_output><jupyter_text>Getting status of all parametersUsing *print_long* on 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.<jupyter_code>src.show_parameters()<jupyter_output><empty_output><jupyter_text>Adding an instrument parameter to control wavelengthControlling 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 current instrument parameters can be viewed with the *show_parameters* method on the isntrument object.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.<jupyter_code>instrument.add_parameter("wavelength", value=5.0, comment="Wavelength in [Ang]")
instrument.add_parameter("int", "order", value=1, comment="Monochromator order, integer")
instrument.show_parameters()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>src.lambda0="wavelength"
src.dlambda="0.01*wavelength"
src.print_long()<jupyter_output><empty_output><jupyter_text>Using keyword arguments when adding a componentWhen adding a component, several keyword arguments are available, for example for setting the position of the component.- AT set position with list of x,y,z coordinates- AT_RELATIVE set reference point for position (name of component instance or object)- ROTATED set rotation around x,y,z axis- ROTATED_RELATIVE set reference rotation (name of component instance or object)- RELATIVE set both reference position and rotation (name of component instance or object)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.<jupyter_code>guide = instrument.add_component("guide", "Guide_gravity", AT=[0,0,2], RELATIVE="source")<jupyter_output><empty_output><jupyter_text>Next we set the parameters for our guide component.<jupyter_code>guide.w1 = 0.05
guide.w2 = 0.05
guide.h1 = 0.05
guide.h2 = 0.05
guide.l = 8.0
guide.m = 3.5
guide.G = -9.82

guide.print_long()<jupyter_output><empty_output><jupyter_text>Adding calculations to an instrument fileOne 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*.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.*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.<jupyter_code>instrument.add_declare_var("double", "mono_Q", value=1.714) # Q for Ge 311
instrument.add_declare_var("double", "wavevector")
instrument.append_initialize("wavevector = 2.0*PI/wavelength;")

instrument.add_declare_var("double", "mono_rotation")
instrument.append_initialize("mono_rotation = asin(mono_Q/(2.0*wavevector))*RAD2DEG;")
instrument.append_initialize('printf("monochromator rotation = %g deg\\n", mono_rotation);')<jupyter_output><empty_output><jupyter_text>Adding the monochromatorHere 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.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.<jupyter_code>mono = instrument.add_component("mono", "Monochromator_flat")
mono.zwidth = 0.05
mono.yheight = 0.08
mono.Q = "mono_Q"
mono.set_AT([0, 0, 8.5], RELATIVE=guide)
mono.set_ROTATED([0, "mono_rotation", 0], RELATIVE="guide")
mono.print_long()<jupyter_output><empty_output><jupyter_text>Using an arm to define the beam directionAs 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*.<jupyter_code>beam_direction = instrument.add_component("beam_dir", "Arm", AT_RELATIVE="mono")
beam_direction.set_ROTATED([0, "mono_rotation", 0], RELATIVE="mono")<jupyter_output><empty_output><jupyter_text>Adding a sampleWe 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.<jupyter_code>sample = instrument.add_component("sample", "PowderN", AT=[0,0,1.1], RELATIVE="beam_dir")
sample.radius = 0.015
sample.yheight = 0.05
sample.reflections = '"Na2Ca3Al2F14.laz"'
sample.print_long()<jupyter_output><empty_output><jupyter_text>Adding a cylindrical monitorThe 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.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"'.<jupyter_code>banana = instrument.add_component("banana", "Monitor_nD", RELATIVE=sample)
banana.xwidth = 2.0
banana.yheight = 0.3
banana.restore_neutron = 1
banana.filename = '"banana.dat"'
banana.options = '"theta limits=[5 175] bins=150, banana"'<jupyter_output><empty_output><jupyter_text>Adding a psd monitorWe also add a simple PSD (position sensitive detector) monitor to see the transmitted beam.<jupyter_code>mon = instrument.add_component("monitor", "PSD_monitor")
mon.nx = 100
mon.ny = 100
mon.filename = '"psd.dat"'
mon.xwidth = 0.05
mon.yheight = 0.08
mon.restore_neutron = 1
mon.set_AT([0,0,0.1], RELATIVE=sample)<jupyter_output><empty_output><jupyter_text>Print the components contained in an instrumentBefore 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.<jupyter_code>instrument.print_components()<jupyter_output><empty_output><jupyter_text>Running the simulationThe instrument object has a method called *run_full_instrument* to execute the simulation and return the data. A number of keyword arguments are available to control the execution of the simulation.- ncount sets the number of rays- mpi sets the number of CPU cores used for execution (requires mpi installed)- foldername sets the name of the output folder- increment_folder_name if set to True, automatically changes the foldername if it already exists.- parameters allows setting instrument parameters using a python dictionary<jupyter_code>data = instrument.run_full_instrument(ncount=5E6, foldername="data_folder/mcstas_basics",
                                      increment_folder_name=True,
                                      parameters={"wavelength" : 2.8})<jupyter_output><empty_output><jupyter_text>Plotting the dataThe *run_full_instrument* method returned a list of McStasData objects which can be plotted by the McStasScript plotter module.<jupyter_code>plotter.make_sub_plot(data)<jupyter_output><empty_output><jupyter_text>Adjusting plotsThe 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.The following plot options are often useful:- log [True or False] For plotting on logarithmic axis- orders_of_mag [number] When using logarithmic plotting, limits the maximum orders of magnitudes shown- left_lim [number] lower limit of plot x axis- right_lim [number] upper limit of plot x axis- bottom_lim [number] lower limit of plot y axis- top_lim [number] upper limit of plot y axis- colormap [string] name of matplotlib colormap to use<jupyter_code>functions.name_plot_options("monitor", data, log=True)
functions.name_plot_options("banana", data, left_lim=90, right_lim=150)
plotter.make_sub_plot(data, fontsize=16)<jupyter_output><empty_output><jupyter_text>Behind the scenes 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.<jupyter_code>with open("run_folder/python_tutorial.instr") as file:
    data = file.read()
    print(data)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.058454855724480295, 'fraction_numerical': 0.008687558175612783, 'mean_word_length': 5.91970802919708, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 93, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '490812', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4307, 'n_tokens_neox': 4121, 'n_words': 1948, 'path': 'tutorial/McStasScript_tutorial_1_the_basics.ipynb', 'repo_name': 'PaNOSC-ViNYL/McStasScript', 'size': '16115'}