Logo ROOT   6.10/09
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 
43 class TH1D;
44 class TH2D;
45 
46 class TSVDUnfold : public TObject {
47 
48 public:
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 
106 private:
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
Bool_t fMatToyMode
Internal switch for covariance matrix propagation.
Definition: TSVDUnfold.h:147
TH2D * fXtau
Distribution of singular values.
Definition: TSVDUnfold.h:133
SVD Approach to Data Unfolding.
Definition: TSVDUnfold.h:46
const TH1D * fXini
Definition: TSVDUnfold.h:140
static void H2Verr(const TH1D *histo, TVectorD &vec)
Fill 1D histogram errors into vector.
Definition: TSVDUnfold.cxx:634
Bool_t fToyMode
Toy MC detector response matrix.
Definition: TSVDUnfold.h:146
static Double_t GetCurvature(const TVectorD &vec, const TMatrixD &curv)
Compute curvature of vector.
Definition: TSVDUnfold.cxx:718
TH2D * fBcov
Definition: TSVDUnfold.h:138
TH1D * fSVHist
Distribution of d (for checking regularization)
Definition: TSVDUnfold.h:132
TH1D * Unfold(Int_t kreg)
Perform the unfolding with regularisation parameter kreg.
Definition: TSVDUnfold.cxx:243
Int_t fKReg
Normalize unfolded spectrum to 1.
Definition: TSVDUnfold.h:130
void FillCurvatureMatrix(TMatrixD &tCurv, TMatrixD &tC) const
Definition: TSVDUnfold.cxx:725
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
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TH2D * GetXinv() const
Returns the computed inverse of the covariance matrix.
Definition: TSVDUnfold.cxx:610
static TMatrixD MatDivVec(const TMatrixD &mat, const TVectorD &vec, Int_t zero=0)
Divide matrix entries by vector.
Definition: TSVDUnfold.cxx:690
Bool_t fNormalize
Derivative for curvature matrix.
Definition: TSVDUnfold.h:129
TH2D * GetBCov() const
Returns the covariance matrix.
Definition: TSVDUnfold.cxx:618
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
TH2D * fXinv
Computed regularized covariance matrix.
Definition: TSVDUnfold.h:134
TH2D * GetXtau() const
Returns the computed regularized covariance matrix corresponding to total uncertainties on measured s...
Definition: TSVDUnfold.cxx:602
TH1D * GetD() const
Returns d vector (for choosing appropriate regularisation)
Definition: TSVDUnfold.cxx:582
void SetNormalize(Bool_t normalize)
Definition: TSVDUnfold.h:66
TH1D * GetSV() const
Returns singular values vector.
Definition: TSVDUnfold.cxx:593
#define ClassDef(name, id)
Definition: Rtypes.h:297
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
Int_t fDdim
Truth and reconstructed dimensions.
Definition: TSVDUnfold.h:128
static void RegularisedSymMatInvert(TMatrixDSym &mat, Double_t eps=1e-3)
naive regularised inversion cuts off small elements
Definition: TSVDUnfold.cxx:833
static void H2V(const TH1D *histo, TVectorD &vec)
Fill 1D histogram into vector.
Definition: TSVDUnfold.cxx:626
const TH1D * fBdat
Computed inverse of covariance matrix.
Definition: TSVDUnfold.h:137
TH2D * fToymat
Toy MC histogram.
Definition: TSVDUnfold.h:145
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
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
Int_t fNdim
Definition: TSVDUnfold.h:127
static TVectorD VecDiv(const TVectorD &vec1, const TVectorD &vec2, Int_t zero=0)
Divide entries of two vectors.
Definition: TSVDUnfold.cxx:674
Int_t GetKReg() const
Definition: TSVDUnfold.h:86
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:594
static TVectorD CompProd(const TVectorD &vec1, const TVectorD &vec2)
Multiply entries of two vectors.
Definition: TSVDUnfold.cxx:708
const TH1D * fBini
Definition: TSVDUnfold.h:139
double Double_t
Definition: RtypesCore.h:55
TH1D * fDHist
Regularisation parameter.
Definition: TSVDUnfold.h:131
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
void InitHistos()
Definition: TSVDUnfold.cxx:813
Mother of all ROOT objects.
Definition: TObject.h:37
TH1D * fToyhisto
Definition: TSVDUnfold.h:144
static void M2H(const TMatrixD &mat, TH2D &histo)
Fill 2D histogram into matrix.
Definition: TSVDUnfold.cxx:662
virtual ~TSVDUnfold()
Destructor.
Definition: TSVDUnfold.cxx:202
const TH2D * fAdet
Definition: TSVDUnfold.h:141
tomato 2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:290