Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TUnfoldBinning.cxx
Go to the documentation of this file.
1// Author: Stefan Schmitt
2// DESY, 10/08/11
3
4// Version 17.9, parallel to changes in TUnfold
5//
6// History:
7// Version 17.8, bug fix in GetNonemptyNode() and non-const access of tree
8// Version 17.7, bug fix in ExtractHistogram
9// Version 17.6, bug fix to avoid possible crash in method
10// CreateHistogramOfMigrations(). Bug fix with NaN in GetGlobalBinNumber()
11// Version 17.5, in parallel to changes in TUnfold
12// Version 17.4, bug fix with error handling
13// Version 17.3, bug fix with underflow/overflow bins
14// Version 17.2, with XML support, bug fix with bin map creation,
15// isPeriodic option for neighbour bins
16// Version 17.1, in parallel to changes in TUnfold
17// Version 17.0, initial version, numbered in parallel to TUnfold
18
19/** \class TUnfoldBinning
20Binning schemes for use with the unfolding algorithm TUnfoldDensity.
21
22Binning schemes are used to map analysis bins on a single histogram
23axis and back. The analysis bins may include unconnected bins (e.g
24nuisances for background normalisation) or various multidimensional
25histograms (signal bins, differential background normalisation bins, etc).
26<br/>
27If you use this software, please consider the following citation
28<br/>
29<b>S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]</b>
30<br/>
31Detailed documentation and updates are available on
32http://www.desy.de/~sschmitt
33
34<h3>Functionality</h3>
35
36The TUnfoldBinning objects are connected by a tree-like structure.
37The structure does not hold any data, but is only responsible for
38arranging the analysis bins in the proper order.
39Each node of the tree is responsible for a group of bins. That group
40may consist of
41<ul>
42<li> several unconnected bins, each with a dedicated name.</li>
43<li> bins organized in a multidimensional distribution, defined by a
44set of axes. The axes are defined by a number of bins N and by (N+1)
45bin borders. In addition to the N bins inside there may be an underflow and an
46overflow bin</li>
47</ul>
48Each bin has a "global" bin number, which can be found using the
49GetGlobalBinNumber() methods. The global bin number 0 is reserved and
50corresponds to the case where no bin is found in the
51TUnfoldBinning tree.
52
53<h3>Use in the analysis</h3>
54Booking histograms:
55<ul>
56<li>Define binning schemes on detector level and on truth level. This
57can be done using the XML language, use the class TUnfoldBinningXML to
58read the binning scheme. The TUnfoldBinning objects can be written to
59a root file, preferentially together with the corresponding histograms.</li>
60<li>For Monte Carlo, book histograms for the response matrix (detector
61vs truth level) using the
62method CreateHistogramOfMigrations()</li>
63<li>For data and background, book histograms using the
64"detector level" binning scheme and the method CreateHistogram()</li>
65<li>(if required) for the data covarianve matrix, book a histogram using the
66"detector level" binning scheme and the method CreateErrorMatrixHistogram()</li>
67<li>For truth histograms, book histograms using the
68"truth level" binning scheme and the method CreateHistogram()</li>
69</ul>
70The histograms which are booked have all analysis bins arranged on one
71axis (global bin number). TUnfoldBinning provides methods to locate
72the global bin number:
73<ul>
74<li>Use the method FindNode() to locate a group of bins (e.g. signal,
75control distribution, etc) by their name, then:</li>
76<li>Use the method GetGlobalBinNumber() to locate a bin in a
77distribution, then:</li>
78<li>Use the TH1::Fill() method and the bin number to fill the
79appropriate bin in one of the histograms booked above.</li>
80</ul>
81Unfolding: Specify the response matrix and the binning schemes when
82constructing a TUnfoldDensity object. Tell TUnfoldDensity about the
83data, bakcground, systematic error histograms using the corresponding
84methods of class TUnfoldDensity. Then run the unfolding. Use the
85GetXXX() methods to retreive the unfolding results into properly
86binned multidimensional histograms.
87*/
88
89/*
90 This file is part of TUnfold.
91
92 TUnfold is free software: you can redistribute it and/or modify
93 it under the terms of the GNU General Public License as published by
94 the Free Software Foundation, either version 3 of the License, or
95 (at your option) any later version.
96
97 TUnfold is distributed in the hope that it will be useful,
98 but WITHOUT ANY WARRANTY; without even the implied warranty of
99 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
100 GNU General Public License for more details.
101
102 You should have received a copy of the GNU General Public License
103 along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
104*/
105
106
107#include "TUnfoldBinningXML.h"
108#include <TVectorD.h>
109#include <TAxis.h>
110#include <TString.h>
111#include <TMath.h>
112#include <TF1.h>
113#include <TH1D.h>
114#include <TH2D.h>
115#include <TH3D.h>
116#include <TIterator.h>
117#include <iomanip>
118
119// #define DEBUG
120
121using std::ostream;
122
124
125////////////////////////////////////////////////////////////////////////
127{
128 // delete all children
129 while(childNode) delete childNode;
130 // remove this node from the tree
131 if(GetParentNode() && (GetParentNode()->GetChildNode()==this)) {
132 parentNode->childNode=nextNode;
133 }
134 if(GetPrevNode()) prevNode->nextNode=nextNode;
135 if(GetNextNode()) nextNode->prevNode=prevNode;
136 delete fAxisList;
137 delete fAxisLabelList;
138 if(fBinFactorFunction) {
139 if(!dynamic_cast<TF1 *>(fBinFactorFunction))
140 delete fBinFactorFunction;
141 }
142}
143
144/********************* setup **************************/
145
146////////////////////////////////////////////////////////////////////////
147/// initialize variables for a given number of bins
149{
150 parentNode=nullptr;
151 childNode=nullptr;
152 nextNode=nullptr;
153 prevNode=nullptr;
154 fAxisList=new TObjArray();
159 fHasOverflow=0;
160 fDistributionSize=nBins;
161 fBinFactorFunction=nullptr;
163}
164
165////////////////////////////////////////////////////////////////////////
166/// update fFirstBin and fLastBin members of this node and its children
167///
168/// \param[in] startWithRootNode if true, start the update with the root node
170{
171 if(startWithRootNode) {
173 }
174 if(GetPrevNode()) {
175 // if this is not the first node in a sequence,
176 // start with the end bin of the previous node
178 } else if(GetParentNode()) {
179 // if this is the first node in a sequence but has a parent,
180 // start with the end bin of the parent's distribution
183 } else {
184 // if this is the top level node, the first bin number is 1
185 fFirstBin=1;
186 // ... unless the top level node is the only node
187 // ... with dimension=1
188 // ... and there are no child nodes
189 // ... and there is an underflow bin
190 if((!GetChildNode())&&(GetDistributionDimension()==1)&&
191 (fHasUnderflow==1)) {
192 fFirstBin=0;
193 }
194 }
196 // now update count for all children
197 for(TUnfoldBinning *node=childNode;node;node=node->nextNode) {
198 fLastBin=node->UpdateFirstLastBin(kFALSE);
199 }
200 return fLastBin;
201}
202
203////////////////////////////////////////////////////////////////////////
204/// create a new node without axis
205///
206/// \param[in] name identifier of the node
207/// \param[in] nBin number of unconnected bins (could be zero)
208/// \param[in] binNames (optional) names of the bins separated by ';'
210(const char *name,Int_t nBins,const char *binNames)
211 : TNamed(name ? name : "",name ? name : "")
212{
213 Initialize(nBins);
214 if(binNames) {
215 TString nameString(binNames);
216 delete fAxisLabelList;
217 fAxisLabelList=nameString.Tokenize(";");
218 }
220}
221
222////////////////////////////////////////////////////////////////////////
223/// create a new node containing a distribution with one axis
224///
225/// \param[in] axis the axis to represent
226/// \param[in] includeUnderflow true if underflow bin should be included
227/// \param[in] includeOverflow true if overflow bin should be included
229(const TAxis &axis,Int_t includeUnderflow,Int_t includeOverflow)
230 : TNamed(axis.GetName(),axis.GetTitle())
231{
232 Initialize(0);
233 AddAxis(axis,includeUnderflow,includeOverflow);
235}
236
237////////////////////////////////////////////////////////////////////////
238/// add a new binning node as last last child of this node
239///
240/// \param[in] name name of the node
241/// \param[in] nBin number of extra bins
242/// \param[in] binNames (optional) names of the bins sepatared by ';'
243///
244/// this is a shortcut for AddBinning(new TUnfoldBinning(name,nBins,binNames))
246(const char *name,Int_t nBins,const char *binNames)
247{
248 return AddBinning(new TUnfoldBinning(name,nBins,binNames));
249}
250
251////////////////////////////////////////////////////////////////////////
252/// add a TUnfoldBinning as the last child of this node
253///
254/// \param[in] binning the new binning to be added
255///
256/// return value: if succeeded, return "binning"
257/// otherwise return 0
259{
260 TUnfoldBinning *r=nullptr;
261 if(binning->GetParentNode()) {
262 Error("AddBinning",
263 "binning \"%s\" already has parent \"%s\", can not be added to %s",
264 (char *)binning->GetName(),
265 (char *)binning->GetParentNode()->GetName(),
266 (char *)GetName());
267 } else if(binning->GetPrevNode()) {
268 Error("AddBinning",
269 "binning \"%s\" has previous node \"%s\", can not be added to %s",
270 (char *)binning->GetName(),
271 (char *)binning->GetPrevNode()->GetName(),
272 (char *)GetName());
273 } else if(binning->GetNextNode()) {
274 Error("AddBinning",
275 "binning \"%s\" has next node \"%s\", can not be added to %s",
276 (char *)binning->GetName(),
277 (char *)binning->GetNextNode()->GetName(),
278 (char *)GetName());
279 } else {
280 r=binning;
281 binning->parentNode=this;
282 if(childNode) {
284 // find last child
285 while(child->nextNode) {
286 child=child->nextNode;
287 }
288 // add as last child
289 child->nextNode=r;
290 r->prevNode=child;
291 } else {
292 childNode=r;
293 }
295 r=binning;
296 }
297 return r;
298}
299
300////////////////////////////////////////////////////////////////////////
301/// add an axis with equidistant bins
302///
303/// \param[in] name name of the axis
304/// \param[in] nBin number of bins
305/// \param[in] xMin lower edge of the first bin
306/// \param[in] xMax upper edge of the last bin
307/// \param[in] hasUnderflow decide whether the axis has an underflow bin
308/// \param[in] hasOverflow decide whether the axis has an overflow bin
309///
310/// returns true if the axis has been added
312(const char *name,Int_t nBin,Double_t xMin,Double_t xMax,
313 Bool_t hasUnderflow,Bool_t hasOverflow)
314{
316 if(nBin<=0) {
317 Fatal("AddAxis","number of bins %d is not positive",
318 nBin);
319 } else if((!TMath::Finite(xMin))||(!TMath::Finite(xMax))||
320 (xMin>=xMax)) {
321 Fatal("AddAxis","xmin=%f required to be smaller than xmax=%f",
322 xMin,xMax);
323 } else {
324 Double_t *binBorders=new Double_t[nBin+1];
325 Double_t x=xMin;
326 Double_t dx=(xMax-xMin)/nBin;
327 for(Int_t i=0;i<=nBin;i++) {
328 binBorders[i]=x+i*dx;
329 }
330 r=AddAxis(name,nBin,binBorders,hasUnderflow,hasOverflow);
331 delete [] binBorders;
332 }
333 return r;
334}
335
336////////////////////////////////////////////////////////////////////////
337/// add an axis to the distribution, using the TAxis as blueprint
338///
339/// \param[in] axis blueprint of the axis
340/// \param[in] hasUnderflow decide whether the underflow bin should be included
341/// \param[in] hasOverflow decide whether the overflow bin should be included
342///
343/// returns true if the axis has been added
344/// </br/>
345/// Note: axis labels are not imported
347(const TAxis &axis,Bool_t hasUnderflow,Bool_t hasOverflow)
348{
349 Int_t nBin=axis.GetNbins();
350 Double_t *binBorders=new Double_t[nBin+1];
351 for(Int_t i=0;i<nBin;i++) {
352 binBorders[i]=axis.GetBinLowEdge(i+1);
353 }
354 binBorders[nBin]=axis.GetBinUpEdge(nBin);
355 Bool_t r=AddAxis(axis.GetTitle(),nBin,binBorders,hasUnderflow,hasOverflow);
356 delete [] binBorders;
357 return r;
358}
359
360////////////////////////////////////////////////////////////////////////
361/// add an axis with the specified bin borders
362///
363/// \param[in] name name of the axis
364/// \param[in] nBin number of bins
365/// \param[in] binBorders array of bin borders, with nBin+1 elements
366/// \param[in] hasUnderflow decide whether the axis has an underflow bin
367/// \param[in] hasOverflow decide whether the axis has an overflow bin
368///
369/// returns true if the axis has been added
370
372(const char *name,Int_t nBin,const Double_t *binBorders,
373 Bool_t hasUnderflow,Bool_t hasOverflow)
374{
376 if(HasUnconnectedBins()) {
377 Fatal("AddAxis","node already has %d bins without axis",
379 } else if(nBin<=0) {
380 Fatal("AddAxis","number of bins %d is not positive",
381 nBin);
382 } else {
383 TVectorD *bins=new TVectorD(nBin+1);
384 r=kTRUE;
385 for(Int_t i=0;i<=nBin;i++) {
386 (*bins)(i)=binBorders[i];
387 if(!TMath::Finite((*bins)(i))) {
388 Fatal("AddAxis","bin border %d is not finite",i);
389 r=kFALSE;
390 } else if((i>0)&&((*bins)(i)<=(*bins)(i-1))) {
391 Fatal("AddAxis","bins not in order x[%d]=%f <= %f=x[%d]",
392 i,(*bins)(i),(*bins)(i-1),i-1);
393 r=kFALSE;
394 }
395 }
396 if(r) {
398 Int_t bitMask=1<<axis;
399 Int_t nBinUO=nBin;
400 if(hasUnderflow) {
401 fHasUnderflow |= bitMask;
402 nBinUO++;
403 } else {
404 fHasUnderflow &= ~bitMask;
405 }
406 if(hasOverflow) {
407 fHasOverflow |= bitMask;
408 nBinUO++;
409 } else {
410 fHasOverflow &= ~bitMask;
411 }
412 fAxisList->AddLast(bins);
415 fDistributionSize *= nBinUO;
417 }
418 }
419 return r;
420}
421
422////////////////////////////////////////////////////////////////////////
423/// print some information about this binning tree
424///
425/// \param[out] out stream to write to
426/// \param[in] indent initial indentation (sub-trees have indent+1)
427/// \param[in] debug if debug&gt;0 print more information
428
429void TUnfoldBinning::PrintStream(ostream &out,Int_t indent,int debug)
430 const {
431 for(Int_t i=0;i<indent;i++) out<<" ";
432 out<<"TUnfoldBinning \""<<GetName()<<"\" has ";
433 Int_t nBin=GetEndBin()-GetStartBin();
434 if(nBin==1) {
435 out<<"1 bin";
436 } else {
437 out<<nBin<<" bins";
438 }
439 out<<" ["
440 <<GetStartBin()<<","<<GetEndBin()<<"] nTH1x="
442 <<"\n";
444 for(Int_t i=0;i<indent;i++) out<<" ";
445 out<<" distribution: "<<GetDistributionNumberOfBins()<<" bins\n";
447 /* for(Int_t i=nullptr;i<indent;i++) out<<" ";
448 out<<" axes:\n"; */
449 for(Int_t axis=0;axis<GetDistributionDimension();axis++) {
450 for(Int_t i=0;i<indent;i++) out<<" ";
451 out<<" \""
453 <<"\" nbin="<<GetDistributionBinning(axis)->GetNrows()-1;
454 if(HasUnderflow(axis)) out<<" plus underflow";
455 if(HasOverflow(axis)) out<<" plus overflow";
456 out<<"\n";
457 }
458 } else {
459 for(Int_t i=0;i<indent;i++) out<<" ";
460 out<<" no axis\n";
461 for(Int_t i=0;i<indent;i++) out<<" ";
462 out<<" names: ";
463 for(Int_t ibin=0;(ibin<GetDistributionNumberOfBins())&&
464 (ibin<fAxisLabelList->GetEntriesFast());ibin++) {
465 if(ibin) out<<";";
466 if(GetDistributionAxisLabel(ibin)) {
467 out<<GetDistributionAxisLabel(ibin);
468 }
469 }
470 out<<"\n";
471 }
472 if(debug>0) {
473 // print all bins with full name, size, status, user factor
474 for(int iBin=GetStartBin();iBin<GetEndBin();iBin++) {
475 for(Int_t i=0;i<indent;i++) out<<" ";
476 out<<GetBinName(iBin)
477 <<" size="<<GetBinSize(iBin)
478 <<" factor="<<GetBinFactor(iBin);
479 out<<"\n";
480 }
481 }
482 }
484 if(child) {
485 while(child) {
486 child->PrintStream(out,indent+1,debug);
487 child=child->GetNextNode();
488 }
489 }
490}
491
492////////////////////////////////////////////////////////////////////////
493/// set normalisation factors which are used in calls to GetBinFactor()
494///
495/// \param[in] normalisation normalisation factor
496/// \param[in] binfactor object which defines a factor for each bin
497///
498/// IN the presnet implementation, <b>binfactor</b> can be a TF1 or a
499/// TVectorD. The TF1 is evaluated a the bin centes of the
500/// relevant axes. The TVectorD is indexed by the global bin number
501/// minus the start bin number of this node.
503(Double_t normalisation,TObject *binfactor) {
504 fBinFactorConstant=normalisation;
506 if(!dynamic_cast<TF1 *>(fBinFactorFunction))
507 delete fBinFactorFunction;
508 }
509 fBinFactorFunction=binfactor;
510}
511
512////////////////////////////////////////////////////////////////////////
513/// set normalisation factor and function which are used in calls to GetBinFactor()
514///
515/// \param[in] normalisation normalisation factor
516/// \param[in] userFunc function evaluated at the (multi-dimensional)
517/// bin centres
519(Double_t normalisation,TF1 *userFunc) {
520 SetBinFactor(normalisation,userFunc);
521}
522
523/********************* Navigation **********************/
524
525////////////////////////////////////////////////////////////////////////
526/// traverse the tree and return the first node which matches the given name
527///
528/// \param[in] name the identifier of the node to find (zero matches
529/// the first node)
530///
531/// returns the node found or zero
533{
534 TUnfoldBinning const *r=nullptr;
535 if((!name)||(!TString(GetName()).CompareTo(name))) {
536 r=this;
537 }
538 for(TUnfoldBinning const *child=GetChildNode();
539 (!r) && child;child=child->GetNextNode()) {
540 r=child->FindNode(name);
541 }
542 return r;
543}
544
545////////////////////////////////////////////////////////////////////////
547{
548 // return root node
549 TUnfoldBinning *node=this;
550 while(node->GetParentNode()) node=node->parentNode;
551 return node;
552}
553
554////////////////////////////////////////////////////////////////////////
556{
557 // return root node
558 TUnfoldBinning const *node=this;
559 while(node->GetParentNode()) node=node->GetParentNode();
560 return node;
561}
562
563/********************* Create THxx histograms **********/
564
565////////////////////////////////////////////////////////////////////////
566/// construct a title
567///
568/// \param[in] histogramName distribution name
569/// \param[in] histogramTitle default title
570/// \param[in] axisList array indicating wqhich axis of this node is
571/// mapped to which histogram axis
572///
573/// if histogramTitle!=nullptr thsi title is used. Otherwise, the title is
574/// composed as:
575/// histogramName;axisname[axisList[0]];axisname[axisList[1]];...
577(const char *histogramName,const char *histogramTitle,Int_t const *axisList)
578 const
579{
580 TString r;
581 if(histogramTitle) {
582 r=histogramTitle;
583 } else {
584 r=histogramName;
585 Int_t iEnd;
586 for(iEnd=2;iEnd>0;iEnd--) {
587 if(axisList[iEnd]>=0) break;
588 }
589 for(Int_t i=0;i<=iEnd;i++) {
590 r += ";";
591 if(axisList[i]<0) {
592 r += GetName();
593 } else {
594 r += GetNonemptyNode()->GetDistributionAxisLabel(axisList[i]);
595 }
596 }
597 }
598 return r;
599}
600
601////////////////////////////////////////////////////////////////////////
602/// construct a histogram title for a 2D histogram with different
603/// binning schemes on x and y axis
604///
605/// \param[in] histogramName distribution name
606/// \param[in] histogramTitle default title
607/// \param[in] xAxis indicates which x-axis name to use
608/// \param[in] yAxisBinning binning scheme for y-axis
609/// \param[in] yAxis indicates which y-axis name to use
610///
612(const char *histogramName,const char *histogramTitle,
613 Int_t xAxis,const TUnfoldBinning *yAxisBinning,Int_t yAxis) const
614{
615 // build a title
616 // input:
617 // histogramTitle : if this is non-zero, use that title
618 // otherwise:
619 // title=histogramName;x;y
620 // xAxis : -1 no title for this axis
621 // >=nullptr use name of the corresponding axis
622 TString r;
623 if(histogramTitle) {
624 r=histogramTitle;
625 } else {
626 r=histogramName;
627 r += ";";
628 if(xAxis==-1) {
629 r += GetName();
630 } else if(xAxis>=0) {
632 }
633 r+= ";";
634 if(yAxis==-1) {
635 r += yAxisBinning->GetName();
636 } else if(yAxis>=0) {
637 r += yAxisBinning->GetNonemptyNode()->GetDistributionAxisLabel(yAxis);
638 }
639
640 }
641 return r;
642}
643
644////////////////////////////////////////////////////////////////////////
645/// return the number of histogram bins required when storing
646/// this binning in a one-dimensional histogram
647///
648/// \param[in] originalAxisBinning if true, try to have the histogram
649/// axis reflect precisely the relevant axis of the binnnig scheme
650/// \param[in] axisSteering steering to integrate over axis and/or
651/// skip underflow and overflow bins
652///
653/// returns the number of bins of the TH1, where the underflow/overflow
654/// are not used, unless the distribution has only one axis and
655/// originalAxisBinning=true)
656/// <br/>
657/// axisSteering is a string as follows:
658/// "axis[options];axis[options];..."
659/// where: axis = name or * is an identifier of an axis (* matches all)
660/// and: options is any combination of the letters C,U,O (other
661/// letters are ignored).
662/// <br>
663/// The letter C means that the corresponding axis is collapsed into
664/// one bin, i.e. one dimension is removed from the counting.
665/// The letters U,O remove for the matching axis the underflow.overflow
666/// bins from the counting
668(Bool_t originalAxisBinning,const char *axisSteering) const
669{
670 Int_t axisBins[3],axisList[3];
671 GetTHxxBinning(originalAxisBinning ? 1 : 0,axisBins,axisList,
672 axisSteering);
673 return axisBins[0];
674}
675
676////////////////////////////////////////////////////////////////////////
677/// create a THxx histogram capable to hold the bins of this binning
678/// node and its children
679///
680/// \param[in] histogramName name of the histogram which is created
681/// \param[in] originalAxisBinning if true, try to preserve the axis binning
682/// \param[out] (default=nullptr) binMap mapping of global bins to histogram bins.
683/// if(binMap==nullptr), no binMap is created
684/// \param[in] (default=nullptr) histogramTitle title o fthe histogram. If zero, a title
685/// is selected automatically
686/// \param[in] (default=nullptr) axisSteering steer the handling of underflow/overflow
687/// and projections
688///
689/// returns a new histogram (TH1D, TH2D or TH3D)
690/// <br>
691/// if the parameter <b>originalAxisBinning</b> parameter is true, the
692/// resulting histogram has bin widths and histogram dimension (TH1D,
693/// TH2D, TH3D) in parallel to this binning node, if possible.
694/// <br/>
695/// The <b>binMap</b> is an array which translates global bin numbers to bin
696/// numbers in the histogram returned by this method. The global bin
697/// numbers correspond to the bin numbers in a histogram created by
698/// calling GetRootNode()->CreateHistogram(name,false,0,0,0)
699/// <br/>
700/// The <b>axisSteering</b> is a string to steer whether underflow and
701/// overflow bins are included in the bin map. Furthermore, it is
702/// possible to "collapse" axes, such that their content is summed and
703/// the axis does not show up in the created histogram.
704/// <br/>
705/// The string looks like this: "axis[options];axis[options];..." where
706/// <ul>
707/// <li>axis is the name of an axis or equal to *, the latter matches
708/// all axes</li>
709/// <li>options is a combination of characters chosen from
710/// OUC0123456789
711/// <ul>
712/// <li>if O is included, the overflow bin of that axis is discarded</li>
713/// <li>if U is included, the underflow bin of that axis is discarded</li>
714/// <li>if C is included, the bins on that axes are collapsed,
715/// i.e. the corresponding histogram axis is not present in the output.
716/// The corrsponding bin contents are added up
717/// (projected onto the remaining axes). Using the characters O and U
718/// one can decide to exclude underflow or overflow from the
719/// projection. Using a selection of the characters 0123456789 one can
720/// restrict the sum further to only include the corresponding
721/// bins. In this counting, the first non-underflow bin corresponds to
722/// the character 0. This obviously only works for up to ten
723/// bins.</li>
724/// </ul>
725/// </li>
726/// </ul>
728(const char *histogramName,Bool_t originalAxisBinning,Int_t **binMap,
729 const char *histogramTitle,const char *axisSteering) const
730{
731 Int_t nBin[3],axisList[3];
732 Int_t nDim=GetTHxxBinning(originalAxisBinning ? 3 : 0,nBin,axisList,
733 axisSteering);
734 const TUnfoldBinning *neNode=GetNonemptyNode();
735 TString title=BuildHistogramTitle(histogramName,histogramTitle,axisList);
736 TH1 *r=nullptr;
737 if(nDim>0) {
738 const TVectorD *axisBinsX=
739 neNode->GetDistributionBinning(axisList[0]);
740 if(nDim>1) {
741 const TVectorD *axisBinsY=
742 neNode->GetDistributionBinning(axisList[1]);
743 if(nDim>2) {
744 const TVectorD *axisBinsZ=
745 neNode->GetDistributionBinning(axisList[2]);
746 r=new TH3D(histogramName,title,
747 nBin[0],axisBinsX->GetMatrixArray(),
748 nBin[1],axisBinsY->GetMatrixArray(),
749 nBin[2],axisBinsZ->GetMatrixArray());
750 } else {
751 r=new TH2D(histogramName,title,
752 nBin[0],axisBinsX->GetMatrixArray(),
753 nBin[1],axisBinsY->GetMatrixArray());
754 }
755 } else {
756 r=new TH1D(histogramName,title,nBin[0],axisBinsX->GetMatrixArray());
757 }
758 } else {
759 if(originalAxisBinning) {
760 Warning("CreateHistogram",
761 "Original binning can not be represented as THxx");
762 }
763 r=new TH1D(histogramName,title,nBin[0],0.5,nBin[0]+0.5);
764 nDim=0;
765 }
766 if(binMap) {
767 *binMap=CreateBinMap(r,nDim,axisList,axisSteering);
768 }
769 return r;
770}
771
772////////////////////////////////////////////////////////////////////////
773/// create a TH2D histogram capable to hold a covariance matrix
774///
775///
776/// \param[in] histogramName name of the histogram which is created
777/// \param[in] originalAxisBinning if true, try to preserve the axis binning
778/// \param[out] (default=nullptr) binMap mapping of global bins to histogram bins.
779/// if(binMap==nullptr), no binMap is created
780/// \param[in] (default=nullptr) histogramTitle title o fthe histogram. If zero, a title
781/// is selected automatically
782/// \param[in] (default=nullptr) axisSteering steer the handling of underflow/overflow
783/// and projections
784///
785/// returns a new TH2D. The options are described in greater detail
786/// with the CreateHistogram() method.
787
789(const char *histogramName,Bool_t originalAxisBinning,Int_t **binMap,
790 const char *histogramTitle,const char *axisSteering) const
791{
792 Int_t nBin[3],axisList[3];
793 Int_t nDim=GetTHxxBinning(originalAxisBinning ? 1 : 0,nBin,axisList,
794 axisSteering);
795 TString title=BuildHistogramTitle(histogramName,histogramTitle,axisList);
796 TH2D *r=nullptr;
797 if(nDim==1) {
798 const TVectorD *axisBinsX=(TVectorD const *)
799 GetNonemptyNode()->fAxisList->At(axisList[0]);
800 r=new TH2D(histogramName,title,nBin[0],axisBinsX->GetMatrixArray(),
801 nBin[0],axisBinsX->GetMatrixArray());
802 } else {
803 if(originalAxisBinning) {
804 Info("CreateErrorMatrixHistogram",
805 "Original binning can not be represented on one axis");
806 }
807 r=new TH2D(histogramName,title,nBin[0],0.5,nBin[0]+0.5,
808 nBin[0],0.5,nBin[0]+0.5);
809 nDim=0;
810 }
811 if(binMap) {
812 *binMap=CreateBinMap(r,nDim,axisList,axisSteering);
813 }
814 return r;
815}
816
817////////////////////////////////////////////////////////////////////////
818/// create a TH2D histogram capable to hold the bins of the two
819/// input binning schemes on the x and y axes, respectively
820///
821/// \paran[in] xAxis binning scheme for the x axis
822/// \param[in] yAxis binning scheme for the y axis
823/// \param[in] histogramName name of the histogram which is created
824/// \param[in] originalXAxisBinning preserve x-axis bin widths if possible
825/// \param[in] originalXAxisBinning preserve y-axis bin widths if possible
826/// \param[in] histogramTitle if is non-zero, it is taken as histogram title
827/// otherwise, the title is created automatically
828///
829/// returns a new TH2D.
830
832(TUnfoldBinning const *xAxis,TUnfoldBinning const *yAxis,
833 char const *histogramName,Bool_t originalXAxisBinning,
834 Bool_t originalYAxisBinning,char const *histogramTitle)
835{
836 Int_t nBinX[3],axisListX[3];
837 Int_t nDimX=
838 xAxis->GetTHxxBinning(originalXAxisBinning ? 1 : 0,nBinX,axisListX,nullptr);
839 const TUnfoldBinning *neNodeX=xAxis->GetNonemptyNode();
840 Int_t nBinY[3],axisListY[3];
841 Int_t nDimY=
842 yAxis->GetTHxxBinning(originalYAxisBinning ? 1 : 0,nBinY,axisListY,nullptr);
843 const TUnfoldBinning *neNodeY=yAxis->GetNonemptyNode();
844 TString title=xAxis->BuildHistogramTitle2D
845 (histogramName,histogramTitle,axisListX[0],yAxis,axisListY[0]);
846 if(nDimX==1) {
847 const TVectorD *axisBinsX=(TVectorD const *)
848 neNodeX->fAxisList->At(axisListX[0]);
849 if(nDimY==1) {
850 const TVectorD *axisBinsY=(TVectorD const *)
851 neNodeY->fAxisList->At(axisListY[0]);
852 return new TH2D(histogramName,title,
853 nBinX[0],axisBinsX->GetMatrixArray(),
854 nBinY[0],axisBinsY->GetMatrixArray());
855 } else {
856 return new TH2D(histogramName,title,
857 nBinX[0],axisBinsX->GetMatrixArray(),
858 nBinY[0],0.5,0.5+nBinY[0]);
859 }
860 } else {
861 if(nDimY==1) {
862 const TVectorD *axisBinsY=(TVectorD const *)
863 neNodeY->fAxisList->At(axisListY[0]);
864 return new TH2D(histogramName,title,
865 nBinX[0],0.5,0.5+nBinX[0],
866 nBinY[0],axisBinsY->GetMatrixArray());
867 } else {
868 return new TH2D(histogramName,title,
869 nBinX[0],0.5,0.5+nBinX[0],
870 nBinY[0],0.5,0.5+nBinY[0]);
871 }
872 }
873}
874
875////////////////////////////////////////////////////////////////////////
876/// calculate properties of a THxx histogram to store this binning
877///
878/// \param[in] maxDim maximum dimension of the THxx (0 or 1..3)
879/// maxDim==nullptr is used to indicate that the histogram should be
880/// dimensional with all bins mapped on one axis,
881/// bin centers equal to bin numbers
882/// \param[in] axisSteering see method CreateHistogram()
883/// \param[out] axisBins[3] number of bins on the THxx axes
884/// \param[out] axisList[3] TUnfoldBinning axis number corresponding
885/// to the THxx axis
886///
887/// returns 1-3 dimension of THxx or 0 for 1-dim THxx with equidistant bins
889(Int_t maxDim,Int_t *axisBins,Int_t *axisList,
890 const char *axisSteering) const
891{
892 for(Int_t i=0;i<3;i++) {
893 axisBins[i]=0;
894 axisList[i]=-1;
895 }
896 const TUnfoldBinning *theNode=GetNonemptyNode();
897 if(theNode) {
899 (maxDim,axisBins,axisList,axisSteering);
900 /* cout<<"GetTHxxBinning single node "<<r
901 <<" "<<axisBins[0]
902 <<" "<<axisBins[1]
903 <<" "<<axisBins[2]
904 <<"\n"; */
905 return r;
906 } else {
907 axisBins[0]=GetTHxxBinsRecursive(axisSteering);
908 //cout<<"GetTHxxBinning recursive "<<axisBins[0]<<"\n";
909 return 0;
910 }
911}
912
913////////////////////////////////////////////////////////////////////////
914/// find a node which has non-empty distributions
915/// if there is none or if there are many, return zero
917{
918 int count=0;
919 const TUnfoldBinning *r=GetNonemptyNode_r(count);
920 if(count!=1) r=nullptr;
921 return r;
922}
923
925{
926 const TUnfoldBinning *r=nullptr;
928 r=this;
929 count++;
930 }
932 child=child->GetNextNode()) {
933 const TUnfoldBinning *c=child->GetNonemptyNode_r(count);
934 if(!r) r=c;
935 }
936 return r;
937}
938
939////////////////////////////////////////////////////////////////////////
940/// get the properties of a histogram capable to hold the distribution
941/// attached to this node
942///
943/// \param[in] maxDim maximum dimension of the THxx (0 or 1..3)
944/// maxDim==nullptr is used to indicate that the histogram should
945/// 1-dimensional with all bins mapped on one axis
946/// \param[out] axisBins[3] number of bins on the THxx axes
947/// \param[out] axisList[3] TUnfoldBinning axis numbers
948/// corresponding to the THxx axis
949/// \param[in] axisSteering see method CreateHistogram()
950/// and projection
951///
952/// returns 1-3 dimension of THxx or use 1-dim THxx, binning structure
953/// is not preserved
955(Int_t maxDim,Int_t *axisBins,Int_t *axisList,const char *axisSteering) const
956{
957 // decode axisSteering
958 // isOptionGiven[0] ('C'): bit vector which axes to collapse
959 // isOptionGiven[1] ('U'): bit vector to discard underflow bins
960 // isOptionGiven[2] ('O'): bit vector to discard overflow bins
961 Int_t isOptionGiven[3];
962 DecodeAxisSteering(axisSteering,"CUO",isOptionGiven);
963 // count number of axes after projecting
964 Int_t numDimension=GetDistributionDimension();
965 Int_t r=0;
966 for(Int_t i=0;i<numDimension;i++) {
967 if(isOptionGiven[0] & (1<<i)) continue;
968 r++;
969 }
970 if((r>0)&&(r<=maxDim)) {
971 // 0<r<=maxDim
972 //
973 // -> preserve the original binning
974 // axisList[] and axisBins[] are overwritten
975 r=0;
976 for(Int_t i=0;i<numDimension;i++) {
977 if(isOptionGiven[0] & (1<<i)) continue;
978 axisList[r]=i;
979 axisBins[r]=GetDistributionBinning(i)->GetNrows()-1;
980 r++;
981 }
982 } else {
983 // map everything on one axis
984 // axisBins[0] is the number of bins
986 axisBins[0] = GetDistributionNumberOfBins();
987 } else {
988 Int_t nBin=1;
989 for(Int_t i=0;i<numDimension;i++) {
990 Int_t mask=(1<<i);
991 if(isOptionGiven[0] & mask) continue;
993 if((fHasUnderflow & mask)&& !(isOptionGiven[1] & mask)) nBinI++;
994 if((fHasOverflow & mask)&& !(isOptionGiven[2] & mask)) nBinI++;
995 nBin *= nBinI;
996 }
997 axisBins[0] = nBin;
998 }
999 r=0;
1000 }
1001 return r;
1002}
1003
1004////////////////////////////////////////////////////////////////////////
1005/// calculate number of bins required to store this binning with teh
1006/// given axisSteering
1007///
1008/// \param[in] axisSteering see method CreateHistogram()
1009///
1010/// returns the number of bins
1011Int_t TUnfoldBinning::GetTHxxBinsRecursive(const char *axisSteering) const
1012{
1013
1014 Int_t r=0;
1016 child=child->GetNextNode()) {
1017 r +=child->GetTHxxBinsRecursive(axisSteering);
1018 }
1019 // here: process distribution of this node
1020 Int_t axisBins[3],axisList[3];
1021 GetTHxxBinningSingleNode(0,axisBins,axisList,axisSteering);
1022 r += axisBins[0];
1023 // cout<<"GetTHxxBinsRecursive returns "<<r<<"\n";
1024 return r;
1025}
1026
1027////////////////////////////////////////////////////////////////////////
1028/// create an empty bin map, useful together with the getter methods of
1029/// class TUnfold and TUnfoldSys
1030///
1031/// returns: a new Int array of the proper size, all eleemnts set to -1
1033 // create empty bin map which can be manipulated by
1034 // MapGlobalBin()
1035 Int_t nMax=GetRootNode()->GetEndBin()+1;
1036 Int_t *r=new Int_t[nMax];
1037 for(Int_t i=0;i<nMax;i++) {
1038 r[i]=-1;
1039 }
1040 return r;
1041}
1042
1043////////////////////////////////////////////////////////////////////////
1044/// set one entry in a bin map
1045///
1046/// \param[out] binMap to be used with TUnfoldSys::GetOutput() etc
1047/// \param[in] source bin, global bin number in this binning scheme
1048/// \param[in] destination bin in the output histogram
1049///
1050
1052(Int_t *binMap,Int_t globalBin,Int_t destBin) const {
1053 Int_t nMax=GetRootNode()->GetEndBin()+1;
1054 if((globalBin<0)||(globalBin>=nMax)) {
1055 Error("SetBinMapEntry","global bin number %d outside range (max=%d)",
1056 globalBin,nMax);
1057 } else {
1058 binMap[globalBin]=destBin;
1059 }
1060}
1061
1062////////////////////////////////////////////////////////////////////////
1063/// map all global bins referenced by this node to the one-dimensional
1064/// histogram destHist, starting with bin firstBinX
1065///
1066/// \param[out] binMap to be used with TUnfoldSys::GetOutput() etc
1067/// \param[in] axisSteering steering for underflow/overflow/projections
1068/// \param[in] firstBinX first bin of destination histogram to be filled
1069///
1070/// returns: highest bin number in destibation histogram plus 1
1071/// <br/>The parameter <b>axisSteering</b> is explained with the
1072/// method CreateHistogram()
1074(Int_t *binMap,const char *axisSteering,Int_t firstBinX) const {
1075 Int_t r=firstBinX;
1076 Int_t axisBins[3],axisList[3];
1077 Int_t nDim=GetTHxxBinningSingleNode(3,axisBins,axisList,axisSteering);
1078 if((nDim==1)|| !GetDistributionDimension()) {
1079 r+=FillBinMapSingleNode(nullptr,r,0,nullptr,axisSteering,binMap);
1080 } else {
1081 Error("FillBinMap1D","distribution %s with steering=%s is not 1D",
1082 (char *)GetName(),axisSteering);
1083 }
1085 child=child->GetNextNode()) {
1086 r =child->FillBinMap1D(binMap,axisSteering,r);
1087 }
1088 return r;
1089}
1090
1091////////////////////////////////////////////////////////////////////////
1092/// create mapping from global bin number to a histogram for this node
1093///
1094/// \param[in] hist destination histogram
1095/// \param[in] nDim target dimension
1096/// \param[in] axisList map axes in the binning scheme to histogram axes
1097/// \param[in] axisSteering steering for underflow/overflow/projections
1098///
1099/// The <b>axisSteering</b> is explained with the method CreateHistogram()
1101(const TH1 *hist,Int_t nDim,const Int_t *axisList,const char *axisSteering)
1102 const
1103{
1104 // create mapping from global bin number to a histogram for this node
1105 // global bins are the bins of the root node binning scheme
1106 // when projecting them on a TH1 histogram "hRootNode" without special
1107 // axis steering and without attempting to preserve the axis binnings
1108 //
1109 // The bin map is an array of size hRootNode->GetNbinsX()+2
1110 // For each bin of the "hRootNode" histogram it holds the target bin in
1111 // "hist" or the number -1 if the corresponding "hRootNode" bin is not
1112 // represented in "hist"
1113 //
1114 // input
1115 // hist : the histogram (to calculate root bin numbers)
1116 // nDim : target dimension of the TUnfoldBinning
1117 // if(nDim==nullptr) all bins are mapped linearly
1118 // axisSteering:
1119 // "pattern1;pattern2;...;patternN"
1120 // patternI = axis[mode]
1121 // axis = name or *
1122 // mode = C|U|O
1123 // C: collapse axis into one bin
1124 // U: discard underflow bin
1125 // O: discard overflow bin
1126 //
1127 // input used only if nDim>0:
1128 // axisList : for each THxx axis give the TUnfoldBinning axis number
1129 //
1130 // return value:
1131 // an new array which holds the bin mapping
1132 // r[0] : to which THxx bin to map global bin number 0
1133 // r[1] : to which THxx bin to map global bin number 1
1134 // ...
1135 // r[nmax]
1136 // where nmax=GetRootNode()->GetEndBin()+1
1138 Int_t startBin=GetRootNode()->GetStartBin();
1139 if(nDim>0) {
1140 const TUnfoldBinning *nonemptyNode=GetNonemptyNode();
1141 if(nonemptyNode) {
1142 nonemptyNode->
1143 FillBinMapSingleNode(hist,startBin,nDim,axisList,axisSteering,r);
1144 } else {
1145 Fatal("CreateBinMap","called with nDim=%d but GetNonemptyNode()=nullptr",
1146 nDim);
1147 }
1148 } else {
1149 FillBinMapRecursive(startBin,axisSteering,r);
1150 }
1151 return r;
1152}
1153
1154////////////////////////////////////////////////////////////////////////
1155/// recursively fill bin map
1156///
1157/// \param[in] startBin first histogram bin
1158/// \param[in] axisSteering see CreateHistogram() method
1159/// \param[out] binMap the bin mapping which is to be filled
1160///
1161/// the positions
1162/// binMap[GetStartBin()]...binMap[GetEndBin()-1]
1163/// are filled
1165(Int_t startBin,const char *axisSteering,Int_t *binMap) const
1166{
1167 Int_t nbin=0;
1168 nbin = FillBinMapSingleNode(nullptr,startBin,0,nullptr,axisSteering,binMap);
1170 child=child->GetNextNode()) {
1171 nbin += child->FillBinMapRecursive(startBin+nbin,axisSteering,binMap);
1172 }
1173 return nbin;
1174}
1175
1176////////////////////////////////////////////////////////////////////////
1177/// fill bin map for a single node
1178///
1179/// \param[in] hist the histogram representing this node (used if nDim>0)
1180/// \param[in] startBin start bin in the bin map
1181/// \param[in] nDim number of dimensions to resolve
1182/// \param[in] axisList[3] TUnfoldBinning axis numbers corresponding
1183/// to the axes of <b>hist</b>
1184/// \param[in] axisSteering see documentation of CreateHistogram()
1185/// \param[out] binMap the bin map to fill
1186///
1187/// returns the number of bins mapped.
1188/// <br>
1189/// The result depends on the parameter <b>nDim</b> as follows
1190/// <ul>
1191/// <li> nDim==nullptr: bins are mapped in linear order, ignore hist and
1192/// axisList</li>
1193/// <li> nDim==hist->GetDimension():
1194/// bins are mapped to "hist" bin numbers
1195/// the corresponding TUnfoldBinning axes are taken from
1196/// axisList[]</li>
1197/// <li>nDim=1 and hist->GetDimension()>1:
1198/// bins are mapped to the x-axis of "hist"
1199/// the corresponding TUnfoldBinning axis is taken from
1200/// axisList[0]</li>
1201/// </ul>
1203(const TH1 *hist,Int_t startBin,Int_t nDim,const Int_t *axisList,
1204 const char *axisSteering,Int_t *binMap) const
1205{
1206 // first, decode axisSteering
1207 // isOptionGiven[0] ('C'): bit vector which axes to collapse
1208 // isOptionGiven[1] ('U'): bit vector to discard underflow bins
1209 // isOptionGiven[2] ('O'): bit vector to discard overflow bins
1210 Int_t isOptionGiven[3+10];
1211 DecodeAxisSteering(axisSteering,"CUO0123456789",isOptionGiven);
1212 Int_t haveSelectedBin=0;
1213 for(Int_t i=3;i<3+10;i++) {
1214 haveSelectedBin |= isOptionGiven[i];
1215 }
1216
1217 Int_t axisBins[MAXDIM];
1218 Int_t dimension=GetDistributionDimension();
1219 Int_t axisNbin[MAXDIM];
1220 for(Int_t i=0;i<dimension;i++) {
1221 const TVectorD *binning=GetDistributionBinning(i);
1222 axisNbin[i]=binning->GetNrows()-1;
1223 };
1224 for(Int_t i=0;i<GetDistributionNumberOfBins();i++) {
1225 Int_t globalBin=GetStartBin()+i;
1226 const TUnfoldBinning *dest=ToAxisBins(globalBin,axisBins);
1227 if(dest!=this) {
1228 if(!dest) {
1229 Fatal("FillBinMapSingleNode",
1230 "bin %d outside binning scheme",
1231 globalBin);
1232 } else {
1233 Fatal("FillBinMapSingleNode",
1234 "bin %d located in %s %d-%d rather than %s %d=%d",
1235 i,(const char *)dest->GetName(),
1236 dest->GetStartBin(),dest->GetEndBin(),
1237 (const char *)GetName(),GetStartBin(),GetEndBin());
1238 }
1239 }
1240 // check whether this bin has to be skipped
1241 Bool_t skip=kFALSE;
1242 for(Int_t axis=0;axis<dimension;axis++) {
1243 Int_t mask=(1<<axis);
1244 // underflow/overflow excluded by steering
1245 if(((axisBins[axis]<0)&&(isOptionGiven[1] & mask))||
1246 ((axisBins[axis]>=axisNbin[axis])&&(isOptionGiven[2] & mask)))
1247 skip=kTRUE;
1248 // only certain bins selected by steering
1249 if((axisBins[axis]>=0)&&(axisBins[axis]<axisNbin[axis])&&
1250 (haveSelectedBin & mask)) {
1251 if(!(isOptionGiven[3+axisBins[axis]] & mask)) skip=kTRUE;
1252 }
1253 }
1254 if(skip) {
1255 continue;
1256 }
1257
1258 if(nDim>0) {
1259 // get bin number from THxx function(s)
1260 if(nDim==hist->GetDimension()) {
1261 Int_t ibin[3];
1262 ibin[0]=ibin[1]=ibin[2]=0;
1263 for(Int_t hdim=0;hdim<nDim;hdim++) {
1264 Int_t axis=axisList[hdim];
1265 ibin[hdim]=axisBins[axis]+1;
1266 }
1267 binMap[globalBin]=hist->GetBin(ibin[0],ibin[1],ibin[2]);
1268 //cout<<globalBin<<"->"<<binMap[globalBin]
1269 // <<" ("<<ibin[0]<<","<<ibin[1]<<","<<ibin[2]<<")\n";
1270 } else if(nDim==1) {
1271 // histogram has more dimensions than the binning scheme
1272 // and the binning scheme has one axis only
1273 //
1274 // special case: error histogram is 2-d
1275 // create nor error if ndim==1 && hist->GetDimension()==2
1276 if((nDim!=1)||( hist->GetDimension()!=2)) {
1277 // -> use the first valid axis only
1278 Error("FillBinMapSingleNode","inconsistent dimensions %d %d",nDim,
1279 hist->GetDimension());
1280 }
1281 for(Int_t ii=0;ii<hist->GetDimension();ii++) {
1282 if(axisList[ii]>=0) {
1283 binMap[globalBin]=axisBins[axisList[ii]]+1;
1284 break;
1285 }
1286 }
1287 } else {
1288 Fatal("FillBinMapSingleNode","inconsistent dimensions %d %d",nDim,
1289 hist->GetDimension());
1290 }
1291 } else {
1292 // order all bins in sequence
1293 // calculation in parallel to ToGlobalBin()
1294 // but take care of
1295 // startBin,collapseAxis,discardeUnderflow,discardeOverflow
1296 if(dimension>0) {
1297 Int_t r=0;
1298 for(Int_t axis=dimension-1;axis>=0;axis--) {
1299 Int_t mask=(1<<axis);
1300 if(isOptionGiven[0] & mask) {
1301 // bins on this axis are integrated over
1302 continue;
1303 }
1304 Int_t iBin=axisBins[axis];
1305 Int_t nMax=axisNbin[axis];
1306 if((fHasUnderflow & ~isOptionGiven[1]) & mask) {
1307 nMax +=1;
1308 iBin +=1;
1309 }
1310 if((fHasOverflow & ~isOptionGiven[2]) & mask) {
1311 nMax += 1;
1312 }
1313 r = r*nMax +iBin;
1314 }
1315 binMap[globalBin] = startBin + r;
1316 } else {
1317 binMap[globalBin] = startBin + axisBins[0];
1318 }
1319 }
1320 }
1321 Int_t nbin;
1322 if(dimension>0) {
1323 nbin=1;
1324 for(Int_t axis=dimension-1;axis>=0;axis--) {
1325 Int_t mask=(1<<axis);
1326 if(isOptionGiven[0] & mask) {
1327 // bins on this axis are integrated over
1328 continue;
1329 }
1330 Int_t nMax=axisNbin[axis];
1331 if((fHasUnderflow & ~isOptionGiven[1]) & mask) {
1332 nMax +=1;
1333 }
1334 if((fHasOverflow & ~isOptionGiven[2]) & mask) {
1335 nMax += 1;
1336 }
1337 nbin = nbin*nMax;
1338 }
1339 } else {
1341 }
1342 return nbin;
1343}
1344
1345////////////////////////////////////////////////////////////////////////
1347(const char *histogramName,const TH1 *globalBins,
1348 const TH2 *globalBinsEmatrix,Bool_t originalAxisBinning,
1349 const char *axisSteering) const
1350{
1351 // extract a distribution from the given set of global bins
1352 // input:
1353 // histogramName : name of the histogram which is created
1354 // globalBins : histogram with all bins
1355 // globalBinsEmatrix : corresponding error matrix
1356 // if this pointer is zero, only diagonal errors
1357 // are considered
1358 // originalAxisBinning : extract histogram with proper binning
1359 // (if possible)
1360 // axisSteering
1361 // "pattern1;pattern2;...;patternN"
1362 // patternI = axis[mode]
1363 // axis = name or *
1364 // mode = C|U|O
1365 // C: collapse axis into one bin
1366 // U: discard underflow bin
1367 // O: discard overflow bin
1368 Int_t *binMap=nullptr;
1369 TH1 *r=CreateHistogram(histogramName,originalAxisBinning,&binMap,nullptr,
1370 axisSteering);
1371 if(!r) return nullptr;
1372 TUnfoldBinning const *root=GetRootNode();
1373 Int_t nMax=-1;
1374 for(Int_t iSrc=root->GetStartBin();iSrc<root->GetEndBin();iSrc++) {
1375 if(binMap[iSrc]>nMax) nMax=binMap[iSrc];
1376 }
1377 if(nMax<0) {
1378 delete r;
1379 r=nullptr;
1380 } else {
1381 TVectorD eSquared(nMax+1);
1382 for(Int_t iSrc=root->GetStartBin();iSrc<root->GetEndBin();iSrc++) {
1383 Int_t iDest=binMap[iSrc];
1384 if(iDest>=0) {
1385 Double_t c=r->GetBinContent(iDest);
1386 r->SetBinContent(iDest,c+globalBins->GetBinContent(iSrc));
1387 if(!globalBinsEmatrix) {
1388 eSquared(iDest)+=TMath::Power(globalBins->GetBinError(iSrc),2.);
1389 } else {
1390 for(Int_t jSrc=root->GetStartBin();jSrc<root->GetEndBin();
1391 jSrc++) {
1392 if(binMap[jSrc]==iDest) {
1393 eSquared(iDest) +=
1394 globalBinsEmatrix->GetBinContent(iSrc,jSrc);
1395 }
1396 }
1397 }
1398 }
1399 }
1400 for(Int_t i=0;i<=nMax;i++) {
1401 Double_t e2=eSquared(i);
1402 if(e2>0.0) {
1403 r->SetBinError(i,TMath::Sqrt(e2));
1404 }
1405 }
1406 }
1407 delete binMap;
1408 return r;
1409}
1410
1411/********************* Calculate global bin number ******/
1412
1413////////////////////////////////////////////////////////////////////////
1414/// locate a bin in a one-dimensional distribution
1415///
1416/// \param[in] x coordinate
1417///
1418/// returns the global bin number within the distribution attached to
1419/// this node. THe global bin number is valid for the root node of the
1420/// binning scheme
1422{
1423 if(GetDistributionDimension()!=1) {
1424 Fatal("GetBinNumber",
1425 "called with 1 argument for %d dimensional distribution",
1427 }
1428 return GetGlobalBinNumber(&x);
1429}
1430
1431////////////////////////////////////////////////////////////////////////
1432/// locate a bin in a two-dimensional distribution
1433///
1434/// \param[in] x coordinate on first axis
1435/// \param[in] y coordinate on second axis
1436///
1437/// returns the global bin number within the distribution attached to
1438/// this node. The global bin number is valid for the root node of the
1439/// binning scheme
1441{
1442 if(GetDistributionDimension()!=2) {
1443 Fatal("GetBinNumber",
1444 "called with 2 arguments for %d dimensional distribution",
1446 }
1447 Double_t xx[2];
1448 xx[0]=x;
1449 xx[1]=y;
1450 return GetGlobalBinNumber(xx);
1451}
1452
1453////////////////////////////////////////////////////////////////////////
1454/// locate a bin in a three-dimensional distribution
1455///
1456/// \param[in] x coordinate on first axis
1457/// \param[in] y coordinate on second axis
1458/// \param[in] z coordinate on third axis
1459///
1460/// returns the global bin number within the distribution attached to
1461/// this node. The global bin number is valid for the root node of the
1462/// binning scheme
1464(Double_t x,Double_t y,Double_t z) const
1465{
1466 // locate bin on a three-dimensional distribution
1467 // input
1468 // x,y,z: coordinates to locate
1469 if(GetDistributionDimension()!=3) {
1470 Fatal("GetBinNumber",
1471 "called with 3 arguments for %d dimensional distribution",
1473 }
1474 Double_t xx[3];
1475 xx[0]=x;
1476 xx[1]=y;
1477 xx[2]=z;
1478 return GetGlobalBinNumber(xx);
1479}
1480
1481////////////////////////////////////////////////////////////////////////
1482/// locate a bin in a four-dimensional distribution
1483///
1484/// \param[in] x0 coordinate on first axis
1485/// \param[in] x1 coordinate on second axis
1486/// \param[in] x2 coordinate on third axis
1487/// \param[in] x3 coordinate on fourth axis
1488///
1489/// returns the global bin number within the distribution attached to
1490/// this node. The global bin number is valid for the root node of the
1491/// binning scheme
1493(Double_t x0,Double_t x1,Double_t x2,Double_t x3) const
1494{
1495 // locate bin on a four-dimensional distribution
1496 // input
1497 // x0,x1,x2,x3: coordinates to locate
1498 if(GetDistributionDimension()!=4) {
1499 Fatal("GetBinNumber",
1500 "called with 4 arguments for %d dimensional distribution",
1502 }
1503 Double_t xx[4];
1504 xx[0]=x0;
1505 xx[1]=x1;
1506 xx[2]=x2;
1507 xx[3]=x3;
1508 return GetGlobalBinNumber(xx);
1509}
1510
1511////////////////////////////////////////////////////////////////////////
1512/// locate a bin in a five-dimensional distribution
1513///
1514/// \param[in] x0 coordinate on first axis
1515/// \param[in] x1 coordinate on second axis
1516/// \param[in] x2 coordinate on third axis
1517/// \param[in] x3 coordinate on fourth axis
1518/// \param[in] x4 coordinate on fifth axis
1519///
1520/// returns the global bin number within the distribution attached to
1521/// this node. The global bin number is valid for the root node of the
1522/// binning scheme
1525{
1526 // locate bin on a five-dimensional distribution
1527 // input
1528 // x0,x1,x2,x3,x4: coordinates to locate
1529 if(GetDistributionDimension()!=5) {
1530 Fatal("GetBinNumber",
1531 "called with 5 arguments for %d dimensional distribution",
1533 }
1534 Double_t xx[5];
1535 xx[0]=x0;
1536 xx[1]=x1;
1537 xx[2]=x2;
1538 xx[3]=x3;
1539 xx[4]=x4;
1540 return GetGlobalBinNumber(xx);
1541}
1542
1543////////////////////////////////////////////////////////////////////////
1544/// locate a bin in a six-dimensional distribution
1545///
1546/// \param[in] x0 coordinate on first axis
1547/// \param[in] x1 coordinate on second axis
1548/// \param[in] x2 coordinate on third axis
1549/// \param[in] x3 coordinate on fourth axis
1550/// \param[in] x4 coordinate on fifth axis
1551/// \param[in] x5 coordinate on sixth axis
1552///
1553/// returns the global bin number within the distribution attached to
1554/// this node. The global bin number is valid for the root node of the
1555/// binning scheme
1558{
1559 // locate bin on a five-dimensional distribution
1560 // input
1561 // x0,x1,x2,x3,x4,x5: coordinates to locate
1562 if(GetDistributionDimension()!=6) {
1563 Fatal("GetBinNumber",
1564 "called with 6 arguments for %d dimensional distribution",
1566 }
1567 Double_t xx[6];
1568 xx[0]=x0;
1569 xx[1]=x1;
1570 xx[2]=x2;
1571 xx[3]=x3;
1572 xx[4]=x4;
1573 xx[5]=x5;
1574 return GetGlobalBinNumber(xx);
1575}
1576
1577////////////////////////////////////////////////////////////////////////
1578/// locate a bin in an N-dimensional distribution
1579///
1580/// \param[in] x array of coordinates
1581/// \param[out] isBelow pointer to an integer (bit vector) to indicate
1582/// coordinates which do not fit in the binning scheme
1583/// \param[out] isAbove pointer to an integer (bit vector) to indicate
1584/// coordinates which do not fit in the binning scheme
1585///
1586/// returns the global bin number within the distribution attached to
1587/// this node. The global bin number is valid for the root node of the
1588/// binning scheme. If some coordinates do not fit, zero is returned.
1589/// The integers pointed to by isBelow and isAbove are set to zero.
1590/// However, if coordinate i is below
1591/// the lowest bin border and there is no underflow bin, the bin i is
1592/// set in (*isBelow). Overflows are handled in a similar manner with
1593/// (*isAbove).
1594/// <br>
1595/// If a coordinate is NaN, the result is undefined for TUnfold
1596/// Version<17.6. As of version 17.6, NaN is expected to end up in the
1597/// underflow or by setting the corresponding bit in (*isBelow).
1598
1600(const Double_t *x,Int_t *isBelow,Int_t *isAbove) const
1601{
1602 // locate bin on a n-dimensional distribution
1603 // input
1604 // x[]: coordinates to locate
1605 // output:
1606 // isBelow,isAbove: bit vectors,
1607 // indicating which cut on which axis failed
1609 Fatal("GetBinNumber",
1610 "no axes are defined for node %s",
1611 (char const *)GetName());
1612 }
1613 Int_t iAxisBins[MAXDIM];
1614 for(Int_t dim=0;dim<GetDistributionDimension();dim++) {
1615 TVectorD const *bins=(TVectorD const *) fAxisList->At(dim);
1616 Int_t i0=0;
1617 Int_t i1=bins->GetNrows()-1;
1618 Int_t iBin= 0;
1619 if(!(x[dim]>=(*bins)[i0])) {
1620 // underflow or NaN
1621 iBin += i0-1;
1622 } else if(!(x[dim]<(*bins)[i1])) {
1623 // overflow
1624 iBin += i1;
1625 } else {
1626 while(i1-i0>1) {
1627 Int_t i2=(i0+i1)/2;
1628 if(x[dim]<(*bins)[i2]) {
1629 i1=i2;
1630 } else {
1631 i0=i2;
1632 }
1633 }
1634 iBin += i0;
1635 }
1636 iAxisBins[dim]=iBin;
1637 }
1638 Int_t r=ToGlobalBin(iAxisBins,isBelow,isAbove);
1639 if(r<0) r=0;
1640 return r;
1641}
1642
1643/********************* access by global bin number ******/
1644
1645////////////////////////////////////////////////////////////////////////
1646/// get the name of a bin
1647///
1648/// \param[in] iBin global bin number
1649///
1650/// returns a string describing the bin
1652{
1653 // get the name of a bin in the given tree
1654 // iBin: bin number
1655 Int_t axisBins[MAXDIM];
1656 TString r=TString::Format("#%d",iBin);
1657 TUnfoldBinning const *distribution=ToAxisBins(iBin,axisBins);
1658 if(distribution) {
1659 r +=" (";
1660 r += distribution->GetName();
1661 Int_t dimension=distribution->GetDistributionDimension();
1662 if(dimension>0) {
1663 TString axisString;
1664 for(Int_t axis=0;axis<dimension;axis++) {
1665 TString thisAxisString=
1666 distribution->GetDistributionAxisLabel(axis);
1667 TVectorD const *bins=distribution->GetDistributionBinning(axis);
1668 Int_t i=axisBins[axis];
1669 if(i<0) thisAxisString += "[ufl]";
1670 else if(i>=bins->GetNrows()-1) thisAxisString += "[ofl]";
1671 else {
1672 thisAxisString +=
1673 TString::Format("[%.3g,%.3g]",(*bins)[i],(*bins)[i+1]);
1674 }
1675 axisString = ":"+thisAxisString+axisString;
1676 }
1677 r += axisString;
1678 } else {
1679 // extra bins
1680 Int_t i=axisBins[0];
1681 if((i>=0)&&(i<distribution->fAxisLabelList->GetEntriesFast())) {
1682 r += distribution->GetDistributionAxisLabel(i);
1683 } else {
1684 r += TString::Format(" %d",i);
1685 }
1686 }
1687 r +=")";
1688 }
1689 return r;
1690}
1691
1692////////////////////////////////////////////////////////////////////////
1693/// get N-dimensional bin size
1694///
1695/// \param[in] iBin global bin number
1697{
1698 // includeUO : include underflow/overflow bins or not
1699 Int_t axisBins[MAXDIM];
1700 TUnfoldBinning const *distribution=ToAxisBins(iBin,axisBins);
1701 Double_t r=0.0;
1702 if(distribution) {
1703 if(distribution->GetDistributionDimension()>0) r=1.0;
1704 for(Int_t axis=0;axis<distribution->GetDistributionDimension();axis++) {
1705 TVectorD const *bins=distribution->GetDistributionBinning(axis);
1706 Int_t pos=axisBins[axis];
1707 if(pos<0) {
1708 r *= distribution->GetDistributionUnderflowBinWidth(axis);
1709 } else if(pos>=bins->GetNrows()-1) {
1710 r *= distribution->GetDistributionOverflowBinWidth(axis);
1711 } else {
1712 r *= (*bins)(pos+1)-(*bins)(pos);
1713 }
1714 if(r<=0.) break;
1715 }
1716 }
1717 return r;
1718}
1719
1720////////////////////////////////////////////////////////////////////////
1721/// check whether there is only a global scaling factor for this node
1723 return fBinFactorFunction ? kFALSE : kTRUE;
1724}
1725
1726////////////////////////////////////////////////////////////////////////
1727/// return global scaling factor for this node
1729 return fBinFactorConstant;
1730}
1731
1732////////////////////////////////////////////////////////////////////////
1733/// return scaling factor for the given global bin number
1734///
1735/// \param[in] iBin global bin number
1736///
1737/// returns the scaling factor for this bin.
1738/// The scaling factors can be set using the method SetBinFactorFunction()
1740{
1741 // return user factor for a bin
1742 // iBin : global bin number
1743 Int_t axisBins[MAXDIM];
1744 TUnfoldBinning const *distribution=ToAxisBins(iBin,axisBins);
1745 Double_t r=distribution->fBinFactorConstant;
1746 if((r!=0.0) && distribution->fBinFactorFunction) {
1747 TF1 *function=dynamic_cast<TF1 *>(distribution->fBinFactorFunction);
1748 if(function) {
1749 Double_t x[MAXDIM];
1750 Int_t dimension=distribution->GetDistributionDimension();
1751 if(dimension>0) {
1752 for(Int_t axis=0;axis<dimension;axis++) {
1753 x[axis]=distribution->GetDistributionBinCenter
1754 (axis,axisBins[axis]);
1755 }
1756 r *= function->EvalPar(x,function->GetParameters());
1757 } else {
1758 x[0]=axisBins[0];
1759 r *= function->Eval(x[0]);
1760 }
1761 } else {
1762 TVectorD *vect=dynamic_cast<TVectorD *>
1763 (distribution->fBinFactorFunction);
1764 if(vect) {
1765 r=(*vect)[iBin-GetStartBin()];
1766 } else {
1767 Error("GetBinFactor",
1768 "internal error: user function is neiter TF1 or TVectorD");
1769 }
1770 }
1771 }
1772 return r;
1773}
1774
1775////////////////////////////////////////////////////////////////////////
1776/// get neighbour bins along the specified axis
1777///
1778/// \param[in] bin global bin number
1779/// \param[in] axis axis number of interest
1780/// \param[out] prev bin number of previous bin or -1 if not existing
1781/// \param[out] distPrev distance between bin centres
1782/// \param[out] next bin number of next bin or -1 if not existing
1783/// \param[out] distNext distance between bin centres
1784/// \param[in] isPeriodic (default=false) if true, the first bin is counted as neighbour of the last bin
1785///
1786/// return code
1787/// <ul>
1788/// <li>0 everything is fine</li>
1789/// <li>1,2,3 isPeriodic option was reset to false, because underflow/overflow
1790/// bins are present
1791/// <ul>
1792/// <li>+1 invalid isPeriodic option was specified with underflow bin</li>
1793/// <li>+2 invalid isPeriodic option was specified with overflow bin</li>
1794/// </ul>
1795/// </li>
1796/// </ul>
1797
1799(Int_t bin,Int_t axis,Int_t *prev,Double_t *distPrev,
1800 Int_t *next,Double_t *distNext,Bool_t isPeriodic) const
1801{
1802 Int_t axisBins[MAXDIM];
1803 TUnfoldBinning const *distribution=ToAxisBins(bin,axisBins);
1804 Int_t dimension=distribution->GetDistributionDimension();
1805 *prev=-1;
1806 *next=-1;
1807 *distPrev=0.;
1808 *distNext=0.;
1809 Int_t r=0;
1810 if((axis>=0)&&(axis<dimension)) {
1811 //TVectorD const *bins=distribution->GetDistributionBinning(axis);
1812 //Int_t nBin=bins->GetNrows()-1;
1813 Int_t nMax=GetDistributionBinning(axis)->GetNrows()-1;
1814 Int_t centerBin= axisBins[axis];
1815 axisBins[axis] =centerBin-1;
1816 if(isPeriodic) {
1817 if(HasUnderflow(axis)) {
1818 r +=1;
1819 } else if((axisBins[axis]<0)&&(nMax>=3)) {
1820 axisBins[axis]=nMax-1;
1821 }
1822 }
1823 *prev=ToGlobalBin(axisBins);
1824 if(*prev>=0) {
1825 *distPrev=distribution->GetDistributionBinCenter(axis,axisBins[axis])-
1826 distribution->GetDistributionBinCenter(axis,centerBin);
1827 }
1828 axisBins[axis] =centerBin+1;
1829 if(isPeriodic) {
1830 if(HasOverflow(axis)) {
1831 r +=2;
1832 } else if((axisBins[axis]==nMax)&&(nMax>=3)) {
1833 axisBins[axis]=0;
1834 }
1835 }
1836 *next=ToGlobalBin(axisBins);
1837 if(*next>=0) {
1838 *distNext=distribution->GetDistributionBinCenter(axis,axisBins[axis])-
1839 distribution->GetDistributionBinCenter(axis,centerBin);
1840 }
1841 }
1842 return r;
1843}
1844
1845////////////////////////////////////////////////////////////////////////
1846/// return bit maps indicating underflow and overflow status
1847///
1848/// \param[in] iBin global bin number
1849/// \param[out] uStatus bit map indicating whether the bin is underflow
1850/// \param[out] oStatus bit map indicating whether the bin is overflow
1852(Int_t iBin,Int_t *uStatus,Int_t *oStatus) const
1853{
1854 Int_t axisBins[MAXDIM];
1855 TUnfoldBinning const *distribution=ToAxisBins(iBin,axisBins);
1856 Int_t dimension=distribution->GetDistributionDimension();
1857 *uStatus=0;
1858 *oStatus=0;
1859 for(Int_t axis=0;axis<dimension;axis++) {
1860 TVectorD const *bins=distribution->GetDistributionBinning(axis);
1861 Int_t nBin=bins->GetNrows()-1;
1862 if(axisBins[axis]<0) *uStatus |= (1<<axis);
1863 if(axisBins[axis]>=nBin) *oStatus |= (1<<axis);
1864 }
1865}
1866
1867////////////////////////////////////////////////////////////////////////
1868/// check whether there are bins but no axis
1870{
1872}
1873
1874////////////////////////////////////////////////////////////////////////
1875/// return the bin names of unconnected bins
1876///
1877/// \param[in] bin local bin number
1879 TObjString const *r = nullptr;
1880 if(HasUnconnectedBins()) {
1881 if(bin<fAxisLabelList->GetEntriesFast()) {
1882 r = ((TObjString const *)fAxisLabelList->At(bin));
1883 }
1884 }
1885 return r;
1886}
1887
1888
1889////////////////////////////////////////////////////////////////////////
1891/// get average bin size on the specified axis
1892///
1893/// \param[in] axis axis number
1894/// \param[in] includeUnderflow whether to include the underflow bin
1895/// \param[in] includeOverflow whether to include the overflow bin
1896(Int_t axis,Bool_t includeUnderflow,Bool_t includeOverflow) const
1897{
1898 Double_t r=0.0;
1899 if((axis>=0)&&(axis<GetDistributionDimension())) {
1900 TVectorD const *bins=GetDistributionBinning(axis);
1901 Double_t d=(*bins)[bins->GetNrows()-1]-(*bins)[0];
1902 Double_t nBins=bins->GetNrows()-1;
1903 if(includeUnderflow && HasUnderflow(axis)) {
1904 Double_t w=GetDistributionUnderflowBinWidth(axis);
1905 if(w>0) {
1906 nBins++;
1907 d += w;
1908 }
1909 }
1910 if(includeOverflow && HasOverflow(axis)) {
1911 Double_t w=GetDistributionOverflowBinWidth(axis);
1912 if(w>0.0) {
1913 nBins++;
1914 d += w;
1915 }
1916 }
1917 if(nBins>0) {
1918 r=d/nBins;
1919 }
1920 } else {
1921 Error("GetDistributionAverageBinSize","axis %d does not exist",axis);
1922 }
1923 return r;
1924}
1925
1926////////////////////////////////////////////////////////////////////////
1927/// return bin width assigned to the underflow bin
1928///
1929/// \param[in] axis axis number
1930///
1931/// the bin width of the first bin is returned.
1932/// The method is virtual, so this behaviour can be adjusted.
1934{
1935 // return width of the underflow bin
1936 // axis: axis number
1937 TVectorD const *bins=GetDistributionBinning(axis);
1938 return (*bins)[1]-(*bins)[0];
1939}
1940
1941////////////////////////////////////////////////////////////////////////
1942/// return bin width assigned to the overflow bin
1943///
1944/// \param[in] axis axis number
1945///
1946/// the bin width of the last bin is returned.
1947/// The method is virtual, so this behaviour can be adjusted.
1949{
1950 // return width of the underflow bin
1951 // axis: axis number
1952 TVectorD const *bins=GetDistributionBinning(axis);
1953 return (*bins)[bins->GetNrows()-1]-(*bins)[bins->GetNrows()-2];
1954}
1955
1956////////////////////////////////////////////////////////////////////////
1957/// return bin center for a given axis and bin number
1958///
1959/// \param[in] axis axis number
1960/// \param[in] bin local bin number on the specified axis
1961///
1962/// returns the geometrical bin center.
1963/// for underflow and overflow, the calculation is using the
1964/// GetDistributionUnderflowBinWidth() and
1965/// GetDistributionOverflowBinWidth() methods.
1967(Int_t axis,Int_t bin) const
1968{
1969 // position of the bin center
1970 // input
1971 // axis : axis number
1972 // bin : bin number on the axis
1973 TVectorD const *bins=GetDistributionBinning(axis);
1974 Double_t r=0.0;
1975 if(bin<0) {
1976 // underflow bin
1977 r=(*bins)[0]-0.5*GetDistributionUnderflowBinWidth(axis);
1978 } else if(bin>=bins->GetNrows()-1) {
1979 // overflow bin
1980 r=(*bins)[bins->GetNrows()-1]+0.5*GetDistributionOverflowBinWidth(axis);
1981 } else {
1982 r=0.5*((*bins)[bin+1]+(*bins)[bin]);
1983 }
1984 return r;
1985}
1986
1987////////////////////////////////////////////////////////////////////////
1988/// get global bin number, given axis bin numbers
1989///
1990/// \param[in] axisBins[] bin numbers on each axis
1991/// \param[out] isBelow indicates bins are in underflow but there is
1992/// no undeflow bin
1993/// \param[out] isAbove indicates bins are in overflow but there is
1994/// no overflow bin
1995///
1996/// return: global bin nuber or -1 if not matched.
1998(Int_t const *axisBins,Int_t *isBelow,Int_t *isAbove) const
1999{
2000 Int_t dimension=GetDistributionDimension();
2001 Int_t r=0;
2002 if(isBelow) *isBelow=0;
2003 if(isAbove) *isAbove=0;
2004 if(dimension>0) {
2005 for(Int_t axis=dimension-1;axis>=0;axis--) {
2006 Int_t nMax=GetDistributionBinning(axis)->GetNrows()-1;
2007 Int_t i=axisBins[axis];
2008 if(HasUnderflow(axis)) {
2009 nMax +=1;
2010 i +=1;
2011 }
2012 if(HasOverflow(axis)) nMax +=1;
2013 if((i>=0)&&(i<nMax)) {
2014 if(r>=0) r = r*nMax +i;
2015 } else {
2016 r=-1;
2017 if((i<0)&&(isBelow)) *isBelow |= 1<<axis;
2018 if((i>=nMax)&&(isAbove)) *isAbove |= 1<<axis;
2019 }
2020 }
2021 if(r>=0) {
2022 r += GetStartBin();
2023 }
2024 } else {
2025 if((axisBins[0]>=0)&&(axisBins[0]<GetDistributionNumberOfBins()))
2026 r=GetStartBin()+axisBins[0];
2027 else
2028 Fatal("ToGlobalBin","bad input %d for dimensionless binning %s %d",
2029 axisBins[0],(const char *)GetName(),
2031 }
2032 return r;
2033}
2034
2035////////////////////////////////////////////////////////////////////////
2036/// return distribution in which the bin is located
2037/// and bin numbers on the corresponding axes
2038///
2039/// \param[in] globalBin global bin number
2040/// \param[out] local bin numbers of the distribution's axes
2041///
2042/// returns the distribution in which the globalBin is located
2043/// or 0 if the globalBin is outside this node and its children
2045(Int_t globalBin,Int_t *axisBins) const
2046{
2047 TUnfoldBinning const *r=nullptr;
2048 if((globalBin>=GetStartBin())&&(globalBin<GetEndBin())) {
2049 TUnfoldBinning const *node;
2050 for(node=GetChildNode();node && !r; node=node->GetNextNode()) {
2051 r=node->ToAxisBins(globalBin,axisBins);
2052 }
2053 if(!r) {
2054 r=this;
2055 Int_t i=globalBin-GetStartBin();
2056 Int_t dimension=GetDistributionDimension();
2057 if(dimension>0) {
2058 for(int axis=0;axis<dimension;axis++) {
2059 Int_t nMax=GetDistributionBinning(axis)->GetNrows()-1;
2060 axisBins[axis]=0;
2061 if(HasUnderflow(axis)) {
2062 axisBins[axis] =-1;
2063 nMax += 1;
2064 }
2065 if(HasOverflow(axis)) nMax +=1;
2066 axisBins[axis] += i % nMax;
2067 i /= nMax;
2068 }
2069 } else {
2070 axisBins[0]=i;
2071 }
2072 }
2073 }
2074 return r;
2075}
2076
2077////////////////////////////////////////////////////////////////////////
2078/// decode axis steering
2079///
2080/// \param[in] axisSteering the steering to decode
2081/// \param[in] options the allowed options to extract
2082/// \param[out] isOptionGiven array of decoded steering options,
2083/// the dimension equal to the number of
2084/// characters in <b>options</b>
2085///
2086/// the axis steering is given in the form
2087/// "axis[option];axis[option];..."
2088/// <br/> axis : the name of the axis for which the optionlist is relevant
2089/// the character * matches all axes
2090/// <br/>option : a list of characters taken from <b>options</b>
2091/// <br/>for each match the corresponding bit number corresponding to
2092/// the axis number is set in
2093/// <b>isOptionGiven</b>[i], where i is the position of the matching option
2094/// character in <b>options</b>
2096(const char *axisSteering,const char *options,Int_t *isOptionGiven) const
2097{
2098 Int_t nOpt=TString(options).Length();
2099 for(Int_t i=0;i<nOpt;i++) isOptionGiven[i]=0;
2100 if(axisSteering) {
2101 TObjArray *patterns=TString(axisSteering).Tokenize(";");
2102 Int_t nPattern=patterns->GetEntries();
2104 for(Int_t i=0;i<nPattern;i++) {
2105 TString const &pattern=((TObjString const *)patterns->At(i))
2106 ->GetString();
2107 Int_t bracketBegin=pattern.Last('[');
2108 Int_t len=pattern.Length();
2109 if((bracketBegin>0)&&(pattern[len-1]==']')) {
2110 TString axisId=pattern(0,bracketBegin);
2111 Int_t mask=0;
2112 if((axisId[0]=='*')&&(axisId.Length()==1)) {
2113 // turn all bins on
2114 mask=(1<<nAxis)-1;
2115 } else {
2116 // if axis is there, turn its bit on
2117 for(Int_t j=0;j<nAxis;j++) {
2118 if(!axisId.CompareTo(GetDistributionAxisLabel(j))) {
2119 mask|= (1<<j);
2120 }
2121 }
2122 }
2123 for(Int_t o=0;o<nOpt;o++) {
2124 if(pattern.Last(options[o])>bracketBegin) {
2125 isOptionGiven[o] |= mask;
2126 }
2127 }
2128 } else {
2129 Error("DecodeAxisSteering",
2130 "steering \"%s\" does not end with [options]",
2131 (const char *)pattern);
2132 }
2133 }
2134 }
2135}
2136
#define d(i)
Definition RSha256.hxx:102
#define c(i)
Definition RSha256.hxx:101
constexpr Bool_t kFALSE
Definition RtypesCore.h:94
constexpr Bool_t kTRUE
Definition RtypesCore.h:93
#define ClassImp(name)
Definition Rtypes.h:382
static void indent(ostringstream &buf, int indent_level)
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t mask
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t dest
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t child
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t UChar_t len
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
char name[80]
Definition TGX11.cxx:110
TVectorT< Double_t > TVectorD
Definition TVectorDfwd.h:23
Class to manage histogram axis.
Definition TAxis.h:31
const char * GetTitle() const override
Returns title of object.
Definition TAxis.h:135
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition TAxis.cxx:518
Int_t GetNbins() const
Definition TAxis.h:125
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition TAxis.cxx:528
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
1-Dim function class
Definition TF1.h:233
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:671
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:9096
virtual Int_t GetDimension() const
Definition TH1.h:284
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:4990
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5090
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:358
Service class for 2-D histogram classes.
Definition TH2.h:30
Double_t GetBinContent(Int_t binx, Int_t biny) const override
Definition TH2.h:94
3-D histogram with a double per channel (see TH1 documentation)
Definition TH3.h:364
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
An array of TObjects.
Definition TObjArray.h:31
Int_t GetEntriesFast() const
Definition TObjArray.h:58
Int_t GetEntries() const override
Return the number of objects in array (i.e.
TObject * At(Int_t idx) const override
Definition TObjArray.h:164
void AddLast(TObject *obj) override
Add object in the next empty slot in the array.
Collectable string class.
Definition TObjString.h:28
Mother of all ROOT objects.
Definition TObject.h:41
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:991
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1005
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition TObject.cxx:1033
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:979
Basic string class.
Definition TString.h:139
Ssiz_t Length() const
Definition TString.h:417
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition TString.cxx:457
Ssiz_t Last(char c) const
Find last occurrence of a character c.
Definition TString.cxx:931
TObjArray * Tokenize(const TString &delim) const
This function is used to isolate sequential tokens in a TString.
Definition TString.cxx:2264
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:2378
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
void PrintStream(std::ostream &out, Int_t indent=0, int debug=0) const
print some information about this binning tree
Int_t ToGlobalBin(Int_t const *axisBins, Int_t *isBelow=nullptr, Int_t *isAbove=nullptr) const
get global bin number, given axis bin numbers
Bool_t HasOverflow(int axis) const
check whether the axis has an overflow bin
TH1 * CreateHistogram(const char *histogramName, Bool_t originalAxisBinning=kFALSE, Int_t **binMap=nullptr, const char *histogramTitle=nullptr, const char *axisSteering=nullptr) const
create a THxx histogram capable to hold the bins of this binning node and its children
Int_t FillBinMapSingleNode(const TH1 *hist, Int_t startBin, Int_t nDim, const Int_t *axisList, const char *axisSteering, Int_t *binMap) const
fill bin map for a single node
virtual Double_t GetDistributionOverflowBinWidth(Int_t axis) const
return bin width assigned to the overflow bin
Int_t fFirstBin
global bin number of the first bin
void SetBinFactorFunction(Double_t normalisation, TF1 *userFunc=nullptr)
set normalisation factor and function which are used in calls to GetBinFactor()
TH1 * ExtractHistogram(const char *histogramName, const TH1 *globalBins, const TH2 *globalBinsEmatrix=nullptr, Bool_t originalAxisBinning=kTRUE, const char *axisSteering=nullptr) const
Int_t UpdateFirstLastBin(Bool_t startWithRootNode=kTRUE)
update fFirstBin and fLastBin members of this node and its children
TString BuildHistogramTitle(const char *histogramName, const char *histogramTitle, Int_t const *axisList) const
construct a title
Int_t FillBinMapRecursive(Int_t startBin, const char *axisSteering, Int_t *binMap) const
recursively fill bin map
Int_t FillBinMap1D(Int_t *binMap, const char *axisSteering, Int_t firstBinX) const
map all global bins referenced by this node to the one-dimensional histogram destHist,...
TUnfoldBinning * parentNode
mother node
Int_t GetDistributionDimension(void) const
query dimension of this node's distribution
TUnfoldBinning const * GetPrevNode(void) const
previous sister node
void SetBinFactor(Double_t normalisation, TObject *factors)
set normalisation factors which are used in calls to GetBinFactor()
TUnfoldBinning * childNode
first daughter node
TString GetBinName(Int_t iBin) const
get the name of a bin
Int_t GetTH1xNumberOfBins(Bool_t originalAxisBinning=kTRUE, const char *axisSteering=nullptr) const
return the number of histogram bins required when storing this binning in a one-dimensional histogram
virtual Double_t GetDistributionUnderflowBinWidth(Int_t axis) const
return bin width assigned to the underflow bin
Double_t GetBinSize(Int_t iBin) const
get N-dimensional bin size
Int_t GetDistributionNumberOfBins(void) const
number of bins in the distribution possibly including under/overflow
Double_t GetGlobalFactor(void) const
return global scaling factor for this node
Int_t GetGlobalBinNumber(Double_t x) const
locate a bin in a one-dimensional distribution
Bool_t HasUnconnectedBins(void) const
check whether there are bins but no axis
virtual Double_t GetBinFactor(Int_t iBin) const
return scaling factor for the given global bin number
TUnfoldBinning * nextNode
next sister
TObject * fBinFactorFunction
function to calculate a scale factor from bin centres (may be a TF1 or a TVectorD
Int_t GetTHxxBinning(Int_t maxDim, Int_t *axisBins, Int_t *axisList, const char *axisSteering) const
calculate properties of a THxx histogram to store this binning
const TUnfoldBinning * GetNonemptyNode_r(int &count) const
Int_t fHasUnderflow
bit fields indicating whether there are underflow bins on the axes
Int_t GetEndBin(void) const
last+1 bin of this node (includes children)
TObjArray * fAxisLabelList
for each axis its name (TObjString), or names of unconnected bins
Double_t fBinFactorConstant
common scale factor for all bins of this node
TUnfoldBinning const * GetNextNode(void) const
next sister node
Bool_t HasUnderflow(int axis) const
check whether an axis has an underflow bin
Int_t * CreateEmptyBinMap(void) const
create an empty bin map, useful together with the getter methods of class TUnfold and TUnfoldSys
TUnfoldBinning const * GetParentNode(void) const
mother node
@ MAXDIM
maximum numner of axes per distribution
static TH2D * CreateHistogramOfMigrations(TUnfoldBinning const *xAxis, TUnfoldBinning const *yAxis, char const *histogramName, Bool_t originalXAxisBinning=kFALSE, Bool_t originalYAxisBinning=kFALSE, char const *histogramTitle=nullptr)
create a TH2D histogram capable to hold the bins of the two input binning schemes on the x and y axes...
TObjArray * fAxisList
for each axis the bin borders (TVectorD)
Int_t fLastBin
global bin number of the last(+1) bin, including daughters
void GetBinUnderflowOverflowStatus(Int_t iBin, Int_t *uStatus, Int_t *oStatus) const
return bit maps indicating underflow and overflow status
virtual Double_t GetDistributionBinCenter(Int_t axis, Int_t bin) const
return bin center for a given axis and bin number
Int_t fDistributionSize
number of bins in this node's distribution
Int_t fHasOverflow
bit fields indicating whether there are overflow bins on the axes
const TUnfoldBinning * GetNonemptyNode(void) const
find a node which has non-empty distributions if there is none or if there are many,...
void Initialize(Int_t nBins)
initialize variables for a given number of bins
void SetBinMapEntry(Int_t *binMap, Int_t globalBin, Int_t destBin) const
set one entry in a bin map
virtual Double_t GetDistributionAverageBinSize(Int_t axis, Bool_t includeUnderflow, Bool_t includeOverflow) const
get average bin size on the specified axis
void DecodeAxisSteering(const char *axisSteering, const char *options, Int_t *isOptionGiven) const
decode axis steering
virtual Bool_t IsBinFactorGlobal(void) const
check whether there is only a global scaling factor for this node
TString GetDistributionAxisLabel(Int_t axis) const
get name of an axis
TUnfoldBinning * AddBinning(TUnfoldBinning *binning)
add a TUnfoldBinning as the last child of this node
const TObjString * GetUnconnectedBinName(Int_t bin) const
return the bin names of unconnected bins
Int_t GetBinNeighbours(Int_t globalBin, Int_t axis, Int_t *prev, Double_t *distPrev, Int_t *next, Double_t *distNext, Bool_t isPeriodic=kFALSE) const
get neighbour bins along the specified axis
TVectorD const * GetDistributionBinning(Int_t axis) const
get vector of bin borders for one axis
Int_t * CreateBinMap(const TH1 *hist, Int_t nDim, const Int_t *axisList, const char *axisSteering) const
create mapping from global bin number to a histogram for this node
TH2D * CreateErrorMatrixHistogram(const char *histogramName, Bool_t originalAxisBinning, Int_t **binMap=nullptr, const char *histogramTitle=nullptr, const char *axisSteering=nullptr) const
create a TH2D histogram capable to hold a covariance matrix
Bool_t AddAxis(const char *name, Int_t nBins, const Double_t *binBorders, Bool_t hasUnderflow, Bool_t hasOverflow)
add an axis with the specified bin borders
TUnfoldBinning * prevNode
previous sister
TUnfoldBinning const * GetRootNode(void) const
return root node of the binnig scheme
Int_t GetStartBin(void) const
first bin of this node
TUnfoldBinning const * ToAxisBins(Int_t globalBin, Int_t *axisBins) const
return distribution in which the bin is located and bin numbers on the corresponding axes
TUnfoldBinning const * FindNode(char const *name) const
traverse the tree and return the first node which matches the given name
TUnfoldBinning(const char *name=nullptr, Int_t nBins=0, const char *binNames=nullptr)
create a new node without axis
Int_t GetTHxxBinsRecursive(const char *axisSteering) const
calculate number of bins required to store this binning with teh given axisSteering
TUnfoldBinning const * GetChildNode(void) const
first daughter node
Int_t GetTHxxBinningSingleNode(Int_t maxDim, Int_t *axisBins, Int_t *axisList, const char *axisSteering) const
get the properties of a histogram capable to hold the distribution attached to this node
TString BuildHistogramTitle2D(const char *histogramName, const char *histogramTitle, Int_t xAxis, const TUnfoldBinning *yAxisBinning, Int_t yAxis) const
construct a histogram title for a 2D histogram with different binning schemes on x and y axis
Int_t GetNrows() const
Definition TVectorT.h:73
Element * GetMatrixArray()
Definition TVectorT.h:76
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
Int_t Finite(Double_t x)
Check if it is finite with a mask in order to be consistent in presence of fast math.
Definition TMath.h:774
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:666
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:725