%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% WORKING WITH THL SCAN DATA
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% In Pixelman, THL scans are used to find spectra.
% The script loads a series of images taken during a THL scan, in order to
% find the total counts at each THL value.
%
% The basic version just plots differential count rates vs. the THL value.
% More sophisticated work can involve curve fitting, and using data from
% X-ray tests and noise tests to calibrate the results rather than just
% using the THL value.
% We can also calculated the total no of hits, and no of charge shared
% events
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Open the data files
% Here, we assume we already have the following data:
% Mask - 256*256 array with pixels we want to ignore set to 0, others set to 1.
% FilenamePrefix - Initial part of the filename
% THLStart - Range of THL values in scan
% THLStop - As above
% NStop - No of separate scans done
kMax=THLStop-THLStart+1;
% Loop over the no. of data sets taken
for N=1:NStop
% Loop over the THL values
for k=1:kMax
% Get correct THL value, and insert it into the vector form
THL=THLStart+k-1;
THLvalues(k)=THL;
% Use concatenation to make file name - in our tests, we have
% 3D1_15keV_THLscan_10__THL_399.txt
if (N<10)
Filename=[FilenamePrefix '0' num2str(N) '_THL_' num2str(THL) '.txt'];
else
Filename=[FilenamePrefix num2str(N) '_THL_' num2str(THL) '.txt'];
end
Image=load(Filename);
% Apply mask, find the total no of counts. Put this into
% appropriate entry of data
Data_counts(k,N)=sum(Image(Mask));
end
end
% So, we now have the Data_counts with all the data sets as columns, and
% the THLvalues vector.
% Can call individual data sets as Data_counts(:,n)
% Construct the differential count rate from this
% The diff function will work down columns on our data sets
% Note that the length of the data set will be reduced by 1; need to
% shorten the THL vector
Data_counts_diff = diff(Data_counts);
THLvalues_short = THLvalues;
THLvalues_short(end) = [];
% Can now find the average counts and the average differential count
% Mean function works down columns, so need to use double-transpose
% approach to average rows
Data_counts_mean = (mean(Data_counts.')).';
Data_counts_diff_mean = (mean(Data_counts_diff.')).';
% Likewise, can calculate the std deviation, and the standard error (std.
% dev divided by the square root of the no of measurements)
Data_counts_diff_stddev = (std(Data_counts_diff.')).';
Data_counts_diff_stderror = Data_counts_diff_stddev/(NStop.^0.5);
% Try a basic plot
figure;
% Spectrum from differential count rate
errorbar(THLvalues_short,Data_counts_diff_mean,Data_counts_diff_stddev,'*');
title('Differential count rate vs THL');
xlabel('THL setting');
ylabel('Rate of change in counts');
%axis([340 400 0 3e5]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CURVE FITTING - SELECTING THE DATA
% When finding the peak, we want to be able to ignore both the noise peak
% and the charge-shared signal
% We will take the point on the falling edge of the peak with the steepest
% gradient as our cut-off point
% So, start by finding the second derivative of our data, smoothing as
% appropriate, then finding the minimum.
% Take the derivative of the spectrum
Data_counts_diff2_mean = diff(Data_counts_diff_mean);
THLvalues_short2 = THLvalues_short;
THLvalues_short2(end) = [];
%figure;
%plot(THLvalues_short2,Data_counts_diff2_mean);
% At this point, the values fluctuate a lot. With better data, there might
% be less fluctuation. Will do a moving average, with a
% span of about 11. This will have the additional advantage of tending to
% keep us away from both the good signal edge and the noise edge.
% The moving average smooth uses similar syntax to lowess. The "span"
% argument here says how many data points to use - must be odd (e.g. span 3
% means we use x(n-1), x(n), x(n+1)
Falling_edge_smooth_span=11;
Data_counts_diff2_mean_smoothed = smooth(Data_counts_diff2_mean,Falling_edge_smooth_span,'moving');
figure;
plot(THLvalues_short2,Data_counts_diff2_mean_smoothed);
% Next, we find the minimum value of this, and the corresponding THL, using
% the index
% When doing this, might want to ignore part of the range, to avoid noise
UpperPeakEdgeLimit=400;
LowerPeakEdgeLimit=200;
% For finding harmonics, need to select range carefully
%UpperPeakEdgeLimit=180;
%LowerPeakEdgeLimit=100;
Data_counts_diff2_select=((THLvalues_short2LowerPeakEdgeLimit));
% So, set the appropriate range of the double-diff data to zero
Data_counts_diff2_mean_smoothed=Data_counts_diff2_mean_smoothed.*(Data_counts_diff2_select.');
[Min_Data_counts_diff2_mean_smoothed,PeakEdgeIndex]=min(Data_counts_diff2_mean_smoothed);
% This grabs the THL value of the steepest-falling point on the edge.
PeakEdgeTHL=THLvalues_short2(PeakEdgeIndex)
% If this all goes wrong, can insert the upper range by hand
%PeakEdgeTHL=400;
% Then, we can apply a Gaussian fit to the data, with a cut-off point set
% by this edge
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CURVE FIT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We want to exclude the data at THL values above the cut-off we've chosen.
% After saying we want default options for Gauss, we supply an "exclude"
% vector where the elements corresponding to the points we want to mask are
% set to 1.
% Create a vector containing 0s corresponding to the points where THL is less than or
% equal to the cut-off value, and 1s in the regions we want to exclude from
% the fit.
% If we have a lot of empty space below the peak, or perhaps harmonics,
% start at a certain point
DataStart=300;
SignalPeakFitExclude=((THLvalues_short>PeakEdgeTHL)|(THLvalues_shortNoiseEdgeUpperLimit)=0;
[NoiseEdgeCounts,NoiseEdgeIndex]=min(abs(Data_counts_NoiseEdge-20e6));
% This grabs the THL value of the steepest-falling point on the edge.
NoiseEdgeTHL=THLvalues_short(NoiseEdgeIndex)
%{
% The following section didn't work very well - ignore.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOISE PEAK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assuming the original noise peak is a Gaussian, then when we find the
% "spectrum" the new "peak" will actually be the derivative of a Gaussian.
% So, I need to be careful. It makes more sense for me to fit a Gaussian to
% the original data set, and then differentiate the result.
% For the current version, I'm not sure how the fit will cope if we have
% higher statistics in the peak, or if we have a large no of entries with
% charge sharing. The effect of the charge sharing and the remenants of the
% peak will carry different weight depending the size of the noise peak.
% We want to work with the original data, minus the actual signal
% First, we grab the values of the Gaussian fitted to the signal peak.
% Then, after we insert a leading zero (to compensate for how diff works)
% we subtract the results from the original data. This effectively allows
% us to integrate the signal peak and remove it.
% Using this, get an idea of how many counts are in the Gaussian (i.e. total no
% of counts from source)
TotalSignalCounts=CumulativeSignalPeakFit(end);
%Subtract the signal from the original data
Peaksub_Data_counts_mean=Data_counts_mean-CumulativeSignalPeakFit;
%figure;
%plot(THLvalues,Data_counts_mean,THLvalues,Peaksub_Data_counts_mean,THLvalues,cumsum(SignalPeakFitEval));
% Now that we've done this, we can actually make the fit to the noise
% When we're doing this, we don't want the incredibly high noise count
% region to dominate the results. So, when making the fit, we cut off the
% noise signal where the noise count rate exceeds the no of signal counts
% by a suitably high factor.
CountRateCutoff=5;
NoisePeakFitExclude=(Peaksub_Data_counts_mean>(TotalSignalCounts*CountRateCutoff))|(THLvalues(E0/2))&(keVvalues_short<(E0*3))))
TotalCountsAtHalfE0=sum(Data_counts_diff_mean((keVvalues_short>(E0/2))&(keVvalues_short<(E0*1.33))))
ChargeSharedCounts=TotalCountsAtHalfE0-TotalCountsSignalPeak
ProportionShared=ChargeSharedCounts/TotalCountsAtHalfE0
% Could potentially subtract the peak from the data, but no real need given
% that charge-sharing is uniform.