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.15 by kiesel, Thu Apr 18 13:38:55 2013 UTC vs.
Revision 1.20 by kiesel, Wed Apr 24 07:16:49 2013 UTC

#	Line 16 | Line 16 | TreeWriter::TreeWriter( TChain* inputTre
16		Init( outputName, loggingVerbosity_ );
17		}
18
19	–
19		void TreeWriter::Init( std::string outputName, int loggingVerbosity_ ) {
20
21		if (loggingVerbosity_ > 0)
#	Line 57 | Line 56 | TreeWriter::~TreeWriter() {
56		delete tree;
57		}
58
59	+	template<class Type>
60	+	std::vector<Type>* getVectorFromMap( map<TString, vector<Type> > myMap, TString search ){
61	+	// if nothing is found, return a empty vector (which has to be deleted never the less)
62	+	std::vector<Type>* vec;
63	+	typename map<TString, vector<Type> >::iterator mapIt = myMap.find( search );
64	+	if( mapIt == myMap.end() ) {
65	+	cout << "ERROR: Collection \"" << search << "\" not found!" << endl;
66	+	vec = new std::vector<Type>;
67	+	} else
68	+	vec = &mapIt->second;
69	+	cout << vec->size() << endl;
70	+	return vec;
71	+	}
72	+
73		// useful functions
74	<	float deltaR( TLorentzVector v1, TLorentzVector v2 ) {
74	>	float deltaR( const TLorentzVector& v1, const TLorentzVector& v2 ) {
75	>	// deltaR  = sqrt ( deltaEta^2 + deltaPhi^2 )
76		return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
77		}
78
79		float effectiveAreaElectron( float eta ) {
80	+	// needed by calculating the isolation for electrons
81		// see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection
82		// only for Delta R = 0.3 on 2012 Data
83		eta = fabs( eta );
#	Line 78 | Line 93 | float effectiveAreaElectron( float eta )
93		}
94
95		// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
96	<	float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
96	>	float chargedHadronIso_corrected(const susy::Photon& gamma, float rho) {
97		float eta = fabs(gamma.caloPosition.Eta());
98		float ea;
99
#	Line 96 | Line 111 | float chargedHadronIso_corrected(susy::P
111		return iso;
112		}
113
114	<	float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
114	>	float neutralHadronIso_corrected(const susy::Photon& gamma, float rho) {
115		float eta = fabs(gamma.caloPosition.Eta());
116		float ea;
117
#	Line 114 | Line 129 | float neutralHadronIso_corrected(susy::P
129		return iso;
130		}
131
132	<	float photonIso_corrected(susy::Photon gamma, float rho) {
132	>	float photonIso_corrected(const susy::Photon& gamma, float rho) {
133		float eta = fabs(gamma.caloPosition.Eta());
134		float ea;
135
#	Line 133 | Line 148 | float photonIso_corrected(susy::Photon g
148		}
149
150		float d0correction( const susy::Electron& electron, const susy::Event& event ) {
151	+	// copied from Brian Francis
152		TVector3 beamspot = event.vertices[0].position;
153		susy::Track track = event.tracks[electron.gsfTrackIndex];
154		float d0 = track.d0() - beamspot.X()*sin(track.phi()) + beamspot.Y()*cos(track.phi());
#	Line 140 | Line 156 | float d0correction( const susy::Electron
156		}
157
158		float dZcorrection( const susy::Electron& electron, const susy::Event& event ) {
159	+	// copied from Brian Francis
160		TVector3 beamspot = event.vertices[0].position;
161		susy::Track track = event.tracks[electron.gsfTrackIndex];
162
#	Line 148 | Line 165 | float dZcorrection( const susy::Electron
165		return dz;
166		}
167
168	<	float getPtFromMatchedJet( susy::Photon myPhoton, susy::PFJetCollection jetColl, int loggingVerbosity = 0 ) {
168	>	float getPtFromMatchedJet( const susy::Photon& myPhoton, const susy::PFJetCollection& jetColl, int loggingVerbosity = 0 ) {
169		/**
170		* \brief Takes jet p_T as photon p_T
171		*
#	Line 160 | Line 177 | float getPtFromMatchedJet( susy::Photon
177		std::vector<susy::PFJet> nearJets;
178		nearJets.clear();
179
180	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
180	>	for(std::vector<susy::PFJet>::const_iterator it = jetColl.begin();
181		it != jetColl.end(); ++it) {
182	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
182	>	float scale = 1.;
183	>	std::map<TString,Float_t>::const_iterator s_it = it->jecScaleFactors.find("L2L3");
184		if (s_it == it->jecScaleFactors.end()) {
185		std::cout << "JEC is not available for this jet!!!" << std::endl;
186		continue;
187	+	} else {
188	+	scale = s_it->second;
189		}
170	–	float scale = s_it->second;
190		TLorentzVector corrP4 = scale * it->momentum;
191		float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
192		if (deltaR_ > 0.3) continue;
#	Line 202 | Line 221 | float getPtFromMatchedJet( susy::Photon
221
222
223		void TreeWriter::Loop() {
224	+	/**
225	+	* \brief Loops over input chain and fills tree
226	+	*
227	+	* This is the major function of treeWriter, which initialize the output, loops
228	+	* over all events and fill the tree. In the end, the tree is saved to the
229	+	* output File
230	+	*/
231
232		// here the event loop is implemented and the tree is filled
233		if (inputTree == 0) return;
#	Line 209 | Line 235 | void TreeWriter::Loop() {
235		// get number of events to be proceeded
236		Long64_t nentries = inputTree->GetEntries();
237		// store them in histo
238	<	eventNumbers->Fill( "Number of generated Events", nentries );
238	>	eventNumbers->Fill( "Number of generated events", nentries );
239		if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;
240
241		if( loggingVerbosity > 0 )
#	Line 229 | Line 255 | void TreeWriter::Loop() {
255		tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
256
257
258	<	for (long jentry=0; jentry < processNEvents; ++jentry) {
259	<	if ( loggingVerbosity>2 || jentry%reportEvery==0 )
258	>	for (unsigned long jentry=0; jentry < processNEvents; ++jentry) {
259	>	if ( loggingVerbosity>1 || jentry%reportEvery==0 )
260		std::cout << jentry << " / " << processNEvents << std :: endl;
261	+	inputTree->LoadTree( jentry );
262	+	inputTree->GetEntry( jentry );
263
264		photon.clear();
265		jet.clear();
#	Line 253 | Line 281 | void TreeWriter::Loop() {
281		}
282
283		// get ak5 jets
284	<	susy::PFJetCollection jetVector;
257	<	if(event->pfJets.count("ak5") == 0)
258	<	cout << "ERROR: Jet collection 'ak5' not found" << endl;
259	<	else
260	<	jetVector = event->pfJets.find("ak5")->second;
284	>	std::vector<susy::PFJet> jetVector = event->pfJets["ak5"];
285
286		// photons
287	<	if( loggingVerbosity > 1 )
264	<	std::cout << "Process photons" << std::endl;
265	<	susy::PhotonCollection photonVector;
266	<	if(event->photons.count("photons") == 0)
267	<	cout << "ERROR: Photon collection 'photons' not found" << endl;
268	<	else
269	<	photonVector = event->photons.find("photons")->second;
287	>	std::vector<susy::Photon> photonVector = event->photons["photons"];
288
289	<	for(susy::PhotonCollection :: iterator it = photonVector.begin();
289	>	for(std::vector<susy::Photon>::iterator it = photonVector.begin();
290		it != photonVector.end(); ++it ) {
291	<	if( !(it->isEB() || it->isEE()) && skim )
291	>	if( !it->isEB() && skim )
292		continue;
293		tree::Photon thisphoton;
294		thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
#	Line 287 | Line 305 | void TreeWriter::Loop() {
305		&& thisphoton.chargedIso<2.6
306		&& thisphoton.neutralIso<3.5+0.04*thisphoton.pt
307		&& thisphoton.photonIso<1.3+0.005*thisphoton.pt;
308	<	bool loose_photon_endcap = thisphoton.pt > 20
308	>	/*bool loose_photon_endcap = thisphoton.pt > 20
309		&& it->isEE()
310		&& it->passelectronveto
311		&& it->hadTowOverEm<0.05
312		&& it->sigmaIetaIeta<0.034
313		&& thisphoton.chargedIso<2.3
314		&& thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
315	<	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim )
315	>	*/
316	>	if(!(loose_photon_barrel || thisphoton.pt > 75 ) && skim )
317		continue;
318		thisphoton.eta = it->momentum.Eta();
319		thisphoton.phi = it->momentum.Phi();
#	Line 315 | Line 334 | void TreeWriter::Loop() {
334		std::cout << "Found " << photon.size() << " photons" << std::endl;
335
336		// jets
318	–	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
319	–	if(pfJets_it == event->pfJets.end()){
320	–	if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
321	–	} else {
337
323	–	susy::PFJetCollection& jetColl = pfJets_it->second;
338
339	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
340	<	it != jetColl.end(); ++it) {
341	<	tree::Jet thisjet;
342	<	/*std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
343	<	if (s_it == it->jecScaleFactors.end()) {
344	<	std::cout << "JEC is not available for this jet!!!" << std::endl;
345	<	continue;
346	<	}
347	<	float scale = s_it->second;
348	<	*/
349	<	float scale = 1.;
350	<	if(it->jecScaleFactors.count("L2L3") == 0)
351	<	std::cout << "ERROR: JEC is not available for this jet" << std::endl;
352	<	else
353	<	scale = it->jecScaleFactors.find("L2L3")->second;
354	<	TLorentzVector corrP4 = scale * it->momentum;
341	<
342	<	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
343	<	if(corrP4.Pt() < 30 && skim ) continue;
344	<	thisjet.pt = corrP4.Pt();
345	<	thisjet.eta = corrP4.Eta();
346	<	thisjet.phi = corrP4.Phi();
347	<	thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
348	<	// jet composition
349	<	thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
350	<	thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
351	<	thisjet.photonEnergy = it->photonEnergy;
352	<	thisjet.electronEnergy = it->electronEnergy;
353	<	thisjet.muonEnergy = it->muonEnergy;
354	<	thisjet.HFHadronEnergy = it->HFHadronEnergy;
355	<	thisjet.HFEMEnergy = it->HFEMEnergy;
356	<	thisjet.chargedEmEnergy = it->chargedEmEnergy;
357	<	thisjet.chargedMuEnergy = it->chargedMuEnergy;
358	<	thisjet.neutralEmEnergy = it->neutralEmEnergy;
339	>	for(std::vector<susy::PFJet>::iterator it = jetVector.begin();
340	>	it != jetVector.end(); ++it) {
341	>	tree::Jet thisjet;
342	>
343	>	// scale with JEC
344	>	float scale = 1.;
345	>	if(it->jecScaleFactors.count("L2L3") == 0)
346	>	std::cout << "ERROR: JEC is not available for this jet" << std::endl;
347	>	else
348	>	scale = it->jecScaleFactors.find("L2L3")->second;
349	>	TLorentzVector corrP4 = scale * it->momentum;
350	>
351	>	// Calculate HT.
352	>	// The definiton differs from the saved jet, since trigger is described better
353	>	if( std::abs( corrP4.Eta() ) < 3 && corrP4.Pt > 40 )
354	>	ht += thisjet.pt;
355
356	<	if( loggingVerbosity > 2 )
357	<	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
356	>	if(std::abs(corrP4.Eta()) > 2.6 && skim ) continue;
357	>	if(corrP4.Pt() < 30 && skim ) continue;
358	>	thisjet.pt = corrP4.Pt();
359	>	thisjet.eta = corrP4.Eta();
360	>	thisjet.phi = corrP4.Phi();
361	>	thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
362	>	// jet composition
363	>	thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
364	>	thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
365	>	thisjet.photonEnergy = it->photonEnergy;
366	>	thisjet.electronEnergy = it->electronEnergy;
367	>	thisjet.muonEnergy = it->muonEnergy;
368	>	thisjet.HFHadronEnergy = it->HFHadronEnergy;
369	>	thisjet.HFEMEnergy = it->HFEMEnergy;
370	>	thisjet.chargedEmEnergy = it->chargedEmEnergy;
371	>	thisjet.chargedMuEnergy = it->chargedMuEnergy;
372	>	thisjet.neutralEmEnergy = it->neutralEmEnergy;
373
374	<	jet.push_back( thisjet );
375	<	ht += thisjet.pt;
376	<	}// for jet
377	<	}// if, else
374	>	if( loggingVerbosity > 2 )
375	>	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
376	>
377	>	jet.push_back( thisjet );
378	>	}// for jet
379		if( jet.size() < 2 && skim )
380		continue;
381		std::sort( jet.begin(), jet.end(), tree::EtGreater);
#	Line 387 | Line 399 | void TreeWriter::Loop() {
399		std::cout << " type1met = " << type1met << std::endl;
400
401		// electrons
402	<	std::vector<susy::Electron>* eVector;
403	<	if( event->electrons.count("gsfElectrons") == 0) {
392	<	cout << "EROR: no electron collection found" << endl;
393	<	continue;
394	<	} else
395	<	eVector = &event->electrons.find("gsfElectrons")->second;
396	<
397	<	for(vector<susy::Electron>::iterator it = eVector->begin(); it < eVector->end(); ++it) {
402	>	std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
403	>	for(std::vector<susy::Electron>::iterator it = eVector.begin(); it < eVector.end(); ++it) {
404		tree::Particle thiselectron;
405		if( loggingVerbosity > 2 )
406		cout << " electron pt = " << it->momentum.Pt() << endl;
#	Line 403 | Line 409 | void TreeWriter::Loop() {
409		if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
410		continue; // spike rejection
411		float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
412	<	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
412	>	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (float)0. )
413		) / it->momentum.Pt();
408	–	cout << iso << endl;
414		float d0 = d0correction( *it, *event );
410	–	cout << d0 << endl;
415		float dZ = std::abs( dZcorrection( *it, *event ) );
412	–	cout << dZ << endl;
416		if ( it->isEB() ){
414	–	cout << "electron in barrel" << endl;
417		if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
418		|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
419		|| it->sigmaIetaIeta > 0.01
420		|| it->hcalOverEcalBc > 0.15
421		|| d0 > 0.04
422		|| dZ > 0.2
423	<	|| iso > 0.15 ){
422	<	cout << " no real electron" << endl;
423	>	|| iso > 0.15 )
424		continue;
425		}
425	–	}
426		else if( it->isEE() ) {
427	–	cout << "electron in endcap" << endl;
427		if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
428		|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
429		|| it->sigmaIetaIeta > 0.03
430		|| d0 > 0.04
431		|| dZ > 0.2
432	<	|| iso > 0.15 ){
434	<	cout << "no real electron" << endl;
432	>	|| iso > 0.15 )
433		continue;
434		}
435	<	}
438	<	else{ // not in barrel nor in endcap
439	<	cout << " electron in gap" << endl;
435	>	else // not in barrel nor in endcap
436		continue;
437	<	// TODO: conversion rejection information not implemented yet, see twiki for more details
442	<	}
437	>
438		thiselectron.pt = it->momentum.Pt();
439		if( loggingVerbosity > 2 )
440		std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
441		thiselectron.eta = it->momentum.Eta();
442		thiselectron.phi = it->momentum.Phi();
448	–	cout << " adde electron" << endl;
443		electron.push_back( thiselectron );
444		}
445		if( loggingVerbosity > 1 )
446		std::cout << "Found " << electron.size() << " electrons" << std::endl;
447
448		// muons
455	–	std::vector<susy::Muon> mVector = event->muons["muons"];
449		tree::Particle thismuon;
450	+	std::vector<susy::Muon> mVector = event->muons["muons"];
451		for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
452		if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
453		continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
#	Line 467 | Line 461 | void TreeWriter::Loop() {
461		if( loggingVerbosity > 1 )
462		std::cout << "Found " << muon.size() << " muons" << std::endl;
463
470	–
464		// vertices
465		nVertex = event->vertices.size();
466

Diff Legend

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