Hello everyone.
I am trying to read data from a ROOT file in order to draw some histograms.
What I have is the following. Inside of the ROOT file there is a TTree called summary, inside the tree there are two branches I need to read, injEnergy and nucEnergy. I need to read all data from each branch into a histogram. I am not sure how to do so. Below there is my code.
#include "TFile.h"
#include "TTree.h"
void ratioplot() {
// open the file
TFile *f = TFile::Open("/home/oscar/Downloads/v2r4/SimProp_N200000.root");
if (f == 0)
{
// if we cannot open the file, print an error message and return immediatly
printf("Error: cannot open file!\n");
}
TTreeReader myReader("summary", f);
// Define two histograms. Note the X and Y title are defined
// at booking time using the convention "Hist_title ; X_title ; Y_title"
TH1F *h1 = new TH1F("h1", "Initial and final energy spectrum for the arriving particles;x title; h1 and h2 gaussian histograms", 100, 0, 5e+19);
TH1F *h2 = new TH1F("h2", "h2", 100, -5, 5);
//h1->FillRandom("gaus"); //Here I think I should read the data from the ROOT file into the histograms
//h2->FillRandom("gaus");
// Define the Canvas
TCanvas *c = new TCanvas("c", "canvas", 800, 800);
// Upper plot will be in pad1
TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
pad1->SetBottomMargin(0); // Upper and lower plot are joined
pad1->SetGridx(); // Vertical grid
pad1->Draw(); // Draw the upper pad: pad1
pad1->cd(); // pad1 becomes the current pad
h1->SetStats(0); // No statistics on upper plot
h1->Draw(); // Draw h1
h2->Draw("same"); // Draw h2 on top of h1
// Do not draw the Y axis label on the upper plot and redraw a small
// axis instead, in order to avoid the first label (0) to be clipped.
h1->GetYaxis()->SetLabelSize(0.);
TGaxis *axis = new TGaxis( -5, 20, -5, 220, 20,220,510,"");
axis->SetLabelFont(43); // Absolute font size in pixel (precision 3)
axis->SetLabelSize(15);
axis->Draw();
// lower plot will be in pad
c->cd(); // Go back to the main canvas before defining pad2
TPad *pad2 = new TPad("pad2", "pad2", 0, 0.05, 1, 0.3);
pad2->SetTopMargin(0);
pad2->SetBottomMargin(0.2);
pad2->SetGridx(); // vertical grid
pad2->Draw();
pad2->cd(); // pad2 becomes the current pad
// Define the ratio plot
TH1F *h3 = (TH1F*)h1->Clone("h3");
h3->SetLineColor(kBlack);
h3->SetMinimum(0.8); // Define Y ..
h3->SetMaximum(1.35); // .. range
h3->Sumw2();
h3->SetStats(0); // No statistics on lower plot
h3->Divide(h2);
h3->SetMarkerStyle(21);
h3->Draw("ep"); // Draw the ratio plot
// h1 settings
h1->SetLineColor(kBlue+1);
h1->SetLineWidth(2);
// Y axis h1 plot settings
h1->GetYaxis()->SetTitleSize(20);
h1->GetYaxis()->SetTitleFont(43);
h1->GetYaxis()->SetTitleOffset(1.55);
// h2 settings
h2->SetLineColor(kRed);
h2->SetLineWidth(2);
// Ratio plot (h3) settings
h3->SetTitle(""); // Remove the ratio title
// Y axis ratio plot settings
h3->GetYaxis()->SetTitle("Ratio of final and initial energy spectrum ");
h3->GetYaxis()->SetNdivisions(505);
h3->GetYaxis()->SetTitleSize(20);
h3->GetYaxis()->SetTitleFont(43);
h3->GetYaxis()->SetTitleOffset(1.55);
h3->GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
h3->GetYaxis()->SetLabelSize(15);
// X axis ratio plot settings
h3->GetXaxis()->SetTitleSize(20);
h3->GetXaxis()->SetTitleFont(43);
h3->GetXaxis()->SetTitleOffset(4.);
h3->GetXaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
h3->GetXaxis()->SetLabelSize(15);
}
I would use the Project Function instead of the TTreeReader:
TTree* tree = (TTree*) f->Get("summary");
// Histogram declarations
tree->Project("h1", "injEnergy");
tree->Project("h2", "nucEnergy");
Hi @JohnyBoy
You almost have it, once you have created the TTreeReader you need to select the branches you want to read, this can be done with TTreeReaderValue:
// Create a TTreeReader for the tree
TTreeReader myReader("summary", f);
// If the branch "injEnergy" contains floats; access them as myInjEnergy.
TTreeReaderValue<Float_t> myInjEnergy(myReader, "injEnergy");
// If the branch "nucEnergy" contains floats, too; access those as myNucEnergy.
TTreeReaderValue<Float_t> myNucEnergy(myReader, "nucEnergy");
Now that we have access to the branch values we can fill the histograms with the values corresponding to each event. The following code fills h1 with values from injEnergy and h2 with values from nucEnergy:
// Loop over all entries of the TTree
while (myReader.Next()) {
// Just access the data as if myInjEnergy and myNucEnergy were iterators
// (note the '*' in front of them):
h1->Fill(*myInjEnergy);
h2->Fill(*myNucEnergy);
}
Alternatively, you can use the new RDataFrame interface rather than TTreeReader, which makes your code less verbose and simpler:
// Create an RDataFrame for the tree "summary" and file "f'.
ROOT::RDataFrame rdf("summary", f);
// There is no need to create an specific reader, you can just create a
// histogram out of each branch (note that histogram models can also be specified):
auto h1 = rdf.Histo1D({"h1", "Initial and final...", 100, 0, 5e+19}, "injEnergy")
auto h2 = rdf.Histo1D({"h2", "h2", 100, -5, 5}, "nucEnergy")
// Rest of your code does not change
Thank you all for your help. I have solved it.
@JohnyBoy
Is it possible that you could post your code with the changes ? I would like to see an example 
Sorry. Here is my code.
// Example displaying two histograms and their ratio.
// Author: Olivier Couet.
#include "TFile.h"
#include "TTree.h"
void ratioplot() {
// open the file
TFile *f = TFile::Open("/home/oscar/File2.root");
if (f == 0)
{
// if we cannot open the file, print an error message and return immediatly
printf("Error: cannot open http://lcg-heppkg.web.cern.ch/lcg-heppkg/ROOT/eventdata.root!\n");
}
TH1F *h1 = new TH1F("h1", "Injection and arrival energy spectrum and their ratio; Energy (eV); Number of particles", 1000, 1e+18, 1e+21);
TH1F *h2 = new TH1F("h2", "h2", 1000, 1e+18, 1e+21);
TTree* tree = (TTree*) f->Get("summary");
// Histogram declarations
tree->Project("h1", "nucEnergy");
tree->Project("h2", "injEnergy");
// Define the Canvas
TCanvas *c = new TCanvas("c", "canvas", 800, 800);
// Upper plot will be in pad1
TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
pad1->SetLogx();
pad1->SetLogy();
pad1->SetBottomMargin(0); // Upper and lower plot are joined
pad1->SetGridx(); // Vertical grid
pad1->Draw(); // Draw the upper pad: pad1
pad1->cd(); // pad1 becomes the current pad
h1->SetStats(0); // No statistics on upper plot
h1->Draw(); // Draw h1
h2->Draw("same"); // Draw h2 on top of h1
// Do not draw the Y axis label on the upper plot and redraw a small
// axis instead, in order to avoid the first label (0) to be clipped.
h1->GetYaxis()->SetLabelSize(0.);
TGaxis *axis = new TGaxis( -5, 20, -5, 220, 20,220,510,"");
axis->SetLabelFont(43); // Absolute font size in pixel (precision 3)
axis->SetLabelSize(15);
axis->Draw();
TLegend *legend = new TLegend(0.55,0.65,0.76,0.82);
legend->AddEntry(h1,"Arrival energy","f");
legend->AddEntry(h2,"Injection energy","f");
legend->Draw();
// lower plot will be in pad
c->cd(); // Go back to the main canvas before defining pad2
TPad *pad2 = new TPad("pad2", "pad2", 0, 0.05, 1, 1);
pad2->SetLogx();
pad2->SetLogy();
pad2->SetTopMargin(0);
pad2->SetBottomMargin(0.2);
pad2->SetGridx(); // vertical grid
pad2->Draw();
pad2->cd(); // pad2 becomes the current pad
// Define the ratio plot
TH1F *h3 = (TH1F*)h1->Clone("h3");
h3->SetLineColor(kBlack);
h3->SetMinimum(0.001); // Define Y ..
h3->SetMaximum(3.6); // .. range
h3->Sumw2();
h3->SetStats(0); // No statistics on lower plot
h3->Divide(h2);
h3->SetMarkerStyle(21);
h3->Draw("HIST"); // Draw the ratio plot
// h1 settings
h1->SetLineColor(kBlue+1);
h1->SetLineWidth(2);
// Y axis h1 plot settings
h1->GetYaxis()->SetTitleSize(20);
h1->GetYaxis()->SetTitleFont(43);
h1->GetYaxis()->SetTitleOffset(1.55);
// h2 settings
h2->SetLineColor(kRed);
h2->SetLineWidth(2);
// Ratio plot (h3) settings
h3->SetTitle(""); // Remove the ratio title
// Y axis ratio plot settings
h3->GetYaxis()->SetTitle("Energy spectrum ratios");
h3->GetYaxis()->SetNdivisions(505);
h3->GetYaxis()->SetTitleSize(20);
h3->GetYaxis()->SetTitleFont(43);
h3->GetYaxis()->SetTitleOffset(1.55);
h3->GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
h3->GetYaxis()->SetLabelSize(15);
// X axis ratio plot settings
h3->GetXaxis()->SetTitleSize(20);
h3->GetXaxis()->SetTitleFont(43);
h3->GetXaxis()->SetTitleOffset(2.);
h3->GetXaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
h3->GetXaxis()->SetLabelSize(15);
}
1 Like