Collimation
Introduction
In Xsuite, collimation is added to the simulations by the Xcoll package.
The collimators themselves are created as instances of xcoll.EverestCollimator, xcoll.EverestCrystal and xcoll.BlackAbsorber. In addition, we also have xcoll.EverestBlock. The Xcoll package also describes the interaction between the particles and the collimators in the simulation.
Loss maps are created after the simulations to assess the performance of the LHC collimation system. They give information about where the beam losses are located in the LHC. The loss map, after one simulation, is created as an instance of xcoll.LossMap.
Collimator objects
BaseCollimator
Since xcoll.EverestCollimator, xcoll.EverestCrystal and xcoll.BlackAbsorber describes different kinds of collimators they all inherit from the abstract class BaseCollimator which contains all the basic attributes of a collimator. Some of these attributes, or fields, are made accessible to the C tracking code as seen in the code block below.
class BaseCollimator(xt.BeamElement):
_xofields = {
'inactive_front': xo.Float64, # Drift before jaws
'active_length': xo.Float64, # Length of jaws
'inactive_back': xo.Float64, # Drift after jaws
'jaw_L': xo.Float64, # left jaw (distance to ref)
'jaw_R': xo.Float64, # right jaw
'ref_x': xo.Float64, # center of collimator reference frame
'ref_y': xo.Float64,
'sin_zL': xo.Float64, # angle of left jaw
'cos_zL': xo.Float64,
'sin_zR': xo.Float64, # angle of right jaw
'cos_zR': xo.Float64,
'sin_yL': xo.Float64, # tilt of left jaw (around jaw midpoint)
'cos_yL': xo.Float64,
'tan_yL': xo.Float64,
'sin_yR': xo.Float64, # tilt of right jaw (around jaw midpoint)
'cos_yR': xo.Float64,
'tan_yR': xo.Float64,
'_side': xo.Int8, # is it a onesided collimator?
'active': xo.Int8
}
However, to make things more user-friendly the BaseCollimator has more properties defined in the class itself. These properties uses the fields from the code block above to define properties such as jaw_LU(), angle() and tilt().
The collimator jaws are separated into the left and the right jaw, as seen defined in the code block above as well. In the figure below, you can see how ‘jaw_L’ and ‘jaw_R’ are defined. However, sometimes it is necessary to know, or change, the position of the corners of the collimator. The first pair of corners are called jaw_LU() and jaw_RU(), where the ‘U’ stands for upstream. The two remaining corners are defined as jaw_LD() and jaw_RD(), where the ‘D’ stands for downstream. Note that when setting the value of e.g. jaw_RU() then jaw_RD() is kept fixed while both ‘jaw_R’ and the tilt changes.
collimatorObject.jaw
collimatorObject.jaw_LU
collimatorObject.jaw_RU
collimatorObject.jaw_RD
collimatorObject.jaw_RU
In addition, it is convenient to have the angle and the tilt of the jaws. The angle is the angle of the tilt of the jaws in the xy-plane while the tilt is the angle of the tilt of the jaws in the x’z-plane. For clarity, the angle is shown as alpha, and the tilt as theta in the images. Note that the tilt is in radians while the angle is in degrees.
collimatorObject.angle
collimatorObject.angle_L
collimatorObject.angle_R
collimatorObject.tilt
collimatorObject.tilt_L
collimatorObject.tilt_R
Furthermore, it is also possible to get (and set) the length and the side of the collimator.
collimatorObject.side
collimatorObject.length
EverestCollimator
The EverestCollimator contains all the same fields as BaseCollimator, but in addition it also has these:
class EverestCollimator(BaseCollimator):
_xofields = { **BaseCollimator._xofields,
'_material': Material,
'rutherford_rng': xt.RandomRutherford,
'_tracking': xo.Int8
}
The new field of interest for the user is the material of the collimator. This is accessed by the method EverestCollimator.material(). The material itself is an instance of xcoll.materials.Material.
EverestCrystal
EverestCrystal has the same fields as BaseCollimator, but in addition it has some extra fields which are needed to describe the crystal collimator. Note also that the material used for the crystal collimator is not the same as for xcoll.EverestCollimator, but an instance of xcoll.materials.CrystalMaterial.
class EverestCrystal(BaseCollimator):
_xofields = { **BaseCollimator._xofields,
'align_angle': xo.Float64, # = - sqrt(eps/beta)*alpha*nsigma
'_bending_radius': xo.Float64,
'_bending_angle': xo.Float64,
'_critical_angle': xo.Float64,
'xdim': xo.Float64,
'ydim': xo.Float64,
'thick': xo.Float64,
'miscut': xo.Float64,
'_orient': xo.Int8,
'_material': CrystalMaterial,
'rutherford_rng': xt.RandomRutherford,
'_tracking': xo.Int8
}
These new fields are accessed through the methods:
EverestCrystal.critical_angle
EverestCrystal.bending_radius
EverestCrystal.bending_angle
EverestCrystal.material
EverestCrystal.lattice
BlackAbsorber
xcoll.BlackAbsorber is different from xcoll.EverestCollimator and EverestCrystal. The absorber has only one extra field compared to the xcoll.BaseCollimator:
class BlackAbsorber(BaseCollimator):
_xofields = { **BaseCollimator._xofields,
'_tracking': xo.Int8
}
BaseBlock and EverestBlock
xcoll.EverestBlock, which inherit xcoll.BaseBlock, describes a block with an infinite transversal length. This class has the fields:
class EverestBlock(BaseBlock):
_xofields = { **BaseBlock._xofields,
'_material': Material,
'rutherford_rng': xt.RandomRutherford,
'_tracking': xo.Int8,
'_only_mcs': xo.Int8
}
Furthermore, xcoll.EverestBlock needs a material, which is an instance of xcoll.materials.Material, and that can be accessed through
Creating a Collimator or Block object
A collimator (or block) object can be created in two different ways; either directly with the class or by loading from file. For example:
import xcoll as xc
block = xc.EverestBlock(length=1., material=xc.materials.Tungsten)
collimator = xc.EverestCollimator(length=1., material=xc.materials.Tungsten)
collimator_crystal = xc.EverestCrystal(length=1., material=xc.materials.SiliconCrystal)
black_absorber = xc.BlackAbsorber(length=1., material=xc.materials.Graphite)
Or, by using the CollimationManager to load from file:
path_in = xc._pkg_root.parent / 'examples'
coll_manager = xc.CollimatorManager.from_yaml(path_in / 'colldb' / f'lhc_run3.yaml',
line=line, beam=beam, _context=context)
# Install collimators in line as black absorbers
coll_manager.install_everest_collimators(verbose=True)
Generating particles on a collimator
For some collimation studies it is convenient to generate a initial pencil distribution on a collimator. Xcoll has its own function for this xcoll.generate_pencil_on_collimator(). An example is shown below.
import xcoll as xc
import xpart as xp
import xtrack as xt
import numpy as np
import json
beam = 1
plane = 'H'
num_turns = 200
num_particles = 10000
path_in = xc._pkg_root.parent / 'examples'
# Load from json
with open(os.devnull, 'w') as fid:
with contextlib.redirect_stdout(fid):
line = xt.Line.from_json(path_in / 'machines' / f'lhc_run3_b{beam}.json')
# Initialise collmanager
coll_manager = xc.CollimatorManager.from_yaml(path_in / 'colldb' / f'lhc_run3.yaml',
line=line, beam=beam, _context=context)
# Install collimators into line
coll_manager.install_everest_collimators(verbose=True)
# Aperture model check
print('\nAperture model check after introducing collimators:')
with open(os.devnull, 'w') as fid:
with contextlib.redirect_stdout(fid):
df_with_coll = line.check_aperture()
assert not np.any(df_with_coll.has_aperture_problem)
# Build the tracker
coll_manager.build_tracker()
# Set openings
coll_manager.set_openings()
tcp = f"tcp.{'c' if plane=='H' else 'd'}6{'l' if int(beam)==1 else 'r'}7.b{beam}"
emittance = coll_manager.colldb.emittance
beta_gamma_rel = coll_manager.colldb._beta_gamma_rel
part = xc.generate_pencil_on_collimator(line=line, emittance=emittance, beta_gamma_rel=bet,
collimator_name=tcp, num_particles=num_particles)
Lossmaps
Lossmaps are created as instances of xcoll.LossMap. The lossmap itself and its summary are calculated when the object is created. It is also possible to save both the lossmap and summary to file with xcoll.LossMap.to_json() and xcoll.LossMap.save_summary(). For example:
Loss location refinement
In Xtrack simulations particles are lost at defined aperture elements (e.g.
xtrack.LimitRect, xtrack.LimitEllipse, xtrack.LimitRectEllipse,
xtrack.LimitPolygon). A more accurate estimate of the loss locations can be
obtained after the tracking is finished using the
xtrack.LossLocationRefinement tool . The tool builds
an interpolated aperture model between the aperture elements and backtracks the
particles in order to identify the impact point. The following example illustrates
how to use this feature.
See also: xtrack.LossLocationRefinement
import numpy as np
import xtrack as xt
import xobjects as xo
###################
# Build test line #
###################
ctx = xo.context_default
buf = ctx.new_buffer()
# We build a test line having two aperture elements which are shifted and
# rotated w.r.t. the accelerator reference frame.
# Define aper_0
aper_0 = xt.LimitEllipse(_buffer=buf, a=2e-2, b=1e-2)
shift_aper_0 = (1e-2, 0.5e-2)
rot_deg_aper_0 = 10.
# Define aper_1
aper_1 = xt.LimitRect(_buffer=buf, min_x=-1e-2, max_x=1e-2,
min_y=-2e-2, max_y=2e-2)
shift_aper_1 = (-5e-3, 1e-2)
rot_deg_aper_1 = 10.
aper_1.shift_x = shift_aper_1[0]
aper_1.shift_y = shift_aper_1[1]
aper_1.rot_s_rad = np.deg2rad(rot_deg_aper_1)
# aper_0_sandwitch
line_aper_0 = xt.Line(
elements=[xt.XYShift(_buffer=buf, dx=shift_aper_0[0], dy=shift_aper_0[1]),
xt.SRotation(_buffer=buf, angle=rot_deg_aper_0),
aper_0,
xt.Multipole(_buffer=buf, knl=[0.001]),
xt.SRotation(_buffer=buf, angle=-rot_deg_aper_0),
xt.XYShift(_buffer=buf, dx=-shift_aper_0[0], dy=-shift_aper_0[1])])
line_aper_0.build_tracker(_buffer=buf)
# aper_1_sandwitch
line_aper_1 = xt.Line(
elements=[aper_1,
xt.Multipole(_buffer=buf, knl=[0.001])
])
line_aper_1.build_tracker(_buffer=buf)
#################
# Build tracker #
#################
line=xt.Line(
elements = ((xt.Drift(_buffer=buf, length=0.5),)
+ line_aper_0.elements
+ (xt.Drift(_buffer=buf, length=1),
xt.Drift(_buffer=buf, length=1),
xt.Drift(_buffer=buf, length=1.),)
+ line_aper_1.elements))
line.build_tracker(_buffer=buf)
num_elements = len(line.element_names)
# Generate test particles
particles = xt.Particles(_context=ctx,
px=np.random.uniform(-0.01, 0.01, 10000),
py=np.random.uniform(-0.01, 0.01, 10000))
#########
# Track #
#########
line.track(particles)
########################
# Refine loss location #
########################
loss_loc_refinement = xt.LossLocationRefinement(line,
n_theta = 360, # Angular resolution in the polygonal approximation of the aperture
r_max = 0.5, # Maximum transverse aperture in m
dr = 50e-6, # Transverse loss refinement accuracy [m]
ds = 0.1, # Longitudinal loss refinement accuracy [m]
save_refine_lines=True # Diagnostics flag
)
loss_loc_refinement.refine_loss_location(particles)
# Complete source: xtrack/examples/collimation/001_loss_location_refinement.py
Generated transition between the defined apertures. Red dots represent the location of the particle-loss events. See the code generating the image.
Beam interaction (generation of secondary particles)
Xtrack includes an interface to ease the modeling of beam-matter interaction (collimators, beam-gas, collisions with another beam), including the loss of the impacting particles and the production of secondary particles, which need to be tracked together with the surviving beam. Such interface can be used to create a link with other programs for the modeling of these effects, e.g. GEANT, FLUKA, K2, GuineaPig.
The interaction is defined as an object that provides a .interact(particles)
method, which sets to zero or negative the state flag for the particles that are lost and
returns a dictionary with the coordinates of the secondary particles that are
emitted. The interaction process is embedded in one or multiple
xtrack.BeamInteraction beam elements that can be included in Xtrack line.
This is illustrated by the following example:
import numpy as np
import xtrack as xt
######################################
# Create a dummy collimation process #
######################################
class DummyInteractionProcess:
'''
I kill some particles. I kick some others by an given angle
and I generate some secondaries with the opposite angles.
'''
def __init__(self, fraction_lost, fraction_secondary, length, kick_x):
self.fraction_lost = fraction_lost
self.fraction_secondary = fraction_secondary
self.kick_x = kick_x
self.length = length
self.drift = xt.Drift(length=self.length)
def interact(self, particles):
self.drift.track(particles)
n_part = particles._num_active_particles
# Kill some particles
mask_kill = np.random.uniform(size=n_part) < self.fraction_lost
particles.state[:n_part][mask_kill] = -1 # special flag`
# Generate some more particles
mask_secondary = np.random.uniform(size=n_part) < self.fraction_secondary
n_products = np.sum(mask_secondary)
if n_products>0:
products = {
's': particles.s[:n_part][mask_secondary],
'x': particles.x[:n_part][mask_secondary],
'px': particles.px[:n_part][mask_secondary] + self.kick_x,
'y': particles.y[:n_part][mask_secondary],
'py': particles.py[:n_part][mask_secondary],
'zeta': particles.zeta[:n_part][mask_secondary],
'delta': particles.delta[:n_part][mask_secondary],
'mass_ratio': particles.x[:n_part][mask_secondary] *0 + 1.,
'charge_ratio': particles.x[:n_part][mask_secondary] *0 + 1.,
'parent_particle_id': particles.particle_id[:n_part][mask_secondary],
'at_element': particles.at_element[:n_part][mask_secondary],
'at_turn': particles.at_turn[:n_part][mask_secondary],
}
else:
products = None
return products
############################################################
# Create a beam interaction from the process defined above #
############################################################
interaction_process=DummyInteractionProcess(length=1., kick_x=4e-3,
fraction_lost=0.0,
fraction_secondary=0.2)
collimator = xt.BeamInteraction(length=interaction_process.length,
interaction_process=interaction_process)
########################################################
# Create a line including the collimator defined above #
########################################################
line = xt.Line(elements=[
xt.Multipole(knl=[0,0]),
xt.LimitEllipse(a=2e-2, b=2e-2),
xt.Drift(length=1.),
xt.Multipole(knl=[0,0]),
xt.LimitEllipse(a=2e-2, b=2e-2),
xt.Drift(length=2.),
collimator,
xt.Multipole(knl=[0,0]),
xt.LimitEllipse(a=2e-2, b=2e-2),
xt.Drift(length=10.),
xt.LimitEllipse(a=2e-2, b=2e-2),
xt.Drift(length=10.),
])
#################
# Build tracker #
#################
line.build_tracker()
line.config.XTRACK_GLOBAL_XY_LIMIT = 1e3
##########################
# Build particles object #
##########################
# We prepare empty slots to store the product particles that will be
# generated during the tracking.
particles = xt.Particles(
_capacity=200000,
x=np.zeros(100000))
#########
# Track #
#########
line.track(particles)
############################
# Loss location refinement #
############################
loss_loc_refinement = xt.LossLocationRefinement(line,
n_theta = 360,
r_max = 0.5, # m
dr = 50e-6,
ds = 0.05,
save_refine_lines=True)
loss_loc_refinement.refine_loss_location(particles)
# Complete source: xtrack/examples/collimation/003_all_together.py