00001 // Author: Stefan Schmitt 00002 // DESY, 23/01/09 00003 00004 // Version 17.5, bug fixes in TUnfold fix problem with GetEmatrixSysUncorr 00005 // 00006 // History: 00007 // Version 17.4, in parallel to changes in TUnfoldBinning 00008 // Version 17.3, in parallel to changes in TUnfoldBinning 00009 // Version 17.2, add methods to find back systematic and background sources 00010 // Version 17.1, bug fix with background uncertainty 00011 // Version 17.0, possibility to specify an error matrix with SetInput 00012 // Version 16.2, bug-fix with the calculation of background errors 00013 // Version 16.1, parallel to changes in TUnfold 00014 // Version 16.0, parallel to changes in TUnfold 00015 // Version 15, fix bugs with uncorr. uncertainties, add backgnd subtraction 00016 // Version 14, with changes in TUnfoldSys.cxx 00017 // Version 13, support for systematic errors 00018 00019 #ifndef ROOT_TUnfoldSys 00020 #define ROOT_TUnfoldSys 00021 00022 ////////////////////////////////////////////////////////////////////////// 00023 // // 00024 // // 00025 // TUnfoldSys, an extension of the class TUnfold to correct for // 00026 // migration effects. It provides methods for background subtraction // 00027 // and propagation of systematic uncertainties // 00028 // // 00029 // Citation: S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201] // 00030 // // 00031 ////////////////////////////////////////////////////////////////////////// 00032 00033 /* 00034 This file is part of TUnfold. 00035 00036 TUnfold is free software: you can redistribute it and/or modify 00037 it under the terms of the GNU General Public License as published by 00038 the Free Software Foundation, either version 3 of the License, or 00039 (at your option) any later version. 00040 00041 TUnfold is distributed in the hope that it will be useful, 00042 but WITHOUT ANY WARRANTY; without even the implied warranty of 00043 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00044 GNU General Public License for more details. 00045 00046 You should have received a copy of the GNU General Public License 00047 along with TUnfold. If not, see <http://www.gnu.org/licenses/>. 00048 */ 00049 00050 #include <TMap.h> 00051 #include <TSortedList.h> 00052 #include "TUnfold.h" 00053 00054 #define TUnfoldSys TUnfoldSysV17 00055 00056 class TUnfoldSysV17 : public TUnfoldV17 { 00057 private: 00058 void InitTUnfoldSys(void); // initialize all data members 00059 protected: 00060 /// Input: normalized errors from input matrix 00061 TMatrixDSparse *fDAinRelSq; 00062 /// Input: normalized column err.sq. (inp.matr.) 00063 TMatrixD* fDAinColRelSq; 00064 /// Input: underflow/overflow bins 00065 TMatrixD* fAoutside; 00066 /// Input: correlated errors 00067 TMap *fSysIn; 00068 /// Input: size of background sources 00069 TMap *fBgrIn; 00070 /// Input: uncorr error squared from bgr sources 00071 TMap *fBgrErrUncorrInSq; 00072 /// Input: background sources correlated error 00073 TMap *fBgrErrScaleIn; 00074 /// Input: error on tau 00075 Double_t fDtau; 00076 /// Input: fY prior to bgr subtraction 00077 TMatrixD *fYData; 00078 /// Input: error on fY prior to bgr subtraction 00079 TMatrixDSparse *fVyyData; 00080 /// Result: syst.error from fDA2 on fX 00081 TMatrixDSparse *fEmatUncorrX; 00082 /// Result: syst.error from fDA2 on fAx 00083 TMatrixDSparse *fEmatUncorrAx; 00084 /// Result: syst.shift from fSysIn on fX 00085 TMap *fDeltaCorrX; 00086 /// Result: syst.shift from fSysIn on fAx 00087 TMap *fDeltaCorrAx; 00088 /// Result: systematic shift from tau 00089 TMatrixDSparse *fDeltaSysTau; 00090 protected: 00091 virtual void ClearResults(void); // clear all results 00092 virtual void PrepareSysError(void); // common calculations for syst.errors 00093 virtual TMatrixDSparse *PrepareUncorrEmat(const TMatrixDSparse *m1,const TMatrixDSparse *m2); // calculate uncorrelated error matrix 00094 virtual TMatrixDSparse *PrepareCorrEmat(const TMatrixDSparse *m1,const TMatrixDSparse *m2,const TMatrixDSparse *dsys); // calculate correlated error matrix 00095 void ScaleColumnsByVector(TMatrixDSparse *m,const TMatrixTBase<Double_t> *v) const; // scale columns of m by the corresponding rows of v 00096 void VectorMapToHist(TH1 *hist_delta,const TMatrixDSparse *delta,const Int_t *binMap); // map and sum vector delta, save in hist_delta 00097 void GetEmatrixFromVyy(const TMatrixDSparse *vyy,TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat); // propagate error matrix vyy to the result 00098 void DoBackgroundSubtraction(void); 00099 TMatrixDSparse *GetSummedErrorMatrixYY(void); 00100 TMatrixDSparse *GetSummedErrorMatrixXX(void); 00101 public: 00102 /// type of matrix specified with AddSysError() 00103 enum ESysErrMode { 00104 /// matrix is an alternative to the default matrix, the errors are the difference to the original matrix 00105 kSysErrModeMatrix=0, 00106 /// matrix gives the absolute shifts 00107 kSysErrModeShift=1, 00108 /// matrix gives the relative shifts 00109 kSysErrModeRelative=2 00110 }; 00111 TUnfoldSysV17(const TH2 *hist_A, EHistMap histmap, ERegMode regmode = kRegModeSize, 00112 EConstraint constraint=kEConstraintArea); // constructor 00113 TUnfoldSysV17(void); // for derived classes 00114 virtual ~ TUnfoldSysV17(void); // delete data members 00115 void AddSysError(const TH2 *sysError,const char *name, EHistMap histmap, 00116 ESysErrMode mode); // add a systematic error source 00117 void SubtractBackground(const TH1 *hist_bgr,const char *name, 00118 Double_t scale=1.0, 00119 Double_t scale_error=0.0); // subtract background prior to unfolding 00120 virtual Int_t SetInput(const TH1 *hist_y,Double_t scaleBias=0.0,Double_t oneOverZeroError=0.0,const TH2 *hist_vyy=0,const TH2 *hist_vyy_inv=0); // define input consistently in case of background subtraction 00121 void SetTauError(Double_t delta_tau); // set uncertainty on tau 00122 TSortedList *GetBgrSources(void) const; // get names of background sources 00123 TSortedList *GetSysSources(void) const; // get names of systematic sources 00124 void GetBackground(TH1 *bgr,const char *bgrSource=0,const Int_t *binMap=0,Int_t includeError=3,Bool_t clearHist=kTRUE) const; // get background as histogram 00125 void GetEmatrixSysBackgroundUncorr(TH2 *ematrix,const char *source, 00126 const Int_t *binMap=0,Bool_t clearEmat=kTRUE); // get error matrix from uncorrelated error of one background source 00127 void GetEmatrixSysBackgroundScale(TH2 *ematrix,const char *source, 00128 const Int_t *binMap=0,Bool_t clearEmat=kTRUE); // get error matrix from the scale error of one background source 00129 Bool_t GetDeltaSysBackgroundScale(TH1 *delta,const char *source, 00130 const Int_t *binMap=0); // get correlated uncertainty induced by the scale uncertainty of a background source 00131 void GetEmatrixSysUncorr(TH2 *ematrix,const Int_t *binMap=0,Bool_t clearEmat=kTRUE); // get error matrix contribution from uncorrelated errors on the matrix A 00132 void GetEmatrixSysSource(TH2 *ematrix,const char *source, 00133 const Int_t *binMap=0,Bool_t clearEmat=kTRUE); // get error matrix from one systematic source 00134 Bool_t GetDeltaSysSource(TH1 *hist_delta,const char *source, 00135 const Int_t *binMap=0); // get systematic shifts from one systematic source 00136 void GetEmatrixSysTau(TH2 *ematrix, 00137 const Int_t *binMap=0,Bool_t clearEmat=kTRUE); // get error matrix from tau variation 00138 Bool_t GetDeltaSysTau(TH1 *delta,const Int_t *binMap=0); // get correlated uncertainty from varying tau 00139 void GetEmatrixInput(TH2 *ematrix,const Int_t *binMap=0,Bool_t clearEmat=kTRUE); // get error contribution from input vector 00140 void GetEmatrixTotal(TH2 *ematrix,const Int_t *binMap=0); // get total error including systematic,statistical,background,tau errors 00141 void GetRhoItotal(TH1 *rhoi,const Int_t *binMap=0,TH2 *invEmat=0); // get global correlation coefficients including systematic,statistical,background,tau errors 00142 Double_t GetChi2Sys(void); // get total chi**2 including all systematic errors 00143 ClassDef(TUnfoldSysV17, TUnfold_CLASS_VERSION) //Unfolding with support for systematic error propagation 00144 }; 00145 00146 #endif