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