ROOT Tutorials Exercises for GridKa School

• create a TF1 object using a formula expression. In the case of a parametric functions, the two parameters are defined as [0] and [1]
• call TF1::Draw()
• call TF1::Eval(x), TF1::Derivative(x) and TF1::Integral(a,b)

You can also find the available member functions of the class TF1, by using the Tab key on the ROOT prompt, for example

root [0] TF1 * f1 = ....
root [1] f1-><TAB>

root [1] f1->Derivative(<TAB>

#include "TF1.h"

void plotFunction() {

   TF1 * f1 = new TF1("f1","sin(x)/x",0,10);
   f1->Draw();


   TF1 * fp = new TF1("fp","[0]*sin([1]*x)/x",0,10);
   fp->SetParameters(1,2);
   fp->Draw("same");
   fp->SetLineColor(kBlue);

   // to change axis y margins                                                                                                                  
   // (or invert the order of plotting the functions)                                                                                           
   f1->SetMaximum(2);
   f1->SetMinimum(-2);

   std::cout << "Value of f(x) at x = 1 is       " << fp->Eval(1.) << std::endl; 
   std::cout << "Derivative of f(x) at x = 1 is  " << fp->Derivative(1.) << std::endl; 
   std::cout << "Integral of f(x) in [0,3] is    " << fp->Integral(0,4)  << std::endl;
 
}

#include "TF1.h"
#include "Math/DistFunc.h"


void plotGamma() {

   // Note that parameter [0] is called alpha in definition of gamma_pdf or kappa in Wikipedia                                                  
   // Note that parameter [1] is  theta                                                                                                         

   // use range [0,20] as in Wikipedia plot                                                                                                     
   TF1 * f1 = new TF1("f","ROOT::Math::gamma_pdf(x,[0],[1])",0,20);

   f1->SetLineColor(kRed);
   f1->SetParameter(0,1);
   f1->SetParameter(1,2);

   // use DrawClone because we will plot many different copies of same object but with different                                                
   // parameter values                                                                                                                          
   f1->DrawClone();


   // now change parameters and draw at different parameter values                                                                              
   f1->SetLineColor(kGreen);
   f1->SetParameter(0,2);
   f1->DrawClone("SAME");

   f1->SetLineColor(kBlue);
   f1->SetParameter(0,3);
   f1->DrawClone("SAME");

   f1->SetLineColor(kCyan);
   f1->SetParameter(0,5);
   f1->SetParameter(1,1);
   f1->DrawClone("SAME");


   f1->SetLineColor(kOrange);
   f1->SetParameter(0,9);
   f1->SetParameter(1,0.5);
   f1->DrawClone("SAME");


}

#include "TH1.h"
#include "TRandom.h"
#include "TStyle.h"

void plotHistogram() { 

   TH1D * h1 = new TH1D("h1","h1",50,0,10);

   for (int i = 0; i < 10000; ++i) {
      double x = gRandom->Gaus(5,2);
      h1->Fill(x); 
   }

   h1->Draw();

   gStyle->SetOptStat(111111110);

}

• Use TH1::Add to add the two histogram.
• Use TH1::Scale(10000/ 100000) to re-normalise the histogram.
• Use again TH1::Add to subtract the histogram, but with a second coefficient equal to -1. (TH1::Add(h1,h2,1,-1)).
• For getting the right errors you must call TH1::Sumw2 before doing the operations on the histograms (i.e. before scaling and before subtracting them)

#include "TH1.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "TLine.h"
#include <iostream>

void HistoOperations() { 


   TH1D * h1 = new TH1D("h1","flat histogram",100,0,10);
   for (int i = 0; i < 10000; ++i) { 
      h1->Fill( gRandom->Uniform(0,10) ); 
   }

   TH1D * h2 = new TH1D("h2","gaus histogram",100,0,10);
   for (int i = 0; i < 1000; ++i) { 
      h2->Fill( gRandom->Gaus(5,1) ); 
   }

   TH1D * h3 = new TH1D("h3","flat+gaus histogram",100,0,10);
   h3->Add(h1,h2);

   h3->Draw();


   TH1D * h4 = new TH1D("h4","second flat histogram",100,0,10);
   for (int i = 0; i < 100000; ++i) { 
      h4->Fill( gRandom->Uniform(0,10) ); 
   }

   // renormalize histogram
   //!! VERY IMPORTANT - NEED TO CALL Sumw2() to get right ERRORS!!!
   h4->Sumw2();
   h3->Sumw2();
   h4->Scale(0.1);

   // plot in a new Canvas
   new TCanvas("c2","c2");


  TH1D * h5 = new TH1D("h5","(flat+gaus) histogram - flat histogram",100,0,10);
  h5->Add(h3,h4,1.,-1.);
  h5->Draw("E");

  TLine * line = new TLine(0,0,10,0);
  line->Draw();
   

}

new TCanvas("c2","c2");

  TH1 * h6 = h5->Rebin(4,"h6");
  h6->SetTitle("Rebinned histogram");
  h6->Draw();

#include "TH2.h"
#include "TProfile.h"
#include "TRandom.h"
#include "TCanvas.h"
#include <iostream>

void Histogram2D() { 


   TH2D * h2d = new TH2D("h2d","2d histogram",40,-5,5, 40, -5, 5);
   for (int i = 0; i < 100000; ++i) { 
      double u =  gRandom->Gaus(0,1);
      double w =  gRandom->Gaus(0,1);      
      double x = u; 
      double y = w + 0.5 * u;
      h2d->Fill( x,y ); 
   }

   h2d->Draw("COLZ");

   std::cout << "correlation factor " << h2d->GetCorrelationFactor() << std::endl;


}

For dividing the canvas call TCanvas::Divide(1,2) and navigate in the pad contained in the canvas by calling TCanvas::cd(pad_number).

TH1 * hx = h2d->ProjectionX(); 
   TH1 * px = h2d->ProfileX(); 

   TCanvas * c2 = new TCanvas("c2","c2"); 
   // divide in 2 pad in x and one in y
   c2->Divide(2,1);
   c2->cd(1); 
   hx->Draw(); 

   c2->cd(2);
   px->Draw();

#include "TFile.h"
#include "TH1.h"
#include "TRandom.h"

void histogramWrite() { 

   TFile f("histogram.root","RECREATE");

   TH1D * h1 = new TH1D("h1","h1",100,0,10); 
   for (int i = 0; i < 10000; ++i) 
      h1->Fill(gRandom->Exp(5) ); 

   h1->Fit("expo");
   h1->Draw();

   f.Write("h1"); 
   f.Close();
}

void histogramRead() { 

   TFile * file = new TFile("histogram.root");

   TH1 * h1 = 0;
   file->GetObject("h1",h1);
   // you can also use nut you need to cast if you compile the code
   //TH1 * h1 = (TH1*) file->Get("h1");

   h1->Draw();
}

#include <vector>

class Particle {
public:
   Particle() {memset(fTags, 0, sizeof(fTags)); }
   double fPosX,fPosY,fPosZ; // particle position nearest to interaction point
   double fMomentum; // particle momentum
   double fMomentumPhi; // particle direction (phi)
   double fMomentumEta; // particle direction (eta)
   Long64_t fTags[128]; // particle tags
};

class EventData {
public:
   std::vector<Particle> fParticles; // particles of the event
   int fEventSize; // size (in bytes) of the event

   void SetSize() {
      fEventSize = sizeof(EventData) + fParticles.size() * sizeof(Particle);
   }
   void Clear() {
      fParticles.clear();
   }
   void AddParticle(const Particle& p) { fParticles.push_back(p); }

   //ClassDef(EventData,1); // Data for an event
};

#ifdef __MAKECINT__ 
#pragma link C++ class Particle+;
#pragma link C++ class EventData+;
#endif

// Prefer compiled:
#include "TTree.h"
#include "TFile.h"
#include "TRandom.h"
#include "TMath.h"
#include <vector>

#include "EventData.h"

void createTree(ULong64_t numEvents = 200) {
   TFile* f = new TFile("eventdata.root", "RECREATE");
   TTree* tree = new TTree("EventTree", "Tutorial tree");

   EventData* event = new EventData();
   tree->Branch("event", &event);

   Particle p;

   for (ULong64_t i = 0; i < numEvents; ++i) {       event->Clear();

      // generate event information
      int nParticles = 10 * gRandom->Exp(10.);

      // for every particle in the event generate the particle data 
      for (int ip = 0; ip < nParticles; ++ip) {          do {             p.fPosX = gRandom->Gaus(gRandom->PoissonD(0.1), 1.)
               + gRandom->BreitWigner(0.1, 0.1);;
         } while (fabs(p.fPosX) > 10.);
         p.fPosY = gRandom->Gaus(gRandom->PoissonD(0.01), .7);
         p.fPosZ = gRandom->Gaus(gRandom->PoissonD(10), 19.);

         p.fMomentum = gRandom->Exp(12);
         p.fMomentumPhi = gRandom->Uniform(2*TMath::Pi());
         do {
            p.fMomentumEta = gRandom->BreitWigner(0.01, 10.);
         } while (fabs(p.fMomentumEta) > 12.);

         
         event->AddParticle(p);
      }
      event->SetSize();

      tree->Fill();

      if (i % (numEvents/50) == 0) {
         printf("*");
         fflush(stdout);
      }
   }
   printf("\n");
   tree->Write();
   delete f;
}

• Open the file using its file name in TFile::Open() and get the Tree. Remember to check if the file pointer is not null. If it is null means the file is not existing.
• Get then a pointer to the tree.
• Connect a Tree Branch with the Data Member.We have to somehow connect the branch we want to read with the variables used to actually store the data by calling TTree::SetBranchAddress().
• Load the TTree data. For the analysis example we need to access the events' size, which is stored in the variable eventSize. But the TTree first needs to load the data for each event it contains. For that call TBranch::GetEntry(entry) in a loop, passing the TTree entry number from the loop index to GetEntry(). Again TBranch is the class name, but you obviously need to call it on an object. To know how many entries the tree contains, simply call TTree::GetEntries(). The branch is stored in eventSizeBranch.
• In the same loop, compute the total size of all events (simply add the current event size to the total size)
• Without the call to GetEntry(), the variables will not contain data. GetEntry() loads the data into the variables connected with the tree by the call to SetBranchAddress().
• Access the Analysis Result. At the end of the loop, print the sum of all event sizes. This sum shows you the real power of a TTree: even though you can analyze large amounts of data (our example tree with 22MB is tiny!) ROOT needs just a few MB of your RAM, no matter how many events you analyze. Imagine what it would be like if you had to load all data into memory, e.g. using a simple vector<EventData>

#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"

void AnalyzeTree()
{
   // Variables used to store the data
   Int_t     totalSize = 0;        // Sum of data size (in bytes) of all events
   Int_t     eventSize = 0;        // Size of the current event
   TBranch  *eventSizeBranch = 0;  // Pointer to the event.fEventsize branch


   // open the file
   TFile *f = TFile::Open("eventdata.root");
   if (f == 0) {
      // if we cannot open the file, print an error message and return immediatly
      printf("Error: cannot open eventdata.root!\n");
      return;
   }
   // get a pointer to the tree
   TTree *tree = (TTree *)f->Get("EventTree");



   // To use SetBranchAddress() with simple types (e.g. double, int)
   // instead of objects (e.g. std::vector<Particle>).
   tree->SetMakeClass(1);

   // Connect the branch "fEventSize" with the variable 
   // eventSize that we want to contain the data.
   // While we are at it, ask the tree to save the branch 
   // in eventSizeBranch
   tree->SetBranchAddress("fEventSize", &eventSize, &eventSizeBranch);

  

   // First, get the total number of entries
   Long64_t nentries = tree->GetEntries();
   // then loop over all of them
   for (Long64_t i=0;i<nentries;i++) {
      // Load the data for TTree entry number "i" from branch 
      // fEventSize into the connected variable (eventSize):
      eventSizeBranch->GetEntry(i);
      // compute the total size of all events
      totalSize += eventSize;
   }


   Int_t sizeInMB = totalSize/1024/1024;
   printf("Total size of all events: %d MB\n", sizeInMB);
}

You can also download the macro AnalyzeTree.C.

#include "TRandom.h"
#include "TFile.h"
#include "TTree.h"

void SimpleTree(const char * filename= "tree.root") {

   TTree data("tree","Example TTree");
   double x, y, z, t;
   data.Branch("x",&x,"x/D");
   data.Branch("y",&y,"y/D");
   data.Branch("z",&z,"z/D");
   data.Branch("t",&t,"t/D");

// fill it with random data                                                                                                   
   for (int i = 0; i<10000; ++i) {
      x = gRandom->Uniform(-10,10);
      y = gRandom->Gaus(0,5);
      z = gRandom->Exp(10);
      t = gRandom->Landau(0,2);

      data.Fill();
   }
// write in a file                                                                                                            
   TFile f(filename,"RECREATE");
   data.Write();
   f.Close();

}

This are the few lines to create the TChain, that you can run directly from the prompt. You can also use wildcard's to chain many files

TChain chain("tree");
chain.Add("tree*.root")
chain.Draw("t")

#include "TRandom.h"
#include "TFile.h"
#include "TTree.h"

void SimpleTree2() {

   TTree data("tree_2","Example TTree");
   double x, u;
   data.Branch("x",&x,"x/D");
   data.Branch("u",&u,"u/D");


// fill it with random data                                                                                                   
   for (int i = 0; i<10000; ++i) {
      x = gRandom->Exp(100);
      u = gRandom->Uniform(0,10);

      data.Fill();
   }
// write in a file                                                                                                            
   TFile f("tree_2.root","RECREATE");
   data.Write();
   f.Close();

}

Here are the lines of codes to Draw the x of the first tree versus the x of the second tree with a selection depending on u and t. You can run these lines from the ROOT prompt.

TFile f("tree.root");   // to get  the first tree
 tree->AddFriend("tree_2","tree_2.root");         
tree->Draw("x:tree_2.x","t<100 && tree_2.u<6","COLZ");

tree->MakeSelector("MySelector.C");

The file MySelector.h and MySelector.C will be created. Add in MySelector.h, inside the class MySelector, a new data member, the histogram you want to create,

TH1D * h_t;

Edit then the file MySelector.C and add in MySelector::SlaveBegin the booking of the histogram.

h_t = new TH1D("h_t","t",100,0,100);

In MySelector::Process the filling of the histogram after calling TSelector::GetEntry()

GetEntry(entry);
h_t->Fill(t);

In MySelector::Terminate the drawing of the histogram.

After having saved the file run the selection by doing (for example from the ROOT prompt):

TFile f("tree.root");
tree->Process("MySelector.C+");

To get access to the TFitResult object after fitting use option "S", as shown in slide 10 of the lecture.

Use gStyle->SetOptFit(1) to display the fit result in the statistics box.

#include "TF1.h"
#include "TH1.h"
#include "TFitResult.h"
#include "TMatrixDSym.h"
#include "TStyle.h"

void gausFit() { 

   TH1D * h1 = new TH1D("h1","h1",50,-5,5);

   h1->FillRandom("gaus",1000);

   TF1 * f1 = new TF1("f1","gaus"); 
   
   // add also option "Q" (quite) to avoid prrinting two time the result
   TFitResultPtr r = h1->Fit(f1,"S Q");
   
   // print the result
   r->Print(); 

   // get the correlation matrix and print it
   TMatrixDSym corrMatrix = r->GetCorrelationMatrix(); 
   corrMatrix.Print(); 

   gStyle->SetOptFit(1);

  // to get the sigma of the gaussian
  std::cout  << "Gaussian sigma = " << f1->GetParameter("sigma") << "  +/- " << f1->GetParError(f1->GetParNumber("sigma")) << std::endl;

}

TF1 * f1 = new TF1("f1","[0]*ROOT::Math::normal_pdf(x,[2],[1])");  
// set the parameters (needed if not using a pre-defined function as "gauss")
 f1->SetParameters(1,0,1); 

After the fit works, you can compute the number of peak events, by using TF1::Integral on the Gaussian only function. For the error you can use TF1::IntegralError, but you need to extract the correlation matrix from the fit and use the sub-matrix referring to the gaussian part.

#include "TF1.h"
#include "TH1.h"
#include "TFitResult.h"
#include "TMatrixDSym.h"
#include "TRandom.h"
#include "TStyle.h"

void PeakFit() { 

   TH1D * h1 = new TH1D("h1","h1",50,0,10);

   for (int i = 0; i < 100; ++i) {        h1->Fill(gRandom->Gaus(2,0.3)); 
   }
   for (int i = 0; i < 1000; ++i) {        h1->Fill(gRandom->Exp(2)); 
   }

   

   // fit first a single gaussian in range [0,5]
   TFitResultPtr r1 = h1->Fit("gaus","S","",1,3);  // first fit of gaussian
   TFitResultPtr r2 = h1->Fit("expo","S");  // first exponential 

   TF1 * f1 = new TF1("fitFunc","gaus(0)+expo(3)");
   // get parameters and set in global TF1
  
   // parameters of first gaussian
   f1->SetParameter(0,r1->Parameter(0));
   f1->SetParameter(1,r1->Parameter(1));
   f1->SetParameter(2,r1->Parameter(2));
   // parameters of second gaussian
   f1->SetParameter(3,r2->Parameter(0));
   f1->SetParameter(4,r2->Parameter(1));



   TFitResultPtr res = h1->Fit(f1,"SL"); 

   TMatrixDSym cov = res->GetCovarianceMatrix(); 

   // compute number of signal events
   // fitted gaussian function
   TF1 * peakFunc = new TF1("peakFunc","gaus");
   peakFunc->SetParameter(0, f1->GetParameter(0));
   peakFunc->SetParameter(1, f1->GetParameter(1));
   peakFunc->SetParameter(2, f1->GetParameter(2));

   TMatrixDSym covPeak = cov.GetSub(0,2,0,2);

   double nsignal = peakFunc->Integral(0,10) / h1->GetBinWidth(1);
   double err  = peakFunc->IntegralError(0,10,peakFunc->GetParameters(), covPeak.GetMatrixArray())/ h1->GetBinWidth(2); 

   std::cout << "number of signal events = " << nsignal << " +/- " << err << std::endl;

      


   


   

}

• Select the fit panel in the Canvas/Tools menu or by right-clicking on the histogram
• Make the fit function (gaussian plus exponential) using the Operation/Add button
• Set the initial parameters by playing with the sliders
• Press the Fit button

• Use the syntax of the RooWorkspace factory to create first the variables (observables and parameters) of the Gaussian probability density function (p.d.f.), as shown in slide 15.

• You need to define the [value, min, max] of a variable only the first time you create in the factory. Afterwards you can reference the variable by its name.

• After you have created the variable and p.d.f in the workspace, you can access their pointers, by using RooWorkspace::var for variables or =RooWorkspace::pdf for p.d.f.
• You can then do as in slide 12 to generate the data set and plot it.
• To fit the data set with the pdf, you need to call the RooAbsPdf::fitTo. See slide 14 on how to fit and on how to plot the function after the fit.
• To save the model call, RooWorkspace::write(file_name). See slide 19.

// make a simple Gaussian model

#include "RooWorkspace.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooDataSet.h"
#include "RooPlot.h"

#include "TCanvas.h"

using namespace RooFit; 

void GaussianModel(int n = 1000) { 


   RooWorkspace w("w");
   // define a Gaussian pdf
   w.factory("Gaussian:pdf(x[-10,10],mu[1,-1000,1000],s[2,0,1000])");   

   RooAbsPdf * pdf = w.pdf("pdf");   // access object from workspace
   RooRealVar * x = w.var("x");   // access object from workspace
   

   // generate n gaussian measurement for a Gaussian(x, 1, 1);
   RooDataSet * data = pdf->generate( *x, n); 

   data->SetName("data");  

   // RooFit plotting capabilities 
   RooPlot * pl = x->frame(); 
   data->plotOn(pl); 
   pl->Draw(); 


   // remove this line if you want to fit the data
   return; 

   // now fit the data set 
   pdf->fitTo(*data); 
 

   // plot the pdf on the same RooPlot object we have plotted the data 
   pdf->plotOn(pl);

   pl->Draw();

   // import data in workspace (IMPORTANT for saving it )
   w.import(*data); 

   w.Print();

   // write workspace in the file (recreate file if already existing)
   w.writeToFile("GaussianModel.root", true);

   cout << "model written to file " << endl;

}
   

• Open the TFile object and use TFile::Get to get a pointer to the workspace using its name
• Once you have the workspace in memory, retrieve from it the p.d.f. and the data set with their names using RooWorkspace::pdf and RooWorkspace::data.
• Re-issue again the call to RooAbsPdf::fitTo. You can set the fit range using the RooFit::Range(xmim,xmas) command arg option in fitTo. See the reference documentation for all the possible options that you can pass (some are shown in the solution code).

#include "RooWorkspace.h"
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooPlot.h"
#include "RooDataSet.h"
#include "RooFitResult.h"

#include "TFile.h"


// roofit tutorial showing how to fit whatever model we get from a file 

// we assume the name of the workspace is w
// the name of the pdf is pdf 
// the name of the data is data
const char * workspaceName = "w"; 
const char * pdfName = "pdf"; 
const char * dataName = "data"; 


using namespace RooFit; 




void fitModel(const char * filename = "GaussianModel.root" ) { 

   // read file: 
   // following lines are for reading workspace
   // and to check that is fine 

   // Check if example input file exists
   TFile *file = TFile::Open(filename);
   
   // if input file was specified but not found, quit
   if(!file ){
      cout <<"file " << filename << " not found" << endl;
      return;
   } 

   // get the workspace out of the file
   RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
   if(!w){
      cout <<"workspace with name " << workspaceName << " not found" << endl;
      return;
   }

   // fit a pdf from workspace with name pdfName

   RooAbsPdf * pdf = w->pdf(pdfName);
   if (!pdf) { 
      w->Print();
      cout << "pdf with name " << pdfName << " does not exist in workspace " << endl;
      return;
   }

   // get the data out of the file
   RooAbsData* data = w->data(dataName);

   if(!data ){
      w->Print();
      cout << "data " << dataName << " was not found" <<endl;
      return;
   }

   //----------------------------------------------------------------------
   //// real code starts here 

   // get variable x (is the first of the data)
   RooRealVar * x = w->var("x"); 
   RooPlot * plot = x->frame();

   data->plotOn(plot);
   plot->Draw();

   // fit pdf - (example using option: save result and using diffferent minimizer

   // global fit  

   pdf->fitTo( *data );

   // for doing a reduce fit in a Range (plus other options)
   
   RooFitResult * r = pdf->fitTo( *data, Save(true), Minimizer("Minuit2","Migrad"),Range(0.,10.) );

   pdf->plotOn(plot);
   pdf->paramOn(plot);

   plot->Draw();

   // if we have a result we can do 

   r->Print();

   r->correlationMatrix().Print();


}

• Use the syntax of the RooWorkspace factory to create first the Gaussian p.d.f. function with the variables (observables and parameters) and the Exponential
• Create then a RooAddPdf using the Gaussian and the Exponential and the relative number of events. Note that we could instead of the number of events a relative fraction. In this last case we would fit only for the shape and not the normalisation.
• The RooAddPdf is created using the SUM operator of the factory syntax.

// create model
RooWorkspace w("w"); 
 w.factory("Exponential:bkg_pdf(x[0,10], a[-0.5,-2,-0.2])");
 w.factory("Gaussian:sig_pdf(x, mass[2], sigma[0.3])");

// create RooAddPdf
 w.factory("SUM:model(nsig[100,0,10000]*sig_pdf, nbkg[1000,0,10000]*bkg_pdf)");  // for extended model

• After you have created the variable and p.d.f in the workspace, you can access their pointers, by using RooWorkspace::var for variables or RooWorkspace::pdf for p.d.f. Example:

// retrieving model components
RooAbsPdf * pdf = w.pdf("pdf");   // access object from workspace
RooRealVar * x = w.var("x");   // access object from workspace

• We generate the dataset using the generate function passing the observables we want to generate and the number of events. In case of an extended pdf (as the RooAddPdf we have created), if we don't pass the number of events, the number of events will be generated according to a Poisson distribution given the expected events (nsig+nbkg).
• Note that you can generate also binned data sets. In that case you have two options:
  • • use a RooDataHist and generateBinned: RooDataHist * binData = pdf->generateBinned(*x, 1000);.
    • use the option AllBinned() which will generate a weighted data set: RooDataSet * data = pdf->generate(*x, NumEvents(1000), AllBinned());.

// generate an unbinned dataset of 1000 events
RooDataSet * data = pdf->generate( *x, 1000); 

// create a RooPlot object and plot the data
RooPlot * plot = x->frame(); 
data->plotOn(plot); 

// fit the data and save its result. Use the optional =Minuit2= minimiser
RooFitResult * res = pdf->fitTo(*data, Minimizer("Minuit2","Migrad"), Save(true) );

• Plot the resulting fit function in the same plot.
• Plot also the signal and background fit components with different colour and style

pdf->plotOn(plot); 
//draw the two separate pdf's
pdf->plotOn(plot, RooFit::Components("bkg_pdf"), RooFit::LineStyle(kDashed) );
pdf->plotOn(plot, RooFit::Components("sig_pdf"), RooFit::LineColor(kRed), RooFit::LineStyle(kDashed) );
plot->Draw();   // to show the RooPlot in the current ROOT Canvas

w.writeToFile("GausExpModel.root",true);  // true means the file is re-created 

#include "RooWorkspace.h"
#include "RooAbsPdf.h"
#include "RooDataSet.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "RooRealVar.h"
#include "RooRandom.h"

#include "RooStats/ModelConfig.h"

using namespace RooFit; 



void GausExpModel(int nsig = 100,    // number of signal events 
                  int nbkg = 1000 )  // number of background events
{ 

   RooWorkspace w("w"); 
   w.factory("Exponential:bkg_pdf(x[0,10], a[-0.5,-2,-0.2])");
   w.factory("Gaussian:sig_pdf(x, mass[2], sigma[0.3])");

   w.factory("SUM:model(nsig[0,10000]*sig_pdf, nbkg[0,10000]*bkg_pdf)");  // for extended model

   RooAbsPdf * pdf = w.pdf("model");
   RooRealVar * x = w.var("x");  // the observable

   // set the desired value of signal and background events
   w.var("nsig")->setVal(nsig);
   w.var("nbkg")->setVal(nbkg);

   // generate the data

   // use fixed random numbers for reproducibility (use 0 for changing every time)
   RooRandom::randomGenerator()->SetSeed(111);

   // fix number of bins to 50 to plot or to generate data (default is 100 bins) 
   x->setBins(50);

   RooDataSet * data = pdf->generate( *x);  // will generate accordint to total S+B events
   //RooDataSet * data = pdf->generate( *x, AllBinned());  // will generate accordint to total S+B events
   data->SetName("data");
   w.import(*data);

   data->Print(); 


   RooPlot * plot = x->frame(Title("Gaussian Signal over Exponential Background"));
   data->plotOn(plot);
   plot->Draw();

   RooFitResult * r = pdf->fitTo(*data, RooFit::Save(true), RooFit::Minimizer("Minuit2","Migrad"));
   r->Print();

   pdf->plotOn(plot);
   //draw the two separate pdf's
   pdf->plotOn(plot, RooFit::Components("bkg_pdf"), RooFit::LineStyle(kDashed) );
   pdf->plotOn(plot, RooFit::Components("sig_pdf"), RooFit::LineColor(kRed), RooFit::LineStyle(kDashed) );

   pdf->paramOn(plot,Layout(0.5,0.9,0.85));

   plot->Draw();


   RooStats::ModelConfig mc("ModelConfig",&w);
   mc.SetPdf(*pdf);
   mc.SetParametersOfInterest(*w.var("nsig"));
   mc.SetObservables(*w.var("x"));
   // define set of nuisance parameters
   w.defineSet("nuisParams","a,nbkg");

   mc.SetNuisanceParameters(*w.set("nuisParams"));

   // import model in the workspace 
   w.import(mc);

   // write the workspace in the file
   TString fileName = "GausExpModel.root";
   w.writeToFile(fileName,true);
   cout << "model written to file " << fileName << endl;
}

• Gaussian plus Exponential fit result:

Error: (3) can't find GausExpModelFit.png in Main

AsymptoticCalculator   ac(*data, *sbModel, *bModel);
ac.SetOneSidedDiscovery(true);  // for one-side discovery test

// this run the hypothesis test and print the result
HypoTestResult * result = ac.GetHypoTest();
result->Print();

using namespace RooStats;
using namespace RooFit;

void SimpleHypoTest( const char* infile =  "GausExpModel.root", 
                     const char* workspaceName = "w",
                     const char* modelConfigName = "ModelConfig",
                     const char* dataName = "data" )
{

  /////////////////////////////////////////////////////////////
  // First part is just to access the workspace file 
  ////////////////////////////////////////////////////////////

  // open input file 
  TFile *file = TFile::Open(infile);
  if (!file) return;

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);


  // get the data  out of the file
  RooAbsData* data = w->data(dataName);

  // get the modelConfig (S+B) out of the file
  // and create the B model from the S+B model
  ModelConfig*  sbModel = (RooStats::ModelConfig*) w->obj(modelConfigName);
  sbModel->SetName("S+B Model");      
  RooRealVar* poi = (RooRealVar*) sbModel->GetParametersOfInterest()->first();
  poi->setVal(50);  // set POI snapshot in S+B model for expected significance
  sbModel->SetSnapshot(*poi);
  ModelConfig * bModel = (ModelConfig*) sbModel->Clone();
  bModel->SetName("B Model");      
  poi->setVal(0);
  bModel->SetSnapshot( *poi  );

  // create the AsymptoticCalculator from data,alt model, null model
  AsymptoticCalculator  ac(*data, *sbModel, *bModel);
  ac.SetOneSidedDiscovery(true);  // for one-side discovery test
  //ac.SetPrintLevel(-1);  // to suppress print level 

  // run the calculator
  HypoTestResult * asResult = ac.GetHypoTest();
  asResult->Print();


  return;  // comment this line if you want to run the FrequentistCalculator
  std::cout << "\n\nRun now FrequentistCalculator.....\n";   FrequentistCalculator   fc(*data, *sbModel, *bModel);   fc.SetToys(2000,500);    // 2000 for null (B) and 500 for alt (S+B)    // create the test statistics   ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
  // use one-sided profile likelihood
  profll.SetOneSidedDiscovery(true);

  // configure  ToyMCSampler and set the test statistics
  ToyMCSampler *toymcs = (ToyMCSampler*)fc.GetTestStatSampler();
  toymcs->SetTestStatistic(&profll);
  
  if (!sbModel->GetPdf()->canBeExtended())
     toymcs->SetNEventsPerToy(1);
 
  // run the test
  HypoTestResult * fqResult = fc.GetHypoTest();
  fqResult->Print();

  // plot test statistic distributions
  HypoTestPlot * plot = new HypoTestPlot(*fqResult);
  plot->SetLogYaxis(true);
  plot->Draw();

}

