diff --git a/plot_tally.py b/plot_tally.py
index b695f50..63e0303 100644
--- a/plot_tally.py
+++ b/plot_tally.py
@@ -13,7 +13,7 @@ from pylab import *
 from matplotlib.colors import LogNorm
 
 COLORMAP='tab20c'
-SCALING=1.56e16/8.7374e15
+SCALING=1.0
 
 def read_tally(fname, use_tally=None):
   if fname.endswith('.gz'):
@@ -68,9 +68,16 @@ def plot_xy(outfile, res, zidx=5, Imin=1e-6, Imax=1e4, cticks=None):
   #fig.savefig(outfile)
   #pc=p.get_facecolors()
   #pc[:,3]=where(I[:,:,zidx].T.flatten()>0, 0.9-0.25*dI[:,:,zidx].T.flatten(), 0.)
-  ax.add_patch(Circle((0,0), radius=5.50, color='black', fill=False, lw=2))
-  ax.add_patch(Circle((0,0), radius=11.50, color='black', fill=False, lw=2))
-  ax.add_patch(Circle((0,0), radius=15.00, color='black', fill=False, lw=2))
+  ax.add_patch(Circle((-0.151,0.021), radius=5.50, color='black', fill=False, lw=2))
+  ax.add_patch(Circle((-0.151,0.021), radius=11.50, color='black', fill=False, lw=2))
+  ax.add_patch(Circle((-0.151,0.021), radius=15.00, color='black', fill=False, lw=2))
+  # shielding walls outside bunker
+  ax.add_line(Line2D([21.12-0.151, 21.12-0.151], [-1.5,1.5], color='black', lw=2))
+  ax.add_line(Line2D([14.845,19.426], [-0.541,-0.713], color='black', lw=2))
+  ax.add_line(Line2D([19.426,19.405], [-0.713,-1.248], color='black', lw=2))
+  ax.add_line(Line2D([19.405,21.205], [-1.248,-1.460], color='black', lw=2))
+  ax.add_line(Line2D([14.723, 22.740], [2.232, 3.084], color='black', lw=2))
+  ax.add_line(Line2D([14.801, 20.894], [-1.210, -1.724], color='black', lw=2))
   if outfile is not None:
     fig.savefig(outfile, transparent=False, dpi=300)
 
@@ -89,9 +96,12 @@ def plot_xz(outfile, res, yidx=5, Imin=1e-6, Imax=1e4, cticks=None):
   #fig.savefig(outfile)
   #pc=p.get_facecolors()
   #pc[:,3]=where(I[:,yidx].T.flatten()>0, 0.9-0.25*dI[:,yidx].T.flatten(), 0.)
-  ax.plot([sqrt(5.5**2-ypos**2), sqrt(5.5**2-ypos**2)], [z.min()/100., z.max()/100.], color='black', lw=2)
-  ax.plot([sqrt(11.5**2-ypos**2), sqrt(11.5**2-ypos**2)], [z.min()/100., z.max()/100.], color='black', lw=2)
-  ax.plot([sqrt(15.0**2-ypos**2), sqrt(15.0**2-ypos**2)], [z.min()/100., z.max()/100.], color='black', lw=2)
+  ax.plot([sqrt(5.5**2-ypos**2)-0.151, sqrt(5.5**2-ypos**2)-0.151], 
+          [z.min()/100., z.max()/100.], color='black', lw=2)
+  ax.plot([sqrt(11.5**2-ypos**2)-0.151, sqrt(11.5**2-ypos**2)-0.151], 
+          [z.min()/100., z.max()/100.], color='black', lw=2)
+  ax.plot([sqrt(15.0**2-ypos**2)-0.151, sqrt(15.0**2-ypos**2)-0.151], 
+          [z.min()/100., z.max()/100.], color='black', lw=2)
   if outfile is not None:
     fig.savefig(outfile, transparent=False, dpi=300)
 
@@ -112,26 +122,67 @@ def plot_yz(outfile, res, xidx=5, Imin=1e-6, Imax=1e4, cticks=None):
   if outfile is not None:
     fig.savefig(outfile, transparent=False, dpi=300)
 
+def get_name(fname):
+  # extract tag-name and version info from filename
+  stripped=os.path.basename(fname)[9:].split('.')[0]
+  tag, version=stripped.split('-', 1)
+  return tag, version
+  
+
 if __name__=='__main__':
   import sys, os
-  res=read_tally(sys.argv[1], use_tally=1)
+  if '-o' in sys.argv:
+    idx=sys.argv.index('-o')
+    sys.argv.pop(idx)
+    prefix=sys.argv.pop(idx)
+  else:
+    prefix='estia_dosemap'
+  print sys.argv[1]
+  if ':' in sys.argv[1]:
+    fname, tally=sys.argv[1].split(':',2)
+    tally=int(tally)
+  else:
+    fname=sys.argv[1]
+    tally=None
+  res=read_tally(fname, use_tally=tally)
   x,y,z,I,dI=res
   I*=SCALING
-  fpre=os.path.join('.', 'images', os.path.basename(sys.argv[1])[:-4])
+  fpath=os.path.join('.', 'images', get_name(fname)[0]+'-'+get_name(fname)[1])
+  if tally is not None:
+    fpath+='[%s]'%tally
+  if len(sys.argv)>2:
+    for fi in sys.argv[2:]:
+      print fi
+      if ':' in fi:
+        fi, tally=fi.split(':',2)
+        tally=int(tally)
+      else:
+        tally=None
+      resi=read_tally(fi, use_tally=tally)
+      xi,yi,zi,Ii,dIi=resi
+      Ii*=SCALING
+      I+=Ii
+      dI=sqrt(dI**2+dIi**2)
+      res=(x,y,z,I,dI)
+      fpath+='+'+get_name(fi)[1]
+      if tally is not None:
+        fpath+='[%s]'%tally
+  os.mkdir(fpath)
+  fpre=os.path.join(fpath, prefix)
   xi, yi, zi=res[3].shape
   for zidx in range(zi):
     print fpre+'_xy_%02i.png'%zidx
     figure(figsize=(12,8))
     plot_xy(fpre+'_xy_%02i.png'%zidx, res, zidx=zidx, 
-            Imin=1e-8, Imax=1e2, cticks=[1e-8, 1e-6, 1e-4, 1e-2, 1, 100])
+            Imin=1e-7, Imax=1e-2, cticks=[1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2])
   for yidx in range(yi):
     print fpre+'_xz_%02i.png'%yidx
     figure(figsize=(12,8))
     plot_xz(fpre+'_xz_%02i.png'%yidx, res, yidx=yidx, 
-            Imin=1e-8, Imax=1e2, cticks=[1e-8, 1e-6, 1e-4, 1e-2, 1, 100])
+            Imin=1e-7, Imax=1e-2, cticks=[1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2])
   for xidx in range(xi):
     print fpre+'_yz_%02i.png'%xidx
     figure(figsize=(8,10))
     plot_yz(fpre+'_yz_%02i.png'%xidx, res, xidx=xidx, 
-            Imin=1e-8, Imax=1e2, cticks=[1e-8, 1e-6, 1e-4, 1e-2, 1, 100])
+            Imin=1e-7, Imax=1e-2, cticks=[1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2])