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.16 by kiesel, Thu Apr 18 16:26:44 2013 UTC vs.
Revision 1.32 by kiesel, Tue May 7 19:41:48 2013 UTC

#	Line 56 | 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	+	float deltaPhi( float phi1, float phi2) {
74	+	float result = phi1 - phi2;
75	+	while (result > M_PI) result -= 2*M_PI;
76	+	while (result <= -M_PI) result += 2*M_PI;
77	+	return result;
78	+	}
79	+
80		// useful functions
81	<	float deltaR( TLorentzVector v1, TLorentzVector v2 ) {
82	<	return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
81	>	float deltaR( const TLorentzVector& v1, const TLorentzVector& v2 ) {
82	>	// deltaR  = sqrt ( deltaEta^2 + deltaPhi^2 )
83	>	return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(deltaPhi(v1.Phi(),v2.Phi()), 2) );
84		}
85
86		float effectiveAreaElectron( float eta ) {
87	+	// needed by calculating the isolation for electrons
88		// see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection
89		// only for Delta R = 0.3 on 2012 Data
90		eta = fabs( eta );
#	Line 77 | Line 100 | float effectiveAreaElectron( float eta )
100		}
101
102		// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
103	<	float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
103	>	float chargedHadronIso_corrected(const susy::Photon& gamma, float rho) {
104		float eta = fabs(gamma.caloPosition.Eta());
105		float ea;
106
#	Line 95 | Line 118 | float chargedHadronIso_corrected(susy::P
118		return iso;
119		}
120
121	<	float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
121	>	float neutralHadronIso_corrected(const susy::Photon& gamma, float rho) {
122		float eta = fabs(gamma.caloPosition.Eta());
123		float ea;
124
#	Line 113 | Line 136 | float neutralHadronIso_corrected(susy::P
136		return iso;
137		}
138
139	<	float photonIso_corrected(susy::Photon gamma, float rho) {
139	>	float photonIso_corrected(const susy::Photon& gamma, float rho) {
140		float eta = fabs(gamma.caloPosition.Eta());
141		float ea;
142
#	Line 132 | Line 155 | float photonIso_corrected(susy::Photon g
155		}
156
157		float d0correction( const susy::Electron& electron, const susy::Event& event ) {
158	+	// copied from Brian Francis
159		TVector3 beamspot = event.vertices[0].position;
160		susy::Track track = event.tracks[electron.gsfTrackIndex];
161		float d0 = track.d0() - beamspot.X()*sin(track.phi()) + beamspot.Y()*cos(track.phi());
#	Line 139 | Line 163 | float d0correction( const susy::Electron
163		}
164
165		float dZcorrection( const susy::Electron& electron, const susy::Event& event ) {
166	+	// copied from Brian Francis
167		TVector3 beamspot = event.vertices[0].position;
168		susy::Track track = event.tracks[electron.gsfTrackIndex];
169
#	Line 147 | Line 172 | float dZcorrection( const susy::Electron
172		return dz;
173		}
174
175	<	float getPtFromMatchedJet( susy::Photon myPhoton, susy::PFJetCollection jetColl, int loggingVerbosity = 0 ) {
175	>	float getPtFromMatchedJet( const susy::Photon& myPhoton, const susy::PFJetCollection& jetColl, int loggingVerbosity = 0 ) {
176		/**
177		* \brief Takes jet p_T as photon p_T
178		*
#	Line 159 | Line 184 | float getPtFromMatchedJet( susy::Photon
184		std::vector<susy::PFJet> nearJets;
185		nearJets.clear();
186
187	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
187	>	for(std::vector<susy::PFJet>::const_iterator it = jetColl.begin();
188		it != jetColl.end(); ++it) {
189	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
189	>	float scale = 1.;
190	>	std::map<TString,Float_t>::const_iterator s_it = it->jecScaleFactors.find("L2L3");
191		if (s_it == it->jecScaleFactors.end()) {
192		std::cout << "JEC is not available for this jet!!!" << std::endl;
193		continue;
194	+	} else {
195	+	scale = s_it->second;
196		}
169	–	float scale = s_it->second;
197		TLorentzVector corrP4 = scale * it->momentum;
198		float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
199		if (deltaR_ > 0.3) continue;
#	Line 201 | Line 228 | float getPtFromMatchedJet( susy::Photon
228
229
230		void TreeWriter::Loop() {
231	+	/**
232	+	* \brief Loops over input chain and fills tree
233	+	*
234	+	* This is the major function of treeWriter, which initialize the output, loops
235	+	* over all events and fill the tree. In the end, the tree is saved to the
236	+	* output File
237	+	*/
238
239		// here the event loop is implemented and the tree is filled
240		if (inputTree == 0) return;
#	Line 225 | Line 259 | void TreeWriter::Loop() {
259		tree->Branch("type1metPhi", &type1met_phi, "type1metPhi/F");
260		tree->Branch("ht", &ht, "ht/F");
261		tree->Branch("nVertex", &nVertex, "nVertex/I");
262	<	tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
263	<
262	>	tree->Branch("weight", &weight, "weight/D");
263	>	tree->Branch("genElectron", &genElectron);
264	>	tree->Branch("genPhoton", &genPhoton);
265
266		for (long jentry=0; jentry < processNEvents; ++jentry) {
267	<	if ( loggingVerbosity>2 || jentry%reportEvery==0 )
267	>	if ( loggingVerbosity>1 || jentry%reportEvery==0 )
268		std::cout << jentry << " / " << processNEvents << std :: endl;
269		inputTree->LoadTree( jentry );
270		inputTree->GetEntry( jentry );
#	Line 238 | Line 273 | void TreeWriter::Loop() {
273		jet.clear();
274		electron.clear();
275		muon.clear();
276	+	genElectron.clear();
277	+	genPhoton.clear();
278		ht = 0;
279
280		// weights
281		if (pileupHisto == 0) {
282	<	pu_weight = 1.;
282	>	weight = 1.;
283		} else {
284		float trueNumInteractions = -1;
285		for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
#	Line 250 | Line 287 | void TreeWriter::Loop() {
287		if (iBX->BX == 0)
288		trueNumInteractions = iBX->trueNumInteractions;
289		}
290	<	pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
290	>	weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
291		}
292
293		// get ak5 jets
294	<	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;
294	>	std::vector<susy::PFJet> jetVector = event->pfJets["ak5"];
295
296		// photons
297	<	if( loggingVerbosity > 1 )
265	<	std::cout << "Process photons" << std::endl;
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;
297	>	std::vector<susy::Photon> photonVector = event->photons["photons"];
298
299	<	for(susy::PhotonCollection :: iterator it = photonVector.begin();
299	>	for(std::vector<susy::Photon>::iterator it = photonVector.begin();
300		it != photonVector.end(); ++it ) {
301	<	if( !(it->isEB() || it->isEE()) && skim )
301	>	if( !(it->isEE() || it->isEB()) && it->momentum.Pt()<20 && it->isEBEtaGap() && it->isEBPhiGap() && it->isEERingGap() && it->isEEDeeGap() && it->isEBEEGap() && skim )
302		continue;
303		tree::Photon thisphoton;
277	–	thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
304
305		thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
306		thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
307		thisphoton.photonIso = photonIso_corrected(*it, event->rho25);
308
309	<	bool loose_photon_barrel = thisphoton.pt>20
284	<	&& it->isEB()
285	<	&& it->passelectronveto
309	>	bool loose_photon_barrel = it->isEB()
310		&& it->hadTowOverEm<0.05
311		&& it->sigmaIetaIeta<0.012
312		&& thisphoton.chargedIso<2.6
313		&& thisphoton.neutralIso<3.5+0.04*thisphoton.pt
314		&& thisphoton.photonIso<1.3+0.005*thisphoton.pt;
315	<	bool loose_photon_endcap = thisphoton.pt > 20
316	<	&& it->isEE()
293	<	&& it->passelectronveto
315	>
316	>	bool loose_photon_endcap = it->isEE()
317		&& it->hadTowOverEm<0.05
318		&& it->sigmaIetaIeta<0.034
319		&& thisphoton.chargedIso<2.3
320		&& thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
321	<	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim )
321	>
322	>	thisphoton.ptJet = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
323	>	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.ptJet > 75 ) && skim )
324		continue;
325	+	thisphoton.pt = it->momentum.Pt();
326		thisphoton.eta = it->momentum.Eta();
327		thisphoton.phi = it->momentum.Phi();
328		thisphoton.r9 = it->r9;
#	Line 304 | Line 330 | void TreeWriter::Loop() {
330		thisphoton.hadTowOverEm = it->hadTowOverEm;
331		thisphoton.pixelseed = it->nPixelSeeds;
332		thisphoton.conversionSafeVeto = it->passelectronveto;
333	+	thisphoton.genInformation = 0;
334		photon.push_back( thisphoton );
335		if( loggingVerbosity > 2 )
336		std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
#	Line 316 | Line 343 | void TreeWriter::Loop() {
343		std::cout << "Found " << photon.size() << " photons" << std::endl;
344
345		// jets
319	–	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
320	–	if(pfJets_it == event->pfJets.end()){
321	–	if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
322	–	} else {
346
324	–	susy::PFJetCollection& jetColl = pfJets_it->second;
347
348	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
349	<	it != jetColl.end(); ++it) {
350	<	tree::Jet thisjet;
351	<	/*std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
352	<	if (s_it == it->jecScaleFactors.end()) {
353	<	std::cout << "JEC is not available for this jet!!!" << std::endl;
354	<	continue;
355	<	}
356	<	float scale = s_it->second;
357	<	*/
358	<	float scale = 1.;
359	<	if(it->jecScaleFactors.count("L2L3") == 0)
360	<	std::cout << "ERROR: JEC is not available for this jet" << std::endl;
361	<	else
362	<	scale = it->jecScaleFactors.find("L2L3")->second;
363	<	TLorentzVector corrP4 = scale * it->momentum;
342	<
343	<	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
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;
348	>	for(std::vector<susy::PFJet>::iterator it = jetVector.begin();
349	>	it != jetVector.end(); ++it) {
350	>	tree::Jet thisjet;
351	>
352	>	// scale with JEC
353	>	float scale = 1.;
354	>	if(it->jecScaleFactors.count("L2L3") == 0)
355	>	std::cout << "ERROR: JEC is not available for this jet" << std::endl;
356	>	else
357	>	scale = it->jecScaleFactors.find("L2L3")->second;
358	>	TLorentzVector corrP4 = scale * it->momentum;
359	>
360	>	// Calculate HT.
361	>	// The definiton differs from the saved jet, since trigger is described better
362	>	if( std::abs( corrP4.Eta() ) < 3 && corrP4.Pt() > 40 )
363	>	ht += thisjet.pt;
364
365	<	if( loggingVerbosity > 2 )
366	<	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
365	>	if(std::abs(corrP4.Eta()) > 2.6 && skim ) continue;
366	>	if(corrP4.Pt() < 30 && skim ) continue;
367	>	thisjet.pt = corrP4.Pt();
368	>	thisjet.eta = corrP4.Eta();
369	>	thisjet.phi = corrP4.Phi();
370	>	thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
371	>	// jet composition
372	>	thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
373	>	thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
374	>	thisjet.photonEnergy = it->photonEnergy;
375	>	thisjet.electronEnergy = it->electronEnergy;
376	>	thisjet.muonEnergy = it->muonEnergy;
377	>	thisjet.HFHadronEnergy = it->HFHadronEnergy;
378	>	thisjet.HFEMEnergy = it->HFEMEnergy;
379	>	thisjet.chargedEmEnergy = it->chargedEmEnergy;
380	>	thisjet.chargedMuEnergy = it->chargedMuEnergy;
381	>	thisjet.neutralEmEnergy = it->neutralEmEnergy;
382
383	<	jet.push_back( thisjet );
384	<	ht += thisjet.pt;
385	<	}// for jet
386	<	}// if, else
383	>	if( loggingVerbosity > 2 )
384	>	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
385	>
386	>	jet.push_back( thisjet );
387	>	}// for jet
388		if( jet.size() < 2 && skim )
389		continue;
390		std::sort( jet.begin(), jet.end(), tree::EtGreater);
391		if( loggingVerbosity > 1 )
392		std::cout << "Found " << jet.size() << " jets" << std::endl;
393
394	+	if( ht < 450 && skim)
395	+	continue;
396	+
397	+
398
399		// met
400		std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
#	Line 388 | Line 412 | void TreeWriter::Loop() {
412		std::cout << " type1met = " << type1met << std::endl;
413
414		// electrons
415	<	std::vector<susy::Electron>* eVector;
416	<	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	<	for(vector<susy::Electron>::iterator it = eVector->begin(); it < eVector->end(); ++it) {
415	>	std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
416	>	for(std::vector<susy::Electron>::iterator it = eVector.begin(); it < eVector.end(); ++it) {
417		tree::Particle thiselectron;
418		if( loggingVerbosity > 2 )
419		cout << " electron pt = " << it->momentum.Pt() << endl;
#	Line 403 | Line 421 | void TreeWriter::Loop() {
421		// use veto electrons
422		if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
423		continue; // spike rejection
424	<	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
407	<	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
424	>	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (float)0. )
425		) / it->momentum.Pt();
426		float d0 = d0correction( *it, *event );
427		float dZ = std::abs( dZcorrection( *it, *event ) );
#	Line 441 | Line 458 | void TreeWriter::Loop() {
458		std::cout << "Found " << electron.size() << " electrons" << std::endl;
459
460		// muons
444	–	std::vector<susy::Muon> mVector = event->muons["muons"];
461		tree::Particle thismuon;
462	+	std::vector<susy::Muon> mVector = event->muons["muons"];
463		for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
464		if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
465		continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
#	Line 456 | Line 473 | void TreeWriter::Loop() {
473		if( loggingVerbosity > 1 )
474		std::cout << "Found " << muon.size() << " muons" << std::endl;
475
459	–
476		// vertices
477		nVertex = event->vertices.size();
478
479	<	if( ht < 450 && skim)
480	<	continue;
481	<
479	>	tree::Particle thisGenParticle;
480	>	for( std::vector<susy::Particle>::iterator it = event->genParticles.begin(); it != event->genParticles.end(); ++it ) {
481	>	if( it->momentum.Pt() < 20 ) continue;
482	>	thisGenParticle.pt = it->momentum.Pt();
483	>	thisGenParticle.eta = it->momentum.Eta();
484	>	thisGenParticle.phi = it->momentum.Phi();
485	>	switch( std::abs(it->pdgId) ) {
486	>	case 22: // photon
487	>	genPhoton.push_back( thisGenParticle );
488	>	break;
489	>	case 11: // electron
490	>	// Demand a W boson as mother particle of electron
491	>	if( abs(event->genParticles[it->motherIndex].pdgId) == 24 )
492	>	genElectron.push_back( thisGenParticle );
493	>	break;
494	>	}
495	>	}
496
497		tree->Fill();
498		} // for jentry

Diff Legend

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