GUINEA-PIG++ · src/guineapigCPP · src/gridCPP · src/fieldCPP · src/particleBeamCPP

Which slice meets which — and what nt is actually for

GUINEA-PIG chops each bunch into nz longitudinal slices and each encounter into nt sub-steps. They are not two versions of the same knob. Everything below is read off the source in this repository, and the toy collision re-derives it — validated against four analytic limits before it shipped.

Fig. 01

Schedule

Emphasise

two ways to discretise a collision
crossing step k
sub-step i₀
pairs this step
make_step calls so far
beam 1 → (slice 0 = head) ← beam 2 pair being solved now frozen — code does not touch it
crossing schedule i1+i2=k k = 0 … 2nz−2 · guineapigCPP.h:513
meeting point zmeet=(i2i1)2Δz centre of the bin whose lower edge is min_z
push per sub-step step_=Δz2nt=cutznznt gridCPP.cc:146
i1, i2
slice index in beam 1 / beam 2, counted from the head of that bunch
k
crossing step — equivalently, twice the meeting time in units of Δz/2c
i0
sub-step within one k, running 0 … nt−1. The pair list does not change with it
Δz
slice thickness, 2·cutz/nz
Check yourself — do the beams move during the n_t loop?
Set nt = 8 and press Step i₀ repeatedly without changing k. The slabs visibly creep. Yet the pair list, and the min_z bin every luminosity event is written into, do not change. Is the creep real physics or a drawing artefact — and if it is real, why does the bookkeeping ignore it?
It is real, and the bookkeeping ignores it on purpose.

Each sub-step calls advancePosition(step_) on every particle of the two active slices, with step_ = Δz/(2nt). Over the nt sub-steps of one k, each beam therefore advances Δz/2, and since they counter-propagate the relative closing is exactly Δz — one whole slice thickness. That is precisely the displacement that makes the next value of k the correct pairing. The march and the scan are the same motion sampled at two rates.

What the code does not do is re-derive the pairing from that motion. The pair list is recomputed only at the top of the k loop, and min_z = ½(i₂−i₁−1) (guineapigCPP.cc:677) is a function of the indices alone. So within one k the geometry is held piecewise constant while the dynamics runs continuously — a standard operator-split: transport the pairing on the coarse grid, integrate the force on the fine one.

There is a measurable consequence, and the toy in the Convergence tab shows it. With rigid beams the split is exact: the 1/nt in the luminosity estimator (gridCPP.cc:916) cancels the loop and L comes out bit-identical for nt = 1, 2, 4, 8. Give the beams their real angular divergence and it stops being exact — the slabs now change transverse size as they creep, so the nt sub-steps sample σ(z) at nt different points inside one Δz bin instead of one. Measured effect: 0.4 % between nt = 1 and 8. Compare the same quantity's nz dependence over the same range: 10 %. That ratio, 25:1, is the quantitative version of “nz carries the geometry, nt only refines the push”.

The judgement call: whether 0.4 % matters is not something the code can tell you. For a luminosity number it is negligible against the transverse-grid bias. For a z-differential background spectrum, where the Δz binning is already the limiting resolution, it is invisible. It would start to matter only if you were using GUINEA-PIG to study the longitudinal structure of the pinch itself — and at that point you would raise nz, not nt.

Probe

Drag on the map.

Beam

Fill is cost, nz2nt field solves. Contours are constant push length step_ = cutz/(nznt) — lines of roughly equal integration accuracy. They have slope −1 in log–log while cost climbs twice as fast along nz, so the cheapest point on any accuracy contour is its bottom-right end.
cost and push length over the (n_z, n_t) plane
field solves
push length step_
slice Δz
lumi z-bins
what it controlsn_zn_tn_x, n_y
push length step_∝ 1/nz∝ 1/nt
field solves∝ nz2∝ nt
cost of one solve2nx×2ny FFT
particle pushes∝ nz∝ nt
Δz output bin, z-resolution∝ 1/nz
hourglass / pinch geometryresolves it
macroparticles per slice∝ 1/nz
overlap-integral biascells per σ

