Logo ROOT   6.19/01
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"
30 #include "Math/QuantFuncMathCore.h"
31 
33 
34 /** \class TGraphAsymmErrors
35  \ingroup Hist
36 TGraph with asymmetric error bars.
37 
38 The TGraphAsymmErrors painting is performed thanks to the TGraphPainter
39 class. All details about the various painting options are given in this class.
40 <p>
41 The picture below gives an example:
42 
43 Begin_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 }
64 End_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 
137 TGraphAsymmErrors::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 
160 TGraphAsymmErrors::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 
184 TGraphAsymmErrors::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 
214 TGraphAsymmErrors::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 
298 TGraphAsymmErrors::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 
574 void 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.)
638  Double_t (*pBound)(Double_t,Double_t,Double_t,Bool_t) = &TEfficiency::ClopperPearson; // default method
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","");
676  pBound = &TEfficiency::ClopperPearson;
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","");
688  pBound = &TEfficiency::AgrestiCoull;
689  }
690  // Feldman-Cousins interval
691  if(option.Contains("fc")) {
692  option.ReplaceAll("fc","");
693  pBound = &TEfficiency::FeldmanCousins;
694  }
695  // mid-P Lancaster interval
696  if(option.Contains("midp")) {
697  option.ReplaceAll("midp","");
698  pBound = &TEfficiency::MidPInterval;
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 
1232 void 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 
1419 void 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 }
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save primitive as a C++ statement(s) on output stream out.
Definition: TH1.cxx:6753
Int_t fNpoints
Number of points <= fMaxSize.
Definition: TGraph.h:46
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Int_t GetNcells() const
Definition: TH1.h:295
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.
float xmin
Definition: THbookFile.cxx:93
Double_t * fX
[fNpoints] array of X points
Definition: TGraph.h:47
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8554
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:854
virtual void ComputeRange(Double_t &xmin, Double_t &ymin, Double_t &xmax, Double_t &ymax) const
Compute Range.
Double_t GetErrorY(Int_t bin) const
This function is called by GraphFitChisquare.
short Version_t
Definition: RtypesCore.h:61
TLine * line
float Float_t
Definition: RtypesCore.h:53
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.
const char Option_t
Definition: RtypesCore.h:62
Double_t GetErrorX(Int_t bin) const
This function is called by GraphFitChisquare.
float ymin
Definition: THbookFile.cxx:93
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
#define g(i)
Definition: RSha256.hxx:105
Double_t GetErrorYlow(Int_t i) const
Get low error on Y.
Double_t * fEXhigh
[fNpoints] array of X high errors
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
virtual void CopyAndRelease(Double_t **newarrays, Int_t ibegin, Int_t iend, Int_t obegin)
Copy and release.
R__EXTERN TStyle * gStyle
Definition: TStyle.h:407
Int_t GetLwb() const
Definition: TVectorT.h:73
static Bool_t CheckBinning(const TH1 &pass, const TH1 &total)
Checks binning for each axis.
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7389
Bool_t IsFloat() const
Returns kTRUE if string contains a floating point or integer number.
Definition: TString.cxx:1791
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4882
TVectorT.
Definition: TMatrixTBase.h:79
Double_t GetErrorYhigh(Int_t i) const
Get high error on Y.
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
virtual Int_t Merge(TCollection *list)
Adds all graphs with asymmetric errors from the collection to this graph.
#define gROOT
Definition: TROOT.h:414
Basic string class.
Definition: TString.h:131
Int_t GetNrows() const
Definition: TVectorT.h:75
#define f(i)
Definition: RSha256.hxx:104
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TList * fFunctions
Pointer to list of functions (fits and user)
Definition: TGraph.h:49
virtual void SetTitle(const char *title="")
Change (i.e.
Definition: TGraph.cxx:2237
Bool_t CtorAllocate()
Should be called from ctors after fNpoints has been set Note: This function should be called only fro...
TH1F * fHistogram
Pointer to histogram used for drawing axis.
Definition: TGraph.h:50
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition: TH1.cxx:8565
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
double beta(double x, double y)
Calculates the beta function.
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 != ie...
Definition: TGraph.cxx:695
virtual Int_t GetDimension() const
Definition: TH1.h:278
double sqrt(double)
Double_t GetErrorXhigh(Int_t i) const
Get high error on X.
TGraph with asymmetric error bars.
Double_t x[n]
Definition: legend1.C:17
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 Class()
Definition: Class.C:29
TGraph & operator=(const TGraph &)
Equal operator for this graph.
Definition: TGraph.cxx:187
virtual void SetName(const char *name="")
Set graph name.
Definition: TGraph.cxx:2221
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.
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 SetPointEXhigh(Int_t i, Double_t exh)
Set EXhigh for point i.
virtual void SetPointEYlow(Int_t i, Double_t eyl)
Set EYlow for point i.
virtual Bool_t DoMerge(const TGraph *g)
protected function to perform the merge operation of a graph with asymmetric errors ...
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
const Double_t sigma
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
virtual void Apply(TF1 *f)
Apply a function to all data points y = f(x,y)
virtual TArrayD * GetSumw2()
Definition: TH1.h:308
virtual ~TGraphAsymmErrors()
TGraphAsymmErrors default destructor.
static Double_t BetaMean(Double_t alpha, Double_t beta)
Compute the mean (average) of the beta distribution.
char * R__STRTOK_R(char *str, const char *delim, char **saveptr)
Definition: Rtypes.h:486
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...
Bool_t CtorAllocate()
In constructors set fNpoints than call this method.
Definition: TGraph.cxx:721
virtual void BayesDivide(const TH1 *pass, const TH1 *total, Option_t *opt="")
This function is only kept for backward compatibility.
Float_t GetErrorX() const
Definition: TStyle.h:175
Int_t fN
Definition: TArray.h:38
float ymax
Definition: THbookFile.cxx:93
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 void SwapPoints(Int_t pos1, Int_t pos2)
Swap points.
Definition: TGraph.cxx:2404
R__EXTERN TSystem * gSystem
Definition: TSystem.h:560
virtual void SwapPoints(Int_t pos1, Int_t pos2)
Swap points.
auto * a
Definition: textangle.C:12
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 dist...
static void SwapValues(Double_t *arr, Int_t pos1, Int_t pos2)
Swap values.
Definition: TGraph.cxx:2413
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:443
Collection abstract base class.
Definition: TCollection.h:63
virtual Bool_t DoMerge(const TGraph *g)
protected function to perform the merge operation of a graph
Definition: TGraph.cxx:2477
void Form(const char *fmt,...)
Formats a string using a printf style format descriptor.
Definition: TString.cxx:2289
unsigned int UInt_t
Definition: RtypesCore.h:42
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 Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
char * Form(const char *fmt,...)
Ssiz_t Length() const
Definition: TString.h:405
virtual void Divide(const TH1 *pass, const TH1 *total, Option_t *opt="cp")
Fill this TGraphAsymmErrors by dividing two 1-dimensional histograms pass/total.
Int_t GetN() const
Definition: TGraph.h:123
Double_t * fEYlow
[fNpoints] array of Y low errors
float xmax
Definition: THbookFile.cxx:93
TGraphErrors * gr
Definition: legend1.C:25
#define h(i)
Definition: RSha256.hxx:106
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t GetErrorXlow(Int_t i) const
Get low error on X.
virtual void SetName(const char *name)
Change the name of this histogram.
Definition: TH1.cxx:8374
Double_t At(Int_t i) const
Definition: TArrayD.h:79
static unsigned int total
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition: TH1.cxx:8576
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)
TGraphAsymmErrors & operator=(const TGraphAsymmErrors &gr)
TGraphAsymmErrors assignment operator.
#define ClassImp(name)
Definition: Rtypes.h:365
TGraphAsymmErrors()
TGraphAsymmErrors default constructor.
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
The TH1 histogram class.
Definition: TH1.h:56
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
virtual void SetPointEXlow(Int_t i, Double_t exl)
Set EXlow for point i.
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.
Mother of all ROOT objects.
Definition: TObject.h:37
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
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
Double_t * fY
[fNpoints] array of Y points
Definition: TGraph.h:48
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.
auto * l
Definition: textangle.C:4
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
void MakeZombie()
Definition: TObject.h:49
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
#define gPad
Definition: TVirtualPad.h:286
Int_t Atoi() const
Return integer value of string.
Definition: TString.cxx:1921
static Int_t CalculateScanfFields(const char *fmt)
Calculate scan fields.
Double_t Atof() const
Return floating-point value contained in string.
Definition: TString.cxx:1987
Bool_t IsDigit() const
Returns true if all characters in string are digits (0-9) or white spaces, i.e.
Definition: TString.cxx:1763
virtual Bool_t ExpandPathName(TString &path)
Expand a pathname getting rid of special shell characters like ~.
Definition: TSystem.cxx:1257
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save primitive as a C++ statement(s) on output stream out.
Int_t fMaxSize
!Current dimension of arrays fX and fY
Definition: TGraph.h:45
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
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
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:357
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
const Bool_t kTRUE
Definition: RtypesCore.h:87
virtual void Print(Option_t *chopt="") const
Print graph and errors values.
const Int_t n
Definition: legend1.C:16
Double_t * fEXlow
[fNpoints] array of X low errors
virtual void SetPointEYhigh(Int_t i, Double_t eyh)
Set EYhigh for point i.
static Double_t BetaMode(Double_t alpha, Double_t beta)
Compute the mode of the beta distribution.
Double_t * fEYhigh
[fNpoints] array of Y high errors
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
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.
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
const char * Data() const
Definition: TString.h:364
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 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.