<jupyter_start><jupyter_text>Introduction To SixTrackLib Implementation And DesignIn this section, we cover the following topics:- Modelling the Particle State- Lattice And Beam Elements- Complementary and Required Externall Libraries &amp; Modules        - pysixtrack    - sixtracktools    - CObjects Buffer- Logistics of Tracking: Tracking Modes, TrackJob      Modeling the Particle State<jupyter_code>p = st.Particles(num_particles=1, p0c=6.5e12, q0=1)
print( f"""Particle-State:
Six main degrees of freedom: x     = {p.x} m\r\n     
                             y     = {p.y} m\r\n
                             zeta  = {p.zeta} m\r\n
                             px    = {p.px} rad ( px = Px / P0 )\r\n
                             py    = {p.py} rad ( py = Py / P0 )\r\n
                             delta = {p.delta}     ( δ  = ( P - P0 ) / P0 )\r\n
""")
print(f"""
4 Logical coordinates: state       = {p.state} ( 1 == active, 0 == lost )\r\n
                       at_element  = {p.at_element}\r\n
                       at_turn     = {p.at_turn}\r\n
                       particle_id = {p.particle_id}\r\n""")
print(f"""
5 Attributes of the ref particle: q0     = {p.q0} x proton charge\r\n
                                  mass0  = {p.mass0} eV/c²\r\n
                                  beta0  = {p.beta0} ( β0 = v0/c )\r\n
                                  gamma0 = {p.gamma0} ( γ0 = 1 / √( 1 - β0²) )\r\n
                                  p0c    = {p.p0c} eV  ( p0c = P0 * c )\r\n 
""")
print(f"""
6 auxilliariy coordinates: s            = {p.s} m   ( s ~ distance from begin of lattice )\r\n
                           psigma       = {p.psigma}     ( psigma = (E-E0)/(beta0*P0*c) )\r\n
                           rpp          = {p.rpp}     ( rpp = P0 / P )\r\n
                           rvv          = {p.rvv}     ( rvv = β / β0 = (v/c) / β0 )\r\n
                           charge_ratio = {p.charge_ratio}     ( charge_ratio = q/q0 )\r\n
                           chi          = {p.chi}     ( χ = ( q/q0 ) / ( m/mass0 ) )\r\n""")<jupyter_output>6 auxilliariy coordinates: s            = [0.] m   ( s ~ distance from begin of lattice )

                           psigma       = [0.]     ( psigma = (E-E0)/(beta0*P0*c) )

                           rpp          = [1.]     ( rpp = P0 / P )

                           rvv          = [1.]     ( rvv = β / β0 = (v/c) / β0 )

                           charge_ratio = [1.]     ( charge_ratio = q/q0 )

                           chi          = [1.]     ( χ = ( q/q0 ) / ( m/mass0 ) )<jupyter_text>Lattice &amp; Beam ElementsThere are several ways to get a lattice1. manually, element by element2. load from a binary dump3. import from pysixtrack4. import from MAD-X (via pysixtrack)5. import from SixTrack (via pysixtrack)  1. Build Lattice Manually<jupyter_code># We can use a st.Elements() class to organize individual beam-elements into a lattice

lattice = st.Elements()
drift = lattice.Drift(lenth=0.2) # length in [m]
quad  = lattice.Multipole( knl=[0.0, 0.01] ) # ksl ... skew multipole parameters

lattice.get_elements()
# Dump the lattice to a binary file
lattice.to_file( "./demo_lattice.bin" )
del lattice # make sure we start from scratch<jupyter_output><empty_output><jupyter_text>2. Load From Binary Dump<jupyter_code>lattice = st.Elements().fromfile( "./demo_lattice.bin" )

lattice.get_elements()<jupyter_output><empty_output><jupyter_text>3. Import From PySixtrack Question: What is pysixtrack ? - pysixtrack is a minimal &amp; straight-forward particle tracking implementation written purely in Python 3- https://github.com/SixTrack/pysixtrack- Independent of SixTrackLib, easy to understand &amp; extend - "Playground", Area for prototyping, testing and debugging new physics<jupyter_code>import pysixtrack as py6tr

# We can use any iterable object to store a lattice:
seq = [ py6tr.elements.Drift(length=0.2), py6tr.elements.Multipole(knl=[0.0, 0.01]) ]

print( seq )
# Use pysixtrack.Line object to store the sequence of beam-elements
other_lattice = py6tr.Line( seq )
    
# ensure we start from scratch with lattice
del lattice

