Error-prone code because accelerator element properties are energy-dependent

[RemiLehe]

Overview

In Cheetah, when the user creates an accelerator element (e.g., a quadrupole, a dipole), they pass parameters that depend on the element itself, but also on the energy of the beam that will go through that element. (For instance, k1 for a quadrupole and angle for a dipole depend on the beam energy.)

This can be error prone, because the user might then use this accelerator element with a beam that is not at the energy for which the accelerator element was setup.

Example

Consider the following example, where a "naive" user would try to see the influence of a change in beam energy on the final beam size after a quad+drift.

from cheetah import ParticleBeam, Quadrupole, Drift, Segment
import torch

# Setup a quadrupole followed by a drift
# The k1 parameter (100.) of the quadrupole is the nominal one for 1 MeV
# but the syntax does not make this obvious, and a user might forget this.
segment = Segment(
    elements=[
        Quadrupole(length=torch.tensor(2.e-1), k1=torch.tensor(100.), tracking_method="bmadx"),
        Drift(length=torch.tensor(2.e-1))
    ]
)

# Try for a 1 MeV beam (correct energy for this segment)
beam_energy = torch.tensor(1e6)
incoming = ParticleBeam.from_twiss(
        beta_x=torch.tensor(0.1),
        alpha_x=torch.tensor(0.0),
        emittance_x=torch.tensor(1.e-7),
        beta_y=torch.tensor(0.1),
        alpha_y=torch.tensor(0.0),
        emittance_y=torch.tensor(1.e-7),
        sigma_tau=torch.tensor(1e-4), # in m
        sigma_p=torch.tensor(1e-3), # dimensionless
        energy=beam_energy, # in eV
        total_charge=torch.tensor(100e-12), # in C
)
print( 'Beam size at 1 MeV:', segment.track(incoming).sigma_y )


# Try for a 0.7 MeV beam, naively (incorrect energy for this segment)
beam_energy = torch.tensor(0.7e6)
incoming = ParticleBeam.from_twiss(
        beta_x=torch.tensor(0.1),
        alpha_x=torch.tensor(0.0),
        emittance_x=torch.tensor(1.e-7),
        beta_y=torch.tensor(0.1),
        alpha_y=torch.tensor(0.0),
        emittance_y=torch.tensor(1.e-7),
        sigma_tau=torch.tensor(1e-4), # in m
        sigma_p=torch.tensor(1e-3), # dimensionless
        energy=beam_energy, # in eV
        total_charge=torch.tensor(100e-12), # in C
)
print( 'Beam size at 0.7 MeV (wrong):', segment.track(incoming).sigma_y )

# Try for a 0.7 MeV beam, in the correct way
# - Initialize a beam at the energy (1 MeV), for which the quadrupole was setup
beam_energy = torch.tensor(1e6)
incoming = ParticleBeam.from_twiss(
        beta_x=torch.tensor(0.1),
        alpha_x=torch.tensor(0.0),
        emittance_x=torch.tensor(1.e-7),
        beta_y=torch.tensor(0.1),
        alpha_y=torch.tensor(0.0),
        emittance_y=torch.tensor(1.e-7),
        sigma_tau=torch.tensor(1e-4), # in m
        sigma_p=torch.tensor(1e-3), # dimensionless
        energy=beam_energy, # in eV
        total_charge=torch.tensor(100e-12), # in C
)
# - Now change the beam to give it a 0.7 MeV energy
incoming.p -= 0.3
# - Track the beam through the beamline
print( 'Beam size at 0.7 MeV (correct):', segment.track(incoming).sigma_y )

With this example, one gets:

Beam size at 1 MeV: tensor(0.0016)
Beam size at 0.7 MeV (wrong): tensor(0.0016)
Beam size at 0.7 MeV (correct): tensor(0.0026)

In other words, with the first (naive) method, the user would mistakenly think that the beam size does not depend on the beam energy. The second method gives the correct answer, but is convoluted: the user first needs to create a 1 MeV beam (which corresponds to the nominal energy of the element) but then modifies it by hand.

Proposed solution

Define Cheetah element by their intrinsic, beam-independent properties (e.g. the user passes B_gradient instead of k1 for a quad, B instead of angle for a dipole -- and these properties would be the ones that are stored internally in Cheetah).

[cr-xu]

Programs like mad-x and elegant use the normalized field strengths $k_1, k_2,...$, while codes like bmad use field gradient $B$ and derive the $k_1$ from the energy. I think using normalized field strengths in magnets is more convenient if the focus is the lattice design. We can definitely add convenient methods to get the actual field gradient.

In the end isn't it actually a design choice whether to take field gradient or normalized strengths as the independent properties?

[RemiLehe]

Another solution discussed offline would be to pass the design energy when creating an element ; we could then add a check in track that the beam reference energy is the same as the design energy.