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