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