diff --git a/docs/source/usage/python.rst b/docs/source/usage/python.rst
index 9954002173..e3bec183da 100644
--- a/docs/source/usage/python.rst
+++ b/docs/source/usage/python.rst
@@ -167,9 +167,9 @@ Particles
 
    This class stores particles, distributed over MPI ranks.
 
-   .. py:method:: add_n_particles(lev, x, y, z, px, py, pyz, qm, bchchg)
+   .. py:method:: add_n_particles(lev, x, y, t, px, py, pt, qm, bchchg)
 
-      Add new particles to the container
+      Add new particles to the container for fixed s.
 
       Note: This can only be used *after* the initialization (grids) have
             been created, meaning after the call to :py:meth:`ImpactX.init_grids`
@@ -178,10 +178,10 @@ Particles
       :param lev: mesh-refinement level
       :param x: positions in x
       :param y: positions in y
-      :param z: positions in z
+      :param t: positions as time-of-flight in c*t
       :param px: momentum in x
       :param py: momentum in y
-      :param pz: momentum in z
+      :param pt: momentum in t
       :param qm: charge over mass in 1/eV
       :param bchchg: total charge within a bunch in C
 
diff --git a/examples/cfchannel/analysis_cfchannel.py b/examples/cfchannel/analysis_cfchannel.py
index 8ef3d15f42..4e928ff877 100755
--- a/examples/cfchannel/analysis_cfchannel.py
+++ b/examples/cfchannel/analysis_cfchannel.py
@@ -23,7 +23,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -34,7 +34,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/cfchannel/analysis_cfchannel_10nC.py b/examples/cfchannel/analysis_cfchannel_10nC.py
index e383cd972a..df7c05908e 100755
--- a/examples/cfchannel/analysis_cfchannel_10nC.py
+++ b/examples/cfchannel/analysis_cfchannel_10nC.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/chicane/analysis_chicane.py b/examples/chicane/analysis_chicane.py
index d094519ace..07ceb5cdd3 100755
--- a/examples/chicane/analysis_chicane.py
+++ b/examples/chicane/analysis_chicane.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/chicane/plot_chicane.py b/examples/chicane/plot_chicane.py
index f8505fed19..1f87e4ae97 100755
--- a/examples/chicane/plot_chicane.py
+++ b/examples/chicane/plot_chicane.py
@@ -28,7 +28,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -39,7 +39,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
@@ -178,7 +178,7 @@ def read_time_series(file_pattern):
     ax = axs[(0, ncol_ax)]
     beam_at_step = series.iterations[step].particles["beam"].to_df()
     ax.scatter(
-        beam_at_step.position_ct.multiply(millimeter),
+        beam_at_step.position_t.multiply(millimeter),
         beam_at_step.momentum_t.multiply(mrad),
         s=0.01,
     )
@@ -202,7 +202,7 @@ def read_time_series(file_pattern):
     ax = axs[(2, ncol_ax)]
     beam_at_step = series.iterations[step].particles["beam"].to_df()
     ax.scatter(
-        beam_at_step.position_ct.multiply(millimeter),
+        beam_at_step.position_t.multiply(millimeter),
         beam_at_step.position_x.multiply(millimeter),
         s=0.01,
     )
@@ -212,7 +212,7 @@ def read_time_series(file_pattern):
     ax = axs[(3, ncol_ax)]
     beam_at_step = series.iterations[step].particles["beam"].to_df()
     ax.scatter(
-        beam_at_step.position_ct.multiply(millimeter),
+        beam_at_step.position_t.multiply(millimeter),
         beam_at_step.momentum_x.multiply(mrad),
         s=0.01,
     )
diff --git a/examples/distgen/analysis_gaussian.py b/examples/distgen/analysis_gaussian.py
index 70250b3525..c2c582b7df 100755
--- a/examples/distgen/analysis_gaussian.py
+++ b/examples/distgen/analysis_gaussian.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/distgen/analysis_kurth4d.py b/examples/distgen/analysis_kurth4d.py
index 70777b3bcb..40851889bc 100755
--- a/examples/distgen/analysis_kurth4d.py
+++ b/examples/distgen/analysis_kurth4d.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/distgen/analysis_kvdist.py b/examples/distgen/analysis_kvdist.py
index 4b967161b7..2399af61c9 100755
--- a/examples/distgen/analysis_kvdist.py
+++ b/examples/distgen/analysis_kvdist.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/distgen/analysis_semigaussian.py b/examples/distgen/analysis_semigaussian.py
index 70250b3525..c2c582b7df 100755
--- a/examples/distgen/analysis_semigaussian.py
+++ b/examples/distgen/analysis_semigaussian.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/expanding_beam/analysis_expanding.py b/examples/expanding_beam/analysis_expanding.py
index 6577d874b4..cdf5ac0cf6 100755
--- a/examples/expanding_beam/analysis_expanding.py
+++ b/examples/expanding_beam/analysis_expanding.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/fodo/analysis_fodo.py b/examples/fodo/analysis_fodo.py
index 7294bf1f70..d4c4fc29c5 100755
--- a/examples/fodo/analysis_fodo.py
+++ b/examples/fodo/analysis_fodo.py
@@ -23,7 +23,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -34,7 +34,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/fodo/plot_fodo.py b/examples/fodo/plot_fodo.py
index 925da2b1b6..8fa9b8f585 100755
--- a/examples/fodo/plot_fodo.py
+++ b/examples/fodo/plot_fodo.py
@@ -28,7 +28,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -39,7 +39,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/fodo_rf/analysis_fodo_rf.py b/examples/fodo_rf/analysis_fodo_rf.py
index bf6c3a29fe..32ed5f6eed 100755
--- a/examples/fodo_rf/analysis_fodo_rf.py
+++ b/examples/fodo_rf/analysis_fodo_rf.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/iota_lattice/analysis_iotalattice.py b/examples/iota_lattice/analysis_iotalattice.py
index dd184dba1c..376033bf78 100755
--- a/examples/iota_lattice/analysis_iotalattice.py
+++ b/examples/iota_lattice/analysis_iotalattice.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/kurth/analysis_kurth.py b/examples/kurth/analysis_kurth.py
index 30fc0e6b2f..6a59b829ff 100755
--- a/examples/kurth/analysis_kurth.py
+++ b/examples/kurth/analysis_kurth.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/kurth/analysis_kurth_10nC.py b/examples/kurth/analysis_kurth_10nC.py
index 7feb4d9bb3..28622ec4db 100755
--- a/examples/kurth/analysis_kurth_10nC.py
+++ b/examples/kurth/analysis_kurth_10nC.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/kurth/analysis_kurth_10nC_periodic.py b/examples/kurth/analysis_kurth_10nC_periodic.py
index 6c57f4f82e..e1c599ccbe 100755
--- a/examples/kurth/analysis_kurth_10nC_periodic.py
+++ b/examples/kurth/analysis_kurth_10nC_periodic.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/kurth/analysis_kurth_periodic.py b/examples/kurth/analysis_kurth_periodic.py
index 2a11b5f267..88a7b7a5c8 100755
--- a/examples/kurth/analysis_kurth_periodic.py
+++ b/examples/kurth/analysis_kurth_periodic.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/multipole/analysis_multipole.py b/examples/multipole/analysis_multipole.py
index 4b0e27103e..4aa2358b63 100755
--- a/examples/multipole/analysis_multipole.py
+++ b/examples/multipole/analysis_multipole.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/quadrupole_softedge/analysis_quadrupole_softedge.py b/examples/quadrupole_softedge/analysis_quadrupole_softedge.py
index 70250b3525..c2c582b7df 100755
--- a/examples/quadrupole_softedge/analysis_quadrupole_softedge.py
+++ b/examples/quadrupole_softedge/analysis_quadrupole_softedge.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/rfcavity/analysis_rfcavity.py b/examples/rfcavity/analysis_rfcavity.py
index 3b371a9e2f..51a07d4d89 100755
--- a/examples/rfcavity/analysis_rfcavity.py
+++ b/examples/rfcavity/analysis_rfcavity.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/rotation/analysis_rotation.py b/examples/rotation/analysis_rotation.py
index b69341b132..d4f88ce0e4 100755
--- a/examples/rotation/analysis_rotation.py
+++ b/examples/rotation/analysis_rotation.py
@@ -23,7 +23,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -34,7 +34,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/solenoid/analysis_solenoid.py b/examples/solenoid/analysis_solenoid.py
index bb109a16f0..2a686f030f 100755
--- a/examples/solenoid/analysis_solenoid.py
+++ b/examples/solenoid/analysis_solenoid.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/examples/solenoid_softedge/analysis_solenoid_softedge.py b/examples/solenoid_softedge/analysis_solenoid_softedge.py
index ed98710892..2e99c02982 100755
--- a/examples/solenoid_softedge/analysis_solenoid_softedge.py
+++ b/examples/solenoid_softedge/analysis_solenoid_softedge.py
@@ -22,7 +22,7 @@ def get_moments(beam):
     sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
     sigy = moment(beam["position_y"], moment=2) ** 0.5
     sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
