Merge pull request #162 from scipp/fix-estia-detector

fix: shape and size of the estia detector

Author: jokasimr

docs/user-guide/estia/estia-mcstas-reduction.ipynb

@@ -57,7 +57,7 @@
     "wf[Filename[ReferenceRun]] = estia_mcstas_example('reference')\n",
     "\n",
     "wf[YIndexLimits]  = sc.scalar(35), sc.scalar(64)\n",
-    "wf[ZIndexLimits] = sc.scalar(0), sc.scalar(14 * 32)\n",
+    "wf[ZIndexLimits] = sc.scalar(0), sc.scalar(48 * 32)\n",
     "wf[BeamDivergenceLimits] = sc.scalar(-0.75, unit='deg'), sc.scalar(0.75, unit='deg')\n",
     "wf[WavelengthBins] = sc.geomspace('wavelength', 3.5, 12, 2001, unit='angstrom')\n",
     "wf[QBins] = 1000\n",

src/ess/amor/conversions.py

@@ -163,7 +163,7 @@ def add_masks(
     and by wavelength.
     """
     da = da.assign_masks(
-        stripe_range=_not_between(da.coords["stripe"], *ylim),
+        strip_range=_not_between(da.coords["strip"], *ylim),
         z_range=_not_between(da.coords["z_index"], *zlims),
     )
     da = da.bins.assign_masks(

src/ess/amor/geometry.py

@@ -10,8 +10,8 @@ class Detector:
     nBlades = sc.scalar(14)
     # number of wires per blade
     nWires = sc.scalar(32)
-    # number of stripes per blade
-    nStripes = sc.scalar(64)
+    # number of strips per blade
+    nStrips = sc.scalar(64)
     # angle of incidence of the beam on the blades (def: 5.1)
     angle = sc.scalar(5.1, unit="degree").to(unit="rad")
     # height-distance of neighboring pixels on one blade
@@ -34,21 +34,21 @@ def pixel_coordinates_in_detector_system() -> tuple[sc.Variable, sc.Variable]:
             'row',
             1,
             (
-                Detector.nBlades * Detector.nWires * Detector.nStripes + sc.scalar(1)
+                Detector.nBlades * Detector.nWires * Detector.nStrips + sc.scalar(1)
             ).values,
             unit=None,
         ).fold(
             'row',
             sizes={
                 'blade': Detector.nBlades,
                 'wire': Detector.nWires,
-                'stripe': Detector.nStripes,
+                'strip': Detector.nStrips,
             },
         ),
         coords={
             'blade': sc.arange('blade', sc.scalar(0), Detector.nBlades),
             'wire': sc.arange('wire', sc.scalar(0), Detector.nWires),
-            'stripe': sc.arange('stripe', sc.scalar(0), Detector.nStripes),
+            'strip': sc.arange('strip', sc.scalar(0), Detector.nStrips),
         },
     )
     # x position in detector

src/ess/estia/beamline.py

@@ -0,0 +1,24 @@
+import sciline
+
+from ess.reduce.nexus.types import DetectorBankSizes
+from ess.reduce.nexus.workflow import GenericNeXusWorkflow
+from ess.reflectometry.types import ReferenceRun, SampleRun
+
+DETECTOR_BANK_SIZES = {
+    "multiblade_detector": {
+        "blade": 48,
+        "wire": 32,
+        "strip": 64,
+    },
+}
+
+
+def LoadNeXusWorkflow() -> sciline.Pipeline:
+    """
+    Workflow for loading NeXus data.
+    """
+    workflow = GenericNeXusWorkflow(
+        run_types=[SampleRun, ReferenceRun], monitor_types=[]
+    )
+    workflow[DetectorBankSizes] = DETECTOR_BANK_SIZES
+    return workflow

src/ess/estia/load.py

@@ -9,6 +9,7 @@
     RunType,
     SampleRotationOffset,
 )
+from .beamline import DETECTOR_BANK_SIZES
 from .mcstas import parse_events_ascii, parse_events_h5
 
 
@@ -38,11 +39,14 @@ def load_mcstas_events(
         unit=da.coords['sample_rotation'].unit
     )
 
-    xbins = sc.linspace('x', -0.25, 0.25, 14 * 32 + 1)
-    ybins = sc.linspace('y', -0.25, 0.25, 65)
-    da = da.bin(x=xbins, y=ybins).rename_dims({'y': 'stripe'})
-    da.coords['stripe'] = sc.arange('stripe', 0, 64)
-    da.coords['z_index'] = sc.arange('x', 14 * 32 - 1, -1, -1)
+    nblades = DETECTOR_BANK_SIZES['multiblade_detector']['blade']
+    nwires = DETECTOR_BANK_SIZES['multiblade_detector']['wire']
+    nstrips = DETECTOR_BANK_SIZES['multiblade_detector']['strip']
+    xbins = sc.linspace('x', -0.25, 0.25, nblades * nwires + 1)
+    ybins = sc.linspace('y', -0.13, 0.13, nstrips + 1)
+    da = da.bin(y=ybins, x=xbins).rename_dims({'y': 'strip'})
+    da.coords['strip'] = sc.arange('strip', 0, nstrips)
+    da.coords['z_index'] = sc.arange('x', nblades * nwires - 1, -1, -1)
 
     # Information is not available in the mcstas output files, therefore it's hardcoded
     da.coords['sample_position'] = sc.vector([0.264298, -0.427595, 35.0512], unit='m')
@@ -64,7 +68,7 @@ def load_mcstas_events(
         x=sc.midpoints(da.coords['x']) * sc.scalar(1.0, unit='m'),
         y=sc.midpoints(da.coords['y']) * sc.scalar(1.0, unit='m'),
         z=sc.scalar(0.0, unit='m'),
-    )
+    ).transpose(da.dims)
     da.coords['position'] = (
         da.coords['detector_position'] + rotation_by_detector_rotation * position
     )
@@ -82,7 +86,14 @@ def load_mcstas_events(
     )
     da.coords["beam_size"] = sc.scalar(2.0, unit='mm')
 
-    da = da.fold('x', sizes={'blade': 14, 'wire': 32})
+    da = da.fold(
+        'x',
+        sizes={
+            k: v
+            for k, v in DETECTOR_BANK_SIZES['multiblade_detector'].items()
+            if k in ('blade', 'wire')
+        },
+    )
     da.bins.coords.pop('L')
     da.bins.coords.pop('t')
     return DetectorData[RunType](da)

src/ess/estia/maskings.py

@@ -27,7 +27,7 @@ def add_masks(
     and by wavelength.
     """
     da = da.assign_masks(
-        stripe_range=_not_between(da.coords["stripe"], *ylim),
+        strip_range=_not_between(da.coords["strip"], *ylim),
         z_range=_not_between(da.coords["z_index"], *zlims),
         divergence_too_large=_not_between(
             da.coords["divergence_angle"],

src/ess/estia/workflow.py

@@ -4,20 +4,16 @@
 import sciline
 import scipp as sc
 
-from ess.reduce import nexus
-
 from ..reflectometry import providers as reflectometry_providers
 from ..reflectometry import supermirror
 from ..reflectometry.types import (
     BeamDivergenceLimits,
     DetectorSpatialResolution,
     NeXusDetectorName,
-    ReferenceRun,
     RunType,
     SampleRotationOffset,
-    SampleRun,
 )
-from . import conversions, corrections, load, maskings, normalization, orso
+from . import beamline, conversions, corrections, load, maskings, normalization, orso
 
 _general_providers = (
     *reflectometry_providers,
@@ -74,9 +70,7 @@ def EstiaMcStasWorkflow() -> sciline.Pipeline:
 
 def EstiaWorkflow() -> sciline.Pipeline:
     """Workflow for reduction of data for the Estia instrument."""
-    workflow = nexus.GenericNeXusWorkflow(
-        run_types=[SampleRun, ReferenceRun], monitor_types=[]
-    )
+    workflow = beamline.LoadNeXusWorkflow()
     for provider in providers:
         workflow.insert(provider)
     for name, param in default_parameters().items():

src/ess/reflectometry/normalization.py

@@ -40,10 +40,10 @@ def reduce_reference(
     )
     reference = reference.bins.assign_masks(invalid=sc.isnan(R))
     reference = reference / R
-    out = reference.bins.concat(('stripe',)).hist(wavelength=wavelength_bins)
+    out = reference.bins.concat(('strip',)).hist(wavelength=wavelength_bins)
 
     if 'position' in reference.coords:
-        out.coords['position'] = reference.coords['position'].mean('stripe')
+        out.coords['position'] = reference.coords['position'].mean('strip')
     return out
 
 
@@ -96,7 +96,7 @@ def reduce_sample_over_zw(
 
     Returns reflectivity as a function of ``blade``, ``wire`` and :math:`\\wavelength`.
     """
-    return sample.bins.concat(('stripe',)).bin(wavelength=wbins) / sc.values(
+    return sample.bins.concat(('strip',)).bin(wavelength=wbins) / sc.values(
         reference.data
     )
 

tests/estia/mcstas_data_test.py

@@ -37,7 +37,7 @@ def estia_mcstas_pipeline() -> sciline.Pipeline:
     wf[Filename[ReferenceRun]] = estia_mcstas_reference_run()
 
     wf[YIndexLimits] = sc.scalar(35), sc.scalar(64)
-    wf[ZIndexLimits] = sc.scalar(0), sc.scalar(14 * 32)
+    wf[ZIndexLimits] = sc.scalar(0), sc.scalar(48 * 32)
     wf[BeamDivergenceLimits] = sc.scalar(-1.0, unit='deg'), sc.scalar(1.0, unit='deg')
     wf[WavelengthBins] = sc.geomspace('wavelength', 3.5, 12, 2001, unit='angstrom')
     wf[QBins] = sc.geomspace('Q', 0.005, 0.1, 200, unit='1/angstrom')
@@ -79,8 +79,8 @@ def estia_mcstas_pipeline() -> sciline.Pipeline:
 def test_mcstas_compute_reducible_data(estia_mcstas_pipeline: sciline.Pipeline):
     estia_mcstas_pipeline[Filename[SampleRun]] = estia_mcstas_sample_run(11)
     da = estia_mcstas_pipeline.compute(ReducibleData[SampleRun])
-    assert da.dims == ('blade', 'wire', 'stripe')
-    assert da.shape == (14, 32, 64)
+    assert da.dims == ('strip', 'blade', 'wire')
+    assert da.shape == (64, 48, 32)
     assert 'position' in da.coords
     assert 'sample_rotation' in da.coords
     assert 'detector_rotation' in da.coords

tests/reflectometry/normalization_test.py

@@ -15,10 +15,10 @@ def sample(request):
         coords={
             'wavelength': sc.linspace('events', 1, 5, n),
             'wire': sc.array(dims=('events',), values=np.random.randint(0, 5, n)),
-            'stripe': sc.array(dims=('events',), values=np.random.randint(0, 10, n)),
+            'strip': sc.array(dims=('events',), values=np.random.randint(0, 10, n)),
         },
     )
-    return da.group('wire', 'stripe')
+    return da.group('wire', 'strip')
 
 
 @pytest.fixture
@@ -29,7 +29,7 @@ def reference(request):
         coords={
             'wavelength': sc.linspace('events', 1, 5, n),
             'wire': sc.array(dims=('events',), values=np.random.randint(0, 5, n)),
-            'stripe': sc.array(dims=('events',), values=np.random.randint(0, 10, n)),
+            'strip': sc.array(dims=('events',), values=np.random.randint(0, 10, n)),
         },
     )
     return da.group('wire').bin(wavelength=2).bins.sum()