From 7d5c7aec16dd0773c928a7c6a11184194660a996 Mon Sep 17 00:00:00 2001 From: Bertram Kopf <bertram@ep1.rub.de> Date: Tue, 3 Sep 2024 15:00:01 +0200 Subject: [PATCH] added option to start fits with covariance matrix --- AppUtils/AppBase.cc | 56 +++++++++++++++++-- AppUtils/Jamfile | 4 +- FitParams/PwaCovMatrix.cc | 14 ++++- FitParams/PwaCovMatrix.hh | 13 ++++- MinFunctions/MinuitMinimizer.cc | 19 ++++++- MinFunctions/MinuitMinimizer.hh | 7 ++- MinFunctions/PwaFcnServer.cc | 26 ++++----- PwaUtils/GlobalEnv.cc | 17 +++++- PwaUtils/GlobalEnv.hh | 6 +- .../pbarpTopi0pi0eta/pbarpReactionDefault.cfg | 3 + 10 files changed, 138 insertions(+), 27 deletions(-) diff --git a/AppUtils/AppBase.cc b/AppUtils/AppBase.cc index b2ffd8e1..81bd3dee 100644 --- a/AppUtils/AppBase.cc +++ b/AppUtils/AppBase.cc @@ -36,6 +36,7 @@ #include <boost/filesystem/path.hpp> #include "Minuit2/FCNBase.h" #include "Minuit2/FCNGradientBase.h" +#include "Minuit2/MnUserCovariance.h" #include "AppUtils/AppBase.hh" @@ -85,7 +86,7 @@ #include "ggUtils/GGChannelEnv.hh" #include "pipiScatteringUtils/PiPiScatteringChannelEnv.hh" - +using namespace ROOT::Minuit2; AppBase::AppBase() { @@ -805,6 +806,18 @@ void AppBase::fitServerMode(std::shared_ptr<AbsPawianParameters> upar){ } } } + + + std::shared_ptr<MnUserCovariance> theMnUserCov; + if (GlobalEnv::instance()->useCovMatrix()){ + std::shared_ptr<PwaCovMatrix> theCovMatrix=GlobalEnv::instance()->pwaCovMatrix(); + const std::vector<double> dataFromMnCovMat=theCovMatrix->dataFromMnCovMatrix(); + //int nrow=std::sqrt(2.*dataFromMnCovMat.size()+1./4.)-0.5; + unsigned int nrow=theCovMatrix->nRow(); + InfoMsg << "dataFromMnCovMat.size(): " << dataFromMnCovMat.size() << " nrow: " << nrow << endmsg; + theMnUserCov = std::shared_ptr<MnUserCovariance>(new MnUserCovariance(dataFromMnCovMat, nrow)); + } + if(GlobalEnv::instance()->parser()->mode()=="serverQA" || GlobalEnv::instance()->parser()->mode()=="server"){ std::shared_ptr<AbsFcn<FCNBase>> absFcn; std::shared_ptr<NetworkServer> theServer(new NetworkServer(GlobalEnv::instance()->parser()->serverPort(), @@ -825,9 +838,29 @@ void AppBase::fitServerMode(std::shared_ptr<AbsPawianParameters> upar){ absMinimizerPtr->printFitResultQA(evtWeightSumData); return; } - - if(GlobalEnv::instance()->parser()->mode()=="server") - absMinimizerPtr = std::shared_ptr<AbsPawianMinimizer<FCNBase>>(new MinuitMinimizer(absFcn, upar)); + + // std::shared_ptr<MnUserCovariance> theMnUserCov; + // if (GlobalEnv::instance()->useCovMatrix()){ + // std::shared_ptr<PwaCovMatrix> theCovMatrix=GlobalEnv::instance()->pwaCovMatrix(); + // const std::vector<double> dataFromMnCovMat=theCovMatrix->dataFromMnCovMatrix(); + // //int nrow=std::sqrt(2.*dataFromMnCovMat.size()+1./4.)-0.5; + // unsigned int nrow=theCovMatrix->nRow(); + // InfoMsg << "dataFromMnCovMat.size(): " << dataFromMnCovMat.size() << " nrow: " << nrow << endmsg; + // theMnUserCov = std::shared_ptr<MnUserCovariance>(new MnUserCovariance(dataFromMnCovMat, nrow)); + // } + + if(GlobalEnv::instance()->parser()->mode()=="server"){ + if (GlobalEnv::instance()->useCovMatrix()){ + // std::shared_ptr<PwaCovMatrix> theCovMatrix=GlobalEnv::instance()->pwaCovMatrix(); + // const std::vector<double> dataFromMnCovMat=theCovMatrix->dataFromMnCovMatrix(); + // //int nrow=std::sqrt(2.*dataFromMnCovMat.size()+1./4.)-0.5; + // unsigned int nrow=theCovMatrix->nRow(); + // InfoMsg << "dataFromMnCovMat.size(): " << dataFromMnCovMat.size() << " nrow: " << nrow << endmsg; + // std::shared_ptr<MnUserCovariance> theMnUserCov(new MnUserCovariance(dataFromMnCovMat, nrow)); + absMinimizerPtr = std::shared_ptr<AbsPawianMinimizer<FCNBase>>(new MinuitMinimizer<FCNBase>(absFcn, upar, theMnUserCov)); + } + else absMinimizerPtr = std::shared_ptr<AbsPawianMinimizer<FCNBase>>(new MinuitMinimizer<FCNBase>(absFcn, upar)); + } else if (GlobalEnv::instance()->parser()->mode()=="evoserver") absMinimizerPtr = std::shared_ptr<AbsPawianMinimizer<FCNBase>>(new EvoMinimizer(absFcn, upar, @@ -857,8 +890,21 @@ void AppBase::fitServerMode(std::shared_ptr<AbsPawianParameters> upar){ clientNumberWeights())); absFcn=std::shared_ptr<AbsFcn<FCNGradientBase>>(new PwaFcnServer<FCNGradientBase>(theServer)); theServer->WaitForFirstClientLogin(); + + std::shared_ptr<AbsPawianMinimizer<FCNGradientBase>> absMinimizerPtr; + if (GlobalEnv::instance()->useCovMatrix()){ + // std::shared_ptr<PwaCovMatrix> theCovMatrix=GlobalEnv::instance()->pwaCovMatrix(); + // const std::vector<double> dataFromMnCovMat=theCovMatrix->dataFromMnCovMatrix(); + // //int nrow=std::sqrt(2.*dataFromMnCovMat.size()+1./4.)-0.5; + // unsigned int nrow=theCovMatrix->nRow(); + // InfoMsg << "dataFromMnCovMat.size(): " << dataFromMnCovMat.size() << " nrow: " << nrow << endmsg; + // std::shared_ptr<MnUserCovariance> theMnUserCov(new MnUserCovariance(dataFromMnCovMat, nrow)); + absMinimizerPtr = std::shared_ptr<AbsPawianMinimizer<FCNGradientBase>>(new MinuitMinimizer<FCNGradientBase>(absFcn, upar, theMnUserCov)); + } + else absMinimizerPtr = std::shared_ptr<AbsPawianMinimizer<FCNGradientBase>>(new MinuitMinimizer<FCNGradientBase>(absFcn, upar)); + - std::shared_ptr<AbsPawianMinimizer<FCNGradientBase>> absMinimizerPtr(new MinuitMinimizer<FCNGradientBase>(absFcn, upar)); + // std::shared_ptr<AbsPawianMinimizer<FCNGradientBase>> absMinimizerPtr(new MinuitMinimizer<FCNGradientBase>(absFcn, upar)); absMinimizerPtr->minimize(); absMinimizerPtr->printFitResult(evtWeightSumData); diff --git a/AppUtils/Jamfile b/AppUtils/Jamfile index 228343b5..72a51783 100644 --- a/AppUtils/Jamfile +++ b/AppUtils/Jamfile @@ -27,6 +27,6 @@ lib AppUtils : $(TOP)/ggUtils//ggUtils $(TOP)/pipiScatteringUtils//pipiScatteringUtils $(TOP)/qaErrorExtract//qaErrorExtract - : + : <use>$(TOP)//Minuit2 : - : ; + : <library>$(TOP)//Minuit2 ; diff --git a/FitParams/PwaCovMatrix.cc b/FitParams/PwaCovMatrix.cc index d831a83d..76b122ac 100644 --- a/FitParams/PwaCovMatrix.cc +++ b/FitParams/PwaCovMatrix.cc @@ -60,7 +60,7 @@ PwaCovMatrix::PwaCovMatrix(ROOT::Minuit2::MnUserCovariance &theMinuitCovMatrix, _n = theMinuitCovMatrix.Nrow(); unsigned int _nPar = theMinuitParameters.Params().size(); - + int iCov=0; int jCov=0; @@ -97,6 +97,18 @@ PwaCovMatrix::PwaCovMatrix(ROOT::Minuit2::MnUserCovariance &theMinuitCovMatrix, iCov++; jCov=0; } + + //fill data from minuit cov matrix + _dataFromMnCovMatrix=theMinuitCovMatrix.Data(); + + // _dataFromMnCovMatrix.resize(_n*(_n+1)/2); + + //for(unsigned int i=0; i<_nPar; i++){ + // for(unsigned int j=i; j<_nPar; j++){ + // double covValue = theMinuitCovMatrix(iCov, jCov); + // _dataFromMnCovMatrix.push_back(covValue); + // } + //} } diff --git a/FitParams/PwaCovMatrix.hh b/FitParams/PwaCovMatrix.hh index 8866889c..7d99226a 100644 --- a/FitParams/PwaCovMatrix.hh +++ b/FitParams/PwaCovMatrix.hh @@ -27,10 +27,12 @@ #pragma once #include <map> +#include <vector> #include <boost/archive/text_oarchive.hpp> #include <boost/archive/text_iarchive.hpp> #include <boost/serialization/map.hpp> +#include <boost/serialization/vector.hpp> #include "Minuit2/MnUserCovariance.h" #include "Minuit2/MnUserParameters.h" @@ -42,6 +44,7 @@ class PwaCovMatrix friend class boost::serialization::access; unsigned short _n; std::map<std::string, std::map<std::string, double> > _covMatrix; + std::vector<double> _dataFromMnCovMatrix; public: PwaCovMatrix(); @@ -55,12 +58,20 @@ class PwaCovMatrix static bool DiagonalIsValid(const ROOT::Minuit2::MnUserCovariance &theMinuitCovMatrix); bool DiagonalIsValid(); bool CheckCorrelationCoefficients(); - void printElements(); + const std::vector<double>& dataFromMnCovMatrix() const {return _dataFromMnCovMatrix;} + const unsigned int nRow() const {return _n;} + + void printElements(); template<class Archive> void serialize(Archive & ar, const unsigned int version){ ar & _n; ar & _covMatrix; + if(version > 0){ + ar & _dataFromMnCovMatrix; + } } }; + +BOOST_CLASS_VERSION(PwaCovMatrix, 1) diff --git a/MinFunctions/MinuitMinimizer.cc b/MinFunctions/MinuitMinimizer.cc index 058e1419..63972096 100644 --- a/MinFunctions/MinuitMinimizer.cc +++ b/MinFunctions/MinuitMinimizer.cc @@ -46,6 +46,16 @@ template<typename T> MinuitMinimizer<T>::MinuitMinimizer(std::shared_ptr<AbsFcn<T>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar) : AbsPawianMinimizer<T>(theAbsFcnPtr, upar) ,_startMnUserParametersPtr(this->_startPawianParams->mnUserParametersPtr()) + ,_startWithCovMat(false) +{ +} + +template<typename T> +MinuitMinimizer<T>::MinuitMinimizer(std::shared_ptr<AbsFcn<T>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar, std::shared_ptr<MnUserCovariance> mnCovMatrix) : + AbsPawianMinimizer<T>(theAbsFcnPtr, upar) + ,_startMnUserParametersPtr(this->_startPawianParams->mnUserParametersPtr()) + ,_startMnCovMatrix(mnCovMatrix) + ,_startWithCovMat(true) { } @@ -56,12 +66,15 @@ void MinuitMinimizer<T>::minimize(){ // MnMigrad migrad(*_absFcn, _startPawianParams->mnUserParameters()); // MnUserParameters startMnUserP(*_startMnUserParametersPtr); unsigned int stratLevel=GlobalEnv::instance()->parser()->minuitStrategyLevel(); - MnMigrad migrad(*this->_absFcn, *_startMnUserParametersPtr); + std::shared_ptr<MnMigrad> migrad; + if (_startWithCovMat) migrad=std::shared_ptr<MnMigrad>(new MnMigrad(*this->_absFcn, *_startMnUserParametersPtr, *_startMnCovMatrix)); + else migrad=std::shared_ptr<MnMigrad>(new MnMigrad(*this->_absFcn, *_startMnUserParametersPtr)); + // (new MnMigrad(*this->_absFcn, *_startMnUserParametersPtr)); FunctionMinimum* currentFunctionMinimum=0; if(stratLevel==1){ InfoMsg << "start migrad with strategy level " << 1 << endmsg; - currentFunctionMinimum= new FunctionMinimum(migrad(0, GlobalEnv::instance()->parser()->tolerance())); + currentFunctionMinimum= new FunctionMinimum((*migrad)(0, GlobalEnv::instance()->parser()->tolerance())); } else if(stratLevel==2){ MnMigrad migrad2a(*this->_absFcn, *_startMnUserParametersPtr, MnStrategy(2)); @@ -226,3 +239,5 @@ void MinuitMinimizer<T>::dumpFitResult(){ template MinuitMinimizer<FCNBase>::MinuitMinimizer(std::shared_ptr<AbsFcn<FCNBase>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar); template MinuitMinimizer<FCNGradientBase>::MinuitMinimizer(std::shared_ptr<AbsFcn<FCNGradientBase>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar); +template MinuitMinimizer<FCNBase>::MinuitMinimizer(std::shared_ptr<AbsFcn<FCNBase>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar, std::shared_ptr<MnUserCovariance> mnCovMatrix); +template MinuitMinimizer<FCNGradientBase>::MinuitMinimizer(std::shared_ptr<AbsFcn<FCNGradientBase>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar, std::shared_ptr<MnUserCovariance> mnCovMatrix); diff --git a/MinFunctions/MinuitMinimizer.hh b/MinFunctions/MinuitMinimizer.hh index 87b7d242..840226a4 100644 --- a/MinFunctions/MinuitMinimizer.hh +++ b/MinFunctions/MinuitMinimizer.hh @@ -31,6 +31,7 @@ #include "MinFunctions/AbsPawianMinimizer.hh" #include "Minuit2/FunctionMinimum.h" #include "Minuit2/MnUserParameters.h" +#include "Minuit2/MnUserCovariance.h" #include <boost/random/normal_distribution.hpp> @@ -41,7 +42,8 @@ class MinuitMinimizer : public AbsPawianMinimizer<T> { public: MinuitMinimizer(std::shared_ptr<AbsFcn<T>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar); - + MinuitMinimizer(std::shared_ptr<AbsFcn<T>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar, std::shared_ptr<MnUserCovariance> mnCovMatrix); + virtual std::string type() {return "MinuitMinimizer";}; virtual void minimize(); // virtual void printFitResultQA(double evtWeightSumData); @@ -52,6 +54,9 @@ protected: std::shared_ptr<FunctionMinimum> _mnFunctionMinimumFinalPtr; std::shared_ptr<MnUserParameters> _startMnUserParametersPtr; + std::shared_ptr<MnUserCovariance> _startMnCovMatrix; + bool _startWithCovMat; + private: }; diff --git a/MinFunctions/PwaFcnServer.cc b/MinFunctions/PwaFcnServer.cc index 35805978..8c02c25d 100644 --- a/MinFunctions/PwaFcnServer.cc +++ b/MinFunctions/PwaFcnServer.cc @@ -85,18 +85,18 @@ double PwaFcnServer<T>::collectLH() const{ } // Add LLHs of different channels if(lhPrint) output << "current LH = "; - for(auto it = theLHDataMap.begin(); it!=theLHDataMap.end();++it){ - (*it).second.weightSum = _networkServerPtr->weightSum((*it).first); - (*it).second.squaredWeightSum = _networkServerPtr->squaredWeightSum((*it).first); - (*it).second.num_mc = _networkServerPtr->numMCs((*it).first); - double channelLH = AbsLh::mergeLogLhData((*it).second, (*it).first); - result += channelLH; - if(lhPrint) output << std::setprecision(16) << channelLH << "\t"; - if(lhDump) outputLHDump << std::setprecision(16) << channelLH << "\t"; - } - if(lhPrint && theLHDataMap.size() > 1){ - output << "sum = " << result; - } + for(auto it = theLHDataMap.begin(); it!=theLHDataMap.end();++it){ + (*it).second.weightSum = _networkServerPtr->weightSum((*it).first); + (*it).second.squaredWeightSum = _networkServerPtr->squaredWeightSum((*it).first); + (*it).second.num_mc = _networkServerPtr->numMCs((*it).first); + double channelLH = AbsLh::mergeLogLhData((*it).second, (*it).first); + result += channelLH; + if(lhPrint) output << std::setprecision(16) << channelLH << "\t"; + if(lhDump) outputLHDump << std::setprecision(16) << channelLH << "\t"; + } + if(lhPrint && theLHDataMap.size() > 1){ + output << "sum = " << result; + } } if(lhPrint){ @@ -130,7 +130,7 @@ std::vector<double> PwaFcnServer<T>::Gradient(const std::vector<double>& par) co } else if (std::abs(currentVal)<1.e-10) epsilon=_numStepSize*1.e-10; double dx=(currentVal+epsilon)-currentVal; - if(this->_currentPawianParms->HasLimits(i) && std::abs(currentVal-this->_currentPawianParms->UpperLimit(i))<1.e-6){ + if(this->_currentPawianParms->HasLimits(i) && std::abs(currentVal-this->_currentPawianParms->UpperLimit(i))<epsilon){ this->_currentPawianParms->SetValue(i, currentVal-epsilon); ParamDepHandler::instance()->ApplyDependencies(this->_currentPawianParms); double currentLH=collectLH(); diff --git a/PwaUtils/GlobalEnv.cc b/PwaUtils/GlobalEnv.cc index 3b66f194..cb64fb30 100644 --- a/PwaUtils/GlobalEnv.cc +++ b/PwaUtils/GlobalEnv.cc @@ -27,6 +27,7 @@ #include <boost/random.hpp> #include <chrono> #include "PwaUtils/GlobalEnv.hh" +#include "FitParams/PwaCovMatrix.hh" #include "ConfigParser/ParserBase.hh" #include "PwaUtils/AbsLh.hh" #include "Particle/PdtParser.hh" @@ -47,7 +48,8 @@ GlobalEnv* GlobalEnv::instance(){ GlobalEnv::GlobalEnv() : _alreadySetUp(false) , _channelEnvsAlredySetup(false), - _theParser(0) + _theParser(0), + _useCovMatrix(false) //_topDirPath(getenv("TOP_DIR")) //_KMatPath(getenv("KMAT_DIR")) { @@ -355,6 +357,19 @@ void GlobalEnv::setup(ParserBase* theParser){ _outputFileNameSuffix = theParser->outputFileNameSuffix(); _serializationFileName = theParser->serializationFile(); + std::ifstream serializationStream(_serializationFileName.c_str()); + if(!serializationStream.is_open()){ + InfoMsg << "Could not open serialization file." << endmsg; + InfoMsg << "fit/calculation is performed without errors!!!" << endmsg; + _useCovMatrix=false; + } + else{ + _pwaCovMatrix = std::shared_ptr<PwaCovMatrix>(new PwaCovMatrix); + boost::archive::text_iarchive boostInputArchive(serializationStream); + boostInputArchive >> *_pwaCovMatrix; + _useCovMatrix=true; + } + // pdtTable PdtParser pdtParser; std::string theSourcePath=getenv("TOP_DIR"); diff --git a/PwaUtils/GlobalEnv.hh b/PwaUtils/GlobalEnv.hh index df5f5855..f1179a15 100644 --- a/PwaUtils/GlobalEnv.hh +++ b/PwaUtils/GlobalEnv.hh @@ -35,6 +35,7 @@ class ParserBase; class ParticleTable; class AbsPawianParameters; +class PwaCovMatrix; typedef std::vector<std::pair<std::shared_ptr<AbsChannelEnv>, short> > ChannelEnvList; @@ -83,7 +84,8 @@ public: std::string topDirPath() const {return _topDirPath;} std::string KMatrixStorePath() const {return _KMatStorePath;} std::string evtStorePath() const {return _evtStorePath;} - + std::shared_ptr<PwaCovMatrix> pwaCovMatrix() const {return _pwaCovMatrix;} + bool useCovMatrix() const {return _useCovMatrix;} private: GlobalEnv(); static GlobalEnv* _instance; @@ -103,4 +105,6 @@ private: std::map<std::string, std::string> _toBeReplacedSuffixMap; std::map<std::string, std::string> _alreadyReplacedSuffixMap; std::map<std::string, std::string> _fitParamReplacementMap; + std::shared_ptr<PwaCovMatrix> _pwaCovMatrix; + bool _useCovMatrix; }; diff --git a/benchmark/pbarpTopi0pi0eta/pbarpReactionDefault.cfg b/benchmark/pbarpTopi0pi0eta/pbarpReactionDefault.cfg index 486cdce1..ae1deb1f 100644 --- a/benchmark/pbarpTopi0pi0eta/pbarpReactionDefault.cfg +++ b/benchmark/pbarpTopi0pi0eta/pbarpReactionDefault.cfg @@ -1,6 +1,8 @@ errLogMode = trace #errLogMode = notice +#minimumTolerance = 100 + stepSizeLhPrint = 10 stepSizeLhDump = 20 stepSizeTimer = 20 @@ -17,6 +19,7 @@ orderInFile = E Px Py Pz noOfDataEvents = 15000 ratioMcToData = 3 +#ratioMcToData = 1 useDataEventWeight = false useMCEventWeight = false -- GitLab