From 214457e78ceb101908ebc102356d9d3b79469dad Mon Sep 17 00:00:00 2001
From: Julian Pychy <julian@pc17.ep1.rub.de>
Date: Mon, 18 Nov 2013 18:30:13 +0100
Subject: [PATCH] workaround for minuit returning bad covariance matrix,
improved startup performance
---
Examples/epem/epemReactionApp.cc | 6 ++---
Examples/pbarp/pbarpReactionApp.cc | 5 ++--
PwaUtils/AppBase.cc | 42 ++++++++++++++++++++++++++----
PwaUtils/PwaCovMatrix.cc | 16 ++++++++++++
PwaUtils/PwaCovMatrix.hh | 2 ++
5 files changed, 59 insertions(+), 12 deletions(-)
diff --git a/Examples/epem/epemReactionApp.cc b/Examples/epem/epemReactionApp.cc
index ad2b9eae..bd46e63a 100644
--- a/Examples/epem/epemReactionApp.cc
+++ b/Examples/epem/epemReactionApp.cc
@@ -82,9 +82,6 @@ int main(int __argc,char *__argv[]){
Info << welcomeScreen << endmsg;
Info << "Compiled " << __DATE__ << " " << __TIME__ << endmsg;
- // Disable output buffering
- setvbuf(stdout, NULL, _IONBF, 0);
-
// Parse the command line
epemParser* theAppParams=new epemParser(__argc, __argv);
@@ -162,7 +159,8 @@ int main(int __argc,char *__argv[]){
theAppBase.fixParams(upar,fixedParams);
const unsigned int noOfFreeFitParams = upar.VariableParameters();
-
+ // Disable output buffering
+ setvbuf(stdout, NULL, _IONBF, 0);
if(mode == "client"){
diff --git a/Examples/pbarp/pbarpReactionApp.cc b/Examples/pbarp/pbarpReactionApp.cc
index 54f06285..ca6dda92 100644
--- a/Examples/pbarp/pbarpReactionApp.cc
+++ b/Examples/pbarp/pbarpReactionApp.cc
@@ -90,9 +90,6 @@ int main(int __argc,char *__argv[]){
Info << welcomeScreen << endmsg;
Info << "Compiled " << __DATE__ << " " << __TIME__ << endmsg;
- // Disable output buffering
- setvbuf(stdout, NULL, _IONBF, 0);
-
// Parse the command line
pbarpParser* theAppParams=new pbarpParser(__argc, __argv);
@@ -177,6 +174,8 @@ int main(int __argc,char *__argv[]){
theAppBase.fixParams(upar,fixedParams);
const unsigned int noOfFreeFitParams = upar.VariableParameters();
+ // Disable output buffering
+ setvbuf(stdout, NULL, _IONBF, 0);
if(mode == "client"){
diff --git a/PwaUtils/AppBase.cc b/PwaUtils/AppBase.cc
index 6b53aa69..699fd5f8 100644
--- a/PwaUtils/AppBase.cc
+++ b/PwaUtils/AppBase.cc
@@ -294,11 +294,43 @@ FunctionMinimum AppBase::migradDefault(AbsFcn& theFcn, MnUserParameters& upar){
MnMigrad migrad(theFcn, upar);
Info <<"start migrad "<< endmsg;
FunctionMinimum funcMin = migrad();
- if(!funcMin.IsValid()) {
- //try with higher strategy
- Info <<"FM is invalid, try with strategy = 2."<< endmsg;
- MnMigrad migrad2(theFcn, funcMin.UserState(), MnStrategy(2));
- funcMin = migrad2();
+
+ if(funcMin.IsValid()){
+ return funcMin;
+ }
+
+ // Two more tries to get a valid result unsing strategy 2
+ for(int i=0; i<2; i++){
+ Warning <<"FM is invalid, try with strategy = 2."<< endmsg;
+
+ // Check minimum covariance matrix
+ bool badCovarianceDiagonal=false;
+ if(funcMin.HasCovariance()){
+ badCovarianceDiagonal = !PwaCovMatrix::DiagonalIsValid(funcMin.UserCovariance());
+ }
+
+ if(badCovarianceDiagonal){
+ Warning << "Using default errors" << endmsg;
+ MnUserParameters newParams = upar;
+ for(unsigned int i=0; i< funcMin.UserParameters().Params().size();i++){
+ newParams.SetValue(i, funcMin.UserParameters().Params().at(i));
+ }
+ MnMigrad migrad2(theFcn, newParams, MnStrategy(2));
+ funcMin = migrad2();
+ }
+ else{
+ MnUserParameters newParams = upar;
+ for(unsigned int i=0; i< funcMin.UserParameters().Params().size();i++){
+ newParams.SetValue(i, funcMin.UserParameters().Params().at(i));
+ newParams.SetError(i, funcMin.UserParameters().Errors().at(i));
+ }
+ MnMigrad migrad2(theFcn, newParams, MnStrategy(2));
+ funcMin = migrad2();
+ }
+
+ if(funcMin.IsValid()){
+ break;
+ }
}
return funcMin;
diff --git a/PwaUtils/PwaCovMatrix.cc b/PwaUtils/PwaCovMatrix.cc
index e100f998..ae5574e6 100644
--- a/PwaUtils/PwaCovMatrix.cc
+++ b/PwaUtils/PwaCovMatrix.cc
@@ -43,6 +43,8 @@ PwaCovMatrix::PwaCovMatrix(ROOT::Minuit2::MnUserCovariance &theMinuitCovMatrix,
ROOT::Minuit2::MnUserParameters &theMinuitParameters,
fitParams &theFitParams)
{
+ DiagonalIsValid(theMinuitCovMatrix);
+
_n = theMinuitCovMatrix.Nrow();
unsigned int _nPar = theMinuitParameters.Params().size();
@@ -106,3 +108,17 @@ double PwaCovMatrix::GetElement(std::string parameter1, std::string parameter2){
return _covMatrix[parameter1][parameter2];
}
+
+
+bool PwaCovMatrix::DiagonalIsValid(const ROOT::Minuit2::MnUserCovariance &theMinuitCovMatrix){
+
+ bool result = true;
+ for(unsigned int i=0; i<theMinuitCovMatrix.Nrow(); i++){
+ if(theMinuitCovMatrix(i, i) < 0){
+ Warning << "Covariance element (" << i << "," << i << ")"
+ << " = " << theMinuitCovMatrix(i,i) << " < 0" << endmsg;
+ result = false;
+ }
+ }
+ return result;
+}
diff --git a/PwaUtils/PwaCovMatrix.hh b/PwaUtils/PwaCovMatrix.hh
index 8f344acc..f0f5ae18 100644
--- a/PwaUtils/PwaCovMatrix.hh
+++ b/PwaUtils/PwaCovMatrix.hh
@@ -52,6 +52,8 @@ class PwaCovMatrix
ROOT::Minuit2::MnUserParameters &theMinuitParameters,
fitParams &theFitParams);
double GetElement(std::string parameter1, std::string parameter2);
+ static bool DiagonalIsValid(const ROOT::Minuit2::MnUserCovariance &theMinuitCovMatrix);
+
template<class Archive>
void serialize(Archive & ar, const unsigned int version){
--
GitLab