diff --git a/CMakeLists.txt b/CMakeLists.txt index e4b80e6d9a62702af711522d8b45ee4692f0493a..c78032a76b02f94ff002d314f85bcc02dcc710f7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,7 +32,7 @@ execute_process( COMMAND ln -sf ${RooUnfoldHeaders} -t ${CMAKE_CURRENT_BINARY_DI include(cmake/PlainROOT.cmake) -if(${RooUnfoldTests} AND NOT RooUnfoldTests) +if(RooUnfoldTests) include(cmake/testing.cmake) endif() diff --git a/src/RooUnfold.h b/src/RooUnfold.h index b99ada2fac131a70e5a5079db02647befe498ed0..7ddfe15e71bad32f6499e03f60d025689199c75c 100644 --- a/src/RooUnfold.h +++ b/src/RooUnfold.h @@ -57,7 +57,7 @@ public: // Special constructors - RooUnfoldT<Hist,Hist2D> (const RooUnfoldResponseT<Hist,Hist2D>* res, const Hist* meas, const char* name= 0, const char* title= 0); + RooUnfoldT(const RooUnfoldResponseT<Hist,Hist2D>* res, const Hist* meas, const char* name= 0, const char* title= 0); // Set up an existing object diff --git a/src/RooUnfoldBayes.cxx b/src/RooUnfoldBayes.cxx index c3fe3e666e9c2aa275de3608fa0427e92b124876..db8dd55df6629e3d5b4fc272dcb571c415838abb 100644 --- a/src/RooUnfoldBayes.cxx +++ b/src/RooUnfoldBayes.cxx @@ -51,6 +51,29 @@ RooUnfoldBayesT<Hist,Hist2D>::RooUnfoldBayesT (const RooUnfoldResponseT<Hist,His Init(); } +template<class Hist,class Hist2D> +void RooUnfoldBayesT<Hist,Hist2D>::SetPriors(const TH1 *priors) +{ + if(priors == nullptr){ + std::cout << "Error SetPriors: no priors provided" << std::endl; + } + Int_t nbins = priors->GetNbinsX(); + if( nbins != _nc){ + std::cout << "Error SetPriors: user's prior size:" << nbins << " different from response matrix size:" << _nc << std::endl; + } + _P0C.Clear(); + _P0C.ResizeTo(nbins); + for(Int_t i = 1; i< nbins+1; i++){ + _P0C[i-1] = priors->GetBinContent(i); + } + _priors = 0; + if(this->verbose() >= 2) { + std::cout << "Priors set to:" << std::endl; + _P0C.Print(); + } +} + + template<class Hist,class Hist2D> RooUnfolding::Algorithm RooUnfoldBayesT<Hist,Hist2D>::GetAlgorithm () const { @@ -168,10 +191,19 @@ RooUnfoldBayesT<Hist,Hist2D>::setup() const if (this->_dosys) this->_dnCidPjk.ResizeTo(this->_nc,this->_ne*this->_nc); // Initial distribution - this->_N0C= this->_nCi.Sum(); - if (this->_N0C!=0.0) { - this->_P0C= this->_nCi; - this->_P0C *= 1.0/this->_N0C; + if(_priors==1){ + this->_N0C= this->_nCi.Sum(); + if (this->_N0C!=0.0) { + this->_P0C= this->_nCi; + this->_P0C *= 1.0/this->_N0C; + } + if(this->verbose() >=0 ) std::cout << "Using response matrix priors" << std::endl; + } else { + if(this->verbose() >= 0) std::cout << "Using user's priors" << std::endl; + } + if(this->verbose()>=1) { + std::cout << "Priors:" << std::endl; + _P0C.Print(); } } diff --git a/src/RooUnfoldBayes.h b/src/RooUnfoldBayes.h index 1f0a28e66c2967c13eea974dde85bb184b76e1de..582b9fd597d1b0a69c8b92e96cf1f2b1d7a53078 100644 --- a/src/RooUnfoldBayes.h +++ b/src/RooUnfoldBayes.h @@ -40,6 +40,7 @@ public: Int_t GetIterations() const; Int_t GetSmoothing() const; const TMatrixD& UnfoldingMatrix() const; + void SetPriors (const TH1* priors); virtual RooUnfolding::Algorithm GetAlgorithm() const override; virtual void SetRegParm (Double_t parm) override; @@ -71,6 +72,7 @@ protected: mutable int _niter; mutable int _smoothit; mutable int _handleFakes; + mutable int _priors =1; mutable int _nc; // number of causes (same as _nt) mutable int _ne; // number of effects (same as _nm) @@ -92,7 +94,7 @@ protected: mutable TMatrixD _dnCidPjk; // response error propagation matrix (stack j,k into each column) public: - ClassDefOverride (RooUnfoldBayesT, 1) // Bayesian Unfolding + ClassDefOverride (RooUnfoldBayesT, 2) // Bayesian Unfolding }; diff --git a/src/RooUnfoldResponse.cxx b/src/RooUnfoldResponse.cxx index 7539289b1dde4b39e48030d8f1c281c9368bfcf4..1314b9a7954032591d319bc8cf62a22b58f2baaa 100644 --- a/src/RooUnfoldResponse.cxx +++ b/src/RooUnfoldResponse.cxx @@ -66,7 +66,7 @@ using std::sqrt; template<class Hist, class Hist2D> class RooUnfoldFoldingFunction { public: - RooUnfoldFoldingFunction<Hist,Hist2D> (const RooUnfoldResponseT<Hist,Hist2D>* res, TF1* func, Double_t eps=1e-12, bool verbose=false) + RooUnfoldFoldingFunction(const RooUnfoldResponseT<Hist,Hist2D>* res, TF1* func, Double_t eps=1e-12, bool verbose=false) : _res(res), _func(func), _eps(eps), _verbose(verbose), _fvals(_res->GetNbinsMeasured()) { _ndim= dynamic_cast<TF3*>(_func) ? 3 : dynamic_cast<TF2*>(_func) ? 2 : 1;