From f80b78e2931611e25eaf6702ae62d6b6773d3e4b Mon Sep 17 00:00:00 2001 From: Blake Leverington Date: Tue, 12 Sep 2023 16:26:47 +0200 Subject: [PATCH] added hit_analyse_v2 class. Currently produces and fits a fake data set. --- hit2023v2/eventbuilder.cpp | 62 +++-- hit2023v2/eventbuilder.h | 3 + hit2023v2/hit2023v2.pro | 5 +- hit2023v2/hit_analyse_v2.cpp | 154 ++++++++++++ hit2023v2/hit_analyse_v2.h | 72 ++++++ hit2023v2/hitreader.h | 441 +++++++++++++++++++++++++++++++++++ 6 files changed, 723 insertions(+), 14 deletions(-) create mode 100644 hit2023v2/hit_analyse_v2.cpp create mode 100644 hit2023v2/hit_analyse_v2.h create mode 100644 hit2023v2/hitreader.h diff --git a/hit2023v2/eventbuilder.cpp b/hit2023v2/eventbuilder.cpp index 212183c..25c4dad 100644 --- a/hit2023v2/eventbuilder.cpp +++ b/hit2023v2/eventbuilder.cpp @@ -1,5 +1,6 @@ #include "eventbuilder.h" #include "udpserver.h" +#include "hit_analyse_v2.h" EventBuilder::EventBuilder(QObject *parent) : QObject(parent) { @@ -10,6 +11,7 @@ EventBuilder::EventBuilder(QObject *parent) : QObject(parent) connect(this, EventBuilder::sigStartTakingHistos, this, EventBuilder::onStartTakingHistos); connect(this, EventBuilder::sigStopTakingHistos, this, EventBuilder::onStopTakingHistos); + moveToThread(&thread); thread.start(); init(); @@ -52,18 +54,57 @@ void EventBuilder::onNewData(DataReceiver* receiver) } + + + + //************ TODO ************ + //Here we can do something more with the complete frame + // I probably want to find the position and focus with the linear regression algorithm, but first, just send data to the udpserver to test. + + //currentFrame[dev_nr].sensor_data[ch] + //currentFrame is BufferData + // unsigned short* sensor_data; + + + //ToDo: + //1. Background subtraction. + + frame_counter++; + + while (frame_counter<10000){ + for (unsigned int dev_nr = 0; dev_nr < nrReceivers; dev_nr++){ + for (unsigned int ch = 0; ch < channelCounts[dev_nr]; ch++) + backgroundFrame[dev_nr].sensor_data[ch]+= currentFrame[dev_nr].sensor_data[ch]; + + } + } + if (frame_counter==10000){ + for (unsigned int dev_nr = 0; dev_nr < nrReceivers; dev_nr++){ + for (unsigned int ch = 0; ch < channelCounts[dev_nr]; ch++) + backgroundFrame[dev_nr].sensor_data[ch]/= 10000 ; + } + } + if (frame_counter>10000){ + for (unsigned int dev_nr = 0; dev_nr < nrReceivers; dev_nr++){ + for (unsigned int ch = 0; ch < channelCounts[dev_nr]; ch++) + currentFrame[dev_nr].sensor_data[ch]-=backgroundFrame[dev_nr].sensor_data[ch] ; + } + } + + + lastFrameMutex.lock(); if (newDataSemaphore.available() == 0) newDataSemaphore.release(1); lastFrame = currentFrame; lastFrameMutex.unlock(); - //histogram stuff + //histogram stuff if (histogramSamplesToTake) { for (int dev_nr = 0; dev_nr < nrReceivers; dev_nr++) for (int ch = 0; ch < channelCounts[dev_nr]; ch++) - histograms[baseAddresses[dev_nr] + ch].shoot(currentFrame[dev_nr].sensor_data[ch]); + histograms[baseAddresses[dev_nr] + ch].shoot(currentFrame[dev_nr].sensor_data[ch]); if (histogramSamplesToTake != -1) histogramSamplesToTake--; @@ -71,19 +112,13 @@ void EventBuilder::onNewData(DataReceiver* receiver) emit sigHistoCompleted(); } - //log data - if (loggingData) - logDataToFile(); + //log data + if (loggingData) logDataToFile(); + HIT_ANALYSE_V2 hit_analyse_v2;//create the object + QString dataString = hit_analyse_v2.analyseBeamData(currentFrame); - - //************ TODO ************ - //Here we can do something more with the complete frame - // I probably want to find the position and focus with the linear regression algorithm, but first, just send data to the udpserver to test. - intensity+=1.0; - position+=0.1; - focus+=0.01; // Call sendData method of the UDP server - QString dataString = QString::number(intensity) + ',' + QString::number(position) + ',' + QString::number(focus); + //QString dataString = QString::number(intensity) + ',' + QString::number(position) + ',' + QString::number(focus); QByteArray data = dataString.toUtf8(); udpServer.sendData(data); @@ -252,6 +287,7 @@ void EventBuilder::addSource(DataReceiver* source) receivers.append(source); nrReceivers = receivers.length(); currentFrame.resize(nrReceivers); + backgroundFrame.resize(nrReceivers); connect(source, DataReceiver::sigDataReady, this, EventBuilder::onNewData); } diff --git a/hit2023v2/eventbuilder.h b/hit2023v2/eventbuilder.h index f4567c0..17b378d 100644 --- a/hit2023v2/eventbuilder.h +++ b/hit2023v2/eventbuilder.h @@ -64,6 +64,8 @@ protected: QVector receivers; QVector currentFrame; + QVector backgroundFrame; + QVector lastFrame; QVector histograms; @@ -88,6 +90,7 @@ protected slots: void onStartTakingHistos(int sample_count); void onStopTakingHistos(); private: + long unsigned int frame_counter = 0; double intensity = 0.0; double position = 0.0; double focus = 0.0; diff --git a/hit2023v2/hit2023v2.pro b/hit2023v2/hit2023v2.pro index 70b1934..d2c4a5c 100644 --- a/hit2023v2/hit2023v2.pro +++ b/hit2023v2/hit2023v2.pro @@ -15,6 +15,7 @@ TEMPLATE = app SOURCES += main.cpp\ + hit_analyse_v2.cpp \ mainwindow.cpp \ q_debugstream.cpp \ dialoglogsettings.cpp \ @@ -45,6 +46,7 @@ HEADERS += mainwindow.h \ device.h \ dev_commands.h \ datareceiver.h \ + hit_analyse_v2.h \ hw.h \ dialoghostip.h \ dialogtriggersettings.h \ @@ -63,7 +65,8 @@ HEADERS += mainwindow.h \ dialogprofiler.h \ stepper.h \ dialogbeta.h \ - udpserver.h + udpserver.h \ + hitreader.h FORMS += mainwindow.ui \ dialoglogsettings.ui \ diff --git a/hit2023v2/hit_analyse_v2.cpp b/hit2023v2/hit_analyse_v2.cpp new file mode 100644 index 0000000..6a096f4 --- /dev/null +++ b/hit2023v2/hit_analyse_v2.cpp @@ -0,0 +1,154 @@ +#include "hit_analyse_v2.h" +#include + +HIT_ANALYSE_V2::HIT_ANALYSE_V2(QObject *parent) : QObject(parent) +{ + + + +} + +QString HIT_ANALYSE_V2::analyseBeamData(QVector 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 signal_list(vector_length); + std::vector channel_list(vector_length); + + + // Create a random number generator with a Gaussian distribution + std::random_device rd; + std::mt19937 gen(rd()); + std::normal_distribution dist(0.0, 17.0); // Mean of 0 and Sigma of 17 + + // Create a vector to store the generated values + std::vector result(vector_length); + + // Fill the vector with random values + for (int i = 0; i < vector_length; ++i) { + result[i] = static_cast(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() +{ + +} diff --git a/hit2023v2/hit_analyse_v2.h b/hit2023v2/hit_analyse_v2.h new file mode 100644 index 0000000..4aadf1a --- /dev/null +++ b/hit2023v2/hit_analyse_v2.h @@ -0,0 +1,72 @@ +#ifndef HIT_ANALYSE_V2_H +#define HIT_ANALYSE_V2_H + +#include +#include +#include +#include +#include +#include // std::count +#include +#include +#include +#include // std::setprecision +#include +#include +#include +#include +#include "datareceiver.h" +//#include +//#include +//#include +//#include + +class HIT_ANALYSE_V2 : public QObject +{ + Q_OBJECT +public: + explicit HIT_ANALYSE_V2(QObject *parent = nullptr); + ~HIT_ANALYSE_V2(); +private: + struct beamRecon { + + double Position; + double Focus; + double Peak; + double Rsqr; + double Skew; + double Kurtosis; + double Sum; + int n_channels; + }; + + struct bpm_frame_v2 { + double channel_amp[320]; + double channel_position[320]; + double avg_position; + double avg_width; + double integratedsignalamp; + double maxchannel_amp; + int maxchannel; + int nrChannels; + int board_number; + int sma_state; + + }; + + +public slots: + QString analyseBeamData(QVector dataframe); + // void processPendingDatagrams(); + +}; + +class Channel{ +public: + double amplitude;; + double position; + int chnumber; + int last_neighbour; + int next_neighbour; +}; +#endif // HIT_ANALYSE_V2_H diff --git a/hit2023v2/hitreader.h b/hit2023v2/hitreader.h new file mode 100644 index 0000000..8bb465d --- /dev/null +++ b/hit2023v2/hitreader.h @@ -0,0 +1,441 @@ +//This is an object interface for reading HIT data files +//See HIT documentation for details and examples. +/* + THIS DOESN'T WORK! + .L hitreader.c + + Hitdata data; + data.read(“my_file.da2”); //to load whole file at once – forget it! See below. + data.read(“my_file.da2”,1000,100,10) //to read 100 frames starting from frame 1000 and incrementing by 10 (i.e. frame 1000, 1010, 1020, ... 1990 will be read) + //Reading 10 000 frames is reasonable. Reading 100 000 frames made my VM beg for memory. + + data.nrFrames //to see how many frames you have + data.frames[0].nrBoards //to see how many boards you had in the system + data.frames[0].boards[0].nrChannels //to see how many channels you have in board 0 + data.frames[10].boards[0].data[100] //get signal value for frame 10, board 0, channel 100 + data.frames[10].boards[0].syncframe.local_ctr //get the local synchro counter for frame 10, board 0 + //same for .global_ctr, .sma_state, .dummy, .device_nr, .data_ok + + */ + +//*********************** Helper ************************* +#include +#include + +using namespace std; + +//#define debug(str) std::cout << "HIT DEBUG: " << str << endl; +#define debug(str) + +//*********************** Syncframe ************************* +class Syncframe +{ + public: + Syncframe() + { + debug("Syncframe()"); + + local_ctr = global_ctr = 0; + sma_state = dummy = 0; + device_nr = -1; + data_ok = 0; + }; + + ~Syncframe() + { + debug("~Syncframe()"); + }; + + int sizeInFile() + { + return 16; + }; + + int read(std::ifstream* file) + { + char buffer[16]; + file->read(buffer,16); + if (file->fail()) + return 0; + local_ctr = *(unsigned short*)(buffer+0); + global_ctr = *(unsigned short*)(buffer+2); + sma_state = *(unsigned short*)(buffer+4); + dummy = *(unsigned short*)(buffer+6); + device_nr = *(int*)(buffer+8); + data_ok = *(int*)(buffer+12); + // std::cout << "Syncframe:" << local_ctr << " " << global_ctr << " " << sma_state << " " << dummy << " " << device_nr << " " << data_ok << std::endl; + + return 1; + }; + + + unsigned short local_ctr; + unsigned short global_ctr; + unsigned short sma_state; + unsigned short dummy; + int device_nr; + unsigned int data_ok; +}; + +//*********************** Sensorframe ************************* +class Boardframe +{ + public: + Boardframe(int nr_channels = 0) + { + debug("Boardframe()"); + + data = NULL; + resize (nr_channels); + }; + + Boardframe(const Boardframe& in) + { + debug("Boardframe(Boardframe&)"); + + data = NULL; + resize(in.nrChannels); + for (int i = 0; i < nrChannels; i++) + data[i] = in.data[i]; + syncframe = in.syncframe; + }; + + Boardframe& operator=(const Boardframe& in) + { + debug("Boardframe::operator=="); + + resize(in.nrChannels);//creates an array called data of length nrChannels + for (int i = 0; i < nrChannels; i++) + data[i] = in.data[i]; + syncframe = in.syncframe; + return *this; + }; + + ~Boardframe() + { + debug("~Boardframe()"); + + if (data) + delete[] data; + }; + + void resize(int nr_channels) + { + if (data) + delete[] data; + nrChannels = nr_channels; + if (nrChannels) + data = new unsigned short[nrChannels]; + else + data = NULL; + }; + + int sizeInFile() + { + // std::cout << "boardframe.sizeInFile() = " << syncframe.sizeInFile() + nrChannels*2 << std::endl; + return syncframe.sizeInFile() + nrChannels*2; + + }; + + int read(std::ifstream* file) + { + if (syncframe.read(file) == 0)//get the syncframe before the board data + return 0; + //I must be already resized at this point! + file->read((char*)data,2*nrChannels); + if (file->fail()) + return 0; + // std::cout<< "data[" << nrChannels << "]: "; + // for (int i = 0;iread((char*)&nr_boards,2); + if (nr_boards==4) break; + } + if ( nr_boards!=4){ + std::cerr << "Resync failed." << std::endl; + return 0; + } + + } + //std::cout << " nr_boards: " << nr_boards << std::endl; + //Read channel counts + unsigned short* channel_counts = new unsigned short[nr_boards]; + file->read((char*)channel_counts,nr_boards*2); + if (file->fail()) + { + delete[] channel_counts; + return 0; + } + + + + //Read board frames + resize(nr_boards); + for (int board_nr = 0; board_nr < nr_boards; board_nr++) + { + // std::cout << " channel_counts[" << board_nr << "]: "<< channel_counts[board_nr] << std::endl; + + boards[board_nr].resize(channel_counts[board_nr]); + if (boards[board_nr].read(file) == 0)//read the board + { + delete[] channel_counts; + return 0; + } + } + + delete[] channel_counts; + return 1; + }; + + int nrChannels() + { + int result = 0; + for (int board_nr = 0; board_nr < nrBoards; board_nr++) + result += boards[board_nr].nrChannels; + return result; + }; + + unsigned short& operator[] (int index) + { + for (int board_nr = 0; board_nr < nrBoards; board_nr++) + { + if (index >= boards[board_nr].nrChannels) + index -= boards[board_nr].nrChannels; + else + return boards[board_nr][index]; + } + + std::cerr << " ### Fullframe::operator[]: index out of range!" << std::endl; + // return (*NULL); //this will cause crash (intended). + return boards[nrBoards][index]; + }; + + int nrBoards; + Boardframe* boards; +}; + +//*********************** Hitdata ************************* + +class Hitdata +{ + public: + Hitdata(int nr_frames = 0) + { + frames = NULL; + resize(nr_frames); + }; + + Hitdata(const Hitdata& in) + { + frames = NULL; + resize(in.nrFrames); + for (int i = 0; i < nrFrames; i++) + frames[i] = in.frames[i]; + }; + + Hitdata& operator=(const Hitdata& in) + { + resize(in.nrFrames); + for (int i = 0; i < nrFrames; i++) + frames[i] = in.frames[i]; + + return *this; + }; + + ~Hitdata() + { + if (nrFrames) + delete[] frames; + }; + + void resize (int nr_frames) + { + if (nrFrames) + delete[] frames; + nrFrames = nr_frames; + if (nrFrames) + frames = new Fullframe[nrFrames]; + else + frames = NULL; + }; + + //Read data from a given file. + //first_frame is the number of first frame to be read + //nr_frames is the maximum number of frames to be read + //-1 to read all of them + //increment allows you reading once every nth sample + //Return number of frames read or 0 in case of failure + int readFile(char* filename, int first_frame = 1, int nr_frames = -1, int increment = 1) + { + std::ifstream file; + //Open the file + file.open(filename, ios_base::in | ios_base::binary); + if (!file.is_open()) + { + std::cerr << " ### Hitdata: File could not be open!" << std::endl; + return 0; //file could not be opened + } + + //Read first record to find board configuration + Fullframe sampleframe; + if (sampleframe.read(&file) == 0) + { + std::cerr << " ### Hitdata: First frame could not be read!" << std::endl; + file.close(); + return 0; + } + else { + std::cout << "Sample frame size (bytes): " << sampleframe.sizeInFile() << std::endl; + } + + //Check file size + file.seekg(0, std::ios::beg); + std::streamsize fsize = file.tellg(); + file.seekg(0, std::ios::end); + fsize = file.tellg() - fsize; + + //Determine real frames to read + unsigned int max_frames = fsize / sampleframe.sizeInFile(); + if ((max_frames == -1) || (max_frames < nr_frames)) + nr_frames = max_frames; + + std::cout << " Hitdata: Nr frames to be read: " << nr_frames << std::endl; + + //Read! + // resize(nr_frames); //make an array of Fullframes called frames of size nr_frames + file.seekg(first_frame * sampleframe.sizeInFile(), std::ios::beg); + for (int frame_nr = first_frame; frame_nr < nr_frames; frame_nr++) + { + if ((frame_nr%10000) == 0) + std::cout << " Frame " << frame_nr << std::endl; + + file.seekg((frame_nr*increment) * sampleframe.sizeInFile() , std::ios::beg); + if (file.eof()) { + std::cerr<< "end of file reached." << std::endl; + return frame_nr; + } + if ( sampleframe.read(&file) == 0) //read the next frame + { + std::cerr << " ### Hitdata: Frame " << frame_nr << " could not be read!" << std::endl; + file.close(); //read error, finish! + // frames = frame_nr; //Kinky! We decrease nr_frames, but the actual array size remains unchanged! + ///???? I don't know what the above line does. + return frame_nr; + } + // std::cout << frames[frame_nr].nrBoards << std::endl; + } + + //Finished + file.close(); + return nr_frames; + }; + + Fullframe& operator[] (int index) + { + if (index < nrFrames) + return frames[index]; + else + { + std::cerr << " ### Hitdata::operator[]: index out of range!" << std::endl; + // return (*NULL); //this will cause crash (intended). + return frames[index]; + + } + }; + + int nrFrames; + Fullframe* frames; +}; + + + + + +