Logo ROOT  
Reference Guide
TH1.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Author: Rene Brun 26/12/94
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 <cstdlib>
13 #include <cstring>
14 #include <cstdio>
15 #include <cctype>
16 #include <sstream>
17 #include <cmath>
18 #include <iostream>
19 
20 #include "TROOT.h"
21 #include "TBuffer.h"
22 #include "TEnv.h"
23 #include "TClass.h"
24 #include "TMath.h"
25 #include "THashList.h"
26 #include "TH1.h"
27 #include "TH2.h"
28 #include "TH3.h"
29 #include "TF2.h"
30 #include "TF3.h"
31 #include "TPluginManager.h"
32 #include "TVirtualPad.h"
33 #include "TRandom.h"
34 #include "TVirtualFitter.h"
35 #include "THLimitsFinder.h"
36 #include "TProfile.h"
37 #include "TStyle.h"
38 #include "TVectorF.h"
39 #include "TVectorD.h"
40 #include "TBrowser.h"
41 #include "TError.h"
42 #include "TVirtualHistPainter.h"
43 #include "TVirtualFFT.h"
44 #include "TVirtualPaveStats.h"
45 
46 #include "HFitInterface.h"
47 #include "Fit/DataRange.h"
48 #include "Fit/BinData.h"
49 #include "Math/GoFTest.h"
50 #include "Math/MinimizerOptions.h"
51 #include "Math/QuantFuncMathCore.h"
52 
53 #include "TH1Merger.h"
54 
55 /** \addtogroup Hist
56 @{
57 \class TH1C
58 \brief 1-D histogram with a byte per channel (see TH1 documentation)
59 \class TH1S
60 \brief 1-D histogram with a short per channel (see TH1 documentation)
61 \class TH1I
62 \brief 1-D histogram with an int per channel (see TH1 documentation)}
63 \class TH1F
64 \brief 1-D histogram with a float per channel (see TH1 documentation)}
65 \class TH1D
66 \brief 1-D histogram with a double per channel (see TH1 documentation)}
67 @}
68 */
69 
70 /** \class TH1
71 The TH1 histogram class.
72 
73 ### The Histogram classes
74 ROOT supports the following histogram types:
75 
76  - 1-D histograms:
77  - TH1C : histograms with one byte per channel. Maximum bin content = 127
78  - TH1S : histograms with one short per channel. Maximum bin content = 32767
79  - TH1I : histograms with one int per channel. Maximum bin content = 2147483647
80  - TH1F : histograms with one float per channel. Maximum precision 7 digits
81  - TH1D : histograms with one double per channel. Maximum precision 14 digits
82  - 2-D histograms:
83  - TH2C : histograms with one byte per channel. Maximum bin content = 127
84  - TH2S : histograms with one short per channel. Maximum bin content = 32767
85  - TH2I : histograms with one int per channel. Maximum bin content = 2147483647
86  - TH2F : histograms with one float per channel. Maximum precision 7 digits
87  - TH2D : histograms with one double per channel. Maximum precision 14 digits
88  - 3-D histograms:
89  - TH3C : histograms with one byte per channel. Maximum bin content = 127
90  - TH3S : histograms with one short per channel. Maximum bin content = 32767
91  - TH3I : histograms with one int per channel. Maximum bin content = 2147483647
92  - TH3F : histograms with one float per channel. Maximum precision 7 digits
93  - TH3D : histograms with one double per channel. Maximum precision 14 digits
94  - Profile histograms: See classes TProfile, TProfile2D and TProfile3D.
95  Profile histograms are used to display the mean value of Y and its standard deviation
96  for each bin in X. Profile histograms are in many cases an elegant
97  replacement of two-dimensional histograms : the inter-relation of two
98  measured quantities X and Y can always be visualized by a two-dimensional
99  histogram or scatter-plot; If Y is an unknown (but single-valued)
100  approximate function of X, this function is displayed by a profile
101  histogram with much better precision than by a scatter-plot.
102 
103 
104 All histogram classes are derived from the base class TH1
105 ~~~ {.cpp}
106  TH1
107  ^
108  |
109  |
110  |
111  +----------------+-------+------+------+-----+-----+
112  | | | | | | |
113  | | TH1C TH1S TH1I TH1F TH1D
114  | | |
115  | | |
116  | TH2 TProfile
117  | |
118  | |
119  | +-------+------+------+-----+-----+
120  | | | | | |
121  | TH2C TH2S TH2I TH2F TH2D
122  | |
123  TH3 |
124  | TProfile2D
125  |
126  +-------+------+------+------+------+
127  | | | | |
128  TH3C TH3S TH3I TH3F TH3D
129  |
130  |
131  TProfile3D
132 
133  The TH*C classes also inherit from the array class TArrayC.
134  The TH*S classes also inherit from the array class TArrayS.
135  The TH*I classes also inherit from the array class TArrayI.
136  The TH*F classes also inherit from the array class TArrayF.
137  The TH*D classes also inherit from the array class TArrayD.
138 ~~~
139 
140 #### Creating histograms
141 
142 Histograms are created by invoking one of the constructors, e.g.
143 ~~~ {.cpp}
144  TH1F *h1 = new TH1F("h1", "h1 title", 100, 0, 4.4);
145  TH2F *h2 = new TH2F("h2", "h2 title", 40, 0, 4, 30, -3, 3);
146 ~~~
147 Histograms may also be created by:
148 
149  - calling the Clone function, see below
150  - making a projection from a 2-D or 3-D histogram, see below
151  - reading an histogram from a file
152 
153  When an histogram is created, a reference to it is automatically added
154  to the list of in-memory objects for the current file or directory.
155  This default behaviour can be changed by:
156 ~~~ {.cpp}
157  h->SetDirectory(0); for the current histogram h
158  TH1::AddDirectory(kFALSE); sets a global switch disabling the reference
159 ~~~
160  When the histogram is deleted, the reference to it is removed from
161  the list of objects in memory.
162  When a file is closed, all histograms in memory associated with this file
163  are automatically deleted.
164 
165 #### Fix or variable bin size
166 
167  All histogram types support either fix or variable bin sizes.
168  2-D histograms may have fix size bins along X and variable size bins
169  along Y or vice-versa. The functions to fill, manipulate, draw or access
170  histograms are identical in both cases.
171 
172  Each histogram always contains 3 objects TAxis: fXaxis, fYaxis and fZaxis
173  o access the axis parameters, do:
174 ~~~ {.cpp}
175  TAxis *xaxis = h->GetXaxis(); etc.
176  Double_t binCenter = xaxis->GetBinCenter(bin), etc.
177 ~~~
178  See class TAxis for a description of all the access functions.
179  The axis range is always stored internally in double precision.
180 
181 #### Convention for numbering bins
182 
183  For all histogram types: nbins, xlow, xup
184 ~~~ {.cpp}
185  bin = 0; underflow bin
186  bin = 1; first bin with low-edge xlow INCLUDED
187  bin = nbins; last bin with upper-edge xup EXCLUDED
188  bin = nbins+1; overflow bin
189 ~~~
190  In case of 2-D or 3-D histograms, a "global bin" number is defined.
191  For example, assuming a 3-D histogram with (binx, biny, binz), the function
192 ~~~ {.cpp}
193  Int_t gbin = h->GetBin(binx, biny, binz);
194 ~~~
195  returns a global/linearized gbin number. This global gbin is useful
196  to access the bin content/error information independently of the dimension.
197  Note that to access the information other than bin content and errors
198  one should use the TAxis object directly with e.g.:
199 ~~~ {.cpp}
200  Double_t xcenter = h3->GetZaxis()->GetBinCenter(27);
201 ~~~
202  returns the center along z of bin number 27 (not the global bin)
203  in the 3-D histogram h3.
204 
205 #### Alphanumeric Bin Labels
206 
207  By default, an histogram axis is drawn with its numeric bin labels.
208  One can specify alphanumeric labels instead with:
209 
210  - call TAxis::SetBinLabel(bin, label);
211  This can always be done before or after filling.
212  When the histogram is drawn, bin labels will be automatically drawn.
213  See examples labels1.C and labels2.C
214  - call to a Fill function with one of the arguments being a string, e.g.
215 ~~~ {.cpp}
216  hist1->Fill(somename, weight);
217  hist2->Fill(x, somename, weight);
218  hist2->Fill(somename, y, weight);
219  hist2->Fill(somenamex, somenamey, weight);
220 ~~~
221  See examples hlabels1.C and hlabels2.C
222  - via TTree::Draw. see for example cernstaff.C
223 ~~~ {.cpp}
224  tree.Draw("Nation::Division");
225 ~~~
226  where "Nation" and "Division" are two branches of a Tree.
227 
228 When using the options 2 or 3 above, the labels are automatically
229  added to the list (THashList) of labels for a given axis.
230  By default, an axis is drawn with the order of bins corresponding
231  to the filling sequence. It is possible to reorder the axis
232 
233  - alphabetically
234  - by increasing or decreasing values
235 
236  The reordering can be triggered via the TAxis context menu by selecting
237  the menu item "LabelsOption" or by calling directly
238  TH1::LabelsOption(option, axis) where
239 
240  - axis may be "X", "Y" or "Z"
241  - option may be:
242  - "a" sort by alphabetic order
243  - ">" sort by decreasing values
244  - "<" sort by increasing values
245  - "h" draw labels horizontal
246  - "v" draw labels vertical
247  - "u" draw labels up (end of label right adjusted)
248  - "d" draw labels down (start of label left adjusted)
249 
250  When using the option 2 above, new labels are added by doubling the current
251  number of bins in case one label does not exist yet.
252  When the Filling is terminated, it is possible to trim the number
253  of bins to match the number of active labels by calling
254 ~~~ {.cpp}
255  TH1::LabelsDeflate(axis) with axis = "X", "Y" or "Z"
256 ~~~
257  This operation is automatic when using TTree::Draw.
258  Once bin labels have been created, they become persistent if the histogram
259  is written to a file or when generating the C++ code via SavePrimitive.
260 
261 #### Histograms with automatic bins
262 
263  When an histogram is created with an axis lower limit greater or equal
264  to its upper limit, the SetBuffer is automatically called with an
265  argument fBufferSize equal to fgBufferSize (default value=1000).
266  fgBufferSize may be reset via the static function TH1::SetDefaultBufferSize.
267  The axis limits will be automatically computed when the buffer will
268  be full or when the function BufferEmpty is called.
269 
270 #### Filling histograms
271 
272  An histogram is typically filled with statements like:
273 ~~~ {.cpp}
274  h1->Fill(x);
275  h1->Fill(x, w); //fill with weight
276  h2->Fill(x, y)
277  h2->Fill(x, y, w)
278  h3->Fill(x, y, z)
279  h3->Fill(x, y, z, w)
280 ~~~
281  or via one of the Fill functions accepting names described above.
282  The Fill functions compute the bin number corresponding to the given
283  x, y or z argument and increment this bin by the given weight.
284  The Fill functions return the bin number for 1-D histograms or global
285  bin number for 2-D and 3-D histograms.
286  If TH1::Sumw2 has been called before filling, the sum of squares of
287  weights is also stored.
288  One can also increment directly a bin number via TH1::AddBinContent
289  or replace the existing content via TH1::SetBinContent.
290  To access the bin content of a given bin, do:
291 ~~~ {.cpp}
292  Double_t binContent = h->GetBinContent(bin);
293 ~~~
294 
295  By default, the bin number is computed using the current axis ranges.
296  If the automatic binning option has been set via
297 ~~~ {.cpp}
298  h->SetCanExtend(TH1::kAllAxes);
299 ~~~
300  then, the Fill Function will automatically extend the axis range to
301  accomodate the new value specified in the Fill argument. The method
302  used is to double the bin size until the new value fits in the range,
303  merging bins two by two. This automatic binning options is extensively
304  used by the TTree::Draw function when histogramming Tree variables
305  with an unknown range.
306  This automatic binning option is supported for 1-D, 2-D and 3-D histograms.
307 
308  During filling, some statistics parameters are incremented to compute
309  the mean value and Root Mean Square with the maximum precision.
310 
311  In case of histograms of type TH1C, TH1S, TH2C, TH2S, TH3C, TH3S
312  a check is made that the bin contents do not exceed the maximum positive
313  capacity (127 or 32767). Histograms of all types may have positive
314  or/and negative bin contents.
315 
316 #### Rebinning
317  At any time, an histogram can be rebinned via TH1::Rebin. This function
318  returns a new histogram with the rebinned contents.
319  If bin errors were stored, they are recomputed during the rebinning.
320 
321 #### Associated errors
322  By default, for each bin, the sum of weights is computed at fill time.
323  One can also call TH1::Sumw2 to force the storage and computation
324  of the sum of the square of weights per bin.
325  If Sumw2 has been called, the error per bin is computed as the
326  sqrt(sum of squares of weights), otherwise the error is set equal
327  to the sqrt(bin content).
328  To return the error for a given bin number, do:
329 ~~~ {.cpp}
330  Double_t error = h->GetBinError(bin);
331 ~~~
332 
333 #### Associated functions
334  One or more object (typically a TF1*) can be added to the list
335  of functions (fFunctions) associated to each histogram.
336  When TH1::Fit is invoked, the fitted function is added to this list.
337  Given an histogram h, one can retrieve an associated function
338  with:
339 ~~~ {.cpp}
340  TF1 *myfunc = h->GetFunction("myfunc");
341 ~~~
342 
343 #### Operations on histograms
344 
345  Many types of operations are supported on histograms or between histograms
346 
347  - Addition of an histogram to the current histogram.
348  - Additions of two histograms with coefficients and storage into the current
349  histogram.
350  - Multiplications and Divisions are supported in the same way as additions.
351  - The Add, Divide and Multiply functions also exist to add, divide or multiply
352  an histogram by a function.
353 
354  If an histogram has associated error bars (TH1::Sumw2 has been called),
355  the resulting error bars are also computed assuming independent histograms.
356  In case of divisions, Binomial errors are also supported.
357  One can mark a histogram to be an "average" histogram by setting its bit kIsAverage via
358  myhist.SetBit(TH1::kIsAverage);
359  When adding (see TH1::Add) average histograms, the histograms are averaged and not summed.
360 
361 #### Fitting histograms
362 
363  Histograms (1-D, 2-D, 3-D and Profiles) can be fitted with a user
364  specified function via TH1::Fit. When an histogram is fitted, the
365  resulting function with its parameters is added to the list of functions
366  of this histogram. If the histogram is made persistent, the list of
367  associated functions is also persistent. Given a pointer (see above)
368  to an associated function myfunc, one can retrieve the function/fit
369  parameters with calls such as:
370 ~~~ {.cpp}
371  Double_t chi2 = myfunc->GetChisquare();
372  Double_t par0 = myfunc->GetParameter(0); value of 1st parameter
373  Double_t err0 = myfunc->GetParError(0); error on first parameter
374 ~~~
375 
376 #### Projections of histograms
377 
378  One can:
379 
380  - make a 1-D projection of a 2-D histogram or Profile
381  see functions TH2::ProjectionX,Y, TH2::ProfileX,Y, TProfile::ProjectionX
382  - make a 1-D, 2-D or profile out of a 3-D histogram
383  see functions TH3::ProjectionZ, TH3::Project3D.
384 
385  One can fit these projections via:
386 ~~~ {.cpp}
387  TH2::FitSlicesX,Y, TH3::FitSlicesZ.
388 ~~~
389 
390 #### Random Numbers and histograms
391 
392  TH1::FillRandom can be used to randomly fill an histogram using
393  the contents of an existing TF1 function or another
394  TH1 histogram (for all dimensions).
395  For example the following two statements create and fill an histogram
396  10000 times with a default gaussian distribution of mean 0 and sigma 1:
397 ~~~ {.cpp}
398  TH1F h1("h1", "histo from a gaussian", 100, -3, 3);
399  h1.FillRandom("gaus", 10000);
400 ~~~
401  TH1::GetRandom can be used to return a random number distributed
402  according the contents of an histogram.
403 
404 #### Making a copy of an histogram
405  Like for any other ROOT object derived from TObject, one can use
406  the Clone() function. This makes an identical copy of the original
407  histogram including all associated errors and functions, e.g.:
408 ~~~ {.cpp}
409  TH1F *hnew = (TH1F*)h->Clone("hnew");
410 ~~~
411 
412 #### Normalizing histograms
413 
414  One can scale an histogram such that the bins integral is equal to
415  the normalization parameter via TH1::Scale(Double_t norm), where norm
416  is the desired normalization divided by the integral of the histogram.
417 
418 #### Drawing histograms
419 
420  Histograms are drawn via the THistPainter class. Each histogram has
421  a pointer to its own painter (to be usable in a multithreaded program).
422  Many drawing options are supported.
423  See THistPainter::Paint() for more details.
424 
425  The same histogram can be drawn with different options in different pads.
426  When an histogram drawn in a pad is deleted, the histogram is
427  automatically removed from the pad or pads where it was drawn.
428  If an histogram is drawn in a pad, then filled again, the new status
429  of the histogram will be automatically shown in the pad next time
430  the pad is updated. One does not need to redraw the histogram.
431  To draw the current version of an histogram in a pad, one can use
432 ~~~ {.cpp}
433  h->DrawCopy();
434 ~~~
435  This makes a clone (see Clone below) of the histogram. Once the clone
436  is drawn, the original histogram may be modified or deleted without
437  affecting the aspect of the clone.
438 
439  One can use TH1::SetMaximum() and TH1::SetMinimum() to force a particular
440  value for the maximum or the minimum scale on the plot. (For 1-D
441  histograms this means the y-axis, while for 2-D histograms these
442  functions affect the z-axis).
443 
444  TH1::UseCurrentStyle() can be used to change all histogram graphics
445  attributes to correspond to the current selected style.
446  This function must be called for each histogram.
447  In case one reads and draws many histograms from a file, one can force
448  the histograms to inherit automatically the current graphics style
449  by calling before gROOT->ForceStyle().
450 
451 #### Setting Drawing histogram contour levels (2-D hists only)
452 
453  By default contours are automatically generated at equidistant
454  intervals. A default value of 20 levels is used. This can be modified
455  via TH1::SetContour() or TH1::SetContourLevel().
456  the contours level info is used by the drawing options "cont", "surf",
457  and "lego".
458 
459 #### Setting histogram graphics attributes
460 
461  The histogram classes inherit from the attribute classes:
462  TAttLine, TAttFill, and TAttMarker.
463  See the member functions of these classes for the list of options.
464 
465 #### Giving titles to the X, Y and Z axis
466 
467 ~~~ {.cpp}
468  h->GetXaxis()->SetTitle("X axis title");
469  h->GetYaxis()->SetTitle("Y axis title");
470 ~~~
471  The histogram title and the axis titles can be any TLatex string.
472  The titles are part of the persistent histogram.
473  It is also possible to specify the histogram title and the axis
474  titles at creation time. These titles can be given in the "title"
475  parameter. They must be separated by ";":
476 ~~~ {.cpp}
477  TH1F* h=new TH1F("h", "Histogram title;X Axis;Y Axis;Z Axis", 100, 0, 1);
478 ~~~
479  Any title can be omitted:
480 ~~~ {.cpp}
481  TH1F* h=new TH1F("h", "Histogram title;;Y Axis", 100, 0, 1);
482  TH1F* h=new TH1F("h", ";;Y Axis", 100, 0, 1);
483 ~~~
484  The method SetTitle has the same syntax:
485 ~~~ {.cpp}
486  h->SetTitle("Histogram title;Another X title Axis");
487 ~~~
488 
489 #### Saving/Reading histograms to/from a ROOT file
490 
491  The following statements create a ROOT file and store an histogram
492  on the file. Because TH1 derives from TNamed, the key identifier on
493  the file is the histogram name:
494 ~~~ {.cpp}
495  TFile f("histos.root", "new");
496  TH1F h1("hgaus", "histo from a gaussian", 100, -3, 3);
497  h1.FillRandom("gaus", 10000);
498  h1->Write();
499 ~~~
500  To read this histogram in another Root session, do:
501 ~~~ {.cpp}
502  TFile f("histos.root");
503  TH1F *h = (TH1F*)f.Get("hgaus");
504 ~~~
505  One can save all histograms in memory to the file by:
506 ~~~ {.cpp}
507  file->Write();
508 ~~~
509 
510 #### Miscellaneous operations
511 
512 ~~~ {.cpp}
513  TH1::KolmogorovTest(): statistical test of compatibility in shape
514  between two histograms
515  TH1::Smooth() smooths the bin contents of a 1-d histogram
516  TH1::Integral() returns the integral of bin contents in a given bin range
517  TH1::GetMean(int axis) returns the mean value along axis
518  TH1::GetStdDev(int axis) returns the sigma distribution along axis
519  TH1::GetEntries() returns the number of entries
520  TH1::Reset() resets the bin contents and errors of an histogram
521 ~~~
522 */
523 
524 TF1 *gF1=0; //left for back compatibility (use TVirtualFitter::GetUserFunc instead)
525 
526 Int_t TH1::fgBufferSize = 1000;
530 
531 extern void H1InitGaus();
532 extern void H1InitExpo();
533 extern void H1InitPolynom();
534 extern void H1LeastSquareFit(Int_t n, Int_t m, Double_t *a);
535 extern void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail);
536 extern void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b);
537 
538 // Internal exceptions for the CheckConsistency method
539 class DifferentDimension: public std::exception {};
540 class DifferentNumberOfBins: public std::exception {};
541 class DifferentAxisLimits: public std::exception {};
542 class DifferentBinLimits: public std::exception {};
543 class DifferentLabels: public std::exception {};
544 
545 ClassImp(TH1);
546 
547 ////////////////////////////////////////////////////////////////////////////////
548 /// Histogram default constructor.
549 
551 {
552  fDirectory = 0;
553  fFunctions = new TList;
554  fNcells = 0;
555  fIntegral = 0;
556  fPainter = 0;
557  fEntries = 0;
558  fNormFactor = 0;
560  fMaximum = -1111;
561  fMinimum = -1111;
562  fBufferSize = 0;
563  fBuffer = 0;
565  fStatOverflows = EStatOverflows::kNeutral;
566  fXaxis.SetName("xaxis");
567  fYaxis.SetName("yaxis");
568  fZaxis.SetName("zaxis");
569  fXaxis.SetParent(this);
570  fYaxis.SetParent(this);
571  fZaxis.SetParent(this);
572  UseCurrentStyle();
573 }
574 
575 ////////////////////////////////////////////////////////////////////////////////
576 /// Histogram default destructor.
577 
579 {
580  if (!TestBit(kNotDeleted)) {
581  return;
582  }
583  delete[] fIntegral;
584  fIntegral = 0;
585  delete[] fBuffer;
586  fBuffer = 0;
587  if (fFunctions) {
589 
591  TObject* obj = 0;
592  //special logic to support the case where the same object is
593  //added multiple times in fFunctions.
594  //This case happens when the same object is added with different
595  //drawing modes
596  //In the loop below we must be careful with objects (eg TCutG) that may
597  // have been added to the list of functions of several histograms
598  //and may have been already deleted.
599  while ((obj = fFunctions->First())) {
600  while(fFunctions->Remove(obj)) { }
601  if (!obj->TestBit(kNotDeleted)) {
602  break;
603  }
604  delete obj;
605  obj = 0;
606  }
607  delete fFunctions;
608  fFunctions = 0;
609  }
610  if (fDirectory) {
611  fDirectory->Remove(this);
612  fDirectory = 0;
613  }
614  delete fPainter;
615  fPainter = 0;
616 }
617 
618 ////////////////////////////////////////////////////////////////////////////////
619 /// Normal constructor for fix bin size histograms.
620 /// Creates the main histogram structure.
621 ///
622 /// \param[in] name name of histogram (avoid blanks)
623 /// \param[in] title histogram title.
624 /// If title is of the form `stringt;stringx;stringy;stringz`,
625 /// the histogram title is set to `stringt`,
626 /// the x axis title to `stringx`, the y axis title to `stringy`, etc.
627 /// \param[in] nbins number of bins
628 /// \param[in] xlow low edge of first bin
629 /// \param[in] xup upper edge of last bin (not included in last bin)
630 ///
631 /// When an histogram is created, it is automatically added to the list
632 /// of special objects in the current directory.
633 /// To find the pointer to this histogram in the current directory
634 /// by its name, do:
635 /// ~~~ {.cpp}
636 /// TH1F *h1 = (TH1F*)gDirectory->FindObject(name);
637 /// ~~~
638 
639 TH1::TH1(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
640  :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
641 {
642  Build();
643  if (nbins <= 0) {Warning("TH1","nbins is <=0 - set to nbins = 1"); nbins = 1; }
644  fXaxis.Set(nbins,xlow,xup);
645  fNcells = fXaxis.GetNbins()+2;
646 }
647 
648 ////////////////////////////////////////////////////////////////////////////////
649 /// Normal constructor for variable bin size histograms.
650 /// Creates the main histogram structure.
651 ///
652 /// \param[in] name name of histogram (avoid blanks)
653 /// \param[in] title histogram title.
654 /// If title is of the form `stringt;stringx;stringy;stringz`
655 /// the histogram title is set to `stringt`,
656 /// the x axis title to `stringx`, the y axis title to `stringy`, etc.
657 /// \param[in] nbins number of bins
658 /// \param[in] xbins array of low-edges for each bin.
659 /// This is an array of size nbins+1
660 
661 TH1::TH1(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
662  :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
663 {
664  Build();
665  if (nbins <= 0) {Warning("TH1","nbins is <=0 - set to nbins = 1"); nbins = 1; }
666  if (xbins) fXaxis.Set(nbins,xbins);
667  else fXaxis.Set(nbins,0,1);
668  fNcells = fXaxis.GetNbins()+2;
669 }
670 
671 ////////////////////////////////////////////////////////////////////////////////
672 /// Normal constructor for variable bin size histograms.
673 ///
674 /// \param[in] name name of histogram (avoid blanks)
675 /// \param[in] title histogram title.
676 /// If title is of the form `stringt;stringx;stringy;stringz`
677 /// the histogram title is set to `stringt`,
678 /// the x axis title to `stringx`, the y axis title to `stringy`, etc.
679 /// \param[in] nbins number of bins
680 /// \param[in] xbins array of low-edges for each bin.
681 /// This is an array of size nbins+1
682 
683 TH1::TH1(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
684  :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
685 {
686  Build();
687  if (nbins <= 0) {Warning("TH1","nbins is <=0 - set to nbins = 1"); nbins = 1; }
688  if (xbins) fXaxis.Set(nbins,xbins);
689  else fXaxis.Set(nbins,0,1);
690  fNcells = fXaxis.GetNbins()+2;
691 }
692 
693 ////////////////////////////////////////////////////////////////////////////////
694 /// Copy constructor.
695 /// The list of functions is not copied. (Use Clone if needed)
696 
698 {
699  ((TH1&)h).Copy(*this);
700 }
701 
702 ////////////////////////////////////////////////////////////////////////////////
703 /// Static function: cannot be inlined on Windows/NT.
704 
706 {
707  return fgAddDirectory;
708 }
709 
710 ////////////////////////////////////////////////////////////////////////////////
711 /// Browse the Histogram object.
712 
714 {
715  Draw(b ? b->GetDrawOption() : "");
716  gPad->Update();
717 }
718 
719 ////////////////////////////////////////////////////////////////////////////////
720 /// Creates histogram basic data structure.
721 
723 {
724  fDirectory = 0;
725  fPainter = 0;
726  fIntegral = 0;
727  fEntries = 0;
728  fNormFactor = 0;
730  fMaximum = -1111;
731  fMinimum = -1111;
732  fBufferSize = 0;
733  fBuffer = 0;
735  fStatOverflows = EStatOverflows::kNeutral;
736  fXaxis.SetName("xaxis");
737  fYaxis.SetName("yaxis");
738  fZaxis.SetName("zaxis");
739  fYaxis.Set(1,0.,1.);
740  fZaxis.Set(1,0.,1.);
741  fXaxis.SetParent(this);
742  fYaxis.SetParent(this);
743  fZaxis.SetParent(this);
744 
745  SetTitle(fTitle.Data());
746 
747  fFunctions = new TList;
748 
749  UseCurrentStyle();
750 
751  if (TH1::AddDirectoryStatus()) {
753  if (fDirectory) {
755  fDirectory->Append(this,kTRUE);
756  }
757  }
758 }
759 
760 ////////////////////////////////////////////////////////////////////////////////
761 /// Performs the operation: `this = this + c1*f1`
762 /// if errors are defined (see TH1::Sumw2), errors are also recalculated.
763 ///
764 /// By default, the function is computed at the centre of the bin.
765 /// if option "I" is specified (1-d histogram only), the integral of the
766 /// function in each bin is used instead of the value of the function at
767 /// the centre of the bin.
768 ///
769 /// Only bins inside the function range are recomputed.
770 ///
771 /// IMPORTANT NOTE: If you intend to use the errors of this histogram later
772 /// you should call Sumw2 before making this operation.
773 /// This is particularly important if you fit the histogram after TH1::Add
774 ///
775 /// The function return kFALSE if the Add operation failed
776 
778 {
779  if (!f1) {
780  Error("Add","Attempt to add a non-existing function");
781  return kFALSE;
782  }
783 
784  TString opt = option;
785  opt.ToLower();
786  Bool_t integral = kFALSE;
787  if (opt.Contains("i") && fDimension == 1) integral = kTRUE;
788 
789  Int_t ncellsx = GetNbinsX() + 2; // cells = normal bins + underflow bin + overflow bin
790  Int_t ncellsy = GetNbinsY() + 2;
791  Int_t ncellsz = GetNbinsZ() + 2;
792  if (fDimension < 2) ncellsy = 1;
793  if (fDimension < 3) ncellsz = 1;
794 
795  // delete buffer if it is there since it will become invalid
796  if (fBuffer) BufferEmpty(1);
797 
798  // - Add statistics
799  Double_t s1[10];
800  for (Int_t i = 0; i < 10; ++i) s1[i] = 0;
801  PutStats(s1);
802  SetMinimum();
803  SetMaximum();
804 
805  // - Loop on bins (including underflows/overflows)
806  Int_t bin, binx, biny, binz;
807  Double_t cu=0;
808  Double_t xx[3];
809  Double_t *params = 0;
810  f1->InitArgs(xx,params);
811  for (binz = 0; binz < ncellsz; ++binz) {
812  xx[2] = fZaxis.GetBinCenter(binz);
813  for (biny = 0; biny < ncellsy; ++biny) {
814  xx[1] = fYaxis.GetBinCenter(biny);
815  for (binx = 0; binx < ncellsx; ++binx) {
816  xx[0] = fXaxis.GetBinCenter(binx);
817  if (!f1->IsInside(xx)) continue;
819  bin = binx + ncellsx * (biny + ncellsy * binz);
820  if (integral) {
821  cu = c1*f1->Integral(fXaxis.GetBinLowEdge(binx), fXaxis.GetBinUpEdge(binx), 0.) / fXaxis.GetBinWidth(binx);
822  } else {
823  cu = c1*f1->EvalPar(xx);
824  }
825  if (TF1::RejectedPoint()) continue;
826  AddBinContent(bin,cu);
827  }
828  }
829  }
830 
831  return kTRUE;
832 }
833 
834 ////////////////////////////////////////////////////////////////////////////////
835 /// Performs the operation: `this = this + c1*h1`
836 /// If errors are defined (see TH1::Sumw2), errors are also recalculated.
837 ///
838 /// Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
839 /// if not already set.
840 ///
841 /// Note also that adding histogram with labels is not supported, histogram will be
842 /// added merging them by bin number independently of the labels.
843 /// For adding histogram with labels one should use TH1::Merge
844 ///
845 /// SPECIAL CASE (Average/Efficiency histograms)
846 /// For histograms representing averages or efficiencies, one should compute the average
847 /// of the two histograms and not the sum. One can mark a histogram to be an average
848 /// histogram by setting its bit kIsAverage with
849 /// myhist.SetBit(TH1::kIsAverage);
850 /// Note that the two histograms must have their kIsAverage bit set
851 ///
852 /// IMPORTANT NOTE1: If you intend to use the errors of this histogram later
853 /// you should call Sumw2 before making this operation.
854 /// This is particularly important if you fit the histogram after TH1::Add
855 ///
856 /// IMPORTANT NOTE2: if h1 has a normalisation factor, the normalisation factor
857 /// is used , ie this = this + c1*factor*h1
858 /// Use the other TH1::Add function if you do not want this feature
859 ///
860 /// The function return kFALSE if the Add operation failed
861 
863 {
864  if (!h1) {
865  Error("Add","Attempt to add a non-existing histogram");
866  return kFALSE;
867  }
868 
869  // delete buffer if it is there since it will become invalid
870  if (fBuffer) BufferEmpty(1);
871 
872  bool useMerge = (c1 == 1. && !this->TestBit(kIsAverage) && !h1->TestBit(kIsAverage) );
873  try {
874  CheckConsistency(this,h1);
875  useMerge = kFALSE;
876  } catch(DifferentNumberOfBins&) {
877  if (useMerge)
878  Info("Add","Attempt to add histograms with different number of bins - trying to use TH1::Merge");
879  else {
880  Error("Add","Attempt to add histograms with different number of bins : nbins h1 = %d , nbins h2 = %d",GetNbinsX(), h1->GetNbinsX());
881  return kFALSE;
882  }
883  } catch(DifferentAxisLimits&) {
884  if (useMerge)
885  Info("Add","Attempt to add histograms with different axis limits - trying to use TH1::Merge");
886  else
887  Warning("Add","Attempt to add histograms with different axis limits");
888  } catch(DifferentBinLimits&) {
889  if (useMerge)
890  Info("Add","Attempt to add histograms with different bin limits - trying to use TH1::Merge");
891  else
892  Warning("Add","Attempt to add histograms with different bin limits");
893  } catch(DifferentLabels&) {
894  // in case of different labels -
895  if (useMerge)
896  Info("Add","Attempt to add histograms with different labels - trying to use TH1::Merge");
897  else
898  Info("Warning","Attempt to add histograms with different labels");
899  }
900 
901  if (useMerge) {
902  TList l;
903  l.Add(const_cast<TH1*>(h1));
904  auto iret = Merge(&l);
905  return (iret >= 0);
906  }
907 
908  // Create Sumw2 if h1 has Sumw2 set
909  if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
910 
911  // - Add statistics
912  Double_t entries = TMath::Abs( GetEntries() + c1 * h1->GetEntries() );
913 
914  // statistics can be preserved only in case of positive coefficients
915  // otherwise with negative c1 (histogram subtraction) one risks to get negative variances
916  Bool_t resetStats = (c1 < 0);
917  Double_t s1[kNstat] = {0};
918  Double_t s2[kNstat] = {0};
919  if (!resetStats) {
920  // need to initialize to zero s1 and s2 since
921  // GetStats fills only used elements depending on dimension and type
922  GetStats(s1);
923  h1->GetStats(s2);
924  }
925 
926  SetMinimum();
927  SetMaximum();
928 
929  // - Loop on bins (including underflows/overflows)
930  Double_t factor = 1;
931  if (h1->GetNormFactor() != 0) factor = h1->GetNormFactor()/h1->GetSumOfWeights();;
932  Double_t c1sq = c1 * c1;
933  Double_t factsq = factor * factor;
934 
935  for (Int_t bin = 0; bin < fNcells; ++bin) {
936  //special case where histograms have the kIsAverage bit set
937  if (this->TestBit(kIsAverage) && h1->TestBit(kIsAverage)) {
938  Double_t y1 = h1->RetrieveBinContent(bin);
939  Double_t y2 = this->RetrieveBinContent(bin);
940  Double_t e1sq = h1->GetBinErrorSqUnchecked(bin);
941  Double_t e2sq = this->GetBinErrorSqUnchecked(bin);
942  Double_t w1 = 1., w2 = 1.;
943 
944  // consider all special cases when bin errors are zero
945  // see http://root-forum.cern.ch/viewtopic.php?f=3&t=13299
946  if (e1sq) w1 = 1. / e1sq;
947  else if (h1->fSumw2.fN) {
948  w1 = 1.E200; // use an arbitrary huge value
949  if (y1 == 0) {
950  // use an estimated error from the global histogram scale
951  double sf = (s2[0] != 0) ? s2[1]/s2[0] : 1;
952  w1 = 1./(sf*sf);
953  }
954  }
955  if (e2sq) w2 = 1. / e2sq;
956  else if (fSumw2.fN) {
957  w2 = 1.E200; // use an arbitrary huge value
958  if (y2 == 0) {
959  // use an estimated error from the global histogram scale
960  double sf = (s1[0] != 0) ? s1[1]/s1[0] : 1;
961  w2 = 1./(sf*sf);
962  }
963  }
964 
965  double y = (w1*y1 + w2*y2)/(w1 + w2);
966  UpdateBinContent(bin, y);
967  if (fSumw2.fN) {
968  double err2 = 1./(w1 + w2);
969  if (err2 < 1.E-200) err2 = 0; // to remove arbitrary value when e1=0 AND e2=0
970  fSumw2.fArray[bin] = err2;
971  }
972  } else { // normal case of addition between histograms
973  AddBinContent(bin, c1 * factor * h1->RetrieveBinContent(bin));
974  if (fSumw2.fN) fSumw2.fArray[bin] += c1sq * factsq * h1->GetBinErrorSqUnchecked(bin);
975  }
976  }
977 
978  // update statistics (do here to avoid changes by SetBinContent)
979  if (resetStats) {
980  // statistics need to be reset in case coefficient are negative
981  ResetStats();
982  }
983  else {
984  for (Int_t i=0;i<kNstat;i++) {
985  if (i == 1) s1[i] += c1*c1*s2[i];
986  else s1[i] += c1*s2[i];
987  }
988  PutStats(s1);
989  SetEntries(entries);
990  }
991  return kTRUE;
992 }
993 
994 ////////////////////////////////////////////////////////////////////////////////
995 /// Replace contents of this histogram by the addition of h1 and h2.
996 ///
997 /// `this = c1*h1 + c2*h2`
998 /// if errors are defined (see TH1::Sumw2), errors are also recalculated
999 ///
1000 /// Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this
1001 /// if not already set.
1002 ///
1003 /// Note also that adding histogram with labels is not supported, histogram will be
1004 /// added merging them by bin number independently of the labels.
1005 /// For adding histogram ith labels one should use TH1::Merge
1006 ///
1007 /// SPECIAL CASE (Average/Efficiency histograms)
1008 /// For histograms representing averages or efficiencies, one should compute the average
1009 /// of the two histograms and not the sum. One can mark a histogram to be an average
1010 /// histogram by setting its bit kIsAverage with
1011 /// myhist.SetBit(TH1::kIsAverage);
1012 /// Note that the two histograms must have their kIsAverage bit set
1013 ///
1014 /// IMPORTANT NOTE: If you intend to use the errors of this histogram later
1015 /// you should call Sumw2 before making this operation.
1016 /// This is particularly important if you fit the histogram after TH1::Add
1017 ///
1018 /// ANOTHER SPECIAL CASE : h1 = h2 and c2 < 0
1019 /// do a scaling this = c1 * h1 / (bin Volume)
1020 ///
1021 /// The function returns kFALSE if the Add operation failed
1022 
1024 {
1025 
1026  if (!h1 || !h2) {
1027  Error("Add","Attempt to add a non-existing histogram");
1028  return kFALSE;
1029  }
1030 
1031  // delete buffer if it is there since it will become invalid
1032  if (fBuffer) BufferEmpty(1);
1033 
1034  Bool_t normWidth = kFALSE;
1035  if (h1 == h2 && c2 < 0) {c2 = 0; normWidth = kTRUE;}
1036 
1037  if (h1 != h2) {
1038  bool useMerge = (c1 == 1. && c2 == 1. && !this->TestBit(kIsAverage) && !h1->TestBit(kIsAverage) );
1039 
1040  try {
1041  CheckConsistency(h1,h2);
1042  CheckConsistency(this,h1);
1043  useMerge = kFALSE;
1044  } catch(DifferentNumberOfBins&) {
1045  if (useMerge)
1046  Info("Add","Attempt to add histograms with different number of bins - trying to use TH1::Merge");
1047  else {
1048  Error("Add","Attempt to add histograms with different number of bins : nbins h1 = %d , nbins h2 = %d",GetNbinsX(), h1->GetNbinsX());
1049  return kFALSE;
1050  }
1051  } catch(DifferentAxisLimits&) {
1052  if (useMerge)
1053  Info("Add","Attempt to add histograms with different axis limits - trying to use TH1::Merge");
1054  else
1055  Warning("Add","Attempt to add histograms with different axis limits");
1056  } catch(DifferentBinLimits&) {
1057  if (useMerge)
1058  Info("Add","Attempt to add histograms with different bin limits - trying to use TH1::Merge");
1059  else
1060  Warning("Add","Attempt to add histograms with different bin limits");
1061  } catch(DifferentLabels&) {
1062  // in case of different labels -
1063  if (useMerge)
1064  Info("Add","Attempt to add histograms with different labels - trying to use TH1::Merge");
1065  else
1066  Info("Warning","Attempt to add histograms with different labels");
1067  }
1068 
1069  if (useMerge) {
1070  TList l;
1071  // why TList takes non-const pointers ????
1072  l.Add(const_cast<TH1*>(h1));
1073  l.Add(const_cast<TH1*>(h2));
1074  Reset("ICE");
1075  auto iret = Merge(&l);
1076  return (iret >= 0);
1077  }
1078  }
1079 
1080  // Create Sumw2 if h1 or h2 have Sumw2 set
1081  if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
1082 
1083  // - Add statistics
1084  Double_t nEntries = TMath::Abs( c1*h1->GetEntries() + c2*h2->GetEntries() );
1085 
1086  // TODO remove
1087  // statistics can be preserved only in case of positive coefficients
1088  // otherwise with negative c1 (histogram subtraction) one risks to get negative variances
1089  // also in case of scaling with the width we cannot preserve the statistics
1090  Double_t s1[kNstat] = {0};
1091  Double_t s2[kNstat] = {0};
1092  Double_t s3[kNstat];
1093 
1094 
1095  Bool_t resetStats = (c1*c2 < 0) || normWidth;
1096  if (!resetStats) {
1097  // need to initialize to zero s1 and s2 since
1098  // GetStats fills only used elements depending on dimension and type
1099  h1->GetStats(s1);
1100  h2->GetStats(s2);
1101  for (Int_t i=0;i<kNstat;i++) {
1102  if (i == 1) s3[i] = c1*c1*s1[i] + c2*c2*s2[i];
1103  //else s3[i] = TMath::Abs(c1)*s1[i] + TMath::Abs(c2)*s2[i];
1104  else s3[i] = c1*s1[i] + c2*s2[i];
1105  }
1106  }
1107 
1108  SetMinimum();
1109  SetMaximum();
1110 
1111  if (normWidth) { // DEPRECATED CASE: belongs to fitting / drawing modules
1112 
1113  Int_t nbinsx = GetNbinsX() + 2; // normal bins + underflow, overflow
1114  Int_t nbinsy = GetNbinsY() + 2;
1115  Int_t nbinsz = GetNbinsZ() + 2;
1116 
1117  if (fDimension < 2) nbinsy = 1;
1118  if (fDimension < 3) nbinsz = 1;
1119 
1120  Int_t bin, binx, biny, binz;
1121  for (binz = 0; binz < nbinsz; ++binz) {
1122  Double_t wz = h1->GetZaxis()->GetBinWidth(binz);
1123  for (biny = 0; biny < nbinsy; ++biny) {
1124  Double_t wy = h1->GetYaxis()->GetBinWidth(biny);
1125  for (binx = 0; binx < nbinsx; ++binx) {
1126  Double_t wx = h1->GetXaxis()->GetBinWidth(binx);
1127  bin = GetBin(binx, biny, binz);
1128  Double_t w = wx*wy*wz;
1129  UpdateBinContent(bin, c1 * h1->RetrieveBinContent(bin) / w);
1130  if (fSumw2.fN) {
1131  Double_t e1 = h1->GetBinError(bin)/w;
1132  fSumw2.fArray[bin] = c1*c1*e1*e1;
1133  }
1134  }
1135  }
1136  }
1137  } else if (h1->TestBit(kIsAverage) && h2->TestBit(kIsAverage)) {
1138  for (Int_t i = 0; i < fNcells; ++i) { // loop on cells (bins including underflow / overflow)
1139  // special case where histograms have the kIsAverage bit set
1140  Double_t y1 = h1->RetrieveBinContent(i);
1141  Double_t y2 = h2->RetrieveBinContent(i);
1142  Double_t e1sq = h1->GetBinErrorSqUnchecked(i);
1143  Double_t e2sq = h2->GetBinErrorSqUnchecked(i);
1144  Double_t w1 = 1., w2 = 1.;
1145 
1146  // consider all special cases when bin errors are zero
1147  // see http://root-forum.cern.ch/viewtopic.php?f=3&t=13299
1148  if (e1sq) w1 = 1./ e1sq;
1149  else if (h1->fSumw2.fN) {
1150  w1 = 1.E200; // use an arbitrary huge value
1151  if (y1 == 0 ) { // use an estimated error from the global histogram scale
1152  double sf = (s1[0] != 0) ? s1[1]/s1[0] : 1;
1153  w1 = 1./(sf*sf);
1154  }
1155  }
1156  if (e2sq) w2 = 1./ e2sq;
1157  else if (h2->fSumw2.fN) {
1158  w2 = 1.E200; // use an arbitrary huge value
1159  if (y2 == 0) { // use an estimated error from the global histogram scale
1160  double sf = (s2[0] != 0) ? s2[1]/s2[0] : 1;
1161  w2 = 1./(sf*sf);
1162  }
1163  }
1164 
1165  double y = (w1*y1 + w2*y2)/(w1 + w2);
1166  UpdateBinContent(i, y);
1167  if (fSumw2.fN) {
1168  double err2 = 1./(w1 + w2);
1169  if (err2 < 1.E-200) err2 = 0; // to remove arbitrary value when e1=0 AND e2=0
1170  fSumw2.fArray[i] = err2;
1171  }
1172  }
1173  } else { // case of simple histogram addition
1174  Double_t c1sq = c1 * c1;
1175  Double_t c2sq = c2 * c2;
1176  for (Int_t i = 0; i < fNcells; ++i) { // Loop on cells (bins including underflows/overflows)
1178  if (fSumw2.fN) {
1179  fSumw2.fArray[i] = c1sq * h1->GetBinErrorSqUnchecked(i) + c2sq * h2->GetBinErrorSqUnchecked(i);
1180  }
1181  }
1182  }
1183 
1184  if (resetStats) {
1185  // statistics need to be reset in case coefficient are negative
1186  ResetStats();
1187  }
1188  else {
1189  // update statistics (do here to avoid changes by SetBinContent) FIXME remove???
1190  PutStats(s3);
1191  SetEntries(nEntries);
1192  }
1193 
1194  return kTRUE;
1195 }
1196 
1197 ////////////////////////////////////////////////////////////////////////////////
1198 /// Increment bin content by 1.
1199 
1201 {
1202  AbstractMethod("AddBinContent");
1203 }
1204 
1205 ////////////////////////////////////////////////////////////////////////////////
1206 /// Increment bin content by a weight w.
1207 
1209 {
1210  AbstractMethod("AddBinContent");
1211 }
1212 
1213 ////////////////////////////////////////////////////////////////////////////////
1214 /// Sets the flag controlling the automatic add of histograms in memory
1215 ///
1216 /// By default (fAddDirectory = kTRUE), histograms are automatically added
1217 /// to the list of objects in memory.
1218 /// Note that one histogram can be removed from its support directory
1219 /// by calling h->SetDirectory(0) or h->SetDirectory(dir) to add it
1220 /// to the list of objects in the directory dir.
1221 ///
1222 /// NOTE that this is a static function. To call it, use;
1223 /// TH1::AddDirectory
1224 
1226 {
1227  fgAddDirectory = add;
1228 }
1229 
1230 ////////////////////////////////////////////////////////////////////////////////
1231 /// Auxilliary function to get the power of 2 next (larger) or previous (smaller)
1232 /// a given x
1233 ///
1234 /// next = kTRUE : next larger
1235 /// next = kFALSE : previous smaller
1236 ///
1237 /// Used by the autobin power of 2 algorithm
1238 
1240 {
1241  Int_t nn;
1242  Double_t f2 = std::frexp(x, &nn);
1243  return ((next && x > 0.) || (!next && x <= 0.)) ? std::ldexp(std::copysign(1., f2), nn)
1244  : std::ldexp(std::copysign(1., f2), --nn);
1245 }
1246 
1247 ////////////////////////////////////////////////////////////////////////////////
1248 /// Auxilliary function to get the next power of 2 integer value larger then n
1249 ///
1250 /// Used by the autobin power of 2 algorithm
1251 
1253 {
1254  Int_t nn;
1255  Double_t f2 = std::frexp(n, &nn);
1256  if (TMath::Abs(f2 - .5) > 0.001)
1257  return (Int_t)std::ldexp(1., nn);
1258  return n;
1259 }
1260 
1261 ////////////////////////////////////////////////////////////////////////////////
1262 /// Buffer-based estimate of the histogram range using the power of 2 algorithm.
1263 ///
1264 /// Used by the autobin power of 2 algorithm.
1265 ///
1266 /// Works on arguments (min and max from fBuffer) and internal inputs: fXmin,
1267 /// fXmax, NBinsX (from fXaxis), ...
1268 /// Result save internally in fXaxis.
1269 ///
1270 /// Overloaded by TH2 and TH3.
1271 ///
1272 /// Return -1 if internal inputs are incosistent, 0 otherwise.
1273 ///
1274 
1276 {
1277  // We need meaningful raw limits
1278  if (xmi >= xma)
1279  return -1;
1280 
1282  Double_t xhmi = fXaxis.GetXmin();
1283  Double_t xhma = fXaxis.GetXmax();
1284 
1285  // Now adjust
1286  if (TMath::Abs(xhma) > TMath::Abs(xhmi)) {
1287  // Start from the upper limit
1288  xhma = TH1::AutoP2GetPower2(xhma);
1289  xhmi = xhma - TH1::AutoP2GetPower2(xhma - xhmi);
1290  } else {
1291  // Start from the lower limit
1292  xhmi = TH1::AutoP2GetPower2(xhmi, kFALSE);
1293  xhma = xhmi + TH1::AutoP2GetPower2(xhma - xhmi);
1294  }
1295 
1296  // Round the bins to the next power of 2; take into account the possible inflation
1297  // of the range
1298  Double_t rr = (xhma - xhmi) / (xma - xmi);
1299  Int_t nb = TH1::AutoP2GetBins((Int_t)(rr * GetNbinsX()));
1300 
1301  // Adjust using the same bin width and offsets
1302  Double_t bw = (xhma - xhmi) / nb;
1303  // Bins to left free on each side
1304  Double_t autoside = gEnv->GetValue("Hist.Binning.Auto.Side", 0.05);
1305  Int_t nbside = (Int_t)(nb * autoside);
1306 
1307  // Side up
1308  Int_t nbup = (xhma - xma) / bw;
1309  if (nbup % 2 != 0)
1310  nbup++; // Must be even
1311  if (nbup != nbside) {
1312  // Accounts also for both case: larger or smaller
1313  xhma -= bw * (nbup - nbside);
1314  nb -= (nbup - nbside);
1315  }
1316 
1317  // Side low
1318  Int_t nblw = (xmi - xhmi) / bw;
1319  if (nblw % 2 != 0)
1320  nblw++; // Must be even
1321  if (nblw != nbside) {
1322  // Accounts also for both case: larger or smaller
1323  xhmi += bw * (nblw - nbside);
1324  nb -= (nblw - nbside);
1325  }
1326 
1327  // Set everything and project
1328  SetBins(nb, xhmi, xhma);
1329 
1330  // Done
1331  return 0;
1332 }
1333 
1334 /// Fill histogram with all entries in the buffer.
1335 ///
1336 /// - action = -1 histogram is reset and refilled from the buffer (called by THistPainter::Paint)
1337 /// - action = 0 histogram is reset and filled from the buffer. When the histogram is filled from the
1338 /// buffer the value fBuffer[0] is set to a negative number (= - number of entries)
1339 /// When calling with action == 0 the histogram is NOT refilled when fBuffer[0] is < 0
1340 /// While when calling with action = -1 the histogram is reset and ALWAYS refilled independently if
1341 /// the histogram was filled before. This is needed when drawing the histogram
1342 /// - action = 1 histogram is filled and buffer is deleted
1343 /// The buffer is automatically deleted when filling the histogram and the entries is
1344 /// larger than the buffer size
1345 
1347 {
1348  // do we need to compute the bin size?
1349  if (!fBuffer) return 0;
1350  Int_t nbentries = (Int_t)fBuffer[0];
1351 
1352  // nbentries correspond to the number of entries of histogram
1353 
1354  if (nbentries == 0) {
1355  // if action is 1 we delete the buffer
1356  // this will avoid infinite recursion
1357  if (action > 0) {
1358  delete [] fBuffer;
1359  fBuffer = 0;
1360  fBufferSize = 0;
1361  }
1362  return 0;
1363  }
1364  if (nbentries < 0 && action == 0) return 0; // case histogram has been already filled from the buffer
1365 
1366  Double_t *buffer = fBuffer;
1367  if (nbentries < 0) {
1368  nbentries = -nbentries;
1369  // a reset might call BufferEmpty() giving an infinite recursion
1370  // Protect it by setting fBuffer = 0
1371  fBuffer=0;
1372  //do not reset the list of functions
1373  Reset("ICES");
1374  fBuffer = buffer;
1375  }
1376  if (CanExtendAllAxes() || (fXaxis.GetXmax() <= fXaxis.GetXmin())) {
1377  //find min, max of entries in buffer
1378  Double_t xmin = fBuffer[2];
1379  Double_t xmax = xmin;
1380  for (Int_t i=1;i<nbentries;i++) {
1381  Double_t x = fBuffer[2*i+2];
1382  if (x < xmin) xmin = x;
1383  if (x > xmax) xmax = x;
1384  }
1385  if (fXaxis.GetXmax() <= fXaxis.GetXmin()) {
1386  Int_t rc = -1;
1387  if (TestBit(TH1::kAutoBinPTwo)) {
1388  if ((rc = AutoP2FindLimits(xmin, xmax)) < 0)
1389  Warning("BufferEmpty",
1390  "incosistency found by power-of-2 autobin algorithm: fallback to standard method");
1391  }
1392  if (rc < 0)
1394  } else {
1395  fBuffer = 0;
1396  Int_t keep = fBufferSize; fBufferSize = 0;
1397  if (xmin < fXaxis.GetXmin()) ExtendAxis(xmin, &fXaxis);
1398  if (xmax >= fXaxis.GetXmax()) ExtendAxis(xmax, &fXaxis);
1399  fBuffer = buffer;
1400  fBufferSize = keep;
1401  }
1402  }
1403 
1404  // call DoFillN which will not put entries in the buffer as FillN does
1405  // set fBuffer to zero to avoid re-emptying the buffer from functions called
1406  // by DoFillN (e.g Sumw2)
1407  buffer = fBuffer; fBuffer = 0;
1408  DoFillN(nbentries,&buffer[2],&buffer[1],2);
1409  fBuffer = buffer;
1410 
1411  // if action == 1 - delete the buffer
1412  if (action > 0) {
1413  delete [] fBuffer;
1414  fBuffer = 0;
1415  fBufferSize = 0;
1416  } else {
1417  // if number of entries is consistent with buffer - set it negative to avoid
1418  // refilling the histogram every time BufferEmpty(0) is called
1419  // In case it is not consistent, by setting fBuffer[0]=0 is like resetting the buffer
1420  // (it will not be used anymore the next time BufferEmpty is called)
1421  if (nbentries == (Int_t)fEntries)
1422  fBuffer[0] = -nbentries;
1423  else
1424  fBuffer[0] = 0;
1425  }
1426  return nbentries;
1427 }
1428 
1429 ////////////////////////////////////////////////////////////////////////////////
1430 /// accumulate arguments in buffer. When buffer is full, empty the buffer
1431 ///
1432 /// - `fBuffer[0]` = number of entries in buffer
1433 /// - `fBuffer[1]` = w of first entry
1434 /// - `fBuffer[2]` = x of first entry
1435 
1437 {
1438  if (!fBuffer) return -2;
1439  Int_t nbentries = (Int_t)fBuffer[0];
1440 
1441 
1442  if (nbentries < 0) {
1443  // reset nbentries to a positive value so next time BufferEmpty() is called
1444  // the histogram will be refilled
1445  nbentries = -nbentries;
1446  fBuffer[0] = nbentries;
1447  if (fEntries > 0) {
1448  // set fBuffer to zero to avoid calling BufferEmpty in Reset
1449  Double_t *buffer = fBuffer; fBuffer=0;
1450  Reset("ICES"); // do not reset list of functions
1451  fBuffer = buffer;
1452  }
1453  }
1454  if (2*nbentries+2 >= fBufferSize) {
1455  BufferEmpty(1);
1456  if (!fBuffer)
1457  // to avoid infinite recursion Fill->BufferFill->Fill
1458  return Fill(x,w);
1459  // this cannot happen
1460  R__ASSERT(0);
1461  }
1462  fBuffer[2*nbentries+1] = w;
1463  fBuffer[2*nbentries+2] = x;
1464  fBuffer[0] += 1;
1465  return -2;
1466 }
1467 
1468 ////////////////////////////////////////////////////////////////////////////////
1469 /// Check bin limits.
1470 
1471 bool TH1::CheckBinLimits(const TAxis* a1, const TAxis * a2)
1472 {
1473  const TArrayD * h1Array = a1->GetXbins();
1474  const TArrayD * h2Array = a2->GetXbins();
1475  Int_t fN = h1Array->fN;
1476  if ( fN != 0 ) {
1477  if ( h2Array->fN != fN ) {
1478  throw DifferentBinLimits();
1479  return false;
1480  }
1481  else {
1482  for ( int i = 0; i < fN; ++i ) {
1483  // for i==fN (nbin+1) a->GetBinWidth() returns last bin width
1484  // we do not need to exclude that case
1485  double binWidth = a1->GetBinWidth(i);
1486  if ( ! TMath::AreEqualAbs( h1Array->GetAt(i), h2Array->GetAt(i), binWidth*1E-10 ) ) {
1487  throw DifferentBinLimits();
1488  return false;
1489  }
1490  }
1491  }
1492  }
1493 
1494  return true;
1495 }
1496 
1497 ////////////////////////////////////////////////////////////////////////////////
1498 /// Check that axis have same labels.
1499 
1500 bool TH1::CheckBinLabels(const TAxis* a1, const TAxis * a2)
1501 {
1502  THashList *l1 = a1->GetLabels();
1503  THashList *l2 = a2->GetLabels();
1504 
1505  if (!l1 && !l2 )
1506  return true;
1507  if (!l1 || !l2 ) {
1508  throw DifferentLabels();
1509  return false;
1510  }
1511  // check now labels sizes are the same
1512  if (l1->GetSize() != l2->GetSize() ) {
1513  throw DifferentLabels();
1514  return false;
1515  }
1516  for (int i = 1; i <= a1->GetNbins(); ++i) {
1517  TString label1 = a1->GetBinLabel(i);
1518  TString label2 = a2->GetBinLabel(i);
1519  if (label1 != label2) {
1520  throw DifferentLabels();
1521  return false;
1522  }
1523  }
1524 
1525  return true;
1526 }
1527 
1528 ////////////////////////////////////////////////////////////////////////////////
1529 /// Check that the axis limits of the histograms are the same.
1530 /// If a first and last bin is passed the axis is compared between the given range
1531 
1532 bool TH1::CheckAxisLimits(const TAxis *a1, const TAxis *a2 )
1533 {
1534  double firstBin = a1->GetBinWidth(1);
1535  double lastBin = a1->GetBinWidth( a1->GetNbins() );
1536  if ( ! TMath::AreEqualAbs(a1->GetXmin(), a2->GetXmin(), firstBin* 1.E-10) ||
1537  ! TMath::AreEqualAbs(a1->GetXmax(), a2->GetXmax(), lastBin*1.E-10) ) {
1538  throw DifferentAxisLimits();
1539  return false;
1540  }
1541  return true;
1542 }
1543 
1544 ////////////////////////////////////////////////////////////////////////////////
1545 /// Check that the axis are the same
1546 
1547 bool TH1::CheckEqualAxes(const TAxis *a1, const TAxis *a2 )
1548 {
1549  if (a1->GetNbins() != a2->GetNbins() ) {
1550  //throw DifferentNumberOfBins();
1551  ::Info("CheckEqualAxes","Axes have different number of bins : nbin1 = %d nbin2 = %d",a1->GetNbins(),a2->GetNbins() );
1552  return false;
1553  }
1554  try {
1555  CheckAxisLimits(a1,a2);
1556  } catch (DifferentAxisLimits&) {
1557  ::Info("CheckEqualAxes","Axes have different limits");
1558  return false;
1559  }
1560  try {
1561  CheckBinLimits(a1,a2);
1562  } catch (DifferentBinLimits&) {
1563  ::Info("CheckEqualAxes","Axes have different bin limits");
1564  return false;
1565  }
1566 
1567  // check labels
1568  try {
1569  CheckBinLabels(a1,a2);
1570  } catch (DifferentLabels&) {
1571  ::Info("CheckEqualAxes","Axes have different labels");
1572  return false;
1573  }
1574 
1575  return true;
1576 }
1577 
1578 ////////////////////////////////////////////////////////////////////////////////
1579 /// Check that two sub axis are the same.
1580 /// The limits are defined by first bin and last bin
1581 /// N.B. no check is done in this case for variable bins
1582 
1583 bool TH1::CheckConsistentSubAxes(const TAxis *a1, Int_t firstBin1, Int_t lastBin1, const TAxis * a2, Int_t firstBin2, Int_t lastBin2 )
1584 {
1585  // By default is assumed that no bins are given for the second axis
1586  Int_t nbins1 = lastBin1-firstBin1 + 1;
1587  Double_t xmin1 = a1->GetBinLowEdge(firstBin1);
1588  Double_t xmax1 = a1->GetBinUpEdge(lastBin1);
1589 
1590  Int_t nbins2 = a2->GetNbins();
1591  Double_t xmin2 = a2->GetXmin();
1592  Double_t xmax2 = a2->GetXmax();
1593 
1594  if (firstBin2 < lastBin2) {
1595  // in this case assume no bins are given for the second axis
1596  nbins2 = lastBin1-firstBin1 + 1;
1597  xmin2 = a1->GetBinLowEdge(firstBin1);
1598  xmax2 = a1->GetBinUpEdge(lastBin1);
1599  }
1600 
1601  if (nbins1 != nbins2 ) {
1602  ::Info("CheckConsistentSubAxes","Axes have different number of bins");
1603  return false;
1604  }
1605 
1606  Double_t firstBin = a1->GetBinWidth(firstBin1);
1607  Double_t lastBin = a1->GetBinWidth(lastBin1);
1608  if ( ! TMath::AreEqualAbs(xmin1,xmin2,1.E-10 * firstBin) ||
1609  ! TMath::AreEqualAbs(xmax1,xmax2,1.E-10 * lastBin) ) {
1610  ::Info("CheckConsistentSubAxes","Axes have different limits");
1611  return false;
1612  }
1613 
1614  return true;
1615 }
1616 
1617 ////////////////////////////////////////////////////////////////////////////////
1618 /// Check histogram compatibility.
1619 
1620 bool TH1::CheckConsistency(const TH1* h1, const TH1* h2)
1621 {
1622  if (h1 == h2) return true;
1623 
1624  if (h1->GetDimension() != h2->GetDimension() ) {
1625  throw DifferentDimension();
1626  return false;
1627  }
1628  Int_t dim = h1->GetDimension();
1629 
1630  // returns kTRUE if number of bins and bin limits are identical
1631  Int_t nbinsx = h1->GetNbinsX();
1632  Int_t nbinsy = h1->GetNbinsY();
1633  Int_t nbinsz = h1->GetNbinsZ();
1634 
1635  // Check whether the histograms have the same number of bins.
1636  if (nbinsx != h2->GetNbinsX() ||
1637  (dim > 1 && nbinsy != h2->GetNbinsY()) ||
1638  (dim > 2 && nbinsz != h2->GetNbinsZ()) ) {
1639  throw DifferentNumberOfBins();
1640  return false;
1641  }
1642 
1643  bool ret = true;
1644 
1645  // check axis limits
1646  ret &= CheckAxisLimits(h1->GetXaxis(), h2->GetXaxis());
1647  if (dim > 1) ret &= CheckAxisLimits(h1->GetYaxis(), h2->GetYaxis());
1648  if (dim > 2) ret &= CheckAxisLimits(h1->GetZaxis(), h2->GetZaxis());
1649 
1650  // check bin limits
1651  ret &= CheckBinLimits(h1->GetXaxis(), h2->GetXaxis());
1652  if (dim > 1) ret &= CheckBinLimits(h1->GetYaxis(), h2->GetYaxis());
1653  if (dim > 2) ret &= CheckBinLimits(h1->GetZaxis(), h2->GetZaxis());
1654 
1655  // check labels if histograms are both not empty
1656  if ( !h1->IsEmpty() && !h2->IsEmpty() ) {
1657  ret &= CheckBinLabels(h1->GetXaxis(), h2->GetXaxis());
1658  if (dim > 1) ret &= CheckBinLabels(h1->GetYaxis(), h2->GetYaxis());
1659  if (dim > 2) ret &= CheckBinLabels(h1->GetZaxis(), h2->GetZaxis());
1660  }
1661 
1662  return ret;
1663 }
1664 
1665 ////////////////////////////////////////////////////////////////////////////////
1666 /// \f$ \chi^{2} \f$ test for comparing weighted and unweighted histograms
1667 ///
1668 /// Function: Returns p-value. Other return values are specified by the 3rd parameter
1669 ///
1670 /// \param[in] h2 the second histogram
1671 /// \param[in] option
1672 /// - "UU" = experiment experiment comparison (unweighted-unweighted)
1673 /// - "UW" = experiment MC comparison (unweighted-weighted). Note that
1674 /// the first histogram should be unweighted
1675 /// - "WW" = MC MC comparison (weighted-weighted)
1676 /// - "NORM" = to be used when one or both of the histograms is scaled
1677 /// but the histogram originally was unweighted
1678 /// - by default underflows and overflows are not included:
1679 /// * "OF" = overflows included
1680 /// * "UF" = underflows included
1681 /// - "P" = print chi2, ndf, p_value, igood
1682 /// - "CHI2" = returns chi2 instead of p-value
1683 /// - "CHI2/NDF" = returns \f$ \chi^{2} \f$/ndf
1684 /// \param[in] res not empty - computes normalized residuals and returns them in this array
1685 ///
1686 /// The current implementation is based on the papers \f$ \chi^{2} \f$ test for comparison
1687 /// of weighted and unweighted histograms" in Proceedings of PHYSTAT05 and
1688 /// "Comparison weighted and unweighted histograms", arXiv:physics/0605123
1689 /// by N.Gagunashvili. This function has been implemented by Daniel Haertl in August 2006.
1690 ///
1691 /// #### Introduction:
1692 ///
1693 /// A frequently used technique in data analysis is the comparison of
1694 /// histograms. First suggested by Pearson [1] the \f$ \chi^{2} \f$ test of
1695 /// homogeneity is used widely for comparing usual (unweighted) histograms.
1696 /// This paper describes the implementation modified \f$ \chi^{2} \f$ tests
1697 /// for comparison of weighted and unweighted histograms and two weighted
1698 /// histograms [2] as well as usual Pearson's \f$ \chi^{2} \f$ test for
1699 /// comparison two usual (unweighted) histograms.
1700 ///
1701 /// #### Overview:
1702 ///
1703 /// Comparison of two histograms expect hypotheses that two histograms
1704 /// represent identical distributions. To make a decision p-value should
1705 /// be calculated. The hypotheses of identity is rejected if the p-value is
1706 /// lower then some significance level. Traditionally significance levels
1707 /// 0.1, 0.05 and 0.01 are used. The comparison procedure should include an
1708 /// analysis of the residuals which is often helpful in identifying the
1709 /// bins of histograms responsible for a significant overall \f$ \chi^{2} \f$ value.
1710 /// Residuals are the difference between bin contents and expected bin
1711 /// contents. Most convenient for analysis are the normalized residuals. If
1712 /// hypotheses of identity are valid then normalized residuals are
1713 /// approximately independent and identically distributed random variables
1714 /// having N(0,1) distribution. Analysis of residuals expect test of above
1715 /// mentioned properties of residuals. Notice that indirectly the analysis
1716 /// of residuals increase the power of \f$ \chi^{2} \f$ test.
1717 ///
1718 /// #### Methods of comparison:
1719 ///
1720 /// \f$ \chi^{2} \f$ test for comparison two (unweighted) histograms:
1721 /// Let us consider two histograms with the same binning and the number
1722 /// of bins equal to r. Let us denote the number of events in the ith bin
1723 /// in the first histogram as ni and as mi in the second one. The total
1724 /// number of events in the first histogram is equal to:
1725 /// \f[
1726 /// N = \sum_{i=1}^{r} n_{i}
1727 /// \f]
1728 /// and
1729 /// \f[
1730 /// M = \sum_{i=1}^{r} m_{i}
1731 /// \f]
1732 /// in the second histogram. The hypothesis of identity (homogeneity) [3]
1733 /// is that the two histograms represent random values with identical
1734 /// distributions. It is equivalent that there exist r constants p1,...,pr,
1735 /// such that
1736 /// \f[
1737 ///\sum_{i=1}^{r} p_{i}=1
1738 /// \f]
1739 /// and the probability of belonging to the ith bin for some measured value
1740 /// in both experiments is equal to pi. The number of events in the ith
1741 /// bin is a random variable with a distribution approximated by a Poisson
1742 /// probability distribution
1743 /// \f[
1744 ///\frac{e^{-Np_{i}}(Np_{i})^{n_{i}}}{n_{i}!}
1745 /// \f]
1746 ///for the first histogram and with distribution
1747 /// \f[
1748 ///\frac{e^{-Mp_{i}}(Mp_{i})^{m_{i}}}{m_{i}!}
1749 /// \f]
1750 /// for the second histogram. If the hypothesis of homogeneity is valid,
1751 /// then the maximum likelihood estimator of pi, i=1,...,r, is
1752 /// \f[
1753 ///\hat{p}_{i}= \frac{n_{i}+m_{i}}{N+M}
1754 /// \f]
1755 /// and then
1756 /// \f[
1757 /// X^{2} = \sum_{i=1}^{r}\frac{(n_{i}-N\hat{p}_{i})^{2}}{N\hat{p}_{i}} + \sum_{i=1}^{r}\frac{(m_{i}-M\hat{p}_{i})^{2}}{M\hat{p}_{i}} =\frac{1}{MN} \sum_{i=1}^{r}\frac{(Mn_{i}-Nm_{i})^{2}}{n_{i}+m_{i}}
1758 /// \f]
1759 /// has approximately a \f$ \chi^{2}_{(r-1)} \f$ distribution [3].
1760 /// The comparison procedure can include an analysis of the residuals which
1761 /// is often helpful in identifying the bins of histograms responsible for
1762 /// a significant overall \f$ \chi^{2} \f$ value. Most convenient for
1763 /// analysis are the adjusted (normalized) residuals [4]
1764 /// \f[
1765 /// r_{i} = \frac{n_{i}-N\hat{p}_{i}}{\sqrt{N\hat{p}_{i}}\sqrt{(1-N/(N+M))(1-(n_{i}+m_{i})/(N+M))}}
1766 /// \f]
1767 /// If hypotheses of homogeneity are valid then residuals ri are
1768 /// approximately independent and identically distributed random variables
1769 /// having N(0,1) distribution. The application of the \f$ \chi^{2} \f$ test has
1770 /// restrictions related to the value of the expected frequencies Npi,
1771 /// Mpi, i=1,...,r. A conservative rule formulated in [5] is that all the
1772 /// expectations must be 1 or greater for both histograms. In practical
1773 /// cases when expected frequencies are not known the estimated expected
1774 /// frequencies \f$ M\hat{p}_{i}, N\hat{p}_{i}, i=1,...,r \f$ can be used.
1775 ///
1776 /// #### Unweighted and weighted histograms comparison:
1777 ///
1778 /// A simple modification of the ideas described above can be used for the
1779 /// comparison of the usual (unweighted) and weighted histograms. Let us
1780 /// denote the number of events in the ith bin in the unweighted
1781 /// histogram as ni and the common weight of events in the ith bin of the
1782 /// weighted histogram as wi. The total number of events in the
1783 /// unweighted histogram is equal to
1784 ///\f[
1785 /// N = \sum_{i=1}^{r} n_{i}
1786 ///\f]
1787 /// and the total weight of events in the weighted histogram is equal to
1788 ///\f[
1789 /// W = \sum_{i=1}^{r} w_{i}
1790 ///\f]
1791 /// Let us formulate the hypothesis of identity of an unweighted histogram
1792 /// to a weighted histogram so that there exist r constants p1,...,pr, such
1793 /// that
1794 ///\f[
1795 /// \sum_{i=1}^{r} p_{i} = 1
1796 ///\f]
1797 /// for the unweighted histogram. The weight wi is a random variable with a
1798 /// distribution approximated by the normal probability distribution
1799 /// \f$ N(Wp_{i},\sigma_{i}^{2}) \f$ where \f$ \sigma_{i}^{2} \f$ is the variance of the weight wi.
1800 /// If we replace the variance \f$ \sigma_{i}^{2} \f$
1801 /// with estimate \f$ s_{i}^{2} \f$ (sum of squares of weights of
1802 /// events in the ith bin) and the hypothesis of identity is valid, then the
1803 /// maximum likelihood estimator of pi,i=1,...,r, is
1804 ///\f[
1805 /// \hat{p}_{i} = \frac{Ww_{i}-Ns_{i}^{2}+\sqrt{(Ww_{i}-Ns_{i}^{2})^{2}+4W^{2}s_{i}^{2}n_{i}}}{2W^{2}}
1806 ///\f]
1807 /// We may then use the test statistic
1808 ///\f[
1809 /// X^{2} = \sum_{i=1}^{r} \frac{(n_{i}-N\hat{p}_{i})^{2}}{N\hat{p}_{i}} + \sum_{i=1}^{r} \frac{(w_{i}-W\hat{p}_{i})^{2}}{s_{i}^{2}}
1810 ///\f]
1811 /// and it has approximately a \f$ \sigma^{2}_{(r-1)} \f$ distribution [2]. This test, as well
1812 /// as the original one [3], has a restriction on the expected frequencies. The
1813 /// expected frequencies recommended for the weighted histogram is more than 25.
1814 /// The value of the minimal expected frequency can be decreased down to 10 for
1815 /// the case when the weights of the events are close to constant. In the case
1816 /// of a weighted histogram if the number of events is unknown, then we can
1817 /// apply this recommendation for the equivalent number of events as
1818 ///\f[
1819 /// n_{i}^{equiv} = \frac{ w_{i}^{2} }{ s_{i}^{2} }
1820 ///\f]
1821 /// The minimal expected frequency for an unweighted histogram must be 1. Notice
1822 /// that any usual (unweighted) histogram can be considered as a weighted
1823 /// histogram with events that have constant weights equal to 1.
1824 /// The variance \f$ z_{i}^{2} \f$ of the difference between the weight wi
1825 /// and the estimated expectation value of the weight is approximately equal to:
1826 ///\f[
1827 /// z_{i}^{2} = Var(w_{i}-W\hat{p}_{i}) = N\hat{p}_{i}(1-N\hat{p}_{i})\left(\frac{Ws_{i}^{2}}{\sqrt{(Ns_{i}^{2}-w_{i}W)^{2}+4W^{2}s_{i}^{2}n_{i}}}\right)^{2}+\frac{s_{i}^{2}}{4}\left(1+\frac{Ns_{i}^{2}-w_{i}W}{\sqrt{(Ns_{i}^{2}-w_{i}W)^{2}+4W^{2}s_{i}^{2}n_{i}}}\right)^{2}
1828 ///\f]
1829 /// The residuals
1830 ///\f[
1831 /// r_{i} = \frac{w_{i}-W\hat{p}_{i}}{z_{i}}
1832 ///\f]
1833 /// have approximately a normal distribution with mean equal to 0 and standard
1834 /// deviation equal to 1.
1835 ///
1836 /// #### Two weighted histograms comparison:
1837 ///
1838 /// Let us denote the common weight of events of the ith bin in the first
1839 /// histogram as w1i and as w2i in the second one. The total weight of events
1840 /// in the first histogram is equal to
1841 ///\f[
1842 /// W_{1} = \sum_{i=1}^{r} w_{1i}
1843 ///\f]
1844 /// and
1845 ///\f[
1846 /// W_{2} = \sum_{i=1}^{r} w_{2i}
1847 ///\f]
1848 /// in the second histogram. Let us formulate the hypothesis of identity of
1849 /// weighted histograms so that there exist r constants p1,...,pr, such that
1850 ///\f[
1851 /// \sum_{i=1}^{r} p_{i} = 1
1852 ///\f]
1853 /// and also expectation value of weight w1i equal to W1pi and expectation value
1854 /// of weight w2i equal to W2pi. Weights in both the histograms are random
1855 /// variables with distributions which can be approximated by a normal
1856 /// probability distribution \f$ N(W_{1}p_{i},\sigma_{1i}^{2}) \f$ for the first histogram
1857 /// and by a distribution \f$ N(W_{2}p_{i},\sigma_{2i}^{2}) \f$ for the second.
1858 /// Here \f$ \sigma_{1i}^{2} \f$ and \f$ \sigma_{2i}^{2} \f$ are the variances
1859 /// of w1i and w2i with estimators \f$ s_{1i}^{2} \f$ and \f$ s_{2i}^{2} \f$ respectively.
1860 /// If the hypothesis of identity is valid, then the maximum likelihood and
1861 /// Least Square Method estimator of pi,i=1,...,r, is
1862 ///\f[
1863 /// \hat{p}_{i} = \frac{w_{1i}W_{1}/s_{1i}^{2}+w_{2i}W_{2} /s_{2i}^{2}}{W_{1}^{2}/s_{1i}^{2}+W_{2}^{2}/s_{2i}^{2}}
1864 ///\f]
1865 /// We may then use the test statistic
1866 ///\f[
1867 /// X^{2} = \sum_{i=1}^{r} \frac{(w_{1i}-W_{1}\hat{p}_{i})^{2}}{s_{1i}^{2}} + \sum_{i=1}^{r} \frac{(w_{2i}-W_{2}\hat{p}_{i})^{2}}{s_{2i}^{2}} = \sum_{i=1}^{r} \frac{(W_{1}w_{2i}-W_{2}w_{1i})^{2}}{W_{1}^{2}s_{2i}^{2}+W_{2}^{2}s_{1i}^{2}}
1868 ///\f]
1869 /// and it has approximately a \f$ \chi^{2}_{(r-1)} \f$ distribution [2].
1870 /// The normalized or studentised residuals [6]
1871 ///\f[
1872 /// r_{i} = \frac{w_{1i}-W_{1}\hat{p}_{i}}{s_{1i}\sqrt{1 - \frac{1}{(1+W_{2}^{2}s_{1i}^{2}/W_{1}^{2}s_{2i}^{2})}}}
1873 ///\f]
1874 /// have approximately a normal distribution with mean equal to 0 and standard
1875 /// deviation 1. A recommended minimal expected frequency is equal to 10 for
1876 /// the proposed test.
1877 ///
1878 /// #### Numerical examples:
1879 ///
1880 /// The method described herein is now illustrated with an example.
1881 /// We take a distribution
1882 ///\f[
1883 /// \phi(x) = \frac{2}{(x-10)^{2}+1} + \frac{1}{(x-14)^{2}+1} (1)
1884 ///\f]
1885 /// defined on the interval [4,16]. Events distributed according to the formula
1886 /// (1) are simulated to create the unweighted histogram. Uniformly distributed
1887 /// events are simulated for the weighted histogram with weights calculated by
1888 /// formula (1). Each histogram has the same number of bins: 20. Fig.1 shows
1889 /// the result of comparison of the unweighted histogram with 200 events
1890 /// (minimal expected frequency equal to one) and the weighted histogram with
1891 /// 500 events (minimal expected frequency equal to 25)
1892 /// Begin_Macro
1893 /// ../../../tutorials/math/chi2test.C
1894 /// End_Macro
1895 /// Fig 1. An example of comparison of the unweighted histogram with 200 events
1896 /// and the weighted histogram with 500 events:
1897 /// 1. unweighted histogram;
1898 /// 2. weighted histogram;
1899 /// 3. normalized residuals plot;
1900 /// 4. normal Q-Q plot of residuals.
1901 ///
1902 /// The value of the test statistic \f$ \chi^{2} \f$ is equal to
1903 /// 21.09 with p-value equal to 0.33, therefore the hypothesis of identity of
1904 /// the two histograms can be accepted for 0.05 significant level. The behavior
1905 /// of the normalized residuals plot (see Fig. 1c) and the normal Q-Q plot
1906 /// (see Fig. 1d) of residuals are regular and we cannot identify the outliers
1907 /// or bins with a big influence on \f$ \chi^{2} \f$.
1908 ///
1909 /// The second example presents the same two histograms but 17 events was added
1910 /// to content of bin number 15 in unweighted histogram. Fig.2 shows the result
1911 /// of comparison of the unweighted histogram with 217 events (minimal expected
1912 /// frequency equal to one) and the weighted histogram with 500 events (minimal
1913 /// expected frequency equal to 25)
1914 /// Begin_Macro
1915 /// ../../../tutorials/math/chi2test.C(17)
1916 /// End_Macro
1917 /// Fig 2. An example of comparison of the unweighted histogram with 217 events
1918 /// and the weighted histogram with 500 events:
1919 /// 1. unweighted histogram;
1920 /// 2. weighted histogram;
1921 /// 3. normalized residuals plot;
1922 /// 4. normal Q-Q plot of residuals.
1923 ///
1924 /// The value of the test statistic \f$ \chi^{2} \f$ is equal to
1925 /// 32.33 with p-value equal to 0.029, therefore the hypothesis of identity of
1926 /// the two histograms is rejected for 0.05 significant level. The behavior of
1927 /// the normalized residuals plot (see Fig. 2c) and the normal Q-Q plot (see
1928 /// Fig. 2d) of residuals are not regular and we can identify the outlier or
1929 /// bin with a big influence on \f$ \chi^{2} \f$.
1930 ///
1931 /// #### References:
1932 ///
1933 /// - [1] Pearson, K., 1904. On the Theory of Contingency and Its Relation to
1934 /// Association and Normal Correlation. Drapers' Co. Memoirs, Biometric
1935 /// Series No. 1, London.
1936 /// - [2] Gagunashvili, N., 2006. \f$ \sigma^{2} \f$ test for comparison
1937 /// of weighted and unweighted histograms. Statistical Problems in Particle
1938 /// Physics, Astrophysics and Cosmology, Proceedings of PHYSTAT05,
1939 /// Oxford, UK, 12-15 September 2005, Imperial College Press, London, 43-44.
1940 /// Gagunashvili,N., Comparison of weighted and unweighted histograms,
1941 /// arXiv:physics/0605123, 2006.
1942 /// - [3] Cramer, H., 1946. Mathematical methods of statistics.
1943 /// Princeton University Press, Princeton.
1944 /// - [4] Haberman, S.J., 1973. The analysis of residuals in cross-classified tables.
1945 /// Biometrics 29, 205-220.
1946 /// - [5] Lewontin, R.C. and Felsenstein, J., 1965. The robustness of homogeneity
1947 /// test in 2xN tables. Biometrics 21, 19-33.
1948 /// - [6] Seber, G.A.F., Lee, A.J., 2003, Linear Regression Analysis.
1949 /// John Wiley & Sons Inc., New York.
1950 
1951 Double_t TH1::Chi2Test(const TH1* h2, Option_t *option, Double_t *res) const
1952 {
1953  Double_t chi2 = 0;
1954  Int_t ndf = 0, igood = 0;
1955 
1956  TString opt = option;
1957  opt.ToUpper();
1958 
1959  Double_t prob = Chi2TestX(h2,chi2,ndf,igood,option,res);
1960 
1961  if(opt.Contains("P")) {
1962  printf("Chi2 = %f, Prob = %g, NDF = %d, igood = %d\n", chi2,prob,ndf,igood);
1963  }
1964  if(opt.Contains("CHI2/NDF")) {
1965  if (ndf == 0) return 0;
1966  return chi2/ndf;
1967  }
1968  if(opt.Contains("CHI2")) {
1969  return chi2;
1970  }
1971 
1972  return prob;
1973 }
1974 
1975 ////////////////////////////////////////////////////////////////////////////////
1976 /// The computation routine of the Chisquare test. For the method description,
1977 /// see Chi2Test() function.
1978 ///
1979 /// \return p-value
1980 /// \param[in] h2 the second histogram
1981 /// \param[in] option
1982 /// - "UU" = experiment experiment comparison (unweighted-unweighted)
1983 /// - "UW" = experiment MC comparison (unweighted-weighted). Note that the first
1984 /// histogram should be unweighted
1985 /// - "WW" = MC MC comparison (weighted-weighted)
1986 /// - "NORM" = if one or both histograms is scaled
1987 /// - "OF" = overflows included
1988 /// - "UF" = underflows included
1989 /// by default underflows and overflows are not included
1990 /// \param[out] igood test output
1991 /// - igood=0 - no problems
1992 /// - For unweighted unweighted comparison
1993 /// - igood=1'There is a bin in the 1st histogram with less than 1 event'
1994 /// - igood=2'There is a bin in the 2nd histogram with less than 1 event'
1995 /// - igood=3'when the conditions for igood=1 and igood=2 are satisfied'
1996 /// - For unweighted weighted comparison
1997 /// - igood=1'There is a bin in the 1st histogram with less then 1 event'
1998 /// - igood=2'There is a bin in the 2nd histogram with less then 10 effective number of events'
1999 /// - igood=3'when the conditions for igood=1 and igood=2 are satisfied'
2000 /// - For weighted weighted comparison
2001 /// - igood=1'There is a bin in the 1st histogram with less then 10 effective
2002 /// number of events'
2003 /// - igood=2'There is a bin in the 2nd histogram with less then 10 effective
2004 /// number of events'
2005 /// - igood=3'when the conditions for igood=1 and igood=2 are satisfied'
2006 /// \param[out] chi2 chisquare of the test
2007 /// \param[out] ndf number of degrees of freedom (important, when both histograms have the same empty bins)
2008 /// \param[out] res normalized residuals for further analysis
2009 
2010 Double_t TH1::Chi2TestX(const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood, Option_t *option, Double_t *res) const
2011 {
2012 
2013  Int_t i_start, i_end;
2014  Int_t j_start, j_end;
2015  Int_t k_start, k_end;
2016 
2017  Double_t sum1 = 0.0, sumw1 = 0.0;
2018  Double_t sum2 = 0.0, sumw2 = 0.0;
2019 
2020  chi2 = 0.0;
2021  ndf = 0;
2022 
2023  TString opt = option;
2024  opt.ToUpper();
2025 
2026  if (fBuffer) const_cast<TH1*>(this)->BufferEmpty();
2027 
2028  const TAxis *xaxis1 = GetXaxis();
2029  const TAxis *xaxis2 = h2->GetXaxis();
2030  const TAxis *yaxis1 = GetYaxis();
2031  const TAxis *yaxis2 = h2->GetYaxis();
2032  const TAxis *zaxis1 = GetZaxis();
2033  const TAxis *zaxis2 = h2->GetZaxis();
2034 
2035  Int_t nbinx1 = xaxis1->GetNbins();
2036  Int_t nbinx2 = xaxis2->GetNbins();
2037  Int_t nbiny1 = yaxis1->GetNbins();
2038  Int_t nbiny2 = yaxis2->GetNbins();
2039  Int_t nbinz1 = zaxis1->GetNbins();
2040  Int_t nbinz2 = zaxis2->GetNbins();
2041 
2042  //check dimensions
2043  if (this->GetDimension() != h2->GetDimension() ){
2044  Error("Chi2TestX","Histograms have different dimensions.");
2045  return 0.0;
2046  }
2047 
2048  //check number of channels
2049  if (nbinx1 != nbinx2) {
2050  Error("Chi2TestX","different number of x channels");
2051  }
2052  if (nbiny1 != nbiny2) {
2053  Error("Chi2TestX","different number of y channels");
2054  }
2055  if (nbinz1 != nbinz2) {
2056  Error("Chi2TestX","different number of z channels");
2057  }
2058 
2059  //check for ranges
2060  i_start = j_start = k_start = 1;
2061  i_end = nbinx1;
2062  j_end = nbiny1;
2063  k_end = nbinz1;
2064 
2065  if (xaxis1->TestBit(TAxis::kAxisRange)) {
2066  i_start = xaxis1->GetFirst();
2067  i_end = xaxis1->GetLast();
2068  }
2069  if (yaxis1->TestBit(TAxis::kAxisRange)) {
2070  j_start = yaxis1->GetFirst();
2071  j_end = yaxis1->GetLast();
2072  }
2073  if (zaxis1->TestBit(TAxis::kAxisRange)) {
2074  k_start = zaxis1->GetFirst();
2075  k_end = zaxis1->GetLast();
2076  }
2077 
2078 
2079  if (opt.Contains("OF")) {
2080  if (GetDimension() == 3) k_end = ++nbinz1;
2081  if (GetDimension() >= 2) j_end = ++nbiny1;
2082  if (GetDimension() >= 1) i_end = ++nbinx1;
2083  }
2084 
2085  if (opt.Contains("UF")) {
2086  if (GetDimension() == 3) k_start = 0;
2087  if (GetDimension() >= 2) j_start = 0;
2088  if (GetDimension() >= 1) i_start = 0;
2089  }
2090 
2091  ndf = (i_end - i_start + 1) * (j_end - j_start + 1) * (k_end - k_start + 1) - 1;
2092 
2093  Bool_t comparisonUU = opt.Contains("UU");
2094  Bool_t comparisonUW = opt.Contains("UW");
2095  Bool_t comparisonWW = opt.Contains("WW");
2096  Bool_t scaledHistogram = opt.Contains("NORM");
2097 
2098  if (scaledHistogram && !comparisonUU) {
2099  Info("Chi2TestX", "NORM option should be used together with UU option. It is ignored");
2100  }
2101 
2102  // look at histo global bin content and effective entries
2103  Stat_t s[kNstat];
2104  GetStats(s);// s[1] sum of squares of weights, s[0] sum of weights
2105  Double_t sumBinContent1 = s[0];
2106  Double_t effEntries1 = (s[1] ? s[0] * s[0] / s[1] : 0.0);
2107 
2108  h2->GetStats(s);// s[1] sum of squares of weights, s[0] sum of weights
2109  Double_t sumBinContent2 = s[0];
2110  Double_t effEntries2 = (s[1] ? s[0] * s[0] / s[1] : 0.0);
2111 
2112  if (!comparisonUU && !comparisonUW && !comparisonWW ) {
2113  // deduce automatically from type of histogram
2114  if (TMath::Abs(sumBinContent1 - effEntries1) < 1) {
2115  if ( TMath::Abs(sumBinContent2 - effEntries2) < 1) comparisonUU = true;
2116  else comparisonUW = true;
2117  }
2118  else comparisonWW = true;
2119  }
2120  // check unweighted histogram
2121  if (comparisonUW) {
2122  if (TMath::Abs(sumBinContent1 - effEntries1) >= 1) {
2123  Warning("Chi2TestX","First histogram is not unweighted and option UW has been requested");
2124  }
2125  }
2126  if ( (!scaledHistogram && comparisonUU) ) {
2127  if ( ( TMath::Abs(sumBinContent1 - effEntries1) >= 1) || (TMath::Abs(sumBinContent2 - effEntries2) >= 1) ) {
2128  Warning("Chi2TestX","Both histograms are not unweighted and option UU has been requested");
2129  }
2130  }
2131 
2132 
2133  //get number of events in histogram
2134  if (comparisonUU && scaledHistogram) {
2135  for (Int_t i = i_start; i <= i_end; ++i) {
2136  for (Int_t j = j_start; j <= j_end; ++j) {
2137  for (Int_t k = k_start; k <= k_end; ++k) {
2138 
2139  Int_t bin = GetBin(i, j, k);
2140 
2141  Double_t cnt1 = RetrieveBinContent(bin);
2142  Double_t cnt2 = h2->RetrieveBinContent(bin);
2143  Double_t e1sq = GetBinErrorSqUnchecked(bin);
2144  Double_t e2sq = h2->GetBinErrorSqUnchecked(bin);
2145 
2146  if (e1sq > 0.0) cnt1 = TMath::Floor(cnt1 * cnt1 / e1sq + 0.5); // avoid rounding errors
2147  else cnt1 = 0.0;
2148 
2149  if (e2sq > 0.0) cnt2 = TMath::Floor(cnt2 * cnt2 / e2sq + 0.5); // avoid rounding errors
2150  else cnt2 = 0.0;
2151 
2152  // sum contents
2153  sum1 += cnt1;
2154  sum2 += cnt2;
2155  sumw1 += e1sq;
2156  sumw2 += e2sq;
2157  }
2158  }
2159  }
2160  if (sumw1 <= 0.0 || sumw2 <= 0.0) {
2161  Error("Chi2TestX", "Cannot use option NORM when one histogram has all zero errors");
2162  return 0.0;
2163  }
2164 
2165  } else {
2166  for (Int_t i = i_start; i <= i_end; ++i) {
2167  for (Int_t j = j_start; j <= j_end; ++j) {
2168  for (Int_t k = k_start; k <= k_end; ++k) {
2169 
2170  Int_t bin = GetBin(i, j, k);
2171 
2172  sum1 += RetrieveBinContent(bin);
2173  sum2 += h2->RetrieveBinContent(bin);
2174 
2175  if ( comparisonWW ) sumw1 += GetBinErrorSqUnchecked(bin);
2176  if ( comparisonUW || comparisonWW ) sumw2 += h2->GetBinErrorSqUnchecked(bin);
2177  }
2178  }
2179  }
2180  }
2181  //checks that the histograms are not empty
2182  if (sum1 == 0.0 || sum2 == 0.0) {
2183  Error("Chi2TestX","one histogram is empty");
2184  return 0.0;
2185  }
2186 
2187  if ( comparisonWW && ( sumw1 <= 0.0 && sumw2 <= 0.0 ) ){
2188  Error("Chi2TestX","Hist1 and Hist2 have both all zero errors\n");
2189  return 0.0;
2190  }
2191 
2192  //THE TEST
2193  Int_t m = 0, n = 0;
2194 
2195  //Experiment - experiment comparison
2196  if (comparisonUU) {
2197  Double_t sum = sum1 + sum2;
2198  for (Int_t i = i_start; i <= i_end; ++i) {
2199  for (Int_t j = j_start; j <= j_end; ++j) {
2200  for (Int_t k = k_start; k <= k_end; ++k) {
2201 
2202  Int_t bin = GetBin(i, j, k);
2203 
2204  Double_t cnt1 = RetrieveBinContent(bin);
2205  Double_t cnt2 = h2->RetrieveBinContent(bin);
2206 
2207  if (scaledHistogram) {
2208  // scale bin value to effective bin entries
2209  Double_t e1sq = GetBinErrorSqUnchecked(bin);
2210  Double_t e2sq = h2->GetBinErrorSqUnchecked(bin);
2211 
2212  if (e1sq > 0) cnt1 = TMath::Floor(cnt1 * cnt1 / e1sq + 0.5); // avoid rounding errors
2213  else cnt1 = 0;
2214 
2215  if (e2sq > 0) cnt2 = TMath::Floor(cnt2 * cnt2 / e2sq + 0.5); // avoid rounding errors
2216  else cnt2 = 0;
2217  }
2218 
2219  if (Int_t(cnt1) == 0 && Int_t(cnt2) == 0) --ndf; // no data means one degree of freedom less
2220  else {
2221 
2222  Double_t cntsum = cnt1 + cnt2;
2223  Double_t nexp1 = cntsum * sum1 / sum;
2224  //Double_t nexp2 = binsum*sum2/sum;
2225 
2226  if (res) res[i - i_start] = (cnt1 - nexp1) / TMath::Sqrt(nexp1);
2227 
2228  if (cnt1 < 1) ++m;
2229  if (cnt2 < 1) ++n;
2230 
2231  //Habermann correction for residuals
2232  Double_t correc = (1. - sum1 / sum) * (1. - cntsum / sum);
2233  if (res) res[i - i_start] /= TMath::Sqrt(correc);
2234 
2235  Double_t delta = sum2 * cnt1 - sum1 * cnt2;
2236  chi2 += delta * delta / cntsum;
2237  }
2238  }
2239  }
2240  }
2241  chi2 /= sum1 * sum2;
2242 
2243  // flag error only when of the two histogram is zero
2244  if (m) {
2245  igood += 1;
2246  Info("Chi2TestX","There is a bin in h1 with less than 1 event.\n");
2247  }
2248  if (n) {
2249  igood += 2;
2250  Info("Chi2TestX","There is a bin in h2 with less than 1 event.\n");
2251  }
2252 
2253  Double_t prob = TMath::Prob(chi2,ndf);
2254  return prob;
2255 
2256  }
2257 
2258  // unweighted - weighted comparison
2259  // case of error = 0 and content not zero is treated without problems by excluding second chi2 sum
2260  // and can be considered as a data-theory comparison
2261  if ( comparisonUW ) {
2262  for (Int_t i = i_start; i <= i_end; ++i) {
2263  for (Int_t j = j_start; j <= j_end; ++j) {
2264  for (Int_t k = k_start; k <= k_end; ++k) {
2265 
2266  Int_t bin = GetBin(i, j, k);
2267 
2268  Double_t cnt1 = RetrieveBinContent(bin);
2269  Double_t cnt2 = h2->RetrieveBinContent(bin);
2270  Double_t e2sq = h2->GetBinErrorSqUnchecked(bin);
2271 
2272  // case both histogram have zero bin contents
2273  if (cnt1 * cnt1 == 0 && cnt2 * cnt2 == 0) {
2274  --ndf; //no data means one degree of freedom less
2275  continue;
2276  }
2277 
2278  // case weighted histogram has zero bin content and error
2279  if (cnt2 * cnt2 == 0 && e2sq == 0) {
2280  if (sumw2 > 0) {
2281  // use as approximated error as 1 scaled by a scaling ratio
2282  // estimated from the total sum weight and sum weight squared
2283  e2sq = sumw2 / sum2;
2284  }
2285  else {
2286  // return error because infinite discrepancy here:
2287  // bin1 != 0 and bin2 =0 in a histogram with all errors zero
2288  Error("Chi2TestX","Hist2 has in bin (%d,%d,%d) zero content and zero errors\n", i, j, k);
2289  chi2 = 0; return 0;
2290  }
2291  }
2292 
2293  if (cnt1 < 1) m++;
2294  if (e2sq > 0 && cnt2 * cnt2 / e2sq < 10) n++;
2295 
2296  Double_t var1 = sum2 * cnt2 - sum1 * e2sq;
2297  Double_t var2 = var1 * var1 + 4. * sum2 * sum2 * cnt1 * e2sq;
2298 
2299  // if cnt1 is zero and cnt2 = 1 and sum1 = sum2 var1 = 0 && var2 == 0
2300  // approximate by incrementing cnt1
2301  // LM (this need to be fixed for numerical errors)
2302  while (var1 * var1 + cnt1 == 0 || var1 + var2 == 0) {
2303  sum1++;
2304  cnt1++;
2305  var1 = sum2 * cnt2 - sum1 * e2sq;
2306  var2 = var1 * var1 + 4. * sum2 * sum2 * cnt1 * e2sq;
2307  }
2308  var2 = TMath::Sqrt(var2);
2309 
2310  while (var1 + var2 == 0) {
2311  sum1++;
2312  cnt1++;
2313  var1 = sum2 * cnt2 - sum1 * e2sq;
2314  var2 = var1 * var1 + 4. * sum2 * sum2 * cnt1 * e2sq;
2315  while (var1 * var1 + cnt1 == 0 || var1 + var2 == 0) {
2316  sum1++;
2317  cnt1++;
2318  var1 = sum2 * cnt2 - sum1 * e2sq;
2319  var2 = var1 * var1 + 4. * sum2 * sum2 * cnt1 * e2sq;
2320  }
2321  var2 = TMath::Sqrt(var2);
2322  }
2323 
2324  Double_t probb = (var1 + var2) / (2. * sum2 * sum2);
2325 
2326  Double_t nexp1 = probb * sum1;
2327  Double_t nexp2 = probb * sum2;
2328 
2329  Double_t delta1 = cnt1 - nexp1;
2330  Double_t delta2 = cnt2 - nexp2;
2331 
2332  chi2 += delta1 * delta1 / nexp1;
2333 
2334  if (e2sq > 0) {
2335  chi2 += delta2 * delta2 / e2sq;
2336  }
2337 
2338  if (res) {
2339  if (e2sq > 0) {
2340  Double_t temp1 = sum2 * e2sq / var2;
2341  Double_t temp2 = 1.0 + (sum1 * e2sq - sum2 * cnt2) / var2;
2342  temp2 = temp1 * temp1 * sum1 * probb * (1.0 - probb) + temp2 * temp2 * e2sq / 4.0;
2343  // invert sign here
2344  res[i - i_start] = - delta2 / TMath::Sqrt(temp2);
2345  }
2346  else
2347  res[i - i_start] = delta1 / TMath::Sqrt(nexp1);
2348  }
2349  }
2350  }
2351  }
2352 
2353  if (m) {
2354  igood += 1;
2355  Info("Chi2TestX","There is a bin in h1 with less than 1 event.\n");
2356  }
2357  if (n) {
2358  igood += 2;
2359  Info("Chi2TestX","There is a bin in h2 with less than 10 effective events.\n");
2360  }
2361 
2362  Double_t prob = TMath::Prob(chi2, ndf);
2363 
2364  return prob;
2365  }
2366 
2367  // weighted - weighted comparison
2368  if (comparisonWW) {
2369  for (Int_t i = i_start; i <= i_end; ++i) {
2370  for (Int_t j = j_start; j <= j_end; ++j) {
2371  for (Int_t k = k_start; k <= k_end; ++k) {
2372 
2373  Int_t bin = GetBin(i, j, k);
2374  Double_t cnt1 = RetrieveBinContent(bin);
2375  Double_t cnt2 = h2->RetrieveBinContent(bin);
2376  Double_t e1sq = GetBinErrorSqUnchecked(bin);
2377  Double_t e2sq = h2->GetBinErrorSqUnchecked(bin);
2378 
2379  // case both histogram have zero bin contents
2380  // (use square of content to avoid numerical errors)
2381  if (cnt1 * cnt1 == 0 && cnt2 * cnt2 == 0) {
2382  --ndf; //no data means one degree of freedom less
2383  continue;
2384  }
2385 
2386  if (e1sq == 0 && e2sq == 0) {
2387  // cannot treat case of booth histogram have zero zero errors
2388  Error("Chi2TestX","h1 and h2 both have bin %d,%d,%d with all zero errors\n", i,j,k);
2389  chi2 = 0; return 0;
2390  }
2391 
2392  Double_t sigma = sum1 * sum1 * e2sq + sum2 * sum2 * e1sq;
2393  Double_t delta = sum2 * cnt1 - sum1 * cnt2;
2394  chi2 += delta * delta / sigma;
2395 
2396  if (res) {
2397  Double_t temp = cnt1 * sum1 * e2sq + cnt2 * sum2 * e1sq;
2398  Double_t probb = temp / sigma;
2399  Double_t z = 0;
2400  if (e1sq > e2sq) {
2401  Double_t d1 = cnt1 - sum1 * probb;
2402  Double_t s1 = e1sq * ( 1. - e2sq * sum1 * sum1 / sigma );
2403  z = d1 / TMath::Sqrt(s1);
2404  }
2405  else {
2406  Double_t d2 = cnt2 - sum2 * probb;
2407  Double_t s2 = e2sq * ( 1. - e1sq * sum2 * sum2 / sigma );
2408  z = -d2 / TMath::Sqrt(s2);
2409  }
2410  res[i - i_start] = z;
2411  }
2412 
2413  if (e1sq > 0 && cnt1 * cnt1 / e1sq < 10) m++;
2414  if (e2sq > 0 && cnt2 * cnt2 / e2sq < 10) n++;
2415  }
2416  }
2417  }
2418  if (m) {
2419  igood += 1;
2420  Info("Chi2TestX","There is a bin in h1 with less than 10 effective events.\n");
2421  }
2422  if (n) {
2423  igood += 2;
2424  Info("Chi2TestX","There is a bin in h2 with less than 10 effective events.\n");
2425  }
2426  Double_t prob = TMath::Prob(chi2, ndf);
2427  return prob;
2428  }
2429  return 0;
2430 }
2431 ////////////////////////////////////////////////////////////////////////////////
2432 /// Compute and return the chisquare of this histogram with respect to a function
2433 /// The chisquare is computed by weighting each histogram point by the bin error
2434 /// By default the full range of the histogram is used.
2435 /// Use option "R" for restricting the chisquare calculation to the given range of the function
2436 /// Use option "L" for using the chisquare based on the poisson likelihood (Baker-Cousins Chisquare)
2437 
2438 Double_t TH1::Chisquare(TF1 * func, Option_t *option) const
2439 {
2440  if (!func) {
2441  Error("Chisquare","Function pointer is Null - return -1");
2442  return -1;
2443  }
2444 
2445  TString opt(option); opt.ToUpper();
2446  bool useRange = opt.Contains("R");
2447  bool usePL = opt.Contains("L");
2448 
2449  return ROOT::Fit::Chisquare(*this, *func, useRange, usePL);
2450 }
2451 
2452 ////////////////////////////////////////////////////////////////////////////////
2453 /// Remove all the content from the underflow and overflow bins, without changing the number of entries
2454 /// After calling this method, every undeflow and overflow bins will have content 0.0
2455 /// The Sumw2 is also cleared, since there is no more content in the bins
2456 
2458 {
2459  for (Int_t bin = 0; bin < fNcells; ++bin)
2460  if (IsBinUnderflow(bin) || IsBinOverflow(bin)) {
2461  UpdateBinContent(bin, 0.0);
2462  if (fSumw2.fN) fSumw2.fArray[bin] = 0.0;
2463  }
2464 }
2465 
2466 ////////////////////////////////////////////////////////////////////////////////
2467 /// Compute integral (cumulative sum of bins)
2468 /// The result stored in fIntegral is used by the GetRandom functions.
2469 /// This function is automatically called by GetRandom when the fIntegral
2470 /// array does not exist or when the number of entries in the histogram
2471 /// has changed since the previous call to GetRandom.
2472 /// The resulting integral is normalized to 1
2473 /// If the routine is called with the onlyPositive flag set an error will
2474 /// be produced in case of negative bin content and a NaN value returned
2475 
2476 Double_t TH1::ComputeIntegral(Bool_t onlyPositive)
2477 {
2478  if (fBuffer) BufferEmpty();
2479 
2480  // delete previously computed integral (if any)
2481  if (fIntegral) delete [] fIntegral;
2482 
2483  // - Allocate space to store the integral and compute integral
2484  Int_t nbinsx = GetNbinsX();
2485  Int_t nbinsy = GetNbinsY();
2486  Int_t nbinsz = GetNbinsZ();
2487  Int_t nbins = nbinsx * nbinsy * nbinsz;
2488 
2489  fIntegral = new Double_t[nbins + 2];
2490  Int_t ibin = 0; fIntegral[ibin] = 0;
2491 
2492  for (Int_t binz=1; binz <= nbinsz; ++binz) {
2493  for (Int_t biny=1; biny <= nbinsy; ++biny) {
2494  for (Int_t binx=1; binx <= nbinsx; ++binx) {
2495  ++ibin;
2496  Double_t y = RetrieveBinContent(GetBin(binx, biny, binz));
2497  if (onlyPositive && y < 0) {
2498  Error("ComputeIntegral","Bin content is negative - return a NaN value");
2499  fIntegral[nbins] = TMath::QuietNaN();
2500  break;
2501  }
2502  fIntegral[ibin] = fIntegral[ibin - 1] + y;
2503  }
2504  }
2505  }
2506 
2507  // - Normalize integral to 1
2508  if (fIntegral[nbins] == 0 ) {
2509  Error("ComputeIntegral", "Integral = zero"); return 0;
2510  }
2511  for (Int_t bin=1; bin <= nbins; ++bin) fIntegral[bin] /= fIntegral[nbins];
2512  fIntegral[nbins+1] = fEntries;
2513  return fIntegral[nbins];
2514 }
2515 
2516 ////////////////////////////////////////////////////////////////////////////////
2517 /// Return a pointer to the array of bins integral.
2518 /// if the pointer fIntegral is null, TH1::ComputeIntegral is called
2519 /// The array dimension is the number of bins in the histograms
2520 /// including underflow and overflow (fNCells)
2521 /// the last value integral[fNCells] is set to the number of entries of
2522 /// the histogram
2523 
2525 {
2526  if (!fIntegral) ComputeIntegral();
2527  return fIntegral;
2528 }
2529 
2530 ////////////////////////////////////////////////////////////////////////////////
2531 /// Return a pointer to an histogram containing the cumulative The
2532 /// cumulative can be computed both in the forward (default) or backward
2533 /// direction; the name of the new histogram is constructed from
2534 /// the name of this histogram with the suffix suffix appended.
2535 ///
2536 /// The cumulative distribution is formed by filling each bin of the
2537 /// resulting histogram with the sum of that bin and all previous
2538 /// (forward == kTRUE) or following (forward = kFALSE) bins.
2539 ///
2540 /// note: while cumulative distributions make sense in one dimension, you
2541 /// may not be getting what you expect in more than 1D because the concept
2542 /// of a cumulative distribution is much trickier to define; make sure you
2543 /// understand the order of summation before you use this method with
2544 /// histograms of dimension >= 2.
2545 
2546 TH1 *TH1::GetCumulative(Bool_t forward, const char* suffix) const
2547 {
2548  const Int_t nbinsx = GetNbinsX();
2549  const Int_t nbinsy = GetNbinsY();
2550  const Int_t nbinsz = GetNbinsZ();
2551  TH1* hintegrated = (TH1*) Clone(fName + suffix);
2552  hintegrated->Reset();
2553  if (forward) { // Forward computation
2554  Double_t sum = 0.;
2555  for (Int_t binz = 1; binz <= nbinsz; ++binz) {
2556  for (Int_t biny = 1; biny <= nbinsy; ++biny) {
2557  for (Int_t binx = 1; binx <= nbinsx; ++binx) {
2558  const Int_t bin = hintegrated->GetBin(binx, biny, binz);
2559  sum += GetBinContent(bin);
2560  hintegrated->SetBinContent(bin, sum);
2561  }
2562  }
2563  }
2564  } else { // Backward computation
2565  Double_t sum = 0.;
2566  for (Int_t binz = nbinsz; binz >= 1; --binz) {
2567  for (Int_t biny = nbinsy; biny >= 1; --biny) {
2568  for (Int_t binx = nbinsx; binx >= 1; --binx) {
2569  const Int_t bin = hintegrated->GetBin(binx, biny, binz);
2570  sum += GetBinContent(bin);
2571  hintegrated->SetBinContent(bin, sum);
2572  }
2573  }
2574  }
2575  }
2576  return hintegrated;
2577 }
2578 
2579 ////////////////////////////////////////////////////////////////////////////////
2580 /// Copy this histogram structure to newth1.
2581 ///
2582 /// Note that this function does not copy the list of associated functions.
2583 /// Use TObject::Clone to make a full copy of an histogram.
2584 ///
2585 /// Note also that the histogram it will be created in gDirectory (if AddDirectoryStatus()=true)
2586 /// or will not be added to any directory if AddDirectoryStatus()=false
2587 /// independently of the current directory stored in the original histogram
2588 
2589 void TH1::Copy(TObject &obj) const
2590 {
2591  if (((TH1&)obj).fDirectory) {
2592  // We are likely to change the hash value of this object
2593  // with TNamed::Copy, to keep things correct, we need to
2594  // clean up its existing entries.
2595  ((TH1&)obj).fDirectory->Remove(&obj);
2596  ((TH1&)obj).fDirectory = 0;
2597  }
2598  TNamed::Copy(obj);
2599  ((TH1&)obj).fDimension = fDimension;
2600  ((TH1&)obj).fNormFactor= fNormFactor;
2601  ((TH1&)obj).fNcells = fNcells;
2602  ((TH1&)obj).fBarOffset = fBarOffset;
2603  ((TH1&)obj).fBarWidth = fBarWidth;
2604  ((TH1&)obj).fOption = fOption;
2605  ((TH1&)obj).fBinStatErrOpt = fBinStatErrOpt;
2606  ((TH1&)obj).fBufferSize= fBufferSize;
2607  // copy the Buffer
2608  // delete first a previously existing buffer
2609  if (((TH1&)obj).fBuffer != 0) {
2610  delete [] ((TH1&)obj).fBuffer;
2611  ((TH1&)obj).fBuffer = 0;
2612  }
2613  if (fBuffer) {
2614  Double_t *buf = new Double_t[fBufferSize];
2615  for (Int_t i=0;i<fBufferSize;i++) buf[i] = fBuffer[i];
2616  // obj.fBuffer has been deleted before
2617  ((TH1&)obj).fBuffer = buf;
2618  }
2619 
2620 
2621  TArray* a = dynamic_cast<TArray*>(&obj);
2622  if (a) a->Set(fNcells);
2623  for (Int_t i = 0; i < fNcells; i++) ((TH1&)obj).UpdateBinContent(i, RetrieveBinContent(i));
2624 
2625  ((TH1&)obj).fEntries = fEntries;
2626 
2627  // which will call BufferEmpty(0) and set fBuffer[0] to a Maybe one should call
2628  // assignment operator on the TArrayD
2629 
2630  ((TH1&)obj).fTsumw = fTsumw;
2631  ((TH1&)obj).fTsumw2 = fTsumw2;
2632  ((TH1&)obj).fTsumwx = fTsumwx;
2633  ((TH1&)obj).fTsumwx2 = fTsumwx2;
2634  ((TH1&)obj).fMaximum = fMaximum;
2635  ((TH1&)obj).fMinimum = fMinimum;
2636 
2637  TAttLine::Copy(((TH1&)obj));
2638  TAttFill::Copy(((TH1&)obj));
2639  TAttMarker::Copy(((TH1&)obj));
2640  fXaxis.Copy(((TH1&)obj).fXaxis);
2641  fYaxis.Copy(((TH1&)obj).fYaxis);
2642  fZaxis.Copy(((TH1&)obj).fZaxis);
2643  ((TH1&)obj).fXaxis.SetParent(&obj);
2644  ((TH1&)obj).fYaxis.SetParent(&obj);
2645  ((TH1&)obj).fZaxis.SetParent(&obj);
2646  fContour.Copy(((TH1&)obj).fContour);
2647  fSumw2.Copy(((TH1&)obj).fSumw2);
2648  // fFunctions->Copy(((TH1&)obj).fFunctions);
2649  // when copying an histogram if the AddDirectoryStatus() is true it
2650  // will be added to gDirectory independently of the fDirectory stored.
2651  // and if the AddDirectoryStatus() is false it will not be added to
2652  // any directory (fDirectory = 0)
2653  if (fgAddDirectory && gDirectory) {
2654  gDirectory->Append(&obj);
2655  ((TH1&)obj).fFunctions->UseRWLock();
2656  ((TH1&)obj).fDirectory = gDirectory;
2657  } else
2658  ((TH1&)obj).fDirectory = 0;
2659 
2660 }
2661 
2662 ////////////////////////////////////////////////////////////////////////////////
2663 /// Make a complete copy of the underlying object. If 'newname' is set,
2664 /// the copy's name will be set to that name.
2665 
2666 TObject* TH1::Clone(const char* newname) const
2667 {
2668  TH1* obj = (TH1*)IsA()->GetNew()(0);
2669  Copy(*obj);
2670 
2671  // Now handle the parts that Copy doesn't do
2672  if(fFunctions) {
2673  // The Copy above might have published 'obj' to the ListOfCleanups.
2674  // Clone can call RecursiveRemove, for example via TCheckHashRecursiveRemoveConsistency
2675  // when dictionary information is initialized, so we need to
2676  // keep obj->fFunction valid during its execution and
2677  // protect the update with the write lock.
2678 
2679  // Reset stats parent - else cloning the stats will clone this histogram, too.
2680  auto oldstats = dynamic_cast<TVirtualPaveStats*>(fFunctions->FindObject("stats"));
2681  TObject *oldparent = nullptr;
2682  if (oldstats) {
2683  oldparent = oldstats->GetParent();
2684  oldstats->SetParent(nullptr);
2685  }
2686 
2687  auto newlist = (TList*)fFunctions->Clone();
2688 
2689  if (oldstats)
2690  oldstats->SetParent(oldparent);
2691  auto newstats = dynamic_cast<TVirtualPaveStats*>(obj->fFunctions->FindObject("stats"));
2692  if (newstats)
2693  newstats->SetParent(obj);
2694 
2695  auto oldlist = obj->fFunctions;
2696  {
2698  obj->fFunctions = newlist;
2699  }
2700  delete oldlist;
2701  }
2702  if(newname && strlen(newname) ) {
2703  obj->SetName(newname);
2704  }
2705  return obj;
2706 }
2707 
2708 ////////////////////////////////////////////////////////////////////////////////
2709 /// Perform the automatic addition of the histogram to the given directory
2710 ///
2711 /// Note this function is called in place when the semantic requires
2712 /// this object to be added to a directory (I.e. when being read from
2713 /// a TKey or being Cloned)
2714 
2716 {
2717  Bool_t addStatus = TH1::AddDirectoryStatus();
2718  if (addStatus) {
2719  SetDirectory(dir);
2720  if (dir) {
2722  }
2723  }
2724 }
2725 
2726 ////////////////////////////////////////////////////////////////////////////////
2727 /// Compute distance from point px,py to a line.
2728 ///
2729 /// Compute the closest distance of approach from point px,py to elements
2730 /// of an histogram.
2731 /// The distance is computed in pixels units.
2732 ///
2733 /// #### Algorithm:
2734 /// Currently, this simple model computes the distance from the mouse
2735 /// to the histogram contour only.
2736 
2738 {
2739  if (!fPainter) return 9999;
2740  return fPainter->DistancetoPrimitive(px,py);
2741 }
2742 
2743 ////////////////////////////////////////////////////////////////////////////////
2744 /// Performs the operation: `this = this/(c1*f1)`
2745 /// if errors are defined (see TH1::Sumw2), errors are also recalculated.
2746 ///
2747 /// Only bins inside the function range are recomputed.
2748 /// IMPORTANT NOTE: If you intend to use the errors of this histogram later
2749 /// you should call Sumw2 before making this operation.
2750 /// This is particularly important if you fit the histogram after TH1::Divide
2751 ///
2752 /// The function return kFALSE if the divide operation failed
2753 
2755 {
2756  if (!f1) {
2757  Error("Divide","Attempt to divide by a non-existing function");
2758  return kFALSE;
2759  }
2760 
2761  // delete buffer if it is there since it will become invalid
2762  if (fBuffer) BufferEmpty(1);
2763 
2764  Int_t nx = GetNbinsX() + 2; // normal bins + uf / of
2765  Int_t ny = GetNbinsY() + 2;
2766  Int_t nz = GetNbinsZ() + 2;
2767  if (fDimension < 2) ny = 1;
2768  if (fDimension < 3) nz = 1;
2769 
2770 
2771  SetMinimum();
2772  SetMaximum();
2773 
2774  // - Loop on bins (including underflows/overflows)
2775  Int_t bin, binx, biny, binz;
2776  Double_t cu, w;
2777  Double_t xx[3];
2778  Double_t *params = 0;
2779  f1->InitArgs(xx,params);
2780  for (binz = 0; binz < nz; ++binz) {
2781  xx[2] = fZaxis.GetBinCenter(binz);
2782  for (biny = 0; biny < ny; ++biny) {
2783  xx[1] = fYaxis.GetBinCenter(biny);
2784  for (binx = 0; binx < nx; ++binx) {
2785  xx[0] = fXaxis.GetBinCenter(binx);
2786  if (!f1->IsInside(xx)) continue;
2788  bin = binx + nx * (biny + ny * binz);
2789  cu = c1 * f1->EvalPar(xx);
2790  if (TF1::RejectedPoint()) continue;
2791  if (cu) w = RetrieveBinContent(bin) / cu;
2792  else w = 0;
2793  UpdateBinContent(bin, w);
2794  if (fSumw2.fN) {
2795  if (cu != 0) fSumw2.fArray[bin] = GetBinErrorSqUnchecked(bin) / (cu * cu);
2796  else fSumw2.fArray[bin] = 0;
2797  }
2798  }
2799  }
2800  }
2801  ResetStats();
2802  return kTRUE;
2803 }
2804 
2805 ////////////////////////////////////////////////////////////////////////////////
2806 /// Divide this histogram by h1.
2807 ///
2808 /// `this = this/h1`
2809 /// if errors are defined (see TH1::Sumw2), errors are also recalculated.
2810 /// Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
2811 /// if not already set.
2812 /// The resulting errors are calculated assuming uncorrelated histograms.
2813 /// See the other TH1::Divide that gives the possibility to optionally
2814 /// compute binomial errors.
2815 ///
2816 /// IMPORTANT NOTE: If you intend to use the errors of this histogram later
2817 /// you should call Sumw2 before making this operation.
2818 /// This is particularly important if you fit the histogram after TH1::Scale
2819 ///
2820 /// The function return kFALSE if the divide operation failed
2821 
2822 Bool_t TH1::Divide(const TH1 *h1)
2823 {
2824  if (!h1) {
2825  Error("Divide", "Input histogram passed does not exist (NULL).");
2826  return kFALSE;
2827  }
2828 
2829  // delete buffer if it is there since it will become invalid
2830  if (fBuffer) BufferEmpty(1);
2831 
2832  try {
2833  CheckConsistency(this,h1);
2834  } catch(DifferentNumberOfBins&) {
2835  Error("Divide","Cannot divide histograms with different number of bins");
2836  return kFALSE;
2837  } catch(DifferentAxisLimits&) {
2838  Warning("Divide","Dividing histograms with different axis limits");
2839  } catch(DifferentBinLimits&) {
2840  Warning("Divide","Dividing histograms with different bin limits");
2841  } catch(DifferentLabels&) {
2842  Warning("Divide","Dividing histograms with different labels");
2843  }
2844 
2845  // Create Sumw2 if h1 has Sumw2 set
2846  if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
2847 
2848  // - Loop on bins (including underflows/overflows)
2849  for (Int_t i = 0; i < fNcells; ++i) {
2850  Double_t c0 = RetrieveBinContent(i);
2852  if (c1) UpdateBinContent(i, c0 / c1);
2853  else UpdateBinContent(i, 0);
2854 
2855  if(fSumw2.fN) {
2856  if (c1 == 0) { fSumw2.fArray[i] = 0; continue; }
2857  Double_t c1sq = c1 * c1;
2858  fSumw2.fArray[i] = (GetBinErrorSqUnchecked(i) * c1sq + h1->GetBinErrorSqUnchecked(i) * c0 * c0) / (c1sq * c1sq);
2859  }
2860  }
2861  ResetStats();
2862  return kTRUE;
2863 }
2864 
2865 ////////////////////////////////////////////////////////////////////////////////
2866 /// Replace contents of this histogram by the division of h1 by h2.
2867 ///
2868 /// `this = c1*h1/(c2*h2)`
2869 ///
2870 /// If errors are defined (see TH1::Sumw2), errors are also recalculated
2871 /// Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this
2872 /// if not already set.
2873 /// The resulting errors are calculated assuming uncorrelated histograms.
2874 /// However, if option ="B" is specified, Binomial errors are computed.
2875 /// In this case c1 and c2 do not make real sense and they are ignored.
2876 ///
2877 /// IMPORTANT NOTE: If you intend to use the errors of this histogram later
2878 /// you should call Sumw2 before making this operation.
2879 /// This is particularly important if you fit the histogram after TH1::Divide
2880 ///
2881 /// Please note also that in the binomial case errors are calculated using standard
2882 /// binomial statistics, which means when b1 = b2, the error is zero.
2883 /// If you prefer to have efficiency errors not going to zero when the efficiency is 1, you must
2884 /// use the function TGraphAsymmErrors::BayesDivide, which will return an asymmetric and non-zero lower
2885 /// error for the case b1=b2.
2886 ///
2887 /// The function return kFALSE if the divide operation failed
2888 
2889 Bool_t TH1::Divide(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
2890 {
2891 
2892  TString opt = option;
2893  opt.ToLower();
2894  Bool_t binomial = kFALSE;
2895  if (opt.Contains("b")) binomial = kTRUE;
2896  if (!h1 || !h2) {
2897  Error("Divide", "At least one of the input histograms passed does not exist (NULL).");
2898  return kFALSE;
2899  }
2900 
2901  // delete buffer if it is there since it will become invalid
2902  if (fBuffer) BufferEmpty(1);
2903 
2904  try {
2905  CheckConsistency(h1,h2);
2906  CheckConsistency(this,h1);
2907  } catch(DifferentNumberOfBins&) {
2908  Error("Divide","Cannot divide histograms with different number of bins");
2909  return kFALSE;
2910  } catch(DifferentAxisLimits&) {
2911  Warning("Divide","Dividing histograms with different axis limits");
2912  } catch(DifferentBinLimits&) {
2913  Warning("Divide","Dividing histograms with different bin limits");
2914  } catch(DifferentLabels&) {
2915  Warning("Divide","Dividing histograms with different labels");
2916  }
2917 
2918 
2919  if (!c2) {
2920  Error("Divide","Coefficient of dividing histogram cannot be zero");
2921  return kFALSE;
2922  }
2923 
2924  // Create Sumw2 if h1 or h2 have Sumw2 set, or if binomial errors are explicitly requested
2925  if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0 || binomial)) Sumw2();
2926 
2927  SetMinimum();
2928  SetMaximum();
2929 
2930  // - Loop on bins (including underflows/overflows)
2931  for (Int_t i = 0; i < fNcells; ++i) {
2932  Double_t b1 = h1->RetrieveBinContent(i);
2933  Double_t b2 = h2->RetrieveBinContent(i);
2934  if (b2) UpdateBinContent(i, c1 * b1 / (c2 * b2));
2935  else UpdateBinContent(i, 0);
2936 
2937  if (fSumw2.fN) {
2938  if (b2 == 0) { fSumw2.fArray[i] = 0; continue; }
2939  Double_t b1sq = b1 * b1; Double_t b2sq = b2 * b2;
2940  Double_t c1sq = c1 * c1; Double_t c2sq = c2 * c2;
2941  Double_t e1sq = h1->GetBinErrorSqUnchecked(i);
2942  Double_t e2sq = h2->GetBinErrorSqUnchecked(i);
2943  if (binomial) {
2944  if (b1 != b2) {
2945  // in the case of binomial statistics c1 and c2 must be 1 otherwise it does not make sense
2946  // c1 and c2 are ignored
2947  //fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/(c2*b2));//this is the formula in Hbook/Hoper1
2948  //fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/b2); // old formula from G. Flucke
2949  // formula which works also for weighted histogram (see http://root-forum.cern.ch/viewtopic.php?t=3753 )
2950  fSumw2.fArray[i] = TMath::Abs( ( (1. - 2.* b1 / b2) * e1sq + b1sq * e2sq / b2sq ) / b2sq );
2951  } else {
2952  //in case b1=b2 error is zero
2953  //use TGraphAsymmErrors::BayesDivide for getting the asymmetric error not equal to zero
2954  fSumw2.fArray[i] = 0;
2955  }
2956  } else {
2957  fSumw2.fArray[i] = c1sq * c2sq * (e1sq * b2sq + e2sq * b1sq) / (c2sq * c2sq * b2sq * b2sq);
2958  }
2959  }
2960  }
2961  ResetStats();
2962  if (binomial)
2963  // in case of binomial division use denominator for number of entries
2964  SetEntries ( h2->GetEntries() );
2965 
2966  return kTRUE;
2967 }
2968 
2969 ////////////////////////////////////////////////////////////////////////////////
2970 /// Draw this histogram with options.
2971 ///
2972 /// Histograms are drawn via the THistPainter class. Each histogram has
2973 /// a pointer to its own painter (to be usable in a multithreaded program).
2974 /// The same histogram can be drawn with different options in different pads.
2975 /// When an histogram drawn in a pad is deleted, the histogram is
2976 /// automatically removed from the pad or pads where it was drawn.
2977 /// If an histogram is drawn in a pad, then filled again, the new status
2978 /// of the histogram will be automatically shown in the pad next time
2979 /// the pad is updated. One does not need to redraw the histogram.
2980 /// To draw the current version of an histogram in a pad, one can use
2981 /// `h->DrawCopy();`
2982 /// This makes a clone of the histogram. Once the clone is drawn, the original
2983 /// histogram may be modified or deleted without affecting the aspect of the
2984 /// clone.
2985 /// By default, TH1::Draw clears the current pad.
2986 ///
2987 /// One can use TH1::SetMaximum and TH1::SetMinimum to force a particular
2988 /// value for the maximum or the minimum scale on the plot.
2989 ///
2990 /// TH1::UseCurrentStyle can be used to change all histogram graphics
2991 /// attributes to correspond to the current selected style.
2992 /// This function must be called for each histogram.
2993 /// In case one reads and draws many histograms from a file, one can force
2994 /// the histograms to inherit automatically the current graphics style
2995 /// by calling before gROOT->ForceStyle();
2996 ///
2997 /// See the THistPainter class for a description of all the drawing options.
2998 
2999 void TH1::Draw(Option_t *option)
3000 {
3001  TString opt1 = option; opt1.ToLower();
3002  TString opt2 = option;
3003  Int_t index = opt1.Index("same");
3004 
3005  // Check if the string "same" is part of a TCutg name.
3006  if (index>=0) {
3007  Int_t indb = opt1.Index("[");
3008  if (indb>=0) {
3009  Int_t indk = opt1.Index("]");
3010  if (index>indb && index<indk) index = -1;
3011  }
3012  }
3013 
3014  // If there is no pad or an empty pad the "same" option is ignored.
3015  if (gPad) {
3016  if (!gPad->IsEditable()) gROOT->MakeDefCanvas();
3017  if (index>=0) {
3018  if (gPad->GetX1() == 0 && gPad->GetX2() == 1 &&
3019  gPad->GetY1() == 0 && gPad->GetY2() == 1 &&
3020  gPad->GetListOfPrimitives()->GetSize()==0) opt2.Remove(index,4);
3021  } else {
3022  //the following statement is necessary in case one attempts to draw
3023  //a temporary histogram already in the current pad
3024  if (TestBit(kCanDelete)) gPad->GetListOfPrimitives()->Remove(this);
3025  gPad->Clear();
3026  }
3027  gPad->IncrementPaletteColor(1, opt1);
3028  } else {
3029  if (index>=0) opt2.Remove(index,4);
3030  }
3031 
3032  AppendPad(opt2.Data());
3033 }
3034 
3035 ////////////////////////////////////////////////////////////////////////////////
3036 /// Copy this histogram and Draw in the current pad.
3037 ///
3038 /// Once the histogram is drawn into the pad, any further modification
3039 /// using graphics input will be made on the copy of the histogram,
3040 /// and not to the original object.
3041 /// By default a postfix "_copy" is added to the histogram name. Pass an empty postfix in case
3042 /// you want to draw an histogram with the same name
3043 ///
3044 /// See Draw for the list of options
3045 
3046 TH1 *TH1::DrawCopy(Option_t *option, const char * name_postfix) const
3047 {
3048  TString opt = option;
3049  opt.ToLower();
3050  if (gPad && !opt.Contains("same")) gPad->Clear();
3051  TString newName = (name_postfix) ? TString::Format("%s%s",GetName(),name_postfix) : "";
3052  TH1 *newth1 = (TH1 *)Clone(newName);
3053  newth1->SetDirectory(0);
3054  newth1->SetBit(kCanDelete);
3055  if (gPad) gPad->IncrementPaletteColor(1, opt);
3056 
3057  newth1->AppendPad(option);
3058  return newth1;
3059 }
3060 
3061 ////////////////////////////////////////////////////////////////////////////////
3062 /// Draw a normalized copy of this histogram.
3063 ///
3064 /// A clone of this histogram is normalized to norm and drawn with option.
3065 /// A pointer to the normalized histogram is returned.
3066 /// The contents of the histogram copy are scaled such that the new
3067 /// sum of weights (excluding under and overflow) is equal to norm.
3068 /// Note that the returned normalized histogram is not added to the list
3069 /// of histograms in the current directory in memory.
3070 /// It is the user's responsibility to delete this histogram.
3071 /// The kCanDelete bit is set for the returned object. If a pad containing
3072 /// this copy is cleared, the histogram will be automatically deleted.
3073 ///
3074 /// See Draw for the list of options
3075 
3076 TH1 *TH1::DrawNormalized(Option_t *option, Double_t norm) const
3077 {
3079  if (sum == 0) {
3080  Error("DrawNormalized","Sum of weights is null. Cannot normalize histogram: %s",GetName());
3081  return 0;
3082  }
3083  Bool_t addStatus = TH1::AddDirectoryStatus();
3085  TH1 *h = (TH1*)Clone();
3086  h->SetBit(kCanDelete);
3087  // in case of drawing with error options - scale correctly the error
3088  TString opt(option); opt.ToUpper();
3089  if (fSumw2.fN == 0) {
3090  h->Sumw2();
3091  // do not use in this case the "Error option " for drawing which is enabled by default since the normalized histogram has now errors
3092  if (opt.IsNull() || opt == "SAME") opt += "HIST";
3093  }
3094  h->Scale(norm/sum);
3095  if (TMath::Abs(fMaximum+1111) > 1e-3) h->SetMaximum(fMaximum*norm/sum);
3096  if (TMath::Abs(fMinimum+1111) > 1e-3) h->SetMinimum(fMinimum*norm/sum);
3097  h->Draw(opt);
3098  TH1::AddDirectory(addStatus);
3099  return h;
3100 }
3101 
3102 ////////////////////////////////////////////////////////////////////////////////
3103 /// Display a panel with all histogram drawing options.
3104 ///
3105 /// See class TDrawPanelHist for example
3106 
3107 void TH1::DrawPanel()
3108 {
3109  if (!fPainter) {Draw(); if (gPad) gPad->Update();}
3110  if (fPainter) fPainter->DrawPanel();
3111 }
3112 
3113 ////////////////////////////////////////////////////////////////////////////////
3114 /// Evaluate function f1 at the center of bins of this histogram.
3115 ///
3116 /// - If option "R" is specified, the function is evaluated only
3117 /// for the bins included in the function range.
3118 /// - If option "A" is specified, the value of the function is added to the
3119 /// existing bin contents
3120 /// - If option "S" is specified, the value of the function is used to
3121 /// generate a value, distributed according to the Poisson
3122 /// distribution, with f1 as the mean.
3123 
3124 void TH1::Eval(TF1 *f1, Option_t *option)
3125 {
3126  Double_t x[3];
3127  Int_t range, stat, add;
3128  if (!f1) return;
3129 
3130  TString opt = option;
3131  opt.ToLower();
3132  if (opt.Contains("a")) add = 1;
3133  else add = 0;
3134  if (opt.Contains("s")) stat = 1;
3135  else stat = 0;
3136  if (opt.Contains("r")) range = 1;
3137  else range = 0;
3138 
3139  // delete buffer if it is there since it will become invalid
3140  if (fBuffer) BufferEmpty(1);
3141 
3142  Int_t nbinsx = fXaxis.GetNbins();
3143  Int_t nbinsy = fYaxis.GetNbins();
3144  Int_t nbinsz = fZaxis.GetNbins();
3145  if (!add) Reset();
3146 
3147  for (Int_t binz = 1; binz <= nbinsz; ++binz) {
3148  x[2] = fZaxis.GetBinCenter(binz);
3149  for (Int_t biny = 1; biny <= nbinsy; ++biny) {
3150  x[1] = fYaxis.GetBinCenter(biny);
3151  for (Int_t binx = 1; binx <= nbinsx; ++binx) {
3152  Int_t bin = GetBin(binx,biny,binz);
3153  x[0] = fXaxis.GetBinCenter(binx);
3154  if (range && !f1->IsInside(x)) continue;
3155  Double_t fu = f1->Eval(x[0], x[1], x[2]);
3156  if (stat) fu = gRandom->PoissonD(fu);
3157  AddBinContent(bin, fu);
3158  if (fSumw2.fN) fSumw2.fArray[bin] += TMath::Abs(fu);
3159  }
3160  }
3161  }
3162 }
3163 
3164 ////////////////////////////////////////////////////////////////////////////////
3165 /// Execute action corresponding to one event.
3166 ///
3167 /// This member function is called when a histogram is clicked with the locator
3168 ///
3169 /// If Left button clicked on the bin top value, then the content of this bin
3170 /// is modified according to the new position of the mouse when it is released.
3171 
3172 void TH1::ExecuteEvent(Int_t event, Int_t px, Int_t py)
3173 {
3174  if (fPainter) fPainter->ExecuteEvent(event, px, py);
3175 }
3176 
3177 ////////////////////////////////////////////////////////////////////////////////
3178 /// This function allows to do discrete Fourier transforms of TH1 and TH2.
3179 /// Available transform types and flags are described below.
3180 ///
3181 /// To extract more information about the transform, use the function
3182 /// TVirtualFFT::GetCurrentTransform() to get a pointer to the current
3183 /// transform object.
3184 ///
3185 /// \param[out] h_output histogram for the output. If a null pointer is passed, a new histogram is created
3186 /// and returned, otherwise, the provided histogram is used and should be big enough
3187 /// \param[in] option option parameters consists of 3 parts:
3188 /// - option on what to return
3189 /// - "RE" - returns a histogram of the real part of the output
3190 /// - "IM" - returns a histogram of the imaginary part of the output
3191 /// - "MAG"- returns a histogram of the magnitude of the output
3192 /// - "PH" - returns a histogram of the phase of the output
3193 /// - option of transform type
3194 /// - "R2C" - real to complex transforms - default
3195 /// - "R2HC" - real to halfcomplex (special format of storing output data,
3196 /// results the same as for R2C)
3197 /// - "DHT" - discrete Hartley transform
3198 /// real to real transforms (sine and cosine):
3199 /// - "R2R_0", "R2R_1", "R2R_2", "R2R_3" - discrete cosine transforms of types I-IV
3200 /// - "R2R_4", "R2R_5", "R2R_6", "R2R_7" - discrete sine transforms of types I-IV
3201 /// To specify the type of each dimension of a 2-dimensional real to real
3202 /// transform, use options of form "R2R_XX", for example, "R2R_02" for a transform,
3203 /// which is of type "R2R_0" in 1st dimension and "R2R_2" in the 2nd.
3204 /// - option of transform flag
3205 /// - "ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
3206 /// performance
3207 /// - "M" (from "measure") - some time spend in finding the optimal way to do the transform
3208 /// - "P" (from "patient") - more time spend in finding the optimal way to do the transform
3209 /// - "EX" (from "exhaustive") - the most optimal way is found
3210 /// This option should be chosen depending on how many transforms of the same size and
3211 /// type are going to be done. Planning is only done once, for the first transform of this
3212 /// size and type. Default is "ES".
3213 ///
3214 /// Examples of valid options: "Mag R2C M" "Re R2R_11" "Im R2C ES" "PH R2HC EX"
3215 
3216 TH1* TH1::FFT(TH1* h_output, Option_t *option)
3217 {
3218 
3219  Int_t ndim[3];
3220  ndim[0] = this->GetNbinsX();
3221  ndim[1] = this->GetNbinsY();
3222  ndim[2] = this->GetNbinsZ();
3223 
3224  TVirtualFFT *fft;
3225  TString opt = option;
3226  opt.ToUpper();
3227  if (!opt.Contains("2R")){
3228  if (!opt.Contains("2C") && !opt.Contains("2HC") && !opt.Contains("DHT")) {
3229  //no type specified, "R2C" by default
3230  opt.Append("R2C");
3231  }
3232  fft = TVirtualFFT::FFT(this->GetDimension(), ndim, opt.Data());
3233  }
3234  else {
3235  //find the kind of transform
3236  Int_t ind = opt.Index("R2R", 3);
3237  Int_t *kind = new Int_t[2];
3238  char t;
3239  t = opt[ind+4];
3240  kind[0] = atoi(&t);
3241  if (h_output->GetDimension()>1) {
3242  t = opt[ind+5];
3243  kind[1] = atoi(&t);
3244  }
3245  fft = TVirtualFFT::SineCosine(this->GetDimension(), ndim, kind, option);
3246  delete [] kind;
3247  }
3248 
3249  if (!fft) return 0;
3250  Int_t in=0;
3251  for (Int_t binx = 1; binx<=ndim[0]; binx++) {
3252  for (Int_t biny=1; biny<=ndim[1]; biny++) {
3253  for (Int_t binz=1; binz<=ndim[2]; binz++) {
3254  fft->SetPoint(in, this->GetBinContent(binx, biny, binz));
3255  in++;
3256  }
3257  }
3258  }
3259  fft->Transform();
3260  h_output = TransformHisto(fft, h_output, option);
3261  return h_output;
3262 }
3263 
3264 ////////////////////////////////////////////////////////////////////////////////
3265 /// Increment bin with abscissa X by 1.
3266 ///
3267 /// if x is less than the low-edge of the first bin, the Underflow bin is incremented
3268 /// if x is equal to or greater than the upper edge of last bin, the Overflow bin is incremented
3269 ///
3270 /// If the storage of the sum of squares of weights has been triggered,
3271 /// via the function Sumw2, then the sum of the squares of weights is incremented
3272 /// by 1 in the bin corresponding to x.
3273 ///
3274 /// The function returns the corresponding bin number which has its content incremented by 1
3275 
3277 {
3278  if (fBuffer) return BufferFill(x,1);
3279 
3280  Int_t bin;
3281  fEntries++;
3282  bin =fXaxis.FindBin(x);
3283  if (bin <0) return -1;
3284  AddBinContent(bin);
3285  if (fSumw2.fN) ++fSumw2.fArray[bin];
3286  if (bin == 0 || bin > fXaxis.GetNbins()) {
3287  if (!GetStatOverflowsBehaviour()) return -1;
3288  }
3289  ++fTsumw;
3290  ++fTsumw2;
3291  fTsumwx += x;
3292  fTsumwx2 += x*x;
3293  return bin;
3294 }
3295 
3296 ////////////////////////////////////////////////////////////////////////////////
3297 /// Increment bin with abscissa X with a weight w.
3298 ///
3299 /// if x is less than the low-edge of the first bin, the Underflow bin is incremented
3300 /// if x is equal to or greater than the upper edge of last bin, the Overflow bin is incremented
3301 ///
3302 /// If the weight is not equal to 1, the storage of the sum of squares of
3303 /// weights is automatically triggered and the sum of the squares of weights is incremented
3304 /// by \f$ w^2 \f$ in the bin corresponding to x.
3305 ///
3306 /// The function returns the corresponding bin number which has its content incremented by w
3307 
3309 {
3310 
3311  if (fBuffer) return BufferFill(x,w);
3312 
3313  Int_t bin;
3314  fEntries++;
3315  bin =fXaxis.FindBin(x);
3316  if (bin <0) return -1;
3317  if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW) ) Sumw2(); // must be called before AddBinContent
3318  if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
3319  AddBinContent(bin, w);
3320  if (bin == 0 || bin > fXaxis.GetNbins()) {
3321  if (!GetStatOverflowsBehaviour()) return -1;
3322  }
3323  Double_t z= w;
3324  fTsumw += z;
3325  fTsumw2 += z*z;
3326  fTsumwx += z*x;
3327  fTsumwx2 += z*x*x;
3328  return bin;
3329 }
3330 
3331 ////////////////////////////////////////////////////////////////////////////////
3332 /// Increment bin with namex with a weight w
3333 ///
3334 /// if x is less than the low-edge of the first bin, the Underflow bin is incremented
3335 /// if x is equal to or greater than the upper edge of last bin, the Overflow bin is incremented
3336 ///
3337 /// If the weight is not equal to 1, the storage of the sum of squares of
3338 /// weights is automatically triggered and the sum of the squares of weights is incremented
3339 /// by \f$ w^2 \f$ in the bin corresponding to x.
3340 ///
3341 /// The function returns the corresponding bin number which has its content
3342 /// incremented by w.
3343 
3344 Int_t TH1::Fill(const char *namex, Double_t w)
3345 {
3346  Int_t bin;
3347  fEntries++;
3348  bin =fXaxis.FindBin(namex);
3349  if (bin <0) return -1;
3350  if (!fSumw2.fN && w != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2();
3351  if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
3352  AddBinContent(bin, w);
3353  if (bin == 0 || bin > fXaxis.GetNbins()) return -1;
3354  Double_t z= w;
3355  fTsumw += z;
3356  fTsumw2 += z*z;
3357  // this make sense if the histogram is not expanding (no axis can be extended)
3358  if (!CanExtendAllAxes()) {
3359  Double_t x = fXaxis.GetBinCenter(bin);
3360  fTsumwx += z*x;
3361  fTsumwx2 += z*x*x;
3362  }
3363  return bin;
3364 }
3365 
3366 ////////////////////////////////////////////////////////////////////////////////
3367 /// Fill this histogram with an array x and weights w.
3368 ///
3369 /// \param[in] ntimes number of entries in arrays x and w (array size must be ntimes*stride)
3370 /// \param[in] x array of values to be histogrammed
3371 /// \param[in] w array of weighs
3372 /// \param[in] stride step size through arrays x and w
3373 ///
3374 /// If the weight is not equal to 1, the storage of the sum of squares of
3375 /// weights is automatically triggered and the sum of the squares of weights is incremented
3376 /// by \f$ w^2 \f$ in the bin corresponding to x.
3377 /// if w is NULL each entry is assumed a weight=1
3378 
3379 void TH1::FillN(Int_t ntimes, const Double_t *x, const Double_t *w, Int_t stride)
3380 {
3381  //If a buffer is activated, fill buffer
3382  if (fBuffer) {
3383  ntimes *= stride;
3384  Int_t i = 0;
3385  for (i=0;i<ntimes;i+=stride) {
3386  if (!fBuffer) break; // buffer can be deleted in BufferFill when is empty
3387  if (w) BufferFill(x[i],w[i]);
3388  else BufferFill(x[i], 1.);
3389  }
3390  // fill the remaining entries if the buffer has been deleted
3391  if (i < ntimes && fBuffer==0) {
3392  auto weights = w ? &w[i] : nullptr;
3393  DoFillN((ntimes-i)/stride,&x[i],weights,stride);
3394  }
3395  return;
3396  }
3397  // call internal method
3398  DoFillN(ntimes, x, w, stride);
3399 }
3400 
3401 ////////////////////////////////////////////////////////////////////////////////
3402 /// Internal method to fill histogram content from a vector
3403 /// called directly by TH1::BufferEmpty
3404 
3405 void TH1::DoFillN(Int_t ntimes, const Double_t *x, const Double_t *w, Int_t stride)
3406 {
3407  Int_t bin,i;
3408 
3409  fEntries += ntimes;
3410  Double_t ww = 1;
3411  Int_t nbins = fXaxis.GetNbins();
3412  ntimes *= stride;
3413  for (i=0;i<ntimes;i+=stride) {
3414  bin =fXaxis.FindBin(x[i]);
3415  if (bin <0) continue;
3416  if (w) ww = w[i];
3417  if (!fSumw2.fN && ww != 1.0 && !TestBit(TH1::kIsNotW)) Sumw2();
3418  if (fSumw2.fN) fSumw2.fArray[bin] += ww*ww;
3419  AddBinContent(bin, ww);
3420  if (bin == 0 || bin > nbins) {
3421  if (!GetStatOverflowsBehaviour()) continue;
3422  }
3423  Double_t z= ww;
3424  fTsumw += z;
3425  fTsumw2 += z*z;
3426  fTsumwx += z*x[i];
3427  fTsumwx2 += z*x[i]*x[i];
3428  }
3429 }
3430 
3431 ////////////////////////////////////////////////////////////////////////////////
3432 /// Fill histogram following distribution in function fname.
3433 ///
3434 /// The distribution contained in the function fname (TF1) is integrated
3435 /// over the channel contents for the bin range of this histogram.
3436 /// It is normalized to 1.
3437 ///
3438 /// Getting one random number implies:
3439 /// - Generating a random number between 0 and 1 (say r1)
3440 /// - Look in which bin in the normalized integral r1 corresponds to
3441 /// - Fill histogram channel
3442 /// ntimes random numbers are generated
3443 ///
3444 /// One can also call TF1::GetRandom to get a random variate from a function.
3445 
3446 void TH1::FillRandom(const char *fname, Int_t ntimes)
3447 {
3448  Int_t bin, binx, ibin, loop;
3449  Double_t r1, x;
3450  // - Search for fname in the list of ROOT defined functions
3451  TF1 *f1 = (TF1*)gROOT->GetFunction(fname);
3452  if (!f1) { Error("FillRandom", "Unknown function: %s",fname); return; }
3453 
3454  // - Allocate temporary space to store the integral and compute integral
3455 
3456  TAxis * xAxis = &fXaxis;
3457 
3458  // in case axis of histogram is not defined use the function axis
3459  if (fXaxis.GetXmax() <= fXaxis.GetXmin()) {
3460  Double_t xmin,xmax;
3461  f1->GetRange(xmin,xmax);
3462  Info("FillRandom","Using function axis and range [%g,%g]",xmin, xmax);
3463  xAxis = f1->GetHistogram()->GetXaxis();
3464  }
3465 
3466  Int_t first = xAxis->GetFirst();
3467  Int_t last = xAxis->GetLast();
3468  Int_t nbinsx = last-first+1;
3469 
3470  Double_t *integral = new Double_t[nbinsx+1];
3471  integral[0] = 0;
3472  for (binx=1;binx<=nbinsx;binx++) {
3473  Double_t fint = f1->Integral(xAxis->GetBinLowEdge(binx+first-1),xAxis->GetBinUpEdge(binx+first-1), 0.);
3474  integral[binx] = integral[binx-1] + fint;
3475  }
3476 
3477  // - Normalize integral to 1
3478  if (integral[nbinsx] == 0 ) {
3479  delete [] integral;
3480  Error("FillRandom", "Integral = zero"); return;
3481  }
3482  for (bin=1;bin<=nbinsx;bin++) integral[bin] /= integral[nbinsx];
3483 
3484  // --------------Start main loop ntimes
3485  for (loop=0;loop<ntimes;loop++) {
3486  r1 = gRandom->Rndm();
3487  ibin = TMath::BinarySearch(nbinsx,&integral[0],r1);
3488  //binx = 1 + ibin;
3489  //x = xAxis->GetBinCenter(binx); //this is not OK when SetBuffer is used
3490  x = xAxis->GetBinLowEdge(ibin+first)
3491  +xAxis->GetBinWidth(ibin+first)*(r1-integral[ibin])/(integral[ibin+1] - integral[ibin]);
3492  Fill(x);
3493  }
3494  delete [] integral;
3495 }
3496 
3497 ////////////////////////////////////////////////////////////////////////////////
3498 /// Fill histogram following distribution in histogram h.
3499 ///
3500 /// The distribution contained in the histogram h (TH1) is integrated
3501 /// over the channel contents for the bin range of this histogram.
3502 /// It is normalized to 1.
3503 ///
3504 /// Getting one random number implies:
3505 /// - Generating a random number between 0 and 1 (say r1)
3506 /// - Look in which bin in the normalized integral r1 corresponds to
3507 /// - Fill histogram channel ntimes random numbers are generated
3508 ///
3509 /// SPECIAL CASE when the target histogram has the same binning as the source.
3510 /// in this case we simply use a poisson distribution where
3511 /// the mean value per bin = bincontent/integral.
3512 
3513 void TH1::FillRandom(TH1 *h, Int_t ntimes)
3514 {
3515  if (!h) { Error("FillRandom", "Null histogram"); return; }
3516  if (fDimension != h->GetDimension()) {
3517  Error("FillRandom", "Histograms with different dimensions"); return;
3518  }
3519  if (std::isnan(h->ComputeIntegral(true))) {
3520  Error("FillRandom", "Histograms contains negative bins, does not represent probabilities");
3521  return;
3522  }
3523 
3524  //in case the target histogram has the same binning and ntimes much greater
3525  //than the number of bins we can use a fast method
3526  Int_t first = fXaxis.GetFirst();
3527  Int_t last = fXaxis.GetLast();
3528  Int_t nbins = last-first+1;
3529  if (ntimes > 10*nbins) {
3530  try {
3531  CheckConsistency(this,h);
3532  Double_t sumw = h->Integral(first,last);
3533  if (sumw == 0) return;
3534  Double_t sumgen = 0;
3535  for (Int_t bin=first;bin<=last;bin++) {
3536  Double_t mean = h->RetrieveBinContent(bin)*ntimes/sumw;
3537  Double_t cont = (Double_t)gRandom->Poisson(mean);
3538  sumgen += cont;
3539  AddBinContent(bin,cont);
3540  if (fSumw2.fN) fSumw2.fArray[bin] += cont;
3541  }
3542 
3543  // fix for the fluctuations in the total number n
3544  // since we use Poisson instead of multinomial
3545  // add a correction to have ntimes as generated entries
3546  Int_t i;
3547  if (sumgen < ntimes) {
3548  // add missing entries
3549  for (i = Int_t(sumgen+0.5); i < ntimes; ++i)
3550  {
3551  Double_t x = h->GetRandom();
3552  Fill(x);
3553  }
3554  }
3555  else if (sumgen > ntimes) {
3556  // remove extra entries
3557  i = Int_t(sumgen+0.5);
3558  while( i > ntimes) {
3559  Double_t x = h->GetRandom();
3560  Int_t ibin = fXaxis.FindBin(x);
3561  Double_t y = RetrieveBinContent(ibin);
3562  // skip in case bin is empty
3563  if (y > 0) {
3564  SetBinContent(ibin, y-1.);
3565  i--;
3566  }
3567  }
3568  }
3569 
3570  ResetStats();
3571  return;
3572  }
3573  catch(std::exception&) {} // do nothing
3574  }
3575  // case of different axis and not too large ntimes
3576 
3577  if (h->ComputeIntegral() ==0) return;
3578  Int_t loop;
3579  Double_t x;
3580  for (loop=0;loop<ntimes;loop++) {
3581  x = h->GetRandom();
3582  Fill(x);
3583  }
3584 }
3585 
3586 ////////////////////////////////////////////////////////////////////////////////
3587 /// Return Global bin number corresponding to x,y,z
3588 ///
3589 /// 2-D and 3-D histograms are represented with a one dimensional
3590 /// structure. This has the advantage that all existing functions, such as
3591 /// GetBinContent, GetBinError, GetBinFunction work for all dimensions.
3592 /// This function tries to extend the axis if the given point belongs to an
3593 /// under-/overflow bin AND if CanExtendAllAxes() is true.
3594 ///
3595 /// See also TH1::GetBin, TAxis::FindBin and TAxis::FindFixBin
3596 
3598 {
3599  if (GetDimension() < 2) {
3600  return fXaxis.FindBin(x);
3601  }
3602  if (GetDimension() < 3) {
3603  Int_t nx = fXaxis.GetNbins()+2;
3604  Int_t binx = fXaxis.FindBin(x);
3605  Int_t biny = fYaxis.FindBin(y);
3606  return binx + nx*biny;
3607  }
3608  if (GetDimension() < 4) {
3609  Int_t nx = fXaxis.GetNbins()+2;
3610  Int_t ny = fYaxis.GetNbins()+2;
3611  Int_t binx = fXaxis.FindBin(x);
3612  Int_t biny = fYaxis.FindBin(y);
3613  Int_t binz = fZaxis.FindBin(z);
3614  return binx + nx*(biny +ny*binz);
3615  }
3616  return -1;
3617 }
3618 
3619 ////////////////////////////////////////////////////////////////////////////////
3620 /// Return Global bin number corresponding to x,y,z.
3621 ///
3622 /// 2-D and 3-D histograms are represented with a one dimensional
3623 /// structure. This has the advantage that all existing functions, such as
3624 /// GetBinContent, GetBinError, GetBinFunction work for all dimensions.
3625 /// This function DOES NOT try to extend the axis if the given point belongs
3626 /// to an under-/overflow bin.
3627 ///
3628 /// See also TH1::GetBin, TAxis::FindBin and TAxis::FindFixBin
3629 
3631 {
3632  if (GetDimension() < 2) {
3633  return fXaxis.FindFixBin(x);
3634  }
3635  if (GetDimension() < 3) {
3636  Int_t nx = fXaxis.GetNbins()+2;
3637  Int_t binx = fXaxis.FindFixBin(x);
3638  Int_t biny = fYaxis.FindFixBin(y);
3639  return binx + nx*biny;
3640  }
3641  if (GetDimension() < 4) {
3642  Int_t nx = fXaxis.GetNbins()+2;
3643  Int_t ny = fYaxis.GetNbins()+2;
3644  Int_t binx = fXaxis.FindFixBin(x);
3645  Int_t biny = fYaxis.FindFixBin(y);
3646  Int_t binz = fZaxis.FindFixBin(z);
3647  return binx + nx*(biny +ny*binz);
3648  }
3649  return -1;
3650 }
3651 
3652 ////////////////////////////////////////////////////////////////////////////////
3653 /// Find first bin with content > threshold for axis (1=x, 2=y, 3=z)
3654 /// if no bins with content > threshold is found the function returns -1.
3655 /// The search will occur between the specified first and last bin. Specifying
3656 /// the value of the last bin to search to less than zero will search until the
3657 /// last defined bin.
3658 
3659 Int_t TH1::FindFirstBinAbove(Double_t threshold, Int_t axis, Int_t firstBin, Int_t lastBin) const
3660 {
3661  if (fBuffer) ((TH1*)this)->BufferEmpty();
3662 
3663  if (axis < 1 || (axis > 1 && GetDimension() == 1 ) ||
3664  ( axis > 2 && GetDimension() == 2 ) || ( axis > 3 && GetDimension() > 3 ) ) {
3665  Warning("FindFirstBinAbove","Invalid axis number : %d, axis x assumed\n",axis);
3666  axis = 1;
3667  }
3668  if (firstBin < 1) {
3669  firstBin = 1;
3670  }
3671  Int_t nbinsx = fXaxis.GetNbins();
3672  Int_t nbinsy = (GetDimension() > 1 ) ? fYaxis.GetNbins() : 1;
3673  Int_t nbinsz = (GetDimension() > 2 ) ? fZaxis.GetNbins() : 1;
3674 
3675  if (axis == 1) {
3676  if (lastBin < 0 || lastBin > fXaxis.GetNbins()) {
3677  lastBin = fXaxis.GetNbins();
3678  }
3679  for (Int_t binx = firstBin; binx <= lastBin; binx++) {
3680  for (Int_t biny = 1; biny <= nbinsy; biny++) {
3681  for (Int_t binz = 1; binz <= nbinsz; binz++) {
3682  if (RetrieveBinContent(GetBin(binx,biny,binz)) > threshold) return binx;
3683  }
3684  }
3685  }
3686  }
3687  else if (axis == 2) {
3688  if (lastBin < 0 || lastBin > fYaxis.GetNbins()) {
3689  lastBin = fYaxis.GetNbins();
3690  }
3691  for (Int_t biny = firstBin; biny <= lastBin; biny++) {
3692  for (Int_t binx = 1; binx <= nbinsx; binx++) {
3693  for (Int_t binz = 1; binz <= nbinsz; binz++) {
3694  if (RetrieveBinContent(GetBin(binx,biny,binz)) > threshold) return biny;
3695  }
3696  }
3697  }
3698  }
3699  else if (axis == 3) {
3700  if (lastBin < 0 || lastBin > fZaxis.GetNbins()) {
3701  lastBin = fZaxis.GetNbins();
3702  }
3703  for (Int_t binz = firstBin; binz <= lastBin; binz++) {
3704  for (Int_t binx = 1; binx <= nbinsx; binx++) {
3705  for (Int_t biny = 1; biny <= nbinsy; biny++) {
3706  if (RetrieveBinContent(GetBin(binx,biny,binz)) > threshold) return binz;
3707  }
3708  }
3709  }
3710  }
3711 
3712  return -1;
3713 }
3714 
3715 ////////////////////////////////////////////////////////////////////////////////
3716 /// Find last bin with content > threshold for axis (1=x, 2=y, 3=z)
3717 /// if no bins with content > threshold is found the function returns -1.
3718 /// The search will occur between the specified first and last bin. Specifying
3719 /// the value of the last bin to search to less than zero will search until the
3720 /// last defined bin.
3721 
3722 Int_t TH1::FindLastBinAbove(Double_t threshold, Int_t axis, Int_t firstBin, Int_t lastBin) const
3723 {
3724  if (fBuffer) ((TH1*)this)->BufferEmpty();
3725 
3726 
3727  if (axis < 1 || ( axis > 1 && GetDimension() == 1 ) ||
3728  ( axis > 2 && GetDimension() == 2 ) || ( axis > 3 && GetDimension() > 3) ) {
3729  Warning("FindFirstBinAbove","Invalid axis number : %d, axis x assumed\n",axis);
3730  axis = 1;
3731  }
3732  if (firstBin < 1) {
3733  firstBin = 1;
3734  }
3735  Int_t nbinsx = fXaxis.GetNbins();
3736  Int_t nbinsy = (GetDimension() > 1 ) ? fYaxis.GetNbins() : 1;
3737  Int_t nbinsz = (GetDimension() > 2 ) ? fZaxis.GetNbins() : 1;
3738 
3739  if (axis == 1) {
3740  if (lastBin < 0 || lastBin > fXaxis.GetNbins()) {
3741  lastBin = fXaxis.GetNbins();
3742  }
3743  for (Int_t binx = lastBin; binx >= firstBin; binx--) {
3744  for (Int_t biny = 1; biny <= nbinsy; biny++) {
3745  for (Int_t binz = 1; binz <= nbinsz; binz++) {
3746  if (RetrieveBinContent(GetBin(binx, biny, binz)) > threshold) return binx;
3747  }
3748  }
3749  }
3750  }
3751  else if (axis == 2) {
3752  if (lastBin < 0 || lastBin > fYaxis.GetNbins()) {
3753  lastBin = fYaxis.GetNbins();
3754  }
3755  for (Int_t biny = lastBin; biny >= firstBin; biny--) {
3756  for (Int_t binx = 1; binx <= nbinsx; binx++) {
3757  for (Int_t binz = 1; binz <= nbinsz; binz++) {
3758  if (RetrieveBinContent(GetBin(binx, biny, binz)) > threshold) return biny;
3759  }
3760  }
3761  }
3762  }
3763  else if (axis == 3) {
3764  if (lastBin < 0 || lastBin > fZaxis.GetNbins()) {
3765  lastBin = fZaxis.GetNbins();
3766  }
3767  for (Int_t binz = lastBin; binz >= firstBin; binz--) {
3768  for (Int_t binx = 1; binx <= nbinsx; binx++) {
3769  for (Int_t biny = 1; biny <= nbinsy; biny++) {
3770  if (RetrieveBinContent(GetBin(binx, biny, binz)) > threshold) return binz;
3771  }
3772  }
3773  }
3774  }
3775 
3776  return -1;
3777 }
3778 
3779 ////////////////////////////////////////////////////////////////////////////////
3780 /// Search object named name in the list of functions.
3781 
3782 TObject *TH1::FindObject(const char *name) const
3783 {
3784  if (fFunctions) return fFunctions->FindObject(name);
3785  return 0;
3786 }
3787 
3788 ////////////////////////////////////////////////////////////////////////////////
3789 /// Search object obj in the list of functions.
3790 
3791 TObject *TH1::FindObject(const TObject *obj) const
3792 {
3793  if (fFunctions) return fFunctions->FindObject(obj);
3794  return 0;
3795 }
3796 
3797 ////////////////////////////////////////////////////////////////////////////////
3798 /// Fit histogram with function fname.
3799 ///
3800 /// fname is the name of an already predefined function created by TF1 or TF2
3801 /// Predefined functions such as gaus, expo and poln are automatically
3802 /// created by ROOT.
3803 /// fname can also be a formula, accepted by the linear fitter (linear parts divided
3804 /// by "++" sign), for example "x++sin(x)" for fitting "[0]*x+[1]*sin(x)"
3805 ///
3806 /// This function finds a pointer to the TF1 object with name fname
3807 /// and calls TH1::Fit(TF1 *f1,...)
3808 
3809 TFitResultPtr TH1::Fit(const char *fname ,Option_t *option ,Option_t *goption, Double_t xxmin, Double_t xxmax)
3810 {
3811  char *linear;
3812  linear= (char*)strstr(fname, "++");
3813  Int_t ndim=GetDimension();
3814  if (linear){
3815  if (ndim<2){
3816  TF1 f1(fname, fname, xxmin, xxmax);
3817  return Fit(&f1,option,goption,xxmin,xxmax);
3818  }
3819  else if (ndim<3){
3820  TF2 f2(fname, fname);
3821  return Fit(&f2,option,goption,xxmin,xxmax);
3822  }
3823  else{
3824  TF3 f3(fname, fname);
3825  return Fit(&f3,option,goption,xxmin,xxmax);
3826  }
3827  }
3828  else{
3829  TF1 * f1 = (TF1*)gROOT->GetFunction(fname);
3830  if (!f1) { Printf("Unknown function: %s",fname); return -1; }
3831  return Fit(f1,option,goption,xxmin,xxmax);
3832  }
3833 }
3834 
3835 ////////////////////////////////////////////////////////////////////////////////
3836 /// Fit histogram with function f1.
3837 ///
3838 /// \param[in] option fit options is given in parameter option.
3839 /// - "W" Ignore the bin uncertainties when fitting using the default least square (chi2) method but skip empty bins
3840 /// - "WW" Ignore the bin uncertainties when fitting using the default least square (chi2) method and include also the empty bins
3841 /// - "I" Use integral of function in bin, normalized by the bin volume,
3842 /// instead of value at bin center
3843 /// - "L" Use Loglikelihood method (default is chisquare method)
3844 /// - "WL" Use Loglikelihood method and bin contents are not integer,
3845 /// i.e. histogram is weighted (must have Sumw2() set)
3846 /// -"MULTI" Use Loglikelihood method based on multi-nomial distribution.
3847 /// In this case function must be normalized and one fits only the function shape (a not extended binned
3848 /// likelihood fit)
3849 /// - "P" Use Pearson chi2 (using expected errors instead of observed errors)
3850 /// - "U" Use a User specified fitting algorithm (via SetFCN)
3851 /// - "Q" Quiet mode (minimum printing)
3852 /// - "V" Verbose mode (default is between Q and V)
3853 /// - "E" Perform better Errors estimation using Minos technique
3854 /// - "B" User defined parameter settings are used for predefined functions
3855 /// like "gaus", "expo", "poln", "landau".
3856 /// Use this option when you want to fix one or more parameters for these functions.
3857 /// - "M" More. Improve fit results.
3858 /// It uses the IMPROVE command of TMinuit (see TMinuit::mnimpr).
3859 /// This algorithm attempts to improve the found local minimum by searching for a
3860 /// better one.
3861 /// - "R" Use the Range specified in the function range
3862 /// - "N" Do not store the graphics function, do not draw
3863 /// - "0" Do not plot the result of the fit. By default the fitted function
3864 /// is drawn unless the option"N" above is specified.
3865 /// - "+" Add this new fitted function to the list of fitted functions
3866 /// (by default, any previous function is deleted)
3867 /// - "C" In case of linear fitting, don't calculate the chisquare
3868 /// (saves time)
3869 /// - "F" If fitting a polN, switch to minuit fitter
3870 /// - "S" The result of the fit is returned in the TFitResultPtr
3871 /// (see below Access to the Fit Result)
3872 /// \param[in] goption specify a list of graphics options. See TH1::Draw for a complete list of these options.
3873 /// \param[in] xxmin range
3874 /// \param[in] xxmax range
3875 ///
3876 /// In order to use the Range option, one must first create a function
3877 /// with the expression to be fitted. For example, if your histogram
3878 /// has a defined range between -4 and 4 and you want to fit a gaussian
3879 /// only in the interval 1 to 3, you can do:
3880 ///
3881 /// ~~~ {.cpp}
3882 /// TF1 *f1 = new TF1("f1", "gaus", 1, 3);
3883 /// histo->Fit("f1", "R");
3884 /// ~~~
3885 ///
3886 /// ## Setting initial conditions
3887 /// Parameters must be initialized before invoking the Fit function.
3888 /// The setting of the parameter initial values is automatic for the
3889 /// predefined functions : poln, expo, gaus, landau. One can however disable
3890 /// this automatic computation by specifying the option "B".
3891 /// Note that if a predefined function is defined with an argument,
3892 /// eg, gaus(0), expo(1), you must specify the initial values for
3893 /// the parameters.
3894 /// You can specify boundary limits for some or all parameters via
3895 ///
3896 /// ~~~ {.cpp}
3897 /// f1->SetParLimits(p_number, parmin, parmax);
3898 /// ~~~
3899 ///
3900 /// if parmin>=parmax, the parameter is fixed
3901 /// Note that you are not forced to fix the limits for all parameters.
3902 /// For example, if you fit a function with 6 parameters, you can do:
3903 ///
3904 /// ~~~ {.cpp}
3905 /// func->SetParameters(0, 3.1, 1.e-6, -8, 0, 100);
3906 /// func->SetParLimits(3, -10, -4);
3907 /// func->FixParameter(4, 0);
3908 /// func->SetParLimits(5, 1, 1);
3909 /// ~~~
3910 ///
3911 /// With this setup, parameters 0->2 can vary freely
3912 /// Parameter 3 has boundaries [-10,-4] with initial value -8
3913 /// Parameter 4 is fixed to 0
3914 /// Parameter 5 is fixed to 100.
3915 /// When the lower limit and upper limit are equal, the parameter is fixed.
3916 /// However to fix a parameter to 0, one must call the FixParameter function.
3917 ///
3918 ///
3919 /// #### Changing the fitting objective function
3920 ///
3921 /// By default a chi square function is used for fitting. When option "L" (or "LL") is used
3922 /// a Poisson likelihood function (see note below) is used.
3923 /// Using option "MULTI" a multinomial likelihood fit is used. In this case the function normalization is not fitted
3924 /// but only the function shape. Therefore the provided function must be normalized.
3925 /// The functions are defined in the header Fit/Chi2Func.h or Fit/PoissonLikelihoodFCN and they
3926 /// are implemented using the routines FitUtil::EvaluateChi2 or FitUtil::EvaluatePoissonLogL in
3927 /// the file math/mathcore/src/FitUtil.cxx.
3928 /// To specify a User defined fitting function, specify option "U" and
3929 /// call the following functions:
3930 ///
3931 /// ~~~ {.cpp}
3932 /// TVirtualFitter::Fitter(myhist)->SetFCN(MyFittingFunction)
3933 /// ~~~
3934 ///
3935 /// where MyFittingFunction is of type:
3936 ///
3937 /// ~~~ {.cpp}
3938 /// extern void MyFittingFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
3939 /// ~~~
3940 ///
3941 /// #### Chi2 Fits
3942 ///
3943 /// By default a chi2 (least-square) fit is performed on the histogram. The so-called modified least-square method
3944 /// is used where the residual for each bin is computed using as error the observed value (the bin error)
3945 ///
3946 /// \f[
3947 /// Chi2 = \sum{ \left(\frac{y(i) - f(x(i) | p )}{e(i)} \right)^2 }
3948 /// \f]
3949 ///
3950 /// where y(i) is the bin content for each bin i, x(i) is the bin center and e(i) is the bin error (sqrt(y(i) for
3951 /// an un-weighted histogram. Bins with zero errors are excluded from the fit. See also later the note on the treatment
3952 /// of empty bins. When using option "I" the residual is computed not using the function value at the bin center, f
3953 /// (x(i) | p), but the integral of the function in the bin, Integral{ f(x|p)dx } divided by the bin volume
3954 ///
3955 /// #### Likelihood Fits
3956 ///
3957 /// When using option "L" a likelihood fit is used instead of the default chi2 square fit.
3958 /// The likelihood is built assuming a Poisson probability density function for each bin.
3959 /// The negative log-likelihood to be minimized is
3960 ///
3961 /// \f[
3962 /// NLL = \sum{ log Poisson ( y(i) | f(x(i) | p ) ) }
3963 /// \f]
3964 ///
3965 /// The exact likelihood used is the Poisson likelihood described in this paper:
3966 /// S. Baker and R. D. Cousins, “Clarification of the use of chi-square and likelihood functions in fits to histograms,”
3967 /// Nucl. Instrum. Meth. 221 (1984) 437.
3968 ///
3969 /// This method can then be used only when the bin content represents counts (i.e. errors are sqrt(N) ).
3970 /// The likelihood method has the advantage of treating correctly bins with low statistics. In case of high
3971 /// statistics/bin the distribution of the bin content becomes a normal distribution and the likelihood and chi2 fit
3972 /// give the same result.
3973 ///
3974 /// The likelihood method, although a bit slower, it is therefore the recommended method in case of low
3975 /// bin statistics, where the chi2 method may give incorrect results, in particular when there are
3976 /// several empty bins (see also below).
3977 /// In case of a weighted histogram, it is possible to perform a likelihood fit by using the
3978 /// option "WL". Note a weighted histogram is an histogram which has been filled with weights and it
3979 /// contains the sum of the weight square ( TH1::Sumw2() has been called). The bin error for a weighted
3980 /// histogram is the square root of the sum of the weight square.
3981 ///
3982 /// #### Treatment of Empty Bins
3983 ///
3984 /// Empty bins, which have the content equal to zero AND error equal to zero,
3985 /// are excluded by default from the chisquare fit, but they are considered in the likelihood fit.
3986 /// since they affect the likelihood if the function value in these bins is not negligible.
3987 /// When using option "WW" these bins will be considered in the chi2 fit with an error of 1.
3988 /// Note that if the histogram is having bins with zero content and non zero-errors they are considered as
3989 /// any other bins in the fit. Instead bins with zero error and non-zero content are excluded in the chi2 fit.
3990 /// A likelihood fit should also not be performed on such an histogram, since we are assuming a wrong pdf for each bin.
3991 /// In general, one should not fit an histogram with non-empty bins and zero errors, apart if all the bins have zero
3992 /// errors. In this case one could use the option "w", which gives a weight=1 for each bin (unweighted least-square
3993 /// fit).
3994 /// Note that in case of histogram with no errors (chi2 fit with option W or W1) the resulting fitted parameter errors
3995 /// are corrected by the obtained chi2 value using this expression: errorp *= sqrt(chisquare/(ndf-1))
3996 ///
3997 /// #### Fitting a histogram of dimension N with a function of dimension N-1
3998 ///
3999 /// It is possible to fit a TH2 with a TF1 or a TH3 with a TF2.
4000 /// In this case the option "Integral" is not allowed and each cell has
4001 /// equal weight. Also in this case th eobtained parameter error are corrected as in the case when the
4002 /// option "W" is used (see above)
4003 ///
4004 /// #### Associated functions
4005 ///
4006 /// One or more object (typically a TF1*) can be added to the list
4007 /// of functions (fFunctions) associated to each histogram.
4008 /// When TH1::Fit is invoked, the fitted function is added to this list.
4009 /// Given an histogram h, one can retrieve an associated function
4010 /// with:
4011 ///
4012 /// ~~~ {.cpp}
4013 /// TF1 *myfunc = h->GetFunction("myfunc");
4014 /// ~~~
4015 ///
4016 /// #### Access to the fit result
4017 ///
4018 /// The function returns a TFitResultPtr which can hold a pointer to a TFitResult object.
4019 /// By default the TFitResultPtr contains only the status of the fit which is return by an
4020 /// automatic conversion of the TFitResultPtr to an integer. One can write in this case directly:
4021 ///
4022 /// ~~~ {.cpp}
4023 /// Int_t fitStatus = h->Fit(myFunc)
4024 /// ~~~
4025 ///
4026 /// If the option "S" is instead used, TFitResultPtr contains the TFitResult and behaves as a smart
4027 /// pointer to it. For example one can do:
4028 ///
4029 /// ~~~ {.cpp}
4030 /// TFitResultPtr r = h->Fit(myFunc,"S");
4031 /// TMatrixDSym cov = r->GetCovarianceMatrix(); // to access the covariance matrix
4032 /// Double_t chi2 = r->Chi2(); // to retrieve the fit chi2
4033 /// Double_t par0 = r->Parameter(0); // retrieve the value for the parameter 0
4034 /// Double_t err0 = r->ParError(0); // retrieve the error for the parameter 0
4035 /// r->Print("V"); // print full information of fit including covariance matrix
4036 /// r->Write(); // store the result in a file
4037 /// ~~~
4038 ///
4039 /// The fit parameters, error and chi2 (but not covariance matrix) can be retrieved also
4040 /// from the fitted function.
4041 /// If the histogram is made persistent, the list of
4042 /// associated functions is also persistent. Given a pointer (see above)
4043 /// to an associated function myfunc, one can retrieve the function/fit
4044 /// parameters with calls such as:
4045 ///
4046 /// ~~~ {.cpp}
4047 /// Double_t chi2 = myfunc->GetChisquare();
4048 /// Double_t par0 = myfunc->GetParameter(0); //value of 1st parameter
4049 /// Double_t err0 = myfunc->GetParError(0); //error on first parameter
4050 /// ~~~
4051 ///
4052 /// #### Access to the fit status
4053 ///
4054 /// The status of the fit can be obtained converting the TFitResultPtr to an integer
4055 /// independently if the fit option "S" is used or not:
4056 ///
4057 /// ~~~ {.cpp}
4058 /// TFitResultPtr r = h->Fit(myFunc,opt);
4059 /// Int_t fitStatus = r;
4060 /// ~~~
4061 ///
4062 /// The fitStatus is 0 if the fit is OK (i.e no error occurred).
4063 /// The value of the fit status code is negative in case of an error not connected with the
4064 /// minimization procedure, for example when a wrong function is used.
4065 /// Otherwise the return value is the one returned from the minimization procedure.
4066 /// When TMinuit (default case) or Minuit2 are used as minimizer the status returned is :
4067 /// `fitStatus = migradResult + 10*minosResult + 100*hesseResult + 1000*improveResult`.
4068 /// TMinuit will return 0 (for migrad, minos, hesse or improve) in case of success and 4 in
4069 /// case of error (see the documentation of TMinuit::mnexcm). So for example, for an error
4070 /// only in Minos but not in Migrad a fitStatus of 40 will be returned.
4071 /// Minuit2 will return also 0 in case of success and different values in migrad minos or
4072 /// hesse depending on the error. See in this case the documentation of
4073 /// Minuit2Minimizer::Minimize for the migradResult, Minuit2Minimizer::GetMinosError for the
4074 /// minosResult and Minuit2Minimizer::Hesse for the hesseResult.
4075 /// If other minimizers are used see their specific documentation for the status code returned.
4076 /// For example in the case of Fumili, for the status returned see TFumili::Minimize.
4077 ///
4078 /// #### Excluding points
4079 ///
4080 /// Use TF1::RejectPoint inside your fitting function to exclude points
4081 /// within a certain range from the fit. Example:
4082 ///
4083 /// ~~~ {.cpp}
4084 /// Double_t fline(Double_t *x, Double_t *par)
4085 /// {
4086 /// if (x[0] > 2.5 && x[0] < 3.5) {
4087 /// TF1::RejectPoint();
4088 /// return 0;
4089 /// }
4090 /// return par[0] + par[1]*x[0];
4091 /// }
4092 ///
4093 /// void exclude() {
4094 /// TF1 *f1 = new TF1("f1", "[0] +[1]*x +gaus(2)", 0, 5);
4095 /// f1->SetParameters(6, -1,5, 3, 0.2);
4096 /// TH1F *h = new TH1F("h", "background + signal", 100, 0, 5);
4097 /// h->FillRandom("f1", 2000);
4098 /// TF1 *fline = new TF1("fline", fline, 0, 5, 2);
4099 /// fline->SetParameters(2, -1);
4100 /// h->Fit("fline", "l");
4101 /// }
4102 /// ~~~
4103 ///
4104 /// #### Warning when using the option "0"
4105 ///
4106 /// When selecting the option "0", the fitted function is added to
4107 /// the list of functions of the histogram, but it is not drawn.
4108 /// You can undo what you disabled in the following way:
4109 ///
4110 /// ~~~ {.cpp}
4111 /// h.Fit("myFunction", "0"); // fit, store function but do not draw
4112 /// h.Draw(); function is not drawn
4113 /// const Int_t kNotDraw = 1<<9;
4114 /// h.GetFunction("myFunction")->ResetBit(kNotDraw);
4115 /// h.Draw(); // function is visible again
4116 /// ~~~
4117 ///
4118 /// #### Access to the Minimizer information during fitting
4119 ///
4120 /// This function calls, the ROOT::Fit::FitObject function implemented in HFitImpl.cxx
4121 /// which uses the ROOT::Fit::Fitter class. The Fitter class creates the objective function
4122 /// (e.g. chi2 or likelihood) and uses an implementation of the Minimizer interface for minimizing
4123 /// the function.
4124 /// The default minimizer is Minuit (class TMinuitMinimizer which calls TMinuit).
4125 /// The default can be set in the resource file in etc/system.rootrc. For example
4126 ///
4127 /// ~~~ {.cpp}
4128 /// Root.Fitter: Minuit2
4129 /// ~~~
4130 ///
4131 /// A different fitter can also be set via ROOT::Math::MinimizerOptions::SetDefaultMinimizer
4132 /// (or TVirtualFitter::SetDefaultFitter).
4133 /// For example ROOT::Math::MinimizerOptions::SetDefaultMinimizer("GSLMultiMin","BFGS");
4134 /// will set the usage of the BFGS algorithm of the GSL multi-dimensional minimization
4135 /// (implemented in libMathMore). ROOT::Math::MinimizerOptions can be used also to set other
4136 /// default options, like maximum number of function calls, minimization tolerance or print
4137 /// level. See the documentation of this class.
4138 ///
4139 /// For fitting linear functions (containing the "++" sign" and polN functions,
4140 /// the linear fitter is automatically initialized.
4141 
4142 TFitResultPtr TH1::Fit(TF1 *f1 ,Option_t *option ,Option_t *goption, Double_t xxmin, Double_t xxmax)
4143 {
4144  // implementation of Fit method is in file hist/src/HFitImpl.cxx
4145  Foption_t fitOption;
4147 
4148  // create range and minimizer options with default values
4149  ROOT::Fit::DataRange range(xxmin,xxmax);
4150  ROOT::Math::MinimizerOptions minOption;
4151 
4152  // need to empty the buffer before
4153  // (t.b.d. do a ML unbinned fit with buffer data)
4154  if (fBuffer) BufferEmpty();
4155 
4156  return ROOT::Fit::FitObject(this, f1 , fitOption , minOption, goption, range);
4157 }
4158 
4159 ////////////////////////////////////////////////////////////////////////////////
4160 /// Display a panel with all histogram fit options.
4161 ///
4162 /// See class TFitPanel for example
4163 
4164 void TH1::FitPanel()
4165 {
4166  if (!gPad)
4167  gROOT->MakeDefCanvas();
4168 
4169  if (!gPad) {
4170  Error("FitPanel", "Unable to create a default canvas");
4171  return;
4172  }
4173 
4174 
4175  // use plugin manager to create instance of TFitEditor
4176  TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TFitEditor");
4177  if (handler && handler->LoadPlugin() != -1) {
4178  if (handler->ExecPlugin(2, gPad, this) == 0)
4179  Error("FitPanel", "Unable to create the FitPanel");
4180  }
4181  else
4182  Error("FitPanel", "Unable to find the FitPanel plug-in");
4183 }
4184 
4185 ////////////////////////////////////////////////////////////////////////////////
4186 /// Return an histogram containing the asymmetry of this histogram with h2,
4187 /// where the asymmetry is defined as:
4188 ///
4189 /// ~~~ {.cpp}
4190 /// Asymmetry = (h1 - h2)/(h1 + h2) where h1 = this
4191 /// ~~~
4192 ///
4193 /// works for 1D, 2D, etc. histograms
4194 /// c2 is an optional argument that gives a relative weight between the two
4195 /// histograms, and dc2 is the error on this weight. This is useful, for example,
4196 /// when forming an asymmetry between two histograms from 2 different data sets that
4197 /// need to be normalized to each other in some way. The function calculates
4198 /// the errors assuming Poisson statistics on h1 and h2 (that is, dh = sqrt(h)).
4199 ///
4200 /// example: assuming 'h1' and 'h2' are already filled
4201 ///
4202 /// ~~~ {.cpp}
4203 /// h3 = h1->GetAsymmetry(h2)
4204 /// ~~~
4205 ///
4206 /// then 'h3' is created and filled with the asymmetry between 'h1' and 'h2';
4207 /// h1 and h2 are left intact.
4208 ///
4209 /// Note that it is the user's responsibility to manage the created histogram.
4210 /// The name of the returned histogram will be `Asymmetry_nameOfh1-nameOfh2`
4211 ///
4212 /// code proposed by Jason Seely (seely@mit.edu) and adapted by R.Brun
4213 ///
4214 /// clone the histograms so top and bottom will have the
4215 /// correct dimensions:
4216 /// Sumw2 just makes sure the errors will be computed properly
4217 /// when we form sums and ratios below.
4218 
4220 {
4221  TH1 *h1 = this;
4222  TString name = TString::Format("Asymmetry_%s-%s",h1->GetName(),h2->GetName() );
4223  TH1 *asym = (TH1*)Clone(name);
4224 
4225  // set also the title
4226  TString title = TString::Format("(%s - %s)/(%s+%s)",h1->GetName(),h2->GetName(),h1->GetName(),h2->GetName() );
4227  asym->SetTitle(title);
4228 
4229  asym->Sumw2();
4230  Bool_t addStatus = TH1::AddDirectoryStatus();
4232  TH1 *top = (TH1*)asym->Clone();
4233  TH1 *bottom = (TH1*)asym->Clone();
4234  TH1::AddDirectory(addStatus);
4235 
4236  // form the top and bottom of the asymmetry, and then divide:
4237  top->Add(h1,h2,1,-c2);
4238  bottom->Add(h1,h2,1,c2);
4239  asym->Divide(top,bottom);
4240 
4241  Int_t xmax = asym->GetNbinsX();
4242  Int_t ymax = asym->GetNbinsY();
4243  Int_t zmax = asym->GetNbinsZ();
4244 
4245  if (h1->fBuffer) h1->BufferEmpty(1);
4246  if (h2->fBuffer) h2->BufferEmpty(1);
4247  if (bottom->fBuffer) bottom->BufferEmpty(1);
4248 
4249  // now loop over bins to calculate the correct errors
4250  // the reason this error calculation looks complex is because of c2
4251  for(Int_t i=1; i<= xmax; i++){
4252  for(Int_t j=1; j<= ymax; j++){
4253  for(Int_t k=1; k<= zmax; k++){
4254  Int_t bin = GetBin(i, j, k);
4255  // here some bin contents are written into variables to make the error
4256  // calculation a little more legible:
4257  Double_t a = h1->RetrieveBinContent(bin);
4258  Double_t b = h2->RetrieveBinContent(bin);
4259  Double_t bot = bottom->RetrieveBinContent(bin);
4260 
4261  // make sure there are some events, if not, then the errors are set = 0
4262  // automatically.
4263  //if(bot < 1){} was changed to the next line from recommendation of Jason Seely (28 Nov 2005)
4264  if(bot < 1e-6){}
4265  else{
4266  // computation of errors by Christos Leonidopoulos
4267  Double_t dasq = h1->GetBinErrorSqUnchecked(bin);
4268  Double_t dbsq = h2->GetBinErrorSqUnchecked(bin);
4269  Double_t error = 2*TMath::Sqrt(a*a*c2*c2*dbsq + c2*c2*b*b*dasq+a*a*b*b*dc2*dc2)/(bot*bot);
4270  asym->SetBinError(i,j,k,error);
4271  }
4272  }
4273  }
4274  }
4275  delete top;
4276  delete bottom;
4277 
4278  return asym;
4279 }
4280 
4281 ////////////////////////////////////////////////////////////////////////////////
4282 /// Static function
4283 /// return the default buffer size for automatic histograms
4284 /// the parameter fgBufferSize may be changed via SetDefaultBufferSize
4285 
4287 {
4288  return fgBufferSize;
4289 }
4290 
4291 ////////////////////////////////////////////////////////////////////////////////
4292 /// Return kTRUE if TH1::Sumw2 must be called when creating new histograms.
4293 /// see TH1::SetDefaultSumw2.
4294 
4296 {
4297  return fgDefaultSumw2;
4298 }
4299 
4300 ////////////////////////////////////////////////////////////////////////////////
4301 /// Return the current number of entries.
4302 
4303 Double_t TH1::GetEntries() const
4304 {
4305  if (fBuffer) {
4306  Int_t nentries = (Int_t) fBuffer[0];
4307  if (nentries > 0) return nentries;
4308  }
4309 
4310  return fEntries;
4311 }
4312 
4313 ////////////////////////////////////////////////////////////////////////////////
4314 /// Number of effective entries of the histogram.
4315 ///
4316 /// \f[
4317 /// neff = \frac{(\sum Weights )^2}{(\sum Weight^2 )}
4318 /// \f]
4319 ///
4320 /// In case of an unweighted histogram this number is equivalent to the
4321 /// number of entries of the histogram.
4322 /// For a weighted histogram, this number corresponds to the hypothetical number of unweighted entries
4323 /// a histogram would need to have the same statistical power as this weighted histogram.
4324 /// Note: The underflow/overflow are included if one has set the TH1::StatOverFlows flag
4325 /// and if the statistics has been computed at filling time.
4326 /// If a range is set in the histogram the number is computed from the given range.
4327 
4329 {
4330  Stat_t s[kNstat];
4331  this->GetStats(s);// s[1] sum of squares of weights, s[0] sum of weights
4332  return (s[1] ? s[0]*s[0]/s[1] : TMath::Abs(s[0]) );
4333 }
4334 
4335 ////////////////////////////////////////////////////////////////////////////////
4336 /// Set highlight (enable/disable) mode for the histogram
4337 /// by default highlight mode is disable
4338 
4339 void TH1::SetHighlight(Bool_t set)
4340 {
4341  if (IsHighlight() == set) return;
4342  if (fDimension > 2) {
4343  Info("SetHighlight", "Supported only 1-D or 2-D histograms");
4344  return;
4345  }
4346 
4347  if (!fPainter) {
4348  Info("SetHighlight", "Need to draw histogram first");
4349  return;
4350  }
4351  SetBit(kIsHighlight, set);
4353 }
4354 
4355 ////////////////////////////////////////////////////////////////////////////////
4356 /// Redefines TObject::GetObjectInfo.
4357 /// Displays the histogram info (bin number, contents, integral up to bin
4358 /// corresponding to cursor position px,py
4359 
4360 char *TH1::GetObjectInfo(Int_t px, Int_t py) const
4361 {
4362  return ((TH1*)this)->GetPainter()->GetObjectInfo(px,py);
4363 }
4364 
4365 ////////////////////////////////////////////////////////////////////////////////
4366 /// Return pointer to painter.
4367 /// If painter does not exist, it is created
4368 
4370 {
4371  if (!fPainter) {
4372  TString opt = option;
4373  opt.ToLower();
4374  if (opt.Contains("gl") || gStyle->GetCanvasPreferGL()) {
4375  //try to create TGLHistPainter
4376  TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TGLHistPainter");
4377 
4378  if (handler && handler->LoadPlugin() != -1)
4379  fPainter = reinterpret_cast<TVirtualHistPainter *>(handler->ExecPlugin(1, this));
4380  }
4381  }
4382 
4384 
4385  return fPainter;
4386 }
4387 
4388 ////////////////////////////////////////////////////////////////////////////////
4389 /// Compute Quantiles for this histogram
4390 /// Quantile x_q of a probability distribution Function F is defined as
4391 ///
4392 /// ~~~ {.cpp}
4393 /// F(x_q) = q with 0 <= q <= 1.
4394 /// ~~~
4395 ///
4396 /// For instance the median x_0.5 of a distribution is defined as that value
4397 /// of the random variable for which the distribution function equals 0.5:
4398 ///
4399 /// ~~~ {.cpp}
4400 /// F(x_0.5) = Probability(x < x_0.5) = 0.5
4401 /// ~~~
4402 ///
4403 /// code from Eddy Offermann, Renaissance
4404 ///
4405 /// \param[in] nprobSum maximum size of array q and size of array probSum (if given)
4406 /// \param[in] probSum array of positions where quantiles will be computed.
4407 /// - if probSum is null, probSum will be computed internally and will
4408 /// have a size = number of bins + 1 in h. it will correspond to the
4409 /// quantiles calculated at the lowest edge of the histogram (quantile=0) and
4410 /// all the upper edges of the bins.
4411 /// - if probSum is not null, it is assumed to contain at least nprobSum values.
4412 /// \param[out] q array q filled with nq quantiles
4413 /// \return value nq (<=nprobSum) with the number of quantiles computed
4414 ///
4415 /// Note that the Integral of the histogram is automatically recomputed
4416 /// if the number of entries is different of the number of entries when
4417 /// the integral was computed last time. In case you do not use the Fill
4418 /// functions to fill your histogram, but SetBinContent, you must call
4419 /// TH1::ComputeIntegral before calling this function.
4420 ///
4421 /// Getting quantiles q from two histograms and storing results in a TGraph,
4422 /// a so-called QQ-plot
4423 ///
4424 /// ~~~ {.cpp}
4425 /// TGraph *gr = new TGraph(nprob);
4426 /// h1->GetQuantiles(nprob,gr->GetX());
4427 /// h2->GetQuantiles(nprob,gr->GetY());
4428 /// gr->Draw("alp");
4429 /// ~~~
4430 ///
4431 /// Example:
4432 ///
4433 /// ~~~ {.cpp}
4434 /// void quantiles() {
4435 /// // demo for quantiles
4436 /// const Int_t nq = 20;
4437 /// TH1F *h = new TH1F("h","demo quantiles",100,-3,3);
4438 /// h->FillRandom("gaus",5000);
4439 ///
4440 /// Double_t xq[nq]; // position where to compute the quantiles in [0,1]
4441 /// Double_t yq[nq]; // array to contain the quantiles
4442 /// for (Int_t i=0;i<nq;i++) xq[i] = Float_t(i+1)/nq;
4443 /// h->GetQuantiles(nq,yq,xq);
4444 ///
4445 /// //show the original histogram in the top pad
4446 /// TCanvas *c1 = new TCanvas("c1","demo quantiles",10,10,700,900);
4447 /// c1->Divide(1,2);
4448 /// c1->cd(1);
4449 /// h->Draw();
4450 ///
4451 /// // show the quantiles in the bottom pad
4452 /// c1->cd(2);
4453 /// gPad->SetGrid();
4454 /// TGraph *gr = new TGraph(nq,xq,yq);
4455 /// gr->SetMarkerStyle(21);
4456 /// gr->Draw("alp");
4457 /// }
4458 /// ~~~
4459 
4460 Int_t TH1::GetQuantiles(Int_t nprobSum, Double_t *q, const Double_t *probSum)
4461 {
4462  if (GetDimension() > 1) {
4463  Error("GetQuantiles","Only available for 1-d histograms");
4464  return 0;
4465  }
4466 
4467  const Int_t nbins = GetXaxis()->GetNbins();
4468  if (!fIntegral) ComputeIntegral();
4469  if (fIntegral[nbins+1] != fEntries) ComputeIntegral();
4470 
4471  Int_t i, ibin;
4472  Double_t *prob = (Double_t*)probSum;
4473  Int_t nq = nprobSum;
4474  if (probSum == 0) {
4475  nq = nbins+1;
4476  prob = new Double_t[nq];
4477  prob[0] = 0;
4478  for (i=1;i<nq;i++) {
4479  prob[i] = fIntegral[i]/fIntegral[nbins];
4480  }
4481  }
4482 
4483  for (i = 0; i < nq; i++) {
4484  ibin = TMath::BinarySearch(nbins,fIntegral,prob[i]);
4485  while (ibin < nbins-1 && fIntegral[ibin+1] == prob[i]) {
4486  if (fIntegral[ibin+2] == prob[i]) ibin++;
4487  else break;
4488  }
4489  q[i] = GetBinLowEdge(ibin+1);
4490  const Double_t dint = fIntegral[ibin+1]-fIntegral[ibin];
4491  if (dint > 0) q[i] += GetBinWidth(ibin+1)*(prob[i]-fIntegral[ibin])/dint;
4492  }
4493 
4494  if (!probSum) delete [] prob;
4495  return nq;
4496 }
4497 
4498 ////////////////////////////////////////////////////////////////////////////////
4499 /// Decode string choptin and fill fitOption structure.
4500 
4501 Int_t TH1::FitOptionsMake(Option_t *choptin, Foption_t &fitOption)
4502 {
4504  return 1;
4505 }
4506 
4507 ////////////////////////////////////////////////////////////////////////////////
4508 /// Compute Initial values of parameters for a gaussian.
4509 
4510 void H1InitGaus()
4511 {
4512  Double_t allcha, sumx, sumx2, x, val, stddev, mean;
4513  Int_t bin;
4514  const Double_t sqrtpi = 2.506628;
4515 
4516  // - Compute mean value and StdDev of the histogram in the given range
4518  TH1 *curHist = (TH1*)hFitter->GetObjectFit();
4519  Int_t hxfirst = hFitter->GetXfirst();
4520  Int_t hxlast = hFitter->GetXlast();
4521  Double_t valmax = curHist->GetBinContent(hxfirst);
4522  Double_t binwidx = curHist->GetBinWidth(hxfirst);
4523  allcha = sumx = sumx2 = 0;
4524  for (bin=hxfirst;bin<=hxlast;bin++) {
4525  x = curHist->GetBinCenter(bin);
4526  val = TMath::Abs(curHist->GetBinContent(bin));
4527  if (val > valmax) valmax = val;
4528  sumx += val*x;
4529  sumx2 += val*x*x;
4530  allcha += val;
4531  }
4532  if (allcha == 0) return;
4533  mean = sumx/allcha;
4534  stddev = sumx2/allcha - mean*mean;
4535  if (stddev > 0) stddev = TMath::Sqrt(stddev);
4536  else stddev = 0;
4537  if (stddev == 0) stddev = binwidx*(hxlast-hxfirst+1)/4;
4538  //if the distribution is really gaussian, the best approximation
4539  //is binwidx*allcha/(sqrtpi*stddev)
4540  //However, in case of non-gaussian tails, this underestimates
4541  //the normalisation constant. In this case the maximum value
4542  //is a better approximation.
4543  //We take the average of both quantities
4544  Double_t constant = 0.5*(valmax+binwidx*allcha/(sqrtpi*stddev));
4545 
4546  //In case the mean value is outside the histo limits and
4547  //the StdDev is bigger than the range, we take
4548  // mean = center of bins
4549  // stddev = half range
4550  Double_t xmin = curHist->GetXaxis()->GetXmin();
4551  Double_t xmax = curHist->GetXaxis()->GetXmax();
4552  if ((mean < xmin || mean > xmax) && stddev > (xmax-xmin)) {
4553  mean = 0.5*(xmax+xmin);
4554  stddev = 0.5*(xmax-xmin);
4555  }
4556  TF1 *f1 = (TF1*)hFitter->GetUserFunc();
4557  f1->SetParameter(0,constant);
4558  f1->SetParameter(1,mean);
4559  f1->SetParameter(2,stddev);
4560  f1->SetParLimits(2,0,10*stddev);
4561 }
4562 
4563 ////////////////////////////////////////////////////////////////////////////////
4564 /// Compute Initial values of parameters for an exponential.
4565 
4566 void H1InitExpo()
4567 {
4568  Double_t constant, slope;
4569  Int_t ifail;
4571  Int_t hxfirst = hFitter->GetXfirst();
4572  Int_t hxlast = hFitter->GetXlast();
4573  Int_t nchanx = hxlast - hxfirst + 1;
4574 
4575  H1LeastSquareLinearFit(-nchanx, constant, slope, ifail);
4576 
4577  TF1 *f1 = (TF1*)hFitter->GetUserFunc();
4578  f1->SetParameter(0,constant);
4579  f1->SetParameter(1,slope);
4580 
4581 }
4582 
4583 ////////////////////////////////////////////////////////////////////////////////
4584 /// Compute Initial values of parameters for a polynom.
4585 
4586 void H1InitPolynom()
4587 {
4588  Double_t fitpar[25];
4589 
4591  TF1 *f1 = (TF1*)hFitter->GetUserFunc();
4592  Int_t hxfirst = hFitter->GetXfirst();
4593  Int_t hxlast = hFitter->GetXlast();
4594  Int_t nchanx = hxlast - hxfirst + 1;
4595  Int_t npar = f1->GetNpar();
4596 
4597  if (nchanx <=1 || npar == 1) {
4598  TH1 *curHist = (TH1*)hFitter->GetObjectFit();
4599  fitpar[0] = curHist->GetSumOfWeights()/Double_t(nchanx);
4600  } else {
4601  H1LeastSquareFit( nchanx, npar, fitpar);
4602  }
4603  for (Int_t i=0;i<npar;i++) f1->SetParameter(i, fitpar[i]);
4604 }
4605 
4606 ////////////////////////////////////////////////////////////////////////////////
4607 /// Least squares lpolynomial fitting without weights.
4608 ///
4609 /// \param[in] n number of points to fit
4610 /// \param[in] m number of parameters
4611 /// \param[in] a array of parameters
4612 ///
4613 /// based on CERNLIB routine LSQ: Translated to C++ by Rene Brun
4614 /// (E.Keil. revised by B.Schorr, 23.10.1981.)
4615 
4617 {
4618  const Double_t zero = 0.;
4619  const Double_t one = 1.;
4620  const Int_t idim = 20;
4621 
4622  Double_t b[400] /* was [20][20] */;
4623  Int_t i, k, l, ifail;
4624  Double_t power;
4625  Double_t da[20], xk, yk;
4626 
4627  if (m <= 2) {
4628  H1LeastSquareLinearFit(n, a[0], a[1], ifail);
4629  return;
4630  }
4631  if (m > idim || m > n) return;
4632  b[0] = Double_t(n);
4633  da[0] = zero;
4634  for (l = 2; l <= m; ++l) {
4635  b[l-1] = zero;
4636  b[m + l*20 - 21] = zero;
4637  da[l-1] = zero;
4638  }
4640  TH1 *curHist = (TH1*)hFitter->GetObjectFit();
4641  Int_t hxfirst = hFitter->GetXfirst();
4642  Int_t hxlast = hFitter->GetXlast();
4643  for (k = hxfirst; k <= hxlast; ++k) {
4644  xk = curHist->GetBinCenter(k);
4645  yk = curHist->GetBinContent(k);
4646  power = one;
4647  da[0] += yk;
4648  for (l = 2; l <= m; ++l) {
4649  power *= xk;
4650  b[l-1] += power;
4651  da[l-1] += power*yk;
4652  }
4653  for (l = 2; l <= m; ++l) {
4654  power *= xk;
4655  b[m + l*20 - 21] += power;
4656  }
4657  }
4658  for (i = 3; i <= m; ++i) {
4659  for (k = i; k <= m; ++k) {
4660  b[k - 1 + (i-1)*20 - 21] = b[k + (i-2)*20 - 21];
4661  }
4662  }
4663  H1LeastSquareSeqnd(m, b, idim, ifail, 1, da);
4664 
4665  for (i=0; i<m; ++i) a[i] = da[i];
4666 
4667 }
4668 
4669 ////////////////////////////////////////////////////////////////////////////////
4670 /// Least square linear fit without weights.
4671 ///
4672 /// extracted from CERNLIB LLSQ: Translated to C++ by Rene Brun
4673 /// (added to LSQ by B. Schorr, 15.02.1982.)
4674 
4675 void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail)
4676 {
4677  Double_t xbar, ybar, x2bar;
4678  Int_t i, n;
4679  Double_t xybar;
4680  Double_t fn, xk, yk;
4681  Double_t det;
4682 
4683  n = TMath::Abs(ndata);
4684  ifail = -2;
4685  xbar = ybar = x2bar = xybar = 0;
4687  TH1 *curHist = (TH1*)hFitter->GetObjectFit();
4688  Int_t hxfirst = hFitter->GetXfirst();
4689  Int_t hxlast = hFitter->GetXlast();
4690  for (i = hxfirst; i <= hxlast; ++i) {
4691  xk = curHist->GetBinCenter(i);
4692  yk = curHist->GetBinContent(i);
4693  if (ndata < 0) {
4694  if (yk <= 0) yk = 1e-9;
4695  yk = TMath::Log(yk);
4696  }
4697  xbar += xk;
4698  ybar += yk;
4699  x2bar += xk*xk;
4700  xybar += xk*yk;
4701  }
4702  fn = Double_t(n);
4703  det = fn*x2bar - xbar*xbar;
4704  ifail = -1;
4705  if (det <= 0) {
4706  a0 = ybar/fn;
4707  a1 = 0;
4708  return;
4709  }
4710  ifail = 0;
4711  a0 = (x2bar*ybar - xbar*xybar) / det;
4712  a1 = (fn*xybar - xbar*ybar) / det;
4713 
4714 }
4715 
4716 ////////////////////////////////////////////////////////////////////////////////
4717 /// Extracted from CERN Program library routine DSEQN.
4718 ///
4719 /// Translated to C++ by Rene Brun
4720 
4721 void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b)
4722 {
4723  Int_t a_dim1, a_offset, b_dim1, b_offset;
4724  Int_t nmjp1, i, j, l;
4725  Int_t im1, jp1, nm1, nmi;
4726  Double_t s1, s21, s22;
4727  const Double_t one = 1.;
4728 
4729  /* Parameter adjustments */
4730  b_dim1 = idim;
4731  b_offset = b_dim1 + 1;
4732  b -= b_offset;
4733  a_dim1 = idim;
4734  a_offset = a_dim1 + 1;
4735  a -= a_offset;
4736 
4737  if (idim < n) return;
4738 
4739  ifail = 0;
4740  for (j = 1; j <= n; ++j) {
4741  if (a[j + j*a_dim1] <= 0) { ifail = -1; return; }
4742  a[j + j*a_dim1] = one / a[j + j*a_dim1];
4743  if (j == n) continue;
4744  jp1 = j + 1;
4745  for (l = jp1; l <= n; ++l) {
4746  a[j + l*a_dim1] = a[j + j*a_dim1] * a[l + j*a_dim1];
4747  s1 = -a[l + (j+1)*a_dim1];
4748  for (i = 1; i <= j; ++i) { s1 = a[l + i*a_dim1] * a[i + (j+1)*a_dim1] + s1; }
4749  a[l + (j+1)*a_dim1] = -s1;
4750  }
4751  }
4752  if (k <= 0) return;
4753 
4754  for (l = 1; l <= k; ++l) {
4755  b[l*b_dim1 + 1] = a[a_dim1 + 1]*b[l*b_dim1 + 1];
4756  }
4757  if (n == 1) return;
4758  for (l = 1; l <= k; ++l) {
4759  for (i = 2; i <= n; ++i) {
4760  im1 = i - 1;
4761  s21 = -b[i + l*b_dim1];
4762  for (j = 1; j <= im1; ++j) {
4763  s21 = a[i + j*a_dim1]*b[j + l*b_dim1] + s21;
4764  }
4765  b[i + l*b_dim1] = -a[i + i*a_dim1]*s21;
4766  }
4767  nm1 = n - 1;
4768  for (i = 1; i <= nm1; ++i) {
4769  nmi = n - i;
4770  s22 = -b[nmi + l*b_dim1];
4771  for (j = 1; j <= i; ++j) {
4772  nmjp1 = n - j + 1;
4773  s22 = a[nmi + nmjp1*a_dim1]*b[nmjp1 + l*b_dim1] + s22;
4774  }
4775  b[nmi + l*b_dim1] = -s22;
4776  }
4777  }
4778 }
4779 
4780 ////////////////////////////////////////////////////////////////////////////////
4781 /// Return Global bin number corresponding to binx,y,z.
4782 ///
4783 /// 2-D and 3-D histograms are represented with a one dimensional
4784 /// structure.
4785 /// This has the advantage that all existing functions, such as
4786 /// GetBinContent, GetBinError, GetBinFunction work for all dimensions.
4787 ///
4788 /// In case of a TH1x, returns binx directly.
4789 /// see TH1::GetBinXYZ for the inverse transformation.
4790 ///
4791 /// Convention for numbering bins
4792 ///
4793 /// For all histogram types: nbins, xlow, xup
4794 ///
4795 /// - bin = 0; underflow bin
4796 /// - bin = 1; first bin with low-edge xlow INCLUDED
4797 /// - bin = nbins; last bin with upper-edge xup EXCLUDED
4798 /// - bin = nbins+1; overflow bin
4799 ///
4800 /// In case of 2-D or 3-D histograms, a "global bin" number is defined.
4801 /// For example, assuming a 3-D histogram with binx,biny,binz, the function
4802 ///
4803 /// ~~~ {.cpp}
4804 /// Int_t bin = h->GetBin(binx,biny,binz);
4805 /// ~~~
4806 ///
4807 /// returns a global/linearized bin number. This global bin is useful
4808 /// to access the bin information independently of the dimension.
4809 
4810 Int_t TH1::GetBin(Int_t binx, Int_t, Int_t) const
4811 {
4812  Int_t ofx = fXaxis.GetNbins() + 1; // overflow bin
4813  if (binx < 0) binx = 0;
4814  if (binx > ofx) binx = ofx;
4815 
4816  return binx;
4817 }
4818 
4819 ////////////////////////////////////////////////////////////////////////////////
4820 /// Return binx, biny, binz corresponding to the global bin number globalbin
4821 /// see TH1::GetBin function above
4822 
4823 void TH1::GetBinXYZ(Int_t binglobal, Int_t &binx, Int_t &biny, Int_t &binz) const
4824 {
4825  Int_t nx = fXaxis.GetNbins()+2;
4826  Int_t ny = fYaxis.GetNbins()+2;
4827 
4828  if (GetDimension() == 1) {
4829  binx = binglobal%nx;
4830  biny = 0;
4831  binz = 0;
4832  return;
4833  }
4834  if (GetDimension() == 2) {
4835  binx = binglobal%nx;
4836  biny = ((binglobal-binx)/nx)%ny;
4837  binz = 0;
4838  return;
4839  }
4840  if (GetDimension() == 3) {
4841  binx = binglobal%nx;
4842  biny = ((binglobal-binx)/nx)%ny;
4843  binz = ((binglobal-binx)/nx -biny)/ny;
4844  }
4845 }
4846 
4847 ////////////////////////////////////////////////////////////////////////////////
4848 /// Return a random number distributed according the histogram bin contents.
4849 /// This function checks if the bins integral exists. If not, the integral
4850 /// is evaluated, normalized to one.
4851 ///
4852 /// The integral is automatically recomputed if the number of entries
4853 /// is not the same then when the integral was computed.
4854 /// NB Only valid for 1-d histograms. Use GetRandom2 or 3 otherwise.
4855 /// If the histogram has a bin with negative content a NaN is returned
4856 
4857 Double_t TH1::GetRandom() const
4858 {
4859  if (fDimension > 1) {
4860  Error("GetRandom","Function only valid for 1-d histograms");
4861  return 0;
4862  }
4863  Int_t nbinsx = GetNbinsX();
4864  Double_t integral = 0;
4865  // compute integral checking that all bins have positive content (see ROOT-5894)
4866  if (fIntegral) {
4867  if (fIntegral[nbinsx+1] != fEntries) integral = ((TH1*)this)->ComputeIntegral(true);
4868  else integral = fIntegral[nbinsx];
4869  } else {
4870  integral = ((TH1*)this)->ComputeIntegral(true);
4871  }
4872  if (integral == 0) return 0;
4873  // return a NaN in case some bins have negative content
4874  if (integral == TMath::QuietNaN() ) return TMath::QuietNaN();
4875 
4876  Double_t r1 = gRandom->Rndm();
4877  Int_t ibin = TMath::BinarySearch(nbinsx,fIntegral,r1);
4878  Double_t x = GetBinLowEdge(ibin+1);
4879  if (r1 > fIntegral[ibin]) x +=
4880  GetBinWidth(ibin+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1] - fIntegral[ibin]);
4881  return x;
4882 }
4883 
4884 ////////////////////////////////////////////////////////////////////////////////
4885 /// Return content of bin number bin.
4886 ///
4887 /// Implemented in TH1C,S,F,D
4888 ///
4889 /// Convention for numbering bins
4890 ///
4891 /// For all histogram types: nbins, xlow, xup
4892 ///
4893 /// - bin = 0; underflow bin
4894 /// - bin = 1; first bin with low-edge xlow INCLUDED
4895 /// - bin = nbins; last bin with upper-edge xup EXCLUDED
4896 /// - bin = nbins+1; overflow bin
4897 ///
4898 /// In case of 2-D or 3-D histograms, a "global bin" number is defined.
4899 /// For example, assuming a 3-D histogram with binx,biny,binz, the function
4900 ///
4901 /// ~~~ {.cpp}
4902 /// Int_t bin = h->GetBin(binx,biny,binz);
4903 /// ~~~
4904 ///
4905 /// returns a global/linearized bin number. This global bin is useful
4906 /// to access the bin information independently of the dimension.
4907 
4909 {
4910  if (fBuffer) const_cast<TH1*>(this)->BufferEmpty();
4911  if (bin < 0) bin = 0;
4912  if (bin >= fNcells) bin = fNcells-1;
4913 
4914  return RetrieveBinContent(bin);
4915 }
4916 
4917 ////////////////////////////////////////////////////////////////////////////////
4918 /// Compute first binx in the range [firstx,lastx] for which
4919 /// diff = abs(bin_content-c) <= maxdiff
4920 ///
4921 /// In case several bins in the specified range with diff=0 are found
4922 /// the first bin found is returned in binx.
4923 /// In case several bins in the specified range satisfy diff <=maxdiff
4924 /// the bin with the smallest difference is returned in binx.
4925 /// In all cases the function returns the smallest difference.
4926 ///
4927 /// NOTE1: if firstx <= 0, firstx is set to bin 1
4928 /// if (lastx < firstx then firstx is set to the number of bins
4929 /// ie if firstx=0 and lastx=0 (default) the search is on all bins.
4930 ///
4931 /// NOTE2: if maxdiff=0 (default), the first bin with content=c is returned.
4932 
4933 Double_t TH1::GetBinWithContent(Double_t c, Int_t &binx, Int_t firstx, Int_t lastx,Double_t maxdiff) const
4934 {
4935  if (fDimension > 1) {
4936  binx = 0;
4937  Error("GetBinWithContent","function is only valid for 1-D histograms");
4938  return 0;
4939  }
4940 
4941  if (fBuffer) ((TH1*)this)->BufferEmpty();
4942 
4943  if (firstx <= 0) firstx = 1;
4944  if (lastx < firstx) lastx = fXaxis.GetNbins();
4945  Int_t binminx = 0;
4946  Double_t diff, curmax = 1.e240;
4947  for (Int_t i=firstx;i<=lastx;i++) {
4948  diff = TMath::Abs(RetrieveBinContent(i)-c);
4949  if (diff <= 0) {binx = i; return diff;}
4950  if (diff < curmax && diff <= maxdiff) {curmax = diff, binminx=i;}
4951  }
4952  binx = binminx;
4953  return curmax;
4954 }
4955 
4956 ////////////////////////////////////////////////////////////////////////////////
4957 /// Given a point x, approximates the value via linear interpolation
4958 /// based on the two nearest bin centers
4959 ///
4960 /// Andy Mastbaum 10/21/08
4961 
4963 {
4964  if (fBuffer) ((TH1*)this)->BufferEmpty();
4965 
4966  Int_t xbin = fXaxis.FindFixBin(x);
4967  Double_t x0,x1,y0,y1;
4968 
4969  if(x<=GetBinCenter(1)) {
4970  return RetrieveBinContent(1);
4971  } else if(x>=GetBinCenter(GetNbinsX())) {
4972  return RetrieveBinContent(GetNbinsX());
4973  } else {
4974  if(x<=GetBinCenter(xbin)) {
4975  y0 = RetrieveBinContent(xbin-1);
4976  x0 = GetBinCenter(xbin-1);
4977  y1 = RetrieveBinContent(xbin);
4978  x1 = GetBinCenter(xbin);
4979  } else {
4980  y0 = RetrieveBinContent(xbin);
4981  x0 = GetBinCenter(xbin);
4982  y1 = RetrieveBinContent(xbin+1);
4983  x1 = GetBinCenter(xbin+1);
4984  }
4985  return y0 + (x-x0)*((y1-y0)/(x1-x0));
4986  }
4987 }
4988 
4989 ////////////////////////////////////////////////////////////////////////////////
4990 /// 2d Interpolation. Not yet implemented.
4991 
4993 {
4994  Error("Interpolate","This function must be called with 1 argument for a TH1");
4995  return 0;
4996 }
4997 
4998 ////////////////////////////////////////////////////////////////////////////////
4999 /// 3d Interpolation. Not yet implemented.
5000 
5002 {
5003  Error("Interpolate","This function must be called with 1 argument for a TH1");
5004  return 0;
5005 }
5006 
5007 ///////////////////////////////////////////////////////////////////////////////
5008 /// Check if an histogram is empty
5009 /// (this a protected method used mainly by TH1Merger )
5010 
5011 Bool_t TH1::IsEmpty() const
5012 {
5013  // if fTsumw or fentries are not zero histogram is not empty
5014  // need to use GetEntries() instead of fEntries in case of bugger histograms
5015  // so we will flash the buffer
5016  if (fTsumw != 0) return kFALSE;
5017  if (GetEntries() != 0) return kFALSE;
5018  // case fTSumw == 0 amd entries are also zero
5019  // this should not really happening, but if one sets content by hand
5020  // it can happen. a call to ResetStats() should be done in such cases
5021  double sumw = 0;
5022  for (int i = 0; i< GetNcells(); ++i) sumw += RetrieveBinContent(i);
5023  return (sumw != 0) ? kFALSE : kTRUE;
5024 }
5025 
5026 ////////////////////////////////////////////////////////////////////////////////
5027 /// Return true if the bin is overflow.
5028 
5029 Bool_t TH1::IsBinOverflow(Int_t bin, Int_t iaxis) const
5030 {
5031  Int_t binx, biny, binz;
5032  GetBinXYZ(bin, binx, biny, binz);
5033 
5034  if (iaxis == 0) {
5035  if ( fDimension == 1 )
5036  return binx >= GetNbinsX() + 1;
5037  if ( fDimension == 2 )
5038  return (binx >= GetNbinsX() + 1) ||
5039  (biny >= GetNbinsY() + 1);
5040  if ( fDimension == 3 )
5041  return (binx >= GetNbinsX() + 1) ||
5042  (biny >= GetNbinsY() + 1) ||
5043  (binz >= GetNbinsZ() + 1);
5044  return kFALSE;
5045  }
5046  if (iaxis == 1)
5047  return binx >= GetNbinsX() + 1;
5048  if (iaxis == 2)
5049  return biny >= GetNbinsY() + 1;
5050  if (iaxis == 3)
5051  return binz >= GetNbinsZ() + 1;
5052 
5053  Error("IsBinOverflow","Invalid axis value");
5054  return kFALSE;
5055 }
5056 
5057 ////////////////////////////////////////////////////////////////////////////////
5058 /// Return true if the bin is underflow.
5059 /// If iaxis = 0 make OR with all axes otherwise check only for the given axis
5060 
5061 Bool_t TH1::IsBinUnderflow(Int_t bin, Int_t iaxis) const
5062 {
5063  Int_t binx, biny, binz;
5064  GetBinXYZ(bin, binx, biny, binz);
5065 
5066  if (iaxis == 0) {
5067  if ( fDimension == 1 )
5068  return (binx <= 0);
5069  else if ( fDimension == 2 )
5070  return (binx <= 0 || biny <= 0);
5071  else if ( fDimension == 3 )
5072  return (binx <= 0 || biny <= 0 || binz <= 0);
5073  else
5074  return kFALSE;
5075  }
5076  if (iaxis == 1)
5077  return (binx <= 0);
5078  if (iaxis == 2)
5079  return (biny <= 0);
5080  if (iaxis == 3)
5081  return (binz <= 0);
5082 
5083  Error("IsBinUnderflow","Invalid axis value");
5084  return kFALSE;
5085 }
5086 
5087 ////////////////////////////////////////////////////////////////////////////////
5088 /// Reduce the number of bins for the axis passed in the option to the number of bins having a label.
5089 /// The method will remove only the extra bins existing after the last "labeled" bin.
5090 /// Note that if there are "un-labeled" bins present between "labeled" bins they will not be removed
5091 
5092 void TH1::LabelsDeflate(Option_t *ax)
5093 {
5094  Int_t iaxis = AxisChoice(ax);
5095  TAxis *axis = 0;
5096  if (iaxis == 1) axis = GetXaxis();
5097  if (iaxis == 2) axis = GetYaxis();
5098  if (iaxis == 3) axis = GetZaxis();
5099  if (!axis) {
5100  Error("LabelsDeflate","Invalid axis option %s",ax);
5101  return;
5102  }
5103  if (!axis->GetLabels()) return;
5104 
5105  // find bin with last labels
5106  // bin number is object ID in list of labels
5107  // therefore max bin number is number of bins of the deflated histograms
5108  TIter next(axis->GetLabels());
5109  TObject *obj;
5110  Int_t nbins = 0;
5111  while ((obj = next())) {
5112  Int_t ibin = obj->GetUniqueID();
5113  if (ibin > nbins) nbins = ibin;
5114  }
5115  if (nbins < 1) nbins = 1;
5116 
5117  // Do nothing in case it was the last bin
5118  if (nbins==axis->GetNbins()) return;
5119 
5120  TH1 *hold = (TH1*)IsA()->New();
5121  R__ASSERT(hold);
5122  hold->SetDirectory(0);
5123  Copy(*hold);
5124 
5125  Bool_t timedisp = axis->GetTimeDisplay();
5126  Double_t xmin = axis->GetXmin();
5127  Double_t xmax = axis->GetBinUpEdge(nbins);
5128  if (xmax <= xmin) xmax = xmin +nbins;
5129  axis->SetRange(0,0);
5130  axis->Set(nbins,xmin,xmax);
5131  SetBinsLength(-1); // reset the number of cells
5132  Int_t errors = fSumw2.fN;
5133  if (errors) fSumw2.Set(fNcells);
5134  axis->SetTimeDisplay(timedisp);
5135  // reset histogram content
5136  Reset("ICE");
5137 
5138  //now loop on all bins and refill
5139  // NOTE that if the bins without labels have content
5140  // it will be put in the underflow/overflow.
5141  // For this reason we use AddBinContent method
5142  Double_t oldEntries = fEntries;
5143  Int_t bin,binx,biny,binz;
5144  for (bin=0; bin < hold->fNcells; ++bin) {
5145  hold->GetBinXYZ(bin,binx,biny,binz);
5146  Int_t ibin = GetBin(binx,biny,binz);
5147  Double_t cu = hold->RetrieveBinContent(bin);
5148  AddBinContent(ibin,cu);
5149  if (errors) {
5150  fSumw2.fArray[ibin] += hold->fSumw2.fArray[bin];
5151  }
5152  }
5153  fEntries = oldEntries;
5154  delete hold;
5155 }
5156 
5157 ////////////////////////////////////////////////////////////////////////////////
5158 /// Double the number of bins for axis.
5159 /// Refill histogram
5160 /// This function is called by TAxis::FindBin(const char *label)
5161 
5162 void TH1::LabelsInflate(Option_t *ax)
5163 {
5164  Int_t iaxis = AxisChoice(ax);
5165  TAxis *axis = 0;
5166  if (iaxis == 1) axis = GetXaxis();
5167  if (iaxis == 2) axis = GetYaxis();
5168  if (iaxis == 3) axis = GetZaxis();
5169  if (!axis) return;
5170 
5171  TH1 *hold = (TH1*)IsA()->New();;
5172  hold->SetDirectory(0);
5173  Copy(*hold);
5174 
5175  Bool_t timedisp = axis->GetTimeDisplay();
5176  Int_t nbins = axis->GetNbins();
5177  Double_t xmin = axis->GetXmin();
5178  Double_t xmax = axis->GetXmax();
5179  xmax = xmin + 2*(xmax-xmin);
5180  axis->SetRange(0,0);
5181  // double the bins and recompute ncells
5182  axis->Set(2*nbins,xmin,xmax);
5183  SetBinsLength(-1);
5184  Int_t errors = fSumw2.fN;
5185  if (errors) fSumw2.Set(fNcells);
5186  axis->SetTimeDisplay(timedisp);
5187 
5188  Reset("ICE"); // reset content and error
5189 
5190  //now loop on all bins and refill
5191  Double_t oldEntries = fEntries;
5192  Int_t bin,ibin,binx,biny,binz;
5193  for (ibin =0; ibin < hold->fNcells; ibin++) {
5194  // get the binx,y,z values . The x-y-z (axis) bin values will stay the same between new-old after the expanding
5195  hold->GetBinXYZ(ibin,binx,biny,binz);
5196  bin = GetBin(binx,biny,binz);
5197 
5198  // underflow and overflow will be cleaned up because their meaning has been altered
5199  if (hold->IsBinUnderflow(ibin,iaxis) || hold->IsBinOverflow(ibin,iaxis)) {
5200  continue;
5201  }
5202  else {
5203  AddBinContent(bin, hold->RetrieveBinContent(ibin));
5204  if (errors) fSumw2.fArray[bin] += hold->fSumw2.fArray[ibin];
5205  }
5206  }
5207  fEntries = oldEntries;
5208  delete hold;
5209 }
5210 
5211 ////////////////////////////////////////////////////////////////////////////////
5212 /// Set option(s) to draw axis with labels
5213 /// \param[in] option
5214 /// - "a" sort by alphabetic order
5215 /// - ">" sort by decreasing values
5216 /// - "<" sort by increasing values
5217 /// - "h" draw labels horizontal
5218 /// - "v" draw labels vertical
5219 /// - "u" draw labels up (end of label right adjusted)
5220 /// - "d" draw labels down (start of label left adjusted)
5221 /// \param[in] ax axis
5222 
5223 void TH1::LabelsOption(Option_t *option, Option_t *ax)
5224 {
5225  Int_t iaxis = AxisChoice(ax);
5226  TAxis *axis = 0;
5227  if (iaxis == 1) axis = GetXaxis();
5228  if (iaxis == 2) axis = GetYaxis();
5229  if (iaxis == 3) axis = GetZaxis();
5230  if (!axis) return;
5231  THashList *labels = axis->GetLabels();
5232  if (!labels) {
5233  Warning("LabelsOption","Cannot sort. No labels");
5234  return;
5235  }
5236  TString opt = option;
5237  opt.ToLower();
5238  if (opt.Contains("h")) {
5239  axis->SetBit(TAxis::kLabelsHori);
5242  axis->ResetBit(TAxis::kLabelsUp);
5243  }
5244  if (opt.Contains("v")) {
5245  axis->SetBit(TAxis::kLabelsVert);
5248  axis->ResetBit(TAxis::kLabelsUp);
5249  }
5250  if (opt.Contains("u")) {
5251  axis->SetBit(TAxis::kLabelsUp);
5255  }
5256  if (opt.Contains("d")) {
5257  axis->SetBit(TAxis::kLabelsDown);
5260  axis->ResetBit(TAxis::kLabelsUp);
5261  }
5262  Int_t sort = -1;
5263  if (opt.Contains("a")) sort = 0;
5264  if (opt.Contains(">")) sort = 1;
5265  if (opt.Contains("<")) sort = 2;
5266  if (sort < 0) return;
5267  if (sort > 0 && GetDimension() > 2) {
5268  Error("LabelsOption","Sorting by value not implemented for 3-D histograms");
5269  return;
5270  }
5271 
5272  Double_t entries = fEntries;
5273  Int_t n = TMath::Min(axis->GetNbins(), labels->GetSize());
5274  std::vector<Int_t> a(n+2);
5275 
5276  Int_t i,j,k;
5277  std::vector<Double_t> cont;
5278  std::vector<Double_t> errors;
5279  THashList *labold = new THashList(labels->GetSize(),1);
5280  TIter nextold(labels);
5281  TObject *obj;
5282  while ((obj=nextold())) {
5283  labold->Add(obj);
5284  }
5285  labels->Clear();
5286 
5287  // delete buffer if it is there since bins will be reordered.
5288  if (fBuffer) BufferEmpty(1);
5289 
5290  if (sort > 0) {
5291  //---sort by values of bins
5292  if (GetDimension() == 1) {
5293  cont.resize(n);
5294  if (fSumw2.fN) errors.resize(n);
5295  for (i=1;i<=n;i++) {
5296  cont[i-1] = GetBinContent(i);
5297  if (!errors.empty()) errors[i-1] = GetBinError(i);
5298  }
5299  if (sort ==1) TMath::Sort(n,cont.data(),a.data(),kTRUE); //sort by decreasing values
5300  else TMath::Sort(n,cont.data(),a.data(),kFALSE); //sort by increasing values
5301  for (i=1;i<=n;i++) {
5302  SetBinContent(i,cont[a[i-1]]);
5303  if (!errors.empty()) SetBinError(i,errors[a[i-1]]);
5304  }
5305  for (i=1;i<=n;i++) {
5306  obj = labold->At(a[i-1]);
5307  labels->Add(obj);
5308  obj->SetUniqueID(i);
5309  }
5310  } else if (GetDimension()== 2) {
5311  std::vector<Double_t> pcont(n+2);
5312  Int_t nx = fXaxis.GetNbins();
5313  Int_t ny = fYaxis.GetNbins();
5314  cont.resize( (nx+2)*(ny+2));
5315  if (fSumw2.fN) errors.resize( (nx+2)*(ny+2));
5316  for (i=1;i<=nx;i++) {
5317  for (j=1;j<=ny;j++) {
5318  cont[i+nx*j] = GetBinContent(i,j);
5319  if (!errors.empty()) errors[i+nx*j] = GetBinError(i,j);
5320  if (axis == GetXaxis()) k = i;
5321  else k = j;
5322  pcont[k-1] += cont[i+nx*j];
5323  }
5324  }
5325  if (sort ==1) TMath::Sort(n,pcont.data(),a.data(),kTRUE); //sort by decreasing values
5326  else TMath::Sort(n,pcont.data(),a.data(),kFALSE); //sort by increasing values
5327  for (i=0;i<n;i++) {
5328  obj = labold->At(a[i]);
5329  labels->Add(obj);
5330  obj->SetUniqueID(i+1);
5331  }
5332  if (axis == GetXaxis()) {
5333  for (i=1;i<=n;i++) {
5334  for (j=1;j<=ny;j++) {
5335  SetBinContent(i,j,cont[a[i-1]+1+nx*j]);
5336  if (!errors.empty()) SetBinError(i,j,errors[a[i-1]+1+nx*j]);
5337  }
5338  }
5339  }
5340  else {
5341  // using y axis
5342  for (i=1;i<=nx;i++) {
5343  for (j=1;j<=n;j++) {
5344  SetBinContent(i,j,cont[i+nx*(a[j-1]+1)]);
5345  if (!errors.empty()) SetBinError(i,j,errors[i+nx*(a[j-1]+1)]);
5346  }
5347  }
5348  }
5349  } else {
5350  //to be implemented for 3d
5351  }
5352  } else {
5353  //---alphabetic sort
5354  const UInt_t kUsed = 1<<18;
5355  TObject *objk=0;
5356  a[0] = 0;
5357  a[n+1] = n+1;
5358  for (i=1;i<=n;i++) {
5359  const char *label = "zzzzzzzzzzzz";
5360  for (j=1;j<=n;j++) {
5361  obj = labold->At(j-1);
5362  if (!obj) continue;
5363  if (obj->TestBit(kUsed)) continue;
5364  //use strcasecmp for case non-sensitive sort (may be an option)
5365  if (strcmp(label,obj->GetName()) < 0) continue;
5366  objk = obj;
5367  a[i] = j;
5368  label = obj->GetName();
5369  }
5370  if (objk) {
5371  objk->SetUniqueID(i);
5372  labels->Add(objk);
5373  objk->SetBit(kUsed);
5374  }
5375  }
5376  for (i=1;i<=n;i++) {
5377  obj = labels->At(i-1);
5378  if (!obj) continue;
5379  obj->ResetBit(kUsed);
5380  }
5381 
5382  if (GetDimension() == 1) {
5383  cont.resize(n+2);
5384  if (fSumw2.fN) errors.resize(n+2);
5385  for (i=1;i<=n;i++) {
5386  cont[i] = GetBinContent(a[i]);
5387  if (!errors.empty()) errors[i] = GetBinError(a[i]);
5388  }
5389  for (i=1;i<=n;i++) {
5390  SetBinContent(i,cont[i]);
5391  if (!errors.empty()) SetBinError(i,errors[i]);
5392  }
5393  } else if (GetDimension()== 2) {
5394  Int_t nx = fXaxis.GetNbins()+2;
5395  Int_t ny = fYaxis.GetNbins()+2;
5396  cont.resize(nx*ny);
5397  if (fSumw2.fN) errors.resize(nx*ny);
5398  for (i=0;i<nx;i++) {
5399  for (j=0;j<ny;j++) {
5400  cont[i+nx*j] = GetBinContent(i,j);
5401  if (!errors.empty()) errors[i+nx*j] = GetBinError(i,j);
5402  }
5403  }
5404  if (axis == GetXaxis()) {
5405  for (i=1;i<=n;i++) {
5406  for (j=0;j<ny;j++) {
5407  SetBinContent(i,j,cont[a[i]+nx*j]);
5408  if (!errors.empty()) SetBinError(i,j,errors[a[i]+nx*j]);
5409  }
5410  }
5411  } else {
5412  for (i=0;i<nx;i++) {
5413  for (j=1;j<=n;j++) {
5414  SetBinContent(i,j,cont[i+nx*a[j]]);
5415  if (!errors.empty()) SetBinError(i,j,errors[i+nx*a[j]]);
5416  }
5417  }
5418  }
5419  } else {
5420  Int_t nx = fXaxis.GetNbins()+2;
5421  Int_t ny = fYaxis.GetNbins()+2;
5422  Int_t nz = fZaxis.GetNbins()+2;
5423  cont.resize(nx*ny*nz);
5424  if (fSumw2.fN) errors.resize(nx*ny*nz);
5425  for (i=0;i<nx;i++) {
5426  for (j=0;j<ny;j++) {
5427  for (k=0;k<nz;k++) {
5428  cont[i+nx*(j+ny*k)] = GetBinContent(i,j,k);
5429  if (!errors.empty()) errors[i+nx*(j+ny*k)] = GetBinError(i,j,k);
5430  }
5431  }
5432  }
5433  if (axis == GetXaxis()) {
5434  // labels on x axis
5435  for (i=1;i<=n;i++) {
5436  for (j=0;j<ny;j++) {
5437  for (k=0;k<nz;k++) {
5438  SetBinContent(i,j,k,cont[a[i]+nx*(j+ny*k)]);
5439  if (!errors.empty()) SetBinError(i,j,k,errors[a[i]+nx*(j+ny*k)]);
5440  }
5441  }
5442  }
5443  }
5444  else if (axis == GetYaxis()) {
5445  // labels on y axis
5446  for (i=0;i<nx;i++) {
5447  for (j=1;j<=n;j++) {
5448  for (k=0;k<nz;k++) {
5449  SetBinContent(i,j,k,cont[i+nx*(a[j]+ny*k)]);
5450  if (!errors.empty()) SetBinError(i,j,k,errors[i+nx*(a[j]+ny*k)]);
5451  }
5452  }
5453  }
5454  }
5455  else {
5456  // labels on z axis
5457  for (i=0;i<nx;i++) {
5458  for (j=0;j<ny;j++) {
5459  for (k=1;k<=n;k++) {
5460  SetBinContent(i,j,k,cont[i+nx*(j+ny*a[k])]);
5461  if (!errors.empty()) SetBinError(i,j,k,errors[i+nx*(j+ny*a[k])]);
5462  }
5463  }
5464  }
5465  }
5466  }
5467  }
5468  fEntries = entries;
5469  delete labold;
5470 }
5471 
5472 ////////////////////////////////////////////////////////////////////////////////
5473 /// Test if two double are almost equal.
5474 
5475 static inline Bool_t AlmostEqual(Double_t a, Double_t b, Double_t epsilon = 0.00000001)
5476 {
5477  return TMath::Abs(a - b) < epsilon;
5478 }
5479 
5480 ////////////////////////////////////////////////////////////////////////////////
5481 /// Test if a double is almost an integer.
5482 
5483 static inline Bool_t AlmostInteger(