Page MenuHomeHEPForge
Files F11221649

D93.id400.diff
No OneTemporary
Actions

D93.id400.diff
View Options

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<TString>& 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<TString>& 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<TString>& 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<TString> 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<LauParameter*> conLauPars_;
+ //! The values of the LauParameters
+ TVectorD values_;
+ };
//! Const access to the constraints store
const std::vector<StoreConstraints>& constraintsStore() const {return storeCon_;}
@@ -161,6 +195,12 @@
//! Access to the constraints store
std::vector<StoreConstraints>& constraintsStore() {return storeCon_;}
+ //! Const access to the ND constraints store
+ const std::vector<StoreNDConstraints>& NDConstraintsStore() const {return storeNDCon_;}
+
+ //! Access to the ND constraints store
+ std::vector<StoreNDConstraints>& 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<StoreConstraints> storeCon_;
+ //! Store the ND constraints for fit parameters until initialisation is complete
+ std::vector<StoreNDConstraints> 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<StoreNDConstraints>& 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<StoreNDConstraints>& storeNDCon = this->NDConstraintsStore();
+ for ( ULong_t i = 0; i<storeNDCon.size(); ++i ){
+ LauParameterPList& params = storeNDCon[i].conLauPars_;
+
+ for ( ULong_t j = 0; j<params.size(); ++j ){
+ LauParameter* param = params[j];
+ storeNDCon[i].values_[j] = param->unblindValue();
+ }
+ 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<StoreNDConstraints>& 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 <iostream>
+
+#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<TString>& 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<TString>& 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<std::size_t>(means.GetNrows()) ) || ( pars.size() != static_cast<std::size_t>(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<TString>& 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<LauAbsRValue*>::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<StoreNDConstraints>& storeNDCon = this->NDConstraintsStore();
+ for ( ULong_t i = 0; i<storeNDCon.size(); ++i ){
+ std::vector<LauParameter*>& params = storeNDCon[i].conLauPars_;
+
+ for ( ULong_t j = 0; j<params.size(); ++j ){
+ LauParameter* param = params[j];
+ storeNDCon[i].values_[j] = param->unblindValue();
+ }
+ 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<StoreNDConstraints>& 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()

File Metadata

Mime Type
text/plain
Expires
Tue, May 13, 10:41 AM (9 h, 5 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5088214
Default Alt Text
D93.id400.diff (11 KB)

Event Timeline