This interactive test file tests the inference of beam parameters (optics) via the beam coordinates as implemented in the Particles class. One can estimate e.g. the beta function of the machine at the current location from the macro-particle distribution.
The following tests are available:
The first and last slice should show a bigger error than the other slices in the inference of epsn_x. The machine parameters are the ones from the SPS. The following gammas are plotted: [1,20] (for small gammas), [17, 19] (around gamma-transition at 18), [1, 1e4] for the big picture
from __future__ import division, print_function
import numpy as np
np.random.seed(111)
import math as m
from scipy.constants import m_p, c, e
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
%matplotlib inline
# sets the PyHEADTAIL directory etc.
try:
from settings import *
except:
pass
from PyHEADTAIL.trackers.transverse_tracking import TransverseMap
from PyHEADTAIL.trackers.longitudinal_tracking import LinearMap
import PyHEADTAIL.particles.generators as generators
import PyHEADTAIL.cobra_functions.stats as st
from PyHEADTAIL.particles.slicing import UniformBinSlicer
from PyHEADTAIL.general.printers import SilentPrinter
PyHEADTAIL v1.10.5.268
def make_params_plotting_vs_gamma(
ax=[-1.4,-1], ay=[-1.1, 1.1], bx=[35, 81], by=[4, 51],
dx=[-3, 20], dy=[-4, 1000], N=50000):
# global parameters
nsamples = 10
Q_x = 20.13
Q_y = 20.18
Q_s = 0.017
epsx = [3.75e-6, 4e-5]
epsy = [3.75e-6, 3.75e-6]
gx = (1. + np.array(ax)*np.array(ax))/np.array(bx)
C = 6911.
R = C / (2.*np.pi)
alpha_0 = [0.00308]
# for the slicer
n_slices = 5
def _run_for_range(gammas):
"""Run the test for the specified gamma range."""
def _plot_error(values_dp, gammas, title='param', logscale=False, param=ax):
"""Plot the computed values using dp, dE in a nice way."""
plt.figure()
ax1, = plt.plot(gammas, values_dp[:,0], label= title + ' ' + str(param[0]))
ax2, = plt.plot(gammas, values_dp[:,1], label= title + ' ' + str(param[1]))
plt.xlabel('gamma')
plt.ylabel('relative <L1 error> of ' + str(nsamples) + ' samples, N=' + str(N))
plt.title(title)
plt.legend(handles=[ax1, ax2])
if logscale:
plt.yscale('log')
plt.show()
def _plot_error_per_slice(values, gammas, title='param', logscale=False, param=epsx):
slice_n = values.shape[2]
for pidx in range(len(epsx)):
plt.figure()
plt.plot(gammas, values[:,pidx,:])
plt.legend(['slice ' + str(i) for i in range(slice_n)])
plt.xlabel('gamma')
plt.ylabel('relative <L1 error> of ' + str(nsamples) + ' samples, N=' + str(N))
plt.title(title + ' = ' + str(epsx[pidx]))
plt.show()
ax_L1 = np.zeros(shape=(len(gammas),2))
bx_L1 = np.zeros(shape=(len(gammas),2))
dx_L1 = np.zeros(shape=(len(gammas),2))
epsx_L1 = np.zeros(shape=(len(gammas),2))
epsx_L1_per_slice = np.zeros(shape=(len(gammas),2, n_slices))
gx_L1 = np.zeros(shape=(len(gammas),2))
for i, gamma in enumerate(gammas):
for j in range(len(ax)):
for s in range(nsamples):
b = generate_bunch(N, ax[j], ay[j], bx[j], by[j],
LinearMap(alpha_0, C, Q_s), dx[j], dy[j], gamma=gamma,
epsn_x=epsx[j], epsn_y=epsy[j])
b._warningprinter=SilentPrinter() # to supress the small gamma warnings
# prepare the slicer
slicer = UniformBinSlicer(n_slices, z_cuts=(np.min(b.z)-0.1, np.max(b.z)+0.1))
slice_set = slicer.slice(b)
ax_L1[i,j] += abs(ax[j] - b.alpha_Twiss_x())/abs(ax[j])
bx_L1[i,j] += abs(bx[j] - b.beta_Twiss_x())/abs(bx[j])
dx_L1[i,j] += abs(dx[j] - b.dispersion_x())/abs(dx[j])
epsx_L1[i,j] += abs(epsx[j] - b.epsn_x())/abs(epsx[j])
gx_L1[i,j] += abs(gx[j] - b.gamma_Twiss_x())/abs(gx[j])
epsx_per_slice = slicer._epsn_x(slice_set, b)
for sliceidx in xrange(n_slices):
epsx_L1_per_slice[i,j,sliceidx] += abs(
epsx[j] - epsx_per_slice[sliceidx]) / abs(epsx[j])
ax_L1[i,j] /= nsamples
dx_L1[i,j] /= nsamples
bx_L1[i,j] /= nsamples
epsx_L1[i,j] /= nsamples
gx_L1[i,j] /= nsamples
for sliceidx in range(n_slices):
epsx_L1_per_slice[i,j,sliceidx] /= nsamples
_plot_error(ax_L1, gammas, title='alpha_x', param=ax)
_plot_error(bx_L1, gammas, title='beta_x' , param=bx)
_plot_error(gx_L1, gammas, title='gamma_x', param=gx)
_plot_error(dx_L1, gammas, title='disp_x', param=dx, logscale=True)
_plot_error(epsx_L1, gammas, title='epsn_x',param=epsx)
# plot the emittance per slice
_plot_error_per_slice(epsx_L1_per_slice, gammas, title='epsx', param=epsx)
return _run_for_range
run_gamma_range = make_params_plotting_vs_gamma()
def generate_bunch(n_macroparticles, alpha_x, alpha_y, beta_x, beta_y,
linear_map, dispx, dispy, epsn_x=3.75e-6,
epsn_y=3.75e-6, gamma=3730.27):
# generate a alpha=0, beta=1 bunch first
# then transform...
intensity = 1.05e11
sigma_z = 0.23#0.59958
gamma_t = 1. / np.sqrt(linear_map.alpha_array[0])
p0 = np.sqrt(gamma**2 - 1) * m_p * c
beta_z = np.abs((linear_map.eta(dp=0, gamma=gamma) * linear_map.circumference /
(2 * np.pi * linear_map.Q_s)))
epsn_z = 4 * np.pi * sigma_z**2 * p0 / (beta_z * e)
#print ('epsn_z: ' + str(epsn_z))
bunch = generators.generate_Gaussian6DTwiss(
macroparticlenumber=n_macroparticles, intensity=intensity, charge=e,
gamma=gamma, mass=m_p, circumference=linear_map.circumference,
alpha_x=0., beta_x=1., epsn_x=epsn_x,
alpha_y=0., beta_y=1., epsn_y=epsn_y,
beta_z=beta_z, epsn_z=epsn_z)
# Scale to correct beta and alpha
xx = bunch.x.copy()
yy = bunch.y.copy()
bunch.x *= np.sqrt(beta_x)
bunch.xp = -alpha_x/np.sqrt(beta_x) * xx + 1./np.sqrt(beta_x) * bunch.xp
bunch.y *= np.sqrt(beta_y)
bunch.yp = -alpha_y/np.sqrt(beta_y) * yy + 1./np.sqrt(beta_y) * bunch.yp
bunch.x += dispx * bunch.dp
bunch.y += dispy * bunch.dp
return bunch
def test_map(n_macroparticles, linear_map, gamma=3730.27):
'''Test the beam optics extraction functions when tracking the beam
through a transverse and longitudinal map (including dispersion!).
'''
ax = np.array([-1.2, 2., 3.])
ay = np.array([0.3, 2., 3.])
bx = np.array([82., 90., 100.])
by = np.array([45., 30, 80.])
gx = (ax*ax + 1.) / bx
gy = (ay*ay + 1.) / by
dx = np.array([1, 2., 5])
dy = np.array([0.2, -4.1, 2.09])
segments = np.linspace(0., C, num=4)
# generate a map with the given parameters
trans_map = TransverseMap(segments, ax, bx, dx, ay, by, dy, Q_x, Q_y)
map_ = [m for m in trans_map] + [LinearMap(alpha_0, C, Q_s, D_x=1, D_y=0.2)]
# generate a bunch matched to the alpha, beta and d of the first segment
b = generate_bunch(n_macroparticles, ax[0], ay[0], bx[0], by[0],
linear_map, dx[0], dy[0], gamma=gamma)
def print_parameters(b, it):
it = it % 3 #access 1st segment's parameters after last segment
print('alphax: ' + str(ax[it]) + ' <---> ' + str(b.alpha_Twiss_x()))
print('alphay: ' + str(ay[it]) + ' <---> ' + str(b.alpha_Twiss_y()))
print('betax: ' + str(bx[it]) + ' <---> ' + str(b.beta_Twiss_x()))
print('betay: ' + str(by[it]) + ' <---> ' + str(b.beta_Twiss_y()))
print('gammax: ' + str(gx[it]) + ' <---> ' + str(b.gamma_Twiss_x()))
print('gammay: ' + str(gy[it]) + ' <---> ' + str(b.gamma_Twiss_y()))
print('dispx: ' + str(dx[it]) + ' <---> ' + str(b.dispersion_x()))
print('dispy: ' + str(dy[it]) + ' <---> ' + str(b.dispersion_y()))
print('epsnx: ' + str(3.75e-6) + ' <---> ' + str(b.epsn_x()))
print('epsny: ' + str(3.75e-6) + ' <---> ' + str(b.epsn_y()))
# print('epsnz: ' + 'TODO' + ' <---> ' + str(b.epsn_z()))
print('\n')
def plot_params(data, title, yaxis):
data = [1e6*d for d in data]
plt.plot(data)
plt.ylim(min(data), max(data))
plt.xlabel('# turns')
plt.ylabel('value * 1e6')
plt.title(title)
plt.show()
print_parameters(b, 0) # matched to the first segment
#track and check parameters!
nturns = 100
e_x = np.empty(nturns, dtype=np.float64)
e_y = np.empty(nturns, dtype=np.float64)
for n in range(nturns):
for it,m in enumerate(map_):
m.track(b)
#print_parameters(b, it+1)
plt.show()
e_x[n] = b.epsn_x()
e_y[n] = b.epsn_y()
plot_params(data=e_x, title='epsx', yaxis='value')
plot_params(data=e_x, title='epsy', yaxis='value')
#check after the last segment
def plot_phase_space(beam, direction):
plt.figure()
u = getattr(beam, direction[0])
up = getattr(beam, direction[1])
plt.scatter(u, up)
plt.xlim([min(u), max(u)])
plt.ylim([min(up), max(up)])
plt.xlabel(direction[0])
plt.ylabel(direction[1])
plt.show()
def compute_statistical_beam_optics(alphax, alphay, betax, betay,
long_map, dispx, dispy, gamma=3000):
'''Generate a bunch with the given optics and compare them
to the statistically computed parameters.
'''
b = generate_bunch(n_macroparticles, alphax, alphay, betax, betay,
long_map, dispx, dispy, gamma=gamma)
gammax = (1. + alphax**2)/betax
gammay = (1. + alphay**2)/betay
print ('Compare the nominal and the statistically computed quantities from the beam')
print ('nominal: <---> statistical:')
print ('alphax: ' + str(alphax) + ' <---> ' + str(b.alpha_Twiss_x()))
print ('alphay: ' + str(alphay) + ' <---> ' + str(b.alpha_Twiss_y()))
print ('betax: ' + str(betax) + ' <---> ' + str(b.beta_Twiss_x()))
print ('betay: ' + str(betay) + ' <---> ' + str(b.beta_Twiss_y()))
print ('gammax: ' + str(gammax) + ' <---> ' + str(b.gamma_Twiss_x()))
print ('gammay: ' + str(gammay) + ' <---> ' + str(b.gamma_Twiss_y()))
print ('dispx: ' + str(dispx) + ' <---> ' + str(b.dispersion_x()))
print ('dispy: ' + str(dispy) + ' <---> ' + str(b.dispersion_y()))
def plot_histogram(measurements, exact_value, title="Dummy"):
plt.hist(measurements)
plt.axvline(measurements.mean(), color='r', linestyle='dashed', linewidth=2)
plt.axvline(exact_value, color='g', linewidth=2)
plt.title(title)
plt.xlabel('Value')
plt.ylabel('# observations')
m = mpatches.Patch(color='red', label='mean')
ex = mpatches.Patch(color='green', label='input')
plt.legend(handles=[m, ex])
plt.show()
def compute_distribution_beam_optics(alphax, alphay, betax, betay, long_map, dispx, dispy, gamma=3000):
nsamples = 100
stat_betax = np.zeros(nsamples)
stat_betay = np.zeros(nsamples)
stat_alphax = np.zeros(nsamples)
stat_alphay = np.zeros(nsamples)
stat_dispx = np.zeros(nsamples)
stat_dispy = np.zeros(nsamples)
for i in range(nsamples):
b = generate_bunch(n_macroparticles, alphax, alphay, betax, betay, long_map, dispx, dispy,gamma=gamma)
stat_betax[i] = b.beta_Twiss_x()
stat_betay[i] = b.beta_Twiss_y()
stat_dispx[i] = b.dispersion_x()
stat_dispy[i] = b.dispersion_y()
stat_alphax[i] = b.alpha_Twiss_x()
stat_alphay[i] = b.alpha_Twiss_y()
plot_histogram(stat_betax, betax, 'beta_x')
plot_histogram(stat_betay, betay, 'beta_y')
plot_histogram(stat_alphax, alphax, 'alpha_x')
plot_histogram(stat_alphay, alphay, 'alpha_y')
plot_histogram(stat_dispx, dispx, 'disp_x')
plot_histogram(stat_dispy, dispy, 'disp_y')
plt.show()
Infer the lattice functions statistically from the macro-particle distribution and compare against the real values:
# Basic simulation parameters
n_macroparticles = 100000
Q_x = 20.13
Q_y = 20.18
Q_s = 0.017
C = 6911.
R = C / (2.*np.pi)
alpha_0 = [0.00308]
long_map = LinearMap(alpha_0, C, Q_s)
alphax = 0.5
alphay = -1.7
betax = 90
betay = 40
dispx = 4.3
dispy = -3.
gamma = 1000
beam = generate_bunch(n_macroparticles, alphax, alphay, betax, betay, long_map, dispx, dispy, gamma=gamma)
print ('total energy [J]: ' + str(beam.gamma*beam.mass*c*c) + '\n')
compute_statistical_beam_optics(alphax, alphay, betax, betay, long_map, dispx, dispy, gamma)
plot_phase_space(beam, ['x', 'xp'])
plot_phase_space(beam, ['y', 'yp'])
plot_phase_space(beam, ['z', 'dp'])
total energy [J]: 1.5032775929e-07 Compare the nominal and the statistically computed quantities from the beam nominal: <---> statistical: alphax: 0.5 <---> 0.504882093933 alphay: -1.7 <---> -1.70320922969 betax: 90 <---> 89.8278809489 betay: 40 <---> 39.9530974011 gammax: 0.0138888888889 <---> 0.0139701161323 gammay: 0.09725 <---> 0.0976375283482 dispx: 4.3 <---> 4.29792798153 dispy: -3.0 <---> -3.00188080776
Repeat the comparison several times with different random generator seeds for the macro-particle distribution. Plot histograms for each of the inferred optics parameters:
compute_distribution_beam_optics(alphax, alphay, betax, betay, long_map, dispx, dispy, gamma)
The following plots show the deviation (error) from the inferred optics functions (based on the macro-particle distribution with $N=50000$) compared to the defined values in the lattice. For higher macroparticlenumber
s the error would become smaller. This becomes obvious in the plots showing the error over the slices: the central slice 2
of the longitudinal thermal ("Gaussian") distribution hosts most macro-particles particles -- its error is relatively small compared to the head and tail slices 0
and 4
which consist of much less macro-particles.
gammas = np.arange(1.5,19,1)
# gammas = np.arange(17,19,0.1)
# gammas = np.arange(2,10000,100)
run_gamma_range(gammas)
Running a simulation over a few turns and checking the distribution inferred optics parameters
test_map(50000, long_map, 4000)
Non-zero dispersion; ensure the beam has been "blown-up" accordingly upon creation! alphax: -1.2 <---> -1.20269283555 alphay: 0.3 <---> 0.297017573584 betax: 82.0 <---> 81.8275895628 betay: 45.0 <---> 44.8106137866 gammax: 0.029756097561 <---> 0.0298978629305 gammay: 0.0242222222222 <---> 0.0242848590336 dispx: 1.0 <---> 0.999584001566 dispy: 0.2 <---> 0.20048512313 epsnx: 3.75e-06 <---> 3.74801803588e-06 epsny: 3.75e-06 <---> 3.75350957245e-06