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 virtual ~TSVDUnfold();
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 // Class members
127 Int_t fNdim; //! Truth and reconstructed dimensions
128 Int_t fDdim; //! Derivative for curvature matrix
129 Bool_t fNormalize; //! Normalize unfolded spectrum to 1
130 Int_t fKReg; //! Regularisation parameter
131 TH1D* fDHist; //! Distribution of d (for checking regularization)
132 TH1D* fSVHist; //! Distribution of singular values
133 TH2D* fXtau; //! Computed regularized covariance matrix
134 TH2D* fXinv; //! Computed inverse of covariance matrix
135
136 // Input histos
137 const TH1D* fBdat; // measured distribution (data)
138 TH2D* fBcov; // covariance matrix of measured distribution (data)
139 const TH1D* fBini; // reconstructed distribution (MC)
140 const TH1D* fXini; // truth distribution (MC)
141 const TH2D* fAdet; // Detector response matrix
142
143 // Evaluation of covariance matrices
144 TH1D* fToyhisto; //! Toy MC histogram
145 TH2D* fToymat; //! Toy MC detector response matrix
146 Bool_t fToyMode; //! Internal switch for covariance matrix propagation
147 Bool_t fMatToyMode; //! Internal switch for evaluation of statistical uncertainties from response matrix
148
149
150 ClassDef( TSVDUnfold, 0 ) // Data unfolding using Singular Value Decomposition (hep-ph/9509307)
151};
152
153#endif
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
#define ClassDef(name, id)
Definition: Rtypes.h:326
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
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:593
TH2D * GetBCov() const
Returns the covariance matrix.
Definition: TSVDUnfold.cxx:618
TH1D * fSVHist
Distribution of d (for checking regularization)
Definition: TSVDUnfold.h:132
TH2D * GetXtau() const
Returns the computed regularized covariance matrix corresponding to total uncertainties on measured s...
Definition: TSVDUnfold.cxx:602
Bool_t fToyMode
Toy MC detector response matrix.
Definition: TSVDUnfold.h:146
static void V2H(const TVectorD &vec, TH1D &histo)
Fill vector into 1D histogram.
Definition: TSVDUnfold.cxx:642
static void H2M(const TH2D *histo, TMatrixD &mat)
Fill 2D histogram into matrix.
Definition: TSVDUnfold.cxx:650
static void RegularisedSymMatInvert(TMatrixDSym &mat, Double_t eps=1e-3)
naive regularised inversion cuts off small elements
Definition: TSVDUnfold.cxx:833
static TVectorD CompProd(const TVectorD &vec1, const TVectorD &vec2)
Multiply entries of two vectors.
Definition: TSVDUnfold.cxx:708
static void M2H(const TMatrixD &mat, TH2D &histo)
Fill 2D histogram into matrix.
Definition: TSVDUnfold.cxx:662
Bool_t fMatToyMode
Internal switch for covariance matrix propagation.
Definition: TSVDUnfold.h:147
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:79
Bool_t fNormalize
Derivative for curvature matrix.
Definition: TSVDUnfold.h:129
const TH1D * fBini
Definition: TSVDUnfold.h:139
TH1D * Unfold(Int_t kreg)
Perform the unfolding with regularisation parameter kreg.
Definition: TSVDUnfold.cxx:243
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:886
TH2D * fToymat
Toy MC histogram.
Definition: TSVDUnfold.h:145
static Double_t GetCurvature(const TVectorD &vec, const TMatrixD &curv)
Compute curvature of vector.
Definition: TSVDUnfold.cxx:718
const TH1D * fBdat
Computed inverse of covariance matrix.
Definition: TSVDUnfold.h:137
static TMatrixD MatDivVec(const TMatrixD &mat, const TVectorD &vec, Int_t zero=0)
Divide matrix entries by vector.
Definition: TSVDUnfold.cxx:690
void SetNormalize(Bool_t normalize)
Definition: TSVDUnfold.h:66
virtual ~TSVDUnfold()
Destructor.
Definition: TSVDUnfold.cxx:202
TH2D * fXinv
Computed regularized covariance matrix.
Definition: TSVDUnfold.h:134
static TVectorD VecDiv(const TVectorD &vec1, const TVectorD &vec2, Int_t zero=0)
Divide entries of two vectors.
Definition: TSVDUnfold.cxx:674
void FillCurvatureMatrix(TMatrixD &tCurv, TMatrixD &tC) const
Definition: TSVDUnfold.cxx:725
TH1D * GetD() const
Returns d vector (for choosing appropriate regularisation)
Definition: TSVDUnfold.cxx:582
Int_t GetKReg() const
Definition: TSVDUnfold.h:86
TH2D * fXtau
Distribution of singular values.
Definition: TSVDUnfold.h:133
TH1D * fDHist
Regularisation parameter.
Definition: TSVDUnfold.h:131
static void H2Verr(const TH1D *histo, TVectorD &vec)
Fill 1D histogram errors into vector.
Definition: TSVDUnfold.cxx:634
Int_t fDdim
Truth and reconstructed dimensions.
Definition: TSVDUnfold.h:128
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:411
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:517
const TH2D * fAdet
Definition: TSVDUnfold.h:141
TH2D * GetXinv() const
Returns the computed inverse of the covariance matrix.
Definition: TSVDUnfold.cxx:610
static void H2V(const TH1D *histo, TVectorD &vec)
Fill 1D histogram into vector.
Definition: TSVDUnfold.cxx:626
Int_t fKReg
Normalize unfolded spectrum to 1.
Definition: TSVDUnfold.h:130
TH1D * fToyhisto
Definition: TSVDUnfold.h:144
void InitHistos()
Definition: TSVDUnfold.cxx:813
TH2D * fBcov
Definition: TSVDUnfold.h:138
const TH1D * fXini
Definition: TSVDUnfold.h:140
Int_t fNdim
Definition: TSVDUnfold.h:127