diff --git a/python/bidi.py b/python/bidi.py
index 0070778..3e15fad 100644
--- a/python/bidi.py
+++ b/python/bidi.py
@@ -11,7 +11,7 @@ def __init__(self, channels, names):
         # works only from python 2.7 
         # self.reverse = {v: k for k, v in self.dictionary.iteritems()}
         self.reverse = {}
-        for k,v in self.dictionary.iteritems():
+        for k,v in self.dictionary.items():
             self.reverse[v] = k
         #
         self.dictionary.update(self.reverse)
diff --git a/python/grammar.py b/python/grammar.py
index 2bc04d2..a6bdd3f 100644
--- a/python/grammar.py
+++ b/python/grammar.py
@@ -3,7 +3,7 @@
 # (sequencer with variables)
 #
 # Author: Laurent Le Guillou
-
+from __future__ import print_function
 import re
 
 class SeqParser(object):
@@ -1031,7 +1031,7 @@ def m_func_clocks_names_line(self, pos):
     
         pnext = self.m_zmsp(pnext)
         if self.verbose:
-            print pnext
+            print(pnext)
 
         comment = ''
         r = self.m_comment(pnext)
@@ -2034,7 +2034,7 @@ def m_seq(self, pos):
         pnext = self.m_empty_lines(pnext)
 
         if self.verbose:
-            print "[includes]"
+            print("[includes]")
         r = self.m_include_section(pnext)
         if r == None:
             result['includes'] = []  # keep it optional
@@ -2044,7 +2044,7 @@ def m_seq(self, pos):
         pnext = self.m_empty_lines(pnext)
 
         if self.verbose:
-            print "[constants]"
+            print("[constants]")
         r = self.m_constant_section(pnext)
         if r == None:
             return None
@@ -2053,7 +2053,7 @@ def m_seq(self, pos):
         pnext = self.m_empty_lines(pnext)
 
         if self.verbose:
-            print "[clocks]"
+            print("[clocks]")
         r = self.m_clock_section(pnext)
         if r == None:
             return None
@@ -2062,7 +2062,7 @@ def m_seq(self, pos):
         pnext = self.m_empty_lines(pnext)
 
         if self.verbose:
-            print "[pointers]"
+            print("[pointers]")
         r = self.m_ptr_section(pnext)
         if r != None:
             pnext, result['pointers'] = r
@@ -2070,7 +2070,7 @@ def m_seq(self, pos):
         pnext = self.m_empty_lines(pnext)
 
         if self.verbose:
-            print "[functions]"
+            print("[functions]")
         r = self.m_func_section(pnext)
         if r == None:
             return None
@@ -2079,7 +2079,7 @@ def m_seq(self, pos):
         pnext = self.m_empty_lines(pnext)
 
         if self.verbose:
-            print "[subroutines]"
+            print("[subroutines]")
         r = self.m_subr_section(pnext)
         if r != None:
             pnext, result['subroutines'] = r
@@ -2087,7 +2087,7 @@ def m_seq(self, pos):
         pnext = self.m_empty_lines(pnext)
             
         if self.verbose:
-            print "[mains]"
+            print("[mains]")
         r = self.m_subr_section(pnext, main=True)
         if r == None:
             return None
diff --git a/python/multiscope.py b/python/multiscope.py
index c6ba601..a7437c9 100644
--- a/python/multiscope.py
+++ b/python/multiscope.py
@@ -1,4 +1,4 @@
-# !/usr/bin/env python
+#!/usr/bin/env python
 
 # Displays a scope-like view of the CCD output waveforms for the whole raft.
 #
@@ -12,15 +12,15 @@
 #
 # Syntax in a script:
 # import multiscope
-# multiscope.raft_display("tm-scan.fits")
-
+# multiscope.raft_display_allchans("tm-scan.fits")
+from __future__ import print_function
 import sys
 import os.path
 import numpy as np
 from matplotlib import pyplot as plt
 from matplotlib import colors as mplcol
 import scope
-import pyfits
+from astropy.io import fits
 
 
 def get_scandata_raft(inputfile, datadir=''):
