summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore2
-rw-r--r--calib_plots.py32
-rw-r--r--importdata.py101
-rw-r--r--plotting.py8
-rw-r--r--quick_tol.py78
-rw-r--r--stats.py196
-rw-r--r--tof_data.py72
-rw-r--r--tof_mc.py203
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);
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-21 11:27:20 +0000