% MATLAB Script for Signal Analysis and Particle Identification
% Load your signal data from the oscilloscope
% Assuming 'Channel_1' is a matrix of size [num_points x num_signals]

% Example data loading (replace this with actual data loading)
% load('your_data_file.mat'); % Uncomment and replace with your actual file

% For demonstration, use Channel_1 from workspace
signals = Channel_1; % Replace with actual data if necessary

% Parameters
[num_points, num_signals] = size(signals);
time_per_point = 40 / 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);
rise_times = zeros(1, num_signals);
fall_times = zeros(1, num_signals);
delta_T = 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 using trapezoidal rule
    areas(i) = trapz(signal);
    
    % Find width at 10% of maximum
    threshold = 0.1 * amplitudes(i); % 10% of the amplitude
    FWHM = 0.5 * amplitudes(i);
    above_ten_percent = find(signal > threshold);
    if isempty(above_ten_percent)
        width_points = 0;
    else
        width_points = length(above_ten_percent);
    end
    widths_10_percent_ns(i) = width_points * time_per_point;
    
    % Calculate form factor
    calculated_area = amplitudes(i) * widths_10_percent_ns(i);
    if areas(i) == 0
        form_factors(i) = Inf; % Avoid division by zero
    else
        form_factors(i) = calculated_area / areas(i);
    end
    
    % Calculate Rise Time (10% to 90%)
    lower_threshold = 0.1 * amplitudes(i);
    upper_threshold = 0.9 * amplitudes(i);
    
    % Find indices for rise time
    rise_start_idx = find(signal >= lower_threshold, 1, 'first');
    rise_end_idx = find(signal >= upper_threshold, 1, 'first');
    
    if ~isempty(rise_start_idx) && ~isempty(rise_end_idx) && rise_end_idx > rise_start_idx
        rise_times(i) = (rise_end_idx - rise_start_idx) * time_per_point;
    else
        rise_times(i) = NaN; % Assign NaN if not found
    end
    
    % Calculate Fall Time (90% to 10%)
    fall_start_idx = find(signal <= upper_threshold, 1, 'last');
    fall_end_idx = find(signal <= lower_threshold, 1, 'last');
    
    if ~isempty(fall_start_idx) && ~isempty(fall_end_idx) && fall_end_idx < fall_start_idx
        fall_times(i) = (fall_start_idx - fall_end_idx) * time_per_point;
    else
        fall_times(i) = NaN; % Assign NaN if not found
    end
    
    % Calculate Delta T (Rise Time - Fall Time)
    if ~isnan(rise_times(i)) && ~isnan(fall_times(i))
        delta_T(i) = rise_times(i) - fall_times(i);
    else
        delta_T(i) = NaN; % Assign NaN if either rise or fall time is undefined
    end
end

% Classify signals based on form factor
particle_types = cell(1, num_signals);
for i = 1:num_signals
    if form_factors(i) < 1.4 && delta_T(i) >= 10
        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', rise_times', fall_times', delta_T', particle_types', ...
    'VariableNames', {'SignalNumber', 'Amplitude', 'Area', ...
    'Width_10_Percent_ns', 'FormFactor', 'RiseTime_ns', ...
    'FallTime_ns', 'DeltaT_ns', 'ParticleType'});

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

% --- Filter and Plot the Rectangular Signals ---
% Criteria:
% 1. Particle Type: Rectangular Signal
% 2. Width at 10% > 10 ns
% 3. |Delta T| < Threshold (e.g., 5 ns)

% Define threshold for Delta T
delta_T_threshold = 5; % nanoseconds

% Find indices of signals that meet all criteria
valid_rectangular_signals = find(strcmp(particle_types, 'Rectangular Signal (Heavy Ionizing Particle)') & ...
    widths_10_percent_ns > 10 & abs(delta_T) < delta_T_threshold);

% Check if there are any valid rectangular signals
if isempty(valid_rectangular_signals)
    disp('No valid Rectangular Signals found with width > 10 ns and |Delta T| < 5 ns.');
else
    % Plot the valid rectangular signals (waveforms)
    figure;
    hold on;
    num_plots = min(length(valid_rectangular_signals), 50); % Limit to first 50 for clarity
    for i = 1:num_plots
        signal_index = valid_rectangular_signals(i);
        signal = signals(:, signal_index);
        
        % Plot the signal (rectangular pulse)
        plot((1:num_points) * time_per_point, signal, 'DisplayName', ['Signal ' num2str(signal_index)]);
    end
    hold off;
    
    xlabel('Time (ns)');
    ylabel('Amplitude (V)');
    title(['Rectangular Signal Waveforms (Heavy Ionizing Particles) - Width > 10 ns & |ΔT| < ' num2str(delta_T_threshold) ' ns']);
    legend('show', 'Location', 'best');
    grid on;
end

% --- Plotting the Spectrum of Width vs Area ---
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');
hold off;

% --- Plot Histogram of Delta T to Analyze Signal Symmetry ---
figure;
histogram(delta_T, 'BinWidth', 1, 'FaceColor', 'g');
xlabel('Rise Time - Fall Time (ns)');
ylabel('Count');
title('Histogram of Delta T (Rise Time - Fall Time)');
grid on;

% --- Optional: Plotting Histogram of Rise Times and Fall Times ---
figure;
subplot(2,1,1);
histogram(rise_times, 'BinWidth', 1, 'FaceColor', 'b');
xlabel('Rise Time (ns)');
ylabel('Count');
title('Histogram of Rise Times');
grid on;

subplot(2,1,2);
histogram(fall_times, 'BinWidth', 1, 'FaceColor', 'm');
xlabel('Fall Time (ns)');
ylabel('Count');
title('Histogram of Fall Times');
grid on;