# import other_lattice as a SixTrackLib lattice:
lattice = st.Elements().from_line( other_lattice )

print( lattice.get_elements() )<jupyter_output>[<Drift at 128
  length:0.2
>, <Multipole at 136
  order:1
  length:0.0
  hxl:0.0
  hyl:0.0
  bal:[0.   0.   0.01 0.  ]
>]<jupyter_text>4. Import From MAD-X Using pysixtrack (&amp; cpymad)The idea is to keep SixTrackLib as minimal as possible.&rightarrow; No I/O and Import/Export Helpers in SixTrackLib, use pysixtrack as an intermediate layerNote: For users of the C/C++ API of SixTrackLib, use binary dumps to import lattices and particle state from MAD-X, SixTrack, ...Note: Import from MAD-X requires the cpymad cython bindings:https://github.com/hibtc/cpymad<jupyter_code># requires cpymad -> https://github.com/hibtc/cpymad
from cpymad.madx import Madx
import sixtracklib as st 
import pysixtrack as py6tr

from scipy.constants import e, m_p, c
import numpy as np 

# Note: pysixtrack / SixTrackLib and MAD-X use different default units for energies!
p0c = 4.7e9 # p0c = P0 * c ; [p0c] = 1 eV
Etot_in_GeV = np.sqrt( p0c * p0c + ( m_p / e ) ** 2 * c ** 4 ) * 1e-9 # [Etot] = 1 GeV

mad = Madx(stdout=False)
mad.call( file="./demo_lattice.madx" )
mad.command.beam(particle='proton', energy=str(Etot_in_GeV))
mad.use(sequence="DEMO_LATTICE")

# Use the from_madx_sequence() method of pysixtrack's Line to import DEMO_LATTICE 
# Also, use the remove_zero_length_drifts() and merge_consequentive_drifts() 
# helpers to optimize the lattice for SixTrackLib's preferred way of torage :

imported_lattice = py6tr.Line.from_madx_sequence( 
    mad.sequence.DEMO_LATTICE, ).remove_zero_length_drifts(
    inplace=True).merge_consecutive_drifts(inplace=True)
  
# As before, use imported_lattice to build the SixTrackLib lattice
# that we are actually interested in:
lattice = st.Elements().from_line( imported_lattice )
print( lattice.get_elements() )<jupyter_output>[<Drift at 128
  length:0.2
>, <Multipole at 136
  order:1
  length:0.0
  hxl:0.0
  hyl:0.0
  bal:[0.   0.   0.01 0.  ]
>]<jupyter_text>5. Import From SixTrack Using pysixtrack (&amp; sixtracktools)Similar approach, but this time we are importing a more sophisticated Lattice:Note: sixtracktools is a helper library which runs the SixTrack binary and interprets the output, allowing import data into pysixtrackhttps://github.com/SixTrack/sixtracktools<jupyter_code>import sixtracktools 
import pysixtrack as py6tr
import sixtracklib as st

# SixTrack input files + helper script to run the SixTrack binary are 
# in subdirectory here:

!ls -al ./sixtrack_lhc_no_bb_example/
six = sixtracktools.SixInput("./sixtrack_lhc_no_bb_example" )
import_lattice = py6tr.Line.from_sixinput( six ).remove_zero_length_drifts(
    inplace=True).merge_consecutive_drifts(inplace=True)

print( f"import_lattice has {len(import_lattice)} elements" )

# Same procedure as usual -> convert import_lattice to a SixTrackLattice
lattice = st.Elements().from_line( import_lattice )

# Create a binary dump for the machine description:
lattice.to_file( "./lhc_no_bb_lattice.bin" )

# Verify the size of the binary dump:
!ls -alh ./lhc_no_bb_lattice.bin<jupyter_output>import_lattice has 18403 elements
-rw-rw-r-- 1 <NAME> 3,8M Jun 16 23:19 ./lhc_no_bb_lattice.bin<jupyter_text>Tracking Examples  Simple Tracking Example (CPU)<jupyter_code># Create an initial particle distribution:

beam = st.ParticlesSet()
p = beam.Particles(num_particles=10, p0c=6.5e12)
p.x[:] = np.linspace(-1e-6, +1e-6, p.num_particles)

if p.num_particles <= 16:
    print( f"initial transversal displacement for particles: {p.x}\r\n" )

# Load the lattice from the binary dump we crated earlier
lattice = st.Elements().fromfile("./lhc_no_bb_lattice.bin")
print( f"number of elements in lattice: {lattice.cbuffer.n_objects}") 

