Logo ROOT  
Reference Guide
No Matches
rs401d_FeldmanCousins.C File Reference

Detailed Description

View in nbviewer Open in SWAN
Neutrino Oscillation Example from Feldman & Cousins

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:37, CP time 157.810
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:09, CP time 189.140
#include "RooGlobalFunc.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooRealVar.h"
#include "RooConstVar.h"
#include "RooAddition.h"
#include "RooProduct.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "TROOT.h"
#include "RooPolynomial.h"
#include "RooRandom.h"
#include "RooProfileLL.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TTree.h"
#include "TMarker.h"
#include "TStopwatch.h"
#include <iostream>
// PDF class created for this macro
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "../tutorials/roostats/NuMuToNuE_Oscillation.h"
#include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx" // so that it can be executed directly
#include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx+" // so that it can be executed directly
// use this order for safety on library loading
using namespace RooFit;
using namespace RooStats;
void rs401d_FeldmanCousins(bool doFeldmanCousins = false, bool doMCMC = true)
// to time the macro
// Taken from Feldman & Cousins paper, Phys.Rev.D57:3873-3889,1998.
// e-Print: physics/9711021 (see page 13.)
// Quantum mechanics dictates that the probability of such a transformation is given by the formula
// $P (\nu\mu \rightarrow \nu e ) = sin^2 (2\theta) sin^2 (1.27 \Delta m^2 L /E )$
// where P is the probability for a $\nu\mu$ to transform into a $\nu e$ , L is the distance in km between
// the creation of the neutrino from meson decay and its interaction in the detector, E is the
// neutrino energy in GeV, and $\Delta m^2 = |m^2 - m^2 |$ in $(eV/c^2 )^2$ .
// To demonstrate how this works in practice, and how it compares to alternative approaches
// that have been used, we consider a toy model of a typical neutrino oscillation experiment.
// The toy model is defined by the following parameters: Mesons are assumed to decay to
// neutrinos uniformly in a region 600 m to 1000 m from the detector. The expected background
// from conventional $\nu e$ interactions and misidentified $\nu\mu$ interactions is assumed to be 100
// events in each of 5 energy bins which span the region from 10 to 60 GeV. We assume that
// the $\nu\mu$ flux is such that if $P (\nu\mu \rightarrow \nu e ) = 0.01$ averaged over any bin, then that bin
// would
// have an expected additional contribution of 100 events due to $\nu\mu \rightarrow \nu e$ oscillations.
// Make signal model model
RooRealVar E("E", "", 15, 10, 60, "GeV");
RooRealVar L("L", "", .800, .600, 1.0, "km"); // need these units in formula
RooRealVar deltaMSq("deltaMSq", "#Delta m^{2}", 40, 1, 300, "eV/c^{2}");
RooRealVar sinSq2theta("sinSq2theta", "sin^{2}(2#theta)", .006, .0, .02);
// RooRealVar deltaMSq("deltaMSq","#Delta m^{2}",40,20,70,"eV/c^{2}");
// RooRealVar sinSq2theta("sinSq2theta","sin^{2}(2#theta)", .006,.001,.01);
// PDF for oscillation only describes deltaMSq dependence, sinSq2theta goes into sigNorm
// 1) The code for this PDF was created by issuing these commands
// root [0] RooClassFactory x
// root [1] x.makePdf("NuMuToNuE_Oscillation","L,E,deltaMSq","","pow(sin(1.27*deltaMSq*L/E),2)")
NuMuToNuE_Oscillation PnmuTone("PnmuTone", "P(#nu_{#mu} #rightarrow #nu_{e}", L, E, deltaMSq);
// only E is observable, so create the signal model by integrating out L
RooAbsPdf *sigModel = PnmuTone.createProjection(L);
// create $ \int dE' dL' P(E',L' | \Delta m^2)$.
// Given RooFit will renormalize the PDF in the range of the observables,
// the average probability to oscillate in the experiment's acceptance
// needs to be incorporated into the extended term in the likelihood.
// Do this by creating a RooAbsReal representing the integral and divide by
// the area in the E-L plane.
// The integral should be over "primed" observables, so we need
// an independent copy of PnmuTone not to interfere with the original.
// Independent copy for Integral
RooRealVar EPrime("EPrime", "", 15, 10, 60, "GeV");
RooRealVar LPrime("LPrime", "", .800, .600, 1.0, "km"); // need these units in formula
NuMuToNuE_Oscillation PnmuTonePrime("PnmuTonePrime", "P(#nu_{#mu} #rightarrow #nu_{e}", LPrime, EPrime, deltaMSq);
RooAbsReal *intProbToOscInExp = PnmuTonePrime.createIntegral(RooArgSet(EPrime, LPrime));
// Getting the flux is a bit tricky. It is more clear to include a cross section term that is not
// explicitly referred to in the text, eg.
// number events in bin = flux * cross-section for nu_e interaction in E bin * average prob nu_mu osc. to nu_e in bin
// let maxEventsInBin = flux * cross-section for nu_e interaction in E bin
// maxEventsInBin * 1% chance per bin = 100 events / bin
// therefore maxEventsInBin = 10,000.
// for 5 bins, this means maxEventsTot = 50,000
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()));
// $sigNorm = maxEventsTot \cdot \int dE dL \frac{P_{oscillate\ in\ experiment}}{Area} \cdot {sin}^2(2\theta)$
RooProduct sigNorm("sigNorm", "", RooArgSet(maxEventsTot, *intProbToOscInExp, inverseArea, sinSq2theta));
// bkg = 5 bins * 100 events / bin
RooConstVar bkgNorm("bkgNorm", "normalization for background", 500);
// flat background (0th order polynomial, so no arguments for coefficients)
RooPolynomial bkgEShape("bkgEShape", "flat bkg shape", E);
// total model
RooAddPdf model("model", "", RooArgList(*sigModel, bkgEShape), RooArgList(sigNorm, bkgNorm));
// for debugging, check model tree
// model.printCompactTree();
// model.graphVizTree("model.dot");
// turn off some messages
// --------------------------------------
// n events in data to data, simply sum of sig+bkg
Int_t nEventsData = bkgNorm.getVal() + sigNorm.getVal();
cout << "generate toy data with nEvents = " << nEventsData << endl;
// adjust random seed to get a toy dataset similar to one in paper.
// Found by trial and error (3 trials, so not very "fine tuned")
// create a toy dataset
RooDataSet *data = model.generate(RooArgSet(E), nEventsData);
// --------------------------------------
// make some plots
TCanvas *dataCanvas = new TCanvas("dataCanvas");
dataCanvas->Divide(2, 2);
// plot the PDF
TH1 *hh = PnmuTone.createHistogram("hh", E, Binning(40), YVar(L, Binning(40)), Scaling(kFALSE));
hh->SetTitle("True Signal Model");
// plot the data with the best fit
RooPlot *Eframe = E.frame();
model.fitTo(*data, Extended());
model.plotOn(Eframe, Components(*sigModel), LineColor(kRed));
model.plotOn(Eframe, Components(bkgEShape), LineColor(kGreen));
Eframe->SetTitle("toy data with best fit model (and sig+bkg components)");
// plot the likelihood function
std::unique_ptr<RooAbsReal> nll{model.createNLL(*data, Extended(true))};
RooProfileLL pll("pll", "", *nll, RooArgSet(deltaMSq, sinSq2theta));
// TH1* hhh = nll.createHistogram("hhh",sinSq2theta,Binning(40),YVar(deltaMSq,Binning(40))) ;
TH1 *hhh = pll.createHistogram("hhh", sinSq2theta, Binning(40), YVar(deltaMSq, Binning(40)), Scaling(kFALSE));
hhh->SetTitle("Likelihood Function");
// --------------------------------------------------------------
// show use of Feldman-Cousins utility in RooStats
// set the distribution creator, which encodes the test statistic
RooArgSet parameters(deltaMSq, sinSq2theta);
ModelConfig modelConfig;
RooStats::FeldmanCousins fc(*data, modelConfig);
fc.SetTestSize(.1); // set size of test
fc.SetNBins(10); // number of points to test per parameter
// use the Feldman-Cousins tool
ConfInterval *interval = 0;
if (doFeldmanCousins)
interval = fc.GetInterval();
// ---------------------------------------------------------
// show use of ProfileLikeihoodCalculator utility in RooStats
ConfInterval *plcInterval = plc.GetInterval();
// --------------------------------------------
// show use of MCMCCalculator utility in RooStats
MCMCInterval *mcInt = NULL;
if (doMCMC) {
// turn some messages back on
TStopwatch mcmcWatch;
RooArgList axisList(deltaMSq, sinSq2theta);
MCMCCalculator mc(*data, modelConfig);
mc.SetAxes(axisList); // set which is x and y axis in posterior histogram
// mc.SetNumBins(50);
mcInt = (MCMCInterval *)mc.GetInterval();
// -------------------------------
// make plot of resulting interval
// first plot a small dot for every point tested
if (doFeldmanCousins) {
RooDataHist *parameterScan = (RooDataHist *)fc.GetPointsToScan();
TH2F *hist = parameterScan->createHistogram(deltaMSq,sinSq2theta, 30, 30);
// hist->Draw();
TH2F *forContour = (TH2F *)hist->Clone();
// now loop through the points and put a marker if it's in the interval
RooArgSet *tmpPoint;
// loop over points to test
for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
// get a parameter point from the list of points to test.
tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");
if (interval) {
if (interval->IsInInterval(*tmpPoint)) {
hist->FindBin(tmpPoint->getRealValue("sinSq2theta"), tmpPoint->getRealValue("deltaMSq")), 1);
} else {
hist->FindBin(tmpPoint->getRealValue("sinSq2theta"), tmpPoint->getRealValue("deltaMSq")), 0);
delete tmpPoint;
if (interval) {
Double_t level = 0.5;
forContour->SetContour(1, &level);
MCMCIntervalPlot *mcPlot = NULL;
if (mcInt) {
cout << "MCMC actual confidence level: " << mcInt->GetActualConfidenceLevel() << endl;
mcPlot = new MCMCIntervalPlot(*mcInt);
LikelihoodIntervalPlot plotInt((LikelihoodInterval *)plcInterval);
plotInt.SetTitle("90% Confidence Intervals");
if (mcInt)
/// print timing info
int Int_t
Definition RtypesCore.h:45
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
@ kRed
Definition Rtypes.h:66
@ kGreen
Definition Rtypes.h:66
@ kMagenta
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
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.
Definition RooAbsPdf.h:40
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...
Definition RooAbsReal.h:59
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.
Definition RooAddPdf.h:33
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
TObject * clone(const char *newname) const override
Definition RooArgSet.h:148
RooConstVar represent a constant real-valued object.
Definition RooConstVar.h:23
The RooDataHist is a container class to hold N-dimensional binned data.
Definition RooDataHist.h:39
const RooArgSet * get() const override
Get bin centre of current bin.
Definition RooDataHist.h:76
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:57
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.
Definition RooPlot.h:43
void SetTitle(const char *name) override
Set the title of the RooPlot to 'title'.
Definition RooPlot.cxx:1258
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:652
RooPolynomial implements a polynomial p.d.f of the form.
A RooProduct represents the product of a given set of RooAbsReal objects.
Definition RooProduct.h:29
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.
Definition RooRandom.cxx:51
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:37
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...
Definition ModelConfig.h:35
virtual void SetWorkspace(RooWorkspace &ws)
Definition ModelConfig.h:70
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.
Definition ModelConfig.h:87
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.
Definition TAttLine.h:43
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
The Canvas class.
Definition TCanvas.h:23
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
Definition TCanvas.cxx:716
void Update() override
Update canvas pad buffers.
Definition TCanvas.cxx:2475
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
void SetTitle(const char *title) override
Change/set the title.
Definition TH1.cxx:6707
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3067
virtual void SetContour(Int_t nlevels, const Double_t *levels=nullptr)
Set the number and values of contour levels.
Definition TH1.cxx:8400
virtual Int_t FindBin(Double_t x, Double_t y=0, Double_t z=0)
Return Global bin number corresponding to x,y,z.
Definition TH1.cxx:3675
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition TH1.cxx:2734
2-D histogram with a float per channel (see TH1 documentation)}
Definition TH2.h:258
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
Definition TH2.cxx:2554
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.
Definition TPad.cxx:1153
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition TRandom.cxx:608
Stopwatch class.
Definition TStopwatch.h:28
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...
Definition JSONIO.h:26
Namespace for the RooStats classes.
Definition Asimov.h:19
Kyle Cranmer

Definition in file rs401d_FeldmanCousins.C.