diff --git a/doc/source/images/tutorial1a.py b/doc/source/images/tutorial1a.py
index ce906eadcb3f961bd8ca51bc339b616ec5affd33..0e14b6ebb96b3946a02aa802189863651ada42c1 100644
--- a/doc/source/images/tutorial1a.py
+++ b/doc/source/images/tutorial1a.py
@@ -76,14 +76,14 @@ for j in xrange(W):
 sys.attach_lead(lead0)
 sys.attach_lead(lead1)
 
-# finalize the system
+# Plot it, to make sure it's OK
 
-fsys = sys.finalized()
+kwant.plot(sys, "tutorial1a_sys.pdf", width=latex.figwidth_pt)
+kwant.plot(sys, "tutorial1a_sys.png", width=html.figwidth_px)
 
-# and plot it, to make sure it's proper
+# finalize the system
 
-kwant.plot(fsys, "tutorial1a_sys.pdf", width=latex.figwidth_pt)
-kwant.plot(fsys, "tutorial1a_sys.png", width=html.figwidth_px)
+fsys = sys.finalized()
 
 # Now that we have the system, we can compute conductance
 
diff --git a/doc/source/images/tutorial2a.py b/doc/source/images/tutorial2a.py
index 942b6bdf0a2e31e9e3e05a8f1808ada7e77c34ba..d0a1238fbae42caa3b7cb05bb1dd5a7a01e220ca 100644
--- a/doc/source/images/tutorial2a.py
+++ b/doc/source/images/tutorial2a.py
@@ -67,7 +67,7 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
 
-    return sys.finalized()
+    return sys
 
 def plot_conductance(fsys, energies):
     # Compute conductance
@@ -94,7 +94,10 @@ def plot_conductance(fsys, energies):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     # We should see non-monotonic conductance steps.
     plot_conductance(fsys, energies=[0.01 * i - 0.3 for i in xrange(100)])
diff --git a/doc/source/images/tutorial2b.py b/doc/source/images/tutorial2b.py
index de4a8c68af2ce04ad4bddb1e88dde07f86649b6d..99861e3b7e2647ab22101774fe4d1cd8edb45f19 100644
--- a/doc/source/images/tutorial2b.py
+++ b/doc/source/images/tutorial2b.py
@@ -62,7 +62,7 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
 
-    return sys.finalized()
+    return sys
 
 def plot_conductance(fsys, energy, welldepths):
     # We specify that we want to not only read, but also write to a
@@ -96,7 +96,10 @@ def plot_conductance(fsys, energy, welldepths):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     # We should see conductance steps.
     plot_conductance(fsys, energy=0.2,
diff --git a/doc/source/images/tutorial2c.py b/doc/source/images/tutorial2c.py
index 2188d3d4fa516ed885a75b36dea2d3b61347656c..e876d9f343c818d86dc404bad76481a79497ef1b 100644
--- a/doc/source/images/tutorial2c.py
+++ b/doc/source/images/tutorial2c.py
@@ -80,7 +80,7 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
 
-    return sys.finalized()
+    return sys
 
 
 def make_system_note1(a=1, t=1.0, W=10, r1=10, r2=20):
@@ -104,7 +104,7 @@ def make_system_note1(a=1, t=1.0, W=10, r1=10, r2=20):
     lead1 = lead0.reversed()
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
-    return sys.finalized()
+    return sys
 
 
 def make_system_note2(a=1, t=1.0, W=10, r1=10, r2=20):
@@ -128,7 +128,7 @@ def make_system_note2(a=1, t=1.0, W=10, r1=10, r2=20):
     lead1 = lead0.reversed()
     sys.attach_lead(lead0)
     sys.attach_lead(lead1, lat(0, 0))
-    return sys.finalized()
+    return sys
 
 
 def plot_conductance(fsys, energy, fluxes):
@@ -162,23 +162,26 @@ def plot_conductance(fsys, energy, fluxes):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
 
     # Check that the system looks as intended.
-    kwant.plot(fsys, "tutorial2c_sys.pdf", width=latex.figwidth_pt)
-    kwant.plot(fsys, "tutorial2c_sys.png", width=html.figwidth_px)
+    kwant.plot(sys, "tutorial2c_sys.pdf", width=latex.figwidth_pt)
+    kwant.plot(sys, "tutorial2c_sys.png", width=html.figwidth_px)
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     # We should see a conductance that is periodic with the flux quantum
     plot_conductance(fsys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
                                                 for i in xrange(100)])
 
     # Finally, some plots needed for the notes
-    fsys = make_system_note1()
-    kwant.plot(fsys, "tutorial2c_note1.pdf", width=latex.figwidth_small_pt)
-    kwant.plot(fsys, "tutorial2c_note1.png", width=html.figwidth_small_px)
-    fsys = make_system_note2()
-    kwant.plot(fsys, "tutorial2c_note2.pdf", width=latex.figwidth_small_pt)
-    kwant.plot(fsys, "tutorial2c_note2.png", width=html.figwidth_small_px)
+    sys = make_system_note1()
+    kwant.plot(sys, "tutorial2c_note1.pdf", width=latex.figwidth_small_pt)
+    kwant.plot(sys, "tutorial2c_note1.png", width=html.figwidth_small_px)
+    sys = make_system_note2()
+    kwant.plot(sys, "tutorial2c_note2.pdf", width=latex.figwidth_small_pt)
+    kwant.plot(sys, "tutorial2c_note2.png", width=html.figwidth_small_px)
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/images/tutorial3a.py b/doc/source/images/tutorial3a.py
index 2fa54f7a480742f11a541b7ac66b93164e303ed8..843886e405319659e6b25f0a1fdb03f1f9376a11 100644
--- a/doc/source/images/tutorial3a.py
+++ b/doc/source/images/tutorial3a.py
@@ -36,7 +36,7 @@ def make_lead(a=1, t=1.0, W=10):
         lead[(1, j), (0, j)] = - t
 
     # return a finalized lead
