Quickstart Guide

Getting Started

Let's start by constructing a simple FODO cell Beamline, which consists of a quadrupole magnet that focuses the beam in the horizontal plane (defocuses in the vertical), followed by a drift (empty space), then a quadrupole that defocuses in the horizontal plane (focuses in the vertical), and finally another drift.

To do this, we will define four LineElements corresponding to each of these objects. The lengths of each object are specified by L (in meters), and the quadrupole strengths can be set using the property Kn1, where n means the "normal" multipole (s would be "skew") and 1 means 1st order multipole (quadrupole).

using Beamlines

qf = Quadrupole(Kn1=0.36, L=0.5)
d = Drift(L=1.2)
qd = Quadrupole(Kn1=-0.36, L=0.5)

fodo = Beamline([qf, d, qd, d])

Beamline:
 species_ref = Inferred
 p_over_q_ref = Inferred

  Index   Name   Kind         s [m] 
  1              Quadrupole   0
  2              Drift        0.5
  3              Quadrupole   1.7
  4              Drift        2.2

As you can see, because we did not specify a reference particle species species_ref, nor a signed reference magnetic rigidity p_over_q_ref, both of those show Inferred. This means that it will infer those values from either a preceeding Beamline, or simply leave it up to a tracking code to decide. One can also specify E_ref or pc_ref instead of p_over_q_ref for convenience. Let's do that now, and specify electrons with a total reference energy of 18 GeV:

fodo = Beamline([qd, d, qd, d], species_ref=Species("electron"), E_ref=18e9)

Beamline:
 species_ref = electron
 E_ref = 1.8e10

  Index   Name   Kind         s [m] 
  1              Quadrupole   0
  2              Drift        0.5
  3              Quadrupole   1.7
  4              Drift        2.2

Beamlines.jl uses the AtomicAndPhysicalConstants.jl package for specifying particle species, and so any species defined by that package may be provided.

The rest of the output looks ok, except for the fact that the name column is empty! This is because we didn't specify a name property for each element. It would often be convenient if we can make the variable symbols (e.g. qf, d, etc.) automatically fill in the name field for each element. We can do exactly this by wrapping the element definitions in a @elements block:

@elements begin
    qf = Quadrupole(Kn1=0.36, L=0.5)
    d = Drift(L=1.2)
    qd = Quadrupole(Kn1=-0.36, L=0.5)
end

fodo = Beamline([qf, d, qd, d], species_ref=Species("electron"), E_ref=18e9)

Beamline:
 species_ref = electron
 E_ref = 1.8e10

  Index   Name   Kind         s [m] 
  1       qf     Quadrupole   0
  2       d      Drift        0.5
  3       qd     Quadrupole   1.7
  4       d      Drift        2.2

Much better!

Python users may use the dict-based naming function elements instead.

@elements ...;
@elements begin
  ...
end

julia> @elements qf = Quadrupole();

julia> qf.name
"qf"

julia> @elements begin
         my_drift = Drift(L=0.5);
         my_sextupole = Sextupole(Kn2L=3.2);
       end;

julia> my_drift.name
"my_drift"

julia> my_sextupole.name
"my_sextupole"

elements(eles_dict)

Deferred Expressions

Earlier we set qf.Kn1 = 0.36, and qd.Kn1 = -0.36. But what if we want to ensure that qd.Kn1 == -qf.Kn1 always? We can bake-in such an interdependence, common in particle accelerator parameters, using a "deferred expression" - an expression where evaluation is postponed until its result is actually needed, rather than immediately when it is defined.

To do this, let's first define a function that returns the current value of -qf.Kn1. We can do this without giving the function any explicit name using lambda/anonymous functions:

lambdafun = () -> -qf.Kn1
println("Before: ", lambdafun())
qf.Kn1 = 0.1
println("After: ", lambdafun())

Before: -0.36
After: -0.1

Here lambdafun takes no arguments (specified by the empty tuple ()) and returns -qf.Kn1. In the context of programming, lambdafun is specifically called a closure, because it "encloses" qf, and at the time of evaluation gets the Kn1 of that enclosed qf and negates its sign.

Now we just wrap this function in Beamlines's DefExpr type, and we can set any LineElement parameter to be such a deferred expression:

qd.Kn1 = DefExpr(lambdafun)
qd.Kn1

-0.1

Now if we change qf.Kn1, evaluation of qd.Kn1 will always be -qf.Kn1:

qf.Kn1 = 0.7
qd.Kn1

-0.7

Deferred expressions can also be manipulated like any other number:

a = 1
da = DefExpr(()->a)
b = 2
db = DefExpr(()->b)
dc = da + db
println(dc())
a = 4
println(dc())
dd = sin(dc)
println(dd())

