Closed orbit and trajectory correction

Xsuite offers the possibility to correct the beam closed orbit for circular accelerators and the beam trajectory in transfer lines. The correction is performed using a linearized response matrix, which is built from a twiss table. The correction is computed using the singular value decomposition (SVD) of the response matrix (least squares solution) or the MICADO algorithm.

In the case of rings, in order to proceed with the correction with the approach described above, it is in necessary to successfully measure the closed orbit before the correction. In certain cases, when strong lattice perturbations are present (e.g. field errors or large element misalignments), the closed orbit search might fail. In such cases, the user can use a threading capability to perform a first correction of the trajectory, after which the closed orbit search can be performed.

The following sections illustrate the capabilities of the Xsuite trajectory correction module.

Basic usage

See also: xtrack.Line.correct_trajectory()

The following example shows how to correct the closed orbit of a ring using the least squares approach.

import xtrack as xt
import numpy as np

line = xt.Line.from_json(
    '../../test_data/hllhc15_thick/lhc_thick_with_knobs.json')
tt = line.get_table()

# Define elements to be used as monitors for orbit correction
# (for LHC all element names starting by "bpm" and not ending by "_entry" or "_exit")
tt_monitors = tt.rows['bpm.*'].rows['.*(?<!_entry)$'].rows['.*(?<!_exit)$']
line.steering_monitors_x = tt_monitors.name
line.steering_monitors_y = tt_monitors.name

# Define elements to be used as correctors for orbit correction
# (for LHC all element namesstarting by "mcb.", containing "h." or "v.")
tt_h_correctors = tt.rows['mcb.*'].rows['.*h\..*']
line.steering_correctors_x = tt_h_correctors.name
tt_v_correctors = tt.rows['mcb.*'].rows['.*v\..*']
line.steering_correctors_y = tt_v_correctors.name

# Reference twiss (no misalignments)
tw_ref = line.twiss4d()

# Introduce misalignments on all quadrupoles
tt = line.get_table()
tt_quad = tt.rows[tt.element_type == 'Quadrupole']
rgen = np.random.RandomState(1) # fix seed for random number generator
shift_x = rgen.randn(len(tt_quad)) * 0.01e-3 # 0.01 mm rms shift on all quads
shift_y = rgen.randn(len(tt_quad)) * 0.01e-3 # 0.01 mm rms shift on all quads
for nn_quad, sx, sy in zip(tt_quad.name, shift_x, shift_y):
    line.element_refs[nn_quad].shift_x = sx
    line.element_refs[nn_quad].shift_y = sy

# Twiss before correction
tw_before = line.twiss4d()

# Correct orbit
orbit_correction = line.correct_trajectory(twiss_table=tw_ref)
# prints:
#
# Iteration 0, x_rms: 1.65e-03 -> 7.22e-06, y_rms: 2.23e-03 -> 1.54e-04
# Iteration 1, x_rms: 1.11e-05 -> 2.09e-06, y_rms: 1.55e-04 -> 5.54e-06
# Iteration 2, x_rms: 2.09e-06 -> 2.09e-06, y_rms: 5.54e-06 -> 2.15e-06
# Iteration 3, x_rms: 2.09e-06 -> 2.09e-06, y_rms: 2.15e-06 -> 2.13e-06

# Twiss after correction
tw_after = line.twiss4d()

# Extract correction strength
s_x_correctors = orbit_correction.x_correction.s_correctors
s_y_correctors = orbit_correction.y_correction.s_correctors
kicks_x = orbit_correction.x_correction.get_kick_values()
kicks_y = orbit_correction.y_correction.get_kick_values()

# Complete source: xtrack/examples/orbit_and_tracjectory_correction/000_closed_orbit_correction_basic.py
_images/orbit_correction_basic.png

Result of the orbit correction (SVD based) for the LHC ring.

MICADO correction

See also: xtrack.Line.correct_trajectory()

