clc;clear;close all
%% 参数设置
num_events = 1e5; % 模拟事件数
E = 662e3; % γ射线能量（eV）
Y = 38000; % 光产额（photon/MeV，NaI(Tl)）
eta = 0.2; % 量子效率
G = 1e6; % PMT增益
F = 1.5; % 法诺因子
C = 100e-12; % 电容（F）
e = 1.6e-19; % 电子电荷量（C）
sigma_noise_e = 5000; % 电子学噪声（电子数RMS）%% 各环节模拟
% 1. 光子生成（泊松涨落）
mean_photon = Y * E / 1e6; % 转换为MeV计算
N_photon = poissrnd(mean_photon, [num_events, 1]);% 2. 光电子转换（二项涨落）
N_pe = arrayfun(@(n) binornd(n, eta), N_photon);% 3. PMT倍增（法诺涨落，高斯近似）
Q_pmt = normrnd(N_pe*G*e, sqrt(N_pe*G^2*F)*e);% 4. 电子学噪声（高斯噪声）
Q_noise = normrnd(0, sigma_noise_e*e, [num_events, 1]);
Q_total = Q_pmt + Q_noise;% 5. 电压转换
V = Q_total / C;%% 结果可视化
% 直方图与高斯拟合
figure
h = histogram(V, 100, 'Normalization', 'pdf');
hold on% 计算理论高斯参数
mu_V = mean(V);
sigma_V = std(V);
x = linspace(min(V), max(V), 500);
y = normpdf(x, mu_V, sigma_V);
plot(x, y, 'r', 'LineWidth', 2)% 图像标注
xlabel('Pulse Amplitude (V)')
ylabel('Probability Density')
title(['Simulated Spectrum @ ', num2str(E/1e3), ' keV'])
legend('Simulation', 'Gaussian Fit')
grid on% 显示能量分辨率
fwhm = 2.355 * sigma_V;
resolution = fwhm/mu_V;
text(0.05, 0.9, sprintf('Resolution: %.1f%%', resolution*100),...'Units', 'normalized', 'FontSize', 12)%% 关键参数输出
disp('===== 关键参数 =====')
fprintf('理论脉冲幅度: %.2e V\n', Y*eta*G*e*E/(1e6*C))
fprintf('实测脉冲幅度: %.2e V\n', mu_V)
fprintf('实测能量分辨率: %.1f%%\n', resolution*100)
运行结果如下
考虑康普顿沉积的话
clc; clear; close all;%% 参数设置
num_events = 1e5; % 模拟事件数
E_MeV = 0.662; % γ射线能量（MeV）
E = E_MeV * 1e6; % 转换为eV
Y = 38000; % 光产额（photon/MeV，NaI(Tl)）
eta = 0.2; % 量子效率
G = 1e6; % PMT增益
F = 1.5; % 法诺因子
C = 100e-12; % 电容（F）
e = 1.6e-19; % 电子电荷量（C）
sigma_noise_e = 5000; % 电子学噪声（电子数RMS）% 相互作用概率（示例值，需根据实际能量计算）
mu_photo = 0.6; % 光电效应占比
mu_compton = 0.3; % 康普顿散射占比
mu_pair = 0.0; % 电子对效应（E < 1.022 MeV时为0）
total_mu = mu_photo + mu_compton + mu_pair;
prob = [mu_photo, mu_compton, mu_pair] / total_mu;%% 蒙特卡洛模拟能量沉积
deposited_energy = zeros(num_events, 1);
for i = 1:num_eventsE_remaining = E;while E_remaining > 0interaction = string(randsample({'photo', 'compton', 'pair'}, 1, true, prob));switch interactioncase 'photo'% 光电效应 + X射线逃逸if rand() < 0.2E_X = 30e3; % 碘的X射线能量deposited_energy(i) = deposited_energy(i) + (E_remaining - E_X);elsedeposited_energy(i) = deposited_energy(i) + E_remaining;endE_remaining = 0;case 'compton'% 康普顿散射cos_theta = 2*rand() - 1;E_deposit = E_remaining * (1 - 1/(1 + (E_remaining/511e3)*(1 - cos_theta)));deposited_energy(i) = deposited_energy(i) + E_deposit;E_remaining = E_remaining - E_deposit;if rand() < 0.5 % 50%概率逃逸E_remaining = 0;endendend
end%% 探测过程模拟
% 1. 光子生成（基于实际沉积能量）
mean_photon = Y * deposited_energy / 1e6; % deposited_energy单位为eV
N_photon = poissrnd(mean_photon);% 2-5. 后续步骤与原代码相同
N_pe = arrayfun(@(n) binornd(n, eta), N_photon);
Q_pmt = normrnd(N_pe*G*e, sqrt(N_pe*G^2*F)*e);
Q_noise = normrnd(0, sigma_noise_e*e, [num_events, 1]);
Q_total = Q_pmt + Q_noise;
V = Q_total / C;%% 结果可视化
figure;
histogram(V, 1024, 'BinLimits', [0, 1.2*max(V)], 'EdgeColor', 'none');
xlabel('Pulse Amplitude (V)');
ylabel('Counts');
title(sprintf('Simulated Spectrum with Compton Edge @ %.1f MeV', E_MeV));
grid on;% 标注特征峰
hold on;
% 光电峰
line([Y*eta*G*e*E/(1e6*C), Y*eta*G*e*E/(1e6*C)], ylim, 'Color', 'r', 'LineStyle', '--');
% 康普顿边缘
E_compton_edge = E_MeV * (1 - 1/(1 + 2*E_MeV/0.511)); % 康普顿最大能量沉积
V_compton = Y*eta*G*e*(E_compton_edge*1e6)/(1e6*C);
line([V_compton, V_compton], ylim, 'Color', 'g', 'LineStyle', '--');
legend('Simulated', 'Photopeak', 'Compton Edge');
模拟Co、Ba、Cs三种放射性核素的能谱，其活度比例为1:2:3，并包含各核素的特征γ射线能量：
clc; clear; close all;%% 参数设置
num_events_total = 6e4; % 总事件数（1+2+3=6份）
e = 1.6e-19; % 电子电荷量（C）% 定义探测器参数（与闪烁体相关）
Y = 38000; % 光产额（photon/MeV，NaI(Tl)）
eta = 0.2; % 量子效率
G = 1e6; % PMT增益
F = 1.5; % 法诺因子
C = 100e-12; % 电容（F）
sigma_noise_e = 5000; % 电子学噪声（电子数RMS）% 定义放射源特性（核素名称、能量[MeV]、分支比、活度权重）
sources = struct(...'name', {'Co60', 'Ba133', 'Cs137'},...'energy', {[1.173, 1.332], 0.356, 0.662},... % γ射线能量（MeV）'branch', {[0.5, 0.5], 1.0, 1.0},... % 分支比（概率和=1）'activity', {1, 2, 3}... % 活度比例
);%% 生成各核素事件
all_E = []; % 存储所有事件的能量% 遍历每个放射源
for i = 1:length(sources)% 计算该核素的事件数 (按活度比例分配)num_events = round(num_events_total * sources(i).activity / sum([sources.activity]));% 生成该核素的γ射线能量抽样if length(sources(i).energy) > 1% 多能峰：按分支比随机选择能量E_choices = sources(i).energy;probs = sources(i).branch;% 累积概率法cum_probs = cumsum(probs) / sum(probs); % 确保分支比总和为1rand_vals = rand(num_events, 1); % 生成随机数selected_E = zeros(num_events, 1);for k = 1:length(E_choices)if k == 1selected_E(rand_vals <= cum_probs(k)) = E_choices(k);elseselected_E((rand_vals > cum_probs(k-1)) & (rand_vals <= cum_probs(k))) = E_choices(k);endendelse% 单能峰：直接重复能量selected_E = repmat(sources(i).energy, num_events, 1);end% 转换为eV并存储all_E = [all_E; selected_E * 1e6]; % 转换为eV
end% 定义本底参数
background_intensity = 0.1; % 本底强度（相对于总事件数的比例）
num_background_events = round(num_events_total * background_intensity); % 本底事件数% 生成本底能量分布（假设为均匀分布，范围覆盖所有核素能量）
E_min = min([sources.energy]); % 最小能量 (MeV)
E_max = max([sources.energy]); % 最大能量 (MeV)
background_E = E_min + (E_max - E_min) * rand(num_background_events, 1); % 本底能量 (MeV)
background_E = background_E * 1e6; % 转换为 eV% 合并本底事件与核素事件
all_E = [all_E; background_E]; % 将本底能量添加到总事件中
%% 模拟探测过程（向量化计算）
% 1. 光子生成（泊松涨落）
N_photon = poissrnd(Y * all_E / 1e6); % 输入能量已为eV，/1e6转换为MeV% % 2. 光电子转换（二项涨落）
% N_pe = binornd(N_photon, eta);% 使用正态分布近似二项分布
mu = N_photon * eta; % 均值
sigma = sqrt(N_photon * eta * (1 - eta)); % 标准差
N_pe = normrnd(mu, sigma);% 3. PMT倍增（法诺涨落，高斯近似）
Q_pmt = normrnd(N_pe*G*e, sqrt(N_pe*G^2*F)*e);% 4. 电子学噪声（高斯噪声）
Q_noise = normrnd(0, sigma_noise_e*e, size(Q_pmt));
Q_total = Q_pmt + Q_noise;% 5. 电压转换
V = Q_total / C;%% 模拟本底探测过程
% 1. 光子生成（泊松涨落）
N_photon_background = poissrnd(Y * background_E / 1e6); % 计算本底光子数% 使用正态分布近似二项分布
mu_background = N_photon_background * eta; % 均值
sigma_background = sqrt(N_photon_background * eta * (1 - eta)); % 标准差
N_pe_background = normrnd(mu_background, sigma_background);% 2. PMT倍增（法诺涨落，高斯近似）
Q_pmt_background = normrnd(N_pe_background * G * e, sqrt(N_pe_background * G^2 * F) * e);% 3. 电子学噪声（高斯噪声）
Q_noise_background = normrnd(0, sigma_noise_e * e, size(Q_pmt_background));
Q_total_background = Q_pmt_background + Q_noise_background;% 4. 电压转换
V_background = Q_total_background / C;% 合并本底信号与核素信号
V = [V; V_background];%% 结果可视化
figure;
hold on;% 绘制总谱线
histogram(V, 2048, 'DisplayName', 'Total Spectrum', 'EdgeColor', 'none');% 绘制本底谱线
histogram(V_background, 2048, 'FaceColor', [0.8, 0.8, 0.8], 'FaceAlpha', 0.5,...'DisplayName', 'Background Spectrum');% % 单独绘制各核素谱线（透明叠加）
colors = lines(length(sources));
% for i = 1:length(sources)
% % 提取该核素对应的事件范围
% start_idx = sum([sources(1:i-1).activity])/sum([sources.activity]) * num_events_total + 1;
% end_idx = sum([sources(1:i).activity])/sum([sources.activity]) * num_events_total;
% range = round(start_idx):round(end_idx);
%
% histogram(V(range), 500,...
% 'FaceColor', colors(i,:), 'FaceAlpha', 0.5,...
% 'DisplayName', sources(i).name);
% end% 理论峰位标注
for i = 1:length(sources)for j = 1:length(sources(i).energy)E_MeV = sources(i).energy(j);mu_V = Y * eta * G * e * (E_MeV*1e6) / (1e6 * C);line([mu_V, mu_V], ylim, 'Color', colors(i,:),...'LineStyle', '--', 'LineWidth', 1.5,...'DisplayName', sprintf('%s %.3f MeV', sources(i).name, E_MeV));end
end% 图形修饰
xlabel('Pulse Amplitude (V)');
ylabel('Counts');
title('Simulated Gamma Spectrum (Co:Ba:Cs = 1:2:3)');
legend('show', 'Location', 'northwest');
grid on;
% set(gca, 'YScale', 'log'); % 对数坐标更好展示低计数区域%% 关键参数输出
disp('===== 模拟参数 =====');
fprintf('总事件数: %d\n', num_events_total);
fprintf('本底事件数: %d\n', num_background_events);
disp('各核素理论峰位:');
for i = 1:length(sources)for j = 1:length(sources(i).energy)E_MeV = sources(i).energy(j);mu_V = Y * eta * G * e * (E_MeV * 1e6) / (1e6 * C);fprintf('%s %.3f MeV → %.2e V\n',...sources(i).name, E_MeV, mu_V);end
end
