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