From c1f4d8889b425fb842efd50756dc61b38f986b25 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Thu, 4 May 2023 10:40:46 +0200 Subject: [PATCH 1/3] TRestAnalysisTree::GetObservableAverage. Changing the way average is calculated --- source/framework/core/src/TRestAnalysisTree.cxx | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/source/framework/core/src/TRestAnalysisTree.cxx b/source/framework/core/src/TRestAnalysisTree.cxx index eea161e65..139f22192 100644 --- a/source/framework/core/src/TRestAnalysisTree.cxx +++ b/source/framework/core/src/TRestAnalysisTree.cxx @@ -1008,13 +1008,16 @@ Double_t TRestAnalysisTree::GetObservableIntegral(const TString& obsName, Double /// Double_t TRestAnalysisTree::GetObservableAverage(const TString& obsName, Double_t xLow, Double_t xHigh, Int_t nBins) { - TString histDefinition = Form("havg(%5d,%lf,%lf)", nBins, xLow, xHigh); - if (xHigh == -1) - this->Draw(obsName + ">>havg"); - else - this->Draw(obsName + ">>" + histDefinition); - TH1F* htemp = (TH1F*)gPad->GetPrimitive("havg"); - return htemp->GetMean(); + Int_t id = GetObservableID((std::string)obsName); + + Double_t sum = 0; + for (Int_t n = 0; n < TTree::GetEntries(); n++) { + TTree::GetEntry(n); + sum += GetDblObservableValue(id); + } + if (TTree::GetEntries() <= 0) return 0; + + return sum / TTree::GetEntries(); } /////////////////////////////////////////////// From a5af002d6bc60b2cf3b28f57e6d46a0046bdf4ca Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Thu, 4 May 2023 22:34:04 +0200 Subject: [PATCH 2/3] TRestAnalysisTree::GetObservableAverage. Now it implements a range --- source/framework/core/inc/TRestAnalysisTree.h | 7 +++---- source/framework/core/src/TRestAnalysisTree.cxx | 13 +++++++++---- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/source/framework/core/inc/TRestAnalysisTree.h b/source/framework/core/inc/TRestAnalysisTree.h index 0be736b3a..a243c6973 100644 --- a/source/framework/core/inc/TRestAnalysisTree.h +++ b/source/framework/core/inc/TRestAnalysisTree.h @@ -284,8 +284,8 @@ class TRestAnalysisTree : public TTree { return GetObservableIntegral(obsName, xLow, xHigh, nBins); } - Double_t GetAverage(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, Int_t nBins = 1000) { - return GetObservableAverage(obsName, xLow, xHigh, nBins); + Double_t GetAverage(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1) { + return GetObservableAverage(obsName, xLow, xHigh); } Double_t GetRMS(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, Int_t nBins = 1000) { @@ -308,8 +308,7 @@ class TRestAnalysisTree : public TTree { Double_t GetObservableIntegral(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, Int_t nBins = 1000); - Double_t GetObservableAverage(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, - Int_t nBins = 1000); + Double_t GetObservableAverage(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1); Double_t GetObservableRMS(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, Int_t nBins = 1000); diff --git a/source/framework/core/src/TRestAnalysisTree.cxx b/source/framework/core/src/TRestAnalysisTree.cxx index 139f22192..f20a53a8c 100644 --- a/source/framework/core/src/TRestAnalysisTree.cxx +++ b/source/framework/core/src/TRestAnalysisTree.cxx @@ -1006,18 +1006,23 @@ Double_t TRestAnalysisTree::GetObservableIntegral(const TString& obsName, Double /// \brief It returns the average of the observable considering the given range. If no range is given /// the full histogram range will be considered. /// -Double_t TRestAnalysisTree::GetObservableAverage(const TString& obsName, Double_t xLow, Double_t xHigh, - Int_t nBins) { +Double_t TRestAnalysisTree::GetObservableAverage(const TString& obsName, Double_t xLow, Double_t xHigh) { Int_t id = GetObservableID((std::string)obsName); Double_t sum = 0; + Int_t N = 0; for (Int_t n = 0; n < TTree::GetEntries(); n++) { TTree::GetEntry(n); + Double_t value = GetDblObservableValue(id); + + if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue; + N++; sum += GetDblObservableValue(id); } - if (TTree::GetEntries() <= 0) return 0; - return sum / TTree::GetEntries(); + if (N <= 0) return 0; + + return sum / N; } /////////////////////////////////////////////// From 8ed3882dc3d67cf495dd06c0162f10284f392b04 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Fri, 5 May 2023 08:49:44 +0200 Subject: [PATCH 3/3] TRestAnalysisTree. Updating GetObservableIntegral,Maximum,Minimum,RMS --- source/framework/core/inc/TRestAnalysisTree.h | 29 ++--- .../framework/core/src/TRestAnalysisTree.cxx | 122 ++++++++++++------ 2 files changed, 96 insertions(+), 55 deletions(-) diff --git a/source/framework/core/inc/TRestAnalysisTree.h b/source/framework/core/inc/TRestAnalysisTree.h index a243c6973..69975da7e 100644 --- a/source/framework/core/inc/TRestAnalysisTree.h +++ b/source/framework/core/inc/TRestAnalysisTree.h @@ -279,25 +279,24 @@ class TRestAnalysisTree : public TTree { void EnableQuickObservableValueSetting(); void DisableQuickObservableValueSetting(); - Double_t GetIntegral(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, - Int_t nBins = 1000) { - return GetObservableIntegral(obsName, xLow, xHigh, nBins); + Double_t GetIntegral(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1) { + return GetObservableIntegral(obsName, xLow, xHigh); } Double_t GetAverage(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1) { return GetObservableAverage(obsName, xLow, xHigh); } - Double_t GetRMS(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, Int_t nBins = 1000) { - return GetObservableRMS(obsName, xLow, xHigh, nBins); + Double_t GetRMS(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1) { + return GetObservableRMS(obsName, xLow, xHigh); } - Double_t GetMinimum(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, Int_t nBins = 1000) { - return GetObservableMinimum(obsName, xLow, xHigh, nBins); + Double_t GetMinimum(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1) { + return GetObservableMinimum(obsName, xLow, xHigh); } - Double_t GetMaximum(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, Int_t nBins = 1000) { - return GetObservableMaximum(obsName, xLow, xHigh, nBins); + Double_t GetMaximum(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1) { + return GetObservableMaximum(obsName, xLow, xHigh); } Double_t GetContour(const TString& obsName, const TString& obsIndexer, Double_t level = 0.5, @@ -305,18 +304,14 @@ class TRestAnalysisTree : public TTree { return GetObservableContour(obsName, obsIndexer, level, nBins, xLow, xHigh); } - Double_t GetObservableIntegral(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, - Int_t nBins = 1000); + Double_t GetObservableIntegral(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1); Double_t GetObservableAverage(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1); - Double_t GetObservableRMS(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, - Int_t nBins = 1000); + Double_t GetObservableRMS(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1); - Double_t GetObservableMinimum(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, - Int_t nBins = 1000); - Double_t GetObservableMaximum(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1, - Int_t nBins = 1000); + Double_t GetObservableMinimum(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1); + Double_t GetObservableMaximum(const TString& obsName, Double_t xLow = -1, Double_t xHigh = -1); Double_t GetObservableContour(const TString& obsName, const TString& obsIndexer, Double_t level = 0.5, Int_t nBins = -1, Double_t xLow = -1, Double_t xHigh = -1); diff --git a/source/framework/core/src/TRestAnalysisTree.cxx b/source/framework/core/src/TRestAnalysisTree.cxx index f20a53a8c..01d2c6622 100644 --- a/source/framework/core/src/TRestAnalysisTree.cxx +++ b/source/framework/core/src/TRestAnalysisTree.cxx @@ -990,16 +990,25 @@ void TRestAnalysisTree::DisableQuickObservableValueSetting() { this->fQuickSetOb /// \brief It returns the integral of the observable considering the given range. If no range is given /// the full histogram range will be considered. /// -Double_t TRestAnalysisTree::GetObservableIntegral(const TString& obsName, Double_t xLow, Double_t xHigh, - Int_t nBins) { - TString histDefinition = Form("hint(%5d,%lf,%lf)", nBins, xLow, xHigh); - if (xHigh == -1) - this->Draw(obsName + ">>hint", obsName); - else - this->Draw(obsName + ">>" + histDefinition, obsName); +Double_t TRestAnalysisTree::GetObservableIntegral(const TString& obsName, Double_t xLow, Double_t xHigh) { + Int_t id = GetObservableID((std::string)obsName); + + if (id < 0) { + RESTError << "TRestAnalysisTree::GetObservableIntegral. Observable not found! : " << obsName + << RESTendl; + return 0; + } + + Double_t sum = 0; + for (Int_t n = 0; n < TTree::GetEntries(); n++) { + TTree::GetEntry(n); + Double_t value = GetDblObservableValue(id); + + if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue; + sum += value; + } - TH1F* htemp = (TH1F*)gPad->GetPrimitive("hint"); - return htemp->Integral(); + return sum; } /////////////////////////////////////////////// @@ -1009,6 +1018,12 @@ Double_t TRestAnalysisTree::GetObservableIntegral(const TString& obsName, Double Double_t TRestAnalysisTree::GetObservableAverage(const TString& obsName, Double_t xLow, Double_t xHigh) { Int_t id = GetObservableID((std::string)obsName); + if (id < 0) { + RESTError << "TRestAnalysisTree::GetObservableAverage. Observable not found! : " << obsName + << RESTendl; + return 0; + } + Double_t sum = 0; Int_t N = 0; for (Int_t n = 0; n < TTree::GetEntries(); n++) { @@ -1017,7 +1032,7 @@ Double_t TRestAnalysisTree::GetObservableAverage(const TString& obsName, Double_ if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue; N++; - sum += GetDblObservableValue(id); + sum += value; } if (N <= 0) return 0; @@ -1029,45 +1044,76 @@ Double_t TRestAnalysisTree::GetObservableAverage(const TString& obsName, Double_ /// \brief It returns the RMS of the observable considering the given range. If no range is given /// the full histogram range will be considered. /// -Double_t TRestAnalysisTree::GetObservableRMS(const TString& obsName, Double_t xLow, Double_t xHigh, - Int_t nBins) { - TString histDefinition = Form("hrms(%5d,%lf,%lf)", nBins, xLow, xHigh); - if (xHigh == -1) - this->Draw(obsName + ">>hrms"); - else - this->Draw(obsName + ">>" + histDefinition); - TH1F* htemp = (TH1F*)gPad->GetPrimitive("hrms"); - return htemp->GetRMS(); +Double_t TRestAnalysisTree::GetObservableRMS(const TString& obsName, Double_t xLow, Double_t xHigh) { + Double_t mean = GetObservableAverage(obsName, xLow, xHigh); + + Int_t id = GetObservableID((std::string)obsName); + + Double_t sum = 0; + Int_t N = 0; + for (Int_t n = 0; n < TTree::GetEntries(); n++) { + TTree::GetEntry(n); + Double_t value = GetDblObservableValue(id); + + if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue; + N++; + + sum += (value - mean) * (value - mean); + } + + if (N <= 0) return 0; + + return TMath::Sqrt(sum / N); } /////////////////////////////////////////////// /// \brief It returns the maximum value of obsName considering the given range. If no range is given /// the full histogram range will be considered. /// -Double_t TRestAnalysisTree::GetObservableMaximum(const TString& obsName, Double_t xLow, Double_t xHigh, - Int_t nBins) { - TString histDefinition = Form("hmax(%5d,%lf,%lf)", nBins, xLow, xHigh); - if (xHigh == -1) - this->Draw(obsName + ">>hmax"); - else - this->Draw(obsName + ">>" + histDefinition); - TH1F* htemp = (TH1F*)gPad->GetPrimitive("hmax"); - return htemp->GetMaximumStored(); +Double_t TRestAnalysisTree::GetObservableMaximum(const TString& obsName, Double_t xLow, Double_t xHigh) { + Int_t id = GetObservableID((std::string)obsName); + + if (id < 0) { + RESTError << "TRestAnalysisTree::GetObservableMaximum. Observable not found! : " << obsName + << RESTendl; + return 0; + } + + Double_t max = -DBL_MAX; + for (Int_t n = 0; n < TTree::GetEntries(); n++) { + TTree::GetEntry(n); + Double_t value = GetDblObservableValue(id); + + if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue; + if (value > max) max = value; + } + + return max; } /////////////////////////////////////////////// /// \brief It returns the minimum value of obsName considering the given range. If no range is given /// the full histogram range will be considered. /// -Double_t TRestAnalysisTree::GetObservableMinimum(const TString& obsName, Double_t xLow, Double_t xHigh, - Int_t nBins) { - TString histDefinition = Form("hmin(%5d,%lf,%lf)", nBins, xLow, xHigh); - if (xHigh == -1) - this->Draw(obsName + ">>hmin"); - else - this->Draw(obsName + ">>" + histDefinition); - TH1F* htemp = (TH1F*)gPad->GetPrimitive("hmin"); - return htemp->GetMinimumStored(); +Double_t TRestAnalysisTree::GetObservableMinimum(const TString& obsName, Double_t xLow, Double_t xHigh) { + Int_t id = GetObservableID((std::string)obsName); + + if (id < 0) { + RESTError << "TRestAnalysisTree::GetObservableMinimum. Observable not found! : " << obsName + << RESTendl; + return 0; + } + + Double_t min = DBL_MAX; + for (Int_t n = 0; n < TTree::GetEntries(); n++) { + TTree::GetEntry(n); + Double_t value = GetDblObservableValue(id); + + if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue; + if (value < min) min = value; + } + + return min; } /////////////////////////////////////////////// @@ -1102,7 +1148,7 @@ Double_t TRestAnalysisTree::GetObservableContour(const TString& obsName, const T return 0; } - Double_t integral = this->GetIntegral(obsWeight, xLow, xHigh, nBins); + Double_t integral = this->GetIntegral(obsWeight, xLow, xHigh); TString histDefinition = Form("hc(%5d,%lf,%lf)", nBins, xLow, xHigh); if (xHigh == -1)