%% Program to read spectrometric profiles in .csv format and output graphs.
% Written by Olivia Leavitt and Luke Brown
% Nuclear Science IQP, Fukushima I Team: Radiation in PPG Project

%% Read all .csv files from the folder of the directory and place each in a cell array
files = dir('*.csv');
datacell = cell(10, length(files));
for fileindex = 1:length(files)
    file = files(fileindex).name;
    data = csvread(file,7,0);
    mass = csvread(file,2,2, [2 2 2 2]);
    [~,category,~] = xlsread(file, 1, 'C4');
    datacell{7,fileindex} = mass;
    datacell{10,fileindex} = category;
    data(:, 4) = [];
    datacell{1,fileindex} = file;
    datacell{2,fileindex} = data;
    figure
    semilogy(data(:,2), data(:,3));
    title_str = datacell{1, fileindex};
    title_str = regexprep(title_str, '_', ' ');
    title(title_str);
    xlabel('Energy (kEV)');
    ylabel('Counts');
    grid
end
grid
title('Raw Data');
xlabel('Energy (kEV)');
ylabel('Number of Counts');

%% Determine whether each file name contains the label 24 hours and assigns divisors (number of seconds of count) in the cell array
for fileindex = 1:length(files)
    is_24hr = strfind(datacell{1,fileindex}, '24hr');
    if isempty(is_24hr)
        datacell{3,fileindex} = 2*3600;
    else
        datacell{3,fileindex} = 24*3600;
    end
end

%% Normalize each spectrum by time and plot the normalized spectra
figure
for fileindex = 1:length(files)
    energies = datacell{2, fileindex};
    raw_counts = datacell{2, fileindex};
    counts_per_second = raw_counts(:,3) ./ datacell{3, fileindex};
    raw_counts(:,3) = counts_per_second;
    datacell{4,fileindex} = raw_counts;
    semilogy(energies(:,2), raw_counts(:,3));
    hold on;
end
title('Time-normalized Data');
xlabel('Energy (kEV)')
ylabel('Counts/second')
grid

%% Create "control" spectrum from all blank counts
for fileindex = 1:length(files)
    if isempty(strfind(datacell{1,fileindex}, 'control')) == 0
        is_control = 1;
    elseif isempty(strfind(datacell{1,fileindex}, 'blank')) == 0
        is_control = 1;
    elseif isempty(strfind(datacell{1,fileindex}, 'background')) == 0
        is_control = 1;
    elseif isempty(strfind(datacell{1,fileindex}, 'Background')) == 0
        is_control = 1;
    elseif isempty(strfind(datacell{1,fileindex}, 'Blank')) == 0
        is_control = 1;
    else
        is_control = 0;
    end
    datacell{5, fileindex} = is_control;  
end

num_controls = 0;
for fileindex = 1:length(files)
    is_control = datacell{5,fileindex};
    if is_control == 1
        num_controls = num_controls + 1;
    end
end
control_matrix = zeros(8192, num_controls + 2);
control_matrix(:, [1 2]) = energies(:, [1 2]);

counter = 2;
for fileindex = 1:length(files)
    is_control = datacell{5, fileindex};
    if is_control == 1
        counter = counter + 1;
        data = datacell{4, fileindex};
        control_matrix(:, counter) = data(:,3);
    end
end

control_average_matrix = zeros(8192, 4);
control_average_matrix(:, [1 2]) = control_matrix(:, [1 2]);
values = control_matrix(:, 3:counter);
control_average_matrix(:, 3) = mean(values, 2);
control_average_matrix(:, 4) = std(values, 0, 2);

figure
semilogy(control_average_matrix(:,2), (control_average_matrix(:,3)-control_average_matrix(:,4)),'color', 'r');
hold on
semilogy(control_average_matrix(:,2), (control_average_matrix(:,3)+control_average_matrix(:,4)), 'color', 'r');
semilogy(control_average_matrix(:,2), control_average_matrix(:,3), 'color', 'g');
title('Background Data');
xlabel('Energy (kEV)');
ylabel('Counts/s');
grid

%% Plot each sample spectrum against the background spectrum and add the relative counts/second to the data cell
% then plot each spectrum's values above background

