Page MenuHomeHEPForge
Files F19252253

D93.1759122925.diff
No OneTemporary
Actions

Size
11 KB
Referenced Files
None
Subscribers
None

D93.1759122925.diff
View Options

diff --git a/examples/GenFitDpipi.cc b/examples/GenFitDpipi.cc
--- a/examples/GenFitDpipi.cc
+++ b/examples/GenFitDpipi.cc
@@ -176,6 +176,24 @@
LauParameter* nSigEvents = new LauParameter(sigEventsName,nSig,-2.0*nSig,2.0*nSig,kTRUE);
fitModel->setNSigEvents(nSigEvents);
+ // Test 2D Gaussian constraint
+ // Names of LauParameters in the fit
+ std::vector<TString> params;
+ params.push_back("A0_A");
+ params.push_back("A0_Delta");
+ // Mean values for the constraints
+ TVectorD means(2);
+ means[0] = 0.55;
+ means[1] = 1.26;
+ // Covariance matrix (not correlation!)
+ TMatrixD covMat(2,2);
+ covMat(0,0) = 0.0001345;
+ covMat(1,0) = 0.0001183;
+ covMat(0,1) = 0.0001183;
+ covMat(1,1) = 0.001078;
+
+ fitModel->addNDConstraint(params, means, covMat);
+
// Set the number of experiments to generate or fit and which
// experiment to start with
fitModel->setNExpts( nExpt, firstExpt );
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,8 @@
#include "TMatrixD.h"
#include "TObject.h"
#include "TString.h"
+#include "TVectorT.h"
+#include "TSystem.h"
#include "LauAbsFitter.hh"
@@ -135,6 +137,20 @@
*/
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
+ */
+ virtual void checkRepetition(const std::vector<TString>& names);
+
protected:
//! Constructor
LauFitObject();
@@ -154,6 +170,22 @@
//! 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_;
+ };
//! Const access to the constraints store
const std::vector<StoreConstraints>& constraintsStore() const {return storeCon_;}
@@ -161,6 +193,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 +271,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 +320,3 @@
};
#endif
-
diff --git a/src/LauAbsFitModel.cc b/src/LauAbsFitModel.cc
--- a/src/LauAbsFitModel.cc
+++ b/src/LauAbsFitModel.cc
@@ -36,6 +36,7 @@
#include "TSocket.h"
#include "TSystem.h"
#include "TVirtualFitter.h"
+#include "TVectorT.h"
#include "LauAbsFitModel.hh"
#include "LauAbsFitter.hh"
@@ -744,7 +745,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,13 +758,30 @@
{
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();
+ if ( ! conVars_.empty() ){
+ 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();
+
+ Double_t term = ( val - mean )*( val - mean );
+ penalty += term/( 2*width*width );
+ }
+ }
+
+ const std::vector<StoreNDConstraints>& storeNDCon = this->NDConstraintsStore();
+ if ( ! storeNDCon.empty() ){
+ for ( ULong_t i = 0; i<storeNDCon.size(); ++i ){
+ TVectorD paramValues(storeNDCon[i].conLauPars_.size());
+ LauParameterPList params = storeNDCon[i].conLauPars_;
- Double_t term = ( val - mean )*( val - mean );
- penalty += term/( 2*width*width );
+ for ( ULong_t j = 0; j<params.size(); ++j ){
+ LauParameter* param = params[j];
+ paramValues[j] = param->unblindValue();
+ }
+ TVectorD diff = paramValues-storeNDCon[i].means_;
+ penalty += 0.5*storeNDCon[i].invCovMat_.Similarity(diff);
+ }
}
return penalty;
@@ -914,8 +933,31 @@
}
}
+ // Add n-dimensional constraints
+ std::vector<StoreNDConstraints>& storeNDCon = this->NDConstraintsStore();
+ std::vector<LauParameter*> params;
+ for ( auto& itervec : storeNDCon ){
+ for ( auto& iternames : itervec.conPars_ ){
+ for ( auto& iterfit : fitVars_ ){
+ if ( iternames == iterfit->name() ){
+ // Check parameters do not have a 1D Gaussian constraint applied
+ if ( iterfit->gaussConstraint() ){
+ std::cerr << "ERROR in LauAbsFitModel::addConParameters: parameter in n-dimsensional constraint already has a 1d constraint applied" << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ params.push_back(iterfit);
+ }
+ }
+ }
+ if ( params.size() != itervec.conPars_.size() ){
+ std::cerr << "Error in LauAbsFitModel::addConParameters : Could not match parameter names for n-dimensional constraint" << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ itervec.conLauPars_ = params;
+ }
}
+
void LauAbsFitModel::updateFitParameters(LauPdfList& pdfList)
{
for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
diff --git a/src/LauFitObject.cc b/src/LauFitObject.cc
--- a/src/LauFitObject.cc
+++ b/src/LauFitObject.cc
@@ -26,6 +26,10 @@
\brief File containing implementation of LauFitObject class.
*/
+#include <iostream>
+
+#include "TMatrixDfwd.h"
+#include "TVectorT.h"
#include "LauFitObject.hh"
@@ -90,6 +94,8 @@
void LauFitObject::addConstraint(const TString& formula, const std::vector<TString>& pars, const Double_t mean, const Double_t width)
{
+ this->checkRepetition(pars);
+
StoreConstraints newCon;
newCon.formula_ = formula;
newCon.conPars_ = pars;
@@ -97,5 +103,64 @@
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() != means.GetNrows()) || ( pars.size() != covMat.GetNcols() ) ){
+ std::cerr << "ERROR in LauFitObject::addNDConstraint : Different number of elements in vectors/covariance matrix!" << std::endl;
+ gSystem->Exit( EXIT_FAILURE );
+ }
+
+ this->checkRepetition(pars);
+
+ 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 );
+ }
+
+ StoreNDConstraints newCon { pars, means, invCovMat, {} };
+ storeNDCon_.emplace_back( newCon );
+
+ std::cout << "INFO in LauFitObject::addNDConstraint : Added list of parameters to constrain" << std::endl;
}
+void LauFitObject::checkRepetition( const std::vector<TString>& names )
+{
+ if ( storeCon_.size()==0 && storeNDCon_.size()==0 ) {
+ return;
+ }
+
+ //Check in formula constraints
+ for ( auto& itervec : storeCon_ ){
+ for ( auto& itername : itervec.conPars_ ){
+ for ( auto& newname : names ){
+ if ( itername == newname ){
+ std::cerr << "ERROR in LauFitObject::checkRepetition: named parameter " << newname << " already used in a constraint" << std::endl;
+ gSystem->Exit( EXIT_FAILURE );
+ }
+ }
+ }
+ }
+ //Check in ND constraints
+ for ( auto& itervec : storeNDCon_ ){
+ for ( auto& itername : itervec.conPars_ ){
+ for ( auto& newname : names ){
+ if ( itername == newname ){
+ std::cerr << "ERROR in LauFitObject::checkRepetition: named parameter " << newname << " already used in a constraint" << std::endl;
+ gSystem->Exit( EXIT_FAILURE );
+ }
+ }
+ }
+ }
+}
diff --git a/src/LauSimFitCoordinator.cc b/src/LauSimFitCoordinator.cc
--- a/src/LauSimFitCoordinator.cc
+++ b/src/LauSimFitCoordinator.cc
@@ -785,6 +785,21 @@
penalty += term/( 2*width*width );
}
+ const std::vector<StoreNDConstraints>& storeNDCon = this->NDConstraintsStore();
+ if ( ! storeNDCon.empty() ){
+ for ( ULong_t i = 0; i<storeNDCon.size(); ++i ){
+ TVectorD paramValues(storeNDCon[i].conLauPars_.size());
+ std::vector<LauParameter*> params = storeNDCon[i].conLauPars_;
+
+ for ( ULong_t j = 0; j<params.size(); ++j ){
+ LauParameter* param = params[j];
+ paramValues[j] = param->unblindValue();
+ }
+ TVectorD diff = paramValues-storeNDCon[i].means_;
+ penalty += 0.5*storeNDCon[i].invCovMat_.Similarity(diff);
+ }
+ }
+
return penalty;
}
@@ -829,7 +844,29 @@
std::cout << " : Parameter: " << (*iterparam)->name() << std::endl;
}
}
-
+
+ // Add n-dimensional constraints
+ std::vector<StoreNDConstraints>& storeNDCon = this->NDConstraintsStore();
+ std::vector<LauParameter*> params;
+ for ( auto& itervec : storeNDCon ){
+ for ( auto& iternames : itervec.conPars_ ){
+ for ( auto& iterfit : params_ ){
+ if ( iternames == iterfit->name() ){
+ // Check parameters do not have a 1D Gaussian constraint applied
+ if ( iterfit->gaussConstraint() ){
+ std::cerr << "ERROR in LauAbsFitModel::addConParameters: parameter in n-dimsensional constraint already has a 1d constraint applied" << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ params.push_back(iterfit);
+ }
+ }
+ }
+ if ( params.size() != itervec.conPars_.size() ){
+ std::cerr << "Error in LauAbsFitModel::addConParameters : Could not match parameter names for n-dimensional constraint" << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ itervec.conLauPars_ = params;
+ }
}
Bool_t LauSimFitCoordinator::finalise()

File Metadata

Mime Type
text/plain
Expires
Mon, Sep 29, 6:15 AM (10 h, 8 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6557007
Default Alt Text
D93.1759122925.diff (11 KB)

Event Timeline

HEPForge.org · .