-- WilliamBreadenMadden - 2011-08-16

Higgs analysis at ATLAS using RooStats

* THIS PAGE IS UNDER CONSTRUCTION *

This page contains basic information on getting started with Higgs analysis at ATLAS using RooStats.

A note on code and formatting

Example code is generally indented and given highlighting. Explanatory code and text directly relating to it has yellow highlighting while example code that one might run directly in ROOT is given a classical green-text-on-black-background console style. Scripts are given grey highlighting. File contents are given the golden verbatim colouring. So, in a nutshell, code segments are in yellow, full ROOT scripts are in terminal, general scripts are in grey and file contents are in gold. In code examples given, user-created objects are generally prefaced with "my", for example, "myData", for the purposes of clarity.

What is RooStats?

RooStats is a project to create statistical tools built on top of the RooFit library, which is a data modelling toolkit. It is distributed in ROOT. Specifically, it has been distributed in the ROOT release since version 5.22 (December 2008). The latest version of ROOT is recommended as the RooStats project is developing quickly.

Using the appropriate version of ROOT at Glasgow

There are instructions on how to use the different versions of ROOT at Glasgow here.

Setting up RooStats

There are three main options available for acquiring ROOT with RooStats included.

Option 1: Download the latest ROOT release binaries.

The latest ROOT binaries for various operating systems are accessible here.

Option 2: Build the ROOT trunk from source.

Follow the appropriate instructions here to build the ROOT trunk.

Shell script: building ROOT with RooFit and RooStats

#!/bin/bash # This script builds the latest version of ROOT with RooFit and RooStats in Ubuntu. # First, the ROOT prerequisites are installed, # then, the most common ROOT optional packages are installed. # Next, the latest version of ROOT in the CERN Subversion repository is checked out. # Finally, ROOT is compiled. # Install the ROOT prerequisites. sudo apt-get install subversion sudo apt-get install make sudo apt-get install g++ sudo apt-get install gcc sudo apt-get install binutils sudo apt-get install libx11-dev sudo apt-get install libxpm-dev sudo apt-get install libxft-dev sudo apt-get install libxext-dev # Install the optional ROOT packages. sudo apt-get install gfortran sudo apt-get install ncurses-dev sudo apt-get install libpcre3-dev sudo apt-get install xlibmesa-glu-dev sudo apt-get install libglew1.5-dev sudo apt-get install libftgl-dev sudo apt-get install libmysqlclient-dev sudo apt-get install libfftw3-dev sudo apt-get install cfitsio-dev sudo apt-get install graphviz-dev sudo apt-get install libavahi-compat-libdnssd-dev sudo apt-get install libldap-dev sudo apt-get install python-dev sudo apt-get install libxml2-dev sudo apt-get install libssl-dev sudo apt-get install libgsl0-dev # Check out latest ROOT trunk. svn co http://root.cern.ch/svn/root/trunk ~/root # The configuration for the build is set. cd ~/root # Run this to define the system architecture and to enable building of the libRooFit advanced fitting package: ./configure linuxx8664gcc --enable-roofit # See other possible configurations using the following command: ./configure --help # Start compiling. make # Upon completion, ROOT is run by executing ~/root/bin/root. # The following line could be added to the ~/.bashrc file: # export PATH=$PATH:/home/wbm/root/bin

Option 3: Build the RooStats branch.

Do this if you want the latest development of RooStats (that has not yet been incorporated into a ROOT version).

The necessary instructions can be found here.

RooFit

General description

The RooFit library provides a toolkit for modelling the expected distribution of events in a physics analysis. Models can be used to perform unbinned maximum likelihood fits, produce plots and generate "toy Monte Carlo" samples for various studies.

The core functionality of RooFit is to enable the modelling of 'event data' distributions, where each event is a discrete occurrence in time and has one or more measured observables associated with it. Experiments of this nature result in datasets obeying Poisson (or binomial) statistics. The natural modeling language for such distributions are probability density functions (PDFs), F(x;p), that describe the probability density of the distribution of observables x in terms of the function parameter p.

In RooFit, every variable, data point, function and PDF is represented in a C++ object. So, for example, in constructing a RooFit model, the mathematical components of the model map to separate C++ objects. Objects are classified by the data or function type that they represent, not by their respective role in a particular setup. All objects are self-documenting. The name of an object is a unique identified for the object while the title of an object is a more elaborate description of the object.

Here are a few examples of mathematical concepts that correspond to various RooFit classes:

Mathematical concept RooFit class
variable RooRealVar
function RooAbsReal
PDF RooAbsPdf
space point (set of parameters) RooArgSet
list of space points RooAbsData
integral RooRealIntegral

Composite functions correspond to composite objects. The ArgSet class is dependent on argument order while the ArgList class is not.

Example code: defining a RooFit variable

General form for defining a RooFit variable: RooRealVar x(<object name>, <object title>, <value>, <minimum value>, <maximum value>) Specific example for defining a RooFit variable x with the value 5: RooRealVar x("x", "x observable", 5, -10, 10)

RooPlot

A RooPlot is essentially an empty frame that is capable of holding anything plotted verses its variable.

PDFs

One of the things that makes the RooFit PDFs nice and flexible (but perhaps counterintuitive) is that they have no idea what is considered the observable. For example, a Gaussian has x, the mean and sigma while the PDF is always normalised to unity. What is one integrating in order to get 1? For the Gaussian, it is x, by convention; the mean and sigma are parameters of the model. RooFit doesn't know that x is special; x, the mean and sigma are all on equal footing. You can tell RooFit that x is the variable to normalise over.

So, RooGaussian has no intrinsic notion of distinction between observables and parameters. The choice of observables (for unit normalisation) is always passed to gauss.getVal().

What is the value of this? In a nutshell, it allows one to do Bayesian stuff very easily.

Bayes' theorem holds that the probability of b given a is related to the probability of a given b. Normally, one might say that the mean of a Gaussian is 1 and that then gives a distribution for x. However, if one had a dataset for x, one might want to know the posterior for the mean. The RooFit ability to switch what the probability density function is of is very useful for Bayesian stuff. Because RooFit allows this composition, the thing that acts as x in the Gaussian could actually be a function of, say, x^2. A Jacobian factor is picked up in going from x to x^2, so the normalisation no longer makes sense. Sometimes RooFit can figure out what the Jacobian factor should be and sometimes it resorts to numerical integration.

Example code: create a Gaussian PDF using RooStats and plot it using the RooPlot class

{ // Build a Gaussian PDF. RooRealVar x("x", "x", -10, 10); RooRealVar mean("mean", "mean of Gaussian", 0, -10, 10); RooRealVar sigma("sigma", "width of Gaussian", 3); RooGaussian gauss("gauss", "Gaussian PDF", x, mean, sigma); // Plot the PDF. RooPlot* xframe = x.frame(); gauss.plotOn(xframe); xframe->Draw(); }

Example code: telling a RooFit PDF what to normalise over

Not normalised (i.e., this is not a PDF): gauss.getVal(); Hey, RooFit! This is the thing to normalise over (i.e., guarantees Int[xmin, xmax] Gauss(x, m, s)dx == 1): gauss.getVal(x); What is the value if sigma is considered the observable? (i.e., guarantees Int[smin, smax] Gauss(x, m, s)ds == 1):

The RooFit workspace

General description

The RooFit "workspace" provides the ability to store the full likelihood model, any desired priors and the minimal data necessary to reproduce the likelihood function in a ROOT file. Thus, the workspace is needed for combinations and has potential for digitally publishing results (PhyStats agreed to publish likelihood functions). The RooFit workspace can be used for this.

One might create a Gaussian PDF and then import it into a workspace. All of the dependencies of the Gaussian shall be drawn in and owned by the workspace. There are no nightmarish ownership problems. Alternatively, one might create simply the Gaussian inside the workspace using the "Workspace Factory".

Example code: using the Workspace Factory to create a Gaussian PDF

// Create a Gaussian PDF using the Workspace Factory (this is essentially the shorthand for creating a Gaussian). RooWorkspace* myWorkspace = new RooWorkspace("myWorkspace"); myWorkspace->factory("Gaussian::g(x[-5, 5], mu[0], sigma[1]");

What's in the RooFit workspace?

Example code: What's in the workspace?

// Open the appropriate ROOT file. root -l myFile.root // Import the workspace. myWorkspace = (RooWorkspace*) _file0->Get("myWorkspace"); // Print the workspace contents. myWorkspace.Print(); // Example printout: variables --------- (x,m,s) p.d.f.s ------- RooGaussian::g[ x=x mean=m sigma=s ] = 0 datasets -------- RooDataSet::d(x) // Import the variable saved as x. RooRealVar* myVariable = myWorkspace.var("x"); // Import the PDF saved as g. RooAbsPdf* myPDF = myWorkspace.pdf("g"); // Import the data saved as d. RooAbsData* myData = myWorkspace.data("d"); // Import the ModelConfig saved as m. ModelConfig* myModelConfig = (ModelConfig*) myWorkspace.obj("m");

Visual representations of the model/PDF contents

Graphviz

Graphviz consists of a graph description language called the DOT language and a set of tools that can generate and/or process DOT files.

Example code: examining PDFs and creating graphical representations of them

// Create variables and a PDF using those variables. RooRealVar mu("mu", "mu", 150); RooRealVar sigma("sigma", "sigma", 5, 0, 20); RooGaussian myGaussianPDF("myGaussianPDF", "Gaussian PDF", x, mu, sigma); // Create a Graphviz DOT file with a representation of the object tree. myGaussianPDF.graphVizTree("myGaussianPDFTree.dot"); // This produced DOT file can be converted to some graphical representation: // Convert the DOT file to a 'top-to-bottom graph' using UNIX commands: // dot -Tgif -o myGaussianPDF_top-to-bottom_graph.gif myGaussianPDFTree.dot // Convert the DOT file to a 'spring-model graph' using UNIX commands: // fdp -Tgif -o myGaussianPDF_spring-model_graph.gif myGaussianPDFTree.dot // Print the PDF contents. myGaussianPDF.Print(); // Example output: // RooGaussian::G[ x=x mean=mu sigma=sigma ] = 1 // sigma.Print() // RooRealVar::sigma = 5 L(0 - 20) // Print the PDF contents in a detailed manner. myGaussianPDF.Print("verbose"); // Print the PDF contents to stdout. myGaussianPDF.Print("t"); // Example output: // 0x166eab0 RooGaussian::G = 1 [Auto] // 0x15f7fe0/V- RooRealVar::x = 150 // 0x1487090/V- RooRealVar::mu = 150 // 0x1487bc0/V- RooRealVar::sigma = 5 // Print the PDF contents to a file. myGaussianPDF.printCompactTree("", "myGaussianPDFTree.txt") // Example output file contents: // 0x166eab0 RooGaussian::G = 1 [Auto] // 0x15f7fe0/V- RooRealVar::x = 150 // 0x1487090/V- RooRealVar::mu = 150 // 0x1487bc0/V- RooRealVar::sigma = 5

Accessing the RooFit workspace

Example code: accessing the workspace

// Open the appropriate ROOT file. root -l BR5_MSSM_signal90_combined_datastat_model.root // Alternatively, you could open the file in a manner such as the following: myFileName = "BR5_MSSM_signal90_combined_datastat_model.root" TFile *myFile = TFile::Open(myFileName); // Import the workspace. RooWorkspace* myWorkspace = (RooWorkspace*) _file0->Get("combined"); // Print the workspace contents. myWorkspace->Print(); // Import the PDF. RooAbsPdf* myPDF = myWorkspace->pdf("model_BR5_MSSM_signal90"); // Import the variable representing the observable. RooRealVar* myObservable = myWorkspace->var("obs"); // Create a RooPlot frame using the imported variable.. RooPlot* myFrame = myObservable.frame(); // Plot the PDF on the created RooPlot frame. myPDF.plotOn(myFrame); // Draw the RooPlot. myFrame->Draw();

Example code: accessing both data and PDF from a workspace stored in a file

// Note that the following code is independent of actual PDF in the file. So, for example, a full Higgs combination could work with identical code. // Open a file and import the workspace. TFile myFile("myResults.root") ; RooWorkspace* myWorkspace = f.Get("myWorkspace") ; // Plot the data and PDF RooPlot* xframe = w->var("x")->frame() ; w->data("d")->plotOn(xframe) ; w->pdf("g")->plotOn(xframe) ; // Construct a likelihood and profile likelihood RooNLLVar nll("nll","nll",*myWorkspace->pdf("g"),*w->data("d")) ; RooProfileLL pll("pll","pll", nll,*myWorkspace->var("m")) ; RooPlot* myFrame = w->var("m")->frame(-1,1) ; pll.plotOn(myFrame) ; myFrame->Draw()

RooStats


General description

RooStats provides tools for high-level statistics questions in ROOT. It builds on RooFit, which provides basic building blocks for statistical questions.

The main goal is to standardise the interface for major statistical procedures so that they can work on an arbitrary RooFit model and dataset and handle any parameters of interest and nuisance parameters. Another goal is to implement most accepted techniques from frequentist, Bayesian and likelihood based approaches. A further goal is to provide utilities to perform combined measurements.

Datasets

In RooFit, data can be stored in an unbinned or binned manner. Unbinned data is stored using the RooDataSet class while binned data is stored using the RooDataHist class. The user can define how many bins there are in a variable. For the purposes of plotting using the RooPlot class, a RooDataSet is binned into a histogram.

Example code: generating toy Monte Carlo, storing it as unbinned data and then plotting it

} // Create a RooDataSet and fill it with generated toy Monte Carlo data: RooDataSet* myData = gauss.generate(x, 100); // Plot the dataset. RooPlot* myFrame = x.frame(); myData.plotOn()(myFrame); myFrame.Draw(); }

Importing data (how to populate datasets from histograms and TTrees)

Importing data from ROOT trees

* Put stuff here, Will. *

Importing data from ROOT TH histogram objects (take a histogram and map it to a binned data set)

Example code: import a ROOT histogram into a RooDataHist (a RooFit binned dataset)

{ // Access the file. TFile* myFile = new TFile("myFile.root"); // Load the histogram. TH1* myHistogram = myFile.Get("myHistogram"); // Draw the loaded histogram. myHistogram.Draw(); // Declare an observable x. RooRealVar x("x", "x", -1, 2); // Create a binned dataset that imports the contents of TH1 and associates its contents to observable 'x'. RooDataHist myData("myData", "myData", RooArgList(x), myHistogram); // Plot the imported dataset. RooPlot* myFrame = x.frame(); myData.plotOn(myFrame) myFrame.Draw() }

Fitting

Fitting a PDF to unbinned data

Example code: fit a Gaussian PDF to data

// Fit gauss to unbinned data gauss.fitTo(*myData);

Example code: create a simple model using the RooFit Workspace Factory. Specify parts of the model using ModelConfig. Create a simple dataset. Complete a confidence interval test using the ProfileLikelihoodCalculator of RooStats

{ // In this script, a simple model is created using the Workspace Factory in RooFit. // ModelConfig is used to specify the parts of the model necessary for the statistical tools of RooStats. // A 95% confidence interval test is run using the ProfileLikelihoodCalculator of RooStats. // Define a RooFit random seed in order to produce reproducible results. RooRandom::randomGenerator()->SetSeed(271); // Make a simple model using the Workspace Factory. // Create a new workspace. RooWorkspace* myWorkspace = new RooWorkspace(); // Create the PDF G(x|mu,1) and the variables x, mu and sigma in one command using the factory syntax. myWorkspace->factory("Gaussian::normal(x[-10,10], mu[-1,1], sigma[1])"); // Define parameter sets for observables and parameters of interest. myWorkspace->defineSet("poi","mu"); myWorkspace->defineSet("obs","x"); // Print the workspace contents. myWorkspace->Print() ; // Specify for the statistical tools the components of the defined model. // Create a new ModelConfig. ModelConfig* myModelConfig = new ModelConfig("my G(x|mu,1)"); // Specify the workspace. myModelConfig->SetWorkspace(*myWorkspace); // Specify the PDF. myModelConfig->SetPdf(*myWorkspace->pdf("normal")); // Specify the parameters of interest. myModelConfig->SetParametersOfInterest(*myWorkspace->set("poi")); // Specify the observables. myModelConfig->SetObservables(*myWorkspace->set("obs")); // Create a toy dataset. // Create a toy dataset of 100 measurements of the observables (x). RooDataSet* myData = myWorkspace->pdf("normal")->generate(*myWorkspace->set("obs"), 100); //myData->print(); // Use the ProfileLikelihoodCalculator to obtain a 95% confidence interval. // Specify the confidence level required. double confidenceLevel = 0.95; // Create an instance of the ProfileLikelihoodCalculator, specifying the data and the ModelConfig for it. ProfileLikelihoodCalculator myProfileLikelihoodCalculator(*myData, *myModelConfig); // Set the confidence level. myProfileLikelihoodCalculator.SetConfidenceLevel(confidenceLevel); // Obtain the resulting interval. LikelihoodInterval* myProfileLikelihoodInterval = myProfileLikelihoodCalculator.GetInterval(); // Use this interval result. In this case, it makes sense to say what the lower and upper limits are. // Define the object variables for the purposes of the confidence interval. RooRealVar* x = myWorkspace->var("x"); RooRealVar* mu = myWorkspace->var("mu"); cout << "The profile likelihood calculator interval is [ "<< myProfileLikelihoodInterval->LowerLimit(*mu) << ", " << myProfileLikelihoodInterval->UpperLimit(*mu) << "] " << endl; // Set mu equal to zero. mu->setVal(0); // Is mu in the interval? cout << "Is mu = 0 in the interval?" << endl; if (myProfileLikelihoodInterval->IsInInterval(*mu) == 1){ cout << "Yes" << endl; } else{ cout << "No" << endl; } }

ModelConfig

The ModelConfig RooStats class encapsulates the configuration of a model to define a particular hypothesis. It is now used extensively by the calculator tools. ModelConfig always contains a reference to an external workspace that manages all of the objects that are a part of the model (PDFs and parameter sets). So, in order to use ModelConfig, the user must specify a workspace pointer before creating the various objects of the model.

HistFactory

General description

The HistFactory can be used to run an analysis without being a RooFit/RooStats expert. In a nutshell, ROOT files containing input histograms are set up and XML files are set up for those input ROOT files. The XML files specify details on the histograms and specify how RooFit should interpret the information in the files. An little executable contained in ROOT called hist2workspace is used to import the histograms into a RooFit workspace appropriately.

prepareHistFactory

The ROOT release ships with a script called prepareHistFactory and a binary file called hist2workspace in the $ROOTSYS/bin directories. prepareHistFactory prepares a working area. It creates results/, data/ and config/ directories. It copies the HistFactorySchema.dtd and example XML files into the config/ directory. Also, it copies a ROOT file into the data/ directory for use with the examples.

HistFactorySchema.dtd

HistFactorySchema.dtd: This file is located in $ROOTSYS/etc/ specifies the XML schema. It is typically placed in the config/ directory of a working area together with the top-level XML file and the individual channel XML files. The user should not modify this file. The HistFactorySchema.dtd is commented to specify exactly the meaning of the various options.

hist2workspace

The hist2workspace executable is used in the following manner: hist2workspace input.xml

Top-Level XML file

General description

This file specifies a top level 'Combination' that is composed of:

  • several 'Channels', which are described in separate XML files.
  • several 'Measurements' (corresponding to a full fit of the model), each of which specifies
    • a name for this measurement to be used in tables and files.
    • the luminosity associated with the measurement in picobarns.
    • which bins of the histogram should be used.
    • what the relative uncertainty on the luminosity is.
    • what is(/are) the parameter(/s) of interest that will be measured.
    • which parameter(/s) should be fixed/floating (e.g., nuisance parameters)
    • which type of constriants are desired:
      • default: Gaussian
      • supported: Gamma, LogNormal, Uniform
    • whether the tool should export the model only and skip the default fit.

Specific instructions

There are several parts in the top-level XML file. Initially, the XML schema is specified. The output is specified next. Specifically, the prefix for any output files (ROOT files containing workspaces) is given. Now, the channel XML files are given for each measurement (e.g., signal, background, systematics information) using the Input tag. Next, the details for each channel are defined using the Measurement tag. Which bins to use are specified inclusively using the Measurement tag attributes BinLow and BinHigh. Within the Measurement tag, the POI tag is used to specify the point of interest for the channel, i.e. the test statistic.

Example file: $ROOTSYS/tutorials/histfactory/example.xml


<!DOCTYPE Combination  SYSTEM 'HistFactorySchema.dtd'>

<Combination OutputFilePrefix="./results/example" Mode="comb" >

  <Input>./config/example_channel.xml</Input>

  <Measurement Name="GaussExample" Lumi="1." LumiRelErr="0.1" BinLow="0" BinHigh="2" Mode="comb" >
    <POI>SigXsecOverSM</POI>
    <ParamSetting Const="True">Lumi alpha_syst1</ParamSetting>
    <!-- don't need <ConstraintTerm> default is Gaussian-->
  </Measurement>

  <Measurement Name="GammaExample" Lumi="1." LumiRelErr="0.1" BinLow="0" BinHigh="2" Mode="comb" >
    <POI>SigXsecOverSM</POI>
    <ParamSetting Const="True">Lumi alpha_syst1</ParamSetting>
    <ConstraintTerm Type="Gamma" RelativeUncertainty=".3">syst2</ConstraintTerm>
  </Measurement>

  <Measurement Name="LogNormExample" Lumi="1." LumiRelErr="0.1" BinLow="0" BinHigh="2" Mode="comb" >
    <POI>SigXsecOverSM</POI>
    <ParamSetting Const="True">Lumi alpha_syst1</ParamSetting>
    <ConstraintTerm Type="LogNormal" RelativeUncertainty=".3">syst2</ConstraintTerm>
  </Measurement>

  <Measurement Name="ConstExample" Lumi="1." LumiRelErr="0.1" BinLow="0" BinHigh="2" Mode="comb" ExportOnly="True">
    <POI>SigXsecOverSM</POI>
    <ParamSetting Const="True">Lumi alpha_syst1</ParamSetting>
  </Measurement>


</Combination>

Channel XML files

General description

These files specify for each channel

  • observed data (if absent, the tool will use the expectation, which is useful for expected sensitivity)
  • several 'Samples' (e.g., signal, bkg1, bkg2 etc.), each of which specifies
    • a name.
    • whether the sample is normalized by theory (e.g., N = L * sigma) or whether the sample is data driven.
    • a nominal expectation histogram.
    • a named 'Normalization Factor' (which can be fixed or allowed to float in a fit).
    • several 'Overall Systematics' in normalization with
      • a name.
      • +/- 1 sigma variations (e.g., 1.05 and 0.95 for a 5% uncertainty).
    • several 'Histogram Systematics' in shape with
      • a name (which can be shared with the OverallSyst if correlated).
      • +/- 1 sigma variational histograms.

Specific instructions

* text *

Example file: $ROOTSYS/tutorials/histfactory/example_channel.xml


<!DOCTYPE Channel  SYSTEM 'HistFactorySchema.dtd'>

  <Channel Name="channel1" InputFile="./data/example.root" HistoName="" >
    <Data HistoName="data" HistoPath="" />
    <Sample Name="signal" HistoPath="" HistoName="signal">
      <OverallSys Name="syst1" High="1.05" Low="0.95"/>
      <NormFactor Name="SigXsecOverSM" Val="1" Low="0." High="3." Const="True" />
    </Sample>
    <Sample Name="background1" HistoPath="" NormalizeByTheory="True" HistoName="background1">
      <OverallSys Name="syst2" Low="0.95" High="1.05"/>
    </Sample>
    <Sample Name="background2" HistoPath="" NormalizeByTheory="True" HistoName="background2">
      <OverallSys Name="syst3" Low="0.95" High="1.05"/>
      <!-- <HistoSys Name="syst4" HistoPathHigh="" HistoPathLow="histForSyst4"/>-->
    </Sample>
  </Channel>

Further information

ROOT links:

ROOT User's Guide

RooFit links

User's Manual

Tutorials

RooStats links

Wiki

RooStats User's Guide

Tutorials

HistFactory manual

E-group for support with ATLAS-sensitive information: atlas-phys-stat-root@cern.ch

E-mail support for software issues, bugs etc.: roostats-development@cern.ch

-- WilliamBreadenMadden - 2010-10-29

Edit | Attach | Print version | History: r51 | r24 < r23 < r22 < r21 | Backlinks | Raw View | Raw edit | More topic actions...
Topic revision: r22 - 2011-08-16 - WilliamBreadenMadden
 
This site is powered by the TWiki collaboration platform Powered by PerlCopyright © 2008-2020 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback