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