Logo ROOT   6.10/09
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) {
101  TGraph::operator=(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->GetBinError(i+1);
250  fEYhigh[i] = fEYlow[i];
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  while (std::getline(infile, line, '\n')) {
383  if (line != "") {
384  if (line[line.size() - 1] == char(13)) { // removing DOS CR character
385  line.erase(line.end() - 1, line.end()) ;
386  }
387  token = strtok(const_cast<char*>(line.c_str()), option) ;
388  while (token != NULL && value_idx < ntokensToBeSaved) {
389  if (isTokenToBeSaved[token_idx]) {
390  token_str = TString(token) ;
391  token_str.ReplaceAll("\t", "") ;
392  if (!token_str.IsFloat()) {
393  isLineToBeSkipped = kTRUE ;
394  break ;
395  } else {
396  value[value_idx] = token_str.Atof() ;
397  value_idx++ ;
398  }
399  }
400  token = strtok(NULL, option) ; //next token
401  token_idx++ ;
402  }
403  if (!isLineToBeSkipped && value_idx > 1) { //i.e. 2,3 or 4
404  x = value[0] ;
405  y = value[1] ;
406  exl = value[2] ;
407  exh = value[3] ;
408  eyl = value[4] ;
409  eyh = value[5] ;
410  SetPoint(np, x, y) ;
411  SetPointError(np, exl, exh, eyl, eyh);
412  np++ ;
413  }
414  }
415  isLineToBeSkipped = kFALSE ;
416  token = NULL ;
417  token_idx = 0 ;
418  value_idx = 0 ;
419  }
420  Set(np) ;
421 
422  // Cleaning
423  delete [] isTokenToBeSaved ;
424  delete [] value ;
425  delete token ;
426  }
427  infile.close();
428 }
429 
430 ////////////////////////////////////////////////////////////////////////////////
431 /// TGraphAsymmErrors default destructor.
432 
434 {
435  if(fEXlow) delete [] fEXlow;
436  if(fEXhigh) delete [] fEXhigh;
437  if(fEYlow) delete [] fEYlow;
438  if(fEYhigh) delete [] fEYhigh;
439 }
440 
441 
442 ////////////////////////////////////////////////////////////////////////////////
443 /// Apply a function to all data points `y = f(x,y)`
444 ///
445 /// Errors are calculated as `eyh = f(x,y+eyh)-f(x,y)` and
446 /// `eyl = f(x,y)-f(x,y-eyl)`
447 ///
448 /// Special treatment has to be applied for the functions where the
449 /// role of "up" and "down" is reversed.
450 /// function suggested/implemented by Miroslav Helbich <helbich@mail.desy.de>
451 
453 {
454  Double_t x,y,exl,exh,eyl,eyh,eyl_new,eyh_new,fxy;
455 
456  if (fHistogram) {
457  delete fHistogram;
458  fHistogram = 0;
459  }
460  for (Int_t i=0;i<GetN();i++) {
461  GetPoint(i,x,y);
462  exl=GetErrorXlow(i);
463  exh=GetErrorXhigh(i);
464  eyl=GetErrorYlow(i);
465  eyh=GetErrorYhigh(i);
466 
467  fxy = f->Eval(x,y);
468  SetPoint(i,x,fxy);
469 
470  // in the case of the functions like y-> -1*y the roles of the
471  // upper and lower error bars is reversed
472  if (f->Eval(x,y-eyl)<f->Eval(x,y+eyh)) {
473  eyl_new = TMath::Abs(fxy - f->Eval(x,y-eyl));
474  eyh_new = TMath::Abs(f->Eval(x,y+eyh) - fxy);
475  }
476  else {
477  eyh_new = TMath::Abs(fxy - f->Eval(x,y-eyl));
478  eyl_new = TMath::Abs(f->Eval(x,y+eyh) - fxy);
479  }
480 
481  //error on x doesn't change
482  SetPointError(i,exl,exh,eyl_new,eyh_new);
483  }
484  if (gPad) gPad->Modified();
485 }
486 
487 ////////////////////////////////////////////////////////////////////////////////
488 ///This function is only kept for backward compatibility.
489 ///You should rather use the Divide method.
490 ///It calls `Divide(pass,total,"cl=0.683 b(1,1) mode")` which is equivalent to the
491 ///former BayesDivide method.
492 
494 {
495  Divide(pass,total,"cl=0.683 b(1,1) mode");
496 }
497 
498 ////////////////////////////////////////////////////////////////////////////////
499 /// Fill this TGraphAsymmErrors by dividing two 1-dimensional histograms pass/total
500 ///
501 /// This method serves two purposes:
502 ///
503 /// ### 1) calculating efficiencies:
504 ///
505 /// The assumption is that the entries in "pass" are a subset of those in
506 /// "total". That is, we create an "efficiency" graph, where each entry is
507 /// between 0 and 1, inclusive.
508 ///
509 /// If the histograms are not filled with unit weights, the number of effective
510 /// entries is used to normalise the bin contents which might lead to wrong results.
511 /// \f[
512 /// \text{effective entries} = \frac{(\sum w_{i})^{2}}{\sum w_{i}^{2}}
513 /// \f]
514 /// The points are assigned a x value at the center of each histogram bin.
515 /// The y values are \f$\text{eff} = \frac{\text{pass}}{\text{total}}\f$
516 /// for all options except for the
517 /// bayesian methods where the result depends on the chosen option.
518 ///
519 /// If the denominator becomes 0 or pass > total, the corresponding bin is
520 /// skipped.
521 ///
522 /// ### 2) calculating ratios of two Poisson means (option 'pois'):
523 ///
524 /// The two histograms are interpreted as independent Poisson processes and the ratio
525 /// \f[
526 /// \tau = \frac{n_{1}}{n_{2}} = \frac{\varepsilon}{1 - \varepsilon}
527 /// \f]
528 /// with \f$\varepsilon = \frac{n_{1}}{n_{1} + n_{2}}\f$.
529 /// The histogram 'pass' is interpreted as \f$n_{1}\f$ and the total histogram
530 /// is used for \f$n_{2}\f$.
531 ///
532 /// The (asymmetric) uncertainties of the Poisson ratio are linked to the uncertainties
533 /// of efficiency by a parameter transformation:
534 /// \f[
535 /// \Delta \tau_{low/up} = \frac{1}{(1 - \varepsilon)^{2}} \Delta \varepsilon_{low/up}
536 /// \f]
537 ///
538 /// The x errors span each histogram bin (lowedge ... lowedge+width)
539 /// The y errors depend on the chosen statistic methode which can be determined
540 /// by the options given below. For a detailed description of the used statistic
541 /// calculations please have a look at the corresponding functions!
542 ///
543 /// Options:
544 /// - v : verbose mode: prints information about the number of used bins
545 /// and calculated efficiencies with their errors
546 /// - cl=x : determine the used confidence level (0<x<1) (default is 0.683)
547 /// - cp : Clopper-Pearson interval (see TEfficiency::ClopperPearson)
548 /// - w : Wilson interval (see TEfficiency::Wilson)
549 /// - n : normal approximation propagation (see TEfficiency::Normal)
550 /// - ac : Agresti-Coull interval (see TEfficiency::AgrestiCoull)
551 /// - fc : Feldman-Cousins interval (see TEfficiency::FeldmanCousinsInterval)
552 /// - midp : Lancaster mid-P interval (see TEfficiency::MidPInterval)
553 /// - b(a,b): bayesian interval using a prior probability ~Beta(a,b); a,b > 0
554 /// (see TEfficiency::Bayesian)
555 /// - mode : use mode of posterior for Bayesian interval (default is mean)
556 /// - shortest: use shortest interval (done by default if mode is set)
557 /// - central: use central interval (done by default if mode is NOT set)
558 /// - pois: interpret histograms as poisson ratio instead of efficiency
559 /// - e0 : plot (in Bayesian case) efficiency and interval for bins where total=0
560 /// (default is to skip them)
561 ///
562 /// Note:
563 /// Unfortunately there is no straightforward approach for determining a confidence
564 /// interval for a given confidence level. The actual coverage probability of the
565 /// confidence interval oscillates significantly according to the total number of
566 /// events and the true efficiency. In order to decrease the impact of this
567 /// oscillation on the actual coverage probability a couple of approximations and
568 /// methodes has been developed. For a detailed discussion, please have a look at
569 /// this statistical paper:
570 /// http://www-stat.wharton.upenn.edu/~tcai/paper/Binomial-StatSci.pdf
571 
572 
573 void TGraphAsymmErrors::Divide(const TH1* pass, const TH1* total, Option_t *opt)
574 {
575  //check pointers
576  if(!pass || !total) {
577  Error("Divide","one of the passed pointers is zero");
578  return;
579  }
580 
581  //check dimension of histograms; only 1-dimensional ones are accepted
582  if((pass->GetDimension() > 1) || (total->GetDimension() > 1)) {
583  Error("Divide","passed histograms are not one-dimensional");
584  return;
585  }
586 
587  //check whether histograms are filled with weights -> use number of effective
588  //entries
589  Bool_t bEffective = false;
590  //compare sum of weights with sum of squares of weights
591  // re-compute here to be sure to get the right values
592  Double_t psumw = 0;
593  Double_t psumw2 = 0;
594  if (pass->GetSumw2()->fN > 0) {
595  for (int i = 0; i < pass->GetNcells(); ++i) {
596  psumw += pass->GetBinContent(i);
597  psumw2 += pass->GetSumw2()->At(i);
598  }
599  }
600  else {
601  psumw = pass->GetSumOfWeights();
602  psumw2 = psumw;
603  }
604  if (TMath::Abs(psumw - psumw2) > 1e-6)
605  bEffective = true;
606 
607  Double_t tsumw = 0;
608  Double_t tsumw2 = 0;
609  if (total->GetSumw2()->fN > 0) {
610  for (int i = 0; i < total->GetNcells(); ++i) {
611  tsumw += total->GetBinContent(i);
612  tsumw2 += total->GetSumw2()->At(i);
613  }
614  }
615  else {
616  tsumw = total->GetSumOfWeights();
617  tsumw2 = tsumw;
618  }
619  if (TMath::Abs(tsumw - tsumw2) > 1e-6)
620  bEffective = true;
621 
622 
623 
624  // we do not want to ignore the weights
625  // if (bEffective && (pass->GetSumw2()->fN == 0 || total->GetSumw2()->fN == 0) ) {
626  // 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");
627  // bEffective = false;
628  // }
629 
630  //parse option
631  TString option = opt;
632  option.ToLower();
633 
634  Bool_t bVerbose = false;
635  //pointer to function returning the boundaries of the confidence interval
636  //(is only used in the frequentist cases.)
637  Double_t (*pBound)(Double_t,Double_t,Double_t,Bool_t) = &TEfficiency::ClopperPearson; // default method
638  //confidence level
639  Double_t conf = 0.682689492137;
640  //values for bayesian statistics
641  Bool_t bIsBayesian = false;
642  Double_t alpha = 1;
643  Double_t beta = 1;
644 
645  //verbose mode
646  if(option.Contains("v")) {
647  option.ReplaceAll("v","");
648  bVerbose = true;
649  if (bEffective)
650  Info("Divide","weight will be considered in the Histogram Ratio");
651  }
652 
653 
654  //confidence level
655  if(option.Contains("cl=")) {
656  Double_t level = -1;
657  // coverity [secure_coding : FALSE]
658  sscanf(strstr(option.Data(),"cl="),"cl=%lf",&level);
659  if((level > 0) && (level < 1))
660  conf = level;
661  else
662  Warning("Divide","given confidence level %.3lf is invalid",level);
663  option.ReplaceAll("cl=","");
664  }
665 
666  //normal approximation
667  if(option.Contains("n")) {
668  option.ReplaceAll("n","");
669  pBound = &TEfficiency::Normal;
670  }
671 
672  //clopper pearson interval
673  if(option.Contains("cp")) {
674  option.ReplaceAll("cp","");
675  pBound = &TEfficiency::ClopperPearson;
676  }
677 
678  //wilson interval
679  if(option.Contains("w")) {
680  option.ReplaceAll("w","");
681  pBound = &TEfficiency::Wilson;
682  }
683 
684  //agresti coull interval
685  if(option.Contains("ac")) {
686  option.ReplaceAll("ac","");
687  pBound = &TEfficiency::AgrestiCoull;
688  }
689  // Feldman-Cousins interval
690  if(option.Contains("fc")) {
691  option.ReplaceAll("fc","");
692  pBound = &TEfficiency::FeldmanCousins;
693  }
694  // mid-P Lancaster interval
695  if(option.Contains("midp")) {
696  option.ReplaceAll("midp","");
697  pBound = &TEfficiency::MidPInterval;
698  }
699 
700  //bayesian with prior
701  if(option.Contains("b(")) {
702  Double_t a = 0;
703  Double_t b = 0;
704  sscanf(strstr(option.Data(),"b("),"b(%lf,%lf)",&a,&b);
705  if(a > 0)
706  alpha = a;
707  else
708  Warning("Divide","given shape parameter for alpha %.2lf is invalid",a);
709  if(b > 0)
710  beta = b;
711  else
712  Warning("Divide","given shape parameter for beta %.2lf is invalid",b);
713  option.ReplaceAll("b(","");
714  bIsBayesian = true;
715  }
716 
717  // use posterior mode
718  Bool_t usePosteriorMode = false;
719  if(bIsBayesian && option.Contains("mode") ) {
720  usePosteriorMode = true;
721  option.ReplaceAll("mode","");
722  }
723 
724  Bool_t plot0Bins = false;
725  if(option.Contains("e0") ) {
726  plot0Bins = true;
727  option.ReplaceAll("e0","");
728  }
729 
730  Bool_t useShortestInterval = false;
731  if (bIsBayesian && ( option.Contains("sh") || (usePosteriorMode && !option.Contains("cen") ) ) ) {
732  useShortestInterval = true;
733  }
734 
735  // interpret as Poisson ratio
736  Bool_t bPoissonRatio = false;
737  if(option.Contains("pois") ) {
738  bPoissonRatio = true;
739  option.ReplaceAll("pois","");
740  }
741 
742  // weights works only in case of Normal approximation or Bayesian for binomial interval
743  // in case of Poisson ratio we can use weights by rescaling the obtained results using the effective entries
744  if ( ( bEffective && !bPoissonRatio) && !bIsBayesian && pBound != &TEfficiency::Normal ) {
745  Warning("Divide","Histograms have weights: only Normal or Bayesian error calculation is supported");
746  Info("Divide","Using now the Normal approximation for weighted histograms");
747  }
748 
749  if(bPoissonRatio)
750  {
751  if(pass->GetDimension() != total->GetDimension()) {
752  Error("Divide","passed histograms are not of the same dimension");
753  return;
754  }
755 
756  if(!TEfficiency::CheckBinning(*pass,*total)) {
757  Error("Divide","passed histograms are not consistent");
758  return;
759  }
760  }
761  else
762  {
763  //check consistency of histograms, allowing weights
764  if(!TEfficiency::CheckConsistency(*pass,*total,"w")) {
765  Error("Divide","passed histograms are not consistent");
766  return;
767  }
768  }
769 
770  //Set the graph to have a number of points equal to the number of histogram
771  //bins
772  Int_t nbins = pass->GetNbinsX();
773  Set(nbins);
774 
775  // Ok, now set the points for each bin
776  // (Note: the TH1 bin content is shifted to the right by one:
777  // bin=0 is underflow, bin=nbins+1 is overflow.)
778 
779  //efficiency with lower and upper boundary of confidence interval
780  double eff, low, upper;
781  //this keeps track of the number of points added to the graph
782  int npoint=0;
783  //number of total and passed events
784  Double_t t = 0 , p = 0;
785  Double_t tw = 0, tw2 = 0, pw = 0, pw2 = 0, wratio = 1; // for the case of weights
786  //loop over all bins and fill the graph
787  for (Int_t b=1; b<=nbins; ++b) {
788 
789  // default value when total =0;
790  eff = 0;
791  low = 0;
792  upper = 0;
793 
794  // special case in case of weights we have to consider the sum of weights and the sum of weight squares
795  if(bEffective) {
796  tw = total->GetBinContent(b);
797  tw2 = (total->GetSumw2()->fN > 0) ? total->GetSumw2()->At(b) : tw;
798  pw = pass->GetBinContent(b);
799  pw2 = (pass->GetSumw2()->fN > 0) ? pass->GetSumw2()->At(b) : pw;
800 
801  if(bPoissonRatio)
802  {
803  // tw += pw;
804  // tw2 += pw2;
805  // compute ratio on the effective entries ( p and t)
806  // special case is when (pw=0, pw2=0) in this case we cannot get the bin weight.
807  // we use then the overall weight of the full histogram
808  if (pw == 0 && pw2 == 0)
809  p = 0;
810  else
811  p = (pw*pw)/pw2;
812 
813  if (tw == 0 && tw2 == 0)
814  t = 0;
815  else
816  t = (tw*tw)/tw2;
817 
818  if (pw > 0 && tw > 0)
819  // this is the ratio of the two bin weights ( pw/p / t/tw )
820  wratio = (pw*t)/(p * tw);
821  else if (pw == 0 && tw > 0)
822  // case p histogram has zero compute the weights from all the histogram
823  // weight of histogram - sumw2/sumw
824  wratio = (psumw2 * t) /(psumw * tw );
825  else if (tw == 0 && pw > 0)
826  // case t histogram has zero compute the weights from all the histogram
827  // weight of histogram - sumw2/sumw
828  wratio = (pw * tsumw) /(p * tsumw2 );
829  else if (p > 0)
830  wratio = pw/p; //not sure if needed
831  else {
832  // case both pw and tw are zero - we skip these bins
833  if (!plot0Bins) continue; // skip bins with total <= 0
834  }
835 
836  t += p;
837  //std::cout << p << " " << t << " " << wratio << std::endl;
838  }
839  else
840  if (tw <= 0 && !plot0Bins) continue; // skip bins with total <= 0
841 
842  // in the case of weights have the formula only for
843  // the normal and bayesian statistics (see below)
844 
845  }
846 
847  //use bin contents
848  else {
849  t = int( total->GetBinContent(b) + 0.5);
850  p = int(pass->GetBinContent(b) + 0.5);
851 
852  if(bPoissonRatio) t += p;
853 
854  if (!t && !plot0Bins) continue; // skip bins with total = 0
855  }
856 
857 
858  //using bayesian statistics
859  if(bIsBayesian) {
860  double aa,bb;
861 
862  if ((bEffective && !bPoissonRatio) && tw2 <= 0) {
863  // case of bins with zero errors
864  eff = pw/tw;
865  low = eff; upper = eff;
866  }
867  else {
868 
869  if (bEffective && !bPoissonRatio) {
870  // tw/tw2 re-normalize the weights
871  double norm = tw/tw2; // case of tw2 = 0 is treated above
872  aa = pw * norm + alpha;
873  bb = (tw - pw) * norm + beta;
874  }
875  else {
876  aa = double(p) + alpha;
877  bb = double(t-p) + beta;
878  }
879  if (usePosteriorMode)
880  eff = TEfficiency::BetaMode(aa,bb);
881  else
882  eff = TEfficiency::BetaMean(aa,bb);
883 
884  if (useShortestInterval) {
885  TEfficiency::BetaShortestInterval(conf,aa,bb,low,upper);
886  }
887  else {
888  low = TEfficiency::BetaCentralInterval(conf,aa,bb,false);
889  upper = TEfficiency::BetaCentralInterval(conf,aa,bb,true);
890  }
891  }
892  }
893  // case of non-bayesian statistics
894  else {
895  if (bEffective && !bPoissonRatio) {
896 
897  if (tw > 0) {
898 
899  eff = pw/tw;
900 
901  // use normal error calculation using variance of MLE with weights (F.James 8.5.2)
902  // this is the same formula used in ROOT for TH1::Divide("B")
903 
904  double variance = ( pw2 * (1. - 2 * eff) + tw2 * eff *eff ) / ( tw * tw) ;
905  double sigma = sqrt(variance);
906 
907  double prob = 0.5 * (1.-conf);
908  double delta = ROOT::Math::normal_quantile_c(prob, sigma);
909  low = eff - delta;
910  upper = eff + delta;
911  if (low < 0) low = 0;
912  if (upper > 1) upper = 1.;
913  }
914  }
915  else {
916  // when not using weights (all cases) or in case of Poisson ratio with weights
917  if(t)
918  eff = ((Double_t)p)/t;
919 
920  low = pBound(t,p,conf,false);
921  upper = pBound(t,p,conf,true);
922  }
923  }
924  // treat as Poisson ratio
925  if(bPoissonRatio)
926  {
927  Double_t ratio = eff/(1 - eff);
928  // take the intervals in eff as intervals in the Poisson ratio
929  low = low/(1. - low);
930  upper = upper/(1.-upper);
931  eff = ratio;
932  if (bEffective) {
933  //scale result by the ratio of the weight
934  eff *= wratio;
935  low *= wratio;
936  upper *= wratio;
937  }
938  }
939  //Set the point center and its errors
940  if (TMath::Finite(eff)) {
941  SetPoint(npoint,pass->GetBinCenter(b),eff);
942  SetPointError(npoint,
943  pass->GetBinCenter(b)-pass->GetBinLowEdge(b),
944  pass->GetBinLowEdge(b)-pass->GetBinCenter(b)+pass->GetBinWidth(b),
945  eff-low,upper-eff);
946  npoint++;//we have added a point to the graph
947  }
948  }
949 
950  Set(npoint);//tell the graph how many points we've really added
951  if (npoint < nbins)
952  Warning("Divide","Number of graph points is different than histogram bins - %d points have been skipped",nbins-npoint);
953 
954 
955  if (bVerbose) {
956  Info("Divide","made a graph with %d points from %d bins",npoint,nbins);
957  Info("Divide","used confidence level: %.2lf\n",conf);
958  if(bIsBayesian)
959  Info("Divide","used prior probability ~ beta(%.2lf,%.2lf)",alpha,beta);
960  Print();
961  }
962 }
963 
964 ////////////////////////////////////////////////////////////////////////////////
965 /// Compute Range
966 
968 {
969  TGraph::ComputeRange(xmin,ymin,xmax,ymax);
970 
971  for (Int_t i=0;i<fNpoints;i++) {
972  if (fX[i] -fEXlow[i] < xmin) {
973  if (gPad && gPad->GetLogx()) {
974  if (fEXlow[i] < fX[i]) xmin = fX[i]-fEXlow[i];
975  else xmin = TMath::Min(xmin,fX[i]/3);
976  } else {
977  xmin = fX[i]-fEXlow[i];
978  }
979  }
980  if (fX[i] +fEXhigh[i] > xmax) xmax = fX[i]+fEXhigh[i];
981  if (fY[i] -fEYlow[i] < ymin) {
982  if (gPad && gPad->GetLogy()) {
983  if (fEYlow[i] < fY[i]) ymin = fY[i]-fEYlow[i];
984  else ymin = TMath::Min(ymin,fY[i]/3);
985  } else {
986  ymin = fY[i]-fEYlow[i];
987  }
988  }
989  if (fY[i] +fEYhigh[i] > ymax) ymax = fY[i]+fEYhigh[i];
990  }
991 }
992 
993 
994 ////////////////////////////////////////////////////////////////////////////////
995 /// Copy and release.
996 
998  Int_t ibegin, Int_t iend, Int_t obegin)
999 {
1000  CopyPoints(newarrays, ibegin, iend, obegin);
1001  if (newarrays) {
1002  delete[] fEXlow;
1003  fEXlow = newarrays[0];
1004  delete[] fEXhigh;
1005  fEXhigh = newarrays[1];
1006  delete[] fEYlow;
1007  fEYlow = newarrays[2];
1008  delete[] fEYhigh;
1009  fEYhigh = newarrays[3];
1010  delete[] fX;
1011  fX = newarrays[4];
1012  delete[] fY;
1013  fY = newarrays[5];
1014  delete[] newarrays;
1015  }
1016 }
1017 
1018 
1019 ////////////////////////////////////////////////////////////////////////////////
1020 /// Copy errors from fE*** to arrays[***]
1021 /// or to f*** Copy points.
1022 
1024  Int_t ibegin, Int_t iend, Int_t obegin)
1025 {
1026  if (TGraph::CopyPoints(arrays ? arrays+4 : 0, ibegin, iend, obegin)) {
1027  Int_t n = (iend - ibegin)*sizeof(Double_t);
1028  if (arrays) {
1029  memmove(&arrays[0][obegin], &fEXlow[ibegin], n);
1030  memmove(&arrays[1][obegin], &fEXhigh[ibegin], n);
1031  memmove(&arrays[2][obegin], &fEYlow[ibegin], n);
1032  memmove(&arrays[3][obegin], &fEYhigh[ibegin], n);
1033  } else {
1034  memmove(&fEXlow[obegin], &fEXlow[ibegin], n);
1035  memmove(&fEXhigh[obegin], &fEXhigh[ibegin], n);
1036  memmove(&fEYlow[obegin], &fEYlow[ibegin], n);
1037  memmove(&fEYhigh[obegin], &fEYhigh[ibegin], n);
1038  }
1039  return kTRUE;
1040  } else {
1041  return kFALSE;
1042  }
1043 }
1044 
1045 
1046 ////////////////////////////////////////////////////////////////////////////////
1047 /// Should be called from ctors after fNpoints has been set
1048 /// Note: This function should be called only from the constructor
1049 /// since it does not delete previously existing arrays
1050 
1052 {
1053  if (!fNpoints) {
1054  fEXlow = fEYlow = fEXhigh = fEYhigh = 0;
1055  return kFALSE;
1056  }
1057  fEXlow = new Double_t[fMaxSize];
1058  fEYlow = new Double_t[fMaxSize];
1059  fEXhigh = new Double_t[fMaxSize];
1060  fEYhigh = new Double_t[fMaxSize];
1061  return kTRUE;
1062 }
1063 
1064 ////////////////////////////////////////////////////////////////////////////////
1065 /// protected function to perform the merge operation of a graph with asymmetric errors
1066 
1068 {
1069  if (g->GetN() == 0) return kFALSE;
1070 
1071  Double_t * exl = g->GetEXlow();
1072  Double_t * exh = g->GetEXhigh();
1073  Double_t * eyl = g->GetEYlow();
1074  Double_t * eyh = g->GetEYhigh();
1075  if (exl == 0 || exh == 0 || eyl == 0 || eyh == 0) {
1076  if (g->IsA() != TGraph::Class() )
1077  Warning("DoMerge","Merging a %s is not compatible with a TGraphAsymmErrors - errors will be ignored",g->IsA()->GetName());
1078  return TGraph::DoMerge(g);
1079  }
1080  for (Int_t i = 0 ; i < g->GetN(); i++) {
1081  Int_t ipoint = GetN();
1082  Double_t x = g->GetX()[i];
1083  Double_t y = g->GetY()[i];
1084  SetPoint(ipoint, x, y);
1085  SetPointError(ipoint, exl[i], exh[i], eyl[i], eyh[i] );
1086  }
1087 
1088  return kTRUE;
1089 }
1090 
1091 ////////////////////////////////////////////////////////////////////////////////
1092 /// Set zero values for point arrays in the range [begin, end)
1093 
1095  Bool_t from_ctor)
1096 {
1097  if (!from_ctor) {
1098  TGraph::FillZero(begin, end, from_ctor);
1099  }
1100  Int_t n = (end - begin)*sizeof(Double_t);
1101  memset(fEXlow + begin, 0, n);
1102  memset(fEXhigh + begin, 0, n);
1103  memset(fEYlow + begin, 0, n);
1104  memset(fEYhigh + begin, 0, n);
1105 }
1106 
1107 
1108 ////////////////////////////////////////////////////////////////////////////////
1109 /// This function is called by GraphFitChisquare.
1110 /// It returns the error along X at point i.
1111 
1113 {
1114  if (i < 0 || i >= fNpoints) return -1;
1115  if (!fEXlow && !fEXhigh) return -1;
1116  Double_t elow=0, ehigh=0;
1117  if (fEXlow) elow = fEXlow[i];
1118  if (fEXhigh) ehigh = fEXhigh[i];
1119  return TMath::Sqrt(0.5*(elow*elow + ehigh*ehigh));
1120 }
1121 
1122 
1123 ////////////////////////////////////////////////////////////////////////////////
1124 /// This function is called by GraphFitChisquare.
1125 /// It returns the error along Y at point i.
1126 
1128 {
1129  if (i < 0 || i >= fNpoints) return -1;
1130  if (!fEYlow && !fEYhigh) return -1;
1131  Double_t elow=0, ehigh=0;
1132  if (fEYlow) elow = fEYlow[i];
1133  if (fEYhigh) ehigh = fEYhigh[i];
1134  return TMath::Sqrt(0.5*(elow*elow + ehigh*ehigh));
1135 }
1136 
1137 
1138 ////////////////////////////////////////////////////////////////////////////////
1139 /// Get high error on X.
1140 
1142 {
1143  if (i<0 || i>fNpoints) return -1;
1144  if (fEXhigh) return fEXhigh[i];
1145  return -1;
1146 }
1147 
1148 
1149 ////////////////////////////////////////////////////////////////////////////////
1150 /// Get low error on X.
1151 
1153 {
1154  if (i<0 || i>fNpoints) return -1;
1155  if (fEXlow) return fEXlow[i];
1156  return -1;
1157 }
1158 
1159 
1160 ////////////////////////////////////////////////////////////////////////////////
1161 /// Get high error on Y.
1162 
1164 {
1165  if (i<0 || i>fNpoints) return -1;
1166  if (fEYhigh) return fEYhigh[i];
1167  return -1;
1168 }
1169 
1170 
1171 ////////////////////////////////////////////////////////////////////////////////
1172 /// Get low error on Y.
1173 
1175 {
1176  if (i<0 || i>fNpoints) return -1;
1177  if (fEYlow) return fEYlow[i];
1178  return -1;
1179 }
1180 
1181 
1182 ////////////////////////////////////////////////////////////////////////////////
1183 /// Adds all graphs with asymmetric errors from the collection to this graph.
1184 /// Returns the total number of poins in the result or -1 in case of an error.
1185 
1187 {
1188  TIter next(li);
1189  while (TObject* o = next()) {
1190  TGraph *g = dynamic_cast<TGraph*>(o);
1191  if (!g) {
1192  Error("Merge",
1193  "Cannot merge - an object which doesn't inherit from TGraph found in the list");
1194  return -1;
1195  }
1196  int n0 = GetN();
1197  int n1 = n0+g->GetN();
1198  Set(n1);
1199  Double_t * x = g->GetX();
1200  Double_t * y = g->GetY();
1201  Double_t * exlow = g->GetEXlow();
1202  Double_t * exhigh = g->GetEXhigh();
1203  Double_t * eylow = g->GetEYlow();
1204  Double_t * eyhigh = g->GetEYhigh();
1205  for (Int_t i = 0 ; i < g->GetN(); i++) {
1206  SetPoint(n0+i, x[i], y[i]);
1207  if (exlow) fEXlow[n0+i] = exlow[i];
1208  if (exhigh) fEXhigh[n0+i] = exhigh[i];
1209  if (eylow) fEYlow[n0+i] = eylow[i];
1210  if (eyhigh) fEYhigh[n0+i] = eyhigh[i];
1211  }
1212  }
1213  return GetN();
1214 }
1215 
1216 ////////////////////////////////////////////////////////////////////////////////
1217 /// Print graph and errors values.
1218 
1220 {
1221  for (Int_t i=0;i<fNpoints;i++) {
1222  printf("x[%d]=%g, y[%d]=%g, exl[%d]=%g, exh[%d]=%g, eyl[%d]=%g, eyh[%d]=%g\n"
1223  ,i,fX[i],i,fY[i],i,fEXlow[i],i,fEXhigh[i],i,fEYlow[i],i,fEYhigh[i]);
1224  }
1225 }
1226 
1227 
1228 ////////////////////////////////////////////////////////////////////////////////
1229 /// Save primitive as a C++ statement(s) on output stream out
1230 
1231 void TGraphAsymmErrors::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1232 {
1233  char quote = '"';
1234  out << " " << std::endl;
1235  static Int_t frameNumber = 3000;
1236  frameNumber++;
1237 
1238  Int_t i;
1239  TString fXName = TString(GetName()) + Form("_fx%d",frameNumber);
1240  TString fYName = TString(GetName()) + Form("_fy%d",frameNumber);
1241  TString fElXName = TString(GetName()) + Form("_felx%d",frameNumber);
1242  TString fElYName = TString(GetName()) + Form("_fely%d",frameNumber);
1243  TString fEhXName = TString(GetName()) + Form("_fehx%d",frameNumber);
1244  TString fEhYName = TString(GetName()) + Form("_fehy%d",frameNumber);
1245  out << " Double_t " << fXName << "[" << fNpoints << "] = {" << std::endl;
1246  for (i = 0; i < fNpoints-1; i++) out << " " << fX[i] << "," << std::endl;
1247  out << " " << fX[fNpoints-1] << "};" << std::endl;
1248  out << " Double_t " << fYName << "[" << fNpoints << "] = {" << std::endl;
1249  for (i = 0; i < fNpoints-1; i++) out << " " << fY[i] << "," << std::endl;
1250  out << " " << fY[fNpoints-1] << "};" << std::endl;
1251  out << " Double_t " << fElXName << "[" << fNpoints << "] = {" << std::endl;
1252  for (i = 0; i < fNpoints-1; i++) out << " " << fEXlow[i] << "," << std::endl;
1253  out << " " << fEXlow[fNpoints-1] << "};" << std::endl;
1254  out << " Double_t " << fElYName << "[" << fNpoints << "] = {" << std::endl;
1255  for (i = 0; i < fNpoints-1; i++) out << " " << fEYlow[i] << "," << std::endl;
1256  out << " " << fEYlow[fNpoints-1] << "};" << std::endl;
1257  out << " Double_t " << fEhXName << "[" << fNpoints << "] = {" << std::endl;
1258  for (i = 0; i < fNpoints-1; i++) out << " " << fEXhigh[i] << "," << std::endl;
1259  out << " " << fEXhigh[fNpoints-1] << "};" << std::endl;
1260  out << " Double_t " << fEhYName << "[" << fNpoints << "] = {" << std::endl;
1261  for (i = 0; i < fNpoints-1; i++) out << " " << fEYhigh[i] << "," << std::endl;
1262  out << " " << fEYhigh[fNpoints-1] << "};" << std::endl;
1263 
1264  if (gROOT->ClassSaved(TGraphAsymmErrors::Class())) out<<" ";
1265  else out << " TGraphAsymmErrors *";
1266  out << "grae = new TGraphAsymmErrors("<< fNpoints << ","
1267  << fXName << "," << fYName << ","
1268  << fElXName << "," << fEhXName << ","
1269  << fElYName << "," << fEhYName << ");"
1270  << std::endl;
1271 
1272  out << " grae->SetName(" << quote << GetName() << quote << ");" << std::endl;
1273  out << " grae->SetTitle(" << quote << GetTitle() << quote << ");" << std::endl;
1274 
1275  SaveFillAttributes(out, "grae", 0, 1001);
1276  SaveLineAttributes(out, "grae", 1, 1, 1);
1277  SaveMarkerAttributes(out, "grae", 1, 1, 1);
1278 
1279  if (fHistogram) {
1280  TString hname = fHistogram->GetName();
1281  hname += frameNumber;
1282  fHistogram->SetName(Form("Graph_%s",hname.Data()));
1283  fHistogram->SavePrimitive(out,"nodraw");
1284  out<<" grae->SetHistogram("<<fHistogram->GetName()<<");"<<std::endl;
1285  out<<" "<<std::endl;
1286  }
1287 
1288  // save list of functions
1289  TIter next(fFunctions);
1290  TObject *obj;
1291  while ((obj = next())) {
1292  obj->SavePrimitive(out, Form("nodraw #%d\n",++frameNumber));
1293  if (obj->InheritsFrom("TPaveStats")) {
1294  out << " grae->GetListOfFunctions()->Add(ptstats);" << std::endl;
1295  out << " ptstats->SetParent(grae->GetListOfFunctions());" << std::endl;
1296  } else {
1297  out << " grae->GetListOfFunctions()->Add("
1298  << Form("%s%d",obj->GetName(),frameNumber) << ");" << std::endl;
1299  }
1300  }
1301 
1302  const char *l = strstr(option,"multigraph");
1303  if (l) {
1304  out<<" multigraph->Add(grae,"<<quote<<l+10<<quote<<");"<<std::endl;
1305  } else {
1306  out<<" grae->Draw("<<quote<<option<<quote<<");"<<std::endl;
1307  }
1308 }
1309 
1310 ////////////////////////////////////////////////////////////////////////////////
1311 /// Set ex and ey values for point pointed by the mouse.
1312 
1314 {
1315  Int_t px = gPad->GetEventX();
1316  Int_t py = gPad->GetEventY();
1317 
1318  //localize point to be deleted
1319  Int_t ipoint = -2;
1320  Int_t i;
1321  // start with a small window (in case the mouse is very close to one point)
1322  for (i=0;i<fNpoints;i++) {
1323  Int_t dpx = px - gPad->XtoAbsPixel(gPad->XtoPad(fX[i]));
1324  Int_t dpy = py - gPad->YtoAbsPixel(gPad->YtoPad(fY[i]));
1325  if (dpx*dpx+dpy*dpy < 25) {ipoint = i; break;}
1326  }
1327  if (ipoint == -2) return;
1328 
1329  fEXlow[ipoint] = exl;
1330  fEYlow[ipoint] = eyl;
1331  fEXhigh[ipoint] = exh;
1332  fEYhigh[ipoint] = eyh;
1333  gPad->Modified();
1334 }
1335 
1336 
1337 ////////////////////////////////////////////////////////////////////////////////
1338 /// Set ex and ey values for point number i.
1339 
1341 {
1342  if (i < 0) return;
1343  if (i >= fNpoints) {
1344  // re-allocate the object
1346  }
1347  fEXlow[i] = exl;
1348  fEYlow[i] = eyl;
1349  fEXhigh[i] = exh;
1350  fEYhigh[i] = eyh;
1351 }
1352 
1353 
1354 ////////////////////////////////////////////////////////////////////////////////
1355 /// Set EXlow for point i
1356 
1358 {
1359  if (i < 0) return;
1360  if (i >= fNpoints) {
1361  // re-allocate the object
1363  }
1364  fEXlow[i] = exl;
1365 }
1366 
1367 
1368 ////////////////////////////////////////////////////////////////////////////////
1369 /// Set EXhigh for point i
1370 
1372 {
1373  if (i < 0) return;
1374  if (i >= fNpoints) {
1375  // re-allocate the object
1377  }
1378  fEXhigh[i] = exh;
1379 }
1380 
1381 
1382 ////////////////////////////////////////////////////////////////////////////////
1383 /// Set EYlow for point i
1384 
1386 {
1387  if (i < 0) return;
1388  if (i >= fNpoints) {
1389  // re-allocate the object
1391  }
1392  fEYlow[i] = eyl;
1393 }
1394 
1395 
1396 ////////////////////////////////////////////////////////////////////////////////
1397 /// Set EYhigh for point i
1398 
1400 {
1401  if (i < 0) return;
1402  if (i >= fNpoints) {
1403  // re-allocate the object
1405  }
1406  fEYhigh[i] = eyh;
1407 }
1408 
1409 
1410 ////////////////////////////////////////////////////////////////////////////////
1411 /// Stream an object of class TGraphAsymmErrors.
1412 
1413 void TGraphAsymmErrors::Streamer(TBuffer &b)
1414 {
1415  if (b.IsReading()) {
1416  UInt_t R__s, R__c;
1417  Version_t R__v = b.ReadVersion(&R__s, &R__c);
1418  if (R__v > 2) {
1419  b.ReadClassBuffer(TGraphAsymmErrors::Class(), this, R__v, R__s, R__c);
1420  return;
1421  }
1422  //====process old versions before automatic schema evolution
1423  TGraph::Streamer(b);
1424  fEXlow = new Double_t[fNpoints];
1425  fEYlow = new Double_t[fNpoints];
1426  fEXhigh = new Double_t[fNpoints];
1427  fEYhigh = new Double_t[fNpoints];
1428  if (R__v < 2) {
1429  Float_t *exlow = new Float_t[fNpoints];
1430  Float_t *eylow = new Float_t[fNpoints];
1431  Float_t *exhigh = new Float_t[fNpoints];
1432  Float_t *eyhigh = new Float_t[fNpoints];
1433  b.ReadFastArray(exlow,fNpoints);
1434  b.ReadFastArray(eylow,fNpoints);
1435  b.ReadFastArray(exhigh,fNpoints);
1436  b.ReadFastArray(eyhigh,fNpoints);
1437  for (Int_t i=0;i<fNpoints;i++) {
1438  fEXlow[i] = exlow[i];
1439  fEYlow[i] = eylow[i];
1440  fEXhigh[i] = exhigh[i];
1441  fEYhigh[i] = eyhigh[i];
1442  }
1443  delete [] eylow;
1444  delete [] exlow;
1445  delete [] eyhigh;
1446  delete [] exhigh;
1447  } else {
1452  }
1453  b.CheckByteCount(R__s, R__c, TGraphAsymmErrors::IsA());
1454  //====end of old versions
1455 
1456  } else {
1458  }
1459 }
1460 
1461 
1462 ////////////////////////////////////////////////////////////////////////////////
1463 /// Swap points.
1464 
1466 {
1467  SwapValues(fEXlow, pos1, pos2);
1468  SwapValues(fEXhigh, pos1, pos2);
1469  SwapValues(fEYlow, pos1, pos2);
1470  SwapValues(fEYhigh, pos1, pos2);
1471  TGraph::SwapPoints(pos1, pos2);
1472 }
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save primitive as a C++ statement(s) on output stream out.
Definition: TH1.cxx:6472
Int_t fNpoints
Number of points <= fMaxSize.
Definition: TGraph.h:46
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Bool_t IsReading() const
Definition: TBuffer.h:81
virtual Int_t GetNcells() const
Definition: TH1.h:280
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 Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8253
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:847
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...
virtual Double_t * GetEYlow() const
Definition: TGraph.h:136
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:640
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:402
Int_t GetLwb() const
Definition: TVectorT.h:73
static Bool_t CheckBinning(const TH1 &pass, const TH1 &total)
Checks binning for each axis.
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:131
TH1 * h
Definition: legend2.C:5
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7103
Bool_t IsFloat() const
Returns kTRUE if string contains a floating point or integer number.
Definition: TString.cxx:1845
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4639
TVectorT.
Definition: TMatrixTBase.h:77
Double_t GetErrorYhigh(Int_t i) const
Get high error on Y.
Buffer base class used for serializing objects.
Definition: TBuffer.h:40
virtual Int_t Merge(TCollection *list)
Adds all graphs with asymmetric errors from the collection to this graph.
virtual Int_t CheckByteCount(UInt_t startpos, UInt_t bcnt, const TClass *clss)=0
#define gROOT
Definition: TROOT.h:375
Basic string class.
Definition: TString.h:129
Int_t GetNrows() const
Definition: TVectorT.h:75
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1099
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
TList * fFunctions
Pointer to list of functions (fits and user)
Definition: TGraph.h:49
virtual void SetTitle(const char *title="")
Set graph title.
Definition: TGraph.cxx:2180
Bool_t CtorAllocate()
Should be called from ctors after fNpoints has been set Note: This function should be called only fro...
int nbins[3]
TH1F * fHistogram
Pointer to histogram used for drawing axis.
Definition: TGraph.h:50
#define NULL
Definition: RtypesCore.h:88
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition: TH1.cxx:8264
static std::string format(double x, double y, int digits, int width)
Float_t delta
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
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:693
virtual Int_t GetDimension() const
Definition: TH1.h:263
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:2345
void Class()
Definition: Class.C:29
TGraph & operator=(const TGraph &)
Equal operator for this graph.
Definition: TGraph.cxx:187
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:260
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)
Definition: TMath.h:655
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:244
virtual void Apply(TF1 *f)
Apply a function to all data points y = f(x,y)
virtual TArrayD * GetSumw2()
Definition: TH1.h:293
virtual ~TGraphAsymmErrors()
TGraphAsymmErrors default destructor.
static Double_t BetaMean(Double_t alpha, Double_t beta)
Compute the mean (average) of the 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...
Bool_t CtorAllocate()
In constructors set fNpoints than call this method.
Definition: TGraph.cxx:719
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:172
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:2327
R__EXTERN TSystem * gSystem
Definition: TSystem.h:539
virtual void SwapPoints(Int_t pos1, Int_t pos2)
Swap points.
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 void SwapValues(Double_t *arr, Int_t pos1, Int_t pos2)
Swap values.
Definition: TGraph.cxx:2336
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:436
Collection abstract base class.
Definition: TCollection.h:42
virtual Bool_t DoMerge(const TGraph *g)
protected function to perform the merge operation of a graph
Definition: TGraph.cxx:2400
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:232
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:873
char * Form(const char *fmt,...)
Ssiz_t Length() const
Definition: TString.h:388
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:122
Double_t * fEYlow
[fNpoints] array of Y low errors
TLine * l
Definition: textangle.C:4
const std::string sname
Definition: testIO.cxx:45
float xmax
Definition: THbookFile.cxx:93
virtual void ReadFastArray(Bool_t *b, Int_t n)=0
TGraphErrors * gr
Definition: legend1.C:25
const Bool_t kFALSE
Definition: RtypesCore.h:92
Double_t * GetX() const
Definition: TGraph.h:129
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:8073
Double_t At(Int_t i) const
Definition: TArrayD.h:79
static unsigned int total
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
virtual Double_t Eval(Double_t x, Double_t y=0, Double_t z=0, Double_t t=0) const
Evaluate this function.
Definition: TF1.cxx:1197
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition: TH1.cxx:8275
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:336
double f(double x)
TGraphAsymmErrors()
TGraphAsymmErrors default constructor.
double Double_t
Definition: RtypesCore.h:55
virtual Double_t * GetEXlow() const
Definition: TGraph.h:134
Double_t y[n]
Definition: legend1.C:17
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:572
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:1014
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2156
Double_t * fY
[fNpoints] array of Y points
Definition: TGraph.h:48
Double_t * GetY() const
Definition: TGraph.h:130
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 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:150
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:284
Int_t Atoi() const
Return integer value of string.
Definition: TString.cxx:1975
static Int_t CalculateScanfFields(const char *fmt)
Calculate scan fields.
Double_t Atof() const
Return floating-point value contained in string.
Definition: TString.cxx:2041
Bool_t IsDigit() const
Returns true if all characters in string are digits (0-9) or white spaces, i.e.
Definition: TString.cxx:1817
virtual Bool_t ExpandPathName(TString &path)
Expand a pathname getting rid of special shell characters like ~.
Definition: TSystem.cxx:1250
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:277
Double_t Sqrt(Double_t x)
Definition: TMath.h:591
virtual Double_t * GetEYhigh() const
Definition: TGraph.h:135
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:645
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:364
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:2105
const Bool_t kTRUE
Definition: RtypesCore.h:91
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
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.
virtual Double_t * GetEXhigh() const
Definition: TGraph.h:133
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:859
static Double_t Wilson(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries for the frequentist Wilson interval.
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
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
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8175
const char * Data() const
Definition: TString.h:347
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TObject.cxx:657
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.