Logo ROOT  
Reference Guide
MCMCIntervalPlot.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Authors: Kevin Belasco 17/06/2009
3 // Authors: Kyle Cranmer 17/06/2009
4 /*************************************************************************
5  * Project: RooStats *
6  * Package: RooFit/RooStats *
7  *************************************************************************
8  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
9  * All rights reserved. *
10  * *
11  * For the licensing terms see $ROOTSYS/LICENSE. *
12  * For the list of contributors see $ROOTSYS/README/CREDITS. *
13  *************************************************************************/
14 
15 /** \class RooStats::MCMCIntervalPlot
16  \ingroup Roostats
17 
18 This class provides simple and straightforward utilities to plot a MCMCInterval
19 object. Basic use only requires a few lines once you have an MCMCInterval*:
20 
21 ~~~ {.cpp}
22  MCMCIntervalPlot plot(*interval);
23  plot.Draw();
24 ~~~
25 
26 The standard Draw() function will currently draw the confidence interval
27 range with bars if 1-D and a contour if 2-D. The MCMC posterior will also be
28 plotted for the 1-D case.
29 
30 */
31 
33 #include "RooStats/MCMCInterval.h"
34 #include "RooStats/MarkovChain.h"
35 
36 #include "TLine.h"
37 #include "TList.h"
38 #include "TGraph.h"
39 #include "RooRealVar.h"
40 #include "RooPlot.h"
41 #include "TH2.h"
42 #include "TH1F.h"
43 #include "RooArgList.h"
44 #include "TAxis.h"
45 #include "RooGlobalFunc.h"
46 
47 #include <iostream>
48 
49 // Extra draw commands
50 //static const char* POSTERIOR_HIST = "posterior_hist";
51 //static const char* POSTERIOR_KEYS_PDF = "posterior_keys_pdf";
52 //static const char* POSTERIOR_KEYS_PRODUCT = "posterior_keys_product";
53 //static const char* HIST_INTERVAL = "hist_interval";
54 //static const char* KEYS_PDF_INTERVAL = "keys_pdf_interval";
55 //static const char* TAIL_FRACTION_INTERVAL = "tail_fraction_interval";
56 //static const char* OPTION_SEP = ":";
57 
59 
60 using namespace std;
61 using namespace RooStats;
62 
63 ////////////////////////////////////////////////////////////////////////////////
64 
65 MCMCIntervalPlot::MCMCIntervalPlot()
66 {
67  fInterval = NULL;
68  fParameters = NULL;
69  fPosteriorHist = NULL;
70  fPosteriorKeysPdf = NULL;
71  fPosteriorKeysProduct = NULL;
72  fDimension = 0;
73  fLineColor = kBlack;
74  fShadeColor = kGray;
75  fLineWidth = 1;
76  //fContourColor = kBlack;
77  fShowBurnIn = kTRUE;
78  fWalk = NULL;
79  fBurnIn = NULL;
80  fFirst = NULL;
81  fParamGraph = NULL;
82  fNLLGraph = NULL;
83  fNLLHist = NULL;
84  fWeightHist = NULL;
85  fPosteriorHistHistCopy = NULL;
86  fPosteriorHistTFCopy = NULL;
87 }
88 
89 ////////////////////////////////////////////////////////////////////////////////
90 
91 MCMCIntervalPlot::MCMCIntervalPlot(MCMCInterval& interval)
92 {
93  SetMCMCInterval(interval);
94  fPosteriorHist = NULL;
95  fPosteriorKeysPdf = NULL;
96  fPosteriorKeysProduct = NULL;
97  fLineColor = kBlack;
98  fShadeColor = kGray;
99  fLineWidth = 1;
100  //fContourColor = kBlack;
101  fShowBurnIn = kTRUE;
102  fWalk = NULL;
103  fBurnIn = NULL;
104  fFirst = NULL;
105  fParamGraph = NULL;
106  fNLLGraph = NULL;
107  fNLLHist = NULL;
108  fWeightHist = NULL;
109  fPosteriorHistHistCopy = NULL;
110  fPosteriorHistTFCopy = NULL;
111 }
112 
113 ////////////////////////////////////////////////////////////////////////////////
114 
115 MCMCIntervalPlot::~MCMCIntervalPlot()
116 {
117  delete fParameters;
118  // kbelasco: why does deleting fPosteriorHist remove the graphics
119  // but deleting TGraphs doesn't?
120  //delete fPosteriorHist;
121  // can we delete fNLLHist and fWeightHist?
122  //delete fNLLHist;
123  //delete fWeightHist;
124 
125  // kbelasco: should we delete fPosteriorKeysPdf and fPosteriorKeysProduct?
126  delete fPosteriorKeysPdf;
127  delete fPosteriorKeysProduct;
128 
129  delete fWalk;
130  delete fBurnIn;
131  delete fFirst;
132  delete fParamGraph;
133  delete fNLLGraph;
134 }
135 
136 ////////////////////////////////////////////////////////////////////////////////
137 
138 void MCMCIntervalPlot::SetMCMCInterval(MCMCInterval& interval)
139 {
140  fInterval = &interval;
141  fDimension = fInterval->GetDimension();
142  fParameters = fInterval->GetParameters();
143 }
144 
145 ////////////////////////////////////////////////////////////////////////////////
146 
147 void MCMCIntervalPlot::Draw(const Option_t* options)
148 {
149  DrawInterval(options);
150 }
151 
152 ////////////////////////////////////////////////////////////////////////////////
153 
154 void MCMCIntervalPlot::DrawPosterior(const Option_t* options)
155 {
156  if (fInterval->GetUseKeys())
157  DrawPosteriorKeysPdf(options);
158  else
159  DrawPosteriorHist(options);
160 }
161 
162 ////////////////////////////////////////////////////////////////////////////////
163 
164 void* MCMCIntervalPlot::DrawPosteriorHist(const Option_t* /*options*/,
165  const char* title, Bool_t scale)
166 {
167  if (fPosteriorHist == NULL)
168  fPosteriorHist = fInterval->GetPosteriorHist();
169 
170  if (fPosteriorHist == NULL) {
171  coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorHist: "
172  << "Couldn't get posterior histogram." << endl;
173  return NULL;
174  }
175 
176  // kbelasco: annoying hack because histogram drawing fails when it sees
177  // an un-recognized option like POSTERIOR_HIST, etc.
178  //const Option_t* myOpt = NULL;
179 
180  //TString tmpOpt(options);
181  //if (tmpOpt.Contains("same"))
182  // myOpt = "same";
183 
184  // scale so highest bin has height 1
185  if (scale)
186  fPosteriorHist->Scale(1/fPosteriorHist->GetBinContent(
187  fPosteriorHist->GetMaximumBin()));
188 
189  TString ourTitle(GetTitle());
190  if (ourTitle.CompareTo("") == 0) {
191  if (title)
192  fPosteriorHist->SetTitle(title);
193  } else
194  fPosteriorHist->SetTitle(GetTitle());
195 
196  //fPosteriorHist->Draw(myOpt);
197 
198  return (void*)fPosteriorHist;
199 }
200 
201 ////////////////////////////////////////////////////////////////////////////////
202 
203 void* MCMCIntervalPlot::DrawPosteriorKeysPdf(const Option_t* options)
204 {
205  if (fPosteriorKeysPdf == NULL)
206  fPosteriorKeysPdf = fInterval->GetPosteriorKeysPdf();
207 
208  if (fPosteriorKeysPdf == NULL) {
209  coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
210  << "Couldn't get posterior Keys PDF." << endl;
211  return NULL;
212  }
213 
214  TString title(GetTitle());
215  Bool_t isEmpty = (title.CompareTo("") == 0);
216 
217  if (fDimension == 1) {
218  RooRealVar* v = (RooRealVar*)fParameters->first();
219  RooPlot* frame = v->frame();
220  if (frame == NULL) {
221  coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
222  << "Invalid parameter" << endl;
223  return NULL;
224  }
225  if (isEmpty)
226  frame->SetTitle(Form("Posterior Keys PDF for %s", v->GetName()));
227  else
228  frame->SetTitle(GetTitle());
229  //fPosteriorKeysPdf->plotOn(frame);
230  //fPosteriorKeysPdf->plotOn(frame,
231  // RooFit::Normalization(1, RooAbsReal::Raw));
232  //frame->Draw(options);
233  return (void*)frame;
234  } else if (fDimension == 2) {
235  RooArgList* axes = fInterval->GetAxes();
236  RooRealVar* xVar = (RooRealVar*)axes->at(0);
237  RooRealVar* yVar = (RooRealVar*)axes->at(1);
238  TH2F* keysHist = (TH2F*)fPosteriorKeysPdf->createHistogram(
239  "keysPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(kFALSE));
240  if (isEmpty)
241  keysHist->SetTitle(
242  Form("MCMC histogram of posterior Keys PDF for %s, %s",
243  axes->at(0)->GetName(), axes->at(1)->GetName()));
244  else
245  keysHist->SetTitle(GetTitle());
246 
247  keysHist->Draw(options);
248  delete axes;
249  return NULL;
250  }
251  return NULL;
252 }
253 
254 ////////////////////////////////////////////////////////////////////////////////
255 
256 void MCMCIntervalPlot::DrawInterval(const Option_t* options)
257 {
258  switch (fInterval->GetIntervalType()) {
259  case MCMCInterval::kShortest:
260  DrawShortestInterval(options);
261  break;
262  case MCMCInterval::kTailFraction:
263  DrawTailFractionInterval(options);
264  break;
265  default:
266  coutE(InputArguments) << "MCMCIntervalPlot::DrawInterval(): " <<
267  "Interval type not supported" << endl;
268  break;
269  }
270 }
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 
274 void MCMCIntervalPlot::DrawShortestInterval(const Option_t* options)
275 {
276  if (fInterval->GetUseKeys())
277  DrawKeysPdfInterval(options);
278  else
279  DrawHistInterval(options);
280 }
281 
282 ////////////////////////////////////////////////////////////////////////////////
283 
284 void MCMCIntervalPlot::DrawKeysPdfInterval(const Option_t* options)
285 {
286  TString title(GetTitle());
287  Bool_t isEmpty = (title.CompareTo("") == 0);
288 
289  if (fDimension == 1) {
290  // Draw the posterior keys PDF as well so the user can see where the
291  // limit bars line up
292  // fDimension == 1, so we know we will receive a RooPlot
293  RooPlot* frame = (RooPlot*)DrawPosteriorKeysPdf(options);
294 
295  //Double_t height = 1;
296  //Double_t height = 2.0 * fInterval->GetKeysPdfCutoff();
297  Double_t height = fInterval->GetKeysMax();
298 
299  RooRealVar* p = (RooRealVar*)fParameters->first();
300  Double_t ul = fInterval->UpperLimitByKeys(*p);
301  Double_t ll = fInterval->LowerLimitByKeys(*p);
302 
303  if (frame != NULL && fPosteriorKeysPdf != NULL) {
304  // draw shading in interval
305  if (isEmpty)
306  frame->SetTitle(NULL);
307  else
308  frame->SetTitle(GetTitle());
309  frame->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
310  p->GetName()));
311  fPosteriorKeysPdf->plotOn(frame,
313  RooFit::Range(ll, ul, kFALSE),
314  RooFit::VLines(),
315  RooFit::DrawOption("F"),
317  RooFit::FillColor(fShadeColor));
318 
319  // hack - this is drawn twice now:
320  // once by DrawPosteriorKeysPdf (which also configures things and sets
321  // the title), and once again here so the shading shows up behind.
322  fPosteriorKeysPdf->plotOn(frame,
324  }
325  if (frame) {
326  frame->Draw(options);
327  }
328 
329  TLine* llLine = new TLine(ll, 0, ll, height);
330  TLine* ulLine = new TLine(ul, 0, ul, height);
331  llLine->SetLineColor(fLineColor);
332  ulLine->SetLineColor(fLineColor);
333  llLine->SetLineWidth(fLineWidth);
334  ulLine->SetLineWidth(fLineWidth);
335  llLine->Draw(options);
336  ulLine->Draw(options);
337  } else if (fDimension == 2) {
338  if (fPosteriorKeysPdf == NULL)
339  fPosteriorKeysPdf = fInterval->GetPosteriorKeysPdf();
340 
341  if (fPosteriorKeysPdf == NULL) {
342  coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
343  << "Couldn't get posterior Keys PDF." << endl;
344  return;
345  }
346 
347  RooArgList* axes = fInterval->GetAxes();
348  RooRealVar* xVar = (RooRealVar*)axes->at(0);
349  RooRealVar* yVar = (RooRealVar*)axes->at(1);
350  TH2F* contHist = (TH2F*)fPosteriorKeysPdf->createHistogram(
351  "keysContour2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(kFALSE));
352  //if (isEmpty)
353  // contHist->SetTitle(Form("MCMC Keys conf. interval for %s, %s",
354  // axes->at(0)->GetName(), axes->at(1)->GetName()));
355  //else
356  // contHist->SetTitle(GetTitle());
357  if (!isEmpty)
358  contHist->SetTitle(GetTitle());
359  else
360  contHist->SetTitle(NULL);
361 
362  contHist->SetStats(kFALSE);
363 
364  TString tmpOpt(options);
365  if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
366 
367  Double_t cutoff = fInterval->GetKeysPdfCutoff();
368  contHist->SetContour(1, &cutoff);
369  contHist->SetLineColor(fLineColor);
370  contHist->SetLineWidth(fLineWidth);
371  contHist->Draw(tmpOpt.Data());
372  delete axes;
373  } else {
374  coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
375  << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
376  }
377 }
378 
379 ////////////////////////////////////////////////////////////////////////////////
380 
381 void MCMCIntervalPlot::DrawHistInterval(const Option_t* options)
382 {
383  TString title(GetTitle());
384  Bool_t isEmpty = (title.CompareTo("") == 0);
385 
386  if (fDimension == 1) {
387  // draw lower and upper limits
388  RooRealVar* p = (RooRealVar*)fParameters->first();
389  Double_t ul = fInterval->UpperLimitByHist(*p);
390  Double_t ll = fInterval->LowerLimitByHist(*p);
391 
392  // Draw the posterior histogram as well so the user can see where the
393  // limit bars line up
394  // fDimension == 1, so we know will get a TH1F*
395  TH1F* hist = (TH1F*)DrawPosteriorHist(options, NULL, false);
396  if (hist == NULL) return;
397  if (isEmpty)
398  hist->SetTitle(NULL);
399  else
400  hist->SetTitle(GetTitle());
401  hist->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
402  p->GetName()));
403  hist->SetStats(kFALSE);
404  TH1F* copy = (TH1F*)hist->Clone(Form("%s_copy", hist->GetTitle()));
405  Double_t histCutoff = fInterval->GetHistCutoff();
406 
407  Int_t i;
408  Int_t nBins = copy->GetNbinsX();
409  Double_t height;
410  for (i = 1; i <= nBins; i++) {
411  // remove bins with height < cutoff
412  height = copy->GetBinContent(i);
413  if (height < histCutoff) {
414  copy->SetBinContent(i, 0);
415  copy->SetBinError(i, 0);
416  }
417  }
418 
419  hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
420  copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
421 
422  copy->SetFillStyle(1001);
423  copy->SetFillColor(fShadeColor);
424  hist->Draw(options);
425  // to show the interval area
426  copy->Draw("HIST SAME");
427 
428  fPosteriorHistHistCopy = copy;
429 
430  TLine* llLine = new TLine(ll, 0, ll, 1);
431  TLine* ulLine = new TLine(ul, 0, ul, 1);
432  llLine->SetLineColor(fLineColor);
433  ulLine->SetLineColor(fLineColor);
434  llLine->SetLineWidth(fLineWidth);
435  ulLine->SetLineWidth(fLineWidth);
436  llLine->Draw(options);
437  ulLine->Draw(options);
438 
439  } else if (fDimension == 2) {
440  if (fPosteriorHist == NULL)
441  fPosteriorHist = fInterval->GetPosteriorHist();
442 
443  if (fPosteriorHist == NULL) {
444  coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
445  << "Couldn't get posterior histogram." << endl;
446  return;
447  }
448 
449  RooArgList* axes = fInterval->GetAxes();
450  //if (isEmpty)
451  // fPosteriorHist->SetTitle(
452  // Form("MCMC histogram conf. interval for %s, %s",
453  // axes->at(0)->GetName(), axes->at(1)->GetName()));
454  //else
455  // fPosteriorHist->SetTitle(GetTitle());
456  if (!isEmpty)
457  fPosteriorHist->SetTitle(GetTitle());
458  else
459  fPosteriorHist->SetTitle(NULL);
460  delete axes;
461 
462  fPosteriorHist->SetStats(kFALSE);
463 
464  TString tmpOpt(options);
465  if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
466 
467  Double_t cutoff = fInterval->GetHistCutoff();
468  fPosteriorHist->SetContour(1, &cutoff);
469  fPosteriorHist->SetLineColor(fLineColor);
470  fPosteriorHist->SetLineWidth(fLineWidth);
471  fPosteriorHist->Draw(tmpOpt.Data());
472  } else {
473  coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
474  << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
475  }
476 }
477 
478 ////////////////////////////////////////////////////////////////////////////////
479 
480 void MCMCIntervalPlot::DrawTailFractionInterval(const Option_t* options)
481 {
482  TString title(GetTitle());
483  Bool_t isEmpty = (title.CompareTo("") == 0);
484 
485  if (fDimension == 1) {
486  // Draw the posterior histogram as well so the user can see where the
487  // limit bars line up
488  RooRealVar* p = (RooRealVar*)fParameters->first();
489  Double_t ul = fInterval->UpperLimitTailFraction(*p);
490  Double_t ll = fInterval->LowerLimitTailFraction(*p);
491 
492  TH1F* hist = (TH1F*)DrawPosteriorHist(options, NULL, false);
493  if (hist == NULL) return;
494  if (isEmpty)
495  hist->SetTitle(NULL);
496  else
497  hist->SetTitle(GetTitle());
498  hist->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
499  p->GetName()));
500  hist->SetStats(kFALSE);
501  TH1F* copy = (TH1F*)hist->Clone(Form("%s_copy", hist->GetTitle()));
502 
503  Int_t i;
504  Int_t nBins = copy->GetNbinsX();
505  Double_t center;
506  for (i = 1; i <= nBins; i++) {
507  // remove bins outside interval
508  center = copy->GetBinCenter(i);
509  if (center < ll || center > ul) {
510  copy->SetBinContent(i, 0);
511  copy->SetBinError(i, 0);
512  }
513  }
514 
515  hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
516  copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
517 
518  copy->SetFillStyle(1001);
519  copy->SetFillColor(fShadeColor);
520  hist->Draw(options);
521  copy->Draw("hist same");
522 
523  // draw lower and upper limits
524  TLine* llLine = new TLine(ll, 0, ll, 1);
525  TLine* ulLine = new TLine(ul, 0, ul, 1);
526  llLine->SetLineColor(fLineColor);
527  ulLine->SetLineColor(fLineColor);
528  llLine->SetLineWidth(fLineWidth);
529  ulLine->SetLineWidth(fLineWidth);
530  llLine->Draw(options);
531  ulLine->Draw(options);
532  } else {
533  coutE(InputArguments) << "MCMCIntervalPlot::DrawTailFractionInterval: "
534  << " Sorry: " << fDimension << "-D plots not currently supported"
535  << endl;
536  }
537 }
538 
539 ////////////////////////////////////////////////////////////////////////////////
540 
541 void* MCMCIntervalPlot::DrawPosteriorKeysProduct(const Option_t* options)
542 {
543  if (fPosteriorKeysProduct == NULL)
544  fPosteriorKeysProduct = fInterval->GetPosteriorKeysProduct();
545 
546  if (fPosteriorKeysProduct == NULL) {
547  coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysProduct: "
548  << "Couldn't get posterior Keys product." << endl;
549  return NULL;
550  }
551 
552  RooArgList* axes = fInterval->GetAxes();
553 
554  TString title(GetTitle());
555  Bool_t isEmpty = (title.CompareTo("") == 0);
556 
557  if (fDimension == 1) {
558  RooPlot* frame = ((RooRealVar*)fParameters->first())->frame();
559  if (!frame) return NULL;
560  if (isEmpty)
561  frame->SetTitle(Form("Posterior Keys PDF * Heaviside product for %s",
562  axes->at(0)->GetName()));
563  else
564  frame->SetTitle(GetTitle());
565  //fPosteriorKeysProduct->plotOn(frame);
566  fPosteriorKeysProduct->plotOn(frame,
568  frame->Draw(options);
569  return (void*)frame;
570  } else if (fDimension == 2) {
571  RooRealVar* xVar = (RooRealVar*)axes->at(0);
572  RooRealVar* yVar = (RooRealVar*)axes->at(1);
573  TH2F* productHist = (TH2F*)fPosteriorKeysProduct->createHistogram(
574  "prodPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(kFALSE));
575  if (isEmpty)
576  productHist->SetTitle(
577  Form("MCMC Posterior Keys Product Hist. for %s, %s",
578  axes->at(0)->GetName(), axes->at(1)->GetName()));
579  else
580  productHist->SetTitle(GetTitle());
581  productHist->Draw(options);
582  return NULL;
583  }
584  delete axes;
585  return NULL;
586 }
587 
588 ////////////////////////////////////////////////////////////////////////////////
589 
590 void MCMCIntervalPlot::DrawChainScatter(RooRealVar& xVar, RooRealVar& yVar)
591 {
592  const MarkovChain* markovChain = fInterval->GetChain();
593 
594  Int_t size = markovChain->Size();
595  Int_t burnInSteps;
596  if (fShowBurnIn)
597  burnInSteps = fInterval->GetNumBurnInSteps();
598  else
599  burnInSteps = 0;
600 
601  Double_t* x = new Double_t[size - burnInSteps];
602  Double_t* y = new Double_t[size - burnInSteps];
603  Double_t* burnInX = NULL;
604  Double_t* burnInY = NULL;
605  if (burnInSteps > 0) {
606  burnInX = new Double_t[burnInSteps];
607  burnInY = new Double_t[burnInSteps];
608  }
609  Double_t firstX;
610  Double_t firstY;
611 
612  for (Int_t i = burnInSteps; i < size; i++) {
613  x[i - burnInSteps] = markovChain->Get(i)->getRealValue(xVar.GetName());
614  y[i - burnInSteps] = markovChain->Get(i)->getRealValue(yVar.GetName());
615  }
616 
617  for (Int_t i = 0; i < burnInSteps; i++) {
618  burnInX[i] = markovChain->Get(i)->getRealValue(xVar.GetName());
619  burnInY[i] = markovChain->Get(i)->getRealValue(yVar.GetName());
620  }
621 
622  firstX = markovChain->Get(0)->getRealValue(xVar.GetName());
623  firstY = markovChain->Get(0)->getRealValue(yVar.GetName());
624 
625  TString title(GetTitle());
626  Bool_t isEmpty = (title.CompareTo("") == 0);
627 
628  TGraph* walk = new TGraph(size - burnInSteps, x, y);
629  if (isEmpty)
630  walk->SetTitle(Form("2-D Scatter Plot of Markov chain for %s, %s",
631  xVar.GetName(), yVar.GetName()));
632  else
633  walk->SetTitle(GetTitle());
634  // kbelasco: figure out how to set TGraph variable ranges
635  walk->GetXaxis()->Set(xVar.numBins(), xVar.getMin(), xVar.getMax());
636  walk->GetXaxis()->SetTitle(xVar.GetName());
637  walk->GetYaxis()->Set(yVar.numBins(), yVar.getMin(), yVar.getMax());
638  walk->GetYaxis()->SetTitle(yVar.GetName());
639  walk->SetLineColor(kGray);
640  walk->SetMarkerStyle(6);
641  walk->SetMarkerColor(kViolet);
642  walk->Draw("A,L,P,same");
643 
644  TGraph* burnIn = NULL;
645  if (burnInX != NULL && burnInY != NULL) {
646  burnIn = new TGraph(burnInSteps - 1, burnInX, burnInY);
647  burnIn->SetLineColor(kPink);
648  burnIn->SetMarkerStyle(6);
649  burnIn->SetMarkerColor(kPink);
650  burnIn->Draw("L,P,same");
651  }
652 
653  TGraph* first = new TGraph(1, &firstX, &firstY);
654  first->SetLineColor(kGreen);
655  first->SetMarkerStyle(3);
656  first->SetMarkerSize(2);
657  first->SetMarkerColor(kGreen);
658  first->Draw("L,P,same");
659 
660  //walkCanvas->Update();
661  delete [] x;
662  delete [] y;
663  if (burnInX != NULL) delete [] burnInX;
664  if (burnInY != NULL) delete [] burnInY;
665  //delete walk;
666  //delete burnIn;
667  //delete first;
668 }
669 
670 ////////////////////////////////////////////////////////////////////////////////
671 
672 void MCMCIntervalPlot::DrawParameterVsTime(RooRealVar& param)
673 {
674  const MarkovChain* markovChain = fInterval->GetChain();
675  Int_t size = markovChain->Size();
676  Int_t numEntries = 2 * size;
677  Double_t* value = new Double_t[numEntries];
678  Double_t* time = new Double_t[numEntries];
679  Double_t val;
680  Int_t weight;
681  Int_t t = 0;
682  for (Int_t i = 0; i < size; i++) {
683  val = markovChain->Get(i)->getRealValue(param.GetName());
684  weight = (Int_t)markovChain->Weight();
685  value[2*i] = val;
686  value[2*i + 1] = val;
687  time[2*i] = t;
688  t += weight;
689  time[2*i + 1] = t;
690  }
691 
692  TString title(GetTitle());
693  Bool_t isEmpty = (title.CompareTo("") == 0);
694 
695  TGraph* paramGraph = new TGraph(numEntries, time, value);
696  if (isEmpty)
697  paramGraph->SetTitle(Form("%s vs. time in Markov chain",param.GetName()));
698  else
699  paramGraph->SetTitle(GetTitle());
700  paramGraph->GetXaxis()->SetTitle("Time (discrete steps)");
701  paramGraph->GetYaxis()->SetTitle(param.GetName());
702  //paramGraph->SetLineColor(fLineColor);
703  paramGraph->Draw("A,L,same");
704  delete [] value;
705  delete [] time;
706  //gPad->Update();
707 }
708 
709 ////////////////////////////////////////////////////////////////////////////////
710 
711 void MCMCIntervalPlot::DrawNLLVsTime()
712 {
713  const MarkovChain* markovChain = fInterval->GetChain();
714  Int_t size = markovChain->Size();
715  Int_t numEntries = 2 * size;
716  Double_t* nllValue = new Double_t[numEntries];
717  Double_t* time = new Double_t[numEntries];
718  Double_t nll;
719  Int_t weight;
720  Int_t t = 0;
721  for (Int_t i = 0; i < size; i++) {
722  nll = markovChain->NLL(i);
723  weight = (Int_t)markovChain->Weight();
724  nllValue[2*i] = nll;
725  nllValue[2*i + 1] = nll;
726  time[2*i] = t;
727  t += weight;
728  time[2*i + 1] = t;
729  }
730 
731  TString title(GetTitle());
732  Bool_t isEmpty = (title.CompareTo("") == 0);
733 
734  TGraph* nllGraph = new TGraph(numEntries, time, nllValue);
735  if (isEmpty)
736  nllGraph->SetTitle("NLL value vs. time in Markov chain");
737  else
738  nllGraph->SetTitle(GetTitle());
739  nllGraph->GetXaxis()->SetTitle("Time (discrete steps)");
740  nllGraph->GetYaxis()->SetTitle("NLL (-log(likelihood))");
741  //nllGraph->SetLineColor(fLineColor);
742  nllGraph->Draw("A,L,same");
743  //gPad->Update();
744  delete [] nllValue;
745  delete [] time;
746 }
747 
748 ////////////////////////////////////////////////////////////////////////////////
749 
750 void MCMCIntervalPlot::DrawNLLHist(const Option_t* options)
751 {
752  if (fNLLHist == NULL) {
753  const MarkovChain* markovChain = fInterval->GetChain();
754  // find the max NLL value
755  Double_t maxNLL = 0;
756  Int_t size = markovChain->Size();
757  for (Int_t i = 0; i < size; i++)
758  if (markovChain->NLL(i) > maxNLL)
759  maxNLL = markovChain->NLL(i);
760  RooRealVar* nllVar = fInterval->GetNLLVar();
761  fNLLHist = new TH1F("mcmc_nll_hist", "MCMC NLL Histogram",
762  nllVar->getBins(), 0, maxNLL);
763  TString title(GetTitle());
764  Bool_t isEmpty = (title.CompareTo("") == 0);
765  if (!isEmpty)
766  fNLLHist->SetTitle(GetTitle());
767  fNLLHist->GetXaxis()->SetTitle("-log(likelihood)");
768  for (Int_t i = 0; i < size; i++)
769  fNLLHist->Fill(markovChain->NLL(i), markovChain->Weight());
770  }
771  fNLLHist->Draw(options);
772 }
773 
774 ////////////////////////////////////////////////////////////////////////////////
775 
776 void MCMCIntervalPlot::DrawWeightHist(const Option_t* options)
777 {
778  if (fWeightHist == NULL) {
779  const MarkovChain* markovChain = fInterval->GetChain();
780  // find the max weight value
781  Double_t maxWeight = 0;
782  Int_t size = markovChain->Size();
783  for (Int_t i = 0; i < size; i++)
784  if (markovChain->Weight(i) > maxWeight)
785  maxWeight = markovChain->Weight(i);
786  fWeightHist = new TH1F("mcmc_weight_hist", "MCMC Weight Histogram",
787  (Int_t)(maxWeight + 1), 0, maxWeight * 1.02);
788  for (Int_t i = 0; i < size; i++)
789  fWeightHist->Fill(markovChain->Weight(i));
790  }
791  fWeightHist->Draw(options);
792 }
793 
794 /*
795 /////////////////////////////////////////////////////////////////////
796  // 3-d plot of the parameter points
797  dataCanvas->cd(2);
798  // also plot the points in the markov chain
799  RooDataSet* markovChainData = ((MCMCInterval*)mcmcint)->GetChainAsDataSet();
800 
801  TTree& chain = ((RooTreeDataStore*) markovChainData->store())->tree();
802  chain.SetMarkerStyle(6);
803  chain.SetMarkerColor(kRed);
804  chain.Draw("s:ratioSigEff:ratioBkgEff","","box"); // 3-d box proportional to posterior
805 
806  // the points used in the profile construction
807  TTree& parameterScan = ((RooTreeDataStore*) fc.GetPointsToScan()->store())->tree();
808  parameterScan.SetMarkerStyle(24);
809  parameterScan.Draw("s:ratioSigEff:ratioBkgEff","","same");
810 
811  chain.SetMarkerStyle(6);
812  chain.SetMarkerColor(kRed);
813  //chain.Draw("s:ratioSigEff:ratioBkgEff", "_MarkovChain_local_nll","box");
814  //chain.Draw("_MarkovChain_local_nll");
815 ////////////////////////////////////////////////////////////////////
816 */
first
Definition: first.py:1
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
TH1F.h
RooStats::MCMCInterval
Definition: MCMCInterval.h:32
TLine.h
TGraph::SetTitle
virtual void SetTitle(const char *title="")
Change (i.e.
Definition: TGraph.cxx:2324
TGraph::Draw
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:760
TH1::SetContour
virtual void SetContour(Int_t nlevels, const Double_t *levels=0)
Set the number and values of contour levels.
Definition: TH1.cxx:7946
TAxis::Set
virtual void Set(Int_t nbins, Double_t xmin, Double_t xmax)
Initialize axis with fix bins.
Definition: TAxis.cxx:731
TLine
Definition: TLine.h:22
RooAbsRealLValue::getMax
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
Definition: RooAbsRealLValue.h:89
RooArgSet::getRealValue
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:474
RooFit::InputArguments
@ InputArguments
Definition: RooGlobalFunc.h:68
kGreen
@ kGreen
Definition: Rtypes.h:66
TH2F
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
TString::Data
const char * Data() const
Definition: TString.h:369
RooStats::MarkovChain::Size
virtual Int_t Size() const
get the number of steps in the chain
Definition: MarkovChain.h:61
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
Form
char * Form(const char *fmt,...)
TNamed::GetTitle
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:54
TGraph.h
kViolet
@ kViolet
Definition: Rtypes.h:67
coutE
#define coutE(a)
Definition: RooMsgService.h:33
RooStats::MarkovChain::Get
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
Definition: MarkovChain.h:63
RooArgList
Definition: RooArgList.h:21
RooStats::MarkovChain
Definition: MarkovChain.h:36
kPink
@ kPink
Definition: Rtypes.h:67
Int_t
int Int_t
Definition: RtypesCore.h:45
RooFit::Normalization
RooCmdArg Normalization(Double_t scaleFactor)
Definition: RooGlobalFunc.cxx:51
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
TH1::SetBinContent
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8677
x
Double_t x[n]
Definition: legend1.C:17
TList.h
RooArgList::at
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:88
RooFit::YVar
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
Definition: RooGlobalFunc.cxx:240
RooFit::FillColor
RooCmdArg FillColor(Color_t color)
Definition: RooGlobalFunc.cxx:59
TString
Definition: TString.h:136
RooFit::MoveToBack
RooCmdArg MoveToBack()
Definition: RooGlobalFunc.cxx:66
TGraph::GetYaxis
TAxis * GetYaxis() const
Get y axis of the graph.
Definition: TGraph.cxx:1636
RooStats::MCMCIntervalPlot
Definition: MCMCIntervalPlot.h:37
RooStats::MCMCInterval::GetDimension
virtual Int_t GetDimension() const
Get the number of parameters of interest in this interval.
Definition: MCMCInterval.h:193
v
@ v
Definition: rootcling_impl.cxx:3635
bool
RooAbsRealLValue::getMin
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
Definition: RooAbsRealLValue.h:86
TH1::SetTitle
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6344
TH1::Clone
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2664
kBlack
@ kBlack
Definition: Rtypes.h:65
TGraph::GetXaxis
TAxis * GetXaxis() const
Get x axis of the graph.
Definition: TGraph.cxx:1626
MarkovChain.h
Option_t
const typedef char Option_t
Definition: RtypesCore.h:66
RooFit::DrawOption
RooCmdArg DrawOption(const char *opt)
Definition: RooGlobalFunc.cxx:38
RooAbsRealLValue::numBins
virtual Int_t numBins(const char *rangeName=0) const
Definition: RooAbsRealLValue.h:54
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4906
TH1::GetYaxis
TAxis * GetYaxis()
Definition: TH1.h:318
TH1::GetMaximumBin
virtual Int_t GetMaximumBin() const
Return location of bin with maximum value in the range.
Definition: TH1.cxx:8035
MCMCIntervalPlot.h
Draw
th1 Draw()
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
TString::Append
TString & Append(const char *cs)
Definition: TString.h:564
MCMCInterval.h
TH1::GetBinCenter
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8596
RooPlot.h
RooPlot::GetYaxis
TAxis * GetYaxis() const
Definition: RooPlot.cxx:1258
RooPlot
Definition: RooPlot.h:44
y
Double_t y[n]
Definition: legend1.C:17
RooStats::MarkovChain::Weight
virtual Double_t Weight() const
get the weight of the current (last indexed) entry
Definition: MarkovChain.cxx:279
RooRealVar.h
TNamed::SetTitle
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
TH1::SetStats
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8445
TH2.h
RooGlobalFunc.h
RooAbsReal::Raw
@ Raw
Definition: RooAbsReal.h:253
RooStats::MarkovChain::NLL
virtual Double_t NLL(Int_t i) const
get the NLL value of entry at position i
Definition: MarkovChain.cxx:263
TObject::Draw
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:197
TString::CompareTo
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition: TString.cxx:418
Double_t
double Double_t
Definition: RtypesCore.h:59
TGraph
Definition: TGraph.h:41
RooPlot::SetTitle
void SetTitle(const char *name)
Set the title of the RooPlot to 'title'.
Definition: RooPlot.cxx:1237
TH1::SetBinError
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition: TH1.cxx:8661
RooStats
Definition: Asimov.h:19
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:572
TAxis.h
RooFit::Range
RooCmdArg Range(const char *rangeName, Bool_t adjustNorm=kTRUE)
Definition: RooGlobalFunc.cxx:52
fFirst
T1 fFirst
Definition: X11Events.mm:86
kGray
@ kGray
Definition: Rtypes.h:65
RooAbsRealLValue::getBins
virtual Int_t getBins(const char *name=0) const
Get number of bins of currently defined range.
Definition: RooAbsRealLValue.h:83
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooRealVar
Definition: RooRealVar.h:36
RooFit::VLines
RooCmdArg VLines()
Definition: RooGlobalFunc.cxx:55
RooArgList.h
RooFit::Scaling
RooCmdArg Scaling(Bool_t flag)
Definition: RooGlobalFunc.cxx:243
TH1::Scale
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6245
TH1::GetNbinsX
virtual Int_t GetNbinsX() const
Definition: TH1.h:293
int
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2997