The following example shows how to correct the closed orbit of a ring using the MICADO algorithm, to achive the best possible correction while using an assigned number of correctors.

import xtrack as xt
import numpy as np

line = xt.Line.from_json(
    '../../test_data/hllhc15_thick/lhc_thick_with_knobs.json')
tt = line.get_table()

# Define elements to be used as monitors for orbit correction
# (for LHC all element names starting by "bpm" and not ending by "_entry" or "_exit")
tt_monitors = tt.rows['bpm.*'].rows['.*(?<!_entry)$'].rows['.*(?<!_exit)$']
line.steering_monitors_x = tt_monitors.name
line.steering_monitors_y = tt_monitors.name

# Define elements to be used as correctors for orbit correction
# (for LHC all element namesstarting by "mcb.", containing "h." or "v.")
tt_h_correctors = tt.rows['mcb.*'].rows['.*h\..*']
line.steering_correctors_x = tt_h_correctors.name
tt_v_correctors = tt.rows['mcb.*'].rows['.*v\..*']
line.steering_correctors_y = tt_v_correctors.name

# Reference twiss (no misalignments)
tw_ref = line.twiss4d()

# Introduce misalignments on all quadrupoles
tt = line.get_table()
tt_quad = tt.rows[tt.element_type == 'Quadrupole']
rgen = np.random.RandomState(1) # fix seed for random number generator
shift_x = rgen.randn(len(tt_quad)) * 0.01e-3 # 0.01 mm rms shift on all quads
shift_y = rgen.randn(len(tt_quad)) * 0.01e-3 # 0.01 mm rms shift on all quads
for nn_quad, sx, sy in zip(tt_quad.name, shift_x, shift_y):
    line.element_refs[nn_quad].shift_x = sx
    line.element_refs[nn_quad].shift_y = sy

# Twiss before correction
tw_before = line.twiss4d()

# Correct orbit using 5 correctors in each plane
orbit_correction = line.correct_trajectory(twiss_table=tw_ref, n_micado=5,
                                           n_iter=1)
# prints:
#
# Iteration 0, x_rms: 1.65e-03 -> 1.06e-04, y_rms: 2.25e-03 -> 1.76e-04

# Twiss after correction
tw_after = line.twiss4d()

# Extract correction strength
s_x_correctors = orbit_correction.x_correction.s_correctors
s_y_correctors = orbit_correction.y_correction.s_correctors
kicks_x = orbit_correction.x_correction.get_kick_values()
kicks_y = orbit_correction.y_correction.get_kick_values()

# Complete source: xtrack/examples/orbit_and_tracjectory_correction/001_closed_orbit_correction_micado.py
_images/orbit_correction_micado.png

Result of the orbit correction (MICADO based) for the LHC ring (5 correctors per plane are used).

Customized correction

See also: xtrack.Line.correct_trajectory()

In order to customize the correction, the user can generate a correction object without actually running the correction. This allows advanced control of the correction process. In the following example, we illustrate how to undo an existing correction and apply a new one performed with a selected number of singular values.

import xtrack as xt
import numpy as np

line = xt.Line.from_json(
    '../../test_data/hllhc15_thick/lhc_thick_with_knobs.json')
tt = line.get_table()

# Define elements to be used as monitors for orbit correction
# (for LHC all element names starting by "bpm" and not ending by "_entry" or "_exit")
tt_monitors = tt.rows['bpm.*'].rows['.*(?<!_entry)$'].rows['.*(?<!_exit)$']
line.steering_monitors_x = tt_monitors.name
line.steering_monitors_y = tt_monitors.name

# Define elements to be used as correctors for orbit correction
# (for LHC all element namesstarting by "mcb.", containing "h." or "v.")
tt_h_correctors = tt.rows['mcb.*'].rows['.*h\..*']
line.steering_correctors_x = tt_h_correctors.name
tt_v_correctors = tt.rows['mcb.*'].rows['.*v\..*']
line.steering_correctors_y = tt_v_correctors.name

