worker.cpp 18.90 KiB
#include "worker.h"
#include "TSystem.h"
#include <QFile>
#include <QString>
#include <QTextStream>
#include <QDataStream>
#include <QDebug>
#include <QDir>
bool Worker::processCoincidence(const std::vector<MyPulse>& pulses, Cone& newCone){
// fill energy spectrum
// TODO
if (pulses.size() != 2)
{
return false;
}
// energy cut
double Etot = pulses[1].height + pulses[0].height ;
if (config->energyCut)
{
if (Etot < config->energyLow || Etot > config->energyUp)
{
return false;
}
}
// if (pulses[0].height < 0.01 || pulses[1].height < 0.01)
// return false;
double cosHalfangle = 1 - pulses[0].height * 511 / (pulses[1].height * Etot);
if (std::abs(cosHalfangle)>=1)
{
return false;
}
Vector3D axis = pulses[1].pos - pulses[0].pos;
if (axis * axis == 0){
return false;
}
newCone = Cone(pulses[0].pos, axis, cosHalfangle,
Etot, pulses[0].height);
return true;
}
Worker::Worker(QObject *parent, const Setup* config_, PlotData* plotdata_) :
QThread(parent),
config(config_),
plotdata(plotdata_)
{
// stopped=true;
stopped = false;
exitted=false;
enabledChs = EnabledChannel(config->channelSettings);
infile.setFileName(config->filePath);
// qDebug()<< outfp << '\n';
if (config->inputFormat == "BIN")
{
if (!infile.open(QIODevice::ReadOnly))
{
throw std::invalid_argument("Cannot open file: " + config->filePath.toStdString());
}
}
else if(config->inputFormat == "CSV FULL" || config->inputFormat == "CSV DECODE EVENTS")
{
if (!infile.open(QIODevice::ReadOnly | QIODevice::Text))
{
throw std::invalid_argument("Cannot open file: " + config->filePath.toStdString());
}
}
QString outfp = config->outputDir.absolutePath() + "/cones.txt";
outfile.setFileName(outfp);
if (!outfile.open(QIODevice::WriteOnly | QIODevice::Text))
{
throw std::invalid_argument("Cannot open file: " + outfp.toStdString());
}
else {
QTextStream out(&outfile);
out << qSetFieldWidth(8)
<< "apex.x(cm)" // << " "
<< "apex.y(cm)" // << " "
<< "apex.z(cm)" // << " "
<< "axis.x(cm)" // << " "
<< "axis.y(cm)" // << " "
<< "axis.z(cm)" // << " "
<< qSetFieldWidth(13)
<< "CosHalfAngle"
<< "Etot(keV)"
<< "E_1(keV)" << '\n';
}
}
void Worker::run()
{
QTextStream out(&outfile);
if (config->inputFormat == "BIN")
{
QDataStream in(&infile);
while(!exitted)
{
if (!stopped)
{
// read some cones
std::vector<Cone> cones;
readBIN(cones, in);
if (exitted)
break;
// update image
plotdata->updateImage(cones.cbegin(), cones.cend(), true);
// save cones
saveCones(cones, out);
if (exitted || plotdata->counts >= config->maxN)
break;
}
}
}
else if (config->inputFormat == "CSV FULL")
{
QTextStream in(&infile);
while(!exitted)
{
if (!stopped)
{
// read some cones
std::vector<Cone> cones;
readCSVFULL(cones, in);
if (exitted)
break;
// update image
plotdata->updateImage(cones.cbegin(), cones.cend(), true);
// save cones
saveCones(cones, out);
if (exitted || plotdata->counts >= config->maxN)
break;
}
}
}
else if (config->inputFormat == "CSV DECODE EVENTS") {
QTextStream in(&infile);
while(!exitted)
{
if (!stopped)
{
// read some cones
std::vector<Cone> cones;
readCSVDECODEEVENTS(cones, in);
if (exitted)
break;
// update image
plotdata->updateImage(cones.cbegin(), cones.cend(), true);
// save cones
saveCones(cones, out);
if (exitted || plotdata->counts >= config->maxN)
break;
}
}
}
// save image & spectrum before quitting
qDebug() << "Save results before exiting worker thread";
plotdata->saveOnExit();
qDebug() << "Worker thread exited.";
}
void Worker::readCSVFULL(std::vector<Cone> &cones, QTextStream &in)
{
QString line;
static std::vector<DataCITIROC> packets;
static double prevTT=0;
double currTT=0;
// skip header
in.readLineInto(&line);
while (!exitted && cones.size() < config->chuckSize) {
while (!exitted && cones.size() < config->chuckSize && in.readLineInto(&line)) {
// split line into data
QVector<QStringRef> data = line.splitRef(";");
if(data.size() < 4)
break;
int numPackets = data[3].toInt();
int expectedLen = 102*data[3].toInt() + 4;
if(data.size() != expectedLen)
break;
currTT = data[1].toDouble();
if (packets.size()>0 && std::abs(currTT - prevTT) >= config->timeWindow){
processPackets(packets, cones);
packets.clear();
}
prevTT = currTT;
// parse data
int startIndex = 4;
for (int j=0;j<numPackets;j++)
{
DataCITIROC newPacket(data.cbegin() + startIndex);
startIndex += 102;
// update LG HG Spectra
plotdata->updateLGSpectra(newPacket.AsicID, newPacket.chargeLG);
plotdata->updateHGSpectra(newPacket.AsicID, newPacket.chargeHG);
packets.push_back(std::move(newPacket));
}
// qDebug() << newPacket.EventCounter << newPacket.RunEventTimecode << '\n';
}
// wait for new data
if (in.atEnd())
msleep(1000);
}}
void Worker::readCSVDECODEEVENTS(std::vector<Cone> &cones, QTextStream &in)
{
QString line;
static std::vector<DataCITIROC> packets;
static double prevTT=0;
double currTT=0;
// skip header
in.readLineInto(&line);
while (!exitted && cones.size() < config->chuckSize) {
while (!exitted && cones.size() < config->chuckSize && in.readLineInto(&line)) {
// split line into data
QVector<QStringRef> data = line.splitRef(";");
// int expectedLen = 103;
if(data.size() != 103)
break;
currTT = data[4].toDouble();
if (packets.size()>0 && std::abs(currTT - prevTT) >= config->timeWindow){
processPackets(packets, cones);
packets.clear();
}
prevTT = currTT;
// parse data
DataCITIROC newPacket(data.cbegin() + 1);
// update LG HG Spectra
plotdata->updateLGSpectra(newPacket.AsicID, newPacket.chargeLG);
plotdata->updateHGSpectra(newPacket.AsicID, newPacket.chargeHG);
packets.push_back(std::move(newPacket));
// qDebug() << newPacket.EventCounter << newPacket.RunEventTimecode << '\n';
}
// wait for new data
if (in.atEnd())
msleep(1000);
}
}
//void Worker::readCSV(std::vector<Cone> &cones, QTextStream& in)
//{
// QString line;
// std::vector<ulong> histNums;
// static std::vector<MyPulse> pulses;
// static double prevTT=0;
// int i=0;
// while (!exitted && i < config->chuckSize && in.readLineInto(&line))
// {
// MyPulse newPulse = MyPulse(line);
// auto it = enabledChs.indexOf(newPulse.cellNo);
// if (it < 0){
// // cell not found
// // go to next line
// continue;
// }
// else {
// newPulse.pos.X = config->channelSettings[it].x;
// newPulse.pos.Y = config->channelSettings[it].y;
// newPulse.height *= config->channelSettings[it].caliCoef;
// // TO DO: calculate z
// }
// if (pulses.size()==0 || std::abs(newPulse.timeStamp - prevTT) < config->timeWindow){
// prevTT = newPulse.timeStamp;
// pulses.push_back(std::move(newPulse));
// }
// else {
// // process coincident events
// Cone newCone;
// if (processCoincidence(pulses, newCone)){
// histNums.push_back(pulses[0].histNo);
// cones.push_back(std::move(newCone));// i++;
// }
// pulses.clear();
// prevTT = newPulse.timeStamp;
// pulses.push_back(std::move(newPulse));
// }
// }
// // save cones to file
// for (int i=0;i<cones.size();i++)
// {
// outfile << std::fixed << std::setprecision(2)
// << std::setw(8) << cones[i].apex.X // << " "
// << std::setw(8) << cones[i].apex.Y // << " "
// << std::setw(8) << cones[i].apex.Z // << " "
// << std::setw(8) << cones[i].axis.X // << " "
// << std::setw(8) << cones[i].axis.Y // << " "
// << std::setw(8) << cones[i].axis.Z // << " "
// << std::fixed << std::setprecision(8)
// << std::setw(13) << cones[i].cosHalfAngle
// << std::setw(13) << cones[i].E0
// << std::setw(13) << cones[i].Edpst
// << std::setw(13) << histNums[i] << '\n';
// }
//}
void Worker::readBIN(std::vector<Cone> &cones, QDataStream &in)
{
// hold a single event
char buffer[eventSize];
static std::vector<DataCITIROC> packets;
static double prevTT=0;
// int i = 0;
// discard the garbage data at the begninning
const uint32_t *p_ui32 = (uint32_t*)&(buffer[0]);
while (true)
{
in.readRawData(buffer, 4);
uint32_t bword = ((uint32_t)(p_ui32[0]));
if (((bword >> 4) & 0xc000000) == 0x8000000)
break;
}
int len = 4;
while (!exitted && cones.size() < config->chuckSize) {
while (!exitted && cones.size() < config->chuckSize) {
// read data into buffer
if (len != eventSize)
{
int tmp=in.readRawData(buffer + len, eventSize-len);
if(tmp < 0)
break;
len += tmp;
}
if(len != eventSize)
break;
else
len = 0;
// parse data
DataCITIROC newPacket;
DecodeCITIROCRowEvents(buffer, eventSize, newPacket);
// update LG HG Spectra
plotdata->updateLGSpectra(newPacket.AsicID, newPacket.chargeLG);
plotdata->updateHGSpectra(newPacket.AsicID, newPacket.chargeHG);
if (packets.size()>0 && std::abs(newPacket.RunEventTimecode_ns - prevTT) >= config->timeWindow){
processPackets(packets, cones);
packets.clear();
} prevTT = newPacket.RunEventTimecode_ns;
packets.push_back(std::move(newPacket));
// qDebug() << newPacket.EventCounter << newPacket.RunEventTimecode << '\n';
}
// wait for new data
if (in.atEnd())
msleep(1000);
}
}
void Worker::processPackets(const std::vector<DataCITIROC>& packets, std::vector<Cone> &cones)
{
// if(packets.size() > 2)
// qDebug() << "Coincidence";
// process coincident events
std::vector<Event> events;
if (extractEvents(packets, events))
{
// sort by channel id
std::sort(events.begin(), events.end(),
[](const Event & a, const Event & b) -> bool
{
return a.id < b.id;
});
std::vector<MyPulse> pulses;
if(event2Pulses(events, pulses))
{
// upddate pair spectra
plotdata->updatePairSpectra(pulses);
if(config->coincidenceEnabled)
{
Cone newCone;
if (processCoincidence(pulses, newCone)){
cones.push_back(std::move(newCone));
}
}
}
}
}
int Worker::DecodeCITIROCRowEvents(const char* buf, const int bufSize, DataCITIROC& DataCITIROCA)
{
//Convert char buffer in uin32
const uint32_t* bufferA = (uint32_t*)&(buf[0]);
const uint32_t valid_wordA = bufSize/4;
int DecodedPackets = 0;
double Time =0;
double minTime = 0;
uint32_t i, j, t, s;
uint datarow[97];
t= 0;
s=0;
while (t < valid_wordA)
{
switch (s) {
case 0:
if(((bufferA[t] >> 4) & 0xc000000) == 0x8000000)
{
s = 1;
DataCITIROCA.AsicID = (uint16_t) (bufferA[t] & 0xF);
DataCITIROCA.EventTimecode = ((uint64_t)bufferA[t + 1]);
DataCITIROCA.RunEventTimecode = (((uint64_t)bufferA[t+ 2]) ) + (((uint64_t)bufferA[t + 3]) << 32);
DataCITIROCA.EventCounter = ((uint64_t)bufferA[t + 4]);
DataCITIROCA.EventTimecode_ns = DataCITIROCA.EventTimecode * 25;
DataCITIROCA.RunEventTimecode_ns =DataCITIROCA.RunEventTimecode * 25;
t = t + 5; minTime = 100000000000000;
}
else
t++;
break;
case 1:
for(i = 0; i < 32; i++)
{
uint32_t bword = bufferA[t];
datarow[i * 3 + 0] = (bword >> 0) & 0x3FFF;
datarow[i * 3 + 1] = (bword >> 14) & 0x3FFF;
datarow[i * 3 + 2] = (bword >> 28) & 0x1;
t++;
}
for(i=0;i<32;i++)
{
DataCITIROCA.hit[31-i] = (bool)((datarow[i*3+2] & 0x1) == 1 ? true: false);
int dataHG = (int)datarow[(i * 3) + 0];
int dataLG= (int)datarow[(i * 3) + 1];
// if(dataHG > ThresholdSoftware)
// DataCITIROCA.chargeHG[i] = (ushort)dataHG;
// else
// DataCITIROCA.chargeHG[i] = 0;
// if(dataHG > ThresholdSoftware)
// DataCITIROCA.chargeLG[i] = (ushort)dataLG;
// else
// DataCITIROCA.chargeLG[i] = 0;
DataCITIROCA.chargeHG[31-i] = (ushort)dataHG;
DataCITIROCA.chargeLG[31-i] = (ushort)dataLG;
}
s = 2;
break;
case 2:
if((bufferA[t] & 0xc0000000) == 0xc0000000)
{
// pC.Enqueue(DataCITIROCA);
DecodedPackets++;
}
t++;
s = 0;
break;
}
}
return DecodedPackets;
}
bool Worker::extractEvents(const std::vector<DataCITIROC>& packets, std::vector<Event>& events)
{
// pick channels with hit = 1 && value > threshold
for (auto it = packets.cbegin(); it!=packets.cend();it++)
{
for (int i=0;i<32;i++)
{
if(it->hit[i] == 0)
continue;
int id = enabledChs.getID(qMakePair(it->AsicID, i));
if (id < 0)
continue;
Event newEvent;
double energy;
if(config->channelSettings[id].source == "LG")
{
energy = it->chargeLG[i] * config->channelSettings[id].caliCoef;
if(energy < config->channelSettings[id].threshold)
continue;
newEvent.energyLG = energy;
newEvent.energyHG = it->chargeHG[i]; }
else if (config->channelSettings[id].source == "HG")
{
energy = it->chargeHG[i] * config->channelSettings[id].caliCoef;
if(energy < config->channelSettings[id].threshold)
continue;
newEvent.energyLG = it->chargeLG[i];
newEvent.energyHG = energy;
}
newEvent.id = id;
newEvent.AsicID = it->AsicID;
newEvent.channelNum = i;
newEvent.EventCounter = it->EventCounter;
newEvent.RunEventTimecode_ns = it->RunEventTimecode_ns;
events.push_back(std::move(newEvent));
}
}
if (events.size() == 0)
return false;
return true;
}
bool Worker::event2Pulses(const std::vector<Event>& events, std::vector<MyPulse>& pulses)
{
// find pair and construct pulses
for (auto it = events.begin(); it != events.end(); it++)
{
auto next = it+1;
if (it->id % 2 == 1)
continue;
if (next == events.end())
continue;
if (next->id != it->id + 1)
continue;
// we have a pair (it, next)
MyPulse newPulse;
double e1 = (config->channelSettings[it->id].source == "LG") ? it->energyLG : it->energyHG;
double e2 = (config->channelSettings[next->id].source == "LG") ? next->energyLG : next->energyHG;
double z = findZ(e1, e2, config->channelSettings[it->id].z, config->channelSettings[next->id].z);
newPulse.height = e1 + e2;
newPulse.pos = Vector3D(config->channelSettings[it->id].x, config->channelSettings[it->id].y, z);
newPulse.timeStamp = it->RunEventTimecode_ns;
newPulse.cellNo = it->id / 2;
newPulse.histNo = it->EventCounter;
pulses.push_back(std::move(newPulse));
}
if(pulses.size() == 0)
return false;
std::sort(pulses.begin(), pulses.end(),
[](const MyPulse & a, const MyPulse & b) -> bool
{
return a.timeStamp < b.timeStamp;
});
return true;
}
void Worker::saveCones(const std::vector<Cone>& cones, QTextStream& out)
{
for (const Cone& p : cones)
{
out << qSetFieldWidth(8)
<< QString::number(p.apex.X, 'f', 3) // << " "
<< QString::number(p.apex.Y, 'f', 3) // << " "
<< QString::number(p.apex.Z, 'f', 3) // << " "
<< QString::number(p.axis.X, 'f', 3) // << " "
<< QString::number(p.axis.Y, 'f', 3) // << " "
<< QString::number(p.axis.Z, 'f', 3) // << " "
<< qSetFieldWidth(13)
<< QString::number(p.cosHalfAngle, 'f', 5) << QString::number(p.E0, 'f', 5)
<< QString::number(p.Edpst, 'f', 5) << '\n';
}
}
double Worker::findZ(const double & E1, const double & E2, const double & Z1, const double & Z2)
{
// linear interpolation
return (E1 * Z1 + E2 * Z2) / (E1 + E2);
}
void Worker::handleStart()
{
stopped=false;
}
void Worker::handlePause()
{
stopped=true;
}
void Worker::handleStop()
{
exitted = true;
qDebug() << "Exiting thread signal received.";
}