Logo ROOT   6.18/05
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"
34
35
36/** \class THnBase
37 \ingroup Hist
38Multidimensional histogram base.
39Defines 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
51THnBase::THnBase(const char* name, const char* title, Int_t dim,
52 const Int_t* nbins, const Double_t* xmin, const Double_t* xmax):
53TNamed(name, title), fNdimensions(dim), fAxes(dim), fBrowsables(dim),
54fEntries(0), fTsumw(0), fTsumw2(-1.), fTsumwx(dim), fTsumwx2(dim),
55fIntegral(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);
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
80THnBase* 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
96void 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
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
144TH1* 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
197THnBase* 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) \
224if (sparse) { \
225s = new _NAME2_(THnSparse,TAG)(name, title, ndim, nbins, \
226minRange, maxRange, chunkSize); \
227} else { \
228s = new _NAME2_(THn,TAG)(name, title, ndim, nbins, \
229minRange, maxRange); \
230} \
231break;
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
267THnBase* 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
331void 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
394void THnBase::GetRandom(Double_t *rand, Bool_t subBinRandom /* = kTRUE */)
395{
396 // check whether the integral array is valid
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{
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
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
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
777 Double_t v1 = GetBinContent(i);
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
871 Double_t v1 = GetBinContent(i);
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
905void 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
987Bool_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
1005void 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
1021void 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) {
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
1152void 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) {
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)
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
1227void 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
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
1285void 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
1335void 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
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
1434fHist(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
void Class()
Definition: Class.C:29
SVector< double, 2 > v
Definition: Dict.h:5
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
#define h(i)
Definition: RSha256.hxx:106
const Ssiz_t kNPOS
Definition: RtypesCore.h:111
int Int_t
Definition: RtypesCore.h:41
int Ssiz_t
Definition: RtypesCore.h:63
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
long long Long64_t
Definition: RtypesCore.h:69
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
int type
Definition: TGX11.cxx:120
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
#define R__THNBCASE(TAG)
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
char * Form(const char *fmt,...)
void Printf(const char *fmt,...)
#define gPad
Definition: TVirtualPad.h:286
class describing the range in the coordinates it supports multiple range in a coordinate.
Definition: DataRange.h:34
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
virtual ~THnBaseBinIter()
Destruct a bin iterator.
Definition: THnBase.cxx:1405
TBrowser helper for THnBase.
Definition: THnBase.h:270
~THnBaseBrowsable()
Destruct a THnBaseBrowsable.
Definition: THnBase.cxx:1449
THnBaseBrowsable(THnBase *hist, Int_t axis)
Construct a THnBaseBrowsable.
Definition: THnBase.cxx:1433
void Browse(TBrowser *b)
Browse an axis of a THnBase, i.e. draw its projection.
Definition: THnBase.cxx:1457
Double_t At(Int_t i) const
Definition: TArrayD.h:79
const Double_t * GetArray() const
Definition: TArrayD.h:43
Abstract array base class.
Definition: TArray.h:31
Int_t GetSize() const
Definition: TArray.h:47
Class to manage histogram axis.
Definition: TAxis.h:30
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
const TArrayD * GetXbins() const
Definition: TAxis.h:130
Double_t GetXmax() const
Definition: TAxis.h:134
@ kAxisRange
Definition: TAxis.h:61
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:504
virtual void Set(Int_t nbins, Double_t xmin, Double_t xmax)
Initialize axis with fix bins.
Definition: TAxis.cxx:717
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:455
Double_t GetXmin() const
Definition: TAxis.h:133
Int_t GetNbins() const
Definition: TAxis.h:121
const char * GetTitle() const
Returns title of object.
Definition: TAxis.h:129
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 Double_t GetBinWidth(Int_t bin) const
Return bin width.
Definition: TAxis.cxx:526
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:514
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:444
Using a TBrowser one can browse all ROOT objects.
Definition: TBrowser.h:37
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition: TClass.h:75
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:2895
Collection abstract base class.
Definition: TCollection.h:63
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
virtual Bool_t IsEmpty() const
Definition: TCollection.h:186
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition: TCollection.h:182
1-Dim function class
Definition: TF1.h:211
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
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
Definition: TFitResultPtr.h:31
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
The TH1 histogram class.
Definition: TH1.h:56
virtual Double_t GetEffectiveEntries() const
Number of effective entries of the histogram.
Definition: TH1.cxx:4302
virtual void Rebuild(Option_t *option="")
Using the current bin info, recompute the arrays for contents and errors.
Definition: TH1.cxx:6684
TAxis * GetZaxis()
Definition: TH1.h:318
virtual void AddBinContent(Int_t bin)
Increment bin content by 1.
Definition: TH1.cxx:1200
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8476
static Int_t FitOptionsMake(Option_t *option, Foption_t &Foption)
Decode string choptin and fill fitOption structure.
Definition: TH1.cxx:4475
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
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
virtual Int_t GetNcells() const
Definition: TH1.h:295
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 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:8619
TAxis * GetYaxis()
Definition: TH1.h:317
virtual Double_t GetEntries() const
Return the current number of entries.
Definition: TH1.cxx:4277
virtual void ResetStats()
Reset the statistics including the number of entries and replace with values calculates from bin cont...
Definition: TH1.cxx:7374
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4882
virtual void SetEntries(Double_t n)
Definition: TH1.h:381
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:289
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:302
Multidimensional histogram base.
Definition: THnBase.h:43
virtual void Sumw2()=0
void SetEntries(Double_t entries)
Definition: THnBase.h:165
virtual void SetFilledBins(Long64_t)
Definition: THnBase.h:100
virtual void InitStorage(Int_t *nbins, Int_t chunkSize)=0
Double_t GetBinError(const Int_t *idx) const
Definition: THnBase.h:160
enum THnBase::@77 fIntegralStatus
array with bin weight sums
void SetBinError(const Int_t *idx, Double_t e)
Definition: THnBase.h:162
Double_t fEntries
browser-helpers for each axis
Definition: THnBase.h:48
Double_t GetSumw2() const
Definition: THnBase.h:185
Bool_t IsInRange(Int_t *coord) const
Check whether bin coord is in range, as defined by TAxis::SetRange().
Definition: THnBase.cxx:429
TFitResultPtr Fit(TF1 *f1, Option_t *option="", Option_t *goption="")
Fit a THnSparse with function f.
Definition: THnBase.cxx:365
void AddBinContent(const Int_t *x, Double_t v=1.)
Definition: THnBase.h:164
void Scale(Double_t c)
Scale contents and errors of this histogram by c: this = this * c It does not modify the histogram's ...
Definition: THnBase.cxx:611
virtual void SetBinError2(Long64_t bin, Double_t e2)=0
TObjArray * GetListOfAxes()
Definition: THnBase.h:123
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
@ kNoInt
Definition: THnBase.h:55
@ kValidInt
Definition: THnBase.h:56
Bool_t PrintBin(Long64_t idx, Int_t *coord, Option_t *options) const
Print one bin.
Definition: THnBase.cxx:1241
Double_t GetSumwx(Int_t dim) const
Definition: THnBase.h:186
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
Int_t GetNdimensions() const
Definition: THnBase.h:135
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 ResetBase(Option_t *option="")
Clear the histogram.
Definition: THnBase.cxx:1152
virtual Long64_t GetNbins() const =0
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
TObjArray fAxes
Definition: THnBase.h:46
void Browse(TBrowser *b)
Browse a THnSparse: create an entry (ROOT::THnSparseBrowsable) for each dimension.
Definition: THnBase.cxx:1382
void SetBinEdges(Int_t idim, const Double_t *bins)
Set the axis # of bins and bin limits on dimension idim.
Definition: THnBase.cxx:1005
Bool_t GetCalculateErrors() const
Definition: THnBase.h:136
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
Double_t GetBinContent(const Int_t *idx) const
Definition: THnBase.h:168
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 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
Double_t fTsumw2
Definition: THnBase.h:50
virtual void Reserve(Long64_t)
Definition: THnBase.h:99
virtual ~THnBase()
Destruct a THnBase.
Definition: THnBase.cxx:69
THnBase()
Definition: THnBase.h:65
Double_t GetEntries() const
Definition: THnBase.h:133
virtual Double_t GetBinError2(Long64_t linidx) const =0
virtual void Reset(Option_t *option="")=0
void SetBinContent(const Int_t *idx, Double_t v)
Definition: THnBase.h:178
void Multiply(const THnBase *h)
Multiply this histogram by histogram h this = this * h Note that if h has Sumw2 set,...
Definition: THnBase.cxx:757
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 GetSumwx2(Int_t dim) const
Definition: THnBase.h:187
TObjArray fBrowsables
Definition: THnBase.h:47
void Print(Option_t *option="") const
Print a THnBase.
Definition: THnBase.cxx:1335
Double_t GetSumw() const
Definition: THnBase.h:184
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 !...
Definition: THnBase.cxx:455
void GetRandom(Double_t *rand, Bool_t subBinRandom=kTRUE)
Generate an n-dimensional random tuple based on the histogrammed distribution.
Definition: THnBase.cxx:394
TAxis * GetAxis(Int_t dim) const
Definition: THnBase.h:125
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
Long64_t Merge(TCollection *list)
Merge this with a list of THnBase's.
Definition: THnBase.cxx:722
Double_t fTsumw
Definition: THnBase.h:49
virtual Long64_t GetBin(const Int_t *idx) const =0
Double_t ComputeIntegral()
Calculate the integral of the histogram.
Definition: THnBase.cxx:1166
Double_t * fIntegral
Definition: THnBase.h:53
virtual void AddBinError2(Long64_t bin, Double_t e2)=0
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
void SetTitle(const char *title)
Change (i.e.
Definition: THnBase.cxx:1021
Int_t fNdimensions
Definition: THnBase.h:45
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
Iterator over THnBase bins.
Definition: THnBase.h:303
Long64_t Next(Int_t *coord=0)
Return the next bin's index.
Definition: THnBase.h:313
Bool_t HaveSkippedBin() const
Definition: THnBase.h:318
Int_t GetCoord(Int_t dim) const
Definition: THnBase.h:317
virtual ~THnIter()
Definition: THnBase.cxx:1417
Efficient multidimensional histogram.
Definition: THnSparse.h:36
Double_t GetSparseFractionBins() const
Return the amount of filled bins over all bins.
Definition: THnSparse.cxx:855
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
Int_t GetNChunks() const
Definition: THnSparse.h:88
Int_t GetChunkSize() const
Definition: THnSparse.h:87
void Reset()
Definition: TCollection.h:252
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
TString fTitle
Definition: TNamed.h:33
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual void SetNameTitle(const char *name, const char *title)
Set all the TNamed parameters (name and title).
Definition: TNamed.cxx:154
An array of TObjects.
Definition: TObjArray.h:37
Int_t GetEntriesFast() const
Definition: TObjArray.h:64
virtual void AddAtAndExpand(TObject *obj, Int_t idx)
Add object at position idx.
Definition: TObjArray.cxx:234
Bool_t IsEmpty() const
Definition: TObjArray.h:71
virtual void AddAt(TObject *obj, Int_t idx)
Add object at position ids.
Definition: TObjArray.cxx:253
TObject * At(Int_t idx) const
Definition: TObjArray.h:166
Mother of all ROOT objects.
Definition: TObject.h:37
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:357
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:443
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
void ResetBit(UInt_t f)
Definition: TObject.h:171
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
Basic string class.
Definition: TString.h:131
Ssiz_t Length() const
Definition: TString.h:405
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
TString & Insert(Ssiz_t pos, const char *s)
Definition: TString.h:644
Ssiz_t First(char c) const
Find first occurrence of a character c.
Definition: TString.cxx:499
const char * Data() const
Definition: TString.h:364
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
Bool_t IsNull() const
Definition: TString.h:402
TString & Remove(Ssiz_t pos)
Definition: TString.h:668
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
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition: TString.h:634
return c1
Definition: legend1.C:41
Double_t x[n]
Definition: legend1.C:17
TH1F * h1
Definition: legend1.C:5
return c2
Definition: legend2.C:14
#define F(x, y, z)
#define I(x, y, z)
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:964
static double C[]
RooArgSet S(const RooAbsArg &v1)
static constexpr double s
Double_t Floor(Double_t x)
Definition: TMath.h:691
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
Int_t CeilNint(Double_t x)
Definition: TMath.h:687
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMathBase.h:278
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
int Like
Definition: Foption.h:34
int Nostore
Definition: Foption.h:41
int Chi2
Definition: Foption.h:32