StatMonitor

Introduction

StatMonitor is a beam statistics monitor that records bunch statistics turn-by-turn at a specified longitudinal position, including centroid position, beam size, emittance, Twiss parameters, higher-order moments, and beam loss. It is the core tool for evaluating beam quality evolution and diagnosing beam dynamics behavior.

Code location: PASS/commands/monitor/statistic.py
Class name: StatMonitor, registered name "statmonitor"
Key features:
- Computes bunch statistics in 6D phase space turn-by-turn (first through fourth order moments);
- Derives emittance and Twiss parameters (beta, alpha, gamma) from second-order moments;
- Records beam loss count and loss percentage;
- CPU uses numpy vectorized computation, GPU uses CUDA kernel functions + warp reduction;
- Appends data to CSV each turn, converts to TFS format uniformly on the final turn;
- Only surviving particles (tag > 0) are counted; lost particles are excluded.

Working Principle

Statistics Computation

For \(N\) surviving particles in the bunch (\(\text{tag} > 0\)), the moments of each order are defined as:

First-order moment (centroid):

\[\langle x \rangle = \frac{1}{N} \sum_{i=1}^{N} x_i\]

Second-order moment:

\[\langle x^2 \rangle = \frac{1}{N} \sum_{i=1}^{N} x_i^2\]

Covariance:

\[\langle x \, p_x \rangle = \frac{1}{N} \sum_{i=1}^{N} x_i \, p_{x,i}\]

Beam size (standard deviation):

\[\sigma_x = \sqrt{\langle x^2 \rangle - \langle x \rangle^2}\]

Similarly, \(\sigma_{p_x}\), \(\sigma_y\), \(\sigma_{p_y}\), \(\sigma_z\), \(\sigma_{\delta}\) are computed.

Longitudinal statistics use the bunch-relative coordinate \(z_{\mathrm{rel}}\). Before calculating moments, the program folds a temporary copy into \([-C/2,C/2)\) to prevent equivalent ring-period representations from inflating the mean and variance. This statistical fold is not written back to the tracked particle coordinates. Under the ring-coordinate convention, the laboratory centroid can be reconstructed from \(z_{\mathrm{center}}+\langle z_{\mathrm{rel}}\rangle\).

Emittance and Twiss Parameters

The 2D emittance is derived from second-order moments:

\[\varepsilon_x = \sqrt{\sigma_x^2 \, \sigma_{p_x}^2 - \sigma_{x,p_x}^2}\]

where \(\sigma_{x,p_x} = \langle x \, p_x \rangle - \langle x \rangle \langle p_x \rangle\) is the covariance.

Twiss parameters:

\[\beta_x = \frac{\sigma_x^2}{\varepsilon_x}\]

\[\alpha_x = -\frac{\sigma_{x,p_x}}{\varepsilon_x}\]

\[\gamma_x = \frac{\sigma_{p_x}^2}{\varepsilon_x}\]

Invariant verification:

\[\gamma_x \, \beta_x - \alpha_x^2 = 1\]

The formulas for the vertical direction (y) are identical in form; simply replace the subscript x with y.

Higher-Order Moments

Skewness (third standardized moment):

\[S_x = \frac{\langle x^3 \rangle - 3 \langle x \rangle \sigma_x^2 - \langle x \rangle^3}{\sigma_x^3}\]

Kurtosis (fourth standardized moment):

\[K_x = \frac{\langle x^4 \rangle - 4 \langle x \rangle \langle x^3 \rangle + 2 \langle x \rangle^2 \langle x^2 \rangle + 4 \langle x \rangle^2 \sigma_x^2 + \langle x \rangle^4}{\sigma_x^4}\]

Beam Loss

\[N_{\text{loss}} = N_{\text{total}} - N_{\text{alive}}\]

\[\text{loss\%} = \frac{N_{\text{loss}}}{N_{\text{total}}} \times 100\%\]

where \(N_{\text{total}}\) is the initial number of macro particles in the bunch, and \(N_{\text{alive}}\) is the current number of surviving particles.

GPU Implementation

The GPU version uses the CUDA kernel function calc_all_stats, employing a grid stride loop to traverse particles. Each thread accumulates 22 statistics in registers, then through warp reduction (__shfl_down_sync) and block reduction, writes the global result via atomicAdd. The maximum number of blocks is 512 (due to atomicAdd contention overhead).

Interface Parameters

Property	JSON key	Type	Default	Description
`s`	`"S (m)"`	float	Required	Longitudinal position of the monitor in the beamline
`cmd_name`	`"name"`	str	Required	Monitor name (automatically filled from the sequence key name)
`command`	`"Command"`	str	`"StatMonitor"`	Command type identifier

Note

StatMonitor has no additional configuration parameters. The statistics target all surviving particles in the bunch at that position (tag > 0); no particle indices need to be specified.