• Result of the significance test on the Gaussian Signal plus Background model (nsig=50, nbckg=1000):

Results HypoTestAsymptotic_result: 
 - Null p-value = 0.00118838
 - Significance = 3.0386
 - CL_b: 0.00118838
 - CL_s+b: 0.502432
 - CL_s: 422.787

• set the test statistic we want to use
• set the number of toys for the null and alternate model

FrequentistCalculator   fc(*data, *sbModel, *bModel);
fc.SetToys(2000,500);    // 2000 for null (B) and 500 for alt (S+B) 

// create the test statistics
ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
// use one-sided profile likelihood
profll.SetOneSidedDiscovery(true);

// configure  ToyMCSampler and set the test statistics
ToyMCSampler *toymcs = (ToyMCSampler*)fc.GetTestStatSampler();
toymcs->SetTestStatistic(&profll);
  
if (!sbModel->GetPdf()->canBeExtended())
   toymcs->SetNEventsPerToy(1);
 

HypoTestResult * fqResult = fc.GetHypoTest();
fqResult->Print();

// plot test statistic distributions
HypoTestPlot * plot = new HypoTestPlot(*fqResult);
plot->SetLogYaxis();
plot->Draw();

• Result of the significance test on the Gaussian Signal plus Background model (nsig=50, nbckg=1000) using the Frequentist calculator:

