general imports

In [1]:
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
In [2]:
# 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 .

In [3]:
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

In [4]:
from PyHEADTAIL.trackers.transverse_tracking import TransverseMap
from PyHEADTAIL.trackers.detuners import Chromaticity, AmplitudeDetuning
import PyHEADTAIL.particles.generators as generators
PyHEADTAIL v1.10.5.269


Setting up the machine and functions

In [5]:
# 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)
In [6]:
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:

In [7]:
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)
analysing particle spectra ... this may take some time.
particle 0
particle 100
particle 200
particle 300
particle 400
std dev. Qx 0.000493943454511
std dev. Qy 0.000504673628496

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:

In [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], 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)
analysing particle spectra ... this may take some time.
particle 0
particle 100
particle 200
particle 300
particle 400
std dev. Qx 0.000561082708596
std dev. Qy 0.00028045654874

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$:

In [9]:
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)
analysing particle spectra ... this may take some time.
particle 0
particle 100
particle 200
particle 300
particle 400
std dev. Qx 0.000568932172386
std dev. Qy 0.000322448107092

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:

In [10]:
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()