Notebook

This notebook was created by Sergey Tomin (sergey.tomin@desy.de). Source and license info is on GitHub. April 2020. Introduced a few little examples with new features, June 2023, Sergey

An Introduction to Ocelot¶

Ocelot is a multiphysics simulation toolkit designed for studying Free Electron Lasers (FEL) and storage ring-based light sources. Implemented in Python, Ocelot caters to researchers seeking the flexibility provided by high-level languages like Matlab and Python. Its core principle revolves around scripting beam physics simulations in Python, utilizing Ocelot's modules and extensive collection of Python libraries.

Users developing high-level control applications can accelerate development by using physics models from Ocelot and Python graphics libraries such as PyQt and PyQtGraph to create a GUI.

Developing machine learning (ML) applications for accelerators can also benefit from using Ocelot, as many popular ML frameworks are written in Python. Ocelot provides a seamless connection between physics and ML methods, making it easier to integrate physical accelerator simulators with machine learning algorithms.

Ocelot main modules:¶

Charged particle beam dynamics module (CPBD)
- optics
- tracking
- matching
- collective effects (description can be found here and here)
  - Space Charge (3D Laplace solver)
  - CSR (Coherent Synchrotron Radiation) (1D model with arbitrary number of dipoles).
  - Wakefields (Taylor expansion up to second order for arbitrary geometry).
- MOGA (Multi Objective Genetics Algorithm) ref.
Native module for spontaneous radiation calculation (some details can be found here and here)
FEL calculations: interface to GENESIS and pre/post-processing
Modules for online beam control and online optimization of accelerator performances. ref1, ref2, ref3, ref4.
- This module is being developed in collaboration with SLAC. The module has been migrated to a separate repository (in ocelot-collab organization) for ease of collaborative development.

Ocelot extensively uses Python's NumPy (Numerical Python) and SciPy (Scientific Python) libraries, which enable efficient in-core numerical and scientific computation within Python and give you access to various mathematical and optimization techniques and algorithms. To produce high quality figures Python's matplotlib library is used.

It is an open source project and it is being developed by physicists from The European XFEL, DESY (Germany), NRC Kurchatov Institute (Russia).

We still have no documentation but you can find a lot of examples in /demos/ folder and jupyter tutorials

Ocelot installation¶

Requirements¶

Python 3.6 - 3.10
numpy version 1.8 or higher: http://www.numpy.org/
scipy version 0.15 or higher: http://www.scipy.org/
matplotlib version 1.5 or higher: http://matplotlib.org/
h5py version 3.10 or higher, https://www.h5py.org

Optional, but highly recommended for speeding up calculations

numexpr (version 2.6.1 or higher)
pyfftw (version 0.10 or higher)
numba

Orbit Correction module is required

pandas

Installation¶

GitHub (for advanced python users)¶

Clone OCELOT from GitHub:

$ git clone https://github.com/ocelot-collab/ocelot.git

or download last release zip file. Now you can install OCELOT from the source:

$ python setup.py install

Anaconda Cloud (recommended)¶

The easiest way to install OCELOT is to use Anaconda cloud. In that case use command:

$ conda install -c ocelot-collab ocelot

PythonPath¶

Another way is download ocelot from GitHub

you have to download from GitHub zip file.
Unzip ocelot-master.zip to your working folder /your_working_dir/.
Add ../your_working_dir/ocelot-master to PYTHONPATH
- Windows 7: go to Control Panel -> System and Security -> System -> Advance System Settings -> Environment Variables.
and in User variables add /your_working_dir/ocelot-master/ to PYTHONPATH. If variable PYTHONPATH does not exist, create it

Variable name: PYTHONPATH

Variable value: ../your_working_dir/ocelot-master/
- Linux:
```
$ export PYTHONPATH=/your_working_dir/ocelot-master:$PYTHONPATH
```

In [1]:

from IPython.display import Image
# Image(filename='gui_example.png')

Tutorials¶

