Logo ROOT   6.12/07
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
6 Binning schemes for use with the unfolding algorithm TUnfoldDensity.
7 
8 Binning schemes are used to map analysis bins on a single histogram
9 axis and back. The analysis bins may include unconnected bins (e.g
10 nuisances for background normalisation) or various multidimensional
11 histograms (signal bins, differential background normalisation bins, etc).
12 
13 If you use this software, please consider the following citation
14 
15 <b>S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]</b>
16 
17 Detailed documentation and updates are available on
18 http://www.desy.de/~sschmitt
19 
20 ### Functionality
21 
22 The TUnfoldBinning objects are connected by a tree-like structure.
23 The structure does not hold any data, but is only responsible for
24 arranging the analysis bins in the proper order.
25 Each node of the tree is responsible for a group of bins. That group
26 may consist of
27 
28  - several unconnected bins, each with a dedicated name.
29  - bins organized in a multidimensional distribution, defined by a
30 set of axes. The axes are defined by a number of bins N and by (N+1)
31 bin borders. In addition to the N bins inside there may be an underflow and an
32 overflow bin
33 
34 Each bin has a "global" bin number, which can be found using the
35 GetGlobalBinNumber() methods. The global bin number 0 is reserved and
36 corresponds to the case where no bin is found in the
37 TUnfoldBinning tree.
38 
39 ### Use in the analysis
40 Booking histograms:
41 
42  - Define binning schemes on detector level and on truth level. This
43 can be done using the XML language, use the class TUnfoldBinningXML to
44 read the binning scheme. The TUnfoldBinning objects can be written to
45 a root file, preferentially together with the corresponding histograms.
46  - For Monte Carlo, book histograms for the response matrix (detector
47 vs truth level) using the
48 method 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 
56 The histograms which are booked have all analysis bins arranged on one
57 axis (global bin number). TUnfoldBinning provides methods to locate
58 the global bin number:
59 
60  - Use the method FindNode() to locate a group of bins (e.g. signal,
61 control distribution, etc) by their name, then:
62  - Use the method GetGlobalBinNumber() to locate a bin in a
63 distribution, then:
64  - Use the TH1::Fill() method and the bin number to fill the
65 appropriate bin in one of the histograms booked above.
66 
67 Unfolding: Specify the response matrix and the binning schemes when
68 constructing a TUnfoldDensity object. Tell TUnfoldDensity about the
69 data, background, systematic error histograms using the corresponding
70 methods of class TUnfoldDensity. Then run the unfolding. Use the
71 GetXXX() methods to retrieve the unfolding results into properly
72 binned 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
92 CreateHistogramOfMigrations(). 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 
118 using 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)) {
131  parentNode->childNode=nextNode;
132  }
133  if(GetPrevNode()) prevNode->nextNode=nextNode;
134  if(GetNextNode()) nextNode->prevNode=prevNode;
135  delete fAxisList;
136  delete fAxisLabelList;
137  if(fBinFactorFunction) {
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();
155  fAxisLabelList=new TObjArray();
156  fAxisList->SetOwner();
157  fAxisLabelList->SetOwner();
158  fHasUnderflow=0;
159  fHasOverflow=0;
160  fDistributionSize=nBins;
161  fBinFactorFunction=0;
162  fBinFactorConstant=1.0;
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) {
173  return GetRootNode()->UpdateFirstLastBin(kFALSE);
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
178  fFirstBin=GetPrevNode()->GetEndBin();
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
182  fFirstBin=GetParentNode()->GetStartBin()+
183  GetParentNode()->GetDistributionNumberOfBins();
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  }
196  fLastBin=fFirstBin+fDistributionSize;
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  }
221  UpdateFirstLastBin();
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);
237  UpdateFirstLastBin();
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) {
288  TUnfoldBinning *child=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  }
299  UpdateFirstLastBin();
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 {
321  Bool_t r=kFALSE;
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 {
382  Bool_t r=kFALSE;
383  if(HasUnconnectedBins()) {
384  Fatal("AddAxis","node already has %d bins without axis",
385  GetDistributionNumberOfBins());
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) {
404  Int_t axis=fAxisList->GetEntriesFast();
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);
420  fAxisLabelList->AddLast(new TObjString(name));
421  if(!fDistributionSize) fDistributionSize=1;
422  fDistributionSize *= nBinUO;
423  UpdateFirstLastBin();
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 
436 void 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="
448  <<GetTH1xNumberOfBins()
449  <<"\n";
450  if(GetDistributionNumberOfBins()) {
451  for(Int_t i=0;i<indent;i++) out<<" ";
452  out<<" distribution: "<<GetDistributionNumberOfBins()<<" bins\n";
453  if(fAxisList->GetEntriesFast()) {
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<<" \""
459  <<GetDistributionAxisLabel(axis)
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;
513  if(fBinFactorFunction) {
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 
542 TUnfoldBinning const *TUnfoldBinning::FindNode(char const *name) const
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) {
646  r += GetNonemptyNode()->GetDistributionAxisLabel(xAxis);
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
993  if(HasUnconnectedBins() || (GetDistributionNumberOfBins()<=0)) {
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;
1000  Int_t nBinI=GetDistributionBinning(i)->GetNrows()-1;
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 
1020 Int_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 {
1147  Int_t *r=CreateEmptyBinMap();
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 {
1349  nbin=GetDistributionNumberOfBins();
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",
1438  GetDistributionDimension());
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",
1458  GetDistributionDimension());
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 
1483 {
1484  if(GetDistributionDimension()!=3) {
1485  Fatal("GetBinNumber",
1486  "called with 3 arguments for %d dimensional distribution",
1487  GetDistributionDimension());
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",
1518  GetDistributionDimension());
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",
1551  GetDistributionDimension());
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",
1586  GetDistributionDimension());
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 {
1630  if(!GetDistributionDimension()) {
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 {
1900  return (!GetDistributionDimension())&&(GetDistributionNumberOfBins()>0);
1901 }
1902 
1903 ////////////////////////////////////////////////////////////////////////////////
1904 /// Return the bin names of unconnected bins.
1905 ///
1906 /// \param[in] bin local bin number
1907 
1909  TObjString *r=0;
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)) {
1935  Double_t w=GetDistributionUnderflowBinWidth(axis);
1936  if(w>0) {
1937  nBins++;
1938  d += w;
1939  }
1940  }
1941  if(includeOverflow && HasOverflow(axis)) {
1942  Double_t w=GetDistributionOverflowBinWidth(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(),
2068  GetDistributionNumberOfBins());
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();
2143  Int_t nAxis=fAxisLabelList->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 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
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, starting with bin firstBinX.
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:145
const TObjString * GetUnconnectedBinName(Int_t bin) const
Return the bin names of unconnected bins.
An array of TObjects.
Definition: TObjArray.h:37
TUnfoldBinning * parentNode
mother node
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.
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.
virtual Double_t GetDistributionBinCenter(Int_t axis, Int_t bin) const
return bin center for a given axis and bin number
const TUnfoldBinning * GetNonemptyNode(void) const
Find a node which has non-empty distributions if there is none or if there are many, return zero.
void Fatal(const char *location, const char *msgfmt,...)
Collectable string class.
Definition: TObjString.h:28
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:504
virtual Double_t GetDistributionAverageBinSize(Int_t axis, Bool_t includeUnderflow, Bool_t includeOverflow) const
Get average bin size on the specified 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.
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...
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...
Double_t GetBinSize(Int_t iBin) const
Get N-dimensional bin size.
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4763
TUnfoldBinning * nextNode
next sister
void GetBinUnderflowOverflowStatus(Int_t iBin, Int_t *uStatus, Int_t *oStatus) const
Return bit maps indicating underflow and overflow status.
Int_t GetNrows() const
Definition: TVectorT.h:75
void SetBinFactor(Double_t normalisation, TObject *factors)
Set normalisation factors which are used in calls to GetBinFactor().
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
STL namespace.
void PrintStream(std::ostream &out, Int_t indent=0, int debug=0) const
Print some information about this binning tree.
TObject * At(Int_t idx) const
Definition: TObjArray.h:165
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.
virtual ~TUnfoldBinning(void)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:627
TUnfoldBinning const * GetNextNode(void) const
next sister node
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:514
Int_t UpdateFirstLastBin(Bool_t startWithRootNode=kTRUE)
Update fFirstBin and fLastBin members of this node and its children.
Int_t ToGlobalBin(Int_t const *axisBins, Int_t *isBelow=0, Int_t *isAbove=0) const
Get global bin number, given axis bin numbers.
Double_t fBinFactorConstant
common scale factor for all bins of this node
virtual Int_t GetDimension() const
Definition: TH1.h:277
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
Int_t GetEndBin(void) const
last+1 bin of this node (includes children)
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:2365
void SetBinMapEntry(Int_t *binMap, Int_t globalBin, Int_t destBin) const
Set one entry in a bin map.
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
Double_t GetGlobalFactor(void) const
Return global scaling factor for this node.
void Initialize(Int_t nBins)
Initialize variables for a given number of bins.
static const double x4[22]
void Info(const char *location, const char *msgfmt,...)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:146
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.
Int_t Finite(Double_t x)
Definition: TMath.h:654
Int_t GetStartBin(void) const
first bin of this node
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition: TH1.cxx:8461
void Error(const char *location, const char *msgfmt,...)
TVectorT< Double_t > TVectorD
Definition: TVectorDfwd.h:22
TVectorD const * GetDistributionBinning(Int_t axis) const
get vector of bin borders for one axis
Int_t GetDistributionDimension(void) const
query dimension of this node&#39;s distribution
Element * GetMatrixArray()
Definition: TVectorT.h:78
TUnfoldBinning const * GetParentNode(void) const
mother node
THist< 3, double, THistStatContent, THistStatUncertainty > TH3D
Definition: THist.hxx:296
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.
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.
void DecodeAxisSteering(const char *axisSteering, const char *options, Int_t *isOptionGiven) const
Decode axis steering.
void SetBinFactorFunction(Double_t normalisation, TF1 *userFunc=0)
Set normalisation factor and function which are used in calls to GetBinFactor().
ROOT::R::TRInterface & r
Definition: Object.C:4
const char * GetTitle() const
Returns title of object.
Definition: TAxis.h:129
Service class for 2-Dim histogram classes.
Definition: TH2.h:30
Class to manage histogram axis.
Definition: TAxis.h:30
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.cxx:176
virtual Double_t GetDistributionUnderflowBinWidth(Int_t axis) const
Return bin width assigned to the underflow bin.
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8477
virtual Bool_t IsBinFactorGlobal(void) const
Check whether there is only a global scaling factor for this node.
TUnfoldBinning const * GetPrevNode(void) const
previous sister node
void Warning(const char *location, const char *msgfmt,...)
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:4665
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...
const Bool_t kFALSE
Definition: RtypesCore.h:88
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
TString GetDistributionAxisLabel(Int_t axis) const
get name of an axis
static const double x1[5]
#define ClassImp(name)
Definition: Rtypes.h:359
TObject * fBinFactorFunction
function to calculate a scale factor from bin centres (may be a TF1 or a TVectorD ...
double Double_t
Definition: RtypesCore.h:55
TUnfoldBinning * prevNode
previous sister
Double_t y[n]
Definition: legend1.C:17
The TH1 histogram class.
Definition: TH1.h:56
Int_t * CreateEmptyBinMap(void) const
Create an empty bin map, useful together with the getter methods of class TUnfold and TUnfoldSys...
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
Definition: THist.hxx:290
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.
Int_t GetGlobalBinNumber(Double_t x) const
Locate a bin in a one-dimensional distribution.
TString GetBinName(Int_t iBin) const
Get the name of a bin.
Mother of all ROOT objects.
Definition: TObject.h:37
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
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 FillBinMapRecursive(Int_t startBin, const char *axisSteering, Int_t *binMap) const
Recursively fill bin map.
Bool_t HasUnconnectedBins(void) const
Check whether there are bins but no axis.
TUnfoldBinning const * GetRootNode(void) const
return root node of the binnig scheme
TObjArray * fAxisList
for each axis the bin borders (TVectorD)
virtual Double_t GetDistributionOverflowBinWidth(Int_t axis) const
Return bin width assigned to the overflow bin.
#define dest(otri, vertexptr)
Definition: triangle.c:1040
1-Dim function class
Definition: TF1.h:211
TUnfoldBinning(const char *name=0, Int_t nBins=0, const char *binNames=0)
Create a new node without axis.
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:284
Int_t GetEntries() const
Return the number of objects in array (i.e.
Definition: TObjArray.cxx:522
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...
Double_t Sqrt(Double_t x)
Definition: TMath.h:590
Int_t GetNbins() const
Definition: TAxis.h:121
TString BuildHistogramTitle(const char *histogramName, const char *histogramTitle, Int_t const *axisList) const
Construct a title.
Int_t GetTHxxBinsRecursive(const char *axisSteering) const
Calculate number of bins required to store this binning with the given axisSteering.
const Bool_t kTRUE
Definition: RtypesCore.h:87
TUnfoldBinning * AddBinning(TUnfoldBinning *binning)
Add a TUnfoldBinning as the last child of this node.
static char * skip(char **buf, const char *delimiters)
Definition: civetweb.c:2039
TUnfoldBinning const * FindNode(char const *name) const
Traverse the tree and return the first node which matches the given name.
char name[80]
Definition: TGX11.cxx:109
virtual Double_t GetBinFactor(Int_t iBin) const
Return scaling factor for the given global bin number.
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8319
static const double x3[11]
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:290