diff --git a/inc/TRestDetectorSignal.h b/inc/TRestDetectorSignal.h index 2c1fbe67..461c1ab6 100644 --- a/inc/TRestDetectorSignal.h +++ b/inc/TRestDetectorSignal.h @@ -138,7 +138,7 @@ class TRestDetectorSignal { void ExponentialConvolution(Double_t fromTime, Double_t decayTime, Double_t offset = 0); void SignalAddition(TRestDetectorSignal* inSgnl); - Bool_t isSorted(); + Bool_t isSorted() const; void Sort(); void GetDifferentialSignal(TRestDetectorSignal* diffSgnl, Int_t smearPoints = 5); @@ -157,7 +157,7 @@ class TRestDetectorSignal { fSignalCharge.clear(); } - void WriteSignalToTextFile(const TString& filename); + void WriteSignalToTextFile(const TString& filename) const; void Print() const; TGraph* GetGraph(Int_t color = 1); diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index 6e7a2f09..0248cd5a 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -291,31 +291,48 @@ TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the pe 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. + // Define fit limits + Double_t threshold = GetData(maxRaw) * 0.9; // 90% of the maximum value + + Double_t lowerLimit = maxRawTime, upperLimit = maxRawTime; + + // Find the lower limit: time when signal drops to 90% of the max before the peak + for (int i = maxRaw; i >= 0; --i) { + if (GetData(i) <= threshold) { + lowerLimit = GetTime(i); + break; + } + } + + // Find the upper limit: time when signal drops to 90% of the max after the peak + for (int i = maxRaw; i < GetNumberOfPoints(); ++i) { + if (GetData(i) <= threshold) { + upperLimit = GetTime(i); + break; + } + } + + TF1 gaus("gaus", "gaus", lowerLimit, upperLimit); + TH1F h("h", "h", GetNumberOfPoints(), GetTime(0), GetTime(GetNumberOfPoints() - 1)); // copying the signal peak to a histogram for (int i = 0; i < GetNumberOfPoints(); i++) { - h1->Fill(GetTime(i), GetData(i)); + h.SetBinContent(i + 1, GetData(i)); } /* TCanvas* c = new TCanvas("c", "Signal fit", 200, 10, 1280, 720); - h1->GetXaxis()->SetTitle("Time (us)"); - h1->GetYaxis()->SetTitle("Amplitude"); - h1->Draw(); + h->GetXaxis()->SetTitle("Time (us)"); + h->GetYaxis()->SetTitle("Amplitude"); + h->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 + TFitResultPtr fitResult = h.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); + energy = gaus.GetParameter(0); + time = gaus.GetParameter(1); } else { // the fit failed, return -1 to indicate failure energy = -1; @@ -324,24 +341,19 @@ TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the pe << "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" + << "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(); + h->Draw(); c2->Update(); getchar(); delete c2; */ } - TVector2 fitParam(time, energy); - - delete h1; - delete gaus; - - return fitParam; + return {time, energy}; } // z position by landau fit @@ -353,24 +365,42 @@ TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the p 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. + // Define fit limits + Double_t threshold = GetData(maxRaw) * 0.9; // 90% of the maximum value + + Double_t lowerLimit = maxRawTime, upperLimit = maxRawTime; + + // Find the lower limit: time when signal drops to 90% of the max before the peak + for (int i = maxRaw; i >= 0; --i) { + if (GetData(i) <= threshold) { + lowerLimit = GetTime(i); + break; + } + } + + // Find the upper limit: time when signal drops to 90% of the max after the peak + for (int i = maxRaw; i < GetNumberOfPoints(); ++i) { + if (GetData(i) <= threshold) { + upperLimit = GetTime(i); + break; + } + } + + TF1 landau("landau", "landau", lowerLimit, upperLimit); + TH1F h("h", "h", GetNumberOfPoints(), GetTime(0), GetTime(GetNumberOfPoints() - 1)); // copying the signal peak to a histogram for (int i = 0; i < GetNumberOfPoints(); i++) { - h1->Fill(GetTime(i), GetData(i)); + h.SetBinContent(i + 1, 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 + h.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); + energy = landau.GetParameter(0); + time = landau.GetParameter(1); } else { // the fit failed, return -1 to indicate failure energy = -1; @@ -379,24 +409,19 @@ TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the p << "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" + << "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(); + h->Draw(); c2->Update(); getchar(); delete c2; */ } - TVector2 fitParam(time, energy); - - delete h1; - delete landau; - - return fitParam; + return {time, energy}; } // z position by aget fit @@ -420,27 +445,43 @@ TRestDetectorSignal::GetMaxAget() // returns a 2vector with the time of the pea 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); + // Define fit limits + Double_t threshold = GetData(maxRaw) * 0.9; // 90% of the maximum value + + Double_t lowerLimit = maxRawTime, upperLimit = maxRawTime; + + // Find the lower limit: time when signal drops to 90% of the max before the peak + for (int i = maxRaw; i >= 0; --i) { + if (GetData(i) <= threshold) { + lowerLimit = GetTime(i); + break; + } + } + + // Find the upper limit: time when signal drops to 90% of the max after the peak + for (int i = maxRaw; i < GetNumberOfPoints(); ++i) { + if (GetData(i) <= threshold) { + upperLimit = GetTime(i); + break; + } + } + + TF1 aget("aget", agetResponseFunction, lowerLimit, upperLimit, 3); // + TH1F h("h", "h", GetNumberOfPoints(), GetTime(0), GetTime(GetNumberOfPoints() - 1)); + 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)); + h.SetBinContent(i + 1, 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 + TFitResultPtr fitResult = h.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); + energy = aget.GetParameter(0); + time = aget.GetParameter(1); } else { // the fit failed, return -1 to indicate failure energy = -1; @@ -449,18 +490,14 @@ TRestDetectorSignal::GetMaxAget() // returns a 2vector with the time of the pea << "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" + << "Failed fit parameters = " << aget.GetParameter(0) << " || " << aget.GetParameter(1) << " || " + << aget.GetParameter(2) << "\n" << "Assigned fit parameters : energy = " << energy << ", time = " << time << endl; } - TVector2 fitParam(time, energy); - - delete h1; - delete aget; - - return fitParam; + return {time, energy}; } + Double_t TRestDetectorSignal::GetMaxPeakTime(Int_t from, Int_t to) { return GetTime(GetMaxIndex(from, to)); } Double_t TRestDetectorSignal::GetMinPeakValue() { return GetData(GetMinIndex()); } @@ -518,7 +555,7 @@ Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) { return -1; } -Bool_t TRestDetectorSignal::isSorted() { +Bool_t TRestDetectorSignal::isSorted() const { for (int i = 0; i < GetNumberOfPoints() - 1; i++) { if (GetTime(i + 1) < GetTime(i)) { return false; @@ -718,7 +755,7 @@ void TRestDetectorSignal::GetSignalGaussianConvolution(TRestDetectorSignal* conv cout << "Final charge of the pulse " << totChargeFinal << endl; } -void TRestDetectorSignal::WriteSignalToTextFile(const TString& filename) { +void TRestDetectorSignal::WriteSignalToTextFile(const TString& filename) const { FILE* fff = fopen(filename.Data(), "w"); for (int i = 0; i < GetNumberOfPoints(); i++) { fprintf(fff, "%e\t%e\n", GetTime(i), GetData(i));