lass_j/AMBER
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 4 Wiki Activity

Add initial documentation
Gitea Actions Demo / build_test (3.13.3) (push) Successful in 11s
Details
Gitea Actions Demo / release-test (push) Successful in 1s
Details

Browse Source 
Add synthetic data tutorial
Add read the docs setup
This commit is contained in:
Jakob Lass committed 2025-04-11 14:02:35 +02:00
1 parent 1dee498b38
commit 682c9aeac9
9 files changed
+833 -9
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

.gitignore
+2
View File
@@ -3,3 +3,5 @@ tests/__pycache__/*
build/*
src/AMBER.egg-info/*
source/tutorials/.ipynb_checkpoints/*
make.bat
.readthedocs.yaml
+19
View File
@@ -0,0 +1,19 @@
# Read the Docs configuration file
# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
# Required
version: 2
# Set the OS, Python version, and other tools you might need
build:
os: ubuntu
tools:
python: "3.13"
# Build documentation in the "docs/" directory with Sphinx
sphinx:
configuration: source/conf.py
python:
install:
- requirements: docs/requirements.txt
Makefile
+20
View File
@@ -0,0 +1,20 @@
# Minimal makefile for Sphinx documentation
#
# You can set these variables from the command line, and also
# from the environment for the first two.
SPHINXOPTS ?=
SPHINXBUILD ?= sphinx-build
SOURCEDIR = source
BUILDDIR = build
# Put it first so that "make" without argument is like "make help".
help:
@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
.PHONY: help Makefile
# Catch-all target: route all unknown targets to Sphinx using the new
# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
%: Makefile
@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
README.md
+2 -9
View File
@@ -1,14 +1,7 @@
AMBER
=====
[![Test status](https://gitea.psi.ch/lass_j/AMBER/actions/workflows/test.yaml/badge.svg)]()
General introduction
## Installation
## Quick Installation
```
pip install AMBER
@@ -22,7 +15,7 @@ python3 -m pip install AMBER
## Documentation and Tutorials
The documentation for AMBER can be found at
## Contribute
source/conf.py
+37
View File
@@ -0,0 +1,37 @@
# Configuration file for the Sphinx documentation builder.
#
# For the full list of built-in configuration values, see the documentation:
# https://www.sphinx-doc.org/en/master/usage/configuration.html
# -- Project information -----------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information
project = 'AMBER'
copyright = '2025, J. Lass, V. Cohen, B. B. Haro, & D. G. Mazzone'
author = 'J. Lass, V. Cohen, B. B. Haro, & D. G. Mazzone'
release = '0.0.2'
# -- General configuration ---------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
extensions = [
'nbsphinx', # allows notebooks
'myst_parser' # allows markdown
]
source_suffix = {
'.rst': 'restructuredtext',
'.txt': 'markdown',
'.md': 'markdown',
}
templates_path = ['_templates']
exclude_patterns = []
# -- Options for HTML output -------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
html_theme = 'pyramid'
html_static_path = ['_static']
output_dir = ''
source/index.rst
+94
View File
@@ -0,0 +1,94 @@
.. AMBER documentation master file, created by
sphinx-quickstart on Fri Apr 11 10:48:14 2025.
You can adapt this file completely to your liking, but it should at least
contain the root `toctree` directive.
AMBER documentation
===================
This is the documentation for the **A** lgorithm for **M** ultiplexing spectrometer
**B** ackground **E** stimation with **R** otation-independence **AMBER** which is a method to estimate and background within 3D data sets.
This software is the result of the collaboration between the instrument responsible at the `CAMEA <https://www.psi.ch/en/sinq/camea>`_ beamline at SINQ, at the
`Paul Scherrer Institut <https://psi.ch/>`_ with the `Swiss Data Science Center <https://www.datascience.ch/>`_ through the "Data Science for Multiplexing Spectrometers" (DS4MS) project supported by
the fifth call for Collaborative Research Data Science Projects.
Badges
------
.. image:: https://gitea.psi.ch/lass_j/AMBER/actions/workflows/test.yaml/badge.svg
:width: 200
:alt: Test Status
Installation
------------
Installation can be done through PyPI using
.. code-block::
python -m pip install AMBER
Alternatively, the newest version can be installed directly from the source repository through
.. code-block::
python -m pip install git https://gitea.psi.ch/lass_j/AMBER.git@master
General Introduction
--------------------
**AMBER** has been developed with multiplexing neutron spectrometers in mind, i.e. a novel version of the triple axis instrument which acquires multiple data points simultaneously.
However, the algorithm is more general and can be utilized on data with other origins.
The main requirements of the data are
#. Rotation independence of the background in two out of three directions.
#. Smooth change of background along all directions
#. The signal is sparse but continuous in all directions
**AMBER** minimizes the following cost function for the signal and background given the set of parameters :math:`\lambda`, :math:`\beta`, and :math:`\mu` on a voxelated data set:
.. math::
\min_{X,b} \frac{1}{2}\lVert Y-X-\mathcal{R}b\rVert_{2}^2+\lambda\vert| X |\vert_{1} +\frac{\beta}{2} \mathrm{Tr} \left( b^T L_{b} b \right) +\frac{\mu}{2} \boldsymbol{1}_{n_x}^TX^T L_{\omega} X\boldsymbol{1}_{n_y}
where :math:`Y`, :math:`X`, and :math:`b` are the measured signal, the subtracted signal and the rotation independent background, respectively. :math:`\mathcal{R}` is a rotation operator taking
the radial :math:`b` to the full space. These three parameters constitute the fitting term while the rest are denoted penalty terms. :math:`\lambda` controls the weighing of sparcity of the subtracted signal gauge through the L-1 norm
:math:`\vert| X |\vert_{1}`, :math:`\beta` controls the smoothness of the background while :math:`\mu` the smoothness of the subtracted signal.
For further details, see the the paper.
Documentation & Tutorials
-------------------------
For **AMBER** the following tutorials have been created
.. toctree::
:maxdepth: 2
:caption: Contents:
tutorials/index.rst
Contribute
----------
Contact
-------
source/requirements.txt
+4
View File
@@ -0,0 +1,4 @@
docutils<0.18
ipykernel
nbsphinx
sphinx>5
source/tutorials/index.md
+11
View File
@@ -0,0 +1,11 @@
# Tutorials
Below is a list of tutorials utilizing __AMBER__ to perform a background determination and background subtraction. Some of these require additional packages and data to function. For those notebooks, links to the and installation requirements are highlighted.
```{toctree}
---
maxdepth: 1
caption: Tutorials:
---
simulated_data.ipynb
```
source/tutorials/simulated_data.ipynb
+644
View File
@@ -0,0 +1,644 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "709c49bd",
"metadata": {},
"source": [
"# Showcasing of AMBER utilizing synthetic neutron scattering data"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "d8b17015",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from AMBER.background import background\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"id": "63aaf7ad",
"metadata": {},
"source": [
"## Generate signal data\n",
"The data for this tutorial is generated using the expected neutron scattering signal in an inelastic exmperiment measuring MnF$_2$. The dispersion relation, i.e. the energy at a give (H,K,L) position is given by the analytical formula [Yamany et al. 2010](https://www.doi.org/10.1139/P10-081):\n",
"\n",
"$$\\omega(H,K,L) = \\sqrt{\\left(2S z_2 J_2+D+2 z_1 S J_1 \\sin{L \\pi}^2\\right)^2-\\left(2 S z_2 J_2 \\cos(H \\pi)\\cos(K \\pi) \\cos(L \\pi)\\right)^2}$$\n",
"\n",
"In the above $J_1$, $J_2$, and $D$ are the magnetic coupling strengths and single-ion anisotropy which determine the amplitude of the dispersion. The parameters $S = 5/2$, $z_1 = 2$, and $z_2 = 8$ corresponds the the size of the magnetic spin, and the number of nearest and next-nearest neighbours in the MnF$_2$ lattice. \n",
"\n",
"From the above equation, the imporant fact is that it depends on the 3D vector $Q = (H,K,L)$. For simplicity, we assume that $K$ is zero and thus reduce the data to 3D, i.e. $(H,L,\\omega)$\n",
"\n",
"Due to how neutron scattering experiments function, the above dispersion will not be infinitely thin but rather extended. Here, this is replicated by a simple gaussian smearing along $\\omega$\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "4c67ae0c-684b-4d85-92ba-44907594cce6",
"metadata": {},
"outputs": [],
"source": [
"# Define parameters, dispersion, and Smearing\n",
"S = 5/2\n",
"z1 = 2\n",
"z2 = 8\n",
"\n",
"def SpinWave(Q,J1,J2,D):\n",
" return np.sqrt(np.power(2*S*z2*J2+D+2*z1*S*J1*np.sin(Q[2]*np.pi)**2,2.0)-\n",
" \n",
" np.power(2*S*z2*J2*np.cos(Q[0]*np.pi)*np.cos(Q[1]*np.pi)*np.cos(Q[2]*np.pi),2.0))\n",
"\n",
"def Intensity(H,K,L,E):\n",
" sigmaE = 0.25\n",
" omega = SpinWave(np.array([H,K,L]), J1=0.0354,J2=0.1499,D=0.131 )\n",
" I = np.exp(-np.power(omega-E,2.0)/(2*sigmaE**2))\n",
" \n",
" return I"
]
},
{
"cell_type": "markdown",
"id": "b8cbd985",
"metadata": {},
"source": [
"## Define 3D grid domain\n",
"\n",
"The input to AMBER is a 3D cube of intensities which is define below"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "b1e107a4-306c-49dd-95d9-732ba4284c1e",
"metadata": {},
"outputs": [],
"source": [
"# Data will ge simulated along H and L but with K = 0\n",
"\n",
"h = np.linspace(-0.1,2.1,101)\n",
"k = 0\n",
"l = np.linspace(-0.1,2.1,101)\n",
"\n",
"# Choose a sufficient energy range\n",
"e = np.linspace(0.5,8,61)\n",
"\n",
"# Generate the grid upon which the dispersion is calculated\n",
"\n",
"H,L,E = np.meshgrid(h,l,e)\n",
"K = np.zeros_like(H)\n",
"\n",
"# Intensities are calculated with the smearing and scaled\n",
"I = 30*Intensity(H,K,L,E)\n",
"\n",
"# Poisson noice is added to the intensity\n",
"I = np.random.poisson(I).astype(float)"
]
},
{
"cell_type": "markdown",
"id": "5c087a9d-e273-4be6-8151-87169990074b",
"metadata": {},
"source": [
"### Introduction of NaN-values\n",
"\n",
"Triple axis instruments cannot measure all lengths of $Q$. this we can mimic by exchaning the intensity at these points by NaNs"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "9180b624-f76b-4d62-b87c-bda5388ef3f9",
"metadata": {},
"outputs": [],
"source": [
"QLength = np.linalg.norm([H,L],axis=0)\n",
"\n",
"I[QLength<0.35] = np.nan\n",
"I[QLength>2.5] = np.nan\n"
]
},
{
"cell_type": "markdown",
"id": "e1bfcdd4-2525-44f9-9f75-3b9e81c2491f",
"metadata": {},
"source": [
"### Add background\n",
"The main feature of AMBER is the background segmentation, which requires the data to have a background. As describe in the [article](???), background is assumed to:\n",
"\n",
"- Rotation independence of the background.\n",
"- Smooth change of background along energy and $|\\vec{Q}|$.\n",
"- The signal is sparse but continuous in energy and $|\\vec{Q}|$.\n",
"\n",
"In the following, a background is generated with a higher amplitude for low and high $|\\vec{Q}|$ mimicing instrumental background artefacts from triple axis neutron experiments - these corresponds to the direct beam contribution and to increased background at larger scattering angles. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "99d07fc8-250c-4870-bfa8-a68dd212ad95",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGxCAYAAAB4AFyyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABVQElEQVR4nO3dd3hT9eIG8Ddpm3Qm3Qs6aWmBthSQUfYeIrJUEC9D8SoKIgJ6xfsTQVFc4Bb1qiAigoMhIHuDZXUIpVDaUtrSPWjSPZLz+6NQraympD0Z7+d58jw2SdOXY2hezvkOiSAIAoiIiIiMkFTsAERERETNxSJDRERERotFhoiIiIwWiwwREREZLRYZIiIiMlosMkRERGS0WGSIiIjIaLHIEBERkdGyFDtAS9NqtcjOzoaDgwMkEonYcYiIiKgJBEFAaWkpvL29IZXe/ryLyReZ7Oxs+Pj4iB2DiIiImiEzMxNt27a97eMmX2QcHBwA1B8IhUIhchoiIiJqCrVaDR8fn4bP8dsx+SJz43KSQqFgkSEiIjIydxsWwsG+REREZLRYZIiIiMhoscgQERGR0WKRISIiIqPFIkNERERGi0WGiIiIjBaLDBERERktFhkiIiIyWiwyREREZLRYZIiIiMhoscgQERGR0WKRISIiIqPFItNM2/7MxqzvYxCdWiR2FCIiIrPFItNMx5ILset8Lg5czBM7ChERkdlikWmmvsGuAICjyYUiJyEiIjJfLDLN1CfIFRIJcDG3FPmlVWLHISIiMkssMs3kbCdDmLcSAHA8hWdliIiIxMAicw94eYmIiEhcLDL3oF9QfZE5llwIQRBETkNERGR+WGTuQTd/J1hbSZFfWo1LeWVixyEiIjI7LDL3QG5pgZ4BLgCAo8kFIqchIiIyPywy96gfx8kQERGJhkXmHvULdgMAnEwrQnWdRuQ0RERE5oVF5h6197CHu4McVbVaxKRfEzsOERGRWWGRuUcSiQR9g3h5iYiISAwsMnrQr/1f07CJiIio9bDI6EGf62dkErJVKC6vETkNERGR+WCR0QN3B2uEejpAELhdARERUWtikdGTG9OweXmJiIio9bDI6Enf69OwjyYXcLsCIiKiVsIioyc9/J0hs5AiW1WFy4XlYschIiIyCywyemIjs0D3ACcAvLxERETUWlhk9Khv0F+Xl4iIiKjliVpkVq1ahYiICCgUCigUCkRFRWHnzp0Nj1dVVWH27NlwcXGBvb09Jk6ciLy8PBET39mNAb8nLhejVqMVOQ0REZHpE7XItG3bFm+//TZiYmJw5swZDB48GGPHjsX58+cBAC+88AK2bduGn3/+GYcPH0Z2djYmTJggZuQ76uilgLOdDGXVdYjPLBE7DhERkcmTCAY2xcbZ2RnvvfceHnroIbi5uWH9+vV46KGHAAAXL15Ehw4dEB0djV69ejXp9dRqNZRKJVQqFRQKRUtGBwA892Mctv2ZjbmDgzB/eEiL/zwiIiJT1NTPb4MZI6PRaLBhwwaUl5cjKioKMTExqK2txdChQxueExoaCl9fX0RHR9/2daqrq6FWqxvdWlO/G/sucWE8IiKiFid6kTl37hzs7e0hl8sxa9YsbN68GR07dkRubi5kMhkcHR0bPd/DwwO5ubm3fb3ly5dDqVQ23Hx8fFr4T9BY3+vjZP7MLIGqsrZVfzYREZG5Eb3IhISEID4+HidPnsQzzzyD6dOnIzExsdmvt2jRIqhUqoZbZmamHtPenbejDdq52UErANGpPCtDRETUkizFDiCTyRAUFAQA6NatG06fPo2PPvoIkyZNQk1NDUpKShqdlcnLy4Onp+dtX08ul0Mul7d07DvqF+yG1IJyHE0uxMgwL1GzEBERmTLRz8j8k1arRXV1Nbp16wYrKyvs37+/4bGkpCRkZGQgKipKxIR317DvEsfJEBERtShRz8gsWrQIo0aNgq+vL0pLS7F+/XocOnQIu3fvhlKpxMyZMzF//nw4OztDoVDgueeeQ1RUVJNnLImlZ6ALLKUSpBdVIKOoAr4utmJHIiIiMkmiFpn8/HxMmzYNOTk5UCqViIiIwO7duzFs2DAAwAcffACpVIqJEyeiuroaI0aMwOeffy5m5Caxl1uiq68TTl0pxtGUAjzm4id2JCIiIpNkcOvI6FtrryNzwyf7k7Fi7yWM7OSJL6Z2a7WfS0REZAqMbh0ZU3NjGvYfqYXQaE26KxIREYmGRaaFRLR1hMLaEuqqOpy9WiJ2HCIiIpPEItNCLKQS9Lmxym8yZy8RERG1BBaZFnTj8tIxFhkiIqIWwSLTgvoFuQEAYjOuoay6TuQ0REREpodFpgX5utjCz8UWdVoBJ1KLxI5DRERkclhkWljfIK7yS0RE1FJYZFrYje0KjiYXiJyEiIjI9LDItLCodq6QSoDUgnJkl1SKHYeIiMiksMi0MKWNFTr7OALg7CUiIiJ9Y5FpBf2C62cvHeU4GSIiIr1ikWkFN8bJHE8phJbbFRAREekNi0wriPRxhL3cEsXlNUjMUYsdh4iIyGSwyLQCKwspegW6AOB2BURERPrEItNKOA2biIhI/1hkWsmNfZfOXLmGyhqNyGmIiIhMA4tMKwl0tUMbRxvUaLQ4mcbtCoiIiPSBRaaVSCSSv7Yr4DgZIiIivWCRaUU3Li9x3yUiIiL9YJFpRX2CXCGRABdzS5GvrhI7DhERkdFjkWlFznYyhHkrAfCsDBERkT6wyLSyhstLHCdDRER0z1hkWlnDejIphRAEbldARER0L1hkWlk3PyfYWFmgoLQaSXmlYschIiIyaiwyrUxuaYEeAc4AeHmJiIjoXrHIiODG5aUjLDJERET3hEVGBP2C3QAAp9KKUFXL7QqIiIiai0VGBO097OHuIEdVrRax6dfEjkNERGS0WGREIJFIGqZh8/ISERFR87HIiKRfw3YFBSInISIiMl4sMiLpc30DyYQsNYrKqkVOQ0REZJxYZETi7mCNUE8HAMDx1CKR0xARERknFhkRNVxeSublJSIiouZgkRFR3+vTsI8mc7sCIiKi5mCREVEPf2fILKXIUVUhtaBc7DhERERGh0VGRDYyC3T3dwLAy0tERETNwSIjsr5Bf11eIiIiIt2wyIjsxoDfE5eLUKvRipyGiIjIuLDIiKyjlwIudjKU12gQl1EidhwiIiKjwiIjMqlUgt7XF8c7ynEyREREOmGRMQA3Li9xnAwREZFuWGQMwI0ic/ZqCVQVtSKnISIiMh4sMgbAS2mDdm520ArAH6k8K0NERNRULDIGot+NVX5TWGSIiIiaikXGQPRvX3956XBSAbcrICIiaiIWGQPRK9AFMgspskoqkVpQJnYcIiIio8AiYyBsZZboGegMADiUxGnYRERETcEiY0AGhrgDYJEhIiJqKhYZAzIwpH7A76m0YpRX14mchoiIyPCJWmSWL1+O7t27w8HBAe7u7hg3bhySkpIaPWfgwIGQSCSNbrNmzRIpccsKdLWDj7MNajRa/JFaJHYcIiIigydqkTl8+DBmz56NEydOYO/evaitrcXw4cNRXl7e6Hn//ve/kZOT03B79913RUrcsiQSCQY1XF7KFzkNERGR4bMU84fv2rWr0ddr1qyBu7s7YmJi0L9//4b7bW1t4enp2drxRDEwxA1ro9Nx6Po0bIlEInYkIiIig2VQY2RUKhUAwNnZudH9P/zwA1xdXREWFoZFixahoqLitq9RXV0NtVrd6GZMogJdIbOsn4adks9p2ERERHci6hmZv9NqtZg3bx769OmDsLCwhvunTJkCPz8/eHt74+zZs/jPf/6DpKQkbNq06Zavs3z5cixdurS1YuudjcwCvQJdcORSAQ4lFSDYw0HsSERERAZLIhjIMrLPPPMMdu7ciWPHjqFt27a3fd6BAwcwZMgQpKSkoF27djc9Xl1djerq6oav1Wo1fHx8oFKpoFAoWiS7vn17LA2vb09E73YuWP/vXmLHISIianVqtRpKpfKun98GcWlpzpw52L59Ow4ePHjHEgMAPXv2BACkpKTc8nG5XA6FQtHoZmxuTMM+faUYZZyGTUREdFuiFhlBEDBnzhxs3rwZBw4cQEBAwF2/Jz4+HgDg5eXVwunEE+BqBz8XW9RqBPzBTSSJiIhuS9QiM3v2bKxbtw7r16+Hg4MDcnNzkZubi8rKSgBAamoq3njjDcTExODKlSv47bffMG3aNPTv3x8RERFiRm9REokEA9vXn5U5yFV+iYiIbkvUIrNq1SqoVCoMHDgQXl5eDbeNGzcCAGQyGfbt24fhw4cjNDQUCxYswMSJE7Ft2zYxY7eKgaH168kcTsrnbthERES3Ieqspbt9QPv4+ODw4cOtlMawRAW6QG4pRbaqCpfyyhDiydlLRERE/2QQg33pZtZWFohq5wIAOHCRq/wSERHdCouMARty/fLSgYt5IichIiIyTCwyBmzQ9SITk34N18prRE5DRERkeFhkDFhbJ1uEejpAKwCHL3H2EhER0T+xyBi4wdfPyuznOBkiIqKb6FxkLl++3BI56DaGdPhrGnatRityGiIiIsOic5EJCgrCoEGDsG7dOlRVVbVEJvqbSB8nONvJoK6qQ0z6NbHjEBERGRSdi0xsbCwiIiIwf/58eHp64umnn8apU6daIhsBsJBKGvZe4jRsIiKixnQuMpGRkfjoo4+QnZ2Nb7/9Fjk5Oejbty/CwsKwcuVKFBRwUKq+DQn1AADsv8Bp2ERERH/X7MG+lpaWmDBhAn7++We88847SElJwcKFC+Hj44Np06YhJydHnznNWr/2rrCUSpBaUI4rheVixyEiIjIYzS4yZ86cwbPPPgsvLy+sXLkSCxcuRGpqKvbu3Yvs7GyMHTtWnznNmsLaCj0CnAHw8hIREdHf6VxkVq5cifDwcPTu3RvZ2dlYu3Yt0tPTsWzZMgQEBKBfv35Ys2YNYmNjWyKv2fprGjYvLxEREd2gc5FZtWoVpkyZgvTLine truncated
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Background definition\n",
"def background_simulation(q,amplitude,gamma,mu,amplitude2,gamma2):\n",
" return amplitude*((gamma/(q**2+gamma**2))+amplitude2*np.exp(-np.power(q-mu,2.0)/(2*gamma2**2)))\n",
"\n",
"\n",
"Q = np.linalg.norm([H,K,L],axis=0)\n",
"\n",
"# Choose suitable valies\n",
"gamma = 0.5\n",
"mu = 3.0\n",
"amplitude = 20\n",
"amplitude2 = 1.0\n",
"gamma2 = 0.5\n",
"\n",
"# Generate an example of the background for visual inspection\n",
"q = np.linspace(0.25,Q.max(),201)\n",
"bg_test = background_simulation(q,amplitude,gamma,mu,amplitude2,gamma2)\n",
"\n",
"# display background amplitude\n",
"fig,ax = plt.subplots()\n",
"ax.plot(q,bg_test)\n",
"ax.set_xlabel('|Q|')\n",
"ax.set_ylabel('Background intensity')\n",
"\n",
"\n",
"## Add background to data\n",
"bg_tmp = background_simulation(Q,amplitude,gamma,mu,amplitude2,gamma2)\n",
"\n",
"I+=np.random.poisson(bg_tmp)"
]
},
{
"cell_type": "markdown",
"id": "0d7f0e5e-6f8f-450d-b1ff-afce5699f64a",
"metadata": {},
"source": [
"## Plot a constant energy plot\n",
"The now background affected data is plotted to show a \"before\" picture. This is done by choosing a specific energy ($\\omega \\sim 4.5$ meV) and plotting the intensity as a color map"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "456b2e68-9d32-4e78-9575-e6a2bcf7e9e5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'L [r.l.u.]')"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAHHCAYAAABXx+fLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAACxGklEQVR4nO29eXxW1bX/v54pT+aBhExMCVOYQVAiTqBGgVqvtNcBrwNah16vfK8Wqy391bm3WOvYW65UqyJtEbVOvdViFcERQUFEVMKUACEJISHz8OQZzu8PL8Fz1idhMybwfN6vFy/NYp9z9j5nn5PNXp+1lsuyLEsIIYQQQqIId3d3gBBCCCHkWMMFECGEEEKiDi6ACCGEEBJ1cAFECCGEkKiDCyBCCCGERB1cABFCCCEk6uACiBBCCCFRBxdAhBBCCIk6uAAihBBCSNTBBRAhhBBCog4ugAghioULF4rL5er0zyeffNJtfTvvvPPE5XLJ7NmzjdpPmTIFjmHatGmqbSAQkJ/97GeSm5srcXFxUlhYKG+//faRHsIh88gjj4jL5ZJ33nmn0zZPPfWUuFwu+dvf/nYMe0bI8Ye3uztACOm53HfffZKfn6/sgwcP7obeiLzyyiuycuXKgz6ub9++Mm/ePJstNzdXtbvmmmvkr3/9q9x6660yZMgQWbhwoXzve9+T5cuXyxlnnHHI/T5SzJw5U26//XZZvHixFBUVwTaLFy+W9PR0mT59+jHuHSHHF1wAEUI6Zfr06XLyySd3dzdERKStrU1uu+02+dnPfiZ33XXXQR2bkpIiV155ZZdtVq9eLUuWLJHf/va38tOf/lRERK6++moZNWqU3HHHHfLxxx8fct+PFLm5uXL22WfLK6+8Ik888YT4/X7b3+/atUvef/99ufHGG8Xn83VTLwk5PqALjBByXPDggw9KJBLpWJwcLKFQSJqamjr9+7/+9a/i8Xjkxhtv7LDFxsbKddddJytXrpSdO3d2ef4pU6bIqFGjZP369TJ58mSJj4+XwYMHy1//+lcREXnvvfeksLBQ4uLipKCgALqxdu3aJT/60Y8kKytL/H6/jBw5Up555hlbmyuvvFLq6+vljTfeUMcvWbJEIpGIXHHFFV32lRDCBRAhpAvq6+ulurra9qempuaAxzU1Nanj0J/6+nqjfuzYsUMeeOAB+c1vfiNxcXEHPY5NmzZJQkKCJCUlSXZ2ttx5550SDAZtbT7//HMZOnSoJCcn2+wTJ04UEZF169Yd8Dq1tbXy/e9/XwoLC+XBBx8Uv98vM2fOlBdeeEFmzpwp3/ve9+SBBx6Q5uZmufjii6WxsbHj2N27d8upp54q77zzjsyePVsef/xxGTx4sFx33XXy2GOPdbT74Q9/KLGxsbJ48WJ1/cWLF8uAAQPk9NNPP4i7Q0iUYhFCiINnn33WEhH4x+/3H/D4WbNmdXr8d/9MnjzZqD8XX3yxddppp3X8LCLWzTffbHTsj370I+uee+6xXn75ZWvRokXWv/zLv1giYl166aW2diNHjrTOOeccdfxXX31liYi1YMGCLq8zefJkS0SsxYsXd9g2btxoiYjldrutTz75pMP+1ltvWSJiPfvssx226667zsrJybGqq6tt5505c6aVkpJitbS0dNguueQSKzY21qqvr1fXmjt3btc3hBBiWZZlUQNECOmU+fPny9ChQ202j8dzwOPuuOOOA2puRETS0tIO2Gb58uXy8ssvy6pVqw7YFvH000/bfr7qqqvkxhtvlKeeekp+8pOfyKmnnioiIq2trUpTI/KtG2zf3x+IxMREmTlzZsfPBQUFkpqaKn369JHCwsIO+77/37Ztm4iIWJYlL7/8slx66aViWZZUV1d3tJ06daosWbJE1q5d27Gzc+WVV8pLL70kr7zyilxzzTUiIh07QnR/EWIGF0CEkE6ZOHHiIYmgR4wYISNGjDjs64dCIfnP//xPueqqq+SUU0457PPt47bbbpOnnnpK3nnnnY4FUFxcnAQCAdW2ra2t4+8PRN++fcXlctlsKSkp0q9fP2UT+dZlJiKyZ88eqaurkyeffFKefPJJeO6qqqqO/58+fbr06tVLFi9e3LEAev7552Xs2LEycuTIA/aTEMIFECHkKFBfX2+0YxITEyO9evXq9O8XLVokxcXF8oc//EFKS0ttf9fY2CilpaWSmZkp8fHxB9W/fQuSvXv3dthycnJk165dqm1FRYWI4LB5J53tjnVmtyxLREQikYiIfLuzM2vWLNh2zJgxHf/v8/nk0ksvlaeeekp2794tO3bskM2bN8uDDz54wD4SQr6FCyBCyBHnlltukeeee+6A7SZPniwrVqzo9O937NghwWAQinoXLVokixYtkldffVVmzJhxUP3b53rq3bt3h23cuHGyfPlyaWhosAmh97nexo0bd1DXOBh69+4tSUlJEg6HO83v4+SKK66QBQsWyAsvvCAlJSXicrnk8ssvP2p9JOREgwsgQsgR50hpgGbOnAkXHj/4wQ/ke9/7ntxwww02bc3GjRslPj5e+vfvLyIiDQ0N4vf7bdoey7LkV7/6lYh8q6/Zx8UXXywPPfSQPPnkkx2h9oFAQJ599lkpLCxUbqwjicfjkX/913+VxYsXy4YNG2TUqFG2v9+zZ49tsSYicvrpp0teXp78+c9/lrKyMpk8ebL07dv3qPWRkBMNLoAIIZ3yj3/8QzZu3Kjsp512mgwcOLDT446UBmjYsGEybNgw+Hf5+flq52f48OG2XaW1a9fK5ZdfLpdffrkMHjxYWltb5dVXX5WPPvpIbrzxRhk/fnzHsYWFhXLJJZfI3LlzpaqqSgYPHizPPfeclJaWKiH10eCBBx6Q5cuXS2Fhodxwww0yYsQI2bt3r6xdu1beeecdm7tORMTlcsm//du/ya9//WsR+TZrNyHEHC6ACCGd0lnG5WeffbbLBVBPYcCAAXLmmWfKq6++KpWVleJ2u2X48OGyYMECW8LDfSxatEjuvPNO+dOf/iS1tbUyZswY+fvf/y5nnXXWUe9rVlaWrF69Wu677z555ZVX5H/+538kPT1dRo4cKb/5zW/gMVdccYX8+te/Fr/fLxdffPFR7yMhJxIua58KjxBCCCEkSmAmaEIIIYREHVwAEUIIISTq4AKIEEIIIVEHF0CEEEIIiTq4ACKEEEJI1MEFECGEEEKiDuYBAkQiESkvL5ekpCRV2JAQQgghPRPLsqSxsVFyc3PF7e56j4cLIEB5eflRTXtPCCGEkKPHzp07D1gahgsgQFJSkoiITJwyV7ze2A67ryWk2vq279En8McYXSeYnWI/V2W9btTYrExWr1RlC6XratjofJGEON2P9Fhl89W0KduuohRlu+GKN5Xt9/+Yrmy9vrb/nLZ2r2rjCgaVDfXXvVePy0pMUDZTIon6ebnaI/oaMfpfE+6mdkcb/Uq5m3VVdNPn4Ab9iIB+xFQ0KVvzIP283O0672n9YJ+ynfFva2w/t4T1PVr7wmhly/lnpbI557mISFMfv7Ilb9NzPZCm70koQVdWj9sTUDZ3QL+vEb/3gG08O6qUremUAfr84F7G7m5RtrYss0r16PtSWajnSUuebvdq0Xxlm/Xft9h+Ttqhj0NjQP0Ixut5/Y+Xb1U2QrqbhoYG6devX8fv8a7gAgiwz+3l9caK17f/A+z16g+D160/5OI2WwBZXvvH3evWiw5x60WB5QHX9OpfFOh8EXCssx8iIl6P/jB6/LpdXCL4hR+r23kct8QL+uEK61/sqL9ucM/hPTEk4kV9CetrePUvXrfH7iK1POB+uMEixvA5uCPgWK++T16Pniffnbsd5wOJ3z1+vQDyJ9ptIbAAQvMBvQ9wfvl0O69X3/MwGIP49HPwerWr2h0GCyCv94BtPOD9Nb2XXo8eAzoWgb4v6B6743S7xCQ9J5zHen1gAYTGAPph+fS8Tk5OVjZCegom8hWKoAkhhBASdbAWGKChoUFSUlKkqM+/2/5FG87tpdq6G/XWO3KDhJP1v3h9xWW2nyM5vVWbYIbeAveXaveRFaP/FQ/dMTXafYTGFUrS/woOAtfDzguUSdytel2dsMNuS90C/hX7FXAnBtqV6XCeg7sJuKMS9T12tev+NY7Q101etdN+zSQzV1woQ7dzzgcRESs9TdvAuFB/I0l6zgV66R2F7d/Xz+vHU961/bzkyfNUm6xPGpXNFNQPhH8v2BUFeGq0C9DknUDvVyBVHwc2wCApXzcoW8uARGULJoDdTp/+F2viLj3/d52ln+s3NzyhbCva7Of71Q3X6r5l6rEGUnXfEnfp+YX44LXbjdoRcrTY9/u7vr7+gLuU3AEihBBCSNTBBRAhhBBCog6KoLsgOKC3TcDpK9HRIZF0HeHiKt2lbN7cLGVrGZ9n+zm2QkfBIHdXMDdV2VA0iysARKXAfYTcB+5ifd36i0cqW/YKLdI992cfKdtrz59pPz+IPgmna1cBcm2hcblAtJwAdxRydyG3IHqu/jotNA4UZNt+9lVrFxtyf3qrdX+DBTpkE0V8+Yt1pJWpGxNFUEm8njtPrrU/r17g9qJ5g55hcz/9HBJ26hM6I7QOhoZx+v1C13DOnaY++h4ll4DgAfAckKus/NxUZctdVqdsll8/B+QWrM/TcycCbtNfGvXzf7NmjO3nmpEg8q4UuE6BKw69r5WT9PhPv+RhZUvcot2Cb31+r7IRcqzhDhAhhBBCog4ugAghhBASdXABRAghhJCogxqgLvDWtNgSmwXzM1UbFLoMM0EDjUpshf32o1BmFMoLs9eWa80O6kdEQAZipCkCWo7eH1crW8nMDGXr469VtpBBMlykWYkBGiB3BQiXT9L9RWHg6D4hvQ/SaMSU67Bvp/YIncvToMeA+oZ0Qa52rTtCID1OKEnrQhDpWVqjcfGAz20/v/pOkWoD9VTgeSXsVCao4/IZaorQsUhnUj9Ch8D6a+3H9vq0RrVBaQYiIFwLine truncated
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig,ax = plt.subplots()\n",
"\n",
"# Find the closest energy slice\n",
"energy = 4.5\n",
"EIdx = np.argmin(np.abs(E-energy))\n",
"ax.pcolormesh(H[:, :, EIdx], L[:, :, EIdx], I[:, :, EIdx], vmin=0, vmax = 50)\n",
"\n",
"ax.set_title('E = {:.2f} meV'.format(e[EIdx]))\n",
"ax.set_xlabel('H [r.l.u.]')\n",
"ax.set_ylabel('L [r.l.u.]')"
]
},
{
"cell_type": "markdown",
"id": "19859951-3686-46a3-b11e-0cc03b5e97fc",
"metadata": {},
"source": [
"## Plot a constant H map"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "5484ba62-d027-408b-9e34-e0f0e0ab4528",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'E [meV]')"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGzCAYAAAABsTylAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABtcElEQVR4nO2deZhdVZnuvzPXPKeSqlTmEQIJhAAdBgGZQS60z6VFgw0iTh275XLVJn1tkcch0NdWbC+NaCOhWwniANg0gqIERIEQkkASQmVOqlJVVKqSmqtOnXP2vn/QlO79viE7IcnZSb2/56lHz8fae6+1h5NVp9bvvBHXdV0TQgghhAgh0Xx3QAghhBBif2iiIoQQQojQoomKEEIIIUKLJipCCCGECC2aqAghhBAitGiiIoQQQojQoomKEEIIIUKLJipCCCGECC2aqAghhBAitGiiIoQQQojQEs/nwXO5nH3lK1+xH/3oR9bW1mb19fV244032pe+9CWLRCIH3N5xHGtpabHS0tJA7YUQQgiRf1zXtd7eXquvr7do9N0/M8nrROWuu+6ye++91x588EGbM2eOrVq1yj72sY9ZeXm5/d3f/d0Bt29pabEJEyYchZ4KIYQQ4nDT1NRkDQ0N79omrxOVP/7xj3b11VfblVdeaWZmkydPtuXLl9vKlSsDbV9aWmpmbw+0rKzsiPXzeGXu9757SNsVtOHsd2B8DmqJnhjUMtVZqEWHgv0FMtmJ7XIF2C61F2tRPKyV7vIWB2vwceia40Ctem2wT+/iA5j3mUtiu0Q/HiPVlYHa3hMLoVa6A9v11ye8++rG/Zds68WOEIbGFmHf9g1BLda2D2p7LpwIteo13VBzk3hdI8PY52j/IHawt9/z0hlXg/va1YK16kqoDUythlrRa7vwmKXFB+yHmVlmRj3U4ht24LbjaqHkJr33YnQvnjcbxmvvVlVAzSnBmy62D/vr7G6DWrSsFI/LYOeEkR72vk5h39xEAmqRvV1Qc3rwHmb9Ze0YdKxsXP4xmJlTVe7dV1tHoGPmJuK1j+1qx4ZJPCfs+rP+svNpbeQY5D6M9JH7etIY3DYAic3e5zDrDNuKzv8Y+Xf83cjrROWss86y73//+7Zp0yabOXOmvfbaa/bCCy/Yt771Ldo+nU5bOp0eed3b+/YNWFZWponKIRAtIP/KByCWwn9YooU4UYkO40QlWkgmKgGXSrHjumQIsRTWotgViye8fYkl8XGIFuA/mLFkwIlKBicqEfKeEU/gMeJx7HAsiYONJ1i7hK8N2X8M32wZ8QQ5JnnXiEXxpNP+xnCS45KxRnLkforiOCzqfbN2yMWPRPAfwwjpLx1rlMwsybb+fpiZuXGyP9IXdsO6Md9EhR4T70OX7MuJk2tDZu4O6VuUjZ/B+kfb+fpM9u/GyESFtAvaX9aOdi3wtcbz7r/vgp63CLlHYrQfrEbeh0h/2fm0gPdhhNwn7L4OAn2WzAIt28jrROW2226znp4emz17tsViMcvlcvb1r3/dFi1aRNsvXbrU7rjjjqPcSyGEEELki7xaP4888oj9+Mc/toceeshWr15tDz74oH3zm9+0Bx98kLZfsmSJdXd3j/w0NTUd5R4LIYQQ4miS109UvvCFL9htt91m1113nZmZnXzyybZz505bunSp3XDDDdA+lUpZKhXwI0YhhBBCHPPkdaIyMDAAWlIsFjPHIX+LFocdtwrXKUT24t8RS7d5r1HvVLJuYyDYh3PF2/CWq9iC+xuoJQt2xwU6RGAyxd5jlOzGdQZ9E/HvuznyJ98hXMNpxa1YK9mN5zxdgTvcNxMXzo5dgQv03jofD1y5Ke15nUuQtT0pXBeSrsK/PSf68W/UrF1BuhxqlZtw8WumBseV6CCLZAmZ+gqoRdMlntexFlxJnZuFi3oZg7V4bxaSRbeRTlw47NThAsNYTxpqVj8Wty0la0h843B7+6CNO4T7j5aWQC3WSRZdkgWhNnsq1vrw2rgdeI7dOryukR27sebrX66a9JdcQ0uS9SgNddiOjCs6BhdJU4bJItkSMi4y/kiQxcTs2rCxknbsnLN7KbdxE9Ti06ZAzWH3DlkQzoh3eNtla3Ds7N73PyNOLm22J9Ah8ztRueqqq+zrX/+6TZw40ebMmWNr1qyxb33rW3bTTTfls1tCCCGECAl5nah897vftX/8x3+0v/mbv7H29narr6+3T33qU/blL385n90SQgghREjI60SltLTU7r77brv77rvz2Q0hhBBChBRl/QghhBAitGiiIoQQQojQktc//Yj8UrYK7Q2H3BHO+7u8r9/CVemJftwwV4w2T38dWiRDp2Ct7Hn86vZ0NX5b6aQn8BjZYjRa2FfXZ4q88/T+Opy3V76J+++diO2YucQYqEXDp2wrrrYvbMfzufd0tBfKtqPR4T8G+4r+aC+uymfif880XNGf2ofXYWAS3hNFO9FU2begAmq1jfjV7Zkp5KvFmUXjwyX2BbWUWvGcV/8Ovy6f9SPe8hbUIsPkGzzJNx0zw4ddC/84hmeh8pZah/11yVetM0vJZXYM+zp7ZriQdkYMH2Yl5bq9UQDRPZ24L2LzOM2o0FGbh/WNGFPMIsq2o4ISJ+1cYttAG2aLMauGmEYMvy319rb47DPDh5lQ7rwZ2I4ZSOR89s+o8LwuWY3Xnl4HXz+izoGf55G2gVsKIYQQQhxlNFERQgghRGjRREUIIYQQoUUTFSGEEEKEFk1UhBBCCBFaZP2MYrpnor3BMnvcHWWe15WN2KZ7FpolkeEI1IobyWpww9q+03Gleul6bNe2kGTWEDso0U1MIJ+VVLku2Lx97Cpcrb5vJtociX4XamXbh6DWPwFNlUwxnjuWE9Q3Hs9J1StekyKofeKksF2MSAnpSjyXzATKleExalZ1Ybv6KjwuMXwCmTXEyEmQLBKGU415RbQfzI4h+4t2dmM7YiUFyfpJdOC1YWYJu9Z4J3GLhGUpJZgxQsZAzSJi/cRrvXkvbAxG7Bhm+NB+BMyrYTBjhmbsEIsI7glq6ZCHiVhFlIB2EIPlREV78dyxdk495mQVb+7ybseuAzGS3gv6REUIIYQQoUUTFSGEEEKEFk1UhBBCCBFaNFERQgghRGjRREUIIYQQoUXWzygmNW4AaplmYiUUeu2Yoepg89uSZvQN+hqYH4FU/wFXw/dMxXZFKHlYfACtlOJWPC7L9vHjJMgYiGlTvX4QavtmFkKtZwrmzkQzwc5JEMPHDLM4osO4/1gGLa1kSy/U4sX4FlHchOZSrBONgVw1MUtq8JzEe9FoYIYPw9+u/7xZ0CbRj0aSESMp0dgMNaduDNT2XTgNahVr0Q5hFgUzhqJ42jEr5c1t2GYqGhnMNGIw0yrR0hVoW2ZWMcuFZvH4YJbOe7JIiB3DriGDWlqTx2NDci3guhJbihLUBGKZQ4fZejKSYcSea8iTIvdDZi6+Wcc7vH1zc8GnH/pERQghhBChRRMVIYQQQoQWTVSEEEIIEVo0URFCCCFEaNFERQghhBChRdbPKCb5x1KoZcehIRINMJ8tJ/k/A+OwXVEbWjTZImw3UIc1J0H6lsH9RcldzcyaZI93265Z2Kb+BTRG2s7EA+QSaLMUdKFZkwsY7VH5GhoILRdWQG24Hq9hpth7LYr60ZjIJfB6+W0hM7PUPtyWZROVpIlZQ0hXJALViptw2yBmETN8Eh1oZEV27MYD1KAJE+3DbUua0RhiGTtRYvjQPB1m2/iMIXfeDGgSadyF27GMGYI/S8iMWyRuCd7XgX+7DWC0BDVXnD1ot0UK8DpQ+4hcQwbrS7R1D7Yjx/WbRUGPGbQf1Hoito1LzB1mmhnZn0OuNctwMl8+U4Q8N4nt7bidf98OZkHtD32iIoQQQojQoomKEEIIIUKLJipCCCGECC2aqAghhBAitGiiIoQQQojQIutnFNM7FQ2JhfM3QW3TPu+K9s54JbTJFmK+TpwsfE9e2IHF39ZAaejAMSFmZtbXgLVMDebEDFXjre7PCapoRIPIIXbMuJdx/0U70UjZN68casw+KttK8k5SeD7HvoInNJrGvji+TCCWpZMm5g7DSeL448SsYRk+0eFg1lP5Gz2B9meGRkfEZxvFDcfK9pUcRjvCWE4KMSb8mSX7g2XM0G1JLox/W2Y8MUuJWUWxlRuhljvjBKhR+4jl+pDjuh1oETHwCcP7kOUQZWehQph4YQPuP6D1xHKoosSiYilcuW408uJ93j4zgyZKsp9Yhg8zfNj5pTlEBGakMVjWUfpkzJPyW3RB84X8NpOyfoQQQghxXKCJihBCCCFCiyYqQgghhAgtmqgIIYQQIrRooiKEEEKI0KKJihBCCCFCi/TkUUxq3ADUXtoyBWrxpgLPa4yQMyvAzDDrn4h6asETqCIPk/DCTDkqsJUbUNntmYICYfE2EhBHMr1SXd5tE/1MRkS6puP+0xWoIqe6cAztp+G2FWsxnGvPWeQ84SGs6g0cWNn2Ic/r3RdgcOH4Z3uhFtvRhgcgMGUxtYPok0m8U8qITs0Iuj9/kGJhGyrcqXUkvI9ooUyLZeF9qNhyfTSzYDoetgXPu1NNNPYAoXbsfDD9OTKGuP49JBCOaNJGQu7YU8JU2SBjYGor++05Re5Ndyreh0wxz27dDrV4L6rjbFwscC/OQv580LEzdZopyyTMkZ1f9hUGTLFm15Bp5w45BnsOoc/kWWLBiu8FfaIihBBCiNCiiYoQQgghQosmKkIIIYQILXmdqEyePNkikQj8LF68OJ/dEkIIIURIyOti2ldeecVyuT8tOFy/fr1dfPHFdu211+axV0IIIYQIC3mdqIwZ4115feedd9q0adPsvPPOy1OPRhdOI9ogLgn08xs4qU5cbc4Mn0Q3fmCX6sJ23bOwb4luPEb3LNw2MowOxnAZ7i9ThyvpBzq81kRFY7APGItbsR+pfWj4JPrxXNa8jsdovRgNH2YpjXmNmBqEvoaLine truncated
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig,ax = plt.subplots()\n",
"\n",
"h_value = 0.25\n",
"HIdx = np.argmin(np.abs(h-h_value))\n",
"ax.pcolormesh(L[:, HIdx, :],E[:, HIdx, :],I[:, HIdx, :], vmin=0, vmax=50)\n",
"\n",
"ax.set_xlabel('L [r.l.u.]')\n",
"ax.set_ylabel('E [meV]')\n"
]
},
{
"cell_type": "markdown",
"id": "2f7f37f2",
"metadata": {},
"source": [
"## Run denoising algorithm"
]
},
{
"cell_type": "code",
"execution_count": 42,
"id": "3a58ab49",
"metadata": {},
"outputs": [],
"source": [
"# Initialize the AMBER background object\n",
"\n",
"AMBER = background(dtype=np.float32)\n",
"\n",
"# Set the grid sizes\n",
"AMBER.set_gridcell_size(dqx = 0.022, dqy = 0.022, dE = 0.125)\n",
"\n",
"# Alternatively these can be set from the h, l, and e arrays like\n",
"# AMBER.set_volume_from_limits([h[0],l[0],e[0]],[h[-1],l[-1],e[-1]],)\n",
"\n",
"# Input the data\n",
"AMBER.set_binned_data(h, l, e, I)\n",
"\n",
"bins = int((q.max()-q.min())/0.022)\n",
"\n",
"# define maximum radius and number of bins\n",
"AMBER.set_radial_bins(q.max(),n_bins=bins)"
]
},
{
"cell_type": "markdown",
"id": "cfb11870",
"metadata": {},
"source": [
"### Set algorithm parameters ($\\lambda, \\beta, \\mu$)\n",
"\n",
"lambda and mu will be determined as described in the paper and we select beta using cross validation\n",
"\n",
"beta_ is selected using cross validation, i.e. mask out the top q quantile of intensity. The beta value for the lowest Root-Mean-Square-Error is then chosen"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "60e2fbf2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Test - ( 0.1 )\n",
"RMSE - ( 4.4478 0.1 55.50635141061015 ) : 3.4844415\n",
"Test - ( 1.0 )\n",
"RMSE - ( 4.4478 1.0 55.50635141061015 ) : 3.4835355\n",
"Test - ( 10.0 )\n",
"RMSE - ( 4.4478 10.0 55.50635141061015 ) : 3.4762373\n",
"Test - ( 100.0 )\n",
"RMSE - ( 4.4478 100.0 55.50635141061015 ) : 3.4504814\n",
"Test - ( 200.0 )\n",
"RMSE - ( 4.4478 200.0 55.50635141061015 ) : 3.443771\n",
"Test - ( 300.0 )\n",
"RMSE - ( 4.4478 300.0 55.50635141061015 ) : 3.4432797\n",
"Test - ( 400.0 )\n",
"RMSE - ( 4.4478 400.0 55.50635141061015 ) : 3.4455476\n",
"Test - ( 500.0 )\n",
"RMSE - ( 4.4478 500.0 55.50635141061015 ) : 3.44926\n"
]
}
],
"source": [
"\n",
"lambda_tmp = AMBER.MAD_lambda()\n",
"mu_tmp = AMBER.mu_estimator()\n",
"\n",
"\n",
"beta_range_tmp = np.array([0.1,1.0,10.0,100.0,200.0,300.0,400.0,500.0])\n",
"rmse = AMBER.cross_validation(q=0.3,beta_range= beta_range_tmp, lambda_=lambda_tmp, mu_=mu_tmp,n_epochs=15,verbose=False)\n",
"\n",
"\n",
"beta_tmp = beta_range_tmp[np.argmin(rmse)]"
]
},
{
"cell_type": "markdown",
"id": "e81bce2a",
"metadata": {},
"source": [
"### Run the denoising algorithm using the parameters obtained using the heuristic"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "1f9a8679",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Iteration 1\n",
" Loss function: 24608190.0\n",
" Iteration 2\n",
" Loss function: 18320822.0\n",
" Iteration 3\n",
" Loss function: 15589522.0\n",
" Iteration 4\n",
" Loss function: 14332520.0\n",
" Iteration 5\n",
" Loss function: 13731277.0\n",
" Iteration 6\n",
" Loss function: 13459373.0\n",
" Iteration 7\n",
" Loss function: 13346343.0\n",
" Iteration 8\n",
" Loss function: 13302241.0\n",
" Iteration 9\n",
" Loss function: 13285422.0\n",
" Iteration 10\n",
" Loss function: 13278990.0\n",
" Iteration 11\n",
" Loss function: 13276512.0\n",
" Iteration 12\n",
" Loss function: 13275553.0\n",
" Iteration 13\n",
" Loss function: 13275184.0\n",
" Iteration 14\n",
" Loss function: 13275036.0\n",
" Iteration 15\n",
" Loss function: 13274984.0\n",
" Iteration 16\n",
" Loss function: 13274959.0\n",
" Iteration 17\n",
" Loss function: 13274950.0\n",
" Iteration 18\n",
" Loss function: 13274948.0\n"
]
}
],
"source": [
"# set number of epochs\n",
"n_epochs = 20\n",
"\n",
"AMBER.denoising(AMBER.Ygrid,lambda_tmp,beta_tmp,mu_tmp,n_epochs,verbose=True) "
]
},
{
"cell_type": "markdown",
"id": "e909ad23",
"metadata": {},
"source": [
"## Compute substracted signal"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "15d23f11",
"metadata": {},
"outputs": [],
"source": [
"# The subtracted signal is given by\n",
"\n",
"Y_sub = AMBER.Ygrid - AMBER.b_grid\n",
"Y_back = AMBER.b_grid\n",
"\n",
"# reshape the data to fit \n",
"Y_sub = Y_sub.reshape(AMBER.E_size,AMBER.Qx_size,AMBER.Qy_size).T\n",
"Y_back = Y_back.reshape(AMBER.E_size,AMBER.Qx_size,AMBER.Qy_size).T\n",
"\n",
"# reshape the observation\n",
"Y_obs = AMBER.Ygrid.reshape(AMBER.E_size,AMBER.Qx_size,AMBER.Qy_size).T"
]
},
{
"cell_type": "markdown",
"id": "bf1032b5",
"metadata": {},
"source": [
"## Display substracted signal"
]
},
{
"cell_type": "code",
"execution_count": 46,
"id": "e9e74383",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAB38AAAN6CAYAAACOo4O+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdeXwlVZ3w/29V3dyb5SY3naWT9JZe6Y1NdhQBEaZlcEERhfFxWHXGAR8VxefB8SfgzNCj84y4IbiDzDgCyuCow9KyyCCL7DQNNL2nO52k09lvlntzq+r3B9MZY+ec7+mQ7r6Ez3te/RrJ+dY5p6rOOXVOVarixXEcCwAAAAAAAAAAAADgDc0/2BUAAAAAAAAAAAAAALx+PPwFAAAAAAAAAAAAgGmAh78AAAAAAAAAAAAAMA3w8BcAAAAAAAAAAAAApgEe/gIAAAAAAAAAAADANMDDXwAAAAAAAAAAAACYBnj4CwAAAAAAAAAAAADTAA9/AQAAAAAAAAAAAGAa4OEvAAAAAAAAAAAAAEwDPPwFAKAIXXPNNeJ5nuzevftgV+WA27PvAAAAADDdvJnXen/swgsvlPnz509JXp7nyTXXXDMleb1eDz30kHieJw899NDBrgoA4E2Mh78AABxA69atk//1v/6XzJ49W1KplMyaNUs+8pGPyLp16w521Q6ooaEhueaaa1gQAwAAAChqa9eulQ9+8IPS3NwspaWlMnv2bDnjjDPkW9/61n4v+7rrrpO77rprv5dj8tJLL8k111wjW7duPWh1AAAA+46HvwAAHCB33nmnHHXUUXL//ffLRRddJN/5znfkkksukQcffFCOOuoo+fd///eDXcUDZmhoSK699toJH/5+8YtflOHh4QNfKQAAAAD4I48++qgcc8wx8vzzz8vHPvYx+fa3vy2XXnqp+L4v3/jGN/Z7+cXw8Pfaa68t+oe/w8PD8sUvfvFgVwMAgKKRONgVAADgzWDTpk3y0Y9+VBYuXCgPP/yw1NfXj6V96lOfkre//e3y0Y9+VF544QVZuHDhQazpeFEUST6fl9LS0gNWZiKRkESCKQoAAACAg+sf/uEfJJPJyJNPPinV1dXj0nbt2nVwKmUwODgoFRUVB7saB8WBXK8CAPBGwJu/AAAcAP/0T/8kQ0ND8r3vfW/cg18Rkbq6Ovnud78rg4OD8tWvfnVc2u7du+VDH/qQVFVVSW1trXzqU5+SkZGRcTFr1qyRk046SaqrqyWdTsvSpUvlC1/4wriYXC4nV199tSxevFhSqZTMnTtXPv/5z0sulxsX53meXH755fKv//qvsnLlSkmlUvKrX/1Kampq5KKLLtprv/r7+6W0tFQ+97nPiYhIPp+XL33pS3L00UdLJpORiooKefvb3y4PPvjg2DZbt24dOwbXXnuteJ437m80TfQ3fwuFgvzd3/2dLFq0SFKplMyfP1++8IUv7FX/+fPny7vf/W555JFH5LjjjpPS0lJZuHCh/OQnPxkXNzo6Ktdee60sWbJESktLpba2Vk466SRZs2bNXvsIAAAA4M1p06ZNsnLlyr0e/IqIzJw5c+x/b926VTzPk5tvvnmvONPfo9XWep7nyeDgoNxyyy1ja6YLL7xQRP5nzfTSSy/JX/zFX8iMGTPkpJNOEhGRF154QS688EJZuHChlJaWSmNjo1x88cXS1dW1Vx1aW1vlkksukVmzZkkqlZIFCxbIJz7xCcnn83LzzTfLueeeKyIi73jHO8bq8Mdfb7r77rvl7W9/u1RUVEhlZaWcddZZE/5Jo7vuuksOPfRQKS0tlUMPPXSfvnr11FNPyapVq6Surk7KyspkwYIFcvHFF6vH+KGHHpJjjjlGSktLZdGiRfLd7353wrXmnjXwnjqmUilZuXKl3HPPPePitm3bJn/zN38jS5culbKyMqmtrZVzzz236N+KBgC8OfFaDQAAB8CvfvUrmT9/vrz97W+fMP3kk0+W+fPny29+85txP//Qhz4k8+fPl9WrV8vjjz8u3/zmN6Wnp2fsYea6devk3e9+txx++OHy5S9/WVKplGzcuFF+//vfj+URRZG8973vlUceeUQ+/vGPy/Lly2Xt2rVy/fXXy6uvvrrXZ8QeeOABuf322+Xyyy+Xuro6WbJkibz//e+XO++8U7773e9KMpkci73rrrskl8vJeeedJyKvPQz+wQ9+IOeff7587GMfk4GBAfnhD38oq1atkj/84Q9y5JFHSn19vdx4443yiU98Qt7//vfLBz7wAREROfzww43H79JLL5VbbrlFPvjBD8pnP/tZeeKJJ2T16tXy8ssv73XjYOPGjfLBD35QLrnkErngggvkRz/6kVx44YVy9NFHy8qVK0XktZslq1evlksvvVSOO+446e/vl6eeekqeeeYZOeOMM2ynEgAAAMCbRHNzszz22GPy4osvyqGHHjqleWtrvVtvvXVsvfLxj39cREQWLVo0Lo9zzz1XlixZItddd53EcSwir/1y8ObNm+Wiiy6SxsZGWbdunXzve9+TdevWyeOPPz728HPnzp1y3HHHSW9vr3z84x+XZcuWSWtrq/z85z+XoaEhOfnkk+V//+//Ld/85jflC1/4gixfvlxEZOz/33rrrXLBBRfIqlWr5Ctf+YoMDQ3JjTfeKCeddJI8++yzMn/+fBERue++++Scc86RFStWyOrVq6Wrq0suuugimTNnjnqMdu3aJX/2Z38m9fX18n//7/+V6upq2bp1q9x5553W7Z599ll517veJU1NTXLttddKGIby5S9/ea9fxN7jkUcekTvvvFP+5m/+RiorK+Wb3/ymnHPOOdLS0iK1tbUiIvLkk0/Ko48+Kuedd57MmTNHtm7dKjfeeKOceuqp8tJLL0l5ebm6PwAAHDAxAADYr3p7e2MRid/3vvdZ49773vfGIhL39/fHV199dSwi8Xvf+95xMX/zN38Ti0j8/PPPx3Ecx9dff30sInFnZ6cx31tvvTX2fT/+r//6r3E/v+mmm2IRiX//+9+P/UxEYt/343Xr1o2Lvffee2MRiX/1q1+N+/mf//mfxwsXLhz770KhEOdyuXExPT09cUNDQ3zxxReP/ayzszMWkfjqq6/eq7579n2P5557LhaR+NJLLx0X97nPfS4WkfiBBx4Y+1lzc3MsIvHDDz889rNdu3bFqVQq/uxnPzv2syOOOCI+66yz9iobAAAAAPa477774iAI4iAI4hNPPDH+/Oc/H997771xPp8fF7dly5ZYROIf//jHe+Xxp+se17VeHMdxRUVFfMEFF+yV5548zj///L3ShoaG9vrZv/3bv+21TvrLv/zL2Pf9+Mknn9wrPoqiOI7j+I477ohFJH7wwQfHpQ8MDMTV1dXxxz72sXE/b29vjzOZzLifH3nkkXFTU1Pc29s79rP77rsvFpG4ubl5r7L/2L//+7/HIjJhHf/Ynx7j97znPXF5eXnc2to69rMNGzbEiUQi/tPb4SISJ5PJeOPGjWM/e/7552MRib/1rW+N/Wyi4/rYY4/FIhL/5Cc/GfvZgw8+OOExAwDgQOKzzwAA7GcDAwMiIlJZWWmN25Pe398/9rPLLrtsXMwnP/lJERH5z//8TxGRsc+P/fKXv5QoiibM94477pDly5fLsmXLZPfu3WP/TjvtNBGRcZ9kFhE55ZRTZMWKFeN+dtppp0ldXZ3cdtttYz/r6emRNWvWyIc//OGxnwVBMPZmcBRF0t3dLYVCQY455hh55plnrPtvsmdfr7jiinE//+xnPysistfb0itWrBj3hnV9fb0sXbpUNm/ePPaz6upqWbdunWzYsGFSdQIAAAAw/Z1xxhny2GOPyXvf+155/vnn5atf/aqsWrVKZs+eLf/xH//xuvLW1nou/vqv/3qvn5WVlY3975GREdm9e7eccMIJIiJja7IoiuSuu+6S97znPXLMMcfslceffhr5T61Zs0Z6e3vl/PPPH7fGDIJAjj/++LE1Zltbmzz33HNywQUXSCaTGdv+jDPO2GvNOZE9691f//rXMjo6qsaLiIRhKL/97W/l7LPPllmzZo39fPHixXLmmWdOuM3pp58+7q3qww8/XKqqqsatIf/4uI6OjkpXV5csXrxYqqurJ73WBQBgf+HhLwAA+9meh7p7HgKbTPSQeMmSJeNiFi1aJL7vj/1doQ9/+MPytre9TS699FJpaGiQ8847T26//fZxD4I3bNgg69atk/r6+nH/DjnkEBF57VNaf2zBggV71S2RSMg555wjv/zlL8f+zu6dd94po6Oj4x7+iojccsstcvjhh4/9Ld36+nr5zW9+I319fdb9N9m2bZv4vi+LFy8e9/PGxkaprq6Wbdu2jfv5vHnz9spjxowZ0tPTM/bfX/7yl6W3t1cOOeQQOeyww+TKK6+UF154YVL1AwAAADB9HXvssXLnnXdKT0+P/OEPf5CrrrpKBgYG5IMf/KC89NJLk85XW+u5mGjt1t3dLZ/61KekoaFBysrKpL6+fixuz5qss7NT+vv7J/0p6z2/RHvaaafttc687777xtaYe9Zqf7qvIiJLly5VyznllFPknHPOkWuvvVbq6urkfe97n/z4xz8eW5NOZNeuXTI8PLzX+lFEJvyZiNsacnh4WL70pS/J3LlzJZVKSV1dndTX10tvb++k17oAAOwv/M1fAAD2s0wmI01NTerDxRdeeEFmz54tVVVVxpg//Q3ssrIyefjhh+XBBx+U3/zmN3LPPffIbbfdJqeddprcd999EgSBRFEkhx12mHzta1+bMM+5c+fuledEzjvvPPnud78rd999t5x99tly++23y7Jly+SII44Yi/mXf/kXufDCC+Xss8+WK6+8UmbOnClBEMjq1atl06ZN1v3XaL99vkcQBBP+PP7vv4El8trfWN60aZP88pe/lPvuu09+8IMfyPXXXy833XSTXHrppa+rngAAAACmn2QyKccee6wce+yxcsghh8hFF10Line truncated
"text/plain": [
"<Figure size 2000x900 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# reshape the observation\n",
"Y_obs = AMBER.Ygrid.reshape(AMBER.E_size,AMBER.Qx_size,AMBER.Qy_size).T\n",
"\n",
"\n",
"fig0 = plt.figure(figsize=(20, 9))\n",
"\n",
"# Plot 1: Observations Y\n",
"ax0 = fig0.add_subplot(1, 2, 1)\n",
"\n",
"energy = 4.5\n",
"EIdx = np.argmin(np.abs(E-energy))\n",
"ax0.pcolormesh(H[:,:,EIdx],L[:,:,EIdx],Y_obs[:,:,EIdx],vmin=-2,vmax=50)\n",
"\n",
"ax0.set_xlabel('H [r.l.u.]')\n",
"ax0.set_ylabel('L [r.l.u.]')\n",
"ax0.set_title('Observations')\n",
"\n",
"# Plot 2: Subtracted Y\n",
"ax1 = fig0.add_subplot(1, 2, 2)\n",
"\n",
"energy = 4.5\n",
"EIdx = np.argmin(np.abs(E-energy))\n",
"p = ax1.pcolormesh(H[:,:,EIdx],L[:,:,EIdx],Y_sub[:,:,EIdx],vmin=-2,vmax=50)\n",
"\n",
"ax1.set_xlabel('H [r.l.u.]')\n",
"ax1.set_ylabel('L [r.l.u.]')\n",
"ax1.set_title('Subtracted signal')\n",
"\n",
"fig0.colorbar(p)\n",
"plt.tight_layout()\n",
"plt.show()\n"
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "da0f7346",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAB3wAAAN6CAYAAABllDi9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdeZwdVZ3w/2/dtW/37SXd6U53Z+msZCUsQUhkUwTRAUWMKAyjrK6Iu/OIM4+AMxL1eY04g4iKPMjgoIIP6DCKbLKIEJZAWEII2TqdrbuT3te7VNXvD6R/NOlzviehTd8bP+/XKy+lz/cst+rUqXOqbtX1wjAMBQAAAAAAAAAAAABQdCIT3QAAAAAAAAAAAAAAwIHhhi8AAAAAAAAAAAAAFClu+AIAAAAAAAAAAABAkeKGLwAAAAAAAAAAAAAUKW74AgAAAAAAAAAAAECR4oYvAAAAAAAAAAAAABQpbvgCAAAAAAAAAAAAQJHihi8AAAAAAAAAAAAAFClu+AIAAAAAAAAAAABAkeKGLwAAB8lVV10lnufJ3r17J7opB93rnx0AAAAADhV/y2u8N7rwwgtl5syZ41KW53ly1VVXjUtZb9XDDz8snufJww8/PNFNAQBAxQ1fAADeonXr1sk//MM/yNSpUyWZTEpjY6Ocf/75sm7duolu2kE1ODgoV111FYthAAAAAAXpxRdflA996EPS1NQkJSUlMnXqVDnttNPkuuuu+6vXfc0118hvfvObv3o9Ji+//LJcddVV0tzcPGFtAAAAfz3c8AUA4C2488475eijj5YHH3xQLrroIvnhD38ol1xyiTz00ENy9NFHy1133TXRTTxoBgcH5eqrrx7zhu8///M/y9DQ0MFvFAAAAACIyOOPPy7HHHOMPP/88/Lxj39cfvCDH8ill14qkUhE/v3f//2vXn8h3PC9+uqrC/6G79DQkPzzP//zRDcDAICiE5voBgAAUKw2b94sH/3oR2X27Nny6KOPSm1t7Uja5z//eTnxxBPlox/9qLzwwgsye/bsCWzpaEEQSDablZKSkoNWZywWk1iMaQcAAACAifGtb31LKisr5emnn5aqqqpRae3t7RPTKIOBgQEpKyub6GZMiIO5TgUA4FDCE74AAByg//N//o8MDg7KT37yk1E3e0VEJk+eLD/+8Y9lYGBAvvvd745K27t3r3z4wx+WiooKqampkc9//vMyPDw8Kub++++XE044QaqqqiSdTsv8+fPl61//+qiYTCYjV155pcydO1eSyaRMnz5d/vEf/1EymcyoOM/z5LOf/az813/9lyxevFiSyaTcfffdUl1dLRdddNE+n6u3t1dKSkrkK1/5ioiIZLNZ+cY3viHLli2TyspKKSsrkxNPPFEeeuihkTzNzc0j2+Dqq68Wz/NG/fbSWL/hm8/n5V/+5V9kzpw5kkwmZebMmfL1r399n/bPnDlTzjzzTHnsscfk2GOPlZKSEpk9e7b853/+56i4XC4nV199tcybN09KSkqkpqZGTjjhBLn//vv3+YwAAAAA/rZs3rxZFi9evM/NXhGRurq6kf/f3NwsnufJz372s33iTL8vq63xPM+TgYEBueWWW0bWShdeeKGI/P9rpZdffln+/u//XiZNmiQnnHCCiIi88MILcuGFF8rs2bOlpKRE6uvr5eKLL5aOjo592rBz50655JJLpLGxUZLJpMyaNUs+/elPSzablZ/97GdyzjnniIjIO9/5zpE2vPHtTPfcc4+ceOKJUlZWJuXl5XLGGWeM+TNFv/nNb2TJkiVSUlIiS5Ys2a+3Wj3zzDNy+umny+TJkyWVSsmsWbPk4osvVrfxww8/LMccc4yUlJTInDlz5Mc//vGYa8zX176vtzGZTMrixYvlD3/4w6i4bdu2yWc+8xmZP3++pFIpqampkXPOOafgn34GAMCGR20AADhAd999t8ycOVNOPPHEMdNPOukkmTlzpvzud78b9fcPf/jDMnPmTFm1apWsXr1a/uM//kO6urpGbmCuW7dOzjzzTFm6dKl885vflGQyKZs2bZI///nPI2UEQSDvf//75bHHHpNPfOITsnDhQnnxxRfl2muvlVdffXWfV4X98Y9/lNtvv10++9nPyuTJk2XevHly9tlny5133ik//vGPJZFIjMT+5je/kUwmI+eee66IvHYD+Kc//amcd9558vGPf1z6+vrkpptuktNPP12eeuopOfLII6W2tlZuuOEG+fSnPy1nn322fPCDHxQRkaVLlxq336WXXiq33HKLfOhDH5Ivf/nL8uSTT8qqVatk/fr1+1w02LRpk3zoQx+SSy65RC644AL5v//3/8qFF14oy5Ytk8WLF4vIaxdKVq1aJZdeeqkce+yx0tvbK88884w8++yzctppp9l2JQAAAIBDXFNTkzzxxBPy0ksvyZIlS8a1bG2Nd+utt46sUz7xiU+IiMicOXNGlXHOOefIvHnz5JprrpEwDEXktS8Cb9myRS666CKpr6+XdevWyU9+8hNZt26drF69euSG565du+TYY4+V7u5u+cQnPiELFiyQnTt3yq9//WsZHByUk046ST73uc/Jf/zHf8jXv/51WbhwoYjIyP/eeuutcsEFF8jpp58u3/nOd2RwcFBuuOEGOeGEE+S5556TmTNniojIfffdJytXrpRFixbJqlWrpKOjQy666CKZNm2auo3a29vl3e9+t9TW1srXvvY1qaqqkubmZrnzzjut+Z577jl5z3veIw0NDXL11VeL7/vyzW9+c58vXb/usccekzvvvFM+85nPSHl5ufzHf/yHrFy5UlpaWqSmpkZERJ5++ml5/PHH5dxzz5Vp06ZJc3Oz3HDDDfKOd7xDXn75ZSktLVU/DwAABScEAAD7rbu7OxSR8KyzzrLGvf/97w9FJOzt7Q2vvPLKUETC97///aNiPvOZz4QiEj7//PNhGIbhtddeG4pIuGfPHmO5t956axiJRMI//elPo/7+ox/9KBSR8M9//vPI30QkjEQi4bp160bF3nvvvaGIhHffffeov//d3/1dOHv27JH/zufzYSaTGRXT1dUVTpkyJbz44otH/rZnz55QRMIrr7xyn/a+/tlft3bt2lBEwksvvXRU3Fe+8pVQRMI//vGPI39ramoKRSR89NFHR/7W3t4eJpPJ8Mtf/vLI34444ojwjDPO2KduAAAAALjvvvvCaDQaRqPRcMWKFeE//uM/hvfee2+YzWZHxW3dujUUkfDmm2/ep4w3r3dc13hhGIZlZWXhBRdcsE+Zr5dx3nnn7ZM2ODi4z99+8Ytf7LM++tjHPhZGIpHw6aef3ic+CIIwDMPwjjvuCEUkfOihh0al9/X1hVVVVeHHP/7xUX9vbW0NKysrR/39yCOPDBsaGsLu7u6Rv913332hiIRNTU371P1Gd911VygiY7bxjd68jd/3vveFpaWl4c6dO0f+tnHjxjAWi4VvvrQtImEikQg3bdo08rfnn38+FJHwuuuuG/nbWNv1iSeeCEUk/M///M+Rvz300ENjbjMAAAoRr3QGAOAA9PX1iYhIeXm5Ne719N7e3pG/XXbZZaNiLr/8chER+f3vfy8iMvKKsd/+9rcSBMGY5d5xxx2ycOFCWbBggezdu3fk3ymnnCIiMup1yyIiJ598sixatGjU30455RSZPHmy/OpXvxr5W1dXl9x///3ykY98ZORv0Wh05AngIAiks7NT8vm8HHPMMfLss89aP7/J65/1S1/60qi/f/nLXxYR2eep6EWLFo16krq2tlbmz58vW7ZsGflbVVWVrFu3TjZu3HhAbQIAAABw6DrttNPkiSeekPe///3y/PPPy3e/+105/fTTZerUqfLf//3fb6lsbY3n4lOf+tQ+f0ulUiP/f3h4WPbu3SvLly8XERlZiwVBIL/5zW/kfe97nxxzzDH7lPHm1x6/2f333y/d3d1y3nnnjVpbRqNROe6440bWlrt375a1a9fKBRdcIJWVlSP5TzvttH3WmmN5fZ37P//zP5LL5dR4ERHf9+WBBx6QD3zgA9LY2Djy97lz58p73/veMfOceuqpo56eXrp0qVRUVIxaO75xu+ZyOeno6JC5c+dKVVXVAa9xAQCYaNzwBQDgALx+I/f1G78mY90Ynjdv3qiYOXPmSCQSGfm9oI985CNy/PHHy6WXXipTpkyRc889V26//fZRN383btwo69atk9ra2lH/DjvsMBF57XVZbzRr1qx92haLxWTlypXy29/+duR3c++8807J5XKjbviKiNxyyy2ydOnSkd/Gra2tld/97nfS09Nj/fwm27Ztk0gkInPnzh319/r6eqmqqpJt27aN+vuMGTP2KWPSpEnS1dU18t/f/OY3pbu7Ww477DA5/PDD5atf/aq88MILB9Q+AAAAAIeet73tbXLnnXdKV1eXPPXUU3LFFVdIX1+ffOhDH5KXX375gMvV1nguxlqzdXZ2yuc//3mZMmWKpFIpqa2tHYl7fS22Z88e6e3tPeDXVL/+hdlTTjlln/XlfffdN7K2fH2N9ubPKiIyf/58tZ6TTz5ZVq5cKVdffbVMnjxZzjrrLLn55ptH1qJjaW9vl6GhoX3WjSIy5t9E3NaOQ0ND8o1vfEOmT58uyWRSJk+eLLW1tdLd3X3Aa1wAACYav+ELAMABqKyslIaGBvWG4gsvvCBTp06ViooKY8ybv3GdSqXk0UcflYceekh+97vfyR/+8Af51a9+Jaeccorcd999Eo1GJQgCOfzww+V73/vemGVOnz59nzLHcu6558qPf/xjueeee+QDH/iA3H777bJgwQI54ogjRmJ+/vOfy4UXXigf+MAH5Ktf/arU1dVJNBqVVatWyebNm62fX6N92/x10Wh0zL+Hf/ltK5HXfjN58+bN8tvf/lbuu+8++elPfyrXXnut/OhHP5JLL730LbUTAAAAwKEjkUjI2972Nnnb294mhx12mFx00UVyxx13yJVXXmlLine truncated
"text/plain": [
"<Figure size 2000x900 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig0 = plt.figure(figsize=(20, 9))\n",
"\n",
"# Plot 1: Observations Y\n",
"ax0 = fig0.add_subplot(1, 2, 1)\n",
"\n",
"h_value = 0.25\n",
"HIdx = np.argmin(np.abs(h-h_value))\n",
"ax0.pcolormesh(L[:,HIdx,:],E[:,HIdx,:],Y_obs[:,HIdx,:],vmin=-2,vmax=50)\n",
"\n",
"ax0.set_xlabel('H [r.l.u.]')\n",
"ax0.set_ylabel('L [r.l.u.]')\n",
"ax0.set_title('Observations')\n",
"\n",
"# Plot 2: Subtracted Y\n",
"ax1 = fig0.add_subplot(1, 2, 2)\n",
"\n",
"\n",
"HIdx = np.argmin(np.abs(h-h_value))\n",
"p = ax1.pcolormesh(L[:,HIdx,:],E[:,HIdx,:],Y_sub[:,HIdx,:],vmin=-2,vmax=50)\n",
"\n",
"ax1.set_xlabel('L [r.l.u.]')\n",
"ax1.set_ylabel('E [meV]')\n",
"ax1.set_title('Subtracted signal')\n",
"\n",
"fig0.colorbar(p)\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "71d4ed7e-dec7-4c84-842a-d03030e42db3",
"metadata": {},
"source": [
"## Compare the signal with the subtracted signal in a 1D cut"
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "12829c7d-fbfd-4e4e-b3a1-c759d83f834e",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAACFPElEQVR4nO3de3wTZboH8F/SK4U25dI2LUIpyK1yE+RSLgpatgVEUFyFRQFldY+CCyJaWEREdpeCCrrHy+56ATxeEFQUhC0CAioXywJVsIAIFRB6wSJtobTQZM4fIaFpMslMMklmkt/3fDhup5PJ22k6efLO+zyPThAEAURERESkefpAD4CIiIiIlMHAjoiIiChIMLAjIiIiChIM7IiIiIiCBAM7IiIioiDBwI6IiIgoSDCwIyIiIgoSDOyIiIiIgkR4oAegNmazGWfOnEFsbCx0Ol2gh0NEREQhThAEVFVVISUlBXq96zk5BnYNnDlzBq1atQr0MIiIiIjsnDp1Ctddd53LfRjYNRAbGwvAcvLi4uICPBoiIiIKdZWVlWjVqpUtRnGFgV0D1tuvcXFxDOyIiIhINaQsEWPyBBEREVGQYGBHREREFCQY2BEREREFCa6xIyKioGUymXDlypVAD4PIpYiICISFhSlyLAZ2REQUdARBQElJCc6fPx/ooRBJEh8fD6PR6HUNXQZ2REQUdKxBXWJiImJiYlhwnlRLEARUV1ejrKwMAJCcnOzV8RjYERFRUDGZTLagrnnz5oEeDpFbjRo1AgCUlZUhMTHRq9uyTJ4gIqKgYl1TFxMTE+CREElnfb16uyaUgR0REQUl3n4lLVHq9crAjoiIiChIMLBTCZNZwK5j5fis4DR2HSuHySwEekhERKRSbdq0wUsvvRToYShm27Zt0Ol0zGJWAAM7Fcg7WIyBi77EuDd2Y9rKAox7YzcGLvoSeQeLAz00IiLys1OnTuHBBx9ESkoKIiMjkZqaimnTpqG8vDzQQ1PE4MGDMX36dLtt/fv3R3FxMQwGQ2AGFUQY2AVY3sFiPPLuPhRX1NhtL6mowSPv7mNwR0QUQP6+m3L8+HHcdNNNOHr0KD744AP89NNP+Oc//4ktW7YgIyMD586d8+nzizGZTDCbzT47fmRkpCI13IiBXUCZzALmryuEs8uEddv8dYW8LUtEFACBuJsyZcoUREZG4osvvsAtt9yC1q1bY9iwYdi8eTNOnz6NOXPm2PatqqrCuHHj0LhxY7Rs2RKvvvqq7XuCIODZZ59F69atERUVhZSUFPz5z3+2fb+2thYzZ85Ey5Yt0bhxY/Tt2xfbtm2zfX/58uWIj4/H2rVrkZ6ejqioKLz55puIjo52uF06bdo03HrrrQCA8vJyjBs3Di1btkRMTAy6du2KDz74wLbvpEmTsH37drz88svQ6XTQ6XT4+eefnd6K/fjjj3HDDTcgKioKbdq0wYsvvmj3vG3atMHf//53PPjgg4iNjUXr1q3x73//2/b9y5cvY+rUqUhOTkZ0dDRSU1OxcOFCj34vWsLALoDyi845zNTVJwAorqhBflFgPqEREYWqQNxNOXfuHDZu3IhHH33UVtfMymg0Yvz48fjwww8hCJYP+88//zy6d++O/fv3Y9asWZg2bRo2bdoEwBIULV26FP/6179w9OhRfPrpp+jatavteFOnTsWuXbuwcuVKfP/99/j973+P7OxsHD161LZPdXU1Fi1ahDfffBM//PADxo8fj/j4eHz88ce2fUwmEz788EOMHz8eAFBTU4NevXph/fr1OHjwIB5++GHcf//9yM/PBwC8/PLLyMjIwEMPPYTi4mIUFxejVatWDudi7969uOeeezB27FgcOHAAzz77LObOnYvly5fb7ffiiy/ipptuwv79+/Hoo4/ikUcewZEjRwAA//jHP7B27VqsWrUKR44cwXvvvYc2bdp4+NvRDhYoDqCyKvGgzpP9iIjIe+7upuhguZsyNN2IML1ytw6PHj0KQRDQuXNnp9/v3LkzfvvtN5w9exYAMGDAAMyaNQsA0KFDB+zYsQNLly7F0KFDcfLkSRiNRmRmZiIiIgKtW7dGnz59AAAnT57EsmXLcPLkSaSkpAAAZs6ciby8PCxbtgx///vfAVjqqb322mvo3r27bQxjx47F+++/j8mTJwMAtmzZgvPnz2PMmDEAgJYtW2LmzJm2/R977DFs3LgRq1atQp8+fWAwGBAZGYmYmBgYjUbRc7FkyRLcdtttmDt3ru3nKywsxPPPP49JkybZ9hs+fDgeffRRAEBOTg6WLl2KrVu3omPHjjh58iTat2+PgQMHQqfTITU1VeJvQts4YxdAibHRiu5HRETeC/TdFOuMnDsZGRkOXx86dAgA8Pvf/x6XLl1C27Zt8dBDD2HNmjWoq6sDABw4cAAmkwkdOnRAkyZNbP+2b9+OY8eO2Y4XGRmJbt262T3H+PHjsW3bNpw5cwYA8N5772HEiBGIj48HYJnBW7BgAbp27YpmzZqhSZMm2LhxI06ePCnrHBw6dAgDBgyw2zZgwAAcPXoUJpPJtq3++HQ6HYxGo60116RJk1BQUICOHTviz3/+M7744gtZY9AqBnYB1CetGZIN0RD7vKcDkGyIRp+0Zv4cFhFRSAvU3ZTrr78eOp3OFpw1dOjQITRt2hQJCQluj9WqVSscOXIEr732Gho1aoRHH30UN998M65cuYILFy4gLCwMe/fuRUFBge3foUOH8PLLL9uO0ahRI4dkht69e6Ndu3ZYuXIlLl26hDVr1thuwwKW28Mvv/wycnJysHXrVhQUFCArKwuXL1/28Ky4FhERYfe1TqezJXn07NkTRUVFWLBgAS5duoR77rkHd999t0/GoSa8FRtAYXod5o1MxyPv7oMOsJv2t/4pzRuZruhUPxERuRaouynNmzfH0KFD8dprr+Hxxx+3W2dXUlKC9957DxMmTLAFW7t377Z7/O7du+1u4zZq1AgjR47EyJEjMWXKFHTq1AkHDhzAjTfeCJPJhLKyMgwaNEj2OMePH4/33nsP1113HfR6PUaMGGH73o4dOzBq1Cjcd999AACz2Ywff/wR6enptn0iIyPtZt2c6dy5M3bs2GG3bceOHejQoYOsPqpxcXG49957ce+99+Luu+9GdnY2zp07h2bNgnfChDN2AZbdJRmv39cTRoP9BcJoiMbr9/VEdpfkAI2MiCg0BfJuyiuvvILa2lpkZWXhq6++wqlTp5CXl4ehQ4eiZcuW+Nvf/mbbd8eOHVi8eDF+/PFHvPrqq1i9ejWmTZsGwJLV+tZbb+HgwYM4fvw43n33XTRq1Aipqano0KEDxo8fjwkTJuCTTz5BUVER8vPzsXDhQqxfv97tGMePH499+/bhb3/7G+6++25ERUXZvte+fXts2rQJO3fuxKFDh/CnP/0JpaWldo9v06YNvv32W/z888/49ddfnZZReeKJJ7BlyxYsWLAAP/74I1asWIFXXnnFbv2eO0uWLMEHH3yAw4cP48cff8Tq1athNBptt42DFWfsVCC7SzKGphuRX3QOZVU1SIy1XDA4U0dE5H+BvJvSvn17/Pe//8W8efNwzz334Ny5czAajRg9ejTmzZtnN9P0xBNP4L///S/mz5+PuLg4LFmyBFlZWQCA+Ph45ObmYsaMGTCZTOjatSvWrVuH5s2bAwCWLVuGv/71r3jiiSdw+vRptGjRAv369cPtt9/udozXX389+vTpg/z8fIfuF08//TSOHz+OrKwsxMTE4OGHH8bo0aNRUVFh22fmzJmYOHEi0tPTcenSJRQVFTk8R8+ePbFq1So888wzWLBgAZKTk/Hcc8/ZJU64Exsbi8WLF+Po0aMICwtD7969sWHDBuj1wT2npROkrtIMEZWVlTAYDKioqEBcXFygh0NERDLV1NSgqKgIaWlpiI72/HZp3sFizF9XaJdIkWyIxryR6bybQopz9bqVE5twxo6IiMgJ3k0hLWJgR0REJCJMr0NGu+aBHgaRZMF9o5mIiIgohDCwIyIiIgoSDOyIiIiIggQDOyIiIqIgwcCOiIiIKEgwsCMiIiIKEgzsiIiIgsi2bdug0+lw/vz5QA/FJ9q0aePQ8SIQx/DE4MGDMX36dJ8+BwM7IiIilTh79iweeeQRtG7dGlFRUTAajcjKysKOHTsUfZ7ly5f7tWeqPwOp6upqzJ49G+3atUN0dDQSEhJwyy234LPPPrPts2fPHjz88MN+GY+/sUAxERGRSowZMwaXL1/GihUr0LZtW5SWlmLLli0oLy8PyHguX76MyMjIgDy3p/7nf/4H3377Lf73f/8X6enpKC8vx86dO+3OYUJCQgBH6FucsSMiIhJjNgFFXwMHPrL812zy2VOdP38eX3/9NRYtWoQhQ4YgNTUVffr0wezZs3HHHXcAAH7++WfodDoUFBTYPU6n02Hbtm12x9uxYwe6deuG6Oho9OvXDwcPHgRguVX7wAMPoKKiAjqdDjqdDs8++ywAy8zaggULMGHCBMTFxdlmtXJyctChQwfExMSgbdu2mDt3Lq5cuWL3fOvWrUPv3r0RHR2NFi1a4M477wRguf144sQJPP7447bns/rmm28waNAgNGrUCK1atcKf//xnXLx40fb9srIyjBw5Eo0aNUJaWhree+89t+dx7dq1+Mtf/oLhw4ejTZs26NWrFx577DE8+OCDtn0aziAePnwYAwcORHR0NNLT07F582bodDp8+umnduf9k08+wZAhQxATE4Pu3btj165dtmOUl5dj3LhxaNmyJWJiYtC1a1d88MEHbserNAZ2REREzhSuBV7qAqy4Hfh4suW/L3WxbPeBJk2aoEmTJvj0009RW1vr9fGefPJJvPjii9izZw8SEhIwcuRIXLlyBf3798dLL72EuLg4FBcXo7i4GDNnzrQ97oUXXkD37t2xf/9+zJ07FwAQGxuL5cuXo7CwEC+//DLeeOMNLF261PaY9evX484778TLine truncated
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig0 = plt.figure()\n",
"\n",
"# Plot 1: Observations Y\n",
"ax0 = fig0.add_subplot(1, 1, 1)\n",
"\n",
"h_value = 0.25\n",
"e_value = 5.5 # meV\n",
"EIdx = np.argmin(np.abs(e-e_value))\n",
"HIdx = np.argmin(np.abs(h-h_value))\n",
"ax0.scatter(L[:,HIdx,EIdx],Y_obs[:,HIdx,EIdx],label='Observations')\n",
"ax0.scatter(L[:,HIdx,EIdx],Y_sub[:,HIdx,EIdx],label='Subtracted Signal')\n",
"ax0.scatter(L[:,HIdx,EIdx],Y_back[:,HIdx,EIdx],label='Background')\n",
"\n",
"\n",
"\n",
"ax0.set_xlabel('L [r.l.u.]')\n",
"ax0.set_ylabel('I [Arb. unit]')\n",
"ax0.legend()\n",
"\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
}
],
"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",
"version": "3.11.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}