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