ROOT 6.07/09 Reference Guide |
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
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
A: migration matrix (dimension ny x nx)
x: unknown underlying distribution (dimension nx) The error matrix of x, V_xx, is also determined
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
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:
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:
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
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
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
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
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
Public Types | |
enum | EConstraint { kEConstraintNone =0, kEConstraintArea =1 } |
enum | EHistMap { kHistMapOutputHoriz = 0, kHistMapOutputVert = 1 } |
enum | ERegMode { kRegModeNone = 0, kRegModeSize = 1, kRegModeDerivative = 2, kRegModeCurvature = 3, kRegModeMixed = 4 } |
Public Types inherited from TObject | |
enum | { kIsOnHeap = 0x01000000, kNotDeleted = 0x02000000, kZombie = 0x04000000, kBitMask = 0x00ffffff } |
enum | { kSingleKey = BIT(0), kOverwrite = BIT(1), kWriteDelete = BIT(2) } |
enum | EStatusBits { kCanDelete = BIT(0), kMustCleanup = BIT(3), kObjInCanvas = BIT(3), kIsReferenced = BIT(4), kHasUUID = BIT(5), kCannotPick = BIT(6), kNoContextMenu = BIT(8), kInvalidObject = BIT(13) } |
Public Member Functions | |
TUnfold (const TH2 *hist_A, EHistMap histmap, ERegMode regmode=kRegModeSize, EConstraint constraint=kEConstraintArea) | |
virtual | ~TUnfold (void) |
Double_t | DoUnfold (Double_t tau, const TH1 *hist_y, Double_t scaleBias=0.0) |
virtual Double_t | DoUnfold (Double_t tau) |
void | GetBias (TH1 *bias, const Int_t *binMap=0) const |
Double_t | GetChi2A (void) const |
Double_t | GetChi2L (void) const |
void | GetEmatrix (TH2 *ematrix, const Int_t *binMap=0) const |
Double_t | GetEpsMatrix (void) const |
void | GetFoldedOutput (TH1 *folded, const Int_t *binMap=0) const |
void | GetInput (TH1 *inputData, const Int_t *binMap=0) const |
void | GetInputInverseEmatrix (TH2 *ematrix) |
void | GetL (TH2 *l) const |
virtual Double_t | GetLcurveX (void) const |
virtual Double_t | GetLcurveY (void) const |
void | GetLsquared (TH2 *lsquared) const |
Int_t | GetNdf (void) const |
void | GetNormalisationVector (TH1 *s, const Int_t *binMap=0) const |
Int_t | GetNpar (void) const |
Int_t | GetNr (void) const |
void | GetOutput (TH1 *output, const Int_t *binMap=0) const |
void | GetProbabilityMatrix (TH2 *A, EHistMap histmap) const |
Double_t | GetRhoAvg (void) const |
Double_t | GetRhoI (TH1 *rhoi, const Int_t *binMap=0, TH2 *invEmat=0) const |
void | GetRhoIJ (TH2 *rhoij, const Int_t *binMap=0) const |
Double_t | GetRhoMax (void) const |
Double_t | GetTau (void) const |
Int_t | RegularizeBins (int start, int step, int nbin, ERegMode regmode) |
Int_t | RegularizeBins2D (int start_bin, int step1, int nbin1, int step2, int nbin2, ERegMode regmode) |
Int_t | RegularizeCurvature (int left_bin, int center_bin, int right_bin, Double_t scale_left=1.0, Double_t scale_right=1.0) |
Int_t | RegularizeDerivative (int left_bin, int right_bin, Double_t scale=1.0) |
Int_t | RegularizeSize (int bin, Double_t scale=1.0) |
virtual Int_t | ScanLcurve (Int_t nPoint, Double_t tauMin, Double_t tauMax, TGraph **lCurve, TSpline **logTauX=0, TSpline **logTauY=0) |
void | SetBias (const TH1 *bias) |
void | SetConstraint (EConstraint constraint) |
void | SetEpsMatrix (Double_t eps) |
virtual Int_t | SetInput (const TH1 *hist_y, Double_t scaleBias=0.0, Double_t oneOverZeroError=0.0, const TH2 *hist_vyy=0, const TH2 *hist_vyy_inv=0) |
Public Member Functions inherited from TObject | |
TObject () | |
TObject constructor. More... | |
TObject (const TObject &object) | |
TObject copy ctor. More... | |
virtual | ~TObject () |
TObject destructor. More... | |
void | AbstractMethod (const char *method) const |
Use this method to implement an "abstract" method that you don't want to leave purely abstract. More... | |
virtual void | AppendPad (Option_t *option="") |
Append graphics object to current pad. More... | |
virtual void | Browse (TBrowser *b) |
Browse object. May be overridden for another default action. More... | |
virtual const char * | ClassName () const |
Returns name of class to which the object belongs. More... | |
virtual void | Clear (Option_t *="") |
virtual TObject * | Clone (const char *newname="") const |
Make a clone of an object using the Streamer facility. More... | |
virtual Int_t | Compare (const TObject *obj) const |
Compare abstract method. More... | |
virtual void | Copy (TObject &object) const |
Copy this to obj. More... | |
virtual void | Delete (Option_t *option="") |
Delete this object. More... | |
virtual Int_t | DistancetoPrimitive (Int_t px, Int_t py) |
Computes distance from point (px,py) to the object. More... | |
virtual void | Draw (Option_t *option="") |
Default Draw method for all objects. More... | |
virtual void | DrawClass () const |
Draw class inheritance tree of the class to which this object belongs. More... | |
virtual TObject * | DrawClone (Option_t *option="") const |
Draw a clone of this object in the current pad. More... | |
virtual void | Dump () const |
Dump contents of object on stdout. More... | |
virtual void | Error (const char *method, const char *msgfmt,...) const |
Issue error message. More... | |
virtual void | Execute (const char *method, const char *params, Int_t *error=0) |
Execute method on this object with the given parameter string, e.g. More... | |
virtual void | Execute (TMethod *method, TObjArray *params, Int_t *error=0) |
Execute method on this object with parameters stored in the TObjArray. More... | |
virtual void | ExecuteEvent (Int_t event, Int_t px, Int_t py) |
Execute action corresponding to an event at (px,py). More... | |
virtual void | Fatal (const char *method, const char *msgfmt,...) const |
Issue fatal error message. More... | |
virtual TObject * | FindObject (const char *name) const |
Must be redefined in derived classes. More... | |
virtual TObject * | FindObject (const TObject *obj) const |
Must be redefined in derived classes. More... | |
virtual Option_t * | GetDrawOption () const |
Get option used by the graphics system to draw this object. More... | |
virtual const char * | GetIconName () const |
Returns mime type name of object. More... | |
virtual const char * | GetName () const |
Returns name of object. More... | |
virtual char * | GetObjectInfo (Int_t px, Int_t py) const |
Returns string containing info about the object at position (px,py). More... | |
virtual Option_t * | GetOption () const |
virtual const char * | GetTitle () const |
Returns title of object. More... | |
virtual UInt_t | GetUniqueID () const |
Return the unique object id. More... | |
virtual Bool_t | HandleTimer (TTimer *timer) |
Execute action in response of a timer timing out. More... | |
virtual ULong_t | Hash () const |
Return hash value for this object. More... | |
virtual void | Info (const char *method, const char *msgfmt,...) const |
Issue info message. More... | |
virtual Bool_t | InheritsFrom (const char *classname) const |
Returns kTRUE if object inherits from class "classname". More... | |
virtual Bool_t | InheritsFrom (const TClass *cl) const |
Returns kTRUE if object inherits from TClass cl. More... | |
virtual void | Inspect () const |
Dump contents of this object in a graphics canvas. More... | |
void | InvertBit (UInt_t f) |
virtual Bool_t | IsEqual (const TObject *obj) const |
Default equal comparison (objects are equal if they have the same address in memory). More... | |
virtual Bool_t | IsFolder () const |
Returns kTRUE in case object contains browsable objects (like containers or lists of other objects). More... | |
Bool_t | IsOnHeap () const |
virtual Bool_t | IsSortable () const |
Bool_t | IsZombie () const |
virtual void | ls (Option_t *option="") const |
The ls function lists the contents of a class on stdout. More... | |
void | MayNotUse (const char *method) const |
Use this method to signal that a method (defined in a base class) may not be called in a derived class (in principle against good design since a child class should not provide less functionality than its parent, however, sometimes it is necessary). More... | |
virtual Bool_t | Notify () |
This method must be overridden to handle object notification. More... | |
void | Obsolete (const char *method, const char *asOfVers, const char *removedFromVers) const |
Use this method to declare a method obsolete. More... | |
void | operator delete (void *ptr) |
Operator delete. More... | |
void | operator delete[] (void *ptr) |
Operator delete []. More... | |
void * | operator new (size_t sz) |
void * | operator new (size_t sz, void *vp) |
void * | operator new[] (size_t sz) |
void * | operator new[] (size_t sz, void *vp) |
TObject & | operator= (const TObject &rhs) |
TObject assignment operator. More... | |
virtual void | Paint (Option_t *option="") |
This method must be overridden if a class wants to paint itself. More... | |
virtual void | Pop () |
Pop on object drawn in a pad to the top of the display list. More... | |
virtual void | Print (Option_t *option="") const |
This method must be overridden when a class wants to print itself. More... | |
virtual Int_t | Read (const char *name) |
Read contents of object with specified name from the current directory. More... | |
virtual void | RecursiveRemove (TObject *obj) |
Recursively remove this object from a list. More... | |
void | ResetBit (UInt_t f) |
virtual void | SaveAs (const char *filename="", Option_t *option="") const |
Save this object in the file specified by filename. More... | |
virtual void | SavePrimitive (std::ostream &out, Option_t *option="") |
Save a primitive as a C++ statement(s) on output stream "out". More... | |
void | SetBit (UInt_t f, Bool_t set) |
Set or unset the user status bits as specified in f. More... | |
void | SetBit (UInt_t f) |
virtual void | SetDrawOption (Option_t *option="") |
Set drawing option for object. More... | |
virtual void | SetUniqueID (UInt_t uid) |
Set the unique object id. More... | |
virtual void | SysError (const char *method, const char *msgfmt,...) const |
Issue system error message. More... | |
Bool_t | TestBit (UInt_t f) const |
Int_t | TestBits (UInt_t f) const |
virtual void | UseCurrentStyle () |
Set current style settings in this object This function is called when either TCanvas::UseCurrentStyle or TROOT::ForceStyle have been invoked. More... | |
virtual void | Warning (const char *method, const char *msgfmt,...) const |
Issue warning message. More... | |
virtual Int_t | Write (const char *name=0, Int_t option=0, Int_t bufsize=0) |
Write this object to the current directory. More... | |
virtual Int_t | Write (const char *name=0, Int_t option=0, Int_t bufsize=0) const |
Write this object to the current directory. More... | |
Static Public Member Functions | |
static const char * | GetTUnfoldVersion (void) |
Static Public Member Functions inherited from TObject | |
static Long_t | GetDtorOnly () |
Return destructor only flag. More... | |
static Bool_t | GetObjectStat () |
Get status of object stat flag. More... | |
static void | SetDtorOnly (void *obj) |
Set destructor only flag. More... | |
static void | SetObjectStat (Bool_t stat) |
Turn on/off tracking of objects in the TObjectTable. More... | |
Static Protected Member Functions | |
static void | DeleteMatrix (TMatrixD **m) |
static void | DeleteMatrix (TMatrixDSparse **m) |
Protected Attributes | |
TMatrixDSparse * | fA |
Double_t | fBiasScale |
EConstraint | fConstraint |
TArrayI | fHistToX |
TMatrixDSparse * | fL |
ERegMode | fRegMode |
TArrayD | fSumOverY |
Double_t | fTauSquared |
TMatrixDSparse * | fVyy |
TMatrixD * | fX0 |
TArrayI | fXToHist |
TMatrixD * | fY |
Private Member Functions | |
void | InitTUnfold (void) |
Private Attributes | |
TMatrixDSparse * | fAx |
Double_t | fChi2A |
TMatrixDSparse * | fDXDAM [2] |
TMatrixDSparse * | fDXDAZ [2] |
TMatrixDSparse * | fDXDtauSquared |
TMatrixDSparse * | fDXDY |
TMatrixDSparse * | fE |
TMatrixDSparse * | fEinv |
Double_t | fEpsMatrix |
Int_t | fIgnoredBins |
Double_t | fLXsquared |
Int_t | fNdf |
Double_t | fRhoAvg |
Double_t | fRhoMax |
TMatrixDSparse * | fVxx |
TMatrixDSparse * | fVxxInv |
TMatrixDSparse * | fVyyInv |
TMatrixD * | fX |
#include <TUnfold.h>
enum TUnfold::EConstraint |
enum TUnfold::EHistMap |
enum TUnfold::ERegMode |
|
protected |
Definition at line 372 of file TUnfold.cxx.
TUnfold::TUnfold | ( | const TH2 * | hist_A, |
EHistMap | histmap, | ||
ERegMode | regmode = kRegModeSize , |
||
EConstraint | constraint = kEConstraintArea |
||
) |
Definition at line 1745 of file TUnfold.cxx.
|
virtual |
Definition at line 1948 of file TUnfold.cxx.
|
protected |
Definition at line 1001 of file TUnfold.cxx.
|
protected |
Definition at line 1973 of file TUnfold.cxx.
|
protected |
Definition at line 2007 of file TUnfold.cxx.
Definition at line 3394 of file TUnfold.cxx.
Reimplemented in TUnfoldSys.
Definition at line 345 of file TUnfold.cxx.
|
protected |
Definition at line 698 of file TUnfold.cxx.
Definition at line 333 of file TUnfold.cxx.
|
staticprotected |
Definition at line 339 of file TUnfold.cxx.
Definition at line 378 of file TUnfold.cxx.
Definition at line 2186 of file TUnfold.cxx.
Definition at line 2407 of file TUnfold.cxx.
|
protected |
Definition at line 3109 of file TUnfold.cxx.
|
inlineprotected |
Definition at line 2794 of file TUnfold.cxx.
|
inlineprotected |
Definition at line 3009 of file TUnfold.cxx.
|
inlineprotected |
|
inlineprotected |
|
inlineprotected |
|
inlineprotected |
|
inlineprotected |
|
inlineprotected |
Definition at line 3175 of file TUnfold.cxx.
Definition at line 2815 of file TUnfold.cxx.
Definition at line 2884 of file TUnfold.cxx.
Definition at line 2917 of file TUnfold.cxx.
Definition at line 2977 of file TUnfold.cxx.
Definition at line 3021 of file TUnfold.cxx.
Definition at line 3027 of file TUnfold.cxx.
Definition at line 2958 of file TUnfold.cxx.
Definition at line 2774 of file TUnfold.cxx.
Definition at line 3015 of file TUnfold.cxx.
Definition at line 3033 of file TUnfold.cxx.
Reimplemented in TUnfoldDensity.
Definition at line 1737 of file TUnfold.cxx.
Definition at line 2863 of file TUnfold.cxx.
Definition at line 3216 of file TUnfold.cxx.
|
protected |
Definition at line 3275 of file TUnfold.cxx.
Definition at line 3188 of file TUnfold.cxx.
|
inlineprotected |
Definition at line 3003 of file TUnfold.cxx.
|
static |
Definition at line 289 of file TUnfold.cxx.
|
inlineprotected |
|
inlineprotected |
|
inlineprotected |
Definition at line 295 of file TUnfold.cxx.
|
protected |
Definition at line 1064 of file TUnfold.cxx.
|
protected |
Definition at line 858 of file TUnfold.cxx.
|
protected |
Definition at line 711 of file TUnfold.cxx.
|
protected |
Definition at line 911 of file TUnfold.cxx.
|
protected |
Definition at line 780 of file TUnfold.cxx.
Definition at line 2123 of file TUnfold.cxx.
Int_t TUnfold::RegularizeBins2D | ( | int | start_bin, |
int | step1, | ||
int | nbin1, | ||
int | step2, | ||
int | nbin2, | ||
ERegMode | regmode | ||
) |
Definition at line 2163 of file TUnfold.cxx.
Int_t TUnfold::RegularizeCurvature | ( | int | left_bin, |
int | center_bin, | ||
int | right_bin, | ||
Double_t | scale_left = 1.0 , |
||
Double_t | scale_right = 1.0 |
||
) |
Definition at line 2099 of file TUnfold.cxx.
Definition at line 2083 of file TUnfold.cxx.
Definition at line 2069 of file TUnfold.cxx.
|
virtual |
Definition at line 2424 of file TUnfold.cxx.
Definition at line 1961 of file TUnfold.cxx.
void TUnfold::SetConstraint | ( | EConstraint | constraint | ) |
Definition at line 2995 of file TUnfold.cxx.
Definition at line 3418 of file TUnfold.cxx.
|
virtual |
Reimplemented in TUnfoldSys.
Definition at line 2208 of file TUnfold.cxx.
|
protected |
|
private |
|
protected |
|
private |
|
private |
|
private |
|
private |
|
private |
|
private |
|
protected |
|
private |
|
private |
|
protected |
|
private |