From 9a817d4295bb72bfc38aac0459ba4f868c023edc Mon Sep 17 00:00:00 2001 From: juanan Date: Thu, 16 Feb 2023 18:26:05 +0100 Subject: [PATCH 01/10] Implementing generic signal processing functions inside detectorlib --- inc/TRestDetectorSignal.h | 13 ++------ inc/TRestDetectorSignalEvent.h | 4 +-- src/TRestDetectorSignal.cxx | 52 ++++++-------------------------- src/TRestDetectorSignalEvent.cxx | 26 +++++++--------- 4 files changed, 25 insertions(+), 70 deletions(-) diff --git a/inc/TRestDetectorSignal.h b/inc/TRestDetectorSignal.h index e381e460..f5ab573b 100644 --- a/inc/TRestDetectorSignal.h +++ b/inc/TRestDetectorSignal.h @@ -43,13 +43,8 @@ class TRestDetectorSignal { void SetPoint(TVector2 p); public: -#ifndef __CINT__ - TGraph* fGraph; //! - std::vector fPointsOverThreshold; //! -#endif - // TODO other objects should probably skip using GetMaxIndex direclty Int_t GetMaxIndex(Int_t from = 0, Int_t to = 0); @@ -97,8 +92,6 @@ Double_t GetIntegralWithThreshold(Int_t from, Int_t to, Double_t baseline, Doubl void Normalize(Double_t scale = 1.); - std::vector GetPointsOverThreshold() { return fPointsOverThreshold; } - Double_t GetAverage(Int_t start = 0, Int_t end = 0); Int_t GetMaxPeakWidth(); Double_t GetMaxPeakWithTime(Double_t startTime, Double_t endTime); @@ -129,11 +122,9 @@ Double_t GetIntegralWithThreshold(Int_t from, Int_t to, Double_t baseline, Doubl void SetPoint(Double_t t, Double_t d); void SetPoint(Int_t index, Double_t t, Double_t d); - Double_t GetStandardDeviation(Int_t startBin, Int_t endBin); - Double_t GetBaseLine(Int_t startBin, Int_t endBin); - Double_t GetBaseLineSigma(Int_t startBin, Int_t endBin, Double_t baseline = 0); + void CalculateBaseLineAndSigma(Int_t startBin, Int_t endBin, Double_t& baseline, Double_t& baseLineSigma); - Double_t SubstractBaseline(Int_t startBin, Int_t endBin); + void SubstractBaseline(Int_t startBin, Int_t endBin); void AddOffset(Double_t offset); void MultiplySignalBy(Double_t factor); diff --git a/inc/TRestDetectorSignalEvent.h b/inc/TRestDetectorSignalEvent.h index f329e9ee..5069abb6 100644 --- a/inc/TRestDetectorSignalEvent.h +++ b/inc/TRestDetectorSignalEvent.h @@ -75,8 +75,8 @@ class TRestDetectorSignalEvent : public TRestEvent { Int_t GetSignalIndex(Int_t signalID); - Double_t GetBaseLineAverage(Int_t startBin, Int_t endBin); - Double_t GetBaseLineSigmaAverage(Int_t startBin, Int_t endBin); + void CalculateBaseLineAndSigmaAverage(Int_t startBin, Int_t endBin, Double_t& baseLineAvg, + Double_t& baseLineSigmaAvg); Double_t GetIntegral(Int_t startBin = 0, Int_t endBin = 0); Double_t GetMaxValue(); Double_t GetMinValue(); diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index 2fa54a53..abce7cf7 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -24,8 +24,8 @@ ///_______________________________________________________________________________ #include "TRestDetectorSignal.h" - #include "TFitResult.h" +#include "TRestSignalAnalysis.h" using namespace std; #include @@ -43,8 +43,6 @@ TRestDetectorSignal::TRestDetectorSignal() { fSignalID = -1; fSignalTime.clear(); fSignalCharge.clear(); - - fPointsOverThreshold.clear(); } TRestDetectorSignal::~TRestDetectorSignal() { @@ -531,52 +529,22 @@ void TRestDetectorSignal::GetSignalDelayed(TRestDetectorSignal* delayedSignal, I void TRestDetectorSignal::GetSignalSmoothed(TRestDetectorSignal* smthSignal, Int_t averagingPoints) { this->Sort(); - averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints - - Double_t sum = GetIntegral(0, averagingPoints); - for (int i = 0; i <= averagingPoints / 2; i++) smthSignal->AddPoint(GetTime(i), sum / averagingPoints); + auto smoothed = TRestSignalAnalysis::GetSignalSmoothed(fSignalCharge, averagingPoints); - for (int i = averagingPoints / 2 + 1; i < GetNumberOfPoints() - averagingPoints / 2; i++) { - sum -= this->GetData(i - (averagingPoints / 2 + 1)); - sum += this->GetData(i + averagingPoints / 2); - smthSignal->AddPoint(this->GetTime(i), sum / averagingPoints); - } - - for (int i = GetNumberOfPoints() - averagingPoints / 2; i < GetNumberOfPoints(); i++) - smthSignal->AddPoint(GetTime(i), sum / averagingPoints); + for (int i = 0; i < GetNumberOfPoints(); i++) smthSignal->AddPoint(GetTime(i), smoothed[i]); } -Double_t TRestDetectorSignal::GetBaseLine(Int_t startBin, Int_t endBin) { - if (endBin - startBin <= 0) return 0.; - - Double_t baseLine = 0; - for (int i = startBin; i < endBin; i++) baseLine += fSignalCharge[i]; - - return baseLine / (endBin - startBin); +void TRestDetectorSignal::CalculateBaseLineAndSigma(Int_t startBin, Int_t endBin, Double_t& baseLine, + Double_t& baseLineSigma) { + TRestSignalAnalysis::CalculateBaselineAndSigmaSD(fSignalCharge, startBin, endBin, baseLine, + baseLineSigma); } -Double_t TRestDetectorSignal::GetStandardDeviation(Int_t startBin, Int_t endBin) { - Double_t bL = GetBaseLine(startBin, endBin); - return GetBaseLineSigma(startBin, endBin, bL); -} - -Double_t TRestDetectorSignal::GetBaseLineSigma(Int_t startBin, Int_t endBin, Double_t baseline) { - Double_t bL = baseline; - if (bL == 0) bL = GetBaseLine(startBin, endBin); - - Double_t baseLineSigma = 0; - for (int i = startBin; i < endBin; i++) - baseLineSigma += (bL - fSignalCharge[i]) * (bL - fSignalCharge[i]); - - return TMath::Sqrt(baseLineSigma / (endBin - startBin)); -} - -Double_t TRestDetectorSignal::SubstractBaseline(Int_t startBin, Int_t endBin) { - Double_t bL = GetBaseLine(startBin, endBin); +void TRestDetectorSignal::SubstractBaseline(Int_t startBin, Int_t endBin) { + Double_t bL, bLS; + CalculateBaseLineAndSigma(startBin, endBin, bL, bLS); AddOffset(-bL); - - return bL; } void TRestDetectorSignal::AddOffset(Double_t offset) { diff --git a/src/TRestDetectorSignalEvent.cxx b/src/TRestDetectorSignalEvent.cxx index 46571715..6ef97ee1 100644 --- a/src/TRestDetectorSignalEvent.cxx +++ b/src/TRestDetectorSignalEvent.cxx @@ -91,26 +91,22 @@ Double_t TRestDetectorSignalEvent::GetIntegralWithThreshold(Int_t from, Int_t to } */ -Double_t TRestDetectorSignalEvent::GetBaseLineAverage(Int_t startBin, Int_t endBin) { - Double_t baseLineMean = 0; +void TRestDetectorSignalEvent::CalculateBaseLineAndSigmaAverage(Int_t startBin, Int_t endBin, + Double_t& baseLineAvg, + Double_t& baseLineSigmaAvg) { + baseLineAvg = 0; + baseLineSigmaAvg = 0; - for (int signal = 0; signal < GetNumberOfSignals(); signal++) { - Double_t baseline = GetSignal(signal)->GetBaseLine(startBin, endBin); - baseLineMean += baseline; - } - - return baseLineMean / GetNumberOfSignals(); -} - -Double_t TRestDetectorSignalEvent::GetBaseLineSigmaAverage(Int_t startBin, Int_t endBin) { - Double_t baseLineSigmaMean = 0; + Double_t bL, bLS; for (int signal = 0; signal < GetNumberOfSignals(); signal++) { - Double_t baselineSigma = GetSignal(signal)->GetBaseLineSigma(startBin, endBin); - baseLineSigmaMean += baselineSigma; + GetSignal(signal)->CalculateBaseLineAndSigma(startBin, endBin, bL, bLS); + baseLineAvg += bL; + baseLineSigmaAvg += bLS; } - return baseLineSigmaMean / GetNumberOfSignals(); + baseLineAvg /= GetNumberOfSignals(); + baseLineSigmaAvg /= GetNumberOfSignals(); } void TRestDetectorSignalEvent::SubstractBaselines(Int_t startBin, Int_t endBin) { From c4ae8747195d256fa3cc92f90406cf20196c2f1b Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 10:29:10 +0200 Subject: [PATCH 02/10] Implementing generic signal processing methods inside detectorlib --- inc/TRestDetectorSignal.h | 8 +- inc/TRestDetectorSignalToHitsProcess.h | 8 + src/TRestDetectorSignal.cxx | 343 +++-------------------- src/TRestDetectorSignalToHitsProcess.cxx | 149 +++------- 4 files changed, 95 insertions(+), 413 deletions(-) diff --git a/inc/TRestDetectorSignal.h b/inc/TRestDetectorSignal.h index c1e90cc4..7ea58dea 100644 --- a/inc/TRestDetectorSignal.h +++ b/inc/TRestDetectorSignal.h @@ -42,8 +42,7 @@ class TRestDetectorSignal { public: TGraph* fGraph; //! - // TODO other objects should probably skip using GetMaxIndex direclty - Int_t GetMaxIndex(Int_t from = 0, Int_t to = 0); + Int_t GetMaxIndex(); TVector2 GetMaxGauss(); TVector2 GetMaxLandau(); @@ -84,12 +83,11 @@ class TRestDetectorSignal { Double_t GetAverage(Int_t start = 0, Int_t end = 0); Int_t GetMaxPeakWidth(); - Double_t GetMaxPeakWithTime(Double_t startTime, Double_t endTime); Double_t GetMaxPeakValue(); Double_t GetMinPeakValue(); - Double_t GetMaxPeakTime(Int_t from = 0, Int_t to = 0); + Double_t GetMaxPeakTime(); Double_t GetMaxValue() { return GetMaxPeakValue(); } Double_t GetMinValue() { return GetMinPeakValue(); } @@ -105,9 +103,11 @@ class TRestDetectorSignal { void SetID(Int_t sID) { fSignalID = sID; } void NewPoint(Float_t time, Float_t data); + void IncreaseAmplitude(const TVector2& p); void IncreaseAmplitude(Double_t t, Double_t d); void SetPoint(Double_t t, Double_t d); + void SetPoint(const TVector2& p); void SetPoint(Int_t index, Double_t t, Double_t d); void CalculateBaseLineAndSigma(Int_t startBin, Int_t endBin, Double_t& baseline, Double_t& baseLineSigma); diff --git a/inc/TRestDetectorSignalToHitsProcess.h b/inc/TRestDetectorSignalToHitsProcess.h index 57a708f2..bdc39b74 100644 --- a/inc/TRestDetectorSignalToHitsProcess.h +++ b/inc/TRestDetectorSignalToHitsProcess.h @@ -77,6 +77,14 @@ class TRestDetectorSignalToHitsProcess : public TRestEventProcess { void LoadConfig(const std::string& configFilename, const std::string& name = ""); + /// Returns the Z coordinate from the hitTime, drifVelocity, fieldZDirection and zPositon (relative to the + /// readout) + inline Double_t GetHitZCoordinate(Double_t hitTime, Double_t driftVelocity, Double_t fieldZDirection, + Double_t zPosition) const { + const Double_t distanceToPlane = hitTime * fDriftVelocity; + return zPosition + fieldZDirection * distanceToPlane; + } + /// It prints out the process parameters stored in the metadata structure void PrintMetadata() override { BeginPrintProcess(); diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index 0dcc3519..e7a9c5ba 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -69,17 +69,15 @@ void TRestDetectorSignal::IncreaseAmplitude(Double_t t, Double_t d) { /// /// The input vector should contain a physical time and an amplitude. /// -void TRestDetectorSignal::IncreaseAmplitude(TVector2 p) { +void TRestDetectorSignal::IncreaseAmplitude(const TVector2& p) { Int_t index = GetTimeIndex(p.X()); - Float_t x = p.X(); - Float_t y = p.Y(); if (index >= 0) { - fSignalTime[index] = x; - fSignalCharge[index] += y; + fSignalTime[index] = p.X(); + fSignalCharge[index] += p.Y(); } else { - fSignalTime.push_back(x); - fSignalCharge.push_back(y); + fSignalTime.push_back(p.X()); + fSignalCharge.push_back(p.Y()); } } @@ -90,17 +88,15 @@ void TRestDetectorSignal::IncreaseAmplitude(TVector2 p) { /// /// The input vector should contain a physical time and an amplitude. /// -void TRestDetectorSignal::SetPoint(TVector2 p) { +void TRestDetectorSignal::SetPoint(const TVector2& p) { Int_t index = GetTimeIndex(p.X()); - Float_t x = p.X(); - Float_t y = p.Y(); if (index >= 0) { - fSignalTime[index] = x; - fSignalCharge[index] = y; + fSignalTime[index] = p.X(); + fSignalCharge[index] = p.Y(); } else { - fSignalTime.push_back(x); - fSignalCharge.push_back(y); + fSignalTime.push_back(p.X()); + fSignalCharge.push_back(p.Y()); } } @@ -129,10 +125,7 @@ Double_t TRestDetectorSignal::GetIntegral(Int_t startBin, Int_t endBin) { if (startBin < 0) startBin = 0; if (endBin <= 0 || endBin > GetNumberOfPoints()) endBin = GetNumberOfPoints(); - Double_t sum = 0; - for (int i = startBin; i < endBin; i++) sum += GetData(i); - - return sum; + return TRestSignalAnalysis::GetIntegral(fSignalCharge, startBin, endBin); } void TRestDetectorSignal::Normalize(Double_t scale) { @@ -142,22 +135,10 @@ void TRestDetectorSignal::Normalize(Double_t scale) { } Double_t TRestDetectorSignal::GetIntegralWithTime(Double_t startTime, Double_t endTime) { - Double_t sum = 0; - for (int i = 0; i < GetNumberOfPoints(); i++) - if (GetTime(i) >= startTime && GetTime(i) < endTime) sum += GetData(i); + int startBin = GetTimeIndex(startTime); + int endBin = GetTimeIndex(endTime); - return sum; -} - -Double_t TRestDetectorSignal::GetMaxPeakWithTime(Double_t startTime, Double_t endTime) { - Double_t max = -1E10; - - for (int i = 0; i < GetNumberOfPoints(); i++) - if (GetTime(i) >= startTime && GetTime(i) < endTime) { - if (this->GetData(i) > max) max = GetData(i); - } - - return max; + return GetIntegral(startBin, endBin); } /* {{{ @@ -223,118 +204,24 @@ Double_t TRestDetectorSignal::GetIntegralWithThreshold(Int_t from, Int_t to, Dou Double_t TRestDetectorSignal::GetAverage(Int_t start, Int_t end) { this->Sort(); - if (end == 0) end = this->GetNumberOfPoints(); + if (end <= 0) end = this->GetNumberOfPoints(); - Double_t sum = 0; - for (int i = start; i <= end; i++) { - sum += this->GetData(i); - } - return sum / (end - start + 1); + return TRestSignalAnalysis::GetAverage(fSignalCharge, start, end); } Int_t TRestDetectorSignal::GetMaxPeakWidth() { this->Sort(); - Int_t mIndex = this->GetMaxIndex(); - Double_t maxValue = this->GetData(mIndex); - - Double_t value = maxValue; - Int_t rightIndex = mIndex; - while (value > maxValue / 2) { - value = this->GetData(rightIndex); - rightIndex++; - } - Int_t leftIndex = mIndex; - value = maxValue; - while (value > maxValue / 2) { - value = this->GetData(leftIndex); - leftIndex--; - } - - return rightIndex - leftIndex; + return TRestSignalAnalysis::GetMaxPeakWidth(fSignalCharge); } Double_t TRestDetectorSignal::GetMaxPeakValue() { return GetData(GetMaxIndex()); } -Int_t TRestDetectorSignal::GetMaxIndex(Int_t from, Int_t to) { - Double_t max = -1E10; - Int_t index = 0; - - if (from < 0) from = 0; - if (to > GetNumberOfPoints()) to = GetNumberOfPoints(); - - if (to == 0) to = GetNumberOfPoints(); - - for (int i = from; i < to; i++) { - if (this->GetData(i) > max) { - max = GetData(i); - index = i; - } - } - - return index; -} - -// z position by gaussian fit - TVector2 TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the peak time in us and the energy { - Int_t maxRaw = GetMaxIndex(); // The bin where the maximum of the raw signal is found - Double_t maxRawTime = - GetTime(maxRaw); // The time of the bin where the maximum of the raw signal is found - Double_t energy = 0, time = 0; - Double_t lowerLimit = maxRawTime - 0.2; // us - Double_t upperLimit = maxRawTime + 0.4; // us - - TF1* gaus = new TF1("gaus", "gaus", lowerLimit, upperLimit); - TH1F* h1 = new TH1F("h1", "h1", 1000, 0, - 10); // Histogram to store the signal. For now the number of bins is fixed. - - // copying the signal peak to a histogram - for (int i = 0; i < GetNumberOfPoints(); i++) { - h1->Fill(GetTime(i), GetData(i)); - } - /* - TCanvas* c = new TCanvas("c", "Signal fit", 200, 10, 1280, 720); - h1->GetXaxis()->SetTitle("Time (us)"); - h1->GetYaxis()->SetTitle("Amplitude"); - h1->Draw(); - */ - - TFitResultPtr fitResult = - h1->Fit(gaus, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in the function range; S - // = save and return the fit result - - if (fitResult->IsValid()) { - energy = gaus->GetParameter(0); - time = gaus->GetParameter(1); - } else { - // the fit failed, return -1 to indicate failure - energy = -1; - time = -1; - cout << endl - << "WARNING: bad fit to signal with ID " << GetID() << " with maximum at time = " << maxRawTime - << " ns " - << "\n" - << "Failed fit parameters = " << gaus->GetParameter(0) << " || " << gaus->GetParameter(1) - << " || " << gaus->GetParameter(2) << "\n" - << "Assigned fit parameters : energy = " << energy << ", time = " << time << endl; - /* - TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720); - h1->Draw(); - c2->Update(); - getchar(); - delete c2; - */ - } - - TVector2 fitParam(time, energy); - - delete h1; - delete gaus; - - return fitParam; + auto gr = GetGraph(); + return TRestSignalAnalysis::GetMaxGauss(gr); } // z position by landau fit @@ -342,185 +229,50 @@ TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the pe TVector2 TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the peak time in us and the energy { - Int_t maxRaw = GetMaxIndex(); // The bin where the maximum of the raw signal is found - Double_t maxRawTime = - GetTime(maxRaw); // The time of the bin where the maximum of the raw signal is found - Double_t energy = 0, time = 0; - Double_t lowerLimit = maxRawTime - 0.2; // us - Double_t upperLimit = maxRawTime + 0.4; // us - - TF1* landau = new TF1("landau", "landau", lowerLimit, upperLimit); - TH1F* h1 = new TH1F("h1", "h1", 1000, 0, - 10); // Histogram to store the signal. For now the number of bins is fixed. - - // copying the signal peak to a histogram - for (int i = 0; i < GetNumberOfPoints(); i++) { - h1->Fill(GetTime(i), GetData(i)); - } - - TFitResultPtr fitResult = - h1->Fit(landau, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in the function range; - // S = save and return the fit result - if (fitResult->IsValid()) { - energy = landau->GetParameter(0); - time = landau->GetParameter(1); - } else { - // the fit failed, return -1 to indicate failure - energy = -1; - time = -1; - cout << endl - << "WARNING: bad fit to signal with ID " << GetID() << " with maximum at time = " << maxRawTime - << " ns " - << "\n" - << "Failed fit parameters = " << landau->GetParameter(0) << " || " << landau->GetParameter(1) - << " || " << landau->GetParameter(2) << "\n" - << "Assigned fit parameters : energy = " << energy << ", time = " << time << endl; - /* - TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720); - h1->Draw(); - c2->Update(); - getchar(); - delete c2; - */ - } - - TVector2 fitParam(time, energy); - - delete h1; - delete landau; - - return fitParam; -} - -// z position by aget fit - -Double_t agetResponseFunction(Double_t* x, Double_t* par) { // x contains as many elements as the number of - // dimensions (in this case 1, i.e. x[0]), and - // par contains as many elements as the number of free parameters in my function. - - Double_t arg = - (x[0] - par[1] + 1.1664) / - par[2]; // 1.1664 is the x value where the maximum of the base function (i.e. without parameters) is. - Double_t f = par[0] / 0.0440895 * exp(-3 * (arg)) * (arg) * (arg) * - (arg)*sin(arg); // to rescale the Y axis and get amplitude. - return f; + auto gr = GetGraph(); + return TRestSignalAnalysis::GetMaxLandau(gr); } TVector2 TRestDetectorSignal::GetMaxAget() // returns a 2vector with the time of the peak time in us and the energy { - Int_t maxRaw = GetMaxIndex(); // The bin where the maximum of the raw signal is found - Double_t maxRawTime = - GetTime(maxRaw); // The time of the bin where the maximum of the raw signal is found - Double_t energy = 0, time = 0; - // The intervals below are small because otherwise the function doesn't fit anymore. - Double_t lowerLimit = maxRawTime - 0.2; // us - Double_t upperLimit = maxRawTime + 0.7; // us - - TF1* aget = new TF1("aget", agetResponseFunction, lowerLimit, upperLimit, 3); // - TH1F* h1 = new TH1F("h1", "h1", 1000, 0, - 10); // Histogram to store the signal. For now the number of bins is fixed. - aget->SetParameters(500, maxRawTime, 1.2); - - // copying the signal peak to a histogram - for (int i = 0; i < GetNumberOfPoints(); i++) { - h1->Fill(GetTime(i), GetData(i)); - } - - TFitResultPtr fitResult = - h1->Fit(aget, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in - // the function range; S = save and return the fit result - - if (fitResult->IsValid()) { - energy = aget->GetParameter(0); - time = aget->GetParameter(1); - } else { - // the fit failed, return -1 to indicate failure - energy = -1; - time = -1; - cout << endl - << "WARNING: bad fit to signal with ID " << GetID() << " with maximum at time = " << maxRawTime - << " ns " - << "\n" - << "Failed fit parameters = " << aget->GetParameter(0) << " || " << aget->GetParameter(1) - << " || " << aget->GetParameter(2) << "\n" - << "Assigned fit parameters : energy = " << energy << ", time = " << time << endl; - /* - TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720); - h1->Draw(); - c2->Update(); - getchar(); - delete c2; - */ - } - - TVector2 fitParam(time, energy); - - delete h1; - delete aget; - - return fitParam; + auto gr = GetGraph(); + return TRestSignalAnalysis::GetMaxAget(gr); } -Double_t TRestDetectorSignal::GetMaxPeakTime(Int_t from, Int_t to) { return GetTime(GetMaxIndex(from, to)); } -Double_t TRestDetectorSignal::GetMinPeakValue() { return GetData(GetMinIndex()); } +Double_t TRestDetectorSignal::GetMaxPeakTime() { return GetTime(GetMaxIndex()); } -Int_t TRestDetectorSignal::GetMinIndex() { - Double_t min = 1E10; - Int_t index = 0; +Double_t TRestDetectorSignal::GetMinPeakValue() { return GetData(GetMinIndex()); } - for (int i = 0; i < GetNumberOfPoints(); i++) { - if (this->GetData(i) < min) { - min = GetData(i); - index = i; - } - } +Int_t TRestDetectorSignal::GetMaxIndex() { return TRestSignalAnalysis::GetMaxBin(fSignalCharge); } - return index; -} +Int_t TRestDetectorSignal::GetMinIndex() { return TRestSignalAnalysis::GetMinBin(fSignalCharge); } Double_t TRestDetectorSignal::GetMinTime() { - Double_t minTime = 1E10; - for (int n = 0; n < GetNumberOfPoints(); n++) - if (minTime > fSignalTime[n]) minTime = fSignalTime[n]; - - return minTime; + int index = TRestSignalAnalysis::GetMinBin(fSignalTime); + return GetTime(index); } Double_t TRestDetectorSignal::GetMaxTime() { - Double_t maxTime = -1E10; - for (int n = 0; n < GetNumberOfPoints(); n++) - if (maxTime < fSignalTime[n]) maxTime = fSignalTime[n]; - - return maxTime; + int index = TRestSignalAnalysis::GetMaxBin(fSignalTime); + return GetTime(index); } Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) { - Float_t time = t; + auto it = find(fSignalTime.begin(), fSignalTime.end(), t); + if (it != fSignalTime.end()) { + return *it; + } - for (int n = 0; n < GetNumberOfPoints(); n++) - if (time == fSignalTime[n]) return n; return -1; } -Bool_t TRestDetectorSignal::isSorted() { - for (int i = 0; i < GetNumberOfPoints() - 1; i++) { - if (GetTime(i + 1) < GetTime(i)) return false; - } - return true; -} +Bool_t TRestDetectorSignal::isSorted() { return is_sorted(fSignalTime.begin(), fSignalTime.end()); } void TRestDetectorSignal::Sort() { - while (!isSorted()) { - for (int i = 0; i < GetNumberOfPoints(); i++) { - for (int j = i; j < GetNumberOfPoints(); j++) { - if (GetTime(i) > GetTime(j)) { - iter_swap(fSignalTime.begin() + i, fSignalTime.begin() + j); - iter_swap(fSignalCharge.begin() + i, fSignalCharge.begin() + j); - } - } - } - } + sort(fSignalCharge.begin(), fSignalCharge.end(), + [&](size_t i, size_t j) { return fSignalTime[i] > fSignalTime[j]; }); + sort(fSignalTime.begin(), fSignalTime.end()); } void TRestDetectorSignal::GetDifferentialSignal(TRestDetectorSignal* diffSgnl, Int_t smearPoints) { @@ -554,7 +306,7 @@ void TRestDetectorSignal::GetSignalSmoothed(TRestDetectorSignal* smthSignal, Int auto smoothed = TRestSignalAnalysis::GetSignalSmoothed(fSignalCharge, averagingPoints); - for (int i = 0; i < GetNumberOfPoints(); i++) smthSignal->AddPoint(GetTime(i), smoothed[i]); + for (int i = 0; i < GetNumberOfPoints(); i++) smthSignal->IncreaseAmplitude(GetTime(i), smoothed[i]); } void TRestDetectorSignal::CalculateBaseLineAndSigma(Int_t startBin, Int_t endBin, Double_t& baseLine, @@ -624,11 +376,8 @@ void TRestDetectorSignal::GetWhiteNoiseSignal(TRestDetectorSignal* noiseSgnl, Do this->Sort(); for (int i = 0; i < GetNumberOfPoints(); i++) { - TRandom3* fRandom = new TRandom3(0); - - noiseSgnl->IncreaseAmplitude(GetTime(i), GetData(i) + fRandom->Gaus(0, noiseLevel)); - - delete fRandom; + TRandom3 fRandom(0); + noiseSgnl->IncreaseAmplitude(GetTime(i), GetData(i) + fRandom.Gaus(0, noiseLevel)); } } @@ -684,7 +433,7 @@ TGraph* TRestDetectorSignal::GetGraph(Int_t color) { fGraph = nullptr; } - fGraph = new TGraph(); + fGraph = new TGraph(GetNumberOfPoints(), fSignalTime.data(), fSignalCharge.data()); // cout << "Signal ID " << this->GetSignalID( ) << " points " << // this->GetNumberOfPoints() << endl; @@ -693,11 +442,5 @@ TGraph* TRestDetectorSignal::GetGraph(Int_t color) { fGraph->SetLineColor(color); fGraph->SetMarkerStyle(7); - int points = 0; - for (int n = 0; n < GetNumberOfPoints(); n++) { - fGraph->SetPoint(points, GetTime(n), GetData(n)); - points++; - } - return fGraph; } diff --git a/src/TRestDetectorSignalToHitsProcess.cxx b/src/TRestDetectorSignalToHitsProcess.cxx index 6e229cf3..59f23ed1 100644 --- a/src/TRestDetectorSignalToHitsProcess.cxx +++ b/src/TRestDetectorSignalToHitsProcess.cxx @@ -125,6 +125,8 @@ /// #include "TRestDetectorSignalToHitsProcess.h" +#include "TRestSignalAnalysis.h" + #include "TRestDetectorSetup.h" using namespace std; @@ -278,13 +280,7 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven if (fMethod == "onlyMax") { Double_t hitTime = sgnl->GetMaxPeakTime(); - Double_t distanceToPlane = hitTime * fDriftVelocity; - - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) - cout << "Distance to plane : " << distanceToPlane << endl; - - Double_t z = zPosition + fieldZDirection * distanceToPlane; - + Double_t z = GetHitZCoordinate(hitTime, fDriftVelocity, fieldZDirection, zPosition); Double_t energy = sgnl->GetMaxPeakValue(); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) @@ -293,76 +289,29 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven fHitsEvent->AddHit(x, y, z, energy, 0, type); } else if (fMethod == "tripleMax") { - Int_t bin = sgnl->GetMaxIndex(); - int binprev = (bin - 1) < 0 ? bin : bin - 1; - int binnext = (bin + 1) > sgnl->GetNumberOfPoints() - 1 ? bin : bin + 1; - - Double_t hitTime = sgnl->GetTime(bin); - Double_t energy = sgnl->GetData(bin); - - Double_t distanceToPlane = hitTime * fDriftVelocity; - Double_t z = zPosition + fieldZDirection * distanceToPlane; - - fHitsEvent->AddHit(x, y, z, energy, 0, type); + auto gr = sgnl->GetGraph(); + auto tripleMax = TRestSignalAnalysis::GetTripleMax(gr); - hitTime = sgnl->GetTime(binprev); - energy = sgnl->GetData(binprev); - - distanceToPlane = hitTime * fDriftVelocity; - z = zPosition + fieldZDirection * distanceToPlane; - - fHitsEvent->AddHit(x, y, z, energy, 0, type); - - hitTime = sgnl->GetTime(binnext); - energy = sgnl->GetData(binnext); - - distanceToPlane = hitTime * fDriftVelocity; - z = zPosition + fieldZDirection * distanceToPlane; - - fHitsEvent->AddHit(x, y, z, energy, 0, type); + for (const auto& [hitTime, energy] : tripleMax) { + const double z = GetHitZCoordinate(hitTime, fDriftVelocity, fieldZDirection, zPosition); + fHitsEvent->AddHit(x, y, z, energy, 0, type); - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { - cout << "Distance to plane : " << distanceToPlane << endl; - cout << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << z - << " Energy : " << energy << endl; + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { + cout << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << z + << " Energy : " << energy << endl; + } } - } else if (fMethod == "tripleMaxAverage") { - Int_t bin = sgnl->GetMaxIndex(); - int binprev = (bin - 1) < 0 ? bin : bin - 1; - int binnext = (bin + 1) > sgnl->GetNumberOfPoints() - 1 ? bin : bin + 1; - - Double_t hitTime = sgnl->GetTime(bin); - Double_t energy1 = sgnl->GetData(bin); - - Double_t distanceToPlane = hitTime * fDriftVelocity; - Double_t z1 = zPosition + fieldZDirection * distanceToPlane; - - hitTime = sgnl->GetTime(binprev); - Double_t energy2 = sgnl->GetData(binprev); - - distanceToPlane = hitTime * fDriftVelocity; - Double_t z2 = zPosition + fieldZDirection * distanceToPlane; - hitTime = sgnl->GetTime(binnext); - Double_t energy3 = sgnl->GetData(binnext); - - distanceToPlane = hitTime * fDriftVelocity; - Double_t z3 = zPosition + fieldZDirection * distanceToPlane; - - Double_t eTot = energy1 + energy2 + energy3; - - Double_t zAvg = ((z1 * energy1) + (z2 * energy2) + (z3 * energy3)) / eTot; - // Double_t zAvg = (z1 + z2 + z3) / 3.0; - Double_t eAvg = eTot / 3.0; + } else if (fMethod == "tripleMaxAverage") { + auto gr = sgnl->GetGraph(); + const TVector2 tripleMaxAvg = TRestSignalAnalysis::GetTripleMaxAverage(gr); + const double z = GetHitZCoordinate(tripleMaxAvg.X(), fDriftVelocity, fieldZDirection, zPosition); - fHitsEvent->AddHit(x, y, zAvg, eAvg, 0, type); + fHitsEvent->AddHit(x, y, z, tripleMaxAvg.Y(), 0, type); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { - cout << "Distance to plane: " << distanceToPlane << endl; - cout << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << zAvg - << " Energy : " << eAvg << endl; - cout << "z1, z2, z3 = " << z1 << ", " << z2 << ", " << z3 << endl; - cout << "E1, E2, E3 = " << energy1 << ", " << energy2 << ", " << energy3 << endl; + cout << "Adding hit. Time : " << tripleMaxAvg.X() << " x : " << x << " y : " << y + << " z : " << z << " Energy : " << tripleMaxAvg.Y() << endl; } } else if (fMethod == "gaussFit") { @@ -372,8 +321,7 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven Double_t z = -1.0; if (gaussFit.X() >= 0.0) { hitTime = gaussFit.X(); - Double_t distanceToPlane = hitTime * fDriftVelocity; - z = zPosition + fieldZDirection * distanceToPlane; + z = GetHitZCoordinate(hitTime, fDriftVelocity, fieldZDirection, zPosition); } Double_t energy = -1.0; if (gaussFit.Y() >= 0.0) { @@ -400,8 +348,7 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven Double_t z = -1.0; if (landauFit.X() >= 0.0) { hitTime = landauFit.X(); - Double_t distanceToPlane = hitTime * fDriftVelocity; - z = zPosition + fieldZDirection * distanceToPlane; + z = GetHitZCoordinate(hitTime, fDriftVelocity, fieldZDirection, zPosition); } Double_t energy = -1.0; if (landauFit.Y() >= 0.0) { @@ -428,8 +375,7 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven Double_t z = -1.0; if (agetFit.X() >= 0.0) { hitTime = agetFit.X(); - Double_t distanceToPlane = hitTime * fDriftVelocity; - z = zPosition + fieldZDirection * distanceToPlane; + z = GetHitZCoordinate(hitTime, fDriftVelocity, fieldZDirection, zPosition); } Double_t energy = -1.0; if (agetFit.Y() >= 0.0) { @@ -450,30 +396,29 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven fHitsEvent->AddHit(x, y, z, energy, 0, type); } else if (fMethod == "qCenter") { - Double_t energy_signal = 0; - Double_t distanceToPlane = 0; + Double_t energy = 0; + Double_t hitTime = 0; for (int j = 0; j < sgnl->GetNumberOfPoints(); j++) { - Double_t energy_point = sgnl->GetData(j); - energy_signal += energy_point; - distanceToPlane += sgnl->GetTime(j) * fDriftVelocity * energy_point; + energy += sgnl->GetData(j); + hitTime += sgnl->GetTime(j) * sgnl->GetData(j); } - Double_t energy = energy_signal / sgnl->GetNumberOfPoints(); - Double_t z = zPosition + fieldZDirection * (distanceToPlane / energy_signal); + if (energy != 0) hitTime /= energy; + energy /= (double)sgnl->GetNumberOfPoints(); + + const Double_t z = GetHitZCoordinate(hitTime, fDriftVelocity, fieldZDirection, zPosition); fHitsEvent->AddHit(x, y, z, energy, 0, type); } else if (fMethod == "all") { for (int j = 0; j < sgnl->GetNumberOfPoints(); j++) { - Double_t energy = sgnl->GetData(j); - - Double_t distanceToPlane = sgnl->GetTime(j) * fDriftVelocity; + const Double_t energy = sgnl->GetData(j); + const Double_t hitTime = sgnl->GetTime(j); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { cout << "Time : " << sgnl->GetTime(j) << " Drift velocity : " << fDriftVelocity << endl; - cout << "Distance to plane : " << distanceToPlane << endl; } - Double_t z = zPosition + fieldZDirection * distanceToPlane; + Double_t z = GetHitZCoordinate(hitTime, fDriftVelocity, fieldZDirection, zPosition); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) cout << "Adding hit. Time : " << sgnl->GetTime(j) << " x : " << x << " y : " << y @@ -482,30 +427,16 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven fHitsEvent->AddHit(x, y, z, energy, 0, type); } } else if (fMethod == "intwindow") { - Int_t nPoints = sgnl->GetNumberOfPoints(); - std::map > windowMap; - for (int j = 0; j < nPoints; j++) { - int index = sgnl->GetTime(j) / fIntWindow; - auto it = windowMap.find(index); - if (it != windowMap.end()) { - it->second.first++; - it->second.second += sgnl->GetData(j); - } else { - windowMap[index] = std::make_pair(1, sgnl->GetData(j)); - } - } + auto gr = sgnl->GetGraph(); + auto intWindow = TRestSignalAnalysis::GetIntWindow(gr, fIntWindow); - for (const auto& [index, pair] : windowMap) { - Double_t hitTime = index * fIntWindow + fIntWindow / 2.; - Double_t energy = pair.second / pair.first; + for (const auto& [hitTime, energy] : intWindow) { if (energy < fThreshold) continue; - RESTDebug << "TimeBin " << index << " Time " << hitTime << " Charge: " << energy - << " Thr: " << (fThreshold) << RESTendl; - Double_t distanceToPlane = hitTime * fDriftVelocity; - Double_t z = zPosition + fieldZDirection * distanceToPlane; + RESTDebug << " Time " << hitTime << " Charge: " << energy << " Thr: " << (fThreshold) + << RESTendl; + Double_t z = GetHitZCoordinate(hitTime, fDriftVelocity, fieldZDirection, zPosition); - RESTDebug << "Time : " << hitTime << " Drift velocity : " << fDriftVelocity - << "\nDistance to plane : " << distanceToPlane << RESTendl; + RESTDebug << "Time : " << hitTime << " Drift velocity : " << fDriftVelocity << RESTendl; RESTDebug << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << z << " type " << type << RESTendl; From dcf6d3b76ed8a34a9fef2f51e515bbe189e8bcfd Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 31 Mar 2023 08:42:19 +0000 Subject: [PATCH 03/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/TRestDetectorSignal.cxx | 1 + src/TRestDetectorSignalToHitsProcess.cxx | 3 +-- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index e7a9c5ba..ab5ecd26 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -24,6 +24,7 @@ ///_______________________________________________________________________________ #include "TRestDetectorSignal.h" + #include "TFitResult.h" #include "TRestSignalAnalysis.h" using namespace std; diff --git a/src/TRestDetectorSignalToHitsProcess.cxx b/src/TRestDetectorSignalToHitsProcess.cxx index 59f23ed1..c72f51ea 100644 --- a/src/TRestDetectorSignalToHitsProcess.cxx +++ b/src/TRestDetectorSignalToHitsProcess.cxx @@ -125,9 +125,8 @@ /// #include "TRestDetectorSignalToHitsProcess.h" -#include "TRestSignalAnalysis.h" - #include "TRestDetectorSetup.h" +#include "TRestSignalAnalysis.h" using namespace std; From fc304681aca69a403a11a87249db246f440afe75 Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 13:49:02 +0200 Subject: [PATCH 04/10] Adding protection at TRestDetectorSignal::Sort() in case both vectors have different size --- src/TRestDetectorSignal.cxx | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index ab5ecd26..efe9004a 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -271,6 +271,11 @@ Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) { Bool_t TRestDetectorSignal::isSorted() { return is_sorted(fSignalTime.begin(), fSignalTime.end()); } void TRestDetectorSignal::Sort() { + if (fSignalCharge.size() != fSignalTime.size()) { + cout << "Time and charge vectors have different size " << fSignalCharge.size() << " " + << fSignalTime.size() << endl; + exit(0); + } sort(fSignalCharge.begin(), fSignalCharge.end(), [&](size_t i, size_t j) { return fSignalTime[i] > fSignalTime[j]; }); sort(fSignalTime.begin(), fSignalTime.end()); From 9c910cf2aa3c3d8dbaa0d529527c5b7de4d1cf0b Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 15:54:04 +0200 Subject: [PATCH 05/10] Reverting Sort function since it is buggy in new implementation --- src/TRestDetectorSignal.cxx | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index efe9004a..72829593 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -271,14 +271,16 @@ Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) { Bool_t TRestDetectorSignal::isSorted() { return is_sorted(fSignalTime.begin(), fSignalTime.end()); } void TRestDetectorSignal::Sort() { - if (fSignalCharge.size() != fSignalTime.size()) { - cout << "Time and charge vectors have different size " << fSignalCharge.size() << " " - << fSignalTime.size() << endl; - exit(0); + while (!isSorted()) { + for (int i = 0; i < GetNumberOfPoints(); i++) { + for (int j = i; j < GetNumberOfPoints(); j++) { + if (GetTime(i) > GetTime(j)) { + iter_swap(fSignalTime.begin() + i, fSignalTime.begin() + j); + iter_swap(fSignalCharge.begin() + i, fSignalCharge.begin() + j); + } + } + } } - sort(fSignalCharge.begin(), fSignalCharge.end(), - [&](size_t i, size_t j) { return fSignalTime[i] > fSignalTime[j]; }); - sort(fSignalTime.begin(), fSignalTime.end()); } void TRestDetectorSignal::GetDifferentialSignal(TRestDetectorSignal* diffSgnl, Int_t smearPoints) { From 42bbebfefd065961f7ef14619e93c5452ac56430 Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 18:22:13 +0200 Subject: [PATCH 06/10] Addressing bug which made the pipeline fail --- src/TRestDetectorSignal.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index 72829593..75e64d53 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -262,7 +262,7 @@ Double_t TRestDetectorSignal::GetMaxTime() { Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) { auto it = find(fSignalTime.begin(), fSignalTime.end(), t); if (it != fSignalTime.end()) { - return *it; + return distance(fSignalTime.begin(), it); } return -1; From 5156e0e9912678bc068b41da0096c45ae53c9758 Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 20:21:50 +0200 Subject: [PATCH 07/10] Addressing bugs --- inc/TRestDetectorSignalToHitsProcess.h | 2 +- src/TRestDetectorSignal.cxx | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/inc/TRestDetectorSignalToHitsProcess.h b/inc/TRestDetectorSignalToHitsProcess.h index bdc39b74..90309cd1 100644 --- a/inc/TRestDetectorSignalToHitsProcess.h +++ b/inc/TRestDetectorSignalToHitsProcess.h @@ -81,7 +81,7 @@ class TRestDetectorSignalToHitsProcess : public TRestEventProcess { /// readout) inline Double_t GetHitZCoordinate(Double_t hitTime, Double_t driftVelocity, Double_t fieldZDirection, Double_t zPosition) const { - const Double_t distanceToPlane = hitTime * fDriftVelocity; + const Double_t distanceToPlane = hitTime * driftVelocity; return zPosition + fieldZDirection * distanceToPlane; } diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index 75e64d53..0b9de4e9 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -136,10 +136,11 @@ void TRestDetectorSignal::Normalize(Double_t scale) { } Double_t TRestDetectorSignal::GetIntegralWithTime(Double_t startTime, Double_t endTime) { - int startBin = GetTimeIndex(startTime); - int endBin = GetTimeIndex(endTime); + Double_t sum = 0; + for (int i = 0; i < GetNumberOfPoints(); i++) + if (GetTime(i) >= startTime && GetTime(i) < endTime) sum += GetData(i); - return GetIntegral(startBin, endBin); + return sum; } /* {{{ From c6de5fc80945c8e6b79feef73f7f9468774bae84 Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 23:37:11 +0200 Subject: [PATCH 08/10] Addressing pipeline failure, although the functions should be doing the same --- src/TRestDetectorSignal.cxx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index 0b9de4e9..6f085f65 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -261,11 +261,10 @@ Double_t TRestDetectorSignal::GetMaxTime() { } Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) { - auto it = find(fSignalTime.begin(), fSignalTime.end(), t); - if (it != fSignalTime.end()) { - return distance(fSignalTime.begin(), it); - } + Float_t time = t; + for (int n = 0; n < GetNumberOfPoints(); n++) + if (time == fSignalTime[n]) return n; return -1; } From dd34daed22263ec7d8bbb3766ae45f845453ff92 Mon Sep 17 00:00:00 2001 From: juanan Date: Sat, 1 Apr 2023 17:51:52 +0200 Subject: [PATCH 09/10] Removing baseline related methods --- inc/TRestDetectorSignal.h | 3 -- src/TRestDetectorSignal.cxx | 73 ------------------------------------- 2 files changed, 76 deletions(-) diff --git a/inc/TRestDetectorSignal.h b/inc/TRestDetectorSignal.h index 7ea58dea..99a69867 100644 --- a/inc/TRestDetectorSignal.h +++ b/inc/TRestDetectorSignal.h @@ -110,9 +110,6 @@ class TRestDetectorSignal { void SetPoint(const TVector2& p); void SetPoint(Int_t index, Double_t t, Double_t d); - void CalculateBaseLineAndSigma(Int_t startBin, Int_t endBin, Double_t& baseline, Double_t& baseLineSigma); - - void SubstractBaseline(Int_t startBin, Int_t endBin); void AddOffset(Double_t offset); void MultiplySignalBy(Double_t factor); diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index 6f085f65..9285c195 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -143,66 +143,6 @@ Double_t TRestDetectorSignal::GetIntegralWithTime(Double_t startTime, Double_t e return sum; } -/* {{{ -Double_t TRestDetectorSignal::GetIntegralWithThreshold(Int_t from, Int_t to, Int_t startBaseline, Int_t -endBaseline, - Double_t nSigmas, Int_t nPointsOverThreshold, - Double_t nMinSigmas) { - if (startBaseline < 0) startBaseline = 0; - if (endBaseline <= 0 || endBaseline > GetNumberOfPoints()) endBaseline = GetNumberOfPoints(); - - Double_t baseLine = GetBaseLine(startBaseline, endBaseline); - - Double_t pointThreshold = nSigmas * GetBaseLineSigma(startBaseline, endBaseline); - Double_t signalThreshold = nMinSigmas * GetBaseLineSigma(startBaseline, endBaseline); - - return GetIntegralWithThreshold(from, to, baseLine, pointThreshold, nPointsOverThreshold, - signalThreshold); -} - -Double_t TRestDetectorSignal::GetIntegralWithThreshold(Int_t from, Int_t to, Double_t baseline, - Double_t pointThreshold, Int_t nPointsOverThreshold, - Double_t signalThreshold) { - Double_t sum = 0; - Int_t nPoints = 0; - fPointsOverThreshold.clear(); - - if (to > GetNumberOfPoints()) to = GetNumberOfPoints(); - - Float_t maxValue = 0; - for (int i = from; i < to; i++) { - if (GetData(i) > baseline + pointThreshold) { - if (GetData(i) > maxValue) maxValue = GetData(i); - nPoints++; - } else { - if (nPoints >= nPointsOverThreshold) { - Double_t sig = GetStandardDeviation(i - nPoints, i); - if (sig > signalThreshold) { - for (int j = i - nPoints; j < i; j++) { - sum += this->GetData(j); - fPointsOverThreshold.push_back(j); - } - } - } - nPoints = 0; - maxValue = 0; - } - } - - if (nPoints >= nPointsOverThreshold) { - Double_t sig = GetStandardDeviation(to - nPoints, to); - if (sig > signalThreshold) { - for (int j = to - nPoints; j < to; j++) { - sum += this->GetData(j); - fPointsOverThreshold.push_back(j); - } - } - } - - return sum; -} -}}} */ - Double_t TRestDetectorSignal::GetAverage(Int_t start, Int_t end) { this->Sort(); @@ -317,19 +257,6 @@ void TRestDetectorSignal::GetSignalSmoothed(TRestDetectorSignal* smthSignal, Int for (int i = 0; i < GetNumberOfPoints(); i++) smthSignal->IncreaseAmplitude(GetTime(i), smoothed[i]); } -void TRestDetectorSignal::CalculateBaseLineAndSigma(Int_t startBin, Int_t endBin, Double_t& baseLine, - Double_t& baseLineSigma) { - TRestSignalAnalysis::CalculateBaselineAndSigmaSD(fSignalCharge, startBin, endBin, baseLine, - baseLineSigma); -} - -void TRestDetectorSignal::SubstractBaseline(Int_t startBin, Int_t endBin) { - Double_t bL, bLS; - CalculateBaseLineAndSigma(startBin, endBin, bL, bLS); - - AddOffset(-bL); -} - void TRestDetectorSignal::AddOffset(Double_t offset) { for (int i = 0; i < GetNumberOfPoints(); i++) fSignalCharge[i] = fSignalCharge[i] + offset; } From 87f846075dd6763912dae34e9c7cff5483b6454f Mon Sep 17 00:00:00 2001 From: juanan Date: Mon, 3 Apr 2023 20:29:16 +0200 Subject: [PATCH 10/10] Renaming TRestSignalAnalysis to TRestPulseShapeAnalysis --- src/TRestDetectorSignal.cxx | 24 ++++++++++++------------ src/TRestDetectorSignalToHitsProcess.cxx | 8 ++++---- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index 9285c195..129036cc 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -26,7 +26,7 @@ #include "TRestDetectorSignal.h" #include "TFitResult.h" -#include "TRestSignalAnalysis.h" +#include "TRestPulseShapeAnalysis.h" using namespace std; #include @@ -126,7 +126,7 @@ Double_t TRestDetectorSignal::GetIntegral(Int_t startBin, Int_t endBin) { if (startBin < 0) startBin = 0; if (endBin <= 0 || endBin > GetNumberOfPoints()) endBin = GetNumberOfPoints(); - return TRestSignalAnalysis::GetIntegral(fSignalCharge, startBin, endBin); + return TRestPulseShapeAnalysis::GetIntegral(fSignalCharge, startBin, endBin); } void TRestDetectorSignal::Normalize(Double_t scale) { @@ -148,13 +148,13 @@ Double_t TRestDetectorSignal::GetAverage(Int_t start, Int_t end) { if (end <= 0) end = this->GetNumberOfPoints(); - return TRestSignalAnalysis::GetAverage(fSignalCharge, start, end); + return TRestPulseShapeAnalysis::GetAverage(fSignalCharge, start, end); } Int_t TRestDetectorSignal::GetMaxPeakWidth() { this->Sort(); - return TRestSignalAnalysis::GetMaxPeakWidth(fSignalCharge); + return TRestPulseShapeAnalysis::GetMaxPeakWidth(fSignalCharge); } Double_t TRestDetectorSignal::GetMaxPeakValue() { return GetData(GetMaxIndex()); } @@ -163,7 +163,7 @@ TVector2 TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the peak time in us and the energy { auto gr = GetGraph(); - return TRestSignalAnalysis::GetMaxGauss(gr); + return TRestPulseShapeAnalysis::GetMaxGauss(gr); } // z position by landau fit @@ -172,31 +172,31 @@ TVector2 TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the peak time in us and the energy { auto gr = GetGraph(); - return TRestSignalAnalysis::GetMaxLandau(gr); + return TRestPulseShapeAnalysis::GetMaxLandau(gr); } TVector2 TRestDetectorSignal::GetMaxAget() // returns a 2vector with the time of the peak time in us and the energy { auto gr = GetGraph(); - return TRestSignalAnalysis::GetMaxAget(gr); + return TRestPulseShapeAnalysis::GetMaxAget(gr); } Double_t TRestDetectorSignal::GetMaxPeakTime() { return GetTime(GetMaxIndex()); } Double_t TRestDetectorSignal::GetMinPeakValue() { return GetData(GetMinIndex()); } -Int_t TRestDetectorSignal::GetMaxIndex() { return TRestSignalAnalysis::GetMaxBin(fSignalCharge); } +Int_t TRestDetectorSignal::GetMaxIndex() { return TRestPulseShapeAnalysis::GetMaxBin(fSignalCharge); } -Int_t TRestDetectorSignal::GetMinIndex() { return TRestSignalAnalysis::GetMinBin(fSignalCharge); } +Int_t TRestDetectorSignal::GetMinIndex() { return TRestPulseShapeAnalysis::GetMinBin(fSignalCharge); } Double_t TRestDetectorSignal::GetMinTime() { - int index = TRestSignalAnalysis::GetMinBin(fSignalTime); + int index = TRestPulseShapeAnalysis::GetMinBin(fSignalTime); return GetTime(index); } Double_t TRestDetectorSignal::GetMaxTime() { - int index = TRestSignalAnalysis::GetMaxBin(fSignalTime); + int index = TRestPulseShapeAnalysis::GetMaxBin(fSignalTime); return GetTime(index); } @@ -252,7 +252,7 @@ void TRestDetectorSignal::GetSignalDelayed(TRestDetectorSignal* delayedSignal, I void TRestDetectorSignal::GetSignalSmoothed(TRestDetectorSignal* smthSignal, Int_t averagingPoints) { this->Sort(); - auto smoothed = TRestSignalAnalysis::GetSignalSmoothed(fSignalCharge, averagingPoints); + auto smoothed = TRestPulseShapeAnalysis::GetSignalSmoothed(fSignalCharge, averagingPoints); for (int i = 0; i < GetNumberOfPoints(); i++) smthSignal->IncreaseAmplitude(GetTime(i), smoothed[i]); } diff --git a/src/TRestDetectorSignalToHitsProcess.cxx b/src/TRestDetectorSignalToHitsProcess.cxx index c72f51ea..071ffbc7 100644 --- a/src/TRestDetectorSignalToHitsProcess.cxx +++ b/src/TRestDetectorSignalToHitsProcess.cxx @@ -126,7 +126,7 @@ #include "TRestDetectorSignalToHitsProcess.h" #include "TRestDetectorSetup.h" -#include "TRestSignalAnalysis.h" +#include "TRestPulseShapeAnalysis.h" using namespace std; @@ -289,7 +289,7 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven fHitsEvent->AddHit(x, y, z, energy, 0, type); } else if (fMethod == "tripleMax") { auto gr = sgnl->GetGraph(); - auto tripleMax = TRestSignalAnalysis::GetTripleMax(gr); + auto tripleMax = TRestPulseShapeAnalysis::GetTripleMax(gr); for (const auto& [hitTime, energy] : tripleMax) { const double z = GetHitZCoordinate(hitTime, fDriftVelocity, fieldZDirection, zPosition); @@ -303,7 +303,7 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven } else if (fMethod == "tripleMaxAverage") { auto gr = sgnl->GetGraph(); - const TVector2 tripleMaxAvg = TRestSignalAnalysis::GetTripleMaxAverage(gr); + const TVector2 tripleMaxAvg = TRestPulseShapeAnalysis::GetTripleMaxAverage(gr); const double z = GetHitZCoordinate(tripleMaxAvg.X(), fDriftVelocity, fieldZDirection, zPosition); fHitsEvent->AddHit(x, y, z, tripleMaxAvg.Y(), 0, type); @@ -427,7 +427,7 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven } } else if (fMethod == "intwindow") { auto gr = sgnl->GetGraph(); - auto intWindow = TRestSignalAnalysis::GetIntWindow(gr, fIntWindow); + auto intWindow = TRestPulseShapeAnalysis::GetIntWindow(gr, fIntWindow); for (const auto& [hitTime, energy] : intWindow) { if (energy < fThreshold) continue;