ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/kiesel/TreeWriter/treeWriter.cc

Comparing UserCode/kiesel/TreeWriter/treeWriter.cc (file contents):
Revision 1.4 by kiesel, Mon Mar 25 12:00:06 2013 UTC vs.
Revision 1.9 by kiesel, Sat Apr 13 09:53:02 2013 UTC

#	Line 1 | Line 1
1		#include<iostream>
2		#include<math.h>
3	+	#include<string>
4
5	<	#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	<	};
5	>	#include "TSystem.h"
6
7	<	TreeWriter::TreeWriter(TString inputName, TString outputName) {
7	>	#include "treeWriter.h"
8	>
9	>	using namespace std;
10	>
11	>	TreeWriter::TreeWriter( TString inputName, TString outputName, int loggingVerbosity_ ) {
12		// read the input file
13	<	inputFile = new TFile( inputName, "read" );
14	<	inputTree = (TTree*) inputFile->Get("susyTree");
13	>	inputTree = new TChain("susyTree");
14	>	if (loggingVerbosity_ > 0)
15	>	std::cout << "Add files to chain" << std::endl;
16	>	inputTree->Add( inputName );
17	>
18	>	if (loggingVerbosity_ > 0)
19	>	std::cout << "Set Branch Address of susy::Event" << std::endl;
20		event = new susy::Event;
21		inputTree->SetBranchAddress("susyEvent", &event);
22
23		// open the output file
24	+	if (loggingVerbosity_>0)
25	+	std::cout << "Open file " << outputName << " for writing." << std::endl;
26		outFile = new TFile( outputName, "recreate" );
27		tree = new TTree("susyTree","Tree for single photon analysis");
28
29		// set default parameter
30		processNEvents = -1;
31		reportEvery = 1000;
32	<	loggingVerbosity = 0;
32	>	loggingVerbosity = loggingVerbosity_;
33	>	skim = true;
34	>	pileupHisto = 0;
35	>	}
36
37	+	void TreeWriter::PileUpWeightFile( string pileupFileName ) {
38	+	TFile *puFile = new TFile( pileupFileName.c_str() );
39	+	pileupHisto = (TH1F*) puFile->Get("pileup");
40		}
41
42	+
43	+
44		TreeWriter::~TreeWriter() {
45	<	if (!inputTree) return;
46	<	delete inputTree->GetCurrentFile();
45	>	if (pileupHisto != 0 )
46	>	delete pileupHisto;
47	>	inputTree->GetCurrentFile()->Close();
48		}
49
50	+	// useful functions
51		float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
52		return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
53		}
54
55	<	float TreeWriter::getPtFromMatchedJet( susy::Photon photon, susy::Event event ) {
55	>	// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
56	>	float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
57	>	float eta = fabs(gamma.caloPosition.Eta());
58	>	float ea;
59	>
60	>	if(eta < 1.0) ea = 0.012;
61	>	else if(eta < 1.479) ea = 0.010;
62	>	else if(eta < 2.0) ea = 0.014;
63	>	else if(eta < 2.2) ea = 0.012;
64	>	else if(eta < 2.3) ea = 0.016;
65	>	else if(eta < 2.4) ea = 0.020;
66	>	else ea = 0.012;
67	>
68	>	float iso = gamma.chargedHadronIso;
69	>	iso = max(iso - rho*ea, (float)0.);
70	>
71	>	return iso;
72	>	}
73	>
74	>	float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
75	>	float eta = fabs(gamma.caloPosition.Eta());
76	>	float ea;
77	>
78	>	if(eta < 1.0) ea = 0.030;
79	>	else if(eta < 1.479) ea = 0.057;
80	>	else if(eta < 2.0) ea = 0.039;
81	>	else if(eta < 2.2) ea = 0.015;
82	>	else if(eta < 2.3) ea = 0.024;
83	>	else if(eta < 2.4) ea = 0.039;
84	>	else ea = 0.072;
85	>
86	>	float iso = gamma.neutralHadronIso;
87	>	iso = max(iso - rho*ea, (float)0.);
88	>
89	>	return iso;
90	>	}
91	>
92	>	float photonIso_corrected(susy::Photon gamma, float rho) {
93	>	float eta = fabs(gamma.caloPosition.Eta());
94	>	float ea;
95	>
96	>	if(eta < 1.0) ea = 0.148;
97	>	else if(eta < 1.479) ea = 0.130;
98	>	else if(eta < 2.0) ea = 0.112;
99	>	else if(eta < 2.2) ea = 0.216;
100	>	else if(eta < 2.3) ea = 0.262;
101	>	else if(eta < 2.4) ea = 0.260;
102	>	else ea = 0.266;
103	>
104	>	float iso = gamma.photonIso;
105	>	iso = max(iso - rho*ea, (float)0.);
106	>
107	>	return iso;
108	>	}
109	>
110	>	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
111		/**
112		* \brief Takes jet p_T as photon p_T
113		*
114		* At first all jets with DeltaR < 0.3 (isolation cone) are searched.
115		* If several jets are found, take the one with the minimal pt difference
116	<	* compared to the photon. If no such jets are found, return 0
116	>	* compared to the photon. If no such jets are found, keep the photon_pt
117		* TODO: remove photon matched jet from jet-selection?
118		*/
119		std::vector<susy::PFJet> nearJets;
120		nearJets.clear();
121
122	<	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event.pfJets.find("ak5");
123	<	if(pfJets_it == event.pfJets.end()){
124	<	if(event.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
122	>	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
123	>	if(pfJets_it == myEvent.pfJets.end()){
124	>	if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
125		} else {
126		susy::PFJetCollection& jetColl = pfJets_it->second;
127		for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
128	<	it != jetColl.end(); it++) {
128	>	it != jetColl.end(); ++it) {
129		std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
130		if (s_it == it->jecScaleFactors.end()) {
131		std::cout << "JEC is not available for this jet!!!" << std::endl;
#	Line 95 | Line 133 | float TreeWriter::getPtFromMatchedJet( s
133		}
134		float scale = s_it->second;
135		TLorentzVector corrP4 = scale * it->momentum;
136	<	float deltaR_ = deltaR(photon.momentum, corrP4 );
136	>	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
137		if (deltaR_ > 0.3) continue;
138	+	if( loggingVerbosity > 0 )
139	+	std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
140		nearJets.push_back( *it );
141		}// for jet
142		}// if, else
143
144		if ( nearJets.size() == 0 ) {
145	<	std::cout << "No jet with ΔR < .3 found, set p_T to 0" << std::endl;
146	<	return 0;
145	>	//std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
146	>	return myPhoton.momentum.Et();
147		}
148
149		float pt = 0;
150		float minPtDifferenz = 1E20; // should be very high
151		for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
152	<	it != jetEnd; it++ ) {
153	<	float ptDiff = fabs(photon.momentum.Et() - it->momentum.Et());
152	>	it != jetEnd; ++it ) {
153	>	float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
154		if (  ptDiff < minPtDifferenz ) {
155		minPtDifferenz = ptDiff;
156		pt = it->momentum.Et();
#	Line 118 | Line 158 | float TreeWriter::getPtFromMatchedJet( s
158		}
159
160		// testing
161	<	if( nearJets.size() > 1 )
161	>	if( nearJets.size() > 1 && loggingVerbosity > 0 )
162		std::cout << "There are several jets matching to this photon. "
163	<	<< "Please check if the upper assumtion is reasonable" << std::endl;
163	>	<< "Please check if jet-matching is correct." << std::endl;
164		return pt;
165		}
166
167
168		void TreeWriter::Loop() {
129	–	// only for testing
130	–	bool skim = true;
169
170		// here the event loop is implemented and the tree is filled
171		if (inputTree == 0) return;
#	Line 145 | Line 183 | void TreeWriter::Loop() {
183
184		tree->Branch("photon", &photon);
185		tree->Branch("jet", &jet);
186	+	tree->Branch("electron", &electron);
187	+	tree->Branch("muon", &muon);
188		tree->Branch("met", &met, "met/F");
189	+	tree->Branch("ht", &ht, "ht/F");
190		tree->Branch("nVertex", &nVertex, "nVertex/I");
191	<	tree->Branch("weigth", &weight, "weight/F");
151	<	tree->Branch("nElectron", &nElectron, "nElectron/I");
191	>	tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
192
193
194	<	for (long jentry=0; jentry < processNEvents; jentry++) {
194	>	for (long jentry=0; jentry < processNEvents; ++jentry) {
195		if ( loggingVerbosity>0 && jentry%reportEvery==0 )
196		std::cout << jentry << " / " << processNEvents << std :: endl;
197		inputTree->GetEntry(jentry);
198
199		photon.clear();
200		jet.clear();
201	+	electron.clear();
202	+	muon.clear();
203	+	ht = 0;
204	+
205	+	// weights
206	+	if (pileupHisto == 0) {
207	+	pu_weight = 1.;
208	+	} else {
209	+	float trueNumInteractions = -1;
210	+	for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
211	+	iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) {
212	+	if (iBX->BX == 0)
213	+	trueNumInteractions = iBX->trueNumInteractions;
214	+	}
215	+	pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
216	+	}
217	+
218
219		// photons
220	+	if( loggingVerbosity > 1 )
221	+	std::cout << "Process photons" << std::endl;
222		std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
223		for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
224	<	it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
224	>	it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
225		if( ! it->isEB() && skim )
226		continue;
227		thisphoton->pt = getPtFromMatchedJet( *it, *event );
#	Line 170 | Line 229 | void TreeWriter::Loop() {
229		continue;
230		thisphoton->eta = it->momentum.Eta();
231		thisphoton->phi = it->momentum.Phi();
232	<	thisphoton->chargedIso = it->chargedHadronIso;
233	<	thisphoton->neutralIso = it->neutralHadronIso;
234	<	thisphoton->photonIso = it->photonIso;
235	<	if ( it->r9 > 1 && skim ) // if == 1 ?
232	>	thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
233	>	thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
234	>	thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
235	>	if ( it->r9 >= 1 && skim )
236		continue;
237		thisphoton->r9 = it->r9;
238		thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
239		thisphoton->hadTowOverEm = it->hadTowOverEm;
240		thisphoton->pixelseed = it->nPixelSeeds;
241		photon.push_back( *thisphoton );
242	+	if( loggingVerbosity > 2 )
243	+	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
244		}
245	+
246		if( photon.size() == 0 && skim )
247		continue;
248		std::sort( photon.begin(), photon.end(), tree::EtGreater);
249	+	if( loggingVerbosity > 1 )
250	+	std::cout << "Found " << photon.size() << " photons" << std::endl;
251
252		// jets
253		std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
#	Line 194 | Line 258 | void TreeWriter::Loop() {
258		susy::PFJetCollection& jetColl = pfJets_it->second;
259
260		for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
261	<	it != jetColl.end(); it++) {
261	>	it != jetColl.end(); ++it) {
262		std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
263		if (s_it == it->jecScaleFactors.end()) {
264		std::cout << "JEC is not available for this jet!!!" << std::endl;
#	Line 204 | Line 268 | void TreeWriter::Loop() {
268		TLorentzVector corrP4 = scale * it->momentum;
269
270		if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
271	+	if(corrP4.Et() < 30 && skim ) continue;
272		thisjet->pt = corrP4.Et();
273		thisjet->eta = corrP4.Eta();
274		thisjet->phi = corrP4.Phi();
275	+	thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
276	+	if( loggingVerbosity > 2 )
277	+	std::cout << " p_T, jet = " << thisjet->pt << std::endl;
278	+
279		jet.push_back( *thisjet );
280	+	ht += thisjet->pt;
281		}// for jet
282		}// if, else
283		if( jet.size() < 2 && skim )
284		continue;
285		std::sort( jet.begin(), jet.end(), tree::EtGreater);
286	+	if( loggingVerbosity > 1 )
287	+	std::cout << "Found " << jet.size() << " jets" << std::endl;
288	+
289
290		// met
291		std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
292		susy::MET* metobj = &(met_it->second);
293		met = metobj->met();
294	+	met_phi = metobj->mEt.Phi();
295	+	if( loggingVerbosity > 2 )
296	+	std::cout << " met = " << met << std::endl;
297	+
298	+	std::map<TString, susy::MET>::iterator type1met_it = event->metMap.find("pfType1CorrectedMet");
299	+	susy::MET* type1metobj = &(type1met_it->second);
300	+	type1met = type1metobj->met();
301	+	type1met_phi = type1metobj->mEt.Phi();
302	+	if( loggingVerbosity > 2 )
303	+	std::cout << " type1met = " << type1met << std::endl;
304
305		// electrons
306	<	std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
224	<	nElectron = eVector.size();
306	>	tree::Particle* thiselectron = new tree::Particle();
307
308	<	// vertices
308	>	map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
309	>	if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
310	>	cout << "gsfElectrons not found!" << endl;
311	>	} else {
312	>	for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
313	>	thiselectron->pt = it->momentum.Et();
314	>	if( thiselectron->pt < 20 )
315	>	continue;
316	>	if( loggingVerbosity > 2 )
317	>	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
318	>	thiselectron->eta = it->momentum.Eta();
319	>	thiselectron->phi = it->momentum.Phi();
320	>	electron.push_back( *thiselectron );
321	>	}
322	>	}
323	>	if( loggingVerbosity > 1 )
324	>	std::cout << "Found " << electron.size() << " electrons" << std::endl;
325	>
326	>	// this seems not to work yet, where is the bug?
327	>	/*
328	>	std::vector<susy::Electron> eVector = event->electrons["gsfElectronsx"];
329	>	for( std::vector<susy::Electron>::iterator it = eVector.begin(); it != eVector.end(); ++it) {
330	>	thiselectron->pt = it->momentum.Et();
331	>	if( thiselectron->pt < 20 )
332	>	continue;
333	>	if( loggingVerbosity > 2 )
334	>	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
335	>	thiselectron->eta = it->momentum.Eta();
336	>	thiselectron->phi = it->momentum.Phi();
337	>	electron.push_back( *thiselectron );
338	>	}
339	>	*/
340
341	+	// muons
342	+	std::vector<susy::Muon> mVector = event->muons["muons"];
343	+	tree::Particle* thismuon = new tree::Particle();
344	+	for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
345	+	thismuon->pt = it->momentum.Et();
346	+	if( thismuon->pt < 20 )
347	+	continue;
348	+	thismuon->eta = it->momentum.Eta();
349	+	thismuon->phi = it->momentum.Phi();
350	+	muon.push_back( *thismuon );
351	+	}
352	+	if( loggingVerbosity > 1 )
353	+	std::cout << "Found " << muon.size() << " muons" << std::endl;
354	+
355	+
356	+	// vertices
357		nVertex = event->vertices.size();
229	–	weight = 1;
358
359	<	// bjets
359	>	if( ht < 450 && skim)
360	>	continue;
361
362		tree->Fill();
363		} // for jentry
364
365
237	–
238	–	tree->Write();
366		outFile->cd();
367	+	tree->Write();
368		outFile->Write();
369		outFile->Close();
370		}
371
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	–

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines