Skip to content

Commit

Permalink
Use average BaseLine and BaseLineSigma and small changes in peak pos (#…
Browse files Browse the repository at this point in the history 
…143)

* Use average BaseLine and BaseLineSigma and small changes in peak position

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

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

* Solve warning that is being treated as error

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

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

* add check just in case

* typo in rml variable name

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

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

* emplace_back

* Set threshold variable of tpc as input

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

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

* bad initialization defaults

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Luis Antonio Obis Aparicio <[email protected]>
  • Loading branch information
@JPorron @pre-commit-ci @lobis
3 people authored Nov 11, 2024
1 parent f1b939d commit 2f69a48
Show file tree
Hide file tree
Showing 3 changed files with 122 additions and 45 deletions.
4 changes: 3 additions & 1 deletion inc/TRestRawPeaksFinderProcess.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -58,7 +60,7 @@ class TRestRawPeaksFinderProcess : public TRestEventProcess {
TRestRawPeaksFinderProcess() = default;
~TRestRawPeaksFinderProcess() = default;

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

#endif // REST_TRESTRAWPEAKSFINDERPROCESS_H
67 changes: 59 additions & 8 deletions src/TRestRawPeaksFinderProcess.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,39 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) {

vector<tuple<UShort_t, UShort_t, double>> 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();
Expand All @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand All @@ -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);
}
Expand All @@ -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;
Expand Down
96 changes: 60 additions & 36 deletions src/TRestRawSignal.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
}
}
}
}
Expand All @@ -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);
}
}
}
Expand All @@ -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;
}

Expand All @@ -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;
}

Expand All @@ -328,14 +347,15 @@ 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 "
"called first!"
<< endl;
fShowWarnings = false;
}
}
return fThresholdIntegral;
}

Expand Down Expand Up @@ -917,7 +937,9 @@ vector<pair<UShort_t, double>> 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<double> smoothedValues(numPoints, 0.0);
Expand Down Expand Up @@ -977,30 +999,32 @@ vector<pair<UShort_t, double>> 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<double>::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<UShort_t>(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);
}
}
}
Expand Down

0 comments on commit 2f69a48

Please sign in to comment.