diff --git a/studyprocessor.py b/studyprocessor.py
index c86de6c..e0ce27c 100644
--- a/studyprocessor.py
+++ b/studyprocessor.py
@@ -438,6 +438,34 @@ class StudyProcessor(Aggregator, HDF5Container):
                 series[key][run] = series[key][run][ind_sort]
         return series
 
+    def get_coldens0(self, run):
+        mp =  1.4 * 1.66 * 10**(-24) * U.g
+        z0 = self.get_nml("galbox_params/height0", run) * U.pc
+        n0 = self.get_nml("galbox_params/dens0", run) * U.cm**(-3)
+        return (np.sqrt(2 * np.pi) * mp * z0 * n0).express(U.coldens)
+
+    def total_mass(self):
+        """
+        Sum of gas plus sink mass
+        """
+        
+        time_gas = self.get_value("/series/coarse_step_from_log/time")        
+        mass_gas = self.get_value("/series/coarse_step_from_log/mcons")
+        mass_sink = self.get_value("/series/sinks_from_log/mass_sink")
+
+        total_mass = dict()
+        for run in self.runs:
+
+            # A bit specific ... needs to be generalized (TODO)
+            info = self.snaps[run][self.nums[run][0]].info
+            surface = (info["unit_length"].express(U.pc)) ** 2
+            m0 = self.get_coldens0(run) * surface    # Initial mass in  Msun
+            offset = time_gas[run].size - mass_sink[run].size
+            mass_gas[run] = m0 + m0*mass_gas[run] # convert in  Msun
+            total_mass[run] = mass_gas[run].copy()
+            total_mass[run][offset:] = mass_gas[run][offset:] + mass_sink[run] # re add sink_mass            
+        return time_gas, total_mass, mass_gas
+
     def _ssfr_from_mass_sink(self, avg_window=None):
         """
         avg_window in year