272 lines
8.7 KiB
C++
272 lines
8.7 KiB
C++
#include <iostream>
|
||
#include <QDebug>
|
||
#include "CoincidenceSpectrumProcess.h"
|
||
|
||
using namespace arma;
|
||
|
||
|
||
namespace CoincidenceSpectrum {
|
||
|
||
mat EventFilter2Order(const mat& raw_spec_data, double coincidence_event_time_win)
|
||
{
|
||
// 原始谱数据排序
|
||
auto raw_spec_data_size = size(raw_spec_data);
|
||
mat temp_data(raw_spec_data_size.n_rows, 3, fill::zeros);
|
||
temp_data.col(0) = raw_spec_data.col(3);
|
||
temp_data.col(1) = raw_spec_data.col(2);
|
||
temp_data.col(2) = 4 * raw_spec_data.col(0) + raw_spec_data.col(1) + 1;
|
||
uvec indices = sort_index(temp_data.col(0), "ascend");
|
||
mat raw_spec_sorted_data = temp_data.rows(indices);
|
||
|
||
// 整理次级谱数据
|
||
temp_data.clear();
|
||
SizeMat raw_spec_sorted_data_size = size(raw_spec_sorted_data);
|
||
uword rows_count = raw_spec_sorted_data_size.n_rows;
|
||
uword cols_count = raw_spec_sorted_data_size.n_cols;
|
||
uword last_row_num = rows_count - 1;
|
||
uword last_col_num = cols_count - 1;
|
||
temp_data = raw_spec_sorted_data.submat(1, 0, last_row_num, last_col_num);
|
||
temp_data.insert_rows(last_row_num, raw_spec_sorted_data.row(last_row_num));
|
||
|
||
// 计算得到时间差谱数据
|
||
mat coincidence_spec_data(rows_count, 6);
|
||
coincidence_spec_data.col(0) = temp_data.col(0) - raw_spec_sorted_data.col(0);
|
||
coincidence_spec_data.col(1) = raw_spec_sorted_data.col(1);
|
||
coincidence_spec_data.col(2) = raw_spec_sorted_data.col(2);
|
||
coincidence_spec_data.col(3) = temp_data.col(1);
|
||
coincidence_spec_data.col(4) = temp_data.col(2);
|
||
coincidence_spec_data.col(5) = temp_data.col(0);
|
||
raw_spec_sorted_data.clear();
|
||
temp_data.clear();
|
||
|
||
// 筛选出符合谱数据
|
||
indices = sort_index(coincidence_spec_data.col(0), "ascend");
|
||
coincidence_spec_data = coincidence_spec_data.rows(indices);
|
||
auto time_difference_col = coincidence_spec_data.col(0);
|
||
indices = find(time_difference_col <= coincidence_event_time_win);
|
||
coincidence_spec_data = coincidence_spec_data.rows(indices);
|
||
|
||
return coincidence_spec_data;
|
||
}
|
||
|
||
// 计算二维直方图
|
||
void Hist3(const mat& data, int nx, int ny, mat& counts, vec& x_edges, vec& y_edges)
|
||
{
|
||
// 检查输入数据是否为N×2矩阵
|
||
if (data.n_cols != 2) {
|
||
cerr << "输入数据必须是N×2矩阵" << endl;
|
||
return;
|
||
}
|
||
uword n = data.n_rows;
|
||
if (n == 0) {
|
||
counts = zeros<mat>(nx, ny);
|
||
x_edges = zeros<vec>(nx + 1);
|
||
y_edges = zeros<vec>(ny + 1);
|
||
return;
|
||
}
|
||
|
||
// 提取x和y数据
|
||
vec x = data.col(0); // 对应MATLAB的第2列
|
||
vec y = data.col(1); // 对应MATLAB的第3列
|
||
|
||
// 计算x和y的分箱边界（等宽分箱）
|
||
double x_min = x.min();
|
||
double x_max = x.max();
|
||
double y_min = y.min();
|
||
double y_max = y.max();
|
||
|
||
// 处理x或y所有值相同的特殊情况（避免除以0）
|
||
if (x_min == x_max) {
|
||
x_edges = linspace<vec>(x_min - 1e-6, x_max + 1e-6, nx + 1); // 微小范围
|
||
} else {
|
||
x_edges = linspace<vec>(x_min, x_max, nx + 1); // nx+1个边界，nx个分箱
|
||
}
|
||
if (y_min == y_max) {
|
||
y_edges = linspace<vec>(y_min - 1e-6, y_max + 1e-6, ny + 1);
|
||
} else {
|
||
y_edges = linspace<vec>(y_min, y_max, ny + 1);
|
||
}
|
||
|
||
// 初始化频数矩阵（nx行，ny列）
|
||
counts = zeros<mat>(nx, ny);
|
||
|
||
// 遍历每个数据点，统计频数
|
||
for (int i = 0; i < n; ++i) {
|
||
double xi = x(i);
|
||
double yi = y(i);
|
||
// 找到xi所在的x分箱索引（0-based）
|
||
int bin_x = -1;
|
||
for (int j = 0; j < nx; ++j) {
|
||
if (xi >= x_edges(j) && xi < x_edges(j + 1)) {
|
||
bin_x = j;
|
||
break;
|
||
}
|
||
}
|
||
// 处理刚好等于最大值的点（放入最后一个分箱）
|
||
if (bin_x == -1 && xi == x_edges(nx)) {
|
||
bin_x = nx - 1;
|
||
}
|
||
// 找到yi所在的y分箱索引（0-based）
|
||
int bin_y = -1;
|
||
for (int j = 0; j < ny; ++j) {
|
||
if (yi >= y_edges(j) && yi < y_edges(j + 1)) {
|
||
bin_y = j;
|
||
break;
|
||
}
|
||
}
|
||
if (bin_y == -1 && yi == y_edges(ny)) {
|
||
bin_y = ny - 1;
|
||
}
|
||
// 累加频数（忽略超出范围的点）
|
||
if (bin_x >= 0 && bin_x < nx && bin_y >= 0 && bin_y < ny) {
|
||
counts(bin_x, bin_y)++;
|
||
}
|
||
}
|
||
}
|
||
|
||
namespace F2t9Order {
|
||
|
||
// 读取CSV数据文件
|
||
std::vector<SpectrumData> ReadCsv(const std::string& filename)
|
||
{
|
||
std::vector<SpectrumData> data;
|
||
std::ifstream file(filename);
|
||
|
||
if (!file.is_open()) {
|
||
std::cerr << "错误：无法打开文件 " << filename << std::endl;
|
||
return data;
|
||
}
|
||
|
||
std::string line;
|
||
// 跳过表头
|
||
std::getline(file, line);
|
||
|
||
while (std::getline(file, line)) {
|
||
std::stringstream ss(line);
|
||
std::string token;
|
||
std::vector<std::string> tokens;
|
||
|
||
while (std::getline(ss, token, ',')) {
|
||
tokens.push_back(token);
|
||
}
|
||
|
||
if (tokens.size() == 4) {
|
||
SpectrumData entry;
|
||
entry.board_id = std::stoi(tokens[0]);
|
||
entry.channel_id = std::stoi(tokens[1]);
|
||
entry.energy = std::stod(tokens[2]);
|
||
entry.timestamp = std::stoull(tokens[3]);
|
||
data.push_back(entry);
|
||
}
|
||
}
|
||
|
||
file.close();
|
||
std::cout << "成功读取 " << data.size() << " 条有效数据" << std::endl;
|
||
return data;
|
||
}
|
||
|
||
// 按时间戳排序数据
|
||
void SortDataByTimestamp(std::vector<SpectrumData>& data)
|
||
{
|
||
std::sort(data.begin(), data.end(),
|
||
[](const SpectrumData& a, const SpectrumData& b) {
|
||
return a.timestamp < b.timestamp;
|
||
});
|
||
}
|
||
|
||
// 2-9次能谱符合处理
|
||
bool ProcessCoincidence(
|
||
std::vector<SpectrumData>& data,
|
||
std::vector<CoincidenceEvent>& events,
|
||
unsigned int time_window,
|
||
int min_order,
|
||
int max_order
|
||
)
|
||
{
|
||
int n = data.size();
|
||
if (n < min_order) {
|
||
return false;
|
||
}
|
||
// 使用滑动窗口寻找符合事件
|
||
for (int i = 0; i < n; ++i) {
|
||
std::vector<SpectrumData> current_coincidence;
|
||
current_coincidence.push_back(data[i]);
|
||
// 寻找时间窗口内的其他事件
|
||
for (int j = i + 1; j < n; ++j) {
|
||
unsigned long long time_diff = data[j].timestamp - data[i].timestamp;
|
||
if (time_diff > time_window) {
|
||
break; // 由于已排序，后续事件时间差会更大
|
||
}
|
||
current_coincidence.push_back(data[j]);
|
||
// 达到最大符合次数时停止当前窗口
|
||
if (current_coincidence.size() >= max_order) {
|
||
break;
|
||
}
|
||
}
|
||
// 保存符合条件的事件
|
||
int coincidence_order = current_coincidence.size();
|
||
if (coincidence_order >= min_order && coincidence_order <= max_order) {
|
||
CoincidenceEvent result;
|
||
result.coincidence_order = coincidence_order;
|
||
result.events = current_coincidence;
|
||
result.time_window = time_window;
|
||
events.push_back(result);
|
||
}
|
||
}
|
||
return true;
|
||
}
|
||
|
||
// 统计符合事件结果
|
||
arma::mat AnalyzeCoincidenceResults(const std::vector<CoincidenceEvent>& results, int max_order)
|
||
{
|
||
// 统计各粒子数符合事件的数量和能量信息
|
||
arma::mat stats(max_order + 1, 2, arma::fill::zeros);
|
||
// 列：符合粒子数、事件数量
|
||
for (const auto& event : results) {
|
||
int order = event.coincidence_order;
|
||
stats(order, 0) = order;
|
||
stats(order, 1) += 1;
|
||
}
|
||
return stats;
|
||
}
|
||
|
||
// 保存结果到文件
|
||
void SaveResults(
|
||
const std::vector<CoincidenceEvent>& results,
|
||
const arma::mat& stats,
|
||
const std::string& base_filename)
|
||
{
|
||
// 保存详细符合事件
|
||
std::ofstream detail_file(base_filename + "_detail.csv");
|
||
detail_file << "符合次数,时间窗口(ns),板卡号_通道号,能量,时间戳\n";
|
||
|
||
for (const auto& event : results) {
|
||
for (size_t i = 0; i < event.events.size(); ++i) {
|
||
const auto& data = event.events[i];
|
||
detail_file << event.coincidence_order << ","
|
||
<< event.time_window << ","
|
||
<< data.board_id << "_" << data.channel_id << ","
|
||
<< data.energy << ","
|
||
<< std::fixed << data.timestamp << "\n";
|
||
}
|
||
}
|
||
detail_file.close();
|
||
|
||
// 保存统计结果
|
||
std::ofstream stats_file(base_filename + "_stats.csv");
|
||
stats_file << "符合次数,事件数量\n";
|
||
|
||
for (int i = 2; i < stats.n_rows; ++i) {
|
||
if (stats(i, 1) > 0) {
|
||
stats_file << (int)stats(i, 0) << "," << (int)stats(i, 1) << "\n";
|
||
}
|
||
}
|
||
stats_file.close();
|
||
|
||
std::cout << "结果已保存到：" << base_filename << "_detail.csv 和 " << base_filename << "_stats.csv" << std::endl;
|
||
}
|
||
|
||
}
|
||
|
||
}
|