ROOT logo
ROOT » HIST » HIST » TSVDUnfold

class TSVDUnfold: public TObject


SVD Approach to Data Unfolding

Reference: Nucl. Instrum. Meth. A372, 469 (1996) [hep-ph/9509307]

TSVDUnfold implements the singular value decomposition based unfolding method (see reference). Currently, the unfolding of one-dimensional histograms is supported, with the same number of bins for the measured and the unfolded spectrum.

The unfolding procedure is based on singular value decomposition of the response matrix. The regularisation of the unfolding is implemented via a discrete minimum-curvature condition.

Monte Carlo inputs:

  • xini: true underlying spectrum (TH1D, n bins)
  • bini: reconstructed spectrum (TH1D, n bins)
  • Adet: response matrix (TH2D, nxn bins)
Consider the unfolding of a measured spectrum bdat with covariance matrix Bcov (if not passed explicitly, a diagonal covariance will be built given the errors of bdat). The corresponding spectrum in the Monte Carlo is given by bini, with the true underlying spectrum given by xini. The detector response is described by Adet, with Adet filled with events (not probabilities) with the true observable on the y-axis and the reconstructed observable on the x-axis.

The measured distribution can be unfolded for any combination of resolution, efficiency and acceptance effects, provided an appropriate definition of xini and Adet.

The unfolding can be performed by

     TSVDUnfold *tsvdunf = new TSVDUnfold( bdat, Bcov, bini, xini, Adet );
     TH1D* unfresult = tsvdunf->Unfold( kreg );
     
where kreg determines the regularisation of the unfolding. In general, overregularisation (too small kreg) will bias the unfolded spectrum towards the Monte Carlo input, while underregularisation (too large kreg) will lead to large fluctuations in the unfolded spectrum. The optimal regularisation can be determined following guidelines in Nucl. Instrum. Meth. A372, 469 (1996) [hep-ph/9509307] using the distribution of the |d_i|<\tt> that can be obtained by tsvdunf->GetD() and/or using pseudo-experiments.

Covariance matrices on the measured spectrum (for either the total uncertainties or individual sources of uncertainties) can be propagated to covariance matrices using the GetUnfoldCovMatrix method, which uses pseudo experiments for the propagation. In addition, GetAdetCovMatrix allows for the propagation of the statistical uncertainties on the response matrix using pseudo experiments. The covariance matrix corresponding to Bcov is also computed as described in Nucl. Instrum. Meth. A372, 469 (1996) [hep-ph/9509307] and can be obtained from tsvdunf->GetXtau() and its (regularisation independent) inverse from tsvdunf->GetXinv(). The distribution of singular values can be retrieved using tsvdunf->GetSV().

See also the tutorial for a toy example.

 

Function Members (Methods)

public:
TSVDUnfold(const TSVDUnfold& other)
TSVDUnfold(const TH1D* bdat, const TH1D* bini, const TH1D* xini, const TH2D* Adet)
TSVDUnfold(const TH1D* bdat, TH2D* Bcov, const TH1D* bini, const TH1D* xini, const TH2D* Adet)
virtual~TSVDUnfold()
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
Double_tComputeChiSquared(const TH1D& truspec, const TH1D& unfspec)
virtual voidTObject::Copy(TObject& object) const
virtual voidTObject::Delete(Option_t* option = "")MENU
virtual Int_tTObject::DistancetoPrimitive(Int_t px, Int_t py)
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
TH2D*GetAdetCovMatrix(Int_t ntoys, Int_t seed = 1)
TH2D*GetBCov() const
TH1D*GetD() const
virtual Option_t*TObject::GetDrawOption() const
static Long_tTObject::GetDtorOnly()
virtual const char*TObject::GetIconName() const
Int_tGetKReg() const
virtual const char*TObject::GetName() const
virtual char*TObject::GetObjectInfo(Int_t px, Int_t py) const
static Bool_tTObject::GetObjectStat()
virtual Option_t*TObject::GetOption() const
TH1D*GetSV() const
virtual const char*TObject::GetTitle() const
TH2D*GetUnfoldCovMatrix(const TH2D* cov, Int_t ntoys, Int_t seed = 1)
virtual UInt_tTObject::GetUniqueID() const
TH2D*GetXinv() const
TH2D*GetXtau() 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
static voidTObject::operator delete(void* ptr)
static voidTObject::operator delete(void* ptr, void* vp)
static voidTObject::operator delete[](void* ptr)
static 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)
TObject&TObject::operator=(const TObject& rhs)
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)
voidTObject::ResetBit(UInt_t f)
virtual voidTObject::SaveAs(const char* filename = "", Option_t* option = "") constMENU
virtual voidTObject::SavePrimitive(ostream& out, Option_t* option = "")
voidTObject::SetBit(UInt_t f)
voidTObject::SetBit(UInt_t f, Bool_t set)
virtual voidTObject::SetDrawOption(Option_t* option = "")MENU
static voidTObject::SetDtorOnly(void* obj)
voidSetNormalize(Bool_t normalize)
static voidTObject::SetObjectStat(Bool_t stat)
virtual voidTObject::SetUniqueID(UInt_t uid)
virtual voidShowMembers(TMemberInspector&)
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
TH1D*Unfold(Int_t kreg)
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:
virtual voidTObject::DoError(int level, const char* location, const char* fmt, va_list va) const
voidTObject::MakeZombie()
private:
static TVectorDCompProd(const TVectorD& vec1, const TVectorD& vec2)
voidFillCurvatureMatrix(TMatrixD& tCurv, TMatrixD& tC) const
static Double_tGetCurvature(const TVectorD& vec, const TMatrixD& curv)
static voidH2M(const TH2D* histo, TMatrixD& mat)
static voidH2V(const TH1D* histo, TVectorD& vec)
static voidH2Verr(const TH1D* histo, TVectorD& vec)
voidInitHistos()
static voidM2H(const TMatrixD& mat, TH2D& histo)
static TMatrixDMatDivVec(const TMatrixD& mat, const TVectorD& vec, Int_t zero = 0)
static voidRegularisedSymMatInvert(TMatrixDSym& mat, Double_t eps = 1e-3)
static voidV2H(const TVectorD& vec, TH1D& histo)
static TVectorDVecDiv(const TVectorD& vec1, const TVectorD& vec2, Int_t zero = 0)

