00001 // Author: Kerstin Tackmann, Andreas Hoecker, Heiko Lacker 00002 00003 /********************************************************************************** 00004 * * 00005 * Project: TSVDUnfold - data unfolding based on Singular Value Decomposition * 00006 * Package: ROOT * 00007 * Class : TSVDUnfold * 00008 * * 00009 * Description: * 00010 * Single class implementation of SVD data unfolding based on: * 00011 * A. Hoecker, V. Kartvelishvili, * 00012 * "SVD approach to data unfolding" * 00013 * NIM A372, 469 (1996) [hep-ph/9509307] * 00014 * * 00015 * Authors: * 00016 * Kerstin Tackmann <Kerstin.Tackmann@cern.ch> - CERN, Switzerland * 00017 * Andreas Hoecker <Andreas.Hoecker@cern.ch> - CERN, Switzerland * 00018 * Heiko Lacker <lacker@physik.hu-berlin.de> - Humboldt U, Germany * 00019 * * 00020 * Copyright (c) 2010: * 00021 * CERN, Switzerland * 00022 * Humboldt University, Germany * 00023 * * 00024 **********************************************************************************/ 00025 00026 ////////////////////////////////////////////////////////////////////////// 00027 // // 00028 // TSVDUnfold // 00029 // // 00030 // Data unfolding using Singular Value Decomposition (hep-ph/9509307) // 00031 // Authors: Kerstin Tackmann, Andreas Hoecker, Heiko Lacker // 00032 // // 00033 ////////////////////////////////////////////////////////////////////////// 00034 00035 #ifndef TSVDUNFOLD_H 00036 #define TSVDUNFOLD_H 00037 00038 #ifndef ROOT_TObject 00039 #include "TObject.h" 00040 #endif 00041 #ifndef ROOT_TMatrixD 00042 #include "TMatrixD.h" 00043 #endif 00044 #ifndef ROOT_TVectorD 00045 #include "TVectorD.h" 00046 #endif 00047 #ifndef ROOT_TMatrixDSym 00048 #include "TMatrixDSym.h" 00049 #endif 00050 00051 class TH1D; 00052 class TH2D; 00053 00054 class TSVDUnfold : public TObject { 00055 00056 public: 00057 00058 // Constructor 00059 // Initialisation of unfolding 00060 // "bdat" - measured data distribution (number of events) 00061 // "bini" - reconstructed MC distribution (number of events) 00062 // "xini" - truth MC distribution (number of events) 00063 // "Adet" - detector response matrix (number of events) 00064 TSVDUnfold( const TH1D* bdat, const TH1D* bini, const TH1D* xini, const TH2D* Adet ); 00065 TSVDUnfold( const TSVDUnfold& other ); 00066 00067 // Destructor 00068 virtual ~TSVDUnfold(); 00069 00070 // Set option to normalize unfolded spectrum to unit area 00071 // "normalize" - switch 00072 void SetNormalize ( Bool_t normalize ) { fNormalize = normalize; } 00073 00074 // Do the unfolding 00075 // "kreg" - number of singular values used (regularisation) 00076 TH1D* Unfold ( Int_t kreg ); 00077 00078 // Determine for given input error matrix covariance matrix of unfolded 00079 // spectrum from toy simulation 00080 // "cov" - covariance matrix on the measured spectrum, to be propagated 00081 // "ntoys" - number of pseudo experiments used for the propagation 00082 // "seed" - seed for pseudo experiments 00083 TH2D* GetUnfoldCovMatrix( const TH2D* cov, Int_t ntoys, Int_t seed = 1 ); 00084 00085 // Determine covariance matrix of unfolded spectrum from finite statistics in 00086 // response matrix 00087 // "ntoys" - number of pseudo experiments used for the propagation 00088 // "seed" - seed for pseudo experiments 00089 TH2D* GetAdetCovMatrix( Int_t ntoys, Int_t seed=1 ); 00090 00091 // Regularisation parameter 00092 Int_t GetKReg() const { return fKReg; } 00093 00094 // Obtain the distribution of |d| (for determining the regularization) 00095 TH1D* GetD() const; 00096 00097 // Obtain the distribution of singular values 00098 TH1D* GetSV() const; 00099 00100 // Helper functions 00101 static Double_t ComputeChiSquared( const TH1D& truspec, const TH1D& unfspec, const TH2D& covmat, Double_t regpar = 0.01 ); 00102 00103 private: 00104 00105 // Helper functions for vector and matrix operations 00106 void FillCurvatureMatrix( TMatrixD& tCurv, TMatrixD& tC ) const; 00107 static Double_t GetCurvature ( const TVectorD& vec, const TMatrixD& curv ); 00108 00109 void InitHistos ( ); 00110 00111 // Helper functions 00112 static void H2V ( const TH1D* histo, TVectorD& vec ); 00113 static void H2Verr ( const TH1D* histo, TVectorD& vec ); 00114 static void V2H ( const TVectorD& vec, TH1D& histo ); 00115 static void H2M ( const TH2D* histo, TMatrixD& mat ); 00116 static TMatrixD MatDivVec( const TMatrixD& mat, const TVectorD& vec, Int_t zero=0 ); 00117 static TVectorD CompProd ( const TVectorD& vec1, const TVectorD& vec2 ); 00118 00119 static TVectorD VecDiv ( const TVectorD& vec1, const TVectorD& vec2, Int_t zero = 0 ); 00120 static void RegularisedSymMatInvert( TMatrixDSym& mat, Double_t eps = 1e-3 ); 00121 00122 // Class members 00123 Int_t fNdim; //! Truth and reconstructed dimensions 00124 Int_t fDdim; //! Derivative for curvature matrix 00125 Bool_t fNormalize; //! Normalize unfolded spectrum to 1 00126 Int_t fKReg; //! Regularisation parameter 00127 TH1D* fDHist; // Distribution of d (for checking regularization) 00128 TH1D* fSVHist; // Distribution of singular values 00129 00130 // Input histos 00131 const TH1D* fBdat; // measured distribution (data) 00132 const TH1D* fBini; // reconstructed distribution (MC) 00133 const TH1D* fXini; // truth distribution (MC) 00134 const TH2D* fAdet; // Detector response matrix 00135 00136 // Evaluation of covariance matrices 00137 TH1D* fToyhisto; //! Toy MC histogram 00138 TH2D* fToymat; //! Toy MC detector response matrix 00139 Bool_t fToyMode; //! Internal switch for covariance matrix propagation 00140 Bool_t fMatToyMode; //! Internal switch for evaluation of statistical uncertainties from response matrix 00141 00142 00143 ClassDef( TSVDUnfold, 0 ) // Data unfolding using Singular Value Decomposition (hep-ph/9509307) 00144 }; 00145 00146 #endif
1.5.1