ROOT Tutorials Exercises for Desy C++ School

Topics

This tutorial aims to provide a complement to the lecture shown at the Desy C++ School, November 2013. This tutorial is focus on data analysis in ROOT using Tree's and Proof. It does not cover other features of ROOT such as histogramming or fitting techniques. The main features of ROOT are presented: histogramming, data analysis using trees and advanced fitting techniques.

Material for the course

The slides of the lectures are available in electronic form from the school Agenda page. As complementary material, one can look at
  • ROOT Primer Guide, available in pdf, html or epub format. This introductory guide illustrates the main features of ROOT, relevant for the typical problems of data analysis: input and plotting of data from measurements and fitting of analytical functions.
  • ROOT user guide. It can be downloaded in various format (or only individual chatters) from here.
  • RooFit User Guide, available in pdf format. A coincide RooFit quick start guide is also available here.
  • A tar file with all the exercise solutions (all the ROOT macro required) is available here

Introduction

We will focus first on simpler exercises. The solution of the exercise is often to write a running ROOT macro. Two levels of help will be provided: a hint and a full solution. Try not to jump straight to the solution even if you experience some frustration. The help is organised as follow:

Here the hint is shown.

Here the solution is shown.

Some points linked to the exercises are optional and marked with a Pointing hand icon: they are useful to scrutinise in more detail some aspects. Try to solve them if you have the time.

Setting up ROOT

ROOT is installed in the virtual machines prepared for the course. It should be automatically available from the Terminal. If not, one should run the fallowing lines (or put in the default login shell script):

export ROOTSYS=/usr/local/ROOT
export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
. $ROOTSYS/bin/thisroot.sh

The version of ROOT available in the school virtual machines which should be 5.34.11.

Setting up ROOT pre-release version 6

Some exercises require some new functionality of ROOT, which will be available only in the new version 6, expected for the end of the year. In order to solve for this exercise and try the new version, you need to set-up a different environment variable for $ROOTSYS.

ROOT I/O and Trees

Set of Exercises working with the Trees in ROOT. First will start with an exercise on the I/O of ROOT by storing and reading an histogram from a file

Exercise 1: Writing and Reading Histogram from a file

Open a file then create a simple histogram, for example an histogram generated with exponential distribution. Fit it and write it in a file. Why the ROOT Canvas does not show the histogram ?

Use TFile::Open to open the file or just create a TFile object. Call TH1::Write to write the histogram in the file after having filled it.

#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();
}
The histogram is not shown, becausem when the file is close, it is automatically deleted.

Now read the histogram from the file and plot it.

Create a file object (or call TFile::Open) and then TFile::Get

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

Pointing hand You can also use the TBrowser to open the file and display the histogram.

Pointing hand What is going to happen if you delete the file after having retrieved the histogram from the file ?

Exercise 2: Creating a ROOT Tree

Create a simple ROOT tree containing 4 variables (for example x,y,z,t). Fill the tree with data (for example 10000 events) where x is generated according to a uniform distribution, y a gaussian and z an exponential and t a Landau distribution. Write also the tree in the file.

Create the Tree using simple variables as described in the lecture slides. Alternatively you can use also the TNtuple class.

Afterwards having saved the file, re-open the file and get the tree. Plot each single variable and also one variable versus another one (for example x versus y) using TTree::Draw. You can also use the TBrowser

--+++ Exercise 3: Creating a ROOT Tree containing a collection of LorentzVector's

