diff --git a/doc/ReleaseNotes.md b/doc/ReleaseNotes.md --- a/doc/ReleaseNotes.md +++ b/doc/ReleaseNotes.md @@ -1,5 +1,9 @@ # Laura++ release notes +9th June 2023 Mark Whitehead +* Added functionality to include n-dimensional Gaussian constraints + - see https://phab.hepforge.org/T214 + 21st November 2022 Andy Morris * Use Minuit to automatically determine ASq max for signal in LauIsobarDynamics - see https://phab.hepforge.org/D85 diff --git a/inc/LauAbsFitModel.hh b/inc/LauAbsFitModel.hh --- a/inc/LauAbsFitModel.hh +++ b/inc/LauAbsFitModel.hh @@ -768,7 +768,7 @@ //! Extra variables that aren't in the fit but are stored in the ntuple LauParameterList extraVars_; - //! Internal vectors of Gaussian parameters + //! Internal vectors of Gaussian parameters LauAbsRValuePList conVars_; // Input data and output ntuple diff --git a/inc/LauFitObject.hh b/inc/LauFitObject.hh --- a/inc/LauFitObject.hh +++ b/inc/LauFitObject.hh @@ -34,6 +34,7 @@ #include "TMatrixD.h" #include "TObject.h" #include "TString.h" +#include "TVectorD.h" #include "LauAbsFitter.hh" @@ -135,6 +136,21 @@ */ virtual void addConstraint(const TString& formula, const std::vector& pars, const Double_t mean, const Double_t width); + //! Store n-dimensional constraint information for fit parameters + /*! + \param [in] pars a vector of LauParameter names to be used in the constraint + \param [in] means the values of the means of the Gaussian constraint + \param [in] covMat the covariance matrix of the parameters of the Gaussian constraint + */ + virtual void addNDConstraint(const std::vector& pars, const TVectorD& means, const TMatrixD& covMat); + + //! Check if parameters names for constraints have already been used elsewhere + /*! + \param [in] names a vector of parameter names + \return kTRUE if no repetitions found, kFALSE if one or more repetitions found + */ + virtual Bool_t checkRepetition(const std::vector& names); + protected: //! Constructor LauFitObject(); @@ -154,6 +170,24 @@ //! The width of the Gaussian constraint to be applied Double_t width_; }; + + // Setup a struct to store information on n-dimensional constrained fit parameters + /*! + \struct StoreNDConstraints + \brief Struct to store n-dimensional constraint information until the fit is run + */ + struct StoreNDConstraints { + //! The list of LauParameter names to be used in the constraint + std::vector conPars_; + //! The mean value of the Gaussian constraints to be applied + TVectorD means_; + //! The inverse covariance matrix of the parameters of the Gaussian constraint to be applied + TMatrixD invCovMat_; + //! The LauParameters used in the constraints + std::vector conLauPars_; + //! The values of the LauParameters + TVectorD values_; + }; //! Const access to the constraints store const std::vector& constraintsStore() const {return storeCon_;} @@ -161,6 +195,12 @@ //! Access to the constraints store std::vector& constraintsStore() {return storeCon_;} + //! Const access to the ND constraints store + const std::vector& NDConstraintsStore() const {return storeNDCon_;} + + //! Access to the ND constraints store + std::vector& NDConstraintsStore() {return storeNDCon_;} + //! Reset the good/bad fit counters void resetFitCounters(); @@ -233,6 +273,9 @@ //! Store the constraints for fit parameters until initialisation is complete std::vector storeCon_; + //! Store the ND constraints for fit parameters until initialisation is complete + std::vector storeNDCon_; + //! Option to perform a two stage fit Bool_t twoStageFit_; @@ -279,4 +322,3 @@ }; #endif - diff --git a/src/LauAbsFitModel.cc b/src/LauAbsFitModel.cc --- a/src/LauAbsFitModel.cc +++ b/src/LauAbsFitModel.cc @@ -744,7 +744,8 @@ } // Calculate any penalty terms from Gaussian constrained variables - if ( ! conVars_.empty() ){ + const std::vector& storeNDCon = this->NDConstraintsStore(); + if ( ! conVars_.empty() || ! storeNDCon.empty() ){ logLike -= this->getLogLikelihoodPenalty(); } @@ -756,15 +757,27 @@ { Double_t penalty(0.0); - for ( LauAbsRValuePList::const_iterator iter = conVars_.begin(); iter != conVars_.end(); ++iter ) { - Double_t val = (*iter)->unblindValue(); - Double_t mean = (*iter)->constraintMean(); - Double_t width = (*iter)->constraintWidth(); + for ( LauAbsRValue* par : conVars_ ) { + Double_t val = par->unblindValue(); + Double_t mean = par->constraintMean(); + Double_t width = par->constraintWidth(); Double_t term = ( val - mean )*( val - mean ); penalty += term/( 2*width*width ); } + std::vector& storeNDCon = this->NDConstraintsStore(); + for ( ULong_t i = 0; iunblindValue(); + } + TVectorD diff = storeNDCon[i].values_-storeNDCon[i].means_; + penalty += 0.5*storeNDCon[i].invCovMat_.Similarity(diff); + } + return penalty; } @@ -914,6 +927,26 @@ } } + // Add n-dimensional constraints + std::vector& storeNDCon = this->NDConstraintsStore(); + for ( auto& constraint : storeNDCon ){ + for ( auto& parname : constraint.conPars_ ){ + for ( auto& fitPar : fitVars_ ){ + if ( parname == fitPar->name() ){ + // Check parameters do not have a 1D Gaussian constraint applied + if ( fitPar->gaussConstraint() ){ + std::cerr << "ERROR in LauAbsFitModel::addConParameters: parameter in n-dimensional constraint already has a 1d constraint applied" << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + constraint.conLauPars_.push_back(fitPar); + } + } + } + if ( constraint.conLauPars_.size() != constraint.conPars_.size() ){ + std::cerr << "Error in LauAbsFitModel::addConParameters : Could not match parameter names for n-dimensional constraint" << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + } } void LauAbsFitModel::updateFitParameters(LauPdfList& pdfList) diff --git a/src/LauFitObject.cc b/src/LauFitObject.cc --- a/src/LauFitObject.cc +++ b/src/LauFitObject.cc @@ -26,6 +26,11 @@ \brief File containing implementation of LauFitObject class. */ +#include + +#include "TMatrixD.h" +#include "TSystem.h" +#include "TVectorD.h" #include "LauFitObject.hh" @@ -90,6 +95,11 @@ void LauFitObject::addConstraint(const TString& formula, const std::vector& pars, const Double_t mean, const Double_t width) { + if ( ! this->checkRepetition(pars) ){ + std::cerr << "ERROR in LauFitObject::addConstraint : Parameter(s) added to multiple constraints!" << std::endl; + gSystem->Exit( EXIT_FAILURE ); + } + StoreConstraints newCon; newCon.formula_ = formula; newCon.conPars_ = pars; @@ -97,5 +107,70 @@ newCon.width_ = width; storeCon_.push_back(newCon); + + std::cout << "INFO in LauFitObject::addConstraint : Added formula to constrain" << std::endl; +} + +void LauFitObject::addNDConstraint( const std::vector& pars, const TVectorD& means, const TMatrixD& covMat) +{ + if ( covMat.GetNcols() != covMat.GetNrows() ){ + std::cerr << "ERROR in LauFitObject::addNDConstraint : Covariance matrix is not square!" << std::endl; + gSystem->Exit( EXIT_FAILURE ); + } + + if ( ( pars.size() != static_cast(means.GetNrows()) ) || ( pars.size() != static_cast(covMat.GetNcols()) ) ){ + std::cerr << "ERROR in LauFitObject::addNDConstraint : Different number of elements in vectors/covariance matrix!" << std::endl; + gSystem->Exit( EXIT_FAILURE ); + } + + if ( ! this->checkRepetition(pars) ){ + std::cerr << "ERROR in LauFitObject::addNDConstraint : Parameter(s) added to multiple constraints!" << std::endl; + gSystem->Exit( EXIT_FAILURE ); + } + + TMatrixD invCovMat {TMatrixD::kInverted, covMat}; + // Check invertion was successful + if ( invCovMat == covMat ){ + std::cerr << "ERROR in LauFitObject::addNDConstraint : covariance matrix inversion failed, check your input!" << std::endl; + gSystem->Exit( EXIT_FAILURE ); + } + + TVectorD values(means.GetNrows()); + StoreNDConstraints newCon { pars, means, invCovMat, {}, values }; + storeNDCon_.emplace_back( newCon ); + + std::cout << "INFO in LauFitObject::addNDConstraint : Added list of parameters to constrain" << std::endl; } +Bool_t LauFitObject::checkRepetition( const std::vector& names ) +{ + Bool_t allOK(kTRUE); + + if ( storeCon_.size()==0 && storeNDCon_.size()==0 ) { + return allOK; + } + + //Check in formula constraints + for ( auto& constraint : storeCon_ ){ + for ( auto& parname : constraint.conPars_ ){ + for ( auto& newname : names ){ + if ( parname == newname ){ + std::cerr << "WARNING in LauFitObject::checkRepetition: named parameter " << newname << " already used in a constraint" << std::endl; + allOK = kFALSE; + } + } + } + } + //Check in ND constraints + for ( auto& constraint : storeNDCon_ ){ + for ( auto& parname : constraint.conPars_ ){ + for ( auto& newname : names ){ + if ( parname == newname ){ + std::cerr << "WARNING in LauFitObject::checkRepetition: named parameter " << newname << " already used in a constraint" << std::endl; + allOK = kFALSE; + } + } + } + } + return allOK; +} diff --git a/src/LauSimFitCoordinator.cc b/src/LauSimFitCoordinator.cc --- a/src/LauSimFitCoordinator.cc +++ b/src/LauSimFitCoordinator.cc @@ -776,15 +776,27 @@ { Double_t penalty(0.0); - for ( std::vector::const_iterator iter = conVars_.begin(); iter != conVars_.end(); ++iter ) { - Double_t val = (*iter)->unblindValue(); - Double_t mean = (*iter)->constraintMean(); - Double_t width = (*iter)->constraintWidth(); + for ( LauAbsRValue* par : conVars_ ) { + Double_t val = par->unblindValue(); + Double_t mean = par->constraintMean(); + Double_t width = par->constraintWidth(); Double_t term = ( val - mean )*( val - mean ); penalty += term/( 2*width*width ); } + std::vector& storeNDCon = this->NDConstraintsStore(); + for ( ULong_t i = 0; i& params = storeNDCon[i].conLauPars_; + + for ( ULong_t j = 0; junblindValue(); + } + TVectorD diff = storeNDCon[i].values_-storeNDCon[i].means_; + penalty += 0.5*storeNDCon[i].invCovMat_.Similarity(diff); + } + return penalty; } @@ -829,7 +841,27 @@ std::cout << " : Parameter: " << (*iterparam)->name() << std::endl; } } - + + // Add n-dimensional constraints + std::vector& storeNDCon = this->NDConstraintsStore(); + for ( auto& constraint : storeNDCon ){ + for ( auto& parname : constraint.conPars_ ){ + for ( auto& fitPar : params_ ){ + if ( parname == fitPar->name() ){ + // Check parameters do not have a 1D Gaussian constraint applied + if ( fitPar->gaussConstraint() ){ + std::cerr << "ERROR in LauAbsFitModel::addConParameters: parameter in n-dimensional constraint already has a 1d constraint applied" << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + constraint.conLauPars_.push_back(fitPar); + } + } + } + if ( constraint.conLauPars_.size() != constraint.conPars_.size() ){ + std::cerr << "Error in LauAbsFitModel::addConParameters : Could not match parameter names for n-dimensional constraint" << std::endl; + gSystem->Exit(EXIT_FAILURE); + } + } } Bool_t LauSimFitCoordinator::finalise()