From 81a0c260d321da151e244548188d13a963bf111f Mon Sep 17 00:00:00 2001 From: Bryan Newbold Date: Mon, 11 Feb 2008 21:40:22 -0500 Subject: initialization and import of old files from 8.13 in fall 2006 --- .gitignore | 2 + calib_plots.py | 32 +++++++++ importdata.py | 101 ++++++++++++++++++++++++++++ plotting.py | 8 +++ quick_tol.py | 78 ++++++++++++++++++++++ stats.py | 196 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ tof_data.py | 72 ++++++++++++++++++++ tof_mc.py | 203 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 692 insertions(+) create mode 100644 .gitignore create mode 100644 calib_plots.py create mode 100644 importdata.py create mode 100644 plotting.py create mode 100644 quick_tol.py create mode 100644 stats.py create mode 100644 tof_data.py create mode 100644 tof_mc.py 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); -- cgit v1.2.3