diff --git a/mbtrack2/utilities/beamloading.py b/mbtrack2/utilities/beamloading.py
index cfd3548defab3c657001f484d15dbf4740dc7997..d8d64b1e842ddf6c5d7cc26016f6aaf17224aedf 100644
--- a/mbtrack2/utilities/beamloading.py
+++ b/mbtrack2/utilities/beamloading.py
@@ -69,7 +69,7 @@ class BeamLoadingEquilibrium():
         self.N = N
         self.z0 = np.linspace(self.B1, self.B2, self.N)
         self.tau0 = self.z0/c
-        self.rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0)
+        self.rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0)/c
 
         # Define constants for scaled potential u(z)
         self.u0 = self.ring.U0 / (self.ring.ac * self.ring.sigma_delta**2 *
@@ -394,3 +394,9 @@ class BeamLoadingEquilibrium():
                 (2 * np.pi * HC.m**2 * self.F[1] * self.ring.h * I0 * f))
 
         return (eta, RQth, f)
+    
+    @property
+    def R_factor(self):
+        rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0)/c
+        R = trapz(rho0**2, self.tau0)/trapz(self.rho0**2, self.tau0)
+        return R