Skip to content

Commit

Permalink
Update peak finder process (#141)
Browse files Browse the repository at this point in the history 
* add observables

* update window observables

* more window obs

* parse map

* convert to physical time

* adc to energy using map

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* only one energy parameter can be defined

* bug

* fix bad observables

* electronics mode

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
  • Loading branch information
@lobis @pre-commit-ci
lobis and pre-commit-ci[bot] authored Nov 5, 2024
1 parent 52aeac4 commit f1b939d
Show file tree
Hide file tree
Showing 2 changed files with 192 additions and 51 deletions.
21 changes: 15 additions & 6 deletions inc/TRestRawPeaksFinderProcess.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,14 +19,23 @@ class TRestRawPeaksFinderProcess : public TRestEventProcess {
Double_t fThresholdOverBaseline = 2.0;
/// \brief range of samples to calculate baseline for peak finding
TVector2 fBaselineRange = {0, 10};
/// \brief distance between two peaks to consider them as different
/// \brief distance between two peaks to consider them as different (ADC units)
UShort_t fDistance = 10;
/// \brief window size to calculate the peak amplitude
/// \brief window size to calculate the peak amplitude (time bins)
UShort_t fWindow = 10;
/// \brief option to remove all veto signals after finding the peaks
Bool_t fRemoveAllVetos = false;
/// \brief option to remove peakless veto signals after finding the peaks
Bool_t fRemovePeaklessVetos = false;
Bool_t fRemoveAllVetoes = false;
/// \brief option to remove peak-less veto signals after finding the peaks
Bool_t fRemovePeaklessVetoes = false;

Double_t fTimeBinToTimeFactorMultiplier = 0.0;
Double_t fTimeBinToTimeFactorOffset = 0.0;
Double_t fTimeBinToTimeFactorOffsetTCM = 0.0;

Bool_t fTimeConversionElectronics = false;

Double_t fADCtoEnergyFactor = 0.0;
std::map<UShort_t, Double_t> fChannelIDToADCtoEnergyFactor = {};

std::set<std::string> fChannelTypes = {}; // this process will only be applied to selected channel types

Expand All @@ -49,7 +58,7 @@ class TRestRawPeaksFinderProcess : public TRestEventProcess {
TRestRawPeaksFinderProcess() = default;
~TRestRawPeaksFinderProcess() = default;

ClassDefOverride(TRestRawPeaksFinderProcess, 4);
ClassDefOverride(TRestRawPeaksFinderProcess, 5);
};

#endif // REST_TRESTRAWPEAKSFINDERPROCESS_H
222 changes: 177 additions & 45 deletions src/TRestRawPeaksFinderProcess.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) {
exit(1);
}

vector<tuple<UShort_t, UShort_t, double>> eventPeaks;
vector<tuple<UShort_t, UShort_t, double>> eventPeaks; // signalId, time, amplitude

for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) {
const auto signal = fInputEvent->GetSignal(signalIndex);
Expand Down Expand Up @@ -71,91 +71,174 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) {
return tie(get<1>(a), get<0>(a)) < tie(get<1>(b), get<0>(b));
});

// SetObservableValue("peaks", eventPeaks); // problems with dictionaries
vector<UShort_t> peaksChannelId;
vector<UShort_t> peaksTime;
vector<double> peaksAmplitude;
double peaksEnergy;
double amplitudeTotal = 0.0;

double peaksEnergy = 0.0;
UShort_t peaksCount = 0;
UShort_t peaksCountUnique = 0;

set<UShort_t> uniquePeaks;
for (const auto& [channelId, time, amplitude] : eventPeaks) {
peaksChannelId.push_back(channelId);
peaksTime.push_back(time);
peaksAmplitude.push_back(amplitude);

amplitudeTotal += amplitude;
}
peaksEnergy += amplitude;
peaksCount++;

peaksEnergy = amplitudeTotal;
if (uniquePeaks.find(channelId) == uniquePeaks.end()) {
uniquePeaks.insert(channelId);
peaksCountUnique++;
}
}

SetObservableValue("peaksChannelId", peaksChannelId);
SetObservableValue("peaksTime", peaksTime);
SetObservableValue("peaksAmplitude", peaksAmplitude);
SetObservableValue("totalPeaksEnergy", peaksEnergy);
SetObservableValue("peaksTimeBin", peaksTime);
SetObservableValue("peaksAmplitudeADC", peaksAmplitude);

SetObservableValue("peaksAmplitudeADCSum", peaksEnergy);
SetObservableValue("peaksCount", peaksCount);
SetObservableValue("peaksCountUnique", peaksCountUnique);

vector<UShort_t> windowIndex(eventPeaks.size(), 0); // Initialize with zeros
vector<UShort_t> windowCenter; // for each different window, the center of the window
vector<UShort_t> windowPeakIndex(eventPeaks.size(), std::numeric_limits<UShort_t>::max());
vector<UShort_t> windowPeakCenter(eventPeaks.size(), 0);
vector<UShort_t> windowPeakMultiplicity(eventPeaks.size(), 0);

UShort_t window_index = 0;
for (size_t peakIndex = 0; peakIndex < eventPeaks.size(); peakIndex++) {
const auto& [channelId, time, amplitude] = eventPeaks[peakIndex];
const auto windowTime = time - fWindow / 2;
const auto windowEnd = time + fWindow / 2;
const auto windowTimeStart = time - fWindow / 2;
const auto windowTimeEnd = time + fWindow / 2;

// check if the peak is already in a window
if (windowIndex[peakIndex] != 0) {
if (windowPeakIndex[peakIndex] != std::numeric_limits<UShort_t>::max()) {
continue;
}

// create a new window
windowCenter.push_back(time);
UShort_t window_center_time = time;

set<UShort_t> windowPeaksIndex;

// add the peak to the window
windowIndex[peakIndex] = windowCenter.size();
windowPeaksIndex.insert(peakIndex);

// add the peaks that are in the window
for (size_t otherPeakIndex = peakIndex + 1; otherPeakIndex < eventPeaks.size(); otherPeakIndex++) {
const auto& [otherChannelId, otherTime, otherAmplitude] = eventPeaks[otherPeakIndex];

if (otherTime < windowTime) {
if (otherTime < windowTimeStart) {
continue;
}

if (otherTime > windowEnd) {
if (otherTime > windowTimeEnd) {
break;
}

windowIndex[otherPeakIndex] = windowCenter.size();
windowPeaksIndex.insert(otherPeakIndex);
}

for (const auto& index : windowPeaksIndex) {
windowPeakIndex[index] = window_index;
windowPeakCenter[index] = window_center_time;
windowPeakMultiplicity[index] = windowPeaksIndex.size();
}

window_index++;
}

// subtract 1 from windowIndex so it starts on 0
for (auto& index : windowIndex) {
index--;
SetObservableValue("windowPeakIndex", windowPeakIndex);
SetObservableValue("windowPeakTimeBin", windowPeakCenter);
SetObservableValue("windowPeakMultiplicity", windowPeakMultiplicity);

vector<UShort_t> windowCenter;
vector<UShort_t> windowMultiplicity;

set<UShort_t> windowPeaksIndexInserted;
for (size_t peakIndex = 0; peakIndex < eventPeaks.size(); peakIndex++) {
const auto& window_peak_index = windowPeakIndex[peakIndex];
if (windowPeaksIndexInserted.find(window_peak_index) != windowPeaksIndexInserted.end()) {
continue;
}

windowPeaksIndexInserted.insert(window_peak_index);

windowCenter.push_back(windowPeakCenter[peakIndex]);
windowMultiplicity.push_back(windowPeakMultiplicity[peakIndex]);
}

// validation
// check only values from 0 ... windowCenter.size() -1 are in windowIndex
// ALL values in this range should appear at least once
for (size_t index = 0; index < windowCenter.size(); index++) {
bool found = false;
for (const auto& window : windowIndex) {
if (window == index) {
found = true;
break;
SetObservableValue("windowTimeBin", windowCenter);
SetObservableValue("windowMultiplicity", windowMultiplicity);

if (fTimeBinToTimeFactorMultiplier != 0.0) {
vector<Double_t> windowCenterTime(windowCenter.size(), 0.0);
vector<Double_t> windowPeakCenterTime(windowPeakCenter.size(), 0.0);
vector<Double_t> peaksTimePhysical(peaksTime.size(), 0.0);

// lambda to convert time bin to time
auto timeBinToTime = [this](UShort_t timeBin) {
if (!fTimeConversionElectronics) {
return fTimeBinToTimeFactorMultiplier * timeBin - fTimeBinToTimeFactorOffset;
} else {
// @jporron
double zero = -1 + (512 * fTimeBinToTimeFactorMultiplier - fTimeBinToTimeFactorOffset -
fTimeBinToTimeFactorOffsetTCM) /
fTimeBinToTimeFactorMultiplier;
return fTimeBinToTimeFactorMultiplier * (timeBin - zero);
}
};

for (size_t i = 0; i < windowCenter.size(); i++) {
windowCenterTime[i] = timeBinToTime(windowCenter[i]);
}
if (!found) {
cerr << "TRestRawPeaksFinderProcess::ProcessEvent: window index " << index << " not found"
<< endl;
exit(1);

for (size_t i = 0; i < windowPeakCenter.size(); i++) {
windowPeakCenterTime[i] = timeBinToTime(windowPeakCenter[i]);
}

for (size_t i = 0; i < peaksTime.size(); i++) {
peaksTimePhysical[i] = timeBinToTime(peaksTime[i]);
}

SetObservableValue("windowTime", windowCenterTime);
SetObservableValue("windowPeakTime", windowPeakCenterTime);
SetObservableValue("peaksTime", peaksTimePhysical);
}

SetObservableValue("windowIndex", windowIndex);
SetObservableValue("windowCenter", windowCenter);
if (fADCtoEnergyFactor != 0.0 || !fChannelIDToADCtoEnergyFactor.empty()) {
vector<Double_t> peaksEnergyPhysical(peaksAmplitude.size(), 0.0);
Double_t peaksEnergySum = 0.0;

if (fADCtoEnergyFactor != 0.0) {
// same factor for all peaks
for (size_t i = 0; i < peaksAmplitude.size(); i++) {
peaksEnergyPhysical[i] = peaksAmplitude[i] * fADCtoEnergyFactor;
peaksEnergySum += peaksEnergyPhysical[i];
}
} else {
// use map to get factor for each channel
for (size_t i = 0; i < peaksAmplitude.size(); i++) {
const auto channelId = peaksChannelId[i];
// check if the channel is in the map
if (fChannelIDToADCtoEnergyFactor.find(channelId) == fChannelIDToADCtoEnergyFactor.end()) {
cerr << "TRestRawPeaksFinderProcess::ProcessEvent: channel id " << channelId
<< " not found in the map of ADC to energy factors" << endl;
exit(1);
}
const auto factor = fChannelIDToADCtoEnergyFactor[channelId];

peaksEnergyPhysical[i] = peaksAmplitude[i] * factor;
peaksEnergySum += peaksEnergyPhysical[i];
}
}

// Remove peakless veto signals after the peak finding if chosen
if (fRemovePeaklessVetos && !fRemoveAllVetos) {
SetObservableValue("peaksEnergy", peaksEnergyPhysical);
SetObservableValue("peaksEnergySum", peaksEnergySum);
}

// Remove peak-less veto signals after the peak finding if chosen
if (fRemovePeaklessVetoes && !fRemoveAllVetoes) {
set<UShort_t> peakSignalIds;
for (const auto& [channelId, time, amplitude] : eventPeaks) {
peakSignalIds.insert(channelId);
Expand All @@ -179,7 +262,7 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) {
}

// Remove all veto signals after the peak finding if chosen
if (fRemoveAllVetos) {
if (fRemoveAllVetoes) {
vector<UShort_t> signalsToRemove;
for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) {
const auto signal = fInputEvent->GetSignal(signalIndex);
Expand All @@ -200,6 +283,34 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) {
return inputEvent;
}

std::map<UShort_t, double> parseStringToMap(const std::string& input) {
std::map<UShort_t, double> my_map;
std::string cleaned_input;

// Step 1: Clean the input to remove curly braces and spaces
for (char c : input) {
if (!std::isspace(c) && c != '{' && c != '}') {
cleaned_input += c;
}
}

// Step 2: Parse each key-value pair
std::stringstream ss(cleaned_input);
std::string pair;
while (std::getline(ss, pair, ',')) {
// Find the colon to split key and value
size_t colon_pos = pair.find(':');
if (colon_pos != std::string::npos) {
// Extract and convert key and value
int key = std::stoi(pair.substr(0, colon_pos));
double value = std::stod(pair.substr(colon_pos + 1));
my_map[key] = value;
}
}

return my_map;
}

void TRestRawPeaksFinderProcess::InitFromConfigFile() {
const auto filterType = GetParameter("channelType", "");
if (!filterType.empty()) {
Expand All @@ -214,8 +325,29 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() {
fBaselineRange = Get2DVectorParameterWithUnits("baselineRange", fBaselineRange);
fDistance = StringToDouble(GetParameter("distance", fDistance));
fWindow = StringToDouble(GetParameter("window", fWindow));
fRemoveAllVetos = StringToBool(GetParameter("removeAllVetos", fRemoveAllVetos));
fRemovePeaklessVetos = StringToBool(GetParameter("removePeaklessVetos", fRemovePeaklessVetos));
fRemoveAllVetoes = StringToBool(GetParameter("removeAllVetos", fRemoveAllVetoes));
fRemovePeaklessVetoes = StringToBool(GetParameter("removePeaklessVetos", fRemovePeaklessVetoes));

fTimeBinToTimeFactorMultiplier = GetDblParameterWithUnits("sampling", fTimeBinToTimeFactorMultiplier);
fTimeBinToTimeFactorOffset = GetDblParameterWithUnits("trigDelay", fTimeBinToTimeFactorOffset);
fTimeBinToTimeFactorOffsetTCM = GetDblParameterWithUnits("trigDelayTCM", fTimeBinToTimeFactorOffsetTCM);
fTimeConversionElectronics =
StringToBool(GetParameter("trigDelayElectronics", fTimeConversionElectronics));

fADCtoEnergyFactor = GetDblParameterWithUnits("adcToEnergyFactor", fADCtoEnergyFactor);
const string fChannelIDToADCtoEnergyFactorAsString = GetParameter("channelIDToADCtoEnergyFactor", "");
if (!fChannelIDToADCtoEnergyFactorAsString.empty()) {
// map should be in the format: "{channelId1: factor1, channelId2: factor2, ...}" (spaces are allowed
// but not required)
fChannelIDToADCtoEnergyFactor = parseStringToMap(fChannelIDToADCtoEnergyFactorAsString);
}

if (fADCtoEnergyFactor != 0 && !fChannelIDToADCtoEnergyFactor.empty()) {
cerr << "TRestRawPeaksFinderProcess::InitFromConfigFile: both adcToEnergyFactor and "
"channelIDToADCtoEnergyFactor are defined. Please, remove one of them."
<< endl;
exit(1);
}

if (fBaselineRange.X() > fBaselineRange.Y() || fBaselineRange.X() < 0 || fBaselineRange.Y() < 0) {
cerr << "TRestRawPeaksFinderProcess::InitProcess: baseline range is not sorted or < 0" << endl;
Expand All @@ -237,7 +369,7 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() {
exit(1);
}

if (filterType != "veto" && fRemovePeaklessVetos) {
if (filterType != "veto" && fRemovePeaklessVetoes) {
cerr << "TRestRawPeaksFinderProcess::InitProcess: removing veto signals only makes sense when the "
"process is applied to veto signals. Remove \"removePeaklessVetos\" parameter"
<< endl;
Expand Down

0 comments on commit f1b939d

Please sign in to comment.