ROOT   6.21/01 Reference Guide
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"
18 #include "RooStats/ConfInterval.h"
25 #include "RooStats/MCMCInterval.h"
26
27 #include "RooDataSet.h"
28 #include "RooDataHist.h"
29 #include "RooRealVar.h"
30 #include "RooConstVar.h"
32 #include "RooProduct.h"
33 #include "RooProdPdf.h"
35
36 #include "TROOT.h"
37 #include "RooPolynomial.h"
38 #include "RooRandom.h"
39
40 #include "RooNLLVar.h"
41 #include "RooProfileLL.h"
42
43 #include "RooPlot.h"
44
45 #include "TCanvas.h"
46 #include "TH1F.h"
47 #include "TH2F.h"
48 #include "TTree.h"
49 #include "TMarker.h"
50 #include "TStopwatch.h"
51
52 #include <iostream>
53
54 // PDF class created for this macro
55 #if !defined(__CINT__) || defined(__MAKECINT__)
56 #include "../tutorials/roostats/NuMuToNuE_Oscillation.h"
57 #include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx" // so that it can be executed directly
58 #else
59 #include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx+" // so that it can be executed directly
60 #endif
61
63 using namespace RooFit;
64 using namespace RooStats;
65
66 void rs401d_FeldmanCousins(bool doFeldmanCousins = false, bool doMCMC = true)
67 {
68
69  // to time the macro
70  TStopwatch t;
71  t.Start();
72
73  // Taken from Feldman & Cousins paper, Phys.Rev.D57:3873-3889,1998.
74  // e-Print: physics/9711021 (see page 13.)
75  //
76  // Quantum mechanics dictates that the probability of such a transformation is given by the formula
77  // $P (\nu\mu \rightarrow \nu e ) = sin^2 (2\theta) sin^2 (1.27 \Delta m^2 L /E )$
78  // where P is the probability for a $\nu\mu$ to transform into a $\nu e$ , L is the distance in km between
79  // the creation of the neutrino from meson decay and its interaction in the detector, E is the
80  // neutrino energy in GeV, and $\Delta m^2 = |m^2 - m^2 |$ in $(eV/c^2 )^2$ .
81  //
82  // To demonstrate how this works in practice, and how it compares to alternative approaches
83  // that have been used, we consider a toy model of a typical neutrino oscillation experiment.
84  // The toy model is defined by the following parameters: Mesons are assumed to decay to
85  // neutrinos uniformly in a region 600 m to 1000 m from the detector. The expected background
86  // from conventional $\nu e$ interactions and misidentified $\nu\mu$ interactions is assumed to be 100
87  // events in each of 5 energy bins which span the region from 10 to 60 GeV. We assume that
88  // the $\nu\mu$ flux is such that if $P (\nu\mu \rightarrow \nu e ) = 0.01$ averaged over any bin, then that bin
89  // would
90  // have an expected additional contribution of 100 events due to $\nu\mu \rightarrow \nu e$ oscillations.
91
92  // Make signal model model
93  RooRealVar E("E", "", 15, 10, 60, "GeV");
94  RooRealVar L("L", "", .800, .600, 1.0, "km"); // need these units in formula
95  RooRealVar deltaMSq("deltaMSq", "#Delta m^{2}", 40, 1, 300, "eV/c^{2}");
96  RooRealVar sinSq2theta("sinSq2theta", "sin^{2}(2#theta)", .006, .0, .02);
97  // RooRealVar deltaMSq("deltaMSq","#Delta m^{2}",40,20,70,"eV/c^{2}");
98  // RooRealVar sinSq2theta("sinSq2theta","sin^{2}(2#theta)", .006,.001,.01);
99  // PDF for oscillation only describes deltaMSq dependence, sinSq2theta goes into sigNorm
100  // 1) The code for this PDF was created by issuing these commands
101  // root [0] RooClassFactory x
102  // root [1] x.makePdf("NuMuToNuE_Oscillation","L,E,deltaMSq","","pow(sin(1.27*deltaMSq*L/E),2)")
103  NuMuToNuE_Oscillation PnmuTone("PnmuTone", "P(#nu_{#mu} #rightarrow #nu_{e}", L, E, deltaMSq);
104
105  // only E is observable, so create the signal model by integrating out L
106  RooAbsPdf *sigModel = PnmuTone.createProjection(L);
107
108  // create $\int dE' dL' P(E',L' | \Delta m^2)$.
109  // Given RooFit will renormalize the PDF in the range of the observables,
110  // the average probability to oscillate in the experiment's acceptance
111  // needs to be incorporated into the extended term in the likelihood.
112  // Do this by creating a RooAbsReal representing the integral and divide by
113  // the area in the E-L plane.
114  // The integral should be over "primed" observables, so we need
115  // an independent copy of PnmuTone not to interfere with the original.
116
117  // Independent copy for Integral
118  RooRealVar EPrime("EPrime", "", 15, 10, 60, "GeV");
119  RooRealVar LPrime("LPrime", "", .800, .600, 1.0, "km"); // need these units in formula
120  NuMuToNuE_Oscillation PnmuTonePrime("PnmuTonePrime", "P(#nu_{#mu} #rightarrow #nu_{e}", LPrime, EPrime, deltaMSq);
121  RooAbsReal *intProbToOscInExp = PnmuTonePrime.createIntegral(RooArgSet(EPrime, LPrime));
122
123  // Getting the flux is a bit tricky. It is more clear to include a cross section term that is not
124  // explicitly referred to in the text, eg.
125  // number events in bin = flux * cross-section for nu_e interaction in E bin * average prob nu_mu osc. to nu_e in bin
126  // let maxEventsInBin = flux * cross-section for nu_e interaction in E bin
127  // maxEventsInBin * 1% chance per bin = 100 events / bin
128  // therefore maxEventsInBin = 10,000.
129  // for 5 bins, this means maxEventsTot = 50,000
130  RooConstVar maxEventsTot("maxEventsTot", "maximum number of sinal events", 50000);
131  RooConstVar inverseArea("inverseArea", "1/(#Delta E #Delta L)",
132  1. / (EPrime.getMax() - EPrime.getMin()) / (LPrime.getMax() - LPrime.getMin()));
133
134  // $sigNorm = maxEventsTot \cdot \int dE dL \frac{P_{oscillate\ in\ experiment}}{Area} \cdot {sin}^2(2\theta)$
135  RooProduct sigNorm("sigNorm", "", RooArgSet(maxEventsTot, *intProbToOscInExp, inverseArea, sinSq2theta));
136  // bkg = 5 bins * 100 events / bin
137  RooConstVar bkgNorm("bkgNorm", "normalization for background", 500);
138
139  // flat background (0th order polynomial, so no arguments for coefficients)
140  RooPolynomial bkgEShape("bkgEShape", "flat bkg shape", E);
141
142  // total model
143  RooAddPdf model("model", "", RooArgList(*sigModel, bkgEShape), RooArgList(sigNorm, bkgNorm));
144
145  // for debugging, check model tree
146  // model.printCompactTree();
147  // model.graphVizTree("model.dot");
148
149  // turn off some messages
153
154  // --------------------------------------
155  // n events in data to data, simply sum of sig+bkg
156  Int_t nEventsData = bkgNorm.getVal() + sigNorm.getVal();
157  cout << "generate toy data with nEvents = " << nEventsData << endl;
158  // adjust random seed to get a toy dataset similar to one in paper.
159  // Found by trial and error (3 trials, so not very "fine tuned")
161  // create a toy dataset
162  RooDataSet *data = model.generate(RooArgSet(E), nEventsData);
163
164  // --------------------------------------
165  // make some plots
166  TCanvas *dataCanvas = new TCanvas("dataCanvas");
167  dataCanvas->Divide(2, 2);
168
169  // plot the PDF
170  dataCanvas->cd(1);
171  TH1 *hh = PnmuTone.createHistogram("hh", E, Binning(40), YVar(L, Binning(40)), Scaling(kFALSE));
172  hh->SetLineColor(kBlue);
173  hh->SetTitle("True Signal Model");
174  hh->Draw("surf");
175
176  // plot the data with the best fit
177  dataCanvas->cd(2);
178  RooPlot *Eframe = E.frame();
179  data->plotOn(Eframe);
180  model.fitTo(*data, Extended());
181  model.plotOn(Eframe);
182  model.plotOn(Eframe, Components(*sigModel), LineColor(kRed));
183  model.plotOn(Eframe, Components(bkgEShape), LineColor(kGreen));
184  model.plotOn(Eframe);
185  Eframe->SetTitle("toy data with best fit model (and sig+bkg components)");
186  Eframe->Draw();
187
188  // plot the likelihood function
189  dataCanvas->cd(3);
190  RooNLLVar nll("nll", "nll", model, *data, Extended());
191  RooProfileLL pll("pll", "", nll, RooArgSet(deltaMSq, sinSq2theta));
192  // TH1* hhh = nll.createHistogram("hhh",sinSq2theta,Binning(40),YVar(deltaMSq,Binning(40))) ;
193  TH1 *hhh = pll.createHistogram("hhh", sinSq2theta, Binning(40), YVar(deltaMSq, Binning(40)), Scaling(kFALSE));
194  hhh->SetLineColor(kBlue);
195  hhh->SetTitle("Likelihood Function");
196  hhh->Draw("surf");
197
198  dataCanvas->Update();
199
200  // --------------------------------------------------------------
201  // show use of Feldman-Cousins utility in RooStats
202  // set the distribution creator, which encodes the test statistic
203  RooArgSet parameters(deltaMSq, sinSq2theta);
204  RooWorkspace *w = new RooWorkspace();
205
206  ModelConfig modelConfig;
207  modelConfig.SetWorkspace(*w);
208  modelConfig.SetPdf(model);
209  modelConfig.SetParametersOfInterest(parameters);
210
211  RooStats::FeldmanCousins fc(*data, modelConfig);
212  fc.SetTestSize(.1); // set size of test
214  fc.SetNBins(10); // number of points to test per parameter
215
216  // use the Feldman-Cousins tool
217  ConfInterval *interval = 0;
218  if (doFeldmanCousins)
219  interval = fc.GetInterval();
220
221  // ---------------------------------------------------------
222  // show use of ProfileLikeihoodCalculator utility in RooStats
223  RooStats::ProfileLikelihoodCalculator plc(*data, modelConfig);
224  plc.SetTestSize(.1);
225
226  ConfInterval *plcInterval = plc.GetInterval();
227
228  // --------------------------------------------
229  // show use of MCMCCalculator utility in RooStats
230  MCMCInterval *mcInt = NULL;
231
232  if (doMCMC) {
233  // turn some messages back on
236
237  TStopwatch mcmcWatch;
238  mcmcWatch.Start();
239
240  RooArgList axisList(deltaMSq, sinSq2theta);
241  MCMCCalculator mc(*data, modelConfig);
242  mc.SetNumIters(5000);
243  mc.SetNumBurnInSteps(100);
244  mc.SetUseKeys(true);
245  mc.SetTestSize(.1);
246  mc.SetAxes(axisList); // set which is x and y axis in posterior histogram
247  // mc.SetNumBins(50);
248  mcInt = (MCMCInterval *)mc.GetInterval();
249
250  mcmcWatch.Stop();
251  mcmcWatch.Print();
252  }
253  // -------------------------------
254  // make plot of resulting interval
255
256  dataCanvas->cd(4);
257
258  // first plot a small dot for every point tested
259  if (doFeldmanCousins) {
260  RooDataHist *parameterScan = (RooDataHist *)fc.GetPointsToScan();
261  TH2F *hist = (TH2F *)parameterScan->createHistogram("sinSq2theta:deltaMSq", 30, 30);
262  // hist->Draw();
263  TH2F *forContour = (TH2F *)hist->Clone();
264
265  // now loop through the points and put a marker if it's in the interval
266  RooArgSet *tmpPoint;
267  // loop over points to test
268  for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
269  // get a parameter point from the list of points to test.
270  tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");
271
272  if (interval) {
273  if (interval->IsInInterval(*tmpPoint)) {
274  forContour->SetBinContent(
275  hist->FindBin(tmpPoint->getRealValue("sinSq2theta"), tmpPoint->getRealValue("deltaMSq")), 1);
276  } else {
277  forContour->SetBinContent(
278  hist->FindBin(tmpPoint->getRealValue("sinSq2theta"), tmpPoint->getRealValue("deltaMSq")), 0);
279  }
280  }
281
282  delete tmpPoint;
283  }
284
285  if (interval) {
286  Double_t level = 0.5;
287  forContour->SetContour(1, &level);
288  forContour->SetLineWidth(2);
289  forContour->SetLineColor(kRed);
290  forContour->Draw("cont2,same");
291  }
292  }
293
294  MCMCIntervalPlot *mcPlot = NULL;
295  if (mcInt) {
296  cout << "MCMC actual confidence level: " << mcInt->GetActualConfidenceLevel() << endl;
297  mcPlot = new MCMCIntervalPlot(*mcInt);
298  mcPlot->SetLineColor(kMagenta);
299  mcPlot->Draw();
300  }
301  dataCanvas->Update();
302
303  LikelihoodIntervalPlot plotInt((LikelihoodInterval *)plcInterval);
304  plotInt.SetTitle("90% Confidence Intervals");
305  if (mcInt)
306  plotInt.Draw("same");
307  else
308  plotInt.Draw();
309  dataCanvas->Update();
310
311  /// print timing info
312  t.Stop();
313  t.Print();
314 }
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
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:3579
RooAddPdf is an efficient implementation of a sum of PDFs of the form.
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
virtual const RooArgSet * get() const
Definition: RooDataHist.h:79
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
Calls RooPlot* plotOn(RooPlot* frame, const RooLinkedList& cmdList) const ;.
Definition: RooAbsData.cxx:549
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
This class provides simple and straightforward utilities to plot a MCMCInterval object.
Class RooProfileLL implements the profile likelihood estimator for a given likelihood and set of para...
Definition: RooProfileLL.h:26
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
RooCmdArg LineColor(Color_t color)
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual void SetWorkspace(RooWorkspace &ws)
Definition: ModelConfig.h:66
Definition: Rtypes.h:64
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
virtual void SetContour(Int_t nlevels, const Double_t *levels=0)
Set the number and values of contour levels.
Definition: TH1.cxx:7913
Definition: TCanvas.cxx:696
This class provides simple and straightforward utilities to plot a LikelihoodInterval object...
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
virtual Bool_t IsInInterval(const RooArgSet &) const =0
check if given point is in the interval
int Int_t
Definition: RtypesCore.h:41
void Draw(const Option_t *options=NULL)
void SetTitle(const char *name)
Set the title of the RooPlot to &#39;title&#39;.
Definition: RooPlot.cxx:1258
Definition: Rtypes.h:64
static RooMsgService & instance()
Return reference to singleton instance.
RooCmdArg Extended(Bool_t flag=kTRUE)
void setStreamStatus(Int_t id, Bool_t active)
(De)Activate stream with given unique ID
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:597
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the the workspace if not already there.
Definition: ModelConfig.h:81
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:54
RooConstVar represent a constant real-valued object.
Definition: RooConstVar.h:25
virtual RooAbsPdf * createProjection(const RooArgSet &iset)
Return a p.d.f that represent a projection of this p.d.f integrated over given observables.
Definition: RooAbsPdf.cxx:3332
virtual TObject * clone(const char *newname) const
Definition: RooArgSet.h:84
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
RooAbsReal * createIntegral(const RooArgSet &iset, 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 an object that represents the integral of the function over one or more observables listed in ...
Definition: RooAbsReal.cxx:562
Class RooNLLVar implements a -log(likelihood) calculation from a dataset and a PDF.
Definition: RooNLLVar.h:26
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2981
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
void SetLineColor(Color_t color)
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
A RooProduct represents the product of a given set of RooAbsReal objects.
Definition: RooProduct.h:32
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:97
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:44
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
const Bool_t kFALSE
Definition: RtypesCore.h:88
virtual Double_t GetActualConfidenceLevel()
virtual Double_t GetKeysPdfCutoff() { return fKeysCutoff; }
The Canvas class.
Definition: TCanvas.h:31
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
Namespace for the RooStats classes.
Definition: Asimov.h:20
virtual Int_t numEntries() const
Return the number of bins.
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
static constexpr double L
RooCmdArg Components(const RooArgSet &compSet)
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
The TH1 histogram class.
Definition: TH1.h:56
RooCmdArg Scaling(Bool_t flag)
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Definition: RooArgSet.cxx:472
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2664
RooPolynomial implements a polynomial p.d.f of the form By default, the coefficient is chosen to be...
Definition: RooPolynomial.h:28
virtual void SetParametersOfInterest(const RooArgSet &set)
Specify parameters of interest.
Definition: ModelConfig.h:100
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
Definition: Rtypes.h:64
virtual void SetTitle(const char *title)
Definition: TH1.cxx:6316
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2441
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
virtual void Update()
Definition: TCanvas.cxx:2339
TH1 * createHistogram(const char *name, const RooAbsRealLValue &xvar, 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
Calls createHistogram(const char *name, const RooAbsRealLValue& xvar, const RooLinkedList& argList) c...
Definition: RooAbsData.cxx:631
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
const Bool_t kTRUE
Definition: RtypesCore.h:87
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712
RooCmdArg Binning(const RooAbsBinning &binning)
Stopwatch class.
Definition: TStopwatch.h:28