McStasRead demo#

[1]:

%matplotlib widget

[2]:

from read_example import make_instrument, plot

[3]:

instr = make_instrument()
instr.show_parameters()

 wavelength        = 2.0  //
 delta_wavelength  = 0.1  //

[4]:

instr.set_parameters(wavelength=1.8, delta_wavelength=1.3)
instr.settings(ncount=1E6)

[5]:

instr.show_diagram()

Run simulation#

[6]:

data = instr.backengine()


INFO: Using directory: "/home/runner/work/McStasToX/McStasToX/docs/user-guide/test"
INFO: Regenerating c-file: test.c
CFLAGS=

WARNING:
 The parameter format of sample is initialized
 using a static {,,,} vector.
  -> Such static vectors support literal numbers ONLY.
  -> Any vector use of variables or defines must happen via a
     DECLARE/INITIALIZE pointer.


-----------------------------------------------------------

Generating single GPU kernel or single CPU section layout:

-----------------------------------------------------------

Generating GPU/CPU -DFUNNEL layout:

-----------------------------------------------------------
INFO: Recompiling: ./test.out
lto-wrapper: warning: using serial compilation of 3 LTRANS jobs
lto-wrapper: note: see the '-flto' option documentation for more information
/home/runner/work/McStasToX/McStasToX/.pixi/envs/docs/bin/../libexec/gcc/x86_64-conda-linux-gnu/14.3.0/ld: /tmp/cc52STJs.ltrans1.ltrans.o: in function `read_line_data.part.0.constprop.0':
/home/runner/work/McStasToX/McStasToX/docs/user-guide/./test.c:10224:(.text.read_line_data.part.0.constprop.0+0xea): warning: the use of `tmpnam' is dangerous, better use `mkstemp'
INFO: ===
Opening input file '/home/runner/work/McStasToX/McStasToX/.pixi/envs/docs/share/mcstas/resources/data/Cu.laz' (Table_Read_Offset)
Table from file 'Cu.laz' (block 1) is 19 x 18 (x=1:6), constant step. interpolation: linear
  '# TITLE *-Cu-[FM3-M] Otte, H.M.[1961];# CELL 3.615050 3.615050 3.615050 90. ...'
PowderN: sample: Reading 19 rows from Cu.laz
PowderN: sample: Read 19 reflections from file 'Cu.laz'
PowderN: sample: Vc=47.24 [Angs] sigma_abs=15.12 [barn] sigma_inc=2.2 [barn] reflections=Cu.laz
*** TRACE end ***
Events:   "direct_event_square_signal_dat_list.p.x.y.n.id.t.L"
Events:   "scattered_event_square_signal_dat_list.p.x.y.n.id.t"
Events:   "direct_event_banana_signal_dat_list.p.th.y.n.id.t"
Events:   "scattered_event_banana_signal_dat_list.p.th.y.n.id.t"
PowderN: sample: Info: you may highly improve the computation efficiency by using
    SPLIT 8 COMPONENT sample=PowderN(...)
  in the instrument description test.instr.
INFO: Placing instr file copy test.instr in dataset /home/runner/work/McStasToX/McStasToX/docs/user-guide/test
INFO: Placing generated c-code copy test.c in dataset /home/runner/work/McStasToX/McStasToX/docs/user-guide/test


[7]:

data

[7]:

[
 McStasDataEvent: Banana_1 with 5543 events. Variables: p th y n id t,

 McStasDataEvent: Square_1 with 4175 events. Variables: p x y n id t L,

 McStasDataEvent: Banana_2 with 18799 events. Variables: p th y n id t,

 McStasDataEvent: Square_2 with 7932 events. Variables: p x y n id t]

Histogram data to display#

[8]:

hist_data = [mon.make_2d(mon.variables[1], "y") for mon in data]

[9]:

import mcstasscript as ms
ms.make_sub_plot(hist_data)

[10]:

file_path = data[0].original_data_location
print(file_path)

/home/runner/work/McStasToX/McStasToX/docs/user-guide/test

Show list of components found in NeXus file#

