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 of the neutron to the sample position. Then the time at that moment 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 a small range of known and unknown samples. 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(input_path="run_folder")
The following components are found in the work directory / input_path:
     Union_abs_logger_1D_space_event.comp
These definitions will be used instead of the installed versions.

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()
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
string sample_choice      = "Elastic"  // Choice of sample type
double gamma_ueV          = 10         // [ueV] Energy width of known 
                                           quasi-elastic sample 
double analyzer_distance  = 3.0        // [m] Sample analyzer distance

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 5 \(\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=5, 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 folder.

instrument.settings(ncount=5e7, mpi=4, 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.

data = instrument.backengine()
data
[
 McStasData: signal_tof_all type: 2D  I:3.5467 E:0.0133054 N:369174.0,
 
 McStasData: signal_tof type: 2D  I:3.5467 E:0.0133054 N:369174.0,
 
 McStasDataEvent: signal_tof_event with 565220 events. Variables: p x y n id t,
 
 McStasData: signal_space type: 1D  I:3.54686 E:0.0133054 N:565222.0,
 
 McStasData: signal_time type: 1D  I:3.5467 E:0.0133054 N:369174.0]

Visualize the data#

The data objects in the returned list can be plotted with the McStasScript function make_sub_plot.

ms.make_sub_plot(data, figsize=(9.5, 8))
Skipped plotting signal_tof_event as it contains event data.
../_images/22e41cbfe8e2af67b395fd799747c4b61aafa27aeac802b8300dd0ea9664be53.png

Explanation of data#

The data is all from a He3 tube which can only be reached from the sample by being scattered almost backwards from a Si analyzer crystal. Two 1D monitors show the spatial and time distribution of the signal respectively. The two 2D monitors show the correlation between time and position, though one will zoom out to capture multiple pulses if n_pulse is set to more than one.

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/ae9c600803bb4efe54261156031f6d2f1b71abd6fd5135affc7022da81d1ae52.png ../_images/f54807bc203d1aafe35ad3796a67f840d51b70a448df5d7b4a9d0430afe70f6b.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  enable_elastic                                               
double  enable_known_quasi                                           
double  enable_unknown_quasi                                         
double  gamma_meV                                                    
double  analyzer_direction                       30                  
double  analyzer_d                               3.135               
double  analyzer_q                               2.0042058396107136  
double  t_min                                                        
double  t_max                                                        
double  t_max_pulses                                                 
double  chopper_distance                         6.5                 
double  chopper_delay                                                
double  chopper_frequency                                            
double  chopper_theta                            12                  
double  chopper_time_width                                           

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 added correctly, provide your instrument object to the question_7 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/ae9c600803bb4efe54261156031f6d2f1b71abd6fd5135affc7022da81d1ae52.png ../_images/3db5b51b75fa152c05cc4f5533dbeb7127c0d12915e89aaac6fed6b6db274a02.png

Run improved instrument#

Run the improved instrument with the following parameters:

  • sample_distance: 150 m

  • energy_width_ueV: 5 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 and number of pulses currently used.

Store the resulting data in a variable called data_improved.

Hide code cell content
instrument.settings(output_path="QENS_elastic_1_pulse")
instrument.set_parameters(energy_width_ueV=5, 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, 8))
Skipped plotting signal_tof_event as it contains event data.
../_images/d76083f56f34661c4b72c4a45d7c1814d34b6f8eb411b9948da9c5190bd56423.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.248503-0.248258)
Correct!
Yes!

Run with known calibration sample#

Next we will run a simulation with a known calibration sample, its energy width can be adjusted with the gamma_ueV (HWHM). Run with the following parameters:

  • sample_choice: ‘“Known_quasi-elastic”’

  • gamma_ueV: 12 ueV

  • energy_width_ueV: 150 ueV

Again set a descriptive name, call the produced data data_known.

Hide code cell content
instrument.settings(output_path="QENS_known_quasi_elastic_1_pulse")
instrument.set_parameters(energy_width_ueV=150, sample_choice='"Known_quasi-elastic"',
                          gamma_ueV=12, frequency_multiplier=10)

data_known = instrument.backengine()

Plot the data with make_sub_plot.

Hide code cell content
ms.make_sub_plot(data_known, figsize=(9.5, 8))
Skipped plotting signal_tof_event as it contains event data.
../_images/7a8a28d01ef8c7690a2591bada86f9bd10c9cb7b1704c55f5d8fff6a6ddefcf4.png

What is the time width when using a known sample with 12 ueV broadening? Insert the answer with units of seconds.

quiz.question_10()
Hide code cell content
quiz.question_10(0.249004-0.247830)
Correct!
Yes!

Run with unknown sample#

Finally we can run a virtual experiment with an unknown sample.

  • sample_choice: "Unknown_quasi-elastic"

  • energy_width_ueV: 150 ueV

Again set a descriptive name, call the produced data data_unknown.

Hide code cell content
instrument.settings(
    output_path="QENS_unknown_quasi_elastic_1_pulse",
)

instrument.set_parameters(
    energy_width_ueV=150,
    sample_choice='"Unknown_quasi-elastic"',
    n_pulses=1,
    frequency_multiplier=10,
)

data_unknown = instrument.backengine()

Plot the data using make_sub_plot.

Hide code cell content
ms.make_sub_plot(data_unknown, figsize=(9.5, 8))
Skipped plotting signal_tof_event as it contains event data.
../_images/d5d20e04b67c69dc998f4bd20ac1aa96b67c8d76b17fa3e9a47bfc04d75e9686.png

Estimate unknown sample energy width#

By comparing the data from the known and unknown sample, provide an estimate of the energy width of the unknown sample. Insert the result below in ueV (HWHM). This questions is optional, if short on time go to the next question.

quiz.question_11()
Hide code cell content
quiz.question_11((0.248782 - 0.247976)/(0.249004-0.247830)*12.0)
Correct!
Yes! Lets see how this early estimate holds up with more thorough analysis over
the next days.

Increase the number of pulses#

Your final task is to re-run the simulations for the 3 different samples, using 3 pulses instead of 1. We will use this data in the exercises for the rest of the week.

Hints:

  • Change the destination folder so that you don’t overwrite the results from the 1-pulse simulations.

  • Remember to adjust the ncount accordingly, we would like 3 times more rays now that we use 3 pulses.

Solution:

Hide code cell content
# Elastic sample
instrument.settings(
    ncount=1.5e8,
    output_path="QENS_elastic_3_pulse",
)
instrument.set_parameters(
    energy_width_ueV=5,
    sample_choice='"Elastic"',
    n_pulses=3,
    frequency_multiplier=10,
)
elastic_3_pulses = instrument.backengine()


# Known calibration sample
instrument.settings(
    ncount=1.5e8,
    output_path="QENS_known_quasi_elastic_3_pulse",
)
instrument.set_parameters(
    energy_width_ueV=150,
    sample_choice='"Known_quasi-elastic"',
    gamma_ueV=12,
    n_pulses=3,
    frequency_multiplier=10,
)
known_3_pulses = instrument.backengine()


# Unknown sample
instrument.settings(
    ncount=1.5e8,
    output_path="QENS_unknown_quasi_elastic_3_pulse",
)
instrument.set_parameters(
    energy_width_ueV=150,
    sample_choice='"Unknown_quasi-elastic"',
    n_pulses=3,
    frequency_multiplier=10,
)
unknown_3_pulses = instrument.backengine()