Reading vector<double> from tree

ROOT
1

Hello,

I’m trying to learn how to use vectors in trees. I found some exemplary code and I tried to introduce it into my code. I was successful in writing vectors into trees and than saving it into root file. But I have problems reading tree.

Below is all code I’m using (although it is part of something bigger). I compile it and run as standalone executable.

void Sample::ElectronPreselection( const string & dir ) const {
	cout << "Entering ElectronPreselection() ..." << endl;
	gROOT->ProcessLine("#include <vector>");
	TFile * file = new TFile( GetPath().c_str() );
	TTree * tree = ( TTree * ) file->Get( "HZZ2l2nuAnalysis" );
	Event ev( tree );

	string outputFile = dir + GetName() + "_preselected.root";
	cout << outputFile << endl;
	TFile * out = new TFile( outputFile.c_str(), "recreate" );
	unsigned run;
	unsigned lumi;
	unsigned event;
	double pfmet;
	int nele;
	int nmu;
	int nsoftmu;
	double l1pt;
	double l1eta;
	double l2pt;
	double l2eta;
	double zmass;
	double zpt;
	double zeta;
	double mt;
	int njet;
	double maxJetBTag;
	double minDeltaPhiJetMet;
	int nvtx;
	int ni;

	TTree * smallTree = new TTree("HZZ2l2nuAnalysis", "HZZ2l2nu Analysis Tree");
	smallTree->Branch( "Run", &run, "Run/i" );
	smallTree->Branch( "Lumi", &lumi, "Lumi/i" );
	smallTree->Branch( "Event", &event, "Event/i" );
	smallTree->Branch( "PFMET", &pfmet, "PFMET/D" );
	smallTree->Branch( "NELE", &nele, "NELE/I" );
	smallTree->Branch( "NMU", &nmu, "NMU/I" );
	smallTree->Branch( "NSOFTMU", &nsoftmu, "NSOFTMU/I" );
	smallTree->Branch( "L1PT", &l1pt, "L1PT/D" );
	smallTree->Branch( "L1ETA", &l1eta, "L1ETA/D" );
	smallTree->Branch( "L2PT", &l2pt, "L2PT/D" );
	smallTree->Branch( "L2ETA", &l2eta, "L2ETA/D" );
	smallTree->Branch( "ZMASS", &zmass, "ZMASS/D" );
	smallTree->Branch( "ZPT", &zpt, "ZPT/D" );
	smallTree->Branch( "ZETA", &zeta, "ZETA/D" );
	smallTree->Branch( "MT", &mt, "MT/D" );
	smallTree->Branch( "NJET", &njet, "NJET/I" );
	smallTree->Branch( "MAXJETBTAG", &maxJetBTag, "MAXJETBTAG/D" );
	smallTree->Branch( "MINDPJETMET", &minDeltaPhiJetMet, "MINDPJETMET/D" );
	smallTree->Branch( "NVTX", &nvtx, "NVTX/I" );
	smallTree->Branch( "nInter" , &ni, "nInter/I" );

	unsigned long nentries = tree->GetEntries();
	for ( unsigned long i = 0; i < nentries; i++ ) {
		cout << "TEST 0, Loop " << i << endl;
		tree->GetEntry( i );
		cout << "TEST 1, Loop " << i << endl;

		run = ev.getUnsignedVariableValue("Run");
		lumi = ev.getUnsignedVariableValue("LumiSection");
		event = ev.getUnsignedVariableValue("Event");
		cout << "Event: " << event << endl;
		pfmet = ev.getDoubleVariableValue("PFMET");

		vector<unsigned> electrons20;
		selectElectrons( ev, electrons20, 20 );
		if ( electrons20.size() < 2 )
			continue;

		vector<unsigned> electrons10;
		selectElectrons( ev, electrons10 );
		nele = electrons10.size();

		vector<unsigned> muons10;
		selectMuons( ev, muons10 );
		nmu = muons10.size();

		vector<unsigned> softMuons;
		selectSoftMuons( ev, softMuons );
		nsoftmu = softMuons.size();

		const int * lepCharge = ev.getIntLAArrayAdr("Electrons_CHARGE", 0);
		if ( lepCharge[0] == lepCharge[1] )
			continue;

		const double * lepE = ev.getDoubleLAArrayAdr("Electrons_E", 0);
		const double * lepPt = ev.getDoubleLAArrayAdr("Electrons_PT", 0);
		const double * lepEta = ev.getDoubleLAArrayAdr("Electrons_ETA", 0);
		const double * lepPhi = ev.getDoubleLAArrayAdr("Electrons_PHI", 0);
		TLorentzVector lep1;
		lep1.SetPtEtaPhiE( lepPt[0], lepEta[0], lepPhi[0], lepE[0]);
		l1pt = lep1.Pt();
		l1eta = lep1.Eta();

		TLorentzVector lep2;
		lep2.SetPtEtaPhiE( lepPt[1], lepEta[1], lepPhi[1], lepE[1]);
		l2pt = lep2.Pt();
		l2eta = lep2.Eta();

		TLorentzVector Zcand = lep1 + lep2;
		zmass = Zcand.M();
		zpt = Zcand.Pt();
		zeta = Zcand.Eta();
		if (zmass < 40)
			continue;

		double metphi = ev.getDoubleVariableValue("PFMET_PHI");
		double metPx = pfmet * cos( metphi );
		double metPy = pfmet * sin( metphi );
		double px = metPx + Zcand.Px();
		double py = metPy + Zcand.Py();
		double pt2 = px * px + py * py;
		double e = sqrt(zpt * zpt + zmass * zmass) + sqrt(pfmet * pfmet + zmass * zmass);
		double mt2 = e * e - pt2;
		mt = (mt2 > 0) ? sqrt(mt2) : 0;

		vector<TLorentzVector> leptons;
		leptons.push_back( lep1 );
		leptons.push_back( lep2 );
		vector<unsigned> jets;
		selectJets( ev, jets, leptons );

		njet = jets.size();

		const double * btag = ev.getDoubleLAArrayAdr("Jets_BTAG_TrackCountingHighEff", 0);
		const double * jetEta = ev.getDoubleLAArrayAdr("Jets_ETA", 0);
		maxJetBTag = -999;
		for ( int j = 0; j < njet; ++j ) {
			if ( btag[j] > maxJetBTag && fabs(jetEta[j]) < 2.4 )
				maxJetBTag = btag[j];
		}

		const double * delPhiJetMet = ev.getDoubleLAArrayAdr("Jets_DELTAPHIMET", 0);
		minDeltaPhiJetMet = 999;
		for ( int j = 0; j < njet; ++j ) {
			if ( delPhiJetMet[j] > -1 && delPhiJetMet[j] < minDeltaPhiJetMet )
				minDeltaPhiJetMet = delPhiJetMet[j];
		}

		nvtx = ev.getIntVariableValue("NVTX");

		ni = nvtx;
		const int * bx = ev.getInt5ArrayAdr("PUP_BunchCrossing", 0);
		const int * ninter = ev.getInt5ArrayAdr("PUP_NumInteractions", 0);
		for ( int j = 0; j < 5; ++j )
			if ( bx[j] == 0 )
				ni = ninter[j];
		
		smallTree->Fill();
	}
	cout << smallTree->GetEntries() << endl;
	delete file;
	smallTree->Write("", TObject::kOverwrite);
	delete smallTree;
	delete out;
}

