Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TGraphAsymmErrors.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: Rene Brun 03/03/99
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12
13#include "TEfficiency.h"
14#include "TROOT.h"
15#include "TBuffer.h"
16#include "TGraphAsymmErrors.h"
17#include "TGraphErrors.h"
18#include "TStyle.h"
19#include "TMath.h"
20#include "TVirtualPad.h"
21#include "TF1.h"
22#include "TH1.h"
23#include "TVector.h"
24#include "TVectorD.h"
25#include "TSystem.h"
27#include "strtok.h"
28
29#include <cstring>
30#include <iostream>
31#include <fstream>
32
33
35
36/** \class TGraphAsymmErrors
37 \ingroup Graphs
38TGraph with asymmetric error bars.
39
40The TGraphAsymmErrors painting is performed thanks to the TGraphPainter
41class. All details about the various painting options are given in this class.
42
43The picture below gives an example:
44
45Begin_Macro(source)
46{
47 auto c1 = new TCanvas("c1","A Simple Graph with asymmetric error bars",200,10,700,500);
48 c1->SetFillColor(42);
49 c1->SetGrid();
50 c1->GetFrame()->SetFillColor(21);
51 c1->GetFrame()->SetBorderSize(12);
52 const Int_t n = 10;
53 Double_t x[n] = {-0.22, 0.05, 0.25, 0.35, 0.5, 0.61,0.7,0.85,0.89,0.95};
54 Double_t y[n] = {1,2.9,5.6,7.4,9,9.6,8.7,6.3,4.5,1};
55 Double_t exl[n] = {.05,.1,.07,.07,.04,.05,.06,.07,.08,.05};
56 Double_t eyl[n] = {.8,.7,.6,.5,.4,.4,.5,.6,.7,.8};
57 Double_t exh[n] = {.02,.08,.05,.05,.03,.03,.04,.05,.06,.03};
58 Double_t eyh[n] = {.6,.5,.4,.3,.2,.2,.3,.4,.5,.6};
59 auto gr = new TGraphAsymmErrors(n,x,y,exl,exh,eyl,eyh);
60 gr->SetTitle("TGraphAsymmErrors Example");
61 gr->SetMarkerColor(4);
62 gr->SetMarkerStyle(21);
63 gr->Draw("ALP");
64}
65End_Macro
66*/
67
68
69////////////////////////////////////////////////////////////////////////////////
70/// TGraphAsymmErrors default constructor.
71
73{
74}
75
76
77////////////////////////////////////////////////////////////////////////////////
78/// TGraphAsymmErrors copy constructor
79
81 : TGraph(gr)
82{
83 if (!CtorAllocate()) return;
84 Int_t n = fNpoints*sizeof(Double_t);
85 memcpy(fEXlow, gr.fEXlow, n);
86 memcpy(fEYlow, gr.fEYlow, n);
87 memcpy(fEXhigh, gr.fEXhigh, n);
88 memcpy(fEYhigh, gr.fEYhigh, n);
89}
90
91
92////////////////////////////////////////////////////////////////////////////////
93/// TGraphAsymmErrors assignment operator.
94
96{
97 if(this!=&gr) {
99 // delete arrays
100 if (fEXlow) delete [] fEXlow;
101 if (fEYlow) delete [] fEYlow;
102 if (fEXhigh) delete [] fEXhigh;
103 if (fEYhigh) delete [] fEYhigh;
104
105 if (!CtorAllocate()) return *this;
106 Int_t n = fNpoints*sizeof(Double_t);
107 memcpy(fEXlow, gr.fEXlow, n);
108 memcpy(fEYlow, gr.fEYlow, n);
109 memcpy(fEXhigh, gr.fEXhigh, n);
110 memcpy(fEYhigh, gr.fEYhigh, n);
111 }
112 return *this;
113}
114
115
116////////////////////////////////////////////////////////////////////////////////
117/// TGraphAsymmErrors normal constructor.
118///
119/// the arrays are preset to zero
120
122 : TGraph(n)
123{
124 if (!CtorAllocate()) return;
125 FillZero(0, fNpoints);
126}
127
128
129////////////////////////////////////////////////////////////////////////////////
130/// TGraphAsymmErrors normal constructor.
131///
132/// if exl,h or eyl,h are null, the corresponding arrays are preset to zero
133
134TGraphAsymmErrors::TGraphAsymmErrors(Int_t n, const Float_t *x, const Float_t *y, const Float_t *exl, const Float_t *exh, const Float_t *eyl, const Float_t *eyh)
135 : TGraph(n,x,y)
136{
137 if (!CtorAllocate()) return;
138
139 for (Int_t i=0;i<n;i++) {
140 if (exl) fEXlow[i] = exl[i];
141 else fEXlow[i] = 0;
142 if (exh) fEXhigh[i] = exh[i];
143 else fEXhigh[i] = 0;
144 if (eyl) fEYlow[i] = eyl[i];
145 else fEYlow[i] = 0;
146 if (eyh) fEYhigh[i] = eyh[i];
147 else fEYhigh[i] = 0;
148 }
149}
150
151
152////////////////////////////////////////////////////////////////////////////////
153/// TGraphAsymmErrors normal constructor.
154///
155/// if exl,h or eyl,h are null, the corresponding arrays are preset to zero
156
157TGraphAsymmErrors::TGraphAsymmErrors(Int_t n, const Double_t *x, const Double_t *y, const Double_t *exl, const Double_t *exh, const Double_t *eyl, const Double_t *eyh)
158 : TGraph(n,x,y)
159{
160 if (!CtorAllocate()) return;
161
162 n = fNpoints*sizeof(Double_t);
163 if (exl)
164 memcpy(fEXlow, exl, n);
165 else
166 memset(fEXlow, 0, n);
167 if (exh)
168 memcpy(fEXhigh, exh, n);
169 else
170 memset(fEXhigh, 0, n);
171 if (eyl)
172 memcpy(fEYlow, eyl, n);
173 else
174 memset(fEYlow, 0, n);
175 if (eyh)
176 memcpy(fEYhigh, eyh, n);
177 else
178 memset(fEYhigh, 0, n);
179}
180
181
182////////////////////////////////////////////////////////////////////////////////
183/// Constructor with six vectors of floats in input
184/// A grapherrors is built with the X coordinates taken from vx and Y coord from vy
185/// and the errors from vectors vexl/h and veyl/h.
186/// The number of points in the graph is the minimum of number of points
187/// in vx and vy.
188
189TGraphAsymmErrors::TGraphAsymmErrors(const TVectorF &vx, const TVectorF &vy, const TVectorF &vexl, const TVectorF &vexh, const TVectorF &veyl, const TVectorF &veyh)
190 :TGraph()
191{
192 fNpoints = TMath::Min(vx.GetNrows(), vy.GetNrows());
193 if (!TGraph::CtorAllocate()) return;
194 if (!CtorAllocate()) return;
195 Int_t ivxlow = vx.GetLwb();
196 Int_t ivylow = vy.GetLwb();
197 Int_t ivexllow = vexl.GetLwb();
198 Int_t ivexhlow = vexh.GetLwb();
199 Int_t iveyllow = veyl.GetLwb();
200 Int_t iveyhlow = veyh.GetLwb();
201 for (Int_t i=0;i<fNpoints;i++) {
202 fX[i] = vx(i+ivxlow);
203 fY[i] = vy(i+ivylow);
204 fEXlow[i] = vexl(i+ivexllow);
205 fEYlow[i] = veyl(i+iveyllow);
206 fEXhigh[i] = vexh(i+ivexhlow);
207 fEYhigh[i] = veyh(i+iveyhlow);
208 }
209}
210
211
212////////////////////////////////////////////////////////////////////////////////
213/// Constructor with six vectors of doubles in input
214/// A grapherrors is built with the X coordinates taken from vx and Y coord from vy
215/// and the errors from vectors vexl/h and veyl/h.
216/// The number of points in the graph is the minimum of number of points
217/// in vx and vy.
218
219TGraphAsymmErrors::TGraphAsymmErrors(const TVectorD &vx, const TVectorD &vy, const TVectorD &vexl, const TVectorD &vexh, const TVectorD &veyl, const TVectorD &veyh)
220 :TGraph()
221{
222 fNpoints = TMath::Min(vx.GetNrows(), vy.GetNrows());
223 if (!TGraph::CtorAllocate()) return;
224 if (!CtorAllocate()) return;
225 Int_t ivxlow = vx.GetLwb();
226 Int_t ivylow = vy.GetLwb();
227 Int_t ivexllow = vexl.GetLwb();
228 Int_t ivexhlow = vexh.GetLwb();
229 Int_t iveyllow = veyl.GetLwb();
230 Int_t iveyhlow = veyh.GetLwb();
231 for (Int_t i=0;i<fNpoints;i++) {
232 fX[i] = vx(i+ivxlow);
233 fY[i] = vy(i+ivylow);
234 fEXlow[i] = vexl(i+ivexllow);
235 fEYlow[i] = veyl(i+iveyllow);
236 fEXhigh[i] = vexh(i+ivexhlow);
237 fEYhigh[i] = veyh(i+iveyhlow);
238 }
239}
240
241
242////////////////////////////////////////////////////////////////////////////////
243/// TGraphAsymmErrors constructor importing its parameters from the TH1 object passed as argument
244/// the low and high errors are set to the bin error of the histogram.
245
247 : TGraph(h)
248{
249 if (!CtorAllocate()) return;
250
251 for (Int_t i = 0; i < fNpoints; i++) {
252 fEXlow[i] = h->GetBinWidth(i+1)*gStyle->GetErrorX();
253 fEXhigh[i] = fEXlow[i];
254 fEYlow[i] = h->GetBinErrorLow(i+1);
255 fEYhigh[i] = h->GetBinErrorUp(i+1);;
256 }
257}
258
259
260////////////////////////////////////////////////////////////////////////////////
261/// Creates a TGraphAsymmErrors by dividing two input TH1 histograms:
262/// pass/total. (see TGraphAsymmErrors::Divide)
263
265 : TGraph((pass)?pass->GetNbinsX():0)
266{
267 if (!pass || !total) {
268 Error("TGraphAsymmErrors","Invalid histogram pointers");
269 return;
270 }
271 if (!CtorAllocate()) return;
272
273 std::string sname = "divide_" + std::string(pass->GetName()) + "_by_" +
274 std::string(total->GetName());
275 SetName(sname.c_str());
276 SetTitle(pass->GetTitle());
277
278 //copy style from pass
279 pass->TAttLine::Copy(*this);
280 pass->TAttFill::Copy(*this);
281 pass->TAttMarker::Copy(*this);
282
283 Divide(pass, total, option);
284}
285
286
287////////////////////////////////////////////////////////////////////////////////
288/// TGraphAsymmErrors constructor reading input from filename
289/// filename is assumed to contain at least 2 columns of numbers
290///
291/// convention for format (default=`"%lg %lg %lg %lg %lg %lg"`)
292/// - format = `"%lg %lg"` read only 2 first columns into X, Y
293/// - format = `"%lg %lg %lg %lg"` read only 4 first columns into X, Y, ELY, EHY
294/// - format = `"%lg %lg %lg %lg %lg %lg"` read only 6 first columns into X, Y, EXL, EYH, EYL, EHY
295///
296/// For files separated by a specific delimiter different from `' '` and `'\\t'` (e.g. `';'` in csv files)
297/// you can avoid using `%*s` to bypass this delimiter by explicitly specify the `"option" argument,
298/// e.g. `option=" \\t,;"` for columns of figures separated by any of these characters `(' ', '\\t', ',', ';')`
299/// used once `(e.g. "1;1")` or in a combined way `(" 1;,;; 1")`.
300/// Note in that case, the instantiation is about 2 times slower.
301/// In case a delimiter is specified, the format `"%lg %lg %lg"` will read X,Y,EX.
302
304 : TGraph(100)
305{
306 if (!CtorAllocate()) return;
307 Double_t x, y, exl, exh, eyl, eyh;
308 TString fname = filename;
309 gSystem->ExpandPathName(fname);
310 std::ifstream infile(fname.Data());
311 if (!infile.good()) {
312 MakeZombie();
313 Error("TGraphAsymmErrors", "Cannot open file: %s, TGraphAsymmErrors is Zombie", filename);
314 fNpoints = 0;
315 return;
316 }
317 std::string line;
318 Int_t np = 0;
319
320 if (strcmp(option, "") == 0) { // No delimiters specified (standard constructor).
321
322 Int_t ncol = TGraphErrors::CalculateScanfFields(format); //count number of columns in format
323 Int_t res;
324 while (std::getline(infile, line, '\n')) {
325 exl = exh = eyl = eyh = 0.;
326 if (ncol < 3) {
327 res = sscanf(line.c_str(), format, &x, &y);
328 } else if (ncol < 5) {
329 res = sscanf(line.c_str(), format, &x, &y, &eyl, &eyh);
330 } else {
331 res = sscanf(line.c_str(), format, &x, &y, &exl, &exh, &eyl, &eyh);
332 }
333 if (res < 2) {
334 continue; //skip empty and ill-formed lines
335 }
336 SetPoint(np, x, y);
337 SetPointError(np, exl, exh, eyl, eyh);
338 np++;
339 }
340 Set(np);
341
342 } else { // A delimiter has been specified in "option"
343
344 // Checking format and creating its boolean equivalent
345 TString format_ = TString(format) ;
346 format_.ReplaceAll(" ", "") ;
347 format_.ReplaceAll("\t", "") ;
348 format_.ReplaceAll("lg", "") ;
349 format_.ReplaceAll("s", "") ;
350 format_.ReplaceAll("%*", "0") ;
351 format_.ReplaceAll("%", "1") ;
352 if (!format_.IsDigit()) {
353 Error("TGraphAsymmErrors", "Incorrect input format! Allowed format tags are {\"%%lg\",\"%%*lg\" or \"%%*s\"}");
354 return ;
355 }
356 Int_t ntokens = format_.Length() ;
357 if (ntokens < 2) {
358 Error("TGraphAsymmErrors", "Incorrect input format! Only %d tag(s) in format whereas at least 2 \"%%lg\" tags are expected!", ntokens);
359 return ;
360 }
361 Int_t ntokensToBeSaved = 0;
362 Bool_t * isTokenToBeSaved = new Bool_t[ntokens];
363 for (Int_t idx = 0; idx < ntokens; idx++) {
364 isTokenToBeSaved[idx] = TString::Format("%c", format_[idx]).Atoi(); //atoi(&format_[idx]) does not work for some reason...
365 if (isTokenToBeSaved[idx] == 1) {
366 ntokensToBeSaved++ ;
367 }
368 }
369 if (ntokens >= 2 && (ntokensToBeSaved < 2 || ntokensToBeSaved > 4)) { //first condition not to repeat the previous error message
370 Error("TGraphAsymmErrors", "Incorrect input format! There are %d \"%%lg\" tag(s) in format whereas 2,3 or 4 are expected!", ntokensToBeSaved);
371 delete [] isTokenToBeSaved;
372 return ;
373 }
374
375 // Initializing loop variables
376 Bool_t isLineToBeSkipped = kFALSE; //empty and ill-formed lines
377 char *token = nullptr;
378 TString token_str = "";
379 Int_t token_idx = 0;
380 Double_t value[6]; //x,y,exl, exh, eyl, eyh buffers
381 for (Int_t k = 0; k < 6; k++)
382 value[k] = 0.;
383 Int_t value_idx = 0;
384
385 // Looping
386 char *rest;
387 while (std::getline(infile, line, '\n')) {
388 if (!line.empty()) {
389 if (line[line.size() - 1] == char(13)) { // removing DOS CR character
390 line.erase(line.end() - 1, line.end()) ;
391 }
392 token = R__STRTOK_R(const_cast<char*>(line.c_str()), option, &rest) ;
393 while (token != nullptr && value_idx < ntokensToBeSaved) {
394 if (isTokenToBeSaved[token_idx]) {
395 token_str = TString(token) ;
396 token_str.ReplaceAll("\t", "") ;
397 if (!token_str.IsFloat()) {
398 isLineToBeSkipped = kTRUE ;
399 break ;
400 } else {
401 value[value_idx] = token_str.Atof() ;
402 value_idx++ ;
403 }
404 }
405 token = R__STRTOK_R(nullptr, option, &rest); // next token
406 token_idx++ ;
407 }
408 if (!isLineToBeSkipped && value_idx > 1) { //i.e. 2,3 or 4
409 x = value[0];
410 y = value[1];
411 exl = value[2];
412 exh = value[3];
413 eyl = value[4];
414 eyh = value[5];
415 SetPoint(np, x, y);
416 SetPointError(np, exl, exh, eyl, eyh);
417 np++ ;
418 }
419 }
420 isLineToBeSkipped = kFALSE;
421 token = nullptr;
422 token_idx = 0;
423 value_idx = 0;
424 }
425 Set(np) ;
426
427 // Cleaning
428 delete [] isTokenToBeSaved;
429 delete token;
430 }
431 infile.close();
432}
433
434////////////////////////////////////////////////////////////////////////////////
435/// TGraphAsymmErrors default destructor.
436
438{
439 if(fEXlow) delete [] fEXlow;
440 if(fEXhigh) delete [] fEXhigh;
441 if(fEYlow) delete [] fEYlow;
442 if(fEYhigh) delete [] fEYhigh;
443}
444
445////////////////////////////////////////////////////////////////////////////////
446/// Allocate internal data structures for `size` points.
447
449 return AllocateArrays(6, size);
450}
451
452////////////////////////////////////////////////////////////////////////////////
453/// Apply a function to all data points \f$ y = f(x,y) \f$
454///
455/// Errors are calculated as \f$ eyh = f(x,y+eyh)-f(x,y) \f$ and
456/// \f$ eyl = f(x,y)-f(x,y-eyl) \f$
457///
458/// Special treatment has to be applied for the functions where the
459/// role of "up" and "down" is reversed.
460///
461/// Function suggested/implemented by Miroslav Helbich <helbich@mail.desy.de>
462
464{
465 Double_t x,y,exl,exh,eyl,eyh,eyl_new,eyh_new,fxy;
466
467 if (fHistogram) {
468 delete fHistogram;
469 fHistogram = nullptr;
470 }
471 for (Int_t i=0;i<GetN();i++) {
472 GetPoint(i,x,y);
473 exl = GetErrorXlow(i);
474 exh = GetErrorXhigh(i);
475 eyl = GetErrorYlow(i);
476 eyh = GetErrorYhigh(i);
477
478 fxy = f->Eval(x,y);
479 SetPoint(i,x,fxy);
480
481 // in the case of the functions like y-> -1*y the roles of the
482 // upper and lower error bars is reversed
483 if (f->Eval(x,y-eyl)<f->Eval(x,y+eyh)) {
484 eyl_new = TMath::Abs(fxy - f->Eval(x,y-eyl));
485 eyh_new = TMath::Abs(f->Eval(x,y+eyh) - fxy);
486 } else {
487 eyh_new = TMath::Abs(fxy - f->Eval(x,y-eyl));
488 eyl_new = TMath::Abs(f->Eval(x,y+eyh) - fxy);
489 }
490
491 //error on x doesn't change
492 SetPointError(i,exl,exh,eyl_new,eyh_new);
493 }
494 if (gPad) gPad->Modified();
495}
496
497////////////////////////////////////////////////////////////////////////////////
498///This function is only kept for backward compatibility.
499///You should rather use the Divide method.
500///It calls `Divide(pass,total,"cl=0.683 b(1,1) mode")` which is equivalent to the
501///former BayesDivide method.
502
504{
505 Divide(pass,total,"cl=0.683 b(1,1) mode");
506}
507
508////////////////////////////////////////////////////////////////////////////////
509/// Fill this TGraphAsymmErrors by dividing two 1-dimensional histograms pass/total
510///
511/// This method serves two purposes:
512///
513/// ### 1) calculating efficiencies:
514///
515/// The assumption is that the entries in "pass" are a subset of those in
516/// "total". That is, we create an "efficiency" graph, where each entry is
517/// between 0 and 1, inclusive.
518///
519/// If the histograms are not filled with unit weights, the number of effective
520/// entries is used to normalise the bin contents which might lead to wrong results.
521/// \f[
522/// \text{effective entries} = \frac{(\sum w_{i})^{2}}{\sum w_{i}^{2}}
523/// \f]
524/// The points are assigned a x value at the center of each histogram bin.
525/// The y values are \f$\text{eff} = \frac{\text{pass}}{\text{total}}\f$
526/// for all options except for the
527/// bayesian methods where the result depends on the chosen option.
528///
529/// If the denominator becomes 0 or pass > total, the corresponding bin is
530/// skipped.
531///
532/// ### 2) calculating ratios of two Poisson means (option 'pois'):
533///
534/// The two histograms are interpreted as independent Poisson processes and the ratio
535/// \f[
536/// \tau = \frac{n_{1}}{n_{2}} = \frac{\varepsilon}{1 - \varepsilon}
537/// \f]
538/// with \f$\varepsilon = \frac{n_{1}}{n_{1} + n_{2}}\f$.
539/// The histogram 'pass' is interpreted as \f$n_{1}\f$ and the total histogram
540/// is used for \f$n_{2}\f$.
541///
542/// The (asymmetric) uncertainties of the Poisson ratio are linked to the uncertainties
543/// of efficiency by a parameter transformation:
544/// \f[
545/// \Delta \tau_{low/up} = \frac{1}{(1 - \varepsilon)^{2}} \Delta \varepsilon_{low/up}
546/// \f]
547///
548/// The x errors span each histogram bin (lowedge ... lowedge+width)
549/// The y errors depend on the chosen statistic methode which can be determined
550/// by the options given below. For a detailed description of the used statistic
551/// calculations please have a look at the corresponding functions!
552///
553/// Options:
554/// - v : verbose mode: prints information about the number of used bins
555/// and calculated efficiencies with their errors
556/// - cl=x : determine the used confidence level (0<x<1) (default is 0.683)
557/// - cp : Clopper-Pearson interval (see TEfficiency::ClopperPearson)
558/// - w : Wilson interval (see TEfficiency::Wilson)
559/// - n : normal approximation propagation (see TEfficiency::Normal)
560/// - ac : Agresti-Coull interval (see TEfficiency::AgrestiCoull)
561/// - fc : Feldman-Cousins interval (see TEfficiency::FeldmanCousinsInterval)
562/// - midp : Lancaster mid-P interval (see TEfficiency::MidPInterval)
563/// - b(a,b): bayesian interval using a prior probability ~Beta(a,b); a,b > 0
564/// (see TEfficiency::Bayesian)
565/// - mode : use mode of posterior for Bayesian interval (default is mean)
566/// - shortest: use shortest interval (done by default if mode is set)
567/// - central: use central interval (done by default if mode is NOT set)
568/// - pois: interpret histograms as poisson ratio instead of efficiency
569/// - e0 : plot efficiency and interval for bins where total=0
570/// (default is to skip them)
571///
572/// Note:
573/// Unfortunately there is no straightforward approach for determining a confidence
574/// interval for a given confidence level. The actual coverage probability of the
575/// confidence interval oscillates significantly according to the total number of
576/// events and the true efficiency. In order to decrease the impact of this
577/// oscillation on the actual coverage probability a couple of approximations and
578/// methodes has been developed. For a detailed discussion, please have a look at
579/// this statistical paper:
580/// http://www-stat.wharton.upenn.edu/~tcai/paper/Binomial-StatSci.pdf
581
582
583void TGraphAsymmErrors::Divide(const TH1* pass, const TH1* total, Option_t *opt)
584{
585 //check pointers
586 if(!pass || !total) {
587 Error("Divide","one of the passed pointers is zero");
588 return;
589 }
590
591 //check dimension of histograms; only 1-dimensional ones are accepted
592 if((pass->GetDimension() > 1) || (total->GetDimension() > 1)) {
593 Error("Divide","passed histograms are not one-dimensional");
594 return;
595 }
596
597 //check whether histograms are filled with weights -> use number of effective
598 //entries
599 Bool_t bEffective = false;
600 //compare sum of weights with sum of squares of weights
601 // re-compute here to be sure to get the right values
602 Double_t psumw = 0;
603 Double_t psumw2 = 0;
604 if (pass->GetSumw2()->fN > 0) {
605 for (int i = 0; i < pass->GetNcells(); ++i) {
606 psumw += pass->GetBinContent(i);
607 psumw2 += pass->GetSumw2()->At(i);
608 }
609 }
610 else {
611 psumw = pass->GetSumOfWeights();
612 psumw2 = psumw;
613 }
614 if (TMath::Abs(psumw - psumw2) > 1e-6)
615 bEffective = true;
616
617 Double_t tsumw = 0;
618 Double_t tsumw2 = 0;
619 if (total->GetSumw2()->fN > 0) {
620 for (int i = 0; i < total->GetNcells(); ++i) {
621 tsumw += total->GetBinContent(i);
622 tsumw2 += total->GetSumw2()->At(i);
623 }
624 }
625 else {
626 tsumw = total->GetSumOfWeights();
627 tsumw2 = tsumw;
628 }
629 if (TMath::Abs(tsumw - tsumw2) > 1e-6)
630 bEffective = true;
631
632
633
634 // we do not want to ignore the weights
635 // if (bEffective && (pass->GetSumw2()->fN == 0 || total->GetSumw2()->fN == 0) ) {
636 // 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");
637 // bEffective = false;
638 // }
639
640 //parse option
641 TString option = opt;
642 option.ToLower();
643
644 Bool_t bVerbose = false;
645 //pointer to function returning the boundaries of the confidence interval
646 //(is only used in the frequentist cases.)
648 //confidence level
649 Double_t conf = 0.682689492137;
650 //values for bayesian statistics
651 Bool_t bIsBayesian = false;
652 Double_t alpha = 1;
653 Double_t beta = 1;
654
655 //verbose mode
656 if(option.Contains("v")) {
657 option.ReplaceAll("v","");
658 bVerbose = true;
659 if (bEffective)
660 Info("Divide","weight will be considered in the Histogram Ratio");
661 }
662
663
664 //confidence level
665 if(option.Contains("cl=")) {
666 Double_t level = -1;
667 // coverity [secure_coding : FALSE]
668 sscanf(strstr(option.Data(),"cl="),"cl=%lf",&level);
669 if((level > 0) && (level < 1))
670 conf = level;
671 else
672 Warning("Divide","given confidence level %.3lf is invalid",level);
673 option.ReplaceAll("cl=","");
674 }
675
676 // look for statistic options
677 // check first Bayesian case
678 // bayesian with prior
679 Bool_t usePosteriorMode = false;
680 Bool_t useShortestInterval = false;
681 if (option.Contains("b(")) {
682 Double_t a = 0;
683 Double_t b = 0;
684 sscanf(strstr(option.Data(), "b("), "b(%lf,%lf)", &a, &b);
685 if (a > 0)
686 alpha = a;
687 else
688 Warning("Divide", "given shape parameter for alpha %.2lf is invalid", a);
689 if (b > 0)
690 beta = b;
691 else
692 Warning("Divide", "given shape parameter for beta %.2lf is invalid", b);
693 option.ReplaceAll("b(", "");
694 bIsBayesian = true;
695
696 // look for specific bayesian options
697
698 // use posterior mode
699
700 if (option.Contains("mode")) {
701 usePosteriorMode = true;
702 option.ReplaceAll("mode", "");
703 }
704 if (option.Contains("sh") || (usePosteriorMode && !option.Contains("cen"))) {
705 useShortestInterval = true;
706 }
707 }
708 // normal approximation
709 else if (option.Contains("n")) {
710 option.ReplaceAll("n", "");
711 pBound = &TEfficiency::Normal;
712 }
713 // clopper pearson interval
714 else if (option.Contains("cp")) {
715 option.ReplaceAll("cp", "");
717 }
718 // wilson interval
719 else if (option.Contains("w")) {
720 option.ReplaceAll("w", "");
721 pBound = &TEfficiency::Wilson;
722 }
723 // agresti coull interval
724 else if (option.Contains("ac")) {
725 option.ReplaceAll("ac", "");
727 }
728 // Feldman-Cousins interval
729 else if (option.Contains("fc")) {
730 option.ReplaceAll("fc", "");
732 }
733 // mid-P Lancaster interval
734 else if (option.Contains("midp")) {
735 option.ReplaceAll("midp", "");
737 }
738
739 // interpret as Poisson ratio
740 Bool_t bPoissonRatio = false;
741 if (option.Contains("pois")) {
742 bPoissonRatio = true;
743 option.ReplaceAll("pois", "");
744 }
745 Bool_t plot0Bins = false;
746 if (option.Contains("e0")) {
747 plot0Bins = true;
748 option.ReplaceAll("e0", "");
749 }
750
751 // weights works only in case of Normal approximation or Bayesian for binomial interval
752 // in case of Poisson ratio we can use weights by rescaling the obtained results using the effective entries
753 if ((bEffective && !bPoissonRatio) && !bIsBayesian && pBound != &TEfficiency::Normal) {
754 Warning("Divide", "Histograms have weights: only Normal or Bayesian error calculation is supported");
755 Info("Divide", "Using now the Normal approximation for weighted histograms");
756 }
757
758 if (bPoissonRatio) {
759 if (pass->GetDimension() != total->GetDimension()) {
760 Error("Divide", "passed histograms are not of the same dimension");
761 return;
762 }
763
764 if (!TEfficiency::CheckBinning(*pass, *total)) {
765 Error("Divide", "passed histograms are not consistent");
766 return;
767 }
768 } else {
769 // check consistency of histograms, allowing weights
770 if (!TEfficiency::CheckConsistency(*pass, *total, "w")) {
771 Error("Divide", "passed histograms are not consistent");
772 return;
773 }
774 }
775
776 // Set the graph to have a number of points equal to the number of histogram
777 // bins
778 Int_t nbins = pass->GetNbinsX();
779 Set(nbins);
780
781 // Ok, now set the points for each bin
782 // (Note: the TH1 bin content is shifted to the right by one:
783 // bin=0 is underflow, bin=nbins+1 is overflow.)
784
785 //efficiency with lower and upper boundary of confidence interval
786 double eff, low, upper;
787 //this keeps track of the number of points added to the graph
788 int npoint=0;
789 //number of total and passed events
790 Double_t t = 0 , p = 0;
791 Double_t tw = 0, tw2 = 0, pw = 0, pw2 = 0, wratio = 1; // for the case of weights
792 //loop over all bins and fill the graph
793 for (Int_t b=1; b<=nbins; ++b) {
794
795 // default value when total =0;
796 eff = 0;
797 low = 0;
798 upper = 0;
799
800 // special case in case of weights we have to consider the sum of weights and the sum of weight squares
801 if (bEffective) {
802 tw = total->GetBinContent(b);
803 tw2 = (total->GetSumw2()->fN > 0) ? total->GetSumw2()->At(b) : tw;
804 pw = pass->GetBinContent(b);
805 pw2 = (pass->GetSumw2()->fN > 0) ? pass->GetSumw2()->At(b) : pw;
806
807 if (bPoissonRatio) {
808 // tw += pw;
809 // tw2 += pw2;
810 // compute ratio on the effective entries ( p and t)
811 // special case is when (pw=0, pw2=0) in this case we cannot get the bin weight.
812 // we use then the overall weight of the full histogram
813 if (pw == 0 && pw2 == 0)
814 p = 0;
815 else
816 p = (pw * pw) / pw2;
817
818 if (tw == 0 && tw2 == 0)
819 t = 0;
820 else
821 t = (tw * tw) / tw2;
822
823 if (pw > 0 && tw > 0)
824 // this is the ratio of the two bin weights ( pw/p / t/tw )
825 wratio = (pw * t) / (p * tw);
826 else if (pw == 0 && tw > 0)
827 // case p histogram has zero compute the weights from all the histogram
828 // weight of histogram - sumw2/sumw
829 wratio = (psumw2 * t) / (psumw * tw);
830 else if (tw == 0 && pw > 0)
831 // case t histogram has zero compute the weights from all the histogram
832 // weight of histogram - sumw2/sumw
833 wratio = (pw * tsumw) / (p * tsumw2);
834 else if (p > 0)
835 wratio = pw / p; // not sure if needed
836 else {
837 // case both pw and tw are zero - we skip these bins
838 if (!plot0Bins) continue; // skip bins with total <= 0
839 }
840
841 t += p;
842 // std::cout << p << " " << t << " " << wratio << std::endl;
843 } else if (tw <= 0 && !plot0Bins)
844 continue; // skip bins with total <= 0
845
846 // in the case of weights have the formula only for
847 // the normal and bayesian statistics (see below)
848
849 }
850
851 // use bin contents
852 else {
853 t = std::round(total->GetBinContent(b));
854 p = std::round(pass->GetBinContent(b));
855
856 if (bPoissonRatio)
857 t += p;
858
859 if (t == 0.0 && !plot0Bins)
860 continue; // skip bins with total = 0
861 }
862
863 //using bayesian statistics
864 if(bIsBayesian) {
865 double aa,bb;
866
867 if ((bEffective && !bPoissonRatio) && tw2 <= 0) {
868 // case of bins with zero errors
869 eff = pw/tw;
870 low = eff; upper = eff;
871 }
872 else {
873
874 if (bEffective && !bPoissonRatio) {
875 // tw/tw2 re-normalize the weights
876 double norm = tw/tw2; // case of tw2 = 0 is treated above
877 aa = pw * norm + alpha;
878 bb = (tw - pw) * norm + beta;
879 }
880 else {
881 aa = double(p) + alpha;
882 bb = double(t-p) + beta;
883 }
884 if (usePosteriorMode)
885 eff = TEfficiency::BetaMode(aa,bb);
886 else
887 eff = TEfficiency::BetaMean(aa,bb);
888
889 if (useShortestInterval) {
890 TEfficiency::BetaShortestInterval(conf,aa,bb,low,upper);
891 }
892 else {
893 low = TEfficiency::BetaCentralInterval(conf,aa,bb,false);
894 upper = TEfficiency::BetaCentralInterval(conf,aa,bb,true);
895 }
896 }
897 }
898 // case of non-bayesian statistics
899 else {
900 if (bEffective && !bPoissonRatio) {
901
902 if (tw > 0) {
903
904 eff = pw/tw;
905
906 // use normal error calculation using variance of MLE with weights (F.James 8.5.2)
907 // this is the same formula used in ROOT for TH1::Divide("B")
908
909 double variance = ( pw2 * (1. - 2 * eff) + tw2 * eff *eff ) / ( tw * tw) ;
910 double sigma = sqrt(variance);
911
912 double prob = 0.5 * (1.-conf);
913 double delta = ROOT::Math::normal_quantile_c(prob, sigma);
914 low = eff - delta;
915 upper = eff + delta;
916 if (low < 0) low = 0;
917 if (upper > 1) upper = 1.;
918 }
919 }
920 else {
921 // when not using weights (all cases) or in case of Poisson ratio with weights
922 if(t != 0.0)
923 eff = ((Double_t)p)/t;
924
925 low = pBound(t,p,conf,false);
926 upper = pBound(t,p,conf,true);
927 }
928 }
929 // treat as Poisson ratio
930 if(bPoissonRatio)
931 {
932 Double_t ratio = eff/(1 - eff);
933 // take the intervals in eff as intervals in the Poisson ratio
934 low = low/(1. - low);
935 upper = upper/(1.-upper);
936 eff = ratio;
937 if (bEffective) {
938 //scale result by the ratio of the weight
939 eff *= wratio;
940 low *= wratio;
941 upper *= wratio;
942 }
943 }
944 //Set the point center and its errors
945 if (TMath::Finite(eff)) {
946 SetPoint(npoint,pass->GetBinCenter(b),eff);
947 SetPointError(npoint,
948 pass->GetBinCenter(b)-pass->GetBinLowEdge(b),
949 pass->GetBinLowEdge(b)-pass->GetBinCenter(b)+pass->GetBinWidth(b),
950 eff-low,upper-eff);
951 npoint++;//we have added a point to the graph
952 }
953 }
954
955 Set(npoint);//tell the graph how many points we've really added
956 if (npoint < nbins)
957 Warning("Divide","Number of graph points is different than histogram bins - %d points have been skipped",nbins-npoint);
958
959
960 if (bVerbose) {
961 Info("Divide","made a graph with %d points from %d bins",npoint,nbins);
962 Info("Divide","used confidence level: %.2lf\n",conf);
963 if(bIsBayesian)
964 Info("Divide","used prior probability ~ beta(%.2lf,%.2lf)",alpha,beta);
965 Print();
966 }
967}
968
969////////////////////////////////////////////////////////////////////////////////
970/// Compute Range.
971
973{
975
976 for (Int_t i=0;i<fNpoints;i++) {
977 if (fX[i] -fEXlow[i] < xmin) {
978 if (gPad && gPad->GetLogx()) {
979 if (fEXlow[i] < fX[i]) xmin = fX[i]-fEXlow[i];
980 else xmin = TMath::Min(xmin,fX[i]/3);
981 } else {
982 xmin = fX[i]-fEXlow[i];
983 }
984 }
985 if (fX[i] +fEXhigh[i] > xmax) xmax = fX[i]+fEXhigh[i];
986 if (fY[i] -fEYlow[i] < ymin) {
987 if (gPad && gPad->GetLogy()) {
988 if (fEYlow[i] < fY[i]) ymin = fY[i]-fEYlow[i];
989 else ymin = TMath::Min(ymin,fY[i]/3);
990 } else {
991 ymin = fY[i]-fEYlow[i];
992 }
993 }
994 if (fY[i] +fEYhigh[i] > ymax) ymax = fY[i]+fEYhigh[i];
995 }
996}
997
998
999////////////////////////////////////////////////////////////////////////////////
1000/// Copy and release.
1001
1003 Int_t ibegin, Int_t iend, Int_t obegin)
1004{
1005 CopyPoints(newarrays, ibegin, iend, obegin);
1006 if (newarrays) {
1007 delete[] fEXlow;
1008 fEXlow = newarrays[0];
1009 delete[] fEXhigh;
1010 fEXhigh = newarrays[1];
1011 delete[] fEYlow;
1012 fEYlow = newarrays[2];
1013 delete[] fEYhigh;
1014 fEYhigh = newarrays[3];
1015 delete[] fX;
1016 fX = newarrays[4];
1017 delete[] fY;
1018 fY = newarrays[5];
1019 delete[] newarrays;
1020 }
1021}
1022
1023
1024////////////////////////////////////////////////////////////////////////////////
1025/// Copy errors from `fE***` to `arrays[***]`
1026/// or to `f***` Copy points.
1027
1029 Int_t ibegin, Int_t iend, Int_t obegin)
1030{
1031 if (TGraph::CopyPoints(arrays ? arrays+4 : nullptr, ibegin, iend, obegin)) {
1032 Int_t n = (iend - ibegin)*sizeof(Double_t);
1033 if (arrays) {
1034 memmove(&arrays[0][obegin], &fEXlow[ibegin], n);
1035 memmove(&arrays[1][obegin], &fEXhigh[ibegin], n);
1036 memmove(&arrays[2][obegin], &fEYlow[ibegin], n);
1037 memmove(&arrays[3][obegin], &fEYhigh[ibegin], n);
1038 } else {
1039 memmove(&fEXlow[obegin], &fEXlow[ibegin], n);
1040 memmove(&fEXhigh[obegin], &fEXhigh[ibegin], n);
1041 memmove(&fEYlow[obegin], &fEYlow[ibegin], n);
1042 memmove(&fEYhigh[obegin], &fEYhigh[ibegin], n);
1043 }
1044 return kTRUE;
1045 } else {
1046 return kFALSE;
1047 }
1048}
1049
1050
1051////////////////////////////////////////////////////////////////////////////////
1052/// Should be called from ctors after `fNpoints` has been set.
1053/// Note: This function should be called only from the constructor
1054/// since it does not delete previously existing arrays
1055
1057{
1058 if (!fNpoints) {
1059 fEXlow = fEYlow = fEXhigh = fEYhigh = nullptr;
1060 return kFALSE;
1061 }
1062 fEXlow = new Double_t[fMaxSize];
1063 fEYlow = new Double_t[fMaxSize];
1064 fEXhigh = new Double_t[fMaxSize];
1065 fEYhigh = new Double_t[fMaxSize];
1066 return kTRUE;
1067}
1068
1069////////////////////////////////////////////////////////////////////////////////
1070/// Protected function to perform the merge operation of a graph with asymmetric errors.
1071
1073{
1074 if (g->GetN() == 0) return kFALSE;
1075
1076 Double_t * exl = g->GetEXlow();
1077 Double_t * exh = g->GetEXhigh();
1078 Double_t * eyl = g->GetEYlow();
1079 Double_t * eyh = g->GetEYhigh();
1080 if (exl == nullptr || exh == nullptr || eyl == nullptr || eyh == nullptr) {
1081 if (g->IsA() != TGraph::Class() )
1082 Warning("DoMerge","Merging a %s is not compatible with a TGraphAsymmErrors - errors will be ignored",g->IsA()->GetName());
1083 return TGraph::DoMerge(g);
1084 }
1085 for (Int_t i = 0 ; i < g->GetN(); i++) {
1086 Int_t ipoint = GetN();
1087 Double_t x = g->GetX()[i];
1088 Double_t y = g->GetY()[i];
1089 SetPoint(ipoint, x, y);
1090 SetPointError(ipoint, exl[i], exh[i], eyl[i], eyh[i] );
1091 }
1092
1093 return kTRUE;
1094}
1095
1096////////////////////////////////////////////////////////////////////////////////
1097/// Set zero values for point arrays in the range `[begin, end]`
1098
1100 Bool_t from_ctor)
1101{
1102 if (!from_ctor) {
1103 TGraph::FillZero(begin, end, from_ctor);
1104 }
1105 Int_t n = (end - begin)*sizeof(Double_t);
1106 memset(fEXlow + begin, 0, n);
1107 memset(fEXhigh + begin, 0, n);
1108 memset(fEYlow + begin, 0, n);
1109 memset(fEYhigh + begin, 0, n);
1110}
1111
1112
1113////////////////////////////////////////////////////////////////////////////////
1114/// Returns the combined error along X at point i by computing the average
1115/// of the lower and upper variance.
1116
1118{
1119 if (i < 0 || i >= fNpoints) return -1;
1120 if (!fEXlow && !fEXhigh) return -1;
1121 Double_t elow=0, ehigh=0;
1122 if (fEXlow) elow = fEXlow[i];
1123 if (fEXhigh) ehigh = fEXhigh[i];
1124 return TMath::Sqrt(0.5*(elow*elow + ehigh*ehigh));
1125}
1126
1127
1128////////////////////////////////////////////////////////////////////////////////
1129/// Returns the combined error along Y at point i by computing the average
1130/// of the lower and upper variance.
1131
1133{
1134 if (i < 0 || i >= fNpoints) return -1;
1135 if (!fEYlow && !fEYhigh) return -1;
1136 Double_t elow=0, ehigh=0;
1137 if (fEYlow) elow = fEYlow[i];
1138 if (fEYhigh) ehigh = fEYhigh[i];
1139 return TMath::Sqrt(0.5*(elow*elow + ehigh*ehigh));
1140}
1141
1142
1143////////////////////////////////////////////////////////////////////////////////
1144/// Get high error on X.
1145
1147{
1148 if (i<0 || i>fNpoints) return -1;
1149 if (fEXhigh) return fEXhigh[i];
1150 return -1;
1151}
1152
1153
1154////////////////////////////////////////////////////////////////////////////////
1155/// Get low error on X.
1156
1158{
1159 if (i<0 || i>fNpoints) return -1;
1160 if (fEXlow) return fEXlow[i];
1161 return -1;
1162}
1163
1164
1165////////////////////////////////////////////////////////////////////////////////
1166/// Get high error on Y.
1167
1169{
1170 if (i<0 || i>fNpoints) return -1;
1171 if (fEYhigh) return fEYhigh[i];
1172 return -1;
1173}
1174
1175
1176////////////////////////////////////////////////////////////////////////////////
1177/// Get low error on Y.
1178
1180{
1181 if (i<0 || i>fNpoints) return -1;
1182 if (fEYlow) return fEYlow[i];
1183 return -1;
1184}
1185
1186
1187////////////////////////////////////////////////////////////////////////////////
1188/// Adds all graphs with asymmetric errors from the collection to this graph.
1189/// Returns the total number of points in the result or -1 in case of an error.
1190
1192{
1193 TIter next(li);
1194 while (TObject* o = next()) {
1195 TGraph *g = dynamic_cast<TGraph*>(o);
1196 if (!g) {
1197 Error("Merge",
1198 "Cannot merge - an object which doesn't inherit from TGraph found in the list");
1199 return -1;
1200 }
1201 int n0 = GetN();
1202 int n1 = n0+g->GetN();
1203 Set(n1);
1204 Double_t * x = g->GetX();
1205 Double_t * y = g->GetY();
1206 Double_t * exlow = g->GetEXlow();
1207 Double_t * exhigh = g->GetEXhigh();
1208 Double_t * eylow = g->GetEYlow();
1209 Double_t * eyhigh = g->GetEYhigh();
1210 for (Int_t i = 0 ; i < g->GetN(); i++) {
1211 SetPoint(n0+i, x[i], y[i]);
1212 if (exlow) fEXlow[n0+i] = exlow[i];
1213 if (exhigh) fEXhigh[n0+i] = exhigh[i];
1214 if (eylow) fEYlow[n0+i] = eylow[i];
1215 if (eyhigh) fEYhigh[n0+i] = eyhigh[i];
1216 }
1217 }
1218 return GetN();
1219}
1220
1221////////////////////////////////////////////////////////////////////////////////
1222/// Print graph and errors values.
1223
1225{
1226 for (Int_t i=0;i<fNpoints;i++) {
1227 printf("x[%d]=%g, y[%d]=%g, exl[%d]=%g, exh[%d]=%g, eyl[%d]=%g, eyh[%d]=%g\n"
1228 ,i,fX[i],i,fY[i],i,fEXlow[i],i,fEXhigh[i],i,fEYlow[i],i,fEYhigh[i]);
1229 }
1230}
1231
1232
1233////////////////////////////////////////////////////////////////////////////////
1234/// Save primitive as a C++ statement(s) on output stream out.
1235
1236void TGraphAsymmErrors::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1237{
1238 out << " " << std::endl;
1239 static Int_t frameNumber = 3000;
1240 frameNumber++;
1241
1242 auto fXName = SaveArray(out, "fx", frameNumber, fX);
1243 auto fYName = SaveArray(out, "fy", frameNumber, fY);
1244 auto fElXName = SaveArray(out, "felx", frameNumber, fEXlow);
1245 auto fElYName = SaveArray(out, "fely", frameNumber, fEYlow);
1246 auto fEhXName = SaveArray(out, "fehx", frameNumber, fEXhigh);
1247 auto fEhYName = SaveArray(out, "fehy", frameNumber, fEYhigh);
1248
1249 if (gROOT->ClassSaved(TGraphAsymmErrors::Class()))
1250 out<<" ";
1251 else
1252 out << " TGraphAsymmErrors *";
1253 out << "grae = new TGraphAsymmErrors("<< fNpoints << ","
1254 << fXName << "," << fYName << ","
1255 << fElXName << "," << fEhXName << ","
1256 << fElYName << "," << fEhYName << ");"
1257 << std::endl;
1258
1259 SaveHistogramAndFunctions(out, "grae", frameNumber, option);
1260}
1261
1262////////////////////////////////////////////////////////////////////////////////
1263/// Multiply the values and errors of a TGraphAsymmErrors by a constant c1.
1264///
1265/// If option contains "x" the x values and errors are scaled
1266/// If option contains "y" the y values and errors are scaled
1267/// If option contains "xy" both x and y values and errors are scaled
1268
1270{
1272 TString opt = option; opt.ToLower();
1273 if (opt.Contains("x") && GetEXlow()) {
1274 for (Int_t i=0; i<GetN(); i++)
1275 GetEXlow()[i] *= c1;
1276 }
1277 if (opt.Contains("x") && GetEXhigh()) {
1278 for (Int_t i=0; i<GetN(); i++)
1279 GetEXhigh()[i] *= c1;
1280 }
1281 if (opt.Contains("y") && GetEYlow()) {
1282 for (Int_t i=0; i<GetN(); i++)
1283 GetEYlow()[i] *= c1;
1284 }
1285 if (opt.Contains("y") && GetEYhigh()) {
1286 for (Int_t i=0; i<GetN(); i++)
1287 GetEYhigh()[i] *= c1;
1288 }
1289}
1290
1291////////////////////////////////////////////////////////////////////////////////
1292/// Set ex and ey values for point pointed by the mouse.
1293
1295{
1296 if (!gPad) {
1297 Error("SetPointError", "Cannot be used without gPad, requires last mouse position");
1298 return;
1299 }
1300
1301 Int_t px = gPad->GetEventX();
1302 Int_t py = gPad->GetEventY();
1303
1304 //localize point to be deleted
1305 Int_t ipoint = -2;
1306 Int_t i;
1307 // start with a small window (in case the mouse is very close to one point)
1308 for (i=0;i<fNpoints;i++) {
1309 Int_t dpx = px - gPad->XtoAbsPixel(gPad->XtoPad(fX[i]));
1310 Int_t dpy = py - gPad->YtoAbsPixel(gPad->YtoPad(fY[i]));
1311 if (dpx*dpx+dpy*dpy < 25) {ipoint = i; break;}
1312 }
1313 if (ipoint == -2) return;
1314
1315 fEXlow[ipoint] = exl;
1316 fEYlow[ipoint] = eyl;
1317 fEXhigh[ipoint] = exh;
1318 fEYhigh[ipoint] = eyh;
1319 gPad->Modified();
1320}
1321
1322
1323////////////////////////////////////////////////////////////////////////////////
1324/// Set ex and ey values for point number i.
1325
1327{
1328 if (i < 0) return;
1329 if (i >= fNpoints) {
1330 // re-allocate the object
1332 }
1333 fEXlow[i] = exl;
1334 fEYlow[i] = eyl;
1335 fEXhigh[i] = exh;
1336 fEYhigh[i] = eyh;
1337}
1338
1339
1340////////////////////////////////////////////////////////////////////////////////
1341/// Set EXlow for point `i`.
1342
1344{
1345 if (i < 0) return;
1346 if (i >= fNpoints) {
1347 // re-allocate the object
1349 }
1350 fEXlow[i] = exl;
1351}
1352
1353
1354////////////////////////////////////////////////////////////////////////////////
1355/// Set EXhigh for point `i`.
1356
1358{
1359 if (i < 0) return;
1360 if (i >= fNpoints) {
1361 // re-allocate the object
1363 }
1364 fEXhigh[i] = exh;
1365}
1366
1367
1368////////////////////////////////////////////////////////////////////////////////
1369/// Set EYlow for point `i`.
1370
1372{
1373 if (i < 0) return;
1374 if (i >= fNpoints) {
1375 // re-allocate the object
1377 }
1378 fEYlow[i] = eyl;
1379}
1380
1381
1382////////////////////////////////////////////////////////////////////////////////
1383/// Set EYhigh for point `i`.
1384
1386{
1387 if (i < 0) return;
1388 if (i >= fNpoints) {
1389 // re-allocate the object
1391 }
1392 fEYhigh[i] = eyh;
1393}
1394
1395
1396////////////////////////////////////////////////////////////////////////////////
1397/// Stream an object of class TGraphAsymmErrors.
1398
1400{
1401 if (b.IsReading()) {
1402 UInt_t R__s, R__c;
1403 Version_t R__v = b.ReadVersion(&R__s, &R__c);
1404 if (R__v > 2) {
1405 b.ReadClassBuffer(TGraphAsymmErrors::Class(), this, R__v, R__s, R__c);
1406 return;
1407 }
1408 //====process old versions before automatic schema evolution
1410 fEXlow = new Double_t[fNpoints];
1411 fEYlow = new Double_t[fNpoints];
1412 fEXhigh = new Double_t[fNpoints];
1413 fEYhigh = new Double_t[fNpoints];
1414 if (R__v < 2) {
1415 Float_t *exlow = new Float_t[fNpoints];
1416 Float_t *eylow = new Float_t[fNpoints];
1417 Float_t *exhigh = new Float_t[fNpoints];
1418 Float_t *eyhigh = new Float_t[fNpoints];
1419 b.ReadFastArray(exlow,fNpoints);
1420 b.ReadFastArray(eylow,fNpoints);
1421 b.ReadFastArray(exhigh,fNpoints);
1422 b.ReadFastArray(eyhigh,fNpoints);
1423 for (Int_t i=0;i<fNpoints;i++) {
1424 fEXlow[i] = exlow[i];
1425 fEYlow[i] = eylow[i];
1426 fEXhigh[i] = exhigh[i];
1427 fEYhigh[i] = eyhigh[i];
1428 }
1429 delete [] eylow;
1430 delete [] exlow;
1431 delete [] eyhigh;
1432 delete [] exhigh;
1433 } else {
1434 b.ReadFastArray(fEXlow,fNpoints);
1435 b.ReadFastArray(fEYlow,fNpoints);
1436 b.ReadFastArray(fEXhigh,fNpoints);
1437 b.ReadFastArray(fEYhigh,fNpoints);
1438 }
1439 b.CheckByteCount(R__s, R__c, TGraphAsymmErrors::IsA());
1440 //====end of old versions
1441
1442 } else {
1443 b.WriteClassBuffer(TGraphAsymmErrors::Class(),this);
1444 }
1445}
1446
1447
1448////////////////////////////////////////////////////////////////////////////////
1449/// Swap points.
1450
1452{
1453 SwapValues(fEXlow, pos1, pos2);
1454 SwapValues(fEXhigh, pos1, pos2);
1455 SwapValues(fEYlow, pos1, pos2);
1456 SwapValues(fEYhigh, pos1, pos2);
1457 TGraph::SwapPoints(pos1, pos2);
1458}
1459
1460////////////////////////////////////////////////////////////////////////////////
1461/// Update the fX, fY, fEXlow, fEXhigh, fEYlow and fEYhigh arrays with the sorted values.
1462
1463void TGraphAsymmErrors::UpdateArrays(const std::vector<Int_t> &sorting_indices, Int_t numSortedPoints, Int_t low)
1464{
1465 std::vector<Double_t> fEXlowSorted(numSortedPoints);
1466 std::vector<Double_t> fEXhighSorted(numSortedPoints);
1467 std::vector<Double_t> fEYlowSorted(numSortedPoints);
1468 std::vector<Double_t> fEYhighSorted(numSortedPoints);
1469
1470 // Fill the sorted X and Y error values based on the sorted indices
1471 std::generate(fEXlowSorted.begin(), fEXlowSorted.end(),
1472 [begin = low, &sorting_indices, this]() mutable { return fEXlow[sorting_indices[begin++]]; });
1473 std::generate(fEXhighSorted.begin(), fEXhighSorted.end(),
1474 [begin = low, &sorting_indices, this]() mutable { return fEXhigh[sorting_indices[begin++]]; });
1475 std::generate(fEYlowSorted.begin(), fEYlowSorted.end(),
1476 [begin = low, &sorting_indices, this]() mutable { return fEYlow[sorting_indices[begin++]]; });
1477 std::generate(fEYhighSorted.begin(), fEYhighSorted.end(),
1478 [begin = low, &sorting_indices, this]() mutable { return fEYhigh[sorting_indices[begin++]]; });
1479
1480 // Copy the sorted X and Y error values back to the original arrays
1481 std::copy(fEXlowSorted.begin(), fEXlowSorted.end(), fEXlow + low);
1482 std::copy(fEXhighSorted.begin(), fEXhighSorted.end(), fEXhigh + low);
1483 std::copy(fEYlowSorted.begin(), fEYlowSorted.end(), fEYlow + low);
1484 std::copy(fEYhighSorted.begin(), fEYhighSorted.end(), fEYhigh + low);
1485
1486 TGraph::UpdateArrays(sorting_indices, numSortedPoints, low);
1487}
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
bool Bool_t
Definition RtypesCore.h:63
short Version_t
Definition RtypesCore.h:65
float Float_t
Definition RtypesCore.h:57
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
static unsigned int total
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t format
float xmin
float ymin
float xmax
float ymax
#define gROOT
Definition TROOT.h:406
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
R__EXTERN TSystem * gSystem
Definition TSystem.h:555
#define gPad
Double_t At(Int_t i) const
Definition TArrayD.h:79
Int_t fN
Definition TArray.h:38
Buffer base class used for serializing objects.
Definition TBuffer.h:43
Collection abstract base class.
Definition TCollection.h:65
static Double_t BetaMode(Double_t alpha, Double_t beta)
Compute the mode of the beta distribution.
static Bool_t BetaShortestInterval(Double_t level, Double_t alpha, Double_t beta, Double_t &lower, Double_t &upper)
Calculates the boundaries for a shortest confidence interval for a Beta distribution.
static Double_t BetaMean(Double_t alpha, Double_t beta)
Compute the mean (average) of the beta distribution.
static Double_t AgrestiCoull(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries for the frequentist Agresti-Coull interval.
static Double_t FeldmanCousins(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries for the frequentist Feldman-Cousins interval.
static Bool_t CheckBinning(const TH1 &pass, const TH1 &total)
Checks binning for each axis.
static Double_t BetaCentralInterval(Double_t level, Double_t alpha, Double_t beta, Bool_t bUpper)
Calculates the boundaries for a central confidence interval for a Beta distribution.
static Double_t MidPInterval(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries using the mid-P binomial interval (Lancaster method) from B.
static Double_t Normal(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Returns the confidence limits for the efficiency supposing that the efficiency follows a normal distr...
static Double_t Wilson(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries for the frequentist Wilson interval.
static Bool_t CheckConsistency(const TH1 &pass, const TH1 &total, Option_t *opt="")
Checks the consistence of the given histograms.
static Double_t ClopperPearson(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries for the frequentist Clopper-Pearson interval.
1-Dim function class
Definition TF1.h:233
TGraph with asymmetric error bars.
Double_t * GetEXlow() const override
virtual void SetPointEYlow(Int_t i, Double_t eyl)
Set EYlow for point i.
Double_t * fEXhigh
[fNpoints] array of X high errors
Double_t GetErrorY(Int_t bin) const override
Returns the combined error along Y at point i by computing the average of the lower and upper varianc...
void FillZero(Int_t begin, Int_t end, Bool_t from_ctor=kTRUE) override
Set zero values for point arrays in the range [begin, end]
virtual void Divide(const TH1 *pass, const TH1 *total, Option_t *opt="cp")
Fill this TGraphAsymmErrors by dividing two 1-dimensional histograms pass/total.
Bool_t CtorAllocate()
Should be called from ctors after fNpoints has been set.
static TClass * Class()
virtual void SetPointError(Double_t exl, Double_t exh, Double_t eyl, Double_t eyh)
Set ex and ey values for point pointed by the mouse.
Double_t GetErrorXhigh(Int_t i) const override
Get high error on X.
Double_t * fEYhigh
[fNpoints] array of Y high errors
Bool_t CopyPoints(Double_t **arrays, Int_t ibegin, Int_t iend, Int_t obegin) override
Copy errors from fE*** to arrays[***] or to f*** Copy points.
void ComputeRange(Double_t &xmin, Double_t &ymin, Double_t &xmax, Double_t &ymax) const override
Compute Range.
virtual void SetPointEXlow(Int_t i, Double_t exl)
Set EXlow for point i.
void Streamer(TBuffer &) override
Stream an object of class TGraphAsymmErrors.
virtual void SetPointEYhigh(Int_t i, Double_t eyh)
Set EYhigh for point i.
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save primitive as a C++ statement(s) on output stream out.
virtual void BayesDivide(const TH1 *pass, const TH1 *total, Option_t *opt="")
This function is only kept for backward compatibility.
Double_t * GetEYhigh() const override
void Print(Option_t *chopt="") const override
Print graph and errors values.
Double_t GetErrorYhigh(Int_t i) const override
Get high error on Y.
Int_t Merge(TCollection *list) override
Adds all graphs with asymmetric errors from the collection to this graph.
Double_t GetErrorXlow(Int_t i) const override
Get low error on X.
Double_t * fEYlow
[fNpoints] array of Y low errors
void UpdateArrays(const std::vector< Int_t > &sorting_indices, Int_t numSortedPoints, Int_t low) override
Update the fX, fY, fEXlow, fEXhigh, fEYlow and fEYhigh arrays with the sorted values.
void SwapPoints(Int_t pos1, Int_t pos2) override
Swap points.
Double_t * GetEXhigh() const override
Double_t * fEXlow
[fNpoints] array of X low errors
void Apply(TF1 *f) override
Apply a function to all data points .
Double_t ** Allocate(Int_t size) override
Allocate internal data structures for size points.
Double_t * GetEYlow() const override
Bool_t DoMerge(const TGraph *g) override
Protected function to perform the merge operation of a graph with asymmetric errors.
void CopyAndRelease(Double_t **newarrays, Int_t ibegin, Int_t iend, Int_t obegin) override
Copy and release.
TClass * IsA() const override
virtual void SetPointEXhigh(Int_t i, Double_t exh)
Set EXhigh for point i.
Double_t GetErrorYlow(Int_t i) const override
Get low error on Y.
Double_t GetErrorX(Int_t bin) const override
Returns the combined error along X at point i by computing the average of the lower and upper varianc...
TGraphAsymmErrors()
TGraphAsymmErrors default constructor.
TGraphAsymmErrors & operator=(const TGraphAsymmErrors &gr)
TGraphAsymmErrors assignment operator.
~TGraphAsymmErrors() override
TGraphAsymmErrors default destructor.
void Scale(Double_t c1=1., Option_t *option="y") override
Multiply the values and errors of a TGraphAsymmErrors by a constant c1.
static Int_t CalculateScanfFields(const char *fmt)
Calculate scan fields.
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
static TClass * Class()
Int_t fNpoints
Number of points <= fMaxSize.
Definition TGraph.h:46
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2325
Int_t fMaxSize
!Current dimension of arrays fX and fY
Definition TGraph.h:45
void SaveHistogramAndFunctions(std::ostream &out, const char *varname, Int_t &frameNumber, Option_t *option)
Save histogram and list of functions of TGraph as C++ statement Used in all TGraph-derived classes.
Definition TGraph.cxx:2182
TH1F * fHistogram
Pointer to histogram used for drawing axis.
Definition TGraph.h:50
virtual void UpdateArrays(const std::vector< Int_t > &sorting_indices, Int_t numSortedPoints, Int_t low)
Update the fX and fY arrays with the sorted values.
Definition TGraph.cxx:2575
Int_t GetN() const
Definition TGraph.h:131
Double_t * fY
[fNpoints] array of Y points
Definition TGraph.h:48
Bool_t CtorAllocate()
In constructors set fNpoints than call this method.
Definition TGraph.cxx:788
TString SaveArray(std::ostream &out, const char *suffix, Int_t frameNumber, Double_t *arr)
Save array as C++ code Returns name of created array.
Definition TGraph.cxx:2158
virtual void ComputeRange(Double_t &xmin, Double_t &ymin, Double_t &xmax, Double_t &ymax) const
Compute the x/y range of the points in this graph.
Definition TGraph.cxx:714
Double_t ** AllocateArrays(Int_t Narrays, Int_t arraySize)
Allocate arrays.
Definition TGraph.cxx:598
virtual void Scale(Double_t c1=1., Option_t *option="y")
Multiply the values of a TGraph by a constant c1.
Definition TGraph.cxx:2242
static void SwapValues(Double_t *arr, Int_t pos1, Int_t pos2)
Swap values.
Definition TGraph.cxx:2594
void Streamer(TBuffer &) override
Stream an object of class TGraph.
Definition TGraph.cxx:2499
virtual Bool_t DoMerge(const TGraph *g)
protected function to perform the merge operation of a graph
Definition TGraph.cxx:2659
void SetName(const char *name="") override
Set graph name.
Definition TGraph.cxx:2364
virtual void SwapPoints(Int_t pos1, Int_t pos2)
Swap points.
Definition TGraph.cxx:2566
virtual void FillZero(Int_t begin, Int_t end, Bool_t from_ctor=kTRUE)
Set zero values for point arrays in the range [begin, end) Should be redefined in descendant classes.
Definition TGraph.cxx:1085
Double_t * fX
[fNpoints] array of X points
Definition TGraph.h:47
void SetTitle(const char *title="") override
Change (i.e.
Definition TGraph.cxx:2380
virtual void Set(Int_t n)
Set number of points in the graph Existing coordinates are preserved New coordinates above fNpoints a...
Definition TGraph.cxx:2260
virtual Int_t GetPoint(Int_t i, Double_t &x, Double_t &y) const
Get x and y values for point number i.
Definition TGraph.cxx:1516
virtual Bool_t CopyPoints(Double_t **newarrays, Int_t ibegin, Int_t iend, Int_t obegin)
Copy points from fX and fY to arrays[0] and arrays[1] or to fX and fY if arrays == 0 and ibegin !...
Definition TGraph.cxx:762
TGraph & operator=(const TGraph &)
Equal operator for this graph.
Definition TGraph.cxx:232
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition TH1.cxx:9109
virtual Int_t GetDimension() const
Definition TH1.h:283
virtual Int_t GetNcells() const
Definition TH1.h:300
virtual Int_t GetNbinsX() const
Definition TH1.h:297
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition TH1.cxx:9120
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5029
virtual TArrayD * GetSumw2()
Definition TH1.h:313
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition TH1.cxx:9131
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition TH1.cxx:7885
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition TNamed.h:48
Mother of all ROOT objects.
Definition TObject.h:41
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:973
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:987
void MakeZombie()
Definition TObject.h:53
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:961
Basic string class.
Definition TString.h:139
Ssiz_t Length() const
Definition TString.h:417
void ToLower()
Change string to lower-case.
Definition TString.cxx:1182
Int_t Atoi() const
Return integer value of string.
Definition TString.cxx:1988
Double_t Atof() const
Return floating-point value contained in string.
Definition TString.cxx:2054
Bool_t IsFloat() const
Returns kTRUE if string contains a floating point or integer number.
Definition TString.cxx:1858
const char * Data() const
Definition TString.h:376
Bool_t IsDigit() const
Returns true if all characters in string are digits (0-9) or white spaces, i.e.
Definition TString.cxx:1830
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition TString.h:704
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2378
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
Float_t GetErrorX() const
Definition TStyle.h:185
virtual Bool_t ExpandPathName(TString &path)
Expand a pathname getting rid of special shell characters like ~.
Definition TSystem.cxx:1274
TVectorT.
Definition TVectorT.h:27
Int_t GetNrows() const
Definition TVectorT.h:73
Int_t GetLwb() const
Definition TVectorT.h:71
TLine * line
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
const Double_t sigma
Double_t y[n]
Definition legend1.C:17
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TGraphErrors * gr
Definition legend1.C:25
Int_t Finite(Double_t x)
Check if it is finite with a mask in order to be consistent in presence of fast math.
Definition TMath.h:770
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:198
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123