Merge pull request #180 from scipp/wrapped-tof-lookup

Make a wrapped lookup table to time-of-flight

Author: nvaytet

docs/user-guide/tof/dream.ipynb

@@ -239,7 +239,7 @@
     "raw_data = ess_beamline.get_monitor(\"detector\")[0]\n",
     "\n",
     "# Visualize\n",
-    "raw_data.hist(event_time_offset=300).sum(\"pulse\").plot()"
+    "raw_data.hist(event_time_offset=300).squeeze().plot()"
    ]
   },
   {
@@ -290,8 +290,10 @@
     "    time_of_flight.providers(), params=time_of_flight.default_parameters()\n",
     ")\n",
     "workflow[time_of_flight.RawData] = raw_data\n",
-    "workflow[time_of_flight.LtotalRange] = (sc.scalar(75.5, unit='m'),\n",
-    "                                        sc.scalar(78.0, unit='m'))\n",
+    "workflow[time_of_flight.LtotalRange] = (\n",
+    "    sc.scalar(75.5, unit=\"m\"),\n",
+    "    sc.scalar(78.0, unit=\"m\"),\n",
+    ")\n",
     "\n",
     "workflow.visualize(time_of_flight.TofData)"
    ]
@@ -314,8 +316,7 @@
    "outputs": [],
    "source": [
     "workflow[time_of_flight.SimulationResults] = time_of_flight.simulate_beamline(\n",
-    "    choppers=disk_choppers,\n",
-    "    neutrons=2_000_000\n",
+    "    choppers=disk_choppers, neutrons=2_000_000\n",
     ")"
    ]
   },
@@ -343,17 +344,24 @@
    "outputs": [],
    "source": [
     "sim = workflow.compute(time_of_flight.SimulationResults)\n",
-    "# Compute time-of-arrival at the detector\n",
-    "tarrival = sim.time_of_arrival + ((Ltotal - sim.distance) / sim.speed).to(unit=\"us\")\n",
-    "# Compute time-of-flight at the detector\n",
-    "tflight = (Ltotal / sim.speed).to(unit=\"us\")\n",
-    "\n",
-    "events = sc.DataArray(\n",
-    "    data=sim.weight,\n",
-    "    coords={\"wavelength\": sim.wavelength, \"toa\": tarrival, \"tof\": tflight},\n",
-    ")\n",
-    "fig1 = events.hist(wavelength=300, toa=300).plot(norm=\"log\")\n",
-    "fig2 = events.hist(tof=300, toa=300).plot(norm=\"log\")\n",
+    "\n",
+    "\n",
+    "def to_event_time_offset(sim):\n",
+    "    # Compute event_time_offset at the detector\n",
+    "    eto = (\n",
+    "        sim.time_of_arrival + ((Ltotal - sim.distance) / sim.speed).to(unit=\"us\")\n",
+    "    ) % sc.scalar(1e6 / 14.0, unit=\"us\")\n",
+    "    # Compute time-of-flight at the detector\n",
+    "    tof = (Ltotal / sim.speed).to(unit=\"us\")\n",
+    "    return sc.DataArray(\n",
+    "        data=sim.weight,\n",
+    "        coords={\"wavelength\": sim.wavelength, \"event_time_offset\": eto, \"tof\": tof},\n",
+    "    )\n",
+    "\n",
+    "\n",
+    "events = to_event_time_offset(sim)\n",
+    "fig1 = events.hist(wavelength=300, event_time_offset=300).plot(norm=\"log\")\n",
+    "fig2 = events.hist(tof=300, event_time_offset=300).plot(norm=\"log\")\n",
     "fig1 + fig2"
    ]
   },
@@ -374,10 +382,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "table = workflow.compute(time_of_flight.TimeOfFlightLookupTable)\n",
+    "table = workflow.compute(time_of_flight.TimeOfFlightLookupTable).squeeze()\n",
     "\n",
     "# Overlay mean on the figure above\n",
