QENS exercise

import mcstasscript as ms
import make_QENS_instrument
import quizlib

quiz = quizlib.QENS_Quiz()

QENS exercise#

This notebook contains code and questions for a McStas simulation of a simplified backscattering instrument that can investigate quasi-elastic scattering from samples.

Quasi-elastic scattering is inelastic scattering with small energy transfers and typically seen as a broadening of the elastic signal. At ESS the backscattering instrument under construction is called MIRACLES and uses the indirect geometry time of flight technique, here neutrons are scattered of the sample and some hit an analyzer crystal afterwards. This analyzer is angled such that if the neutron is scattered almost backwards, and due to Braggs law this will happen with a certain energy.

It turns out the precision of that energy is highest when the neutron is scattered back in the direction it came from, but most instruments choose a slightly lower angle to avoid hitting the sample a second time. The detector is then placed slightly above or below the sample.

Since the analyzer choose a specific energy, the final energy of the neutrons being recorded in the detector is known, this can be used to calculate the time between the neutron scattering in the sample and being detected, as both the distance and speed is known. The scattering time and the known pulse time can be used to calculate the time-of-flight, which with the known distance from source to sample gives the speed and thus energy before scattering in the sample. The difference between the known initial and final energy provide the energy transfer, which for backscattering can be down to \(\mu\)eV, where most other inelastic techniques look at meV.

In this notebook you will get a McStas model of simplified backscattering instrument which can be used to perform virtual experiments. You will run a few and answer questions about the results. You will also get to improve it and run experiments with an elastic calibration sample and a sample with an unknown QENS signal. We will use the Python McStas API McStasScript to work with the instrument, you can find documentation here.

Get the instrument object#

First we need the McStas instrument object. Here it is retrieved from a local python function that generates it.

instrument = make_QENS_instrument.make()

Investigate instrument#

First investigate the instrument object instrument using some of the available methods. All the methods that help do that start with the word show. In particular, look at what parameters are available and take a look at the instrument geometry.

Hide code cell content

 
instrument.show_parameters()

instrument.show_diagram()

double integration_time            = 21600      // [s] Time span of experiment
double energy_width_ueV            = 3          // Simulated energy range in 
                                                    micro eV 
double n_pulses                    = 3.0        // Number of pulses from source
double sample_distance             = 8.0        // [m] Source Sample distance
double sample_analyzer_distance    = 3.0        // [m] Sample analyzer distance
double analyzer_detector_distance  = 2.8        // [m] Sample analyzer distance
string sample_choice               = "Elastic"  // Choice of sample type, 
                                                    Elastic or QENS_sample
../_images/91eff03f8ebe45d718f836c3ec71b9acbcd58f405706d517e0e6b87da153d3de.png ../_images/93dd031ffa18722032949cde67d2610b86512fb506426762898c3a2325b895ce.png

Set parameters#

Before running the instrument we need to set some parameters on the instrument. The most important one is the sample_distance parameter describing the distance between the source and the sample. Given the need for high precision in determining the energy of the neutron, which of the following instrument lengths should be chosen?

  • A: 30 m

  • B: 60 m

  • C: 150 m

quiz.question_1()

Hide code cell content

 
quiz.question_1("C")

Correct!
Yes, a long instrument inherently provides a greater accuracy when determining
the neutron flight time and thus energy.

Set the sample_distance corresponding to the answer above and set the simulated energy width to 20 \(\mu\)eV. Keep the sample type to Elastic and number of pulses to 1. Use the set_parameters method on the instrument object.

Hide code cell content

 
instrument.set_parameters(sample_distance=150, energy_width_ueV=20, sample_choice='"Elastic"', n_pulses=1)

# Validate the instrument by giving it to the question_2
quiz.question_2(instrument)

Hide code cell output

 
Correct!
The parameters of the instrument were correctly set!

Instrument settings#

Before performing a simulation a few settings pertaining to the technical side should be set. These use a different method to clearly distinguish them from the instrument parameters. One important parameter is called output_path which sets the name of the generated data folder.

instrument.settings(ncount=5e6, mpi="auto", suppress_output=True, NeXus=True, output_path="first_run")

Run the simulation#

Now the simulation can be executed with the backengine method on the instrument object. Store the returned data in a python variable called data.

instrument.settings(suppress_output=True)

data = instrument.backengine()
data

