Logo ROOT   6.21/01
Reference Guide
THnBase.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Author: Axel Naumann (2011-12-20)
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2012, 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 "THnBase.h"
13 
14 #include "TAxis.h"
15 #include "TBrowser.h"
16 #include "TBuffer.h"
17 #include "TError.h"
18 #include "TClass.h"
19 #include "TF1.h"
20 #include "TH1D.h"
21 #include "TH2D.h"
22 #include "TH3D.h"
23 #include "THn.h"
24 #include "THnSparse.h"
25 #include "TMath.h"
26 #include "TRandom.h"
27 #include "TVirtualPad.h"
28 
29 #include "HFitInterface.h"
30 #include "Fit/DataRange.h"
31 #include "Fit/SparseData.h"
32 #include "Math/MinimizerOptions.h"
33 #include "Math/WrappedMultiTF1.h"
34 
35 
36 /** \class THnBase
37  \ingroup Hist
38 Multidimensional histogram base.
39 Defines common functionality and interfaces for THn, THnSparse.
40 */
41 
43 
44 ////////////////////////////////////////////////////////////////////////////////
45 /// Construct a THnBase with "dim" dimensions,
46 /// "nbins" holds the number of bins for each dimension;
47 /// "xmin" and "xmax" the minimal and maximal value for each dimension.
48 /// The arrays "xmin" and "xmax" can be NULL; in that case SetBinEdges()
49 /// must be called for each dimension.
50 
51 THnBase::THnBase(const char* name, const char* title, Int_t dim,
52  const Int_t* nbins, const Double_t* xmin, const Double_t* xmax):
53 TNamed(name, title), fNdimensions(dim), fAxes(dim), fBrowsables(dim),
54 fEntries(0), fTsumw(0), fTsumw2(-1.), fTsumwx(dim), fTsumwx2(dim),
55 fIntegral(0), fIntegralStatus(kNoInt)
56 {
57  for (Int_t i = 0; i < fNdimensions; ++i) {
58  TAxis* axis = new TAxis(nbins[i], xmin ? xmin[i] : 0., xmax ? xmax[i] : 1.);
59  axis->SetName(TString::Format("axis%d", i));
60  fAxes.AddAtAndExpand(axis, i);
61  }
62  SetTitle(title);
63  fAxes.SetOwner();
64 }
65 
66 ////////////////////////////////////////////////////////////////////////////////
67 /// Destruct a THnBase
68 
70  if (fIntegralStatus != kNoInt) delete [] fIntegral;
71 }
72 
73 
74 ////////////////////////////////////////////////////////////////////////////////
75 /// Create a new THnBase object that is of the same type as *this,
76 /// but with dimensions and bins given by axes.
77 /// If keepTargetAxis is true, the axes will keep their original xmin / xmax,
78 /// else they will be restricted to the range selected (first / last).
79 
80 THnBase* THnBase::CloneEmpty(const char* name, const char* title,
81  const TObjArray* axes, Bool_t keepTargetAxis) const
82 {
83  THnBase* ret = (THnBase*)IsA()->New();
84  Int_t chunkSize = 1024 * 16;
86  chunkSize = ((const THnSparse*)this)->GetChunkSize();
87  }
88  ret->Init(name, title, axes, keepTargetAxis, chunkSize);
89  return ret;
90 }
91 
92 
93 ////////////////////////////////////////////////////////////////////////////////
94 /// Initialize axes and name.
95 
96 void THnBase::Init(const char* name, const char* title,
97  const TObjArray* axes, Bool_t keepTargetAxis,
98  Int_t chunkSize /*= 1024 * 16*/)
99 {
100  SetNameTitle(name, title);
101 
102  TIter iAxis(axes);
103  const TAxis* axis = 0;
104  Int_t pos = 0;
105  Int_t *nbins = new Int_t[axes->GetEntriesFast()];
106  while ((axis = (TAxis*)iAxis())) {
107  TAxis* reqaxis = new TAxis(*axis);
108  if (!keepTargetAxis && axis->TestBit(TAxis::kAxisRange)) {
109  Int_t binFirst = axis->GetFirst();
110  // The lowest egde of the underflow is meaningless.
111  if (binFirst == 0)
112  binFirst = 1;
113  Int_t binLast = axis->GetLast();
114  // The overflow edge is implicit.
115  if (binLast > axis->GetNbins())
116  binLast = axis->GetNbins();
117  Int_t nBins = binLast - binFirst + 1;
118  if (axis->GetXbins()->GetSize()) {
119  // non-uniform bins:
120  reqaxis->Set(nBins, axis->GetXbins()->GetArray() + binFirst - 1);
121  } else {
122  // uniform bins:
123  reqaxis->Set(nBins, axis->GetBinLowEdge(binFirst), axis->GetBinUpEdge(binLast));
124  }
125  reqaxis->ResetBit(TAxis::kAxisRange);
126  }
127 
128  nbins[pos] = reqaxis->GetNbins();
129  fAxes.AddAtAndExpand(new TAxis(*reqaxis), pos++);
130  }
131  fAxes.SetOwner();
132 
133  fNdimensions = axes->GetEntriesFast();
134  InitStorage(nbins, chunkSize);
135  delete [] nbins;
136 }
137 
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 /// Create an empty histogram with name and title with a given
141 /// set of axes. Create a TH1D/TH2D/TH3D, depending on the number
142 /// of elements in axes.
143 
144 TH1* THnBase::CreateHist(const char* name, const char* title,
145  const TObjArray* axes,
146  Bool_t keepTargetAxis ) const {
147  const int ndim = axes->GetSize();
148 
149  TH1* hist = 0;
150  // create hist with dummy axes, we fix them later.
151  if (ndim == 1)
152  hist = new TH1D(name, title, 1, 0., 1.);
153  else if (ndim == 2)
154  hist = new TH2D(name, title, 1, 0., 1., 1, 0., 1.);
155  else if (ndim == 3)
156  hist = new TH3D(name, title, 1, 0., 1., 1, 0., 1., 1, 0., 1.);
157  else {
158  Error("CreateHist", "Cannot create histogram %s with %d dimensions!", name, ndim);
159  return 0;
160  }
161 
162  TAxis* hax[3] = {hist->GetXaxis(), hist->GetYaxis(), hist->GetZaxis()};
163  for (Int_t d = 0; d < ndim; ++d) {
164  TAxis* reqaxis = (TAxis*)(*axes)[d];
165  hax[d]->SetTitle(reqaxis->GetTitle());
166  if (!keepTargetAxis && reqaxis->TestBit(TAxis::kAxisRange)) {
167  // axis cannot extend to underflow/overflows (fix ROOT-8781)
168  Int_t binFirst = std::max(reqaxis->GetFirst(),1);
169  Int_t binLast = std::min(reqaxis->GetLast(), reqaxis->GetNbins() );
170  Int_t nBins = binLast - binFirst + 1;
171  if (reqaxis->GetXbins()->GetSize()) {
172  // non-uniform bins:
173  hax[d]->Set(nBins, reqaxis->GetXbins()->GetArray() + binFirst - 1);
174  } else {
175  // uniform bins:
176  hax[d]->Set(nBins, reqaxis->GetBinLowEdge(binFirst), reqaxis->GetBinUpEdge(binLast));
177  }
178  } else {
179  if (reqaxis->GetXbins()->GetSize()) {
180  // non-uniform bins:
181  hax[d]->Set(reqaxis->GetNbins(), reqaxis->GetXbins()->GetArray());
182  } else {
183  // uniform bins:
184  hax[d]->Set(reqaxis->GetNbins(), reqaxis->GetXmin(), reqaxis->GetXmax());
185  }
186  }
187  }
188 
189  hist->Rebuild();
190 
191  return hist;
192 }
193 
194 ////////////////////////////////////////////////////////////////////////////////
195 /// Create a THn / THnSparse object from a histogram deriving from TH1.
196 
197 THnBase* THnBase::CreateHnAny(const char* name, const char* title,
198  const TH1* h, Bool_t sparse, Int_t chunkSize)
199 {
200  // Get the dimension of the TH1
201  int ndim = h->GetDimension();
202 
203  // Axis properties
204  int nbins[3] = {0,0,0};
205  double minRange[3] = {0.,0.,0.};
206  double maxRange[3] = {0.,0.,0.};
207  const TAxis* axis[3] = { h->GetXaxis(), h->GetYaxis(), h->GetZaxis() };
208  for (int i = 0; i < ndim; ++i) {
209  nbins[i] = axis[i]->GetNbins();
210  minRange[i] = axis[i]->GetXmin();
211  maxRange[i] = axis[i]->GetXmax();
212  }
213 
214  // Create the corresponding THnSparse, depending on the storage
215  // type of the TH1. The class name will be "TH??\0" where the first
216  // ? is 1,2 or 3 and the second ? indicates the storage as C, S,
217  // I, F or D.
218  THnBase* s = 0;
219  const char* cname( h->ClassName() );
220  if (cname[0] == 'T' && cname[1] == 'H'
221  && cname[2] >= '1' && cname[2] <= '3' && cname[4] == 0) {
222 
223 #define R__THNBCASE(TAG) \
224 if (sparse) { \
225 s = new _NAME2_(THnSparse,TAG)(name, title, ndim, nbins, \
226 minRange, maxRange, chunkSize); \
227 } else { \
228 s = new _NAME2_(THn,TAG)(name, title, ndim, nbins, \
229 minRange, maxRange); \
230 } \
231 break;
232 
233  switch (cname[3]) {
234  case 'F': R__THNBCASE(F);
235  case 'D': R__THNBCASE(D);
236  case 'I': R__THNBCASE(I);
237  case 'S': R__THNBCASE(S);
238  case 'C': R__THNBCASE(C);
239  }
240 #undef R__THNBCASE
241  }
242  if (!s) {
243  ::Warning("THnSparse::CreateHnAny", "Unknown Type of Histogram");
244  return 0;
245  }
246 
247  for (int i = 0; i < ndim; ++i) {
248  s->GetAxis(i)->SetTitle(axis[i]->GetTitle());
249  }
250 
251  // Get the array to know the number of entries of the TH1
252  const TArray *array = dynamic_cast<const TArray*>(h);
253  if (!array) {
254  ::Warning("THnSparse::CreateHnAny", "Unknown Type of Histogram");
255  return 0;
256  }
257 
258  s->Add(h);
259  return s;
260 }
261 
262 
263 ////////////////////////////////////////////////////////////////////////////////
264 /// Create a THnSparse (if "sparse") or THn from "hn", possibly
265 /// converting THn <-> THnSparse.
266 
267 THnBase* THnBase::CreateHnAny(const char* name, const char* title,
268  const THnBase* hn, Bool_t sparse,
269  Int_t chunkSize /*= 1024 * 16*/)
270 {
271  TClass* type = 0;
272  if (hn->InheritsFrom(THnSparse::Class())) {
273  if (sparse) type = hn->IsA();
274  else {
275  char bintype;
276  if (hn->InheritsFrom(THnSparseD::Class())) bintype = 'D';
277  else if (hn->InheritsFrom(THnSparseF::Class())) bintype = 'F';
278  else if (hn->InheritsFrom(THnSparseL::Class())) bintype = 'L';
279  else if (hn->InheritsFrom(THnSparseI::Class())) bintype = 'I';
280  else if (hn->InheritsFrom(THnSparseS::Class())) bintype = 'S';
281  else if (hn->InheritsFrom(THnSparseC::Class())) bintype = 'C';
282  else {
283  hn->Error("CreateHnAny", "Type %s not implemented; please inform the ROOT team!",
284  hn->IsA()->GetName());
285  return 0;
286  }
287  type = TClass::GetClass(TString::Format("THn%c", bintype));
288  }
289  } else if (hn->InheritsFrom(THn::Class())) {
290  if (!sparse) type = hn->IsA();
291  else {
292  char bintype = 0;
293  if (hn->InheritsFrom(THnD::Class())) bintype = 'D';
294  else if (hn->InheritsFrom(THnF::Class())) bintype = 'F';
295  else if (hn->InheritsFrom(THnC::Class())) bintype = 'C';
296  else if (hn->InheritsFrom(THnS::Class())) bintype = 'S';
297  else if (hn->InheritsFrom(THnI::Class())) bintype = 'I';
298  else if (hn->InheritsFrom(THnL::Class())) bintype = 'L';
299  else if (hn->InheritsFrom(THnL64::Class())) {
300  hn->Error("CreateHnAny", "Type THnSparse with Long64_t bins is not available!");
301  return 0;
302  }
303  if (bintype) {
304  type = TClass::GetClass(TString::Format("THnSparse%c", bintype));
305  }
306  }
307  } else {
308  hn->Error("CreateHnAny", "Unhandled type %s, not deriving from THn nor THnSparse!",
309  hn->IsA()->GetName());
310  return 0;
311  }
312  if (!type) {
313  hn->Error("CreateHnAny", "Unhandled type %s, please inform the ROOT team!",
314  hn->IsA()->GetName());
315  return 0;
316  }
317 
318  THnBase* ret = (THnBase*)type->New();
319  ret->Init(name, title, hn->GetListOfAxes(),
320  kFALSE /*keepTargetAxes*/, chunkSize);
321 
322  ret->Add(hn);
323  return ret;
324 }
325 
326 
327 ////////////////////////////////////////////////////////////////////////////////
328 /// Fill the THnBase with the bins of hist that have content
329 /// or error != 0.
330 
331 void THnBase::Add(const TH1* hist, Double_t c /*=1.*/)
332 {
333  Long64_t nbins = hist->GetNcells();
334  int x[3] = {0,0,0};
335  for (int i = 0; i < nbins; ++i) {
336  double value = hist->GetBinContent(i);
337  double error = hist->GetBinError(i);
338  if (!value && !error) continue;
339  hist->GetBinXYZ(i, x[0], x[1], x[2]);
340  SetBinContent(x, value * c);
341  SetBinError(x, error * c);
342  }
343 }
344 
345 ////////////////////////////////////////////////////////////////////////////////
346 /// Fit a THnSparse with function f
347 ///
348 /// since the data is sparse by default a likelihood fit is performed
349 /// merging all the regions with empty bins for better performance efficiency
350 ///
351 /// Since the THnSparse is not drawn no graphics options are passed
352 /// Here is the list of possible options
353 ///
354 /// = "I" Use integral of function in bin instead of value at bin center
355 /// = "X" Use chi2 method (default is log-likelihood method)
356 /// = "U" Use a User specified fitting algorithm (via SetFCN)
357 /// = "Q" Quiet mode (minimum printing)
358 /// = "V" Verbose mode (default is between Q and V)
359 /// = "E" Perform better Errors estimation using Minos technique
360 /// = "B" Use this option when you want to fix one or more parameters
361 /// and the fitting function is like "gaus", "expo", "poln", "landau".
362 /// = "M" More. Improve fit results
363 /// = "R" Use the Range specified in the function range
364 
366 {
367 
368  Foption_t fitOption;
369 
370  if (!TH1::FitOptionsMake(option,fitOption)) return 0;
371 
372  // The function used to fit cannot be stored in a THnSparse. It
373  // cannot be drawn either. Perhaps in the future.
374  fitOption.Nostore = true;
375  // Use likelihood fit if not specified
376  if (!fitOption.Chi2) fitOption.Like = true;
377  // create range and minimizer options with default values
379  for ( int i = 0; i < GetNdimensions(); ++i ) {
380  TAxis *axis = GetAxis(i);
381  range.AddRange(i, axis->GetXmin(), axis->GetXmax());
382  }
384 
385  return ROOT::Fit::FitObject(this, f , fitOption , minOption, goption, range);
386 }
387 
388 ////////////////////////////////////////////////////////////////////////////////
389 /// Generate an n-dimensional random tuple based on the histogrammed
390 /// distribution. If subBinRandom, the returned tuple will be additionally
391 /// randomly distributed within the randomized bin, using a flat
392 /// distribution.
393 
394 void THnBase::GetRandom(Double_t *rand, Bool_t subBinRandom /* = kTRUE */)
395 {
396  // check whether the integral array is valid
397  if (fIntegralStatus != kValidInt)
398  ComputeIntegral();
399 
400  // generate a random bin
401  Double_t p = gRandom->Rndm();
403  const Int_t nStaticBins = 40;
404  Int_t bin[nStaticBins];
405  Int_t* pBin = bin;
406  if (GetNdimensions() > nStaticBins) {
407  pBin = new Int_t[GetNdimensions()];
408  }
409  GetBinContent(idx, pBin);
410 
411  // convert bin coordinates to real values
412  for (Int_t i = 0; i < fNdimensions; i++) {
413  rand[i] = GetAxis(i)->GetBinCenter(pBin[i]);
414 
415  // randomize the vector within a bin
416  if (subBinRandom)
417  rand[i] += (gRandom->Rndm() - 0.5) * GetAxis(i)->GetBinWidth(pBin[i]);
418  }
419  if (pBin != bin) {
420  delete [] pBin;
421  }
422 
423  return;
424 }
425 
426 ////////////////////////////////////////////////////////////////////////////////
427 /// Check whether bin coord is in range, as defined by TAxis::SetRange().
428 
430 {
431  Int_t min = 0;
432  Int_t max = 0;
433  for (Int_t i = 0; i < fNdimensions; ++i) {
434  TAxis *axis = GetAxis(i);
435  if (!axis->TestBit(TAxis::kAxisRange)) continue;
436  min = axis->GetFirst();
437  max = axis->GetLast();
438  if (coord[i] < min || coord[i] > max)
439  return kFALSE;
440  }
441  return kTRUE;
442 }
443 
444 ////////////////////////////////////////////////////////////////////////////////
445 /// Project all bins into a ndim-dimensional THn / THnSparse (whatever
446 /// *this is) or if (ndim < 4 and !wantNDim) a TH1/2/3 histogram,
447 /// keeping only axes in dim (specifying ndim dimensions).
448 /// If "option" contains:
449 /// - "E" errors will be calculated.
450 /// - "A" ranges of the target axes will be ignored.
451 /// - "O" original axis range of the target axes will be
452 /// kept, but only bins inside the selected range
453 /// will be filled.
454 
456  Bool_t wantNDim,
457  Option_t* option /*= ""*/) const
458 {
459  TString name(GetName());
460  name +="_proj";
461 
462  for (Int_t d = 0; d < ndim; ++d) {
463  name += "_";
464  name += dim[d];
465  }
466 
467  TString title(GetTitle());
468  Ssiz_t posInsert = title.First(';');
469  if (posInsert == kNPOS) {
470  title += " projection ";
471  for (Int_t d = 0; d < ndim; ++d)
472  title += GetAxis(dim[d])->GetTitle();
473  } else {
474  for (Int_t d = ndim - 1; d >= 0; --d) {
475  title.Insert(posInsert, GetAxis(d)->GetTitle());
476  if (d)
477  title.Insert(posInsert, ", ");
478  }
479  title.Insert(posInsert, " projection ");
480  }
481 
482  TObjArray newaxes(ndim);
483  for (Int_t d = 0; d < ndim; ++d) {
484  newaxes.AddAt(GetAxis(dim[d]),d);
485  }
486 
487  THnBase* hn = 0;
488  TH1* hist = 0;
489  TObject* ret = 0;
490 
491  Bool_t* hadRange = 0;
492  Bool_t ignoreTargetRange = (option && (strchr(option, 'A') || strchr(option, 'a')));
493  Bool_t keepTargetAxis = ignoreTargetRange || (option && (strchr(option, 'O') || strchr(option, 'o')));
494  if (ignoreTargetRange) {
495  hadRange = new Bool_t[ndim];
496  for (Int_t d = 0; d < ndim; ++d){
497  TAxis *axis = GetAxis(dim[d]);
498  hadRange[d] = axis->TestBit(TAxis::kAxisRange);
500  }
501  }
502 
503  if (wantNDim)
504  ret = hn = CloneEmpty(name, title, &newaxes, keepTargetAxis);
505  else
506  ret = hist = CreateHist(name, title, &newaxes, keepTargetAxis);
507 
508  if (keepTargetAxis) {
509  // make the whole axes visible, i.e. unset the range
510  if (wantNDim) {
511  for (Int_t d = 0; d < ndim; ++d) {
512  hn->GetAxis(d)->SetRange(0, 0);
513  }
514  } else {
515  hist->GetXaxis()->SetRange(0, 0);
516  hist->GetYaxis()->SetRange(0, 0);
517  hist->GetZaxis()->SetRange(0, 0);
518  }
519  }
520 
521  Bool_t haveErrors = GetCalculateErrors();
522  Bool_t wantErrors = haveErrors || (option && (strchr(option, 'E') || strchr(option, 'e')));
523 
524  Int_t* bins = new Int_t[ndim];
525  Long64_t myLinBin = 0;
526 
527  THnIter iter(this, kTRUE /*use axis range*/);
528 
529  while ((myLinBin = iter.Next()) >= 0) {
530  Double_t v = GetBinContent(myLinBin);
531 
532  for (Int_t d = 0; d < ndim; ++d) {
533  bins[d] = iter.GetCoord(dim[d]);
534  if (!keepTargetAxis && GetAxis(dim[d])->TestBit(TAxis::kAxisRange)) {
535  Int_t binOffset = GetAxis(dim[d])->GetFirst();
536  // Don't subtract even more if underflow is alreday included:
537  if (binOffset > 0) --binOffset;
538  bins[d] -= binOffset;
539  }
540  }
541 
542  Long64_t targetLinBin = -1;
543  if (!wantNDim) {
544  if (ndim == 1) targetLinBin = bins[0];
545  else if (ndim == 2) targetLinBin = hist->GetBin(bins[0], bins[1]);
546  else if (ndim == 3) targetLinBin = hist->GetBin(bins[0], bins[1], bins[2]);
547  } else {
548  targetLinBin = hn->GetBin(bins, kTRUE /*allocate*/);
549  }
550 
551  if (wantErrors) {
552  Double_t err2 = 0.;
553  if (haveErrors) {
554  err2 = GetBinError2(myLinBin);
555  } else {
556  err2 = v;
557  }
558  if (wantNDim) {
559  hn->AddBinError2(targetLinBin, err2);
560  } else {
561  Double_t preverr = hist->GetBinError(targetLinBin);
562  hist->SetBinError(targetLinBin, TMath::Sqrt(preverr * preverr + err2));
563  }
564  }
565 
566  // only _after_ error calculation, or sqrt(v) is taken into account!
567  if (wantNDim)
568  hn->AddBinContent(targetLinBin, v);
569  else
570  hist->AddBinContent(targetLinBin, v);
571  }
572 
573  delete [] bins;
574 
575  if (wantNDim) {
576  hn->SetEntries(fEntries);
577  } else {
578  if (!iter.HaveSkippedBin()) {
579  hist->SetEntries(fEntries);
580  } else {
581  // re-compute the entries
582  // in case of error calculation (i.e. when Sumw2() is set)
583  // use the effective entries for the entries
584  // since this is the only way to estimate them
585  hist->ResetStats();
586  Double_t entries = hist->GetEffectiveEntries();
587  if (!wantErrors) {
588  // to avoid numerical rounding
589  entries = TMath::Floor(entries + 0.5);
590  }
591  hist->SetEntries(entries);
592  }
593  }
594 
595  if (hadRange) {
596  // reset kAxisRange bit:
597  for (Int_t d = 0; d < ndim; ++d)
598  GetAxis(dim[d])->SetBit(TAxis::kAxisRange, hadRange[d]);
599 
600  delete [] hadRange;
601  }
602 
603  return ret;
604 }
605 
606 ////////////////////////////////////////////////////////////////////////////////
607 /// Scale contents and errors of this histogram by c:
608 /// this = this * c
609 /// It does not modify the histogram's number of entries.
610 
612 {
613 
614  Double_t nEntries = GetEntries();
615  // Scale the contents & errors
616  Bool_t haveErrors = GetCalculateErrors();
617  Long64_t i = 0;
618  THnIter iter(this);
619  while ((i = iter.Next()) >= 0) {
620  // Get the content of the bin from the current histogram
621  Double_t v = GetBinContent(i);
622  SetBinContent(i, c * v);
623  if (haveErrors) {
624  Double_t err2 = GetBinError2(i);
625  SetBinError2(i, c * c * err2);
626  }
627  }
628  SetEntries(nEntries);
629 }
630 
631 ////////////////////////////////////////////////////////////////////////////////
632 /// Add() implementation for both rebinned histograms and those with identical
633 /// binning. See THnBase::Add().
634 
636 {
637  if (fNdimensions != h->GetNdimensions()) {
638  Warning("RebinnedAdd", "Different number of dimensions, cannot carry out operation on the histograms");
639  return;
640  }
641 
642  // Trigger error calculation if h has it
643  if (!GetCalculateErrors() && h->GetCalculateErrors())
644  Sumw2();
645  Bool_t haveErrors = GetCalculateErrors();
646 
647  Double_t* x = 0;
648  if (rebinned) {
649  x = new Double_t[fNdimensions];
650  }
651  Int_t* coord = new Int_t[fNdimensions];
652 
653  // Expand the exmap if needed, to reduce collisions
654  Long64_t numTargetBins = GetNbins() + h->GetNbins();
655  Reserve(numTargetBins);
656 
657  Long64_t i = 0;
658  THnIter iter(h);
659  // Add to this whatever is found inside the other histogram
660  while ((i = iter.Next(coord)) >= 0) {
661  // Get the content of the bin from the second histogram
662  Double_t v = h->GetBinContent(i);
663 
664  Long64_t mybinidx = -1;
665  if (rebinned) {
666  // Get the bin center given a coord
667  for (Int_t j = 0; j < fNdimensions; ++j)
668  x[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
669 
670  mybinidx = GetBin(x, kTRUE /* allocate*/);
671  } else {
672  mybinidx = GetBin(coord, kTRUE /*allocate*/);
673  }
674 
675  if (haveErrors) {
676  Double_t err2 = h->GetBinError2(i) * c * c;
677  AddBinError2(mybinidx, err2);
678  }
679  // only _after_ error calculation, or sqrt(v) is taken into account!
680  AddBinContent(mybinidx, c * v);
681  }
682 
683  delete [] coord;
684  delete [] x;
685 
686  Double_t nEntries = GetEntries() + c * h->GetEntries();
687  SetEntries(nEntries);
688 }
689 
690 ////////////////////////////////////////////////////////////////////////////////
691 /// Add contents of h scaled by c to this histogram:
692 /// this = this + c * h
693 /// Note that if h has Sumw2 set, Sumw2 is automatically called for this
694 /// if not already set.
695 
697 {
698  // Check consistency of the input
699  if (!CheckConsistency(h, "Add")) return;
700 
701  AddInternal(h, c, kFALSE);
702 }
703 
704 ////////////////////////////////////////////////////////////////////////////////
705 /// Add contents of h scaled by c to this histogram:
706 /// this = this + c * h
707 /// Note that if h has Sumw2 set, Sumw2 is automatically called for this
708 /// if not already set.
709 /// In contrast to Add(), RebinnedAdd() does not require consistent binning of
710 /// this and h; instead, each bin's center is used to determine the target bin.
711 
713 {
714  AddInternal(h, c, kTRUE);
715 }
716 
717 
718 ////////////////////////////////////////////////////////////////////////////////
719 /// Merge this with a list of THnBase's. All THnBase's provided
720 /// in the list must have the same bin layout!
721 
723 {
724  if (!list) return 0;
725  if (list->IsEmpty()) return (Long64_t)GetEntries();
726 
727  Long64_t sumNbins = GetNbins();
728  TIter iter(list);
729  const TObject* addMeObj = 0;
730  while ((addMeObj = iter())) {
731  const THnBase* addMe = dynamic_cast<const THnBase*>(addMeObj);
732  if (addMe) {
733  sumNbins += addMe->GetNbins();
734  }
735  }
736  Reserve(sumNbins);
737 
738  iter.Reset();
739  while ((addMeObj = iter())) {
740  const THnBase* addMe = dynamic_cast<const THnBase*>(addMeObj);
741  if (!addMe)
742  Error("Merge", "Object named %s is not THnBase! Skipping it.",
743  addMeObj->GetName());
744  else
745  Add(addMe);
746  }
747  return (Long64_t)GetEntries();
748 }
749 
750 
751 ////////////////////////////////////////////////////////////////////////////////
752 /// Multiply this histogram by histogram h
753 /// this = this * h
754 /// Note that if h has Sumw2 set, Sumw2 is automatically called for this
755 /// if not already set.
756 
758 {
759  // Check consistency of the input
760  if(!CheckConsistency(h, "Multiply"))return;
761 
762  // Trigger error calculation if h has it
763  Bool_t wantErrors = kFALSE;
764  if (GetCalculateErrors() || h->GetCalculateErrors())
765  wantErrors = kTRUE;
766 
767  if (wantErrors) Sumw2();
768 
769  Double_t nEntries = GetEntries();
770  // Now multiply the contents: in this case we have the intersection of the sets of bins
771  Int_t* coord = new Int_t[fNdimensions];
772  Long64_t i = 0;
773  THnIter iter(this);
774  // Add to this whatever is found inside the other histogram
775  while ((i = iter.Next(coord)) >= 0) {
776  // Get the content of the bin from the current histogram
778  // Now look at the bin with the same coordinates in h
779  Long64_t idxh = h->GetBin(coord);
780  Double_t v2 = 0.;
781  if (idxh >= 0) v2 = h->GetBinContent(idxh);
782  SetBinContent(i, v1 * v2);
783  if (wantErrors) {
784  Double_t err1 = GetBinError(i) * v2;
785  Double_t err2 = 0.;
786  if (idxh >= 0) err2 = h->GetBinError(idxh) * v1;
787  SetBinError(i, TMath::Sqrt((err2 * err2 + err1 * err1)));
788  }
789  }
790  SetEntries(nEntries);
791 
792  delete [] coord;
793 }
794 
795 ////////////////////////////////////////////////////////////////////////////////
796 /// Performs the operation: this = this*c*f1
797 /// if errors are defined, errors are also recalculated.
798 ///
799 /// Only bins inside the function range are recomputed.
800 /// IMPORTANT NOTE: If you intend to use the errors of this histogram later
801 /// you should call Sumw2 before making this operation.
802 /// This is particularly important if you fit the histogram after
803 /// calling Multiply()
804 
806 {
807  Int_t* coord = new Int_t[fNdimensions];
809 
810  Bool_t wantErrors = GetCalculateErrors();
811  if (wantErrors) Sumw2();
812 
813  Long64_t i = 0;
814  THnIter iter(this);
815  // Add to this whatever is found inside the other histogram
816  while ((i = iter.Next(coord)) >= 0) {
817  Double_t value = GetBinContent(i);
818 
819  // Get the bin coordinates given an index array
820  for (Int_t j = 0; j < fNdimensions; ++j)
821  x[j] = GetAxis(j)->GetBinCenter(coord[j]);
822 
823  if (!f->IsInside(x))
824  continue;
826 
827  // Evaluate function at points
828  Double_t fvalue = f->EvalPar(x, NULL);
829 
830  SetBinContent(i, c * fvalue * value);
831  if (wantErrors) {
832  Double_t error = GetBinError(i);
833  SetBinError(i, c * fvalue * error);
834  }
835  }
836 
837  delete [] x;
838  delete [] coord;
839 }
840 
841 ////////////////////////////////////////////////////////////////////////////////
842 /// Divide this histogram by h
843 /// this = this/(h)
844 /// Note that if h has Sumw2 set, Sumw2 is automatically called for
845 /// this if not already set.
846 /// The resulting errors are calculated assuming uncorrelated content.
847 
849 {
850  // Check consistency of the input
851  if (!CheckConsistency(h, "Divide"))return;
852 
853  // Trigger error calculation if h has it
854  Bool_t wantErrors=GetCalculateErrors();
855  if (!GetCalculateErrors() && h->GetCalculateErrors())
856  wantErrors=kTRUE;
857 
858  // Remember original histogram statistics
859  Double_t nEntries = fEntries;
860 
861  if (wantErrors) Sumw2();
862  Bool_t didWarn = kFALSE;
863 
864  // Now divide the contents: also in this case we have the intersection of the sets of bins
865  Int_t* coord = new Int_t[fNdimensions];
866  Long64_t i = 0;
867  THnIter iter(this);
868  // Add to this whatever is found inside the other histogram
869  while ((i = iter.Next(coord)) >= 0) {
870  // Get the content of the bin from the first histogram
872  // Now look at the bin with the same coordinates in h
873  Long64_t hbin = h->GetBin(coord);
874  Double_t v2 = h->GetBinContent(hbin);
875  if (!v2) {
876  v1 = 0.;
877  v2 = 1.;
878  if (!didWarn) {
879  Warning("Divide(h)", "Histogram h has empty bins - division by zero! Setting bin to 0.");
880  didWarn = kTRUE;
881  }
882  }
883  SetBinContent(i, v1 / v2);
884  if (wantErrors) {
885  Double_t err1 = GetBinError(i) * v2;
886  Double_t err2 = h->GetBinError(hbin) * v1;
887  Double_t b22 = v2 * v2;
888  Double_t err = (err1 * err1 + err2 * err2) / (b22 * b22);
889  SetBinError2(i, err);
890  }
891  }
892  delete [] coord;
893  SetEntries(nEntries);
894 }
895 
896 ////////////////////////////////////////////////////////////////////////////////
897 /// Replace contents of this histogram by multiplication of h1 by h2
898 /// this = (c1*h1)/(c2*h2)
899 /// Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for
900 /// this if not already set.
901 /// The resulting errors are calculated assuming uncorrelated content.
902 /// However, if option ="B" is specified, Binomial errors are computed.
903 /// In this case c1 and c2 do not make real sense and they are ignored.
904 
905 void THnBase::Divide(const THnBase *h1, const THnBase *h2, Double_t c1, Double_t c2, Option_t *option)
906 {
907 
908  TString opt = option;
909  opt.ToLower();
910  Bool_t binomial = kFALSE;
911  if (opt.Contains("b")) binomial = kTRUE;
912 
913  // Check consistency of the input
914  if (!CheckConsistency(h1, "Divide") || !CheckConsistency(h2, "Divide"))return;
915  if (!c2) {
916  Error("Divide","Coefficient of dividing histogram cannot be zero");
917  return;
918  }
919 
920  Reset();
921 
922  // Trigger error calculation if h1 or h2 have it
923  if (!GetCalculateErrors() && (h1->GetCalculateErrors()|| h2->GetCalculateErrors() != 0))
924  Sumw2();
925 
926  // Count filled bins
927  Long64_t nFilledBins=0;
928 
929  // Now divide the contents: we have the intersection of the sets of bins
930 
931  Int_t* coord = new Int_t[fNdimensions];
932  memset(coord, 0, sizeof(Int_t) * fNdimensions);
933  Bool_t didWarn = kFALSE;
934 
935  Long64_t i = 0;
936  THnIter iter(h1);
937  // Add to this whatever is found inside the other histogram
938  while ((i = iter.Next(coord)) >= 0) {
939  // Get the content of the bin from the first histogram
940  Double_t v1 = h1->GetBinContent(i);
941  // Now look at the bin with the same coordinates in h2
942  Long64_t h2bin = h2->GetBin(coord);
943  Double_t v2 = h2->GetBinContent(h2bin);
944  if (!v2) {
945  v1 = 0.;
946  v2 = 1.;
947  if (!didWarn) {
948  Warning("Divide(h1, h2)", "Histogram h2 has empty bins - division by zero! Setting bin to 0.");
949  didWarn = kTRUE;
950  }
951  }
952  nFilledBins++;
953  Long64_t myBin = GetBin(coord);
954  SetBinContent(myBin, c1 * v1 / c2 / v2);
955  if(GetCalculateErrors()){
956  Double_t err1 = h1->GetBinError(i);
957  Double_t err2 = h2->GetBinError(h2bin);
958  Double_t errSq = 0.;
959  if (binomial) {
960  if (v1 != v2) {
961  Double_t w = v1 / v2;
962  err2 *= w;
963  errSq = TMath::Abs( ( (1. - 2.*w) * err1 * err1 + err2 * err2 ) / (v2 * v2) );
964  }
965  } else {
966  c1 *= c1;
967  c2 *= c2;
968  Double_t b22 = v2 * v2 * c2;
969  err1 *= v2;
970  err2 *= v1;
971  errSq = c1 * c2 * (err1 * err1 + err2 * err2) / (b22 * b22);
972  }
973  SetBinError2(myBin, errSq);
974  }
975  }
976 
977  delete [] coord;
978  SetFilledBins(nFilledBins);
979 
980  // Set as entries in the result histogram the entries in the numerator
982 }
983 
984 ////////////////////////////////////////////////////////////////////////////////
985 /// Consistency check on (some of) the parameters of two histograms (for operations).
986 
987 Bool_t THnBase::CheckConsistency(const THnBase *h, const char *tag) const
988 {
989  if (fNdimensions != h->GetNdimensions()) {
990  Warning(tag, "Different number of dimensions, cannot carry out operation on the histograms");
991  return kFALSE;
992  }
993  for (Int_t dim = 0; dim < fNdimensions; dim++){
994  if (GetAxis(dim)->GetNbins() != h->GetAxis(dim)->GetNbins()) {
995  Warning(tag, "Different number of bins on axis %i, cannot carry out operation on the histograms", dim);
996  return kFALSE;
997  }
998  }
999  return kTRUE;
1000 }
1001 
1002 ////////////////////////////////////////////////////////////////////////////////
1003 /// Set the axis # of bins and bin limits on dimension idim
1004 
1005 void THnBase::SetBinEdges(Int_t idim, const Double_t* bins)
1006 {
1007  TAxis* axis = (TAxis*) fAxes[idim];
1008  axis->Set(axis->GetNbins(), bins);
1009 }
1010 
1011 ////////////////////////////////////////////////////////////////////////////////
1012 /// Change (i.e. set) the title.
1013 ///
1014 /// If title is in the form "stringt;string0;string1;string2 ..."
1015 /// the histogram title is set to stringt, the title of axis0 to string0,
1016 /// of axis1 to string1, of axis2 to string2, etc, just like it is done
1017 /// for TH1/TH2/TH3.
1018 /// To insert the character ";" in one of the titles, one should use "#;"
1019 /// or "#semicolon".
1020 
1021 void THnBase::SetTitle(const char *title)
1022 {
1023  fTitle = title;
1024  fTitle.ReplaceAll("#;",2,"#semicolon",10);
1025 
1026  Int_t endHistTitle = fTitle.First(';');
1027  if (endHistTitle >= 0) {
1028  // title contains a ';' so parse the axis titles
1029  Int_t posTitle = endHistTitle + 1;
1030  Int_t lenTitle = fTitle.Length();
1031  Int_t dim = 0;
1032  while (posTitle > 0 && posTitle < lenTitle && dim < fNdimensions){
1033  Int_t endTitle = fTitle.Index(";", posTitle);
1034  TString axisTitle = fTitle(posTitle, endTitle - posTitle);
1035  axisTitle.ReplaceAll("#semicolon", 10, ";", 1);
1036  GetAxis(dim)->SetTitle(axisTitle);
1037  dim++;
1038  if (endTitle > 0)
1039  posTitle = endTitle + 1;
1040  else
1041  posTitle = -1;
1042  }
1043  // Remove axis titles from histogram title
1044  fTitle.Remove(endHistTitle, lenTitle - endHistTitle);
1045  }
1046 
1047  fTitle.ReplaceAll("#semicolon", 10, ";", 1);
1048 
1049 }
1050 
1051 ////////////////////////////////////////////////////////////////////////////////
1052 /// Combine the content of "group" neighboring bins into
1053 /// a new bin and return the resulting THnBase.
1054 /// For group=2 and a 3 dimensional histogram, all "blocks"
1055 /// of 2*2*2 bins will be put into a bin.
1056 
1058 {
1059  Int_t* ngroup = new Int_t[GetNdimensions()];
1060  for (Int_t d = 0; d < GetNdimensions(); ++d)
1061  ngroup[d] = group;
1062  THnBase* ret = RebinBase(ngroup);
1063  delete [] ngroup;
1064  return ret;
1065 }
1066 
1067 ////////////////////////////////////////////////////////////////////////////////
1068 /// Combine the content of "group" neighboring bins for each dimension
1069 /// into a new bin and return the resulting THnBase.
1070 /// For group={2,1,1} and a 3 dimensional histogram, pairs of x-bins
1071 /// will be grouped.
1072 
1074 {
1075  Int_t ndim = GetNdimensions();
1076  TString name(GetName());
1077  for (Int_t d = 0; d < ndim; ++d)
1078  name += Form("_%d", group[d]);
1079 
1080 
1081  TString title(GetTitle());
1082  Ssiz_t posInsert = title.First(';');
1083  if (posInsert == kNPOS) {
1084  title += " rebin ";
1085  for (Int_t d = 0; d < ndim; ++d)
1086  title += Form("{%d}", group[d]);
1087  } else {
1088  for (Int_t d = ndim - 1; d >= 0; --d)
1089  title.Insert(posInsert, Form("{%d}", group[d]));
1090  title.Insert(posInsert, " rebin ");
1091  }
1092 
1093  TObjArray newaxes(ndim);
1094  newaxes.SetOwner();
1095  for (Int_t d = 0; d < ndim; ++d) {
1096  newaxes.AddAt(new TAxis(*GetAxis(d) ),d);
1097  if (group[d] > 1) {
1098  TAxis* newaxis = (TAxis*) newaxes.At(d);
1099  Int_t newbins = (newaxis->GetNbins() + group[d] - 1) / group[d];
1100  if (newaxis->GetXbins() && newaxis->GetXbins()->GetSize()) {
1101  // variable bins
1102  Double_t *edges = new Double_t[newbins + 1];
1103  for (Int_t i = 0; i < newbins + 1; ++i)
1104  if (group[d] * i <= newaxis->GetNbins())
1105  edges[i] = newaxis->GetXbins()->At(group[d] * i);
1106  else edges[i] = newaxis->GetXmax();
1107  newaxis->Set(newbins, edges);
1108  delete [] edges;
1109  } else {
1110  newaxis->Set(newbins, newaxis->GetXmin(), newaxis->GetXmax());
1111  }
1112  }
1113  }
1114 
1115  THnBase* h = CloneEmpty(name.Data(), title.Data(), &newaxes, kTRUE);
1116  Bool_t haveErrors = GetCalculateErrors();
1117  Bool_t wantErrors = haveErrors;
1118 
1119  Int_t* bins = new Int_t[ndim];
1120  Int_t* coord = new Int_t[fNdimensions];
1121 
1122  Long64_t i = 0;
1123  THnIter iter(this);
1124  while ((i = iter.Next(coord)) >= 0) {
1125  Double_t v = GetBinContent(i);
1126  for (Int_t d = 0; d < ndim; ++d) {
1127  bins[d] = TMath::CeilNint( (double) coord[d]/group[d] );
1128  }
1129  Long64_t idxh = h->GetBin(bins, kTRUE /*allocate*/);
1130 
1131  if (wantErrors) {
1132  // wantErrors == haveErrors, thus:
1133  h->AddBinError2(idxh, GetBinError2(i));
1134  }
1135 
1136  // only _after_ error calculation, or sqrt(v) is taken into account!
1137  h->AddBinContent(idxh, v);
1138  }
1139 
1140  delete [] bins;
1141  delete [] coord;
1142 
1143  h->SetEntries(fEntries);
1144 
1145  return h;
1146 
1147 }
1148 
1149 ////////////////////////////////////////////////////////////////////////////////
1150 /// Clear the histogram
1151 
1152 void THnBase::ResetBase(Option_t * /*option = ""*/)
1153 {
1154  fEntries = 0.;
1155  fTsumw = 0.;
1156  fTsumw2 = -1.;
1157  if (fIntegralStatus != kNoInt) {
1158  delete [] fIntegral;
1160  }
1161 }
1162 
1163 ////////////////////////////////////////////////////////////////////////////////
1164 /// Calculate the integral of the histogram
1165 
1167 {
1168  // delete old integral
1169  if (fIntegralStatus != kNoInt) {
1170  delete [] fIntegral;
1172  }
1173 
1174  // check number of bins
1175  if (GetNbins() == 0) {
1176  Error("ComputeIntegral", "The histogram must have at least one bin.");
1177  return 0.;
1178  }
1179 
1180  // allocate integral array
1181  fIntegral = new Double_t [GetNbins() + 1];
1182  fIntegral[0] = 0.;
1183 
1184  // fill integral array with contents of regular bins (non over/underflow)
1185  Int_t* coord = new Int_t[fNdimensions];
1186  Long64_t i = 0;
1187  THnIter iter(this);
1188  while ((i = iter.Next(coord)) >= 0) {
1189  Double_t v = GetBinContent(i);
1190 
1191  // check whether the bin is regular
1192  bool regularBin = true;
1193  for (Int_t dim = 0; dim < fNdimensions; dim++) {
1194  if (coord[dim] < 1 || coord[dim] > GetAxis(dim)->GetNbins()) {
1195  regularBin = false;
1196  break;
1197  }
1198  }
1199 
1200  // if outlayer, count it with zero weight
1201  if (!regularBin) v = 0.;
1202 
1203  fIntegral[i + 1] = fIntegral[i] + v;
1204  }
1205  delete [] coord;
1206 
1207  // check sum of weights
1208  if (fIntegral[GetNbins()] == 0.) {
1209  Error("ComputeIntegral", "No hits in regular bins (non over/underflow).");
1210  delete [] fIntegral;
1211  return 0.;
1212  }
1213 
1214  // normalize the integral array
1215  for (Long64_t j = 0; j <= GetNbins(); ++j)
1216  fIntegral[j] = fIntegral[j] / fIntegral[GetNbins()];
1217 
1218  // set status to valid
1220  return fIntegral[GetNbins()];
1221 }
1222 
1223 ////////////////////////////////////////////////////////////////////////////////
1224 /// Print bin with linex index "idx".
1225 /// For valid options see PrintBin(Long64_t idx, Int_t* bin, Option_t* options).
1226 
1227 void THnBase::PrintBin(Long64_t idx, Option_t* options) const
1228 {
1229  Int_t* coord = new Int_t[fNdimensions];
1230  PrintBin(idx, coord, options);
1231  delete [] coord;
1232 }
1233 
1234 ////////////////////////////////////////////////////////////////////////////////
1235 /// Print one bin. If "idx" is != -1 use that to determine the bin,
1236 /// otherwise (if "idx" == -1) use the coordinate in "bin".
1237 /// If "options" contains:
1238 /// - '0': only print bins with an error or content != 0
1239 /// Return whether the bin was printed (depends on options)
1240 
1241 Bool_t THnBase::PrintBin(Long64_t idx, Int_t* bin, Option_t* options) const
1242 {
1243  Double_t v = -42;
1244  if (idx == -1) {
1245  idx = GetBin(bin);
1246  v = GetBinContent(idx);
1247  } else {
1248  v = GetBinContent(idx, bin);
1249  }
1250 
1251  Double_t err = 0.;
1252  if (GetCalculateErrors()) {
1253  if (idx != -1) {
1254  err = GetBinError(idx);
1255  }
1256  }
1257 
1258  if (v == 0. && err == 0. && options && strchr(options, '0')) {
1259  // suppress zeros, and we have one.
1260  return kFALSE;
1261  }
1262 
1263  TString coord;
1264  for (Int_t dim = 0; dim < fNdimensions; ++dim) {
1265  coord += bin[dim];
1266  coord += ',';
1267  }
1268  coord.Remove(coord.Length() - 1);
1269 
1270  if (GetCalculateErrors()) {
1271  Printf("Bin at (%s) = %g (+/- %g)", coord.Data(), v, err);
1272  } else {
1273  Printf("Bin at (%s) = %g", coord.Data(), v);
1274  }
1275 
1276  return kTRUE;
1277 }
1278 
1279 ////////////////////////////////////////////////////////////////////////////////
1280 /// Print "howmany" entries starting at "from". If "howmany" is -1, print all.
1281 /// If "options" contains:
1282 /// - 'x': print in the order of axis bins, i.e. (0,0,...,0), (0,0,...,1),...
1283 /// - '0': only print bins with content != 0
1284 
1285 void THnBase::PrintEntries(Long64_t from /*=0*/, Long64_t howmany /*=-1*/,
1286  Option_t* options /*=0*/) const
1287 {
1288  if (from < 0) from = 0;
1289  if (howmany == -1) howmany = GetNbins();
1290 
1291  Int_t* bin = new Int_t[fNdimensions];
1292 
1293  if (options && (strchr(options, 'x') || strchr(options, 'X'))) {
1294  Int_t* nbins = new Int_t[fNdimensions];
1295  for (Int_t dim = fNdimensions - 1; dim >= 0; --dim) {
1296  nbins[dim] = GetAxis(dim)->GetNbins();
1297  bin[dim] = from % nbins[dim];
1298  from /= nbins[dim];
1299  }
1300 
1301  for (Long64_t i = 0; i < howmany; ++i) {
1302  if (!PrintBin(-1, bin, options))
1303  ++howmany;
1304  // Advance to next bin:
1305  ++bin[fNdimensions - 1];
1306  for (Int_t dim = fNdimensions - 1; dim >= 0; --dim) {
1307  if (bin[dim] >= nbins[dim]) {
1308  bin[dim] = 0;
1309  if (dim > 0) {
1310  ++bin[dim - 1];
1311  } else {
1312  howmany = -1; // aka "global break"
1313  }
1314  }
1315  }
1316  }
1317  delete [] nbins;
1318  } else {
1319  for (Long64_t i = from; i < from + howmany; ++i) {
1320  if (!PrintBin(i, bin, options))
1321  ++howmany;
1322  }
1323  }
1324  delete [] bin;
1325 }
1326 
1327 ////////////////////////////////////////////////////////////////////////////////
1328 /// Print a THnBase. If "option" contains:
1329 /// - 'a': print axis details
1330 /// - 'm': print memory usage
1331 /// - 's': print statistics
1332 /// - 'c': print its content, too (this can generate a LOT of output!)
1333 /// Other options are forwarded to PrintEntries().
1334 
1335 void THnBase::Print(Option_t* options) const
1336 {
1337  Bool_t optAxis = options && (strchr(options, 'A') || (strchr(options, 'a')));
1338  Bool_t optMem = options && (strchr(options, 'M') || (strchr(options, 'm')));
1339  Bool_t optStat = options && (strchr(options, 'S') || (strchr(options, 's')));
1340  Bool_t optContent = options && (strchr(options, 'C') || (strchr(options, 'c')));
1341 
1342  Printf("%s (*0x%lx): \"%s\" \"%s\"", IsA()->GetName(), (unsigned long)this, GetName(), GetTitle());
1343  Printf(" %d dimensions, %g entries in %lld filled bins", GetNdimensions(), GetEntries(), GetNbins());
1344 
1345  if (optAxis) {
1346  for (Int_t dim = 0; dim < fNdimensions; ++dim) {
1347  TAxis* axis = GetAxis(dim);
1348  Printf(" axis %d \"%s\": %d bins (%g..%g), %s bin sizes", dim,
1349  axis->GetTitle(), axis->GetNbins(), axis->GetXmin(), axis->GetXmax(),
1350  (axis->GetXbins() ? "variable" : "fixed"));
1351  }
1352  }
1353 
1354  if (optStat) {
1355  Printf(" %s error calculation", (GetCalculateErrors() ? "with" : "without"));
1356  if (GetCalculateErrors()) {
1357  Printf(" Sum(w)=%g, Sum(w^2)=%g", GetSumw(), GetSumw2());
1358  for (Int_t dim = 0; dim < fNdimensions; ++dim) {
1359  Printf(" axis %d: Sum(w*x)=%g, Sum(w*x^2)=%g", dim, GetSumwx(dim), GetSumwx2(dim));
1360  }
1361  }
1362  }
1363 
1364  if (optMem && InheritsFrom(THnSparse::Class())) {
1365  const THnSparse* hsparse = dynamic_cast<const THnSparse*>(this);
1366  Printf(" coordinates stored in %d chunks of %d entries\n %g of bins filled using %g of memory compared to an array",
1367  hsparse->GetNChunks(), hsparse->GetChunkSize(),
1368  hsparse->GetSparseFractionBins(), hsparse->GetSparseFractionMem());
1369  }
1370 
1371  if (optContent) {
1372  Printf(" BIN CONTENT:");
1373  PrintEntries(0, -1, options);
1374  }
1375 }
1376 
1377 
1378 ////////////////////////////////////////////////////////////////////////////////
1379 /// Browse a THnSparse: create an entry (ROOT::THnSparseBrowsable) for each
1380 /// dimension.
1381 
1383 {
1384  if (fBrowsables.IsEmpty()) {
1385  for (Int_t dim = 0; dim < fNdimensions; ++dim) {
1387  }
1389  }
1390 
1391  for (Int_t dim = 0; dim < fNdimensions; ++dim) {
1392  b->Add(fBrowsables[dim]);
1393  }
1394 }
1395 
1396 
1397 
1398 
1399 /** \class ROOT::Internal::THnBaseBinIter
1400  Iterator over THnBase bins (internal implementation).
1401 */
1402 
1403 /// Destruct a bin iterator.
1404 
1406  // Not much to do, but pin vtable
1407 }
1408 
1409 
1410 
1411 /** \class THnIter
1412  Iterator over THnBase bins
1413 */
1414 
1415 ClassImp(THnIter);
1416 
1418  // Destruct a bin iterator.
1419  delete fIter;
1420 }
1421 
1422 
1423 /** \class ROOT::Internal::THnBaseBrowsable
1424  TBrowser helper for THnBase.
1425 */
1426 
1427 
1429 
1430 ////////////////////////////////////////////////////////////////////////////////
1431 /// Construct a THnBaseBrowsable.
1432 
1434 fHist(hist), fAxis(axis), fProj(0)
1435 {
1436  TString axisName = hist->GetAxis(axis)->GetName();
1437  if (axisName.IsNull()) {
1438  axisName = TString::Format("axis%d", axis);
1439  }
1440 
1441  SetNameTitle(axisName,
1442  TString::Format("Projection on %s of %s", axisName.Data(),
1443  hist->IsA()->GetName()).Data());
1444 }
1445 
1446 ////////////////////////////////////////////////////////////////////////////////
1447 /// Destruct a THnBaseBrowsable.
1448 
1450 {
1451  delete fProj;
1452 }
1453 
1454 ////////////////////////////////////////////////////////////////////////////////
1455 /// Browse an axis of a THnBase, i.e. draw its projection.
1456 
1458 {
1459  if (!fProj) {
1460  fProj = fHist->Projection(fAxis);
1461  }
1462  fProj->Draw(b ? b->GetDrawOption() : "");
1463  gPad->Update();
1464 }
1465 
TString fTitle
Definition: TNamed.h:33
Abstract array base class.
Definition: TArray.h:31
Double_t ComputeIntegral()
Calculate the integral of the histogram.
Definition: THnBase.cxx:1166
Double_t GetBinError(const Int_t *idx) const
Definition: THnBase.h:160
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Double_t GetSparseFractionMem() const
Return the amount of used memory over memory that would be used by a non-sparse n-dimensional histogr...
Definition: THnSparse.cxx:866
virtual Int_t GetNcells() const
Definition: TH1.h:295
virtual Double_t GetEffectiveEntries() const
Number of effective entries of the histogram.
Definition: TH1.cxx:4302
An array of TObjects.
Definition: TObjArray.h:37
float xmin
Definition: THbookFile.cxx:93
virtual void Reset(Option_t *option="")=0
Bool_t GetCalculateErrors() const
Definition: THnBase.h:136
void Divide(const THnBase *h)
Divide this histogram by h this = this/(h) Note that if h has Sumw2 set, Sumw2 is automatically calle...
Definition: THnBase.cxx:848
Double_t Floor(Double_t x)
Definition: TMath.h:693
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:444
long long Long64_t
Definition: RtypesCore.h:69
Double_t * fIntegral
Definition: THnBase.h:53
void SetEntries(Double_t entries)
Definition: THnBase.h:165
TFitResultPtr Fit(TF1 *f1, Option_t *option="", Option_t *goption="")
Fit a THnSparse with function f.
Definition: THnBase.cxx:365
Double_t fTsumw2
Definition: THnBase.h:50
THnBase * RebinBase(Int_t group) const
Combine the content of "group" neighboring bins into a new bin and return the resulting THnBase...
Definition: THnBase.cxx:1057
enum THnBase::@73 fIntegralStatus
array with bin weight sums
virtual Long64_t GetNbins() const =0
const char Option_t
Definition: RtypesCore.h:62
void Multiply(const THnBase *h)
Multiply this histogram by histogram h this = this * h Note that if h has Sumw2 set, Sumw2 is automatically called for this if not already set.
Definition: THnBase.cxx:757
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:504
return c1
Definition: legend1.C:41
const Ssiz_t kNPOS
Definition: RtypesCore.h:111
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
const Double_t * GetArray() const
Definition: TArrayD.h:43
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
virtual ~THnIter()
Definition: THnBase.cxx:1417
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4882
Double_t GetBinContent(const Int_t *idx) const
Definition: THnBase.h:168
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition: TString.h:634
void Scale(Double_t c)
Scale contents and errors of this histogram by c: this = this * c It does not modify the histogram&#39;s ...
Definition: THnBase.cxx:611
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
virtual void GetBinXYZ(Int_t binglobal, Int_t &binx, Int_t &biny, Int_t &binz) const
Return binx, biny, binz corresponding to the global bin number globalbin see TH1::GetBin function abo...
Definition: TH1.cxx:4797
static THnBase * CreateHnAny(const char *name, const char *title, const TH1 *h1, Bool_t sparse, Int_t chunkSize=1024 *16)
Create a THn / THnSparse object from a histogram deriving from TH1.
Definition: THnBase.cxx:197
void GetRandom(Double_t *rand, Bool_t subBinRandom=kTRUE)
Generate an n-dimensional random tuple based on the histogrammed distribution.
Definition: THnBase.cxx:394
Basic string class.
Definition: TString.h:131
#define f(i)
Definition: RSha256.hxx:104
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
Bool_t IsEmpty() const
Definition: TObjArray.h:71
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TObject * At(Int_t idx) const
Definition: TObjArray.h:166
Int_t GetNChunks() const
Definition: THnSparse.h:88
virtual void SetNameTitle(const char *name, const char *title)
Set all the TNamed parameters (name and title).
Definition: TNamed.cxx:154
TString & Insert(Ssiz_t pos, const char *s)
Definition: TString.h:644
void AddBinContent(const Int_t *x, Double_t v=1.)
Definition: THnBase.h:164
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
TObjArray fAxes
Definition: THnBase.h:46
int Nostore
Definition: Foption.h:41
void Reset()
Definition: TCollection.h:252
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
Double_t fTsumw
Definition: THnBase.h:49
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:514
Double_t fEntries
browser-helpers for each axis
Definition: THnBase.h:48
THnBaseBrowsable(THnBase *hist, Int_t axis)
Construct a THnBaseBrowsable.
Definition: THnBase.cxx:1433
Double_t GetXmin() const
Definition: TAxis.h:133
Double_t x[n]
Definition: legend1.C:17
TBrowser helper for THnBase.
Definition: THnBase.h:270
static Int_t FitOptionsMake(Option_t *option, Foption_t &Foption)
Decode string choptin and fill fitOption structure.
Definition: TH1.cxx:4475
Efficient multidimensional histogram.
Definition: THnSparse.h:36
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
Definition: TString.cxx:2311
void Class()
Definition: Class.C:29
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
int Chi2
Definition: Foption.h:32
THnBase * CloneEmpty(const char *name, const char *title, const TObjArray *axes, Bool_t keepTargetAxis) const
Create a new THnBase object that is of the same type as *this, but with dimensions and bins given by ...
Definition: THnBase.cxx:80
Int_t GetCoord(Int_t dim) const
Definition: THnBase.h:317
virtual void Reserve(Long64_t)
Definition: THnBase.h:99
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
void PrintEntries(Long64_t from=0, Long64_t howmany=-1, Option_t *options=0) const
Print "howmany" entries starting at "from".
Definition: THnBase.cxx:1285
void SetBinContent(const Int_t *idx, Double_t v)
Definition: THnBase.h:178
Double_t GetSumw2() const
Definition: THnBase.h:185
virtual ~THnBaseBinIter()
Destruct a bin iterator.
Definition: THnBase.cxx:1405
Double_t GetSumw() const
Definition: THnBase.h:184
TObjArray fBrowsables
Definition: THnBase.h:47
Ssiz_t First(char c) const
Find first occurrence of a character c.
Definition: TString.cxx:499
Bool_t CheckConsistency(const THnBase *h, const char *tag) const
Consistency check on (some of) the parameters of two histograms (for operations). ...
Definition: THnBase.cxx:987
static constexpr double s
void Browse(TBrowser *b)
Browse an axis of a THnBase, i.e. draw its projection.
Definition: THnBase.cxx:1457
TH1F * h1
Definition: legend1.C:5
virtual void AddBinContent(Int_t bin)
Increment bin content by 1.
Definition: TH1.cxx:1200
virtual void ResetStats()
Reset the statistics including the number of entries and replace with values calculates from bin cont...
Definition: TH1.cxx:7382
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition: TH1.cxx:8628
virtual Double_t GetBinError2(Long64_t linidx) const =0
Bool_t PrintBin(Long64_t idx, Int_t *coord, Option_t *options) const
Print one bin.
Definition: THnBase.cxx:1241
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
Using a TBrowser one can browse all ROOT objects.
Definition: TBrowser.h:37
RooArgSet S(const RooAbsArg &v1)
void AddInternal(const THnBase *h, Double_t c, Bool_t rebinned)
Add() implementation for both rebinned histograms and those with identical binning.
Definition: THnBase.cxx:635
#define F(x, y, z)
virtual void SetRange(Int_t first=0, Int_t last=0)
Set the viewing range for the axis from bin first to last.
Definition: TAxis.cxx:903
virtual void InitStorage(Int_t *nbins, Int_t chunkSize)=0
virtual void SetFilledBins(Long64_t)
Definition: THnBase.h:100
static double C[]
TAxis * GetAxis(Int_t dim) const
Definition: THnBase.h:125
#define R__THNBCASE(TAG)
virtual void AddAtAndExpand(TObject *obj, Int_t idx)
Add object at position idx.
Definition: TObjArray.cxx:234
void Print(Option_t *option="") const
Print a THnBase.
Definition: THnBase.cxx:1335
const char * GetTitle() const
Returns title of object.
Definition: TAxis.h:129
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:455
Class to manage histogram axis.
Definition: TAxis.h:30
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:304
Int_t GetSize() const
Definition: TArray.h:47
Long64_t Next(Int_t *coord=0)
Return the next bin&#39;s index.
Definition: THnBase.h:313
virtual Long64_t GetBin(const Int_t *idx) const =0
void Init(const char *name, const char *title, const TObjArray *axes, Bool_t keepTargetAxis, Int_t chunkSize=1024 *16)
Initialize axes and name.
Definition: THnBase.cxx:96
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
Definition: TFitResultPtr.h:31
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:443
Collection abstract base class.
Definition: TCollection.h:63
virtual ~THnBase()
Destruct a THnBase.
Definition: THnBase.cxx:69
Int_t GetEntriesFast() const
Definition: TObjArray.h:64
Double_t GetSumwx2(Int_t dim) const
Definition: THnBase.h:187
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
char * Form(const char *fmt,...)
Ssiz_t Length() const
Definition: TString.h:405
static void RejectPoint(Bool_t reject=kTRUE)
Static function to set the global flag to reject points the fgRejectPoint global flag is tested by al...
Definition: TF1.cxx:3650
int Like
Definition: Foption.h:34
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition: TClass.h:75
TAxis * GetYaxis()
Definition: TH1.h:317
float xmax
Definition: THbookFile.cxx:93
void SetBinEdges(Int_t idim, const Double_t *bins)
Set the axis # of bins and bin limits on dimension idim.
Definition: THnBase.cxx:1005
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
virtual Int_t GetBin(Int_t binx, Int_t biny=0, Int_t binz=0) const
Return Global bin number corresponding to binx,y,z.
Definition: TH1.cxx:4784
virtual void Rebuild(Option_t *option="")
Using the current bin info, recompute the arrays for contents and errors.
Definition: TH1.cxx:6691
virtual void AddAt(TObject *obj, Int_t idx)
Add object at position ids.
Definition: TObjArray.cxx:253
Double_t GetEntries() const
Definition: THnBase.h:133
#define h(i)
Definition: RSha256.hxx:106
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t At(Int_t i) const
Definition: TArrayD.h:79
Double_t GetSumwx(Int_t dim) const
Definition: THnBase.h:186
TString & Remove(Ssiz_t pos)
Definition: TString.h:668
void AddRange(unsigned int icoord, double xmin, double xmax)
add a range [xmin,xmax] for the new coordinate icoord Adding a range does not delete existing one...
Definition: DataRange.cxx:94
int Ssiz_t
Definition: RtypesCore.h:63
#define d(i)
Definition: RSha256.hxx:102
void SetBinError(const Int_t *idx, Double_t e)
Definition: THnBase.h:162
virtual Bool_t IsEmpty() const
Definition: TCollection.h:186
return c2
Definition: legend2.C:14
virtual void SetBinError2(Long64_t bin, Double_t e2)=0
#define ClassImp(name)
Definition: Rtypes.h:365
class describing the range in the coordinates it supports multiple range in a coordinate.
Definition: DataRange.h:34
double Double_t
Definition: RtypesCore.h:55
Bool_t IsInRange(Int_t *coord) const
Check whether bin coord is in range, as defined by TAxis::SetRange().
Definition: THnBase.cxx:429
void Add(const THnBase *h, Double_t c=1.)
Add contents of h scaled by c to this histogram: this = this + c * h Note that if h has Sumw2 set...
Definition: THnBase.cxx:696
void Printf(const char *fmt,...)
Long64_t Merge(TCollection *list)
Merge this with a list of THnBase&#39;s.
Definition: THnBase.cxx:722
int type
Definition: TGX11.cxx:120
THnBase()
Definition: THnBase.h:65
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
The TH1 histogram class.
Definition: TH1.h:56
virtual Double_t GetEntries() const
Return the current number of entries.
Definition: TH1.cxx:4277
~THnBaseBrowsable()
Destruct a THnBaseBrowsable.
Definition: THnBase.cxx:1449
virtual void AddBinError2(Long64_t bin, Double_t e2)=0
Double_t GetSparseFractionBins() const
Return the amount of filled bins over all bins.
Definition: THnSparse.cxx:855
virtual void Sumw2()=0
static TClass * GetClass(const char *name, Bool_t load=kTRUE, Bool_t silent=kFALSE)
Static method returning pointer to TClass of the specified class name.
Definition: TClass.cxx:2906
TObject * ProjectionAny(Int_t ndim, const Int_t *dim, Bool_t wantNDim, Option_t *option="") const
Project all bins into a ndim-dimensional THn / THnSparse (whatever *this is) or if (ndim < 4 and !wan...
Definition: THnBase.cxx:455
Bool_t IsNull() const
Definition: TString.h:402
Int_t GetNdimensions() const
Definition: THnBase.h:135
TAxis * GetZaxis()
Definition: TH1.h:318
void Browse(TBrowser *b)
Browse a THnSparse: create an entry (ROOT::THnSparseBrowsable) for each dimension.
Definition: THnBase.cxx:1382
Mother of all ROOT objects.
Definition: TObject.h:37
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width.
Definition: TAxis.cxx:526
Int_t fNdimensions
Definition: THnBase.h:45
Bool_t HaveSkippedBin() const
Definition: THnBase.h:318
Iterator over THnBase bins.
Definition: THnBase.h:303
1-Dim function class
Definition: TF1.h:211
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
#define gPad
Definition: TVirtualPad.h:286
TObjArray * GetListOfAxes()
Definition: THnBase.h:123
#define c(i)
Definition: RSha256.hxx:101
virtual void SetEntries(Double_t n)
Definition: TH1.h:381
Multidimensional histogram base.
Definition: THnBase.h:43
void ResetBit(UInt_t f)
Definition: TObject.h:171
TFitResultPtr FitObject(TH1 *h1, TF1 *f1, Foption_t &option, const ROOT::Math::MinimizerOptions &moption, const char *goption, ROOT::Fit::DataRange &range)
fitting function for a TH1 (called from TH1::Fit)
Definition: HFitImpl.cxx:966
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Int_t GetNbins() const
Definition: TAxis.h:121
#define I(x, y, z)
TH1 * CreateHist(const char *name, const char *title, const TObjArray *axes, Bool_t keepTargetAxis) const
Create an empty histogram with name and title with a given set of axes.
Definition: THnBase.cxx:144
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:357
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition: TCollection.h:182
virtual void Set(Int_t nbins, Double_t xmin, Double_t xmax)
Initialize axis with fix bins.
Definition: TAxis.cxx:717
Int_t GetChunkSize() const
Definition: THnSparse.h:87
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
const Bool_t kTRUE
Definition: RtypesCore.h:87
Double_t GetXmax() const
Definition: TAxis.h:134
void RebinnedAdd(const THnBase *h, Double_t c=1.)
Add contents of h scaled by c to this histogram: this = this + c * h Note that if h has Sumw2 set...
Definition: THnBase.cxx:712
Int_t CeilNint(Double_t x)
Definition: TMath.h:689
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMathBase.h:278
char name[80]
Definition: TGX11.cxx:109
const TArrayD * GetXbins() const
Definition: TAxis.h:130
void ResetBase(Option_t *option="")
Clear the histogram.
Definition: THnBase.cxx:1152
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
void SetTitle(const char *title)
Change (i.e.
Definition: THnBase.cxx:1021
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8485
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
const char * Data() const
Definition: TString.h:364