-    "table[\"distance\", 13].plot(ax=fig2.ax, color=\"C1\", ls='-', marker=None)"
+    "table[\"distance\", 13].plot(ax=fig2.ax, color=\"C1\", ls=\"-\", marker=None)"
    ]
   },
   {
@@ -447,7 +455,8 @@
     "# Define wavelength bin edges\n",
     "wavs = sc.linspace(\"wavelength\", 0.8, 4.6, 201, unit=\"angstrom\")\n",
     "\n",
-    "wav_wfm.hist(wavelength=wavs).sum(\"pulse\").plot()"
+    "histogrammed = wav_wfm.hist(wavelength=wavs).squeeze()\n",
+    "histogrammed.plot()"
    ]
   },
   {
@@ -473,7 +482,7 @@
     "\n",
     "pp.plot(\n",
     "    {\n",
-    "        \"wfm\": wav_wfm.hist(wavelength=wavs).sum(\"pulse\"),\n",
+    "        \"wfm\": histogrammed,\n",
     "        \"ground_truth\": ground_truth.hist(wavelength=wavs),\n",
     "    }\n",
     ")"
@@ -530,16 +539,12 @@
    "outputs": [],
    "source": [
     "raw_data = sc.concat(\n",
-    "    [ess_beamline.get_monitor(key)[0] for key in monitors.keys()],\n",
+    "    [ess_beamline.get_monitor(key)[0].squeeze() for key in monitors.keys()],\n",
     "    dim=\"detector_number\",\n",
     ")\n",
     "\n",
     "# Visualize\n",
-    "pp.plot(\n",
-    "    sc.collapse(\n",
-    "        raw_data.hist(event_time_offset=300).sum(\"pulse\"), keep=\"event_time_offset\"\n",
-    "    )\n",
-    ")"
+    "pp.plot(sc.collapse(raw_data.hist(event_time_offset=300), keep=\"event_time_offset\"))"
    ]
   },
   {
@@ -655,22 +660,14 @@
    "source": [
     "# Update workflow\n",
     "workflow[time_of_flight.SimulationResults] = time_of_flight.simulate_beamline(\n",
-    "    choppers=disk_choppers,\n",
-    "    neutrons=2_000_000\n",
+    "    choppers=disk_choppers, neutrons=2_000_000\n",
     ")\n",
     "workflow[time_of_flight.RawData] = ess_beamline.get_monitor(\"detector\")[0]\n",
     "\n",
     "sim = workflow.compute(time_of_flight.SimulationResults)\n",
-    "# Compute time-of-arrival at the detector\n",
-    "tarrival = sim.time_of_arrival + ((Ltotal - sim.distance) / sim.speed).to(unit=\"us\")\n",
-    "# Compute time-of-flight at the detector\n",
-    "tflight = (Ltotal / sim.speed).to(unit=\"us\")\n",
-    "\n",
-    "events = sc.DataArray(\n",
-    "    data=sim.weight,\n",
-    "    coords={\"wavelength\": sim.wavelength, \"toa\": tarrival, \"tof\": tflight},\n",
-    ")\n",
-    "events.hist(wavelength=300, toa=300).plot(norm=\"log\")"
+    "\n",
+    "events = to_event_time_offset(sim)\n",
+    "events.hist(wavelength=300, event_time_offset=300).plot(norm=\"log\")"
    ]
   },
   {
@@ -696,7 +693,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "table = workflow.compute(time_of_flight.TimeOfFlightLookupTable)\n",
+    "table = workflow.compute(time_of_flight.TimeOfFlightLookupTable).squeeze()\n",
     "table.plot() / (sc.stddevs(table) / sc.values(table)).plot(norm=\"log\")"
    ]
   },
@@ -720,9 +717,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "workflow[time_of_flight.LookupTableRelativeErrorThreshold] = 1.0e-2\n",
+    "workflow[time_of_flight.LookupTableRelativeErrorThreshold] = 0.01\n",
     "\n",
-    "workflow.compute(time_of_flight.MaskedTimeOfFlightLookupTable).plot()"
+    "workflow.compute(time_of_flight.TimeOfFlightLookupTable).squeeze().plot()"
    ]
   },
   {
@@ -757,7 +754,7 @@
     "\n",
     "pp.plot(\n",
     "    {\n",
-    "        \"wfm\": wav_wfm.hist(wavelength=wavs).sum(\"pulse\"),\n",
+    "        \"wfm\": wav_wfm.hist(wavelength=wavs).squeeze(),\n",
     "        \"ground_truth\": ground_truth.hist(wavelength=wavs),\n",
     "    }\n",
     ")"