StandardProfileInspectorDemo.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roostats
/// \notebook -js
/// Standard demo of the ProfileInspector class
/// StandardProfileInspectorDemo
///
/// This is a standard demo that can be used with any ROOT file
/// prepared in the standard way.  You specify:
///  - name for input ROOT file
///  - name of workspace inside ROOT file that holds model and data
///  - name of ModelConfig that specifies details for calculator tools
///  - name of dataset
///
/// With default parameters the macro will attempt to run the
/// standard hist2workspace example and read the ROOT file
/// that it produces.
///
/// The actual heart of the demo is only about 10 lines long.
///
/// The ProfileInspector plots the conditional maximum likelihood estimate
/// of each nuisance parameter in the model vs. the parameter of interest.
/// (aka. profiled value of nuisance parameter vs. parameter of interest)
/// (aka. best fit nuisance parameter with p.o.i fixed vs. parameter of interest)
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \author Kyle Cranmer
 
#include "TFile.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TList.h"
#include "TMath.h"
#include "TSystem.h"
#include "RooWorkspace.h"
#include "RooAbsData.h"
 
#include "RooStats/ModelConfig.h"
#include "RooStats/ProfileInspector.h"
 
using namespace RooFit;
using namespace RooStats;
 
void StandardProfileInspectorDemo(const char *infile = "", const char *workspaceName = "combined",
                                  const char *modelConfigName = "ModelConfig", const char *dataName = "obsData")
{
 
   // -------------------------------------------------------
   // First part is just to access a user-defined file
   // or create the standard example file if it doesn't exist
 
   const char *filename = "";
   if (!strcmp(infile, "")) {
      filename = "results/example_combined_GaussExample_model.root";
      bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
      // if file does not exists generate with histfactory
      if (!fileExist) {
#ifdef _WIN32
         cout << "HistFactory file cannot be generated on Windows - exit" << endl;
         return;
#endif
         // Normally this would be run on the command line
         cout << "will run standard hist2workspace example" << endl;
         gROOT->ProcessLine(".! prepareHistFactory .");
         gROOT->ProcessLine(".! hist2workspace config/example.xml");
         cout << "\n\n---------------------" << endl;
         cout << "Done creating example input" << endl;
         cout << "---------------------\n\n" << endl;
      }
 
   } else
      filename = infile;
 
   // Try to open the file
   TFile *file = TFile::Open(filename);
 
   // if input file was specified byt not found, quit
   if (!file) {
      cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
      return;
   }
 
   // -------------------------------------------------------
   // Tutorial starts here
   // -------------------------------------------------------
 
   // get the workspace out of the file
   RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
   if (!w) {
      cout << "workspace not found" << endl;
      return;
   }
 
   // get the modelConfig out of the file
   ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
 
   // get the modelConfig out of the file
   RooAbsData *data = w->data(dataName);
 
   // make sure ingredients are found
   if (!data || !mc) {
      w->Print();
      cout << "data or ModelConfig was not found" << endl;
      return;
   }
 
   // -----------------------------
   // now use the profile inspector
   ProfileInspector p;
   TList *list = p.GetListOfProfilePlots(*data, mc);
 
   // now make plots
   TCanvas *c1 = new TCanvas("c1", "ProfileInspectorDemo", 800, 200);
   if (list->GetSize() > 4) {
      double n = list->GetSize();
      int nx = (int)sqrt(n);
      int ny = TMath::CeilNint(n / nx);
      nx = TMath::CeilNint(sqrt(n));
      c1->Divide(ny, nx);
   } else
      c1->Divide(list->GetSize());
   for (int i = 0; i < list->GetSize(); ++i) {
      c1->cd(i + 1);
      list->At(i)->Draw("al");
   }
 
   cout << endl;
}