From 33f9ab616cd5c3983faebbe0464757be2f3daaa1 Mon Sep 17 00:00:00 2001
From: Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
Date: Tue, 29 Sep 2009 09:01:04 +0000
Subject: [PATCH] Merging r30520 through r30523 for
https://root.cern.ch/svn/root/branches/dev/roostats
- fix an unitialized variable in BayesianCalculator
- remove log messages in LikelihoodInterval
git-svn-id: http://root.cern.ch/svn/root/trunk@30524 27541ba8-7e3a-0410-8455-c3a389f83636
---
roofit/roostats/inc/BayesianCalculator.h | 2 -
roofit/roostats/src/BayesianCalculator.cxx | 19 ++---
roofit/roostats/src/LikelihoodInterval.cxx | 71 +++++--------------
tutorials/roostats/rs701_BayesianCalculator.C | 1 +
4 files changed, 30 insertions(+), 63 deletions(-)
diff --git a/roofit/roostats/inc/BayesianCalculator.h b/roofit/roostats/inc/BayesianCalculator.h
index e58d4f68156..2cb0e6238b0 100644
--- a/roofit/roostats/inc/BayesianCalculator.h
+++ b/roofit/roostats/inc/BayesianCalculator.h
@@ -91,8 +91,6 @@ namespace RooStats {
double fSize;
- mutable Double_t fLowerLimit;
- mutable Double_t fUpperLimit;
protected:
diff --git a/roofit/roostats/src/BayesianCalculator.cxx b/roofit/roostats/src/BayesianCalculator.cxx
index 15a84732749..72763ecc26f 100644
--- a/roofit/roostats/src/BayesianCalculator.cxx
+++ b/roofit/roostats/src/BayesianCalculator.cxx
@@ -37,7 +37,8 @@ BayesianCalculator::BayesianCalculator() :
fData(0),
fPdf(0),
fPriorPOI(0),
- fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0)
+ fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
+ fSize(0.05)
{
// default constructor
}
@@ -54,11 +55,10 @@ BayesianCalculator::BayesianCalculator( /* const char* name, const char* title,
fPOI(POI),
fPriorPOI(&priorPOI),
fNuisanceParameters(*nuisanceParameters),
- fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0)
+ fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
+ fSize(0.05)
{
// constructor
- fLowerLimit = -999;
- fUpperLimit = +999;
}
BayesianCalculator::BayesianCalculator( RooAbsData& data,
@@ -66,7 +66,8 @@ BayesianCalculator::BayesianCalculator( RooAbsData& data,
fData(&data),
fPdf(model.GetPdf()),
fPriorPOI( model.GetPriorPdf()),
- fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0)
+ fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
+ fSize(0.05)
{
// constructor from Model Config
SetModel(model);
@@ -171,17 +172,19 @@ SimpleInterval* BayesianCalculator::GetInterval() const
assert(poi);
double y = fSize;
- brf.findRoot(fLowerLimit,poi->getMin(),poi->getMax(),y);
+ double lowerLimit = 0;
+ double upperLimit = 0;
+ brf.findRoot(lowerLimit,poi->getMin(),poi->getMax(),y);
y=1-fSize;
- brf.findRoot(fUpperLimit,poi->getMin(),poi->getMax(),y);
+ brf.findRoot(upperLimit,poi->getMin(),poi->getMax(),y);
delete cdf_bind;
delete cdf;
delete nll;
TString interval_name = TString("BayesianInterval_a") + TString(this->GetName());
- SimpleInterval* interval = new SimpleInterval(interval_name,"SimpleInterval from BayesianCalculator",poi,fLowerLimit,fUpperLimit);
+ SimpleInterval* interval = new SimpleInterval(interval_name,"SimpleInterval from BayesianCalculator",poi,lowerLimit,upperLimit);
return interval;
}
diff --git a/roofit/roostats/src/LikelihoodInterval.cxx b/roofit/roostats/src/LikelihoodInterval.cxx
index 027d4c64d5a..42cab8c0345 100644
--- a/roofit/roostats/src/LikelihoodInterval.cxx
+++ b/roofit/roostats/src/LikelihoodInterval.cxx
@@ -263,10 +263,9 @@ Double_t LikelihoodInterval::LowerLimit(RooRealVar& param)
step=thisArgVal+step-x_app;
if (fabs(step)<1E-5) {
- std::cout<<"WARNING lower limit is outside the parameters bounderies. Abort!"<<std::endl;
+ ccoutW(Eval) << "Lower limit is outside the parameters bounderies. Abort!" <<std::endl;
delete newProfile;
double ret = myarg->getMin();
- //delete myarg;
return ret;
}
}
@@ -277,7 +276,7 @@ Double_t LikelihoodInterval::LowerLimit(RooRealVar& param)
// If L is below target and we are at boundary stop here
if ((fabs(myarg->getVal()-myarg->getMin())<1e-5) && (L<target)) {
- std::cout<<"WARNING lower limit is outside the parameters bounderies (L at lower bound of " << myarg->GetName() << " is " << L
+ ccoutW(Eval) <<"WARNING lower limit is outside the parameters bounderies (L at lower bound of " << myarg->GetName() << " is " << L
<< ", which is less than target value of " << target << "). Abort!"<<std::endl;
delete newProfile;
double ret = myarg->getMin();
@@ -287,7 +286,7 @@ Double_t LikelihoodInterval::LowerLimit(RooRealVar& param)
//Compute the linear function
a=(L-L_old)/(step);
if (fabs(a)<1E-3) {
- std::cout<<"WARNING: the slope of the Likelihood linear approximation is close to zero" <<std::endl;
+ ccoutD(Eval) <<"WARNING: the slope of the Likelihood linear approximation is close to zero" <<std::endl;
}
b=L-a*thisArgVal;
//approximate the position of the desired L value
@@ -298,19 +297,17 @@ Double_t LikelihoodInterval::LowerLimit(RooRealVar& param)
diff=L-target;
if(a>0) {
- std::cout<<"WARNING: you are on the right of the MLE. Step towards MLE"<<std::endl;
+ ccoutD(Eval) <<"WARNING: you are on the right of the MLE. Step towards MLE"<<std::endl;
step=-(myarg->getMax()-myarg->getMin())/100; //Do a constant step, typical for the dimenions of the problem towards the minimum
}
}
-
-
- std::cout<<"LL search Iterations:"<< nIterations<<std::endl;
+
+ ccoutD(Eval) <<"LL search Iterations:"<< nIterations<<std::endl;
delete newProfile;
double ret=myarg->getVal();
- //delete myarg;
return ret;
}
@@ -332,7 +329,6 @@ Double_t LikelihoodInterval::UpperLimit(RooRealVar& param)
Double_t maxStep = (myarg->getMax()-myarg->getMin())/10 ;
- bool toggleplot=false; //For debugging purposes, could think of interfacing a Filepattern to outside to have controlplots stored
double step=1;
double L= newProfile->getVal();
double L_old=0;
@@ -346,41 +342,29 @@ Double_t LikelihoodInterval::UpperLimit(RooRealVar& param)
L_old=L;
thisArgVal=thisArgVal+step;
if (thisArgVal>myarg->getMax())
- {
- std::cout<<"WARNING: near the upper boundery"<<std::endl;
+ {
+ ccoutD(Eval) <<"WARNING: near the upper boundery"<<std::endl;
thisArgVal=myarg->getMax();
step=thisArgVal+step-x_app;
//std::cout<<"DEBUG: step:"<<step<<" thistArgVal:"<<thisArgVal<<std::endl;
if (fabs(step)<1E-5) {
- toggleplot=true;
- std::cout<<"WARNING upper limit is outside the parameters bounderies. Abort!"<<std::endl;
- // TCanvas c1;
- // c1.cd();
- // RooPlot* frame=myarg->frame();
- // newProfile->plotOn(frame);
- // frame->Draw();
- // TString fname("");
- // fname+=MLE;
- // fname+="_"; fname+="abort.ps";
- // c1.Print(fname);
- // delete frame;
+ ccoutW(Eval) <<"Upper limit is outside the parameters bounderies. Abort!"<<std::endl;
delete newProfile;
double ret=myarg->getMax();
- //delete myarg;
return ret;
}
- }
+ }
myarg->setVal( thisArgVal );
L=newProfile->getVal();
// If L is below target and we are at boundary stop here
if ((fabs(myarg->getVal()-myarg->getMax())<1e-5) && (L<target)) {
- std::cout<<"WARNING upper limit is outside the parameters bounderies (L at upper bound of " << myarg->GetName() << " is " << L
+ ccoutW(Eval) <<"WARNING upper limit is outside the parameters bounderies (L at upper bound of " << myarg->GetName() << " is " << L
<< ", which is less than target value of " << target << "). Abort!"<<std::endl;
- delete newProfile;
- double ret = myarg->getMax();
- return ret;
+ delete newProfile;
+ double ret = myarg->getMax();
+ return ret;
}
@@ -388,8 +372,7 @@ Double_t LikelihoodInterval::UpperLimit(RooRealVar& param)
a=(L-L_old)/(step);
if (fabs(a)<1E-3){
- std::cout<<"WARNING: the slope of the Likelihood linear approximation is close to zero."<<std::endl;
- toggleplot=true;
+ ccoutD(Eval) <<"WARNING: the slope of the Likelihood linear approximation is close to zero."<<std::endl;
}
b=L-a*thisArgVal;
//approximate the position of the desired L value
@@ -401,38 +384,20 @@ Double_t LikelihoodInterval::UpperLimit(RooRealVar& param)
//If slope is negative you are below the minimum
if(a<0) {
- std::cout<<"WARNING: you are on the left of the MLE. Step towards MLE"<<std::endl;
+ ccoutD(Eval)<<"WARNING: you are on the left of the MLE. Step towards MLE"<<std::endl;
step=(myarg->getMax()-myarg->getMin())/100; //Do a constant step, typical for the dimenions of the problem towards the minimum
//L_old=0;
}
}
-
- //________________________________________
- //Controlplots for debugging
- if (toggleplot) {
- // TCanvas c1;
- // c1.cd();
- // RooPlot* frame=myarg->frame();
- // newProfile->plotOn(frame);
- // frame->Draw();
- // TString fname("/afs/cern.ch/user/r/ruthmann/development/c++/2009.08.14_higgs_combination_7tev/limit_controlplots/");
- // fname+=MLE;
- // fname+="_"; fname+=myarg->getVal(); fname+=".ps";
- // c1.Print(fname);
- // delete frame;
- }
- //_______________________________________
-
- std::cout<<"UL search Iterations:"<< nIterations<<std::endl;
+
+ ccoutD(Eval) <<"UL search Iterations:"<< nIterations<<std::endl;
delete newProfile;
double ret=myarg->getVal();
- //delete myarg;
return ret;
- //return x_app;
RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;
}
diff --git a/tutorials/roostats/rs701_BayesianCalculator.C b/tutorials/roostats/rs701_BayesianCalculator.C
index e6f18ae6ad0..0ec9f62c095 100644
--- a/tutorials/roostats/rs701_BayesianCalculator.C
+++ b/tutorials/roostats/rs701_BayesianCalculator.C
@@ -40,6 +40,7 @@ void rs701_BayesianCalculator()
data.add(RooArgSet(*(w->var("x"))),w->var("n")->getVal());
BayesianCalculator bcalc(data,*model,RooArgSet(*POI),*priorPOI,&nuisanceParameters);
+ bcalc.SetTestSize(0.05);
SimpleInterval* interval = bcalc.GetInterval();
std::cout << "90% CL interval: [ " << interval->LowerLimit() << " - " << interval->UpperLimit() << " ] or 95% CL limits\n";
bcalc.PlotPosterior();
--
GitLab