Read examples¶
from beamphysics import ParticleGroup
elegant hdf5 format¶
from beamphysics.interfaces.elegant import elegant_h5_to_data
This will convert to a data dict
data = elegant_h5_to_data("data/elegant_raw.h5")
data
{'x': array([-1.12390181e-04, -7.20268439e-05, -1.14543953e-04, ...,
1.05737263e-04, -7.24786314e-05, 5.78264247e-05], shape=(50000,)),
'y': array([-1.23229455e-04, -1.08820261e-04, -9.27723036e-05, ...,
8.55980575e-05, 8.58175111e-05, 9.07622257e-05], shape=(50000,)),
'z': array([0, 0, 0, ..., 0, 0, 0], shape=(50000,)),
'px': array([ -3083.62426208, 2225.50916089, 2725.71118497, ...,
-15364.10019502, 15821.42535919, -4733.37387114], shape=(50000,)),
'py': array([ 83545.59793901, 83389.57159329, 68275.02707993, ...,
-70430.7890369 , -70600.14629199, -69547.74384872], shape=(50000,)),
'pz': array([8.00396788e+09, 8.00007259e+09, 8.00139351e+09, ...,
7.99685995e+09, 7.99660959e+09, 7.99818066e+09], shape=(50000,)),
't': array([5.11804566e-06, 5.11804568e-06, 5.11804567e-06, ...,
5.11804569e-06, 5.11804569e-06, 5.11804569e-06], shape=(50000,)),
'status': array([1, 1, 1, ..., 1, 1, 1], shape=(50000,)),
'species': 'electron',
'weight': array([2.e-17, 2.e-17, 2.e-17, ..., 2.e-17, 2.e-17, 2.e-17],
shape=(50000,)),
'id': array([ 1, 2, 3, ..., 49998, 49999, 50000],
shape=(50000,), dtype=int32)}
Create ParticleGroup and plot:
P = ParticleGroup(data=data)
P.plot("delta_t", "delta_pz")
Genesis4 HDF5 format¶
Genesis4 stores particles in z slices in its HDF5 file format. Each slice has a z length equal to the reference wavelength lambda0. In order to speed up the simulation, Genesis4 also has an integer parameter sample that essentially skips adjacent slices, defined in the input file:
# Input
!cat data/genesis4/genesis4.in
# Lattice
!cat data/genesis4/genesis4.lat
&setup rootname = genesis4 lattice = genesis4.lat beamline = LAT gamma0 = 1000.0 lambda0 = 1e-07 delz = 0.026 seed = 123456 npart = 128 &end &time slen = 7.175993210500245e-05 sample = 100 &end # A negligible-power field is declared so Genesis4 sets the radiation # wavenumber xks = 2*pi/lambda0 in the longitudinal ODE. Without a field, # xks defaults to 1 and the energy-dependent slippage (velocity bunching / # R56) is switched off. In a pure drift (aw=0) the field is dynamically # inert -- it exchanges no energy with the beam -- so this only activates # the correct kinematic longitudinal physics. &field power = 1e-6 dgrid = 1e-3 ngrid = 51 &end # 100 pC &profile_gauss label = beamcurrent c0 = 1000.0 s0 = 3.5879966052501225e-05 sig = 1.1959988684167075e-05 &end &beam gamma = 1000.0 delgam = 1.0 current = @beamcurrent &end &write beam=beginning &end &track zstop = 1.0 &end &write beam=end &end
D1: drift = {l=1.0};
LAT: LINE = {D1};
Running Genesis4 on this input (genesis4 genesis4.in) produces the raw particle file end.par.h5, which contains the particle data for each slice. This file and the one4one.par.h5 file used later in this notebook are small and are committed to the repository. The data directory provides a run.sh script that regenerates them by running Genesis4 on genesis4.in and the one4one input; run it if you have the genesis4 executable and want to reproduce the data.
parfile = "data/genesis4/end.par.h5"
By default ParticleGroup.from_genesis4 will load these particles relative to their slice position. Because sample=100 we see a comb-like distribution:
P_default = ParticleGroup.from_genesis4(parfile)
P_default.plot("z", "pz", bins=200)
P_default
<ParticleGroup with 896 particles at 0x7fd114b28910>
ParticleGroup.from_genesis4 has a smear option that shifts the particles in each slice by a random integer multiple of the slice spacing, so that the overall distribution is smoother while the bunching characteristics are preserved. The smearing uses random numbers; pass an integer seed or a numpy.random.Generator as the rng argument to make the result reproducible.
P_smear = ParticleGroup.from_genesis4(parfile, smear=True)
P_smear.plot("z", "pz", bins=100)
P_smear
<ParticleGroup with 896 particles at 0x7fd11bb30b90>
Looking closer:
P = P_default
P[(P.z > 30e-6) & (P.z < 35e-6)].plot("z", "pz", bins=100)
P
<ParticleGroup with 896 particles at 0x7fd114b28910>
P = P_smear
P[(P.z > 30e-6) & (P.z < 35e-6)].plot("z", "pz", bins=100)
P
<ParticleGroup with 896 particles at 0x7fd11bb30b90>
Normally Genesis4 populates each slice with the same number of particles, and each slice carries its own weight (internally, the current). In some cases it is simpler to work with particles that all have the same weight, so ParticleGroup.from_genesis4 also has an equal_weights option that resamples the particles within each slice according to the relative slice weights. Only existing particles are used — the beam is never upsampled — so there are fewer particles in this case, but the projected densities are roughly the same as above.
P_equal = ParticleGroup.from_genesis4(parfile, smear=True, equal_weights=True)
P_equal.plot("z", "pz", bins=100)
P_equal
<ParticleGroup with 402 particles at 0x7fd114dd8290>
Alternatively, the beam can be thinned to an exact number of particles with the n_particle option. Unlike equal_weights, this retains the per-slice weights, which are rescaled once at the end so that the total charge is conserved. Each slice is thinned in proportion to its particle count, so the longitudinal charge profile is preserved. Because the two options control the output in different ways, n_particle and equal_weights are mutually exclusive, and requesting both raises a ValueError.
Here the beam is thinned from its 896 particles to exactly 500:
P_thin = ParticleGroup.from_genesis4(parfile, smear=True, n_particle=500, rng=0)
P_thin.plot("z", "pz", bins=100)
P_thin
<ParticleGroup with 500 particles at 0x7fd11bb908d0>
one4one mode¶
Genesis4 also has a one4one mode, in which each macroparticle represents a single real electron rather than a quiet-loaded sample. The number of particles in a slice is then the integer electron count round(ne), so low-current slices in the tails of the distribution can contain very few — or even zero — particles. ParticleGroup.from_genesis4 handles these variable-size and empty slices automatically:
P = ParticleGroup.from_genesis4("data/genesis4/one4one.par.h5")
set(P.weight)
{np.float64(1.602176634e-19)}
In one4one mode every macroparticle is a single real electron, so ParticleGroup.from_genesis4 assigns each one a weight of exactly the elementary charge qe — set(P.weight) is a single value.
This is worth a note because of how Genesis4 stores a one4one beam internally: each slice keeps its smooth requested current but is populated with round(ne) integer particles, where ne = current · slicespacing / (c · qe). Reconstructing the weight from the current as q1 = ne·qe/round(ne) would scatter it around qe (most in the low-count tail slices). Since each macroparticle physically represents one electron, the reader instead assigns qe directly. This is exactly what Genesis4's &sort namelist does — it resets each slice current to npart · qe — so a sorted file reconstructs to the same uniform weights.
Because the per-particle weight is fixed at qe, the cutoff argument (which drops sub-electron macroparticles for quiet-loaded beams) does not apply to one4one files and the full beam is preserved — here all 935 macroparticles, totalling 935 electrons.
Subsampling on read¶
The n_particle option shown above is especially useful for one4one files, which can be very large. The downsampling happens immediately, as the file is read: the particles to keep are chosen up front from cheap HDF5 metadata and are apportioned across the slices in proportion to their sizes, and only the selected particles are ever read from the file, so the full beam is never loaded into memory.
To decide how many particles to keep, genesis4_parfile_n_particle cheaply reports the file's total count by reading only the dataset shapes, with no bulk particle data:
from beamphysics.interfaces.genesis import genesis4_parfile_n_particle
genesis4_parfile_n_particle("data/genesis4/one4one.par.h5")
935
Here the beam is thinned from its 935 particles to exactly 200, and the total charge is conserved:
P_full = ParticleGroup.from_genesis4("data/genesis4/one4one.par.h5")
P_sub = ParticleGroup.from_genesis4(
"data/genesis4/one4one.par.h5", n_particle=200, rng=0
)
P_sub.n_particle, P_sub.charge, P_full.charge
(200, np.float64(1.4980351527899998e-16), np.float64(1.4980351527899996e-16))
Requesting n_particle greater than or equal to the file's total count is a no-op, and all of the particles are returned:
ParticleGroup.from_genesis4("data/genesis4/one4one.par.h5", n_particle=9350)
<ParticleGroup with 935 particles at 0x7fd114ee7d40>