Code breaks ( *** Break *** segmentation violation) on line
tree->GetEntry( i );

I also add definition and constructor of Event class:

#ifndef EVENT_H
#define EVENT_H

// ROOT Libraries
#include <TTree.h>
// Standard Libraries
#include <string>
#include <utility>
#include <vector>

class Event {
	public :
		const static int nLA = 25;
		//const static int nJET = 40;
		enum VARTYPE { UNSIGN, INT, DOUBLE, UNSIGN10, DOUBLE5, DOUBLE3, INT5, DOUBLELA };
		Event(TTree *tree);
		~Event();
		unsigned getUnsignedVariableValue(const std::string & name) const;
		unsigned getUnsigned10ArrayValue(const std::string & name, unsigned n) const;
		const unsigned * getUnsignedVariableAdr(const std::string & name) const;
		const unsigned * getUnsigned10ArrayAdr(const std::string & name, unsigned n) const;
		int getIntVariableValue(const std::string & name) const;
		int getInt5ArrayValue(const std::string & name, unsigned n) const;
		const int * getIntVariableAdr(const std::string & name) const;
		const int * getInt5ArrayAdr(const std::string & name, unsigned n) const;
		const int * getIntLAArrayAdr(const std::string & name, unsigned n) const;
		double getDoubleVariableValue(const std::string & name) const;
		double getDouble3ArrayValue(const std::string & name, unsigned n) const;
		double getDouble5ArrayValue(const std::string & name, unsigned n) const;
		double getDouble10ArrayValue(const std::string & name, unsigned n) const;
		double getDoubleLAArrayValue(const std::string & name, unsigned n) const;
		const double * getDoubleVariableAdr(const std::string & name) const;
		const double * getDouble3ArrayAdr(const std::string & name, unsigned n) const;
		const double * getDouble5ArrayAdr(const std::string & name, unsigned n) const;
		const double * getDouble10ArrayAdr(const std::string & name, unsigned n) const;
		const double * getDoubleLAArrayAdr(const std::string & name, unsigned n) const;
		bool findVariable(const std::string & name, VARTYPE & vt) const;
		double getValue(const std::string & name) const;
	private :
		std::vector<std::pair<std::string, unsigned *> > unsignedTypes;
		std::vector<std::pair<std::string, int *> > intTypes;
		std::vector<std::pair<std::string, double *> > doubleTypes;
		std::vector<std::pair<std::string, unsigned *> > unsigned10Types;
		std::vector<std::pair<std::string, int *> > int5Types;
		std::vector<std::pair<std::string, int *> > intLATypes;
		std::vector<std::pair<std::string, double *> > double3Types;
		std::vector<std::pair<std::string, double *> > double5Types;
		std::vector<std::pair<std::string, double *> > double10Types;
		std::vector<std::pair<std::string, double *> > doubleLATypes;
		std::vector<std::pair<std::string, std::vector<double> *> > VectorsDouble;

};

