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, nb=200, maxdiv=1.5, ndiv=200,
                          filename='"acceptance.dat"', restore_neutron=1)

instrument.set_parameters(wavelength=2.8)
instrument.settings(ncount=5E6, suppress_output=True, mpi="auto")
data = instrument.backengine()

ms.make_sub_plot(data)

../_images/cb3cd1d52ce74fc6f7e60ab7566c8211b9a01f83c7aea409da969268087164a1.png

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, nb=200, maxdiv=1.5, ndiv=200, restore_neutron=1)
    
    print(acceptance)
    print("")

COMPONENT acceptance_0 = DivPos_monitor(
  nb = 200, // [1]
  ndiv = 200, // [1]
  filename = "acceptance_0.dat", // [string]
  xwidth = 0.08, // [m]
  yheight = 0.05, // [m]
  maxdiv = 1.5, // [degrees]
  restore_neutron = 1 // [1]
) WHEN (n_reflections == 0)
AT (0, 0, 15.1) RELATIVE guide

COMPONENT acceptance_1 = DivPos_monitor(
  nb = 200, // [1]
  ndiv = 200, // [1]
  filename = "acceptance_1.dat", // [string]
  xwidth = 0.08, // [m]
  yheight = 0.05, // [m]
  maxdiv = 1.5, // [degrees]
  restore_neutron = 1 // [1]
) WHEN (n_reflections == 1)
AT (0, 0, 15.1) RELATIVE guide

COMPONENT acceptance_2 = DivPos_monitor(
  nb = 200, // [1]
  ndiv = 200, // [1]
  filename = "acceptance_2.dat", // [string]
  xwidth = 0.08, // [m]
  yheight = 0.05, // [m]
  maxdiv = 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()

../_images/5b9d74d8aad674f1ec5a51bd83da79dcd69b5ec090a2a0955c7e7caba11a9d70.png

../_images/fbef2c84aed79ec8f41251fe2baafaa450a7c62a2a9a63680c7d2f876e473393.png

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
  mpi:              auto 
  output_path:      data_folder/mcstas_EXTEND_WHEN
  run_path:         run_folder
  package_path:     /home/runner/micromamba/envs/mcstasscript-docs/share/mcstas/resources
  executable_path:  /home/runner/micromamba/envs/mcstasscript-docs/bin
  executable:       mcrun
  force_compile:    True
  NeXus:            False
  save_comp_pars:   False

ms.make_sub_plot(new_data)

../_images/3276da05fe1b1f2bbd075b3463da000a7a25fee1a38443c4faedd1c239a69d02.png

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: 18:39:04 on July 24, 2026
* Origin: ESS DMSC
* %INSTRUMENT_SITE: Generated_instruments
* 
* !!Please write a short instrument description (1 line) here!!
* 
* %Description
* Please write a longer instrument description here!
* 
* 
* %Parameters
* wavelength: [unit] Wavelength in [Ang]
* 
* %Link 
* 
* %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 = 10000000000000.0)
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.0,
 mright = 4.0, mtop = 0.0,
 mbottom = 0.0, G = -9.82)
AT (0, 0, 1.5) RELATIVE source
EXTEND %{
n_reflections = SCATTERED - 2;
%}

COMPONENT acceptance = DivPos_monitor(
 nb = 200, ndiv = 200,
 filename = "acceptance.dat", xwidth = 0.08,
 yheight = 0.05, maxdiv = 1.5,
 restore_neutron = 1)
AT (0, 0, 15.1) RELATIVE guide

COMPONENT acceptance_0 = DivPos_monitor(
 nb = 200, ndiv = 200,
 filename = "acceptance_0.dat", xwidth = 0.08,
 yheight = 0.05, maxdiv = 1.5,
 restore_neutron = 1)
WHEN (n_reflections == 0)
AT (0, 0, 15.1) RELATIVE guide

COMPONENT acceptance_1 = DivPos_monitor(
 nb = 200, ndiv = 200,
 filename = "acceptance_1.dat", xwidth = 0.08,
 yheight = 0.05, maxdiv = 1.5,
 restore_neutron = 1)
WHEN (n_reflections == 1)
AT (0, 0, 15.1) RELATIVE guide

COMPONENT acceptance_2 = DivPos_monitor(
 nb = 200, ndiv = 200,
 filename = "acceptance_2.dat", xwidth = 0.08,
 yheight = 0.05, maxdiv = 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

Advanced McStas features: EXTEND and WHEN

Contents