Higgs analysis at ATLAS using RooStats

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

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)

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: 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]");

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();
}

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

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