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