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