diff --git a/AppUtils/AppBase.cc b/AppUtils/AppBase.cc
index b2ffd8e163c13a419b038e72c2a5e33b8ab1a120..81bd3dee18681f611d52a1591b155e98049046fa 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 228343b5db4f0cd6a31822c993e5f4498aa38e66..72a517835afb67dda29a925c5ae46f6147881058 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 d831a83d01f199a0db81737293d18187a01537b8..76b122ac7b7f44953c249deb247827ed976faed7 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 8866889c464d3e152fca0f04fee5caf1b9d2d6c7..7d99226ae17d8dc4987c85d2978cb2ef83ba6f57 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 058e141910afa84954220ed4d08d578d34160625..639720964c2b53975d83172dffd83eda450212a2 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 87b7d24214d4dabb05bcaf8e55e7c3aa29e2f35b..840226a4f0764eb0a800b596ae65c91fbad8e1e1 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 35805978bd0465f4537ee0c633aa6b8113931193..8c02c25d488592325ec1338953546be8233d1475 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 3b66f194f5c5ad0a75d6a3299340541d1de26a35..cb64fb30d69686d92faac8f24d7948eee9199ee5 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 df5f5855dc08e7f5ce6b134c9370ab8cfc36461f..f1179a150d48472ad25366da21cfdeb34914ae32 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 486cdce1fea76f063271b9f9f55a089d82b244b0..ae1deb1f6891add0523ca7390bf1e5665073cfdc 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