diff --git a/analysis.py b/analysis.py
index 9a79085..bf7b07f 100755
--- a/analysis.py
+++ b/analysis.py
@@ -7,6 +7,7 @@
 from StringIO import StringIO
 import sys
 from astropy.io import fits
+from astropy.table import Table, Column
 
 
 
@@ -19,7 +20,7 @@ def main(folder, quiet=0):
 Parameters
 ----------
 folder  String   Specifiecs the folder where the files are
-quiet   boolean  =1 surpresses all standard output
+quiet   boolean  =1 suppresses all standard output
 
 Returns
 -------
@@ -58,8 +59,19 @@ def main(folder, quiet=0):
             x = -2.5*(np.log10(data['c%s' % color1]/64130) - np.log10(data['c%s' % color2]/7140))
             y = -2.5*(np.log10(data['c%s' % color2]/7140))
 
-            n = sum(np.logical_and( (y > -10./3. * (x-1.) + 10.), np.logical_and(max_mag < y, y < min_mag)))
+            
+            sel = np.logical_and( (y > -10./3. * (x-1.) + 10.), np.logical_and(max_mag < y, y < min_mag))
+            n = sum(sel)
+            t = Table(hdulist[1].data)
+            if 'sel' in t.columns:
+                t.remove_column('sel')
+            t.add_column(Column(name='sel', data=sel.astype('int')))
+            
+            hdulist[1].data = np.array(t)
             tmp, av, apera, age = fil.split('_')
+            #hdulist = fits.open('%s/%s' %(folder,fil), 'write')
+            #hdulist.writeto('%s/%s' %(folder,fil), clobber=True)
+            fits.update('%s/%s' %(folder,fil), np.array(t), ext = 1, clobber=True)
             out.append([Decimal(av), Decimal(apera), Decimal(age), n])
 
     #writing obtained data to "folder/__expected_number"
diff --git a/execute.py b/execute.py
index 3cb49eb..1c662e1 100644
--- a/execute.py
+++ b/execute.py
@@ -20,7 +20,7 @@ def main(quiet = False):
 
 Parameters
 ----------
-quiet   boolean  if true surpresses all standard output
+quiet   boolean  if true suppresses all standard output
 
 
 Returns
diff --git a/plot.py b/plot.py
index a4b119b..e6f0d26 100644
--- a/plot.py
+++ b/plot.py
@@ -232,3 +232,31 @@ def init():
 
     return avs, aperas, ages, sfr, headers, data
 
+
+def histogram(folder, av, apera, age, color1 = "I4", color2 = "M1"):
+    
+    hdulist = fits.open('%s/%s' %(folder,'sim_%03d_%06d_%09d'%(av,apera,age)))
+    data = hdulist[1].data
+
+
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+
+
+    pdf, bins, patches = ax.hist(data['mass'], range=(-.5,1.5))
+    dat = data['mass'][data['sel'] == 1]
+    pdf, bins, patches = ax.hist(dat, range=(-.5,1.5))
+    #ax.set_yscale('log')
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/setup.py b/setup.py
index e9fc211..8ba4993 100644
--- a/setup.py
+++ b/setup.py
@@ -23,7 +23,7 @@ def main(quiet=False):
 
 Parameter
 ---------
-quiet  boolean  if true surpresses all standard output
+quiet  boolean  if true suppresses all standard output
 '''
     if quiet:
         output_stream = StringIO()
diff --git a/starformation.py b/starformation.py
index bcd4a28..e16494b 100755
--- a/starformation.py
+++ b/starformation.py
@@ -28,7 +28,7 @@ def main(massfunction = 0, starformationhistory = 0, A_v = 10.0, sfr = .01, aper
 maxage    float     age of the star formation site, sfr is assumed to be constant
 distance  float     distance to the simulated starformation site
 appendix  String    sets the outputfilename, default is the starting time (via time.time())
-quiet     boolean   if true (=1) surpresses all standard output
+quiet     boolean   if true (=1) suppresses all standard output
 precise   boolean   if true (=1) sample single star till expected mass reached, else calculate
                     expected number and sample as an array
 
@@ -178,18 +178,6 @@ def g(x):
     newfluxes = np.asarray(newfluxes)
     output = np.vstack([np.asarray(output).transpose(), matches, fluxes, newfluxes]).transpose()
 
-    # creating the output file
-    #head = ['#', 'age', 'mass', 'model']
-    #for mod in models:
-        #head.append('flux %s' % mod)
-    #for mod in models:
-        #head.append('corrected_flux %s' % mod)
-    #f = open('out/%s' % appendix, 'w')
-    #f.write( ','.join(head)+'\n' )
-    ##np.savetxt(f, output)
-    #np.savetxt(f, output)
-    #f.close()
-
     # create table
     # data table
     t = Table()
@@ -200,6 +188,7 @@ def g(x):
         t.add_column(Column(name='%s' % models[i], data=output[:,4+i]))
     for i in range(len(models)):
         t.add_column(Column(name='c%s' % models[i], data=output[:,4+len(models)+i]))
+   # t.add_column(Column(name='sel', data=zeros(len(models[0]))))
     # head table
     header = Table()
     header.add_column(Column(name='AV', data = [A_v]))