Logo ROOT  
Reference Guide
TSVDUnfold.h
Go to the documentation of this file.
1// Author: Kerstin Tackmann, Andreas Hoecker, Heiko Lacker
2
3/**********************************************************************************
4 * *
5 * Project: TSVDUnfold - data unfolding based on Singular Value Decomposition *
6 * Package: ROOT *
7 * Class : TSVDUnfold *
8 * *
9 * Description: *
10 * Single class implementation of SVD data unfolding based on: *
11 * A. Hoecker, V. Kartvelishvili, *
12 * "SVD approach to data unfolding" *
13 * NIM A372, 469 (1996) [hep-ph/9509307] *
14 * *
15 * Authors: *
16 * Kerstin Tackmann <Kerstin.Tackmann@cern.ch> - CERN, Switzerland *
17 * Andreas Hoecker <Andreas.Hoecker@cern.ch> - CERN, Switzerland *
18 * Heiko Lacker <lacker@physik.hu-berlin.de> - Humboldt U, Germany *
19 * *
20 * Copyright (c) 2010: *
21 * CERN, Switzerland *
22 * Humboldt University, Germany *
23 * *
24 **********************************************************************************/
25
26//////////////////////////////////////////////////////////////////////////
27// //
28// TSVDUnfold //
29// //
30// Data unfolding using Singular Value Decomposition (hep-ph/9509307) //
31// Authors: Kerstin Tackmann, Andreas Hoecker, Heiko Lacker //
32// //
33//////////////////////////////////////////////////////////////////////////
34
35#ifndef TSVDUNFOLD_H
36#define TSVDUNFOLD_H
37
38#include "TObject.h"
39#include "TMatrixD.h"
40#include "TVectorD.h"
41#include "TMatrixDSym.h"
42
43class TH1D;
44class TH2D;
45
46class TSVDUnfold : public TObject {
47
48public:
49
50 // Constructor
51 // Initialisation of unfolding
52 // "bdat" - measured data distribution (number of events)
53 // "Bcov" - covariance matrix for measured data distribution
54 // "bini" - reconstructed MC distribution (number of events)
55 // "xini" - truth MC distribution (number of events)
56 // "Adet" - detector response matrix (number of events)
57 TSVDUnfold( const TH1D* bdat, const TH1D* bini, const TH1D* xini, const TH2D* Adet );
58 TSVDUnfold( const TH1D* bdat, TH2D* Bcov, const TH1D* bini, const TH1D* xini, const TH2D* Adet );
59 TSVDUnfold( const TSVDUnfold& other );
60
61 // Destructor
62 ~TSVDUnfold() override;
63
64 // Set option to normalize unfolded spectrum to unit area
65 // "normalize" - switch
66 void SetNormalize ( Bool_t normalize ) { fNormalize = normalize; }
67
68 // Do the unfolding
69 // "kreg" - number of singular values used (regularisation)
70 TH1D* Unfold ( Int_t kreg );
71
72 // Determine for given input error matrix covariance matrix of unfolded
73 // spectrum from toy simulation
74 // "cov" - covariance matrix on the measured spectrum, to be propagated
75 // "ntoys" - number of pseudo experiments used for the propagation
76 // "seed" - seed for pseudo experiments
77 TH2D* GetUnfoldCovMatrix( const TH2D* cov, Int_t ntoys, Int_t seed = 1 );
78
79 // Determine covariance matrix of unfolded spectrum from finite statistics in
80 // response matrix
81 // "ntoys" - number of pseudo experiments used for the propagation
82 // "seed" - seed for pseudo experiments
83 TH2D* GetAdetCovMatrix( Int_t ntoys, Int_t seed=1 );
84
85 // Regularisation parameter
86 Int_t GetKReg() const { return fKReg; }
87
88 // Obtain the distribution of |d| (for determining the regularization)
89 TH1D* GetD() const;
90
91 // Obtain the distribution of singular values
92 TH1D* GetSV() const;
93
94 // Obtain the computed regularized covariance matrix
95 TH2D* GetXtau() const;
96
97 // Obtain the computed inverse of the covariance matrix
98 TH2D* GetXinv() const;
99
100 //Obtain the covariance matrix on the data
101 TH2D* GetBCov() const;
102
103 // Helper functions
104 Double_t ComputeChiSquared( const TH1D& truspec, const TH1D& unfspec );
105
106private:
107
108 // Helper functions for vector and matrix operations
109 void FillCurvatureMatrix( TMatrixD& tCurv, TMatrixD& tC ) const;
110 static Double_t GetCurvature ( const TVectorD& vec, const TMatrixD& curv );
111
112 void InitHistos ( );
113
114 // Helper functions
115 static void H2V ( const TH1D* histo, TVectorD& vec );
116 static void H2Verr ( const TH1D* histo, TVectorD& vec );
117 static void V2H ( const TVectorD& vec, TH1D& histo );
118 static void H2M ( const TH2D* histo, TMatrixD& mat );
119 static void M2H ( const TMatrixD& mat, TH2D& histo );
120 static TMatrixD MatDivVec( const TMatrixD& mat, const TVectorD& vec, Int_t zero=0 );
121 static TVectorD CompProd ( const TVectorD& vec1, const TVectorD& vec2 );
122
123 static TVectorD VecDiv ( const TVectorD& vec1, const TVectorD& vec2, Int_t zero = 0 );
124 static void RegularisedSymMatInvert( TMatrixDSym& mat, Double_t eps = 1e-3 );
125
126 /// @name Class members
127 ///@{
128 Int_t fNdim; ///<! Truth and reconstructed dimensions
129 Int_t fDdim; ///<! Derivative for curvature matrix
130 Bool_t fNormalize; ///<! Normalize unfolded spectrum to 1
131 Int_t fKReg; ///<! Regularisation parameter
132 TH1D* fDHist; ///<! Distribution of d (for checking regularization)
133 TH1D* fSVHist; ///<! Distribution of singular values
134 TH2D* fXtau; ///<! Computed regularized covariance matrix
135 TH2D* fXinv; ///<! Computed inverse of covariance matrix
136 ///@}
137
138 /// @name Input histos
139 ///@{
140 const TH1D* fBdat; ///< Measured distribution (data)
141 TH2D* fBcov; ///< Covariance matrix of measured distribution (data)
142 const TH1D* fBini; ///< Reconstructed distribution (MC)
143 const TH1D* fXini; ///< Truth distribution (MC)
144 const TH2D* fAdet; ///< Detector response matrix
145 ///@}
146
147 /// @name Evaluation of covariance matrices
148 ///@{
149 TH1D* fToyhisto; ///<! Toy MC histogram
150 TH2D* fToymat; ///<! Toy MC detector response matrix
151 Bool_t fToyMode; ///<! Internal switch for covariance matrix propagation
152 Bool_t fMatToyMode; ///<! Internal switch for evaluation of statistical uncertainties from response matrix
153 ///@}
154
155
156 ClassDefOverride( TSVDUnfold, 0 ) // Data unfolding using Singular Value Decomposition (hep-ph/9509307)
157};
158
159#endif
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:45
double Double_t
Definition: RtypesCore.h:59
#define ClassDefOverride(name, id)
Definition: Rtypes.h:339
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:617
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:300
Mother of all ROOT objects.
Definition: TObject.h:37
SVD Approach to Data Unfolding.
Definition: TSVDUnfold.h:46
TH1D * GetSV() const
Returns singular values vector.
Definition: TSVDUnfold.cxx:591
TH2D * GetBCov() const
Returns the covariance matrix.
Definition: TSVDUnfold.cxx:616
TH1D * fSVHist
! Distribution of singular values
Definition: TSVDUnfold.h:133
TH2D * GetXtau() const
Returns the computed regularized covariance matrix corresponding to total uncertainties on measured s...
Definition: TSVDUnfold.cxx:600
Bool_t fToyMode
! Internal switch for covariance matrix propagation
Definition: TSVDUnfold.h:151
static void V2H(const TVectorD &vec, TH1D &histo)
Fill vector into 1D histogram.
Definition: TSVDUnfold.cxx:640
static void H2M(const TH2D *histo, TMatrixD &mat)
Fill 2D histogram into matrix.
Definition: TSVDUnfold.cxx:648
static void RegularisedSymMatInvert(TMatrixDSym &mat, Double_t eps=1e-3)
naive regularised inversion cuts off small elements
Definition: TSVDUnfold.cxx:831
static TVectorD CompProd(const TVectorD &vec1, const TVectorD &vec2)
Multiply entries of two vectors.
Definition: TSVDUnfold.cxx:706
static void M2H(const TMatrixD &mat, TH2D &histo)
Fill 2D histogram into matrix.
Definition: TSVDUnfold.cxx:660
Bool_t fMatToyMode
! Internal switch for evaluation of statistical uncertainties from response matrix
Definition: TSVDUnfold.h:152
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,...
Definition: TSVDUnfold.cxx:77
Bool_t fNormalize
! Normalize unfolded spectrum to 1
Definition: TSVDUnfold.h:130
const TH1D * fBini
Reconstructed distribution (MC)
Definition: TSVDUnfold.h:142
TH1D * Unfold(Int_t kreg)
Perform the unfolding with regularisation parameter kreg.
Definition: TSVDUnfold.cxx:241
Double_t ComputeChiSquared(const TH1D &truspec, const TH1D &unfspec)
Helper routine to compute chi-squared between distributions using the computed inverse of the covaria...
Definition: TSVDUnfold.cxx:884
TH2D * fToymat
! Toy MC detector response matrix
Definition: TSVDUnfold.h:150
static Double_t GetCurvature(const TVectorD &vec, const TMatrixD &curv)
Compute curvature of vector.
Definition: TSVDUnfold.cxx:716
const TH1D * fBdat
Measured distribution (data)
Definition: TSVDUnfold.h:140
~TSVDUnfold() override
Destructor.
Definition: TSVDUnfold.cxx:200
static TMatrixD MatDivVec(const TMatrixD &mat, const TVectorD &vec, Int_t zero=0)
Divide matrix entries by vector.
Definition: TSVDUnfold.cxx:688
void SetNormalize(Bool_t normalize)
Definition: TSVDUnfold.h:66
TH2D * fXinv
! Computed inverse of covariance matrix
Definition: TSVDUnfold.h:135
static TVectorD VecDiv(const TVectorD &vec1, const TVectorD &vec2, Int_t zero=0)
Divide entries of two vectors.
Definition: TSVDUnfold.cxx:672
void FillCurvatureMatrix(TMatrixD &tCurv, TMatrixD &tC) const
Definition: TSVDUnfold.cxx:723
TH1D * GetD() const
Returns d vector (for choosing appropriate regularisation)
Definition: TSVDUnfold.cxx:580
Int_t GetKReg() const
Definition: TSVDUnfold.h:86
TH2D * fXtau
! Computed regularized covariance matrix
Definition: TSVDUnfold.h:134
TH1D * fDHist
! Distribution of d (for checking regularization)
Definition: TSVDUnfold.h:132
static void H2Verr(const TH1D *histo, TVectorD &vec)
Fill 1D histogram errors into vector.
Definition: TSVDUnfold.cxx:632
Int_t fDdim
! Derivative for curvature matrix
Definition: TSVDUnfold.h:129
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 giv...
Definition: TSVDUnfold.cxx:409
TH2D * GetAdetCovMatrix(Int_t ntoys, Int_t seed=1)
Determine covariance matrix of unfolded spectrum from finite statistics in response matrix using pseu...
Definition: TSVDUnfold.cxx:515
const TH2D * fAdet
Detector response matrix.
Definition: TSVDUnfold.h:144
TH2D * GetXinv() const
Returns the computed inverse of the covariance matrix.
Definition: TSVDUnfold.cxx:608
static void H2V(const TH1D *histo, TVectorD &vec)
Fill 1D histogram into vector.
Definition: TSVDUnfold.cxx:624
Int_t fKReg
! Regularisation parameter
Definition: TSVDUnfold.h:131
TH1D * fToyhisto
! Toy MC histogram
Definition: TSVDUnfold.h:149
void InitHistos()
Definition: TSVDUnfold.cxx:811
TH2D * fBcov
Covariance matrix of measured distribution (data)
Definition: TSVDUnfold.h:141
const TH1D * fXini
Truth distribution (MC)
Definition: TSVDUnfold.h:143
Int_t fNdim
! Truth and reconstructed dimensions
Definition: TSVDUnfold.h:128
Definition: civetweb.c:1856