general imports¶
from __future__ import division, print_function
import numpy as np
np.random.seed(42)
from scipy.constants import m_p, c, e
import matplotlib.pyplot as plt
%matplotlib inline
# sets the PyHEADTAIL directory etc.
try:
from settings import *
except:
pass
PySussix imports¶
Sussix is a tune analysis tool: https://cds.cern.ch/record/702438/files/sl-note-98-017.pdf . If you do not have PySussix, its python interface, find it under https://github.com/PyCOMPLETE/PySussix .
try:
from PySussix import Sussix
except ImportError:
print ('ERROR: This interactive test needs the PySussix module. Trying PySUSSIX')
from PySUSSIX import Sussix
print('PySUSSIX found')
PyHEADTAIL imports¶
from PyHEADTAIL.trackers.transverse_tracking import TransverseMap
from PyHEADTAIL.trackers.detuners import Chromaticity, AmplitudeDetuning
import PyHEADTAIL.particles.generators as generators
Setting up the machine and functions¶
# Basic parameters.
n_turns = 1024
n_segments = 1
n_macroparticles = 500
Q_x = 64.28
Q_y = 59.31
Q_s = 0.0020443
C = 26658.883
R = C / (2.*np.pi)
alpha_x_inj = 0.
alpha_y_inj = 0.
beta_x_inj = 66.0064
beta_y_inj = 71.5376
alpha_0 = 0.0003225
# Parameters for transverse map.
s = np.arange(0, n_segments + 1) * C / n_segments
alpha_x = alpha_x_inj * np.ones(n_segments)
beta_x = beta_x_inj * np.ones(n_segments)
D_x = np.zeros(n_segments)
alpha_y = alpha_y_inj * np.ones(n_segments)
beta_y = beta_y_inj * np.ones(n_segments)
D_y = np.zeros(n_segments)
def calc_sussix_spec(x, xp, y, yp, p_idx, turn, window_width, q_x, q_y, n_lines=10):
# Initialise Sussix object
SX = Sussix()
SX.sussix_inp(nt1=1, nt2=window_width, idam=2, ir=0, tunex=q_x, tuney=q_y)
tunes_x = np.empty(n_lines)
tunes_y = np.empty(n_lines)
SX.sussix(x[p_idx,turn:turn+window_width], xp[p_idx,turn:turn+window_width],
y[p_idx,turn:turn+window_width], yp[p_idx,turn:turn+window_width],
x[p_idx,turn:turn+window_width], xp[p_idx,turn:turn+window_width]) # this line is not used by sussix!
return SX.ox[:n_lines], SX.oy[:n_lines]
def track_n_save(bunch, map_):
n_particles = bunch.macroparticlenumber
x_i = np.empty((n_particles, n_turns))
xp_i = np.empty((n_particles, n_turns))
y_i = np.empty((n_particles, n_turns))
yp_i = np.empty((n_particles, n_turns))
for i in xrange(n_turns):
x_i[:,i] = bunch.x[:]
xp_i[:,i] = bunch.xp[:]
y_i[:,i] = bunch.y[:]
yp_i[:,i] = bunch.yp[:]
for m_ in map_:
m_.track(bunch)
return x_i, xp_i, y_i, yp_i
def analyse_n_plot(bunch, x_i, xp_i, y_i, yp_i):
ox = np.empty(bunch.macroparticlenumber)
oy = np.empty(bunch.macroparticlenumber)
print ('analysing particle spectra ... this may take some time.')
for p_idx in range(bunch.macroparticlenumber):
ox[p_idx], oy[p_idx] = calc_sussix_spec(x_i, xp_i, y_i, yp_i, p_idx,
turn=0, window_width=512, q_x=Q_x%1, q_y=Q_y%1, n_lines=1)
if (p_idx)%100 == 0:
print ('particle', p_idx)
fig = plt.figure(figsize=(20,20))
ax1 = fig.add_subplot(311)
ax2 = fig.add_subplot(312, sharex=ax1)
ax3 = fig.add_subplot(313)
ax1.scatter(ox, oy)
ax1.set_ylabel(r'$Q_y$', fontsize=20)
ax1.set_xlabel(r'$Q_x$', fontsize=20)
ax2.hist(ox, bins=50, color='blue')
ax2.set_xlabel(r'$Q_x$', fontsize=20)
ax3.hist(oy, bins=50, color='red')
ax3.set_xlabel(r'$Q_y$', fontsize=20)
print ('std dev. Qx', np.std(ox))
print ('std dev. Qy', np.std(oy))
plt.show()
os.remove('sussix.inp')
def generate_bunch(n_macroparticles, alpha_x, alpha_y, beta_x, beta_y, alpha_0, Q_s, R):
intensity = 1.05e11
sigma_z = 0.059958
gamma = 3730.26
eta = alpha_0 - 1. / gamma**2
gamma_t = 1. / np.sqrt(alpha_0)
p0 = np.sqrt(gamma**2 - 1) * m_p * c
beta_z = eta * R / Q_s
epsn_x = 3.75e-6 # [m rad]
epsn_y = 3.75e-6 # [m rad]
epsn_z = 4 * np.pi * sigma_z**2 * p0 / (beta_z * e)
bunch = generators.generate_Gaussian6DTwiss(
macroparticlenumber=n_macroparticles, intensity=intensity, charge=e,
gamma=gamma, mass=m_p, circumference=C,
alpha_x=alpha_x, beta_x=beta_x, epsn_x=epsn_x,
alpha_y=alpha_y, beta_y=beta_y, epsn_y=epsn_y,
beta_z=beta_z, epsn_z=epsn_z)
return bunch
Let's go¶
Amplitude detuning¶
Amplitude detuning (such as from octupole magnets) can be modelled with the AmplitudeDetuning
detuner class. Check out this LHC example with currents of $I=400~A$ where we expect $\approx 5\times 10^{-4}$ detuning in both transverse planes:
bunch = generate_bunch(
n_macroparticles, alpha_x_inj, alpha_y_inj, beta_x_inj, beta_y_inj,
alpha_0, Q_s, R)
ampl_det = AmplitudeDetuning.from_octupole_currents_LHC(i_focusing=400, i_defocusing=-400)
trans_map = TransverseMap(
s, alpha_x, beta_x, D_x, alpha_y, beta_y, D_y, Q_x, Q_y, [ampl_det])
map_ = [m for m in trans_map]
x_i, xp_i, y_i, yp_i = track_n_save(bunch, map_)
analyse_n_plot(bunch, x_i, xp_i, y_i, yp_i)
Chromatic detuning¶
Chromaticity can be modelled with the Chromaticity
detuner class. First-order chromaticity for values of $Q'_x=6$ and $Q'_y=3$ with LHC tunes yields:
bunch = generate_bunch(
n_macroparticles, alpha_x_inj, alpha_y_inj, beta_x_inj, beta_y_inj,
alpha_0, Q_s, R)
chroma = Chromaticity(Qp_x=[6], Qp_y=[3]) # note the lists!
trans_map = TransverseMap(
s, alpha_x, beta_x, D_x, alpha_y, beta_y, D_y, Q_x, Q_y, [chroma])
map_ = [m for m in trans_map]
x_i, xp_i, y_i, yp_i = track_n_save(bunch, map_)
analyse_n_plot(bunch, x_i, xp_i, y_i, yp_i)
Higher-order Chromatic Detuning¶
The Chromaticity
detuner class supports higher orders in chromatic detuning. Check out this LHC example with horizontal chromaticities of $Q'_x=6, Q''_x = 4\times10^4$ and vertical chromaticities of $Q'_y=3, Q''_y=0, Q'''_y=2\times10^8$:
bunch = generate_bunch(
n_macroparticles, alpha_x_inj, alpha_y_inj, beta_x_inj, beta_y_inj,
alpha_0, Q_s, R)
chroma = Chromaticity(Qp_x=[6., 4e4], Qp_y=[3., 0., 2e8])
# here the lists are extended beyond first order! ^^^
trans_map = TransverseMap(
s, alpha_x, beta_x, D_x, alpha_y, beta_y, D_y, Q_x, Q_y, [chroma])
map_ = [m for m in trans_map]
x_i, xp_i, y_i, yp_i = track_n_save(bunch, map_)
analyse_n_plot(bunch, x_i, xp_i, y_i, yp_i)
The corresponding detuning curves versus $\delta = \frac{p-p_0}{p_0}$ show a parabola in the horizontal and a cubic polynomial in the vertical:
fig = plt.figure(figsize=(12,8))
chroma = Chromaticity(Qp_x=[6., 4e4], Qp_y=[3., 0., 2e8])
chroma.generate_segment_detuner(dmu_x=1, dmu_y=1)
dp = np.linspace(-5e-4, 5e-4, 500)
dQx = chroma.segment_detuners[0].calc_detuning_x(dp)
dQy = chroma.segment_detuners[0].calc_detuning_y(dp)
plt.plot(dp, dQx, label='Detuning x')
plt.plot(dp, dQy, label='Detuning y')
plt.grid('on')
plt.xlabel('dp')
plt.ylabel('Detuning')
plt.legend(loc='upper left')
plt.show()