ROOT tutorial II

In the previous ROOT tutorial, we learned how to use basic ROOT classes such as TF1, TGraph, TRandom3, TH1, TH2, and TFile. In this lecture, we move to one of the most important classes in ROOT:

TTree

A TTree is used to store many events of the same type.
Each entry of the tree usually corresponds to one event.
Each branch stores one variable, or one group of closely related variables.

This is the standard data structure used in ROOT-based event analysis.

---

1. A first look at TTree

In a simple detector event, we may want to store:

• event ID
• measured energy
• time
• hit positions

Some variables are scalars, such as eventID or energy.

Some variables may naturally appear as arrays, such as signals from several detector channels.

---

2. Create a simple tree

In the following example, each event contains:

• evt: event number
• energy: measured total energy
• time: measured time
• adc[4]: ADC values from 4 channels
• x[2]: two reconstructed x positions

%jsroot on
TCanvas *c1 = new TCanvas("c1","c1",800,600);

TFile *fout = new TFile("tree_basic.root","recreate");
TTree *tree = new TTree("tree","simple event tree");

Int_t evt;
Float_t energy, t;
Int_t adc[4];

Now create branches:

tree->Branch("evt",    &evt,    "evt/I");
tree->Branch("energy", &energy, "energy/F");
tree->Branch("t",      &t,       "t/F");
tree->Branch("adc",    adc,     "adc[4]/I");

Now fill the tree event by event.

TRandom3 r1(0);

for(evt=0; evt<5000; evt++) {

    if(r1.Uniform()<0.65) energy = r1.Gaus(950, 70);
    else                  energy = r1.Gaus(1350, 90);

    t = 20 + 0.02*energy + r1.Gaus(0,4);

    adc[0] = r1.Poisson(0.08*energy);
    adc[1] = r1.Poisson(0.09*energy);
    adc[2] = r1.Poisson(0.07*energy);
    adc[3] = r1.Poisson(0.06*energy);

    tree->Fill();
}

Finally write the tree into the file:

tree->Write();
fout->Close();

3. Reopen the file and inspect the tree

Now reopen the file:

TFile *fin = new TFile("tree_basic.root");
fin->ls(); // use .ls in ROOT terminal

TFile**		tree_basic.root	
 TFile*		tree_basic.root	
  KEY: TTree	tree;1	simple event tree

Get the tree:

TTree *tree = (TTree*)fin->Get("tree"); //can be skiped in ROOT terminal

Print the tree structure:

• This is the first thing you should do after opening a tree. It shows the branch names and their types.

tree->Print();

******************************************************************************
*Tree    :tree      : simple event tree                                      *
*Entries :     5000 : Total =          142724 bytes  File  Size =      69907 *
*        :          : Tree compression factor =   2.03                       *
******************************************************************************
*Br    0 :evt       : evt/I                                                  *
*Entries :     5000 : Total  Size=      20547 bytes  File Size  =       7080 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=   2.83     *
*............................................................................*
*Br    1 :energy    : energy/F                                               *
*Entries :     5000 : Total  Size=      20562 bytes  File Size  =      17141 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=   1.17     *
*............................................................................*
*Br    2 :t         : t/F                                                    *
*Entries :     5000 : Total  Size=      20537 bytes  File Size  =      17184 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=   1.17     *
*............................................................................*
*Br    3 :adc       : adc[4]/I                                               *
*Entries :     5000 : Total  Size=      80693 bytes  File Size  =      27809 *
*Baskets :        3 : Basket Size=      32000 bytes  Compression=   2.88     *
*............................................................................*

Check the number of entries:

cout << "Entries = " << tree->GetEntriesFast() << endl;

Entries = 5000

Show one event:

tree->Show(10);// 11th event

======> EVENT:10
 evt             = 10
 energy          = 859.16
 t               = 34.4011
 adc             = 65, 
                  84, 75, 43

Scan several events:

tree->Scan("evt:energy:t:adc[0]:adc[1]","","",10,1); //show from 2nd event to 11th events

************************************************************************
*    Row   *       evt *    energy *         t *    adc[0] *    adc[1] *
************************************************************************
*        1 *         1 * 1023.4143 * 43.048313 *        76 *        90 *
*        2 *         2 * 939.31671 * 40.579158 *        61 *        83 *
*        3 *         3 * 943.59063 * 29.387554 *        80 *       103 *
*        4 *         4 * 887.30120 * 33.671905 *        65 *        83 *
*        5 *         5 * 956.72790 * 41.347721 *        74 *        90 *
*        6 *         6 * 1270.0981 * 39.712917 *       106 *       136 *
*        7 *         7 * 771.20013 * 35.892658 *        67 *        76 *
*        8 *         8 * 1261.6889 * 41.559322 *        94 *       140 *
*        9 *         9 * 875.31134 * 29.064889 *        72 *        79 *
*       10 *        10 * 859.15979 * 34.401138 *        65 *        84 *
************************************************************************