# Reference twiss (no misalignments)
tw_ref = line.twiss4d()

# Introduce misalignments on all quadrupoles
tt = line.get_table()
tt_quad = tt.rows[tt.element_type == 'Quadrupole']
rgen = np.random.RandomState(1) # fix seed for random number generator
shift_x = rgen.randn(len(tt_quad)) * 0.01e-3 # 0.01 mm rms shift on all quads
shift_y = rgen.randn(len(tt_quad)) * 0.01e-3 # 0.01 mm rms shift on all quads
for nn_quad, sx, sy in zip(tt_quad.name, shift_x, shift_y):
    line.element_refs[nn_quad].shift_x = sx
    line.element_refs[nn_quad].shift_y = sy

# Create orbit correction object without running the correction
orbit_correction = line.correct_trajectory(twiss_table=tw_ref, run=False)

# Inspect singular values of the response matrices
x_sv = orbit_correction.x_correction.singular_values
y_sv = orbit_correction.y_correction.singular_values

# Twiss before correction
tw_before = line.twiss4d()

# Correct
orbit_correction.correct()
# prints:
# Iteration 0, x_rms: 1.65e-03 -> 7.22e-06, y_rms: 2.23e-03 -> 1.54e-04
# Iteration 1, x_rms: 1.11e-05 -> 2.09e-06, y_rms: 1.55e-04 -> 5.54e-06
# Iteration 2, x_rms: 2.09e-06 -> 2.09e-06, y_rms: 5.54e-06 -> 2.15e-06
# Iteration 3, x_rms: 2.09e-06 -> 2.09e-06, y_rms: 2.15e-06 -> 2.13e-06

# Remove applied correction
orbit_correction.clear_correction_knobs()

# Correct with a customized number of singular values
orbit_correction.correct(n_singular_values=(200, 210))
# prints:
#
# Iteration 0, x_rms: 1.65e-03 -> 9.38e-06, y_rms: 2.23e-03 -> 1.55e-04
# Iteration 1, x_rms: 1.24e-05 -> 6.42e-06, y_rms: 1.56e-04 -> 7.40e-06
# Iteration 2, x_rms: 6.42e-06 -> 6.42e-06, y_rms: 7.40e-06 -> 5.22e-06
# Iteration 3, x_rms: 6.42e-06 -> 6.42e-06, y_rms: 5.22e-06 -> 5.22e-06

# Twiss after correction
tw_after = line.twiss4d()

# Extract correction strength
s_x_correctors = orbit_correction.x_correction.s_correctors
s_y_correctors = orbit_correction.y_correction.s_correctors
kicks_x = orbit_correction.x_correction.get_kick_values()
kicks_y = orbit_correction.y_correction.get_kick_values()

# Complete source: xtrack/examples/orbit_and_tracjectory_correction/002_closed_orbit_correction_customize.py
_images/orbit_correction_svalues.png

Singular values of the response matrix for the LHC ring.

_images/orbit_correction_custom.png

Result of the customized orbit correction for the LHC ring.

Threading

See also: xtrack.Line.correct_trajectory()

In order to proceed with the correction with the approach described above, it is in necessary to successfully measure the closed orbit before the correction. In certain cases, when strong lattice perturbations are present (e.g. field errors or large element misalignments), the closed orbit search might fail. In such cases, the user can use a threading capability to perform a first correction of the trajectory, after which the closed orbit search can be performed. This is illustrated in the following example.

import xtrack as xt
import numpy as np

line = xt.Line.from_json(
    '../../test_data/hllhc15_thick/lhc_thick_with_knobs.json')
tt = line.get_table()

# Define elements to be used as monitors for orbit correction
# (for LHC all element names starting by "bpm" and not ending by "_entry" or "_exit")
tt_monitors = tt.rows['bpm.*'].rows['.*(?<!_entry)$'].rows['.*(?<!_exit)$']
line.steering_monitors_x = tt_monitors.name
line.steering_monitors_y = tt_monitors.name

