rs102_hypotestwithshapes.C File Reference

Detailed Description

rs102_hypotestwithshapes for RooStats project

This tutorial macro shows a typical search for a new particle by studying an invariant mass distribution. The macro creates a simple signal model and two background models, which are added to a RooWorkspace. The macro creates a toy dataset, and then uses a RooStats ProfileLikleihoodCalculator to do a hypothesis test of the background-only and signal+background hypotheses. In this example, shape uncertainties are not taken into account, but normalization uncertainties are.

Processing /mnt/build/workspace/root-makedoc-v612/rootspi/rdoc/src/v6-12-00-patches/tutorials/roostats/rs102_hypotestwithshapes.C...

RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby 
                Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
                All rights reserved, please read http://roofit.sourceforge.net/license.txt

[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooAddPdf::model
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooGaussian::sigModel
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::invMass
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::mH
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::sigma1
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooProduct::fsig
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::mu
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::ratioSigEff
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::fsigExpected
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooGaussian::zjjModel
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::mZ
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::sigma1_z
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::fzjj
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooChebychev::qcdModel
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::a0
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::a1
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing RooRealVar::a2
[#1] INFO:ObjectHandling -- RooWorkspace::import(myWS) importing dataset modelData
[#1] INFO:ObjectHandling -- RooWorkSpace::import(myWS) changing name of dataset from  modelData to data
[#1] INFO:Minization -- createNLL: caching constraint set under name CONSTR_OF_PDF_model_FOR_OBS_invMass with 0 entries
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE 
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit / Migrad with strategy 1
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization --  The following expressions have been identified as constant and will be precalculated and cached: (sigModel,zjjModel,qcdModel)
[#1] INFO:Minization -- 
  RooFitResult: minimized FCN value: 717.039, estimated distance to minimum: 8.90226e-10
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                  fzjj    3.1152e-01 +/-  5.03e-02
                    mu    1.0968e+00 +/-  3.03e-01

[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::GetHypoTest - do conditional fit 
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit / Migrad with strategy 1
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization -- 
  RooFitResult: minimized FCN value: 723.97, estimated distance to minimum: 2.09863e-09
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                  fzjj    2.6213e-01 +/-  5.18e-02

-------------------------------------------------
The p-value for the null is 9.83108e-05
Corresponding to a significance of 3.72332
-------------------------------------------------


[#1] INFO:Minization -- createNLL picked up cached consraints from workspace with 0 entries
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization --  The following expressions have been identified as constant and will be precalculated and cached: (sigModel,zjjModel,qcdModel)
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (sigModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (zjjModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (qcdModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Minization -- createNLL picked up cached consraints from workspace with 0 entries
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization --  The following expressions have been identified as constant and will be precalculated and cached: (sigModel,zjjModel,qcdModel)
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (zjjModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (qcdModel)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()

#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooAddition.h"
#include "RooProduct.h"
#include "TCanvas.h"
#include "RooChebychev.h"
#include "RooAbsPdf.h"
#include "RooFit.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "RooAbsArg.h"
#include "RooWorkspace.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/HypoTestResult.h"
#include <string>


// use this order for safety on library loading
using namespace RooFit;
using namespace RooStats;

// see below for implementation
void AddModel(RooWorkspace*);
void AddData(RooWorkspace*);
void DoHypothesisTest(RooWorkspace*);
void MakePlots(RooWorkspace*);

//____________________________________
void rs102_hypotestwithshapes() {

   // The main macro.

   // Create a workspace to manage the project.
   RooWorkspace* wspace = new RooWorkspace("myWS");

   // add the signal and background models to the workspace
   AddModel(wspace);

   // add some toy data to the workspace
   AddData(wspace);

   // inspect the workspace if you wish
   //  wspace->Print();

   // do the hypothesis test
   DoHypothesisTest(wspace);

   // make some plots
   MakePlots(wspace);

   // cleanup
   delete wspace;
}

//____________________________________
void AddModel(RooWorkspace* wks){

   // Make models for signal (Higgs) and background (Z+jets and QCD)
   // In real life, this part requires an intelligent modeling
   // of signal and background -- this is only an example.

   // set range of observable
   Double_t lowRange = 60, highRange = 200;

   // make a RooRealVar for the observable
   RooRealVar invMass("invMass", "M_{inv}", lowRange, highRange,"GeV");


   // --------------------------------------
   // make a simple signal model.
   RooRealVar mH("mH","Higgs Mass",130,90,160) ;
   RooRealVar sigma1("sigma1","Width of Gaussian",12.,2,100)  ;
   RooGaussian sigModel("sigModel", "Signal Model", invMass, mH, sigma1);
   // we will test this specific mass point for the signal
   mH.setConstant();
   // and we assume we know the mass resolution
   sigma1.setConstant();

   // --------------------------------------
   // make zjj model.  Just like signal model
   RooRealVar mZ("mZ", "Z Mass", 91.2, 0, 100);
   RooRealVar sigma1_z("sigma1_z","Width of Gaussian",10.,6,100)  ;
   RooGaussian zjjModel("zjjModel", "Z+jets Model", invMass, mZ, sigma1_z);
   // we know Z mass
   mZ.setConstant();
   // assume we know resolution too
   sigma1_z.setConstant();


   // --------------------------------------
   // make QCD model
   RooRealVar a0("a0","a0",0.26,-1,1) ;
   RooRealVar a1("a1","a1",-0.17596,-1,1) ;
   RooRealVar a2("a2","a2",0.018437,-1,1) ;
   RooRealVar a3("a3","a3",0.02,-1,1) ;
   RooChebychev qcdModel("qcdModel","A  Polynomial for QCD",invMass,RooArgList(a0,a1,a2)) ;

   // let's assume this shape is known, but the normalization is not
   a0.setConstant();
   a1.setConstant();
   a2.setConstant();

   // --------------------------------------
   // combined model

   // Setting the fraction of Zjj to be 40% for initial guess.
   RooRealVar fzjj("fzjj","fraction of zjj background events",.4,0.,1) ;

   // Set the expected fraction of signal to 20%.
   RooRealVar fsigExpected("fsigExpected","expected fraction of signal events",.2,0.,1) ;
   fsigExpected.setConstant(); // use mu as main parameter, so fix this.

   // Introduce mu: the signal strength in units of the expectation.
   // eg. mu = 1 is the SM, mu = 0 is no signal, mu=2 is 2x the SM
   RooRealVar mu("mu","signal strength in units of SM expectation",1,0.,2) ;

   // Introduce ratio of signal efficiency to nominal signal efficiency.
   // This is useful if you want to do limits on cross section.
   RooRealVar ratioSigEff("ratioSigEff","ratio of signal efficiency to nominal signal efficiency",1. ,0.,2) ;
   ratioSigEff.setConstant(kTRUE);

   // finally the signal fraction is the product of the terms above.
   RooProduct fsig("fsig","fraction of signal events",RooArgSet(mu,ratioSigEff,fsigExpected)) ;

   // full model
   RooAddPdf model("model","sig+zjj+qcd background shapes",RooArgList(sigModel,zjjModel, qcdModel),RooArgList(fsig,fzjj)) ;

   // interesting for debugging and visualizing the model
   //  model.printCompactTree("","fullModel.txt");
   //  model.graphVizTree("fullModel.dot");

   wks->import(model);
}

//____________________________________
void AddData(RooWorkspace* wks){
   // Add a toy dataset

   Int_t nEvents = 150;
   RooAbsPdf* model = wks->pdf("model");
   RooRealVar* invMass = wks->var("invMass");

   RooDataSet* data = model->generate(*invMass,nEvents);

   wks->import(*data, Rename("data"));

}

//____________________________________
void DoHypothesisTest(RooWorkspace* wks){

   // Use a RooStats ProfileLikleihoodCalculator to do the hypothesis test.
   ModelConfig model;
   model.SetWorkspace(*wks);
   model.SetPdf("model");

   //plc.SetData("data");

   ProfileLikelihoodCalculator plc;
   plc.SetData( *(wks->data("data") ));

   // here we explicitly set the value of the parameters for the null.
   // We want no signal contribution, eg. mu = 0
   RooRealVar* mu = wks->var("mu");
   //   RooArgSet* nullParams = new RooArgSet("nullParams");
   //   nullParams->addClone(*mu);
   RooArgSet poi(*mu);
   RooArgSet * nullParams = (RooArgSet*) poi.snapshot();
   nullParams->setRealValue("mu",0);


   //plc.SetNullParameters(*nullParams);
   plc.SetModel(model);
   // NOTE: using snapshot will import nullparams
   // in the WS and merge with existing "mu"
   // model.SetSnapshot(*nullParams);

   //use instead setNuisanceParameters
   plc.SetNullParameters( *nullParams);



   // We get a HypoTestResult out of the calculator, and we can query it.
   HypoTestResult* htr = plc.GetHypoTest();
   cout << "-------------------------------------------------" << endl;
   cout << "The p-value for the null is " << htr->NullPValue() << endl;
   cout << "Corresponding to a significance of " << htr->Significance() << endl;
   cout << "-------------------------------------------------\n\n" << endl;


}

//____________________________________
void MakePlots(RooWorkspace* wks) {

   // Make plots of the data and the best fit model in two cases:
   // first the signal+background case
   // second the background-only case.

   // get some things out of workspace
   RooAbsPdf* model = wks->pdf("model");
   RooAbsPdf* sigModel = wks->pdf("sigModel");
   RooAbsPdf* zjjModel = wks->pdf("zjjModel");
   RooAbsPdf* qcdModel = wks->pdf("qcdModel");

   RooRealVar* mu = wks->var("mu");
   RooRealVar* invMass = wks->var("invMass");
   RooAbsData* data = wks->data("data");


   // --------------------------------------
   // Make plots for the Alternate hypothesis, eg. let mu float

   mu->setConstant(kFALSE);

   model->fitTo(*data,Save(kTRUE),Minos(kFALSE), Hesse(kFALSE),PrintLevel(-1));

   //plot sig candidates, full model, and individual components
   new TCanvas();
   RooPlot* frame = invMass->frame() ;
   data->plotOn(frame ) ;
   model->plotOn(frame) ;
   model->plotOn(frame,Components(*sigModel),LineStyle(kDashed), LineColor(kRed)) ;
   model->plotOn(frame,Components(*zjjModel),LineStyle(kDashed),LineColor(kBlack)) ;
   model->plotOn(frame,Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;

   frame->SetTitle("An example fit to the signal + background model");
   frame->Draw() ;
   //  cdata->SaveAs("alternateFit.gif");

   // --------------------------------------
   // Do Fit to the Null hypothesis.  Eg. fix mu=0

   mu->setVal(0); // set signal fraction to 0
   mu->setConstant(kTRUE); // set constant

   model->fitTo(*data, Save(kTRUE), Minos(kFALSE), Hesse(kFALSE),PrintLevel(-1));

   // plot signal candidates with background model and components
   new TCanvas();
   RooPlot* xframe2 = invMass->frame() ;
   data->plotOn(xframe2, DataError(RooAbsData::SumW2)) ;
   model->plotOn(xframe2) ;
   model->plotOn(xframe2, Components(*zjjModel),LineStyle(kDashed),LineColor(kBlack)) ;
   model->plotOn(xframe2, Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;

   xframe2->SetTitle("An example fit to the background-only model");
   xframe2->Draw() ;
   //  cbkgonly->SaveAs("nullFit.gif");

}

Author

Kyle Cranmer

Definition in file rs102_hypotestwithshapes.C.