diff --git a/eventdisplay_anasum/spectral_analysis.C b/eventdisplay_anasum/spectral_analysis.C
index 1a79fd6..1885557 100644
--- a/eventdisplay_anasum/spectral_analysis.C
+++ b/eventdisplay_anasum/spectral_analysis.C
@@ -119,6 +119,17 @@ void spectral_analysis(
         fFitFunction->Draw("same");
     }
     a.plotFitValues();
+    if (c)
+    {
+        c->Print((output_file + ".pdf").c_str());
+    }
+
+    TCanvas *cCounting = a.plotCountingHistograms();
+    if (cCounting)
+    {
+        cCounting->Print((output_file + "_counting.pdf").c_str());
+    }
+
 
     // write fit results as yaml file
     string yaml_file = output_file + ".yaml";
@@ -145,11 +156,6 @@ void spectral_analysis(
     }
     yaml_out.close();
 
-    if (c)
-    {
-        c->Print((output_file + ".pdf").c_str());
-    }
-
     // write spectral points
     a.writeSpectralPointsToCSVFile((output_file + ".ecsv").c_str());
 }
diff --git a/eventdisplay_anasum/spectral_analysis_config.txt b/eventdisplay_anasum/spectral_analysis_config.txt
index 749c916..6188f2d 100644
--- a/eventdisplay_anasum/spectral_analysis_config.txt
+++ b/eventdisplay_anasum/spectral_analysis_config.txt
@@ -6,7 +6,7 @@ ENERGY_MIN: 0.09
 ENERGY_MAX: 100.
 SIGNIFICANCE_FIT_MIN: 0.0
 EXCESS_EVENTS_FIT_MIN: 0.0
-SIGNIFICANCE_PLOT_MIN: 1.0
-EXCESS_EVENTS_PLOT_MIN: 1.0
+SIGNIFICANCE_PLOT_MIN: 0.5
+EXCESS_EVENTS_PLOT_MIN: 0.5
 PLOT_UPPER_LIMITS: 1
 PLOT_EVENT_NUMBERS: 1
diff --git a/v2dl5/light_curves/binary_plotting.py b/v2dl5/light_curves/binary_plotting.py
index 02fe8d2..8f35a4a 100644
--- a/v2dl5/light_curves/binary_plotting.py
+++ b/v2dl5/light_curves/binary_plotting.py
@@ -555,339 +555,3 @@ def plot_distribution(
             file_type=file_type,
             figure_dir=figure_dir,
         )