Output Files

A pair of files is generated for each bunch at each monitor position:

CSV (appended turn-by-turn): {hms}_stat_beam{bid}_bunch{bid}_Np_{Np}_s_{s:.4f}.csv
TFS (converted from CSV on the final turn): {hms}_stat_beam{bid}_bunch{bid}_Np_{Np}_s_{s:.4f}.tfs

The output directory is output_dir_stat.

TFS file header:

@ Name             PASS Statistic Data
@ Time             2026-07-14 00:11:03

Output columns (35 columns total):

Column name	Group	Description
`turn`	Basic	Turn number
`xAverage`	Centroid	Horizontal position mean \(\langle x \rangle\)
`pxAverage`	Centroid	Horizontal momentum mean \(\langle p_x \rangle\)
`sigmaX`	Beam size	Horizontal position standard deviation \(\sigma_x\)
`sigmaPx`	Beam size	Horizontal momentum standard deviation \(\sigma_{p_x}\)
`yAverage`	Centroid	Vertical position mean
`pyAverage`	Centroid	Vertical momentum mean
`sigmaY`	Beam size	Vertical position standard deviation
`sigmaPy`	Beam size	Vertical momentum standard deviation
`zAverage`	Centroid	Mean of the folded bunch-relative coordinate \(\langle z_{\mathrm{rel}}\rangle\)
`dpAverage`	Centroid	Momentum deviation mean
`sigmaZ`	Beam size	Standard deviation of the folded bunch-relative coordinate
`sigmadp`	Beam size	Momentum deviation standard deviation
`xEmittance`	Emittance	Horizontal 2D emittance \(\varepsilon_x\)
`yEmittance`	Emittance	Vertical 2D emittance \(\varepsilon_y\)
`betax`	Twiss	Horizontal beta function
`betay`	Twiss	Vertical beta function
`alphax`	Twiss	Horizontal alpha function
`alphay`	Twiss	Vertical alpha function
`gammax`	Twiss	Horizontal gamma function
`gammay`	Twiss	Vertical gamma function
`invariantx`	Verification	Horizontal invariant \(\gamma_x \beta_x - \alpha_x^2\) (should equal 1)
`invarianty`	Verification	Vertical invariant (should equal 1)
`zCenter`	Longitudinal reference	Laboratory longitudinal center of the bunch, \(z_{\mathrm{center}}\)
`xzAverage`	Correlation	\(\langle x \, z \rangle\)
`xyAverage`	Correlation	\(\langle x \, y \rangle\)
`yzAverage`	Correlation	\(\langle y \, z \rangle\)
`xzDevideSigmaxSigmaz`	Correlation	\(\langle x \, z \rangle / (\sigma_x \, \sigma_z)\) normalized correlation
`beamLossTotal`	Loss	Number of lost particles
`lossPercent`	Loss	Loss percentage
`xSkewness`	Higher-order moments	Horizontal skewness
`xKurtosis`	Higher-order moments	Horizontal kurtosis
`ySkewness`	Higher-order moments	Vertical skewness
`yKurtosis`	Higher-order moments	Vertical kurtosis
`Ek`	Energy	Bunch kinetic energy

Usage Example

The following JSON snippet places a statistics monitor at \(s = 0.0\) m:

"SM1": {
    "S (m)": 0.0,
    "Command": "StatMonitor"
}

The statistics monitor requires no additional parameters; only the position and command type need to be specified. During simulation, the bunch statistics at that position are recorded turn-by-turn.

Multi-position Monitoring

Multiple statistics monitors can be placed at different positions to compare the variation of bunch statistics along the beamline:

"SM_start": {
    "S (m)": 0.0,
    "Command": "StatMonitor"
},
"SM_mid": {
    "S (m)": 100.0,
    "Command": "StatMonitor"
},
"SM_end": {
    "S (m)": 250.0,
    "Command": "StatMonitor"
}

Application Scenarios

Beam quality assessment: Monitor the evolution of emittance, beam size, and centroid position turn-by-turn to evaluate whether beam quality is stable or degrading
Emittance measurement: Compute emittance and Twiss parameters from second-order moments, and compare with design values for verification
Beam loss diagnostics: Monitor beam loss rate through beamLossTotal and lossPercent, identifying the turn and position where losses occur
Nonlinear effect identification: Use higher-order moment information (skewness and kurtosis) to determine the degree to which the beam distribution deviates from Gaussian, identifying nonlinear resonances or dispersion coupling
Momentum-spread monitoring: sigmadp and sigmaZ describe longitudinal beam quality within a bunch; use zCenter as well when comparing absolute azimuths across bunches
Correlation diagnostics: Correlation quantities such as xzAverage can be used to diagnose dispersion coupling or transverse-longitudinal coupling