From 32e8f57747d8a752cb6ab6ae0dcfc9262a4333df Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 1 Sep 2026 22:55:30 -0700 Subject: [PATCH 1/5] Add an ImpactX interface ImpactX is an s-based beam dynamics code, the successor of IMPACT-Z, and openPMD-beamphysics had no interface for it. This adds one alongside the existing astra/bmad/elegant/gpt/impact interfaces: - beamphysics/interfaces/impactx.py: converters in both directions, plus a reader for ImpactX' openPMD BeamMonitor output that needs no ImpactX installation (openpmd-api is imported lazily). - ParticleGroup.to_impactx() / .from_impactx(), mirroring to_bmad/from_bmad. - Tests, and a small BeamMonitor fixture written by ImpactX 26.08 with the generate.py script beside it, so the reader is covered in CI without ImpactX. Coordinates: ImpactX holds particles at fixed s, so the result is in z-coordinates -- all z equal, the bunch length a spread in t -- as bmad.py does. The transverse coordinates stay in the local frame relative to the reference particle, since adding x_ref would be wrong wherever the reference orbit bends, while t is absolute lab time, which is unambiguous. A few ImpactX specifics are handled rather than papered over: t is c*t in metres there; qm is reported in different units depending on how the particles were inserted, so it is never trusted on read; the openPMD id record holds AMReX' packed idcpu, whose per-rank id field alone is not unique, so the whole value becomes the ParticleGroup id with only the validity bit moved to status; and every particle in a particles_lost file is lost by construction even though AMReX marks them valid, so they are not reported as alive. Spin and runtime per-particle components have no ParticleGroup representation and are refused rather than silently dropped, with strict=False to opt out. Verified against a live ImpactX 26.08 run: a bunch injected through add_n_particles and read back out of the monitor file matches to 2e-16 of the spread in x, y and t and 2.4e-13 in pz, with identical emittances; and the reader reproduces ImpactX' own reduced beam characteristics to 1e-12. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01DhbdaEDSBqFN974t9QcXem --- beamphysics/interfaces/impactx.py | 907 +++++++++++++++++++ beamphysics/particles.py | 87 +- docs/examples/data/impactx/generate.py | 76 ++ docs/examples/data/impactx/monitor.h5 | Bin 0 -> 130410 bytes docs/examples/data/impactx/particles_lost.h5 | Bin 0 -> 46655 bytes environment.yml | 1 + pyproject.toml | 9 +- tests/test_impactx.py | 592 ++++++++++++ 8 files changed, 1670 insertions(+), 2 deletions(-) create mode 100644 beamphysics/interfaces/impactx.py create mode 100755 docs/examples/data/impactx/generate.py create mode 100644 docs/examples/data/impactx/monitor.h5 create mode 100644 docs/examples/data/impactx/particles_lost.h5 create mode 100644 tests/test_impactx.py diff --git a/beamphysics/interfaces/impactx.py b/beamphysics/interfaces/impactx.py new file mode 100644 index 00000000..35f4a0a3 --- /dev/null +++ b/beamphysics/interfaces/impactx.py @@ -0,0 +1,907 @@ +"""ImpactX beam data <-> ParticleGroup. + +[ImpactX](https://impactx.readthedocs.io) is an s-based beam dynamics code, the +successor of IMPACT-Z. Its particles are held at a common ``s`` with a spread in +arrival time, which is z-coordinates on this side -- all ``z`` equal, ``t`` varying -- +so the conversion is a direct algebraic map, like the Bmad interface and unlike the +time-based ones. + +Coordinates and frames +---------------------- +ImpactX describes each particle at fixed ``s`` by ``(x, y, t, px, py, pt)``: + +- ``x``, ``y`` [m] are the transverse displacement from the reference particle, in the + local (curvilinear) frame that follows the reference orbit. +- ``t`` [m] is ``c`` times the difference between the particle's and the reference + particle's arrival time, i.e. a length, not a time. +- ``px``, ``py``, ``pt`` are dimensionless, normalized by the magnitude of the + reference momentum: ``px = Delta(beta_x gamma) / (beta_0 gamma_0)`` and + ``pt = -Delta(gamma) / (beta_0 gamma_0)``. + +`ParticleGroup` is a lab-frame container, so the mapping has to choose a frame: + +- The transverse coordinates stay in the **local frame**: ``x`` and ``y`` are the + displacement from the reference particle and ``z`` is zero, the reference plane. + Adding ``x_ref``/``z_ref`` would be wrong wherever the reference orbit is bent, + because local ``x`` is then not lab ``x``. Use `ImpactXRefPart` (``x``, ``y``, ``z``, + ``s``, ``px``, ``py``, ``pz``) if you need to place the bunch in the lab. +- The time is **absolute**: ``t = t_ref + position_t / c``. That one is unambiguous, it + is what openPMD's ``position/t + positionOffset/t`` means in ImpactX output, and it + keeps quantities like `ParticleGroup.average_current` meaningful. + +Originally developed in lume-impactx. +""" + +from __future__ import annotations + +import pathlib +import warnings +from dataclasses import dataclass +from typing import Any + +import numpy as np + +from ..particles import ParticleGroup +from ..species import charge_of, e_charge, mass_of +from ..status import ParticleStatus +from ..units import c_light + + +__all__ = [ + "ImpactXRefPart", + "PARTICLE_STATUS_LOST", + "UnrepresentableParticleData", + "beam_monitor_iterations", + "impactx_to_particlegroup_data", + "particle_id_from_idcpu", + "particlegroup_to_impactx", + "pmd_species_of", + "read_beam_monitor", + "read_beam_monitor_data", + "refpart_from_openpmd", +] + + +#: ImpactX's built-in species names mapped to openPMD-beamphysics names. +#: ImpactX only knows these four; anything else needs explicit mass and charge. +IMPACTX_TO_PMD_SPECIES = { + "electron": "electron", + "positron": "positron", + "proton": "proton", + "Hminus": "H-", +} +PMD_TO_IMPACTX_SPECIES = {v: k for k, v in IMPACTX_TO_PMD_SPECIES.items()} + + +# -------------------------------------------------------------------------------------- +# Reference particle +# -------------------------------------------------------------------------------------- + + +@dataclass(frozen=True) +class ImpactXRefPart: + """An ImpactX reference particle, detached from any live ``ImpactX`` session. + + Holding this as a plain dataclass rather than wrapping ``impactx.RefPart`` is what + lets the converters run with no ImpactX object in the process -- which the openPMD + reader and the whole test suite rely on. + + Attributes + ---------- + x, y, z : float + Lab-frame position of the reference particle, in metres. + t : float + ``c * t`` of the reference particle, in **metres** (ImpactX convention). + px, py, pz : float + Lab-frame momenta normalized by ``m * c``, i.e. ``beta_i * gamma``. + Dimensionless. + pt : float + ``-gamma`` of the reference particle. Dimensionless. + mass_MeV : float + Rest mass in MeV. + charge_qe : float + Charge in units of the elementary charge, e.g. -1 for an electron. + s : float + Integrated path length along the reference orbit, in metres. + gyromagnetic_anomaly : float + Anomalous magnetic moment, dimensionless. Carried for round-tripping only; + `ParticleGroup` has no spin. + """ + + x: float + y: float + z: float + t: float + px: float + py: float + pz: float + pt: float + mass_MeV: float + charge_qe: float + s: float = 0.0 + gyromagnetic_anomaly: float = 0.0 + + @property + def mass_eV(self) -> float: + """Rest mass in eV.""" + return self.mass_MeV * 1.0e6 + + @property + def gamma(self) -> float: + """Relativistic gamma of the reference particle.""" + return -self.pt + + @property + def beta_gamma(self) -> float: + """Magnitude of the normalized reference momentum.""" + return float(np.sqrt(self.px**2 + self.py**2 + self.pz**2)) + + @property + def qm_eV(self) -> float: + """Charge over mass in 1/eV, the form ``add_n_particles`` expects.""" + return self.charge_qe / self.mass_eV + + @property + def qm_SI(self) -> float: + """Charge over mass in C/kg, the form ``qm`` is written with in openPMD output.""" + return self.qm_eV * c_light**2 + + +def _reference_is_physical(ref: ImpactXRefPart) -> bool: + """True when ``ref`` has a positive mass and is actually moving.""" + return bool(ref.mass_MeV > 0.0 and ref.gamma >= 1.0 and ref.beta_gamma > 0.0) + + +def _check_reference_particle(ref: ImpactXRefPart) -> None: + """Refuse a reference particle the conversion cannot use. + + A default-constructed ``RefPart`` has zero mass and zero energy, which makes + ``beta_gamma`` zero and every converted momentum infinite. ImpactX writes exactly + that into its ``particles_lost`` output, so this is not a hypothetical. + + Parameters + ---------- + ref : ImpactXRefPart + The reference particle to check. + + Raises + ------ + ValueError + If the mass is not positive or the particle is not moving. + """ + if _reference_is_physical(ref): + return + raise ValueError( + f"The reference particle is not physical: mass_MeV={ref.mass_MeV}, " + f"gamma={ref.gamma}, beta_gamma={ref.beta_gamma}. A mass and an energy are " + "both required to convert ImpactX' normalized coordinates." + ) + + +def pmd_species_of(ref: ImpactXRefPart, rtol: float = 1e-6) -> str: + """Infer the openPMD-beamphysics species name from a reference particle. + + Parameters + ---------- + ref : ImpactXRefPart + The reference particle. + rtol : float + Relative tolerance for the mass match. ImpactX and openPMD-beamphysics carry + electron masses that differ in the 9th digit, so this cannot be exact. + + Returns + ------- + str + A species name such as ``"electron"``. + + Raises + ------ + ValueError + If no known species matches; pass ``species=`` explicitly in that case. + """ + for pmd_name in IMPACTX_TO_PMD_SPECIES.values(): + charge_matches = np.isclose( + ref.charge_qe, charge_of(pmd_name) / e_charge, rtol=rtol + ) + mass_matches = np.isclose(ref.mass_eV, mass_of(pmd_name), rtol=rtol) + if charge_matches and mass_matches: + return pmd_name + raise ValueError( + f"Cannot infer a species from charge_qe={ref.charge_qe} and " + f"mass_MeV={ref.mass_MeV}. Pass species= explicitly." + ) + + +def _check_species_matches_reference( + species: str, ref: ImpactXRefPart, rtol: float = 1e-6 +) -> None: + """Refuse a species that is not the one the reference particle describes. + + The momenta are un-normalized with the *reference* mass, and ``ParticleGroup`` then + reads them back with the *species* mass. If the two disagree nothing raises on its + own: the bunch is silently relabelled and comes back at the wrong energy -- a 2 GeV + electron beam read as protons reports gamma = 2.35. + + Parameters + ---------- + species : str + openPMD-beamphysics species name. + ref : ImpactXRefPart + The reference particle the coordinates are relative to. + rtol : float + Relative tolerance. ImpactX and openPMD-beamphysics carry electron masses that + differ in the 9th digit, so this cannot be exact. + + Raises + ------ + ValueError + If the species' rest mass or charge does not match the reference particle's. + """ + # named for the reference particle, not just "matches": lume-impactx, which shares + # this module verbatim, has its own _check_species_matches(pg, ref) for injection + mass_matches = np.isclose(mass_of(species), ref.mass_eV, rtol=rtol) + charge_matches = np.isclose(charge_of(species) / e_charge, ref.charge_qe, rtol=rtol) + if mass_matches and charge_matches: + return + raise ValueError( + f"species={species!r} (mass {mass_of(species):.6e} eV, charge " + f"{charge_of(species) / e_charge:+.3f} e) is not the species the reference " + f"particle describes (mass {ref.mass_eV:.6e} eV, charge {ref.charge_qe:+.3f} " + "e). The momenta are normalized by the reference mass, so relabelling alone " + "would return the bunch at the wrong energy. Pass species= only to name a " + "species that pmd_species_of() cannot infer, not to convert between species." + ) + + +# -------------------------------------------------------------------------------------- +# Coordinate conversion +# -------------------------------------------------------------------------------------- + + +def particlegroup_to_impactx(pg: ParticleGroup, ref: ImpactXRefPart) -> dict: + """Convert a ``ParticleGroup`` to ImpactX fixed-s beam arrays. + + Parameters + ---------- + pg : ParticleGroup + The bunch to convert. Unless it already sits exactly on one plane it is copied + and drifted to its own mean ``z``, so the input is never mutated. + + ``pg.x`` and ``pg.y`` are taken to be *already relative to the reference + particle*, which is the frame the reader returns and the frame ImpactX works + in; ``ref.x`` and ``ref.y`` are not subtracted. A lab-frame bunch around a + reference orbit that is off-axis must be shifted first. ``pg.z`` does not enter + the result at all: in the local frame the plane *is* the reference particle's + location. + + ``pg.status`` and ``pg.id`` are not carried -- ImpactX has no equivalent of the + first, and ``add_n_particles`` assigns its own ids -- so filter dead particles + beforehand if that matters. + ref : ImpactXRefPart + The reference particle the ImpactX coordinates are relative to. + + Returns + ------- + dict + Keys ``position_x``, ``position_y``, ``position_t``, ``momentum_x``, + ``momentum_y``, ``momentum_t`` (arrays), ``weighting`` (array, real particles + per macroparticle), ``qm`` (scalar, 1/eV) and ``species`` (str) -- the names + ``ImpactXParticleContainer.to_df()`` uses, which is also what + :func:`impactx_to_particlegroup_data` reads back. + + ``ImpactXParticleContainer.add_n_particles`` takes the same quantities under + shorter names and has no species argument, so feed it as + ``add_n_particles(x=data["position_x"], ..., px=data["momentum_x"], ..., + qm=data["qm"], w=data["weighting"])``. + """ + # The bunch must occupy a single plane. A t-coordinate bunch (spread in z) is + # drifted to its own mean z on a copy, so the input is never mutated. + _check_reference_particle(ref) + + if not pg.in_z_coordinates: + pg = pg.copy() + pg.drift_to_z() + + mass_eV = ref.mass_eV + beta_gamma = ref.beta_gamma + + position_x = pg.x + position_y = pg.y + position_t = c_light * pg.t - ref.t + + momentum_x = pg.px / mass_eV / beta_gamma + momentum_y = pg.py / mass_eV / beta_gamma + + # gamma - gamma_ref through the algebraic identity + # (gamma^2 - gamma_ref^2) / (gamma + gamma_ref), which avoids subtracting two + # numbers that are both ~4000 for a 2 GeV beam and differ by ~1e-3. The float64 + # representation of pg.p already sets the accuracy floor, so the gain over the + # plain difference is a small constant factor rather than orders of magnitude -- + # but it is free. + p_mc2 = (pg.p / mass_eV) ** 2 + gamma = np.sqrt(1.0 + p_mc2) + dgamma = (p_mc2 - beta_gamma**2) / (gamma + ref.gamma) + momentum_t = -dgamma / beta_gamma + + return { + "position_x": position_x, + "position_y": position_y, + "position_t": position_t, + "momentum_x": momentum_x, + "momentum_y": momentum_y, + "momentum_t": momentum_t, + "weighting": pg.weight / abs(charge_of(pg.species)), + "qm": ref.qm_eV, + "species": pg.species, + } + + +def impactx_to_particlegroup_data( + data: dict, + ref: ImpactXRefPart, + species: str | None = None, +) -> dict: + """Convert ImpactX fixed-s beam arrays to ``ParticleGroup`` data. + + The inverse of :func:`particlegroup_to_impactx`. + + Parameters + ---------- + data : dict + Arrays keyed as ``ImpactXParticleContainer.to_df()`` names them: + ``position_x/y/t``, ``momentum_x/y/t``, ``weighting``. Optional ``id`` and + ``status`` keys are carried through. + ref : ImpactXRefPart + The reference particle the ImpactX coordinates are relative to. + species : str, optional + openPMD-beamphysics species name. Inferred from ``ref`` when omitted. + + Returns + ------- + dict + Suitable for ``ParticleGroup(data=...)``: ``x``, ``y``, ``z`` in metres, + ``px``, ``py``, ``pz`` in eV/c, ``t`` in seconds, ``weight`` in Coulomb, + ``status`` and ``species``. The result is in z-coordinates: every ``z`` is + zero, the reference plane, and the bunch length shows up as a spread in ``t``. + + Raises + ------ + ValueError + If a particle's transverse momentum exceeds its total momentum, which makes + ``pz`` imaginary. That means ``data`` and ``ref`` do not belong together. + """ + _check_reference_particle(ref) + + if species is None: + species = pmd_species_of(ref) + else: + _check_species_matches_reference(species, ref) + + mass_eV = ref.mass_eV + beta_gamma = ref.beta_gamma + n = len(np.asarray(data["position_x"])) + + # In the local frame the reference particle is the origin: its transverse momentum + # is zero by construction and its longitudinal momentum is |p_ref|. ref.px / ref.pz + # are *lab* components and must not be mixed in here. + px_mc = beta_gamma * np.asarray(data["momentum_x"], dtype=float) + py_mc = beta_gamma * np.asarray(data["momentum_y"], dtype=float) + # ref.gamma is -ref.pt, so this is gamma_ref + (gamma - gamma_ref) = gamma. + gamma = ref.gamma - beta_gamma * np.asarray(data["momentum_t"], dtype=float) + + pz_mc2 = gamma**2 - 1.0 - px_mc**2 - py_mc**2 + n_bad = int(np.count_nonzero(pz_mc2 < 0.0)) + if n_bad: + raise ValueError( + f"{n_bad} of {n} particles have a transverse momentum larger than their " + "total momentum, so pz would be imaginary. The reference particle most " + "likely does not belong to this bunch: check mass_MeV, pt and the " + "momentum normalization." + ) + pz_mc = np.sqrt(pz_mc2) + + if "status" in data: + status = np.asarray(data["status"]) + else: + status = np.full(n, int(ParticleStatus.ALIVE)) + + pg_data = { + "x": np.asarray(data["position_x"], dtype=float), + "y": np.asarray(data["position_y"], dtype=float), + "z": np.zeros(n), # Zero by definition in z-coordinates + "px": px_mc * mass_eV, + "py": py_mc * mass_eV, + "pz": pz_mc * mass_eV, + "t": (ref.t + np.asarray(data["position_t"], dtype=float)) / c_light, + "weight": np.asarray(data["weighting"], dtype=float) * abs(charge_of(species)), + "status": status, + "species": species, + } + if "id" in data: + pg_data["id"] = np.asarray(data["id"]) + return pg_data + + +# -------------------------------------------------------------------------------------- +# Per-particle data ParticleGroup cannot hold +# -------------------------------------------------------------------------------------- + +#: The spin components, which ImpactX always allocates and always writes. They stay at +#: exactly zero unless the beam was seeded with a spin distribution -- ``sim.spin = +#: True`` alone is not enough, the gate is the ``spin_distr`` argument to +#: ``add_particles``. So testing for "any non-zero" is exact: zero means there is +#: genuinely nothing to lose. +SPIN_COLUMNS = ("spin_x", "spin_y", "spin_z") + + +class UnrepresentableParticleData(NotImplementedError): + """Raised when a bunch carries per-particle data ``ParticleGroup`` cannot hold. + + Converting anyway would return a bunch that looks right and has silently lost + physics, so the conversion refuses instead. Pass ``strict=False`` to + :func:`read_beam_monitor` to drop the extra data deliberately. + """ + + +def _unrepresentable_in(columns: dict) -> list[str]: + """Name the per-particle data in ``columns`` that ``ParticleGroup`` cannot hold. + + Parameters + ---------- + columns : dict + Extra per-particle arrays keyed by ImpactX SoA name, i.e. everything beyond + the coordinates, weighting and id that this module maps. + + Returns + ------- + list of str + Human-readable descriptions, empty when the bunch converts losslessly. + """ + carried = [] + if any( + name in columns and np.any(np.asarray(columns[name]) != 0.0) + for name in SPIN_COLUMNS + ): + carried.append("spin (spin_x/y/z)") + runtime = sorted(name for name in columns if name not in SPIN_COLUMNS) + if runtime: + carried.append(f"runtime components {runtime}") + return carried + + +def _check_representable(columns: dict) -> None: + """Refuse to convert a bunch whose extra per-particle data would be dropped. + + Parameters + ---------- + columns : dict + Extra per-particle arrays keyed by ImpactX SoA name. + + Raises + ------ + UnrepresentableParticleData + If any spin component is non-zero, or any runtime component is present. + """ + carried = _unrepresentable_in(columns) + if not carried: + return + raise UnrepresentableParticleData( + f"This bunch carries {' and '.join(carried)}, which ParticleGroup cannot " + "represent. Converting would silently drop it. Pass strict=False to drop it " + "deliberately, or work with the ImpactX particle container directly. Note " + "that ImpactX's particles_lost output always carries a runtime 's_lost' " + "component." + ) + + +# -------------------------------------------------------------------------------------- +# openPMD BeamMonitor reader +# +# ImpactX's elements.BeamMonitor writes standard openPMD with the species always named +# "beam" -- also in the particles_lost output, where only the *file* is named +# differently. It records position/{x,y,t}, momentum/{x,y,t}, positionOffset/{x,y,t}, +# weighting, qm, spin/{x,y,z} and id, plus the reference particle and the reduced beam +# characteristics as per-iteration species attributes. Verified against ImpactX 26.08. +# +# This path needs no live ImpactX object, which is what makes ImpactXRefPart a plain +# dataclass rather than a wrapper around impactx.RefPart. +# -------------------------------------------------------------------------------------- + +#: openPMD records this reader consumes or deliberately ignores. Everything else in a +#: species is per-particle data ParticleGroup has no place for -- ImpactX's ``spin``, +#: or a runtime component such as the ``s/lost`` of the particles_lost output. +#: +#: ``positionOffset`` holds ``(x_ref, y_ref, t_ref)``, not zeros: the longitudinal part +#: is applied via the ``t_ref`` attribute, the transverse part is deliberately not, see +#: the module docstring. ``qm`` is redundant with the reference particle's mass and +#: charge, which are what this reader uses. +_CONSUMED_RECORDS = frozenset( + { + "position", + "positionOffset", + "momentum", + "weighting", + "qm", + "id", + } +) + +#: ImpactX writes this zero-extent placeholder instead of particle records when a +#: BeamMonitor is configured with ``particles=False`` (moments-only output). +_EMPTY_PLACEHOLDER_RECORD = "empty" + +#: The runtime component ImpactX attaches to every particle in its ``particles_lost`` +#: output: the path length ``s`` at which the particle was lost. +_S_LOST_COLUMN = "s_lost" + +#: ``status`` for a particle ImpactX has lost. +#: +#: `ParticleStatus` defines only ``CATHODE = 0`` and ``ALIVE = 1``; `ParticleGroup` +#: counts ``status == 1`` as alive and everything else as dead, and each interface +#: passes its own source code's value straight through -- `beamphysics.interfaces.bmad` +#: hands Bmad's ``state`` over as ``status`` verbatim, loss codes and all. There is no +#: universal "lost" value to reach for, so this is a choice. +#: +#: The one value it must not be is ``0``: that is a positive claim that the particle is +#: sitting at the source, and `beamphysics.interfaces.astra` writes ``status == 0`` back +#: out as Astra's ``-1``, "at the cathode". ``2`` is outside Bmad's loss-direction range +#: and carries no such meaning -- it reads as simply "not alive". +PARTICLE_STATUS_LOST = 2 + +# AMReX packs a particle's identity into one uint64 ``idcpu``: bit 63 marks the +# particle valid, bits 24-62 hold the id, bits 0-23 the originating MPI rank. The id +# counter is per-rank, so the id on its own repeats across the ranks of a parallel run +# -- only the packed value is globally unique, and that is what becomes the +# ``ParticleGroup`` id. +_AMREX_VALID_BIT = np.uint64(1) << np.uint64(63) + + +def particle_id_from_idcpu(idcpu) -> tuple[np.ndarray, np.ndarray]: + """Split AMReX's packed ``idcpu`` into a ``ParticleGroup`` id and a validity flag. + + ImpactX writes the raw ``idcpu`` as the openPMD ``id`` record. Only the validity + bit is stripped off here, for two reasons: with it set the value exceeds the range + of a signed 64-bit integer, which is what ``ParticleGroup`` stores ids in, and + aliveness belongs in ``status`` rather than in the id. Everything identifying the + particle -- AMReX' per-rank id *and* the rank it came from -- is kept, so the + result is unique across a parallel run where the id alone would not be. + + The original value is recovered as ``np.uint64(id) | (np.uint64(1) << 63)`` for a + live particle; AMReX' own id and rank are ``id >> 24`` and ``id & 0xFFFFFF``. + + Parameters + ---------- + idcpu : array_like of uint64 + The ``id`` record as stored. + + Returns + ------- + ids : np.ndarray of int64 + ``idcpu`` with the validity bit cleared. + valid : np.ndarray of bool + True where the particle is marked valid. + """ + idcpu = np.asarray(idcpu, dtype=np.uint64) + valid = (idcpu & _AMREX_VALID_BIT).astype(bool) + return (idcpu & ~_AMREX_VALID_BIT).astype(np.int64), valid + + +def refpart_from_openpmd(species: Any) -> ImpactXRefPart: + """Rebuild a reference particle from a BeamMonitor species' attributes. + + Parameters + ---------- + species : openpmd_api.ParticleSpecies + A species from a BeamMonitor iteration. + + Returns + ------- + ImpactXRefPart + ``mass_ref`` is stored in kg and ``charge_ref`` in Coulomb, so both are + converted here to the MeV / elementary-charge units the converters use. + + Raises + ------ + KeyError + If the species carries no reference particle attributes, i.e. it was not + written by an ImpactX BeamMonitor. + """ + attributes = set(species.attributes) + required = { + "x_ref", "y_ref", "z_ref", "t_ref", + "px_ref", "py_ref", "pz_ref", "pt_ref", + "mass_ref", "charge_ref", "s_ref", + } # fmt: skip + missing = sorted(required - attributes) + if missing: + raise KeyError( + f"Species is missing the reference particle attributes {missing}. " + "This does not look like ImpactX BeamMonitor output." + ) + + get = species.get_attribute + return ImpactXRefPart( + x=get("x_ref"), + y=get("y_ref"), + z=get("z_ref"), + t=get("t_ref"), + px=get("px_ref"), + py=get("py_ref"), + pz=get("pz_ref"), + pt=get("pt_ref"), + mass_MeV=get("mass_ref") * c_light**2 / e_charge / 1.0e6, + charge_qe=get("charge_ref") / e_charge, + s=get("s_ref"), + gyromagnetic_anomaly=( + get("gyromagnetic_anomaly_ref") + if "gyromagnetic_anomaly_ref" in attributes + else 0.0 + ), + ) + + +def _import_openpmd_api(): + """Import ``openpmd_api``, with an actionable message when it is missing.""" + try: + import openpmd_api + except ImportError as exc: # pragma: no cover - depends on the install + raise ImportError( + "Reading ImpactX BeamMonitor output needs the openpmd-api Python package: " + "conda install -c conda-forge openpmd-api, or pip install openpmd-api." + ) from exc + return openpmd_api + + +def read_beam_monitor_data( + path: str | pathlib.Path, + iteration: int | None = None, + species_name: str = "beam", + species: str | None = None, + strict: bool = True, + ref: ImpactXRefPart | None = None, +) -> dict: + """Read an ImpactX ``BeamMonitor`` openPMD file into ``ParticleGroup`` data. + + See :func:`read_beam_monitor` for the parameters; this is the same reader without + the final ``ParticleGroup`` construction. + + Returns + ------- + dict + Suitable for ``ParticleGroup(data=...)``. + """ + io = _import_openpmd_api() + + series = io.Series(str(path), io.Access.read_only) + try: + iterations = list(series.iterations) + if not iterations: + raise KeyError(f"No iterations in {str(path)!r}.") + if iteration is None: + iteration = iterations[-1] + elif iteration not in iterations: + raise KeyError( + f"Iteration {iteration} not in {str(path)!r}; have {iterations}." + ) + + particles = series.iterations[iteration].particles + if species_name not in particles: + raise KeyError( + f"Species {species_name!r} not in {str(path)!r}; have " + f"{list(particles)}. ImpactX always names it 'beam', including in " + "particles_lost output." + ) + beam = particles[species_name] + + records = list(beam) + if _EMPTY_PLACEHOLDER_RECORD in records: + raise KeyError( + f"{str(path)!r} iteration {iteration} holds no particles: this " + "BeamMonitor was configured with particles=False and wrote only the " + "reduced beam characteristics, which are on the species' attributes." + ) + + def _load(record_name: str, component: str): + record_component = beam[record_name][component] + return record_component, record_component.load_chunk() + + components = { + "position_x": _load("position", "x"), + "position_y": _load("position", "y"), + "position_t": _load("position", "t"), + "momentum_x": _load("momentum", "x"), + "momentum_y": _load("momentum", "y"), + "momentum_t": _load("momentum", "t"), + "weighting": _load("weighting", io.Record_Component.SCALAR), + } + idcpu = beam["id"][io.Record_Component.SCALAR].load_chunk() + + # ParticleGroup cannot hold spin or runtime components, so collect them -- + # either to refuse loudly, or to report what strict=False threw away. + extras = {} + for name in records: + if name in _CONSUMED_RECORDS: + continue + for component, record_component in beam[name].items(): + key = ( + name + if component == io.Record_Component.SCALAR + else f"{name}_{component}" + ) + extras[key] = record_component.load_chunk() + + file_ref = refpart_from_openpmd(beam) + series.flush() + + # openPMD stores values that must be multiplied by unitSI to reach the unit + # the record declares, and ImpactX writes 1.0 throughout. The positions are + # genuinely lengths, so the factor is applied. The momenta and the weighting + # are ImpactX' own dimensionless quantities -- normalized by the reference + # momentum, and a count of real particles -- so a factor other than 1 would + # mean the file no longer follows the convention decoded below. Refuse rather + # than scale them into silent nonsense. + data = {} + for key, (record_component, chunk) in components.items(): + unit_si = record_component.unit_SI + if key.startswith("position_"): + data[key] = np.asarray(chunk) * unit_si + elif unit_si != 1.0: + raise ValueError( + f"{key} in {str(path)!r} has unitSI={unit_si}, but this reader " + "decodes ImpactX' dimensionless convention, in which the momenta " + "are normalized by the reference momentum and the weighting " + "counts real particles." + ) + else: + data[key] = np.asarray(chunk) + extras = {key: np.asarray(chunk) for key, chunk in extras.items()} + ids, valid = particle_id_from_idcpu(idcpu) + finally: + series.close() + + if ref is None: + ref = file_ref + if not _reference_is_physical(ref): + raise ValueError( + f"{str(path)!r} iteration {iteration} carries a zeroed reference " + f"particle (mass_ref={ref.mass_MeV} MeV, gamma_ref={ref.gamma}), so " + "its normalized coordinates cannot be converted. ImpactX writes this " + "for its particles_lost output. Pass ref= with the ImpactXRefPart of " + "the BeamMonitor iteration the particles were lost at, which " + "refpart_from_openpmd() reads from the monitor file." + ) + + if strict: + _check_representable(extras) + else: + dropped = _unrepresentable_in(extras) + if dropped: + warnings.warn( + f"Dropping {' and '.join(dropped)} from {str(path)!r} iteration " + f"{iteration}: ParticleGroup cannot represent it.", + stacklevel=2, + ) + + data["id"] = ids + # AMReX' validity bit is not aliveness: ImpactX copies lost particles into a + # separate container and marks them valid *there*, so the bit is True for every + # particle in a particles_lost file and taking it at face value would report a + # bunch that is entirely alive. + # + # The s_lost runtime component identifies such a file: CollectLost always adds it. + # The zeroed reference particle is a second signature, but only because ImpactX + # fails to set one -- keyed on that alone, this would silently go back to reporting + # lost particles as alive the day ImpactX fixes it. + from_lost_file = _S_LOST_COLUMN in extras or not _reference_is_physical(file_ref) + alive = valid & (not from_lost_file) + data["status"] = np.where(alive, int(ParticleStatus.ALIVE), PARTICLE_STATUS_LOST) + + return impactx_to_particlegroup_data(data, ref, species=species) + + +def read_beam_monitor( + path: str | pathlib.Path, + iteration: int | None = None, + species_name: str = "beam", + species: str | None = None, + strict: bool = True, + ref: ImpactXRefPart | None = None, +) -> ParticleGroup: + """Read an ImpactX ``BeamMonitor`` openPMD file into a ``ParticleGroup``. + + Parameters + ---------- + path : str or pathlib.Path + Path to the file ImpactX wrote, e.g. ``diags/openPMD/monitor.h5``. Any backend + openpmd-api supports works, including ``.bp`` and a ``%T``-templated file-based + series. + iteration : int, optional + Which iteration to read. The last one when omitted. + species_name : str + openPMD species to read. ImpactX always writes ``"beam"`` -- its lost-particle + output differs in the *file* name (``particles_lost.*``), not the species name. + species : str, optional + openPMD-beamphysics species name; inferred from the reference particle when + omitted. + strict : bool + When True (the default), refuse to read a bunch that carries per-particle data + `ParticleGroup` cannot hold -- non-zero spin, or a runtime component such as + the ``s_lost`` that ImpactX's ``particles_lost`` output always carries. Set it + to False to drop that data and read the bunch anyway. + ref : ImpactXRefPart, optional + Reference particle to interpret the coordinates against. Taken from the file + when omitted, which is what you want for a BeamMonitor. ImpactX's + ``particles_lost`` output carries a *zeroed* reference particle, so reading it + requires passing one, e.g. + ``refpart_from_openpmd(series.iterations[n].particles["beam"])`` from the + monitor file. + + Be aware that this is an approximation for lost particles: they were lost at + whatever ``s`` the file's own ``s_lost`` record says, not at the monitor's, so + they are un-normalized with the wrong reference momentum and given the wrong + reference time unless the two happen to coincide. It is exact only for + particles lost at the monitor. + + Returns + ------- + ParticleGroup + In z-coordinates: every ``z`` is zero, the bunch length is a spread in ``t``, + and the transverse coordinates are relative to the reference particle. See the + module docstring for the frame conventions. + + ``status`` follows openPMD-beamphysics, where ``1`` is alive and anything else + is not -- ``pg.n_alive`` and ``pg.n_dead`` split on exactly that. Every particle + in a ``particles_lost`` file is marked :data:`PARTICLE_STATUS_LOST`, because + they are all lost by construction. + + Raises + ------ + KeyError + If the requested iteration or species is not in the file, if the species is + not ImpactX BeamMonitor output, or if the monitor recorded no particles. + UnrepresentableParticleData + If ``strict`` and the monitor recorded spin or runtime components. + ValueError + If neither ``ref`` nor the file provides a usable reference particle, or if a + record carries a ``unitSI`` this reader cannot honour. + ImportError + If the ``openpmd-api`` package is not installed. + + Examples + -------- + >>> from beamphysics.interfaces.impactx import read_beam_monitor + >>> P = read_beam_monitor("diags/openPMD/monitor.h5") # doctest: +SKIP + >>> P.norm_emit_x # doctest: +SKIP + """ + return ParticleGroup( + data=read_beam_monitor_data( + path, + iteration=iteration, + species_name=species_name, + species=species, + strict=strict, + ref=ref, + ) + ) + + +def beam_monitor_iterations(path: str | pathlib.Path) -> list[int]: + """List the iterations available in an ImpactX ``BeamMonitor`` file. + + Parameters + ---------- + path : str or pathlib.Path + Path to the file ImpactX wrote. + + Returns + ------- + list of int + The iteration numbers, in file order. + """ + io = _import_openpmd_api() + + series = io.Series(str(path), io.Access.read_only) + try: + return list(series.iterations) + finally: + series.close() diff --git a/beamphysics/particles.py b/beamphysics/particles.py index 64aba263..a6118b78 100644 --- a/beamphysics/particles.py +++ b/beamphysics/particles.py @@ -3,7 +3,7 @@ import os import pathlib from copy import deepcopy -from typing import Union, Optional, Sequence +from typing import TYPE_CHECKING, Union, Optional, Sequence import numpy as np from h5py import File, Group @@ -45,6 +45,9 @@ from .wakefields import WakefieldBase from .writers import pmd_init, write_pmd_bunch +if TYPE_CHECKING: + from .interfaces.impactx import ImpactXRefPart + __all__ = [ "load_bunch_data", # Re-exported for backwards compatibility ] @@ -1073,6 +1076,29 @@ def to_bmad(self, p0c=None, tref=None): """ return bmad.particlegroup_to_bmad(self, p0c=p0c, tref=tref) + def to_impactx(self, ref: ImpactXRefPart) -> dict: + """ + Convert to ImpactX fixed-s beam arrays. + + See `beamphysics.interfaces.impactx` for the coordinate conventions; + note in particular that `x` and `y` are taken to be relative to the + reference particle already. + + Parameters + ---------- + ref : beamphysics.interfaces.impactx.ImpactXRefPart + The reference particle the ImpactX coordinates are relative to. + + Returns + ------- + dict + Arrays keyed as `ImpactXParticleContainer.to_df()` names them. + """ + # imported here: interfaces.impactx imports ParticleGroup at module level + from .interfaces.impactx import particlegroup_to_impactx + + return particlegroup_to_impactx(self, ref) + @classmethod def from_bmad(cls, bmad_dict): """ @@ -1124,6 +1150,65 @@ def from_hdf5( data = load_species_data(group, include_time_offset=include_time_offset) return cls(data=data) + @classmethod + def from_impactx( + cls, + path: str | pathlib.Path, + iteration: int | None = None, + species_name: str = "beam", + species: str | None = None, + strict: bool = True, + ref: ImpactXRefPart | None = None, + ) -> ParticleGroup: + """ + Load an ImpactX `BeamMonitor` openPMD file. + + ImpactX is an s-based code, so the result is in z-coordinates: every `z` + is zero, the bunch length is a spread in `t`, and the transverse + coordinates are relative to the reference particle. See + `beamphysics.interfaces.impactx` for the frame conventions and for the + reference particle, which this does not return. + + Parameters + ---------- + path : str or pathlib.Path + Path to the file ImpactX wrote, e.g. `diags/openPMD/monitor.h5`. + iteration : int, optional + Which iteration to read. The last one when omitted. + species_name : str, optional + openPMD species to read. ImpactX always writes `"beam"`. + species : str, optional + openPMD-beamphysics species name; inferred from the reference + particle when omitted. + strict : bool, optional + Refuse to read a bunch carrying per-particle data ParticleGroup + cannot hold (non-zero spin, or a runtime component such as the + `s_lost` of ImpactX's `particles_lost` output). Default is True. + ref : beamphysics.interfaces.impactx.ImpactXRefPart, optional + Reference particle to interpret the coordinates against. Taken from + the file when omitted. ImpactX's `particles_lost` output carries a + zeroed reference particle, so reading it needs one; see + `beamphysics.interfaces.impactx.read_beam_monitor`, which this + delegates to, for the caveats. + + Returns + ------- + ParticleGroup + """ + # imported here: interfaces.impactx imports ParticleGroup at module level + from .interfaces.impactx import read_beam_monitor_data + + return cls( + data=read_beam_monitor_data( + path, + iteration=iteration, + species_name=species_name, + species=species, + strict=strict, + ref=ref, + ) + ) + @classmethod def from_genesis4( cls, diff --git a/docs/examples/data/impactx/generate.py b/docs/examples/data/impactx/generate.py new file mode 100755 index 00000000..05664b16 --- /dev/null +++ b/docs/examples/data/impactx/generate.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 +"""Regenerate the ImpactX BeamMonitor test data in this directory. + +Run with an environment that has ImpactX installed, e.g.:: + + conda install -c conda-forge impactx + python generate.py + +It writes ``monitor.h5`` (a two-iteration beam monitor) and +``particles_lost.h5`` (the lost-particle output, which carries the runtime +``s_lost`` component). Both are read by ``tests/test_impactx.py``. +""" + +import os +import pathlib +import shutil +import tempfile + +from impactx import ImpactX, distribution, elements + +HERE = pathlib.Path(__file__).resolve().parent + +# ImpactX writes its diags/ tree into the working directory, so run somewhere +# disposable rather than leaving it next to the data files. +_cwd = pathlib.Path.cwd() +_tmp = tempfile.TemporaryDirectory() +os.chdir(_tmp.name) + +sim = ImpactX() +sim.space_charge = False +sim.slice_step_diagnostics = False +sim.particle_lost_diagnostics_backend = "h5" +sim.init_grids() + +# 2 GeV electron beam, as in the ImpactX FODO example +sim.beam.ref.set_species("electron").set_kin_energy_MeV(2.0e3) +distr = distribution.Waterbag( + lambdaX=3.9984884770e-5, + lambdaY=3.9984884770e-5, + lambdaT=1.0e-3, + lambdaPx=2.6623538760e-5, + lambdaPy=2.6623538760e-5, + lambdaPt=2.0e-3, + muxpx=-0.846574929020762, + muypy=0.846574929020762, +) +sim.add_particles(1.0e-9, distr, 200) + +monitor = elements.BeamMonitor("monitor", backend="h5") +sim.lattice.extend( + [ + monitor, + elements.Drift(name="d1", ds=0.25, nslice=1), + # a collimator tight enough to scrape the tails, so that + # particles_lost.h5 has a few particles in it + elements.Aperture( + name="ap", + aperture_x=1.55e-4, + aperture_y=1.55e-4, + shape="rectangular", + action="transmit", + ), + elements.Quad(name="q1", ds=1.0, k=1.0, nslice=1), + monitor, + ] +) + +sim.track_particles() +sim.finalize() + +for name in ("monitor.h5", "particles_lost.h5"): + shutil.copy(pathlib.Path("diags") / "openPMD" / name, HERE / name) + print(f"wrote {HERE / name}") + +os.chdir(_cwd) +_tmp.cleanup() diff --git a/docs/examples/data/impactx/monitor.h5 b/docs/examples/data/impactx/monitor.h5 new file mode 100644 index 0000000000000000000000000000000000000000..4229a091beb52fcbca0d7782fcdf32817f7defed GIT binary patch literal 130410 zcmeFa2VB<5zV{n?Zx3C16$J&Qh=3VDs(|z&9i=KwX`&!tM^sQLc0s^~qM#_Y41!p& zfr6;0sGy>tQWcPL6Q1EL*1l_<{kqTH=e*wY44>rt%w%RV^UEapXUI&_t;{U=xTkS5 zb8=#nZXCJgqf8S)T4=*eyW^kit<0R5PcB$o(leIw~eCHY_S~of!(2q+;_=2HOV*#{|Wm zUWqo+^IjQI6NdumvcyD(Mea^PLG*){5EK>?8XFcFQgxahjtQZoh`bDt3f$_h@3wHV za=@r74l_o!yL^^4d%V6eP`K#mcF2pcvB@#IR$l*vSBtY1dZ{4(u%1bSGoDsomnhmT~^OgfE+o|8-(B{S)yg)eM1rT=8Ce>(-x`EN4iyeqlhQ;aE5 z#~U?>;HH`q>;~1WP~#PoO*@AfQ!Qj3y@fp18NP(YlUB@Rl88l*uV5Ba12SJ-#rmi` zP^OvL(N&nmlxFALT2hO>ozj~7Yt?UKo2b(aCAp!mUWZu-k^$&Mw9#H(kBL#4qVxYj zmd6XMVoH)Z#p*4X0(A{Rc?{jgUSfh%YI|3p_$Agzl>=qka%pl*c7M!uiikRIqn2WpR18$TviNPwFy|>PX$g(q!~~9z>A$u|Ia3mUAN# zn?+q_QAzRC(@|LV6q(I(XT)QRDG|}7%e-PN35%yRL&snF6jL%bcgoe>&R}{fradJO zzpF7@FtI5~2*ur7u*NA(rl{k`R_qv6Qgk39H*{>nHcUBo0jG>Ju%;<7%rEQgz*4E> zfRgxnH1EbPPbq29t|fU`@DvG!=M)~mET)u$-MRP>c4Nw2rpR&V5Juib|D$r&@$2(3 zmMN7hI4)m+J)uNK)&Fq#L7{ZqMJf}z`)eZ@V`3t09A+jaraLI!Bw2VdCMrLw)nS-1 z#zsDQv2bvb_AU!E^C6UPD!GY@@YfFyykuEeFdj^VvWjBEB7*1!FD4`}hwLPIOq>%F zPa$)l8Vu#A2aX4_^xqrk!07pK(DR9l42yNPCi4LQQ8XPMOlO=_crp(R>jsz#O8n1I zb}w_4i<7xIY1WhN_s{MRYwpjvJearykYQB2Ctd2q`j6GgtAG5aex)KWT}$L)fS~|G z0fqt$1sDo26!_f~pdQy#$%V2xC{hLCcr!)Su!W5EQ;hmZ-9tI^0vqWj>+-q9LLY;sQh*5g-~P znjo4YE<`j#G)J^Rv_!N*v_@QnXoF~rXot8M(H_wO(Gk%J(HYSN(G}4Rk?LGbB^M9G zC5WDgUWnd^K8U`Eeu(~v0fl=-453R3cIAn=bs9KE9cvnTCej2EtM-Da4+`&Y9sM1f;Mh2O51Mi^j#AW{h{a1 zKv6JoC=uynKiGp;d@s9Xqhbl!a$YMm?+gIiIjjp$J1>KWecg`iUM~oZyjKRxTwe&j z=~}Jw(OU%yUKVQV4Fy8X?DEqcQqE9kb&=~URROTC6P}KX7KEt}Z46|$2ZQdDPh|rh zJHyK%Y&6a39#HTs=g5OoO2BkO;Io@LK_Gprj?TEABXnjrPOjMK4-}@qxs>{W zr$vVkHj7%rn)X5Y^pj!0eVO5=JueqRGuuyn1uufZ!*>opB5ni%funDY6f#}mjp=I6 zBB#OtxRz9SP}UAKTpe7Zy4?#}R7qCu+2{nWG&_bmFLeXwmOiZSSQQ9q18ZMo@3aT5 z&hbCO&zOMPFAZu-rZ0nYR^jJA7+b)x9mMr5O{O5Y*70$Itvfg@eb{i5xDR+TmKk_C zb17Iq`+T6Zuq3FMd3>WnqdL?rTd9uLO{)GI|FyN$a_Bs8P4v>DS z)zg$71oka!F3j`^26fDBX>6S@X0zgjj3eC?UJi@%K29xdzARf0plk{8$$eF2wf>)4R4E9l7Hqh$BW5{TVTvZ@wy2G^~& z8F_jIg6UqIRcmqsKxIx+uJEQX7_jM7DAWi7Y+I6q<^&Rmr#4SaD-#u(YgBY z;AwNO03&~JyuxU`%n?hdM7!xP-mn7B*Hp&le6#{*E_RQk`I|AXm2QK3l zfdPO!>a%~}3|Gj#?ESU;m7ZXO`A$>O@DSkOcxy4=*dqAA3X!Qa}f`?KvmB zz|9-@-tC=r;b$PwO6(}sZ&(UNzUww>#Vi1Jw#QWzdjY_H zp_L6Uc|~C%Z_gs+A4YXbn2{YEkul_d`~JFbg6v|UYyiqllJ5KU_25@y+kfTXtfXy0 z{xdqZjk=eU4*fHFwI+2h{==WRmbxdB?yoy_s1BOj)ZLDBm*^nZyyC_mN71vYr0AQ2$a75p|3D?{H_;RD*xHhiMur87rNs|4`7A@XrC0+bM|{ zj-dcU0fqt$1sDo26!<+A__x?nHdJQ|=~Bt%H?^f?o-bDTRaK`PLK?U)5qGK5Qw-!v-w|5$ ztpuMPx$JDoXFhmzv(2W4lplCrz8srk(c4Dqg_( z&Z3+{l?OV~u6$x4#9-T%0v+LbckpE8uHxf2xMA#x(wrSPnZRM+J&I=^i-7R%t~Q0! zvM@D&^LSw27d&dc_X(%nJRoD~x*?trR-iuHzR+Mu1aw||M0|DOfC?su59{u0!3~s! zoW9uzfxD|TRR{e=p~|XRcfD6WCCc|kloxPz5j$qf81dg7BJQfs_lrCx2ktImlJU>t z0zdNVEe_g$#S`?rlnv*J0=v<+9SSU5AUNtxKzNNH98=%Wp=~b?VuDNh9DerW6}IsC z`hc&5ep=I|H{wii!xp7)%y)UAlu+<47x^AM{c)`49z7;-+k3Nm^&W08^YYcGi%(b} zIH*4m$SwqXxIIRvy%z!>ogNOdWfT)K$A{g-+c|-5+uTcp3=2qWe4%pDSqzK?t(&JPNr_o&Z1!vod$J)N^OSV8$hrWFH&qeRU11-)8)y+lR*foy{u4zOb` z-}}|FJ;W*Yy8&+HgZR1k?cLt$(jbjL-HN~EBwot>uuwO19RH#HP_oNb2poB{igh>Z zN1``d>s|H@Vc^xC`X*|FAdt&dtc!c}o#=Fmy&V`iM$o=54)43sNE|%ds|ccc@Rz}T z*X?U~fOh;$pC#_>fQ6@AOF2{o#IBxg8oZ|gzpz?KBG}+N;h`IEDXY#4cH=9!x!2O* zWeZ13y|P^d)-L{_T7?xhaxU7=9zKSL>f{c}vx)!?X6?{D`WTdW9cQ{PrUI9=Wp1eT z7KMWxjW!o%3ByAwmkvE$aEBP(y;Oncn*c1C`z_ITP9M&4Hqkssg&BU9qm}B`4ifOm zi&%FDRv^c0d^0I~lnCkdt-5IW0%y{QB{GcJ!Q4->-RjXopoR6SP3Jc;Ab*a<;CZ$L z6f73FxzV?i020nmE+(=;M~<*HUTR(V@gchn5~ui~d%M%l3^O)hIb1!*NVJn!0&|Q; zz4{4LRZQ{Ts0jF?scW98#s_C)z8F$E{t|z4N0BLWuP_jQRr+Kdn*gwMs&sE&NdwvS zdsg0!lY+4(N3EDoio(oHZ+9*YQNR;tDg7kz2i_lGG-^5LJ>hoJwu)<$5De_A^5aT= zOAwvwzU9fh#uJWoKkF`Fh0c79f%o&?;F|9v684?!Ct8+le0uhj6r8crzyEkM2hjHF z3plh<5E$tWzgeQh3s}R9F1Me5jMHq_*U!4i2c>flIt#IO5}6k{T7ITA;d`DJGwnNw z!AI)0Nm8SIgm}i4=B@T@;GoxmlxUVu#N~NvMMLLAfZ^lEUw3i|g4GrymF?B{h|yue ztzv6M;DEfJq~2~0pk-a3DP185)L);uGChS8UV=YXr#>6NlL}&*9aqVN;k{~SOE=J< zQRn*e`yPoxt2-|TJb4Ad(N(FJTa%bz@MYFw7x}Mv$E#C3DQavmb%pNu24{B25qi22 zFX|wg#nMLhkFdj&9%c$B9Rxv|u2n0i4=0#k&-J9q@iQ@_YZ}j*W)@&@J~4H9vlL_z z7_gE$(@UJKTD$bL=?%ip*zyXiULVm_|90(3Q+}}If>fR8IVR{GD1gaTv=OPZ3g7v9 z&|qNBx8^#%0b+CfnZp*j-|;f@Ge7OT_~F}}@>cO54GfPMRa?y7++dSGIE!Z-qcWr^Kn} zPGk6E`|rnN^8`UBt$I!a7Y)d|%G@ZCYQuAu_-vdXA^~)aqo;MI58(&Gwm53rH4v$5 zF2q+y*c|Kx+Px|s-Wo{6G=69};2slb5K7_P{eT9en*(}+?0BG^ z0Z&*%ECzR!7_NzUao9VMN5;-KzzUgm-qFNs^qn{qDiZ6lJMbZjzT zOTcc0l2?L~`S^78$9uTGh`}#A_ExE7$-?dz^~r4mEWpS4)wbLS0VpSZ>7D2GH^c|y zVSz|rPMG+v*rUBn0D4<-R@{aSkD&-5`>%13^xg1Vg+S$d$c~J^%7mZQL{akGs7^I({;vP z9r)Z(8s9-yQCOWd&we4_cl_Vdmy$=lFhaUia`{dBQnHydoPW75g?gWca(~yCa+|(+ za-E8PF#1v^6=9H}07C(W0t^Kh3jBTwOuYF-zJvR#!@B+l-XR!R5A+fvs>zE=KQtkPN3>&B!@C70jSDt&(- z?C7srr9RYK(v0U6K%lxSS$EhQo=D`>yA$jTzU3ZFuVtSPY2f@fpy~|N&xqy=u5*P|_Mb;h)p5Yi zwL>azSOYkJv=)0%xCG=&9bfA1X$+g&W3N;FeYcBYH#_B`|yD2=j-%GrLm=0c@4~y(9_JOh%@)8-5PC#>ghdg$~ z1Jrpj(e;!1~HP^#*JSV zgKZZJmvfg{z|ZTSKe}+y1$=c%NZDjP2NZEFS9CuH;1bJ6M{X4wg7Z(){J0w|;I(z! zPOf9l;M}gPoo&im@aG%*)pdmqFha#zb0B^Z*x+`7U8BSR)W5%D!ky~?tj-FhMWs5x z(bXU0YF~Q8YZ|Xbu4aG79SW10oX^_ACGW+)@OaxmCOd4^v4!qXGSzzTz10g~+db>x zL{2u?;(I*w>strF6346Ob=MU#Em+XD@tGfZEw;t1<;irIv^J+KL(K!e8ozqJI@$sT zJh-xLca8zD?l^dSEJ+0hHQ9GuQnCXY`mcwLGF-p`_49+)1NPv{^gexq@x?G++*o1Z z8ynz1dQ;Zj$`!`D0f!y$P2rLg4TUo`MPNctrbfJ)2{`ieU z{mq;k=CI#7?eJ`OMW8VI4ENf(M!=>1OHF3me6Trv`lm1X#=z-Dm`VB?Ul?%y`9kcW z7R-p8q21wc0k%EcY&sMt3H(GC2~})#0ClEC;c+T*@T&ARXl*eYnwCd7c`BN~hKyBg zG;>$*zNL5WoGrFcKRtZoV1Pb+kZw@_@Q^7?bK-;Add5(mD|FP+X6PE z#~Co5>@5MuJFmjy>IZDLB1#o)#s z@gyc`5qKt+pQk=&1aEm%8-Cn<0AK9cay-&R323+#@wK*y0IiMl>~4v1f`{L?Op|St z23wT+rpXP8gRgfBv<5Z$i7TIDw>*mx0*lU>UT^s}iucwZD*Rc>4DMxoI~Mrl4X!M_ zf04!_CCIbD)wVkWgZIBE931Zw1-$jAnlIY7(u~L@VZ?@Aw zmyFTsdt0Z0#k$uPT)p2%@I}p|WpD@p=EHUCmguN}y!`#X--cK~pkqSNUJqqpQ`|Dr zm7xO8$upml{W?mtvL(DM+Q9_D5R3XEK4K z{C(H@56FZ4S3U;{UA#wJj(GXdE}0h=MlJH@KKT;o*V-_z+(`oRXQXcC^A-csyl>uW zy(J2FS=m-~l+u9V7ugrxxWU&Y13asvxmIxe zARMl>-Dy8Gj?;owAb2kW272{NvZ{o@<#p#YlFrJ3ZGJy)W`)W_9lMHVlNkf}=n|z_ zdrhSv&+Z);HfOLy7SWC7+AZQx=4D}M=M5T&Q`3HIVJrqWUGw`c+VBy#y4P{q~h=-#ocaML!sQDU*sY$WVZx07C(W0t^LyKL!3Rwp6V@U5Hc?zo{)H z8|u9IS8b_*PGHj_o>%#Cc6ksU8P=KKTLxJB(f!!M)eLGTpgV?tD zc=jMc_75JA5?p>x2!3-O|GIM}J6KtCGR-N-%#171n;`6ni)DS54hK_ zdj6qQ1O)iM$dIlQggd48=9Tky<2tDc^?9EKz>NJ40<6lSpsh}IdRB`#Xu;h*7M!jp zqEwYl7p{Ga8@%+*lWY}-zRH}Z%6X?j<8j^M^qJkbTp8Pq&9g^{@9TG6IVLX)4&3_O zyE=^(b|mi6S8Hd5yQ&{I+<(?h+}{83w#XO_F1Fv!{7qLF%vr|vu8fBdmWpz2ncgP{ zx=#xU7HYo3Qx}8samJ5oQ&v?h=HbmqPbn4q?NP?(!g`02I@xrFg*?SWA@xxvE`4Pj8 z%%Ez0mcBAMQpdV)-U-ehGSDpahh_Y-ZhWzvfB|QoJY+GC<} z@=O`tcMZPuba(%;tB@&WdP9mN@4EFdXb!q@DUFf@t0w1D}} zARc4JY-%5NgLrngZ5Jmd2@D(cn0GRF;4b>dPPYUIz=DAl%tCwk!6uf|I`%E1!0^+d z@5Tmw!~wnpi-R+SU~H7DlJ>MZ;n)`cQ!Bf?N6EByBH)+Ssd{k`pd*QbQW z*(dO5wGw!4B-nGA9I;7EFRE_jlm;JNR39ep-GwWR28>nRkO0RlWIWb<<_6p$Iq4yo zD4e(a%AEwecX;``97mbV9^6&><8^C)F?d~Lw&%C$b-2~j5V2)Jyx{&b_ddy8FY&z! z4bOepB>-Euyw#ZnSMY!j?=;F1&8it>=)hLmdFRoS@KQ5<#)Jk$2TjbvP4A>Jp7?kQH5)kPCV7YM z{g=L6lJ)#(7o)o52O9tXdvvMkWB(@yQT!J!3;(arX#OW(5c)rv-?f#Rpqq55SDe}~+>^xePQGfd5MPPxDPFMBh6^W-`e{b2MA zPb$J7Lji^Y3rehQ*Z0%?e{IZBh8fpsdW=7w z@(!yXUiEKz_rN5LzN=Ic|IEVt#cxdhcYQ&ifpcNeuUeQZ)9KazTMM%VZAVB&WLTJg z)`j6?D8NvFp#Vbxh5`%)7z!{HU?{*)fS~|G0fqwqH!JXO=}XzMn_lHq62EC*O0;ei z&oB3->_A8Q@BUlYinb%9A~O0?{;Uhb$54Qw07C(W0t^LyCk6g3w$$?jbO};P{HC^) z`5Di9ziLZy6`>>lceYd;+K!Nl$grjUtP8`(P=KKTLji^Y39T@H4f8TlVKRzb**DsR)@dcs(uLr-qJVP%`Dv5vA1;j%ivKku}5o8}6 z91|2vj$4&OFZar0;+&W`HP(%po*_oPW6gV!?)kS4e%p|hT=~mgKqpWs|E>$@BifFT zipb~!`m-(!A436#0t^Kh3NRG-ofP=LVelI@1%n1PuLld=t_q~=DKKG$VVE(d!-MIJ zv;1LCVOTfdolV7O`@;_VCn-n^OovJ_#d33v42yNP{wL|W0qcwNV}Csb$m;)W#uIjO zei9~(D>`6RY}m5!pqQIzblu+%!^=HFpRwfVWb+Z5{PB;Jp?O*VD@78Zip07!C?Miw zC(09jQ1c)6@XLN%@5`m8E>>pbK_ksxA#y)~YFoSM*)U`7D8JqUi3JwU&gPi=WqJsc z4>MLcQSFdLa+4o3uyC?+m^k>%7->0Ay2;67$P^Q9GKfNL`0KOt%NBGY$r z${01W^dj%+2Lr=0yStuFtr?A5a7dE+<1W)}77S3e-vB~mHrk~_wF<3DNY&E7v zdg)t|E9qOLP=rY)CtBV#b+NS|m3T~TBIn8(AIC7Q-XR&xiAt(U4@gloHpNy`0qXRC zCodlRV9Hbgd5%y?Di5fnV`RxNYOhp)Z!+6HY@&6MF(_vf-d2j`Q@`fYW5kUPR$>#= zVvuJImBckMrJ9^wgp8p?4=@pb@ouGH$WbxE+ytn=+MLb|PjXW+76eiBrjidfAlX9>i?K9K>A2JjA_-`w*#7 z$f@K)y=X!4AmSm!!-z)^^AQUWsjFftxfCK6Ar>PZLp+Xn0`Vl`DZ~;)YRq#gxlm92 z6iX4$BA!DmLo7!;k9YyG0+D*YPbC-XGEb2@t0_`f3W`?{sS}3sUqif(NFYMQD#ROz z)rd8SHxX+QZz0}Byn{%B;eJ;Ie7^dAs^}AjOq#(PfL18%%?}pKw3!bYYWiDbzdM0l zeY5Wej{1Q8VYC;f0zp7#13VOyXa?^;6!AMJxB||>&K7sBRfg@k90fKGPH?X+Q;(#x z58NANd(S72{AS*K*}aEl{6T`|>*DvdOX2x9WrCWnhV`c2lzU6wBzF3V4zr0ej_==5-wG;C|XkG3fIo&g6orf!1W`NExEpSAh1b&YwinY zc>R>3c=MeQDDhNaxtg;pB#4u;I=ccveuh<&VuRN9sD84Z>=?68V=9R+#zab?*yEm z$}X1PI1dbKNan-h>EL&5_l$DW6^L1#?r_wQ705n4H2N`#1D3ymw?3rEz|0E0%2)Y( z@JXUk@;P1_XbsrE$Y3`+d|f3N&NMnm9B_#{A+V1F)IE*tOtJ4EmT$T6k{jcL4@dL( ztt*}q=L8Quw8!|Nv!7Dv;?qOK_Sb>+HqSZW9It~SS!k++>skDHVT1O3S2;HhvkoZlHZ8~5|t;?bugeJ6!*yMwkSBjE8!PC zIz%{1JW!!2^b(iWhmSAa!VLTa9OPgU8+2aDaZIdS2yh2P-%~ll18OSIZ`Hld1YPgx zs$C=J7H-j2WEsj62BvCbcFh|^;ee9a-KLj9FmPY%iz74Gp|U{a+rpEA&{v1Ay-K5x zc=1@k*yQ6~{MkF0=66mKUP%bzvbrY*A~)-3U)vxKY}J+-Y`cnq=tXyVM}#DRXzzjb zL6Yp?NS3p}`3o%Y*uB1zB;W7&{+rKB2{AsX{7ok`&7hk&*nUzvkf#?vE$}Jt2I$6@ z==SeQEoX)z}0R=L+hTK)%aJ177Zq(+DZ4^|9LqX~nNK;ks>90@3Rl=~vn8XoY~^=a+0 zm7~PmEe;|gfdatlhW(eUWC?Kd;;DGa0tsmA73;C8?gy^I;#bt7^Pada__5_d*5}(PQHBLk|+#Kz+eMoejNm;x% zrG*fDGw3$sLJtma<$vLLDhLh>SDeo?!pLuAdJRf*y&)PRRUV1@R1--J51i8)eiBX3 zQ?5Pv%ma(xX?^(k_B}E4l(@r^?@VAI;YBbHtrxfax#efuks(}4FtL+;FDrzfmA@~L z9m4ItEnJiO8Uq>wZ`4Q1#__8nC-w`o4&t<@g?x?$w{0k0o9-6INCXR#sP;0sHRl&#o_a0K9Gw6Uy9t;J)&0vu1fa!u6_2 z_6x%Nfz|z?PkX(mL!}3SDb^L9-~h4PV}HP6n4oiPe#fvYUYzB=NFu z+%e4=TwS-hUH+6KFt08+7&mMJ3*N`r5;Zp9&Kd9;Bw4}lCo@jc!X3aVQ|l1l4LWd3 z;Pu$16cwPmYn#M|S9ZX3HMjOs16~-hqq3yELJPzNDh{m8QUKeWpE?%LFauio2U>OS zn?o4Z(3ic<9~8)zewo#223VY%OC)vO80Y_QeC|_49~w#1&RH)8j(s$e%UUo*yvffP zw%}8Qfof-KZs~s{ZuXekKH97V*f!^C!lqGTn5R-_E6;nJ?Ri(D#&a6bmUtK6`+@~1 zb!s)_$$ld4UsB}!@md6~bCotCpXs4$=&`OtPbA>o139p*krUQvc#lP0k_No@Sv>-F z%fOqex0d!k6@p2PQ3=_dG+^<`{1fwb8a({5D<^Ct24q$I6g$32g6S;%BV47du;f#G zSarN0kU3UvSYs>?UUR36>nV-kiw`W<9el(BEQ)`~ylLzve(rW(xtg?(RbTPzrach^ z*4nO%OPH8IN%O0%K|6SWsZ~e`cl%G`!Z%@#*Mc(OCe!+tW9EFo_}jW^?{-N6E%Pr& zg=FOb=avtj@9f4Pb4A7-DTg*9qE^0OrX(xS6u)ystceXOv(xSk^C<#-Z9!J~98qAg z?CZg!snYNuYhKktAp!7VBPjIpmjHE@;qlMi_+gB7iS45nSr`f$j7H+e2-ks)u3=mG z;hu;`rLUDaK^$wT|8<)YqT6b{!1GR3sQ+f?#rwDAf%x#JW9re0K<;V7>Sf+HaYYk3 zJ0PtLlEsgfx0uO+Y29Cg#eRyzthL~Yu~QErw)TE!-8FHz`C|THrRXs6VV2u;z6(qs z^KiEiw(2MDzmJu3@F71uWL!-|=P3aXyvW_yas(IfF!Ep0!wTkCcpq!}DglExuo~fA zoS=J6FJ4{B4$l~`PkMP+68f7ah3yp-Vm$x<=I6iWc!W=**8tIIE>>0$CJY0653aaX z^aOAB`Ru;CTORzJ@kP!iVhrzEXrC+W+e19k;0$ZWq+su97+O+8p7&{N{`~1e&{p?E z0Vc)&T{?4XuLUV#jsYWA-DsUi?#$V6W;^g~#v?XM=qITv1qS8RBKA z_znN0U_KHa$p(sF2Zi>^NW&2OBH;rE_<=%)-W}a#a&V{rQC_<*(m+j3L$o}E2|9;; zI+(IR0=O{qusxRHg=`rDk@g={VZLGPc5ICjsLo5wDt$UY6rA)>F5_bXP6}r$RT^Z# zos4fExB7~Mi~Me!tyxlFm+twR01jywsYer4EFU0vrH@Xp*Wv?ft?kW|1BQq-;;udX zUwMJ*%=LG4Zn43#g==rAmdd~#f(knxU*!VA$?^vIr`UkVtM(_)hF%ffC-O3cnA?eI z&6|u3blE}SZFqe{i743rye@3>7B2W$#3H-B=N^%gK3&K$rxus@aq(2x^&B5La`g>s zEHhjOr(O5w?;<({JyniYaDm3AXm``TVS-29y?k5M4}wSj@chUMF_1eas?2_e6uf5; zcVve&7r5hQU_?_>h5ca$wY9vg01j}-Y#S8-w@zrNx_!sMw=Z6eWut=dP3R-hF>g+| zA|!sf=Vli0e)lDwfRobzSK5=VcJg}in>YXe^~Gm=tBUchs{hKjs{YrXKnwvG3NRF4 zD8NvFp#Vbxh5`%)7z!{HU?}jLD8Lw(&ls2gUm2JGn>2GqG5U>dv02P}xCORx~UnmNq4~I$!3JVF14T}u9okI`D zgb?YaVyr0zHl=O|`OgPI3qJPDv7QagG4iMAb!yzT2}qKS*L{C^hux0mRcn zX{?{h6cwyLo?j6Y6Z|8f^H94I#zqw#?e?^JstUG{$`mC@z`azl03k8})rCDLMh%Oj z@<0J9pG;<8%G`hCap{5VY|L+pggQZ?v$3El61IO3(V~yCL!-O0mT2l?c~fEp)s`4y z!c)qz6HhR~Vx|O~*#5u?`#B}a*PeAQ*oG;phv(P1VYXyL!nEmnw7cqf9+_2YL5BQ| zGj*|EOb%}o<0Q8u0W4Uoys)F6VJYD%yhRI@^jS4=kT9A-?lka_eL@>pm15*E+~)2$BKl zM6}UfUXO`UnWFRmL6*l0tYS)%ImPNNm;!YTL3s?_#$IBAQ)+uxp!g-$NRXIwQ(T~MXDcHf2#CPl%&mVc%U2Xn>6;FvVx0Q1k zHe|Xw457}*evQYi%N>8o{qw@r^sxUJ0l)jOo@mtUFH>INmx9k89M&Tr2bM}52b9Fuqj@)Wc}hu(b}h-nf~QC*Jg4vgW-+B4?9Ro9up3kEGDVIfA$-CvuxdSPZ_VnVj-#8^ry z1q(06MCC`dI%;A7HfD^8g@cncSy`Bwi+)UmVyPsk2!B1$o|i1kpJoC0Kb*g*k%PJ! zV%Ua>H94EWL_WXxiCPtu_@7~H)_<(@^XV0s44K?c*K#tH zL52bh1sDo26ksU8P~bOHfV#n_k_%<%Qltz#imZr~%}n_zYmy>myiue~0E(0yN|Ca@ zDN^PiMaq<Ao-9nLGt-@g5*Qi1m`0fAQ~brKwN}Ky|$Z5F4l-vh^B}O5zP?I z5iJlc5p57{5$zBcBibW6AUYyCAvzDf6fq2OIpPY$aKs42NW>__XvCF>s}N%lV-e#J;}H`OS0k=RT!Xk4 zF$pmlaUEg`Vk%-9Viw{?#7xA^h+7c1BBoPxeAlS;7W#s-DLMM-!A?Nf*TSo%!4Y0B z@(MCs>jkT~*p6Ir@&_8p^B%Y8c|fJS(~l}nxI>wbj+ciU{GqA0zW$RMYgneAeW~?2Pta1i@&Wg951=*@&mw5!2BWm?#!lZg0hdGAXqwYKVCq8~ z1KI7tpyGSkB^wn>$d>b3p?PNj(9U69c-naxJnZXsWcPYOXym;zSmydd@J-ihm5<&k zQ1G%)Q*S5`VrG}0?vQeZI;)FZU#SX!eVy=hT(lq*Jj*%q;FJRcC zuIC7y*^QGcHu?jF=`Z=~AKHNVs^O*kj=4kGmKui?H&1ws1~pt04S@SH!%cf$E`%1X zn=iaOvj}KQ>11w=^?+w@-r;g7^#f0f4j*h5wT3nAgYxMo!+@FXr@n$0!QkOLhaV9) z0)fEMw?+z?t}wXP@o|E!I{?>`3J=QKfrhJtOH{XeL5nKM$~_yMz?EjlQ0Jv?;M~%O z^&P7MA#Gsoi|n2Dz|}eaNB9{NQ2V7pZOQayaLy|H{0CzTIJSehzNN_&9F{(8xJldx zeC?UJi@%KV!p9O@PWzk=UKuiHA4Z00t^Kh3NRG-dj;0K8VvF&s3YE{jB^KF;3g;b zRX$)AE(DEsu{jGS^TO=np4iKMEa3I_#)zZ~-NameOuXK*fjBGjT0Y-U2%axl{PTtz z3w*F~RchFakNCs1>lfBHh=9XNTV>*k`tV4PPdjg&t|0W!ZZ!0LB?Qj-m7kg(&kS3J zmbVUn%4q5 zi0Dg;BV&fgiOO4V!awMXz<9P4_G`<^2;XOw3&Z0CLFs4JFGtUe5a9-Ehh%G+pnQ7r zd>vCExH9%gp#YB{JhqIhdTUN2VeeDhIPK?8;<$#B^;NT8qPC#W#MFQp-l=grdfTg) zC|+sU;=Q$l2$er@>m0;Dn%JPkYJCAXrhUoFa+^3*W9Q1s7ZL`UEK;63FEE3gPl@uX zi-w3M$FtwFzVblkDE(OH=bwq^a-q3VU%K%{L;i0y-!MZKx$r8lCw;`UeGhu@WEuzt zw>;rZak5ugaa+%;`?$KI%~+-2Z9>T4ikqhPcY=AabWgKV3*J{DgW9 z1&Z!tv*VxAVE6r9_6Zn2cv$29G&h_Z^1n#x>l%LqOoqVs4a zO;H5gb~tjQBY+=#8ZFkaILZw4PAC`hENa6a2_4~Cn#&K{M!N1KhOoocL(iMFYsT<1 z@2ZQRJQ0PVZhLOsGit+=9^{lgj~XOa70t26lSEr4uzRd9E!0TQE31Jv^?&k`T zG-2S&&wJF~={u2iVpfWKkTBTgu;YTV4I5Bi=~BtV$q%&JgOr(fa6%8Z+dcN=M9(&I z+vdEv)=j(%k<{BbPJ{2q?}>`obAZ=bmMg=JzT-1R43|7s{Z5o!j5#rvmmhLi(Kf$y z=YzM?R$t5&Cnu(|pVr`ff*BOZ&2hiJbC~FGXMPb=HcaFm-n;cu*k|I}(Y4PTpKycC z+m9U|yU~Y>9>FF3!`MN|hmmPJYQ(|r&d?Vf$;>eAU6#AT7#C;_6?Etwc#7X=DPhUY z!Jv=s>lqg%MIep4x3zW7Tl}&iuUCU2JG|pTn$!XV75H?jwTJieBS$)!z<8CSp*7*0+x+ZzNu- z6)sY2V*y7-x&8L87Jvp0my^od25|ZJwHK6}hj76GOrs!A6n>sxcHqgqulVNuGM<7~ zBZNgbxALr`)1dtO%6%E*7?^iXY`?}D4#<+I=xwOcMucYLv`1GlWz?|GWcvltM!=mJF zJT>!JhY_~~h+dN2Bm8QRIJgBmK2cx;ZF>$j8(-!FMk5Oq%FaoEPg-XL+O$Q$f#Rc; zY{hQ~FZR;;2Sb6Q%(lEYFlWpDB{TRG*M_c4u=MZk#ucu^nnGHVd=6S_y*hH*PQ;XB~ zc#qp`d&5;8$p?++51SrI|A<>zu4-Fb^&Gb-O&iT`;fD^z%AXX@@d7nFtoqy7cVfX} zyZ!5pXb?;vNj1Ixli-LMd0sof0u%k0ztwl(fCKXvohS^tMZ~E1i-rwJLY3tuikIc` z@qW*Q;cO*782c>Q-)^@M{F)~*wBzv`+_dTQ?RhPNFgtI%{lM)vL{Ho*hnu?v$@4RA zZT8k0T-YPHDWFCK8lM;TldEbWDzoMG9B83Imz!Io^lor~Yho#DYVPqv+YYa*J)=*cOA5&6{ zp%4i0S@_dTYlw(=A(q2(rI#Si>b&J!aEsWhs`oZn_%`wWgsX5apuscTr@Sxuf5+$k z=5#dLZuw9rA&oPH^_9Lx zTk<*?Ieg&j0>}y&2#*ZwJJkJEQnrg;g(F6eRoO0L2LuL5u>&q9GN3+1= zAGfZM&*&yBEMyEef8v32-&QRX7^o-ePnHe$9Q^tJw09-&RBhkiOy+qkQz;4|L!q-x z8Ol6_44H>AWhg_%Yt$S{v!XPp(4>LRMtV&mDQPgIsEi>*$^Tpqmsj`7`~BW~zuy1c zK3R8I>+HS0d+oi}*=z5!k-uavC1vo?r(b$5MSiu_)d|m~NQ^U=a*miwp&W@4zYquz z2oMMm2oMPTng|T8d{UbDwk~#T{U$xOn3*5@M-zb^25aDtJqsZWd+NcD&xFM`q>)xk zlU~Dg=+ngiZu5M+YA^)+bE~C<6^9}ke05Z$7Z~$7D%*JI^N$xMkycB6bsXX#ua>Ia zm__@SM5U^;hl+~BVERQfd>QQ1E=@El_0kwpo8oGzZ2WiJdPwk)Uc_ptAsmRm1OfyC z1OfyC1OmSz0?fMlDJoX3;P$l`v!zc|;Y{t$y=r~VFp(jM`>u}-&@6AMC}NrieKitF zz;qjMWKKjr_gXs`dg%TFy;%s@tQx;FxlbOPGkh=h@VFy5JzML!{SqB`rMBjJ^wov1 zV!7|MR0Sho=r-+=SA`<{6uB@;E8Z0xv?}Nh_p^k3#;R`yQf#52%X-%GDdtePHD*)B zGgs*8Y|QAd;R^Tkdhy!iSOfOlG|3}w24Ed4FXEh_f*m!)C-#9!4<4^lHoT0^0AUyT z1Dn!a;B}X4g4=?u!Lq)~Te1z!p+@U2GxGsmV7q1OaxZODFny;dvT2S8V3to(FIBRH zt?33?ZBdrM(O$G>6H}k-U5=Z~mV>&?D{=Qsd82Jbq6yzcXO5(>FB*p+l|CNU`OG-IW;|;bmX4RjzuP;vG$zc zUMq9^4eaRObqBLbCn<`<$^J!Yiq&O^#iWK$iLPQWL_rdAtmX!eyk;ko^M$~dDmt}5 z(X zl)vn7KbY%tL=&9D;yT0Zf!=U7dKV)pk$P#^{@b45+c zdg$P>_W5nK2Rc#z`v+>1lX!sH>op&|FMUSErZ(l;b8v&j1JP?$nx}#8no}X3g4`fD z2BG5utdMnb18d`jw@9L77;DLrTS#^d8ae%=G?1QpvP5wnJrrR*9Z+P?0vUElc{OJ| zLKb=0ugtgof=Fo{Jml^q1Q;)tDIF8y0yE~ODQ;ZB0gl%oRc=ZWux7(r$vP$$;LvqB zZT&$3sM(S=Yu82|uv~EJ(&d}Q0Z-$yXRU`>;mNcCg|#AdFqgi3hNkRPXuw_KSG_#G-FV`|OZr}%(Seab-PjkWoS7|=2`bKnS{yUR( ziPGTmuDw}OrNVHr;ASZ`u4iaP_14VJXQE(P+4>+?{>fl&-P1KZ_7cG8beZy*7#V16 zrWPWQD*@$ow)e*dOaaVC)vncUU;@>0cP%2q9-({MwxquangXq>U-q~2^`H^@r?%KX z?L@E0b!c4MAO=^wZ!0?bkOpMZqTADY*kRN95aXBv4k)B=pI)+t4@jjowVW`(j?9|+ zf~_3Oq*vz%EPT~Z1Dc0 zd~rIz{9H=>y-?N(&!uRMGnaCim`nK;V`QSv1OfyC1OfyC1On6$_)GFq#U{UeETvL( z_V7ewsjpV!#8TIYSc;l062B1$5C{+m5C{+m{E7%%-J@B4_5>R&c{af%y-pYp(|TnCbFLopb;BjHqOEk6kJ?C#m|59gqfPTU72EsS zpv7jT`p;rQaC>=`@7*|Ypl7E4;=C9)xO}TCajz#I$gbNnz>(35CYcxV8(}*Z{F_9= zx8*Ux(hw>i%R z*7L7j?WijX*ZTB6SI*!Eyh~)RnhJEH&tLCPFn%)?HgZ1srp+e|-geXn%~KGA5A*J- zUVA(V#(aI#k+Y)#Wp?*#nzfY|B!2{j@71}%q0-pGb5<+>DZX7UGk+4;`;a5(>;(bP z-)G@+ai;)G_B#FB7fpJ&%j!mwr(G{{tK`XKp)cH^J+Ig!7SMtz9%p5;4IU#_wubkeCFo5Td9aECjSs{J;uGCow1N41# zu8&Vl0y?bHGdOx*0xknfj}#nd23p(GW}NWmg%yn-6~fKwVCL53hJ71(;MvSAkHrFM zU{-15eW#3iG_cY(!m6GNEUjbWD$RO}1~D-$v2){tYR3(Hgu6c@1CDFl&-*tZhz0NL zTVDJy?}%aTZf-H4y^p!bJVO{Jc7e~wPKdxuwtRAF3fV|`hpGPiCj+R(rjEV0Yq??E z+CGULlfEH{)T_68C%K{CeWjNJmW~eWUH0#t1j=cBq~a$qjE^ zdhnbp?IT+MYzrTK7B6g^yi=g+Mn5`#>d74)l1%XGsdUeTX~KX>Z{e46pX=z?lh>RM z1_*%IzG=GIRg!R%uk{NJRxaSnC~|U+pD<8Z@P1zRgIeU2F{gcrV?UC0&pcA+>1)JI zG@xzWhZ>|IRYFP1QxpoX&pGYztqYm`N#lr3-fJYqgW-H`EF0LoP0IP~9dWqD#PUHy zf;70JhBB1iU;t~Ie3CBQ;Q*?6)15vw@qNQQ*?W_3Bn8FSy~b_|BDM?I?>d&z0cyJ?NU&bzQN$M4_1bO{p7=?Z~47 zzXMSrEhurMJaMG_e{!TeVL!jN{ruyrr9YA$Z#n*~E?o3XMQ~ zO&~xZKp;RMKp^l7An;3{AP_0=i7UN%&jILs`5(LJ&J)hFC08ev9;4yjFXPi#K^JKkGJ@Q_|aggJx*@s~h= zK!8AiK!8AiK!8AiK!8AiK!8AiK!8BtS3uw|nM>Id^UKerFunY|e8O`nE615jDI(@l zeuX|k)R91dK!8AiK!8ACA_)8?u~hAbUp|&PtNi5NL}MxL9pl7OcZpbPBDP8J5(p3o z5C{+m5C{+m5C{+m5D5G)BT)J2(_ZFVYjg=p;6T{H|G$3lkMD_{_~1Blg3y0^Y4V%o zpPvF^8``oOj-8ui;o}qL9ggi=mB1JGO410i&tCZ>QUZl#0-fdGL3fdGL3fx!PZ0{;t}-$;*OkS9IYgMKKk@*>4kVC115 zs%$i>1N7v0iiYMRU{@sZGm_)5aU`%Ohp3VyMn!IR*vXeRCgaF`1e1m531CgQ75g8D z{ojVu?Ks_m)15d?z-b~*lW@8Rr+aaljMIHMO~L5_oF2sKA)Fq@=@FbB#pyAers6aW zr^j)60;lOXJ&Ds(I6aNiGdRt_=~o_HC zHjKmm=i)REr#Eqm;uPXEAE$S5T8PtII4#2IZJhp&(>plD2=O8i_;U!nO=eDAWN8ie zo6@$6*to#`!KM#f_Sr*`dE)z1rnrGsv#ZnVZ+pOgA#E;G>t*1+on+!J-WBlSrZboO zAL#*3+ILlpb}xen%I3NIh;^r~~G5&hz&JcM} z`@waV@E31zY=Q5q!tNlD-_ur673%>9A_aHJUT}h^W_&YzDd_-Ozoos;QSkxNS1uN; z@il@TQw>iz=G(y-#Yu48Y8P}*#tx$B z8FAGFFK}#|ahd%CSD@72wC%Sm#vo~hQ-f5J8>|-++UQzg1FW1Xd9!BdL&cI}-;Hy8 zz#R=b;kxaPV4n5K-sMJaKyUfk>}1J6FxUSY&&XPJ~yQLch$7*itcm%-g znL7l`EUbY|g}C|T4H^@E@YbKR5<=?Violy$GgH!}U!(qr_f$S%9$?{iBm2r$K_D?} z{?V@6yx_R^<)}RzeTY9|){zs>04@Oyz*or*4b9x6>qS|>Ij`Sdoe-jdC9Y0dr|LMN zuH4JCX-xx2YFBhhj4LbL$~XTV>!warwPUvyhguKXR9(b*x0w#wxo+O4vbz%r%*|~Q zc*73wsvb=Yxj_%#6?yG^(=!PstbGx)%dZlN>6a7wX2B27#~%CL)`%TSEO*+Jx``jQ zeR=1UGr$XrUwSQ&v+6|LCU=;c&L}`4JWlMPInafspS+y8cy~8awMdBl?KWYkvCQT} zs2UrbleT?b0jB^+C}xuF)T4!aceXAu6lRAjT=h~wMmLhRjydEyFF(k6{?KD9V*?VK zs`z}8cR!-yEfrW2%?x?tFXq3l{fx*M>|1=Qg&VxH*W|Lk#t+Y?JMkuF2*PO#Uvj1u zzeb{Oa#@}E)P-D-EDG*)m?#K5>p=tR4BT$ywxe#&d`1BU^xy{T zE8F+<+%Wr-y7mg1XXta2xLw(9^w9Rjk&MFmJRo6?=9(&g2ACbYX!DU-^dRF4ARtmsn`#z>SS9GBrb23goo;4Zz1uFp=gL-6A-bb&MUUX2`$6f$ru>qeBUzz07 z9cZGMiYJ=Zh8{SO?i`hqfxLbRt=Z1!zhDzO}CP?|LfH!Y!B9#llXMK6BNy1-qW3%!sOL zt5Yo^-uZ1-F!t>3H_QuP*A$i^W=iWAZ#zFn16JrXRpd@4`v2e5PyVgN6sO)Yw!Ufu zUbr4%T&cbgtjXSEwN=*)%%h^@YUpIu0jKt9^rHXJ4|La$&6)JEOYnt-Yyan*k64{+! zy=TMLndW`Dm-L~R(W*@T1r{KeSI=fcqBXp6?8uxC?{opX^K00n$_p!(Jo(o0XbJH1 zxV~_rSGoP>iW2C|*ncf@w-m5>yY=%4IRhZw&=sv- zya1jua-i9})eRhaR@AT1WI*))e{0A5oAhFrt!CPq!v%GobIR!0OoFOcV@IuogGJG2{h2UrS?zoK=2Z-;0Qa)SxVaPqzI^png#Cgsww)!wh zusus+u0!27RCQpA(93*Apuhaetam4wzze!dm5=s)L=Je%r@xaj+JIrG`8d0C$>Kw8SJgA}%}GOpH^w;XExCbVIKgsP1(PdL!}{xqo`yO}bJ+ z5WhOs!=802us$3&S0J4cN-EEPUE9=+-fWm{EAWXO*jkLFwe_NGxi9=rAFpFh{`l@P}xV^Crb?UueU1{ zVx|Quaw#I&pWh?A8oCTin8o09yQU=0_I~u%$I!kxZgc?ct`~K__zg)kP+jvNj0du+ z20qs|YDZ$ty&8GtdBB#4!wx+$?~#2R=`YK-2t&r1a`iPLG(gr}qk1~SXSAb3`=&Fa zIQ*N;gec};ia;Avl1A6DYBWa4)Ohz*ZYWu@-%l=!8Lqw_bdw)DGn?-Im3Q{FY&s(VDr~GR#D9B$oDfuU7LjBGaM4%W({+Mie(4@3ri0Re z_fUL*TlEJt!h5-0XFVV6&{>wS`mHED?Q%(=)Qvw;otQ(0;HDo`^Fo4+CJ_u#fV1uPsE9Dcs9-&_H7FTr; zX22FPV0iAS7_c#2wI(2<1Ks{9;8b3gB(R^w=B*PZ0#(``YT5MNLio?kV$1aDM*^A| z^F5aLBIhq$dw$wi8182AxmJ_kjW#3&GoHWMhFBlHB zQ_Wv{jhBh_} zXl64GRf85gh_suT^q!WNhNczUn1H>LrlG0k8hkg4hUO{T;JXkyn%z94ciAr+%my!^ z-b(t2u}}CC3C#YbhQz5B2t;cy+u|*bWQ)IM=A9f$otM*@q|QM=Fo8Q>%RT zl%x8a5{uB$m#0$s;j`{Xu}`S-(M!i?|3qXJh8vxre+vR*`#(EYW@FXi9V1;R%-?ru z2<4qqy;)x^f|HhlvHX;~sNSPX3f+@g{=gqau}hz(&2gUDGzsI!M*vcNC|75VrEN%Z zf#q0s%C&i>$K2EXmVYce<;sLj+h!73!i}*{8X6q2$^^Rq96VzsC|8LL&{zq|wIV|_ zMk0)I8N+t<2R*_lmqAG)nqtF|C0Isd^)pR3Mu4&wPSu+XFIyWi zZH%3oH*wRWN37z-$i{^6Q}HH>y&SZ~s3#n;>~N~?Vfd)>5q8YcQgMZg2cp-iG*9~p zZ!)(3lg@^l&ek0`qBN$$(h$mxFucZpEI;Ko7+xDZq6X>ul#Sv$kVx#5Yx7-rL1L#| z8L~sigPjoXQ2$`BpRnm+e(Dy~re+-Kq8^G_1a6@r z6y?J(8zS>li!sz~UFk~pX6BRmsaw(6hOeSKZV8Q2V4QG};f?%!?K`ztL^WI?>ha)r z<`WW0wHJ^_mEr4Cj|hjCdXrmVoRDz0ToR*;`u$M~U;o4;z*NItQ zi`@mV4oUgC2V#zq8W6ZCVv332B2l5if$qLR-r@e9OWlL8_W@BV;aZH8UIi0P27PPsDVn?T*2vDvPxigLy_^H5Vi~!{l z$lY;_z|Vv=!%a*=iE^8fgM!iglrbXTPmkuOtN{5Ic9@^!uPN7O*kcSUPD3-c|C7G& z;FzyFGGMv~g!s8rQSe}6X|aiR0p@~m1cwADYslWHT6cZT&in@cd4Ffu{ z61av&3J>pa_od^HFyx^VA0GzgsEtYVCaYTq2HLQy;b0C*PzJN+P{H{i>DPbYr>uJ9 zj>I-WBX(@e8-0Dsn$KeE$~HEp!T50vJpJ55eZ80JQ;bmd#fd-hJ|`3=4$9{DY}RevK!zRX#jegem_6G0HkUNzmqaYrF&#!yU}1 z(kLT$GQQA{iD^Fkkrd^!V!xd0`*h9^D~o+dGqN+I+^CK#9THx9tnKA}JSD{VsWN;*Kgz>N~9{L*ek7>bprsQV6-<55F@s zHW@YlNBS169eOw$`U zKgf4XZG1z6BSQ4B;=M2l?mtKv28RZ^hY#}qs8TFL?7<=4L6&9)gKSKAy_6M|6-jmA z!I}_5184D$4EGKqbvvvi8q)XJOLF9SeWWHpNmnOc1Ok5z0)zd3)sPL-4BGR}qxtG1j&GdX$TNq|O5VwD!G<_I z(xrYE`O$bC*Yl?7n0FXqr{?g6H6CH7eh)dBQMYox&va~P$G!GAw~&>xyhUMRspQI2 zzk_7l#ua4IFsky@?|QzgJoWp>VRnKN2>hK0jP3uVtitX_t+0_UG)`6_ZCUZIeNsQ4 zRXFbZphrlPjod{CuTMD|lFKM{)m31{`0)XQdRBqMA|`iX+)rl}8oe$?NJag0R)H^s zK1;)e539_W`Fkoiab$SLN>Gk8$pDR&_!)@CNc_~z_gD$aWegkSA2bj9ncHvtZgVID zLpho!H{BQk%1DtLZ;SwCCCJT(+YH@40Uw*+ zNSzahWBDn!LCcQ)`Ler7{J0vV>r*z0??584Q`Y*s@PfonxiVyjj>ilCbXGx9XXi>9 zX$i8)kQVX zr5+D{XFmVstitNluxD1GqZJq@B+S2`@!mQy4<#!=JtidE^RU7{8x#)nQ;!PCQ)9#a zk7@%Cj}jB0WBotrK46tcIol8BRH$SXlH}7rd|gw6HT(bTtOB{a;Eg+&g`hnCkgaog zFr=DQAUmW{?34$?VHZFiYN=)w$Tv@TUmq{4Kpw2d2vDvNd2kpbKp7u$XB;i?Q-RGG z0m>x|_mhzUZcyo;3TuX&n0i)$92AV^r;HK#etI-NWd+E$u*3W$e@)q;4ttDY#qq4d zi2o;jzqL$DXatb)vI$ij+0XO v-)Th|3e@|x{>msxN?C<5{hz;gC~Lsqh0}JNw&3)`5T%e*7+AOCX~q8m!=Xy3 literal 0 HcmV?d00001 diff --git a/docs/examples/data/impactx/particles_lost.h5 b/docs/examples/data/impactx/particles_lost.h5 new file mode 100644 index 0000000000000000000000000000000000000000..ba0357aa7c0ebc7863177910be7b68989ef25ba3 GIT binary patch literal 46655 zcmeHQ30xFM)_*Vvf+%=FOcX&uP@{n288eLtBZx!-NEAgJVSq6(jLgWP7;%j`jPd$N zG>S*!kq9Q{@<>c<(TJK|qMujdkwlSoz2kkw$W~2NGd)TVGiHekTlM?Vue*tY(3{~Tr^Sqm*->)HURc6 zfQVg+Qm@=eW!km^JAgh*$>nZJZCl}Mp&o&PKNYf^R7TxDA)!+3bx=i}v7$Bs@knVy z3b(wjQ}V>Bv4mZ0H61WDjx_Mxd|r^Q{9*E zfBrZsk>;GBQf3y|(lWyyF8eXB)81I1F{tDv_9J}$tBtvqBXHy*pTfzp38D7QF?hD> zlo{v*C;QeoqiJSVrb?rqoEhEGBwL%MQ|r}PnpXzV>~@&Dg^)cWHC3n5uOCfslCmdd z?o6Z^#CdcJ)95vUNvBn7{IY2lY0O-eIxStV)}+PHBl%z>I2^tHUnB#z%=}Fe6CMqC z#IXhB$wCsmdr=FL<4>CxJShi))sQHCiC)%4khf5W!R?HbAEG`>h@35a={IC*NL+X* z(g^u^OLQ;AZ-#?#SEyxQ?{Zf~EaHNrv^NOq36x0S%Az3;)gjcjoA%TpltG|L=Qlt| z+;9b|cRn^b5RZuf^A?>4*o!^67s>E6_cE%xu|{9rOJgJ$@g8{8GVos4B{uJGsyNg@b?ROvlIC4#W(Oq&2$ zg=?Xh&L8dH2NsDqRJ%%_gVrKF=E`#igHDw3(VU7H_+2D&`RXfiFjItj#>(&GVL0}+ z-UN@TkwFP7dVHclIT9%(F5(5${%*K4PyiIN@Ud_jquPg*^7hT{Yi3X9Kucf=-l*pTJ?CQQ7X5nr| zMauej{~i>sD1&yy$leq8gBvcC_Wz3cqmIEwkx1a0?#ICkk0G>-ySw28xQeXZ>$2tq z9Kw3gLXVfv{ssOi>fP_n*PVjiqJF1&!Mzf`6OA&b{5MvDA^f3(O3$TURUi{_WbWK_ z7VK~tw3cV0Ub+P9&^7cY5Uzai;x%|mhBBxv@47kMfEEXO|N46B@%7_Q}{M*-RZ@1#=}Svh(V{aCcse4LCyOtXLd5g ziag}D|7;b!E1D33w!D`P({an7p#jej-}Eb(2)S4_YMxElwu!J>)R^PX?VAX0BA4l( zSdoe5F=Mqe-YL^S4?N6LP348_v!F->bLhKWbKnIGA|1NIHa;kT9E=%lfA12TY0yhF zx&IWl8W`a9G3< z_Epf^P>9=sii{ie<3jjKr0MDfqZdJ{2*TEZtCoU7q=$X%nq^QXnq^kMcy}3~S+vHu zh@8Kz1ewUV_doP}9}Z%WsrhG?SFCD3c@r*#q}qdE6eyw?$JWNi2Fv0kP?<9r^03zd z1Va<_eE@s9IwLd&KBVAwUw!E^CAau3$FM4N6H_G+}KI8u7+Qa=FeaG|4{9rL~j4l9TU{|J9E1#VicFh zB=XhDxjFusN&ND0|JP~5fi_&mm5Zmn)_d@KyAMdxtnJ9ne7*B;wtJoV8HbMfY>IW` zZukD7-NoT8xobrMl|^N(xiM9RC$na{a>qx%dSpY%GyLZ#`SLTem7re#7>u6V`{6Z97^N_?B0M6W^&@>fi3iJ9670=JqtZ6LLHMNh7mt9QjkjgL^n8 zp5C@7&_3bYHot7!7vH?{kG^HY4gLM~-%qc& zE9cL9-+JETu6)+!HQ$wgDC5%yb&+5CrXT;3qr<(tub$%~?Oc0CcjfLXjnM+#V@W2 zMZs|qvDxduKi!d3 zX-}^6{JWZ5=TAk8ZaCGxsDl@1c0WfrHWPx5%Vux-!Ri&JHegfRGLFegQ?)0l1IJr+ ztu$G=LZl)k*ZXmN+^KkLSD2#QP1VBYQjM{TFq;?x3;~7!Lx3S*7J+(+r52$d-s04NV<}BoV5+sTRMb<% zY~je(r354xWe6|?7y=9dhCt&*!0_fX!1|wcKf$-9ZzUhd{S$9XTRu-(Q+QjNNWZb7 z*ZJ#4ko(eaOK*-M9LBe$gEqa~zFwlz9{A^GC^g`y^r`k^6ReF&D_B%Y+MOL40t^9$ z07HNw@IQn=y~I-aLkWTnIF?$M)3U9#u~bw$Q^#Q|e5QD?D~13=fFZyTUpeXy)dqgs#q=1;sW zF&42x21RavnOj;=}_G}_?m9$RbB0jwq=unJ=sev zN0~8MRr}|oS@JAGXVpw=&8*MW?2wL)e_pjqf&#~mm(E6POUD8PDe`+T|ytR+H)H!)p zsI^aqQBRweC=3|v`)*TQ*cC&7A;1t|2rvX169Nr%2O+etHTu6h8Cg2LPdaH1=?A$< zlL$x6gWNM_8*uGYHe^$y|8#^oSWW7+B*OeJ-9hL^v)$J3QOzPuv;M|zG6Wa`3;~7! zL!dDt@WkIeSP$x-`%)`!Yc{D;_%h3(z~p*LEZ$=qCFziL#azNs6S6|lx$%m6`9;i+ zI#N{wj!N$g<4dew_k8FffUSEHH(*DG07HNwzz|>vG-d?qC6?+ujX>6bW2w!t$_LiQ zQWY$gYRp}Q*~SoH2rvW~0t|s#5vZ3~D(f{;R|Afv7F4t?wKkTDYHtcLTbDA$gIzHM z7y=9dh5$pLF(S~A-r%g6`nDf&9TPQSuA_r&~72)`cv59&H30Gxv! z`t*1t`tQ}T5p~F2;p~om@jww98wUPO2#_{T0JjYP%-W_Ez-1ZEwdn-#3;Ik1oplE| z+{|#MwT1aD@!6-x&qNuDxU>{sBBn%JSY`NXzWb~#VT$mG5XiM7w+u(_b8^9bLTL%5 zYbpJJ(1$gRMO+~0Cv1u!P@z0K zh=x4yLzm=^(-Lw=c&#naiv`nLg~Sb4pn4Br4zz_A zK%R#=4jwy@2;MjY^%0!0WvS1bfHVVKE6x}KY^9_tQX;!{KuUw zkbyH$ufZZC^QZ4X-@K!Y&Lm^7`?1PNu4!8_4@>ogfFe9J1p3fH!%#<8^*(zUv_MXQ zoK~yM)~k~bE@W0--_KSilm`uoITz9xhZ z<){pW;ZLb>VZ8!+_3aTbphrNUXRlrZd-oX_80gs}pkF`$l4(U`)OuC6Qm@X^glUqq zQq-C>6K&G6vnFdpQ1vND!u25uMOJpEQg7g|O{u1YJxi<7L`Q`h*zBoZfjt9zVjWt7 zGe(25sPgnGjSi_|7(7ruN}dKg@L3HSx~x=vt}>fgiuWJDmPTZ1mC5>1MlA)dcvZHc zNQ|*ieXqVf1Nx(DmKXx{gn*&{*F*+*W$NUSqq}85KiJGm%2X*ek`YXZw)|toC>O3~ z;~&JYmS&F(`?|Wys~az3&sL>MmM)lF!=9g%smzmP5Ax|(w)cWNQidvz*v%;)J)m^Q z!j(Ja;^NIHzkFHWZnnl7YMrzL5n zpOv_%(dRn%9c}D#27c)-di7-orNce9hxqd)yY%?`;S&x{u|fQFV=-L5v^wJyDW^W1 z6Jg3OtxbfJ+tK9%eN5SEo10C}cKSZ)OYaVRjjJMlGqXw5V3+>imnMkWt(r~HL~CX$JW_VyPC~ISCp3dJ zHSPa6O08VnX-JLH$b=-VbR%S3xp>-Zy$7;_f@I*AZiAymfp2+5I2oZ}wCd|Kc)??Gk0ez9#z)SL`GU$ikiOx58+qdIVBL zBp|B5j1N^S=jQllCK=J!>Me>fjOCm6qS(@*H^!{+6d6!QIJj_0-eINh0~ApS2AGcr zA8{XTo|Oiej|ju}oVc0b>?$&#j*xKsZ&ja1xPK-x0123B0i!<)23r~w3jCHug*AhU z&LVuIyZa+|O#456mR-MCy-5bmiOZ@o)q1^BldMY0mk#Smp?9if;6$X-l$mH^v|v0d zltDWxlAoQGsZ7(T&~ihPQiIMjEQONJdzng|W$4s+CmKP`&ECfTOB>6;FFpJiyNj4# zdi*hZC!LUR`O<@-aY7QYOKT%c0X43JOP5y0`0zvrD?HhmL&rE+nFvTLWE>n!1f*3m zc1E#)MS+cpfOHK)Z>Sle4MyuM32THV#wba*nK393^Gjo7e4ZBbOB-N(gcbPld@Zey zFk=XY1DN)IoPTNi^^zJulo{G|r5OWVK3X^IHh0Q6!V&>#3s3JT(zxZiZhxo-X%mz8 zU2Tg%P!F@K=)sRhaD%Iq{%(hv(G6#DZaDp9x9c-gWA!8b*%LxB)p zi!>Mc_0QQgZrm0m_}2ZaWa!`)ZCHFiX6A#2;h+!P;bc$_()aKR{1AReeGmekJ_rrd#tLY Up@`fG_WS?kME{AgvG)V~8(VX2DF6Tf literal 0 HcmV?d00001 diff --git a/environment.yml b/environment.yml index be31c899..548f1688 100644 --- a/environment.yml +++ b/environment.yml @@ -16,6 +16,7 @@ dependencies: - pytest-cov - pyyaml - pint # test-only: independent oracle for beamphysics.units + - openpmd-api >=0.17 # test-only: reads ImpactX BeamMonitor output - jupyterlab>=3 - nbconvert # Documentation diff --git a/pyproject.toml b/pyproject.toml index dd49cb21..5fb0b595 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,14 @@ readme = { file = "README.md", content-type = "text/markdown" } requires-python = ">=3.10" [project.optional-dependencies] -dev = ["pytest", "pytest-cov", "pytest-benchmark", "pyyaml", "pint"] +dev = [ + "pytest", + "pytest-cov", + "pytest-benchmark", + "pyyaml", + "pint", # test-only: independent oracle for beamphysics.units + "openpmd-api>=0.17", # test-only: reads ImpactX BeamMonitor output +] doc = [ "mkdocs", "mkdocs-jupyter", diff --git a/tests/test_impactx.py b/tests/test_impactx.py new file mode 100644 index 00000000..cbb4fa3c --- /dev/null +++ b/tests/test_impactx.py @@ -0,0 +1,592 @@ +"""Tests for the ImpactX interface. + +The pure conversion tests need nothing but numpy. The reader tests run against +``docs/examples/data/impactx/*.h5``, written by ImpactX 26.08 with the +``generate.py`` script next to them, and need ``openpmd-api`` but not ImpactX. +""" + +from __future__ import annotations + +import pathlib + +import numpy as np +import pytest + +from beamphysics import ParticleGroup +from beamphysics.interfaces.impactx import ( + IMPACTX_TO_PMD_SPECIES, + PARTICLE_STATUS_LOST, + ImpactXRefPart, + UnrepresentableParticleData, + _check_representable, + particle_id_from_idcpu, + beam_monitor_iterations, + impactx_to_particlegroup_data, + particlegroup_to_impactx, + pmd_species_of, + read_beam_monitor, + refpart_from_openpmd, +) +from beamphysics.species import charge_of, e_charge, mass_of +from beamphysics.status import ParticleStatus +from beamphysics.units import c_light + +DATA_DIR = ( + pathlib.Path(__file__).resolve().parent.parent + / "docs" + / "examples" + / "data" + / "impactx" +) +MONITOR = DATA_DIR / "monitor.h5" +PARTICLES_LOST = DATA_DIR / "particles_lost.h5" + +try: + import openpmd_api # noqa: F401 + + HAVE_OPENPMD_API = True +except ImportError: + HAVE_OPENPMD_API = False + +requires_data = pytest.mark.skipif( + not (HAVE_OPENPMD_API and MONITOR.exists()), + reason="needs openpmd-api and the ImpactX test data", +) + +# The reference energy of both the fixtures and the committed test data. +KIN_ENERGY_MeV = 2.0e3 + + +def make_ref(species: str = "electron", kin_energy_MeV: float = KIN_ENERGY_MeV): + """An on-axis reference particle at s = 1.5 m, using beamphysics' own constants.""" + mass_MeV = mass_of(species) / 1.0e6 + gamma = 1.0 + kin_energy_MeV / mass_MeV + return ImpactXRefPart( + x=0.0, + y=0.0, + z=1.5, + t=1.5, + px=0.0, + py=0.0, + pz=np.sqrt(gamma**2 - 1.0), + pt=-gamma, + mass_MeV=mass_MeV, + charge_qe=charge_of(species) / e_charge, + s=1.5, + ) + + +@pytest.fixture +def electron_ref() -> ImpactXRefPart: + return make_ref() + + +@pytest.fixture +def bunch(electron_ref: ImpactXRefPart) -> ParticleGroup: + """A 2 GeV electron bunch in z-coordinates, matched to `electron_ref`.""" + rng = np.random.default_rng(42) + n = 500 + mass_eV = electron_ref.mass_eV + p0 = electron_ref.beta_gamma * mass_eV # eV/c + + px = p0 * rng.normal(0.0, 1.0e-5, n) + py = p0 * rng.normal(0.0, 1.0e-5, n) + pz = p0 * (1.0 + rng.normal(0.0, 2.0e-3, n)) + + return ParticleGroup( + data={ + "x": rng.normal(0.0, 1.0e-4, n), + "y": rng.normal(0.0, 1.0e-4, n), + "z": np.full(n, electron_ref.z), + "px": px, + "py": py, + "pz": pz, + "t": electron_ref.t / c_light + rng.normal(0.0, 3.0e-12, n), + "weight": np.full(n, 1.0e-9 / n), + "status": np.full(n, int(ParticleStatus.ALIVE)), + "species": "electron", + } + ) + + +# --------------------------------------------------------------------------- +# Reference particle +# --------------------------------------------------------------------------- + + +def test_qm_units_relate_by_c_squared(electron_ref): + assert electron_ref.qm_SI == pytest.approx(electron_ref.qm_eV * c_light**2) + # -e / m_e in C/kg + assert electron_ref.qm_SI == pytest.approx(-1.75882e11, rel=1e-5) + + +def test_gamma_and_beta_gamma_agree(electron_ref): + assert electron_ref.gamma == pytest.approx(-electron_ref.pt) + assert electron_ref.beta_gamma == pytest.approx( + np.sqrt(electron_ref.gamma**2 - 1.0) + ) + + +@pytest.mark.parametrize("pmd_name", sorted(IMPACTX_TO_PMD_SPECIES.values())) +def test_species_inference(pmd_name): + assert pmd_species_of(make_ref(pmd_name)) == pmd_name + + +def test_species_inference_refuses_the_unknown(electron_ref): + from dataclasses import replace + + with pytest.raises(ValueError, match="Pass species= explicitly"): + pmd_species_of(replace(electron_ref, mass_MeV=123.456)) + + +# --------------------------------------------------------------------------- +# Coordinate conversion +# --------------------------------------------------------------------------- + + +def test_roundtrip_is_exact(bunch, electron_ref): + data = particlegroup_to_impactx(bunch, electron_ref) + back = ParticleGroup(data=impactx_to_particlegroup_data(data, electron_ref)) + + for key in ("x", "y", "px", "py", "pz", "t", "weight"): + np.testing.assert_allclose( + back[key], bunch[key], rtol=1e-11, atol=0, err_msg=key + ) + assert back.species == bunch.species + + +def test_roundtrip_preserves_beam_statistics(bunch, electron_ref): + data = particlegroup_to_impactx(bunch, electron_ref) + back = ParticleGroup(data=impactx_to_particlegroup_data(data, electron_ref)) + + assert back.charge == pytest.approx(bunch.charge, rel=1e-12) + assert back.norm_emit_x == pytest.approx(bunch.norm_emit_x, rel=1e-9) + assert back.norm_emit_y == pytest.approx(bunch.norm_emit_y, rel=1e-9) + assert back.std("t") == pytest.approx(bunch.std("t"), rel=1e-9) + assert back.avg("energy") == pytest.approx(bunch.avg("energy"), rel=1e-12) + + +def test_result_is_in_z_coordinates(bunch, electron_ref): + """The bunch comes back on the reference plane, with the length in t.""" + data = particlegroup_to_impactx(bunch, electron_ref) + back = ParticleGroup(data=impactx_to_particlegroup_data(data, electron_ref)) + + assert back.in_z_coordinates + np.testing.assert_array_equal(np.unique(back.z), [0.0]) + assert back.std("t") > 0.0 + + +def test_weighting_counts_real_particles(bunch, electron_ref): + data = particlegroup_to_impactx(bunch, electron_ref) + assert data["weighting"].sum() * e_charge == pytest.approx(bunch.charge, rel=1e-12) + + +def test_momentum_t_is_accurate(bunch, electron_ref): + """momentum_t must resolve gamma - gamma_ref, which is 1e-3 of gamma at 2 GeV. + + Errors are measured against the spread of momentum_t itself: a particle sitting + exactly on the reference momentum has momentum_t ~ 0, where a per-particle + relative error is meaningless. + """ + data = particlegroup_to_impactx(bunch, electron_ref) + + mass_eV = np.longdouble(electron_ref.mass_eV) + gamma_ref = np.longdouble(electron_ref.gamma) + beta_gamma = np.longdouble(electron_ref.beta_gamma) + gamma = np.sqrt(1.0 + (np.longdouble(bunch.p) / mass_eV) ** 2) + exact = -(gamma - gamma_ref) / beta_gamma + + naive_gamma = np.sqrt(1.0 + (bunch.p / electron_ref.mass_eV) ** 2) + naive = -(naive_gamma - electron_ref.gamma) / electron_ref.beta_gamma + + scale = float(np.std(exact)) + ours_error = float(np.max(np.abs(data["momentum_t"] - exact))) / scale + naive_error = float(np.max(np.abs(naive - exact))) / scale + assert ours_error < 1e-13 + # the algebraic form must never be worse than the plain difference + assert ours_error <= naive_error + + +def test_input_particlegroup_is_not_mutated(electron_ref): + """A t-coordinate bunch is drifted on a copy, never in place.""" + n = 32 + rng = np.random.default_rng(0) + p0 = electron_ref.beta_gamma * electron_ref.mass_eV + pg = ParticleGroup( + data={ + "x": np.zeros(n), + "y": np.zeros(n), + "z": rng.normal(0.0, 1.0e-3, n), + "px": np.zeros(n), + "py": np.zeros(n), + "pz": np.full(n, p0), + "t": np.zeros(n), + "weight": np.full(n, 1.0e-12), + "status": np.ones(n, dtype=int), + "species": "electron", + } + ) + assert pg.in_t_coordinates + z_before, t_before = pg.z.copy(), pg.t.copy() + + particlegroup_to_impactx(pg, electron_ref) + + np.testing.assert_array_equal(pg.z, z_before) + np.testing.assert_array_equal(pg.t, t_before) + + +def test_mismatched_reference_is_refused(bunch, electron_ref): + """A reference particle from another beam makes pz imaginary; say so.""" + data = particlegroup_to_impactx(bunch, electron_ref) + data["momentum_x"] = np.full(bunch.n_particle, 10.0) + + with pytest.raises(ValueError, match="transverse momentum"): + impactx_to_particlegroup_data(data, electron_ref) + + +def test_status_and_id_are_carried_through(bunch, electron_ref): + data = particlegroup_to_impactx(bunch, electron_ref) + data["id"] = np.arange(bunch.n_particle) + data["status"] = np.zeros(bunch.n_particle, dtype=int) + + back = ParticleGroup(data=impactx_to_particlegroup_data(data, electron_ref)) + np.testing.assert_array_equal(back.id, np.arange(bunch.n_particle)) + assert back.n_alive == 0 + + +# --------------------------------------------------------------------------- +# AMReX id packing +# --------------------------------------------------------------------------- + + +VALID_BIT = np.uint64(1) << np.uint64(63) + + +def pack_idcpu(ids, cpus, valid=True): + """Pack ids and ranks the way AMReX does, for the tests below.""" + ids = np.asarray(ids, dtype=np.uint64) + cpus = np.asarray(cpus, dtype=np.uint64) + packed = (ids << np.uint64(24)) | cpus + return (packed | VALID_BIT) if valid else packed + + +def test_particle_id_keeps_the_whole_idcpu(): + """Only the validity bit comes off; id and rank both stay in the value.""" + packed = pack_idcpu([1, 2, 12345], [0, 3, 7]) + + got, valid = particle_id_from_idcpu(packed) + np.testing.assert_array_equal(valid, [True, True, True]) + np.testing.assert_array_equal(got >> 24, [1, 2, 12345]) # AMReX' own id + np.testing.assert_array_equal(got & 0xFFFFFF, [0, 3, 7]) # originating rank + # int64 is what ParticleGroup stores ids in, and the value must fit + assert got.dtype == np.int64 + assert got.min() >= 0 + # and the original is recoverable + np.testing.assert_array_equal(got.astype(np.uint64) | VALID_BIT, packed) + + +def test_particle_id_survives_a_particlegroup(bunch): + """The raw idcpu overflows the int64 ParticleGroup stores ids in; the id must not. + + Stripping the validity bit is what makes it fit, so this checks the whole path, + not just the dtype: ParticleGroup coerces ids with `_round_to_int_array`, which + wraps silently rather than raising. + """ + n = bunch.n_particle + packed = pack_idcpu(np.arange(1, n + 1), np.zeros(n, dtype=np.uint64)) + assert packed.max() > np.iinfo(np.int64).max + + ids, _ = particle_id_from_idcpu(packed) + data = particlegroup_to_impactx(bunch, make_ref()) + data["id"] = ids + stored = ParticleGroup(data=impactx_to_particlegroup_data(data, make_ref())).id + + np.testing.assert_array_equal(stored, ids) + assert stored.min() > 0 + + +def test_particle_id_is_unique_across_ranks(): + """AMReX counts ids per rank, so the id alone collides on a parallel run.""" + packed = pack_idcpu([1, 1, 2, 2], [0, 1, 0, 1]) + + got, _ = particle_id_from_idcpu(packed) + assert len(np.unique(got)) == 4 + assert len(np.unique(got >> 24)) == 2 # what the AMReX id alone would give + + +def test_particle_id_reports_invalid_particles(): + packed = pack_idcpu([1, 2], [0, 0], valid=False) + + got, valid = particle_id_from_idcpu(packed) + np.testing.assert_array_equal(valid, [False, False]) + np.testing.assert_array_equal(got >> 24, [1, 2]) + + +# --------------------------------------------------------------------------- +# Unrepresentable per-particle data +# --------------------------------------------------------------------------- + + +def test_zero_spin_is_representable(): + _check_representable({name: np.zeros(4) for name in ("spin_x", "spin_y", "spin_z")}) + + +def test_nonzero_spin_refuses_loudly(): + columns = {name: np.zeros(4) for name in ("spin_x", "spin_y", "spin_z")} + columns["spin_z"][2] = 1.0 + with pytest.raises(UnrepresentableParticleData, match="spin"): + _check_representable(columns) + + +def test_runtime_component_refuses_loudly(): + with pytest.raises(UnrepresentableParticleData, match="s_lost"): + _check_representable({"s_lost": np.full(4, 0.25)}) + + +# --------------------------------------------------------------------------- +# Reference particle from openPMD attributes +# --------------------------------------------------------------------------- + + +class _FakeSpecies: + def __init__(self, attributes: dict): + self._attributes = attributes + + @property + def attributes(self): + return list(self._attributes) + + def get_attribute(self, name): + return self._attributes[name] + + +def test_refpart_from_openpmd_rejects_a_foreign_species(): + with pytest.raises(KeyError, match="ImpactX BeamMonitor"): + refpart_from_openpmd(_FakeSpecies({"x_ref": 0.0})) + + +# --------------------------------------------------------------------------- +# BeamMonitor reader +# --------------------------------------------------------------------------- + + +@requires_data +def test_beam_monitor_iterations(): + assert beam_monitor_iterations(MONITOR) == [1, 5] + + +@requires_data +def test_read_beam_monitor_matches_impactx_moments(): + """ImpactX's own reduced beam characteristics must come back out.""" + import openpmd_api as io + + P = read_beam_monitor(MONITOR) + + series = io.Series(str(MONITOR), io.Access.read_only) + beam = series.iterations[5].particles["beam"] + attrs = {name: beam.get_attribute(name) for name in beam.attributes} + series.close() + + assert P.species == "electron" + assert P.n_particle == 184 + assert P.n_alive == P.n_particle + assert P.charge == pytest.approx(abs(attrs["charge_C"]), rel=1e-12) + + assert P.std("x") == pytest.approx(attrs["sigma_x"], rel=1e-12) + assert P.std("y") == pytest.approx(attrs["sigma_y"], rel=1e-12) + assert P.std("t") * c_light == pytest.approx(attrs["sigma_t"], rel=1e-12) + assert P.avg("x") == pytest.approx(attrs["mean_x"], rel=1e-9) + assert P.avg("y") == pytest.approx(attrs["mean_y"], rel=1e-9) + assert P.avg("t") * c_light - attrs["t_ref"] == pytest.approx( + attrs["mean_t"], rel=1e-9 + ) + # beamphysics' covariance is bias-corrected and ImpactX' is not, and with equal + # weights that is exactly a factor n/(n-1) on the emittance + n = P.n_particle + bias = n / (n - 1) + assert P.norm_emit_x == pytest.approx(attrs["emittance_xn"] * bias, rel=1e-9) + assert P.norm_emit_y == pytest.approx(attrs["emittance_yn"] * bias, rel=1e-9) + + +@requires_data +def test_conversion_to_impactx_matches_the_file_it_came_from(): + """Check the write direction against real ImpactX output, not synthetic data. + + Every other conversion test starts from a bunch this module itself built, so a + paired sign or normalization error present in *both* directions would round-trip + cleanly and pass. This one converts a bunch read from the file back to ImpactX + arrays and compares against what ImpactX actually wrote. + """ + import openpmd_api as io + + series = io.Series(str(MONITOR), io.Access.read_only) + beam = series.iterations[5].particles["beam"] + ref = refpart_from_openpmd(beam) + raw = { + "position_x": beam["position"]["x"].load_chunk(), + "position_y": beam["position"]["y"].load_chunk(), + "position_t": beam["position"]["t"].load_chunk(), + "momentum_x": beam["momentum"]["x"].load_chunk(), + "momentum_y": beam["momentum"]["y"].load_chunk(), + "momentum_t": beam["momentum"]["t"].load_chunk(), + "weighting": beam["weighting"][io.Record_Component.SCALAR].load_chunk(), + } + series.flush() + series.close() + + data = particlegroup_to_impactx(read_beam_monitor(MONITOR), ref) + + for key, expected in raw.items(): + expected = np.asarray(expected) + # measure against the spread, or the level for a uniform column like weighting + scale = np.std(expected) or np.abs(np.mean(expected)) or 1.0 + error = np.max(np.abs(data[key] - expected)) / scale + assert error < 1e-12, f"{key}: {error:e}" + + +@requires_data +def test_read_beam_monitor_reference_particle(): + import openpmd_api as io + + series = io.Series(str(MONITOR), io.Access.read_only) + ref = refpart_from_openpmd(series.iterations[5].particles["beam"]) + series.close() + + assert pmd_species_of(ref) == "electron" + assert ref.charge_qe == pytest.approx(-1.0) + assert ref.mass_MeV == pytest.approx(mass_of("electron") / 1e6, rel=1e-8) + assert ref.gamma == pytest.approx(1.0 + KIN_ENERGY_MeV / ref.mass_MeV, rel=1e-8) + assert ref.s == pytest.approx(1.25) + assert ref.t == pytest.approx(ref.s, rel=1e-6) # ultrarelativistic + + +@requires_data +def test_read_beam_monitor_carries_ids(): + P = read_beam_monitor(MONITOR) + assert P.id.dtype.kind == "i" + assert len(np.unique(P.id)) == P.n_particle + # the test data is a serial run, so every rank field is 0 and the AMReX ids are + # 1..200 for the 200 macroparticles the run started with + np.testing.assert_array_equal(P.id & 0xFFFFFF, 0) + assert (P.id >> 24).min() >= 1 + assert (P.id >> 24).max() <= 200 + + +@requires_data +def test_read_beam_monitor_iteration_selection(): + first = read_beam_monitor(MONITOR, iteration=1) + last = read_beam_monitor(MONITOR, iteration=5) + + assert first.n_particle == 200 # before the collimator + assert last.n_particle == 184 + assert first.std("x") != last.std("x") + # the default is the last iteration + assert read_beam_monitor(MONITOR).n_particle == last.n_particle + + +@requires_data +def test_read_beam_monitor_rejects_bad_arguments(): + with pytest.raises(KeyError, match="Iteration 99"): + read_beam_monitor(MONITOR, iteration=99) + with pytest.raises(KeyError, match="always names it 'beam'"): + read_beam_monitor(MONITOR, species_name="particles_lost") + + +@requires_data +def test_particlegroup_from_impactx_classmethod(): + from beamphysics.testing import assert_pg_close + + assert_pg_close(ParticleGroup.from_impactx(MONITOR), read_beam_monitor(MONITOR)) + + +@requires_data +@pytest.mark.parametrize("species", ["proton", "positron"]) +def test_species_cannot_relabel_across_species(species): + """The momenta are normalized by the *reference* mass, so a relabel is not free. + + Read as protons, this 2 GeV electron beam would come back at gamma = 2.35 and no + exception. `positron` is the sharp case: same mass as the reference, opposite + charge, so only checking the mass would let it through. + """ + with pytest.raises(ValueError, match="not the species the reference particle"): + read_beam_monitor(MONITOR, species=species) + + +def test_species_can_name_what_inference_cannot(bunch): + """`species=` is for species ImpactX has no name for, not for converting.""" + from dataclasses import replace + + muon_ref = replace( + make_ref(), + mass_MeV=mass_of("muon") / 1.0e6, + charge_qe=charge_of("muon") / e_charge, + ) + with pytest.raises(ValueError, match="Pass species= explicitly"): + pmd_species_of(muon_ref) + + data = particlegroup_to_impactx(bunch, make_ref()) + back = ParticleGroup( + data=impactx_to_particlegroup_data(data, muon_ref, species="muon") + ) + assert back.species == "muon" + + +requires_lost_data = pytest.mark.skipif( + not (HAVE_OPENPMD_API and PARTICLES_LOST.exists()), + reason="needs openpmd-api and the ImpactX particles_lost test data", +) + + +@requires_lost_data +def test_particles_lost_has_no_reference_particle(): + """ImpactX writes a default-constructed RefPart into particles_lost output.""" + import openpmd_api as io + + series = io.Series(str(PARTICLES_LOST), io.Access.read_only) + ref = refpart_from_openpmd(series.iterations[0].particles["beam"]) + series.close() + + assert ref.mass_MeV == 0.0 + assert ref.gamma == 0.0 + + with pytest.raises(ValueError, match="zeroed reference particle"): + read_beam_monitor(PARTICLES_LOST) + + +@requires_lost_data +def test_particles_lost_carries_an_unrepresentable_component(): + """ImpactX' lost-particle output always has a runtime `s_lost` column.""" + import openpmd_api as io + + series = io.Series(str(MONITOR), io.Access.read_only) + ref = refpart_from_openpmd(series.iterations[5].particles["beam"]) + series.close() + + with pytest.raises(UnrepresentableParticleData, match=r"s_lost"): + read_beam_monitor(PARTICLES_LOST, ref=ref) + + with pytest.warns(UserWarning, match=r"s_lost"): + lost = read_beam_monitor(PARTICLES_LOST, ref=ref, strict=False) + kept = read_beam_monitor(MONITOR, iteration=5) + + assert lost.n_particle == 16 + assert lost.n_particle + kept.n_particle == 200 + # AMReX marks these valid -- they are live entries of the loss container -- so the + # validity bit alone would report them alive, which is exactly backwards + assert lost.n_alive == 0 + assert lost.n_dead == 16 + np.testing.assert_array_equal(lost.status, PARTICLE_STATUS_LOST) + assert kept.n_alive == kept.n_particle + # the same species is used for lost particles; ids partition the original beam + assert set(lost.id).isdisjoint(set(kept.id)) + assert len(set(lost.id) | set(kept.id)) == 200 + + +def test_unphysical_reference_particle_is_refused(bunch): + zero_ref = ImpactXRefPart( + x=0.0, y=0.0, z=0.0, t=0.0, + px=0.0, py=0.0, pz=0.0, pt=0.0, + mass_MeV=0.0, charge_qe=0.0, + ) # fmt: skip + with pytest.raises(ValueError, match="not physical"): + particlegroup_to_impactx(bunch, zero_ref) From efa63c4ac9315f3fa97d6c76bd0e660d4f938a62 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 1 Sep 2026 22:59:01 -0700 Subject: [PATCH 2/5] Point at the real ImpactX and lume-impactx repositories The module docstring had no link for lume-impactx (I had removed a guessed one) and named only the ImpactX documentation site. Add both source repositories: ImpactX is BLAST-ImpactX/impactx and lume-impactx is ax3l/lume-impactx. impactx.readthedocs.io stays as the documentation link, confirmed canonical by ImpactX' own README. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01DhbdaEDSBqFN974t9QcXem --- beamphysics/interfaces/impactx.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/beamphysics/interfaces/impactx.py b/beamphysics/interfaces/impactx.py index 35f4a0a3..0b7ce43d 100644 --- a/beamphysics/interfaces/impactx.py +++ b/beamphysics/interfaces/impactx.py @@ -29,7 +29,8 @@ is what openPMD's ``position/t + positionOffset/t`` means in ImpactX output, and it keeps quantities like `ParticleGroup.average_current` meaningful. -Originally developed in lume-impactx. +ImpactX source: . This interface was +originally developed in [lume-impactx](https://github.com/ax3l/lume-impactx). """ from __future__ import annotations From 25700863fff2a18769d2166ba534669ff5ef43e9 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 1 Sep 2026 23:06:26 -0700 Subject: [PATCH 3/5] Stay correct once ImpactX stores the lost-particle reference particle BLAST-ImpactX/impactx#1647 hands the beam's live reference particle to the lost-particle container, so particles_lost.* will carry a usable one instead of a default-constructed RefPart. Three consequences here: - The "zeroed reference particle" error now says which ImpactX versions it applies to, rather than presenting it as how ImpactX behaves. - The ref= caveat generalizes. It was written as a cost of *supplying* a reference particle from a monitor; it is really intrinsic. The momenta are normalized by beta_gamma at the reference particle's own s, so wherever the reference energy changes between the loss point and the reference particle in hand, the conversion is approximate. That is the same limitation ImpactX documents for the reference particle it now stores: the invariants (mass, charge) describe lost particles exactly, the kinematic attributes are end-of-tracking values. - The test asserted the fixture has a zeroed reference particle, which would have failed the day someone regenerated it against a newer ImpactX. It now branches on what the file actually holds and asserts the right behaviour for either. The lost-file detection needs no change: it keys on the s_lost record, which CollectLost always adds, and only falls back to the zeroed reference particle. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01DhbdaEDSBqFN974t9QcXem --- beamphysics/interfaces/impactx.py | 28 ++++++++++++++++------------ tests/test_impactx.py | 25 ++++++++++++++++++------- 2 files changed, 34 insertions(+), 19 deletions(-) diff --git a/beamphysics/interfaces/impactx.py b/beamphysics/interfaces/impactx.py index 0b7ce43d..26d545d8 100644 --- a/beamphysics/interfaces/impactx.py +++ b/beamphysics/interfaces/impactx.py @@ -766,10 +766,11 @@ def _load(record_name: str, component: str): raise ValueError( f"{str(path)!r} iteration {iteration} carries a zeroed reference " f"particle (mass_ref={ref.mass_MeV} MeV, gamma_ref={ref.gamma}), so " - "its normalized coordinates cannot be converted. ImpactX writes this " - "for its particles_lost output. Pass ref= with the ImpactXRefPart of " - "the BeamMonitor iteration the particles were lost at, which " - "refpart_from_openpmd() reads from the monitor file." + "its normalized coordinates cannot be converted. ImpactX wrote this " + "into its particles_lost output before BLAST-ImpactX/impactx#1647; " + "newer files carry a usable one. Pass ref= with the ImpactXRefPart of " + "a BeamMonitor iteration, which refpart_from_openpmd() reads from the " + "monitor file." ) if strict: @@ -831,17 +832,20 @@ def read_beam_monitor( to False to drop that data and read the bunch anyway. ref : ImpactXRefPart, optional Reference particle to interpret the coordinates against. Taken from the file - when omitted, which is what you want for a BeamMonitor. ImpactX's - ``particles_lost`` output carries a *zeroed* reference particle, so reading it - requires passing one, e.g. + when omitted, which is what you want for a BeamMonitor. ImpactX wrote a + *zeroed* reference particle into its ``particles_lost`` output before + BLAST-ImpactX/impactx#1647, so reading such a file requires passing one, e.g. ``refpart_from_openpmd(series.iterations[n].particles["beam"])`` from the monitor file. - Be aware that this is an approximation for lost particles: they were lost at - whatever ``s`` the file's own ``s_lost`` record says, not at the monitor's, so - they are un-normalized with the wrong reference momentum and given the wrong - reference time unless the two happen to coincide. It is exact only for - particles lost at the monitor. + Either way, a lost particle is only converted exactly if the reference energy + did not change between where it was lost -- the file's own ``s_lost`` record -- + and the reference particle in hand. The momenta are normalized by ``beta_gamma`` + at the reference particle's own ``s``, so through an RF cavity or any other + accelerating element the two differ, and so does the reference time. This is + the same limitation ImpactX documents for the reference particle it now stores + alongside its lost particles: the invariants (mass, charge) describe them + exactly, the kinematic attributes are end-of-tracking values. Returns ------- diff --git a/tests/test_impactx.py b/tests/test_impactx.py index cbb4fa3c..90573811 100644 --- a/tests/test_impactx.py +++ b/tests/test_impactx.py @@ -538,19 +538,30 @@ def test_species_can_name_what_inference_cannot(bunch): @requires_lost_data -def test_particles_lost_has_no_reference_particle(): - """ImpactX writes a default-constructed RefPart into particles_lost output.""" +def test_particles_lost_reference_particle_either_way(): + """ImpactX wrote a default-constructed RefPart into particles_lost output before + BLAST-ImpactX/impactx#1647; newer versions store a usable one. + + The committed fixture predates that fix, so it takes the first branch. Assert + whichever the file actually holds rather than the version it was made with, so + regenerating the data against a newer ImpactX does not look like a regression. + """ import openpmd_api as io series = io.Series(str(PARTICLES_LOST), io.Access.read_only) ref = refpart_from_openpmd(series.iterations[0].particles["beam"]) series.close() - assert ref.mass_MeV == 0.0 - assert ref.gamma == 0.0 - - with pytest.raises(ValueError, match="zeroed reference particle"): - read_beam_monitor(PARTICLES_LOST) + if ref.mass_MeV == 0.0: + assert ref.gamma == 0.0 + with pytest.raises(ValueError, match="zeroed reference particle"): + read_beam_monitor(PARTICLES_LOST) + else: + # a usable reference particle in the file means no ref= is needed + with pytest.warns(UserWarning, match="s_lost"): + lost = read_beam_monitor(PARTICLES_LOST, strict=False) + assert lost.n_alive == 0 + assert pmd_species_of(ref) == "electron" @requires_lost_data From f1d273acd276ad8ed467cb9b28aa92e938bbd4e4 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 1 Sep 2026 23:13:37 -0700 Subject: [PATCH 4/5] Do not document an unsettled ImpactX design as settled The ref= docstring claimed that ImpactX stores mass and charge as invariants and the kinematic attributes as end-of-tracking values for lost particles. That is the current shape of BLAST-ImpactX/impactx#1647, which is still in flux -- what to key on, and whether final state is even the right thing to store, is undecided. Keep only the part that does not depend on how that lands: a lost particle converts exactly when the reference energy did not change between the loss point and the reference particle in hand, and this reader applies no correction of its own. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01DhbdaEDSBqFN974t9QcXem --- beamphysics/interfaces/impactx.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/beamphysics/interfaces/impactx.py b/beamphysics/interfaces/impactx.py index 26d545d8..fd27c185 100644 --- a/beamphysics/interfaces/impactx.py +++ b/beamphysics/interfaces/impactx.py @@ -842,10 +842,10 @@ def read_beam_monitor( did not change between where it was lost -- the file's own ``s_lost`` record -- and the reference particle in hand. The momenta are normalized by ``beta_gamma`` at the reference particle's own ``s``, so through an RF cavity or any other - accelerating element the two differ, and so does the reference time. This is - the same limitation ImpactX documents for the reference particle it now stores - alongside its lost particles: the invariants (mass, charge) describe them - exactly, the kinematic attributes are end-of-tracking values. + accelerating element the two differ, and so does the reference time. Which + reference particle ImpactX stores alongside lost particles is still settling + upstream; this reader applies no correction of its own -- it uses whatever it + is given -- so pass ``ref=`` explicitly when it matters which one that is. Returns ------- From f50d003335725e87e8502f08ad5551a2c94ea45a Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Wed, 2 Sep 2026 14:48:17 -0700 Subject: [PATCH 5/5] improve docs --- beamphysics/interfaces/impactx.py | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/beamphysics/interfaces/impactx.py b/beamphysics/interfaces/impactx.py index fd27c185..8857c151 100644 --- a/beamphysics/interfaces/impactx.py +++ b/beamphysics/interfaces/impactx.py @@ -1,10 +1,15 @@ """ImpactX beam data <-> ParticleGroup. -[ImpactX](https://impactx.readthedocs.io) is an s-based beam dynamics code, the -successor of IMPACT-Z. Its particles are held at a common ``s`` with a spread in -arrival time, which is z-coordinates on this side -- all ``z`` equal, ``t`` varying -- -so the conversion is a direct algebraic map, like the Bmad interface and unlike the -time-based ones. +[ImpactX](https://impactx.readthedocs.io) enables high-performance modeling of beam +dynamics in particle accelerators with collective effects. This is the next +generation of the IMPACT-Z code. ImpactX runs on modern GPUs or CPUs alike, provides +user-friendly interfaces suitable for AI/ML workflows, has many benchmarks to ensure +its correctness, and an extensive documentation. + +ImpactX models particle beams with respect to a common ``s`` (in this repo called ``z``) +variable for the reference trajectory, with a spread in arrival time: all ``z`` equal, +``t`` varying, so the conversion is a direct algebraic map, like the Bmad interface and +unlike time-integrating codes (e.g., WarpX, ASTRA, etc.). Coordinates and frames ---------------------- @@ -18,6 +23,8 @@ reference momentum: ``px = Delta(beta_x gamma) / (beta_0 gamma_0)`` and ``pt = -Delta(gamma) / (beta_0 gamma_0)``. +See: https://impactx.readthedocs.io/en/latest/theory/coordinates_units.html + `ParticleGroup` is a lab-frame container, so the mapping has to choose a frame: - The transverse coordinates stay in the **local frame**: ``x`` and ``y`` are the @@ -29,8 +36,10 @@ is what openPMD's ``position/t + positionOffset/t`` means in ImpactX output, and it keeps quantities like `ParticleGroup.average_current` meaningful. -ImpactX source: . This interface was -originally developed in [lume-impactx](https://github.com/ax3l/lume-impactx). +See also: + +- ImpactX source: https://github.com/BLAST-ImpactX/impactx +- ImpactX manual: https://impactx.readthedocs.io """ from __future__ import annotations