The third column is the point about nx and ny: the transform is always 2nx × 2ny (fieldCPP.h:40–41) whatever nz and nt are. They set how many transforms, never how big.

Beam

Grid

Physics switches

cumulative L/L₀ through the crossing
L / L₀
σ_y head → tail
outside grid
solves · wall time
L(k)/L₀ and dL/dz beam 1 σ_y, final beam 2 σ_y, final σ_y initial
luminosity estimator L=stepscellsρ1ρ2ΔxΔynt gridCPP.cc:916
The 1/nt is what makes nt luminosity-neutral: each pair is scored nt times and divided by nt. Switch the force off and the readout is bit-identical for nt = 1, 2, 4, 8 — verified to six decimals. Switch divergence off too and it must land on 0.9946, the ±3σz truncation factor erf(3/√2)²; whatever is left is grid bias.
Check yourself — why is L below 1 even with the force off?
Turn the beam-beam force off and leave divergence on. L/L₀ settles near 0.87, not 1. Nothing is deflecting anything. Where did 13 % of the luminosity go?
Into the hourglass, and it is the one thing in this file that tests backSlices().

L₀ = N²/(4πσ*xσ*y) is the luminosity two bunches would have if every pair of slices met exactly at the waist. They do not. Slice i1 meets slice i2 at z = (i2−i1)Δz/2, and at that longitudinal position each beam has expanded to σ(z) = σ*√(1+z²/β*²). For C³-250, σz/β*y = 0.83, so the tails collide with a vertical size several times the waist value.

Summing that over the toy's own slicing gives an analytic prediction with no free parameters: H = ∑ij wiwj / √((1+u²/β*x²)(1+u²/β*y²)) with u the meeting point above — 0.8834 at nz = 8, converging to 0.8868. The toy returns 0.8732 at nz = 8, a ratio of 0.9885; the missing 1.15 % is the −0.64 % transverse-grid bias plus shot noise. That is the validation.

Why it tests backSlices. Nothing in the toy tells a slice how big it should be when it collides. Each slice is generated at the waist and then pre-drifted backwards by distk = (k + ½(1−1/nt))·Δz/2 (particleBeamCPP.h:232), after which it simply drifts forward one step_ at a time. The correct σ(z) at the correct moment is an emergent property of that offset being right. Change the Δz/2 increment to Δz and the hourglass factor comes out wrong; the luminosity closure test above would still pass, because it does not involve v at all.

Note also the nz dependence of H itself: 0.808 at nz = 2, 0.873 at 4, 0.883 at 8, 0.887 in the limit. Two slices cannot represent an hourglass. This is pure geometry — nt does not appear.

Sweep grid

nz ∈ {2,4,8,16} × nt ∈ {1,2,4,8}, sixteen collisions at fixed transverse grid.
L/L₀ against the product n_z·n_t — the thing that sets step_
the same points against what they cost

If the push length were the whole story the four curves would coincide. They do not, and the readout above quantifies the spread for the run you just did. Everything at nz ≥ 8 collapses onto one curve; everything below it is a different physical model. The absolute scale depends on the transverse grid, so raise ny before reading values off the axis — the shape is what this panel is for.

Check yourself — n_z = 4 overshoots. Why?
Run the sweep. The nz = 2 family lands well below the nz = 16 answer, but nz = 4 ends up above it. The error is not monotonic in nz. What are the two effects, and why do they have opposite signs? (Reference values quoted below — 1.29, 1.80 and 1.69 for nz = 2, 4 and 16 — were measured on a 64×128 grid with nt up to 16. The in-browser sweep defaults to a coarser grid and to nt ≤ 8, so its numbers sit a few percent lower.)
Two errors of opposite sign that cross over between nz = 2 and 8.

Downward: the hourglass. Coarse slicing puts too much of the charge at |z| near zero. With nz = 2 the only available meeting points are z = 0 and ±Δz/2 = ±150 µm, and the discrete hourglass factor is 0.808 against 0.887 in the limit — a 9 % deficit sitting directly on the luminosity. This dominates at nz = 2.

