Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rs401d_FeldmanCousins.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook
4/// Neutrino Oscillation Example from Feldman & Cousins
5///
6/// This tutorial shows a more complex example using the FeldmanCousins utility
7/// to create a confidence interval for a toy neutrino oscillation experiment.
8/// The example attempts to faithfully reproduce the toy example described in Feldman & Cousins'
9/// original paper, Phys.Rev.D57:3873-3889,1998.
10///
11/// \macro_image
12/// \macro_output
13/// \macro_code
14///
15/// \author Kyle Cranmer
16
17#include "RooGlobalFunc.h"
26
27#include "RooDataSet.h"
28#include "RooDataHist.h"
29#include "RooRealVar.h"
30#include "RooConstVar.h"
31#include "RooAddition.h"
32#include "RooProduct.h"
33#include "RooProdPdf.h"
34#include "RooAddPdf.h"
35
36#include "TROOT.h"
37#include "RooPolynomial.h"
38#include "RooRandom.h"
39
40#include "RooProfileLL.h"
41
42#include "RooPlot.h"
43
44#include "TCanvas.h"
45#include "TH1F.h"
46#include "TH2F.h"
47#include "TTree.h"
48#include "TMarker.h"
49#include "TStopwatch.h"
50
51#include <iostream>
52
53// PDF class created for this macro
54#if !defined(__CINT__) || defined(__MAKECINT__)
55#include "../tutorials/roostats/NuMuToNuE_Oscillation.h"
56#include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx" // so that it can be executed directly
57#else
58#include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx+" // so that it can be executed directly
59#endif
60
61// use this order for safety on library loading
62using namespace RooFit;
63using namespace RooStats;
64
65void rs401d_FeldmanCousins(bool doFeldmanCousins = false, bool doMCMC = true)
66{
67
68 // to time the macro
69 TStopwatch t;
70 t.Start();
71
72 // Taken from Feldman & Cousins paper, Phys.Rev.D57:3873-3889,1998.
73 // e-Print: physics/9711021 (see page 13.)
74 //
75 // Quantum mechanics dictates that the probability of such a transformation is given by the formula
76 // $P (\nu\mu \rightarrow \nu e ) = sin^2 (2\theta) sin^2 (1.27 \Delta m^2 L /E )$
77 // where P is the probability for a $\nu\mu$ to transform into a $\nu e$ , L is the distance in km between
78 // the creation of the neutrino from meson decay and its interaction in the detector, E is the
79 // neutrino energy in GeV, and $\Delta m^2 = |m^2 - m^2 |$ in $(eV/c^2 )^2$ .
80 //
81 // To demonstrate how this works in practice, and how it compares to alternative approaches
82 // that have been used, we consider a toy model of a typical neutrino oscillation experiment.
83 // The toy model is defined by the following parameters: Mesons are assumed to decay to
84 // neutrinos uniformly in a region 600 m to 1000 m from the detector. The expected background
85 // from conventional $\nu e$ interactions and misidentified $\nu\mu$ interactions is assumed to be 100
86 // events in each of 5 energy bins which span the region from 10 to 60 GeV. We assume that
87 // the $\nu\mu$ flux is such that if $P (\nu\mu \rightarrow \nu e ) = 0.01$ averaged over any bin, then that bin
88 // would
89 // have an expected additional contribution of 100 events due to $\nu\mu \rightarrow \nu e$ oscillations.
90
91 // Make signal model model
92 RooRealVar E("E", "", 15, 10, 60, "GeV");
93 RooRealVar L("L", "", .800, .600, 1.0, "km"); // need these units in formula
94 RooRealVar deltaMSq("deltaMSq", "#Delta m^{2}", 40, 1, 300, "eV/c^{2}");
95 RooRealVar sinSq2theta("sinSq2theta", "sin^{2}(2#theta)", .006, .0, .02);
96 // RooRealVar deltaMSq("deltaMSq","#Delta m^{2}",40,20,70,"eV/c^{2}");
97 // RooRealVar sinSq2theta("sinSq2theta","sin^{2}(2#theta)", .006,.001,.01);
98 // PDF for oscillation only describes deltaMSq dependence, sinSq2theta goes into sigNorm
99 // 1) The code for this PDF was created by issuing these commands
100 // root [0] RooClassFactory x
101 // root [1] x.makePdf("NuMuToNuE_Oscillation","L,E,deltaMSq","","pow(sin(1.27*deltaMSq*L/E),2)")
102 NuMuToNuE_Oscillation PnmuTone("PnmuTone", "P(#nu_{#mu} #rightarrow #nu_{e}", L, E, deltaMSq);
103
104 // only E is observable, so create the signal model by integrating out L
105 RooAbsPdf *sigModel = PnmuTone.createProjection(L);
106
107 // create $ \int dE' dL' P(E',L' | \Delta m^2)$.
108 // Given RooFit will renormalize the PDF in the range of the observables,
109 // the average probability to oscillate in the experiment's acceptance
110 // needs to be incorporated into the extended term in the likelihood.
111 // Do this by creating a RooAbsReal representing the integral and divide by
112 // the area in the E-L plane.
113 // The integral should be over "primed" observables, so we need
114 // an independent copy of PnmuTone not to interfere with the original.
115
116 // Independent copy for Integral
117 RooRealVar EPrime("EPrime", "", 15, 10, 60, "GeV");
118 RooRealVar LPrime("LPrime", "", .800, .600, 1.0, "km"); // need these units in formula
119 NuMuToNuE_Oscillation PnmuTonePrime("PnmuTonePrime", "P(#nu_{#mu} #rightarrow #nu_{e}", LPrime, EPrime, deltaMSq);
120 RooAbsReal *intProbToOscInExp = PnmuTonePrime.createIntegral(RooArgSet(EPrime, LPrime));
121
122 // Getting the flux is a bit tricky. It is more clear to include a cross section term that is not
123 // explicitly referred to in the text, eg.
124 // number events in bin = flux * cross-section for nu_e interaction in E bin * average prob nu_mu osc. to nu_e in bin
125 // let maxEventsInBin = flux * cross-section for nu_e interaction in E bin
126 // maxEventsInBin * 1% chance per bin = 100 events / bin
127 // therefore maxEventsInBin = 10,000.
128 // for 5 bins, this means maxEventsTot = 50,000
129 RooConstVar maxEventsTot("maxEventsTot", "maximum number of sinal events", 50000);
130 RooConstVar inverseArea("inverseArea", "1/(#Delta E #Delta L)",
131 1. / (EPrime.getMax() - EPrime.getMin()) / (LPrime.getMax() - LPrime.getMin()));
132
133 // $sigNorm = maxEventsTot \cdot \int dE dL \frac{P_{oscillate\ in\ experiment}}{Area} \cdot {sin}^2(2\theta)$
134 RooProduct sigNorm("sigNorm", "", RooArgSet(maxEventsTot, *intProbToOscInExp, inverseArea, sinSq2theta));
135 // bkg = 5 bins * 100 events / bin
136 RooConstVar bkgNorm("bkgNorm", "normalization for background", 500);
137
138 // flat background (0th order polynomial, so no arguments for coefficients)
139 RooPolynomial bkgEShape("bkgEShape", "flat bkg shape", E);
140
141 // total model
142 RooAddPdf model("model", "", RooArgList(*sigModel, bkgEShape), RooArgList(sigNorm, bkgNorm));
143
144 // for debugging, check model tree
145 // model.printCompactTree();
146 // model.graphVizTree("model.dot");
147
148 // turn off some messages
152
153 // --------------------------------------
154 // n events in data to data, simply sum of sig+bkg
155 Int_t nEventsData = bkgNorm.getVal() + sigNorm.getVal();
156 cout << "generate toy data with nEvents = " << nEventsData << endl;
157 // adjust random seed to get a toy dataset similar to one in paper.
158 // Found by trial and error (3 trials, so not very "fine tuned")
160 // create a toy dataset
161 RooDataSet *data = model.generate(RooArgSet(E), nEventsData);
162
163 // --------------------------------------
164 // make some plots
165 TCanvas *dataCanvas = new TCanvas("dataCanvas");
166 dataCanvas->Divide(2, 2);
167
168 // plot the PDF
169 dataCanvas->cd(1);
170 TH1 *hh = PnmuTone.createHistogram("hh", E, Binning(40), YVar(L, Binning(40)), Scaling(kFALSE));
171 hh->SetLineColor(kBlue);
172 hh->SetTitle("True Signal Model");
173 hh->Draw("surf");
174
175 // plot the data with the best fit
176 dataCanvas->cd(2);
177 RooPlot *Eframe = E.frame();
178 data->plotOn(Eframe);
179 model.fitTo(*data, Extended());
180 model.plotOn(Eframe);
181 model.plotOn(Eframe, Components(*sigModel), LineColor(kRed));
182 model.plotOn(Eframe, Components(bkgEShape), LineColor(kGreen));
183 model.plotOn(Eframe);
184 Eframe->SetTitle("toy data with best fit model (and sig+bkg components)");
185 Eframe->Draw();
186
187 // plot the likelihood function
188 dataCanvas->cd(3);
189 std::unique_ptr<RooAbsReal> nll{model.createNLL(*data, Extended(true))};
190 RooProfileLL pll("pll", "", *nll, RooArgSet(deltaMSq, sinSq2theta));
191 // TH1* hhh = nll.createHistogram("hhh",sinSq2theta,Binning(40),YVar(deltaMSq,Binning(40))) ;
192 TH1 *hhh = pll.createHistogram("hhh", sinSq2theta, Binning(40), YVar(deltaMSq, Binning(40)), Scaling(kFALSE));
193 hhh->SetLineColor(kBlue);
194 hhh->SetTitle("Likelihood Function");
195 hhh->Draw("surf");
196
197 dataCanvas->Update();
198
199 // --------------------------------------------------------------
200 // show use of Feldman-Cousins utility in RooStats
201 // set the distribution creator, which encodes the test statistic
202 RooArgSet parameters(deltaMSq, sinSq2theta);
203 RooWorkspace *w = new RooWorkspace();
204
205 ModelConfig modelConfig;
206 modelConfig.SetWorkspace(*w);
207 modelConfig.SetPdf(model);
208 modelConfig.SetParametersOfInterest(parameters);
209
210 RooStats::FeldmanCousins fc(*data, modelConfig);
211 fc.SetTestSize(.1); // set size of test
212 fc.UseAdaptiveSampling(true);
213 fc.SetNBins(10); // number of points to test per parameter
214
215 // use the Feldman-Cousins tool
216 ConfInterval *interval = 0;
217 if (doFeldmanCousins)
218 interval = fc.GetInterval();
219
220 // ---------------------------------------------------------
221 // show use of ProfileLikeihoodCalculator utility in RooStats
223 plc.SetTestSize(.1);
224
225 ConfInterval *plcInterval = plc.GetInterval();
226
227 // --------------------------------------------
228 // show use of MCMCCalculator utility in RooStats
229 MCMCInterval *mcInt = NULL;
230
231 if (doMCMC) {
232 // turn some messages back on
235
236 TStopwatch mcmcWatch;
237 mcmcWatch.Start();
238
239 RooArgList axisList(deltaMSq, sinSq2theta);
240 MCMCCalculator mc(*data, modelConfig);
241 mc.SetNumIters(5000);
242 mc.SetNumBurnInSteps(100);
243 mc.SetUseKeys(true);
244 mc.SetTestSize(.1);
245 mc.SetAxes(axisList); // set which is x and y axis in posterior histogram
246 // mc.SetNumBins(50);
247 mcInt = (MCMCInterval *)mc.GetInterval();
248
249 mcmcWatch.Stop();
250 mcmcWatch.Print();
251 }
252 // -------------------------------
253 // make plot of resulting interval
254
255 dataCanvas->cd(4);
256
257 // first plot a small dot for every point tested
258 if (doFeldmanCousins) {
259 RooDataHist *parameterScan = (RooDataHist *)fc.GetPointsToScan();
260 TH2F *hist = parameterScan->createHistogram(deltaMSq,sinSq2theta, 30, 30);
261 // hist->Draw();
262 TH2F *forContour = (TH2F *)hist->Clone();
263
264 // now loop through the points and put a marker if it's in the interval
265 RooArgSet *tmpPoint;
266 // loop over points to test
267 for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
268 // get a parameter point from the list of points to test.
269 tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");
270
271 if (interval) {
272 if (interval->IsInInterval(*tmpPoint)) {
273 forContour->SetBinContent(
274 hist->FindBin(tmpPoint->getRealValue("sinSq2theta"), tmpPoint->getRealValue("deltaMSq")), 1);
275 } else {
276 forContour->SetBinContent(
277 hist->FindBin(tmpPoint->getRealValue("sinSq2theta"), tmpPoint->getRealValue("deltaMSq")), 0);
278 }
279 }
280
281 delete tmpPoint;
282 }
283
284 if (interval) {
285 Double_t level = 0.5;
286 forContour->SetContour(1, &level);
287 forContour->SetLineWidth(2);
288 forContour->SetLineColor(kRed);
289 forContour->Draw("cont2,same");
290 }
291 }
292
293 MCMCIntervalPlot *mcPlot = NULL;
294 if (mcInt) {
295 cout << "MCMC actual confidence level: " << mcInt->GetActualConfidenceLevel() << endl;
296 mcPlot = new MCMCIntervalPlot(*mcInt);
297 mcPlot->SetLineColor(kMagenta);
298 mcPlot->Draw();
299 }
300 dataCanvas->Update();
301
302 LikelihoodIntervalPlot plotInt((LikelihoodInterval *)plcInterval);
303 plotInt.SetTitle("90% Confidence Intervals");
304 if (mcInt)
305 plotInt.Draw("same");
306 else
307 plotInt.Draw();
308 dataCanvas->Update();
309
310 /// print timing info
311 t.Stop();
312 t.Print();
313}
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.
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)
RooArgList L(Args_t &&... args)
Definition RooArgList.h:156
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
constexpr Double_t E()
Base of natural log: .
Definition TMath.h:93