From 8cdfdd33564094d9200e631bf77989532f94cc67 Mon Sep 17 00:00:00 2001
From: Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
Date: Mon, 7 May 2012 16:41:44 +0000
Subject: [PATCH] merge changes from roostats development branch up to
revision 44155
Modifications by Jakub Ademek and Max Baak :
ensure that ToyMCSampler can return TTree with detailed fit output.
See for a detailed description:
https://indico.cern.ch/getFile.py/access?contribId=0&resId=3&materialId=slides&confId=182211
Modifications by Sven Kreiss:
ProfileLikelihoodTestStat:
- added option to add errors and pulls to detailed output
Revision 44128 - Directory Listing
Modified Fri May 4 21:33:08 2012 UTC (2 days, 17 hours ago) by sven
ToyMCSampler, ToyMCStudy, SamplingDistribution, DetailedOutputAggregator:
- bugfixes for the naming of the result
- return toys with weight=1 when generating unweighted toys (not -1)
- make it optional to return errors and pulls with detailed output
- from Tim Adye: catch NaNs at the right place
ProfileLikelihoodTestStat:
- added option to add errors and pulls to detailed output
FrequentistCalculator (from Tim Adye):
- skip the profiling if there are no free nuisance parameters
--This line, and those below, wil
_M roostats
_M roostats/src
M roostats/src/HypoTestResult.cxx
M roostats/src/FrequentistCalculator.cxx
M roostats/src/HypoTestCalculatorGeneric.cxx
M roostats/src/RooStatsUtils.cxx
M roostats/src/ToyMCSampler.cxx
M roostats/src/ToyMCStudy.cxx
M roostats/src/ProfileLikelihoodTestStat.cxx
M roostats/src/SamplingDistribution.cxx
_M roostats/inc
M roostats/inc/LinkDef.h
M roostats/inc/RooStats/ProfileLikelihoodTestStat.h
_M roostats/inc/RooStats/ToyMCSamplerOld.h
M roostats/inc/RooStats/FrequentistCalculator.h
M roostats/inc/RooStats/HypoTestCalculatorGeneric.h
A + roostats/inc/RooStats/DetailedOutputAggregator.h
M roostats/inc/RooStats/RatioOfProfiledLikelihoodsTestStat.h
M roostats/inc/RooStats/ToyMCStudy.h
M roostats/inc/RooStats/TestStatistic.h
M roostats/inc/RooStats/SamplingDistribution.h
M roostats/inc/RooStats/HypoTestResult.h
A + roostats/inc/RooStats/MinNLLTestStat.h
M roostats/inc/RooStats/RooStatsUtils.h
M roostats/inc/RooStats/ToyMCSampler.h
git-svn-id: http://root.cern.ch/svn/root/trunk@44157 27541ba8-7e3a-0410-8455-c3a389f83636
---
roofit/roostats/inc/LinkDef.h | 4 +
.../inc/RooStats/DetailedOutputAggregator.h | 172 ++++++++++++++++++
.../inc/RooStats/FrequentistCalculator.h | 28 ++-
.../inc/RooStats/HypoTestCalculatorGeneric.h | 3 +
roofit/roostats/inc/RooStats/HypoTestResult.h | 3 +
roofit/roostats/inc/RooStats/MinNLLTestStat.h | 107 +++++++++++
.../inc/RooStats/ProfileLikelihoodTestStat.h | 36 ++--
.../RatioOfProfiledLikelihoodsTestStat.h | 57 +++++-
roofit/roostats/inc/RooStats/RooStatsUtils.h | 8 +
.../inc/RooStats/SamplingDistribution.h | 2 +-
roofit/roostats/inc/RooStats/TestStatistic.h | 2 +-
roofit/roostats/inc/RooStats/ToyMCSampler.h | 9 +-
roofit/roostats/inc/RooStats/ToyMCStudy.h | 1 +
roofit/roostats/src/FrequentistCalculator.cxx | 117 ++++++++----
.../src/HypoTestCalculatorGeneric.cxx | 9 +
roofit/roostats/src/HypoTestResult.cxx | 12 +-
.../src/ProfileLikelihoodTestStat.cxx | 60 +++---
roofit/roostats/src/RooStatsUtils.cxx | 80 ++++++++
roofit/roostats/src/SamplingDistribution.cxx | 29 ++-
roofit/roostats/src/ToyMCSampler.cxx | 91 +++------
roofit/roostats/src/ToyMCStudy.cxx | 2 +-
21 files changed, 662 insertions(+), 170 deletions(-)
create mode 100755 roofit/roostats/inc/RooStats/DetailedOutputAggregator.h
create mode 100755 roofit/roostats/inc/RooStats/MinNLLTestStat.h
diff --git a/roofit/roostats/inc/LinkDef.h b/roofit/roostats/inc/LinkDef.h
index b2a2cef3ee4..21425817a40 100644
--- a/roofit/roostats/inc/LinkDef.h
+++ b/roofit/roostats/inc/LinkDef.h
@@ -61,6 +61,8 @@
#pragma link C++ class RooStats::HybridPlot+;
#pragma link C++ class RooStats::HybridResult+;
+#pragma link C++ class RooStats::DetailedOutputAggregator+;
+
#pragma link C++ class RooStats::TestStatSampler+; // interface, not concrete
#pragma link C++ class RooStats::DebuggingSampler+;
#pragma link C++ class RooStats::ToyMCSampler+;
@@ -77,6 +79,7 @@
#pragma link C++ class RooStats::NumEventsTestStat+;
#pragma link C++ class RooStats::SimpleLikelihoodRatioTestStat+;
#pragma link C++ class RooStats::MaxLikelihoodEstimateTestStat+;
+#pragma link C++ class RooStats::MinNLLTestStat+;
#pragma link C++ class RooStats::SamplingDistribution+;
#pragma link C++ class RooStats::NeymanConstruction+;
@@ -121,6 +124,7 @@
#pragma link C++ function RooStats::SetAllConstant(const RooAbsCollection & , bool );
#pragma link C++ function RooStats::MakeNuisancePdf(RooAbsPdf & , const RooArgSet &, const char * );
#pragma link C++ function RooStats::MakeNuisancePdf(const RooStats::ModelConfig & , const char * );
+#pragma link C++ function RooStats::GetAsTTree(TString, TString, const RooDataSet&);
// need for auto_ptr object in Likelihoodinterval since they are forwd declared
diff --git a/roofit/roostats/inc/RooStats/DetailedOutputAggregator.h b/roofit/roostats/inc/RooStats/DetailedOutputAggregator.h
new file mode 100755
index 00000000000..ab06791fe5d
--- /dev/null
+++ b/roofit/roostats/inc/RooStats/DetailedOutputAggregator.h
@@ -0,0 +1,172 @@
+// @(#)root/roostats:$Id: DetailedOutputAggregator.h 37084 2010-11-29 21:37:13Z moneta $
+// Author: Sven Kreiss, Kyle Cranmer Nov 2010
+/*************************************************************************
+ * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
+ * All rights reserved. *
+ * *
+ * For the licensing terms see $ROOTSYS/LICENSE. *
+ * For the list of contributors see $ROOTSYS/README/CREDITS. *
+ *************************************************************************/
+
+#ifndef ROOSTATS_DetailedOutputAggregator
+#define ROOSTATS_DetailedOutputAggregator
+
+//_________________________________________________
+/*
+ BEGIN_HTML
+ <p>
+ This class is designed to aid in the construction of RooDataSets and RooArgSets,
+ particularly those naturally arising in fitting operations.
+
+ Typically, the usage of this class is as follows:
+ <ol>
+ <li> create DetailedOutputAggregator instance </li>
+ <li> use AppendArgSet to add value sets to be stored as one row of the dataset </li>
+ <li> call CommitSet when an entire row's worth of values has been added </li>
+ <li> repeat steps 2 and 3 until all rows have been added </li>
+ <li> call GetAsDataSet to extract result RooDataSet </li>
+ </ol>
+
+ </p>
+ END_HTML
+ */
+//
+
+#include "RooFitResult.h"
+#include "RooPullVar.h"
+#include "RooDataSet.h"
+
+namespace RooStats {
+
+ class DetailedOutputAggregator {
+
+ public:
+ // Translate the given fit result to a RooArgSet in a generic way.
+ // Prefix is prepended to all variable names.
+ static RooArgSet *GetAsArgSet(RooFitResult *result, TString prefix="", bool withErrorsAndPulls=false) {
+ RooArgSet *detailedOutput = new RooArgSet;
+ const RooArgSet &detOut = result->floatParsFinal();
+ const RooArgSet &truthSet = result->floatParsInit();
+ RooRealVar* var(0);
+ TIterator *it = detOut.createIterator();
+ for(;(var = dynamic_cast<RooRealVar*>(it->Next()));) {
+ RooAbsArg* clone = var->cloneTree(TString().Append(prefix).Append(var->GetName()));
+ clone->SetTitle( TString().Append(prefix).Append(var->GetTitle()) );
+ detailedOutput->addOwned(*clone);
+
+ if( withErrorsAndPulls ) {
+ TString pullname = TString().Append(prefix).Append(TString::Format("%s_pull", var->GetName()));
+ // TString pulldesc = TString::Format("%s pull for fit %u", var->GetTitle(), fitNumber);
+ RooRealVar* truth = dynamic_cast<RooRealVar*>(truthSet.find(var->GetName()));
+ RooPullVar pulltemp("temppull", "temppull", *var, *truth);
+ RooRealVar* pull = new RooRealVar(pullname, pullname, pulltemp.getVal());
+ detailedOutput->addOwned(*pull);
+
+ TString errloname = TString().Append(prefix).Append(TString::Format("%s_errlo", var->GetName()));
+ // TString errlodesc = TString::Format("%s low error for fit %u", var->GetTitle(), fitNumber);
+ RooRealVar* errlo = new RooRealVar(errloname, errloname, var->getErrorLo());
+ detailedOutput->addOwned(*errlo);
+
+ TString errhiname = TString().Append(prefix).Append(TString::Format("%s_errhi", var->GetName()));
+ // TString errhidesc = TString::Format("%s high error for fit %u", var->GetTitle(), fitNumber);
+ RooRealVar* errhi = new RooRealVar(errhiname, errhiname, var->getErrorHi());
+ detailedOutput->addOwned(*errhi);
+ }
+ }
+ delete it;
+
+ // monitor a few more variables
+ detailedOutput->addOwned( *new RooRealVar(TString().Append(prefix).Append("minNLL"), TString().Append(prefix).Append("minNLL"), result->minNll() ) );
+ detailedOutput->addOwned( *new RooRealVar(TString().Append(prefix).Append("fitStatus"), TString().Append(prefix).Append("fitStatus"), result->status() ) );
+ detailedOutput->addOwned( *new RooRealVar(TString().Append(prefix).Append("covQual"), TString().Append(prefix).Append("covQual"), result->covQual() ) );
+ detailedOutput->addOwned( *new RooRealVar(TString().Append(prefix).Append("numInvalidNLLEval"), TString().Append(prefix).Append("numInvalidNLLEval"), result->numInvalidNLL() ) );
+ return detailedOutput;
+ }
+
+ DetailedOutputAggregator() {
+ result = NULL;
+ builtSet = NULL;
+ }
+
+ // For each variable in aset, prepend prefix to its name and add
+ // to the internal store. Note this will not appear in the produced
+ // dataset unless CommitSet is called.
+ void AppendArgSet(const RooArgSet *aset, TString prefix="") {
+ if (aset == NULL) {
+ // silently ignore
+ //std::cout << "Attempted to append NULL" << endl;
+ return;
+ }
+ if (builtSet == NULL) {
+ builtSet = new RooArgSet();
+ }
+ TIterator* iter = aset->createIterator();
+ while(RooRealVar* v = dynamic_cast<RooRealVar*>( iter->Next() ) ) {
+ TString renamed(TString::Format("%s%s", prefix.Data(), v->GetName()));
+ if (result != NULL) {
+ // we already commited an argset once, so we expect all columns to already be in the set
+ builtSet->setRealValue(renamed, v->getVal());
+ }
+ else {
+ // we never commited, so by default all columns are expected to not exist
+ RooRealVar *var = new RooRealVar(renamed, v->GetTitle(), v->getVal());
+ if (!builtSet->addOwned(*var)) {
+ delete var;
+ builtSet->setRealValue(renamed, v->getVal());
+ }
+ }
+ }
+ delete iter;
+ }
+
+ // Commit to the result RooDataSet.
+ void CommitSet(double weight=1.0) {
+ if (result == NULL) {
+ result = new RooDataSet("", "",
+ RooArgSet( *(new RooRealVar("weight","weight",1.0)), "tmpSet" ), "weight");
+ InitializeColumns(result, builtSet);
+ }
+ result->add(*builtSet, weight);
+ TIterator* iter = builtSet->createIterator();
+ while(RooAbsArg* v = dynamic_cast<RooAbsArg*>( iter->Next() ) ) {
+ builtSet->setRealValue(v->GetName(), -999.0);
+ }
+ delete iter;
+ }
+
+ RooDataSet *GetAsDataSet(TString name, TString title) {
+ RooDataSet* temp = NULL;
+ if( result ) {
+ temp = new RooDataSet( *result );
+ temp->SetNameTitle( name.Data(), title.Data() );
+ }else{
+ temp = new RooDataSet(name.Data(), title.Data(),
+ RooArgSet( *(new RooRealVar("weight","weight",1.0)), "tmpSet" ), "weight");
+ }
+
+ return temp;
+ }
+
+ virtual ~DetailedOutputAggregator() {
+ if (result != NULL) delete result;
+ if (builtSet != NULL) delete builtSet;
+ }
+
+ private:
+ RooDataSet *result;
+ RooArgSet *builtSet;
+
+ void InitializeColumns(RooDataSet *dset, RooArgSet *aset) {
+ TIterator* iter = aset->createIterator();
+ while(RooAbsArg* v = dynamic_cast<RooAbsArg*>( iter->Next() ) ) {
+ dset->addColumn( *(new RooRealVar(v->GetName(), v->GetTitle(), -1.0)));
+ }
+ delete iter;
+ }
+
+ protected:
+ ClassDef(DetailedOutputAggregator,1)
+ };
+}
+
+#endif
diff --git a/roofit/roostats/inc/RooStats/FrequentistCalculator.h b/roofit/roostats/inc/RooStats/FrequentistCalculator.h
index d5594577d21..eabdfcad58b 100644
--- a/roofit/roostats/inc/RooStats/FrequentistCalculator.h
+++ b/roofit/roostats/inc/RooStats/FrequentistCalculator.h
@@ -31,6 +31,12 @@ END_HTML
#include "RooStats/ToyMCSampler.h"
#endif
+#ifndef ROOSTATS_DetailedOutputAggregator
+#include "RooStats/DetailedOutputAggregator.h"
+#endif
+
+#include "RooFitResult.h"
+
namespace RooStats {
class FrequentistCalculator : public HypoTestCalculatorGeneric {
@@ -48,13 +54,16 @@ namespace RooStats {
fNToysNull(-1),
fNToysAlt(-1),
fNToysNullTail(0),
- fNToysAltTail(0)
+ fNToysAltTail(0),
+ fFitInfo(NULL),
+ fStoreFitInfo(false)
{
}
~FrequentistCalculator() {
if( fConditionalMLEsNull ) delete fConditionalMLEsNull;
- if( fConditionalMLEsAlt ) delete fConditionalMLEsAlt;
+ if( fConditionalMLEsAlt ) delete fConditionalMLEsAlt;
+ if( fFitInfo ) delete fFitInfo;
}
@@ -82,6 +91,14 @@ namespace RooStats {
else fConditionalMLEsAlt = NULL;
}
+ void StoreFitInfo(bool val = true) {
+ fStoreFitInfo = val;
+ }
+
+ const RooArgSet* GetFitInfo() const {
+ return fFitInfo;
+ }
+
protected:
// configure TestStatSampler for the Null run
int PreNullHook(RooArgSet *parameterPoint, double obsTestStat) const;
@@ -89,6 +106,9 @@ namespace RooStats {
// configure TestStatSampler for the Alt run
int PreAltHook(RooArgSet *parameterPoint, double obsTestStat) const;
+ void PreHook() const;
+ void PostHook() const;
+
protected:
// MLE inputs
const RooArgSet* fConditionalMLEsNull;
@@ -102,6 +122,10 @@ namespace RooStats {
int fNToysNullTail;
int fNToysAltTail;
+ private:
+ mutable RooArgSet* fFitInfo;
+ bool fStoreFitInfo;
+
protected:
ClassDef(FrequentistCalculator,1)
};
diff --git a/roofit/roostats/inc/RooStats/HypoTestCalculatorGeneric.h b/roofit/roostats/inc/RooStats/HypoTestCalculatorGeneric.h
index 51ea7a2db67..92006f5a18d 100644
--- a/roofit/roostats/inc/RooStats/HypoTestCalculatorGeneric.h
+++ b/roofit/roostats/inc/RooStats/HypoTestCalculatorGeneric.h
@@ -65,6 +65,7 @@ namespace RooStats {
virtual void SetNullModel(const ModelConfig &nullModel) { fNullModel = &nullModel; }
const RooAbsData * GetData(void) const { return fData; }
const ModelConfig* GetNullModel(void) const { return fNullModel; }
+ virtual const RooArgSet* GetFitInfo() const { return NULL; }
// set the model for the alternate hypothesis (S+B)
virtual void SetAlternateModel(const ModelConfig &altModel) { fAltModel = &altModel; }
const ModelConfig* GetAlternateModel(void) const { return fAltModel; }
@@ -80,6 +81,8 @@ namespace RooStats {
virtual int CheckHook(void) const { return 0; }
virtual int PreNullHook(RooArgSet* /*parameterPoint*/, double /*obsTestStat*/) const { return 0; }
virtual int PreAltHook(RooArgSet* /*parameterPoint*/, double /*obsTestStat*/) const { return 0; }
+ virtual void PreHook() const { }
+ virtual void PostHook() const { }
protected:
const ModelConfig *fAltModel;
diff --git a/roofit/roostats/inc/RooStats/HypoTestResult.h b/roofit/roostats/inc/RooStats/HypoTestResult.h
index 5cec5734b62..b7a7e19e3c8 100644
--- a/roofit/roostats/inc/RooStats/HypoTestResult.h
+++ b/roofit/roostats/inc/RooStats/HypoTestResult.h
@@ -109,6 +109,7 @@ namespace RooStats {
SamplingDistribution* GetAltDistribution(void) const { return fAltDistr; }
RooDataSet* GetNullDetailedOutput(void) const { return fNullDetailedOutput; }
RooDataSet* GetAltDetailedOutput(void) const { return fAltDetailedOutput; }
+ RooDataSet* GetFitInfo(void) const { return fFitInfo; }
Double_t GetTestStatisticData(void) const { return fTestStatisticData; }
const RooArgList* GetAllTestStatisticsData(void) const { return fAllTestStatisticsData; }
Bool_t HasTestStatisticData(void) const;
@@ -117,6 +118,7 @@ namespace RooStats {
void SetNullDistribution(SamplingDistribution *null);
void SetAltDetailedOutput(RooDataSet* d) { fAltDetailedOutput = d; }
void SetNullDetailedOutput(RooDataSet* d) { fNullDetailedOutput = d; }
+ void SetFitInfo(RooDataSet* d) { fFitInfo = d; }
void SetTestStatisticData(const Double_t tsd);
void SetAllTestStatisticsData(const RooArgList* tsd);
@@ -156,6 +158,7 @@ namespace RooStats {
SamplingDistribution *fAltDistr;
RooDataSet* fNullDetailedOutput;
RooDataSet* fAltDetailedOutput;
+ RooDataSet* fFitInfo;
Bool_t fPValueIsRightTail;
Bool_t fBackgroundIsAlt;
diff --git a/roofit/roostats/inc/RooStats/MinNLLTestStat.h b/roofit/roostats/inc/RooStats/MinNLLTestStat.h
new file mode 100755
index 00000000000..8ce6b8d095b
--- /dev/null
+++ b/roofit/roostats/inc/RooStats/MinNLLTestStat.h
@@ -0,0 +1,107 @@
+// @(#)root/roostats:$Id: MinNLLTestStat.h 43035 2012-02-16 16:48:57Z sven $
+// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
+// Additional Contributions: Giovanni Petrucciani
+/*************************************************************************
+ * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
+ * All rights reserved. *
+ * *
+ * For the licensing terms see $ROOTSYS/LICENSE. *
+ * For the list of contributors see $ROOTSYS/README/CREDITS. *
+ *************************************************************************/
+
+#ifndef ROOSTATS_MinNLLTestStat
+#define ROOSTATS_MinNLLTestStat
+
+//_________________________________________________
+/*
+BEGIN_HTML
+<p>
+MinNLLTestStat is an implementation of the TestStatistic interface that calculates the minimum value of the negative log likelihood
+function and returns it as a test statistic.
+Internaly it operates by delegating to a MinNLLTestStat object.
+</p>
+END_HTML
+*/
+//
+
+#ifndef ROOT_Rtypes
+#include "Rtypes.h"
+#endif
+
+#include <vector>
+
+#include "RooStats/RooStatsUtils.h"
+
+//#include "RooStats/DistributionCreator.h"
+#include "RooStats/SamplingDistribution.h"
+#include "RooStats/TestStatistic.h"
+
+#include "RooStats/RooStatsUtils.h"
+
+#include "RooRealVar.h"
+#include "RooProfileLL.h"
+#include "RooNLLVar.h"
+#include "RooMsgService.h"
+
+#include "RooMinuit.h"
+#include "RooMinimizer.h"
+#include "Math/MinimizerOptions.h"
+#include "TStopwatch.h"
+#include "ProfileLikelihoodTestStat.h"
+
+namespace RooStats {
+
+ class MinNLLTestStat : public TestStatistic{
+
+ public:
+ MinNLLTestStat() {
+ // Proof constructor. Do not use.
+ proflts = 0;
+ }
+ MinNLLTestStat(RooAbsPdf& pdf) {
+ proflts = new ProfileLikelihoodTestStat(pdf);
+ }
+
+ virtual ~MinNLLTestStat() {
+ delete proflts;
+ }
+
+ void SetOneSided(Bool_t flag=true) {proflts->SetOneSided(flag);}
+ void SetOneSidedDiscovery(Bool_t flag=true) {proflts->SetOneSidedDiscovery(flag);}
+ void SetReuseNLL(Bool_t flag) { proflts->SetReuseNLL(flag); }
+ void SetMinimizer(const char* minimizer){ proflts->SetMinimizer(minimizer); }
+ void SetStrategy(Int_t strategy){ proflts->SetStrategy(strategy); }
+ void SetTolerance(double tol){ proflts->SetTolerance(tol); }
+ void SetPrintLevel(Int_t printlevel){ proflts->SetPrintLevel(printlevel); }
+
+ // Main interface to evaluate the test statistic on a dataset
+ virtual Double_t Evaluate(RooAbsData& data, RooArgSet& paramsOfInterest) {
+ return proflts->EvaluateProfileLikelihood(1, data, paramsOfInterest); //find unconditional NLL minimum
+ }
+
+ virtual void EnableDetailedOutput( bool e=true ) { proflts->EnableDetailedOutput(e); }
+
+ virtual const RooArgSet* GetDetailedOutput(void) const {
+ // Returns detailed output. The value returned by this function is updated after each call to Evaluate().
+ // The returned RooArgSet contains the following:
+ // <ul>
+ // <li> the minimum nll, fitstatus and convergence quality for each fit </li>
+ // <li> for all non-constant parameters their value, error and pull </li>
+ // </ul>
+ return proflts->GetDetailedOutput();
+ }
+
+ virtual void SetVarName(const char* name) { proflts->SetVarName(name); }
+
+ virtual const TString GetVarName() const { return proflts->GetVarName(); }
+
+ private:
+ ProfileLikelihoodTestStat* proflts;
+
+ protected:
+ ClassDef(MinNLLTestStat,1) // implements the minimum NLL as a test statistic to be used with several tools
+ };
+}
+
+
+#endif
diff --git a/roofit/roostats/inc/RooStats/ProfileLikelihoodTestStat.h b/roofit/roostats/inc/RooStats/ProfileLikelihoodTestStat.h
index fdd0820e4c6..c97071c47bc 100644
--- a/roofit/roostats/inc/RooStats/ProfileLikelihoodTestStat.h
+++ b/roofit/roostats/inc/RooStats/ProfileLikelihoodTestStat.h
@@ -69,11 +69,6 @@ namespace RooStats {
fDetailedOutputEnabled = false;
fDetailedOutput = NULL;
- fUncML = new RooRealVar("uncondML","unconditional ML", 0.0);
- fFitStatus = new RooRealVar("fitStatus","fit status", 0.0);
- fCovQual = new RooRealVar("covQual","quality of covariance matrix", 0.0);
- fNumInvalidNLLEval = new RooRealVar("numInvalidNLLEval","number of invalid NLL evaluations", 0.0);
-
fVarName = "Profile Likelihood Ratio";
fReuseNll = false;
fMinimizer=::ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str();
@@ -92,11 +87,6 @@ namespace RooStats {
fDetailedOutputEnabled = false;
fDetailedOutput = NULL;
- fUncML = new RooRealVar("uncondML","unconditional ML", 0.0);
- fFitStatus = new RooRealVar("fitStatus","fit status", 0.0);
- fCovQual = new RooRealVar("covQual","quality of covariance matrix", 0.0);
- fNumInvalidNLLEval = new RooRealVar("numInvalidNLLEval","number of invalid NLL evaluations", 0.0);
-
fVarName = "Profile Likelihood Ratio";
fReuseNll = false;
fMinimizer=::ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str();
@@ -135,8 +125,20 @@ namespace RooStats {
// evaluate the profile likelihood ratio (type = 0) or the minimum of likelihood (type=1) or the conditional LL (type = 2)
virtual Double_t EvaluateProfileLikelihood(int type, RooAbsData &data, RooArgSet & paramsOfInterest);
- virtual void EnableDetailedOutput( bool e=true ) { fDetailedOutputEnabled = e; fDetailedOutput = NULL; }
- virtual const RooArgSet* GetDetailedOutput(void) const { return fDetailedOutput; }
+ virtual void EnableDetailedOutput( bool e=true, bool withErrorsAndPulls=false ) {
+ fDetailedOutputEnabled = e;
+ fDetailedOutputWithErrorsAndPulls = withErrorsAndPulls;
+ fDetailedOutput = NULL;
+ }
+ virtual const RooArgSet* GetDetailedOutput(void) const {
+ // Returns detailed output. The value returned by this function is updated after each call to Evaluate().
+ // The returned RooArgSet contains the following:
+ // <ul>
+ // <li> the minimum nll, fitstatus and convergence quality for each fit </li>
+ // <li> for each fit and for each non-constant parameter, the value, error and pull of the parameter are stored </li>
+ // </ul>
+ return fDetailedOutput;
+ }
virtual void SetVarName(const char* name) { fVarName = name; }
virtual const TString GetVarName() const {return fVarName;}
@@ -145,7 +147,7 @@ namespace RooStats {
private:
- double GetMinNLL(int& status);
+ RooFitResult* GetMinNLL();
private:
@@ -160,12 +162,8 @@ namespace RooStats {
// this will store a snapshot of the unconditional nuisance
// parameter fit.
bool fDetailedOutputEnabled;
- const RooArgSet* fDetailedOutput; //!
-
- RooRealVar* fUncML; //!
- RooRealVar* fFitStatus; //!
- RooRealVar* fCovQual; //!
- RooRealVar* fNumInvalidNLLEval; //!
+ bool fDetailedOutputWithErrorsAndPulls;
+ RooArgSet* fDetailedOutput; //!
TString fVarName;
diff --git a/roofit/roostats/inc/RooStats/RatioOfProfiledLikelihoodsTestStat.h b/roofit/roostats/inc/RooStats/RatioOfProfiledLikelihoodsTestStat.h
index 38d1ae225fe..4a13af3e6ea 100644
--- a/roofit/roostats/inc/RooStats/RatioOfProfiledLikelihoodsTestStat.h
+++ b/roofit/roostats/inc/RooStats/RatioOfProfiledLikelihoodsTestStat.h
@@ -47,7 +47,9 @@ class RatioOfProfiledLikelihoodsTestStat: public TestStatistic {
fNullProfile(),
fAltProfile(),
fAltPOI(NULL),
- fSubtractMLE(true)
+ fSubtractMLE(true),
+ fDetailedOutputEnabled(false),
+ fDetailedOutput(NULL)
{
// Proof constructor. Don't use.
}
@@ -56,7 +58,9 @@ class RatioOfProfiledLikelihoodsTestStat: public TestStatistic {
const RooArgSet* altPOI=0) :
fNullProfile(nullPdf),
fAltProfile(altPdf),
- fSubtractMLE(true)
+ fSubtractMLE(true),
+ fDetailedOutputEnabled(false),
+ fDetailedOutput(NULL)
{
/*
Calculates the ratio of profiled likelihoods.
@@ -94,6 +98,7 @@ class RatioOfProfiledLikelihoodsTestStat: public TestStatistic {
//__________________________________________
~RatioOfProfiledLikelihoodsTestStat(void) {
if(fAltPOI) delete fAltPOI;
+ if(fDetailedOutput) delete fDetailedOutput;
}
//__________________________________________
@@ -125,13 +130,34 @@ class RatioOfProfiledLikelihoodsTestStat: public TestStatistic {
int type = (fSubtractMLE) ? 0 : 2;
-
+
// null
double nullNLL = fNullProfile.EvaluateProfileLikelihood(type, data, nullParamsOfInterest);
+ const RooArgSet *nullset = fNullProfile.GetDetailedOutput();
// alt
double altNLL = fAltProfile.EvaluateProfileLikelihood(type, data, *fAltPOI);
-
+ const RooArgSet *altset = fAltProfile.GetDetailedOutput();
+
+ if (fDetailedOutput != NULL) {
+ delete fDetailedOutput;
+ fDetailedOutput = NULL;
+ }
+ if (fDetailedOutputEnabled) {
+ fDetailedOutput = new RooArgSet();
+ RooRealVar* var(0);
+ for(TIterator *it = nullset->createIterator();(var = dynamic_cast<RooRealVar*>(it->Next()));) {
+ RooRealVar* cloneVar = new RooRealVar(TString::Format("nullprof_%s", var->GetName()),
+ TString::Format("%s for null", var->GetTitle()), var->getVal());
+ fDetailedOutput->addOwned(*cloneVar);
+ }
+ for(TIterator *it = altset->createIterator();(var = dynamic_cast<RooRealVar*>(it->Next()));) {
+ RooRealVar* cloneVar = new RooRealVar(TString::Format("altprof_%s", var->GetName()),
+ TString::Format("%s for null", var->GetTitle()), var->getVal());
+ fDetailedOutput->addOwned(*cloneVar);
+ }
+ }
+
/*
// set variables back to where they were
nullParamsOfInterest = *saveNullPOI;
@@ -143,8 +169,13 @@ class RatioOfProfiledLikelihoodsTestStat: public TestStatistic {
return nullNLL -altNLL;
}
-
- static void SetAlwaysReuseNLL(Bool_t flag) { fgAlwaysReuseNll = flag ; }
+ virtual void EnableDetailedOutput( bool e=true ) {
+ fDetailedOutputEnabled = e;
+ fNullProfile.EnableDetailedOutput(fDetailedOutputEnabled);
+ fAltProfile.EnableDetailedOutput(fDetailedOutputEnabled);
+ }
+
+ static void SetAlwaysReuseNLL(Bool_t flag) { fgAlwaysReuseNll = flag ; }
void SetReuseNLL(Bool_t flag) {
fNullProfile.SetReuseNLL(flag);
@@ -168,6 +199,16 @@ class RatioOfProfiledLikelihoodsTestStat: public TestStatistic {
fAltProfile.SetPrintLevel(printLevel);
}
+ virtual const RooArgSet* GetDetailedOutput(void) const {
+ // Returns detailed output. The value returned by this function is updated after each call to Evaluate().
+ // The returned RooArgSet contains the following for the alternative and null hypotheses:
+ // <ul>
+ // <li> the minimum nll, fitstatus and convergence quality for each fit </li>
+ // <li> for each fit and for each non-constant parameter, the value, error and pull of the parameter are stored </li>
+ // </ul>
+ return fDetailedOutput;
+ }
+
virtual const TString GetVarName() const { return "log(L(#mu_{1},#hat{#nu}_{1}) / L(#mu_{0},#hat{#nu}_{0}))"; }
@@ -178,10 +219,14 @@ class RatioOfProfiledLikelihoodsTestStat: public TestStatistic {
private:
ProfileLikelihoodTestStat fNullProfile;
ProfileLikelihoodTestStat fAltProfile;
+
RooArgSet* fAltPOI;
Bool_t fSubtractMLE;
static Bool_t fgAlwaysReuseNll ;
+ bool fDetailedOutputEnabled;
+ RooArgSet* fDetailedOutput;
+
protected:
ClassDef(RatioOfProfiledLikelihoodsTestStat,3)
diff --git a/roofit/roostats/inc/RooStats/RooStatsUtils.h b/roofit/roostats/inc/RooStats/RooStatsUtils.h
index ca6bc804272..39875ea4e4d 100644
--- a/roofit/roostats/inc/RooStats/RooStatsUtils.h
+++ b/roofit/roostats/inc/RooStats/RooStatsUtils.h
@@ -16,6 +16,10 @@
#include "TMath.h"
#endif
+#ifndef ROOT_TTree
+#include "TTree.h"
+#endif
+
#ifndef ROOT_Math_DistFuncMathCore
#include"Math/DistFuncMathCore.h"
#endif
@@ -26,6 +30,7 @@
#include "TIterator.h"
#include "RooStats/ModelConfig.h"
#include "RooProdPdf.h"
+#include "RooDataSet.h"
namespace RooStats {
@@ -101,6 +106,9 @@ namespace RooStats {
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name);
RooAbsPdf * MakeNuisancePdf(const RooStats::ModelConfig &model, const char *name);
+
+ // Create a TTree with the given name and description. All RooRealVars in the RooDataSet are represented as branches that contain values of type Double_t.
+ TTree* GetAsTTree(TString name, TString desc, const RooDataSet& data);
}
diff --git a/roofit/roostats/inc/RooStats/SamplingDistribution.h b/roofit/roostats/inc/RooStats/SamplingDistribution.h
index 96c71fdc9f2..5261852abf0 100644
--- a/roofit/roostats/inc/RooStats/SamplingDistribution.h
+++ b/roofit/roostats/inc/RooStats/SamplingDistribution.h
@@ -41,7 +41,7 @@ namespace RooStats {
SamplingDistribution(const char *name,const char *title, const char * varName = 0);
- SamplingDistribution(const char *name,const char *title, RooDataSet& dataSet, const char * varName = 0);
+ SamplingDistribution(const char *name,const char *title, RooDataSet& dataSet, const char * columnName = 0, const char * varName = 0);
// Default constructor for SamplingDistribution
SamplingDistribution();
diff --git a/roofit/roostats/inc/RooStats/TestStatistic.h b/roofit/roostats/inc/RooStats/TestStatistic.h
index ac0591979ce..495db265635 100644
--- a/roofit/roostats/inc/RooStats/TestStatistic.h
+++ b/roofit/roostats/inc/RooStats/TestStatistic.h
@@ -47,7 +47,7 @@ class TestStatistic {
// Defines the sign convention of the test statistic. Overwrite function if necessary.
virtual bool PValueIsRightTail(void) const { return true; }
- // return detailed output: for fits this can be pulls, processing time, ...
+ // return detailed output: for fits this can be pulls, processing time, ... The returned pointer will not loose validity until another call to Evaluate.
virtual const RooArgSet* GetDetailedOutput() const { return NULL; }
protected:
diff --git a/roofit/roostats/inc/RooStats/ToyMCSampler.h b/roofit/roostats/inc/RooStats/ToyMCSampler.h
index f06b7e64429..31460eb18e1 100644
--- a/roofit/roostats/inc/RooStats/ToyMCSampler.h
+++ b/roofit/roostats/inc/RooStats/ToyMCSampler.h
@@ -188,7 +188,7 @@ class ToyMCSampler: public TestStatSampler {
if( fParametersForTestStat ) delete fParametersForTestStat;
fParametersForTestStat = (const RooArgSet*)nullpoi.snapshot();
}
-
+
virtual void SetPdf(RooAbsPdf& pdf) { fPdf = &pdf; ClearCache(); }
// How to randomize the prior. Set to NULL to deactivate randomization.
@@ -208,11 +208,14 @@ class ToyMCSampler: public TestStatSampler {
// Set the TestStatistic (want the argument to be a function of the data & parameter points
virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i) {
- if( fTestStatistics.size() <= i ) {
+ if( fTestStatistics.size() < i ) {
oocoutE((TObject*)NULL,InputArguments) << "Cannot set test statistic for this index." << endl;
return;
}
- fTestStatistics[i] = testStatistic;
+ if( fTestStatistics.size() == i)
+ fTestStatistics.push_back(testStatistic);
+ else
+ fTestStatistics[i] = testStatistic;
}
virtual void SetTestStatistic(TestStatistic *t) { return SetTestStatistic(t,0); }
diff --git a/roofit/roostats/inc/RooStats/ToyMCStudy.h b/roofit/roostats/inc/RooStats/ToyMCStudy.h
index 876eef9f802..d872f892487 100644
--- a/roofit/roostats/inc/RooStats/ToyMCStudy.h
+++ b/roofit/roostats/inc/RooStats/ToyMCStudy.h
@@ -87,6 +87,7 @@ class ToyMCPayload : public TNamed {
ToyMCPayload() {
// proof constructor, do not use
+ fDataSet = NULL;
}
ToyMCPayload(RooDataSet* sd)
diff --git a/roofit/roostats/src/FrequentistCalculator.cxx b/roofit/roostats/src/FrequentistCalculator.cxx
index 4bc3118b522..b5b095ed990 100644
--- a/roofit/roostats/src/FrequentistCalculator.cxx
+++ b/roofit/roostats/src/FrequentistCalculator.cxx
@@ -15,13 +15,26 @@ MLEs.
#include "RooStats/FrequentistCalculator.h"
#include "RooStats/ToyMCSampler.h"
+#include "RooMinuit.h"
+#include "RooProfileLL.h"
ClassImp(RooStats::FrequentistCalculator)
using namespace RooStats;
+void FrequentistCalculator::PreHook() const {
+ if (fFitInfo != NULL) {
+ delete fFitInfo;
+ fFitInfo = NULL;
+ }
+ if (fStoreFitInfo) {
+ fFitInfo = new RooArgSet();
+ }
+}
+void FrequentistCalculator::PostHook() const {
+}
int FrequentistCalculator::PreNullHook(RooArgSet *parameterPoint, double obsTestStat) const {
@@ -31,35 +44,55 @@ int FrequentistCalculator::PreNullHook(RooArgSet *parameterPoint, double obsTest
RooArgSet * allParams = fNullModel->GetPdf()->getParameters(*fData);
RemoveConstantParameters(allParams);
- oocoutI((TObject*)0,InputArguments) << "Profiling conditional MLEs for Null." << endl;
// note: making nll or profile class variables can only be done in the constructor
// as all other hooks are const (which has to be because GetHypoTest is const). However,
// when setting it only in constructor, they would have to be changed every time SetNullModel
// or SetAltModel is called. Simply put, converting them into class variables breaks
// encapsulation.
+ bool doProfile = true;
RooArgSet allButNuisance(*allParams);
- if( fNullModel->GetNuisanceParameters() ) allButNuisance.remove(*fNullModel->GetNuisanceParameters());
- if( fConditionalMLEsNull ) {
- oocoutI((TObject*)0,InputArguments) << "Using given conditional MLEs for Null." << endl;
- *allParams = *fConditionalMLEsNull;
- allButNuisance.add( *fConditionalMLEsNull );
+ if( fNullModel->GetNuisanceParameters() ) {
+ allButNuisance.remove(*fNullModel->GetNuisanceParameters());
+ if( fConditionalMLEsNull ) {
+ oocoutI((TObject*)0,InputArguments) << "Using given conditional MLEs for Null." << endl;
+ *allParams = *fConditionalMLEsNull;
+ allButNuisance.add( *fConditionalMLEsNull );
+ if (fNullModel->GetNuisanceParameters()) {
+ RooArgSet remain(*fNullModel->GetNuisanceParameters());
+ remain.remove(*fConditionalMLEsNull,true,true);
+ if( remain.getSize() == 0 ) doProfile = false;
+ }
+ }
+ }else{
+ doProfile = false;
}
- RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
- RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
-
- RooAbsReal* nll = fNullModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(kFALSE), RooFit::Constrain(*allParams));
- RooAbsReal* profile = nll->createProfile(allButNuisance);
- profile->getVal(); // this will do fit and set nuisance parameters to profiled values
+ if (doProfile) {
+ oocoutI((TObject*)0,InputArguments) << "Profiling conditional MLEs for Null." << endl;
+ RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
+ RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
+
+ RooAbsReal* nll = fNullModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(kFALSE), RooFit::Constrain(*allParams));
+ RooProfileLL* profile = dynamic_cast<RooProfileLL*>(nll->createProfile(allButNuisance));
+ profile->getVal(); // this will do fit and set nuisance parameters to profiled values
+
+ // Hack to extract a RooFitResult
+ if (fStoreFitInfo) {
+ RooFitResult *result = profile->minuit()->save();
+ fFitInfo->addOwned(*DetailedOutputAggregator::GetAsArgSet(result, "fitNull_"));
+ delete result;
+ }
+
+ delete profile;
+ delete nll;
+ RooMsgService::instance().setGlobalKillBelow(msglevel);
+ }
+
// add nuisance parameters to parameter point
if(fNullModel->GetNuisanceParameters())
parameterPoint->add(*fNullModel->GetNuisanceParameters());
- delete profile;
- delete nll;
-
delete allParams;
- RooMsgService::instance().setGlobalKillBelow(msglevel);
@@ -103,29 +136,49 @@ int FrequentistCalculator::PreAltHook(RooArgSet *parameterPoint, double obsTestS
RooArgSet * allParams = fAltModel->GetPdf()->getParameters(*fData);
RemoveConstantParameters(allParams);
- oocoutI((TObject*)0,InputArguments) << "Profiling conditional MLEs for Alt." << endl;
-
+ bool doProfile = true;
RooArgSet allButNuisance(*allParams);
- if( fAltModel->GetNuisanceParameters() ) allButNuisance.remove(*fAltModel->GetNuisanceParameters());
- if( fConditionalMLEsAlt ) {
- oocoutI((TObject*)0,InputArguments) << "Using given conditional MLEs for Alt." << endl;
- *allParams = *fConditionalMLEsAlt;
- allButNuisance.add( *fConditionalMLEsAlt );
+ if( fAltModel->GetNuisanceParameters() ) {
+ allButNuisance.remove(*fAltModel->GetNuisanceParameters());
+ if( fConditionalMLEsAlt ) {
+ oocoutI((TObject*)0,InputArguments) << "Using given conditional MLEs for Alt." << endl;
+ *allParams = *fConditionalMLEsAlt;
+ allButNuisance.add( *fConditionalMLEsAlt );
+ if (fAltModel->GetNuisanceParameters()) {
+ RooArgSet remain(*fAltModel->GetNuisanceParameters());
+ remain.remove(*fConditionalMLEsAlt,true,true);
+ if( remain.getSize() == 0 ) doProfile = false;
+ }
+ }
+ }else{
+ doProfile = false;
+ }
+ if (doProfile) {
+ oocoutI((TObject*)0,InputArguments) << "Profiling conditional MLEs for Alt." << endl;
+ RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
+ RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
+
+ RooAbsReal* nll = fAltModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(kFALSE), RooFit::Constrain(*allParams));
+ RooProfileLL* profile = dynamic_cast<RooProfileLL*>(nll->createProfile(allButNuisance));
+ profile->getVal(); // this will do fit and set nuisance parameters to profiled values
+
+ // Hack to extract a RooFitResult
+ if (fStoreFitInfo) {
+ RooFitResult *result = profile->minuit()->save();
+ fFitInfo->addOwned(*DetailedOutputAggregator::GetAsArgSet(result, "fitAlt_"));
+ delete result;
+ }
+
+ delete profile;
+ delete nll;
+ RooMsgService::instance().setGlobalKillBelow(msglevel);
}
- RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
- RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
- RooAbsReal* nll = fAltModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(kFALSE), RooFit::Constrain(*allParams));
- RooAbsReal* profile = nll->createProfile(allButNuisance);
- profile->getVal(); // this will do fit and set nuisance parameters to profiled values
// add nuisance parameters to parameter point
if(fAltModel->GetNuisanceParameters())
parameterPoint->add(*fAltModel->GetNuisanceParameters());
- delete profile;
- delete nll;
-
+
delete allParams;
- RooMsgService::instance().setGlobalKillBelow(msglevel);
diff --git a/roofit/roostats/src/HypoTestCalculatorGeneric.cxx b/roofit/roostats/src/HypoTestCalculatorGeneric.cxx
index 772ac84a02b..ecb39309b05 100644
--- a/roofit/roostats/src/HypoTestCalculatorGeneric.cxx
+++ b/roofit/roostats/src/HypoTestCalculatorGeneric.cxx
@@ -92,6 +92,7 @@ HypoTestResult* HypoTestCalculatorGeneric::GetHypoTest() const {
// - nuisance parameters are floating, so they do float in test statistic
// initial setup
+ PreHook();
const_cast<ModelConfig*>(fNullModel)->GuessObsAndNuisance(*fData);
const_cast<ModelConfig*>(fAltModel)->GuessObsAndNuisance(*fData);
@@ -188,12 +189,20 @@ HypoTestResult* HypoTestCalculatorGeneric::GetHypoTest() const {
res->SetNullDetailedOutput( detOut_null );
res->SetAllTestStatisticsData( allTS );
+ const RooArgSet *aset = GetFitInfo();
+ if (aset != NULL) {
+ RooDataSet *dset = new RooDataSet(NULL, NULL, *aset);
+ dset->add(*aset);
+ res->SetFitInfo( dset );
+ }
+
*bothParams = *saveAll;
delete bothParams;
delete saveAll;
delete altParams;
delete nullParams;
delete nullSnapshot;
+ PostHook();
return res;
}
diff --git a/roofit/roostats/src/HypoTestResult.cxx b/roofit/roostats/src/HypoTestResult.cxx
index 09f59465576..f11b045c399 100644
--- a/roofit/roostats/src/HypoTestResult.cxx
+++ b/roofit/roostats/src/HypoTestResult.cxx
@@ -60,7 +60,7 @@ HypoTestResult::HypoTestResult(const char* name) :
fTestStatisticData(NaN),
fAllTestStatisticsData(NULL),
fNullDistr(NULL), fAltDistr(NULL),
- fNullDetailedOutput(NULL), fAltDetailedOutput(NULL),
+ fNullDetailedOutput(NULL), fAltDetailedOutput(NULL), fFitInfo(NULL),
fPValueIsRightTail(kTRUE),
fBackgroundIsAlt(kFALSE)
{
@@ -76,7 +76,7 @@ HypoTestResult::HypoTestResult(const char* name, Double_t nullp, Double_t altp)
fTestStatisticData(NaN),
fAllTestStatisticsData(NULL),
fNullDistr(NULL), fAltDistr(NULL),
- fNullDetailedOutput(NULL), fAltDetailedOutput(NULL),
+ fNullDetailedOutput(NULL), fAltDetailedOutput(NULL), fFitInfo(NULL),
fPValueIsRightTail(kTRUE),
fBackgroundIsAlt(kFALSE)
{
@@ -91,7 +91,7 @@ HypoTestResult::HypoTestResult(const HypoTestResult& other) :
fTestStatisticData(NaN),
fAllTestStatisticsData(NULL),
fNullDistr(NULL), fAltDistr(NULL),
- fNullDetailedOutput(NULL), fAltDetailedOutput(NULL),
+ fNullDetailedOutput(NULL), fAltDetailedOutput(NULL), fFitInfo(NULL),
fPValueIsRightTail( other.GetPValueIsRightTail() ),
fBackgroundIsAlt( other.GetBackGroundIsAlt() )
{
@@ -142,6 +142,12 @@ void HypoTestResult::Append(const HypoTestResult* other) {
if( other->GetAltDetailedOutput() ) fAltDetailedOutput = new RooDataSet( *other->GetAltDetailedOutput() );
}
+ if( fFitInfo ) {
+ if( other->GetFitInfo() ) fFitInfo->append( *other->GetFitInfo() );
+ }else{
+ if( other->GetFitInfo() ) fFitInfo = new RooDataSet( *other->GetFitInfo() );
+ }
+
// if no data is present use the other HypoTestResult's data
if(IsNaN(fTestStatisticData)) fTestStatisticData = other->GetTestStatisticData();
diff --git a/roofit/roostats/src/ProfileLikelihoodTestStat.cxx b/roofit/roostats/src/ProfileLikelihoodTestStat.cxx
index 736d5a0b075..4a21253f51d 100644
--- a/roofit/roostats/src/ProfileLikelihoodTestStat.cxx
+++ b/roofit/roostats/src/ProfileLikelihoodTestStat.cxx
@@ -10,13 +10,16 @@
*************************************************************************/
#include "RooStats/ProfileLikelihoodTestStat.h"
+#include "RooFitResult.h"
+#include "RooPullVar.h"
+#include "RooStats/DetailedOutputAggregator.h"
#include "RooProfileLL.h"
#include "RooNLLVar.h"
#include "RooMsgService.h"
#include "RooMinimizer.h"
#include "RooArgSet.h"
-#include "RooAbsData.h"
+#include "RooDataSet.h"
#include "TStopwatch.h"
#include "RooStats/RooStatsUtils.h"
@@ -33,6 +36,13 @@ Double_t RooStats::ProfileLikelihoodTestStat::EvaluateProfileLikelihood(int type
cout << "problem with data" << endl;
return 0 ;
}
+ if( fDetailedOutputEnabled and fDetailedOutput ) {
+ delete fDetailedOutput;
+ fDetailedOutput = 0;
+ }
+ if( fDetailedOutputEnabled and !fDetailedOutput ) {
+ fDetailedOutput = new RooArgSet();
+ }
//data.Print("V");
@@ -98,33 +108,17 @@ Double_t RooStats::ProfileLikelihoodTestStat::EvaluateProfileLikelihood(int type
if (type != 2) {
// minimize and count eval errors
fNll->clearEvalErrorLog();
- uncondML = GetMinNLL(statusD);
- int invalidNLLEval = fNll->numEvalErrors();
+ RooFitResult* result = GetMinNLL();
+ uncondML = result->minNll();
+ statusD = result->status();
// get best fit value for one-sided interval
- if (firstPOI) fit_favored_mu = attachedSet->getRealValue(firstPOI->GetName()) ;
-
- // save this snapshot
- if( fDetailedOutputEnabled ) {
- if( fDetailedOutput ) delete fDetailedOutput;
-
- RooArgSet detOut( *attachedSet );
- RooStats::RemoveConstantParameters( &detOut );
-
- // monitor a few more variables
- fUncML->setVal( uncondML ); detOut.add( *fUncML );
- fFitStatus->setVal( statusD ); detOut.add( *fFitStatus );
- //fCovQual->setVal( covQual ); detOut.add( *fCovQual );
- fNumInvalidNLLEval->setVal( invalidNLLEval ); detOut.add( *fNumInvalidNLLEval );
-
- // store it
- fDetailedOutput = (const RooArgSet*)detOut.snapshot();
-
- cout << endl << "STORING THIS AS DETAILED OUTPUT:" << endl;
- fDetailedOutput->Print("v");
- cout << endl;
- }
+ if (firstPOI) fit_favored_mu = attachedSet->getRealValue(firstPOI->GetName()) ;
+ // save this snapshot
+ if( fDetailedOutputEnabled )
+ fDetailedOutput->addOwned(*DetailedOutputAggregator::GetAsArgSet(result, "fitUncond_", fDetailedOutputWithErrorsAndPulls));
+ delete result;
}
tsw.Stop();
double fitTime1 = tsw.CpuTime();
@@ -171,7 +165,13 @@ Double_t RooStats::ProfileLikelihoodTestStat::EvaluateProfileLikelihood(int type
condML = fNll->getVal();
}
else {
- condML = GetMinNLL(statusN);
+ fNll->clearEvalErrorLog();
+ RooFitResult* result = GetMinNLL();
+ condML = result->minNll();
+ statusN = result->status();
+ if( fDetailedOutputEnabled )
+ fDetailedOutput->addOwned(*DetailedOutputAggregator::GetAsArgSet(result, "fitCond_", fDetailedOutputWithErrorsAndPulls));
+ delete result;
}
}
@@ -227,7 +227,7 @@ Double_t RooStats::ProfileLikelihoodTestStat::EvaluateProfileLikelihood(int type
}
-double RooStats::ProfileLikelihoodTestStat::GetMinNLL(int& status) {
+RooFitResult* RooStats::ProfileLikelihoodTestStat::GetMinNLL() {
//find minimum of NLL using RooMinimizer
RooMinimizer minim(*fNll);
@@ -241,6 +241,7 @@ double RooStats::ProfileLikelihoodTestStat::GetMinNLL(int& status) {
TString minimizer = fMinimizer;
TString algorithm = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo();
if (algorithm == "Migrad") algorithm = "Minimize"; // prefer to use Minimize instead of Migrad
+ int status;
for (int tries = 1, maxtries = 4; tries <= maxtries; ++tries) {
status = minim.minimize(minimizer,algorithm);
if (status%1000 == 0) { // ignore erros from Improve
@@ -264,8 +265,7 @@ double RooStats::ProfileLikelihoodTestStat::GetMinNLL(int& status) {
}
}
- double val = fNll->getVal();
+ //how to get cov quality faster?
+ return minim.save();
//minim.optimizeConst(false);
-
- return val;
}
diff --git a/roofit/roostats/src/RooStatsUtils.cxx b/roofit/roostats/src/RooStatsUtils.cxx
index 9f40cca2b5a..5510148a0cb 100644
--- a/roofit/roostats/src/RooStatsUtils.cxx
+++ b/roofit/roostats/src/RooStatsUtils.cxx
@@ -20,6 +20,7 @@
NamespaceImp(RooStats)
#endif
+#include "TTree.h"
#include "RooProdPdf.h"
#include "RooSimultaneous.h"
@@ -88,9 +89,88 @@ namespace RooStats {
return MakeNuisancePdf(*model.GetPdf(), *model.GetObservables(), name);
}
+ // Helper class for GetAsTTree
+ class BranchStore {
+ public:
+ std::map<TString, Double_t> varVals;
+ double inval;
+ BranchStore(const vector <TString> ¶ms = vector <TString>(), double _inval = -999.) {
+ inval = _inval;
+ for(unsigned int i = 0;i<params.size();i++)
+ varVals[params[i]] = _inval;
+ }
+ void AssignToTTree(TTree &myTree) {
+ for(std::map<TString, Double_t>::iterator it = varVals.begin();it!=varVals.end();it++) {
+ const TString& name = it->first;
+ myTree.Branch( name, &varVals[name], TString::Format("%s/D", name.Data()));
+ }
+ }
+ void ResetValues() {
+ for(std::map<TString, Double_t>::iterator it = varVals.begin();it!=varVals.end();it++) {
+ const TString& name = it->first;
+ varVals[name] = inval;
+ }
+ }
+ };
+ BranchStore* CreateBranchStore(const RooDataSet& data) {
+ if (data.numEntries() == 0) {
+ return new BranchStore;
+ }
+ vector <TString> V;
+ const RooArgSet* aset = data.get(0);
+ RooAbsArg *arg(0);
+ TIterator *it = aset->createIterator();
+ for(;(arg = dynamic_cast<RooAbsArg*>(it->Next()));) {
+ RooRealVar *rvar = dynamic_cast<RooRealVar*>(arg);
+ if (rvar == NULL)
+ continue;
+ V.push_back(rvar->GetName());
+ if (rvar->hasAsymError()) {
+ V.push_back(TString::Format("%s_errlo", rvar->GetName()));
+ V.push_back(TString::Format("%s_errhi", rvar->GetName()));
+ }
+ else if (rvar->hasError()) {
+ V.push_back(TString::Format("%s_err", rvar->GetName()));
+ }
+ }
+ delete it;
+ return new BranchStore(V);
+ }
+
+ void FillTree(TTree &myTree, const RooDataSet &data) {
+ BranchStore *bs = CreateBranchStore(data);
+ bs->AssignToTTree(myTree);
+
+ for(int entry = 0;entry<data.numEntries();entry++) {
+ bs->ResetValues();
+ const RooArgSet* aset = data.get(entry);
+ RooAbsArg *arg(0);
+ TIterator *it = aset->createIterator();
+ for(;(arg = dynamic_cast<RooAbsArg*>(it->Next()));) {
+ RooRealVar *rvar = dynamic_cast<RooRealVar*>(arg);
+ if (rvar == NULL)
+ continue;
+ bs->varVals[rvar->GetName()] = rvar->getValV();
+ if (rvar->hasAsymError()) {
+ bs->varVals[TString::Format("%s_errlo", rvar->GetName())] = rvar->getAsymErrorLo();
+ bs->varVals[TString::Format("%s_errhi", rvar->GetName())] = rvar->getAsymErrorHi();
+ }
+ else if (rvar->hasError()) {
+ bs->varVals[TString::Format("%s_err", rvar->GetName())] = rvar->getError();
+ }
+ }
+ delete it;
+ myTree.Fill();
+ }
+ }
+ TTree * GetAsTTree(TString name, TString desc, const RooDataSet& data) {
+ TTree* myTree = new TTree(name, desc);
+ FillTree(*myTree, data);
+ return myTree;
+ }
}
diff --git a/roofit/roostats/src/SamplingDistribution.cxx b/roofit/roostats/src/SamplingDistribution.cxx
index aa5e9e933f5..aab87c31da7 100644
--- a/roofit/roostats/src/SamplingDistribution.cxx
+++ b/roofit/roostats/src/SamplingDistribution.cxx
@@ -23,6 +23,10 @@ a weighted set of points (eg. for the FFT method).
The class supports merging.
*/
+#ifndef ROO_MSG_SERVICE
+#include "RooMsgService.h"
+#endif
+
#include "RooStats/SamplingDistribution.h"
#include "RooNumber.h"
#include "TMath.h"
@@ -80,6 +84,7 @@ SamplingDistribution::SamplingDistribution(
const char *name,
const char *title,
RooDataSet& dataSet,
+ const char * _columnName,
const char * varName
) : TNamed(name, title) {
// Creates a SamplingDistribution from a RooDataSet for debugging
@@ -93,14 +98,30 @@ SamplingDistribution::SamplingDistribution(
// If varName is not given, the first variable will be used.
// This is useful mostly for RooDataSets with only one observable.
- fVarName = varName;
- if(fVarName.Length() == 0) {
+ // check there are any meaningful entries in the given dataset
+ if( dataSet.numEntries() == 0 || !dataSet.get()->first() ) {
+ if( varName ) fVarName = varName;
+ return;
+ }
+
+ TString columnName( _columnName );
+
+ if( !columnName.Length() ) {
+ columnName.Form( "%s_TS0", name );
+ if( !dataSet.get()->find(columnName) ) {
+ columnName = dataSet.get()->first()->GetName();
+ }
+ }
+
+ if( !varName ) {
// no leak. none of these transfers ownership.
- fVarName = dataSet.get()->first()->GetName();
+ fVarName = (*dataSet.get())[columnName].GetTitle();
+ }else{
+ fVarName = varName;
}
for(Int_t i=0; i < dataSet.numEntries(); i++) {
- fSamplingDist.push_back(dataSet.get(i)->getRealValue(fVarName));
+ fSamplingDist.push_back(dataSet.get(i)->getRealValue(columnName));
fSampleWeights.push_back(dataSet.weight());
}
}
diff --git a/roofit/roostats/src/ToyMCSampler.cxx b/roofit/roostats/src/ToyMCSampler.cxx
index a038d2db5f6..d415b719ed1 100644
--- a/roofit/roostats/src/ToyMCSampler.cxx
+++ b/roofit/roostats/src/ToyMCSampler.cxx
@@ -29,6 +29,7 @@
#include "RooStudyManager.h"
#include "RooStats/ToyMCStudy.h"
+#include "RooStats/DetailedOutputAggregator.h"
#include "RooSimultaneous.h"
#include "TMath.h"
@@ -126,7 +127,7 @@ Bool_t ToyMCSampler::fgAlwaysUseMultiGen = kFALSE ;
-ToyMCSampler::ToyMCSampler() : fSamplingDistName("temp"), fNToys(1)
+ToyMCSampler::ToyMCSampler() : fSamplingDistName("SD"), fNToys(1)
{
// Proof constructor. Do not use.
@@ -315,57 +316,17 @@ RooDataSet* ToyMCSampler::GetSamplingDistributionsSingleWorker(RooArgSet& paramP
RooArgSet *saveAll = (RooArgSet*) allVars->snapshot();
- // bug workaround: evaluate all test statistics once so that they have their detailed
- // output generated ... but only do that if there is no detailedOutput present already
- if( fTestStatistics[0] && !fTestStatistics[0]->GetDetailedOutput() ) {
- RooAbsData* dummyData = GenerateToyData(*paramPoint);
- EvaluateAllTestStatistics(*dummyData, *paramPoint);
- }
+ DetailedOutputAggregator detOutAgg;
- RooDataSet* outputs = new RooDataSet(
- fSamplingDistName.c_str(),fSamplingDistName.c_str(),
- RooArgSet( *(new RooRealVar("weight","weight",1.0)), "tmpSet" ),
- "weight"
- );
vector<TString> namesOfTSColumns;
+ RooArgSet* allVarsToBeSaved = new RooArgSet("tsValues");
for( unsigned int i=0; i < fTestStatistics.size(); i++ ) {
TString name( TString::Format("%s_TS%d", fSamplingDistName.c_str(), i) );
namesOfTSColumns.push_back( name );
-
- // Add a column for the test statistic value.
- // Its name is "NameOfSamplingDistribution_i" where is the index of
- // the snapshot/test statistic.
- // Title comes from the TestStatistic.
- outputs->addColumn( *(new RooRealVar(name, fTestStatistics[i]->GetVarName(), -1.0)) );
-
- // add columns for detailed output if this TS supplies them
- const RooArgSet* detOut = fTestStatistics[i]->GetDetailedOutput();
- if( detOut ) {
- // rename all variables to "NameOfSamplingDistribution_i_detailedOutputVariable"
- // for "namespacing" in the columns of the data set.
- // Titles remain the same.
- //cout << "Adding det out for TS " << i << endl;
- TIterator* iter = detOut->createIterator();
- while(RooAbsArg* v = dynamic_cast<RooAbsArg*>( iter->Next() ) ) {
- TString origName( v->GetName() );
- v->SetName( TString::Format( "%s_%s", name.Data(), v->GetName() ) );
- //cout << "adding variable: " << v->GetName() << endl;
- outputs->addColumn( *v );
- v->SetName( origName );
- }
- }
+ allVarsToBeSaved->addOwned( *new RooRealVar(name, fTestStatistics[i]->GetVarName(), -1) );
}
-
- // retrieve the set of variables to be stored and print some info
- RooArgSet* allVarsToBeSaved = new RooArgSet( *outputs->get() );
-// oocoutI((TObject*)NULL, InputArguments) << endl;
-// oocoutI((TObject*)NULL, InputArguments) << "The variables that will be stored for each toy:" << endl;
-// allVarsToBeSaved->Print();
-// oocoutI((TObject*)NULL, InputArguments) << endl;
-
-
-
+
// counts the number of toys in the limits set for adaptive sampling
// (taking weights into account; always on first test statistic)
Double_t toysInTails = 0.0;
@@ -385,53 +346,46 @@ RooDataSet* ToyMCSampler::GetSamplingDistributionsSingleWorker(RooArgSet& paramP
// TODO: change this treatment to keep track of all values so that the threshold
// for adaptive sampling is counted for all distributions and not just the
// first one.
- Double_t valueFirst = -1.0, weight = -1.0;
+ Double_t valueFirst = -999.0, weight = 1.0;
// set variables to requested parameter point
*allVars = *saveAll; // important for example for SimpleLikelihoodRatioTestStat
*allVars = *fParametersForTestStat;
-
+
RooAbsData* toydata = GenerateToyData(*paramPoint, weight);
RooArgSet* saveVarsWithGlobObsSet = (RooArgSet*)allVars->snapshot();
+ detOutAgg.AppendArgSet( fGlobalObservables, "globObs_" );
// evaluate all test statistics
for( unsigned int tsi = 0; tsi < fTestStatistics.size(); tsi++ ) {
*allVars = *saveVarsWithGlobObsSet;
*allVars = *fParametersForTestStat;
-
+
// evaluate test statistic; only depends on null POI
- Double_t value = fTestStatistics[tsi]->Evaluate(*toydata, *(RooArgSet*)fParametersForTestStat);
+ RooArgSet* parForTS = (RooArgSet*)fParametersForTestStat->snapshot();
+ Double_t value = fTestStatistics[tsi]->Evaluate(*toydata, *parForTS);
+ delete parForTS;
allVarsToBeSaved->setRealValue( namesOfTSColumns[tsi], value );
// get detailed output, construct name in dataset, store
- const RooArgSet* detOut = fTestStatistics[tsi]->GetDetailedOutput();
- if( detOut ) {
- TIterator* iter = detOut->createIterator();
- while(RooAbsArg* v = dynamic_cast<RooAbsArg*>( iter->Next() ) ) {
- TString varName( TString::Format( "%s_%s", namesOfTSColumns[tsi].Data(), v->GetName() ) );
- ((RooRealVar&)((*allVarsToBeSaved)[varName])).setVal( ((RooRealVar*)v)->getVal() );
- }
- //cout << "detOut: " << endl;
- //detOut->Print("v");
- //cout << "allVarsToSave: " << endl;
- //allVarsToBeSaved->Print("v");
- }
+ const RooArgSet *ndetout = fTestStatistics[tsi]->GetDetailedOutput();
+ if (ndetout != NULL)
+ detOutAgg.AppendArgSet(ndetout, TString(namesOfTSColumns[tsi]).Append("_"));
if(valueFirst < 0.0) { valueFirst = value; }
}
delete saveVarsWithGlobObsSet;
delete toydata;
- //cout << "weight for this data entry: " << weight << endl;
- outputs->add( *allVarsToBeSaved, weight );
-
-
// check for nan
if(valueFirst != valueFirst) {
oocoutW((TObject*)NULL, Generation) << "skip: " << valueFirst << ", " << weight << endl;
continue;
}
+ detOutAgg.AppendArgSet(allVarsToBeSaved);
+ detOutAgg.CommitSet(weight);
+
// adaptive sampling checks
if (valueFirst <= fAdaptiveLowLimit || valueFirst >= fAdaptiveHighLimit) {
if(weight >= 0.) toysInTails += weight;
@@ -439,6 +393,7 @@ RooDataSet* ToyMCSampler::GetSamplingDistributionsSingleWorker(RooArgSet& paramP
}
}
+ delete allVarsToBeSaved;
// clean up
*allVars = *saveAll;
@@ -446,7 +401,7 @@ RooDataSet* ToyMCSampler::GetSamplingDistributionsSingleWorker(RooArgSet& paramP
delete allVars;
delete paramPoint;
- return outputs;
+ return detOutAgg.GetAsDataSet(fSamplingDistName, fSamplingDistName);
}
void ToyMCSampler::GenerateGlobalObservables(RooAbsPdf& pdf) const {
@@ -546,7 +501,7 @@ RooAbsData* ToyMCSampler::GenerateToyData(RooArgSet& paramPoint, double& weight,
// get nuisance parameter point and weight
fNuisanceParametersSampler->NextPoint(allVarsMinusParamPoint, weight);
}else{
- weight = -1.0;
+ weight = 1.0;
}
RooAbsData *data = Generate(pdf, observables);
@@ -662,7 +617,7 @@ SamplingDistribution* ToyMCSampler::AppendSamplingDistribution(
void ToyMCSampler::ClearCache() {
// clear the cache obtained from the pdf used for speeding the toy and global observables generation
- // needs to be called every time the model pdf (fPdf) changes
+ // needs to be called every time the model pdf (fPdf) changes
if (_gs1) delete _gs1;
diff --git a/roofit/roostats/src/ToyMCStudy.cxx b/roofit/roostats/src/ToyMCStudy.cxx
index d5cec61fda0..b4a48676e19 100644
--- a/roofit/roostats/src/ToyMCStudy.cxx
+++ b/roofit/roostats/src/ToyMCStudy.cxx
@@ -84,7 +84,7 @@ RooDataSet* ToyMCStudy::merge() {
continue;
}
- if( !samplingOutput ) samplingOutput = new RooDataSet(*oneWorker->GetSamplingDistributions(), "proofOutput");
+ if( !samplingOutput ) samplingOutput = new RooDataSet(*oneWorker->GetSamplingDistributions());
else samplingOutput->append( *oneWorker->GetSamplingDistributions() );
//delete oneWorker;
--
GitLab