diff --git a/hist/hist/inc/TUnfold.h b/hist/hist/inc/TUnfold.h
index 779e4fadc915a2a9f28647798888e0f045257e56..df338532a5bb2c65071dbd893e093e7463b594ea 100644
--- a/hist/hist/inc/TUnfold.h
+++ b/hist/hist/inc/TUnfold.h
@@ -1,9 +1,14 @@
// Author: Stefan Schmitt
// DESY, 13/10/08
-// Version 6, completely remove definition of class XY
+// Version 11, regularisation methods have return values
//
// History:
+// Version 10, with bug-fix in TUnfold.cxx
+// Version 9, implements method for optimized inversion of sparse matrix
+// Version 8, replace all TMatrixSparse matrix operations by private code
+// Version 7, fix problem with TMatrixDSparse,TMatrixD multiplication
+// Version 6, completely remove definition of class XY
// Version 5, move definition of class XY from TUnfold.C to this file
// Version 4, with bug-fix in TUnfold.C
// Version 3, with bug-fix in TUnfold.C
@@ -75,11 +80,16 @@ class TUnfold:public TObject {
Double_t fChi2L; // Result: chi**2 contribution from tau(x-s*x0)Lsquared(x-s*x0)
Double_t fRhoMax; // Result: maximum global correlation
Double_t fRhoAvg; // Result: average global correlation
+ Int_t fNdf; // Result: number of degrees of freedom
protected:
TUnfold(void); // for derived classes
virtual Double_t DoUnfold(void); // the unfolding algorithm
virtual void CalculateChi2Rho(void); // supplementory calculations
- static TMatrixDSparse *CreateSparseMatrix(Int_t nr, Int_t nc, Int_t * row, Int_t * col, Double_t const *data); // create matrices
+ static TMatrixDSparse *MultiplyMSparseM(TMatrixDSparse const &a,TMatrixD const &b); // multiply sparse and non-sparse matrix
+ static TMatrixDSparse *MultiplyMSparseMSparse(TMatrixDSparse const &a,TMatrixDSparse const &b); // multiply sparse and sparse matrix
+ static TMatrixDSparse *MultiplyMSparseTranspMSparse(TMatrixDSparse const &a,TMatrixDSparse const &b); // multiply transposed sparse and sparse matrix
+ static Double_t MultiplyVecMSparseVec(TMatrixDSparse const &a,TMatrixD const &v); // scalar product of v and Av
+ static TMatrixD *InvertMSparse(TMatrixDSparse const &A); // invert sparse matrix
inline Int_t GetNx(void) const {
return fA->GetNcols();
} // number of non-zero output bins
@@ -96,19 +106,19 @@ class TUnfold:public TObject {
kRegModeNone = 0, // no regularisation
kRegModeSize = 1, // regularise the size of the output
kRegModeDerivative = 2, // regularize the 1st derivative of the output
- kRegModeCurvature = 3 // regularize the 2nd derivative of the output
+ kRegModeCurvature = 3, // regularize the 2nd derivative of the output
};
TUnfold(TH2 const *hist_A, EHistMap histmap, ERegMode regmode = kRegModeSize); // constructor
virtual ~ TUnfold(void); // delete data members
void SetBias(TH1 const *bias); // set alternative bias
- void RegularizeSize(int bin, Double_t const &scale = 1.0); // regularise the size of one output bin
- void RegularizeDerivative(int left_bin, int right_bin, Double_t const &scale = 1.0); // regularize difference of two output bins (1st derivative)
- void RegularizeCurvature(int left_bin, int center_bin, int right_bin, Double_t const &scale_left = 1.0, Double_t const &scale_right = 1.0); // regularize curvature of three output bins (2nd derivative)
- void RegularizeBins(int start, int step, int nbin, ERegMode regmode); // regularize a 1-dimensional curve
- void RegularizeBins2D(int start_bin, int step1, int nbin1, int step2, int nbin2, ERegMode regmode); // regularize a 2-dimensional grid
+ Int_t RegularizeSize(int bin, Double_t const &scale = 1.0); // regularise the size of one output bin
+ Int_t RegularizeDerivative(int left_bin, int right_bin, Double_t const &scale = 1.0); // regularize difference of two output bins (1st derivative)
+ Int_t RegularizeCurvature(int left_bin, int center_bin, int right_bin, Double_t const &scale_left = 1.0, Double_t const &scale_right = 1.0); // regularize curvature of three output bins (2nd derivative)
+ Int_t RegularizeBins(int start, int step, int nbin, ERegMode regmode); // regularize a 1-dimensional curve
+ Int_t RegularizeBins2D(int start_bin, int step1, int nbin1, int step2, int nbin2, ERegMode regmode); // regularize a 2-dimensional grid
Double_t DoUnfold(Double_t const &tau,
TH1 const *hist_y, Double_t const &scaleBias=0.0); // do the unfolding
- void SetInput(TH1 const *hist_y, Double_t const &scaleBias=0.0); // define input distribution for ScanLCurve
+ Int_t SetInput(TH1 const *hist_y, Double_t const &scaleBias=0.0,Double_t oneOverZeroError=0.0); // define input distribution for ScanLCurve
virtual Double_t DoUnfold(Double_t const &tau); // Unfold with given choice of tau
virtual Int_t ScanLcurve(Int_t nPoint,Double_t const &tauMin,
Double_t const &tauMax,TGraph **lCurve,
@@ -134,8 +144,10 @@ class TUnfold:public TObject {
Double_t const &GetChi2L(void) const; // chi**2 contribution from L
Double_t GetLcurveX(void) const; // x axis of L curve
Double_t GetLcurveY(void) const; // y axis of L curve
+ Int_t GetNdf(void) const; // number of degrees of freedom
+ Int_t GetNpar(void) const; // number of parameters
- ClassDef(TUnfold, 0)
+ ClassDef(TUnfold, 0) //Unfolding with support for L-curve analysis
};
#endif
diff --git a/hist/hist/src/TUnfold.cxx b/hist/hist/src/TUnfold.cxx
index e3706daecc5f745d72ece04a531fdf651551bd9a..b70b0fc429bc4e17bead3eb87247809904305898 100644
--- a/hist/hist/src/TUnfold.cxx
+++ b/hist/hist/src/TUnfold.cxx
@@ -1,10 +1,14 @@
// Author: Stefan Schmitt
// DESY, 13/10/08
-
-// Version 6, replace class XY by std::pair
+// Version 11, reduce the amount of printout
//
// History:
+// Version 10, more correct definition of the L curve, update references
+// Version 9, faster matrix inversion and skip edge points for L-curve scan
+// Version 8, replace all TMatrixSparse matrix operations by private code
+// Version 7, fix problem with TMatrixDSparse,TMatrixD multiplication
+// Version 6, replace class XY by std::pair
// Version 5, replace run-time dynamic arrays by new and delete[]
// Version 4, fix new bug from V3 with initial regularisation condition
// Version 3, fix bug with initial regularisation condition
@@ -55,9 +59,19 @@
//
// References:
// ===========
-// Talk by V. Blobel, Terascale Statistics school
+// A nice overview of the method is given in:
+// The L-curve and Its Use in the Numerical Treatment of Inverse Problems
+// (2000) by P. C. Hansen, in Computational Inverse Problems in
+// Electrocardiology, ed. P. Johnston,
+// Advances in Computational Bioengineering
+// http://www.imm.dtu.dk/~pch/TR/Lcurve.ps
+// The relevant equations are (1), (2) for the unfolding
+// and (14) for the L-curve curvature definition
+//
+// Related literature on unfolding:
+// Talk by V. Blobel, Terascale Statistics school
// https://indico.desy.de/contributionDisplay.py?contribId=23&confId=1149
-// References quoted in Blobel's talk:
+// References quoted in Blobel's talk:
// Per Chistian Hansen, Rank-Deficient and Discrete Ill-posed Problems,
// Siam (1998)
// Jari Kaipio and Erkki Somersalo, Statistical and Computational
@@ -274,6 +288,7 @@
// RegularizeDerivative() regularize the slope given by two bins
// RegularizeCurvature() regularize the curvature given by three bins
+
#include <iostream>
#include <TMatrixD.h>
#include <TMatrixDSparse.h>
@@ -282,6 +297,22 @@
#include "TUnfold.h"
#include <map>
+#include <vector>
+
+
+// this option enables the use of Schur complement for matrix inversion
+// with large diagonal sub-matrices
+#define SCHUR_COMPLEMENT_MATRIX_INVERSION
+
+// this option saves the spline of the L curve curvature to a file
+// named splinec.ps for debugging
+
+// #define DEBUG_LCURVE
+
+
+#ifdef DEBUG_LCURVE
+#include <TCanvas.h>
+#endif
ClassImp(TUnfold)
//______________________________________________________________________________
@@ -355,48 +386,48 @@ Double_t TUnfold::DoUnfold(void)
// T
// fA fV = mAt_V
//
- TMatrixDSparse mAt_V(TMatrixDSparse(TMatrixDSparse::kTransposed, *fA),
- TMatrixDSparse::kMult, *fV);
+ TMatrixDSparse *mAt_V=MultiplyMSparseTranspMSparse(*fA,*fV);
//
// get T
// fA fV fY + fTau fBiasScale Lsquared fX0 = rhs
//
- TMatrixDSparse rhs(mAt_V, TMatrixDSparse::kMult, *fY);
+ TMatrixDSparse *rhs=MultiplyMSparseM(*mAt_V,*fY);
if (fBiasScale != 0.0) {
- rhs +=
- (fTau * fBiasScale) * TMatrixDSparse(*fLsquared,
- TMatrixDSparse::kMult,
- *fX0);
+ TMatrixDSparse *rhs2=MultiplyMSparseM(*fLsquared,*fX0);
+ (*rhs) += (fTau * fBiasScale) * (*rhs2);
+ delete rhs2;
}
//
// get matrix
// T
// (fA fV)fA + fTau*fLsquared = fEinv
//
- fEinv = new TMatrixDSparse(mAt_V, TMatrixDSparse::kMult, *fA);
- *fEinv += fTau * (*fLsquared);
+ TMatrixDSparse *tmp=MultiplyMSparseMSparse(*mAt_V,*fA);
+ fEinv=new TMatrixDSparse(*tmp,TMatrixDSparse::kPlus,fTau * (*fLsquared));
+ delete tmp;
//
// get matrix
// -1
// fEinv = fE
//
-#ifdef MEASURE_TIMING
- struct timespec t0, t1;
- clockid_t type;
- clock_getcpuclockid(0, &type);
- clock_gettime(type, &t0);
-#endif
- fE = new TMatrixD(TMatrixD::kInverted, *fEinv);
+ fE = InvertMSparse(*fEinv);
//
// get result
// fE rhs = x
//
- fX = new TMatrixD(*fE, TMatrixD::kMult, rhs);
+ fX = new TMatrixD(*fE, TMatrixD::kMult, *rhs);
//
// get result
// fA x = fAx
//
- fAx = new TMatrixDSparse(*fA, TMatrixDSparse::kMult, *fX);
+ fAx = MultiplyMSparseM(*fA,*fX);
+
+ //
+ // clean up
+ //
+ delete rhs;
+ delete mAt_V;
+
//
// calculate chi**2 etc
//
@@ -414,24 +445,10 @@ void TUnfold::CalculateChi2Rho(void)
// fChi2A,fChi2L,fRhoMax,fRhoAvg
// chi**2 contribution from (y-Ax)V(y-Ax)
- fChi2A = 0.0;
TMatrixD dy(*fY, TMatrixD::kMinus, *fAx);
- // Vd is made sparse, because it is made from sparse matrix V
- TMatrixDSparse Vd(*fV, TMatrixDSparse::kMult, dy);
- for (int iy = 0; iy < GetNy(); iy++) {
- fChi2A += dy(iy, 0) * Vd(iy, 0);
- }
-
- // chi**2 contribution from tau(x-s*x0)Lsquared(x-s*x0)
- fChi2L = 0.0;
+ fChi2A = MultiplyVecMSparseVec(*fV,dy);
TMatrixD dx(*fX, TMatrixD::kMinus, fBiasScale * (*fX0));
- // LsquaredDx is made sparse, because it is made from sparse matrix
- // fLsquared
- TMatrixDSparse LsquaredDx(*fLsquared, TMatrixDSparse::kMult, dx);
- for (int ix = 0; ix < GetNx(); ix++) {
- fChi2L += fTau * dx(ix, 0) * LsquaredDx(ix, 0);
- }
-
+ fChi2L = fTau*MultiplyVecMSparseVec(*fLsquared,dx);
// maximum global correlation coefficient
Double_t rho_squared_max = 0.0;
Double_t rho_sum = 0.0;
@@ -448,30 +465,423 @@ void TUnfold::CalculateChi2Rho(void)
}
fRhoMax = TMath::Sqrt(rho_squared_max);
fRhoAvg = (n_rho>0) ? (rho_sum/n_rho) : -1.0;
+
}
-TMatrixDSparse *TUnfold::CreateSparseMatrix(Int_t nr, Int_t nc,
- Int_t * row, Int_t * col,
- Double_t const *data)
+Double_t TUnfold::MultiplyVecMSparseVec(TMatrixDSparse const &a,TMatrixD const &v)
{
- // Create a sparse matrix, to set up fA or fV.
- // nr,nc: number of rows and columns
- // row: row index array
- // col: column index array
- // data: non-zero matrix elements
- // Return value:
- // a new matrix
- //
- // This implementation cirumvents certain problems with TMatrixDSparse
- // constructors. Eventually calls to this method should be replaced
- // by something like
- // new TMatrixDSparse( some_arguments )
- TMatrixDSparse *m = new TMatrixDSparse(nr, nc);
- m->SetSparseIndex(row[nr]);
- m->SetRowIndexArray(row);
- m->SetColIndexArray(col);
- m->SetMatrixArray(data);
- return m;
+ // calculate the product v# a v
+ // where v# is the transposed vector v
+ // a: a matrix
+ // v: a vector
+ Double_t r=0.0;
+ if((a.GetNrows()!=a.GetNcols())||
+ (v.GetNrows()!=a.GetNrows())||
+ (v.GetNcols()!=1)) {
+ std::cout<<"TUnfold::MultiplyVecMSparseVec inconsistent row/col numbers "
+ <<" a("<<a.GetNrows()<<","<<a.GetNcols()
+ <<") v("<<v.GetNrows()<<","<<v.GetNcols()<<")\n";
+ }
+ const Int_t *a_rows=a.GetRowIndexArray();
+ const Int_t *a_cols=a.GetColIndexArray();
+ const Double_t *a_data=a.GetMatrixArray();
+ for (Int_t irow = 0; irow < a.GetNrows(); irow++) {
+ for(Int_t i=a_rows[irow];i<a_rows[irow+1];i++) {
+ r += v(irow,0) * a_data[i] * v(a_cols[i],0);
+ }
+ }
+ return r;
+}
+
+TMatrixDSparse *TUnfold::MultiplyMSparseMSparse(TMatrixDSparse const &a,TMatrixDSparse const &b)
+{
+ // calculate the product of two sparse matrices
+ // a,b: two sparse matrices, where a.GetNcols()==b.GetNrows()
+ // this is a replacement for the call
+ // new TMatrixDSparse(a,TMatrixDSparse::kMult,b);
+ if(a.GetNcols()!=b.GetNrows()) {
+ std::cout<<"TUnfold::MultiplyMSparseMSparse inconsistent row/col number "
+ <<a.GetNcols()<<" "<<b.GetNrows()<<"\n";
+ }
+
+ TMatrixDSparse *r=new TMatrixDSparse(a.GetNrows(),b.GetNcols());
+ const Int_t *a_rows=a.GetRowIndexArray();
+ const Int_t *a_cols=a.GetColIndexArray();
+ const Double_t *a_data=a.GetMatrixArray();
+ const Int_t *b_rows=b.GetRowIndexArray();
+ const Int_t *b_cols=b.GetColIndexArray();
+ const Double_t *b_data=b.GetMatrixArray();
+ // maximum size of the output matrix
+ int nMax=0;
+ for (Int_t irow = 0; irow < a.GetNrows(); irow++) {
+ if(a_rows[irow+1]>a_rows[irow]) nMax += b.GetNcols();
+ }
+ if(nMax>0) {
+ Int_t *r_rows=new Int_t[nMax];
+ Int_t *r_cols=new Int_t[nMax];
+ Double_t *r_data=new Double_t[nMax];
+ Double_t *row_data=new Double_t[b.GetNcols()];
+ Int_t n=0;
+ for (Int_t irow = 0; irow < a.GetNrows(); irow++) {
+ if(a_rows[irow+1]<=a_rows[irow]) continue;
+ // clear row data
+ for(Int_t icol=0;icol<b.GetNcols();icol++) {
+ row_data[icol]=0.0;
+ }
+ // loop over a-columns in this a-row
+ for(Int_t ia=a_rows[irow];ia<a_rows[irow+1];ia++) {
+ Int_t k=a_cols[ia];
+ // loop over b-columns in b-row k
+ for(Int_t ib=b_rows[k];ib<b_rows[k+1];ib++) {
+ row_data[b_cols[ib]] += a_data[ia]*b_data[ib];
+ }
+ }
+ // store nonzero elements
+ for(Int_t icol=0;icol<b.GetNcols();icol++) {
+ if(row_data[icol] != 0.0) {
+ r_rows[n]=irow;
+ r_cols[n]=icol;
+ r_data[n]=row_data[icol];
+ n++;
+ }
+ }
+ }
+ r->SetMatrixArray(n,r_rows,r_cols,r_data);
+ delete [] r_rows;
+ delete [] r_cols;
+ delete [] r_data;
+ delete [] row_data;
+ }
+
+ return r;
+}
+
+
+TMatrixDSparse *TUnfold::MultiplyMSparseTranspMSparse(TMatrixDSparse const &a,TMatrixDSparse const &b)
+{
+ // multiply a transposed Sparse matrix with another Sparse matrix
+ // a: sparse matrix (to be transposed
+ // b: sparse matrix
+ // this is a replacement for the call
+ // new TMatrixDSparse(TMatrixDSparse(TMatrixDSparse::kTransposed,a),
+ // TMatrixDSparse::kMult,b)
+ if(a.GetNrows() != b.GetNrows()) {
+ std::cout<<"TUnfold::MultiplyMSparseTranspMSparse "
+ <<"inconsistent row numbers "
+ <<a.GetNrows()<<" "<<b.GetNrows()<<"\n";
+ }
+
+ TMatrixDSparse *r=new TMatrixDSparse(a.GetNcols(),b.GetNcols());
+ const Int_t *a_rows=a.GetRowIndexArray();
+ const Int_t *a_cols=a.GetColIndexArray();
+ const Double_t *a_data=a.GetMatrixArray();
+ const Int_t *b_rows=b.GetRowIndexArray();
+ const Int_t *b_cols=b.GetColIndexArray();
+ const Double_t *b_data=b.GetMatrixArray();
+ // maximum size of the output matrix
+
+ // matrix multiplication
+ typedef std::map<Int_t,Double_t> MMatrixRow_t;
+ typedef std::map<Int_t, MMatrixRow_t > MMatrix_t;
+ MMatrix_t matrix;
+
+
+ for(Int_t iRowAB=0;iRowAB<a.GetNrows();iRowAB++) {
+ for(Int_t ia=a_rows[iRowAB];ia<a_rows[iRowAB+1];ia++) {
+ for(Int_t ib=b_rows[iRowAB];ib<b_rows[iRowAB+1];ib++) {
+ MMatrixRow_t &row=matrix[a_cols[ia]]; // creates a new row if necessary
+ MMatrixRow_t::iterator icol=row.find(b_cols[ib]);
+ if(icol!=row.end()) {
+ // update existing row
+ (*icol).second += a_data[ia]*b_data[ib];
+ } else {
+ // create new row
+ row[b_cols[ib]] = a_data[ia]*b_data[ib];
+ }
+ }
+ }
+ }
+
+ Int_t n=0;
+ for(MMatrix_t::const_iterator irow=matrix.begin();
+ irow!=matrix.end();irow++) {
+ n += (*irow).second.size();
+ }
+ if(n>0) {
+ // pack matrix into arrays
+ Int_t *r_rows=new Int_t[n];
+ Int_t *r_cols=new Int_t[n];
+ Double_t *r_data=new Double_t[n];
+ n=0;
+ for(MMatrix_t::const_iterator irow=matrix.begin();
+ irow!=matrix.end();irow++) {
+ for(MMatrixRow_t::const_iterator icol=(*irow).second.begin();
+ icol!=(*irow).second.end();icol++) {
+ r_rows[n]=(*irow).first;
+ r_cols[n]=(*icol).first;
+ r_data[n]=(*icol).second;
+ n++;
+ }
+ }
+ // pack arrays into TMatrixDSparse
+ r->SetMatrixArray(n,r_rows,r_cols,r_data);
+ delete [] r_rows;
+ delete [] r_cols;
+ delete [] r_data;
+ }
+
+ return r;
+}
+
+TMatrixDSparse *TUnfold::MultiplyMSparseM(TMatrixDSparse const &a,TMatrixD const &b)
+{
+ // multiply a Sparse matrix with a non-sparse matrix
+ // a: sparse matrix
+ // b: non-sparse matrix
+ // this is a replacement for the call
+ // new TMatrixDSparse(a,TMatrixDSparse::kMult,b);
+ if(a.GetNcols()!=b.GetNrows()) {
+ std::cout<<"TUnfold::MultiplyMSparseM inconsistent row/col number "
+ <<a.GetNcols()<<" "<<b.GetNrows()<<"\n";
+ }
+
+ TMatrixDSparse *r=new TMatrixDSparse(a.GetNrows(),b.GetNcols());
+ const Int_t *a_rows=a.GetRowIndexArray();
+ const Int_t *a_cols=a.GetColIndexArray();
+ const Double_t *a_data=a.GetMatrixArray();
+ // maximum size of the output matrix
+ int nMax=0;
+ for (Int_t irow = 0; irow < a.GetNrows(); irow++) {
+ if(a_rows[irow+1]-a_rows[irow]>0) nMax += b.GetNcols();
+ }
+ if(nMax>0) {
+ Int_t *r_rows=new Int_t[nMax];
+ Int_t *r_cols=new Int_t[nMax];
+ Double_t *r_data=new Double_t[nMax];
+
+ Int_t n=0;
+ // fill matrix r
+ for (Int_t irow = 0; irow < a.GetNrows(); irow++) {
+ if(a_rows[irow+1]-a_rows[irow]<=0) continue;
+ for(Int_t icol=0;icol<b.GetNcols();icol++) {
+ r_rows[n]=irow;
+ r_cols[n]=icol;
+ r_data[n]=0.0;
+ for(Int_t i=a_rows[irow];i<a_rows[irow+1];i++) {
+ Int_t j=a_cols[i];
+ r_data[n] += a_data[i]*b(j,icol);
+ }
+ if(r_data[n]!=0.0) n++;
+ }
+ }
+
+ r->SetMatrixArray(n,r_rows,r_cols,r_data);
+ delete [] r_rows;
+ delete [] r_cols;
+ delete [] r_data;
+ }
+ return r;
+}
+
+TMatrixD *TUnfold::InvertMSparse(TMatrixDSparse const &A) {
+ // get the inverse of a sparse matrix
+ // A: the original matrix
+ // this is a replacement of the call
+ // new TMatrixD(TMatrixD::kInverted, a);
+ // the matrix inversion is optimized for the case
+ // where a large submatrix of A is diagonal
+
+ if(A.GetNcols()!=A.GetNrows()) {
+ std::cout<<"TUnfold::InvertMSparse inconsistent row/col number "
+ <<A.GetNcols()<<" "<<A.GetNrows()<<"\n";
+ }
+
+#ifdef SCHUR_COMPLEMENT_MATRIX_INVERSION
+
+ const Int_t *a_rows=A.GetRowIndexArray();
+ const Int_t *a_cols=A.GetColIndexArray();
+ const Double_t *a_data=A.GetMatrixArray();
+
+ Int_t nmin=0,nmax=0;
+ // find largest diagonal submatrix
+ for(Int_t imin=0;imin<A.GetNrows();imin++) {
+ Int_t imax=A.GetNrows();
+ for(Int_t i2=imin;i2<imax;i2++) {
+ for(Int_t i=a_rows[i2];i<a_rows[i2+1];i++) {
+ if(a_data[i]==0.0) continue;
+ Int_t icol=a_cols[i];
+ if(icol<imin) continue; // ignore first part of the matrix
+ if(icol==i2) continue; // ignore diagonals
+ if(icol<i2) {
+ // this row spoils the matrix, so do not use
+ imax=i2;
+ break;
+ } else {
+ // this entry limits imax
+ imax=icol;
+ break;
+ }
+ }
+ }
+ if(imax-imin>nmax-nmin) {
+ nmin=imin;
+ nmax=imax;
+ }
+ }
+ // special case: complete matrix is diagonal
+ // make sure the "non-diagonal" submatrix dimension is nonzero
+ if((nmin==0)&&(nmax==A.GetNrows())) {
+ nmin++;
+ }
+ // if the diagonal part has size zero or one, use standard matrix inversion
+ if(nmin>=nmax-1) return new TMatrixD(TMatrixD::kInverted, A);
+ // A B
+ // ( )
+ // C D
+ //
+ // get inverse of diagonal part
+ std::vector<Double_t> Dinv;
+ Dinv.resize(nmax-nmin);
+ for(Int_t irow=nmin;irow<nmax;irow++) {
+ for(Int_t i=a_rows[irow];i<a_rows[irow+1];i++) {
+ if(a_cols[i]==irow) {
+ Dinv[irow-nmin]=1./a_data[i];
+ break;
+ }
+ }
+ }
+ // B*Dinv and C
+ Int_t nBDinv=0,nC=0;
+ for(Int_t irow_a=0;irow_a<A.GetNrows();irow_a++) {
+ if((irow_a<nmin)||(irow_a>=nmax)) {
+ for(Int_t i=a_rows[irow_a];i<a_rows[irow_a+1];i++) {
+ Int_t icol=a_cols[i];
+ if((icol>=nmin)&&(icol<nmax)) nBDinv++;
+ }
+ } else {
+ for(Int_t i=a_rows[irow_a];i<a_rows[irow_a+1];i++) {
+ Int_t icol=a_cols[i];
+ if((icol<nmin)||(icol>=nmax)) nC++;
+ }
+ }
+ }
+ Int_t *row_BDinv=new Int_t[nBDinv+1];
+ Int_t *col_BDinv=new Int_t[nBDinv+1];
+ Double_t *data_BDinv=new Double_t[nBDinv+1];
+
+ Int_t *row_C=new Int_t[nC+1];
+ Int_t *col_C=new Int_t[nC+1];
+ Double_t *data_C=new Double_t[nC+1];
+
+ TMatrixD Aschur(A.GetNrows()-(nmax-nmin),A.GetNcols()-(nmax-nmin));
+
+ nBDinv=0;
+ nC=0;
+ for(Int_t irow_a=0;irow_a<A.GetNrows();irow_a++) {
+ if((irow_a<nmin)||(irow_a>=nmax)) {
+ Int_t row=(irow_a<nmin) ? irow_a : (irow_a-(nmax-nmin));
+ for(Int_t i=a_rows[irow_a];i<a_rows[irow_a+1];i++) {
+ Int_t icol_a=a_cols[i];
+ if(icol_a<nmin) {
+ Aschur(row,icol_a)=a_data[i];
+ } else if(icol_a>=nmax) {
+ Aschur(row,icol_a-(nmax-nmin))=a_data[i];
+ } else {
+ row_BDinv[nBDinv]=row;
+ col_BDinv[nBDinv]=icol_a-nmin;
+ data_BDinv[nBDinv]=a_data[i]*Dinv[icol_a-nmin];
+ nBDinv++;
+ }
+ }
+ } else {
+ for(Int_t i=a_rows[irow_a];i<a_rows[irow_a+1];i++) {
+ Int_t icol_a=a_cols[i];
+ if(icol_a<nmin) {
+ row_C[nC]=irow_a-nmin;
+ col_C[nC]=icol_a;
+ data_C[nC]=a_data[i];
+ nC++;
+ } else if(icol_a>=nmax) {
+ row_C[nC]=irow_a-nmin;
+ col_C[nC]=icol_a-(nmax-nmin);
+ data_C[nC]=a_data[i];
+ nC++;
+ }
+ }
+ }
+ }
+ TMatrixDSparse BDinv(A.GetNrows()-(nmax-nmin),nmax-nmin);
+ BDinv.SetMatrixArray(nBDinv,row_BDinv,col_BDinv,data_BDinv);
+ delete[] row_BDinv;
+ delete[] col_BDinv;
+ delete[] data_BDinv;
+
+ TMatrixDSparse C(nmax-nmin,A.GetNcols()-(nmax-nmin));
+ C.SetMatrixArray(nC,row_C,col_C,data_C);
+ delete[] row_C;
+ delete[] col_C;
+ delete[] data_C;
+
+ TMatrixDSparse *BDinvC=MultiplyMSparseMSparse(BDinv,C);
+ Aschur -= *BDinvC;
+ Aschur.Invert();
+ delete BDinvC;
+
+ TMatrixD *r=new TMatrixD(A.GetNrows(),A.GetNcols());
+
+ for(Int_t row_a=0;row_a<Aschur.GetNrows();row_a++) {
+ for(Int_t col_a=0;col_a<Aschur.GetNcols();col_a++) {
+ (*r)((row_a<nmin) ? row_a : (row_a+nmax-nmin),
+ (col_a<nmin) ? col_a : (col_a+nmax-nmin))=Aschur(row_a,col_a);
+ }
+ }
+
+ TMatrixDSparse *CAschur=MultiplyMSparseM(C,Aschur);
+ TMatrixDSparse *CAschurBDinv=MultiplyMSparseMSparse(*CAschur,BDinv);
+
+ const Int_t *CAschurBDinv_row=CAschurBDinv->GetRowIndexArray();
+ const Int_t *CAschurBDinv_col=CAschurBDinv->GetColIndexArray();
+ const Double_t *CAschurBDinv_data=CAschurBDinv->GetMatrixArray();
+ for(Int_t row=0;row<CAschurBDinv->GetNrows();row++) {
+ bool has_diagonal=false;
+ for(Int_t i=CAschurBDinv_row[row];i<CAschurBDinv_row[row+1];i++) {
+ Int_t col=CAschurBDinv_col[i];
+ (*r)(row+nmin,col+nmin)=CAschurBDinv_data[i]*Dinv[row];
+ }
+ (*r)(row+nmin,row+nmin) += Dinv[row];
+ }
+
+ delete CAschurBDinv;
+
+ const Int_t *CAschur_row=CAschur->GetRowIndexArray();
+ const Int_t *CAschur_col=CAschur->GetColIndexArray();
+ const Double_t *CAschur_data=CAschur->GetMatrixArray();
+ for(Int_t row=0;row<CAschur->GetNrows();row++) {
+ for(Int_t i=CAschur_row[row];i<CAschur_row[row+1];i++) {
+ Int_t col=CAschur_col[i];
+ (*r)(row+nmin,
+ (col<nmin) ? col : (col+nmax-nmin))= -CAschur_data[i]*Dinv[row];
+ }
+ }
+ delete CAschur;
+
+ const Int_t *BDinv_row=BDinv.GetRowIndexArray();
+ const Int_t *BDinv_col=BDinv.GetColIndexArray();
+ const Double_t *BDinv_data=BDinv.GetMatrixArray();
+ for(Int_t row_aschur=0;row_aschur<Aschur.GetNrows();row_aschur++) {
+ Int_t row=(row_aschur<nmin) ? row_aschur : (row_aschur+nmax-nmin);
+ for(Int_t row_bdinv=0;row_bdinv<BDinv.GetNrows();row_bdinv++) {
+ for(Int_t i=BDinv_row[row_bdinv];i<BDinv_row[row_bdinv+1];i++) {
+ (*r)(row,BDinv_col[i]+nmin) -= Aschur(row_aschur,row_bdinv)*
+ BDinv_data[i];
+ }
+ }
+ }
+
+ return r;
+#else
+ return new TMatrixD(TMatrixD::kInverted, A);
+#endif
}
TUnfold::TUnfold(TH2 const *hist_A, EHistMap histmap, ERegMode regmode)
@@ -511,7 +921,7 @@ TUnfold::TUnfold(TH2 const *hist_A, EHistMap histmap, ERegMode regmode)
TArrayD sum_over_y(nx0);
fXToHist.Set(nx0);
fHistToX.Set(nx0);
- TArrayI colCount_A(ny);
+ Int_t nonzeroA=0;
// loop
// - calculate bias distribution
// sum_over_y
@@ -520,8 +930,9 @@ TUnfold::TUnfold(TH2 const *hist_A, EHistMap histmap, ERegMode regmode)
// - histogram bins are added to the lookup-tables
// fHistToX[] and fXToHist[]
// these convert histogram bins to matrix indices and vice versa
- // - the number of columns per row of the final matrix is counted
- // colCount_A
+ // - the number of non-zero elements is counted
+ // nonzeroA
+ Int_t nskipped=0;
for (Int_t ix = 0; ix < nx0; ix++) {
// calculate sum over all detector bins
// excluding the overflow bins
@@ -534,7 +945,7 @@ TUnfold::TUnfold(TH2 const *hist_A, EHistMap histmap, ERegMode regmode)
z = hist_A->GetBinContent(iy + 1, ix);
}
if (z > 0.0) {
- colCount_A[iy]++;
+ nonzeroA++;
sum += z;
}
}
@@ -556,27 +967,46 @@ TUnfold::TUnfold(TH2 const *hist_A, EHistMap histmap, ERegMode regmode)
}
nx++;
} else {
+ nskipped++;
// histogram bin pointing to -1 (non-existing matrix entry)
fHistToX[ix] = -1;
- std::cout << "TUnfold: bin " << ix
- << " has no influence on the data -> skipped!\n";
}
}
+ if(nskipped) {
+ std::cout << "TUnfold: the following output bins "
+ <<"are not connected to the input side\n";
+ Int_t nprint=0;
+ Int_t ixfirst=-1,ixlast=-1;
+ for (Int_t ix = 0; ix < nx0; ix++) {
+ if(fHistToX[ix]<0) {
+ nprint++;
+ if(ixlast<0) {
+ std::cout<<" "<<ix;
+ ixfirst=ix;
+ }
+ ixlast=ix;
+ }
+ if(((fHistToX[ix]>=0)&&(ixlast>=0))||
+ (nprint==nskipped)) {
+ if(ixlast>ixfirst) std::cout<<"-"<<ixlast;
+ ixfirst=-1;
+ ixlast=-1;
+ }
+ if(nprint==nskipped) {
+ std::cout<<"\n";
+ break;
+ }
+ }
+ }
// store bias as matrix
fX0 = new TMatrixD(nx, 1, sum_over_y.GetArray());
// store non-zero elements in sparse matrix fA
// construct sparse matrix fA
- // row structure
- Int_t *row_A = new Int_t[ny + 1];
- row_A[0] = 0;
- for (Int_t iy = 0; iy < ny; iy++) {
- row_A[iy + 1] = row_A[iy] + colCount_A[iy];
- }
- // column structure and data points
- Int_t *col_A= new Int_t[row_A[ny]];
- Double_t *data_A=new Double_t[row_A[ny]];
+ Int_t *rowA = new Int_t[nonzeroA];
+ Int_t *colA = new Int_t[nonzeroA];
+ Double_t *dataA = new Double_t[nonzeroA];
+ Int_t index=0;
for (Int_t iy = 0; iy < ny; iy++) {
- Int_t index = row_A[iy];
for (Int_t ix = 0; ix < nx; ix++) {
Int_t ibinx = fXToHist[ix];
Double_t z;
@@ -586,20 +1016,29 @@ TUnfold::TUnfold(TH2 const *hist_A, EHistMap histmap, ERegMode regmode)
z = hist_A->GetBinContent(iy + 1, ibinx);
}
if (z > 0.0) {
- col_A[index] = ix;
- data_A[index] = z / sum_over_y[ix];
- index++;
+ rowA[index]=iy;
+ colA[index] = ix;
+ dataA[index] = z / sum_over_y[ix];
+ index++;
}
}
}
- fA = CreateSparseMatrix(ny, nx, row_A, col_A, data_A);
- delete[] row_A;
- delete[] col_A;
- delete[] data_A;
+ fA = new TMatrixDSparse(ny,nx);
+ fA->SetMatrixArray(index,rowA,colA,dataA);
+ delete[] rowA;
+ delete[] colA;
+ delete[] dataA;
// regularisation conditions squared
fLsquared = new TMatrixDSparse(GetNx(), GetNx());
if (regmode != kRegModeNone) {
- RegularizeBins(0, 1, nx0, regmode);
+ Int_t nError=RegularizeBins(0, 1, nx0, regmode);
+ if(nError>0) {
+ if(nError>1) {
+ std::cout<<"TUnfold: "<<nError<<" regularisation conditions have been skipped\n";
+ } else {
+ std::cout<<"TUnfold: "<<nError<<" regularisation condition has been skipped\n";
+ }
+ }
}
}
@@ -640,44 +1079,45 @@ void TUnfold::SetBias(TH1 const *bias)
}
}
-void TUnfold::RegularizeSize(int bin, Double_t const &scale)
+Int_t TUnfold::RegularizeSize(int bin, Double_t const &scale)
{
// add regularisation on the size of bin i
// bin: bin number
// scale: size of the regularisation
+ // return value: number of conditions which have been skipped
// modifies data member fLsquared
Int_t i = fHistToX[bin];
if (i < 0) {
- std::cout << "TUnfold::RegularizeSize skip bin " << bin << "\n";
- return;
+ return 1;
}
(*fLsquared) (i, i) = scale * scale;
+ return 0;
}
-void TUnfold::RegularizeDerivative(int left_bin, int right_bin,
+Int_t TUnfold::RegularizeDerivative(int left_bin, int right_bin,
Double_t const &scale)
{
// add regularisation on the difference of two bins
// left_bin: 1st bin
// right_bin: 2nd bin
// scale: size of the regularisation
+ // return value: number of conditions which have been skipped
// modifies data member fLsquared
Int_t il = fHistToX[left_bin];
Int_t ir = fHistToX[right_bin];
if ((il < 0) || (ir < 0)) {
- std::cout << "TUnfold::RegularizeDerivative skip bins "
- << left_bin << "," << right_bin << "\n";
- return;
+ return 1;
}
Double_t scale_squared = scale * scale;
(*fLsquared) (il, il) += scale_squared;
(*fLsquared) (il, ir) -= scale_squared;
(*fLsquared) (ir, il) -= scale_squared;
(*fLsquared) (ir, ir) += scale_squared;
+ return 0;
}
-void TUnfold::RegularizeCurvature(int left_bin, int center_bin,
+Int_t TUnfold::RegularizeCurvature(int left_bin, int center_bin,
int right_bin,
Double_t const &scale_left,
Double_t const &scale_right)
@@ -688,6 +1128,7 @@ void TUnfold::RegularizeCurvature(int left_bin, int center_bin,
// right_bin: 3rd bin
// scale_left: scale factor on center-left difference
// scale_right: scale factor on right-center difference
+ // return value: number of conditions which have been skipped
// modifies data member fLsquared
Int_t il, ic, ir;
@@ -695,9 +1136,7 @@ void TUnfold::RegularizeCurvature(int left_bin, int center_bin,
ic = fHistToX[center_bin];
ir = fHistToX[right_bin];
if ((il < 0) || (ic < 0) || (ir < 0)) {
- std::cout << "TUnfold::RegularizeCurvature skip bins "
- << left_bin << "," << center_bin << "," << right_bin << "\n";
- return;
+ return 1;
}
Double_t r[3];
r[0] = -scale_left;
@@ -713,9 +1152,10 @@ void TUnfold::RegularizeCurvature(int left_bin, int center_bin,
(*fLsquared) (ir, il) += r[2] * r[0];
(*fLsquared) (ir, ic) += r[2] * r[1];
(*fLsquared) (ir, ir) += r[2] * r[2];
+ return 0;
}
-void TUnfold::RegularizeBins(int start, int step, int nbin,
+Int_t TUnfold::RegularizeBins(int start, int step, int nbin,
ERegMode regmode)
{
// set regulatisation on a 1-dimensional curve
@@ -723,6 +1163,8 @@ void TUnfold::RegularizeBins(int start, int step, int nbin,
// step: distance between neighbouring bins
// nbin: total number of bins
// regmode: regularisation mode
+ // return value:
+ // number of errors (i.e. conditions which have been skipped)
// modifies data member fLsquared
Int_t i0, i1, i2;
@@ -730,25 +1172,27 @@ void TUnfold::RegularizeBins(int start, int step, int nbin,
i1 = i0 + step;
i2 = i1 + step;
Int_t nSkip = 0;
+ Int_t nError= 0;
if (regmode == kRegModeDerivative)
nSkip = 1;
else if (regmode == kRegModeCurvature)
nSkip = 2;
for (Int_t i = nSkip; i < nbin; i++) {
if (regmode == kRegModeSize) {
- RegularizeSize(i0);
+ nError += RegularizeSize(i0);
} else if (regmode == kRegModeDerivative) {
- RegularizeDerivative(i0, i1);
+ nError += RegularizeDerivative(i0, i1);
} else if (regmode == kRegModeCurvature) {
- RegularizeCurvature(i0, i1, i2);
+ nError += RegularizeCurvature(i0, i1, i2);
}
i0 = i1;
i1 = i2;
i2 += step;
}
+ return nError;
}
-void TUnfold::RegularizeBins2D(int start_bin, int step1, int nbin1,
+Int_t TUnfold::RegularizeBins2D(int start_bin, int step1, int nbin1,
int step2, int nbin2, ERegMode regmode)
{
// set regularisation on a 2-dimensional grid of bins
@@ -757,14 +1201,18 @@ void TUnfold::RegularizeBins2D(int start_bin, int step1, int nbin1,
// nbin1: number of bins in 1st direction
// step2: distance between bins in 2nd direction
// nbin2: number of bins in 2nd direction
- // modifies data member fLsquared
+ // return value:
+ // number of errors (i.e. conditions which have been skipped)
+ // modifies data member fLsquared
+ Int_t nError = 0;
for (Int_t i1 = 0; i1 < nbin1; i1++) {
- RegularizeBins(start_bin + step1 * i1, step2, nbin2, regmode);
+ nError += RegularizeBins(start_bin + step1 * i1, step2, nbin2, regmode);
}
for (Int_t i2 = 0; i2 < nbin2; i2++) {
- RegularizeBins(start_bin + step2 * i2, step1, nbin1, regmode);
+ nError += RegularizeBins(start_bin + step2 * i2, step1, nbin1, regmode);
}
+ return nError;
}
Double_t TUnfold::DoUnfold(Double_t const &tau_reg, TH1 const *input,
@@ -789,35 +1237,61 @@ Double_t TUnfold::DoUnfold(Double_t const &tau_reg, TH1 const *input,
return DoUnfold(tau_reg);
}
-void TUnfold::SetInput(TH1 const *input, Double_t const &scaleBias) {
+Int_t TUnfold::SetInput(TH1 const *input, Double_t const &scaleBias,
+ Double_t oneOverZeroError) {
// Define the input data for subsequent calls to DoUnfold(Double_t)
// input: input distribution with errors
// scaleBias: scale factor applied to the bias
+ // oneOverZeroError: for bins with zero error, this number defines 1/error.
+ // Return value: number of bins with bad error
// Data members modified:
- // fY, fV, fBiasScale
+ // fY, fV, fBiasScale, fNdf
fBiasScale = scaleBias;
+
// construct inverted error matrix of measured quantities
// from errors of input histogram
- Int_t *row_V=new Int_t[GetNy() + 1];
- Int_t *col_V=new Int_t[GetNy()];
- Double_t *data_V=new Double_t[GetNy()];
- // fV is set up in increasing row/column order
- // to avoid too much memory management
+ // and count number of degrees of freedom
+ fNdf = -GetNpar();
+ Int_t *rowColV=new Int_t[GetNy()];
+ Double_t *dataV=new Double_t[GetNy()];
+ Int_t nError=0;
for (Int_t iy = 0; iy < GetNy(); iy++) {
Double_t dy = input->GetBinError(iy + 1);
- if (dy <= 0.0)
- dy = 1.0;
- row_V[iy] = iy;
- col_V[iy] = iy;
- data_V[iy] = 1 / dy / dy;
+ rowColV[iy] = iy;
+ if (dy <= 0.0) {
+ nError++;
+ dataV[iy] = oneOverZeroError*oneOverZeroError;
+ } else {
+ dataV[iy] = 1. / dy / dy;
+ }
+ if( dataV[iy]>0.0) fNdf ++;
+ }
+ if(nError>0) {
+ if(nError>1) {
+ std::cout
+ <<"TUnfold::SetInput "<<nError<<" input bins have zero error. ";
+ } else {
+ std::cout
+ <<"TUnfold::SetInput "<<nError<<" input bin has zero error. ";
+ }
+ if(oneOverZeroError !=0.0) {
+ std::cout<<"1/error is set to "<<oneOverZeroError<<"\n";
+ } else {
+ if(nError>1) {
+ std::cout
+ <<"These bins are ignored.\n";
+ } else {
+ std::cout
+ <<"This bin is ignored.\n";
+ }
+ }
}
- row_V[GetNy()] = GetNy();
if (fV)
delete fV;
- fV = CreateSparseMatrix(GetNy(), GetNy(), row_V, col_V, data_V);
- delete[] row_V;
- delete[] col_V;
- delete[] data_V;
+ fV = new TMatrixDSparse(GetNy(),GetNy());
+ fV->SetMatrixArray(GetNy(),rowColV,rowColV, dataV);
+ delete[] rowColV;
+ delete[] dataV;
//
// get input vector
if (fY)
@@ -825,7 +1299,8 @@ void TUnfold::SetInput(TH1 const *input, Double_t const &scaleBias) {
fY = new TMatrixD(GetNy(), 1);
for (Int_t i = 0; i < GetNy(); i++) {
(*fY) (i, 0) = input->GetBinContent(i + 1);
- }
+ }
+ return nError;
}
Double_t TUnfold::DoUnfold(Double_t const &tau) {
@@ -944,12 +1419,19 @@ Int_t TUnfold::ScanLcurve(Int_t nPoint,
}
// create curvature Spline
TSpline3 *splineC=new TSpline3("L curve curvature",cTi,cCi,n-1);
- Double_t cCmax=cCi[0];
- Double_t cTmax=cTi[0];
// find the maximum of the curvature
- for(Int_t i=0;i<n-2;i++) {
+ // if the parameter iskip is non-zero, then iskip points are
+ // ignored when looking for the largest curvature
+ // (there are problems with the curvature determined from the first
+ // few points of splineX,splineY in teh algorithm above)
+ Int_t iskip=0;
+ if(n>4) iskip=1;
+ if(n>7) iskip=2;
+ Double_t cCmax=cCi[iskip];
+ Double_t cTmax=cTi[iskip];
+ for(Int_t i=iskip;i<n-2-iskip;i++) {
// find maximum on this spline section
- // check boundary conditions for x[i] and x[i+1]
+ // check boundary conditions for x[i+1]
Double_t xMax=cTi[i+1];
Double_t yMax=cCi[i+1];
if(cCi[i]>yMax) {
@@ -1004,6 +1486,15 @@ Int_t TUnfold::ScanLcurve(Int_t nPoint,
cTmax=xMax;
}
}
+#ifdef DEBUG_LCURVE
+ {
+ TCanvas lcc;
+ lcc.Divide(1,1);
+ lcc.cd(1);
+ splineC->Draw();
+ lcc.SaveAs("splinec.ps");
+ }
+#endif
delete splineC;
delete[] cTi;
delete[] cCi;
@@ -1054,14 +1545,7 @@ TH1D *TUnfold::GetOutput(char const *name, char const *title,
xmax = nbins;
}
TH1D *out = new TH1D(name, title, nbins, xmin, xmax);
-#ifdef UNUSED
- for (Int_t i = 0; i < GetNx(); i++) {
- out->SetBinContent(fXToHist[i], (*fX) (i, 0));
- out->SetBinError(fXToHist[i], TMath::Sqrt((*fE) (i, i)));
- }
-#else
GetOutput(out);
-#endif
return out;
}
@@ -1132,11 +1616,13 @@ TH1D *TUnfold::GetInput(char const *name, char const *title,
y1 = GetNy();
}
TH1D *out = new TH1D(name, title, GetNy(), y0, y1);
- TMatrixD Vinv(TMatrixD::kInverted, *fV);
+
+ TMatrixD *Vinv=InvertMSparse(*fV);
for (Int_t i = 0; i < GetNy(); i++) {
out->SetBinContent(i + 1, (*fY) (i, 0));
- out->SetBinError(i + 1, TMath::Sqrt(Vinv(i, i)));
+ out->SetBinError(i + 1, TMath::Sqrt((*Vinv)(i, i)));
}
+ delete Vinv;
return out;
}
@@ -1156,19 +1642,7 @@ TH2D *TUnfold::GetRhoIJ(char const *name, char const *title,
xmax = nbins;
}
TH2D *out = new TH2D(name, title, nbins, xmin, xmax, nbins, xmin, xmax);
-#ifdef UNUSED
- for (Int_t i = 0; i < GetNx(); i++) {
- for (Int_t j = 0; j < GetNx(); j++) {
- out->SetBinContent(fXToHist[i], fXToHist[j],
- (*fE) (i,
- j) / TMath::Sqrt((*fE) (i,
- i) * (*fE) (j,
- j)));
- }
- }
-#else
GetRhoIJ(out);
-#endif
return out;
}
@@ -1187,15 +1661,7 @@ TH2D *TUnfold::GetEmatrix(char const *name, char const *title,
xmax = nbins;
}
TH2D *out = new TH2D(name, title, nbins, xmin, xmax, nbins, xmin, xmax);
-#ifdef UNUSED
- for (Int_t i = 0; i < GetNx(); i++) {
- for (Int_t j = 0; j < GetNx(); j++) {
- out->SetBinContent(fXToHist[i], fXToHist[j], (*fE) (i, j));
- }
- }
-#else
GetEmatrix(out);
-#endif
return out;
}
@@ -1295,9 +1761,21 @@ Double_t const &TUnfold::GetChi2L(void) const
return fChi2L;
}
+Int_t TUnfold::GetNdf(void) const
+{
+ // return number of degrees of freedom
+ return fNdf;
+}
+
+Int_t TUnfold::GetNpar(void) const
+{
+ // return number of parameters
+ return GetNx();
+}
+
Double_t TUnfold::GetLcurveX(void) const {
// return value on x axis of L curve
- return TMath::Log10(fChi2A+fChi2L);
+ return TMath::Log10(fChi2A);
}
Double_t TUnfold::GetLcurveY(void) const {
diff --git a/tutorials/math/testUnfold1.C b/tutorials/math/testUnfold1.C
index 1032dde2d1c07d2e0871a5070a39c6311e668ddb..17ac7fdc46d04db05811e5ad711351d01cdb5c74 100644
--- a/tutorials/math/testUnfold1.C
+++ b/tutorials/math/testUnfold1.C
@@ -1,10 +1,14 @@
-// Test program for the class TUnfold
// Author: Stefan Schmitt
// DESY, 14.10.2008
-// Version 6a, fix problem with dynamic array allocation under windows
+// Version 11, print chi**2 and number of degrees of freedom
//
// History:
+// Version 10, with bug-fix in TUnfold.cxx
+// Version 9, with bug-fix in TUnfold.cxx and TUnfold.h
+// Version 8, with bug-fix in TUnfold.cxx and TUnfold.h
+// Version 7, with bug-fix in TUnfold.cxx and TUnfold.h
+// Version 6a, fix problem with dynamic array allocation under windows
// Version 6, bug-fixes in TUnfold.C
// Version 5, replace main() by testUnfold1()
// Version 4, with bug-fix in TUnfold.C
@@ -253,11 +257,12 @@ int testUnfold1()
Int_t iBest;
TSpline *logTauX,*logTauY;
TGraph *lCurve;
-
// this method scans the parameter tau and finds the kink in the L curve
// finally, the unfolding is done for the best choice of tau
iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
std::cout<<"tau="<<unfold.GetTau()<<"\n";
+ std::cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
+ <<" / "<<unfold.GetNdf()<<"\n";
// save graphs with one point to visualize best choice of tau
Double_t t[1],x[1],y[1];
@@ -302,7 +307,6 @@ int testUnfold1()
// fit Breit-Wigner shape to unfolded data, using the full error matrix
// here we use a "user" chi**2 function to take into account
// the full covariance matrix
- //TVirtualFitter::SetDefaultFitter("TMinuit");
gFitter=TVirtualFitter::Fitter(histMunfold);
gFitter->SetFCN(chisquare_corr);
@@ -372,6 +376,7 @@ int testUnfold1()
bestLcurve->Draw("*");
output.SaveAs("c1.ps");
+
return 0;
}
diff --git a/tutorials/math/testUnfold2.C b/tutorials/math/testUnfold2.C
index 61c69dec7ffef7c0914b88d8325d09104ce7ca85..7477c93a95007a05f0afc350019877ac596f4188 100644
--- a/tutorials/math/testUnfold2.C
+++ b/tutorials/math/testUnfold2.C
@@ -1,10 +1,14 @@
-// Test program for the class MyUnfold, derived from TUnfold
// Author: Stefan Schmitt
// DESY, 14.10.2008
-// Version 6a, fix problem with dynamic array allocation under windows
+// Version 11, print chi**2 and number of degrees of freedom
//
// History:
+// Version 10, with bug-fix in TUnfold.cxx
+// Version 9, with bug-fix in TUnfold.cxx, TUnfold.h
+// Version 8, with bug-fix in TUnfold.cxx, TUnfold.h
+// Version 7, with bug-fix in TUnfold.cxx, TUnfold.h
+// Version 6a, fix problem with dynamic array allocation under windows
// Version 6, re-include class MyUnfold in the example
// Version 5, move class MyUnfold to seperate files
// Version 4, with bug-fix in TUnfold.C
@@ -67,7 +71,7 @@ protected:
Int_t const *fBinMap; // bin mapping to extract the global correlation
Double_t fTauBest; // tau with the smallest correlation
Double_t fRhoMin; // smallest correlation
- ClassDef(MyUnfold,0);
+ //ClassDef(MyUnfold,0);
};
//ClassImp(MyUnfold)
@@ -350,6 +354,8 @@ int testUnfold2()
// finally, the unfolding is done for the best choice of tau
iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
std::cout<<"tau="<<unfold.GetTau()<<"\n";
+ std::cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
+ <<" / "<<unfold.GetNdf()<<"\n";
Double_t t[1],x[1],y[1];
logTauX->GetKnot(iBest,t[0],x[0]);
logTauY->GetKnot(iBest,t[0],y[0]);