Physics307L:People/Trujillo/LAB NOTEBOOK/071003

From OpenWetWare

Jump to: navigation, search

firstweek

Raw Data Files

10ms

Oct 03 2007     03:07:41 am     Elt: 000000 Seconds.  Real Time: 000000

ID: No spectrum identifier defined

Memory Size: 16384 Chls  Conversion Gain: 1024  Adc Offset: 0000 Chls



 Chn    Counts  ROI

   0,        0, 000

   1,        0, 000

   2,        0, 000

   3,        0, 000

   4,        0, 000

   5,        0, 000

   6,        0, 000

   7,        0, 000

   8,        0, 000

   9,        0, 000

  10,        0, 000

  11,        0, 000

  12,        0, 000

  13,        0, 000

  14,        0, 000

  15,        0, 000

  16,        0, 000

  17,        0, 000

  18,        0, 000

  19,        0, 000

  20,        0, 000

  21,        0, 000

  22,        0, 000

  23,        0, 000

  24,        0, 000

  25,        0, 000

  26,        0, 000

  27,        0, 000

  28,        0, 000

  29,        0, 000

  30,        0, 000

  31,        0, 000

  32,        0, 000

  33,        0, 000

  34,        0, 000

  35,        0, 000

  36,        0, 000

  37,        0, 000

  38,        0, 000

  39,        0, 000

  40,        0, 000

  41,        0, 000

  42,        0, 000

  43,        0, 000

  44,        0, 000

  45,        0, 000

  46,        0, 000

  47,        0, 000

  48,        0, 000

  49,        0, 000

  50,        0, 000

  51,        0, 000

  52,        0, 000

  53,        0, 000

  54,        0, 000

  55,        0, 000

  56,        0, 000

  57,        0, 000

  58,        0, 000

  59,        0, 000

  60,        0, 000

  61,        0, 000

  62,        0, 000

  63,        0, 000

  64,        1, 000

  65,        0, 000

  66,        0, 000

  67,        0, 000

  68,        0, 000

  69,        0, 000

  70,        0, 000

  71,        0, 000

  72,        0, 000

  73,        0, 000

  74,        0, 000

  75,        0, 000

  76,        0, 000

  77,        0, 000

  78,        0, 000

  79,        0, 000

  80,        0, 000

  81,        0, 000

  82,        0, 000

  83,        0, 000

  84,        0, 000

  85,        0, 000

  86,        0, 000

  87,        0, 000

  88,        0, 000

  89,        0, 000

  90,        0, 000

  91,        0, 000

  92,        0, 000

  93,        0, 000

  94,        0, 000

  95,        0, 000

  96,        0, 000

  97,        0, 000

  98,        0, 000

  99,        0, 000

 100,        0, 000

 101,        0, 000

 102,        0, 000

 103,        0, 000

 104,        0, 000

 105,        0, 000

 106,        1, 000

 107,        0, 000

 108,        0, 000

 109,        0, 000

 110,        0, 000

 111,        0, 000

 112,        0, 000

 113,        0, 000

 114,        0, 000

 115,        0, 000

 116,        0, 000

 117,        0, 000

 118,        0, 000

 119,        0, 000

 120,        0, 000

 121,        0, 000

 122,        0, 000

 123,        0, 000

 124,        0, 000

 125,        0, 000

 126,        0, 000

 127,        0, 000

 128,        0, 000

 129,        0, 000

 130,        0, 000

 131,        0, 000

 132,        0, 000

 133,        0, 000

 134,        0, 000

 135,        0, 000

 136,        0, 000

 137,        0, 000

 138,        0, 000

 139,        0, 000

 140,        0, 000

 141,        0, 000

 142,        0, 000

 143,        0, 000

 144,        0, 000

 145,        0, 000

 146,        1, 000

 147,        0, 000

 148,        0, 000

 149,        0, 000

 150,        0, 000

 151,        0, 000

 152,        0, 000

 153,        0, 000

 154,        0, 000

 155,        0, 000

 156,        0, 000

 157,        0, 000

 158,        0, 000

 159,        0, 000

 160,        0, 000

 161,        0, 000

 162,        0, 000

 163,        0, 000

 164,        0, 000

 165,        0, 000

 166,        0, 000

 167,        1, 000

 168,        0, 000

 169,        0, 000

 170,        0, 000

 171,        0, 000

 172,        0, 000

 173,        0, 000

 174,        0, 000

 175,        0, 000

 176,        0, 000

 177,        0, 000

 178,        0, 000

 179,        0, 000

 180,        0, 000

 181,        0, 000

 182,        0, 000

 183,        0, 000

 184,        0, 000

 185,        0, 000

 186,        0, 000

 187,        0, 000

 188,        0, 000

 189,        0, 000

 190,        0, 000

 191,        0, 000

 192,        0, 000

 193,        0, 000

 194,        0, 000

 195,        0, 000

 196,        0, 000

 197,        0, 000

 198,        0, 000

 199,        0, 000

 200,        0, 000

 201,        0, 000

 202,        0, 000

 203,        0, 000

 204,        0, 000

 205,        0, 000

 206,        0, 000

 207,        0, 000

 208,        0, 000

 209,        0, 000

 210,        0, 000

 211,        0, 000

 212,        0, 000

 213,        0, 000

 214,        2, 000

 215,        0, 000

 216,        0, 000

 217,        0, 000

 218,        1, 000

 219,        0, 000

 220,        0, 000

 221,        0, 000

 222,        0, 000

 223,        0, 000

 224,        2, 000

 225,        0, 000

 226,        0, 000

 227,        0, 000

 228,        0, 000

 229,        0, 000

 230,        0, 000

 231,        0, 000

 232,        0, 000

 233,        0, 000

 234,        0, 000

 235,        0, 000

 236,        0, 000

 237,        0, 000

 238,        0, 000

 239,        0, 000

 240,        0, 000

 241,        0, 000

 242,        0, 000

 243,        0, 000

 244,        0, 000

 245,        0, 000

 246,        0, 000

 247,        0, 000

 248,        0, 000

 249,        0, 000

 250,        0, 000

 251,        2, 000

 252,        0, 000

 253,        0, 000

 254,        0, 000

 255,        0, 000

==Octave script for analysis & plotting==
{{SJK comment|label=01:21, 22 October 2007 (CDT)|comment=I had never heard of Octave before, so found [http://en.wikipedia.org/wiki/GNU_Octave a description on wikipedia].  Are you happy with it?  Is it easy to use if you don't know matlab?  In any case, thanks for posting your analysis code!  I hope to have some time to check it out...}}
chop off the first ten lines of the data files before use

rename variable and filenames as necessary. Probably only needs minor tweeks so that it can run in matlab.
<pre><nowiki>#octave
#a sort of diary of commands to automate the data processing day to day
#some of octave's functions do not work in matlab,and vice 
#versa. Some functions may work, but not in the same way
#Most things work, though.

#Tomas Mondragon

load LAB3_D10.ASC;     #The load function is intended
load LAB3_D20.ASC;     #to read files saved by octave.
load LAB3_D40.ASC;     #Load can also read .mat files,
load LAB3_D80.ASC;     #but only versions saved by matlab
load LAB3_D100.ASC;    #version 4. It can also read text files.
load LAB3_D200.ASC;    #if no variable name info is in the file
load LAB3_D400.ASC;    #it names the variable after the file name
load LAB3_D800.ASC;
load LAB3_D1S.ASC;
load LAB3_D10S.ASC;

dwell10ms=LAB3_D10(:,2);     #load put what was on LAB3_D10.ASC into the 
dwell20ms=LAB3_D20(:,2);     #variable LAB3_D10. I only care about the second
dwell40ms=LAB3_D40(:,2);     #column of the data, so here I extract it to an
dwell80ms=LAB3_D80(:,2);     #appropriately named variable. 
dwell100ms=LAB3_D100(:,2);   
dwell200ms=LAB3_D200(:,2);   #The first column of LAB3_D10 is the channel number
dwell400ms=LAB3_D400(:,2);   #but a more appropriate name would be the bin index.
dwell800ms=LAB3_D800(:,2);   #The second column is a count of how many events fell
dwell1s=LAB3_D1S(:,2);       #into the bin. The third column isn't important for this
dwell10s=LAB3_D10S(:,2);     #experiment, but knowing what it is will allow one to take 
                             #advantage PCA3's region of interest feature. It just marks
                             #where one has marked ROI's. for example, 4 indicates the bin or channel 
                             #lies within ROI 4

#idea:use max(dwell10s) as max bin for all so grapfh are same x scale

#[y,x]=hist(data,bincenters) is a a function that sets up bins whose values are centered 
#at the values given by the vector bincenter and counts the number of times the values in 
#data fall into each of the bins. y contains the frequncy counts and x contains the 
#corresponding bin index. bar(x,y) will plot the histogram of data. In general, x=bincenters,
#so one can shorten [y,x]=hist(data,bincenters) to y=hist(data,bincenters) and use 
#bar(bincenters,y) to plot the same thing.

bincenters=0:max(dwell10s);

freq10ms=hist(dwell10ms,bincenters);
freq20ms=hist(dwell20ms,bincenters);
freq40ms=hist(dwell40ms,bincenters);
freq80ms=hist(dwell80ms,bincenters);
freq100ms=hist(dwell100ms,bincenters);
freq200ms=hist(dwell200ms,bincenters);
freq400ms=hist(dwell400ms,bincenters);
freq800ms=hist(dwell800ms,bincenters);
freq1s=hist(dwell1s,bincenters);
freq10s=hist(dwell10s,bincenters);

#plots. press any key to move to next plot
bar(bincenters,freq10ms)
title("frequency counts for dwell time=10ms")
xlabel("number of events occuring during dwell time")
ylabel("frequency count")
pause
replot
bar(bincenters,freq20ms)
title("frequency counts for dwell time=20ms")
%xlabel("number of events occuring during dwell time")
%ylabel("frequency count")
pause
replot
bar(bincenters,freq40ms)
title("frequency counts for dwell time=40ms")
%xlabel("number of events occuring during dwell time")
%ylabel("frequency count")
pause
replot
bar(bincenters,freq80ms)
title("frequency counts for dwell time=80ms")
%xlabel("number of events occuring during dwell time")
%ylabel("frequency count")
pause
replot
bar(bincenters,freq100ms)
title("frequency counts for dwell time=100ms")
%xlabel("number of events occuring during dwell time")
%ylabel("frequency count")
pause
replot
bar(bincenters,freq200ms)
title("frequency counts for dwell time=200ms")
%xlabel("number of events occuring during dwell time")
%ylabel("frequency count")
pause
replot
bar(bincenters,freq400ms)
title("frequency counts for dwell time=400ms")
%xlabel("number of events occuring during dwell time")
%ylabel("frequency count")
pause
replot
bar(bincenters,freq800ms)
title("frequency counts for dwell time=800ms")
%xlabel("number of events occuring during dwell time")
%ylabel("frequency count")
pause
replot
bar(bincenters,freq1s)
title("frequency counts for dwell time=1s")
%xlabel("number of events occuring during dwell time")
%ylabel("frequency count")
pause
replot
bar(bincenters,freq10s)
title("frequency counts for dwell time=10s")
%xlabel("number of events occuring during dwell time")
%ylabel("frequency count")

Compilation of plots

The plot that the script above produced with the data I obtained are shown in sequence below. As the dwell time increases, the most often occurring event count increases and the frequency vs. event count plots drift from something resembling a poisson distibution of a rarely occuring event to a the poisson distribution of a more common event, which resembles a gaussian.SJK 01:40, 22 October 2007 (CDT)
01:40, 22 October 2007 (CDT)These graphs are very snazzy, and your data certainly are great...however, it's a bit tough to look at the data, because I don't have a way to "play" and "pause" the animated GIF.  So, making the individual images available will definitely be necessary in your formal writeup.
01:40, 22 October 2007 (CDT)
These graphs are very snazzy, and your data certainly are great...however, it's a bit tough to look at the data, because I don't have a way to "play" and "pause" the animated GIF. So, making the individual images available will definitely be necessary in your formal writeup.

Image:Poissongraphsanimation.gif

I should point out that something interesting happened while we we counting events with a dwell of 200 ms. Our equipment recorded one instance where 27 events happened in that small dwell time. The mode of the data for dwell time=200ms was 0 events, and were were rarely getting counts more than or equal to 5.SJK 01:27, 22 October 2007 (CDT)
01:27, 22 October 2007 (CDT)Very interesting!  Bradley and Nik left the setup going from last Wednesday until today: 100 second bins, going like 5 days or something.  I wonder if they will also record a similar occurance (though 27 events may not be obvious in a 100 second bin)?By the way, you can easily calculate the probability of 27 events happening in a 200 ms window, given your Poisson fit.  You can then use this probability to argue about throwing out that data point when fitting.  (Also, you should directly state whether you throw the data out or not, and why.)
01:27, 22 October 2007 (CDT)
Very interesting! Bradley and Nik left the setup going from last Wednesday until today: 100 second bins, going like 5 days or something. I wonder if they will also record a similar occurance (though 27 events may not be obvious in a 100 second bin)?

By the way, you can easily calculate the probability of 27 events happening in a 200 ms window, given your Poisson fit. You can then use this probability to argue about throwing out that data point when fitting. (Also, you should directly state whether you throw the data out or not, and why.)

In Dr. Gold's lab manual, the poisson distribution is a model for the results of experiments that count random events that occur at a definite average rate and the frequency that one should expect to get a certain count number. If the expected average count is λ, the probability that one will obtain the result of k after performing the counting experiment is

f(k;\lambda)=\frac{\lambda^k e^{-\lambda}}{k!}

So, to fit my data, I will have to normalize my data and find a suitable λ for each data set.

First off, normalize. If the frequency count of an event count k is xk, find some normalizing constant N so that \sum_{k=\operatorname{min}\ k}^{\operatorname{max}\ k} \frac{x_k}{N}=1

Then, find λ such that\sum_{k=\operatorname{min}\ k}^{\operatorname{max}\ k}\left(f(k,\lambda)-\frac{x_k}{N} \right)^2 is a minimum

Duurrrrr... N = 256 because 256 tests were performed each time.

SJK 01:38, 22 October 2007 (CDT)
01:38, 22 October 2007 (CDT)It appears what you are doing here is a least squares fit of the Poisson distribution to the histogram, which is a great idea.  However, I'd have to think some more, but I think maybe you would want a weighting factor in the denominator, perhaps  (because the uncertainty in #counts is the sqrt(#counts)).  Also, though, via talking with Bradley about this lab, I found that in the wikipedia article on Poisson distribution, they show that the maximum likelihood fit is when lambda is set to equal the average counts for the data.  In your formal report, you should compare what you did with the much easier maximum likelihood method.
01:38, 22 October 2007 (CDT)
It appears what you are doing here is a least squares fit of the Poisson distribution to the histogram, which is a great idea. However, I'd have to think some more, but I think maybe you would want a weighting factor in the denominator, perhaps \sqrt{x_k} (because the uncertainty in #counts is the sqrt(#counts)). Also, though, via talking with Bradley about this lab, I found that in the wikipedia article on Poisson distribution, they show that the maximum likelihood fit is when lambda is set to equal the average counts for the data. In your formal report, you should compare what you did with the much easier maximum likelihood method.

Through an iterative process, I determine λ to be 0.0261205 for dwell=10ms. Ack, I won't go any further, I'll just stop at 6 significant digits for the others too.

  • 10ms,0.0261205
  • 20ms,0.0736898
  • 40ms,0.185892
  • 80ms,0.457154
  • 100ms,0.468385
  • 200ms,0.937453
  • 400ms,2.58578
  • 800ms,5.74562
  • 1s,7.27233
  • 10s,73.9022
Image:Poissongraphswithfits.gifSJK 01:44, 22 October 2007 (CDT)
01:44, 22 October 2007 (CDT)From my glimpses of the moving GIFs, the fits do look good.  In addition to being able to look at the individual images, the fits would probably also be easier to see if they were not also bars, but instead were data points, or curves...
01:44, 22 October 2007 (CDT)
From my glimpses of the moving GIFs, the fits do look good. In addition to being able to look at the individual images, the fits would probably also be easier to see if they were not also bars, but instead were data points, or curves...

how good are these fits, since I decided to stop at 6 sig figs (quite arbitrarily). I suppose I shall quantify this by taking the average of how far off the actual data is from the fit.

I think the formula is Error=\frac{\sqrt{\sum_{k=\operatorname{min}\ k}^{\operatorname{max}\ k}\left(f(k,\lambda)-\frac{x_k}{N} \right)^2}}{\sqrt{N}\sqrt{N-1}}

  • 10ms, 0.0261205, 5.45905*10^-6
  • 20ms, 0.0736898, 1.19310*10^-4
  • 40ms, 0.185892, 3.28259*10^-4
  • 80ms, 0.457154, 5.95464*10^-4
  • 100ms, 0.468385, 3.15888*10^-4
  • 200ms, 0.937453, 5.23153*10^-4
  • 400ms, 2.58578, 5.22450*10^-4
  • 800ms, 5.74562, 2.97039*10^-4
  • 1s, 7.27233, 2.59775*10^-4
  • 10s, 73.9022, 2.96150*10^-4
The poisson distribution fits count data of random events that happen at a definite average rate. Therefore \frac{\lambda}{dwell time}=some constant

A plot of λ vs. dwell time should fit on a line. to do this, I use Dr. Gold's linefit matlab function, modified for use with Octave. (www-hep.phys.unm.edu is down at the moment, maybe provide link when it comes back up or upload modded program here?) SJK 01:38, 22 October 2007 (CDT)

01:38, 22 October 2007 (CDT)Really cool way to analyze your data (checking lambda versus window width to be linear)!
01:38, 22 October 2007 (CDT)
Really cool way to analyze your data (checking lambda versus window width to be linear)!

Image:Lambdatimeplot.png

slope=7.39233+/-0.00003
intercept= -0.048099+/-0.000005
cov(m,b)=0.000000
correlation=0.999983
point error estimate=0.231612
chisq/ndf=265341
So, with our equipment set up the way it was (sorry, no data on this, oops!) we were observing events that occured at an average rate of 7.39 events per second.

Personal tools