-
-
-
-
-
-# Temporary stuff - probably not needed
-#
-# def plotLightCurve_fluxvsPhase_inOrbits(
-#         fDataDict,
-#         PlotInstruments,
-#         F,
-#         lc_spline_bin_centers,
-#         lc_spline_sv,
-#         lc_spline_sv_err,
-#         orbital_period_BL=315.,
-#         orbital_periodBins=20,
-#         plot_variable=None):
-#     """plot flux vs orbital phase (separated in orbits, all in one plot)
-#     """
-#
-#     MJD_firstOrbit, N_orbits = getNumberOfOrbits(
-#         fDataDict, PlotInstruments,
-#         orbital_period_BL, False)
-#
-#     lightCurvePlottingUtilities.paper_figures(4, 4)
-#     colors = lightCurvePlottingUtilities.getColorList(N_orbits)
-#     if N_orbits < 6:
-#         markers = lightCurvePlottingUtilities.getMarkerList()
-#     else:
-#         markers = ['o']*N_orbits
-#
-#     if len(PlotInstruments) < 1:
-#         return
-#
-#     # plot average and interpolated light curves
-#     glabel = PlotInstruments[0] + " (average)"
-#     if lc_spline_bin_centers and len(lc_spline_bin_centers) > 0:
-#         plt.plot(
-#             lc_spline_bin_centers,
-#             lc_spline_sv,
-#             color='tab:gray',
-#             linestyle='--',
-#             linewidth=lightCurvePlottingUtilities.getLineWidth())
-#         plt.fill_between(lc_spline_bin_centers,
-#                         lc_spline_sv - lc_spline_sv_err,
-#                         lc_spline_sv + lc_spline_sv_err,
-#                         color='tab:gray', linestyle='--',
-#                         linewidth=lightCurvePlottingUtilities.getLineWidth(),
-#                         alpha=0.3)
-#
-#     # plot one light curve per orbital period
-#     for i in range(0, N_orbits):
-#         x = []
-#         y = []
-#         ex = []
-#         mj = []
-#         orbit_min = MJD_firstOrbit + i * orbital_period_BL
-#         orbit_max = MJD_firstOrbit + (i + 1) * orbital_period_BL - 1.
-#         for j in range(len(PlotInstruments)):
-#             i_plotValue, i_plotError = lightCurvePlottingUtilities.get_plotting_variable(
-#                 plot_variable, j)
-#
-#             for p in range(len(fDataDict[PlotInstruments[j]]['phaseN'])):
-#                 if fDataDict[PlotInstruments[j]]['MJD'][p] >= orbit_min \
-#                         and fDataDict[PlotInstruments[j]]['MJD'][p] < orbit_max + 1.:
-#
-#                     x.append(fDataDict[PlotInstruments[j]]['phase'][p])
-#                     y.append(fDataDict[PlotInstruments[j]][i_plotValue][p])
-#                     ex.append(fDataDict[PlotInstruments[j]][i_plotError][p])
-#                     mj.append(fDataDict[PlotInstruments[j]]['MJD'][p])
-#
-#         if len(mj):
-#             OrbitPhStr = "MJD %d - %d" % (
-#                 orbit_min, orbit_max)
-#             plt.errorbar(
-#                 x,
-#                 y,
-#                 ex,
-#                 None,
-#                 color=colors[i],
-#                 marker=markers[i],
-#                 linestyle='none',
-#                 label=OrbitPhStr,
-#                 linewidth=lightCurvePlottingUtilities.getLineWidth(),
-#                 markersize=lightCurvePlottingUtilities.getMarkerSize())
-#
-#     plt.xlabel(
-#         lightCurvePlottingUtilities.get_orbital_phase_axis_string(orbital_period_BL) )
-#     plt.ylabel(
-#         lightCurvePlottingUtilities.getFluxAxisString(
-#         PlotInstruments[0],plot_variable) )
-#     if PlotInstruments[0].find('Optical') < 0:
-#         plt.axhline(y=0, linestyle=':')
-#     if getPrintInstrumentName(PlotInstruments).find( "fwhm" ) < 0 and \
-# getPrintInstrumentName(PlotInstruments).find( "ew" ) < 0:
-#         plt.legend(prop={'size': 10}, framealpha=0.1)
-#
-#     lightCurvePlottingUtilities.printFigure(
-#         getPrintInstrumentName(PlotInstruments) +
-#         "-HESSJ0632p057-LC-phaseFolded-%dd-Orbits" %
-#         orbital_period_BL)
-#
-#
-#
-# def plotAverageLightCurve_fluxvsPhase(
-#         fDataDict,
-#         PlotInstruments,
-#         orbital_period_BL=315.,
-#         orbital_periodBins=20):
-#     """
-#     plot average flux vs orbital phase
-#
-#     calculates also average light curves
-#     """
-#
-#     print("Plot phase binned averaged light curve:")
-#     print("\t orbital phase (%.1f d)" % orbital_period_BL)
-#     print("\t number of phase bins (%d)" % orbital_periodBins)
-#     print("\t calculating average light curve for ", PlotInstruments)
-#     if len(PlotInstruments) < 1:
-#         return
-#
-#     # copy all data into one set of arrays
-#     # averaged light curve is calculated from
-#     # light curve bins of all data
-#     i_MJD = []
-#     i_flux = []
-#     i_flux_err = []
-#     for I in PlotInstruments:
-#         i_MJD.extend(fDataDict[I]['MJD'])
-#         i_flux.extend(fDataDict[I]['flux'])
-#         i_flux_err.extend(fDataDict[I]['flux_err'])
-#
-#     lightCurvePlottingUtilities.paper_figures(4, 4)
-#
-#     # calculate phase binned average light curve
-#     lc_bincenters, lc_mean, lc_std = lightCurveAverageing.calculateAverageLightCurve(
-#         i_MJD, i_flux, i_flux_err, orbital_period_BL, orbital_periodBins)
-#
-#     # bin width
-#     lc_binw = np.repeat(0.5 / orbital_periodBins, len(lc_bincenters))
-#     plt.errorbar(lc_bincenters, lc_mean, lc_std, lc_binw,
-#                  color='b', linestyle='none',
-#                  linewidth=lightCurvePlottingUtilities.getLineWidth(),
-#                  markersize=lightCurvePlottingUtilities.getMarkerSize(),
-#                  marker='o', label='average')
-#
-#     # calculate cubic spline smoothed average light curve
-#     F, lc_spline_bin_centers, lc_spline_sv = \
-#         lightCurveAverageing.smoothCubeSplineLightCurve(
-#             i_MJD, i_flux, i_flux_err,
-#             orbital_period_BL, orbital_periodBins,
-#             True, 0.)
-#     # get errors on the same
-#     F_temp, lc_spline_bin_centers, lc_spline_sv_err = \
-#         lightCurveAverageing.smoothCubeSplineLightCurve(
-#             i_MJD, i_flux, i_flux_err,
-#             orbital_period_BL, orbital_periodBins,
-#             True, 1.)
-#     lc_spline_sv_err = list(
-#         np.array(lc_spline_sv_err) -
-#         np.array(lc_spline_sv))
-#
-#     plt.plot(lc_spline_bin_centers, lc_spline_sv, color='r', linestyle='--',
-#              linewidth=lightCurvePlottingUtilities.getLineWidth(),
-#              label='spline average')
-#
-#     # data set used for averaging
-#     wP = []
-#     wP = [lightCurveAnalysisorbital_period.getOrbitalPhase(
-#         x, orbital_period_BL) for x in i_MJD]
-#     plt.errorbar(wP, i_flux, i_flux_err, None,
-#                  color='g', linestyle='none',
-#                  linewidth=lightCurvePlottingUtilities.getLineWidth(),
-#                  markersize=lightCurvePlottingUtilities.getMarkerSize(),
-#                  marker='o', label='average')
-#
-#     plt.fill_between(lc_spline_bin_centers,
-#                      lc_spline_sv - lc_spline_sv_err,
-#                      lc_spline_sv + lc_spline_sv_err,
-#                      color='r', linestyle='--',
-#                      linewidth=lightCurvePlottingUtilities.getLineWidth(),
-#                      alpha=0.3)
-#
-#     plt.xlabel(
-#         lightCurvePlottingUtilities.get_orbital_phase_axis_string(orbital_period_BL))
-#     plt.ylabel(lightCurvePlottingUtilities.getFluxAxisString(PlotInstruments[0]))
-#     plt.legend()
-#     plt.axhline(y=lc_spline_sv[0], linestyle=':')
-#
-#     lightCurvePlottingUtilities.printFigure(
-#         getPrintInstrumentName(PlotInstruments) +
-#         "-HESSJ0632p057-LC-phaseFolded-%dd-Average" %
-#         orbital_period_BL)
-#
-#     return F, lc_spline_bin_centers, lc_spline_sv, lc_spline_sv_err
-#
-#
-# def plotLightCurve_fluxvsMJD_XandGray(
-#         fDataDict, mjd_min, mjd_max,
-#         icrc2019Plots=False,
-#         yaxis_min = -0.9, yaxis_max=7.5,
-#         fColorDict=None,
-#         plot_title=None,
-#         convert_erg=False):
-#     """plot flux vs MJD with two different y-axis
-#
-#     note: fixed limits for the y-axis
-#     """
-#
-#     lightCurvePlottingUtilities.paper_figures(4, 4, 1, False)
-#     # quick and dirty fix to get consistent colors and markers
-#     if len(fColorDict) == 5:
-#         colors = lightCurvePlottingUtilities.getColorList(3)
-#         colors[-1] = 'black'
-#         colors.append( 'gray' )
-#     else:
-#         colors = lightCurvePlottingUtilities.getColorList(len(fColorDict))
-#     markers = lightCurvePlottingUtilities.getMarkerList(icrc2019Plots)
-#     markers[3]='+'
-#     markers[4]='x'
-#
-#     fig, ax1 = plt.subplots()
-#     ax1_y = ax1.twinx()
-#     ax1_x = ax1.twiny()
-#     if convert_erg:
-#         ax1.set_ylim(ymin=yaxis_min,ymax=yaxis_max*1.8)
-#     else:
-#         ax1.set_ylim(ymin=yaxis_min,ymax=yaxis_max)
-#     ax1_y.set_ylim(ymin=yaxis_min,ymax=yaxis_max)
-#
-#     c = 0
-#     for key, fData in fDataDict.items():
-#         print('Plotting %s (%d data points)' % (key, len(fData)))
-#
-#         # scale all data by 1.e-12
-#         # (scale is added to the legend text)
-#         y = np.asarray(fData['flux']) / 1.e-12
-#         y_err = np.asarray(fData['flux_err']) / 1.e-12
-#
-#         if key in fColorDict:
-#             c = fColorDict[key]
-#
-#         if len(y) == 0:
-#             c += 1
-#             continue
-#
-#         # left labelled data
-#         if "VERITAS" in key or \
-#                 "HESS" in key or \
-#                 "MAGIC" in key:
-#             ln1 = ax1.errorbar(
-#                 fData['MJD'],
-#                 y,
-#                 y_err,
-#                 fData['MJD_err'],
-#                 color=colors[c],
-#                 marker=markers[c],
-#                 label=key,
-#                 linestyle='none',
-#                 fillstyle='full',
-#                 linewidth=lightCurvePlottingUtilities.getLineWidth(),
-#                 markersize=lightCurvePlottingUtilities.getMarkerSize())
-#             ax1_y.plot(np.nan, '-r',
-#                      color=colors[c], marker=markers[c], label=key,
-#                      linestyle='none', fillstyle='full',
-#                      linewidth=lightCurvePlottingUtilities.getLineWidth(),
-#                      markersize=lightCurvePlottingUtilities.getMarkerSize())
-#         # right labelled data
-#         else:
-#             pLabel=key
-#             if key.find('XRT')>=0:
-#                 pLabel="$\it{Swift}$-XRT"
-#             ln2 = ax1_y.errorbar(
-#                 fData['MJD'],
-#                 y,
-#                 y_err,
-#                 fData['MJD_err'],
-#                 color=colors[c],
-#                 marker=markers[c],
-#                 label=pLabel,
-#                 linestyle='none',
-#                 fillstyle='full',
-#                 linewidth=lightCurvePlottingUtilities.getLineWidth(),
-#                 markersize=lightCurvePlottingUtilities.getMarkerSize())
-#         c += 1
-#
-#     ax1.axis(xmin=mjd_min,xmax=mjd_max)
-#     ax1.locator_params(axis='x', tight=True, nbins=4)
-#     ax1.set_xlabel('Modified Julian Day (MJD)')
-#     # left label: assume always gamma rays
-#     ax1.set_ylabel(lightCurvePlottingUtilities.getFluxAxisString("VERITAS",
-#                                                                  None,
-#                                                                  "10$^{-12} \\times$",
-#                                                                  convert_erg))
-#
-#     if plot_title:
-#         ty=0.9*yaxis_max
-#         if convert_erg:
-#             ty*=1.8
-#         ax1.text(
-#               mjd_min+0.75*(mjd_max-mjd_min),
-#               ty,
-#               plot_title,
-#               fontsize=10,
-#               weight='bold',
-#               )
-#
-#     # right label: assume always X-rays
-#     if 'Swift XRT' in fDataDict or 'NuSTAR' in fDataDict:
-#         ax1_y.set_ylabel(lightCurvePlottingUtilities.getFluxAxisString(
-#   "Swift XRT", None, "10$^{-12} \\times$"))
-#
-#     ax1_y.legend(loc=2)
-#     if yaxis_min < 0.:
-#         plt.axhline(y=0, linestyle=':')
-#     lightCurvePlottingUtilities.align_yaxis(ax1, 0, ax1_y, 0)
-#
-#     # orbital phase axis
-#     # (note: use default orbital period here
-#     # as defined in lightCurveAnalysisorbital_period
-#     phase_ticks = []
-#     label_format = '%.2f'
-#     for t in ax1_x.get_xticks():
-#         mjd=mjd_min+t*(mjd_max-mjd_min)
-#         phase_ticks.append(label_format % (
-#             lightCurveAnalysisorbital_period.getOrbitalPhase(mjd),))
-#     # remove first and last label from plotting
-#     if len(phase_ticks)>0:
-#         phase_ticks[0]=None
-#         phase_ticks[-1]=None
-#     ax1_x.set_xticklabels(phase_ticks)
-#     ax1_x.set_xlabel('orbital phase')
-#
-#     lightCurvePlottingUtilities.printFigure('XG-HESSJ0632p057-LC')
-#
diff --git a/v2dl5/light_curves/data_reader.py b/v2dl5/light_curves/data_reader.py
index e21b9aa..dfcd68a 100644
--- a/v2dl5/light_curves/data_reader.py
+++ b/v2dl5/light_curves/data_reader.py
@@ -8,6 +8,7 @@
 import yaml
 from astropy.table import Table
 