for fileindex = 1:length(files)
    if datacell{5, fileindex} ~= 1
        figure
        data = datacell{4, fileindex};
        data(:,3) = data(:,3) - control_average_matrix(:,3);
        semilogy(data(:,2), data(:,3));
        title_str = datacell{1, fileindex};
        title_str = regexprep(title_str, '_', ' ');
        title(title_str);
        xlabel('Energy (kEV)');
        ylabel('Counts/second over background');
        grid
        datacell{6, fileindex} = data;
    end
end

%% Divide each sample's data by the sample mass
for fileindex = 1:length(files)
    if datacell{5, fileindex} ~= 1
        data = datacell{6, fileindex};
        mass = datacell{7, fileindex};
        data(:, 3) = data(:,3)/mass;
        datacell{8, fileindex} = data;
        figure
        semilogy(data(:,2), data(:,3));
        title_str = datacell{1, fileindex};
        title_str = regexprep(title_str, '_', ' ');
        title(title_str);
        xlabel('Energy (kEV)');
        ylabel('Counts/second/g over background');
        grid
    end
end

%% Assign efficiency curve from BEGE Detector
data = datacell{2, 1};
energies = data(:, 2);
efficiency_curve = exp(1.475e2 - 1.499e2*log(energies) + (57.61*log(energies).^2) - (10.77*log(energies).^3) + (.9782*log(energies).^4) - (3.476e-2*log(energies).^5));
start_index = find(efficiency_curve<1, 1);
plot(energies(start_index:end), efficiency_curve(start_index:end));
title('MATLAB Calculated Efficiency Curve');
xlabel('Energy (kEV)');
ylabel('Efficiency');
ylim([0, 0.014]);
grid;

%% Divide each non-background count by the detector efficiency at that energy
% Assume, since each sample was on the detector face, that the detector gets 1/2 of the counts in the sample.

figure
for fileindex = 1: length(files)
    if datacell{5,fileindex} ~= 1
        data = datacell{8, fileindex};
        data(:,3) = data(:,3)./((1/2) .* efficiency_curve);
        datacell{9,fileindex} = data;
        %figure
        semilogy(data(:,2), data(:,3));
        hold on
        title_str = datacell{1, fileindex};
        title_str = regexprep(title_str, '_', ' ');
        title(title_str);
        xlabel('Energy (kEV)');
        ylim([0.01 1000]);
        ylabel('Bq/g over background');
        grid
    end
end
title('All Samples Data over Background');
        
%% For each category of experimental sample, create an average curve for the Bq/g over background of that data
% Include standard deviation curves

air_counter = 0;
feet_counter = 0;
gloves_counter = 0;
pants_counter = 0;
vest_counter = 0;
face_counter = 0;

for fileindex = 1:length(files)
    is_control = datacell{5, fileindex};
    if is_control ~= 1
        if strcmp(datacell{10, fileindex}, 'air') == 1
            air_counter = air_counter + 1;
        elseif strcmp(datacell{10,fileindex}, 'feet') == 1
            feet_counter = feet_counter + 1;
        elseif strcmp(datacell{10,fileindex}, 'gloves') == 1
            gloves_counter = gloves_counter + 1;
        elseif strcmp(datacell{10, fileindex}, 'pants') == 1
            pants_counter = pants_counter + 1;
        elseif strcmp(datacell{10, fileindex}, 'vest') == 1
            vest_counter = vest_counter + 1;
        elseif strcmp(datacell{10, fileindex}, 'face') == 1
            face_counter = face_counter + 1;
        end
    end
end
% Create matrices for each category in which to place data
energies = data(:, [1,2]);
air_matrix = zeros(8192, air_counter + 2);
air_matrix(:, [1,2]) = energies;
feet_matrix = zeros(8192, feet_counter + 2);
feet_matrix(:,[1,2]) = energies;
gloves_matrix = zeros(8192, gloves_counter + 2);
gloves_matrix(:, [1,2]) = energies;
pants_matrix = zeros(8192, pants_counter + 2);
pants_matrix(:, [1,2]) = energies;
vest_matrix = zeros(8192, vest_counter + 2);
vest_matrix(:, [1,2]) = energies;
face_matrix = zeros(8192, face_counter + 2);
face_matrix(:, [1,2]) = energies;

