diff --git a/inc/TRestRawPeaksFinderProcess.h b/inc/TRestRawPeaksFinderProcess.h index db59d91..85693f9 100644 --- a/inc/TRestRawPeaksFinderProcess.h +++ b/inc/TRestRawPeaksFinderProcess.h @@ -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 fChannelIDToADCtoEnergyFactor = {}; std::set fChannelTypes = {}; // this process will only be applied to selected channel types @@ -49,7 +58,7 @@ class TRestRawPeaksFinderProcess : public TRestEventProcess { TRestRawPeaksFinderProcess() = default; ~TRestRawPeaksFinderProcess() = default; - ClassDefOverride(TRestRawPeaksFinderProcess, 4); + ClassDefOverride(TRestRawPeaksFinderProcess, 5); }; #endif // REST_TRESTRAWPEAKSFINDERPROCESS_H diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index 6f16c43..112fe1e 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -28,7 +28,7 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { exit(1); } - vector> eventPeaks; + vector> eventPeaks; // signalId, time, amplitude for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) { const auto signal = fInputEvent->GetSignal(signalIndex); @@ -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 peaksChannelId; vector peaksTime; vector peaksAmplitude; - double peaksEnergy; - double amplitudeTotal = 0.0; + double peaksEnergy = 0.0; + UShort_t peaksCount = 0; + UShort_t peaksCountUnique = 0; + + set 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 windowIndex(eventPeaks.size(), 0); // Initialize with zeros - vector windowCenter; // for each different window, the center of the window + vector windowPeakIndex(eventPeaks.size(), std::numeric_limits::max()); + vector windowPeakCenter(eventPeaks.size(), 0); + vector 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::max()) { continue; } - // create a new window - windowCenter.push_back(time); + UShort_t window_center_time = time; + + set 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 windowCenter; + vector windowMultiplicity; + + set 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 windowCenterTime(windowCenter.size(), 0.0); + vector windowPeakCenterTime(windowPeakCenter.size(), 0.0); + vector 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 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 peakSignalIds; for (const auto& [channelId, time, amplitude] : eventPeaks) { peakSignalIds.insert(channelId); @@ -179,7 +262,7 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { } // Remove all veto signals after the peak finding if chosen - if (fRemoveAllVetos) { + if (fRemoveAllVetoes) { vector signalsToRemove; for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) { const auto signal = fInputEvent->GetSignal(signalIndex); @@ -200,6 +283,34 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { return inputEvent; } +std::map parseStringToMap(const std::string& input) { + std::map 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()) { @@ -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; @@ -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;