diff --git a/pyfstat/core.py b/pyfstat/core.py
index 0fd551788f93140cd11852f7f5c4f4a6f4975578..130039dca0b9aa92996c3d0ef2782637df9152a1 100755
--- a/pyfstat/core.py
+++ b/pyfstat/core.py
@@ -116,6 +116,7 @@ class BaseSearchClass(object):
                 post_theta_at_ith_glitch = pre_theta_at_ith_glitch + dt
                 thetas.append(self.shift_coefficients(
                     post_theta_at_ith_glitch, self.tref - tbounds[i+1]))
+        self.thetas_at_tref = thetas
         return thetas
 
     def generate_loudest(self):
diff --git a/pyfstat/mcmc_based_searches.py b/pyfstat/mcmc_based_searches.py
index 8b3ddab4d6115a8fb2b13632c12af15c6da18c76..a32046a54317c5d092f35a28be91582e9a39a65e 100644
--- a/pyfstat/mcmc_based_searches.py
+++ b/pyfstat/mcmc_based_searches.py
@@ -851,6 +851,26 @@ class MCMCSearch(BaseSearchClass):
             d[k+'_std'] = np.std(s)
         return d
 
+    def check_if_samples_are_railing(self, threshold=0.01):
+        return_flag = False
+        for s, k in zip(self.samples.T, self.theta_keys):
+            prior = self.theta_prior[k]
+            if prior['type'] == 'unif':
+                prior_range = prior['upper'] - prior['lower']
+                edges = []
+                fracs = []
+                for l in ['lower', 'upper']:
+                    bools = np.abs(s - prior[l])/prior_range < threshold
+                    if np.any(bools):
+                        edges.append(l)
+                        fracs.append(str(100*float(np.sum(bools))/len(bools)))
+                if len(edges) > 0:
+                    logging.warning(
+                        '{}% of the {} posterior is railing on the {} edges'
+                        .format('% & '.join(fracs), k, ' & '.join(edges)))
+                    return_flag = True
+        return return_flag
+
     def write_par(self, method='med'):
         """ Writes a .par of the best-fit params with an estimated std """
         logging.info('Writing {}/{}.par using the {} method'.format(