3
6
-0.27941549819892586

One can really "go crazy" with deferred expressions if they want to. They can be infinitely nested, and you can write any function that the programming language allows, for example file I/O, or even control system gets/puts with a real accelerator for a digital twin.

Parameters

Beamlines.jl supports a continually-growing list of parameters to define accelerator elements. To see a full list of the parameters you can set, look at the docstring for the LineElement type. This can be retrieved in a Julia session using Doc.docs(LineElement).

LineElement

  AlignmentParams         ApertureParams               BMultipoleParams
   x_offset                x1_limit                     KnX
   y_offset                x2_limit                     KsX
   z_offset                y1_limit                     BnX
   x_rot                   y2_limit                     BsX
   y_rot                   aperture_shape               KnXL
   tilt                    aperture_at                  KsXL
                           aperture_shifts_with_body    BnXL
                           aperture_active              BsXL
                                                        tiltX

  BeamlineParams          BendParams                   FourPotentialParams
   beamline                g_ref                        four_potential
   beamline_index          tilt_ref                     four_potential_params
   s                       e1                           four_potential_normalized
   s_downstream            e2
   line                    edge1_int
   branch                  edge2_int
   branch_index

  InitialBeamlineParams   MapParams                    MetaParams
   species_ref             transport_map                alias
   E_ref                   transport_map_params         label
   pc_ref                                               description
   p_over_q_ref
   dE_ref
   dpc_ref
   dp_over_q_ref

  PatchParams             RFParams                     UniversalParams
   dt                      rf_frequency                 kind
   dx                      harmon                       name
   dy                      phi0                         L
   dz                      zero_phase                   tracking_method
   dx_rot                  traveling_wave
   dy_rot                  is_crabcavity
   dz_rot

ele = LineElement()
ele.L = 1   # Length 
ele.Kn1 = 2 # Normal quadrupole strength
ele.rf_frequency = 1e6 

For a more detailed description of each parameter, see the docstrings for individual parameter groups. These are shown in the Parameter Groups section of the documentation.

Multiple LineElements in (multiple) Beamlines

Let's go back to the earlier example,

@elements begin
    qf = Quadrupole(Kn1=0.36, L=0.5)
    d = Drift(L=1.2)
    qd = Quadrupole(Kn1=-0.36, L=0.5)
end

fodo = Beamline([qf, d, qd, d], species_ref=Species("electron"), E_ref=18e9);

Beamline:
 species_ref = electron
 E_ref = 1.8e10

  Index   Name   Kind         s [m] 
  1       qf     Quadrupole   0
  2       d      Drift        0.5
  3       qd     Quadrupole   1.7
  4       d      Drift        2.2

Here, fodo contains two instances of the line element d. Therefore, if the length of d is changed, then both instances of d will see this new, changed length:

d.L = 2.0
println(fodo.line[2].L)
println(fodo.line[4].L)

2.0
2.0

However, both drifts in fodo are unique elements. We can check this using the === operator:

println(fodo.line[2] === fodo.line[4])
println(d === fodo.line[2])
println(d === fodo.line[4])

false
false
false

Under the hood, when an element is placed in a Beamline, a shallow copy of that element is created that points to the "parent" element, from which it inherits its parameters. So, in this above example when the "get" fodo.line[2].L is executed, the code goes to the parent element d and returns d.L. "Sets", such as fodo.line[2].L = 10, will also pass through from the child to the parent:

fodo.line[2].L = 3.0
println(d)
println(fodo.line[4].L)
println(d === fodo.line[4].L)

LineElement:
  UniversalParams
   kind            = "Drift"
   name            = "d"
   L               = 3.0
   tracking_method = SciBmadStandard(
     radiation_damping_on      = false,
     radiation_fluctuations_on = false,
     ibs_damping_on            = false,
     ibs_fluctuations_on       = false,
   )


3.0
false

The only case where a child element can have parameters different from its parent is when a given parameter group is contained within the child. For example, fodo.line[2] and fodo.line[4] both have their own instance of BeamlineParams, from which we can extract things like beamline_index, s, and s_downstream. On the other hand, the parent element d does not have a BeamlineParams.

println("beamline_index:")
println(fodo.line[2].beamline_index)
println(fodo.line[4].beamline_index)

println("s_downstream:")
println(fodo.line[2].s_downstream)
println(fodo.line[4].s_downstream)

# This will error:
try
d.beamline_index
catch err
println(err)
end

beamline_index:
2
4
s_downstream:
3.5
7.0
ErrorException("  Unable to get key beamline_index from LineElement: element is not in a Beamline. \n  If you placed this element in a Beamline, use `findchildren` to find \n  the child instances of this element in a given Beamline.\n")

