kiesel/TreeWriter/treeWriter.cc

#include<iostream>

#include "TFile.h"
#include "TTree.h"
#include "TString.h"

#include "SusyEvent.h"
#include "TreeObjects.h"


class TreeWriter {
        public :
                TreeWriter(TString inputName, TString outputName );
                virtual ~TreeWriter();
                virtual void Loop();

                void SetProcessNEvents(int nEvents) { processNEvents = nEvents; }
                void SetReportEvents(int nEvents) { reportEvery = nEvents; }

                TFile *inputFile;
                TTree *inputTree;
                susy::Event *event;

                TFile *outFile;
                TTree *tree;

        private:
                int processNEvents; // number of events to be processed
                int reportEvery;
                int loggingVerbosity;

                // variables which will be stored in the tree
                std::vector<tree::Photon> photon;
                std::vector<tree::Jet> jet;
                float met;
                int nVertex;
                float weight;
                float leadingPhotonPt;
};


TreeWriter::TreeWriter(TString inputName, TString outputName) {
        // read the input file
        inputFile = new TFile( inputName, "read" );
        inputTree = (TTree*) inputFile->Get("susyTree");
        event = new susy::Event;
        inputTree->SetBranchAddress("susyEvent", &event);

        // open the output file
        outFile = new TFile( outputName, "recreate" );
        tree = new TTree("susyTree","Tree for single photon analysis");

        // set default parameter
        processNEvents = -1;
        reportEvery = 1000;
        loggingVerbosity = 0;

}

TreeWriter::~TreeWriter() {
        if (!inputTree) return;
        delete inputTree->GetCurrentFile();
}

void TreeWriter::Loop() {
        // here the event loop is implemented and the tree is filled
        if (inputTree == 0) return;

        // get number of events to be proceeded
        Long64_t nentries = inputTree->GetEntries();
        if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;

        if( loggingVerbosity > 0 )
                std::cout << "Processing " << processNEvents << " ouf of "
                        << nentries << " events. " << std::endl;

        tree::Photon *thisphoton = new tree::Photon();
        tree::Jet *thisjet = new tree::Jet();

        tree->Branch("photon", &photon);
        tree->Branch("jet", &jet);
        tree->Branch("met", &met, "met/F");
        tree->Branch("nVertex", &nVertex, "nVertex/I");
        tree->Branch("weigth", &weight, "weight/F");
        tree->Branch("photonPt", &leadingPhotonPt, "photonPt/F");


        for (long jentry=0; jentry < processNEvents; jentry++) {
                if (jentry%reportEvery==0)
                        std::cout << jentry << " / " << processNEvents << std :: endl;
                inputTree->GetEntry(jentry);

                photon.clear();
                jet.clear();

                // photons
                std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
                for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
                                it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
                        if( ! it->isEB() )
                                continue;
                        thisphoton->pt = it->momentum.Et();
                        if( thisphoton->pt < 80 )
                                continue;
                        thisphoton->eta = it->momentum.Eta();
                        thisphoton->chargedIso = it->chargedHadronIso;
                        thisphoton->neutralIso = it->neutralHadronIso;
                        thisphoton->photonIso = it->photonIso;
                        if ( it->r9 > 1 ) // if == 1 ?
                                continue;
                        thisphoton->r9 = it->r9;
                        thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
                        thisphoton->hadTowOverEm = it->hadTowOverEm;
                        thisphoton->pixelseed = it->nPixelSeeds;
                        photon.push_back( *thisphoton );
                }
                if( photon.size() == 0 )
                        continue;
                std::sort( photon.begin(), photon.end(), tree::EtGreater);
                //leadingPhotonPt = photon.at(0)->pt;


                // jets
                std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
                if(pfJets_it == event->pfJets.end()){
                        if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
                } else {

                        susy::PFJetCollection& jetColl = pfJets_it->second;

                        for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
                                        it != jetColl.end(); it++) {
                                std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
                                if (s_it == it->jecScaleFactors.end()) {
                                        std::cout << "JEC is not available for this jet!!!" << std::endl;
                                        continue;
                                }
                                float scale = s_it->second;
                                TLorentzVector corrP4 = scale * it->momentum;

                                if(std::abs(corrP4.Eta()) > 3.0) continue;
                                thisjet->pt = corrP4.Et();
                                thisjet->eta = corrP4.Eta();
                                jet.push_back( *thisjet );
                        }// for jet
                }// if, else
                if( jet.size() == 0 )
                        std::cout << "error, no jets found " << std::endl;
                else
                        std::sort( jet.begin(), jet.end(), tree::EtGreater);

                // met
                std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
                susy::MET* metobj = &(met_it->second);
                met = metobj->met();

                // vertices

                nVertex = event->vertices.size();
                weight = 1;

                tree->Fill();
        } // for jentry


        tree->Write();
        outFile->cd();
        outFile->Write();
        outFile->Close();
}

int main(int argc, char** argv) {
        TreeWriter *tw = new TreeWriter("qcd-1000-nTuple-test.root", "myQCDTree.root");
        tw->SetProcessNEvents(10);
        tw->Loop();
}

Revision:	1.3
Committed:	Thu Mar 21 09:01:14 2013 UTC (12 years, 1 month ago) by kiesel
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+6 -28 lines
Log Message:	Added own tree objects to treewriter
#	User	Rev	Content
1	kiesel	1.2	#include<iostream>
2
3	kiesel	1.1	#include "TFile.h"
4			#include "TTree.h"
5			#include "TString.h"
6
7			#include "SusyEvent.h"
8	kiesel	1.3	#include "TreeObjects.h"
9	kiesel	1.1
10
11			class TreeWriter {
12			public :
13			TreeWriter(TString inputName, TString outputName );
14			virtual ~TreeWriter();
15			virtual void Loop();
16
17			void SetProcessNEvents(int nEvents) { processNEvents = nEvents; }
18			void SetReportEvents(int nEvents) { reportEvery = nEvents; }
19
20			TFile *inputFile;
21			TTree *inputTree;
22			susy::Event *event;
23
24			TFile *outFile;
25			TTree *tree;
26
27			private:
28			int processNEvents; // number of events to be processed
29			int reportEvery;
30			int loggingVerbosity;
31
32			// variables which will be stored in the tree
33			std::vector<tree::Photon> photon;
34			std::vector<tree::Jet> jet;
35			float met;
36			int nVertex;
37			float weight;
38	kiesel	1.3	float leadingPhotonPt;
39	kiesel	1.1	};
40
41
42			TreeWriter::TreeWriter(TString inputName, TString outputName) {
43			// read the input file
44			inputFile = new TFile( inputName, "read" );
45			inputTree = (TTree*) inputFile->Get("susyTree");
46			event = new susy::Event;
47			inputTree->SetBranchAddress("susyEvent", &event);
48
49			// open the output file
50			outFile = new TFile( outputName, "recreate" );
51			tree = new TTree("susyTree","Tree for single photon analysis");
52
53			// set default parameter
54			processNEvents = -1;
55			reportEvery = 1000;
56			loggingVerbosity = 0;
57
58			}
59
60			TreeWriter::~TreeWriter() {
61			if (!inputTree) return;
62			delete inputTree->GetCurrentFile();
63			}
64
65			void TreeWriter::Loop() {
66			// here the event loop is implemented and the tree is filled
67			if (inputTree == 0) return;
68
69			// get number of events to be proceeded
70			Long64_t nentries = inputTree->GetEntries();
71			if(processNEvents <= 0 \|\| processNEvents > nentries) processNEvents = nentries;
72
73			if( loggingVerbosity > 0 )
74			std::cout << "Processing " << processNEvents << " ouf of "
75			<< nentries << " events. " << std::endl;
76
77			tree::Photon *thisphoton = new tree::Photon();
78			tree::Jet *thisjet = new tree::Jet();
79
80			tree->Branch("photon", &photon);
81			tree->Branch("jet", &jet);
82			tree->Branch("met", &met, "met/F");
83			tree->Branch("nVertex", &nVertex, "nVertex/I");
84			tree->Branch("weigth", &weight, "weight/F");
85	kiesel	1.3	tree->Branch("photonPt", &leadingPhotonPt, "photonPt/F");
86	kiesel	1.1
87
88			for (long jentry=0; jentry < processNEvents; jentry++) {
89			if (jentry%reportEvery==0)
90			std::cout << jentry << " / " << processNEvents << std :: endl;
91			inputTree->GetEntry(jentry);
92
93			photon.clear();
94			jet.clear();
95
96			// photons
97			std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
98			for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
99			it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
100			if( ! it->isEB() )
101			continue;
102			thisphoton->pt = it->momentum.Et();
103			if( thisphoton->pt < 80 )
104			continue;
105			thisphoton->eta = it->momentum.Eta();
106			thisphoton->chargedIso = it->chargedHadronIso;
107			thisphoton->neutralIso = it->neutralHadronIso;
108			thisphoton->photonIso = it->photonIso;
109			if ( it->r9 > 1 ) // if == 1 ?
110			continue;
111			thisphoton->r9 = it->r9;
112			thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
113			thisphoton->hadTowOverEm = it->hadTowOverEm;
114			thisphoton->pixelseed = it->nPixelSeeds;
115			photon.push_back( *thisphoton );
116			}
117			if( photon.size() == 0 )
118			continue;
119			std::sort( photon.begin(), photon.end(), tree::EtGreater);
120	kiesel	1.3	//leadingPhotonPt = photon.at(0)->pt;
121	kiesel	1.1
122
123			// jets
124			std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
125			if(pfJets_it == event->pfJets.end()){
126			if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
127			} else {
128
129			susy::PFJetCollection& jetColl = pfJets_it->second;
130
131			for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
132			it != jetColl.end(); it++) {
133			std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
134			if (s_it == it->jecScaleFactors.end()) {
135			std::cout << "JEC is not available for this jet!!!" << std::endl;
136			continue;
137			}
138			float scale = s_it->second;
139			TLorentzVector corrP4 = scale * it->momentum;
140
141			if(std::abs(corrP4.Eta()) > 3.0) continue;
142			thisjet->pt = corrP4.Et();
143			thisjet->eta = corrP4.Eta();
144			jet.push_back( *thisjet );
145			}// for jet
146			}// if, else
147			if( jet.size() == 0 )
148			std::cout << "error, no jets found " << std::endl;
149			else
150			std::sort( jet.begin(), jet.end(), tree::EtGreater);
151
152			// met
153			std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
154			susy::MET* metobj = &(met_it->second);
155			met = metobj->met();
156
157			// vertices
158
159			nVertex = event->vertices.size();
160			weight = 1;
161
162			tree->Fill();
163			} // for jentry
164
165
166
167			tree->Write();
168			outFile->cd();
169			outFile->Write();
170			outFile->Close();
171			}
172	kiesel	1.3
173	kiesel	1.1	int main(int argc, char** argv) {
174			TreeWriter *tw = new TreeWriter("qcd-1000-nTuple-test.root", "myQCDTree.root");
175			tw->SetProcessNEvents(10);
176			tw->Loop();
177			}
178	kiesel	1.3