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