diff --git a/inc/TRestDetectorSignal.h b/inc/TRestDetectorSignal.h index 11eee7bb..99a69867 100644 --- a/inc/TRestDetectorSignal.h +++ b/inc/TRestDetectorSignal.h @@ -42,13 +42,7 @@ class TRestDetectorSignal { public: TGraph* fGraph; //! - std::vector fPointsOverThreshold; //! - - void IncreaseAmplitude(TVector2 p); - void SetPoint(TVector2 p); - - // 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(); @@ -87,16 +81,13 @@ class TRestDetectorSignal { 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); 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(); } @@ -112,16 +103,13 @@ 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); - 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); - - Double_t SubstractBaseline(Int_t startBin, Int_t endBin); void AddOffset(Double_t offset); void MultiplySignalBy(Double_t factor); diff --git a/inc/TRestDetectorSignalToHitsProcess.h b/inc/TRestDetectorSignalToHitsProcess.h index b59d6cdf..96ca9e39 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 * driftVelocity; + 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 b26b0f14..129036cc 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -26,6 +26,7 @@ #include "TRestDetectorSignal.h" #include "TFitResult.h" +#include "TRestPulseShapeAnalysis.h" using namespace std; #include @@ -43,8 +44,6 @@ TRestDetectorSignal::TRestDetectorSignal() { fSignalID = -1; fSignalTime.clear(); fSignalCharge.clear(); - - fPointsOverThreshold.clear(); } TRestDetectorSignal::~TRestDetectorSignal() { @@ -71,17 +70,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()); } } @@ -92,17 +89,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()); } } @@ -131,10 +126,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 TRestPulseShapeAnalysis::GetIntegral(fSignalCharge, startBin, endBin); } void TRestDetectorSignal::Normalize(Double_t scale) { @@ -151,192 +143,27 @@ Double_t TRestDetectorSignal::GetIntegralWithTime(Double_t startTime, Double_t e 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; -} - -/* {{{ -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(); - 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 TRestPulseShapeAnalysis::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 TRestPulseShapeAnalysis::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 TRestPulseShapeAnalysis::GetMaxGauss(gr); } // z position by landau fit @@ -344,157 +171,33 @@ 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 TRestPulseShapeAnalysis::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 TRestPulseShapeAnalysis::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 TRestPulseShapeAnalysis::GetMaxBin(fSignalCharge); } - return index; -} +Int_t TRestDetectorSignal::GetMinIndex() { return TRestPulseShapeAnalysis::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 = TRestPulseShapeAnalysis::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 = TRestPulseShapeAnalysis::GetMaxBin(fSignalTime); + return GetTime(index); } Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) { @@ -505,12 +208,7 @@ Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) { 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()) { @@ -554,53 +252,9 @@ 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->IncreaseAmplitude(GetTime(i), sum / 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->IncreaseAmplitude(this->GetTime(i), sum / averagingPoints); - } - - for (int i = GetNumberOfPoints() - averagingPoints / 2; i < GetNumberOfPoints(); i++) - smthSignal->IncreaseAmplitude(GetTime(i), sum / averagingPoints); -} - -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); -} + auto smoothed = TRestPulseShapeAnalysis::GetSignalSmoothed(fSignalCharge, averagingPoints); -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); - - AddOffset(-bL); - - return bL; + for (int i = 0; i < GetNumberOfPoints(); i++) smthSignal->IncreaseAmplitude(GetTime(i), smoothed[i]); } void TRestDetectorSignal::AddOffset(Double_t offset) { @@ -657,11 +311,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)); } } @@ -717,7 +368,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; @@ -726,11 +377,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 98149350..ad4365ac 100644 --- a/src/TRestDetectorSignalToHitsProcess.cxx +++ b/src/TRestDetectorSignalToHitsProcess.cxx @@ -126,6 +126,7 @@ #include "TRestDetectorSignalToHitsProcess.h" #include "TRestDetectorSetup.h" +#include "TRestPulseShapeAnalysis.h" using namespace std; @@ -276,13 +277,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) @@ -291,76 +286,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); - - hitTime = sgnl->GetTime(binprev); - energy = sgnl->GetData(binprev); + auto gr = sgnl->GetGraph(); + auto tripleMax = TRestPulseShapeAnalysis::GetTripleMax(gr); - 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 = TRestPulseShapeAnalysis::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") { @@ -370,8 +318,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) { @@ -398,8 +345,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) { @@ -426,8 +372,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) { @@ -448,30 +393,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 @@ -480,30 +424,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 = TRestPulseShapeAnalysis::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;