diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py
index 58d7170..44893f1 100755
--- a/bin/rabbit_fit.py
+++ b/bin/rabbit_fit.py
@@ -149,6 +149,12 @@ def make_parser():
         type=str,
         help="Specify result from external postfit file",
     )
+    parser.add_argument(
+        "--noPostfitProfileBB",
+        default=False,
+        action="store_true",
+        help="Do not profile the bin-by-bin parameters in the postfit (e.g. if parameters are loaded from another fit using --externalPostfit and/or no data is available to be profiled).",
+    )
     parser.add_argument(
         "--doImpacts",
         default=False,
@@ -276,41 +282,46 @@ def fit(args, fitter, ws, dofit=True):
 
     if args.externalPostfit is not None:
         fitter.load_fitresult(args.externalPostfit, args.externalPostfitResult)
-    else:
-        if dofit:
-            cb = fitter.minimize()
 
-            if cb is not None:
-                ws.add_loss_time_hist(cb.loss_history, cb.time_history)
+    if dofit:
+        cb = fitter.minimize()
+
+        # force profiling of beta with final parameter values
+        # TODO avoid the extra calculation and jitting if possible since the relevant calculation
+        # usually would have been done during the minimization
+        if fitter.binByBinStat and not args.noPostfitProfileBB:
+            fitter._profile_beta()
 
-        if not args.noHessian:
-            # compute the covariance matrix and estimated distance to minimum
+        if cb is not None:
+            ws.add_loss_time_hist(cb.loss_history, cb.time_history)
 
-            val, grad, hess = fitter.loss_val_grad_hess()
-            edmval, cov = edmval_cov(grad, hess)
+    if not args.noHessian:
+        # compute the covariance matrix and estimated distance to minimum
 
-            logger.info(f"edmval: {edmval}")
+        val, grad, hess = fitter.loss_val_grad_hess()
+        edmval, cov = edmval_cov(grad, hess)
+        logger.info(f"edmval: {edmval}")
 
-            fitter.cov.assign(cov)
+        fitter.cov.assign(cov)
 
-            del cov
+        del cov
 
-            if fitter.binByBinStat and fitter.diagnostics:
-                # This is the estimated distance to minimum with respect to variations of
-                # the implicit binByBinStat nuisances beta at fixed parameter values.
-                # It should be near-zero by construction as long as the analytic profiling is
-                # correct
-                _, gradbeta, hessbeta = fitter.loss_val_grad_hess_beta()
-                edmvalbeta, covbeta = edmval_cov(gradbeta, hessbeta)
-                logger.info(f"edmvalbeta: {edmvalbeta}")
+        if fitter.binByBinStat and fitter.diagnostics:
+            # This is the estimated distance to minimum with respect to variations of
+            # the implicit binByBinStat nuisances beta at fixed parameter values.
+            # It should be near-zero by construction as long as the analytic profiling is
+            # correct
+            _, gradbeta, hessbeta = fitter.loss_val_grad_hess_beta()
+            edmvalbeta, covbeta = edmval_cov(gradbeta, hessbeta)
+            logger.info(f"edmvalbeta: {edmvalbeta}")
 
-            if args.doImpacts:
-                ws.add_impacts_hists(*fitter.impacts_parms(hess))
+        if args.doImpacts:
+            ws.add_impacts_hists(*fitter.impacts_parms(hess))
 
-            del hess
+        del hess
 
-            if args.globalImpacts:
-                ws.add_global_impacts_hists(*fitter.global_impacts_parms())
+        if args.globalImpacts:
+            ws.add_global_impacts_hists(*fitter.global_impacts_parms())
 
     nllvalreduced = fitter.reduced_nll().numpy()
 
@@ -331,7 +342,7 @@ def fit(args, fitter, ws, dofit=True):
             "nllvalreduced": nllvalreduced,
             "ndfsat": ndfsat,
             "edmval": edmval,
-            "postfit_profile": args.externalPostfit is None,
+            "postfit_profile": not args.noPostfitProfileBB,
         }
     )
 
