220 for(
Int_t i=0;i<2;i++) {
322 if(rank !=
GetNx()) {
323 Warning(
"DoUnfold",
"rank of matrix E %d expect %d",rank,
GetNx());
347 epsilonNosparse(A_cols[i],0) += A_data[i];
359 Fatal(
"Unfold",
"epsilon#Eepsilon has dimension %d != 1",
366 y_minus_epsx += (*fY)(iy,0);
370 y_minus_epsx -= epsilonNosparse(ix,0) * (*fX)(ix,0);
373 lambda_half=y_minus_epsx*one_over_epsEeps;
378 if(EEpsilon_rows[ix]<EEpsilon_rows[ix+1]) {
379 (*fX)(ix,0) += lambda_half * EEpsilon_data[EEpsilon_rows[ix]];
398 data[i]=one_over_epsEeps;
408 AddMSparse(temp, -one_over_epsEeps, epsilonB);
445 if(VyyinvDy_rows[iy]<VyyinvDy_rows[iy+1]) {
446 fChi2A += VyyinvDy_data[VyyinvDy_rows[iy]]*dy(iy,0);
455 if(LsquaredDx_rows[ix]<LsquaredDx_rows[ix+1]) {
456 fLXsquared += LsquaredDx_data[LsquaredDx_rows[ix]]*dx(ix,0);
471 "epsilon#fDXDtauSquared has dimension %d != 1",
494 (Eepsilon,Eepsilon,
nullptr);
519 if(rank !=
GetNx()) {
520 Warning(
"DoUnfold",
"rank of output covariance is %d expect %d",
529 for(
int ik=VxxInv_rows[ix];ik<VxxInv_rows[ix+1];ik++) {
530 if(ix==VxxInv_cols[ik]) {
531 VxxInvDiag(ix)=VxxInv_data[ik];
545 for(
int ik=Vxx_rows[ix];ik<Vxx_rows[ix+1];ik++) {
546 if(ix==Vxx_cols[ik]) {
548 1. - 1. / (VxxInvDiag(ix) * Vxx_data[ik]);
549 if (rho_squared > rho_squared_max)
550 rho_squared_max = rho_squared;
551 if(rho_squared>0.0) {
560 fRhoAvg = (n_rho>0) ? (rho_sum/n_rho) : -1.0;
606 if(
a->GetNcols()!=
b->GetNrows()) {
607 Fatal(
"MultiplyMSparseMSparse",
608 "inconsistent matrix col/ matrix row %d !=%d",
609 a->GetNcols(),
b->GetNrows());
613 const Int_t *a_rows=
a->GetRowIndexArray();
614 const Int_t *a_cols=
a->GetColIndexArray();
615 const Double_t *a_data=
a->GetMatrixArray();
616 const Int_t *b_rows=
b->GetRowIndexArray();
617 const Int_t *b_cols=
b->GetColIndexArray();
618 const Double_t *b_data=
b->GetMatrixArray();
621 for (
Int_t irow = 0; irow <
a->GetNrows(); irow++) {
622 if(a_rows[irow+1]>a_rows[irow]) nMax +=
b->GetNcols();
624 if((nMax>0)&&(a_cols)&&(b_cols)) {
630 for (
Int_t irow = 0; irow <
a->GetNrows(); irow++) {
631 if(a_rows[irow+1]<=a_rows[irow])
continue;
633 for(
Int_t icol=0;icol<
b->GetNcols();icol++) {
637 for(
Int_t ia=a_rows[irow];ia<a_rows[irow+1];ia++) {
640 for(
Int_t ib=b_rows[k];ib<b_rows[k+1];ib++) {
641 row_data[b_cols[ib]] += a_data[ia]*b_data[ib];
645 for(
Int_t icol=0;icol<
b->GetNcols();icol++) {
646 if(row_data[icol] != 0.0) {
649 r_data[
n]=row_data[icol];
655 r->SetMatrixArray(
n,r_rows,r_cols,r_data);
680 if(
a->GetNrows() !=
b->GetNrows()) {
681 Fatal(
"MultiplyMSparseTranspMSparse",
682 "inconsistent matrix row numbers %d!=%d",
683 a->GetNrows(),
b->GetNrows());
687 const Int_t *a_rows=
a->GetRowIndexArray();
688 const Int_t *a_cols=
a->GetColIndexArray();
689 const Double_t *a_data=
a->GetMatrixArray();
690 const Int_t *b_rows=
b->GetRowIndexArray();
691 const Int_t *b_cols=
b->GetColIndexArray();
692 const Double_t *b_data=
b->GetMatrixArray();
696 typedef std::map<Int_t,Double_t> MMatrixRow_t;
697 typedef std::map<Int_t, MMatrixRow_t > MMatrix_t;
700 for(
Int_t iRowAB=0;iRowAB<
a->GetNrows();iRowAB++) {
701 for(
Int_t ia=a_rows[iRowAB];ia<a_rows[iRowAB+1];ia++) {
702 for(
Int_t ib=b_rows[iRowAB];ib<b_rows[iRowAB+1];ib++) {
704 MMatrixRow_t &row=matrix[a_cols[ia]];
705 MMatrixRow_t::iterator icol=row.find(b_cols[ib]);
706 if(icol!=row.end()) {
708 (*icol).second += a_data[ia]*b_data[ib];
711 row[b_cols[ib]] = a_data[ia]*b_data[ib];
718 for(MMatrix_t::const_iterator irow=matrix.begin();
719 irow!=matrix.end();irow++) {
720 n += (*irow).second.size();
728 for(MMatrix_t::const_iterator irow=matrix.begin();
729 irow!=matrix.end();irow++) {
730 for(MMatrixRow_t::const_iterator icol=(*irow).second.begin();
731 icol!=(*irow).second.end();icol++) {
732 r_rows[
n]=(*irow).first;
733 r_cols[
n]=(*icol).first;
734 r_data[
n]=(*icol).second;
740 r->SetMatrixArray(
n,r_rows,r_cols,r_data);
763 if(
a->GetNcols()!=
b->GetNrows()) {
764 Fatal(
"MultiplyMSparseM",
"inconsistent matrix col /matrix row %d!=%d",
765 a->GetNcols(),
b->GetNrows());
769 const Int_t *a_rows=
a->GetRowIndexArray();
770 const Int_t *a_cols=
a->GetColIndexArray();
771 const Double_t *a_data=
a->GetMatrixArray();
774 for (
Int_t irow = 0; irow <
a->GetNrows(); irow++) {
775 if(a_rows[irow+1]-a_rows[irow]>0) nMax +=
b->GetNcols();
784 for (
Int_t irow = 0; irow <
a->GetNrows(); irow++) {
785 if(a_rows[irow+1]-a_rows[irow]<=0)
continue;
786 for(
Int_t icol=0;icol<
b->GetNcols();icol++) {
790 for(
Int_t i=a_rows[irow];i<a_rows[irow+1];i++) {
792 r_data[
n] += a_data[i]*(*b)(j,icol);
794 if(r_data[
n]!=0.0)
n++;
798 r->SetMatrixArray(
n,r_rows,r_cols,r_data);
824 (
v && ((m1->
GetNcols()!=
v->GetNrows())||(
v->GetNcols()!=1)))) {
826 Fatal(
"MultiplyMSparseMSparseTranspVector",
827 "matrix cols/vector rows %d!=%d!=%d or vector rows %d!=1\n",
830 Fatal(
"MultiplyMSparseMSparseTranspVector",
839 if(rows_m1[i]<rows_m1[i+1]) num_m1++;
846 if(rows_m2[j]<rows_m2[j+1]) num_m2++;
849 const Int_t *v_rows=
nullptr;
855 Int_t num_r=num_m1*num_m2+1;
863 Int_t index_m1=rows_m1[i];
864 Int_t index_m2=rows_m2[j];
865 while((index_m1<rows_m1[i+1])&&(index_m2<rows_m2[j+1])) {
866 Int_t k1=cols_m1[index_m1];
867 Int_t k2=cols_m2[index_m2];
874 Int_t v_index=v_rows[k1];
875 if(v_index<v_rows[k1+1]) {
876 data_r[num_r] += data_m1[index_m1] * data_m2[index_m2]
880 data_r[num_r] += data_m1[index_m1] * data_m2[index_m2]
883 data_r[num_r] += data_m1[index_m1] * data_m2[index_m2];
889 if(data_r[num_r] !=0.0) {
897 num_r,row_r,col_r,data_r);
918 const Int_t *dest_rows=
dest->GetRowIndexArray();
919 const Int_t *dest_cols=
dest->GetColIndexArray();
921 const Int_t *src_rows=
src->GetRowIndexArray();
922 const Int_t *src_cols=
src->GetColIndexArray();
925 if((
dest->GetNrows()!=
src->GetNrows())||
926 (
dest->GetNcols()!=
src->GetNcols())) {
927 Fatal(
"AddMSparse",
"inconsistent matrix rows %d!=%d OR cols %d!=%d",
928 src->GetNrows(),
dest->GetNrows(),
src->GetNcols(),
dest->GetNcols());
935 for(
Int_t row=0;row<
dest->GetNrows();row++) {
936 Int_t i_dest=dest_rows[row];
937 Int_t i_src=src_rows[row];
938 while((i_dest<dest_rows[row+1])||(i_src<src_rows[row+1])) {
939 Int_t col_dest=(i_dest<dest_rows[row+1]) ?
940 dest_cols[i_dest] :
dest->GetNcols();
941 Int_t col_src =(i_src <src_rows[row+1] ) ?
942 src_cols [i_src] :
src->GetNcols();
944 if(col_dest<col_src) {
945 result_cols[
n]=col_dest;
946 result_data[
n]=dest_data[i_dest++];
947 }
else if(col_dest>col_src) {
948 result_cols[
n]=col_src;
949 result_data[
n]=
f*src_data[i_src++];
951 result_cols[
n]=col_dest;
952 result_data[
n]=dest_data[i_dest++]+
f*src_data[i_src++];
954 if(result_data[
n] !=0.0) {
956 Fatal(
"AddMSparse",
"Nan detected %d %d %d",
969 dest->SetMatrixArray(
n,result_rows,result_cols,result_data);
970 delete[] result_data;
971 delete[] result_rows;
972 delete[] result_cols;
997 Fatal(
"InvertMSparseSymmPos",
"inconsistent matrix row/col %d!=%d",
1014 for(
Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
1015 Int_t jA=a_cols[indexA];
1017 if(!(a_data[indexA]>=0.0)) nError++;
1018 aII(iA)=a_data[indexA];
1023 Fatal(
"InvertMSparseSymmPos",
1024 "Matrix has %d negative elements on the diagonal", nError);
1039 for(
Int_t iStart=0;iStart<iBlock;iStart++) {
1041 for(
Int_t i=iStart;i<iBlock;i++) {
1045 Int_t tmp=swap[iStart];
1046 swap[iStart]=swap[i];
1052 Int_t iA=swap[iStart];
1053 for(
Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
1054 Int_t jA=a_cols[indexA];
1063 for(
Int_t i=iStart+1;i<iBlock;) {
1064 if(isZero[swap[i]]) {
1068 Int_t tmp=swap[iBlock];
1069 swap[iBlock]=swap[i];
1081#ifdef CONDITION_BLOCK_PART
1083 for(
int inc=(nn+1)/2;inc;inc /= 2) {
1084 for(
int i=inc;i<nn;i++) {
1085 int itmp=swap[i+iBlock];
1087 for(j=i;(j>=inc)&&(aII(swap[j-inc+iBlock])<aII(itmp));j -=inc) {
1088 swap[j+iBlock]=swap[j-inc+iBlock];
1090 swap[j+iBlock]=itmp;
1097 swapBack[swap[i]]=i;
1101 std::cout<<
" "<<i<<
" "<<swap[i]<<
" "<<swapBack[i]<<
"\n";
1103 std::cout<<
"after sorting\n";
1105 if(i==iDiagonal) std::cout<<
"iDiagonal="<<i<<
"\n";
1106 if(i==iBlock) std::cout<<
"iBlock="<<i<<
"\n";
1107 std::cout<<
" "<<swap[i]<<
" "<<aII(swap[i])<<
"\n";
1121 for(
Int_t indexA=a_rows[row];indexA<a_rows[row+1];indexA++) {
1122 Int_t j=swapBack[a_cols[indexA]];
1124 else if((i<iDiagonal)||(j<iDiagonal)) nError1++;
1125 else if((i<iBlock)&&(j<iBlock)) nError2++;
1126 else if((i<iBlock)||(j<iBlock)) nBlock++;
1130 if(nError1+nError2>0) {
1131 Fatal(
"InvertMSparseSymmPos",
"sparse matrix analysis failed %d %d %d %d %d",
1132 iDiagonal,iBlock,A->
GetNrows(),nError1,nError2);
1137 Info(
"InvertMSparseSymmPos",
"iDiagonal=%d iBlock=%d nRow=%d",
1197 for(
Int_t i=0;i<iDiagonal;++i) {
1202 rEl_data[rNumEl]=1./aII(iA);
1207 if((!rankPtr)&&(rankD1!=iDiagonal)) {
1208 Fatal(
"InvertMSparseSymmPos",
1209 "diagonal part 1 has rank %d != %d, matrix can not be inverted",
1217 Int_t nD2=iBlock-iDiagonal;
1219 if((rNumEl>=0)&&(nD2>0)) {
1224 for(
Int_t i=0;i<nD2;++i) {
1225 Int_t iA=swap[i+iDiagonal];
1227 D2inv_col[D2invNumEl]=i;
1228 D2inv_row[D2invNumEl]=i;
1229 D2inv_data[D2invNumEl]=1./aII(iA);
1233 if(D2invNumEl==nD2) {
1236 }
else if(!rankPtr) {
1237 Fatal(
"InvertMSparseSymmPos",
1238 "diagonal part 2 has rank %d != %d, matrix can not be inverted",
1242 delete [] D2inv_data;
1243 delete [] D2inv_col;
1244 delete [] D2inv_row;
1254 if((rNumEl>=0)&&(nF>0)&&((nD2==0)||D2inv)) {
1263 Int_t nBmax=nF*(nD2+1);
1269 for(
Int_t i=0;i<nF;i++) {
1270 Int_t iA=swap[i+iBlock];
1271 for(
Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
1272 Int_t jA=a_cols[indexA];
1273 Int_t j0=swapBack[jA];
1278 F_data[FnumEl]=a_data[indexA];
1280 }
else if(j0>=iDiagonal) {
1281 Int_t j=j0-iDiagonal;
1284 B_data[BnumEl]=a_data[indexA];
1291#ifndef FORCE_EIGENVALUE_DECOMPOSITION
1311 for(
Int_t i=0;i<mbd2_nMax;i++) {
1312 mbd2_data[i] = - mbd2_data[i];
1316 if(minusBD2inv &&
F) {
1325 const Int_t *f_rows=
F->GetRowIndexArray();
1326 const Int_t *f_cols=
F->GetColIndexArray();
1327 const Double_t *f_data=
F->GetMatrixArray();
1331 for(
Int_t i=0;i<nF;i++) {
1332 for(
Int_t indexF=f_rows[i];indexF<f_rows[i+1];indexF++) {
1333 if(f_cols[indexF]>=i)
c(f_cols[indexF],i)=f_data[indexF];
1337 for(
Int_t j=0;j<i;j++) {
1348 for(
Int_t j=i+1;j<nF;j++) {
1350 for(
Int_t k=0;k<i;k++) {
1351 c_ji -=
c(i,k)*
c(j,k);
1360 for(
Int_t i=1;i<nF;i++) {
1365 std::cout<<
"dmin,dmax: "<<dmin<<
" "<<dmax<<
"\n";
1367 if(dmin/dmax<epsilonF2*nF) nErrorF=2;
1373 for(
Int_t i=0;i<nF;i++) {
1374 cinv(i,i)=1./
c(i,i);
1376 for(
Int_t i=0;i<nF;i++) {
1377 for(
Int_t j=i+1;j<nF;j++) {
1379 for(
Int_t k=i+1;k<j;k++) {
1380 tmp -= cinv(k,i)*
c(j,k);
1382 cinv(j,i)=tmp*cinv(j,j);
1387 (&cInvSparse,&cInvSparse);
1398 Int_t rankD2Block=0;
1400 if( ((nD2==0)||D2inv) && ((nF==0)||Finv) ) {
1410 if(Finv && minusBD2inv) {
1415 if(
G && minusBD2inv) {
1422 E=minusBD2invTransG;
1427 const Int_t *e_rows=E->GetRowIndexArray();
1428 const Int_t *e_cols=E->GetColIndexArray();
1429 const Double_t *e_data=E->GetMatrixArray();
1430 for(
Int_t iE=0;iE<E->GetNrows();iE++) {
1431 Int_t iA=swap[iE+iDiagonal];
1432 for(
Int_t indexE=e_rows[iE];indexE<e_rows[iE+1];++indexE) {
1433 Int_t jE=e_cols[indexE];
1434 Int_t jA=swap[jE+iDiagonal];
1437 rEl_data[rNumEl]=e_data[indexE];
1445 const Int_t *g_rows=
G->GetRowIndexArray();
1446 const Int_t *g_cols=
G->GetColIndexArray();
1447 const Double_t *g_data=
G->GetMatrixArray();
1448 for(
Int_t iG=0;iG<
G->GetNrows();iG++) {
1449 Int_t iA=swap[iG+iBlock];
1450 for(
Int_t indexG=g_rows[iG];indexG<g_rows[iG+1];++indexG) {
1451 Int_t jG=g_cols[indexG];
1452 Int_t jA=swap[jG+iDiagonal];
1456 rEl_data[rNumEl]=g_data[indexG];
1461 rEl_data[rNumEl]=g_data[indexG];
1473 Int_t iA=swap[iF+iBlock];
1474 for(
Int_t indexF=finv_rows[iF];indexF<finv_rows[iF+1];++indexF) {
1475 Int_t jF=finv_cols[indexF];
1476 Int_t jA=swap[jF+iBlock];
1479 rEl_data[rNumEl]=finv_data[indexF];
1485 }
else if(!rankPtr) {
1487 Fatal(
"InvertMSparseSymmPos",
1488 "non-trivial part has rank < %d, matrix can not be inverted",
1495 Info(
"InvertMSparseSymmPos",
1496 "cholesky-decomposition failed, try eigenvalue analysis");
1498 std::cout<<
"nEV="<<nEV<<
" iDiagonal="<<iDiagonal<<
"\n";
1501 for(
Int_t i=0;i<nEV;i++) {
1502 Int_t iA=swap[i+iDiagonal];
1503 for(
Int_t indexA=a_rows[iA];indexA<a_rows[iA+1];indexA++) {
1504 Int_t jA=a_cols[indexA];
1505 Int_t j0=swapBack[jA];
1507 Int_t j=j0-iDiagonal;
1509 if((i<0)||(j<0)||(i>=nEV)||(j>=nEV)) {
1510 std::cout<<
" error "<<nEV<<
" "<<i<<
" "<<j<<
"\n";
1514 Fatal(
"InvertMSparseSymmPos",
1515 "non-finite number detected element %d %d\n",
1518 EV(i,j)=a_data[indexA];
1525 std::cout<<
"Eigenvalues\n";
1535 if(condition<-epsilonEV2) {
1540 Error(
"InvertMSparseSymmPos",
1541 "Largest Eigenvalue is negative %f",
1544 Error(
"InvertMSparseSymmPos",
1545 "Some Eigenvalues are negative (EV%d/EV0=%g epsilon=%g)",
1555 for(
Int_t i=0;i<nEV;i++) {
1558 inverseEV(i,0)=1./
x;
1571 Int_t iA=swap[iVDVt+iDiagonal];
1572 for(
Int_t indexVDVt=vdvt_rows[iVDVt];
1573 indexVDVt<vdvt_rows[iVDVt+1];++indexVDVt) {
1574 Int_t jVDVt=vdvt_cols[indexVDVt];
1575 Int_t jA=swap[jVDVt+iDiagonal];
1578 rEl_data[rNumEl]=vdvt_data[indexVDVt];
1586 *rankPtr=rankD1+rankD2Block;
1601 rEl_row,rEl_col,rEl_data) :
nullptr;
1624 for(
Int_t i=
nullptr;i<
a.GetNrows();i++) {
1625 for(
Int_t j=
nullptr;j<
a.GetNcols();j++) {
1629 std::cout<<
"Ar is not symmetric Ar("<<i<<
","<<j<<
")="<<ar(i,j)
1630 <<
" Ar("<<j<<
","<<i<<
")="<<ar(j,i)<<
"\n";
1635 std::cout<<
"ArA is not equal A ArA("<<i<<
","<<j<<
")="<<ara(i,j)
1636 <<
" A("<<i<<
","<<j<<
")="<<
a(i,j)<<
"\n";
1641 std::cout<<
"rAr is not equal r rAr("<<i<<
","<<j<<
")="<<rar(i,j)
1642 <<
" r("<<i<<
","<<j<<
")="<<
R(i,j)<<
"\n";
1646 if(rankPtr) std::cout<<
"rank="<<*rankPtr<<
"\n";
1648 std::cout<<
"Matrix is not positive\n";
1738 for (
Int_t ix = 0; ix < nx0; ix++) {
1743 for (
Int_t iy = 0; iy < ny; iy++) {
1781 Int_t underflowBin=0,overflowBin=0;
1782 for (
Int_t ix = 0; ix < nx; ix++) {
1784 if(ibinx<1) underflowBin=1;
1786 if(ibinx>hist_A->
GetNbinsX()) overflowBin=1;
1788 if(ibinx>hist_A->
GetNbinsY()) overflowBin=1;
1793 Int_t ixfirst=-1,ixlast=-1;
1795 int nDisconnected=0;
1796 for (
Int_t ix = 0; ix < nx0; ix++) {
1807 if(((
fHistToX[ix]>=0)&&(ixlast>=0))||
1808 (nprint==nskipped)) {
1809 if(ixlast>ixfirst) {
1816 if(nprint==nskipped) {
1820 if(nskipped==(2-underflowBin-overflowBin)) {
1821 Info(
"TUnfold",
"underflow and overflow bin "
1822 "do not depend on the input data");
1824 Warning(
"TUnfold",
"%d output bins "
1825 "do not depend on the input data %s",nDisconnected,
1826 static_cast<const char *
>(binlist));
1837 for (
Int_t iy = 0; iy < ny; iy++) {
1838 for (
Int_t ix = 0; ix < nx; ix++) {
1854 if(underflowBin && overflowBin) {
1855 Info(
"TUnfold",
"%d input bins and %d output bins (includes 2 underflow/overflow bins)",ny,nx);
1856 }
else if(underflowBin) {
1857 Info(
"TUnfold",
"%d input bins and %d output bins (includes 1 underflow bin)",ny,nx);
1858 }
else if(overflowBin) {
1859 Info(
"TUnfold",
"%d input bins and %d output bins (includes 1 overflow bin)",ny,nx);
1861 Info(
"TUnfold",
"%d input bins and %d output bins",ny,nx);
1865 Error(
"TUnfold",
"too few (ny=%d) input bins for nx=%d output bins",ny,nx);
1867 Warning(
"TUnfold",
"too few (ny=%d) input bins for nx=%d output bins",ny,nx);
1879 "%d regularisation conditions have been skipped",nError);
1882 "One regularisation condition has been skipped");
1969 for(
Int_t k=l0_rows[row];k<l0_rows[row+1];k++) {
1971 l_col[nF]=l0_cols[k];
1972 l_data[nF]=l0_data[k];
1984 for(
Int_t i=0;i<nEle;i++) {
1992 l_data[nF]=rowData[i];
2116 (left_bin,-scale_left,
2117 center_bin,scale_left+scale_right,
2118 right_bin,-scale_right)
2166 Error(
"RegularizeBins",
"regmode = %d is not valid",regmode);
2168 for (
Int_t i = nSkip; i < nbin; i++) {
2205 int step2,
int nbin2,
ERegMode regmode)
2218 for (
Int_t i1 = 0; i1 < nbin1; i1++) {
2219 nError +=
RegularizeBins(start_bin + step1 * i1, step2, nbin2, regmode);
2221 for (
Int_t i2 = 0; i2 < nbin2; i2++) {
2222 nError +=
RegularizeBins(start_bin + step2 * i2, step1, nbin1, regmode);
2276 const TH2 *hist_vyy_inv)
2304 Int_t nVarianceZero=0;
2305 Int_t nVarianceForced=0;
2315 if(oneOverZeroError>0.0) {
2316 dy2 = 1./ ( oneOverZeroError*oneOverZeroError);
2328 rowVyyN[nVyyN] = iy;
2329 colVyyN[nVyyN] = iy;
2330 rowVyy1[nVyy1] = iy;
2332 dataVyyDiag[iy] = dy2;
2334 dataVyyN[nVyyN++] = dy2;
2335 dataVyy1[nVyy1++] = dy2;
2340 Int_t nOffDiagNonzero=0;
2343 if(dataVyyDiag[iy]<=0.0) {
2353 if(iy==jy)
continue;
2355 if(dataVyyDiag[jy]<=0.0)
continue;
2357 rowVyyN[nVyyN] = iy;
2358 colVyyN[nVyyN] = jy;
2360 if(dataVyyN[nVyyN] == 0.0)
continue;
2366 "inverse of input covariance is taken from user input");
2373 rowVyyInv[nVyyInv] = iy;
2374 colVyyInv[nVyyInv] = jy;
2375 dataVyyInv[nVyyInv]= hist_vyy_inv->
GetBinContent(iy+1,jy+1);
2376 if(dataVyyInv[nVyyInv] == 0.0)
continue;
2381 (
GetNy(),
GetNy(),nVyyInv,rowVyyInv,colVyyInv,dataVyyInv);
2382 delete [] rowVyyInv;
2383 delete [] colVyyInv;
2384 delete [] dataVyyInv;
2386 if(nOffDiagNonzero) {
2388 "input covariance has elements C(X,Y)!=nullptr where V(X)==0");
2394 (
GetNy(),
GetNy(),nVyyN,rowVyyN,colVyyN,dataVyyN);
2401 (
GetNy(),1,nVyy1,rowVyy1,colVyy1, dataVyy1);
2412 (*fY) (i, 0) =
input->GetBinContent(i + 1);
2426 if(nVarianceForced) {
2427 if(nVarianceForced>1) {
2428 Warning(
"SetInput",
"%d/%d input bins have zero error,"
2429 " 1/error set to %lf.",
2430 nVarianceForced,
GetNy(),oneOverZeroError);
2432 Warning(
"SetInput",
"One input bin has zero error,"
2433 " 1/error set to %lf.",oneOverZeroError);
2437 if(nVarianceZero>1) {
2438 Warning(
"SetInput",
"%d/%d input bins have zero error,"
2439 " and are ignored.",nVarianceZero,
GetNy());
2441 Warning(
"SetInput",
"One input bin has zero error,"
2442 " and is ignored.");
2448 if(oneOverZeroError<=0.0) {
2454 TString binlist(
"no data to constrain output bin ");
2465 Warning(
"SetInput",
"%s",(
char const *)binlist);
2470 Error(
"SetInput",
"%d/%d output bins are not constrained by any data.",
2473 Error(
"SetInput",
"One output bin [%d] is not constrained by any data.",
2479 delete[] dataVyyDiag;
2481 return nVarianceForced+nVarianceZero+10000*nError2;
2515 for(
Int_t irow=0;irow<nRow;irow++) {
2516 for(
Int_t indexE=rows[irow];indexE<rows[irow+1];indexE++) {
2517 Int_t icol=cols[indexE];
2518 eSym(irow,icol)=
data[indexE];
2523 for(
int i=0;i<
r.GetNrows();i++) {
2563 typedef std::map<Double_t,std::pair<Double_t,Double_t> > XYtau_t;
2588 if((tauMin<=0)||(tauMax<=0.0)||(tauMin>=tauMax)) {
2598 Error(
"ScanLcurve",
"too few input bins, NDF<=nullptr %d",
GetNdf());
2603 Info(
"ScanLcurve",
"logtau=-Infinity X=%lf Y=%lf",x0,y0);
2605 Fatal(
"ScanLcurve",
"problem (too few input bins?) X=%f",x0);
2608 Fatal(
"ScanLcurve",
"problem (missing regularisation?) Y=%f",y0);
2618 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2622 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2625 if((*curve.begin()).second.first<x0) {
2634 Double_t logTau=(*curve.begin()).first-0.5;
2637 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2641 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2644 while(((
int)curve.size()<(nPoint-1)/2)&&
2645 ((*curve.begin()).second.first<x0));
2652 while(((
int)curve.size()<nPoint-1)&&
2653 (((*curve.begin()).second.first-x0>0.00432)||
2654 ((*curve.begin()).second.second-y0>0.00432)||
2655 (curve.size()<2))) {
2656 Double_t logTau=(*curve.begin()).first-0.5;
2659 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2663 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2674 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2677 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2684 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2687 Info(
"ScanLcurve",
"logtau=%lf X=%lf Y=%lf",
2697 while(
int(curve.size())<nPoint-1) {
2700 XYtau_t::const_iterator i0,i1;
2705 for(i1++;i1!=curve.end();i1++) {
2706 const std::pair<Double_t,Double_t> &xy0=(*i0).second;
2707 const std::pair<Double_t,Double_t> &xy1=(*i1).second;
2713 logTau=0.5*((*i0).first+(*i1).first);
2719 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2731 XYtau_t::const_iterator i0,i1;
2736 if( ((
int)curve.size())<nPoint) {
2745 for( XYtau_t::const_iterator i=curve.begin();i!=curve.end();i++) {
2746 lXi[
n]=(*i).second.first;
2747 lYi[
n]=(*i).second.second;
2755 for(
Int_t i=0;i<
n-1;i++) {
2758 Double_t tauBar=0.5*(lTi[i]+lTi[i+1]);
2759 Double_t dTau=0.5*(lTi[i+1]-lTi[i]);
2760 Double_t dx1=bi+dTau*(2.*ci+3.*di*dTau);
2763 Double_t dy1=bi+dTau*(2.*ci+3.*di*dTau);
2766 cCi[i]=(dy2*dx1-dy1*dx2)/
TMath::Power(dx1*dx1+dy1*dy1,1.5);
2786 for(
Int_t i=iskip;i<
n-2-iskip;i++) {
2809 xx = m_p_half + discr;
2811 xx = m_p_half - discr;
2815 if((xx>0.0)&&(xx<dx)) {
2829 if((xx>0.0)&&(xx<dx)) {
2843 if(logTauCurvature) {
2844 *logTauCurvature=splineC;
2853 Fatal(
"ScanLcurve",
"problem (missing regularisation?) X=%f Y=%f",
2856 Info(
"ScanLcurve",
"Result logtau=%lf X=%lf Y=%lf",
2866 Int_t bestChoice=-1;
2867 if(!curve.empty()) {
2872 for( XYtau_t::const_iterator i=curve.begin();i!=curve.end();i++) {
2873 if(logTauFin==(*i).first) {
2876 x[
n]=(*i).second.first;
2877 y[
n]=(*i).second.second;
2883 (*lCurve)->SetTitle(
"L curve");
2885 if(logTauX) (*logTauX)=
new TSpline3(
"log(chi**2)%log(tau)",logT,
x,
n);
2886 if(logTauY) (*logTauY)=
new TSpline3(
"log(reg.cond)%log(tau)",logT,
y,
n);
2943 out->GetBinContent(
dest));
2976 Int_t destI = binMap ? binMap[i+1] : i+1;
2977 if(destI<0)
continue;
2979 out->SetBinContent(destI, (*
fAx) (i, 0)+ out->GetBinContent(destI));
2981 Int_t index_a=rows_A[i];
2982 Int_t index_av=rows_AVxx[i];
2983 while((index_a<rows_A[i+1])&&(index_av<rows_AVxx[i+1])) {
2984 Int_t j_a=cols_A[index_a];
2985 Int_t j_av=cols_AVxx[index_av];
2988 }
else if(j_a>j_av) {
2991 e2 += data_AVxx[index_av] * data_A[index_a];
3019 for(
Int_t indexA=rows_A[iy];indexA<rows_A[iy+1];indexA++) {
3020 Int_t ix = cols_A[indexA];
3046 for(
Int_t indexDXDY=rows_dxdy[ix];indexDXDY<rows_dxdy[ix+1];indexDXDY++) {
3047 Int_t iy = cols_dxdy[indexDXDY];
3078 Int_t destI=binMap ? binMap[i+1] : i+1;
3079 if(destI<0)
continue;
3081 out->SetBinContent(destI, (*
fY) (i, 0)+out->GetBinContent(destI));
3085 if(cols_Vyy[
index]==i) {
3089 out->SetBinError(destI,
e);
3109 Warning(
"GetInputInverseEmatrix",
"input covariance matrix has rank %d expect %d",
3113 Error(
"GetInputInverseEmatrix",
"number of parameters %d > %d (rank of input covariance). Problem can not be solved",
GetNpar(),rank);
3114 }
else if(
fNdf==0) {
3115 Warning(
"GetInputInverseEmatrix",
"number of parameters %d = input rank %d. Problem is ill posed",
GetNpar(),rank);
3124 for(
int i=0;i<=out->GetNbinsX()+1;i++) {
3125 for(
int j=0;j<=out->GetNbinsY()+1;j++) {
3126 out->SetBinContent(i,j,0.);
3133 out->SetBinContent(i+1,j+1,data_VyyInv[
index]);
3201 out->SetBinContent(indexH,row+1,
data[
cindex]);
3309 std::map<Int_t,Double_t> e2;
3316 for(
Int_t i=0;i<binMapSize;i++) {
3317 Int_t destBinI=binMap ? binMap[i] : i;
3319 if((destBinI>=0)&&(srcBinI>=0)) {
3321 (destBinI, (*
fX)(srcBinI,0)+
output->GetBinContent(destBinI));
3327 Int_t index_vxx=rows_Vxx[srcBinI];
3329 while((j<binMapSize)&&(index_vxx<rows_Vxx[srcBinI+1])) {
3330 Int_t destBinJ=binMap ? binMap[j] : j;
3331 if(destBinI!=destBinJ) {
3340 if(cols_Vxx[index_vxx]<srcBinJ) {
3343 }
else if(cols_Vxx[index_vxx]>srcBinJ) {
3348 e2[destBinI] += data_Vxx[index_vxx];
3358 for(std::map<Int_t,Double_t>::const_iterator i=e2.begin();
3384 for(
Int_t i=0;i<nbin+2;i++) {
3385 for(
Int_t j=0;j<nbin+2;j++) {
3398 for(
Int_t i=0;i<binMapSize;i++) {
3399 Int_t destBinI=binMap ? binMap[i] : i;
3401 if((destBinI>=0)&&(destBinI<nbin+2)&&(srcBinI>=0)) {
3407 Int_t index_vxx=rows_emat[srcBinI];
3408 while((j<binMapSize)&&(index_vxx<rows_emat[srcBinI+1])) {
3409 Int_t destBinJ=binMap ? binMap[j] : j;
3411 if((destBinJ<0)||(destBinJ>=nbin+2)||(srcBinJ<0)) {
3415 if(cols_emat[index_vxx]<srcBinJ) {
3418 }
else if(cols_emat[index_vxx]>srcBinJ) {
3424 + data_emat[index_vxx];
3466 for(
Int_t i=0;i<nbin+2;i++) {
3469 for(
Int_t i=0;i<nbin+2;i++) {
3470 for(
Int_t j=0;j<nbin+2;j++) {
3471 if((
e[i]>0.0)&&(
e[j]>0.0)) {
3527 for(
int index_vxx=rows_Vxx[i];index_vxx<rows_Vxx[i+1];index_vxx++) {
3528 if(cols_Vxx[index_vxx]==i) {
3529 e_ii=data_Vxx[index_vxx];
3533 for(
int index_vxxinv=rows_VxxInv[i];index_vxxinv<rows_VxxInv[i+1];
3535 if(cols_VxxInv[index_vxxinv]==i) {
3536 einv_ii=data_VxxInv[index_vxxinv];
3541 if((einv_ii>0.0)&&(e_ii>0.0)) rho=1.-1./(einv_ii*e_ii);
3554 const Int_t *binMap,
TH2 *invEmat)
const
3577 std::map<int,int> histToLocalBin;
3579 for(
Int_t i=0;i<binMapSize;i++) {
3580 Int_t mapped_i=binMap ? binMap[i] : i;
3582 if(histToLocalBin.find(mapped_i)==histToLocalBin.end()) {
3583 histToLocalBin[mapped_i]=nBin;
3593 for(std::map<int,int>::const_iterator i=histToLocalBin.begin();
3594 i!=histToLocalBin.end();i++) {
3595 localBinToHist[(*i).second]=(*i).first;
3612 Int_t destI=binMap ? binMap[origI] : origI;
3613 if(destI<0)
continue;
3614 Int_t ematBinI=histToLocalBin[destI];
3615 for(
int index_eOrig=rows_eOrig[i];index_eOrig<rows_eOrig[i+1];
3618 Int_t j=cols_eOrig[index_eOrig];
3622 Int_t destJ=binMap ? binMap[origJ] : origJ;
3623 if(destJ<0)
continue;
3624 Int_t ematBinJ=histToLocalBin[destJ];
3625 eRemap(ematBinI,ematBinJ) += data_eOrig[index_eOrig];
3633 Warning(
"GetRhoIFormMatrix",
"Covariance matrix has rank %d expect %d",
3644 for(
Int_t index_einv=rows_eInv[i];index_einv<rows_eInv[i+1];
3646 Int_t j=cols_eInv[index_einv];
3648 einv_ii=data_eInv[index_einv];
3652 data_eInv[index_einv]);
3656 if((einv_ii>0.0)&&(e_ii>0.0)) rho=1.-1./(einv_ii*e_ii);
3667 delete [] localBinToHist;
3683 nxyz[0]=
h->GetNbinsX()+1;
3684 nxyz[1]=
h->GetNbinsY()+1;
3685 nxyz[2]=
h->GetNbinsZ()+1;
3686 for(
int i=
h->GetDimension();i<3;i++) nxyz[i]=0;
3688 for(
int i=0;i<3;i++) ixyz[i]=0;
3689 while((ixyz[0]<=nxyz[0])&&
3690 (ixyz[1]<=nxyz[1])&&
3691 (ixyz[2]<=nxyz[2])){
3692 Int_t ibin=
h->GetBin(ixyz[0],ixyz[1],ixyz[2]);
3693 h->SetBinContent(ibin,
x);
3694 h->SetBinError(ibin,0.0);
3695 for(
Int_t i=0;i<3;i++) {
3697 if(ixyz[i]<=nxyz[i])
break;
3756 for(
Int_t indexA=rows_A[iy];indexA<rows_A[iy+1];indexA++) {
3757 Int_t ix = cols_A[indexA];
3758 r += data_A[indexA]*(*fDXDY)(ix,iy);
3790 double DF,SURE,chi2A,
x,
y;
3792 std::map<double,ScanResult > scan;
3797 while((
int)scan.size()<nPoint) {
3802 if((tauMin>=0.)&&(tauMax>tauMin)) {
3808 if((tauMin>=0.)&&(tauMax>tauMin)) {
3813 if(tau<=tauMin) tau=tauMin+1.;
3817 std::vector<double> t,s;
3820 for(std::map<double,ScanResult>::const_iterator i=scan.begin();
3821 i!=scan.end();i++) {
3822 t.push_back((*i).first);
3823 s.push_back((*i).second.SURE);
3825 if((
n<nPoint-1)||(
n<3)) {
3832 double s0=0.,
s1=0.,s2=0.;
3834 for(
size_t i=0;i<t.size()-1;i++) {
3847 double diffSq=s2-mean*mean*
s0;
3848 if(diffSq<0.5*(rMax-mean)) {
3853 if((where==0)&&(t[where]<=0.0)) {
3864 for(
size_t i=2;i<t.size()-1;i++) {
3865 if(s[i]<s[where]) where=i;
3870 for(
int i=-1;i<2;i++) t[where+i]=
TMath::Log10(t[where+i]);
3872 Info(
"ScanSURE",
"minimum near: [%f,%f,%f] -> [%f,%f,%f}",
3873 t[where-1],t[where],t[where+1],s[where-1],s[where],s[where+1]);
3874 double s21=s[where+1]-s[where];
3875 double s01=s[where-1]-s[where];
3876 double t21=t[where+1]-t[where];
3877 double t01=t[where-1]-t[where];
3878 double t20=t[where+1]-t[where-1];
3879 double ct0t2=(s21*t01-s01*t21)/t20;
3880 double bt0t2=(s01*t21*t21-s21*t01*t01)/t20;
3882 tau=t[where]-0.5*bt0t2/ct0t2;
3887 tau=t[where]+0.5*t21;
3889 tau=t[where]+0.5*t01;
3897 ScanResult &
r=scan[tau];
3904 if((tau<=0.)&&(
GetNdf()<=0)) {
3905 Error(
"ScanSURE",
"too few input bins, NDF<=nullptr %d",
GetNdf());
3908 Info(
"ScanSURE",
"logtau=-Infinity Chi2A=%lf SURE=%lf DF=%lf X=%lf Y=%lf",
3909 r.chi2A,
r.SURE,
r.DF,
r.x,
r.y);
3911 Fatal(
"ScanSURE",
"problem (too few input bins?) x=%f",
r.x);
3914 Fatal(
"ScanSURE",
"problem (missing regularisation?) y=%f",
r.y);
3917 Info(
"ScanSURE",
"logtau=%lf Chi2A=%lf SURE=%lf DF=%lf X=%lf Y=%lf",
3925 std::vector<double> logTau,SURE,DF,chi2A,
X,Y;
3926 for(std::map<double,ScanResult>::const_iterator i=scan.begin();
3927 i!=scan.end();i++) {
3930 if((*i).first>0.0) {
3935 double s=(*i).second.SURE;
3936 if((iBest<0)||(SURE[iBest]>s)) {
3939 logTau.push_back(log10tau);
3941 DF.push_back((*i).second.DF);
3942 chi2A.push_back((*i).second.chi2A);
3943 X.push_back((*i).second.x);
3944 Y.push_back((*i).second.y);
3948 int n=logTau.size();
3952 *logTauSURE=
new TGraph(
n,logTau.data(),SURE.data());
3955 *df_chi2A=
new TGraph(
n,DF.data(),chi2A.data());
3958 *lCurve=
new TGraph(
n,
X.data(),Y.data());
#define R(a, b, c, d, e, f, g, h, i)
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void input
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t dest
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint xy
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t src
TMatrixTSparse< Double_t > TMatrixDSparse
TMatrixT< Double_t > TMatrixD
void Set(Int_t n) override
Set size of this array to n doubles.
const Double_t * GetArray() const
void Set(Int_t n) override
Set size of this array to n ints.
A TGraph is an object made of two arrays X and Y with npoints each.
TH1 is the base class of all histogram classes in ROOT.
virtual Int_t GetNbinsY() const
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-D histogram classes.
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
Double_t GetBinContent(Int_t binx, Int_t biny) const override
const TVectorD & GetEigenValues() const
const TMatrixD & GetEigenVectors() const
TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *="") override
Copy array data to matrix .
const Int_t * GetRowIndexArray() const override
const Int_t * GetColIndexArray() const override
const Element * GetMatrixArray() const override
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...
void GetCoeff(Int_t i, Double_t &x, Double_t &y, Double_t &b, Double_t &c, Double_t &d) const
Double_t Eval(Double_t x) const override
Eval this spline at x.
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
double GetDF(void) const
return the effecive number of degrees of freedom See e.g.
TMatrixDSparse * fE
matrix E
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
double GetSURE(void) const
return Stein's unbiased risk estimator See e.g.
void GetBias(TH1 *bias, const Int_t *binMap=nullptr) const
get bias vector including bias scale
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 M1*V*M2T 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 χ2L 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 outpt bin.
Double_t fBiasScale
scale factor for the bias
virtual Int_t ScanSURE(Int_t nPoint, Double_t tauMin, Double_t tauMax, TGraph **logTauSURE=nullptr, TGraph **df_chi2A=nullptr, TGraph **lCurve=nullptr)
minimize Stein's unbiased risk estimator "SURE" using successive calls to DoUnfold at various tau.
virtual Int_t ScanLcurve(Int_t nPoint, Double_t tauMin, Double_t tauMax, TGraph **lCurve, TSpline **logTauX=nullptr, TSpline **logTauY=nullptr, TSpline **logTauCurvature=nullptr)
scan the L curve, determine tau and unfold at the final value of tau
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 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
void GetFoldedOutput(TH1 *folded, const Int_t *binMap=nullptr) const
get unfolding result on detector level
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
Double_t GetRhoI(TH1 *rhoi, const Int_t *binMap=nullptr, TH2 *invEmat=nullptr) const
get global correlation coefficiencts, possibly cumulated over several bins
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
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=nullptr) const
Input vector of measurements.
void SetEpsMatrix(Double_t eps)
set numerical accuracy for Eigenvalue analysis when inverting matrices with rank problems
const TMatrixDSparse * GetE(void) const
matrix E, using internal bin counting
TVectorD GetSqrtEvEmatrix(void) const
void GetOutput(TH1 *output, const Int_t *binMap=nullptr) const
get output distribution, possibly cumulated over several bins
void GetRhoIJ(TH2 *rhoij, const Int_t *binMap=nullptr) const
get correlation coefficiencts, 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
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
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
void GetEmatrix(TH2 *ematrix, const Int_t *binMap=nullptr) const
get output covariance matrix, possibly cumulated over several bins
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
void GetNormalisationVector(TH1 *s, const Int_t *binMap=nullptr) const
histogram of truth bins, determined from suming over the response 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
void InitTUnfold(void)
initialize data menbers, for use in constructors
Double_t GetTau(void) const
return regularisation parameter
Double_t GetChi2A(void) const
get χ2A contribution determined in recent unfolding
Int_t RegularizeBins2D(int start_bin, int step1, int nbin1, int step2, int nbin2, ERegMode regmode)
add regularisation conditions for 2d unfolding
const TMatrixDSparse * GetDXDY(void) const
matrix of derivatives dx/dy
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
TMatrixDSparse * fL
regularisation conditions L
virtual Int_t SetInput(const TH1 *hist_y, Double_t scaleBias=0.0, Double_t oneOverZeroError=0.0, const TH2 *hist_vyy=nullptr, const TH2 *hist_vyy_inv=nullptr)
Define input data for subsequent calls to DoUnfold(tau)
Int_t fIgnoredBins
number of input bins which are dropped because they have error=nullptr
Int_t GetNdf(void) const
get number of degrees of freedom determined in recent unfolding
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and 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.
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Double_t Sqrt(Double_t x)
Returns the square root of x.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Double_t Log10(Double_t x)
Returns the common (base-10) logarithm of x.
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
static uint64_t sum(uint64_t i)