-    sigt = moment(beam["position_ct"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
     sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
 
     epstrms = beam.cov(ddof=0)
@@ -33,7 +33,7 @@ def get_moments(beam):
         sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2
     ) ** 0.5
     emittance_t = (
-        sigt**2 * sigpt**2 - epstrms["position_ct"]["momentum_t"] ** 2
+        sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2
     ) ** 0.5
 
     return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
diff --git a/src/particles/ImpactXParticleContainer.H b/src/particles/ImpactXParticleContainer.H
index 40eaa2938d..dd83041154 100644
--- a/src/particles/ImpactXParticleContainer.H
+++ b/src/particles/ImpactXParticleContainer.H
@@ -39,18 +39,23 @@ namespace impactx
     {
         enum
         {
-            x,  ///< position in x [m] (at fixed t OR fixed s)
-            y,  ///< position in y [m] (at fixed t OR fixed s)
-            z,  ///< position in z [m] (at fixed t) OR time-of-flight ct [m] (at fixed s)
+            x,  ///< position in x [m] (at fixed s OR fixed t)
+            y,  ///< position in y [m] (at fixed s OR fixed t)
+            t,  ///< time-of-flight c*t [m] (at fixed s)
             nattribs ///< the number of attributes above (always last)
         };
 
+        // at fixed t, the third component represents the position z
+        enum {
+            z = t  ///< position in z [m] (at fixed t)
+        };
+
+        //! named labels for fixed s
+        static constexpr auto names_s = { "position_x", "position_y", "position_t" };
         //! named labels for fixed t
         static constexpr auto names_t = { "position_x", "position_y", "position_z" };
-        //! named labels for fixed s
-        static constexpr auto names_s = { "position_x", "position_y", "position_ct" };
-        static_assert(names_t.size() == nattribs);
         static_assert(names_s.size() == nattribs);
+        static_assert(names_t.size() == nattribs);
     };
 
     /** This struct indexes the additional Real attributes
@@ -60,20 +65,25 @@ namespace impactx
     {
         enum
         {
-            ux,  ///< momentum in x, scaled by the magnitude of the reference momentum [unitless] (at fixed t or s)
-            uy,  ///< momentum in y, scaled by the magnitude of the reference momentum [unitless] (at fixed t or s)
-            pt,  ///< momentum in z, scaled by the magnitude of the reference momentum [unitless] (at fixed t) OR energy deviation, scaled by speed of light * the magnitude of the reference momentum [unitless] (at fixed s)
-            m_qm, ///< charge to mass ratio, in q_e/m_e (q_e/eV) TODO: rename to qm_m
-            w,   ///< particle weight, unitless
+            ux,  ///< momentum in x, scaled by the magnitude of the reference momentum [unitless] (at fixed s or t)
+            uy,  ///< momentum in y, scaled by the magnitude of the reference momentum [unitless] (at fixed s or t)
+            pt,  ///< energy deviation, scaled by speed of light * the magnitude of the reference momentum [unitless] (at fixed s)
+            qm,  ///< charge to mass ratio, in q_e/m_e [q_e/eV]
+            w,   ///< particle weight, number of real particles represented by this macroparticle [unitless]
             nattribs ///< the number of attributes above (always last)
         };
 
-        //! named labels for fixed t
-        static constexpr auto names_t = { "momentum_x", "momentum_y", "momentum_z", "qmm", "weighting" };
+        // at fixed t, the third component represents the momentum in z
+        enum {
+            pz = pt  ///< momentum in z, scaled by the magnitude of the reference momentum [unitless] (at fixed t)
+        };
+
         //! named labels for fixed s
-        static constexpr auto names_s = { "momentum_x", "momentum_y", "momentum_t", "qmm", "weighting" };
-        static_assert(names_t.size() == nattribs);
+        static constexpr auto names_s = { "momentum_x", "momentum_y", "momentum_t", "qm", "weighting" };
+        //! named labels for fixed t
+        static constexpr auto names_t = { "momentum_x", "momentum_y", "momentum_z", "qm", "weighting" };
         static_assert(names_s.size() == nattribs);
+        static_assert(names_t.size() == nattribs);
     };
 
     /** This struct indexes the additional Integer attributes
@@ -139,7 +149,7 @@ namespace impactx
         //! Destruct a particle container
         virtual ~ImpactXParticleContainer() = default;
 
-        /** Add new particles to the container
+        /** Add new particles to the container for fixed s.
          *
          * Note: This can only be used *after* the initialization (grids) have
          *       been created, meaning after the call to AmrCore::InitFromScratch
@@ -149,10 +159,10 @@ namespace impactx
          * @param lev mesh-refinement level
          * @param x positions in x
          * @param y positions in y
-         * @param z positions in z
+         * @param t positions as time-of-flight in c*t
          * @param px momentum in x
          * @param py momentum in y
-         * @param pz momentum in z
+         * @param pt momentum in t
          * @param qm charge over mass in 1/eV
          * @param bchchg total charge within a bunch in C
          */
@@ -160,10 +170,10 @@ namespace impactx
         AddNParticles (int lev,
                        amrex::Vector<amrex::ParticleReal> const & x,
                        amrex::Vector<amrex::ParticleReal> const & y,
-                       amrex::Vector<amrex::ParticleReal> const & z,
+                       amrex::Vector<amrex::ParticleReal> const & t,
                        amrex::Vector<amrex::ParticleReal> const & px,
                        amrex::Vector<amrex::ParticleReal> const & py,
-                       amrex::Vector<amrex::ParticleReal> const & pz,
+                       amrex::Vector<amrex::ParticleReal> const & pt,
                        amrex::ParticleReal const & qm,
                        amrex::ParticleReal const & bchchg);
 
diff --git a/src/particles/ImpactXParticleContainer.cpp b/src/particles/ImpactXParticleContainer.cpp
index 6056183d6f..165b1295b0 100644
--- a/src/particles/ImpactXParticleContainer.cpp
+++ b/src/particles/ImpactXParticleContainer.cpp
@@ -81,10 +81,10 @@ namespace impactx
     ImpactXParticleContainer::AddNParticles (int lev,
                                              amrex::Vector<amrex::ParticleReal> const & x,
                                              amrex::Vector<amrex::ParticleReal> const & y,
-                                             amrex::Vector<amrex::ParticleReal> const & z,
+                                             amrex::Vector<amrex::ParticleReal> const & t,
                                              amrex::Vector<amrex::ParticleReal> const & px,
                                              amrex::Vector<amrex::ParticleReal> const & py,
-                                             amrex::Vector<amrex::ParticleReal> const & pz,
+                                             amrex::Vector<amrex::ParticleReal> const & pt,
                                              amrex::ParticleReal const & qm,
                                              amrex::ParticleReal const & bchchg)
     {
@@ -92,10 +92,10 @@ namespace impactx
 
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lev == 0, "AddNParticles: only lev=0 is supported yet.");
         AMREX_ALWAYS_ASSERT(x.size() == y.size());
-        AMREX_ALWAYS_ASSERT(x.size() == z.size());
+        AMREX_ALWAYS_ASSERT(x.size() == t.size());
         AMREX_ALWAYS_ASSERT(x.size() == px.size());
         AMREX_ALWAYS_ASSERT(x.size() == py.size());
-        AMREX_ALWAYS_ASSERT(x.size() == pz.size());
+        AMREX_ALWAYS_ASSERT(x.size() == pt.size());
 
         // number of particles to add
         int const np = x.size();
@@ -124,9 +124,9 @@ namespace impactx
             ParticleType p;
             p.id() = ParticleType::NextID();
             p.cpu() = amrex::ParallelDescriptor::MyProc();
-            p.pos(0) = x[i];
-            p.pos(1) = y[i];
-            p.pos(2) = z[i];
+            p.pos(RealAoS::x) = x[i];
+            p.pos(RealAoS::y) = y[i];
+            p.pos(RealAoS::t) = t[i];
             // write position, creating cpu id, and particle id
             pinned_tile.push_back(p);
         }