@@ -33,7 +33,7 @@ def get_scandata_raft(inputfile, datadir=''):
     raftarrays = []
     if os.path.splitext(inputfile)[1] in [".fits", ".fz"]:
         # starts with 00 through 22
-        seglist = ["%d%d" % (i, j) for i in range(3) for j in range(3)]
+        seglist = slot_ids()
         # when REB2 data is missing
         #seglist = ["%d%d" % (i, j) for i in range(2) for j in range(3)]
 
@@ -58,7 +58,7 @@ def get_scandata_raft(inputfile, datadir=''):
             #print fullreb.shape
             raftarrays.extend([a for a in np.split(fullreb, 3, axis=0)])  # splits REB data into 3 CCDs
     elif inputfile == '':
-        seglist = ["%d%d" % (i, j) for i in range(3) for j in range(3)]
+        seglist = slot_ids()
         for segstr in seglist:
             d = os.path.join(datadir, 'S%s' % segstr)
             for f in os.listdir(d):
@@ -72,6 +72,7 @@ def get_scandata_raft(inputfile, datadir=''):
 
     return raftarrays, seglist
 
+
 def raft_display_allchans(inputfile, datadir='', suptitle=''):
     """
     Builds up data from all raft files and display scans.
@@ -81,7 +82,36 @@ def raft_display_allchans(inputfile, datadir='', suptitle=''):
     :return:
     """
     raftarrays, seglist = get_scandata_raft(inputfile, datadir)
+    if inputfile:
+        dataname = scope.get_rootfile(inputfile)
+    else:
+        dataname = os.path.split(datadir)[-1]
+    if not suptitle:
+        suptitle = dataname
+    outfile = os.path.join(datadir, "multiscope-%s.png" % dataname)
+    plot_raft_allchans(raftarrays, seglist, suptitle)
+    plt.savefig(outfile)
+    plt.show()
+
+
+def slot_ids():
+    return ["%d%d" % (i, j) for i in range(3) for j in range(3)]
+
+
+def plot_raft_allchans(raftarrays, seglist, suptitle=''):
+    """
+    Builds up data from all raft files and display scans.
 
+    Parameters
+    ----------
+    raftarrays: list
+        list of channel data for a single CCD.
+    seglist: list
+        list of CCD slot ids, e.g.,
+        ['00', '01', '02', '10', '11', '12', '20', '21', '22']
+    suptitle: str ['']
+        Title of the plot.
+    """
     fig, axes = plt.subplots(nrows = 3, ncols = 3, figsize=(15, 9))
     # when REB2 data is missing
     # fig, axes = plt.subplots(nrows = 2, ncols = 3, figsize=(15, 9))
@@ -91,7 +121,7 @@ def raft_display_allchans(inputfile, datadir='', suptitle=''):
     maxplot = 0
     minplot = 25000
     for num,tmscope in enumerate(raftarrays):
-        ax = axes[num  / 3, num  % 3 ]
+        ax = axes[num // 3, num  % 3 ]
 
         # single CCD plot
         for c in range(16):
@@ -109,30 +139,19 @@ def raft_display_allchans(inputfile, datadir='', suptitle=''):
 
         if num%3 == 0:
             ax.set_ylabel('Scan (ADU)')
-        if num/3 == 2:
+        if num//3 == 2:
             ax.set_xlabel('Time increment (10 ns)')
         ax.grid(True)
 
     for num in range(len(raftarrays)):
-        ax = axes[num / 3, num % 3]
+        ax = axes[num // 3, num % 3]
         # common scale for all subplots
         ax.set_ylim(0, maxplot * 1.02)
         #ax.set_ylim(10000, 40000)
         #ax.set_ylim(minplot * 0.95, maxplot * 1.02)
 
     ax.legend(bbox_to_anchor=(1.05, 0), loc='lower left', borderaxespad=0.)
-
-    if inputfile:
-        dataname = scope.get_rootfile(inputfile)
-    else:
-        dataname = os.path.split(datadir)[-1]
-
-    if suptitle:
-        plt.suptitle(suptitle, fontsize='x-large')
-    else:
-        plt.suptitle(dataname)
-    plt.savefig(os.path.join(datadir, "multiscope-%s.png" % dataname))
-    plt.show()
+    plt.suptitle(suptitle, fontsize='x-large')
 
 
 def raft_display_longscan(inputfile, niter, datadir='', suptitle=''):
@@ -153,7 +172,7 @@ def raft_display_longscan(inputfile, niter, datadir='', suptitle=''):
 
     # plot all channels, with one subplot per CCD
     for num in range(len(raftfits)):
-        ax = axes[num  / 3, num  % 3 ]
+        ax = axes[num // 3, num  % 3 ]
         tmscope = scope.stitch_long_scan(raftfits[num], niter, datadir)
 
         # single CCD plot
@@ -228,19 +247,19 @@ def reframe_scope_frames(listfits):
     :return:
     """
     for f in listfits:
-        h = pyfits.open(f)
+        h = fits.open(f)
         imgcols = h[1].header['naxis1']
         imglines = h[1].header['naxis2']
         detstring = '[1:%d,1:%d]' % (imgcols * 8, 2 * imglines)
         # Create HDU list
-        primaryhdu = pyfits.PrimaryHDU()
+        primaryhdu = fits.PrimaryHDU()
         primaryhdu.header = h[0].header.copy()
         primaryhdu.header['DETSIZE'] = (detstring, 'NOAO MOSAIC keywords')
-        hcopy = pyfits.HDUList([primaryhdu])
+        hcopy = fits.HDUList([primaryhdu])
 
         #print h.info()
         for extnum in range(1,17,1):
-            exthdu = pyfits.ImageHDU(data=h[extnum].data)
+            exthdu = fits.ImageHDU(data=h[extnum].data)
             get_extension_header(imgcols, imgcols, extnum-1, exthdu, detstring)
             hcopy.append(exthdu)
 
@@ -260,13 +279,13 @@ def corrcoef_raftscope(raftsfits, ROIrows, ROIcols, norm=True):
     stackh = []
 
     for fl in raftsfits:
-        h = pyfits.open(fl)
+        h = fits.open(fl)
         for i in range(1, 17):
             stackh.append(h[i].data[ROIrows, ROIcols])
         h.close()
         del h
     stackh = np.stack(stackh)
-    print stackh.shape
+    print(stackh.shape)
 
     a = []
 
@@ -276,7 +295,7 @@ def corrcoef_raftscope(raftsfits, ROIrows, ROIcols, norm=True):
         else:
             a.append(np.cov(stackh[:, :, numcol]))
     a = np.stack(a)
-    print a.shape
+    print(a.shape)
 
     return a.mean(axis=0)
 
diff --git a/python/rebtxt.py b/python/rebtxt.py
index ca12768..a1f1d7b 100644
--- a/python/rebtxt.py
+++ b/python/rebtxt.py
@@ -4,7 +4,7 @@
 # Python minimal interface for the REB FPGA
 # TXT assembler-style IO
 #
-
+from __future__ import print_function
 from sequencer import *
 import grammar
 
@@ -148,7 +148,7 @@ def parse_functions(self, functions_node):
             self.functions_desc[name] = dict(func_dict)
 
             if self.verbose:
-                print name, fullname
+                print(name, fullname)
 
             function = Function(name=name,
                                 fullname=fullname,
@@ -173,7 +173,7 @@ def parse_functions(self, functions_node):
                 cpos += 1
 
             if self.verbose:
-                print channel_position
+                print(channel_position)
             self.functions_desc[name]['clocks'] = func['clocks']  # list of active channels
 
             # self.timelengths = {0: 12, 1: 14}
@@ -208,7 +208,7 @@ def parse_functions(self, functions_node):
                         cval = values[cpos]
                         output |= (cval << crank)
                 if self.verbose:
-                    print bin(output)
+                    print(bin(output))
 
                 outputs[islice] = output
 
@@ -331,7 +331,7 @@ def parse_subroutine(self, sub_node):
         subname = sub_node['name']
         fullname = sub_node['comment']
         if self.verbose:
-            print "   name = ", subname
+            print("   name = ", subname)
             # print "   fullname = ", fullname
 
         sub = Subroutine()
@@ -385,12 +385,12 @@ def parse_result(self, result):
         # Parse all subroutines
         self.parse_subroutines(result['subroutines'])
         if self.verbose:
-            print "SUBS", self.subroutines_names
+            print("SUBS", self.subroutines_names)
 
         # Parse all 'mains' 
         self.parse_mains(result['mains'])
         if self.verbose:
-            print "MAINS", self.mains_names
+            print("MAINS", self.mains_names)
         # we will use 0x340000 (now in pointers) to point to the right one
         # initializing Main pointer if does not exist in file
         if 'Main' not in self.pointers:
diff --git a/python/scope.py b/python/scope.py
index f1190b2..476bb8e 100644
--- a/python/scope.py
+++ b/python/scope.py
@@ -1,4 +1,4 @@
-# !/usr/bin/env python
+#! /usr/bin/env python
 
 # Displays a scope-like view of the CCD output waveform, matched with the clock states of the sequencer function
 # used to acquire it.
@@ -19,10 +19,10 @@
 # Syntax in a script:
 # import scope
 # scope.combined_scope_display("dsi-scan.fits", "tm-scan.fits", "seq-newflush.txt", "Segment00", "Bias")
-
+from __future__ import print_function
 import sys
 import os.path
-import pyfits
+from astropy.io import fits
 import numpy as np
 from matplotlib import pyplot as plt
 
@@ -69,7 +69,7 @@ def get_scandata_fromfile(inputfile, datadir='', selectchannels=None):
         else:
             displayamps = selectchannels
 
-        hdulist = pyfits.open(os.path.join(datadir, inputfile))
+        hdulist = fits.open(os.path.join(datadir, inputfile))
         imgdata = []
         for i in displayamps:
             imgdata.append(hdulist[i + 1].data)
@@ -356,7 +356,7 @@ def compare_scope_display(scanlist, labellist, datadir='', displayamps=range(16)
         except:
             ndisplay -= 1
             continue
-    print "Found %d scan files to display" % ndisplay
+    print("Found %d scan files to display" % ndisplay)
     # plot
     fig, axes = plt.subplots(nrows=(len(displayamps) + 3)/4, ncols=4, figsize=(14, 9))
 
@@ -460,36 +460,43 @@ def cut_scan_plot(scanfile, cutcolumns=[180], datadir='', polynomfit=True, displ
         figY.savefig(os.path.join(datadir, 'plotscanfit' + rootname + '.png'))
 
     else:
-        fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 12))
-
-        # plots mean and standard deviation along line direction
-        axes[0].set_title(rootname)
-        axes[0].set_xlim(0, Nbins)
-        axes[0].set_ylim(0, img[:, 5:-2].max())
-        for c in range(Nchan):
-            axes[0].plot(img[c].mean(axis=0), color=color_idx[c])
-        axes[0].set_xlabel('Time increment (10 ns)')
-        axes[0].set_ylabel('Average (ADU)')
-        axes[0].grid(True)
-
-        axes[1].set_xlim(0, Nbins)
-        maxstd = 0
-        for c in range(Nchan):
-            stdscan = img[c].std(axis=0)
-            maxstd = max(maxstd, stdscan[5:-1].max())
-            axes[1].plot(stdscan, color=color_idx[c], label="Ch%02d" % c)
-        axes[1].set_ylim(0, maxstd)
-        axes[1].set_xlabel('Time increment (10 ns)')
-        axes[1].set_ylabel('Dispersion (ADU)')
-        axes[1].grid(True)
-        # single legend and title
-        axes[1].legend(bbox_to_anchor=(1.02, 0), loc='lower left', borderaxespad=0.)
+        plot_scan_dispersion(img, len(displayamps), title=rootname)
         plt.suptitle('Statistics for %s' % rootname, fontsize='large')
         plt.savefig(os.path.join(datadir, 'scanstats' + rootname + '.png'))
 
     plt.show()
 
 
+def plot_scan_dispersion(img, Nchan=16, title=''):
+    Nlines = img.shape[1]
+    Nbins = img.shape[2]
+    color_idx = [plt.cm.jet(i) for i in np.linspace(0, 1, Nchan)]
+    fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 12))
+
+    # plots mean and standard deviation along line direction
+    axes[0].set_title(title)
+    axes[0].set_xlim(0, Nbins)
+    axes[0].set_ylim(0, img[:, 5:-2].max())
+    for c in range(Nchan):
+        axes[0].plot(img[c].mean(axis=0), color=color_idx[c])
+    axes[0].set_xlabel('Time increment (10 ns)')
+    axes[0].set_ylabel('Average (ADU)')
+    axes[0].grid(True)
+
+    axes[1].set_xlim(0, Nbins)
+    maxstd = 0
+    for c in range(Nchan):
+        stdscan = img[c].std(axis=0)
+        maxstd = max(maxstd, stdscan[5:-1].max())
+        axes[1].plot(stdscan, color=color_idx[c], label="Ch%02d" % c)
+    axes[1].set_ylim(0, maxstd)
+    axes[1].set_xlabel('Time increment (10 ns)')
+    axes[1].set_ylabel('Dispersion (ADU)')
+    axes[1].grid(True)
+    # single legend and title
+    axes[1].legend(bbox_to_anchor=(1.02, 0), loc='lower left', borderaxespad=0.)
+
+
 def stats_scan_plot(scanfile, datadir='', basecols=slice(70, 90), signalcols=slice(140, 160)):
     """
     Computes average and standard deviation of scan plots in two bands.
@@ -579,7 +586,7 @@ def stitch_long_scan(scanfile, niter, datadir='', displayamps=range(16)):
     """
     nchan = len(displayamps)
 
-    itm = pyfits.open(os.path.join(datadir,scanfile))
+    itm = fits.open(os.path.join(datadir,scanfile))
     tmscope = np.zeros((nchan, 256 * niter))
 
     for chan in range(nchan):
@@ -645,7 +652,7 @@ def stats_long_scan(scanfile, niter, datadir='', displayamps=range(16)):
 
     nchan = len(displayamps)
 
-    itm = pyfits.open(os.path.join(datadir, scanfile))
+    itm = fits.open(os.path.join(datadir, scanfile))
     nsplitlines = itm[1].header["NAXIS2"] / niter
     tmscope = np.zeros((nchan, nsplitlines, 256 * niter))
 
diff --git a/python/sequencer.py b/python/sequencer.py
index 58619d5..2ca3904 100644
--- a/python/sequencer.py
+++ b/python/sequencer.py
@@ -4,7 +4,7 @@
 # Python objects for the FPGA sequencer
 #
 # Authors: Laurent Le Guillou, Claire Juramy
-
+from __future__ import print_function
 import re
 import bidi
 import numpy as np
@@ -412,7 +412,7 @@ def fromstring(cls, s):
         # JSR addr
         m = cls.pattern_JSR_addr.match(s)
         if m is not None:
-            print m.groups()
+            print(m.groups())
             address = int(m.group(1), base=16)
             repeat = int(m.group(3))
             return Instruction(opcode="JSR",
@@ -421,7 +421,7 @@ def fromstring(cls, s):
 
         # JSR name
         m = cls.pattern_JSR_name.match(s)
-        print m, s
+        print(m, s)
         if m is not None:
             subroutine = m.group(1)
             repeat = int(m.group(2))
@@ -569,12 +569,12 @@ def fromstring(cls, s):
 
         prg = Program_UnAssembled()
 
-        print lines
+        print(lines)
 
         for iline in xrange(nlines):
-            print iline + 1
+            print(iline + 1)
             line = lines[iline]
-            print line
+            print(line)
             elts = line.split()
 
             if len(elts) < 1:
@@ -597,7 +597,7 @@ def fromstring(cls, s):
             s = " ".join(elts)
 
             instr = Instruction.fromstring(s)
-            print "INSTR = ", instr
+            print("INSTR = ", instr)
             if instr is None:
                 continue