@@ -559,7 +570,7 @@ def main():
                             ifitter,
                             ws,
                             prefit=False,
-                            profile=args.externalPostfit is None,
+                            profile=not args.noPostfitProfileBB,
                         )
                 else:
                     fit_time.append(time.time())
diff --git a/bin/rabbit_plot_hists.py b/bin/rabbit_plot_hists.py
index bb66e0a..9d92e8d 100755
--- a/bin/rabbit_plot_hists.py
+++ b/bin/rabbit_plot_hists.py
@@ -817,16 +817,16 @@ def make_plot(
             h_den = h_inclusive
 
         if diff:
-            h1 = hh.addHists(h_inclusive, h_den, scale2=-1)
+            h0 = hh.addHists(h_num, h_num, scale2=-1)
             h2 = hh.addHists(h_data, h_den, scale2=-1)
             if h_data_stat is not None:
                 h2_stat = hh.divideHists(
                     h_data_stat, h_den, cutoff=cutoff, rel_unc=True
                 )
         else:
-            h1 = hh.divideHists(
-                h_inclusive,
-                h_den,
+            h0 = hh.divideHists(
+                h_num,
+                h_num,
                 cutoff=1e-8,
                 rel_unc=True,
                 flow=False,
@@ -839,7 +839,7 @@ def make_plot(
                 )
 
         hep.histplot(
-            h1,
+            h0,
             histtype="step",
             color="grey",
             alpha=0.5,
@@ -1311,14 +1311,7 @@ def make_plots(
                 h_data_stat = None
 
             if hists_up is not None:
-                hup = [
-                    (
-                        h[{k.replace("Sig", ""): v for k, v in idxs.items()}]
-                        if h is not None
-                        else None
-                    )
-                    for h in hists_up
-                ]
+                hup = [h[idxs] if h is not None else None for h in hists_up]
             else:
                 hup = None
 
diff --git a/bin/rabbit_plot_pulls_and_impacts.py b/bin/rabbit_plot_pulls_and_impacts.py
index d70b50b..d109aba 100755
--- a/bin/rabbit_plot_pulls_and_impacts.py
+++ b/bin/rabbit_plot_pulls_and_impacts.py
@@ -206,6 +206,9 @@ def plotImpacts(
         else:
             text_on_bars = True
 
+    name_suffix = f" ({name})" if name else ""
+    ref_name_suffix = f" ({ref_name})" if ref_name else ""
+
     if impacts:
 
         def make_bar(
@@ -234,7 +237,7 @@ def make_bar(
             )
 
         sign = "+/-" if (oneSidedImpacts and not any(df["impact_down"] > 0)) else "+"
-        label = f"{sign}1σ impact ({name})" if name else f"{sign}1σ impact"
+        label = f"{sign}1σ impact{name_suffix}"
         fig.add_trace(
             make_bar(
                 key="impact_up",
@@ -249,7 +252,7 @@ def make_bar(
             fig.add_trace(
                 make_bar(
                     key="impact_up_ref",
-                    name=f"{sign}1σ impact ({ref_name})",
+                    name=f"{sign}1σ impact{ref_name_suffix}",
                     filled=False,
                 ),
                 row=1,
@@ -260,7 +263,7 @@ def make_bar(
             fig.add_trace(
                 make_bar(
                     key="impact_down",
-                    name=f"-1σ impact ({name})" if name else "-1σ impact",
+                    name=f"-1σ impact{name_suffix}",
                     color="#e41a1c",
                     text_on_bars=text_on_bars,
                     opacity=0.5 if include_ref else 1,
@@ -272,7 +275,7 @@ def make_bar(
                 fig.add_trace(
                     make_bar(
                         key="impact_down_ref",
-                        name=f"-1σ impact ({ref_name})",
+                        name=f"-1σ impact{ref_name_suffix}",
                         color="#e41a1c",
                         filled=False,
                     ),
@@ -325,7 +328,7 @@ def make_bar(
                     size=8,
                 ),
                 error_x=error_x,
-                name="Pulls ± Constraints",
+                name=f"Pulls ± Constraints{name_suffix}",
                 showlegend=include_ref,
             ),
             row=1,
@@ -346,7 +349,7 @@ def make_bar(
                     y=labels,
                     orientation="h",
                     **get_marker(filled=False, color="black"),
-                    name=f"Pulls ± Constraints ({ref_name})",
+                    name=f"Pulls ± Constraints{ref_name_suffix}",
                     showlegend=True,
                 ),
                 row=1,
@@ -386,7 +389,7 @@ def make_bar(
                                 width=1
                             ),  # Adjust the thickness of the marker lines
                         ),
-                        name=f"Asym. pulls ({ref_name})",
+                        name=f"Asym. pulls{ref_name_suffix}",
                         showlegend=include_ref,
                     ),
                     row=1,
@@ -760,11 +763,13 @@ def parseArgs():
     parser.add_argument(
         "--refName",
         type=str,
+        default="ref.",
         help="Name of reference input for legend",
     )
     parser.add_argument(
         "--name",
         type=str,
+        default=None,
         help="Name of input for legend",
     )
     parser.add_argument(
diff --git a/rabbit/fitter.py b/rabbit/fitter.py
index faba62a..bb5e46a 100644
--- a/rabbit/fitter.py
+++ b/rabbit/fitter.py
@@ -288,12 +288,14 @@ def __init__(
 
     def load_fitresult(self, fitresult_file, fitresult_key):
         # load results from external fit and set postfit value and covariance elements for common parameters
+        cov_ext = None
         with h5py.File(fitresult_file, "r") as fext:
             if "x" in fext.keys():
                 # fitresult from combinetf
                 x_ext = fext["x"][...]
                 parms_ext = fext["parms"][...].astype(str)
-                cov_ext = fext["cov"][...]
+                if "cov" in fext.keys():
+                    cov_ext = fext["cov"][...]
             else:
                 # fitresult from rabbit
                 h5results_ext = io_tools.get_fitresult(fext, fitresult_key)
@@ -301,10 +303,10 @@ def load_fitresult(self, fitresult_file, fitresult_key):
 
                 x_ext = h_parms_ext.values()
                 parms_ext = np.array(h_parms_ext.axes["parms"])
-                cov_ext = h5results_ext["cov"].get().values()
+                if "cov" in h5results_ext.keys():
+                    cov_ext = h5results_ext["cov"].get().values()
 
         xvals = self.x.numpy()
-        covval = self.cov.numpy()
         parms = self.parms.astype(str)
 
         # Find common elements with their matching indices
@@ -312,10 +314,13 @@ def load_fitresult(self, fitresult_file, fitresult_key):
             parms, parms_ext, assume_unique=True, return_indices=True
         )
         xvals[idxs] = x_ext[idxs_ext]
-        covval[np.ix_(idxs, idxs)] = cov_ext[np.ix_(idxs_ext, idxs_ext)]
 
         self.x.assign(xvals)
-        self.cov.assign(tf.constant(covval))
+
+        if cov_ext is not None:
+            covval = self.cov.numpy()
+            covval[np.ix_(idxs, idxs)] = cov_ext[np.ix_(idxs_ext, idxs_ext)]
+            self.cov.assign(tf.constant(covval))
 
     def update_frozen_params(self):
         new_mask_np = np.isin(self.parms, self.frozen_params)
@@ -2152,12 +2157,6 @@ def scipy_hess(xval):
 
             self.x.assign(xval)
 
-        # force profiling of beta with final parameter values
-        # TODO avoid the extra calculation and jitting if possible since the relevant calculation
-        # usually would have been done during the minimization
-        if self.binByBinStat:
-            self._profile_beta()
-
         return callback
 
     def nll_scan(self, param, scan_range, scan_points, use_prefit=False):