# What's cbuffer -> Cf. BE Seminar talk for details!

# Setup a track-job instance:
job = st.TrackJob( lattice, beam )

# Print particle state before tracking:
if p.num_particles <= 16:    
    print( f"at_element before tracking: {p.at_element}" )
    print( f"at_turn    before tracking: {p.at_turn}" )
    print( f"state      before tracking: {p.state}" )
    print( f"x          before tracking: {p.x}" )
# Track particles <b>until</b> they are in turn 100
job.track_until( 100 )

# Print particle state after tracking for 100 turns:
if p.num_particles <= 16:    
    print( f"at_element after tracking: {p.at_element}" )
    print( f"at_turn    after tracking: {p.at_turn}" )
    print( f"state      after tracking: {p.state}" )
    print( f"x          after tracking: {p.x}" )
    
# Note: the command is called <tt>track_until</tt> - if we call it again,
# it will have no effect because all particles are alreaedy at turn 100!
# -> we would have to call <tt>track_until( 200 )</tt> to get the desired effect
# Let's Loose a particle
p.state[0] = 0
if p.num_particles <= 16:        
    print( f"state after manually loosing a particle: {p.state}\r\n" )
    
# Track until turn 200 and verify the result:
job.track_until( 200 )

if p.num_particles <= 16:    
    print( f"at_element after 200 turns : {p.at_element}" )
    print( f"at_turn    after 200 turns : {p.at_turn}" )
    print( f"state      after 200 turns : {p.state}" )
    print( f"x          after 200 turns : {p.x}" )<jupyter_output>state after manually loosing a particle: [0 1 1 1 1 1 1 1 1 1]

at_element after 200 turns : [0 0 0 0 0 0 0 0 0 0]
at_turn    after 200 turns : [100 200 200 200 200 200 200 200 200 200]
state      after 200 turns : [0 1 1 1 1 1 1 1 1 1]
x          after 200 turns : [-9.99845051e-07 -7.77491911e-07 -5.55303462e-07 -3.33123094e-07
 -1.10949224e-07  1.11219787e-07  3.33385573e-07  5.55549726e-07
  7.77713872e-07  9.99879665e-07]<jupyter_text>Simple Tracking Example (OpenCL, GPU)First, check whether we have any OpenCL devices and whether SixTrackLib has been compiled with OpenCL support:<jupyter_code>print( f"SixTrackLib has OpenCL support enabled: {st.config.SIXTRACKLIB_MODULES['opencl']}")
!clinfo -l 
# Re-create the initial state:
beam = st.ParticlesSet()
p = beam.Particles(num_particles=10, p0c=6.5e12)
p.x[:] = np.linspace(-1e-6, +1e-6, p.num_particles)
lattice = st.Elements().fromfile("./lhc_no_bb_lattice.bin")

# Again, create a TrackJob. But this time, we pass the "device" string
opencl_job = st.TrackJob( lattice, beam, device="opencl:1.0" )

print( f"track job instance has architecture {opencl_job.arch_str}" )
# Again: track until turn 100
opencl_job.track_until( 100 )

# Print particle state after tracking for 100 turns:
if p.num_particles <= 16:    
    print( f"at_element after tracking: {p.at_element}" )
    print( f"at_turn    after tracking: {p.at_turn}" )
    print( f"state      after tracking: {p.state}" )
    print( f"x          after tracking: {p.x}" )
    
# Spoiler: this does not seem to work. Why?<jupyter_output>at_element after tracking: [0 0 0 0 0 0 0 0 0 0]
at_turn    after tracking: [0 0 0 0 0 0 0 0 0 0]
state      after tracking: [1 1 1 1 1 1 1 1 1 1]
x          after tracking: [-1.00000000e-06 -7.77777778e-07 -5.55555556e-07 -3.33333333e-07
 -1.11111111e-07  1.11111111e-07  3.33333333e-07  5.55555556e-07
  7.77777778e-07  1.00000000e-06]<jupyter_text>Remember From Introduction To GPU Programming Talk:<jupyter_code># Copying the data back from the device to the host can be a costly operation
# Thus, it is not done automatically -> we have to "collect" the results
opencl_job.collect_particles()

# If we print the the particle data now, we should have the exprected particle state
if p.num_particles <= 16:    
    print( f"at_element after tracking: {p.at_element}" )
    print( f"at_turn    after tracking: {p.at_turn}" )
    print( f"state      after tracking: {p.state}" )
    print( f"x          after tracking: {p.x}" )<jupyter_output>at_element after tracking: [0 0 0 0 0 0 0 0 0 0]
at_turn    after tracking: [100 100 100 100 100 100 100 100 100 100]
state      after tracking: [1 1 1 1 1 1 1 1 1 1]
x          after tracking: [-9.99845052e-07 -7.77634846e-07 -5.55429505e-07 -3.33228214e-07
 -1.11030165e-07  1.11165447e-07  3.33359448e-07  5.55552645e-07
  7.77745826e-07  9.99939829e-07]<jupyter_text>When To Call collect_* And push_*:  Whenever we need the current content of a dataset (e.g. particles, beam elements, ...) on the host, we have to call collect_*  Whenever we want to send the current content of a dataset (again: particles, beam elements, ... ) to the device, we have to call push_* Host To Device: &nbsp;push_*<jupyter_code># Repeat the experiment with the explicitly "lost" particle:
p.state[ 0 ] = 0

if p.num_particles <= 16:        
    print( f"state after manually loosing a particle: {p.state}\r\n" )

# It should not be a big surprise, that we need something equivalent to "collect"
# but working in the other direction, i.e. from Host -> Device
opencl_job.push_particles()
    
# Track until turn 200:
opencl_job.track_until( 200 )

# everytime we need the particle state on the host side, we
# have to collect the data:
opencl_job.collect_particles()

# Now, we expect the same output as before:
if p.num_particles <= 16:    
    print( f"at_element after 200 turns : {p.at_element}" )
    print( f"at_turn    after 200 turns : {p.at_turn}" )
    print( f"state      after 200 turns : {p.state}" )
    print( f"x          after 200 turns : {p.x}" )<jupyter_output>state after manually loosing a particle: [0 1 1 1 1 1 1 1 1 1]

at_element after 200 turns : [0 0 0 0 0 0 0 0 0 0]
at_turn    after 200 turns : [100 200 200 200 200 200 200 200 200 200]
state      after 200 turns : [0 1 1 1 1 1 1 1 1 1]
x          after 200 turns : [-9.99845052e-07 -7.77491913e-07 -5.55303460e-07 -3.33123095e-07
 -1.10949223e-07  1.11219786e-07  3.33385572e-07  5.55549725e-07
  7.77713871e-07  9.99879664e-07]<jupyter_text>collect_* and push_* are potentially expensive calls (band-width for transfer, latency, waiting for all running kernels    They contribute to $t_s$!!!    $\rightarrow$With the exception of the initial push_* when setting up the track job, these are not performed automatically!     Calling push_* and collect_* has (almost) no negative run-time-cost effect on a CPU track-job     $\Longrightarrow$ If you call them also with a CPU track-job, your code works on the GPU with just changing the setup line of the track-job!  Generic Tracking Program (Works On CPUs &amp; GPUs With Minimal Changes)<jupyter_code>beam = st.ParticlesSet()
p = beam.Particles(num_particles=10, p0c=6.5e12)
p.x[:] = np.linspace(-1e-6, +1e-6, p.num_particles)
lattice = st.Elements().fromfile("./lhc_no_bb_lattice.bin")

#device=None  # Or: 
device="opencl:0.0" #for GPU

job = st.TrackJob( lattice, beam, device=device )
print( f"Architecture of the track job: {job.arch_str}")

job.track_until( 100 )
job.collect_particles()

p.state[0] = 0 # Mark particle 0 explicitly as lost
job.push_particles()

job.track_until( 200 )
job.collect_particles()

if p.num_particles <= 16:    
    print( f"at_element after 200 turns : {p.at_element}" )
    print( f"at_turn    after 200 turns : {p.at_turn}" )
    print( f"state      after 200 turns : {p.state}" )<jupyter_output>Architecture of the track job: opencl
at_element after 200 turns : [0 0 0 0 0 0 0 0 0 0]
at_turn    after 200 turns : [100 200 200 200 200 200 200 200 200 200]
state      after 200 turns : [0 1 1 1 1 1 1 1 1 1]
{'fraction_non_alphanumeric': 0.10964572019653478, 'fraction_numerical': 0.06633048875096974, 'mean_word_length': 2.9766066838046275, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 53, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 10, 'word_list_counts': {'cursed_substrings.json': 38, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '431869', 'license': 'MIT', 'n_tokens_mistral': 6109, 'n_tokens_neox': 5393, 'n_words': 1730, 'path': '20200617_hss_section_meeting_sixtracklib/notebook/introduction_to_sixtracklib.ipynb', 'repo_name': 'martinschwinzerl/presentations', 'size': '15475'}