% MATLAB Script for Signal Analysis and Particle Identification
% Load your signal data from the oscilloscope
% Assuming 'signals' is a matrix of size 400 x 10000

% Example data loading (replace this with actual data loading)
% signals = load('your_data_file.mat'); % Assuming a MAT-file

% For demonstration, generate random data
% Comment this out when using actual data
signals = Channel_1; % Replace with actual data

% Parameters
[num_points, num_signals] = size(signals);
thresholds = 0.1 * max(signals); % Calculate thresholds for each signal
time_per_point = 50 / num_points; % Time per point in nanoseconds

% Preallocate arrays for performance
amplitudes = zeros(1, num_signals);
areas = zeros(1, num_signals);
widths_10_percent_ns = zeros(1, num_signals);
form_factors = zeros(1, num_signals);

% Analyze each signal
for i = 1:num_signals
    % Extract the signal
    signal = signals(:, i);
    
    % Calculate amplitude as the peak voltage value
    amplitudes(i) = max(signal);
    
    % Calculate area under the curve
    areas(i) = trapz(signal);
    
    % Find width at 10% of maximum
    threshold = 0.1 * amplitudes(i); % 10% of the amplitude
    above_ten_percent = find(signal > threshold);
    width_points = length(above_ten_percent);
    widths_10_percent_ns(i) = width_points * time_per_point;
    
    % Calculate form factor
    calculated_area = amplitudes(i) * widths_10_percent_ns(i);
    form_factors(i) = calculated_area / areas(i);
end

% Classify signals based on form factor
particle_types = cell(1, num_signals);
for i = 1:num_signals
    if areas(i) > 2.3
        particle_types{i} = 'Rectangular Signal (Heavy Ionizing Particle)';
    elseif form_factors(i) >= 1.4 && form_factors(i) <= 2.0
        particle_types{i} = 'Triangular Signal (Traversing Particle)';
    else
        particle_types{i} = 'Unknown Signal Shape';
    end
end

% Store the parameters in a table named signal_iden
signal_iden = table((1:num_signals)', amplitudes', areas', widths_10_percent_ns', form_factors', particle_types', ...
    'VariableNames', {'SignalNumber', 'Amplitude', 'Area', 'Width_10_Percent_ns', 'FormFactor', 'ParticleType'});

% Display the first few rows of the table
disp(signal_iden(1:1000, :));

% Plotting the spectrum of wb versus area with logarithmic color scale
figure;
scatter(areas, widths_10_percent_ns, 20, log10(1:num_signals), 'filled'); % Ensuring the index matches
colormap(jet);
colorbar;
xlabel('Area (pVs)');
ylabel('Width at 10% of max (ns)');
set(gca, 'ColorScale', 'log');
grid on;

% Highlight specific regions
hold on;
% Mark region for fast neutrons
fast_neutrons_area = areas > 300 & widths_10_percent_ns > 5 & widths_10_percent_ns < 7;
scatter(areas(fast_neutrons_area), widths_10_percent_ns(fast_neutrons_area), 40, 'r', 'o');

% Mark region for 6Li(n,α)3H reaction products
reaction_products_area = areas <= 300 & widths_10_percent_ns <= 5;
scatter(areas(reaction_products_area), widths_10_percent_ns(reaction_products_area), 40, 'b', 'o');

legend('Logarithmic Count Rate', 'Fast Neutrons', '6Li(n,α)3H Reaction Products');