StandardProfileLikelihoodDemo.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roostats
/// \notebook
/// Standard demo of the Profile Likelihood calculator
/// StandardProfileLikelihoodDemo
///
/// 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 ProfileLikelihoodCalculator is based on Wilks's theorem
/// and the asymptotic properties of the profile likelihood ratio
/// (eg. that it is chi-square distributed for the true value).
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \author Kyle Cranmer

#include "TFile.h"
#include "TROOT.h"
#include "TSystem.h"
#include "RooWorkspace.h"
#include "RooAbsData.h"
#include "RooRealVar.h"

#include "RooStats/ModelConfig.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodInterval.h"
#include "RooStats/LikelihoodIntervalPlot.h"

using namespace RooFit;
using namespace RooStats;

struct ProfileLikelihoodOptions {

double confLevel = 0.95;
int nScanPoints = 50;
bool plotAsTF1 = false;
double poiMinPlot = 1;
double poiMaxPlot = 0;
   bool doHypoTest = false;
   double nullValue = 0;

};

ProfileLikelihoodOptions optPL;

void StandardProfileLikelihoodDemo(const char* infile = "",
                                   const char* workspaceName = "combined",
                                   const char* modelConfigName = "ModelConfig",
                                   const char* dataName = "obsData"){

      double confLevel = optPL.confLevel;
      double nScanPoints = optPL.nScanPoints;
      bool plotAsTF1 = optPL.plotAsTF1;
      double poiXMin = optPL.poiMinPlot;
      double poiXMax = optPL.poiMaxPlot;
      bool doHypoTest = optPL.doHypoTest;
      double nullParamValue = optPL.nullValue;

   // -------------------------------------------------------
   // 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;
   }

   // ---------------------------------------------
   // create and use the ProfileLikelihoodCalculator
   // to find and plot the 95% confidence interval
   // on the parameter of interest as specified
   // in the model config
   ProfileLikelihoodCalculator pl(*data,*mc);
   pl.SetConfidenceLevel(confLevel); // 95% interval
   LikelihoodInterval* interval = pl.GetInterval();

   // print out the interval on the first Parameter of Interest
   RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
   cout << "\n>>>> RESULT : " << confLevel*100 << "% interval on " <<firstPOI->GetName()<<" is : ["<<
      interval->LowerLimit(*firstPOI) << ", "<<
      interval->UpperLimit(*firstPOI) <<"]\n "<<endl;

   // make a plot

   cout << "making a plot of the profile likelihood function ....(if it is taking a lot of time use less points or the TF1 drawing option)\n";
   LikelihoodIntervalPlot plot(interval);
   plot.SetNPoints(nScanPoints);  // do not use too many points, it could become very slow for some models
   if (poiXMin < poiXMax) plot.SetRange(poiXMin, poiXMax);
   TString opt;
   if (plotAsTF1) opt += TString("tf1");
   plot.Draw(opt);  // use option TF1 if too slow (plot.Draw("tf1")


   // if requested perform also an hypothesis test for the significance
   if (doHypoTest) {
      RooArgSet nullparams("nullparams");
      nullparams.addClone(*firstPOI);
      nullparams.setRealValue(firstPOI->GetName(), nullParamValue);
      pl.SetNullParameters(nullparams);
      std::cout << "Perform Test of Hypothesis : null Hypothesis is " << firstPOI->GetName() << nullParamValue << std::endl;
      auto result = pl.GetHypoTest();
      std::cout << "\n>>>> Hypotheis Test Result \n";
      result->Print();
   }


}