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 #include <iostream>
00108 #include <TMatrixD.h>
00109 #include <TMatrixDSparse.h>
00110 #include <TMatrixDSym.h>
00111 #include <TMatrixDSymEigen.h>
00112 #include <TMath.h>
00113 #include "TUnfold.h"
00114
00115 #include <map>
00116 #include <vector>
00117
00118
00119
00120
00121
00122 ClassImp(TUnfoldV17)
00123
00124 TUnfoldV17::~TUnfoldV17(void)
00125 {
00126
00127
00128 DeleteMatrix(&fA);
00129 DeleteMatrix(&fL);
00130 DeleteMatrix(&fVyy);
00131 DeleteMatrix(&fY);
00132 DeleteMatrix(&fX0);
00133 DeleteMatrix(&fVyyInv);
00134
00135 ClearResults();
00136 }
00137
00138
00139
00140
00141 void TUnfoldV17::InitTUnfold(void)
00142 {
00143
00144 fXToHist.Set(0);
00145 fHistToX.Set(0);
00146 fSumOverY.Set(0);
00147 fA = 0;
00148 fL = 0;
00149 fVyy = 0;
00150 fY = 0;
00151 fX0 = 0;
00152 fTauSquared = 0.0;
00153 fBiasScale = 0.0;
00154 fConstraint = kEConstraintNone;
00155 fRegMode = kRegModeNone;
00156
00157 fX = 0;
00158 fVyyInv = 0;
00159 fVxx = 0;
00160 fVxxInv = 0;
00161 fAx = 0;
00162 fChi2A = 0.0;
00163 fLXsquared = 0.0;
00164 fRhoMax = 999.0;
00165 fRhoAvg = -1.0;
00166 fNdf = 0;
00167 fDXDAM[0] = 0;
00168 fDXDAZ[0] = 0;
00169 fDXDAM[1] = 0;
00170 fDXDAZ[1] = 0;
00171 fDXDtauSquared = 0;
00172 fDXDY = 0;
00173 fEinv = 0;
00174 fE = 0;
00175 fEpsMatrix=1.E-13;
00176 fIgnoredBins=0;
00177 }
00178
00179
00180
00181
00182
00183
00184
00185 void TUnfoldV17::DeleteMatrix(TMatrixD **m)
00186 {
00187 if(*m) delete *m;
00188 *m=0;
00189 }
00190
00191
00192
00193
00194
00195
00196
00197 void TUnfoldV17::DeleteMatrix(TMatrixDSparse **m)
00198 {
00199 if(*m) delete *m;
00200 *m=0;
00201 }
00202
00203
00204
00205 void TUnfoldV17::ClearResults(void)
00206 {
00207
00208
00209
00210
00211
00212
00213
00214 DeleteMatrix(&fVxx);
00215 DeleteMatrix(&fX);
00216 DeleteMatrix(&fAx);
00217 for(Int_t i=0;i<2;i++) {
00218 DeleteMatrix(fDXDAM+i);
00219 DeleteMatrix(fDXDAZ+i);
00220 }
00221 DeleteMatrix(&fDXDtauSquared);
00222 DeleteMatrix(&fDXDY);
00223 DeleteMatrix(&fEinv);
00224 DeleteMatrix(&fE);
00225 DeleteMatrix(&fVxxInv);
00226 fChi2A = 0.0;
00227 fLXsquared = 0.0;
00228 fRhoMax = 999.0;
00229 fRhoAvg = -1.0;
00230 }
00231
00232
00233
00234
00235 TUnfoldV17::TUnfoldV17(void)
00236 {
00237
00238 InitTUnfold();
00239 }
00240
00241
00242
00243 Double_t TUnfoldV17::DoUnfold(void)
00244 {
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 ClearResults();
00278
00279
00280 if(!fVyyInv) {
00281 GetInputInverseEmatrix(0);
00282 if(fConstraint != kEConstraintNone) {
00283 fNdf--;
00284 }
00285 }
00286
00287
00288
00289
00290
00291 TMatrixDSparse *AtVyyinv=MultiplyMSparseTranspMSparse(fA,fVyyInv);
00292
00293
00294
00295
00296
00297 TMatrixDSparse *rhs=MultiplyMSparseM(AtVyyinv,fY);
00298 TMatrixDSparse *lSquared=MultiplyMSparseTranspMSparse(fL,fL);
00299 if (fBiasScale != 0.0) {
00300 TMatrixDSparse *rhs2=MultiplyMSparseM(lSquared,fX0);
00301 AddMSparse(rhs, fTauSquared * fBiasScale ,rhs2);
00302 DeleteMatrix(&rhs2);
00303 }
00304
00305
00306
00307
00308
00309 fEinv=MultiplyMSparseMSparse(AtVyyinv,fA);
00310 AddMSparse(fEinv,fTauSquared,lSquared);
00311
00312
00313
00314
00315
00316
00317 Int_t rank=0;
00318 fE = InvertMSparseSymmPos(fEinv,&rank);
00319 if(rank != GetNx()) {
00320 Warning("DoUnfold","rank of matrix E %d expect %d",rank,GetNx());
00321 }
00322
00323
00324
00325
00326
00327 TMatrixDSparse *xSparse=MultiplyMSparseMSparse(fE,rhs);
00328 fX = new TMatrixD(*xSparse);
00329 DeleteMatrix(&rhs);
00330 DeleteMatrix(&xSparse);
00331
00332
00333 Double_t lambda_half=0.0;
00334 Double_t one_over_epsEeps=0.0;
00335 TMatrixDSparse *epsilon=0;
00336 TMatrixDSparse *Eepsilon=0;
00337 if(fConstraint != kEConstraintNone) {
00338
00339 const Int_t *A_rows=fA->GetRowIndexArray();
00340 const Int_t *A_cols=fA->GetColIndexArray();
00341 const Double_t *A_data=fA->GetMatrixArray();
00342 TMatrixD epsilonNosparse(fA->GetNcols(),1);
00343 for(Int_t i=0;i<A_rows[fA->GetNrows()];i++) {
00344 epsilonNosparse(A_cols[i],0) += A_data[i];
00345 }
00346 epsilon=new TMatrixDSparse(epsilonNosparse);
00347
00348 Eepsilon=MultiplyMSparseMSparse(fE,epsilon);
00349
00350 TMatrixDSparse *epsilonEepsilon=MultiplyMSparseTranspMSparse(epsilon,
00351 Eepsilon);
00352
00353 if(epsilonEepsilon->GetRowIndexArray()[1]==1) {
00354 one_over_epsEeps=1./epsilonEepsilon->GetMatrixArray()[0];
00355 } else {
00356 Fatal("Unfold","epsilon#Eepsilon has dimension %d != 1",
00357 epsilonEepsilon->GetRowIndexArray()[1]);
00358 }
00359 DeleteMatrix(&epsilonEepsilon);
00360
00361 Double_t y_minus_epsx=0.0;
00362 for(Int_t iy=0;iy<fY->GetNrows();iy++) {
00363 y_minus_epsx += (*fY)(iy,0);
00364 }
00365
00366 for(Int_t ix=0;ix<epsilonNosparse.GetNrows();ix++) {
00367 y_minus_epsx -= epsilonNosparse(ix,0) * (*fX)(ix,0);
00368 }
00369
00370 lambda_half=y_minus_epsx*one_over_epsEeps;
00371
00372 const Int_t *EEpsilon_rows=Eepsilon->GetRowIndexArray();
00373 const Double_t *EEpsilon_data=Eepsilon->GetMatrixArray();
00374 for(Int_t ix=0;ix<Eepsilon->GetNrows();ix++) {
00375 if(EEpsilon_rows[ix]<EEpsilon_rows[ix+1]) {
00376 (*fX)(ix,0) += lambda_half * EEpsilon_data[EEpsilon_rows[ix]];
00377 }
00378 }
00379 }
00380
00381
00382
00383
00384 fDXDY = MultiplyMSparseMSparse(fE,AtVyyinv);
00385
00386
00387 if(fConstraint != kEConstraintNone) {
00388
00389 Int_t *rows=new Int_t[GetNy()];
00390 Int_t *cols=new Int_t[GetNy()];
00391 Double_t *data=new Double_t[GetNy()];
00392 for(Int_t i=0;i<GetNy();i++) {
00393 rows[i]=0;
00394 cols[i]=i;
00395 data[i]=one_over_epsEeps;
00396 }
00397 TMatrixDSparse *temp=CreateSparseMatrix
00398 (1,GetNy(),GetNy(),rows,cols,data);
00399 delete[] data;
00400 delete[] rows;
00401 delete[] cols;
00402
00403 TMatrixDSparse *epsilonB=MultiplyMSparseTranspMSparse(epsilon,fDXDY);
00404
00405 AddMSparse(temp, -one_over_epsEeps, epsilonB);
00406 DeleteMatrix(&epsilonB);
00407
00408 TMatrixDSparse *corr=MultiplyMSparseMSparse(Eepsilon,temp);
00409 DeleteMatrix(&temp);
00410
00411 AddMSparse(fDXDY,1.0,corr);
00412 DeleteMatrix(&corr);
00413 }
00414
00415 DeleteMatrix(&AtVyyinv);
00416
00417
00418
00419
00420 TMatrixDSparse *fDXDYVyy = MultiplyMSparseMSparse(fDXDY,fVyy);
00421 fVxx = MultiplyMSparseMSparseTranspVector(fDXDYVyy,fDXDY,0);
00422
00423 DeleteMatrix(&fDXDYVyy);
00424
00425
00426
00427
00428
00429 fAx = MultiplyMSparseM(fA,fX);
00430
00431
00432
00433
00434
00435 TMatrixD dy(*fY, TMatrixD::kMinus, *fAx);
00436 TMatrixDSparse *VyyinvDy=MultiplyMSparseM(fVyyInv,&dy);
00437
00438 const Int_t *VyyinvDy_rows=VyyinvDy->GetRowIndexArray();
00439 const Double_t *VyyinvDy_data=VyyinvDy->GetMatrixArray();
00440 fChi2A=0.0;
00441 for(Int_t iy=0;iy<VyyinvDy->GetNrows();iy++) {
00442 if(VyyinvDy_rows[iy]<VyyinvDy_rows[iy+1]) {
00443 fChi2A += VyyinvDy_data[VyyinvDy_rows[iy]]*dy(iy,0);
00444 }
00445 }
00446 TMatrixD dx( fBiasScale * (*fX0), TMatrixD::kMinus,(*fX));
00447 TMatrixDSparse *LsquaredDx=MultiplyMSparseM(lSquared,&dx);
00448 const Int_t *LsquaredDx_rows=LsquaredDx->GetRowIndexArray();
00449 const Double_t *LsquaredDx_data=LsquaredDx->GetMatrixArray();
00450 fLXsquared = 0.0;
00451 for(Int_t ix=0;ix<LsquaredDx->GetNrows();ix++) {
00452 if(LsquaredDx_rows[ix]<LsquaredDx_rows[ix+1]) {
00453 fLXsquared += LsquaredDx_data[LsquaredDx_rows[ix]]*dx(ix,0);
00454 }
00455 }
00456
00457
00458
00459 fDXDtauSquared=MultiplyMSparseMSparse(fE,LsquaredDx);
00460
00461 if(fConstraint != kEConstraintNone) {
00462 TMatrixDSparse *temp=MultiplyMSparseTranspMSparse(epsilon,fDXDtauSquared);
00463 Double_t f=0.0;
00464 if(temp->GetRowIndexArray()[1]==1) {
00465 f=temp->GetMatrixArray()[0]*one_over_epsEeps;
00466 } else if(temp->GetRowIndexArray()[1]>1) {
00467 Fatal("Unfold",
00468 "epsilon#fDXDtauSquared has dimension %d != 1",
00469 temp->GetRowIndexArray()[1]);
00470 }
00471 if(f!=0.0) {
00472 AddMSparse(fDXDtauSquared, -f,Eepsilon);
00473 }
00474 DeleteMatrix(&temp);
00475 }
00476 DeleteMatrix(&epsilon);
00477
00478 DeleteMatrix(&LsquaredDx);
00479 DeleteMatrix(&lSquared);
00480
00481
00482 fDXDAM[0]=new TMatrixDSparse(*fE);
00483 fDXDAM[1]=new TMatrixDSparse(*fDXDY);
00484 fDXDAZ[0]=VyyinvDy;
00485 VyyinvDy=0;
00486 fDXDAZ[1]=new TMatrixDSparse(*fX);
00487
00488 if(fConstraint != kEConstraintNone) {
00489
00490 TMatrixDSparse *temp1=MultiplyMSparseMSparseTranspVector
00491 (Eepsilon,Eepsilon,0);
00492 AddMSparse(fDXDAM[0], -one_over_epsEeps,temp1);
00493 DeleteMatrix(&temp1);
00494
00495 Int_t *rows=new Int_t[GetNy()];
00496 Int_t *cols=new Int_t[GetNy()];
00497 Double_t *data=new Double_t[GetNy()];
00498 for(Int_t i=0;i<GetNy();i++) {
00499 rows[i]=i;
00500 cols[i]=0;
00501 data[i]=lambda_half;
00502 }
00503 TMatrixDSparse *temp2=CreateSparseMatrix
00504 (GetNy(),1,GetNy(),rows,cols,data);
00505 delete[] data;
00506 delete[] rows;
00507 delete[] cols;
00508 AddMSparse(fDXDAZ[0],1.0,temp2);
00509 DeleteMatrix(&temp2);
00510 }
00511
00512 DeleteMatrix(&Eepsilon);
00513
00514 rank=0;
00515 fVxxInv = InvertMSparseSymmPos(fVxx,&rank);
00516 if(rank != GetNx()) {
00517 Warning("DoUnfold","rank of output covariance is %d expect %d",
00518 rank,GetNx());
00519 }
00520
00521 TVectorD VxxInvDiag(fVxxInv->GetNrows());
00522 const Int_t *VxxInv_rows=fVxxInv->GetRowIndexArray();
00523 const Int_t *VxxInv_cols=fVxxInv->GetColIndexArray();
00524 const Double_t *VxxInv_data=fVxxInv->GetMatrixArray();
00525 for (int ix = 0; ix < fVxxInv->GetNrows(); ix++) {
00526 for(int ik=VxxInv_rows[ix];ik<VxxInv_rows[ix+1];ik++) {
00527 if(ix==VxxInv_cols[ik]) {
00528 VxxInvDiag(ix)=VxxInv_data[ik];
00529 }
00530 }
00531 }
00532
00533
00534 const Int_t *Vxx_rows=fVxx->GetRowIndexArray();
00535 const Int_t *Vxx_cols=fVxx->GetColIndexArray();
00536 const Double_t *Vxx_data=fVxx->GetMatrixArray();
00537
00538 Double_t rho_squared_max = 0.0;
00539 Double_t rho_sum = 0.0;
00540 Int_t n_rho=0;
00541 for (int ix = 0; ix < fVxx->GetNrows(); ix++) {
00542 for(int ik=Vxx_rows[ix];ik<Vxx_rows[ix+1];ik++) {
00543 if(ix==Vxx_cols[ik]) {
00544 Double_t rho_squared =
00545 1. - 1. / (VxxInvDiag(ix) * Vxx_data[ik]);
00546 if (rho_squared > rho_squared_max)
00547 rho_squared_max = rho_squared;
00548 if(rho_squared>0.0) {
00549 rho_sum += TMath::Sqrt(rho_squared);
00550 n_rho++;
00551 }
00552 break;
00553 }
00554 }
00555 }
00556 fRhoMax = TMath::Sqrt(rho_squared_max);
00557 fRhoAvg = (n_rho>0) ? (rho_sum/n_rho) : -1.0;
00558
00559 return fRhoMax;
00560 }
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575 TMatrixDSparse *TUnfoldV17::CreateSparseMatrix
00576 (Int_t nrow,Int_t ncol,Int_t nel,Int_t *row,Int_t *col,Double_t *data) const
00577 {
00578
00579
00580
00581
00582 TMatrixDSparse *A=new TMatrixDSparse(nrow,ncol);
00583 if(nel>0) {
00584 A->SetMatrixArray(nel,row,col,data);
00585 }
00586 return A;
00587 }
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600 TMatrixDSparse *TUnfoldV17::MultiplyMSparseMSparse(const TMatrixDSparse *a,
00601 const TMatrixDSparse *b) const
00602 {
00603 if(a->GetNcols()!=b->GetNrows()) {
00604 Fatal("MultiplyMSparseMSparse",
00605 "inconsistent matrix col/ matrix row %d !=%d",
00606 a->GetNcols(),b->GetNrows());
00607 }
00608
00609 TMatrixDSparse *r=new TMatrixDSparse(a->GetNrows(),b->GetNcols());
00610 const Int_t *a_rows=a->GetRowIndexArray();
00611 const Int_t *a_cols=a->GetColIndexArray();
00612 const Double_t *a_data=a->GetMatrixArray();
00613 const Int_t *b_rows=b->GetRowIndexArray();
00614 const Int_t *b_cols=b->GetColIndexArray();
00615 const Double_t *b_data=b->GetMatrixArray();
00616
00617 int nMax=0;
00618 for (Int_t irow = 0; irow < a->GetNrows(); irow++) {
00619 if(a_rows[irow+1]>a_rows[irow]) nMax += b->GetNcols();
00620 }
00621 if((nMax>0)&&(a_cols)&&(b_cols)) {
00622 Int_t *r_rows=new Int_t[nMax];
00623 Int_t *r_cols=new Int_t[nMax];
00624 Double_t *r_data=new Double_t[nMax];
00625 Double_t *row_data=new Double_t[b->GetNcols()];
00626 Int_t n=0;
00627 for (Int_t irow = 0; irow < a->GetNrows(); irow++) {
00628 if(a_rows[irow+1]<=a_rows[irow]) continue;
00629
00630 for(Int_t icol=0;icol<b->GetNcols();icol++) {
00631 row_data[icol]=0.0;
00632 }
00633
00634 for(Int_t ia=a_rows[irow];ia<a_rows[irow+1];ia++) {
00635 Int_t k=a_cols[ia];
00636
00637 for(Int_t ib=b_rows[k];ib<b_rows[k+1];ib++) {
00638 row_data[b_cols[ib]] += a_data[ia]*b_data[ib];
00639 }
00640 }
00641
00642 for(Int_t icol=0;icol<b->GetNcols();icol++) {
00643 if(row_data[icol] != 0.0) {
00644 r_rows[n]=irow;
00645 r_cols[n]=icol;
00646 r_data[n]=row_data[icol];
00647 n++;
00648 }
00649 }
00650 }
00651 if(n>0) {
00652 r->SetMatrixArray(n,r_rows,r_cols,r_data);
00653 }
00654 delete[] r_rows;
00655 delete[] r_cols;
00656 delete[] r_data;
00657 delete[] row_data;
00658 }
00659
00660 return r;
00661 }
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674 TMatrixDSparse *TUnfoldV17::MultiplyMSparseTranspMSparse
00675 (const TMatrixDSparse *a,const TMatrixDSparse *b) const
00676 {
00677 if(a->GetNrows() != b->GetNrows()) {
00678 Fatal("MultiplyMSparseTranspMSparse",
00679 "inconsistent matrix row numbers %d!=%d",
00680 a->GetNrows(),b->GetNrows());
00681 }
00682
00683 TMatrixDSparse *r=new TMatrixDSparse(a->GetNcols(),b->GetNcols());
00684 const Int_t *a_rows=a->GetRowIndexArray();
00685 const Int_t *a_cols=a->GetColIndexArray();
00686 const Double_t *a_data=a->GetMatrixArray();
00687 const Int_t *b_rows=b->GetRowIndexArray();
00688 const Int_t *b_cols=b->GetColIndexArray();
00689 const Double_t *b_data=b->GetMatrixArray();
00690
00691
00692
00693 typedef std::map<Int_t,Double_t> MMatrixRow_t;
00694 typedef std::map<Int_t, MMatrixRow_t > MMatrix_t;
00695 MMatrix_t matrix;
00696
00697 for(Int_t iRowAB=0;iRowAB<a->GetNrows();iRowAB++) {
00698 for(Int_t ia=a_rows[iRowAB];ia<a_rows[iRowAB+1];ia++) {
00699 for(Int_t ib=b_rows[iRowAB];ib<b_rows[iRowAB+1];ib++) {
00700
00701 MMatrixRow_t &row=matrix[a_cols[ia]];
00702 MMatrixRow_t::iterator icol=row.find(b_cols[ib]);
00703 if(icol!=row.end()) {
00704
00705 (*icol).second += a_data[ia]*b_data[ib];
00706 } else {
00707
00708 row[b_cols[ib]] = a_data[ia]*b_data[ib];
00709 }
00710 }
00711 }
00712 }
00713
00714 Int_t n=0;
00715 for(MMatrix_t::const_iterator irow=matrix.begin();
00716 irow!=matrix.end();irow++) {
00717 n += (*irow).second.size();
00718 }
00719 if(n>0) {
00720
00721 Int_t *r_rows=new Int_t[n];
00722 Int_t *r_cols=new Int_t[n];
00723 Double_t *r_data=new Double_t[n];
00724 n=0;
00725 for(MMatrix_t::const_iterator irow=matrix.begin();
00726 irow!=matrix.end();irow++) {
00727 for(MMatrixRow_t::const_iterator icol=(*irow).second.begin();
00728 icol!=(*irow).second.end();icol++) {
00729 r_rows[n]=(*irow).first;
00730 r_cols[n]=(*icol).first;
00731 r_data[n]=(*icol).second;
00732 n++;
00733 }
00734 }
00735
00736 if(n>0) {
00737 r->SetMatrixArray(n,r_rows,r_cols,r_data);
00738 }
00739 delete[] r_rows;
00740 delete[] r_cols;
00741 delete[] r_data;
00742 }
00743
00744 return r;
00745 }
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757 TMatrixDSparse *TUnfoldV17::MultiplyMSparseM(const TMatrixDSparse *a,
00758 const TMatrixD *b) const
00759 {
00760 if(a->GetNcols()!=b->GetNrows()) {
00761 Fatal("MultiplyMSparseM","inconsistent matrix col /matrix row %d!=%d",
00762 a->GetNcols(),b->GetNrows());
00763 }
00764
00765 TMatrixDSparse *r=new TMatrixDSparse(a->GetNrows(),b->GetNcols());
00766 const Int_t *a_rows=a->GetRowIndexArray();
00767 const Int_t *a_cols=a->GetColIndexArray();
00768 const Double_t *a_data=a->GetMatrixArray();
00769
00770 int nMax=0;
00771 for (Int_t irow = 0; irow < a->GetNrows(); irow++) {
00772 if(a_rows[irow+1]-a_rows[irow]>0) nMax += b->GetNcols();
00773 }
00774 if(nMax>0) {
00775 Int_t *r_rows=new Int_t[nMax];
00776 Int_t *r_cols=new Int_t[nMax];
00777 Double_t *r_data=new Double_t[nMax];
00778
00779 Int_t n=0;
00780
00781 for (Int_t irow = 0; irow < a->GetNrows(); irow++) {
00782 if(a_rows[irow+1]-a_rows[irow]<=0) continue;
00783 for(Int_t icol=0;icol<b->GetNcols();icol++) {
00784 r_rows[n]=irow;
00785 r_cols[n]=icol;
00786 r_data[n]=0.0;
00787 for(Int_t i=a_rows[irow];i<a_rows[irow+1];i++) {
00788 Int_t j=a_cols[i];
00789 r_data[n] += a_data[i]*(*b)(j,icol);
00790 }
00791 if(r_data[n]!=0.0) n++;
00792 }
00793 }
00794 if(n>0) {
00795 r->SetMatrixArray(n,r_rows,r_cols,r_data);
00796 }
00797 delete[] r_rows;
00798 delete[] r_cols;
00799 delete[] r_data;
00800 }
00801 return r;
00802 }
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816 TMatrixDSparse *TUnfoldV17::MultiplyMSparseMSparseTranspVector
00817 (const TMatrixDSparse *m1,const TMatrixDSparse *m2,
00818 const TMatrixTBase<Double_t> *v) const
00819 {
00820 if((m1->GetNcols() != m2->GetNcols())||
00821 (v && ((m1->GetNcols()!=v->GetNrows())||(v->GetNcols()!=1)))) {
00822 if(v) {
00823 Fatal("MultiplyMSparseMSparseTranspVector",
00824 "matrix cols/vector rows %d!=%d!=%d or vector rows %d!=1\n",
00825 m1->GetNcols(),m2->GetNcols(),v->GetNrows(),v->GetNcols());
00826 } else {
00827 Fatal("MultiplyMSparseMSparseTranspVector",
00828 "matrix cols %d!=%d\n",m1->GetNcols(),m2->GetNcols());
00829 }
00830 }
00831 const Int_t *rows_m1=m1->GetRowIndexArray();
00832 const Int_t *cols_m1=m1->GetColIndexArray();
00833 const Double_t *data_m1=m1->GetMatrixArray();
00834 Int_t num_m1=0;
00835 for(Int_t i=0;i<m1->GetNrows();i++) {
00836 if(rows_m1[i]<rows_m1[i+1]) num_m1++;
00837 }
00838 const Int_t *rows_m2=m2->GetRowIndexArray();
00839 const Int_t *cols_m2=m2->GetColIndexArray();
00840 const Double_t *data_m2=m2->GetMatrixArray();
00841 Int_t num_m2=0;
00842 for(Int_t j=0;j<m2->GetNrows();j++) {
00843 if(rows_m2[j]<rows_m2[j+1]) num_m2++;
00844 }
00845 const TMatrixDSparse *v_sparse=dynamic_cast<const TMatrixDSparse *>(v);
00846 const Int_t *v_rows=0;
00847 const Double_t *v_data=0;
00848 if(v_sparse) {
00849 v_rows=v_sparse->GetRowIndexArray();
00850 v_data=v_sparse->GetMatrixArray();
00851 }
00852 Int_t num_r=num_m1*num_m2+1;
00853 Int_t *row_r=new Int_t[num_r];
00854 Int_t *col_r=new Int_t[num_r];
00855 Double_t *data_r=new Double_t[num_r];
00856 num_r=0;
00857 for(Int_t i=0;i<m1->GetNrows();i++) {
00858 for(Int_t j=0;j<m2->GetNrows();j++) {
00859 data_r[num_r]=0.0;
00860 Int_t index_m1=rows_m1[i];
00861 Int_t index_m2=rows_m2[j];
00862 while((index_m1<rows_m1[i+1])&&(index_m2<rows_m2[j+1])) {
00863 Int_t k1=cols_m1[index_m1];
00864 Int_t k2=cols_m2[index_m2];
00865 if(k1<k2) {
00866 index_m1++;
00867 } else if(k1>k2) {
00868 index_m2++;
00869 } else {
00870 if(v_sparse) {
00871 Int_t v_index=v_rows[k1];
00872 if(v_index<v_rows[k1+1]) {
00873 data_r[num_r] += data_m1[index_m1] * data_m2[index_m2]
00874 * v_data[v_index];
00875 }
00876 } else if(v) {
00877 data_r[num_r] += data_m1[index_m1] * data_m2[index_m2]
00878 * (*v)(k1,0);
00879 } else {
00880 data_r[num_r] += data_m1[index_m1] * data_m2[index_m2];
00881 }
00882 index_m1++;
00883 index_m2++;
00884 }
00885 }
00886 if(data_r[num_r] !=0.0) {
00887 row_r[num_r]=i;
00888 col_r[num_r]=j;
00889 num_r++;
00890 }
00891 }
00892 }
00893 TMatrixDSparse *r=CreateSparseMatrix(m1->GetNrows(),m2->GetNrows(),
00894 num_r,row_r,col_r,data_r);
00895 delete[] row_r;
00896 delete[] col_r;
00897 delete[] data_r;
00898 return r;
00899 }
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912 void TUnfoldV17::AddMSparse(TMatrixDSparse *dest,Double_t f,
00913 const TMatrixDSparse *src) const
00914 {
00915 const Int_t *dest_rows=dest->GetRowIndexArray();
00916 const Int_t *dest_cols=dest->GetColIndexArray();
00917 const Double_t *dest_data=dest->GetMatrixArray();
00918 const Int_t *src_rows=src->GetRowIndexArray();
00919 const Int_t *src_cols=src->GetColIndexArray();
00920 const Double_t *src_data=src->GetMatrixArray();
00921
00922 if((dest->GetNrows()!=src->GetNrows())||
00923 (dest->GetNcols()!=src->GetNcols())) {
00924 Fatal("AddMSparse","inconsistent matrix rows %d!=%d OR cols %d!=%d",
00925 src->GetNrows(),dest->GetNrows(),src->GetNcols(),dest->GetNcols());
00926 }
00927 Int_t nmax=dest->GetNrows()*dest->GetNcols();
00928 Double_t *result_data=new Double_t[nmax];
00929 Int_t *result_rows=new Int_t[nmax];
00930 Int_t *result_cols=new Int_t[nmax];
00931 Int_t n=0;
00932 for(Int_t row=0;row<dest->GetNrows();row++) {
00933 Int_t i_dest=dest_rows[row];
00934 Int_t i_src=src_rows[row];
00935 while((i_dest<dest_rows[row+1])||(i_src<src_rows[row+1])) {
00936 Int_t col_dest=(i_dest<dest_rows[row+1]) ?
00937 dest_cols[i_dest] : dest->GetNcols();
00938 Int_t col_src =(i_src <src_rows[row+1] ) ?
00939 src_cols [i_src] : src->GetNcols();
00940 result_rows[n]=row;
00941 if(col_dest<col_src) {
00942 result_cols[n]=col_dest;
00943 result_data[n]=dest_data[i_dest++];
00944 } else if(col_dest>col_src) {
00945 result_cols[n]=col_src;
00946 result_data[n]=f*src_data[i_src++];
00947 } else {
00948 result_cols[n]=col_dest;
00949 result_data[n]=dest_data[i_dest++]+f*src_data[i_src++];
00950 }
00951 if(result_data[n] !=0.0) {
00952 if(!TMath::Finite(result_data[n])) {
00953 Fatal("AddMSparse","Nan detected %d %d %d",
00954 row,i_dest,i_src);
00955 }
00956 n++;
00957 }
00958 }
00959 }
00960 if(n<=0) {
00961 n=1;
00962 result_rows[0]=0;
00963 result_cols[0]=0;
00964 result_data[0]=0.0;
00965 }
00966 dest->SetMatrixArray(n,result_rows,result_cols,result_data);
00967 delete[] result_data;
00968 delete[] result_rows;
00969 delete[] result_cols;
00970 }
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989 TMatrixDSparse *TUnfoldV17::InvertMSparseSymmPos
00990 (const TMatrixDSparse *A,Int_t *rankPtr) const
00991 {
00992
00993 if(A->GetNcols()!=A->GetNrows()) {
00994 Fatal("InvertMSparseSymmPos","inconsistent matrix row/col %d!=%d",
00995 A->GetNcols(),A->GetNrows());
00996 }
00997
00998 Bool_t *isZero=new Bool_t[A->GetNrows()];
00999 const Int_t *a_rows=A->GetRowIndexArray();
01000 const Int_t *a_cols=A->GetColIndexArray();
01001 const Double_t *a_data=A->GetMatrixArray();
01002
01003
01004
01005 Int_t iDiagonal=0;
01006 Int_t iBlock=A->GetNrows();
01007 Bool_t isDiagonal=kTRUE;
01008 TVectorD aII(A->GetNrows());
01009 Int_t nError=0;
01010 for(Int_t iA=0;iA<A->GetNrows();iA++) {
01011 for(Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
01012 Int_t jA=a_cols[indexA];
01013 if(iA==jA) {
01014 if(!(a_data[indexA]>=0.0)) nError++;
01015 aII(iA)=a_data[indexA];
01016 }
01017 }
01018 }
01019 if(nError>0) {
01020 Fatal("InvertMSparseSymmPos",
01021 "Matrix has %d negative elements on the diagonal", nError);
01022 return 0;
01023 }
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034 Int_t *swap=new Int_t[A->GetNrows()];
01035 for(Int_t j=0;j<A->GetNrows();j++) swap[j]=j;
01036 for(Int_t iStart=0;iStart<iBlock;iStart++) {
01037 Int_t nZero=0;
01038 for(Int_t i=iStart;i<iBlock;i++) {
01039 Int_t iA=swap[i];
01040 Int_t n=A->GetNrows()-(a_rows[iA+1]-a_rows[iA]);
01041 if(n>nZero) {
01042 Int_t tmp=swap[iStart];
01043 swap[iStart]=swap[i];
01044 swap[i]=tmp;
01045 nZero=n;
01046 }
01047 }
01048 for(Int_t i=0;i<A->GetNrows();i++) isZero[i]=kTRUE;
01049 Int_t iA=swap[iStart];
01050 for(Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
01051 Int_t jA=a_cols[indexA];
01052 isZero[jA]=kFALSE;
01053 if(iA!=jA) {
01054 isDiagonal=kFALSE;
01055 }
01056 }
01057 if(isDiagonal) {
01058 iDiagonal=iStart+1;
01059 } else {
01060 for(Int_t i=iStart+1;i<iBlock;) {
01061 if(isZero[swap[i]]) {
01062 i++;
01063 } else {
01064 iBlock--;
01065 Int_t tmp=swap[iBlock];
01066 swap[iBlock]=swap[i];
01067 swap[i]=tmp;
01068 }
01069 }
01070 }
01071 }
01072
01073
01074
01075
01076
01077
01078 #ifdef CONDITION_BLOCK_PART
01079 Int_t nn=A->GetNrows()-iBlock;
01080 for(int inc=(nn+1)/2;inc;inc /= 2) {
01081 for(int i=inc;i<nn;i++) {
01082 int itmp=swap[i+iBlock];
01083 int j;
01084 for(j=i;(j>=inc)&&(aII(swap[j-inc+iBlock])<aII(itmp));j -=inc) {
01085 swap[j+iBlock]=swap[j-inc+iBlock];
01086 }
01087 swap[j+iBlock]=itmp;
01088 }
01089 }
01090 #endif
01091
01092 Int_t *swapBack=new Int_t[A->GetNrows()];
01093 for(Int_t i=0;i<A->GetNrows();i++) {
01094 swapBack[swap[i]]=i;
01095 }
01096 #ifdef DEBUG_DETAIL
01097 for(Int_t i=0;i<A->GetNrows();i++) {
01098 std::cout<<" "<<i<<" "<<swap[i]<<" "<<swapBack[i]<<"\n";
01099 }
01100 std::cout<<"after sorting\n";
01101 for(Int_t i=0;i<A->GetNrows();i++) {
01102 if(i==iDiagonal) std::cout<<"iDiagonal="<<i<<"\n";
01103 if(i==iBlock) std::cout<<"iBlock="<<i<<"\n";
01104 std::cout<<" "<<swap[i]<<" "<<aII(swap[i])<<"\n";
01105 }
01106 {
01107
01108
01109
01110
01111 Int_t nDiag=0;
01112 Int_t nError1=0;
01113 Int_t nError2=0;
01114 Int_t nBlock=0;
01115 Int_t nNonzero=0;
01116 for(Int_t i=0;i<A->GetNrows();i++) {
01117 Int_t row=swap[i];
01118 for(Int_t indexA=a_rows[row];indexA<a_rows[row+1];indexA++) {
01119 Int_t j=swapBack[a_cols[indexA]];
01120 if(i==j) nDiag++;
01121 else if((i<iDiagonal)||(j<iDiagonal)) nError1++;
01122 else if((i<iBlock)&&(j<iBlock)) nError2++;
01123 else if((i<iBlock)||(j<iBlock)) nBlock++;
01124 else nNonzero++;
01125 }
01126 }
01127 if(nError1+nError2>0) {
01128 Fatal("InvertMSparseSymmPos","sparse matrix analysis failed %d %d %d %d %d",
01129 iDiagonal,iBlock,A->GetNrows(),nError1,nError2);
01130 }
01131 }
01132 #endif
01133 #ifdef DEBUG
01134 Info("InvertMSparseSymmPos","iDiagonal=%d iBlock=%d nRow=%d",
01135 iDiagonal,iBlock,A->GetNrows());
01136 #endif
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
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 Double_t *rEl_data=new Double_t[A->GetNrows()*A->GetNrows()];
01187 Int_t *rEl_col=new Int_t[A->GetNrows()*A->GetNrows()];
01188 Int_t *rEl_row=new Int_t[A->GetNrows()*A->GetNrows()];
01189 Int_t rNumEl=0;
01190
01191
01192
01193 Int_t rankD1=0;
01194 for(Int_t i=0;i<iDiagonal;++i) {
01195 Int_t iA=swap[i];
01196 if(aII(iA)>0.0) {
01197 rEl_col[rNumEl]=iA;
01198 rEl_row[rNumEl]=iA;
01199 rEl_data[rNumEl]=1./aII(iA);
01200 ++rankD1;
01201 ++rNumEl;
01202 }
01203 }
01204 if((!rankPtr)&&(rankD1!=iDiagonal)) {
01205 Fatal("InvertMSparseSymmPos",
01206 "diagonal part 1 has rank %d != %d, matrix can not be inverted",
01207 rankD1,iDiagonal);
01208 rNumEl=-1;
01209 }
01210
01211
01212
01213
01214 Int_t nD2=iBlock-iDiagonal;
01215 TMatrixDSparse *D2inv=0;
01216 if((rNumEl>=0)&&(nD2>0)) {
01217 Double_t *D2inv_data=new Double_t[nD2];
01218 Int_t *D2inv_col=new Int_t[nD2];
01219 Int_t *D2inv_row=new Int_t[nD2];
01220 Int_t D2invNumEl=0;
01221 for(Int_t i=0;i<nD2;++i) {
01222 Int_t iA=swap[i+iDiagonal];
01223 if(aII(iA)>0.0) {
01224 D2inv_col[D2invNumEl]=i;
01225 D2inv_row[D2invNumEl]=i;
01226 D2inv_data[D2invNumEl]=1./aII(iA);
01227 ++D2invNumEl;
01228 }
01229 }
01230 if(D2invNumEl==nD2) {
01231 D2inv=CreateSparseMatrix(nD2,nD2,D2invNumEl,D2inv_row,D2inv_col,
01232 D2inv_data);
01233 } else if(!rankPtr) {
01234 Fatal("InvertMSparseSymmPos",
01235 "diagonal part 2 has rank %d != %d, matrix can not be inverted",
01236 D2invNumEl,nD2);
01237 rNumEl=-2;
01238 }
01239 delete [] D2inv_data;
01240 delete [] D2inv_col;
01241 delete [] D2inv_row;
01242 }
01243
01244
01245
01246
01247 Int_t nF=A->GetNrows()-iBlock;
01248 TMatrixDSparse *Finv=0;
01249 TMatrixDSparse *B=0;
01250 TMatrixDSparse *minusBD2inv=0;
01251 if((rNumEl>=0)&&(nF>0)&&((nD2==0)||D2inv)) {
01252
01253 Int_t nFmax=nF*nF;
01254 Double_t epsilonF2=fEpsMatrix;
01255 Double_t *F_data=new Double_t[nFmax];
01256 Int_t *F_col=new Int_t[nFmax];
01257 Int_t *F_row=new Int_t[nFmax];
01258 Int_t FnumEl=0;
01259
01260 Int_t nBmax=nF*(nD2+1);
01261 Double_t *B_data=new Double_t[nBmax];
01262 Int_t *B_col=new Int_t[nBmax];
01263 Int_t *B_row=new Int_t[nBmax];
01264 Int_t BnumEl=0;
01265
01266 for(Int_t i=0;i<nF;i++) {
01267 Int_t iA=swap[i+iBlock];
01268 for(Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
01269 Int_t jA=a_cols[indexA];
01270 Int_t j0=swapBack[jA];
01271 if(j0>=iBlock) {
01272 Int_t j=j0-iBlock;
01273 F_row[FnumEl]=i;
01274 F_col[FnumEl]=j;
01275 F_data[FnumEl]=a_data[indexA];
01276 FnumEl++;
01277 } else if(j0>=iDiagonal) {
01278 Int_t j=j0-iDiagonal;
01279 B_row[BnumEl]=i;
01280 B_col[BnumEl]=j;
01281 B_data[BnumEl]=a_data[indexA];
01282 BnumEl++;
01283 }
01284 }
01285 }
01286 TMatrixDSparse *F=0;
01287 if(FnumEl) {
01288 #ifndef FORCE_EIGENVALUE_DECOMPOSITION
01289 F=CreateSparseMatrix(nF,nF,FnumEl,F_row,F_col,F_data);
01290 #endif
01291 }
01292 delete [] F_data;
01293 delete [] F_col;
01294 delete [] F_row;
01295 if(BnumEl) {
01296 B=CreateSparseMatrix(nF,nD2,BnumEl,B_row,B_col,B_data);
01297 }
01298 delete [] B_data;
01299 delete [] B_col;
01300 delete [] B_row;
01301
01302 if(B && D2inv) {
01303 minusBD2inv=MultiplyMSparseMSparse(B,D2inv);
01304 if(minusBD2inv) {
01305 Int_t mbd2_nMax=minusBD2inv->GetRowIndexArray()
01306 [minusBD2inv->GetNrows()];
01307 Double_t *mbd2_data=minusBD2inv->GetMatrixArray();
01308 for(Int_t i=0;i<mbd2_nMax;i++) {
01309 mbd2_data[i] = - mbd2_data[i];
01310 }
01311 }
01312 }
01313 if(minusBD2inv && F) {
01314 TMatrixDSparse *minusBD2invBt=
01315 MultiplyMSparseMSparseTranspVector(minusBD2inv,B,0);
01316 AddMSparse(F,1.,minusBD2invBt);
01317 DeleteMatrix(&minusBD2invBt);
01318 }
01319
01320 if(F) {
01321
01322 const Int_t *f_rows=F->GetRowIndexArray();
01323 const Int_t *f_cols=F->GetColIndexArray();
01324 const Double_t *f_data=F->GetMatrixArray();
01325
01326 TMatrixD c(nF,nF);
01327 Int_t nErrorF=0;
01328 for(Int_t i=0;i<nF;i++) {
01329 for(Int_t indexF=f_rows[i];indexF<f_rows[i+1];indexF++) {
01330 if(f_cols[indexF]>=i) c(f_cols[indexF],i)=f_data[indexF];
01331 }
01332
01333 Double_t c_ii=c(i,i);
01334 for(Int_t j=0;j<i;j++) {
01335 Double_t c_ij=c(i,j);
01336 c_ii -= c_ij*c_ij;
01337 }
01338 if(c_ii<=0.0) {
01339 nErrorF++;
01340 break;
01341 }
01342 c_ii=TMath::Sqrt(c_ii);
01343 c(i,i)=c_ii;
01344
01345 for(Int_t j=i+1;j<nF;j++) {
01346 Double_t c_ji=c(j,i);
01347 for(Int_t k=0;k<i;k++) {
01348 c_ji -= c(i,k)*c(j,k);
01349 }
01350 c(j,i) = c_ji/c_ii;
01351 }
01352 }
01353
01354 if(!nErrorF) {
01355 Double_t dmin=c(0,0);
01356 Double_t dmax=dmin;
01357 for(Int_t i=1;i<nF;i++) {
01358 dmin=TMath::Min(dmin,c(i,i));
01359 dmax=TMath::Max(dmax,c(i,i));
01360 }
01361 #ifdef DEBUG
01362 std::cout<<"dmin,dmax: "<<dmin<<" "<<dmax<<"\n";
01363 #endif
01364 if(dmin/dmax<epsilonF2*nF) nErrorF=2;
01365 }
01366 if(!nErrorF) {
01367
01368
01369 TMatrixD cinv(nF,nF);
01370 for(Int_t i=0;i<nF;i++) {
01371 cinv(i,i)=1./c(i,i);
01372 }
01373 for(Int_t i=0;i<nF;i++) {
01374 for(Int_t j=i+1;j<nF;j++) {
01375 Double_t tmp=-c(j,i)*cinv(i,i);
01376 for(Int_t k=i+1;k<j;k++) {
01377 tmp -= cinv(k,i)*c(j,k);
01378 }
01379 cinv(j,i)=tmp*cinv(j,j);
01380 }
01381 }
01382 TMatrixDSparse cInvSparse(cinv);
01383 Finv=MultiplyMSparseTranspMSparse
01384 (&cInvSparse,&cInvSparse);
01385 }
01386 DeleteMatrix(&F);
01387 }
01388 }
01389
01390
01391
01392
01393
01394
01395 Int_t rankD2Block=0;
01396 if(rNumEl>=0) {
01397 if( ((nD2==0)||D2inv) && ((nF==0)||Finv) ) {
01398
01399
01400
01401
01402
01403
01404
01405
01406 TMatrixDSparse *G=0;
01407 if(Finv && minusBD2inv) {
01408 G=MultiplyMSparseMSparse(Finv,minusBD2inv);
01409 }
01410 TMatrixDSparse *E=0;
01411 if(D2inv) E=new TMatrixDSparse(*D2inv);
01412 if(G && minusBD2inv) {
01413 TMatrixDSparse *minusBD2invTransG=
01414 MultiplyMSparseTranspMSparse(minusBD2inv,G);
01415 if(E) {
01416 AddMSparse(E,1.,minusBD2invTransG);
01417 DeleteMatrix(&minusBD2invTransG);
01418 } else {
01419 E=minusBD2invTransG;
01420 }
01421 }
01422
01423 if(E) {
01424 const Int_t *e_rows=E->GetRowIndexArray();
01425 const Int_t *e_cols=E->GetColIndexArray();
01426 const Double_t *e_data=E->GetMatrixArray();
01427 for(Int_t iE=0;iE<E->GetNrows();iE++) {
01428 Int_t iA=swap[iE+iDiagonal];
01429 for(Int_t indexE=e_rows[iE];indexE<e_rows[iE+1];++indexE) {
01430 Int_t jE=e_cols[indexE];
01431 Int_t jA=swap[jE+iDiagonal];
01432 rEl_col[rNumEl]=iA;
01433 rEl_row[rNumEl]=jA;
01434 rEl_data[rNumEl]=e_data[indexE];
01435 ++rNumEl;
01436 }
01437 }
01438 DeleteMatrix(&E);
01439 }
01440
01441 if(G) {
01442 const Int_t *g_rows=G->GetRowIndexArray();
01443 const Int_t *g_cols=G->GetColIndexArray();
01444 const Double_t *g_data=G->GetMatrixArray();
01445 for(Int_t iG=0;iG<G->GetNrows();iG++) {
01446 Int_t iA=swap[iG+iBlock];
01447 for(Int_t indexG=g_rows[iG];indexG<g_rows[iG+1];++indexG) {
01448 Int_t jG=g_cols[indexG];
01449 Int_t jA=swap[jG+iDiagonal];
01450
01451 rEl_col[rNumEl]=iA;
01452 rEl_row[rNumEl]=jA;
01453 rEl_data[rNumEl]=g_data[indexG];
01454 ++rNumEl;
01455
01456 rEl_col[rNumEl]=jA;
01457 rEl_row[rNumEl]=iA;
01458 rEl_data[rNumEl]=g_data[indexG];
01459 ++rNumEl;
01460 }
01461 }
01462 DeleteMatrix(&G);
01463 }
01464 if(Finv) {
01465
01466 const Int_t *finv_rows=Finv->GetRowIndexArray();
01467 const Int_t *finv_cols=Finv->GetColIndexArray();
01468 const Double_t *finv_data=Finv->GetMatrixArray();
01469 for(Int_t iF=0;iF<Finv->GetNrows();iF++) {
01470 Int_t iA=swap[iF+iBlock];
01471 for(Int_t indexF=finv_rows[iF];indexF<finv_rows[iF+1];++indexF) {
01472 Int_t jF=finv_cols[indexF];
01473 Int_t jA=swap[jF+iBlock];
01474 rEl_col[rNumEl]=iA;
01475 rEl_row[rNumEl]=jA;
01476 rEl_data[rNumEl]=finv_data[indexF];
01477 ++rNumEl;
01478 }
01479 }
01480 }
01481 rankD2Block=nD2+nF;
01482 } else if(!rankPtr) {
01483 rNumEl=-3;
01484 Fatal("InvertMSparseSymmPos",
01485 "non-trivial part has rank < %d, matrix can not be inverted",
01486 nF);
01487 } else {
01488
01489
01490 Int_t nEV=nD2+nF;
01491 Double_t epsilonEV2=fEpsMatrix ;
01492 Info("InvertMSparseSymmPos",
01493 "cholesky-decomposition failed, try eigenvalue analysis");
01494 #ifdef DEBUG
01495 std::cout<<"nEV="<<nEV<<" iDiagonal="<<iDiagonal<<"\n";
01496 #endif
01497 TMatrixDSym EV(nEV);
01498 for(Int_t i=0;i<nEV;i++) {
01499 Int_t iA=swap[i+iDiagonal];
01500 for(Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
01501 Int_t jA=a_cols[indexA];
01502 Int_t j0=swapBack[jA];
01503 if(j0>=iDiagonal) {
01504 Int_t j=j0-iDiagonal;
01505 #ifdef DEBUG
01506 if((i<0)||(j<0)||(i>=nEV)||(j>=nEV)) {
01507 std::cout<<" error "<<nEV<<" "<<i<<" "<<j<<"\n";
01508 }
01509 #endif
01510 if(!TMath::Finite(a_data[indexA])) {
01511 Fatal("InvertMSparseSymmPos",
01512 "non-finite number detected element %d %d\n",
01513 iA,jA);
01514 }
01515 EV(i,j)=a_data[indexA];
01516 }
01517 }
01518 }
01519
01520 TMatrixDSymEigen Eigen(EV);
01521 #ifdef DEBUG
01522 std::cout<<"Eigenvalues\n";
01523
01524 #endif
01525 Int_t errorEV=0;
01526 Double_t condition=1.0;
01527 if(Eigen.GetEigenValues()(0)<0.0) {
01528 errorEV=1;
01529 } else if(Eigen.GetEigenValues()(0)>0.0) {
01530 condition=Eigen.GetEigenValues()(nEV-1)/Eigen.GetEigenValues()(0);
01531 }
01532 if(condition<-epsilonEV2) {
01533 errorEV=2;
01534 }
01535 if(errorEV) {
01536 if(errorEV==1) {
01537 Error("InvertMSparseSymmPos",
01538 "Largest Eigenvalue is negative %f",
01539 Eigen.GetEigenValues()(0));
01540 } else {
01541 Error("InvertMSparseSymmPos",
01542 "Some Eigenvalues are negative (EV%d/EV0=%g epsilon=%g)",
01543 nEV-1,
01544 Eigen.GetEigenValues()(nEV-1)/Eigen.GetEigenValues()(0),
01545 epsilonEV2);
01546 }
01547 }
01548
01549 rankD2Block=0;
01550 Double_t setToZero=epsilonEV2*Eigen.GetEigenValues()(0);
01551 TMatrixD inverseEV(nEV,1);
01552 for(Int_t i=0;i<nEV;i++) {
01553 Double_t x=Eigen.GetEigenValues()(i);
01554 if(x>setToZero) {
01555 inverseEV(i,0)=1./x;
01556 ++rankD2Block;
01557 }
01558 }
01559 TMatrixDSparse V(Eigen.GetEigenVectors());
01560 TMatrixDSparse *VDVt=MultiplyMSparseMSparseTranspVector
01561 (&V,&V,&inverseEV);
01562
01563
01564 const Int_t *vdvt_rows=VDVt->GetRowIndexArray();
01565 const Int_t *vdvt_cols=VDVt->GetColIndexArray();
01566 const Double_t *vdvt_data=VDVt->GetMatrixArray();
01567 for(Int_t iVDVt=0;iVDVt<VDVt->GetNrows();iVDVt++) {
01568 Int_t iA=swap[iVDVt+iDiagonal];
01569 for(Int_t indexVDVt=vdvt_rows[iVDVt];
01570 indexVDVt<vdvt_rows[iVDVt+1];++indexVDVt) {
01571 Int_t jVDVt=vdvt_cols[indexVDVt];
01572 Int_t jA=swap[jVDVt+iDiagonal];
01573 rEl_col[rNumEl]=iA;
01574 rEl_row[rNumEl]=jA;
01575 rEl_data[rNumEl]=vdvt_data[indexVDVt];
01576 ++rNumEl;
01577 }
01578 }
01579 DeleteMatrix(&VDVt);
01580 }
01581 }
01582 if(rankPtr) {
01583 *rankPtr=rankD1+rankD2Block;
01584 }
01585
01586
01587 DeleteMatrix(&D2inv);
01588 DeleteMatrix(&Finv);
01589 DeleteMatrix(&B);
01590 DeleteMatrix(&minusBD2inv);
01591
01592 delete [] swap;
01593 delete [] swapBack;
01594 delete [] isZero;
01595
01596 TMatrixDSparse *r=(rNumEl>=0) ?
01597 CreateSparseMatrix(A->GetNrows(),A->GetNrows(),rNumEl,
01598 rEl_row,rEl_col,rEl_data) : 0;
01599 delete [] rEl_data;
01600 delete [] rEl_col;
01601 delete [] rEl_row;
01602
01603 #ifdef DEBUG_DETAIL
01604
01605 if(r) {
01606 TMatrixDSparse *Ar=MultiplyMSparseMSparse(A,r);
01607 TMatrixDSparse *ArA=MultiplyMSparseMSparse(Ar,A);
01608 TMatrixDSparse *rAr=MultiplyMSparseMSparse(r,Ar);
01609
01610 TMatrixD ar(*Ar);
01611 TMatrixD a(*A);
01612 TMatrixD ara(*ArA);
01613 TMatrixD R(*r);
01614 TMatrixD rar(*rAr);
01615
01616 DeleteMatrix(&Ar);
01617 DeleteMatrix(&ArA);
01618 DeleteMatrix(&rAr);
01619
01620 Double_t epsilonA2=fEpsMatrix ;
01621 for(Int_t i=0;i<a.GetNrows();i++) {
01622 for(Int_t j=0;j<a.GetNcols();j++) {
01623
01624 if(TMath::Abs(ar(i,j)-ar(j,i))>
01625 epsilonA2*(TMath::Abs(ar(i,j))+TMath::Abs(ar(j,i)))) {
01626 std::cout<<"Ar is not symmetric Ar("<<i<<","<<j<<")="<<ar(i,j)
01627 <<" Ar("<<j<<","<<i<<")="<<ar(j,i)<<"\n";
01628 }
01629
01630 if(TMath::Abs(ara(i,j)-a(i,j))>
01631 epsilonA2*(TMath::Abs(ara(i,j))+TMath::Abs(a(i,j)))) {
01632 std::cout<<"ArA is not equal A ArA("<<i<<","<<j<<")="<<ara(i,j)
01633 <<" A("<<i<<","<<j<<")="<<a(i,j)<<"\n";
01634 }
01635
01636 if(TMath::Abs(rar(i,j)-R(i,j))>
01637 epsilonA2*(TMath::Abs(rar(i,j))+TMath::Abs(R(i,j)))) {
01638 std::cout<<"rAr is not equal r rAr("<<i<<","<<j<<")="<<rar(i,j)
01639 <<" r("<<i<<","<<j<<")="<<R(i,j)<<"\n";
01640 }
01641 }
01642 }
01643 if(rankPtr) std::cout<<"rank="<<*rankPtr<<"\n";
01644 } else {
01645 std::cout<<"Matrix is not positive\n";
01646 }
01647 #endif
01648 return r;
01649
01650 }
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664 TString TUnfoldV17::GetOutputBinName(Int_t iBinX) const
01665 {
01666 return TString::Format("#%d",iBinX);
01667 }
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696 TUnfoldV17::TUnfoldV17(const TH2 *hist_A, EHistMap histmap, ERegMode regmode,
01697 EConstraint constraint)
01698 {
01699
01700
01701
01702
01703
01704 InitTUnfold();
01705 SetConstraint(constraint);
01706 Int_t nx0, nx, ny;
01707 if (histmap == kHistMapOutputHoriz) {
01708
01709 nx0 = hist_A->GetNbinsX()+2;
01710 ny = hist_A->GetNbinsY();
01711 } else {
01712
01713 nx0 = hist_A->GetNbinsY()+2;
01714 ny = hist_A->GetNbinsX();
01715 }
01716 nx = 0;
01717
01718
01719
01720 fSumOverY.Set(nx0);
01721 fXToHist.Set(nx0);
01722 fHistToX.Set(nx0);
01723 Int_t nonzeroA=0;
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734 Int_t nskipped=0;
01735 for (Int_t ix = 0; ix < nx0; ix++) {
01736
01737
01738 Double_t sum = 0.0;
01739 Int_t nonZeroY = 0;
01740 for (Int_t iy = 0; iy < ny; iy++) {
01741 Double_t z;
01742 if (histmap == kHistMapOutputHoriz) {
01743 z = hist_A->GetBinContent(ix, iy + 1);
01744 } else {
01745 z = hist_A->GetBinContent(iy + 1, ix);
01746 }
01747 if (z != 0.0) {
01748 nonzeroA++;
01749 nonZeroY++;
01750 sum += z;
01751 }
01752 }
01753
01754
01755 if (nonZeroY) {
01756
01757 fXToHist[nx] = ix;
01758 fHistToX[ix] = nx;
01759
01760
01761 fSumOverY[nx] = sum;
01762 if (histmap == kHistMapOutputHoriz) {
01763 fSumOverY[nx] +=
01764 hist_A->GetBinContent(ix, 0) +
01765 hist_A->GetBinContent(ix, ny + 1);
01766 } else {
01767 fSumOverY[nx] +=
01768 hist_A->GetBinContent(0, ix) +
01769 hist_A->GetBinContent(ny + 1, ix);
01770 }
01771 nx++;
01772 } else {
01773 nskipped++;
01774
01775 fHistToX[ix] = -1;
01776 }
01777 }
01778 Int_t underflowBin=0,overflowBin=0;
01779 for (Int_t ix = 0; ix < nx; ix++) {
01780 Int_t ibinx = fXToHist[ix];
01781 if(ibinx<1) underflowBin=1;
01782 if (histmap == kHistMapOutputHoriz) {
01783 if(ibinx>hist_A->GetNbinsX()) overflowBin=1;
01784 } else {
01785 if(ibinx>hist_A->GetNbinsY()) overflowBin=1;
01786 }
01787 }
01788 if(nskipped) {
01789 Int_t nprint=0;
01790 Int_t ixfirst=-1,ixlast=-1;
01791 TString binlist;
01792 int nDisconnected=0;
01793 for (Int_t ix = 0; ix < nx0; ix++) {
01794 if(fHistToX[ix]<0) {
01795 nprint++;
01796 if(ixlast<0) {
01797 binlist +=" ";
01798 binlist +=ix;
01799 ixfirst=ix;
01800 }
01801 ixlast=ix;
01802 nDisconnected++;
01803 }
01804 if(((fHistToX[ix]>=0)&&(ixlast>=0))||
01805 (nprint==nskipped)) {
01806 if(ixlast>ixfirst) {
01807 binlist += "-";
01808 binlist += ixlast;
01809 }
01810 ixfirst=-1;
01811 ixlast=-1;
01812 }
01813 if(nprint==nskipped) {
01814 break;
01815 }
01816 }
01817 if(nskipped==(2-underflowBin-overflowBin)) {
01818 Info("TUnfold","underflow and overflow bin "
01819 "do not depend on the input data");
01820 } else {
01821 Warning("TUnfold","%d output bins "
01822 "do not depend on the input data %s",nDisconnected,
01823 static_cast<const char *>(binlist));
01824 }
01825 }
01826
01827 fX0 = new TMatrixD(nx, 1, fSumOverY.GetArray());
01828
01829
01830 Int_t *rowA = new Int_t[nonzeroA];
01831 Int_t *colA = new Int_t[nonzeroA];
01832 Double_t *dataA = new Double_t[nonzeroA];
01833 Int_t index=0;
01834 for (Int_t iy = 0; iy < ny; iy++) {
01835 for (Int_t ix = 0; ix < nx; ix++) {
01836 Int_t ibinx = fXToHist[ix];
01837 Double_t z;
01838 if (histmap == kHistMapOutputHoriz) {
01839 z = hist_A->GetBinContent(ibinx, iy + 1);
01840 } else {
01841 z = hist_A->GetBinContent(iy + 1, ibinx);
01842 }
01843 if (z != 0.0) {
01844 rowA[index]=iy;
01845 colA[index] = ix;
01846 dataA[index] = z / fSumOverY[ix];
01847 index++;
01848 }
01849 }
01850 }
01851 if(underflowBin && overflowBin) {
01852 Info("TUnfold","%d input bins and %d output bins (includes 2 underflow/overflow bins)",ny,nx);
01853 } else if(underflowBin) {
01854 Info("TUnfold","%d input bins and %d output bins (includes 1 underflow bin)",ny,nx);
01855 } else if(overflowBin) {
01856 Info("TUnfold","%d input bins and %d output bins (includes 1 overflow bin)",ny,nx);
01857 } else {
01858 Info("TUnfold","%d input bins and %d output bins",ny,nx);
01859 }
01860 fA = CreateSparseMatrix(ny,nx,index,rowA,colA,dataA);
01861 if(ny<nx) {
01862 Error("TUnfold","too few (ny=%d) input bins for nx=%d output bins",ny,nx);
01863 } else if(ny==nx) {
01864 Warning("TUnfold","too few (ny=%d) input bins for nx=%d output bins",ny,nx);
01865 }
01866 delete[] rowA;
01867 delete[] colA;
01868 delete[] dataA;
01869
01870 fL = new TMatrixDSparse(1, GetNx());
01871 if (regmode != kRegModeNone) {
01872 Int_t nError=RegularizeBins(0, 1, nx0, regmode);
01873 if(nError>0) {
01874 if(nError>1) {
01875 Info("TUnfold",
01876 "%d regularisation conditions have been skipped",nError);
01877 } else {
01878 Info("TUnfold",
01879 "One regularisation condition has been skipped");
01880 }
01881 }
01882 }
01883 }
01884
01885
01886
01887
01888
01889
01890
01891
01892 void TUnfoldV17::SetBias(const TH1 *bias)
01893 {
01894 DeleteMatrix(&fX0);
01895 fX0 = new TMatrixD(GetNx(), 1);
01896 for (Int_t i = 0; i < GetNx(); i++) {
01897 (*fX0) (i, 0) = bias->GetBinContent(fXToHist[i]);
01898 }
01899 }
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914 Bool_t TUnfoldV17::AddRegularisationCondition
01915 (Int_t i0,Double_t f0,Int_t i1,Double_t f1,Int_t i2,Double_t f2)
01916 {
01917 Int_t indices[3];
01918 Double_t data[3];
01919 Int_t nEle=0;
01920
01921 if(i2>=0) {
01922 data[nEle]=f2;
01923 indices[nEle]=i2;
01924 nEle++;
01925 }
01926 if(i1>=0) {
01927 data[nEle]=f1;
01928 indices[nEle]=i1;
01929 nEle++;
01930 }
01931 if(i0>=0) {
01932 data[nEle]=f0;
01933 indices[nEle]=i0;
01934 nEle++;
01935 }
01936 return AddRegularisationCondition(nEle,indices,data);
01937 }
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951 Bool_t TUnfoldV17::AddRegularisationCondition
01952 (Int_t nEle,const Int_t *indices,const Double_t *rowData)
01953 {
01954 Bool_t r=kTRUE;
01955 const Int_t *l0_rows=fL->GetRowIndexArray();
01956 const Int_t *l0_cols=fL->GetColIndexArray();
01957 const Double_t *l0_data=fL->GetMatrixArray();
01958
01959 Int_t nF=l0_rows[fL->GetNrows()]+nEle;
01960 Int_t *l_row=new Int_t[nF];
01961 Int_t *l_col=new Int_t[nF];
01962 Double_t *l_data=new Double_t[nF];
01963
01964 nF=0;
01965 for(Int_t row=0;row<fL->GetNrows();row++) {
01966 for(Int_t k=l0_rows[row];k<l0_rows[row+1];k++) {
01967 l_row[nF]=row;
01968 l_col[nF]=l0_cols[k];
01969 l_data[nF]=l0_data[k];
01970 nF++;
01971 }
01972 }
01973
01974
01975 Int_t rowMax=0;
01976 if(nF>0) {
01977 rowMax=fL->GetNrows();
01978 }
01979
01980
01981 for(Int_t i=0;i<nEle;i++) {
01982 Int_t col=fHistToX[indices[i]];
01983 if(col<0) {
01984 r=kFALSE;
01985 break;
01986 }
01987 l_row[nF]=rowMax;
01988 l_col[nF]=col;
01989 l_data[nF]=rowData[i];
01990 nF++;
01991 }
01992
01993
01994 if(r) {
01995 DeleteMatrix(&fL);
01996 fL=CreateSparseMatrix(rowMax+1,GetNx(),nF,l_row,l_col,l_data);
01997 }
01998 delete [] l_row;
01999 delete [] l_col;
02000 delete [] l_data;
02001 return r;
02002 }
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022 Int_t TUnfoldV17::RegularizeSize(int bin, Double_t scale)
02023 {
02024
02025
02026
02027
02028
02029
02030 if(fRegMode==kRegModeNone) fRegMode=kRegModeSize;
02031 if(fRegMode!=kRegModeSize) fRegMode=kRegModeMixed;
02032
02033 return AddRegularisationCondition(bin,scale) ? 0 : 1;
02034 }
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056 Int_t TUnfoldV17::RegularizeDerivative(int left_bin, int right_bin,
02057 Double_t scale)
02058 {
02059
02060
02061
02062
02063
02064
02065
02066 if(fRegMode==kRegModeNone) fRegMode=kRegModeDerivative;
02067 if(fRegMode!=kRegModeDerivative) fRegMode=kRegModeMixed;
02068
02069 return AddRegularisationCondition(left_bin,-scale,right_bin,scale) ? 0 : 1;
02070 }
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095 Int_t TUnfoldV17::RegularizeCurvature(int left_bin, int center_bin,
02096 int right_bin,
02097 Double_t scale_left,
02098 Double_t scale_right)
02099 {
02100
02101
02102
02103
02104
02105
02106
02107
02108
02109 if(fRegMode==kRegModeNone) fRegMode=kRegModeCurvature;
02110 if(fRegMode!=kRegModeCurvature) fRegMode=kRegModeMixed;
02111
02112 return AddRegularisationCondition
02113 (left_bin,-scale_left,
02114 center_bin,scale_left+scale_right,
02115 right_bin,-scale_right)
02116 ? 0 : 1;
02117 }
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140 Int_t TUnfoldV17::RegularizeBins(int start, int step, int nbin,
02141 ERegMode regmode)
02142 {
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152 Int_t i0, i1, i2;
02153 i0 = start;
02154 i1 = i0 + step;
02155 i2 = i1 + step;
02156 Int_t nSkip = 0;
02157 Int_t nError= 0;
02158 if (regmode == kRegModeDerivative) {
02159 nSkip = 1;
02160 } else if (regmode == kRegModeCurvature) {
02161 nSkip = 2;
02162 } else if(regmode != kRegModeSize) {
02163 Error("RegularizeBins","regmode = %d is not valid",regmode);
02164 }
02165 for (Int_t i = nSkip; i < nbin; i++) {
02166 if (regmode == kRegModeSize) {
02167 nError += RegularizeSize(i0);
02168 } else if (regmode == kRegModeDerivative) {
02169 nError += RegularizeDerivative(i0, i1);
02170 } else if (regmode == kRegModeCurvature) {
02171 nError += RegularizeCurvature(i0, i1, i2);
02172 }
02173 i0 = i1;
02174 i1 = i2;
02175 i2 += step;
02176 }
02177 return nError;
02178 }
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191
02192
02193
02194
02195
02196
02197
02198
02199
02200
02201 Int_t TUnfoldV17::RegularizeBins2D(int start_bin, int step1, int nbin1,
02202 int step2, int nbin2, ERegMode regmode)
02203 {
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214 Int_t nError = 0;
02215 for (Int_t i1 = 0; i1 < nbin1; i1++) {
02216 nError += RegularizeBins(start_bin + step1 * i1, step2, nbin2, regmode);
02217 }
02218 for (Int_t i2 = 0; i2 < nbin2; i2++) {
02219 nError += RegularizeBins(start_bin + step2 * i2, step1, nbin1, regmode);
02220 }
02221 return nError;
02222 }
02223
02224
02225
02226
02227
02228
02229
02230
02231
02232 Double_t TUnfoldV17::DoUnfold(Double_t tau_reg,const TH1 *input,
02233 Double_t scaleBias)
02234 {
02235
02236
02237
02238
02239
02240
02241
02242
02243
02244
02245 SetInput(input,scaleBias);
02246 return DoUnfold(tau_reg);
02247 }
02248
02249
02250
02251
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264
02265
02266
02267
02268
02269
02270
02271 Int_t TUnfoldV17::SetInput(const TH1 *input, Double_t scaleBias,
02272 Double_t oneOverZeroError,const TH2 *hist_vyy,
02273 const TH2 *hist_vyy_inv)
02274 {
02275
02276
02277
02278
02279
02280
02281 DeleteMatrix(&fVyyInv);
02282 fNdf=0;
02283
02284 fBiasScale = scaleBias;
02285
02286
02287 ClearResults();
02288
02289
02290
02291
02292 Int_t *rowVyyN=new Int_t[GetNy()*GetNy()+1];
02293 Int_t *colVyyN=new Int_t[GetNy()*GetNy()+1];
02294 Double_t *dataVyyN=new Double_t[GetNy()*GetNy()+1];
02295
02296 Int_t *rowVyy1=new Int_t[GetNy()];
02297 Int_t *colVyy1=new Int_t[GetNy()];
02298 Double_t *dataVyy1=new Double_t[GetNy()];
02299 Double_t *dataVyyDiag=new Double_t[GetNy()];
02300
02301 Int_t nVarianceZero=0;
02302 Int_t nVarianceForced=0;
02303 Int_t nVyyN=0;
02304 Int_t nVyy1=0;
02305 for (Int_t iy = 0; iy < GetNy(); iy++) {
02306
02307 Double_t dy2;
02308 if(!hist_vyy) {
02309 Double_t dy = input->GetBinError(iy + 1);
02310 dy2=dy*dy;
02311 if (dy2 <= 0.0) {
02312 if(oneOverZeroError>0.0) {
02313 dy2 = 1./ ( oneOverZeroError*oneOverZeroError);
02314 nVarianceForced++;
02315 } else {
02316 nVarianceZero++;
02317 }
02318 }
02319 } else {
02320 dy2 = hist_vyy->GetBinContent(iy+1,iy+1);
02321 if (dy2 <= 0.0) {
02322 nVarianceZero++;
02323 }
02324 }
02325 rowVyyN[nVyyN] = iy;
02326 colVyyN[nVyyN] = iy;
02327 rowVyy1[nVyy1] = iy;
02328 colVyy1[nVyy1] = 0;
02329 dataVyyDiag[iy] = dy2;
02330 if(dy2>0.0) {
02331 dataVyyN[nVyyN++] = dy2;
02332 dataVyy1[nVyy1++] = dy2;
02333 }
02334 }
02335 if(hist_vyy) {
02336
02337 Int_t nOffDiagNonzero=0;
02338 for (Int_t iy = 0; iy < GetNy(); iy++) {
02339
02340 if(dataVyyDiag[iy]<=0.0) {
02341 for (Int_t jy = 0; jy < GetNy(); jy++) {
02342 if(hist_vyy->GetBinContent(iy+1,jy+1)!=0.0) {
02343 nOffDiagNonzero++;
02344 }
02345 }
02346 continue;
02347 }
02348 for (Int_t jy = 0; jy < GetNy(); jy++) {
02349
02350 if(iy==jy) continue;
02351
02352 if(dataVyyDiag[jy]<=0.0) continue;
02353
02354 rowVyyN[nVyyN] = iy;
02355 colVyyN[nVyyN] = jy;
02356 dataVyyN[nVyyN]= hist_vyy->GetBinContent(iy+1,jy+1);
02357 if(dataVyyN[nVyyN] == 0.0) continue;
02358 nVyyN ++;
02359 }
02360 }
02361 if(hist_vyy_inv) {
02362 Warning("SetInput",
02363 "inverse of input covariance is taken from user input");
02364 Int_t *rowVyyInv=new Int_t[GetNy()*GetNy()+1];
02365 Int_t *colVyyInv=new Int_t[GetNy()*GetNy()+1];
02366 Double_t *dataVyyInv=new Double_t[GetNy()*GetNy()+1];
02367 Int_t nVyyInv=0;
02368 for (Int_t iy = 0; iy < GetNy(); iy++) {
02369 for (Int_t jy = 0; jy < GetNy(); jy++) {
02370 rowVyyInv[nVyyInv] = iy;
02371 colVyyInv[nVyyInv] = jy;
02372 dataVyyInv[nVyyInv]= hist_vyy_inv->GetBinContent(iy+1,jy+1);
02373 if(dataVyyInv[nVyyInv] == 0.0) continue;
02374 nVyyInv ++;
02375 }
02376 }
02377 fVyyInv=CreateSparseMatrix
02378 (GetNy(),GetNy(),nVyyInv,rowVyyInv,colVyyInv,dataVyyInv);
02379 delete [] rowVyyInv;
02380 delete [] colVyyInv;
02381 delete [] dataVyyInv;
02382 } else {
02383 if(nOffDiagNonzero) {
02384 Error("SetInput",
02385 "input covariance has elements C(X,Y)!=0 where V(X)==0");
02386 }
02387 }
02388 }
02389 DeleteMatrix(&fVyy);
02390 fVyy = CreateSparseMatrix
02391 (GetNy(),GetNy(),nVyyN,rowVyyN,colVyyN,dataVyyN);
02392
02393 delete[] rowVyyN;
02394 delete[] colVyyN;
02395 delete[] dataVyyN;
02396
02397 TMatrixDSparse *vecV=CreateSparseMatrix
02398 (GetNy(),1,nVyy1,rowVyy1,colVyy1, dataVyy1);
02399
02400 delete[] rowVyy1;
02401 delete[] colVyy1;
02402 delete[] dataVyy1;
02403
02404
02405
02406 DeleteMatrix(&fY);
02407 fY = new TMatrixD(GetNy(), 1);
02408 for (Int_t i = 0; i < GetNy(); i++) {
02409 (*fY) (i, 0) = input->GetBinContent(i + 1);
02410 }
02411
02412 TMatrixDSparse *mAtV=MultiplyMSparseTranspMSparse(fA,vecV);
02413 DeleteMatrix(&vecV);
02414 Int_t nError2=0;
02415 for (Int_t i = 0; i <mAtV->GetNrows();i++) {
02416 if(mAtV->GetRowIndexArray()[i]==
02417 mAtV->GetRowIndexArray()[i+1]) {
02418 nError2 ++;
02419 }
02420 }
02421 if(nVarianceForced) {
02422 if(nVarianceForced>1) {
02423 Warning("SetInput","%d/%d input bins have zero error,"
02424 " 1/error set to %lf.",
02425 nVarianceForced,GetNy(),oneOverZeroError);
02426 } else {
02427 Warning("SetInput","One input bin has zero error,"
02428 " 1/error set to %lf.",oneOverZeroError);
02429 }
02430 }
02431 if(nVarianceZero) {
02432 if(nVarianceZero>1) {
02433 Warning("SetInput","%d/%d input bins have zero error,"
02434 " and are ignored.",nVarianceZero,GetNy());
02435 } else {
02436 Warning("SetInput","One input bin has zero error,"
02437 " and is ignored.");
02438 }
02439 fIgnoredBins=nVarianceZero;
02440 }
02441 if(nError2>0) {
02442
02443 if(oneOverZeroError<=0.0) {
02444
02445
02446 for (Int_t col = 0; col <mAtV->GetNrows();col++) {
02447 if(mAtV->GetRowIndexArray()[col]==
02448 mAtV->GetRowIndexArray()[col+1]) {
02449 TString binlist("no data to constrain output bin ");
02450 binlist += GetOutputBinName(fXToHist[col]);
02451
02452
02453
02454
02455
02456
02457
02458
02459
02460 Warning("SetInput",binlist);
02461 }
02462 }
02463 }
02464 if(nError2>1) {
02465 Error("SetInput","%d/%d output bins are not constrained by any data.",
02466 nError2,mAtV->GetNrows());
02467 } else {
02468 Error("SetInput","One output bins is not constrained by any data.");
02469 }
02470 }
02471 DeleteMatrix(&mAtV);
02472
02473 delete[] dataVyyDiag;
02474
02475 return nVarianceForced+nVarianceZero+10000*nError2;
02476 }
02477
02478
02479
02480
02481
02482
02483
02484
02485 Double_t TUnfoldV17::DoUnfold(Double_t tau)
02486 {
02487
02488
02489
02490
02491
02492
02493
02494
02495
02496 fTauSquared=tau*tau;
02497 return DoUnfold();
02498 }
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526 Int_t TUnfoldV17::ScanLcurve(Int_t nPoint,
02527 Double_t tauMin,Double_t tauMax,
02528 TGraph **lCurve,TSpline **logTauX,
02529 TSpline **logTauY,TSpline **logTauCurvature)
02530 {
02531 typedef std::map<Double_t,std::pair<Double_t,Double_t> > XYtau_t;
02532 XYtau_t curve;
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556 if((tauMin<=0)||(tauMax<=0.0)||(tauMin>=tauMax)) {
02557
02558
02559
02560
02561
02562 DoUnfold(0.0);
02563
02564
02565 if(GetNdf()<=0) {
02566 Error("ScanLcurve","too few input bins, NDF<=0 %d",GetNdf());
02567 }
02568
02569 Double_t x0=GetLcurveX();
02570 Double_t y0=GetLcurveY();
02571 Info("ScanLcurve","logtau=-Infinity X=%lf Y=%lf",x0,y0);
02572 if(!finite(x0)) {
02573 Fatal("ScanLcurve","problem (too few input bins?) X=%f",x0);
02574 }
02575 if(!finite(y0)) {
02576 Fatal("ScanLcurve","problem (missing regularisation?) Y=%f",y0);
02577 }
02578 {
02579
02580 Double_t logTau=
02581 0.5*(TMath::Log10(fChi2A+3.*TMath::Sqrt(GetNdf()+1.0))
02582 -GetLcurveY());
02583 DoUnfold(TMath::Power(10.,logTau));
02584 if((!finite(GetLcurveX())) ||(!finite(GetLcurveY()))) {
02585 Fatal("ScanLcurve","problem (missing regularisation?) X=%f Y=%f",
02586 GetLcurveX(),GetLcurveY());
02587 }
02588 curve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
02589 Info("ScanLcurve","logtau=%lf X=%lf Y=%lf",
02590 logTau,GetLcurveX(),GetLcurveY());
02591 }
02592 if((*curve.begin()).second.first<x0) {
02593
02594
02595
02596
02597
02598
02599 do {
02600 x0=GetLcurveX();
02601 Double_t logTau=(*curve.begin()).first-0.5;
02602 DoUnfold(TMath::Power(10.,logTau));
02603 if((!finite(GetLcurveX())) ||(!finite(GetLcurveY()))) {
02604 Fatal("ScanLcurve","problem (missing regularisation?) X=%f Y=%f",
02605 GetLcurveX(),GetLcurveY());
02606 }
02607 curve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
02608 Info("ScanLcurve","logtau=%lf X=%lf Y=%lf",
02609 logTau,GetLcurveX(),GetLcurveY());
02610 }
02611 while(((int)curve.size()<(nPoint-1)/2)&&
02612 ((*curve.begin()).second.first<x0));
02613 } else {
02614
02615
02616
02617
02618
02619 while(((int)curve.size()<nPoint-1)&&
02620 (((*curve.begin()).second.first-x0>0.00432)||
02621 ((*curve.begin()).second.second-y0>0.00432)||
02622 (curve.size()<2))) {
02623 Double_t logTau=(*curve.begin()).first-0.5;
02624 DoUnfold(TMath::Power(10.,logTau));
02625 if((!finite(GetLcurveX())) ||(!finite(GetLcurveY()))) {
02626 Fatal("ScanLcurve","problem (missing regularisation?) X=%f Y=%f",
02627 GetLcurveX(),GetLcurveY());
02628 }
02629 curve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
02630 Info("ScanLcurve","logtau=%lf X=%lf Y=%lf",
02631 logTau,GetLcurveX(),GetLcurveY());
02632 }
02633 }
02634 } else {
02635 Double_t logTauMin=TMath::Log10(tauMin);
02636 Double_t logTauMax=TMath::Log10(tauMax);
02637 if(nPoint>1) {
02638
02639 DoUnfold(TMath::Power(10.,logTauMax));
02640 if((!finite(GetLcurveX())) ||(!finite(GetLcurveY()))) {
02641 Fatal("ScanLcurve","problem (missing regularisation?) X=%f Y=%f",
02642 GetLcurveX(),GetLcurveY());
02643 }
02644 Info("ScanLcurve","logtau=%lf X=%lf Y=%lf",
02645 logTauMax,GetLcurveX(),GetLcurveY());
02646 curve[logTauMax]=std::make_pair(GetLcurveX(),GetLcurveY());
02647 }
02648
02649 DoUnfold(TMath::Power(10.,logTauMin));
02650 if((!finite(GetLcurveX())) ||(!finite(GetLcurveY()))) {
02651 Fatal("ScanLcurve","problem (missing regularisation?) X=%f Y=%f",
02652 GetLcurveX(),GetLcurveY());
02653 }
02654 Info("ScanLcurve","logtau=%lf X=%lf Y=%lf",
02655 logTauMin,GetLcurveX(),GetLcurveY());
02656 curve[logTauMin]=std::make_pair(GetLcurveX(),GetLcurveY());
02657 }
02658
02659
02660
02661
02662
02663
02664 while(int(curve.size())<nPoint-1) {
02665
02666
02667 XYtau_t::const_iterator i0,i1;
02668 i0=curve.begin();
02669 i1=i0;
02670 Double_t logTau=(*i0).first;
02671 Double_t distMax=0.0;
02672 for(i1++;i1!=curve.end();i1++) {
02673 const std::pair<Double_t,Double_t> &xy0=(*i0).second;
02674 const std::pair<Double_t,Double_t> &xy1=(*i1).second;
02675 Double_t dx=xy1.first-xy0.first;
02676 Double_t dy=xy1.second-xy0.second;
02677 Double_t d=TMath::Sqrt(dx*dx+dy*dy);
02678 if(d>=distMax) {
02679 distMax=d;
02680 logTau=0.5*((*i0).first+(*i1).first);
02681 }
02682 i0=i1;
02683 }
02684 DoUnfold(TMath::Power(10.,logTau));
02685 if((!finite(GetLcurveX())) ||(!finite(GetLcurveY()))) {
02686 Fatal("ScanLcurve","problem (missing regularisation?) X=%f Y=%f",
02687 GetLcurveX(),GetLcurveY());
02688 }
02689 Info("ScanLcurve","logtau=%lf X=%lf Y=%lf",logTau,GetLcurveX(),GetLcurveY());
02690 curve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
02691 }
02692
02693
02694
02695
02696
02697
02698 XYtau_t::const_iterator i0,i1;
02699 i0=curve.begin();
02700 i1=i0;
02701 i1++;
02702 Double_t logTauFin=(*i0).first;
02703 if( ((int)curve.size())<nPoint) {
02704
02705 Double_t *cTi=new Double_t[curve.size()-1];
02706 Double_t *cCi=new Double_t[curve.size()-1];
02707 Int_t n=0;
02708 {
02709 Double_t *lXi=new Double_t[curve.size()];
02710 Double_t *lYi=new Double_t[curve.size()];
02711 Double_t *lTi=new Double_t[curve.size()];
02712 for( XYtau_t::const_iterator i=curve.begin();i!=curve.end();i++) {
02713 lXi[n]=(*i).second.first;
02714 lYi[n]=(*i).second.second;
02715 lTi[n]=(*i).first;
02716 n++;
02717 }
02718 TSpline3 *splineX=new TSpline3("x vs tau",lTi,lXi,n);
02719 TSpline3 *splineY=new TSpline3("y vs tau",lTi,lYi,n);
02720
02721
02722 for(Int_t i=0;i<n-1;i++) {
02723 Double_t ltau,xy,bi,ci,di;
02724 splineX->GetCoeff(i,ltau,xy,bi,ci,di);
02725 Double_t tauBar=0.5*(lTi[i]+lTi[i+1]);
02726 Double_t dTau=0.5*(lTi[i+1]-lTi[i]);
02727 Double_t dx1=bi+dTau*(2.*ci+3.*di*dTau);
02728 Double_t dx2=2.*ci+6.*di*dTau;
02729 splineY->GetCoeff(i,ltau,xy,bi,ci,di);
02730 Double_t dy1=bi+dTau*(2.*ci+3.*di*dTau);
02731 Double_t dy2=2.*ci+6.*di*dTau;
02732 cTi[i]=tauBar;
02733 cCi[i]=(dy2*dx1-dy1*dx2)/TMath::Power(dx1*dx1+dy1*dy1,1.5);
02734 }
02735 delete splineX;
02736 delete splineY;
02737 delete[] lXi;
02738 delete[] lYi;
02739 delete[] lTi;
02740 }
02741
02742 TSpline3 *splineC=new TSpline3("L curve curvature",cTi,cCi,n-1);
02743
02744
02745
02746
02747
02748 Int_t iskip=0;
02749 if(n>4) iskip=1;
02750 if(n>7) iskip=2;
02751 Double_t cCmax=cCi[iskip];
02752 Double_t cTmax=cTi[iskip];
02753 for(Int_t i=iskip;i<n-2-iskip;i++) {
02754
02755
02756 Double_t xMax=cTi[i+1];
02757 Double_t yMax=cCi[i+1];
02758 if(cCi[i]>yMax) {
02759 yMax=cCi[i];
02760 xMax=cTi[i];
02761 }
02762
02763
02764
02765 Double_t x,y,b,c,d;
02766 splineC->GetCoeff(i,x,y,b,c,d);
02767
02768 Double_t m_p_half=-c/(3.*d);
02769 Double_t q=b/(3.*d);
02770 Double_t discr=m_p_half*m_p_half-q;
02771 if(discr>=0.0) {
02772
02773 discr=TMath::Sqrt(discr);
02774 Double_t xx;
02775 if(m_p_half>0.0) {
02776 xx = m_p_half + discr;
02777 } else {
02778 xx = m_p_half - discr;
02779 }
02780 Double_t dx=cTi[i+1]-x;
02781
02782 if((xx>0.0)&&(xx<dx)) {
02783 y=splineC->Eval(x+xx);
02784 if(y>yMax) {
02785 yMax=y;
02786 xMax=x+xx;
02787 }
02788 }
02789
02790 if(xx !=0.0) {
02791 xx= q/xx;
02792 } else {
02793 xx=0.0;
02794 }
02795
02796 if((xx>0.0)&&(xx<dx)) {
02797 y=splineC->Eval(x+xx);
02798 if(y>yMax) {
02799 yMax=y;
02800 xMax=x+xx;
02801 }
02802 }
02803 }
02804
02805 if(yMax>cCmax) {
02806 cCmax=yMax;
02807 cTmax=xMax;
02808 }
02809 }
02810 if(logTauCurvature) {
02811 *logTauCurvature=splineC;
02812 } else {
02813 delete splineC;
02814 }
02815 delete[] cTi;
02816 delete[] cCi;
02817 logTauFin=cTmax;
02818 DoUnfold(TMath::Power(10.,logTauFin));
02819 if((!finite(GetLcurveX())) ||(!finite(GetLcurveY()))) {
02820 Fatal("ScanLcurve","problem (missing regularisation?) X=%f Y=%f",
02821 GetLcurveX(),GetLcurveY());
02822 }
02823 Info("ScanLcurve","Result logtau=%lf X=%lf Y=%lf",
02824 logTauFin,GetLcurveX(),GetLcurveY());
02825 curve[logTauFin]=std::make_pair(GetLcurveX(),GetLcurveY());
02826 }
02827
02828
02829
02830
02831
02832
02833 Int_t bestChoice=-1;
02834 if(curve.size()>0) {
02835 Double_t *x=new Double_t[curve.size()];
02836 Double_t *y=new Double_t[curve.size()];
02837 Double_t *logT=new Double_t[curve.size()];
02838 int n=0;
02839 for( XYtau_t::const_iterator i=curve.begin();i!=curve.end();i++) {
02840 if(logTauFin==(*i).first) {
02841 bestChoice=n;
02842 }
02843 x[n]=(*i).second.first;
02844 y[n]=(*i).second.second;
02845 logT[n]=(*i).first;
02846 n++;
02847 }
02848 if(lCurve) {
02849 (*lCurve)=new TGraph(n,x,y);
02850 (*lCurve)->SetTitle("L curve");
02851 }
02852 if(logTauX) (*logTauX)=new TSpline3("log(chi**2)%log(tau)",logT,x,n);
02853 if(logTauY) (*logTauY)=new TSpline3("log(reg.cond)%log(tau)",logT,y,n);
02854 delete[] x;
02855 delete[] y;
02856 delete[] logT;
02857 }
02858
02859 return bestChoice;
02860 }
02861
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877 void TUnfoldV17::GetNormalisationVector(TH1 *out,const Int_t *binMap) const
02878 {
02879
02880 ClearHistogram(out);
02881 for (Int_t i = 0; i < GetNx(); i++) {
02882 int dest=binMap ? binMap[fXToHist[i]] : fXToHist[i];
02883 if(dest>=0) {
02884 out->SetBinContent(dest, fSumOverY[i] + out->GetBinContent(dest));
02885 }
02886 }
02887 }
02888
02889
02890
02891
02892
02893
02894
02895
02896
02897
02898
02899
02900
02901
02902 void TUnfoldV17::GetBias(TH1 *out,const Int_t *binMap) const
02903 {
02904
02905 ClearHistogram(out);
02906 for (Int_t i = 0; i < GetNx(); i++) {
02907 int dest=binMap ? binMap[fXToHist[i]] : fXToHist[i];
02908 if(dest>=0) {
02909 out->SetBinContent(dest, fBiasScale*(*fX0) (i, 0) +
02910 out->GetBinContent(dest));
02911 }
02912 }
02913 }
02914
02915
02916
02917
02918
02919
02920
02921
02922
02923
02924
02925
02926
02927
02928
02929 void TUnfoldV17::GetFoldedOutput(TH1 *out,const Int_t *binMap) const
02930 {
02931 ClearHistogram(out);
02932
02933 TMatrixDSparse *AVxx=MultiplyMSparseMSparse(fA,fVxx);
02934
02935 const Int_t *rows_A=fA->GetRowIndexArray();
02936 const Int_t *cols_A=fA->GetColIndexArray();
02937 const Double_t *data_A=fA->GetMatrixArray();
02938 const Int_t *rows_AVxx=AVxx->GetRowIndexArray();
02939 const Int_t *cols_AVxx=AVxx->GetColIndexArray();
02940 const Double_t *data_AVxx=AVxx->GetMatrixArray();
02941
02942 for (Int_t i = 0; i < GetNy(); i++) {
02943 Int_t destI = binMap ? binMap[i+1] : i+1;
02944 if(destI<0) continue;
02945
02946 out->SetBinContent(destI, (*fAx) (i, 0)+ out->GetBinContent(destI));
02947 Double_t e2=0.0;
02948 Int_t index_a=rows_A[i];
02949 Int_t index_av=rows_AVxx[i];
02950 while((index_a<rows_A[i+1])&&(index_av<rows_AVxx[i+1])) {
02951 Int_t j_a=cols_A[index_a];
02952 Int_t j_av=cols_AVxx[index_av];
02953 if(j_a<j_av) {
02954 index_a++;
02955 } else if(j_a>j_av) {
02956 index_av++;
02957 } else {
02958 e2 += data_AVxx[index_av] * data_A[index_a];
02959 index_a++;
02960 index_av++;
02961 }
02962 }
02963 out->SetBinError(destI,TMath::Sqrt(e2));
02964 }
02965 DeleteMatrix(&AVxx);
02966 }
02967
02968
02969
02970
02971
02972
02973
02974
02975
02976
02977 void TUnfoldV17::GetProbabilityMatrix(TH2 *A,EHistMap histmap) const
02978 {
02979
02980
02981
02982 const Int_t *rows_A=fA->GetRowIndexArray();
02983 const Int_t *cols_A=fA->GetColIndexArray();
02984 const Double_t *data_A=fA->GetMatrixArray();
02985 for (Int_t iy = 0; iy <fA->GetNrows(); iy++) {
02986 for(Int_t indexA=rows_A[iy];indexA<rows_A[iy+1];indexA++) {
02987 Int_t ix = cols_A[indexA];
02988 Int_t ih=fXToHist[ix];
02989 if (histmap == kHistMapOutputHoriz) {
02990 A->SetBinContent(ih, iy+1,data_A[indexA]);
02991 } else {
02992 A->SetBinContent(iy+1, ih,data_A[indexA]);
02993 }
02994 }
02995 }
02996 }
02997
02998
02999
03000
03001
03002
03003
03004
03005
03006
03007
03008
03009
03010
03011
03012
03013 void TUnfoldV17::GetInput(TH1 *out,const Int_t *binMap) const
03014 {
03015 ClearHistogram(out);
03016
03017 const Int_t *rows_Vyy=fVyy->GetRowIndexArray();
03018 const Int_t *cols_Vyy=fVyy->GetColIndexArray();
03019 const Double_t *data_Vyy=fVyy->GetMatrixArray();
03020
03021 for (Int_t i = 0; i < GetNy(); i++) {
03022 Int_t destI=binMap ? binMap[i+1] : i+1;
03023 if(destI<0) continue;
03024
03025 out->SetBinContent(destI, (*fY) (i, 0)+out->GetBinContent(destI));
03026
03027 Double_t e=0.0;
03028 for(int index=rows_Vyy[i];index<rows_Vyy[i+1];index++) {
03029 if(cols_Vyy[index]==i) {
03030 e=TMath::Sqrt(data_Vyy[index]);
03031 }
03032 }
03033 out->SetBinError(destI, e);
03034 }
03035 }
03036
03037
03038
03039
03040
03041
03042 void TUnfoldV17::GetInputInverseEmatrix(TH2 *out)
03043 {
03044
03045
03046 if(!fVyyInv) {
03047 Int_t rank=0;
03048 fVyyInv=InvertMSparseSymmPos(fVyy,&rank);
03049
03050 fNdf = rank-GetNpar();
03051
03052 if(rank<GetNy()-fIgnoredBins) {
03053 Warning("GetInputInverseEmatrix","input covariance matrix has rank %d expect %d",
03054 rank,GetNy());
03055 }
03056 if(fNdf<0) {
03057 Error("GetInputInverseEmatrix","number of parameters %d > %d (rank of input covariance). Problem can not be solved",GetNpar(),rank);
03058 } else if(fNdf==0) {
03059 Warning("GetInputInverseEmatrix","number of parameters %d = input rank %d. Problem is ill posed",GetNpar(),rank);
03060 }
03061 }
03062 if(out) {
03063
03064 const Int_t *rows_VyyInv=fVyyInv->GetRowIndexArray();
03065 const Int_t *cols_VyyInv=fVyyInv->GetColIndexArray();
03066 const Double_t *data_VyyInv=fVyyInv->GetMatrixArray();
03067
03068 for(int i=0;i<=out->GetNbinsX()+1;i++) {
03069 for(int j=0;j<=out->GetNbinsY()+1;j++) {
03070 out->SetBinContent(i,j,0.);
03071 }
03072 }
03073
03074 for (Int_t i = 0; i < fVyyInv->GetNrows(); i++) {
03075 for(int index=rows_VyyInv[i];index<rows_VyyInv[i+1];index++) {
03076 Int_t j=cols_VyyInv[index];
03077 out->SetBinContent(i+1,j+1,data_VyyInv[index]);
03078 }
03079 }
03080 }
03081 }
03082
03083
03084
03085
03086
03087
03088
03089
03090
03091
03092
03093
03094
03095 void TUnfoldV17::GetLsquared(TH2 *out) const
03096 {
03097
03098
03099
03100 TMatrixDSparse *lSquared=MultiplyMSparseTranspMSparse(fL,fL);
03101
03102 const Int_t *rows=lSquared->GetRowIndexArray();
03103 const Int_t *cols=lSquared->GetColIndexArray();
03104 const Double_t *data=lSquared->GetMatrixArray();
03105 for (Int_t i = 0; i < GetNx(); i++) {
03106 for (Int_t cindex = rows[i]; cindex < rows[i+1]; cindex++) {
03107 Int_t j=cols[cindex];
03108 out->SetBinContent(fXToHist[i], fXToHist[j],data[cindex]);
03109 }
03110 }
03111 DeleteMatrix(&lSquared);
03112 }
03113
03114
03115
03116
03117
03118
03119
03120 Int_t TUnfoldV17::GetNr(void) const {
03121 return fL->GetNrows();
03122 }
03123
03124
03125
03126
03127
03128
03129
03130
03131
03132
03133
03134
03135 void TUnfoldV17::GetL(TH2 *out) const
03136 {
03137
03138 const Int_t *rows=fL->GetRowIndexArray();
03139 const Int_t *cols=fL->GetColIndexArray();
03140 const Double_t *data=fL->GetMatrixArray();
03141 for (Int_t row = 0; row < GetNr(); row++) {
03142 for (Int_t cindex = rows[row]; cindex < rows[row+1]; cindex++) {
03143 Int_t col=cols[cindex];
03144 Int_t indexH=fXToHist[col];
03145 out->SetBinContent(indexH,row+1,data[cindex]);
03146 }
03147 }
03148 }
03149
03150
03151
03152
03153
03154
03155 void TUnfoldV17::SetConstraint(EConstraint constraint)
03156 {
03157
03158 if(fConstraint !=constraint) ClearResults();
03159 fConstraint=constraint;
03160 Info("SetConstraint","fConstraint=%d",fConstraint);
03161 }
03162
03163
03164
03165
03166
03167 Double_t TUnfoldV17::GetTau(void) const
03168 {
03169
03170 return TMath::Sqrt(fTauSquared);
03171 }
03172
03173
03174
03175 Double_t TUnfoldV17::GetChi2L(void) const
03176 {
03177
03178 return fLXsquared*fTauSquared;
03179 }
03180
03181
03182
03183
03184
03185
03186 Int_t TUnfoldV17::GetNpar(void) const
03187 {
03188 return GetNx();
03189 }
03190
03191
03192
03193
03194
03195 Double_t TUnfoldV17::GetLcurveX(void) const
03196 {
03197 return TMath::Log10(fChi2A);
03198 }
03199
03200
03201
03202
03203
03204 Double_t TUnfoldV17::GetLcurveY(void) const
03205 {
03206 return TMath::Log10(fLXsquared);
03207 }
03208
03209
03210
03211
03212
03213
03214
03215
03216
03217
03218
03219
03220
03221
03222
03223
03224
03225
03226
03227
03228
03229
03230
03231
03232
03233 void TUnfoldV17::GetOutput(TH1 *output,const Int_t *binMap) const
03234 {
03235
03236
03237
03238
03239
03240
03241
03242
03243
03244 ClearHistogram(output);
03245
03246
03247
03248
03249
03250
03251
03252
03253 std::map<Int_t,Double_t> e2;
03254
03255 const Int_t *rows_Vxx=fVxx->GetRowIndexArray();
03256 const Int_t *cols_Vxx=fVxx->GetColIndexArray();
03257 const Double_t *data_Vxx=fVxx->GetMatrixArray();
03258
03259 Int_t binMapSize = fHistToX.GetSize();
03260 for(Int_t i=0;i<binMapSize;i++) {
03261 Int_t destBinI=binMap ? binMap[i] : i;
03262 Int_t srcBinI=fHistToX[i];
03263 if((destBinI>=0)&&(srcBinI>=0)) {
03264 output->SetBinContent
03265 (destBinI, (*fX)(srcBinI,0)+ output->GetBinContent(destBinI));
03266
03267
03268
03269
03270 Int_t j=0;
03271 Int_t index_vxx=rows_Vxx[srcBinI];
03272
03273 while((j<binMapSize)&&(index_vxx<rows_Vxx[srcBinI+1])) {
03274 Int_t destBinJ=binMap ? binMap[j] : j;
03275 if(destBinI!=destBinJ) {
03276
03277 j++;
03278 } else {
03279 Int_t srcBinJ=fHistToX[j];
03280 if(srcBinJ<0) {
03281
03282 j++;
03283 } else {
03284 if(cols_Vxx[index_vxx]<srcBinJ) {
03285
03286 index_vxx++;
03287 } else if(cols_Vxx[index_vxx]>srcBinJ) {
03288
03289 j++;
03290 } else {
03291
03292 e2[destBinI] += data_Vxx[index_vxx];
03293 j++;
03294 index_vxx++;
03295 }
03296 }
03297 }
03298 }
03299
03300 }
03301 }
03302 for(std::map<Int_t,Double_t>::const_iterator i=e2.begin();
03303 i!=e2.end();i++) {
03304
03305 output->SetBinError((*i).first,TMath::Sqrt((*i).second));
03306 }
03307 }
03308
03309
03310
03311
03312
03313
03314
03315
03316
03317
03318
03319
03320
03321
03322
03323 void TUnfoldV17::ErrorMatrixToHist(TH2 *ematrix,const TMatrixDSparse *emat,
03324 const Int_t *binMap,Bool_t doClear) const
03325 {
03326 Int_t nbin=ematrix->GetNbinsX();
03327 if(doClear) {
03328 for(Int_t i=0;i<nbin+2;i++) {
03329 for(Int_t j=0;j<nbin+2;j++) {
03330 ematrix->SetBinContent(i,j,0.0);
03331 ematrix->SetBinError(i,j,0.0);
03332 }
03333 }
03334 }
03335
03336 if(emat) {
03337 const Int_t *rows_emat=emat->GetRowIndexArray();
03338 const Int_t *cols_emat=emat->GetColIndexArray();
03339 const Double_t *data_emat=emat->GetMatrixArray();
03340
03341 Int_t binMapSize = fHistToX.GetSize();
03342 for(Int_t i=0;i<binMapSize;i++) {
03343 Int_t destBinI=binMap ? binMap[i] : i;
03344 Int_t srcBinI=fHistToX[i];
03345 if((destBinI>=0)&&(destBinI<nbin+2)&&(srcBinI>=0)) {
03346
03347
03348
03349
03350 Int_t j=0;
03351 Int_t index_vxx=rows_emat[srcBinI];
03352 while((j<binMapSize)&&(index_vxx<rows_emat[srcBinI+1])) {
03353 Int_t destBinJ=binMap ? binMap[j] : j;
03354 Int_t srcBinJ=fHistToX[j];
03355 if((destBinJ<0)||(destBinJ>=nbin+2)||(srcBinJ<0)) {
03356
03357 j++;
03358 } else {
03359 if(cols_emat[index_vxx]<srcBinJ) {
03360
03361 index_vxx++;
03362 } else if(cols_emat[index_vxx]>srcBinJ) {
03363
03364 j++;
03365 } else {
03366
03367 Double_t e2= ematrix->GetBinContent(destBinI,destBinJ)
03368 + data_emat[index_vxx];
03369 ematrix->SetBinContent(destBinI,destBinJ,e2);
03370 j++;
03371 index_vxx++;
03372 }
03373 }
03374 }
03375 }
03376 }
03377 }
03378 }
03379
03380
03381
03382
03383
03384
03385
03386
03387
03388
03389
03390 void TUnfoldV17::GetEmatrix(TH2 *ematrix,const Int_t *binMap) const
03391 {
03392 ErrorMatrixToHist(ematrix,fVxx,binMap,kTRUE);
03393 }
03394
03395
03396
03397
03398
03399
03400
03401
03402
03403
03404
03405 void TUnfoldV17::GetRhoIJ(TH2 *rhoij,const Int_t *binMap) const
03406 {
03407 GetEmatrix(rhoij,binMap);
03408 Int_t nbin=rhoij->GetNbinsX();
03409 Double_t *e=new Double_t[nbin+2];
03410 for(Int_t i=0;i<nbin+2;i++) {
03411 e[i]=TMath::Sqrt(rhoij->GetBinContent(i,i));
03412 }
03413 for(Int_t i=0;i<nbin+2;i++) {
03414 for(Int_t j=0;j<nbin+2;j++) {
03415 if((e[i]>0.0)&&(e[j]>0.0)) {
03416 rhoij->SetBinContent(i,j,rhoij->GetBinContent(i,j)/e[i]/e[j]);
03417 } else {
03418 rhoij->SetBinContent(i,j,0.0);
03419 }
03420 }
03421 }
03422 delete[] e;
03423 }
03424
03425
03426
03427
03428
03429
03430
03431
03432
03433
03434
03435
03436
03437
03438
03439
03440
03441
03442
03443
03444
03445
03446
03447
03448 Double_t TUnfoldV17::GetRhoI(TH1 *rhoi,const Int_t *binMap,TH2 *invEmat) const
03449 {
03450 ClearHistogram(rhoi,-1.);
03451
03452 if(binMap) {
03453
03454
03455
03456 return GetRhoIFromMatrix(rhoi,fVxx,binMap,invEmat);
03457 } else {
03458 Double_t rhoMax=0.0;
03459
03460 const Int_t *rows_Vxx=fVxx->GetRowIndexArray();
03461 const Int_t *cols_Vxx=fVxx->GetColIndexArray();
03462 const Double_t *data_Vxx=fVxx->GetMatrixArray();
03463
03464 const Int_t *rows_VxxInv=fVxxInv->GetRowIndexArray();
03465 const Int_t *cols_VxxInv=fVxxInv->GetColIndexArray();
03466 const Double_t *data_VxxInv=fVxxInv->GetMatrixArray();
03467
03468 for(Int_t i=0;i<GetNx();i++) {
03469 Int_t destI=fXToHist[i];
03470 Double_t e_ii=0.0,einv_ii=0.0;
03471 for(int index_vxx=rows_Vxx[i];index_vxx<rows_Vxx[i+1];index_vxx++) {
03472 if(cols_Vxx[index_vxx]==i) {
03473 e_ii=data_Vxx[index_vxx];
03474 break;
03475 }
03476 }
03477 for(int index_vxxinv=rows_VxxInv[i];index_vxxinv<rows_VxxInv[i+1];
03478 index_vxxinv++) {
03479 if(cols_VxxInv[index_vxxinv]==i) {
03480 einv_ii=data_VxxInv[index_vxxinv];
03481 break;
03482 }
03483 }
03484 Double_t rho=1.0;
03485 if((einv_ii>0.0)&&(e_ii>0.0)) rho=1.-1./(einv_ii*e_ii);
03486 if(rho>=0.0) rho=TMath::Sqrt(rho);
03487 else rho=-TMath::Sqrt(-rho);
03488 if(rho>rhoMax) {
03489 rhoMax = rho;
03490 }
03491 rhoi->SetBinContent(destI,rho);
03492 }
03493 return rhoMax;
03494 }
03495 }
03496
03497 Double_t TUnfoldV17::GetRhoIFromMatrix(TH1 *rhoi,const TMatrixDSparse *eOrig,
03498 const Int_t *binMap,TH2 *invEmat) const
03499 {
03500
03501
03502
03503
03504
03505
03506
03507
03508
03509
03510
03511 Double_t rhoMax=0.;
03512
03513
03514
03515
03516 Int_t binMapSize = fHistToX.GetSize();
03517
03518
03519
03520
03521 std::map<int,int> histToLocalBin;
03522 Int_t nBin=0;
03523 for(Int_t i=0;i<binMapSize;i++) {
03524 Int_t mapped_i=binMap ? binMap[i] : i;
03525 if(mapped_i>=0) {
03526 if(histToLocalBin.find(mapped_i)==histToLocalBin.end()) {
03527 histToLocalBin[mapped_i]=nBin;
03528 nBin++;
03529 }
03530 }
03531 }
03532
03533 if(nBin>0) {
03534
03535
03536 Int_t *localBinToHist=new Int_t[nBin];
03537 for(std::map<int,int>::const_iterator i=histToLocalBin.begin();
03538 i!=histToLocalBin.end();i++) {
03539 localBinToHist[(*i).second]=(*i).first;
03540 }
03541
03542 const Int_t *rows_eOrig=eOrig->GetRowIndexArray();
03543 const Int_t *cols_eOrig=eOrig->GetColIndexArray();
03544 const Double_t *data_eOrig=eOrig->GetMatrixArray();
03545
03546
03547
03548
03549
03550 TMatrixD eRemap(nBin,nBin);
03551
03552 for(Int_t i=0;i<GetNx();i++) {
03553
03554 Int_t origI=fXToHist[i];
03555
03556 Int_t destI=binMap ? binMap[origI] : origI;
03557 if(destI<0) continue;
03558 Int_t ematBinI=histToLocalBin[destI];
03559 for(int index_eOrig=rows_eOrig[i];index_eOrig<rows_eOrig[i+1];
03560 index_eOrig++) {
03561
03562 Int_t j=cols_eOrig[index_eOrig];
03563
03564 Int_t origJ=fXToHist[j];
03565
03566 Int_t destJ=binMap ? binMap[origJ] : origJ;
03567 if(destJ<0) continue;
03568 Int_t ematBinJ=histToLocalBin[destJ];
03569 eRemap(ematBinI,ematBinJ) += data_eOrig[index_eOrig];
03570 }
03571 }
03572
03573 TMatrixDSparse eSparse(eRemap);
03574 Int_t rank=0;
03575 TMatrixDSparse *einvSparse=InvertMSparseSymmPos(&eSparse,&rank);
03576 if(rank!=einvSparse->GetNrows()) {
03577 Warning("GetRhoIFormMatrix","Covariance matrix has rank %d expect %d",
03578 rank,einvSparse->GetNrows());
03579 }
03580
03581 const Int_t *rows_eInv=einvSparse->GetRowIndexArray();
03582 const Int_t *cols_eInv=einvSparse->GetColIndexArray();
03583 const Double_t *data_eInv=einvSparse->GetMatrixArray();
03584
03585 for(Int_t i=0;i<einvSparse->GetNrows();i++) {
03586 Double_t e_ii=eRemap(i,i);
03587 Double_t einv_ii=0.0;
03588 for(Int_t index_einv=rows_eInv[i];index_einv<rows_eInv[i+1];
03589 index_einv++) {
03590 Int_t j=cols_eInv[index_einv];
03591 if(j==i) {
03592 einv_ii=data_eInv[index_einv];
03593 }
03594 if(invEmat) {
03595 invEmat->SetBinContent(localBinToHist[i],localBinToHist[j],
03596 data_eInv[index_einv]);
03597 }
03598 }
03599 Double_t rho=1.0;
03600 if((einv_ii>0.0)&&(e_ii>0.0)) rho=1.-1./(einv_ii*e_ii);
03601 if(rho>=0.0) rho=TMath::Sqrt(rho);
03602 else rho=-TMath::Sqrt(-rho);
03603 if(rho>rhoMax) {
03604 rhoMax = rho;
03605 }
03606
03607 rhoi->SetBinContent(localBinToHist[i],rho);
03608 }
03609
03610 DeleteMatrix(&einvSparse);
03611 delete [] localBinToHist;
03612 }
03613 return rhoMax;
03614 }
03615
03616
03617
03618
03619
03620
03621
03622
03623
03624 void TUnfoldV17::ClearHistogram(TH1 *h,Double_t x) const
03625 {
03626 Int_t nxyz[3];
03627 nxyz[0]=h->GetNbinsX()+1;
03628 nxyz[1]=h->GetNbinsY()+1;
03629 nxyz[2]=h->GetNbinsZ()+1;
03630 for(int i=h->GetDimension();i<3;i++) nxyz[i]=0;
03631 Int_t ixyz[3];
03632 for(int i=0;i<3;i++) ixyz[i]=0;
03633 while((ixyz[0]<=nxyz[0])&&
03634 (ixyz[1]<=nxyz[1])&&
03635 (ixyz[2]<=nxyz[2])){
03636 Int_t ibin=h->GetBin(ixyz[0],ixyz[1],ixyz[2]);
03637 h->SetBinContent(ibin,x);
03638 h->SetBinError(ibin,0.0);
03639 for(Int_t i=0;i<3;i++) {
03640 ixyz[i] += 1;
03641 if(ixyz[i]<=nxyz[i]) break;
03642 if(i<2) ixyz[i]=0;
03643 }
03644 }
03645 }
03646
03647 void TUnfoldV17::SetEpsMatrix(Double_t eps) {
03648
03649 if((eps>0.0)&&(eps<1.0)) fEpsMatrix=eps;
03650 }
03651
03652
03653
03654
03655
03656
03657
03658
03659
03660
03661 const char *TUnfoldV17::GetTUnfoldVersion(void)
03662 {
03663 return TUnfold_VERSION;
03664 }
03665