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