Create a TTree containing a collection of 4D LorentzVectors. For example one could generated a list of pions (let's suppose 20/event in average) with an exponential distribution in pt and a uniform distribution in phi and Eta.

Measure the time to write a TTree with 100000 events.

Try to generate the tree in split mode (default) and no-split mode. Did you see a difference in performances in writing ?

<span class='twikiAlert'>
      Failed to include URL https://twiki.cern.ch/twiki/pub/Main/ROOTDesyTutorial2013vectorCollection.C Not Found
</span>
</div></div> <!--/twistyPlugin-->

Try now to use a TClonesArray and a TObjArray and measure the time to create the tree. Did you see an increase/decrease in performances ? 

--+++ Exercise 4: Read a !ROOT Tree containing a collection of LorentzVector's

Using TTree::Draw plot: 

 - the number of tracks per event
 - a profile plot showing the pt of the number of tracks vs the event number
 - a candle plot with the pt of the first 5 tracks when the number of tracks is larger or equal 5. 

Using some C/C++ code plot the Px  distribution of the tracks. Measure the time to read the tree and compare in case the tree was generated with splitting or not splitting.

Plot also the invariant mass of all tracks combinations. This you cannot do with TTree;:Draw. 
 

- try to do  this also using a Proxy function (see  TTree::MakeProxy) and/or  or TTree::MakeClass or the Selector. 

---+++ Exercise 2: Creating a !ROOT Tree containing an EventData object

Create a !ROOT tree which contains an !EventData object. The Event data object is defined in the file below
<div class="twistyPlugin twikiMakeVisibleInline">  <span id="twistyIdMainROOTDesyTutorial20139show" class="twistyTrigger twikiUnvisited twistyHidden twistyInited"><a href="#"><span class="twikiLinkLabel twikiUnvisited">Show&nbsp;</span></a> </span> <span id="twistyIdMainROOTDesyTutorial20139hide" class="twistyTrigger twikiUnvisited twistyHidden twistyInited"><a href="#"><span class="twikiLinkLabel twikiUnvisited">Hide&nbsp;</span></a> </span>  </div><!--/twistyPlugin twikiMakeVisibleInline--> <div class="twistyPlugin"><div id="twistyIdMainROOTDesyTutorial20139toggle" class="twistyContent twikiMakeHidden twistyInited">
%CODE{"cpp" style="background: yellow;"  }%
#ifndef EventData_h
#define EventData_h


#include <vector>
#include "Math/Vector4D.h"
#include "Math/Vector3D.h"


class Particle {
public:
   Particle() { //memset(fTags, 0, sizeof(fTags));
 }
  
   
   
   ROOT::Math::XYZVector  fPosition;   // vertex position
   ROOT::Math::PtEtaPhiMVector fVector;     // particle vector 
   int  fCharge;        // particle charge
   int fType;           // particle type (0=pho, 1 ele, 2...)
};

class TF1; 

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

   void Generate();

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


#ifdef __MAKECINT__ 
#pragma link C++ class Particle+;
#pragma link C++ class EventData+;
#pragma link C++ class std::vector<Particle>+;
#endif

#endif
   

The EventData object provides a method (Generate() ) to generate one event and it is implemented in the file below

#include "EventData.h"
#include "TRandom.h"
#include "TGenPhaseSpace.h"
#include "Math/Vector4D.h"
#include "TLorentzVector.h"

// generate the events

// parameters

int NTRACKS = 40; 
// parameter for pt distribution
double PtAvg = 20;   
int NTYPES = 8;
double masses[] = { 0, 0.0005, 0.105, 0.135, 0.139, 0.4937, 0.4976, 3.096 }; 
int charge [] = { 0, 1, 1, 0, 1, 1, 0, 0 }; 
double fractions[] = { 0.1, 0.12, 0.08, 0.05, 0.5, 0.1, 0.049 , 0.01 };  
double sigmax = 10.E-6; 
double sigmay = 10.E-6; 
double sigmaz = 5.; 

using namespace ROOT::Math;

ROOT::Math::PtEtaPhiMVector SmearVector(const ROOT::Math::XYZTVector & v) { 
   double x = v.X()*(1. + gRandom->Gaus(0, 0.1) );
   double y = v.Y()*(1. + gRandom->Gaus(0, 0.1) );
   double z = v.Z()*(1. + gRandom->Gaus(0, 0.1) );
   ROOT::Math::PxPyPzMVector tmp(x,y,z,v.M() );
   return PtEtaPhiMVector(tmp);
}


void EventData::Generate()  { 


   // get expected value for each type
   for (int i = 0; i < NTYPES; ++i) {        double nexp = fractions[i] * NTRACKS;        int np = gRandom->Poisson(nexp); 
      for (int j = 0; j < np; ++j) {           Particle p;           p.fPosition = XYZVector( gRandom->Gaus(0,sigmax), gRandom->Gaus(0, sigmay), gRandom->Gaus(0, sigmaz) ); 
         double pt = gRandom->Exp(PtAvg); 
         double eta = gRandom->Uniform(-3,3);
         double phi = gRandom->Uniform(-TMath::Pi(), TMath::Pi() ); 
         double mass = masses[i];
         p.fVector = PtEtaPhiMVector(pt, eta, phi, mass); 
         p.fType = i; 
         p.fCharge = charge[i];
         if (p.fCharge) { 
            int tmp = gRandom->Integer(2); 
            if (tmp == 0) p.fCharge = -1; 
         }
         // special case for decays
         if (i == 3 ) { 
            // pi0
            TGenPhaseSpace evt; 
            double m[2] = {0,0};
            TLorentzVector W( p.fVector.X(), p.fVector.Y(), p.fVector.Z(), p.fVector.E() );
            evt.SetDecay(W, 2, m);
            evt.Generate();
            TLorentzVector * v1 = evt.GetDecay(0); 
            TLorentzVector * v2 = evt.GetDecay(1); 
            Particle p1; 
            Particle p2; 
            p1.fPosition = p.fPosition; 
            p2.fPosition = p.fPosition; 
            p1.fCharge = 0;
            p2.fCharge = 0;
            p1.fType = 0; 
            p2.fType = 0;

            p1.fVector = SmearVector(XYZTVector(v1->X(), v1->Y(), v1->Z(), v1->E() )); 
            p2.fVector = SmearVector(XYZTVector(v2->X(), v2->Y(), v2->Z(), v2->E() )); 
            AddParticle(p1);
            AddParticle(p2);
         }
         if (i == 6) { 
            // Ks
            TGenPhaseSpace evt; 
            double m[2] = {0.139,0.139};
            TLorentzVector W( p.fVector.X(), p.fVector.Y(), p.fVector.Z(), p.fVector.E() );
            evt.SetDecay(W, 2, m);
            evt.Generate();
            TLorentzVector * v1 = evt.GetDecay(0); 
            TLorentzVector * v2 = evt.GetDecay(1); 
            double ctau = 2.6844 * p.fVector.Gamma(); 
            double disp = gRandom->Exp(ctau); 
            double dispX = disp * p.fVector.X()/p.fVector.P(); 
            double dispY = disp * p.fVector.Y()/p.fVector.P(); 
            double dispZ = disp * p.fVector.Z()/p.fVector.P(); 
            Particle p1; 
            Particle p2; 
            p1.fPosition = XYZVector( p.fPosition.X() + dispX,  p.fPosition.Y() + dispY,  p.fPosition.Z() + dispZ);
            p2.fPosition = p1.fPosition; 
            p1.fCharge = 1;
            p2.fCharge = -1;
            p1.fType = 4; 
            p2.fType = 4;

            p1.fVector = SmearVector(XYZTVector(v1->X(), v1->Y(), v1->Z(), v1->E() )); 
            p2.fVector = SmearVector(XYZTVector(v2->X(), v2->Y(), v2->Z(), v2->E() )); 
            AddParticle(p1);
            AddParticle(p2);
         }
         if (i == 7 ) { 
            // J/psi
            TGenPhaseSpace evt; 
            double m1[2] = {0.0005,0.0005};
            double m2[2] = {0.105,0.105};
            int tmp = gRandom->Integer(2);
            TLorentzVector W( p.fVector.X(), p.fVector.Y(), p.fVector.Z(), p.fVector.E() );
            if (tmp == 0)
               evt.SetDecay(W, 2, m1);
            else 
               evt.SetDecay(W, 2, m2);

            evt.Generate();

            TLorentzVector * v1 = evt.GetDecay(0); 
            TLorentzVector * v2 = evt.GetDecay(1); 
            Particle p1; 
            Particle p2; 
            p1.fPosition = p.fPosition; 
            p2.fPosition = p.fPosition; 
            p1.fCharge = 1;
            p2.fCharge = -1;
            p1.fType = 1; 
            p2.fType = 1;
            if (tmp == 1) {
               p1.fType = 2; 
               p2.fType = 2;
            }

            p1.fVector = SmearVector(XYZTVector(v1->X(), v1->Y(), v1->Z(), v1->E() )); 
            p2.fVector = SmearVector(XYZTVector(v2->X(), v2->Y(), v2->Z(), v2->E() )); 
            AddParticle(p1);
            AddParticle(p2);
         }
         else
            AddParticle(p);

      }
   }
}

Try to write macro creating the tree and filling with them with the events containing the EventData object.

<span class='twikiAlert'>
     Error: File attachment at https://twiki.cern.ch/twiki/pub/Main/ROOTDesyTutorial2013/CreateTree.C,%PARAM2% does not exist 
</span>

Look at the macro and try to understand. Run the macro to create and write the tree in the file.

Exercise 6: Reading a ROOT Tree

Read the tree from the file and by using the TBrowser or TTree::Draw plot for example the momentum of all the particle or the event size.

Now read the Tree with a macro and calculate the sum of all event sizes.

  • 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.

Exercise 7: Chaining ROOT Trees.

Create a TTree now containing for each event simple double variables (like a tuple), for example, x, y, z and t.  Run the macro few times, but changing always the name of the file where the tree is stored. Use then the =TChain class to merge the trees

You can look at the tutorial tree/tree1.C as example, in particular how TTree:Branch is used to define the tree branches containing the variables. Then see the example in the lecture slide on how to use TTree::Chain.

Here is for example the macro to create several trees
#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")

Exercise 8: Using Tree Friends

Suppose we have the Tree with the vector of LorentzVector. We want now to add a new branch to this tree containing an std::vector< std::vector< double> > which contains the mass invariant of each track with respect all the other tracks.

Make a new Tree. Read from the file the Tree used before and add as a friend to this tree. Plot the x variable of the first tree versus the x variable of the second one.

Here is the macro to create a second tree, containing x, u:
#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");

Exercise 9: Using the TSelector class for analysing a TTree

Create the Selector for the EventData TTree made before. Use TTree::MakeSelector to create your own Selector class. The aim is to plot the invariant masses for: - the photons, - the opposite charged particles - for the muon and electrons

Inside the code of your Selector do the following:

  • book the histograms in the initialisation routine
  • fill the histogram in the Process function
  • draw the histogram in the Terminate function

Here is what you need to do, after having opened the file with the tree
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+");

See the attached file MySelector.h and MySelector.C.


This topic: Main > WebPreferences > RootGridKaTutorial2013 > ROOTDesyTutorial2013
Topic revision: r2 - 2013-11-13 - LorenzoMoneta
 
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