diff --git a/examples/bcdi_pipeline_example_xfel.ipynb b/examples/bcdi_pipeline_example_xfel.ipynb new file mode 100644 index 00000000..3f5adbbb --- /dev/null +++ b/examples/bcdi_pipeline_example_xfel.ipynb @@ -0,0 +1,673 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# **BCDI Pipeline** \n", + "### A Notebook to Run the `BcdiPipeline` Instance \n", + "\n", + "This notebook provides a structured workflow for running a **Bragg Coherent Diffraction Imaging (BCDI) pipeline**. \n", + "\n", + "The `BcdiPipeline` class handles the entire process, including: \n", + "- **Pre-processing** → Data preparation and corrections. \n", + "- **Phase retrieval** → Running PyNX algorithms to reconstruct the phase. \n", + "- **Post-processing** → Refining, analysing (get the strain!), and visualising results. \n", + "\n", + "You can provide **either**: \n", + "- A **YAML parameter file** for full automation. \n", + "- A **Python dictionary** for interactive control in this notebook. \n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "execution": { + "iopub.execute_input": "2026-05-29T14:27:33.509559Z", + "iopub.status.busy": "2026-05-29T14:27:33.509033Z", + "iopub.status.idle": "2026-05-29T14:27:34.305653Z", + "shell.execute_reply": "2026-05-29T14:27:34.304508Z", + "shell.execute_reply.started": "2026-05-29T14:27:33.509534Z" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "# import required packages\n", + "import os\n", + "\n", + "import cdiutils # core library for BCDI processing\n", + "import os, shutil\n", + "\n", + "# os.environ[\"PATH\"] = (\n", + "# \"/home/azakaria/.conda/envs/cdiutils_abd313/bin:\"\n", + "# + os.environ[\"PATH\"]\n", + "# )\n", + "\n", + "# print(shutil.which(\"nvcc\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **General Parameters**\n", + "Here, define the key parameters for **accessing and saving data** before running the pipeline. \n", + "- **These parameters must be set manually by the user** before execution. \n", + "- The output data will be saved in a structured directory format based on `sample_name` and `scan`. However, you can change the directory path if you like.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "execution": { + "iopub.execute_input": "2026-05-29T14:27:34.314583Z", + "iopub.status.busy": "2026-05-29T14:27:34.314412Z", + "iopub.status.idle": "2026-05-29T14:27:37.280075Z", + "shell.execute_reply": "2026-05-29T14:27:37.279286Z", + "shell.execute_reply.started": "2026-05-29T14:27:34.314560Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[INFO] BcdiPipeline initialised.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Dump directory already exists, results will be saved in:\n", + "/home/azakaria/analysis/testing/results//S221/.\n" + ] + } + ], + "source": [ + "# define the key parameters (must be filled in by the user)\n", + "beamline_setup: str = \"xfel\" # example: \"ID01\" (provide the beamline setup)\n", + "experiment_file_path: str = (\n", + " \"/gpfs/exfel/exp/MID/202601/p010402/scratch/\" # example: \"/path/to/experiment/file.h5\"\n", + ")\n", + "sample_name: str = \"\" # example: \"Sample_Pt\" (specify the sample name)\n", + "scan: int = 221 # example: 42 (specify the scan number)\n", + "hkl=[1, 1, 1]\n", + "# choose where to save the results (default: current working directory)\n", + "dump_dir = os.getcwd() + f\"/results/{sample_name}/S{scan}/\"\n", + "detector_name = \"agipd\"\n", + "# load the parameters and parse them into the BcdiPipeline class instance\n", + "params = cdiutils.pipeline.get_params_from_variables(dir(), globals())\n", + "bcdi_pipeline = cdiutils.BcdiPipeline(params=params)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Pre-Processing** \n", + "\n", + "If you need to update specific parameters, you can **pass them directly** into the `preprocess` method. \n", + "\n", + "### **Main Parameters**\n", + "- `preprocess_shape` → The shape of the cropped window used throughout the processes. \n", + " - Can be a **tuple of 2 or 3 values**. \n", + " - If only **2 values**, the entire rocking curve is used. \n", + "\n", + "- `voxel_reference_methods` → A `list` (or a single value) defining how to centre the data. \n", + " - Can include `\"com\"`, `\"max\"`, or a `tuple` of `int` (specific voxel position). \n", + " - Example:\n", + " ```python\n", + " voxel_reference_methods = [(70, 200, 200), \"com\", \"com\"]\n", + " ```\n", + " - This centres a box of size `preprocess_shape` around `(70, 200, 200)`, then iteratively refines it using `\"com\"` (only computed within this box).\n", + " - Useful when `\"com\"` fails due to artifacts or `\"max\"` fails due to hot pixels. \n", + " - Default: `[\"max\", \"com\", \"com\"]`. \n", + "\n", + "- `rocking_angle_binning` → If you want to bin in the **rocking curve direction**, provide a binning factor (ex.: `2`). \n", + "\n", + "- `light_loading` → If `True`, loads only the **ROI of the data** based on `voxel_reference_methods` and `preprocess_output_shape`. \n", + "\n", + "- `hot_pixel_filter` → Removes isolated hot pixels. \n", + " - Default: `False`. \n", + "\n", + "- `background_level` → Sets the background intensity to be removed. \n", + " - Example: `3`. \n", + " - Default: `None`. \n", + "\n", + "- `hkl` → Defines the **Bragg reflection** measured to extend *d*-spacing values to the lattice parameter. \n", + " - Default: `[1, 1, 1]`. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "execution_failed": "2026-05-29T14:27:25.901Z", + "iopub.execute_input": "2026-05-29T14:27:25.466360Z", + "iopub.status.busy": "2026-05-29T14:27:25.466142Z" + }, + "scrolled": true + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[INFO] \n", + "*******************************************************************************\n", + "* Starting process: preprocess *\n", + "*******************************************************************************\n", + "\n", + "[INFO] Additional parameters provided, will update the current dictionary of parameters.\n", + "[INFO] \n", + "PyNX needs a different shape to comply with gpu-based FFT, requested shape\n", + "(151, 200, 200) will be cropped to (150, 200, 200).\n", + "\n" + ] + } + ], + "source": [ + "%%time\n", + "bcdi_pipeline.preprocess(\n", + " preprocess_shape=(151, 200, 200), # define cropped window size\n", + " voxel_reference_methods=['max',], # centring method sequence\n", + " hot_pixel_filter=True, # remove isolated hot pixels\n", + " background_level=None, # background intensity level to remove\n", + " light_loading=False, # use light loading to save memory upon loading\n", + " # det_calib_params={\n", + " # \"distance\": 1.82169, # direct beam position vertical,\n", + " # \"cch1\": 1055.74, # horizontal\n", + " # \"cch2\": 1332.94, # detector pixel size in m\n", + " # \"pwidth1\": 7.5e-05, # detector pixel size in m\n", + " # \"pwidth2\": 7.5e-05, # sample to detector distance in m\n", + " # \"detrot\": 0,\n", + " # \"tilt\": 0,\n", + " # \"tiltazimuth\": 0,\n", + " # },\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **[PyNX](https://pynx.esrf.fr/en/latest/index.html) Phase Retrieval**\n", + "See the [pynx.cdi](https://pynx.esrf.fr/en/latest/scripts/pynx-cdi-id01.html) documentation for details on the phasing algorithms used here. \n", + "\n", + "**Algorithm recipe**\n", + "\n", + "You can either: \n", + "- provide the exact chain of algorithms. \n", + "- or specify the number of iterations for **RAAR**, **HIO**, and **ER**. \n", + "\n", + "```python\n", + "algorithm = None # ex: \"(Sup * (ER**20)) ** 10, (Sup*(HIO**20)) ** 15, (Sup*(RAAR**20)) ** 25\"\n", + "nb_raar = 500\n", + "nb_hio = 300\n", + "nb_er = 200\n", + "psf = \"pseudo-voigt,1,0.05,20\"\n", + "```\n", + "**Support-related parameters**\n", + "```python\n", + "support = \"auto\" # ex: bcdi_pipeline.pynx_phasing_dir + \"support.cxi\" (path to an existing support)\n", + "```\n", + ">_Note: If strain seems to large, don't use \"auto\" (autocorrelation) but use \"circle\" or \"square\", in combination with \"support_size\"_ \n", + "```python\n", + "support_threshold = \"0.15, 0.40\" # must be a string\n", + "support_update_period = 20\n", + "support_only_shrink = False\n", + "support_post_expand = None # ex: \"-1,1\" or \"-1,2,-1\"\n", + "support_update_border_n = None\n", + "support_smooth_width_begin = 2\n", + "support_smooth_width_end = 0.5\n", + "```\n", + "**Other parameters**\n", + "```python\n", + "positivity = False\n", + "beta = 0.9 # β parameter in HIO and RAAR\n", + "detwin = True\n", + "rebin = \"1, 1, 1\" # must be a string\n", + "```\n", + "**Number of Runs & Reconstructions to Keep**\n", + "```python\n", + "nb_run = 20 # total number of runs\n", + "nb_run_keep = 10 # number of reconstructions to keep\n", + "```\n", + "\n", + "**Override defaults in `phase_retrieval`**\n", + "\n", + "You can override any default parameter directly in the phase_retrieval method:\n", + "```python\n", + "bcdi_pipeline.phase_retrieval(nb_run=50, nb_run_keep=25)\n", + "```\n", + "If a parameter is not provided, the default value is used.\n", + "\n", + "### **Phase Retrieval GUI**\n", + "You can also launch a **Graphical User Interface (GUI)** to interactively set parameters and run phase retrieval. \n", + "```python\n", + "bcdi_pipeline.phase_retrieval_gui()\n", + "```\n", + "In that case, you can take care of the the result analysis, the selection of the best reconstructions, and the mode decomposition. Then, simply jump to the post-processing step cell." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "execution": { + "iopub.execute_input": "2026-05-29T14:28:16.281826Z", + "iopub.status.busy": "2026-05-29T14:28:16.281587Z", + "iopub.status.idle": "2026-05-29T14:28:19.016392Z", + "shell.execute_reply": "2026-05-29T14:28:19.015648Z", + "shell.execute_reply.started": "2026-05-29T14:28:16.281803Z" + }, + "scrolled": true + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[INFO] \n", + "*******************************************************************************\n", + "* Starting process: phase_retrieval *\n", + "*******************************************************************************\n", + "\n", + "[INFO] Removing former results.\n", + "\n", + "[INFO] Assuming the current machine is running PyNX. Will run the provided command.\n", + "[INFO] Using parameters file: pynx-cdi-inputs.txt\n", + "[INFO] Loading data: /home/azakaria/analysis/testing/results//S221//pynx_phasing/S221_pynx_input_data.npz\n", + "[INFO] Finished loading iobs data, with size: 6000000\n", + "[INFO] CDI runner: preparing processing unit\n", + "[INFO] Computing speed for available CUDA GPU [ranking by global memory bandwidth]:\n", + "[INFO] ERROR occurred on: max-p3ag012.desy.de [Linux-5.14.0-611.54.1.el9_7.x86_64-x86_64-with-glibc2.34]\n", + "[INFO] \n", + "[INFO] \n", + "[INFO] Caught exception for scan None:\n", + "[INFO] Failed initialising GPU. Please check GPU name or CUDA/OpenCL installation\n", + "[INFO] original error: error invoking 'nvcc --preprocess -arch sm_60 -I/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/cuda /tmp/tmpzzmmxe5h.cu --compiler-options -P': [Errno 2] No such file or directory: 'nvcc'\n", + "[INFO] on max-p3ag012.desy.de [Linux-5.14.0-611.54.1.el9_7.x86_64-x86_64-with-glibc2.34]\n", + "[INFO] \n", + "[INFO] Traceback (most recent call last):\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pytools/prefork.py\", line 101, in call_capture_output\n", + "[INFO] popen = Popen(cmdline, cwd=cwd, stdin=PIPE, stdout=PIPE,\n", + "[INFO] stderr=PIPE)\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/subprocess.py\", line 1039, in __init__\n", + "[INFO] self._execute_child(args, executable, preexec_fn, close_fds,\n", + "[INFO] ~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] pass_fds, cwd, env,\n", + "[INFO] ^^^^^^^^^^^^^^^^^^^\n", + "[INFO] ...<5 lines>...\n", + "[INFO] gid, gids, uid, umask,\n", + "[INFO] ^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] start_new_session, process_group)\n", + "[INFO] ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/subprocess.py\", line 1991, in _execute_child\n", + "[INFO] raise child_exception_type(errno_num, err_msg, err_filename)\n", + "[INFO] FileNotFoundError: [Errno 2] No such file or directory: 'nvcc'\n", + "[INFO] \n", + "[INFO] The above exception was the direct cause of the following exception:\n", + "[INFO] \n", + "[INFO] Traceback (most recent call last):\n", + "[INFO] File \"/home/azakaria/packages/PyNX/pynx/cdi/runner/runner.py\", line 1132, in prepare_processing_unit\n", + "[INFO] default_processing_unit.select_gpu(gpu_name=self.params['gpu'], verbose=True)\n", + "[INFO] ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] File \"/home/azakaria/packages/PyNX/pynx/processing_unit/cu_processing_unit.py\", line 77, in select_gpu\n", + "[INFO] super().select_gpu(gpu_name=gpu_name, gpu_rank=gpu_rank, ranking=ranking, language=None, verbose=verbose)\n", + "[INFO] ~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] File \"/home/azakaria/packages/PyNX/pynx/processing_unit/__init__.py\", line 328, in select_gpu\n", + "[INFO] benchmark_results += cuda_device.available_gpu_speed(fft_shape=fft_shape,\n", + "[INFO] ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] batch=self.benchmark_fft_batch, min_gpu_mem=None,\n", + "[INFO] ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] verbose=verbose, gpu_name=gpu_name,\n", + "[INFO] ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] ranking=ranking)\n", + "[INFO] ^^^^^^^^^^^^^^^^\n", + "[INFO] File \"/home/azakaria/packages/PyNX/pynx/processing_unit/cuda_device.py\", line 231, in available_gpu_speed\n", + "[INFO] gbps = int(round(cuda_device_global_mem_bandwidth(d, measured=True)))\n", + "[INFO] ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^\n", + "[INFO] File \"/home/azakaria/packages/PyNX/pynx/processing_unit/cuda_device.py\", line 136, in cuda_device_global_mem_bandwidth\n", + "[INFO] a = curandom.rand(shape=(512, 512, 512), dtype=np.float32)\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/curandom.py\", line 195, in rand\n", + "[INFO] func = get_elwise_kernel(\n", + "[INFO] \"float *dest, unsigned int seed\",\n", + "[INFO] ...<11 lines>...\n", + "[INFO] \"md5_rng_float\",\n", + "[INFO] )\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/elementwise.py\", line 193, in get_elwise_kernel\n", + "[INFO] _mod, func, arguments = get_elwise_kernel_and_types(\n", + "[INFO] ~~~~~~~~~~~~~~~~~~~~~~~~~~~^\n", + "[INFO] arguments, operation, name, keep, options, **kwargs\n", + "[INFO] ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] )\n", + "[INFO] ^\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/elementwise.py\", line 179, in get_elwise_kernel_and_types\n", + "[INFO] mod = module_builder(arguments, operation, name, keep, options, **kwargs)\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/elementwise.py\", line 51, in get_elwise_module\n", + "[INFO] return SourceModule(\n", + "[INFO] \"\"\"\n", + "[INFO] ...<33 lines>...\n", + "[INFO] no_extern_c=True,\n", + "[INFO] )\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/compiler.py\", line 348, in __init__\n", + "[INFO] cubin = compile(\n", + "[INFO] source,\n", + "[INFO] ...<7 lines>...\n", + "[INFO] include_dirs,\n", + "[INFO] )\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/compiler.py\", line 295, in compile\n", + "[INFO] return compile_plain(source, options, keep, nvcc, cache_dir, target)\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/compiler.py\", line 90, in compile_plain\n", + "[INFO] checksum.update(preprocess_source(source, options, nvcc).encode(\"utf-8\"))\n", + "[INFO] ~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/compiler.py\", line 53, in preprocess_source\n", + "[INFO] result, stdout, stderr = call_capture_output(cmdline, error_on_nonzero=False)\n", + "[INFO] ~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pytools/prefork.py\", line 306, in call_capture_output\n", + "[INFO] return forker.call_capture_output(cmdline, cwd, error_on_nonzero)\n", + "[INFO] ~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n", + "[INFO] File \"/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pytools/prefork.py\", line 113, in call_capture_output\n", + "[INFO] raise ExecError(\n", + "[INFO] \"error invoking '{}': {}\".format(\" \".join(cmdline), e)) from e\n", + "[INFO] pytools.prefork.ExecError: error invoking 'nvcc --preprocess -arch sm_60 -I/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/cuda /tmp/tmpzzmmxe5h.cu --compiler-options -P': [Errno 2] No such file or directory: 'nvcc'\n", + "[INFO] \n", + "[INFO] During handling of the above exception, another exception occurred:\n", + "[INFO] \n", + "[INFO] Traceback (most recent call last):\n", + "[INFO] File \"/home/azakaria/packages/PyNX/pynx/cdi/runner/runner.py\", line 2211, in process_scans\n", + "[INFO] self.ws.prepare_processing_unit()\n", + "[INFO] ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^\n", + "[INFO] File \"/home/azakaria/packages/PyNX/pynx/cdi/runner/runner.py\", line 1139, in prepare_processing_unit\n", + "[INFO] raise CDIRunnerException(\n", + "[INFO] \"Failed initialising GPU. Please check GPU name or CUDA/OpenCL installation\" + s0)\n", + "[INFO] pynx.cdi.runner.runner.CDIRunnerException: Failed initialising GPU. Please check GPU name or CUDA/OpenCL installation\n", + "[INFO] original error: error invoking 'nvcc --preprocess -arch sm_60 -I/home/azakaria/.conda/envs/cdiutils_abd313/lib/python3.13/site-packages/pycuda/cuda /tmp/tmpzzmmxe5h.cu --compiler-options -P': [Errno 2] No such file or directory: 'nvcc'\n", + "[INFO] \n", + "[ERROR] PyNX phasing process failed with return code 1\n", + "[ERROR] \n", + "Error occurred in the 'phase_retrieval' process:\n", + "PyNX script failed due to the following error:\n", + "Process exited with code 1 (no stderr output captured)\n" + ] + }, + { + "ename": "PyNXScriptError", + "evalue": "PyNX script failed due to the following error:\nProcess exited with code 1 (no stderr output captured)", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mPyNXScriptError\u001b[0m Traceback (most recent call last)", + "File \u001b[0;32m:1\u001b[0m\n", + "File \u001b[0;32m~/packages/cdiutils/src/cdiutils/pipeline/base.py:276\u001b[0m, in \u001b[0;36mPipeline.process..wrapper\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 273\u001b[0m sys\u001b[38;5;241m.\u001b[39mstdout \u001b[38;5;241m=\u001b[39m LoggerWriter(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlogger, logging\u001b[38;5;241m.\u001b[39mINFO)\n\u001b[1;32m 275\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 276\u001b[0m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 277\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlogger\u001b[38;5;241m.\u001b[39minfo(\n\u001b[1;32m 278\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mProcess \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mfunc\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m completed successfully.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 279\u001b[0m )\n\u001b[1;32m 280\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m e:\n", + "File \u001b[0;32m~/packages/cdiutils/src/cdiutils/pipeline/bcdi.py:1166\u001b[0m, in \u001b[0;36mBcdiPipeline.phase_retrieval\u001b[0;34m(self, jump_to_cluster, pynx_slurm_file_template, clear_former_results, cmd, search_pattern, **pynx_params)\u001b[0m\n\u001b[1;32m 1162\u001b[0m cmd \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpynx-cdi-id01 pynx-cdi-inputs.txt\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 1163\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlogger\u001b[38;5;241m.\u001b[39minfo(\n\u001b[1;32m 1164\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNo command provided. Will use the default: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mcmd\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 1165\u001b[0m )\n\u001b[0;32m-> 1166\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_run_cmd\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcmd\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpynx_phasing_dir\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/packages/cdiutils/src/cdiutils/pipeline/bcdi.py:1240\u001b[0m, in \u001b[0;36mBcdiPipeline._run_cmd\u001b[0;34m(self, cmd, cwd)\u001b[0m\n\u001b[1;32m 1228\u001b[0m stderr_text \u001b[38;5;241m=\u001b[39m (\n\u001b[1;32m 1229\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mjoin(stderr_lines)\n\u001b[1;32m 1230\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m stderr_lines\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 1234\u001b[0m )\n\u001b[1;32m 1235\u001b[0m )\n\u001b[1;32m 1236\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlogger\u001b[38;5;241m.\u001b[39merror(\n\u001b[1;32m 1237\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPyNX phasing process failed with return code \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 1238\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mproc\u001b[38;5;241m.\u001b[39mreturncode\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 1239\u001b[0m )\n\u001b[0;32m-> 1240\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m PyNXScriptError(stderr_text)\n\u001b[1;32m 1241\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m subprocess\u001b[38;5;241m.\u001b[39mCalledProcessError \u001b[38;5;28;01mas\u001b[39;00m e:\n\u001b[1;32m 1242\u001b[0m \u001b[38;5;66;03m# log the error if the job submission fails\u001b[39;00m\n\u001b[1;32m 1243\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlogger\u001b[38;5;241m.\u001b[39merror(\n\u001b[1;32m 1244\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSubprocess process failed with return code \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 1245\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00me\u001b[38;5;241m.\u001b[39mreturncode\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00me\u001b[38;5;241m.\u001b[39mstderr\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 1246\u001b[0m )\n", + "\u001b[0;31mPyNXScriptError\u001b[0m: PyNX script failed due to the following error:\nProcess exited with code 1 (no stderr output captured)" + ] + } + ], + "source": [ + "%%time\n", + "bcdi_pipeline.phase_retrieval(\n", + " cmd=r'~/.conda/envs/cdiutils_abd313/bin/pynx-cdi-id01 pynx-cdi-inputs.txt',\n", + " clear_former_results=True,\n", + " nb_run=20,\n", + " nb_run_keep=10,\n", + " # support=bcdi_pipeline.pynx_phasing_dir + \"support.cxi\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Analyse the phasing results**\n", + "\n", + "This step evaluates the quality of the phase retrieval results by sorting reconstructions based on a `sorting_criterion`. \n", + "\n", + "##### **Available Sorting Criteria**\n", + "- `\"mean_to_max\"` → Difference between the mean of the **Gaussian fit of the amplitude histogram** and its maximum value. A **smaller difference** indicates a more homogeneous reconstruction. \n", + "- `\"sharpness\"` → Sum of the amplitude within the support raised to the power of 4. **Lower values** indicate greater homogeneity. \n", + "- `\"std\"` → **Standard deviation** of the amplitude. \n", + "- `\"llk\"` → **Log-likelihood** of the reconstruction. \n", + "- `\"llkf\"` → **Free log-likelihood** of the reconstruction. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "execution_failed": "2026-05-29T14:27:25.901Z" + } + }, + "outputs": [], + "source": [ + "%%time\n", + "bcdi_pipeline.analyse_phasing_results(\n", + " sorting_criterion=\"mean_to_max\", # selects the sorting method\n", + " # Optional parameters\n", + " # plot_phasing_results=False, # uncomment to disable plotting\n", + " # plot_phase=True, # uncomment to enable phase plotting\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Optionally, generate a support for further phasing attempts** \n", + "\n", + "##### **Parameters**\n", + "- `run` → set to either: \n", + " - `\"best\"` to use the best reconstruction. \n", + " - an **integer** corresponding to the specific run you want. \n", + "- `output_path` → the location to save the generated support. By default, it will be saved in the `pynx_phasing` folder. \n", + "- `fill` → whether to fill the support if it contains holes. \n", + " - Default: `False`.\n", + "- `verbose` → whether to print logs and display a plot of the support. \n", + "\n", + "" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "execution_failed": "2026-05-29T14:27:25.901Z" + } + }, + "outputs": [], + "source": [ + "# bcdi_pipeline.generate_support_from(\"best\", fill=False) # uncomment to generate a support" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Selection of the best reconstructions & mode decomposition**\n", + "\n", + "You can select the best reconstructions based on a **sorting criterion** and keep a specified number of top candidates. \n", + "\n", + "##### **Parameters**\n", + "- `nb_of_best_sorted_runs` → the number of best reconstructions to keep, selected based on the `sorting_criterion` used in the `analyse_phasing_results` method above. \n", + "- `best_runs` → instead of selecting based on sorting, you can manually specify a list of reconstruction numbers.\n", + "\n", + "By default, the **best reconstructions** are automatically selected. \n", + "\n", + "Once the best candidates are chosen, `mode_decomposition` analyses them to extract dominant features. \n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "execution_failed": "2026-05-29T14:27:25.901Z" + } + }, + "outputs": [], + "source": [ + "%%time\n", + "# define how many of the best candidates to keep\n", + "number_of_best_candidates: int = 5\n", + "\n", + "# select the best reconstructions based on the sorting criterion\n", + "bcdi_pipeline.select_best_candidates(\n", + " nb_of_best_sorted_runs=number_of_best_candidates\n", + " # best_runs=[10] # uncomment to manually select a specific run\n", + ")\n", + "\n", + "# perform mode decomposition on the selected reconstructions\n", + "bcdi_pipeline.mode_decomposition()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Post-processing**\n", + "\n", + "This stage includes several key operations: \n", + "- **orthogonalisation** of the reconstructed data. \n", + "- **phase manipulation**: \n", + " - phase unwrapping \n", + " - phase ramp removal \n", + "- **computation of physical properties**: \n", + " - displacement field \n", + " - strain \n", + " - d-spacing \n", + "- **visualisation**: Generate multiple plots for analysis. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "execution_failed": "2026-05-29T14:27:25.902Z" + } + }, + "outputs": [], + "source": [ + "%%time\n", + "bcdi_pipeline.postprocess(\n", + " isosurface=0.5, # threshold for isosurface\n", + " voxel_size=None, # use default voxel size if not provided\n", + " flip=True, # whether to flip the reconstruction if you got the twin image (enantiomorph)\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **3D interactive plot**\n", + "\n", + "Display an interactive 3D isosurface of the final reconstruction and explore different quantities for colouring.\n", + "\n", + "What you can do\n", + "- Visualise an isosurface of the reconstructed object (amplitude / support).\n", + "- Colour the surface by different quantities: amplitude, phase, displacement, strain, d-spacing, etc.\n", + "- Interactively adjust:\n", + " - isosurface threshold (isosurface level)\n", + " - colormap and value range\n", + "- Rotate, zoom and pan the scene with the mouse; use the toolbar to reset view or save a screenshot.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "execution_failed": "2026-05-29T14:27:25.902Z" + } + }, + "outputs": [], + "source": [ + "bcdi_pipeline.show_3d_final_result()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Feedback & Issue Reporting** \n", + "\n", + "If you have **comments, suggestions, or encounter any issues**, please reach out: \n", + "\n", + "📧 **Email:** [clement.atlan@esrf.fr](mailto:clement.atlan@esrf.fr?subject=cdiutils) \n", + "🐙 **GitHub Issues:** [Report an issue](https://github.com/clatlan/cdiutils/issues) \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Credits\n", + "This notebook was created by Clément Atlan, ESRF, 2025. It is part of the `cdiutils` package, which provides tools for BCDI data analysis and visualisation.\n", + "If you have used this notebook or the `cdiutils` package in your research, please consider citing the package https://github.com/clatlan/cdiutils/\n", + "You'll find the citation information in the `cdiutils` package documentation.\n", + "\n", + "```bibtex\n", + "@software{Atlan_Cdiutils_A_python,\n", + "author = {Atlan, Clement},\n", + "doi = {10.5281/zenodo.7656853},\n", + "license = {MIT},\n", + "title = {{Cdiutils: A python package for Bragg Coherent Diffraction Imaging processing, analysis and visualisation workflows}},\n", + "url = {https://github.com/clatlan/cdiutils},\n", + "version = {0.2.0}\n", + "}\n", + "```\n", + "\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (cdiutils_abd313)", + "language": "python", + "name": "cdiutils_abd313" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.13" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/pyproject.toml b/pyproject.toml index 9f0407d3..797307c4 100755 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,6 +43,8 @@ dependencies = [ "tabulate", "xrayutilities>=1.7.3", "pyyaml", + "damnit", + "EXtra-data", ] [project.optional-dependencies] diff --git a/src/cdiutils/geometry.py b/src/cdiutils/geometry.py index 4d7e629c..701f87e3 100644 --- a/src/cdiutils/geometry.py +++ b/src/cdiutils/geometry.py @@ -287,6 +287,17 @@ def from_setup( sample_surface_normal=[0, 1, 0], # default sample facing up name="ID27", ) + if beamline.lower() == "xfel": + geometry = cls( + sample_circles=["y-", "x+", "z+"], # theta, chi, phi + detector_circles=["y-"], # twotheta + detector_axis0_orientation="y-", + detector_axis1_orientation="z-", + beam_direction=[1, 0, 0], + sample_surface_normal=[0, 0, 1], + name="XFEL", + is_cxi=False, + ) if geometry is None: raise NotImplementedError( f"The beamline name {beamline} is not valid. Available:\n" diff --git a/src/cdiutils/io/__init__.py b/src/cdiutils/io/__init__.py index 569c9ae3..37203f36 100755 --- a/src/cdiutils/io/__init__.py +++ b/src/cdiutils/io/__init__.py @@ -14,6 +14,7 @@ from .p10 import P10Loader from .sixs import SIXSLoader from .vtk import save_as_vti +from .xfel import XFELLoader, load_xfel __all__ = [ "Loader", @@ -25,6 +26,8 @@ "SIXSLoader", "CristalLoader", "NanoMAXLoader", + "XFELLoader", + "load_xfel", "CXIFile", "CXIExplorer", "save_as_cxi", diff --git a/src/cdiutils/io/loader.py b/src/cdiutils/io/loader.py index 0d6fb99d..5e9407db 100644 --- a/src/cdiutils/io/loader.py +++ b/src/cdiutils/io/loader.py @@ -212,7 +212,10 @@ def from_setup(cls, beamline_setup: str, **metadata) -> "Loader": "Only 2019 and 2022 versions are available for now. Specify " "the version in the beamline_setup." ) + if beamline_setup.lower() == "xfel": + from . import XFELLoader + return XFELLoader(**metadata) if "p10" in beamline_setup.lower(): from . import P10Loader @@ -802,6 +805,8 @@ def get_mask( mask = np.zeros(shape=(512, 1028)) elif detector_name.lower() == "merlin": mask = np.zeros(shape=(512, 512)) + elif detector_name.lower() == "agipd": + mask = np.zeros(shape=(1325, 1196)) else: raise ValueError(f"Invalid detector name: {detector_name}") if channel: diff --git a/src/cdiutils/io/xfel.py b/src/cdiutils/io/xfel.py new file mode 100644 index 00000000..74261f6d --- /dev/null +++ b/src/cdiutils/io/xfel.py @@ -0,0 +1,496 @@ +"""European XFEL I/O utilities.""" + +# TODO: +# Verify mapping between XFEL theta/chi/phi/twotheta +# and cdiutils sample/detector angle conventions. +from __future__ import annotations + +from pathlib import Path + +import numpy as np + +from cdiutils.io.loader import H5TypeLoader + +try: + from extra_data import RunDirectory +except ImportError as exc: + raise ImportError( + "XFELLoader requires EXtra-data. " + "Install/use it in the European XFEL analysis environment." + ) from exc +try: + from damnit import Damnit +except ImportError as exc: + raise ImportError( + "XFELLoader requires DAMNIT to read processed XFEL variables." + ) from exc + + +def xfel_safe_load(func): + """Safe loader wrapper for XFEL loading.""" + + def wrap(self, *args, **kwargs): + if not self.experiment_file_path.exists(): + raise FileNotFoundError( + f"XFEL experiment path does not exist: {self.experiment_file_path}" + ) + + if not self.experiment_file_path.is_dir(): + raise NotADirectoryError( + f"XFEL experiment path should be a directory: {self.experiment_file_path}" + ) + + return func(self, *args, **kwargs) + + return wrap + + +class XFELLoader(H5TypeLoader): + """Loader for European XFEL DAMNIT/EXtra-data processed datasets.""" + + angle_names = { + "sample_outofplane_angle": "mu", + "sample_inplane_angle": "omega", + "detector_outofplane_angle": "gamma", + "detector_inplane_angle": "delta", + } + + authorised_detector_names = ("agipd", "jungfrau", "epix") + + def __init__( + self, + experiment_file_path: str, + scan: int = None, + sample_name: str = None, + detector_name: str = None, + flat_field: np.ndarray | str = None, + alien_mask: np.ndarray | str = None, + aliases_file_name: str = "extra-data-aliases.yml", + data_key: str = "peak_images", + pulse_dimension: str = "pulseIndex", + pulse_reduction: str = "mean", + **kwargs, + ) -> None: + """ + Initialise an XFEL loader. + + Args: + experiment_file_path: Path to the XFEL scratch directory. + scan: DAMNIT run/scan number. + sample_name: Optional sample name. + detector_name: Detector name, e.g. "agipd". + flat_field: Optional flat-field correction. + alien_mask: Optional detector mask. + aliases_file_name: Alias file name inside the run directory. + data_key: DAMNIT variable containing detector images. + pulse_dimension: Xarray pulse dimension name. + pulse_reduction: Reduction over pulse dimension: + "mean", "sum", or None. + """ + self.aliases_file_name = aliases_file_name + self.data_key = data_key + self.pulse_dimension = pulse_dimension + self.pulse_reduction = pulse_reduction + + super().__init__( + experiment_file_path, + scan, + sample_name, + detector_name, + flat_field, + alien_mask, + ) + + self.experiment_file_path = Path(self.experiment_file_path) + + def _get_run(self): + """Return the EXtra-data run object.""" + run_dir_name = self.sample_name or self.run_dir_name + run_dir = self.experiment_file_path.parent / run_dir_name + aliases_file = run_dir / self.aliases_file_name + + run = RunDirectory(run_dir) + if aliases_file.exists(): + run = run.with_aliases(aliases_file) + + return run + + def _get_run_vars(self, scan: int = None): + """Return DAMNIT variables for one XFEL run/scan.""" + scan, _ = self._check_scan_sample(scan, None) + + damnit_path = self.experiment_file_path + db = Damnit(damnit_path) + + return db[scan] + + def _read_images(self, scan: int = None): + """Read the detector image variable from DAMNIT.""" + run_vars = self._get_run_vars(scan) + available_keys = list(run_vars.keys()) + + try: + return run_vars[self.data_key].read() + + except KeyError as exc: + + # Preferred fallback: xarray dataset containing detector images + if "peak_dataset" in available_keys: + ds = run_vars["peak_dataset"].read() + + if "images" in ds: + print( + f"\nDAMNIT variable {self.data_key!r} was not found." + "\nUsing fallback dataset: 'peak_dataset/images'" + ) + return ds["images"] + + candidates = [] + + for key in available_keys: + try: + data = run_vars[key].read() + shape = np.shape(data) + + if len(shape) in (3, 4): + candidates.append((key, shape)) + + except Exception: + continue + + if candidates: + print( + f"\nDAMNIT variable {self.data_key!r} was not found.\n" + "Available 3D/4D datasets:" + ) + + for key, shape in candidates: + print(f" {key:30s} shape={shape}") + + raise KeyError( + f"DAMNIT variable {self.data_key!r} was not found " + f"for scan/run {scan}." + ) from exc + def _reduce_pulses(self, images, pulse_reduction: str = None): + """Reduce the XFEL pulse dimension if present.""" + reduction = ( + self.pulse_reduction + if pulse_reduction is None + else pulse_reduction + ) + + if reduction is None: + return images + + if not hasattr(images, "dims"): + return images + + pulse_dim = self.pulse_dimension + + if pulse_dim not in images.dims: + if "pulseId" in images.dims: + pulse_dim = "pulseId" + elif "pulseIndex" in images.dims: + pulse_dim = "pulseIndex" + else: + return images + + if reduction == "mean": + return images.mean(pulse_dim) + + if reduction == "sum": + return images.sum(pulse_dim) + + raise ValueError("pulse_reduction should be 'mean', 'sum', or None.") + + @xfel_safe_load + def load_detector_data( + self, + scan: int = None, + sample_name: str = None, + roi: tuple[slice] = None, + rocking_angle_binning: int = None, + binning_method: str = "sum", + pulse_reduction: str = None, + ) -> np.ndarray: + """ + Load XFEL detector data as a cdiutils-compatible 3D array. + + Returns: + np.ndarray with shape: + (rocking_position, detector_y, detector_x) + """ + scan, sample_name = self._check_scan_sample(scan, sample_name) + + images = self._read_images(scan) + images = self._reduce_pulses(images, pulse_reduction) + + data = np.asarray(images) + data = np.nan_to_num(data, nan=0.0, posinf=0.0, neginf=0.0) + + roi = self._check_roi(roi) + data = data[roi] + + return self.bin_flat_mask( + data, + roi, + self.flat_field, + self.alien_mask, + rocking_angle_binning, + binning_method, + ) + + def _detect_scan_positions( + self, + motor_values, + n_steps=None, + ): + """Detect scan-step positions from train-resolved motor values.""" + values = np.asarray(motor_values, dtype=float) + values = values[np.isfinite(values)] + + if values.size == 0: + return values + + if n_steps is None: + raise ValueError( + "n_steps is required for XFEL scan-position detection." + ) + + chunks = np.array_split(values, n_steps) + + return np.asarray( + [np.nanmean(chunk) for chunk in chunks if chunk.size > 0] + ) + + @xfel_safe_load + def load_motor_positions( + self, + scan: int = None, + sample_name: str = None, + roi: tuple[slice] = None, + rocking_angle_binning: int = None, + theta=None, + chi=None, + phi=None, + twotheta=None, + theta_offset: float = 0.0, + chi_offset: float = 0.0, + phi_offset: float = 0.0, + twotheta_offset: float = 0.0, + ) -> dict: + """Load XFEL diffractometer angles.""" + scan, sample_name = self._check_scan_sample(scan, sample_name) + + angles, scanned_motor_name = self._load_geometry_angles( + theta=theta, + chi=chi, + phi=phi, + twotheta=twotheta, + theta_offset=theta_offset, + chi_offset=chi_offset, + phi_offset=phi_offset, + twotheta_offset=twotheta_offset, + ) + + self.rocking_angle = { + "theta": "sample_outofplane_angle", + "chi": "sample_inplane_angle", + "phi": "sample_inplane_angle", + "twotheta": "detector_inplane_angle", + }[scanned_motor_name] + + # Apply ROI/binning only to the scanned array angle. + scanned_angle = angles[scanned_motor_name] + + if rocking_angle_binning and np.ndim(scanned_angle) > 0: + scanned_angle = self.bin_rocking_angle_values( + scanned_angle, + rocking_angle_binning, + ) + + if roi is not None and np.ndim(scanned_angle) > 0: + if isinstance(roi, (tuple, list)) and len(roi) == 3: + scanned_angle = scanned_angle[roi[0]] + elif isinstance(roi, slice): + scanned_angle = scanned_angle[roi] + + angles[scanned_motor_name] = scanned_angle + + return { + "sample_outofplane_angle": angles["theta"], + "sample_inplane_angle": angles["chi"], + "detector_outofplane_angle": angles["phi"], + "detector_inplane_angle": angles["twotheta"], + } + + @xfel_safe_load + def load_det_calib_params(self) -> dict: + """Load XFEL detector calibration parameters.""" + run = self._get_run() + sdd = np.nanmedian(run.alias["sdd"].ndarray()) + + detector_shape = self.load_detector_shape() + cch1 = detector_shape[0] // 2 + cch2 = detector_shape[1] // 2 + + return { + "distance": float(sdd), # mm + "pwidth1": 0.2, # mm + "pwidth2": 0.2, # mm + "cch1": cch1, + "cch2": cch2, + "detrot": 0, + "tilt": 0, + "tiltazimuth": 0, + } + + @xfel_safe_load + def load_energy(self, scan: int = None) -> float: + """Load photon energy in eV.""" + run = self._get_run() + + # Direct energy aliases + for alias in ("energy-kev", "undulator-energy"): + if alias in run._aliases: + energy = run.alias[alias].ndarray() + energy = float(np.nanmean(energy)) + + # If value looks like keV, convert to eV + if energy < 100: + energy *= 1e3 + + return energy + + # Wavelength fallback + if "xgm-wavelength" in run._aliases: + wavelength = run.alias["xgm-wavelength"].ndarray() + wavelength = float(np.nanmean(wavelength)) + + # Usually XFEL wavelength may be in m; convert m -> Å + if wavelength < 1e-6: + wavelength *= 1e10 + + # E[eV] = hc / lambda[Å] + return 12398.419843320026 / wavelength + + raise RuntimeError( + "Could not load photon energy. Tried aliases: " + "'energy-kev', 'undulator-energy', 'xgm-wavelength'.\n" + f"Available aliases are:\n{sorted(run._aliases.keys())}" + ) + + @xfel_safe_load + def load_detector_shape(self, scan: int = None) -> tuple: + """Return detector image shape after pulse reduction.""" + data = self.load_detector_data(scan=scan) + + if data.ndim != 3: + raise ValueError( + f"Expected detector data with 3 dimensions, got shape {data.shape}." + ) + + return data.shape[1:] + + def _infer_scanned_motor_name(self, run, angles): + """Infer scanned motor from theta/chi/phi/twotheta aliases.""" + candidates = {} + + for name in ("theta", "chi", "phi", "twotheta"): + if angles[name] is not None: + continue + + values = np.asarray(run.alias[name].ndarray(), dtype=float) + diffs = np.diff(values) + movement = np.nanmax(values) - np.nanmin(values) + nonzero_steps = np.count_nonzero(np.abs(diffs) > 0) + + candidates[name] = (movement, nonzero_steps) + if not candidates: + raise ValueError( + "Cannot infer scanned motor because all angles " + "were provided explicitly." + ) + scanned_motor_name = max( + candidates, key=lambda key: candidates[key][0] + ) + + if candidates[scanned_motor_name][0] == 0: + raise RuntimeError( + "Could not infer scanned motor: theta, chi, phi, and twotheta " + "appear constant." + ) + + return scanned_motor_name + + def _load_geometry_angles( + self, + theta=None, + chi=None, + phi=None, + twotheta=None, + theta_offset=0.0, + chi_offset=0.0, + phi_offset=0.0, + twotheta_offset=0.0, + ): + """Load theta, chi, phi, and twotheta from EXtra-data aliases.""" + + run = self._get_run() + + angles = { + "theta": theta, + "chi": chi, + "phi": phi, + "twotheta": twotheta, + } + + images = self._read_images() + + if "position" in images.coords and angles["theta"] is None: + scanned_motor_name = "theta" + angles["theta"] = np.asarray(images.coords["position"].values) + else: + scanned_motor_name = self._infer_scanned_motor_name(run, angles) + + if angles[scanned_motor_name] is None: + motor_values = run.alias[scanned_motor_name].ndarray() + n_steps = np.asarray(images).shape[0] + angles[scanned_motor_name] = self._detect_scan_positions( + motor_values, + n_steps=n_steps, + ) + + for key, value in angles.items(): + if value is None: + angles[key] = run.alias[key].as_single_value() + + offsets = { + "theta": theta_offset, + "chi": chi_offset, + "phi": phi_offset, + "twotheta": twotheta_offset, + } + + for key in angles: + angles[key] = np.asarray(angles[key]) + offsets[key] + + return angles, scanned_motor_name + + +def load_xfel( + experiment_file_path: str, + scan: int = None, + sample_name: str = None, + detector_name: str = "agipd", + **kwargs, +) -> np.ndarray: + """Load XFEL detector data using XFELLoader.""" + loader = XFELLoader( + experiment_file_path=experiment_file_path, + scan=scan, + sample_name=sample_name, + detector_name=detector_name, + **kwargs, + ) + return loader.load_detector_data() diff --git a/src/cdiutils/pipeline/bcdi.py b/src/cdiutils/pipeline/bcdi.py index 870f77d2..15e619d5 100755 --- a/src/cdiutils/pipeline/bcdi.py +++ b/src/cdiutils/pipeline/bcdi.py @@ -623,6 +623,14 @@ def _load(self, roi: tuple[slice] = None) -> None: detector_name=self.params["detector_name"], roi=(slice(None), roi[1], roi[2]) if roi else None, ) + + if self.mask.shape != self.detector_data.shape: + self.logger.warning( + "Loaded detector mask shape does not match detector data shape. " + f"mask={self.mask.shape}, data={self.detector_data.shape}. " + "Using an empty mask instead." + ) + self.mask = np.zeros_like(self.detector_data, dtype=np.uint8) if loader.get_alien_mask() is not None: self.logger.info("Alien mask provided. Will update detector mask.") alien_mask = loader.get_alien_mask(roi)