00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 #include "TUnfoldBinningXML.h"
00105 #include <TVectorD.h>
00106 #include <TAxis.h>
00107 #include <TString.h>
00108 #include <TMath.h>
00109 #include <TF1.h>
00110 #include <TH1D.h>
00111 #include <TH2D.h>
00112 #include <TH3D.h>
00113 #include <TIterator.h>
00114 #include <iomanip>
00115
00116
00117
00118 using namespace std;
00119
00120 ClassImp(TUnfoldBinningV17)
00121
00122
00123 TUnfoldBinningV17::~TUnfoldBinningV17(void)
00124 {
00125
00126 while(childNode) delete childNode;
00127
00128 if(GetParentNode() && (GetParentNode()->GetChildNode()==this)) {
00129 parentNode->childNode=nextNode;
00130 }
00131 if(GetPrevNode()) prevNode->nextNode=nextNode;
00132 if(GetNextNode()) nextNode->prevNode=prevNode;
00133 delete fAxisList;
00134 delete fAxisLabelList;
00135 if(fBinFactorFunction) {
00136 if(!dynamic_cast<TF1 *>(fBinFactorFunction))
00137 delete fBinFactorFunction;
00138 }
00139 }
00140
00141
00142
00143
00144
00145 void TUnfoldBinningV17::Initialize(Int_t nBins)
00146 {
00147 parentNode=0;
00148 childNode=0;
00149 nextNode=0;
00150 prevNode=0;
00151 fAxisList=new TObjArray();
00152 fAxisLabelList=new TObjArray();
00153 fAxisList->SetOwner();
00154 fAxisLabelList->SetOwner();
00155 fHasUnderflow=0;
00156 fHasOverflow=0;
00157 fDistributionSize=nBins;
00158 fBinFactorFunction=0;
00159 fBinFactorConstant=1.0;
00160 }
00161
00162
00163
00164
00165
00166 Int_t TUnfoldBinningV17::UpdateFirstLastBin(Bool_t startWithRootNode)
00167 {
00168 if(startWithRootNode) {
00169 return GetRootNode()->UpdateFirstLastBin(kFALSE);
00170 }
00171 if(GetPrevNode()) {
00172
00173
00174 fFirstBin=GetPrevNode()->GetEndBin();
00175 } else if(GetParentNode()) {
00176
00177
00178 fFirstBin=GetParentNode()->GetStartBin()+
00179 GetParentNode()->GetDistributionNumberOfBins();
00180 } else {
00181
00182 fFirstBin=1;
00183
00184
00185
00186
00187 if((!GetChildNode())&&(GetDistributionDimension()==1)&&
00188 (fHasUnderflow==1)) {
00189 fFirstBin=0;
00190 }
00191 }
00192 fLastBin=fFirstBin+fDistributionSize;
00193
00194 for(TUnfoldBinningV17 *node=childNode;node;node=node->nextNode) {
00195 fLastBin=node->UpdateFirstLastBin(kFALSE);
00196 }
00197 return fLastBin;
00198 }
00199
00200
00201
00202
00203
00204
00205
00206 TUnfoldBinningV17::TUnfoldBinningV17
00207 (const char *name,Int_t nBins,const char *binNames)
00208 : TNamed(name ? name : "",name ? name : "")
00209 {
00210 Initialize(nBins);
00211 if(binNames) {
00212 TString nameString(binNames);
00213 delete fAxisLabelList;
00214 fAxisLabelList=nameString.Tokenize(";");
00215 }
00216 UpdateFirstLastBin();
00217 }
00218
00219
00220
00221
00222
00223
00224
00225 TUnfoldBinningV17::TUnfoldBinningV17
00226 (const TAxis &axis,Int_t includeUnderflow,Int_t includeOverflow)
00227 : TNamed(axis.GetName(),axis.GetTitle())
00228 {
00229 Initialize(0);
00230 AddAxis(axis,includeUnderflow,includeOverflow);
00231 UpdateFirstLastBin();
00232 }
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242 TUnfoldBinningV17 *TUnfoldBinningV17::AddBinning
00243 (const char *name,Int_t nBins,const char *binNames)
00244 {
00245 return AddBinning(new TUnfoldBinning(name,nBins,binNames));
00246 }
00247
00248
00249
00250
00251
00252
00253
00254
00255 TUnfoldBinningV17 *TUnfoldBinningV17::AddBinning(TUnfoldBinningV17 *binning)
00256 {
00257 TUnfoldBinning *r=0;
00258 if(binning->GetParentNode()) {
00259 Error("AddBinning",
00260 "binning \"%s\" already has parent \"%s\", can not be added to %s",
00261 (char *)binning->GetName(),
00262 (char *)binning->GetParentNode()->GetName(),
00263 (char *)GetName());
00264 } else if(binning->GetPrevNode()) {
00265 Error("AddBinning",
00266 "binning \"%s\" has previous node \"%s\", can not be added to %s",
00267 (char *)binning->GetName(),
00268 (char *)binning->GetPrevNode()->GetName(),
00269 (char *)GetName());
00270 } else if(binning->GetNextNode()) {
00271 Error("AddBinning",
00272 "binning \"%s\" has next node \"%s\", can not be added to %s",
00273 (char *)binning->GetName(),
00274 (char *)binning->GetNextNode()->GetName(),
00275 (char *)GetName());
00276 } else {
00277 r=binning;
00278 binning->parentNode=this;
00279 if(childNode) {
00280 TUnfoldBinningV17 *child=childNode;
00281
00282 while(child->nextNode) {
00283 child=child->nextNode;
00284 }
00285
00286 child->nextNode=r;
00287 r->prevNode=child;
00288 } else {
00289 childNode=r;
00290 }
00291 UpdateFirstLastBin();
00292 r=binning;
00293 }
00294 return r;
00295 }
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308 Bool_t TUnfoldBinningV17::AddAxis
00309 (const char *name,Int_t nBin,Double_t xMin,Double_t xMax,
00310 Bool_t hasUnderflow,Bool_t hasOverflow)
00311 {
00312 Bool_t r=kFALSE;
00313 if(nBin<=0) {
00314 Fatal("AddAxis","number of bins %d is not positive",
00315 nBin);
00316 } else if((!TMath::Finite(xMin))||(!TMath::Finite(xMax))||
00317 (xMin>=xMax)) {
00318 Fatal("AddAxis","xmin=%f required to be smaller than xmax=%f",
00319 xMin,xMax);
00320 } else {
00321 Double_t *binBorders=new Double_t[nBin+1];
00322 Double_t x=xMin;
00323 Double_t dx=(xMax-xMin)/nBin;
00324 for(Int_t i=0;i<=nBin;i++) {
00325 binBorders[i]=x+i*dx;
00326 }
00327 r=AddAxis(name,nBin,binBorders,hasUnderflow,hasOverflow);
00328 delete [] binBorders;
00329 }
00330 return r;
00331 }
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343 Bool_t TUnfoldBinningV17::AddAxis
00344 (const TAxis &axis,Bool_t hasUnderflow,Bool_t hasOverflow)
00345 {
00346 Int_t nBin=axis.GetNbins();
00347 Double_t *binBorders=new Double_t[nBin+1];
00348 for(Int_t i=0;i<nBin;i++) {
00349 binBorders[i]=axis.GetBinLowEdge(i+1);
00350 }
00351 binBorders[nBin]=axis.GetBinUpEdge(nBin);
00352 Bool_t r=AddAxis(axis.GetTitle(),nBin,binBorders,hasUnderflow,hasOverflow);
00353 delete [] binBorders;
00354 return r;
00355 }
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368 Bool_t TUnfoldBinningV17::AddAxis
00369 (const char *name,Int_t nBin,const Double_t *binBorders,
00370 Bool_t hasUnderflow,Bool_t hasOverflow)
00371 {
00372 Bool_t r=kFALSE;
00373 if(HasUnconnectedBins()) {
00374 Fatal("AddAxis","node already has %d bins without axis",
00375 GetDistributionNumberOfBins());
00376 } else if(nBin<=0) {
00377 Fatal("AddAxis","number of bins %d is not positive",
00378 nBin);
00379 } else {
00380 TVectorD *bins=new TVectorD(nBin+1);
00381 r=kTRUE;
00382 for(Int_t i=0;i<=nBin;i++) {
00383 (*bins)(i)=binBorders[i];
00384 if(!TMath::Finite((*bins)(i))) {
00385 Fatal("AddAxis","bin border %d is not finite",i);
00386 r=kFALSE;
00387 } else if((i>0)&&((*bins)(i)<=(*bins)(i-1))) {
00388 Fatal("AddAxis","bins not in order x[%d]=%f <= %f=x[%d]",
00389 i,(*bins)(i),(*bins)(i-1),i-1);
00390 r=kFALSE;
00391 }
00392 }
00393 if(r) {
00394 Int_t axis=fAxisList->GetEntriesFast();
00395 Int_t bitMask=1<<axis;
00396 Int_t nBinUO=nBin;
00397 if(hasUnderflow) {
00398 fHasUnderflow |= bitMask;
00399 nBinUO++;
00400 } else {
00401 fHasUnderflow &= ~bitMask;
00402 }
00403 if(hasOverflow) {
00404 fHasOverflow |= bitMask;
00405 nBinUO++;
00406 } else {
00407 fHasOverflow &= ~bitMask;
00408 }
00409 fAxisList->AddLast(bins);
00410 fAxisLabelList->AddLast(new TObjString(name));
00411 if(!fDistributionSize) fDistributionSize=1;
00412 fDistributionSize *= nBinUO;
00413 UpdateFirstLastBin();
00414 }
00415 }
00416 return r;
00417 }
00418
00419
00420
00421
00422
00423
00424
00425
00426 void TUnfoldBinningV17::PrintStream(ostream &out,Int_t indent,int debug)
00427 const {
00428 for(Int_t i=0;i<indent;i++) out<<" ";
00429 out<<"TUnfoldBinning \""<<GetName()<<"\" has ";
00430 Int_t nBin=GetEndBin()-GetStartBin();
00431 if(nBin==1) {
00432 out<<"1 bin";
00433 } else {
00434 out<<nBin<<" bins";
00435 }
00436 out<<" ["
00437 <<GetStartBin()<<","<<GetEndBin()<<"] nTH1x="
00438 <<GetTH1xNumberOfBins()
00439 <<"\n";
00440 if(GetDistributionNumberOfBins()) {
00441 for(Int_t i=0;i<indent;i++) out<<" ";
00442 out<<" distribution: "<<GetDistributionNumberOfBins()<<" bins\n";
00443 if(fAxisList->GetEntriesFast()) {
00444
00445
00446 for(Int_t axis=0;axis<GetDistributionDimension();axis++) {
00447 for(Int_t i=0;i<indent;i++) out<<" ";
00448 out<<" \""
00449 <<GetDistributionAxisLabel(axis)
00450 <<"\" nbin="<<GetDistributionBinning(axis)->GetNrows()-1;
00451 if(HasUnderflow(axis)) out<<" plus underflow";
00452 if(HasOverflow(axis)) out<<" plus overflow";
00453 out<<"\n";
00454 }
00455 } else {
00456 for(Int_t i=0;i<indent;i++) out<<" ";
00457 out<<" no axis\n";
00458 for(Int_t i=0;i<indent;i++) out<<" ";
00459 out<<" names: ";
00460 for(Int_t ibin=0;(ibin<GetDistributionNumberOfBins())&&
00461 (ibin<fAxisLabelList->GetEntriesFast());ibin++) {
00462 if(ibin) out<<";";
00463 if(GetDistributionAxisLabel(ibin)) {
00464 out<<GetDistributionAxisLabel(ibin);
00465 }
00466 }
00467 out<<"\n";
00468 }
00469 if(debug>0) {
00470
00471 for(int iBin=GetStartBin();iBin<GetEndBin();iBin++) {
00472 for(Int_t i=0;i<indent;i++) out<<" ";
00473 out<<GetBinName(iBin)
00474 <<" size="<<GetBinSize(iBin)
00475 <<" factor="<<GetBinFactor(iBin);
00476 out<<"\n";
00477 }
00478 }
00479 }
00480 TUnfoldBinningV17 const *child=GetChildNode();
00481 if(child) {
00482 while(child) {
00483 child->PrintStream(out,indent+1,debug);
00484 child=child->GetNextNode();
00485 }
00486 }
00487 }
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499 void TUnfoldBinningV17::SetBinFactor
00500 (Double_t normalisation,TObject *binfactor) {
00501 fBinFactorConstant=normalisation;
00502 if(fBinFactorFunction) {
00503 if(!dynamic_cast<TF1 *>(fBinFactorFunction))
00504 delete fBinFactorFunction;
00505 }
00506 fBinFactorFunction=binfactor;
00507 }
00508
00509
00510
00511
00512
00513
00514
00515 void TUnfoldBinningV17::SetBinFactorFunction
00516 (Double_t normalisation,TF1 *userFunc) {
00517 SetBinFactor(normalisation,userFunc);
00518 }
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529 TUnfoldBinningV17 const *TUnfoldBinningV17::FindNode(char const *name) const
00530 {
00531 TUnfoldBinningV17 const *r=0;
00532 if((!name)||(!TString(GetName()).CompareTo(name))) {
00533 r=this;
00534 }
00535 for(TUnfoldBinningV17 const *child=GetChildNode();
00536 (!r) && child;child=child->GetNextNode()) {
00537 r=child->FindNode(name);
00538 }
00539 return r;
00540 }
00541
00542
00543 TUnfoldBinningV17 *TUnfoldBinningV17::GetRootNode(void)
00544 {
00545
00546 TUnfoldBinningV17 *node=this;
00547 while(node->GetParentNode()) node=node->parentNode;
00548 return node;
00549 }
00550
00551
00552 TUnfoldBinningV17 const *TUnfoldBinningV17::GetRootNode(void) const
00553 {
00554
00555 TUnfoldBinningV17 const *node=this;
00556 while(node->GetParentNode()) node=node->GetParentNode();
00557 return node;
00558 }
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573 TString TUnfoldBinningV17::BuildHistogramTitle
00574 (const char *histogramName,const char *histogramTitle,Int_t const *axisList)
00575 const
00576 {
00577 TString r;
00578 if(histogramTitle) {
00579 r=histogramTitle;
00580 } else {
00581 r=histogramName;
00582 Int_t iEnd;
00583 for(iEnd=2;iEnd>0;iEnd--) {
00584 if(axisList[iEnd]>=0) break;
00585 }
00586 for(Int_t i=0;i<=iEnd;i++) {
00587 r += ";";
00588 if(axisList[i]<0) {
00589 r += GetName();
00590 } else {
00591 r += GetNonemptyNode()->GetDistributionAxisLabel(axisList[i]);
00592 }
00593 }
00594 }
00595 return r;
00596 }
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608 TString TUnfoldBinningV17::BuildHistogramTitle2D
00609 (const char *histogramName,const char *histogramTitle,
00610 Int_t xAxis,const TUnfoldBinningV17 *yAxisBinning,Int_t yAxis) const
00611 {
00612
00613
00614
00615
00616
00617
00618
00619 TString r;
00620 if(histogramTitle) {
00621 r=histogramTitle;
00622 } else {
00623 r=histogramName;
00624 r += ";";
00625 if(xAxis==-1) {
00626 r += GetName();
00627 } else if(xAxis>=0) {
00628 r += GetNonemptyNode()->GetDistributionAxisLabel(xAxis);
00629 }
00630 r+= ";";
00631 if(yAxis==-1) {
00632 r += yAxisBinning->GetName();
00633 } else if(yAxis>=0) {
00634 r += yAxisBinning->GetNonemptyNode()->GetDistributionAxisLabel(yAxis);
00635 }
00636
00637 }
00638 return r;
00639 }
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664 Int_t TUnfoldBinningV17::GetTH1xNumberOfBins
00665 (Bool_t originalAxisBinning,const char *axisSteering) const
00666 {
00667 Int_t axisBins[3],axisList[3];
00668 GetTHxxBinning(originalAxisBinning ? 1 : 0,axisBins,axisList,
00669 axisSteering);
00670 return axisBins[0];
00671 }
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724 TH1 *TUnfoldBinningV17::CreateHistogram
00725 (const char *histogramName,Bool_t originalAxisBinning,Int_t **binMap,
00726 const char *histogramTitle,const char *axisSteering) const
00727 {
00728 Int_t nBin[3],axisList[3];
00729 Int_t nDim=GetTHxxBinning(originalAxisBinning ? 3 : 0,nBin,axisList,
00730 axisSteering);
00731 const TUnfoldBinningV17 *neNode=GetNonemptyNode();
00732 TString title=BuildHistogramTitle(histogramName,histogramTitle,axisList);
00733 TH1 *r=0;
00734 if(nDim>0) {
00735 const TVectorD *axisBinsX=
00736 neNode->GetDistributionBinning(axisList[0]);
00737 if(nDim>1) {
00738 const TVectorD *axisBinsY=
00739 neNode->GetDistributionBinning(axisList[1]);
00740 if(nDim>2) {
00741 const TVectorD *axisBinsZ=
00742 neNode->GetDistributionBinning(axisList[2]);
00743 r=new TH3D(histogramName,title,
00744 nBin[0],axisBinsX->GetMatrixArray(),
00745 nBin[1],axisBinsY->GetMatrixArray(),
00746 nBin[2],axisBinsZ->GetMatrixArray());
00747 } else {
00748 r=new TH2D(histogramName,title,
00749 nBin[0],axisBinsX->GetMatrixArray(),
00750 nBin[1],axisBinsY->GetMatrixArray());
00751 }
00752 } else {
00753 r=new TH1D(histogramName,title,nBin[0],axisBinsX->GetMatrixArray());
00754 }
00755 } else {
00756 if(originalAxisBinning) {
00757 Warning("CreateHistogram",
00758 "Original binning can not be represented as THxx");
00759 }
00760 r=new TH1D(histogramName,title,nBin[0],0.5,nBin[0]+0.5);
00761 nDim=0;
00762 }
00763 if(binMap) {
00764 *binMap=CreateBinMap(r,nDim,axisList,axisSteering);
00765 }
00766 return r;
00767 }
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785 TH2D *TUnfoldBinningV17::CreateErrorMatrixHistogram
00786 (const char *histogramName,Bool_t originalAxisBinning,Int_t **binMap,
00787 const char *histogramTitle,const char *axisSteering) const
00788 {
00789 Int_t nBin[3],axisList[3];
00790 Int_t nDim=GetTHxxBinning(originalAxisBinning ? 1 : 0,nBin,axisList,
00791 axisSteering);
00792 TString title=BuildHistogramTitle(histogramName,histogramTitle,axisList);
00793 TH2D *r=0;
00794 if(nDim==1) {
00795 const TVectorD *axisBinsX=(TVectorD const *)
00796 GetNonemptyNode()->fAxisList->At(axisList[0]);
00797 r=new TH2D(histogramName,title,nBin[0],axisBinsX->GetMatrixArray(),
00798 nBin[0],axisBinsX->GetMatrixArray());
00799 } else {
00800 if(originalAxisBinning) {
00801 Info("CreateErrorMatrixHistogram",
00802 "Original binning can not be represented on one axis");
00803 }
00804 r=new TH2D(histogramName,title,nBin[0],0.5,nBin[0]+0.5,
00805 nBin[0],0.5,nBin[0]+0.5);
00806 nDim=0;
00807 }
00808 if(binMap) {
00809 *binMap=CreateBinMap(r,nDim,axisList,axisSteering);
00810 }
00811 return r;
00812 }
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828 TH2D *TUnfoldBinningV17::CreateHistogramOfMigrations
00829 (TUnfoldBinningV17 const *xAxis,TUnfoldBinningV17 const *yAxis,
00830 char const *histogramName,Bool_t originalXAxisBinning,
00831 Bool_t originalYAxisBinning,char const *histogramTitle)
00832 {
00833 Int_t nBinX[3],axisListX[3];
00834 Int_t nDimX=
00835 xAxis->GetTHxxBinning(originalXAxisBinning ? 1 : 0,nBinX,axisListX,0);
00836 const TUnfoldBinningV17 *neNodeX=xAxis->GetNonemptyNode();
00837 Int_t nBinY[3],axisListY[3];
00838 Int_t nDimY=
00839 yAxis->GetTHxxBinning(originalYAxisBinning ? 1 : 0,nBinY,axisListY,0);
00840 const TUnfoldBinningV17 *neNodeY=yAxis->GetNonemptyNode();
00841 TString title=xAxis->BuildHistogramTitle2D
00842 (histogramName,histogramTitle,axisListX[0],yAxis,axisListY[0]);
00843 if(nDimX==1) {
00844 const TVectorD *axisBinsX=(TVectorD const *)
00845 neNodeX->fAxisList->At(axisListX[0]);
00846 if(nDimY==1) {
00847 const TVectorD *axisBinsY=(TVectorD const *)
00848 neNodeY->fAxisList->At(axisListY[0]);
00849 return new TH2D(histogramName,title,
00850 nBinX[0],axisBinsX->GetMatrixArray(),
00851 nBinY[0],axisBinsY->GetMatrixArray());
00852 } else {
00853 return new TH2D(histogramName,title,
00854 nBinX[0],axisBinsX->GetMatrixArray(),
00855 nBinY[0],0.5,0.5+nBinY[0]);
00856 }
00857 } else {
00858 if(nDimY==1) {
00859 const TVectorD *axisBinsY=(TVectorD const *)
00860 neNodeY->fAxisList->At(axisListY[0]);
00861 return new TH2D(histogramName,title,
00862 nBinX[0],0.5,0.5+nBinX[0],
00863 nBinY[0],axisBinsY->GetMatrixArray());
00864 } else {
00865 return new TH2D(histogramName,title,
00866 nBinX[0],0.5,0.5+nBinX[0],
00867 nBinY[0],0.5,0.5+nBinY[0]);
00868 }
00869 }
00870 }
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885 Int_t TUnfoldBinningV17::GetTHxxBinning
00886 (Int_t maxDim,Int_t *axisBins,Int_t *axisList,
00887 const char *axisSteering) const
00888 {
00889 for(Int_t i=0;i<3;i++) {
00890 axisBins[i]=0;
00891 axisList[i]=-1;
00892 }
00893 const TUnfoldBinningV17 *theNode=GetNonemptyNode();
00894 if(theNode) {
00895 Int_t r=theNode->GetTHxxBinningSingleNode
00896 (maxDim,axisBins,axisList,axisSteering);
00897 return r;
00898 } else {
00899 axisBins[0]=GetTHxxBinsRecursive(axisSteering);
00900 return 0;
00901 }
00902 }
00903
00904
00905
00906
00907 const TUnfoldBinningV17 *TUnfoldBinningV17::GetNonemptyNode(void) const
00908 {
00909 const TUnfoldBinningV17 *r=GetDistributionNumberOfBins()>0 ? this : 0;
00910 for(TUnfoldBinningV17 const *child=GetChildNode();child;
00911 child=child->GetNextNode()) {
00912 const TUnfoldBinningV17 *c=child->GetNonemptyNode();
00913 if(!r) {
00914
00915 r=c;
00916 } else {
00917 if(c) {
00918
00919 r=0;
00920 break;
00921 }
00922 }
00923 }
00924 return r;
00925 }
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942 Int_t TUnfoldBinningV17::GetTHxxBinningSingleNode
00943 (Int_t maxDim,Int_t *axisBins,Int_t *axisList,const char *axisSteering) const
00944 {
00945
00946
00947
00948
00949 Int_t isOptionGiven[3];
00950 DecodeAxisSteering(axisSteering,"CUO",isOptionGiven);
00951
00952 Int_t numDimension=GetDistributionDimension();
00953 Int_t r=0;
00954 for(Int_t i=0;i<numDimension;i++) {
00955 if(isOptionGiven[0] & (1<<i)) continue;
00956 r++;
00957 }
00958 if((r>0)&&(r<=maxDim)) {
00959
00960
00961
00962
00963 r=0;
00964 for(Int_t i=0;i<numDimension;i++) {
00965 if(isOptionGiven[0] & (1<<i)) continue;
00966 axisList[r]=i;
00967 axisBins[r]=GetDistributionBinning(i)->GetNrows()-1;
00968 r++;
00969 }
00970 } else {
00971
00972
00973 if(HasUnconnectedBins() || (GetDistributionNumberOfBins()<=0)) {
00974 axisBins[0] = GetDistributionNumberOfBins();
00975 } else {
00976 Int_t nBin=1;
00977 for(Int_t i=0;i<numDimension;i++) {
00978 Int_t mask=(1<<i);
00979 if(isOptionGiven[0] & mask) continue;
00980 Int_t nBinI=GetDistributionBinning(i)->GetNrows()-1;
00981 if((fHasUnderflow & mask)&& !(isOptionGiven[1] & mask)) nBinI++;
00982 if((fHasOverflow & mask)&& !(isOptionGiven[2] & mask)) nBinI++;
00983 nBin *= nBinI;
00984 }
00985 axisBins[0] = nBin;
00986 }
00987 r=0;
00988 }
00989 return r;
00990 }
00991
00992
00993
00994
00995
00996
00997
00998
00999 Int_t TUnfoldBinningV17::GetTHxxBinsRecursive(const char *axisSteering) const
01000 {
01001
01002 Int_t r=0;
01003 for(TUnfoldBinningV17 const *child=GetChildNode();child;
01004 child=child->GetNextNode()) {
01005 r +=child->GetTHxxBinsRecursive(axisSteering);
01006 }
01007
01008 Int_t axisBins[3],axisList[3];
01009 GetTHxxBinningSingleNode(0,axisBins,axisList,axisSteering);
01010 r += axisBins[0];
01011 return r;
01012 }
01013
01014
01015
01016
01017
01018
01019 Int_t *TUnfoldBinningV17::CreateEmptyBinMap(void) const {
01020
01021
01022 Int_t nMax=GetRootNode()->GetEndBin()+1;
01023 Int_t *r=new Int_t[nMax];
01024 for(Int_t i=0;i<nMax;i++) {
01025 r[i]=-1;
01026 }
01027 return r;
01028 }
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038 void TUnfoldBinningV17::SetBinMapEntry
01039 (Int_t *binMap,Int_t globalBin,Int_t destBin) const {
01040 Int_t nMax=GetRootNode()->GetEndBin()+1;
01041 if((globalBin<0)||(globalBin>=nMax)) {
01042 Error("SetBinMapEntry","global bin number %d outside range (max=%d)",
01043 globalBin,nMax);
01044 } else {
01045 binMap[globalBin]=destBin;
01046 }
01047 }
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060 Int_t TUnfoldBinningV17::FillBinMap1D
01061 (Int_t *binMap,const char *axisSteering,Int_t firstBinX) const {
01062 Int_t r=firstBinX;
01063 Int_t axisBins[3],axisList[3];
01064 Int_t nDim=GetTHxxBinningSingleNode(3,axisBins,axisList,axisSteering);
01065 if((nDim==1)|| !GetDistributionDimension()) {
01066 r+=FillBinMapSingleNode(0,r,0,0,axisSteering,binMap);
01067 } else {
01068 Error("FillBinMap1D","distribution %s with steering=%s is not 1D",
01069 (char *)GetName(),axisSteering);
01070 }
01071 for(TUnfoldBinningV17 const *child=GetChildNode();child;
01072 child=child->GetNextNode()) {
01073 r =child->FillBinMap1D(binMap,axisSteering,r);
01074 }
01075 return r;
01076 }
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087 Int_t *TUnfoldBinningV17::CreateBinMap
01088 (const TH1 *hist,Int_t nDim,const Int_t *axisList,const char *axisSteering)
01089 const
01090 {
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124 Int_t *r=CreateEmptyBinMap();
01125 Int_t startBin=GetRootNode()->GetStartBin();
01126 if(nDim>0) {
01127 const TUnfoldBinning *nonemptyNode=GetNonemptyNode();
01128 if(nonemptyNode) {
01129 nonemptyNode->
01130 FillBinMapSingleNode(hist,startBin,nDim,axisList,axisSteering,r);
01131 } else {
01132 Fatal("CreateBinMap","called with nDim=%d but GetNonemptyNode()=0",
01133 nDim);
01134 }
01135 } else {
01136 FillBinMapRecursive(startBin,axisSteering,r);
01137 }
01138 return r;
01139 }
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151 Int_t TUnfoldBinningV17::FillBinMapRecursive
01152 (Int_t startBin,const char *axisSteering,Int_t *binMap) const
01153 {
01154 Int_t nbin=0;
01155 nbin = FillBinMapSingleNode(0,startBin,0,0,axisSteering,binMap);
01156 for(TUnfoldBinningV17 const *child=GetChildNode();child;
01157 child=child->GetNextNode()) {
01158 nbin += child->FillBinMapRecursive(startBin+nbin,axisSteering,binMap);
01159 }
01160 return nbin;
01161 }
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189 Int_t TUnfoldBinningV17::FillBinMapSingleNode
01190 (const TH1 *hist,Int_t startBin,Int_t nDim,const Int_t *axisList,
01191 const char *axisSteering,Int_t *binMap) const
01192 {
01193
01194
01195
01196
01197 Int_t isOptionGiven[3+10];
01198 DecodeAxisSteering(axisSteering,"CUO0123456789",isOptionGiven);
01199 Int_t haveSelectedBin=0;
01200 for(Int_t i=3;i<3+10;i++) {
01201 haveSelectedBin |= isOptionGiven[i];
01202 }
01203
01204 Int_t axisBins[MAXDIM];
01205 Int_t dimension=GetDistributionDimension();
01206 Int_t axisNbin[MAXDIM];
01207 for(Int_t i=0;i<dimension;i++) {
01208 const TVectorD *binning=GetDistributionBinning(i);
01209 axisNbin[i]=binning->GetNrows()-1;
01210 };
01211 for(Int_t i=0;i<GetDistributionNumberOfBins();i++) {
01212 Int_t globalBin=GetStartBin()+i;
01213 const TUnfoldBinningV17 *dest=ToAxisBins(globalBin,axisBins);
01214 if(dest!=this) {
01215 if(!dest) {
01216 Fatal("FillBinMapSingleNode",
01217 "bin %d outside binning scheme",
01218 globalBin);
01219 } else {
01220 Fatal("FillBinMapSingleNode",
01221 "bin %d located in %s %d-%d rather than %s %d=%d",
01222 i,(const char *)dest->GetName(),
01223 dest->GetStartBin(),dest->GetEndBin(),
01224 (const char *)GetName(),GetStartBin(),GetEndBin());
01225 }
01226 }
01227
01228 Bool_t skip=kFALSE;
01229 for(Int_t axis=0;axis<dimension;axis++) {
01230 Int_t mask=(1<<axis);
01231
01232 if(((axisBins[axis]<0)&&(isOptionGiven[1] & mask))||
01233 ((axisBins[axis]>=axisNbin[axis])&&(isOptionGiven[2] & mask)))
01234 skip=kTRUE;
01235
01236 if((axisBins[axis]>=0)&&(axisBins[axis]<axisNbin[axis])&&
01237 (haveSelectedBin & mask)) {
01238 if(!(isOptionGiven[3+axisBins[axis]] & mask)) skip=kTRUE;
01239 }
01240 }
01241 if(skip) {
01242 continue;
01243 }
01244
01245 if(nDim>0) {
01246
01247 if(nDim==hist->GetDimension()) {
01248 Int_t ibin[3];
01249 ibin[0]=ibin[1]=ibin[2]=0;
01250 for(Int_t hdim=0;hdim<nDim;hdim++) {
01251 Int_t axis=axisList[hdim];
01252 ibin[hdim]=axisBins[axis]+1;
01253 }
01254 binMap[globalBin]=hist->GetBin(ibin[0],ibin[1],ibin[2]);
01255 } else if(nDim==1) {
01256
01257
01258
01259
01260
01261 if((nDim!=1)||( hist->GetDimension()!=2)) {
01262
01263 Error("FillBinMapSingleNode","inconsistent dimensions %d %d",nDim,
01264 hist->GetDimension());
01265 }
01266 for(Int_t ii=0;ii<hist->GetDimension();ii++) {
01267 if(axisList[ii]>=0) {
01268 binMap[globalBin]=axisBins[axisList[ii]]+1;
01269 break;
01270 }
01271 }
01272 } else {
01273 Fatal("FillBinMapSingleNode","inconsistent dimensions %d %d",nDim,
01274 hist->GetDimension());
01275 }
01276 } else {
01277
01278
01279
01280
01281 if(dimension>0) {
01282 Int_t r=0;
01283 for(Int_t axis=dimension-1;axis>=0;axis--) {
01284 Int_t mask=(1<<axis);
01285 if(isOptionGiven[0] & mask) {
01286
01287 continue;
01288 }
01289 Int_t iBin=axisBins[axis];
01290 Int_t nMax=axisNbin[axis];
01291 if((fHasUnderflow & ~isOptionGiven[1]) & mask) {
01292 nMax +=1;
01293 iBin +=1;
01294 }
01295 if((fHasOverflow & ~isOptionGiven[2]) & mask) {
01296 nMax += 1;
01297 }
01298 r = r*nMax +iBin;
01299 }
01300 binMap[globalBin] = startBin + r;
01301 } else {
01302 binMap[globalBin] = startBin + axisBins[0];
01303 }
01304 }
01305 }
01306 Int_t nbin;
01307 if(dimension>0) {
01308 nbin=1;
01309 for(Int_t axis=dimension-1;axis>=0;axis--) {
01310 Int_t mask=(1<<axis);
01311 if(isOptionGiven[0] & mask) {
01312
01313 continue;
01314 }
01315 Int_t nMax=axisNbin[axis];
01316 if((fHasUnderflow & ~isOptionGiven[1]) & mask) {
01317 nMax +=1;
01318 }
01319 if((fHasOverflow & ~isOptionGiven[2]) & mask) {
01320 nMax += 1;
01321 }
01322 nbin = nbin*nMax;
01323 }
01324 } else {
01325 nbin=GetDistributionNumberOfBins();
01326 }
01327 return nbin;
01328 }
01329
01330
01331 TH1 *TUnfoldBinningV17::ExtractHistogram
01332 (const char *histogramName,const TH1 *globalBins,
01333 const TH2 *globalBinsEmatrix,Bool_t originalAxisBinning,
01334 const char *axisSteering) const
01335 {
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353 Int_t *binMap=0;
01354 TH1 *r=CreateHistogram(histogramName,originalAxisBinning,&binMap,0,
01355 axisSteering);
01356 if(!r) return 0;
01357 TUnfoldBinningV17 const *root=GetRootNode();
01358 Int_t nMax=-1;
01359 for(Int_t iSrc=root->GetStartBin();iSrc<root->GetEndBin();iSrc++) {
01360 if(binMap[iSrc]>nMax) nMax=binMap[iSrc];
01361 }
01362 if(nMax<0) {
01363 delete r;
01364 r=0;
01365 } else {
01366 TVectorD eSquared(nMax+1);
01367 for(Int_t iSrc=root->GetStartBin();iSrc<root->GetEndBin();iSrc++) {
01368 Int_t iDest=binMap[iSrc];
01369 if(iDest>=0) {
01370 Double_t c=r->GetBinContent(iDest);
01371 r->SetBinContent(iDest,c+globalBins->GetBinContent(iSrc));
01372 if(!globalBinsEmatrix) {
01373 eSquared(iDest)+=TMath::Power(globalBins->GetBinError(iSrc),2.);
01374 } else {
01375 for(Int_t jSrc=root->GetStartBin();jSrc<root->GetEndBin();
01376 jSrc++) {
01377 if(binMap[jSrc]==iDest) {
01378 eSquared(iDest) +=
01379 TMath::Power(globalBins->GetBinError(jSrc),2.);
01380 }
01381 }
01382 }
01383 }
01384 }
01385 for(Int_t i=0;i<nMax;i++) {
01386 Double_t e2=eSquared(i);
01387 if(e2>0.0) {
01388 r->SetBinError(i,TMath::Sqrt(e2));
01389 }
01390 }
01391 }
01392 delete binMap;
01393 return r;
01394 }
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406 Int_t TUnfoldBinningV17::GetGlobalBinNumber(Double_t x) const
01407 {
01408 if(GetDistributionDimension()!=1) {
01409 Fatal("GetBinNumber",
01410 "called with 1 argument for %d dimensional distribution",
01411 GetDistributionDimension());
01412 }
01413 return GetGlobalBinNumber(&x);
01414 }
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425 Int_t TUnfoldBinningV17::GetGlobalBinNumber(Double_t x,Double_t y) const
01426 {
01427 if(GetDistributionDimension()!=2) {
01428 Fatal("GetBinNumber",
01429 "called with 2 arguments for %d dimensional distribution",
01430 GetDistributionDimension());
01431 }
01432 Double_t xx[2];
01433 xx[0]=x;
01434 xx[1]=y;
01435 return GetGlobalBinNumber(xx);
01436 }
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448 Int_t TUnfoldBinningV17::GetGlobalBinNumber
01449 (Double_t x,Double_t y,Double_t z) const
01450 {
01451
01452
01453
01454 if(GetDistributionDimension()!=3) {
01455 Fatal("GetBinNumber",
01456 "called with 3 arguments for %d dimensional distribution",
01457 GetDistributionDimension());
01458 }
01459 Double_t xx[3];
01460 xx[0]=x;
01461 xx[1]=y;
01462 xx[2]=z;
01463 return GetGlobalBinNumber(xx);
01464 }
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477 Int_t TUnfoldBinningV17::GetGlobalBinNumber
01478 (Double_t x0,Double_t x1,Double_t x2,Double_t x3) const
01479 {
01480
01481
01482
01483 if(GetDistributionDimension()!=4) {
01484 Fatal("GetBinNumber",
01485 "called with 4 arguments for %d dimensional distribution",
01486 GetDistributionDimension());
01487 }
01488 Double_t xx[4];
01489 xx[0]=x0;
01490 xx[1]=x1;
01491 xx[2]=x2;
01492 xx[3]=x3;
01493 return GetGlobalBinNumber(xx);
01494 }
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508 Int_t TUnfoldBinningV17::GetGlobalBinNumber
01509 (Double_t x0,Double_t x1,Double_t x2,Double_t x3,Double_t x4) const
01510 {
01511
01512
01513
01514 if(GetDistributionDimension()!=5) {
01515 Fatal("GetBinNumber",
01516 "called with 5 arguments for %d dimensional distribution",
01517 GetDistributionDimension());
01518 }
01519 Double_t xx[5];
01520 xx[0]=x0;
01521 xx[1]=x1;
01522 xx[2]=x2;
01523 xx[3]=x3;
01524 xx[4]=x4;
01525 return GetGlobalBinNumber(xx);
01526 }
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541 Int_t TUnfoldBinningV17::GetGlobalBinNumber
01542 (Double_t x0,Double_t x1,Double_t x2,Double_t x3,Double_t x4,Double_t x5) const
01543 {
01544
01545
01546
01547 if(GetDistributionDimension()!=6) {
01548 Fatal("GetBinNumber",
01549 "called with 6 arguments for %d dimensional distribution",
01550 GetDistributionDimension());
01551 }
01552 Double_t xx[6];
01553 xx[0]=x0;
01554 xx[1]=x1;
01555 xx[2]=x2;
01556 xx[3]=x3;
01557 xx[4]=x4;
01558 xx[5]=x5;
01559 return GetGlobalBinNumber(xx);
01560 }
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584 Int_t TUnfoldBinningV17::GetGlobalBinNumber
01585 (const Double_t *x,Int_t *isBelow,Int_t *isAbove) const
01586 {
01587
01588
01589
01590
01591
01592
01593 if(!GetDistributionDimension()) {
01594 Fatal("GetBinNumber",
01595 "no axes are defined for node %s",
01596 (char const *)GetName());
01597 }
01598 Int_t iAxisBins[MAXDIM];
01599 for(Int_t dim=0;dim<GetDistributionDimension();dim++) {
01600 TVectorD const *bins=(TVectorD const *) fAxisList->At(dim);
01601 Int_t i0=0;
01602 Int_t i1=bins->GetNrows()-1;
01603 Int_t iBin= 0;
01604 if(!(x[dim]>=(*bins)[i0])) {
01605
01606 iBin += i0-1;
01607 } else if(!(x[dim]<(*bins)[i1])) {
01608
01609 iBin += i1;
01610 } else {
01611 while(i1-i0>1) {
01612 Int_t i2=(i0+i1)/2;
01613 if(x[dim]<(*bins)[i2]) {
01614 i1=i2;
01615 } else {
01616 i0=i2;
01617 }
01618 }
01619 iBin += i0;
01620 }
01621 iAxisBins[dim]=iBin;
01622 }
01623 Int_t r=ToGlobalBin(iAxisBins,isBelow,isAbove);
01624 if(r<0) r=0;
01625 return r;
01626 }
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636 TString TUnfoldBinningV17::GetBinName(Int_t iBin) const
01637 {
01638
01639
01640 Int_t axisBins[MAXDIM];
01641 TString r=TString::Format("#%d",iBin);
01642 TUnfoldBinningV17 const *distribution=ToAxisBins(iBin,axisBins);
01643 if(distribution) {
01644 r +=" (";
01645 r += distribution->GetName();
01646 Int_t dimension=distribution->GetDistributionDimension();
01647 if(dimension>0) {
01648 TString axisString;
01649 for(Int_t axis=0;axis<dimension;axis++) {
01650 TString thisAxisString=
01651 distribution->GetDistributionAxisLabel(axis);
01652 TVectorD const *bins=distribution->GetDistributionBinning(axis);
01653 Int_t i=axisBins[axis];
01654 if(i<0) thisAxisString += "[ufl]";
01655 else if(i>=bins->GetNrows()-1) thisAxisString += "[ofl]";
01656 else {
01657 thisAxisString +=
01658 TString::Format("[%.3g,%.3g]",(*bins)[i],(*bins)[i+1]);
01659 }
01660 axisString = ":"+thisAxisString+axisString;
01661 }
01662 r += axisString;
01663 } else {
01664
01665 Int_t i=axisBins[0];
01666 if((i>=0)&&(i<distribution->fAxisLabelList->GetEntriesFast())) {
01667 r += distribution->GetDistributionAxisLabel(i);
01668 } else {
01669 r += TString::Format(" %d",i);
01670 }
01671 }
01672 r +=")";
01673 }
01674 return r;
01675 }
01676
01677
01678
01679
01680
01681 Double_t TUnfoldBinningV17::GetBinSize(Int_t iBin) const
01682 {
01683
01684 Int_t axisBins[MAXDIM];
01685 TUnfoldBinningV17 const *distribution=ToAxisBins(iBin,axisBins);
01686 Double_t r=0.0;
01687 if(distribution) {
01688 if(distribution->GetDistributionDimension()>0) r=1.0;
01689 for(Int_t axis=0;axis<distribution->GetDistributionDimension();axis++) {
01690 TVectorD const *bins=distribution->GetDistributionBinning(axis);
01691 Int_t pos=axisBins[axis];
01692 if(pos<0) {
01693 r *= distribution->GetDistributionUnderflowBinWidth(axis);
01694 } else if(pos>=bins->GetNrows()-1) {
01695 r *= distribution->GetDistributionOverflowBinWidth(axis);
01696 } else {
01697 r *= (*bins)(pos+1)-(*bins)(pos);
01698 }
01699 if(r<=0.) break;
01700 }
01701 }
01702 return r;
01703 }
01704
01705
01706
01707 Bool_t TUnfoldBinningV17::IsBinFactorGlobal(void) const {
01708 return fBinFactorFunction ? kFALSE : kTRUE;
01709 }
01710
01711
01712
01713 Double_t TUnfoldBinningV17::GetGlobalFactor(void) const {
01714 return fBinFactorConstant;
01715 }
01716
01717
01718
01719
01720
01721
01722
01723
01724 Double_t TUnfoldBinningV17::GetBinFactor(Int_t iBin) const
01725 {
01726
01727
01728 Int_t axisBins[MAXDIM];
01729 TUnfoldBinningV17 const *distribution=ToAxisBins(iBin,axisBins);
01730 Double_t r=distribution->fBinFactorConstant;
01731 if((r!=0.0) && distribution->fBinFactorFunction) {
01732 TF1 *function=dynamic_cast<TF1 *>(distribution->fBinFactorFunction);
01733 if(function) {
01734 Double_t x[MAXDIM];
01735 Int_t dimension=distribution->GetDistributionDimension();
01736 if(dimension>0) {
01737 for(Int_t axis=0;axis<dimension;axis++) {
01738 x[axis]=distribution->GetDistributionBinCenter
01739 (axis,axisBins[axis]);
01740 }
01741 r *= function->EvalPar(x,function->GetParameters());
01742 } else {
01743 x[0]=axisBins[0];
01744 r *= function->Eval(x[0]);
01745 }
01746 } else {
01747 TVectorD *vect=dynamic_cast<TVectorD *>
01748 (distribution->fBinFactorFunction);
01749 if(vect) {
01750 r=(*vect)[iBin-GetStartBin()];
01751 } else {
01752 Error("GetBinFactor",
01753 "internal error: user function is neiter TF1 or TVectorD");
01754 }
01755 }
01756 }
01757 return r;
01758 }
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783 Int_t TUnfoldBinningV17::GetBinNeighbours
01784 (Int_t bin,Int_t axis,Int_t *prev,Double_t *distPrev,
01785 Int_t *next,Double_t *distNext,Bool_t isPeriodic) const
01786 {
01787 Int_t axisBins[MAXDIM];
01788 TUnfoldBinningV17 const *distribution=ToAxisBins(bin,axisBins);
01789 Int_t dimension=distribution->GetDistributionDimension();
01790 *prev=-1;
01791 *next=-1;
01792 *distPrev=0.;
01793 *distNext=0.;
01794 Int_t r=0;
01795 if((axis>=0)&&(axis<dimension)) {
01796
01797
01798 Int_t nMax=GetDistributionBinning(axis)->GetNrows()-1;
01799 Int_t centerBin= axisBins[axis];
01800 axisBins[axis] =centerBin-1;
01801 if(isPeriodic) {
01802 if(HasUnderflow(axis)) {
01803 r +=1;
01804 } else if((axisBins[axis]<0)&&(nMax>=3)) {
01805 axisBins[axis]=nMax-1;
01806 }
01807 }
01808 *prev=ToGlobalBin(axisBins);
01809 if(*prev>=0) {
01810 *distPrev=distribution->GetDistributionBinCenter(axis,axisBins[axis])-
01811 distribution->GetDistributionBinCenter(axis,centerBin);
01812 }
01813 axisBins[axis] =centerBin+1;
01814 if(isPeriodic) {
01815 if(HasOverflow(axis)) {
01816 r +=2;
01817 } else if((axisBins[axis]==nMax)&&(nMax>=3)) {
01818 axisBins[axis]=0;
01819 }
01820 }
01821 *next=ToGlobalBin(axisBins);
01822 if(*next>=0) {
01823 *distNext=distribution->GetDistributionBinCenter(axis,axisBins[axis])-
01824 distribution->GetDistributionBinCenter(axis,centerBin);
01825 }
01826 }
01827 return r;
01828 }
01829
01830
01831
01832
01833
01834
01835
01836 void TUnfoldBinningV17::GetBinUnderflowOverflowStatus
01837 (Int_t iBin,Int_t *uStatus,Int_t *oStatus) const
01838 {
01839 Int_t axisBins[MAXDIM];
01840 TUnfoldBinningV17 const *distribution=ToAxisBins(iBin,axisBins);
01841 Int_t dimension=distribution->GetDistributionDimension();
01842 *uStatus=0;
01843 *oStatus=0;
01844 for(Int_t axis=0;axis<dimension;axis++) {
01845 TVectorD const *bins=distribution->GetDistributionBinning(axis);
01846 Int_t nBin=bins->GetNrows()-1;
01847 if(axisBins[axis]<0) *uStatus |= (1<<axis);
01848 if(axisBins[axis]>=nBin) *oStatus |= (1<<axis);
01849 }
01850 }
01851
01852
01853
01854 Bool_t TUnfoldBinningV17::HasUnconnectedBins(void) const
01855 {
01856 return (!GetDistributionDimension())&&(GetDistributionNumberOfBins()>0);
01857 }
01858
01859
01860
01861
01862
01863 const TObjString *TUnfoldBinningV17::GetUnconnectedBinName(Int_t bin) const {
01864 TObjString *r=0;
01865 if(HasUnconnectedBins()) {
01866 if(bin<fAxisLabelList->GetEntriesFast()) {
01867 r=((TObjString * const)fAxisLabelList->At(bin));
01868 }
01869 }
01870 return r;
01871 }
01872
01873
01874
01875 Double_t TUnfoldBinningV17::GetDistributionAverageBinSize
01876
01877
01878
01879
01880
01881 (Int_t axis,Bool_t includeUnderflow,Bool_t includeOverflow) const
01882 {
01883 Double_t r=0.0;
01884 if((axis>=0)&&(axis<GetDistributionDimension())) {
01885 TVectorD const *bins=GetDistributionBinning(axis);
01886 Double_t d=(*bins)[bins->GetNrows()-1]-(*bins)[0];
01887 Double_t nBins=bins->GetNrows()-1;
01888 if(includeUnderflow && HasUnderflow(axis)) {
01889 Double_t w=GetDistributionUnderflowBinWidth(axis);
01890 if(w>0) {
01891 nBins++;
01892 d += w;
01893 }
01894 }
01895 if(includeOverflow && HasOverflow(axis)) {
01896 Double_t w=GetDistributionOverflowBinWidth(axis);
01897 if(w>0.0) {
01898 nBins++;
01899 d += w;
01900 }
01901 }
01902 if(nBins>0) {
01903 r=d/nBins;
01904 }
01905 } else {
01906 Error("GetDistributionAverageBinSize","axis %d does not exist",axis);
01907 }
01908 return r;
01909 }
01910
01911
01912
01913
01914
01915
01916
01917
01918 Double_t TUnfoldBinningV17::GetDistributionUnderflowBinWidth(Int_t axis) const
01919 {
01920
01921
01922 TVectorD const *bins=GetDistributionBinning(axis);
01923 return (*bins)[1]-(*bins)[0];
01924 }
01925
01926
01927
01928
01929
01930
01931
01932
01933 Double_t TUnfoldBinningV17::GetDistributionOverflowBinWidth(Int_t axis) const
01934 {
01935
01936
01937 TVectorD const *bins=GetDistributionBinning(axis);
01938 return (*bins)[bins->GetNrows()-1]-(*bins)[bins->GetNrows()-2];
01939 }
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951 Double_t TUnfoldBinningV17::GetDistributionBinCenter
01952 (Int_t axis,Int_t bin) const
01953 {
01954
01955
01956
01957
01958 TVectorD const *bins=GetDistributionBinning(axis);
01959 Double_t r=0.0;
01960 if(bin<0) {
01961
01962 r=(*bins)[0]-0.5*GetDistributionUnderflowBinWidth(axis);
01963 } else if(bin>=bins->GetNrows()-1) {
01964
01965 r=(*bins)[bins->GetNrows()-1]+0.5*GetDistributionOverflowBinWidth(axis);
01966 } else {
01967 r=0.5*((*bins)[bin+1]+(*bins)[bin]);
01968 }
01969 return r;
01970 }
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982 Int_t TUnfoldBinningV17::ToGlobalBin
01983 (Int_t const *axisBins,Int_t *isBelow,Int_t *isAbove) const
01984 {
01985 Int_t dimension=GetDistributionDimension();
01986 Int_t r=0;
01987 if(isBelow) *isBelow=0;
01988 if(isAbove) *isAbove=0;
01989 if(dimension>0) {
01990 for(Int_t axis=dimension-1;axis>=0;axis--) {
01991 Int_t nMax=GetDistributionBinning(axis)->GetNrows()-1;
01992 Int_t i=axisBins[axis];
01993 if(HasUnderflow(axis)) {
01994 nMax +=1;
01995 i +=1;
01996 }
01997 if(HasOverflow(axis)) nMax +=1;
01998 if((i>=0)&&(i<nMax)) {
01999 if(r>=0) r = r*nMax +i;
02000 } else {
02001 r=-1;
02002 if((i<0)&&(isBelow)) *isBelow |= 1<<axis;
02003 if((i>=nMax)&&(isAbove)) *isAbove |= 1<<axis;
02004 }
02005 }
02006 if(r>=0) {
02007 r += GetStartBin();
02008 }
02009 } else {
02010 if((axisBins[0]>=0)&&(axisBins[0]<GetDistributionNumberOfBins()))
02011 r=GetStartBin()+axisBins[0];
02012 else
02013 Fatal("ToGlobalBin","bad input %d for dimensionless binning %s %d",
02014 axisBins[0],(const char *)GetName(),
02015 GetDistributionNumberOfBins());
02016 }
02017 return r;
02018 }
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029 TUnfoldBinningV17 const *TUnfoldBinningV17::ToAxisBins
02030 (Int_t globalBin,Int_t *axisBins) const
02031 {
02032 TUnfoldBinningV17 const *r=0;
02033 if((globalBin>=GetStartBin())&&(globalBin<GetEndBin())) {
02034 TUnfoldBinningV17 const *node;
02035 for(node=GetChildNode();node && !r; node=node->GetNextNode()) {
02036 r=node->ToAxisBins(globalBin,axisBins);
02037 }
02038 if(!r) {
02039 r=this;
02040 Int_t i=globalBin-GetStartBin();
02041 Int_t dimension=GetDistributionDimension();
02042 if(dimension>0) {
02043 for(int axis=0;axis<dimension;axis++) {
02044 Int_t nMax=GetDistributionBinning(axis)->GetNrows()-1;
02045 axisBins[axis]=0;
02046 if(HasUnderflow(axis)) {
02047 axisBins[axis] =-1;
02048 nMax += 1;
02049 }
02050 if(HasOverflow(axis)) nMax +=1;
02051 axisBins[axis] += i % nMax;
02052 i /= nMax;
02053 }
02054 } else {
02055 axisBins[0]=i;
02056 }
02057 }
02058 }
02059 return r;
02060 }
02061
02062
02063
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080 void TUnfoldBinningV17::DecodeAxisSteering
02081 (const char *axisSteering,const char *options,Int_t *isOptionGiven) const
02082 {
02083 Int_t nOpt=TString(options).Length();
02084 for(Int_t i=0;i<nOpt;i++) isOptionGiven[i]=0;
02085 if(axisSteering) {
02086 TObjArray *patterns=TString(axisSteering).Tokenize(";");
02087 Int_t nPattern=patterns->GetEntries();
02088 Int_t nAxis=fAxisLabelList->GetEntries();
02089 for(Int_t i=0;i<nPattern;i++) {
02090 TString const &pattern=((TObjString * const)patterns->At(i))
02091 ->GetString();
02092 Int_t bracketBegin=pattern.Last('[');
02093 Int_t len=pattern.Length();
02094 if((bracketBegin>0)&&(pattern[len-1]==']')) {
02095 TString axisId=pattern(0,bracketBegin);
02096 Int_t mask=0;
02097 if((axisId[0]=='*')&&(axisId.Length()==1)) {
02098
02099 mask=(1<<nAxis)-1;
02100 } else {
02101
02102 for(Int_t j=0;j<nAxis;j++) {
02103 if(!axisId.CompareTo(GetDistributionAxisLabel(j))) {
02104 mask|= (1<<j);
02105 }
02106 }
02107 }
02108 for(Int_t o=0;o<nOpt;o++) {
02109 if(pattern.Last(options[o])>bracketBegin) {
02110 isOptionGiven[o] |= mask;
02111 }
02112 }
02113 } else {
02114 Error("DecodeAxisSteering",
02115 "steering \"%s\" does not end with [options]",
02116 (const char *)pattern);
02117 }
02118 }
02119 }
02120 }
02121