diff --git a/docs/source/dataanalysis/dataanalysis.rst b/docs/source/dataanalysis/dataanalysis.rst
index cfae75ce2..67d1e2d8f 100644
--- a/docs/source/dataanalysis/dataanalysis.rst
+++ b/docs/source/dataanalysis/dataanalysis.rst
@@ -69,6 +69,8 @@ The code writes out the values in an ASCII file prefixed ``reduced_beam_characte
 
 * ``step``
     Iteration within the simulation
+* ``s``
+    Reference particle coordinate ``s`` (unit: meter)
 * ``x_mean/min/max``, ``y_mean/min/max``, ``t_mean/min/max``
     Average / minimum / maximum beam particle position in the dimensions of ``x``, ``y`` (transverse coordinates, unit: meter), and ``t`` (normalized time difference :math:`ct`, unit: meter)
 * ``sig_x``, ``sig_y``, ``sig_t``
diff --git a/src/particles/diagnostics/DiagnosticOutput.cpp b/src/particles/diagnostics/DiagnosticOutput.cpp
index 9ea12082d..1b0c6190b 100644
--- a/src/particles/diagnostics/DiagnosticOutput.cpp
+++ b/src/particles/diagnostics/DiagnosticOutput.cpp
@@ -43,7 +43,7 @@ namespace impactx::diagnostics
             if (otype == OutputType::PrintRefParticle) {
                 file_handler << "step s beta gamma beta_gamma x y z t px py pz pt\n";
             } else if (otype == OutputType::PrintReducedBeamCharacteristics) {
-                file_handler << "step" << " "
+                file_handler << "step" << " " << "s" << " "
                              << "x_mean" << " " << "x_min" << " " << "x_max" << " "
                              << "y_mean" << " " << "y_min" << " " << "y_max" << " "
                              << "t_mean" << " " << "t_min" << " " << "t_max" << " "
@@ -88,7 +88,9 @@ namespace impactx::diagnostics
             std::unordered_map<std::string, amrex::ParticleReal> const rbc =
                 diagnostics::reduced_beam_characteristics(pc);
 
-            file_handler << step << " "
+            amrex::ParticleReal const s = pc.GetRefParticle().s;
+
+            file_handler << step << " " << s << " "
                          << rbc.at("x_mean") << " " << rbc.at("x_min") << " " << rbc.at("x_max") << " "
                          << rbc.at("y_mean") << " " << rbc.at("y_min") << " " << rbc.at("y_max") << " "
                          << rbc.at("t_mean") << " " << rbc.at("t_min") << " " << rbc.at("t_max") << " "