155 lines
4.8 KiB
C++
155 lines
4.8 KiB
C++
#include "hit_analyse_v2.h"
|
|
#include <random>
|
|
|
|
HIT_ANALYSE_V2::HIT_ANALYSE_V2(QObject *parent) : QObject(parent)
|
|
{
|
|
|
|
|
|
|
|
}
|
|
|
|
QString HIT_ANALYSE_V2::analyseBeamData(QVector<BufferData> dataframe){
|
|
|
|
double position=0.1;
|
|
double focus=8;
|
|
double intensity=10000.0;
|
|
QString dataString;
|
|
|
|
|
|
const int vector_length = 300; // Replace with the actual length of your vectors
|
|
|
|
std::vector<double> signal_list(vector_length);
|
|
std::vector<double> channel_list(vector_length);
|
|
|
|
|
|
// Create a random number generator with a Gaussian distribution
|
|
std::random_device rd;
|
|
std::mt19937 gen(rd());
|
|
std::normal_distribution<double> dist(0.0, 17.0); // Mean of 0 and Sigma of 17
|
|
|
|
// Create a vector to store the generated values
|
|
std::vector<short int> result(vector_length);
|
|
|
|
// Fill the vector with random values
|
|
for (int i = 0; i < vector_length; ++i) {
|
|
result[i] = static_cast<short int>(dist(gen));
|
|
signal_list.push_back(result[i]);
|
|
channel_list.push_back(i);
|
|
}
|
|
//add a gaussian profile, focus is FWHM, position is random between 50 and 250
|
|
position = 50 + (rand() % (int)(250 - 50 + 1));
|
|
for (int i = 0; i < vector_length; ++i) {
|
|
signal_list[i] += intensity*exp(-4*log(2)*pow((channel_list[i]-position)/focus,2));
|
|
}
|
|
|
|
|
|
|
|
|
|
// Fill signal_list and channel_list with your data
|
|
|
|
double SumT = 0.0, SumS = 0.0, SumS2 = 0.0, SumST = 0.0, SumT2 = 0.0, SumY = 0.0, SumYS = 0.0, SumYT = 0.0;
|
|
double b_den = 0.0, b_num = 0.0, b = 0.0, p = 0.0, c = 0.0, SumYYP = 0.0, SumYYM = 0.0, MeanY = 0.0;
|
|
double S[vector_length];
|
|
double T[vector_length];
|
|
|
|
for (int k = 0; k < vector_length; k++) {
|
|
if (k == 0) {
|
|
S[k] = 0.0;
|
|
T[k] = 0.0;
|
|
} else {
|
|
S[k] = S[k - 1] + 0.5 * (signal_list[k] + signal_list[k - 1]) * (channel_list[k] - channel_list[k - 1]);
|
|
T[k] = T[k - 1] + 0.5 * (channel_list[k] * signal_list[k] + channel_list[k - 1] * signal_list[k - 1]) *
|
|
(channel_list[k] - channel_list[k - 1]);
|
|
}
|
|
SumS += S[k];
|
|
SumT += T[k];
|
|
SumY += signal_list[k];
|
|
SumS2 += S[k] * S[k];
|
|
SumST += S[k] * T[k];
|
|
SumT2 += T[k] * T[k];
|
|
SumYS += signal_list[k] * S[k];
|
|
SumYT += signal_list[k] * T[k];
|
|
MeanY += signal_list[k];
|
|
}
|
|
MeanY /= vector_length;
|
|
|
|
// Calculate M1 matrix elements
|
|
double M1_00 = SumT2;
|
|
double M1_01 = SumST;
|
|
double M1_02 = SumT;
|
|
double M1_10 = SumST;
|
|
double M1_11 = SumS2;
|
|
double M1_12 = SumS;
|
|
double M1_20 = SumT;
|
|
double M1_21 = SumS;
|
|
double M1_22 = vector_length;
|
|
|
|
// Calculate M2 vector elements
|
|
double M2_0 = SumYT;
|
|
double M2_1 = SumYS;
|
|
double M2_2 = SumY;
|
|
|
|
// Calculate the inverse of M1
|
|
double detM1 = M1_00 * (M1_11 * M1_22 - M1_12 * M1_21) -
|
|
M1_01 * (M1_10 * M1_22 - M1_12 * M1_20) +
|
|
M1_02 * (M1_10 * M1_21 - M1_11 * M1_20);
|
|
|
|
if (detM1 == 0.0) {
|
|
std::cerr << "M1 is not invertible." << std::endl;
|
|
//return 1;
|
|
}
|
|
|
|
double invM1_00 = (M1_11 * M1_22 - M1_12 * M1_21) / detM1;
|
|
double invM1_01 = (M1_02 * M1_21 - M1_01 * M1_22) / detM1;
|
|
double invM1_02 = (M1_01 * M1_12 - M1_02 * M1_11) / detM1;
|
|
double invM1_10 = (M1_12 * M1_20 - M1_10 * M1_22) / detM1;
|
|
double invM1_11 = (M1_00 * M1_22 - M1_02 * M1_20) / detM1;
|
|
double invM1_12 = (M1_02 * M1_10 - M1_00 * M1_12) / detM1;
|
|
double invM1_20 = (M1_10 * M1_21 - M1_11 * M1_20) / detM1;
|
|
double invM1_21 = (M1_01 * M1_20 - M1_00 * M1_21) / detM1;
|
|
double invM1_22 = (M1_00 * M1_11 - M1_01 * M1_10) / detM1;
|
|
|
|
// Calculate ABC vector
|
|
double ABC_0 = invM1_00 * M2_0 + invM1_01 * M2_1 + invM1_02 * M2_2;
|
|
double ABC_1 = invM1_10 * M2_0 + invM1_11 * M2_1 + invM1_12 * M2_2;
|
|
double ABC_2 = invM1_20 * M2_0 + invM1_21 * M2_1 + invM1_22 * M2_2;
|
|
|
|
// Calculate b, p, and c
|
|
p = -ABC_0 / 2.0;
|
|
c = -ABC_1 / ABC_0;
|
|
|
|
for (int k = 0; k < vector_length; k++) {
|
|
double exp_term = exp(-p * (channel_list[k] - c) * (channel_list[k] - c));
|
|
b_num += exp_term * signal_list[k];
|
|
b_den += exp_term;
|
|
}
|
|
|
|
b = b_num / b_den;
|
|
|
|
for (int k = 0; k < vector_length; k++) {
|
|
double y_pred = b * exp(-p * (channel_list[k] - c) * (channel_list[k] - c));
|
|
SumYYM += (signal_list[k] - MeanY) * (signal_list[k] - MeanY);
|
|
SumYYP += (y_pred - MeanY) * (y_pred - MeanY);
|
|
}
|
|
|
|
double R_squared = SumYYP / SumYYM;
|
|
//std::cout << "R-squared = " << R_squared << endl;
|
|
position = -ABC_1/ ABC_0;
|
|
//sigma = sqrt(1.0 / (2.0 * ABC_0));
|
|
focus = 2.3548/sqrt(2*p);
|
|
intensity = b;
|
|
|
|
dataString += QString::number(intensity) + ',' + QString::number(position) + ',' + QString::number(focus)
|
|
+ ',' + QString::number(R_squared);
|
|
|
|
|
|
return dataString;
|
|
|
|
|
|
}
|
|
|
|
HIT_ANALYSE_V2::~HIT_ANALYSE_V2()
|
|
{
|
|
|
|
}
|