Minimal working example

Author: ewaagaard

data/acc-models-sps

@@ -0,0 +1,142 @@
+"""
+Minimalistic example to load an xsuite line, generate a particle distribution,
+install space charge, IBS, and track the particles.
+"""
+import xtrack as xt
+import xpart as xp
+import xfields as xf
+import xobjects as xo   
+from records import Records
+import time
+
+# Tracking parameters and context
+n_turns = 20
+n_particles = 20
+context = xo.ContextCpu(omp_num_threads='auto')
+
+# Beam parameters
+exn = 1.0e-6  # horizontal emittance in m
+eyn = 1.0e-6  # vertical emittance in m
+sigma_z = 0.2  # bunch length in m
+Nb = 1e8  # number of particles in the bunch
+
+# Space charge and IBS parameters
+q_val = 1.0  # q-value of longitudinal profile, for space charge
+num_spacecharge_interactions = 1080
+ibs_step = 100  # number of turns between recomputing the IBS coefficients
+
+# Load the line from a file
+line = xt.Line.from_json('SPS_2021_Pb_nominal_deferred_exp.json')
+harmonic_nb = 4653
+
+# Set RF voltage to correct value
+line['actcse.31632' ].voltage = 3.0e6 
+print('RF voltage set to {:.3e} V\n'.format(line['actcse.31632'].voltage))
+
+# Twiss command, inspect tunes and reference particle
+tw = line.twiss() 
+print('Tunes: Qx = {:.6f}, Qy = {:.6f}'.format(tw.qx, tw.qy))
+print('Reference particle: {}'.format(line.particle_ref.show()))
+
+# Add longitudinal limit rectangle - to kill particles that fall out of bucket
+bucket_length = line.get_length()/harmonic_nb
+print('\nBucket length is {:.4f} m'.format(bucket_length))
+line.unfreeze() # if you had already build the tracker
+line.append_element(element=xt.LongitudinalLimitRect(min_zeta=-bucket_length/2, max_zeta=bucket_length/2), name='long_limit')
+line.build_tracker(_context=context)
+
+# Generate a particle distribution
+particles = xp.generate_matched_gaussian_bunch(_context=context,
+    num_particles=n_particles, 
+    total_intensity_particles=Nb,
+    nemitt_x=exn, 
+    nemitt_y=eyn, 
+    sigma_z=sigma_z,
+    particle_ref=line.particle_ref, 
+    line=line)
+
+# Initialize the dataclasses to store particle values
+tbt = Records.init_zeroes(n_turns)  # only emittances and bunch intensity
+tbt.update_at_turn(0, particles, tw)
+tbt.store_initial_particles(particles)
+tbt.store_twiss(tw.to_pandas())
+
+######### Frozen space charge #########
+
+# Store the initial buffer of the line
+_buffer = line._buffer
+line.discard_tracker()
+
+# Install space charge
+lprofile = xf.LongitudinalProfileQGaussian(
+        number_of_particles = Nb,
+        sigma_z = sigma_z,
+        z0=0.,
+        q_parameter=q_val)
+
+# Install frozen space charge as base 
+xf.install_spacecharge_frozen(line = line,
+                    particle_ref = line.particle_ref,
+                    longitudinal_profile = lprofile,
+                    nemitt_x = exn, nemitt_y = eyn,
+                    sigma_z = sigma_z,
+                    num_spacecharge_interactions = num_spacecharge_interactions)
+line.build_tracker(_buffer=_buffer)
+
+######### IBS kinetic kicks #########
+
+#  friction and diffusion terms of the kinetic theory of gases
+ibs_kick = xf.IBSKineticKick(num_slices=50)
+
+### Install the IBS kinetic kick element ###
+#line.configure_intrabeam_scattering(
+#    element=ibs_kick, name="ibskick", index=-1, update_every=ibs_step
+#)
+
+# THESE LINES ABOVE WILL NOT WORK if space charge is already installed
+# Instead, follow manual steps Felix Soubelet's tips
+# Directly copy steps from https://github.com/xsuite/xfields/blob/6882e0d03bb6772f873ce57ef6cf2592e5779359/xfields/ibs/_api.py
+_buffer = line._buffer
+line.discard_tracker()
+line.insert_element(element=ibs_kick, name="ibskick", index=-1)
+line.build_tracker(_buffer=_buffer)
+
+line_sc_off = line.filter_elements(exclude_types_starting_with='SpaceCh')
+twiss_no_sc = line_sc_off.twiss(method="4d")
+
+# Figure out the IBS kick element and its name in the line
+only_ibs_kicks = {name: element for name, element in line.element_dict.items() if isinstance(element, xf.ibs._kicks.IBSKick)}
+assert len(only_ibs_kicks) == 1, "Only one 'IBSKick' element should be present in the line"
+name, element = only_ibs_kicks.popitem()
+
+# Set necessary (private) attributes for the kick to function
+element.update_every = ibs_step
+element._name = name
+element._twiss = twiss_no_sc
+element._scale_strength = 1  # element is now ON, will track
+
+print('\nFixed IBS coefficient recomputation at interval = {} steps\n'.format(ibs_step))
+
+#### Track the particles ####
+time00 = time.time()
+
+for turn in range(1, n_turns):
+    
+    # Print out info at specified intervals
+    if turn % 5 == 0:
+        print('\nTracking turn {}'.format(turn))        
+
+    # ----- Track and update records for tracked particles ----- #
+    line.track(particles, num_turns=1)
+
+    tbt.update_at_turn(turn, particles, tw)
+
+time01 = time.time()
+dt0 = time01-time00
+print('\nTracking time: {:.1f} s = {:.1f} h'.format(dt0, dt0/3600))
+
+# Final turn-by-turn records
+tbt.store_final_particles(particles)
+tbt.to_dict(convert_to_numpy=True)
+
+# then save tbt dict if desired

