diff --git a/inc/TRestRawPeaksFinderProcess.h b/inc/TRestRawPeaksFinderProcess.h index 85693f9..7e99b12 100644 --- a/inc/TRestRawPeaksFinderProcess.h +++ b/inc/TRestRawPeaksFinderProcess.h @@ -17,6 +17,8 @@ class TRestRawPeaksFinderProcess : public TRestEventProcess { /// \brief threshold over baseline to consider a peak Double_t fThresholdOverBaseline = 2.0; + /// \brief choose times the sigma of the baseline must be overcome to consider a peak + Double_t fSigmaOverBaseline = 10.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 (ADC units) @@ -58,7 +60,7 @@ class TRestRawPeaksFinderProcess : public TRestEventProcess { TRestRawPeaksFinderProcess() = default; ~TRestRawPeaksFinderProcess() = default; - ClassDefOverride(TRestRawPeaksFinderProcess, 5); + ClassDefOverride(TRestRawPeaksFinderProcess, 6); }; #endif // REST_TRESTRAWPEAKSFINDERPROCESS_H diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index 112fe1e..92c2107 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -30,6 +30,39 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { vector> eventPeaks; // signalId, time, amplitude + // Calculate average baseline and sigma of all the TPC signals + double BaseLineMean = 0.0; + double BaseLineSigmaMean = 0.0; + unsigned int countTPC = 0; + + for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) { + const auto signal = fInputEvent->GetSignal(signalIndex); + const UShort_t signalId = signal->GetSignalID(); + + const string channelType = fReadoutMetadata->GetTypeForChannelDaqId(signalId); + const string channelName = fReadoutMetadata->GetNameForChannelDaqId(signalId); + + // Check if channel type is in the list of selected channel types + if (fChannelTypes.find(channelType) == fChannelTypes.end()) { + continue; + } + + // Choose appropriate function based on channel type + if (channelType == "tpc") { + signal->CalculateBaseLine(fBaselineRange.X(), fBaselineRange.Y()); + // Accumulate the baseline and sigma values + BaseLineMean += signal->GetBaseLine(); + BaseLineSigmaMean += signal->GetBaseLineSigma(); + countTPC++; // Count the signals considered + } + } + + // Calculate the average if there were any matching signals + if (countTPC > 0) { + BaseLineMean /= countTPC; + BaseLineSigmaMean /= countTPC; + } + for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) { const auto signal = fInputEvent->GetSignal(signalIndex); const UShort_t signalId = signal->GetSignalID(); @@ -45,8 +78,18 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { // Choose appropriate function based on channel type if (channelType == "tpc") { signal->CalculateBaseLine(fBaselineRange.X(), fBaselineRange.Y()); - const auto peaks = - signal->GetPeaks(signal->GetBaseLine() + 5 * signal->GetBaseLineSigma(), fDistance); + + // I think count will never be 0, just in case + const double threshold = + (countTPC > 0) ? BaseLineMean + fSigmaOverBaseline * BaseLineSigmaMean + : signal->GetBaseLine() + fSigmaOverBaseline * signal->GetBaseLineSigma(); + if (countTPC <= 0) { + cerr << "TRestRawPeaksFinderProcess::ProcessEvent: TPC count is 0 in TPC loop, this should " + "not happen" + << endl; + exit(1); + } + const auto peaks = signal->GetPeaks(threshold, fDistance); for (const auto& [time, amplitude] : peaks) { eventPeaks.emplace_back(signalId, time, amplitude); @@ -321,12 +364,13 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() { // if no channel type is specified, use all channel types } - fThresholdOverBaseline = StringToDouble(GetParameter("thresholdOverBaseline", fThresholdOverBaseline)); + fThresholdOverBaseline = GetDblParameterWithUnits("thresholdOverBaseline", fThresholdOverBaseline); + fSigmaOverBaseline = GetDblParameterWithUnits("sigmaOverBaseline", fSigmaOverBaseline); fBaselineRange = Get2DVectorParameterWithUnits("baselineRange", fBaselineRange); - fDistance = StringToDouble(GetParameter("distance", fDistance)); - fWindow = StringToDouble(GetParameter("window", fWindow)); - fRemoveAllVetoes = StringToBool(GetParameter("removeAllVetos", fRemoveAllVetoes)); - fRemovePeaklessVetoes = StringToBool(GetParameter("removePeaklessVetos", fRemovePeaklessVetoes)); + fDistance = UShort_t(GetDblParameterWithUnits("distance", fDistance)); + fWindow = UShort_t(GetDblParameterWithUnits("window", fWindow)); + fRemoveAllVetoes = StringToBool(GetParameter("removeAllVetoes", fRemoveAllVetoes)); + fRemovePeaklessVetoes = StringToBool(GetParameter("removePeaklessVetoes", fRemovePeaklessVetoes)); fTimeBinToTimeFactorMultiplier = GetDblParameterWithUnits("sampling", fTimeBinToTimeFactorMultiplier); fTimeBinToTimeFactorOffset = GetDblParameterWithUnits("trigDelay", fTimeBinToTimeFactorOffset); @@ -359,6 +403,11 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() { exit(1); } + if (fSigmaOverBaseline < 0) { + cerr << "TRestRawPeaksFinderProcess::InitProcess: sigma over baseline is < 0" << endl; + exit(1); + } + if (fDistance <= 0) { cerr << "TRestRawPeaksFinderProcess::InitProcess: distance is < 0" << endl; exit(1); @@ -371,7 +420,7 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() { if (filterType != "veto" && fRemovePeaklessVetoes) { cerr << "TRestRawPeaksFinderProcess::InitProcess: removing veto signals only makes sense when the " - "process is applied to veto signals. Remove \"removePeaklessVetos\" parameter" + "process is applied to veto signals. Remove \"removePeaklessVetoes\" parameter" << endl; exit(1); } @@ -388,6 +437,8 @@ void TRestRawPeaksFinderProcess::PrintMetadata() { RESTMetadata << "Baseline range for tpc signals: " << fBaselineRange.X() << " - " << fBaselineRange.Y() << RESTendl; + RESTMetadata << "Sigma over baseline for tpc signals: " << fSigmaOverBaseline << RESTendl; + RESTMetadata << "Threshold over baseline for veto signals: " << fThresholdOverBaseline << RESTendl; RESTMetadata << "Distance: " << fDistance << RESTendl; diff --git a/src/TRestRawSignal.cxx b/src/TRestRawSignal.cxx index bf2c4a6..72747e3 100644 --- a/src/TRestRawSignal.cxx +++ b/src/TRestRawSignal.cxx @@ -264,12 +264,15 @@ void TRestRawSignal::InitializePointsOverThreshold(const TVector2& thrPar, Int_t if (pulse.size() >= (unsigned int)nPointsOver) { // auto stdev = TMath::StdDev(begin(pulse), end(pulse)); // calculate stdev - double mean = std::accumulate(pulse.begin(), pulse.end(), 0.0) / pulse.size(); + double mean = std::accumulate(pulse.begin(), pulse.end(), 0.0) / double(pulse.size()); double sq_sum = std::inner_product(pulse.begin(), pulse.end(), pulse.begin(), 0.0); - double stdev = std::sqrt(sq_sum / pulse.size() - mean * mean); + double stdev = std::sqrt(sq_sum / double(pulse.size()) - mean * mean); - if (stdev > signalTh * fBaseLineSigma) - for (int j = pos; j < i; j++) fPointsOverThreshold.push_back(j); + if (stdev > signalTh * fBaseLineSigma) { + for (int j = pos; j < i; j++) { + fPointsOverThreshold.push_back(j); + } + } } } } @@ -283,14 +286,18 @@ void TRestRawSignal::InitializePointsOverThreshold(const TVector2& thrPar, Int_t /// of fThresholdIntegral. This method is only used internally. /// void TRestRawSignal::CalculateThresholdIntegral() { - if (fRange.X() < 0) fRange.SetX(0); - if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) fRange.SetY(GetNumberOfPoints()); + if (fRange.X() < 0) { + fRange.SetX(0); + } + if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) { + fRange.SetY(GetNumberOfPoints()); + } fThresholdIntegral = 0; - for (unsigned int n = 0; n < fPointsOverThreshold.size(); n++) { - if (fPointsOverThreshold[n] >= fRange.X() && fPointsOverThreshold[n] < fRange.Y()) { - fThresholdIntegral += GetData(fPointsOverThreshold[n]); + for (int n : fPointsOverThreshold) { + if (n >= fRange.X() && n < fRange.Y()) { + fThresholdIntegral += GetData(n); } } } @@ -301,11 +308,17 @@ void TRestRawSignal::CalculateThresholdIntegral() { /// the integral is calculated in the full range. /// Double_t TRestRawSignal::GetIntegral() { - if (fRange.X() < 0) fRange.SetX(0); - if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) fRange.SetY(GetNumberOfPoints()); + if (fRange.X() < 0) { + fRange.SetX(0); + } + if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) { + fRange.SetY(GetNumberOfPoints()); + } Double_t sum = 0; - for (int i = fRange.X(); i < fRange.Y(); i++) sum += GetData(i); + for (int i = fRange.X(); i < fRange.Y(); i++) { + sum += GetData(i); + } return sum; } @@ -314,11 +327,17 @@ Double_t TRestRawSignal::GetIntegral() { /// by (startBin,endBin). /// Double_t TRestRawSignal::GetIntegralInRange(Int_t startBin, Int_t endBin) { - if (startBin < 0) startBin = 0; - if (endBin <= 0 || endBin > GetNumberOfPoints()) endBin = GetNumberOfPoints(); + if (startBin < 0) { + startBin = 0; + } + if (endBin <= 0 || endBin > GetNumberOfPoints()) { + endBin = GetNumberOfPoints(); + } Double_t sum = 0; - for (int i = startBin; i < endBin; i++) sum += GetRawData(i); + for (int i = startBin; i < endBin; i++) { + sum += GetRawData(i); + } return sum; } @@ -328,7 +347,7 @@ Double_t TRestRawSignal::GetIntegralInRange(Int_t startBin, Int_t endBin) { /// have been called first. /// Double_t TRestRawSignal::GetThresholdIntegral() { - if (fThresholdIntegral == -1) + if (fThresholdIntegral == -1) { if (fShowWarnings) { std::cout << "TRestRawSignal::GetThresholdIntegral. " "InitializePointsOverThreshold should be " @@ -336,6 +355,7 @@ Double_t TRestRawSignal::GetThresholdIntegral() { << endl; fShowWarnings = false; } + } return fThresholdIntegral; } @@ -917,7 +937,9 @@ vector> TRestRawSignal::GetPeaks(double threshold, UShort 10; // Region to compare for peak/no peak classification. 10 means 5 bins to each side const size_t numPoints = GetNumberOfPoints(); - if (numPoints == 0) return peaks; + if (numPoints == 0) { + return peaks; + } // Pre-calculate smoothed values for all bins using a rolling sum vector smoothedValues(numPoints, 0.0); @@ -977,30 +999,32 @@ vector> TRestRawSignal::GetPeaks(double threshold, UShort } // If it's a peak and it´s above the threshold and further than distance to the previous peak, add - // to peaks + // to peaks the biggest amplitude bin within the next "distance" bins and as amplitude the + // TripleMaxAverage. This is because for flat regions the detected peak is more to the left than + // the actual one. if (isPeak && smoothedValue > threshold) { if (peaks.empty() || i - peaks.back().first >= distance) { - double fitMinRange = i - 20; - double fitMaxRange = i + 20; - - // Create a Gaussian fit function - TF1 fitFunction("gaussianFit", "gaus", fitMinRange, fitMaxRange); - // Fit the data with the Gaussian function - fitFunction.SetRange(fitMinRange, fitMaxRange); // Initial parameters - - // Create histogram with the values to fit - TH1D histogram("hist", "hist", 40, fitMinRange, fitMaxRange); - for (int k = i - 20; k <= i + 20; ++k) { - histogram.SetBinContent(k - (i - 20) + 1, GetRawData(k)); // Set bin content + // Initialize variables to find the max amplitude within the next "distance" bins + int maxBin = i; + double maxAmplitude = smoothedValues[i]; + + // Look ahead within the specified distance to find the bin with the maximum amplitude + for (std::vector::size_type j = i + 1; + j <= i + distance && j < smoothedValues.size(); ++j) { + if (smoothedValues[j] > maxAmplitude) { + maxAmplitude = smoothedValues[j]; + maxBin = j; + } } - histogram.Fit(&fitFunction, "RQ"); - // Get peak position and amplitude from the fit - double peakPosition = fitFunction.GetParameter(1); - UShort_t formattedPeakPosition = static_cast(peakPosition); - double peakAmplitude = GetRawData(formattedPeakPosition); + // Calculate the peak amplitude as the average of maxBin and its two neighbors + double amplitude1 = GetRawData(maxBin - 1); + double amplitude2 = GetRawData(maxBin); + double amplitude3 = GetRawData(maxBin + 1); + double peakAmplitude = (amplitude1 + amplitude2 + amplitude3) / 3.0; - peaks.push_back(std::make_pair(formattedPeakPosition, peakAmplitude)); + // Store the peak position and amplitude + peaks.emplace_back(maxBin, peakAmplitude); } } }