Skip to content

Commit

Permalink
limit time overlap instead of wavelength overlap
Browse files Browse the repository at this point in the history
  • Loading branch information
@nvaytet
nvaytet committed Oct 2, 2024
1 parent 3a414c9 commit 9ac463a
Showing 1 changed file with 36 additions and 33 deletions.
69 changes: 36 additions & 33 deletions docs/wfm-stitching.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -344,7 +344,8 @@
"source": [
"## Reducing overlap between frames\n",
"\n",
"There is a significant amount of overlap between the two frames, as illustrated by this figure"
"There is a small amount of time overlap between some of the frames, as illustrated by this figure\n",
"(for example, some of the neutrons that belong to <span style=\"color: #2ca02c;\">frame 2</span> are bleeding into the range that belongs to <span style=\"color: #ff7f0e;\">frame 1</span>)."
]
},
{
Expand All @@ -355,49 +356,47 @@
"outputs": [],
"source": [
"# Define a common binning that covers the entire range\n",
"full_wavelengths = (\n",
"full_tofs = (\n",
" full_beamline_results.detectors[\"detector\"]\n",
" .wavelengths.visible.data['pulse:0']\n",
" .coords['wavelength']\n",
")\n",
"bins = sc.linspace(\n",
" 'wavelength', full_wavelengths.min(), full_wavelengths.max(), num=301\n",
" .tofs.visible.data['pulse:0']\n",
" .coords['tof']\n",
")\n",
"bins = sc.linspace('tof', full_tofs.min(), full_tofs.max(), num=301)\n",
"\n",
"pp.plot(\n",
" {\n",
" f\"Frame-{i}\": res.detectors[\"detector\"]\n",
" .wavelengths.visible.data['pulse:0']\n",
" .hist(wavelength=bins)\n",
" .tofs.visible.data['pulse:0']\n",
" .hist(tof=bins)\n",
" for i, res in enumerate(results)\n",
" }\n",
")"
]
},
{
"cell_type": "markdown",
"id": "542d589d-0a46-4503-a37b-ce2435403aad",
"id": "8f76007b-da07-4fa7-8ce2-98d57bcbd2c4",
"metadata": {},
"source": [
"Reducing or removing overlap between frames is beneficial as it reduces the appearance of artifacts between frames in the final wavelength spectrum.\n",
"Regions where time of flight overlaps between frames are regions where the wavelength of a neutron cannot accurately be determined,\n",
"as neutrons have the same arrival time at the detector,\n",
"but different start times at the source end of the beamline.\n",
"\n",
"To achieve this as follows:\n",
"They are usually a sign of a faulty design in the chopper cascade,\n",
"and we thus need to discard neutrons from any overlapping regions in our final result.\n",
"\n",
"1. we use percentiles (3% and 97%) to remove wavelength outliers in each frame.\n",
"1. we find the maximum wavelength in frame $i$ and the minimum wavelength in frame $i+1$.\n",
"1. we compute the mean of these two wavelengths.\n",
"1. we remove all neutrons in frame $i$ with wavelengths above the mean.\n",
"1. we remove all neutrons in frame $i+1$ with wavelengths below the mean."
"One way to achieve this is to remove the outliers from the edges of the frames, until the frames no longer overlap.\n",
"This is done below by gradually increasing percentile threshold on the `tof` reading at the detector until overlap is gone."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "63564c46-37ae-4d62-9ee5-01b9c5fa5125",
"id": "8eb2487f-8a15-473b-bd21-ef5142be857c",
"metadata": {},
"outputs": [],
"source": [
"percentile = 3.0\n",
"percentile_step = 0.1 # Step by which to increase percentile thresholds\n",
"\n",
"non_overlapping = []\n",
"for res in results:\n",
Expand All @@ -406,15 +405,20 @@
"\n",
"for i in range(len(results) - 1):\n",
" # Find mean wavelength between frames\n",
" wavs1 = non_overlapping[i].coords['wavelength']\n",
" wavs2 = non_overlapping[i + 1].coords['wavelength']\n",
" # Remove wavelength outliers\n",
" wmax = np.percentile(wavs1.values, 100 - percentile)\n",
" wmin = np.percentile(wavs2.values, percentile)\n",
" wmean = 0.5 * (wmin + wmax) * wavs1.unit\n",
" # Filter on wavelengths\n",
" non_overlapping[i] = non_overlapping[i][wavs1 <= wmean]\n",
" non_overlapping[i + 1] = non_overlapping[i + 1][wavs2 >= wmean]\n",
" tofs1 = non_overlapping[i].coords['tof']\n",
" tofs2 = non_overlapping[i + 1].coords['tof']\n",
" p = 0.0\n",
" overlap = True\n",
" while overlap:\n",
" tofmin = np.percentile(tofs2.values, p)\n",
" tofmax = np.percentile(tofs1.values, 100 - p)\n",
" overlap = tofmin < tofmax\n",
" p += percentile_step\n",
" # Filter on tofs\n",
" non_overlapping[i] = non_overlapping[i][tofs1 <= sc.scalar(tofmax, unit=tofs1.unit)]\n",
" non_overlapping[i + 1] = non_overlapping[i + 1][\n",
" tofs2 >= sc.scalar(tofmin, unit=tofs2.unit)\n",
" ]\n",
"\n",
"for frame in non_overlapping:\n",
" w = frame.coords[\"wavelength\"]\n",
Expand Down Expand Up @@ -544,8 +548,7 @@
"id": "54b85305-c89b-4077-97c2-10c447793725",
"metadata": {},
"source": [
"It is important to note here that, because we removed the wavelength overlap between frames,\n",
"the max wavelength of frame $i$ should be equal (or very close) to the min wavelength of frame $i+1$."
"It is important to note here that the maximum wavelength of frame $i$ should be close to the minimum wavelength of frame $i+1$."
]
},
{
Expand Down Expand Up @@ -714,11 +717,11 @@
"source": [
"pp.plot(\n",
" {\n",
" \"naive\": wav_naive.hist(wavelength=bins),\n",
" \"wfm\": wav_wfm.hist(wavelength=bins),\n",
" \"naive\": wav_naive.hist(wavelength=300),\n",
" \"wfm\": wav_wfm.hist(wavelength=300),\n",
" \"truth\": full_beamline_results[\"detector\"]\n",
" .wavelengths.data[\"visible\"][\"pulse:0\"]\n",
" .hist(wavelength=bins),\n",
" .hist(wavelength=300),\n",
" }\n",
")"
]
Expand Down

0 comments on commit 9ac463a

Please sign in to comment.