#endif // EVENT_H
Event::Event( TTree *tree )
{
	gROOT->ProcessLine("#include <vector>");
	TObjArray * treeLeaves = tree->GetListOfLeaves();
	int nEntries = treeLeaves->GetEntries();
	TObjArray::Iterator_t iter( treeLeaves );
	for( int i = 0; i < nEntries; ++i ) {
		iter.Next();
		TLeaf * leafPointer = dynamic_cast<TLeaf *>( *iter );
		//cout << leafPointer->GetTypeName() << endl;
		//cout << leafPointer->GetLen() << endl;
		//cout << leafPointer->GetBranch()->GetName() << endl;
		if( leafPointer ) {
			if ( leafPointer->GetTypeName() == string("UInt_t") && leafPointer->GetLen() == 1 ) {
				pair<string, unsigned *> temp(leafPointer->GetBranch()->GetName(), new unsigned);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				unsignedTypes.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("Int_t") && leafPointer->GetLen() == 1 ) {
				pair<string, int *> temp(leafPointer->GetBranch()->GetName(), new int);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				intTypes.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("Double_t") && leafPointer->GetLen() == 1 ) {
				pair<string, double *> temp(leafPointer->GetBranch()->GetName(), new double);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				doubleTypes.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("UInt_t") && leafPointer->GetLen() == 10 ) {
				pair<string, unsigned *> temp(leafPointer->GetBranch()->GetName(), new unsigned[10]);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				unsigned10Types.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("Int_t") && leafPointer->GetLen() == 5 ) {
				pair<string, int *> temp(leafPointer->GetBranch()->GetName(), new int[5]);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				int5Types.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("Int_t") && leafPointer->GetLen() == Event::nLA ) {
				pair<string, int *> temp(leafPointer->GetBranch()->GetName(), new int[Event::nLA]);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				intLATypes.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("Double_t") && leafPointer->GetLen() == 3 ) {
				pair<string, double *> temp(leafPointer->GetBranch()->GetName(), new double[3]);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				double3Types.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("Double_t") && leafPointer->GetLen() == 5 ) {
				pair<string, double *> temp(leafPointer->GetBranch()->GetName(), new double[5]);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				double5Types.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("Double_t") && leafPointer->GetLen() == 10 ) {
				pair<string, double *> temp(leafPointer->GetBranch()->GetName(), new double[10]);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				double10Types.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("Double_t") && leafPointer->GetLen() == Event::nLA ) {
				pair<string, double *> temp(leafPointer->GetBranch()->GetName(), new double[Event::nLA]);
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), temp.second );
				doubleLATypes.push_back(temp);
			} else if ( leafPointer->GetTypeName() == string("vector<double>") && leafPointer->GetLen() == 1 ) {
				pair<string, vector<double> *> temp(leafPointer->GetBranch()->GetName(), 0 );
				temp.second = 0;
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), &(temp.second) );
				VectorsDouble.push_back(temp);
			} else {
				cout << "ERROR: Don't know what to do with this type: " << leafPointer->GetTypeName() << '[' << leafPointer->GetLen() << "]!" << endl;
				exit (EXIT_FAILURE);
			}
		}
	}
}

I would be very grateful on understanding what might be wrong.

2

Hi,

pair<string, vector<double> *> temp(leafPointer->GetBranch()->GetName(), 0 ); temp.second = 0; tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), &(temp.second) ); The address passed to SetBranchAddress is recorded in the TTree. However in the code above the address being passed is the address of a local variable which goes aways at the end of the code block (and thus GetEntry will use this address and cause a segmentation fault).

Cheers,
Philippe.

3

Philippe,

I was also thinking that this might be the issue, so I changed code to:

            } else if ( leafPointer->GetTypeName() == string("vector<double>") && leafPointer->GetLen() == 1 ) {
                    pair<string, vector<double> *> temp(leafPointer->GetBranch()->GetName(), 0 );
                    temp.second = 0;
                    VectorsDouble.push_back(temp);
                    tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), &(VectorsDouble[VectorsDouble.size() - 1].second) );

but this also doesn’t work.

By the way, why in case of custom objects SetBranchAddress needs address of the pointer. That seems a bit strange.

4

[quote]By the way, why in case of custom objects SetBranchAddress needs address of the pointer. That seems a bit strange[/quote]That is because in this case the TTree needs the ability to create the object itself and returns it back to the user (We are planning upgrading SetBranchAddress to also support passing the address of an existing object in the next year).

Cheers,
Philippe.

5

Hi,

[quote]but this also doesn’t work.[/quote]humm … I am not sure. The best thing is to run your example with valgrind to pinpoint the (memory) problem.

Cheers,
Philippe.

6

I tried using valgrind, but I’m not expert and I couldn’t find out reason. If something it looks like error comes from ROOT libraries.

I attach below fully functional code with input sample. I would appreciate if you could look at this.
macro.tgz (514 KB)

7

Hi,

The location of objects in a std::vector is not stable. If the vector needs to grows, it will move the object into a larger area of memory. One possible solution is to call VectorsDouble.reserve(large_enough_value); which, if you are really able to pick a value of ‘large_enough_value’ that is guaranteed to be larger than the maximum size of the vector, will insure the object in the vector are not moved.

Cheers,
Philippe.

8

Philippe,

yesterday evening I was thinking about this problem and I also came up with idea you are suggesting. Yet this is not very robust as you have know how big vector you will need. I found more efficient solution. I give it here in case someone else will encounter similar problem.

			} else if ( leafPointer->GetTypeName() == string("vector<double>") && leafPointer->GetLen() == 1 ) {
				vector<double> ** tmpPtr = new vector<double> *;
				*tmpPtr = 0;
				tree->SetBranchAddress( leafPointer->GetBranch()->GetName(), tmpPtr );
				pair<string, vector<double> **> temp(leafPointer->GetBranch()->GetName(), tmpPtr );
				VectorsDouble.push_back(temp);

P.S. I think that this quite weird “feature” of SetBranchAddress with custom objects should be well documented, since unaware users can have real problems with that.