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