Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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