Beam dynamics¶

Introduction. Tutorial N1. Linear optics.. Web version.
- Linear optics. Double Bend Achromat (DBA). Simple example of usage OCELOT functions to get periodic solution for a storage ring cell.
Tutorial N2. Tracking.. Web version.
- Linear optics of the European XFEL Injector.
- Tracking. First and second order.
- Artificial beam matching - BeamTransform
Tutorial N3. Space Charge.. Web version.
- Tracking through RF cavities with SC effects and RF focusing.
Tutorial N4. Wakefields.. Web version.
- Tracking through corrugated structure (energy chirper) with Wakefields
Tutorial N5. CSR.. Web version.
- Tracking trough bunch compressor with CSR effect.
Tutorial N6. RF Coupler Kick.. Web version.
- Coupler Kick. Example of RF coupler kick influence on trajectory and optics.
Tutorial N7. Lattice design.. Web version.
- Lattice design, twiss matching, twiss backtracking
Tutorial N8. Physics process addition. Laser heater. Web version.
- Theory of Laser Heater, implementation of new Physics Process, track particles w/o laser heater effect.
Tutorial N9. Simple accelerator based THz source. Web version.
- A simple accelerator with the electron beam formation system and an undulator to generate THz radiation.
Tutorial N10. Simple accelerator based THz source. Web verstion
- In this tutorial, a few examples for tracking with parallel-plate corrugated structures are shown. The wakefields model are based on analytical wakefield formulas for flat corrugated structures.

Photon field simulation¶

PFS tutorial N1. Synchrotron radiation module. Web version.
- Simple examples how to calculate synchrotron radiation with OCELOT Synchrotron Radiation Module.
PFS tutorial N2. Coherent radiation module and RadiationField object. Web version.
PFS tutorial N3. Reflection from imperfect highly polished mirror. Web version.
PFS tutorial N4. Converting synchrotron radiation Screen object to RadiationField object for viewing and propagation. Web version.
PFS tutorial N5: SASE estimation and imitation. Web version.

Appendixes¶

Undulator matching. Web version.
- brief theory and example in OCELOT
Some useful OCELOT functions. Web version
- Aperture
- Losses along accelerator lattice
- RK tracking
- Dump the beam distribution at a specific location of the lattice
- Energy jitter. Or simulation of the jitter in the RF parameters.
- Get Twiss paremeters from the beam slice
Example of an accelerator section optimization. Web version
- A simple demo of accelerator section optimization with a standard scipy numerical optimization method.

Checking your installation¶

You can run the following code to check the versions of the packages on your system:

(in IPython notebook, press shift and return together to execute the contents of a cell)

In [1]:

import IPython
print('IPython:', IPython.__version__)

import numpy
print('numpy:', numpy.__version__)

import scipy
print('scipy:', scipy.__version__)

import matplotlib
print('matplotlib:', matplotlib.__version__)

import ocelot
print('ocelot:', ocelot.__version__)

IPython: 8.20.0
numpy: 1.26.3
scipy: 1.11.4
matplotlib: 3.8.2
initializing ocelot...
ocelot: 24.03.0

Tutorial N1. Double Bend Achromat.¶

We designed a simple lattice to demonstrate the basic concepts and syntax of the optics functions calculation. Also, we chose DBA to demonstrate the periodic solution for the optical functions calculation.

In [2]:

from __future__ import print_function

# the output of plotting commands is displayed inline within frontends, 
# directly below the code cell that produced it
%matplotlib inline

# import from Ocelot main modules and functions
from ocelot import *

# import from Ocelot graphical modules
from ocelot.gui.accelerator import *

Creating lattice¶

Ocelot has following elements: Drift, Quadrupole, Sextupole, Octupole, Bend, SBend, RBend, Edge, Multipole, Hcor, Vcor, Solenoid, Cavity, Monitor, Marker, Undulator.

In [3]:

