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