# Define elements to be used as correctors for orbit correction
# (for LHC all element namesstarting by "mcb.", containing "h." or "v.")
tt_h_correctors = tt.rows['mcb.*'].rows['.*h\..*']
line.steering_correctors_x = tt_h_correctors.name
tt_v_correctors = tt.rows['mcb.*'].rows['.*v\..*']
line.steering_correctors_y = tt_v_correctors.name

# Reference twiss (no misalignments)
tw_ref = line.twiss4d()

# Introduce misalignments on all quadrupoles
tt = line.get_table()
tt_quad = tt.rows[tt.element_type == 'Quadrupole']
rgen = np.random.RandomState(2) # fix seed for random number generator
shift_x = rgen.randn(len(tt_quad)) * 1e-3 # 1. mm rms shift on all quads
shift_y = rgen.randn(len(tt_quad)) * 1e-3 # 1. mm rms shift on all quads
for nn_quad, sx, sy in zip(tt_quad.name, shift_x, shift_y):
    line.element_refs[nn_quad].shift_x = sx
    line.element_refs[nn_quad].shift_y = sy

# Closed twiss fails (closed orbit is not found)
# line.twiss4d()

# Create orbit correction object without running the correction
orbit_correction = line.correct_trajectory(twiss_table=tw_ref, run=False)

# Thread
threader = orbit_correction.thread(ds_thread=500., # correct in sections of 500 m
                                   rcond_short=1e-4, rcond_long=1e-4)
# prints:
#
# Stop at s=500.0, local rms  = [x: 1.44e-03 -> 3.04e-05, y: 3.47e-03 -> 6.91e-08]
# Stop at s=500.0, global rms = [x: 3.04e-05 -> 3.04e-05, y: 6.91e-08 -> 1.90e-12]
# Stop at s=1000.0, local rms  = [x: 6.44e-03 -> 9.02e-05, y: 3.45e-03 -> 6.44e-05]
# Stop at s=1000.0, global rms = [x: 1.37e-03 -> 2.19e-05, y: 3.59e-03 -> 7.25e-06]
# Stop at s=1500.0, local rms  = [x: 3.24e-03 -> 5.33e-04, y: 1.37e-03 -> 6.29e-06]
# Stop at s=1500.0, global rms = [x: 5.58e-04 -> 1.89e-05, y: 2.17e-03 -> 5.66e-06]
# ...
# Stop at s=26000.0, local rms  = [x: 7.07e-03 -> 1.66e-04, y: 5.25e-03 -> 3.62e-04]
# Stop at s=26000.0, global rms = [x: 5.46e-04 -> 5.82e-06, y: 2.53e-04 -> 2.40e-07]
# Stop at s=26500.0, local rms  = [x: 4.30e-03 -> 1.16e-03, y: 1.73e-03 -> 3.11e-06]
# Stop at s=26500.0, global rms = [x: 3.19e-04 -> 5.09e-06, y: 4.06e-04 -> 2.81e-06]
# Stop at s=26658.88, local rms  = [x: 8.09e-04 -> 3.74e-05, y: 1.11e-03 -> 2.86e-04]
# Stop at s=26658.88, global rms = [x: 3.07e-05 -> 5.07e-06, y: 5.45e-05 -> 3.52e-09]

# Closed twiss after threading (closed orbit is found)
tw_after_thread = line.twiss4d()

# Correct (with custom number of singular values)
orbit_correction.correct(n_singular_values=200)
# prints:
#
# Iteration 0, x_rms: 1.30e-03 -> 2.33e-04, y_rms: 6.91e-04 -> 2.20e-04
# Iteration 1, x_rms: 2.33e-04 -> 2.33e-04, y_rms: 2.20e-04 -> 2.16e-04

# Twiss after closed orbit correction
tw_after = line.twiss4d()

