Track
Tracking particles with Xsuite
The tracking of particles through a beam line can be simulated using the
xtrack.Line.track()
method of the xtrack.Line
class. This is
illustrated in the following example:
import numpy as np
import xobjects as xo
import xtrack as xt
## Generate a simple line
line = xt.Line(
elements=[xt.Drift(length=2.),
xt.Multipole(knl=[0, 0.5], ksl=[0,0]),
xt.Drift(length=1.),
xt.Multipole(knl=[0, -0.5], ksl=[0,0])],
element_names=['drift_0', 'quad_0', 'drift_1', 'quad_1'])
## Attach a reference particle to the line (optional)
## (defines the reference mass, charge and energy)
line.particle_ref = xt.Particles(p0c=6500e9, #eV
q0=1, mass0=xt.PROTON_MASS_EV)
## Choose a context
context = xo.ContextCpu() # For CPU
# context = xo.ContextCupy() # For CUDA GPUs
# context = xo.ContextPyopencl() # For OpenCL GPUs
## Transfer lattice on context and compile tracking code
line.build_tracker(_context=context)
## Build particle object on context
n_part = 20
particles = line.build_particles(
x=np.random.uniform(-1e-3, 1e-3, n_part),
px=np.random.uniform(-1e-5, 1e-5, n_part),
y=np.random.uniform(-2e-3, 2e-3, n_part),
py=np.random.uniform(-3e-5, 3e-5, n_part),
zeta=np.random.uniform(-1e-2, 1e-2, n_part),
delta=np.random.uniform(-1e-4, 1e-4, n_part))
# Reference mass, charge, energy are taken from the reference particle.
# Particles are allocated on the context chosen for the line.
## Track (no saving of turn-by-turn data)
n_turns = 100
line.track(particles, num_turns=n_turns)
particles.state # > 0 for particles still alive
particles.at_turn # turn number (for lost particles, it is the turn of loss)
particles.x # x position after tracking
particles.px # x momentum after tracking
# etc...
## Track (saving turn-by-turn data)
particles = line.build_particles( # fresh particles
x=np.random.uniform(-1e-3, 1e-3, n_part),
px=np.random.uniform(-1e-5, 1e-5, n_part),
y=np.random.uniform(-2e-3, 2e-3, n_part),
py=np.random.uniform(-3e-5, 3e-5, n_part),
zeta=np.random.uniform(-1e-2, 1e-2, n_part),
delta=np.random.uniform(-1e-4, 1e-4, n_part))
n_turns = 100
line.track(particles, num_turns=n_turns,
turn_by_turn_monitor=True)
## Turn-by-turn data is available at:
line.record_last_track.x
line.record_last_track.px
# etc...
# Complete source: xtrack/examples/tracker/000_track.py
Dynamic aperture
The following example illustrate how to use Xsuite to study the Dynamic Aperture of a ring.
import numpy as np
import xtrack as xt
import xpart as xp
import xobjects as xo
line = xt.Line.from_json(
'../../test_data/hllhc14_no_errors_with_coupling_knobs/line_b1.json')
line.build_tracker()
# Switch on LHC octupole circuits to have a smaller dynamic aperture
for arc in ['12', '23', '34', '45', '56', '67', '78', '81']:
line.vars[f'kod.a{arc}b1'] = 2.0
line.vars[f'kof.a{arc}b1'] = 2.0
# Generate normalized particle coordinates on a polar grid
n_r = 50
n_theta = 60
x_normalized, y_normalized, r_xy, theta_xy = xp.generate_2D_polar_grid(
r_range=(0, 40.), # beam sigmas
theta_range=(0, np.pi/2),
nr=n_r, ntheta=n_theta)
# Set initial momentum deviation
delta_init = 0 # In case off-momentum DA is needed
# Match particles to the machine optics and orbit
particles = line.build_particles(
x_norm=x_normalized, px_norm=0,
y_norm=y_normalized, py_norm=0,
nemitt_x=3e-6, nemitt_y=3e-6, # normalized emittances
delta=delta_init)
# # Optional: activate multi-core CPU parallelization
# line.discard_tracker()
# line.build_tracker(_context=xo.ContextCpu(omp_num_threads='auto'))
# Track
line.track(particles, num_turns=200, time=True, with_progress=5)
print(f'Tracked in {line.time_last_track} seconds')
# Sort particles to get the initial order
# (during tracking lost particles are moved to the end)
particles.sort(interleave_lost_particles=True)
# Plot result using scatter or pcolormesh
import matplotlib.pyplot as plt
plt.close('all')
plt.figure(1)
plt.scatter(x_normalized, y_normalized, c=particles.at_turn)
plt.xlabel(r'$A_x [\sigma]$')
plt.ylabel(r'$A_y [\sigma]$')
cb = plt.colorbar()
cb.set_label('Lost at turn')
plt.figure(2)
plt.pcolormesh(
x_normalized.reshape(n_r, n_theta), y_normalized.reshape(n_r, n_theta),
particles.at_turn.reshape(n_r, n_theta), shading='gouraud')
plt.xlabel(r'$A_x [\sigma]$')
plt.ylabel(r'$A_y [\sigma]$')
ax = plt.colorbar()
ax.set_label('Lost at turn')
plt.show()
# Complete source: xtrack/examples/dynamic_aperture/000_tracking_for_da.py
Monitors
See also: xtrack.ParticlesMonitor
, xtrack.LastTurnsMonitor
, xtrack.BeamPositionMonitor
, xtrack.BeamProfileMonitor
, xtrack.BeamSizeMonitor
.
The easy way
When starting a tracking simulation with the Xtrack Line object, the easiest
way of logging the coordinates of all particles for all turns is to enable the
default turn-by-turn monitor, as illustrated by the following example.
Note: this mode requires that particles.at_turn
is 0
for all particles
at the beginning of the simulation.
import json
import xtrack as xt
import xpart as xp
import xobjects as xo
context = xo.ContextCpu()
with open('../../test_data/hllhc15_noerrors_nobb/line_and_particle.json') as f:
dct = json.load(f)
line = xt.Line.from_dict(dct['line'])
line.particle_ref = xt.Particles.from_dict(dct['particle'])
line.build_tracker()
num_particles = 50
particles = xp.generate_matched_gaussian_bunch(line=line,
num_particles=num_particles,
nemitt_x=2.5e-6,
nemitt_y=2.5e-6,
sigma_z=9e-2)
num_turns = 30
line.track(particles, num_turns=num_turns,
turn_by_turn_monitor=True # <--
)
# line.record_last_track contains the measured data. For example,
# line.record_last_track.x contains the x coordinate for all particles
# and all turns, e.g. line.record_last_track.x[3, 5] for the particle
# having particle_id = 3 and for the turn number 5.
# Monitor objects can be dumped to a dictionary and loaded back
mon = line.record_last_track
dct = mon.to_dict()
mon2 = xt.ParticlesMonitor.from_dict(dct)
# Complete source: xtrack/examples/monitor/000_example_quick_monitor.py
Custom monitors
In order to customize the monitor’s behaviour,
a custom monitor object can be built and passed to the Line.track
function.
Particles coordinates can be recorded only in a selected range of turns
by specifying start_at_turn
and stop_at_turn
.
The monitoring can also be limited to a selected range of particles IDs,
by using the argument particle_id_range
of the ParticlesMonitor
class
to provide a tuple defining the range to be recorded. In that case the
num_particles
input of the monitor is omitted.
The example above is changed as follows:
...
monitor = xt.ParticlesMonitor(_context=context,
particle_id_range=(5, 42),
start_at_turn=5, # <-- first turn to monitor (including)
stop_at_turn=15, # <-- last turn to monitor (excluding)
)
line.track(particles, num_turns=num_turns,
turn_by_turn_monitor=monitor, # <-- pass the monitor here
)
Now, line.record_last_track.x[3, 5]
gives the x coordinates for the
particle 3 (which has the id 8) and the
recorded turn 5 (which is turn number 10)
The particle ids that are recorded can be inspected in line.record_last_track.particle_id
and the turn indeces in line.record_last_track.at_turn
.
Multi-frame particles monitor
The particles monitor can also record periodically spaced intervals of turns (frames)
This feature can be activated by providing the arguments n_repetitions
and
repetition_period
when creating the monitor.
In the following example, we record turns in range 5 to 10 (first frame),
range 25 to 30 (second frame) and range 45 to 50 (third frame).
Note that each frame consists of 5 turns since stop_at_turn
is excluding.
monitor = xt.ParticlesMonitor(_context=context,
num_particles=num_particles,
start_at_turn=5,
stop_at_turn=10,
n_repetitions=3, # <--
repetition_period=20, # <--
)
Now, the measured data are 3D array with the first index being the frame number.
For example, line.record_last_track.x[0, :, :]
contains the recorded
x position for the first frame (turns 5, 6, 7, 8 and 9) and
line.record_last_track.x[-1, :, 0]
refers to the last frame and the first turn within,
which is turn number turn 25.
As before, the turn numbers recorded can be inspected with line.record_last_track.at_turn
.
Particles monitor as beam elements
Particles monitors can be used as regular beam element to record the particle coordinates at specific locations in the beam line. For this purpose they can be inserted in the line, as illustrated in the following example.
import json
import xtrack as xt
import xpart as xp
import xobjects as xo
context = xo.ContextCpu()
with open('../../test_data/hllhc15_noerrors_nobb/line_and_particle.json') as f:
dct = json.load(f)
line = xt.Line.from_dict(dct['line'])
line.particle_ref = xt.Particles.from_dict(dct['particle'])
num_particles = 50
monitor_ip5 = xt.ParticlesMonitor(start_at_turn=5, stop_at_turn=15,
num_particles=num_particles)
monitor_ip8 = xt.ParticlesMonitor(start_at_turn=5, stop_at_turn=15,
num_particles=num_particles)
line.insert_element(index='ip5', element=monitor_ip5, name='mymon5')
line.insert_element(index='ip8', element=monitor_ip8, name='mymon8')
line.build_tracker()
particles = xp.generate_matched_gaussian_bunch(line=line,
num_particles=num_particles,
nemitt_x=2.5e-6,
nemitt_y=2.5e-6,
sigma_z=9e-2)
num_turns = 30
monitor = xt.ParticlesMonitor(_context=context,
start_at_turn=5, stop_at_turn=15,
num_particles=num_particles)
line.track(particles, num_turns=num_turns)
# monitor_ip5 contains the data recorded in before the element 'ip5', while
# monitor_ip8 contains the data recorded in before the element 'ip8'
# The element index at which the recording is made can be inspected in
# monitor_ip5.at_element.
# Complete source: xtrack/examples/monitor/003_monitors_as_beam_elements.py
As all Xtrack elements, the Particles Monitor has a track method and can be used in stand-alone mode as illustrated in the following example.
# line.track(particles, num_turns=num_turns)
for iturn in range(num_turns):
monitor.track(particles)
line.track(particles)
Last turns monitor
The xtrack.LastTurnsMonitor
records particle data in the last turns before respective particle loss
(or the end of tracking).
The idea is to use a rolling buffer instead of saving all the turns. This saves a lot of memory resources
when the interest lies only in the last few turns.
For each particle, the recorded data will cover up to n_last_turns*every_n_turns
turns before it is lost (or the tracking ends).
monitor = LastTurnsMonitor(
particle_id_range=(0, 5),
n_last_turns=5, # amount of turns to store
every_n_turns=3, # only consider turns which are a multiples of this
)
... # track
monitor.at_turn[:,-1] # turn number of each particle before it is lost (last turn alive)
monitor.x[3,-2] # x coordinate of particle 3 in one but last turn (-2)
The monitor provides the following data as 2D array of shape (num_particles, n_last_turns)
,
where the first index corresponds to the particle in particle_id_range
and the second index corresponds to the turn (or every_n_turns) before the respective particle is lost:
particle_id
, at_turn
, x
, px
, y
, py
, delta
, zeta
Beam position monitor
The xtrack.BeamPositionMonitor
records transverse beam positions,
i.e. it stores the x and y centroid positions of particles.
This can be useful for tune or beam-transfer-function diagnostics
as well as transverse schottky spectra.
The monitor allows for arbitrary sampling frequencies and can thus not only be used for bunch positions, but also coasting beam positions. Higher sampling frequencies give access to transverse beam oscillations at higher harmonics, which is especially useful for schottky diagnostics. Internally, the particle arrival time is used when determining the record index. For coasting beams this ensures, that the centroid is computed considering all particles which arrive at the monitor at the same time (as in a real-world measurement device), even if some particles might have made more or less turns than the synchronous particle due to a non-negligible momentum deviation.
where
\(f_{samp}\) is the sampling frequency,
\(f_{rev}\) is the revolution frequency,
\(n\) is the current turn number and \(n_0\) is the first turn recorded,
\(\zeta=(s-\beta_0\cdot c_0\cdot t)\) is the longitudinal zeta
coordinate of the particle,
\(\beta_0\) is the relativistic beta factor of the particle
and \(c_0\) is the speed of light.
For non-circular lines \(n\) is always zero and \(f_{rev}\) can be omitted.
Note that the index is rounded, i.e. the result array represents data of particles equally distributed around the reference particle, which is useful for bunched beams. For example, if the sampling frequency is twice the revolution frequency, the first item contains data from particles in the range zeta/circumference = -0.25 .. 0.25, the second item in the range 0.25 .. 0.75 and so on.
monitor = xt.BeamPositionMonitor(
#particle_id_range=(5, 42), # optional, defaults to all particles if not given
start_at_turn=5, stop_at_turn=10, # turn refers to the synchronous particle (at zeta=0)
frev=1e6, # revolution frequency (only for circular lines)
sampling_frequency=2e6, # sampling frequency
)
... # track
print(monitor.count) # waveform of number of particles (intensity)
print(monitor.x_mean) # waveform of horizontal centroid positions (alias monitor.x_cen)
print(monitor.y_mean) # waveform of vertical centroid positions (alias monitor.y_cen)
The result arrays can be understood as waveforms recorded at the specified sampling frequency. In the special case where sampling frequency was set to the same value as the revolution frequency, the indices are identical to the recorded turn numbers (of the synchronous particle).
Beam size monitor
The xtrack.BeamSizeMonitor
records transverse beam sizes,
i.e. it stores the standard deviation of the particles x and y positions.
Like the Beam position monitor also the beam size monitor is based on particle arrival time and an arbitrary sampling frequency.
monitor = xt.BeamSizeMonitor(
#particle_id_range=(5, 42), # optional, defaults to all particles if not given
start_at_turn=5, stop_at_turn=10, # turn refers to the synchronous particle (at zeta=0)
frev=1e6, # revolution frequency (only for circular lines)
sampling_frequency=2e6, # sampling frequency
)
... # track
print(monitor.count) # waveform of number of particles (intensity)
print(monitor.x_mean) # waveform of horizontal centroid positions
print(monitor.y_std) # waveform of vertical position standard deviation (i.e. beam size)
print(monitor.x_var) # waveform of horizontal position variances
Beam profile monitor
The xtrack.BeamProfileMonitor
records transverse beam profiles,
i.e. it stores the number of particles on a defined raster (like a histogram).
Like the Beam position monitor also the beam profile monitor is based on particle arrival time and an arbitrary sampling frequency.
monitor = xt.BeamProfileMonitor(
#particle_id_range=(5, 42), # optional, defaults to all particles if not given
start_at_turn=5, stop_at_turn=10, # turn refers to the synchronous particle (at zeta=0)
frev=1e6, # revolution frequency (only for circular lines)
sampling_frequency=2e6, # sampling frequency
n=100, # number of bins in the profile (can also specify nx and ny separately)
x_range=(-4,2), # save horizontal profile extending from -4 to 2
y_range=5, # shorthand for (-2.5, 2.5)
)
... # track
print(monitor.x_edges) # the bin edges
print(monitor.x_grid) # the bin midpoints
print(monitor.x_intensity) # the actual profile (particle count per bin)
The recorded profiles are 2D arrays of shape (sample_size, n)
where sample_size = round(( stop_at_turn - start_at_turn ) * sampling_frequency / frev)
.
I.e. monitor.x_intensity[0,:]
is the first recorded profile and monitor.x_intensity[-1,:]
the last.
Start/stop tracking at specific elements
It is possible to start and/or stop the tracking at specific elements of the beam line. This is illustrated in the following example:
import xtrack as xt
### Behavior:
# For num_turns = 1 (default):
# - line.track(particles, ele_start=3, ele_stop=5) tracks once from element 3 to
# element 5 e. particles.at_turn is not incremented as the particles never pass through
# the end of the beam line.
# - line.track(particles, ele_start=7, ele_stop=3) tracks once from element 5 to
# the end of the beamline and then from element 0 to element 3 excluded.
# particles.at_turn is incremented to one as the particles pass once
# the end of the beam line.
# When indicating num_turns = N with N > 1, additional (N-1) full turns are added to logic
# above. Therefore:
# - if ele_stop < ele_start: stops at element ele_stop when particles.at_turn = N - 1
# - if ele_stop >= ele_start: stops at element ele_stop when particles.at_turn = N
elements = [xt.Drift(length=1.) for _ in range(10)]
line=xt.Line(elements=elements)
line.build_tracker()
line.reset_s_at_end_turn = False
# Standard mode
p = xt.Particles(x=[1e-3, 2e-3, 3e-3], p0c=7e12)
line.track(p, num_turns=4)
p.at_turn # is [4, 4, 4]
# ele_start > 0
p = xt.Particles(x=[1e-3, 2e-3, 3e-3], p0c=7e12, s=4., at_element=2)
line.track(p, num_turns=3, ele_start=2)
p.at_turn # is [3, 3, 3]
# 0 <= ele_start < ele_stop
p = xt.Particles(x=[1e-3, 2e-3, 3e-3], p0c=7e12, s=5 * 2., at_element=5)
line.track(p, num_turns=4, ele_start=5, ele_stop=8)
p.at_turn # is [3, 3, 3]
# 0 <= ele_start < ele_stop -- num_turns = 1
p = xt.Particles(x=[1e-3, 2e-3, 3e-3], p0c=7e12, s=5 * 2., at_element=5)
line.track(p, num_turns=1, ele_start=5, ele_stop=2)
p.at_turn # is [1, 1, 1]
# 0 <= ele_start < ele_stop -- num_turns > 1
p = xt.Particles(x=[1e-3, 2e-3, 3e-3], p0c=7e12, s=5 * 2., at_element=5)
line.track(p, num_turns=4, ele_start=5, ele_stop=2)
p.at_turn # is [4, 4, 4]
# Use ele_start and num_turns
p = xt.Particles(x=[1e-3, 2e-3, 3e-3], p0c=7e12, s=5 * 2., at_element=5)
line.track(p, ele_start=5, num_elements=7)
p.at_turn # is [1, 1, 1]
# Complete source: xtrack/examples/tracker/001_tracker_start_stop.py
Backtracking
It is possible to track particles backwards through a beam line, provided that all elements included in the line support backtracking. The following example illustrates how backtrack for a full turn or between specified elements:
import xtrack as xt
line = xt.Line.from_json(
'../../test_data/hllhc15_noerrors_nobb/line_w_knobs_and_particle.json')
line.build_tracker()
p = xt.Particles(
p0c=7000e9, x=1e-4, px=1e-6, y=2e-4, py=3e-6, zeta=0.01, delta=1e-4)
p.get_table().show(digits=3)
# prints:
#
# particle_id s x px y py zeta delta chi ...
# 0 0 0.0001 1e-06 0.0002 3e-06 0.01 0.0001 1
# Track one turn
line.track(p)
p.get_table().show()
# prints:
#
# particle_id s x px y py zeta delta chi ...
# 0 0 0.000324 -1.07e-05 1.21e-05 -1.42e-06 0.00909 0.0001 1
# Track back one turn
line.track(p, backtrack=True)
# The particle is back to the initial coordinates
p.get_table().show(digits=3)
# prints:
#
# particle_id s x px y py zeta delta chi ...
# 0 0 0.0001 1e-06 0.0002 3e-06 0.01 0.0001 1
# It is also possible to backtrack with a specified start/stop elements
# Track three elements
line.track(p, ele_start=0, ele_stop=3)
p.get_table().cols['x px y py zeta delta at_element'].show()
# prints:
#
# particle_id x px y py zeta delta at_element
# 0 0.000121028 1e-06 0.000263084 3e-06 0.01 0.0001 3
# Track back three elements
line.track(p, ele_start=0, ele_stop=3, backtrack=True)
p.get_table().cols['x px y py zeta delta at_element'].show()
# prints:
#
# particle_id x px y py zeta delta at_element
# 0 0.0001 1e-06 0.0002 3e-06 0.01 0.0001 0
# Complete source: xtrack/examples/tracker/002_backtrack.py
Freeze longitudinal coordinates
In certain studies, it is convenient to track particles updating only the transverse coordinates, while keeping the longitudinal coordinates fixed (frozen). Xsuite offers the possibility to freeze the longitudinal coordinates within a single method or changing the state of the Line, as illustrated in the following sections.
Freezing longitudinal when calling methods
The Line.twiss
and Line.track
can work with frozen longitudinal
coordinates. This is done by setting the freeze_longitudinal
argument to
True
, as shown in the following example:
import json
import xtrack as xt
# import a line and add reference particle
line = xt.Line.from_json(
'../../test_data/hllhc15_noerrors_nobb/line_w_knobs_and_particle.json')
line.particle_ref = xt.Particles(mass0=xt.PROTON_MASS_EV, p0c=7e12)
# Build the tracker
line.build_tracker()
# Track some particles with frozen longitudinal coordinates
particles = line.build_particles(delta=1e-3, x=[-1e-3, 0, 1e-3])
line.track(particles, num_turns=10, freeze_longitudinal=True)
print(particles.delta) # gives [0.001 0.001 0.001], same as initial value
# Twiss with frozen longitudinal coordinates (needs to be 4d)
twiss = line.twiss(method='4d', freeze_longitudinal=True)
print(twiss.slip_factor) # gives 0 (no longitudinal motion)
# Track some particles with unfrozen longitudinal coordinates
particles = line.build_particles(delta=1e-3, x=[-1e-3, 0, 1e-3])
line.track(particles, num_turns=10)
print(particles.delta) # gives [0.00099218, ...], different from initial value
# Twiss with unfrozen longitudinal coordinates (can be 6d)
twiss = line.twiss(method='6d')
print(twiss.slip_factor) # gives 0.00032151, from longitudinal motion
# Complete source: xtrack/examples/freeze_longitudinal/002_freeze_individual_methods.py
Freezing longitudinal coordinates within a with
block
A context manager is also available to freeze the longitudinal coordinates within
a with
block. The normal tracking mode, updating the longitudinal
coordinates, is automatically restored when exiting the with
block, as it is
illustrated in the following example:
import json
import xtrack as xt
# import a line and add reference particle
# import a line and add reference particle
line = xt.Line.from_json(
'../../test_data/hllhc15_noerrors_nobb/line_w_knobs_and_particle.json')
line.particle_ref = xt.Particles(mass0=xt.PROTON_MASS_EV, p0c=7e12)
# Build the tracker
line.build_tracker()
# Perform a set of operations with frozen longitudinal coordinates
with xt.freeze_longitudinal(line):
# Track some particles with frozen longitudinal coordinates
particles = line.build_particles(delta=1e-3, x=[-1e-3, 0, 1e-3])
line.track(particles, num_turns=10)
print(particles.delta) # gives [0.001 0.001 0.001], same as initial value
# Twiss with frozen longitudinal coordinates (needs to be 4d)
twiss = line.twiss(method='4d')
print(twiss.slip_factor) # gives 0 (no longitudinal motion)
# 6d tracking is automatically restored when the with block is exited
# Track some particles with unfrozen longitudinal coordinates
particles = line.build_particles(delta=1e-3, x=[-1e-3, 0, 1e-3])
line.track(particles, num_turns=10)
print(particles.delta) # gives [0.00099218, ...], different from initial value
# Twiss with unfrozen longitudinal coordinates (can be 6d)
twiss = line.twiss(method='6d')
print(twiss.slip_factor) # gives 0.00032151, from longitudinal motion
# Complete source: xtrack/examples/freeze_longitudinal/001_freeze_context_manager.py
Freezing longitudinal coordinates with Line method
The xtrack.Line
class provides a method called freeze_longitudinal()
to explicitly freeze the longitudinal coordinates. The normal tracking mode,
updating the longitudinal coordinates, can be restored by calling
freeze_longitudinal(False)
. This is illustrated in the following example:
import xtrack as xt
# import a line and add reference particle
line = xt.Line.from_json(
'../../test_data/hllhc15_noerrors_nobb/line_w_knobs_and_particle.json')
line.particle_ref = xt.Particles(mass0=xt.PROTON_MASS_EV, p0c=7e12)
# Build the tracker
line.build_tracker()
# Freeze longitudinal coordinates
line.freeze_longitudinal()
# Track some particles with frozen longitudinal coordinates
particles = line.build_particles(delta=1e-3, x=[-1e-3, 0, 1e-3])
line.track(particles, num_turns=10)
print(particles.delta) # gives [0.001 0.001 0.001], same as initial value
# Twiss with frozen longitudinal coordinates (needs to be 4d)
twiss = line.twiss(method='4d')
print(twiss.slip_factor) # gives 0 (no longitudinal motion)
# Unfreeze longitudinal coordinates
line.freeze_longitudinal(False)
# Track some particles with unfrozen longitudinal coordinates
particles = line.build_particles(delta=1e-3, x=[-1e-3, 0, 1e-3])
line.track(particles, num_turns=10)
print(particles.delta) # gives [0.00099218, ...], different from initial value
# Twiss with unfrozen longitudinal coordinates (can be 6d)
twiss = line.twiss(method='6d')
print(twiss.slip_factor) # gives 0.00032151, from longitudinal motion
# Complete source: xtrack/examples/freeze_longitudinal/000_freeze_unfreeze_explicit.py
Time-dependent line properties
Time dependent knobs
To simulate the effect of time-changing properties of the beam-line it is possible
to control any lattice element attribute with a time-dependent function.
For this purpose, the variable t_turn_s
provides the time in seconds since
the start of the simulation and is updated automatically during tracking.
The simulation of an orbit bump driven by a sinusoidal
function is shown in the following example.
import numpy as np
import matplotlib.pyplot as plt
import xtrack as xt
# Build a simple ring
env = xt.Environment()
pi = np.pi
lbend = 3
line = env.new_line(components=[
# Three dipoles to make a closed orbit bump
env.new('d0.1', xt.Drift, length=0.05),
env.new('bumper_0', xt.Bend, length=0.05, k0=0, h=0),
env.new('d0.2', xt.Drift, length=0.3),
env.new('bumper_1', xt.Bend, length=0.05, k0=0, h=0),
env.new('d0.3', xt.Drift, length=0.3),
env.new('bumper_2', xt.Bend, length=0.05, k0=0, h=0),
env.new('d0.4', xt.Drift, length=0.2),
# Simple ring with two FODO cells
env.new('mqf.1', xt.Quadrupole, length=0.3, k1=0.1),
env.new('d1.1', xt.Drift, length=1),
env.new('mb1.1', xt.Bend, length=lbend, k0=pi / 2 / lbend, h=pi / 2 / lbend),
env.new('d2.1', xt.Drift, length=1),
env.new('mqd.1', xt.Quadrupole, length=0.3, k1=-0.7),
env.new('d3.1', xt.Drift, length=1),
env.new('mb2.1', xt.Bend, length=lbend, k0=pi / 2 / lbend, h=pi / 2 / lbend),
env.new('d3.4', xt.Drift, length=1),
env.new('mqf.2', xt.Quadrupole, length=0.3, k1=0.1),
env.new('d1.2', xt.Drift, length=1),
env.new('mb1.2', xt.Bend, length=lbend, k0=pi / 2 / lbend, h=pi / 2 / lbend),
env.new('d2.2', xt.Drift, length=1),
env.new('mqd.2', xt.Quadrupole, length=0.3, k1=-0.7),
env.new('d3.2', xt.Drift, length=1),
env.new('mb2.2', xt.Bend, length=lbend, k0=pi / 2 / lbend, h=pi / 2 / lbend),
])
kin_energy_0 = 50e6 # 50 MeV
line.particle_ref = xt.Particles(energy0=kin_energy_0 + xt.PROTON_MASS_EV, # total energy
mass0=xt.PROTON_MASS_EV)
# Twiss
tw = line.twiss(method='4d')
# Power the correctors to make a closed orbit bump
line.vars['bumper_strength'] = 0.
line.element_refs['bumper_0'].k0 = -line.vars['bumper_strength']
line.element_refs['bumper_1'].k0 = 2 * line.vars['bumper_strength']
line.element_refs['bumper_2'].k0 = -line.vars['bumper_strength']
# Drive the correctors with a sinusoidal function
T_sin = 100e-6
sin = line.functions.sin
line.vars['bumper_strength'] = (0.1 * sin(2 * np.pi / T_sin * line.vars['t_turn_s']))
# --- Probe behavior with twiss at different t_turn_s ---
t_test = np.linspace(0, 100e-6, 15)
tw_list = []
bumper_0_list = []
for tt in t_test:
line.vars['t_turn_s'] = tt
bumper_0_list.append(line.element_refs['bumper_0'].k0) # Inspect bumper
tw_list.append(line.twiss(method='4d')) # Twiss
# Plot
plt.close('all')
plt.figure(1, figsize=(6.4*1.2, 4.8*0.85))
ax1 = plt.subplot(1, 2, 1, xlabel='Time [us]', ylabel='x at bump center [mm]')
ax2 = plt.subplot(1, 2, 2, xlabel='s [m]', ylabel='x [mm]')
colors = plt.cm.jet(np.linspace(0, 1, len(t_test)))
for ii, tt in enumerate(t_test):
tw_tt = tw_list[ii]
ax1.plot(tt * 1e6, tw_tt['x', 'bumper_1'] * 1e3, 'o', color=colors[ii])
ax2.plot(tw_tt.s, tw_tt.x * 1e3, color=colors[ii], label=f'{tt * 1e6:.1f} us')
ax2.set_xlim(0, tw_tt['s', 'bumper_2'] + 1)
plt.subplots_adjust(left=.1, right=.97, top=.92, wspace=.27)
ax2.legend()
# --- Track particles with time-dependent bumpers ---
num_particles = 100
num_turns = 1000
# Install monitor in the middle of the bump
monitor = xt.ParticlesMonitor(num_particles=num_particles,
start_at_turn=0,
stop_at_turn =num_turns)
line.discard_tracker()
line.insert_element('monitor', monitor, index='bumper_1')
line.build_tracker()
# Generate particles
particles = line.build_particles(
x=np.random.uniform(-0.1e-3, 0.1e-3, num_particles),
px=np.random.uniform(-0.1e-3, 0.1e-3, num_particles))
# Enable time-dependent variables
line.enable_time_dependent_vars = True
# Track
line.track(particles, num_turns=num_turns, with_progress=True)
# Plot
plt.figure(2, figsize=(6.4*0.8, 4.8*0.8))
plt.plot(1e6 * monitor.at_turn.T * tw.T_rev0, 1e3 * monitor.x.T,
lw=1, color='r', alpha=0.05)
plt.plot(1e6 * monitor.at_turn[0, :] * tw.T_rev0, 1e3 * monitor.x.mean(axis=0),
lw=2, color='k')
plt.xlabel('Time [us]')
plt.ylabel('x [mm]')
plt.show()
# Complete source: xtrack/examples/toy_ring/006a_dynamic_bump_sin.py
# Control the bumpers with a piece-wise linear function of time
# (outside the defined time interval, first and last point are held constant)
line.functions['my_fun'] = xt.FunctionPieceWiseLinear(
x=np.array([10, 40, 70, 100]) * 1e-6, # time in s
y=np.array([0, 0.5, 0.5, 1.0]) # value
)
line.vars['bumper_strength'] = 0.1 * line.functions['my_fun'](line.vars['t_turn_s'])
# Complete source: xtrack/examples/toy_ring/006b_dynamic_bump_piece_wise_linear.py
Instead of built-in functions, it is also possible to use piece-wise linear functions defined by the user. For example, the function in the figure below can be obtained with the following code:
It is also possible to combine piece-wise linear functions with built-in functions in the same expression. For example, the following code can be used to generate a modulated sinusoidal function:
# Control the bumpers with a modulated sinusoidal function of time
line.functions['pulse'] = xt.FunctionPieceWiseLinear(
x=np.array([10, 50, 150, 190]) * 1e-6, # s
y=np.array([0, 1, 1, 0]) # knob value
)
T_sin = 10e-6
sin = line.functions.sin
line.vars['bumper_strength'] = (0.1 * line.functions['pulse'](line.vars['t_turn_s'])
* sin(2 * np.pi / T_sin * line.vars['t_turn_s']))
# Complete source: xtrack/examples/toy_ring/006c_dynamic_bump_sin_env.py
Multisetters
When the number of beam elements to be changed is very large, this can introduce
a significant overhead on the simulation time. For these cases it is possible to use
the xtrack.MultiSetter
class, to perform the changes in a more efficient way,
bypassing the expressions and acting directly on the element properties.
Furthermore, the MultiSetter
stores the
memory addresses of the quantities to be changed and performs the changes with a single
compiled kernel, using multithreading when allowed by the context.
The following example shows how to use the MultiSetter
to apply a sinusoidal
ripple to the strength of several quadropoles of a synchrotron.
See also: xtrack.MultiSetter
import numpy as np
from cpymad.madx import Madx
import xobjects as xo
import xtrack as xt
import xpart as xp
# Choose a context
ctx = xo.ContextCpu()
# Import SPS lattice and build a tracker
mad = Madx()
seq_name = 'sps'
mad.call('../../test_data/sps_w_spacecharge/sps_thin.seq')
mad.use(seq_name)
madtw = mad.twiss()
line = xt.Line.from_madx_sequence(mad.sequence[seq_name])
line.particle_ref = xt.Particles(p0c=400e9, mass0=xt.PROTON_MASS_EV)
line.build_tracker(_context=ctx)
# Switch on RF and twiss
line['acta.31637'].voltage = 7e6
line['acta.31637'].lag = 180.
twxt = line.twiss()
# Get revolution period
T_rev = twxt['T_rev0']
# Extract list of elements to trim (all focusing quads)
tt = line.get_table()
elements_to_trim = tt.rows[tt.element_type == 'Multipole'].rows['qf.*'].name
# => contains ['qf.52010', 'qf.52210', 'qf.52410', 'qf.52610', 'qf.52810',
# 'qf.53010', 'qf.53210', 'qf.53410', 'qf.60010', 'qf.60210', ...]
# Build a custom setter
qf_setter = xt.MultiSetter(line, elements_to_trim,
field='knl', index=1 # we want to change knl[1]
)
# Get the initial values of the quad strength
k1l_0 = qf_setter.get_values()
# Generate particles to be tracked
# (we choose to match the distribution without accounting for spacecharge)
particles = xp.generate_matched_gaussian_bunch(_context=ctx,
num_particles=100, total_intensity_particles=1e10,
nemitt_x=3e-6, nemitt_y=3e-6, sigma_z=15e-2,
line=line)
# Define amplitude and phase of the quadrupole ripple
f_quad = 50. # Hz
A_quad = 0.01 # relative amplitude
# Track the particles and apply turn-by-turn change to all selected quads
num_turns = 2000
check_trim = []
for ii in range(num_turns):
if ii % 100 == 0: print(f'Turn {ii} of {num_turns}')
# Change the strength of the quads
k1l = k1l_0 * (1 + A_quad * np.sin(2*np.pi*f_quad*ii*T_rev))
qf_setter.set_values(k1l)
# Track one turn
line.track(particles)
# Log the strength of one quad to check
check_trim.append(ctx.nparray_from_context_array(line['qf.52010'].knl)[1])
# Plot the evolution of the quad strength
import matplotlib.pyplot as plt
plt.plot(check_trim)
plt.show()
# Complete source: xtrack/examples/multisetter/000_sps_50hz_ripple.py
Exciter beam element
The beam element xtrack.Exciter
provides a model for a transverse exciter as a time-dependent thin multipole.
By providing an array of samples, and the sampling frequency, the element can provide an arbitrary excitation waveform.
This can be used for RFKO slow extraction, excitation tune measurement, power supply ripples, etc.
The given multipole components knl and ksl (normal and skew respectively) are multiplied according to an array of samples which allows for arbitrary time dependance:
To provide for an arbitrary frequency spectrum, the variations are not assumed to be slow compared to the revolution frequency \(f_{rev}\), and the particle arrival time is taken into account when determining the sample index
where \(\zeta=(s-\beta_0\cdot c_0\cdot t)\) is the longitudinal coordinate of the particle, \(\beta_0\) is the relativistic beta factor of the particle, \(c_0\) is the speed of light, \(n\) is the current turn number, \(f_{rev}\) is the revolution frequency, and \(f_{samp}\) is the sample frequency.
The excitation starts with the first sample when the reference particle arrives at the element in \(n_0\)
The samples can therefore be considered as a waveform sampled at the sampling frequency. To compute a sample for a sinusoidal excitation at frequency f_ex using NumPy:
samples[i] = np.sin(2*np.pi*f_ex*i/sampling_freq)
To generate the full samples array at sampling_freq sampling frequency, for n_turns turns at f_rev revolution frequency:
total_time = n_turns / f_rev
time = np.arange(0, total_time, 1/sampling_freq)
samples = np.sin(2*np.pi*f_ex*time)
To generate a chirp array at sampling_freq, between frequencies f_start and f_stop, lasting n_turns turns at f_rev revolution frequency:
chirp_time = n_turns / f_rev
time = np.arange(0, chirp_time, 1/sampling_freq)
samples = sp.signal.chirp(time, f_start, chirp_time, f_stop)
To then define an Exciter element with the custom waveform (array of samples at sampling frequency sampling freq) and normal and skew components KNL and KSL:
# Create beam element
exciter = xt.Exciter(_context = ctx,
samples = samples,
sampling_frequency = sampling_freq,
duration = None, # defaults to waveform duration
frev = f_rev,
start_turn = 0, # default, seconds
knl = KNL,
ksl = KSL,
)
# Add it to the line for tracking as usual
line.insert_element(
element = exciter,
name = 'RF_KO_EXCITER',
index = 42,
)
The optional parameter duration (seconds) may be used to repeat (or truncate) the excitation waveform. It defaults to len(samples)/sampling_freq, the duration of samples.
The element also provides the read-only parameter order, the multipole order, equal to the order of the largest non-zero multipole component knl or ksl.
Acceleration
An energy ramp with an arbitrary time dependence of the beam energy can be
simulated by attaching and xtrack.EnergyProgram
object to the line, as
shown in the example below.
import numpy as np
from cpymad.madx import Madx
import xtrack as xt
# Import a line and build a tracker
line = xt.Line.from_json(
'../../test_data/psb_injection/line_and_particle.json')
e_kin_start_eV = 160e6
line.particle_ref = xt.Particles(mass0=xt.PROTON_MASS_EV, q0=1.,
energy0=xt.PROTON_MASS_EV + e_kin_start_eV)
line.build_tracker()
# User-defined energy ramp
t_s = np.array([0., 0.0006, 0.0008, 0.001 , 0.0012, 0.0014, 0.0016, 0.0018,
0.002 , 0.0022, 0.0024, 0.0026, 0.0028, 0.003, 0.01])
E_kin_GeV = np.array([0.16000000,0.16000000,
0.16000437, 0.16001673, 0.16003748, 0.16006596, 0.16010243, 0.16014637,
0.16019791, 0.16025666, 0.16032262, 0.16039552, 0.16047524, 0.16056165,
0.163586])
# Attach energy program to the line
line.energy_program = xt.EnergyProgram(t_s=t_s, kinetic_energy0=E_kin_GeV*1e9)
# Plot energy and revolution frequency vs time
t_plot = np.linspace(0, 10e-3, 20)
E_kin_plot = line.energy_program.get_kinetic_energy0_at_t_s(t_plot)
f_rev_plot = line.energy_program.get_frev_at_t_s(t_plot)
import matplotlib.pyplot as plt
plt.close('all')
plt.figure(1, figsize=(6.4 * 1.5, 4.8))
ax1 = plt.subplot(2,2,1)
plt.plot(t_plot * 1e3, E_kin_plot * 1e-6)
plt.ylabel(r'$E_{kin}$ [MeV]')
ax2 = plt.subplot(2,2,3, sharex=ax1)
plt.plot(t_plot * 1e3, f_rev_plot * 1e-3)
plt.ylabel(r'$f_{rev}$ [kHz]')
plt.xlabel('t [ms]')
# Setup frequency of the RF cavity to stay on the second harmonic of the
# revolution frequency during the acceleration
t_rf = np.linspace(0, 3e-3, 100) # time samples for the frequency program
f_rev = line.energy_program.get_frev_at_t_s(t_rf)
h_rf = 2 # harmonic number
f_rf = h_rf * f_rev # frequency program
# Build a function with these samples and link it to the cavity
line.functions['fun_f_rf'] = xt.FunctionPieceWiseLinear(x=t_rf, y=f_rf)
line.element_refs['br.c02'].frequency = line.functions['fun_f_rf'](
line.vars['t_turn_s'])
# Setup voltage and lag
line.element_refs['br.c02'].voltage = 3000 # V
line.element_refs['br.c02'].lag = 0 # degrees (below transition energy)
# When setting line.vars['t_turn_s'] the reference energy and the rf frequency
# are updated automatically
line.vars['t_turn_s'] = 0
line.particle_ref.kinetic_energy0 # is 160.00000 MeV
line['br.c02'].frequency # is 1983931.935 Hz
line.vars['t_turn_s'] = 3e-3
line.particle_ref.kinetic_energy0 # is 160.56165 MeV
line['br.c02'].frequency # is 1986669.0559674294
# Back to zero for tracking!
line.vars['t_turn_s'] = 0
# Track a few particles to visualize the longitudinal phase space
p_test = line.build_particles(x_norm=0, zeta=np.linspace(0, line.get_length(), 101))
# Enable time-dependent variables (t_turn_s and all the variables that depend on
# it are automatically updated at each turn)
line.enable_time_dependent_vars = True
# Track
line.track(p_test, num_turns=9000, turn_by_turn_monitor=True, with_progress=True)
mon = line.record_last_track
# Plot
plt.subplot2grid((2,2), (0,1), rowspan=2)
plt.plot(mon.zeta[:, -2000:].T, mon.delta[:, -2000:].T, color='C0')
plt.xlabel(r'$\zeta$ [m]')
plt.ylabel('$\delta$')
plt.xlim(-40, 30)
plt.ylim(-0.0025, 0.0025)
plt.title('Last 2000 turns')
plt.subplots_adjust(left=0.08, right=0.95, wspace=0.26)
plt.show()
# Complete source: xtrack/examples/acceleration/001_energy_ramp.py
Off-momentum from RF frequency change
The xtrack default RF are synchronized with the reference particle (simply
because of the coordinate choice). For this reason a change of the RF frequency
does not result in a change in the revolution frequency. To obtain a change in
the revolution frequency (and hence in the momentum) it is necessary to
introduce explicitly a time delay, using the element xtrack.ZetaShift
.
This is illustrated in the following example:
import numpy as np
import xtrack as xt
# Load a line and build a tracker
line = xt.Line.from_json(
'../../test_data/hllhc15_noerrors_nobb/line_and_particle.json')
line.particle_ref = xt.Particles(mass0=xt.PROTON_MASS_EV, q0=1, energy0=7e12)
line.build_tracker()
# Frequency trim to make
df_hz = 50 # Frequency trim
# Twiss
tw0 = line.twiss()
# Compute corresponding delay to be introduced in the line:
#
# T_rev = h_rf/f_rf
# dt = h_rf/(f_rf + df_hz) - h_rf/f_rf = h_rf/f_rf (1/(1+df_hz/f_rf) - 1)
# ~= h_rf/f_rf * (1 - df_hz/f_rf -1)
# = -h_rf/(f_rf^2) * df_hz
# = -T_rev / f_rf * df_hz
# dzeta = -beta0 * clight * dt = circumference * df_hz / f_rf
h_rf = 35640
f_rf = h_rf/tw0.T_rev0
beta0 = line.particle_ref.beta0[0]
dzeta = tw0.circumference * df_hz / f_rf
# Append delay element to the line
line.unfreeze()
line.append_element(element=xt.ZetaShift(dzeta=dzeta), name='zeta_shift')
line.build_tracker()
# Twiss
tw1 = line.twiss()
# Expected momentum from slip factor (eta = -df_rev / f_rev / delta)
f_rev = 1/tw0.T_rev0
df_rev = df_hz / h_rf
eta = tw0.slip_factor
delta_expected = -df_rev / f_rev / eta
tw1.delta[0] # is -0.00035789
delta_expected # is -0.00035803
# Check that particle stays off momentum in multi-turn tracking
p = tw1.particle_on_co.copy()
p_test = line.build_particles(x_norm=0,
delta=np.linspace(delta_expected - 8e-4, delta_expected + 8e-4, 60))
# p_test = line.build_particles(x_norm=0, delta=delta_expected,
# zeta=np.linspace(0, 0.5, 101))
# Track
line.track(p_test, num_turns=1000, turn_by_turn_monitor=True, with_progress=True)
mon = line.record_last_track
# Plot
import matplotlib.pyplot as plt
plt.close('all')
plt.figure()
plt.plot(mon.zeta[:, :].T, mon.delta[:, :].T * 1e4, color='C0')
plt.xlabel(r'$\zeta$ [m]')
plt.ylabel('$\delta$ [$10^{-4}$]')
plt.xlim(-.8, .8)
plt.ylim(delta_expected * 1e4 - 4, delta_expected * 1e4 + 4)
plt.axhline(delta_expected * 1e4, color='r', linestyle='--',
label='$\delta$ expected')
plt.show()
# Complete source: xtrack/examples/radial_steering/000_radial_steering.py
Optimize line for multi-turn tracking speed
When the tracker is used for multi-turn tracking, the line and tracker can be
optimized to increase the tracking speed. This can be done with the
Line.optimize_for_tracking()
method, as shown in the following example.
See also: Line.optimize_for_tracking()
import json
import time
import numpy as np
import xtrack as xt
import xpart as xp
#################################
# Load a line and build tracker #
#################################
line = xt.Line.from_json(
'../../test_data/hllhc15_noerrors_nobb/line_w_knobs_and_particle.json')
line.particle_ref = xt.Particles(mass0=xt.PROTON_MASS_EV, p0c=7e12)
line.build_tracker()
###########################
# Generate some particles #
###########################
particles = line.build_particles(
x_norm=np.linspace(-2, 2, 1000), y_norm=0.1, delta=3e-4,
nemitt_x=2.5e-6, nemitt_y=2.5e-6)
p_no_optimized = particles.copy()
p_optimized = particles.copy()
##############################
# Track without optimization #
##############################
num_turns = 10
line.track(p_no_optimized, num_turns=num_turns, time=True)
t_not_optimized = line.time_last_track
####################
# Optimize tracker #
####################
line.optimize_for_tracking()
# This performs the following actions (physics model is unchanged):
# - Disables xdeps expressions
# - Removes inactive multipoles
# - Merges consecutive multipoles
# - Removes drifts with zero length
# - Merges consecutive drifts
###########################
# Track with optimization #
###########################
line.track(p_optimized, num_turns=num_turns, time=True)
t_optimized = line.time_last_track
#################
# Compare times #
#################
num_particles = len(p_no_optimized.x)
print(f'Time not optimized {t_not_optimized*1e6/num_particles/num_turns:.1f} us/part/turn')
print(f'Time optimized {t_optimized*1e6/num_particles/num_turns:.1f} us/part/turn')
###################################
# Check that result are identical #
###################################
assert np.all(p_no_optimized.state == 1)
assert np.all(p_optimized.state == 1)
assert np.allclose(p_no_optimized.x, p_optimized.x, rtol=0, atol=1e-14)
assert np.allclose(p_no_optimized.y, p_optimized.y, rtol=0, atol=1e-14)
assert np.allclose(p_no_optimized.px, p_optimized.px, rtol=0, atol=1e-14)
assert np.allclose(p_no_optimized.py, p_optimized.py, rtol=0, atol=1e-14)
assert np.allclose(p_no_optimized.zeta, p_optimized.zeta, rtol=0, atol=1e-14)
assert np.allclose(p_no_optimized.delta, p_optimized.delta, rtol=0, atol=1e-14)
# Complete source: xtrack/examples/optimized_tracker/000_optimized_tracker.py