# defining of the drifts
D1 = Drift(l=2.)
D2 = Drift(l=0.6)
D3 = Drift(l=0.3)
D4 = Drift(l=0.7)
D5 = Drift(l=0.9)
D6 = Drift(l=0.2)

# defining of the quads
Q1 = Quadrupole(l=0.4, k1=-1.3)
Q2 = Quadrupole(l=0.8, k1=1.4)
Q3 = Quadrupole(l=0.4, k1=-1.7)
Q4 = Quadrupole(l=0.5, k1=1.3)

# defining of the bending magnet
B = Bend(l=2.7, k1=-.06, angle=2*pi/16., e1=pi/16., e2=pi/16.)

# defining of the sextupoles
SF = Sextupole(l=0.01, k2=1.5) #random value
SD = Sextupole(l=0.01, k2=-1.5) #random value

# cell creating
cell = (D1, Q1, D2, Q2, D3, Q3, D4, B, D5, SD, D5, SF, D6, Q4, D6,
        SF, D5, SD, D5, B, D4, Q3, D3, Q2, D2, Q1, D1)

In [4]:

cell

Out[4]:

(<Drift: name=ID_41638795_ at 0x2815f6aa0>,
 <Quadrupole: name=ID_48212357_ at 0x2815f5d50>,
 <Drift: name=ID_92167354_ at 0x2815f6260>,
 <Quadrupole: name=ID_62144763_ at 0x2815f5de0>,
 <Drift: name=ID_65974384_ at 0x2815f69e0>,
 <Quadrupole: name=ID_10916876_ at 0x2815f5db0>,
 <Drift: name=ID_95483235_ at 0x2815f4850>,
 <Bend: name=ID_90868229_ at 0x2815f5f60>,
 <Drift: name=ID_55934688_ at 0x2815f4ac0>,
 <Sextupole: name=ID_22464080_ at 0x11818fee0>,
 <Drift: name=ID_55934688_ at 0x2815f4ac0>,
 <Sextupole: name=ID_30803941_ at 0x2815f5ff0>,
 <Drift: name=ID_49578218_ at 0x2815f4eb0>,
 <Quadrupole: name=ID_27836471_ at 0x2815f5ed0>,
 <Drift: name=ID_49578218_ at 0x2815f4eb0>,
 <Sextupole: name=ID_30803941_ at 0x2815f5ff0>,
 <Drift: name=ID_55934688_ at 0x2815f4ac0>,
 <Sextupole: name=ID_22464080_ at 0x11818fee0>,
 <Drift: name=ID_55934688_ at 0x2815f4ac0>,
 <Bend: name=ID_90868229_ at 0x2815f5f60>,
 <Drift: name=ID_95483235_ at 0x2815f4850>,
 <Quadrupole: name=ID_10916876_ at 0x2815f5db0>,
 <Drift: name=ID_65974384_ at 0x2815f69e0>,
 <Quadrupole: name=ID_62144763_ at 0x2815f5de0>,
 <Drift: name=ID_92167354_ at 0x2815f6260>,
 <Quadrupole: name=ID_48212357_ at 0x2815f5d50>,
 <Drift: name=ID_41638795_ at 0x2815f6aa0>)

hint: to see a simple description of the function put cursor inside () and press Shift-Tab* or you can type sign ? before function. To extend dialog window press + * Also, one can get more info about element just using print(element)

In [5]:

# all infro about an element can be seen with 
print(B)

Bend(l=2.70000, angle=3.926991e-01, e1=1.963495e-01, e2=1.963495e-01, eid="ID_90868229_")

The cell is a list of the simple objects which contain a physical information of lattice elements such as length, strength, voltage and so on. In order to create a transport map for every element and bind it with lattice object we have to create new Ocelot object - MagneticLattice() which makes these things automatically.

MagneticLattice(sequence, start=None, stop=None, method={"global": TransferMap}):

sequence - list of the elements,

other parameters we will consider in tutorial N2.