% For each set of data, input the experimental values into columns of data
% matrices
air_count = 0;
feet_count = 0;
gloves_count = 0;
pants_count = 0;
vest_count = 0;
face_count = 0;

for fileindex = 1:length(files)
    is_control = datacell{5, fileindex};
    if is_control ~= 1
        if strcmp(datacell{10, fileindex}, 'air') == 1
            air_count = air_count + 1;
            data = datacell{9, fileindex};
            data_column = data(:,3);
            air_matrix(:, air_count + 2) = data_column;
        elseif strcmp(datacell{10,fileindex}, 'feet') == 1
            data = datacell{9, fileindex};
            data_column = data(:,3);
            feet_count = feet_count + 1;
            feet_matrix(:, feet_count + 2) = data_column;
        elseif strcmp(datacell{10,fileindex}, 'gloves') == 1
            data = datacell{9, fileindex};
            data_column = data(:,3);
            gloves_count = gloves_count + 1;
            gloves_matrix(:, gloves_count + 2) = data_column;
        elseif strcmp(datacell{10, fileindex}, 'pants') == 1
            data = datacell{9, fileindex};
            data_column = data(:,3);
            pants_count = pants_count + 1;
            pants_matrix(:, pants_count + 2) = data_column;
        elseif strcmp(datacell{10, fileindex}, 'vest') == 1
            data = datacell{9, fileindex};
            data_column = data(:,3);
            vest_count = vest_count + 1;
            vest_matrix(:, vest_count + 2) = data_column;
        elseif strcmp(datacell{10, fileindex}, 'face') == 1
            data = datacell{9, fileindex};
            data_column = data(:,3);
            face_count = face_count + 1;
            face_matrix(:, face_count + 2) = data_column;
        end
    end
end

% Create matrices for the average and standard deviations of each set

air_average_matrix = zeros(8192, 4);
air_average_matrix(:, [1 2]) = energies;
values = air_matrix(:, 3:air_count);
air_average_matrix(:, 3) = mean(values, 2);
air_average_matrix(:, 4) = std(values, 0, 2);

feet_average_matrix = zeros(8192, 4);
feet_average_matrix(:, [1 2]) = energies;
values = feet_matrix(:, 3:gloves_count);
feet_average_matrix(:, 3) = mean(values, 2);
feet_average_matrix(:, 4) = std(values, 0, 2);

gloves_average_matrix = zeros(8192, 4);
gloves_average_matrix(:, [1 2]) = energies;
values = gloves_matrix(:, 3:gloves_count);
gloves_average_matrix(:, 3) = mean(values, 2);
gloves_average_matrix(:, 4) = std(values, 0, 2);

pants_average_matrix = zeros(8192, 4);
pants_average_matrix(:, [1 2]) = energies;
values = pants_matrix(:, 3:pants_count);
pants_average_matrix(:, 3) = mean(values, 2);
pants_average_matrix(:, 4) = std(values, 0, 2);

vest_average_matrix = zeros(8192, 4);
vest_average_matrix(:, [1 2]) = energies;
values = vest_matrix(:, 3:vest_count);
vest_average_matrix(:, 3) = mean(values, 2);
vest_average_matrix(:, 4) = std(values, 0, 2);

face_average_matrix = zeros(8192, 4);
face_average_matrix(:, [1 2]) = energies;
values = face_matrix(:, 3:face_count);
face_average_matrix(:, 3) = mean(values, 2);
face_average_matrix(:, 4) = std(values, 0, 2);

% Plot average Bq/g above background with standard deviations for each
% category

figure
semilogy(air_average_matrix(:,2), (air_average_matrix(:,3) - air_average_matrix(:, 4)), 'color', 'c');
hold on
semilogy(air_average_matrix(:,2), (air_average_matrix(:,3) + air_average_matrix(:, 4)), 'color', 'c');
semilogy(air_average_matrix(:, 2), air_average_matrix(:,3), 'color', 'b');
ylim([.001 1000]);
ylabel('Average Bq/g');
xlabel('Energy');
title('Average Spectrum for Air Samples');

