Advanced McStas features: EXTEND and WHEN
Contents
Advanced McStas features: EXTEND and WHEN¶
In this tutorial we will look at two advanced features in McStas, the EXTEND block and WHEN condition. Here we will use them to flag certain neutrons with EXTEND, and only record them in monitors when the flag is set using a WHEN condition.
import mcstasscript as ms
instrument = ms.McStas_instr("python_tutorial", input_path="run_folder",
output_path="data_folder/mcstas_EXTEND_WHEN")
Set up an example McStas instrument¶
First we set up an example instrument conisiting of a source, a guide and a position/divergence monitor. The guide is set up such that it only has mirrors on the left and right side, and absorbs neutrons if they hit the top or bottom. This is done to look at the horizontal behavior independently from the vertical, as this is easier to analyze.
src = instrument.add_component("source", "Source_simple")
src.focus_xw = guide_opening_w = 0.05
src.focus_yh = guide_opening_h = 0.06
src.dist = source_guide_distance = 1.5
wavelength = instrument.add_parameter("wavelength", value=5.0, comment="Wavelength in [Ang]")
src.set_parameters(xwidth=0.02, yheight=0.02, flux=1E13, lambda0=wavelength, dlambda="0.001*wavelength")
guide = instrument.add_component("guide", "Guide_gravity")
guide.set_AT([0, 0, source_guide_distance], RELATIVE=src)
guide.set_parameters(w1=guide_opening_w, w2=guide_opening_w,
h1=guide_opening_h, h2=guide_opening_h,
l=15, G=-9.82, mleft=4.0, mright=4.0, mtop=0.0, mbottom=0.0)
acceptance = instrument.add_component("acceptance", "DivPos_monitor")
acceptance.set_AT([0,0, guide.l + 0.1], RELATIVE=guide)
acceptance.set_parameters(xwidth=0.08, yheight=0.05, nh=200, maxdiv_h=1.5, ndiv=200,
filename='"acceptance.dat"', restore_neutron=1)
instrument.set_parameters(wavelength=2.8)
instrument.settings(ncount=5E6)
data = instrument.backengine()
INFO: Using directory: "/Users/madsbertelsen/PaNOSC/McStasScript/github/McStasScript/docs/source/tutorial/data_folder/mcstas_EXTEND_WHEN_1"
INFO: Regenerating c-file: python_tutorial.c
CFLAGS=
INFO: Recompiling: ./python_tutorial.out
mccode-r.c:2837:3: warning: expression result unused [-Wunused-value]
*t0;
^~~
1 warning generated.
INFO: ===
Detector: acceptance_I=1.58935e+07 acceptance_ERR=28258.1 acceptance_N=375738 "acceptance.dat"
INFO: Placing instr file copy python_tutorial.instr in dataset /Users/madsbertelsen/PaNOSC/McStasScript/github/McStasScript/docs/source/tutorial/data_folder/mcstas_EXTEND_WHEN_1
loading system configuration
ms.make_sub_plot(data)
Interpreting the data¶
Here we see an acceptance monitor, with position along the x-axis and divergence along the y-axis. The guide is under illuminated by the small source, so there are gaps in the acceptance diagram. We see the position and divergence of the beam consist of a large number of stripes, the ones with lowest divergence has the largest intensity.
Add an flag¶
A flag is just a name for a variable that records some information on the neutron during the simulation, and can be used later to make a decision. Here we could check how many times the ray was reflected in the guide.
We use an EXTEND block after a component to access variables internal to the component in the instrument scope. We declare a variable in the instrument scope called n_reflections. In the component scope, one can use the SCATTERED variable which contains the number of times the ray has encountered the SCATTER keyword within the component. Usually this is done when entering and leaving, and under each scattering / reflection, so the number of reflections is SCATTERED - 2.
instrument.add_declare_var("int", "n_reflections")
guide.append_EXTEND("n_reflections = SCATTERED - 2;")
print(guide)
COMPONENT guide = Guide_gravity(
w1 = 0.05, // [m]
h1 = 0.06, // [m]
w2 = 0.05, // [m]
h2 = 0.06, // [m]
l = 15, // [m]
mleft = 4.0, // [1]
mright = 4.0, // [1]
mtop = 0.0, // [1]
mbottom = 0.0, // [1]
G = -9.82 // [m/s2]
)
AT (0, 0, 1.5) RELATIVE source
EXTEND %{
n_reflections = SCATTERED - 2;
%}
Use the flag to limit what is recorded in a monitor¶
A WHEN statement can be used to activate / deactivate a component when some condition is true / false. For example we could require 0 reflection in our guide. We add a few monitors similar to the original, with the only difference being WHEN statements requiring 0, 1 or 2 reflections in the guide for the component to be active. We use a for loop to add the similar components, only changing the component instance name, filename and WHEN statement between each.
reflection_numbers = [0, 1, 2]
for reflections in reflection_numbers:
reflections_string = str(reflections)
component_name = "acceptance_" + reflections_string
acceptance = instrument.add_component(component_name, "DivPos_monitor")
acceptance.filename = '"acceptance_' + reflections_string + '.dat"'
acceptance.set_WHEN("n_reflections == " + reflections_string)
acceptance.set_AT([0,0, guide.l + 0.1], RELATIVE=guide)
acceptance.set_parameters(xwidth=0.08, yheight=0.05, nh=200, maxdiv_h=1.5, ndiv=200, restore_neutron=1)
print(acceptance)
print("")
COMPONENT acceptance_0 = DivPos_monitor(
nh = 200, // [1]
ndiv = 200, // [1]
filename = "acceptance_0.dat", // [string]
xwidth = 0.08, // [m]
yheight = 0.05, // [m]
maxdiv_h = 1.5, // [degrees]
restore_neutron = 1 // [1]
) WHEN (n_reflections == 0)
AT (0, 0, 15.1) RELATIVE guide
COMPONENT acceptance_1 = DivPos_monitor(
nh = 200, // [1]
ndiv = 200, // [1]
filename = "acceptance_1.dat", // [string]
xwidth = 0.08, // [m]
yheight = 0.05, // [m]
maxdiv_h = 1.5, // [degrees]
restore_neutron = 1 // [1]
) WHEN (n_reflections == 1)
AT (0, 0, 15.1) RELATIVE guide
COMPONENT acceptance_2 = DivPos_monitor(
nh = 200, // [1]
ndiv = 200, // [1]
filename = "acceptance_2.dat", // [string]
xwidth = 0.08, // [m]
yheight = 0.05, // [m]
maxdiv_h = 1.5, // [degrees]
restore_neutron = 1 // [1]
) WHEN (n_reflections == 2)
AT (0, 0, 15.1) RELATIVE guide
The instrument diagram shows the use of EXTEND and WHEN keywords by decorating the corresponding component boxes with a heavier outline or dashed outline respectively.
instrument.show_diagram()
Running the simulation¶
We now run the simulation with the new monitors to see how they differ from the original version.
instrument.set_parameters(wavelength=2.8)
instrument.settings(ncount=5E6)
instrument.show_settings()
new_data = instrument.backengine()
Instrument settings:
ncount: 5.00e+06
output_path: data_folder/mcstas_EXTEND_WHEN
run_path: run_folder
package_path: /Applications/McStas-2.7.1.app/Contents/Resources/mcstas/2.7.1
executable_path: /Applications/McStas-2.7.1.app/Contents/Resources/mcstas/2.7.1/bin/
executable: mcrun
force_compile: True
INFO: Using directory: "/Users/madsbertelsen/PaNOSC/McStasScript/github/McStasScript/docs/source/tutorial/data_folder/mcstas_EXTEND_WHEN_2"
INFO: Regenerating c-file: python_tutorial.c
CFLAGS=
INFO: Recompiling: ./python_tutorial.out
mccode-r.c:2837:3: warning: expression result unused [-Wunused-value]
*t0;
^~~
./python_tutorial.c:9457:21: warning: equality comparison with extraneous parentheses [-Wparentheses-equality]
if (( n_reflections == 0 ))
~~~~~~~~~~~~~~^~~~
./python_tutorial.c:9457:21: note: remove extraneous parentheses around the comparison to silence this warning
if (( n_reflections == 0 ))
~ ^ ~
./python_tutorial.c:9457:21: note: use '=' to turn this equality comparison into an assignment
if (( n_reflections == 0 ))
^~
=
./python_tutorial.c:9615:21: warning: equality comparison with extraneous parentheses [-Wparentheses-equality]
if (( n_reflections == 1 ))
~~~~~~~~~~~~~~^~~~
./python_tutorial.c:9615:21: note: remove extraneous parentheses around the comparison to silence this warning
if (( n_reflections == 1 ))
~ ^ ~
./python_tutorial.c:9615:21: note: use '=' to turn this equality comparison into an assignment
if (( n_reflections == 1 ))
^~
=
./python_tutorial.c:9773:21: warning: equality comparison with extraneous parentheses [-Wparentheses-equality]
if (( n_reflections == 2 ))
~~~~~~~~~~~~~~^~~~
./python_tutorial.c:9773:21: note: remove extraneous parentheses around the comparison to silence this warning
if (( n_reflections == 2 ))
~ ^ ~
./python_tutorial.c:9773:21: note: use '=' to turn this equality comparison into an assignment
if (( n_reflections == 2 ))
^~
=
4 warnings generated.
INFO: ===
Detector: acceptance_I=1.59141e+07 acceptance_ERR=28279.2 acceptance_N=375873 "acceptance.dat"
Detector: acceptance_0_I=2.05685e+06 acceptance_0_ERR=11084.3 acceptance_0_N=34434 "acceptance_0.dat"
Detector: acceptance_1_I=4.06628e+06 acceptance_1_ERR=15545.6 acceptance_1_N=68419 "acceptance_1.dat"
Detector: acceptance_2_I=3.85377e+06 acceptance_2_ERR=14739.7 acceptance_2_N=68445 "acceptance_2.dat"
INFO: Placing instr file copy python_tutorial.instr in dataset /Users/madsbertelsen/PaNOSC/McStasScript/github/McStasScript/docs/source/tutorial/data_folder/mcstas_EXTEND_WHEN_2
loading system configuration
ms.make_sub_plot(new_data)
Interpretation of the data¶
The original monitor is unchanged as it was not modified. On the monitors with different numbers of reflections, we see the middle line correspond to zero reflections, the two lines around those are for one reflection and so forth. This explains why the lines further from the center has lower intensity, as they underwent more reflections while also having a larger angle of incidence on the mirrors.
The McStas instrument file¶
We here show the generated McStas instrument file in order to clarify how this would be accomplished without the McStasScript API.
instrument.show_instrument_file()
/********************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright (C) 1997-2008, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* This file was written by McStasScript, which is a
* python based McStas instrument generator written by
* Mads Bertelsen in 2019 while employed at the
* European Spallation Source Data Management and
* Software Centre
*
* Instrument python_tutorial
*
* %Identification
* Written by: Python McStas Instrument Generator
* Date: 14:58:55 on April 25, 2023
* Origin: ESS DMSC
* %INSTRUMENT_SITE: Generated_instruments
*
*
* %Parameters
*
* %End
********************************************************************************/
DEFINE INSTRUMENT python_tutorial (
wavelength = 2.8 // Wavelength in [Ang]
)
DECLARE
%{
int n_reflections;
%}
INITIALIZE
%{
// Start of initialize for generated python_tutorial
%}
TRACE
COMPONENT source = Source_simple(
yheight = 0.02, xwidth = 0.02,
dist = 1.5, focus_xw = 0.05,
focus_yh = 0.06, lambda0 = wavelength,
dlambda = 0.001*wavelength, flux = 1E+13)
AT (0,0,0) ABSOLUTE
COMPONENT guide = Guide_gravity(
w1 = 0.05, h1 = 0.06,
w2 = 0.05, h2 = 0.06,
l = 15, mleft = 4,
mright = 4, mtop = 0,
mbottom = 0, G = -9.82)
AT (0,0,1.5) RELATIVE source
EXTEND %{
n_reflections = SCATTERED - 2;
%}
COMPONENT acceptance = DivPos_monitor(
nh = 200, ndiv = 200,
filename = "acceptance.dat", xwidth = 0.08,
yheight = 0.05, maxdiv_h = 1.5,
restore_neutron = 1)
AT (0,0,15.1) RELATIVE guide
COMPONENT acceptance_0 = DivPos_monitor(
nh = 200, ndiv = 200,
filename = "acceptance_0.dat", xwidth = 0.08,
yheight = 0.05, maxdiv_h = 1.5,
restore_neutron = 1)
WHEN (n_reflections == 0)
AT (0,0,15.1) RELATIVE guide
COMPONENT acceptance_1 = DivPos_monitor(
nh = 200, ndiv = 200,
filename = "acceptance_1.dat", xwidth = 0.08,
yheight = 0.05, maxdiv_h = 1.5,
restore_neutron = 1)
WHEN (n_reflections == 1)
AT (0,0,15.1) RELATIVE guide
COMPONENT acceptance_2 = DivPos_monitor(
nh = 200, ndiv = 200,
filename = "acceptance_2.dat", xwidth = 0.08,
yheight = 0.05, maxdiv_h = 1.5,
restore_neutron = 1)
WHEN (n_reflections == 2)
AT (0,0,15.1) RELATIVE guide
FINALLY
%{
// Start of finally for generated python_tutorial
%}
END