[
 McStasData: signal_tof_1 type: 2D  I:2.66624 E:0.036036 N:26728.0,
 
 McStasData: signal_tof_3 type: 2D  I:2.99611 E:0.0400794 N:27038.0,
 
 McStasDataEvent: signal_tof_event_0 with 39606 events. Variables: p x y n id t,
 
 McStasDataEvent: signal_tof_event_1 with 41034 events. Variables: p x y n id t,
 
 McStasDataEvent: signal_tof_event_2 with 41434 events. Variables: p x y n id t,
 
 McStasDataEvent: signal_tof_event_3 with 41201 events. Variables: p x y n id t,
 
 McStasDataEvent: signal_tof_event_4 with 41480 events. Variables: p x y n id t,
 
 McStasDataEvent: signal_tof_event_5 with 41776 events. Variables: p x y n id t,
 
 McStasDataEvent: signal_tof_event_6 with 40980 events. Variables: p x y n id t,
 
 McStasDataEvent: signal_tof_event_7 with 39395 events. Variables: p x y n id t,
 
 McStasData: signal_space_1 type: 1D  I:2.66636 E:0.0360359 N:41034.0,
 
 McStasData: signal_space_3 type: 1D  I:2.99622 E:0.0400792 N:41201.0,
 
 McStasData: signal_time_1 type: 1D  I:2.66624 E:0.036036 N:26728.0,
 
 McStasData: signal_time_3 type: 1D  I:2.99611 E:0.0400794 N:27038.0]

Visualize the data#

The data objects in the returned list can be plotted with the McStasScript function make_sub_plot. Two types of data are returned, histogrammed data (McStasData) and event data (McStasDataEvent), only the histograms can be directly visualized.

ms.make_sub_plot(data, figsize=(9.5, 6))

Skipped plotting signal_tof_event_0 as it contains event data.
Skipped plotting signal_tof_event_1 as it contains event data.
Skipped plotting signal_tof_event_2 as it contains event data.
Skipped plotting signal_tof_event_3 as it contains event data.
Skipped plotting signal_tof_event_4 as it contains event data.
Skipped plotting signal_tof_event_5 as it contains event data.
Skipped plotting signal_tof_event_6 as it contains event data.
Skipped plotting signal_tof_event_7 as it contains event data.
../_images/be02d51a48786cf53753b435b39081d62b657066abcca2bcb6ddac512aa8f1af.png

Explanation of data#

The data is from two of the eight He3 tube detectors in the instrument, these can only be reached from the sample by being scattered almost backwards from a Si analyzer crystal. The detected signal from each of these two detectors are shown in three figures, two 1D and a 2D. The 1D monitors show the spatial and time distribution of the signal respectively, while the 2D monitor show the correlation between time and position.

Questions#

Look at the time distribution of the signal, which statement about this data is true?

  • A: The data looks like a typical inelastic signal

  • B: The data looks like the ESS pulse structure

  • C: The data looks like a typical elastic signal

  • D: The data looks like the analyzer selected to broad an energy range

quiz.question_3()

Hide code cell content

 
quiz.question_3("B")

Correct!
Yes, the time distribution is dominated by the ESS pulse structure

Is this a problem for a backscattering instrument?

  • A: Yes, the low time resolution means low energy resolution

  • B: No, the low time resolution is not necessary for high energy resolution

quiz.question_4()

Hide code cell content

 
quiz.question_4("A")

Correct!
Exactly, the time of flight is used to calculate the initial neutron energy

How can the instrument be improved?

  • A: Add a chopper to control the time aspect

  • B: Add a slit before sample to reduce the illuminated area

  • C: Add a slit before analyzer to ensure same angle

  • D: Add a spin polarizer to select spin state

quiz.question_5()

Hide code cell content

 
quiz.question_5("A")

Correct!
Yes, that would limit the uncertainty on the time of flight

Improve the instrument#

In order to improve the performance of the instrument, we will add a McStas component. The first aspect to consider when doing so is where to place it, both in the component sequence and its physical location. We start by investigating where to put the component in the component sequence, as that needs to be specified when adding the component, and the physical location can be updated later.

McStas sequence#

Use either the show_diagram or show_components method on the instrument object to get an overview of the component sequence in the instrument. Where would you place the new component?

  • A: After the source

  • B: Before the sample position

  • C: After the sample position

  • D: After the analyzer

Hide code cell content

 
instrument.show_diagram()
../_images/5e150711108476795d56b818379138aab76696f0fe57b026e2d82b9a7a1ee648.png ../_images/52029e5d053776c19bb8b8cb02d293a742ff785f167af88297026875424968ed.png 
quiz.question_6()

Hide code cell content

 
quiz.question_6("A")

