Skip to content
Snippets Groups Projects
Unverified Commit ea77c599 authored by tegenolf's avatar tegenolf Committed by GitHub
Browse files

Upload notebook for week 9

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