From 4f9bc46ffee944702c1e7af895f7d8e66f89cf46 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Thu, 21 Aug 2025 16:05:57 +0200 Subject: [PATCH 01/37] Add `compt_scale` as a user parameter. `compt_scale` is a multiplicative factor applied to the total Compton cross-section. Increasing this factor scales up photon generation in the radiative Bhabha model. This can be useful for testing or for producing a larger number of radiative Bhabha scattering events. --- xfields/beam_elements/beambeam3d.py | 10 ++++++++-- xfields/beam_elements/beambeam_src/beambeam3d.h | 4 +++- xfields/headers/bhabha_spectrum.h | 2 +- 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/xfields/beam_elements/beambeam3d.py b/xfields/beam_elements/beambeam3d.py index 46dbbaa6..e32103e9 100644 --- a/xfields/beam_elements/beambeam3d.py +++ b/xfields/beam_elements/beambeam3d.py @@ -181,6 +181,7 @@ class BeamBeamBiGaussian3D(xt.BeamElement): #bhabha 'flag_bhabha': xo.Int64, 'compt_x_min': xo.Float64, + 'compt_scale': xo.Float64, 'flag_beamsize_effect': xo.Int64, #lumi @@ -251,6 +252,7 @@ def __init__(self, flag_bhabha=0, compt_x_min=1e-4, + compt_scale=1., flag_beamsize_effect=1, flag_luminosity=0, @@ -330,6 +332,7 @@ def __init__(self, slices_other_beam_sqrtSigma_{135}{135}_beamstrahlung (float array): Array storing the per-slice standard deviations of x (=1), y (=3) and zeta (=5) of the opposing bunch, in the unboosted accelerator frame. Length of the array is the number of longitudinal slices. Used for beamstrahlung only. flag_bhabha (int): Flag to simulate small angle radiative Bhabha scattering. 1: ON (quantum), 0: OFF compt_x_min (float): Low energy cut on virtual photon spectrum, used for Bhabha, in units of [gamma^-2] where gamma is the rel. Lorentz factor. + compt_scale (float): Scaling factor to scale up photon generation, used for Bhabha. flag_beamsize_effect (int): Flag to simulate beamsize effect, used for Bhabha. 1: ON, 0: OFF. Results in ~factor 2 reduction in cross section. flag_luminosity (int): Flag to record soft-Gaussian luminosity per bunch crossing in a buffer. Luminosity will be in units of [m^-2]. slices_other_beam_{x,px,y,py,zeta,pzeta}_center_star (float array): Array storing the per-slice centroid variables of the opposing bunch, in the unboosted accelerator frame. Length of the array is the number of longitudinal slices. @@ -465,7 +468,7 @@ def __init__(self, ) # initialize bhabha - self._init_bhabha(flag_bhabha, compt_x_min, flag_beamsize_effect) + self._init_bhabha(flag_bhabha, compt_x_min, compt_scale, flag_beamsize_effect) self._init_luminosity(flag_luminosity) @@ -635,9 +638,10 @@ def flag_beamstrahlung(self, flag_beamstrahlung): 'needs to be correctly set') self._flag_beamstrahlung = flag_beamstrahlung - def _init_bhabha(self, flag_bhabha, compt_x_min, flag_beamsize_effect): + def _init_bhabha(self, flag_bhabha, compt_x_min, compt_scale, flag_beamsize_effect): self.flag_beamsize_effect = flag_beamsize_effect self.compt_x_min = compt_x_min + self.compt_scale = compt_scale self.flag_bhabha = flag_bhabha # Trigger property setter @property @@ -649,6 +653,8 @@ def flag_bhabha(self, flag_bhabha): if flag_bhabha == 1: if self.compt_x_min <= 0: raise ValueError('compt_x_min must be larger than 0') + if self.compt_scale <= 0: + raise ValueError('compt_scale must be larger than 0') self._flag_bhabha = flag_bhabha def _init_luminosity(self, flag_luminosity): diff --git a/xfields/beam_elements/beambeam_src/beambeam3d.h b/xfields/beam_elements/beambeam_src/beambeam3d.h index e028c020..482d9c36 100644 --- a/xfields/beam_elements/beambeam_src/beambeam3d.h +++ b/xfields/beam_elements/beambeam_src/beambeam3d.h @@ -103,9 +103,11 @@ void do_bhabha(BeamBeamBiGaussian3DData el, LocalParticle *part, py_photon = py_slice_star; pzeta_photon = pzeta_slice_star; + const double compt_scale = BeamBeamBiGaussian3DData_get_compt_scale(el); // [1] can be used to scale up photon generation for tests + // for each virtual photon get compton scatterings; updates pzeta and energy vars inside compt_do(part, bhabha_record, bhabha_table_index, bhabha_table, - e_photon, compt_x_min, q2, + e_photon, compt_x_min, compt_scale, q2, x_photon, y_photon, S, px_photon, py_photon, pzeta_photon, *wgt, px_star, py_star, pzeta_star, q0); diff --git a/xfields/headers/bhabha_spectrum.h b/xfields/headers/bhabha_spectrum.h index 0e2f9c11..a732842d 100644 --- a/xfields/headers/bhabha_spectrum.h +++ b/xfields/headers/bhabha_spectrum.h @@ -250,6 +250,7 @@ double compt_select(LocalParticle *part, void compt_do(LocalParticle *part, BeamBeamBiGaussian3DRecordData bhabha_record, RecordIndex bhabha_table_index, BhabhaTableData bhabha_table, double e_photon, // [GeV] single equivalent virtual photon energy before Compton scattering const double compt_x_min, // [1] scaling factor in the minimum energy cutoff + double compt_scale, // [1] can be used to scale up photon generation for tests double q2, // [GeV^2] single equivalent virtual photon virtuality double x_photon, double y_photon, double z_photon, // [m] (boosted) coords of the virtual photon double vx_photon, // [1] transverse x momentum component of virtual photon (vx = dx/ds/p0) @@ -284,7 +285,6 @@ void compt_do(LocalParticle *part, BeamBeamBiGaussian3DRecordData bhabha_record, double e_photon_prime, px_photon_prime, py_photon_prime; // [GeV, GeV/c] scattered (real) Compton photon momenta double e_loss_primary; // [GeV] energy lost from one photon emission double eps = 0.0; // 1e-5 in guinea, used due to limited resolution of RNG in the past - const double compt_scale = 1; // [1] can be used to scale up photon generation for tests const double compt_emax = 200; // [GeV] upper cutoff on e_e_prime from guineapig const double pair_ecut = 0.005; // [GeV] lower cutoff on e_e_prime from guineapig double r1, r2; // [1] uniform random numbers From 49b473e6e9651108a651f9da51003be14d3dc28e Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Thu, 4 Sep 2025 16:55:15 +0200 Subject: [PATCH 02/37] Touschek! --- xfields/__init__.py | 3 + xfields/beam_elements/touschek.py | 201 ++++++++ xfields/beam_elements/touschek_src/touschek.h | 488 ++++++++++++++++++ xfields/headers/constants.h | 4 + xfields/touschek/manager.py | 305 +++++++++++ 5 files changed, 1001 insertions(+) create mode 100644 xfields/beam_elements/touschek.py create mode 100644 xfields/beam_elements/touschek_src/touschek.h create mode 100644 xfields/touschek/manager.py diff --git a/xfields/__init__.py b/xfields/__init__.py index e68c3215..e059e685 100644 --- a/xfields/__init__.py +++ b/xfields/__init__.py @@ -16,6 +16,8 @@ from .solvers.fftsolvers import FFTSolver3D +from .touschek.manager import TouschekManager + from .beam_elements.spacecharge import SpaceCharge3D, SpaceChargeBiGaussian from .beam_elements.beambeam2d import BeamBeamBiGaussian2D from .beam_elements.beambeam2d import ConfigForUpdateBeamBeamBiGaussian2D @@ -25,6 +27,7 @@ from .beam_elements.temp_slicer import TempSlicer from .beam_elements.electroncloud import ElectronCloud from .beam_elements.electronlens_interpolated import ElectronLensInterpolated +from .beam_elements.touschek import TouschekScattering from .general import _pkg_root from .config_tools import replace_spacecharge_with_quasi_frozen diff --git a/xfields/beam_elements/touschek.py b/xfields/beam_elements/touschek.py new file mode 100644 index 00000000..2f5a4b9d --- /dev/null +++ b/xfields/beam_elements/touschek.py @@ -0,0 +1,201 @@ +import xobjects as xo +import xtrack as xt + +import numpy as np +from scipy.integrate import quad +from scipy.special import i0 +from scipy.constants import physical_constants + +from ..general import _pkg_root + +# =========================================== +# Constants +# =========================================== +ELECTRON_MASS_EV = xt.ELECTRON_MASS_EV +C_LIGHT_VACUUM = physical_constants['speed of light in vacuum'][0] +CLASSICAL_ELECTRON_RADIUS = physical_constants['classical electron radius'][0] + +class TouschekScattering(xt.BeamElement): + + _xofields = { + # All these are handled by the manager for now + + # Everything that need to be accessed with TouschekElementData_get in C + '_p0c': xo.Float64, + '_bunch_population': xo.Float64, + '_gemitt_x': xo.Float64, + '_gemitt_y': xo.Float64, + '_alfx': xo.Float64, + '_betx': xo.Float64, + '_alfy': xo.Float64, + '_bety': xo.Float64, + '_dx': xo.Float64, + '_dpx': xo.Float64, + '_dy': xo.Float64, + '_dpy': xo.Float64, + '_deltaN': xo.Float64, + '_deltaP': xo.Float64, + '_sigma_z': xo.Float64, + '_sigma_delta': xo.Float64, + '_n_simulated': xo.Int64, + '_nx': xo.Float64, + '_ny': xo.Float64, + '_nz': xo.Float64, + '_ignored_portion': xo.Float64, + '_integrated_piwinski_rate': xo.Float64, + } + + # allow_track = False + _depends_on = [xt.RandomUniformAccurate] + + _extra_c_sources = [ + _pkg_root.joinpath('headers/constants.h'), + _pkg_root.joinpath('beam_elements/touschek_src/touschek.h'), + ] + + _per_particle_kernels = { + '_scatter': xo.Kernel( + c_name='TouschekScatter', + args=[ + # xo.Arg(xo.Int64, name='n_simulated'), + # xo.Arg(xo.Float64, name='nx'), + # xo.Arg(xo.Float64, name='ny'), + # xo.Arg(xo.Float64, name='nz'), + # xo.Arg(xo.Float64, name='ignored_portion'), + + # These are variables that we will get out from the C kernel + # We will need to pre-allocate memory for them in Python + # Then the C kernel will fill them in + xo.Arg(xo.Float64, name='x_out', pointer=True), + xo.Arg(xo.Float64, name='px_out', pointer=True), + xo.Arg(xo.Float64, name='y_out', pointer=True), + xo.Arg(xo.Float64, name='py_out', pointer=True), + xo.Arg(xo.Float64, name='zeta_out', pointer=True), + xo.Arg(xo.Float64, name='delta_out', pointer=True), + xo.Arg(xo.Float64, name='weight_out', pointer=True), + xo.Arg(xo.Float64, name='totalMCRate_out', pointer=True), + xo.Arg(xo.Int64, name='n_selected_out', pointer=True), + ], + ), + } + + def __init__(self, + s=0.0, + particle_ref=xt.Particles(), + bunch_population=0.0, + alfx=0.0, betx=0.0, alfy=0.0, bety=0.0, + dx=0.0, dpx=0.0, dy=0.0, dpy=0.0, + deltaN=0.0, deltaP=0.0, + gemitt_x=0.0, gemitt_y=0.0, + sigma_z=0.0, sigma_delta=0.0, + n_simulated=0, nx=0.0, ny=0.0, nz=0.0, + piwinski_rate=0.0, + ignored_portion=0.0, + integrated_piwinski_rate=0.0, + **kwargs): + + # This gives AttributeError: 'TouschekScattering' object has no attribute '_xobject' + # if not isinstance(self._context, xo.ContextCpu) or self._context.openmp_enabled: + # raise ValueError('TouschekScattering only enabled on CPU.') + + if '_xobject' in kwargs.keys(): + self.xoinitialize(**kwargs) + return + + super().__init__(**kwargs) + + self._s = s + self._particle_ref = particle_ref + self._bunch_population = bunch_population + self._alfx = alfx + self._betx = betx + self._alfy = alfy + self._bety = bety + self._dx = dx + self._dpx = dpx + self._dy = dy + self._dpy = dpy + self._deltaN = deltaN + self._deltaP = deltaP + self._gemitt_x = gemitt_x + self._gemitt_y = gemitt_y + self._sigma_z = sigma_z + self._sigma_delta = sigma_delta + self._n_simulated = n_simulated + self._nx = nx + self._ny = ny + self._nz = nz + self._ignored_portion = ignored_portion + self._integrated_piwinski_rate = integrated_piwinski_rate + + self.piwinski_rate = piwinski_rate + + + def _configure(self, **kwargs): + config_allowed = { + "_particle_ref", "_bunch_population", + "_gemitt_x", "_gemitt_y", + "_alfx", "_betx", "_alfy", "_bety", + "_dx", "_dpx", "_dy", "_dpy", + "_deltaN", "_deltaP", + "_sigma_z", "_sigma_delta", + "_n_simulated", "_nx", "_ny", "_nz", + "_ignored_portion", "piwinski_rate", + "_integrated_piwinski_rate" + } + + unknown = set(kwargs) - config_allowed + if unknown: + bad = ", ".join(sorted(unknown)) + raise KeyError(f"Unsupported configure() keys: {bad}") + + for kk, vv in kwargs.items(): + setattr(self, kk, vv) + if kk == "_particle_ref": + self._p0c = self._particle_ref.p0c[0] + + + def scatter(self): + context = self._context + particles = xt.Particles(_context=context) + + if not particles._has_valid_rng_state(): + particles._init_random_number_generator() + + x_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) + px_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) + y_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) + py_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) + zeta_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) + delta_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) + weight_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) + totalMCRate_out = context.zeros(shape=(1,), dtype=np.float64) + n_selected_out = context.zeros(shape=(1,), dtype=np.int64) + + self._scatter(particles=particles, + x_out=x_out, px_out=px_out, + y_out=y_out, py_out=py_out, + zeta_out=zeta_out, delta_out=delta_out, + weight_out=weight_out, + totalMCRate_out=totalMCRate_out, + n_selected_out=n_selected_out) + + n = n_selected_out[0] + part = xt.Particles(_capacity=2*n, + p0c=self._p0c, + mass0=self._particle_ref.mass0, + q0=self._particle_ref.q0, + pdg_id=self._particle_ref.pdg_id, + x=x_out[:n], px=px_out[:n], + y=y_out[:n], py=py_out[:n], + zeta=zeta_out[:n], delta=delta_out[:n], + weight=weight_out[:n], + s=getattr(self, '_s', 0.0)) + + self.particles = part + self.total_mc_rate = totalMCRate_out + self.ignored_rate = self._ignored_portion * self.total_mc_rate + + + def track(self, particles): + super().track(particles) \ No newline at end of file diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h new file mode 100644 index 00000000..b9090573 --- /dev/null +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -0,0 +1,488 @@ +#ifndef XTRACK_TOUSCHEK_H +#define XTRACK_TOUSCHEK_H + +#include +#include +#include +#include + +static inline double sqr(double x){ return x*x; } + +/*gpufun*/ +void TouschekScattering_track_local_particle(TouschekScatteringData el, LocalParticle* part0) { + (void)el; (void)part0; + return; +} + +/*gpufun*/ +void shuffle_double_array(double* a, int n, LocalParticle* part) { + for (int i = n - 1; i > 0; --i) { + // Option A: use floating RNG + int j = (int)floor(RandomUniformAccurate_generate(part) * (i + 1)); + + double tmp = a[i]; + a[i] = a[j]; + a[j] = tmp; + } +} + + +void selectPartGauss(double *p1, double *p2, + double *dens1, double *dens2, + const double *ran1, + const double *range, + const double *alfa, + const double *beta, + const double *disp, + const double *gemitt) { + int i; + double U[3], V1[3], V2[3], densa[3], densb[3]; + + /* Select random particle coordinates in normalized phase space */ + for (i = 0; i < 3; i++) { + U[i] = (ran1[i] - 0.5) * range[i] * sqrt(gemitt[i]); + V1[i] = (ran1[i + 3] - 0.5) * range[i] * sqrt(gemitt[i]); + V2[i] = (ran1[i + 6] - 0.5) * range[i] * sqrt(gemitt[i]); + densa[i] = exp(-0.5 * (U[i] * U[i] + V1[i] * V1[i]) / gemitt[i]); + } + densb[2] = exp(-0.5 * (U[2] * U[2] + V2[2] * V2[2]) / gemitt[2]); + /* Transform particle coordinates from normalized to real phase space */ + for (i = 0; i < 3; i++) { + p1[i] = p2[i] = sqrt(beta[i]) * U[i]; + p1[i + 3] = (V1[i] - alfa[i] * U[i]) / sqrt(beta[i]); + p2[i + 3] = (V2[i] - alfa[i] * U[i]) / sqrt(beta[i]); + } + /* Dispersion correction */ + p1[0] = p1[0] + p1[5] * disp[0]; + p1[1] = p1[1] + p1[5] * disp[1]; + p1[3] = p1[3] + p1[5] * disp[2]; + p1[4] = p1[4] + p1[5] * disp[3]; + + p2[0] = p1[0] - p2[5] * disp[0]; + p2[1] = p1[1] - p2[5] * disp[1]; + U[0] = p2[0] / sqrt(beta[0]); + U[1] = p2[1] / sqrt(beta[1]); + p2[3] = (V2[0] - alfa[0] * U[0]) / sqrt(beta[0]); + p2[4] = (V2[1] - alfa[1] * U[1]) / sqrt(beta[1]); + densb[0] = exp(-0.5 * (U[0] * U[0] + V2[0] * V2[0]) / gemitt[0]); + densb[1] = exp(-0.5 * (U[1] * U[1] + V2[1] * V2[1]) / gemitt[1]); + + p2[0] = p1[0]; + p2[1] = p1[1]; + p2[3] = p2[3] + p2[5] * disp[2]; + p2[4] = p2[4] + p2[5] * disp[3]; + + *dens1 = densa[0] * densa[1] * densa[2]; + *dens2 = densb[0] * densb[1] * densb[2]; + + return; +} + +void bunch2cm(double *p1, double *p2, double *q, double *beta, double *gamma) { + double pp1, pp2, e1, e2, ee; + int i; + double bb, betap1, factor; + + pp1 = 0.0; + pp2 = 0.0; + for (i = 3; i < 6; i++) { + pp1 = pp1 + sqr(p1[i]); + pp2 = pp2 + sqr(p2[i]); + } + e1 = sqrt(MELECTRON_EV * MELECTRON_EV + pp1); + e2 = sqrt(MELECTRON_EV * MELECTRON_EV + pp2); + ee = e1 + e2; + + betap1 = 0.0; + bb = 0.0; + for (i = 0; i < 3; i++) { + beta[i] = (p1[i + 3] + p2[i + 3]) / ee; + betap1 = betap1 + beta[i] * p1[i + 3]; + bb = bb + beta[i] * beta[i]; + } + + *gamma = 1. / sqrt(1. - bb); + factor = ((*gamma) - 1.) * betap1 / bb; + + for (i = 0; i < 3; i++) { + q[i] = p1[i + 3] + factor * beta[i] - (*gamma) * e1 * beta[i]; + } + + return; +} + +/* Rotate scattered p in c.o.m system */ +void eulertrans(double *v0, double theta, double phi, double *v1, double *v) { + double th, ph, s1, s2, c1, c2; + double x0, y0, z0; + + *v = sqrt(v0[0] * v0[0] + v0[1] * v0[1] + v0[2] * v0[2]); + th = acos(v0[2] / (*v)); + ph = atan2(v0[1], v0[0]); + + s1 = sin(th); + s2 = sin(ph); + c1 = cos(th); + c2 = cos(ph); + + x0 = cos(theta); + y0 = sin(theta) * cos(phi); + z0 = sin(theta) * sin(phi); + + v1[0] = (*v) * (s1 * c2 * x0 - s2 * y0 - c1 * c2 * z0); + v1[1] = (*v) * (s1 * s2 * x0 + c2 * y0 - c1 * s2 * z0); + v1[2] = (*v) * (c1 * x0 + s1 * z0); + + return; +} + +void cm2bunch(double *p1, double *p2, double *q, double *beta, double *gamma) { + int i; + double pq, e, betaq, bb, factor; + + pq = 0.0; + for (i = 0; i < 3; i++) { + pq = pq + q[i] * q[i]; + } + + e = sqrt(MELECTRON_EV * MELECTRON_EV + pq); + + betaq = 0.0; + bb = 0.0; + for (i = 0; i < 3; i++) { + betaq = betaq + beta[i] * q[i]; + bb = bb + beta[i] * beta[i]; + } + + factor = ((*gamma) - 1) * betaq / bb; + for (i = 0; i < 3; i++) { + p1[i + 3] = q[i] + (*gamma) * beta[i] * e + factor * beta[i]; + p2[i + 3] = -q[i] + (*gamma) * beta[i] * e - factor * beta[i]; + } + + return; +} + +double moeller(double beta0, double theta) { + double cross; + double beta2, st2; + + beta2 = beta0 * beta0; + st2 = sqr(sin(theta)); + + cross = (1. - beta2) * (sqr(1. + 1. / beta2) * (4. / st2 / st2 - 3. / st2) + 1. + 4. / st2); + + return cross; +} + +void pickPart(double *weight, long *index, long start, long end, + long *iTotal, double *wTotal, double weight_limit, double weight_ave) { + long i, i1, i2, N; + double w1, w2; + long *index1, *index2; + double *weight1, *weight2; + + i1 = i2 = 0; + w1 = w2 = 0.; + N = end - start; + if (N < 3) + return; /* scattered particles normally appear in pair */ + index2 = (long *)malloc(sizeof(long) * N); + weight2 = (double *)malloc(sizeof(double) * N); + index1 = (long *)malloc(sizeof(long) * N); + weight1 = (double *)malloc(sizeof(double) * N); + + for (i = start; i < end; i++) { + if (weight[i] > weight_ave) { + weight2[i2] = weight[i]; + index2[i2++] = index[i]; + w2 += weight[i]; + } else { + weight1[i1] = weight[i]; + index1[i1++] = index[i]; + w1 += weight[i]; + } + } + if ((w2 + (*wTotal)) > weight_limit) { + weight_ave = w2 / (double)i2; + for (i = 0; i < i2; i++) { + index[start + i] = index2[i]; + weight[start + i] = weight2[i]; + } + free(weight1); + free(index1); + free(weight2); + free(index2); + pickPart(weight, index, start, start + i2, + iTotal, wTotal, weight_limit, weight_ave); + return; + } + + *iTotal += i2; + *wTotal += w2; + weight_ave = w1 / (double)i1; + for (i = 0; i < i2; i++) { + index[start + i] = index2[i]; + weight[start + i] = weight2[i]; + } + for (i = 0; i < i1; i++) { + index[start + i2 + i] = index1[i]; + weight[start + i2 + i] = weight1[i]; + } + free(weight1); + free(index1); + free(weight2); + free(index2); + pickPart(weight, index, i2 + start, end, + iTotal, wTotal, weight_limit, weight_ave); + return; +} + +/*gpufun*/ +void TouschekScatter(TouschekScatteringData el, + LocalParticle* part0, + double* x_out, + double* px_out, + double* y_out, + double* py_out, + double* zeta_out, + double* delta_out, + double* weight_out, + double* totalMCRate_out, + int64_t* n_selected_out){ + + const double p0c = TouschekScatteringData_get__p0c(el); + const double bunch_population = TouschekScatteringData_get__bunch_population(el); + const double gemitt_x = TouschekScatteringData_get__gemitt_x(el); + const double gemitt_y = TouschekScatteringData_get__gemitt_y(el); + const double alfx = TouschekScatteringData_get__alfx(el); + const double betx = TouschekScatteringData_get__betx(el); + const double alfy = TouschekScatteringData_get__alfy(el); + const double bety = TouschekScatteringData_get__bety(el); + const double dx = TouschekScatteringData_get__dx(el); + const double dpx = TouschekScatteringData_get__dpx(el); + const double dy = TouschekScatteringData_get__dy(el); + const double dpy = TouschekScatteringData_get__dpy(el); + const double deltaN = TouschekScatteringData_get__deltaN(el); + const double deltaP = TouschekScatteringData_get__deltaP(el); + const double sigma_z = TouschekScatteringData_get__sigma_z(el); + const double sigma_delta = TouschekScatteringData_get__sigma_delta(el); + const double n_simulated = TouschekScatteringData_get__n_simulated(el); + const double nx = TouschekScatteringData_get__nx(el); + const double ny = TouschekScatteringData_get__ny(el); + const double nz = TouschekScatteringData_get__nz(el); + const double ignoredPortion = TouschekScatteringData_get__ignored_portion(el); + const double integrated_piwinski_rate = TouschekScatteringData_get__integrated_piwinski_rate(el); + + long i, j, total_event, simuCount, iTotal; + double ran1[11]; + + long *index = NULL; + double *weight = (double*)malloc(sizeof(double) * n_simulated); + + const double twissAlpha[3] = { alfx, alfy, 0.0 }; + + const double bets = sigma_z/sigma_delta; + const double twissBeta[3] = { betx, bety, bets }; + + const double twissDisp[4] = { dx, dy, dpx, dpy }; + + const double gemitt_z = sigma_z*sigma_delta; + const double gemitt[3] = { gemitt_x, gemitt_y, gemitt_z }; + const double range[3] = { 2.0*nx, 2.0*ny, 2.0*nz }; + + double totalWeight = 0.0; + double totalMCRate = 0.0; + double ignoredRate = 0.0; + + double pTemp[6], p1[6], p2[6], dens1, dens2; + double theta, phi, qa[3], qb[3], beta[3], qabs, gamma; + double beta0, cross, temp; + + double weight_limit, weight_ave, wTotal; + + const double sigxyz = sqrt(twissBeta[0]*gemitt[0]) * sqrt(twissBeta[1]*gemitt[1]) * sigma_z; + temp = sqr(bunch_population) * sqr(PI) * sqr(RE) * C_LIGHT / 4.; + double factor = temp * pow(range[0], 3.0) * pow(range[1], 3.0) * pow(range[2], 3.0) / pow(2 * PI, 6.0) / sigxyz; + + double *xtemp = (double*)malloc(sizeof(double) * n_simulated); + double *pxtemp = (double*)malloc(sizeof(double) * n_simulated); + double *ytemp = (double*)malloc(sizeof(double) * n_simulated); + double *pytemp = (double*)malloc(sizeof(double) * n_simulated); + double *zetatemp = (double*)malloc(sizeof(double) * n_simulated); + double *deltatemp = (double*)malloc(sizeof(double) * n_simulated); + + i = 0; + j = 0; + total_event = 0; + simuCount = 0; + + while (1) { + if (i >= n_simulated) + break; + + /* Select 11 random numbers, then mix them. */ + for (j = 0; j < 11; j++) { + ran1[j] = RandomUniformAccurate_generate(part0); // replaces random_1_elegant + } + // shuffle_double_array(ran1, 11, part0); // may replace ELEGANT's randomizeOrder ? + + total_event++; + + selectPartGauss(p1, p2, &dens1, &dens2, ran1, range, twissAlpha, twissBeta, twissDisp, gemitt); + + if (!dens1 || !dens2) { + continue; + } + /* Change from slope to momentum */ + for (j = 3; j < 5; j++) { + p1[j] *= p0c; + p2[j] *= p0c; + } + p1[5] = (p1[5] + 1) * p0c; + p2[5] = (p2[5] + 1) * p0c; + + bunch2cm(p1, p2, qa, beta, &gamma); + + theta = (ran1[9] * 0.9999 + 0.00005) * PI; + phi = ran1[10] * PI; + + temp = dens1 * dens2 * sin(theta); + eulertrans(qa, theta, phi, qb, &qabs); + cm2bunch(p1, p2, qb, beta, &gamma); + p1[5] = (p1[5] - p0c) / p0c; + p2[5] = (p2[5] - p0c) / p0c; + + if (p1[5] > p2[5]) { + for (j = 0; j < 6; j++) { + pTemp[j] = p2[j]; + p2[j] = p1[j]; + p1[j] = pTemp[j]; + } + } + + if (p1[5] < deltaN || p2[5] > deltaP) { + beta0 = qabs / sqrt(qabs * qabs + MELECTRON_EV * MELECTRON_EV); + cross = moeller(beta0, theta); + temp *= cross * beta0 / gamma / gamma; + + if (p1[5] < deltaN) { + totalWeight += temp; + p1[3] /= p0c; + p1[4] /= p0c; + simuCount++; + + xtemp[i] = p1[0]; + pxtemp[i] = p1[3]; + ytemp[i] = p1[1]; + pytemp[i] = p1[4]; + zetatemp[i] = p1[2]; + deltatemp[i] = p1[5]; + weight[i] = temp; + i++; + } + + if (i >= n_simulated) + break; + + if (p2[5] > deltaP) { + totalWeight += temp; + p2[3] /= p0c; + p2[4] /= p0c; + simuCount++; + + xtemp[i] = p2[0]; + pxtemp[i] = p2[3]; + ytemp[i] = p2[1]; + pytemp[i] = p2[4]; + zetatemp[i] = p2[2]; + deltatemp[i] = p2[5]; + weight[i] = temp; + i++; + } + } + // ELEGANT warnings + // if (total_event * 11 > (long)2e9) { + // printWarning("TouschekScatter: the total number of random numbers used exceeded 2e9.", + // "Use smaller n_simulated or smaller delta."); + // break; + // } + // } + // if (verbosity) + // report_stats(stdout, "After particle generation: "); + // if (simuCount == 0) + // bombElegant("It appears that the Touschek lifetime is extremely wrong and there is no need to perform Touschek simulation.\n If you think this is a wrong statement, send input to developers for evaluation.\n", NULL); + + // if (total_event / tsptr->simuCount > 20) { + // if (distribution_cutoff[0] < 5 || distribution_cutoff[1] < 5) + // printWarning("touschek_scatter: scattering rate is low.", + // "Please use >=5 sigma beam for better simulation."); + // else + // printWarning("touschek_scatter: Sscattering rate is very low.", + // "Please ignore the rate from Monte Carlo simulation. Use Piwinski's rate only."); + } + factor = factor / (double)(total_event); + totalMCRate = totalWeight * factor; + ignoredRate = totalMCRate * ignoredPortion; + + /* Pick tracking particles from the simulated scattered particles */ + index = (long *)malloc(sizeof(long) * simuCount); + for (i = 0; i < simuCount; i++) { + index[i] = i; + } + + if (ignoredPortion <= 1e-9) { + iTotal = simuCount; + wTotal = totalWeight; + for (long k = 0; k < iTotal; ++k) { + x_out[k] = xtemp[k]; + px_out[k] = pxtemp[k]; + y_out[k] = ytemp[k]; + py_out[k] = pytemp[k]; + zeta_out[k] = zetatemp[k]; + delta_out[k] = deltatemp[k]; + weight_out[k] = weight[k]; + } + } else { + iTotal = 0; + wTotal = 0.; + weight_limit = totalWeight * (1 - ignoredPortion); + weight_ave = totalWeight / simuCount; + + pickPart(weight, index, 0, simuCount, + &iTotal, &wTotal, weight_limit, weight_ave); + + for (long k = 0; k < iTotal; ++k) { + long src = index[k]; + x_out[k] = xtemp[src]; + px_out[k] = pxtemp[src]; + y_out[k] = ytemp[src]; + py_out[k] = pytemp[src]; + zeta_out[k] = zetatemp[src]; + delta_out[k] = deltatemp[src]; + weight_out[k] = weight[src]; + } + } + + printf("%ld of %ld particles selected for tracking\n", iTotal, simuCount); + fflush(stdout); + + // Update weight_out to match ELEGANT + for (long k = 0; k < iTotal; ++k) { + weight_out[k] *= (factor / totalMCRate) * integrated_piwinski_rate; + } + + *n_selected_out = iTotal; + *totalMCRate_out = totalMCRate; + + free(index); + free(weight); + free(xtemp); + free(pxtemp); + free(ytemp); + free(pytemp); + free(zetatemp); + free(deltatemp); +} + +#endif // XTRACK_TOUSCHEK_H \ No newline at end of file diff --git a/xfields/headers/constants.h b/xfields/headers/constants.h index 22f0af9f..43a84af2 100644 --- a/xfields/headers/constants.h +++ b/xfields/headers/constants.h @@ -51,6 +51,10 @@ #define MELECTRON_GEV (0.00051099895000) #endif +#if !defined( MELECTRON_EV ) + #define MELECTRON_EV (510998.95) +#endif + #if !defined( MELECTRON_KG ) #define MELECTRON_KG (9.1093837015e-31) #endif diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py new file mode 100644 index 00000000..a89a06c7 --- /dev/null +++ b/xfields/touschek/manager.py @@ -0,0 +1,305 @@ +import xtrack as xt +import xfields as xf + +import numpy as np +from scipy.integrate import quad +from scipy.special import i0 +from scipy.constants import physical_constants + +# =========================================== +# Constants +# =========================================== +ELECTRON_MASS_EV = xt.ELECTRON_MASS_EV +C_LIGHT_VACUUM = physical_constants['speed of light in vacuum'][0] +CLASSICAL_ELECTRON_RADIUS = physical_constants['classical electron radius'][0] + +class TouschekCalculator: + def __init__(self, manager): + self.manager = manager + self.twiss = None + + def _compute_piwinski_integral(self, tm, B1, B2): + km = np.arctan(np.sqrt(tm)) + + def int_piwinski(k, km, B1, B2): + t = np.tan(k) ** 2 + tm = np.tan(km) ** 2 + fact = ( + (2*t + 1)**2 * (t/tm / (1+t) - 1) / t + t - np.sqrt(t*tm * (1 + t)) + - (2 + 1 / (2*t)) * np.log(t/tm / (1+t)) + ) + if B2 * t < 500: + intp = fact * np.exp(-B1*t) * i0(B2*t) * np.sqrt(1+t) + else: + intp = ( + fact + * np.exp(B2*t - B1*t) + / np.sqrt(2*np.pi * B2*t) + * np.sqrt(1+t) + ) + return intp + + args = (km, B1, B2) + val, _ = quad( + int_piwinski, + km, + np.pi / 2, + args=args, + epsabs=1e-16, + epsrel=1e-12 + ) + + return val + + def _compute_piwinski_scattering_rate(self, element): + p0c = self.manager.particle_ref.p0c[0] + bunch_population = self.manager.bunch_population + gemitt_x = self.manager.gemitt_x + gemitt_y = self.manager.gemitt_y + alfx = self.twiss['alfx', element] + betx = self.twiss['betx', element] + alfy = self.twiss['alfy', element] + bety = self.twiss['bety', element] + sigma_z = self.manager.sigma_z + sigma_delta = self.manager.sigma_delta + delta = self.twiss['delta', element] + dx = self.twiss['dx', element] + dpx = self.twiss['dpx', element] + dxt = alfx * dx + betx * dpx # dxt: dx tilde + dy = self.twiss['dy', element] + dpy = self.twiss['dpy', element] + dyt = alfy * dy + bety * dpy # dyt: dy tilde + + deltaN = self.manager.momentum_aperture.at[element, "deltaN"] + deltaP = self.manager.momentum_aperture.at[element, "deltaP"] + + sigmab_x = np.sqrt(gemitt_x * betx) # Horizontal betatron beam size + sigma_x = np.sqrt(gemitt_x * betx + dx**2 * sigma_delta**2) # Horizontal beam size + + sigmab_y = np.sqrt(gemitt_y * bety) # Vertical betatron beam size + sigma_y = np.sqrt(gemitt_y * bety + dy**2 * sigma_delta**2) # Vertical beam size + + sigma_h = (sigma_delta**-2 + (dx**2 + dxt**2)/sigmab_x**2 + (dy**2 + dyt**2)/sigmab_y**2)**(-0.5) + + p = p0c * (1 + delta) + gamma = np.sqrt(1 + p**2 / ELECTRON_MASS_EV**2) + beta = np.sqrt(1 - gamma**-2) + + B1 = betx**2 / (2 * beta**2 * gamma**2 * sigmab_x**2) * (1 - sigma_h**2 * dxt**2 / sigmab_x**2) \ + + bety**2 / (2 * beta**2 * gamma**2 * sigmab_y**2) * (1 - sigma_h**2 * dyt**2 / sigmab_y**2) + + B2 = np.sqrt(B1**2 - betx**2 * bety**2 * sigma_h**2 / (beta**4 * gamma**4 * sigmab_x**4 * sigmab_y**4 * sigma_delta**2) \ + * (sigma_x**2 * sigma_y**2 - sigma_delta**4 * dx**2 * dy**2)) + + tmN = beta**2 * (deltaN**2) + tmP = beta**2 * (deltaP**2) + + piwinski_integralN = self._compute_piwinski_integral(tmN, B1, B2) + piwinski_integralP = self._compute_piwinski_integral(tmP, B1, B2) + + rateN = CLASSICAL_ELECTRON_RADIUS**2 * C_LIGHT_VACUUM * bunch_population**2 \ + / (8*np.pi * gamma**2 * sigma_z * np.sqrt(sigma_x**2 * sigma_y**2 - sigma_delta**4 * dx**2 * dy**2)) \ + * 2 * np.sqrt(np.pi * (B1**2 - B2**2)) * piwinski_integralN + + rateP = CLASSICAL_ELECTRON_RADIUS**2 * C_LIGHT_VACUUM * bunch_population**2 \ + / (8*np.pi * gamma**2 * sigma_z * np.sqrt(sigma_x**2 * sigma_y**2 - sigma_delta**4 * dx**2 * dy**2)) \ + * 2 * np.sqrt(np.pi * (B1**2 - B2**2)) * piwinski_integralP + + rate = (rateN + rateP) / 2 + + return rate + + def _compute_integrated_piwinski_rates(self): + """ + Integrate Piwinski Touschek scattering rate along s using trapezoidal rule and + store the per-TouschekElelement the integrated rate per bunch. + """ + line = self.manager.line + tab = line.get_table() + T_rev0 = float(self.twiss.T_rev0) + + # Indexes of the TouschekScatterings + ii_t = [ii for ii, nn in enumerate(tab.name[:-1]) if isinstance(line[nn], xf.TouschekScattering)] + + integrated = 0.0 + s0 = 0.0 + r0 = self._compute_piwinski_scattering_rate(tab.name[0]) + + s_before = s0 + rate_before = r0 + ii_current = 0 + + for ii, nn in enumerate(tab.name): + s = tab.rows[nn].s[0] + ds = s - s_before + if ds > 0.0: + rate = self._compute_piwinski_scattering_rate(nn) + integrated += 0.5 * (rate_before + rate) * ds + s_before = s + rate_before = rate + + if ii_current < len(ii_t) and ii == ii_t[ii_current]: + # divide by c and by T_rev0 -> per-bunch rate + integrated_piwinski_rate = integrated / C_LIGHT_VACUUM / T_rev0 + elem = line[nn] # xf.TouschekScattering + elem._configure(_integrated_piwinski_rate=integrated_piwinski_rate) + integrated = 0.0 + ii_current += 1 + if ii_current >= len(ii_t): + break + +class TouschekManager: + def __init__(self, line, momentum_aperture, nemitt_x=None, nemitt_y=None, + sigma_z=None, sigma_delta=None, bunch_population=None, + n_simulated=None, gemitt_x=None, gemitt_y=None, + momentum_aperture_scale=0.85, ignored_portion=0.01, + nx=3, ny=3, nz=3, **kwargs): + + # Input validation + if line is None: + raise ValueError("`line` is required.") + if not hasattr(line, "particle_ref"): + raise ValueError("`line` must have a `particle_ref`.") + if momentum_aperture is None: + raise ValueError("`momentum_aperture` is required.") + if sigma_z is None: + raise ValueError("`sigma_z` is required.") + if sigma_delta is None: + raise ValueError("`sigma_delta` is required.") + if bunch_population is None: + raise ValueError("`bunch_population` is required.") + if n_simulated is None: + raise ValueError("`n_simulated` is required.") + + # Momentum aperture validation + required_cols = {"s", "name", "deltaN", "deltaP"} + if not hasattr(momentum_aperture, "columns") or not hasattr(momentum_aperture, "__getitem__"): + raise TypeError("`momentum_aperture` must be a DataFrame-like object with columns " + "'s', 'name', 'deltaN', 'deltaP'.") + missing = required_cols - set(momentum_aperture.columns) + if missing: + raise ValueError(f"`momentum_aperture` missing columns: {sorted(missing)}") + + for col in ("s", "deltaN", "deltaP"): + try: + vals = momentum_aperture[col].astype(float) + except Exception: + raise TypeError(f"`{col}` column must be numeric (cannot coerce to float).") + if not vals.notna().all(): + bad = list(vals.index[~vals.notna()][:5]) + raise ValueError(f"`{col}` contains NaN at rows {bad}.") + if (abs(vals) == float("inf")).any(): + bad = list(vals.index[(abs(vals) == float("inf"))][:5]) + raise ValueError(f"`{col}` contains inf at rows {bad}.") + + self.line = line + self.particle_ref = line.particle_ref + + momentum_aperture = momentum_aperture.copy() + momentum_aperture['deltaN'] *= momentum_aperture_scale + momentum_aperture['deltaP'] *= momentum_aperture_scale + self.momentum_aperture= momentum_aperture.set_index("name") + + self.sigma_z = sigma_z + self.sigma_delta = sigma_delta + self.bunch_population = bunch_population + self.n_simulated = n_simulated + self.momentum_aperture_scale = momentum_aperture_scale + self.ignored_portion = ignored_portion + self.nx = nx + self.ny = ny + self.nz = nz + + # Emittance validation + nemitt_given = nemitt_x is not None and nemitt_y is not None + gemitt_given = gemitt_x is not None and gemitt_y is not None + + if nemitt_given and gemitt_given: + raise ValueError("Provide either normalized emittances (nemitt_x, nemitt_y) " + "OR geometric emittances (gemitt_x, gemitt_y), not both.") + if not (nemitt_given or gemitt_given): + raise ValueError("You must provide either both normalized emittances (nemitt_x, nemitt_y) " + "OR both geometric emittances (gemitt_x, gemitt_y).") + + if nemitt_given: + beta0 = line.particle_ref.beta0[0] + gamma0 = line.particle_ref.gamma0[0] + self.gemitt_x = nemitt_x / (beta0 * gamma0) + self.gemitt_y = nemitt_y / (beta0 * gamma0) + else: + self.gemitt_x = gemitt_x + self.gemitt_y = gemitt_y + + self.kwargs = kwargs + + self.touschek = TouschekCalculator(self) + + # Check that the line contains TouschekScatterings + tab = self.line.get_table() + try: + has = "TouschekScattering" in set(np.unique(tab.element_type)) + except Exception: + has = "TouschekScattering" in set(getattr(tab, "element_type", [])) + if not has: + raise ValueError("The line does not contain any TouschekScattering. " + "Please add them before initializing the TouschekManager.") + + + def initialise_touschek(self, element=None): + line = self.line + tab = line.get_table() + + twiss_method = self.kwargs.get("method", "6d") + twiss = self.line.twiss(method=twiss_method) + # Pass the twiss to the TouschekCalculator + self.touschek.twiss = twiss + + self.touschek._compute_integrated_piwinski_rates() + + # Helper to assign all fields to a single TouschekScattering + def _assign(nn: str): + alfx = twiss["alfx", nn]; betx = twiss["betx", nn] + alfy = twiss["alfy", nn]; bety = twiss["bety", nn] + dx = twiss["dx", nn]; dpx = twiss["dpx", nn] + dy = twiss["dy", nn]; dpy = twiss["dpy", nn] + dN = self.momentum_aperture.at[nn, "deltaN"] + dP = self.momentum_aperture.at[nn, "deltaP"] + + piwinski_rate = self.touschek._compute_piwinski_scattering_rate(nn) + + elem = line[nn] # xf.TouschekScattering + + elem._configure( + _particle_ref=self.particle_ref, + _bunch_population=self.bunch_population, + _gemitt_x=self.gemitt_x, + _gemitt_y=self.gemitt_y, + _alfx=alfx, _betx=betx, + _alfy=alfy, _bety=bety, + _dx=dx, _dpx=dpx, + _dy=dy, _dpy=dpy, + _deltaN=dN, _deltaP=dP, + _sigma_z=self.sigma_z, + _sigma_delta=self.sigma_delta, + _n_simulated=self.n_simulated, + _nx=self.nx, _ny=self.ny, _nz=self.nz, + _ignored_portion=self.ignored_portion, + piwinski_rate=piwinski_rate, + ) + + if element is None: + for nn in tab.name[:-1]: # Avoid the last tab.name which is _end_point + if isinstance(line[nn], xf.TouschekScattering): + _assign(nn) + else: + if not isinstance(element, str): + raise TypeError(f"`element` must be a string (got {type(element).__name__}).") + if element not in set(tab.name): + raise ValueError( + f"`element='{element}'` is not present in the line provided to the TouschekManager." + ) + if not isinstance(line[element], xf.TouschekScattering): + raise TypeError( + f"`line['{element}']` is not a TouschekScattering (got {type(line[element]).__name__})." + ) + _assign(element) \ No newline at end of file From d366006ca2bfde6e57562752179ab6112d2fc380 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Thu, 4 Sep 2025 17:12:18 +0200 Subject: [PATCH 03/37] Return particle object. --- xfields/beam_elements/touschek.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/xfields/beam_elements/touschek.py b/xfields/beam_elements/touschek.py index 2f5a4b9d..80e9eefa 100644 --- a/xfields/beam_elements/touschek.py +++ b/xfields/beam_elements/touschek.py @@ -192,10 +192,11 @@ def scatter(self): weight=weight_out[:n], s=getattr(self, '_s', 0.0)) - self.particles = part self.total_mc_rate = totalMCRate_out self.ignored_rate = self._ignored_portion * self.total_mc_rate + return part + def track(self, particles): super().track(particles) \ No newline at end of file From 43c8b87d8bd91444c01f999f4f6dd6a05b412410 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Thu, 4 Sep 2025 17:59:26 +0200 Subject: [PATCH 04/37] Return a scalar. --- xfields/beam_elements/touschek.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/xfields/beam_elements/touschek.py b/xfields/beam_elements/touschek.py index 80e9eefa..fb273868 100644 --- a/xfields/beam_elements/touschek.py +++ b/xfields/beam_elements/touschek.py @@ -192,11 +192,11 @@ def scatter(self): weight=weight_out[:n], s=getattr(self, '_s', 0.0)) - self.total_mc_rate = totalMCRate_out + self.total_mc_rate = totalMCRate_out[0] self.ignored_rate = self._ignored_portion * self.total_mc_rate return part - + def track(self, particles): super().track(particles) \ No newline at end of file From 830052090c7e43f99fa1bf7b5c70c4f4700a8edb Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Thu, 4 Sep 2025 18:33:32 +0200 Subject: [PATCH 05/37] Configure correct s-location. --- xfields/beam_elements/touschek.py | 2 +- xfields/touschek/manager.py | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/xfields/beam_elements/touschek.py b/xfields/beam_elements/touschek.py index fb273868..b5345833 100644 --- a/xfields/beam_elements/touschek.py +++ b/xfields/beam_elements/touschek.py @@ -133,7 +133,7 @@ def __init__(self, def _configure(self, **kwargs): config_allowed = { - "_particle_ref", "_bunch_population", + "_s", "_particle_ref", "_bunch_population", "_gemitt_x", "_gemitt_y", "_alfx", "_betx", "_alfy", "_bety", "_dx", "_dpx", "_dy", "_dpy", diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index a89a06c7..0803b0e4 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -258,6 +258,7 @@ def initialise_touschek(self, element=None): # Helper to assign all fields to a single TouschekScattering def _assign(nn: str): + s = tab.rows[nn].s[0] alfx = twiss["alfx", nn]; betx = twiss["betx", nn] alfy = twiss["alfy", nn]; bety = twiss["bety", nn] dx = twiss["dx", nn]; dpx = twiss["dpx", nn] @@ -270,6 +271,7 @@ def _assign(nn: str): elem = line[nn] # xf.TouschekScattering elem._configure( + _s=s, _particle_ref=self.particle_ref, _bunch_population=self.bunch_population, _gemitt_x=self.gemitt_x, From 0b6707038ed8597d737992af9a2781f84bc026fc Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Fri, 5 Sep 2025 19:58:10 +0200 Subject: [PATCH 06/37] Use ELEGANT's random number generator. --- xfields/beam_elements/touschek.py | 24 +- xfields/beam_elements/touschek_src/touschek.h | 62 +++-- xfields/headers/elegant_rng.h | 221 ++++++++++++++++++ xfields/touschek/manager.py | 12 +- 4 files changed, 265 insertions(+), 54 deletions(-) create mode 100644 xfields/headers/elegant_rng.h diff --git a/xfields/beam_elements/touschek.py b/xfields/beam_elements/touschek.py index b5345833..b81cddc8 100644 --- a/xfields/beam_elements/touschek.py +++ b/xfields/beam_elements/touschek.py @@ -2,19 +2,9 @@ import xtrack as xt import numpy as np -from scipy.integrate import quad -from scipy.special import i0 -from scipy.constants import physical_constants from ..general import _pkg_root -# =========================================== -# Constants -# =========================================== -ELECTRON_MASS_EV = xt.ELECTRON_MASS_EV -C_LIGHT_VACUUM = physical_constants['speed of light in vacuum'][0] -CLASSICAL_ELECTRON_RADIUS = physical_constants['classical electron radius'][0] - class TouschekScattering(xt.BeamElement): _xofields = { @@ -43,6 +33,8 @@ class TouschekScattering(xt.BeamElement): '_nz': xo.Float64, '_ignored_portion': xo.Float64, '_integrated_piwinski_rate': xo.Float64, + '_seed': xo.Int64, + '_inhibit_permute': xo.Int64 } # allow_track = False @@ -50,6 +42,7 @@ class TouschekScattering(xt.BeamElement): _extra_c_sources = [ _pkg_root.joinpath('headers/constants.h'), + _pkg_root.joinpath('headers/elegant_rng.h'), _pkg_root.joinpath('beam_elements/touschek_src/touschek.h'), ] @@ -79,8 +72,7 @@ class TouschekScattering(xt.BeamElement): ), } - def __init__(self, - s=0.0, + def __init__(self, s=0.0, particle_ref=xt.Particles(), bunch_population=0.0, alfx=0.0, betx=0.0, alfy=0.0, bety=0.0, @@ -92,6 +84,8 @@ def __init__(self, piwinski_rate=0.0, ignored_portion=0.0, integrated_piwinski_rate=0.0, + seed=1997, + inhibit_permute=0, **kwargs): # This gives AttributeError: 'TouschekScattering' object has no attribute '_xobject' @@ -127,8 +121,9 @@ def __init__(self, self._nz = nz self._ignored_portion = ignored_portion self._integrated_piwinski_rate = integrated_piwinski_rate - self.piwinski_rate = piwinski_rate + self._seed = seed + self._inhibit_permute = inhibit_permute def _configure(self, **kwargs): @@ -141,7 +136,8 @@ def _configure(self, **kwargs): "_sigma_z", "_sigma_delta", "_n_simulated", "_nx", "_ny", "_nz", "_ignored_portion", "piwinski_rate", - "_integrated_piwinski_rate" + "_integrated_piwinski_rate", + "_seed", "_inhibit_permute" } unknown = set(kwargs) - config_allowed diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h index b9090573..a2425d9b 100644 --- a/xfields/beam_elements/touschek_src/touschek.h +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -5,6 +5,7 @@ #include #include #include +#include static inline double sqr(double x){ return x*x; } @@ -14,19 +15,6 @@ void TouschekScattering_track_local_particle(TouschekScatteringData el, LocalPar return; } -/*gpufun*/ -void shuffle_double_array(double* a, int n, LocalParticle* part) { - for (int i = n - 1; i > 0; --i) { - // Option A: use floating RNG - int j = (int)floor(RandomUniformAccurate_generate(part) * (i + 1)); - - double tmp = a[i]; - a[i] = a[j]; - a[j] = tmp; - } -} - - void selectPartGauss(double *p1, double *p2, double *dens1, double *dens2, const double *ran1, @@ -312,6 +300,14 @@ void TouschekScatter(TouschekScatteringData el, double *zetatemp = (double*)malloc(sizeof(double) * n_simulated); double *deltatemp = (double*)malloc(sizeof(double) * n_simulated); + static int seeded_once = 0; + if (!seeded_once){ + long seed = TouschekScatteringData_get__seed(el); + short inhibit = (short)TouschekScatteringData_get__inhibit_permute(el); + seedElegantRandomNumbers(seed, inhibit); + seeded_once = 1; + } + i = 0; j = 0; total_event = 0; @@ -322,10 +318,23 @@ void TouschekScatter(TouschekScatteringData el, break; /* Select 11 random numbers, then mix them. */ + + // These 11 random numbers are assigned to: + // particle 1 (p1) as: { x, y, px, py, zeta, delta} + // particle 2 (p2) as: { -, -, px, py, -, delta} + // scattering angles in the cm frame: theta and phi + + // In ELEGANT the 11 random numbers are assigned to: + // particle 1 (p1) as: { x, y, xp, yp, zeta, delta} + // particle 2 (p2) as: { -, -, xp, yp, -, delta} + // scattering angles in the cm frame: theta and phi + + // NOTE: ELEGANT uses slopes xp=dx/ds, yp=dy/ds instead of the normalized momentum components px=Px/p0c, py=Py/p0c for (j = 0; j < 11; j++) { - ran1[j] = RandomUniformAccurate_generate(part0); // replaces random_1_elegant + // ran1[j] = RandomUniformAccurate_generate(part0); // Does not match with ELEGANT + ran1[j] = random_1_elegant(1); } - // shuffle_double_array(ran1, 11, part0); // may replace ELEGANT's randomizeOrder ? + randomizeOrder((char*)ran1, sizeof(ran1[0]), 11, 0, random_4); // like ELEGANT total_event++; @@ -334,7 +343,9 @@ void TouschekScatter(TouschekScatteringData el, if (!dens1 || !dens2) { continue; } - /* Change from slope to momentum */ + /* Here ELEGANT changes from slopes to momentum components */ + // Since we use already the normalized momentum components {px, py} instead of the slopes {xp, yp} + // here we just unormalize the momentum components: Px=px*p0c, Py=py*p0c for (j = 3; j < 5; j++) { p1[j] *= p0c; p2[j] *= p0c; @@ -401,25 +412,6 @@ void TouschekScatter(TouschekScatteringData el, i++; } } - // ELEGANT warnings - // if (total_event * 11 > (long)2e9) { - // printWarning("TouschekScatter: the total number of random numbers used exceeded 2e9.", - // "Use smaller n_simulated or smaller delta."); - // break; - // } - // } - // if (verbosity) - // report_stats(stdout, "After particle generation: "); - // if (simuCount == 0) - // bombElegant("It appears that the Touschek lifetime is extremely wrong and there is no need to perform Touschek simulation.\n If you think this is a wrong statement, send input to developers for evaluation.\n", NULL); - - // if (total_event / tsptr->simuCount > 20) { - // if (distribution_cutoff[0] < 5 || distribution_cutoff[1] < 5) - // printWarning("touschek_scatter: scattering rate is low.", - // "Please use >=5 sigma beam for better simulation."); - // else - // printWarning("touschek_scatter: Sscattering rate is very low.", - // "Please ignore the rate from Monte Carlo simulation. Use Piwinski's rate only."); } factor = factor / (double)(total_event); totalMCRate = totalWeight * factor; diff --git a/xfields/headers/elegant_rng.h b/xfields/headers/elegant_rng.h new file mode 100644 index 00000000..ca278a4a --- /dev/null +++ b/xfields/headers/elegant_rng.h @@ -0,0 +1,221 @@ +// elegant_rng.h (C99) + +#ifndef ELEGANT_RNG_H +#define ELEGANT_RNG_H + +#include +#include +#include +#include + +static inline double dlaran_core(int32_t iseed[4]) { + // Seed chunks: iseed[0..3], iseed[3] must be odd + int32_t it1, it2, it3, it4; + it4 = iseed[3] * 2549; + it3 = it4 / 4096; + it4 -= (it3 << 12); + it3 = it3 + iseed[2] * 2549 + iseed[3] * 2508; + it2 = it3 / 4096; + it3 -= (it2 << 12); + it2 = it2 + iseed[1] * 2549 + iseed[2] * 2508 + iseed[3] * 322; + it1 = it2 / 4096; + it2 -= (it1 << 12); + it1 = it1 + iseed[0] * 2549 + iseed[1] * 2508 + iseed[2] * 322 + iseed[3] * 494; + it1 %= 4096; + + iseed[0] = it1; + iseed[1] = it2; + iseed[2] = it3; + iseed[3] = it4; + + const double twoneg12 = 2.44140625e-4; // 2^-12 + return ( (double)it1 + + ( (double)it2 + ( (double)it3 + (double)it4 * twoneg12 ) * twoneg12 ) * twoneg12 + ) * twoneg12; // in (0,1) +} + +/* --------- seed permutation control ---------- */ +static short g_inhibitPermute = 0; +static inline short inhibitRandomSeedPermutation(short state){ + if (state >= 0) g_inhibitPermute = state; + return g_inhibitPermute; +} + +static inline uint32_t permuteSeedBitOrder(uint32_t input0){ + if (g_inhibitPermute) return input0; + uint32_t input = input0; + uint32_t newValue = 0u; + uint32_t offset = input0 % 1000u; + static const uint32_t bitMask[32] = { + 0x00000001u,0x00000002u,0x00000004u,0x00000008u, + 0x00000010u,0x00000020u,0x00000040u,0x00000080u, + 0x00000100u,0x00000200u,0x00000400u,0x00000800u, + 0x00001000u,0x00002000u,0x00004000u,0x00008000u, + 0x00010000u,0x00020000u,0x00040000u,0x00080000u, + 0x00100000u,0x00200000u,0x00400000u,0x00800000u, + 0x01000000u,0x02000000u,0x04000000u,0x08000000u, + 0x10000000u,0x20000000u,0x40000000u,0x80000000u + }; + for (int i=0;i<31;i++) + if (input & bitMask[i]) newValue |= bitMask[(i + offset) % 31]; + if (newValue == input){ + offset++; + newValue = 0u; + for (int i=0;i<31;i++) + if (input & bitMask[i]) newValue |= bitMask[(i + offset) % 31]; + } + return newValue; +} + +/* --------- RNG streams 1..6 (LAPACK DLARAN core) ---------- */ +static inline void seed_from_long(int32_t seed[4], long iseed_in, int force_odd_last){ + uint32_t s = (uint32_t)(iseed_in < 0 ? -iseed_in : iseed_in); + s = permuteSeedBitOrder(s); + // pack into 4x12-bit chunks, last must be odd + seed[3] = (int32_t)(s & 4095u); s >>= 12; + if (force_odd_last) seed[3] = (seed[3] | 1); // ensure odd + seed[2] = (int32_t)(s & 4095u); s >>= 12; + seed[1] = (int32_t)(s & 4095u); s >>= 12; + seed[0] = (int32_t)(s & 4095u); +} + +double random_2(long iseed); +double random_3(long iseed); +double random_4(long iseed); +double random_5(long iseed); +double random_6(long iseed); + +double random_1(long iseed){ + static int initialized = 0; + static int32_t seed[4] = {0,0,0,0}; + if (!initialized || iseed < 0){ + long base = (iseed < 0) ? -iseed : iseed; // abs + base = (long)permuteSeedBitOrder((uint32_t)base); // permute + // reseed degli altri stream come in SDDS (argomento NEGATIVO) + random_2(-(base + 2)); + random_3(-(base + 4)); + random_4(-(base + 6)); + random_5(-(base + 8)); + random_6(-(base + 10)); + // forza odd e spezza in 4×12 bit + base = (base/2)*2 + 1; + uint32_t s = (uint32_t)base; + seed[3] = (int32_t)(s & 4095u); s >>= 12; + seed[2] = (int32_t)(s & 4095u); s >>= 12; + seed[1] = (int32_t)(s & 4095u); s >>= 12; + seed[0] = (int32_t)(s & 4095u); + initialized = 1; + } + return dlaran_core(seed); +} + +double random_2(long iseed){ + static int initialized = 0; + static int32_t seed[4] = {0,0,0,0}; + if (!initialized || iseed < 0){ + if (iseed >= 0) iseed = -1; + seed_from_long(seed, -iseed, /*force_odd_last=*/1); + initialized = 1; + } + return dlaran_core(seed); +} +double random_3(long iseed){ + static int initialized = 0; + static int32_t seed[4] = {0,0,0,0}; + if (!initialized || iseed < 0){ + if (iseed >= 0) iseed = -1; + seed_from_long(seed, -iseed, /*force_odd_last=*/1); + initialized = 1; + } + return dlaran_core(seed); +} +double random_4(long iseed){ + static int initialized = 0; + static int32_t seed[4] = {0,0,0,0}; + if (!initialized || iseed < 0){ + if (iseed >= 0) iseed = -1; + seed_from_long(seed, -iseed, /*force_odd_last=*/1); + initialized = 1; + } + return dlaran_core(seed); +} +double random_5(long iseed){ + static int initialized = 0; + static int32_t seed[4] = {0,0,0,0}; + if (!initialized || iseed < 0){ + if (iseed >= 0) iseed = -1; + seed_from_long(seed, -iseed, /*force_odd_last=*/1); + initialized = 1; + } + return dlaran_core(seed); +} +double random_6(long iseed){ + static int initialized = 0; + static int32_t seed[4] = {0,0,0,0}; + if (!initialized || iseed < 0){ + if (iseed >= 0) iseed = -1; + seed_from_long(seed, -iseed, /*force_odd_last=*/1); + initialized = 1; + } + return dlaran_core(seed); +} + +/* comodo alias se nel tuo codice usi random_1_elegant */ +static inline double random_1_elegant(long iseed){ return random_1(iseed); } + +/* --------- randomizeOrder stile Elegant (qsort) ---------- */ +typedef struct RANDOMIZATION_HOLDER_ { + void* buffer; + double randomValue; +} RANDOMIZATION_HOLDER; + +static int randomizeOrderCmp(const void *p1, const void *p2) { + const RANDOMIZATION_HOLDER *rh1 = (const RANDOMIZATION_HOLDER *)p1; + const RANDOMIZATION_HOLDER *rh2 = (const RANDOMIZATION_HOLDER *)p2; + if (rh1->randomValue > rh2->randomValue) return 1; + if (rh1->randomValue < rh2->randomValue) return -1; + return 0; +} + +static long randomizeOrder(char *ptr, long size, long length, + long iseed, double (*urandom)(long iseed1)) { + if (!ptr || size<=0 || length<=1 || !urandom) return 0; + if (iseed < 0) urandom(iseed); + RANDOMIZATION_HOLDER *rh = (RANDOMIZATION_HOLDER*)malloc(sizeof(*rh)* (size_t)length); + if (!rh) return 0; + for (long i=0;i Date: Sun, 7 Sep 2025 09:43:47 +0200 Subject: [PATCH 07/37] Clean-up. --- xfields/beam_elements/touschek.py | 23 ++++++------------- xfields/touschek/manager.py | 37 ++++++++++++++++++++++++++----- 2 files changed, 38 insertions(+), 22 deletions(-) diff --git a/xfields/beam_elements/touschek.py b/xfields/beam_elements/touschek.py index b81cddc8..63ca3b47 100644 --- a/xfields/beam_elements/touschek.py +++ b/xfields/beam_elements/touschek.py @@ -1,6 +1,10 @@ +# copyright ################################# # +# This file is part of the Xfields Package. # +# Copyright (c) CERN, 2021. # +# ########################################### # + import xobjects as xo import xtrack as xt - import numpy as np from ..general import _pkg_root @@ -8,9 +12,6 @@ class TouschekScattering(xt.BeamElement): _xofields = { - # All these are handled by the manager for now - - # Everything that need to be accessed with TouschekElementData_get in C '_p0c': xo.Float64, '_bunch_population': xo.Float64, '_gemitt_x': xo.Float64, @@ -50,15 +51,6 @@ class TouschekScattering(xt.BeamElement): '_scatter': xo.Kernel( c_name='TouschekScatter', args=[ - # xo.Arg(xo.Int64, name='n_simulated'), - # xo.Arg(xo.Float64, name='nx'), - # xo.Arg(xo.Float64, name='ny'), - # xo.Arg(xo.Float64, name='nz'), - # xo.Arg(xo.Float64, name='ignored_portion'), - - # These are variables that we will get out from the C kernel - # We will need to pre-allocate memory for them in Python - # Then the C kernel will fill them in xo.Arg(xo.Float64, name='x_out', pointer=True), xo.Arg(xo.Float64, name='px_out', pointer=True), xo.Arg(xo.Float64, name='y_out', pointer=True), @@ -125,7 +117,6 @@ def __init__(self, s=0.0, self._seed = seed self._inhibit_permute = inhibit_permute - def _configure(self, **kwargs): config_allowed = { "_s", "_particle_ref", "_bunch_population", @@ -150,7 +141,6 @@ def _configure(self, **kwargs): if kk == "_particle_ref": self._p0c = self._particle_ref.p0c[0] - def scatter(self): context = self._context particles = xt.Particles(_context=context) @@ -177,6 +167,8 @@ def scatter(self): n_selected_out=n_selected_out) n = n_selected_out[0] + # Create particle object for tracking + # TODO: add at_element, start_tracking_at_element, ... part = xt.Particles(_capacity=2*n, p0c=self._p0c, mass0=self._particle_ref.mass0, @@ -193,6 +185,5 @@ def scatter(self): return part - def track(self, particles): super().track(particles) \ No newline at end of file diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index bec922cd..674faaeb 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -1,3 +1,8 @@ +# copyright ################################# # +# This file is part of the Xfields Package. # +# Copyright (c) CERN, 2021. # +# ########################################### # + import xtrack as xt import xfields as xf @@ -6,9 +11,6 @@ from scipy.special import i0 from scipy.constants import physical_constants -# =========================================== -# Constants -# =========================================== ELECTRON_MASS_EV = xt.ELECTRON_MASS_EV C_LIGHT_VACUUM = physical_constants['speed of light in vacuum'][0] CLASSICAL_ELECTRON_RADIUS = physical_constants['classical electron radius'][0] @@ -19,6 +21,15 @@ def __init__(self, manager): self.twiss = None def _compute_piwinski_integral(self, tm, B1, B2): + """ + Compute Piwinski integral for Touschek scattering rate calculation. + + Reference: + A. Piwinski, + "The Touschek Effect in Strong Focusing Storage Rings", + arXiv:physics/9903034, 1999. + URL: https://arxiv.org/abs/physics/9903034 + """ km = np.arctan(np.sqrt(tm)) def int_piwinski(k, km, B1, B2): @@ -52,6 +63,15 @@ def int_piwinski(k, km, B1, B2): return val def _compute_piwinski_scattering_rate(self, element): + """ + Compute Piwinski Touschek scattering rate. + + Reference: + A. Piwinski, + "The Touschek Effect in Strong Focusing Storage Rings", + arXiv:physics/9903034, 1999. + URL: https://arxiv.org/abs/physics/9903034 + """ p0c = self.manager.particle_ref.p0c[0] bunch_population = self.manager.bunch_population gemitt_x = self.manager.gemitt_x @@ -111,8 +131,13 @@ def _compute_piwinski_scattering_rate(self, element): def _compute_integrated_piwinski_rates(self): """ - Integrate Piwinski Touschek scattering rate along s using trapezoidal rule and - store the per-TouschekElelement the integrated rate per bunch. + Integrate the Piwinski Touschek scattering rate along the line using + the trapezoidal rule, between successive TouschekScattering elements. + + For each TouschekScattering element, the method stores the integrated + rate per bunch over the preceding section of the line. This per-bunch + rate is later used to assign the correct weights to Touschek-scattered + particles at the corresponding element. """ line = self.manager.line tab = line.get_table() @@ -139,7 +164,7 @@ def _compute_integrated_piwinski_rates(self): rate_before = rate if ii_current < len(ii_t) and ii == ii_t[ii_current]: - # divide by c and by T_rev0 -> per-bunch rate + # divide by c and by T_rev0 --> per-bunch rate integrated_piwinski_rate = integrated / C_LIGHT_VACUUM / T_rev0 elem = line[nn] # xf.TouschekScattering elem._configure(_integrated_piwinski_rate=integrated_piwinski_rate) From 204d05dcbdb4cf3282f932c682b656916d652b17 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 09:45:25 +0200 Subject: [PATCH 08/37] Include ELEGANT and SDDS credits. --- THIRD_PARTY_NOTICES | 3 + third_party/SDDS/LICENSE | 66 ++++++++ third_party/elegant/LICENSE | 66 ++++++++ xfields/beam_elements/touschek_src/touschek.h | 49 +++++- xfields/headers/elegant_rng.h | 156 ++++++++++++++++-- 5 files changed, 321 insertions(+), 19 deletions(-) create mode 100644 THIRD_PARTY_NOTICES create mode 100644 third_party/SDDS/LICENSE create mode 100644 third_party/elegant/LICENSE diff --git a/THIRD_PARTY_NOTICES b/THIRD_PARTY_NOTICES new file mode 100644 index 00000000..6ccdbad0 --- /dev/null +++ b/THIRD_PARTY_NOTICES @@ -0,0 +1,3 @@ +This package includes code adapted from Elegant and SDDS. +© 2002 The University of Chicago; © 2002 The Regents of the University of California. +Distributed subject to their Software License Agreements. See third_party/elegant/LICENSE and third_party/SDDS/LICENSE. \ No newline at end of file diff --git a/third_party/SDDS/LICENSE b/third_party/SDDS/LICENSE new file mode 100644 index 00000000..c14ef88f --- /dev/null +++ b/third_party/SDDS/LICENSE @@ -0,0 +1,66 @@ +Copyright (c) 2002 University of Chicago. All rights reserved. + +SDDS ToolKit is distributed subject to the following license conditions: + + SOFTWARE LICENSE AGREEMENT + Software: SDDS ToolKit + + 1. The "Software", below, refers to SDDS ToolKit (in either source code, or + binary form and accompanying documentation). Each licensee is + addressed as "you" or "Licensee." + + 2. The copyright holders shown above and their third-party licensors + hereby grant Licensee a royalty-free nonexclusive license, subject to + the limitations stated herein and U.S. Government license rights. + + 3. You may modify and make a copy or copies of the Software for use + within your organization, if you meet the following conditions: + a. Copies in source code must include the copyright notice and this + Software License Agreement. + b. Copies in binary form must include the copyright notice and this + Software License Agreement in the documentation and/or other + materials provided with the copy. + + 4. You may modify a copy or copies of the Software or any portion of it, + thus forming a work based on the Software, and distribute copies of + such work outside your organization, if you meet all of the following + conditions: + a. Copies in source code must include the copyright notice and this + Software License Agreement; + b. Copies in binary form must include the copyright notice and this + Software License Agreement in the documentation and/or other + materials provided with the copy; + c. Modified copies and works based on the Software must carry + prominent notices stating that you changed specified portions of + the Software. + + 5. Portions of the Software resulted from work developed under a U.S. + Government contract and are subject to the following license: the + Government is granted for itself and others acting on its behalf a + paid-up, nonexclusive, irrevocable worldwide license in this computer + software to reproduce, prepare derivative works, and perform publicly + and display publicly. + + 6. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" WITHOUT + WARRANTY OF ANY KIND. THE COPYRIGHT HOLDERS, THEIR THIRD PARTY + LICENSORS, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF + ENERGY, AND THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED + WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, + TITLE OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY + OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR USEFULNESS + OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF THE SOFTWARE + WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) DO NOT WARRANT + THAT THE SOFTWARE WILL FUNCTION UNINTERRUPTED, THAT IT IS + ERROR-FREE OR THAT ANY ERRORS WILL BE CORRECTED. + + 7. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT HOLDERS, + THEIR THIRD PARTY LICENSORS, THE UNITED STATES, THE UNITED + STATES DEPARTMENT OF ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR + ANY INDIRECT, INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE + DAMAGES OF ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS + OF PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER + SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT + (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, EVEN + IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE POSSIBILITY OF + SUCH LOSS OR DAMAGES. \ No newline at end of file diff --git a/third_party/elegant/LICENSE b/third_party/elegant/LICENSE new file mode 100644 index 00000000..bd9e10ea --- /dev/null +++ b/third_party/elegant/LICENSE @@ -0,0 +1,66 @@ +Copyright (c) 2002 University of Chicago. All rights reserved. + +elegant is distributed subject to the following license conditions: + + SOFTWARE LICENSE AGREEMENT + Software: elegant + + 1. The "Software", below, refers to elegant (in either source code, or + binary form and accompanying documentation). Each licensee is + addressed as "you" or "Licensee." + + 2. The copyright holders shown above and their third-party licensors + hereby grant Licensee a royalty-free nonexclusive license, subject to + the limitations stated herein and U.S. Government license rights. + + 3. You may modify and make a copy or copies of the Software for use + within your organization, if you meet the following conditions: + a. Copies in source code must include the copyright notice and this + Software License Agreement. + b. Copies in binary form must include the copyright notice and this + Software License Agreement in the documentation and/or other + materials provided with the copy. + + 4. You may modify a copy or copies of the Software or any portion of it, + thus forming a work based on the Software, and distribute copies of + such work outside your organization, if you meet all of the following + conditions: + a. Copies in source code must include the copyright notice and this + Software License Agreement; + b. Copies in binary form must include the copyright notice and this + Software License Agreement in the documentation and/or other + materials provided with the copy; + c. Modified copies and works based on the Software must carry + prominent notices stating that you changed specified portions of + the Software. + + 5. Portions of the Software resulted from work developed under a U.S. + Government contract and are subject to the following license: the + Government is granted for itself and others acting on its behalf a + paid-up, nonexclusive, irrevocable worldwide license in this computer + software to reproduce, prepare derivative works, and perform publicly + and display publicly. + + 6. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" WITHOUT + WARRANTY OF ANY KIND. THE COPYRIGHT HOLDERS, THEIR THIRD PARTY + LICENSORS, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF + ENERGY, AND THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED + WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, + TITLE OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY + OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR USEFULNESS + OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF THE SOFTWARE + WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) DO NOT WARRANT + THAT THE SOFTWARE WILL FUNCTION UNINTERRUPTED, THAT IT IS + ERROR-FREE OR THAT ANY ERRORS WILL BE CORRECTED. + + 7. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT HOLDERS, + THEIR THIRD PARTY LICENSORS, THE UNITED STATES, THE UNITED + STATES DEPARTMENT OF ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR + ANY INDIRECT, INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE + DAMAGES OF ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS + OF PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER + SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT + (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, EVEN + IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE POSSIBILITY OF + SUCH LOSS OR DAMAGES. \ No newline at end of file diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h index a2425d9b..3b977b8a 100644 --- a/xfields/beam_elements/touschek_src/touschek.h +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -1,3 +1,43 @@ +/* touschek.h — Touschek scattering routine (C99, header-only kernel) + + Portions adapted from Elegant/SDDS. + + Original notice (preserved as required): + ---------------------------------------------------------------------- + Copyright (c) 2002 The University of Chicago, as Operator of Argonne + National Laboratory. + Copyright (c) 2002 The Regents of the University of California, as + Operator of Los Alamos National Laboratory. + This file is distributed subject to a Software License Agreement found + in the file LICENSE that is included with this distribution. + ---------------------------------------------------------------------- + + This derivative file is distributed with the same notice; see + third_party/elegant/LICENSE (Elegant) + third_party/sdds/LICENSE (SDDS) + included in this source tree. + + Modifications (c) 2025 Giacomo Broggi / CERN. + Prominent changes from Elegant’s `touschekScatter.c`: + - Converted to a header-only C99 kernel and simplified API (no SDDS I/O). + Which has been made compatible with `xobjects` via the `xobjects` API. + - Uses `elegant_rng.h` for RNG with Elegant-identical streams: + draws via `random_1_elegant` and shuffling via `random_4` + `randomizeOrder`, + matching Elegant’s RNG consumption. + - Works in terms of normalized momentum (px,py) and then un-normalizes to eV, + documenting the slope (xp,yp) vs momentum difference used in Elegant. + - Small safety/cleanup changes (bounds checks, allocations, comments). + - Kept physics and selection logic identical. + + Attribution / citation: + If you publish results produced with this routines, please also cite: + M. Borland, “elegant: A Flexible SDDS-Compliant Code for Accelerator Simulation,” + Advanced Photon Source LS-287, September 2000. + + SPDX (license identifiers for scanners): + SPDX-License-Identifier: LicenseRef-ELEGANT + SPDX-License-Identifier: LicenseRef-SDDS +*/ #ifndef XTRACK_TOUSCHEK_H #define XTRACK_TOUSCHEK_H @@ -15,6 +55,7 @@ void TouschekScattering_track_local_particle(TouschekScatteringData el, LocalPar return; } +/* Adapted from Elegant touschekScatter.c: selectPartGauss (logic unchanged). */ void selectPartGauss(double *p1, double *p2, double *dens1, double *dens2, const double *ran1, @@ -66,6 +107,7 @@ void selectPartGauss(double *p1, double *p2, return; } +/* From Elegant touschekScatter.c: bunch2cm */ void bunch2cm(double *p1, double *p2, double *q, double *beta, double *gamma) { double pp1, pp2, e1, e2, ee; int i; @@ -99,7 +141,9 @@ void bunch2cm(double *p1, double *p2, double *q, double *beta, double *gamma) { return; } + /* Rotate scattered p in c.o.m system */ +/* From Elegant touschekScatter.c: eulertrans*/ void eulertrans(double *v0, double theta, double phi, double *v1, double *v) { double th, ph, s1, s2, c1, c2; double x0, y0, z0; @@ -124,6 +168,7 @@ void eulertrans(double *v0, double theta, double phi, double *v1, double *v) { return; } +/* From Elegant touschekScatter.c: cm2bunch*/ void cm2bunch(double *p1, double *p2, double *q, double *beta, double *gamma) { int i; double pq, e, betaq, bb, factor; @@ -151,6 +196,7 @@ void cm2bunch(double *p1, double *p2, double *q, double *beta, double *gamma) { return; } +/* From Elegant touschekScatter.c: moeller */ double moeller(double beta0, double theta) { double cross; double beta2, st2; @@ -163,6 +209,7 @@ double moeller(double beta0, double theta) { return cross; } +/* From Elegant touschekScatter.c: pickPart */ void pickPart(double *weight, long *index, long start, long end, long *iTotal, double *wTotal, double weight_limit, double weight_ave) { long i, i1, i2, N; @@ -226,7 +273,7 @@ void pickPart(double *weight, long *index, long start, long end, return; } -/*gpufun*/ +/* Adapted from Elegant touschekScatter.c: TouschekDistribution (logic unchanged) */ void TouschekScatter(TouschekScatteringData el, LocalParticle* part0, double* x_out, diff --git a/xfields/headers/elegant_rng.h b/xfields/headers/elegant_rng.h index ca278a4a..1c41dfde 100644 --- a/xfields/headers/elegant_rng.h +++ b/xfields/headers/elegant_rng.h @@ -1,4 +1,49 @@ // elegant_rng.h (C99) +/* + * Portions adapted from Elegant/SDDS. + * Copyright (c) 2002 University of Chicago. All rights reserved. + * Modifications: C99 refactor, static-state consolidation, MPI modes omitted, + * randomizeOrder safety checks, bit-mask fix, and comments. + * (c) 2025 . All rights reserved. + */ +// +// This file mirrors algorithms and seeding conventions from: +// - Elegant: M. Borland, “elegant: A Flexible SDDS-Compliant Code for Accelerator Simulation,” +// Advanced Photon Source LS-287, September 2000. +// - SDDS: see the SDDS and Elegant source distributions and their LICENSE files. +// +// Please retain upstream copyright and license notices where code/logic +// has been adapted, and also cite LS-287 when publishing results produced with this RNG. +// +// ----------------------------------------------------------------------------- +// Purpose +// ------- +// Re-implements the random-number utilities used by Elegant/SDDS so that +// C-kernels compiled via xobjects can reproduce *bitwise-identical* +// sequences to Elegant. This file exposes: +// +// - A LAPACK-compatible 48-bit LCG core (DLARAN) working on 4×12-bit chunks +// - Six RNG streams random_1..random_6 with Elegant-compatible seeding rules +// - Seed-bit permutation with Elegant’s “inhibit” switch and special seed +// behavior (987654321) +// - randomizeOrder() that shuffles buffers using the same qsort+random-key +// scheme used by Elegant to match sequence consumption +// +// Notes +// ----- +// * Based on the DLARAN routine from LAPACK and the RNG/seeding conventions from +// SDDS/Elegant (e.g., drand.c and friends in the Elegant/SDDS sources). +// * This is a single-threaded, static-state implementation. +// If you call these from multiple threads, guard access yourself. +// * MPI-related seed diversification is (for now) intentionally omitted here. +// * “Negative seed” calls *reinitialize* the stream; “non-negative” consume. +// +// Attribution +// ----------- +// - DLARAN: LAPACK auxiliary RNG (48-bit LCG with multiplier 33952834046453). +// - Seeding conventions and the 987654321 “inhibit permutation” behavior mirror +// SDDS/Elegant (random_1..random_6, randomizeOrder). +// ----------------------------------------------------------------------------- #ifndef ELEGANT_RNG_H #define ELEGANT_RNG_H @@ -8,6 +53,14 @@ #include #include +// ----------------------------------------------------------------------------- +// LAPACK-style DLARAN core +// ------------------------ +// This is the 48-bit multiplicative LCG used by DLARAN, with the 48-bit state +// stored in 4 integers of 12 bits each. The last chunk (iseed[3]) must be odd. +// Returns a uniform double in (0,1). The recurrence and normalization constants +// reproduce LAPACK DLARAN bit-for-bit. +// ----------------------------------------------------------------------------- static inline double dlaran_core(int32_t iseed[4]) { // Seed chunks: iseed[0..3], iseed[3] must be odd int32_t it1, it2, it3, it4; @@ -34,7 +87,16 @@ static inline double dlaran_core(int32_t iseed[4]) { ) * twoneg12; // in (0,1) } -/* --------- seed permutation control ---------- */ +/* ----------------------------------------------------------------------------- + Seed permutation control (Elegant-compatible) + -------------------------------------------- + Elegant permutes the bit order of integer seeds before packing 12-bit chunks. + This helps decorrelate “close” seeds. A global flag can inhibit this step. + + - inhibitRandomSeedPermutation(state>=0) sets the global flag. + - permuteSeedBitOrder(x) permutes bit positions unless inhibited. + - Elegant disables permutation when random_number_seed == 987654321. +----------------------------------------------------------------------------- */ static short g_inhibitPermute = 0; static inline short inhibitRandomSeedPermutation(short state){ if (state >= 0) g_inhibitPermute = state; @@ -67,7 +129,10 @@ static inline uint32_t permuteSeedBitOrder(uint32_t input0){ return newValue; } -/* --------- RNG streams 1..6 (LAPACK DLARAN core) ---------- */ +/* ----------------------------------------------------------------------------- + Packing helper: split a 32-bit integer into 4×12-bit chunks (DLARAN format). + The last chunk must be odd for DLARAN to work correctly (Elegant behavior). +----------------------------------------------------------------------------- */ static inline void seed_from_long(int32_t seed[4], long iseed_in, int force_odd_last){ uint32_t s = (uint32_t)(iseed_in < 0 ? -iseed_in : iseed_in); s = permuteSeedBitOrder(s); @@ -79,6 +144,20 @@ static inline void seed_from_long(int32_t seed[4], long iseed_in, int force_odd_ seed[0] = (int32_t)(s & 4095u); } +/* ----------------------------------------------------------------------------- + RNG streams random_1 .. random_6 + -------------------------------- + These reproduce the SDDS/Elegant API and seeding semantics: + + - Calling with a *negative* iseed re-initializes that stream from |iseed|. + - Calling with a non-negative iseed consumes the next variate. + - random_1 (on (re)seed) also re-seeds streams 2..6 using |base|+{2,4,6,8,10}. + - All streams use the same DLARAN core, independent 48-bit states. + + Important: + * This file intentionally omits MPI diversification (modes 1..4 in Elegant). + * Keep static state: not thread-safe by design (matches Elegant). +----------------------------------------------------------------------------- */ double random_2(long iseed); double random_3(long iseed); double random_4(long iseed); @@ -91,13 +170,13 @@ double random_1(long iseed){ if (!initialized || iseed < 0){ long base = (iseed < 0) ? -iseed : iseed; // abs base = (long)permuteSeedBitOrder((uint32_t)base); // permute - // reseed degli altri stream come in SDDS (argomento NEGATIVO) + // SDDS-like reseed random_2(-(base + 2)); random_3(-(base + 4)); random_4(-(base + 6)); random_5(-(base + 8)); random_6(-(base + 10)); - // forza odd e spezza in 4×12 bit + // force odd and split in 4×12 bit base = (base/2)*2 + 1; uint32_t s = (uint32_t)base; seed[3] = (int32_t)(s & 4095u); s >>= 12; @@ -160,10 +239,35 @@ double random_6(long iseed){ return dlaran_core(seed); } -/* comodo alias se nel tuo codice usi random_1_elegant */ -static inline double random_1_elegant(long iseed){ return random_1(iseed); } +static inline double random_1_elegant(long iseed) { + static int initialized = 0; + static int32_t seed[4] = {0,0,0,0}; + if (!initialized || iseed < 0){ + long base = (iseed < 0) ? -iseed : iseed; // abs + // Respect current inhibit flag (seedElegantRandomNumbers sets it) + base = (long)permuteSeedBitOrder((uint32_t)base); // no-op if inhibited + base = (base/2)*2 + 1; // force odd + uint32_t s = (uint32_t)base; + seed[3] = (int32_t)(s & 4095u); s >>= 12; + seed[2] = (int32_t)(s & 4095u); s >>= 12; + seed[1] = (int32_t)(s & 4095u); s >>= 12; + seed[0] = (int32_t)(s & 4095u); + initialized = 1; + } + return dlaran_core(seed); +} + +/* ----------------------------------------------------------------------------- + randomizeOrder + ------------------------------ + Elegant shuffles arrays by: + 1) creating an array of (buffer copy, random key) pairs, + 2) sorting by the random key (qsort), + 3) copying buffers back in sorted order. -/* --------- randomizeOrder stile Elegant (qsort) ---------- */ + This consumes RNG like Elegant does and matches its order exactly. + It allocates O(N*size) memory. +----------------------------------------------------------------------------- */ typedef struct RANDOMIZATION_HOLDER_ { void* buffer; double randomValue; @@ -179,22 +283,28 @@ static int randomizeOrderCmp(const void *p1, const void *p2) { static long randomizeOrder(char *ptr, long size, long length, long iseed, double (*urandom)(long iseed1)) { - if (!ptr || size<=0 || length<=1 || !urandom) return 0; + if (!ptr || size<=0 || !urandom) return 0; + if (length < 2) return 1; if (iseed < 0) urandom(iseed); - RANDOMIZATION_HOLDER *rh = (RANDOMIZATION_HOLDER*)malloc(sizeof(*rh)* (size_t)length); + + RANDOMIZATION_HOLDER *rh = + (RANDOMIZATION_HOLDER*)malloc(sizeof(*rh) * (size_t)length); if (!rh) return 0; - for (long i=0;i inhibit permutation globally. + + Call this exactly once (per process) before any RNG use to mimic Elegant’s + &run_setup random_number_seed behavior. +----------------------------------------------------------------------------- */ static inline void seedElegantRandomNumbers(long seed, short inhibit_permute){ - long s0 = labs(seed); - long s1 = labs(seed + 2); - long s2 = labs(seed + 4); - long s3 = labs(seed + 6); + long s0 = labs(seed), s1 = labs(seed + 2), s2 = labs(seed + 4), s3 = labs(seed + 6); - /* ELEGANT: if seed==987654321 inhibit seed permutation */ if (s0 == 987654321) inhibitRandomSeedPermutation(1); else inhibitRandomSeedPermutation(inhibit_permute ? 1 : 0); From c828fbe3e7b2d4bc111b058b0b95ff7d2a8cddef Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 09:51:40 +0200 Subject: [PATCH 09/37] Move third party notice to `contributors_and_credits.txt`. --- THIRD_PARTY_NOTICES | 3 --- contributors_and_credits.txt | 4 ++++ 2 files changed, 4 insertions(+), 3 deletions(-) delete mode 100644 THIRD_PARTY_NOTICES diff --git a/THIRD_PARTY_NOTICES b/THIRD_PARTY_NOTICES deleted file mode 100644 index 6ccdbad0..00000000 --- a/THIRD_PARTY_NOTICES +++ /dev/null @@ -1,3 +0,0 @@ -This package includes code adapted from Elegant and SDDS. -© 2002 The University of Chicago; © 2002 The Regents of the University of California. -Distributed subject to their Software License Agreements. See third_party/elegant/LICENSE and third_party/SDDS/LICENSE. \ No newline at end of file diff --git a/contributors_and_credits.txt b/contributors_and_credits.txt index bf72b891..64044355 100644 --- a/contributors_and_credits.txt +++ b/contributors_and_credits.txt @@ -18,3 +18,7 @@ External contributors: The Particle In Cell implementation is largely based on the PyPIC package (https://github.com/pycomplete/pypic) + +The Touschek scattering routine in this package contains code adapted from Elegant and SDDS. +© 2002 The University of Chicago; © 2002 The Regents of the University of California. +Distributed subject to their Software License Agreements. See third_party/elegant/LICENSE and third_party/SDDS/LICENSE. From 0cbc05336005fc1ec5708c18bf6f06b1adac4300 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 09:52:41 +0200 Subject: [PATCH 10/37] Add myself as contributor. --- contributors_and_credits.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/contributors_and_credits.txt b/contributors_and_credits.txt index 64044355..e8253569 100644 --- a/contributors_and_credits.txt +++ b/contributors_and_credits.txt @@ -2,6 +2,7 @@ The following people contributed to the development of this package: CERN contributors: - H. Bartosik + - G. Broggi - X. Buffat - R. De Maria - L. Giacomel From ea07b3c0e27a3742e9a153b7e520e6c15db7e74c Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 10:16:06 +0200 Subject: [PATCH 11/37] Add README for touschek. --- xfields/beam_elements/touschek_src/touschek.h | 4 ++-- xfields/touschek/README.md | 22 +++++++++++++++++++ 2 files changed, 24 insertions(+), 2 deletions(-) create mode 100644 xfields/touschek/README.md diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h index 3b977b8a..fb468780 100644 --- a/xfields/beam_elements/touschek_src/touschek.h +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -18,7 +18,7 @@ included in this source tree. Modifications (c) 2025 Giacomo Broggi / CERN. - Prominent changes from Elegant’s `touschekScatter.c`: + Changes from Elegant’s `touschekScatter.c`: - Converted to a header-only C99 kernel and simplified API (no SDDS I/O). Which has been made compatible with `xobjects` via the `xobjects` API. - Uses `elegant_rng.h` for RNG with Elegant-identical streams: @@ -30,7 +30,7 @@ - Kept physics and selection logic identical. Attribution / citation: - If you publish results produced with this routines, please also cite: + If you publish results produced with this routine, please also cite: M. Borland, “elegant: A Flexible SDDS-Compliant Code for Accelerator Simulation,” Advanced Photon Source LS-287, September 2000. diff --git a/xfields/touschek/README.md b/xfields/touschek/README.md new file mode 100644 index 00000000..a7ff3b55 --- /dev/null +++ b/xfields/touschek/README.md @@ -0,0 +1,22 @@ +# Monte Carlo simulation of Touschek scattering + +The Monte-Carlo Touschek scattering routine in this package is based on the work by A. Xiao and M. Borland developed for **ELEGANT**. + +If you publish results obtained with this routine, please cite: + +- M. Borland, “elegant: A Flexible SDDS-Compliant Code for Accelerator Simulation,” APS LS-287 (2000). + +- A. Xiao and M. Borland, “Monte Carlo simulation of Touschek effect,” *Phys. Rev. ST Accel. Beams* **13**, 074201 (2010). + DOI: 10.1103/PhysRevSTAB.13.074201 + +## RNG compatibility with ELEGANT + +For reproducibility with ELEGANT, this routine uses the same RNG conventions; see `xfields/xfields/headers/elegant_rng.h`. + +## Third-party notice and licenses + +This routine contains code adapted from **ELEGANT** and **SDDS**. +© 2002 The University of Chicago; © 2002 The Regents of the University of California. +Distributed subject to their Software License Agreements. See: +- `third_party/elegant/LICENSE` +- `third_party/SDDS/LICENSE` \ No newline at end of file From 7069a69d3ab519f95a9a4fea52cb92a6355d1a7f Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 10:21:11 +0200 Subject: [PATCH 12/37] Add link to ELEGANT. --- xfields/touschek/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xfields/touschek/README.md b/xfields/touschek/README.md index a7ff3b55..97c046fe 100644 --- a/xfields/touschek/README.md +++ b/xfields/touschek/README.md @@ -1,6 +1,6 @@ # Monte Carlo simulation of Touschek scattering -The Monte-Carlo Touschek scattering routine in this package is based on the work by A. Xiao and M. Borland developed for **ELEGANT**. +The Monte-Carlo Touschek scattering routine in this package is based on the work by A. Xiao and M. Borland developed for [ELEGANT](https://github.com/rtsoliday/elegant). If you publish results obtained with this routine, please cite: From 30cffd1893ec9559345f77d99180a66bf395dd7b Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 10:22:36 +0200 Subject: [PATCH 13/37] Fix. --- xfields/touschek/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xfields/touschek/README.md b/xfields/touschek/README.md index 97c046fe..43e3e81c 100644 --- a/xfields/touschek/README.md +++ b/xfields/touschek/README.md @@ -1,6 +1,6 @@ # Monte Carlo simulation of Touschek scattering -The Monte-Carlo Touschek scattering routine in this package is based on the work by A. Xiao and M. Borland developed for [ELEGANT](https://github.com/rtsoliday/elegant). +The Monte-Carlo Touschek scattering routine in this package is based on the implementation by A. Xiao and M. Borland developed for [ELEGANT](https://github.com/rtsoliday/elegant). If you publish results obtained with this routine, please cite: From 690393633de27a6006c14a041003d44939667695 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 12:28:12 +0200 Subject: [PATCH 14/37] Cleanup headers. --- xfields/beam_elements/touschek_src/touschek.h | 109 +++++++++++------- xfields/headers/elegant_rng.h | 99 ++++++++-------- 2 files changed, 119 insertions(+), 89 deletions(-) diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h index fb468780..4e4c2bcc 100644 --- a/xfields/beam_elements/touschek_src/touschek.h +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -1,43 +1,58 @@ -/* touschek.h — Touschek scattering routine (C99, header-only kernel) - - Portions adapted from Elegant/SDDS. - - Original notice (preserved as required): - ---------------------------------------------------------------------- - Copyright (c) 2002 The University of Chicago, as Operator of Argonne - National Laboratory. - Copyright (c) 2002 The Regents of the University of California, as - Operator of Los Alamos National Laboratory. - This file is distributed subject to a Software License Agreement found - in the file LICENSE that is included with this distribution. - ---------------------------------------------------------------------- - - This derivative file is distributed with the same notice; see - third_party/elegant/LICENSE (Elegant) - third_party/sdds/LICENSE (SDDS) - included in this source tree. - - Modifications (c) 2025 Giacomo Broggi / CERN. - Changes from Elegant’s `touschekScatter.c`: - - Converted to a header-only C99 kernel and simplified API (no SDDS I/O). - Which has been made compatible with `xobjects` via the `xobjects` API. - - Uses `elegant_rng.h` for RNG with Elegant-identical streams: - draws via `random_1_elegant` and shuffling via `random_4` + `randomizeOrder`, - matching Elegant’s RNG consumption. - - Works in terms of normalized momentum (px,py) and then un-normalizes to eV, - documenting the slope (xp,yp) vs momentum difference used in Elegant. - - Small safety/cleanup changes (bounds checks, allocations, comments). - - Kept physics and selection logic identical. - - Attribution / citation: - If you publish results produced with this routine, please also cite: - M. Borland, “elegant: A Flexible SDDS-Compliant Code for Accelerator Simulation,” - Advanced Photon Source LS-287, September 2000. - - SPDX (license identifiers for scanners): - SPDX-License-Identifier: LicenseRef-ELEGANT - SPDX-License-Identifier: LicenseRef-SDDS -*/ +/*************************************************************************\ +* Copyright (c) 2002 The University of Chicago, as Operator of Argonne +* National Laboratory. +* Copyright (c) 2002 The Regents of the University of California, as +* Operator of Los Alamos National Laboratory. +* This file is distributed subject to a Software License Agreement found +* in the file LICENSE that is included with this distribution. +\*************************************************************************/ +/* + * touschek.h — Touschek scattering kernel (C99, header-only) + * + * Overview + * -------- + * Header-only C99 implementation of the Monte Carlo Touschek scattering + * routine for Xsuite/xfields. The physics and selection logic follow the + * Elegant implementation by A. Xiao and M. Borland (PRSTAB 13, 074201, 2010), + * with a simplified API suitable for xobjects; SDDS I/O is omitted. + * + * Provenance (portions adapted from) + * ---------------------------------- + * - Elegant: src/touschekScatter.c (A. Xiao, M. Borland) + * - SDDS Toolkit utilities: mdbmth/drand.c, mdbmth/dlaran.c + * - LAPACK DLARAN (48-bit LCG RNG core) + * + * Licenses & Notices + * ------------------ + * - Upstream notice preserved above (Elegant/SDDS). + * - See license texts in: + * xfields/third_party/elegant/LICENSE + * xfields/third_party/SDDS/LICENSE + * xfields/third_party/lapack/LICENSE + * - This file is a derivative work. + * Modifications © 2025 Giacomo Broggi / CERN. + * + * RNG compatibility + * ----------------- + * RNG streams are chosen to reproduce Elegant’s sequences. + * See: xfields/xfields/headers/elegant_rng.h + * (uses random_1_elegant, random_4, randomizeOrder, and DLARAN). + * + * Summary of modifications vs Elegant + * ----------------------------------- + * - Converted to a header-only C99 kernel; simplified API; no SDDS I/O. + * - Integrated with xobjects; clarified slope (xp, yp) vs momentum usage. + * - Minor safety/cleanup (bounds checks, allocations, comments). + * - Physics and selection logic preserved. + * + * Citation + * -------- + * If you publish results obtained with this routine, please cite: + * - M. Borland, “elegant: A Flexible SDDS-Compliant Code for Accelerator + * Simulation,” APS LS-287 (2000). + * - A. Xiao and M. Borland, “Monte Carlo simulation of Touschek effect,” + * Phys. Rev. ST Accel. Beams 13, 074201 (2010). DOI: 10.1103/PhysRevSTAB.13.074201 + */ #ifndef XTRACK_TOUSCHEK_H #define XTRACK_TOUSCHEK_H @@ -397,8 +412,11 @@ void TouschekScatter(TouschekScatteringData el, p1[j] *= p0c; p2[j] *= p0c; } - p1[5] = (p1[5] + 1) * p0c; - p2[5] = (p2[5] + 1) * p0c; + // p1[5] = (p1[5] + 1) * p0c; + // p2[5] = (p2[5] + 1) * p0c; + // Use exact formula to compute p1[5]=Pz1 and p2[5]=Pz2 + p1[5] = sqrt(sqr(p0c)*sqr(1. + p1[5]) - sqr(p1[3]) - sqr(p1[4])); + p2[5] = sqrt(sqr(p0c)*sqr(1. + p2[5]) - sqr(p2[3]) - sqr(p2[4])); bunch2cm(p1, p2, qa, beta, &gamma); @@ -408,8 +426,11 @@ void TouschekScatter(TouschekScatteringData el, temp = dens1 * dens2 * sin(theta); eulertrans(qa, theta, phi, qb, &qabs); cm2bunch(p1, p2, qb, beta, &gamma); - p1[5] = (p1[5] - p0c) / p0c; - p2[5] = (p2[5] - p0c) / p0c; + // p1[5] = (p1[5] - p0c) / p0c; + // p2[5] = (p2[5] - p0c) / p0c; + // Use exact formula to compute p1[5]=delta1 and p2[5]=delta2 + p1[5] = (sqrt(sqr(p1[3]) + sqr(p1[4]) + sqr(p1[5])) - p0c) / p0c; + p2[5] = (sqrt(sqr(p2[3]) + sqr(p2[4]) + sqr(p2[5])) - p0c) / p0c; if (p1[5] > p2[5]) { for (j = 0; j < 6; j++) { diff --git a/xfields/headers/elegant_rng.h b/xfields/headers/elegant_rng.h index 1c41dfde..02a141d9 100644 --- a/xfields/headers/elegant_rng.h +++ b/xfields/headers/elegant_rng.h @@ -1,50 +1,59 @@ -// elegant_rng.h (C99) +/*************************************************************************\ +* Portions adapted from Elegant and the SDDS Toolkit. +* Copyright (c) 2002 The University of Chicago. +* Copyright (c) 2002 The Regents of the University of California. +* This file is distributed subject to a Software License Agreement found +* in the file LICENSE that is included with this distribution. +\*************************************************************************/ /* - * Portions adapted from Elegant/SDDS. - * Copyright (c) 2002 University of Chicago. All rights reserved. - * Modifications: C99 refactor, static-state consolidation, MPI modes omitted, - * randomizeOrder safety checks, bit-mask fix, and comments. - * (c) 2025 . All rights reserved. + * elegant_rng.h — Elegant-compatible RNG utilities (C99) + * + * Overview + * -------- + * Re-implements the random-number utilities used by Elegant/SDDS so that + * kernels compiled via xobjects can reproduce (where applicable) the same + * sequences in Xsuite/xfields. Includes the LAPACK DLARAN core (48-bit LCG) + * and the streams random_1..random_6 with Elegant-compatible seeding rules. + * Provides random_1_elegant, randomizeOrder, and related helpers. + * + * Provenance (portions adapted from) + * ---------------------------------- + * - SDDS: mdbmth/drand.c (random_* streams, randomizeOrder, seeding) + * - SDDS: mdbmth/dlaran.c (C translation of LAPACK's DLARAN, via f2c) + * - Elegant: src/drand_oag.c (random_1_elegant and seed behavior) + * - LAPACK: DLARAN (48-bit LCG RNG core) + * + * Licenses & Notices + * ------------------ + * - Upstream notice preserved above (Elegant/SDDS). + * - License texts: + * xfields/third_party/elegant/LICENSE + * xfields/third_party/SDDS/LICENSE + * xfields/third_party/lapack/LICENSE + * - This file is a derivative work. + * Modifications © 2025 Giacomo Broggi / CERN. + * + * Purpose / Exposed API + * --------------------- + * - LAPACK-compatible DLARAN core (48-bit, 4×12-bit seed) + * - Streams: random_1 .. random_6 (Elegant-style seeding conventions) + * - random_1_elegant() and seedElegantRandomNumbers() (for Elegant compatibility) + * - permuteSeedBitOrder(), inhibitRandomSeedPermutation() + * - randomizeOrder() (qsort + random keys to match Elegant’s consumption) + * + * Usage Notes + * ----------- + * - Single-threaded, static state; synchronize if used from multiple threads. + * - Calls with negative seeds reinitialize the stream; non-negative seeds consume. + * - MPI seed diversification is intentionally omitted here. + * - Special behavior for seed 987654321 (permute inhibition) is preserved. + * - Goal: reproduce Elegant sequences (bitwise where possible). + * + * References + * ---------- + * - M. Borland, “elegant: A Flexible SDDS-Compliant Code for Accelerator Simulation,” + * APS LS-287 (2000). */ -// -// This file mirrors algorithms and seeding conventions from: -// - Elegant: M. Borland, “elegant: A Flexible SDDS-Compliant Code for Accelerator Simulation,” -// Advanced Photon Source LS-287, September 2000. -// - SDDS: see the SDDS and Elegant source distributions and their LICENSE files. -// -// Please retain upstream copyright and license notices where code/logic -// has been adapted, and also cite LS-287 when publishing results produced with this RNG. -// -// ----------------------------------------------------------------------------- -// Purpose -// ------- -// Re-implements the random-number utilities used by Elegant/SDDS so that -// C-kernels compiled via xobjects can reproduce *bitwise-identical* -// sequences to Elegant. This file exposes: -// -// - A LAPACK-compatible 48-bit LCG core (DLARAN) working on 4×12-bit chunks -// - Six RNG streams random_1..random_6 with Elegant-compatible seeding rules -// - Seed-bit permutation with Elegant’s “inhibit” switch and special seed -// behavior (987654321) -// - randomizeOrder() that shuffles buffers using the same qsort+random-key -// scheme used by Elegant to match sequence consumption -// -// Notes -// ----- -// * Based on the DLARAN routine from LAPACK and the RNG/seeding conventions from -// SDDS/Elegant (e.g., drand.c and friends in the Elegant/SDDS sources). -// * This is a single-threaded, static-state implementation. -// If you call these from multiple threads, guard access yourself. -// * MPI-related seed diversification is (for now) intentionally omitted here. -// * “Negative seed” calls *reinitialize* the stream; “non-negative” consume. -// -// Attribution -// ----------- -// - DLARAN: LAPACK auxiliary RNG (48-bit LCG with multiplier 33952834046453). -// - Seeding conventions and the 987654321 “inhibit permutation” behavior mirror -// SDDS/Elegant (random_1..random_6, randomizeOrder). -// ----------------------------------------------------------------------------- - #ifndef ELEGANT_RNG_H #define ELEGANT_RNG_H From 5b59ddbe0665828977d93f398a0047d66ec771ba Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 12:28:31 +0200 Subject: [PATCH 15/37] Cleanup contributors_and_credits.txt --- contributors_and_credits.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/contributors_and_credits.txt b/contributors_and_credits.txt index e8253569..3af431d3 100644 --- a/contributors_and_credits.txt +++ b/contributors_and_credits.txt @@ -20,6 +20,7 @@ External contributors: The Particle In Cell implementation is largely based on the PyPIC package (https://github.com/pycomplete/pypic) -The Touschek scattering routine in this package contains code adapted from Elegant and SDDS. +The Touschek scattering routine contains portions adapted from Elegant and the SDDS Toolkit. © 2002 The University of Chicago; © 2002 The Regents of the University of California. Distributed subject to their Software License Agreements. See third_party/elegant/LICENSE and third_party/SDDS/LICENSE. +The RNG core used in the Touschek scattering routine incorporates LAPACK’s DLARAN (modified BSD license). See third_party/lapack/LICENSE. \ No newline at end of file From a411d27de915925d4c45864c5eef33ea7532f239 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 12:28:48 +0200 Subject: [PATCH 16/37] Add NOTICE --- MANIFEST.in | 8 ++++++-- NOTICE | 7 +++++++ 2 files changed, 13 insertions(+), 2 deletions(-) create mode 100644 NOTICE diff --git a/MANIFEST.in b/MANIFEST.in index 829d9280..e1fa24b4 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -3,8 +3,12 @@ include pyproject.toml # Include the README include *.md -# Include the license file +# Include the license files include LICENSE.txt +recursive-include xfields/third_party *LICENSE* + +# Include the notice +include NOTICE recursive-include xfields *.h -recursive-include xfields *.clh +recursive-include xfields *.clh \ No newline at end of file diff --git a/NOTICE b/NOTICE new file mode 100644 index 00000000..8c0c5460 --- /dev/null +++ b/NOTICE @@ -0,0 +1,7 @@ +This package incorporates portions adapted from Argonne National Laboratory’s “Elegant” +and from the SDDS Toolkit. © 2002 The University of Chicago; © 2002 The Regents of the +University of California. Distributed under their respective Software License Agreements. +See: xfields/third_party/elegant/LICENSE and xfields/third_party/SDDS/LICENSE. + +This package also incorporates LAPACK’s DLARAN algorithm (modified BSD license). +See: xfields/third_party/lapack/LICENSE. \ No newline at end of file From 11a3cbf3354340573a3dc0d291ea7487764ac5ec Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 12:29:01 +0200 Subject: [PATCH 17/37] Add LAPACK license in third_party --- third_party/LAPACK/LICENSE | 48 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 third_party/LAPACK/LICENSE diff --git a/third_party/LAPACK/LICENSE b/third_party/LAPACK/LICENSE new file mode 100644 index 00000000..c9ae626c --- /dev/null +++ b/third_party/LAPACK/LICENSE @@ -0,0 +1,48 @@ +Copyright (c) 1992-2013 The University of Tennessee and The University + of Tennessee Research Foundation. All rights + reserved. +Copyright (c) 2000-2013 The University of California Berkeley. All + rights reserved. +Copyright (c) 2006-2013 The University of Colorado Denver. All rights + reserved. + +$COPYRIGHT$ + +Additional copyrights may follow + +$HEADER$ + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +- Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +- Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer listed + in this license in the documentation and/or other materials + provided with the distribution. + +- Neither the name of the copyright holders nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +The copyright holders provide no reassurances that the source code +provided does not infringe any patent, copyright, or any other +intellectual property rights of third parties. The copyright holders +disclaim any liability to any recipient for claims brought against +recipient by any third party for infringement of that parties +intellectual property rights. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. \ No newline at end of file From 91b7eccab54700e3765bc3a653d15d8f1e3e0aa3 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sun, 7 Sep 2025 12:29:16 +0200 Subject: [PATCH 18/37] Cleanup touschek readme. --- xfields/touschek/README.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/xfields/touschek/README.md b/xfields/touschek/README.md index 43e3e81c..aa5570b2 100644 --- a/xfields/touschek/README.md +++ b/xfields/touschek/README.md @@ -1,6 +1,6 @@ # Monte Carlo simulation of Touschek scattering -The Monte-Carlo Touschek scattering routine in this package is based on the implementation by A. Xiao and M. Borland developed for [ELEGANT](https://github.com/rtsoliday/elegant). +The Monte Carlo Touschek scattering routine in this package is based on the implementation by A. Xiao and M. Borland developed for [ELEGANT](https://github.com/rtsoliday/elegant). If you publish results obtained with this routine, please cite: @@ -13,6 +13,9 @@ If you publish results obtained with this routine, please cite: For reproducibility with ELEGANT, this routine uses the same RNG conventions; see `xfields/xfields/headers/elegant_rng.h`. +The uniform RNG core is LAPACK’s DLARAN (48-bit LCG, modified BSD license) as used by Elegant. +See `xfields/third_party/LAPACK/LICENSE`. + ## Third-party notice and licenses This routine contains code adapted from **ELEGANT** and **SDDS**. From 74b85a2af4ccb00d0d406818413108c86b417cc3 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Wed, 10 Sep 2025 10:20:17 +0200 Subject: [PATCH 19/37] Take momentum aperture based on s instead of element_name. --- xfields/touschek/manager.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index 674faaeb..be0a6eeb 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -90,8 +90,9 @@ def _compute_piwinski_scattering_rate(self, element): dpy = self.twiss['dpy', element] dyt = alfy * dy + bety * dpy # dyt: dy tilde - deltaN = self.manager.momentum_aperture.at[element, "deltaN"] - deltaP = self.manager.momentum_aperture.at[element, "deltaP"] + s = self.twiss.rows[element].s[0] + deltaN = np.interp(s, self.manager.momentum_aperture.s, self.manager.momentum_aperture.deltan) + deltaP = np.interp(s, self.manager.momentum_aperture.s, self.manager.momentum_aperture.deltap) sigmab_x = np.sqrt(gemitt_x * betx) # Horizontal betatron beam size sigma_x = np.sqrt(gemitt_x * betx + dx**2 * sigma_delta**2) # Horizontal beam size @@ -197,15 +198,15 @@ def __init__(self, line, momentum_aperture, nemitt_x=None, nemitt_y=None, raise ValueError("`n_simulated` is required.") # Momentum aperture validation - required_cols = {"s", "name", "deltaN", "deltaP"} + required_cols = {"s", "deltan", "deltap"} if not hasattr(momentum_aperture, "columns") or not hasattr(momentum_aperture, "__getitem__"): raise TypeError("`momentum_aperture` must be a DataFrame-like object with columns " - "'s', 'name', 'deltaN', 'deltaP'.") + "'s', 'deltan', 'deltap'.") missing = required_cols - set(momentum_aperture.columns) if missing: raise ValueError(f"`momentum_aperture` missing columns: {sorted(missing)}") - for col in ("s", "deltaN", "deltaP"): + for col in ("s", "deltan", "deltap"): try: vals = momentum_aperture[col].astype(float) except Exception: @@ -221,9 +222,8 @@ def __init__(self, line, momentum_aperture, nemitt_x=None, nemitt_y=None, self.particle_ref = line.particle_ref momentum_aperture = momentum_aperture.copy() - momentum_aperture['deltaN'] *= momentum_aperture_scale - momentum_aperture['deltaP'] *= momentum_aperture_scale - self.momentum_aperture= momentum_aperture.set_index("name") + momentum_aperture['deltan'] *= momentum_aperture_scale + momentum_aperture['deltap'] *= momentum_aperture_scale self.sigma_z = sigma_z self.sigma_delta = sigma_delta @@ -289,8 +289,8 @@ def _config(nn): alfy = twiss["alfy", nn]; bety = twiss["bety", nn] dx = twiss["dx", nn]; dpx = twiss["dpx", nn] dy = twiss["dy", nn]; dpy = twiss["dpy", nn] - dN = self.momentum_aperture.at[nn, "deltaN"] - dP = self.momentum_aperture.at[nn, "deltaP"] + dN = np.interp(s, self.momentum_aperture.s, self.momentum_aperture.deltan) + dP = np.interp(s, self.momentum_aperture.s, self.momentum_aperture.deltap) piwinski_rate = self.touschek._compute_piwinski_scattering_rate(nn) From c345044b7a9ef30bcc4fb8d9bd8816f8165f2b69 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Wed, 10 Sep 2025 11:27:16 +0200 Subject: [PATCH 20/37] Add back momentum_aperture attribute. --- xfields/touschek/manager.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index be0a6eeb..77ed4b04 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -90,7 +90,10 @@ def _compute_piwinski_scattering_rate(self, element): dpy = self.twiss['dpy', element] dyt = alfy * dy + bety * dpy # dyt: dy tilde - s = self.twiss.rows[element].s[0] + s = self.manager.line.get_s_position(element) + # The following would be better (or line table): + # s = self.twiss.rows[element].s[0] + # However, sequencename$start and sequencename$end do not have s in the tables! deltaN = np.interp(s, self.manager.momentum_aperture.s, self.manager.momentum_aperture.deltan) deltaP = np.interp(s, self.manager.momentum_aperture.s, self.manager.momentum_aperture.deltap) @@ -224,6 +227,7 @@ def __init__(self, line, momentum_aperture, nemitt_x=None, nemitt_y=None, momentum_aperture = momentum_aperture.copy() momentum_aperture['deltan'] *= momentum_aperture_scale momentum_aperture['deltap'] *= momentum_aperture_scale + self.momentum_aperture = momentum_aperture self.sigma_z = sigma_z self.sigma_delta = sigma_delta From 60ad9bc01d5d56b2a95f27c5b4fa1d839f5b2f95 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Wed, 10 Sep 2025 11:31:30 +0200 Subject: [PATCH 21/37] Line table not suitable when lines come from MAD-X conversion. When a line is generated through MAD-X sequence conversion, the line end up having markers sequencename$start and sequencename$end. In the line table these are present, but they do not have an associated s position and this causes an error! --- xfields/touschek/manager.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index 77ed4b04..92a81c83 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -159,7 +159,7 @@ def _compute_integrated_piwinski_rates(self): ii_current = 0 for ii, nn in enumerate(tab.name): - s = tab.rows[nn].s[0] + s = self.manager.line.get_s_position(nn) ds = s - s_before if ds > 0.0: rate = self._compute_piwinski_scattering_rate(nn) @@ -288,7 +288,7 @@ def initialise_touschek(self, element=None): # Helper to config all fields to a single TouschekScattering def _config(nn): - s = tab.rows[nn].s[0] + s = self.line.get_s_position(nn) alfx = twiss["alfx", nn]; betx = twiss["betx", nn] alfy = twiss["alfy", nn]; bety = twiss["bety", nn] dx = twiss["dx", nn]; dpx = twiss["dpx", nn] From a8483e4e5098bbfa10dc85bfd747678df8d4eddf Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Mon, 15 Sep 2025 18:22:40 +0200 Subject: [PATCH 22/37] Speed up Piwinski integral. --- xfields/touschek/manager.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index 92a81c83..7d839352 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -30,34 +30,34 @@ def _compute_piwinski_integral(self, tm, B1, B2): arXiv:physics/9903034, 1999. URL: https://arxiv.org/abs/physics/9903034 """ - km = np.arctan(np.sqrt(tm)) + from math import atan, tan, sqrt, exp, log, pi - def int_piwinski(k, km, B1, B2): + km = atan(sqrt(tm)) + tm + + def int_piwinski(k): t = np.tan(k) ** 2 - tm = np.tan(km) ** 2 fact = ( - (2*t + 1)**2 * (t/tm / (1+t) - 1) / t + t - np.sqrt(t*tm * (1 + t)) - - (2 + 1 / (2*t)) * np.log(t/tm / (1+t)) + (2*t + 1)**2 * (t/tm / (1+t) - 1) / t + t - sqrt(t*tm * (1 + t)) + - (2 + 1 / (2*t)) * log(t/tm / (1+t)) ) if B2 * t < 500: - intp = fact * np.exp(-B1*t) * i0(B2*t) * np.sqrt(1+t) + intp = fact * exp(-B1*t) * i0(B2*t) * sqrt(1+t) else: intp = ( fact - * np.exp(B2*t - B1*t) - / np.sqrt(2*np.pi * B2*t) - * np.sqrt(1+t) + * exp(B2*t - B1*t) + / sqrt(2*pi * B2*t) + * sqrt(1+t) ) return intp - args = (km, B1, B2) val, _ = quad( int_piwinski, km, - np.pi / 2, - args=args, - epsabs=1e-16, - epsrel=1e-12 + pi / 2, + epsabs=1e-12, + epsrel=1e-8 ) return val From 2da379bc2a02b3813a000eb09bd546b341d141ae Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Mon, 15 Sep 2025 18:23:27 +0200 Subject: [PATCH 23/37] Momentum aperture selection based on s instead of element_name. --- xfields/touschek/manager.py | 180 ++++++++++++++++++++++++------------ 1 file changed, 122 insertions(+), 58 deletions(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index 7d839352..c614c3d2 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -7,6 +7,7 @@ import xfields as xf import numpy as np +import pandas as pd from scipy.integrate import quad from scipy.special import i0 from scipy.constants import physical_constants @@ -76,26 +77,28 @@ def _compute_piwinski_scattering_rate(self, element): bunch_population = self.manager.bunch_population gemitt_x = self.manager.gemitt_x gemitt_y = self.manager.gemitt_y - alfx = self.twiss['alfx', element] - betx = self.twiss['betx', element] - alfy = self.twiss['alfy', element] - bety = self.twiss['bety', element] + twiss = self.twiss + alfx = twiss['alfx', element] + betx = twiss['betx', element] + alfy = twiss['alfy', element] + bety = twiss['bety', element] sigma_z = self.manager.sigma_z sigma_delta = self.manager.sigma_delta - delta = self.twiss['delta', element] - dx = self.twiss['dx', element] - dpx = self.twiss['dpx', element] + delta = twiss['delta', element] + dx = twiss['dx', element] + dpx = twiss['dpx', element] dxt = alfx * dx + betx * dpx # dxt: dx tilde - dy = self.twiss['dy', element] - dpy = self.twiss['dpy', element] + dy = twiss['dy', element] + dpy = twiss['dpy', element] dyt = alfy * dy + bety * dpy # dyt: dy tilde - s = self.manager.line.get_s_position(element) - # The following would be better (or line table): - # s = self.twiss.rows[element].s[0] - # However, sequencename$start and sequencename$end do not have s in the tables! - deltaN = np.interp(s, self.manager.momentum_aperture.s, self.manager.momentum_aperture.deltan) - deltaP = np.interp(s, self.manager.momentum_aperture.s, self.manager.momentum_aperture.deltap) + try: + s = self.twiss.rows[element].s[0] + except: + s = self.manager.line.get_s_position(element) + + deltaN = self.manager.momentum_aperture.at[s, "deltan"] + deltaP = self.manager.momentum_aperture.at[s, "deltap"] sigmab_x = np.sqrt(gemitt_x * betx) # Horizontal betatron beam size sigma_x = np.sqrt(gemitt_x * betx + dx**2 * sigma_delta**2) # Horizontal beam size @@ -133,7 +136,7 @@ def _compute_piwinski_scattering_rate(self, element): return rate - def _compute_integrated_piwinski_rates(self): + def _compute_integrated_piwinski_rates(self, element): """ Integrate the Piwinski Touschek scattering rate along the line using the trapezoidal rule, between successive TouschekScattering elements. @@ -143,6 +146,22 @@ def _compute_integrated_piwinski_rates(self): rate is later used to assign the correct weights to Touschek-scattered particles at the corresponding element. """ + def _get_s(name): + try: + return tab.rows[name].s[0] + except (KeyError, AttributeError, IndexError, TypeError): + return self.manager.line.get_s_position(name) + + def _step(name, s_before, rate_before, integrated): + s = _get_s(name) + ds = s - s_before + if ds > 0.0: + rate = self._compute_piwinski_scattering_rate(name) + integrated += 0.5 * (rate_before + rate) * ds + return s, rate, integrated + else: + return s_before, rate_before, integrated + line = self.manager.line tab = line.get_table() T_rev0 = float(self.twiss.T_rev0) @@ -151,34 +170,54 @@ def _compute_integrated_piwinski_rates(self): ii_t = [ii for ii, nn in enumerate(tab.name[:-1]) if isinstance(line[nn], xf.TouschekScattering)] integrated = 0.0 - s0 = 0.0 - r0 = self._compute_piwinski_scattering_rate(tab.name[0]) + + if element is None: + ii_current = 0 + s0 = 0.0 + r0 = self._compute_piwinski_scattering_rate(tab.name[0]) + else: + import re + ii_current = int(re.search(r'\d+', element).group()) + tscatter_before = tab.name[ii_t[ii_current - 1]] if ii_current != 0 else tab.name[0] + s0 = _get_s(tscatter_before) + r0 = self._compute_piwinski_scattering_rate(tscatter_before) s_before = s0 rate_before = r0 - ii_current = 0 - for ii, nn in enumerate(tab.name): - s = self.manager.line.get_s_position(nn) - ds = s - s_before - if ds > 0.0: - rate = self._compute_piwinski_scattering_rate(nn) - integrated += 0.5 * (rate_before + rate) * ds - s_before = s - rate_before = rate - - if ii_current < len(ii_t) and ii == ii_t[ii_current]: - # divide by c and by T_rev0 --> per-bunch rate - integrated_piwinski_rate = integrated / C_LIGHT_VACUUM / T_rev0 - elem = line[nn] # xf.TouschekScattering - elem._configure(_integrated_piwinski_rate=integrated_piwinski_rate) - integrated = 0.0 - ii_current += 1 - if ii_current >= len(ii_t): + if element is None: + # Configure all the TouschekScattering elements + for ii, nn in enumerate(tab.name): + s_before, rate_before, integrated = _step(nn, s_before, rate_before, integrated) + + if ii == ii_t[ii_current]: + # divide by c and by T_rev0 --> per-bunch rate + integrated_piwinski_rate = integrated / C_LIGHT_VACUUM / T_rev0 + elem = line[nn] # xf.TouschekScattering + # print(f'Integrated Piwinski rate at {nn}: {integrated_piwinski_rate*1e-3} [kHz]') + elem._configure(_integrated_piwinski_rate=integrated_piwinski_rate) + integrated = 0.0 + ii_current += 1 + if ii_current == len(ii_t): + break + else: + # Configure only the TouschekScattering element named `element` + subtab = tab.rows[tscatter_before:element] + for nn in subtab.name: + s_before, rate_before, integrated = _step(nn, s_before, rate_before, integrated) + + if nn == element: + # divide by c and by T_rev0 --> per-bunch rate + integrated_piwinski_rate = integrated / C_LIGHT_VACUUM / T_rev0 + elem = line[nn] # xf.TouschekScattering + # print(f'Integrated Piwinski rate at {nn}: {integrated_piwinski_rate*1e-3} [kHz]') + elem._configure(_integrated_piwinski_rate=integrated_piwinski_rate) break + class TouschekManager: - def __init__(self, line, momentum_aperture, nemitt_x=None, nemitt_y=None, + def __init__(self, line=None, twiss=None, momentum_aperture=None, + nemitt_x=None, nemitt_y=None, sigma_z=None, sigma_delta=None, bunch_population=None, n_simulated=None, gemitt_x=None, gemitt_y=None, momentum_aperture_scale=0.85, ignored_portion=0.01, @@ -223,10 +262,34 @@ def __init__(self, line, momentum_aperture, nemitt_x=None, nemitt_y=None, self.line = line self.particle_ref = line.particle_ref + self.twiss = twiss - momentum_aperture = momentum_aperture.copy() - momentum_aperture['deltan'] *= momentum_aperture_scale - momentum_aperture['deltap'] *= momentum_aperture_scale + # Check that the line contains TouschekScatterings + tab = line.get_table() + try: + has = "TouschekScattering" in set(np.unique(tab.element_type)) + except Exception: + has = "TouschekScattering" in set(getattr(tab, "element_type", [])) + if not has: + raise ValueError("The line does not contain any TouschekScattering. " + "Please add them before initializing the TouschekManager.") + + # Momentum aperture + ap = momentum_aperture.copy() + ap['deltan'] *= momentum_aperture_scale + ap['deltap'] *= momentum_aperture_scale + + deltan = np.interp(tab.s, ap.s, ap.deltan) + deltap = np.interp(tab.s, ap.s, ap.deltap) + + momentum_aperture = pd.DataFrame({ + "s": tab.s, + "deltan": deltan, + "deltap": deltap, + }, index=tab.s) + + # Remove duplicates + momentum_aperture = momentum_aperture[~momentum_aperture.index.duplicated(keep="first")] self.momentum_aperture = momentum_aperture self.sigma_z = sigma_z @@ -264,37 +327,36 @@ def __init__(self, line, momentum_aperture, nemitt_x=None, nemitt_y=None, self.touschek = TouschekCalculator(self) - # Check that the line contains TouschekScatterings - tab = self.line.get_table() - try: - has = "TouschekScattering" in set(np.unique(tab.element_type)) - except Exception: - has = "TouschekScattering" in set(getattr(tab, "element_type", [])) - if not has: - raise ValueError("The line does not contain any TouschekScattering. " - "Please add them before initializing the TouschekManager.") - - def initialise_touschek(self, element=None): line = self.line tab = line.get_table() - twiss_method = self.kwargs.get("method", "6d") - twiss = self.line.twiss(method=twiss_method) + if self.twiss is None: + twiss_method = self.kwargs.get("method", "6d") + self.twiss = self.line.twiss(method=twiss_method) + # Pass the twiss to the TouschekCalculator - self.touschek.twiss = twiss + self.touschek.twiss = self.twiss - self.touschek._compute_integrated_piwinski_rates() + # import time + # t0 = time.time() + self.touschek._compute_integrated_piwinski_rates(element) + # print(f"Computed integrated piwinski rates in {time.time() - t0:.2f} s.") # Helper to config all fields to a single TouschekScattering def _config(nn): - s = self.line.get_s_position(nn) + try: + s = tab.rows[nn].s[0] + except Exception: + s = self.line.get_s_position(nn) + + twiss = self.twiss alfx = twiss["alfx", nn]; betx = twiss["betx", nn] alfy = twiss["alfy", nn]; bety = twiss["bety", nn] dx = twiss["dx", nn]; dpx = twiss["dpx", nn] dy = twiss["dy", nn]; dpy = twiss["dpy", nn] - dN = np.interp(s, self.momentum_aperture.s, self.momentum_aperture.deltan) - dP = np.interp(s, self.momentum_aperture.s, self.momentum_aperture.deltap) + dN = self.momentum_aperture.at[s, "deltan"] + dP = self.momentum_aperture.at[s, "deltap"] piwinski_rate = self.touschek._compute_piwinski_scattering_rate(nn) @@ -323,6 +385,7 @@ def _config(nn): if element is None: for nn in tab.name[:-1]: # Avoid the last tab.name which is _end_point if isinstance(line[nn], xf.TouschekScattering): + print(f'Initialising TouschekScattering for {nn}') _config(nn) else: if not isinstance(element, str): @@ -335,4 +398,5 @@ def _config(nn): raise TypeError( f"`line['{element}']` is not a TouschekScattering (got {type(line[element]).__name__})." ) + print(f'Initialising TouschekScattering for {element}') _config(element) \ No newline at end of file From ca69ce373ef29394b05e88f85a71ae81ce9b849c Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Mon, 15 Sep 2025 18:23:55 +0200 Subject: [PATCH 24/37] Inherit some warning prints from ELEGANT for safety. --- xfields/beam_elements/touschek_src/touschek.h | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h index 4e4c2bcc..c81bdca5 100644 --- a/xfields/beam_elements/touschek_src/touschek.h +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -524,6 +524,27 @@ void TouschekScatter(TouschekScatteringData el, } } + // printf("Total number of random numbers used: %ld\n", total_event * 11); + // fflush(stdout); + if (total_event * 11 > (long)2e9) { + printf("The total number of random numbers used exceeded 2e9. Use smaller n_simulated or smaller delta."); + fflush(stdout); + } + + // if (simuCount == 0) { + // printf("It appears that the Touschek lifetime is extremely wrong and there is no need to perform Touschek simulation.\n If you think this is a wrong statement, send input to developers for evaluation.\n"); + // } + + if (total_event / simuCount > 20) { + if (nx < 5 || ny < 5) { + printf("The Touschek scattering rate is low. Please use >=5 sigma beam for better simulation."); + fflush(stdout); + } else { + printf("The Touschek scattering rate is very low. Please ignore the rate from Monte Carlo simulation. Use Piwinski rate only."); + fflush(stdout); + } + } + printf("%ld of %ld particles selected for tracking\n", iTotal, simuCount); fflush(stdout); From 9d6e53ee7f076a8cc24f9c8ac43322b1a2479d44 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sat, 20 Sep 2025 16:37:00 +0200 Subject: [PATCH 25/37] Remove typo. --- xfields/touschek/manager.py | 1 - 1 file changed, 1 deletion(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index c614c3d2..420c6643 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -34,7 +34,6 @@ def _compute_piwinski_integral(self, tm, B1, B2): from math import atan, tan, sqrt, exp, log, pi km = atan(sqrt(tm)) - tm def int_piwinski(k): t = np.tan(k) ** 2 From 35843b1578196f8a539385ea5ad25ef0389a4596 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sat, 20 Sep 2025 16:42:10 +0200 Subject: [PATCH 26/37] Change momentum aperture handling. --- xfields/touschek/manager.py | 32 ++++++++++++-------------------- 1 file changed, 12 insertions(+), 20 deletions(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index 420c6643..028de4d8 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -74,6 +74,7 @@ def _compute_piwinski_scattering_rate(self, element): """ p0c = self.manager.particle_ref.p0c[0] bunch_population = self.manager.bunch_population + momentum_aperture = self.manager.momentum_aperture gemitt_x = self.manager.gemitt_x gemitt_y = self.manager.gemitt_y twiss = self.twiss @@ -95,9 +96,9 @@ def _compute_piwinski_scattering_rate(self, element): s = self.twiss.rows[element].s[0] except: s = self.manager.line.get_s_position(element) - - deltaN = self.manager.momentum_aperture.at[s, "deltan"] - deltaP = self.manager.momentum_aperture.at[s, "deltap"] + + deltaN = np.interp(s, momentum_aperture.s, momentum_aperture.deltan) + deltaP = np.interp(s, momentum_aperture.s, momentum_aperture.deltap) sigmab_x = np.sqrt(gemitt_x * betx) # Horizontal betatron beam size sigma_x = np.sqrt(gemitt_x * betx + dx**2 * sigma_delta**2) # Horizontal beam size @@ -277,18 +278,7 @@ def __init__(self, line=None, twiss=None, momentum_aperture=None, ap = momentum_aperture.copy() ap['deltan'] *= momentum_aperture_scale ap['deltap'] *= momentum_aperture_scale - - deltan = np.interp(tab.s, ap.s, ap.deltan) - deltap = np.interp(tab.s, ap.s, ap.deltap) - - momentum_aperture = pd.DataFrame({ - "s": tab.s, - "deltan": deltan, - "deltap": deltap, - }, index=tab.s) - - # Remove duplicates - momentum_aperture = momentum_aperture[~momentum_aperture.index.duplicated(keep="first")] + momentum_aperture = ap self.momentum_aperture = momentum_aperture self.sigma_z = sigma_z @@ -330,6 +320,8 @@ def initialise_touschek(self, element=None): line = self.line tab = line.get_table() + momentum_aperture = self.momentum_aperture + if self.twiss is None: twiss_method = self.kwargs.get("method", "6d") self.twiss = self.line.twiss(method=twiss_method) @@ -337,10 +329,10 @@ def initialise_touschek(self, element=None): # Pass the twiss to the TouschekCalculator self.touschek.twiss = self.twiss - # import time - # t0 = time.time() + import time + t0 = time.time() self.touschek._compute_integrated_piwinski_rates(element) - # print(f"Computed integrated piwinski rates in {time.time() - t0:.2f} s.") + print(f"Computed integrated piwinski rates in {time.time() - t0:.2f} s.") # Helper to config all fields to a single TouschekScattering def _config(nn): @@ -354,8 +346,8 @@ def _config(nn): alfy = twiss["alfy", nn]; bety = twiss["bety", nn] dx = twiss["dx", nn]; dpx = twiss["dpx", nn] dy = twiss["dy", nn]; dpy = twiss["dpy", nn] - dN = self.momentum_aperture.at[s, "deltan"] - dP = self.momentum_aperture.at[s, "deltap"] + dN = np.interp(s, momentum_aperture.s, momentum_aperture.deltan) + dP = np.interp(s, momentum_aperture.s, momentum_aperture.deltap) piwinski_rate = self.touschek._compute_piwinski_scattering_rate(nn) From cdb7b6a4de8cc42c944bc3abe6d9cbd09efc7f5e Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Sat, 20 Sep 2025 16:43:28 +0200 Subject: [PATCH 27/37] Introduce theta_min and theta_max as user_defined parameters. --- xfields/beam_elements/touschek.py | 12 ++++++++++++ xfields/beam_elements/touschek_src/touschek.h | 11 ++++++++++- xfields/touschek/manager.py | 7 ++++++- 3 files changed, 28 insertions(+), 2 deletions(-) diff --git a/xfields/beam_elements/touschek.py b/xfields/beam_elements/touschek.py index 63ca3b47..00a3b644 100644 --- a/xfields/beam_elements/touschek.py +++ b/xfields/beam_elements/touschek.py @@ -32,6 +32,8 @@ class TouschekScattering(xt.BeamElement): '_nx': xo.Float64, '_ny': xo.Float64, '_nz': xo.Float64, + '_theta_min': xo.Float64, + '_theta_max': xo.Float64, '_ignored_portion': xo.Float64, '_integrated_piwinski_rate': xo.Float64, '_seed': xo.Int64, @@ -57,6 +59,7 @@ class TouschekScattering(xt.BeamElement): xo.Arg(xo.Float64, name='py_out', pointer=True), xo.Arg(xo.Float64, name='zeta_out', pointer=True), xo.Arg(xo.Float64, name='delta_out', pointer=True), + xo.Arg(xo.Float64, name='theta_out', pointer=True), xo.Arg(xo.Float64, name='weight_out', pointer=True), xo.Arg(xo.Float64, name='totalMCRate_out', pointer=True), xo.Arg(xo.Int64, name='n_selected_out', pointer=True), @@ -73,6 +76,7 @@ def __init__(self, s=0.0, gemitt_x=0.0, gemitt_y=0.0, sigma_z=0.0, sigma_delta=0.0, n_simulated=0, nx=0.0, ny=0.0, nz=0.0, + theta_min=0.0, theta_max=0.0, piwinski_rate=0.0, ignored_portion=0.0, integrated_piwinski_rate=0.0, @@ -111,6 +115,8 @@ def __init__(self, s=0.0, self._nx = nx self._ny = ny self._nz = nz + self._theta_min = theta_min + self._theta_max = theta_max self._ignored_portion = ignored_portion self._integrated_piwinski_rate = integrated_piwinski_rate self.piwinski_rate = piwinski_rate @@ -126,6 +132,7 @@ def _configure(self, **kwargs): "_deltaN", "_deltaP", "_sigma_z", "_sigma_delta", "_n_simulated", "_nx", "_ny", "_nz", + "_theta_min", "_theta_max", "_ignored_portion", "piwinski_rate", "_integrated_piwinski_rate", "_seed", "_inhibit_permute" @@ -154,6 +161,7 @@ def scatter(self): py_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) zeta_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) delta_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) + theta_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) weight_out = context.zeros(shape=(self._n_simulated,), dtype=np.float64) totalMCRate_out = context.zeros(shape=(1,), dtype=np.float64) n_selected_out = context.zeros(shape=(1,), dtype=np.int64) @@ -162,6 +170,7 @@ def scatter(self): x_out=x_out, px_out=px_out, y_out=y_out, py_out=py_out, zeta_out=zeta_out, delta_out=delta_out, + theta_out=theta_out, weight_out=weight_out, totalMCRate_out=totalMCRate_out, n_selected_out=n_selected_out) @@ -180,6 +189,9 @@ def scatter(self): weight=weight_out[:n], s=getattr(self, '_s', 0.0)) + part_ids = part.filter(part.state == 1).particle_id + self.theta_log = dict(zip(part_ids.astype(int), theta_out[:n].astype(float))) + self.total_mc_rate = totalMCRate_out[0] self.ignored_rate = self._ignored_portion * self.total_mc_rate diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h index c81bdca5..e00fc04d 100644 --- a/xfields/beam_elements/touschek_src/touschek.h +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -297,6 +297,7 @@ void TouschekScatter(TouschekScatteringData el, double* py_out, double* zeta_out, double* delta_out, + double* theta_out, double* weight_out, double* totalMCRate_out, int64_t* n_selected_out){ @@ -321,6 +322,8 @@ void TouschekScatter(TouschekScatteringData el, const double nx = TouschekScatteringData_get__nx(el); const double ny = TouschekScatteringData_get__ny(el); const double nz = TouschekScatteringData_get__nz(el); + const double theta_min = TouschekScatteringData_get__theta_min(el); + const double theta_max = TouschekScatteringData_get__theta_max(el); const double ignoredPortion = TouschekScatteringData_get__ignored_portion(el); const double integrated_piwinski_rate = TouschekScatteringData_get__integrated_piwinski_rate(el); @@ -361,6 +364,7 @@ void TouschekScatter(TouschekScatteringData el, double *pytemp = (double*)malloc(sizeof(double) * n_simulated); double *zetatemp = (double*)malloc(sizeof(double) * n_simulated); double *deltatemp = (double*)malloc(sizeof(double) * n_simulated); + double *thetatemp = (double*)malloc(sizeof(double) * n_simulated); static int seeded_once = 0; if (!seeded_once){ @@ -420,7 +424,7 @@ void TouschekScatter(TouschekScatteringData el, bunch2cm(p1, p2, qa, beta, &gamma); - theta = (ran1[9] * 0.9999 + 0.00005) * PI; + theta = theta_min + (theta_max - theta_min) * ran1[9]; phi = ran1[10] * PI; temp = dens1 * dens2 * sin(theta); @@ -457,6 +461,7 @@ void TouschekScatter(TouschekScatteringData el, pytemp[i] = p1[4]; zetatemp[i] = p1[2]; deltatemp[i] = p1[5]; + thetatemp[i] = theta; weight[i] = temp; i++; } @@ -476,6 +481,7 @@ void TouschekScatter(TouschekScatteringData el, pytemp[i] = p2[4]; zetatemp[i] = p2[2]; deltatemp[i] = p2[5]; + thetatemp[i] = theta; weight[i] = temp; i++; } @@ -501,6 +507,7 @@ void TouschekScatter(TouschekScatteringData el, py_out[k] = pytemp[k]; zeta_out[k] = zetatemp[k]; delta_out[k] = deltatemp[k]; + theta_out[k] = thetatemp[k]; weight_out[k] = weight[k]; } } else { @@ -520,6 +527,7 @@ void TouschekScatter(TouschekScatteringData el, py_out[k] = pytemp[src]; zeta_out[k] = zetatemp[src]; delta_out[k] = deltatemp[src]; + theta_out[k] = thetatemp[src]; weight_out[k] = weight[src]; } } @@ -557,6 +565,7 @@ void TouschekScatter(TouschekScatteringData el, *totalMCRate_out = totalMCRate; free(index); + free(thetatemp); free(weight); free(xtemp); free(pxtemp); diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index 028de4d8..d5dc7c34 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -221,7 +221,8 @@ def __init__(self, line=None, twiss=None, momentum_aperture=None, sigma_z=None, sigma_delta=None, bunch_population=None, n_simulated=None, gemitt_x=None, gemitt_y=None, momentum_aperture_scale=0.85, ignored_portion=0.01, - seed=1997, nx=3, ny=3, nz=3, **kwargs): + seed=1997, nx=3, ny=3, nz=3, + theta_min=0.00005*np.pi, theta_max=0.99995*np.pi, **kwargs): # Input validation if line is None: @@ -292,6 +293,9 @@ def __init__(self, line=None, twiss=None, momentum_aperture=None, self.ny = ny self.nz = nz + self.theta_min = theta_min + self.theta_max = theta_max + # Emittance validation nemitt_given = nemitt_x is not None and nemitt_y is not None gemitt_given = gemitt_x is not None and gemitt_y is not None @@ -368,6 +372,7 @@ def _config(nn): _sigma_delta=self.sigma_delta, _n_simulated=self.n_simulated, _nx=self.nx, _ny=self.ny, _nz=self.nz, + _theta_min=self.theta_min, _theta_max=self.theta_max, _ignored_portion=self.ignored_portion, piwinski_rate=piwinski_rate, _seed=self.seed, _inhibit_permute=0 From 3a01e6803625bfdbfa813175f637ef5b721e17b6 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Mon, 22 Sep 2025 13:49:55 +0200 Subject: [PATCH 28/37] theta_min and theta_max must not be set by the user. --- xfields/touschek/manager.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index d5dc7c34..c3721e5c 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -221,8 +221,7 @@ def __init__(self, line=None, twiss=None, momentum_aperture=None, sigma_z=None, sigma_delta=None, bunch_population=None, n_simulated=None, gemitt_x=None, gemitt_y=None, momentum_aperture_scale=0.85, ignored_portion=0.01, - seed=1997, nx=3, ny=3, nz=3, - theta_min=0.00005*np.pi, theta_max=0.99995*np.pi, **kwargs): + seed=1997, nx=3, ny=3, nz=3, **kwargs): # Input validation if line is None: @@ -286,15 +285,15 @@ def __init__(self, line=None, twiss=None, momentum_aperture=None, self.sigma_delta = sigma_delta self.bunch_population = bunch_population self.n_simulated = n_simulated - self.momentum_aperture_scale = momentum_aperture_scale self.ignored_portion = ignored_portion self.seed = seed self.nx = nx self.ny = ny self.nz = nz - self.theta_min = theta_min - self.theta_max = theta_max + # Limits from ELEGANT + self._theta_min = 0.00005*np.pi + self._theta_max = 0.99995*np.pi # Emittance validation nemitt_given = nemitt_x is not None and nemitt_y is not None From e2329eb0317990d63151cf20151c88407eeadf57 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Mon, 22 Sep 2025 13:56:22 +0200 Subject: [PATCH 29/37] =?UTF-8?q?Touschek:=20clamp=20=CE=B4=20sampling=20t?= =?UTF-8?q?o=20LMA=20(k=3D0.9).?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Prevented sampling of particles with |δ|>LMA in Touschek MC by reducing the longitudinal cutoff per element where needed: nz_eff = min(nz, 0.9 * min(|δN|, δP) / σδ) This ensures sampled δ stays strictly within the LMA, avoiding pathological small-angle events with huge weights that distorted RMC/RP and Touschek lifetime estimates. - Adjustment is local: only elements with tight acceptance are clamped - Prints a warning when nz_eff < nz - Restores stable RMC/RP≈1 and meaningful lifetime results --- xfields/touschek/manager.py | 56 +++++++++++++++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 2 deletions(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index c3721e5c..9df1cc1b 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -352,6 +352,58 @@ def _config(nn): dN = np.interp(s, momentum_aperture.s, momentum_aperture.deltan) dP = np.interp(s, momentum_aperture.s, momentum_aperture.deltap) + # Adjust the effective longitudinal sampling cutoff (nz_eff) to prevent + # generation of initial particles that are already outside + # the local momentum aperture (LMA) before Touschek scattering. + # + # Background: + # In the Touschek Monte Carlo routine, initial particle coordinates + # are drawn from a truncated Gaussian distribution with cutoffs + # {nx, ny, nz} in the transverse and longitudinal planes. The longitudinal + # cutoff nz sets the maximum |δ| ≈ nz*σδ. If nz*σδ exceeds the local + # momentum acceptance at this location, some particles + # are sampled outside the LMA (|δ| > LMA) even before any scattering event occurs. + # + # These particles create a pathological situation: + # • For very small scattering angles (θ* --> 0 in the CM frame), + # such particles are flagged as "candidates for loss" and passed to tracking, + # even though their state is essentially unchanged by scattering. + # • Because the Møller differential cross-section diverges at + # θ* --> 0, the corresponding particle weights become extremely large. + # • The pickPart routine then tends to select a handful of these + # high-weight pathological particles, distorting both the local + # scattering rate (RMC/RP diverges) and the overall Touschek + # lifetime estimate. + # + # Mitigation: + # To eliminate these spurious contributions, we dynamically reduce + # the longitudinal cutoff at each TouschekScattering element: + # + # nz_eff = min(nz, 0.9 * min(|δN|, δP) / σδ) + # + # where δN, δP are the negative/positive momentum aperture limits + # (scaled by momentum_aperture_scale). This ensures that the sampled + # longitudinal range ±nz_eff*σδ always lies strictly inside the local + # momentum aperture, with a small safety factor (0.9). As a result, + # only pathological large-weight events are avoided, and the Monte Carlo + # rate remains consistent with the Piwinski formula. + # + # NOTE: nz_eff is determined independently at each scattering element, + # so tighter cutoffs are applied only where the local momentum aperture + # is restrictive, while wider cutoffs are retained elsewhere. + min_dNdP = min(abs(dN), dP) + nz_eff = min(self.nz, 0.9 * min_dNdP / self.sigma_delta) + + if nz_eff < self.nz: + print(f""" + *********************************************************************************************** + [TouschekManager] Warning: longitudinal cutoff reduced at element '{nn}' (s={s:.2f} m). + + Using nz_eff={nz_eff:.2f} instead of nz={self.nz:.2f}. + This ensures that particles are sampled strictly within the local momentum aperture. + *********************************************************************************************** + """) + piwinski_rate = self.touschek._compute_piwinski_scattering_rate(nn) elem = line[nn] # xf.TouschekScattering @@ -370,8 +422,8 @@ def _config(nn): _sigma_z=self.sigma_z, _sigma_delta=self.sigma_delta, _n_simulated=self.n_simulated, - _nx=self.nx, _ny=self.ny, _nz=self.nz, - _theta_min=self.theta_min, _theta_max=self.theta_max, + _nx=self.nx, _ny=self.ny, _nz=nz_eff, + _theta_min=self._theta_min, _theta_max=self._theta_max, _ignored_portion=self.ignored_portion, piwinski_rate=piwinski_rate, _seed=self.seed, _inhibit_permute=0 From 6212f506e06acaf2780204069b77ea389e0fed68 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Mon, 22 Sep 2025 15:27:57 +0200 Subject: [PATCH 30/37] Touschek: safety factor from LMA 0.9 --> 0.85. --- xfields/touschek/manager.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index 9df1cc1b..3e169e07 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -379,12 +379,12 @@ def _config(nn): # To eliminate these spurious contributions, we dynamically reduce # the longitudinal cutoff at each TouschekScattering element: # - # nz_eff = min(nz, 0.9 * min(|δN|, δP) / σδ) + # nz_eff = min(nz, 0.85 * min(|δN|, δP) / σδ) # # where δN, δP are the negative/positive momentum aperture limits # (scaled by momentum_aperture_scale). This ensures that the sampled # longitudinal range ±nz_eff*σδ always lies strictly inside the local - # momentum aperture, with a small safety factor (0.9). As a result, + # momentum aperture, with a small safety factor (0.85). As a result, # only pathological large-weight events are avoided, and the Monte Carlo # rate remains consistent with the Piwinski formula. # @@ -392,7 +392,7 @@ def _config(nn): # so tighter cutoffs are applied only where the local momentum aperture # is restrictive, while wider cutoffs are retained elsewhere. min_dNdP = min(abs(dN), dP) - nz_eff = min(self.nz, 0.9 * min_dNdP / self.sigma_delta) + nz_eff = min(self.nz, 0.85 * min_dNdP / self.sigma_delta) if nz_eff < self.nz: print(f""" From b9044f5ed9ea7394fac2e73f47f1cf415ff0eadf Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Mon, 22 Sep 2025 17:48:39 +0200 Subject: [PATCH 31/37] Touschek: center Touschek scattered particles around the closed orbit. --- xfields/beam_elements/touschek.py | 28 +++++++++++++++++++++++++++- xfields/touschek/manager.py | 9 +++++++++ 2 files changed, 36 insertions(+), 1 deletion(-) diff --git a/xfields/beam_elements/touschek.py b/xfields/beam_elements/touschek.py index 00a3b644..6e80ee2c 100644 --- a/xfields/beam_elements/touschek.py +++ b/xfields/beam_elements/touschek.py @@ -69,9 +69,12 @@ class TouschekScattering(xt.BeamElement): def __init__(self, s=0.0, particle_ref=xt.Particles(), + element_index=0, bunch_population=0.0, alfx=0.0, betx=0.0, alfy=0.0, bety=0.0, dx=0.0, dpx=0.0, dy=0.0, dpy=0.0, + x_co=0.0, px_co=0.0, y_co=0.0, py_co=0.0, + zeta_co=0.0, delta_co=0.0, deltaN=0.0, deltaP=0.0, gemitt_x=0.0, gemitt_y=0.0, sigma_z=0.0, sigma_delta=0.0, @@ -96,6 +99,7 @@ def __init__(self, s=0.0, self._s = s self._particle_ref = particle_ref + self._element_index = element_index self._bunch_population = bunch_population self._alfx = alfx self._betx = betx @@ -105,6 +109,12 @@ def __init__(self, s=0.0, self._dpx = dpx self._dy = dy self._dpy = dpy + self._x_co = x_co + self._px_co = px_co + self._y_co = y_co + self._py_co = py_co + self._zeta_co = zeta_co + self._delta_co = delta_co self._deltaN = deltaN self._deltaP = deltaP self._gemitt_x = gemitt_x @@ -125,10 +135,13 @@ def __init__(self, s=0.0, def _configure(self, **kwargs): config_allowed = { - "_s", "_particle_ref", "_bunch_population", + "_s", "_particle_ref", "_element_index", + "_bunch_population", "_gemitt_x", "_gemitt_y", "_alfx", "_betx", "_alfy", "_bety", "_dx", "_dpx", "_dy", "_dpy", + "_x_co", "_px_co", "_y_co", "_py_co", + "_zeta_co", "_delta_co", "_deltaN", "_deltaP", "_sigma_z", "_sigma_delta", "_n_simulated", "_nx", "_ny", "_nz", @@ -189,6 +202,19 @@ def scatter(self): weight=weight_out[:n], s=getattr(self, '_s', 0.0)) + # Shift Touschek scattered particles around the closed orbit + part.x[:n] += self._x_co + part.px[:n] += self._px_co + part.y[:n] += self._y_co + part.py[:n] += self._py_co + part.zeta[:n] += self._zeta_co + + delta_temp = part.delta.copy() + delta_temp[:n] += self._delta_co + part.update_delta(delta_temp) + + part.at_element = self._element_index + part_ids = part.filter(part.state == 1).particle_id self.theta_log = dict(zip(part_ids.astype(int), theta_out[:n].astype(float))) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index 3e169e07..18039f05 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -352,6 +352,10 @@ def _config(nn): dN = np.interp(s, momentum_aperture.s, momentum_aperture.deltan) dP = np.interp(s, momentum_aperture.s, momentum_aperture.deltap) + x_co = twiss["x", nn]; px_co = twiss["px", nn] + y_co = twiss["y", nn]; py_co = twiss["py", nn] + zeta_co = twiss["zeta", nn]; delta_co = twiss["delta", nn] + # Adjust the effective longitudinal sampling cutoff (nz_eff) to prevent # generation of initial particles that are already outside # the local momentum aperture (LMA) before Touschek scattering. @@ -407,10 +411,12 @@ def _config(nn): piwinski_rate = self.touschek._compute_piwinski_scattering_rate(nn) elem = line[nn] # xf.TouschekScattering + element_index = line.element_names.index(nn) elem._configure( _s=s, _particle_ref=self.particle_ref, + _element_index=element_index, _bunch_population=self.bunch_population, _gemitt_x=self.gemitt_x, _gemitt_y=self.gemitt_y, @@ -418,6 +424,9 @@ def _config(nn): _alfy=alfy, _bety=bety, _dx=dx, _dpx=dpx, _dy=dy, _dpy=dpy, + _x_co=x_co, _px_co=px_co, + _y_co=y_co, _py_co=py_co, + _zeta_co=zeta_co, _delta_co=delta_co, _deltaN=dN, _deltaP=dP, _sigma_z=self.sigma_z, _sigma_delta=self.sigma_delta, From 60d333ff9e94ebd0d7bbb9b06f1de30b44bde637 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Mon, 22 Sep 2025 19:35:18 +0200 Subject: [PATCH 32/37] Touschek: remove ignoredRate. --- xfields/beam_elements/touschek_src/touschek.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h index e00fc04d..daf89fdd 100644 --- a/xfields/beam_elements/touschek_src/touschek.h +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -346,7 +346,6 @@ void TouschekScatter(TouschekScatteringData el, double totalWeight = 0.0; double totalMCRate = 0.0; - double ignoredRate = 0.0; double pTemp[6], p1[6], p2[6], dens1, dens2; double theta, phi, qa[3], qb[3], beta[3], qabs, gamma; @@ -489,7 +488,6 @@ void TouschekScatter(TouschekScatteringData el, } factor = factor / (double)(total_event); totalMCRate = totalWeight * factor; - ignoredRate = totalMCRate * ignoredPortion; /* Pick tracking particles from the simulated scattered particles */ index = (long *)malloc(sizeof(long) * simuCount); From 48568a31973563924bf24d2ab543dffa5b8158c5 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Wed, 24 Sep 2025 14:32:38 +0200 Subject: [PATCH 33/37] Restore original epsabs and epsrel in Piwinski integral. --- xfields/touschek/manager.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py index 18039f05..b073b486 100644 --- a/xfields/touschek/manager.py +++ b/xfields/touschek/manager.py @@ -56,8 +56,8 @@ def int_piwinski(k): int_piwinski, km, pi / 2, - epsabs=1e-12, - epsrel=1e-8 + epsabs=1e-16, + epsrel=1e-12 ) return val From 0de6fa2f0f6623040bdece521ffe0a821bdff671 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Fri, 13 Feb 2026 17:23:57 +0100 Subject: [PATCH 34/37] Use the include header facility of Xobjects --- xfields/beam_elements/touschek.py | 6 +----- xfields/beam_elements/touschek_src/touschek.h | 13 ++++++++----- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/xfields/beam_elements/touschek.py b/xfields/beam_elements/touschek.py index 6e80ee2c..801b379b 100644 --- a/xfields/beam_elements/touschek.py +++ b/xfields/beam_elements/touschek.py @@ -7,8 +7,6 @@ import xtrack as xt import numpy as np -from ..general import _pkg_root - class TouschekScattering(xt.BeamElement): _xofields = { @@ -44,9 +42,7 @@ class TouschekScattering(xt.BeamElement): _depends_on = [xt.RandomUniformAccurate] _extra_c_sources = [ - _pkg_root.joinpath('headers/constants.h'), - _pkg_root.joinpath('headers/elegant_rng.h'), - _pkg_root.joinpath('beam_elements/touschek_src/touschek.h'), + '#include "xfields/beam_elements/touschek_src/touschek.h"' ] _per_particle_kernels = { diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h index daf89fdd..2567d24e 100644 --- a/xfields/beam_elements/touschek_src/touschek.h +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -56,6 +56,9 @@ #ifndef XTRACK_TOUSCHEK_H #define XTRACK_TOUSCHEK_H +#include "xtrack/headers/track.h" +#include "xfields/headers/elegant_rng.h" + #include #include #include @@ -134,8 +137,8 @@ void bunch2cm(double *p1, double *p2, double *q, double *beta, double *gamma) { pp1 = pp1 + sqr(p1[i]); pp2 = pp2 + sqr(p2[i]); } - e1 = sqrt(MELECTRON_EV * MELECTRON_EV + pp1); - e2 = sqrt(MELECTRON_EV * MELECTRON_EV + pp2); + e1 = sqrt(ELECTRON_MASS_EV * ELECTRON_MASS_EV + pp1); + e2 = sqrt(ELECTRON_MASS_EV * ELECTRON_MASS_EV + pp2); ee = e1 + e2; betap1 = 0.0; @@ -193,7 +196,7 @@ void cm2bunch(double *p1, double *p2, double *q, double *beta, double *gamma) { pq = pq + q[i] * q[i]; } - e = sqrt(MELECTRON_EV * MELECTRON_EV + pq); + e = sqrt(ELECTRON_MASS_EV * ELECTRON_MASS_EV + pq); betaq = 0.0; bb = 0.0; @@ -354,7 +357,7 @@ void TouschekScatter(TouschekScatteringData el, double weight_limit, weight_ave, wTotal; const double sigxyz = sqrt(twissBeta[0]*gemitt[0]) * sqrt(twissBeta[1]*gemitt[1]) * sigma_z; - temp = sqr(bunch_population) * sqr(PI) * sqr(RE) * C_LIGHT / 4.; + temp = sqr(bunch_population) * sqr(PI) * sqr(RADIUS_ELECTRON) * C_LIGHT / 4.; double factor = temp * pow(range[0], 3.0) * pow(range[1], 3.0) * pow(range[2], 3.0) / pow(2 * PI, 6.0) / sigxyz; double *xtemp = (double*)malloc(sizeof(double) * n_simulated); @@ -444,7 +447,7 @@ void TouschekScatter(TouschekScatteringData el, } if (p1[5] < deltaN || p2[5] > deltaP) { - beta0 = qabs / sqrt(qabs * qabs + MELECTRON_EV * MELECTRON_EV); + beta0 = qabs / sqrt(qabs * qabs + ELECTRON_MASS_EV * ELECTRON_MASS_EV); cross = moeller(beta0, theta); temp *= cross * beta0 / gamma / gamma; From 929efcd5fff62fd92c015b986c5f60f079538b1c Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Fri, 13 Feb 2026 17:24:19 +0100 Subject: [PATCH 35/37] Touschek example! --- .../006_touschek/000_touschek_toy_ring.py | 228 ++++++++++++++++++ 1 file changed, 228 insertions(+) create mode 100644 examples/006_touschek/000_touschek_toy_ring.py diff --git a/examples/006_touschek/000_touschek_toy_ring.py b/examples/006_touschek/000_touschek_toy_ring.py new file mode 100644 index 00000000..3751957b --- /dev/null +++ b/examples/006_touschek/000_touschek_toy_ring.py @@ -0,0 +1,228 @@ +# copyright ############################### # +# This file is part of the Xtrack Package. # +# Copyright (c) CERN, 2025. # +# ######################################### # +import numpy as np +import matplotlib.pyplot as plt + +import xobjects as xo +import xtrack as xt +import xfields as xf + +###################################################### +# Beam parameters +###################################################### +nemitt_x = 1e-5 +nemitt_y = 1e-7 + +sigma_z = 4e-3 +sigma_delta = 1e-3 + +bunch_population = 4e9 + +###################################################### +# Build a toy ring +###################################################### +lbend = 3 +angle = np.pi / 2 + +lquad = 0.3 +k1qf = 0.1 +k1qd = 0.7 + +# Create environment +env = xt.Environment() + +# Define the line (toy ring) +line = env.new_line(components=[ + env.new('mqf.1', xt.Quadrupole, length=lquad, k1=k1qf), + env.new('d1.1', xt.Drift, length=1), + env.new('mb1.1', xt.Bend, length=lbend, angle=angle), + env.new('d2.1', xt.Drift, length=1), + + env.new('mqd.1', xt.Quadrupole, length=lquad, k1=-k1qd), + env.new('d3.1', xt.Drift, length=1), + env.new('mb2.1', xt.Bend, length=lbend, angle=angle), + env.new('d4.1', xt.Drift, length=1), + + env.new('mqf.2', xt.Quadrupole, length=lquad, k1=k1qf), + env.new('d1.2', xt.Drift, length=1), + env.new('mb1.2', xt.Bend, length=lbend, angle=angle), + env.new('d2.2', xt.Drift, length=1), + + env.new('mqd.2', xt.Quadrupole, length=lquad, k1=-k1qd), + env.new('d3.2', xt.Drift, length=1), + env.new('mb2.2', xt.Bend, length=lbend, angle=angle), + env.new('d4.2', xt.Drift, length=1), +]) + +# Set the reference particle +line.set_particle_ref('electron', p0c=1e9) + +# Configure the bend model +line.configure_bend_model(core='full', edge=None) + +###################################################### +# Insert Touschek scattering centers +###################################################### +# We insert Touschek scattering centers in the middle of each magnet +# to have good coverage of variations of the optical functions +tab = line.get_table() +tab_bends_quads = tab.rows[(tab.element_type == 'Bend') | (tab.element_type == 'Quadrupole')] + +# Would be good to have env.new for xf.TouschekScattering +for ii, nn in enumerate(tab_bends_quads.name): + tscatter_name = f'TScatter_{ii}' + env.elements[tscatter_name] = xf.TouschekScattering() + line.insert(tscatter_name, at=0.0, from_=nn) + +# The last TouschekScattering element has to be placed at the end of the line +tscatter_name = f'TScatter_{ii+1}' +env.elements[tscatter_name] = xf.TouschekScattering() +line.insert(tscatter_name, at=tab.s[-1]) + + +###################################################### +# Install apertures +###################################################### +tab = line.get_table() +needs_aperture = np.unique(tab.element_type)[ + ~np.isin(np.unique(tab.element_type), ["", "Drift", "Marker"]) +] + +aper_size = 0.040 # m + +placements = [] +for nn, ee in zip(tab.name, tab.element_type): + if ee not in needs_aperture: + continue + + env.new( + f'{nn}_aper_entry', xt.LimitRect, + min_x=-aper_size, max_x=aper_size, + min_y=-aper_size, max_y=aper_size + ) + placements.append(env.place(f'{nn}_aper_entry', at=f'{nn}@start')) + + env.new( + f'{nn}_aper_exit', xt.LimitRect, + min_x=-aper_size, max_x=aper_size, + min_y=-aper_size, max_y=aper_size + ) + placements.append(env.place(f'{nn}_aper_exit', at=f'{nn}@end')) + +line.insert(placements) + +###################################################### +# Evaluate momentum aperture profile +###################################################### +# Norlamized emittance +nemitt_x = 1e-5 +nemitt_y = 1e-7 + +# Evaluate local momentum aperture at the touschek scattering centers +momentum_aperture = line.momentum_aperture( + # twiss=tw, + include_type_pattern="TouschekScattering", + nemitt_x=nemitt_x, + nemitt_y=nemitt_y, + y_offset=1e-9, + delta_negative_limit=-0.012, + delta_positive_limit=0.012, + delta_step_size=1e-4, + n_turns=1000, + method="4d" +) + +df_momentum_aperture = momentum_aperture.to_pandas() + +###################################################### +# Plot +###################################################### +plt.plot(momentum_aperture.s, momentum_aperture.deltan*100, c='r') +plt.plot(momentum_aperture.s, momentum_aperture.deltap*100, c='r') +plt.title('Toy ring: local momentum aperture profile') +plt.xlabel('s [m]') +plt.ylabel(r'$\delta$ [%]') +plt.grid() +plt.show() + +###################################################### +# Touschek simulation +###################################################### +# Parameters +momentum_aperture_scale = 0.85 # scaling factor for momentum aperture +n_simulated = 5e6 # number of simulated scattering events with delta > delta_min +nturns = 1000 # number of turns to simulate + +touschek_manager = xf.TouschekManager( + line, + momentum_aperture=df_momentum_aperture, + momentum_aperture_scale=momentum_aperture_scale, + nemitt_x=nemitt_x, + nemitt_y=nemitt_y, + sigma_z=sigma_z, + sigma_delta=sigma_delta, + bunch_population=bunch_population, + n_simulated=n_simulated, + nx=3, ny=3, nz=3, + ignored_portion=0.01, + seed=1997, + method='4d' +) + +touschek_manager.initialise_touschek() + +touschek_elements = tab.rows[tab.element_type == 'TouschekScattering'].name + +line.discard_tracker() +line.build_tracker(_context=xo.ContextCpu(omp_num_threads='auto')) + +particles_list = [] +for ii in range(len(touschek_elements)): + element = touschek_elements[ii] # xf.TouschekScattering + s_start_elem = tab.rows[tab.name == element].s[0] + + # Touschek! + particles = line[element].scatter() + + # Track! + print(f"\nTrack particles scattered at {element} at s = {s_start_elem}") + line.track(particles, ele_start=element, ele_stop=element, num_turns=nturns, with_progress=1) + + particles_list.append(particles) + +particles = xt.Particles.merge(particles_list) + +# Refine loss location +loss_loc_refinement = xt.LossLocationRefinement(line, + n_theta = 360, # Angular resolution in the polygonal approximation of the aperture + r_max = 0.5, # Maximum transverse aperture in m + dr = 50e-6, # Transverse loss refinement accuracy [m] + ds = 0.1, # Longitudinal loss refinement accuracy [m] + ) + +loss_loc_refinement.refine_loss_location(particles) + +###################################################### +# Compute Touschek lifetime +###################################################### +# Keep lost particles only +particles = particles.filter(particles.state == 0) +# Compute total Touschek loss rate +loss_rate = sum(particles.weight) +# Compute Touschek lifetime +touschek_lifetime = bunch_population / loss_rate + +###################################################### +# Plot: Toy ring Touschek loss map +###################################################### +circumference = line.get_length() +binwidth = 0.1 # m + +plt.title(f'Toy ring Touschek loss map (Touschek lifetime: {touschek_lifetime/60:.2f} min)') +plt.hist(particles.s, bins=np.arange(0, circumference + binwidth, binwidth), weights=particles.weight*1e-3) +plt.xlabel('s [m]') +plt.ylabel('Loss rate [kHz]') +plt.grid() +plt.show() \ No newline at end of file From f1a473db07fc0d2aae0915c8396e8351a6345f10 Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Fri, 13 Feb 2026 17:27:07 +0100 Subject: [PATCH 36/37] Clean up Touschek example. --- examples/006_touschek/000_touschek_toy_ring.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/examples/006_touschek/000_touschek_toy_ring.py b/examples/006_touschek/000_touschek_toy_ring.py index 3751957b..c53fa702 100644 --- a/examples/006_touschek/000_touschek_toy_ring.py +++ b/examples/006_touschek/000_touschek_toy_ring.py @@ -70,7 +70,6 @@ tab = line.get_table() tab_bends_quads = tab.rows[(tab.element_type == 'Bend') | (tab.element_type == 'Quadrupole')] -# Would be good to have env.new for xf.TouschekScattering for ii, nn in enumerate(tab_bends_quads.name): tscatter_name = f'TScatter_{ii}' env.elements[tscatter_name] = xf.TouschekScattering() @@ -81,7 +80,6 @@ env.elements[tscatter_name] = xf.TouschekScattering() line.insert(tscatter_name, at=tab.s[-1]) - ###################################################### # Install apertures ###################################################### @@ -116,10 +114,6 @@ ###################################################### # Evaluate momentum aperture profile ###################################################### -# Norlamized emittance -nemitt_x = 1e-5 -nemitt_y = 1e-7 - # Evaluate local momentum aperture at the touschek scattering centers momentum_aperture = line.momentum_aperture( # twiss=tw, From 728a7be9c03cc45f771e957625ef524d95cb533f Mon Sep 17 00:00:00 2001 From: gbrogginess Date: Fri, 13 Feb 2026 17:53:53 +0100 Subject: [PATCH 37/37] Touschek test draft. --- tests/test_touschek.py | 263 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 263 insertions(+) create mode 100644 tests/test_touschek.py diff --git a/tests/test_touschek.py b/tests/test_touschek.py new file mode 100644 index 00000000..c0c8d062 --- /dev/null +++ b/tests/test_touschek.py @@ -0,0 +1,263 @@ +# copyright ################################# # +# This file is part of the Xfields Package. # +# Copyright (c) CERN, 2025. # +# ########################################### # +import numpy as np +import pytest + +import xobjects as xo +import xtrack as xt +import xfields as xf + +# --- Helpers ----------------------------------------------------------------- +def _build_toy_ring_with_touschek_and_apertures( + *, + aper_size=0.040, + install_apertures=True, +): + lbend = 3.0 + angle = np.pi / 2 + + lquad = 0.3 + k1qf = 0.1 + k1qd = 0.7 + + env = xt.Environment() + line = env.new_line(components=[ + env.new('mqf.1', xt.Quadrupole, length=lquad, k1=k1qf), + env.new('d1.1', xt.Drift, length=1.0), + env.new('mb1.1', xt.Bend, length=lbend, angle=angle), + env.new('d2.1', xt.Drift, length=1.0), + + env.new('mqd.1', xt.Quadrupole, length=lquad, k1=-k1qd), + env.new('d3.1', xt.Drift, length=1.0), + env.new('mb2.1', xt.Bend, length=lbend, angle=angle), + env.new('d4.1', xt.Drift, length=1.0), + + env.new('mqf.2', xt.Quadrupole, length=lquad, k1=k1qf), + env.new('d1.2', xt.Drift, length=1.0), + env.new('mb1.2', xt.Bend, length=lbend, angle=angle), + env.new('d2.2', xt.Drift, length=1.0), + + env.new('mqd.2', xt.Quadrupole, length=lquad, k1=-k1qd), + env.new('d3.2', xt.Drift, length=1.0), + env.new('mb2.2', xt.Bend, length=lbend, angle=angle), + env.new('d4.2', xt.Drift, length=1.0), + ]) + + line.set_particle_ref('electron', p0c=1e9) + line.configure_bend_model(core='full', edge=None) + + # Insert Touschek scattering centers (one per magnet, plus one at end) + tab0 = line.get_table() + tab_magnets = tab0.rows[(tab0.element_type == 'Bend') | (tab0.element_type == 'Quadrupole')] + + for ii, nn in enumerate(tab_magnets.name): + name = f'TScatter_{ii}' + env.elements[name] = xf.TouschekScattering() + line.insert(name, at=0.0, from_=nn) + + # Last TouschekScattering at end of the line + name_last = f'TScatter_{ii + 1}' + env.elements[name_last] = xf.TouschekScattering() + line.insert(name_last, at=float(line.get_length())) + + # Install rectangular apertures + if install_apertures: + tab = line.get_table() + needs_aperture = np.unique(tab.element_type)[ + ~np.isin(np.unique(tab.element_type), ["", "Drift", "Marker"]) + ] + + placements = [] + for nn, ee in zip(tab.name, tab.element_type): + if ee not in needs_aperture: + continue + + env.new( + f'{nn}_aper_entry', xt.LimitRect, + min_x=-aper_size, max_x=aper_size, + min_y=-aper_size, max_y=aper_size + ) + placements.append(env.place(f'{nn}_aper_entry', at=f'{nn}@start')) + + env.new( + f'{nn}_aper_exit', xt.LimitRect, + min_x=-aper_size, max_x=aper_size, + min_y=-aper_size, max_y=aper_size + ) + placements.append(env.place(f'{nn}_aper_exit', at=f'{nn}@end')) + + line.insert(placements) + + return env, line + + +def _compute_fast_lma_at_touschek_centers(line, *, nemitt_x, nemitt_y): + # Fast-ish settings for tests + return line.momentum_aperture( + include_type_pattern="TouschekScattering", + nemitt_x=nemitt_x, + nemitt_y=nemitt_y, + y_offset=1e-12, + delta_negative_limit=-0.012, + delta_positive_limit=+0.012, + delta_step_size=0.001, + n_turns=64, + method="4d", + with_progress=False, + verbose=False, + ).to_pandas() + + +# --- Tests ------------------------------------------------------------------- +def test_touschek_manager_initialise_configures_elements(): + """ + Check that TouschekManager.initialise_touschek() configures all TouschekScattering + elements with consistent optics, momentum aperture, and positive rates. + """ + nemitt_x = 1e-5 + nemitt_y = 1e-7 + sigma_z = 4e-3 + sigma_delta = 1e-3 + bunch_population = 4e9 + + _, line = _build_toy_ring_with_touschek_and_apertures() + + # LMA at Touschek centers (DataFrame-like as expected by TouschekManager) + df_ma = _compute_fast_lma_at_touschek_centers(line, nemitt_x=nemitt_x, nemitt_y=nemitt_y) + + # Keep runtime reasonable + tm = xf.TouschekManager( + line, + momentum_aperture=df_ma, + momentum_aperture_scale=0.85, + nemitt_x=nemitt_x, + nemitt_y=nemitt_y, + sigma_z=sigma_z, + sigma_delta=sigma_delta, + bunch_population=bunch_population, + n_simulated=1e6, + nx=3, ny=3, nz=3, + ignored_portion=0.01, + seed=1997, + method="4d", + ) + + tm.initialise_touschek() + + tab = line.get_table() + tnames = tab.rows[tab.element_type == "TouschekScattering"].name + assert len(tnames) > 0 + + for nn in tnames: + el = line[nn] + assert isinstance(el, xf.TouschekScattering) + + # Config fields should be set + assert np.isfinite(el._p0c) and el._p0c > 0 + assert np.isfinite(el._integrated_piwinski_rate) and el._integrated_piwinski_rate >= 0 + assert np.isfinite(el.piwinski_rate) and el.piwinski_rate >= 0 + + # Acceptance should have correct sign convention + assert el._deltaN <= 0 + assert el._deltaP >= 0 + + # Consistency + assert el._bunch_population == pytest.approx(bunch_population) + assert el._sigma_z == pytest.approx(sigma_z) + assert el._sigma_delta == pytest.approx(sigma_delta) + + # nz may be reduced by nz_eff logic, but never increased + assert el._nz <= 3.0 + 1e-15 + + +def test_touschek_scattering_scatter_returns_particles(): + """ + Run scatter() on one TouschekScattering element and check that: + - returned object is xt.Particles, + - some particles are selected, + - weights and deltas are finite, + - at_element points to the TouschekScattering location. + """ + nemitt_x = 1e-5 + nemitt_y = 1e-7 + sigma_z = 4e-3 + sigma_delta = 1e-3 + bunch_population = 4e9 + + _, line = _build_toy_ring_with_touschek_and_apertures() + + df_ma = _compute_fast_lma_at_touschek_centers(line, nemitt_x=nemitt_x, nemitt_y=nemitt_y) + + tm = xf.TouschekManager( + line, + momentum_aperture=df_ma, + momentum_aperture_scale=0.85, + nemitt_x=nemitt_x, + nemitt_y=nemitt_y, + sigma_z=sigma_z, + sigma_delta=sigma_delta, + bunch_population=bunch_population, + n_simulated=1e6, + nx=3, ny=3, nz=3, + ignored_portion=0.01, + seed=1997, + method="4d", + ) + tm.initialise_touschek() + + tab = line.get_table() + tnames = tab.rows[tab.element_type == "TouschekScattering"].name + nn = tnames[0] + el = line[nn] + + parts = el.scatter() + assert isinstance(parts, xt.Particles) + + # "Selected" particles are those with state==1 immediately after scatter + alive = parts.filter(parts.state == 1) + assert alive._capacity > 0 + assert len(alive.x) > 0 + + # Sanity of generated coordinates and weights + assert np.all(np.isfinite(alive.x)) + assert np.all(np.isfinite(alive.px)) + assert np.all(np.isfinite(alive.y)) + assert np.all(np.isfinite(alive.py)) + assert np.all(np.isfinite(alive.zeta)) + assert np.all(np.isfinite(alive.delta)) + assert np.all(np.isfinite(alive.weight)) + assert np.all(alive.weight >= 0) + + # Particles should be located at the element index + expected_idx = line.element_names.index(nn) + assert int(alive.at_element[0]) == expected_idx + + # Element should have recorded MC totals + assert np.isfinite(el.total_mc_rate) + assert el.total_mc_rate >= 0 + + +def test_touschek_manager_validates_momentum_aperture_columns(): + """ + TouschekManager should reject malformed momentum_aperture objects. + """ + _, line = _build_toy_ring_with_touschek_and_apertures() + + # Minimal fake DataFrame-like with missing columns + import pandas as pd + bad = pd.DataFrame({"s": [0.0], "deltan": [-0.01]}) # missing 'deltap' + + with pytest.raises(ValueError, match=r"missing columns"): + xf.TouschekManager( + line, + momentum_aperture=bad, + sigma_z=4e-3, + sigma_delta=1e-3, + bunch_population=1.0, + n_simulated=10, + nemitt_x=1e-6, + nemitt_y=1e-6, + method="4d", + )