This tutorial shows a more complex example using the FeldmanCousins utility to create a confidence interval for a toy neutrino oscillation experiment. The example attempts to faithfully reproduce the toy example described in Feldman & Cousins' original paper, Phys.Rev.D57:3873-3889,1998.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(PnmuTone_Int[L]) using numeric integrator RooRombergIntegrator to calculate Int(L)
generate toy data with nEvents = 692
Minuit2Minimizer: Minimize with max-calls 1000 convergence for edm < 1 strategy 1
Minuit2Minimizer : Valid minimum - status = 0
FVAL = -1131.34224753934996
Edm = 4.26659709371356039e-08
Nfcn = 33
deltaMSq = 37.5389 +/- 4.12974 (limited)
sinSq2theta = 0.00629097 +/- 0.000861732 (limited)
[#1] INFO:Minimization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(PnmuTonePrime_Int[EPrime,LPrime]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(LPrime,EPrime)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(PnmuTone_Int[E,L]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(L,E)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(PnmuTone_Int[L]_Norm[E,L]) using numeric integrator RooRombergIntegrator to calculate Int(L)
Metropolis-Hastings progress: ....................................................................................................
[#1] INFO:Eval -- Proposal acceptance rate: 3.3%
[#1] INFO:Eval -- Number of steps in chain: 165
[#1] INFO:NumericIntegration -- RooRealIntegral::init(product_Int[deltaMSq,sinSq2theta]_Norm[deltaMSq,sinSq2theta]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(deltaMSq,sinSq2theta)
[#0] WARNING:NumericIntegration -- RooAdaptiveIntegratorND::dtor(product) WARNING: Number of suppressed warningings about integral evaluations where target precision was not reached is 1
[#1] INFO:NumericIntegration -- RooRealIntegral::init(product_Int[deltaMSq,sinSq2theta]_Norm[deltaMSq,sinSq2theta]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(deltaMSq,sinSq2theta)
[#1] INFO:Eval -- cutoff = 0.166573, conf = 0.904333
[#0] WARNING:NumericIntegration -- RooAdaptiveIntegratorND::dtor(product) WARNING: Number of suppressed warningings about integral evaluations where target precision was not reached is 1
[#0] WARNING:NumericIntegration -- RooAdaptiveIntegratorND::dtor(PnmuTone) WARNING: Number of suppressed warningings about integral evaluations where target precision was not reached is 628
[#0] WARNING:NumericIntegration -- RooAdaptiveIntegratorND::dtor(PnmuTonePrime) WARNING: Number of suppressed warningings about integral evaluations where target precision was not reached is 628
Real time 0:02:33, CP time 153.430
MCMC actual confidence level: 0.904333
[#1] INFO:Minimization -- RooProfileLL::evaluate(nll_model_modelData_Profile[deltaMSq,sinSq2theta]) Creating instance of MINUIT
[#1] INFO:Minimization -- RooProfileLL::evaluate(nll_model_modelData_Profile[deltaMSq,sinSq2theta]) determining minimum likelihood for current configurations w.r.t all observable
[#1] INFO:Minimization -- RooProfileLL::evaluate(nll_model_modelData_Profile[deltaMSq,sinSq2theta]) minimum found at (deltaMSq=37.5376, sinSq2theta=0.00629099)
..[#1] INFO:Minimization -- LikelihoodInterval - Finding the contour of deltaMSq ( 0 ) and sinSq2theta ( 1 )
Real time 0:03:02, CP time 181.990
#include <iostream>
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "../tutorials/roostats/NuMuToNuE_Oscillation.h"
#include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx"
#else
#include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx+"
#endif
void rs401d_FeldmanCousins(bool doFeldmanCousins = false, bool doMCMC = true)
{
RooRealVar deltaMSq(
"deltaMSq",
"#Delta m^{2}", 40, 1, 300,
"eV/c^{2}");
RooRealVar sinSq2theta(
"sinSq2theta",
"sin^{2}(2#theta)", .006, .0, .02);
RooRealVar EPrime(
"EPrime",
"", 15, 10, 60,
"GeV");
RooRealVar LPrime(
"LPrime",
"", .800, .600, 1.0,
"km");
NuMuToNuE_Oscillation PnmuTonePrime(
"PnmuTonePrime",
"P(#nu_{#mu} #rightarrow #nu_{e}", LPrime, EPrime, deltaMSq);
RooConstVar maxEventsTot(
"maxEventsTot",
"maximum number of sinal events", 50000);
RooConstVar inverseArea(
"inverseArea",
"1/(#Delta E #Delta L)",
1. / (EPrime.getMax() - EPrime.getMin()) / (LPrime.getMax() - LPrime.getMin()));
RooProduct sigNorm(
"sigNorm",
"",
RooArgSet(maxEventsTot, *intProbToOscInExp, inverseArea, sinSq2theta));
RooConstVar bkgNorm(
"bkgNorm",
"normalization for background", 500);
Int_t nEventsData = bkgNorm.getVal() + sigNorm.getVal();
cout << "generate toy data with nEvents = " << nEventsData << endl;
TH1 *hh = PnmuTone.createHistogram(
"hh", E, Binning(40), YVar(L, Binning(40)), Scaling(
kFALSE));
model.fitTo(*
data, Extended());
model.plotOn(Eframe);
model.plotOn(Eframe, Components(*sigModel), LineColor(
kRed));
model.plotOn(Eframe, Components(bkgEShape), LineColor(
kGreen));
model.plotOn(Eframe);
Eframe->
SetTitle(
"toy data with best fit model (and sig+bkg components)");
std::unique_ptr<RooAbsReal> nll{model.createNLL(*
data, Extended(
true))};
TH1 *hhh = pll.createHistogram(
"hhh", sinSq2theta, Binning(40), YVar(deltaMSq, Binning(40)), Scaling(
kFALSE));
fc.SetTestSize(.1);
fc.UseAdaptiveSampling(true);
fc.SetNBins(10);
if (doFeldmanCousins)
interval = fc.GetInterval();
plc.SetTestSize(.1);
if (doMCMC) {
mc.SetNumIters(5000);
mc.SetNumBurnInSteps(100);
mc.SetUseKeys(true);
mc.SetTestSize(.1);
mc.SetAxes(axisList);
}
if (doFeldmanCousins) {
if (interval) {
} else {
}
}
delete tmpPoint;
}
if (interval) {
forContour->
Draw(
"cont2,same");
}
}
if (mcInt) {
}
plotInt.SetTitle("90% Confidence Intervals");
if (mcInt)
plotInt.Draw("same");
else
plotInt.Draw();
}
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
TH1 * createHistogram(const char *name, const RooAbsRealLValue &xvar, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Calls createHistogram(const char *name, const RooAbsRealLValue& xvar, const RooLinkedList& argList) c...
Abstract interface for all probability density functions.
virtual RooAbsPdf * createProjection(const RooArgSet &iset)
Return a p.d.f that represent a projection of this p.d.f integrated over given observables.
Abstract base class for objects that represent a real value and implements functionality common to al...
RooFit::OwningPtr< RooAbsReal > createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create an object that represents the integral of the function over one or more observables listed in ...
Efficient implementation of a sum of PDFs of the form.
RooArgList is a container object that can hold multiple RooAbsArg objects.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
TObject * clone(const char *newname) const override
RooConstVar represent a constant real-valued object.
The RooDataHist is a container class to hold N-dimensional binned data.
const RooArgSet * get() const override
Get bin centre of current bin.
RooDataSet is a container class to hold unbinned data.
static RooMsgService & instance()
Return reference to singleton instance.
void setStreamStatus(Int_t id, bool active)
(De)Activate stream with given unique ID
A RooPlot is a plot frame and a container for graphics objects within that frame.
void SetTitle(const char *name) override
Set the title of the RooPlot to 'title'.
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
RooPolynomial implements a polynomial p.d.f of the form.
A RooProduct represents the product of a given set of RooAbsReal objects.
Class RooProfileLL implements the profile likelihood estimator for a given likelihood and set of para...
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
RooRealVar represents a variable that can be changed from the outside.
ConfInterval is an interface class for a generic interval in the RooStats framework.
virtual bool IsInInterval(const RooArgSet &) const =0
check if given point is in the interval
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
This class provides simple and straightforward utilities to plot a LikelihoodInterval object.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
This class provides simple and straightforward utilities to plot a MCMCInterval object.
void SetLineColor(Color_t color)
void Draw(const Option_t *options=nullptr) override
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual double GetActualConfidenceLevel()
virtual double GetKeysPdfCutoff() { return fKeysCutoff; }
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
virtual void SetWorkspace(RooWorkspace &ws)
virtual void SetParametersOfInterest(const RooArgSet &set)
Specify parameters of interest.
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the workspace if not already there.
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
Persistable container for RooFit projects.
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
virtual void SetLineColor(Color_t lcolor)
Set the line color.
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
void Update() override
Update canvas pad buffers.
TH1 is the base class of all histogram classes in ROOT.
void SetTitle(const char *title) override
Change/set the title.
void Draw(Option_t *option="") override
Draw this histogram with options.
virtual void SetContour(Int_t nlevels, const Double_t *levels=nullptr)
Set the number and values of contour levels.
virtual Int_t FindBin(Double_t x, Double_t y=0, Double_t z=0)
Return Global bin number corresponding to x,y,z.
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
2-D histogram with a float per channel (see TH1 documentation)}
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void Stop()
Stop the stopwatch.
void Print(Option_t *option="") const override
Print the real and cpu time passed between the start and stop events.
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.