diff --git a/examples/findPeaks.py b/examples/findPeaks.py
index be3b4a339e975fff1adb5e1e57aaad3cf7433919..8333da8cf36a97e5c3538e29909206beeb72e77e 100755
--- a/examples/findPeaks.py
+++ b/examples/findPeaks.py
@@ -1,5 +1,12 @@
 #!/usr/bin/env python3
 """findPeaks: find and fit peaks in gamma spectrum
+
+   algorithm: perform a tow-stage determination of peak posisions:
+
+   1. scipi.signal.find_peaks() is run on a smoothed histogram
+      of channel counts.
+   1. the output is used as a starting point for histogram fits
+      with kafe2 based on the binned neg-log-likelhood cost function.
 """
 
 import argparse
@@ -163,17 +170,28 @@ ax1.set_xlabel('Channel #')
 ax1.grid(linestyle='dotted', which='both')
 
 if plot:
+    # show spectrum and result of find_peaks
     xhst = np.linspace(0, len(hst), len(hst), endpoint=False) +0.5
-    ax0.errorbar(xhst, hst, yerr=np.sqrt(hst), fmt='.', color='grey', markersize=2, linewidth=2)
-    ax1.errorbar(xhst, hst, yerr=np.sqrt(hst), fmt='.', color='grey', alpha=0.1)
-    colors = ('cadetblue', 'orange', 'olive', 'orchid', 'turquoise', 'tomato', 'green', 'pink', 'salmon', 'yellowgreen' )
+    ax0.plot(xhst, hst_s, 'b-', linewidth=3, label = 'smoothed channel counts')
+    ax0.errorbar(xhst, hst, yerr=np.sqrt(hst), label='channel counts',
+                 fmt='.', color='grey', markersize=2, linewidth=2, alpha=0.5)
+    ax0.plot(peaks, hst_s[peaks], '*', color='red', markersize=10, label='result of find_peaks()')
+    ax0.legend(loc = 'best')
+
+    # show fitted peaks
+    ax1.errorbar(xhst, hst, yerr=np.sqrt(hst), fmt='.', color='grey', alpha=0.1,
+                 label = 'channel counts')
+    # select colors for peaks 
+    pcolors = ('cadetblue', 'orange', 'olive', 'orchid', 'turquoise',
+              'tomato', 'green', 'pink', 'salmon', 'yellowgreen' )
     for i, fit in enumerate(fit_results):
         # plot fitted peak in fit range
-        colr = colors[i%10]
+        colr = pcolors[i%10]
         xplt = np.linspace(plot_ranges[i][0], plot_ranges[i][1],
                     10*int((plot_ranges[i][1]-plot_ranges[i][0])))
         ax1.plot(xplt, gauss_plus_bkg(xplt, *fit.parameter_values),
-             linestyle = 'solid', linewidth=3, color= colr)
+                 linestyle = 'solid', linewidth=3, color= colr,
+                 label='peak ' + str(i+1) + '@' + str(int(10*fit.parameter_values[1])/10.) )
         # plot fitted peak near fit region
         xplt2 = np.linspace(plot_ranges2[i][0], plot_ranges2[i][1],
                     10*int((plot_ranges[i][1]-plot_ranges[i][0])))
@@ -190,7 +208,8 @@ if plot:
         ax1.vlines(mu, 0, mx - h, 
                linewidth=1, linestyle='dashed', color='goldenrod')
         ax1.hlines(mx-h/2, mu-fwhm/2, mu+fwhm/2,
-               linewidth=2, color='gold')    
+               linewidth=2, color='gold')
+        ax1.legend(loc = 'best')
 
     if kafe2_plots:
     # kafe2 plots