Difference: HiggsAnalysisAtATLASUsingRooStats (50 vs. 51)

Revision 512016-05-06 - WilliamBreadenMadden

Line: 1 to 1
 
META TOPICPARENT name="WebHome"
Deleted:
<
<		-- Will Breaden Madden - 2016-01-25
 

Higgs analysis at ATLAS using RooStats

Changed:
<
<		This page contains basic information on getting started with Higgs analysis at ATLAS using RooStats.
>
>		This page contains basic information on getting started with Higgs analysis at ATLAS using RooStats. Only the most meager of attempts is made to keep this documentation current.
 
Line: 14 to 12
  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 develops quickly.
Deleted:
<
<

using ROOT on the Glasgow PPELX network

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

Execute commands something like the following in order to set up ROOT version 5.32.00 on PPELX:

%CODE{"bash"}% export ROOTSYS=/data/ppe01/sl5x/x86_64/root/5.32.00 export PATH=$ROOTSYS/bin:$PATH export LD_LIBRARY_PATH=$ROOTSYS/lib/root:$LD_LIBRARY_PATH %ENDCODE%

using ROOT on the CERN LXPLUS network

Execute commands something like the following in order to set up ROOT version 5.32.00 on LXPLUS:

%CODE{"bash"}% . /afs/cern.ch/sw/lcg/external/gcc/4.3.2/x86_64-slc5/setup.sh cd /afs/cern.ch/sw/lcg/app/releases/ROOT/5.32.00/x86_64-slc5-gcc43-opt/root/ . bin/thisroot.sh %ENDCODE%

 

setting up RooStats

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

Line: 208 to 184
  %CODE{"c++"}% // general form for defining a RooFit variable:
Changed:
<
<		RooRealVar x(<object name>, <object title>, <value>, <minimum value>, <maximum value>)
>
>		RooRealVar x(, , , , )
 // specific example for defining a RooFit variable x with the value 5: RooRealVar x("x", "x observable", 5, -10, 10) %ENDCODE%
Line: 225 to 201
  What is the value of this? In a nutshell, it allows one to do Bayesian stuff very easily.
Changed:
<
<		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.
>
>		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

Line: 240 to 216
  // Plot the PDF. RooPlot* xframe = x.frame(); gauss.plotOn(xframe);
Changed:
<
<		 xframe->Draw();
>
>		 xframe->Draw();
 } %ENDCODE%
Line: 284 to 260
 %CODE{"c++"}% RooPlot* myFrame = x.frame() ; myData.plotOn(myFrame, Binning(25)) ;
Changed:
<
<		myFrame->Draw()
>
>		myFrame->Draw()
 %ENDCODE%

importing data from ROOT trees (how to populate RooDataSets from TTrees)
Line: 306 to 282
  // Access the file. TFile* myFile = new TFile("myFile.root"); // Load the histogram.
Changed:
<
<		 TH1* myHistogram = (TH1*) myFile->Get("myHistogram");
>
>		 TH1* myHistogram = (TH1*) myFile->Get("myHistogram");
  // Draw the loaded histogram. myHistogram.Draw();
Line: 351 to 327
 %CODE{"c++"}% // Create a Gaussian PDF using the Workspace Factory (this is essentially the shorthand for creating a Gaussian). RooWorkspace* myWorkspace = new RooWorkspace("myWorkspace");
Changed:
<
<		 myWorkspace->factory("Gaussian::g(x[-5, 5], mu[0], sigma[1]");
>
>		 myWorkspace->factory("Gaussian::g(x[-5, 5], mu[0], sigma[1]");
 %ENDCODE%

What's in the RooFit workspace?

Line: 362 to 338
 // Open the appropriate ROOT file. root -l myFile.root // Import the workspace.
Changed:
<
<		myWorkspace = (RooWorkspace*) _file0->Get("myWorkspace");
>
>		myWorkspace = (RooWorkspace*) _file0->Get("myWorkspace");
 // Print the workspace contents. myWorkspace.Print(); // Example printout:
Line: 478 to 454
 // myFileName = "BR5_MSSM_signal90_combined_datastat_model.root" // TFile *myFile = TFile::Open(myFileName); // Import the workspace.
Changed:
<
<		RooWorkspace* myWorkspace = (RooWorkspace*) _file0->Get("combined");
>
>		RooWorkspace* myWorkspace = (RooWorkspace*) _file0->Get("combined");
 // Print the workspace contents.
Changed:
<
<		myWorkspace->Print();
>
>		myWorkspace->Print();
 // Import the PDF.
Changed:
<
<		RooAbsPdf* myPDF = myWorkspace->pdf("model_BR5_MSSM_signal90");
>
>		RooAbsPdf* myPDF = myWorkspace->pdf("model_BR5_MSSM_signal90");
 // Import the variable representing the observable.
Changed:
<
<		RooRealVar* myObservable = myWorkspace->var("obs");
>
>		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.
Changed:
<
<		myFrame->Draw();
>
>		myFrame->Draw();
 %ENDCODE%

example code: accessing both data and PDF from a workspace stored in a file
Line: 502 to 478
 TFile myFile("myResults.root") ; RooWorkspace* myWorkspace = f.Get("myWorkspace") ; // Plot the data and PDF
Changed:
<
<		RooPlot* xframe = w->var("x")->frame() ; w->data("d")->plotOn(xframe) ; w->pdf("g")->plotOn(xframe) ;
>
>		RooPlot* xframe = w->var("x")->frame() ; w->data("d")->plotOn(xframe) ; w->pdf("g")->plotOn(xframe) ;
 // Construct a likelihood and profile likelihood
