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.
#include <iostream>
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);
auto oscillationFormula = "std::pow(std::sin(1.27 * x[2] * x[0] / x[1]), 2)";
RooGenericPdf PnmuTone(
"PnmuTone",
"P(#nu_{#mu} #rightarrow #nu_{e}", oscillationFormula, {L, E, deltaMSq});
RooRealVar EPrime(
"EPrime",
"", 15, 10, 60,
"GeV");
RooRealVar LPrime(
"LPrime",
"", .800, .600, 1.0,
"km");
RooGenericPdf PnmuTonePrime(
"PnmuTonePrime",
"P(#nu_{#mu} #rightarrow #nu_{e}", oscillationFormula, {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;
model.plotOn(Eframe);
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))};
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
Represents a constant real-valued object.
Container class to hold N-dimensional binned data.
const RooArgSet * get() const override
Get bin centre of current bin.
Container class to hold unbinned data.
Implementation of a probability density function that takes a RooArgList of servers and a C++ express...
static RooMsgService & instance()
Return reference to singleton instance.
void setStreamStatus(Int_t id, bool active)
(De)Activate stream with given unique ID
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.
Represents the product of a given set of RooAbsReal objects.
Implements the profile likelihood estimator for a given likelihood and set of parameters of interest.
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
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.
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg={})
RooCmdArg Scaling(bool flag)
RooCmdArg Extended(bool flag=true)
RooCmdArg Components(Args_t &&... argsOrArgSet)
RooCmdArg Binning(const RooAbsBinning &binning)
RooCmdArg LineColor(Color_t color)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
constexpr Double_t E()
Base of natural log: .