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: 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):

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

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

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

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

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

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>

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>