diff --git a/xray/azav.py b/xray/azav.py
index 699aef4..0f05c0e 100644
--- a/xray/azav.py
+++ b/xray/azav.py
@@ -128,8 +128,7 @@ def doFolder(folder,files='*.edf*',nQ = 1500,force=False,mask=None,
   # which poni file to use:
   ai = _getAI(poni,folder)
 
-  files = glob.glob("%s/%s"%(folder,files))
-  files.sort()
+  files = utils.getFiles(folder,files)
   if saved is not None:
     files = [f for f in files if utils.getBasename(f) not in saved["files"]]
 
@@ -140,8 +139,8 @@ def doFolder(folder,files='*.edf*',nQ = 1500,force=False,mask=None,
     elif mask is not None:
       mask = pyFAIread(mask).astype(bool)
       
-    data   = np.empty( (len(files),nQ),dtype=np.float32 )
-    err    = np.empty( (len(files),nQ),dtype=np.float32 )
+    data   = np.empty( (len(files),nQ) )
+    err    = np.empty( (len(files),nQ) )
     for ifname,fname in enumerate(files):
       img = pyFAIread(fname)
       q,i,e = pyFAI1d(ai,img,mask=mask,npt_radial=nQ)
@@ -149,7 +148,7 @@ def doFolder(folder,files='*.edf*',nQ = 1500,force=False,mask=None,
       err[ifname]  = e
       if saveChi:
         chi_fname = utils.removeExt(fname) + ".chi"
-        utils.saveTxt(chi_fname,q,i,e,info=pyFAI_dict(ai),overwrite=True)
+        utils.saveTxt(chi_fname,q,i,e,info=pyFAI_info(ai),overwrite=True)
 
     files = [ utils.getBasename(f) for f in files ]
     files = np.asarray(files)
@@ -223,30 +222,6 @@ def pyFAI_find_center(img,psize=100e-6,dist=0.1,wavelength=0.8e-10,**kwargs):
   print("Final values: (in pixels) %.3f %.3f"%(xc,yc))
   return ai
 
-def plotdata(q,data,plot=True,showTrend=True,title=None,clim='auto'):
-  if not (plot or showTrend): return
-  if clim == 'auto': clim = np.nanpercentile(data,(1.5,98.5))
-  one_plot = showTrend or plot
-  two_plot = showTrend and plot
-  if one_plot and not two_plot:
-    fig,ax = plt.subplots(1,1)
-  if two_plot:
-    fig,ax = plt.subplots(1,2,sharey=True)
-  ax = np.atleast_1d(ax)
-  if showTrend:
-    plt.sca(ax[0])
-    plt.pcolormesh(np.arange(data.shape[0]),q,data.T)
-    plt.xlabel("image number, 0 being older")
-    plt.ylabel(r"q ($\AA^{-1}$)")
-    plt.clim( *clim )
-  if plot:
-    if showTrend:
-      ax[1].plot(data.mean(axis=0),q)
-    else:
-      ax[0].plot(q,data.mean(axis=0))
-  if (plot or showTrend) and title is not None:
-    plt.title(title)
-
 def average(fileOrFolder,delays=slice(None),scale=1,norm=None,returnAll=False,plot=False,
   showTrend=False):
   data = utils.data_storage(fileOrFolder)
@@ -268,7 +243,7 @@ def average(fileOrFolder,delays=slice(None),scale=1,norm=None,returnAll=False,pl
   if isinstance(norm,np.ndarray):
     i = i/norm[:,np.newaxis]
   title = "%s %s" % (fileOrFolder,str(delays))
-  plotdata(q,i*scale,showTrend=showTrend,plot=plot,title=title)
+  utils.plotdata(q,i*scale,showTrend=showTrend,plot=plot,title=title)
   if returnAll:
     return q,i.mean(axis=0)*scale,i
   else:
diff --git a/xray/dataReduction.py b/xray/dataReduction.py
index 8e6fbdb..dad29f0 100644
--- a/xray/dataReduction.py
+++ b/xray/dataReduction.py
@@ -38,6 +38,7 @@ def subtractReferences(i,idx_ref, useRatio = False):
       # normal reference for an on chi, the weighted average
       iref[_i]      = weight_before*ref_before + weight_after*ref_after
     if _i>=idx_ref_after: _ref += 1
+    log.debug("For image %d : %d-%d"%(_i,idx_ref_before,idx_ref_after))
   if useRatio:
     i /= iref
   else:
@@ -59,9 +60,10 @@ def averageScanPoints(scan,data,isRef=None,lpower=None,useRatio=False,\
       funcForEveraging: is usually np.nanmean or np.nanmedian. it can be any 
          function that support axis=0 as keyword argument
 """
-
+  data = data.astype(np.float)
   if isRef is None: isRef = np.zeros( data.shape[0], dtype=bool )
   assert data.shape[0] == isRef.shape[0]
+
   # subtract reference only is there is at least one
   if isRef.sum()>0:
     data = subtractReferences(data,np.argwhere(isRef), useRatio=useRatio)
@@ -158,7 +160,7 @@ def errorFilter(data,threshold=4):
     idx = data.err[iscan] > threshold*median_err
     data.data[iscan][idx] = \
       np.nanmedian( data.dataInScanPoint[iscan][:,idx],axis=0 )
-#    data.err[iscan][idx] = median_err
+    data.err[iscan][idx] = median_err
   return data
   
 def saveTxt(folder,data,delayToStr=True,info="",**kw):
@@ -173,7 +175,7 @@ def saveTxt(folder,data,delayToStr=True,info="",**kw):
     # try retreiving info on chi2
     try:
       chi2_0  = data.chi2_0[iscan]
-      info_delay = [ " # rep_num , chi2_0 ", ]
+      info_delay = [ "# rep_num : chi2_0 ", ]
       for irep,value in enumerate(chi2_0):
         info_delay.append( "# %d : %.3f" % (irep,value))
       info_delay = "\n".join(info_delay)
diff --git a/id9.py b/xray/id9.py
similarity index 71%
rename from id9.py
rename to xray/id9.py
index ba21b02..c51965c 100644
--- a/id9.py
+++ b/xray/id9.py
@@ -4,9 +4,9 @@ log.basicConfig(level=log.INFO)
 import os
 import collections
 import numpy as np
-from .xray import azav
-from .xray import dataReduction
-from .xray import utils
+from . import azav
+from . import dataReduction
+from . import utils
 
 storage_extension = ".npz"
 
@@ -44,15 +44,25 @@ def doFolder_azav(folder,nQ=1500,force=False,mask=None,saveChi=True,
   return azav.doFolder(folder,files="*.edf*",nQ=nQ,force=force,mask=mask,
     saveChi=saveChi,poni=poni,storageFile=storageFile,diagnostic=diag)
 
-def doFolder_dataRed(folder,monitor=None,funcForEveraging=np.nanmean,errFilter=True):
-  data = utils.data_storage(folder)
-  scan = data.delays
-  q    = data.q
-  i    = data.data
-  diffs = dataReduction.calcTimeResolvedSignal(scan,i,q=q,reference="min",
-          monitor=monitor,funcForEveraging=funcForEveraging)
+def doFolder_dataRed(folder,storageFile='auto',monitor=None,
+                     funcForEveraging=np.nanmean,errFilter=True):
+
+  if storageFile == 'auto' : storageFile = folder +  "/" + "pyfai_1d" + storage_extension
+
+  # read azimuthal averaged curves
+  data = utils.data_storage(storageFile)
+
+  
+  # calculate differences
+  diffs = dataReduction.calcTimeResolvedSignal(data.delays,data.data,q=data.q,\
+          reference="min",monitor=monitor,funcForEveraging=funcForEveraging)
+
+  # filter if asked so
   if funcForEveraging == np.nanmean and errFilter: 
     diffs = dataReduction.errorFilter(diffs)
+  
+  # save txt and npz file
   dataReduction.saveTxt(folder,diffs,info=data.pyfai_info)
   diffs.save(folder + "/" + "diffs" + storage_extension)
-  return diffs
+
+  return data,diffs
diff --git a/xray/storage.py b/xray/storage.py
index 4106039..ff49b04 100644
--- a/xray/storage.py
+++ b/xray/storage.py
@@ -72,6 +72,7 @@ def dictToNpz(npzFile,d): np.savez(npzFile,**d)
 
 def read(fname):
   extension = os.path.splitext(fname)[1]
+  log.info("Reading storage file %s"%fname)
   if extension == ".npz": 
     return npzToDict(fname)
   elif extension == ".h5":
@@ -81,6 +82,7 @@ def read(fname):
 
 def save(fname,d):
   extension = os.path.splitext(fname)[1]
+  log.info("Saving storage file %s"%fname)
   if extension == ".npz": 
     return dictToNpz(fname,d)
   elif extension == ".h5":
diff --git a/xray/utils.py b/xray/utils.py
index 28116d1..0f23755 100644
--- a/xray/utils.py
+++ b/xray/utils.py
@@ -19,6 +19,14 @@ except ImportError:
 _time_regex =      re.compile( "(-?\d+\.?\d*(?:ps|ns|us|ms)?)")
 _timeInStr_regex = re.compile("_(-?\d+\.?\d*(?:ps|ns|us|ms)?)")
 
+def getFiles(folder,basename="*.edf*"):
+  files = glob.glob(folder + "/" + basename)
+  files.sort()
+  return files
+
+def getEdfFiles(folder):
+  return getFiles(folder,basename="*.edf*")
+
 def getDelayFromString(string) :
   match = _timeInStr_regex_regex.search(string)
   return match and match.group(1) or None
@@ -51,6 +59,8 @@ def timeToStr(delay,fmt="%+.0f"):
     ret = fmt % (delay*1e12) + "fs"
   elif 1e-18 <= a_delay < 1e-15: 
     ret = fmt % (delay*1e12) + "as"
+  else:
+    ret = str(delay) +"s"
   return ret
 
 def removeExt(fname):
@@ -61,6 +71,80 @@ def removeExt(fname):
 def getBasename(fname):
   return removeExt(os.path.basename(fname));
 
+def plotdata(q,data,t=None,plot=True,showTrend=True,title=None,clim='auto'):
+  if not (plot or showTrend): return
+  if t is None: t = np.arange(data.shape[0])
+  if clim == 'auto': clim = np.nanpercentile(data,(1.5,98.5))
+  one_plot = showTrend or plot
+  two_plot = showTrend and plot
+  if one_plot and not two_plot:
+    fig,ax = plt.subplots(1,1)
+  if two_plot:
+    fig,ax = plt.subplots(1,2,sharey=True)
+  ax = np.atleast_1d(ax)
+  if showTrend:
+    plt.sca(ax[0])
+    plt.pcolormesh(t,q,data.T)
+    plt.xlabel("image number, 0 being older")
+    plt.ylabel(r"q ($\AA^{-1}$)")
+    plt.clim( *clim )
+  if plot:
+    if showTrend:
+      ax[1].plot(data.mean(axis=0),q)
+    else:
+      ax[0].plot(q,data.mean(axis=0))
+  if (plot or showTrend) and title is not None:
+    plt.title(title)
+
+
+def plotdiffs(q,diffs,t,select=None,err=None,showErr=False,cmap=plt.cm.jet):
+  # this selection trick done in this way allows to keep the same colors when 
+  # subselecting (because I do not change the size of diffs)
+  if select is not None:
+    indices = range(*select.indices(t.shape[0]))
+  else:
+    indices = range(len(t))
+  lines = []
+  for idiff in indices:
+    color = cmap(idiff/(len(diffs)-1))
+    label = timeToStr(t[idiff])
+    kw = dict( color = color, label = label )
+    if err is not None and showErr:
+      line = plt.errorbar(q,diffs[idiff],err[idiff],**kw)[0]
+    else:
+      line = plt.plot(q,diffs[idiff],**kw)[0]
+    lines.append(line)
+
+  fig    = plt.gcf() 
+  legend = plt.legend()
+  plt.grid()
+  plt.xlabel(r"q ($\AA^{-1}$)")
+  # we will set up a dict mapping legend line to orig line, and enable
+  # picking on the legend line
+  lined = dict()
+  for legline, origline in zip(legend.get_lines(), lines):
+    legline.set_picker(5)  # 5 pts tolerance
+    lined[legline] = origline
+
+  def onpick(event):
+    # on the pick event, find the orig line corresponding to the
+    # legend proxy line, and toggle the visibility
+    legline = event.artist
+    origline = lined[legline]
+    vis = not origline.get_visible()
+    origline.set_visible(vis)
+    # Change the alpha on the line in the legend so we can see what lines
+    # have been toggled
+    if vis:
+        legline.set_alpha(1.0)
+    else:
+        legline.set_alpha(0.2)
+    fig    = plt.gcf() 
+    fig.canvas.draw()
+
+  fig.canvas.mpl_connect('pick_event', onpick)
+
+
 def saveTxt(fname,q,i,e=None,headerv=None,info=None,overwrite=True):
   """ Write data to file 'fname' in text format.
       Inputs: