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

#	Line 14 | Line 14 | TreeWriter::TreeWriter( TString inputNam
14		if (loggingVerbosity_ > 0)
15		std::cout << "Add files to chain" << std::endl;
16		inputTree->Add( inputName );
17	+	Init( outputName, loggingVerbosity_ );
18	+	}
19	+
20	+	TreeWriter::TreeWriter( TChain* inputTree_, TString outputName, int loggingVerbosity_ ) {
21	+	inputTree = inputTree_;
22	+	Init( outputName, loggingVerbosity_ );
23	+	}
24	+
25	+
26	+	void TreeWriter::Init( TString outputName, int loggingVerbosity_ ) {
27
28		if (loggingVerbosity_ > 0)
29		std::cout << "Set Branch Address of susy::Event" << std::endl;
#	Line 39 | Line 49 | void TreeWriter::PileUpWeightFile( strin
49		pileupHisto = (TH1F*) puFile->Get("pileup");
50		}
51
42	–
43	–
52		TreeWriter::~TreeWriter() {
53		if (pileupHisto != 0 )
54		delete pileupHisto;
#	Line 52 | Line 60 | float TreeWriter::deltaR( TLorentzVector
60		return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
61		}
62
63	+	float effectiveAreaElectron( float eta ) {
64	+	// see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection
65	+	// only for Delta R = 0.3 on 2012 Data
66	+	eta = fabs( eta );
67	+	float ea;
68	+	if( eta < 1.0 ) ea = 0.13;
69	+	else if( eta < 1.479 ) ea = 0.14;
70	+	else if( eta < 2.0 ) ea = 0.07;
71	+	else if( eta < 2.2 ) ea = 0.09;
72	+	else if( eta < 2.3 ) ea = 0.11;
73	+	else if( eta < 2.4 ) ea = 0.11;
74	+	else ea = 0.14;
75	+	return ea;
76	+	}
77	+
78		// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
79		float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
80	<	float eta = fabs(gamma.caloPosition.Eta());
81	<	float ea;
80	>	float eta = fabs(gamma.caloPosition.Eta());
81	>	float ea;
82
83	<	if(eta < 1.0) ea = 0.012;
84	<	else if(eta < 1.479) ea = 0.010;
85	<	else if(eta < 2.0) ea = 0.014;
86	<	else if(eta < 2.2) ea = 0.012;
87	<	else if(eta < 2.3) ea = 0.016;
88	<	else if(eta < 2.4) ea = 0.020;
89	<	else ea = 0.012;
83	>	if(eta < 1.0) ea = 0.012;
84	>	else if(eta < 1.479) ea = 0.010;
85	>	else if(eta < 2.0) ea = 0.014;
86	>	else if(eta < 2.2) ea = 0.012;
87	>	else if(eta < 2.3) ea = 0.016;
88	>	else if(eta < 2.4) ea = 0.020;
89	>	else ea = 0.012;
90
91	<	float iso = gamma.chargedHadronIso;
92	<	iso = max(iso - rho*ea, (float)0.);
91	>	float iso = gamma.chargedHadronIso;
92	>	iso = max(iso - rho*ea, (float)0.);
93
94	<	return iso;
94	>	return iso;
95		}
96
97		float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
98	<	float eta = fabs(gamma.caloPosition.Eta());
99	<	float ea;
98	>	float eta = fabs(gamma.caloPosition.Eta());
99	>	float ea;
100
101	<	if(eta < 1.0) ea = 0.030;
102	<	else if(eta < 1.479) ea = 0.057;
103	<	else if(eta < 2.0) ea = 0.039;
104	<	else if(eta < 2.2) ea = 0.015;
105	<	else if(eta < 2.3) ea = 0.024;
106	<	else if(eta < 2.4) ea = 0.039;
107	<	else ea = 0.072;
101	>	if(eta < 1.0) ea = 0.030;
102	>	else if(eta < 1.479) ea = 0.057;
103	>	else if(eta < 2.0) ea = 0.039;
104	>	else if(eta < 2.2) ea = 0.015;
105	>	else if(eta < 2.3) ea = 0.024;
106	>	else if(eta < 2.4) ea = 0.039;
107	>	else ea = 0.072;
108
109	<	float iso = gamma.neutralHadronIso;
110	<	iso = max(iso - rho*ea, (float)0.);
109	>	float iso = gamma.neutralHadronIso;
110	>	iso = max(iso - rho*ea, (float)0.);
111
112	<	return iso;
112	>	return iso;
113		}
114
115		float photonIso_corrected(susy::Photon gamma, float rho) {
116	<	float eta = fabs(gamma.caloPosition.Eta());
117	<	float ea;
116	>	float eta = fabs(gamma.caloPosition.Eta());
117	>	float ea;
118
119	<	if(eta < 1.0) ea = 0.148;
120	<	else if(eta < 1.479) ea = 0.130;
121	<	else if(eta < 2.0) ea = 0.112;
122	<	else if(eta < 2.2) ea = 0.216;
123	<	else if(eta < 2.3) ea = 0.262;
124	<	else if(eta < 2.4) ea = 0.260;
125	<	else ea = 0.266;
119	>	if(eta < 1.0) ea = 0.148;
120	>	else if(eta < 1.479) ea = 0.130;
121	>	else if(eta < 2.0) ea = 0.112;
122	>	else if(eta < 2.2) ea = 0.216;
123	>	else if(eta < 2.3) ea = 0.262;
124	>	else if(eta < 2.4) ea = 0.260;
125	>	else ea = 0.266;
126
127	<	float iso = gamma.photonIso;
128	<	iso = max(iso - rho*ea, (float)0.);
127	>	float iso = gamma.photonIso;
128	>	iso = max(iso - rho*ea, (float)0.);
129
130	<	return iso;
130	>	return iso;
131		}
132
133		float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
#	Line 142 | Line 165 | float TreeWriter::getPtFromMatchedJet( s
165		}// if, else
166
167		if ( nearJets.size() == 0 ) {
168	<	//std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
168	>	if( loggingVerbosity > 1 )
169	>	std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
170		return myPhoton.momentum.Et();
171		}
172
#	Line 186 | Line 210 | void TreeWriter::Loop() {
210		tree->Branch("electron", &electron);
211		tree->Branch("muon", &muon);
212		tree->Branch("met", &met, "met/F");
213	+	tree->Branch("metPhi", &met_phi, "metPhi/F");
214	+	tree->Branch("type1met", &type1met, "type1met/F");
215	+	tree->Branch("type1metPhi", &type1met_phi, "type1metPhi/F");
216		tree->Branch("ht", &ht, "ht/F");
217		tree->Branch("nVertex", &nVertex, "nVertex/I");
218		tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
#	Line 222 | Line 249 | void TreeWriter::Loop() {
249		std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
250		for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
251		it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
252	<	if( ! it->isEB() && skim )
252	>	if( !(it->isEB() || it->isEE()) && skim )
253		continue;
254		thisphoton->pt = getPtFromMatchedJet( *it, *event );
255	<	if( thisphoton->pt < 80 && skim )
229	<	continue;
230	<	thisphoton->eta = it->momentum.Eta();
231	<	thisphoton->phi = it->momentum.Phi();
255	>
256		thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
257		thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
258		thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
259	<	if ( it->r9 >= 1 && skim )
259	>
260	>	bool loose_photon_barrel = thisphoton->pt>20
261	>	&& it->isEB()
262	>	&& it->passelectronveto
263	>	&& it->hadTowOverEm<0.05
264	>	&& it->sigmaIetaIeta<0.012
265	>	&& thisphoton->chargedIso<2.6
266	>	&& thisphoton->neutralIso<3.5+0.04*thisphoton->pt
267	>	&& thisphoton->photonIso<1.3+0.005*thisphoton->pt;
268	>	bool loose_photon_endcap = thisphoton->pt > 20
269	>	&& it->isEE()
270	>	&& it->passelectronveto
271	>	&& it->hadTowOverEm<0.05
272	>	&& it->sigmaIetaIeta<0.034
273	>	&& thisphoton->chargedIso<2.3
274	>	&& thisphoton->neutralIso<2.9+0.04*thisphoton->pt;
275	>	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton->pt > 75 ) && skim )
276		continue;
277	+	thisphoton->eta = it->momentum.Eta();
278	+	thisphoton->phi = it->momentum.Phi();
279		thisphoton->r9 = it->r9;
280		thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
281		thisphoton->hadTowOverEm = it->hadTowOverEm;
282		thisphoton->pixelseed = it->nPixelSeeds;
283	+	thisphoton->conversionSafeVeto = it->passelectronveto;
284		photon.push_back( *thisphoton );
285		if( loggingVerbosity > 2 )
286		std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
#	Line 273 | Line 316 | void TreeWriter::Loop() {
316		thisjet->eta = corrP4.Eta();
317		thisjet->phi = corrP4.Phi();
318		thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
319	+	// jet composition
320	+	thisjet->chargedHadronEnergy = it->chargedHadronEnergy;
321	+	thisjet->neutralHadronEnergy = it->neutralHadronEnergy;
322	+	thisjet->photonEnergy = it->photonEnergy;
323	+	thisjet->electronEnergy = it->electronEnergy;
324	+	thisjet->muonEnergy = it->muonEnergy;
325	+	thisjet->HFHadronEnergy = it->HFHadronEnergy;
326	+	thisjet->HFEMEnergy = it->HFEMEnergy;
327	+	thisjet->chargedEmEnergy = it->chargedEmEnergy;
328	+	thisjet->chargedMuEnergy = it->chargedMuEnergy;
329	+	thisjet->neutralEmEnergy = it->neutralEmEnergy;
330	+
331		if( loggingVerbosity > 2 )
332		std::cout << " p_T, jet = " << thisjet->pt << std::endl;
333
#	Line 304 | Line 359 | void TreeWriter::Loop() {
359
360		// electrons
361		tree::Particle* thiselectron = new tree::Particle();
307	–
362		map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
363		if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
364		cout << "gsfElectrons not found!" << endl;
365		} else {
366		for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
367	<	thiselectron->pt = it->momentum.Et();
367	>	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
368	>	if( loggingVerbosity > 0)
369	>	std::cout << "TODO: piAtVtx is not calculated correctly" << std::endl;
370	>	//float pAtVtx = it->trackMomentums.find("AtVtxWithConstraint")->second.P();
371	>	float pAtVtx = it->ecalEnergy;
372	>	if( it->momentum.Pt() < 1  || it->momentum.Pt() > 1e6 || pAtVtx == 0 )
373	>	continue;
374	>	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
375	>	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
376	>	) / it->momentum.Pt();
377	>	if ( it->isEE() ){
378	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
379	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.15
380	>	|| it->sigmaIetaIeta > 0.01
381	>	|| it->hcalOverEcalBc > 0.12
382	>	|| it->vertex.Perp() > 0.02
383	>	|| it->vertex.Z() > 0.2
384	>	|| fabs(1./it->ecalEnergy - 1./pAtVtx ) > 0.05
385	>	|| iso > 0.15 )
386	>	continue;
387	>	}
388	>	else if( it->isEB() ) {
389	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.009
390	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.10
391	>	|| it->sigmaIetaIeta > 0.03
392	>	|| it->hcalOverEcalBc > 0.10
393	>	|| it->vertex.Perp() > 0.02
394	>	|| it->vertex.Z() > 0.2
395	>	|| fabs(1./it->ecalEnergy - 1./pAtVtx) > 0.05
396	>	|| iso > 0.15 )
397	>	continue;
398	>	}
399	>	else // not in barrel nor in endcap
400	>	continue;
401	>	// TODO: conversion rejection information not implemented yet, see twiki for more details
402	>
403	>	thiselectron->pt = it->momentum.Pt();
404		if( thiselectron->pt < 20 )
405		continue;
406		if( loggingVerbosity > 2 )
#	Line 323 | Line 413 | void TreeWriter::Loop() {
413		if( loggingVerbosity > 1 )
414		std::cout << "Found " << electron.size() << " electrons" << std::endl;
415
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	–
416		// muons
417		std::vector<susy::Muon> mVector = event->muons["muons"];
418		tree::Particle* thismuon = new tree::Particle();
419		for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
420	+	if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
421	+	continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
422		thismuon->pt = it->momentum.Et();
423		if( thismuon->pt < 20 )
424		continue;

Diff Legend

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