Page MenuHomeHEPForge

D93.1759178818.diff
No OneTemporary

Size
7 KB
Referenced Files
None
Subscribers
None

D93.1759178818.diff

diff --git a/examples/GenFitDpipi.cc b/examples/GenFitDpipi.cc
--- a/examples/GenFitDpipi.cc
+++ b/examples/GenFitDpipi.cc
@@ -133,8 +133,10 @@
mipw->floatKnotsSecondStage(kFALSE);
//Set magnitude and phase for the knots (including the end points here)
- mipw->setKnotAmp(0, 0.12, -2.82,kFALSE,kFALSE);
- mipw->setKnotAmp(1, 0.58, -1.56,kFALSE,kFALSE);
+ //mipw->setKnotAmp(0, 0.12, -2.82,kFALSE,kFALSE);
+ //mipw->setKnotAmp(1, 0.58, -1.56,kFALSE,kFALSE);
+ mipw->setKnotAmp(0, 0.5, -2.82,kFALSE,kFALSE);
+ mipw->setKnotAmp(1, 1.0, -1.56,kFALSE,kFALSE);
mipw->setKnotAmp(2, 0.73, -1.00,kFALSE,kFALSE);
mipw->setKnotAmp(3, 0.68, -0.42,kFALSE,kFALSE);
mipw->setKnotAmp(4, 0.5, 0.0,kTRUE,kTRUE);
@@ -176,6 +178,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("Dst0_0_A0");
+ params.push_back("Dst0_0_A1");
+ // Mean values for the constraints
+ std::vector<Double_t> means;
+ means.push_back(0.12);
+ means.push_back(0.58);
+ // Covariance matrix (not correlation!)
+ TMatrixD covMat(2,2);
+ covMat(0,0) = 0.0001736;
+ covMat(1,0) = 0.0000755;
+ covMat(0,1) = 0.0000755;
+ covMat(1,1) = 0.0001104;
+
+ 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
@@ -294,6 +294,14 @@
*/
void setParametersFileFallback(const TString& fileName, const TString& treeName, const std::map<TString, Double_t>& parameters, const Bool_t fix);
+ //! Add an n-dimensional Gaussian constraint for a set of parameters including correlations.
+ /*!
+ \param [in] paramNames vector of LauParameter names to constrain
+ \param [in] meanVals vector of corresponding mean values to use in the constraint
+ \param [in] covMat the covariance matrix relating the set of parameters in the constraint
+ */
+ void addnDConstraint(std::vector<TString> paramNames, std::vector<Double_t> meanVals, TMatrixD covMat);
+
protected:
// Some typedefs
@@ -768,9 +776,21 @@
//! 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_;
+ //! Interal vector of n-dimensional Gaussian constraint parameters
+ std::vector<LauParameterPList> nDConParams_;
+
+ //! Interal vector of n-dimensional Gaussian constraint inverse covariance matrices
+ std::vector<TMatrixD> nDConInvCovMats_;
+
+ //! Interal vector of n-dimensional Gaussian constraint central values
+ std::vector<std::vector<Double_t>> nDConMuVals_;
+
+ //! Interal vector of n-dimensional Gaussian constraint parameter names
+ std::vector<std::vector<TString>> nDConNames_;
+
// Input data and output ntuple
//! The input data
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,7 @@
}
// Calculate any penalty terms from Gaussian constrained variables
- if ( ! conVars_.empty() ){
+ if ( ! conVars_.empty() || ! nDConMuVals_.empty() ){
logLike -= this->getLogLikelihoodPenalty();
}
@@ -756,13 +757,34 @@
{
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 );
+ Double_t term = ( val - mean )*( val - mean );
+ penalty += term/( 2*width*width );
+ }
+ }
+
+ if ( ! nDConMuVals_.empty() ){
+
+ for ( ULong_t i = 0; i<nDConMuVals_.size(); ++i ){
+ TVectorD paramValues(nDConParams_[i].size());
+ TVectorD meanValues(nDConMuVals_[i].size());
+ LauParameterPList params = nDConParams_[i];
+ std::vector<Double_t> means = nDConMuVals_[i];
+
+ for ( ULong_t j = 0; j<means.size(); ++j ){
+ LauParameter* param = params[j];
+ paramValues[j] = param->unblindValue();
+ meanValues[j] = means[j];
+ }
+
+ TVectorD diff = paramValues-meanValues;
+ penalty += 0.5*nDConInvCovMats_[i].Similarity(diff);
+ }
}
return penalty;
@@ -914,6 +936,62 @@
}
}
+ // Add n-dimensional constraints
+ std::vector<LauParameter*> params;
+ for ( std::vector<std::vector<TString>>::const_iterator itervec = nDConNames_.begin(); itervec != nDConNames_.end(); ++itervec ){
+ for ( std::vector<TString>::const_iterator iternames = itervec->begin(); iternames != itervec->end(); ++iternames ) {
+ for ( LauParameterPList::const_iterator iterfit = fitVars_.begin(); iterfit != fitVars_.end(); ++iterfit ) {
+ 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->size() ){
+ std::cerr << "Error in LauAbsFitModel::addConParameters : Could not match parameter names for n-dimensional constraint" << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ nDConParams_.push_back(params);
+ }
+
+
+}
+
+void LauAbsFitModel::addnDConstraint(std::vector<TString> paramNames, std::vector<Double_t> meanVals, TMatrixD covMat) {
+
+ if ( covMat.GetNcols() != covMat.GetNrows() ){
+ std::cerr << "ERROR in LauAbsFitModel::addnDConstraint : Covariance matrix is not square!" << std::endl;
+ gSystem->Exit( EXIT_FAILURE );
+ }
+
+ if ( (paramNames.size() != meanVals.size()) || ( paramNames.size() != covMat.GetNcols() ) ){
+ std::cerr << "ERROR in LauAbsFitModel::addnDConstraint : Different number of elements in vectors/covariance matrix!" << std::endl;
+ gSystem->Exit( EXIT_FAILURE );
+ }
+
+ // Check if a parameter name has already been used in another contraint
+ if ( nDConNames_.size() != 0 ){
+ for ( std::vector<std::vector<TString>>::const_iterator itervec = nDConNames_.begin(); itervec != nDConNames_.end(); ++itervec ){
+ for ( std::vector<TString>::const_iterator itername = itervec->begin(); itername != itervec->end(); ++itername ) {
+ for ( std::vector<TString>::const_iterator newname = paramNames.begin(); newname != paramNames.end(); ++newname ){
+ if ( (*itername) == (*newname) ){
+ std::cerr << "ERROR in LauAbsFitModel::addnDConstraint : named parameter " << (*newname) << " already used in n-dimensional constraint" << std::endl;
+ gSystem->Exit( EXIT_FAILURE );
+ }
+ }
+ }
+ }
+ }
+
+ nDConNames_.push_back(paramNames);
+ nDConMuVals_.push_back(meanVals);
+ nDConInvCovMats_.push_back(covMat.Invert());
+
+ std::cout << "INFO in LauAbsFitModel::addnDConstraint : Added list of parameters to constrain" << std::endl;
}
void LauAbsFitModel::updateFitParameters(LauPdfList& pdfList)

File Metadata

Mime Type
text/plain
Expires
Mon, Sep 29, 9:46 PM (22 h, 34 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6571702
Default Alt Text
D93.1759178818.diff (7 KB)

Event Timeline