-    return lead.finalized()
+    return lead
 
 
 def plot_bandstructure(flead, momenta):
@@ -62,7 +62,7 @@ def plot_bandstructure(flead, momenta):
 
 
 def main():
-    flead = make_lead()
+    flead = make_lead().finalized()
 
     # list of momenta at which the bands should be computed
     momenta = np.arange(-pi, pi + .01, 0.02 * pi)
diff --git a/doc/source/images/tutorial3b.py b/doc/source/images/tutorial3b.py
index c0ac9db3e3d0bbfd2654228d776ae31f8f7faaa6..bac129f2decfae627388395ee0f74fc0a3908805 100644
--- a/doc/source/images/tutorial3b.py
+++ b/doc/source/images/tutorial3b.py
@@ -47,7 +47,7 @@ def make_system(a=1, t=1.0, r=10):
     sys[sys.possible_hoppings((0, 1), lat, lat)] = - t
 
     # It's a closed system for a change, so no leads
-    return sys.finalized()
+    return sys
 
 
 def plot_spectrum(fsys, Bfields):
@@ -89,7 +89,10 @@ def plot_spectrum(fsys, Bfields):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     # We should observe energy levels that flow towards Landau
     # level energies with increasing magnetic field
diff --git a/doc/source/images/tutorial4.py b/doc/source/images/tutorial4.py
index 59795c5a582ee0b3fde68928548c9b20107bc3fd..e93014920ddbcac207e97e07e38c5a6dbadaf3ad 100644
--- a/doc/source/images/tutorial4.py
+++ b/doc/source/images/tutorial4.py
@@ -56,10 +56,6 @@ def make_system(r=10, w=2.0, pot=0.1):
     del sys[a(0,0)]
     sys[a(-2,1), b(2, 2)] = -1
 
-    # Keep a copy of the closed system without leads, for
-    # eigenvalue computations
-    closed_fsys = sys.finalized()
-
     #### Define the leads. ####
     # left lead
     sym0 = kwant.TranslationalSymmetry([graphene.vec((-1, 0))])
@@ -73,7 +69,7 @@ def make_system(r=10, w=2.0, pot=0.1):
     for hopping in hoppings:
         lead0[lead0.possible_hoppings(*hopping)] = - 1
 
-    # The second lead, going ot the top right
+    # The second lead, going to the top right
     sym1 = kwant.TranslationalSymmetry([graphene.vec((0, 1))])
 
     def lead1_shape(pos):
@@ -87,11 +83,7 @@ def make_system(r=10, w=2.0, pot=0.1):
     for hopping in hoppings:
         lead1[lead1.possible_hoppings(*hopping)] = - 1
 
-    # Attach the leads
-    sys.attach_lead(lead0)
-    sys.attach_lead(lead1)
-
-    return sys.finalized(), closed_fsys, lead0.finalized()
+    return sys, [lead0, lead1]
 
 
 def compute_evs(sys):
@@ -158,9 +150,7 @@ def plot_bandstructure(flead, momenta):
 
 def main():
     pot = 0.1
-    fsys, closed_fsys, flead = make_system(pot=pot)
-
-    # First, plot the closed system, and compute some eigenvalues
+    sys, leads = make_system(pot=pot)
 
     # To highlight the two sublattices of graphene, we plot one with
     # a filled, and the other one with an open circle:
@@ -169,26 +159,33 @@ def main():
                                                fcol=kwant.plotter.white,
                                                lcol=kwant.plotter.black)}
 
-    kwant.plot(closed_fsys, a=1./sqrt(3.), symbols=plotter_symbols,
-               filename="tutorial4_sys1.pdf",
-               width=latex.figwidth_pt)
-    kwant.plot(closed_fsys, a=1./sqrt(3.), symbols=plotter_symbols,
-               filename="tutorial4_sys1.png",
-               width=html.figwidth_px)
+    # Plot the closed system without leads.
+    kwant.plot(sys, symbols=plotter_symbols,
+               filename="tutorial4_sys1.pdf", width=latex.figwidth_pt)
+    kwant.plot(sys, symbols=plotter_symbols,
+               filename="tutorial4_sys1.png", width=html.figwidth_px)
+
+    # Compute some eigenvalues.
+    compute_evs(sys.finalized())
+
+    # Attach the leads to the system.
+    for lead in leads:
+        sys.attach_lead(lead)
 
-    # Then, plot the system with leads and compute the band structure
-    # of one of the (zigzag) leads, as well as the conductance
+    # Then, plot the system with leads.
+    kwant.plot(sys, symbols=plotter_symbols,
+               filename="tutorial4_sys2.pdf", width=latex.figwidth_pt)
+    kwant.plot(sys, symbols=plotter_symbols,
+               filename="tutorial4_sys2.png", width=html.figwidth_px)
 
-    kwant.plot(fsys, a=1/sqrt(3.), symbols=plotter_symbols,
-               filename="tutorial4_sys2.png",
-               width=html.figwidth_px)
-    kwant.plot(fsys, a=1/sqrt(3.), symbols=plotter_symbols,
-               filename="tutorial4_sys2.pdf",
-               width=latex.figwidth_pt)
+    # Finalize the system.
+    fsys = sys.finalized()
 
+    # Compute the band structure of lead 0.
     momenta = np.arange(-pi, pi + .01, 0.1 * pi)
-    plot_bandstructure(flead, momenta)
+    plot_bandstructure(fsys.leads[0], momenta)
 
+    # Plot conductance.
     energies = np.arange(-2 * pot, 2 * pot, pot / 10.5)
     plot_conductance(fsys, energies)
 