Upward: over-concentrated deflectors. A slice's whole charge acts as one transverse sheet. Coarse slicing therefore delivers the pinch in a few large, perfectly-timed kicks instead of many small ones, and because the focusing is nonlinear — a particle kicked inward early sees a denser core for the rest of the crossing — lumping the kicks over-focuses. This is the same error a leapfrog integrator makes with too long a step, except that here the step is a spatial one and nt cannot fix it: nt subdivides the push through one sheet, it does not make the sheet thinner.

At nz = 2 the hourglass deficit (−9 %) beats the over-focusing; at nz = 4 the hourglass deficit has fallen to −1.5 % while the over-focusing is still worth about +8 %, so the sum goes positive; by nz = 8 both are at the percent level and the curve settles. You can watch the first effect alone by switching the beam-beam force off in the Toy collision tab and stepping nz from 2 to 16 — that isolates the geometry with no dynamics at all.

Judgement: the practical reading is that a non-monotonic error curve makes “I doubled nz and the answer barely moved” an unsafe convergence test near nz ≈ 4–8, because you can be sitting near the crossover where the two errors cancel. Converge in nz from above, not from below.

The loop nest, verbatim

// guineapigCPP.h:460  — the crossing schedule, 2n_z−1 steps
for (k=0; k<n_slice; k++)           iteration_on_overlaping_slices(0, k, sor);
for (k=n_slice; k<2*n_slice-1; k++) iteration_on_overlaping_slices(k-n_slice+1, n_slice-1, sor);

// guineapigCPP.cc:819 — n_t sub-steps; the pair list is frozen inside this loop
for (i0=0; i0<grid_.get_timestep(); i0++)
    make_time_step_on_slices(firstSliceOfBeam1, lastSliceOfBeam2, sor);

// guineapigCPP.h:513 — the pairs of one step: the indices run in OPPOSITE directions
slice_beam_2 = lastSliceOfBeam2;
for (slice_beam_1 = firstSliceOfBeam1; slice_beam_1 <= lastSliceOfBeam2; slice_beam_1++) {
    make_step(slice_beam_1, slice_beam_2, sor);   // deposit → FFT solve → interp → kick + drift
    slice_beam_2--;                               // => i1 + i2 = k is invariant
}

Source map

Critique

1 · the n_t loop re-solves the field but not the pairing

Within one k, make_step re-deposits and re-solves nt times while min_z stays fixed. That is a defensible operator split, but it means the cost of raising nt is a full FFT per sub-step per pair, when the pairing — and therefore the source geometry at slice granularity — has not changed. A scheme that froze the field and sub-cycled only the push would get most of the integration accuracy for a fraction of the cost. GUINEA-PIG does not offer that option; nt buys both together whether you need both or not.

2 · cost is quadratic in n_z, but only the linear part is physics

Field solves scale as nz2nt because every slice must meet every other slice. Particle pushes scale only as nzntnm, since each make_step touches just two slices. So at large nz the run is dominated by field solves whose source distributions differ only slightly from one another — the transverse charge profile of slice i and slice i+1 are nearly the same object. Nothing in the code exploits that.

3 · the z resolution of every differential output is set by a parameter chosen for cost

min_z takes only 2nz−1 discrete values and each event is smeared uniformly over one Δz (meshCPP.h:34). Every z-differential distribution GUINEA-PIG writes — pair vertices, hadron vertices, luminosity against z — inherits that binning. nz is normally chosen by how long the job may take, which means an output resolution is being set by a compute budget without that trade being stated anywhere in the parameter file.

4 · the source tree has an encoding trap

src/particleBeamCPP.cc is ISO-8859, not UTF-8. GNU grep classifies it as binary and prints nothing rather than an error, so a search for backstep across src/ comes back empty and the function looks as though it does not exist. grep -a, or an encoding normalisation, avoids an hour of confusion.

Reference values used here