Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

TRestAnalysisTree fixing observable estimators #409

Merged
merged 3 commits into from
May 8, 2023
Merged

TRestAnalysisTree fixing observable estimators #409

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c1f4d88
TRestAnalysisTree::GetObservableAverage. Changing the way average is …
jgalan May 4, 2023
a5af002
TRestAnalysisTree::GetObservableAverage. Now it implements a range
jgalan May 4, 2023
8ed3882
TRestAnalysisTree. Updating GetObservableIntegral,Maximum,Minimum,RMS
jgalan May 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 15 additions & 21 deletions source/framework/core/inc/TRestAnalysisTree.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -279,45 +279,39 @@ 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, 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) {
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,
Int_t nBins = -1, Double_t xLow = -1, Double_t xHigh = -1) {
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,
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);
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);
Expand Down
146 changes: 100 additions & 46 deletions source/framework/core/src/TRestAnalysisTree.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -990,76 +990,130 @@ 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;
}

///////////////////////////////////////////////
/// \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) {
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();
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++) {
TTree::GetEntry(n);
Double_t value = GetDblObservableValue(id);

if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue;
N++;
sum += value;
}

if (N <= 0) return 0;

return sum / N;
}

///////////////////////////////////////////////
/// \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;
}

///////////////////////////////////////////////
Expand Down Expand Up @@ -1094,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)
Expand Down