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