diff --git a/PeakOTron.py b/PeakOTron.py
index 445084ae53a8e20fd962e33d8b896dc73beb747b..0e3d8576b9a662b71ab94a7619b393d5315dd2e8 100644
--- a/PeakOTron.py
+++ b/PeakOTron.py
@@ -1139,8 +1139,100 @@ class PeakOTron:
         else:
             plt.close(fig)
         
-        
-        
+    def PlotFitSumlab(self,figsize=(10.0,10.0),labelsize=None,ticksize=None,titlesize=None,
+                legendsize=None,title=None,xlabel = None,scaled=False,save_directory=None,
+                label=None,y_limits = [1, None],residual_limits = [-5, 5],residualticksize=None,
+                linewidth_main = 5,  x_limits = [None, None],display=True,prefit=False,plot_in_bins=False):
+        large = 25
+        med = 20
+        small = 12
+        if labelsize is None: labelsize = large
+        if ticksize is None: ticksize = large
+        if titlesize is None: titlesize = large
+        if legendsize is None: legendsize = med
+        if residualticksize is None: residualticksize = med
+        if(plot_in_bins):
+            x_all = self._hist_data["bin_numbers"]
+            bw = 1
+        else:
+            x_all = self._hist_data["bin_centres"]
+            bw = self._hist_data["bw"] 
+        y_all = self._hist_data["counts"]
+        y_err_all = self._hist_data["counts_error"]
+        Q_0_prefit = self._prefit_values_bin["Q_0"]
+        sigma_0_prefit = self._prefit_values_bin["sigma_0"]
+        if(self._hist_data["truncate_nsigma0_do_min"] is not None):
+            lim_nsigma0 =  Q_0_prefit - self._hist_data["truncate_nsigma0_do_min"]*sigma_0_prefit
+            cond_sel_nsigma0 = (self._hist_data["bin_numbers"]>lim_nsigma0)
+        else:
+            cond_sel_nsigma0 = np.full(len(x_all), True).astype(bool)
+        x = x_all[cond_sel_nsigma0]
+        y = y_all[cond_sel_nsigma0]
+        y_err = y_err_all[cond_sel_nsigma0]        
+        y_hat = self.GetModel(x, prefit=prefit, bin_units=plot_in_bins)*bw
+        y_hat_err = np.sqrt(y_hat)
+        chi2, ndf = self.GetChi2(prefit=prefit)
+
+        if(xlabel is None):
+            if(plot_in_bins): xlabel="Bin"
+            else: xlabel = "Q"
+        fig = plt.figure(figsize=figsize)
+        gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
+        ax0 = plt.subplot(gs[0])
+        ax0.tick_params(axis="x", direction='in', top=True, labeltop=False, length=small, width=1)
+        ax0.tick_params(axis="y", direction='in', right=True, labelright=False, length=small, width=1)
+        if(label is None): label = "Histogram" 
+        label = "{:s}".format(label)
+        dusty = (0.3, 0.6, 1.0)
+        brown = (0.35, 0, 0)
+        ax0.step(x_all, y_all, label=label, lw=0.8, color="grey", where="mid")
+        ax0.plot(x, y_hat, linestyle="--", label="Fit", lw=2, color='red')
+        ax0.plot([], [], ' ', label="$\\chi^{{2}}$/NDF = {:s}".format(LatexFormat(chi2/ndf)))
+        ax0.set_ylabel("Entries".format(xlabel, xlabel), fontsize=labelsize, loc='top')
+        ax0.tick_params(axis="y", labelsize=ticksize)
+        if(y_limits[0] is None): y_limits[0] = np.min(y[y>0])
+        ax0.set_ylim(y_limits)
+        ax0.set_xlim(x_limits)
+        if(title is not None): ax0.set_title(title, fontsize=titlesize)
+        ax0.grid(which="both", axis="both", linewidth=0.2, color="grey", linestyle=":")
+        ax0.locator_params(axis='y', nbins=8)
+        #ax0.set_yscale("log")
+        ax0.legend(fontsize=legendsize, loc='upper right', frameon=True, edgecolor='None')
+
+        ax1 = plt.subplot(gs[1], sharex = ax0)
+        ax1.tick_params(axis="x", direction='in', top=True, labeltop=False, length=small, width=1)
+        ax1.tick_params(axis="y", direction='in', right=True, labelright=False, length=small, width=1)
+        ax1.scatter(x, (y - y_hat)/(y_hat_err), color="black", s=1) 
+        ax1.tick_params(axis="x", labelsize=ticksize)
+        ax1.tick_params(axis="y", labelsize=ticksize)
+        ax1.set_ylabel("Pull", fontsize=labelsize, loc='top')
+        ax1.set_xlabel(xlabel, fontsize=labelsize, loc='right')
+        ax1.set_ylim(residual_limits)
+        ax1.set_yticks(np.arange(int(np.floor(np.min(residual_limits))), int(np.ceil(np.max(residual_limits)))))
+        ax1.axhline(y=-1.0, color=dusty, lw=2, linestyle="--")
+        ax1.axhline(y=1.0, color=dusty, lw=2, linestyle="--")
+        ax1.tick_params(axis="x", labelsize=ticksize)
+        ax1.tick_params(axis="y", labelsize=ticksize)
+        ax1.grid(which="both", axis="both", linewidth=0.2, color="grey", linestyle=":")
+        ax1.locator_params(axis='y', nbins=6)
+
+        fig.tight_layout()
+        if(not plot_in_bins):
+            mf = ScalarFormatter(useMathText=True)
+            mf.set_powerlimits((-2,2))
+            plt.gca().xaxis.set_major_formatter(mf)
+            offset = ax1.xaxis.get_offset_text()
+            offset.set_size(int(0.8*ticksize))
+
+        fig.subplots_adjust(hspace=.0)
+        if(save_directory is not None):
+            print("Saving figure to {0}...".format(save_directory))
+            fig.savefig(save_directory)
+        if(display):
+            plt.pause(0.01)
+            fig.show()
+        else:
+            plt.close(fig)        
         
     ########################################################################################################
     ## Fit(prefit=False)
diff --git a/user/sumlab_auto.py b/user/sumlab_auto.py
index 1db86beeb75273412ae1459185c11507db4a7459..2f71407ed10e9c8f3bbb2c73fcc5a95cee243802 100644
--- a/user/sumlab_auto.py
+++ b/user/sumlab_auto.py
@@ -121,6 +121,8 @@ for i, (file, path) in enumerate(files_to_fit):
             fit_out["d_{:s}".format(key)] = value
         f_data.PlotFit(plot_in_bins=True, display=False,
                        save_directory=f"{folder}/{file[:-4]}_fit.png")
+        f_data.PlotFitSumlab(plot_in_bins=True, display=False,
+                       save_directory=f"{folder}/{file[:-4]}_fit_sumlab.png")
 
     df = pd.DataFrame.from_dict([fit_out])
     df.to_csv("{}/fit_results_{:s}.csv".format(folder, file[:-4]))