Results HypoTestCalculator_result: 
 - Null p-value = 0.001 +/- 0.000706753
 - Significance = 3.09023 +/- 0.2099 sigma
 - Number of Alt toys: 500
 - Number of Null toys: 2000
 - Test statistic evaluated on data: 4.61654
 - CL_b: 0.001 +/- 0.000706753
 - CL_s+b: 0.496 +/- 0.02236
 - CL_s: 496 +/- 351.262

Test statistics distribution for the B model (red) and S+B model obtained with the Frequentist calculator

using namespace RooStats;
using namespace RooFit;

void p0Plot() { 

  const char* infile =  "GausExpModel.root"; 
  const char* workspaceName = "w";
  const char* modelConfigName = "ModelConfig";
  const char* dataName = "data";

  /////////////////////////////////////////////////////////////
  // First part is just to access the workspace file 
  ////////////////////////////////////////////////////////////

  // Check if example input file exists
  TFile *file = TFile::Open(infile);

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);


  // get the modelConfig out of the file
  RooStats::ModelConfig* mc = (RooStats::ModelConfig*) w->obj(modelConfigName);

  // get the modelConfig out of the file
  RooAbsData* data = w->data(dataName);

  // get the modelConfig (S+B) out of the file
  // and create the B model from the S+B model
  ModelConfig*  sbModel = (RooStats::ModelConfig*) w->obj(modelConfigName);
  sbModel->SetName("S+B Model");      
  RooRealVar* poi = (RooRealVar*) sbModel->GetParametersOfInterest()->first();
  poi->setVal(50);  // set POI snapshot in S+B model for expected significance
  sbModel->SetSnapshot(*poi);
  ModelConfig * bModel = (ModelConfig*) sbModel->Clone();
  bModel->SetName("B Model");      
  poi->setVal(0);
  bModel->SetSnapshot( *poi  );

  vector<double> masses;
  vector<double> p0values;
  vector<double> p0valuesExpected;

  double massMin = 0.5;
  double massMax = 8.5;

  // loop on the mass values 
 
  for( double mass=massMin; mass<=massMax; mass += (massMax-massMin)/40.0 )   {                    cout << endl << endl << "Running for mass: " << mass << endl << endl;      w->var("mass")->setVal( mass );

     AsymptoticCalculator *  ac  = new AsymptoticCalculator(*data, *sbModel, *bModel);
     ac->SetOneSidedDiscovery(true);  // for one-side discovery test                                      
     AsymptoticCalculator::SetPrintLevel(-1);


     HypoTestResult* asymCalcResult = ac->GetHypoTest();
 
     asymCalcResult->Print();
     
     masses.push_back( mass );
     p0values.push_back( asymCalcResult->NullPValue() );
  
     double expectedP0 = AsymptoticCalculator::GetExpectedPValues(  asymCalcResult->NullPValue(),  asymCalcResult->AlternatePValue(), 0, false);
     p0valuesExpected.push_back( expectedP0 );
     std::cout << "expected p0 = " << expectedP0 << std::endl;   }   TGraph * graph1  = new TGraph(masses.size(),&masses[0],&p0values[0]);   TGraph * graph2  = new TGraph(masses.size(),&masses[0],&p0valuesExpected[0]);   graph1->SetMarkerStyle(20);
  graph1->Draw("APC");
  graph2->SetLineStyle(2);
  graph2->Draw("C");
  graph1->GetXaxis()->SetTitle("mass");
  graph1->GetYaxis()->SetTitle("p0 value");
  graph1->SetTitle("Significance vs Mass");
  graph1->SetMinimum(graph2->GetMinimum());
  graph1->SetLineColor(kBlue);
  graph2->SetLineColor(kRed);
  gPad->SetLogy(true);
  
     
}

p value obtained as function of the signal mass in the Gaussian plus exponential background model