Correct!
Yes, as we want to reduce the uncertainty in when the neutron was at the source,
placing the chopper as close to the source as possible is optimal.

Which component#

Now we need to select what type of component to add to the instrument, here we will need the DiskChopper component. Use the component_help method on the instrument object to learn more about this component.

Hide code cell content

 
instrument.component_help("DiskChopper")

 ___ Help DiskChopper _______________________________________________________________
|optional parameter|required parameter|default value|user specified value|
theta_0 = 0.0 [deg] // Angular width of the slits.
radius = 0.5 [m] // Radius of the disc
yheight [m] // Slit height (if = 0, equal to radius). Auto centering of beam at 
               half height. 
nu [Hz] // Frequency of the Chopper, omega=2*PI*nu (algebraic sign defines the 
           direction of rotation) 
nslit = 3.0 [1] // Number of slits, regularly arranged around the disk
jitter = 0.0 [s] // Jitter in the time phase
delay = 0.0 [s] // Time 'delay'
isfirst = 0.0 [0/1] // Set it to 1 for the first chopper position in a cw 
                       source (it then spreads the neutron time distribution) 
n_pulse = 1.0 [1] // Number of pulses (Only if isfirst)
abs_out = 1.0 [0/1] // Absorb neutrons hitting outside of chopper radius?
phase = 0.0 [deg] // Angular 'delay' (overrides delay)
xwidth = 0.0 [m] // Horizontal slit width opening at beam center
verbose = 0.0 [1] // Set to 1 to display Disk chopper configuration
-------------------------------------------------------------------------------------

Chopper calculations#

When adding a chopper one need to perform calculations to obtain the appropriate rotation frequency and phase. For this exercise, those calculations can be added to the instrument using a function included in the given Python file.

make_QENS_instrument.add_chopper_code(instrument)

To see what variables are used in the instrument, one can use the show_variables method like below. This will included the ones added by the above function.

Hide code cell content

 
instrument.show_variables()

DECLARE VARIABLES 
type    variable name                   array length  value               
------------------------------------------------------------------------
double  detector_offset                               0.25                
double  analyzer_angle                                                    
double  backscattering_wavelength                                         
double  backscattering_energy                                             
double  energy_width_meV                                                  
double  min_energy                                                        
double  max_energy                                                        
double  min_wavelength                                                    
double  max_wavelength                                                    
double  t_min_sample                                                      
double  t_max_sample                                                      
double  t_max_sample_pulses                                               
double  t_min                                                             
double  t_max                                                             
double  t_max_pulses                                                      
double  total_sample_detector_distance                                    
double  analyzer_d                                    3.135               
double  analyzer_q                                    2.0042058396107136  
double  hbar                                                              
double  sample_D                                                          
double  sample_tau                                                        
double  use_inelastic                                                     
double  Q_0                                                               
double  Q_1                                                               
double  Q_2                                                               
double  Q_3                                                               
double  Q_4                                                               
double  Q_5                                                               
double  Q_6                                                               
double  Q_7                                                               
double  chopper_distance                              6.5                 
double  chopper_delay                                                     
double  chopper_frequency                                                 
double  chopper_theta                                 12                  
double  chopper_time_width                                                

USER VARIABLES (per neutron, only use in EXTEND)
type    variable name                   array length  value               
------------------------------------------------------------------------
int     channel_index

Add chopper component and set parameters#

Use the add_component method on the instrument object to add a chopper. Place it in the component sequence according to your answer in question 6 by using either the before or after keyword argument.

The add_component method returns a component object, store that in a python variable.

Set the parameters:

  • yheight: 0.05 m

  • radius: 0.7 m

  • nslit: 1.0

  • nu, delay and theta_0: To the variables calculated in the instrument (use quotation marks)

Hide code cell content

 
chopper = instrument.add_component("chopper", "DiskChopper", after="source")
chopper.set_parameters(
    yheight=0.05,
    radius=0.7,
    nu="chopper_frequency",
    nslit=1.0,
    delay="chopper_delay",
    theta_0="chopper_theta",
)

To check if the component was added correctly, provide your instrument object to the question_7 function below.

# Validate the instrument again
quiz.question_7(instrument)

Hide code cell output

 
Correct!
The DiskChopper was added at the right point in the component sequence and with
correct parameters!

Placing the component in space#

The next task is to specify the physical location of the component, this is done using the set_AT method on the component object. This method takes a list of 3 numbers, corresponding to the x, y and z coordinates of the component.