Data Members

private:
const TH2D*fAdetDetector response matrix
TH2D*fBcovcovariance matrix of measured distribution (data)
const TH1D*fBdatmeasured distribution (data)
const TH1D*fBinireconstructed distribution (MC)
TH1D*fDHist! Distribution of d (for checking regularization)
Int_tfDdim! Derivative for curvature matrix
Int_tfKReg! Regularisation parameter
Bool_tfMatToyMode! Internal switch for evaluation of statistical uncertainties from response matrix
Int_tfNdim! Truth and reconstructed dimensions
Bool_tfNormalize! Normalize unfolded spectrum to 1
TH1D*fSVHist! Distribution of singular values
Bool_tfToyMode! Internal switch for covariance matrix propagation
TH1D*fToyhisto! Toy MC histogram
TH2D*fToymat! Toy MC detector response matrix
const TH1D*fXinitruth distribution (MC)
TH2D*fXinv! Computed inverse of covariance matrix
TH2D*fXtau! Computed regularized covariance matrix

Class Charts

Class Charts

Function documentation

TSVDUnfold(const TH1D* bdat, const TH1D* bini, const TH1D* xini, const TH2D* Adet)
 Alternative constructor
 User provides data and MC test spectra, as well as detector response matrix, diagonal covariance matrix of measured spectrum built from the uncertainties on measured spectrum
TSVDUnfold(const TH1D* bdat, TH2D* Bcov, const TH1D* bini, const TH1D* xini, const TH2D* Adet)
 Default constructor
 Initialisation of TSVDUnfold
 User provides data and MC test spectra, as well as detector response matrix and the covariance matrix of the measured distribution
TSVDUnfold(const TSVDUnfold& other)
 Copy constructor
~TSVDUnfold()
 Destructor
TH1D* Unfold(Int_t kreg)
 Perform the unfolding with regularisation parameter kreg
TH2D* GetUnfoldCovMatrix(const TH2D* cov, Int_t ntoys, Int_t seed = 1)
 Determine for given input error matrix covariance matrix of unfolded
 spectrum from toy simulation given the passed covariance matrix on measured spectrum
 "cov"    - covariance matrix on the measured spectrum, to be propagated
 "ntoys"  - number of pseudo experiments used for the propagation
 "seed"   - seed for pseudo experiments
 Note that this covariance matrix will contain effects of forced normalisation if spectrum is normalised to unit area.
TH2D* GetAdetCovMatrix(Int_t ntoys, Int_t seed = 1)
 Determine covariance matrix of unfolded spectrum from finite statistics in
 response matrix using pseudo experiments
 "ntoys"  - number of pseudo experiments used for the propagation
 "seed"   - seed for pseudo experiments
TH1D* GetD() const
 Returns d vector (for choosing appropriate regularisation)
TH1D* GetSV() const
 Returns singular values vector
TH2D* GetXtau() const
 Returns the computed regularized covariance matrix corresponding to total uncertainties on measured spectrum as passed in the constructor.
 Note that this covariance matrix will not contain the effects of forced normalization if spectrum is normalized to unit area.
TH2D* GetXinv() const
 Returns the computed inverse of the covariance matrix
TH2D* GetBCov() const
 Returns the covariance matrix
void H2V(const TH1D* histo, TVectorD& vec)
 Fill 1D histogram into vector
void H2Verr(const TH1D* histo, TVectorD& vec)
 Fill 1D histogram errors into vector
void V2H(const TVectorD& vec, TH1D& histo)
 Fill vector into 1D histogram
void H2M(const TH2D* histo, TMatrixD& mat)
 Fill 2D histogram into matrix
void M2H(const TMatrixD& mat, TH2D& histo)
 Fill 2D histogram into matrix
TVectorD VecDiv(const TVectorD& vec1, const TVectorD& vec2, Int_t zero = 0)
 Divide entries of two vectors
TMatrixD MatDivVec(const TMatrixD& mat, const TVectorD& vec, Int_t zero = 0)
 Divide matrix entries by vector
TVectorD CompProd(const TVectorD& vec1, const TVectorD& vec2)
 Multiply entries of two vectors
Double_t GetCurvature(const TVectorD& vec, const TMatrixD& curv)
 Compute curvature of vector
void FillCurvatureMatrix(TMatrixD& tCurv, TMatrixD& tC) const
void InitHistos()
void RegularisedSymMatInvert(TMatrixDSym& mat, Double_t eps = 1e-3)
 naive regularised inversion cuts off small elements
Double_t ComputeChiSquared(const TH1D& truspec, const TH1D& unfspec)
 Helper routine to compute chi-squared between distributions using the computed inverse of the covariance matrix for the unfolded spectrum as given in paper.
void SetNormalize(Bool_t normalize)
 Set option to normalize unfolded spectrum to unit area
 "normalize" - switch
{ fNormalize = normalize; }
Int_t GetKReg() const
 Regularisation parameter
{ return fKReg; }