The parent element can be retrieved using parent:

fodo.line[2].parent

LineElement:
  UniversalParams
   kind            = "Drift"
   name            = "d"
   L               = 3.0
   tracking_method = SciBmadStandard(
     radiation_damping_on      = false,
     radiation_fluctuations_on = false,
     ibs_damping_on            = false,
     ibs_fluctuations_on       = false,
   )

Sometimes it can be a pain to find exactly where in a beamline are the corresponding child elements. As such, the findchildren function is provided for convenience:

children = findchildren(d, fodo);

2-element Vector{LineElement}:
 LineElement:
  BeamlineParams          InheritParams
   beamline_index = 2      parent = LineElement:
   s              = 0.5     UniversalParams
   s_downstream   = 3.5      kind            = "Drift"
                             name            = "d"
                             L               = 3.0
                             tracking_method = SciBmadStandard(
                               radiation_damping_on      = false,
                               radiation_fluctuations_on = false,
                               ibs_damping_on            = false,
                               ibs_fluctuations_on       = false,
                             )




 LineElement:
  BeamlineParams          InheritParams
   beamline_index = 4      parent = LineElement:
   s              = 4.0     UniversalParams
   s_downstream   = 7.0      kind            = "Drift"
                             name            = "d"
                             L               = 3.0
                             tracking_method = SciBmadStandard(
                               radiation_damping_on      = false,
                               radiation_fluctuations_on = false,
                               ibs_damping_on            = false,
                               ibs_fluctuations_on       = false,
                             )

Finally, elements in a beamline allow one to "get" parameters that may only be defined when said element is in a beamline. We showed the s and s_downstream, but another example would be the unnormalized magnetic field, if the normalied magnetic field is stored as an independent variable:

ele = Quadrupole(Kn1=2, L=2)
bl = Beamline([ele], p_over_q_ref=3)
println(bl.line[1].Bn1) # Returns Kn1 * p_over_q_ref = 2 * 3

6.0

The last parameter "set" will always define what the independent variable is. So if we then set the unnormalized quadrupole strength Bn1, that will be the independent variable:

ele.Bn1 = 10
println(bl.line[1].Kn1) # Returns Bn1 / p_over_q_ref = 10 / 3

3.3333333

Now, if we then change the reference energy of the beamline, Bn1 will remain constant but Kn1 will change:

bl.p_over_q_ref = 4
println(bl.line[1].Bn1) # == 10
println(bl.line[1].Kn1) # Now equals 10 / 4

10.0
2.5

Polymorphism/Differentiability

To enable full auto-differentiability of all accelerator parameters, Beamlines.jl is fully polymorphic. Full polymorphism this means that you can set any parameter to be any type that you want. For auto-differentiability, a special number type that propagates the partial derivative(s) with actual value must be used in-place of the regular 64-bit floating point numbers – polymorphism allows that. Differentiable codes in accelerator physics are not new: an early example of such a differentiable code is the Polymorphic Tracking Code (PTC), written in Fortran back in the 90s.

As an example, let's see how to compute the derivative of the total length of the beamline w.r.t. a particular element length. We will use the GTPSA.jl package to do so.

using GTPSA
d1 = Descriptor(1, 1) # 1 variable, 1st order

@elements begin
    qf = Quadrupole(Kn1=0.36, L=0.5)
    d = Drift(L=1.2)
    qd = Quadrupole(Kn1=-0.36, L=0.5)
end

fodo = Beamline([qf, d, qd, d], species_ref=Species("electron"), E_ref=18e9);
println(fodo.line[end].s_downstream)

3.4000000000000004

Now we just need to update L to be a differential-algebra variable,

ΔL = vars(d1)[1] # get the first differential

d.L += ΔL
println(fodo.line[end].s_downstream)

GTPSA.TPS64{GTPSA.Dynamic}:
Descriptor(NV=1, MO=1)
 COEFFICIENT             ORDER   EXPONENTS
 3.4000000000000004E+00    0     0
 2.0000000000000000E+00    1     1

This output shows that the total length of fodo is equal to $3.4 + 2L$, which is exactly what we'd expect given that there are two drifts.

Here we just showed the length, but **any*** accelerator parameters defined in Beamlines.jl can be set to any number type (fully polymorphic), so that derivatives can be computed using any automatic-differentiation package in Julia.

After computing derivatives, e.g. during an optimization, one might want to restore all number types back to their primitive values ( Float64, Float32, etc). This can be done using the scalarize! function:

scalarize!(fodo)
println(fodo.line[end].s_downstream)

3.4000000000000004

scalarize!(ele::LineElement)

scalarize!(bl::Beamline)

scalarize!(branch::Branch)