ROOT » HIST » HIST » TUnfold

class TUnfold: public TObject


  TUnfold is used to decompose a measurement y into several sources x
  given the measurement uncertainties and a matrix of migrations A


  * For most applications, it is better to use TUnfoldDensity *
  * instead of using TUnfoldSys or TUnfold                    *


  If you use this software, please consider the following citation
       S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]

  More documentation and updates are available on
      http://www.desy.de/~sschmitt


  A short summary of the algorithm is given in the following:

  the "best" x matching the measurement y within errors
  is determined by minimizing the function
    L1+L2+L3

  where
   L1 = (y-Ax)# Vyy^-1 (y-Ax)
   L2 = tau^2 (L(x-x0))# L(x-x0)
   L3 = lambda sum_i(y_i -(Ax)_i)

  [the notation # means that the matrix is transposed ]

  The term L1 is familiar from a least-square minimisation
  The term L2 defines the regularisation (smootheness condition on x),
    where the parameter tau^2 gives the strength of teh regularisation
  The term L3 is an optional area constraint with Lagrangian parameter
    lambda, ensuring that that the normalisation of x is consistent with the
    normalisation of y

  The method can be applied to a very large number of problems,
  where the measured distribution y is a linear superposition
  of several Monte Carlo shapes

 Input from measurement:

   y: vector of measured quantities  (dimension ny)
   Vyy: covariance matrix for y (dimension ny x ny)
        in many cases V is diagonal and calculated from the errors of y

 From simulation:

   A: migration matrix               (dimension ny x nx)

 Result

   x: unknown underlying distribution (dimension nx)
      The error matrix of x, V_xx, is also determined

 Regularisation

   tau: parameter, defining the regularisation strength
   L: matrix of regularisation conditions (dimension nl x nx)
        depends on the structure of the input data
   x0: bias distribution, from simulation

 Preservation of the area

   lambda: lagrangian multiplier
   y_i: one component of the vector y
   (Ax)_i: one component of the vector Ax


 Determination of the unfolding result x:

   (a) not constrained: minimisation is performed as a function of x
         for fixed lambda=0
  or
   (b) constrained: stationary point is found as a function of x and lambda

 The constraint can be useful to reduce biases on the result x
 in cases where the vector y follows non-Gaussian probability densities
 (example: Poisson statistics at counting experiments in particle physics)

 Some random examples:

   (1) measure a cross-section as a function of, say, E_T(detector)
        and unfold it to obtain the underlying distribution E_T(generator)
   (2) measure a lifetime distribution and unfold the contributions from
        different flavours
   (3) measure the transverse mass and decay angle
        and unfold for the true mass distribution plus background

 Documentation

 Some technical documentation is available here:
   http://www.desy.de/~sschmitt

 Note:
   For most applications it is better to use the derived class
     TUnfoldSys or even better TUnfoldDensity

  TUnfoldSys extends the functionality of TUnfold
             such that systematic errors are propagated to teh result
             and that the unfolding can be done with proper background
             subtraction

  TUnfoldDensity extends further the functionality of TUnfoldSys
                 complex binning schemes are supported
                 The binning of input histograms is handeld consistently:
                  (1) the regularisation may be done by density,
                      i.e respecting the bin widths
                  (2) methods are provided which preserve the proper binning
                      of the result histograms

 Implementation

 The result of the unfolding is calculated as follows:

    Lsquared = L#L            regularisation conditions squared

    epsilon_j = sum_i A_ij    vector of efficiencies

    E^-1  = ((A# Vyy^-1 A)+tau^2 Lsquared)

    x = E (A# Vyy^-1 y + tau^2 Lsquared x0 +lambda/2 * epsilon) is the result

 The derivatives
    dx_k/dy_i
    dx_k/dA_ij
    dx_k/d(tau^2)
 are calculated for further usage.

 The covariance matrix V_xx is calculated as:
    Vxx_ij = sum_kl dx_i/dy_k Vyy_kl dx_j/dy_l

 Warning:

  The algorithm is based on "standard" matrix inversion, with the
  known limitations in numerical accuracy and computing cost for
  matrices with large dimensions.

  Thus the algorithm should not used for large dimensions of x and y
    nx should not be much larger than 200
    ny should not be much larger than 1000


 Proper choice of tau

 One of the difficult questions is about the choice of tau.
 The method implemented in TUnfold is the L-curve method:
 a two-dimensional curve is plotted
   x-axis: log10(chisquare)
   y-axis: log10(regularisation condition)
 In many cases this curve has an L-shape. The best choice of tau is in the
 kink of the L

 Within TUnfold a simple version of the L-curve analysis is available.
 It tests a given number of points in a predefined tau-range and searches
 for the maximum of the curvature in the L-curve (kink position).
 if no tau range is given, the range of the scan is determined automatically

  A nice overview of the L-curve 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

 Alternative Regularisation conditions

 Regularisation is needed for most unfolding problems, in order to avoid
 large oscillations and large correlations on the output bins.
 It means that some extra conditions are applied on the output bins

 Within TUnfold these conditions are posed on the difference (x-x0), where
    x:  unfolding output
    x0: the bias distribution, by default calculated from
        the input matrix A. There is a method SetBias() to change the
        bias distribution.
        The 3rd argument to DoUnfold() is a scale factor applied to the bias
          bias_default[j] = sum_i A[i][j]
          x0[j] = scaleBias*bias[j]
        The scale factor can be used to
         (a) completely suppress the bias by setting it to zero
         (b) compensate differences in the normalisation between data
             and Monte Carlo

 If the regularisation is strong, i.e. large parameter tau,
 then the distribution x or its derivatives will look like the bias
 distribution. If the parameter tau is small, the distribution x is
 independent of the bias.

 Three basic types of regularisation are implemented in TUnfold

    condition            regularisation

    kRegModeNone         none
    kRegModeSize         minimize the size of (x-x0)
    kRegModeDerivative   minimize the 1st derivative of (x-x0)
    kRegModeCurvature    minimize the 2nd derivative of (x-x0)

 kRegModeSize is the regularisation scheme which often is found in
 literature. The second derivative is often named curvature.
 Sometimes the bias is not discussed,  equivalent to a bias scale factor
 of zero.

 The regularisation schemes kRegModeDerivative and
 kRegModeCurvature have the nice feature that they create correlations
 between x-bins, whereas the non-regularized unfolding tends to create
 negative correlations between bins. For these regularisation schemes the
 parameter tau could be tuned such that the correlations are smallest,
 as an alternative to the L-curve method.

 If kRegModeSize is chosen or if x is a smooth function through all bins,
 the regularisation condition can be set on all bins together by giving
 the appropriate argument in the constructor (see examples above).

 If x is composed of independent groups of bins (for example,
 signal and background binning in two variables), it may be necessary to
 set regularisation conditions for the individual groups of bins.
 In this case,  give  kRegModeNone  in the constructor and specify
 the bin grouping with calls to
          RegularizeBins()   specify a 1-dimensional group of bins
          RegularizeBins2D() specify a 2-dimensional group of bins

 Note, the class TUnfoldDensity provides an automatic setup of complex
 regularisation schemes

 For ultimate flexibility, the regularisation condition can be set on each
 bin individually
  -> give  kRegModeNone  in the constructor and use
      RegularizeSize()        regularize one bin
      RegularizeDerivative()  regularize the slope given by two bins
      RegularizeCurvature()   regularize the curvature given by three bins
      AddRegularisationCondition()
                              define an arbitrary regulatisation condition


Function Members (Methods)

public:
virtual~TUnfold()
voidTObject::AbstractMethod(const char* method) const
virtual voidTObject::AppendPad(Option_t* option = "")
virtual voidTObject::Browse(TBrowser* b)
static TClass*Class()
virtual const char*TObject::ClassName() const
virtual voidTObject::Clear(Option_t* = "")
virtual TObject*TObject::Clone(const char* newname = "") const
virtual Int_tTObject::Compare(const TObject* obj) const
virtual voidTObject::Copy(TObject& object) const
virtual voidTObject::Delete(Option_t* option = "")MENU
virtual Int_tTObject::DistancetoPrimitive(Int_t px, Int_t py)
virtual Double_tDoUnfold(Double_t tau)
Double_tDoUnfold(Double_t tau, const TH1* hist_y, Double_t scaleBias = 0.)
virtual voidTObject::Draw(Option_t* option = "")
virtual voidTObject::DrawClass() constMENU
virtual TObject*TObject::DrawClone(Option_t* option = "") constMENU
virtual voidTObject::Dump() constMENU
virtual voidTObject::Error(const char* method, const char* msgfmt) const
virtual voidTObject::Execute(const char* method, const char* params, Int_t* error = 0)
virtual voidTObject::Execute(TMethod* method, TObjArray* params, Int_t* error = 0)
virtual voidTObject::ExecuteEvent(Int_t event, Int_t px, Int_t py)
virtual voidTObject::Fatal(const char* method, const char* msgfmt) const
virtual TObject*TObject::FindObject(const char* name) const
virtual TObject*TObject::FindObject(const TObject* obj) const
voidGetBias(TH1* bias, const Int_t* binMap = 0) const
Double_tGetChi2A() const
Double_tGetChi2L() const
virtual Option_t*TObject::GetDrawOption() const
static Long_tTObject::GetDtorOnly()
voidGetEmatrix(TH2* ematrix, const Int_t* binMap = 0) const
Double_tGetEpsMatrix() const
voidGetFoldedOutput(TH1* folded, const Int_t* binMap = 0) const
virtual const char*TObject::GetIconName() const
voidGetInput(TH1* inputData, const Int_t* binMap = 0) const
voidGetInputInverseEmatrix(TH2* ematrix)
voidGetL(TH2* l) const
virtual Double_tGetLcurveX() const
virtual Double_tGetLcurveY() const
voidGetLsquared(TH2* lsquared) const
virtual const char*TObject::GetName() const
Int_tGetNdf() const
voidGetNormalisationVector(TH1* s, const Int_t* binMap = 0) const
Int_tGetNpar() const
Int_tGetNr() const
virtual char*TObject::GetObjectInfo(Int_t px, Int_t py) const
static Bool_tTObject::GetObjectStat()
virtual Option_t*TObject::GetOption() const
voidGetOutput(TH1* output, const Int_t* binMap = 0) const
voidGetProbabilityMatrix(TH2* A, TUnfold::EHistMap histmap) const
Double_tGetRhoAvg() const
Double_tGetRhoI(TH1* rhoi, const Int_t* binMap = 0, TH2* invEmat = 0) const
voidGetRhoIJ(TH2* rhoij, const Int_t* binMap = 0) const
Double_tGetRhoMax() const
Double_tGetTau() const
virtual const char*TObject::GetTitle() const
static const char*GetTUnfoldVersion()
virtual UInt_tTObject::GetUniqueID() const
virtual Bool_tTObject::HandleTimer(TTimer* timer)
virtual ULong_tTObject::Hash() const
virtual voidTObject::Info(const char* method, const char* msgfmt) const
virtual Bool_tTObject::InheritsFrom(const char* classname) const
virtual Bool_tTObject::InheritsFrom(const TClass* cl) const
virtual voidTObject::Inspect() constMENU
voidTObject::InvertBit(UInt_t f)
virtual TClass*IsA() const
virtual Bool_tTObject::IsEqual(const TObject* obj) const
virtual Bool_tTObject::IsFolder() const
Bool_tTObject::IsOnHeap() const
virtual Bool_tTObject::IsSortable() const
Bool_tTObject::IsZombie() const
virtual voidTObject::ls(Option_t* option = "") const
voidTObject::MayNotUse(const char* method) const
virtual Bool_tTObject::Notify()
voidTObject::Obsolete(const char* method, const char* asOfVers, const char* removedFromVers) const
voidTObject::operator delete(void* ptr)
voidTObject::operator delete(void* ptr, void* vp)
voidTObject::operator delete[](void* ptr)
voidTObject::operator delete[](void* ptr, void* vp)
void*TObject::operator new(size_t sz)
void*TObject::operator new(size_t sz, void* vp)
void*TObject::operator new[](size_t sz)
void*TObject::operator new[](size_t sz, void* vp)
TUnfold&operator=(const TUnfold&)
virtual voidTObject::Paint(Option_t* option = "")
virtual voidTObject::Pop()
virtual voidTObject::Print(Option_t* option = "") const
virtual Int_tTObject::Read(const char* name)
virtual voidTObject::RecursiveRemove(TObject* obj)
Int_tRegularizeBins(int start, int step, int nbin, TUnfold::ERegMode regmode)
Int_tRegularizeBins2D(int start_bin, int step1, int nbin1, int step2, int nbin2, TUnfold::ERegMode regmode)
Int_tRegularizeCurvature(int left_bin, int center_bin, int right_bin, Double_t scale_left = 1., Double_t scale_right = 1.)
Int_tRegularizeDerivative(int left_bin, int right_bin, Double_t scale = 1.)
Int_tRegularizeSize(int bin, Double_t scale = 1.)
voidTObject::ResetBit(UInt_t f)
virtual voidTObject::SaveAs(const char* filename = "", Option_t* option = "") constMENU
virtual voidTObject::SavePrimitive(ostream& out, Option_t* option = "")
virtual Int_tScanLcurve(Int_t nPoint, Double_t tauMin, Double_t tauMax, TGraph** lCurve, TSpline** logTauX = 0, TSpline** logTauY = 0)
voidSetBias(const TH1* bias)
voidTObject::SetBit(UInt_t f)
voidTObject::SetBit(UInt_t f, Bool_t set)
voidSetConstraint(TUnfold::EConstraint constraint)
virtual voidTObject::SetDrawOption(Option_t* option = "")MENU
static voidTObject::SetDtorOnly(void* obj)
voidSetEpsMatrix(Double_t eps)
virtual Int_tSetInput(const TH1* hist_y, Double_t scaleBias = 0., Double_t oneOverZeroError = 0., const TH2* hist_vyy = 0, const TH2* hist_vyy_inv = 0)
static voidTObject::SetObjectStat(Bool_t stat)
virtual voidTObject::SetUniqueID(UInt_t uid)
virtual voidShowMembers(TMemberInspector& insp) const
virtual voidStreamer(TBuffer&)
voidStreamerNVirtual(TBuffer& ClassDef_StreamerNVirtual_b)
virtual voidTObject::SysError(const char* method, const char* msgfmt) const
Bool_tTObject::TestBit(UInt_t f) const
Int_tTObject::TestBits(UInt_t f) const
TUnfold(const TUnfold&)
TUnfold(const TH2* hist_A, TUnfold::EHistMap histmap, TUnfold::ERegMode regmode = kRegModeSize, TUnfold::EConstraint constraint = kEConstraintArea)
virtual voidTObject::UseCurrentStyle()
virtual voidTObject::Warning(const char* method, const char* msgfmt) const
virtual Int_tTObject::Write(const char* name = 0, Int_t option = 0, Int_t bufsize = 0)
virtual Int_tTObject::Write(const char* name = 0, Int_t option = 0, Int_t bufsize = 0) const
protected:
voidAddMSparse(TMatrixDSparse* dest, Double_t f, const TMatrixDSparse* src) const
Bool_tAddRegularisationCondition(Int_t nEle, const Int_t* indices, const Double_t* rowData)
Bool_tAddRegularisationCondition(Int_t i0, Double_t f0, Int_t i1 = -1, Double_t f1 = 0., Int_t i2 = -1, Double_t f2 = 0.)
voidClearHistogram(TH1* h, Double_t x = 0.) const
virtual voidClearResults()
TMatrixDSparse*CreateSparseMatrix(Int_t nrow, Int_t ncol, Int_t nele, Int_t* row, Int_t* col, Double_t* data) const
static voidDeleteMatrix(TMatrixD** m)
static voidDeleteMatrix(TMatrixDSparse** m)
virtual voidTObject::DoError(int level, const char* location, const char* fmt, va_list va) const
virtual Double_tDoUnfold()
voidErrorMatrixToHist(TH2* ematrix, const TMatrixDSparse* emat, const Int_t* binMap, Bool_t doClear) const
const TMatrixDSparse*GetAx() const
Int_tGetBinFromRow(int ix) const
const TMatrixDSparse*GetDXDAM(int i) const
const TMatrixDSparse*GetDXDAZ(int i) const
const TMatrixDSparse*GetDXDtauSquared() const
const TMatrixDSparse*GetDXDY() const
const TMatrixDSparse*GetE() const
const TMatrixDSparse*GetEinv() const
Int_tGetNx() const
Int_tGetNy() const
virtual TStringGetOutputBinName(Int_t iBinX) const
Double_tGetRhoIFromMatrix(TH1* rhoi, const TMatrixDSparse* eOrig, const Int_t* binMap, TH2* invEmat) const
Int_tGetRowFromBin(int ix) const
const TMatrixDSparse*GetVxx() const
const TMatrixDSparse*GetVxxInv() const
const TMatrixDSparse*GetVyyInv() const
const TMatrixD*GetX() const
TMatrixDSparse*InvertMSparseSymmPos(const TMatrixDSparse* A, Int_t* rank) const
voidTObject::MakeZombie()
TMatrixDSparse*MultiplyMSparseM(const TMatrixDSparse* a, const TMatrixD* b) const
TMatrixDSparse*MultiplyMSparseMSparse(const TMatrixDSparse* a, const TMatrixDSparse* b) const
TMatrixDSparse*MultiplyMSparseMSparseTranspVector(const TMatrixDSparse* m1, const TMatrixDSparse* m2, const TMatrixTBase<Double_t>* v) const
TMatrixDSparse*MultiplyMSparseTranspMSparse(const TMatrixDSparse* a, const TMatrixDSparse* b) const
TUnfold()
private:
voidInitTUnfold()

Data Members

public:
static TObject::<anonymous>TObject::kBitMask
static TObject::EStatusBitsTObject::kCanDelete
static TObject::EStatusBitsTObject::kCannotPick
static TUnfold::EConstraintkEConstraintArea
static TUnfold::EConstraintkEConstraintNone
static TObject::EStatusBitsTObject::kHasUUID
static TUnfold::EHistMapkHistMapOutputHoriz
static TUnfold::EHistMapkHistMapOutputVert
static TObject::EStatusBitsTObject::kInvalidObject
static TObject::<anonymous>TObject::kIsOnHeap
static TObject::EStatusBitsTObject::kIsReferenced
static TObject::EStatusBitsTObject::kMustCleanup
static TObject::EStatusBitsTObject::kNoContextMenu
static TObject::<anonymous>TObject::kNotDeleted
static TObject::EStatusBitsTObject::kObjInCanvas
static TObject::<anonymous>TObject::kOverwrite
static TUnfold::ERegModekRegModeCurvature
static TUnfold::ERegModekRegModeDerivative
static TUnfold::ERegModekRegModeMixed
static TUnfold::ERegModekRegModeNone
static TUnfold::ERegModekRegModeSize
static TObject::<anonymous>TObject::kSingleKey
static TObject::<anonymous>TObject::kWriteDelete
static TObject::<anonymous>TObject::kZombie
protected:
TMatrixDSparse*fAInput: matrix
Double_tfBiasScaleInput: scale factor for the bias
TUnfold::EConstraintfConstraintInput: type of constraint to use
TArrayIfHistToXInput: histogram bins -> matrix indices
TMatrixDSparse*fLInput: regularisation conditions
TUnfold::ERegModefRegModeInput: type of regularisation
TArrayDfSumOverYInput: sum of all columns
Double_tfTauSquaredInput: regularisation parameter
TMatrixDSparse*fVyyInput: covariance matrix for y
TMatrixD*fX0Input: x0
TArrayIfXToHistInput: matrix indices -> histogram bins
TMatrixD*fYInput: y
private:
TMatrixDSparse*fAxResult: Ax
Double_tfChi2AResult: chi**2 contribution from (y-Ax)V(y-Ax)
TMatrixDSparse*fDXDAM[2]Result: part of derivative dx_k/dA_ij
TMatrixDSparse*fDXDAZ[2]Result: part of derivative dx_k/dA_ij
TMatrixDSparse*fDXDYResult: derivative dx/dy
TMatrixDSparse*fDXDtauSquaredResult: derivative dx/dtau
TMatrixDSparse*fEResult: matrix E
TMatrixDSparse*fEinvResult: matrix E^(-1)
Double_tfEpsMatrixmachine accuracy for eingenvalue analysis
Int_tfIgnoredBinsnumber of input bins which are dropped because they have error=0
Double_tfLXsquaredResult: chi**2 contribution from (x-s*x0)Lsquared(x-s*x0)
Int_tfNdfResult: number of degrees of freedom
Double_tfRhoAvgResult: average global correlation
Double_tfRhoMaxResult: maximum global correlation
TMatrixDSparse*fVxxResult: covariance matrix on x
TMatrixDSparse*fVxxInvResult: inverse of covariance matrix on x
TMatrixDSparse*fVyyInvResult: inverse of covariance matrix on y
TMatrixD*fXResult: x

Class Charts

Inheritance Inherited Members Includes Libraries
Class Charts

Function documentation

const char * GetTUnfoldVersion(void)
void InitTUnfold(void)
 reset all data members
void DeleteMatrix(TMatrixD** m)
void DeleteMatrix(TMatrixDSparse** m)
void ClearResults(void)
 delete old results (if any)
 this function is virtual, so derived classes may implement their own
 method to flag results as non-valid
TUnfold(const TUnfold& )
 set all matrix pointers to zero
Double_t DoUnfold(Double_t tau)
 main unfolding algorithm. Declared virtual, because other algorithms
 could be implemented

 Purpose: unfold y -> x
 Data members required:
     fA:  matrix to relate x and y
     fY:  measured data points
     fX0: bias on x
     fBiasScale: scale factor for fX0
     fVyy:  covariance matrix for y
     fL: regularisation conditions
     fTauSquared: regularisation strength
     fConstraint: whether the constraint is applied
 Data members modified:
     fVyyInv: inverse of input data covariance matrix
     fNdf: number of dgerres of freedom
     fEinv: inverse of the matrix needed for unfolding calculations
     fE:    the matrix needed for unfolding calculations
     fX:    unfolded data points
     fDXDY: derivative of x wrt y (for error propagation)
     fVxx:  error matrix (covariance matrix) on x
     fAx:   estimate of distribution y from unfolded data
     fChi2A:  contribution to chi**2 from y-Ax
     fChi2L:  contribution to chi**2 from L*(x-x0)
     fDXDtauSquared: derivative of x wrt tau
     fDXDAM[0,1]: matrix parts of derivative x wrt A
     fDXDAZ[0,1]: vector parts of derivative x wrt A
     fRhoMax: maximum global correlation coefficient
     fRhoAvg: average global correlation coefficient
 return code:
     fRhoMax   if(fRhoMax>=1.0) then the unfolding has failed!
TMatrixDSparse * MultiplyMSparseMSparse(const TMatrixDSparse* a, const TMatrixDSparse* b) const
 calculate the product of two sparse matrices
    a,b: pointers to sparse matrices, where a->GetNcols()==b->GetNrows()
 this is a replacement for the call
    new TMatrixDSparse(*a,TMatrixDSparse::kMult,*b);
TMatrixDSparse * MultiplyMSparseM(const TMatrixDSparse* a, const TMatrixD* b) const
 multiply a Sparse matrix with a non-sparse matrix
    a:  pointer to sparse matrix
    b:  pointer to non-sparse matrix
 this is a replacement for the call
    new TMatrixDSparse(*a,TMatrixDSparse::kMult,*b);
void AddMSparse(TMatrixDSparse* dest, Double_t f, const TMatrixDSparse* src) const
 a replacement for
     (*dest) += f*(*src)
TString GetOutputBinName(Int_t iBinX) const
 given a bin number, return the name of the output bin
 this method makes more sense for the class TUnfoldDnesity
 where it gets overwritten
TUnfold(const TH2* hist_A, TUnfold::EHistMap histmap, TUnfold::ERegMode regmode = kRegModeSize, TUnfold::EConstraint constraint = kEConstraintArea)
 set up unfolding matrix and initial regularisation scheme
    hist_A:  matrix that describes the migrations
    histmap: mapping of the histogram axes to the unfolding output
    regmode: global regularisation mode
    constraint: type of constraint to use
 data members initialized to something different from zero:
    fA: filled from hist_A
    fDA: filled from hist_A
    fX0: filled from hist_A
    fL: filled depending on the regularisation scheme
 Treatment of overflow bins
    Bins where the unfolding input (Detector level) is in overflow
    are used for the efficiency correction. They have to be filled
    properly!
    Bins where the unfolding output (Generator level) is in overflow
    are treated as a part of the generator level distribution.
    I.e. the unfolding output could have non-zero overflow bins if the
    input matrix does have such bins.
~TUnfold(void)
 delete all data members
void SetBias(const TH1* bias)
 initialize alternative bias from histogram
 modifies data member fX0
Int_t RegularizeSize(int bin, Double_t scale = 1.)
 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 fL
Int_t RegularizeDerivative(int left_bin, int right_bin, Double_t scale = 1.)
 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 fL
Int_t RegularizeCurvature(int left_bin, int center_bin, int right_bin, Double_t scale_left = 1., Double_t scale_right = 1.)
 add regularisation on the curvature through 3 bins (2nd derivative)
   left_bin: 1st bin
   center_bin: 2nd 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 fL
Int_t RegularizeBins(int start, int step, int nbin, TUnfold::ERegMode regmode)
 set regulatisation on a 1-dimensional curve
   start: first bin
   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 fL
Int_t RegularizeBins2D(int start_bin, int step1, int nbin1, int step2, int nbin2, TUnfold::ERegMode regmode)
 set regularisation on a 2-dimensional grid of bins
     start: first bin
     step1: distance between bins in 1st direction
     nbin1: number of bins in 1st direction
     step2: distance between bins in 2nd direction
     nbin2: number of bins in 2nd direction
 return value:
   number of errors (i.e. conditions which have been skipped)
 modifies data member fL
Double_t DoUnfold(Double_t tau, const TH1* hist_y, Double_t scaleBias = 0.)
 Do unfolding of an input histogram
   tau_reg: regularisation parameter
   input:   input distribution with errors
   scaleBias:  scale factor applied to the bias
 Data members required:
   fA, fX0, fL
 Data members modified:
   those documented in SetInput()
   and those documented in DoUnfold(Double_t)
 Return value:
   maximum global correlation coefficient
   NOTE!!! return value >=1.0 means error, and the result is junk

 Overflow bins of the input distribution are ignored!
Int_t SetInput(const TH1* hist_y, Double_t scaleBias = 0., Double_t oneOverZeroError = 0., const TH2* hist_vyy = 0, const TH2* hist_vyy_inv = 0)
 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.
   hist_vyy: if non-zero, defines the data covariance matrix
             otherwise it is calculated from the data errors
   hist_vyy_inv: if non-zero and if hist_vyy  is set, defines the inverse of the data covariance matrix
 Return value: number of bins with bad error
                 +10000*number of unconstrained output bins
         Note: return values>=10000 are fatal errors,
               for the given input, the unfolding can not be done!
 Data members modified:
   fY, fVyy, , fBiasScale
 Data members cleared
   fVyyInv, fNdf
   + see ClearResults
Double_t DoUnfold(Double_t tau)
 Unfold with given value of regularisation parameter tau
     tau: new tau parameter
 required data members:
     fA:  matrix to relate x and y
     fY:  measured data points
     fX0: bias on x
     fBiasScale: scale factor for fX0
     fV:  inverse of covariance matrix for y
     fL: regularisation conditions
 modified data members:
     fTauSquared and those documented in DoUnfold(void)
Int_t ScanLcurve(Int_t nPoint, Double_t tauMin, Double_t tauMax, TGraph** lCurve, TSpline** logTauX = 0, TSpline** logTauY = 0)
 scan the L curve
   nPoint: number of points on the resulting curve
   tauMin: smallest tau value to study
   tauMax: largest tau value to study
   lCurve: the L curve as graph
   logTauX: output spline of x-coordinates vs tau for the L curve
   logTauY: output spline of y-coordinates vs tau for the L curve
 return value: the coordinate number (0..nPoint-1) with the "best" choice
     of tau
void GetNormalisationVector(TH1* s, const Int_t* binMap = 0) const
 get vector of normalisation factors
   out: output histogram
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

void GetBias(TH1* bias, const Int_t* binMap = 0) const
 get bias distribution, possibly with bin remapping
   out: output histogram
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

void GetFoldedOutput(TH1* folded, const Int_t* binMap = 0) const
 get unfolding result, folded back trough the matrix
   out: output histogram
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

void GetProbabilityMatrix(TH2* A, TUnfold::EHistMap histmap) const
 retreive matrix of probabilities
    histmap: on which axis to arrange the input/output vector
    A: histogram to store the probability matrix
void GetInput(TH1* inputData, const Int_t* binMap = 0) const
 retreive input distribution
   out: output histogram
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

void GetInputInverseEmatrix(TH2* ematrix)
 calculate the inverse of the contribution to the error matrix
 corresponding to the input data
void GetLsquared(TH2* lsquared) const
 retreive matrix of regularisation conditions squared
   out: prebooked matrix
void GetL(TH2* l) const
 retreive matrix of regularisation conditions
   out: prebooked matrix
void SetConstraint(TUnfold::EConstraint constraint)
 set type of constraint for the next unfolding
Double_t GetTau(void)
 return regularisation parameter
Double_t GetChi2L(void)
 return chi**2 contribution from regularisation conditions
Int_t GetNpar(void)
 return number of parameters
Double_t GetLcurveX(void)
 return value on x axis of L curve
Double_t GetLcurveY(void)
 return value on y axis of L curve
void GetOutput(TH1* output, const Int_t* binMap = 0) const
 get output distribution, cumulated over several bins
   output: output histogram
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

void ErrorMatrixToHist(TH2* ematrix, const TMatrixDSparse* emat, const Int_t* binMap, Bool_t doClear) const
 get an error matrix, cumulated over several bins
   ematrix: output error matrix histogram
   emat: error matrix
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

void GetEmatrix(TH2* ematrix, const Int_t* binMap = 0) const
 get output error matrix, cumulated over several bins
   ematrix: output error matrix histogram
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

void GetRhoIJ(TH2* rhoij, const Int_t* binMap = 0) const
 get correlation coefficient matrix, cumulated over several bins
   rhoij:  correlation coefficient matrix histogram
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

Double_t GetRhoI(TH1* rhoi, const Int_t* binMap = 0, TH2* invEmat = 0) const
 get global correlation coefficients with arbitrary min map
   rhoi: global correlation histogram
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

 return value: maximum global correlation
Double_t GetRhoIFromMatrix(TH1* rhoi, const TMatrixDSparse* eOrig, const Int_t* binMap, TH2* invEmat) const
 get global correlation coefficients with arbitrary min map
   rhoi: global correlation histogram
   emat: error matrix
   binMap: for each bin of the original output distribution
           specify the destination bin. A value of -1 means that the bin
           is discarded. 0 means underflow bin, 1 first bin, ...
        binMap[0] : destination of underflow bin
        binMap[1] : destination of first bin

 return value: maximum global correlation
void ClearHistogram(TH1* h, Double_t x = 0.) const
 clear histogram contents and error
void SetEpsMatrix(Double_t eps)
 set accuracy for matrix inversion
TUnfold(const TUnfold& )
 Int_t IsNotSymmetric(TMatrixDSparse const &m) const;
TMatrixDSparse * MultiplyMSparseTranspMSparse(const TMatrixDSparse* a, const TMatrixDSparse* b) const
TMatrixDSparse * InvertMSparseSymmPos(const TMatrixDSparse* A, Int_t* rank) const
TMatrixDSparse * CreateSparseMatrix(Int_t nrow, Int_t ncol, Int_t nele, Int_t* row, Int_t* col, Double_t* data) const
Int_t GetNx(void)
Int_t GetNy(void)
const TMatrixDSparse * GetDXDY(void)
{ return fDXDY; }
const TMatrixDSparse * GetDXDAM(int i) const
{ return fDXDAM[i]; }
const TMatrixDSparse * GetDXDAZ(int i) const
{ return fDXDAZ[i]; }
const TMatrixDSparse * GetDXDtauSquared(void)
{ return fDXDtauSquared; }
const TMatrixDSparse * GetAx(void)
{ return fAx; }
const TMatrixDSparse * GetEinv(void)
{ return fEinv; }
const TMatrixDSparse * GetE(void)
{ return fE; }
const TMatrixDSparse * GetVxx(void)
{ return fVxx; }
const TMatrixDSparse * GetVxxInv(void)
{ return fVxxInv; }
const TMatrixDSparse * GetVyyInv(void)
{ return fVyyInv; }
const TMatrixD * GetX(void)
{ return fX; }
Int_t GetRowFromBin(int ix) const
{ return fHistToX[ix]; }
Int_t GetBinFromRow(int ix) const
{ return fXToHist[ix]; }
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.)
Bool_t AddRegularisationCondition(Int_t nEle, const Int_t* indices, const Double_t* rowData)
Int_t GetNr(void)
{ return fL->GetNrows(); }
Double_t GetRhoMax(void)
{ return fRhoMax; }
Double_t GetRhoAvg(void)
{ return fRhoAvg; }
Double_t GetChi2A(void)
{ return fChi2A; }
Int_t GetNdf(void)
{ return fNdf; }
Double_t GetEpsMatrix(void)
{ return fEpsMatrix; }