From 05eb74d5891eab30b15356fc15cb85ca53c249be Mon Sep 17 00:00:00 2001
From: Fons Rademakers <Fons.Rademakers@cern.ch>
Date: Tue, 29 Jul 2008 11:58:09 +0000
Subject: [PATCH] first import of the RooStats code. Still many coding
violations to be fixed.
git-svn-id: http://root.cern.ch/svn/root/trunk@24970 27541ba8-7e3a-0410-8455-c3a389f83636
---
Makefile | 4 +-
config/Makefile.depend | 59 +-
roofit/roostats/Module.mk | 70 ++
roofit/roostats/inc/LinkDef.h | 41 +
roofit/roostats/inc/NumberCountingUtils.h | 61 ++
roofit/roostats/inc/SPlot.h | 55 ++
roofit/roostats/src/NumberCountingUtils.cxx | 335 ++++++++
roofit/roostats/src/SPlot.cxx | 841 ++++++++++++++++++++
8 files changed, 1445 insertions(+), 21 deletions(-)
create mode 100644 roofit/roostats/Module.mk
create mode 100644 roofit/roostats/inc/LinkDef.h
create mode 100644 roofit/roostats/inc/NumberCountingUtils.h
create mode 100644 roofit/roostats/inc/SPlot.h
create mode 100644 roofit/roostats/src/NumberCountingUtils.cxx
create mode 100644 roofit/roostats/src/SPlot.cxx
diff --git a/Makefile b/Makefile
index 12de714718a..1cac5d3ee3b 100644
--- a/Makefile
+++ b/Makefile
@@ -199,7 +199,7 @@ ifeq ($(BUILDCINTEX),yes)
MODULES += cint/cintex
endif
ifeq ($(BUILDROOFIT),yes)
-MODULES += roofit/roofitcore roofit/roofit
+MODULES += roofit/roofitcore roofit/roofit roofit/roostats
endif
ifeq ($(BUILDGDML),yes)
MODULES += geom/gdml
@@ -264,7 +264,7 @@ MODULES += core/unix core/winnt graf2d/x11 graf2d/x11ttf \
graf2d/qt gui/qtroot gui/qtgsi net/xrootd net/netx net/alien \
proof/proofd proof/proofx proof/clarens proof/peac \
sql/oracle io/xmlparser math/mathmore cint/reflex cint/cintex \
- cint/cint7 roofit/roofitcore roofit/roofit \
+ cint/cint7 roofit/roofitcore roofit/roofit roofit/roostats \
math/minuit2 net/monalisa math/fftw sql/odbc math/unuran \
geom/gdml graf3d/eve montecarlo/g4root net/glite misc/memstat
MODULES := $(sort $(MODULES)) # removes duplicates
diff --git a/config/Makefile.depend b/config/Makefile.depend
index 35b59951da7..0503f973eed 100644
--- a/config/Makefile.depend
+++ b/config/Makefile.depend
@@ -72,7 +72,10 @@ ALIENLIBDEPM = $(XMLLIB) $(NETXLIB) $(TREELIB) $(PROOFPLAYERLIB) \
$(NETLIB) $(IOLIB)
ROOFITCORELIBDEPM = $(HISTLIB) $(GRAFLIB) $(MATRIXLIB) $(TREELIB) \
$(MINUITLIB) $(IOLIB) $(MATHCORELIB)
-ROOFITLIBDEPM = $(ROOFITCORELIB) $(TREELIB) $(IOLIB) $(MATRIXLIB) $(MATHCORELIB)
+ROOFITLIBDEPM = $(ROOFITCORELIB) $(TREELIB) $(IOLIB) $(MATRIXLIB) \
+ $(MATHCORELIB)
+ROOSTATSLIBDEPM = $(ROOFITLIB) $(ROOFITCORELIB) $(TREELIB) $(IOLIB) \
+ $(MATRIXLIB) $(HISTLIB) $(MATHCORELIB)
CINT7LIBDEPM = $(REFLEXLIB)
CINTEXLIBDEPM = $(REFLEXLIB)
REFLEXDICTLIBDEPM = $(REFLEXLIB)
@@ -162,6 +165,7 @@ FOAMLIBDEP = $(FOAMLIBDEPM)
ALIENLIBDEP = $(ALIENLIBDEPM)
ROOFITCORELIBDEP = $(ROOFITCORELIBDEPM)
ROOFITLIBDEP = $(ROOFITLIBDEPM)
+ROOSTATSLIBDEP = $(ROOSTATSLIBDEPM)
CINT7LIBDEP = $(CINT7LIBDEPM)
CINTEXLIBDEP = $(CINTEXLIBDEPM)
REFEXDICTLIBDEP = $(REFLEXDICTLIBDEPM)
@@ -231,11 +235,13 @@ PROOFPLAYERLIBEXTRA = lib/libProof.lib lib/libHist.lib lib/libTree.lib \
lib/libRIO.lib lib/libNet.lib lib/libThread.lib \
lib/libMathCore.lib
PROOFDRAWLIBEXTRA = lib/libTreePlayer.lib lib/libGraf3d.lib \
- lib/libGraf.lib lib/libGpad.lib lib/libProofPlayer.lib \
- lib/libHist.lib lib/libTree.lib lib/libProof.lib
+ lib/libGraf.lib lib/libGpad.lib \
+ lib/libProofPlayer.lib lib/libHist.lib \
+ lib/libTree.lib lib/libProof.lib
PROOFXLIBEXTRA = lib/libNet.lib lib/libProof.lib lib/libThread.lib
-SESSIONVIEWERLIBEXTRA = lib/libProof.lib lib/libGui.lib lib/libHist.lib lib/libGpad.lib \
- lib/libGraf.lib lib/libTree.lib lib/libMathCore.lib
+SESSIONVIEWERLIBEXTRA = lib/libProof.lib lib/libGui.lib lib/libHist.lib \
+ lib/libGpad.lib lib/libGraf.lib lib/libTree.lib \
+ lib/libMathCore.lib
PEACLIBEXTRA = lib/libProof.lib lib/libClarens.lib
PEACGUILIBEXTRA = lib/libGui.lib
EGLIBEXTRA = lib/libGraf3d.lib lib/libGraf.lib lib/libGpad.lib \
@@ -250,15 +256,16 @@ TABLELIBEXTRA = lib/libTree.lib lib/libGpad.lib lib/libRIO.lib \
lib/libGraf3d.lib lib/libGraf.lib lib/libHist.lib \
lib/libMathCore.lib
MLPLIBEXTRA = lib/libHist.lib lib/libMatrix.lib lib/libTree.lib \
- lib/libGraf.lib lib/libGpad.lib lib/libTreePlayer.lib \
- lib/libMathCore.lib
+ lib/libGraf.lib lib/libGpad.lib \
+ lib/libTreePlayer.lib lib/libMathCore.lib
SPECTRUMLIBEXTRA = lib/libHist.lib lib/libMatrix.lib
TMVALIBEXTRA = lib/libRIO.lib lib/libHist.lib lib/libMatrix.lib \
lib/libTree.lib lib/libGraf.lib lib/libGpad.lib \
lib/libTreePlayer.lib lib/libMLP.lib \
lib/libMinuit.lib lib/libMathCore.lib
SPLOTLIBEXTRA = lib/libMatrix.lib lib/libHist.lib lib/libTree.lib \
- lib/libTreePlayer.lib lib/libGraf3d.lib lib/libGraf.lib
+ lib/libTreePlayer.lib lib/libGraf3d.lib \
+ lib/libGraf.lib
QUADPLIBEXTRA = lib/libMatrix.lib
GLLIBEXTRA = lib/libGraf3d.lib lib/libGui.lib lib/libGraf.lib \
lib/libHist.lib lib/libGed.lib lib/libMathCore.lib
@@ -284,7 +291,11 @@ ALIENLIBEXTRA = lib/libXMLIO.lib lib/libNetx.lib lib/libTree.lib \
ROOFITCORELIBEXTRA = lib/libHist.lib lib/libGraf.lib lib/libMatrix.lib \
lib/libTree.lib lib/libMinuit.lib lib/libRIO.lib \
lib/libMathCore.lib
-ROOFITLIBEXTRA = lib/libRooFitCore.lib lib/libTree.lib lib/libRIO.lib lib/libMatrix.lib lib/libMathCore.lib
+ROOFITLIBEXTRA = lib/libRooFitCore.lib lib/libTree.lib lib/libRIO.lib \
+ lib/libMatrix.lib lib/libMathCore.lib
+ROOSTATSLIBEXTRA = lib/libRooFit.lib lib/libRooFitCore.lib \
+ lib/libTree.lib lib/libRIO.lib lib/libHist.lib \
+ lib/libMatrix.lib lib/libMathCore.lib
CINT7LIBEXTRA = lib/libReflex.lib
CINTEXLIBEXTRA = lib/libReflex.lib
REFLEXDICTLIBEXTRA = lib/libReflex.lib
@@ -334,15 +345,18 @@ MINUITLIBEXTRA = -Llib -lGraf -lHist -lMatrix -lMathCore
MINUIT2LIBEXTRA = -Llib -lGraf -lHist -lMatrix -lMathCore
FUMILILIBEXTRA = -Llib -lGraf -lHist -lMathCore
TREELIBEXTRA = -Llib -lNet -lRIO -lThread
-TREEPLAYERLIBEXTRA = -Llib -lTree -lGraf3d -lGraf -lHist -lGpad -lRIO -lMathCore
+TREEPLAYERLIBEXTRA = -Llib -lTree -lGraf3d -lGraf -lHist -lGpad -lRIO \
+ -lMathCore
TREEVIEWERLIBEXTRA = -Llib -lTree -lGpad -lGraf -lHist -lGui -lTreePlayer \
-lGed -lRIO -lMathCore
PROOFLIBEXTRA = -Llib -lNet -lTree -lThread -lRIO -lMathCore
-PROOFPLAYERLIBEXTRA = -Llib -lProof -lHist -lTree -lRIO -lNet -lThread -lMathCore
-PROOFDRAWLIBEXTRA = -Llib -lTreePlayer -lGraf3d -lGraf -lGpad -lProofPlayer \
- -lHist -lTree -lProof
+PROOFPLAYERLIBEXTRA = -Llib -lProof -lHist -lTree -lRIO -lNet -lThread \
+ -lMathCore
+PROOFDRAWLIBEXTRA = -Llib -lTreePlayer -lGraf3d -lGraf -lGpad \
+ -lProofPlayer -lHist -lTree -lProof
PROOFXLIBEXTRA = -Llib -lNet -lProof -lThread
-SESSIONVIEWERLIBEXTRA = -Llib -lProof -lGui -lHist -lGpad -lGraf -lTree -lMathCore
+SESSIONVIEWERLIBEXTRA = -Llib -lProof -lGui -lHist -lGpad -lGraf -lTree \
+ -lMathCore
PEACLIBEXTRA = -Llib -lProof -lClarens
PEACGUILIBEXTRA = -Llib -lGui
X3DLIBEXTRA = -Llib -lGraf3d -lGui
@@ -354,12 +368,15 @@ PYTHIA6LIBEXTRA = -Llib -lEG -lGraf -lVMC -lPhysics $(FPYTHIA6LIBDIR) \
PYTHIA8LIBEXTRA = -Llib -lEG -lGraf -lVMC -lPhysics $(FPYTHIA8LIBDIR) \
$(FPYTHIA8LIB)
X11TTFLIBEXTRA = -Llib -lGX11 -lGraf
-TABLELIBEXTRA = -Llib -lTree -lGpad -lGraf3d -lGraf -lHist -lRIO -lMathCore
-GLLIBEXTRA = -Llib -lGpad -lGraf3d -lGui -lGraf -lHist -lGed -lMathCore
+TABLELIBEXTRA = -Llib -lTree -lGpad -lGraf3d -lGraf -lHist -lRIO \
+ -lMathCore
+GLLIBEXTRA = -Llib -lGpad -lGraf3d -lGui -lGraf -lHist -lGed \
+ -lMathCore
HBOOKLIBEXTRA = -Llib -lHist -lMatrix -lTree -lGraf -lTreePlayer \
-lRIO -lminicern
GEOMLIBEXTRA = -Llib -lRIO -lMathCore
-GEOMPAINTERLIBEXTRA = -Llib -lGeom -lTree -lGraf3d -lHist -lGpad -lRIO -lMathCore
+GEOMPAINTERLIBEXTRA = -Llib -lGeom -lTree -lGraf3d -lHist -lGpad -lRIO \
+ -lMathCore
GEOMBUILDERLIBEXTRA = -Llib -lGeom -lGraf3d -lGpad -lGraf -lGui -lGed
ASIMAGELIBEXTRA = -Llib -lGraf -lMathCore
ASIMAGEGUILIBEXTRA = -Llib -lGraf -lHist -lGui -lASImage -lRIO
@@ -369,7 +386,8 @@ MLPLIBEXTRA = -Llib -lHist -lMatrix -lTree -lGraf -lGpad \
SPECTRUMLIBEXTRA = -Llib -lHist -lMatrix
TMVALIBEXTRA = -Llib -lRIO -lHist -lMatrix -lTree -lGraf -lGpad \
-lTreePlayer -lMLP -lMinuit -lMathCore
-SPLOTLIBEXTRA = -Llib -lMatrix -lHist -lTree -lTreePlayer -lGraf3d -lGraf
+SPLOTLIBEXTRA = -Llib -lMatrix -lHist -lTree -lTreePlayer -lGraf3d \
+ -lGraf
QTROOTLIBEXTRA = -Llib -lGui -lGQt
GQTLIBEXTRA = -Llib -lGui -lGpad -lGraf -lRint
QUADPLIBEXTRA = -Llib -lMatrix
@@ -379,8 +397,11 @@ GUIHTMLLIBEXTRA = -Llib -lGui -lGraf -lNet
XMLLIBEXTRA = -Llib -lRIO
FOAMLIBEXTRA = -Llib -lHist
ALIENLIBEXTRA = -Llib -lXMLIO -lNetx -lTree -lProofPlayer -lNet -lRIO
-ROOFITCORELIBEXTRA = -Llib -lHist -lGraf -lMatrix -lTree -lMinuit -lRIO -lMathCore
+ROOFITCORELIBEXTRA = -Llib -lHist -lGraf -lMatrix -lTree -lMinuit -lRIO \
+ -lMathCore
ROOFITLIBEXTRA = -Llib -lRooFitCore -lTree -lRIO -lMatrix -lMathCore
+ROOSTATSLIBEXTRA = -Llib -lRooFit -lRooFitCore -lTree -lRIO -lMatrix \
+ -lHist -lMathCore
CINT7LIBEXTRA = -Llib -lReflex
CINTEXLIBEXTRA = -Llib -lReflex
REFLEXDICTLIBEXTRA = -Llib -lReflex
diff --git a/roofit/roostats/Module.mk b/roofit/roostats/Module.mk
new file mode 100644
index 00000000000..83ddc786c1a
--- /dev/null
+++ b/roofit/roostats/Module.mk
@@ -0,0 +1,70 @@
+# Module.mk for roostats module
+# Copyright (c) 2008 Rene Brun and Fons Rademakers
+#
+# Author: Kyle Cranmer
+
+MODNAME := roostats
+MODDIR := roofit/$(MODNAME)
+MODDIRS := $(MODDIR)/src
+MODDIRI := $(MODDIR)/inc
+
+ROOSTATSDIR := $(MODDIR)
+ROOSTATSDIRS := $(ROOSTATSDIR)/src
+ROOSTATSDIRI := $(ROOSTATSDIR)/inc
+
+##### libRooStats #####
+ROOSTATSL := $(MODDIRI)/LinkDef.h
+ROOSTATSDS := $(MODDIRS)/G__RooStats.cxx
+ROOSTATSDO := $(ROOSTATSDS:.cxx=.o)
+ROOSTATSDH := $(ROOSTATSDS:.cxx=.h)
+
+ROOSTATSH := $(filter-out $(MODDIRI)/LinkDef%,$(wildcard $(MODDIRI)/*.h))
+ROOSTATSS := $(filter-out $(MODDIRS)/G__%,$(wildcard $(MODDIRS)/*.cxx))
+ROOSTATSO := $(ROOSTATSS:.cxx=.o)
+
+ROOSTATSDEP := $(ROOSTATSO:.o=.d) $(ROOSTATSDO:.o=.d)
+
+ROOSTATSLIB := $(LPATH)/libRooStats.$(SOEXT)
+ROOSTATSMAP := $(ROOSTATSLIB:.$(SOEXT)=.rootmap)
+
+# used in the main Makefile
+ALLHDRS += $(patsubst $(MODDIRI)/%.h,include/%.h,$(ROOSTATSH))
+ALLLIBS += $(ROOSTATSLIB)
+ALLMAPS += $(ROOSTATSMAP)
+
+# include all dependency files
+INCLUDEFILES += $(ROOSTATSDEP)
+
+##### local rules #####
+.PHONY: all-$(MODNAME) clean-$(MODNAME) distclean-$(MODNAME)
+
+include/%.h: $(ROOSTATSDIRI)/%.h
+ cp $< $@
+
+$(ROOSTATSLIB): $(ROOSTATSO) $(ROOSTATSDO) $(ORDER_) $(MAINLIBS) \
+ $(ROOSTATSLIBDEP)
+ @$(MAKELIB) $(PLATFORM) $(LD) "$(LDFLAGS)" \
+ "$(SOFLAGS)" libRooStats.$(SOEXT) $@ \
+ "$(ROOSTATSO) $(ROOSTATSDO)" \
+ "$(ROOSTATSLIBEXTRA)"
+
+$(ROOSTATSDS): $(ROOSTATSH) $(ROOSTATSL) $(ROOTCINTTMPDEP)
+ @echo "Generating dictionary $@..."
+ $(ROOTCINTTMP) -f $@ -c $(ROOSTATSH) $(ROOSTATSL)
+
+$(ROOSTATSMAP): $(RLIBMAP) $(MAKEFILEDEP) $(ROOSTATSL)
+ $(RLIBMAP) -o $(ROOSTATSMAP) -l $(ROOSTATSLIB) \
+ -d $(ROOSTATSLIBDEPM) -c $(ROOSTATSL)
+
+all-$(MODNAME): $(ROOSTATSLIB) $(ROOSTATSMAP)
+
+clean-$(MODNAME):
+ @rm -f $(ROOSTATSO) $(ROOSTATSDO)
+
+clean:: clean-$(MODNAME)
+
+distclean-$(MODNAME): clean-$(MODNAME)
+ @rm -rf $(ROOSTATSDEP) $(ROOSTATSLIB) $(ROOSTATSMAP) \
+ $(ROOSTATSDS) $(ROOSTATSDH)
+
+distclean:: distclean-$(MODNAME)
diff --git a/roofit/roostats/inc/LinkDef.h b/roofit/roostats/inc/LinkDef.h
new file mode 100644
index 00000000000..274b4c586f0
--- /dev/null
+++ b/roofit/roostats/inc/LinkDef.h
@@ -0,0 +1,41 @@
+/* @(#)root/roostats:$Id$ */
+
+/*************************************************************************
+ * 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. *
+ *************************************************************************/
+
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ namespace NumberCountingUtils;
+#pragma link C++ namespace Statistics;
+
+// for auto-loading namespaces
+#ifdef USE_FOR_AUTLOADING
+#pragma link C++ class NumberCountingUtils;
+#pragma link C++ class Statistics;
+#endif
+
+#pragma link C++ class SPlot+;
+
+#pragma link C++ function NumberCountingUtils::BinomialExpZ(Double_t,Double_t,Double_t);
+#pragma link C++ function NumberCountingUtils::BinomialWithTauExpZ(Double_t,Double_t,Double_t);
+#pragma link C++ function NumberCountingUtils::BinomialObsZ(Double_t,Double_t,Double_t);
+#pragma link C++ function NumberCountingUtils::BinomialWithTauObsZ(Double_t,Double_t,Double_t);
+#pragma link C++ function NumberCountingUtils::BinomialExpP(Double_t,Double_t,Double_t);
+#pragma link C++ function NumberCountingUtils::BinomialWithTauExpP(Double_t,Double_t,Double_t);
+#pragma link C++ function NumberCountingUtils::BinomialObsP(Double_t,Double_t,Double_t);
+#pragma link C++ function NumberCountingUtils::BinomialWithTauObsP(Double_t,Double_t,Double_t);
+
+#pragma link C++ function NumberCountingUtils::ProfileCombinationExpZ(Double_t*,Double_t*,Double_t*,Int_t);
+
+#pragma link C++ function Statistics::PValueToSignificance(Double_t);
+
+#endif
diff --git a/roofit/roostats/inc/NumberCountingUtils.h b/roofit/roostats/inc/NumberCountingUtils.h
new file mode 100644
index 00000000000..ae587de901c
--- /dev/null
+++ b/roofit/roostats/inc/NumberCountingUtils.h
@@ -0,0 +1,61 @@
+// @(#)root/roostats:$Id$
+// Author: Kyle Cranmer 28/07/2008
+
+/*************************************************************************
+ * 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 ROOT_NumberCountingUtils
+#define ROOT_NumberCountingUtils
+
+#ifndef ROOT_TMath
+#include "TMath.h"
+#endif
+
+namespace Statistics {
+ inline Double_t PValueToSignificance(Double_t pvalue){
+// return sqrt(2.)*TMath::ErfInverse(1 - 2.*pvalue);
+ return TMath::Abs(TMath::NormQuantile(pvalue) );
+ }
+}
+
+/////////////////////////////////////////
+// NumberCountingUtils
+//
+// Encapsulates common number counting utilities
+/////////////////////////////////////////
+
+namespace NumberCountingUtils {
+
+ ///////////////////////////////////
+ // Standalone Functions.
+ // Naming conventions:
+ // Exp = Expected
+ // Obs = Observed
+ // P = p-value
+ // Z = Z-value or significance in Sigma (one-sided convention)
+ //////////////////////////////////
+
+
+ Double_t BinomialExpZ(Double_t, Double_t, Double_t);
+ Double_t BinomialWithTauExpZ(Double_t, Double_t, Double_t);
+ Double_t BinomialObsZ(Double_t, Double_t, Double_t);
+ Double_t BinomialWithTauObsZ(Double_t, Double_t, Double_t);
+
+ Double_t BinomialExpP(Double_t, Double_t, Double_t);
+ Double_t BinomialWithTauExpP(Double_t, Double_t, Double_t);
+ Double_t BinomialObsP(Double_t, Double_t, Double_t);
+ Double_t BinomialWithTauObsP(Double_t, Double_t, Double_t);
+
+ ///////////////////////////////////
+ // RooFit based Functions
+ //////////////////////////////////
+ Double_t ProfileCombinationExpZ(Double_t*, Double_t*, Double_t*, Int_t );
+
+}
+
+#endif
diff --git a/roofit/roostats/inc/SPlot.h b/roofit/roostats/inc/SPlot.h
new file mode 100644
index 00000000000..9c2278610ba
--- /dev/null
+++ b/roofit/roostats/inc/SPlot.h
@@ -0,0 +1,55 @@
+// @(#)root/roostats:$Id$
+// Author: Kyle Cranmer 21/07/2008
+
+/*************************************************************************
+ * 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 ROOT_SPlot
+#define ROOT_SPlot
+
+class RooAbsReal;
+class RooAbsPdf;
+class RooFitResult;
+class RooRealVar;
+class RooSimultaneous;
+
+#ifndef ROOT_TH1
+#include "TH1.h"
+#endif
+#include "RooFitResult.h"
+#include "RooRealVar.h"
+#include "RooHist.h"
+#include "RooPlot.h"
+#include "RooDataSet.h"
+
+class SPlot : public TH1F {
+public:
+ SPlot();
+ SPlot(const SPlot &other);
+ SPlot(const char* name, const char *title, Int_t nbins, Double_t xmin, Double_t xmax);
+
+ static RooDataSet*
+ AddSWeightToData(const RooSimultaneous* pdf, const RooArgList &yieldsTmp,
+ RooDataSet &data, const RooArgSet &projDeps=RooArgSet()) ;
+
+
+ void FillSPlot(const RooDataSet &data, TString varname, TString weightname);
+
+ void FillSPlot(const RooAbsReal &x, RooAbsReal &nstar, RooDataSet data, const RooFitResult &fitRes, const RooArgList &pdfList, const RooArgList &yields, Bool_t doErrors, const RooArgSet &projDeps=RooArgSet() );
+ void FillSPlot(const RooAbsReal &x, RooAbsReal &nstar, RooDataSet data, const RooFitResult &fitRes, const RooArgList &pdfList, const RooArgList &yields, RooAbsPdf &totalPdf, Bool_t doErrors, const RooArgSet &projDeps=RooArgSet() );
+
+ void FillSPlot(const RooAbsReal &x, RooAbsReal &nstar, RooDataSet data, const RooFitResult &fitRes, RooAbsPdf &totalPdf, RooArgList &yields, Bool_t doErrors, const RooArgSet &projDeps=RooArgSet() );
+
+ Double_t GetComponentValue(RooAbsPdf &pdf, RooArgList &yieldsTmp, Int_t igood, RooArgSet &normSet);
+
+protected:
+
+ ClassDef(SPlot,1) // Make sPlot
+};
+
+#endif
diff --git a/roofit/roostats/src/NumberCountingUtils.cxx b/roofit/roostats/src/NumberCountingUtils.cxx
new file mode 100644
index 00000000000..d9dd9cf8266
--- /dev/null
+++ b/roofit/roostats/src/NumberCountingUtils.cxx
@@ -0,0 +1,335 @@
+// @(#)root/roostats:$Id$
+// Author: Kyle Cranmer 28/07/2008
+
+/*************************************************************************
+ * 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. *
+ *************************************************************************/
+
+/////////////////////////////////////////
+// NumberCountingUtils
+//
+// Encapsulates common number counting utilities
+/////////////////////////////////////////
+ ///////////////////////////////////
+ // Standalone Functions.
+ // Naming conventions:
+ // Exp = Expected
+ // Obs = Observed
+ // P = p-value
+ // Z = Z-value or significance in Sigma (one-sided convention)
+ //////////////////////////////////
+
+#include "NumberCountingUtils.h"
+
+
+// Without this macro the THtml doc for TMath can not be generated
+#if !defined(R__ALPHA) && !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD)
+NamespaceImp(NumberCountingUtils)
+#endif
+
+Double_t NumberCountingUtils::BinomialExpP(Double_t signalExp, Double_t backgroundExp, Double_t relativeBkgUncert){
+ // Expected P-value for s=0 in a ratio of Poisson means.
+ // Here the background and its uncertainty are provided directly and
+ // assumed to be from the double Poisson counting setup described in the
+ // BinomialWithTau functions.
+ // Normally one would know tau directly, but here it is determiend from
+ // the background uncertainty. This is not strictly correct, but a useful
+ // approximation.
+ //
+ //
+ // This is based on code and comments from Bob Cousins
+ // based on the following papers:
+ //
+ // Statistical Challenges for Searches for New Physics at the LHC
+ // Authors: Kyle Cranmer
+ // http://arxiv.org/abs/physics/0511028
+ //
+ // Measures of Significance in HEP and Astrophysics
+ // Authors: J. T. Linnemann
+ // http://arxiv.org/abs/physics/0312059
+ //
+ // In short, this is the exact frequentist solution to the problem of
+ // a main measurement x distributed as a Poisson around s+b and a sideband or
+ // auxiliary measurement y distributed as a Poisson around \taub. Eg.
+ // L(x,y|s,b,\tau) = Pois(x|s+b)Pois(y|\tau b)
+
+
+ //SIDE BAND EXAMPLE
+ //See Eqn. (19) of Cranmer and pp. 36-37 of Linnemann.
+ //150 total events in signalExp region, 100 in sideband of equal size
+ Double_t mainInf = signalExp+backgroundExp; //Given
+ Double_t tau = 1./backgroundExp/(relativeBkgUncert*relativeBkgUncert);
+ Double_t auxiliaryInf = backgroundExp*tau; //Given
+
+ Double_t P_Bi = TMath::BetaIncomplete(1./(1.+tau),mainInf,auxiliaryInf+1);
+ return P_Bi;
+
+/*
+Now, if instead the mean background level b in the signal region is
+specified, along with Gaussian rms sigb, then one can fake a Poisson
+sideband (see Linnemann, p. 35, converted to Cranmer's notation) by
+letting tau = b/(sigb*sigb) and y = b*tau. Thus, for example, if one
+has x=150 and b = 100 +/- 10, one then derives tau and y. Then one
+has the same two lines of ROOT calling BetaIncomplete and ErfInverse.
+Since I chose these numbers to revert to the previous example, we get
+the same answer:
+*/
+/*
+//GAUSSIAN FAKED AS POISSON EXAMPLE
+x = 150. //Given
+b = 100. //Given
+sigb = 10. //Given
+tau = b/(sigb*sigb)
+y = tau*b
+Z_Bi = TMath::BetaIncomplete(1./(1.+tau),x,y+1)
+S = sqrt(2)*TMath::ErfInverse(1 - 2*Z_Bi)
+
+*/
+
+}
+
+
+Double_t NumberCountingUtils::BinomialWithTauExpP(Double_t signalExp, Double_t backgroundExp, Double_t tau){
+ // Expected P-value for s=0 in a ratio of Poisson means.
+ // Based on two expectations, a main measurement that might have signal
+ // and an auxiliarly measurement for the background that is signal free.
+ // The expected background in the auxiliary measurement is a factor
+ // tau larger than in the main measurement.
+
+ Double_t mainInf = signalExp+backgroundExp; //Given
+ Double_t auxiliaryInf = backgroundExp*tau; //Given
+
+ Double_t P_Bi = TMath::BetaIncomplete(1./(1.+tau),mainInf,auxiliaryInf+1);
+
+ return P_Bi;
+
+}
+
+Double_t NumberCountingUtils::BinomialObsP(Double_t mainObs, Double_t backgroundObs, Double_t relativeBkgUncert){
+ // P-value for s=0 in a ratio of Poisson means.
+ // Here the background and its uncertainty are provided directly and
+ // assumed to be from the double Poisson counting setup.
+ // Normally one would know tau directly, but here it is determiend from
+ // the background uncertainty. This is not strictly correct, but a useful
+ // approximation.
+
+ Double_t tau = 1./backgroundObs/(relativeBkgUncert*relativeBkgUncert);
+ Double_t auxiliaryInf = backgroundObs*tau; //Given
+
+
+ //SIDE BAND EXAMPLE
+ //See Eqn. (19) of Cranmer and pp. 36-37 of Linnemann.
+ Double_t P_Bi = TMath::BetaIncomplete(1./(1.+tau),mainObs,auxiliaryInf+1);
+
+ return P_Bi;
+
+}
+
+
+Double_t NumberCountingUtils::BinomialWithTauObsP(Double_t mainObs, Double_t auxiliaryObs, Double_t tau){
+ // P-value for s=0 in a ratio of Poisson means.
+ // Based on two observations, a main measurement that might have signal
+ // and an auxiliarly measurement for the background that is signal free.
+ // The expected background in the auxiliary measurement is a factor
+ // tau larger than in the main measurement.
+
+ //SIDE BAND EXAMPLE
+ //See Eqn. (19) of Cranmer and pp. 36-37 of Linnemann.
+ Double_t P_Bi = TMath::BetaIncomplete(1./(1.+tau),mainObs,auxiliaryObs+1);
+
+ return P_Bi;
+
+}
+
+Double_t NumberCountingUtils::BinomialExpZ(Double_t signalExp, Double_t backgroundExp, Double_t relativeBkgUncert) {
+ // See BinomialExpP
+ return Statistics::PValueToSignificance( BinomialExpP(signalExp,backgroundExp,relativeBkgUncert) ) ;
+ }
+
+Double_t NumberCountingUtils::BinomialWithTauExpZ(Double_t signalExp, Double_t backgroundExp, Double_t tau){
+ // See BinomialWithTauExpP
+ return Statistics::PValueToSignificance( BinomialWithTauExpP(signalExp,backgroundExp,tau) ) ;
+}
+
+
+Double_t NumberCountingUtils::BinomialObsZ(Double_t mainObs, Double_t backgroundObs, Double_t relativeBkgUncert){
+ // See BinomialObsZ
+ return Statistics::PValueToSignificance( BinomialObsP(mainObs,backgroundObs,relativeBkgUncert) ) ;
+}
+
+Double_t NumberCountingUtils::BinomialWithTauObsZ(Double_t mainObs, Double_t auxiliaryObs, Double_t tau){
+ // See BinomialWithTauObsZ
+ return Statistics::PValueToSignificance( BinomialWithTauObsZ(mainObs,auxiliaryObs,tau) ) ;
+}
+
+/////////////////////////////////////////////////////////////
+//
+// RooFit based Functions
+//
+/////////////////////////////////////////////////////////////
+#include "RooRealVar.h"
+#include "RooConstVar.h"
+#include "RooAddition.h"
+#include "RooProduct.h"
+#include "RooDataSet.h"
+#include "RooProdPdf.h"
+#include "RooPlot.h"
+#include "TMath.h"
+#include "TCanvas.h"
+#include "RooFitResult.h"
+#include "RooNLLVar.h"
+#include "RooPoisson.h"
+#include "RooGaussian.h"
+#include "RooMinuit.h"
+#include "RooGlobalFunc.h"
+#include "RooCmdArg.h"
+#include "TH2F.h"
+#include "TTree.h"
+#include <sstream>
+
+
+Double_t NumberCountingUtils::ProfileCombinationExpZ(Double_t* sig,
+ Double_t* back,
+ Double_t* back_syst,
+ Int_t nbins){
+
+ // A number counting combination for N channels with uncorrelated background
+ // uncertainty.
+ // Background uncertainty taken into account via the profile likelihood ratio.
+ // Arguements are an array of expected signal, expected background, and relative
+ // background uncertainty (eg. 0.1 for 10% uncertainty), and the number of channels.
+
+ using namespace RooFit;
+ using std::vector;
+
+ vector<RooRealVar*> backVec, tauVec, xVec, yVec;
+ vector<RooProduct*> sigVec;
+ vector<RooFormulaVar*> splusbVec;
+ vector<RooPoisson*> sigRegions, sidebands;
+ TList likelihoodFactors;
+ TList observablesCollection;
+
+ TTree* tree = new TTree();
+ Double_t* xForTree = new Double_t[nbins];
+ Double_t* yForTree = new Double_t[nbins];
+
+ Double_t MaxSigma = 8; // Needed to set ranges for varaibles.
+
+ RooRealVar* masterSignal =
+ new RooRealVar("masterSignal","masterSignal",1., 0., 3.);
+ for(Int_t i=0; i<nbins; ++i){
+ std::stringstream str;
+ str<<"_"<<i;
+ RooRealVar* expectedSignal =
+ new RooRealVar(("expected_s"+str.str()).c_str(),("expected_s"+str.str()).c_str(),sig[i], 0., 2*sig[i]);
+ expectedSignal->setConstant(kTRUE);
+
+ RooProduct* s =
+ new RooProduct(("s"+str.str()).c_str(),("s"+str.str()).c_str(), RooArgSet(*masterSignal, *expectedSignal));
+
+ RooRealVar* b =
+ new RooRealVar(("b"+str.str()).c_str(),("b"+str.str()).c_str(),back[i], 0., 1.2*back[i]+MaxSigma*(sqrt(back[i])+back[i]*back_syst[i]));
+ b->Print();
+ Double_t _tau = 1./back[i]/back_syst[i]/back_syst[i];
+ RooRealVar* tau =
+ new RooRealVar(("tau"+str.str()).c_str(),("tau"+str.str()).c_str(),_tau,0,2*_tau);
+ tau->setConstant(kTRUE);
+
+ RooAddition* splusb =
+ new RooAddition(("splusb"+str.str()).c_str(),("s"+str.str()+"+"+"b"+str.str()).c_str(),
+ RooArgSet(*s,*b));
+ RooProduct* bTau =
+ new RooProduct(("bTau"+str.str()).c_str(),("b*tau"+str.str()).c_str(), RooArgSet(*b, *tau));
+ RooRealVar* x =
+ new RooRealVar(("x"+str.str()).c_str(),("x"+str.str()).c_str(), sig[i]+back[i], 0., 1.2*sig[i]+back[i]+MaxSigma*sqrt(sig[i]+back[i]));
+ RooRealVar* y =
+ new RooRealVar(("y"+str.str()).c_str(),("y"+str.str()).c_str(), back[i]*_tau, 0., 1.2*back[i]*_tau+MaxSigma*sqrt(back[i]*_tau));
+
+
+ RooPoisson* sigRegion =
+ new RooPoisson(("sigRegion"+str.str()).c_str(),("sigRegion"+str.str()).c_str(), *x,*splusb);
+ RooPoisson* sideband =
+ new RooPoisson(("sideband"+str.str()).c_str(),("sideband"+str.str()).c_str(), *y,*bTau);
+
+ sigVec.push_back(s);
+ backVec.push_back(b);
+ tauVec.push_back(tau);
+ xVec.push_back(x);
+ yVec.push_back(y);
+ sigRegions.push_back(sigRegion);
+ sidebands.push_back(sideband);
+
+ likelihoodFactors.Add(sigRegion);
+ likelihoodFactors.Add(sideband);
+ observablesCollection.Add(x);
+ observablesCollection.Add(y);
+
+ // print to see range on variables
+ // x->Print();
+ // y->Print();
+ // b->Print();
+
+
+ xForTree[i] = sig[i]+back[i];
+ yForTree[i] = back[i]*_tau;
+ tree->Branch(("x"+str.str()).c_str(), xForTree+i ,("x"+str.str()+"/D").c_str());
+ tree->Branch(("y"+str.str()).c_str(), yForTree+i ,("y"+str.str()+"/D").c_str());
+ }
+ tree->Fill();
+ // tree->Print();
+ // tree->Scan();
+
+ RooArgSet likelihoodFactorSet(likelihoodFactors);
+ RooProdPdf joint("joint","joint", likelihoodFactorSet );
+ // likelihoodFactorSet.Print();
+
+ // cout << "\n print model" << endl;
+ // joint.Print();
+ // joint.printCompactTree();
+
+ // RooArgSet* observableSet = new RooArgSet(observablesCollection);
+ RooArgList* observableList = new RooArgList(observablesCollection);
+
+ // observableSet->Print();
+ // observableList->Print();
+
+ // cout << "Make hypothetical dataset:" << endl;
+ RooDataSet* toyMC = new RooDataSet("data","data", tree, *observableList); // one experiment
+ toyMC->Scan();
+
+ // cout << "about to do fit \n\n" << endl;
+ RooFitResult* fit = joint.fitTo(*toyMC,Extended(kFALSE),Strategy(0),Hesse(kFALSE),Save(kTRUE),PrintLevel(-1));
+
+ //RooFitResult* fit = joint.fitTo(*toyMC,"sr");
+ // fit->Print();
+
+ // joint.Print("v");
+
+ ////////////////////////////////////////
+ /// Calculate significance
+ //////////////////////////////
+ // cout << "\nFit to signal plus background:" << endl;
+ masterSignal->Print();
+ for(Int_t i=0; i<nbins; ++i) backVec.at(i)->Print();
+ fit->Print();
+ Double_t NLLatMLE= fit->minNll();
+
+
+
+ // cout << "\nFit to background only:" << endl;
+ masterSignal->setVal(0);
+ masterSignal->setConstant();
+ RooFitResult* fit2 = joint.fitTo(*toyMC,Extended(kFALSE),Hesse(kFALSE),Strategy(0), Minos(kFALSE), Save(kTRUE),PrintLevel(-1));
+
+ masterSignal->Print();
+ for(Int_t i=0; i<nbins; ++i) backVec.at(i)->Print();
+ Double_t NLLatCondMLE= fit2->minNll();
+ fit2->Print();
+
+ return sqrt( 2*(NLLatCondMLE-NLLatMLE));
+
+}
diff --git a/roofit/roostats/src/SPlot.cxx b/roofit/roostats/src/SPlot.cxx
new file mode 100644
index 00000000000..edb3a77df25
--- /dev/null
+++ b/roofit/roostats/src/SPlot.cxx
@@ -0,0 +1,841 @@
+// @(#)root/roostats:$Id$
+// Author: Kyle Cranmer 28/07/2008
+
+/*************************************************************************
+ * 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. *
+ *************************************************************************/
+
+/*****************************************************************************
+ * Project: RooStats
+ * Package: RooFit/RooStats
+ *
+ * Authors:
+ * Original code from M. Pivk as part of MLFit package from BaBar.
+ * Modifications:
+ * Giacinto Piacquadio, Maurizio Pierini: modifications for new RooFit version
+ * George H. Lewis, Kyle Cranmer: generalized for weighted events
+ *
+ * Porting to RooStats (with permission) by Kyle Cranmer, July 2008
+ * documentation for the multiple versions of fillSplot are needed.
+ *
+ *****************************************************************************/
+
+
+
+// This class calculates sWeights used to create an sPlot.
+// The code is based on
+// ``SPlot: A statistical tool to unfold data distributions,''
+// Nucl. Instrum. Meth. A 555, 356 (2005)
+// [arXiv:physics/0402083].
+//
+// An SPlot gives us the distribution of some variable, x in our
+// data sample for a given species (eg. signal or background).
+// The result is similar to a likelihood projection plot, but no cuts are made,
+// so every event contributes to the distribution.
+//
+// [Usage]
+// To use this class, you first must perform your fit twice:
+// The first time perform your nominal fit.
+// For the second fit, fix your parameters at the minimum, float only your yields
+// (normalizations), and remove any PDFs correlated with the variable of interest.
+// Be sure to save the RooFitResult.
+
+#include <vector>
+#include <map>
+
+#include "SPlot.h"
+#include "RooAbsPdf.h"
+#include "RooDataSet.h"
+#include "RooRealVar.h"
+#include "RooSimultaneous.h"
+
+#include "TMatrixD.h"
+
+ClassImp(SPlot) ;
+
+
+//____________________________________________________________________
+SPlot::SPlot() :
+ TH1F()
+{
+ // Default constructor
+}
+
+//____________________________________________________________________
+SPlot::SPlot(const SPlot &other) :
+ TH1F(other)
+{
+}
+
+//____________________________________________________________________
+SPlot::SPlot(const char* name, const char* title, Int_t nbins, Double_t xmin, Double_t xmax) :
+ TH1F(name, title, nbins, xmin, xmax)
+{
+ // Constructor
+}
+
+
+//____________________________________________________________________
+RooDataSet* SPlot::AddSWeightToData(const RooSimultaneous* pdf,
+ const RooArgList &yieldsTmp,
+ RooDataSet &data, const RooArgSet &projDeps)
+{
+
+ // Method which adds the sWeights to the dataset.
+ // input is the PDF, a RooArgList of the yields (floating)
+ // the dataset to which the sWeights should be added,
+ // and a RooArgSet of the projDeps (needs better description).
+
+ Int_t nspec = yieldsTmp.getSize();
+ RooArgList yields = *(RooArgList*)yieldsTmp.snapshot(kFALSE);
+
+ // The list of variables to normalize over when calculating PDF values.
+ RooArgSet vars(*data.get());
+ vars.remove(projDeps, kTRUE, kTRUE);
+
+ // Attach data set
+ const_cast<RooSimultaneous*>(pdf)->attachDataSet(data);
+
+ // first calculate the pdf values for all species and all events
+ std::vector<RooRealVar*> yieldvars ;
+ RooArgSet* parameters = pdf->getParameters(&data) ;
+ //parameters->Print("V") ;
+
+ std::vector<Double_t> yieldvalues ;
+ for (Int_t k = 0; k < nspec; ++k) {
+ RooRealVar* thisyield = dynamic_cast<RooRealVar*>(yields.at(k)) ;
+ RooRealVar* yieldinpdf = dynamic_cast<RooRealVar*>(parameters->find(thisyield->GetName() )) ;
+ std::cout << "yield in pdf: " << yieldinpdf << " " << thisyield->getVal() << std::endl ;
+ yieldvars.push_back(yieldinpdf) ;
+ yieldvalues.push_back(thisyield->getVal()) ;
+ }
+
+ Int_t numevents = data.numEntries() ;
+
+ std::vector<std::vector<Double_t> > pdfvalues(numevents,std::vector<Double_t>(nspec,0)) ;
+
+ // set all yield to zero
+ for(Int_t m=0; m<nspec; ++m) yieldvars[m]->setVal(0) ;
+
+ for (Int_t ievt = 0; ievt <numevents; ievt++) {
+ if (ievt % 100 == 0) {
+ std::cout << ".";
+ std::cout.flush();
+ }
+ RooArgSet row(*data.get(ievt));
+ for(Int_t k = 0; k < nspec; ++k) {
+ // set this yield to 1
+ yieldvars[k]->setVal( 1 ) ;
+ // evaluate the pdf
+ Double_t f_k = pdf->getVal(&vars) ;
+ pdfvalues[ievt][k] = f_k ;
+ if( !(f_k>1 || f_k<1) ) std::cout << "Strange pdf value: " << ievt << " " << k << " " << f_k << std::endl ;
+ yieldvars[k]->setVal( 0 ) ;
+ }
+ }
+
+ // check that the likelihood normalization is fine
+ std::vector<Double_t> norm(nspec,0) ;
+ for (Int_t ievt = 0; ievt <numevents ; ievt++) {
+ Double_t dnorm(0) ;
+ for(Int_t k=0; k<nspec; ++k) dnorm += yieldvalues[k] * pdfvalues[ievt][k] ;
+ for(Int_t j=0; j<nspec; ++j) norm[j] += pdfvalues[ievt][j]/dnorm ;
+ }
+
+ std::cout << "likelihood norms: " ;
+ for(Int_t k=0; k<nspec; ++k) std::cout << norm[k] << " " ;
+ std::cout << std::endl ;
+
+ // Make a TMatrixD to hold the covariance matrix.
+ TMatrixD covInv(nspec, nspec);
+ for (Int_t i = 0; i < nspec; i++) for (Int_t j = 0; j < nspec; j++) covInv(i,j) = 0;
+
+ std::cout << "Calculating covariance matrix";
+ for (Int_t ievt = 0; ievt < numevents; ++ievt) {
+
+ // Calculate contribution to the inverse of the covariance
+ // matrix. See BAD 509 V2 eqn. 15
+
+ // Sum for the denominator
+ Double_t dsum(0);
+ for(Int_t k = 0; k < nspec; ++k)
+ dsum += pdfvalues[ievt][k] * yieldvalues[k] ;
+
+ for(Int_t n=0; n<nspec; ++n)
+ for(Int_t j=0; j<nspec; ++j)
+ covInv(n,j) += pdfvalues[ievt][n]*pdfvalues[ievt][j]/(dsum*dsum);
+ }
+ // Covariance inverse should now be computed!
+ covInv.Print();
+
+ // Invert to get the covariance matrix
+ if (covInv.Determinant() <=0) {
+ std::cout << "SPlot Error: covariance matrix is singular; I can't invert it!" << std::endl;
+ covInv.Print();
+ return 0 ;
+ }
+
+ TMatrixD covMatrix(TMatrixD::kInverted,covInv);
+ covMatrix.Print() ;
+
+ //check cov normalization
+ for(Int_t k=0; k<nspec; ++k) {
+ Double_t covnorm(0) ;
+ for(Int_t m=0; m<nspec; ++m) covnorm += covInv[k][m]*yieldvalues[m] ;
+ Double_t sumrow(0) ;
+ for(Int_t m = 0; m < nspec; ++m) sumrow += covMatrix[k][m] ;
+ std::cout << yieldvalues[k] << " " << sumrow << " " << covnorm << std::endl ;
+ }
+
+ // calculate for each event the sWeight (BAD 509 V2 eq. 21)
+ std::cout << "Calculating sWeight";
+ std::vector<RooRealVar*> sweightvec ;
+ std::vector<RooRealVar*> pdfvec ;
+ RooArgSet sweightset ;
+
+ char wname[256] ;
+ for(Int_t k=0; k<nspec; ++k) {
+ sprintf(wname,"%s_sw", yieldvars[k]->GetName()) ;
+ RooRealVar* var = new RooRealVar(wname,wname,0) ;
+ sweightvec.push_back( var) ;
+ sweightset.add(*var) ;
+ sprintf(wname,"L_%s", yieldvars[k]->GetName()) ;
+ var = new RooRealVar(wname,wname,0) ;
+ pdfvec.push_back( var) ;
+ sweightset.add(*var) ;
+ }
+ sweightset.add(*data.get()) ;
+ RooDataSet* sWeightData = new RooDataSet("dataset", "dataset with sWeights", sweightset);
+
+ for(Int_t ievt = 0; ievt < numevents; ++ievt) {
+
+ data.get(ievt) ;
+
+ // sum for denominator
+ Double_t dsum(0);
+ for(Int_t k = 0; k < nspec; ++k) dsum += pdfvalues[ievt][k] * yieldvalues[k] ;
+
+ // covariance weighted pdf for each specief
+ Double_t sweightsum(0) ;
+ for(Int_t n=0; n<nspec; ++n) {
+ Double_t nsum(0) ;
+ for(Int_t j=0; j<nspec; ++j) nsum += covMatrix(n,j) * pdfvalues[ievt][j] ;
+ sweightvec[n]->setVal(nsum/dsum) ;
+ pdfvec[n]->setVal( pdfvalues[ievt][n] ) ;
+ sweightsum+= nsum/dsum ;
+ if( !(fabs(nsum/dsum)>=0 ) ) {
+ std::cout << "error: " << nsum/dsum << std::endl ;
+ return 0 ;
+ }
+
+ //std::cout << nsum/dsum << " " ;
+ }
+ sWeightData->add(sweightset) ;
+ //std::cout << "sum : " << sweightsum << std::endl ;
+ }
+ std::cout << "number of entries in new dataset: " << data.numEntries() << " "
+ << sWeightData->numEntries() << std::endl ;
+
+ //RooDataSet mergeddata;
+ //mergeddata.merge(&data,&sWeightData) ;
+ //data.merge(&sWeightData);
+ //std::cout << "number of entries in final dataset: " << data.numEntries() << std::endl ;
+ return sWeightData ;
+}
+
+//____________________________________________________________________
+void SPlot::FillSPlot(const RooDataSet &data, TString varname, TString weightname)
+{
+ // Method to fill an SPlot for a given variable varname
+ if (data.get()->find(varname) == NULL) {
+ std::cout << "Can't find variable " << varname << " in data set!" << std::endl;
+ return;
+ }
+
+ if (data.get()->find(weightname) == NULL){
+ std::cout << "Can't find weight " << weightname << " in data set!" << std::endl;
+ return;
+ }
+
+
+ for (Int_t ievt = 0; ievt < data.numEntries(); ievt++) {
+ RooArgList row(*data.get(ievt));
+ Double_t xval = ((RooAbsReal*)row.find(varname))->getVal();
+
+ Double_t p = ((RooAbsReal*)row.find(weightname))->getVal();
+
+ Fill(xval,p);
+ }
+
+
+}
+
+
+
+//____________________________________________________________________
+void SPlot::FillSPlot(const RooAbsReal &x, RooAbsReal &nstar, RooDataSet data, const RooFitResult &fitRes, const RooArgList &pdfListTmp, const RooArgList &yieldsTmp, RooAbsPdf &totalPdf, Bool_t doErrors, const RooArgSet &projDeps)
+{
+ // Alternate method to fill an SPlot for the variable x. Better description of this method is needed.
+
+ Bool_t verbose = kTRUE;
+ if (verbose) {
+ std::cout << "yieldsTmp:" << std::endl;
+ yieldsTmp.Print("V");
+ std::cout << "pdfListTmp:" << std::endl;
+ pdfListTmp.Print();
+ }
+
+ // Size of bins...
+ Double_t xmin = GetXaxis()->GetXmin();
+ Double_t xmax = GetXaxis()->GetXmax();
+ Double_t nbins = GetNbinsX();
+ Double_t binSize = (xmax - xmin)/nbins;
+
+ if (verbose) {
+ std::cout << "Bins: " << xmin << " to " << xmax << " with " << nbins << " bins." << std::endl;
+ std::cout << "binSize = " << binSize << std::endl;
+ }
+
+ // Number of species in this fit.
+ Int_t nspec = yieldsTmp.getSize();
+
+ // The list of parameters (with their final fit values)
+ RooArgList finalPars = fitRes.floatParsFinal();
+
+ // Number of parameters in the fit result.
+ Int_t npars = finalPars.getSize();
+
+ RooArgList pdfList = *(RooArgList*)pdfListTmp.snapshot();
+ RooArgList yields = *(RooArgList*)yieldsTmp.snapshot(kFALSE);
+ //RooAbsPdf totalPdf = *(RooAbsPdf*)totalPdfTmp.Clone();
+
+ if (verbose) {
+ std::cout << "Yields I will use in calculation:" << std::endl;
+ yields.Print("V");
+ std::cout << "pdfList:" << std::endl;
+ pdfList.Print();
+ }
+
+ // The list of variables to normalize over when calculating PDF values.
+ RooArgSet vars(*data.get());
+ vars.remove(projDeps, kTRUE, kTRUE);
+
+ // Make a TMatrixD to hold the covariance matrix.
+ TMatrixD covMatrix(npars, npars);
+
+ // Loop over all the parameters to make the covariance matrix.
+ for (Int_t i = 0; i < npars; i++) {
+ for (Int_t j = 0; j < npars; j++) {
+ const RooRealVar *rowVar= (const RooRealVar*)finalPars.at(i);
+ const RooRealVar *colVar= (const RooRealVar*)finalPars.at(j);
+ assert(0 != rowVar && 0 != colVar);
+ covMatrix(i,j) = rowVar->getError()*colVar->getError()*fitRes.correlation(rowVar->GetName(),colVar->GetName());
+ }
+ }
+
+ // Get the inverse of the covariance matrix
+ // First check if it's singular
+ if (covMatrix.Determinant() == 0) {
+ std::cout << "SPlot Error: covariance matrix is singular; I can't invert it!" << std::endl;
+ covMatrix.Print();
+ return;
+ }
+ TMatrixD covInv(TMatrixD::kInverted,covMatrix);
+
+ if (verbose) {
+ std::cout << "Covariance matrix:" << std::endl;
+ covMatrix.Print();
+ std::cout << "Inverse of covariance matrix:" << std::endl;
+ covInv.Print();
+ }
+
+ // Make a matrix to hold V(i,j) inverse.
+ TMatrixD vinv(nspec, nspec);
+
+ // And fill it with the correct numbers...
+ Int_t istar(0);
+ Int_t vi = 0;
+ for (Int_t ci = 0; ci < npars; ci++) {
+ // If this parameter isn't a yield, move to the next row
+ if (yields.find(finalPars.at(ci)->GetName()) == 0) continue;
+
+ // If this parameter is the one of interest (nstar), then remember its index.
+ TString name = ((RooRealVar*)finalPars.at(ci))->GetName();
+ if (!name.CompareTo(nstar.GetName())) {
+ istar = vi;
+ }
+
+ Int_t vj = 0;
+ for (Int_t cj = 0; cj < npars; cj++) {
+
+ // If this parameter isn't a yield, move to the next column
+ if (!yields.contains(*finalPars.at(cj))) continue;
+
+ // This element's row and column correspond to yield parameters. Put it in V inverse.
+ vinv(vi, vj) = covInv(ci, cj);
+ vj++;
+ }
+ vi++;
+ }
+
+ // Now invert V(i, j) inverse to get V(i, j)
+ if (vinv.Determinant() == 0) {
+ std::cout << "SPlot Error: Yield covariance matrix V inverse is singular and can't be inverted!" << std::endl;
+ vinv.Print();
+ return;
+ }
+
+ TMatrixD v(TMatrixD::kInverted,vinv);
+
+ if (verbose) {
+ std::cout << "V inverse:" << std::endl;
+ vinv.Print();
+ std::cout << "V:" << std::endl;
+ v.Print();
+
+ Double_t sum = 0;
+ for (Int_t j = 0; j < nspec; j++) {
+ sum += v(j,istar);
+ }
+ std::cout << "Sum of star column in V: " << sum << std::endl;
+
+ }
+
+
+
+ Double_t sum = 0;
+ // Double_t sumtmp = 0;
+
+ // This forces the error in a bin to be calculated as the sqrt of the sum of the squares of
+ // weights, as they should be.
+ if (doErrors) Sumw2();
+
+ totalPdf.attachDataSet(data);
+ for (Int_t ievt = 0; ievt < data.numEntries(); ievt++) {
+
+ if (ievt % 100 == 0) std::cout << "Event: " << ievt << std::endl;
+ // Read this event and find the value of x for this event.
+ const RooArgSet *row = data.get(ievt);
+ Double_t xval = ((RooAbsReal*)row->find(x.GetName()))->getVal();
+
+ // Loop over the species and calculate P.
+ Double_t numerator = 0;
+ Double_t denominator = 0;
+ for (Int_t i = 0; i < nspec ; i++){
+
+ RooAbsPdf *pdf = (RooAbsPdf*)pdfList.at(i);
+ pdf->attachDataSet(data);
+ Double_t pdfval = pdf->getVal(&vars);
+ numerator += v(istar, i)*pdfval;
+
+ //Double_t ni = ((RooAbsReal*)yields.at(i))->getVal();
+ //denominator += ni*pdfval;
+ }
+ denominator = totalPdf.getVal(&vars);
+
+ Double_t p = 1/nstar.getVal()*(numerator/denominator);
+ if (xval > xmin && xval < xmax) sum += p*nstar.getVal();
+
+ Fill(xval, p*nstar.getVal());
+
+ }// end event loop
+
+ SetEntries(sum*Double_t(binSize));
+ if (verbose) std::cout << "Entries should be: " << sum*Double_t(binSize) << " (" << sum << " events)" << std::endl;
+
+}
+
+//____________________________________________________________________
+void SPlot::FillSPlot(const RooAbsReal &x, RooAbsReal &nstar, RooDataSet data, const RooFitResult &fitRes, const RooArgList &pdfListTmp, const RooArgList &yieldsTmp, Bool_t doErrors, const RooArgSet &projDeps)
+{
+ // Alternate method to fill an SPlot for the variable x. Better description of this method is needed.
+
+ Bool_t verbose = kTRUE;
+ if (verbose) {
+ std::cout << "nstar: " << std::endl;
+ nstar.Print();
+ std::cout << "yieldsTmp:" << std::endl;
+ yieldsTmp.Print("V");
+ std::cout << "pdfListTmp:" << std::endl;
+ pdfListTmp.Print();
+ }
+
+ // Size of bins...
+ Double_t xmin = GetXaxis()->GetXmin();
+ Double_t xmax = GetXaxis()->GetXmax();
+ Double_t nbins = GetNbinsX();
+ Double_t binSize = (xmax - xmin)/nbins;
+
+ if (verbose) {
+ std::cout << "Bins: " << xmin << " to " << xmax << " with " << nbins << " bins." << std::endl;
+ std::cout << "binSize = " << binSize << std::endl;
+ }
+
+ // Number of species in this fit.
+ Int_t nspec = yieldsTmp.getSize();
+
+ // The list of parameters (with their final fit values)
+ RooArgList finalPars = fitRes.floatParsFinal();
+
+ // Number of parameters in the fit result.
+ Int_t npars = finalPars.getSize();
+
+ RooArgList pdfList = *(RooArgList*)pdfListTmp.snapshot();
+ RooArgList yields = *(RooArgList*)yieldsTmp.snapshot(kFALSE);
+
+ if (verbose) {
+ std::cout << "Yields I will use in calculation:" << std::endl;
+ yields.Print("V");
+ std::cout << "pdfList:" << std::endl;
+ pdfList.Print();
+ }
+
+ // The list of variables to normalize over when calculating PDF values.
+ RooArgSet vars(*data.get());
+ vars.remove(projDeps, kTRUE, kTRUE);
+
+ // Make a TMatrixD to hold the covariance matrix.
+ TMatrixD covMatrix(npars, npars);
+
+ // Loop over all the parameters to make the covariance matrix.
+ for (Int_t i = 0; i < npars; i++) {
+ for (Int_t j = 0; j < npars; j++) {
+ const RooRealVar *rowVar= (const RooRealVar*)finalPars.at(i);
+ const RooRealVar *colVar= (const RooRealVar*)finalPars.at(j);
+ assert(0 != rowVar && 0 != colVar);
+ covMatrix(i,j) = rowVar->getError()*colVar->getError()*fitRes.correlation(rowVar->GetName(),colVar->GetName());
+ }
+ }
+
+ // Get the inverse of the covariance matrix
+ // First check if it's singular
+ if (covMatrix.Determinant() == 0) {
+ std::cout << "SPlot Error: covariance matrix is singular; I can't invert it!" << std::endl;
+ covMatrix.Print();
+ return;
+ }
+ TMatrixD covInv(TMatrixD::kInverted,covMatrix);
+
+ if (verbose) {
+ std::cout << "Covariance matrix:" << std::endl;
+ covMatrix.Print();
+ std::cout << "Inverse of covariance matrix:" << std::endl;
+ covInv.Print();
+ }
+
+ // Make a matrix to hold V(i,j) inverse.
+ TMatrixD vinv(nspec, nspec);
+
+ // And fill it with the correct numbers...
+ Int_t istar(0);
+ Int_t vi = 0;
+ for (Int_t ci = 0; ci < npars; ci++) {
+ // If this parameter isn't a yield, move to the next row
+ if (yields.find(finalPars.at(ci)->GetName()) == 0) continue;
+
+ // If this parameter is the one of interest (nstar), then remember its index.
+ TString name = ((RooRealVar*)finalPars.at(ci))->GetName();
+ if (!name.CompareTo(nstar.GetName())) {
+ istar = vi;
+ }
+
+ Int_t vj = 0;
+ for (Int_t cj = 0; cj < npars; cj++) {
+
+ // If this parameter isn't a yield, move to the next column
+ if (!yields.contains(*finalPars.at(cj))) continue;
+
+ // This element's row and column correspond to yield parameters. Put it in V inverse.
+ vinv(vi, vj) = covInv(ci, cj);
+ vj++;
+ }
+ vi++;
+ }
+
+ // Now invert V(i, j) inverse to get V(i, j)
+ if (vinv.Determinant() == 0) {
+ std::cout << "SPlot Error: Yield covariance matrix V inverse is singular and can't be inverted!" << std::endl;
+ vinv.Print();
+ return;
+ }
+
+ TMatrixD v(TMatrixD::kInverted,vinv);
+
+ if (verbose) {
+ std::cout << "V inverse:" << std::endl;
+ vinv.Print();
+ std::cout << "V:" << std::endl;
+ v.Print();
+
+ Double_t sum = 0;
+ for (Int_t j = 0; j < nspec; j++) {
+ sum += v(j,istar);
+ }
+ std::cout << "Sum of star column in V: " << sum << std::endl;
+
+ }
+
+
+
+ Double_t sum = 0;
+ Double_t sumtmp = 0;
+
+ // This forces the error in a bin to be calculated as the sqrt of the sum of the squares of
+ // weights, as they should be.
+ if (doErrors) Sumw2();
+
+ for (Int_t i = 0; i < nspec; i++) {
+ RooAbsPdf *pdf = (RooAbsPdf*)pdfList.at(i);
+ pdf->attachDataSet(data);
+ }
+
+ for (Int_t ievt = 0; ievt < data.numEntries(); ievt++) {
+
+ if (ievt % 100 == 0) std::cout << ".";
+ // Read this event and find the value of x for this event.
+ const RooArgSet *row = data.get(ievt);
+ Double_t xval = ((RooAbsReal*)row->find(x.GetName()))->getVal();
+
+ // Loop over the species and calculate P.
+ Double_t numerator = 0;
+ Double_t denominator = 0;
+ for (Int_t i = 0; i < nspec ; i++){
+
+ RooAbsPdf *pdf = (RooAbsPdf*)pdfList.at(i);
+ //pdf->attachDataSet(data);
+ Double_t pdfval = pdf->getVal(&vars);
+ numerator += v(istar, i)*pdfval;
+
+ Double_t ni = ((RooAbsReal*)yields.at(i))->getVal();
+ denominator += ni*pdfval;
+ }
+
+ Double_t p = 1/nstar.getVal()*(numerator/denominator);
+ sumtmp += ((RooAbsPdf*)pdfList.at(istar))->getVal(&vars)/denominator;
+
+ //if (xval > xmin && xval < xmax)
+ sum += p*nstar.getVal();
+
+ Fill(xval, p*nstar.getVal());
+
+ }// end event loop
+
+ SetEntries(sum*Double_t(binSize));
+
+ std::cout << std::endl;
+ if (verbose) std::cout << "Entries should be: " << sum*Double_t(binSize) << " (" << sum << " events)" << std::endl;
+ if (verbose) std::cout << "Sum of likelihood ratios for nstar: " << sumtmp << std::endl;
+
+}
+
+
+//____________________________________________________________________
+void SPlot::FillSPlot(const RooAbsReal &x, RooAbsReal &nstar, RooDataSet data, const RooFitResult &fitRes, RooAbsPdf &totalPdf, RooArgList &yields, Bool_t doErrors, const RooArgSet &projDeps)
+{
+ // Alternate method to fill an SPlot for the variable x. Better description of this method is needed.
+
+ Bool_t verbose = kTRUE;
+ if (verbose) {
+ //std::cout << "yieldsTmp:" << std::endl;
+ //yieldsTmp.Print("V");
+ //std::cout << "pdfListTmp:" << std::endl;
+ //pdfListTmp.Print();
+ }
+
+ // Size of bins...
+ Double_t xmin = GetXaxis()->GetXmin();
+ Double_t xmax = GetXaxis()->GetXmax();
+ Double_t nbins = GetNbinsX();
+ Double_t binSize = (xmax - xmin)/nbins;
+
+ if (verbose) {
+ std::cout << "Bins: " << xmin << " to " << xmax << " with " << nbins << " bins." << std::endl;
+ std::cout << "binSize = " << binSize << std::endl;
+ }
+
+ // Number of species in this fit.
+ Int_t nspec = yields.getSize();
+
+ // The list of parameters (with their final fit values)
+ RooArgList finalPars = fitRes.floatParsFinal();
+
+ // Number of parameters in the fit result.
+ Int_t npars = finalPars.getSize();
+
+ //RooArgList pdfList = *(RooArgList*)pdfListTmp.snapshot();
+ //RooArgList yields = *(RooArgList*)yieldsTmp.snapshot(kFALSE);
+ //RooAbsPdf totalPdf = *(RooAbsPdf*)totalPdfTmp.Clone();
+
+ if (verbose) {
+ std::cout << "Yields I will use in calculation:" << std::endl;
+ yields.Print("V");
+ //std::cout << "pdfList:" << std::endl;
+ //pdfList.Print();
+ }
+
+ // The list of variables to normalize over when calculating PDF values.
+ RooArgSet vars(*data.get());
+ vars.remove(projDeps, kTRUE, kTRUE);
+
+ // Make a TMatrixD to hold the covariance matrix.
+ TMatrixD covMatrix(npars, npars);
+
+ // Loop over all the parameters to make the covariance matrix.
+ for (Int_t i = 0; i < npars; i++) {
+ for (Int_t j = 0; j < npars; j++) {
+ const RooRealVar *rowVar= (const RooRealVar*)finalPars.at(i);
+ const RooRealVar *colVar= (const RooRealVar*)finalPars.at(j);
+ assert(0 != rowVar && 0 != colVar);
+ covMatrix(i,j) = rowVar->getError()*colVar->getError()*fitRes.correlation(rowVar->GetName(),colVar->GetName());
+ }
+ }
+
+ // Get the inverse of the covariance matrix
+ // First check if it's singular
+ if (covMatrix.Determinant() == 0) {
+ std::cout << "SPlot Error: covariance matrix is singular; I can't invert it!" << std::endl;
+ covMatrix.Print();
+ return;
+ }
+ TMatrixD covInv(TMatrixD::kInverted,covMatrix);
+
+ if (verbose) {
+ std::cout << "Covariance matrix:" << std::endl;
+ covMatrix.Print();
+ std::cout << "Inverse of covariance matrix:" << std::endl;
+ covInv.Print();
+ }
+
+ // Make a matrix to hold V(i,j) inverse.
+ TMatrixD vinv(nspec, nspec);
+
+ // And fill it with the correct numbers...
+ Int_t istar(0);
+ Int_t vi = 0;
+ for (Int_t ci = 0; ci < npars; ci++) {
+ // If this parameter isn't a yield, move to the next row
+ if (yields.find(finalPars.at(ci)->GetName()) == 0) continue;
+
+ // If this parameter is the one of interest (nstar), then remember its index.
+ TString name = ((RooRealVar*)finalPars.at(ci))->GetName();
+ if (!name.CompareTo(nstar.GetName())) {
+ istar = vi;
+ }
+
+ Int_t vj = 0;
+ for (Int_t cj = 0; cj < npars; cj++) {
+
+ // If this parameter isn't a yield, move to the next column
+ if (!yields.contains(*finalPars.at(cj))) continue;
+
+ // This element's row and column correspond to yield parameters. Put it in V inverse.
+ vinv(vi, vj) = covInv(ci, cj);
+ vj++;
+ }
+ vi++;
+ }
+
+ // Now invert V(i, j) inverse to get V(i, j)
+ if (vinv.Determinant() == 0) {
+ std::cout << "SPlot Error: Yield covariance matrix V inverse is singular and can't be inverted!" << std::endl;
+ vinv.Print();
+ return;
+ }
+
+ TMatrixD v(TMatrixD::kInverted,vinv);
+
+ if (verbose) {
+ std::cout << "V inverse:" << std::endl;
+ vinv.Print();
+ std::cout << "V:" << std::endl;
+ v.Print();
+
+ Double_t sum = 0;
+ for (Int_t j = 0; j < nspec; j++) {
+ sum += v(j,istar);
+ }
+ std::cout << "Sum of star column in V: " << sum << std::endl;
+
+ }
+
+
+
+ Double_t sum = 0;
+ //Double_t sumtmp = 0;
+
+ // This forces the error in a bin to be calculated as the sqrt of the sum of the squares of
+ // weights, as they should be.
+ if (doErrors) Sumw2();
+
+ totalPdf.attachDataSet(data);
+ for (Int_t ievt = 0; ievt < data.numEntries(); ievt++) {
+
+ if (ievt % 100 == 0) std::cout << "Event: " << ievt << std::endl;
+ // Read this event and find the value of x for this event.
+ const RooArgSet *row = data.get(ievt);
+ Double_t xval = ((RooAbsReal*)row->find(x.GetName()))->getVal();
+
+ // Loop over the species and calculate P.
+ Double_t numerator = 0;
+ Double_t denominator = 0;
+ for (Int_t i = 0; i < nspec ; i++){
+
+ //RooAbsPdf *pdf = (RooAbsPdf*)pdfList.at(i);
+ //pdf->attachDataSet(data);
+ //Double_t pdfval = pdf->getVal(&vars);
+ Double_t pdfval = GetComponentValue(totalPdf, yields, i, vars);
+ numerator += v(istar, i)*pdfval;
+
+ //Double_t ni = ((RooAbsReal*)yields.at(i))->getVal();
+ //denominator += ni*pdfval;
+ }
+ denominator = totalPdf.getVal(&vars);
+
+ Double_t p = 1/nstar.getVal()*(numerator/denominator);
+ if (xval > xmin && xval < xmax) sum += p*nstar.getVal();
+
+ Fill(xval, p*nstar.getVal());
+
+ }// end event loop
+
+ SetEntries(sum*Double_t(binSize));
+ if (verbose) std::cout << "Entries should be: " << sum*Double_t(binSize) << " (" << sum << " events)" << std::endl;
+
+}
+
+
+//____________________________________________________________________
+Double_t SPlot::GetComponentValue(RooAbsPdf &pdf, RooArgList &yieldsTmp, Int_t igood, RooArgSet &normSet)
+{
+ // Alternate method to fill an SPlot for the variable x. Better description of this method is needed.
+
+ Int_t i=0;
+ Int_t nspec = yieldsTmp.getSize();
+ std::vector<Double_t> yields(nspec);
+ // std::cout << "Before: " << pdf.getVal(&normSet) << std::endl;
+ for (i=0; i < nspec; i++) {
+ yields[i] = ((RooRealVar*)yieldsTmp.at(i))->getVal();
+ ((RooRealVar*)yieldsTmp.at(i))->setVal(0);
+ }
+
+ ((RooRealVar*)yieldsTmp.at(igood))->setVal(1);
+
+ Double_t result = pdf.getVal(&normSet);
+ // std::cout << "During: " << result << std::endl;
+
+ for (i=0; i < nspec; i++) {
+ ((RooRealVar*)yieldsTmp.at(i))->setVal(yields[i]);
+ }
+ // std::cout << "After: " << pdf.getVal(&normSet) << std::endl;
+ return result;
+
+}
+
+
+
+
+
--
GitLab