examples/minimal_working_examples/000_SPS_Pb_with_space_charge_ibs_and_track.py

@@ -0,0 +1,43 @@
+asttokens @ file:///home/conda/feedstock_root/build_artifacts/asttokens_1733250440834/work
+certifi==2025.8.3
+cffi==1.17.1
+charset-normalizer==3.4.3
+Cython==3.1.2
+decorator @ file:///home/conda/feedstock_root/build_artifacts/decorator_1740384970518/work
+exceptiongroup @ file:///home/conda/feedstock_root/build_artifacts/exceptiongroup_1746947292760/work
+executing @ file:///home/conda/feedstock_root/build_artifacts/executing_1745502089858/work
+idna==3.10
+ipython @ file:///home/conda/feedstock_root/build_artifacts/bld/rattler-build_ipython_1751465044/work
+ipython_pygments_lexers @ file:///home/conda/feedstock_root/build_artifacts/ipython_pygments_lexers_1737123620466/work
+jedi @ file:///home/conda/feedstock_root/build_artifacts/jedi_1733300866624/work
+lark==1.2.2
+matplotlib-inline @ file:///home/conda/feedstock_root/build_artifacts/matplotlib-inline_1733416936468/work
+numpy @ file:///home/conda/feedstock_root/build_artifacts/bld/rattler-build_numpy_1753401560/work/dist/numpy-2.3.2-cp311-cp311-linux_x86_64.whl#sha256=469440884ff65cdd5601d2ff9e1c5b9bcc2d8176c37e37ec78a0666703a54157
+pandas==2.3.1
+parso @ file:///home/conda/feedstock_root/build_artifacts/parso_1733271261340/work
+pexpect @ file:///home/conda/feedstock_root/build_artifacts/pexpect_1733301927746/work
+pickleshare @ file:///home/conda/feedstock_root/build_artifacts/pickleshare_1733327343728/work
+prompt_toolkit @ file:///home/conda/feedstock_root/build_artifacts/prompt-toolkit_1744724089886/work
+ptyprocess @ file:///home/conda/feedstock_root/build_artifacts/ptyprocess_1733302279685/work/dist/ptyprocess-0.7.0-py2.py3-none-any.whl#sha256=92c32ff62b5fd8cf325bec5ab90d7be3d2a8ca8c8a3813ff487a8d2002630d1f
+pure_eval @ file:///home/conda/feedstock_root/build_artifacts/pure_eval_1733569405015/work
+pycparser==2.22
+Pygments @ file:///home/conda/feedstock_root/build_artifacts/pygments_1750615794071/work
+python-dateutil==2.9.0.post0
+pytz==2025.2
+requests==2.32.4
+scipy @ file:///home/conda/feedstock_root/build_artifacts/scipy-split_1751147407503/work/dist/scipy-1.16.0-cp311-cp311-linux_x86_64.whl#sha256=af2237dd5dad564373927317264a4269e1a7afc6dc576b23ef2fd3e6f93926c5
+six==1.17.0
+stack_data @ file:///home/conda/feedstock_root/build_artifacts/stack_data_1733569443808/work
+tqdm==4.67.1
+traitlets @ file:///home/conda/feedstock_root/build_artifacts/traitlets_1733367359838/work
+typing_extensions @ file:///home/conda/feedstock_root/build_artifacts/bld/rattler-build_typing_extensions_1751643513/work
+tzdata==2025.2
+urllib3==2.5.0
+wcwidth @ file:///home/conda/feedstock_root/build_artifacts/wcwidth_1733231326287/work
+xcoll==0.6.2
+xdeps==0.10.5
+xfields==0.25.1
+xobjects==0.5.2
+xpart==0.23.1
+xsuite==0.36.2
+xtrack==0.88.3

examples/minimal_working_examples/SPS_2021_Pb_nominal_deferred_exp.json

