SANS exercise

import mcstasscript as ms
import make_SANS_instrument
import quizlib

quiz = quizlib.SANS_Quiz()

SANS exercise#

In this notebook you will work with a McStas model of a simplified SANS instrument. You will have to answer questions in the notebook by working with this model, both by running simulations and expanding the model. We will use the Python McStas API McStasScript to work with the instrument, you can find documentation here.

SANS is an abbreviation for Small Angle Neutron Scattering, and as the name suggests is concerned with neutrons scattered at very small angles. Here we will look at a sample in solution composed of some simple geometry, which will cause an interesting scattering pattern on the detector.

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_SANS_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#

The first task is to investigate the instrument object instrument using some of the available methods available on that object. Each method that show something about the instrument starts with the word show, so you can use tab to autocomplete in the cell to see the relevant methods.

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   = 1    // [s] Time span of experiment
double wavelength         = 6.0  // [AA]  Mean wavelength of neutrons
double d_wavelength       = 3.0  // [AA]  Wavelength spread of neutrons
double n_pulses           = 1.0  // [1] Number of pulses from source
double sample_distance    = 8.0  // [m] Source Sample distance
double enable_sample      = 0    // [1] 0 for nothing, 1 for SANS sample
double detector_distance  = 2.0  // [m] Sample_detector_distance

Set parameters#

Before we run a simulation using the instrument, we need to set some parameters to the desired values. The most important one is the detector_distance parameter describing the distance between the sample and the detector. Given the need for high angular precision in determining the scattering angle of the neutron, which of these would be best?

  • A: 1 m

  • B: 2 m

  • C: 3 m

quiz.question_1()
Hide code cell content 
quiz.question_1("C")

Correct!
Yes, a large distance between sample and detector provides high angular
precision

Set the parameters of the instrument using the set_parameters method. In addition to the detector distance from the previous question, the parameters should be:

  • sample_distance: 150 m

  • wavelength: 6 Å

  • wavelength band: 1.5 Å

  • enable_sample: 0

  • n_pulses: 1

Hide code cell content 
instrument.set_parameters(
    sample_distance=150,
    wavelength=6,
    d_wavelength=1.5,
    enable_sample=0,
    detector_distance=3,
)

# Test your instrument by giving it to the question_2 function
quiz.question_2(instrument)
Hide code cell output 
Correct!
The parameters of the instrument were correctly set!

Instrument settings#

Before running the simulation, a few settings pertaining to computing options need to be specified. This is done with a different method to clearly distinguish these from the instrument parameters. One important setting is called output_path which sets the name of the generated folder with simulation output.

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

Run the instrument#

Now the simulation can be executed with the backengine method. This method returns a list of data object. Store this returned data in a python variable named data.

data = instrument.backengine()
data

[
 McStasData: signal type: 1D  I:45.0135 E:0.184418 N:3325220.0,
 
 McStasDataEvent: signal_tof_event with 3325219 events. Variables: p x y n id t,
 
 McStasData: signal_tof_all type: 2D  I:45.0135 E:0.184418 N:3325220.0,
 
 McStasData: signal_tof type: 2D  I:45.0135 E:0.184418 N:3325220.0]

Visualize the data#

The data objects in the returned list can be plotted with the McStasScript function make_plot. The plots can be customized, use these keyword arguments:

  • log : True

  • orders_of_mag : 5 (maximum orders of magnitudes plotted when using logarithmic plotting)

Plot overview#

The function should plot three graphs:

1D absorption logger#

This shows intensity along our single detector tube along its height y. Its coordinate system is such that 0 is 25 cm above the center of the beam.

1D TOF absorption logger#

Shows the same intensity as a function of the detector height, though here also as a function of when the neutron was detected. This one has the time axis restricted to show only one pulse.

1D TOF absorption logger#

Same as above, though will show all pulses if several are simulated.

ms.make_sub_plot(data, log=True, orders_of_mag=5)

Skipped plotting signal_tof_event as it contains event data.
../_images/a1358453b55b784afe7d8066e6f7251e17d238bc5124e6ba1e1a90aa12114af3.png

Interpretation of the data#

The detector is a He3 tube centered 25 cm above the beam height and with a metal casing.

What does the signal look like without sample?

  • A: Most of the signal close to the direct beam

  • B: Flat signal over detector height

  • C: Most of the signal is far away from the direct beam

quiz.question_3()
Hide code cell content 
quiz.question_3("A")

Correct!
Yes, since the detector is above the beam, -0.25 m would be beam height, and the
majority of the counts are seen at the lowest y values

Is this a problem for a SANS instrument?

  • A: Yes

  • B: No

quiz.question_4()
Hide code cell content 
quiz.question_4("A")

Correct!
Exactly, since this region correspond to the lowest measured angles, crucial for
small angle scattering

How can it be improved?

  • A: By adding a Velocity selector

  • B: By adding a Chopper

  • C: By adding a Beamstop

  • D: By adding a Slit

quiz.question_5()
Hide code cell content 
quiz.question_5("C")

Correct!
Yes, the issue arise from the direct beam scattering in the casing, a beamstop
can prevent this

Improve the instrument#

In order to improve the performance of the instrument, we will add a McStas component corresponding to the answer of question 5. The first aspect to consider when adding a component is where to place it, both in the component sequence and its physical location. The position in the component sequence needs to be specified when adding the component, so this will be the first decision.

McStas component 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: Before the detector position

instrument.show_diagram()
../_images/9b2dc003e98cfab50b9b1b3ac488d59a6ab46939bf99ec5f73bb9e075ed6e308.png ../_images/540993608c6f5376225f9aafadc312dba8dfc1b70b9e5342f476b653a80fb60f.png 
quiz.question_6()
Hide code cell content 
quiz.question_6("D")

Correct!
Yes, that would remove the part of the beam that did not scatter in the sample

Which component#

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

Hide code cell content 
instrument.component_help("Beamstop")

 ___ Help Beamstop __________________________________________________________________
|optional parameter|required parameter|default value|user specified value|
xmin = -0.05 [m] // Lower x bound
xmax = 0.05 [m] // Upper x bound
ymin = -0.05 [m] // Lower y bound
ymax = 0.05 [m] // Upper y bound
xwidth = 0.0 [m] // Width of beamstop (x). Overrides xmin, xmax.
yheight = 0.0 [m] // Height of beamstop (y). Overrides ymin, ymax.
radius = 0.0 [m] // radius of the beam stop in the z=0 plane, centered at Origo
-------------------------------------------------------------------------------------

Add beamstop component and set parameters#

Use the add_component method on the instrument to add a beamstop. Place it in the component sequence according to your answer in question 6 by using either the before or after keyword arguments in add_component. The add_component method returns the component object, save that in a Python variable.

The component should have these parameters set:

  • xwidth: 0.1 m

  • yheight: 0.02 m

After adding the component to the instrument object, provide the instrument object to the question 7 so it can be confirmed that the component was added correctly.

Hide code cell content 
beamstop = instrument.add_component("beamstop", "Beamstop", before="detector_position")
beamstop.set_parameters(xwidth=0.1, yheight=0.02)

# Validate instrument again
quiz.question_7(instrument)
Hide code cell output 
Correct!
The beamstop was added at the right point in the component sequence and with
right parameters!

Placing the component in space#

The next decision is the physical location of the component, this is done using the set_AT method on 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, this can be that of any preceding component. Use the RELATIVE keyword to work in the sample_position coordinate system. The instrument has a parameter called detector_distance, use this to place the beamstop 90% of the way from the sample to the detector.

After setting the position of the beamstop component, provide the instrument object to question 8 so it can be checked it was updated correctly.

Hide code cell content 
beamstop.set_AT([0, 0, "0.9*detector_distance"], RELATIVE="sample_position")

quiz.question_8(instrument)
Hide code cell output 
Correct!
The beamstop was added at the right point in the component sequence and space!

Verify new component#

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

instrument.show_diagram()
../_images/9b2dc003e98cfab50b9b1b3ac488d59a6ab46939bf99ec5f73bb9e075ed6e308.png ../_images/000d09db8279f4f0fd0073e498e5306addf646653dc39e2606ca70454ccf48e8.png

Run improved instrument#

Its time to run the improved instrument. We will use the same parameters as earlier, but now also set integration_time to \(5\times 10^4\), corresponding to 13.8 hours. This is supposed to emulate a long measurement of the instrument background only done rarely. We use the integration time in the continued workflow to estimate the error on the measured signal.

Store the returned data in a variable called background_data. Name the generated data folder “SANS_without_sample_1_pulse” using the output_path argument in the settings method of the instrument.

instrument.set_parameters(integration_time=5e4)
instrument.settings(output_path="SANS_without_sample_1_pulse")

background_data = instrument.backengine()
background_data

[
 McStasData: signal type: 1D  I:2.68061 E:0.0097524 N:2075290.0,
 
 McStasDataEvent: signal_tof_event with 2075286 events. Variables: p x y n id t,
 
 McStasData: signal_tof_all type: 2D  I:2.68061 E:0.0097524 N:2075290.0,
 
 McStasData: signal_tof type: 2D  I:2.68061 E:0.0097524 N:2075290.0]

Plot the resulting data.

ms.make_sub_plot(background_data, log=True, orders_of_mag=5)

Skipped plotting signal_tof_event as it contains event data.
../_images/b951a5c1617cd7eb72ffea7755e8df5de8e65dde9dcc8ab821bc6c7005830e54.png

Do you see an improvement compared to earlier results?

  • A: Yes

  • B: No

quiz.question_9()
Hide code cell content 
quiz.question_9("A")

Correct!
The large signal at low y has disappeared!

Run with sample#

Its time to add the sample, this can be done by updating the enable_sample parameter to 1 and calling the backengine method again. Here the integration time should be smaller, 500 s, as this run correspond to one of many measurements of different samples. Save the data in a variable called sample_data and use the settings method to provide a reasonable name to the run using the output_path argument.

instrument.settings(output_path="SANS_with_sample_1_pulse")
instrument.set_parameters(enable_sample=1, integration_time=500)

sample_data = instrument.backengine()

Plot the resulting data.

ms.make_sub_plot(sample_data, log=True, orders_of_mag=5)

Skipped plotting signal_tof_event as it contains event data.
../_images/57360c38f53ecd3b5018d06d5a8220c8898a72c54fdc7d15dce574dc01528d7b.png

Compare the results with and without sample. Where on the detector is the difference largest?

  • A: Lowest part of the detector

  • B: Middle of the detector

  • C: Top of the detector

quiz.question_10()
Hide code cell content 
quiz.question_10("A")

Correct!
Yes, there should be around 500 times larger count rate at the lowest part of
the detector

Increase the number of pulses#

Your final task is to re-run the simulations with and without sample, using 3 pulses from the source instead of 1 which was used so far. This is controlled using the n_pulses parameter on the instrument object. Set integration time to 5E4 s when running without sample and 500 s when running with sample. We will use this data throughout the exercises for the rest of the school.

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.

Hide code cell content 
# Without sample
instrument.settings(ncount=1.5e7, output_path="SANS_without_sample_3_pulse")
instrument.set_parameters(enable_sample=0, n_pulses=3, integration_time=5e4)
background_3_pulses = instrument.backengine()

# With sample
instrument.settings(ncount=1.5e7, output_path="SANS_with_sample_3_pulse")
instrument.set_parameters(enable_sample=1, n_pulses=3, integration_time=500)
sample_3_pulses = instrument.backengine()

Compare results with shorter runs#

By saving the results in different python variables, you can compare the two runs by plotting each. Use of the make_sub_plot function might make it easier to view each monitor side by side.

Hide code cell content 
def compare_monitors(monitor_name):
    sample_low_ncount = ms.name_search(monitor_name, sample_data)
    sample_high_ncount = ms.name_search(monitor_name, sample_3_pulses)
    ms.make_sub_plot([sample_low_ncount, sample_high_ncount], log=True, orders_of_mag=5)    

compare_monitors("signal")
../_images/f191b1db6d55f605719e69c308d0de4573335375f8b104d0e06b87bcddbced4e.png
Hide code cell content 
compare_monitors("signal_tof_all")
../_images/fe608bad778e05d1b488437f45867326cbc2cb92504ad8cbbbd3e5add2479772.png
Hide code cell content 
compare_monitors("signal_tof")
../_images/420a94304efbaeaa9596a1146146e4588bd2f73d6f90876efd605ab5eb3b2b9b.png
Contents