Fortunately this has a simple explanation, although I have to admit it's far from obvious...
Short version: increase the number of k-points for the DOS to 50 or 100 in the Z direction and all is fine. You may want this anyway, since the accuracy at 9 points is highly questionable.
Long version: ATK chooses a Gaussian broadening by default if the number of k-points is very low. This is good for molecules, but probably not very good for periodic structure. The breakpoint is 10 points - with time-reversal taken into account, meaning above 1,1,21 k-points it looks fine, since it uses the tetrahedron method instead. You could force the use of the tetrahedron method (at least for the log file, or custom plotting) by adding to the script
dos_tetra = density_of_states.tetrahedronSpectrum()
if processIsMaster:
for e,dos in zip(density_of_states.energies(), dos_tetra):
print e.inUnitsOf(eV),dos.inUnitsOf(eV**-1)
For up to 20 k-points in C, this will be different from nlprint(density_of_states), for 21 and above it will be the same.
See http://quantumwise.com/documents/manuals/latest/ReferenceManual/index.html/ref.densityofstates.html