@@ -0,0 +1,136 @@
+"""
+Container of simplistic record class to store values during tracking
+"""
+import numpy as np
+import xpart as xp
+import xtrack as xt
+import xobjects as xo
+from dataclasses import dataclass
+from typing import Self
+import pandas as pd
+import json
+
+####### Helper functions for bunch length, momentum spread and geometric emittances #######
+def _bunch_length(parts: xp.Particles) -> float:
+    return np.std(parts.zeta[parts.state > 0])
+
+
+def _sigma_delta(parts: xp.Particles) -> float:
+    return np.std(parts.delta[parts.state > 0])
+
+
+def _geom_epsx(parts: xp.Particles, twiss: xt.TwissTable) -> float:
+    """
+    We index dx and betx at 0 which corresponds to the beginning / end of
+    the line, since this is where / when we will be applying the kicks.
+    """
+    delta = parts.delta[parts.state > 0]
+    sigma_x = np.std(parts.x[parts.state > 0])
+    sig_delta = _sigma_delta(parts)
+    return (sigma_x**2 - (twiss["dx"][0] * sig_delta) ** 2) / twiss["betx"][0]
+
+
+def _geom_epsy(parts: xp.Particles, twiss: xt.TwissTable) -> float:
+    """
+    We index dy and bety at 0 which corresponds to the beginning / end of
+    the line, since this is where / when we will be applying the kicks.
+    """
+    delta = parts.delta[parts.state > 0]
+    sigma_y =  np.std(parts.y[parts.state > 0])
+    sig_delta = _sigma_delta(parts)
+    return (sigma_y**2 - (twiss["dy"][0] * sig_delta) ** 2) / twiss["bety"][0]
+
+#############################################################################
+
+@dataclass
+class Records:
+    """
+    Data class to store numpy.ndarray of results during flat-bottom tracking 
+    """
+    exn: np.ndarray
+    eyn: np.ndarray
+    sigma_delta: np.ndarray
+    bunch_length: np.ndarray
+    Nb: np.ndarray
+    turns: np.ndarray
+    twiss: dict
+    particles_i: dict
+    particles_f: dict
+
+    def update_at_turn(self, turn: int, parts: xp.Particles, twiss: xt.TwissTable):
+        """Automatically update the records at given turn from the xpart.Particles."""
+
+        # Store the particle ensemble quantities
+        self.exn[turn] = _geom_epsx(parts, twiss) * parts.beta0[0] * parts.gamma0[0]
+        self.eyn[turn] = _geom_epsy(parts, twiss) * parts.beta0[0] * parts.gamma0[0]
+        self.Nb[turn] = parts.weight[parts.state > 0][0]*len(parts.x[parts.state > 0])
+        self.sigma_delta[turn] = _sigma_delta(parts)
+        self.bunch_length[turn] = _bunch_length(parts)
+
+    @classmethod
+    def init_zeroes(cls, n_turns: int) -> Self:  # noqa: F821
+        """Initialize the dataclass with arrays of zeroes."""
+        return cls(
+            exn=np.zeros(n_turns, dtype=float),
+            eyn=np.zeros(n_turns, dtype=float),
+            Nb=np.zeros(n_turns, dtype=float),
+            sigma_delta=np.zeros(n_turns, dtype=float),
+            bunch_length=np.zeros(n_turns, dtype=float),
+            turns=np.arange(n_turns, dtype=int),
+            twiss={},
+            particles_i={},
+            particles_f={}
+        )
+    
+    def store_twiss(self, df_twiss: pd.DataFrame):
+        """Store twiss table, before collective elements"""
+        self.twiss = df_twiss.to_dict()
+        self.includes_particle_and_twiss_data = True
+
+    def store_initial_particles(self, parts: xp.Particles):
+        """Store initial particle object"""
+        self.particles_i = parts.to_dict()
+        self.includes_particle_and_twiss_data = True
+
+    def store_final_particles(self, parts: xp.Particles):
+        """Store final particle object"""
+        self.particles_f = parts.to_dict()
+        self.includes_particle_and_twiss_data = True
+
+    def to_dict(self, convert_to_numpy=True):
+        """
+        Convert data arrays to dictionary, possible also beam profile monitor data
+        Convert lists to numpy format if desired, but typically not if data is saved to json
+        """
+        data = {
+            'exn': self.exn.tolist(),
+            'eyn': self.eyn.tolist(),
+            'sigma_delta': self.sigma_delta.tolist(),
+            'bunch_length': self.bunch_length.tolist(),
+            'Nb' : self.Nb.tolist(),
+            'Turns': self.turns.tolist()
+        }
+        data['twiss'] = self.twiss
+        data['particles_i'] = self.particles_i
+        data['particles_f'] = self.particles_f 
+
+        # Convert lists to numpy arrays if desired
+        if convert_to_numpy:
+            for key, value in data.items():
+                if isinstance(data[key], list):
+                    data[key] = np.array(data[key])
+                    
+        return data
+
+
+    def to_json(self, file_path=None):
+        """
+        Save the data to a JSON file.
+        """
+        if file_path is None:
+            file_path = './'
+
+        data = self.to_dict()
+
+        with open('{}tbt.json'.format(file_path), 'w') as f:
+            json.dump(data, f, cls=xo.JEncoder)