[11]:

import mcstastox

with mcstastox.Read(file_path) as loaded_data:
    loaded_data.show_components()

All components in file:
2026-04-30 13:35:42,917 INFO:All components in file:
source
2026-04-30 13:35:42,918 INFO:source
sample_position
2026-04-30 13:35:42,919 INFO:sample_position
sample
2026-04-30 13:35:42,919 INFO:sample
detector_direction_square_1
2026-04-30 13:35:42,920 INFO:detector_direction_square_1
Square_1
2026-04-30 13:35:42,921 INFO:Square_1
detector_direction_square_2
2026-04-30 13:35:42,921 INFO:detector_direction_square_2
Square_2
2026-04-30 13:35:42,922 INFO:Square_2
detector_direction_banana_1
2026-04-30 13:35:42,922 INFO:detector_direction_banana_1
Banana_1
2026-04-30 13:35:42,923 INFO:Banana_1
detector_direction_banana_2
2026-04-30 13:35:42,923 INFO:detector_direction_banana_2
Banana_2
2026-04-30 13:35:42,924 INFO:Banana_2

Show list of these with available data#

[12]:

with mcstastox.Read(file_path) as loaded_data:
    loaded_data.show_components_with_data()

All components with data in file:
All components with data in file:
2026-04-30 13:35:42,933 INFO:All components with data in file:
Square_1
Square_1
2026-04-30 13:35:42,934 INFO:Square_1
Square_2
Square_2
2026-04-30 13:35:42,935 INFO:Square_2
Banana_1
Banana_1
2026-04-30 13:35:42,936 INFO:Banana_1
Banana_2
Banana_2
2026-04-30 13:35:42,937 INFO:Banana_2

[13]:

with mcstastox.Read(file_path) as loaded_data:
    variables = loaded_data.get_component_variables("Square_1")
    print(variables)

p x y n id t L

Show list of these with geometry information#

[14]:

with mcstastox.Read(file_path) as loaded_data:
    loaded_data.show_components_with_geometry()

All components with geometry information in file:
All components with geometry information in file:
All components with geometry information in file:
All components with geometry information in file:
2026-04-30 13:35:42,953 INFO:All components with geometry information in file:
Square_1
Square_1
Square_1
Square_1
2026-04-30 13:35:42,955 INFO:Square_1
Square_2
Square_2
Square_2
Square_2
2026-04-30 13:35:42,957 INFO:Square_2
Banana_1
Banana_1
Banana_1
Banana_1
2026-04-30 13:35:42,959 INFO:Banana_1
Banana_2
Banana_2
Banana_2
Banana_2
2026-04-30 13:35:42,961 INFO:Banana_2

Show monitors with pixel ID’s#

[15]:

with mcstastox.Read(file_path) as loaded_data:
    loaded_data.show_components_with_ids()

All components with pixel id information in file:
All components with pixel id information in file:
All components with pixel id information in file:
All components with pixel id information in file:
All components with pixel id information in file:
2026-04-30 13:35:42,970 INFO:All components with pixel id information in file:
Square_1
Square_1
Square_1
Square_1
Square_1
2026-04-30 13:35:42,973 INFO:Square_1
Square_2
Square_2
Square_2
Square_2
Square_2
2026-04-30 13:35:42,974 INFO:Square_2
Banana_1
Banana_1
Banana_1
Banana_1
Banana_1
2026-04-30 13:35:42,977 INFO:Banana_1
Banana_2
Banana_2
Banana_2
Banana_2
Banana_2
2026-04-30 13:35:42,979 INFO:Banana_2

[16]:

with mcstastox.Read(file_path) as loaded_data:
    print(loaded_data.get_components_with_ids())

['Square_1', 'Square_2', 'Banana_1', 'Banana_2']

Get local pixel coordinates for a given component#

[17]:

with mcstastox.Read(file_path) as loaded_data:
    coordinates = loaded_data.get_component_local("Square_1")

plot(coordinates)

Get global coordinates for a given component#

[18]:

with mcstastox.Read(file_path) as loaded_data:
    coordinates = loaded_data.get_component_global("Square_1")

plot(coordinates)

Get sample position as well#

[19]:

with mcstastox.Read(file_path) as loaded_data:
    coordinates = loaded_data.get_component_global("Square_2")
    sample_pos = loaded_data.get_global_component_coordinates("sample_position")

plot(coordinates, points=[sample_pos])

Get pixel coordinates for list of monitors#

[20]:

with mcstastox.Read(file_path) as loaded_data:
    coordinates = loaded_data.get_id_to_global_coordinates(component_name=["Square_2", "Banana_1", "Banana_2"])
    sample_pos = loaded_data.get_global_component_coordinates("sample_position")

plot(coordinates, points=[sample_pos])

Include source and get all pixel coordinates#

[21]:

with mcstastox.Read(file_path) as loaded_data:
    coordinates = loaded_data.get_id_to_global_coordinates()
    sample_pos = loaded_data.get_global_component_coordinates("sample_position")
    source_pos = loaded_data.get_global_component_coordinates("source")

plot(coordinates, points=[sample_pos, source_pos])

Export to (simple) Scipp object#

[22]:

with mcstastox.Read(file_path) as loaded_data:
    scipp_object = loaded_data.export_scipp_simple(source_name="source",
                                                   sample_name="sample_position",
                                                   component_name="Banana_2")

[23]:

scipp_object

[23]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (736.17 KB)
    • events: 18799
    • position
      (events)
      vector3
      m
      [1.10723017 0.5501102 1.14614031], [1.08596947 0.67992982 1.2199218 ], ..., [0.98253646 0.77851101 1.17162874], [1.32632086 0.21235023 1.32292054]
      Values:
      array([[1.10723017, 0.5501102 , 1.14614031],
       [1.08596947, 0.67992982, 1.2199218 ],
       [0.92053824, 0.88841645, 1.21435899],
       ...,
       [1.14435646, 0.42467068, 1.12093269],
       [0.98253646, 0.77851101, 1.17162874],
       [1.32632086, 0.21235023, 1.32292054]], shape=(18799, 3))
    • sample_position
      ()
      vector3
      m
      [0. 0. 2.]
      Values:
      array([0., 0., 2.])
    • source_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • t
      (events)
      float64
      s
      0.001, 0.002, ..., 0.002, 0.002
      Values:
      array([0.00115242, 0.00161145, 0.00161684, ..., 0.00116698, 0.00162522,
       0.00164246], shape=(18799,))
    • (events)
      float64
      counts
      1.918e-11, 6.342e-12, ..., 3.986e-12, 3.280e-11
      Values:
      array([1.91779856e-11, 6.34159333e-12, 3.14056539e-12, ...,
       1.18500788e-11, 3.98586682e-12, 3.27998577e-11], shape=(18799,))

Plot event positions#

[24]:

import plopp as pp

pp.scatter3d(scipp_object[0::3], pos='position', size=0.02, cbar=True, norm="linear")

[24]:

Perform coordinate transforms#

[25]:

from scippneutron.conversion.graph.beamline import beamline
from scippneutron.conversion.graph.tof import elastic

# McStas provides absolute time, not time of flight
scipp_object.coords["tof"] = scipp_object.coords["t"]

graph = {**beamline(scatter=True), **elastic("tof")}

[26]:

scipp_object = scipp_object.transform_coords("dspacing", graph=graph)

[27]:

scipp_object.hist(dspacing=150).plot(norm="log")

[27]:

Load all data in simple scipp format#

[28]:

import mcstastox
with mcstastox.Read(file_path) as loaded_data:
    all_data = loaded_data.export_scipp_simple(source_name="source",
                                                   sample_name="sample_position")

all_data

[28]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (1.39 MB)
    • events: 36449
    • position
      (events)
      vector3
      m
      [-0.24540352 -0.02 1.74902545], [-0.21476174 -0.01333333 1.72331394], ..., [0.98253646 0.77851101 1.17162874], [1.32632086 0.21235023 1.32292054]
      Values:
      array([[-0.24540352, -0.02      ,  1.74902545],
       [-0.21476174, -0.01333333,  1.72331394],
       [-0.18922692,  0.04      ,  1.70188769],
       ...,
       [ 1.14435646,  0.42467068,  1.12093269],
       [ 0.98253646,  0.77851101,  1.17162874],
       [ 1.32632086,  0.21235023,  1.32292054]], shape=(36449, 3))
    • sample_position
      ()
      vector3
      m
      [0. 0. 2.]
      Values:
      array([0., 0., 2.])
    • source_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • t
      (events)
      float64
      s
      0.001, 0.001, ..., 0.002, 0.002
      Values:
      array([0.00076284, 0.00116951, 0.00083817, ..., 0.00116698, 0.00162522,
       0.00164246], shape=(36449,))
    • (events)
      float64
      counts
      3.505e-12, 1.322e-11, ..., 3.986e-12, 3.280e-11
      Values:
      array([3.50488533e-12, 1.32208569e-11, 1.60747884e-11, ...,
       1.18500788e-11, 3.98586682e-12, 3.27998577e-11], shape=(36449,))
[29]:

pp.scatter3d(all_data[0::3], pos='position', size=0.02, cbar=True, norm="log")

[29]:

Convert to d-space and compare#

Lets compare data from all the detectors with the previous scipp object that only contained one.

[30]:

all_data.coords["tof"] = all_data.coords["t"]
all_data = all_data.transform_coords("dspacing", graph=graph)

pp.plot({"Only banana 2": scipp_object.hist(dspacing=200),
         "All data": all_data.hist(dspacing=200)}, norm="log")

[30]:

Also access to data directly#

Its possible to get access to the data directly if necessary, this is demonstrated here.

[31]:

with mcstastox.Read(file_path) as loaded_data:
    raw_event_data = loaded_data.get_event_data(variables=["id", "t"], component_name="Square_2")
    pos, rot = loaded_data.get_component_placement(component_name="Square_2")

print("Raw event data")
print(raw_event_data["id"])
print(raw_event_data["t"])

print("Component placement")
print(pos)
print(rot)

Raw event data
[555. 357. 228. ... 355. 377. 594.]
[0.00122248 0.00078535 0.00114255 ... 0.00108291 0.00110442 0.00116819]
Component placement
[ 0.41933528 -0.09313876  2.25589664]
[[ 0.54463904  0.28684223 -0.78809254]
 [-0.          0.93969262  0.34202014]
 [ 0.83867057 -0.18627752  0.51179328]]

[32]:

with mcstastox.Read(file_path) as loaded_data:
    raw_event_data = loaded_data.get_event_data(variables=["id", "t", "L"], component_name="Square_1")

Add additional data to a scipp object#

We can add data read directly to the scipp object, here the actual wavelength recorded by McStas.

[33]:

import scipp as sc

with mcstastox.Read(file_path) as loaded_data:
    sq1 = loaded_data.export_scipp_simple(source_name="source", sample_name="sample_position", component_name="Square_1")

sq1.coords["sim_wavelength"] = sc.array(dims=["events"], values=raw_event_data["L"], unit="Å")

sq1.coords["tof"] = sq1.coords["t"]
sq1 = sq1.transform_coords("wavelength", graph=graph)
sq1.coords["wavelength_ratio"] = (sq1.coords["sim_wavelength"] - sq1.coords["wavelength"])/sq1.coords["sim_wavelength"]

[34]:

sq1

[34]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (460.88 KB)
    • events: 4175
    • L1
      ()
      float64
      m
      2.0
      Values:
      array(2.)
    • L2
      (events)
      float64
      m
      0.352, 0.351, ..., 0.353, 0.351
      Values:
      array([0.35158372, 0.35050757, 0.35535585, ..., 0.35215211, 0.35259357,
       0.35057096], shape=(4175,))
    • Ltotal
      (events)
      float64
      m
      2.352, 2.351, ..., 2.353, 2.351
      Values:
      array([2.35158372, 2.35050757, 2.35535585, ..., 2.35215211, 2.35259357,
       2.35057096], shape=(4175,))
    • incident_beam
      ()
      vector3
      m
      [0. 0. 2.]
      Values:
      array([0., 0., 2.])
    • position
      (events)
      vector3
      m
      [-0.24540352 -0.02 1.74902545], [-0.21476174 -0.01333333 1.72331394], ..., [-0.25051048 -0.02666667 1.7533107 ], [-0.22497566 0.02 1.73188444]
      Values:
      array([[-0.24540352, -0.02      ,  1.74902545],
       [-0.21476174, -0.01333333,  1.72331394],
       [-0.18922692,  0.04      ,  1.70188769],
       ...,
       [-0.25051048,  0.02      ,  1.7533107 ],
       [-0.25051048, -0.02666667,  1.7533107 ],
       [-0.22497566,  0.02      ,  1.73188444]], shape=(4175, 3))
    • sample_position
      ()
      vector3
      m
      [0. 0. 2.]
      Values:
      array([0., 0., 2.])
    • scattered_beam
      (events)
      vector3
      m
      [-0.24540352 -0.02 -0.25097455], [-0.21476174 -0.01333333 -0.27668606], ..., [-0.25051048 -0.02666667 -0.2466893 ], [-0.22497566 0.02 -0.26811556]
      Values:
      array([[-0.24540352, -0.02      , -0.25097455],
       [-0.21476174, -0.01333333, -0.27668606],
       [-0.18922692,  0.04      , -0.29811231],
       ...,
       [-0.25051048,  0.02      , -0.2466893 ],
       [-0.25051048, -0.02666667, -0.2466893 ],
       [-0.22497566,  0.02      , -0.26811556]], shape=(4175, 3))
    • sim_wavelength
      (events)
      float64
      Å
      1.285, 1.974, ..., 2.006, 1.521
      Values:
      array([1.28487592, 1.97383291, 1.41095847, ..., 2.94379194, 2.00560265,
       1.52111715], shape=(4175,))
    • source_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • t
      (events)
      float64
      s
      0.001, 0.001, ..., 0.001, 0.001
      Values:
      array([0.00076284, 0.00116951, 0.00083817, ..., 0.00174814, 0.00118902,
       0.00090694], shape=(4175,))
    • tof
      (events)
      float64
      s
      0.001, 0.001, ..., 0.001, 0.001
      Values:
      array([0.00076284, 0.00116951, 0.00083817, ..., 0.00174814, 0.00118902,
       0.00090694], shape=(4175,))
    • wavelength
      (events)
      float64
      Å
      1.283, 1.968, ..., 1.999, 1.526
      Values:
      array([1.28331066, 1.96834739, 1.40778214, ..., 2.94016238, 1.99941086,
       1.52639268], shape=(4175,))
    • wavelength_ratio
      (events)
      float64
      𝟙
      0.001, 0.003, ..., 0.003, -0.003
      Values:
      array([ 0.00121821,  0.00277912,  0.00225119, ...,  0.00123295,
        0.00308725, -0.0034682 ], shape=(4175,))
    • (events)
      float64
      counts
      3.505e-12, 1.322e-11, ..., 5.346e-11, 8.130e-12
      Values:
      array([3.50488533e-12, 1.32208569e-11, 1.60747884e-11, ...,
       1.21539814e-11, 5.34598947e-11, 8.13010044e-12], shape=(4175,))

Plotting comparison between actual and calculated wavelength#

[35]:

w_start = 0.5
w_end=3.2
w_steps=401
wavs = sc.linspace("wavelength", w_start, w_end, w_steps, unit="angstrom")
sim_wavs = sc.linspace("sim_wavelength", w_start, w_end, w_steps, unit="angstrom")

pp.plot({
    "calculated wavelength": sq1.hist(wavelength=wavs),
    "true wavelength": sq1.hist(sim_wavelength=sim_wavs).rename(sim_wavelength="wavelength")
}) / sq1.hist(wavelength_ratio=200).plot()

[35]: