StandardTestStatDistributionDemo.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roostats
/// \notebook
/// StandardTestStatDistributionDemo.C
///
/// This simple script plots the sampling distribution of the profile likelihood
/// ratio test statistic based on the input Model File.  To do this one needs to
/// specify the value of the parameter of interest that will be used for evaluating
/// the test statistic and the value of the parameters used for generating the toy data.
/// In this case, it uses the upper-limit estimated from the ProfileLikleihoodCalculator,
/// which assumes the asymptotic chi-square distribution for -2 log profile likelihood ratio.
/// Thus, the script is handy for checking to see if the asymptotic approximations are valid.
/// To aid, that comparison, the script overlays a chi-square distribution as well.
/// The most common parameter of interest is a parameter proportional to the signal rate,
/// and often that has a lower-limit of 0, which breaks the standard chi-square distribution.
/// Thus the script allows the parameter to be negative so that the overlay chi-square is
/// the correct asymptotic distribution.
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \author Kyle Cranmer
 
#include "TFile.h"
#include "TROOT.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TSystem.h"
#include "TF1.h"
#include "TSystem.h"
 
#include "RooWorkspace.h"
#include "RooAbsData.h"
 
#include "RooStats/ModelConfig.h"
#include "RooStats/FeldmanCousins.h"
#include "RooStats/ToyMCSampler.h"
#include "RooStats/PointSetInterval.h"
#include "RooStats/ConfidenceBelt.h"
 
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodInterval.h"
#include "RooStats/ProfileLikelihoodTestStat.h"
#include "RooStats/SamplingDistribution.h"
#include "RooStats/SamplingDistPlot.h"
 
using namespace RooFit;
using namespace RooStats;
 
 
bool useProof = false;  // flag to control whether to use Proof
int nworkers = 0;   // number of workers (default use all available cores)
 
// -------------------------------------------------------
// The actual macro
 
void StandardTestStatDistributionDemo(const char* infile = "",
                                      const char* workspaceName = "combined",
                                      const char* modelConfigName = "ModelConfig",
                                      const char* dataName = "obsData"){
 
 
   // the number of toy MC used to generate the distribution
   int nToyMC = 1000;
   // The parameter below is needed for asymptotic distribution to be chi-square,
   // but set to false if your model is not numerically stable if mu<0
   bool allowNegativeMu=true;
 
 
   // -------------------------------------------------------
   // 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;
      }
 
 
   // -------------------------------------------------------
   // Now get the data and workspace
 
   // 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;
   }
 
   mc->Print();
   // -------------------------------------------------------
   // Now find the upper limit based on the asymptotic results
   RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
   ProfileLikelihoodCalculator plc(*data,*mc);
   LikelihoodInterval* interval = plc.GetInterval();
   double plcUpperLimit = interval->UpperLimit(*firstPOI);
   delete interval;
   cout << "\n\n--------------------------------------"<<endl;
   cout <<"Will generate sampling distribution at " << firstPOI->GetName() << " = " << plcUpperLimit <<endl;
   int nPOI = mc->GetParametersOfInterest()->getSize();
   if(nPOI>1){
      cout <<"not sure what to do with other parameters of interest, but here are their values"<<endl;
      mc->GetParametersOfInterest()->Print("v");
   }
 
   // -------------------------------------------------------
   // create the test stat sampler
   ProfileLikelihoodTestStat ts(*mc->GetPdf());
 
   // to avoid effects from boundary and simplify asymptotic comparison, set min=-max
   if(allowNegativeMu)
      firstPOI->setMin(-1*firstPOI->getMax());
 
   // temporary RooArgSet
   RooArgSet poi;
   poi.add(*mc->GetParametersOfInterest());
 
   // create and configure the ToyMCSampler
   ToyMCSampler sampler(ts,nToyMC);
   sampler.SetPdf(*mc->GetPdf());
   sampler.SetObservables(*mc->GetObservables());
   sampler.SetGlobalObservables(*mc->GetGlobalObservables());
   if(!mc->GetPdf()->canBeExtended() && (data->numEntries()==1)){
      cout << "tell it to use 1 event"<<endl;
      sampler.SetNEventsPerToy(1);
   }
   firstPOI->setVal(plcUpperLimit); // set POI value for generation
   sampler.SetParametersForTestStat(*mc->GetParametersOfInterest()); // set POI value for evaluation
 
   if (useProof) {
      ProofConfig pc(*w, nworkers, "",false);
      sampler.SetProofConfig(&pc); // enable proof
   }
 
   firstPOI->setVal(plcUpperLimit);
   RooArgSet allParameters;
   allParameters.add(*mc->GetParametersOfInterest());
   allParameters.add(*mc->GetNuisanceParameters());
   allParameters.Print("v");
 
   SamplingDistribution* sampDist = sampler.GetSamplingDistribution(allParameters);
   SamplingDistPlot plot;
   plot.AddSamplingDistribution(sampDist);
   plot.GetTH1F(sampDist)->GetYaxis()->SetTitle(Form("f(-log #lambda(#mu=%.2f) | #mu=%.2f)",plcUpperLimit,plcUpperLimit));
   plot.SetAxisTitle(Form("-log #lambda(#mu=%.2f)",plcUpperLimit));
 
   TCanvas* c1 = new TCanvas("c1");
   c1->SetLogy();
   plot.Draw();
   double min = plot.GetTH1F(sampDist)->GetXaxis()->GetXmin();
   double max = plot.GetTH1F(sampDist)->GetXaxis()->GetXmax();
 
   TF1* f = new TF1("f",Form("2*ROOT::Math::chisquared_pdf(2*x,%d,0)",nPOI),min,max);
   f->Draw("same");
   c1->SaveAs("standard_test_stat_distribution.pdf");
 
}