User:Alexander T. J. Barron/Notebook/PHYC 307L Lab Notebook/Poisson notes
Poisson Statistics Lab
Partner: Justin Muehlmeyer
Set UpSJK 09:45, 20 October 2008 (EDT)
Searching for the cosmic radiation signals:
Using a Tektronix 485 Oscilloscope (analog) and ORTEC model # 266 PM (photo-multiplier) base.
Dynode output on PMT (photo-multiplier tube) to oscilloscope Input Voltage -1500 V PMT outside of lead brick cave
- 5 mV/division
- 1 ms/division
- AC triggering: negative slope (internal)
We observe occasional peaks per second that are distinguishable from noise. They change so fast that we can't actually measure the amplitude of their peaks. This is our presumed cosmic radiation.
Using an amplifier:
- We addded a ORTEC amplifier to our board, and connected the output of the diode to the input of the amplifyer. We wanted to test the signal coming out of the amplifier with the oscilloscope before connecting it into the computer. We ran into many difficulties even finding our cosmic radiation signal with the amplifier, and then again while bypassing it.
We had extensive difficulties in trying to get a signal of the cosmic radiation using the amplifier.
New amplifier: Canberra PAD 814.
Switched from dynode to the anode on the PMT, with no effect.
Switched to new PMT - Tracor Northern Model # TN1222. Way clearer, greater magnitude cosmic radiation spikes with just the oscilloscope. Input ~-1400 V. Positive V output.
Still no results from the amplifier, so...
New Oscilloscope: Tektronix TDS 1002 (digital).
"Derivative signal" from the amplifier.
Set up: PMT to PREAMP IN, PREAMP switch to IN, AMP OUT to oscilloscope (DC coupling, rising slope, 36.0 mV threshold)
Hooking Multichannel Analyzer Card (MCA) up to the software:
Software: PCA3 in MSDOS
We decided to first take our dwell time to be 1 second.
Other passes: dwell time of 1, 10, 100 milliseconds, and finally 2 seconds.
Saved to ASCII files individually-- the .txt files are here:
- Image:JAND1MS.txt - 1 ms dwell
- Image:JAND10MS.txt - 10 ms dwell
- Image:JAND100M.txt - 100 ms dwell
- Image:JAND1S.txt - 1 s dwell
- Image:JAND2S.txt - 2 s dwell
The goal of this lab is to see whether the frequencies of cosmic ray interactions with a scintillator over several different dwell times follows a Poisson distribution. They should by definition, assuming that the average rate of cosmic ray detection does not change over time. To test this hypothesis, I will first create separate histograms of event counts, calculate the mean of each, and use said means to fit a Poisson distribution over the histogram for each dwell time. Through suggestion from Dr. Koch, I will compare vs. a gaussian function for each case as well.
Here are the histograms:
Here is the MATLAB code:
%% Poisson Distribution Analysis % Alexander Barron % Junior Lab Fall 2008 close all, clear all; %% Load .txt files M1ms_bin = dlmread('JAND1MS.txt',',',[0 0 511 0]); % bin vector 1ms dwell M1ms_val = dlmread('JAND1MS.txt',',',[0 1 511 1]); % value vector 1ms dwell M10ms_bin = dlmread('JAND10MS.txt',',',[0 0 511 0]); M10ms_val = dlmread('JAND10MS.txt',',',[0 1 511 1]); M100ms_bin = dlmread('JAND100M.txt',',',[0 0 511 0]); M100ms_val = dlmread('JAND100M.txt',',',[0 1 511 1]); M1s_bin = dlmread('JAND1S.txt',',',[0 0 511 0]); M1s_val = dlmread('JAND1S.txt',',',[0 1 511 1]); M2s_bin = dlmread('JAND2S.txt',',',[0 0 511 0]); M2s_val = dlmread('JAND2S.txt',',',[0 1 511 1]); %% sort into frequencies for histograms d1ms_count = zeros(1,50); % Here I am creating empty d10ms_count = zeros(1,50); % vectors in which to put d100ms_count = zeros(1,50); % my manipulated data. d1s_count = zeros(1,50); d2s_count = zeros(1,50); for c=1:5; switch c; case 1; valmat = M1ms_val; % generic value and countmat = d1ms_count; % count matrices case 2; valmat = M10ms_val; countmat = d10ms_count; case 3; valmat = M100ms_val; countmat = d100ms_count; case 4; valmat = M1s_val; countmat = d1s_count; case 5; valmat = M2s_val; countmat = d2s_count; end; for j=1:512; % Here I scan my data for % desired frequencies for k=0:50; % and place them in % generic matrices. if valmat(j) == k; countmat(k+1) = countmat(k+1) + 1; end end end switch c; case 1; d1ms_count = countmat; % transformation of case 2; % generic matrices d10ms_count = countmat; % into appropriately- case 3; % named ones d100ms_count = countmat; case 4; d1s_count = countmat; case 5; d2s_count = countmat; end; clear valmat, clear countmat; end; %% Plot freq = linspace(0,49,50); scrsz = get(0,'ScreenSize'); figure('Position',[1 scrsz(4)/1.5 scrsz(3)/1.25 scrsz(4)/1.4]); for f=1:5; switch f; case 1; data = d1ms_count; titl = '1 ms Dwell Time'; case 2; data = d10ms_count; titl = '10 ms Dwell Time'; case 3; data = d100ms_count; titl = '100 ms Dwell Time'; case 4; data = d1s_count; titl = '1 s Dwell Time'; case 5; data = d2s_count; titl = '2 s Dwell Time'; end; subplot(1,5,f), bar(freq,data,'c'); ylim([0 525]); xlim([-2 50]); set(gca,'XTick',[0 10 20 30 40 50]); ylabel('Frequency of Events'); xlabel('Events'); title(titl); clear data; end;
As we increase the dwell time, each histogram begins to look more like a Gaussian function. Perhaps at very large dwell times, a gaussian will fit the data perfectly. It would have been nice to take more data to test this... maybe as my final lab report...
Calculating the Mean
A good method of determining the mean can be found from the Thermodynamics class notes of Dr. Geremia which I post here:
Where we determine our mean λ to be the rate multiplied by the time. This comes from the definition of the mean, which is the summation of the (# of events) X (Probability of those events) over all values. Geremia shows that when simplified, this summation becomes (average rate per dwell time) X ( dwell time).
I will use this formula to calculate the average for each dwell time:
We took 512 passes per dwell time. ki is the count per bin from the MCA data. Add this code to calculate the mean per dwell time:
%% Calculate Mean (lambda) n = 512; format long; d1mslambda = 1/n * sum(M1ms_val) d10mslambda = 1/n * sum(M10ms_val) d100mslambda = 1/n * sum(M100ms_val) d1slambda = 1/n * sum(M1s_val) d2slambda = 1/n * sum(M2s_val)
|Dwell Time Δt||Frequency Count λ|
|1 ms||λ = 0.013671875000000|
|10 ms||λ = 0.126953125000000|
|100 ms||λ = 1.294921875000000|
|1 s||λ = 13.488281250000000|
|2 s||λ = 25.361328125000000|
This makes sense, as our dwell time increases, our average number of events increases. These values for λ seem to roughly match the means for the above histograms. To make a probability density bar graph, just divide each histogram bin by the number of iterations.
Add this code at the end of the last block:
%% Transform to Probability Density Graph probd1ms_count = d1ms_count./512; probd10ms_count = d10ms_count./512; probd100ms_count = d100ms_count./512; probd1s_count = d1s_count./512; probd2s_count = d2s_count./512; %% Plot Probability Density Graph scrsz = get(0,'ScreenSize'); figure('Position',[1 scrsz(4)/1.5 scrsz(3)/1.25 scrsz(4)/1.4]); for f=1:5; switch f; case 1; data = probd1ms_count; titl = '1 ms Dwell Time'; case 2; data = probd10ms_count; titl = '10 ms Dwell Time'; case 3; data = probd100ms_count; titl = '100 ms Dwell Time'; case 4; data = probd1s_count; titl = '1 s Dwell Time'; case 5; data = probd2s_count; titl = '2 s Dwell Time'; end; subplot(1,5,f), bar(freq,data,'c'); ylim([0 1]); xlim([-2 50]); set(gca,'XTick',[0 10 20 30 40 50]); ylabel('Frequency of Events'); xlabel('Events'); title(titl); clear data; end;
I tested to see if these functions integrated to 1. They do. They are, indeed, probability density graphs.
Now for the error on lambda. I got this formula for the error on lambda from Tomas Mondragon's notebook entry last year:
Where the frequency count of an event count k is xk. N is 512 in this case.
So, once again, tack this code on the end of the MATLAB m-file:
%% Error on Lambda for j=1:5; switch j; case 1; lambda = d1mslambda; count = d1ms_count; t0 = 0; % minimum event count tf = 1; % maximum event count case 2; lambda = d1mslambda; count = d10ms_count; t0 = 0; tf = 2; case 3; lambda = d100mslambda; count = d100ms_count; t0 = 0; tf = 4; case 4; lambda = d1slambda; count = d1s_count; t0 = 7; tf = 20; case 5; lambda = d2slambda; count = d2s_count; t0 = 17; tf = 34; end; errorsum = 0; for t=t0:tf; errorsum = errorsum + ((lambda^t*exp(-lambda)/factorial(t) - count(t+1)/512))^2; end; error = ((errorsum))^.5/(512*511)^.5; switch j; case 1; d1ms_lerror = error case 2; d10ms_lerror = error case 3; d100ms_lerror = error case 4; d1s_lerror = error case 5; d2s_lerror = error end; clear error, clear errorsum, clear count, clear t0, clear tf; end;
So now my full results with error are here:
|Dwell Time Δt||Frequency Count λ||λerror|
|1 ms||λ = 0.013671875000000||4.059741317917446e-07|
|10 ms||λ = 0.126953125000000||2.975294357739427e-04|
|100 ms||λ = 1.294921875000000||4.620595463752644e-04|
|1 s||λ = 13.488281250000000||3.259048960726763e-04|
|2 s||λ = 25.361328125000000||2.425961354924050e-04|
This error seems plausible.SJK 09:20, 20 October 2008 (EDT)
Fitting with Poisson
Here I use the lambda values to create Poisson plots over my data distributions. I also created a "range" overlay of Poisson to contrast the extremes in error I found, but the plots with error accounted for are indistinguishable from those without (at least at this scope). It appears that the Poisson distribution is pretty spot-on for 1ms and 10ms dwell, then matches the mean of the data without modeling its true shape. The data range of counts is decidedly narrower than the Poisson range at 1-2s.
Here is the code:
%% Fitting with Poisson figure('Position',[1 scrsz(4)/1.5 scrsz(3)/1.25 scrsz(4)/1.4]); for j=1:5; switch j; case 1; lambda = d1mslambda; % <- here I define simple lambda w/out error case 2; lambda = d10mslambda; case 3; lambda = d100mslambda; case 4; lambda = d1slambda; case 5; lambda = d2slambda; end; for t=1:50; Poissonvec(t) = lambda^(t-1)*exp(-lambda)/factorial(t-1); end; switch j; case 1; Pd1msvec = Poissonvec; data = probd1ms_count; titl = '1 ms Dwell Time'; case 2; Pd10msvec = Poissonvec; data = probd10ms_count; titl = '10 ms Dwell Time'; case 3; Pd100msvec = Poissonvec; data = probd100ms_count; titl = '100 ms Dwell Time'; case 4; Pd1svec = Poissonvec; data = probd1s_count; titl = '1 s Dwell Time'; case 5; Pd2svec = Poissonvec; data = probd2s_count; titl = '2 s Dwell Time'; end; subplot(1,5,j), bar(freq,data,'c'), hold on, plot(freq,Poissonvec,'.r'); ylim([0 1]); xlim([-2 50]); set(gca,'XTick',[0 10 20 30 40 50]); legend('Data','Poisson'); ylabel('Frequency of Events'); xlabel('Events'); title(titl); clear lambda, clear Poissonvec, clear data, clear titl; end;
Fitting with Gaussian
This is the Gaussian Function, according to the Normal Distribution Wikipedia Article:
where σ is the standard deviation, and μ is the expected value, or mean.
This is the relationship between the standard deviation and the mean, according to the Standard Deviation Wikipedia Article:
- Alexander T. J. Barron 18:36, 26 October 2008 (EDT): OK, looking back, I don't know what I was doing. Above is the sample standard deviation. The standard error would be the sample standard dev. divided by root N.
I tried fitting with the Gaussian function manually, as before with Poisson, but it was a disaster. I obtained these results with the Distribution Fitting Tool in MATLAB:
SJK 09:24, 20 October 2008 (EDT)
The tool obviously uses different bin widths and count ranges from mine when generating its histograms, but it's easy to see the same trend as I obtained before. Clearly, the Gaussian is a horrible fit for low dwell time. Now I'll fit both Gaussian and Poisson together with the tool:
|Gaussian vs. Poisson Fit