4. `tree->Draw(...)

Before writing explicit event loops, ROOT gives a very convenient tool:

tree->Draw(...)

• give a branch name or an expression
• ROOT loops over the tree automatically
• a histogram or graph is produced

Draw the energy distribution:

tree->Draw("energy");
c1->Draw();

• ROOT automatically creates a temporary histogram

Draw into a named histogram

Draw the energy into a histogram named hE:

tree->Draw("energy>>hE(100,800,1600)");
c1->Draw();

Now hE exists and can be used later:

hE->Draw();
c1->Draw();

Use selection cuts in Draw()

Very often we only want to draw events satisfying a condition.

For example:

tree->Draw("energy","t>50");
c1->Draw();

Two-dimensional Draw()

For a 2D plot, use:

tree->Draw("y:x");

• The variable after the colon is the x-axis.
• The variable before the colon is the y-axis.

tree->Draw("t:energy");
c1->Draw();

Draw into a 2D histogram

Just like in 1D, we can create a named histogram directly:

tree->Draw("t:energy>>hET(100,800,1600,100,20,80)","","colz");
c1->Draw();

5. Variable-length arrays

So far we have used fixed arrays like adc[4].

In real event data, however, the number of objects may vary from event to event.

For example:

• the number of hits in a detector
• the number of reconstructed tracks
• the number of clusters
• the number of jets in collider events

ROOT supports variable-length arrays.

The usual pattern is:

• one integer branch stores the number of items
• another branch stores an array whose valid length depends on that integer

For example:

Int_t nhit;
Float_t xpos[10];
Float_t ypos[10];

tree->Branch("nhit", &nhit, "nhit/I");
tree->Branch("xpos", xpos,  "xpos[nhit]/F");
tree->Branch("ypos", ypos,  "ypos[nhit]/F");

This means:

• for each event, nhit tells ROOT how many elements are valid
• only the first nhit elements of xpos and ypos are used

Example: a tree with variable-length hit arrays

TFile *fhit = new TFile("tree_hits.root","recreate");
TTree *thit = new TTree("thit","event tree with variable-length hit arrays");

Int_t evt;
Int_t nhit;
Float_t xpos[10], ypos[10], e[10];

thit->Branch("evt",  &evt,  "evt/I");
thit->Branch("nhit", &nhit, "nhit/I");
thit->Branch("xpos", xpos,  "xpos[nhit]/F");
thit->Branch("ypos", ypos,  "ypos[nhit]/F");
thit->Branch("e", e,  "e[nhit]/F");

TRandom3 r3(0);

for(evt=0; evt<3000; evt++) {

    nhit = r3.Integer(5);   // 0,1,2,3,4 hits

    for(int i=0; i<nhit; i++) {
        xpos[i] = r3.Gaus(0,10);
        ypos[i] = r3.Gaus(0,8);
        e[i] = r3.Exp(2.0);
    }

    thit->Fill();
}

thit->Write();
fhit->Close();

Now inspect it:

TFile *fhitin = new TFile("tree_hits.root");
TTree *thit = (TTree*)fhitin->Get("thit");

thit->Print();

******************************************************************************
*Tree    :thit      : event tree with variable-length hit arrays             *
*Entries :     3000 : Total =          134898 bytes  File  Size =      89054 *
*        :          : Tree compression factor =   1.50                       *
******************************************************************************
*Br    0 :evt       : evt/I                                                  *
*Entries :     3000 : Total  Size=      12547 bytes  File Size  =       4296 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=   2.81     *
*............................................................................*
*Br    1 :nhit      : nhit/I                                                 *
*Entries :     3000 : Total  Size=      12552 bytes  File Size  =       2153 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=   5.61     *
*............................................................................*
*Br    2 :xpos      : xpos[nhit]/F                                           *
*Entries :     3000 : Total  Size=      36534 bytes  File Size  =      27421 *
*Baskets :        2 : Basket Size=      32000 bytes  Compression=   1.31     *
*............................................................................*
*Br    3 :ypos      : ypos[nhit]/F                                           *
*Entries :     3000 : Total  Size=      36534 bytes  File Size  =      27415 *
*Baskets :        2 : Basket Size=      32000 bytes  Compression=   1.31     *
*............................................................................*
*Br    4 :e         : e[nhit]/F                                              *
*Entries :     3000 : Total  Size=      36516 bytes  File Size  =      26940 *
*Baskets :        2 : Basket Size=      32000 bytes  Compression=   1.33     *
*............................................................................*

thit->Scan("evt:nhit:ypos","","",10,1);

***********************************************************
*    Row   * Instance *       evt *      nhit *      ypos *
***********************************************************
*        1 *        0 *         1 *         4 * 4.5903854 *
*        1 *        1 *         1 *         4 * 12.433490 *
*        1 *        2 *         1 *         4 * -4.218780 *
*        1 *        3 *         1 *         4 * 9.8480567 *
*        2 *        0 *         2 *         4 * 14.526621 *
*        2 *        1 *         2 *         4 * 0.5912082 *
*        2 *        2 *         2 *         4 * 0.2238916 *
*        2 *        3 *         2 *         4 * 4.3442978 *
*        3 *        0 *         3 *         4 * -11.62974 *
*        3 *        1 *         3 *         4 * 12.614666 *
*        3 *        2 *         3 *         4 * -4.749458 *
*        3 *        3 *         3 *         4 * 12.738591 *
*        4 *        0 *         4 *         2 * 1.3567740 *
*        4 *        1 *         4 *         2 * 14.820048 *
*        5 *        0 *         5 *         3 * 10.654076 *
*        5 *        1 *         5 *         3 * -3.364059 *
*        5 *        2 *         5 *         3 * -1.159192 *
*        6 *        0 *         6 *         0 *           *
*        7 *        0 *         7 *         2 * -9.866870 *
*        7 *        1 *         7 *         2 * 1.4436848 *
*        8 *        0 *         8 *         4 * 2.2805223 *
*        8 *        1 *         8 *         4 * -1.107305 *
*        8 *        2 *         8 *         4 * 5.5729947 *
*        8 *        3 *         8 *         4 * 0.3025611 *
*        9 *        0 *         9 *         2 * 3.9707925 *
*        9 *        1 *         9 *         2 * -11.77236 *
*       10 *        0 *        10 *         1 * -1.187129 *
***********************************************************

Type <CR> to continue or q to quit ==> 

If nhit=0, the array values are not meaningful for that event, so selection is often useful.

Draw the multiplicity:

thit->Draw("nhit");
c1->Draw();

6. Reading a tree explicitly: SetBranchAddress() and GetEntry()

Draw() is excellent for quick checks and simple plots. But for more flexible analysis, we need to read entries explicitly.

The standard pattern is:

1. declare C++ variables
2. connect them to branches with SetBranchAddress()
3. loop over entries
4. call GetEntry(i)
5. do calculations inside the loop

Example using tree:

TFile *fin = new TFile("tree_basic.root");

Int_t evt_r;
Float_t energy_r, t_r;
Int_t adc_r[4];

tree->SetBranchAddress("evt",    &evt_r);
tree->SetBranchAddress("energy", &energy_r);
tree->SetBranchAddress("t",      &t_r);
tree->SetBranchAddress("adc",    adc_r);

Create a histogram:

TH1F *hsel = new TH1F("hsel","Selected energy;Energy;Counts",100,700,1600);

Loop over entries:

Long64_t nentries = tree->GetEntriesFast();

for(Long64_t i=0; i<nentries; i++) {
    tree->GetEntry(i);

    if(t_r > 45 && adc_r[0] > 70) {
        hsel->Fill(energy_r);
    }
}

Draw the result:

hsel->Draw();
c1->Draw();

Now the difference between Draw() and an explicit loop becomes clear:

• Draw() is fast and convenient
• explicit reading is more flexible

7. A more realistic event example

Now we use ROOT’s physics classes to generate events that look more like real particle or nuclear reaction data.

In this example, we generate events for:

$$ pp \rightarrow pp\pi^+\pi^- $$ using TLorentzVector and TGenPhaseSpace.

We will store the final-state particles in a tree and later read them back.

Create a tree for generated final-state particles

Each event will contain 4 final-state particles:

• proton 1
• proton 2
• pion +
• pion -

For each particle, we store:

• px
• py
• pz
• E

We also store the PDG ID.

TFile *fgen = new TFile("pp2pppipi.root","recreate");
TTree *tgen = new TTree("tgen","generated pp -> pp pi+ pi- events");

Int_t evt;
Int_t npart;
Int_t pid[10];
Float_t px[10], py[10], pz[10], E[10];

tgen->Branch("evt",   &evt,   "evt/I");
tgen->Branch("npart", &npart, "npart/I");
tgen->Branch("pid",   pid,    "pid[npart]/I");
tgen->Branch("px",    px,     "px[npart]/F");
tgen->Branch("py",    py,     "py[npart]/F");
tgen->Branch("pz",    pz,     "pz[npart]/F");
tgen->Branch("E",     E,      "E[npart]/F");

TRandom3 rng4(0);

const double mp  = 0.93827;
const double mpi = 0.13957;
const double sqrts = 3.0;

TLorentzVector pinit(0,0,0,sqrts);
TGenPhaseSpace event;

double masses[4] = {mp, mp, mpi, mpi};

for(evt=0; evt<20000; evt++) {

    if(!event.SetDecay(pinit,4,masses)) continue;
    event.Generate();

    npart = 4;

    Int_t pidtmp[4] = {2212,2212,211,-211};

    for(int i=0; i<npart; i++) {
        TLorentzVector *p = event.GetDecay(i);
        pid[i] = pidtmp[i];
        px[i] = p->Px();
        py[i] = p->Py();
        pz[i] = p->Pz();
        E[i]  = p->E();
    }

    tgen->Fill();
}

tgen->Write();
fgen->Close();

Inspection of the generated event tree

Open it:

TFile *fgenin = new TFile("pp2pppipi.root");
TTree *tgen = (TTree*)fgenin->Get("tgen");

Quick checks with Draw()

Now use Draw() immediately:

Draw pion momentum components:

TCanvas *c2 = new TCanvas("c2","c2", 900,300);
c2->Clear();
c2->Divide(3,1);
c2->cd(1);
tgen->Draw("px[2]");
c2->cd(2);
tgen->Draw("py[2]");
c2->cd(3);
tgen->Draw("pz[2]");
c2->Draw();

Physics Analysis: Reconstructing the $\pi^+\pi^-$ Invariant Mass

In this final step, we treat the stored data as a real physics experiment. The goal is to calculate the invariant mass of the two pions ($M_{\pi^+\pi^-}$) to see the kinematic distribution of the pair.

How we build it:

1. Identify Particles: We loop through the pid array to find the indices of the positive pion (211) and the negative pion (-211).
2. Assign Four-Momenta: We use the stored px, py, pz, E to create two TLorentzVector objects, one for each pion.
3. Vector Addition: We sum the two vectors to get the total four-momentum of the system:
   $P_{total} = P_{\pi^+} + P_{\pi^-}$
4. Calculate Mass: We use the .M() method of the total vector, which calculates:
   $M = \sqrt{E_{total}^2 - |\vec{p}_{total}|^2}$

This value is filled into the histogram hm. In a real experiment, if the two pions came from the decay of a parent particle (like a $\rho^0$ or $K_S^0$), you would see a sharp peak in this mass distribution.

Int_t evt_g, npart_g;
Int_t pid_g[10];
Float_t px_g[10], py_g[10], pz_g[10], E_g[10];

tgen->SetBranchAddress("evt",   &evt_g);
tgen->SetBranchAddress("npart", &npart_g);
tgen->SetBranchAddress("pid",   pid_g);
tgen->SetBranchAddress("px",    px_g);
tgen->SetBranchAddress("py",    py_g);
tgen->SetBranchAddress("pz",    pz_g);
tgen->SetBranchAddress("E",     E_g);

TH1F *hm = new TH1F("hm","M(#pi^{+}#pi^{-});M (GeV/c^{2});Counts",100,0.2,1.5);

Long64_t nentries = tgen->GetEntriesFast();

for(Long64_t i=0; i<nentries; i++) {
    tgen->GetEntry(i);

    int ipip = -1, ipim = -1;

    for(int j=0; j<npart_g; j++) {
        if(pid_g[j] ==  211) ipip = j;
        if(pid_g[j] == -211) ipim = j;
    }

    if(ipip<0 || ipim<0) continue;

    TLorentzVector ppip(px_g[ipip], py_g[ipip], pz_g[ipip], E_g[ipip]);
    TLorentzVector ppim(px_g[ipim], py_g[ipim], pz_g[ipim], E_g[ipim]);

    hm->Fill((ppip + ppim).M());
}

TCanvas *c4 = new TCanvas("c4","c4", 800,600);
hm->Draw();
c4->Draw();

Reference Values for Verification:

• $K_S^0$ (K-Short Meson):
  • Mass: $\approx 0.4976$ GeV/$c^2$ ($497.6$ MeV/$c^2$)
  • Note: This would appear as a very sharp, narrow peak in your $M_{\pi^+\pi^-}$ histogram.
• $\rho^0$ (Rho Meson):
  • Mass: $\approx 0.775$ GeV/$c^2$ ($775.3$ MeV/$c^2$)
  • Note: Unlike the $K_S^0$, the $\rho^0$ is a "resonance" with a very short lifetime, so it would appear as a much broader peak (about $0.15$ GeV wide).

Note on the result:

The resulting histogram shows a smooth continuum. This is because TGenPhaseSpace simulates a “pure phase space” distribution where particles do not interact with each other.