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 #include <iostream>
00151 #include <TMap.h>
00152 #include <TMath.h>
00153 #include <TObjString.h>
00154 #include <TSortedList.h>
00155 #include <RVersion.h>
00156 #include <cmath>
00157
00158 #include "TUnfoldSys.h"
00159
00160 ClassImp(TUnfoldSysV17)
00161
00162 TUnfoldSysV17::~TUnfoldSysV17(void)
00163 {
00164
00165 DeleteMatrix(&fDAinRelSq);
00166 DeleteMatrix(&fDAinColRelSq);
00167 delete fBgrIn;
00168 delete fBgrErrUncorrInSq;
00169 delete fBgrErrScaleIn;
00170 delete fSysIn;
00171 ClearResults();
00172 delete fDeltaCorrX;
00173 delete fDeltaCorrAx;
00174 DeleteMatrix(&fYData);
00175 DeleteMatrix(&fVyyData);
00176 }
00177
00178
00179
00180
00181 TUnfoldSysV17::TUnfoldSysV17(void)
00182 {
00183
00184 InitTUnfoldSys();
00185 }
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 TUnfoldSysV17::TUnfoldSysV17
00200 (const TH2 *hist_A, EHistMap histmap, ERegMode regmode,EConstraint constraint)
00201 : TUnfold(hist_A,histmap,regmode,constraint)
00202 {
00203
00204
00205
00206
00207 InitTUnfoldSys();
00208
00209
00210 fAoutside = new TMatrixD(GetNx(),2);
00211
00212
00213 fDAinColRelSq = new TMatrixD(GetNx(),1);
00214
00215 Int_t nmax=GetNx()*GetNy();
00216 Int_t *rowDAinRelSq = new Int_t[nmax];
00217 Int_t *colDAinRelSq = new Int_t[nmax];
00218 Double_t *dataDAinRelSq = new Double_t[nmax];
00219
00220 Int_t da_nonzero=0;
00221 for (Int_t ix = 0; ix < GetNx(); ix++) {
00222 Int_t ibinx = fXToHist[ix];
00223 Double_t sum_sq= fSumOverY[ix]*fSumOverY[ix];
00224 for (Int_t ibiny = 0; ibiny <= GetNy()+1; ibiny++) {
00225 Double_t dz;
00226 if (histmap == kHistMapOutputHoriz) {
00227 dz = hist_A->GetBinError(ibinx, ibiny);
00228 } else {
00229 dz = hist_A->GetBinError(ibiny, ibinx);
00230 }
00231 Double_t normerr_sq=dz*dz/sum_sq;
00232
00233
00234 (*fDAinColRelSq)(ix,0) += normerr_sq;
00235
00236 if(ibiny==0) {
00237
00238 if (histmap == kHistMapOutputHoriz) {
00239 (*fAoutside)(ix,0)=hist_A->GetBinContent(ibinx, ibiny);
00240 } else {
00241 (*fAoutside)(ix,0)=hist_A->GetBinContent(ibiny, ibinx);
00242 }
00243 } else if(ibiny==GetNy()+1) {
00244
00245 if (histmap == kHistMapOutputHoriz) {
00246 (*fAoutside)(ix,1)=hist_A->GetBinContent(ibinx, ibiny);
00247 } else {
00248 (*fAoutside)(ix,1)=hist_A->GetBinContent(ibiny, ibinx);
00249 }
00250 } else {
00251
00252 rowDAinRelSq[da_nonzero]=ibiny-1;
00253 colDAinRelSq[da_nonzero] = ix;
00254 dataDAinRelSq[da_nonzero] = normerr_sq;
00255 if(dataDAinRelSq[da_nonzero]>0.0) da_nonzero++;
00256 }
00257 }
00258 }
00259 if(da_nonzero) {
00260 fDAinRelSq = CreateSparseMatrix(GetNy(),GetNx(),da_nonzero,
00261 rowDAinRelSq,colDAinRelSq,dataDAinRelSq);
00262 } else {
00263 DeleteMatrix(&fDAinColRelSq);
00264 }
00265 delete[] rowDAinRelSq;
00266 delete[] colDAinRelSq;
00267 delete[] dataDAinRelSq;
00268 }
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288 void TUnfoldSysV17::AddSysError
00289 (const TH2 *sysError,const char *name,EHistMap histmap,ESysErrMode mode)
00290 {
00291
00292 if(fSysIn->FindObject(name)) {
00293 Error("AddSysError","Source %s given twice, ignoring 2nd call.\n",name);
00294 } else {
00295
00296
00297
00298
00299 TMatrixD aCopy(*fA);
00300
00301 Int_t nmax= GetNx()*GetNy();
00302 Double_t *data=new Double_t[nmax];
00303 Int_t *cols=new Int_t[nmax];
00304 Int_t *rows=new Int_t[nmax];
00305 nmax=0;
00306 for (Int_t ix = 0; ix < GetNx(); ix++) {
00307 Int_t ibinx = fXToHist[ix];
00308 Double_t sum=0.0;
00309 for(Int_t loop=0;loop<2;loop++) {
00310 for (Int_t ibiny = 0; ibiny <= GetNy()+1; ibiny++) {
00311 Double_t z;
00312
00313 if (histmap == kHistMapOutputHoriz) {
00314 z = sysError->GetBinContent(ibinx, ibiny);
00315 } else {
00316 z = sysError->GetBinContent(ibiny, ibinx);
00317 }
00318
00319 if(mode!=kSysErrModeMatrix) {
00320 Double_t z0;
00321 if((ibiny>0)&&(ibiny<=GetNy())) {
00322 z0=aCopy(ibiny-1,ix)*fSumOverY[ix];
00323 } else if(ibiny==0) {
00324 z0=(*fAoutside)(ix,0);
00325 } else {
00326 z0=(*fAoutside)(ix,1);
00327 }
00328 if(mode==kSysErrModeShift) {
00329 z += z0;
00330 } else if(mode==kSysErrModeRelative) {
00331 z = z0*(1.+z);
00332 }
00333 }
00334 if(loop==0) {
00335
00336 sum += z;
00337 } else {
00338 if((ibiny>0)&&(ibiny<=GetNy())) {
00339
00340
00341 rows[nmax]=ibiny-1;
00342 cols[nmax]=ix;
00343 if(sum>0.0) {
00344 data[nmax]=z/sum-aCopy(ibiny-1,ix);
00345 } else {
00346 data[nmax]=0.0;
00347 }
00348 if(data[nmax] != 0.0) nmax++;
00349 }
00350 }
00351 }
00352 }
00353 }
00354 if(nmax==0) {
00355 Error("AddSysError",
00356 "source %s has no influence and has not been added.\n",name);
00357 } else {
00358 TMatrixDSparse *dsys=CreateSparseMatrix(GetNy(),GetNx(),
00359 nmax,rows,cols,data);
00360 fSysIn->Add(new TObjString(name),dsys);
00361 }
00362 delete[] data;
00363 delete[] rows;
00364 delete[] cols;
00365 }
00366 }
00367
00368
00369
00370
00371
00372
00373 void TUnfoldSysV17::DoBackgroundSubtraction(void)
00374 {
00375
00376
00377
00378
00379 if(fYData) {
00380 DeleteMatrix(&fY);
00381 DeleteMatrix(&fVyy);
00382 if(fBgrIn->GetEntries()<=0) {
00383
00384 fY=new TMatrixD(*fYData);
00385 fVyy=new TMatrixDSparse(*fVyyData);
00386 } else {
00387 if(GetVyyInv()) {
00388 Warning("DoBackgroundSubtraction",
00389 "inverse error matrix from user input,"
00390 " not corrected for background");
00391 }
00392
00393 fY=new TMatrixD(*fYData);
00394
00395 const TObject *key;
00396 {
00397 TMapIter bgrPtr(fBgrIn);
00398 for(key=bgrPtr.Next();key;key=bgrPtr.Next()) {
00399 const TMatrixD *bgr=(const TMatrixD *)
00400 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
00401 ((const TPair *)*bgrPtr)->Value()
00402 #else
00403 fBgrIn->GetValue(((const TObjString *)key)->GetString())
00404 #endif
00405 ;
00406 for(Int_t i=0;i<GetNy();i++) {
00407 (*fY)(i,0) -= (*bgr)(i,0);
00408 }
00409 }
00410 }
00411
00412 TMatrixD vyy(*fVyyData);
00413
00414 Int_t ny=fVyyData->GetNrows();
00415 const Int_t *vyydata_rows=fVyyData->GetRowIndexArray();
00416 const Int_t *vyydata_cols=fVyyData->GetColIndexArray();
00417 const Double_t *vyydata_data=fVyyData->GetMatrixArray();
00418 Int_t *usedBin=new Int_t[ny];
00419 for(Int_t i=0;i<ny;i++) {
00420 usedBin[i]=0;
00421 }
00422 for(Int_t i=0;i<ny;i++) {
00423 for(Int_t k=vyydata_rows[i];k<vyydata_rows[i+1];k++) {
00424 if(vyydata_data[k]>0.0) {
00425 usedBin[i]++;
00426 usedBin[vyydata_cols[k]]++;
00427 }
00428 }
00429 }
00430
00431 {
00432 TMapIter bgrErrUncorrSqPtr(fBgrErrUncorrInSq);
00433 for(key=bgrErrUncorrSqPtr.Next();key;
00434 key=bgrErrUncorrSqPtr.Next()) {
00435 const TMatrixD *bgrerruncorrSquared=(TMatrixD const *)
00436 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
00437 ((const TPair *)*bgrErrUncorrSqPtr)->Value()
00438 #else
00439 fBgrErrUncorrInSq->GetValue(((const TObjString *)key)
00440 ->GetString())
00441 #endif
00442 ;
00443 for(Int_t yi=0;yi<ny;yi++) {
00444 if(!usedBin[yi]) continue;
00445 vyy(yi,yi) +=(*bgrerruncorrSquared)(yi,0);
00446 }
00447 }
00448 }
00449
00450 {
00451 TMapIter bgrErrScalePtr(fBgrErrScaleIn);
00452 for(key=bgrErrScalePtr.Next();key;key=bgrErrScalePtr.Next()) {
00453 const TMatrixD *bgrerrscale=(const TMatrixD *)
00454 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
00455 ((const TPair *)*bgrErrScalePtr)->Value()
00456 #else
00457 fBgrErrScaleIn->GetValue(((const TObjString *)key)
00458 ->GetString())
00459 #endif
00460 ;
00461 for(Int_t yi=0;yi<ny;yi++) {
00462 if(!usedBin[yi]) continue;
00463 for(Int_t yj=0;yj<ny;yj++) {
00464 if(!usedBin[yj]) continue;
00465 vyy(yi,yj) +=(*bgrerrscale)(yi,0)* (*bgrerrscale)(yj,0);
00466 }
00467 }
00468 }
00469 }
00470 delete[] usedBin;
00471 usedBin=0;
00472
00473
00474 fVyy=new TMatrixDSparse(vyy);
00475
00476 }
00477 } else {
00478 Fatal("DoBackgroundSubtraction","No input vector defined");
00479 }
00480 }
00481
00482 Int_t TUnfoldSysV17::SetInput(const TH1 *hist_y,Double_t scaleBias,
00483 Double_t oneOverZeroError,const TH2 *hist_vyy,
00484 const TH2 *hist_vyy_inv)
00485 {
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501 Int_t r=TUnfold::SetInput(hist_y,scaleBias,oneOverZeroError,hist_vyy,
00502 hist_vyy_inv);
00503 fYData=fY;
00504 fY=0;
00505 fVyyData=fVyy;
00506 fVyy=0;
00507 DoBackgroundSubtraction();
00508
00509 return r;
00510 }
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531 void TUnfoldSysV17::SubtractBackground
00532 (const TH1 *bgr,const char *name,Double_t scale,Double_t scale_error)
00533 {
00534
00535
00536
00537
00538
00539 if(fBgrIn->FindObject(name)) {
00540 Error("SubtractBackground","Source %s given twice, ignoring 2nd call.\n",
00541 name);
00542 } else {
00543 TMatrixD *bgrScaled=new TMatrixD(GetNy(),1);
00544 TMatrixD *bgrErrUncSq=new TMatrixD(GetNy(),1);
00545 TMatrixD *bgrErrCorr=new TMatrixD(GetNy(),1);
00546 for(Int_t row=0;row<GetNy();row++) {
00547 (*bgrScaled)(row,0) = scale*bgr->GetBinContent(row+1);
00548 (*bgrErrUncSq)(row,0) =
00549 TMath::Power(scale*bgr->GetBinError(row+1),2.);
00550 (*bgrErrCorr)(row,0) = scale_error*bgr->GetBinContent(row+1);
00551 }
00552 fBgrIn->Add(new TObjString(name),bgrScaled);
00553 fBgrErrUncorrInSq->Add(new TObjString(name),bgrErrUncSq);
00554 fBgrErrScaleIn->Add(new TObjString(name),bgrErrCorr);
00555 if(fYData) {
00556 DoBackgroundSubtraction();
00557 } else {
00558 Info("SubtractBackground",
00559 "Background subtraction prior to setting input data");
00560 }
00561 }
00562 }
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581 void TUnfoldSysV17::GetBackground
00582 (TH1 *bgrHist,const char *bgrSource,const Int_t *binMap,
00583 Int_t includeError,Bool_t clearHist) const
00584 {
00585 if(clearHist) {
00586 ClearHistogram(bgrHist);
00587 }
00588
00589 const TObject *key;
00590 {
00591 TMapIter bgrPtr(fBgrIn);
00592 for(key=bgrPtr.Next();key;key=bgrPtr.Next()) {
00593 TString bgrName=((const TObjString *)key)->GetString();
00594 if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
00595 const TMatrixD *bgr=(const TMatrixD *)
00596 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
00597 ((const TPair *)*bgrPtr)->Value()
00598 #else
00599 fBgrIn->GetValue(bgrName)
00600 #endif
00601 ;
00602 for(Int_t i=0;i<GetNy();i++) {
00603 Int_t destBin=binMap[i];
00604 bgrHist->SetBinContent(destBin,bgrHist->GetBinContent(destBin)+
00605 (*bgr)(i,0));
00606 }
00607 }
00608 }
00609
00610 if(includeError &1) {
00611 TMapIter bgrErrUncorrSqPtr(fBgrErrUncorrInSq);
00612 for(key=bgrErrUncorrSqPtr.Next();key;key=bgrErrUncorrSqPtr.Next()) {
00613 TString bgrName=((const TObjString *)key)->GetString();
00614 if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
00615 const TMatrixD *bgrerruncorrSquared=(TMatrixD const *)
00616 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
00617 ((const TPair *)*bgrErrUncorrSqPtr)->Value()
00618 #else
00619 fBgrErrUncorrInSq->GetValue(((const TObjString *)key)
00620 ->GetString())
00621 #endif
00622 ;
00623 for(Int_t i=0;i<GetNy();i++) {
00624 Int_t destBin=binMap[i];
00625 bgrHist->SetBinError
00626 (destBin,TMath::Sqrt
00627 ((*bgrerruncorrSquared)(i,0)+
00628 TMath::Power(bgrHist->GetBinError(destBin),2.)));
00629 }
00630 }
00631 }
00632 if(includeError & 2) {
00633 TMapIter bgrErrScalePtr(fBgrErrScaleIn);
00634 for(key=bgrErrScalePtr.Next();key;key=bgrErrScalePtr.Next()) {
00635 TString bgrName=((const TObjString *)key)->GetString();
00636 if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
00637 const TMatrixD *bgrerrscale=(TMatrixD const *)
00638 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
00639 ((const TPair *)*bgrErrScalePtr)->Value()
00640 #else
00641 fBgrErrScaleIn->GetValue(((const TObjString *)key)->GetString())
00642 #endif
00643 ;
00644 for(Int_t i=0;i<GetNy();i++) {
00645 Int_t destBin=binMap[i];
00646 bgrHist->SetBinError(destBin,hypot((*bgrerrscale)(i,0),
00647 bgrHist->GetBinError(destBin)));
00648 }
00649 }
00650 }
00651 }
00652
00653 void TUnfoldSysV17::InitTUnfoldSys(void)
00654 {
00655
00656
00657
00658 fDAinRelSq = 0;
00659 fDAinColRelSq = 0;
00660 fAoutside = 0;
00661 fBgrIn = new TMap();
00662 fBgrErrUncorrInSq = new TMap();
00663 fBgrErrScaleIn = new TMap();
00664 fSysIn = new TMap();
00665 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
00666 fBgrIn->SetOwnerKeyValue();
00667 fBgrErrUncorrInSq->SetOwnerKeyValue();
00668 fBgrErrScaleIn->SetOwnerKeyValue();
00669 fSysIn->SetOwnerKeyValue();
00670 #else
00671 fBgrIn->SetOwner();
00672 fBgrErrUncorrInSq->SetOwner();
00673 fBgrErrScaleIn->SetOwner();
00674 fSysIn->SetOwner();
00675 #endif
00676
00677 fEmatUncorrX = 0;
00678 fEmatUncorrAx = 0;
00679 fDeltaCorrX = new TMap();
00680 fDeltaCorrAx = new TMap();
00681 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
00682 fDeltaCorrX->SetOwnerKeyValue();
00683 fDeltaCorrAx->SetOwnerKeyValue();
00684 #else
00685 fDeltaCorrX->SetOwner();
00686 fDeltaCorrAx->SetOwner();
00687 #endif
00688 fDeltaSysTau = 0;
00689 fDtau=0.0;
00690 fYData=0;
00691 fVyyData=0;
00692 }
00693
00694 void TUnfoldSysV17::ClearResults(void)
00695 {
00696
00697 TUnfold::ClearResults();
00698 DeleteMatrix(&fEmatUncorrX);
00699 DeleteMatrix(&fEmatUncorrAx);
00700 fDeltaCorrX->Clear();
00701 fDeltaCorrAx->Clear();
00702 DeleteMatrix(&fDeltaSysTau);
00703 }
00704
00705
00706
00707 void TUnfoldSysV17::PrepareSysError(void)
00708 {
00709
00710
00711 if(!fEmatUncorrX) {
00712 fEmatUncorrX=PrepareUncorrEmat(GetDXDAM(0),GetDXDAM(1));
00713 }
00714 TMatrixDSparse *AM0=0,*AM1=0;
00715 if(!fEmatUncorrAx) {
00716 if(!AM0) AM0=MultiplyMSparseMSparse(fA,GetDXDAM(0));
00717 if(!AM1) {
00718 AM1=MultiplyMSparseMSparse(fA,GetDXDAM(1));
00719 Int_t *rows_cols=new Int_t[GetNy()];
00720 Double_t *data=new Double_t[GetNy()];
00721 for(Int_t i=0;i<GetNy();i++) {
00722 rows_cols[i]=i;
00723 data[i]=1.0;
00724 }
00725 TMatrixDSparse *one=CreateSparseMatrix
00726 (GetNy(),GetNy(),GetNy(),rows_cols, rows_cols,data);
00727 delete[] data;
00728 delete[] rows_cols;
00729 AddMSparse(AM1,-1.,one);
00730 DeleteMatrix(&one);
00731 fEmatUncorrAx=PrepareUncorrEmat(AM0,AM1);
00732 }
00733 }
00734 if((!fDeltaSysTau )&&(fDtau>0.0)) {
00735 fDeltaSysTau=new TMatrixDSparse(*GetDXDtauSquared());
00736 Double_t scale=2.*TMath::Sqrt(fTauSquared)*fDtau;
00737 Int_t n=fDeltaSysTau->GetRowIndexArray() [fDeltaSysTau->GetNrows()];
00738 Double_t *data=fDeltaSysTau->GetMatrixArray();
00739 for(Int_t i=0;i<n;i++) {
00740 data[i] *= scale;
00741 }
00742 }
00743
00744 TMapIter sysErrIn(fSysIn);
00745 const TObjString *key;
00746
00747
00748 for(key=(const TObjString *)sysErrIn.Next();key;
00749 key=(const TObjString *)sysErrIn.Next()) {
00750 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
00751 const TMatrixDSparse *dsys=
00752 (const TMatrixDSparse *)((const TPair *)*sysErrIn)->Value();
00753 #else
00754 const TMatrixDSparse *dsys=
00755 (const TMatrixDSparse *)(fSysIn->GetValue(key->GetString()));
00756 #endif
00757 const TPair *named_emat=(const TPair *)
00758 fDeltaCorrX->FindObject(key->GetString());
00759 if(!named_emat) {
00760 TMatrixDSparse *emat=PrepareCorrEmat(GetDXDAM(0),GetDXDAM(1),dsys);
00761 fDeltaCorrX->Add(new TObjString(*key),emat);
00762 }
00763 named_emat=(const TPair *)fDeltaCorrAx->FindObject(key->GetString());
00764 if(!named_emat) {
00765 if(!AM0) AM0=MultiplyMSparseMSparse(fA,GetDXDAM(0));
00766 if(!AM1) {
00767 AM1=MultiplyMSparseMSparse(fA,GetDXDAM(1));
00768 Int_t *rows_cols=new Int_t[GetNy()];
00769 Double_t *data=new Double_t[GetNy()];
00770 for(Int_t i=0;i<GetNy();i++) {
00771 rows_cols[i]=i;
00772 data[i]=1.0;
00773 }
00774 TMatrixDSparse *one=CreateSparseMatrix
00775 (GetNy(),GetNy(),GetNy(),rows_cols, rows_cols,data);
00776 delete[] data;
00777 delete[] rows_cols;
00778 AddMSparse(AM1,-1.,one);
00779 DeleteMatrix(&one);
00780 fEmatUncorrAx=PrepareUncorrEmat(AM0,AM1);
00781 }
00782 TMatrixDSparse *emat=PrepareCorrEmat(AM0,AM1,dsys);
00783 fDeltaCorrAx->Add(new TObjString(*key),emat);
00784 }
00785 }
00786 DeleteMatrix(&AM0);
00787 DeleteMatrix(&AM1);
00788 }
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809 void TUnfoldSysV17::GetEmatrixSysUncorr
00810 (TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
00811 {
00812
00813
00814 PrepareSysError();
00815 ErrorMatrixToHist(ematrix,fEmatUncorrX,binMap,clearEmat);
00816 }
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826 TMatrixDSparse *TUnfoldSysV17::PrepareUncorrEmat
00827 (const TMatrixDSparse *m_0,const TMatrixDSparse *m_1)
00828 {
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921 TMatrixDSparse *r=0;
00922 if(fDAinColRelSq && fDAinRelSq) {
00923
00924 TMatrixDSparse *M1A_Z1=MultiplyMSparseMSparse(m_1,fA);
00925 ScaleColumnsByVector(M1A_Z1,GetDXDAZ(1));
00926 TMatrixDSparse *M1Rsq_Z1=MultiplyMSparseMSparse(m_1,fDAinRelSq);
00927 ScaleColumnsByVector(M1Rsq_Z1,GetDXDAZ(1));
00928
00929 TMatrixDSparse *AtZ0 = MultiplyMSparseTranspMSparse(fA,GetDXDAZ(0));
00930 TMatrixDSparse *RsqZ0=
00931 MultiplyMSparseTranspMSparse(fDAinRelSq,GetDXDAZ(0));
00932
00933
00934 TMatrixDSparse *F=new TMatrixDSparse(*m_0);
00935 ScaleColumnsByVector(F,AtZ0);
00936 AddMSparse(F,-1.0,M1A_Z1);
00937
00938
00939 TMatrixDSparse *G=new TMatrixDSparse(*m_0);
00940 ScaleColumnsByVector(G,RsqZ0);
00941 AddMSparse(G,-1.0,M1Rsq_Z1);
00942 DeleteMatrix(&M1A_Z1);
00943 DeleteMatrix(&M1Rsq_Z1);
00944 DeleteMatrix(&AtZ0);
00945 DeleteMatrix(&RsqZ0);
00946
00947 r=MultiplyMSparseMSparseTranspVector(F,F,fDAinColRelSq);
00948
00949 TMatrixDSparse *r1=MultiplyMSparseMSparseTranspVector(F,G,0);
00950
00951 TMatrixDSparse *r2=MultiplyMSparseMSparseTranspVector(G,F,0);
00952
00953 AddMSparse(r,-1.0,r1);
00954 AddMSparse(r,-1.0,r2);
00955 DeleteMatrix(&r1);
00956 DeleteMatrix(&r2);
00957 DeleteMatrix(&F);
00958 DeleteMatrix(&G);
00959 }
00960
00961
00962
00963
00964 if(fDAinRelSq) {
00965
00966 TMatrixDSparse Z0sq(*GetDXDAZ(0));
00967 const Int_t *Z0sq_rows=Z0sq.GetRowIndexArray();
00968 Double_t *Z0sq_data=Z0sq.GetMatrixArray();
00969 for(int index=0;index<Z0sq_rows[Z0sq.GetNrows()];index++) {
00970 Z0sq_data[index] *= Z0sq_data[index];
00971 }
00972
00973 TMatrixDSparse *Z0sqRsq=MultiplyMSparseTranspMSparse(fDAinRelSq,&Z0sq);
00974
00975 TMatrixDSparse *r3=MultiplyMSparseMSparseTranspVector(m_0,m_0,Z0sqRsq);
00976 DeleteMatrix(&Z0sqRsq);
00977
00978
00979 TMatrixDSparse Z1sq(*GetDXDAZ(1));
00980 const Int_t *Z1sq_rows=Z1sq.GetRowIndexArray();
00981 Double_t *Z1sq_data=Z1sq.GetMatrixArray();
00982 for(int index=0;index<Z1sq_rows[Z1sq.GetNrows()];index++) {
00983 Z1sq_data[index] *= Z1sq_data[index];
00984 }
00985
00986 TMatrixDSparse *Z1sqRsq=MultiplyMSparseMSparse(fDAinRelSq,&Z1sq);
00987
00988 TMatrixDSparse *r4=MultiplyMSparseMSparseTranspVector(m_1,m_1,Z1sqRsq);
00989 DeleteMatrix(&Z1sqRsq);
00990
00991
00992 TMatrixDSparse *H=MultiplyMSparseMSparseTranspVector
00993 (m_0,fDAinRelSq,GetDXDAZ(1));
00994
00995 ScaleColumnsByVector(H,GetDXDAZ(0));
00996
00997 TMatrixDSparse *r5=MultiplyMSparseMSparseTranspVector(m_1,H,0);
00998
00999 TMatrixDSparse *r6=MultiplyMSparseMSparseTranspVector(H,m_1,0);
01000 DeleteMatrix(&H);
01001
01002 if(r) {
01003 AddMSparse(r,1.0,r3);
01004 DeleteMatrix(&r3);
01005 } else {
01006 r=r3;
01007 r3=0;
01008 }
01009 AddMSparse(r,1.0,r4);
01010 AddMSparse(r,-1.0,r5);
01011 AddMSparse(r,-1.0,r6);
01012 DeleteMatrix(&r4);
01013 DeleteMatrix(&r5);
01014 DeleteMatrix(&r6);
01015 }
01016 return r;
01017 }
01018
01019
01020
01021
01022
01023
01024
01025
01026 TMatrixDSparse *TUnfoldSysV17::PrepareCorrEmat
01027 (const TMatrixDSparse *m1,const TMatrixDSparse *m2,const TMatrixDSparse *dsys)
01028 {
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039 TMatrixDSparse *dsysT_VYAx = MultiplyMSparseTranspMSparse(dsys,GetDXDAZ(0));
01040 TMatrixDSparse *delta = MultiplyMSparseMSparse(m1,dsysT_VYAx);
01041 DeleteMatrix(&dsysT_VYAx);
01042 TMatrixDSparse *dsys_X = MultiplyMSparseMSparse(dsys,GetDXDAZ(1));
01043 TMatrixDSparse *delta2 = MultiplyMSparseMSparse(m2,dsys_X);
01044 DeleteMatrix(&dsys_X);
01045 AddMSparse(delta,-1.0,delta2);
01046 DeleteMatrix(&delta2);
01047 return delta;
01048 }
01049
01050
01051
01052
01053
01054
01055
01056 void TUnfoldSysV17::SetTauError(Double_t delta_tau)
01057 {
01058
01059 fDtau=delta_tau;
01060 DeleteMatrix(&fDeltaSysTau);
01061 }
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076 Bool_t TUnfoldSysV17::GetDeltaSysSource(TH1 *hist_delta,const char *name,
01077 const Int_t *binMap)
01078 {
01079 PrepareSysError();
01080 const TPair *named_emat=(const TPair *)fDeltaCorrX->FindObject(name);
01081 const TMatrixDSparse *delta=0;
01082 if(named_emat) {
01083 delta=(TMatrixDSparse *)named_emat->Value();
01084 }
01085 VectorMapToHist(hist_delta,delta,binMap);
01086 return delta !=0;
01087 }
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103 Bool_t TUnfoldSysV17::GetDeltaSysBackgroundScale
01104 (TH1 *hist_delta,const char *source,const Int_t *binMap)
01105 {
01106 PrepareSysError();
01107 const TPair *named_err=(const TPair *)fBgrErrScaleIn->FindObject(source);
01108 TMatrixDSparse *dx=0;
01109 if(named_err) {
01110 const TMatrixD *dy=(TMatrixD *)named_err->Value();
01111 dx=MultiplyMSparseM(GetDXDY(),dy);
01112 }
01113 VectorMapToHist(hist_delta,dx,binMap);
01114 if(dx!=0) {
01115 DeleteMatrix(&dx);
01116 return kTRUE;
01117 }
01118 return kFALSE;
01119 }
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135 Bool_t TUnfoldSysV17::GetDeltaSysTau(TH1 *hist_delta,const Int_t *binMap)
01136 {
01137
01138
01139
01140 PrepareSysError();
01141 VectorMapToHist(hist_delta,fDeltaSysTau,binMap);
01142 return fDeltaSysTau !=0;
01143 }
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162 void TUnfoldSysV17::GetEmatrixSysSource
01163 (TH2 *ematrix,const char *name,const Int_t *binMap,Bool_t clearEmat)
01164 {
01165 PrepareSysError();
01166 const TPair *named_emat=(const TPair *)fDeltaCorrX->FindObject(name);
01167 TMatrixDSparse *emat=0;
01168 if(named_emat) {
01169 const TMatrixDSparse *delta=(TMatrixDSparse *)named_emat->Value();
01170 emat=MultiplyMSparseMSparseTranspVector(delta,delta,0);
01171 }
01172 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
01173 DeleteMatrix(&emat);
01174 }
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193 void TUnfoldSysV17::GetEmatrixSysBackgroundScale
01194 (TH2 *ematrix,const char *name,const Int_t *binMap,Bool_t clearEmat)
01195 {
01196 PrepareSysError();
01197 const TPair *named_err=(const TPair *)fBgrErrScaleIn->FindObject(name);
01198 TMatrixDSparse *emat=0;
01199 if(named_err) {
01200 const TMatrixD *dy=(TMatrixD *)named_err->Value();
01201 TMatrixDSparse *dx=MultiplyMSparseM(GetDXDY(),dy);
01202 emat=MultiplyMSparseMSparseTranspVector(dx,dx,0);
01203 DeleteMatrix(&dx);
01204 }
01205 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
01206 DeleteMatrix(&emat);
01207 }
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225 void TUnfoldSysV17::GetEmatrixSysTau
01226 (TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
01227 {
01228
01229
01230
01231
01232 PrepareSysError();
01233 TMatrixDSparse *emat=0;
01234 if(fDeltaSysTau) {
01235 emat=MultiplyMSparseMSparseTranspVector(fDeltaSysTau,fDeltaSysTau,0);
01236 }
01237 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
01238 DeleteMatrix(&emat);
01239 }
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256 void TUnfoldSysV17::GetEmatrixInput
01257 (TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
01258 {
01259 GetEmatrixFromVyy(fVyyData,ematrix,binMap,clearEmat);
01260 }
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278 void TUnfoldSysV17::GetEmatrixSysBackgroundUncorr
01279 (TH2 *ematrix,const char *source,const Int_t *binMap,Bool_t clearEmat)
01280 {
01281 const TPair *named_err=(const TPair *)fBgrErrUncorrInSq->FindObject(source);
01282 TMatrixDSparse *emat=0;
01283 if(named_err) {
01284 TMatrixD const *dySquared=(TMatrixD const *)named_err->Value();
01285 emat=MultiplyMSparseMSparseTranspVector(GetDXDY(),GetDXDY(),dySquared);
01286 }
01287 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
01288 DeleteMatrix(&emat);
01289 }
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299 void TUnfoldSysV17::GetEmatrixFromVyy
01300 (const TMatrixDSparse *vyy,TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
01301 {
01302
01303
01304
01305
01306
01307 PrepareSysError();
01308 TMatrixDSparse *em=0;
01309 if(vyy) {
01310 TMatrixDSparse *dxdyVyy=MultiplyMSparseMSparse(GetDXDY(),vyy);
01311 em=MultiplyMSparseMSparseTranspVector(dxdyVyy,GetDXDY(),0);
01312 DeleteMatrix(&dxdyVyy);
01313 }
01314 ErrorMatrixToHist(ematrix,em,binMap,clearEmat);
01315 DeleteMatrix(&em);
01316 }
01317
01318
01319
01320
01321
01322
01323
01324
01325 void TUnfoldSysV17::GetEmatrixTotal(TH2 *ematrix,const Int_t *binMap)
01326 {
01327
01328
01329
01330 GetEmatrix(ematrix,binMap);
01331 GetEmatrixSysUncorr(ematrix,binMap,kFALSE);
01332 TMapIter sysErrPtr(fDeltaCorrX);
01333 const TObject *key;
01334
01335 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
01336 GetEmatrixSysSource(ematrix,
01337 ((const TObjString *)key)->GetString(),
01338 binMap,kFALSE);
01339 }
01340 GetEmatrixSysTau(ematrix,binMap,kFALSE);
01341 }
01342
01343
01344
01345 TMatrixDSparse *TUnfoldSysV17::GetSummedErrorMatrixYY(void)
01346 {
01347 PrepareSysError();
01348
01349
01350 TMatrixDSparse *emat_sum=new TMatrixDSparse(*fVyy);
01351
01352
01353 if(fEmatUncorrAx) {
01354 AddMSparse(emat_sum,1.0,fEmatUncorrAx);
01355 }
01356 TMapIter sysErrPtr(fDeltaCorrAx);
01357 const TObject *key;
01358
01359
01360 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
01361 const TMatrixDSparse *delta=(TMatrixDSparse *)
01362 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
01363 ((const TPair *)*sysErrPtr)->Value()
01364 #else
01365 fDeltaCorrAx->GetValue(((const TObjString *)key)
01366 ->GetString())
01367 #endif
01368 ;
01369 TMatrixDSparse *emat=MultiplyMSparseMSparseTranspVector(delta,delta,0);
01370 AddMSparse(emat_sum,1.0,emat);
01371 DeleteMatrix(&emat);
01372 }
01373
01374 if(fDeltaSysTau) {
01375 TMatrixDSparse *Adx_tau=MultiplyMSparseMSparse(fA,fDeltaSysTau);
01376 TMatrixDSparse *emat_tau=
01377 MultiplyMSparseMSparseTranspVector(Adx_tau,Adx_tau,0);
01378 DeleteMatrix(&Adx_tau);
01379 AddMSparse(emat_sum,1.0,emat_tau);
01380 DeleteMatrix(&emat_tau);
01381 }
01382 return emat_sum;
01383 }
01384
01385
01386
01387 TMatrixDSparse *TUnfoldSysV17::GetSummedErrorMatrixXX(void)
01388 {
01389 PrepareSysError();
01390
01391
01392 TMatrixDSparse *emat_sum=new TMatrixDSparse(*GetVxx());
01393
01394
01395 if(fEmatUncorrX) {
01396 AddMSparse(emat_sum,1.0,fEmatUncorrX);
01397 }
01398 TMapIter sysErrPtr(fDeltaCorrX);
01399 const TObject *key;
01400
01401
01402 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
01403 const TMatrixDSparse *delta=(TMatrixDSparse *)
01404 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
01405 ((const TPair *)*sysErrPtr)->Value()
01406 #else
01407 fDeltaCorrX->GetValue(((const TObjString *)key)
01408 ->GetString())
01409 #endif
01410 ;
01411 TMatrixDSparse *emat=MultiplyMSparseMSparseTranspVector(delta,delta,0);
01412 AddMSparse(emat_sum,1.0,emat);
01413 DeleteMatrix(&emat);
01414 }
01415
01416 if(fDeltaSysTau) {
01417 TMatrixDSparse *emat_tau=
01418 MultiplyMSparseMSparseTranspVector(fDeltaSysTau,fDeltaSysTau,0);
01419 AddMSparse(emat_sum,1.0,emat_tau);
01420 DeleteMatrix(&emat_tau);
01421 }
01422 return emat_sum;
01423 }
01424
01425
01426
01427
01428
01429 Double_t TUnfoldSysV17::GetChi2Sys(void)
01430 {
01431
01432 TMatrixDSparse *emat_sum=GetSummedErrorMatrixYY();
01433
01434 Int_t rank=0;
01435 TMatrixDSparse *v=InvertMSparseSymmPos(emat_sum,&rank);
01436 TMatrixD dy(*fY, TMatrixD::kMinus, *GetAx());
01437 TMatrixDSparse *vdy=MultiplyMSparseM(v,&dy);
01438 DeleteMatrix(&v);
01439 Double_t r=0.0;
01440 const Int_t *vdy_rows=vdy->GetRowIndexArray();
01441 const Double_t *vdy_data=vdy->GetMatrixArray();
01442 for(Int_t i=0;i<vdy->GetNrows();i++) {
01443 if(vdy_rows[i+1]>vdy_rows[i]) {
01444 r += vdy_data[vdy_rows[i]] * dy(i,0);
01445 }
01446 }
01447 DeleteMatrix(&vdy);
01448 DeleteMatrix(&emat_sum);
01449 return r;
01450 }
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464 void TUnfoldSysV17::GetRhoItotal(TH1 *rhoi,const Int_t *binMap,TH2 *invEmat)
01465 {
01466
01467
01468
01469
01470
01471
01472 ClearHistogram(rhoi,-1.);
01473 TMatrixDSparse *emat_sum=GetSummedErrorMatrixXX();
01474 GetRhoIFromMatrix(rhoi,emat_sum,binMap,invEmat);
01475
01476 DeleteMatrix(&emat_sum);
01477 }
01478
01479
01480
01481
01482
01483
01484
01485
01486 void TUnfoldSysV17::ScaleColumnsByVector
01487 (TMatrixDSparse *m,const TMatrixTBase<Double_t> *v) const
01488 {
01489
01490
01491
01492
01493 if((m->GetNcols() != v->GetNrows())||(v->GetNcols()!=1)) {
01494 Fatal("ScaleColumnsByVector error",
01495 "matrix cols/vector rows %d!=%d OR vector cols %d !=1\n",
01496 m->GetNcols(),v->GetNrows(),v->GetNcols());
01497 }
01498 const Int_t *rows_m=m->GetRowIndexArray();
01499 const Int_t *cols_m=m->GetColIndexArray();
01500 Double_t *data_m=m->GetMatrixArray();
01501 const TMatrixDSparse *v_sparse=dynamic_cast<const TMatrixDSparse *>(v);
01502 if(v_sparse) {
01503 const Int_t *rows_v=v_sparse->GetRowIndexArray();
01504 const Double_t *data_v=v_sparse->GetMatrixArray();
01505 for(Int_t i=0;i<m->GetNrows();i++) {
01506 for(Int_t index_m=rows_m[i];index_m<rows_m[i+1];index_m++) {
01507 Int_t j=cols_m[index_m];
01508 Int_t index_v=rows_v[j];
01509 if(index_v<rows_v[j+1]) {
01510 data_m[index_m] *= data_v[index_v];
01511 } else {
01512 data_m[index_m] =0.0;
01513 }
01514 }
01515 }
01516 } else {
01517 for(Int_t i=0;i<m->GetNrows();i++) {
01518 for(Int_t index_m=rows_m[i];index_m<rows_m[i+1];index_m++) {
01519 Int_t j=cols_m[index_m];
01520 data_m[index_m] *= (*v)(j,0);
01521 }
01522 }
01523 }
01524 }
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538 void TUnfoldSysV17::VectorMapToHist
01539 (TH1 *hist_delta,const TMatrixDSparse *delta,const Int_t *binMap)
01540 {
01541
01542
01543
01544 Int_t nbin=hist_delta->GetNbinsX();
01545 Double_t *c=new Double_t[nbin+2];
01546 for(Int_t i=0;i<nbin+2;i++) {
01547 c[i]=0.0;
01548 }
01549 if(delta) {
01550 Int_t binMapSize = fHistToX.GetSize();
01551 const Double_t *delta_data=delta->GetMatrixArray();
01552 const Int_t *delta_rows=delta->GetRowIndexArray();
01553 for(Int_t i=0;i<binMapSize;i++) {
01554 Int_t destBinI=binMap ? binMap[i] : i;
01555 Int_t srcBinI=fHistToX[i];
01556 if((destBinI>=0)&&(destBinI<nbin+2)&&(srcBinI>=0)) {
01557 Int_t index=delta_rows[srcBinI];
01558 if(index<delta_rows[srcBinI+1]) {
01559 c[destBinI]+=delta_data[index];
01560 }
01561 }
01562 }
01563 }
01564 for(Int_t i=0;i<nbin+2;i++) {
01565 hist_delta->SetBinContent(i,c[i]);
01566 hist_delta->SetBinError(i,0.0);
01567 }
01568 delete[] c;
01569 }
01570
01571
01572
01573
01574
01575 TSortedList *TUnfoldSysV17::GetSysSources(void) const {
01576
01577 TSortedList *r=new TSortedList();
01578 TMapIter i(fSysIn);
01579 for(const TObject *key=i.Next();key;key=i.Next()) {
01580 r->Add(((TObjString *)key)->Clone());
01581 }
01582 return r;
01583 }
01584
01585
01586
01587
01588
01589 TSortedList *TUnfoldSysV17::GetBgrSources(void) const {
01590
01591 TSortedList *r=new TSortedList();
01592 TMapIter i(fBgrIn);
01593 for(const TObject *key=i.Next();key;key=i.Next()) {
01594 r->Add(((TObjString *)key)->Clone());
01595 }
01596 return r;
01597 }