Add MnF2 tutorial

This commit is contained in:
Jakob Lass committed 2025-04-17 10:16:35 +02:00
1 parent e018f965e7
commit 7d8b713ae0
1 file changed
+760
+760
View File
@@ -0,0 +1,760 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "8f9ac288-dfea-438b-a91c-971393102585",
"metadata": {},
"source": [
"# Simulate Neutron Scattering Data"
]
},
{
"cell_type": "markdown",
"id": "256c1739-3586-45cf-8fe3-cb89fdd1a935",
"metadata": {},
"source": [
"> **_NOTE:_** This tutorial requires the installation of the [MJOLNIR package](https://mjolnir.readthedocs.io/). However, no data is required."
]
},
{
"cell_type": "markdown",
"id": "be13ace3-32c1-4985-859f-8f85d48d2394",
"metadata": {},
"source": [
"In this tutorial we will generate synthetic data in the exact same form as measured on the CAMEA instrumnet. We will utilize the case of measuring the magnetic dispersion of single crystal MnF$_2$ as described by [Yamany et al. 2010](https://www.doi.org/10.1139/P10-081). The synthetic data will consist of 12 sample rotation scans as measured on CAMEA. \n",
"\n",
"The dispersion relation of MnF$_2$ can be approximated analytically by:\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 the specfic crystal orientation, $(H,0,L)$ defines our scattering plane meaning $K$ is zero.\n",
"\n",
"The instrument resolution/response function of CAMEA is not perfect and thus data acquired get smeard. For simplicy, this smearing will be approximated by a 1D Gaussian smearing along the energy direction.\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "d10ed368-5d6d-4507-bfee-602da4366810",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from MJOLNIR.Data import DataSet, Sample, Mask, DataFile, BackgroundModel\n",
"from MJOLNIR import TasUBlibDEG\n",
"from MJOLNIR.Geometry import Instrument\n",
"import copy\n",
"from scipy.optimize import minimize\n",
"from scipy.interpolate import CubicSpline\n"
]
},
{
"cell_type": "markdown",
"id": "f22c9d8a-76af-48df-a0e9-732412b9650b",
"metadata": {},
"source": [
"### Define needed functions for dispersion and SpinWave smearing"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "d21fe066-3051-4872-a201-5f5f884a6390",
"metadata": {},
"outputs": [],
"source": [
"S = 5/2\n",
"z1 = 2\n",
"z2 = 8\n",
"\n",
"def MnF2(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",
" 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 SpinWave(H,K,L,E):\n",
" #ani = 1.2\n",
" #amp = 0.5\n",
" sigmaE = 0.25\n",
" omega = MnF2(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": "040d416e-5e5a-409e-8756-236c46dbd93e",
"metadata": {},
"source": [
"## Generate CAMEA Data files\n",
"\n",
"For this tutorial we will generate 12 data files with the following settings\n",
"\n",
"- $E_i$ = 5, 5.13, 7, 7.13, 9, and 9.13 meV\n",
"- $2\\theta$ = -40 and -44 deg\n",
"- $A_3$ range are [ -45 , 45 ] deg at 5 and 5.13 meV\n",
"- $A_3$ range are [ -37 , 37 ] deg at 7 and 7.13 meV\n",
"- $A_3$ range are [ -29 , 29 ] deg at 9 and 9.13 meV\n",
"\n",
"The lattice parameters of MnF$_2$, needed to generate the correct scattering vectors are \n",
"\n",
"- $a$ = $b$ = 4.873 Å\n",
"- $\\alpha$ = $\\beta$ = $\\gamma$ = 3.3107 deg\n",
"\n",
"For the alignment of the sample, two peaks in the scattering plane (which is $(H,0,L)$) are needed. These are\n",
"\n",
"- r1: [1,0,+1] at $A_3$ = 35.74 deg, $A_4$ = -95.21 deg, and incoming and outoging energy equal at $E_i$ = $E_f$ = 5.0 meV\n",
"- r2: [1,0,-1] at $A_3$ =-32.35 deg, $A_4$ = -95.21 deg, and incoming and outoging energy equal at $E_i$ = $E_f$ = 5.0 meV\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "16297260-ace0-4bb7-bbf8-9e070ace8814",
"metadata": {},
"outputs": [],
"source": [
"dfs = [] # Container for generated DataFiles\n",
"\n",
"for i,ei in enumerate([5,7,9]):\n",
" df = DataFile.DataFile()\n",
"\n",
" # Define sample and initialize sample - All of this is MJOLNIR specific\n",
" s = Sample.Sample(a=4.873,b=4.873,c=3.3107,alpha=90,beta=90,gamma=90)\n",
"\n",
" r1 = np.array([1.0,0.0,1.0,35.74,-95.21,0.0,0,0.0,5.0,5.0])\n",
" r2 = np.array([-1.0,0,1.0,-32.35,-94.64,0.0,0,0.0,5.0,5.0])\n",
" s.plane_vector1 = r1\n",
" s.plane_vector2 = r2\n",
"\n",
" # Vectors defining the main axes\n",
" s.projectionVector1 = np.array([1.0,0.0,0.0])\n",
" s.projectionVector2 = np.array([0.0,0.0,1.0])\n",
" s.planeNormal = np.asarray([0.0,1.0,0.0])\n",
"\n",
" s.plane_vector1 = np.delete(s.plane_vector1,3)\n",
" s.plane_vector1[5:7] = 0.0\n",
" s.plane_vector2 = np.delete(s.plane_vector2,3)\n",
" s.plane_vector2[5:7] = 0.0\n",
" s.orientationMatrix = TasUBlibDEG.calcTasUBFromTwoReflections(s.cell, s.plane_vector1, s.plane_vector2)\n",
"\n",
" s.initialize()\n",
" s.calculateProjections()\n",
"\n",
" s.UB = s.orientationMatrix\n",
"\n",
" # Start to populate the parameters for the data file\n",
" df.A3 = np.linspace(-45-8*i,45-8*i,91)\n",
" df.A4 = np.asarray([-40])\n",
" df.Ei = np.asarray([ei])\n",
" df.Monitor = np.asarray([125000]*len(df.A3))\n",
" df.fileLocation = None\n",
" df.name = 'Test'\n",
"\n",
" \n",
" df.I = np.zeros([len(df.A3),104,1024])\n",
" df.sample = s\n",
"\n",
" # For simplicity, df2 is a deep copy of df1 with only the A4 value changed\n",
" df2 = copy.deepcopy(df)\n",
" df2.A4 = np.asarray([-44])\n",
"\n",
" df3 = copy.deepcopy(df)\n",
" df3.Ei = np.asarray([ei+0.13])\n",
"\n",
" df4 = copy.deepcopy(df2)\n",
" df4.Ei = np.asarray([ei+0.13])\n",
"\n",
" dfs.append(df)\n",
" dfs.append(df2)\n",
" \n",
" dfs.append(df3)\n",
" dfs.append(df4)\n",
"\n",
"# Collect all of the 12 data files into a DataSet\n",
"ds = DataSet.DataSet(dfs)\n",
"\n",
"# Convert utilizing the standard binning 8\n",
"ds.convertDataFile(binning=8)\n",
"\n",
"# To allow a comparison between treated and ground truth data, a compy of the initial DataSet is made\n",
"ds_calculated = copy.deepcopy(ds)\n"
]
},
{
"cell_type": "markdown",
"id": "50a47452-5217-4aac-b071-2c57c7b83a5d",
"metadata": {},
"source": [
"### Addition of noise and background\n",
"As described in the documentation, AMBER expects the background signal to be smoothly varying in all 3 directions as well as to be independent of rotations around the origin, i.e. $(0,0,0)$. In the following, such a background, mimicing real world background, is generated and added to the intensity in the DataFiles. In addition, due to the Poisson nature of counting statistics, which neutron scattering falls under, the final neutron count will be found from the calcualted intensity subject to Poisson noice"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "be9bb08d-c82c-4718-92ed-6a5dc0260ff9",
"metadata": {},
"outputs": [],
"source": [
"\n",
"def A4BG(a4,A):\n",
" gamma = 10\n",
" return 2*A*((gamma/(a4**2+gamma**2))+0.02*np.exp(-np.power(a4+120,2.0)/(2*20**2)))\n",
"\n",
"# Containers to retain the different contributions\n",
"trueForeground = []\n",
"trueRandomBackground = []\n",
"trueA3IndependentBackground = []\n",
"\n",
"for d,d_calcualted in zip(ds,ds_calculated):\n",
" trueForeground.append(np.round(SpinWave(d.h,d.k,d.l,d.energy)*(d.Monitor*d.Norm)))\n",
" trueRandomBackground.append(np.ones(d.I.shape)*0.01)\n",
" trueA3IndependentBackground.append(A4BG(d.instrumentCalibrationA4+d.A4,A=50*10).reshape(1,104,64))\n",
" d.I =np.random.poisson(trueForeground[-1])\n",
" d.I+=np.random.poisson(trueRandomBackground[-1]) # constant Background\n",
" d.I+=np.random.poisson(trueA3IndependentBackground[-1])\n",
"\n",
" d_calcualted.I =np.random.poisson(trueForeground[-1])"
]
},
{
"cell_type": "markdown",
"id": "012f2897-84af-45c1-af64-35755d735933",
"metadata": {},
"source": [
"## Utilizing some MJOLNIR features to showcase the data\n",
"Following are a couple of plots of the data using MJOLNIR functions"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "73f245e4-a0e1-4db7-bb0a-ec28730c3e57",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAHVCAYAAAATheoJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAACFkklEQVR4nO29e5hcVZX3v86pe3dX3zvduXQuJKBADOESiJABFDUgogbwhsbMgIoMEBQIiL7KTWVUGMRB8XFgAg4/CUJAXkcQ5x2u4RIkQMAk3EwgnaQ7Sacv1d11r7N/f4T00Omzvqmu6qruqv5+nqeeJ12r9jn7rL1PZdf+nrWWZYwxQgghhBBSZthj3QFCCCGEkELARQ4hhBBCyhIucgghhBBSlnCRQwghhJCyhIscQgghhJQlXOQQQgghpCzhIocQQgghZQkXOYQQQggpS7jIIYQQQkhZwkUOIYQQQsYtDz74oFiWNeR19tlnZ9XWW+C+EUIIIYTkzMaNG+WMM86Q3/zmN4PvBYPBrNpykUMIIYSQccumTZtk7ty50tLSMuK2lKsIIYQQMm7ZuHGjHHLIITm1tUq1CrnjOLJjxw4Jh8NiWdZYd4cQQggpGsYY6evrkylTpohtF2e/Ih6PSzKZzPs4xphh/28HAgEJBAKun62qqpJPf/rTsm7dOslkMvK5z31OrrvuOvH7/Qc8V8nKVTt27JDW1tax7gYhhBAyZrS1tcm0adMKfp54PC6zZlRJx65M3seqqqqS/v7+Ie9dffXVcs011wz77NatWyUajUogEJDf//73smXLFlm+fLnEYjG55ZZbDniukt3J6e3tldraWjnu5KvE6x3+AJLx6Ls76YBuM7ZuC0TSqs1y3N0YeLNDbSNKm702R7chUDu/L7djBoevrgfJgPMZ3WbCla7vW9E4OJd+c2WaalSbpzuqt6urUG1OQP8NYMDmoQWGFbUzHv3XWLLGvS+emO7jVNij2hyf3pHOI1WT2C26L884+G+q7dpJr+nHbH5ZPyEhxJVIJCKtra3S09MjNTX6999onq+mpka2rJsh1eHcd44ifY7MOvpdaWtrk+rq6sH3tZ0cEZGuri6pq6sb3P1ZvXq1fOUrX5H+/n7xePTvOZES3snZd7Feb9B9keMF/5uAL3i0yPF6R77I8dpggSBoIZPjIge1sw+8tefeDlwDWMjARY7H/ZiWjVYI+iLHcpkD+/B4cmvneMfPIsfxKYucNPCxT7/5M369IzYIWrAr9PMFqvRFdDVYcNnv+6IjhIyMYj+uUVm195Urmfe+H6urq4cschD19fVD/j700EMlHo9LV1eXNDU1wbYlu8jZh512xHb5j92K6//TOF79y9iA3RXHr/8nFNgdU23quWL6roUBuqeFVq7Api3EDoQZ0H+9C1gUWsqqXETE2tXlfq54Qm9Tq98QnnfAjhnwiSeV0m1oBysNtmy9+vlMAMw9sPvo63W3ZSr1has3pvcxUa/3Y/Ia1SSJWn3n66HXj1dtq6s+rNr+fql+PkLIxObRRx+Vc845R9ra2qSiYu/3zyuvvCINDQ0HXOCIMLqKEEIIIVniiMn7NRKOP/54CYVC8rWvfU3eeOMNeeSRR2TFihVyxRVXZNW+5HdyCCGEEFIcHHFyfphiX/uREA6H5dFHH5Vvfetbcswxx0g4HJbzzz9fVqxYkVX7kl/kWGkjlsvKMFGnb8enqvQNLCsD5CrwLE96mvs2fhV4ZkIytarJHgAP4AIN1vjBcxiVuvxiQAiid0+/akPnM0AmcoLuUy9Vo/exf4o+puFturzXfYh+zJrNulzVfYh+Pg+IokyGdZtHV+MkCvJcpRrdnwfz9ug+Tlfrc696SkS1fWL6G6ptQeVm1fa5ql7VhqFeRQjROfzww+W///u/c2pb8oscQgghhBSHjDGSySMoO5+2ucBFDiGEEEKyIpfnavZvX0xKfpHj7Uu4BrPYKX2r3k7rESnpYG7PYvsj7pEsSAby9PTpBwSRVyjfjdWrt/MM6JExKCLICesxxXY/0F8yIELM5x6DGNzao7bxd+rj5oR0W83f9X74+vW0APWb9JvRk9Ajl1JV+m2VrtD9XAkCxFJKO1vvvvQcAqTEN+tU2/1H6IlyHvQfodp+09yp2v5tzr2q7YMderp2u+VN1UYIKT6OGMmU0CKH0VWEEEIIKUtKfieHEEIIIcWBclWRsdKOWC6Zde0BXaIIgAefrFogiYBkgKIEPFkpkDQuDbQGUJ7B9OnRThbI0iugHZKr7ChIdIgeIqvU5TG7fY+7AZRusBMh3barR7V5QyCFL/AzvDlA5JgX+FJAwj9EqsG9DAaSZUO79StwAvpcDu0Bci4oidLeqNeSO33HRarNDOj9fOcbqokQMgaU2oPHlKsIIYQQUpaU/E4OIYQQQoqDI7lXVtzXvpiU/CLH+D1iPMMvw0roSd6sFJATQC2mdFC3pSrdI57shJ6Izmf0rHEWiLyyQrpsAyUwVIsJ1IxCEpIAKcipBvJSLlXWQaSaAYkOraTe/+T0etXmjeg+6Zutj50WaScikgSFKvungtpoPe5bvMmwPicT+qVJulLfMjateoSe16/Pr2l1Part+oP+oNoWgmrvItllNSWEFIdMntFV+bTNBcpVhBBCCClLSn4nhxBCCCHFIWP2vvJpX0zKdpGD6ybpNk9clxqcWr1dqspdNsjM0uWc6i2qSbxefZPN8evDZoFaWXZSlxoMkJZSdfo1eAb0Y1rgKfroIU2u7we6dYkoXaEnQUxXgminAZC4D8hH8Xr9fLvn6+NT83ddQupv1W2eqGqSrrnu72fCuv89NXqEYVWlLkm1hHWp9Mi6bartx5NeVW1l/FVDyISCz+QQQgghpCxxxJKMljMly/bFhM/kEEIIIaQsKfmdHLu9S2zbJXmZD11atWrxDOitqkDitUSDewI1Twwkm+vVpRl7T69uS4BIKJQMEEQ0IWkp0AVqXiEcXSaq6HKPhjIgkZ4vrksz/t0gAszSfzn4Irq/MkAeO2i1LgU5IFqoZjOQSmP69aXC7vMr1qSfyxfVI87SICpuZ4MelvXHihmq7aHACaot8UE9oeQnPrBJtf1KWNeKkPGEY/a+8mlfTEp+kUMIIYSQ4pDJU67Kp20uUK4ihBBCSFlS+js5jvKsN0huZ4NEe+LX6/bYAV2+8PW7yyXpihzrGIGkhDApHpKyABaIRpMBoOGBvkDprN/9mBaod2X1gnFDSRDDVeCYej0v2wY+AdiodhWq2wIkVnvAfVw9cRAVF9bnqz8NorwSoI/AtOso/bqDG/V+vt7crNqO7zxLtT3foveFEFIYSm0np/QXOYQQQggpCo6xxDF5RFfl0TYXKFcRQgghpCwp/Z0cTa4CshO0AZkIRSD1TXOXBvx9ept0ld4PXwLINmk9ksiqqlRtktAjgsQD1rtIYkGPyqNj+t3lCwMik1DCQrtPT26XbtLrTCHshC6BRafq42MndZ94QbLJvlZQf0s5pL9Pj5hDtdbidfrYJPXgQ3HArWOASnfoJ99SbZ9o2KDavl7Trh9UbgA2QkghoFxFCCGEkLIkI7Zk8hCBQLKPgsBFDiGEEEKywuT5TI4p8jM5pb/IqawQsYdv85tKsPUPZAinUk+ShmonhbrcZYOMXx/QNIh+8faC+lQgkZuAyB6nukY/Zq5V01KgLlSTLp1ZisxlgcImjhdEBKGIM/CjwwkAf4HzoX5q1yYikgrr4+pJ6e0SYfeLSILaW1E9aEmsuXqkmsejX9z02h7Vdmz9u6rtq3VrVdtsrx79hnA6mCiQEIIp/UUOIYQQQooCn8khhBBCSFmSMbZk0Pb4AduPYmeyoOQXOaY3IsYaHvJh9enb+CajSyx2So/E8cdSqk1LyhadDKQlMNp2Qj+XxPRIIiRX2RE98Z0An4hPl9VQPSz/rj16O79yTFRzDEXFgfpUnt0gqgwcE0V6+Xfq/nKC+jHTtfp8CO4Ctbnq3eVXb0zvR92buk+c54AsG9K/wHoq9dCr/wrrda1WtZ6o2iYfu0O1fbBml2r79bTnVBshhIiUwSKHEEIIIcXBEUucPKKrHJQ2vQBwkUMIIYSQrOAzOcXGslwT+DmxmN4E1VtK6jKRZentbL+7Kz2pkSd4ExExPl12QkkJZUC/bgPqO1mgzpRJRvXzASxNkhLR5THgf8mAkCaUeBAkT5QU8IkDpEaQWNEG5/P3g/EB/gr1u0uUBtUH8+o+SYd1SQ1Fh3mjwJeW3n8b1MraunmSfsyDdNMHtnxVtb31Ob0dIWTiUPqLHEIIIYQUhfwfPKZcRQghhJBxyN5ncvIo0Em5amRYHo9Y9nBpxwqDBGNIrkJRRkDuEY/7wFkgwRuKYknX6nWafNGE3g9ES4NuG9CPaRldJjIVej9lAMhcWkJDVAsL/AJIT9KjflANKiuuy06ZsH5tVlr3iRPSbyvHD8a8Qpcok1Xu7bwx3Sf9U/Xjhfbo/e+Zo/fR36OapPdD4P4Aw/rwqT9XbSnwizHi6FKwyP8BNkJIrjh5lnUo9oPHrEJOCCGEkLKk5HdyCCGEEFIc+ExOsfF6RWyXy0DSEsKjb/EjucTudZdmvFV6FEsGJI1LVeuRKj49dxqU26wekAwQJNNzi14bbIYkKXBMTQrydEbUNk69LkF6d/bq/UD1vECNM0+fHgmVnKQnjfT16O1SzXo9r+BOEBGYdpf3jCKTiohUb9XvgXRA/5Kq2gZkwaB+vop39a+T6GRdHvvRjk+qNi8oEnZ4lX4jnKDUtWJNK0LywxG7pPLkUK4ihBBCSFlS+js5hBBCCCkKGWNJxuSRDDCPtrlQ+oscxxERly3tAIi8SIJaRkjm0iKCRNSILTupb7ejZIBWCiS+S4GEeUFw3SByyUR12clCvkSRakAm8mjJ7ap0H9s7u/VzAUy9HnmFJEgrCaQ/FAUGhs4fAckmQcSWv8d9ziYa9bFxvPqXiS+mnytVpY9boFdLine truncated
"text/plain": [
"<Figure size 640x480 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAF2CAYAAAB3QMMiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABoX0lEQVR4nO29fXhU5Z3//z4zmUxm8pyACU8NqGXLSr3SCN0W+VojUh9at0qxFJYaFbdCXZFaQ2m3++XBh9UWUbv4cFmuAtaVh1XQ/bVgpF+KigorWtiWUPpAgBAIkASSkGQmkznn90cy534fMofMTJLJTPJ5Xde5OHPmnPvc58zk5p73+XzeH80wDAOCIAhC0uEY6A4IgiAIsSEDuCAIQpIiA7ggCEKSIgO4IAhCkiIDuCAIQpIiA7ggCEKSIgO4IAhCkiIDuCAIQpIiA7ggCEKSIgO4IAhCHNi6dSs0TbMsM2fO7FWbKX3UN0EQBOESVFZW4rbbbsPLL79sbktLS+tVmzKAC4IgxIFDhw5h4sSJKCws7LM2RUIRBEGIA5WVlRg/fnyftqmJG6E9uq7j5MmTyMzMhKZpA90dQRgQDMNAc3MzRo4cCYdjcM35fD4f2tvbYz7eMIxuY4Pb7Ybb7e62X0ZGBv7xH/8Rn3zyCYLBIO68806sWLECqampMZ9fJJRLcPLkSYwZM2agu5F0fNU9x1x/x//aAPZE6Euqq6sxevToge5Gn+Hz+TCuKAO1Z4Ixt5GRkYELFy5Yti1duhTLli2zbDt+/DhaW1vhdruxefNmVFVVYeHChWhra8Nzzz0X8/llAL8EmZmZADq/uFlZWf1+voaGBowbNw5VVVXIy8vr9/P1F3detsBcb2xsjOiYwXLtsZDo197U1IQxY8aYfw+Dhfb2dtSeCaLqkyJkZUb/y6KpWce4a451Gx8unn0DQFFREerr65GbmwtN01BcXAxd1zF37lysWrUKTqczpmuQAfwShH4aZWVlxWUADwQCADr/44jH+eJBpNcxGK89UpLl2gerjJiV6YhpADePj3B8uPg/5wkTJsDn86GhoQHDhw+P6dwygAt9ws1Z96gXKeprNe2GfwcA/L+dP4p3lwQhIoKGjmAMTwKDhh7xvhUVFZgzZw6qq6vh9XoBAPv370d+fn7MgzcgUSiCIAxxdBgxL5EyZcoUeDwe3HfffTh8+DC2b9+O8vJyLF68uFd9lxl4AuF2uzFr1qywGtpgR659aF57IqBDR+RzaetxkZKZmYmKigosWrQIkyZNQmZmJu6//36Ul5fHcGaFDOAJhNvtxuzZs5PmD5llEyPQYa47sjLMdeeFQERtJdu19yVD+doTgaBhIBhDNHW0x1x11VXYsWNH1Oe5FDKAC4IwpIlWDuHjBhrRwAVBEJIUmYELUXHLBBVNonk9ap1DzDLSzVVnbUNc+iUIsaLDQDBJZ+AygAuCMKRJZglFBnBBEIY08XqI2R/IAC70yLV3Pm2ue/NVhIkjTZnwaO0qCkVrbDbXjeZOn4gvzV1lbtvz6sP90k9BiAW9a4nluIFGBnBBEIY0wRg18FiO6WskCkUQBCFJkRm4YHJLgXIRRH6uuWpMzDfXz31ORZ6kn3SZ696j59WxAZW8o2V3mvzkVvy5D3sqCH1H0ECMXih935dokQFcEIQhjWjggiAISYoODUFEb5Wrx3BMXyMD+BDnq65vm+uODBVholGIVNbh8+Z6yzjlexz0qEco7QVk9k/rqVVnOlfcsZeNEoT+RDc6l1iOG2hkABcEYUgTjHEGHssxfY1EoQiCICQpMgMfgrBsolFFbPY2QcP5sPt4Xap2n+8yr7neMkpZoWYcbzPX9WHZAABHq6r8/aV/UolBe/7zB9F2XxD6lGSegcsALgjCkEY3NOhGDA8xYzimr0kKCcXn82HevHnIycnBiBEj8PTTT4fdb926ddA0rdvicKjL3Lp1a7f3Z86cGa9LEQQhwQjNwGNZBpqkmIGXl5dj37592LlzJ44dO4aysjIUFRV1G3hnzZqFm2++2XwdCARwww034Otf/7q5rbKyErfddhtefvllc1taWlr/X0QCoaWoBBwHWb+i3aZ6DlnFam1qn7TTLeZ6SquSUIJp6mvlz+/c7uhQj+zdTYkQQSsInQThQDCGuWywH/oSLQk/gLe0tGDNmjXYvn07SkpKUFJSgoMHD2L16tXdBnCPxwOPR+m4//7v/w7DMPDkk0+a2w4dOoSJEyeisLAwbtcgCELiYsQooRgiofTMgQMHEAgEMGXKFHPb1KlTsXfvXui6/UyuoaEBTz31FJ588klLrcHKykqMHz++X/ssCIIQDxJ+Bn7q1CkMGzYMqRQJUVBQAJ/Ph/r6egwfPjzscS+++CJGjhxpmaUbhoHDhw+joqICTzzxBILBIO68806sWLHC0v5gxJKwk50VfqcUFWFiXJZnrmtnVFUdrbmFtvvMdVeh8kthycWR1xmp0pGuvmoa/b979cPPmOuPfuXyS12CIPQLEoXSC9ra2lBTUxP2vREjRqC1tbVbte7Qa7/fH/Y4wzCwZs0aLF682LL9+PHjZnubN29GVVUVFi5ciLa2Njz33HO2fWxoaECADJrcbne/VBAPnYPPNdQYiteeaJ+73++3/G01NzdfYu/kJ2g4EDRi0MAlExPYu3cvSktLw763detWpKWldRuoQ6+9Xm+4w7Bv3z6cOHEC3/72ty3bi4qKUF9fj9zcXGiahuLiYui6jrlz52LVqlVwOp1h2xs3bpzl9axZszB79uyIri8WduzY0W9tJzpy7QPPhg0bsGnTpoHuRtzQoUGPQU2WkmoArr/+ehiXKE304Ycfoq6uDh0dHUhJ6exubW0tPB4PcnJywh7z9ttv47rrrkNubm639/Ly8iyvJ0yYAJ/Ph4aGBls5pqqqCpmZyt+jP2fgO3bswPTp0+FyuXo+oAe+5v2Ouc4+JxZS6CvA/4EdVb+KDJfaJ1RhB7BKMVptvTo2Td2b1PNd+19ZoE7Tpqr3jNyp2lsx/X/M9T03Lgvf30FIX3/uvWXatGl4/vnnzdfNzc3dJjGDCZFQ+pHi4mK4XC7s2bMHU6dOBQDs3r0bkydPtsR3M3v37sW1117bbXtFRQXmzJmD6upqc/a+f/9+5Ofn2w7eQOegn5Vloxv3Ay6XKyH+kAeSoXj9ifK5u1wuZNB/+InQp/4kdglFZuA94vV6UVZWhvnz52Pt2rWoqanBypUrsXbtWgCd+rTT6UR2drZ5zB//+EfMnTu3W1tTpkyBx+PBfffdh6VLl+LIkSMoLy/vppUnM5Y0eYr31jjeu0PNgBFU0ax6Y5O57himfqkYVONSoz9m3ZJuT+fih6HDcgAAqYdOqHPyf7y0nvOUmqXjFgiC0AMJP4ADwKpVq7BgwQKUlpYiOzsby5cvx4wZMwAAM2bMwNixY7Fu3Tpz/9OnT4eVTzIzM1FRUYFFixZh0qRJyMzMxP3334/y8vJ4XYogCAlGpwYufuD9htfrxfr167F+/fpu7+3atavbtra2tm7bQlx11VUJ87BIEISBR48xE1MeYgp9gkU2oYerjhFKkjDqVCw3nDZf1s+OVfv/7XjYXQxKnnIUDFNvBCm4myQao6pLOqEHpFqukrvgUxFGqYfVg9PpUx8313fs/tfw/RWEPkA0cEEQhCRFh0PCCAVBEJKRoKEhGIOvSSzH9DUygA8CNJInHHnq4a1xpk7txBIGx35nqfh2409H1Dq3T3HgjkwVXqaTLOMgEzHDIJklFP1CPzf5OI5q0ahuZsrZwZ39Jwh9gQzggiAMaWK3kxUJRRAEYUDRDQf0GB5i6vIQU4gVO3dBTnXXSB4ByROgffSTtWo7RadwHUyOKtE5lX44ORBSNAn8KjnIaFOOheZxed1j9AEAbOZ0QbkeTp/ymLm+48OfhD9WEGJEZuCCIAhJio7YHkgmQl2phC/oIAiCIIRHZuBJCifsaBQBYpEhKOmGI1I0iiThKBD2RWHZRPOq9lmWMc4qB0KDzuvIzVHttLd3/kvnMVpbVXvsgMjRMdT3lL+pBJ9/+M4qc33vrx6GIPSW2OPAB37+KwO4IAhDmtgzMWUAFwRBGFDEzEqIGzel3wUAcHA1opBMAcCwk0EKlN+5fuKU2sehvoQcecL2sAbbz54+q9Y5wSdNyTJGU/ckHI3qZFokH4YiWfg6NE+auZ6393T4YwUhRmQGLgiCkKTEHkY48AP4wPdAEARBiAmZgScBN3motmVWp1RhkRh0SsBJU3KDJQGHEnZYHmGPFIOkGEZLJ7kmk+pjnlcVfHQfJRCxFJPdFbXCCT0kp7D8gwx1Hq2domnoWtkW9+bi/2uuv71/Rdi+C0JP6IYGPZY4cDGzEgRBGFhiL+gw8AKGDOCCIAxpYvdCkQFciAAtTSXtGF2RGlpejtqBquEYVJiYCxmzrGGRSkiesPMoMc43hu8XSSsWi1qKbMGFrqSdbPJlITnF4pVCsglby4IlGU78qVXJSdNu+Hdz/f/t/FHY/gpCOILQEIwhJDCWY/oaGcAFQRjSJPMMfOB7IAiCIMSEzMATlFDCDnBR0k5IWrig/EQMv0qAiShhh31UOAqFPEqgK6tM9kLhCBKLtSzJNRZPk1Db7MXC3iosBZGEYrkm9l228W5J/cMxc/2W8T8017f/+alufREEJojY5JBgz7v0OzIDFwRhSBOSUGJZYuFrX/sa7r777j7puwzggiAMaUKp9LEs0bJx40Zs27atz/ouEkoC8TVvZ8LO89gA5zBV7cYgjxDzhx5bwlJijm3CDkOSBPucaDlU2YeSdCxRLhy1QraxRouNRWxXpIpG57Gcs7VN7cvSCnu6sMxDlXo4OgcpdE6uDiQIPWDEaGZlRHlMQ0MDysvLMXny5KjPZYcM4IIgDGniZWb1yCOP4Dvf+Q5OnjwZ9bnskAFcEAShFzQ1NVleu91uuClQAAB27tyJ9957D3/4wx+wYMGCPjuLine truncated
"text/plain": [
"<Figure size 700x400 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"\n",
"# Generate 2D cut in momentum transfer and energy\n",
"q1 = np.array([-1.5, 0.0, 0]) # r.l.u.\n",
"q2 = np.array([0, 0.0, 0.0]) # r.l.u.\n",
"EMin = 0.0 # meV\n",
"EMax = 6.0 # meV\n",
"dE = 0.05 # meV\n",
"width = 0.03 # 1/Å orthogonal to cutting direction\n",
"minPixel = 0.03 # 1/Å bin size along cutting direction\n",
"\n",
"ax, *d = ds.plotCutQE(q1=q1, q2=q2, EMin=EMin, EMax=EMax, dE = dE, \n",
" width=width, minPixel=minPixel, vmin=0.0, vmax=5,\n",
" colorbar = True)\n",
"\n",
"# Generate a constant energy cut\n",
"EMin = 4.5 # meV\n",
"EMax = 4.55 # meV\n",
"dqx = 0.03 # 1/Å bin size\n",
"dqy = 0.03 # 1/Å bin size\n",
"\n",
"\n",
"ax, *d = ds.plotQPlane(EMin=EMin,EMax=EMax, xBinTolerance=dqx, yBinTolerance=dqy, \n",
" vmin=0.0, vmax=5,colorbar = True)"
]
},
{
"cell_type": "markdown",
"id": "6bd0bbba-3255-45ac-861a-236085bc28b3",
"metadata": {},
"source": [
"## Initialize AMBER background"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "f7ea4236-8f1e-4edd-90d0-8b90ad7f1469",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Anaconda\\envs\\python311\\Lib\\site-packages\\MJOLNIR\\Data\\BackgroundModel.py:309: RuntimeWarning: invalid value encountered in divide\n",
" self.I = np.divide(self.data[0] * self.data[-1], self.data[1] * self.data[2])\n"
]
}
],
"source": [
"AB = BackgroundModel.AMBERBackground(ds,dQx = 0.03, dQy = 0.03, dE = 0.05, beta=100)\n",
"AB.set_radial_bins()"
]
},
{
"cell_type": "markdown",
"id": "bf3bb3ac-b1c3-4e9f-9926-87d43d6ac4a4",
"metadata": {},
"source": [
"## Minimize RMSQ for beta\n",
"\n",
"> **_NOTE:_** Running the cross-normalization will take quite some time. > 5-10 minutes!\n",
"\n",
"To perform a cross validation of the $\\beta$ value, on needs to set the qantile level where only background is left. That is, byt masking out higher intensity regions one can assume that the original minimization problem only contain the data and the background, i.e.\n",
"\n",
"$\n",
"\\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}\n",
"$\n",
"\n",
"becomes\n",
"\n",
"$\n",
"\\min_{b} \\frac{1}{2}\\lVert Y-\\mathcal{R}b\\rVert_{2}^2+\\frac{\\beta}{2} \\mathrm{Tr} \\left( b^T L_{b} b \\right) \n",
"$\n",
"\n",
"For the above data set, a value of 0.75 (the default setting) is suitable. \n",
"\n",
"### Minimization methods\n",
"\n",
"One can do two things to find the optimal $\\beta$ value\n",
"\n",
"1. Calculate the RMSE across multiple decades and interpolate\n",
"2. Perform a minimization\n",
"\n",
"\n",
"It is adviced to first perform a coarse calculation for then to perform the minimization as the latter can be quite time-consuming."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "af82019c-84fd-47ff-844c-81c2becb4a7f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"q=0.75, beta=array([ 1, 3, 7, 10, 30, 70, 100, 300, 700, 1000]), l=0.2665642114325497, mu=0.2427590851339083, n_epochs=20, verbose=False\n",
"Test - ( 1 )\n",
"RMSE - ( 0.2665642114325497 1 0.2427590851339083 ) : 0.10329160145964443\n",
"Test - ( 3 )\n",
"RMSE - ( 0.2665642114325497 3 0.2427590851339083 ) : 0.1028468402733826\n",
"Test - ( 7 )\n",
"RMSE - ( 0.2665642114325497 7 0.2427590851339083 ) : 0.10225770010882033\n",
"Test - ( 10 )\n",
"RMSE - ( 0.2665642114325497 10 0.2427590851339083 ) : 0.1019720100920882\n",
"Test - ( 30 )\n",
"RMSE - ( 0.2665642114325497 30 0.2427590851339083 ) : 0.10156595368428951\n",
"Test - ( 70 )\n",
"RMSE - ( 0.2665642114325497 70 0.2427590851339083 ) : 0.10227250242411776\n",
"Test - ( 100 )\n",
"RMSE - ( 0.2665642114325497 100 0.2427590851339083 ) : 0.10292055384172495\n",
"Test - ( 300 )\n",
"RMSE - ( 0.2665642114325497 300 0.2427590851339083 ) : 0.10781518872596527\n",
"Test - ( 700 )\n",
"RMSE - ( 0.2665642114325497 700 0.2427590851339083 ) : 0.11480259052941344\n",
"Test - ( 1000 )\n",
"RMSE - ( 0.2665642114325497 1000 0.2427590851339083 ) : 0.11836036118710512\n"
]
}
],
"source": [
"# We define a list of reasonable beta values distributed over 3 orders of magnitude\n",
"\n",
"beta = [1,3,7,10,30,70,100,300,700,1000]\n",
"RMSE = AB.cross_validation(q=0.75, beta = beta,verbose=False,n_epochs=20)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "ef25cca1-bc0b-4842-a7a3-3dbfd173d052",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x211ecd07490>"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAGwCAYAAAC5ACFFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABWZklEQVR4nO3de1xUdf4/8NcMt2EQBgGFAVEBL0iKiGJoeCtrbUs3zKxtVSw1NVddTSu7Gb9yXTdv9bXM1kLTNM1Ny1bNTK1QAW8QN++oiIwKxAy34TLz+f3BMuvERZSBMwOv5+Mxj93zObf3OebMy88553NkQggBIiIiIqqTXOoCiIiIiKwZwxIRERFRAxiWiIiIiBrAsERERETUAIYlIiIiogYwLBERERE1gGGJiIiIqAH2UhfQGhiNRly/fh2urq6QyWRSl0NERESNIIRAUVERfH19IZfX33/EsGQB169fh7+/v9RlEBER0T3Izs5Gp06d6p3PsGQBrq6uAKpPtpubm8TVEBERUWPodDr4+/ubfsfrw7BkATWX3tzc3BiWiIiIbMydbqHhDd5EREREDWBYIiIiImoAwxIRERFRAxiWiIiIiBrAsERERETUAIYlIiIiogYwLBERERE1gGGJiIiIqAEMS0REREQN4AjeREREZJUMRoGkrALcLNKjo6sCAwM8YCdv+RfWMywRERGR1dmXlovY3RnI1epNbWqVAotHh2BUb3WL1sLLcERERGRV9qXlYubmU2ZBCQA0Wj1mbj6FfWm5LVoPwxIRERFZDYNRIHZ3BkQd82raYndnwGCsa4nmwbBEREREViMpq6BWj9LtBIBcrR5JWQUtVhPDEhEREVmNm0X1B6V7Wc4SeIM3ERERSc5oFPg+XYPMXF2jlu/oqmjmiv6HYYmIiIgkI4TAj5k3sfKHc8jI1UHl7ABvNyfc1JXXed+SDICPqnoYgZbCsEREREQtTgiBn8/nYeUP55CSXQgAcHG0w6RBXRDUwQXztqVABpgFppoRlhaPDmnR8ZYYloiIiKhFZebq8NY3aTh++TcAgLODHSYN7oLpQ4Pg4eIIAFA42NUaZ8lHonGWGJaIiIioRTnYyXHyym9wtJdjYmQXzBgWhA6uTmbLjOqtxsMhPhzBm4iIiFq/5OxCnLhcgKlDAgEA3Tq2w3vj+uKBbl7wUdV/o7adXIZBQZ4tVWa9GJaIiIioWaRf12LVD+dwIPMm5DJgRHBHBHVoBwB4sn8niatrPIYlIiIisqhzN4qw6odz2JumAQDIZUB0v05wdrCTuLJ7w7BEREREFnFTp8e7/8nE7l+vQwhAJgPG9PXFnIe6m3qUbBHDEhEREVmEwtEOP527BSGAP/bxwd9G9kAPb1epy2oyhiUiIiK6J9kFpdh5OgezH+wGmUwGN4UD/jG2Dzp7KnGfr0rq8ixG8nfD6fV6TJkyBe7u7lCr1VixYsUd14mPj0dgYGCd85YsWYLJkyebpg8fPgyZTFbn5+rVqwCAuXPn1pq3Zs0aixwfERFRa5NTWIZFX6dixPLDWPnDORw6e9M079E+6lYVlAAr6FlauHAhTpw4gYMHD+LKlSuIiYlBly5dMG7cuDqXT01Nxbhx46BQ1H7UcOvWrVi8eDEmTJhgahs8eDByc3PNlhs/fjw8PT3RuXNnAEBGRgaWLl1qFrLc3NwscHREREStR662DB8euoBtx7NRaageW3tId68WfU+bFCQNSyUlJVi/fj327t2L8PBwhIeHIz09HWvWrKkzLK1btw4LFixAYGAgtFqtqb2qqgqzZ8/Ghg0bEBQUZLaOo6MjfHx8TNNbt25Famoqzp8/b2rLzMzEwoULzZYjIiKiaqUVVfjnvrPYkngVFQYjAOCBbp6YN7IHBnRtuXe0SUXSy3ApKSmorKzE4MGDTW1RUVFITEyE0WistfzevXuxceNGzJs3z6y9uLgYv/76KxITEzFo0KB691dZWYk33ngDr7/+Ory8vAAAOp0OOTk56NGjh4WOioiIqHVR2Nvh6MU8VBiMGBjggS9fiMQXUyPbRFACJO5Zys3NhZeXFxwdHU1t3t7e0Ov1yM/PR4cOHcyW37VrFwBgw4YNZu3u7u44cuTIHfe3fft2FBYWYtasWaa2zMxMyGQyLFmyBHv37oWnpyfmz5+PmJiYerdTXl6O8vJy07ROp7vjvomIiGxFQUkFNhy9jJnDguDsaAe5XIa3x9wHIYDBQZ6QyVr+lSNSkjQslZaWwsnJ/F0wNdO3hxFL+eSTTzB16lQ4Ozub2s6cOQOZTIbg4GDMnj0bP/30E1544QW4ubkhOjq6zu0sXboUsbGxFq+PiIhISoWlFfjXL5ew4chllFQY4Opkj2lDqx+oGhzkJXF10pE0LCkUilqhqGZaqVRadF83b97EL7/8Uuspt0mTJmH06NHw8KjuSgwNDcW5c+ewdu3aesPSokWLMH/+fNO0TqeDv7+/ReslIiJqKdqySnwan4XP4rNQXF4FAOjt54Zgte2PkWQJkoYlPz8/5OXloaqqCvb21aVoNBo4OzvD3d3dovv6/vvvERAQgD59+pi1y2QyU1Cq0atXLxw8eLDebTk5OdXqESMiIrI1QgisOXgBn/xyCUX66pDUS+2GeSO74+EQ7zZ3ua0+kt7gHRYWBgcHByQkJJja4uPjERERAbncsqUlJibigQceqNX+1ltvYeTIkWZtycnJCA4Otuj+iYiIrI1MJkP6dR2K9FXo4d0Oa/8Sjv/MjsIj9/kwKN1G0p4lpVKJmJgYzJgxA3FxccjJycHy5csRFxcHoLqXSaVSmd1jdK/S0tIwatSoWu2jR4/G0qVLsXz5ckRHR2P//v34/PPPcejQoSbvk4iIyJqUVlRh07EreCxUjU7tq293eemRHngsVI3H+qghlzMg1UXyEbxXrlyJ/v37Y8SIEZg1axZiY2MxduxYAIBarca2bdsssp8bN26gffv2tdojIiKwY8cObNq0Cb1798YHH3yALVu2NDgEARERkS3RVxqw/pdLGPrPQ1i69ww+PHTBNK+7tytG9/VlUGqATAghpC7C1ul0OqhUKmi1Wo78TUREVkNfacC249n48NAF3CyqfoCqs4cS8x7ujuh+nSSuTnqN/f2W/HUnREREZHlfncjGyh/OIVerBwD4uTtjzkPdMDa8ExzsJL+wZFMYloiIiFqhy/klyNXqoVYp8NcHu+Gp/v5wtGdIuhcMS0RERDauymDEztM5COzggv5dqofDmTYkEN5uCjwd4Q8nezuJK7RtDEtEREQ2ymAU+DYlB+8fOI/L+aUY2NUD26ZHQiaTwV3piEmDukpdYqvAsERERGRjDEaB/6TmYvWBc7h0qwQA4OHiiJEhHWEUgB0fbLMohiUiIiIb8tO5W1jynwycu1EMAHBXOuCFoYGIGdQVLk78WW8OPKtEREQ2JK+oHOduFMNNYY9pQwIx+YGucFU4SF1Wq8awREREZKWEEIi/kAd9pREPh3gDAP4U5ovfSivw1AB/qJwZkloCwxIREZEVOn65AO99fxZJWQVQqxQY0t0LCgc72NvJMXVIoNTltSkMS0RERFYkJbsQK344h5/P3QIAONrL8WhvNSoNRigcOASAFBiWiIiIrMDFW8X4x94z+CHjBgDAXi7D+Ah/zH6wG9Sqpr9Qnu4dwxIREZEV0JVV4oeMG5DLgCf6+eFvD/VAZ0+l1GURGJaIiIgkkV1QiuTsQozu6wsA6Ne5PRY9GoyHenVEt46uEldHt2NYIiIiakEarR5rDp3HtuPZkMlkiOjqAR+VAgAwfViQxNVRXRiWiIiIWkBecTnWHr6ITQlXUFFlBAAM6e6JskqDxJXRnTAsERERNSNtWSXW/XQRG45eRmlFdTAa2NUDLz3SA/cHekpcHTUGwxIREVEzKq804NP4LJRXGdG3kwovPdITQ7p7QSbjC9xsBcMSERGRBZVVGPDjmRt4PLT6xu2ObgosejQYvu7OeDjEmyHJBjEsERERWUClwYjtJ7LxwY/ncUNXDrVKgf5dPAAAkx8IkLg6agqGJSIioiYwGgW+S83Fyv1ncTm/FADg5+6M4nLeuN1aMCwRERHdAyEEDp+7hff2nUVGrg4A4OniiL8+2A3P3t8ZTvZ8NUlrwbBERER0DyoMRrz+dSqua/VwdbLHtKGBeD4qAO2c+NPa2vBPlIiIqJHO3ShCoJcL7O3kcLK3w8JRPZFxXYeZw7vBw8VR6vKomTAsERER3cHV/FKsOnAOu5JzsOzJUIwf4A8AiO7XCdH9JC6Omh3DEhERUT1uFumx5uAFbE26ikqDAABkXNdJXBW1NIYlIiKi36kZdTvuyGXT60iGdPfCy38IRp9OKomro5bGsERERPQ7L21PxoHMmwCAMH93vDyqJwYHeUlcFUmFYYmIiNq8KoMRVUYBhUP14/4zhgXhakEpXnqkJx7hqNttHsMSERG1WUIIHMi8iWX7zmBkL2+8+mgwAGBAVw/smzsUcjlDEjEsERFRG3X66m9YuucMki4XAACK9VWY93B302CSDEpUg2GJiIjalKy8Erz3/RnsSdUAAJzs5Xg+KgAzhgVx1G2qE8MSERG1Gd8k5+Cl7SmoMgrIZMC48E6Y/0gPqFXOUpdGVoxhiYiI2oyBAR6wk8swtEcHvDIqGD19XKUuiWwAwxIREbVKVQYjtp3IRsZ1HZZE9wEAqFXO+PGlYejUXilxdWRLGJaIiKhVEULg+/Qb+Of3Z3DpVgkA4Mn+nRDeuT0AMCjRXWNYIiKiVuPklQL8fc8ZnLzyGwDAw8URcx7sht6+HHWb7h3DEhER2bz84nK89U06/pOaCwBQOMgxNSoQ04cFwlXhIHF1ZOsYloiIyOYpHe1x8spvkMuAp/r7Y/4jPeDtppC6LGol5FIXoNfrMWXKFLi7u0OtVmPFihV3XCc+Ph6BgYF1zluyZAkmT55s1nb69GnIZDKzz4ABA0zzs7KyMHLkSLi4uCALine truncated
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAG2CAYAAABSy6M+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABq/UlEQVR4nO3deVhUZf8G8HsYlmGdkX3YNzfEDUWBUHG3XEojW1VKK9HSN9PMfpXSW9miqb2mmRaa5p6lmFvmUiiLG4bgLiIgiIAOArLNnN8fxOQEJOjAYbk/18VV5znLfGeOOjfPec5zJIIgCCAiIiIivTAQuwAiIiKiloThioiIiEiPGK6IiIiI9IjhioiIiEiPGK6IiIiI9IjhioiIiEiPGK6IiIiI9IjhioiIiEiPGK6IiIiI9EjUcFVSUoKJEydCoVBAqVRi4cKF990nJiYGXl5e1do3bNgAb29vmJmZYfTo0cjNza1x/6lTpyI0NFSnLTU1FYMGDYK5uTl8fX2xb9++B3o/RERERKKGq1mzZuH48eM4cOAAli1bhsjISGzdurXW7ZOSkhAWFgaNRqPTnpCQgIkTJ2Lu3LmIi4vDrVu3EB4eXm3/o0ePYvny5TptgiDgiSeegKOjI44fP45x48Zh9OjRuHbtml7eIxEREbUuErGeLVhUVARbW1vs3r1b25P04YcfYv/+/Th06FC17VesWIGZM2fCy8sLKpUKV69e1a4bP348DAwMsHr1agBAeno63N3dcfnyZXh6egIAysrK4O/vD7lcDiMjI+1rHDhwAKNGjcKNGzdgbm4OABg0aBBCQkIwb968Bnr3RERE1FKJ1nN1+vRplJeXIzg4WNsWEhKC+Pj4aj1TALB7926sWbMGb7zxRrV1cXFx6Nu3r3bZ1dUVbm5uiIuL07Z98skn6NKlCwYPHlxtX39/f22wqqojNjb2od4fERERtU6GYr1wVlYWbG1tYWxsrG1zcHBASUkJ8vLyYGdnp7P9zz//DADa3ql/HsvJyUmnzcHBARkZGQCAc+fOYdmyZTh9+nS1y4L327cmpaWlKC0t1S5rNBrk5+fDxsYGEomk9jdNRERETYYgCLhz5w6cnJxgYKC//ibRwlVxcTFMTEx02qqW7w0uD3Os0tJSCIKAV155BZGRkXBwcKjXvrWZP38+IiMj61UjERERNU3p6elwcXHR2/FEC1cymaxagKlaNjMz08uxzMzM8M0330CtVuOVV16pdd+8vLwa963NnDlzMGPGDO2ySqWCm5sb0tPTYWVlVa/aiYiISBwFBQVwdXWFpaWlXo8rWrhydnZGbm4uKioqYGhYWUZ2djZMTU2hUCjqfazs7GydtuzsbCiVSnzzzTc4fvy49oMrKyuDWq2GhYUFUlJS4OzsjOTk5Br3rY2JiUm13i4AsLKyYrgiIiJqZvQ9pEe0Ae3dunWDkZGRzqDzmJgYBAQE1Pu6Z2BgIGJiYrTL6enpSE9PR2BgINatW4fk5GQkJiYiMTERkydPRs+ePZGYmAgnJycEBgbi5MmTuHv3rk4dgYGBD/8miYiIqNURrefKzMwMEyZMwOTJkxEVFYXMzEwsWLAAUVFRACp7j+RyOUxNTe97rIiICISGhiIoKAgBAQGYPn06RowYoZ2G4V7W1tYwNTWFj48PAKBfv35wdXXFiy++iPfeew/R0dFISEjQ1kFERERUH6JOIvrFF1+gR48e6N+/P6ZOnYrIyEiMGTMGAKBUKrFp06Y6HScoKAgrVqxAZGQkgoOD0aZNmzqHI6lUiu3btyMrKws9evTAunXr8NNPP8HNze2B3xcRERG1XqJNItqSFBQUQC6XQ6VSccwVEVETp1arUV5eLnYZ1AiMjIwglUprXd9Q39+iXRYkIiJqTIIgIDs7G7dv3xa7FGpECoUCjo6OjToPJcMVERG1ClXByt7eHmZmZpz0uYUTBAHFxcXIyckBgH+dBUDfGK6IiKjFU6vV2mBlY2MjdjnUSKpuisvJyYG9vf2/XiLUJ1EHtBMRETWGqjFW9Z2kmpq/qnPemOPsGK6IiKjV4KXA1keMc85wRURERKRHDFdERETNwK1bt/Dmm2/C09MTZmZm6NixIxYvXgyNRnPffQ8dOqSXHpx58+YhNDS0TtuWlZVh5cqVD/xaq1evhoeHxwPvLyaGKyIioiYuLy8PvXr1wvHjx/Htt98iOTkZ8+bNw8cff4zp06eLXV6NNmzYgI8++kjsMkTBuwWJiIjqQa0RkJCaj5w7JbC3lKGXpzWkBg07ruftt9+GiYkJ9u7dC5lMBgDaHqzHH38cr7/+Otq1a9egNdRXa56jnD1XREREdbTnTBZCPj2AZ1fGYfrGRDy7Mg4hnx7AnjNZDfaapaWl2LhxI1577TVtsKoyYsQI/Pbbb3B3d0dKSgqGDh0KS0tLyGQy9OnTB2fPnq3xmJcuXcKwYcNgYWEBNzc3fPnllwBqvnwYHh6O8PDwGo+zatUqdOjQAcbGxrC1tcXUqVOhVqtx6NAhvPjii0hLS4NEIsHVq1chCAL++9//wsnJCQqFAiNHjsS1a9e0x7p+/ToeffRRmJubw9/fH5cvX36IT01cDFdERER1sOdMFiLWnUSWqkSnPVtVgoh1JxssYF2+fBmFhYUICAiotk4ikaB///4wMjLCyJEj4enpicTERBw9ehQVFRWYPXt2tX1KSkowZMgQWFhYID4+Hl999RXeeecd7Ny5s151HT58GNOmTcPHH3+MCxcu4Ouvv8a3336L7du3Izg4GIsXL4aLiwuysrLg6uqKpUuX4ocffsD69esRFxcHBwcHDBkyRDtFQlhYGNRqNRISEjB79mwsXrz4gT6vpoCXBYmIiO5DrREQGZ2Cmi50CQAkACKjUzDY11HvlwirHtcjl8tr3ebu3buYPHkypkyZAnNzcwCVPU6fffZZtW337duHmzdvIioqCpaWlujUqRP+97//1XuCTQsLC3z77bcYM2YMAMDDwwMLFy5EcnIyxowZA7lcDqlUCkdHRwDAZ599hmXLlmkHxK9YsQJKpRJ79uyBl5cXYmNjkZaWBjc3N3Tq1AknTpzA5s2b61VTU8FwRUREdB8JqfnVeqzuJQDIUpUgITUfQd76nQG+akb5W7du1bqNubk5IiIi8P333+P48eM4d+4cTp48CQcHh2rbnj9/Hu3atYOlpaW27cUXXwRQeVmwrnr06AFTU1PMnTsXycnJSEpKwsWLFzF06NBq2xYWFiIjIwNPP/00DAz+vmh29+5dXLhwASUlJbC2toabm5t2XUBAQLMNV7wsSEREdB85d2oPVg+yXX14e3tDLpfjxIkTNa5//PHHsXPnTgQEBGD9+vXo0KEDIiMj8fnnn9e4vZGRUa2vVdN0DRUVFTVuu3fvXvTo0QPZ2dl49NFHsXXrVjzyyCM1blt1jC1btiAxMVH7c/78eW2w++cAeGNj41rrbOoYroiIiO7D3lJ2/43qsV19GBoa4plnnsHSpUtRVlamsy46Oho7duzAlStXcP36dRw8eBCzZs3CoEGDcO3atRrv2Gvbti0uXbqE4uJibdvMmTMxbdo0baC5c+eOdt2VK1dqrGvlypV46aWXsGLFCkycOBEdO3bE5cuXta95b1BTKBSwt7dHdnY2fHx84OPjAzc3N7z11ls4f/48/Pz8cOvWLVy6dEm7z6lTpx7g02oaGK6IiIjuo5enNZRyGWobTSUBoJRXTsvQEObNm4eCggIMHToUhw8fxuXLl/Htt99iwoQJmD59OgICAlBYWIiff/4ZV69exapVq7B06VKUlpZWO9bQoUPh6OiIV199FefOncOOHTvw9ddfY+jQoejUqRNMTU3x0UcfITU1FZ9//nmtIcfGxgZHjx5FUlISkpOTER4ejqysLO1rmpub49atW7h48SIqKiowY8YM/N///R+io6Nx8eJFTJo0CUeOHEGHDh3QsWNHDBw4EC+99BL+/PNPbN++Hf/73/8a5LNsDAxXRERE9yE1kGDuSF8AqBawqpbnjvRtsPmuHB0dceTIEXh5eeH555+Hn58fFi1ahA8++AALFy5EUFAQ3n//fUyZMgVdunTB6tWr8dVXXyEnJweZmZk6xzI0NMT27dtx/fp1dO/eHdOnT8eCBQswfPhwWFlZYeXKldiwYQM6deqE06dP47XXXquxpnnz5sHe3h6BgYEYPHgwZDIZIiIitGFswIAB8PHxQefOnZGYmIiZM2di0qRJeOWVV9CtWzekpaVh7969aNOmDQBg06ZNsLW1RVBQEObMmdNkJ0etC4nQmmf50pOCggLI5XKoVCpYWVmJXQ4REf1DSUkJUlNT4enpWW2uqPrYcyYLkdEpOoPblXIZ5o70xTA/pT5KJT37t3PfUN/fvFuQiIiojob5KTHY17HRZ2in5oXhioiIqB6kBhK9T7dALQvHXBERERHpEcMVERERkR4xXBERERHpEcMVERERkR4xXBERERHpEcMVERERkR4xXBERERHpEcMVERFRExYaGop58+bdd7uysjKsXLmy4QuqhUQiwaFDh+q0bWJiIo4ePQoAOHTokM5DnlsChisiIqIWYMOGDfjoo4/ELqNORo8ejQsXLgAAgoODkZWVJXJF+sVwRURE1AI0p0cF31ursbExHB0dRaxG/xiuiIiImoHVq1cjNDQUc+fOha2tLRQKBWbMmAFBEHDo0CG8+OKLSEtLg0QiwdWrVyEIAv773//CyckJCoUCI0eOxLVr17THk0gkeP/992Fra4tRo0Zh9erVCAkJwdtvvw1LS0u4ublh1apV2u01Gg0+//xzeHl5wdTUFP3790dSUlKNtWZmZiIsLAxt2rSBiYkJ/P39ceTIEQCVlznT0tLw4osvIjw8vNplwYyMDIwdOxbW1tawtbXFtGnTUFpaet/PoClhuCIiolatuKyi1p+ScrXet30YR48exfnz53HkyBEsXboUS5Yswf79+xEcHIzFixfDxcULine truncated
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# As a guide to the eye, a Cubic Spline will be added\n",
"\n",
"interpolation = CubicSpline(np.asarray(beta),RMSE)\n",
"x = np.linspace(np.min(beta),np.max(beta),201)\n",
"y = interpolation(x)\n",
"\n",
"fig,ax = plt.subplots()\n",
"\n",
"ax.scatter(beta,RMSE,label='Calculated')\n",
"ax.plot(x,y,'--',label='Interpolation')\n",
"\n",
"ax.set_xlabel('Beta')\n",
"ax.set_ylabel('RMSE')\n",
"\n",
"\n",
"x2 = np.linspace(0,100,501)\n",
"y2 = interpolation(x2)\n",
"\n",
"xmin,ymin = x2[np.argmin(y2)],y2[np.argmin(y2)]\n",
"\n",
"fig,ax = plt.subplots()\n",
"ax.scatter(beta,RMSE,label='Calculated')\n",
"ax.plot(x2,y2,'--',label='Interpolation')\n",
"ax.scatter(xmin,ymin,label='minimum',color='r')\n",
"\n",
"ax.set_xlabel('Beta')\n",
"ax.set_ylabel('RMSE')\n",
"ax.set_xlim(0,100)\n",
"ax.set_ylim(0.1010,0.1040)\n",
"ax.legend()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "de8fef6b-5e63-43b7-bbf7-928d2993adda",
"metadata": {},
"outputs": [],
"source": [
"# If the above has been run, utilize that value otherwise set ymin to 23.4\n",
"\n",
"xmin = 23.4\n",
"\n",
"func = lambda beta: AB.cross_validation(q=0.75, beta = beta,verbose=False,n_epochs=20)\n",
"\n",
"result = minimize(func,x0 = [xmin])\n"
]
},
{
"cell_type": "markdown",
"id": "4425261a-3419-499d-ab84-1ac281750b8c",
"metadata": {},
"source": [
"## Generation of AMBER background"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "804af4a9-47e5-49cf-bb7b-7d5aa51022d8",
"metadata": {},
"outputs": [],
"source": [
"AB.beta = xmin"
]
},
{
"cell_type": "code",
"execution_count": 40,
"id": "9bd8bd87-6230-4d1d-a980-3c83c5fc1573",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Iteration 1\n",
" Loss function: 89216.17289327279\n",
" Iteration 2\n",
" Loss function: 68093.90388265462\n",
" Iteration 3\n",
" Loss function: 57929.508628268115\n",
" Iteration 4\n",
" Loss function: 51147.228477046934\n",
" Iteration 5\n",
" Loss function: 46560.03500671341\n",
" Iteration 6\n",
" Loss function: 43283.23323431091\n",
" Iteration 7\n",
" Loss function: 40937.53702364261\n",
" Iteration 8\n",
" Loss function: 39308.08183709762\n",
" Iteration 9\n",
" Loss function: 38236.35126931614\n",
" Iteration 10\n",
" Loss function: 37568.299703508696\n",
" Iteration 11\n",
" Loss function: 37178.28929735461\n",
" Iteration 12\n",
" Loss function: 36963.03718997033\n",
" Iteration 13\n",
" Loss function: 36850.84062922201\n",
" Iteration 14\n",
" Loss function: 36796.66347072205\n",
" Iteration 15\n",
" Loss function: 36772.34667854443\n",
" Iteration 16\n",
" Loss function: 36761.725670802814\n",
" Iteration 17\n",
" Loss function: 36757.366819518116\n",
" Iteration 18\n",
" Loss function: 36755.782935800766\n",
" Iteration 19\n",
" Loss function: 36755.3751834626\n",
" Iteration 20\n",
" Loss function: 36755.43029325381\n",
" Iteration 21\n",
" Loss function: 36755.63232500783\n",
" Iteration 22\n",
" Loss function: 36755.84605170578\n",
" Iteration 23\n",
" Loss function: 36756.02667786971\n"
]
}
],
"source": [
"AB.generateAMBER(n_epochs=25)\n",
"AB.applyAMBER()"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "7f104c8d-50a2-4ae7-80ef-75f31a0083c2",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAHVCAYAAAD/8I8fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAACTlklEQVR4nO39eZhV1ZX/j69z56HmebCYHRBEHEACOLVJcIgxKm1sDcFWEzUqGhMVOx/jmJiO2pq0UX+2CSbt10iMGjvGmHTHEQWMCKhMDkwFVAE1V91bdzz79wdSsajzXtyqoqpuFe+XTz2PnHX3Ofusvc+tXft91lqWMcYIIYQQQshBjmuoO0AIIYQQkg1wUUQIIYQQIlwUEUIIIYSICBdFhBBCCCEiwkURIYQQQoiIcFFECCGEECIiXBQRQgghhIgIF0WEEEIIISLCRREhhBBCiIhwUUQIIYSQLOaTTz6ROXPmSE5OjowaNUruvfde+NmVK1fKCSecIKFQSKZNmyYrVqzo1bW4KCKEEEJIVmLbtpx11llSWloqK1eulEcffVTuvvtueeqpp3p8NhKJyJlnniknnniirFixQmbOnClnnXWWRCKRjK/HRREhhBBCspKdO3fK1KlT5ZFHHpFDDz1UzjzzTDnttNNkyZIlPT67ePFiCQaDcu+998rEiRPlwQcflNzcXHnmmWcyvh4XRYQQQgjJSiorK2Xx4sWSm5srxhh566235I033pBTTjmlx2eXLVsms2fPFsuyRETEsiyZNWuWLF26NOPreQ5Uxwcb27Zlx44dkpub2+UAQggh5GDAGCPt7e1SVVUlLtfg7G/EYjFJJBL9Po8xpsfvbb/fL36/X203ZswY2bp1q3zlK1+R888/v4e9rq5OJk2a1O1YeXm5fPjhhxn3bdguinbs2CE1NTVD3Q1CCCFkyKitrZVDDjlkwK8Ti8Vk7Ogcqd+V7ve5cnJypKOjo9ux2267TW6//Xa13bPPPiv19fVy1VVXyXe/+135+c9/3s0ejUZ7LKz8fr/E4/GM+zZsF0W5ubkiInLEpT8Uty/Qw267cdtkDrYZpV2wAdtcKefjpcuaYBsrZeMTphWbRgp0RERMuKefMiEdxqt3V0p5QAw2xUuCjsd9zTHYxkrjE7aPwYMarsPnjFZgnyTDyg6kZlKGTpuXtheftLPU+bi3w/m4iEi8ANvSQezLf/3yq9C2oHAjPmkfcZWvPODnJGSk09bWJjU1NV2/CweaRCIh9bvSsmnFaMnL7fvOVFu7LWOP2yK1tbWSl5fXdXx/u0QiIscff7yI7FmgXXzxxXLfffeJz+frsgcCgR4LoHg8LqFQKOP+DdtF0d6tN7cvIG5/z19slvbLR/G9UTzi9mEb2r30uPHFLKMtfPq4KFJ+6xqlLxqWBy8cXEZbFOFfvGlwTo/if0vxiceL+6idU2tn+w78okibl5ayKEJD507i87mVNbAJ4LEJ5GCH5eUqN9BHXJ/7YiSE9I7Bfn0knLPnp6/s/ds2Ly+v26IIsXPnTlm6dKl87Wtf6zp25JFHSiKRkLa2NikpKek6Xl1dLfX19d3a19fXS2VlZcb9G7aLor24EyJOX9MevGEiaW3DRPkdn1bWFKGdfVjEtOMwQROJQpvlV1ZnygrAau/bQsvT3I6N2gMZdt4NEhEJRsB2ZksbbGNXg+0SEclftQv3w43/qsnpxPp4OkdZ0CaxL20fvl4q1wttRvnjK7zT2c/xfNzIh10pkSo8br/63Zeh7eFy/IAcM2UTtB2dvw3abq+AJkLIQc6mTZvkvPPOk9raWqmurhYRkRUrVkhpaWm3BZGIyIwZM+QnP/lJ1ztLe1/M/sEPfpDx9Rh9RgghhJCMsMX0+6c3TJs2TY477ji59NJLZe3atfLSSy/JjTfe2LXQqa+vl87OThERmTt3rrS0tMj1118va9euleuvv14ikYhccMEFGV+PiyJCCCGEZIR9AP7rDW63W1544QUJh8PyhS98QS6//HJZsGCBLFiwQET+EbIvskeSe/HFF+XNN9+U4447TpYtWyYvvfSShMPhjK837OUzKyViOSztosqWfKIQr1StZO/f6xARaTrSeX3pTuTDNq40tvmasayjySy2H7/zEStSpBvlVZHwNvyicjqAp5D24nAy7HzBaHkZbNN2KH6YCtZhXzZ9AUce5K3EgxqZgSXMdCf2ZX4Jfvu5XVEiv3zYOmi7oPgdx+MvtU6Fbf4pbw20nR7sf1gtIYQMBlVVVfLcc8852sw+765Onz5d3nvvvT5fa9gvigghhBAyOKSNkbQSRJNJ+2yGiyJCCCGEZERf3gvat302M+wXRcGGlHi8PUPNPHGsB0ViiqzTx+jgIAh+sj34WsGdndDmbsHSja3kDfLWtUKbpxXHURqfIruV4Gg3fzOOCfe240ildMA5Mq3kfeyT3FosWSXysJ/z38X+CrTgB9T9Bs5t4YnidrGiImjzKfPrlU+Phba/5k11PK7JvbUzC6Hte/U4RPWv0x+Ftko31uZdfXxF0a4/DJ+z4qM+nZMQMjDYYiQ9ghdFfNGaEEIIIURGwE4RIYQQQgYHymdZjjthxG33dHKgWcu0jKWiaBrLEVrSRxQR5k4okW4JLcMkjrRy7WrB7ZTEju4GJZufH0tT4Va8oaiV3kgVYvkpZ32z8/lsfN/uKB6A8EYs46VzcTtLK6eiJKa0vdgnuVvx/NIi8izlu6JtlPO4emK40fqth0ObC+fVlJO33Ahtdhg/V9VjcB2cP016Ctr+Hsey7peghRAyFIz0F60pnxFCCCGEyAjYKSKEEELI4GBLnytzdrXPZob9oigVdIk4SBneDrzN7+nUqnli6SOlJMVMFDgf90bx+SzlhAFl5pgCJQIorshISi0vdwQn89NkvnQ+1mE6y7FsFXKQPEVELSJrPMrYFOJrueN4LjRNwr4M7cbtdh6H+xLEKpLEirHNTMRJH5N1zo+qVYITU06sroe2owp2QNvdZaugra8RZiJ4npwWVKRuQkhWke5n9Fl/2g4GlM8IIYQQQmQE7BQRQgghZHBImz0//WmfzQz/RZH57Gcf0n68CZbMwTYvzpko0So8mukiZ4mpMYBlFrNaqRvmz4W2VBDLfy4l2s3TiTU5u1SpAVaJ++lv1XQ+bKqf4RxxFK7H54vn43FLKEkRfTifpcRxfkOJVOKxO+X0VdD213eOgrajjtoCbdvbcf22s05xruUzO2cDbPOloBLdqMINZEKIM3yniBBCCCFERGyxJC3Ke7kZtM9m+CchIYQQQoiMgJ2i3I+axePuKf3YQZyMUESpAaYsE10JpZ5atfP1vDigSEK7sbwR2oITLVrtuD6YloRRJYUjgELbcYSWURIcWml8zvB252SEthf72BvBNvcmvClrXLiPqQZFksvFtnU/whJZRRi327FqHLQFFCnyxeKTHI8/Pcr5uIiIV8nVmcTqrKRqYtBWVIwndFEQa89/OPw5aAtaWLplXTRCsgvb7PnpT/tsZtgvigghhBAyOKT7KZ/1p+1gQPmMEEIIIURGwk5R2vldeFc7lgCCO5T6ZmFcOywVwmtIf7PzObWoKOPG/TAupd6YVxm2Tnzf6r6lB5/T1axogEo/TQBLed7dzudMFWGpzlffDm1WQklaWYK1In+9ImEqUp5GQKmLZiljkA4qUX4tzuf0RrCPO4vx/HLjXJ2S6sSJMDttbKs4ew20nbH2Amh7bdILuDOEkKxipO8UDf9FESGEEEIGBdtYYpt+RJ/1o+1gQPmMEEIIIURGwk6RnRaRnlFOJohlME0iUxexivrUMtE50srfiCWYziK8JvVEQtDmr1NqkRXj4lruKNZMjEdZHyv1yCwlPal2zlSBcy2sVAhPyVgZlm78jbgGWPtoXHdLwxvF0WDNhyqRcIo05YlgW8tEzc/Ox/2NSmSdVquvGneytBxnuywI4sjHsAePASUyQkYGlM8IIYQQQkQkLS5J90Nkyvbyz1wUEUIIISQjTD/fKTJZ/k7RsF8UpQtzxHL3lFWSBTghnCeCI5XiikSj1d4KbXc+nlKUm1gRnhyh3fhadg7uo1EinyJjcBSWS1u+KxKZO44lprYxWKZE17OUwji2MltthwSeezFKEFkyhMfA1qLIFH9ZSsmxeAG2uTtxXxKlzidNl2KHHTcB11l7Ztz/4Y4ME5jYkRByoBn2iyJCCCGEDA58p4gQQgghRETSxiVprR7WftsfwM4MAMN+UeSuaxK3q6dM427A0o0ksXwWihVAm7cDa2HRduckei2H4smjSVbedqzBuDpwgkbLh4c0ZxdOfqjVPjMhLE1ZSdyudKOSSDLoLAFqNetsJaJQe0ZzNiuJHUP4elokXN5GfN/JPHzOaBnW8vK24m+LjkrnvviUvJqblh4Kbcf4sS2h1EVLFCjJJ0uwn0+f/CG0PVy9HF+QEEIGkWG/KCKEEELI4GCLJXY/os9sLbdNFsBFESGEEEIygu8UZTsuy7H+lmlrg00sS5G0lASHHjeWPnxBZ5u7U6uDBU2SDijXUkRZq0mRyJT6YOLDko/VhJPyqQSw7IbkOlcz9r+Vwg7TEkVaMSxFeuLYJ5YSOqglwnTHsbTmb9Dq7uHHMdgA2ijzJO3D1+osxu20yDpfKz5nh+A51JzAmSRfi+FznhLI7r8qCSEji+G/KCKEEELIoND/F62z+w8dLooIIYQQkhF73inqR0FYymcDjMcj4up5G1ZRIW7jVla5ShSWK4FlGLRw1upgaVE+kQo8NN5WXBdNG9DYoeXQ5mtWIsUUtCSZ3iZcJyuV69zOUv6I0OSztvHYJ94O3M7bgcc7VoK96UrgBJqJPKXmW0BJ0JgHTZIocD7uVaLPOo/AY+rZhuWzUSfUQtum+hJoe/wLv4Y27S9LSmSEDB/sfpb5yPYXrft+Z4QQQgghI4jhv1NECCGEkEGB7xRlOz6viNsh6iWpFKCyFE3TqxXYwjJMYLdzhJYLine truncated
"text/plain": [
"<Figure size 640x480 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAHVCAYAAAD/8I8fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB45UlEQVR4nO3deXxU9b0//teZmex7WJIQwu6CIIIKUsDtci0qVatSa7UUr9qqV0Vrq+Ltw+tGa2/Viv1a9eHFYnv9IRSXel2rj+uKLGpYVDaVRQIkQAhJyJ6Z8/n9gURD5v1m5jNLJpPXs488HnU+c875zOecCZ+c93m/P44xxoCIiIiol/N0dweIiIiIEgEnRURERETgpIiIiIgIACdFRERERAA4KSIiIiICwEkREREREQBOioiIiIgAcFJEREREBICTIiIiIiIAnBQRERFRAvvqq68wbdo0ZGdnY9CgQXjggQfE965evRqnnHIKMjMzMX78eJSXl4d1LE6KiIiIKCG5rovp06ejX79+WL16NZ544gnMnTsXCxcu7PLexsZGnHvuuTj11FNRXl6OSZMmYfr06WhsbAz5eJwUERERUULavXs3xo4di8cffxxHHXUUzj33XEydOhVLly7t8t7FixcjIyMDDzzwAEaOHIl58+YhJycHS5YsCfl4nBQRERFRQiopKcHixYuRk5MDYww+/PBDvP/++zjjjDO6vHfFihWYMmUKHMcBADiOg8mTJ2P58uUhH88XrY7Hm+u62LVrF3JycjoGgIiIqDcwxuDAgQMYMGAAPJ743N9oaWlBW1tbxPsxxnT5dzstLQ1paWnqdkOGDMH27dvxgx/8ABdffHGX9srKSowaNarTa0VFRfj8889D7luPnRTt2rULZWVl3d0NIiKiblNRUYGBAwfG/DgtLS0YOjgbVXsCEe8rOzsbDQ0NnV676667cPfdd6vbPf/886iqqsJ1112HX/7yl/jTn/7Uqb2pqanLxCotLQ2tra0h963HTopycnIAAMde+Z/wpqZ3aXe98rbt2XKbUbbLqJbbPP7gr/dbUSNu4/hdeYcBpU3jFzoCwGR1HadQBLLk2bvHr3xBjNzU2jcj6Oup+1vEbZyAvMMDQ+STmlUp77OpWB6T9izlDqTWpJw67bp0U+Sd+oVu9tkgf9kbSlPFNk+73I/zf/WO2Da7YIu8oSVP0eqo75Mo2dXX16OsrKzj38JYa2trQ9WeALaWD0Zujv2dqfoDLoae9DUqKiqQm5vb8fqR7hIBwMknnwzg4ATt8ssvx4MPPojU1G9/z6Wnp3eZALW2tiIzMzPk/vXYSdGhW2/e1HR407r+i+Fo//goY2+UEfHK/8ZAunvp88oHc4w28bGcFCn/6hqlLxrHJ08cPEabFMmTmICwT58y/o4yJr4UuY/aPrXt3NToT4q069JRJkVGOHU+n7yNN1W+YJVuID1bHrDcHG1LO57v/GIkovDE+/GRrOyDP7YO/W2bm5vbaVIk2b17N5YvX44f/vCHHa8dd9xxaGtrQ319Pfr27dvxemlpKaqqqjptX1VVhZKSkpD712MnRYf0W9UIn6/rP8y+/U3iNq0l8okwXuUfmRb5LkzKl7uCvu7fvUfcJhlYTt2QIryu3FxS27JW2fVD+/sh9L8tElO+5XZvLZR/472FEyz3qhzP9iIioqS3detWXHTRRaioqEBpaSkAoLy8HP369es0IQKAiRMn4ve//33HM0uHHsz+zW9+E/LxmH1GREREIXFhIv4Jx/jx43HSSSfhyiuvxPr16/Haa6/h1ltv7ZjoVFVVobm5GQAwY8YM1NbW4uabb8b69etx8803o7GxEZdccknIx+OkiIiIiELiRuF/4fB6vXjppZeQlZWF733ve7j66qsxe/ZszJ49G8C3KfvAwZDcK6+8gg8++AAnnXQSVqxYgddeew1ZWVkhH6/Hh898+5vg83YNn7UVyQ+ftefKH1t7mNdRnsJuGxn86X9v9T55f8ozH+43M99o8g1WsvW05392VoptjvJwnJOutEnb5cnnra1YbkvdFDx8CQC1pw8R2/JekVM12045RmzTeBvlMKtvv1xZ1V+gfHE9Qlh32ZoQe9WZNz9fbGs4Q/7cTf3l70DNifLn3nL+f4fULyKiww0YMAAvvPBC0DZz2L9dEyZMwKpVls9TIAkmRURERBQfAWMQUP6IDmX7RMZJEREREYXE5rmgw7dPZD1/UlRbD3i6hqFSG+TsM+/AvmKbP1PKiwIcpXZQamV90NcDjvzYVixCZBr/1xVW2zleOWTiNsnjDK1NUlklNnk2ypvJgRsge5G8T60Mmff/wltdORRaP+MpUFsrtmX8Y6Xcpuyzj9I2/d6zxbZXPzlabPMUf6HslYjizYVBIIknRXzQmoiIiAjJcKeIiIiI4oLhs0RnTNDMKaOEpryV8tIbKCkUmwLp8nC5mcGzqbwFeeI2/r3KuiEJxAQiX+uGejdXycKcNkAuCMnCjkSJJdkftGb4jIiIiAjJcKeIiIiI4sKF/fJOh7ZPZD1/UpSVCXi6hq7Mnr3iJk6KnGGmCaTLWVjtZcHXi8qqPWB1rGTgCbJQ7yFuq7xyfTx5MuR8qnhnB/oGyIsWBoQ19JIhtKllNxJRYglEmH0WybbxwPAZEREREZLhThERERHFRcAc/Ilk+0TW4ydFpq4exulavFErKqityeXdIWeEObkDxLY2aT21MXJIJPNAg9gWqNkvtsU75OPR1mhra7Nqk9be0ooKaqz7aDle7mnj5L68v1ps8w0ZLLb5t30t71P4fIkUPnNPP1FufH+N2PTlU8p2RJRQ+EwREREREQAXDgIQFqgOcftExmeKiIiIiJAEd4oCtXVwnK7ZZFpGi5MuZ0WhvV1sSt9eK7a53oKgr6dVK1lWbfKxNPHOitLCTyoj3yi1CZM5Pjlr0LqPlrQQmUYLkWni/flseN5bZbXdUVd+Ira5Z3NdNKJE4pqDP5Fsn8h6/KSIiIiI4iMQYfgskm3jgeEzIiIiIiTBnSLH64XjdA2VaVk5gT1yhpknQw6tObk5YltqbfDwRkuxnCmWtTN4wUcAQIOcmZbUHHmebvx24UaKDy2zLlCxQ2z7Yt54sW3YC6eIbdv+PbR+EVH0JPudoh4/KSIiIqL4cI0D10SQfRbBtvHA8BkRERERkuBOkQkEYJSQSzBOivyxAwfktcp8Smht/6lFQV/P3KsU18vOEps8mZlim1aYssdTMtYSibZOmVtbJ7cp587Xr6/Y5q+uCd6QQONlm1l31I0rxLbAmSfJGzJ8RhR3DJ8RERERAQjAg0AEQabEqcEfHCdFREREFBIT4TNFJsGfKerxkyIp+0xaWwsA/Mq6Yt7C4EUYAcD0k9vyNgcv0uimykUk3Vw5M82pjP7jXt68PLHNtLSKbWrhQCV84+vTR97M7w/+unYsR/4yxTuk6N9VGf197pWzIrWsPIkWgm2eOlpsy/pU+WxKUVSjhA2dfPna+/Ln8pqChSfsFdvcKhZ2JKLo6vGTIiIiIooPPlNEREREBCBgPAiYCJ4p4jIfsSVln/n37bPan1tXL7Z5ldBBalvwx8cajpMzijztqWJbihIqshWok8MbsWB7DpKZtiafVnDUJstMCymmvfyR2BY8sBkhJWQ9/Dm5iKlnnnLN2i0/R0Qk6vGTIiIiIooPFw7cCLLPXCT2rSJOioiIiCgkfKYowTm+FDhOSpfXbdfJMq48i9VCa560tKCvO0oA1fjki8MpkLN1oBSYpMSmhsiSmBY29DTKmY/1U4aJbdMu+pnY9tay0PpFRPRdPX5SRERERPER+YPWDJ8RERFREjj4TFEEC8IyfBZbxt+OoAUytWJ3WiaP0qYVMfR4gp9ob6u8v7Ycefh9xXKhSM9eOavLbW6W96kUUwwooUEt5OPNy5X3WVsrtonnx3ItL9/gMrHN1Mj90Na604pdBuqV7XJz5L60y2FdT4ZczBN9g18P/k1fyfs74TixzV27XmxzTj5ebDOffCa3TRkrtvm27BbbNt8nf+72XWITHFcuTklEseFGuMxHoj9oHf2yyUREREQ9UI+/U0RERETxwWeKeirLMIwtf0Xw+/ypBXJ4yR0oF61rzwuezQYAPiVEpolFMUU1RKbwZmcF358SztIymPxfV1j1Q6MWu1TCs7ZFMtX12yzOnRYi02ghMo2zdI3YphWEHH6HHD7zf/Wp2KaFg4cX/S7o65sv/Q+lJ0R0JC48SV2niOEzIiIiIiTznSIiIiKKqoBxEAia3RT69oms50+KHE/QUIZXKX4YUNZh0ngy5WwXJ6VrAUkAgF8OHARS5IvDIye6qaEbT4p8SrUssngXFZSOp2Z8xXntNo2vb6HY5t9bLW9omxXZA/j695MbC/PFJn9h8FAqAHgzjpW3+2yj2PblGW8LLQyfEUUiEGH2WYDhMyIiIqLE1/PvFBEREVFcuMYDN4LsM5fZZzFmXABdww62ITKNmh0k8PWRizD6mpXCjrnyqUnNSBfbbPrYHXpKPyVqiEzTw0NkGv+evWKbOWag2Jayt0He58YvrfoyfeL0oK+/vs1qd0T0jWQPn/X8SRERERHFhYvIHpZO9D8L+UwREREREZLgTpE3OwteJ7XL64EG+Za8RisQqGVoOT4h++yA3A+PXy4+16qEz3KGyqEIrPtCbqOwiOcUB9fck3hz5LXPAg2N8gF7eGhN+9wH+svFSFO27ol6X/zbdwR9/SzPj8Rt3nKXRL0fRMkm8uKNiX0vpsdPioiIiCg+Il/mI7EnRYndOyIiIqI46dY7RS+++CIuuuiiTq9dfPHFeO6550Leh9vcAtfpGtbypHYNqXVs0yZXRrQtYmiznfHID6ulHpCf0A/kymtFaWwLI2phJG+evLZbtNdaU4tnKufbdn02T7oc8gk0yOEzbf02bb2uaI+XN1teW882vKxdQ2ZEmdjWWCz//VX7k8FiW/4WeZ/eVjncuP3i4G1Line truncated
"text/plain": [
"<Figure size 640x480 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAF2CAYAAABtfCOsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB0vElEQVR4nO2de3xU1bn3f3smk7kkk5ncSELAgCgthbacAFYDKogevNSqES9wbIPKq1BfkapBensBb7UVEVu8HMsx0GNFqEI9tUBKD0VBCRUUjxKkhxIghIRLQu6ZyWRmv39Msp9nkxkyM7lNMs+Xz3xYs/faa6+1Z2Zl7d9+LoqqqioEQRCEQY+hvzsgCIIg9A0y4QuCIMQIMuELgiDECDLhC4IgxAgy4QuCIMQIMuELgiDECDLhC4IgxAgy4QuCIMQIMuELgiDECDLhC4Ig9AGHDx/GjBkzkJiYiIsuugjPP/980LqfffYZvvOd78Bms2HSpEnYt29fj/RBJnxBEIRexufz4aabbkJ6ejo+++wzvPbaa3j66afx1ltvdarb1NSEG2+8EVdeeSX27duHvLw83HTTTWhqaup2P2TCFwRB6GVOnTqF8ePH49VXX8Wll16KG2+8EdOnT8euXbs61V2/fj2sViuef/55jBkzBitXroTdbscf/vCHbvdDJnxBEIReJisrC+vXr4fdboeqqvjoo4/w4YcfYurUqZ3qlpSUYMqUKVAUBQCgKAomT56M3bt3d7sfcd1uYRDj8/lw8uRJ2O127eILQqyhqioaGhowdOhQGAyDa43ocrnQ2toa8fGqqnaaG8xmM8xmc9BjRowYgePHj+O73/0ubr/99k77KysrMXbsWN22jIwMfPnllxH3swOZ8C/AyZMnMXz48P7uxoBjhn2OVi5uWNNv/RB6lvLycgwbNqy/u9FjuFwujMxJRNVpb8RtJCYmorGxUbdtyZIlWLp0adBj3n33XVRVVWH+/Pn40Y9+hF//+te6/c3NzZ3+YJjNZrjd7oj72YFM+BfAbrcD8H/Rk5KSev18NTU1GDlyJMrKypCSktLr5+st7hz2iFauq6sL6ZjBMvZIiPax19fXY/jw4drvYbDQ2tqKqtNelO3LQZI9/DuX+gYfRk441ml+uNDqHgAmTpwIwP8H59/+7d+wfPlyxMfHa/stFkunyd3tdsNms4Xdx/ORCf8CdNyqJSUl9cmE7/F4APj/0PTF+fqCUMcxGMceKgNl7INV1kyyGyKa8LXjQ5gfTp06hd27d+PWW2/Vtn3jG99Aa2sr6uvrkZaWpm3Pzs5GVVWV7viqqipkZWVF3McOZMIXeoQbMubTm3iTVpzwwIsAgH2v/6ivuyQIIeFVffBGkPfPq/pCrltWVob8/HyUl5cjOzsbALBv3z6kp6frJnsAuPzyy/Hcc89pzwc6HvL+9Kc/Db+T5zG4nsAIgiCEiQ9qxK9QmTRpEiZMmID77rsPpaWl2Lx5MwoLC7VJvKqqCi0tLQCAmTNnora2FgsXLkRpaSkWLlyIpqYm3Hnnnd0eq0z4UYTZbMZdd93VpQY4GJGxx+bYowFfN/6FitFoxHvvvYeEhARcccUVmDt3LhYsWIAFCxYAILNNwC8Rvf/++9i5cycmTJiAkpISbN68GQkJCd0eq0g6UYTZbMasWbMGzA+fyziqix4yKan04NFyLrQfxUAbe08Sy2OPBryqCq8avqYT7jFDhw7Fxo0bA+5Tz2vrsssuw6effhp2n7pCJnxBEGKacOUZftxAQyQdQRCEGEFW+EJYXHv1s1rZlER22YrToZW9KaQ12v9R2yf9EoRI8UGFN0ZW+DLhC4IQ08SSpCMTviAIMU1fPbSNBmTCF7rka0++qJWdORatbHaka+W4FopHEn+yQSsr1ecAAJf8ito4vEicsITowdf+iuS4gYZM+IIgxDTeCDX8SI7pb8RKRxAEIUaQFb6gccPFj2vl1pxU2nENyTinLqcb2aTD9PVJ39+mlRU3xRdXM/1xQi595Ti1t6hHuisIPYJXRYSxdHq+L72NTPiCIMQ0ouELgiDECD4o8CL80M++CI7pb2TCj3FmWP5NKxvSSMZRfHS/mvWxRyufGU+hjz2J1E7dxRQHRh01VCunfVrrLyR0P3mDIPQGPtX/iuS4gYZM+IIgxDTeCFf4kRzT34iVjiAIQowgK/wYhMs4itVKO5JIozGV12jluESy0kmxOLVy3QijVq4dTaudlAN0r9s83N9mfB1Z8VzyS+aE9YQ4YQn9Syyt8GXCFwQhpvGpCnxqBA9tIzimvxkQko7L5cL9998Pp9OJrKwsvPDCCwHrrVmzBoqidHoZDDTMTZs2ddo/c+bMvhqKIAhRRscKP5LXQGNArPALCwuxd+9ebN++HceOHUNBQQFycnI6TdR33XUXrr/+eu29x+PBNddcg+9+97vattLSUtx88814/fXXtW0WiwWxhMIyKympybSjxRX4APYH09RI0oyjjKSb+HqSdzyJ9ENoHOb/ihld9FWznA2/z4LQW3hhgDeCta+36ypRR9RP+E1NTVi9ejW2bNmC3Nxc5Obm4sCBA1i1alWnCd9qtcLKNOlf/OIXUFUVzz33nLbt4MGDGDduHDIzM/tsDIIgRC9qhJKOKpJOz/P555/D4/EgLy9P2zZlyhTs2bMHPl9wX7eamhr88pe/xHPPPafLFVpaWorRo0f3ap8FQRCikahf4VdWViItLQ3x8fHatoyMDLhcLlRXVyM9PT3gca+++iqGDh2quwtQVRWHDh1CcXExnn32WXi9Xtxxxx148sknde0PRnQOVkPYNePeI+watFxMicitR8hiJ/4MVTcfbtTKcZdmaGWVLSOahvodtVxOWg0ZyY8LX9v4pFb+lfHbXYxCEHoesdLpQ1paWlBRURFwX1ZWFpqbm3UrdADae7fbHfA4VVWxevVqLFqkj9J1/Phxrb0NGzagrKwMCxYsQEtLC1566aWgfaypqYHHQ7OU2Wzu1KeeoOMc/FyxRiyOPdo+d7fbrfttNTQ0XKD2wMerGuBVI9DwxdM2fPbs2YNp06YF3Ldp0yZYLJZOE3vHe5stsLv+3r17ceLECdx999267Tk5OaiurkZycjIURcH48ePh8/lwzz33YMWKFTAajQHbGzlypO79XXfdhVmzZoU0vkjYtm1br7Ud7cjY+59169Zh/fr1/d2NPsMHBb4I1G1JcRgBU6dOhXqBVGEff/wxzp49i7a2NsTF+btbVVUFq9UKp9MZ8JitW7fiqquuQnJycqd9KSkpuvdjxoyBy+VCTU1NUHmorKwMdjsl7O7NFf62bdtw3XXXwWQydX1AF3zXPkcrG5KdtMPAbkUtJOOoJvqDZ/3sKNXhclc1yTsYkkbN/KOK2kmiJOaWKv8P6cxl9FmY6+jZi/M/yELq5v8o1Mre1AOIFXr6c+8u06dPx8svv6y9b2ho6LToGUyIpBNFjB8/HiaTCSUlJZgyZQoAYNeuXZg0aZLOvp6zZ88eTJ48udP24uJizJ49G+Xl5drdwf79+5Gamhp0sgf8fySSkpJ6YDShYTKZouKH35/E4vij5XM3mUxITEzUvR/MRC7pyAq/x7HZbCgoKMC8efNQVFSEiooKLF++HEVFRQD8+rrRaITD4dCO+fLLL3HPPfd0aisvLw9WqxVz587FkiVLcOTIERQWFnbS+gcyurAJ7C5ETacVtuImW3p4WOKSM7R6911ED2ENp85RHfbjV09U0nYbmcMq7I7ANcJ/RzVk52k6p0IrIzWOfmjf+ck8rfzJGxAEoYeJ+gkfAFasWIH58+dj2rRpcDgcWLZsGfLz8wEA+fn5GDFiBNasWaPVP3XqVEA5x263o7i4GAsXLsTEiRNht9vx4IMPorCwsFNdQRBiA7+GL/HwowabzYa1a9di7dq1nfbt2LGj07aWlpagbY0dOzZqHo4JgtD/+CL0tJWHtkK/oJNxEumBqe/S4VrZeJxJKtwaiT3Abb58lFa2lfwzYB2VSUDqqIuoThuThlrJ6dzyWVn7OekH5cui5yWGRvrjnPYR9fE7P1ihlff87lEIQm8hGr4gCEKM4INBzDIFQRBiAa+qwBtBXJxIjulvZMIfBHCrGAyjoHDGMrKiATetM1P9tiFkbmr9oFQr87WLYiVbeSWd5b09cYoqJTInOC+LceRsb5/d/irlQfrF5Cj70eDPYQRBiAyZ8AVBiGkiD48sko4gCMKAwqca4Ivgoa1PHtoKfUXQ6JenqrWiOoTCSPhsJOMYT9dr5bh/nKD68SSvKDwpDAvqpbLQCmpOFvWhkcU7amZyTEN7RE3mbMVlJx5WQ3HReUxnKBLnFbMpw9nutx6DIPQkssIXBEGIEXyI7AFs8Gwc0UvUJ0ARBEEQegZZ4Q9QDCx6p+qkQFdKSyuVPeQApbPY4Xls3UzGaeMxdlhs9iQ6l5dZ9cQdJSsdlYewHjqE2mz2b1eZZZBSy+KrG5gTGIvcyS19nH+nSJyXPE8OWYcLxSFL6D6R2+EPvPWyTPiCIMQ0kXvaDrwJf+D1WBAEoQfpCJ4WySscKioqMHPmTKSkpCA7OxuPPvooXC5XwLq33HILFEXRvd5///1uj1VW+AOM65PnAgAUJ0krSiP70uhkGfpCto3K1srGA0eoDsspoFhYUhcrhTuGm2SiuH+epO2sPo/hgzMUTlnbz0MiswQpujr1TfSGO2/Zqf7I95gFkAQ5FXqAvljhq6qKmTNnIjk5GTt37kRNTQ3uu+8+GI1GPP/8853ql5aW4s0338T06dO1bYEiAIeLTPiCIMQ0kZtlhn7MoUOHUFJSgqqqKmRk+HNNPPnkk3j88cc7TfhutxtlZWWYNGkSMjMzAzUXMSLpCIIg9DKZmZnYunWrNtl3UFdX16nuoUOHoCgKLr744h7vh6zwBwDXO+7Tylosm1ZmRcMMXVQmfyisjvGrY7SdxcbhsWxUlkdAYVKPmkLykWp2amVDJTlLine truncated
"text/plain": [
"<Figure size 700x400 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAF2CAYAAABtfCOsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABgEElEQVR4nO29fXwU5bn//5ndLMkm5IkQkxAwIEqbA1YaiMUQlYh8UaxPEUU41KDwrVC/IrUGsa2/gKLViggtKsfyMuDhgFCFeo4NpHgoCkqogFBJgNYSnkKWh4SHkGSX7M78/gg79zXZnezsZpPsZq+3r3l57+w999yzm7245zPXg6QoigKGYRimx2Pq7gkwDMMwXQMbfIZhmAiBDT7DMEyEwAafYRgmQmCDzzAMEyGwwWcYhokQ2OAzDMNECGzwGYZhIgQ2+AzDMBECG3yGYZgu4LvvvsP48ePRu3dvXHvttXjjjTd0+37zzTf40Y9+hNjYWOTm5mLPnj1BmQMbfIZhmE5GlmXcc889SE1NxTfffIPly5dj4cKFWLNmjUffxsZGTJgwAbfeeiv27NmDvLw83HPPPWhsbOzwPNjgMwzDdDKnT5/G8OHD8e677+KGG27AhAkTMHbsWOzYscOj77p162C1WvHGG28gOzsbS5YsQXx8PP74xz92eB5s8BmGYTqZjIwMrFu3DvHx8VAUBV9++SW++OILjBkzxqNvRUUF8vPzIUkSAECSJIwePRo7d+7s8DyiOjxCD0aWZZw6dQrx8fHqh88wkYaiKGhoaEC/fv1gMvWsNaLdbseVK1cCPl5RFA/bEB0djejoaN1jBg4ciOPHj+PHP/4xHnroIY/3a2trMXToUM2+tLQ0HDhwIOB5umGD3w6nTp3CgAEDunsaYced0kS1/ZnyUTfOhAkmJ06cQP/+/bt7GkHDbrdjUFZv2M64Ah6jd+/euHz5smZfSUkJ5s+fr3vMxx9/DJvNhlmzZuHnP/85fve732neb2pq8vgHIzo6Gg6HI+B5umGD3w7x8fEAWv/QExISOv189fX1GDRoEKqrq9GnT59OP19n8VDSdLV98eJFQ8f0lGsPhFC/9kuXLmHAgAHq76GncOXKFdjOuFC9JwsJ8f7fuVxqkDFoxDEP+9De6h4ARo4cCaD1H5x///d/x6JFi9CrVy/1/ZiYGA/j7nA4EBsb6/cc28IGvx3ct2oJCQldYvBbWloAtP5D0xXn6wqMXkdPvHajhMu191RZMyHeFJDBV483YB9Onz6NnTt34oEHHlD3/du//RuuXLmCS5cuoW/fvur+zMxM2Gw2zfE2mw0ZGRkBz9FNzxLkmG5jnOlhdaPcdv8buO1+fX9jhuluXIoc8GaU6upqFBYWoqamRt23Z88epKamaow9AIwaNQpfffUV3MUI3Q95R40a1eFrZYPPMExEI0MJeDNKbm4uRowYgSeeeAJVVVUoKytDcXExfvWrXwFoXcE3NzcDACZOnIgLFy5gzpw5qKqqwpw5c9DY2IhHHnmkw9fKBj+EiI6OxqRJk3xqgD0RvvbIvPZQQO7Af0Yxm8345JNPEBcXh1tuuQUzZszA7NmzMXv2bADCbRNolYg+/fRTbN++HSNGjEBFRQXKysoQFxfX4WtlDT+EiI6OxuTJk8Pmh99WvlGRxDqi13ljLm/hdu3BJJKvPRRwKQpcivHVOj3OH/r164cNGzZ4fU9pM9bNN9+MvXv3+j0nX7DBZxgmovFXnqHHhRss6TAMw0QIvMJn/OLuAc+IF0S6McUIOcKUInzJlW+ru2ReDBMoMhS4ImSFzwafYZiIJpIkHTb4DMNENF310DYUYIPP+CTnybfUdtL3RbSfJS1ZbZsaRSi4fOKUaDc1AQCGPi/GqHz9550yT4YJBPnqFshx4QYbfIZhIhpXgBp+IMd0N+ylwzAMEyHwCp9RoYFUUdeKNLjyj69V2zW3CW+c5MMiw1/ynrNiIJe42Y1KuwYA0H/pbvH+60GZLsMEBZfSugVyXLjBBp9hmIiGNXyGYZgIQYYEF/xP/SwHcEx3wwY/wtHNh2Myq8307fVq+1yu8My5Ei/+4BuyU9S29H3R7v318dZ9vTue+IlhOgNZad0COS7cYIPPMExE4wpwhR/IMd0Ne+kwDMNECLzCj0D0ZBwTqavpPHZCbZuTE9V2cpzw0rk4RNTYPP89IQH1PeAU42S1eulEnW9S9934rAjC+nYxB2Ex3UskrfDZ4DMME9HIigRZCeChbQDHdDdhIenY7XZMnz4dSUlJyMjIwJtvvum138qVKyFJksdmMonL3Lhxo8f7EydO7KpLYRgmxHCv8APZwo2wWOEXFxdj9+7d2Lp1K44dO4aioiJkZWV5GOpJkybhrrvuUl+3tLTgjjvuwI9//GN1X1VVFe69916899576r6YmJjOv4gwQL4iqlNJZrPXtrmhWW0n/kMcG31BSD0tseIf2Ib+rd45vS5Z1X3WujB0b2B6LC6Y4Apg7evqhLl0NiFv8BsbG7FixQps2rQJOTk5yMnJQWVlJZYtW+Zh8K1WK6xWYVh+85vfQFEUvPbaa+q+gwcPYtiwYUhPT++ya2AYJnRRApR0FJZ0gs/+/fvR0tKCvLw8dV9+fj527doFWdaPdauvr8frr7+O1157TVMrtKqqCkOGDOnUOTMMw4QiIb/Cr62tRd++fdGLeJCkpaXBbrejrq4OqampXo9799130a9fP81dgKIoOHz4MMrLy/Hqq6/C5XLh4YcfxksvvaQZvydipOA4rVolDRygthWS7lg6e160D4tqVnHZg0V/8szEktEq6dhTxJ+auUVIOgOXL1LbSzNuaP8iGKYTYC+dLqS5uRk1NTVe38vIyEBTU5NmhQ5Afe1wOLwdBkVRsGLFCsydO1ez//jx4+p469evR3V1NWbPno3m5mYsXbpUd4719fVoaWnRnL/tnIKB+xz0XJFGJF57qH3vDodD89tqaGjoxtl0Pi7FBJcSgIYfho+iut3g79q1CwUFBV7f27hxI2JiYjwMu/t1bGyst8Owe/dunDx5Eo8++qhmf1ZWFurq6pCcnAxJkjB8+HDIsoypU6di8eLFMJOHk5RBgwZpXk+aNAmTJ082dH2BsGXLlk4bO9Tha+9+1q5di3Xr1nX3NLoMGRLkANRtLnEYAGPGjIHSTqmwr776CufOnYPT6URUVOt0bTYbrFYrkpKSvB6zefNm3HbbbUhOTvZ4r0+fPprX2dnZsNvtqK+v15WHqqurER8fr77uzBX+li1bMG7cOFgslg6PNyF6is8+5jjyjyb5B891kLjgENkHly+L7n3E5ytX/VNtm3r3VtvWY623vcqt4rlJrwtiJfv95Xa1fe/mYnH+lEqfc+8pBPt77yhjx47F22+/rb5uaGjwWPT0JFjSCSGGDx8Oi8WCiooK5OfnAwB27NiB3NxcjX89ZdeuXRg9erTH/vLyckyZMgUnTpxQ7w727duHlJQUXWMPtP4jkZCQEISrMYbFYgmJH353EonXHyrfu8ViQW/yj3YozKkzCVzS4RV+0ImNjUVRURFmzpyJ0tJS1NTUYNGiRSgtLQXQqq+bzWYkJorw/wMHDmDq1KkeY+Xl5cFqtWLGjBkoKSnBkSNHUFxc7KH1hzN6D2fN5A5FuSJW2LJdyGWKU+yPyhBuq07bGa9juurPe92vNAtffWnwQACA9S9/J5MhPy5yV3HHtOlqe9tmr0MzDNMBQt7gA8DixYsxa9YsFBQUIDExEQsWLEBhYSEAoLCwEAMHDsTKlSvV/qdPn/Yq58THx6O8vBxz5szByJEjER8fjyeffBLFxcUefRmGiQxaNXzOhx8yxMbGYtWqVVi1apXHe9u2bfPY10xWmG0ZOnRoyDwcYxim+5EDjLTlh7ZMt6An40Rl9lPbzppTXvtoUijcfKPo/7dvfZ7XTB+au0SguUK8qjQPf93HkYe9VBay/EXUvb3ztlfU9mdf/MrnXBgmUFjDZxiGiRBkmNgtk2EYJhJwKRJcAeTFCeSY7oYNfg/DTLyV9GQcU7TIDmrqnyH6G5BxTCQ5nevCBZ/9pahWlz6FSD563j0Uyz+8R18zDBM4bPAZholoAk+PzJIOwzBMWCErJsgBPLSV+aEt01Xoeea4Ll5U21EpKeINEuwkXxTJsJz/EhkvjSATl9eo1L5kzEuiTQqp0GAuN1QWot499DjnmbNqe+yYV9X2/277pV/zZRhf8AqfYRgmQpAR2ANY/WocoUvIF0BhGIZhggOv8HsAmjw5JOhJbmoSbSLFaAKmHCJbpS4kW2ZUJvHqOendk4bOx3U1lzoN8JLbiYT2humLb9T2TXPeUtv7l/zcr3EYxhuB++GH33qZDT7DMBFN4JG24Wfww2/GDMMwQcSdPC2QzR9qamowceJE9OnTB5mZmXj22Wdht3u/w77//vshSZJm+/TTTzt8rbzCDzO8eedQ6YYGOFEZJeq6gWrbeeSofydVxOMpPRnHRGoCu7yUxNObF5WLvHn0tKXfWpKbZ4nP7gzjk65Y4SuKgokTJyI5ORnbt29HfX09nnjiCZjNZrzxxhse/auqqrB69WqMHTtW3ectA7C/sMFnGCaiCdwt0/gxhw8fRkVFBWw2G9LS0gAAL730Ep577jkPg+9wOFBdXY3c3Fykp6d7Gy5gWNJhGIbpZNLT07F582bV2Lu5SOJm3Bw+fBiSJOG6664L+jx4hR8G6AVZecNMStNRGUVPxqHyCu2vGZN43UhWkYeHBkfRoCmK6WopSSo7aTx2yHFR14gyk3RsivPsObV995Dn1famf7zutT/D+EJWJMiB+OH7cUxSUhLGjx8vjpVlLFu2TCPZuDl48CASExPxk5/8BNu2bcOAAQOwYMEC3H333X7PsS1s8BmGiWgCL4DSesylS5c0+6OjoxEdHd3usXPnzsXevXvx9ddfe7x36NAhNDU1Yfz48Zg3bx42btyIe++9FxUVFRg5cqTf86SLine truncated
"text/plain": [
"<Figure size 700x400 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"\n",
"# Generate 2D cut in momentum transfer and energy\n",
"q1 = np.array([-1.5, 0.0, 0]) # r.l.u.\n",
"q2 = np.array([0, 0.0, 0.0]) # r.l.u.\n",
"EMin = 0.0 # meV\n",
"EMax = 6.0 # meV\n",
"dE = 0.05 # meV\n",
"width = 0.03 # 1/Å orthogonal to cutting direction\n",
"minPixel = 0.03 # 1/Å bin size along cutting direction\n",
"\n",
"ax, *d = ds.plotCutQE(q1=q1, q2=q2, EMin=EMin, EMax=EMax, dE = dE, \n",
" width=width, minPixel=minPixel, vmin=0.0, vmax=3,\n",
" colorbar = True, backgroundSubtraction=False)\n",
"\n",
"ax, *d = ds.plotCutQE(q1=q1, q2=q2, EMin=EMin, EMax=EMax, dE = dE, \n",
" width=width, minPixel=minPixel, vmin=0.0, vmax=3,\n",
" colorbar = True, backgroundSubtraction=True)\n",
"\n",
"# Generate a constant energy cut\n",
"EMin = 4.5 # meV\n",
"EMax = 4.55 # meV\n",
"dqx = 0.03 # 1/Å bin size\n",
"dqy = 0.03 # 1/Å bin size\n",
"\n",
"\n",
"ax, *d = ds.plotQPlane(EMin=EMin,EMax=EMax, xBinTolerance=dqx, yBinTolerance=dqy, \n",
" vmin=0.0, vmax=3,colorbar = True, backgroundSubtraction=False)\n",
"ax, *d = ds.plotQPlane(EMin=EMin,EMax=EMax, xBinTolerance=dqx, yBinTolerance=dqy, \n",
" vmin=0.0, vmax=3,colorbar = True, backgroundSubtraction=True)"
]
},
{
"cell_type": "markdown",
"id": "2e14f805-37bf-434c-ad98-90a701e1b2fd",
"metadata": {},
"source": [
"### Create 1D cuts with overplotting"
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "40c4ef59-1e92-46c9-b27b-bdd8f50cc972",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAADSK0lEQVR4nOzdd3hTZfvA8W+S7pXuCaXsLbOyBEFEBATcuOVVUZwg+nPgBPUFN4p7Ik7kdTBUEEVAlrJn2bsDuvdM8vvj9Jw2TdKm0DZtuT/X1YuMM56UNOfOM+5bZ7FYLAghhBBCiCZP7+oGCCGEEEKIuiGBnRBCCCFEMyGBnRBCCCFEMyGBnRBCCCFEMyGBnRBCCCFEMyGBnRBCCCFEMyGBnRBCCCFEMyGBnRBCCCFEMyGBnRBCCCFEM+HSwC4xMZFrr72W4OBgYmJimDZtGkVFRXa33bZtG/369cPHx4f4+Hi2bNnSwK0VQgghhGjcXBbYWSwWrr32WgoKCvj777/57rvvWLJkCc8884zNtvn5+YwePZrBgwezZcsWBg4cyJgxY8jPz3dBy4UQQgghGiedq2rF7tu3j86dO5OSkkJERAQA3377LY8++iiJiYlW23722We8+OKLHD58GJ1Oh8VioUOHDjz11FNMnDjRBa0XQgghhGh8XNZjFxkZybJly7SgTpWdnW2z7caNG7nooovQ6XQA6HQ6Bg0axIYNGxqkrUIIIYQQTYHLArvAwEBGjhyp3TebzbzzzjsMHz7cZtvk5GSio6OtHouIiODUqVP13k4hhBBCiKbCzdUNUD322GNs3bqVTZs22TxXUFCAp6en1WOenp4UFxc7PF5xcbHV82azmYyMDEJCQrSePyGEEEKIxs5isZCbm0t0dDR6ffV9co0isHv88ceZM2cOCxYsoFu3bjbPe3l52QRxxcXF+Pj4ODzmrFmzmDFjRp23VQghhBDCFU6ePEmLFi2q3cblgd2DDz7I+++/z1dffcU111xjd5uYmBhSUlKsHktJSSEqKsrhcZ988kmmTZum3c/OziY2NpaTJ08SEBBQN40XQgghhKhnOTk5tGzZEn9//xq3dWlgN2PGDD744AO+++47rr32Wofb9e/fn9mzZ2OxWLRVsevWreOpp55yuI+np6fN8C1AQECABHZCCCGEaHKcmUrmssUTCQkJvPDCCzzxxBNcdNFFpKSkaD+g9MgVFhYCcO2115KVlcXUqVPZu3cvU6dOJT8/n+uvv95VzRdCCCGEaHRcFtgtWrQIk8nEiy++SFRUlNUPQFRUFAsWLACUXralS5fy999/06dPHzZu3Mivv/6Kr6+vq5ovhBBCCNHouCxBcUPLycnBaDSSnZ0tQ7FCCNFMmc1mSkpKXN0MIWrF3d0dg8Hg8PnaxDAuXzwhhBBC1IWSkhKOHj2K2Wx2dVOEqLXAwEAiIyPPOSWbBHZCCCGaPIvFQnJyMgaDgZYtW9aY60uIxsJisVBQUMCZM2cAqs344QwJ7IQQQjR5ZWVlFBQUEB0dXW2OUyEaI29vbwDOnDlDeHh4tcOyNZGvNEIIIZo8k8kEgIeHh4tbIsTZUb+QlJaWntNxJLATQgjRbEjJSNFU1dV7VwI7IYQQQohmQgI7IYQQzVZOUSlncooa7Cen6OyH0ebNm4dOp+PTTz+1enzo0KHodDrmz59vs8++ffvQ6XQMHToUgFWrVqHT6bQfNzc34uLiePXVV632i4uLs9qu8s+qVasAmDhxos1zfn5+9OvXjzVr1jh8HWp71R9/f39GjhzJoUOHzvp3U1nlNjY1cXFxzJs3r17PIYsnhBBCNEs5RaXM/fMgGfkNl9cu2NeDB4e3J8DLvdb7fvvtt7Rt25b58+dz5513Wj3n7u7O4sWLue2226we/+mnn+wO4SUnJwNKCpht27Zx6623Ehsby4QJE7Rt5syZY3Vfew3Bwdrt66+/nrfeeku7n5SUxJNPPsn48eM5fvy4w5xqjzzyCI8++igWi4WMjAxeeuklxo8fz+7du2W4vJ5JYCeEEKJZKioxkZFfgqebAR+Ps19l6KyC8vMVlZhqHdidOXOGP//8k88//5zbb7+do0eP0rp1a+35IUOG8Pvvv1NSUmK1QOTnn3+mf//+NseLjIzUbsfGxnLTTTexYMECq0DOaDRabWePt7e31TaRkZF89tlntGjRgr/++ovx48fb3c/Pz0/bLyoqijfeeIOoqCh27drFBRdcUMNvQ5wLGYoVQgjRrPl4GPD1dKv3n3MJHhcuXEhgYCA333wz0dHRNsOuAwcOxMvLi5UrV2qPJSUlcfDgQYYNG1bj8euyBKenpycAbm7O9w1VPX9JSQnTpk0jJiYGd3d34uLi+Oijj7Tn8/PzueeeewgJCSEkJIS7776boqIim+P+/PPPeHt7s3z5cgCOHDnCpZdeio+PD927d+e1114jLi4OUIa6Bw0axFVXXYXRaOTrr7/GbDbz6quv0qZNG7y9vRk2bBi7du3Sjl912HfevHna8VatWkVcXBzvv/8+MTEx+Pr6cuutt1JcXKxt/+GHHxIbG0tAQAAvvvii07+vcyGBnRBCCOFi3333HWPGjEGv1zNu3Djmz59P5Yqfer2eK664gsWLF2uP/fzzz4waNQp39+p7B/fu3ct3333HLbfccs7tzMzM5NFHHyU8PJzBgwc7tU9xcTEvvfQSF1xwAd27dwdg1qxZ/PLLL/zwww/s37+f22+/nQceeIDTp08DcNddd7F27VoWL17MihUrWLt2LU8//bTVcdevX88tt9zCF198wciRIykrK+OKK64gMDCQzZs38+STTzJjxgybfbp27crGjRsZOXIkM2fO5LXXXmPOnDls3bqVVq1acfnll5Ofn+/Ua0tKSuJ///sfy5Yt48cff+SHH37QgvLly5czZcoUXnrpJTZs2MCmTZs4fvy4U8c9FxLYCSGEEC508uRJ1q1bx5VXXgnA1VdfzZEjR1i7dq3VduPHj2fJkiXa/Z9//pmrrrrK7jH9/Pzw8/PDy8uLrl27EhcXx8iRI622mTx5srad+tO1a1erbb7++mvtOV9fX8LDwzl+/DgrVqyotmbpf//7X20/Hx8fXnnlFZ577jltfl2PHj349NNP6d+/P23atGH69OmUlpZy4MABMjMzWbhwIe+++y6DBg2id+/efPjhh7Rq1Uo7/r59+xg7diyvv/46119/PQArV67k5MmTfPbZZ3Tp0oWbbrqJBx54wKpdOp2Op556is6dOxMSEsLcuXN54YUXGDduHJ07d+bjjz/GYDDw1VdfOXxtlZWWlvL222/TvXt3Ro4cyeWXX86mTZsA+OSTT7j55pu59dZb6dq1K5999hleXl5OHfdcSGAnhBBCuNB3332Hl5eXFngNHTqUoKAgvvjiC6vtRowYQXp6Olu3biUrK4uNGzdy+eWX2z3m9u3b2b59Ozt27GD58uUUFRVxxRVXWG0zc+ZMbTv159dff7XaZty4cWzfvp3Nmzczbdo0AgMDeeaZZ2qcJzd58mTtmFu2bOHFF1/khhtu4M8//wTgyiuvpLCwkEceeYQxY8Zow5smk4lDhw5hMpno06ePdrzBgwfz4IMPavenTJlCZmYmsbGx2mM7d+6kQ4cOVgHngAEDrNoVHh5uVeUhIyODfv36ac+7u7vTt29fEhISqn19lbVv3167HRAQoCUY3rt3Lz179tSeCwkJoU2bNk4f92zJ4gkhhBDChb799lsKCwutAhKTycTChQuZO3eu9piPjw8jRoxg8eLFdOjQgaFDh+Ln52f3mO3atdNud+zYEX9/fwYOHMju3bvp1q0boAQ5lbezx9/fX9vmhRdeIDU1lauuuoodO3ZowZg9wcHBVsfu2bMnq1ev5v3332f48OE8/fTTfPzxx/znP//htttu47333tOOV9PQMiiBo4eHBw8++CB79uzB09MTNzc3q+FrwOZ+5R4zR71nJpNJq2RSVVlZmc1jVaudVD5n1fM3RGUU6bETQgghXOTAgQNs27aNt99+26rn7LvvviMnJ4effvrJavvx48ezdOlSFi1a5HAY1h41wHAUsDjrlVdewc/Pj/vuu6/W+1osFu38H3zwAe+88w6zZ89mwoQJ2pw2i8VCmzZtMBgM7NixQ9t30aJF9O7dW7t/1VVX8cwzz5Cfn8/LL78MQNeuXTl48CC5ubnadlu2bHHYHqPRSEREBBs3btQeKy0tZcuWLXTs2BFQArHKxzty5IjTr7dbt27asCxAbm5uneXyq4702AkhhBAu8u233xIcHMzdd9+trTYFJSiYOXOmzXDsFVdcwT333MOhQ4d49913HR43JSVFu52YmMhjjz1Gp06dtMULANnZ2Vbbqfz9/R2uog0ICODVV1/l5ptvZsmSJYwdO9budnl5edqxS0tLWbJkCX/++ac2dy0kJIQlS5bQp08fkpKSmDJlCqAstAgICOD222/noYce4oMPPkCv1zN9+nRGjx5t05ZZs2Zx3333cdtttzF8+HBatmzJpEmTeO6559izZw9vvfWWVV6+qqZNm8azzz5LdHQ07dq14+WXX6aoqEhLCxMfH8/cuXPp1KkTCQkJfP7551b/T9V54IEHuPTSSxkyZAiDBw9mxowZFBQUOLXvuZAeOyGEEM1aQYmJ/OKyev8pKKl9b5i6WtVesHDvvffyxx9/kJiYqD0WHh5Ov3796NmzJ6GhoQ6PGxUVRVRUFNHR0YwYMYIWLVrwyy+/oNdXXPanTp2qbVf5580336y2zTfddBODBw/m4YcftkrtUdnrr7+uHa9du3a8++67vPvuu9x4440AfPbZZ2zfvp2uXbsyceJErr/+ei688EK2bdsGKMmTe/TowYgRIxg1ahTDhg2zmy7k9ttvp3v37jz00EPo9Xp+/PFHEhMT6dmzJy+88AL/+c9/qh3+fOSRR5g0aRKTJk2iT58+nDp1ilWrVhEWFgbA3LlzSU9Pp1u3brzLine truncated
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABaeElEQVR4nO3deXhTZdoG8Psk3ZKWpkJLG8pSKCIWRESLBcoOo8igQ2FEUSkzuIA6gKgMMG7oaBkVBMdxGdEifoqAgIIFmUEsyCaboBRx6cJSGrRlTKUrbd7vj9iQNEmbtElOcnL/uHLRc/Lm5DknafL0XSUhhAARERERBTyV3AEQERERkWcwsSMiIiJSCCZ2RERERArBxI6IiIhIIZjYERERESkEEzsiIiIihWBiR0RERKQQTOyIiIiIFCJE7gB8xWQy4ezZs2jTpg0kSZI7HCIiIiKXCCHw66+/okOHDlCpmq6TC5rE7uzZs+jUqZPcYRARERG1yOnTp9GxY8cmywRNYtemTRsA5osSHR0tczRERERErikvL0enTp0suUxTgiaxa2h+jY6OZmJHREREAceVrmQcPEFERESkEEzsiIiIiBSCiR0RERGRQgRNHzsiIiJyzGQyoba2Vu4wglZoaCjUarVHjsXEjoiIKIjV1taisLAQJpNJ7lCCWkxMDBISElo91y4TOyIioiAlhEBJSQnUajU6derU7OS35HlCCFRWVuKnn34CAOj1+lYdj4kdERFRkKqrq0NlZSU6dOgArVYrdzhBS6PRAAB++ukntG/fvlXNskzNiYiIglR9fT0AICwsTOZIqCGxvnjxYquOw8SOiIgoyHENdfl56jVgYkdEREStUllbh6R5OUial4PK2jq5wwlqTOyIiIgoYK1YsQKSJOGtt96y2T9s2DBIkoSVK1faPebEiROQJAnDhg0DAOTm5kKSJMstJCQESUlJeOGFF2wel5SUZFPO+pabmwsAmDp1qt19UVFRuP7667Fz506vXANrTOyIiIgoYK1atQrJyckOE7jQ0FBs3LjRbv+GDRscNn2WlJSgpKQEBQUFWLZsGZ555hmsXr3apszSpUst5axvAwcOtJS59dZbbe7buXMnYmJicMstt6C8vNwDZ+0cR8USERFRQPrpp5/w2WefITs7G5mZmSgsLETXrl0t9w8ZMgT/+c9/UFtbazNA5KOPPkJaWprd8RISEiw/d+7cGZMnT8bq1asxadIky36dTmdTzhGNRmNTJiEhAW+//TY6duyIzz//HLfcckuLztcVrLEjIiIijzEYq332XGvXrkVMTAzuuOMOdOjQwa7WbuDAgYiIiMD27dst+86ePYsffvgBw4cPb/b4kZGRHos1PDwcABAS4t06NSZ2RERE1CrrDp2x/DxqyQ6sPnDKJ8/7wQcfYOzYsVCpVLj55puxcuVKCCEs96tUKvz+97+3aY796KOPMGbMGISGhjZ57OPHj+ODDz7AnXfe2eo4//e//+GRRx5B+/btMXjw4FYfrylM7IiIiKjFSoxVeHJjnmXbJIAF64+hxFjl1ec9ffo0du/ejT/84Q8AgIyMDBQUFGDXrl025W655RZs2rTJsv3RRx9h/PjxDo8ZFRWFqKgoREREoFevXkhKSsINN9xgU2b69OmWcg23Xr162ZR57733LPdFRkaiffv2OHnyJP773/8iOjraA2fvnF8kdjU1Nejdu7dlRElTioqKEBUV5VJZIiIi8q7C0gqYhO2+eiFQVFrp1ef94IMPEBERYUm8hg0bhssuuwzvvPOOTbnRo0ejrKwMhw8fxi+//IJ9+/bhxhtvdHjMI0eO4MiRIzh69Ci2bt2K6upq/P73v7cp8/TTT1vKNdw2b95sU+bmm2/GkSNHcPDgQcyZMwcxMTF4/PHH0adPHw9eAcdkHzxRXV2NyZMnIy8vr/nCAGbMmIGKigovR0VERESu6BobCZUEm+ROLUlIivXuEmWrVq1CVVWVTQ1YfX091q5di3/+85+WfVqtFqNHj8bGjRvRo0cPDBs2DFFRUQ6P2b17d8vPV1xxBdq0aYOBAwfi2LFj6N27NwCgffv2NuUcadOmjaXMM888g59//hnjx4/H0aNHkZSU1NJTdomsNXbHjx9HWloa8vPzXSr/3nvv4ddff/VyVEREROQqvU6DhTdfaopUScBzGb2h12m89pzff/89vvrqK7z88ss2NWcffPABysvLsWHDBpvyt9xyCz755BN8/PHHTpthHWnor9ew9FpLPf/884iKisL999/fquO4QtbEbseOHRg+fDj27t3bbNmysjLMnTsXb7zxhg8iIyIiIldNuLaj5edtc4ZiUmpnrz7fqlWr0LZtW9x7773o3bu35TZp0iSkpKTYNcf+/ve/tzSvjhs3zulxDQaD5Xbo0CHMnTsXPXv2xFVXXWUpYzQabco13JpqTYyOjsYLL7yALVu22PT38wZZE7sZM2bgpZdesix825Q5c+YgMzPTroMiERER+Y8EXYTXn6NhtGrDFCLWZsyYgW3btqG4uNiyr3379rj++uvRt29fxMbGOj2uXq+HXq9Hhw4dMHr0aHTs2BE5OTlQqS6lS7Nnz7aUs7699NJLTcY8efJkDB48GA899BBqampacNaukb2PnSu2bduGXbt24dixYy4/pqamxubCeXumZyIiIvKNb7/91ul9Dz74IB588EG7/Y1Hyz711FOWn4cNG2YzTYozRUVFzZZZsWKF0/u4pBiAqqoq3HfffXj11Veh0bjeXp+VlQWdTme5derUyYtREhEREcnP7xO7/fv3o6CgABMmTLDMCQMAY8aMwfTp050+bv78+TAajZbb6dOnfRUyERFRUNGGhaBo0VgULRoLbVhANAYqlt9f/f79++OHH36w2Xf55Zdj+fLlGD16tNPHhYeHO2x7JyIiIlIqv03sDAYDdDodNBqNw/liEhMT0b59exkiIyIiIvJPftsUq9frsXr1arnDICIiIgoYflNj13g0SlOjU1wZuUJEREQUbPy2xo6IiIiI3MPEjoiIiEghmNgRERFR69RWAE/pzLda50trecLUqVMhSZLTW25urlef398xsSMiIqKAsWzZMpSUlKCkpARLly5Fx44dLdslJSUYOHCg3CHKym8GTxARERE1p2FFqYaf1Wo1EhISZI7Kf7DGjoiIiBShqKgIkiThmWeewWWXXYYHH3wQTz31FIYNG2ZTLikpybKmqxACzzzzDDp06ICYmBiMGzcOp06d8n3wHsLEjoiIiDyn/KzcEWD37t04ePAgZs2a1WzZV155Be+99x7ef/997Nu3D/Hx8fjd736Hixcv+iBSz2NiR0RERK1z5P1LP/+rP3B4pXyxAJg9ezaSk5Nx+eWXN1v2+eefxwsvvIBhw4ahZ8+eeOONN3D+/Hl8+umnPojU89jHjoiIiFrOWAxsmXtpW5iATbOB5JGALlGWkJKSklwqd+HCBZw5cwaTJk2CSnWprquqqgrff/+9l6LzLiZ2RERE1HLn883JnDVRD5wvkC2xi4iIsPwsSZLd/XV1dTb/r127FldccYVNmbZt23oxQu9hUywRERG1XNtkQGqUTkhqoG03eeJpJCwsDL/++qtl+8KFC/jpp58AADExMWjfvj0MBgO6d++O7t27o3Pnzpg7dy6+++47uUJuFSZ2RERE1HK6RGDM85e2JTUwbqlstXWNpaam4ujRo1i7di2+//573HvvvVCr1Zb758yZg7/97W/YtGkTfvjhB9x9993YvXs3evbsKWPULcemWCIiImqdvpOBzY+Yf37gSyC2+UELvjJy5EjMmTPHktDNmTMHZ89eGrn7yCOP4Ndff8W9996L8vJyXHfdddi6dSsuu+wyGaNuOUkIIeQOwhfKy8uh0+lgNBoRHR0tdzhERESyq66uRmFhIbp27WrTL81ttRXAcx3MPy84C4RFeibAINLUa+FODsOmWCIiIiKFYGJHREREpBDsY0dEREStExYJPGWUOwoCa+yIiIiIFIOJHREREZFCMLEjIiIKckEyQYZf89RrwMSOiIgoSDVM1FtbWytzJFRZWQkACA0NbdVxOHiCiIgoSIWEhECr1eLnn39GaGgoVCrW9/iaEAKVlZX46aefEBMTY7MqRkswsSMiIgpSkiRBr9ejsLAQJ0+elDucoBYTE4OEhIRWH4eJHRERURALCwvD5ZdfzuZYGYWGhra6pq4BEzsiIqIgp1KpWrekGPkNNqYTERERKQQTOyIiIiKFYGJHREREpBBM7IiIiIgUgokdERERkUIwsSMiIiJSCCZ2RERERArBxI6IiIhIIZjYERERESkEEzsiIiIihWBiR0RERKQQTOyIiIiIFMIvEruamhr07t0bubm5Tsvk5OSgb9++iIqKQp8+fbBx40bfBUhEREQUAGRP7Kqrq3H77bcjLy/PaZmvv/4aGRkZ+POf/4wjR47gvvvuw8SJE3H06FEfRkpERETk30LkfPLjx49j8uTJEEI0We7999/HiBEjMHPmTABA9+7dsXHjRqxZswZXX321L0IlIiIi8nuyJnY7duzA8OHD8eyzzyIyMtJpuczMTNTW1trtNxqN3gyPiIiIKKDImtjNmDHDpXJXXnmlzXZeXh4+++wzTJ8+3RthEREREQUkWRO7ligtLcWECRMwaNAg3HLLLU7L1dTUoKamxrJdXl7ui/CIiIiIZCP74Al3nDt3DiNGjIDJZMKHH34Ilcp5+FlZWdDpdJZbp06dfBgpERERke8FTGJXXFyMIUOGoKamBrm5uYiLi2uy/Pz582E0Gi2306dP+yhSIiIiInkERFNsRUUFbrzxRqhUKnz++edISEho9jHh4eEIDw/3QXRERERE/sFvEzuDwQCdTgeNRoPnnnsO+fn5lgmMDQYDAECj0UCn08kYJREREZH/8NumWL1ej9WrVwMA1q1bh6qqKlx//fXQ6/WW26xZs2SOkoiIiMh/+E2NXeNJiq23T5wLine truncated
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Q1 = np.array([-1.0, 0, 0])\n",
"\n",
"ax1D001 = DataSet.generate1DAxisE(Q1,rlu=True,showQ = False)\n",
"ax1D001SUB = DataSet.generate1DAxisE(Q1,rlu=True,showQ = False)\n",
"\n",
"titles = [\"AMBER\",\"True\"]\n",
"fgDone = False\n",
"\n",
"\n",
"\n",
"for dds,title in zip([ds,ds_calculated],titles): \n",
" \n",
" \n",
" #title = bg*'Subtracted'+(1-bg)*'Foreground'#'['+', '.join(['{:.0f}'.format(x) for x in q])+'] '+bg*'Subtracted'\n",
" if dds == ds_calculated:\n",
" ax,data,bins = dds.plotCut1DE(q=Q1, E1 = -0.2, E2 =7.0, width = 0.1, minPixel=0.05,label='Subtracted',ax=None)#,backgroundSubtraction=Fal,plotForeground=True)\n",
" else:\n",
" ax,data,bins = dds.plotCut1DE(q=Q1, E1 = -0.2, E2 =7.0, width = 0.1, minPixel=0.05,label='Subtracted',ax=None,backgroundSubtraction=True,plotForeground=True)\n",
"\n",
" plt.close(ax.get_figure())\n",
" \n",
" pos = data['BinDistance'].to_numpy()\n",
"\n",
" if dds != ds_calculated:\n",
" bgvalue = data['Int_Bg'].to_numpy()\n",
" bgerr = data['Int_Bg_err'].to_numpy()\n",
" fgvalue = data['Int_Fg'].to_numpy()\n",
" fgerr = data['Int_Fg_err'].to_numpy()\n",
" subvalue = data['Int'].to_numpy()\n",
" suberr = data['Int_err'].to_numpy()\n",
" else:\n",
" fgvalue = data['Int']\n",
" fgerr = data['Int_err']\n",
" subvalue = fgvalue\n",
" suberr = fgerr\n",
" \n",
"\n",
" #if not fgDone:\n",
" ax1D001.errorbar(pos,fgvalue,yerr=fgerr,label='Foreground',fmt='.')\n",
" # fgDone = True\n",
" \n",
" c = ax1D001SUB.errorbar(pos,subvalue,yerr=suberr,label=title,fmt='.')[0]\n",
" color = c.get_color()\n",
" \n",
" ax1D001.plot(pos,bgvalue,color='k',zorder=-10)\n",
" ax1D001.fill_between(pos,bgvalue-bgerr,bgvalue+bgerr,color=color,label=title+' Background',alpha=0.5,zorder=-11)\n",
"\n",
"for ax in [ax1D001]:\n",
" ax.set_ylim(-0.3,2)\n",
" ax.legend()\n",
" \n",
" ax.set_ylabel('Int [arb. unit]')\n",
" ax.get_figure().tight_layout()\n",
" \n",
"for ax in [ax1D001SUB]:\n",
" ax.set_ylim(-0.15,1.5)\n",
" ax.legend()\n",
" ax.set_ylabel('Int [arb. unit]')\n",
" ax.get_figure().tight_layout()"
]
}
],
"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.12.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}