+import v2dl5.binaries as binaries
 import v2dl5.orbital_phase as orbit
 
 
@@ -58,7 +59,7 @@ def read_data(self):
         """
         for data_config in self.config:
             self.data_dict[data_config["instrument"]] = self._read_fluxes_from_file(data_config)
-            self._add_orbital_parameters(self.data_dict[data_config["instrument"]])
+            self._add_orbital_parameters(self.data_dict[data_config["instrument"]], data_config)
 
     def _read_fluxes_from_file(self, data_config):
         """Read flux from file."""
@@ -78,19 +79,56 @@ def _read_fluxes_from_file(self, data_config):
 
         return None
 
-    def _add_orbital_parameters(self, data):
+    def _apply_phase_mask(self, data, data_config):
         """
-        Add orbital phase to light-curve data data.
+        Apply phase mask to data.
+
+        Requires 'phase_cut_binary' to be set in data_config.
+        """
+        if not data_config.get("phase_cut_binary"):
+            return data
+
+        try:
+            _binary = binaries.binary_properties()[data_config["phase_cut_binary"]]
+        except KeyError:
+            raise KeyError(f"Binary {data_config['phase_cut_binary']} not found in binaries.py")
+
+        phase_min = data_config.get("phase_min", 0.0)
+        phase_max = data_config.get("phase_max", 1.0)
+
+        phases = [
+            orbit.get_orbital_phase(
+                mjd,
+                orbital_period=_binary["orbital_period"],
+                mjd_0=_binary["mjd_0"],
+                phase_reduce=True
+            )
+            for mjd in data["MJD"]
+        ]
+        phase_mask = [
+            (p >= phase_min) & (p <= phase_max) for p in phases
+        ]
+        for key in data.keys():
+            data[key] = [val for val, mask in zip(data[key], phase_mask) if mask]
+
+        return data
+
+    def _add_orbital_parameters(self, data, data_config):
+        """
+        Add orbital phase to light-curve data and filter by phase range.
 
         Parameters
         ----------
         data: dict
             Data dictionary with light-curve data.