@@ -136,8 +136,8 @@ namespace impactx
 
         pinned_tile.push_back_real(RealSoA::ux, px);
         pinned_tile.push_back_real(RealSoA::uy, py);
-        pinned_tile.push_back_real(RealSoA::pt, pz);
-        pinned_tile.push_back_real(RealSoA::m_qm, np, qm);
+        pinned_tile.push_back_real(RealSoA::pt, pt);
+        pinned_tile.push_back_real(RealSoA::qm, np, qm);
         pinned_tile.push_back_real(RealSoA::w, np, bchchg/ablastr::constant::SI::q_e/np);
 
         /* Redistributes particles to their respective tiles (spatial bucket
diff --git a/src/particles/ReferenceParticle.H b/src/particles/ReferenceParticle.H
index dede525de8..e4ed21b7db 100644
--- a/src/particles/ReferenceParticle.H
+++ b/src/particles/ReferenceParticle.H
@@ -31,7 +31,7 @@ namespace impactx
         amrex::ParticleReal s = 0.0;  ///< integrated orbit path length, in meters
         amrex::ParticleReal x = 0.0;  ///< horizontal position x, in meters
         amrex::ParticleReal y = 0.0;  ///< vertical position y, in meters
-        amrex::ParticleReal z = 0.0;  ///< longitudinal position y, in meters
+        amrex::ParticleReal z = 0.0;  ///< longitudinal position z, in meters
         amrex::ParticleReal t = 0.0;  ///< clock time * c in meters
         amrex::ParticleReal px = 0.0; ///< momentum in x, normalized to proper velocity
         amrex::ParticleReal py = 0.0; ///< momentum in y, normalized to proper velocity
diff --git a/src/particles/diagnostics/DiagnosticOutput.cpp b/src/particles/diagnostics/DiagnosticOutput.cpp
index 0af5f3b7ac..538e1b4985 100644
--- a/src/particles/diagnostics/DiagnosticOutput.cpp
+++ b/src/particles/diagnostics/DiagnosticOutput.cpp
@@ -78,9 +78,9 @@ namespace impactx::diagnostics
 
                         // access AoS data such as positions and cpu/id
                         PType const &p = aos_ptr[i];
-                        amrex::ParticleReal const x = p.pos(0);
-                        amrex::ParticleReal const y = p.pos(1);
-                        amrex::ParticleReal const t = p.pos(2);
+                        amrex::ParticleReal const x = p.pos(RealAoS::x);
+                        amrex::ParticleReal const y = p.pos(RealAoS::y);
+                        amrex::ParticleReal const t = p.pos(RealAoS::t);
                         uint64_t const global_id = ablastr::particles::localIDtoGlobal(p.id(), p.cpu());
 
                         // access SoA Real data
@@ -122,8 +122,8 @@ namespace impactx::diagnostics
 
                         // access AoS data such as positions and cpu/id
                         PType const &p = aos_ptr[i];
-                        amrex::ParticleReal const x = p.pos(0);
-                        amrex::ParticleReal const y = p.pos(1);
+                        amrex::ParticleReal const x = p.pos(RealAoS::x);
+                        amrex::ParticleReal const y = p.pos(RealAoS::y);
                         uint64_t const global_id = ablastr::particles::localIDtoGlobal(p.id(), p.cpu());
 
                         // access SoA Real data
diff --git a/src/particles/elements/ConstF.H b/src/particles/elements/ConstF.H
index accbcc2acf..98ccdf6388 100644
--- a/src/particles/elements/ConstF.H
+++ b/src/particles/elements/ConstF.H
@@ -68,9 +68,9 @@ namespace impactx
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // access reference particle values to find beta*gamma^2
             amrex::ParticleReal const pt_ref = refpart.pt;
@@ -85,13 +85,13 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // advance position and momentum
-            p.pos(0) = cos(m_kx*slice_ds)*x + sin(m_kx*slice_ds)/m_kx*px;
+            p.pos(RealAoS::x) = cos(m_kx*slice_ds)*x + sin(m_kx*slice_ds)/m_kx*px;
             pxout = -m_kx*sin(m_kx*slice_ds)*x + cos(m_kx*slice_ds)*px;
 
-            p.pos(1) = cos(m_ky*slice_ds)*y + sin(m_ky*slice_ds)/m_ky*py;
+            p.pos(RealAoS::y) = cos(m_ky*slice_ds)*y + sin(m_ky*slice_ds)/m_ky*py;
             pyout = -m_ky*sin(m_ky*slice_ds)*y + cos(m_ky*slice_ds)*py;
 
-            p.pos(2) = cos(m_kt*slice_ds)*t + sin(m_kt*slice_ds)/(betgam2*m_kt)*pt;
+            p.pos(RealAoS::t) = cos(m_kt*slice_ds)*t + sin(m_kt*slice_ds)/(betgam2*m_kt)*pt;
             ptout = -(m_kt*betgam2)*sin(m_kt*slice_ds)*t + cos(m_kt*slice_ds)*pt;
 
             // assign updated momenta
diff --git a/src/particles/elements/DipEdge.H b/src/particles/elements/DipEdge.H
index 64e3025938..c86b2ee682 100644
--- a/src/particles/elements/DipEdge.H
+++ b/src/particles/elements/DipEdge.H
@@ -74,8 +74,8 @@ namespace impactx
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
 
             // access reference particle values if needed
 
diff --git a/src/particles/elements/Drift.H b/src/particles/elements/Drift.H
index 80decffe33..f14ba3fb9b 100644
--- a/src/particles/elements/Drift.H
+++ b/src/particles/elements/Drift.H
@@ -63,9 +63,9 @@ namespace impactx
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // initialize output values of momenta
             amrex::ParticleReal pxout = px;
@@ -80,11 +80,11 @@ namespace impactx
             amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt;
 
             // advance position and momentum (drift)
-            p.pos(0) = x + slice_ds * px;
+            p.pos(RealAoS::x) = x + slice_ds * px;
             // pxout = px;
-            p.pos(1) = y + slice_ds * py;
+            p.pos(RealAoS::y) = y + slice_ds * py;
             // pyout = py;
-            p.pos(2) = t + (slice_ds/betgam2) * pt;
+            p.pos(RealAoS::t) = t + (slice_ds/betgam2) * pt;
             // ptout = pt;
 
             // assign updated momenta
diff --git a/src/particles/elements/Multipole.H b/src/particles/elements/Multipole.H
index 80694953df..ef487a8c1d 100644
--- a/src/particles/elements/Multipole.H
+++ b/src/particles/elements/Multipole.H
@@ -76,9 +76,9 @@ namespace impactx
             using Complex = amrex::GpuComplex<amrex::ParticleReal>;
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // access reference particle values to find (beta*gamma)^2
             //amrex::ParticleReal const pt_ref = refpart.pt;
@@ -101,13 +101,13 @@ namespace impactx
             amrex::ParticleReal const dpy = kick.m_imag/m_mfactorial;
 
             // advance position and momentum
-            p.pos(0) = x;
+            p.pos(RealAoS::x) = x;
             pxout = px + dpx;
 
-            p.pos(1) = y;
+            p.pos(RealAoS::y) = y;
             pyout = py + dpy;
 
-            p.pos(2) = t;
+            p.pos(RealAoS::t) = t;
             ptout = pt;
 
             // assign updated momenta
diff --git a/src/particles/elements/NonlinearLens.H b/src/particles/elements/NonlinearLens.H
index 018e43a18e..42f6431ff0 100644
--- a/src/particles/elements/NonlinearLens.H
+++ b/src/particles/elements/NonlinearLens.H
@@ -73,9 +73,9 @@ namespace impactx
             using Complex = amrex::GpuComplex<amrex::ParticleReal>;
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // access reference particle values to find (beta*gamma)^2
             //amrex::ParticleReal const pt_ref = refpart.pt;
@@ -110,13 +110,13 @@ namespace impactx
             amrex::ParticleReal dpy = -kick*dF.m_imag;
 
             // advance position and momentum
-            p.pos(0) = x;
+            p.pos(RealAoS::x) = x;
             pxout = px + dpx;
 
-            p.pos(1) = y;
+            p.pos(RealAoS::y) = y;
             pyout = py + dpy;
 
-            p.pos(2) = t;
+            p.pos(RealAoS::t) = t;
             ptout = pt;
 
             // assign updated momenta
diff --git a/src/particles/elements/PRot.H b/src/particles/elements/PRot.H
index eab45e8949..87a712e7e9 100644
--- a/src/particles/elements/PRot.H
+++ b/src/particles/elements/PRot.H
@@ -71,9 +71,9 @@ namespace impactx
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // access reference particle values to find beta:
             amrex::ParticleReal const beta = refpart.beta();
@@ -91,13 +91,13 @@ namespace impactx
                  sin(m_phi_in))*sin(theta);
 
             // advance position and momentum
-            p.pos(0) = x*pz/pzf;
+            p.pos(RealAoS::x) = x*pz/pzf;
             pxout = px*cos(theta) + (pz - cos(m_phi_in))*sin(theta);
 
-            p.pos(1) = y + py*x*sin(theta)/pzf;
+            p.pos(RealAoS::y) = y + py*x*sin(theta)/pzf;
             pyout = py;
 
-            p.pos(2) = t - (pt - 1.0_prt/beta)*x*sin(theta)/pzf;
+            p.pos(RealAoS::t) = t - (pt - 1.0_prt/beta)*x*sin(theta)/pzf;
             ptout = pt;
 
             // assign updated momenta
diff --git a/src/particles/elements/Quad.H b/src/particles/elements/Quad.H
index 8d90c58b46..30fd2057b3 100644
--- a/src/particles/elements/Quad.H
+++ b/src/particles/elements/Quad.H
@@ -68,9 +68,9 @@ namespace impactx
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -89,23 +89,23 @@ namespace impactx
 
             if(m_k > 0.0) {
                // advance position and momentum (focusing quad)
-               p.pos(0) = cos(omega*slice_ds)*x + sin(omega*slice_ds)/omega*px;
+               p.pos(RealAoS::x) = cos(omega*slice_ds)*x + sin(omega*slice_ds)/omega*px;
                pxout = -omega*sin(omega*slice_ds)*x + cos(omega*slice_ds)*px;
 
-               p.pos(1) = cosh(omega*slice_ds)*y + sinh(omega*slice_ds)/omega*py;
+               p.pos(RealAoS::y) = cosh(omega*slice_ds)*y + sinh(omega*slice_ds)/omega*py;
                pyout = omega*sinh(omega*slice_ds)*y + cosh(omega*slice_ds)*py;
 
-               p.pos(2) = t + (slice_ds/betgam2)*pt;
+               p.pos(RealAoS::t) = t + (slice_ds/betgam2)*pt;
                // ptout = pt;
             } else {
                // advance position and momentum (defocusing quad)
-               p.pos(0) = cosh(omega*slice_ds)*x + sinh(omega*slice_ds)/omega*px;
+               p.pos(RealAoS::x) = cosh(omega*slice_ds)*x + sinh(omega*slice_ds)/omega*px;
                pxout = omega*sinh(omega*slice_ds)*x + cosh(omega*slice_ds)*px;
 
-               p.pos(1) = cos(omega*slice_ds)*y + sin(omega*slice_ds)/omega*py;
+               p.pos(RealAoS::y) = cos(omega*slice_ds)*y + sin(omega*slice_ds)/omega*py;
                pyout = -omega*sin(omega*slice_ds)*y + cos(omega*slice_ds)*py;
 
-               p.pos(2) = t + (slice_ds/betgam2)*pt;
+               p.pos(RealAoS::t) = t + (slice_ds/betgam2)*pt;
                // ptout = pt;
             }
 
diff --git a/src/particles/elements/RFCavity.H b/src/particles/elements/RFCavity.H
index 0825910a0a..ce6b7d8780 100644
--- a/src/particles/elements/RFCavity.H
+++ b/src/particles/elements/RFCavity.H
@@ -186,9 +186,9 @@ namespace RFCavityData
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // initialize output values of momenta
             amrex::ParticleReal pxout = px;
@@ -205,15 +205,15 @@ namespace RFCavityData
             // so that, e.g., R(3,4) = dyf/dpyi.
 
             // push particles using the linear map
-            p.pos(0) = R(1,1)*x + R(1,2)*px + R(1,3)*y
+            p.pos(RealAoS::x) = R(1,1)*x + R(1,2)*px + R(1,3)*y
                      + R(1,4)*py + R(1,5)*t + R(1,6)*pt;
             pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y
                   + R(2,4)*py + R(2,5)*t + R(2,6)*pt;
-            p.pos(1) = R(3,1)*x + R(3,2)*px + R(3,3)*y
+            p.pos(RealAoS::y) = R(3,1)*x + R(3,2)*px + R(3,3)*y
                      + R(3,4)*py + R(3,5)*t + R(3,6)*pt;
             pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y
                   + R(4,4)*py + R(4,5)*t + R(4,6)*pt;
-            p.pos(2) = R(5,1)*x + R(5,2)*px + R(5,3)*y
+            p.pos(RealAoS::t) = R(5,1)*x + R(5,2)*px + R(5,3)*y
                      + R(5,4)*py + R(5,5)*t + R(5,6)*pt;
             ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y
                   + R(6,4)*py + R(6,5)*t + R(6,6)*pt;
diff --git a/src/particles/elements/Sbend.H b/src/particles/elements/Sbend.H
index d0e79dcb69..ac7d49e3ad 100644
--- a/src/particles/elements/Sbend.H
+++ b/src/particles/elements/Sbend.H
@@ -65,9 +65,9 @@ namespace impactx
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // initialize output values of momenta
             amrex::ParticleReal pxout = px;
@@ -89,16 +89,16 @@ namespace impactx
             amrex::ParticleReal const cos_theta = cos(theta);
 
             // advance position and momentum (sector bend)
-            p.pos(0) = cos_theta*x + m_rc*sin_theta*px
+            p.pos(RealAoS::x) = cos_theta*x + m_rc*sin_theta*px
                        - (m_rc/bet)*(1.0_prt - cos_theta)*pt;
 
             pxout = -sin_theta/m_rc*x + cos_theta*px - sin_theta/bet*pt;
 
-            p.pos(1) = y + m_rc*theta*py;
+            p.pos(RealAoS::y) = y + m_rc*theta*py;
 
             // pyout = py;
 
-            p.pos(2) = sin_theta/bet*x + m_rc/bet*(1.0_prt - cos_theta)*px + t
+            p.pos(RealAoS::t) = sin_theta/bet*x + m_rc/bet*(1.0_prt - cos_theta)*px + t
                        + m_rc*(-theta+sin_theta/(bet*bet))*pt;
 
             // ptout = pt;
diff --git a/src/particles/elements/ShortRF.H b/src/particles/elements/ShortRF.H
index aa984674d1..a3ff4645ed 100644
--- a/src/particles/elements/ShortRF.H
+++ b/src/particles/elements/ShortRF.H
@@ -64,9 +64,9 @@ namespace impactx
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // access reference particle values to find (beta*gamma)^2
             amrex::ParticleReal const pt_ref = refpart.pt;
@@ -78,13 +78,13 @@ namespace impactx
             amrex::ParticleReal ptout = pt;
 
             // advance position and momentum
-            p.pos(0) = x;
+            p.pos(RealAoS::x) = x;
             pxout = px + m_k*m_V/(2.0_prt*betgam2)*x;
 
-            p.pos(1) = y;
+            p.pos(RealAoS::y) = y;
             pyout = py + m_k*m_V/(2.0_prt*betgam2)*y;
 
-            p.pos(2) = t;
+            p.pos(RealAoS::t) = t;
             ptout = pt - m_k*m_V*t;
 
             // assign updated momenta
diff --git a/src/particles/elements/SoftQuad.H b/src/particles/elements/SoftQuad.H
index 720ab222bd..885b27af01 100644
--- a/src/particles/elements/SoftQuad.H
+++ b/src/particles/elements/SoftQuad.H
@@ -191,9 +191,9 @@ namespace SoftQuadrupoleData
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // initialize output values of momenta
             amrex::ParticleReal pxout = px;
@@ -210,15 +210,15 @@ namespace SoftQuadrupoleData
             // so that, e.g., R(3,4) = dyf/dpyi.
 
             // push particles using the linear map
-            p.pos(0) = R(1,1)*x + R(1,2)*px + R(1,3)*y
+            p.pos(RealAoS::x) = R(1,1)*x + R(1,2)*px + R(1,3)*y
                      + R(1,4)*py + R(1,5)*t + R(1,6)*pt;
             pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y
                   + R(2,4)*py + R(2,5)*t + R(2,6)*pt;
-            p.pos(1) = R(3,1)*x + R(3,2)*px + R(3,3)*y
+            p.pos(RealAoS::y) = R(3,1)*x + R(3,2)*px + R(3,3)*y
                      + R(3,4)*py + R(3,5)*t + R(3,6)*pt;
             pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y
                   + R(4,4)*py + R(4,5)*t + R(4,6)*pt;
-            p.pos(2) = R(5,1)*x + R(5,2)*px + R(5,3)*y
+            p.pos(RealAoS::t) = R(5,1)*x + R(5,2)*px + R(5,3)*y
                      + R(5,4)*py + R(5,5)*t + R(5,6)*pt;
             ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y
                   + R(6,4)*py + R(6,5)*t + R(6,6)*pt;
diff --git a/src/particles/elements/SoftSol.H b/src/particles/elements/SoftSol.H
index cc6b470c9d..bbd7f4d40b 100644
--- a/src/particles/elements/SoftSol.H
+++ b/src/particles/elements/SoftSol.H
@@ -196,9 +196,9 @@ namespace SoftSolenoidData
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // initialize output values of momenta
             amrex::ParticleReal pxout = px;
@@ -215,15 +215,15 @@ namespace SoftSolenoidData
             // so that, e.g., R(3,4) = dyf/dpyi.
 
             // push particles using the linear map
-            p.pos(0) = R(1,1)*x + R(1,2)*px + R(1,3)*y
+            p.pos(RealAoS::x) = R(1,1)*x + R(1,2)*px + R(1,3)*y
                      + R(1,4)*py + R(1,5)*t + R(1,6)*pt;
             pxout = R(2,1)*x + R(2,2)*px + R(2,3)*y
                   + R(2,4)*py + R(2,5)*t + R(2,6)*pt;
-            p.pos(1) = R(3,1)*x + R(3,2)*px + R(3,3)*y
+            p.pos(RealAoS::y) = R(3,1)*x + R(3,2)*px + R(3,3)*y
                      + R(3,4)*py + R(3,5)*t + R(3,6)*pt;
             pyout = R(4,1)*x + R(4,2)*px + R(4,3)*y
                   + R(4,4)*py + R(4,5)*t + R(4,6)*pt;
-            p.pos(2) = R(5,1)*x + R(5,2)*px + R(5,3)*y
+            p.pos(RealAoS::t) = R(5,1)*x + R(5,2)*px + R(5,3)*y
                      + R(5,4)*py + R(5,5)*t + R(5,6)*pt;
             ptout = R(6,1)*x + R(6,2)*px + R(6,3)*y
                   + R(6,4)*py + R(6,5)*t + R(6,6)*pt;
diff --git a/src/particles/elements/Sol.H b/src/particles/elements/Sol.H
index 84c5a36f9f..0fb93c6f01 100644
--- a/src/particles/elements/Sol.H
+++ b/src/particles/elements/Sol.H
@@ -66,9 +66,9 @@ namespace impactx
             using namespace amrex::literals; // for _rt and _prt
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -106,13 +106,13 @@ namespace impactx
             pt = ptout;
 
             // advance positions and momenta using map for rotation
-            p.pos(0) = cos(theta)*xout + sin(theta)*yout;
+            p.pos(RealAoS::x) = cos(theta)*xout + sin(theta)*yout;
             pxout = cos(theta)*px + sin(theta)*py;
 
-            p.pos(1) = -sin(theta)*xout + cos(theta)*yout;
+            p.pos(RealAoS::y) = -sin(theta)*xout + cos(theta)*yout;
             pyout = -sin(theta)*px + cos(theta)*py;
 
-            p.pos(2) = tout;
+            p.pos(RealAoS::t) = tout;
             ptout = pt;
 
             // assign updated momenta
diff --git a/src/particles/elements/diagnostics/openPMD.cpp b/src/particles/elements/diagnostics/openPMD.cpp
index 27efade837..d61695842a 100644
--- a/src/particles/elements/diagnostics/openPMD.cpp
+++ b/src/particles/elements/diagnostics/openPMD.cpp
@@ -261,8 +261,8 @@ namespace detail
             beam["positionOffset"]["x"].makeConstant(ref_part.x);
             beam["positionOffset"]["y"].resetDataset(d_fl);
             beam["positionOffset"]["y"].makeConstant(ref_part.y);
-            beam["positionOffset"]["ct"].resetDataset(d_fl);
-            beam["positionOffset"]["ct"].makeConstant(ref_part.t);
+            beam["positionOffset"]["t"].resetDataset(d_fl);
+            beam["positionOffset"]["t"].makeConstant(ref_part.t);
         }
 
         // AoS: Int
@@ -374,7 +374,7 @@ namespace detail
         // AoS: position and particle ID
         {
             using vs = std::vector<std::string>;
-            vs const positionComponents{"x", "y", "ct"}; // TODO: generalize
+            vs const positionComponents{"x", "y", "t"}; // TODO: generalize
             for (auto currDim = 0; currDim < AMREX_SPACEDIM; currDim++) {
                 std::shared_ptr<amrex::ParticleReal> curr(
                     new amrex::ParticleReal[numParticleOnTile],
diff --git a/src/particles/spacecharge/GatherAndPush.cpp b/src/particles/spacecharge/GatherAndPush.cpp
index 10a4ce646f..4a85abebfd 100644
--- a/src/particles/spacecharge/GatherAndPush.cpp
+++ b/src/particles/spacecharge/GatherAndPush.cpp
@@ -19,10 +19,10 @@
 namespace impactx::spacecharge
 {
     void GatherAndPush (
-            ImpactXParticleContainer & pc,
-            std::unordered_map<int, std::unordered_map<std::string, amrex::MultiFab> > const & space_charge_field,
-            const amrex::Vector<amrex::Geometry>& geom,
-            amrex::ParticleReal const slice_ds
+        ImpactXParticleContainer & pc,
+        std::unordered_map<int, std::unordered_map<std::string, amrex::MultiFab> > const & space_charge_field,
+        const amrex::Vector<amrex::Geometry>& geom,
+        amrex::ParticleReal const slice_ds
     )
     {
         BL_PROFILE("impactx::spacecharge::GatherAndPush");
@@ -72,7 +72,7 @@ namespace impactx::spacecharge
                 auto& soa_real = pti.GetStructOfArrays().GetRealData();
                 amrex::ParticleReal* const AMREX_RESTRICT part_px = soa_real[RealSoA::ux].dataPtr();
                 amrex::ParticleReal* const AMREX_RESTRICT part_py = soa_real[RealSoA::uy].dataPtr();
-                amrex::ParticleReal* const AMREX_RESTRICT part_pz = soa_real[RealSoA::pt].dataPtr(); // note: currently in z
+                amrex::ParticleReal* const AMREX_RESTRICT part_pz = soa_real[RealSoA::pz].dataPtr(); // note: currently for a fixed t
 
                 // group together constants for the momentum push
                 amrex::ParticleReal const push_consts = dt * charge * inv_gamma2 / pz_ref_SI;
@@ -90,7 +90,7 @@ namespace impactx::spacecharge
                     // force gather
                     amrex::GpuArray<amrex::Real, 3> const field_interp =
                         ablastr::particles::doGatherVectorFieldNodal (
-                            p.pos(0), p.pos(1), p.pos(2),
+                            p.pos(RealAoS::x), p.pos(RealAoS::y), p.pos(RealAoS::t),
                             scf_arr_x, scf_arr_y, scf_arr_z,
                             invdr,
                             prob_lo);
diff --git a/src/particles/transformation/CoordinateTransformation.cpp b/src/particles/transformation/CoordinateTransformation.cpp
index 35187be401..7615354ede 100644
--- a/src/particles/transformation/CoordinateTransformation.cpp
+++ b/src/particles/transformation/CoordinateTransformation.cpp
@@ -52,10 +52,12 @@ namespace transformation {
                 auto &soa_real = pti.GetStructOfArrays().GetRealData();
                 amrex::ParticleReal *const AMREX_RESTRICT part_px = soa_real[RealSoA::ux].dataPtr();
                 amrex::ParticleReal *const AMREX_RESTRICT part_py = soa_real[RealSoA::uy].dataPtr();
-                amrex::ParticleReal *const AMREX_RESTRICT part_pt = soa_real[RealSoA::pt].dataPtr();
 
                 if( direction == Direction::to_fixed_s) {
                     BL_PROFILE("impactx::transformation::CoordinateTransformation::to_fixed_s");
+
+                    amrex::ParticleReal *const AMREX_RESTRICT part_pz = soa_real[RealSoA::pz].dataPtr();
+
                     // Design value of pz/mc = beta*gamma
                     amrex::ParticleReal const pzd = sqrt(pow(pd, 2) - 1.0);
 
@@ -67,12 +69,15 @@ namespace transformation {
                         // access SoA Real data
                         amrex::ParticleReal &px = part_px[i];
                         amrex::ParticleReal &py = part_py[i];
-                        amrex::ParticleReal &pz = part_pt[i];
+                        amrex::ParticleReal &pz = part_pz[i];
 
                         to_s(p, px, py, pz);
                     });
                 } else {
                     BL_PROFILE("impactx::transformation::CoordinateTransformation::to_fixed_t");
+
+                    amrex::ParticleReal *const AMREX_RESTRICT part_pt = soa_real[RealSoA::pt].dataPtr();
+
                     amrex::ParticleReal const ptd = pd;  // Design value of pt/mc2 = -gamma.
                     ToFixedT const to_t(ptd);
                     amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(long i) {
diff --git a/src/particles/transformation/ToFixedS.H b/src/particles/transformation/ToFixedS.H
index 23887501fe..859f08807e 100644
--- a/src/particles/transformation/ToFixedS.H
+++ b/src/particles/transformation/ToFixedS.H
@@ -34,7 +34,7 @@ namespace transformation
          *
          * @param pzd Design value of pz/mc = beta*gamma.
          */