diff --git a/doc/source/tutorial/tutorial1.rst b/doc/source/tutorial/tutorial1.rst
index c35ba6df60f612545ef390c55d1ced26b35be21b..59ff4096511d4b96af717b109db0070a0dbdeb40 100644
--- a/doc/source/tutorial/tutorial1.rst
+++ b/doc/source/tutorial/tutorial1.rst
@@ -102,18 +102,13 @@ attached at the right position using:
 More details about attaching leads can be found in the tutorial
 :ref:`tutorial-abring`.
 
-Now we have finished building our system! We need to finalize it, in
-order to use it for a transport calculation:
+Now we have finished building our system! We plot it, to make sure we didn't
+make any mistakes:
 
 .. literalinclude:: ../../../examples/tutorial1a.py
     :lines: 79
 
-and should plot it, to make sure we didn't make mistakes:
-
-.. literalinclude:: ../../../examples/tutorial1a.py
-    :lines: 83
-
-This command should bring up this picture:
+This should bring up this picture:
 
 .. image:: /images/tutorial1a_sys.*
 
@@ -123,6 +118,11 @@ nonzero hopping element between points there is a line connecting these
 points. From the leads, only a few (default 2) unit cells are shown, with
 fading color.
 
+In order to use our system for a transport calculation, we need to finalize it
+
+.. literalinclude:: ../../../examples/tutorial1a.py
+    :lines: 83
+
 Having successfully created a system, we now can immediately start to compute
 its conductance as a function of energy:
 
@@ -360,14 +360,14 @@ We collect all statements that should be executed in the script
 in a ``main()``-function:
 
 .. literalinclude:: ../../../examples/tutorial1b.py
