 |
ROOT
Reference Guide |
|
| |
ROOT
Reference Guide |
|
#include "TEfficiency.h"#include "TROOT.h"#include "TBuffer.h"#include "TGraphAsymmErrors.h"#include "TGraphErrors.h"#include "TStyle.h"#include "TMath.h"#include "TVirtualPad.h"#include "TF1.h"#include "TH1.h"#include "TVector.h"#include "TVectorD.h"#include "TSystem.h"#include "Math/QuantFuncMathCore.h"#include "strtok.h"#include <cstring>#include <iostream>#include <fstream>Namespaces | |
| namespace | ROOT |
| tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tbb::task_arena without forward declaring tbb::interface7 | |
Functions | |
| for (i=0;i< fNpoints;i++) | |
| else out<< " TGraphAsymmErrors *";out<< "grae = new TGraphAsymmErrors("<< fNpoints<< ","<< fXName<< ","<< fYName<< ","<< fElXName<< ","<< fEhXName<< ","<< fElYName<< ","<< fEhYName<< ");"<< std::endl;SaveHistogramAndFunctions(out, "grae", frameNumber, option);}void TGraphAsymmErrors::Scale(Double_t c1, Option_t *option){ TGraph::Scale(c1, option);TString opt=option;opt.ToLower();if(opt.Contains("x") &&GetEXlow()) { for(Int_t i=0;i< GetN();i++) GetEXlow()[i] *=c1;} if(opt.Contains("x") &&GetEXhigh()) { for(Int_t i=0;i< GetN();i++) GetEXhigh()[i] *=c1;} if(opt.Contains("y") &&GetEYlow()) { for(Int_t i=0;i< GetN();i++) GetEYlow()[i] *=c1;} if(opt.Contains("y") &&GetEYhigh()) { for(Int_t i=0;i< GetN();i++) GetEYhigh()[i] *=c1;}}void TGraphAsymmErrors::SetPointError(Double_t exl, Double_t exh, Double_t eyl, Double_t eyh){ if(!(TVirtualPad::Pad())) { Error("SetPointError", "Cannot be used without gPad, requires last mouse position");return;} Int_t px=(TVirtualPad::Pad()) -> | GetEventX () |
| if ((ROOT::GetROOT()) ->ClassSaved(TGraphAsymmErrors::Class())) out<<" " | |
| if (ipoint==-2) return | |
Variables | |
| auto | fEhXName = SaveArray(out, "fehx", frameNumber, fEXhigh) |
| auto | fEhYName = SaveArray(out, "fehy", frameNumber, fEYhigh) |
| auto | fElXName = SaveArray(out, "felx", frameNumber, fEXlow) |
| auto | fElYName = SaveArray(out, "fely", frameNumber, fEYlow) |
| fEXhigh [ipoint] = exh | |
| fEXlow [ipoint] = exl | |
| fEYhigh [ipoint] = eyh | |
| fEYlow [ipoint] = eyl | |
| static Int_t | frameNumber = 3000 |
| auto | fXName = SaveArray(out, "fx", frameNumber, fX) |
| auto | fYName = SaveArray(out, "fy", frameNumber, fY) |
| Int_t | i |
| Int_t | ipoint = -2 |
| Int_t | py = (TVirtualPad::Pad()) ->GetEventY() |
| for | ( | ) |
Definition at line 1308 of file TGraphAsymmErrors.cxx.
| else out<< " TGraphAsymmErrors *";out<< "grae = new TGraphAsymmErrors("<< fNpoints<< ","<< fXName<< ","<< fYName<< ","<< fElXName<< ","<< fEhXName<< ","<< fElYName<< ","<< fEhYName<< ");"<< std::endl;SaveHistogramAndFunctions(out, "grae", frameNumber, option);}void TGraphAsymmErrors::Scale(Double_t c1, Option_t *option){ TGraph::Scale(c1, option);TString opt=option;opt.ToLower();if(opt.Contains("x") &&GetEXlow()) { for(Int_t i=0;i< GetN();i++) GetEXlow()[i] *=c1;} if(opt.Contains("x") &&GetEXhigh()) { for(Int_t i=0;i< GetN();i++) GetEXhigh()[i] *=c1;} if(opt.Contains("y") &&GetEYlow()) { for(Int_t i=0;i< GetN();i++) GetEYlow()[i] *=c1;} if(opt.Contains("y") &&GetEYhigh()) { for(Int_t i=0;i< GetN();i++) GetEYhigh()[i] *=c1;}}void TGraphAsymmErrors::SetPointError(Double_t exl, Double_t exh, Double_t eyl, Double_t eyh){ if(!(TVirtualPad::Pad())) { Error("SetPointError", "Cannot be used without gPad, requires last mouse position");return;} Int_t px=(TVirtualPad::Pad()) -> GetEventX | ( | ) |
| if | ( | (ROOT::GetROOT()) ->ClassSaved(TGraphAsymmErrors::Class()) | ) |
| if | ( | ipoint | = =-2 | ) |
| auto fEhXName = SaveArray(out, "fehx", frameNumber, fEXhigh) |
Definition at line 1246 of file TGraphAsymmErrors.cxx.
| auto fEhYName = SaveArray(out, "fehy", frameNumber, fEYhigh) |
Definition at line 1247 of file TGraphAsymmErrors.cxx.
| auto fElXName = SaveArray(out, "felx", frameNumber, fEXlow) |
Definition at line 1244 of file TGraphAsymmErrors.cxx.
| auto fElYName = SaveArray(out, "fely", frameNumber, fEYlow) |
Definition at line 1245 of file TGraphAsymmErrors.cxx.
| fEXhigh[ipoint] = exh |
Definition at line 1317 of file TGraphAsymmErrors.cxx.
| fEXlow[ipoint] = exl |
Definition at line 1315 of file TGraphAsymmErrors.cxx.
| fEYhigh[ipoint] = eyh |
Definition at line 1318 of file TGraphAsymmErrors.cxx.
| fEYlow[ipoint] = eyl |
Definition at line 1316 of file TGraphAsymmErrors.cxx.
|
static |
TGraphAsymmErrors constructor reading input from filename filename is assumed to contain at least 2 columns of numbers convention for format (default=`"%lg %lg %lg %lg %lg %lg"`) - format = `"%lg %lg"` read only 2 first columns into X, Y - format = `"%lg %lg %lg %lg"` read only 4 first columns into X, Y, ELY, EHY - format = `"%lg %lg %lg %lg %lg %lg"` read only 6 first columns into X, Y, EXL, EYH, EYL, EHY For files separated by a specific delimiter different from `' '` and `'\\t'` (e.g. `';'` in csv files) you can avoid using `%*s` to bypass this delimiter by explicitly specify the `"option" argument,
/ e.g. option=" \\t,;" for columns of figures separated by any of these characters ‘(’ ', '\t', ',', ';') / used once(e.g. "1;1")or in a combined way(" 1;,;; 1"). / Note in that case, the instantiation is about 2 times slower. / In case a delimiter is specified, the format"%lg %lg %lg"` will read X,Y,EX.
TGraphAsymmErrors::TGraphAsymmErrors(const char *filename, const char *format, Option_t *option) : TGraph(100) { if (!CtorAllocate()) return; Double_t x, y, exl, exh, eyl, eyh; TString fname = filename; gSystem->ExpandPathName(fname); std::ifstream infile(fname.Data()); if (!infile.good()) { MakeZombie(); Error("TGraphAsymmErrors", "Cannot open file: %s, TGraphAsymmErrors is Zombie", filename); fNpoints = 0; return; } std::string line; Int_t np = 0;
if (strcmp(option, "") == 0) { // No delimiters specified (standard constructor).
Int_t ncol = TGraphErrors::CalculateScanfFields(format); //count number of columns in format Int_t res; while (std::getline(infile, line, '
')) { exl = exh = eyl = eyh = 0.; if (ncol < 3) { res = sscanf(line.c_str(), format, &x, &y); } else if (ncol < 5) { res = sscanf(line.c_str(), format, &x, &y, &eyl, &eyh); } else { res = sscanf(line.c_str(), format, &x, &y, &exl, &exh, &eyl, &eyh); } if (res < 2) { continue; //skip empty and ill-formed lines } SetPoint(np, x, y); SetPointError(np, exl, exh, eyl, eyh); np++; } Set(np);
} else { // A delimiter has been specified in "option"
Checking format and creating its boolean equivalent TString format_ = TString(format) ; format_.ReplaceAll(" ", "") ; format_.ReplaceAll("\t", "") ; format_.ReplaceAll("lg", "") ; format_.ReplaceAll("s", "") ; format_.ReplaceAll("%*", "0") ; format_.ReplaceAll("%", "1") ; if (!format_.IsDigit()) { Error("TGraphAsymmErrors", "Incorrect input format! Allowed format tags are {\"%%lg\",\"%%*lg\" or \"%%*s\"}"); return ; } Int_t ntokens = format_.Length() ; if (ntokens < 2) { Error("TGraphAsymmErrors", "Incorrect input format! Only %d tag(s) in format whereas at least 2 \"%%lg\" tags are expected!", ntokens); return ; } Int_t ntokensToBeSaved = 0; Bool_t * isTokenToBeSaved = new Bool_t[ntokens]; for (Int_t idx = 0; idx < ntokens; idx++) { isTokenToBeSaved[idx] = TString::Format("%c", format_[idx]).Atoi(); //atoi(&format_[idx]) does not work for some reason... if (isTokenToBeSaved[idx] == 1) { ntokensToBeSaved++ ; } } if (ntokens >= 2 && (ntokensToBeSaved < 2 || ntokensToBeSaved > 4)) { //first condition not to repeat the previous error message Error("TGraphAsymmErrors", "Incorrect input format! There are %d \"%%lg\" tag(s) in format whereas 2,3 or 4 are expected!", ntokensToBeSaved); delete [] isTokenToBeSaved; return ; }
Initializing loop variables Bool_t isLineToBeSkipped = kFALSE; //empty and ill-formed lines char *token = nullptr; TString token_str = ""; Int_t token_idx = 0; Double_t value[6]; //x,y,exl, exh, eyl, eyh buffers for (Int_t k = 0; k < 6; k++) value[k] = 0.; Int_t value_idx = 0;
Looping char *rest; while (std::getline(infile, line, '
')) { if (!line.empty()) { if (line[line.size() - 1] == char(13)) { // removing DOS CR character line.erase(line.end() - 1, line.end()) ; } token = R__STRTOK_R(const_cast<char*>(line.c_str()), option, &rest) ; while (token != nullptr && value_idx < ntokensToBeSaved) { if (isTokenToBeSaved[token_idx]) { token_str = TString(token) ; token_str.ReplaceAll("\t", "") ; if (!token_str.IsFloat()) { isLineToBeSkipped = kTRUE ; break ; } else { value[value_idx] = token_str.Atof() ; value_idx++ ; } } token = R__STRTOK_R(nullptr, option, &rest); // next token token_idx++ ; } if (!isLineToBeSkipped && value_idx > 1) { //i.e. 2,3 or 4 x = value[0]; y = value[1]; exl = value[2]; exh = value[3]; eyl = value[4]; eyh = value[5]; SetPoint(np, x, y); SetPointError(np, exl, exh, eyl, eyh); np++ ; } } isLineToBeSkipped = kFALSE; token = nullptr; token_idx = 0; value_idx = 0; } Set(np) ;
Cleaning delete [] isTokenToBeSaved; delete token; } infile.close(); }
////////////////////////////////////////////////////////////////////////////// / TGraphAsymmErrors default destructor.
TGraphAsymmErrors::~TGraphAsymmErrors() { if(fEXlow) delete [] fEXlow; if(fEXhigh) delete [] fEXhigh; if(fEYlow) delete [] fEYlow; if(fEYhigh) delete [] fEYhigh; }
////////////////////////////////////////////////////////////////////////////// / Allocate internal data structures for size points.
Double_t** TGraphAsymmErrors::Allocate(Int_t size) { return AllocateArrays(6, size); }
////////////////////////////////////////////////////////////////////////////// / Apply a function to all data points \( y = f(x,y) \) / / Errors are calculated as \( eyh = f(x,y+eyh)-f(x,y) \) and / \( eyl = f(x,y)-f(x,y-eyl) \) / / Special treatment has to be applied for the functions where the / role of "up" and "down" is reversed. / / Function suggested/implemented by Miroslav Helbich helbi.nosp@m.ch@m.nosp@m.ail.d.nosp@m.esy..nosp@m.de
void TGraphAsymmErrors::Apply(TF1 *f) { Double_t x,y,exl,exh,eyl,eyh,eyl_new,eyh_new,fxy;
if (fHistogram) { delete fHistogram; fHistogram = nullptr; } for (Int_t i=0;i<GetN();i++) { GetPoint(i,x,y); exl = GetErrorXlow(i); exh = GetErrorXhigh(i); eyl = GetErrorYlow(i); eyh = GetErrorYhigh(i);
fxy = f->Eval(x,y); SetPoint(i,x,fxy);
in the case of the functions like y-> -1*y the roles of the upper and lower error bars is reversed if (f->Eval(x,y-eyl)<f->Eval(x,y+eyh)) { eyl_new = TMath::Abs(fxy - f->Eval(x,y-eyl)); eyh_new = TMath::Abs(f->Eval(x,y+eyh) - fxy); } else { eyh_new = TMath::Abs(fxy - f->Eval(x,y-eyl)); eyl_new = TMath::Abs(f->Eval(x,y+eyh) - fxy); }
error on x doesn't change SetPointError(i,exl,exh,eyl_new,eyh_new); } if ( (TVirtualPad::Pad()) ) (TVirtualPad::Pad()) ->Modified(); }
////////////////////////////////////////////////////////////////////////////// /This function is only kept for backward compatibility. /You should rather use the Divide method. /It calls Divide(pass,total,"cl=0.683 b(1,1) mode") which is equivalent to the /former BayesDivide method.
void TGraphAsymmErrors::BayesDivide(const TH1* pass, const TH1* total, Option_t *) { Divide(pass,total,"cl=0.683 b(1,1) mode"); }
////////////////////////////////////////////////////////////////////////////// / Fill this TGraphAsymmErrors by dividing two 1-dimensional histograms pass/total / / This method serves two purposes: / / ### 1) calculating efficiencies: / / The assumption is that the entries in "pass" are a subset of those in / "total". That is, we create an "efficiency" graph, where each entry is / between 0 and 1, inclusive. / / If the histograms are not filled with unit weights, the number of effective / entries is used to normalise the bin contents which might lead to wrong results. /
\[/// \text{effective entries} = \frac{(\sum w_{i})^{2}}{\sum w_{i}^{2}} /// \]
/ The points are assigned a x value at the center of each histogram bin. / The y values are \(\text{eff} = \frac{\text{pass}}{\text{total}}\) / for all options except for the / bayesian methods where the result depends on the chosen option. / / If the denominator becomes 0 or pass > total, the corresponding bin is / skipped. / / ### 2) calculating ratios of two Poisson means (option 'pois'): / / The two histograms are interpreted as independent Poisson processes and the ratio /
\[/// \tau = \frac{n_{1}}{n_{2}} = \frac{\varepsilon}{1 - \varepsilon} /// \]
/ with \(\varepsilon = \frac{n_{1}}{n_{1} + n_{2}}\). / The histogram 'pass' is interpreted as \(n_{1}\) and the total histogram / is used for \(n_{2}\). / / The (asymmetric) uncertainties of the Poisson ratio are linked to the uncertainties / of efficiency by a parameter transformation: /
\[/// \Delta \tau_{low/up} = \frac{1}{(1 - \varepsilon)^{2}} \Delta \varepsilon_{low/up} /// \]
/ / The x errors span each histogram bin (lowedge ... lowedge+width) / The y errors depend on the chosen statistic methode which can be determined / by the options given below. For a detailed description of the used statistic / calculations please have a look at the corresponding functions! / / Options: / - v : verbose mode: prints information about the number of used bins / and calculated efficiencies with their errors / - cl=x : determine the used confidence level (0<x<1) (default is 0.683) / - cp : Clopper-Pearson interval (see TEfficiency::ClopperPearson) / - w : Wilson interval (see TEfficiency::Wilson) / - n : normal approximation propagation (see TEfficiency::Normal) / - ac : Agresti-Coull interval (see TEfficiency::AgrestiCoull) / - fc : Feldman-Cousins interval (see TEfficiency::FeldmanCousinsInterval) / - midp : Lancaster mid-P interval (see TEfficiency::MidPInterval) / - b(a,b): bayesian interval using a prior probability ~Beta(a,b); a,b > 0 / (see TEfficiency::Bayesian) / - mode : use mode of posterior for Bayesian interval (default is mean) / - shortest: use shortest interval (done by default if mode is set) / - central: use central interval (done by default if mode is NOT set) / - pois: interpret histograms as poisson ratio instead of efficiency / - e0 : plot efficiency and interval for bins where total=0 / (default is to skip them) / / Note: / Unfortunately there is no straightforward approach for determining a confidence / interval for a given confidence level. The actual coverage probability of the / confidence interval oscillates significantly according to the total number of / events and the true efficiency. In order to decrease the impact of this / oscillation on the actual coverage probability a couple of approximations and / methodes has been developed. For a detailed discussion, please have a look at / this statistical paper: / http://www-stat.wharton.upenn.edu/~tcai/paper/Binomial-StatSci.pdf
void TGraphAsymmErrors::Divide(const TH1* pass, const TH1* total, Option_t *opt) { check pointers if(!pass || !total) { Error("Divide","one of the passed pointers is zero"); return; }
check dimension of histograms; only 1-dimensional ones are accepted if((pass->GetDimension() > 1) || (total->GetDimension() > 1)) { Error("Divide","passed histograms are not one-dimensional"); return; }
check whether histograms are filled with weights -> use number of effective entries Bool_t bEffective = false; compare sum of weights with sum of squares of weights re-compute here to be sure to get the right values Double_t psumw = 0; Double_t psumw2 = 0; if (pass->GetSumw2()->fN > 0) { for (int i = 0; i < pass->GetNcells(); ++i) { psumw += pass->GetBinContent(i); psumw2 += pass->GetSumw2()->At(i); } } else { psumw = pass->GetSumOfWeights(); psumw2 = psumw; } if (TMath::Abs(psumw - psumw2) > 1e-6) bEffective = true;
Double_t tsumw = 0; Double_t tsumw2 = 0; if (total->GetSumw2()->fN > 0) { for (int i = 0; i < total->GetNcells(); ++i) { tsumw += total->GetBinContent(i); tsumw2 += total->GetSumw2()->At(i); } } else { tsumw = total->GetSumOfWeights(); tsumw2 = tsumw; } if (TMath::Abs(tsumw - tsumw2) > 1e-6) bEffective = true;
we do not want to ignore the weights if (bEffective && (pass->GetSumw2()->fN == 0 || total->GetSumw2()->fN == 0) ) { Warning("Divide","histogram have been computed with weights but the sum of weight squares are not stored in the histogram. Error calculation is performed ignoring the weights"); bEffective = false; }
parse option TString option = opt; option.ToLower();
Bool_t bVerbose = false; pointer to function returning the boundaries of the confidence interval (is only used in the frequentist cases.) Double_t (*pBound)(Double_t,Double_t,Double_t,Bool_t) = &TEfficiency::ClopperPearson; // default method confidence level Double_t conf = 0.682689492137; values for bayesian statistics Bool_t bIsBayesian = false; Double_t alpha = 1; Double_t beta = 1;
verbose mode if(option.Contains("v")) { option.ReplaceAll("v",""); bVerbose = true; if (bEffective) Info("Divide","weight will be considered in the Histogram Ratio"); }
confidence level if(option.Contains("cl=")) { Double_t level = -1; coverity [secure_coding : FALSE] sscanf(strstr(option.Data(),"cl="),"cl=%lf",&level); if((level > 0) && (level < 1)) conf = level; else Warning("Divide","given confidence level %.3lf is invalid",level); option.ReplaceAll("cl=",""); }
look for statistic options check first Bayesian case bayesian with prior Bool_t usePosteriorMode = false; Bool_t useShortestInterval = false; if (option.Contains("b(")) { Double_t a = 0; Double_t b = 0; sscanf(strstr(option.Data(), "b("), "b(%lf,%lf)", &a, &b); if (a > 0) alpha = a; else Warning("Divide", "given shape parameter for alpha %.2lf is invalid", a); if (b > 0) beta = b; else Warning("Divide", "given shape parameter for beta %.2lf is invalid", b); option.ReplaceAll("b(", ""); bIsBayesian = true;
look for specific bayesian options
use posterior mode
if (option.Contains("mode")) {
usePosteriorMode = true;
option.ReplaceAll("mode", "");
}
if (option.Contains("sh") || (usePosteriorMode && !option.Contains("cen"))) {
useShortestInterval = true;
}
} normal approximation else if (option.Contains("n")) { option.ReplaceAll("n", ""); pBound = &TEfficiency::Normal; } clopper pearson interval else if (option.Contains("cp")) { option.ReplaceAll("cp", ""); pBound = &TEfficiency::ClopperPearson; } wilson interval else if (option.Contains("w")) { option.ReplaceAll("w", ""); pBound = &TEfficiency::Wilson; } agresti coull interval else if (option.Contains("ac")) { option.ReplaceAll("ac", ""); pBound = &TEfficiency::AgrestiCoull; } Feldman-Cousins interval else if (option.Contains("fc")) { option.ReplaceAll("fc", ""); pBound = &TEfficiency::FeldmanCousins; } mid-P Lancaster interval else if (option.Contains("midp")) { option.ReplaceAll("midp", ""); pBound = &TEfficiency::MidPInterval; }
interpret as Poisson ratio Bool_t bPoissonRatio = false; if (option.Contains("pois")) { bPoissonRatio = true; option.ReplaceAll("pois", ""); } Bool_t plot0Bins = false; if (option.Contains("e0")) { plot0Bins = true; option.ReplaceAll("e0", ""); }
weights works only in case of Normal approximation or Bayesian for binomial interval in case of Poisson ratio we can use weights by rescaling the obtained results using the effective entries if ((bEffective && !bPoissonRatio) && !bIsBayesian && pBound != &TEfficiency::Normal) { Warning("Divide", "Histograms have weights: only Normal or Bayesian error calculation is supported"); Info("Divide", "Using now the Normal approximation for weighted histograms"); }
if (bPoissonRatio) { if (pass->GetDimension() != total->GetDimension()) { Error("Divide", "passed histograms are not of the same dimension"); return; }
if (!TEfficiency::CheckBinning(*pass, *total)) { Error("Divide", "passed histograms are not consistent"); return; } } else { check consistency of histograms, allowing weights if (!TEfficiency::CheckConsistency(*pass, *total, "w")) { Error("Divide", "passed histograms are not consistent"); return; } }
Set the graph to have a number of points equal to the number of histogram bins Int_t nbins = pass->GetNbinsX(); Set(nbins);
Ok, now set the points for each bin (Note: the TH1 bin content is shifted to the right by one: bin=0 is underflow, bin=nbins+1 is overflow.)
efficiency with lower and upper boundary of confidence interval double eff, low, upper; this keeps track of the number of points added to the graph int npoint=0; number of total and passed events Double_t t = 0 , p = 0; Double_t tw = 0, tw2 = 0, pw = 0, pw2 = 0, wratio = 1; // for the case of weights loop over all bins and fill the graph for (Int_t b=1; b<=nbins; ++b) {
default value when total =0; eff = 0; low = 0; upper = 0;
special case in case of weights we have to consider the sum of weights and the sum of weight squares if (bEffective) { tw = total->GetBinContent(b); tw2 = (total->GetSumw2()->fN > 0) ? total->GetSumw2()->At(b) : tw; pw = pass->GetBinContent(b); pw2 = (pass->GetSumw2()->fN > 0) ? pass->GetSumw2()->At(b) : pw;
if (bPoissonRatio) { tw += pw; tw2 += pw2; compute ratio on the effective entries ( p and t) special case is when (pw=0, pw2=0) in this case we cannot get the bin weight. we use then the overall weight of the full histogram if (pw == 0 && pw2 == 0) p = 0; else p = (pw * pw) / pw2;
if (tw == 0 && tw2 == 0) t = 0; else t = (tw * tw) / tw2;
if (pw > 0 && tw > 0) this is the ratio of the two bin weights ( pw/p / t/tw ) wratio = (pw * t) / (p * tw); else if (pw == 0 && tw > 0) case p histogram has zero compute the weights from all the histogram weight of histogram - sumw2/sumw wratio = (psumw2 * t) / (psumw * tw); else if (tw == 0 && pw > 0) case t histogram has zero compute the weights from all the histogram weight of histogram - sumw2/sumw wratio = (pw * tsumw) / (p * tsumw2); else if (p > 0) wratio = pw / p; // not sure if needed else { case both pw and tw are zero - we skip these bins if (!plot0Bins) continue; // skip bins with total <= 0 }
t += p; std::cout << p << " " << t << " " << wratio << std::endl; } else if (tw <= 0 && !plot0Bins) continue; // skip bins with total <= 0
in the case of weights have the formula only for the normal and bayesian statistics (see below)
}
use bin contents else { t = std::round(total->GetBinContent(b)); p = std::round(pass->GetBinContent(b));
if (bPoissonRatio) t += p;
if (t == 0.0 && !plot0Bins) continue; // skip bins with total = 0 }
using bayesian statistics if(bIsBayesian) { double aa,bb;
if ((bEffective && !bPoissonRatio) && tw2 <= 0) { case of bins with zero errors eff = pw/tw; low = eff; upper = eff; } else {
if (bEffective && !bPoissonRatio) { tw/tw2 re-normalize the weights double norm = tw/tw2; // case of tw2 = 0 is treated above aa = pw * norm + alpha; bb = (tw - pw) * norm + beta; } else { aa = double(p) + alpha; bb = double(t-p) + beta; } if (usePosteriorMode) eff = TEfficiency::BetaMode(aa,bb); else eff = TEfficiency::BetaMean(aa,bb);
if (useShortestInterval) { TEfficiency::BetaShortestInterval(conf,aa,bb,low,upper); } else { low = TEfficiency::BetaCentralInterval(conf,aa,bb,false); upper = TEfficiency::BetaCentralInterval(conf,aa,bb,true); } } } case of non-bayesian statistics else { if (bEffective && !bPoissonRatio) {
if (tw > 0) {
eff = pw/tw;
use normal error calculation using variance of MLE with weights (F.James 8.5.2) this is the same formula used in ROOT for TH1::Divide("B")
double variance = ( pw2 * (1. - 2 * eff) + tw2 * eff *eff ) / ( tw * tw) ;
double sigma = sqrt(variance);
double prob = 0.5 * (1.-conf);
double delta = ROOT::Math::normal_quantile_c(prob, sigma);
low = eff - delta;
upper = eff + delta;
if (low < 0) low = 0;
if (upper > 1) upper = 1.;
}
}
else {
when not using weights (all cases) or in case of Poisson ratio with weights if(t != 0.0) eff = ((Double_t)p)/t;
low = pBound(t,p,conf,false); upper = pBound(t,p,conf,true); } } treat as Poisson ratio if(bPoissonRatio) { Double_t ratio = eff/(1 - eff); take the intervals in eff as intervals in the Poisson ratio low = low/(1. - low); upper = upper/(1.-upper); eff = ratio; if (bEffective) { scale result by the ratio of the weight eff *= wratio; low *= wratio; upper *= wratio; } } Set the point center and its errors if (TMath::Finite(eff)) { SetPoint(npoint,pass->GetBinCenter(b),eff); SetPointError(npoint, pass->GetBinCenter(b)-pass->GetBinLowEdge(b), pass->GetBinLowEdge(b)-pass->GetBinCenter(b)+pass->GetBinWidth(b), eff-low,upper-eff); npoint++;//we have added a point to the graph } }
Set(npoint);//tell the graph how many points we've really added if (npoint < nbins) Warning("Divide","Number of graph points is different than histogram bins - %d points have been skipped",nbins-npoint);
if (bVerbose) { Info("Divide","made a graph with %d points from %d bins",npoint,nbins); Info("Divide","used confidence level: %.2lf\n",conf); if(bIsBayesian) Info("Divide","used prior probability ~ beta(%.2lf,%.2lf)",alpha,beta); Print(); } }
////////////////////////////////////////////////////////////////////////////// / Compute Range.
void TGraphAsymmErrors::ComputeRange(Double_t &xmin, Double_t &ymin, Double_t &xmax, Double_t &ymax) const { TGraph::ComputeRange(xmin,ymin,xmax,ymax);
for (Int_t i=0;i<fNpoints;i++) { if (fX[i] -fEXlow[i] < xmin) { if ( (TVirtualPad::Pad()) && (TVirtualPad::Pad()) ->GetLogx()) { if (fEXlow[i] < fX[i]) xmin = fX[i]-fEXlow[i]; else xmin = TMath::Min(xmin,fX[i]/3); } else { xmin = fX[i]-fEXlow[i]; } } if (fX[i] +fEXhigh[i] > xmax) xmax = fX[i]+fEXhigh[i]; if (fY[i] -fEYlow[i] < ymin) { if ( (TVirtualPad::Pad()) && (TVirtualPad::Pad()) ->GetLogy()) { if (fEYlow[i] < fY[i]) ymin = fY[i]-fEYlow[i]; else ymin = TMath::Min(ymin,fY[i]/3); } else { ymin = fY[i]-fEYlow[i]; } } if (fY[i] +fEYhigh[i] > ymax) ymax = fY[i]+fEYhigh[i]; } }
////////////////////////////////////////////////////////////////////////////// / Copy and release.
void TGraphAsymmErrors::CopyAndRelease(Double_t **newarrays, Int_t ibegin, Int_t iend, Int_t obegin) { CopyPoints(newarrays, ibegin, iend, obegin); if (newarrays) { delete[] fEXlow; fEXlow = newarrays[0]; delete[] fEXhigh; fEXhigh = newarrays[1]; delete[] fEYlow; fEYlow = newarrays[2]; delete[] fEYhigh; fEYhigh = newarrays[3]; delete[] fX; fX = newarrays[4]; delete[] fY; fY = newarrays[5]; delete[] newarrays; } }
////////////////////////////////////////////////////////////////////////////// / Copy errors from fE*** to arrays[***] / or to f*** Copy points.
Bool_t TGraphAsymmErrors::CopyPoints(Double_t **arrays, Int_t ibegin, Int_t iend, Int_t obegin) { if (TGraph::CopyPoints(arrays ? arrays+4 : nullptr, ibegin, iend, obegin)) { Int_t n = (iend - ibegin)*sizeof(Double_t); if (arrays) { memmove(&arrays[0][obegin], &fEXlow[ibegin], n); memmove(&arrays[1][obegin], &fEXhigh[ibegin], n); memmove(&arrays[2][obegin], &fEYlow[ibegin], n); memmove(&arrays[3][obegin], &fEYhigh[ibegin], n); } else { memmove(&fEXlow[obegin], &fEXlow[ibegin], n); memmove(&fEXhigh[obegin], &fEXhigh[ibegin], n); memmove(&fEYlow[obegin], &fEYlow[ibegin], n); memmove(&fEYhigh[obegin], &fEYhigh[ibegin], n); } return kTRUE; } else { return kFALSE; } }
////////////////////////////////////////////////////////////////////////////// / Should be called from ctors after fNpoints has been set. / Note: This function should be called only from the constructor / since it does not delete previously existing arrays
Bool_t TGraphAsymmErrors::CtorAllocate() { if (!fNpoints) { fEXlow = fEYlow = fEXhigh = fEYhigh = nullptr; return kFALSE; } fEXlow = new Double_t[fMaxSize]; fEYlow = new Double_t[fMaxSize]; fEXhigh = new Double_t[fMaxSize]; fEYhigh = new Double_t[fMaxSize]; return kTRUE; }
////////////////////////////////////////////////////////////////////////////// / Protected function to perform the merge operation of a graph with asymmetric errors.
Bool_t TGraphAsymmErrors::DoMerge(const TGraph *g) { if (g->GetN() == 0) return kFALSE;
Double_t * exl = g->GetEXlow(); Double_t * exh = g->GetEXhigh(); Double_t * eyl = g->GetEYlow(); Double_t * eyh = g->GetEYhigh(); if (exl == nullptr || exh == nullptr || eyl == nullptr || eyh == nullptr) { if (g->IsA() != TGraph::Class() ) Warning("DoMerge","Merging a %s is not compatible with a TGraphAsymmErrors - errors will be ignored",g->IsA()->GetName()); return TGraph::DoMerge(g); } for (Int_t i = 0 ; i < g->GetN(); i++) { Int_t ipoint = GetN(); Double_t x = g->GetX()[i]; Double_t y = g->GetY()[i]; SetPoint(ipoint, x, y); SetPointError(ipoint, exl[i], exh[i], eyl[i], eyh[i] ); }
return kTRUE; }
////////////////////////////////////////////////////////////////////////////// / Set zero values for point arrays in the range [begin, end]
void TGraphAsymmErrors::FillZero(Int_t begin, Int_t end, Bool_t from_ctor) { if (!from_ctor) { TGraph::FillZero(begin, end, from_ctor); } Int_t n = (end - begin)*sizeof(Double_t); memset(fEXlow + begin, 0, n); memset(fEXhigh + begin, 0, n); memset(fEYlow + begin, 0, n); memset(fEYhigh + begin, 0, n); }
////////////////////////////////////////////////////////////////////////////// / Returns the combined error along X at point i by computing the average / of the lower and upper variance.
Double_t TGraphAsymmErrors::GetErrorX(Int_t i) const { if (i < 0 || i >= fNpoints) return -1; if (!fEXlow && !fEXhigh) return -1; Double_t elow=0, ehigh=0; if (fEXlow) elow = fEXlow[i]; if (fEXhigh) ehigh = fEXhigh[i]; return TMath::Sqrt(0.5*(elow*elow + ehigh*ehigh)); }
////////////////////////////////////////////////////////////////////////////// / Returns the combined error along Y at point i by computing the average / of the lower and upper variance.
Double_t TGraphAsymmErrors::GetErrorY(Int_t i) const { if (i < 0 || i >= fNpoints) return -1; if (!fEYlow && !fEYhigh) return -1; Double_t elow=0, ehigh=0; if (fEYlow) elow = fEYlow[i]; if (fEYhigh) ehigh = fEYhigh[i]; return TMath::Sqrt(0.5*(elow*elow + ehigh*ehigh)); }
////////////////////////////////////////////////////////////////////////////// / Get high error on X.
Double_t TGraphAsymmErrors::GetErrorXhigh(Int_t i) const { if (i<0 || i>fNpoints) return -1; if (fEXhigh) return fEXhigh[i]; return -1; }
////////////////////////////////////////////////////////////////////////////// / Get low error on X.
Double_t TGraphAsymmErrors::GetErrorXlow(Int_t i) const { if (i<0 || i>fNpoints) return -1; if (fEXlow) return fEXlow[i]; return -1; }
////////////////////////////////////////////////////////////////////////////// / Get high error on Y.
Double_t TGraphAsymmErrors::GetErrorYhigh(Int_t i) const { if (i<0 || i>fNpoints) return -1; if (fEYhigh) return fEYhigh[i]; return -1; }
////////////////////////////////////////////////////////////////////////////// / Get low error on Y.
Double_t TGraphAsymmErrors::GetErrorYlow(Int_t i) const { if (i<0 || i>fNpoints) return -1; if (fEYlow) return fEYlow[i]; return -1; }
////////////////////////////////////////////////////////////////////////////// / Adds all graphs with asymmetric errors from the collection to this graph. / Returns the total number of points in the result or -1 in case of an error.
Int_t TGraphAsymmErrors::Merge(TCollection* li) { TIter next(li); while (TObject* o = next()) { TGraph g = dynamic_cast<TGraph>(o); if (!g) { Error("Merge", "Cannot merge - an object which doesn't inherit from TGraph found in the list"); return -1; } int n0 = GetN(); int n1 = n0+g->GetN(); Set(n1); Double_t * x = g->GetX(); Double_t * y = g->GetY(); Double_t * exlow = g->GetEXlow(); Double_t * exhigh = g->GetEXhigh(); Double_t * eylow = g->GetEYlow(); Double_t * eyhigh = g->GetEYhigh(); for (Int_t i = 0 ; i < g->GetN(); i++) { SetPoint(n0+i, x[i], y[i]); if (exlow) fEXlow[n0+i] = exlow[i]; if (exhigh) fEXhigh[n0+i] = exhigh[i]; if (eylow) fEYlow[n0+i] = eylow[i]; if (eyhigh) fEYhigh[n0+i] = eyhigh[i]; } } return GetN(); }
////////////////////////////////////////////////////////////////////////////// / Print graph and errors values.
void TGraphAsymmErrors::Print(Option_t *) const { for (Int_t i=0;i<fNpoints;i++) { printf("x[%d]=%g, y[%d]=%g, exl[%d]=%g, exh[%d]=%g, eyl[%d]=%g, eyh[%d]=%g\n" ,i,fX[i],i,fY[i],i,fEXlow[i],i,fEXhigh[i],i,fEYlow[i],i,fEYhigh[i]); } }
////////////////////////////////////////////////////////////////////////////// / Save primitive as a C++ statement(s) on output stream out.
void TGraphAsymmErrors::SavePrimitive(std::ostream &out, Option_t option /= ""
Definition at line 1239 of file TGraphAsymmErrors.cxx.
| auto fXName = SaveArray(out, "fx", frameNumber, fX) |
Definition at line 1242 of file TGraphAsymmErrors.cxx.
| auto fYName = SaveArray(out, "fy", frameNumber, fY) |
Definition at line 1243 of file TGraphAsymmErrors.cxx.
| Int_t i |
Definition at line 1306 of file TGraphAsymmErrors.cxx.
| Int_t ipoint = -2 |
Definition at line 1305 of file TGraphAsymmErrors.cxx.
| Int_t py = (TVirtualPad::Pad()) ->GetEventY() |
Definition at line 1302 of file TGraphAsymmErrors.cxx.
