kiesel/TreeWriter/treeWriter.cc

#include<iostream>
#include<math.h>

#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;

                float getPtFromMatchedJet( susy::Photon, susy::Event );
                float deltaR( TLorentzVector, TLorentzVector );

        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;
                int nElectron;
                float weight;
};

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();
}

float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
        return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
}

float TreeWriter::getPtFromMatchedJet( susy::Photon photon, susy::Event event ) {
        /**
         * \brief Takes jet p_T as photon p_T
         *
         * At first all jets with DeltaR < 0.3 (isolation cone) are searched.
         * If several jets are found, take the one with the minimal pt difference
         * compared to the photon. If no such jets are found, return 0
         * TODO: remove photon matched jet from jet-selection?
         */
        std::vector<susy::PFJet> nearJets;
        nearJets.clear();

        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;
                        float deltaR_ = deltaR(photon.momentum, corrP4 );
                        if (deltaR_ > 0.3) continue;
                        nearJets.push_back( *it );
                }// for jet
        }// if, else

        if ( nearJets.size() == 0 ) {
                std::cout << "No jet with ΔR < .3 found, set p_T to 0" << std::endl;
                return 0;
        }

        float pt = 0;
        float minPtDifferenz = 1E20; // should be very high
        for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
                        it != jetEnd; it++ ) {
                float ptDiff = fabs(photon.momentum.Et() - it->momentum.Et());
                if (  ptDiff < minPtDifferenz ) {
                        minPtDifferenz = ptDiff;
                        pt = it->momentum.Et();
                }
        }

        // testing
        if( nearJets.size() > 1 )
                std::cout << "There are several jets matching to this photon. "
                                        << "Please check if the upper assumtion is reasonable" << std::endl;
        return pt;
}


void TreeWriter::Loop() {
        // only for testing
        bool skim = true;

        // 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("nElectron", &nElectron, "nElectron/I");


        for (long jentry=0; jentry < processNEvents; jentry++) {
                if ( loggingVerbosity>0 && 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() && skim )
                                continue;
                        thisphoton->pt = getPtFromMatchedJet( *it, *event );
                        if( thisphoton->pt < 80 && skim )
                                continue;
                        thisphoton->eta = it->momentum.Eta();
                        thisphoton->phi = it->momentum.Phi();
                        thisphoton->chargedIso = it->chargedHadronIso;
                        thisphoton->neutralIso = it->neutralHadronIso;
                        thisphoton->photonIso = it->photonIso;
                        if ( it->r9 > 1 && skim ) // 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 && skim )
                        continue;
                std::sort( photon.begin(), photon.end(), tree::EtGreater);

                // 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 && skim ) continue;
                                thisjet->pt = corrP4.Et();
                                thisjet->eta = corrP4.Eta();
                                thisjet->phi = corrP4.Phi();
                                jet.push_back( *thisjet );
                        }// for jet
                }// if, else
                if( jet.size() < 2 && skim )
                        continue;
                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();

                // electrons
                std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
                nElectron = eVector.size();

                // vertices

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

                // bjets

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


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

int main(int argc, char** argv) {
        TreeWriter *tw = new TreeWriter("../root-files/susyEvents_qcd_1000-inf_part.root", "myTree.root");
        tw->SetProcessNEvents(-1);
        tw->Loop();
}

Revision:	1.4
Committed:	Mon Mar 25 12:00:06 2013 UTC (12 years, 1 month ago) by kiesel
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+89 -18 lines
Log Message:	Shared objects are ready to be imported by plot tool
#	Content
1	#include<iostream>
2	#include<math.h>
3
4	#include "TFile.h"
5	#include "TTree.h"
6	#include "TString.h"
7
8	#include "SusyEvent.h"
9	#include "TreeObjects.h"
10
11
12	class TreeWriter {
13	public :
14	TreeWriter(TString inputName, TString outputName );
15	virtual ~TreeWriter();
16	virtual void Loop();
17
18	void SetProcessNEvents(int nEvents) { processNEvents = nEvents; }
19	void SetReportEvents(int nEvents) { reportEvery = nEvents; }
20
21	TFile *inputFile;
22	TTree *inputTree;
23	susy::Event *event;
24
25	TFile *outFile;
26	TTree *tree;
27
28	float getPtFromMatchedJet( susy::Photon, susy::Event );
29	float deltaR( TLorentzVector, TLorentzVector );
30
31	private:
32	int processNEvents; // number of events to be processed
33	int reportEvery;
34	int loggingVerbosity;
35
36	// variables which will be stored in the tree
37	std::vector<tree::Photon> photon;
38	std::vector<tree::Jet> jet;
39	float met;
40	int nVertex;
41	int nElectron;
42	float weight;
43	};
44
45	TreeWriter::TreeWriter(TString inputName, TString outputName) {
46	// read the input file
47	inputFile = new TFile( inputName, "read" );
48	inputTree = (TTree*) inputFile->Get("susyTree");
49	event = new susy::Event;
50	inputTree->SetBranchAddress("susyEvent", &event);
51
52	// open the output file
53	outFile = new TFile( outputName, "recreate" );
54	tree = new TTree("susyTree","Tree for single photon analysis");
55
56	// set default parameter
57	processNEvents = -1;
58	reportEvery = 1000;
59	loggingVerbosity = 0;
60
61	}
62
63	TreeWriter::~TreeWriter() {
64	if (!inputTree) return;
65	delete inputTree->GetCurrentFile();
66	}
67
68	float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
69	return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
70	}
71
72	float TreeWriter::getPtFromMatchedJet( susy::Photon photon, susy::Event event ) {
73	/**
74	* \brief Takes jet p_T as photon p_T
75	*
76	* At first all jets with DeltaR < 0.3 (isolation cone) are searched.
77	* If several jets are found, take the one with the minimal pt difference
78	* compared to the photon. If no such jets are found, return 0
79	* TODO: remove photon matched jet from jet-selection?
80	*/
81	std::vector<susy::PFJet> nearJets;
82	nearJets.clear();
83
84	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event.pfJets.find("ak5");
85	if(pfJets_it == event.pfJets.end()){
86	if(event.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
87	} else {
88	susy::PFJetCollection& jetColl = pfJets_it->second;
89	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
90	it != jetColl.end(); it++) {
91	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
92	if (s_it == it->jecScaleFactors.end()) {
93	std::cout << "JEC is not available for this jet!!!" << std::endl;
94	continue;
95	}
96	float scale = s_it->second;
97	TLorentzVector corrP4 = scale * it->momentum;
98	float deltaR_ = deltaR(photon.momentum, corrP4 );
99	if (deltaR_ > 0.3) continue;
100	nearJets.push_back( *it );
101	}// for jet
102	}// if, else
103
104	if ( nearJets.size() == 0 ) {
105	std::cout << "No jet with ΔR < .3 found, set p_T to 0" << std::endl;
106	return 0;
107	}
108
109	float pt = 0;
110	float minPtDifferenz = 1E20; // should be very high
111	for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
112	it != jetEnd; it++ ) {
113	float ptDiff = fabs(photon.momentum.Et() - it->momentum.Et());
114	if ( ptDiff < minPtDifferenz ) {
115	minPtDifferenz = ptDiff;
116	pt = it->momentum.Et();
117	}
118	}
119
120	// testing
121	if( nearJets.size() > 1 )
122	std::cout << "There are several jets matching to this photon. "
123	<< "Please check if the upper assumtion is reasonable" << std::endl;
124	return pt;
125	}
126
127
128	void TreeWriter::Loop() {
129	// only for testing
130	bool skim = true;
131
132	// here the event loop is implemented and the tree is filled
133	if (inputTree == 0) return;
134
135	// get number of events to be proceeded
136	Long64_t nentries = inputTree->GetEntries();
137	if(processNEvents <= 0 \|\| processNEvents > nentries) processNEvents = nentries;
138
139	if( loggingVerbosity > 0 )
140	std::cout << "Processing " << processNEvents << " ouf of "
141	<< nentries << " events. " << std::endl;
142
143	tree::Photon *thisphoton = new tree::Photon();
144	tree::Jet *thisjet = new tree::Jet();
145
146	tree->Branch("photon", &photon);
147	tree->Branch("jet", &jet);
148	tree->Branch("met", &met, "met/F");
149	tree->Branch("nVertex", &nVertex, "nVertex/I");
150	tree->Branch("weigth", &weight, "weight/F");
151	tree->Branch("nElectron", &nElectron, "nElectron/I");
152
153
154	for (long jentry=0; jentry < processNEvents; jentry++) {
155	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
156	std::cout << jentry << " / " << processNEvents << std :: endl;
157	inputTree->GetEntry(jentry);
158
159	photon.clear();
160	jet.clear();
161
162	// photons
163	std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
164	for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
165	it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
166	if( ! it->isEB() && skim )
167	continue;
168	thisphoton->pt = getPtFromMatchedJet( it, event );
169	if( thisphoton->pt < 80 && skim )
170	continue;
171	thisphoton->eta = it->momentum.Eta();
172	thisphoton->phi = it->momentum.Phi();
173	thisphoton->chargedIso = it->chargedHadronIso;
174	thisphoton->neutralIso = it->neutralHadronIso;
175	thisphoton->photonIso = it->photonIso;
176	if ( it->r9 > 1 && skim ) // if == 1 ?
177	continue;
178	thisphoton->r9 = it->r9;
179	thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
180	thisphoton->hadTowOverEm = it->hadTowOverEm;
181	thisphoton->pixelseed = it->nPixelSeeds;
182	photon.push_back( *thisphoton );
183	}
184	if( photon.size() == 0 && skim )
185	continue;
186	std::sort( photon.begin(), photon.end(), tree::EtGreater);
187
188	// jets
189	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
190	if(pfJets_it == event->pfJets.end()){
191	if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
192	} else {
193
194	susy::PFJetCollection& jetColl = pfJets_it->second;
195
196	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
197	it != jetColl.end(); it++) {
198	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
199	if (s_it == it->jecScaleFactors.end()) {
200	std::cout << "JEC is not available for this jet!!!" << std::endl;
201	continue;
202	}
203	float scale = s_it->second;
204	TLorentzVector corrP4 = scale * it->momentum;
205
206	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
207	thisjet->pt = corrP4.Et();
208	thisjet->eta = corrP4.Eta();
209	thisjet->phi = corrP4.Phi();
210	jet.push_back( *thisjet );
211	}// for jet
212	}// if, else
213	if( jet.size() < 2 && skim )
214	continue;
215	std::sort( jet.begin(), jet.end(), tree::EtGreater);
216
217	// met
218	std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
219	susy::MET* metobj = &(met_it->second);
220	met = metobj->met();
221
222	// electrons
223	std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
224	nElectron = eVector.size();
225
226	// vertices
227
228	nVertex = event->vertices.size();
229	weight = 1;
230
231	// bjets
232
233	tree->Fill();
234	} // for jentry
235
236
237
238	tree->Write();
239	outFile->cd();
240	outFile->Write();
241	outFile->Close();
242	}
243
244	int main(int argc, char** argv) {
245	TreeWriter *tw = new TreeWriter("../root-files/susyEvents_qcd_1000-inf_part.root", "myTree.root");
246	tw->SetProcessNEvents(-1);
247	tw->Loop();
248	}
249