Changed:
<
<		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) ;
>
>		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) ;
Changed:
<
<		myFrame->Draw()
>
>		myFrame->Draw()
 %ENDCODE%

links for RooFit

Line: 535 to 511
  // A 95% confidence interval test is run using the ProfileLikelihoodCalculator of RooStats.

// Define a RooFit random seed in order to produce reproducible results.

Changed:
<
<		 RooRandom::randomGenerator()->SetSeed(271);
>
>		 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.

Changed:
<
<		 myWorkspace->factory("Gaussian::normal(x[-10,10], mu[-1,1], sigma[1])");
>
>		 myWorkspace->factory("Gaussian::normal(x[-10,10], mu[-1,1], sigma[1])");
  // Define parameter sets for observables and parameters of interest.
Changed:
<
<		 myWorkspace->defineSet("poi","mu"); myWorkspace->defineSet("obs","x");
>
>		 myWorkspace->defineSet("poi","mu"); myWorkspace->defineSet("obs","x");
  // Print the workspace contents.
Changed:
<
<		 myWorkspace->Print() ;
>
>		 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.
Changed:
<
<		 myModelConfig->SetWorkspace(*myWorkspace);
>
>		 myModelConfig->SetWorkspace(*myWorkspace);
  // Specify the PDF.
Changed:
<
<		 myModelConfig->SetPdf(*myWorkspace->pdf("normal"));
>
>		 myModelConfig->SetPdf(*myWorkspace->pdf("normal"));
  // Specify the parameters of interest.
Changed:
<
<		 myModelConfig->SetParametersOfInterest(*myWorkspace->set("poi"));
>
>		 myModelConfig->SetParametersOfInterest(*myWorkspace->set("poi"));
  // Specify the observables.
Changed:
<
<		 myModelConfig->SetObservables(*myWorkspace->set("obs"));
>
>		 myModelConfig->SetObservables(*myWorkspace->set("obs"));
  // Create a toy dataset.

// Create a toy dataset of 100 measurements of the observables (x).

Changed:
<
<		 RooDataSet* myData = myWorkspace->pdf("normal")->generate(*myWorkspace->set("obs"), 100); //myData->print();
>
>		 RooDataSet* myData = myWorkspace->pdf("normal")->generate(*myWorkspace->set("obs"), 100); //myData->print();
  // Use the ProfileLikelihoodCalculator to obtain a 95% confidence interval.
Line: 580 to 556
  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.
Changed:
<
<		 RooRealVar* x = myWorkspace->var("x"); RooRealVar* mu = myWorkspace->var("mu"); cout << "The profile likelihood calculator interval is [ "<< myProfileLikelihoodInterval->LowerLimit(*mu) << ", " << myProfileLikelihoodInterval->UpperLimit(*mu) << "] " << endl;
>
>		 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.
Changed:
<
<		 mu->setVal(0);
>
>		 mu->setVal(0);
  // Is mu in the interval?
Changed:
<
<		 cout << "Is mu = 0 in the interval?" << endl; if (myProfileLikelihoodInterval->IsInInterval(*mu) == 1){ cout << "Yes" << endl;
>
>		 cout << "Is mu = 0 in the interval?" << endl; if (myProfileLikelihoodInterval->IsInInterval(*mu) == 1){ cout << "Yes" << endl;
  } else{
Changed:
<
<		 cout << "No" << endl;
>
>		 cout << "No" << endl;
  } } %ENDCODE%
Line: 996 to 972
 
example (early in project development): creation of the Measurement and a Channel and, thence, creation of the channel samples, including signal and backgrounds

%CODE{"c++"}%

Changed:
<
<		// Create the Measurement and a channel
>
>		// Create the Measurement and a channel.
  std::string myInputFile = "./data/myData.root"; std::string myChannel1Path = "";
Line: 1127 to 1103
 
details on the object measurement
Changed:
<
<		A measurement has several methods to configure its options, each of which are equivalent to their XML equivalents.
>
>		A measurement has several methods to configure its options, each of which is equivalent to their XML equivalents.
 
objective code
set the prefix for output files void SetOutputFilePrefix(const std::string& prefix);
Line: 1192 to 1168
 void AddShapeSys(std::string Name, Constraint::Type ConstraintType, std::string HistoName, std::string HistoFile, std::string HistoPath=""); %ENDCODE%
Changed:
<
<		A sample can be included in a channel's bin-by-bin statistical uncertainty fluctuations by "activating" the sample. There are two ways to do this. The first way is to use the default errors that are stored in the histogram's uncertainty array. The second way is to supply the errors using an external histogram (in the case where the desired errors differ from those stored by the HT1 histogram). These can be achieved using thw following methods:
>
>		A sample can be included in a channel's bin-by-bin statistical uncertainty fluctuations by "activating" the sample. There are two ways to do this. The first way is to use the default errors that are stored in the histogram's uncertainty array. The second way is to supply the errors using an external histogram (in the case where the desired errors differ from those stored by the HT1 histogram). These can be achieved using the following methods:
  %CODE{"c++"}% void ActivateStatError();
Line: 1228 to 1204
 

links for HistFactory

Changed:
<
<
>
>
 
View topic | History: r51 < r50 < r49 < r48 | More topic actions...
 
This site is powered by the TWiki collaboration platform Powered by PerlCopyright © 2008-2025 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback