McStasRead demo - with scipp#

[1]:

%matplotlib widget

[2]:

from read_example import make_instrument

[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_0"
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/cchfw6RW.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_0
INFO: Placing generated c-code copy test.c in dataset /home/runner/work/McStasToX/McStasToX/docs/user-guide/test_0


[7]:

data

[7]:

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

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

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

 McStasDataEvent: Square_2 with 7913 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_0

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:56,513 INFO:All components in file:
source
2026-04-30 13:35:56,514 INFO:source
sample_position
2026-04-30 13:35:56,514 INFO:sample_position
sample
2026-04-30 13:35:56,515 INFO:sample
detector_direction_square_1
2026-04-30 13:35:56,516 INFO:detector_direction_square_1
Square_1
2026-04-30 13:35:56,517 INFO:Square_1
detector_direction_square_2
2026-04-30 13:35:56,517 INFO:detector_direction_square_2
Square_2
2026-04-30 13:35:56,518 INFO:Square_2
detector_direction_banana_1
2026-04-30 13:35:56,519 INFO:detector_direction_banana_1
Banana_1
2026-04-30 13:35:56,519 INFO:Banana_1
detector_direction_banana_2
2026-04-30 13:35:56,520 INFO:detector_direction_banana_2
Banana_2
2026-04-30 13:35:56,520 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:56,527 INFO:All components with data in file:
Square_1
Square_1
2026-04-30 13:35:56,528 INFO:Square_1
Square_2
Square_2
2026-04-30 13:35:56,529 INFO:Square_2
Banana_1
Banana_1
2026-04-30 13:35:56,530 INFO:Banana_1
Banana_2
Banana_2
2026-04-30 13:35:56,531 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:56,545 INFO:All components with geometry information in file:
Square_1
Square_1
Square_1
Square_1
2026-04-30 13:35:56,546 INFO:Square_1
Square_2
Square_2
Square_2
Square_2
2026-04-30 13:35:56,548 INFO:Square_2
Banana_1
Banana_1
Banana_1
Banana_1
2026-04-30 13:35:56,549 INFO:Banana_1
Banana_2
Banana_2
Banana_2
Banana_2
2026-04-30 13:35:56,550 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:56,561 INFO:All components with pixel id information in file:
Square_1
Square_1
Square_1
Square_1
Square_1
2026-04-30 13:35:56,563 INFO:Square_1
Square_2
Square_2
Square_2
Square_2
Square_2
2026-04-30 13:35:56,566 INFO:Square_2
Banana_1
Banana_1
Banana_1
Banana_1
Banana_1
2026-04-30 13:35:56,568 INFO:Banana_1
Banana_2
Banana_2
Banana_2
Banana_2
Banana_2
2026-04-30 13:35:56,570 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']

Export to Scipp object#

Using the export_scipp method we get a scipp DataGroup that contains:

  • events : the data

  • positions : positions of the pixel ids

  • bank_ids : pixel id range for each detector bank

  • bank_names : names of the loaded detector banks

This structure requires a little more knowledge to work with than the simple export, but saves space and provides more flexibility.

[17]:

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

[18]:

scipp_object

[18]:
  • events
    scipp
    DataArray
    (pixel_id: 1115)
    float64
    counts
    binned data [len=14, len=17, ..., len=111, len=100]
  • positions
    scipp
    Variable
    (pixel_id: 1115)
    vector3
    m
    [-0.18922692 -0.04666667 1.70188769], [-0.19433389 -0.04666667 1.70617294], ..., [0.92516053 1.01752513 1.39059304], [0.83699271 1.10865008 1.42294328]
  • bank_ids
    scipp
    Variable
    (panel_id: 4, pixel: 2)
    int64
    𝟙
    0, 224, ..., 915, 1114
  • bank_names
    scipp
    Variable
    (panel_id: 4)
    string
    Square_1, Square_2, Banana_1, Banana_2
[19]:

scipp_object["bank_names"]

[19]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (384 Bytes)
    • (panel_id: 4)
      string
      Square_1, Square_2, Banana_1, Banana_2
      Values:
      ["Square_1", "Square_2", "Banana_1", "Banana_2"]
[20]:

scipp_object["bank_ids"]

[20]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (320 Bytes)
    • (panel_id: 4, pixel: 2)
      int64
      𝟙
      0, 224, ..., 915, 1114
      Values:
      array([[   0,  224],
       [ 225,  674],
       [ 675,  914],
       [ 915, 1114]])
[21]:

scipp_object["positions"]

[21]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (26.44 KB)
    • (pixel_id: 1115)
      vector3
      m
      [-0.18922692 -0.04666667 1.70188769], [-0.19433389 -0.04666667 1.70617294], ..., [0.92516053 1.01752513 1.39059304], [0.83699271 1.10865008 1.42294328]
      Values:
      array([[-0.18922692, -0.04666667,  1.70188769],
       [-0.19433389, -0.04666667,  1.70617294],
       [-0.19944085, -0.04666667,  1.71045819],
       ...,
       [ 1.00670892,  0.91880488,  1.36361098],
       [ 0.92516053,  1.01752513,  1.39059304],
       [ 0.83699271,  1.10865008,  1.42294328]], shape=(1115, 3))
[22]:

scipp_object["events"]

[22]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (622.51 KB)
    • pixel_id: 1115
    • pixel_id
      (pixel_id)
      int64
      𝟙
      0, 1, ..., 1113, 1114
      Values:
      array([   0,    1,    2, ..., 1112, 1113, 1114], shape=(1115,))
    • position
      (pixel_id)
      vector3
      m
      [-0.18922692 -0.04666667 1.70188769], [-0.19433389 -0.04666667 1.70617294], ..., [0.92516053 1.01752513 1.39059304], [0.83699271 1.10865008 1.42294328]
      Values:
      array([[-0.18922692, -0.04666667,  1.70188769],
       [-0.19433389, -0.04666667,  1.70617294],
       [-0.19944085, -0.04666667,  1.71045819],
       ...,
       [ 1.00670892,  0.91880488,  1.36361098],
       [ 0.92516053,  1.01752513,  1.39059304],
       [ 0.83699271,  1.10865008,  1.42294328]], shape=(1115, 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.])
    • (pixel_id)
      float64
      counts
      binned data [len=14, len=17, ..., len=111, len=100]
      dim='events',
content=DataArray(
          dims=(events: 36315),
          data=float64[counts],
          coords={'t':float64[s]})

Plot pixels with intensities#

With this setup we can plot the total intensity in each pixel rather than all events individually.

[23]:

import plopp as pp

pp.scatter3d(scipp_object["events"].hist(), pos='position', size=0.015, cbar=True, norm="log")

[23]:

Perform coordinate transforms#

Coordinate transformations can be done almost as before, they just need to be summed over all pixel_id when plotted.

[24]:

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

event_object = scipp_object["events"]

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

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

[25]:

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

Note extra .sum(“pixel_id”) in plot.

[26]:

event_object.hist(dspacing=150).sum("pixel_id").plot(norm="log")

[26]: