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