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
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152 #include "TUnfoldDensity.h"
00153 #include <TMath.h>
00154 #include <TVectorD.h>
00155 #include <TObjString.h>
00156 #include <iostream>
00157 #include <map>
00158
00159
00160
00161 #ifdef DEBUG
00162 using namespace std;
00163 #endif
00164
00165 ClassImp(TUnfoldDensityV17)
00166
00167 TUnfoldDensityV17::~TUnfoldDensityV17(void)
00168 {
00169
00170 if(fOwnedOutputBins) delete fOwnedOutputBins;
00171 if(fOwnedInputBins) delete fOwnedInputBins;
00172 if(fRegularisationConditions) delete fRegularisationConditions;
00173 }
00174
00175
00176
00177
00178 TUnfoldDensityV17::TUnfoldDensityV17(void)
00179 {
00180 fConstOutputBins=0;
00181 fConstInputBins=0;
00182 fOwnedOutputBins=0;
00183 fOwnedInputBins=0;
00184 fRegularisationConditions=0;
00185 }
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 TUnfoldDensityV17::TUnfoldDensityV17
00218 (const TH2 *hist_A, EHistMap histmap,ERegMode regmode,EConstraint constraint,
00219 EDensityMode densityMode,const TUnfoldBinningV17 *outputBins,
00220 const TUnfoldBinningV17 *inputBins,const char *regularisationDistribution,
00221 const char *regularisationAxisSteering) :
00222 TUnfoldSys(hist_A,histmap,kRegModeNone,constraint)
00223 {
00224 fRegularisationConditions=0;
00225
00226 fConstOutputBins = outputBins;
00227 fOwnedOutputBins = 0;
00228 TAxis const *genAxis,*detAxis;
00229 if(histmap==kHistMapOutputHoriz) {
00230 genAxis=hist_A->GetXaxis();
00231 detAxis=hist_A->GetYaxis();
00232 } else {
00233 genAxis=hist_A->GetYaxis();
00234 detAxis=hist_A->GetXaxis();
00235 }
00236 if(!fConstOutputBins) {
00237
00238
00239 fOwnedOutputBins =
00240 new TUnfoldBinningV17(*genAxis,1,1);
00241 fConstOutputBins = fOwnedOutputBins;
00242 }
00243
00244 if(fConstOutputBins->GetParentNode()) {
00245 Error("TUnfoldDensity",
00246 "Invalid output binning scheme (node is not the root node)");
00247 }
00248 fConstInputBins = inputBins;
00249 fOwnedInputBins = 0;
00250 if(!fConstInputBins) {
00251
00252
00253 fOwnedInputBins =
00254 new TUnfoldBinningV17(*detAxis,0,0);
00255 fConstInputBins = fOwnedInputBins;
00256 }
00257 if(fConstInputBins->GetParentNode()) {
00258 Error("TUnfoldDensity",
00259 "Invalid input binning scheme (node is not the root node)");
00260 }
00261
00262
00263 Int_t nOut=genAxis->GetNbins();
00264 Int_t nOutMappedT=TMath::Abs(fConstOutputBins->GetTH1xNumberOfBins(kTRUE));
00265 Int_t nOutMappedF=TMath::Abs(fConstOutputBins->GetTH1xNumberOfBins
00266 (fOwnedOutputBins));
00267 if((nOutMappedT!= nOut)&&(nOutMappedF!=nOut)) {
00268 Error("TUnfoldDensity",
00269 "Output binning incompatible number of bins: axis %d binning scheme %d (%d)",
00270 nOut,nOutMappedT,nOutMappedF);
00271 }
00272
00273 Int_t nInput=detAxis->GetNbins();
00274 Int_t nInputMappedT=TMath::Abs(fConstInputBins->GetTH1xNumberOfBins(kTRUE));
00275 Int_t nInputMappedF=TMath::Abs(fConstInputBins->GetTH1xNumberOfBins
00276 (fOwnedInputBins));
00277 if((nInputMappedT!= nInput)&&(nInputMappedF!= nInput)) {
00278 Error("TUnfoldDensity",
00279 "Input binning incompatible number of bins:axis %d binning scheme %d (%d) ",
00280 nInput,nInputMappedT,nInputMappedF);
00281 }
00282
00283
00284 for (Int_t ix = 0; ix <= nOut+1; ix++) {
00285 if(fHistToX[ix]<0) {
00286 Info("TUnfold","*NOT* unfolding bin "+GetOutputBinName(ix));
00287 }
00288 }
00289
00290
00291 if(regmode !=kRegModeNone) {
00292 RegularizeDistribution
00293 (regmode,densityMode,regularisationDistribution,
00294 regularisationAxisSteering);
00295 }
00296 }
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 TString TUnfoldDensityV17::GetOutputBinName(Int_t iBinX) const {
00308 if(!fConstOutputBins) return TUnfold::GetOutputBinName(iBinX);
00309 else return fConstOutputBins->GetBinName(iBinX);
00310 }
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331 Double_t TUnfoldDensityV17::GetDensityFactor
00332 (EDensityMode densityMode,Int_t iBin) const
00333 {
00334 Double_t factor=1.0;
00335 if((densityMode == kDensityModeBinWidth)||
00336 (densityMode == kDensityModeBinWidthAndUser)) {
00337 Double_t binSize=fConstOutputBins->GetBinSize(iBin);
00338 if(binSize>0.0) factor /= binSize;
00339 else factor=0.0;
00340 }
00341 if((densityMode == kDensityModeUser)||
00342 (densityMode == kDensityModeBinWidthAndUser)) {
00343 factor *= fConstOutputBins->GetBinFactor(iBin);
00344 }
00345 return factor;
00346 }
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382 void TUnfoldDensityV17::RegularizeDistribution
00383 (ERegMode regmode,EDensityMode densityMode,const char *distribution,
00384 const char *axisSteering)
00385 {
00386
00387 RegularizeDistributionRecursive(GetOutputBinning(),regmode,densityMode,
00388 distribution,axisSteering);
00389 }
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399 void TUnfoldDensityV17::RegularizeDistributionRecursive
00400 (const TUnfoldBinningV17 *binning,ERegMode regmode,
00401 EDensityMode densityMode,const char *distribution,const char *axisSteering) {
00402 if((!distribution)|| !TString(distribution).CompareTo(binning->GetName())) {
00403 RegularizeOneDistribution(binning,regmode,densityMode,axisSteering);
00404 }
00405 for(const TUnfoldBinningV17 *child=binning->GetChildNode();child;
00406 child=child->GetNextNode()) {
00407 RegularizeDistributionRecursive(child,regmode,densityMode,distribution,
00408 axisSteering);
00409 }
00410 }
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420 void TUnfoldDensityV17::RegularizeOneDistribution
00421 (const TUnfoldBinningV17 *binning,ERegMode regmode,
00422 EDensityMode densityMode,const char *axisSteering)
00423 {
00424 #ifdef DEBUG
00425 cout<<"TUnfoldDensityV17::RegularizeOneDistribution node="
00426 <<binning->GetName()<<" "<<regmode<<" "<<densityMode
00427 <<" "<<(axisSteering ? axisSteering : "")<<"\n";
00428 #endif
00429 if(!fRegularisationConditions)
00430 fRegularisationConditions=new TUnfoldBinning("regularisation");
00431
00432 TUnfoldBinning *thisRegularisationBinning=
00433 fRegularisationConditions->AddBinning(binning->GetName());
00434
00435
00436 Int_t isOptionGiven[8];
00437 binning->DecodeAxisSteering(axisSteering,"uUoObBpN",isOptionGiven);
00438
00439 isOptionGiven[0] |= isOptionGiven[1];
00440
00441 isOptionGiven[2] |= isOptionGiven[3];
00442
00443 isOptionGiven[4] |= isOptionGiven[5];
00444
00445 for(Int_t i=0;i<7;i++) {
00446 isOptionGiven[7] &= ~isOptionGiven[i];
00447 }
00448
00449 if(isOptionGiven[6] & (isOptionGiven[0]|isOptionGiven[2]) ) {
00450 Error("RegularizeOneDistribution",
00451 "axis steering %s is not valid",axisSteering);
00452 }
00453 #ifdef DEBUG
00454 cout<<" "<<isOptionGiven[0]
00455 <<" "<<isOptionGiven[1]
00456 <<" "<<isOptionGiven[2]
00457 <<" "<<isOptionGiven[3]
00458 <<" "<<isOptionGiven[4]
00459 <<" "<<isOptionGiven[5]
00460 <<" "<<isOptionGiven[6]
00461 <<" "<<isOptionGiven[7]
00462 <<"\n";
00463 #endif
00464 Info("RegularizeOneDistribution","regularizing %s regMode=%d"
00465 " densityMode=%d axisSteering=%s",
00466 binning->GetName(),(Int_t) regmode,(Int_t)densityMode,
00467 axisSteering ? axisSteering : "");
00468 Int_t startBin=binning->GetStartBin();
00469 Int_t endBin=startBin+ binning->GetDistributionNumberOfBins();
00470 std::vector<Double_t> factor(endBin-startBin);
00471 Int_t nbin=0;
00472 for(Int_t bin=startBin;bin<endBin;bin++) {
00473 factor[bin-startBin]=GetDensityFactor(densityMode,bin);
00474 if(factor[bin-startBin] !=0.0) nbin++;
00475 }
00476 #ifdef DEBUG
00477 cout<<"initial number of bins "<<nbin<<"\n";
00478 #endif
00479 Int_t dimension=binning->GetDistributionDimension();
00480
00481
00482 nbin=0;
00483 for(Int_t bin=startBin;bin<endBin;bin++) {
00484 Int_t uStatus,oStatus;
00485 binning->GetBinUnderflowOverflowStatus(bin,&uStatus,&oStatus);
00486 if(uStatus & isOptionGiven[1]) factor[bin-startBin]=0.;
00487 if(oStatus & isOptionGiven[3]) factor[bin-startBin]=0.;
00488 if(factor[bin-startBin] !=0.0) nbin++;
00489 }
00490 #ifdef DEBUG
00491 cout<<"after underflow/overflow bin removal "<<nbin<<"\n";
00492 #endif
00493 if(regmode==kRegModeSize) {
00494 Int_t nRegBins=0;
00495
00496
00497 for(Int_t bin=startBin;bin<endBin;bin++) {
00498 if(factor[bin-startBin]==0.0) continue;
00499 if(AddRegularisationCondition(bin,factor[bin-startBin])) {
00500 nRegBins++;
00501 }
00502 }
00503 if(nRegBins) {
00504 thisRegularisationBinning->AddBinning("size",nRegBins);
00505 }
00506 } else if((regmode==kRegModeDerivative)||(regmode==kRegModeCurvature)) {
00507 for(Int_t direction=0;direction<dimension;direction++) {
00508
00509 Int_t nRegBins=0;
00510 Int_t directionMask=(1<<direction);
00511 if(isOptionGiven[7] & directionMask) {
00512 #ifdef DEBUG
00513 cout<<"skip direction "<<direction<<"\n";
00514 #endif
00515 continue;
00516 }
00517 Double_t binDistanceNormalisation=
00518 (isOptionGiven[5] & directionMask) ?
00519 binning->GetDistributionAverageBinSize
00520 (direction,isOptionGiven[0] & directionMask,
00521 isOptionGiven[2] & directionMask) : 1.0;
00522 for(Int_t bin=startBin;bin<endBin;bin++) {
00523
00524 if(factor[bin-startBin]==0.0) continue;
00525
00526 Int_t iPrev,iNext;
00527 Double_t distPrev,distNext;
00528 Int_t error=binning->GetBinNeighbours
00529 (bin,direction,&iPrev,&distPrev,&iNext,&distNext,
00530 isOptionGiven[6] & directionMask);
00531 if(error) {
00532 Error("RegularizeOneDistribution",
00533 "invalid option %s (isPeriodic) for axis %s"
00534 " (has underflow or overflow)",axisSteering,
00535 binning->GetDistributionAxisLabel(direction).Data());
00536 }
00537 if((regmode==kRegModeDerivative)&&(iNext>=0)) {
00538 Double_t f0 = -factor[bin-startBin];
00539 Double_t f1 = factor[iNext-startBin];
00540 if(isOptionGiven[4] & directionMask) {
00541 if(distNext>0.0) {
00542 f0 *= binDistanceNormalisation/distNext;
00543 f1 *= binDistanceNormalisation/distNext;
00544 } else {
00545 f0=0.;
00546 f1=0.;
00547 }
00548 }
00549 if((f0==0.0)||(f1==0.0)) continue;
00550 if(AddRegularisationCondition(bin,f0,iNext,f1)) {
00551 nRegBins++;
00552 #ifdef DEBUG
00553 std::cout<<"Added Reg: bin "<<bin<<" "<<f0
00554 <<" next: "<<iNext<<" "<<f1<<"\n";
00555 #endif
00556 }
00557 } else if((regmode==kRegModeCurvature)&&(iPrev>=0)&&(iNext>=0)) {
00558 Double_t f0 = factor[iPrev-startBin];
00559 Double_t f1 = -factor[bin-startBin];
00560 Double_t f2 = factor[iNext-startBin];
00561 if(isOptionGiven[4] & directionMask) {
00562 if((distPrev<0.)&&(distNext>0.)) {
00563 distPrev= -distPrev;
00564 Double_t f=TMath::Power(binDistanceNormalisation,2.)/
00565 (distPrev+distNext);
00566 f0 *= f/distPrev;
00567 f1 *= f*(1./distPrev+1./distNext);
00568 f2 *= f/distNext;
00569 } else {
00570 f0=0.;
00571 f1=0.;
00572 f2=0.;
00573 }
00574 }
00575 if((f0==0.0)||(f1==0.0)||(f2==0.0)) continue;
00576 if(AddRegularisationCondition(iPrev,f0,bin,f1,iNext,f2)) {
00577 nRegBins++;
00578 #ifdef DEBUG
00579 std::cout<<"Added Reg: prev "<<iPrev<<" "<<f0
00580 <<" bin: "<<bin<<" "<<f1
00581 <<" next: "<<iNext<<" "<<f2<<"\n";
00582 #endif
00583 }
00584 }
00585 }
00586 if(nRegBins) {
00587 TString name;
00588 if(regmode==kRegModeDerivative) {
00589 name="derivative_";
00590 } else if(regmode==kRegModeCurvature) {
00591 name="curvature_";
00592 }
00593 name += binning->GetDistributionAxisLabel(direction);
00594 thisRegularisationBinning->AddBinning(name,nRegBins);
00595 }
00596 }
00597 }
00598 #ifdef DEBUG
00599
00600 #endif
00601 }
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650 TH1 *TUnfoldDensityV17::GetOutput
00651 (const char *histogramName,const char *histogramTitle,
00652 const char *distributionName,const char *axisSteering,
00653 Bool_t useAxisBinning) const
00654 {
00655 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
00656 Int_t *binMap=0;
00657 TH1 *r=binning->CreateHistogram
00658 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00659 if(r) {
00660 TUnfoldSys::GetOutput(r,binMap);
00661 }
00662 if(binMap) {
00663 delete [] binMap;
00664 }
00665 return r;
00666 }
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681 TH1 *TUnfoldDensityV17::GetBias
00682 (const char *histogramName,const char *histogramTitle,
00683 const char *distributionName,const char *axisSteering,
00684 Bool_t useAxisBinning) const
00685 {
00686 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
00687 Int_t *binMap=0;
00688 TH1 *r=binning->CreateHistogram
00689 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00690 if(r) {
00691 TUnfoldSys::GetBias(r,binMap);
00692 }
00693 if(binMap) delete [] binMap;
00694 return r;
00695 }
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712 TH1 *TUnfoldDensityV17::GetFoldedOutput
00713 (const char *histogramName,const char *histogramTitle,
00714 const char *distributionName,const char *axisSteering,
00715 Bool_t useAxisBinning,Bool_t addBgr) const
00716 {
00717 TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
00718 Int_t *binMap=0;
00719 TH1 *r=binning->CreateHistogram
00720 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00721 if(r) {
00722 TUnfoldSys::GetFoldedOutput(r,binMap);
00723 if(addBgr) {
00724 TUnfoldSys::GetBackground(r,0,binMap,0,kFALSE);
00725 }
00726 }
00727 if(binMap) delete [] binMap;
00728 return r;
00729 }
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747 TH1 *TUnfoldDensityV17::GetBackground
00748 (const char *histogramName,const char *bgrSource,const char *histogramTitle,
00749 const char *distributionName,const char *axisSteering,Bool_t useAxisBinning,
00750 Int_t includeError) const
00751 {
00752 TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
00753 Int_t *binMap=0;
00754 TH1 *r=binning->CreateHistogram
00755 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00756 if(r) {
00757 TUnfoldSys::GetBackground(r,bgrSource,binMap,includeError,kTRUE);
00758 }
00759 if(binMap) delete [] binMap;
00760 return r;
00761 }
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777 TH1 *TUnfoldDensityV17::GetInput
00778 (const char *histogramName,const char *histogramTitle,
00779 const char *distributionName,const char *axisSteering,
00780 Bool_t useAxisBinning) const
00781 {
00782 TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
00783 Int_t *binMap=0;
00784 TH1 *r=binning->CreateHistogram
00785 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00786 if(r) {
00787 TUnfoldSys::GetInput(r,binMap);
00788 }
00789 if(binMap) delete [] binMap;
00790 return r;
00791 }
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809 TH1 *TUnfoldDensityV17::GetRhoItotal
00810 (const char *histogramName,const char *histogramTitle,
00811 const char *distributionName,const char *axisSteering,
00812 Bool_t useAxisBinning,TH2 **ematInv) {
00813 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
00814 Int_t *binMap=0;
00815 TH1 *r=binning->CreateHistogram
00816 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00817 if(r) {
00818 TH2 *invEmat=0;
00819 if(ematInv) {
00820 if(r->GetDimension()==1) {
00821 TString ematName(histogramName);
00822 ematName += "_inverseEMAT";
00823 Int_t *binMap2D=0;
00824 invEmat=binning->CreateErrorMatrixHistogram
00825 (ematName,useAxisBinning,&binMap2D,histogramTitle,
00826 axisSteering);
00827 if(binMap2D) delete [] binMap2D;
00828 } else {
00829 Error("GetRhoItotal",
00830 "can not return inverse of error matrix for this binning");
00831 }
00832 }
00833 TUnfoldSys::GetRhoItotal(r,binMap,invEmat);
00834 if(invEmat) {
00835 *ematInv=invEmat;
00836 }
00837 }
00838 if(binMap) delete [] binMap;
00839 return r;
00840 }
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859 TH1 *TUnfoldDensityV17::GetRhoIstatbgr
00860 (const char *histogramName,const char *histogramTitle,
00861 const char *distributionName,const char *axisSteering,
00862 Bool_t useAxisBinning,TH2 **ematInv) {
00863 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
00864 Int_t *binMap=0;
00865 TH1 *r=binning->CreateHistogram
00866 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00867 if(r) {
00868 TH2 *invEmat=0;
00869 if(ematInv) {
00870 if(r->GetDimension()==1) {
00871 TString ematName(histogramName);
00872 ematName += "_inverseEMAT";
00873 Int_t *binMap2D=0;
00874 invEmat=binning->CreateErrorMatrixHistogram
00875 (ematName,useAxisBinning,&binMap2D,histogramTitle,
00876 axisSteering);
00877 if(binMap2D) delete [] binMap2D;
00878 } else {
00879 Error("GetRhoItotal",
00880 "can not return inverse of error matrix for this binning");
00881 }
00882 }
00883 TUnfoldSys::GetRhoI(r,binMap,invEmat);
00884 if(invEmat) {
00885 *ematInv=invEmat;
00886 }
00887 }
00888 if(binMap) delete [] binMap;
00889 return r;
00890 }
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907 TH1 *TUnfoldDensityV17::GetDeltaSysSource
00908 (const char *source,const char *histogramName,
00909 const char *histogramTitle,const char *distributionName,
00910 const char *axisSteering,Bool_t useAxisBinning) {
00911 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
00912 Int_t *binMap=0;
00913 TH1 *r=binning->CreateHistogram
00914 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00915 if(r) {
00916 if(!TUnfoldSys::GetDeltaSysSource(r,source,binMap)) {
00917 delete r;
00918 r=0;
00919 }
00920 }
00921 if(binMap) delete [] binMap;
00922 return r;
00923 }
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940 TH1 *TUnfoldDensityV17::GetDeltaSysBackgroundScale
00941 (const char *bgrSource,const char *histogramName,
00942 const char *histogramTitle,const char *distributionName,
00943 const char *axisSteering,Bool_t useAxisBinning) {
00944 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
00945 Int_t *binMap=0;
00946 TH1 *r=binning->CreateHistogram
00947 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00948 if(r) {
00949 if(!TUnfoldSys::GetDeltaSysBackgroundScale(r,bgrSource,binMap)) {
00950 delete r;
00951 r=0;
00952 }
00953 }
00954 if(binMap) delete [] binMap;
00955 return r;
00956 }
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972 TH1 *TUnfoldDensityV17::GetDeltaSysTau
00973 (const char *histogramName,const char *histogramTitle,
00974 const char *distributionName,const char *axisSteering,Bool_t useAxisBinning)
00975 {
00976 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
00977 Int_t *binMap=0;
00978 TH1 *r=binning->CreateHistogram
00979 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
00980 if(r) {
00981 if(!TUnfoldSys::GetDeltaSysTau(r,binMap)) {
00982 delete r;
00983 r=0;
00984 }
00985 }
00986 if(binMap) delete [] binMap;
00987 return r;
00988 }
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003 TH2 *TUnfoldDensityV17::GetRhoIJtotal
01004 (const char *histogramName,const char *histogramTitle,
01005 const char *distributionName,const char *axisSteering,
01006 Bool_t useAxisBinning)
01007 {
01008 TH2 *r=GetEmatrixTotal
01009 (histogramName,histogramTitle,distributionName,
01010 axisSteering,useAxisBinning);
01011 if(r) {
01012 for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
01013 Double_t e_i=r->GetBinContent(i,i);
01014 if(e_i>0.0) e_i=TMath::Sqrt(e_i);
01015 else e_i=0.0;
01016 for(Int_t j=0;j<=r->GetNbinsY()+1;j++) {
01017 if(i==j) continue;
01018 Double_t e_j=r->GetBinContent(j,j);
01019 if(e_j>0.0) e_j=TMath::Sqrt(e_j);
01020 else e_j=0.0;
01021 Double_t e_ij=r->GetBinContent(i,j);
01022 if((e_i>0.0)&&(e_j>0.0)) {
01023 r->SetBinContent(i,j,e_ij/e_i/e_j);
01024 } else {
01025 r->SetBinContent(i,j,0.0);
01026 }
01027 }
01028 }
01029 for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
01030 if(r->GetBinContent(i,i)>0.0) {
01031 r->SetBinContent(i,i,1.0);
01032 } else {
01033 r->SetBinContent(i,i,0.0);
01034 }
01035 }
01036 }
01037 return r;
01038 }
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054 TH2 *TUnfoldDensityV17::GetEmatrixSysUncorr
01055 (const char *histogramName,const char *histogramTitle,
01056 const char *distributionName,const char *axisSteering,
01057 Bool_t useAxisBinning)
01058 {
01059 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
01060 Int_t *binMap=0;
01061 TH2 *r=binning->CreateErrorMatrixHistogram
01062 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
01063 if(r) {
01064 TUnfoldSys::GetEmatrixSysUncorr(r,binMap);
01065 }
01066 if(binMap) delete [] binMap;
01067 return r;
01068 }
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084 TH2 *TUnfoldDensityV17::GetEmatrixSysBackgroundUncorr
01085 (const char *bgrSource,const char *histogramName,
01086 const char *histogramTitle,const char *distributionName,
01087 const char *axisSteering,Bool_t useAxisBinning)
01088 {
01089 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
01090 Int_t *binMap=0;
01091 TH2 *r=binning->CreateErrorMatrixHistogram
01092 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
01093 if(r) {
01094 TUnfoldSys::GetEmatrixSysBackgroundUncorr(r,bgrSource,binMap,kFALSE);
01095 }
01096 if(binMap) delete [] binMap;
01097 return r;
01098 }
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115 TH2 *TUnfoldDensityV17::GetEmatrixInput
01116 (const char *histogramName,const char *histogramTitle,
01117 const char *distributionName,const char *axisSteering,
01118 Bool_t useAxisBinning)
01119 {
01120 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
01121 Int_t *binMap=0;
01122 TH2 *r=binning->CreateErrorMatrixHistogram
01123 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
01124 if(r) {
01125 TUnfoldSys::GetEmatrixInput(r,binMap);
01126 }
01127 if(binMap) delete [] binMap;
01128 return r;
01129 }
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141 TH2 *TUnfoldDensityV17::GetProbabilityMatrix
01142 (const char *histogramName,const char *histogramTitle,
01143 Bool_t useAxisBinning) const
01144 {
01145 TH2 *r=TUnfoldBinning::CreateHistogramOfMigrations
01146 (fConstOutputBins,fConstInputBins,histogramName,
01147 useAxisBinning,useAxisBinning,histogramTitle);
01148 TUnfoldV17::GetProbabilityMatrix(r,kHistMapOutputHoriz);
01149 return r;
01150 }
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166 TH2 *TUnfoldDensityV17::GetEmatrixTotal
01167 (const char *histogramName,const char *histogramTitle,
01168 const char *distributionName,const char *axisSteering,
01169 Bool_t useAxisBinning)
01170 {
01171 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
01172 Int_t *binMap=0;
01173 TH2 *r=binning->CreateErrorMatrixHistogram
01174 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
01175 if(r) {
01176 TUnfoldSys::GetEmatrixTotal(r,binMap);
01177 }
01178 if(binMap) delete [] binMap;
01179 return r;
01180 }
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193 TH2 *TUnfoldDensityV17::GetL
01194 (const char *histogramName,const char *histogramTitle,Bool_t useAxisBinning)
01195 {
01196 if(fRegularisationConditions &&
01197 (fRegularisationConditions->GetEndBin()-
01198 fRegularisationConditions->GetStartBin()!= fL->GetNrows())) {
01199 Warning("GetL",
01200 "remove invalid scheme of regularisation conditions %d %d",
01201 fRegularisationConditions->GetEndBin(),fL->GetNrows());
01202 delete fRegularisationConditions;
01203 fRegularisationConditions=0;
01204 }
01205 if(!fRegularisationConditions) {
01206 fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
01207 Warning("GetL","create flat regularisation conditions scheme");
01208 }
01209 TH2 *r=TUnfoldBinning::CreateHistogramOfMigrations
01210 (fConstOutputBins,fRegularisationConditions,histogramName,
01211 useAxisBinning,useAxisBinning,histogramTitle);
01212 TUnfold::GetL(r);
01213 return r;
01214 }
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230 TH1 *TUnfoldDensityV17::GetLxMinusBias
01231 (const char *histogramName,const char *histogramTitle)
01232 {
01233 TMatrixD dx(*GetX(), TMatrixD::kMinus, fBiasScale * (*fX0));
01234 TMatrixDSparse *Ldx=MultiplyMSparseM(fL,&dx);
01235 if(fRegularisationConditions &&
01236 (fRegularisationConditions->GetEndBin()-
01237 fRegularisationConditions->GetStartBin()!= fL->GetNrows())) {
01238 Warning("GetLxMinusBias",
01239 "remove invalid scheme of regularisation conditions %d %d",
01240 fRegularisationConditions->GetEndBin(),fL->GetNrows());
01241 delete fRegularisationConditions;
01242 fRegularisationConditions=0;
01243 }
01244 if(!fRegularisationConditions) {
01245 fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
01246 Warning("GetLxMinusBias","create flat regularisation conditions scheme");
01247 }
01248 TH1 *r=fRegularisationConditions->CreateHistogram
01249 (histogramName,kFALSE,0,histogramTitle);
01250 const Int_t *Ldx_rows=Ldx->GetRowIndexArray();
01251 const Double_t *Ldx_data=Ldx->GetMatrixArray();
01252 for(Int_t row=0;row<Ldx->GetNrows();row++) {
01253 if(Ldx_rows[row]<Ldx_rows[row+1]) {
01254 r->SetBinContent(row+1,Ldx_data[Ldx_rows[row]]);
01255 }
01256 }
01257 delete Ldx;
01258 return r;
01259 }
01260
01261
01262
01263
01264
01265
01266
01267
01268 const TUnfoldBinning *TUnfoldDensityV17::GetInputBinning
01269 (const char *distributionName) const
01270 {
01271
01272
01273 return fConstInputBins->FindNode(distributionName);
01274 }
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284 const TUnfoldBinning *TUnfoldDensityV17::GetOutputBinning
01285 (const char *distributionName) const
01286 {
01287
01288
01289 return fConstOutputBins->FindNode(distributionName);
01290 }
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329 Int_t TUnfoldDensityV17::ScanTau
01330 (Int_t nPoint,Double_t tauMin,Double_t tauMax,TSpline **scanResult,
01331 Int_t mode,const char *distribution,const char *axisSteering,
01332 TGraph **lCurvePlot,TSpline **logTauXPlot,TSpline **logTauYPlot)
01333 {
01334 typedef std::map<Double_t,Double_t> TauScan_t;
01335 typedef std::map<Double_t,std::pair<Double_t,Double_t> > LCurve_t;
01336 TauScan_t curve;
01337 LCurve_t lcurve;
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360 if((tauMin<=0)||(tauMax<=0.0)||(tauMin>=tauMax)) {
01361
01362
01363
01364
01365
01366 DoUnfold(0.0);
01367
01368
01369 if(GetNdf()<=0) {
01370 Error("ScanTau","too few input bins, NDF<=0 %d",GetNdf());
01371 }
01372 Double_t X0=GetLcurveX();
01373 Double_t Y0=GetLcurveY();
01374 Double_t y0=GetScanVariable(mode,distribution,axisSteering);
01375 Info("ScanTau","logtau=-Infinity y=%lf X=%lf Y=%lf",y0,X0,Y0);
01376 {
01377
01378 Double_t logTau=
01379 0.5*(TMath::Log10(GetChi2A()+3.*TMath::Sqrt(GetNdf()+1.0))
01380 -GetLcurveY());
01381 DoUnfold(TMath::Power(10.,logTau));
01382 if((!finite(GetLcurveX())) ||(!finite(GetLcurveY()))) {
01383 Fatal("ScanTau","problem (missing regularisation?) X=%f Y=%f",
01384 GetLcurveX(),GetLcurveY());
01385 }
01386 Double_t y=GetScanVariable(mode,distribution,axisSteering);
01387 curve[logTau]=y;
01388 lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
01389 Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
01390 GetLcurveX(),GetLcurveY());
01391 }
01392
01393
01394
01395 while(((int)curve.size()<nPoint-1)&&
01396 ((TMath::Abs(GetLcurveX()-X0)>0.00432)||
01397 (TMath::Abs(GetLcurveY()-Y0)>0.00432))) {
01398 Double_t logTau=(*curve.begin()).first-0.5;
01399 DoUnfold(TMath::Power(10.,logTau));
01400 Double_t y=GetScanVariable(mode,distribution,axisSteering);
01401 curve[logTau]=y;
01402 lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
01403 Info("ScanTay","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
01404 GetLcurveX(),GetLcurveY());
01405 }
01406 } else {
01407 Double_t logTauMin=TMath::Log10(tauMin);
01408 Double_t logTauMax=TMath::Log10(tauMax);
01409 if(nPoint>1) {
01410
01411 DoUnfold(TMath::Power(10.,logTauMax));
01412 Double_t y=GetScanVariable(mode,distribution,axisSteering);
01413 curve[logTauMax]=y;
01414 lcurve[logTauMax]=std::make_pair(GetLcurveX(),GetLcurveY());
01415 Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMax,y,
01416 GetLcurveX(),GetLcurveY());
01417 }
01418
01419 DoUnfold(TMath::Power(10.,logTauMin));
01420 Double_t y=GetScanVariable(mode,distribution,axisSteering);
01421 curve[logTauMin]=y;
01422 lcurve[logTauMin]=std::make_pair(GetLcurveX(),GetLcurveY());
01423 Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMin,y,
01424 GetLcurveX(),GetLcurveY());
01425 }
01426
01427
01428
01429
01430 while((int)curve.size()<nPoint-1) {
01431
01432
01433
01434
01435
01436 TauScan_t::const_iterator i0,i1;
01437 i0=curve.begin();
01438
01439 Double_t logTauYMin=(*i0).first;
01440 Double_t yMin=(*i0).second;
01441 for(;i0!=curve.end();i0++) {
01442 if((*i0).second<yMin) {
01443 yMin=(*i0).second;
01444 logTauYMin=(*i0).first;
01445 }
01446 }
01447
01448
01449 i0=curve.begin();
01450 i1=i0;
01451 Double_t distMax=0.0;
01452 Double_t logTau=0.0;
01453 for(i1++;i1!=curve.end();i1++) {
01454 Double_t dist;
01455
01456 dist=TMath::Abs((*i0).first-(*i1).first)
01457
01458 +0.25*TMath::Power(0.5*((*i0).first+(*i1).first)-logTauYMin,2.)/
01459 ((*curve.rbegin()).first-(*curve.begin()).first)/nPoint;
01460 if((dist<=0.0)||(dist>distMax)) {
01461 distMax=dist;
01462 logTau=0.5*((*i0).first+(*i1).first);
01463 }
01464 i0=i1;
01465 }
01466 DoUnfold(TMath::Power(10.,logTau));
01467 Double_t y=GetScanVariable(mode,distribution,axisSteering);
01468 curve[logTau]=y;
01469 lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
01470 Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
01471 GetLcurveX(),GetLcurveY());
01472 }
01473
01474
01475
01476
01477
01478
01479
01480 Double_t cTmin=0.0;
01481 {
01482 Double_t *cTi=new Double_t[curve.size()];
01483 Double_t *cCi=new Double_t[curve.size()];
01484 Int_t n=0;
01485 for(TauScan_t::const_iterator i=curve.begin();i!=curve.end();i++) {
01486 cTi[n]=(*i).first;
01487 cCi[n]=(*i).second;
01488 n++;
01489 }
01490
01491 TSpline3 *splineC=new TSpline3("L curve curvature",cTi,cCi,n);
01492
01493
01494
01495
01496
01497 Int_t iskip=0;
01498 if(n>3) iskip=1;
01499 if(n>6) iskip=2;
01500 Double_t cCmin=cCi[iskip];
01501 cTmin=cTi[iskip];
01502 for(Int_t i=iskip;i<n-1-iskip;i++) {
01503
01504
01505 Double_t xMin=cTi[i+1];
01506 Double_t yMin=cCi[i+1];
01507 if(cCi[i]<yMin) {
01508 yMin=cCi[i];
01509 xMin=cTi[i];
01510 }
01511
01512
01513
01514 Double_t x,y,b,c,d;
01515 splineC->GetCoeff(i,x,y,b,c,d);
01516
01517 Double_t m_p_half=-c/(3.*d);
01518 Double_t q=b/(3.*d);
01519 Double_t discr=m_p_half*m_p_half-q;
01520 if(discr>=0.0) {
01521
01522 discr=TMath::Sqrt(discr);
01523 Double_t xx;
01524 if(m_p_half>0.0) {
01525 xx = m_p_half + discr;
01526 } else {
01527 xx = m_p_half - discr;
01528 }
01529 Double_t dx=cTi[i+1]-x;
01530
01531 if((xx>0.0)&&(xx<dx)) {
01532 y=splineC->Eval(x+xx);
01533 if(y<yMin) {
01534 yMin=y;
01535 xMin=x+xx;
01536 }
01537 }
01538
01539 if(xx !=0.0) {
01540 xx= q/xx;
01541 } else {
01542 xx=0.0;
01543 }
01544
01545 if((xx>0.0)&&(xx<dx)) {
01546 y=splineC->Eval(x+xx);
01547 if(y<yMin) {
01548 yMin=y;
01549 xMin=x+xx;
01550 }
01551 }
01552 }
01553
01554 if(yMin<cCmin) {
01555 cCmin=yMin;
01556 cTmin=xMin;
01557 }
01558 }
01559 delete splineC;
01560 delete[] cTi;
01561 delete[] cCi;
01562 }
01563 Double_t logTauFin=cTmin;
01564 DoUnfold(TMath::Power(10.,logTauFin));
01565 {
01566 Double_t y=GetScanVariable(mode,distribution,axisSteering);
01567 curve[logTauFin]=y;
01568 lcurve[logTauFin]=std::make_pair(GetLcurveX(),GetLcurveY());
01569 Info("ScanTau","Result logtau=%lf y=%lf X=%lf Y=%lf",logTauFin,y,
01570 GetLcurveX(),GetLcurveY());
01571 }
01572
01573
01574
01575
01576 Int_t bestChoice=-1;
01577 if(curve.size()>0) {
01578 Double_t *y=new Double_t[curve.size()];
01579 Double_t *logT=new Double_t[curve.size()];
01580 int n=0;
01581 for( TauScan_t::const_iterator i=curve.begin();i!=curve.end();i++) {
01582 if(logTauFin==(*i).first) {
01583 bestChoice=n;
01584 }
01585 y[n]=(*i).second;
01586 logT[n]=(*i).first;
01587 n++;
01588 }
01589 if(scanResult) {
01590 TString name;
01591 name = TString::Format("scan(%d,",mode);
01592 if(distribution) name+= distribution;
01593 name += ",";
01594 if(axisSteering) name += axisSteering;
01595 name +=")";
01596 (*scanResult)=new TSpline3(name+"%log(tau)",logT,y,n);
01597 }
01598 delete[] y;
01599 delete[] logT;
01600 }
01601 if(lcurve.size()) {
01602 Double_t *logT=new Double_t[lcurve.size()];
01603 Double_t *x=new Double_t[lcurve.size()];
01604 Double_t *y=new Double_t[lcurve.size()];
01605 Int_t n=0;
01606 for(LCurve_t::const_iterator i=lcurve.begin();i!=lcurve.end();i++) {
01607 logT[n]=(*i).first;
01608 x[n]=(*i).second.first;
01609 y[n]=(*i).second.second;
01610
01611 n++;
01612 }
01613 if(lCurvePlot) {
01614 *lCurvePlot=new TGraph(n,x,y);
01615 (*lCurvePlot)->SetTitle("L curve");
01616 }
01617 if(logTauXPlot)
01618 *logTauXPlot=new TSpline3("log(chi**2)%log(tau)",logT,x,n);
01619 if(logTauYPlot)
01620 *logTauYPlot=new TSpline3("log(reg.cond)%log(tau)",logT,y,n);
01621 delete [] y;
01622 delete [] x;
01623 delete [] logT;
01624 }
01625 return bestChoice;
01626 }
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655 Double_t TUnfoldDensityV17::GetScanVariable
01656 (Int_t mode,const char *distribution,const char *axisSteering)
01657 {
01658 Double_t r=0.0;
01659 TString name("GetScanVariable(");
01660 name += TString::Format("%d,",mode);
01661 if(distribution) name += distribution;
01662 name += ",";
01663 if(axisSteering) name += axisSteering;
01664 name += ")";
01665 TH1 *rhoi=0;
01666 if((mode==kEScanTauRhoAvg)||(mode==kEScanTauRhoMax)||
01667 (mode==kEScanTauRhoSquareAvg)) {
01668 rhoi=GetRhoIstatbgr(name,0,distribution,axisSteering,kFALSE);
01669 } else if((mode==kEScanTauRhoAvgSys)||(mode==kEScanTauRhoMaxSys)||
01670 (mode==kEScanTauRhoSquareAvgSys)) {
01671 rhoi=GetRhoItotal(name,0,distribution,axisSteering,kFALSE);
01672 }
01673 if(rhoi) {
01674 Double_t sum=0.0;
01675 Double_t sumSquare=0.0;
01676 Double_t rhoMax=0.0;
01677 Int_t n=0;
01678 for(Int_t i=0;i<=rhoi->GetNbinsX()+1;i++) {
01679 Double_t c=rhoi->GetBinContent(i);
01680 if(c>=0.) {
01681 if(c>rhoMax) rhoMax=c;
01682 sum += c;
01683 sumSquare += c*c;
01684 n ++;
01685 }
01686 }
01687 if((mode==kEScanTauRhoAvg)||(mode==kEScanTauRhoAvgSys)) {
01688 r=sum/n;
01689 } else if((mode==kEScanTauRhoSquareAvg)||
01690 (mode==kEScanTauRhoSquareAvgSys)) {
01691 r=sum/n;
01692 } else {
01693 r=rhoMax;
01694 }
01695
01696 delete rhoi;
01697 } else {
01698 Fatal("GetScanVariable","mode %d not implemented",mode);
01699 }
01700 return r;
01701 }