initialization and import of old files from 8.13 in fall 2006

8 files changed, 692 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..52e4e61
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+*.pyc
+*.pyo
diff --git a/calib_plots.py b/calib_plots.py
new file mode 100644
index 0000000..a905d91
--- /dev/null
+++ b/calib_plots.py
@@ -0,0 +1,32 @@
+
+def find_peaks(x,num=8):
+    d = x;
+    #figure();
+    #plot(d);
+
+    peaks = r_[0.:num];
+    i = 0;
+    mmm = max(d)
+    for p in range(0,int(num)):
+        while (d[i] < mmm*.2):
+            i += 1;
+        peaks[p] = 0.;
+        for n in range(-1,10):
+            peaks[p] += r_[i+n]*d[i+n]
+        peaks[p] = peaks[p]/sum(d[i-1:i+9])
+        i += 13
+    return peaks;
+
+
+def fit_peaks(p):
+    l=lambda a,x:a[0]*x+a[1]
+    print p; 
+    plotfit(fittodata_linear(r_[0.:size(p)],p,[25.,1.],err=.5*ones(size(p)),ret_all=True));
+
+def plot_mcacalib():
+    figure();
+    plot(fmca_calib[0:450]);
+    txt = text(0.22,.75,"delta-T/",
+            horizontalalignment='center',
+            verticalalignment='center',
+            transform = ax.transAxes)
diff --git a/importdata.py b/importdata.py
new file mode 100644
index 0000000..a72f71e
--- /dev/null
+++ b/importdata.py
@@ -0,0 +1,101 @@
+# @author: bryan newbold
+# @email: bnewbold@mit.edu
+# @date: nov 2006
+
+import shutil, os, scipy, pylab;
+
+def dothing(arg,thisdir,files):
+    global all_data
+    for fname in files:
+        if os.path.isfile(thisdir + '/' + fname):
+            if (fname.find('.') == -1):
+                try:
+                    dat = pylab.load(thisdir + '/' + fname)
+                    all_data[fname.split('.')[0]] = dat
+                    print "Loaded (no type): %s/%s" % (thisdir,fname);
+                except:
+                    print "Tried, but didn't load (no type):  %s/%s" % (thisdir,fname);
+            elif fname.endswith('.TKA'):
+                try:
+                    dat = pylab.load(thisdir + '/' + fname)
+                    #print all_data
+                    all_data[fname.split('.')[0]] = dat
+                    print "Loaded .TKA: %s/%s" % (thisdir,fname);
+                    #print "all_data has %d" % len(all_data)
+                except:
+                    print "Tried, but didn't load .TKA: %s/%s" % (thisdir,fname);
+            elif fname.endswith('.txt'):
+                try:
+                    dat = pylab.load(thisdir + '/' + fname,skiprows=1)
+                    #print all_data
+                    all_data[fname.split('.')[0]] = dat
+                    print "Loaded .txt: %s/%s" % (thisdir,fname);
+                    #print "all_data has %d" % len(all_data)
+                except:
+                    print "Tried, but didn't load .txt: %s/%s" % (thisdir,fname);
+            elif fname.endswith('.Spe!!!'):
+#                try:
+                    dat = pylab.load(thisdir+'/'+fname,skiprows=12);
+                    # also have to skip 15 from the end...
+                    print "trying: " + fname;
+#                    dat = scipy.io.read_array(thisdir + '/' + fname,lines=((14,-40,),))
+#                    print len(dat)
+#                    print dat
+                    dat = scipy.array(dat);
+                    all_data[fname.split('.')[0]] = dat
+                    print "Loaded .Spe: %s/%s" % (thisdir,fname);
+#                except:
+                    print "Tried, but didn't load .Spe: %s/%s" % (thisdir,fname);
+            else:
+                print "Didn't load %s/%s" % (thisdir,fname);
+        else:
+            print "Wasn't a file: %s/%s" % (thisdir,fname);
+
+def import_data(outfile='data.pydat',root='.',saveit=True,recurse=False):
+    """Imports all that yummah files..."""
+    global all_data
+    if saveit:
+        ofile = file(outfile,'w');
+        ofile.close()
+    all_data = {}
+    if recurse: 
+        os.path.walk(root,dothing,'')
+    else:
+        dothing('',root,os.listdir(os.getcwd()));
+    if saveit: 
+        print "Saving %d objects to %s" % (len(all_data), outfile);
+        scipy.io.objsave(ofile.name,globals(),all_data);
+    g = globals()
+    for a in all_data.keys():
+        g[a.replace('-','_')] = all_data[a]
+
+def load_imported_data(infile='data.pydat'):
+    global all_data
+    scipy.io.objload(infile,globals())
+    g=globals()
+    for a in all_data.keys():
+        if(a.startswith('11')):
+            g['d3' + a.replace('-','_')[9:]]=all_data[a]
+        else:
+            g[a.replace('-','_')]=all_data[a]
+
+def dorenamething(arg,thisdir,files):
+    for fname in files:
+        if os.path.isfile(thisdir + '/' + fname):
+            if fname.endswith('.Spe'):
+                try:
+                    shutil.copy(thisdir+'/'+fname,thisdir + '/' + fname[:-4]);
+                    print "Copied .Spe: %s/%s" % (thisdir,fname);
+                except:
+                    print "Tried, but didn't load .Spe: %s/%s" % (thisdir,fname);
+            else:
+                print "Didn't load %s/%s" % (thisdir,fname);
+        else:
+            print "Wasn't a file: %s/%s" % (thisdir,fname);
+
+def rename_data_files(root='.',ext='.Spe',recurse=False):
+    """Only works for .Spe..."""
+    if recurse: 
+        os.path.walk(root,dorenamething,'')
+    else:
+        dorenamething('',root,os.listdir(os.getcwd()));
diff --git a/plotting.py b/plotting.py
new file mode 100644
index 0000000..5f33dd6
--- /dev/null
+++ b/plotting.py
@@ -0,0 +1,8 @@
+
+import pylab;
+
+justplot = lambda x: pylab.plot(x[:,0],x[:,1]);
+justplotdot = lambda x: pylab.plot(x[:,0],x[:,1],'o');
+
+justplot_x = lambda x: pylab.plot(range(len(x)),x);
+justplotdot_x = lambda x: pylab.plot(range(len(x)),x,'o');
diff --git a/quick_tol.py b/quick_tol.py
new file mode 100644
index 0000000..176ba28
--- /dev/null
+++ b/quick_tol.py
@@ -0,0 +1,78 @@
+
+import pylab
+import scipy
+from pylab import *
+from scipy import *
+
+
+#_ip.magic("cd ~/jlab/muon/scipy/")
+
+
+d_tol_weekend = pylab.load('../data/bnewbold_solar_weekend2/SOLAR000')
+d_tol_weekend += pylab.load('../data/bnewbold_solar_weekend2/SOLAR001')
+d_tol_weekend += pylab.load('../data/bnewbold_solar_weekend2/SOLAR002')
+d_tol_weekend += pylab.load('../data/bnewbold_solar_weekend2/SOLAR003')
+d_tol_weekend += pylab.load('../data/bnewbold_solar_weekend2/SOLAR004')
+d_tol_weekend += pylab.load('../data/bnewbold_solar_weekend2/SOLAR005')
+d_tol_weekend += pylab.load('../data/bnewbold_solar_weekend2/SOLAR006')
+d_tol_weekend += pylab.load('../data/bnewbold_solar_weekend2/SOLAR007')
+d_tol_weekend += pylab.load('../data/bnewbold_solar_weekend2/SOLAR008')
+
+
+fitfunc = lambda p,x: p[0]*e^(-x/p[1])
+errfunc = lambda p,x,y: fitfunc(p,x) - y
+
+p0 = c_[1.0, 20.]
+y = frange(1,1024)
+fitfunc = lambda p,x: p[0]*e**(-x/p[1])
+errfunc = lambda p,x,y: fitfunc(p,x) - y
+p1,success = optimize.leastsq(errfunc, p0.copy(), args=(d_tol_weekend, y))
+p0 = c_[1.0, 20.]
+p0 = c_[50.0, 100.]
+p1,success = optimize.leastsq(errfunc, p0.copy(), args=(d_tol_weekend, y))
+p0 = c_[50.0, 100.]
+success
+p1
+errfunc = lambda p,x,y: fitfunc(p,x) - y
+p1,success = optimize.leastsq(errfunc, p0.copy(), args=(y, d_tol_weekend))
+plot(fitfunc(p1, frange(1,1024)), 'r')
+plot(d_tol_weekend, '.')
+legend(('Exponential Fit', 'Weekend #2 TOL data'))
+txt = text(380, 35, 'x0 = 32.555...\n T=171.1719...',)
+txt.set_bbox({"facecolor":"white","linewidth":1.0,"pad":16.0})
+
+title("Rough Fit of Weekend TOL Data")
+xlabel("Time (in bins...)")
+ylabel("Number of muon decays")
+
+
+fitfunc = lambda p,x: p[0]*e**(-x/p[1])+p[2]
+p0 = c_[50.0, 100.,10.]
+p1,success = optimize.leastsq(errfunc, p0.copy(), args=(y, d_tol_weekend))
+figure()
+plot(d_tol_weekend, '.')
+plot(fitfunc(p1, frange(1,1024)), 'r')
+title("Rough Fit of Weekend TOL Data, with +C")
+xlabel("Time (in bins...)")
+ylabel("Number of muon decays")
+t = text(380, 35, 'x0 = %e\n T=%e\nC=%e' % (p1[0],p1[1],p1[2]))
+t.set_bbox({"facecolor":"white","linewidth":1.0,"pad":16.0})
+success
+
+
+legend(('Weekend #2 TOL data', 'Exponential Fit'))
+
+
+fitfunc = lambda p,x: p[0]*e**(-x/p[1])+p[2]
+p0 = c_[50.0, 100.,10.]
+p1,success = optimize.leastsq(errfunc, p0.copy(), args=(y[13:1024], d_tol_weekend[13:1024]))
+figure()
+plot(d_tol_weekend, '.')
+plot(fitfunc(p1, frange(1,1024)), 'r')
+title("Rough Fit of Weekend TOL Data, with +C, minus outliers")
+xlabel("Time (in bins...)")
+ylabel("Number of muon decays")
+t = text(380, 35, 'x0 = %e\n T=%e\nC=%e' % (p1[0],p1[1],p1[2]))
+t.set_bbox({"facecolor":"white","linewidth":1.0,"pad":16.0})
+legend(('Weekend #2 TOL data', 'Exponential Fit'))
+success
diff --git a/stats.py b/stats.py
new file mode 100644
index 0000000..3315968
--- /dev/null
+++ b/stats.py
@@ -0,0 +1,196 @@
+# bnewbold_stats.py
+# helper functions!
+# @author: bryan newbold
+# @email: bnewbold@mit.edu
+# @date: oct 2006
+
+import scipy, pylab
+
+def floored_val(x,min=1.):
+    if pylab.size(x) > 1:
+        ret = scipy.array(pylab.copy.copy(x));
+        for i in range(0,pylab.size(x)):
+            ret[i] = floored_val(x[i],min=min);
+        return ret;
+    if(abs(x) < min):
+        return min;
+    else: 
+        return abs(x);
+
+def floored_root(x,min=1.):
+    return scipy.sqrt(floored_val(x,min=min));
+
+def rchisqr(tp,ffunc,tx,ty,terr):
+    """finds reduced chi squared for a function compared to data"""
+    
+    if not (pylab.size(tx) == pylab.size(ty) == pylab.size(terr)):
+        print "ERROR: input matrix mismatch... (rchisqr)";
+        return;
+    try:
+        ffunc(tp,tx[-1]);
+    except:
+        print "problem with passed function and/or parameters!";
+    try:
+        eta_err = (ty - ffunc(tp,tx))**2 / (terr**2)
+    except:
+        print "i think you have a zero error in there...";
+        eta_err = (ty - ffunc(tp,tx))**2 / (terr**2)
+    return sum(eta_err)/(pylab.size(tx)-pylab.size(tp));
+
+
+def rchisqr_poisson(tp,ffunc,tx,ty,floor=1.):
+    """finds reduced chi squared for a function compared to data,
+    assuming each data point has a poisson distribution and thus
+    an error of sqrt(val). The error value is floored at 1.0 to 
+    prevent stupid val=0 cases..."""
+    
+    if not (pylab.size(tx) == pylab.size(ty)):
+        print "ERROR: input matrix mismatch... (rchisqr_poisson)";
+        return;
+    terr = floored_root(tx,min=floor);
+    return rchisqr(tp,ffunc,tx,ty,terr);
+
+def fittodata(tx,ty,p0,ffunc,err=None,floor=1.,ret_all=False):
+    """does a fit and returns a dict of shit. 
+    ret_all=True gives more shit."""
+    if err == None:
+        err = floored_root(ty,min=floor);
+    if not (pylab.size(tx) == pylab.size(ty) == pylab.size(err)):
+        print "ERROR: input matrix mismatch... (fittodata)";
+        return;
+    
+    errfunc = lambda n,m,w: ffunc(n,m) - w;
+   
+    tx = scipy.array(tx);
+    #tx = tx.tolist(); tx = scipy.array(tx);
+    ty = scipy.array(ty);
+    #ty = ty.tolist(); ty = scipy.array(ty);
+    err = err.tolist(); err = scipy.array(err);
+    p0 = scipy.array(p0).tolist(); p0 = scipy.array(p0);
+    
+    
+    p1,cov_x,extra1,extra2,success = scipy.optimize.leastsq(errfunc, 
+                    p0.copy(), args=(tx,ty), full_output=1)
+    if size(p1) == 1:
+        p1 = [p1,]
+    print "DEBUG: done with fit"
+    print "p1: " + str(p1)
+    rcs = rchisqr(p1,ffunc,tx,ty,err);    
+    try:
+        p1_err = p1.copy()
+    except:
+        p1_err = [0,]
+    for i in range(scipy.size(p1)):
+        try:
+            scipy.size(cov_x);
+        except:
+            p1_err = array([0.0,]) 
+            p1 = array([p1,]) 
+            break;
+        if scipy.size(cov_x) == 1:
+            p1_err = array([scipy.sqrt(cov_x),]);
+            p1 = array([p1,]);
+        else:
+            p1_err[i] = scipy.sqrt(cov_x[i,i]);
+
+    if ret_all:
+        return {"params":p1, 
+                "param_error":p1_err,
+                "rchisqr":rcs,
+                "ffunc":ffunc,
+                "err":err,
+                "x":tx,
+                "y":ty,
+                "success":success};
+    else:
+        return {"params":p1,
+                "param_error":p1_err,
+                "rchisqr":rcs}
+
+def fittodata_gaussian(tx,ty,p0,err=None,floor=1.,ret_all=False):
+    """runs fittodata for a gaussian (first parameter is mean,
+    second is stddev, third is amplitude)"""
+    #gauss_func = lambda a,b: a[2]/(a[1] * sqrt(2*pi)) * exp(-.5 * ((b - a[0])/a[1])**2)
+    #thing = std(ty);
+    bigsum = scipy.integrate.simps(ty,x=tx)
+    gauss_func = lambda a,b: bigsum/(a[1] * scipy.sqrt(2*pi)) * scipy.exp(-.5 * ((b - a[0])/a[1])**2)
+    return fittodata(tx,ty,p0,gauss_func,err=err,floor=floor,ret_all=ret_all);
+
+def fittodata_linear(tx,ty,p0,err=None,floor=1.,ret_all=False):
+    """runs fittodata for a linear curve (first param is slope,
+    second is vertical offset)"""
+    lin_func = lambda a,b: a[0] * b + a[1]
+    return fittodata(tx,ty,p0,lin_func,err=err,floor=floor,ret_all=ret_all);
+
+def fittodata_poisson(tx,ty,p0,err=None,floor=.5,ret_all=False):
+    """runs fittodata for a poisson (first param is the mean"""
+    bigsum = scipy.integrate.simps(ty,x=tx)
+    if type(p0) == type(0.1):
+        p0 = [p0,]
+    poissfunc = lambda a,b: bigsum * scipy.power(a[0],b) / scipy.factorial(b) * scipy.exp(-1 * a[0])
+    #def poissfunc(a,b):
+    #    print "--------"
+    #    print a
+    #    print b
+    #    print bigsum
+    #    print scipy.power(a[0],b)
+    #    print scipy.exp(a[0])
+    #    print bigsum * scipy.power(a[0],b) / scipy.factorial(b) * scipy.exp(-1 * a[0])
+    #    return bigsum * scipy.power(a[0],b) / scipy.factorial(b) * scipy.exp(-1 * a[0])
+    return fittodata(tx,ty,p0,poissfunc,err=err,floor=floor,ret_all=ret_all);
+
+def fittodata_exp(tx,ty,p0,err=None,floor=1.,ret_all=False):
+    """runs fittodata for an exponential (first param is x0,
+    second param is tau, third param is vertical offset)"""
+    #exp_func = lambda a,b: a[0] * scipy.exp(- b/a[1]) + a[2]
+    exp_func = lambda a,b: a[0] * e ** (- b/a[1]) + a[2]
+    return fittodata(tx,ty,p0,exp_func,err=err,floor=floor,ret_all=ret_all);
+
+def fittodata_lorentz(tx,ty,p0,err=None,floor=1.,ret_all=False):
+    """runs fittodata for a lorentzian (first param is gamma,
+    second param is mean)"""
+    lorentz_func = lambda a,b: a[0]/(2.*pi*((a[1]-b)**2+(a[0]/2.)**2. )) 
+    return fittodata(tx,ty,p0,lorentz_func,err=err,floor=floor,ret_all=ret_all);
+
+def plotfit(indict,newfig=True):
+    if not indict.has_key('x'):
+        print "forgot to use ret_all on the fittodata?";
+        return;
+    p1 = indict['params'];
+    x = indict['x']
+    y = indict['y']
+    err = indict['err']
+    ffunc = indict['ffunc']
+    if newfig: pylab.figure();
+    pylab.errorbar(x,y,fmt='b.',yerr=err)
+    d=pylab.plot(x,y,'k.')
+    fx = x
+    if pylab.size(fx) < 100:
+        fx = r_[min(x):max(x):(max(x)-min(x))/150.]
+    f=pylab.plot(fx,ffunc(p1,fx),'r-')
+    #fy = indict['ffunc'](indict['params'],indict['x'])
+    pylab.title("Rough Fit of Data")
+    pylab.xlabel("X")
+    pylab.ylabel("Y")
+    p1 = indict['params']
+    p1_err = indict['param_error']
+    
+    ax = pylab.gca()
+    extra_txt = 'reduced chi squared: %g'\
+        % (indict['rchisqr'])
+    for i in range(scipy.size(p1)):
+        extra_txt += '\np[%i]: %g +/- %g' \
+            % (i,p1[i],p1_err[i]);
+    
+#    if(pylab.size(p1) == 2):
+#        extra_txt = 'reduced chi squared: %g\np[0]: %g\np[1]: %g' \
+#            % (indict['rchisqr'],p1[0],p1[1])
+#    else:
+#        extra_txt= 'reduced chi squared:%g\np[0]:%g, +/- %g\n p[1]:%g\np[2]:%g' \
+#            % (indict['rchisqr'],p1[0],p1[1],p1[2])
+    
+    t = pylab.text(.77, .75, extra_txt,
+        horizontalalignment='center',
+        verticalalignment='center',
+        transform = ax.transAxes)
+    pylab.legend((d,f),('Data', 'Fit'))
diff --git a/tof_data.py b/tof_data.py
new file mode 100644
index 0000000..5bc51a7
--- /dev/null
+++ b/tof_data.py
@@ -0,0 +1,72 @@
+# tof-data.py
+# Churn and display data from tof muon experiment
+# @author: bryan newbold
+# @email: bnewbold@mit.edu
+# @date: oct 2006
+
+from pylab import *;
+from scipy import *;
+
+data_file = '/home/bnewbold/jlab/muon/data/all_data.pydat';
+
+if not vars().has_key('fsolar2'):
+    io.objload(data_file,globals());
+
+def plot_tof(dat):
+    y = array(dat[0])
+    x = r_[0.:size(y)]
+    dat = array((dat[0],dat[1],dat[2]))
+
+    err = r_[0.:size(y)];
+    dmean = r_[0.:size(y)];
+    dmean_all = 0.;
+    err_all = 0.;
+    for i in range(0,size(y)):
+        i = int(i)
+        err[i] = std((dat[0,i],dat[1,i],dat[2,i]))
+        dmean[i] = mean((dat[0,i],dat[1,i],dat[2,i]))
+        dmean_all += dmean[i]*i
+    ax=gca()
+    errorbar(x,dmean,fmt='b',yerr=err,alpha=0.2)
+    plot(x,dmean,'k.')
+    dmean_all = dmean_all/sum(dmean)
+    err_all = sum(err)/size(err)
+    t = text(.22, .75, 'Mean: %g\nStdDev of Mean: %g' % (dmean_all,err_all),
+        horizontalalignment='center',
+        verticalalignment='center',
+        transform = ax.transAxes)
+    t.set_bbox({"facecolor":"white","linewidth":1.0,"pad":16.0})
+    xlabel('delta-Time (in \"bins units\")"');
+    ylabel('Counts');
+    return (dmean,err,dmean_all,err_all)
+
+quick_calib = lambda aD,aP: plotfit(fittodata_linear(r_[0.:aP],find_peaks(aD,num=aP),[1.,1.],ret_all=True))
+
+
+def plot_tof_more(dat,txt="blah",pos=0,style='b.'):
+    y = array(dat[0])
+    x = r_[0.:size(y)]
+    dat = array((dat[0],dat[1],dat[2]))
+
+    err = r_[0.:size(y)];
+    dmean = r_[0.:size(y)];
+    dmean_all = 0.;
+    err_all = 0.;
+    for i in range(0,size(y)):
+        i = int(i)
+        err[i] = std((dat[0,i],dat[1,i],dat[2,i]))
+        dmean[i] = mean((dat[0,i],dat[1,i],dat[2,i]))
+        dmean_all += dmean[i]*i
+    ax=gca()
+    er=errorbar(x,dmean,fmt=style,yerr=err,alpha=1.,label="_nolegend_")
+    plot(x,dmean,'k.',label="_nolegend_")
+    dmean_all = dmean_all/sum(dmean)
+    err_all = sum(err)/size(err)
+    t = text(.2, .90-.20*pos, '%s\nMean: %g\nStdDev of Mean: %g' % (txt,dmean_all,err_all),
+        horizontalalignment='center',
+        verticalalignment='center',
+        transform = ax.transAxes)
+    t.set_bbox({"facecolor":"white","linewidth":1.0,"pad":16.0})
+    xlabel('delta-Time (in \"bins units\")"');
+    ylabel('Counts');
+    return (dmean,err,dmean_all,err_all)
diff --git a/tof_mc.py b/tof_mc.py
new file mode 100644
index 0000000..2316c3c
--- /dev/null
+++ b/tof_mc.py
@@ -0,0 +1,203 @@
+# tol-mc.py
+# Monte-Carlo Simulation of a basic muon time-of-flight experiment
+# @author: bryan newbold
+# @email: bnewbold@mit.edu
+# @date: oct 2006
+
+from pylab import *;
+from scipy import *;
+
+# default number of tries to run...
+max_tries = 10**2
+
+# default geometry
+h = 100.0 # cm
+h_err = 1.0 # cm
+
+Tw = 15.5*2.54 # inches->cm
+Tw_err = .5*2.54
+Tl = 23.5*2.54
+Tl_err = 1.*2.54
+Tpmtl = 19.*2.54
+Tpmtl_err = 1.*2.54
+
+Bw = 42.
+Bw_err = 4./10 # mm->cm
+Bl = 53.
+Bl_err = 1.
+Bpmtl = 20.
+Bpmtl_err = 1.
+
+n_scint = 1.5 # index of refraction of scintillators
+C_ = 2.99792458*10**10 # speed of light in cm
+
+# shortcuts to random number generators
+ur = lambda: rand();
+gr = lambda: randn();
+
+aTx = lambda: ur()*Tl
+aTy = lambda: ur()*Tw
+#aphi = lambda: arccos( sqrt(ur()))
+aomega = lambda: ur()*2*pi
+aBx = lambda x,p,o,h: x+h*sin(o)*tan(p)
+aBy = lambda y,p,o,h: y+h*cos(o)*tan(p)
+#aP = lambda p: 100.*cos(p)
+aP = lambda p,mc: ((mc/2.99) + gr() * 0.02) * C_
+aTime = lambda x,y,pmt,w: n_scint/(C_)*sqrt((x+pmt)**2 + (w/2 - y)**2) 
+
+if not vars().has_key('phi_table'):
+    phi_table = zeros((1000,1),typecode='f')
+    for i in frange(0., pi/2,0.00001):
+        m = (2*i + sin(2*i))/pi
+        phi_table[int(1000*(m - (m%.001)))] = i
+    phi_table[-1] = pi/2
+
+def aphi():
+    n = ur();
+    return phi_table[int(1000*(n-(n%.001)))][0];
+
+def run(tries=max_tries, h=h, h_err=h_err,printstuff=True,mean_c=2.93):
+    num_events=0;
+    if(printstuff): print "Going to make %d tries..." % tries
+    events=zeros((tries,7),typecode='f')
+    while tries>0:
+        tries = tries-1;
+        Tx,Ty,phi,omega=aTx(),aTy(),aphi(),aomega();
+        Bx = aBx(Tx,phi,omega,h)
+        if((Bx < 0.) or (Bx > Bl)):
+            continue;
+        By = aBy(Ty + abs(Tw-Bw)/2,phi,omega,h)
+        if((By < 0.) or (By > Bw)):
+            continue;
+        P = aP(phi,mean_c)
+        dL = (h/cos(phi)); 
+        dT = -1 *aTime(Tx,Ty + abs(Tw-Bw)/2,Tpmtl,Tw) \
+             + aTime(Bx,By,Bpmtl,Bw) \
+             + dL/P; 
+        events[num_events] = (Tx,Ty,Bx,By,P,dT,dL);
+        num_events = num_events +1;
+    if num_events>0:
+        if(printstuff): print "... %d valid events!" % num_events;
+        return events[0:num_events];
+    else:
+        if(printstuff): print "No valid events!"
+        return [];
+
+def plot_events(events, xres=16,yres=16):
+    if(size(events[0]) != 7):
+        print "Passed matrix is no good...";
+        return;
+    figure();
+    # TOP part
+    subplot(211);
+    x = arange(0.,Tl,Tl/xres);
+    y = arange(0.,Tw,Tw/yres);
+    X,Y = meshgrid(x,y);
+    Z = zeros((xres,yres));
+    for i in events:
+        Z[int(i[1] * yres/Tw),int(i[0] * xres/Tl)] +=1;
+    contourf(X,Y,Z);
+    title("Distribution on Top Paddle (counts per %.3g cm^2)" \
+        % (Tl/xres * Tw/yres));
+    xlabel("Paddle Length (%g cm +/- %g cm)" % (Tl,Tl_err));
+    ylabel("Paddle Width (%g cm +/- %g cm)" % (Tw,Tw_err));
+    colorbar();
+
+    # BOTTOM part
+    subplot(212);
+    x = arange(0.,Bl,Bl/xres);
+    y = arange(0.,Bw,Bw/yres);
+    X,Y = meshgrid(x,y);
+    Z = zeros((xres,yres));
+    for i in events:
+        Z[int(i[3] * yres/Bw),int(i[2] * xres/Bl)] +=1;
+    contourf(X,Y,Z);
+    title("Distribution on Bottom Paddle (counts per %.3g cm^2)" \
+        % (Bl/xres * Tw/yres));
+    xlabel("Paddle Length (%g cm +/- %g)" % (Bl,Bl_err));
+    ylabel("Paddle Width (%g cm +/- %g cm)" % (Bw,Bw_err));
+    colorbar();
+
+def plot_events_3d(events):
+    if(size(events[0] != 6)):
+        print "Passed matrix is no good...";
+        return;
+    figure();
+    return;
+
+def plot_deltaT_h(heights=frange(33.,233,50.),num_events=1000,err=3,c=2.93):
+    
+    points = zeros((size(heights),6),typecode='f');
+    points[:,0] = list(heights[:]);
+    times = frange((int(err)));
+    dL = frange((int(err)));
+    actP = frange((int(err)));
+    for i in points:
+        for r in range(0,size(times)):
+            events = run(int(num_events),h=i[0],printstuff=False,mean_c=c);
+            times[int(r)] = mean(events[:,5]);
+            dL[int(r)] = mean(events[:,6]);
+            actP[int(r)] = mean(events[:,4]);
+        i[1] = mean(times);
+        i[2] = std(times);
+        i[3] = mean(dL);
+        i[4] = mean(actP)
+        i[5] = std(actP)
+    
+    stuff = fittodata_linear(points[:,0],points[:,1],[1.,1.],err=points[:,2],ret_all=True)
+    plotfit(stuff)
+    v = stuff['params']
+    ax = gca()
+    title("Change in delta-T with average h");
+    xlabel("Average h in cm");
+    ylabel("Average delta-T");
+    txt = text(0.65,.22,"delta-T/h: %g \noffset: %g \nmean muon speed: %g cm/s\nmonte carlo input speed:%g +/- %g cm/s" 
+            % (v[0],v[1],1/v[0],mean(points[:,4]),mean(points[:,5])), 
+            horizontalalignment='center',
+            verticalalignment='center',
+            transform = ax.transAxes)
+    
+
+    stuff = fittodata_linear(points[:,3],points[:,1],[1.,1.],err=points[:,2],ret_all=True)
+    plotfit(stuff)
+    vc = stuff['params']
+    ax=gca()
+    title('Change in mean delta-T with change in mean L')
+    xlabel('Mean Length of flight')
+    ylabel('Mean time gap between events')
+    txt = text(0.65,.22,"delta-T/h: %g \noffset: %g \nmean muon speed: %g cm/s\nmean muon speed corrected: %g cm/s\nmonte carlo input speed:%g +/- %g cm/s" 
+            % (v[0],v[1],1/v[0],1/vc[0],mean(points[:,4]),mean(points[:,5])), 
+            horizontalalignment='center',
+            verticalalignment='center',
+            transform = ax.transAxes)
+
+def plot_L_vs_h(heights=r_[18.5:218.5:8.],num_events=400,err=2,c=2.93):
+    points = zeros((size(heights),5),typecode='f');
+    
+    for i in range(0,size(heights)):
+        points[i,0] = heights[i];
+
+    times = r_[0.:int(err)];
+    dL = r_[0.:int(err)];
+    for i in range(0,size(heights)):
+        for r in range(0,size(times)):
+            events = run(int(num_events),h=points[i,0],printstuff=False,mean_c=c);
+            times[int(r)] = mean(events[:,5]);
+            dL[int(r)] = mean(events[:,6]);
+        points[i,1] = mean(times);
+        points[i,2] = std(times);
+        points[i,3] = mean(dL);
+        points[i,4] = std(dL);
+
+    x = points[:,0].tolist(); x = array(x)
+    y = points[:,3].tolist(); y = array(y)
+    err = points[:,4].tolist(); err = array(err)
+
+    plotfit(fittodata_linear(x,y,[3.,-1.],err=err,ret_all=True))
+    
+    title("Change in mean L with h");
+    xlabel("h in cm");
+    ylabel("Mean L");
+    figure();
+    plot(x,y-x);
+    return (x,y,err);