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 "TBuffer.h"
22#include "TEnv.h"
23#include "TClass.h"
24#include "TMath.h"
25#include "THashList.h"
26#include "TH1.h"
27#include "TH2.h"
28#include "TH3.h"
29#include "TF2.h"
30#include "TF3.h"
31#include "TPluginManager.h"
32#include "TVirtualPad.h"
33#include "TRandom.h"
34#include "TVirtualFitter.h"
35#include "THLimitsFinder.h"
36#include "TProfile.h"
37#include "TStyle.h"
38#include "TVectorF.h"
39#include "TVectorD.h"
40#include "TBrowser.h"
41#include "TError.h"
42#include "TVirtualHistPainter.h"
43#include "TVirtualFFT.h"
44#include "TVirtualPaveStats.h"
45#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" Ignore the bin uncertainties when fitting using the default least square (chi2) method but skip empty bins
3841/// - "WW" Ignore the bin uncertainties when fitting using the default least square (chi2) method and include also the empty bins
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/// -"MULTI" Use Loglikelihood method based on multi-nomial distribution.
3848/// In this case function must be normalized and one fits only the function shape (a not extended binned
3849/// likelihood fit)
3850/// - "P" Use Pearson chi2 (using expected errors instead of observed errors)
3851/// - "U" Use a User specified fitting algorithm (via SetFCN)
3852/// - "Q" Quiet mode (minimum printing)
3853/// - "V" Verbose mode (default is between Q and V)
3854/// - "E" Perform better Errors estimation using Minos technique
3855/// - "B" User defined parameter settings are used for predefined functions
3856/// like "gaus", "expo", "poln", "landau".
3857/// Use this option when you want to fix one or more parameters for these functions.
3858/// - "M" More. Improve fit results.
3859/// It uses the IMPROVE command of TMinuit (see TMinuit::mnimpr).
3860/// This algorithm attempts to improve the found local minimum by searching for a
3861/// better one.
3862/// - "R" Use the Range specified in the function range
3863/// - "N" Do not store the graphics function, do not draw
3864/// - "0" Do not plot the result of the fit. By default the fitted function
3865/// is drawn unless the option"N" above is specified.
3866/// - "+" Add this new fitted function to the list of fitted functions
3867/// (by default, any previous function is deleted)
3868/// - "C" In case of linear fitting, don't calculate the chisquare
3869/// (saves time)
3870/// - "F" If fitting a polN, switch to minuit fitter
3871/// - "S" The result of the fit is returned in the TFitResultPtr
3872/// (see below Access to the Fit Result)
3873/// \param[in] goption specify a list of graphics options. See TH1::Draw for a complete list of these options.
3874/// \param[in] xxmin range
3875/// \param[in] xxmax range
3876///
3877/// In order to use the Range option, one must first create a function
3878/// with the expression to be fitted. For example, if your histogram
3879/// has a defined range between -4 and 4 and you want to fit a gaussian
3880/// only in the interval 1 to 3, you can do:
3881///
3882/// ~~~ {.cpp}
3883/// TF1 *f1 = new TF1("f1", "gaus", 1, 3);
3884/// histo->Fit("f1", "R");
3885/// ~~~
3886///
3887/// ## Setting initial conditions
3888/// Parameters must be initialized before invoking the Fit function.
3889/// The setting of the parameter initial values is automatic for the
3890/// predefined functions : poln, expo, gaus, landau. One can however disable
3891/// this automatic computation by specifying the option "B".
3892/// Note that if a predefined function is defined with an argument,
3893/// eg, gaus(0), expo(1), you must specify the initial values for
3894/// the parameters.
3895/// You can specify boundary limits for some or all parameters via
3896///
3897/// ~~~ {.cpp}
3898/// f1->SetParLimits(p_number, parmin, parmax);
3899/// ~~~
3900///
3901/// if parmin>=parmax, the parameter is fixed
3902/// Note that you are not forced to fix the limits for all parameters.
3903/// For example, if you fit a function with 6 parameters, you can do:
3904///
3905/// ~~~ {.cpp}
3906/// func->SetParameters(0, 3.1, 1.e-6, -8, 0, 100);
3907/// func->SetParLimits(3, -10, -4);
3908/// func->FixParameter(4, 0);
3909/// func->SetParLimits(5, 1, 1);
3910/// ~~~
3911///
3912/// With this setup, parameters 0->2 can vary freely
3913/// Parameter 3 has boundaries [-10,-4] with initial value -8
3914/// Parameter 4 is fixed to 0
3915/// Parameter 5 is fixed to 100.
3916/// When the lower limit and upper limit are equal, the parameter is fixed.
3917/// However to fix a parameter to 0, one must call the FixParameter function.
3918///
3919///
3920/// #### Changing the fitting objective function
3921///
3922/// By default a chi square function is used for fitting. When option "L" (or "LL") is used
3923/// a Poisson likelihood function (see note below) is used.
3924/// Using option "MULTI" a multinomial likelihood fit is used. In this case the function normalization is not fitted
3925/// but only the function shape. Therefore the provided function must be normalized.
3926/// The functions are defined in the header Fit/Chi2Func.h or Fit/PoissonLikelihoodFCN and they
3927/// are implemented using the routines FitUtil::EvaluateChi2 or FitUtil::EvaluatePoissonLogL in
3928/// the file math/mathcore/src/FitUtil.cxx.
3929/// To specify a User defined fitting function, specify option "U" and
3930/// call the following functions:
3931///
3932/// ~~~ {.cpp}
3933/// TVirtualFitter::Fitter(myhist)->SetFCN(MyFittingFunction)
3934/// ~~~
3935///
3936/// where MyFittingFunction is of type:
3937///
3938/// ~~~ {.cpp}
3939/// extern void MyFittingFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
3940/// ~~~
3941///
3942/// #### Chi2 Fits
3943///
3944/// By default a chi2 (least-square) fit is performed on the histogram. The so-called modified least-square method
3945/// is used where the residual for each bin is computed using as error the observed value (the bin error)
3946///
3947/// \f[
3948/// Chi2 = \sum{ \left(\frac{y(i) - f(x(i) | p )}{e(i)} \right)^2 }
3949/// \f]
3950///
3951/// 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
3952/// an un-weighted histogram. Bins with zero errors are excluded from the fit. See also later the note on the treatment
3953/// of empty bins. When using option "I" the residual is computed not using the function value at the bin center, f
3954/// (x(i) | p), but the integral of the function in the bin, Integral{ f(x|p)dx } divided by the bin volume
3955///
3956/// #### Likelihood Fits
3957///
3958/// When using option "L" a likelihood fit is used instead of the default chi2 square fit.
3959/// The likelihood is built assuming a Poisson probability density function for each bin.
3960/// The negative log-likelihood to be minimized is
3961///
3962/// \f[
3963/// NLL = \sum{ log Poisson ( y(i) | f(x(i) | p ) ) }
3964/// \f]
3965///
3966/// The exact likelihood used is the Poisson likelihood described in this paper:
3967/// S. Baker and R. D. Cousins, “Clarification of the use of chi-square and likelihood functions in fits to histograms,”
3968/// Nucl. Instrum. Meth. 221 (1984) 437.
3969///
3970/// This method can then be used only when the bin content represents counts (i.e. errors are sqrt(N) ).
3971/// The likelihood method has the advantage of treating correctly bins with low statistics. In case of high
3972/// statistics/bin the distribution of the bin content becomes a normal distribution and the likelihood and chi2 fit
3973/// give the same result.
3974///
3975/// The likelihood method, although a bit slower, it is therefore the recommended method in case of low
3976/// bin statistics, where the chi2 method may give incorrect results, in particular when there are
3977/// several empty bins (see also below).
3978/// In case of a weighted histogram, it is possible to perform a likelihood fit by using the
3979/// option "WL". Note a weighted histogram is an histogram which has been filled with weights and it
3980/// contains the sum of the weight square ( TH1::Sumw2() has been called). The bin error for a weighted
3981/// histogram is the square root of the sum of the weight square.
3982///
3983/// #### Treatment of Empty Bins
3984///
3985/// Empty bins, which have the content equal to zero AND error equal to zero,
3986/// are excluded by default from the chisquare fit, but they are considered in the likelihood fit.
3987/// since they affect the likelihood if the function value in these bins is not negligible.
3988/// When using option "WW" these bins will be considered in the chi2 fit with an error of 1.
3989/// Note that if the histogram is having bins with zero content and non zero-errors they are considered as
3990/// any other bins in the fit. Instead bins with zero error and non-zero content are excluded in the chi2 fit.
3991/// A likelihood fit should also not be performed on such an histogram, since we are assuming a wrong pdf for each bin.
3992/// In general, one should not fit an histogram with non-empty bins and zero errors, apart if all the bins have zero
3993/// errors. In this case one could use the option "w", which gives a weight=1 for each bin (unweighted least-square
3994/// fit).
3995/// Note that in case of histogram with no errors (chi2 fit with option W or W1) the resulting fitted parameter errors
3996/// are corrected by the obtained chi2 value using this expression: errorp *= sqrt(chisquare/(ndf-1))
3997///
3998/// #### Fitting a histogram of dimension N with a function of dimension N-1
3999///
4000/// It is possible to fit a TH2 with a TF1 or a TH3 with a TF2.
4001/// In this case the option "Integral" is not allowed and each cell has
4002/// equal weight. Also in this case th eobtained parameter error are corrected as in the case when the
4003/// option "W" is used (see above)
4004///
4005/// #### Associated functions
4006///
4007/// One or more object (typically a TF1*) can be added to the list
4008/// of functions (fFunctions) associated to each histogram.
4009/// When TH1::Fit is invoked, the fitted function is added to this list.
4010/// Given an histogram h, one can retrieve an associated function
4011/// with:
4012///
4013/// ~~~ {.cpp}
4014/// TF1 *myfunc = h->GetFunction("myfunc");
4015/// ~~~
4016///
4017/// #### Access to the fit result
4018///
4019/// The function returns a TFitResultPtr which can hold a pointer to a TFitResult object.
4020/// By default the TFitResultPtr contains only the status of the fit which is return by an
4021/// automatic conversion of the TFitResultPtr to an integer. One can write in this case directly:
4022///
4023/// ~~~ {.cpp}
4024/// Int_t fitStatus = h->Fit(myFunc)
4025/// ~~~
4026///
4027/// If the option "S" is instead used, TFitResultPtr contains the TFitResult and behaves as a smart
4028/// pointer to it. For example one can do:
4029///
4030/// ~~~ {.cpp}
4031/// TFitResultPtr r = h->Fit(myFunc,"S");
4032/// TMatrixDSym cov = r->GetCovarianceMatrix(); // to access the covariance matrix
4033/// Double_t chi2 = r->Chi2(); // to retrieve the fit chi2
4034/// Double_t par0 = r->Parameter(0); // retrieve the value for the parameter 0
4035/// Double_t err0 = r->ParError(0); // retrieve the error for the parameter 0
4036/// r->Print("V"); // print full information of fit including covariance matrix
4037/// r->Write(); // store the result in a file
4038/// ~~~
4039///
4040/// The fit parameters, error and chi2 (but not covariance matrix) can be retrieved also
4041/// from the fitted function.
4042/// If the histogram is made persistent, the list of
4043/// associated functions is also persistent. Given a pointer (see above)
4044/// to an associated function myfunc, one can retrieve the function/fit
4045/// parameters with calls such as:
4046///
4047/// ~~~ {.cpp}
4048/// Double_t chi2 = myfunc->GetChisquare();
4049/// Double_t par0 = myfunc->GetParameter(0); //value of 1st parameter
4050/// Double_t err0 = myfunc->GetParError(0); //error on first parameter
4051/// ~~~
4052///
4053/// #### Access to the fit status
4054///
4055/// The status of the fit can be obtained converting the TFitResultPtr to an integer
4056/// independently if the fit option "S" is used or not:
4057///
4058/// ~~~ {.cpp}
4059/// TFitResultPtr r = h->Fit(myFunc,opt);
4060/// Int_t fitStatus = r;
4061/// ~~~
4062///
4063/// The fitStatus is 0 if the fit is OK (i.e no error occurred).
4064/// The value of the fit status code is negative in case of an error not connected with the
4065/// minimization procedure, for example when a wrong function is used.
4066/// Otherwise the return value is the one returned from the minimization procedure.
4067/// When TMinuit (default case) or Minuit2 are used as minimizer the status returned is :
4068/// `fitStatus = migradResult + 10*minosResult + 100*hesseResult + 1000*improveResult`.
4069/// TMinuit will return 0 (for migrad, minos, hesse or improve) in case of success and 4 in
4070/// case of error (see the documentation of TMinuit::mnexcm). So for example, for an error
4071/// only in Minos but not in Migrad a fitStatus of 40 will be returned.
4072/// Minuit2 will return also 0 in case of success and different values in migrad minos or
4073/// hesse depending on the error. See in this case the documentation of
4074/// Minuit2Minimizer::Minimize for the migradResult, Minuit2Minimizer::GetMinosError for the
4075/// minosResult and Minuit2Minimizer::Hesse for the hesseResult.
4076/// If other minimizers are used see their specific documentation for the status code returned.
4077/// For example in the case of Fumili, for the status returned see TFumili::Minimize.
4078///
4079/// #### Excluding points
4080///
4081/// Use TF1::RejectPoint inside your fitting function to exclude points
4082/// within a certain range from the fit. Example:
4083///
4084/// ~~~ {.cpp}
4085/// Double_t fline(Double_t *x, Double_t *par)
4086/// {
4087/// if (x[0] > 2.5 && x[0] < 3.5) {
4088/// TF1::RejectPoint();
4089/// return 0;
4090/// }
4091/// return par[0] + par[1]*x[0];
4092/// }
4093///
4094/// void exclude() {
4095/// TF1 *f1 = new TF1("f1", "[0] +[1]*x +gaus(2)", 0, 5);
4096/// f1->SetParameters(6, -1,5, 3, 0.2);
4097/// TH1F *h = new TH1F("h", "background + signal", 100, 0, 5);
4098/// h->FillRandom("f1", 2000);
4099/// TF1 *fline = new TF1("fline", fline, 0, 5, 2);
4100/// fline->SetParameters(2, -1);
4101/// h->Fit("fline", "l");
4102/// }
4103/// ~~~
4104///
4105/// #### Warning when using the option "0"
4106///
4107/// When selecting the option "0", the fitted function is added to
4108/// the list of functions of the histogram, but it is not drawn.
4109/// You can undo what you disabled in the following way:
4110///
4111/// ~~~ {.cpp}
4112/// h.Fit("myFunction", "0"); // fit, store function but do not draw
4113/// h.Draw(); function is not drawn
4114/// const Int_t kNotDraw = 1<<9;
4115/// h.GetFunction("myFunction")->ResetBit(kNotDraw);
4116/// h.Draw(); // function is visible again
4117/// ~~~
4118///
4119/// #### Access to the Minimizer information during fitting
4120///
4121/// This function calls, the ROOT::Fit::FitObject function implemented in HFitImpl.cxx
4122/// which uses the ROOT::Fit::Fitter class. The Fitter class creates the objective function
4123/// (e.g. chi2 or likelihood) and uses an implementation of the Minimizer interface for minimizing
4124/// the function.
4125/// The default minimizer is Minuit (class TMinuitMinimizer which calls TMinuit).
4126/// The default can be set in the resource file in etc/system.rootrc. For example
4127///
4128/// ~~~ {.cpp}
4129/// Root.Fitter: Minuit2
4130/// ~~~
4131///
4132/// A different fitter can also be set via ROOT::Math::MinimizerOptions::SetDefaultMinimizer
4133/// (or TVirtualFitter::SetDefaultFitter).
4134/// For example ROOT::Math::MinimizerOptions::SetDefaultMinimizer("GSLMultiMin","BFGS");
4135/// will set the usage of the BFGS algorithm of the GSL multi-dimensional minimization
4136/// (implemented in libMathMore). ROOT::Math::MinimizerOptions can be used also to set other
4137/// default options, like maximum number of function calls, minimization tolerance or print
4138/// level. See the documentation of this class.
4139///
4140/// For fitting linear functions (containing the "++" sign" and polN functions,
4141/// the linear fitter is automatically initialized.
4142
4143TFitResultPtr TH1::Fit(TF1 *f1 ,Option_t *option ,Option_t *goption, Double_t xxmin, Double_t xxmax)
4144{
4145 // implementation of Fit method is in file hist/src/HFitImpl.cxx
4146 Foption_t fitOption;
4148
4149 // create range and minimizer options with default values
4150 ROOT::Fit::DataRange range(xxmin,xxmax);
4152
4153 // need to empty the buffer before
4154 // (t.b.d. do a ML unbinned fit with buffer data)
4155 if (fBuffer) BufferEmpty();
4156
4157 return ROOT::Fit::FitObject(this, f1 , fitOption , minOption, goption, range);
4158}
4159
4160////////////////////////////////////////////////////////////////////////////////
4161/// Display a panel with all histogram fit options.
4162///
4163/// See class TFitPanel for example
4164
4165void TH1::FitPanel()
4166{
4167 if (!gPad)
4168 gROOT->MakeDefCanvas();
4169
4170 if (!gPad) {
4171 Error("FitPanel", "Unable to create a default canvas");
4172 return;
4173 }
4174
4175
4176 // use plugin manager to create instance of TFitEditor
4177 TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TFitEditor");
4178 if (handler && handler->LoadPlugin() != -1) {
4179 if (handler->ExecPlugin(2, gPad, this) == 0)
4180 Error("FitPanel", "Unable to create the FitPanel");
4181 }
4182 else
4183 Error("FitPanel", "Unable to find the FitPanel plug-in");
4184}
4185
4186////////////////////////////////////////////////////////////////////////////////
4187/// Return an histogram containing the asymmetry of this histogram with h2,
4188/// where the asymmetry is defined as:
4189///
4190/// ~~~ {.cpp}
4191/// Asymmetry = (h1 - h2)/(h1 + h2) where h1 = this
4192/// ~~~
4193///
4194/// works for 1D, 2D, etc. histograms
4195/// c2 is an optional argument that gives a relative weight between the two
4196/// histograms, and dc2 is the error on this weight. This is useful, for example,
4197/// when forming an asymmetry between two histograms from 2 different data sets that
4198/// need to be normalized to each other in some way. The function calculates
4199/// the errors assuming Poisson statistics on h1 and h2 (that is, dh = sqrt(h)).
4200///
4201/// example: assuming 'h1' and 'h2' are already filled
4202///
4203/// ~~~ {.cpp}
4204/// h3 = h1->GetAsymmetry(h2)
4205/// ~~~
4206///
4207/// then 'h3' is created and filled with the asymmetry between 'h1' and 'h2';
4208/// h1 and h2 are left intact.
4209///
4210/// Note that it is the user's responsibility to manage the created histogram.
4211/// The name of the returned histogram will be `Asymmetry_nameOfh1-nameOfh2`
4212///
4213/// code proposed by Jason Seely (seely@mit.edu) and adapted by R.Brun
4214///
4215/// clone the histograms so top and bottom will have the
4216/// correct dimensions:
4217/// Sumw2 just makes sure the errors will be computed properly
4218/// when we form sums and ratios below.
4219
4221{
4222 TH1 *h1 = this;
4223 TString name = TString::Format("Asymmetry_%s-%s",h1->GetName(),h2->GetName() );
4224 TH1 *asym = (TH1*)Clone(name);
4225
4226 // set also the title
4227 TString title = TString::Format("(%s - %s)/(%s+%s)",h1->GetName(),h2->GetName(),h1->GetName(),h2->GetName() );
4228 asym->SetTitle(title);
4229
4230 asym->Sumw2();
4231 Bool_t addStatus = TH1::AddDirectoryStatus();
4233 TH1 *top = (TH1*)asym->Clone();
4234 TH1 *bottom = (TH1*)asym->Clone();
4235 TH1::AddDirectory(addStatus);
4236
4237 // form the top and bottom of the asymmetry, and then divide:
4238 top->Add(h1,h2,1,-c2);
4239 bottom->Add(h1,h2,1,c2);
4240 asym->Divide(top,bottom);
4241
4242 Int_t xmax = asym->GetNbinsX();
4243 Int_t ymax = asym->GetNbinsY();
4244 Int_t zmax = asym->GetNbinsZ();
4245
4246 if (h1->fBuffer) h1->BufferEmpty(1);
4247 if (h2->fBuffer) h2->BufferEmpty(1);
4248 if (bottom->fBuffer) bottom->BufferEmpty(1);
4249
4250 // now loop over bins to calculate the correct errors
4251 // the reason this error calculation looks complex is because of c2
4252 for(Int_t i=1; i<= xmax; i++){
4253 for(Int_t j=1; j<= ymax; j++){
4254 for(Int_t k=1; k<= zmax; k++){
4255 Int_t bin = GetBin(i, j, k);
4256 // here some bin contents are written into variables to make the error
4257 // calculation a little more legible:
4259 Double_t b = h2->RetrieveBinContent(bin);
4260 Double_t bot = bottom->RetrieveBinContent(bin);
4261
4262 // make sure there are some events, if not, then the errors are set = 0
4263 // automatically.
4264 //if(bot < 1){} was changed to the next line from recommendation of Jason Seely (28 Nov 2005)
4265 if(bot < 1e-6){}
4266 else{
4267 // computation of errors by Christos Leonidopoulos
4268 Double_t dasq = h1->GetBinErrorSqUnchecked(bin);
4269 Double_t dbsq = h2->GetBinErrorSqUnchecked(bin);
4270 Double_t error = 2*TMath::Sqrt(a*a*c2*c2*dbsq + c2*c2*b*b*dasq+a*a*b*b*dc2*dc2)/(bot*bot);
4271 asym->SetBinError(i,j,k,error);
4272 }
4273 }
4274 }
4275 }
4276 delete top;
4277 delete bottom;
4278
4279 return asym;
4280}
4281
4282////////////////////////////////////////////////////////////////////////////////
4283/// Static function
4284/// return the default buffer size for automatic histograms
4285/// the parameter fgBufferSize may be changed via SetDefaultBufferSize
4286
4288{
4289 return fgBufferSize;
4290}
4291
4292////////////////////////////////////////////////////////////////////////////////
4293/// Return kTRUE if TH1::Sumw2 must be called when creating new histograms.
4294/// see TH1::SetDefaultSumw2.
4295
4297{
4298 return fgDefaultSumw2;
4299}
4300
4301////////////////////////////////////////////////////////////////////////////////
4302/// Return the current number of entries.
4303
4305{
4306 if (fBuffer) {
4307 Int_t nentries = (Int_t) fBuffer[0];
4308 if (nentries > 0) return nentries;
4309 }
4310
4311 return fEntries;
4312}
4313
4314////////////////////////////////////////////////////////////////////////////////
4315/// Number of effective entries of the histogram.
4316///
4317/// \f[
4318/// neff = \frac{(\sum Weights )^2}{(\sum Weight^2 )}
4319/// \f]
4320///
4321/// In case of an unweighted histogram this number is equivalent to the
4322/// number of entries of the histogram.
4323/// For a weighted histogram, this number corresponds to the hypothetical number of unweighted entries
4324/// a histogram would need to have the same statistical power as this weighted histogram.
4325/// Note: The underflow/overflow are included if one has set the TH1::StatOverFlows flag
4326/// and if the statistics has been computed at filling time.
4327/// If a range is set in the histogram the number is computed from the given range.
4328
4330{
4331 Stat_t s[kNstat];
4332 this->GetStats(s);// s[1] sum of squares of weights, s[0] sum of weights
4333 return (s[1] ? s[0]*s[0]/s[1] : TMath::Abs(s[0]) );
4334}
4335
4336////////////////////////////////////////////////////////////////////////////////
4337/// Set highlight (enable/disable) mode for the histogram
4338/// by default highlight mode is disable
4339
4340void TH1::SetHighlight(Bool_t set)
4341{
4342 if (IsHighlight() == set) return;
4343 if (fDimension > 2) {
4344 Info("SetHighlight", "Supported only 1-D or 2-D histograms");
4345 return;
4346 }
4347
4348 if (!fPainter) {
4349 Info("SetHighlight", "Need to draw histogram first");
4350 return;
4351 }
4352 SetBit(kIsHighlight, set);
4354}
4355
4356////////////////////////////////////////////////////////////////////////////////
4357/// Redefines TObject::GetObjectInfo.
4358/// Displays the histogram info (bin number, contents, integral up to bin
4359/// corresponding to cursor position px,py
4360
4361char *TH1::GetObjectInfo(Int_t px, Int_t py) const
4362{
4363 return ((TH1*)this)->GetPainter()->GetObjectInfo(px,py);
4364}
4365
4366////////////////////////////////////////////////////////////////////////////////
4367/// Return pointer to painter.
4368/// If painter does not exist, it is created
4369
4371{
4372 if (!fPainter) {
4373 TString opt = option;
4374 opt.ToLower();
4375 if (opt.Contains("gl") || gStyle->GetCanvasPreferGL()) {
4376 //try to create TGLHistPainter
4377 TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TGLHistPainter");
4378
4379 if (handler && handler->LoadPlugin() != -1)
4380 fPainter = reinterpret_cast<TVirtualHistPainter *>(handler->ExecPlugin(1, this));
4381 }
4382 }
4383
4385
4386 return fPainter;
4387}
4388
4389////////////////////////////////////////////////////////////////////////////////
4390/// Compute Quantiles for this histogram
4391/// Quantile x_q of a probability distribution Function F is defined as
4392///
4393/// ~~~ {.cpp}
4394/// F(x_q) = q with 0 <= q <= 1.
4395/// ~~~
4396///
4397/// For instance the median x_0.5 of a distribution is defined as that value
4398/// of the random variable for which the distribution function equals 0.5:
4399///
4400/// ~~~ {.cpp}
4401/// F(x_0.5) = Probability(x < x_0.5) = 0.5
4402/// ~~~
4403///
4404/// code from Eddy Offermann, Renaissance
4405///
4406/// \param[in] nprobSum maximum size of array q and size of array probSum (if given)
4407/// \param[in] probSum array of positions where quantiles will be computed.
4408/// - if probSum is null, probSum will be computed internally and will
4409/// have a size = number of bins + 1 in h. it will correspond to the
4410/// quantiles calculated at the lowest edge of the histogram (quantile=0) and
4411/// all the upper edges of the bins.
4412/// - if probSum is not null, it is assumed to contain at least nprobSum values.
4413/// \param[out] q array q filled with nq quantiles
4414/// \return value nq (<=nprobSum) with the number of quantiles computed
4415///
4416/// Note that the Integral of the histogram is automatically recomputed
4417/// if the number of entries is different of the number of entries when
4418/// the integral was computed last time. In case you do not use the Fill
4419/// functions to fill your histogram, but SetBinContent, you must call
4420/// TH1::ComputeIntegral before calling this function.
4421///
4422/// Getting quantiles q from two histograms and storing results in a TGraph,
4423/// a so-called QQ-plot
4424///
4425/// ~~~ {.cpp}
4426/// TGraph *gr = new TGraph(nprob);
4427/// h1->GetQuantiles(nprob,gr->GetX());
4428/// h2->GetQuantiles(nprob,gr->GetY());
4429/// gr->Draw("alp");
4430/// ~~~
4431///
4432/// Example:
4433///
4434/// ~~~ {.cpp}
4435/// void quantiles() {
4436/// // demo for quantiles
4437/// const Int_t nq = 20;
4438/// TH1F *h = new TH1F("h","demo quantiles",100,-3,3);
4439/// h->FillRandom("gaus",5000);
4440///
4441/// Double_t xq[nq]; // position where to compute the quantiles in [0,1]
4442/// Double_t yq[nq]; // array to contain the quantiles
4443/// for (Int_t i=0;i<nq;i++) xq[i] = Float_t(i+1)/nq;
4444/// h->GetQuantiles(nq,yq,xq);
4445///
4446/// //show the original histogram in the top pad
4447/// TCanvas *c1 = new TCanvas("c1","demo quantiles",10,10,700,900);
4448/// c1->Divide(1,2);
4449/// c1->cd(1);
4450/// h->Draw();
4451///
4452/// // show the quantiles in the bottom pad
4453/// c1->cd(2);
4454/// gPad->SetGrid();
4455/// TGraph *gr = new TGraph(nq,xq,yq);
4456/// gr->SetMarkerStyle(21);
4457/// gr->Draw("alp");
4458/// }
4459/// ~~~
4460
4461Int_t TH1::GetQuantiles(Int_t nprobSum, Double_t *q, const Double_t *probSum)
4462{
4463 if (GetDimension() > 1) {
4464 Error("GetQuantiles","Only available for 1-d histograms");
4465 return 0;
4466 }
4467
4468 const Int_t nbins = GetXaxis()->GetNbins();
4469 if (!fIntegral) ComputeIntegral();
4470 if (fIntegral[nbins+1] != fEntries) ComputeIntegral();
4471
4472 Int_t i, ibin;
4473 Double_t *prob = (Double_t*)probSum;
4474 Int_t nq = nprobSum;
4475 if (probSum == 0) {
4476 nq = nbins+1;
4477 prob = new Double_t[nq];
4478 prob[0] = 0;
4479 for (i=1;i<nq;i++) {
4480 prob[i] = fIntegral[i]/fIntegral[nbins];
4481 }
4482 }
4483
4484 for (i = 0; i < nq; i++) {
4485 ibin = TMath::BinarySearch(nbins,fIntegral,prob[i]);
4486 while (ibin < nbins-1 && fIntegral[ibin+1] == prob[i]) {
4487 if (fIntegral[ibin+2] == prob[i]) ibin++;
4488 else break;
4489 }
4490 q[i] = GetBinLowEdge(ibin+1);
4491 const Double_t dint = fIntegral[ibin+1]-fIntegral[ibin];
4492 if (dint > 0) q[i] += GetBinWidth(ibin+1)*(prob[i]-fIntegral[ibin])/dint;
4493 }
4494
4495 if (!probSum) delete [] prob;
4496 return nq;
4497}
4498
4499////////////////////////////////////////////////////////////////////////////////
4500/// Decode string choptin and fill fitOption structure.
4501
4502Int_t TH1::FitOptionsMake(Option_t *choptin, Foption_t &fitOption)
4503{
4505 return 1;
4506}
4507
4508////////////////////////////////////////////////////////////////////////////////
4509/// Compute Initial values of parameters for a gaussian.
4510
4511void H1InitGaus()
4512{
4513 Double_t allcha, sumx, sumx2, x, val, stddev, mean;
4514 Int_t bin;
4515 const Double_t sqrtpi = 2.506628;
4516
4517 // - Compute mean value and StdDev of the histogram in the given range
4519 TH1 *curHist = (TH1*)hFitter->GetObjectFit();
4520 Int_t hxfirst = hFitter->GetXfirst();
4521 Int_t hxlast = hFitter->GetXlast();
4522 Double_t valmax = curHist->GetBinContent(hxfirst);
4523 Double_t binwidx = curHist->GetBinWidth(hxfirst);
4524 allcha = sumx = sumx2 = 0;
4525 for (bin=hxfirst;bin<=hxlast;bin++) {
4526 x = curHist->GetBinCenter(bin);
4527 val = TMath::Abs(curHist->GetBinContent(bin));
4528 if (val > valmax) valmax = val;
4529 sumx += val*x;
4530 sumx2 += val*x*x;
4531 allcha += val;
4532 }
4533 if (allcha == 0) return;
4534 mean = sumx/allcha;
4535 stddev = sumx2/allcha - mean*mean;
4536 if (stddev > 0) stddev = TMath::Sqrt(stddev);
4537 else stddev = 0;
4538 if (stddev == 0) stddev = binwidx*(hxlast-hxfirst+1)/4;
4539 //if the distribution is really gaussian, the best approximation
4540 //is binwidx*allcha/(sqrtpi*stddev)
4541 //However, in case of non-gaussian tails, this underestimates
4542 //the normalisation constant. In this case the maximum value
4543 //is a better approximation.
4544 //We take the average of both quantities
4545 Double_t constant = 0.5*(valmax+binwidx*allcha/(sqrtpi*stddev));
4546
4547 //In case the mean value is outside the histo limits and
4548 //the StdDev is bigger than the range, we take
4549 // mean = center of bins
4550 // stddev = half range
4551 Double_t xmin = curHist->GetXaxis()->GetXmin();
4552 Double_t xmax = curHist->GetXaxis()->GetXmax();
4553 if ((mean < xmin || mean > xmax) && stddev > (xmax-xmin)) {
4554 mean = 0.5*(xmax+xmin);
4555 stddev = 0.5*(xmax-xmin);
4556 }
4557 TF1 *f1 = (TF1*)hFitter->GetUserFunc();
4558 f1->SetParameter(0,constant);
4559 f1->SetParameter(1,mean);
4560 f1->SetParameter(2,stddev);
4561 f1->SetParLimits(2,0,10*stddev);
4562}
4563
4564////////////////////////////////////////////////////////////////////////////////
4565/// Compute Initial values of parameters for an exponential.
4566
4567void H1InitExpo()
4568{
4569 Double_t constant, slope;
4570 Int_t ifail;
4572 Int_t hxfirst = hFitter->GetXfirst();
4573 Int_t hxlast = hFitter->GetXlast();
4574 Int_t nchanx = hxlast - hxfirst + 1;
4575
4576 H1LeastSquareLinearFit(-nchanx, constant, slope, ifail);
4577
4578 TF1 *f1 = (TF1*)hFitter->GetUserFunc();
4579 f1->SetParameter(0,constant);
4580 f1->SetParameter(1,slope);
4581
4582}
4583
4584////////////////////////////////////////////////////////////////////////////////
4585/// Compute Initial values of parameters for a polynom.
4586
4587void H1InitPolynom()
4588{
4589 Double_t fitpar[25];
4590
4592 TF1 *f1 = (TF1*)hFitter->GetUserFunc();
4593 Int_t hxfirst = hFitter->GetXfirst();
4594 Int_t hxlast = hFitter->GetXlast();
4595 Int_t nchanx = hxlast - hxfirst + 1;
4596 Int_t npar = f1->GetNpar();
4597
4598 if (nchanx <=1 || npar == 1) {
4599 TH1 *curHist = (TH1*)hFitter->GetObjectFit();
4600 fitpar[0] = curHist->GetSumOfWeights()/Double_t(nchanx);
4601 } else {
4602 H1LeastSquareFit( nchanx, npar, fitpar);
4603 }
4604 for (Int_t i=0;i<npar;i++) f1->SetParameter(i, fitpar[i]);
4605}
4606
4607////////////////////////////////////////////////////////////////////////////////
4608/// Least squares lpolynomial fitting without weights.
4609///
4610/// \param[in] n number of points to fit
4611/// \param[in] m number of parameters
4612/// \param[in] a array of parameters
4613///
4614/// based on CERNLIB routine LSQ: Translated to C++ by Rene Brun
4615/// (E.Keil. revised by B.Schorr, 23.10.1981.)
4616
4618{
4619 const Double_t zero = 0.;
4620 const Double_t one = 1.;
4621 const Int_t idim = 20;
4622
4623 Double_t b[400] /* was [20][20] */;
4624 Int_t i, k, l, ifail;
4625 Double_t power;
4626 Double_t da[20], xk, yk;
4627
4628 if (m <= 2) {
4629 H1LeastSquareLinearFit(n, a[0], a[1], ifail);
4630 return;
4631 }
4632 if (m > idim || m > n) return;
4633 b[0] = Double_t(n);
4634 da[0] = zero;
4635 for (l = 2; l <= m; ++l) {
4636 b[l-1] = zero;
4637 b[m + l*20 - 21] = zero;
4638 da[l-1] = zero;
4639 }
4641 TH1 *curHist = (TH1*)hFitter->GetObjectFit();
4642 Int_t hxfirst = hFitter->GetXfirst();
4643 Int_t hxlast = hFitter->GetXlast();
4644 for (k = hxfirst; k <= hxlast; ++k) {
4645 xk = curHist->GetBinCenter(k);
4646 yk = curHist->GetBinContent(k);
4647 power = one;
4648 da[0] += yk;
4649 for (l = 2; l <= m; ++l) {
4650 power *= xk;
4651 b[l-1] += power;
4652 da[l-1] += power*yk;
4653 }
4654 for (l = 2; l <= m; ++l) {
4655 power *= xk;
4656 b[m + l*20 - 21] += power;
4657 }
4658 }
4659 for (i = 3; i <= m; ++i) {
4660 for (k = i; k <= m; ++k) {
4661 b[k - 1 + (i-1)*20 - 21] = b[k + (i-2)*20 - 21];
4662 }
4663 }
4664 H1LeastSquareSeqnd(m, b, idim, ifail, 1, da);
4665
4666 for (i=0; i<m; ++i) a[i] = da[i];
4667
4668}
4669
4670////////////////////////////////////////////////////////////////////////////////
4671/// Least square linear fit without weights.
4672///
4673/// extracted from CERNLIB LLSQ: Translated to C++ by Rene Brun
4674/// (added to LSQ by B. Schorr, 15.02.1982.)
4675
4676void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail)
4677{
4678 Double_t xbar, ybar, x2bar;
4679 Int_t i, n;
4680 Double_t xybar;
4681 Double_t fn, xk, yk;
4682 Double_t det;
4683
4684 n = TMath::Abs(ndata);
4685 ifail = -2;
4686 xbar = ybar = x2bar = xybar = 0;
4688 TH1 *curHist = (TH1*)hFitter->GetObjectFit();
4689 Int_t hxfirst = hFitter->GetXfirst();
4690 Int_t hxlast = hFitter->GetXlast();
4691 for (i = hxfirst; i <= hxlast; ++i) {
4692 xk = curHist->GetBinCenter(i);
4693 yk = curHist->GetBinContent(i);
4694 if (ndata < 0) {
4695 if (yk <= 0) yk = 1e-9;
4696 yk = TMath::Log(yk);
4697 }
4698 xbar += xk;
4699 ybar += yk;
4700 x2bar += xk*xk;
4701 xybar += xk*yk;
4702 }
4703 fn = Double_t(n);
4704 det = fn*x2bar - xbar*xbar;
4705 ifail = -1;
4706 if (det <= 0) {
4707 a0 = ybar/fn;
4708 a1 = 0;
4709 return;
4710 }
4711 ifail = 0;
4712 a0 = (x2bar*ybar - xbar*xybar) / det;
4713 a1 = (fn*xybar - xbar*ybar) / det;
4714
4715}
4716
4717////////////////////////////////////////////////////////////////////////////////
4718/// Extracted from CERN Program library routine DSEQN.
4719///
4720/// Translated to C++ by Rene Brun
4721
4722void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b)
4723{
4724 Int_t a_dim1, a_offset, b_dim1, b_offset;
4725 Int_t nmjp1, i, j, l;
4726 Int_t im1, jp1, nm1, nmi;
4727 Double_t s1, s21, s22;
4728 const Double_t one = 1.;
4729
4730 /* Parameter adjustments */
4731 b_dim1 = idim;
4732 b_offset = b_dim1 + 1;
4733 b -= b_offset;
4734 a_dim1 = idim;
4735 a_offset = a_dim1 + 1;
4736 a -= a_offset;
4737
4738 if (idim < n) return;
4739
4740 ifail = 0;
4741 for (j = 1; j <= n; ++j) {
4742 if (a[j + j*a_dim1] <= 0) { ifail = -1; return; }
4743 a[j + j*a_dim1] = one / a[j + j*a_dim1];
4744 if (j == n) continue;
4745 jp1 = j + 1;
4746 for (l = jp1; l <= n; ++l) {
4747 a[j + l*a_dim1] = a[j + j*a_dim1] * a[l + j*a_dim1];
4748 s1 = -a[l + (j+1)*a_dim1];
4749 for (i = 1; i <= j; ++i) { s1 = a[l + i*a_dim1] * a[i + (j+1)*a_dim1] + s1; }
4750 a[l + (j+1)*a_dim1] = -s1;
4751 }
4752 }
4753 if (k <= 0) return;
4754
4755 for (l = 1; l <= k; ++l) {
4756 b[l*b_dim1 + 1] = a[a_dim1 + 1]*b[l*b_dim1 + 1];
4757 }
4758 if (n == 1) return;
4759 for (l = 1; l <= k; ++l) {
4760 for (i = 2; i <= n; ++i) {
4761 im1 = i - 1;
4762 s21 = -b[i + l*b_dim1];
4763 for (j = 1; j <= im1; ++j) {
4764 s21 = a[i + j*a_dim1]*b[j + l*b_dim1] + s21;
4765 }
4766 b[i + l*b_dim1] = -a[i + i*a_dim1]*s21;
4767 }
4768 nm1 = n - 1;
4769 for (i = 1; i <= nm1; ++i) {
4770 nmi = n - i;
4771 s22 = -b[nmi + l*b_dim1];
4772 for (j = 1; j <= i; ++j) {
4773 nmjp1 = n - j + 1;
4774 s22 = a[nmi + nmjp1*a_dim1]*b[nmjp1 + l*b_dim1] + s22;
4775 }
4776 b[nmi + l*b_dim1] = -s22;
4777 }
4778 }
4779}
4780
4781////////////////////////////////////////////////////////////////////////////////
4782/// Return Global bin number corresponding to binx,y,z.
4783///
4784/// 2-D and 3-D histograms are represented with a one dimensional
4785/// structure.
4786/// This has the advantage that all existing functions, such as
4787/// GetBinContent, GetBinError, GetBinFunction work for all dimensions.
4788///
4789/// In case of a TH1x, returns binx directly.
4790/// see TH1::GetBinXYZ for the inverse transformation.
4791///
4792/// Convention for numbering bins
4793///
4794/// For all histogram types: nbins, xlow, xup
4795///
4796/// - bin = 0; underflow bin
4797/// - bin = 1; first bin with low-edge xlow INCLUDED
4798/// - bin = nbins; last bin with upper-edge xup EXCLUDED
4799/// - bin = nbins+1; overflow bin
4800///
4801/// In case of 2-D or 3-D histograms, a "global bin" number is defined.
4802/// For example, assuming a 3-D histogram with binx,biny,binz, the function
4803///
4804/// ~~~ {.cpp}
4805/// Int_t bin = h->GetBin(binx,biny,binz);
4806/// ~~~
4807///
4808/// returns a global/linearized bin number. This global bin is useful
4809/// to access the bin information independently of the dimension.
4810
4811Int_t TH1::GetBin(Int_t binx, Int_t, Int_t) const
4812{
4813 Int_t ofx = fXaxis.GetNbins() + 1; // overflow bin
4814 if (binx < 0) binx = 0;
4815 if (binx > ofx) binx = ofx;
4816
4817 return binx;
4818}
4819
4820////////////////////////////////////////////////////////////////////////////////
4821/// Return binx, biny, binz corresponding to the global bin number globalbin
4822/// see TH1::GetBin function above
4823
4824void TH1::GetBinXYZ(Int_t binglobal, Int_t &binx, Int_t &biny, Int_t &binz) const
4825{
4826 Int_t nx = fXaxis.GetNbins()+2;
4827 Int_t ny = fYaxis.GetNbins()+2;
4828
4829 if (GetDimension() == 1) {
4830 binx = binglobal%nx;
4831 biny = 0;
4832 binz = 0;
4833 return;
4834 }
4835 if (GetDimension() == 2) {
4836 binx = binglobal%nx;
4837 biny = ((binglobal-binx)/nx)%ny;
4838 binz = 0;
4839 return;
4840 }
4841 if (GetDimension() == 3) {
4842 binx = binglobal%nx;
4843 biny = ((binglobal-binx)/nx)%ny;
4844 binz = ((binglobal-binx)/nx -biny)/ny;
4845 }
4846}
4847
4848////////////////////////////////////////////////////////////////////////////////
4849/// Return a random number distributed according the histogram bin contents.
4850/// This function checks if the bins integral exists. If not, the integral
4851/// is evaluated, normalized to one.
4852///
4853/// The integral is automatically recomputed if the number of entries
4854/// is not the same then when the integral was computed.
4855/// NB Only valid for 1-d histograms. Use GetRandom2 or 3 otherwise.
4856/// If the histogram has a bin with negative content a NaN is returned
4857
4859{
4860 if (fDimension > 1) {
4861 Error("GetRandom","Function only valid for 1-d histograms");
4862 return 0;
4863 }
4864 Int_t nbinsx = GetNbinsX();
4865 Double_t integral = 0;
4866 // compute integral checking that all bins have positive content (see ROOT-5894)
4867 if (fIntegral) {
4868 if (fIntegral[nbinsx+1] != fEntries) integral = ((TH1*)this)->ComputeIntegral(true);
4869 else integral = fIntegral[nbinsx];
4870 } else {
4871 integral = ((TH1*)this)->ComputeIntegral(true);
4872 }
4873 if (integral == 0) return 0;
4874 // return a NaN in case some bins have negative content
4875 if (integral == TMath::QuietNaN() ) return TMath::QuietNaN();
4876
4877 Double_t r1 = gRandom->Rndm();
4878 Int_t ibin = TMath::BinarySearch(nbinsx,fIntegral,r1);
4879 Double_t x = GetBinLowEdge(ibin+1);
4880 if (r1 > fIntegral[ibin]) x +=
4881 GetBinWidth(ibin+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1] - fIntegral[ibin]);
4882 return x;
4883}
4884
4885////////////////////////////////////////////////////////////////////////////////
4886/// Return content of bin number bin.
4887///
4888/// Implemented in TH1C,S,F,D
4889///
4890/// Convention for numbering bins
4891///
4892/// For all histogram types: nbins, xlow, xup
4893///
4894/// - bin = 0; underflow bin
4895/// - bin = 1; first bin with low-edge xlow INCLUDED
4896/// - bin = nbins; last bin with upper-edge xup EXCLUDED
4897/// - bin = nbins+1; overflow bin
4898///
4899/// In case of 2-D or 3-D histograms, a "global bin" number is defined.
4900/// For example, assuming a 3-D histogram with binx,biny,binz, the function
4901///
4902/// ~~~ {.cpp}
4903/// Int_t bin = h->GetBin(binx,biny,binz);
4904/// ~~~
4905///
4906/// returns a global/linearized bin number. This global bin is useful
4907/// to access the bin information independently of the dimension.
4908
4910{
4911 if (fBuffer) const_cast<TH1*>(this)->BufferEmpty();
4912 if (bin < 0) bin = 0;
4913 if (bin >= fNcells) bin = fNcells-1;
4914
4915 return RetrieveBinContent(bin);
4916}
4917
4918////////////////////////////////////////////////////////////////////////////////
4919/// Compute first binx in the range [firstx,lastx] for which
4920/// diff = abs(bin_content-c) <= maxdiff
4921///
4922/// In case several bins in the specified range with diff=0 are found
4923/// the first bin found is returned in binx.
4924/// In case several bins in the specified range satisfy diff <=maxdiff
4925/// the bin with the smallest difference is returned in binx.
4926/// In all cases the function returns the smallest difference.
4927///
4928/// NOTE1: if firstx <= 0, firstx is set to bin 1
4929/// if (lastx < firstx then firstx is set to the number of bins
4930/// ie if firstx=0 and lastx=0 (default) the search is on all bins.
4931///
4932/// NOTE2: if maxdiff=0 (default), the first bin with content=c is returned.
4933
4934Double_t TH1::GetBinWithContent(Double_t c, Int_t &binx, Int_t firstx, Int_t lastx,Double_t maxdiff) const
4935{
4936 if (fDimension > 1) {
4937 binx = 0;
4938 Error("GetBinWithContent","function is only valid for 1-D histograms");
4939 return 0;
4940 }
4941
4942 if (fBuffer) ((TH1*)this)->BufferEmpty();
4943
4944 if (firstx <= 0) firstx = 1;
4945 if (lastx < firstx) lastx = fXaxis.GetNbins();
4946 Int_t binminx = 0;
4947 Double_t diff, curmax = 1.e240;
4948 for (Int_t i=firstx;i<=lastx;i++) {
4949 diff = TMath::Abs(RetrieveBinContent(i)-c);
4950 if (diff <= 0) {binx = i; return diff;}
4951 if (diff < curmax && diff <= maxdiff) {curmax = diff, binminx=i;}
4952 }
4953 binx = binminx;
4954 return curmax;
4955}
4956
4957////////////////////////////////////////////////////////////////////////////////
4958/// Given a point x, approximates the value via linear interpolation
4959/// based on the two nearest bin centers
4960///
4961/// Andy Mastbaum 10/21/08
4962
4964{
4965 if (fBuffer) ((TH1*)this)->BufferEmpty();
4966
4967 Int_t xbin = fXaxis.FindFixBin(x);
4968 Double_t x0,x1,y0,y1;
4969
4970 if(x<=GetBinCenter(1)) {
4971 return RetrieveBinContent(1);
4972 } else if(x>=GetBinCenter(GetNbinsX())) {
4973 return RetrieveBinContent(GetNbinsX());
4974 } else {
4975 if(x<=GetBinCenter(xbin)) {
4976 y0 = RetrieveBinContent(xbin-1);
4977 x0 = GetBinCenter(xbin-1);
4978 y1 = RetrieveBinContent(xbin);
4979 x1 = GetBinCenter(xbin);
4980 } else {
4981 y0 = RetrieveBinContent(xbin);
4982 x0 = GetBinCenter(xbin);
4983 y1 = RetrieveBinContent(xbin+1);
4984 x1 = GetBinCenter(xbin+1);
4985 }
4986 return y0 + (x-x0)*((y1-y0)/(x1-x0));
4987 }
4988}
4989
4990////////////////////////////////////////////////////////////////////////////////
4991/// 2d Interpolation. Not yet implemented.
4992
4994{
4995 Error("Interpolate","This function must be called with 1 argument for a TH1");
4996 return 0;
4997}
4998
4999////////////////////////////////////////////////////////////////////////////////
5000/// 3d Interpolation. Not yet implemented.
5001
5003{
5004 Error("Interpolate","This function must be called with 1 argument for a TH1");
5005 return 0;
5006}
5007
5008///////////////////////////////////////////////////////////////////////////////
5009/// Check if an histogram is empty
5010/// (this a protected method used mainly by TH1Merger )
5011
5012Bool_t TH1::IsEmpty() const
5013{
5014 // if fTsumw or fentries are not zero histogram is not empty
5015 // need to use GetEntries() instead of fEntries in case of bugger histograms
5016 // so we will flash the buffer
5017 if (fTsumw != 0) return kFALSE;
5018 if (GetEntries() != 0) return kFALSE;
5019 // case fTSumw == 0 amd entries are also zero
5020 // this should not really happening, but if one sets content by hand
5021 // it can happen. a call to ResetStats() should be done in such cases
5022 double sumw = 0;
5023 for (int i = 0; i< GetNcells(); ++i) sumw += RetrieveBinContent(i);
5024 return (sumw != 0) ? kFALSE : kTRUE;
5025}
5026
5027////////////////////////////////////////////////////////////////////////////////
5028/// Return true if the bin is overflow.
5029
5030Bool_t TH1::IsBinOverflow(Int_t bin, Int_t iaxis) const
5031{
5032 Int_t binx, biny, binz;
5033 GetBinXYZ(bin, binx, biny, binz);
5034
5035 if (iaxis == 0) {
5036 if ( fDimension == 1 )
5037 return binx >= GetNbinsX() + 1;
5038 if ( fDimension == 2 )
5039 return (binx >= GetNbinsX() + 1) ||
5040 (biny >= GetNbinsY() + 1);
5041 if ( fDimension == 3 )
5042 return (binx >= GetNbinsX() + 1) ||
5043 (biny >= GetNbinsY() + 1) ||
5044 (binz >= GetNbinsZ() + 1);
5045 return kFALSE;
5046 }
5047 if (iaxis == 1)
5048 return binx >= GetNbinsX() + 1;
5049 if (iaxis == 2)
5050 return biny >= GetNbinsY() + 1;
5051 if (iaxis == 3)
5052 return binz >= GetNbinsZ() + 1;
5053
5054 Error("IsBinOverflow","Invalid axis value");
5055 return kFALSE;
5056}
5057
5058////////////////////////////////////////////////////////////////////////////////
5059/// Return true if the bin is underflow.
5060/// If iaxis = 0 make OR with all axes otherwise check only for the given axis
5061
5062Bool_t TH1::IsBinUnderflow(Int_t bin, Int_t iaxis) const
5063{
5064 Int_t binx, biny, binz;
5065 GetBinXYZ(bin, binx, biny, binz);
5066
5067 if (iaxis == 0) {
5068 if ( fDimension == 1 )
5069 return (binx <= 0);
5070 else if ( fDimension == 2 )
5071 return (binx <= 0 || biny <= 0);
5072 else if ( fDimension == 3 )
5073 return (binx <= 0 || biny <= 0 || binz <= 0);
5074 else
5075 return kFALSE;
5076 }
5077 if (iaxis == 1)
5078 return (binx <= 0);
5079 if (iaxis == 2)
5080 return (biny <= 0);
5081 if (iaxis == 3)
5082 return (binz <= 0);
5083
5084 Error("IsBinUnderflow","Invalid axis value");
5085 return kFALSE;
5086}
5087
5088////////////////////////////////////////////////////////////////////////////////
5089/// Reduce the number of bins for the axis passed in the option to the number of bins having a label.
5090/// The method will remove only the extra bins existing after the last "labeled" bin.
5091/// Note that if there are "un-labeled" bins present between "labeled" bins they will not be removed
5092
5094{
5095 Int_t iaxis = AxisChoice(ax);
5096 TAxis *axis = 0;
5097 if (iaxis == 1) axis = GetXaxis();
5098 if (iaxis == 2) axis = GetYaxis();
5099 if (iaxis == 3) axis = GetZaxis();
5100 if (!axis) {
5101 Error("LabelsDeflate","Invalid axis option %s",ax);
5102 return;
5103 }
5104 if (!axis->GetLabels()) return;
5105
5106 // find bin with last labels
5107 // bin number is object ID in list of labels
5108 // therefore max bin number is number of bins of the deflated histograms
5109 TIter next(axis->GetLabels());
5110 TObject *obj;
5111 Int_t nbins = 0;
5112 while ((obj = next())) {
5113 Int_t ibin = obj->GetUniqueID();
5114 if (ibin > nbins) nbins = ibin;
5115 }
5116 if (nbins < 1) nbins = 1;
5117
5118 // Do nothing in case it was the last bin
5119 if (nbins==axis->GetNbins()) return;
5120
5121 TH1 *hold = (TH1*)IsA()->New();
5122 R__ASSERT(hold);
5123 hold->SetDirectory(0);
5124 Copy(*hold);
5125
5126 Bool_t timedisp = axis->GetTimeDisplay();
5127 Double_t xmin = axis->GetXmin();
5128 Double_t xmax = axis->GetBinUpEdge(nbins);
5129 if (xmax <= xmin) xmax = xmin +nbins;
5130 axis->SetRange(0,0);
5131 axis->Set(nbins,xmin,xmax);
5132 SetBinsLength(-1); // reset the number of cells
5133 Int_t errors = fSumw2.fN;
5134 if (errors) fSumw2.Set(fNcells);
5135 axis->SetTimeDisplay(timedisp);
5136 // reset histogram content
5137 Reset("ICE");
5138
5139 //now loop on all bins and refill
5140 // NOTE that if the bins without labels have content
5141 // it will be put in the underflow/overflow.
5142 // For this reason we use AddBinContent method
5143 Double_t oldEntries = fEntries;
5144 Int_t bin,binx,biny,binz;
5145 for (bin=0; bin < hold->fNcells; ++bin) {
5146 hold->GetBinXYZ(bin,binx,biny,binz);
5147 Int_t ibin = GetBin(binx,biny,binz);
5148 Double_t cu = hold->RetrieveBinContent(bin);
5149 AddBinContent(ibin,cu);
5150 if (errors) {
5151 fSumw2.fArray[ibin] += hold->fSumw2.fArray[bin];
5152 }
5153 }
5154 fEntries = oldEntries;
5155 delete hold;
5156}
5157
5158////////////////////////////////////////////////////////////////////////////////
5159/// Double the number of bins for axis.
5160/// Refill histogram
5161/// This function is called by TAxis::FindBin(const char *label)
5162
5164{
5165 Int_t iaxis = AxisChoice(ax);
5166 TAxis *axis = 0;
5167 if (iaxis == 1) axis = GetXaxis();
5168 if (iaxis == 2) axis = GetYaxis();
5169 if (iaxis == 3) axis = GetZaxis();
5170 if (!axis) return;
5171
5172 TH1 *hold = (TH1*)IsA()->New();;
5173 hold->SetDirectory(0);
5174 Copy(*hold);
5175
5176 Bool_t timedisp = axis->GetTimeDisplay();
5177 Int_t nbins = axis->GetNbins();
5178 Double_t xmin = axis->GetXmin();
5179 Double_t xmax = axis->GetXmax();
5180 xmax = xmin + 2*(xmax-xmin);
5181 axis->SetRange(0,0);
5182 // double the bins and recompute ncells
5183 axis->Set(2*nbins,xmin,xmax);
5184 SetBinsLength(-1);
5185 Int_t errors = fSumw2.fN;
5186 if (errors) fSumw2.Set(fNcells);
5187 axis->SetTimeDisplay(timedisp);
5188
5189 Reset("ICE"); // reset content and error
5190
5191 //now loop on all bins and refill
5192 Double_t oldEntries = fEntries;
5193 Int_t bin,ibin,binx,biny,binz;
5194 for (ibin =0; ibin < hold->fNcells; ibin++) {
5195 // get the binx,y,z values . The x-y-z (axis) bin values will stay the same between new-old after the expanding
5196 hold->GetBinXYZ(ibin,binx,biny,binz);
5197 bin = GetBin(binx,biny,binz);
5198
5199 // underflow and overflow will be cleaned up because their meaning has been altered
5200 if (hold->IsBinUnderflow(ibin,iaxis) || hold->IsBinOverflow(ibin,iaxis)) {
5201 continue;
5202 }
5203 else {
5204 AddBinContent(bin, hold->RetrieveBinContent(ibin));
5205 if (errors) fSumw2.fArray[bin] += hold->fSumw2.fArray[ibin];
5206 }
5207 }
5208 fEntries = oldEntries;
5209 delete hold;
5210}
5211
5212////////////////////////////////////////////////////////////////////////////////
5213/// Set option(s) to draw axis with labels
5214/// \param[in] option
5215/// - "a" sort by alphabetic order
5216/// - ">" sort by decreasing values
5217/// - "<" sort by increasing values
5218/// - "h" draw labels horizontal
5219/// - "v" draw labels vertical
5220/// - "u" draw labels up (end of label right adjusted)
5221/// - "d" draw labels down (start of label left adjusted)
5222/// \param[in] ax axis
5223
5224void TH1::LabelsOption(Option_t *option, Option_t *ax)
5225{
5226 Int_t iaxis = AxisChoice(ax);
5227 TAxis *axis = 0;
5228 if (iaxis == 1) axis = GetXaxis();
5229 if (iaxis == 2) axis = GetYaxis();
5230 if (iaxis == 3) axis = GetZaxis();
5231 if (!axis) return;
5232 THashList *labels = axis->GetLabels();
5233 if (!labels) {
5234 Warning("LabelsOption","Cannot sort. No labels");
5235 return;
5236 }
5237 TString opt = option;
5238 opt.ToLower();
5239 if (opt.Contains("h")) {
5244 }
5245 if (opt.Contains("v")) {
5250 }
5251 if (opt.Contains("u")) {
5252 axis->SetBit(TAxis::kLabelsUp);
5256 }
5257 if (opt.Contains("d")) {
5262 }
5263 Int_t sort = -1;
5264 if (opt.Contains("a")) sort = 0;
5265 if (opt.Contains(">")) sort = 1;
5266 if (opt.Contains("<")) sort = 2;
5267 if (sort < 0) return;
5268 if (sort > 0 && GetDimension() > 2) {
5269 Error("LabelsOption","Sorting by value not implemented for 3-D histograms");