One can also specify in what coordinate system one wants to work, which can be that of any preceding component. Use the RELATIVE keyword to work in the source coordinate system. The position of the chopper is needed for calculating phase, so it is available as a variable in the instrument, use this variable to set the position.

Hide code cell content

 
chopper.set_AT([0, 0, "chopper_distance"], RELATIVE="source")

To check if the component was placed correctly, provide your instrument object to the question_8 function below.

quiz.question_8(instrument)

Hide code cell output

 
Correct!
The DiskChopper was added at the right point in the space!

Verify new component#

Now that the chopper has been added to the instrument, lets show the component sequence again to verify it was added correctly.

Hide code cell content

 
instrument.show_diagram()
../_images/5e150711108476795d56b818379138aab76696f0fe57b026e2d82b9a7a1ee648.png ../_images/90c9baf03882f573a9712cee0a64b93b701c620816464d652b8a9d3978a301a4.png

Run improved instrument#

Run the improved instrument with the following parameters:

  • sample_distance: 150 m

  • energy_width_ueV: 20 ueV (This controls the simulated energy range, should be larger than inelastic signal width)

  • sample_choice: ‘“Elastic”’

  • frequency_multiplier: 10 (This controls the ratio between chopper and source frequency)

Use the settings method to set a reasonable name, including the sample type.

Store the resulting data in a variable called data_improved.

Hide code cell content

 
instrument.settings(output_path="QENS_elastic")
instrument.set_parameters(energy_width_ueV=20, sample_choice='"Elastic"', n_pulses=1, frequency_multiplier=10)

data_improved = instrument.backengine()

Plot the data using the make_sub_plot function.

Hide code cell content

 
ms.make_sub_plot(data_improved, figsize=(9.5, 6))

Skipped plotting signal_tof_event_0 as it contains event data.
Skipped plotting signal_tof_event_1 as it contains event data.
Skipped plotting signal_tof_event_2 as it contains event data.
Skipped plotting signal_tof_event_3 as it contains event data.
Skipped plotting signal_tof_event_4 as it contains event data.
Skipped plotting signal_tof_event_5 as it contains event data.
Skipped plotting signal_tof_event_6 as it contains event data.
Skipped plotting signal_tof_event_7 as it contains event data.
../_images/c6f444686519058edbb67f4d63e071ab0593a76751c018dcd052148f94914456.png

Time resolution#

Using the plotted data above its possible to get a estimate of the time resolution of the instrument.

Use the cursor to zoom and read of values for the time resolution of the instrument (FWHM). Insert the found time in the question below as a value in seconds.

quiz.question_9()

Hide code cell content

 
quiz.question_9(0.248167-0.247926)

Correct!
Yes!

Run with unknown quasi-elastic sample#

Next we will run a simulation with an interesting sample that actually has some inelastic signal. To cover the inelastic range, we need to simulate neutrons with a wider energy range at the source. Run with the following parameters:

  • sample_choice: ‘“QENS_sample”’

  • energy_width_ueV: 400 ueV

Again set a descriptive name, call the produced data data_unknown. If you have the time, it would be good to increase the ncount to 5e7 to improve the statistics.

Hide code cell content

 
instrument.settings(output_path="QENS_sample", ncount=5e7)
instrument.set_parameters(energy_width_ueV=400, sample_choice='"QENS_sample"',
                          frequency_multiplier=10)

data_unknown = instrument.backengine()

Plot the data with make_sub_plot.

Hide code cell content

 
ms.make_sub_plot(data_unknown, figsize=(9.5, 6), log=False, orders_of_mag=5)

Skipped plotting signal_tof_event_0 as it contains event data.
Skipped plotting signal_tof_event_1 as it contains event data.
Skipped plotting signal_tof_event_2 as it contains event data.
Skipped plotting signal_tof_event_3 as it contains event data.
Skipped plotting signal_tof_event_4 as it contains event data.
Skipped plotting signal_tof_event_5 as it contains event data.
Skipped plotting signal_tof_event_6 as it contains event data.
Skipped plotting signal_tof_event_7 as it contains event data.
../_images/75e24ad330930705ed20c34be55b02d504813909028fea157dce79cb90ab5c37.png

The two visualized channels show different broadening in the time direction. However, one of the channels (the one plotted as signal_time_1) suffers from a detector glitch that displaces part of the signal in time — you can see it in the 2D plot as a band of the tube shifted to later times. Ignore the glitch channel.

What is the time width (FWHM) of the signal in the other channel? Insert the answer with units of seconds.

quiz.question_10()

Hide code cell content

 
quiz.question_10(0.24823-0.24783)

Correct!
Yes!