# Extract correction strength
s_x_correctors = orbit_correction.x_correction.s_correctors
s_y_correctors = orbit_correction.y_correction.s_correctors
kicks_x = orbit_correction.x_correction.get_kick_values()
kicks_y = orbit_correction.y_correction.get_kick_values()

# Complete source: xtrack/examples/orbit_and_tracjectory_correction/003_closed_orbit_correction_thread.py
_images/orbit_correction_thread.png

Closed orbit obtained after threading, followed by conventional SVD correction.

Trajectory correction for transfer lines

See also: xtrack.Line.correct_trajectory()

The following example shows how to correct the beam trajectory in a transfer line.

import numpy as np
from cpymad.madx import Madx
import xtrack as xt

mad_ti2 = Madx()
mad_ti2.call('../../test_data/sps_to_lhc_ti2/ti2.seq')
mad_ti2.call('../../test_data/sps_to_lhc_ti2/ti2_liu.str')
mad_ti2.beam()
mad_ti2.use('ti2')

line = xt.Line.from_madx_sequence(mad_ti2.sequence['ti2'])
line.particle_ref = xt.Particles(p0c=450e9, mass0=xt.PROTON_MASS_EV, q0=1)
tt = line.get_table()

# Define elements to be used as monitors for orbit correction
# (in this case all element names starting by "bpm" and not ending by "_entry" or "_exit")
tt_monitors = tt.rows['bpm.*'].rows['.*(?<!_entry)$'].rows['.*(?<!_exit)$']
line.steering_monitors_x = tt_monitors.name
line.steering_monitors_y = tt_monitors.name

# Define elements to be used as correctors for orbit correction
# (in this case all element names starting by "mci.", containing "h." or "v.")
tt_h_correctors = tt.rows['mci.*'].rows['.*h\..*']
line.steering_correctors_x = tt_h_correctors.name
tt_v_correctors = tt.rows['mci.*'].rows['.*v\..*']
line.steering_correctors_y = tt_v_correctors.name

# Initial conditions from upstream ring
init = xt.TwissInit(betx=27.77906807, bety=120.39920690,
                     alfx=0.63611880, alfy=-2.70621900,
                     dx=-0.59866300, dpx=0.01603536)

# Reference twiss (no misalignments)
tw_ref = line.twiss4d(start='ti2$start', end='ti2$end', init=init)

# Introduce misalignments on all quadrupoles
tt = line.get_table()
tt_quad = tt.rows['mqi.*']
rgen = np.random.RandomState(2) # fix seed for random number generator
shift_x = rgen.randn(len(tt_quad)) * 0.1e-3 # 0.1 mm rms shift on all quads
shift_y = rgen.randn(len(tt_quad)) * 0.1e-3 # 0.1 mm rms shift on all quads
for nn_quad, sx, sy in zip(tt_quad.name, shift_x, shift_y):
    line.element_refs[nn_quad].shift_x = sx
    line.element_refs[nn_quad].shift_y = sy

# Twiss before correction
tw_before = line.twiss4d(start='ti2$start', end='ti2$end', init=init)

# Correct trajectory
correction = line.correct_trajectory(twiss_table=tw_ref, start='ti2$start', end='ti2$end')
# prints:
#
# Iteration 0, x_rms: 2.01e-03 -> 1.80e-04, y_rms: 6.94e-04 -> 1.23e-04
# Iteration 1, x_rms: 1.80e-04 -> 1.80e-04, y_rms: 1.23e-04 -> 1.23e-04

# Twiss after correction
tw_after = line.twiss4d(start='ti2$start', end='ti2$end', init=init)

# Extract correction strength
s_x_correctors = correction.x_correction.s_correctors
s_y_correctors = correction.y_correction.s_correctors
kicks_x = correction.x_correction.get_kick_values()
kicks_y = correction.y_correction.get_kick_values()

# Complete source: xtrack/examples/orbit_and_tracjectory_correction/005_transfer_line_correction.py
_images/orbit_correction_ti2.png

Trajectory correction for one of the LHC injection lines (TI2).