-        ToFixedS( amrex::ParticleReal const pzd )
+        ToFixedS (amrex::ParticleReal const pzd)
         : m_pzd(pzd)
         {
         }
@@ -42,24 +42,24 @@ namespace transformation
         /** This is a t-to-s map, so that a variable of this type can be used like a
          *  t-to-s function.
          *
-         * @param p Particle AoS data for positions and cpu/id
-         * @param px particle momentum in x
-         * @param py particle momentum in y
-         * @param pt particle momentum in t
+         * @param[inout] p Particle AoS data for positions and cpu/id
+         * @param[inout] px particle momentum in x
+         * @param[inout] py particle momentum in y
+         * @param[inout] pz particle momentum in z (in), in t (out)
          */
         AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
         void operator() (
-                PType& p,
-                amrex::ParticleReal & px,
-                amrex::ParticleReal & py,
-                amrex::ParticleReal & pt) const
+            PType& p,
+            amrex::ParticleReal & px,
+            amrex::ParticleReal & py,
+            amrex::ParticleReal & pz) const
         {
             using namespace amrex::literals;
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const z = p.pos(RealAoS::z);
 
             // compute value of reference ptd = -gamma
             amrex::ParticleReal const argd = 1.0_prt + pow(m_pzd, 2);
@@ -68,34 +68,33 @@ namespace transformation
 
             // transform momenta to dynamic units (e.g., so that momenta are
             // normalized by mc):
-            px = px*m_pzd;
-            py = py*m_pzd;
-            pt = pt*m_pzd;
+            px = px * m_pzd;
+            py = py * m_pzd;
+            pz = pz * m_pzd;
 
             // compute value of particle pt = -gamma
-            amrex::ParticleReal const arg = 1.0_prt + pow(m_pzd+pt, 2) + pow(px, 2) + pow(py, 2);
+            amrex::ParticleReal const arg = 1.0_prt + pow(m_pzd + pz, 2) + pow(px, 2) + pow(py, 2);
             AMREX_ASSERT_WITH_MESSAGE(arg > 0.0_prt, "invalid pt arg (<=0)");
             amrex::ParticleReal const ptf = arg > 0.0_prt ? -sqrt(arg) : -1.0_prt;
 
             // transform position and momentum (from fixed t to fixed s)
-
-            p.pos(0) = x - px*t/(m_pzd+pt);
+            p.pos(RealAoS::x) = x - px * z / (m_pzd + pz);
             // px = px;
-            p.pos(1) = y - py*t/(m_pzd+pt);
+            p.pos(RealAoS::y) = y - py * z / (m_pzd + pz);
             // py = py;
-            p.pos(2) = ptf*t/(m_pzd+pt);  // This now represents t.
-            pt = ptf - ptdf;                    // This now represents pt.
+            p.pos(RealAoS::t) = ptf * z / (m_pzd + pz);
+            auto & pt = pz;  // We store pt in the same memory slot as pz.
+            pt = ptf - ptdf;
 
             // transform momenta to static units (eg, so that momenta are
             // normalized by pzd):
-            px = px/m_pzd;
-            py = py/m_pzd;
-            pt = pt/m_pzd;
-
+            px = px / m_pzd;
+            py = py / m_pzd;
+            pt = pz / m_pzd;
         }
 
     private:
