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.12 by kiesel, Tue Apr 16 15:11:20 2013 UTC vs.
Revision 1.16 by kiesel, Thu Apr 18 16:26:44 2013 UTC

#	Line 1 | Line 1
1	–	#include<iostream>
2	–	#include<math.h>
3	–	#include<string>
4	–
5	–	#include "TSystem.h"
6	–
1		#include "treeWriter.h"
2
3		using namespace std;
4
5	<	TreeWriter::TreeWriter( TString inputName, TString outputName, int loggingVerbosity_ ) {
5	>	TreeWriter::TreeWriter( std::string inputName, std::string outputName, int loggingVerbosity_ ) {
6		// read the input file
7		inputTree = new TChain("susyTree");
8		if (loggingVerbosity_ > 0)
9		std::cout << "Add files to chain" << std::endl;
10	<	inputTree->Add( inputName );
10	>	inputTree->Add( inputName.c_str() );
11		Init( outputName, loggingVerbosity_ );
12		}
13
14	<	TreeWriter::TreeWriter( TChain* inputTree_, TString outputName, int loggingVerbosity_ ) {
14	>	TreeWriter::TreeWriter( TChain* inputTree_, std::string outputName, int loggingVerbosity_ ) {
15		inputTree = inputTree_;
16		Init( outputName, loggingVerbosity_ );
17		}
18
19	<
26	<	void TreeWriter::Init( TString outputName, int loggingVerbosity_ ) {
19	>	void TreeWriter::Init( std::string outputName, int loggingVerbosity_ ) {
20
21		if (loggingVerbosity_ > 0)
22		std::cout << "Set Branch Address of susy::Event" << std::endl;
23		event = new susy::Event;
24		inputTree->SetBranchAddress("susyEvent", &event);
25
26	+	// Here the number of proceeded events will be stored. For plotting, simply use L*sigma/eventNumber
27	+	eventNumbers = new TH1F("eventNumbers", "Histogram containing number of generated events", 1, 0, 1);
28	+	eventNumbers->GetXaxis()->SetBinLabel(1,"Number of generated events");
29	+
30		// open the output file
31		if (loggingVerbosity_>0)
32		std::cout << "Open file " << outputName << " for writing." << std::endl;
33	<	outFile = new TFile( outputName, "recreate" );
33	>	outFile = new TFile( outputName.c_str(), "recreate" );
34		tree = new TTree("susyTree","Tree for single photon analysis");
35
36		// set default parameter
#	Line 53 | Line 50 | TreeWriter::~TreeWriter() {
50		if (pileupHisto != 0 )
51		delete pileupHisto;
52		inputTree->GetCurrentFile()->Close();
53	+	delete inputTree;
54	+	delete event;
55	+	delete outFile;
56	+	delete tree;
57		}
58
59		// useful functions
60	<	float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
60	>	float deltaR( TLorentzVector v1, TLorentzVector v2 ) {
61		return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
62		}
63
#	Line 130 | Line 131 | float photonIso_corrected(susy::Photon g
131		return iso;
132		}
133
134	<	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
134	>	float d0correction( const susy::Electron& electron, const susy::Event& event ) {
135	>	TVector3 beamspot = event.vertices[0].position;
136	>	susy::Track track = event.tracks[electron.gsfTrackIndex];
137	>	float d0 = track.d0() - beamspot.X()*sin(track.phi()) + beamspot.Y()*cos(track.phi());
138	>	return d0;
139	>	}
140	>
141	>	float dZcorrection( const susy::Electron& electron, const susy::Event& event ) {
142	>	TVector3 beamspot = event.vertices[0].position;
143	>	susy::Track track = event.tracks[electron.gsfTrackIndex];
144	>
145	>	if(track.momentum.Pt() == 0.) return 1.e6;
146	>	float dz = (track.vertex.Z() - beamspot.Z()) - ((track.vertex.X() - beamspot.X())*track.momentum.Px() + (track.vertex.Y() - beamspot.Y())*track.momentum.Py()) / track.momentum.Pt() * (track.momentum.Pz() / track.momentum.Pt());
147	>	return dz;
148	>	}
149	>
150	>	float getPtFromMatchedJet( susy::Photon myPhoton, susy::PFJetCollection jetColl, int loggingVerbosity = 0 ) {
151		/**
152		* \brief Takes jet p_T as photon p_T
153		*
#	Line 142 | Line 159 | float TreeWriter::getPtFromMatchedJet( s
159		std::vector<susy::PFJet> nearJets;
160		nearJets.clear();
161
162	<	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
163	<	if(pfJets_it == myEvent.pfJets.end()){
164	<	if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
165	<	} else {
166	<	susy::PFJetCollection& jetColl = pfJets_it->second;
167	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
168	<	it != jetColl.end(); ++it) {
169	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
170	<	if (s_it == it->jecScaleFactors.end()) {
171	<	std::cout << "JEC is not available for this jet!!!" << std::endl;
172	<	continue;
173	<	}
174	<	float scale = s_it->second;
175	<	TLorentzVector corrP4 = scale * it->momentum;
176	<	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
160	<	if (deltaR_ > 0.3) continue;
161	<	if( loggingVerbosity > 0 )
162	<	std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
163	<	nearJets.push_back( *it );
164	<	}// for jet
165	<	}// if, else
162	>	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
163	>	it != jetColl.end(); ++it) {
164	>	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
165	>	if (s_it == it->jecScaleFactors.end()) {
166	>	std::cout << "JEC is not available for this jet!!!" << std::endl;
167	>	continue;
168	>	}
169	>	float scale = s_it->second;
170	>	TLorentzVector corrP4 = scale * it->momentum;
171	>	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
172	>	if (deltaR_ > 0.3) continue;
173	>	if( loggingVerbosity > 2 )
174	>	std::cout << " pT_jet / pT_gamma = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
175	>	nearJets.push_back( *it );
176	>	}// for jet
177
178		if ( nearJets.size() == 0 ) {
179		if( loggingVerbosity > 1 )
#	Line 196 | Line 207 | void TreeWriter::Loop() {
207
208		// get number of events to be proceeded
209		Long64_t nentries = inputTree->GetEntries();
210	+	// store them in histo
211	+	eventNumbers->Fill( "Number of generated events", nentries );
212		if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;
213
214		if( loggingVerbosity > 0 )
215		std::cout << "Processing " << processNEvents << " ouf of "
216		<< nentries << " events. " << std::endl;
217
205	–	tree::Photon *thisphoton = new tree::Photon();
206	–	tree::Jet *thisjet = new tree::Jet();
207	–
218		tree->Branch("photon", &photon);
219		tree->Branch("jet", &jet);
220		tree->Branch("electron", &electron);
#	Line 219 | Line 229 | void TreeWriter::Loop() {
229
230
231		for (long jentry=0; jentry < processNEvents; ++jentry) {
232	<	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
232	>	if ( loggingVerbosity>2 || jentry%reportEvery==0 )
233		std::cout << jentry << " / " << processNEvents << std :: endl;
234	<	inputTree->GetEntry(jentry);
234	>	inputTree->LoadTree( jentry );
235	>	inputTree->GetEntry( jentry );
236
237		photon.clear();
238		jet.clear();
#	Line 242 | Line 253 | void TreeWriter::Loop() {
253		pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
254		}
255
256	+	// get ak5 jets
257	+	susy::PFJetCollection jetVector;
258	+	if(event->pfJets.count("ak5") == 0)
259	+	cout << "ERROR: Jet collection 'ak5' not found" << endl;
260	+	else
261	+	jetVector = event->pfJets.find("ak5")->second;
262
263		// photons
264		if( loggingVerbosity > 1 )
265		std::cout << "Process photons" << std::endl;
266	<	std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
267	<	for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
268	<	it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
266	>	susy::PhotonCollection photonVector;
267	>	if(event->photons.count("photons") == 0)
268	>	cout << "ERROR: Photon collection 'photons' not found" << endl;
269	>	else
270	>	photonVector = event->photons.find("photons")->second;
271	>
272	>	for(susy::PhotonCollection :: iterator it = photonVector.begin();
273	>	it != photonVector.end(); ++it ) {
274		if( !(it->isEB() || it->isEE()) && skim )
275		continue;
276	<	thisphoton->pt = getPtFromMatchedJet( *it, *event );
276	>	tree::Photon thisphoton;
277	>	thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
278
279	<	thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
280	<	thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
281	<	thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
279	>	thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
280	>	thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
281	>	thisphoton.photonIso = photonIso_corrected(*it, event->rho25);
282
283	<	bool loose_photon_barrel = thisphoton->pt>20
283	>	bool loose_photon_barrel = thisphoton.pt>20
284		&& it->isEB()
285		&& it->passelectronveto
286		&& it->hadTowOverEm<0.05
287		&& it->sigmaIetaIeta<0.012
288	<	&& thisphoton->chargedIso<2.6
289	<	&& thisphoton->neutralIso<3.5+0.04*thisphoton->pt
290	<	&& thisphoton->photonIso<1.3+0.005*thisphoton->pt;
291	<	bool loose_photon_endcap = thisphoton->pt > 20
288	>	&& thisphoton.chargedIso<2.6
289	>	&& thisphoton.neutralIso<3.5+0.04*thisphoton.pt
290	>	&& thisphoton.photonIso<1.3+0.005*thisphoton.pt;
291	>	bool loose_photon_endcap = thisphoton.pt > 20
292		&& it->isEE()
293		&& it->passelectronveto
294		&& it->hadTowOverEm<0.05
295		&& it->sigmaIetaIeta<0.034
296	<	&& thisphoton->chargedIso<2.3
297	<	&& thisphoton->neutralIso<2.9+0.04*thisphoton->pt;
298	<	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton->pt > 75 ) && skim )
296	>	&& thisphoton.chargedIso<2.3
297	>	&& thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
298	>	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim )
299		continue;
300	<	thisphoton->eta = it->momentum.Eta();
301	<	thisphoton->phi = it->momentum.Phi();
302	<	thisphoton->r9 = it->r9;
303	<	thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
304	<	thisphoton->hadTowOverEm = it->hadTowOverEm;
305	<	thisphoton->pixelseed = it->nPixelSeeds;
306	<	thisphoton->conversionSafeVeto = it->passelectronveto;
307	<	photon.push_back( *thisphoton );
300	>	thisphoton.eta = it->momentum.Eta();
301	>	thisphoton.phi = it->momentum.Phi();
302	>	thisphoton.r9 = it->r9;
303	>	thisphoton.sigmaIetaIeta = it->sigmaIetaIeta;
304	>	thisphoton.hadTowOverEm = it->hadTowOverEm;
305	>	thisphoton.pixelseed = it->nPixelSeeds;
306	>	thisphoton.conversionSafeVeto = it->passelectronveto;
307	>	photon.push_back( thisphoton );
308		if( loggingVerbosity > 2 )
309	<	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
309	>	std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
310		}
311
312		if( photon.size() == 0 && skim )
#	Line 302 | Line 325 | void TreeWriter::Loop() {
325
326		for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
327		it != jetColl.end(); ++it) {
328	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
328	>	tree::Jet thisjet;
329	>	/*std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
330		if (s_it == it->jecScaleFactors.end()) {
331		std::cout << "JEC is not available for this jet!!!" << std::endl;
332		continue;
333		}
334		float scale = s_it->second;
335	+	*/
336	+	float scale = 1.;
337	+	if(it->jecScaleFactors.count("L2L3") == 0)
338	+	std::cout << "ERROR: JEC is not available for this jet" << std::endl;
339	+	else
340	+	scale = it->jecScaleFactors.find("L2L3")->second;
341		TLorentzVector corrP4 = scale * it->momentum;
342
343		if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
344	<	if(corrP4.Et() < 30 && skim ) continue;
345	<	thisjet->pt = corrP4.Et();
346	<	thisjet->eta = corrP4.Eta();
347	<	thisjet->phi = corrP4.Phi();
348	<	thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
344	>	if(corrP4.Pt() < 30 && skim ) continue;
345	>	thisjet.pt = corrP4.Pt();
346	>	thisjet.eta = corrP4.Eta();
347	>	thisjet.phi = corrP4.Phi();
348	>	thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
349		// jet composition
350	<	thisjet->chargedHadronEnergy = it->chargedHadronEnergy;
351	<	thisjet->neutralHadronEnergy = it->neutralHadronEnergy;
352	<	thisjet->photonEnergy = it->photonEnergy;
353	<	thisjet->electronEnergy = it->electronEnergy;
354	<	thisjet->muonEnergy = it->muonEnergy;
355	<	thisjet->HFHadronEnergy = it->HFHadronEnergy;
356	<	thisjet->HFEMEnergy = it->HFEMEnergy;
357	<	thisjet->chargedEmEnergy = it->chargedEmEnergy;
358	<	thisjet->chargedMuEnergy = it->chargedMuEnergy;
359	<	thisjet->neutralEmEnergy = it->neutralEmEnergy;
350	>	thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
351	>	thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
352	>	thisjet.photonEnergy = it->photonEnergy;
353	>	thisjet.electronEnergy = it->electronEnergy;
354	>	thisjet.muonEnergy = it->muonEnergy;
355	>	thisjet.HFHadronEnergy = it->HFHadronEnergy;
356	>	thisjet.HFEMEnergy = it->HFEMEnergy;
357	>	thisjet.chargedEmEnergy = it->chargedEmEnergy;
358	>	thisjet.chargedMuEnergy = it->chargedMuEnergy;
359	>	thisjet.neutralEmEnergy = it->neutralEmEnergy;
360
361		if( loggingVerbosity > 2 )
362	<	std::cout << " p_T, jet = " << thisjet->pt << std::endl;
362	>	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
363
364	<	jet.push_back( *thisjet );
365	<	ht += thisjet->pt;
364	>	jet.push_back( thisjet );
365	>	ht += thisjet.pt;
366		}// for jet
367		}// if, else
368		if( jet.size() < 2 && skim )
#	Line 358 | Line 388 | void TreeWriter::Loop() {
388		std::cout << " type1met = " << type1met << std::endl;
389
390		// electrons
391	<	tree::Particle* thiselectron = new tree::Particle();
392	<	map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
393	<	if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
394	<	cout << "gsfElectrons not found!" << endl;
395	<	} else {
396	<	for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
367	<	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
368	<	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
369	<	-effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
370	<	) / it->momentum.Pt();
371	<	if ( it->isEE() ){
372	<	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
373	<	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.15
374	<	|| it->sigmaIetaIeta > 0.01
375	<	|| it->hcalOverEcalBc > 0.12
376	<	|| it->vertex.Perp() > 0.02
377	<	|| it->vertex.Z() > 0.2
378	<	|| fabs(1./(it->ecalEnergy) - 1./(it->trackMomentums["AtVtx"].P())) > 0.05
379	<	|| it->convFlags() // not official, but perhaps substitude?
380	<	|| iso > 0.15 )
381	<	continue;
382	<	}
383	<	else if( it->isEB() ) {
384	<	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.009
385	<	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.10
386	<	|| it->sigmaIetaIeta > 0.03
387	<	|| it->hcalOverEcalBc > 0.10
388	<	|| it->vertex.Perp() > 0.02
389	<	|| it->vertex.Z() > 0.2
390	<	|| fabs(1./(it->ecalEnergy) - 1./(it->trackMomentums["AtVtx"].P())) > 0.05
391	<	|| it->convFlags() // not official, but perhaps substitude?
392	<	|| iso > 0.15 )
393	<	continue;
394	<	}
395	<	else // not in barrel nor in endcap
396	<	continue;
397	<	// TODO: conversion rejection information not implemented yet, see twiki for more details
391	>	std::vector<susy::Electron>* eVector;
392	>	if( event->electrons.count("gsfElectrons") == 0) {
393	>	cout << "EROR: no electron collection found" << endl;
394	>	continue;
395	>	} else
396	>	eVector = &event->electrons.find("gsfElectrons")->second;
397
398	<	thiselectron->pt = it->momentum.Et();
399	<	if( thiselectron->pt < 20 )
398	>	for(vector<susy::Electron>::iterator it = eVector->begin(); it < eVector->end(); ++it) {
399	>	tree::Particle thiselectron;
400	>	if( loggingVerbosity > 2 )
401	>	cout << " electron pt = " << it->momentum.Pt() << endl;
402	>	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
403	>	// use veto electrons
404	>	if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
405	>	continue; // spike rejection
406	>	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
407	>	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
408	>	) / it->momentum.Pt();
409	>	float d0 = d0correction( *it, *event );
410	>	float dZ = std::abs( dZcorrection( *it, *event ) );
411	>	if ( it->isEB() ){
412	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
413	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
414	>	|| it->sigmaIetaIeta > 0.01
415	>	|| it->hcalOverEcalBc > 0.15
416	>	|| d0 > 0.04
417	>	|| dZ > 0.2
418	>	|| iso > 0.15 )
419		continue;
420	<	if( loggingVerbosity > 2 )
421	<	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
422	<	thiselectron->eta = it->momentum.Eta();
423	<	thiselectron->phi = it->momentum.Phi();
424	<	electron.push_back( *thiselectron );
425	<	}
420	>	}
421	>	else if( it->isEE() ) {
422	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
423	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
424	>	|| it->sigmaIetaIeta > 0.03
425	>	|| d0 > 0.04
426	>	|| dZ > 0.2
427	>	|| iso > 0.15 )
428	>	continue;
429	>	}
430	>	else // not in barrel nor in endcap
431	>	continue;
432	>
433	>	thiselectron.pt = it->momentum.Pt();
434	>	if( loggingVerbosity > 2 )
435	>	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
436	>	thiselectron.eta = it->momentum.Eta();
437	>	thiselectron.phi = it->momentum.Phi();
438	>	electron.push_back( thiselectron );
439		}
440		if( loggingVerbosity > 1 )
441		std::cout << "Found " << electron.size() << " electrons" << std::endl;
442
443		// muons
444		std::vector<susy::Muon> mVector = event->muons["muons"];
445	<	tree::Particle* thismuon = new tree::Particle();
445	>	tree::Particle thismuon;
446		for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
447		if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
448		continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
449	<	thismuon->pt = it->momentum.Et();
450	<	if( thismuon->pt < 20 )
449	>	thismuon.pt = it->momentum.Et();
450	>	if( thismuon.pt < 20 )
451		continue;
452	<	thismuon->eta = it->momentum.Eta();
453	<	thismuon->phi = it->momentum.Phi();
454	<	muon.push_back( *thismuon );
452	>	thismuon.eta = it->momentum.Eta();
453	>	thismuon.phi = it->momentum.Phi();
454	>	muon.push_back( thismuon );
455		}
456		if( loggingVerbosity > 1 )
457		std::cout << "Found " << muon.size() << " muons" << std::endl;
#	Line 432 | Line 463 | void TreeWriter::Loop() {
463		if( ht < 450 && skim)
464		continue;
465
466	+
467		tree->Fill();
468		} // for jentry
469
470
471		outFile->cd();
472	+	eventNumbers->Write();
473		tree->Write();
474		outFile->Write();
475		outFile->Close();

Diff Legend

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