{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bedd5844-1c5c-4fcf-a311-e0625b5e23b7",
   "metadata": {},
   "source": [
    "# Data reduction for POWGEN"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "094942ff-1bde-46c2-ae50-d00177eca3ad",
   "metadata": {},
   "source": [
    "This notebook shows a basic reduction workflow for powder diffraction for the SNS [POWGEN](https://sns.gov/powgen) instrument.\n",
    "It serves mainly to develop and present routines for powder diffraction and will eventually be removed in favor of a workflow for DREAM at ESS.\n",
    "\n",
    "**Note** that we load functions from `external` modules.\n",
    "These modules will be removed when their ESS counterparts exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b58104c6-a196-4576-b3f0-9fb6fb1216e9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import scipp as sc\n",
    "import scippneutron as scn\n",
    "\n",
    "import ess\n",
    "from ess.diffraction.external import load_calibration\n",
    "from ess.diffraction import powder\n",
    "from ess import diffraction\n",
    "from ess.external import powgen"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05194976-6bd6-4d2b-bfa2-c21105608bd2",
   "metadata": {},
   "source": [
    "Initialize the logging system to write logs to a file `powgen.log`.\n",
    "This also displays a widget which shows log messages emitted by the functions that we call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c88d48df-272a-4f83-8b4a-43450a365b1f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ess.logging.configure_workflow('powgen_reduction', filename='powgen.log')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da4040e8-681e-4eb7-a205-af6a88539ffa",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8db8e21f-f19c-4498-a4d1-d34142b2ba02",
   "metadata": {},
   "source": [
    "Load the sample data.\n",
    "We use the dedicated `load` function from `powgen` instead of `scippneutron.load` because the former passes the correct parameters to the underlying algorithm automatically.\n",
    "\n",
    "**Note:** We get the file name from `powgen.data`.\n",
    "This module provides access to managed example files.\n",
    "In the real world, we would need to find the file name in a different way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1d52d95-727c-4d67-b1f9-dc8c4e345d15",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample = powgen.load(powgen.data.sample_file())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "723704af-1dc1-4cce-a396-77a754b19d77",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bab068d8-d02c-4bec-a9a3-2c3b713ab750",
   "metadata": {},
   "source": [
    "## Inspect the raw data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "330aaf4b-465f-47e9-bede-45183971f9f6",
   "metadata": {},
   "source": [
    "We can plot the data array to get an idea of its contents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3848a30b-b872-465c-a1ca-2fad3a5110cb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59585c24-349f-4c5b-ab35-b17772712688",
   "metadata": {},
   "source": [
    "We can see how that data maps onto the detector by using POWGEN's instrument view."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e5a4d3c-cb93-4875-b2be-9e73df22771d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "powgen.instrument_view(sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07912224-f459-4c42-95e0-58d99291a713",
   "metadata": {},
   "source": [
    "## Filter out invalid events"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5121997-9522-4cd7-bbc4-468b38001aee",
   "metadata": {},
   "source": [
    "The file contains events that cannot be associated with a specific pulse.\n",
    "We can get a range of valid time-of-flight values from the instrument characterization file associated with the current run.\n",
    "There is currently no mechanism in `scippneutron` or `ess` to load such a file as it is not clear if ESS will use this approach.\n",
    "The values used below are taken from `PG3_characterization_2011_08_31-HR.txt` which is part of the sample files of Mantid.\n",
    "See, e.g., [PowderDiffractionReduction](https://www.mantidproject.org/PowderDiffractionReduction).\n",
    "\n",
    "We remove all events that have a time-of-flight value outside the valid range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ced73191-082a-4c7c-aa58-6111cd3cf1fd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample = sc.bin(sample, edges=[\n",
    "    sc.array(dims=['tof'], values=[0.0, 16666.67], unit='us')\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2de9ce54-f972-474a-ae82-c9754b80719f",
   "metadata": {},
   "source": [
    "## Normalize by proton charge"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebfdcccc-7638-4dbb-94f4-cf070d9ed28d",
   "metadata": {},
   "source": [
    "Next, we normalize the data by the proton charge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85bd7848-a7db-42c6-974a-1eb5bf3e9860",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample /= sample.attrs['gd_prtn_chrg']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0e9bd7b-897f-41a1-80f7-1a12c8424c82",
   "metadata": {},
   "source": [
    "We can check the unit of the event weight to see that the data was indeed divided by a charge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "063b72ab-7b38-47fa-a180-e4cb659efb1f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample.data.values[0].unit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52e9d8d7-474c-423b-baf4-e9e135f6e46d",
   "metadata": {},
   "source": [
    "## Compute d-spacing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9291753-2b26-4d5d-8396-f4d3a2a9abfc",
   "metadata": {},
   "source": [
    "Here, we compute d-spacing using calibration parameters provided in an example file.\n",
    "\n",
    "First, we load the calibration parameters.\n",
    "We need to specify an instrument definition that Mantid understands.\n",
    "It needs to be selected for the specific data that we are processing.\n",
    "\n",
    "**Note:** ESS instruments will not use instrument definition files (IDF).\n",
    "Instead, the instrument parameters will be encoded in the `ess` module and the instrument can be identified using data in the NeXus file.\n",
    "But for the purposes of this example, we need to provide the IDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ca5395a-4a43-4eaa-85fa-f43e6ac1a576",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cal = load_calibration(\n",
    "    powgen.data.calibration_file(),\n",
    "    instrument_filename='POWGEN_Definition_2011-02-25.xml',\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02c9608f-7b1d-412b-a4f0-df064b340916",
   "metadata": {},
   "source": [
    "The calibration is loaded with a 'detector' dimension.\n",
    "Compute the corresponding spectrum indices using the detector info loaded as part of the sample data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9b964ee-43b3-4bc1-b079-dfede08d3208",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cal = powgen.beamline.map_detector_to_spectrum(\n",
    "    cal, detector_info=sample.coords['detector_info'].value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59566cfb-e175-44ea-bd09-bd2324514405",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63553acb-4f78-4518-b5f8-6dc71dc616a5",
   "metadata": {},
   "source": [
    "Now when can compute d-spacing for the sample using the calibration parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69388346-24ab-4035-b18d-b10bf274e4c4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample_dspacing = powder.to_dspacing_with_calibration(sample, calibration=cal)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc226f10-71fe-4571-8391-db35028c4c76",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Vanadium correction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f670ce66-38cf-44ff-ad7c-d5dd30e3b09e",
   "metadata": {},
   "source": [
    "Before we can process the d-spacing distribution further, we need to normalize the data by a vanadium measurement.\n",
    "`ess.diffraction` provides a helper function to load event data for vanadium and process it in a similar way to what we did above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9494f15-1fb2-4236-95a7-800a04896e2f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vana = powgen.load_and_preprocess_vanadium(\n",
    "    powgen.data.vanadium_file(),\n",
    "    powgen.data.empty_instrument_file())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19cd62a0-94ef-4dd4-b00e-7d3d9d6c9d50",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vana"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8f39532-e0ee-48cc-8c49-40ac7427ded8",
   "metadata": {},
   "source": [
    "Currently, `load_and_preprocess_vanadium` uses a crude, preliminary event filtering mechanism which removed all events if the proton charge is too low.\n",
    "This happens with the empty instrument data in this case as can be seen in the log widget at the top of the notebook.\n",
    "\n",
    "Now, we compute d-spacing using the same calibration parameters as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03105e85-f66c-45b4-9b50-adb9192763a1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vana_dspacing = powder.to_dspacing_with_calibration(vana, calibration=cal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70fd1e76-a16d-4b6d-94d0-fb1284cb26a4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vana_dspacing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0441b09d-9b7e-45e6-a154-a822a21aebb2",
   "metadata": {},
   "source": [
    "## Inspect d-spacing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd2ef840-a5aa-45dd-aaea-c426dea00a3d",
   "metadata": {},
   "source": [
    "We need to histogram the events in order to normalize our sample data by vanadium.\n",
    "For consistency, we use these bin edges for both vanadium and the sample data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba9ab7db-5752-4130-8cfb-ec00c70feed8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "d = vana_dspacing.coords['dspacing']\n",
    "dspacing_edges = sc.linspace('dspacing',\n",
    "                             d.min().value,\n",
    "                             d.max().value,\n",
    "                             200,\n",
    "                             unit=d.unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2278a3db-6b79-4421-a1ff-2628226f4166",
   "metadata": {},
   "source": [
    "### All spectra combined"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84aac5f6-38bd-4246-974b-2a94c2dd80d0",
   "metadata": {},
   "source": [
    "We start simple by combining all spectra using `data.bins.concat('spectrum')`.\n",
    "Then, we can normalize the same data by vanadium to get a d-spacing distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bd287a8-d26e-4134-b1b5-3ed686cca417",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "all_spectra = diffraction.normalize_by_vanadium(\n",
    "    sample_dspacing.bins.concat('spectrum'),\n",
    "    vanadium=vana_dspacing.bins.concat('spectrum'),\n",
    "    edges=dspacing_edges,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ef778f1-e23a-492f-9a08-9a08e889ad83",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.histogram(all_spectra, bins=dspacing_edges).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df5f7a70-ac69-47dc-b0c1-1139728d2878",
   "metadata": {},
   "source": [
    "### Group into $2\\theta$ bins"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb474efb-692f-4f2d-9243-23ff9319dd76",
   "metadata": {},
   "source": [
    "For a better resolution, we now group the sample and vanadium data into a number of bins in the scattering angle $2\\theta$ (see [here](https://scipp.github.io/scippneutron/user-guide/coordinate-transformations.html))\n",
    "and normalize each individually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a412534e-06e5-4723-9856-37e652dfacf7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "two_theta = sc.linspace(dim='two_theta',\n",
    "                        unit='deg',\n",
    "                        start=25.0,\n",
    "                        stop=90.0,\n",
    "                        num=16)\n",
    "sample_by_two_theta = diffraction.group_by_two_theta(sample_dspacing,\n",
    "                                                     edges=two_theta)\n",
    "vana_by_two_theta = diffraction.group_by_two_theta(vana_dspacing,\n",
    "                                                   edges=two_theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0532df2d-a907-40ed-8271-2dc3a71b04ee",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "normalized = diffraction.normalize_by_vanadium(sample_by_two_theta,\n",
    "                                               vanadium=vana_by_two_theta,\n",
    "                                               edges=dspacing_edges)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84c1b8a6-504f-4627-9463-5ff70a6575fd",
   "metadata": {},
   "source": [
    "Histogram the results in order to get a useful binning in the following plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "014d6ef4-7905-42bf-8336-df7095de47dc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "normalized = sc.histogram(normalized, bins=dspacing_edges)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "460c3a5e-2ab6-4dba-8f78-e4cd126960e8",
   "metadata": {},
   "source": [
    "Now we can inspect the d-spacing distribution as a function of $2\\theta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e66a78f7-e91c-4f81-a3ad-e77f6a4c05ae",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "normalized.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "336a7157-1832-4616-b97f-1372fdad3c50",
   "metadata": {},
   "source": [
    "In order to get 1-dimensional plots, we can select some ranges of scattering angles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "694c8d66-58ea-4645-adca-b33d95f81d23",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "angle = sc.midpoints(normalized.coords['two_theta']).values\n",
    "results = {\n",
    "    f'{round(angle[group], 3)} rad':\n",
    "    normalized['two_theta', group]\n",
    "    for group in range(2, 6)\n",
    "}\n",
    "sc.plot(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02bd3242-7e90-4279-8cbf-341f2bd37d0c",
   "metadata": {},
   "source": [
    "Or interactively by plotting with a 1d projection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c41d0965-dfef-40e3-bb9a-5629f51c8b68",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "normalized.plot(projection='1d')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}