-    :lines: 66-73
+    :lines: 66-76
 
 Finally, we use the following python construct [#]_ that executes
 ``main()`` if the program is used as a script (i.e. executed as
 ``python tutorial1b.py``):
 
 .. literalinclude:: ../../../examples/tutorial1b.py
-    :lines: 78-79
+    :lines: 81-82
 
 If the example however is imported using ``import tutorial1b``,
 ``main()`` is not executed automatically. Instead, you can execute it
diff --git a/doc/source/tutorial/tutorial4.rst b/doc/source/tutorial/tutorial4.rst
index 8f9b92858ad364d6dfcf54f0256c5147316d9a56..85ae898e7ee5b2c19048197ca5ec09882cffe7f3 100644
--- a/doc/source/tutorial/tutorial4.rst
+++ b/doc/source/tutorial/tutorial4.rst
@@ -17,11 +17,11 @@ explicitly here to show how to define a new lattice:
 .. literalinclude:: ../../../examples/tutorial4.py
     :lines: 24-27
 
-The first argument to the `make_lattice` function is the list of primitive
-vectors of the lattice; the second one is the coordinates of basis atoms.
-The honeycomb lattice has two basis atoms. Each type of basis atom by itself
-forms a regular lattice of the same type as well, and those *sublattices*
-are referenced as `a` and `b` above.
+The first argument to the `~kwant.lattice.make_lattice` function is the list of
+primitive vectors of the lattice; the second one is the coordinates of basis
+atoms.  The honeycomb lattice has two basis atoms. Each type of basis atom by
+itself forms a regular lattice of the same type as well, and those
+*sublattices* are referenced as `a` and `b` above.
 
 In the next step we define the shape of the scattering region (circle again)
 and add all lattice points using the ``shape()``-functionality:
@@ -71,32 +71,28 @@ as add an additional link:
 Note that the conversion from a tuple `(i,j)` to site
 is done by the sublattices `a` and `b`.
 
-Later, we will compute some eigenvalues of the closed
-scattering region without leads. For that, obtain a finalized
-snapshot of the system:
+The leads are defined as before:
 
 .. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 60
+    :lines: 60-83
 
-Adding leads to the scattering region is done as before:
+Note that the translational vectors ``graphene.vec((-1, 0))`` and
+``graphene.vec((0, 1))`` are *not* orthogonal any more as they would have been
+in a square lattice -- they follow the non-orthogonal primitive vectors defined
+in the beginning.
 
-.. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 64-93
-
-Note here that the translational vectors ``graphene.vec((-1, 0))`` and
-``graphene.vec((0, 1))`` are *not* orthogonal any more as they would
-have been in a square lattice-- they follow the non-orthogonal
-primitive vectors defined in the beginning.
+Later, we will compute some eigenvalues of the closed scattering region without
+leads. This is why we postpone attaching the leads to the system. Instead,
+we return the scattering region and the leads separately.
 
-In the end, we return not only the finalized system with leads, but
-also a finalized copy of the closed system (for eigenvalues)
-as well as a finalized lead (for band structure calculations).
+.. literalinclude:: ../../../examples/tutorial4.py
+    :lines: 85
 
 The computation of some eigenvalues of the closed system is done
 in the following piece of code:
 
 .. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 96-101, 104-107
+    :lines: 88-93, 96-99
 
 Here we use in contrast to the previous example a sparse matrix and
 the sparse linear algebra functionality of scipy (this requires
@@ -111,45 +107,47 @@ Finally, in the `main()` function we make and
 plot the system:
 
 .. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 135-137, 142-147
-
-Here we customize the plotting: `plotter_symbols` is a dictionary
-with the sublattice objects `a` and `b` as keys, and the
-`~kwant.plotter.Circle` objects specify that the sublattice `a` should
-be drawn using a filled black circle, and `b` using a white circle
-with a black outline. The radius of the circle is given in relative
-units: :func:`plot <kwant.plotter.plot>` uses a typical length
-scale as a reference length. By default, the typical length scale is
-the smallest distance between lattice points.  :func:`plot
-<kwant.plotter.plot>` can find this length by itself, but must then go
-through all hoppings. Alternatively, one can specify the typical
-length scale using the argument `a` as in the example (not to be
-confused with the sublattice `a`) which is here set to the distance
-between carbon atoms in the graphene lattice. Specifying ``r=0.3`` in
-`~kwant.plotter.Circle` hence means that the radius of the circle is
-30% of the carbon-carbon distance. Using this relative unit it is
-easy to make good-looking plots where the symbols cover a well-defined
-part of the plot.
+    :lines: 126-129, 138
+
+We customize the plotting: `plotter_symbols` is a dictionary with the
+sublattice objects `a` and `b` as keys, and the `~kwant.plotter.Circle` objects
+specify that the sublattice `a` should be drawn using a filled black circle,
+and `b` using a white circle with a black outline.
+
+.. literalinclude:: ../../../examples/tutorial4.py
+    :lines: 132-135
+
+The radius of the circle is given in relative units: `~kwant.plotter.plot` uses
+a typical length scale as a reference length. By default, the typical length
+scale is the smallest distance between lattice points.  `~kwant.plotter.plot`
+can find this length by itself, but must then go through all
+hoppings. Alternatively, one can specify the typical length scale using the
+argument `a` as in the example (not to be confused with the sublattice `a`)
+which is here set to the distance between carbon atoms in the graphene
+lattice. Specifying ``r=0.3`` in `~kwant.plotter.Circle` hence means that the
+radius of the circle is 30% of the carbon-carbon distance. Using this relative
+unit it is easy to make good-looking plots where the symbols cover a
+well-defined part of the plot.
 
 Plotting the closed system gives this result:
 
 .. image:: ../images/tutorial4_sys1.*
 
-and computing the eigenvalues of largest magnitude,
+Computing the eigenvalues of largest magnitude,
 
 .. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 148
+    :lines: 141
 
-should yield two eigenvalues similar to `[ 3.07869311
-+1.02714523e-17j, -3.06233144 -6.68085759e-18j]` (round-off might
-change the imaginary part which should in exact arithmetics be equal
-to zero).
+should yield two eigenvalues similar to `[ 3.07869311 +1.02714523e-17j,
+-3.06233144 -6.68085759e-18j]` (round-off might change the imaginary part which
+would be equal to zero for exact arithmetics).
 
-The remaining code of `main()` plots the system with leads:
+The remaining code of `main()` attaches the leads to the system and plots it
+again:
 
 .. image:: ../images/tutorial4_sys2.*
 
-It computes the band structure of one of the leads:
+It computes the band structure of one of lead 0:
 
 .. image:: ../images/tutorial4_bs.*
 
@@ -174,7 +172,7 @@ an open quantum dot
    `~kwant.plotter.Polygon`'s as predefined symbols. It is also
    easy to define any custom symbol. Furthermore, plotting offers
    many more options to customize plots. See the documentation of
-   :func:`plot <kwant.plotter.plot>` for more details.
+   `~kwant.plotter.plot` for more details.
 
  - In a lattice with more than one basis atom, you can always act either
    on all sublattices at the same time, or on a single sublattice only.
diff --git a/doc/source/whatsnew/0.2.rst b/doc/source/whatsnew/0.2.rst
index 81f43b600223d40c8dc6322ff2dd96ea0d6fe40b..6a509357e878e096049d0571a71a2b1552e8fb70 100644
--- a/doc/source/whatsnew/0.2.rst
+++ b/doc/source/whatsnew/0.2.rst
@@ -3,6 +3,23 @@ What's New in kwant 0.2
 
 This article explains the user-visible changes in kwant 0.2.
 
+`~kwant.plotter.plot` more useful for low level systems
+-------------------------------------------------------
+The behavior of `~kwant.plotter.plot` has been changed when a `low level system
+<kwant.system.System>` is plotted.  Previously, only low level systems which
+were finalized `builders <kwant.builder.Builder>` were supported and there was
+no difference between plotting a low-level system and a builder.
+
+* Arguments of plot which are functions are given site numbers in place of
+  `~kwant.builder.Site` objects when plotting a low level system.  This
+  provides an easy way to make the appearance of lines and symbols depend on
+  computation results.
+
+* Only the scattering region (without leads) is plotted for low level systems.
+
+* For plotting low level systems, dictionaries with site group keys are no
+  longer supported as plot arguments.
+
 Calculation of the local density of states
 ------------------------------------------
 The new function `kwant.solvers.sparse.ldos` allows the calculation of the
diff --git a/examples/tutorial1a.py b/examples/tutorial1a.py
index 64b3b6f7ec66f884767358ae11d50d5ae2c92523..0919347329eddc85eb66281169e65450b7660c6a 100644
--- a/examples/tutorial1a.py
+++ b/examples/tutorial1a.py
@@ -74,13 +74,13 @@ for j in xrange(W):
 sys.attach_lead(lead0)
 sys.attach_lead(lead1)
 
-# finalize the system
+# Plot it, to make sure it's OK
 
-fsys = sys.finalized()
+kwant.plot(sys)
 
-# and plot it, to make sure it's proper
+# Finalize the system
 
-kwant.plot(fsys)
+fsys = sys.finalized()
 
 # Now that we have the system, we can compute conductance
 
diff --git a/examples/tutorial1b.py b/examples/tutorial1b.py
index cee11853234d23d8267a4b306cbe3f19cc9d55a0..0b8ee7fe48a6d316684ca25e8e7c5de39b8373cf 100644
--- a/examples/tutorial1b.py
+++ b/examples/tutorial1b.py
@@ -44,11 +44,11 @@ def make_system(a=1, t=1.0, W=10, L=30):
     # `lead0` with its direction reversed.
     lead1 = lead0.reversed()
 
-    #### Attach the leads and return the finalized system. ####
+    #### Attach the leads and return the system. ####
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
 
-    return sys.finalized()
+    return sys
 
 def plot_conductance(fsys, energies):
     # Compute conductance
@@ -64,10 +64,13 @@ def plot_conductance(fsys, energies):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
 
     # Check that the system looks as intended.
-    kwant.plot(fsys)
+    kwant.plot(sys)
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     # We should see conductance steps.
     plot_conductance(fsys, energies=[0.01 * i for i in xrange(100)])
diff --git a/examples/tutorial2a.py b/examples/tutorial2a.py
index 30e4ee53ee4275e643b8334e3b4ae2f8f96dac89..d88a1746374233753d31f794055303446f74707f 100644
--- a/examples/tutorial2a.py
+++ b/examples/tutorial2a.py
@@ -65,7 +65,7 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
 
-    return sys.finalized()
+    return sys
 
 def plot_conductance(fsys, energies):
     # Compute conductance
@@ -81,10 +81,13 @@ def plot_conductance(fsys, energies):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
 
     # Check that the system looks as intended.
-    kwant.plot(fsys)
+    kwant.plot(sys)
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     # We should see non-monotonic conductance steps.
     plot_conductance(fsys, energies=[0.01 * i - 0.3 for i in xrange(100)])
diff --git a/examples/tutorial2b.py b/examples/tutorial2b.py
index 4bd284c5e3dac3eaee79dc01a7c1a6e33df8f618..94b778b8117bf29c73a3761f9aa95b8e74f572e2 100644
--- a/examples/tutorial2b.py
+++ b/examples/tutorial2b.py
@@ -60,7 +60,7 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
 
-    return sys.finalized()
+    return sys
 
 def plot_conductance(fsys, energy, welldepths):
     # We specify that we want to not only read, but also write to a
@@ -83,10 +83,13 @@ def plot_conductance(fsys, energy, welldepths):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
 
     # Check that the system looks as intended.
-    kwant.plot(fsys)
+    kwant.plot(sys)
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     # We should see conductance steps.
     plot_conductance(fsys, energy=0.2,
diff --git a/examples/tutorial2c.py b/examples/tutorial2c.py
index 855bbc29d835abaff2fd1999164a1694d2f3b7e8..9a4cdbfde427a15245c6d3ddfa0f6699afa490af 100644
--- a/examples/tutorial2c.py
+++ b/examples/tutorial2c.py
@@ -78,7 +78,7 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
 
-    return sys.finalized()
+    return sys
 
 
 def plot_conductance(fsys, energy, fluxes):
@@ -101,10 +101,13 @@ def plot_conductance(fsys, energy, fluxes):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
 
     # Check that the system looks as intended.
-    kwant.plot(fsys)
+    kwant.plot(sys)
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     # We should see a conductance that is periodic with the flux quantum
     plot_conductance(fsys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
diff --git a/examples/tutorial3a.py b/examples/tutorial3a.py
index 63cf1b1363f5041a2e752e69f62f7768da97852c..058b4e20369710009f3a746e7c29f3a3610684a3 100644
--- a/examples/tutorial3a.py
+++ b/examples/tutorial3a.py
@@ -33,8 +33,7 @@ def make_lead(a=1, t=1.0, W=10):
 
         lead[(1, j), (0, j)] = - t
 
-    # return a finalized lead
-    return lead.finalized()
+    return lead
 
 
 def plot_bandstructure(flead, momenta):
@@ -49,7 +48,7 @@ def plot_bandstructure(flead, momenta):
 
 
 def main():
-    flead = make_lead()
+    flead = make_lead().finalized()
 
     # list of momenta at which the bands should be computed
     momenta = np.arange(-pi, pi + .01, 0.02 * pi)
diff --git a/examples/tutorial3b.py b/examples/tutorial3b.py
index 08f317fe0a262565f4724a7986bde076bd2057e1..273c8d4b423d071000f58300e9fea41b67e46768 100644
--- a/examples/tutorial3b.py
+++ b/examples/tutorial3b.py
@@ -44,7 +44,7 @@ def make_system(a=1, t=1.0, r=10):
     sys[sys.possible_hoppings((0, 1), lat, lat)] = - t
 
     # It's a closed system for a change, so no leads
-    return sys.finalized()
+    return sys
 
 
 def plot_spectrum(fsys, Bfields):
@@ -75,10 +75,13 @@ def plot_spectrum(fsys, Bfields):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
 
     # Check that the system looks as intended.
-    kwant.plot(fsys)
+    kwant.plot(sys)
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     # We should observe energy levels that flow towards Landau
     # level energies with increasing magnetic field
diff --git a/examples/tutorial4.py b/examples/tutorial4.py
index 74ea9e6f74270ab09ac47e5c06bac602cf2e1764..a31f00dde0a5965b07619fc5a4ee4da6818ec391 100644
--- a/examples/tutorial4.py
+++ b/examples/tutorial4.py
@@ -55,10 +55,6 @@ def make_system(r=10, w=2.0, pot=0.1):
     del sys[a(0,0)]
     sys[a(-2,1), b(2, 2)] = -1
 
-    # Keep a copy of the closed system without leads, for
-    # eigenvalue computations
-    closed_fsys = sys.finalized()
-
     #### Define the leads. ####
     # left lead
     sym0 = kwant.TranslationalSymmetry([graphene.vec((-1, 0))])
@@ -72,7 +68,7 @@ def make_system(r=10, w=2.0, pot=0.1):
     for hopping in hoppings:
         lead0[lead0.possible_hoppings(*hopping)] = - 1
 
-    # The second lead, going ot the top right
+    # The second lead, going to the top right
     sym1 = kwant.TranslationalSymmetry([graphene.vec((0, 1))])
 
     def lead1_shape(pos):
@@ -86,11 +82,7 @@ def make_system(r=10, w=2.0, pot=0.1):
     for hopping in hoppings:
         lead1[lead1.possible_hoppings(*hopping)] = - 1
 
-    # Attach the leads
-    sys.attach_lead(lead0)
-    sys.attach_lead(lead1)
-
-    return sys.finalized(), closed_fsys, lead0.finalized()
+    return sys, [lead0, lead1]
 
 
 def compute_evs(sys):
@@ -133,9 +125,7 @@ def plot_bandstructure(flead, momenta):
 
 def main():
     pot = 0.1
-    fsys, closed_fsys, flead = make_system(pot=pot)
-
-    # First, plot the closed system, and compute some eigenvalues
+    sys, leads = make_system(pot=pot)
 
     # To highlight the two sublattices of graphene, we plot one with
     # a filled, and the other one with an open circle:
@@ -144,17 +134,27 @@ def main():
                                                fcol=kwant.plotter.white,
                                                lcol=kwant.plotter.black)}
 
-    kwant.plot(closed_fsys, a=1./sqrt(3.), symbols=plotter_symbols)
-    compute_evs(closed_fsys)
+    # Plot the closed system without leads.
+    kwant.plot(sys, symbols=plotter_symbols)
+
+    # Compute some eigenvalues.
+    compute_evs(sys.finalized())
+
+    # Attach the leads to the system.
+    for lead in leads:
+        sys.attach_lead(lead)
 
-    # Then, plot the system with leads and compute the band structure
-    # of one of the (zigzag) leads, as well as the conductance
+    # Then, plot the system with leads.
+    kwant.plot(sys, symbols=plotter_symbols)
 
-    kwant.plot(fsys, a=1/sqrt(3.), symbols=plotter_symbols)
+    # Finalize the system.
+    fsys = sys.finalized()
 
+    # Compute the band structure of lead 0.
     momenta = np.arange(-pi, pi + .01, 0.1 * pi)
-    plot_bandstructure(flead, momenta)
+    plot_bandstructure(fsys.leads[0], momenta)
 
+    # Plot conductance.
     energies = np.arange(-2 * pot, 2 * pot, pot / 10.5)
     plot_conductance(fsys, energies)
 
diff --git a/kwant/plotter.py b/kwant/plotter.py
index b9bbb6faa89a2c43899bcbd306c16d98e7876989..a69c7e45b848689a3e86db07eeb2db80938e416c 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -275,36 +275,44 @@ def plot(system, filename=defaultname, fmt=None, a=None,
     """Plot two-dimensional systems (or two-dimensional representations
     of a system).
 
-    `plot` can be used to plot both unfinalized systems derived from
-    kwant.builder.Builder, or the corresponding finalized systems.
-
-    The output of `plot` is highly modifyable, as it does not perform
-    any drawing itself, but instead lets objects passed by the user
-    (or as default parameters) do the actual drawing work. `plot`
-    itself does figure out the range of positions occupied by the
-    sites, as well as the smallest distance between two sites which
-    then serves as a reference length, unless the user specifies
-    explicitely a reference length. This reference length is then
-    used so that the sizes of symbols or lines are always given relative
-    to that reference length. This is particularly advantageous for
-    regular lattices, as it makes it easy to specify the area covered
-    by symbols, etc.
-
-    The objects that determine `plot`'s behavior are symbol_like
-    (symbols representing sites), line_like (lines representing
-    hoppings) and color_like (representing colors). The notes below
-    explain in detail how to implement custom classes. In most cases
-    it is enough to use the predefined standard objects:
+    `plot` can be used to plot both unfinalized kwant.builder.Builder
+    instances, and low level systems (i.e. instances of
+    kwant.system.FiniteSystem), including finalized builders.
+
+    This function behaves differently for builders and low-level systems:
+    builders are plotted including those of their leads which are builders
+    themselves.  For the leads, several copies of the lead unit cell are
+    plotted (per default 2), and they are gradually faded towards the
+    background color (at least in the default behavior).  For low-level systems
+    the leads are ignored as there is no general way to recover the necessary
+    information about leads for low level systems.
+
+    When arguments to this function are functions themselves, "sites" will be
+    passed to them as arguments.  The meaning of "site" depends on whether the
+    system to be plotted is a builder or a low level system.  For builders, a
+    site is a kwant.builder.Site object.  For low level systems, a site is an
+    integer -- the site number.
+
+    The output of `plot` is highly modifyable, as it does not perform any
+    drawing itself, but instead lets objects passed by the user (or as default
+    parameters) do the actual drawing work. `plot` itself does figure out the
+    range of positions occupied by the sites, as well as the smallest distance
+    between two sites which then serves as a reference length, unless the user
+    specifies explicitely a reference length. This reference length is then
+    used so that the sizes of symbols or lines are always given relative to
+    that reference length. This is particularly advantageous for regular
+    lattices, as it makes it easy to specify the area covered by symbols, etc.
+
+    The objects that determine `plot`'s behavior are symbol_like (symbols
+    representing sites), line_like (lines representing hoppings) and color_like
+    (representing colors). The notes below explain in detail how to implement
+    custom classes. In most cases it is enough to use the predefined standard
+    objects:
 
     - for symbol_like: `Circle` and `Polygon`
     - for line_like: `Line`
     - for color_like: `Color`.
 
-    `plot` draws both system sites, as well as sites corresponding to the
-    leads. For the leads, several copies of the lead unit cell are plotted
-    (per default 2), and they are gradually faded towards the background
-    color (at least in the default behavior).
-
     Parameters
     ----------
     system : (un)finalized system
@@ -346,44 +354,49 @@ def plot(system, filename=defaultname, fmt=None, a=None,
         maybe [fading to bcol e.g. still makes a white symbol, not a
         transparant symbol], but then again there is no reason for
         having a white box behind everything)
-    pos : callable or None, optional
+    pos : function or None, optional
         When passed a site should return its (2D) position as a sequence of
-        length 2. If None, the method pos() of the site is used.
-        Defaults to None.
-    symbols : {symbol_like, callable, dict, None}, optional
+        length 2. If None, the real space position of the site is used if the
+        system to be plotted is a (finalized) builder.  For other low level
+        systems it is required to specify this argument and an error will be
+        reported if it is missing. Defaults to None.
+    symbols : {symbol_like, function, dict, None}, optional
         Object responsible for drawing the symbols correspodning to sites.
         Either must be a single symbol_like object (the same symbol is drawn
-        for every site, regardless of site group), a callable that
-        returns a symbol_like object when passed a site, a dictionary
-        with site groups as keys and symbol_like as values
-        (allowing to specify different symbols for different site groups),
-        or None (in which case no symbols are drawn). Instead of
-        a symbol_like object the callable or the dict may also return None
-        corresponding to no symbol. Defaults to ``Circle(r=0.3)``.
-
-        The standard symbols available are `Circle` and
-        `Polygon`.
-    lines : {line_like, callable, dict, None}, optional
-        Object responsible for drawing the lines representing the
-        hoppings between sites. Either a single line_like object
-        (the same type of line is drawn for all hoppings), a callable
-        that returns a line_like object when passed two sites,
-        a dictionary with tuples of two site groups as keys and
-        line_like objects as values (allowing to specify different
-        line styles for different hoppings; note that if the hopping
-        (a, b) is specified, (b, a) needs not be included in the dictionary),
-        or None (in which case no hoppings are drawn). Instead of
-        a line_like object the callable or the dict may also return None
+        for every site), a function that returns a symbol_like object when
+        passed a site, or None (in which case no symbols are drawn). Instead of
+        a symbol_like object the function may also return None corresponding to
+        no symbol.
+
+        If the system is a builder, `symbols` may also be a dictionary with
+        site groups as keys and symbol_like as values.  This allows to specify
+        different symbols for different site groups.
+
+        Defaults to ``Circle(r=0.3)``.
+
+        The standard symbols available are `Circle` and `Polygon`.
+    lines : {line_like, function, dict, None}, optional
+        Object responsible for drawing the lines representing the hoppings
+        between sites. Either a single line_like object (the same type of line
+        is drawn for all hoppings), a function that returns a line_like object
+        when passed two sites, or None (in which case no hoppings are
+        drawn). Instead of a line_like object the function may also return None
         corresponding to no line. Defaults to ``Line(lw=0.1)``.
 
+        If the system is a builder, `lines` may also be a dictionary with
+        tuples of two site groups as keys and line_like objects as values.
+        This allows to specify different line styles for different hoppings.
+        Note that if the hopping (a, b) is specified, (b, a) needs not be
+        included in the dictionary.
+
         The standard line available is `Line`.
-    lead_symbols : {symbol_like, callable, dict, -1, None}, optional
+    lead_symbols : {symbol_like, function, dict, -1, None}, optional
         Symbols to be drawn for the sites in the leads. The special
         value -1 indicates that `symbols` (which is used for system sites)
         should be used also for the leads. The other possible values are
         as for the system `symbols`.
         Defaults to -1.
-    lead_lines : {line_like, callable, dict, -1, None}, optional
+    lead_lines : {line_like, function, dict, -1, None}, optional
         Lines to be drawn for the hoppings in the leads. The special
         value -1 indicates that `lines` (which is used for system hoppings)
         should be used also for the leads. The other possible values are
@@ -506,87 +519,57 @@ def plot(system, filename=defaultname, fmt=None, a=None,
                                sym.act(shift + i - 1, site2),
                                i - 1 + shift1, i - 1 + shift2)
 
-    def iterate_system_sites_builder(system):
+    def iterate_scattreg_sites_builder(system):
         for site in system.sites():
             yield site
 
-    def iterate_system_hoppings_builder(system):
+    def iterate_scattreg_hoppings_builder(system):
         for hopping in system.hoppings():
             yield hopping
 
-    def iterate_lead_sites_llsys(system, lead_copies):
-        for ilead in xrange(len(system.leads)):
-            lead = system.leads[ilead]
-            sym = lead.symmetry
-            shift = sym.which(system.site(system.lead_neighbor_seqs[ilead][0]))
-            shift += 1
-
-            for i in xrange(lead_copies):
-                for isite in xrange(lead.slice_size):
-                        yield sym.act(shift + i, lead.site(isite)), i
-
-    def iterate_lead_hoppings_llsys(system, lead_copies):
-        for ilead in xrange(len(system.leads)):
-            lead = system.leads[ilead]
-            sym = lead.symmetry
-            shift = sym.which(system.site(system.lead_neighbor_seqs[ilead][0]))
-            shift += 1
-
-            for i in xrange(lead_copies):
-                for isite in xrange(lead.slice_size):
-                    for jsite in lead.graph.out_neighbors(isite):
-                        # Note: unlike in builder, it is guaranteed
-                        #       in the finalized system that hoppings
-                        #       beyond the unit cell are in the previous
-                        #       "slice"
-                        if jsite < lead.slice_size:
-                            yield (sym.act(shift + i, lead.site(isite)),
-                                   sym.act(shift + i, lead.site(jsite)),
-                                   i, i)
-                        else:
-                            jsite -= lead.slice_size
-                            yield (sym.act(shift + i, lead.site(isite)),
-                                   sym.act(shift + i - 1, lead.site(jsite)),
-                                   i, i - 1)
-
+    def empty_generator(*args, **kwds):
+        return
+        yield
 
-    def iterate_system_sites_llsys(system):
-        for i in xrange(system.graph.num_nodes):
-            yield system.site(i)
+    def iterate_scattreg_sites_llsys(system):
+        return xrange(system.graph.num_nodes)
 
-    def iterate_system_hoppings_llsys(system):
+    def iterate_scattreg_hoppings_llsys(system):
         for i in xrange(system.graph.num_nodes):
-            site1 = system.site(i)
             for j in system.graph.out_neighbors(i):
                 # Only yield half of the hoppings (as builder does)
                 if i < j:
-                    yield (site1, system.site(j))
+                    yield i, j
 
 
     def iterate_all_sites(system, lead_copies=0):
-        for site in iterate_system_sites(system):
+        for site in iterate_scattreg_sites(system):
             yield site
 
         for site, ucindx in iterate_lead_sites(system, lead_copies):
             yield site
 
     def iterate_all_hoppings(system, lead_copies=0):
-        for site1, site2 in iterate_system_hoppings(system):
+        for site1, site2 in iterate_scattreg_hoppings(system):
             yield site1, site2
 
         for site1, site2, i1, i2 in iterate_lead_hoppings(system, lead_copies):
             yield site1, site2
 
-    if isinstance(system, kwant.builder.Builder):
-        iterate_system_sites = iterate_system_sites_builder
-        iterate_system_hoppings = iterate_system_hoppings_builder
+    is_builder = isinstance(system, kwant.builder.Builder)
+    is_lowlevel = isinstance(system, kwant.system.FiniteSystem)
+    if is_builder:
+        iterate_scattreg_sites = iterate_scattreg_sites_builder
+        iterate_scattreg_hoppings = iterate_scattreg_hoppings_builder
         iterate_lead_sites = iterate_lead_sites_builder
         iterate_lead_hoppings = iterate_lead_hoppings_builder
-    elif isinstance(system, kwant.builder.FiniteSystem):
-        iterate_system_sites = iterate_system_sites_llsys
-        iterate_system_hoppings = iterate_system_hoppings_llsys
-        iterate_lead_sites = iterate_lead_sites_llsys
-        iterate_lead_hoppings = iterate_lead_hoppings_llsys
+    elif is_lowlevel:
+        iterate_scattreg_sites = iterate_scattreg_sites_llsys
+        iterate_scattreg_hoppings = iterate_scattreg_hoppings_llsys
+        # We do not plot leads for low level systems, as there is no general
+        # way to do that.
+        iterate_lead_sites = empty_generator
+        iterate_lead_hoppings = empty_generator
     else:
         raise ValueError("Plotting not suported for given system")
 
@@ -602,7 +585,13 @@ def plot(system, filename=defaultname, fmt=None, a=None,
             raise ValueError("The distance a must be >0")
 
     if pos is None:
-        pos = lambda site: site.pos
+        if is_builder:
+            pos = lambda site: site.pos
+        elif is_lowlevel:
+            pos = lambda i: system.site(i).pos
+        else:
+            raise ValueError("`pos` argument needed when plotting"
+                             " systems which are not (finalized) builders")
 
     if fmt is None and filename is not None:
         # Try to figure out the format from the filename
@@ -618,19 +607,19 @@ def plot(system, filename=defaultname, fmt=None, a=None,
         raise ValueError("The requested functionality requires the "
                          "Python Image Library (PIL)")
 
-    # symbols and lines may be constant, functions or dicts
-    # Here they are wrapped with a function
+    # symbols and lines may be constant or functions
+    # Here they are wrapped as a function
 
     if hasattr(symbols, "__call__"):
         fsymbols = symbols
-    elif hasattr(symbols, "__getitem__"):
+    elif is_builder and hasattr(symbols, "__getitem__"):
         fsymbols = lambda x : symbols[x.group]
     else:
         fsymbols = lambda x : symbols
 
     if hasattr(lines, "__call__"):
         flines = lines
-    elif hasattr(lines, "__getitem__"):
+    elif is_builder and hasattr(lines, "__getitem__"):
         flines = lambda x, y : (lines[x.group, y.group] if (x.group, y.group)
                                 in lines else lines[y.group, x.group])
     else:
@@ -640,7 +629,7 @@ def plot(system, filename=defaultname, fmt=None, a=None,
         flsymbols = fsymbols
     elif hasattr(lead_symbols, "__call__"):
         flsymbols = lead_symbols
-    elif hasattr(lead_symbols, "__getitem__"):
+    elif is_builder and hasattr(lead_symbols, "__getitem__"):
         flsymbols = lambda x : lead_symbols[x.group]
     else:
         flsymbols = lambda x : lead_symbols
@@ -649,15 +638,14 @@ def plot(system, filename=defaultname, fmt=None, a=None,
         fllines = flines
     elif hasattr(lead_lines, "__call__"):
         fllines = lead_lines
-    elif hasattr(lines, "__getitem__"):
+    elif is_builder and hasattr(lines, "__getitem__"):
         fllines = lambda x, y : (lead_lines[x.group, y.group]
                                  if (x.group, y.group) in lead_lines
                                  else lead_lines[y.group ,x.group])
     else:
         fllines = lambda x, y : lead_lines
 
-
-    #Figure out the extent of the system
+    # Figure out the extent of the system
     nsites = 0
     first = True
     for site in iterate_all_sites(system, len(lead_fading)):
@@ -804,7 +792,7 @@ def plot(system, filename=defaultname, fmt=None, a=None,
 
     # The lines for the hoppings
 
-    for site1, site2 in iterate_system_hoppings(system):
+    for site1, site2 in iterate_scattreg_hoppings(system):
         line = flines(site1, site2)
 
         if line is not None:
@@ -846,7 +834,7 @@ def plot(system, filename=defaultname, fmt=None, a=None,
                                      0.5 * add(pos(site1), pos(site2)), dist)
     # the symbols for the sites
 
-    for site in iterate_system_sites(system):
+    for site in iterate_scattreg_sites(system):
         symbol = fsymbols(site)
 
         if symbol is not None: