diff --git a/docs/conf.py b/docs/conf.py index 17b45932..37ffb8fa 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -5,35 +5,35 @@ import sys from importlib.metadata import version as get_version -sys.path.insert(0, os.path.abspath('.')) +sys.path.insert(0, os.path.abspath(".")) # General information about the project. -project = 'ESSdiffraction' -copyright = '2024 Scipp contributors' -author = 'Scipp contributors' +project = "ESSdiffraction" +copyright = "2024 Scipp contributors" +author = "Scipp contributors" html_show_sourcelink = True extensions = [ - 'sphinx.ext.autodoc', - 'sphinx.ext.autosummary', - 'sphinx.ext.doctest', - 'sphinx.ext.githubpages', - 'sphinx.ext.intersphinx', - 'sphinx.ext.mathjax', - 'sphinx.ext.napoleon', - 'sphinx.ext.viewcode', - 'sphinx_autodoc_typehints', - 'sphinx_copybutton', - 'sphinx_design', - 'nbsphinx', - 'myst_parser', + "sphinx.ext.autodoc", + "sphinx.ext.autosummary", + "sphinx.ext.doctest", + "sphinx.ext.githubpages", + "sphinx.ext.intersphinx", + "sphinx.ext.mathjax", + "sphinx.ext.napoleon", + "sphinx.ext.viewcode", + "sphinx_autodoc_typehints", + "sphinx_copybutton", + "sphinx_design", + "nbsphinx", + "myst_parser", ] try: import sciline.sphinxext.domain_types # noqa: F401 - extensions.append('sciline.sphinxext.domain_types') + extensions.append("sciline.sphinxext.domain_types") except ModuleNotFoundError: pass @@ -56,13 +56,13 @@ myst_heading_anchors = 3 autodoc_type_aliases = { - 'array_like': 'array_like', + "array_like": "array_like", } intersphinx_mapping = { - 'python': ('https://docs.python.org/3', None), - 'numpy': ('https://numpy.org/doc/stable/', None), - 'scipp': ('https://scipp.github.io/', None), + "python": ("https://docs.python.org/3", None), + "numpy": ("https://numpy.org/doc/stable/", None), + "scipp": ("https://scipp.github.io/", None), } # autodocs includes everything, even irrelevant API internals. autosummary @@ -79,32 +79,32 @@ # objects without namespace: numpy "ndarray": "~numpy.ndarray", } -typehints_defaults = 'comma' +typehints_defaults = "comma" typehints_use_rtype = False -sciline_domain_types_prefix = 'ess.diffraction' +sciline_domain_types_prefix = "ess.diffraction" sciline_domain_types_aliases = { - 'scipp._scipp.core.DataArray': 'scipp.DataArray', - 'scipp._scipp.core.Dataset': 'scipp.Dataset', - 'scipp._scipp.core.DType': 'scipp.DType', - 'scipp._scipp.core.Unit': 'scipp.Unit', - 'scipp._scipp.core.Variable': 'scipp.Variable', - 'scipp.core.data_group.DataGroup': 'scipp.DataGroup', + "scipp._scipp.core.DataArray": "scipp.DataArray", + "scipp._scipp.core.Dataset": "scipp.Dataset", + "scipp._scipp.core.DType": "scipp.DType", + "scipp._scipp.core.Unit": "scipp.Unit", + "scipp._scipp.core.Variable": "scipp.Variable", + "scipp.core.data_group.DataGroup": "scipp.DataGroup", } # Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] +templates_path = ["_templates"] # The suffix(es) of source filenames. # You can specify multiple suffix as a list of string: # -source_suffix = ['.rst', '.md'] -html_sourcelink_suffix = '' # Avoid .ipynb.txt extensions in sources +source_suffix = [".rst", ".md"] +html_sourcelink_suffix = "" # Avoid .ipynb.txt extensions in sources # The master toctree document. -master_doc = 'index' +master_doc = "index" # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the @@ -112,7 +112,7 @@ # release = get_version("essdiffraction") -version = ".".join(release.split('.')[:3]) # CalVer +version = ".".join(release.split(".")[:3]) # CalVer warning_is_error = True @@ -127,15 +127,15 @@ # directories to ignore when looking for source files. # This patterns also effect to html_static_path and html_extra_path exclude_patterns = [ - '_build', - 'Thumbs.db', - '.DS_Store', - '**.ipynb_checkpoints', - 'user-guide/sns-instruments/preprocess_files.ipynb', + "_build", + "Thumbs.db", + ".DS_Store", + "**.ipynb_checkpoints", + "user-guide/sns-instruments/preprocess_files.ipynb", ] # The name of the Pygments (syntax highlighting) style to use. -pygments_style = 'sphinx' +pygments_style = "sphinx" # If true, `todo` and `todoList` produce output, else they produce nothing. todo_include_todos = False @@ -200,14 +200,14 @@ # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] +html_static_path = ["_static"] html_css_files = [] html_js_files = ["anaconda-icon.js"] # -- Options for HTMLHelp output ------------------------------------------ # Output file base name for HTML help builder. -htmlhelp_basename = 'essdiffractiondoc' +htmlhelp_basename = "essdiffractiondoc" # -- Options for Matplotlib in notebooks ---------------------------------- @@ -223,7 +223,7 @@ # In addition, there is no need to make plots in doctest as the documentation # build already tests if those plots can be made. # So we simply disable plots in doctests. -doctest_global_setup = ''' +doctest_global_setup = """ import numpy as np try: @@ -240,7 +240,7 @@ def do_not_plot(*args, **kwargs): except ImportError: # Scipp is not needed by docs if it is not installed. pass -''' +""" # Using normalize whitespace because many __str__ functions in scipp produce # extraneous empty lines and it would look strange to include them in the docs. @@ -255,5 +255,5 @@ def do_not_plot(*args, **kwargs): linkcheck_ignore = [ # Specific lines in Github blobs cannot be found by linkcheck. - r'https?://github\.com/.*?/blob/[a-f0-9]+/.+?#', + r"https?://github\.com/.*?/blob/[a-f0-9]+/.+?#", ] diff --git a/docs/user-guide/dream/dream-data-reduction.ipynb b/docs/user-guide/dream/dream-data-reduction.ipynb index 06cd8586..92bad4de 100644 --- a/docs/user-guide/dream/dream-data-reduction.ipynb +++ b/docs/user-guide/dream/dream-data-reduction.ipynb @@ -23,6 +23,17 @@ "from ess.dream.io.geant4 import providers as geant4_providers" ] }, + { + "cell_type": "markdown", + "id": "1252feab-12d2-46ac-bf74-70b32344473d", + "metadata": {}, + "source": [ + "## Define reduction parameters\n", + "\n", + "We define a dictionary containing the reduction parameters.\n", + "The keys are types defined in [essdiffraction.types](../generated/modules/ess.diffraction.types.rst)." + ] + }, { "cell_type": "code", "execution_count": null, @@ -30,33 +41,42 @@ "metadata": {}, "outputs": [], "source": [ - "providers = (\n", - " *geant4_providers,\n", - " *powder.providers,\n", - " *dream.providers,\n", - ")\n", "params = {\n", - " FilePath[SampleRun]: dream.data.simulated_diamond_sample(),\n", - " FilePath[VanadiumRun]: dream.data.simulated_vanadium_sample(),\n", - " FilePath[EmptyCanRun]: dream.data.simulated_empty_can(),\n", - " DetectorName: DetectorName('mantle'),\n", - "\n", + " Filename[SampleRun]: dream.data.simulated_diamond_sample(),\n", + " Filename[VanadiumRun]: dream.data.simulated_vanadium_sample(),\n", + " Filename[EmptyCanRun]: dream.data.simulated_empty_can(),\n", + " NeXusDetectorName: \"mantle\",\n", " # The upper bounds mode is not yet implemented.\n", " UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,\n", + " # Edges for binning in d-spacing\n", + " DspacingBins: sc.linspace(\"dspacing\", 0.0, 2.3434, 201, unit=\"angstrom\"),\n", + " # Mask in time-of-flight to crop to valid range\n", + " TofMask: lambda x: (x < sc.scalar(0.0, unit=\"ns\"))\n", + " | (x > sc.scalar(86e6, unit=\"ns\")),\n", + "}\n", "\n", - " # Crop data to this range in time-of-flight\n", - " ValidTofRange: sc.array(dims=['tof'], values=[0.0, 86_000.0], unit='us'),\n", + "# Not available in simulated data\n", + "sample = sc.DataGroup(position=sc.vector([0.0, 0.0, 0.0], unit=\"mm\"))\n", + "params[RawSample[SampleRun]] = sample\n", + "params[RawSample[VanadiumRun]] = sample\n", "\n", - " # Edges for binning in d-spacing\n", - " DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 200, unit='angstrom'),\n", + "source = sc.DataGroup(position=sc.vector([-3.478, 0.0, -76550], unit=\"mm\"))\n", + "params[RawSource[SampleRun]] = source\n", + "params[RawSource[VanadiumRun]] = source\n", + "\n", + "charge = sc.scalar(1.0, unit=\"µAh\")\n", + "params[AccumulatedProtonCharge[SampleRun]] = charge\n", + "params[AccumulatedProtonCharge[VanadiumRun]] = charge" + ] + }, + { + "cell_type": "markdown", + "id": "21cb87f2-4ff7-436e-b603-cc8f60c73e7a", + "metadata": {}, + "source": [ + "## Create pipeline using Sciline\n", "\n", - " # Not available in simulated data\n", - " RawSample[SampleRun]: sc.DataGroup(position=sc.vector([0., 0., 0.], unit='mm')),\n", - " RawSample[VanadiumRun]: sc.DataGroup(position=sc.vector([0., 0., 0.], unit='mm')),\n", - " RawSource: sc.DataGroup(position=sc.vector([-3.478, 0.0, -76550], unit='mm')),\n", - " AccumulatedProtonCharge[SampleRun]: sc.scalar(1.0, unit='µAh'),\n", - " AccumulatedProtonCharge[VanadiumRun]: sc.scalar(1.0, unit='µAh'),\n", - "}" + "We use the `powder` and `geant4` providers to build our pipeline." ] }, { @@ -66,6 +86,11 @@ "metadata": {}, "outputs": [], "source": [ + "providers = (\n", + " *geant4_providers,\n", + " *powder.providers,\n", + ")\n", + "\n", "pipeline = sciline.Pipeline(providers, params=params)" ] }, @@ -74,70 +99,171 @@ "id": "21fb4492-e836-41d3-a2d4-9678df43b9f9", "metadata": {}, "source": [ - "We don't have calibration data yet.\n", - "So convert to d-spacing without calibration." + "We can visualize the graph for computing the final normalized result for intensity as a function of d-spacing:" ] }, { "cell_type": "code", "execution_count": null, - "id": "3dd37911-7173-413f-a30b-362da8e7f326", + "id": "93f2ba0d-f256-4b64-b8fa-42cd1081cc41", "metadata": {}, "outputs": [], "source": [ - "from ess.powder.conversion import to_dspacing_with_positions\n", - "pipeline.insert(to_dspacing_with_positions)" + "pipeline.visualize(IofDspacing, graph_attr={\"rankdir\": \"LR\"})" + ] + }, + { + "cell_type": "markdown", + "id": "478f580e-273c-4a46-8dfc-9a2ed31b0cd3", + "metadata": {}, + "source": [ + "We then call `compute()` to compute the result:" ] }, { "cell_type": "code", "execution_count": null, - "id": "93f2ba0d-f256-4b64-b8fa-42cd1081cc41", + "id": "6a17cf1d-0407-41dd-8a84-0975371ac70a", "metadata": {}, "outputs": [], "source": [ - "pipeline.visualize(DspacingHistogram, graph_attr={'rankdir': 'LR'})" + "result = pipeline.compute(IofDspacing)\n", + "result" ] }, { "cell_type": "code", "execution_count": null, - "id": "6a17cf1d-0407-41dd-8a84-0975371ac70a", + "id": "527623d2-0eee-456d-a4dd-764df8618ee4", "metadata": {}, "outputs": [], "source": [ - "result = pipeline.compute(DspacingHistogram)" + "dspacing_histogram = result.hist()\n", + "dspacing_histogram.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "39dc0ea3-b30d-4a88-878c-86cf0e8e76e9", + "metadata": {}, + "source": [ + "We can now save the result to disk:" ] }, { "cell_type": "code", "execution_count": null, - "id": "c05859ef-753b-4ea6-8d90-a17aa4a58c8f", + "id": "f37eba41-6178-4bcb-bf3c-fa177e6112b0", "metadata": {}, "outputs": [], "source": [ - "result" + "dspacing_histogram.coords[\"dspacing\"] = sc.midpoints(\n", + " dspacing_histogram.coords[\"dspacing\"]\n", + ")\n", + "scn.io.save_xye(\"dspacing.xye\", dspacing_histogram)" + ] + }, + { + "cell_type": "markdown", + "id": "0898265e-7edd-4093-9dfb-6e6b2f05c766", + "metadata": {}, + "source": [ + "## Compute intermediate results\n", + "\n", + "For inspection and debugging purposes, we can also compute intermediate results.\n", + "To avoid repeated computation (including costly loading of files), we can request multiple results at once, including the final result, if desired.\n", + "For example:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4bd5a7d5-581f-4d04-bf11-955e55646199", + "metadata": {}, + "outputs": [], + "source": [ + "intermediates = pipeline.compute(\n", + " (\n", + " DataWithScatteringCoordinates[SampleRun],\n", + " MaskedData[SampleRun],\n", + " )\n", + ")\n", + "\n", + "intermediates[DataWithScatteringCoordinates[SampleRun]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "185cb6ee-8e55-43d0-9505-ebe2761785ed", + "metadata": {}, + "outputs": [], + "source": [ + "intermediates[MaskedData[SampleRun]].bins.concat().hist(\n", + " two_theta=300, wavelength=300\n", + ").plot(norm=\"log\")" + ] + }, + { + "cell_type": "markdown", + "id": "7f866ad4-5a0d-4c98-b5ab-a2436f97074d", + "metadata": {}, + "source": [ + "## Grouping by scattering angle\n", + "\n", + "The above pipeline focuses the data by merging all instrument pixels to produce a 1d d-spacing curve.\n", + "If instead we want to group into $2\\theta$ bins, we can alter the pipeline parameters by adding some binning in $2\\theta$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc516311-1ab4-4eb7-9a0b-d83ff61e036b", + "metadata": {}, + "outputs": [], + "source": [ + "pipeline[TwoThetaBins] = sc.linspace(\n", + " dim=\"two_theta\", unit=\"rad\", start=0.8, stop=2.4, num=17\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b19ab368-e0c4-4a03-ac20-61a442e8826e", + "metadata": {}, + "outputs": [], + "source": [ + "grouped_dspacing = pipeline.compute(IofDspacingTwoTheta)\n", + "grouped_dspacing" ] }, { "cell_type": "code", "execution_count": null, - "id": "4b2c884c-776e-4bb0-b517-10f5ddbf8aba", + "id": "fc25ee9c-9395-481f-bd31-3ad139859686", "metadata": {}, "outputs": [], "source": [ - "result.plot()" + "angle = sc.midpoints(grouped_dspacing.coords[\"two_theta\"])\n", + "sc.plot(\n", + " {\n", + " f\"{angle[group].value:.3f} {angle[group].unit}\": grouped_dspacing[\n", + " \"two_theta\", group\n", + " ].hist()\n", + " for group in range(2, 6)\n", + " }\n", + ")" ] }, { "cell_type": "code", "execution_count": null, - "id": "27793efc-ecf2-4e1a-a73f-940ab652007c", + "id": "6e414419-ccf4-420b-899d-a1fa4a6ce4ca", "metadata": {}, "outputs": [], "source": [ - "result.coords['dspacing'] = sc.midpoints(result.coords['dspacing'])\n", - "scn.io.save_xye('dspacing.xye', result)" + "grouped_dspacing.hist().plot(norm=\"log\")" ] } ], @@ -156,8 +282,7 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.13" + "pygments_lexer": "ipython3" } }, "nbformat": 4, diff --git a/docs/user-guide/dream/dream-instrument-view.ipynb b/docs/user-guide/dream/dream-instrument-view.ipynb index 3a3765a2..e1bc9cff 100644 --- a/docs/user-guide/dream/dream-instrument-view.ipynb +++ b/docs/user-guide/dream/dream-instrument-view.ipynb @@ -43,19 +43,15 @@ "metadata": {}, "outputs": [], "source": [ - "dg = dream.io.load_geant4_csv(\n", - " dream.data.get_path('data_dream0_new_hkl_Si_pwd.csv.zip')\n", - ")\n", - "dg = dg['instrument'] # Extract the instrument data\n", + "dg = dream.io.load_geant4_csv(dream.data.get_path(\"data_dream0_new_hkl_Si_pwd.csv.zip\"))\n", + "dg = dg[\"instrument\"] # Extract the instrument data\n", "\n", "# Extract the events from nested data groups\n", - "dg = sc.DataGroup({\n", - " key: detector['events'] for key, detector in dg.items()\n", - "})\n", + "dg = sc.DataGroup({key: detector[\"events\"] for key, detector in dg.items()})\n", "\n", "# Construct the pixel positions from event positions\n", "for da in dg.values():\n", - " da.coords['position'] = da.bins.coords['position'].bins.mean()\n", + " da.coords[\"position\"] = da.bins.coords[\"position\"].bins.mean()\n", "\n", "dg" ] @@ -75,7 +71,7 @@ "outputs": [], "source": [ "# Only plot half of the pixels to reduce html docs size\n", - "dg = dg['counter', 0]" + "dg = dg[\"counter\", 0]" ] }, { @@ -108,7 +104,7 @@ }, "outputs": [], "source": [ - "tof_edges = sc.linspace('tof', 1.0e7, 1.0e8, 51, unit='ns', dtype=int)\n", + "tof_edges = sc.linspace(\"tof\", 1.0e7, 1.0e8, 51, unit=\"ns\", dtype=int)\n", "data = dg.hist(tof=tof_edges)" ] }, @@ -140,7 +136,7 @@ }, "outputs": [], "source": [ - "full_view = dream.instrument_view(data, dim='tof')\n", + "full_view = dream.instrument_view(data, dim=\"tof\")\n", "full_view" ] }, @@ -158,7 +154,7 @@ }, "outputs": [], "source": [ - "full_view[2].controls['tof']['slider'].value = 35" + "full_view[2].controls[\"tof\"][\"slider\"].value = 35" ] }, { @@ -192,7 +188,7 @@ }, "outputs": [], "source": [ - "mantle_view = dream.instrument_view(dg['mantle'].hist(wavelength=50), dim='wavelength')\n", + "mantle_view = dream.instrument_view(dg[\"mantle\"].hist(wavelength=50), dim=\"wavelength\")\n", "mantle_view" ] }, @@ -210,7 +206,7 @@ }, "outputs": [], "source": [ - "mantle_view[1].controls['wavelength']['slider'].value = 43" + "mantle_view[1].controls[\"wavelength\"][\"slider\"].value = 43" ] }, { @@ -231,7 +227,7 @@ "metadata": {}, "outputs": [], "source": [ - "dream.instrument_view(dg['endcap_backward']['module', 0].hist(tof=1))" + "dream.instrument_view(dg[\"endcap_backward\"][\"module\", 0].hist(tof=1))" ] } ], diff --git a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb index f7853190..c6c55efb 100644 --- a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb +++ b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb @@ -50,19 +50,18 @@ "outputs": [], "source": [ "params = {\n", + " NeXusDetectorName: \"powgen_detector\",\n", " # Input data\n", - " Filename[SampleRun]: 'PG3_4844_event.zip',\n", - " Filename[VanadiumRun]: 'PG3_4866_event.zip',\n", - " CalibrationFilename: 'PG3_FERNS_d4832_2011_08_24.zip',\n", - "\n", + " Filename[SampleRun]: powgen.data.powgen_tutorial_sample_file(),\n", + " Filename[VanadiumRun]: powgen.data.powgen_tutorial_vanadium_file(),\n", + " CalibrationFilename: powgen.data.powgen_tutorial_calibration_file(),\n", " # The upper bounds mode is not yet implemented.\n", " UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,\n", - "\n", - " # Crop data to this range in time-of-flight\n", - " ValidTofRange: sc.array(dims=['tof'], values=[0.0, 16666.67], unit='us'),\n", - "\n", " # Edges for binning in d-spacing\n", - " DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 200, unit='angstrom'),\n", + " DspacingBins: sc.linspace(\"dspacing\", 0.0, 2.3434, 201, unit=\"angstrom\"),\n", + " # Mask in time-of-flight to crop to valid range\n", + " TofMask: lambda x: (x < sc.scalar(0.0, unit=\"us\"))\n", + " | (x > sc.scalar(16666.67, unit=\"us\")),\n", "}" ] }, @@ -83,14 +82,8 @@ "metadata": {}, "outputs": [], "source": [ - "providers = [\n", - " *powder.providers,\n", - " *powgen.providers,\n", - "]\n", - "pipeline = sciline.Pipeline(\n", - " providers,\n", - " params=params,\n", - ")" + "providers = [*powder.providers, *powgen.providers]\n", + "pipeline = sciline.Pipeline(providers, params=params)" ] }, { @@ -102,27 +95,17 @@ "\n", "### Compute final result\n", "\n", - "We can get the graph for computing the final d-spacing histogram:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7d72a6ad-b43d-4e05-a70f-7f334ca06b7a", - "metadata": {}, - "outputs": [], - "source": [ - "pipeline.visualize(DspacingHistogram, graph_attr={'rankdir': 'LR'})" + "We can get the graph for computing the final intensity as a function of d-spacing:" ] }, { "cell_type": "code", "execution_count": null, - "id": "d09bc8d5-4251-483f-bb45-aa8cf67f69b2", + "id": "d3069da8-de3c-4cc3-bd14-ff8fcbf58e91", "metadata": {}, "outputs": [], "source": [ - "dspacing_histogram = pipeline.get(DspacingHistogram)" + "pipeline.visualize(IofDspacing, graph_attr={\"rankdir\": \"LR\"})" ] }, { @@ -140,7 +123,7 @@ "metadata": {}, "outputs": [], "source": [ - "result = dspacing_histogram.compute()\n", + "result = pipeline.compute(IofDspacing)\n", "result" ] }, @@ -151,7 +134,8 @@ "metadata": {}, "outputs": [], "source": [ - "result.plot()" + "dspacing_histogram = result.hist()\n", + "dspacing_histogram.plot()" ] }, { @@ -176,12 +160,13 @@ "metadata": {}, "outputs": [], "source": [ - "def save_xye(reduced_data: DspacingHistogram,\n", - " out_filename: OutFilename,\n", + "def save_xye(\n", + " reduced_data: IofDspacing,\n", + " out_filename: OutFilename,\n", ") -> None:\n", - " data = reduced_data.copy(deep=False)\n", - " data.coords['dspacing'] = sc.midpoints(data.coords['dspacing'])\n", - " scn.io.save_xye(out_filename, data, coord='dspacing')" + " data = reduced_data.hist()\n", + " data.coords[\"dspacing\"] = sc.midpoints(data.coords[\"dspacing\"])\n", + " scn.io.save_xye(out_filename, data, coord=\"dspacing\")" ] }, { @@ -200,7 +185,7 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline[OutFilename] = 'reduced.xye'" + "pipeline[OutFilename] = \"reduced.xye\"" ] }, { @@ -241,11 +226,14 @@ "metadata": {}, "outputs": [], "source": [ - "results = pipeline.compute((\n", - " RawDetectorData[SampleRun],\n", - " FilteredData[SampleRun],\n", - " FilteredData[VanadiumRun],\n", - "))" + "results = pipeline.compute(\n", + " (\n", + " ReducibleDetectorData[SampleRun],\n", + " MaskedData[SampleRun],\n", + " FilteredData[SampleRun],\n", + " FilteredData[VanadiumRun],\n", + " )\n", + ")" ] }, { @@ -255,7 +243,7 @@ "metadata": {}, "outputs": [], "source": [ - "results[RawDetectorData[SampleRun]]" + "results[ReducibleDetectorData[SampleRun]]" ] }, { @@ -265,7 +253,7 @@ "metadata": {}, "outputs": [], "source": [ - "scn.instrument_view(results[RawDetectorData[SampleRun]].hist())" + "results[MaskedData[SampleRun]].bins.concat().hist(wavelength=300).plot()" ] }, { @@ -276,8 +264,8 @@ "outputs": [], "source": [ "tof_data = sc.DataGroup(\n", - " sample=results[FilteredData[SampleRun]].bins.concat('spectrum'),\n", - " vanadium=results[FilteredData[VanadiumRun]].bins.concat('spectrum'),\n", + " sample=results[FilteredData[SampleRun]].bins.concat(),\n", + " vanadium=results[FilteredData[VanadiumRun]].bins.concat(),\n", ")\n", "tof_data.hist(tof=100).plot()" ] @@ -290,60 +278,45 @@ "## Group by scattering angle\n", "\n", "The above pipeline focuses the data by merging all instrument pixels to produce a 1d d-spacing curve.\n", - "If instead we want to group into $2\\theta$ bins, we can alter the pipeline by replacing the focussing step:" + "If instead we want to group into $2\\theta$ bins, we can alter the pipeline parameters by adding some binning in $2\\theta$:" ] }, { "cell_type": "code", "execution_count": null, - "id": "d24c8bf8-ad59-4211-ae6a-3ee29b0556a3", + "id": "a4b68853-a70b-42d6-a8cc-58c77e83eaec", "metadata": {}, "outputs": [], "source": [ - "from ess.powder.grouping import group_by_two_theta, merge_all_pixels\n", - "\n", - "grouping_providers = list(providers)\n", - "grouping_providers.remove(merge_all_pixels)\n", - "grouping_providers = (*grouping_providers, group_by_two_theta)" + "pipeline[TwoThetaBins] = sc.linspace(\n", + " dim=\"two_theta\", unit=\"deg\", start=25.0, stop=90.0, num=17\n", + ").to(unit=\"rad\")" ] }, { "cell_type": "markdown", - "id": "3e83ed6c-46fc-4396-8331-3654993def94", - "metadata": {}, - "source": [ - "We also need to specify the grouping with a new parameter:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a4b68853-a70b-42d6-a8cc-58c77e83eaec", + "id": "c8377091-43b5-467a-8f86-b58a96a3dc89", "metadata": {}, - "outputs": [], "source": [ - "params[TwoThetaBins] = sc.linspace(dim='two_theta', unit='deg', start=25.0, stop=90.0, num=16)" + "We then have to request a final result that depends on both d-spacing and $2\\theta$:" ] }, { "cell_type": "code", "execution_count": null, - "id": "ea65b6bd-b9fa-496d-bfa2-f459e33f75fc", + "id": "e88eb36e-9377-4415-8123-a8916fd14793", "metadata": {}, "outputs": [], "source": [ - "grouping_pipeline = sciline.Pipeline(\n", - " grouping_providers,\n", - " params=params\n", - ")" + "pipeline.visualize(IofDspacingTwoTheta, graph_attr={\"rankdir\": \"LR\"})" ] }, { "cell_type": "markdown", - "id": "59ae5b45-8e11-4de1-bef5-c4b58df75804", + "id": "2a8a0853-2829-49ab-b343-31deca3de31c", "metadata": {}, "source": [ - "Inspect the graph to check that the new provider has been inserted:" + "Compute and plot the result:" ] }, { @@ -353,41 +326,44 @@ "metadata": {}, "outputs": [], "source": [ - "grouped_dspacing = grouping_pipeline.get(DspacingHistogram)\n", - "grouped_dspacing.visualize(graph_attr={'rankdir': 'LR'})" + "grouped_dspacing = pipeline.compute(IofDspacingTwoTheta)\n", + "grouped_dspacing" ] }, { - "cell_type": "markdown", - "id": "2a8a0853-2829-49ab-b343-31deca3de31c", + "cell_type": "code", + "execution_count": null, + "id": "7ca04947-b019-4814-9467-4d0519d8d384", "metadata": {}, + "outputs": [], "source": [ - "Compute and plot the result:" + "angle = sc.midpoints(grouped_dspacing.coords[\"two_theta\"])\n", + "sc.plot(\n", + " {\n", + " f\"{angle[group].value:.3f} {angle[group].unit}\": grouped_dspacing[\n", + " \"two_theta\", group\n", + " ].hist()\n", + " for group in range(2, 6)\n", + " }\n", + ")" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "bc73276a-dadb-4c21-9954-e69bebc2ff90", + "cell_type": "markdown", + "id": "3de8061a-7c28-48e2-b569-a08fa09f1295", "metadata": {}, - "outputs": [], "source": [ - "grouped_result = grouped_dspacing.compute()\n", - "grouped_result" + "Or we can view it as a 2D plot, which should display powder peaks as vertical bright lines:" ] }, { "cell_type": "code", "execution_count": null, - "id": "7ca04947-b019-4814-9467-4d0519d8d384", + "id": "4c9ee765-0699-4ca2-9bd8-f8dea27a261c", "metadata": {}, "outputs": [], "source": [ - "angle = sc.midpoints(grouped_result.coords['two_theta'])\n", - "sc.plot({\n", - " f'{angle[group].value:.3f} {angle[group].unit}': grouped_result['two_theta', group]\n", - " for group in range(2, 6)\n", - "})" + "grouped_dspacing.hist().plot()" ] } ], @@ -406,8 +382,7 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.13" + "pygments_lexer": "ipython3" } }, "nbformat": 4, diff --git a/docs/user-guide/sns-instruments/preprocess_files.ipynb b/docs/user-guide/sns-instruments/preprocess_files.ipynb index 4cd44b42..e1a8971c 100644 --- a/docs/user-guide/sns-instruments/preprocess_files.ipynb +++ b/docs/user-guide/sns-instruments/preprocess_files.ipynb @@ -44,7 +44,8 @@ " powgen.data.mantid_sample_file(),\n", " advanced_geometry=True,\n", " load_pulse_times=False,\n", - " mantid_args={\"LoadMonitors\": False})" + " mantid_args={\"LoadMonitors\": False},\n", + ")" ] }, { @@ -64,13 +65,15 @@ "metadata": {}, "outputs": [], "source": [ - "sample_data = sample['data']\n", - "sample_data.coords['gd_prtn_chrg'] = sample['gd_prtn_chrg']\n", - "sample_data.coords.set_aligned('gd_prtn_chrg', False)\n", - "sample_dg = sc.DataGroup({\n", - " 'data': sample_data,\n", - " 'detector_info': sample_data.coords.pop('detector_info').value\n", - "})" + "sample_data = sample[\"data\"]\n", + "sample_data.coords[\"gd_prtn_chrg\"] = sample[\"gd_prtn_chrg\"]\n", + "sample_data.coords.set_aligned(\"gd_prtn_chrg\", False)\n", + "sample_dg = sc.DataGroup(\n", + " {\n", + " \"data\": sample_data,\n", + " \"detector_info\": sample_data.coords.pop(\"detector_info\").value,\n", + " }\n", + ")" ] }, { @@ -90,7 +93,7 @@ "metadata": {}, "outputs": [], "source": [ - "sample_dg.save_hdf5('PG3_4844_event.h5')" + "sample_dg.save_hdf5(\"PG3_4844_event.h5\")" ] }, { @@ -100,7 +103,7 @@ "metadata": {}, "outputs": [], "source": [ - "sc.io.load_hdf5('PG3_4844_event.h5')" + "sc.io.load_hdf5(\"PG3_4844_event.h5\")" ] }, { @@ -122,7 +125,8 @@ " powgen.data.mantid_vanadium_file(),\n", " advanced_geometry=False,\n", " load_pulse_times=True,\n", - " mantid_args={\"LoadMonitors\": False})" + " mantid_args={\"LoadMonitors\": False},\n", + ")" ] }, { @@ -132,13 +136,15 @@ "metadata": {}, "outputs": [], "source": [ - "vana_data = vana['data']\n", - "vana_data.coords['gd_prtn_chrg'] = vana['gd_prtn_chrg']\n", - "vana_data.coords.set_aligned('gd_prtn_chrg', False)\n", - "vana_dg = sc.DataGroup({\n", - " 'data': vana_data,\n", - " 'proton_charge': vana['proton_charge'].rename(time='pulse_time')\n", - "})" + "vana_data = vana[\"data\"]\n", + "vana_data.coords[\"gd_prtn_chrg\"] = vana[\"gd_prtn_chrg\"]\n", + "vana_data.coords.set_aligned(\"gd_prtn_chrg\", False)\n", + "vana_dg = sc.DataGroup(\n", + " {\n", + " \"data\": vana_data,\n", + " \"proton_charge\": vana[\"proton_charge\"].rename(time=\"pulse_time\"),\n", + " }\n", + ")" ] }, { @@ -158,7 +164,7 @@ "metadata": {}, "outputs": [], "source": [ - "vana_dg['data'].bins.constituents['data']" + "vana_dg[\"data\"].bins.constituents[\"data\"]" ] }, { @@ -168,7 +174,7 @@ "metadata": {}, "outputs": [], "source": [ - "vana_dg.save_hdf5('PG3_4866_event.h5')" + "vana_dg.save_hdf5(\"PG3_4866_event.h5\")" ] }, { @@ -178,7 +184,7 @@ "metadata": {}, "outputs": [], "source": [ - "sc.io.load_hdf5('PG3_4866_event.h5')" + "sc.io.load_hdf5(\"PG3_4866_event.h5\")" ] }, { @@ -198,7 +204,7 @@ "source": [ "cal = load_calibration(\n", " powgen.data.mantid_calibration_file(),\n", - " instrument_filename='POWGEN_Definition_2011-02-25.xml',\n", + " instrument_filename=\"POWGEN_Definition_2011-02-25.xml\",\n", ")" ] }, @@ -219,7 +225,7 @@ "metadata": {}, "outputs": [], "source": [ - "cal.save_hdf5('PG3_FERNS_d4832_2011_08_24.h5')" + "cal.save_hdf5(\"PG3_FERNS_d4832_2011_08_24.h5\")" ] }, { @@ -229,7 +235,7 @@ "metadata": {}, "outputs": [], "source": [ - "sc.io.load_hdf5('PG3_FERNS_d4832_2011_08_24.h5')" + "sc.io.load_hdf5(\"PG3_FERNS_d4832_2011_08_24.h5\")" ] } ], diff --git a/src/ess/dream/__init__.py b/src/ess/dream/__init__.py index bb50dc01..020bf577 100644 --- a/src/ess/dream/__init__.py +++ b/src/ess/dream/__init__.py @@ -7,9 +7,9 @@ import importlib.metadata -from . import beamline, data +from . import data from .instrument_view import instrument_view -from .io import fold_nexus_detectors, load_geant4_csv, load_nexus +from .io import load_geant4_csv, nexus try: __version__ = importlib.metadata.version(__package__ or __name__) @@ -18,13 +18,12 @@ del importlib -providers = (*beamline.providers,) +providers = (*nexus.providers,) __all__ = [ 'data', 'beamline', - 'fold_nexus_detectors', 'instrument_view', 'load_geant4_csv', - 'load_nexus', + 'nexus', ] diff --git a/src/ess/dream/beamline.py b/src/ess/dream/beamline.py deleted file mode 100644 index 2b0a5bb1..00000000 --- a/src/ess/dream/beamline.py +++ /dev/null @@ -1,50 +0,0 @@ -# SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) - -"""Beamline tools for DREAM.""" - -from ess.powder.types import DetectorDimensions, RawDetectorData, SampleRun - - -def dream_detector_dimensions(data: RawDetectorData[SampleRun]) -> DetectorDimensions: - """Logical dimensions used by a DREAM detector. - - The dimensions differ between simulated data loaded from GEANT4 CSV files - and measured data loaded from NeXus files. - The dimensions returned by this function match the dimensions found - in the ``data`` argument. - - Parameters - ---------- - data: - Dimensions are deduced based on this data. - - Returns - ------- - : - Logical dimensions used by the given DREAM detector. - """ - geant4_dims = { - 'module', - 'segment', - 'counter', - 'wire', - 'strip', - 'sector', - } - nexus_dims = { - 'wire', - 'mounting_unit', - 'cassette', - 'counter', - 'strip', - 'sector', - 'sumo_cass_ctr', - 'other', - } - dims = (geant4_dims | nexus_dims) & set(data.dims) - return DetectorDimensions(tuple(dims)) - - -providers = (dream_detector_dimensions,) -"""Sciline providers for DREAM detector handling.""" diff --git a/src/ess/dream/data.py b/src/ess/dream/data.py index fe1b694d..00ae7e88 100644 --- a/src/ess/dream/data.py +++ b/src/ess/dream/data.py @@ -2,30 +2,28 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) """Data for tests and documentation with DREAM.""" -from pathlib import Path +_version = "1" -_version = '1' - -__all__ = ['get_path'] +__all__ = ["get_path"] def _make_pooch(): import pooch return pooch.create( - path=pooch.os_cache('ess/dream'), - env='ESS_DATA_DIR', - base_url='https://public.esss.dk/groups/scipp/ess/dream/{version}/', + path=pooch.os_cache("ess/dream"), + env="ESS_DATA_DIR", + base_url="https://public.esss.dk/groups/scipp/ess/dream/{version}/", version=_version, registry={ - 'data_dream_with_sectors.csv.zip': 'md5:52ae6eb3705e5e54306a001bc0ae85d8', - 'data_dream0_new_hkl_Si_pwd.csv.zip': 'md5:d0ae518dc1b943bb817b3d18c354ed01', # noqa: E501 - 'DREAM_nexus_sorted-2023-12-07.nxs': 'md5:22824e14f6eb950d24a720b2a0e2cb66', - 'DREAM_simple_pwd_workflow/data_dream_diamond_vana_container_sample_union.csv.zip': 'md5:33302d0506b36aab74003b8aed4664cc', # noqa: E501 - 'DREAM_simple_pwd_workflow/data_dream_diamond_vana_container_sample_union_run2.csv.zip': 'md5:c7758682f978d162dcb91e47c79abb83', # noqa: E501 - 'DREAM_simple_pwd_workflow/data_dream_vana_container_sample_union.csv.zip': 'md5:1e22917b2bb68b5cacfb506b72700a4d', # noqa: E501 - 'DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip': 'md5:e5addfc06768140c76533946433fa2ec', # noqa: E501 - 'DREAM_simple_pwd_workflow/data_dream_vanadium_inc_coh.csv.zip': 'md5:39d1a44e248b12966b26f7c2f6c602a2', # noqa: E501 + "data_dream_with_sectors.csv.zip": "md5:52ae6eb3705e5e54306a001bc0ae85d8", + "data_dream0_new_hkl_Si_pwd.csv.zip": "md5:d0ae518dc1b943bb817b3d18c354ed01", # noqa: E501 + "DREAM_nexus_sorted-2023-12-07.nxs": "md5:22824e14f6eb950d24a720b2a0e2cb66", + "DREAM_simple_pwd_workflow/data_dream_diamond_vana_container_sample_union.csv.zip": "md5:33302d0506b36aab74003b8aed4664cc", # noqa: E501 + "DREAM_simple_pwd_workflow/data_dream_diamond_vana_container_sample_union_run2.csv.zip": "md5:c7758682f978d162dcb91e47c79abb83", # noqa: E501 + "DREAM_simple_pwd_workflow/data_dream_vana_container_sample_union.csv.zip": "md5:1e22917b2bb68b5cacfb506b72700a4d", # noqa: E501 + "DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip": "md5:e5addfc06768140c76533946433fa2ec", # noqa: E501 + "DREAM_simple_pwd_workflow/data_dream_vanadium_inc_coh.csv.zip": "md5:39d1a44e248b12966b26f7c2f6c602a2", # noqa: E501 }, ) @@ -33,7 +31,7 @@ def _make_pooch(): _pooch = _make_pooch() -def get_path(name: str, unzip: bool = False) -> Path: +def get_path(name: str, unzip: bool = False) -> str: """ Return the path to a data file bundled with ess.dream. @@ -42,10 +40,10 @@ def get_path(name: str, unzip: bool = False) -> Path: """ import pooch - return Path(_pooch.fetch(name, processor=pooch.Unzip() if unzip else None)) + return _pooch.fetch(name, processor=pooch.Unzip() if unzip else None) -def simulated_diamond_sample() -> Path: +def simulated_diamond_sample() -> str: """Path to a GEANT4 CSV file for a diamond sample. **Sample**: @@ -70,12 +68,12 @@ def simulated_diamond_sample() -> Path: absorption of 36.73 1/m """ return get_path( - 'DREAM_simple_pwd_workflow/' - 'data_dream_diamond_vana_container_sample_union.csv.zip' + "DREAM_simple_pwd_workflow/" + "data_dream_diamond_vana_container_sample_union.csv.zip" ) -def simulated_vanadium_sample() -> Path: +def simulated_vanadium_sample() -> str: """Path to a GEANT4 CSV file for a vanadium sample. Contains both coherent and incoherent scattering. @@ -90,10 +88,10 @@ def simulated_vanadium_sample() -> Path: unit cell volume=27.66 angstrom^3 absorption of 36.73 1/m """ - return get_path('DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip') + return get_path("DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip") -def simulated_vanadium_sample_incoherent() -> Path: +def simulated_vanadium_sample_incoherent() -> str: """Path to a GEANT4 CSV file for a vanadium sample with only incoherent scattering. **Sample**: @@ -102,10 +100,10 @@ def simulated_vanadium_sample_incoherent() -> Path: - vertical dimension of sample (along y)=0.01 m - packing factor=1 """ - return get_path('DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip') + return get_path("DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip") -def simulated_empty_can() -> Path: +def simulated_empty_can() -> str: """Path to a GEANT4 CSV file for an empty can measurement. **Container**: @@ -119,5 +117,5 @@ def simulated_empty_can() -> Path: absorption of 36.73 1/m """ return get_path( - 'DREAM_simple_pwd_workflow/data_dream_vana_container_sample_union.csv.zip' + "DREAM_simple_pwd_workflow/data_dream_vana_container_sample_union.csv.zip" ) diff --git a/src/ess/dream/instrument_view.py b/src/ess/dream/instrument_view.py index d1dab39f..0c7c529e 100644 --- a/src/ess/dream/instrument_view.py +++ b/src/ess/dream/instrument_view.py @@ -50,7 +50,7 @@ def instrument_view( def _to_data_group(data: Union[sc.DataArray, sc.DataGroup, dict]) -> sc.DataGroup: if isinstance(data, sc.DataArray): - data = sc.DataGroup({data.name or 'data': data}) + data = sc.DataGroup({data.name or "data": data}) elif isinstance(data, dict): data = sc.DataGroup(data) return data @@ -61,8 +61,8 @@ def _pre_process(da: sc.DataArray, dim: str) -> sc.DataArray: dims = list(da.dims) if dim is not None: dims.remove(dim) - out = da.flatten(dims=dims, to='pixel') - sel = sc.isfinite(out.coords['position']) + out = da.flatten(dims=dims, to="pixel") + sel = sc.isfinite(out.coords["position"]) return out[sel] @@ -87,7 +87,7 @@ def __init__( if dim is not None: self.slider = SliceWidget(next(iter(self.data.values())), dims=[dim]) - self.slider.controls[dim]['slider'].layout = {'width': '600px'} + self.slider.controls[dim]["slider"].layout = {"width": "600px"} self.slider_node = pp.widget_node(self.slider) self.slice_nodes = { key: slice_dims(n, self.slider_node) @@ -101,8 +101,8 @@ def __init__( self.fig = pp.scatter3d( to_scatter, - pos='position', - pixel_size=1.0 * sc.Unit('cm') if pixel_size is None else pixel_size, + pos="position", + pixel_size=1.0 * sc.Unit("cm") if pixel_size is None else pixel_size, **kwargs, ) @@ -137,7 +137,7 @@ def _add_module_control(self): ) for key, ch in self.checkboxes.items(): ch.key = key - ch.observe(self._check_visibility, names='value') + ch.observe(self._check_visibility, names="value") self.children.insert(0, self.modules_widget) def _check_visibility(self, _): diff --git a/src/ess/dream/io/__init__.py b/src/ess/dream/io/__init__.py index 605738ed..2224a078 100644 --- a/src/ess/dream/io/__init__.py +++ b/src/ess/dream/io/__init__.py @@ -3,7 +3,7 @@ """Input/output for DREAM.""" +from . import nexus from .geant4 import load_geant4_csv -from .nexus import fold_nexus_detectors, load_nexus -__all__ = ["fold_nexus_detectors", "load_geant4_csv", "load_nexus"] +__all__ = ["nexus", "load_geant4_csv"] diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index eb602a2c..09bb5a4b 100644 --- a/src/ess/dream/io/geant4.py +++ b/src/ess/dream/io/geant4.py @@ -1,19 +1,24 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) -import os -from io import BytesIO, StringIO -from typing import Dict, Optional, Union +from typing import Dict, Optional import numpy as np import sciline import scipp as sc from ess.powder.types import ( - DetectorName, - FilePath, + Filename, + NeXusDetectorDimensions, + NeXusDetectorName, RawDetector, RawDetectorData, + RawSample, + RawSource, + ReducibleDetectorData, RunType, + SamplePosition, + SampleRun, + SourcePosition, ) from ess.reduce.nexus import extract_detector_data @@ -27,9 +32,7 @@ class AllRawDetectors(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): """Raw data for all detectors.""" -def load_geant4_csv( - file_path: Union[FilePath[RunType], str, StringIO, BytesIO], -) -> AllRawDetectors[RunType]: +def load_geant4_csv(file_path: Filename[RunType]) -> AllRawDetectors[RunType]: """Load a GEANT4 CSV file for DREAM. Parameters @@ -54,15 +57,15 @@ def load_geant4_csv( detectors = _group(detectors) return AllRawDetectors[RunType]( - sc.DataGroup({'instrument': sc.DataGroup(detectors)}) + sc.DataGroup({"instrument": sc.DataGroup(detectors)}) ) def extract_geant4_detector( - detectors: AllRawDetectors[RunType], detector_name: DetectorName + detectors: AllRawDetectors[RunType], detector_name: NeXusDetectorName ) -> RawDetector[RunType]: """Extract a single detector from a loaded GEANT4 simulation.""" - return RawDetector[RunType](detectors['instrument'][detector_name.name]) + return RawDetector[RunType](detectors["instrument"][detector_name]) def extract_geant4_detector_data( @@ -72,42 +75,41 @@ def extract_geant4_detector_data( return RawDetectorData[RunType](extract_detector_data(detector)) -def _load_raw_events( - file_path: Union[str, os.PathLike, StringIO, BytesIO], -) -> sc.DataArray: +def _load_raw_events(file_path: str) -> sc.DataArray: table = sc.io.load_csv( - file_path, sep='\t', header_parser='bracket', data_columns=[] + file_path, sep="\t", header_parser="bracket", data_columns=[] ) - table = table.rename_dims(row='event') + table = table.rename_dims(row="event") return sc.DataArray( - sc.ones(sizes=table.sizes, with_variances=True, unit='counts'), + sc.ones(sizes=table.sizes, with_variances=True, unit="counts"), coords=table.coords, ) def _adjust_coords(da: sc.DataArray) -> None: - da.coords['wavelength'] = da.coords.pop('lambda') - da.coords['position'] = sc.spatial.as_vectors( - da.coords.pop('x_pos'), da.coords.pop('y_pos'), da.coords.pop('z_pos') + da.coords["wavelength"] = da.coords.pop("lambda") + da.coords["wavelength"].unit = "angstrom" + da.coords["position"] = sc.spatial.as_vectors( + da.coords.pop("x_pos"), da.coords.pop("y_pos"), da.coords.pop("z_pos") ) def _group(detectors: Dict[str, sc.DataArray]) -> Dict[str, sc.DataGroup]: - elements = ('module', 'segment', 'counter', 'wire', 'strip') + elements = ("module", "segment", "counter", "wire", "strip") def group(key: str, da: sc.DataArray) -> sc.DataArray: - if key in ['high_resolution', 'sans']: + if key in ["high_resolution", "sans"]: # Only the HR and SANS detectors have sectors. - return da.group('sector', *elements) + return da.group("sector", *elements) res = da.group(*elements) - res.bins.coords.pop('sector', None) + res.bins.coords.pop("sector", None) return res return {key: sc.DataGroup(events=group(key, da)) for key, da in detectors.items()} def _split_detectors( - data: sc.DataArray, detector_id_name: str = 'det ID' + data: sc.DataArray, detector_id_name: str = "det ID" ) -> Dict[str, sc.DataArray]: groups = data.group( sc.concat( @@ -124,15 +126,15 @@ def _split_detectors( if ( mantle := _extract_detector(groups, detector_id_name, MANTLE_DETECTOR_ID) ) is not None: - detectors['mantle'] = mantle.copy() + detectors["mantle"] = mantle.copy() if ( high_res := _extract_detector(groups, detector_id_name, HIGH_RES_DETECTOR_ID) ) is not None: - detectors['high_resolution'] = high_res.copy() + detectors["high_resolution"] = high_res.copy() if ( sans := _extract_detector(groups, detector_id_name, SANS_DETECTOR_ID) ) is not None: - detectors['sans'] = sans.copy() + detectors["sans"] = sans.copy() endcaps_list = [ det @@ -143,13 +145,13 @@ def _split_detectors( endcaps = sc.concat(endcaps_list, data.dim) endcaps = endcaps.bin( z_pos=sc.array( - dims=['z_pos'], + dims=["z_pos"], values=[-np.inf, 0.0, np.inf], - unit=endcaps.coords['z_pos'].unit, + unit=endcaps.coords["z_pos"].unit, ) ) - detectors['endcap_backward'] = endcaps[0].bins.concat().value.copy() - detectors['endcap_forward'] = endcaps[1].bins.concat().value.copy() + detectors["endcap_backward"] = endcaps[0].bins.concat().value.copy() + detectors["endcap_forward"] = endcaps[1].bins.concat().value.copy() return detectors @@ -163,5 +165,44 @@ def _extract_detector( return events -providers = (extract_geant4_detector, extract_geant4_detector_data, load_geant4_csv) +def get_source_position( + raw_source: RawSource[RunType], +) -> SourcePosition[RunType]: + return SourcePosition[RunType](raw_source["position"]) + + +def get_sample_position( + raw_sample: RawSample[RunType], +) -> SamplePosition[RunType]: + return SamplePosition[RunType](raw_sample["position"]) + + +def patch_detector_data( + detector_data: RawDetectorData[RunType], + source_position: SourcePosition[RunType], + sample_position: SamplePosition[RunType], +) -> ReducibleDetectorData[RunType]: + out = detector_data.copy(deep=False) + out.coords["source_position"] = source_position + out.coords["sample_position"] = sample_position + return ReducibleDetectorData[RunType](out) + + +def geant4_detector_dimensions( + data: RawDetectorData[SampleRun], +) -> NeXusDetectorDimensions[NeXusDetectorName]: + # For geant4 data, we group by detector identifier, so the data already has + # logical dimensions, so we simply return the dimensions of the detector. + return NeXusDetectorDimensions[NeXusDetectorName](data.sizes) + + +providers = ( + extract_geant4_detector, + extract_geant4_detector_data, + load_geant4_csv, + get_sample_position, + get_source_position, + patch_detector_data, + geant4_detector_dimensions, +) """Geant4-providers for Sciline pipelines.""" diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index dd485e61..3de86b88 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -1,135 +1,144 @@ # SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +# Copyright (c) 2024 Scipp contributors (https://github.com/scipp) -"""NeXus input/output for DREAM.""" +"""NeXus input/output for DREAM. -import os -from typing import Any, Dict, Union +Notes on the detector dimensions (2024-05-22): + +See https://confluence.esss.lu.se/pages/viewpage.action?pageId=462000005 +and the ICD DREAM interface specification for details. + +- The high-resolution and SANS detectors have a very odd numbering scheme. + The scheme attempts to follows some sort of physical ordering in space (x,y,z), + but it is not possible to reshape the data into all the logical dimensions. +""" import scipp as sc -import scippnexus as snx +from ess.powder.types import ( + Filename, + LoadedNeXusDetector, + NeXusDetectorName, + RawDetectorData, + RawSample, + RawSource, + ReducibleDetectorData, + RunType, + SamplePosition, + SourcePosition, +) +from ess.reduce import nexus + +DETECTOR_BANK_SIZES = { + "endcap_backward_detector": { + "strip": 16, + "wire": 16, + "module": 11, + "segment": 28, + "counter": 2, + }, + "endcap_forward_detector": { + "strip": 16, + "wire": 16, + "module": 5, + "segment": 28, + "counter": 2, + }, + "mantle_detector": { + "wire": 32, + "module": 5, + "segment": 6, + "strip": 256, + "counter": 2, + }, + "high_resolution_detector": { + "strip": 32, + "other": -1, + }, + "sans_detector": lambda x: x.fold( + dim="detector_number", + sizes={ + "strip": 32, + "other": -1, + }, + ), +} + + +def load_nexus_sample(file_path: Filename[RunType]) -> RawSample[RunType]: + return RawSample[RunType](nexus.load_sample(file_path)) -def load_nexus( - path: Union[str, os.PathLike], - *, - load_pixel_shape: bool = False, - entry: str = 'entry', - fold_detectors: bool = True, -) -> sc.DataGroup: +def dummy_load_sample(file_path: Filename[RunType]) -> RawSample[RunType]: """ - Load an unprocessed DREAM NeXus file. - - See https://confluence.esss.lu.se/pages/viewpage.action?pageId=462000005 - and the ICD DREAM interface specification for details. - - Notes (2023-12-06): - - - Mounting-unit, cassette, and counter roughly correspond to the azimuthal angle - in the mantle detector. However, counter is reversed in the current files. This - may also be the case in the other detectors. - - The endcap detectors have a irregular structure that cannot be fully folded. - This is not a problem but note again that the counter may be reversed. It is - currently not clear if this is a bug in the file. - - The high-resolution detector has a very odd numbering scheme. The SANS detector - is using the same, but is not populated in the current files. The scheme - attempts to follows some sort of physical ordering in space (x,y,z), but it - looks partially messed up. - - Parameters - ---------- - path: - Path to the NeXus file. - load_pixel_shape: - If True, load the pixel shape from the file's NXoff_geometry group. This is - often unused by would slow down file loading. Default is False. - entry: - Name of the entry to load. Default is "entry". - fold_detectors: - If True, fold the detector data to (partially) mimic the logical detector - structure. Default is True. - - Returns - ------- - : - A data group with the loaded file contents. + In test files there is not always a sample, so we need a dummy. """ - definitions = snx.base_definitions() - if not load_pixel_shape: - definitions["NXdetector"] = FilteredDetector + return RawSample[RunType]( + sc.DataGroup({'position': sc.vector(value=[0, 0, 0], unit='m')}) + ) - with snx.File(path, definitions=definitions) as f: - dg = f[entry][()] - dg = snx.compute_positions(dg) - return fold_nexus_detectors(dg) if fold_detectors else dg +def load_nexus_source(file_path: Filename[RunType]) -> RawSource[RunType]: + return RawSource[RunType](nexus.load_source(file_path)) -def fold_nexus_detectors(dg: sc.DataGroup) -> sc.DataGroup: - """ - Fold the detector data in a DREAM NeXus file. - The detector banks in the returned data group will have a multi-dimensional shape, - following the logical structure as far as possible. Note that the full structure - cannot be folded, as some dimensions are irregular. - """ - dg = dg.copy() - dg['instrument'] = dg['instrument'].copy() - instrument = dg['instrument'] - mantle = instrument['mantle_detector'] - mantle['mantle_event_data'] = mantle['mantle_event_data'].fold( - dim='detector_number', - sizes={ - 'wire': 32, - 'mounting_unit': 5, - 'cassette': 6, - 'counter': 2, - 'strip': 256, - }, - ) - for direction in ('backward', 'forward'): - endcap = instrument[f'endcap_{direction}_detector'] - endcap[f'endcap_{direction}_event_data'] = endcap[ - f'endcap_{direction}_event_data' - ].fold( - dim='detector_number', - sizes={ - 'strip': 16, - 'wire': 16, - 'sector': 5 if direction == 'forward' else 11, - 'sumo_cass_ctr': -1, # sumo*cassette*counter, irregular, cannot fold - }, - ) - high_resolution = instrument['high_resolution_detector'] - high_resolution['high_resolution_event_data'] = high_resolution[ - 'high_resolution_event_data' - ].fold( - dim='detector_number', - sizes={ - 'strip': 32, - 'other': -1, # some random order that is hard to follow - }, - ) - sans = instrument['sans_detector'] - sans['sans_event_data'] = sans['sans_event_data'].fold( - dim='detector_number', - sizes={ - 'strip': 32, - 'other': -1, # some random order that is hard to follow - }, - ) - return dg +def load_nexus_detector( + file_path: Filename[RunType], detector_name: NeXusDetectorName +) -> LoadedNeXusDetector[RunType]: + out = nexus.load_detector(file_path=file_path, detector_name=detector_name) + out.pop("pixel_shape", None) + return LoadedNeXusDetector[RunType](out) -def _skip(name: str, obj: Union[snx.Field, snx.Group]) -> bool: - skip_classes = (snx.NXoff_geometry,) - return isinstance(obj, snx.Group) and (obj.nx_class in skip_classes) +def get_source_position( + raw_source: RawSource[RunType], +) -> SourcePosition[RunType]: + return SourcePosition[RunType](raw_source["position"]) -class FilteredDetector(snx.NXdetector): - def __init__( - self, attrs: Dict[str, Any], children: Dict[str, Union[snx.Field, snx.Group]] - ): - children = { - name: child for name, child in children.items() if not _skip(name, child) - } - super().__init__(attrs=attrs, children=children) +def get_sample_position( + raw_sample: RawSample[RunType], +) -> SamplePosition[RunType]: + return SamplePosition[RunType](raw_sample["position"]) + + +def get_detector_data( + detector: LoadedNeXusDetector[RunType], + detector_name: NeXusDetectorName, +) -> RawDetectorData[RunType]: + da = nexus.extract_detector_data(detector) + if detector_name in DETECTOR_BANK_SIZES: + da = da.fold(dim="detector_number", sizes=DETECTOR_BANK_SIZES[detector_name]) + return RawDetectorData[RunType](da) + + +def patch_detector_data( + detector_data: RawDetectorData[RunType], + source_position: SourcePosition[RunType], + sample_position: SamplePosition[RunType], +) -> ReducibleDetectorData[RunType]: + """ + Patch a detector data object with source and sample positions. + Also adds variances to the event data if they are missing. + """ + out = detector_data.copy(deep=False) + if out.bins is not None: + content = out.bins.constituents["data"] + if content.variances is None: + content.variances = content.values + out.coords["sample_position"] = sample_position + out.coords["source_position"] = source_position + return ReducibleDetectorData[RunType](out) + + +providers = ( + load_nexus_sample, + load_nexus_source, + load_nexus_detector, + get_source_position, + get_sample_position, + get_detector_data, + patch_detector_data, +) +""" +Providers for loading and processing DREAM NeXus data. +""" diff --git a/src/ess/powder/__init__.py b/src/ess/powder/__init__.py index e2c69c75..240592fc 100644 --- a/src/ess/powder/__init__.py +++ b/src/ess/powder/__init__.py @@ -1,13 +1,20 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) -# flake8: noqa: F401 """ Components for powder diffraction experiments. """ import importlib.metadata -from . import conversion, correction, filtering, grouping, smoothing, uncertainty +from . import ( + conversion, + correction, + filtering, + grouping, + masking, + smoothing, + uncertainty, +) try: __version__ = importlib.metadata.version(__package__ or __name__) @@ -17,9 +24,21 @@ del importlib providers = ( - *conversion.providers, + *conversion.providers_with_positions, *correction.providers, *filtering.providers, *grouping.providers, + *masking.providers, ) """Sciline providers for powder diffraction.""" + +__all__ = [ + "conversion", + "correction", + "filtering", + "grouping", + "masking", + "providers", + "smoothing", + "uncertainty", +] diff --git a/src/ess/powder/conversion.py b/src/ess/powder/conversion.py index 4ab63ee3..982ecf59 100644 --- a/src/ess/powder/conversion.py +++ b/src/ess/powder/conversion.py @@ -13,10 +13,10 @@ from .logging import get_logger from .types import ( CalibrationData, + DataWithScatteringCoordinates, DspacingData, + ElasticCoordTransformGraph, NormalizedByProtonCharge, - RawSample, - RawSource, RunType, ) @@ -130,102 +130,83 @@ def to_dspacing_with_calibration( out = _restore_tof_if_in_wavelength(out) graph = { - 'dspacing': _dspacing_from_diff_calibration, + "dspacing": _dspacing_from_diff_calibration, } # `_dspacing_from_diff_calibration` does not need positions but conceptually, # the conversion maps from positions to d-spacing. # The mechanism with `_tag_positions_consumed` is meant to ensure that, # if positions are present, they are consumed (mad unaligned or dropped) # by the coordinate transform similarly to `to_dspacing_with_positions`. - if 'position' in out.coords or ( - out.bins is not None and 'position' in out.bins.coords + if "position" in out.coords or ( + out.bins is not None and "position" in out.bins.coords ): - graph['_tag_positions_consumed'] = _consume_positions + graph["_tag_positions_consumed"] = _consume_positions else: - graph['_tag_positions_consumed'] = lambda: sc.scalar(0) - out = out.transform_coords('dspacing', graph=graph, keep_intermediate=False) - out.coords.pop('_tag_positions_consumed', None) + graph["_tag_positions_consumed"] = lambda: sc.scalar(0) + out = out.transform_coords("dspacing", graph=graph, keep_intermediate=False) + out.coords.pop("_tag_positions_consumed", None) return DspacingData[RunType](out) -def to_dspacing_with_positions( - data: NormalizedByProtonCharge[RunType], - *, - sample: Optional[RawSample[RunType]] = None, - source: Optional[RawSource] = None, -) -> DspacingData[RunType]: +def powder_coordinate_transformation_graph() -> ElasticCoordTransformGraph: """ - Transform coordinates to d-spacing using detector positions. - - Computes d-spacing from time-of-flight stored in `data`. - - Attention - --------- - `data` may have a wavelength coordinate and dimension, - but those are discarded. - Only the stored time-of-flight is used, that is, any modifications to - the wavelength coordinate after it was computed from time-of-flight are lost. - - Raises - ------ - KeyError - If `data` does not contain a 'tof' coordinate. - - Parameters - ---------- - data: - Input data in tof or wavelength dimension. - Must have a tof coordinate. - sample: - Sample data with a position. - If not given, ``data`` must contain a 'sample_position' coordinate. - source: - Source data with a position. - If not given, ``data`` must contain a 'source_position' coordinate. + Generate a coordinate transformation graph for powder diffraction. Returns ------- : - A DataArray with the same data as the input and a 'dspacing' coordinate. + A dictionary with the graph for the transformation. """ - graph = { - **scn.conversion.graph.beamline.beamline(scatter=True), - **scn.conversion.graph.tof.elastic_dspacing('tof'), - } - if sample is not None: - graph['sample_position'] = lambda: sample['position'] - if source is not None: - graph['source_position'] = lambda: source['position'] - - out = _restore_tof_if_in_wavelength(data) - out = out.transform_coords('dspacing', graph=graph, keep_intermediate=False) - # Add coords to ensure the result is the same whether sample or source are - # coords in the input or separate function arguments. - if sample is not None: - out.coords['sample_position'] = sample['position'] - out.coords.set_aligned('sample_position', False) - if source is not None: - out.coords['source_position'] = source['position'] - out.coords.set_aligned('source_position', False) - - return DspacingData[RunType](out) + return ElasticCoordTransformGraph( + { + **scn.conversion.graph.beamline.beamline(scatter=True), + **scn.conversion.graph.tof.elastic("tof"), + } + ) def _restore_tof_if_in_wavelength(data: sc.DataArray) -> sc.DataArray: out = data.copy(deep=False) - outer = out.coords.pop('wavelength', None) + outer = out.coords.pop("wavelength", None) if out.bins is not None: - binned = out.bins.coords.pop('wavelength', None) + binned = out.bins.coords.pop("wavelength", None) else: binned = None if outer is not None or binned is not None: get_logger().info("Discarded coordinate 'wavelength' in favor of 'tof'.") - if 'wavelength' in out.dims: - out = out.rename_dims(wavelength='tof') + if "wavelength" in out.dims: + out = out.rename_dims(wavelength="tof") return out -providers = (to_dspacing_with_calibration,) +def add_scattering_coordinates_from_positions( + data: NormalizedByProtonCharge[RunType], graph: ElasticCoordTransformGraph +) -> DataWithScatteringCoordinates[RunType]: + """ + Add ``wavelength``, ``two_theta``, and ``dspacing`` coordinates to the data. + The input ``data`` must have a ``tof`` coordinate, as well as the necessary + positions of the beamline components (source, sample, detectors) to compute + the scattering coordinates. + + Parameters + ---------- + data: + Input data with a ``tof`` coordinate. + graph: + Coordinate transformation graph. + """ + out = data.transform_coords( + ["two_theta", "wavelength", "dspacing"], graph=graph, keep_intermediate=False + ) + return DataWithScatteringCoordinates[RunType](out) + + +providers_with_calibration = (to_dspacing_with_calibration,) """Sciline providers for coordinate transformations.""" + +providers_with_positions = ( + powder_coordinate_transformation_graph, + add_scattering_coordinates_from_positions, +) diff --git a/src/ess/powder/correction.py b/src/ess/powder/correction.py index d950f993..46363e7b 100644 --- a/src/ess/powder/correction.py +++ b/src/ess/powder/correction.py @@ -12,11 +12,12 @@ from .smoothing import lowpass from .types import ( AccumulatedProtonCharge, - DspacingBins, FilteredData, - FocussedData, + FocussedDataDspacing, + FocussedDataDspacingTwoTheta, + IofDspacing, + IofDspacingTwoTheta, NormalizedByProtonCharge, - NormalizedByVanadium, RunType, SampleRun, UncertaintyBroadcastMode, @@ -56,9 +57,9 @@ def normalize_by_monitor( : `data` normalized by a monitor. """ - if 'wavelength' not in monitor.coords: + if "wavelength" not in monitor.coords: monitor = monitor.transform_coords( - 'wavelength', + "wavelength", graph={**beamline.beamline(scatter=False), **tof.elastic("tof")}, keep_inputs=False, keep_intermediate=False, @@ -73,18 +74,30 @@ def normalize_by_monitor( "ess.powder.smoothing.lowpass with %s.", smooth_args, ) - monitor = lowpass(monitor, dim='wavelength', **smooth_args) - return data.bins / sc.lookup(func=monitor, dim='wavelength') + monitor = lowpass(monitor, dim="wavelength", **smooth_args) + return data.bins / sc.lookup(func=monitor, dim="wavelength") -def normalize_by_vanadium( - data: FocussedData[SampleRun], - vanadium: FocussedData[VanadiumRun], - edges: DspacingBins, +def _normalize_by_vanadium( + data: sc.DataArray, + vanadium: sc.DataArray, + uncertainty_broadcast_mode: UncertaintyBroadcastMode, +) -> sc.DataArray: + vanadium = broadcast_uncertainties(vanadium, uncertainty_broadcast_mode) + norm = vanadium.hist() + # Converting to unit 'one' because the division might produce a unit + # with a large scale if the proton charges in data and vanadium were + # measured with different units. + return (data / norm).to(unit="one", copy=False) + + +def normalize_by_vanadium_dspacing( + data: FocussedDataDspacing[SampleRun], + vanadium: FocussedDataDspacing[VanadiumRun], uncertainty_broadcast_mode: Optional[UncertaintyBroadcastMode] = None, -) -> NormalizedByVanadium: +) -> IofDspacing: """ - Normalize sample data by a vanadium measurement. + Normalize sample data by a vanadium measurement and return intensity vs d-spacing. Parameters ---------- @@ -92,24 +105,37 @@ def normalize_by_vanadium( Sample data. vanadium: Vanadium data. - edges: - `vanadium` is histogrammed into these bins before dividing the data by it. uncertainty_broadcast_mode: Choose how uncertainties of vanadium are broadcast to the sample data. Defaults to ``UncertaintyBroadcastMode.fail``. + """ + return IofDspacing( + _normalize_by_vanadium(data, vanadium, uncertainty_broadcast_mode) + ) - Returns - ------- - : - `data` normalized by `vanadium`. + +def normalize_by_vanadium_dspacing_and_two_theta( + data: FocussedDataDspacingTwoTheta[SampleRun], + vanadium: FocussedDataDspacingTwoTheta[VanadiumRun], + uncertainty_broadcast_mode: Optional[UncertaintyBroadcastMode] = None, +) -> IofDspacingTwoTheta: """ - vanadium = broadcast_uncertainties(vanadium, uncertainty_broadcast_mode) + Normalize sample data by a vanadium measurement and return intensity vs + (d-spacing, 2theta). - norm = sc.lookup(vanadium.hist({edges.dim: edges}), dim=edges.dim) - # Converting to unit 'one' because the division might produce a unit - # with a large scale if the proton charges in data and vanadium were - # measured with different units. - return (data.bins / norm).to(unit='one', copy=False) + Parameters + ---------- + data: + Sample data. + vanadium: + Vanadium data. + uncertainty_broadcast_mode: + Choose how uncertainties of vanadium are broadcast to the sample data. + Defaults to ``UncertaintyBroadcastMode.fail``. + """ + return IofDspacingTwoTheta( + _normalize_by_vanadium(data, vanadium, uncertainty_broadcast_mode) + ) def normalize_by_proton_charge( @@ -156,22 +182,22 @@ def merge_calibration(*, into: sc.DataArray, calibration: sc.Dataset) -> sc.Data dim = calibration.dim if not sc.identical(into.coords[dim], calibration.coords[dim]): raise ValueError( - f'Coordinate {dim} of calibration and target dataset do not agree.' + f"Coordinate {dim} of calibration and target dataset do not agree." ) out = into.copy(deep=False) - for name in ('difa', 'difc', 'tzero'): + for name in ("difa", "difc", "tzero"): if name in out.coords: raise ValueError( f"Cannot add calibration parameter '{name}' to data, " "there already is metadata with the same name." ) out.coords[name] = calibration[name].data - if 'calibration' in out.masks: + if "calibration" in out.masks: raise ValueError( "Cannot add calibration mask 'calibration' tp data, " "there already is a mask with the same name." ) - out.masks['calibration'] = calibration['mask'].data + out.masks["calibration"] = calibration["mask"].data return out @@ -207,8 +233,8 @@ def apply_lorentz_correction(da: sc.DataArray) -> sc.DataArray: # The implementation is optimized under the assumption that two_theta # is small and dspacing and the data are large. out = _shallow_copy(da) - dspacing = event_or_outer_coord(da, 'dspacing') - two_theta = event_or_outer_coord(da, 'two_theta') + dspacing = event_or_outer_coord(da, "dspacing") + two_theta = event_or_outer_coord(da, "two_theta") theta = 0.5 * two_theta d4 = dspacing.broadcast(sizes=out.sizes) ** 4 @@ -231,5 +257,9 @@ def _shallow_copy(da: sc.DataArray) -> sc.DataArray: return out -providers = (normalize_by_proton_charge, normalize_by_vanadium) +providers = ( + normalize_by_proton_charge, + normalize_by_vanadium_dspacing, + normalize_by_vanadium_dspacing_and_two_theta, +) """Sciline providers for powder diffraction corrections.""" diff --git a/src/ess/powder/external/powgen/beamline.py b/src/ess/powder/external/powgen/beamline.py index b17c1359..3c4e109f 100644 --- a/src/ess/powder/external/powgen/beamline.py +++ b/src/ess/powder/external/powgen/beamline.py @@ -6,9 +6,16 @@ import scipp as sc -from ...types import CalibrationData, DetectorDimensions, RawCalibrationData +from ...types import ( + CalibrationData, + NeXusDetectorDimensions, + NeXusDetectorName, + RawCalibrationData, +) from .types import DetectorInfo +DETECTOR_BANK_SIZES = {"powgen_detector": {"bank": 23, "column": 154, "row": 7}} + def map_detector_to_spectrum( data: sc.Dataset, *, detector_info: sc.Dataset @@ -29,24 +36,24 @@ def map_detector_to_spectrum( `data` with 'detector' coord and dim replaced by 'spectrum'. """ if not sc.identical( - data.coords['detector'].to( - dtype=detector_info.coords['detector'].dtype, copy=False + data.coords["detector"].to( + dtype=detector_info.coords["detector"].dtype, copy=False ), - detector_info.coords['detector'], + detector_info.coords["detector"], ): raise sc.CoordError( "The 'detector' coords of `data` and `detector_info` do not match." ) out = data.copy(deep=False) - del out.coords['detector'] + del out.coords["detector"] # Add 1 because spectrum numbers in the data start at 1 but # detector_info contains spectrum indices which start at 0. - out.coords['spectrum'] = detector_info.coords['spectrum'] + sc.index( - 1, dtype=detector_info.coords['spectrum'].dtype + out.coords["spectrum"] = detector_info.coords["spectrum"] + sc.index( + 1, dtype=detector_info.coords["spectrum"].dtype ) - return out.rename_dims({'detector': 'spectrum'}) + return out.rename_dims({"detector": "spectrum"}) def preprocess_calibration_data( @@ -61,9 +68,13 @@ def preprocess_calibration_data( return CalibrationData(map_detector_to_spectrum(data, detector_info=detector_info)) -def powgen_detector_dimensions() -> DetectorDimensions: +def powgen_detector_dimensions( + detector_name: NeXusDetectorName, +) -> NeXusDetectorDimensions[NeXusDetectorName]: """Dimensions used by POWGEN detectors.""" - return DetectorDimensions(('spectrum',)) + return NeXusDetectorDimensions[NeXusDetectorName]( + DETECTOR_BANK_SIZES[detector_name] + ) providers = (preprocess_calibration_data, powgen_detector_dimensions) diff --git a/src/ess/powder/external/powgen/data.py b/src/ess/powder/external/powgen/data.py index 13197bc8..8be5b5c0 100644 --- a/src/ess/powder/external/powgen/data.py +++ b/src/ess/powder/external/powgen/data.py @@ -9,38 +9,38 @@ AccumulatedProtonCharge, CalibrationFilename, Filename, + NeXusDetectorDimensions, + NeXusDetectorName, ProtonCharge, RawCalibrationData, RawDataAndMetadata, - RawDetectorData, + ReducibleDetectorData, RunType, SampleRun, ) from .types import DetectorInfo -_version = '1' - -__all__ = ['get_path'] +_version = "1" def _make_pooch(): import pooch return pooch.create( - path=pooch.os_cache('ess/powgen'), - env='ESS_DATA_DIR', - base_url='https://public.esss.dk/groups/scipp/ess/powgen/{version}/', + path=pooch.os_cache("ess/powgen"), + env="ESS_DATA_DIR", + base_url="https://public.esss.dk/groups/scipp/ess/powgen/{version}/", version=_version, registry={ # Files loadable with Mantid - 'PG3_4844_event.nxs': 'md5:d5ae38871d0a09a28ae01f85d969de1e', - 'PG3_4866_event.nxs': 'md5:3d543bc6a646e622b3f4542bc3435e7e', - 'PG3_5226_event.nxs': 'md5:58b386ebdfeb728d34fd3ba00a2d4f1e', - 'PG3_FERNS_d4832_2011_08_24.cal': 'md5:c181221ebef9fcf30114954268c7a6b6', + "PG3_4844_event.nxs": "md5:d5ae38871d0a09a28ae01f85d969de1e", + "PG3_4866_event.nxs": "md5:3d543bc6a646e622b3f4542bc3435e7e", + "PG3_5226_event.nxs": "md5:58b386ebdfeb728d34fd3ba00a2d4f1e", + "PG3_FERNS_d4832_2011_08_24.cal": "md5:c181221ebef9fcf30114954268c7a6b6", # Zipped Scipp HDF5 files - 'PG3_4844_event.zip': 'md5:a644c74f5e740385469b67431b690a3e', - 'PG3_4866_event.zip': 'md5:5bc49def987f0faeb212a406b92b548e', - 'PG3_FERNS_d4832_2011_08_24.zip': 'md5:0fef4ed5f61465eaaa3f87a18f5bb80d', + "PG3_4844_event.zip": "md5:a644c74f5e740385469b67431b690a3e", + "PG3_4866_event.zip": "md5:5bc49def987f0faeb212a406b92b548e", + "PG3_FERNS_d4832_2011_08_24.zip": "md5:0fef4ed5f61465eaaa3f87a18f5bb80d", }, ) @@ -48,7 +48,7 @@ def _make_pooch(): _pooch = _make_pooch() -def get_path(name: str, unzip: bool = False) -> str: +def _get_path(name: str) -> str: """ Return the path to a data file bundled with scippneutron. @@ -57,38 +57,39 @@ def get_path(name: str, unzip: bool = False) -> str: """ import pooch - return _pooch.fetch(name, processor=pooch.Unzip() if unzip else None) + if name.endswith(".zip"): + (path,) = _pooch.fetch(name, processor=pooch.Unzip()) + else: + path = _pooch.fetch(name) + return path -def mantid_sample_file() -> str: - return get_path('PG3_4844_event.nxs') +def powgen_tutorial_mantid_sample_file() -> str: + return _get_path("PG3_4844_event.nxs") -def mantid_vanadium_file() -> str: - return get_path('PG3_4866_event.nxs') +def powgen_tutorial_mantid_vanadium_file() -> str: + return _get_path("PG3_4866_event.nxs") -def mantid_empty_instrument_file() -> str: - return get_path('PG3_5226_event.nxs') +def powgen_tutorial_mantid_empty_instrument_file() -> str: + return _get_path("PG3_5226_event.nxs") -def mantid_calibration_file() -> str: - return get_path('PG3_FERNS_d4832_2011_08_24.cal') +def powgen_tutorial_mantid_calibration_file() -> str: + return _get_path("PG3_FERNS_d4832_2011_08_24.cal") -def sample_file() -> str: - (path,) = get_path('PG3_4844_event.zip', unzip=True) - return path +def powgen_tutorial_sample_file() -> str: + return _get_path("PG3_4844_event.zip") -def vanadium_file() -> str: - (path,) = get_path('PG3_4866_event.zip', unzip=True) - return path +def powgen_tutorial_vanadium_file() -> str: + return _get_path("PG3_4866_event.zip") -def calibration_file() -> str: - (path,) = get_path('PG3_FERNS_d4832_2011_08_24.zip', unzip=True) - return path +def powgen_tutorial_calibration_file() -> str: + return _get_path("PG3_FERNS_d4832_2011_08_24.zip") def pooch_load(filename: Filename[RunType]) -> RawDataAndMetadata[RunType]: @@ -99,42 +100,41 @@ def pooch_load(filename: Filename[RunType]) -> RawDataAndMetadata[RunType]: The loaded data holds both the events and any metadata from the file. """ - if filename.endswith('.zip'): - (path,) = get_path(filename, unzip=True) - else: - path = get_path(filename) - return RawDataAndMetadata[RunType](sc.io.load_hdf5(path)) + return RawDataAndMetadata[RunType](sc.io.load_hdf5(filename)) def pooch_load_calibration(filename: CalibrationFilename) -> RawCalibrationData: """Load the calibration data for the POWGEN test data.""" - if filename.endswith('.zip'): - (path,) = get_path(filename, unzip=True) - else: - path = get_path(filename) - return RawCalibrationData(sc.io.load_hdf5(path)) + return RawCalibrationData(sc.io.load_hdf5(filename)) -def extract_raw_data(dg: RawDataAndMetadata[RunType]) -> RawDetectorData[RunType]: +def extract_raw_data( + dg: RawDataAndMetadata[RunType], sizes: NeXusDetectorDimensions[NeXusDetectorName] +) -> ReducibleDetectorData[RunType]: """Return the events from a loaded data group.""" - return RawDetectorData[RunType](dg['data']) + # Remove the tof binning and dimension, as it is not needed and it gets in the way + # of masking. + out = dg["data"].squeeze() + out.coords.pop("tof", None) + out = out.fold(dim="spectrum", sizes=sizes) + return ReducibleDetectorData[RunType](out) def extract_detector_info(dg: RawDataAndMetadata[SampleRun]) -> DetectorInfo: """Return the detector info from a loaded data group.""" - return DetectorInfo(dg['detector_info']) + return DetectorInfo(dg["detector_info"]) def extract_proton_charge(dg: RawDataAndMetadata[RunType]) -> ProtonCharge[RunType]: """Return the proton charge from a loaded data group.""" - return ProtonCharge[RunType](dg['proton_charge']) + return ProtonCharge[RunType](dg["proton_charge"]) def extract_accumulated_proton_charge( - data: RawDetectorData[RunType], + data: ReducibleDetectorData[RunType], ) -> AccumulatedProtonCharge[RunType]: """Return the stored accumulated proton charge from a loaded data group.""" - return AccumulatedProtonCharge[RunType](data.coords['gd_prtn_chrg']) + return AccumulatedProtonCharge[RunType](data.coords["gd_prtn_chrg"]) providers = ( diff --git a/src/ess/powder/filtering.py b/src/ess/powder/filtering.py index 45f9e247..f49ab736 100644 --- a/src/ess/powder/filtering.py +++ b/src/ess/powder/filtering.py @@ -12,15 +12,14 @@ import scipp as sc -from ._util import elem_dtype, elem_unit, event_or_outer_coord -from .types import FilteredData, RawDetectorData, RunType, TofCroppedData, ValidTofRange +from .types import FilteredData, ReducibleDetectorData, RunType def _equivalent_bin_indices(a, b) -> bool: - a_begin = a.bins.constituents['begin'].flatten(to='') - a_end = a.bins.constituents['end'].flatten(to='') - b_begin = b.bins.constituents['begin'].flatten(to='') - b_end = b.bins.constituents['end'].flatten(to='') + a_begin = a.bins.constituents["begin"].flatten(to="") + a_end = a.bins.constituents["end"].flatten(to="") + b_begin = b.bins.constituents["begin"].flatten(to="") + b_end = b.bins.constituents["end"].flatten(to="") non_empty = a_begin != a_end return ( sc.all((a_begin == b_begin)[non_empty]).value @@ -33,10 +32,10 @@ def _temporary_bin_coord(data: sc.DataArray, name: str, coord: sc.Variable) -> N if not _equivalent_bin_indices(data, coord): raise ValueError("data and coord do not have equivalent bin indices") coord = sc.bins( - data=coord.bins.constituents['data'], - begin=data.bins.coords['pulse_time'].bins.constituents['begin'], - end=data.bins.coords['pulse_time'].bins.constituents['end'], - dim=coord.bins.constituents['dim'], + data=coord.bins.constituents["data"], + begin=data.bins.coords["pulse_time"].bins.constituents["begin"], + end=data.bins.coords["pulse_time"].bins.constituents["end"], + dim=coord.bins.constituents["dim"], ) data.bins.coords[name] = coord yield @@ -46,7 +45,7 @@ def _temporary_bin_coord(data: sc.DataArray, name: str, coord: sc.Variable) -> N # TODO non-monotonic proton charge -> raise? def _with_pulse_time_edges(da: sc.DataArray, dim: str) -> sc.DataArray: pulse_time = da.coords[dim] - one = sc.scalar(1, dtype='int64', unit=pulse_time.unit) + one = sc.scalar(1, dtype="int64", unit=pulse_time.unit) lo = pulse_time[0] - one hi = pulse_time[-1] + one mid = sc.midpoints(pulse_time) @@ -64,43 +63,16 @@ def remove_bad_pulses( good_pulse = _with_pulse_time_edges(proton_charge >= min_charge, proton_charge.dim) with _temporary_bin_coord( data, - 'good_pulse', + "good_pulse", sc.lookup(good_pulse, good_pulse.dim)[data.bins.coords[good_pulse.dim]], ): - filtered = data.group(sc.array(dims=['good_pulse'], values=[True])) - filtered = filtered.squeeze('good_pulse').copy(deep=False) - del filtered.coords['good_pulse'] + filtered = data.group(sc.array(dims=["good_pulse"], values=[True])) + filtered = filtered.squeeze("good_pulse").copy(deep=False) + del filtered.coords["good_pulse"] return filtered -def crop_tof( - data: RawDetectorData[RunType], tof_range: ValidTofRange -) -> TofCroppedData[RunType]: - """Remove events outside the specified TOF range. - - Parameters - ---------- - data: - Data to be cropped. - Expected to have a coordinate called `'tof'`. - tof_range: - 1d, len-2 variable containing the lower and upper bounds for - time-of-flight. - - Returns - ------- - : - Cropped data. - """ - tof = event_or_outer_coord(data, 'tof') - tof_unit = elem_unit(tof) - tof_dtype = elem_dtype(tof) - return TofCroppedData[RunType]( - data.bin(tof=tof_range.to(unit=tof_unit, dtype=tof_dtype)) - ) - - -def filter_events(data: TofCroppedData[RunType]) -> FilteredData[RunType]: +def filter_events(data: ReducibleDetectorData[RunType]) -> FilteredData[RunType]: """Remove bad events. Attention @@ -124,8 +96,5 @@ def filter_events(data: TofCroppedData[RunType]) -> FilteredData[RunType]: return FilteredData[RunType](data) -providers = ( - crop_tof, - filter_events, -) +providers = (filter_events,) """Sciline providers for event filtering.""" diff --git a/src/ess/powder/grouping.py b/src/ess/powder/grouping.py index 94d459a4..ecd74092 100644 --- a/src/ess/powder/grouping.py +++ b/src/ess/powder/grouping.py @@ -2,90 +2,34 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) """Grouping and merging of pixels / voxels.""" -from scippneutron.conversion.graph import beamline - from .types import ( - DetectorDimensions, DspacingBins, - DspacingData, - DspacingHistogram, - FocussedData, - NormalizedByVanadium, + FocussedDataDspacing, + FocussedDataDspacingTwoTheta, + MaskedData, RunType, TwoThetaBins, ) -def group_by_two_theta( - data: DspacingData[RunType], edges: TwoThetaBins -) -> FocussedData[RunType]: - """ - Group data into two_theta bins. - - Parameters - ---------- - data: - Input data array with events. Must contain a coord called 'two_theta' - or coords that can be used to compute it. - edges: - Bin edges in two_theta. `data` is grouped into those bins. - - Returns - ------- - : - `data` grouped into two_theta bins. - """ - out = data.transform_coords('two_theta', graph=beamline.beamline(scatter=True)) - return FocussedData[RunType]( - out.bin(two_theta=edges.to(unit=out.coords['two_theta'].unit, copy=False)) - ) - - -def merge_all_pixels( - data: DspacingData[RunType], dims: DetectorDimensions -) -> FocussedData[RunType]: - """Combine all pixels (spectra) of the detector. - - Parameters - ---------- - data: - Input data with pixel dimensions. - dims: - The dimensions to reduce over. - Corresponds to the pixel dimensions of the detector. - - Returns - ------- - : - The input without pixel dimensions. - """ - # Put the dims into the same order as in the data. - # See https://github.com/scipp/scipp/issues/3408 - to_concat = tuple(dim for dim in data.dims if dim in dims) - return FocussedData(data.bins.concat(to_concat)) - - -def finalize_histogram( - data: NormalizedByVanadium, edges: DspacingBins -) -> DspacingHistogram: - """Finalize the d-spacing histogram. - - Histograms the input data into the given d-spacing bins. +def focus_data_dspacing( + data: MaskedData[RunType], + dspacing_bins: DspacingBins, +) -> FocussedDataDspacing[RunType]: + out = data.bins.concat().bin({dspacing_bins.dim: dspacing_bins}) + return FocussedDataDspacing[RunType](out) - Parameters - ---------- - data: - Data to be histogrammed. - edges: - Bin edges in d-spacing. - Returns - ------- - : - Histogrammed data. - """ - return DspacingHistogram(data.hist(dspacing=edges)) +def focus_data_dspacing_and_two_theta( + data: MaskedData[RunType], + dspacing_bins: DspacingBins, + twotheta_bins: TwoThetaBins, +) -> FocussedDataDspacingTwoTheta[RunType]: + bins = {twotheta_bins.dim: twotheta_bins, dspacing_bins.dim: dspacing_bins} + if "two_theta" in data.bins.coords: + data = data.bins.concat() + return FocussedDataDspacingTwoTheta[RunType](data.bin(**bins)) -providers = (merge_all_pixels, finalize_histogram) +providers = (focus_data_dspacing, focus_data_dspacing_and_two_theta) """Sciline providers for grouping pixels.""" diff --git a/src/ess/powder/masking.py b/src/ess/powder/masking.py new file mode 100644 index 00000000..cd51ded8 --- /dev/null +++ b/src/ess/powder/masking.py @@ -0,0 +1,73 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2024 Scipp contributors (https://github.com/scipp) +""" +Masking functions for the powder workflow. +""" + +from typing import Optional + +import numpy as np +import scipp as sc + +from .types import ( + DataWithScatteringCoordinates, + MaskedData, + MaskedDetectorIDs, + PixelMaskFilename, + RunType, + TofMask, + TwoThetaMask, + WavelengthMask, +) + + +def read_pixel_masks( + filename: Optional[PixelMaskFilename] = None, +) -> MaskedDetectorIDs: + """Read a pixel mask from a Scipp hdf5 file. + + Parameters + ---------- + filename: + Path to the hdf5 file. + """ + masked_ids = {} + if filename is not None: + masked_ids = {filename: sc.io.load_hdf5(filename)} + return MaskedDetectorIDs(masked_ids) + + +def apply_masks( + data: DataWithScatteringCoordinates[RunType], + masked_pixel_ids: MaskedDetectorIDs, + tof_mask_func: Optional[TofMask] = None, + wavelength_mask_func: Optional[WavelengthMask] = None, + two_theta_mask_func: Optional[TwoThetaMask] = None, +) -> MaskedData[RunType]: + """ """ + out = data.copy(deep=False) + if len(masked_pixel_ids) > 0: + key = ( + set(out.coords.keys()) & {"detector_number", "detector_id", "spectrum"} + ).pop() + ids = out.coords[key] + for name, masked in masked_pixel_ids.items(): + mask = sc.zeros(sizes=ids.sizes, dtype="bool") + mask.values[np.isin(ids.values, masked.values)] = True + out.masks[name] = mask + + for dim, mask in { + "tof": tof_mask_func, + "wavelength": wavelength_mask_func, + "two_theta": two_theta_mask_func, + }.items(): + if mask is not None: + if dim in out.bins.coords: + out.bins.masks[dim] = mask(out.bins.coords[dim]) + else: + out.masks[dim] = mask(out.coords[dim]) + + return MaskedData[RunType](out) + + +providers = (read_pixel_masks, apply_masks) diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index 5a4310c5..8887af7e 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -8,8 +8,7 @@ """ from enum import Enum -from pathlib import Path -from typing import Any, NewType, TypeVar +from typing import Any, Callable, Dict, NewType, TypeVar import sciline import scipp as sc @@ -17,35 +16,28 @@ # 1 TypeVars used to parametrize the generic parts of the workflow # 1.1 Run types -EmptyCanRun = NewType('EmptyCanRun', int) +EmptyCanRun = NewType("EmptyCanRun", int) """Empty sample can run.""" -EmptyInstrumentRun = NewType('EmptyInstrumentRun', int) +EmptyInstrumentRun = NewType("EmptyInstrumentRun", int) """Empty instrument run.""" -SampleRun = NewType('SampleRun', int) +SampleRun = NewType("SampleRun", int) """Sample run.""" -VanadiumRun = NewType('VanadiumRun', int) +VanadiumRun = NewType("VanadiumRun", int) """Vanadium run.""" -RunType = TypeVar('RunType', EmptyInstrumentRun, SampleRun, VanadiumRun) +RunType = TypeVar("RunType", EmptyInstrumentRun, SampleRun, VanadiumRun) """TypeVar used for specifying the run.""" # 2 Workflow parameters -CalibrationFilename = NewType('CalibrationFilename', str) +CalibrationFilename = NewType("CalibrationFilename", str) """Filename of the instrument calibration file.""" -# In Python 3.11, this can be replaced with a StrEnum -class DetectorName(str, Enum): - """Name of a detector.""" +NeXusDetectorName = NewType("NeXusDetectorName", str) +"""Name of detector entry in NeXus file""" - mantle = 'mantle' - high_resolution = 'high_resolution' - endcap_backward = 'endcap_backward' - endcap_forward = 'endcap_forward' - - -DspacingBins = NewType('DSpacingBins', sc.Variable) +DspacingBins = NewType("DSpacingBins", sc.Variable) """Bin edges for d-spacing.""" @@ -53,14 +45,10 @@ class Filename(sciline.Scope[RunType, str], str): """Name of an input file.""" -class FilePath(sciline.Scope[RunType, Path], Path): - """Path to an input file on disk.""" - - -OutFilename = NewType('OutFilename', str) +OutFilename = NewType("OutFilename", str) """Filename of the output.""" -TwoThetaBins = NewType('TwoThetaBins', sc.Variable) +TwoThetaBins = NewType("TwoThetaBins", sc.Variable) """Bin edges for grouping in 2theta. This is used by an alternative focussing step that groups detector @@ -68,14 +56,14 @@ class FilePath(sciline.Scope[RunType, Path], Path): """ UncertaintyBroadcastMode = Enum( - 'UncertaintyBroadcastMode', ['drop', 'upper_bound', 'fail'] + "UncertaintyBroadcastMode", ["drop", "upper_bound", "fail"] ) """Mode for broadcasting uncertainties. See https://doi.org/10.3233/JNR-220049 for context. """ -ValidTofRange = NewType('ValidTofRange', sc.Variable) +ValidTofRange = NewType("ValidTofRange", sc.Variable) """Min and max tof value of the instrument.""" # 3 Workflow (intermediate) results @@ -90,11 +78,20 @@ def __init__(self, *args: Any, **kwargs: Any) -> None: super().__init__(*args, **kwargs) -CalibrationData = NewType('CalibrationData', sc.Dataset) +CalibrationData = NewType("CalibrationData", sc.Dataset) """Detector calibration data.""" +DataFolder = NewType("DataFolder", str) -class DetectorDimensions(sciline.Scope[DetectorName, tuple[str, ...]], tuple[str, ...]): + +class DataWithScatteringCoordinates(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data with scattering coordinates computed for all events: wavelength, 2theta, + d-spacing.""" + + +class NeXusDetectorDimensions( + sciline.Scope[NeXusDetectorName, Dict[str, int]], Dict[str, int] +): """Logical detector dimensions.""" @@ -106,31 +103,59 @@ class DspacingDataWithoutVariances(sciline.Scope[RunType, sc.DataArray], sc.Data """Data converted to d-spacing where variances where removed.""" -DspacingHistogram = NewType('DspacingHistogram', sc.DataArray) +DspacingHistogram = NewType("DspacingHistogram", sc.DataArray) """Histogrammed intensity vs d-spacing.""" +ElasticCoordTransformGraph = NewType("ElasticCoordTransformGraph", dict) +"""Graph for transforming coordinates in elastic scattering.""" + class FilteredData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Raw data without invalid events.""" -class FocussedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): +class FocussedDataDspacing(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Intensity vs d-spacing after focussing pixels.""" +class FocussedDataDspacingTwoTheta(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Intensity vs (d-spacing, 2theta) after focussing pixels.""" + + +IofDspacing = NewType("IofDspacing", sc.DataArray) +"""Data that has been normalized by a vanadium run.""" + +IofDspacingTwoTheta = NewType("IofDspacingTwoTheta", sc.DataArray) +"""Data that has been normalized by a vanadium run, and grouped into 2theta bins.""" + + +class LoadedNeXusDetector(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): + """Detector data, loaded from a NeXus file, containing not only neutron events + but also pixel shape information, transformations, ...""" + + +class MaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data with masked pixels, tof regions, wavelength regions, 2theta regions, or + dspacing regions.""" + + +MaskedDetectorIDs = NewType("MaskedDetectorIDs", Dict[str, sc.Variable]) +"""1-D variable listing all masked detector IDs.""" + + class NormalizedByProtonCharge(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Data that has been normalized by proton charge.""" -NormalizedByVanadium = NewType('NormalizedByVanadium', sc.DataArray) -"""Data that has been normalized by a vanadium run.""" +PixelMaskFilename = NewType("PixelMaskFilename", str) +"""Filename of a pixel mask.""" class ProtonCharge(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Time-dependent proton charge.""" -RawCalibrationData = NewType('RawCalibrationData', sc.Dataset) +RawCalibrationData = NewType("RawCalibrationData", sc.Dataset) """Calibration data as loaded from file, needs preprocessing before using.""" @@ -150,12 +175,32 @@ class RawSample(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): """Raw data from a loaded sample.""" -RawSource = NewType('RawSource', sc.DataGroup) -"""Raw data from a loaded neutron source.""" +class RawSource(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): + """Raw data from a loaded neutron source.""" + + +class ReducibleDetectorData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data that is in a state ready for reduction.""" + + +class SamplePosition(sciline.Scope[RunType, sc.Variable], sc.Variable): + """Sample position""" + + +class SourcePosition(sciline.Scope[RunType, sc.Variable], sc.Variable): + """Source position""" + + +TofMask = NewType("TofMask", Callable) +"""TofMask is a callable that returns a mask for a given TofData.""" + + +TwoThetaMask = NewType("TwoThetaMask", Callable) +"""TwoThetaMask is a callable that returns a mask for a given TwoThetaData.""" -class TofCroppedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Raw data cropped to the valid TOF range.""" +WavelengthMask = NewType("WavelengthMask", Callable) +"""WavelengthMask is a callable that returns a mask for a given WavelengthData.""" del sc, sciline, NewType, TypeVar diff --git a/tests/dream/geant4_reduction_test.py b/tests/dream/geant4_reduction_test.py new file mode 100644 index 00000000..085fa5a0 --- /dev/null +++ b/tests/dream/geant4_reduction_test.py @@ -0,0 +1,138 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) + +import pytest +import sciline +import scipp as sc +from ess.powder.types import ( + AccumulatedProtonCharge, + DspacingBins, + EmptyCanRun, + Filename, + IofDspacing, + IofDspacingTwoTheta, + MaskedData, + NeXusDetectorName, + NormalizedByProtonCharge, + RawSample, + RawSource, + SampleRun, + TofMask, + TwoThetaBins, + TwoThetaMask, + UncertaintyBroadcastMode, + VanadiumRun, + WavelengthMask, +) + + +@pytest.fixture() +def providers(): + from ess import powder + from ess.dream.io.geant4 import providers as geant4_providers + + return [*powder.providers, *geant4_providers] + + +@pytest.fixture(params=["mantle", "endcap_backward", "endcap_forward"]) +def params(request): + from ess import dream + + # Not available in simulated data + sample = sc.DataGroup(position=sc.vector([0.0, 0.0, 0.0], unit='mm')) + source = sc.DataGroup(position=sc.vector([-3.478, 0.0, -76550], unit='mm')) + charge = sc.scalar(1.0, unit='µAh') + + return { + NeXusDetectorName: request.param, + Filename[SampleRun]: dream.data.simulated_diamond_sample(), + Filename[VanadiumRun]: dream.data.simulated_vanadium_sample(), + Filename[EmptyCanRun]: dream.data.simulated_empty_can(), + UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, + DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 201, unit='angstrom'), + TofMask: lambda x: (x < sc.scalar(0.0, unit='ns')) + | (x > sc.scalar(86e6, unit='ns')), + RawSample[SampleRun]: sample, + RawSample[VanadiumRun]: sample, + RawSource[SampleRun]: source, + RawSource[VanadiumRun]: source, + AccumulatedProtonCharge[SampleRun]: charge, + AccumulatedProtonCharge[VanadiumRun]: charge, + } + + +def test_can_create_pipeline(providers, params): + sciline.Pipeline(providers, params=params) + + +def test_pipeline_can_compute_dspacing_result(providers, params): + pipeline = sciline.Pipeline(providers, params=params) + result = pipeline.compute(IofDspacing) + assert result.sizes == { + 'dspacing': len(params[DspacingBins]) - 1, + } + assert sc.identical(result.coords['dspacing'], params[DspacingBins]) + + +def test_workflow_is_deterministic(providers, params): + pipeline = sciline.Pipeline(providers, params=params) + # This is Sciline's default scheduler, but we want to be explicit here + scheduler = sciline.scheduler.DaskScheduler() + graph = pipeline.get(IofDspacing, scheduler=scheduler) + reference = graph.compute().data + result = graph.compute().data + assert sc.identical(sc.values(result), sc.values(reference)) + + +def test_pipeline_can_compute_intermediate_results(providers, params): + pipeline = sciline.Pipeline(providers, params=params) + result = pipeline.compute(NormalizedByProtonCharge[SampleRun]) + assert set(result.dims) == {'segment', 'wire', 'counter', 'strip', 'module'} + + +def test_pipeline_group_by_two_theta(providers, params): + params[TwoThetaBins] = sc.linspace( + dim='two_theta', unit='rad', start=0.8, stop=2.4, num=17 + ) + pipeline = sciline.Pipeline(providers, params=params) + result = pipeline.compute(IofDspacingTwoTheta) + assert result.sizes == { + 'two_theta': 16, + 'dspacing': len(params[DspacingBins]) - 1, + } + assert sc.identical(result.coords['dspacing'], params[DspacingBins]) + assert sc.allclose(result.coords['two_theta'], params[TwoThetaBins]) + + +def test_pipeline_wavelength_masking(providers, params): + wmin = sc.scalar(0.18, unit="angstrom") + wmax = sc.scalar(0.21, unit="angstrom") + params[WavelengthMask] = lambda x: (x > wmin) & (x < wmax) + pipeline = sciline.Pipeline(providers, params=params) + masked_sample = pipeline.compute(MaskedData[SampleRun]) + assert 'wavelength' in masked_sample.bins.masks + sum_in_masked_region = ( + masked_sample.bin(wavelength=sc.concat([wmin, wmax], dim='wavelength')) + .sum() + .data + ) + assert sc.allclose( + sum_in_masked_region, + sc.scalar(0.0, unit=sum_in_masked_region.unit), + ) + + +def test_pipeline_two_theta_masking(providers, params): + tmin = sc.scalar(1.0, unit="rad") + tmax = sc.scalar(1.2, unit="rad") + params[TwoThetaMask] = lambda x: (x > tmin) & (x < tmax) + pipeline = sciline.Pipeline(providers, params=params) + masked_sample = pipeline.compute(MaskedData[SampleRun]) + assert 'two_theta' in masked_sample.bins.masks + sum_in_masked_region = ( + masked_sample.bin(two_theta=sc.concat([tmin, tmax], dim='two_theta')).sum().data + ) + assert sc.allclose( + sum_in_masked_region, + sc.scalar(0.0, unit=sum_in_masked_region.unit), + ) diff --git a/tests/dream/io/geant4_test.py b/tests/dream/io/geant4_test.py index 5c082f80..5b78cd25 100644 --- a/tests/dream/io/geant4_test.py +++ b/tests/dream/io/geant4_test.py @@ -11,30 +11,30 @@ import scipp as sc import scipp.testing from ess.dream import data, load_geant4_csv -from ess.powder.types import DetectorName, FilePath, RawDetectorData, SampleRun +from ess.powder.types import Filename, NeXusDetectorName, RawDetectorData, SampleRun -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def file_path(): - return data.get_path('data_dream0_new_hkl_Si_pwd.csv.zip') + return data.get_path("data_dream0_new_hkl_Si_pwd.csv.zip") -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def file_path_without_sans(): - return data.get_path('data_dream_with_sectors.csv.zip') + return data.get_path("data_dream_with_sectors.csv.zip") # Load file into memory only once -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def load_file(file_path): - with zipfile.ZipFile(file_path, 'r') as archive: + with zipfile.ZipFile(file_path, "r") as archive: return archive.read(archive.namelist()[0]) # Load file into memory only once -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def load_file_without_sans(file_path_without_sans): - with zipfile.ZipFile(file_path_without_sans, 'r') as archive: + with zipfile.ZipFile(file_path_without_sans, "r") as archive: return archive.read(archive.namelist()[0]) @@ -53,7 +53,7 @@ def assert_index_coord( ) -> None: assert coord.ndim == 1 assert coord.unit is None - assert coord.dtype == 'int64' + assert coord.dtype == "int64" if values is not None: assert set(np.unique(coord.values)) == values @@ -61,120 +61,120 @@ def assert_index_coord( def test_load_geant4_csv_loads_expected_structure(file): loaded = load_geant4_csv(file) assert isinstance(loaded, sc.DataGroup) - assert loaded.keys() == {'instrument'} + assert loaded.keys() == {"instrument"} - instrument = loaded['instrument'] + instrument = loaded["instrument"] assert isinstance(instrument, sc.DataGroup) assert instrument.keys() == { - 'mantle', - 'high_resolution', - 'sans', - 'endcap_forward', - 'endcap_backward', + "mantle", + "high_resolution", + "sans", + "endcap_forward", + "endcap_backward", } def test_load_geant4_csv_loads_expected_structure_without_sans(file_without_sans): loaded = load_geant4_csv(file_without_sans) assert isinstance(loaded, sc.DataGroup) - assert loaded.keys() == {'instrument'} + assert loaded.keys() == {"instrument"} - instrument = loaded['instrument'] + instrument = loaded["instrument"] assert isinstance(instrument, sc.DataGroup) assert instrument.keys() == { - 'mantle', - 'high_resolution', - 'endcap_forward', - 'endcap_backward', + "mantle", + "high_resolution", + "endcap_forward", + "endcap_backward", } @pytest.mark.parametrize( - 'key', ['mantle', 'high_resolution', 'endcap_forward', 'endcap_backward'] + "key", ["mantle", "high_resolution", "endcap_forward", "endcap_backward"] ) def test_load_gean4_csv_set_weights_to_one(file, key): - detector = load_geant4_csv(file)['instrument'][key]['events'] - events = detector.bins.constituents['data'].data + detector = load_geant4_csv(file)["instrument"][key]["events"] + events = detector.bins.constituents["data"].data sc.testing.assert_identical( - events, sc.ones(sizes=events.sizes, with_variances=True, unit='counts') + events, sc.ones(sizes=events.sizes, with_variances=True, unit="counts") ) def test_load_geant4_csv_mantle_has_expected_coords(file): # Only testing ranges that will not change in the future - mantle = load_geant4_csv(file)['instrument']['mantle']['events'] - assert_index_coord(mantle.coords['module']) - assert_index_coord(mantle.coords['segment']) - assert_index_coord(mantle.coords['counter']) - assert_index_coord(mantle.coords['wire'], values=set(range(1, 33))) - assert_index_coord(mantle.coords['strip'], values=set(range(1, 257))) - assert 'sector' not in mantle.coords + mantle = load_geant4_csv(file)["instrument"]["mantle"]["events"] + assert_index_coord(mantle.coords["module"]) + assert_index_coord(mantle.coords["segment"]) + assert_index_coord(mantle.coords["counter"]) + assert_index_coord(mantle.coords["wire"], values=set(range(1, 33))) + assert_index_coord(mantle.coords["strip"], values=set(range(1, 257))) + assert "sector" not in mantle.coords - assert 'sector' not in mantle.bins.coords - assert 'tof' in mantle.bins.coords - assert 'wavelength' in mantle.bins.coords - assert 'position' in mantle.bins.coords + assert "sector" not in mantle.bins.coords + assert "tof" in mantle.bins.coords + assert "wavelength" in mantle.bins.coords + assert "position" in mantle.bins.coords def test_load_geant4_csv_endcap_backward_has_expected_coords(file): - endcap = load_geant4_csv(file)['instrument']['endcap_backward']['events'] - assert_index_coord(endcap.coords['module']) - assert_index_coord(endcap.coords['segment']) - assert_index_coord(endcap.coords['counter']) - assert_index_coord(endcap.coords['wire'], values=set(range(1, 17))) - assert_index_coord(endcap.coords['strip'], values=set(range(1, 17))) - assert 'sector' not in endcap.coords + endcap = load_geant4_csv(file)["instrument"]["endcap_backward"]["events"] + assert_index_coord(endcap.coords["module"]) + assert_index_coord(endcap.coords["segment"]) + assert_index_coord(endcap.coords["counter"]) + assert_index_coord(endcap.coords["wire"], values=set(range(1, 17))) + assert_index_coord(endcap.coords["strip"], values=set(range(1, 17))) + assert "sector" not in endcap.coords - assert 'sector' not in endcap.bins.coords - assert 'tof' in endcap.bins.coords - assert 'wavelength' in endcap.bins.coords - assert 'position' in endcap.bins.coords + assert "sector" not in endcap.bins.coords + assert "tof" in endcap.bins.coords + assert "wavelength" in endcap.bins.coords + assert "position" in endcap.bins.coords def test_load_geant4_csv_endcap_forward_has_expected_coords(file): - endcap = load_geant4_csv(file)['instrument']['endcap_forward']['events'] - assert_index_coord(endcap.coords['module']) - assert_index_coord(endcap.coords['segment']) - assert_index_coord(endcap.coords['counter']) - assert_index_coord(endcap.coords['wire'], values=set(range(1, 17))) - assert_index_coord(endcap.coords['strip'], values=set(range(1, 17))) - assert 'sector' not in endcap.coords + endcap = load_geant4_csv(file)["instrument"]["endcap_forward"]["events"] + assert_index_coord(endcap.coords["module"]) + assert_index_coord(endcap.coords["segment"]) + assert_index_coord(endcap.coords["counter"]) + assert_index_coord(endcap.coords["wire"], values=set(range(1, 17))) + assert_index_coord(endcap.coords["strip"], values=set(range(1, 17))) + assert "sector" not in endcap.coords - assert 'sector' not in endcap.bins.coords - assert 'tof' in endcap.bins.coords - assert 'wavelength' in endcap.bins.coords - assert 'position' in endcap.bins.coords + assert "sector" not in endcap.bins.coords + assert "tof" in endcap.bins.coords + assert "wavelength" in endcap.bins.coords + assert "position" in endcap.bins.coords def test_load_geant4_csv_high_resolution_has_expected_coords(file): - hr = load_geant4_csv(file)['instrument']['high_resolution']['events'] - assert_index_coord(hr.coords['module']) - assert_index_coord(hr.coords['segment']) - assert_index_coord(hr.coords['counter']) - assert_index_coord(hr.coords['wire'], values=set(range(1, 17))) - assert_index_coord(hr.coords['strip'], values=set(range(1, 33))) - assert_index_coord(hr.coords['sector'], values=set(range(1, 5))) + hr = load_geant4_csv(file)["instrument"]["high_resolution"]["events"] + assert_index_coord(hr.coords["module"]) + assert_index_coord(hr.coords["segment"]) + assert_index_coord(hr.coords["counter"]) + assert_index_coord(hr.coords["wire"], values=set(range(1, 17))) + assert_index_coord(hr.coords["strip"], values=set(range(1, 33))) + assert_index_coord(hr.coords["sector"], values=set(range(1, 5))) - assert 'tof' in hr.bins.coords - assert 'wavelength' in hr.bins.coords - assert 'position' in hr.bins.coords + assert "tof" in hr.bins.coords + assert "wavelength" in hr.bins.coords + assert "position" in hr.bins.coords def test_load_geant4_csv_sans_has_expected_coords(file): - sans = load_geant4_csv(file)['instrument']['sans']['events'] - assert_index_coord(sans.coords['module']) - assert_index_coord(sans.coords['segment']) - assert_index_coord(sans.coords['counter']) + sans = load_geant4_csv(file)["instrument"]["sans"]["events"] + assert_index_coord(sans.coords["module"]) + assert_index_coord(sans.coords["segment"]) + assert_index_coord(sans.coords["counter"]) # check ranges only if csv file contains events from SANS detectors - if len(sans.coords['module'].values) > 0: - assert_index_coord(sans.coords['wire'], values=set(range(1, 17))) - assert_index_coord(sans.coords['strip'], values=set(range(1, 33))) - assert_index_coord(sans.coords['sector'], values=set(range(1, 5))) + if len(sans.coords["module"].values) > 0: + assert_index_coord(sans.coords["wire"], values=set(range(1, 17))) + assert_index_coord(sans.coords["strip"], values=set(range(1, 33))) + assert_index_coord(sans.coords["sector"], values=set(range(1, 5))) - assert 'tof' in sans.bins.coords - assert 'wavelength' in sans.bins.coords - assert 'position' in sans.bins.coords + assert "tof" in sans.bins.coords + assert "wavelength" in sans.bins.coords + assert "position" in sans.bins.coords def test_geant4_in_pipeline(file_path, file): @@ -182,8 +182,11 @@ def test_geant4_in_pipeline(file_path, file): pipeline = sciline.Pipeline( providers, - params={FilePath[SampleRun]: file_path, DetectorName: DetectorName('mantle')}, + params={ + Filename[SampleRun]: file_path, + NeXusDetectorName: NeXusDetectorName("mantle"), + }, ) detector = pipeline.compute(RawDetectorData[SampleRun]) - expected = load_geant4_csv(file)['instrument']['mantle']['events'] + expected = load_geant4_csv(file)["instrument"]["mantle"]["events"] sc.testing.assert_identical(detector, expected) diff --git a/tests/dream/io/nexus_test.py b/tests/dream/io/nexus_test.py index 6f7fe5f9..126197ac 100644 --- a/tests/dream/io/nexus_test.py +++ b/tests/dream/io/nexus_test.py @@ -1,76 +1,86 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) import pytest +import sciline from ess import dream +from ess.dream import nexus +from ess.powder.types import ( + Filename, + NeXusDetectorName, + RawDetectorData, + ReducibleDetectorData, + SampleRun, +) - -@pytest.fixture() -def filename(): - return dream.data.get_path('DREAM_nexus_sorted-2023-12-07.nxs') +bank_dims = {'wire', 'module', 'segment', 'strip', 'counter'} +hr_sans_dims = {'strip', 'other'} -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") -def test_load_nexus_loads_file(filename): - dg = dream.load_nexus(filename) - assert 'instrument' in dg - instr = dg['instrument'] - for name in ( - 'mantle', - 'endcap_backward', - 'endcap_forward', - 'high_resolution', - 'sans', - ): - assert f'{name}_detector' in instr - det = instr[f'{name}_detector'] - assert 'pixel_shape' not in det +@pytest.fixture() +def providers(): + return ( + nexus.dummy_load_sample, + nexus.load_nexus_source, + nexus.load_nexus_detector, + nexus.get_sample_position, + nexus.get_source_position, + nexus.get_detector_data, + nexus.patch_detector_data, + ) -def test_load_nexus_fails_if_entry_not_found(filename): - with pytest.raises(KeyError): - dream.load_nexus(filename, entry='foo') +@pytest.fixture( + params=[ + 'mantle_detector', + 'endcap_backward_detector', + 'endcap_forward_detector', + 'high_resolution_detector', + # TODO: the 'sans_detector' is strange in the current files + ] +) +def params(request): + params = { + Filename[SampleRun]: dream.data.get_path('DREAM_nexus_sorted-2023-12-07.nxs'), + NeXusDetectorName: request.param, + } + return params -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") -def test_load_nexus_folds_detectors_by_default(filename): - dg = dream.load_nexus(filename) - instr = dg['instrument'] - # sans_detector is not populated in the current files - for name in ('mantle', 'endcap_backward', 'endcap_forward', 'high_resolution'): - det = instr[f'{name}_detector'] - # There may be other dims, but some are irregular and this may be subject to - # change - assert 'strip' in det.dims +def test_can_load_nexus_detector_data(providers, params): + pipeline = sciline.Pipeline(params=params, providers=providers) + result = pipeline.compute(RawDetectorData[SampleRun]) + assert ( + set(result.dims) == hr_sans_dims + if params[NeXusDetectorName] + in ( + 'high_resolution_detector', + 'sans_detector', + ) + else bank_dims + ) -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") -def test_load_nexus_with_disabled_fold(filename): - dg = dream.load_nexus(filename, fold_detectors=False) - instr = dg['instrument'] - for name in ('mantle', 'endcap_backward', 'endcap_forward', 'high_resolution'): - det = instr[f'{name}_detector'] - assert det.dims == ('detector_number',) +def test_load_fails_with_bad_detector_name(providers): + params = { + Filename[SampleRun]: dream.data.get_path('DREAM_nexus_sorted-2023-12-07.nxs'), + NeXusDetectorName: 'bad_detector', + } + pipeline = sciline.Pipeline(params=params, providers=providers) + with pytest.raises(KeyError, match='bad_detector'): + pipeline.compute(RawDetectorData[SampleRun]) -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") -def test_load_nexus_with_pixel_shape(filename): - dg = dream.load_nexus(filename, load_pixel_shape=True) - assert 'instrument' in dg - instr = dg['instrument'] - # sans_detector is not populated in the current files - for name in ('mantle', 'endcap_backward', 'endcap_forward', 'high_resolution'): - assert f'{name}_detector' in instr - det = instr[f'{name}_detector'] - assert 'pixel_shape' in det +def test_patch_nexus_detector_data(providers, params): + pipeline = sciline.Pipeline(params=params, providers=providers) + result = pipeline.compute(ReducibleDetectorData[SampleRun]) + assert ( + set(result.dims) == hr_sans_dims + if params[NeXusDetectorName] + in ( + 'high_resolution_detector', + 'sans_detector', + ) + else bank_dims + ) + assert "source_position" in result.coords + assert "sample_position" in result.coords diff --git a/tests/powder/conversion_test.py b/tests/powder/conversion_test.py index 08af8a77..281f37f3 100644 --- a/tests/powder/conversion_test.py +++ b/tests/powder/conversion_test.py @@ -6,8 +6,8 @@ import scipp.testing import scippneutron as scn from ess.powder.conversion import ( + add_scattering_coordinates_from_positions, to_dspacing_with_calibration, - to_dspacing_with_positions, ) @@ -138,7 +138,7 @@ def test_dspacing_with_calibration_does_not_use_positions(calibration): ) -def test_dspacing_with_positions(): +def test_add_scattering_coordinates_from_positions(): position = sc.vectors( dims=['spectrum'], values=np.arange(14 * 3).reshape((14, 3)), unit='m' ) @@ -149,20 +149,17 @@ def test_dspacing_with_positions(): coords={ 'position': position, 'tof': sc.linspace('tof', 1.0, 1000.0, 27, unit='us'), + 'sample_position': sample_position, + 'source_position': source_position, }, ) - dspacing = to_dspacing_with_positions( - tof, - sample=sc.DataGroup(position=sample_position), - source=sc.DataGroup(position=source_position), - ) - graph = { **scn.conversion.graph.beamline.beamline(scatter=True), - **scn.conversion.graph.tof.elastic_dspacing('tof'), + **scn.conversion.graph.tof.elastic('tof'), } - tof = tof.assign_coords( - {'sample_position': sample_position, 'source_position': source_position} - ) - expected = tof.transform_coords('dspacing', graph=graph, keep_intermediate=False) - sc.testing.assert_identical(dspacing, expected) + + result = add_scattering_coordinates_from_positions(tof, graph) + + assert 'wavelength' in result.coords + assert 'two_theta' in result.coords + assert 'dspacing' in result.coords diff --git a/tests/powder/external/powgen/powgen_reduction_test.py b/tests/powder/external/powgen/powgen_reduction_test.py index 64bba88f..a723a0e4 100644 --- a/tests/powder/external/powgen/powgen_reduction_test.py +++ b/tests/powder/external/powgen/powgen_reduction_test.py @@ -7,14 +7,19 @@ from ess.powder.types import ( CalibrationFilename, DspacingBins, - DspacingHistogram, Filename, + IofDspacing, + IofDspacingTwoTheta, + MaskedData, + NeXusDetectorName, NormalizedByProtonCharge, SampleRun, + TofMask, TwoThetaBins, + TwoThetaMask, UncertaintyBroadcastMode, - ValidTofRange, VanadiumRun, + WavelengthMask, ) @@ -28,13 +33,17 @@ def providers(): @pytest.fixture() def params(): + from ess.powder.external import powgen + return { - Filename[SampleRun]: 'PG3_4844_event.zip', - Filename[VanadiumRun]: 'PG3_4866_event.zip', - CalibrationFilename: 'PG3_FERNS_d4832_2011_08_24.zip', + NeXusDetectorName: "powgen_detector", + Filename[SampleRun]: powgen.data.powgen_tutorial_sample_file(), + Filename[VanadiumRun]: powgen.data.powgen_tutorial_vanadium_file(), + CalibrationFilename: powgen.data.powgen_tutorial_calibration_file(), UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, - ValidTofRange: sc.array(dims=['tof'], values=[0.0, 16666.67], unit='us'), DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 200, unit='angstrom'), + TofMask: lambda x: (x < sc.scalar(0.0, unit="us")) + | (x > sc.scalar(16666.67, unit="us")), } @@ -42,9 +51,9 @@ def test_can_create_pipeline(providers, params): sciline.Pipeline(providers, params=params) -def test_pipeline_can_compute_dspacing_histogram(providers, params): +def test_pipeline_can_compute_dspacing_result(providers, params): pipeline = sciline.Pipeline(providers, params=params) - result = pipeline.compute(DspacingHistogram) + result = pipeline.compute(IofDspacing) assert result.sizes == { 'dspacing': len(params[DspacingBins]) - 1, } @@ -55,31 +64,64 @@ def test_workflow_is_deterministic(providers, params): pipeline = sciline.Pipeline(providers, params=params) # This is Sciline's default scheduler, but we want to be explicit here scheduler = sciline.scheduler.DaskScheduler() - graph = pipeline.get(DspacingHistogram, scheduler=scheduler) - reference = graph.compute().data - result = graph.compute().data + graph = pipeline.get(IofDspacing, scheduler=scheduler) + reference = graph.compute().hist().data + result = graph.compute().hist().data assert sc.identical(sc.values(result), sc.values(reference)) def test_pipeline_can_compute_intermediate_results(providers, params): pipeline = sciline.Pipeline(providers, params=params) result = pipeline.compute(NormalizedByProtonCharge[SampleRun]) - assert set(result.dims) == {'spectrum', 'tof'} + assert set(result.dims) == {'bank', 'column', 'row'} def test_pipeline_group_by_two_theta(providers, params): - from ess.powder.grouping import group_by_two_theta, merge_all_pixels - - providers.remove(merge_all_pixels) - providers.append(group_by_two_theta) params[TwoThetaBins] = sc.linspace( dim='two_theta', unit='deg', start=25.0, stop=90.0, num=16 - ) + ).to(unit='rad') pipeline = sciline.Pipeline(providers, params=params) - result = pipeline.compute(DspacingHistogram) + result = pipeline.compute(IofDspacingTwoTheta) assert result.sizes == { 'two_theta': 15, 'dspacing': len(params[DspacingBins]) - 1, } assert sc.identical(result.coords['dspacing'], params[DspacingBins]) - assert sc.allclose(result.coords['two_theta'].to(unit='deg'), params[TwoThetaBins]) + assert sc.allclose(result.coords['two_theta'], params[TwoThetaBins]) + + +def test_pipeline_wavelength_masking(providers, params): + wmin = sc.scalar(0.18, unit="angstrom") + wmax = sc.scalar(0.21, unit="angstrom") + params[WavelengthMask] = lambda x: (x > wmin) & (x < wmax) + pipeline = sciline.Pipeline(providers, params=params) + masked_sample = pipeline.compute(MaskedData[SampleRun]) + assert 'wavelength' in masked_sample.bins.masks + sum_in_masked_region = ( + masked_sample.bin(wavelength=sc.concat([wmin, wmax], dim='wavelength')) + .sum() + .data + ) + assert sc.allclose( + sum_in_masked_region, + sc.scalar(0.0, unit=sum_in_masked_region.unit), + ) + + +def test_pipeline_two_theta_masking(providers, params): + tmin = sc.scalar(0.8, unit="rad") + tmax = sc.scalar(1.0, unit="rad") + params[TwoThetaMask] = lambda x: (x > tmin) & (x < tmax) + pipeline = sciline.Pipeline(providers, params=params) + masked_sample = pipeline.compute(MaskedData[SampleRun]) + assert 'two_theta' in masked_sample.masks + sum_in_masked_region = ( + masked_sample.flatten(to='pixel') + .hist(two_theta=sc.concat([tmin, tmax], dim='two_theta')) + .sum() + .data + ) + assert sc.allclose( + sum_in_masked_region, + sc.scalar(0.0, unit=sum_in_masked_region.unit), + )