Skip to content

Commit

Permalink
Merge pull request #409 from rest-for-physics/jgalan_fix_restG4
Browse files Browse the repository at this point in the history 
TRestAnalysisTree fixing observable estimators
  • Loading branch information
@jgalan
jgalan authored May 8, 2023
2 parents 3b4459c + 8ed3882 commit d81f90e
Show file tree
Hide file tree
Showing 2 changed files with 115 additions and 67 deletions.
36 changes: 15 additions & 21 deletions source/framework/core/inc/TRestAnalysisTree.h
View file
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
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

0 comments on commit d81f90e

Please sign in to comment.