ROOT Tree/Ntuple Workflow Comparisons

README.md

There are different ways of handeling ROOT's NTuple / Trees.

I summarize them in the different methods of analyze.C and compare their performances.
generate.C generates real-looking, though fake data needed to run over with analyze.C.

See my blog post on the results and explanation: LINK.

analyze.C

#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TStopwatch.h"
#include <iostream>

// ###############
// The simplest solution of how to fill contents of trees / ntuples into their own histogram: declare every variable explicitly, connect it to a branch, loop over all entries of the tree and, at every point, fill a histogram.
// This is essentiall the same which tree->MakeClass() also generates
// ###############
void explicitVariableDeclaration(TTree * tree) {
	// Create flat float variables for every property that needs to be extracted of tree
	Float_t dPlusMass, dPlusMassSmeared, branchingRatioOfRandomDecayChannel;
	Int_t dPlusPdgCode;
	// Connect branch address to every variable
	tree->SetBranchAddress("dPlusMass", &dPlusMass);
	tree->SetBranchAddress("dPlusMassSmeared", &dPlusMassSmeared);
	tree->SetBranchAddress("branchingRatioOfRandomDecayChannel", &branchingRatioOfRandomDecayChannel);
	tree->SetBranchAddress("dPlusPdgCode", &dPlusPdgCode);

	Int_t numberOfEntries = tree->GetEntries();

	// Create histogram for every variable
	TH1F * histDPlusMass = new TH1F("histDPlusMass", "D^{+} m", 100, 1.7, 2.0);
	TH1F * histDPlusMassSmeared = new TH1F("histDPlusMassSmeared", "D^{+} m (smeared)", 100, 1.5, 2.2);
	TH1F * histDPlusMassDiff = new TH1F("histDPlusMassDiff", "D^{+}: m - m_{smeared}", 100, -0.5, 0.5);
	TH1F * histBranchingRatioOfRandomDecayChannel = new TH1F("histBranchingRatioOfRandomDecayChannel", "Branching Ratio of a random Decay Channel", 100, 0, 0.12);
	TH1F * histDPlusPdgCode = new TH1F("histDPlusPdgCode", "PDG Code of D^{+}", 10, 408, 418);

	// Loop over all entries
	for (int i = 0; i < numberOfEntries; ++i) {
		tree->GetEntry(i); // Get each entry and fill the connected local variables

		// Fill variables into histograms
		histDPlusMass->Fill(dPlusMass);
		histDPlusMassSmeared->Fill(dPlusMassSmeared);
		histDPlusMassDiff->Fill(dPlusMass - dPlusMassSmeared);
		histBranchingRatioOfRandomDecayChannel->Fill(branchingRatioOfRandomDecayChannel);
		histDPlusPdgCode->Fill(dPlusPdgCode);
	}
	// histDPlusMassSmeared->Draw("HIST"); // Commented to not include this in timing
	delete histDPlusMass;
	delete histDPlusMassSmeared;
	delete histDPlusMassDiff;
	delete histBranchingRatioOfRandomDecayChannel;
	delete histDPlusPdgCode;
}

// ###############
// Similar to the first approach, but some more structure - we're using C/C++ after all!
// ###############
// Create a container for all the information in the tree
struct dMesonInfo {
	// For this example, this might only be four, but you can work here with more parameters and even sub-structs and stuff (inheritance!)
	Float_t mass, massSmeared, branchingRatio;
	Int_t pdgCode;
};
// Why not make a function which extracts the boring part of setting branch addresses? Exactly.
void setBranchAddresses(TTree * tree, dMesonInfo & dMeson) {
	TString baseString = "dPlus";
	tree->SetBranchAddress(baseString + "Mass", &(dMeson.mass));
	tree->SetBranchAddress(baseString + "MassSmeared", &(dMeson.massSmeared));
	tree->SetBranchAddress(baseString + "PdgCode", &(dMeson.pdgCode));
	tree->SetBranchAddress("branchingRatioOfRandomDecayChannel", &(dMeson.branchingRatio));
}
void structAndTree(TTree * tree) {
	dMesonInfo dPlus; // Much condensed!
	setBranchAddresses(tree, dPlus); // Very wow!

	Int_t numberOfEntries = tree->GetEntries();

	TH1F * histDPlusMass = new TH1F("histDPlusMass", "D^{+} m", 100, 1.7, 2.0);
	TH1F * histDPlusMassSmeared = new TH1F("histDPlusMassSmeared", "D^{+} m (smeared)", 100, 1.5, 2.2);
	TH1F * histDPlusMassDiff = new TH1F("histDPlusMassDiff", "D^{+}: m - m_{smeared}", 100, -0.5, 0.5);
	TH1F * histBranchingRatioOfRandomDecayChannel = new TH1F("histBranchingRatioOfRandomDecayChannel", "Branching Ratio of a random Decay Channel", 100, 0, 0.12);
	TH1F * histDPlusPdgCode = new TH1F("histDPlusPdgCode", "PDG Code of D^{+}", 10, 408, 418);

	// Loop over all entries
	for (int i = 0; i < numberOfEntries; ++i) {
		tree->GetEntry(i); // Get each entry and fill the connected local variables

		// Fill variables into histograms, this time you can access all parameters neatly as members of the dPlus object
		histDPlusMass->Fill(dPlus.mass);
		histDPlusMassSmeared->Fill(dPlus.massSmeared);
		histDPlusMassDiff->Fill(dPlus.mass - dPlus.massSmeared);
		histBranchingRatioOfRandomDecayChannel->Fill(dPlus.branchingRatio);
		histDPlusPdgCode->Fill(dPlus.pdgCode);
	}
	// histDPlusMassSmeared->Draw("HIST"); // Commented to not include this in timing
	delete histDPlusMass;
	delete histDPlusMassSmeared;
	delete histDPlusMassDiff;
	delete histBranchingRatioOfRandomDecayChannel;
	delete histDPlusPdgCode;
}

// ###############
// ROOT has some built-in functions to shortcut what we did above...
// ###############
void projectOntoHistograms(TTree * tree) {
	// Create histograms as usual
	TH1F * histDPlusMass = new TH1F("histDPlusMass", "D^{+} m", 100, 1.7, 2.0);
	TH1F * histDPlusMassSmeared = new TH1F("histDPlusMassSmeared", "D^{+} m (smeared)", 100, 1.5, 2.2);
	TH1F * histDPlusMassDiff = new TH1F("histDPlusMassDiff", "D^{+}: m - m_{smeared}", 100, -0.5, 0.5);
	TH1F * histBranchingRatioOfRandomDecayChannel = new TH1F("histBranchingRatioOfRandomDecayChannel", "Branching Ratio of a random Decay Channel", 100, 0, 0.12);
	TH1F * histDPlusPdgCode = new TH1F("histDPlusPdgCode", "PDG Code of D^{+}", 10, 408, 418);

	// Now call Project to project onto a histogram a branch from the ROOT tree
	tree->Project("histDPlusMass", "dPlusMass");
	tree->Project("histDPlusMassSmeared", "dPlusMassSmeared");
	tree->Project("histDPlusMassDiff", "dPlusMass - dPlusMassSmeared"); // Note: You can use the usual ROOT maths in-string operators
	tree->Project("histBranchingRatioOfRandomDecayChannel", "branchingRatioOfRandomDecayChannel");
	tree->Project("histDPlusPdgCode", "dPlusPdgCode");

	// histDPlusMassSmeared->Draw("HIST"); // Commented to not include this in timing
	delete histDPlusMass;
	delete histDPlusMassSmeared;
	delete histDPlusMassDiff;
	delete histBranchingRatioOfRandomDecayChannel;
	delete histDPlusPdgCode;
}
// ###############
// Another way to make use of a built-in ROOT function, although less practical (IMHO), as it's running in-string
// ###############
void doubleChevronToHistogram(TTree * tree) {
		// Create histograms as usual
	TH1F * histDPlusMass = new TH1F("histDPlusMass", "D^{+} m", 100, 1.7, 2.0);
	TH1F * histDPlusMassSmeared = new TH1F("histDPlusMassSmeared", "D^{+} m (smeared)", 100, 1.5, 2.2);
	TH1F * histDPlusMassDiff = new TH1F("histDPlusMassDiff", "D^{+}: m - m_{smeared}", 100, -0.5, 0.5);
	TH1F * histBranchingRatioOfRandomDecayChannel = new TH1F("histBranchingRatioOfRandomDecayChannel", "Branching Ratio of a random Decay Channel", 100, 0, 0.12);
	TH1F * histDPlusPdgCode = new TH1F("histDPlusPdgCode", "PDG Code of D^{+}", 10, 408, 418);

	tree->Draw("dPlusMass >> histDPlusMass", "", "goff");
	tree->Draw("dPlusMassSmeared >> histDPlusMassSmeared", "", "goff");
	tree->Draw("dPlusMass - dPlusMassSmeared >> histDPlusMassDiff", "", "goff");
	tree->Draw("branchingRatioOfRandomDecayChannel >> histBranchingRatioOfRandomDecayChannel", "", "goff");
	tree->Draw("dPlusPdgCode >> histDPlusPdgCode", "", "goff");

	// histDPlusMassDiff->Draw("HIST"); // Commented to not include this in timing
	delete histDPlusMass;
	delete histDPlusMassSmeared;
	delete histDPlusMassDiff;
	delete histBranchingRatioOfRandomDecayChannel;
	delete histDPlusPdgCode;
}

// ###############
// Main routine
// ###############
void analyze() {
	TFile * file = new TFile("ntuple.root", "OPEN");
	TTree * tree = (TTree*)file->Get("tree");
	
	TStopwatch timer;

	timer.Start();
	explicitVariableDeclaration(tree);
	timer.Stop();
	std::cout << "## Explicit variables:" << std::endl;
	std::cout << "   Real time: " << timer.RealTime() << ", CPU time " << timer.CpuTime() << std::endl;

	timer.Start();
	structAndTree(tree);
	timer.Stop();
	std::cout << "## Variables, but as a struct:" << std::endl;
	std::cout << "   Real time: " << timer.RealTime() << ", CPU time " << timer.CpuTime() << std::endl;

	timer.Start();
	projectOntoHistograms(tree);
	timer.Stop();
	std::cout << "## Project-ed onto histogram:" << std::endl;
	std::cout << "   Real time: " << timer.RealTime() << ", CPU time " << timer.CpuTime() << std::endl;

	timer.Start();
	doubleChevronToHistogram(tree);
	timer.Stop();
	std::cout << "## Chevron project:" << std::endl;
	std::cout << "   Real time: " << timer.RealTime() << ", CPU time " << timer.CpuTime() << std::endl;
}

int main() {
	analyze();
	return 0;
}

generate.C

#include "TFile.h"
#include "TTree.h"
#include "TRandom2.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"
#include "TDecayChannel.h"


void generate(Int_t numberOfEvents = 100000) {
	// ###############
	// Create infrastructure
	// ###############
	TFile * file = new TFile("ntuple.root", "RECREATE"); // Create File
	TTree * tree = new TTree("tree", "Testing Data."); // Create Tree

	TRandom2 * randomNumberGenerator = new TRandom2(42); // Random generator
	TDatabasePDG * pdgDb = new TDatabasePDG(); // Create instance of ROOTs PDG database
	TParticlePDG * dPlus = pdgDb->GetParticle("D+"); // Get some Particle

	// ###############
	// Create branches
	// ###############
	Float_t dPlusMass, dPlusMassSmeared, branchingRatioOfRandomDecayChannel;
	Int_t dPlusPdgCode;

	tree->Branch("dPlusMass", &dPlusMass, "dPlusMass/F");
	tree->Branch("dPlusMassSmeared", &dPlusMassSmeared, "dPlusMassSmeared/F");
	tree->Branch("dPlusPdgCode", &dPlusPdgCode, "dPlusPdgCode/I");
	tree->Branch("branchingRatioOfRandomDecayChannel", &branchingRatioOfRandomDecayChannel, "branchingRatioOfRandomDecayChannel/F");

	// ###############
	// Now create some data
	// ###############
	dPlusMass = dPlus->Mass();
	dPlusPdgCode = dPlus->PdgCode();
	Int_t numberOfDecayChannels = dPlus->NDecayChannels();
	TDecayChannel * randomDecayChannel;

	for (unsigned int i = 0; i < numberOfEvents; i++) {
		dPlusMassSmeared = randomNumberGenerator->Gaus(dPlusMass, 0.1); // Just to get some randomness

		Int_t numberOfRandomDecayChannel = randomNumberGenerator->Uniform(numberOfDecayChannels);
		randomDecayChannel = dPlus->DecayChannel(numberOfRandomDecayChannel);
		branchingRatioOfRandomDecayChannel = randomDecayChannel->BranchingRatio();

		tree->Fill(); // This fills all values to their connected branches
	}

	// ###############
	// Print, save, and close everything
	// ###############
	tree->Print();
	file->Write();
	file->Close();

}
int main() {
	generate();

	return 0;
}