From ea77c5993190c2d6aa3314b59b319ca3109e29e0 Mon Sep 17 00:00:00 2001 From: tegenolf <138484196+tegenolf@users.noreply.github.com> Date: Wed, 22 Jan 2025 22:21:08 +0100 Subject: [PATCH] Upload notebook for week 9 --- week9/nmap9.ipynb | 1312 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1312 insertions(+) create mode 100644 week9/nmap9.ipynb diff --git a/week9/nmap9.ipynb b/week9/nmap9.ipynb new file mode 100644 index 0000000..ffaedf5 --- /dev/null +++ b/week9/nmap9.ipynb @@ -0,0 +1,1312 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "incomplete-medline", + "metadata": {}, + "source": [ + "<h1 style=\"text-align: center; vertical-align: middle;\">Numerical Methods in Accelerator Physics</h1>\n", + "<h2 style=\"text-align: center; vertical-align: middle;\">Python examples -- Week 9</h2>" + ] + }, + { + "cell_type": "markdown", + "id": "executive-television", + "metadata": {}, + "source": [ + "<h2>Run this first!</h2>\n", + "\n", + "Imports and modules:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "biological-product", + "metadata": {}, + "outputs": [], + "source": [ + "from config9 import (np, plt, radon, iradon, Image, \n", + " tnrange, plot_mp, plot_tomo)\n", + "from scipy.constants import m_p, e, c\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "58404467-c2fd-491c-87f5-42848bd0e7d9", + "metadata": {}, + "source": [ + "<h2>Projection integral: Radon transform</h2>\n", + "\n", + "Load sample image for the tomography:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17ff62e5-9873-49c3-ad6d-239520c31b01", + "metadata": {}, + "outputs": [], + "source": [ + "data = ~np.array(Image.open('src/temf.png').convert('1', dither=False))\n", + "\n", + "plt.imshow(data, cmap='binary');" + ] + }, + { + "cell_type": "markdown", + "id": "d4fe311f-0256-4121-8498-db4126f08a2f", + "metadata": {}, + "source": [ + "Compute the Radon transform at an angle of 0 deg and 90 deg:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0b0dd3e-8c3d-4ed1-b59c-ca8c3ee9e8bc", + "metadata": {}, + "outputs": [], + "source": [ + "Rf_0 = radon(data, [0], circle=False).astype(float)\n", + "Rf_90 = radon(data, [90], circle=False).astype(float)\n", + "\n", + "plt.plot(Rf_0, label='0 deg')\n", + "plt.plot(Rf_90, label='90 deg')\n", + "plt.legend()\n", + "plt.xlabel('$p$')\n", + "plt.ylabel(r'$\\mathcal{R}_{\\theta}(p)f$');" + ] + }, + { + "cell_type": "markdown", + "id": "7e0881e0-c223-4fdd-9ebc-dc568bd7b2c9", + "metadata": {}, + "source": [ + "<h3>Preparing the measurement data</h3>\n", + "\n", + "We take a number of `VIEW` measurements across the angular interval of [0,`ANG`] degrees:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2ba1d99-b3b9-412c-9c69-8077673d28d3", + "metadata": {}, + "outputs": [], + "source": [ + "# parameters\n", + "M = max(data.shape)\n", + "ANG = 180\n", + "VIEW = 180\n", + "THETA = np.linspace(0, ANG, VIEW, endpoint=False)" + ] + }, + { + "cell_type": "markdown", + "id": "316e27f5-6dab-47a2-b329-ab294be233bd", + "metadata": {}, + "source": [ + "`A` is the projection operator, applying the Radon transform along all chosen angles `THETA` to the original object under study, $x$ (`data`):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a7864df-c5d8-4816-ad8e-66871f9b548f", + "metadata": {}, + "outputs": [], + "source": [ + "A = lambda x: radon(x, THETA, circle=False).astype(float)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "102a331c-cc4d-43d3-a8a1-32fcf39580bf", + "metadata": {}, + "outputs": [], + "source": [ + "proj = A(data)" + ] + }, + { + "cell_type": "markdown", + "id": "f9988137-6209-48cb-93a0-45800eae6bf2", + "metadata": {}, + "source": [ + "<h3>Sinogram</h3>\n", + "\n", + "The sinogram represents the measurement, the collection of all projection results:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8efbfb19-cbb0-4574-b3c9-8b538a2fb11c", + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(proj.T)\n", + "\n", + "plt.gca().set_aspect('auto')\n", + "plt.xlabel('Projection location [px]')\n", + "plt.ylabel('Rotation angle [deg]')\n", + "plt.colorbar(label='Intensity');" + ] + }, + { + "cell_type": "markdown", + "id": "07505825-224c-43ca-9718-1930f4d55e51", + "metadata": {}, + "source": [ + "<h2>Using Filtered Back Projection</h2>\n", + "\n", + "Let's apply the back projection as an inverse Radon transform (implemented with a ramp filter):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b67a7197-3eae-4442-8c19-71d74820d72b", + "metadata": {}, + "outputs": [], + "source": [ + "AINV = lambda b: iradon(b, THETA, circle=False, filter_name='ramp', output_size=M).astype(float)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40384134-bcf0-4e3a-b7d3-142c557558d3", + "metadata": {}, + "outputs": [], + "source": [ + "fbp = AINV(proj)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6458f51d-5539-45c3-b332-19d2e8b7f810", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(12, 6))\n", + "plt.imshow(fbp, cmap='binary', vmin=0, vmax=1);" + ] + }, + { + "cell_type": "markdown", + "id": "58261674-9f7a-4c5c-9742-e2d055c56fa4", + "metadata": {}, + "source": [ + "$\\implies$ note the artifacts due to edges in the reconstructed image (high frequency!)." + ] + }, + { + "cell_type": "markdown", + "id": "e76dfda1-1f9c-41e1-9ffc-599e85b1ac6e", + "metadata": {}, + "source": [ + "<h2>Effect of Noise</h2>\n", + "\n", + "In the following, we will compare the FBP with the iterative ART approach on a noisy sinogram.\n", + "\n", + "We add 15% of the maximum sinogram value as Gaussian normal distributed noise:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "46c6c4eb-eae7-423a-8c06-e69031256664", + "metadata": {}, + "outputs": [], + "source": [ + "noise = np.random.normal(0, 0.15 * np.max(proj), size=proj.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a803dc3-eb60-467c-9e42-a86c454eb1bd", + "metadata": {}, + "outputs": [], + "source": [ + "proj_noise = proj + noise" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "332f1aa0-b2ae-4902-819e-d4123a5a0a88", + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(proj_noise.T)\n", + "plt.gca().set_aspect('auto')\n", + "plt.xlabel('Projection location [px]')\n", + "plt.ylabel('Rotation angle [deg]')\n", + "plt.colorbar(label='Intensity');" + ] + }, + { + "cell_type": "markdown", + "id": "66d0a42c-be38-424f-9a9d-fdcda21d33e4", + "metadata": {}, + "source": [ + "<h3>Filtered Back Projection Results</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c21dcf8d-300f-443e-9da0-5bbdfe662fbb", + "metadata": {}, + "outputs": [], + "source": [ + "fbp_noise = AINV(proj_noise)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "afbb978c-4195-4cf2-b908-0195c8f11aba", + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(fbp_noise, cmap='binary', vmin=0, vmax=1, interpolation='None');" + ] + }, + { + "cell_type": "markdown", + "id": "0dd794e4-b190-40a0-af26-856938ff295a", + "metadata": {}, + "source": [ + "<h3>Implement the Algebraic Reconstruction Technique</h3>\n", + "\n", + "Besides the projection operator $A$ (`A`), we need the adjoint $A^\\mathrm{T}$ (`AT`) for the ART implementation:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6164b61-0f72-41b8-a72a-a3fa2d2a5e6e", + "metadata": {}, + "outputs": [], + "source": [ + "AT = lambda b: iradon(b, THETA, circle=False, filter_name=None, output_size=M).astype(float) * 2 * len(THETA) / np.pi" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a6db6a7-f583-4242-a8f0-8571fabda1a5", + "metadata": {}, + "outputs": [], + "source": [ + "def ART(A, AT, b, x, mu=1, niter=10):\n", + " # for all i: ||a_i||^2 = A^T ⋅ A\n", + " ATA = AT(A(np.ones_like(x)))\n", + "\n", + " for i in tnrange(niter):\n", + " x = x + np.divide(mu * AT(b - A(x)), ATA)\n", + "\n", + " # nonlinearity: constrain to >= 0 values\n", + " x[x < 0] = 0\n", + "\n", + " plt.imshow(x, cmap='binary', vmin=0, vmax=1, interpolation='None')\n", + " plt.title(\"%d / %d\" % (i + 1, niter))\n", + " plt.show()\n", + "\n", + " return x\n", + "\n", + "# initialization\n", + "x0 = np.zeros((M, M))\n", + "mu = 1\n", + "niter = 30" + ] + }, + { + "cell_type": "markdown", + "id": "5d21b2e3-fa46-4ee4-8262-7d345e194d61", + "metadata": {}, + "source": [ + "<h3>ART Results</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf1dd707-014e-44d7-85a4-9caeaedceb6a", + "metadata": {}, + "outputs": [], + "source": [ + "x_art = ART(A, AT, proj_noise, x0, mu, niter)" + ] + }, + { + "cell_type": "markdown", + "id": "d4fafe0b-04c1-4a1c-a8a0-d2ca11c2c70c", + "metadata": {}, + "source": [ + "<h3>Comparison of the Results</h3>\n", + "\n", + "ART typically outperforms FBP in terms of reconstruction quality:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18a53b6a-6e3e-46c6-bb6a-90f1837c98d3", + "metadata": {}, + "outputs": [], + "source": [ + "_, ax = plt.subplots(1, 2, figsize=(10, 5))\n", + "plt.subplots_adjust(wspace=0.3)\n", + "\n", + "plt.sca(ax[0])\n", + "plt.imshow(fbp_noise, cmap='binary', vmin=0, vmax=1, interpolation='None')\n", + "plt.title('FBP')\n", + "\n", + "plt.sca(ax[1])\n", + "plt.imshow(x_art, cmap='binary', vmin=0, vmax=1, interpolation='None')\n", + "plt.title('ART')" + ] + }, + { + "cell_type": "markdown", + "id": "67992111-5a01-4f78-9bea-a1232985104d", + "metadata": {}, + "source": [ + "<h2>Phase-space Tomography</h2>\n", + "\n", + "The `longitudinal_tomograpy` package is developed at CERN: control room apps tomographically reconstruct the longitudinal phase-space distribution from bunch profile measurements!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b7a23b9-4a4c-4ece-896a-61fed6f32ded", + "metadata": {}, + "outputs": [], + "source": [ + "import longitudinal_tomography.tomography.tomography as tmo\n", + "from longitudinal_tomography.tracking import particles, machine\n", + "\n", + "import longitudinal_tomography.utils.tomo_input as tomoin\n", + "import longitudinal_tomography.utils.tomo_output as tomoout\n", + "from longitudinal_tomography.tracking import tracking" + ] + }, + { + "cell_type": "markdown", + "id": "dc52aa24-2814-4e0d-8839-1c3c5660873a", + "metadata": {}, + "source": [ + "We work once more with the `PyHEADTAIL` library to generate the longitudinal macro-particle distribution.\n", + "\n", + "If `PyHEADTAIL` is not installed, please run this:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7b7dce29-6d01-4f00-99df-388b0e101513", + "metadata": {}, + "outputs": [], + "source": [ + "!pip install PyHEADTAIL==1.16.4" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f6b61ff-43a2-44cf-8c94-a9ec17110e71", + "metadata": {}, + "outputs": [], + "source": [ + "from PyHEADTAIL.particles.generators import RFBucketMatcher, ThermalDistribution\n", + "from PyHEADTAIL.trackers.rf_bucket import RFBucket" + ] + }, + { + "cell_type": "markdown", + "id": "50599adf-6f5b-4e37-94c8-d2b07fb60a79", + "metadata": {}, + "source": [ + "<h3>Beam Measurements</h3>\n", + "\n", + "Consider a circulating bunch in a synchrotron or storage ring. Longitudinal bunch profiles can be recorded via wall current monitor with a high-bandwidth oscilloscope: store $V_\\mathrm{gap}(t)$ during the bunch passage time $t$.\n", + "\n", + "<h3>Physical Beam Parameters</h3>\n", + "\n", + "During previous lectures we discussed reduced models for the longitudinal plane. We have\n", + "- `voltage`: rf voltage\n", + "- `harmonic`: rf frequency divided by revolution frequency\n", + "- `circumference`: accelerator ring circumference\n", + "- `gamma`: Lorentz gamma of bunch\n", + "- `alpha_c`: momentum compaction of the ring" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c168076-fbb7-49f0-a3bc-23fb77c8ecdb", + "metadata": {}, + "outputs": [], + "source": [ + "voltage = 24e3\n", + "harmonic = 7\n", + "circumference = 100 * 2 * np.pi\n", + "\n", + "gamma = 3.1\n", + "beta = np.sqrt(1 - gamma**-2)\n", + "\n", + "alpha_c = 0.027\n", + "eta = alpha_c - gamma**-2" + ] + }, + { + "cell_type": "markdown", + "id": "89338727-8ec4-4845-9407-7367a8805a65", + "metadata": {}, + "source": [ + "We look at a circuating proton bunch in a CERN Proton Synchrotron like machine. A \"full bunch length\" of $B_L = 180$ns translates to an rms bunch length in meters of $\\sigma_z=B_L/4\\cdot\\beta c$. As before, the following `PyHEADTAIL` `RFBucket` class captures most of the important quantities for computation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "880d296e-a008-4dc9-b28b-144eeff70742", + "metadata": {}, + "outputs": [], + "source": [ + "rfb = RFBucket(circumference, gamma, m_p, e, [alpha_c], 0., [harmonic], [voltage], [np.pi])\n", + "\n", + "sigma_z = 180e-9 * beta * c / 4. # in [m]" + ] + }, + { + "cell_type": "markdown", + "id": "c12cc674-14a3-45a4-a672-9c0876e8cdbe", + "metadata": {}, + "source": [ + "<h3> A. The simulated measurement</h3>\n", + "<h4>Initialisation of particles</h4>\n", + "\n", + "A thermal distribution of $N$ particles is matched into the rf bucket." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6bb10960-61f5-4eae-acd9-95624c43ffd2", + "metadata": {}, + "outputs": [], + "source": [ + "N = 100000" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84720f8d-90d6-4168-9a4a-17238e56ccce", + "metadata": {}, + "outputs": [], + "source": [ + "np.random.seed(12345)\n", + "\n", + "rfb_matcher = RFBucketMatcher(rfb, ThermalDistribution, sigma_z=sigma_z)\n", + "rfb_matcher.integrationmethod = 'cumtrapz'\n", + "\n", + "z, dp, _, _ = rfb_matcher.generate(N)" + ] + }, + { + "cell_type": "markdown", + "id": "0145fe43-fd41-4798-8092-423195364324", + "metadata": {}, + "source": [ + "In case `PyHEADTAIL` is not available, import the generated distribution directly (uncomment one of the following lines):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed3c2dc7-fa68-42ac-8bf0-369b41ed3780", + "metadata": {}, + "outputs": [], + "source": [ + "# z, dp = np.load('data/long-dist-80ns.npy') # short bunch\n", + "# z, dp = np.load('data/long-dist-180ns.npy') # long bunch" + ] + }, + { + "cell_type": "markdown", + "id": "47dad048-2c11-4ab4-affa-91db0bd36d01", + "metadata": {}, + "source": [ + "The distribution in longitudinal phase space ($z$, $\\delta$) looks as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d787a0f8-7997-4795-8f38-ebb4de24a111", + "metadata": {}, + "outputs": [], + "source": [ + "plot_mp(z, dp, rfb);" + ] + }, + { + "cell_type": "markdown", + "id": "dd028fbb-a578-49e2-8c9b-652c98588ecc", + "metadata": {}, + "source": [ + "Let's have a look at the bunch profile as it would appear in a discrete measurement (e.g. via a wall current monitor)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "54f56b56-a024-42d1-a534-1ff657b628ec", + "metadata": {}, + "outputs": [], + "source": [ + "nbins = 100\n", + "z_bins = np.linspace(*rfb.interval, num=nbins-1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7f8f984-c917-441b-88c6-5a67168d6d8e", + "metadata": {}, + "outputs": [], + "source": [ + "plt.hist(z, bins=z_bins, histtype='stepfilled')\n", + "plt.xlabel('$z$ [m]')\n", + "plt.ylabel('Discrete bunch profile');" + ] + }, + { + "cell_type": "markdown", + "id": "986ea0d4-9138-4fb5-8f71-9fd382613635", + "metadata": {}, + "source": [ + "A single bin of the profile measurement amounts to a time (in seconds) of" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc8747e8-cd7b-4c1f-8bdc-ab97ff806a88", + "metadata": {}, + "outputs": [], + "source": [ + "dtbin = (z_bins[1] - z_bins[0]) / (beta * c)\n", + "dtbin" + ] + }, + { + "cell_type": "markdown", + "id": "322e852e-5fe7-4411-a035-e6cc40dca5b8", + "metadata": {}, + "source": [ + "The synchrotron period `T_s`, i.e. the time period of the longitudinal motion, amounts to the following number of turns:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5db3c2ac-8bcb-4ba2-8a6a-159ceb7de70d", + "metadata": {}, + "outputs": [], + "source": [ + "T_s = int(1 / rfb.Q_s)\n", + "T_s" + ] + }, + { + "cell_type": "markdown", + "id": "a67c605b-725f-4339-a39c-06ecd58ead0f", + "metadata": {}, + "source": [ + "<h4>Deliberate mismatch</h4>\n", + "\n", + "To provide some interesting dynamics for the tomographic reconstruction, let's mismatch the distribution in momentum. This will launch a quadrupolar oscillation, i.e. the bunch length and momentum spread will oscillate during the synchrotron period." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2687813-e761-4571-9ee7-4a4730c9d581", + "metadata": {}, + "outputs": [], + "source": [ + "dp *= 0.3" + ] + }, + { + "cell_type": "markdown", + "id": "4bad2836-9110-4f0f-bbe6-40451b9e8f74", + "metadata": {}, + "source": [ + "As one can see, the distribution does not follow the iso-Hamiltonian lines (red) any longer but is slightly compressed along the momentum $\\delta$ axis:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba4eea27-7fa8-4ad0-aebc-dafd4275cd23", + "metadata": {}, + "outputs": [], + "source": [ + "plot_mp(z, dp, rfb);" + ] + }, + { + "cell_type": "markdown", + "id": "5181b045-3d9a-4266-aec7-3798b43505c9", + "metadata": {}, + "source": [ + "<h4>Particle Tracking</h4>\n", + "\n", + "We model the synchrotron motion of the particles due to the rf cavities once more via a second order leap-frog integrator. The `track` function advances the particles by one turn:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "505d1c30-2a1d-48dd-9714-c5e7f1325bb5", + "metadata": {}, + "outputs": [], + "source": [ + "def track(z, dp, rfb=rfb):\n", + " # half drift\n", + " z = z - eta * dp * circumference / 2\n", + " # rf kick\n", + " amplitude = rfb.charge * voltage / (beta * c * rfb.p0)\n", + " phi = harmonic * (2 * np.pi * z / circumference) + rfb.phi_offset_list[0]\n", + " dp += amplitude * np.sin(phi)\n", + " # half drift\n", + " z = z - eta * dp * circumference / 2\n", + " return z, dp" + ] + }, + { + "cell_type": "markdown", + "id": "e71861bb-c215-45fa-afb1-48706102ee5c", + "metadata": {}, + "source": [ + "Let's gather the data for the tomographic reconstruction by recording a bunch profile every few turns during one `T_s`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db897c47-ae66-4dfb-97c8-2bc7680f9efc", + "metadata": {}, + "outputs": [], + "source": [ + "record_every_nturns = 10" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d2d43f6-e358-4050-8fd6-6806ed1773f4", + "metadata": {}, + "outputs": [], + "source": [ + "raw_data = [np.histogram(z, bins=z_bins)[0]]\n", + "\n", + "for i in tnrange(1, T_s + 1):\n", + " z, dp = track(z, dp)\n", + " if not i % record_every_nturns:\n", + " # the discrete WCW measurement:\n", + " raw_data += [np.histogram(z, bins=z_bins)[0]]" + ] + }, + { + "cell_type": "markdown", + "id": "6866540e-3289-4117-83d2-52f8fa0bd58a", + "metadata": {}, + "source": [ + "The quadrupole oscillation is clearly visible:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13d4c0ad-c2bf-49d5-acde-da11290dae65", + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(raw_data)\n", + "plt.xlabel('Profile bin')\n", + "plt.ylabel('#Profile')\n", + "plt.colorbar(label='Density');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "be6d3b02-a46a-498b-96b6-1822da1e7d00", + "metadata": {}, + "outputs": [], + "source": [ + "std=[np.sqrt(np.average((np.array((z_bins[:-1]+z_bins[1:])/2) - np.average(np.array((z_bins[:-1]+z_bins[1:])/2),weights=raw_data[i]))**2, weights=raw_data[i])) for i in range(len(raw_data))]\n", + "plt.plot(std)\n", + "plt.xlabel('#profile')\n", + "plt.ylabel('rms bunch length [m]');" + ] + }, + { + "cell_type": "markdown", + "id": "dad78ee3-fc3c-436a-a309-af065f2bc225", + "metadata": {}, + "source": [ + "$\\implies$ a $\\gtrsim 100\\%$ rms bunch length (quadrupole) oscillation takes place!\n", + "\n", + "<h4>Ready...</h4>\n", + "\n", + "We have now gathered a simulated discrete WCW measurement in the `raw_data` array during one synchrotron period." + ] + }, + { + "cell_type": "markdown", + "id": "246fb7e7-0b6e-4355-8d5d-263709a8d433", + "metadata": {}, + "source": [ + "<h3>B. The longitudinal tomographic reconstruction</h3>\n", + "<h4>Preparations</h4>\n", + "\n", + "The `longitudinal_tomography` package requires a certain inout format for the machine parameters (voltage, frequency, etc.):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b1655988-3fa3-4dbc-bdc4-a9cd7cbff27b", + "metadata": {}, + "outputs": [], + "source": [ + "frame_input_args = {\n", + " 'raw_data_path': './',\n", + " 'framecount': len(raw_data), # Number of frames in input data\n", + " 'skip_frames': 0, # Number of frames to ignore\n", + " 'framelength': len(raw_data[0]), # Number of bins in each frame\n", + " 'dtbin': dtbin, # Width (in s) of each frame bin\n", + " 'skip_bins_start': 0, # Number of frame bins before the lower profile bound to ignore\n", + " 'skip_bins_end': 0, # Number of frame bins after the upper profile bound to ignore\n", + " 'rebin': 1, #Number of frame bins to rebin into one profile bin\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43365e28-886b-43ac-9823-025853a4b8ff", + "metadata": {}, + "outputs": [], + "source": [ + "# computation of dipole B-field from beam rigidity Brho = p/e\n", + "Ekin = (gamma - 1) * m_p * c**2 / e\n", + "p0c = np.sqrt(Ekin**2 + 2*Ekin*m_p/e * c**2)\n", + "brho = p0c / c\n", + "brho / 70.08" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d381a56-b93a-44ec-942d-a48965a066c3", + "metadata": {}, + "outputs": [], + "source": [ + "machine_args = {\n", + " 'output_dir': '/tmp/', # Directory in which to write all output\n", + " 'dtbin': dtbin, # Width (in s) of each frame bin\n", + " 'dturns': record_every_nturns, # Number of machine turns between frames\n", + " 'synch_part_x': frame_input_args['framelength']/2, # Time (in frame bins) \n", + " # from the lower profile bound to the synchronous phase (if <0, \n", + " # a fit is performed) in the \"bunch reference\" frame\n", + " 'demax': -1e6, # noqa - Max energy (in eV) of reconstructed phase space (if >0)\n", + " 'filmstart': 0, # Number of the first profile at which to reconstruct\n", + " 'filmstop': 1, # Number of the last profile at which to reconstruct\n", + " 'filmstep': 1, # Step between reconstructions\n", + " 'niter': 50, # Number of iterations for each reconstruction\n", + " 'snpt': 4, # Square root of the number of test particles to track per cell\n", + " 'full_pp_flag': False, # Flag to extend the region in phase space of map elements (if =1)\n", + " 'beam_ref_frame': 0, # Reference frame for bunch parameters (synchronous phase, baseline, integral)\n", + " 'machine_ref_frame': 0, # Reference frame for machine parameters (RF voltages, B-field)\n", + " 'vrf1': voltage, # Peak RF voltage (in V) of principal RF system\n", + " 'vrf1dot': 0.0, # and its time derivative (in V/s)\n", + " 'vrf2': 0.0, # Peak RF voltage (in V) of higher-harmonic RF system\n", + " 'vrf2dot': 0.0, # and its time derivative (in V/s)\n", + " 'h_num': harmonic, # Harmonic number of principal RF system\n", + " 'h_ratio': 2.0, # Ratio of harmonics between RF systems\n", + " 'phi12': 0, # Phase difference (in radians of the principal harmonic) between RF systems\n", + " 'b0': brho / 70.08, # Dipole magnetic field (in T) -- up to 1.8T\n", + " 'bdot': 0.0, # and its time derivative (in T/s) -- up to 10T/s\n", + " 'mean_orbit_rad': circumference / (2 * np.pi), # Machine radius (in m)\n", + " 'bending_rad': 70.08, # Bending radius (in m)\n", + " 'trans_gamma': alpha_c**-0.5, # Gamma transition\n", + " 'rest_energy': m_p * c**2 / e, # Rest mass (in eV/c**2) of accelerated particle\n", + " 'charge': 1, # Charge state of accelerated particle\n", + " 'self_field_flag': False, # Flag to include self-fields in the tracking (if =1)\n", + " 'g_coupling': 0.0, # Geometrical coupling coefficient\n", + " 'zwall_over_n': 0.0, # Reactive impedance (in Ohms per mode number) over a machine turn\n", + " 'pickup_sensitivity': 1, # Effective pick-up sensitivity (in digitizer units per instantaneous Amp)\n", + " 'nprofiles': frame_input_args['framecount'],\n", + " 'nbins': frame_input_args['framelength'],\n", + " 'min_dt': 0.0,\n", + " 'max_dt': dtbin * frame_input_args['framelength'],\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0878c450-32b2-478c-aa73-91a433f08aa3", + "metadata": {}, + "outputs": [], + "source": [ + "# initialising tomography\n", + "frames = tomoin.Frames(**frame_input_args)\n", + "\n", + "mach = machine.Machine(**machine_args)\n", + "\n", + "mach.values_at_turns()" + ] + }, + { + "cell_type": "markdown", + "id": "0ea910d3-99e5-4c2e-8d15-0621ac87dd1b", + "metadata": {}, + "source": [ + "The tomography package uses the waterfall plot of profiles as measured sinogram input then:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f68167b7-1ea3-4132-a121-ae6c690eb037", + "metadata": {}, + "outputs": [], + "source": [ + "measured_waterfall = frames.to_waterfall(np.array(raw_data, dtype=float).flatten())\n", + "\n", + "plt.imshow(measured_waterfall)\n", + "plt.xlabel('Profile bin')\n", + "plt.ylabel('#Profile')\n", + "plt.colorbar(label='Density');" + ] + }, + { + "cell_type": "markdown", + "id": "19db4412-bdd1-43e7-a39f-c16795065f50", + "metadata": {}, + "source": [ + "<h4>Settings</h4>\n", + "\n", + "Reconstruction will be carried out at the following profile index:\n", + "(choose `0` for the first profile, or `-2` for the last, or any number < `len(measured_waterfall) - 1`)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed0ea569-2995-41d0-b920-0432503a51bf", + "metadata": {}, + "outputs": [], + "source": [ + "reconstr_idx = 0" + ] + }, + { + "cell_type": "markdown", + "id": "afc02db6-d278-4742-8abd-b6f84579307c", + "metadata": {}, + "source": [ + "Number of iterations of ART:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86947745-f690-469f-b4db-f85c6379fe0f", + "metadata": {}, + "outputs": [], + "source": [ + "niterations = 50" + ] + }, + { + "cell_type": "markdown", + "id": "7c68302d-1277-478d-a100-5106a919342f", + "metadata": {}, + "source": [ + "<h4>Projection Matrix</h4>\n", + "\n", + "Establishing the map via tracking (nonlinear synchrotron motion), i.e. matrix $A$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02b77a6e-4afd-4a27-8c1d-6cad50fac001", + "metadata": {}, + "outputs": [], + "source": [ + "tracker = tracking.Tracking(mach)\n", + "\n", + "phip, dEp = tracker.track(reconstr_idx)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69ba16d9-66ed-4520-a9dd-948a3cdf6bde", + "metadata": {}, + "outputs": [], + "source": [ + "# Converting from physical coordinates ([rad], [eV])\n", + "# to internal phase space coordinates.\n", + "if not tracker.self_field_flag:\n", + " phip, dEp = particles.physical_to_coords(\n", + " phip, dEp, mach, tracker.particles.xorigin,\n", + " tracker.particles.dEbin)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2d729ae0-a678-4e65-a8e9-8503bc99ea81", + "metadata": {}, + "outputs": [], + "source": [ + "phip, dEp = particles.ready_for_tomography(phip, dEp, mach.nbins)\n", + "\n", + "profiles = tomoin.raw_data_to_profiles(\n", + " measured_waterfall, mach, frames.rebin, frames.sampling_time)\n", + "profiles.calc_profilecharge()\n", + "\n", + "waterfall = profiles.waterfall" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b0e30199-31dd-4a12-9052-279fbd4d2577", + "metadata": {}, + "outputs": [], + "source": [ + "# further preparations\n", + "nprofs = waterfall.shape[0]\n", + "nbins = waterfall.shape[1]\n", + "nparts = phip.shape[0]\n", + "\n", + "flat_profs = waterfall.copy()\n", + "flat_profs = flat_profs.clip(0.0)\n", + "\n", + "flat_profs /= np.sum(flat_profs, axis=1)[:, None]\n", + "flat_profs = np.ascontiguousarray(flat_profs.flatten()).astype(float)\n", + "\n", + "waterfall /= np.sum(waterfall, axis=1)[:, None]\n", + "\n", + "flat_points = phip.copy()\n", + "for i in range(nprofs):\n", + " flat_points[:, i] += nbins * i" + ] + }, + { + "cell_type": "markdown", + "id": "2bcf5641-f670-42f0-912a-5b100bd28aa1", + "metadata": {}, + "source": [ + "Starting the algebraic reconstruction technique (ART) algorithm:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "123a3147-bdf3-421a-9eca-38d68d262830", + "metadata": {}, + "outputs": [], + "source": [ + "# Initialising arrays with zeros\n", + "weight = np.zeros(nparts)\n", + "rec_wf = np.zeros(waterfall.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9b059c4b-d157-4f07-ab79-8e255e3c9a94", + "metadata": {}, + "outputs": [], + "source": [ + "# Initial estimation of weight factors using (flattened) measured profiles.\n", + "weight = tmo.libtomo.back_project(weight, flat_points, flat_profs,\n", + " nparts, nprofs)\n", + "weight = weight.clip(0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66a3ccd5-07b9-4021-b81b-59dce7e84bff", + "metadata": {}, + "outputs": [], + "source": [ + "diff = []\n", + "for i in range(niterations):\n", + " # Projection from phase space to time projections\n", + " rec_wf = tmo.libtomo.project(rec_wf, flat_points, weight, nparts,\n", + " nprofs, nbins)\n", + "\n", + " # Normalizing reconstructed waterfall\n", + " rec_wf /= np.sum(rec_wf, axis=1)[:, None]\n", + "\n", + " # Finding difference between measured and reconstructed waterfall\n", + " dwaterfall = waterfall - rec_wf\n", + "\n", + " # Setting to zero for next round\n", + " rec_wf[:] = 0.0\n", + "\n", + " # Calculating discrepancy\n", + " diff.append(np.sqrt(np.sum(dwaterfall**2) / (nbins * nprofs)))\n", + "\n", + " # Back projection using the difference between measured and recorded waterfall\n", + " weight = tmo.libtomo.back_project(weight, flat_points, dwaterfall.flatten(),\n", + " nparts, nprofs)\n", + " # non-linearity of ART:\n", + " weight = weight.clip(0.0)\n", + "\n", + " print(f'Iteration: {i:3d}, discrepancy: {diff[-1]:3e}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "58290fc8-2bde-4080-a65a-f9f49e503e28", + "metadata": {}, + "outputs": [], + "source": [ + "# Finding last discrepancy...\n", + "rec_wf = tmo.libtomo.project(rec_wf, flat_points, weight, nparts, nprofs, nbins)\n", + "rec_wf /= np.sum(rec_wf, axis=1)[:, None]\n", + "dwaterfall = waterfall - rec_wf\n", + "diff.append(np.sqrt(np.sum(dwaterfall**2) / (nbins * nprofs)))\n", + "\n", + "print(f'Iteration: {i + 1:3d}, discrepancy: {diff[-1]:3E}')" + ] + }, + { + "cell_type": "markdown", + "id": "84e18024-dcd1-48b0-9039-e8c4716f8890", + "metadata": {}, + "source": [ + "Reconstruction finished, return the reconstructed phase-space distribution:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c91e185-a3c0-4a90-896b-46b34152c4d2", + "metadata": {}, + "outputs": [], + "source": [ + "phasespace = tomoout.create_phase_space_image(\n", + " phip, dEp, weight, nbins, reconstr_idx)" + ] + }, + { + "cell_type": "markdown", + "id": "81dfb247-7596-4cf5-a1c2-318112b824b8", + "metadata": {}, + "source": [ + "<h4>Results from Tomography</h4>\n", + "\n", + "The profile of the reconstructed distribution compared to the input profile:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "679c4c9a-95fa-4a51-b304-8f4117b8f391", + "metadata": {}, + "outputs": [], + "source": [ + "prof_rec = np.sum(phasespace, axis=1)\n", + "z_rec = (0.5 + np.arange(-len(prof_rec)/2, len(prof_rec)/2)) * (-dtbin * beta * c)\n", + "\n", + "plt.plot(z_rec, prof_rec, label='reconstructed')\n", + "plt.plot(z_rec, waterfall[reconstr_idx], label='measured')\n", + "\n", + "plt.xlabel('$z$ [m]')\n", + "plt.ylabel('Histogram')\n", + "plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1));" + ] + }, + { + "cell_type": "markdown", + "id": "442e751f-8916-46b4-973b-deb17017ac04", + "metadata": {}, + "source": [ + "The reconstructed rms bunch length is: $\\sigma_z = \\sqrt{\\cfrac{\\sum_i p(z_i) \\cdot (z_i - \\langle z \\rangle)^2}{\\sum_i p(z_i)}}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e82a46e4-5512-445d-88d1-0cb8b01d8412", + "metadata": {}, + "outputs": [], + "source": [ + "np.sqrt(np.trapz(prof_rec * z_rec**2, z_rec) / np.trapz(prof_rec, z_rec))" + ] + }, + { + "cell_type": "markdown", + "id": "d161c8bb-429a-42a1-8ab8-57e417a06c08", + "metadata": {}, + "source": [ + "The original macro-particle distribution was generated with:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6113c2eb-5ca4-479d-8077-b37fa2de7c0a", + "metadata": {}, + "outputs": [], + "source": [ + "sigma_z" + ] + }, + { + "cell_type": "markdown", + "id": "ffb08752-e4e1-4c94-9aee-4b1610f3c820", + "metadata": {}, + "source": [ + "The momentum distribution of the reconstructed distribution:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88bc333f-df5b-493b-8298-2cc589412960", + "metadata": {}, + "outputs": [], + "source": [ + "Etot = Ekin + m_p/e * c**2 # in eV" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b42e20a-2ccc-48d6-a8b2-bf0dbe2da628", + "metadata": {}, + "outputs": [], + "source": [ + "profdp_rec = np.sum(phasespace, axis=0)\n", + "dp_rec = (0.5 + np.arange(-len(profdp_rec)/2, len(profdp_rec)/2)) * (tracker.particles.dEbin / (Etot) * beta**2)\n", + "\n", + "plt.plot(dp_rec, profdp_rec)\n", + "plt.xlabel('$\\delta$')\n", + "plt.ylabel('Reconstructed histogram');" + ] + }, + { + "cell_type": "markdown", + "id": "83f10c9c-f244-4e63-87d8-ff86db5f6145", + "metadata": {}, + "source": [ + "The reconstructed rms momentum deviation is: $\\sigma_\\delta = \\sqrt{\\cfrac{\\sum_i p(\\delta_i) \\cdot (\\delta_i - \\langle \\delta \\rangle)^2}{\\sum_i p(\\delta_i)}}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d9bc160-c8e6-481d-9006-3cffa6cf1533", + "metadata": {}, + "outputs": [], + "source": [ + "np.sqrt(np.trapz(profdp_rec * dp_rec**2, dp_rec) / np.trapz(profdp_rec, dp_rec))" + ] + }, + { + "cell_type": "markdown", + "id": "edbaf367-4b32-4fb5-8c60-07a6a3e8c153", + "metadata": {}, + "source": [ + "Here is the full reconstructed phase-space distribution:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f6f8ce0-1b2a-4091-ba3e-310dc6306269", + "metadata": {}, + "outputs": [], + "source": [ + "plot_tomo(phasespace, z_rec, dp_rec, rfb);" + ] + }, + { + "cell_type": "markdown", + "id": "da897b01-ba2b-490f-9327-e035356cb45e", + "metadata": {}, + "source": [ + "$\\implies$ Does this fit the initial macro-particle distribution?\n", + "\n", + "$\\implies$ Re-run the tomographic reconstruction for the last profile (see settings part above, start from section B.)! Does it match the final macro-particle distribution?" + ] + }, + { + "cell_type": "markdown", + "id": "3e9f525c-fd4e-407a-b279-b74d0df58f0f", + "metadata": {}, + "source": [ + "The final macro-particle phase-space distribution after the simulation:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9bfba1ed-218b-4de8-9e46-64de032b6a50", + "metadata": {}, + "outputs": [], + "source": [ + "plot_mp(z, dp, rfb);" + ] + } + ], + "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.9.20" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} -- GitLab