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