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 diff --git a/contributors_and_credits.txt b/contributors_and_credits.txt index bf72b891..3af431d3 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 @@ -18,3 +19,8 @@ External contributors: The Particle In Cell implementation is largely based on the PyPIC package (https://github.com/pycomplete/pypic) + +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 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..c53fa702 --- /dev/null +++ b/examples/006_touschek/000_touschek_toy_ring.py @@ -0,0 +1,222 @@ +# 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')] + +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 +###################################################### +# 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 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", + ) 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 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/__init__.py b/xfields/__init__.py index aa9c919c..daaf3db3 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 .config_tools import replace_spacecharge_with_quasi_frozen from .config_tools import replace_spacecharge_with_PIC diff --git a/xfields/beam_elements/beambeam3d.py b/xfields/beam_elements/beambeam3d.py index 0fa2dfbd..69a9495a 100644 --- a/xfields/beam_elements/beambeam3d.py +++ b/xfields/beam_elements/beambeam3d.py @@ -180,6 +180,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 @@ -234,6 +235,7 @@ def __init__(self, flag_bhabha=0, compt_x_min=1e-4, + compt_scale=1., flag_beamsize_effect=1, flag_luminosity=0, @@ -313,6 +315,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. @@ -448,7 +451,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) @@ -618,9 +621,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 @@ -632,6 +636,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 fd444ce1..043206fa 100644 --- a/xfields/beam_elements/beambeam_src/beambeam3d.h +++ b/xfields/beam_elements/beambeam_src/beambeam3d.h @@ -114,9 +114,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/beam_elements/touschek.py b/xfields/beam_elements/touschek.py new file mode 100644 index 00000000..801b379b --- /dev/null +++ b/xfields/beam_elements/touschek.py @@ -0,0 +1,223 @@ +# 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 + +class TouschekScattering(xt.BeamElement): + + _xofields = { + '_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, + '_theta_min': xo.Float64, + '_theta_max': xo.Float64, + '_ignored_portion': xo.Float64, + '_integrated_piwinski_rate': xo.Float64, + '_seed': xo.Int64, + '_inhibit_permute': xo.Int64 + } + + # allow_track = False + _depends_on = [xt.RandomUniformAccurate] + + _extra_c_sources = [ + '#include "xfields/beam_elements/touschek_src/touschek.h"' + ] + + _per_particle_kernels = { + '_scatter': xo.Kernel( + c_name='TouschekScatter', + args=[ + 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='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), + ], + ), + } + + 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, + 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, + seed=1997, + inhibit_permute=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._element_index = element_index + 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._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 + 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._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 + self._seed = seed + self._inhibit_permute = inhibit_permute + + def _configure(self, **kwargs): + config_allowed = { + "_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", + "_theta_min", "_theta_max", + "_ignored_portion", "piwinski_rate", + "_integrated_piwinski_rate", + "_seed", "_inhibit_permute" + } + + 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) + 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) + + 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, + theta_out=theta_out, + weight_out=weight_out, + totalMCRate_out=totalMCRate_out, + 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, + 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)) + + # 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))) + + 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 diff --git a/xfields/beam_elements/touschek_src/touschek.h b/xfields/beam_elements/touschek_src/touschek.h new file mode 100644 index 00000000..2567d24e --- /dev/null +++ b/xfields/beam_elements/touschek_src/touschek.h @@ -0,0 +1,579 @@ +/*************************************************************************\ +* 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 + +#include "xtrack/headers/track.h" +#include "xfields/headers/elegant_rng.h" + +#include +#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; +} + +/* Adapted from Elegant touschekScatter.c: selectPartGauss (logic unchanged). */ +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; +} + +/* 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; + 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(ELECTRON_MASS_EV * ELECTRON_MASS_EV + pp1); + e2 = sqrt(ELECTRON_MASS_EV * ELECTRON_MASS_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 */ +/* 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; + + *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; +} + +/* 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; + + pq = 0.0; + for (i = 0; i < 3; i++) { + pq = pq + q[i] * q[i]; + } + + e = sqrt(ELECTRON_MASS_EV * ELECTRON_MASS_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; +} + +/* From Elegant touschekScatter.c: moeller */ +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; +} + +/* 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; + 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; +} + +/* Adapted from Elegant touschekScatter.c: TouschekDistribution (logic unchanged) */ +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* theta_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 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); + + 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 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(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); + 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); + double *thetatemp = (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; + simuCount = 0; + + while (1) { + if (i >= n_simulated) + 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); // Does not match with ELEGANT + ran1[j] = random_1_elegant(1); + } + randomizeOrder((char*)ran1, sizeof(ran1[0]), 11, 0, random_4); // like ELEGANT + + total_event++; + + selectPartGauss(p1, p2, &dens1, &dens2, ran1, range, twissAlpha, twissBeta, twissDisp, gemitt); + + if (!dens1 || !dens2) { + continue; + } + /* 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; + } + // 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); + + theta = theta_min + (theta_max - theta_min) * ran1[9]; + 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; + // 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++) { + pTemp[j] = p2[j]; + p2[j] = p1[j]; + p1[j] = pTemp[j]; + } + } + + if (p1[5] < deltaN || p2[5] > deltaP) { + beta0 = qabs / sqrt(qabs * qabs + ELECTRON_MASS_EV * ELECTRON_MASS_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]; + thetatemp[i] = theta; + 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]; + thetatemp[i] = theta; + weight[i] = temp; + i++; + } + } + } + factor = factor / (double)(total_event); + totalMCRate = totalWeight * factor; + + /* 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]; + theta_out[k] = thetatemp[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]; + theta_out[k] = thetatemp[src]; + weight_out[k] = weight[src]; + } + } + + // 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); + + // 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(thetatemp); + 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/bhabha_spectrum.h b/xfields/headers/bhabha_spectrum.h index 93358f32..63455fa0 100644 --- a/xfields/headers/bhabha_spectrum.h +++ b/xfields/headers/bhabha_spectrum.h @@ -252,6 +252,7 @@ GPUFUN 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) @@ -286,7 +287,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 diff --git a/xfields/headers/constants.h b/xfields/headers/constants.h index 591e50f5..7f988af9 100644 --- a/xfields/headers/constants.h +++ b/xfields/headers/constants.h @@ -16,6 +16,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/headers/elegant_rng.h b/xfields/headers/elegant_rng.h new file mode 100644 index 00000000..02a141d9 --- /dev/null +++ b/xfields/headers/elegant_rng.h @@ -0,0 +1,350 @@ +/*************************************************************************\ +* 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. +\*************************************************************************/ +/* + * 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). + */ +#ifndef ELEGANT_RNG_H +#define ELEGANT_RNG_H + +#include +#include +#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; + 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 (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; + 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; +} + +/* ----------------------------------------------------------------------------- + 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); + // 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); +} + +/* ----------------------------------------------------------------------------- + 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); +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 + // SDDS-like reseed + random_2(-(base + 2)); + random_3(-(base + 4)); + random_4(-(base + 6)); + random_5(-(base + 8)); + random_6(-(base + 10)); + // 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; + 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); +} + +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. + + 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; +} 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 || !urandom) return 0; + if (length < 2) return 1; + 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 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), s1 = labs(seed + 2), s2 = labs(seed + 4), s3 = labs(seed + 6); + + if (s0 == 987654321) inhibitRandomSeedPermutation(1); + else inhibitRandomSeedPermutation(inhibit_permute ? 1 : 0); + + random_1_elegant(-s0); + random_2(-s1); + random_3(-s2); + random_4(-s3); +} + +#endif // ELEGANT_RNG_H \ No newline at end of file diff --git a/xfields/touschek/README.md b/xfields/touschek/README.md new file mode 100644 index 00000000..aa5570b2 --- /dev/null +++ b/xfields/touschek/README.md @@ -0,0 +1,25 @@ +# 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). + +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`. + +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**. +© 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 diff --git a/xfields/touschek/manager.py b/xfields/touschek/manager.py new file mode 100644 index 00000000..b073b486 --- /dev/null +++ b/xfields/touschek/manager.py @@ -0,0 +1,458 @@ +# copyright ################################# # +# This file is part of the Xfields Package. # +# Copyright (c) CERN, 2021. # +# ########################################### # + +import xtrack as xt +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 + +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): + """ + 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 + """ + from math import atan, tan, sqrt, exp, log, pi + + km = atan(sqrt(tm)) + + def int_piwinski(k): + t = np.tan(k) ** 2 + fact = ( + (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 * exp(-B1*t) * i0(B2*t) * sqrt(1+t) + else: + intp = ( + fact + * exp(B2*t - B1*t) + / sqrt(2*pi * B2*t) + * sqrt(1+t) + ) + return intp + + val, _ = quad( + int_piwinski, + km, + pi / 2, + epsabs=1e-16, + epsrel=1e-12 + ) + + 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 + momentum_aperture = self.manager.momentum_aperture + gemitt_x = self.manager.gemitt_x + gemitt_y = self.manager.gemitt_y + 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 = twiss['delta', element] + dx = twiss['dx', element] + dpx = twiss['dpx', element] + dxt = alfx * dx + betx * dpx # dxt: dx tilde + dy = twiss['dy', element] + dpy = twiss['dpy', element] + dyt = alfy * dy + bety * dpy # dyt: dy tilde + + try: + s = self.twiss.rows[element].s[0] + except: + s = self.manager.line.get_s_position(element) + + 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 + + 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, element): + """ + 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. + """ + 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) + + # Indexes of the TouschekScatterings + ii_t = [ii for ii, nn in enumerate(tab.name[:-1]) if isinstance(line[nn], xf.TouschekScattering)] + + integrated = 0.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 + + 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=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, + seed=1997, 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", "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', '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 + self.twiss = twiss + + # 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 + momentum_aperture = ap + self.momentum_aperture = momentum_aperture + + self.sigma_z = sigma_z + self.sigma_delta = sigma_delta + self.bunch_population = bunch_population + self.n_simulated = n_simulated + self.ignored_portion = ignored_portion + self.seed = seed + self.nx = nx + self.ny = ny + self.nz = nz + + # 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 + 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) + + 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) + + # Pass the twiss to the TouschekCalculator + self.touschek.twiss = self.twiss + + 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): + 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, 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. + # + # 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.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.85). 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.85 * 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 + 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, + _alfx=alfx, _betx=betx, + _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, + _n_simulated=self.n_simulated, + _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 + ) + + 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): + 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__})." + ) + print(f'Initialising TouschekScattering for {element}') + _config(element) \ No newline at end of file