230 for(
Int_t i=0;i<2;i++) {
334 if(rank !=
GetNx()) {
335 Warning(
"DoUnfold",
"rank of matrix E %d expect %d",rank,
GetNx());
359 epsilonNosparse(A_cols[i],0) += A_data[i];
371 Fatal(
"Unfold",
"epsilon#Eepsilon has dimension %d != 1",
378 y_minus_epsx += (*fY)(iy,0);
382 y_minus_epsx -= epsilonNosparse(ix,0) * (*fX)(ix,0);
385 lambda_half=y_minus_epsx*one_over_epsEeps;
390 if(EEpsilon_rows[ix]<EEpsilon_rows[ix+1]) {
391 (*fX)(ix,0) += lambda_half * EEpsilon_data[EEpsilon_rows[ix]];
410 data[i]=one_over_epsEeps;
420 AddMSparse(temp, -one_over_epsEeps, epsilonB);
457 if(VyyinvDy_rows[iy]<VyyinvDy_rows[iy+1]) {
458 fChi2A += VyyinvDy_data[VyyinvDy_rows[iy]]*dy(iy,0);
467 if(LsquaredDx_rows[ix]<LsquaredDx_rows[ix+1]) {
468 fLXsquared += LsquaredDx_data[LsquaredDx_rows[ix]]*dx(ix,0);
483 "epsilon#fDXDtauSquared has dimension %d != 1",
506 (Eepsilon,Eepsilon,0);
531 if(rank !=
GetNx()) {
532 Warning(
"DoUnfold",
"rank of output covariance is %d expect %d",
541 for(
int ik=VxxInv_rows[ix];ik<VxxInv_rows[ix+1];ik++) {
542 if(ix==VxxInv_cols[ik]) {
543 VxxInvDiag(ix)=VxxInv_data[ik];
557 for(
int ik=Vxx_rows[ix];ik<Vxx_rows[ix+1];ik++) {
558 if(ix==Vxx_cols[ik]) {
560 1. - 1. / (VxxInvDiag(ix) * Vxx_data[ik]);
561 if (rho_squared > rho_squared_max)
562 rho_squared_max = rho_squared;
563 if(rho_squared>0.0) {
572 fRhoAvg = (n_rho>0) ? (rho_sum/n_rho) : -1.0;
600 A->SetMatrixArray(nel,row,col,data);
620 if(
a->GetNcols()!=
b->GetNrows()) {
621 Fatal(
"MultiplyMSparseMSparse",
622 "inconsistent matrix col/ matrix row %d !=%d",
623 a->GetNcols(),
b->GetNrows());
627 const Int_t *a_rows=
a->GetRowIndexArray();
628 const Int_t *a_cols=
a->GetColIndexArray();
629 const Double_t *a_data=
a->GetMatrixArray();
630 const Int_t *b_rows=
b->GetRowIndexArray();
631 const Int_t *b_cols=
b->GetColIndexArray();
632 const Double_t *b_data=
b->GetMatrixArray();
635 for (
Int_t irow = 0; irow <
a->GetNrows(); irow++) {
636 if(a_rows[irow+1]>a_rows[irow]) nMax +=
b->GetNcols();
638 if((nMax>0)&&(a_cols)&&(b_cols)) {
644 for (
Int_t irow = 0; irow <
a->GetNrows(); irow++) {
645 if(a_rows[irow+1]<=a_rows[irow])
continue;
647 for(
Int_t icol=0;icol<
b->GetNcols();icol++) {
651 for(
Int_t ia=a_rows[irow];ia<a_rows[irow+1];ia++) {
654 for(
Int_t ib=b_rows[k];ib<b_rows[k+1];ib++) {
655 row_data[b_cols[ib]] += a_data[ia]*b_data[ib];
659 for(
Int_t icol=0;icol<
b->GetNcols();icol++) {
660 if(row_data[icol] != 0.0) {
663 r_data[
n]=row_data[icol];
669 r->SetMatrixArray(
n,r_rows,r_cols,r_data);
695 if(
a->GetNrows() !=
b->GetNrows()) {
696 Fatal(
"MultiplyMSparseTranspMSparse",
697 "inconsistent matrix row numbers %d!=%d",
698 a->GetNrows(),
b->GetNrows());
702 const Int_t *a_rows=
a->GetRowIndexArray();
703 const Int_t *a_cols=
a->GetColIndexArray();
704 const Double_t *a_data=
a->GetMatrixArray();
705 const Int_t *b_rows=
b->GetRowIndexArray();
706 const Int_t *b_cols=
b->GetColIndexArray();
707 const Double_t *b_data=
b->GetMatrixArray();
711 typedef std::map<Int_t,Double_t> MMatrixRow_t;
712 typedef std::map<Int_t, MMatrixRow_t > MMatrix_t;
715 for(
Int_t iRowAB=0;iRowAB<
a->GetNrows();iRowAB++) {
716 for(
Int_t ia=a_rows[iRowAB];ia<a_rows[iRowAB+1];ia++) {
717 for(
Int_t ib=b_rows[iRowAB];ib<b_rows[iRowAB+1];ib++) {
719 MMatrixRow_t &row=matrix[a_cols[ia]];
720 MMatrixRow_t::iterator icol=row.find(b_cols[ib]);
721 if(icol!=row.end()) {
723 (*icol).second += a_data[ia]*b_data[ib];
726 row[b_cols[ib]] = a_data[ia]*b_data[ib];
733 for (MMatrix_t::const_iterator irow = matrix.begin(); irow != matrix.end(); ++irow) {
734 n += (*irow).second.size();
742 for (MMatrix_t::const_iterator irow = matrix.begin(); irow != matrix.end(); ++irow) {
743 for (MMatrixRow_t::const_iterator icol = (*irow).second.begin(); icol != (*irow).second.end(); ++icol) {
744 r_rows[
n]=(*irow).first;
745 r_cols[
n]=(*icol).first;
746 r_data[
n]=(*icol).second;
752 r->SetMatrixArray(
n,r_rows,r_cols,r_data);
776 if(
a->GetNcols()!=
b->GetNrows()) {
777 Fatal(
"MultiplyMSparseM",
"inconsistent matrix col /matrix row %d!=%d",
778 a->GetNcols(),
b->GetNrows());
782 const Int_t *a_rows=
a->GetRowIndexArray();
783 const Int_t *a_cols=
a->GetColIndexArray();
784 const Double_t *a_data=
a->GetMatrixArray();
787 for (
Int_t irow = 0; irow <
a->GetNrows(); irow++) {
788 if(a_rows[irow+1]-a_rows[irow]>0) nMax +=
b->GetNcols();
797 for (
Int_t irow = 0; irow <
a->GetNrows(); irow++) {
798 if(a_rows[irow+1]-a_rows[irow]<=0)
continue;
799 for(
Int_t icol=0;icol<
b->GetNcols();icol++) {
803 for(
Int_t i=a_rows[irow];i<a_rows[irow+1];i++) {
805 r_data[
n] += a_data[i]*(*b)(j,icol);
807 if(r_data[
n]!=0.0)
n++;
811 r->SetMatrixArray(
n,r_rows,r_cols,r_data);
837 (
v && ((m1->
GetNcols()!=
v->GetNrows())||(
v->GetNcols()!=1)))) {
839 Fatal(
"MultiplyMSparseMSparseTranspVector",
840 "matrix cols/vector rows %d!=%d!=%d or vector rows %d!=1\n",
841 m1->
GetNcols(),
m2->GetNcols(),
v->GetNrows(),
v->GetNcols());
843 Fatal(
"MultiplyMSparseMSparseTranspVector",
844 "matrix cols %d!=%d\n",m1->
GetNcols(),
m2->GetNcols());
852 if(rows_m1[i]<rows_m1[i+1]) num_m1++;
854 const Int_t *rows_m2=
m2->GetRowIndexArray();
855 const Int_t *cols_m2=
m2->GetColIndexArray();
858 for(
Int_t j=0;j<
m2->GetNrows();j++) {
859 if(rows_m2[j]<rows_m2[j+1]) num_m2++;
862 const Int_t *v_rows=0;
868 Int_t num_r=num_m1*num_m2+1;
874 for(
Int_t j=0;j<
m2->GetNrows();j++) {
876 Int_t index_m1=rows_m1[i];
877 Int_t index_m2=rows_m2[j];
878 while((index_m1<rows_m1[i+1])&&(index_m2<rows_m2[j+1])) {
879 Int_t k1=cols_m1[index_m1];
880 Int_t k2=cols_m2[index_m2];
887 Int_t v_index=v_rows[k1];
888 if(v_index<v_rows[k1+1]) {
889 data_r[num_r] += data_m1[index_m1] * data_m2[index_m2]
893 data_r[num_r] += data_m1[index_m1] * data_m2[index_m2]
896 data_r[num_r] += data_m1[index_m1] * data_m2[index_m2];
902 if(data_r[num_r] !=0.0) {
910 num_r,row_r,col_r,data_r);
933 const Int_t *dest_rows=
dest->GetRowIndexArray();
934 const Int_t *dest_cols=
dest->GetColIndexArray();
942 Fatal(
"AddMSparse",
"inconsistent matrix rows %d!=%d OR cols %d!=%d",
950 for(
Int_t row=0;row<
dest->GetNrows();row++) {
951 Int_t i_dest=dest_rows[row];
952 Int_t i_src=src_rows[row];
953 while((i_dest<dest_rows[row+1])||(i_src<src_rows[row+1])) {
954 Int_t col_dest=(i_dest<dest_rows[row+1]) ?
955 dest_cols[i_dest] :
dest->GetNcols();
956 Int_t col_src =(i_src <src_rows[row+1] ) ?
959 if(col_dest<col_src) {
960 result_cols[
n]=col_dest;
961 result_data[
n]=dest_data[i_dest++];
962 }
else if(col_dest>col_src) {
963 result_cols[
n]=col_src;
964 result_data[
n]=
f*src_data[i_src++];
966 result_cols[
n]=col_dest;
967 result_data[
n]=dest_data[i_dest++]+
f*src_data[i_src++];
969 if(result_data[
n] !=0.0) {
971 Fatal(
"AddMSparse",
"Nan detected %d %d %d",
984 dest->SetMatrixArray(
n,result_rows,result_cols,result_data);
985 delete[] result_data;
986 delete[] result_rows;
987 delete[] result_cols;
1011 if(
A->GetNcols()!=
A->GetNrows()) {
1012 Fatal(
"InvertMSparseSymmPos",
"inconsistent matrix row/col %d!=%d",
1013 A->GetNcols(),
A->GetNrows());
1017 const Int_t *a_rows=
A->GetRowIndexArray();
1018 const Int_t *a_cols=
A->GetColIndexArray();
1019 const Double_t *a_data=
A->GetMatrixArray();
1024 Int_t iBlock=
A->GetNrows();
1028 for(
Int_t iA=0;iA<
A->GetNrows();iA++) {
1029 for(
Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
1030 Int_t jA=a_cols[indexA];
1032 if(!(a_data[indexA]>=0.0)) nError++;
1033 aII(iA)=a_data[indexA];
1038 Fatal(
"InvertMSparseSymmPos",
1039 "Matrix has %d negative elements on the diagonal", nError);
1054 for(
Int_t j=0;j<
A->GetNrows();j++)
swap[j]=j;
1055 for(
Int_t iStart=0;iStart<iBlock;iStart++) {
1057 for(
Int_t i=iStart;i<iBlock;i++) {
1059 Int_t n=
A->GetNrows()-(a_rows[iA+1]-a_rows[iA]);
1067 for(
Int_t i=0;i<
A->GetNrows();i++) isZero[i]=
kTRUE;
1069 for(
Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
1070 Int_t jA=a_cols[indexA];
1079 for(
Int_t i=iStart+1;i<iBlock;) {
1080 if(isZero[
swap[i]]) {
1097#ifdef CONDITION_BLOCK_PART
1098 Int_t nn=
A->GetNrows()-iBlock;
1099 for(
int inc=(nn+1)/2;inc;inc /= 2) {
1100 for(
int i=inc;i<nn;i++) {
1101 int itmp=
swap[i+iBlock];
1103 for(j=i;(j>=inc)&&(aII(
swap[j-inc+iBlock])<aII(itmp));j -=inc) {
1106 swap[j+iBlock]=itmp;
1112 for(
Int_t i=0;i<
A->GetNrows();i++) {
1113 swapBack[
swap[i]]=i;
1116 for(
Int_t i=0;i<
A->GetNrows();i++) {
1117 std::cout<<
" "<<i<<
" "<<
swap[i]<<
" "<<swapBack[i]<<
"\n";
1119 std::cout<<
"after sorting\n";
1120 for(
Int_t i=0;i<
A->GetNrows();i++) {
1121 if(i==iDiagonal) std::cout<<
"iDiagonal="<<i<<
"\n";
1122 if(i==iBlock) std::cout<<
"iBlock="<<i<<
"\n";
1123 std::cout<<
" "<<
swap[i]<<
" "<<aII(
swap[i])<<
"\n";
1135 for(
Int_t i=0;i<
A->GetNrows();i++) {
1137 for(
Int_t indexA=a_rows[row];indexA<a_rows[row+1];indexA++) {
1138 Int_t j=swapBack[a_cols[indexA]];
1140 else if((i<iDiagonal)||(j<iDiagonal)) nError1++;
1141 else if((i<iBlock)&&(j<iBlock)) nError2++;
1142 else if((i<iBlock)||(j<iBlock)) nBlock++;
1146 if(nError1+nError2>0) {
1147 Fatal(
"InvertMSparseSymmPos",
"sparse matrix analysis failed %d %d %d %d %d",
1148 iDiagonal,iBlock,
A->GetNrows(),nError1,nError2);
1153 Info(
"InvertMSparseSymmPos",
"iDiagonal=%d iBlock=%d nRow=%d",
1154 iDiagonal,iBlock,
A->GetNrows());
1213 for(
Int_t i=0;i<iDiagonal;++i) {
1218 rEl_data[rNumEl]=1./aII(iA);
1223 if((!rankPtr)&&(rankD1!=iDiagonal)) {
1224 Fatal(
"InvertMSparseSymmPos",
1225 "diagonal part 1 has rank %d != %d, matrix can not be inverted",
1233 Int_t nD2=iBlock-iDiagonal;
1235 if((rNumEl>=0)&&(nD2>0)) {
1240 for(
Int_t i=0;i<nD2;++i) {
1243 D2inv_col[D2invNumEl]=i;
1244 D2inv_row[D2invNumEl]=i;
1245 D2inv_data[D2invNumEl]=1./aII(iA);
1249 if(D2invNumEl==nD2) {
1252 }
else if(!rankPtr) {
1253 Fatal(
"InvertMSparseSymmPos",
1254 "diagonal part 2 has rank %d != %d, matrix can not be inverted",
1258 delete [] D2inv_data;
1259 delete [] D2inv_col;
1260 delete [] D2inv_row;
1266 Int_t nF=
A->GetNrows()-iBlock;
1270 if((rNumEl>=0)&&(nF>0)&&((nD2==0)||D2inv)) {
1279 Int_t nBmax=nF*(nD2+1);
1285 for(
Int_t i=0;i<nF;i++) {
1287 for(
Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
1288 Int_t jA=a_cols[indexA];
1289 Int_t j0=swapBack[jA];
1294 F_data[FnumEl]=a_data[indexA];
1296 }
else if(j0>=iDiagonal) {
1297 Int_t j=j0-iDiagonal;
1300 B_data[BnumEl]=a_data[indexA];
1307#ifndef FORCE_EIGENVALUE_DECOMPOSITION
1327 for(
Int_t i=0;i<mbd2_nMax;i++) {
1328 mbd2_data[i] = - mbd2_data[i];
1332 if(minusBD2inv &&
F) {
1341 const Int_t *f_rows=
F->GetRowIndexArray();
1342 const Int_t *f_cols=
F->GetColIndexArray();
1343 const Double_t *f_data=
F->GetMatrixArray();
1347 for(
Int_t i=0;i<nF;i++) {
1348 for(
Int_t indexF=f_rows[i];indexF<f_rows[i+1];indexF++) {
1349 if(f_cols[indexF]>=i)
c(f_cols[indexF],i)=f_data[indexF];
1353 for(
Int_t j=0;j<i;j++) {
1364 for(
Int_t j=i+1;j<nF;j++) {
1366 for(
Int_t k=0;k<i;k++) {
1367 c_ji -=
c(i,k)*
c(j,k);
1376 for(
Int_t i=1;i<nF;i++) {
1381 std::cout<<
"dmin,dmax: "<<dmin<<
" "<<dmax<<
"\n";
1383 if(dmin/dmax<epsilonF2*nF) nErrorF=2;
1389 for(
Int_t i=0;i<nF;i++) {
1390 cinv(i,i)=1./
c(i,i);
1392 for(
Int_t i=0;i<nF;i++) {
1393 for(
Int_t j=i+1;j<nF;j++) {
1395 for(
Int_t k=i+1;k<j;k++) {
1396 tmp -= cinv(k,i)*
c(j,k);
1398 cinv(j,i)=tmp*cinv(j,j);
1403 (&cInvSparse,&cInvSparse);
1414 Int_t rankD2Block=0;
1416 if( ((nD2==0)||D2inv) && ((nF==0)||Finv) ) {
1426 if(Finv && minusBD2inv) {
1431 if(
G && minusBD2inv) {
1438 E=minusBD2invTransG;
1443 const Int_t *e_rows=
E->GetRowIndexArray();
1444 const Int_t *e_cols=
E->GetColIndexArray();
1445 const Double_t *e_data=
E->GetMatrixArray();
1446 for(
Int_t iE=0;iE<
E->GetNrows();iE++) {
1448 for(
Int_t indexE=e_rows[iE];indexE<e_rows[iE+1];++indexE) {
1449 Int_t jE=e_cols[indexE];
1453 rEl_data[rNumEl]=e_data[indexE];
1461 const Int_t *g_rows=
G->GetRowIndexArray();
1462 const Int_t *g_cols=
G->GetColIndexArray();
1463 const Double_t *g_data=
G->GetMatrixArray();
1464 for(
Int_t iG=0;iG<
G->GetNrows();iG++) {
1466 for(
Int_t indexG=g_rows[iG];indexG<g_rows[iG+1];++indexG) {
1467 Int_t jG=g_cols[indexG];
1472 rEl_data[rNumEl]=g_data[indexG];
1477 rEl_data[rNumEl]=g_data[indexG];
1490 for(
Int_t indexF=finv_rows[iF];indexF<finv_rows[iF+1];++indexF) {
1491 Int_t jF=finv_cols[indexF];
1495 rEl_data[rNumEl]=finv_data[indexF];
1501 }
else if(!rankPtr) {
1503 Fatal(
"InvertMSparseSymmPos",
1504 "non-trivial part has rank < %d, matrix can not be inverted",
1511 Info(
"InvertMSparseSymmPos",
1512 "cholesky-decomposition failed, try eigenvalue analysis");
1514 std::cout<<
"nEV="<<nEV<<
" iDiagonal="<<iDiagonal<<
"\n";
1517 for(
Int_t i=0;i<nEV;i++) {
1519 for(
Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
1520 Int_t jA=a_cols[indexA];
1521 Int_t j0=swapBack[jA];
1523 Int_t j=j0-iDiagonal;
1525 if((i<0)||(j<0)||(i>=nEV)||(j>=nEV)) {
1526 std::cout<<
" error "<<nEV<<
" "<<i<<
" "<<j<<
"\n";
1530 Fatal(
"InvertMSparseSymmPos",
1531 "non-finite number detected element %d %d\n",
1534 EV(i,j)=a_data[indexA];
1541 std::cout<<
"Eigenvalues\n";
1551 if(condition<-epsilonEV2) {
1556 Error(
"InvertMSparseSymmPos",
1557 "Largest Eigenvalue is negative %f",
1560 Error(
"InvertMSparseSymmPos",
1561 "Some Eigenvalues are negative (EV%d/EV0=%g epsilon=%g)",
1571 for(
Int_t i=0;i<nEV;i++) {
1574 inverseEV(i,0)=1./
x;
1588 for(
Int_t indexVDVt=vdvt_rows[iVDVt];
1589 indexVDVt<vdvt_rows[iVDVt+1];++indexVDVt) {
1590 Int_t jVDVt=vdvt_cols[indexVDVt];
1594 rEl_data[rNumEl]=vdvt_data[indexVDVt];
1602 *rankPtr=rankD1+rankD2Block;
1617 rEl_row,rEl_col,rEl_data) : 0;
1640 for(
Int_t i=0;i<
a.GetNrows();i++) {
1641 for(
Int_t j=0;j<
a.GetNcols();j++) {
1645 std::cout<<
"Ar is not symmetric Ar("<<i<<
","<<j<<
")="<<ar(i,j)
1646 <<
" Ar("<<j<<
","<<i<<
")="<<ar(j,i)<<
"\n";
1651 std::cout<<
"ArA is not equal A ArA("<<i<<
","<<j<<
")="<<ara(i,j)
1652 <<
" A("<<i<<
","<<j<<
")="<<
a(i,j)<<
"\n";
1657 std::cout<<
"rAr is not equal r rAr("<<i<<
","<<j<<
")="<<rar(i,j)
1658 <<
" r("<<i<<
","<<j<<
")="<<
R(i,j)<<
"\n";
1662 if(rankPtr) std::cout<<
"rank="<<*rankPtr<<
"\n";
1664 std::cout<<
"Matrix is not positive\n";
1754 for (
Int_t ix = 0; ix < nx0; ix++) {
1759 for (
Int_t iy = 0; iy < ny; iy++) {
1797 Int_t underflowBin=0,overflowBin=0;
1798 for (
Int_t ix = 0; ix < nx; ix++) {
1800 if(ibinx<1) underflowBin=1;
1802 if(ibinx>hist_A->
GetNbinsX()) overflowBin=1;
1804 if(ibinx>hist_A->
GetNbinsY()) overflowBin=1;
1809 Int_t ixfirst=-1,ixlast=-1;
1811 int nDisconnected=0;
1812 for (
Int_t ix = 0; ix < nx0; ix++) {
1823 if(((
fHistToX[ix]>=0)&&(ixlast>=0))||
1824 (nprint==nskipped)) {
1825 if(ixlast>ixfirst) {
1832 if(nprint==nskipped) {
1836 if(nskipped==(2-underflowBin-overflowBin)) {
1837 Info(
"TUnfold",
"underflow and overflow bin "
1838 "do not depend on the input data");
1840 Warning(
"TUnfold",
"%d output bins "
1841 "do not depend on the input data %s",nDisconnected,
1842 static_cast<const char *
>(binlist));
1853 for (
Int_t iy = 0; iy < ny; iy++) {
1854 for (
Int_t ix = 0; ix < nx; ix++) {
1870 if(underflowBin && overflowBin) {
1871 Info(
"TUnfold",
"%d input bins and %d output bins (includes 2 underflow/overflow bins)",ny,nx);
1872 }
else if(underflowBin) {
1873 Info(
"TUnfold",
"%d input bins and %d output bins (includes 1 underflow bin)",ny,nx);
1874 }
else if(overflowBin) {
1875 Info(
"TUnfold",
"%d input bins and %d output bins (includes 1 overflow bin)",ny,nx);
1877 Info(
"TUnfold",
"%d input bins and %d output bins",ny,nx);
1881 Error(
"TUnfold",
"too few (ny=%d) input bins for nx=%d output bins",ny,nx);
1883 Warning(
"TUnfold",
"too few (ny=%d) input bins for nx=%d output bins",ny,nx);
1895 "%d regularisation conditions have been skipped",nError);
1898 "One regularisation condition has been skipped");
1988 for(
Int_t k=l0_rows[row];k<l0_rows[row+1];k++) {
1990 l_col[nF]=l0_cols[k];
1991 l_data[nF]=l0_data[k];
2003 for(
Int_t i=0;i<nEle;i++) {
2011 l_data[nF]=rowData[i];
2135 (left_bin,-scale_left,
2136 center_bin,scale_left+scale_right,
2137 right_bin,-scale_right)
2184 Error(
"RegularizeBins",
"regmode = %d is not valid",regmode);
2186 for (
Int_t i = nSkip; i < nbin; i++) {
2223 int step2,
int nbin2,
ERegMode regmode)
2236 for (
Int_t i1 = 0; i1 < nbin1; i1++) {
2237 nError +=
RegularizeBins(start_bin + step1 * i1, step2, nbin2, regmode);
2239 for (
Int_t i2 = 0; i2 < nbin2; i2++) {
2240 nError +=
RegularizeBins(start_bin + step2 * i2, step1, nbin1, regmode);
2302 const TH2 *hist_vyy_inv)
2324 Int_t nVarianceZero=0;
2325 Int_t nVarianceForced=0;
2335 if(oneOverZeroError>0.0) {
2336 dy2 = 1./ ( oneOverZeroError*oneOverZeroError);
2348 rowVyyN[nVyyN] = iy;
2349 colVyyN[nVyyN] = iy;
2350 rowVyy1[nVyy1] = iy;
2352 dataVyyDiag[iy] = dy2;
2354 dataVyyN[nVyyN++] = dy2;
2355 dataVyy1[nVyy1++] = dy2;
2360 Int_t nOffDiagNonzero=0;
2363 if(dataVyyDiag[iy]<=0.0) {
2373 if(iy==jy)
continue;
2375 if(dataVyyDiag[jy]<=0.0)
continue;
2377 rowVyyN[nVyyN] = iy;
2378 colVyyN[nVyyN] = jy;
2380 if(dataVyyN[nVyyN] == 0.0)
continue;
2386 "inverse of input covariance is taken from user input");
2393 rowVyyInv[nVyyInv] = iy;
2394 colVyyInv[nVyyInv] = jy;
2395 dataVyyInv[nVyyInv]= hist_vyy_inv->
GetBinContent(iy+1,jy+1);
2396 if(dataVyyInv[nVyyInv] == 0.0)
continue;
2401 (
GetNy(),
GetNy(),nVyyInv,rowVyyInv,colVyyInv,dataVyyInv);
2402 delete [] rowVyyInv;
2403 delete [] colVyyInv;
2404 delete [] dataVyyInv;
2406 if(nOffDiagNonzero) {
2408 "input covariance has elements C(X,Y)!=0 where V(X)==0");
2414 (
GetNy(),
GetNy(),nVyyN,rowVyyN,colVyyN,dataVyyN);
2421 (
GetNy(),1,nVyy1,rowVyy1,colVyy1, dataVyy1);
2444 if(nVarianceForced) {
2445 if(nVarianceForced>1) {
2446 Warning(
"SetInput",
"%d/%d input bins have zero error,"
2447 " 1/error set to %lf.",
2448 nVarianceForced,
GetNy(),oneOverZeroError);
2450 Warning(
"SetInput",
"One input bin has zero error,"
2451 " 1/error set to %lf.",oneOverZeroError);
2455 if(nVarianceZero>1) {
2456 Warning(
"SetInput",
"%d/%d input bins have zero error,"
2457 " and are ignored.",nVarianceZero,
GetNy());
2459 Warning(
"SetInput",
"One input bin has zero error,"
2460 " and is ignored.");
2466 if(oneOverZeroError<=0.0) {
2472 TString binlist(
"no data to constrain output bin ");
2483 Warning(
"SetInput",
"%s",(
char const *)binlist);
2488 Error(
"SetInput",
"%d/%d output bins are not constrained by any data.",
2491 Error(
"SetInput",
"One output bins is not constrained by any data.");
2496 delete[] dataVyyDiag;
2498 return nVarianceForced+nVarianceZero+10000*nError2;
2553 typedef std::map<Double_t,std::pair<Double_t,Double_t> > XYtau_t;
2578 if((tauMin<=0)||(tauMax<=0.0)||(tauMin>=tauMax)) {
2588 Error(
"ScanLcurve",
"too few input bins, NDF<=0 %d",
GetNdf());
2593 Info(
"ScanLcurve",
"logtau=-Infinity X=%lf Y=%lf",x0,y0);
2595 Fatal(
"ScanLcurve",
"problem (too few input bins?) X=%f",x0);
2598 Fatal(
"ScanLcurve",
"problem (missing regularisation?) Y=%f",y0);
2607 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2611 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2614 if((*curve.begin()).second.first<x0) {
2626 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2630 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2633 while(((
int)curve.size()<(nPoint-1)/2)&&
2634 ((*curve.begin()).second.first<x0));
2641 while(((
int)curve.size()<nPoint-1)&&
2642 (((*curve.begin()).second.first-x0>0.00432)||
2643 ((*curve.begin()).second.second-y0>0.00432)||
2644 (curve.size()<2))) {
2648 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2652 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2663 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2666 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2673 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2676 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2686 while(
int(curve.size())<nPoint-1) {
2689 XYtau_t::const_iterator i0,i1;
2694 for (++i1; i1 != curve.end(); ++i1) {
2695 const std::pair<Double_t, Double_t> &xy0 = (*i0).second;
2696 const std::pair<Double_t, Double_t> &xy1 = (*i1).second;
2697 Double_t dx = xy1.first - xy0.first;
2698 Double_t dy = xy1.second - xy0.second;
2702 logTau = 0.5 * ((*i0).first + (*i1).first);
2708 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2720 XYtau_t::const_iterator i0,i1;
2725 if( ((
int)curve.size())<nPoint) {
2734 for (XYtau_t::const_iterator i = curve.begin(); i != curve.end(); ++i) {
2735 lXi[
n] = (*i).second.first;
2736 lYi[
n] = (*i).second.second;
2737 lTi[
n] = (*i).first;
2744 for(
Int_t i=0;i<
n-1;i++) {
2747 Double_t tauBar=0.5*(lTi[i]+lTi[i+1]);
2748 Double_t dTau=0.5*(lTi[i+1]-lTi[i]);
2749 Double_t dx1=bi+dTau*(2.*ci+3.*di*dTau);
2752 Double_t dy1=bi+dTau*(2.*ci+3.*di*dTau);
2755 cCi[i]=(dy2*dx1-dy1*dx2)/
TMath::Power(dx1*dx1+dy1*dy1,1.5);
2775 for(
Int_t i=iskip;i<
n-2-iskip;i++) {
2798 xx = m_p_half + discr;
2800 xx = m_p_half - discr;
2804 if((xx>0.0)&&(xx<dx)) {
2818 if((xx>0.0)&&(xx<dx)) {
2832 if(logTauCurvature) {
2833 *logTauCurvature=splineC;
2842 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2845 Info(
"ScanLcurve",
"Result logtau=%lf X=%lf Y=%lf",
2855 Int_t bestChoice=-1;
2856 if(curve.size()>0) {
2861 for (XYtau_t::const_iterator i = curve.begin(); i != curve.end(); ++i) {
2862 if (logTauFin == (*i).first) {
2865 x[
n] = (*i).second.first;
2866 y[
n] = (*i).second.second;
2867 logT[
n] = (*i).first;
2872 (*lCurve)->SetTitle(
"L curve");
2874 if(logTauX) (*logTauX)=
new TSpline3(
"log(chi**2)%log(tau)",logT,
x,
n);
2875 if(logTauY) (*logTauY)=
new TSpline3(
"log(reg.cond)%log(tau)",logT,
y,
n);
2962 Int_t destI = binMap ? binMap[i+1] : i+1;
2963 if(destI<0)
continue;
2967 Int_t index_a=rows_A[i];
2968 Int_t index_av=rows_AVxx[i];
2969 while((index_a<rows_A[i+1])&&(index_av<rows_AVxx[i+1])) {
2970 Int_t j_a=cols_A[index_a];
2971 Int_t j_av=cols_AVxx[index_av];
2974 }
else if(j_a>j_av) {
2977 e2 += data_AVxx[index_av] * data_A[index_a];
3005 for(
Int_t indexA=rows_A[iy];indexA<rows_A[iy+1];indexA++) {
3006 Int_t ix = cols_A[indexA];
3009 A->SetBinContent(ih, iy+1,data_A[indexA]);
3011 A->SetBinContent(iy+1, ih,data_A[indexA]);
3040 Int_t destI=binMap ? binMap[i+1] : i+1;
3041 if(destI<0)
continue;
3046 for(
int index=rows_Vyy[i];index<rows_Vyy[i+1];index++) {
3047 if(cols_Vyy[index]==i) {
3071 Warning(
"GetInputInverseEmatrix",
"input covariance matrix has rank %d expect %d",
3075 Error(
"GetInputInverseEmatrix",
"number of parameters %d > %d (rank of input covariance). Problem can not be solved",
GetNpar(),rank);
3076 }
else if(
fNdf==0) {
3077 Warning(
"GetInputInverseEmatrix",
"number of parameters %d = input rank %d. Problem is ill posed",
GetNpar(),rank);
3086 for(
int i=0;i<=out->
GetNbinsX()+1;i++) {
3087 for(
int j=0;j<=out->
GetNbinsY()+1;j++) {
3093 for(
int index=rows_VyyInv[i];index<rows_VyyInv[i+1];index++) {
3094 Int_t j=cols_VyyInv[index];
3124 for (
Int_t cindex = rows[i]; cindex < rows[i+1]; cindex++) {
3125 Int_t j=cols[cindex];
3160 for (
Int_t cindex = rows[row]; cindex < rows[row+1]; cindex++) {
3161 Int_t col=cols[cindex];
3274 std::map<Int_t,Double_t> e2;
3281 for(
Int_t i=0;i<binMapSize;i++) {
3282 Int_t destBinI=binMap ? binMap[i] : i;
3284 if((destBinI>=0)&&(srcBinI>=0)) {
3286 (destBinI, (*
fX)(srcBinI,0)+
output->GetBinContent(destBinI));
3292 Int_t index_vxx=rows_Vxx[srcBinI];
3294 while((j<binMapSize)&&(index_vxx<rows_Vxx[srcBinI+1])) {
3295 Int_t destBinJ=binMap ? binMap[j] : j;
3296 if(destBinI!=destBinJ) {
3305 if(cols_Vxx[index_vxx]<srcBinJ) {
3308 }
else if(cols_Vxx[index_vxx]>srcBinJ) {
3313 e2[destBinI] += data_Vxx[index_vxx];
3323 for (std::map<Int_t, Double_t>::const_iterator i = e2.begin(); i != e2.end(); ++i) {
3348 for(
Int_t i=0;i<nbin+2;i++) {
3349 for(
Int_t j=0;j<nbin+2;j++) {
3362 for(
Int_t i=0;i<binMapSize;i++) {
3363 Int_t destBinI=binMap ? binMap[i] : i;
3365 if((destBinI>=0)&&(destBinI<nbin+2)&&(srcBinI>=0)) {
3371 Int_t index_vxx=rows_emat[srcBinI];
3372 while((j<binMapSize)&&(index_vxx<rows_emat[srcBinI+1])) {
3373 Int_t destBinJ=binMap ? binMap[j] : j;
3375 if((destBinJ<0)||(destBinJ>=nbin+2)||(srcBinJ<0)) {
3379 if(cols_emat[index_vxx]<srcBinJ) {
3382 }
else if(cols_emat[index_vxx]>srcBinJ) {
3388 + data_emat[index_vxx];
3430 for(
Int_t i=0;i<nbin+2;i++) {
3433 for(
Int_t i=0;i<nbin+2;i++) {
3434 for(
Int_t j=0;j<nbin+2;j++) {
3435 if((
e[i]>0.0)&&(
e[j]>0.0)) {
3489 for(
int index_vxx=rows_Vxx[i];index_vxx<rows_Vxx[i+1];index_vxx++) {
3490 if(cols_Vxx[index_vxx]==i) {
3491 e_ii=data_Vxx[index_vxx];
3495 for(
int index_vxxinv=rows_VxxInv[i];index_vxxinv<rows_VxxInv[i+1];
3497 if(cols_VxxInv[index_vxxinv]==i) {
3498 einv_ii=data_VxxInv[index_vxxinv];
3503 if((einv_ii>0.0)&&(e_ii>0.0)) rho=1.-1./(einv_ii*e_ii);
3529 const Int_t *binMap,
TH2 *invEmat)
const
3541 std::map<int,int> histToLocalBin;
3543 for(
Int_t i=0;i<binMapSize;i++) {
3544 Int_t mapped_i=binMap ? binMap[i] : i;
3546 if(histToLocalBin.find(mapped_i)==histToLocalBin.end()) {
3547 histToLocalBin[mapped_i]=nBin;
3557 for (std::map<int, int>::const_iterator i = histToLocalBin.begin(); i != histToLocalBin.end(); ++i) {
3558 localBinToHist[(*i).second]=(*i).first;
3575 Int_t destI=binMap ? binMap[origI] : origI;
3576 if(destI<0)
continue;
3577 Int_t ematBinI=histToLocalBin[destI];
3578 for(
int index_eOrig=rows_eOrig[i];index_eOrig<rows_eOrig[i+1];
3581 Int_t j=cols_eOrig[index_eOrig];
3585 Int_t destJ=binMap ? binMap[origJ] : origJ;
3586 if(destJ<0)
continue;
3587 Int_t ematBinJ=histToLocalBin[destJ];
3588 eRemap(ematBinI,ematBinJ) += data_eOrig[index_eOrig];
3596 Warning(
"GetRhoIFormMatrix",
"Covariance matrix has rank %d expect %d",
3607 for(
Int_t index_einv=rows_eInv[i];index_einv<rows_eInv[i+1];
3609 Int_t j=cols_eInv[index_einv];
3611 einv_ii=data_eInv[index_einv];
3615 data_eInv[index_einv]);
3619 if((einv_ii>0.0)&&(e_ii>0.0)) rho=1.-1./(einv_ii*e_ii);
3630 delete [] localBinToHist;
3646 nxyz[0]=
h->GetNbinsX()+1;
3647 nxyz[1]=
h->GetNbinsY()+1;
3648 nxyz[2]=
h->GetNbinsZ()+1;
3649 for(
int i=
h->GetDimension();i<3;i++) nxyz[i]=0;
3651 for(
int i=0;i<3;i++) ixyz[i]=0;
3652 while((ixyz[0]<=nxyz[0])&&
3653 (ixyz[1]<=nxyz[1])&&
3654 (ixyz[2]<=nxyz[2])){
3655 Int_t ibin=
h->GetBin(ixyz[0],ixyz[1],ixyz[2]);
3656 h->SetBinContent(ibin,
x);
3657 h->SetBinError(ibin,0.0);
3658 for(
Int_t i=0;i<3;i++) {
3660 if(ixyz[i]<=nxyz[i])
break;
#define R(a, b, c, d, e, f, g, h, i)
TMatrixTSparse< Double_t > TMatrixDSparse
TMatrixT< Double_t > TMatrixD
void Set(Int_t n)
Set size of this array to n doubles.
const Double_t * GetArray() const
void Set(Int_t n)
Set size of this array to n ints.
A TGraph is an object made of two arrays X and Y with npoints each.
virtual Int_t GetNbinsY() const
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
virtual Int_t GetNbinsX() const
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Service class for 2-Dim histogram classes.
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
const TVectorD & GetEigenValues() const
const TMatrixD & GetEigenVectors() const
virtual const Int_t * GetRowIndexArray() const
virtual const Element * GetMatrixArray() const
virtual const Int_t * GetColIndexArray() const
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Class to create third splines to interpolate knots Arbitrary conditions can be introduced for first a...
Double_t Eval(Double_t x) const
Eval this spline at x.
void GetCoeff(Int_t i, Double_t &x, Double_t &y, Double_t &b, Double_t &c, Double_t &d)
Base class for spline implementation containing the Draw/Paint methods.
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
An algorithm to unfold distributions from detector to truth level.
TArrayI fHistToX
mapping of histogram bins to matrix indices
void GetRhoIJ(TH2 *rhoij, const Int_t *binMap=0) const
Get correlation coefficients, possibly cumulated over several bins.
TMatrixDSparse * fE
matrix E
void GetBias(TH1 *bias, const Int_t *binMap=0) const
Get bias vector including bias scale.
TMatrixDSparse * fEinv
matrix E^(-1)
virtual Double_t GetLcurveY(void) const
Get value on y-axis of L-curve determined in recent unfolding.
TMatrixDSparse * fAx
result x folded back A*x
TMatrixDSparse * MultiplyMSparseM(const TMatrixDSparse *a, const TMatrixD *b) const
Multiply sparse matrix and a non-sparse matrix.
virtual Double_t DoUnfold(void)
Core unfolding algorithm.
Double_t fChi2A
chi**2 contribution from (y-Ax)Vyy-1(y-Ax)
TMatrixD * fX0
bias vector x0
TMatrixDSparse * MultiplyMSparseTranspMSparse(const TMatrixDSparse *a, const TMatrixDSparse *b) const
Multiply a transposed Sparse matrix with another sparse matrix,.
TMatrixDSparse * MultiplyMSparseMSparseTranspVector(const TMatrixDSparse *m1, const TMatrixDSparse *m2, const TMatrixTBase< Double_t > *v) const
Calculate a sparse matrix product where the diagonal matrix V is given by a vector.
TMatrixDSparse * CreateSparseMatrix(Int_t nrow, Int_t ncol, Int_t nele, Int_t *row, Int_t *col, Double_t *data) const
Create a sparse matrix, given the nonzero elements.
Int_t RegularizeSize(int bin, Double_t scale=1.0)
Add a regularisation condition on the magnitude of a truth bin.
Double_t fEpsMatrix
machine accuracy used to determine matrix rank after eigenvalue analysis
void GetProbabilityMatrix(TH2 *A, EHistMap histmap) const
Get matrix of probabilities.
Double_t GetChi2L(void) const
Get contribution determined in recent unfolding.
TMatrixDSparse * fVxx
covariance matrix Vxx
Int_t GetNy(void) const
returns the number of measurement bins
virtual TString GetOutputBinName(Int_t iBinX) const
Get bin name of an output bin.
Double_t fBiasScale
scale factor for the bias
Double_t fRhoAvg
average global correlation coefficient
TMatrixDSparse * fDXDtauSquared
derivative of the result wrt tau squared
static void DeleteMatrix(TMatrixD **m)
delete matrix and invalidate pointer
void GetNormalisationVector(TH1 *s, const Int_t *binMap=0) const
Histogram of truth bins, determined from summing over the response matrix.
void ClearHistogram(TH1 *h, Double_t x=0.) const
Initialize bin contents and bin errors for a given histogram.
Int_t RegularizeDerivative(int left_bin, int right_bin, Double_t scale=1.0)
Add a regularisation condition on the difference of two truth bin.
Int_t GetNx(void) const
returns internal number of output (truth) matrix rows
static const char * GetTUnfoldVersion(void)
Return a string describing the TUnfold version.
void SetConstraint(EConstraint constraint)
Set type of area constraint.
virtual Int_t ScanLcurve(Int_t nPoint, Double_t tauMin, Double_t tauMax, TGraph **lCurve, TSpline **logTauX=0, TSpline **logTauY=0, TSpline **logTauCurvature=0)
Scan the L curve, determine tau and unfold at the final value of tau.
Int_t RegularizeBins(int start, int step, int nbin, ERegMode regmode)
Add regularisation conditions for a group of bins.
Bool_t AddRegularisationCondition(Int_t i0, Double_t f0, Int_t i1=-1, Double_t f1=0., Int_t i2=-1, Double_t f2=0.)
Add a row of regularisation conditions to the matrix L.
Int_t RegularizeCurvature(int left_bin, int center_bin, int right_bin, Double_t scale_left=1.0, Double_t scale_right=1.0)
Add a regularisation condition on the curvature of three truth bin.
void SetBias(const TH1 *bias)
Set bias vector.
void GetL(TH2 *l) const
Get matrix of regularisation conditions.
ERegMode fRegMode
type of regularisation
Int_t GetNr(void) const
Get number of regularisation conditions.
TMatrixDSparse * fVxxInv
inverse of covariance matrix Vxx-1
TMatrixD * fX
unfolding result x
EConstraint
type of extra constraint
@ kEConstraintNone
use no extra constraint
virtual Double_t GetLcurveX(void) const
Get value on x-axis of L-curve determined in recent unfolding.
TMatrixDSparse * fVyy
covariance matrix Vyy corresponding to y
Int_t fNdf
number of degrees of freedom
TArrayD fSumOverY
truth vector calculated from the non-normalized response matrix
virtual Int_t SetInput(const TH1 *hist_y, Double_t scaleBias=0.0, Double_t oneOverZeroError=0.0, const TH2 *hist_vyy=0, const TH2 *hist_vyy_inv=0)
Define input data for subsequent calls to DoUnfold(tau).
ERegMode
choice of regularisation scheme
@ kRegModeNone
no regularisation, or defined later by RegularizeXXX() methods
@ kRegModeDerivative
regularize the 1st derivative of the output distribution
@ kRegModeSize
regularise the amplitude of the output distribution
@ kRegModeCurvature
regularize the 2nd derivative of the output distribution
@ kRegModeMixed
mixed regularisation pattern
void GetInput(TH1 *inputData, const Int_t *binMap=0) const
Input vector of measurements.
void SetEpsMatrix(Double_t eps)
set numerical accuracy for Eigenvalue analysis when inverting matrices with rank problems
void GetEmatrix(TH2 *ematrix, const Int_t *binMap=0) const
Get output covariance matrix, possibly cumulated over several bins.
void ErrorMatrixToHist(TH2 *ematrix, const TMatrixDSparse *emat, const Int_t *binMap, Bool_t doClear) const
Add up an error matrix, also respecting the bin mapping.
TArrayI fXToHist
mapping of matrix indices to histogram bins
TMatrixDSparse * fDXDY
derivative of the result wrt dx/dy
TMatrixD * fY
input (measured) data y
TMatrixDSparse * InvertMSparseSymmPos(const TMatrixDSparse *A, Int_t *rank) const
Get the inverse or pseudo-inverse of a positive, sparse matrix.
Double_t GetRhoI(TH1 *rhoi, const Int_t *binMap=0, TH2 *invEmat=0) const
Get global correlation coefficients, possibly cumulated over several bins.
TMatrixDSparse * fVyyInv
inverse of the input covariance matrix Vyy-1
Double_t fLXsquared
chi**2 contribution from (x-s*x0)TLTL(x-s*x0)
TMatrixDSparse * fDXDAM[2]
matrix contribution to the of derivative dx_k/dA_ij
Double_t fTauSquared
regularisation parameter tau squared
void GetOutput(TH1 *output, const Int_t *binMap=0) const
Get output distribution, possibly cumulated over several bins.
Int_t GetNpar(void) const
Get number of truth parameters determined in recent unfolding.
virtual void ClearResults(void)
Reset all results.
Double_t fRhoMax
maximum global correlation coefficient
TMatrixDSparse * MultiplyMSparseMSparse(const TMatrixDSparse *a, const TMatrixDSparse *b) const
Multiply two sparse matrices.
EConstraint fConstraint
type of constraint to use for the unfolding
TUnfold(void)
Only for use by root streamer or derived classes.
EHistMap
arrangement of axes for the response matrix (TH2 histogram)
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
void AddMSparse(TMatrixDSparse *dest, Double_t f, const TMatrixDSparse *src) const
Add a sparse matrix, scaled by a factor, to another scaled matrix.
TMatrixDSparse * fDXDAZ[2]
vector contribution to the of derivative dx_k/dA_ij
Double_t GetRhoIFromMatrix(TH1 *rhoi, const TMatrixDSparse *eOrig, const Int_t *binMap, TH2 *invEmat) const
Get global correlation coefficients with arbitrary min map.
void InitTUnfold(void)
Initialize data members, for use in constructors.
Double_t GetTau(void) const
Return regularisation parameter.
Int_t RegularizeBins2D(int start_bin, int step1, int nbin1, int step2, int nbin2, ERegMode regmode)
Add regularisation conditions for 2d unfolding.
void GetLsquared(TH2 *lsquared) const
Get matrix of regularisation conditions squared.
void GetInputInverseEmatrix(TH2 *ematrix)
Get inverse of the measurement's covariance matrix.
TMatrixDSparse * fA
response matrix A
void GetFoldedOutput(TH1 *folded, const Int_t *binMap=0) const
Get unfolding result on detector level.
TMatrixDSparse * fL
regularisation conditions L
Int_t fIgnoredBins
number of input bins which are dropped because they have error=0
Int_t GetNdf(void) const
get number of degrees of freedom determined in recent unfolding
void swap(RDirectoryEntry &e1, RDirectoryEntry &e2) noexcept
static constexpr double m2
Short_t Max(Short_t a, Short_t b)
Int_t Finite(Double_t x)
Check if it is finite with a mask in order to be consistent in presence of fast math.
constexpr Double_t E()
Base of natural log:
Double_t Sqrt(Double_t x)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Short_t Min(Short_t a, Short_t b)
Double_t Log10(Double_t x)
static long int sum(long int i)
#define dest(otri, vertexptr)
static void output(int code)