Note, in the current version of OCELOT, transfer map belongs to element. See example

In [6]:

# R matrix can be printed for any particular element.
print(Q1.R(energy=0))

[array([[ 1.10581521,  0.4140116 ,  0.        ,  0.        ,  0.        ,
         0.        ],
       [ 0.53821508,  1.10581521,  0.        ,  0.        ,  0.        ,
         0.        ],
       [ 0.        ,  0.        ,  0.89779021,  0.38627683,  0.        ,
         0.        ],
       [ 0.        ,  0.        , -0.50215988,  0.89779021,  0.        ,
         0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  1.        ,
         0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         1.        ]])]

In [7]:

# or you can directly get transfer maps 
Q2.tms

Out[7]:

[<ocelot.cpbd.transformations.transfer_map.TransferMap at 0x2815f5ea0>]

In [8]:

lat = MagneticLattice(cell)

# to see total lenth of the lattice 
print("length of the cell: ", lat.totalLen, "m")

# or, for example, you can get R matrix for whole lattice

B, R, T = lat.transfer_maps(energy=0)
print(R)

length of the cell:  20.34 m
[[ 0.68401288  0.38454837  0.          0.          0.          0.05268746]
 [-1.38376969  0.68401288  0.          0.          0.          0.23072876]
 [ 0.          0.          0.81775255 -0.29733817  0.          0.        ]
 [ 0.          0.          1.11415489  0.81775255  0.          0.        ]
 [ 0.23072876  0.05268746  0.          0.          1.          0.02228572]
 [ 0.          0.          0.          0.          0.          1.        ]]

Optical function calculation¶

Uses:

twiss() function and,
Twiss() object contains twiss parameters and other information at one certain position (s) of lattice

To calculate twiss parameters you have to run twiss(lattice, tws0=None, nPoints=None) function. If you want to get a periodic solution leave tws0 by default.

You can change the number of points over the cell, If nPoints=None, then twiss parameters are calculated at the end of each element. twiss() function returns list of Twiss() objects.

You will see the Twiss object contains more information than just twiss parameters.¶

In [9]:

tws = twiss(lat, nPoints=1000)

# to see twiss paraments at the begining of the cell, uncomment next line
# print(tws[0])
print("length = ", len(tws))
# to see twiss paraments at the end of the cell, uncomment next line
print(tws[-1])

length =  1000
emit_x  = 0.0
emit_y  = 0.0
beta_x  = 0.5271613695963895
beta_y  = 0.5165977895295946
alpha_x = -4.440892098500626e-16
alpha_y = 6.661338147750939e-15
gamma_x = 1.8969523521149319
gamma_y = 1.9357419258618653
Dx      = 0.16673927708143915
Dy      = 0.0
Dxp     = 4.440892098500626e-16
Dyp     = 0.0
mux     = 7.100731992120578
muy     = 5.669884351617213
nu_x    = 1.1301165961167512
nu_y    = 0.9023901213192655
E       = 0.0
s        = 20.34

In [10]:

# plot optical functions.
plot_opt_func(lat, tws, top_plot = ["Dx", "Dy"], legend=False, font_size=10)
plt.show()

# you also can use standard matplotlib functions for plotting
#s = [tw.s for tw in tws]
#bx = [tw.beta_x for tw in tws]
#plt.plot(s, bx)
#plt.show()

In [11]:

# you can play with quadrupole strength and try to make achromat
Q4.k1 = 1.18

# to make achromat uncomment next line
# Q4.k1 =  1.18543769836
# To use matching function, please see ocelot/demos/ebeam/dba.py 

# updating transfer maps after changing element parameters. 
#lat.update_transfer_maps() - not needed anymore

# recalculate twiss parameters. Argument nPoints is None by default - Twiss is calculating at the end of each element. 
# If you want smooth twiss functions you can set number of points. 
tws=twiss(lat, nPoints=1000)

plot_opt_func(lat, tws, legend=False)
plt.show()