figure
semilogy(gloves_average_matrix(:,2), (gloves_average_matrix(:,3) - gloves_average_matrix(:, 4)), 'color', 'c');
hold on
semilogy(gloves_average_matrix(:,2), (gloves_average_matrix(:,3) + gloves_average_matrix(:, 4)), 'color', 'c');
semilogy(gloves_average_matrix(:, 2), gloves_average_matrix(:,3), 'color', 'b');
ylim([.001 1000]);
ylabel('Average Bq/g');
xlabel('Energy');
title('Average Spectrum for Hand Covering Samples');

figure
semilogy(feet_average_matrix(:,2), (feet_average_matrix(:,3) - feet_average_matrix(:, 4)), 'color', 'c');
hold on
semilogy(feet_average_matrix(:,2), (feet_average_matrix(:,3) + feet_average_matrix(:, 4)), 'color', 'c');
semilogy(feet_average_matrix(:, 2), feet_average_matrix(:,3), 'color', 'b');
ylim([.001 1000]);
ylabel('Average Bq/g');
xlabel('Energy');
title('Average Spectrum for Foot Covering Samples');

figure
semilogy(pants_average_matrix(:,2), (pants_average_matrix(:,3) - pants_average_matrix(:, 4)), 'color', 'c');
hold on
semilogy(pants_average_matrix(:,2), (pants_average_matrix(:,3) + pants_average_matrix(:, 4)), 'color', 'c');
semilogy(pants_average_matrix(:, 2), pants_average_matrix(:,3), 'color', 'b');
ylim([.001 1000]);
ylabel('Average Bq/g');
xlabel('Energy');
title('Average Spectrum for Leg Covering Samples');

figure
semilogy(vest_average_matrix(:,2), (vest_average_matrix(:,3) - vest_average_matrix(:, 4)), 'color', 'c');
hold on
semilogy(vest_average_matrix(:,2), (vest_average_matrix(:,3) + vest_average_matrix(:, 4)), 'color', 'c');
semilogy(vest_average_matrix(:, 2), vest_average_matrix(:,3), 'color', 'b');
ylim([.001 1000]);
ylabel('Average Bq/g');
xlabel('Energy');
title('Average Spectrum for Chest Covering Samples');

figure
semilogy(face_average_matrix(:,2), (face_average_matrix(:,3) - face_average_matrix(:, 4)), 'color', 'c');
hold on
semilogy(face_average_matrix(:,2), (face_average_matrix(:,3) + face_average_matrix(:, 4)), 'color', 'c');
semilogy(face_average_matrix(:, 2), face_average_matrix(:,3), 'color', 'b');
ylim([.001 1000]);
ylabel('Average Bq/g');
xlabel('Energy');
title('Average Spectrum for Face Covering Samples');

%% Determine the area under the curve for the Cs-137 peaks in each set of samples and the total standard deviation

Cs137_channel_range = [3518 3533];
feet_Cs137_activity = sum(feet_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),3));
feet_Cs137_stdev = sqrt(sum((feet_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),4).^2)));
gloves_Cs137_activity = sum(gloves_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),3));
gloves_Cs137_stdev = sqrt(sum((gloves_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),4).^2)));
air_Cs137_activity = sum(air_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),3));
air_Cs137_stdev = sqrt((sum(air_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),4).^2)));
pants_Cs137_activity = sum(pants_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),3));
pants_Cs137_stdev = sqrt(sum((pants_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),4)).^2));
vest_Cs137_activity = sum(vest_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),3));
vest_Cs137_stdev = sqrt(sum(vest_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),4).^2));
face_Cs137_activity = sum(face_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),3));
face_Cs137_stdev = sqrt(sum(face_average_matrix(Cs137_channel_range(1):Cs137_channel_range(2),4).^2));

% Plot a bar chart of this data

Cs137_activities = [air_Cs137_activity, feet_Cs137_activity, gloves_Cs137_activity, pants_Cs137_activity, face_Cs137_activity];
Cs137_stdev = [air_Cs137_stdev, feet_Cs137_stdev, gloves_Cs137_activity, pants_Cs137_activity, face_Cs137_activity];
figure
p = bar(Cs137_activities, 'b');
hold on
errorbar(Cs137_activities, Cs137_stdev, '.', 'color', 'k');
title('Cesium-137 Activity by Category');