Logo ROOT  
Reference Guide
rs102_hypotestwithshapes.C File Reference

Detailed Description

A typical search for a new particle by studying an invariant mass distribution.

View in nbviewer Open in SWAN The macro creates a simple signal model and two background models, which are added to a RooWorkspace. The macro creates a toy dataset, and then uses a RooStats ProfileLikleihoodCalculator to do a hypothesis test of the background-only and signal+background hypotheses. In this example, shape uncertainties are not taken into account, but normalization uncertainties are.

RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooAddPdf::model
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooGaussian::sigModel
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::invMass
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::mH
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::sigma1
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooProduct::fsig
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::mu
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::ratioSigEff
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::fsigExpected
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooGaussian::zjjModel
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::mZ
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::sigma1_z
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::fzjj
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooChebychev::qcdModel
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::a0
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::a1
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::a2
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing dataset modelData
[#1] INFO:ObjectHandling -- RooWorkSpace::import(myWS) changing name of dataset from modelData to data
[#1] INFO:Minization -- createNLL: caching constraint set under name CONSTR_OF_PDF_model_FOR_OBS_invMass with 0 entries
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit / Migrad with strategy 1
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization -- The following expressions have been identified as constant and will be precalculated and cached: (sigModel,zjjModel,qcdModel)
[#1] INFO:Minization --
RooFitResult: minimized FCN value: 717.039, estimated distance to minimum: 8.90226e-10
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
fzjj 3.1152e-01 +/- 5.03e-02
mu 1.0968e+00 +/- 3.03e-01
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::GetHypoTest - do conditional fit
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit / Migrad with strategy 1
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization --
RooFitResult: minimized FCN value: 723.97, estimated distance to minimum: 2.09862e-09
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
fzjj 2.6213e-01 +/- 5.18e-02
-------------------------------------------------
The p-value for the null is 9.83108e-05
Corresponding to a significance of 3.72332
-------------------------------------------------
[#1] INFO:Minization -- createNLL picked up cached consraints from workspace with 0 entries
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization -- The following expressions have been identified as constant and will be precalculated and cached: (sigModel,zjjModel,qcdModel)
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (sigModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (zjjModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (qcdModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Minization -- createNLL picked up cached consraints from workspace with 0 entries
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization -- The following expressions have been identified as constant and will be precalculated and cached: (sigModel,zjjModel,qcdModel)
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (zjjModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (qcdModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooAddition.h"
#include "RooProduct.h"
#include "TCanvas.h"
#include "RooChebychev.h"
#include "RooAbsPdf.h"
#include "RooFit.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "RooAbsArg.h"
#include "RooWorkspace.h"
#include <string>
// use this order for safety on library loading
using namespace RooFit;
using namespace RooStats;
// see below for implementation
void AddModel(RooWorkspace *);
void AddData(RooWorkspace *);
void DoHypothesisTest(RooWorkspace *);
void MakePlots(RooWorkspace *);
//____________________________________
void rs102_hypotestwithshapes()
{
// The main macro.
// Create a workspace to manage the project.
RooWorkspace *wspace = new RooWorkspace("myWS");
// add the signal and background models to the workspace
AddModel(wspace);
// add some toy data to the workspace
AddData(wspace);
// inspect the workspace if you wish
// wspace->Print();
// do the hypothesis test
DoHypothesisTest(wspace);
// make some plots
MakePlots(wspace);
// cleanup
delete wspace;
}
//____________________________________
void AddModel(RooWorkspace *wks)
{
// Make models for signal (Higgs) and background (Z+jets and QCD)
// In real life, this part requires an intelligent modeling
// of signal and background -- this is only an example.
// set range of observable
Double_t lowRange = 60, highRange = 200;
// make a RooRealVar for the observable
RooRealVar invMass("invMass", "M_{inv}", lowRange, highRange, "GeV");
// --------------------------------------
// make a simple signal model.
RooRealVar mH("mH", "Higgs Mass", 130, 90, 160);
RooRealVar sigma1("sigma1", "Width of Gaussian", 12., 2, 100);
RooGaussian sigModel("sigModel", "Signal Model", invMass, mH, sigma1);
// we will test this specific mass point for the signal
mH.setConstant();
// and we assume we know the mass resolution
sigma1.setConstant();
// --------------------------------------
// make zjj model. Just like signal model
RooRealVar mZ("mZ", "Z Mass", 91.2, 0, 100);
RooRealVar sigma1_z("sigma1_z", "Width of Gaussian", 10., 6, 100);
RooGaussian zjjModel("zjjModel", "Z+jets Model", invMass, mZ, sigma1_z);
// we know Z mass
mZ.setConstant();
// assume we know resolution too
sigma1_z.setConstant();
// --------------------------------------
// make QCD model
RooRealVar a0("a0", "a0", 0.26, -1, 1);
RooRealVar a1("a1", "a1", -0.17596, -1, 1);
RooRealVar a2("a2", "a2", 0.018437, -1, 1);
RooRealVar a3("a3", "a3", 0.02, -1, 1);
RooChebychev qcdModel("qcdModel", "A Polynomial for QCD", invMass, RooArgList(a0, a1, a2));
// let's assume this shape is known, but the normalization is not
a0.setConstant();
a1.setConstant();
a2.setConstant();
// --------------------------------------
// combined model
// Setting the fraction of Zjj to be 40% for initial guess.
RooRealVar fzjj("fzjj", "fraction of zjj background events", .4, 0., 1);
// Set the expected fraction of signal to 20%.
RooRealVar fsigExpected("fsigExpected", "expected fraction of signal events", .2, 0., 1);
fsigExpected.setConstant(); // use mu as main parameter, so fix this.
// Introduce mu: the signal strength in units of the expectation.
// eg. mu = 1 is the SM, mu = 0 is no signal, mu=2 is 2x the SM
RooRealVar mu("mu", "signal strength in units of SM expectation", 1, 0., 2);
// Introduce ratio of signal efficiency to nominal signal efficiency.
// This is useful if you want to do limits on cross section.
RooRealVar ratioSigEff("ratioSigEff", "ratio of signal efficiency to nominal signal efficiency", 1., 0., 2);
ratioSigEff.setConstant(kTRUE);
// finally the signal fraction is the product of the terms above.
RooProduct fsig("fsig", "fraction of signal events", RooArgSet(mu, ratioSigEff, fsigExpected));
// full model
RooAddPdf model("model", "sig+zjj+qcd background shapes", RooArgList(sigModel, zjjModel, qcdModel),
RooArgList(fsig, fzjj));
// interesting for debugging and visualizing the model
// model.printCompactTree("","fullModel.txt");
// model.graphVizTree("fullModel.dot");
wks->import(model);
}
//____________________________________
void AddData(RooWorkspace *wks)
{
// Add a toy dataset
Int_t nEvents = 150;
RooAbsPdf *model = wks->pdf("model");
RooRealVar *invMass = wks->var("invMass");
RooDataSet *data = model->generate(*invMass, nEvents);
wks->import(*data, Rename("data"));
}
//____________________________________
void DoHypothesisTest(RooWorkspace *wks)
{
// Use a RooStats ProfileLikleihoodCalculator to do the hypothesis test.
model.SetWorkspace(*wks);
model.SetPdf("model");
// plc.SetData("data");
plc.SetData(*(wks->data("data")));
// here we explicitly set the value of the parameters for the null.
// We want no signal contribution, eg. mu = 0
RooRealVar *mu = wks->var("mu");
// RooArgSet* nullParams = new RooArgSet("nullParams");
// nullParams->addClone(*mu);
RooArgSet poi(*mu);
RooArgSet *nullParams = (RooArgSet *)poi.snapshot();
nullParams->setRealValue("mu", 0);
// plc.SetNullParameters(*nullParams);
// NOTE: using snapshot will import nullparams
// in the WS and merge with existing "mu"
// model.SetSnapshot(*nullParams);
// use instead setNuisanceParameters
plc.SetNullParameters(*nullParams);
// We get a HypoTestResult out of the calculator, and we can query it.
cout << "-------------------------------------------------" << endl;
cout << "The p-value for the null is " << htr->NullPValue() << endl;
cout << "Corresponding to a significance of " << htr->Significance() << endl;
cout << "-------------------------------------------------\n\n" << endl;
}
//____________________________________
void MakePlots(RooWorkspace *wks)
{
// Make plots of the data and the best fit model in two cases:
// first the signal+background case
// second the background-only case.
// get some things out of workspace
RooAbsPdf *model = wks->pdf("model");
RooAbsPdf *sigModel = wks->pdf("sigModel");
RooAbsPdf *zjjModel = wks->pdf("zjjModel");
RooAbsPdf *qcdModel = wks->pdf("qcdModel");
RooRealVar *mu = wks->var("mu");
RooRealVar *invMass = wks->var("invMass");
RooAbsData *data = wks->data("data");
// --------------------------------------
// Make plots for the Alternate hypothesis, eg. let mu float
model->fitTo(*data, Save(kTRUE), Minos(kFALSE), Hesse(kFALSE), PrintLevel(-1));
// plot sig candidates, full model, and individual components
new TCanvas();
RooPlot *frame = invMass->frame();
data->plotOn(frame);
model->plotOn(frame);
model->plotOn(frame, Components(*sigModel), LineStyle(kDashed), LineColor(kRed));
model->plotOn(frame, Components(*zjjModel), LineStyle(kDashed), LineColor(kBlack));
model->plotOn(frame, Components(*qcdModel), LineStyle(kDashed), LineColor(kGreen));
frame->SetTitle("An example fit to the signal + background model");
frame->Draw();
// cdata->SaveAs("alternateFit.gif");
// --------------------------------------
// Do Fit to the Null hypothesis. Eg. fix mu=0
mu->setVal(0); // set signal fraction to 0
mu->setConstant(kTRUE); // set constant
model->fitTo(*data, Save(kTRUE), Minos(kFALSE), Hesse(kFALSE), PrintLevel(-1));
// plot signal candidates with background model and components
new TCanvas();
RooPlot *xframe2 = invMass->frame();
model->plotOn(xframe2);
model->plotOn(xframe2, Components(*zjjModel), LineStyle(kDashed), LineColor(kBlack));
model->plotOn(xframe2, Components(*qcdModel), LineStyle(kDashed), LineColor(kGreen));
xframe2->SetTitle("An example fit to the background-only model");
xframe2->Draw();
// cbkgonly->SaveAs("nullFit.gif");
}
Author
Kyle Cranmer

Definition in file rs102_hypotestwithshapes.C.

RooWorkspace::data
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1368
RooAbsRealLValue::frame
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
Definition: RooAbsRealLValue.cxx:199
RooStats::HypoTestResult::Significance
virtual Double_t Significance() const
familiar name for the Null p-value in terms of 1-sided Gaussian significance
Definition: HypoTestResult.h:80
RooWorkspace.h
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
RooArgSet::setRealValue
Bool_t setRealValue(const char *name, Double_t newVal=0, Bool_t verbose=kFALSE)
Set value of a RooAbsRealLValye stored in set with given name to newVal No error messages are printed...
Definition: RooArgSet.cxx:495
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:216
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooChebychev.h
RooAddPdf
Definition: RooAddPdf.h:32
RooAbsData
Definition: RooAbsData.h:46
RooAbsData::SumW2
@ SumW2
Definition: RooAbsData.h:96
RooFit.h
kGreen
@ kGreen
Definition: Rtypes.h:66
RooFit::DataError
RooCmdArg DataError(Int_t)
Definition: RooGlobalFunc.cxx:156
RooFit::PrintLevel
RooCmdArg PrintLevel(Int_t code)
Definition: RooGlobalFunc.cxx:189
RooStats::ProfileLikelihoodCalculator
Definition: ProfileLikelihoodCalculator.h:28
RooChebychev
Definition: RooChebychev.h:25
RooStats::HypoTestResult::NullPValue
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
Definition: HypoTestResult.h:57
RooArgList
Definition: RooArgList.h:21
RooGaussian.h
RooFit::Hesse
RooCmdArg Hesse(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:193
RooFit::Minos
RooCmdArg Minos(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:194
RooStats::ProfileLikelihoodCalculator::GetHypoTest
virtual HypoTestResult * GetHypoTest() const
Return the hypothesis test result obtained from the likelihood ratio of the maximum likelihood value ...
Definition: ProfileLikelihoodCalculator.cxx:301
RooAbsArg.h
RooGaussian
Definition: RooGaussian.h:25
RooAddPdf.h
TCanvas.h
HypoTestResult.h
RooDataSet.h
RooWorkspace::import
Bool_t import(const RooAbsArg &arg, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Import a RooAbsArg object, e.g.
Definition: RooWorkspace.cxx:361
RooArgSet::snapshot
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:134
RooProduct
Definition: RooProduct.h:30
kBlack
@ kBlack
Definition: Rtypes.h:65
RooFit::Rename
RooCmdArg Rename(const char *suffix)
Definition: RooGlobalFunc.cxx:315
RooStats::CombinedCalculator::SetModel
virtual void SetModel(const ModelConfig &model)
set the model (in this case can set only the parameters for the null hypothesis)
Definition: CombinedCalculator.h:127
RooProdPdf.h
RooStats::CombinedCalculator::SetNullParameters
virtual void SetNullParameters(const RooArgSet &set)
set parameter values for the null if using a common PDF
Definition: CombinedCalculator.h:153
RooFit
Definition: RooCFunction1Binding.h:29
RooAbsPdf.h
RooAbsData::plotOn
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Definition: RooAbsData.cxx:547
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooPlot.h
RooPlot
Definition: RooPlot.h:44
RooWorkspace::pdf
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1277
RooAddition.h
RooStats::CombinedCalculator::SetData
virtual void SetData(RooAbsData &data)
Set the DataSet, add to the the workspace if not already there.
Definition: CombinedCalculator.h:122
RooRealVar.h
kRed
@ kRed
Definition: Rtypes.h:66
RooProduct.h
RooFitResult.h
RooWorkspace
Definition: RooWorkspace.h:43
RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition: RooGlobalFunc.cxx:56
Double_t
double Double_t
Definition: RtypesCore.h:59
RooPlot::SetTitle
void SetTitle(const char *name)
Set the title of the RooPlot to 'title'.
Definition: RooPlot.cxx:1237
TCanvas
Definition: TCanvas.h:23
RooStats
Definition: Asimov.h:19
ProfileLikelihoodCalculator.h
RooStats::HypoTestResult
Definition: HypoTestResult.h:28
RooAbsRealLValue::setConstant
void setConstant(Bool_t value=kTRUE)
Definition: RooAbsRealLValue.h:115
kDashed
@ kDashed
Definition: TAttLine.h:48
RooWorkspace::var
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1295
RooDataSet
Definition: RooDataSet.h:33
make_cnn_model.model
model
Definition: make_cnn_model.py:6
RooAbsPdf
Definition: RooAbsPdf.h:40
RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition: RooGlobalFunc.cxx:57
RooRealVar
Definition: RooRealVar.h:36
RooFit::Components
RooCmdArg Components(const RooArgSet &compSet)
Definition: RooGlobalFunc.cxx:74
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:187
RooStats::ModelConfig
Definition: ModelConfig.h:36
RooArgSet
Definition: RooArgSet.h:28
int