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