diff --git a/inc/TRestDetectorSignal.h b/inc/TRestDetectorSignal.h index 461c1ab6..c72ee08b 100644 --- a/inc/TRestDetectorSignal.h +++ b/inc/TRestDetectorSignal.h @@ -28,6 +28,7 @@ #include #include +#include class TRestDetectorSignal { private: @@ -55,9 +56,9 @@ class TRestDetectorSignal { // TODO other objects should probably skip using GetMaxIndex directly Int_t GetMaxIndex(Int_t from = 0, Int_t to = 0); - TVector2 GetMaxGauss(); - TVector2 GetMaxLandau(); - TVector2 GetMaxAget(); + std::optional> GetPeakGauss(); + std::optional> GetPeakLandau(); + std::optional> GetPeakAget(); std::string GetSignalName() const { return fName; } std::string GetSignalType() const { return fType; } @@ -66,11 +67,7 @@ class TRestDetectorSignal { void SetSignalType(const std::string& type) { fType = type; } // Getters - TVector2 GetPoint(Int_t n) { - TVector2 vector2(GetTime(n), GetData(n)); - - return vector2; - } + TVector2 GetPoint(Int_t n) { return {GetTime(n), GetData(n)}; } inline Int_t GetSignalID() const { return fSignalID; } inline Int_t GetID() const { return fSignalID; } diff --git a/src/TRestDetectorSignal.cxx b/src/TRestDetectorSignal.cxx index 0248cd5a..362bdfff 100644 --- a/src/TRestDetectorSignal.cxx +++ b/src/TRestDetectorSignal.cxx @@ -32,6 +32,7 @@ #include #include +#include #include using namespace std; @@ -267,10 +268,16 @@ Int_t TRestDetectorSignal::GetMaxIndex(Int_t from, Int_t to) { Double_t max = std::numeric_limits::min(); Int_t index = 0; - if (from < 0) from = 0; - if (to > GetNumberOfPoints()) to = GetNumberOfPoints(); + if (from < 0) { + from = 0; + } + if (to > GetNumberOfPoints()) { + to = GetNumberOfPoints(); + } - if (to == 0) to = GetNumberOfPoints(); + if (to == 0) { + to = GetNumberOfPoints(); + } for (int i = from; i < to; i++) { if (this->GetData(i) > max) { @@ -284,87 +291,60 @@ Int_t TRestDetectorSignal::GetMaxIndex(Int_t from, Int_t to) { // z position by gaussian fit -TVector2 -TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the peak time in us and the energy +optional> +TRestDetectorSignal::GetPeakGauss() // 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; + const auto indexMax = + std::distance(fSignalCharge.begin(), std::max_element(fSignalCharge.begin(), fSignalCharge.end())); + const auto timeMax = fSignalTime[indexMax]; + const auto signalMax = fSignalCharge[indexMax]; // Define fit limits - Double_t threshold = GetData(maxRaw) * 0.9; // 90% of the maximum value + Double_t threshold = signalMax * 0.9; // 90% of the maximum value - Double_t lowerLimit = maxRawTime, upperLimit = maxRawTime; + Double_t lowerLimit = timeMax, upperLimit = timeMax; // 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); + for (auto i = indexMax; i >= 0; --i) { + if (fSignalCharge[i] <= threshold) { + lowerLimit = fSignalTime[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); + for (auto i = indexMax; i < GetNumberOfPoints(); ++i) { + if (fSignalCharge[i] <= threshold) { + lowerLimit = fSignalTime[i]; break; } } - TF1 gaus("gaus", "gaus", lowerLimit, upperLimit); - TH1F h("h", "h", GetNumberOfPoints(), GetTime(0), GetTime(GetNumberOfPoints() - 1)); + TF1 gauss("gaus", "gaus", lowerLimit, upperLimit); - // copying the signal peak to a histogram - for (int i = 0; i < GetNumberOfPoints(); i++) { - h.SetBinContent(i + 1, GetData(i)); + auto signal_graph = std::unique_ptr(GetGraph()); + + TFitResultPtr fitResult = + signal_graph->Fit(&gauss, "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()) { + return nullopt; } - /* - TCanvas* c = new TCanvas("c", "Signal fit", 200, 10, 1280, 720); - h->GetXaxis()->SetTitle("Time (us)"); - h->GetYaxis()->SetTitle("Amplitude"); - h->Draw(); - */ - 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 + double energy = gauss.GetParameter(0); + double time = gauss.GetParameter(1); - 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); - h->Draw(); - c2->Update(); - getchar(); - delete c2; - */ - } - - return {time, energy}; + return make_pair(time, energy); } // z position by landau fit -TVector2 -TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the peak time in us and the energy +optional> +TRestDetectorSignal::GetPeakLandau() // 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; // Define fit limits Double_t threshold = GetData(maxRaw) * 0.9; // 90% of the maximum value @@ -388,40 +368,20 @@ TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the p } 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++) { - h.SetBinContent(i + 1, GetData(i)); - } + auto signal_graph = std::unique_ptr(GetGraph()); TFitResultPtr fitResult = - 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); - } 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); - h->Draw(); - c2->Update(); - getchar(); - delete c2; - */ - } - - return {time, energy}; + signal_graph->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()) { + return nullopt; + } + + double energy = landau.GetParameter(0); + double time = landau.GetParameter(1); + + return make_pair(time, energy); } // z position by aget fit @@ -438,13 +398,12 @@ Double_t agetResponseFunction(Double_t* x, Double_t* par) { // x contains as ma return f; } -TVector2 -TRestDetectorSignal::GetMaxAget() // returns a 2vector with the time of the peak time in us and the energy +optional> +TRestDetectorSignal::GetPeakAget() // 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; // Define fit limits Double_t threshold = GetData(maxRaw) * 0.9; // 90% of the maximum value @@ -468,34 +427,22 @@ TRestDetectorSignal::GetMaxAget() // returns a 2vector with the time of the pea } 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++) { - h.SetBinContent(i + 1, GetData(i)); + auto signal_graph = std::unique_ptr(GetGraph()); + + TFitResultPtr fitResult = + signal_graph->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()) { + return nullopt; } - 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 + double energy = aget.GetParameter(0); + double time = aget.GetParameter(1); - 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; - } - - return {time, energy}; + return make_pair(time, energy); } Double_t TRestDetectorSignal::GetMaxPeakTime(Int_t from, Int_t to) { return GetTime(GetMaxIndex(from, to)); } diff --git a/src/TRestDetectorSignalToHitsProcess.cxx b/src/TRestDetectorSignalToHitsProcess.cxx index 43f132b4..4782aa07 100644 --- a/src/TRestDetectorSignalToHitsProcess.cxx +++ b/src/TRestDetectorSignalToHitsProcess.cxx @@ -378,90 +378,41 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven cout << "E1, E2, E3 = " << energy1 << ", " << energy2 << ", " << energy3 << endl; } - } else if (fMethod == "gaussFit") { - TVector2 gaussFit = signal->GetMaxGauss(); - - Double_t hitTime = 0; - Double_t z = -1.0; - if (gaussFit.X() >= 0.0) { - hitTime = gaussFit.X(); - Double_t distanceToPlane = hitTime * fDriftVelocity; - z = zPosition + fieldZDirection * distanceToPlane; - } - Double_t energy = -1.0; - if (gaussFit.Y() >= 0.0) { - energy = gaussFit.Y(); - } - - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { - cout << "Signal event : " << signal->GetSignalID() - << "--------------------------------------------------------" << endl; - cout << "GausFit : time bin " << gaussFit.X() << " and energy : " << gaussFit.Y() << endl; - cout << "Signal to hit info : zPosition : " << zPosition - << "; fieldZDirection : " << fieldZDirection << " and driftV : " << fDriftVelocity - << endl; - cout << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << z - << " Energy : " << energy << endl; - } - - fHitsEvent->AddHit(x, y, z, energy, 0, type); - - } else if (fMethod == "landauFit") { - TVector2 landauFit = signal->GetMaxLandau(); - - Double_t hitTime = 0; - Double_t z = -1.0; - if (landauFit.X() >= 0.0) { - hitTime = landauFit.X(); - Double_t distanceToPlane = hitTime * fDriftVelocity; - z = zPosition + fieldZDirection * distanceToPlane; + } else if (fMethod == "gaussFit" || fMethod == "landauFit" || fMethod == "agetFit") { + optional> peak; + if (fMethod == "gaussFit") { + peak = signal->GetPeakGauss(); + } else if (fMethod == "landauFit") { + peak = signal->GetPeakLandau(); + } else if (fMethod == "agetFit") { + peak = signal->GetPeakAget(); + } else { + throw std::runtime_error("Invalid method"); } - Double_t energy = -1.0; - if (landauFit.Y() >= 0.0) { - energy = landauFit.Y(); - } - - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { - cout << "Signal event : " << signal->GetSignalID() - << "--------------------------------------------------------" << endl; - cout << "landauFit : time bin " << landauFit.X() << " and energy : " << landauFit.Y() << endl; - cout << "Signal to hit info : zPosition : " << zPosition - << "; fieldZDirection : " << fieldZDirection << " and driftV : " << fDriftVelocity - << endl; - cout << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << z - << " Energy : " << energy << endl; + if (!peak) { + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) { + cout << "Unable to find peak for signal " << signal->GetSignalID() + << " with method: " << fMethod << endl; + } + continue; } + const auto [time, energy] = peak.value(); - fHitsEvent->AddHit(x, y, z, energy, 0, type); - - } else if (fMethod == "agetFit") { - TVector2 agetFit = signal->GetMaxAget(); - - Double_t hitTime = 0; - Double_t z = -1.0; - if (agetFit.X() >= 0.0) { - hitTime = agetFit.X(); - Double_t distanceToPlane = hitTime * fDriftVelocity; - z = zPosition + fieldZDirection * distanceToPlane; - } - Double_t energy = -1.0; - if (agetFit.Y() >= 0.0) { - energy = agetFit.Y(); - } + const auto distanceToPlane = time * fDriftVelocity; + const auto z = zPosition + fieldZDirection * distanceToPlane; if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { cout << "Signal event : " << signal->GetSignalID() << "--------------------------------------------------------" << endl; - cout << "agetFit : time bin " << agetFit.X() << " and energy : " << agetFit.Y() << endl; + cout << "Method: " << fMethod << " : time " << time << " us and energy : " << energy << endl; cout << "Signal to hit info : zPosition : " << zPosition << "; fieldZDirection : " << fieldZDirection << " and driftV : " << fDriftVelocity << endl; - cout << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << z + cout << "Adding hit. Time : " << time << " us x : " << x << " y : " << y << " z : " << z << " Energy : " << energy << endl; } fHitsEvent->AddHit(x, y, z, energy, 0, type); - } else if (fMethod == "qCenter") { Double_t energy_signal = 0; Double_t distanceToPlane = 0; @@ -496,7 +447,7 @@ TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEven } } else if (fMethod == "intwindow") { Int_t nPoints = signal->GetNumberOfPoints(); - std::map > windowMap; + std::map> windowMap; for (int j = 0; j < nPoints; j++) { int index = signal->GetTime(j) / fIntWindow; auto it = windowMap.find(index);