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