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