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