-
+        data_config: dict
+            Data configuration dictionary.
         """
         orbital_period = self.binary["orbital_period"]
         mjd_0 = self.binary["mjd_0"]
 
+        self._apply_phase_mask(data, data_config)
+
         data["phase"] = [
             orbit.get_orbital_phase(
                 mjd, orbital_period=orbital_period, mjd_0=mjd_0, phase_reduce=True
diff --git a/v2dl5/plotting/utilities.py b/v2dl5/plotting/utilities.py
index 15670c6..bcf0c37 100644
--- a/v2dl5/plotting/utilities.py
+++ b/v2dl5/plotting/utilities.py
@@ -111,38 +111,6 @@ def print_figure(print_name, file_type=".pdf", figure_dir="./figures/"):
         plt.savefig(figure_file)
 
 
-def get_orbital_phase_axis_string(orbital_period):
-    """Return a string for the orbital phase axis."""
-    # strOr = "orbital phase (%.1f d)" % orbital_period
-    return "orbital phase"
-
-
-def get_flux_axis_string(instrument, plot_variable=None, scale_factor=None, energy_flux=False):
-    """Flux axis string with the correct unit depending on the type of instrument."""
-    if not scale_factor:
-        scale_factor = ""
-
-    if "VERITAS" in instrument or "HESS" in instrument or "MAGIC" in instrument:
-        if energy_flux:
-            pAxisFluxSTR = "Flux E>350 GeV (" + scale_factor + "erg/(cm$^2$ s))"
-        else:
-            pAxisFluxSTR = "Flux E>350 GeV (" + scale_factor + "1/(cm$^2$s))"
-    elif instrument == "XRT":
-        pAxisFluxSTR = "Swift rate (0.3-10 keV) (counts/s)"
-    elif instrument.find("Optical") >= 0:
-        if plot_variable and len(plot_variable) == 1:
-            if plot_variable[0] == "ew":
-                pAxisFluxSTR = "EW ($\\AA$)"
-            elif plot_variable[0] == "vc":
-                pAxisFluxSTR = "v$_c$ (km/s)"
-            elif plot_variable[0] == "fwhm":
-                pAxisFluxSTR = "FWHM (km/s)"
-    else:
-        pAxisFluxSTR = "Flux 0.3-10 keV (" + scale_factor + "erg/(cm$^2$s))"
-
-    return pAxisFluxSTR
-
-
 def get_paper_figures_parameters(width=None, height=None, xtick_top=True, legend_font_size=8):
     """Return default paper figure parameters."""
     if width is None: