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.7 by kiesel, Tue Apr 9 12:13:55 2013 UTC vs.
Revision 1.9 by kiesel, Sat Apr 13 09:53:02 2013 UTC

#	Line 8 | Line 8
8
9		using namespace std;
10
11	<	TreeWriter::TreeWriter(TString inputName, TString outputName, int loggingVerbosity_ ) {
11	>	TreeWriter::TreeWriter( TString inputName, TString outputName, int loggingVerbosity_ ) {
12		// read the input file
13		inputTree = new TChain("susyTree");
14		if (loggingVerbosity_ > 0)
#	Line 20 | Line 20 | TreeWriter::TreeWriter(TString inputName
20		event = new susy::Event;
21		inputTree->SetBranchAddress("susyEvent", &event);
22
23	<	if (outputName == "") {
24	<	// Here the output filename is generated automaticly
25	<	try{
26	<	string fName = (string) inputName;
27	<	// Search for 'susyEvents_', and get the unique characters afterwards.
28	<	// eg /path/susyEvents_123_Myl.root -> susyTree_123_Mly.root
29	<	outputName = "susyTree_" + fName.substr(fName.find("susyEvents_")+11,-1);
30	<	} catch (int e ) {
31	<	outputName = "susyTree.root";
32	<	}
33	<	}
34	<
23	>	// open the output file
24		if (loggingVerbosity_>0)
25		std::cout << "Open file " << outputName << " for writing." << std::endl;
37	–	// open the output file
26		outFile = new TFile( outputName, "recreate" );
27		tree = new TTree("susyTree","Tree for single photon analysis");
28
#	Line 43 | Line 31 | TreeWriter::TreeWriter(TString inputName
31		reportEvery = 1000;
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	+	// 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
#	Line 73 | Line 125 | float TreeWriter::getPtFromMatchedJet( s
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 97 | Line 149 | float TreeWriter::getPtFromMatchedJet( s
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++ ) {
152	>	it != jetEnd; ++it ) {
153		float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
154		if (  ptDiff < minPtDifferenz ) {
155		minPtDifferenz = ptDiff;
#	Line 131 | 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("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 159 | 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 );
172	–	ht += thisphoton->pt;
242		if( loggingVerbosity > 2 )
243		std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
244		}
#	Line 189 | 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 199 | 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
#	Line 209 | Line 283 | void TreeWriter::Loop() {
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"];
220	<	//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();
358
359		if( ht < 450 && skim)
#	Line 230 | Line 363 | void TreeWriter::Loop() {
363		} // for jentry
364
365
233	–
234	–	tree->Write();
366		outFile->cd();
367	+	tree->Write();
368		outFile->Write();
369		outFile->Close();
370		}

Diff Legend

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