-        amrex::ParticleReal m_pzd;
+        amrex::ParticleReal m_pzd;  ///< Design value of pz/mc = beta*gamma.
     };
 
 } // namespace transformation
diff --git a/src/particles/transformation/ToFixedT.H b/src/particles/transformation/ToFixedT.H
index 0dbedf5df9..1515075b50 100644
--- a/src/particles/transformation/ToFixedT.H
+++ b/src/particles/transformation/ToFixedT.H
@@ -34,7 +34,7 @@ namespace transformation
          *
          * @param ptd Design value of pt/mc2 = -gamma.
          */
-        ToFixedT( amrex::ParticleReal const ptd )
+        ToFixedT (amrex::ParticleReal const ptd)
         : m_ptd(ptd)
         {
         }
@@ -42,60 +42,59 @@ namespace transformation
         /** This is a s-to-t map, so that a variable of this type can be used like a
          *  s-to-t function.
          *
-         * @param p Particle AoS data for positions and cpu/id
-         * @param px particle momentum in x
-         * @param py particle momentum in y
-         * @param pt particle momentum in t
+         * @param[inout] p Particle AoS data for positions and cpu/id
+         * @param[inout] px particle momentum in x
+         * @param[inout] py particle momentum in y
+         * @param[inout] pt particle momentum in t (in), in z (out)
          */
         AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
         void operator() (
-                PType& p,
-                amrex::ParticleReal & px,
-                amrex::ParticleReal & py,
-                amrex::ParticleReal & pt) const
+            PType& p,
+            amrex::ParticleReal & px,
+            amrex::ParticleReal & py,
+            amrex::ParticleReal & pt) const
         {
             using namespace amrex::literals;
 
             // access AoS data such as positions and cpu/id
-            amrex::ParticleReal const x = p.pos(0);
-            amrex::ParticleReal const y = p.pos(1);
-            amrex::ParticleReal const t = p.pos(2);
+            amrex::ParticleReal const x = p.pos(RealAoS::x);
+            amrex::ParticleReal const y = p.pos(RealAoS::y);
+            amrex::ParticleReal const t = p.pos(RealAoS::t);
 
             // compute value of reference pzd = beta*gamma
             amrex::ParticleReal const argd = -1.0_prt + pow(m_ptd, 2);
             AMREX_ASSERT_WITH_MESSAGE(argd > 0.0_prt, "invalid pzd arg (<=0)");
-            amrex::ParticleReal const pzd = argd > 0.0_prt ? sqrt(argd) : 0.0_prt;
+            amrex::ParticleReal const pzdf = argd > 0.0_prt ? sqrt(argd) : 0.0_prt;
 
             // transform momenta to dynamic units (eg, so that momenta are
             // normalized by mc):
-            px = px*pzd;
-            py = py*pzd;
-            pt = pt*pzd;
+            px = px * pzdf;
+            py = py * pzdf;
+            pt = pt * pzdf;
 
             // compute value of particle pz = beta*gamma
             amrex::ParticleReal const arg = -1.0_prt + pow(m_ptd+pt, 2) - pow(px, 2) - pow(py, 2);
             AMREX_ASSERT_WITH_MESSAGE(arg > 0.0_prt, "invalid pz arg (<=0)");
-            amrex::ParticleReal const pz = arg > 0.0_prt ? sqrt(arg) : 0.0_prt;
+            amrex::ParticleReal const pzf = arg > 0.0_prt ? sqrt(arg) : 0.0_prt;
 
             // transform position and momentum (from fixed s to fixed t)
-
-            p.pos(0) = x + px*t/(m_ptd+pt);
+            p.pos(RealAoS::x) = x + px*t/(m_ptd+pt);
             // px = px;
-            p.pos(1) = y + py*t/(m_ptd+pt);
+            p.pos(RealAoS::y) = y + py*t/(m_ptd+pt);
             // py = py;
-            p.pos(2) = pz*t/(m_ptd+pt);  // This now represents z.
-            pt = pz - pzd;                     // This now represents pz.
+            p.pos(RealAoS::z) = pzf * t / (m_ptd + pt);
+            auto & pz = pt;  // We store pz in the same memory slot as pt.
+            pz = pzf - pzdf;
 
             // transform momenta to static units (eg, so that momenta are
-            // normalized by pzd):
-            px = px/pzd;
-            py = py/pzd;
-            pt = pt/pzd;
-
+            // normalized by pzdf):
+            px = px / pzdf;
+            py = py / pzdf;
+            pz = pz / pzdf;
         }
 
     private:
-        amrex::ParticleReal m_ptd;
+        amrex::ParticleReal m_ptd;  ///< Design value of pt/mc2 = -gamma.
     };
 
 } // namespace transformation
diff --git a/src/python/ImpactXParticleContainer.cpp b/src/python/ImpactXParticleContainer.cpp
index 983d3d2066..dae8d2c71d 100644
--- a/src/python/ImpactXParticleContainer.cpp
+++ b/src/python/ImpactXParticleContainer.cpp
@@ -44,20 +44,20 @@ void init_impactxparticlecontainer(py::module& m)
         .def("add_n_particles",
              &ImpactXParticleContainer::AddNParticles,
              py::arg("lev"),
-             py::arg("x"), py::arg("y"), py::arg("z"),
-             py::arg("px"), py::arg("py"), py::arg("pz"),
+             py::arg("x"), py::arg("y"), py::arg("t"),
+             py::arg("px"), py::arg("py"), py::arg("pt"),
              py::arg("qm"), py::arg("bchchg"),
-             "Add new particles to the container\n\n"
+             "Add new particles to the container for fixed s.\n\n"
              "Note: This can only be used *after* the initialization (grids) have\n"
              "      been created, meaning after the call to ImpactX.init_grids\n"
              "      has been made in the ImpactX class.\n\n"
              ":param lev: mesh-refinement level\n"
              ":param x: positions in x\n"
              ":param y: positions in y\n"
-             ":param z: positions in z\n"
+             ":param t: positions as time-of-flight in c*t\n"
              ":param px: momentum in x\n"
              ":param py: momentum in y\n"
-             ":param pz: momentum in z\n"
+             ":param pt: momentum in t\n"
              ":param qm: charge over mass in 1/eV\n"
              ":param bchchg: total charge within a bunch in C"
         )