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.8 by kiesel, Wed Apr 10 20:46:33 2013 UTC vs.
Revision 1.15 by kiesel, Thu Apr 18 13:38:55 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_, std::string outputName, int loggingVerbosity_ ) {
15	>	inputTree = inputTree_;
16	>	Init( outputName, loggingVerbosity_ );
17	>	}
18	>
19	>
20	>	void TreeWriter::Init( std::string outputName, int loggingVerbosity_ ) {
21
22		if (loggingVerbosity_ > 0)
23		std::cout << "Set Branch Address of susy::Event" << std::endl;
24		event = new susy::Event;
25		inputTree->SetBranchAddress("susyEvent", &event);
26
27	<	if (outputName == "") {
28	<	// Here the output filename is generated automaticly
29	<	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	<	}
27	>	// Here the number of proceeded events will be stored. For plotting, simply use L*sigma/eventNumber
28	>	eventNumbers = new TH1F("eventNumbers", "Histogram containing number of generated events", 1, 0, 1);
29	>	eventNumbers->GetXaxis()->SetBinLabel(1,"Number of generated events");
30
31	+	// open the output file
32		if (loggingVerbosity_>0)
33		std::cout << "Open file " << outputName << " for writing." << std::endl;
34	<	// open the output file
38	<	outFile = new TFile( outputName, "recreate" );
34	>	outFile = new TFile( outputName.c_str(), "recreate" );
35		tree = new TTree("susyTree","Tree for single photon analysis");
36
37		// set default parameter
#	Line 43 | Line 39 | TreeWriter::TreeWriter(TString inputName
39		reportEvery = 1000;
40		loggingVerbosity = loggingVerbosity_;
41		skim = true;
42	+	pileupHisto = 0;
43	+	}
44
45	+	void TreeWriter::PileUpWeightFile( string pileupFileName ) {
46	+	TFile *puFile = new TFile( pileupFileName.c_str() );
47	+	pileupHisto = (TH1F*) puFile->Get("pileup");
48		}
49
50		TreeWriter::~TreeWriter() {
51	<	if (!inputTree) return;
52	<	delete inputTree->GetCurrentFile();
51	>	if (pileupHisto != 0 )
52	>	delete pileupHisto;
53	>	inputTree->GetCurrentFile()->Close();
54	>	delete inputTree;
55	>	delete event;
56	>	delete outFile;
57	>	delete tree;
58		}
59
60		// useful functions
61	<	float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
61	>	float deltaR( TLorentzVector v1, TLorentzVector v2 ) {
62		return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
63		}
64
65	+	float effectiveAreaElectron( float eta ) {
66	+	// see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection
67	+	// only for Delta R = 0.3 on 2012 Data
68	+	eta = fabs( eta );
69	+	float ea;
70	+	if( eta < 1.0 ) ea = 0.13;
71	+	else if( eta < 1.479 ) ea = 0.14;
72	+	else if( eta < 2.0 ) ea = 0.07;
73	+	else if( eta < 2.2 ) ea = 0.09;
74	+	else if( eta < 2.3 ) ea = 0.11;
75	+	else if( eta < 2.4 ) ea = 0.11;
76	+	else ea = 0.14;
77	+	return ea;
78	+	}
79	+
80		// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
81		float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
82	<	float eta = fabs(gamma.caloPosition.Eta());
83	<	float ea;
82	>	float eta = fabs(gamma.caloPosition.Eta());
83	>	float ea;
84
85	<	if(eta < 1.0) ea = 0.012;
86	<	else if(eta < 1.479) ea = 0.010;
87	<	else if(eta < 2.0) ea = 0.014;
88	<	else if(eta < 2.2) ea = 0.012;
89	<	else if(eta < 2.3) ea = 0.016;
90	<	else if(eta < 2.4) ea = 0.020;
91	<	else ea = 0.012;
85	>	if(eta < 1.0) ea = 0.012;
86	>	else if(eta < 1.479) ea = 0.010;
87	>	else if(eta < 2.0) ea = 0.014;
88	>	else if(eta < 2.2) ea = 0.012;
89	>	else if(eta < 2.3) ea = 0.016;
90	>	else if(eta < 2.4) ea = 0.020;
91	>	else ea = 0.012;
92
93	<	float iso = gamma.chargedHadronIso;
94	<	iso = max(iso - rho*ea, (float)0.);
93	>	float iso = gamma.chargedHadronIso;
94	>	iso = max(iso - rho*ea, (float)0.);
95
96	<	return iso;
96	>	return iso;
97		}
98
99		float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
100	<	float eta = fabs(gamma.caloPosition.Eta());
101	<	float ea;
100	>	float eta = fabs(gamma.caloPosition.Eta());
101	>	float ea;
102
103	<	if(eta < 1.0) ea = 0.030;
104	<	else if(eta < 1.479) ea = 0.057;
105	<	else if(eta < 2.0) ea = 0.039;
106	<	else if(eta < 2.2) ea = 0.015;
107	<	else if(eta < 2.3) ea = 0.024;
108	<	else if(eta < 2.4) ea = 0.039;
109	<	else ea = 0.072;
103	>	if(eta < 1.0) ea = 0.030;
104	>	else if(eta < 1.479) ea = 0.057;
105	>	else if(eta < 2.0) ea = 0.039;
106	>	else if(eta < 2.2) ea = 0.015;
107	>	else if(eta < 2.3) ea = 0.024;
108	>	else if(eta < 2.4) ea = 0.039;
109	>	else ea = 0.072;
110
111	<	float iso = gamma.neutralHadronIso;
112	<	iso = max(iso - rho*ea, (float)0.);
111	>	float iso = gamma.neutralHadronIso;
112	>	iso = max(iso - rho*ea, (float)0.);
113
114	<	return iso;
114	>	return iso;
115		}
116
117		float photonIso_corrected(susy::Photon gamma, float rho) {
118	<	float eta = fabs(gamma.caloPosition.Eta());
119	<	float ea;
118	>	float eta = fabs(gamma.caloPosition.Eta());
119	>	float ea;
120
121	<	if(eta < 1.0) ea = 0.148;
122	<	else if(eta < 1.479) ea = 0.130;
123	<	else if(eta < 2.0) ea = 0.112;
124	<	else if(eta < 2.2) ea = 0.216;
125	<	else if(eta < 2.3) ea = 0.262;
126	<	else if(eta < 2.4) ea = 0.260;
127	<	else ea = 0.266;
121	>	if(eta < 1.0) ea = 0.148;
122	>	else if(eta < 1.479) ea = 0.130;
123	>	else if(eta < 2.0) ea = 0.112;
124	>	else if(eta < 2.2) ea = 0.216;
125	>	else if(eta < 2.3) ea = 0.262;
126	>	else if(eta < 2.4) ea = 0.260;
127	>	else ea = 0.266;
128
129	<	float iso = gamma.photonIso;
130	<	iso = max(iso - rho*ea, (float)0.);
129	>	float iso = gamma.photonIso;
130	>	iso = max(iso - rho*ea, (float)0.);
131
132	<	return iso;
132	>	return iso;
133		}
134
135	+	float d0correction( const susy::Electron& electron, const susy::Event& event ) {
136	+	TVector3 beamspot = event.vertices[0].position;
137	+	susy::Track track = event.tracks[electron.gsfTrackIndex];
138	+	float d0 = track.d0() - beamspot.X()*sin(track.phi()) + beamspot.Y()*cos(track.phi());
139	+	return d0;
140	+	}
141
142	+	float dZcorrection( const susy::Electron& electron, const susy::Event& event ) {
143	+	TVector3 beamspot = event.vertices[0].position;
144	+	susy::Track track = event.tracks[electron.gsfTrackIndex];
145	+
146	+	if(track.momentum.Pt() == 0.) return 1.e6;
147	+	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());
148	+	return dz;
149	+	}
150
151	<
117	<	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
151	>	float getPtFromMatchedJet( susy::Photon myPhoton, susy::PFJetCollection jetColl, int loggingVerbosity = 0 ) {
152		/**
153		* \brief Takes jet p_T as photon p_T
154		*
#	Line 126 | Line 160 | float TreeWriter::getPtFromMatchedJet( s
160		std::vector<susy::PFJet> nearJets;
161		nearJets.clear();
162
163	<	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
164	<	if(pfJets_it == myEvent.pfJets.end()){
165	<	if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
166	<	} else {
167	<	susy::PFJetCollection& jetColl = pfJets_it->second;
168	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
169	<	it != jetColl.end(); it++) {
170	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
171	<	if (s_it == it->jecScaleFactors.end()) {
172	<	std::cout << "JEC is not available for this jet!!!" << std::endl;
173	<	continue;
174	<	}
175	<	float scale = s_it->second;
176	<	TLorentzVector corrP4 = scale * it->momentum;
177	<	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
144	<	if (deltaR_ > 0.3) continue;
145	<	if( loggingVerbosity > 0 )
146	<	std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
147	<	nearJets.push_back( *it );
148	<	}// for jet
149	<	}// if, else
163	>	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
164	>	it != jetColl.end(); ++it) {
165	>	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
166	>	if (s_it == it->jecScaleFactors.end()) {
167	>	std::cout << "JEC is not available for this jet!!!" << std::endl;
168	>	continue;
169	>	}
170	>	float scale = s_it->second;
171	>	TLorentzVector corrP4 = scale * it->momentum;
172	>	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
173	>	if (deltaR_ > 0.3) continue;
174	>	if( loggingVerbosity > 2 )
175	>	std::cout << " pT_jet / pT_gamma = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
176	>	nearJets.push_back( *it );
177	>	}// for jet
178
179		if ( nearJets.size() == 0 ) {
180	<	//std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
180	>	if( loggingVerbosity > 1 )
181	>	std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
182		return myPhoton.momentum.Et();
183		}
184
185		float pt = 0;
186		float minPtDifferenz = 1E20; // should be very high
187		for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
188	<	it != jetEnd; it++ ) {
188	>	it != jetEnd; ++it ) {
189		float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
190		if (  ptDiff < minPtDifferenz ) {
191		minPtDifferenz = ptDiff;
#	Line 179 | Line 208 | void TreeWriter::Loop() {
208
209		// get number of events to be proceeded
210		Long64_t nentries = inputTree->GetEntries();
211	+	// store them in histo
212	+	eventNumbers->Fill( "Number of generated Events", nentries );
213		if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;
214
215		if( loggingVerbosity > 0 )
216		std::cout << "Processing " << processNEvents << " ouf of "
217		<< nentries << " events. " << std::endl;
218
188	–	tree::Photon *thisphoton = new tree::Photon();
189	–	tree::Jet *thisjet = new tree::Jet();
190	–
219		tree->Branch("photon", &photon);
220		tree->Branch("jet", &jet);
221	+	tree->Branch("electron", &electron);
222	+	tree->Branch("muon", &muon);
223		tree->Branch("met", &met, "met/F");
224	+	tree->Branch("metPhi", &met_phi, "metPhi/F");
225	+	tree->Branch("type1met", &type1met, "type1met/F");
226	+	tree->Branch("type1metPhi", &type1met_phi, "type1metPhi/F");
227		tree->Branch("ht", &ht, "ht/F");
228		tree->Branch("nVertex", &nVertex, "nVertex/I");
229	<	tree->Branch("nElectron", &nElectron, "nElectron/I");
229	>	tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
230
231
232	<	for (long jentry=0; jentry < processNEvents; jentry++) {
233	<	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
232	>	for (long jentry=0; jentry < processNEvents; ++jentry) {
233	>	if ( loggingVerbosity>2 || jentry%reportEvery==0 )
234		std::cout << jentry << " / " << processNEvents << std :: endl;
202	–	inputTree->GetEntry(jentry);
235
236		photon.clear();
237		jet.clear();
238	+	electron.clear();
239	+	muon.clear();
240		ht = 0;
241
242	+	// weights
243	+	if (pileupHisto == 0) {
244	+	pu_weight = 1.;
245	+	} else {
246	+	float trueNumInteractions = -1;
247	+	for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
248	+	iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) {
249	+	if (iBX->BX == 0)
250	+	trueNumInteractions = iBX->trueNumInteractions;
251	+	}
252	+	pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
253	+	}
254	+
255	+	// get ak5 jets
256	+	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;
261	+
262		// photons
263		if( loggingVerbosity > 1 )
264		std::cout << "Process photons" << std::endl;
265	<	std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
266	<	for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
267	<	it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
268	<	if( ! it->isEB() && skim )
269	<	continue;
270	<	thisphoton->pt = getPtFromMatchedJet( *it, *event );
271	<	if( thisphoton->pt < 80 && skim )
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;
270	>
271	>	for(susy::PhotonCollection :: iterator it = photonVector.begin();
272	>	it != photonVector.end(); ++it ) {
273	>	if( !(it->isEB() || it->isEE()) && skim )
274		continue;
275	<	thisphoton->eta = it->momentum.Eta();
276	<	thisphoton->phi = it->momentum.Phi();
277	<	thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
278	<	thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
279	<	thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
280	<	if ( it->r9 >= 1 && skim )
275	>	tree::Photon thisphoton;
276	>	thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
277	>
278	>	thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
279	>	thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
280	>	thisphoton.photonIso = photonIso_corrected(*it, event->rho25);
281	>
282	>	bool loose_photon_barrel = thisphoton.pt>20
283	>	&& it->isEB()
284	>	&& it->passelectronveto
285	>	&& it->hadTowOverEm<0.05
286	>	&& it->sigmaIetaIeta<0.012
287	>	&& thisphoton.chargedIso<2.6
288	>	&& thisphoton.neutralIso<3.5+0.04*thisphoton.pt
289	>	&& thisphoton.photonIso<1.3+0.005*thisphoton.pt;
290	>	bool loose_photon_endcap = thisphoton.pt > 20
291	>	&& it->isEE()
292	>	&& it->passelectronveto
293	>	&& it->hadTowOverEm<0.05
294	>	&& it->sigmaIetaIeta<0.034
295	>	&& thisphoton.chargedIso<2.3
296	>	&& thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
297	>	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim )
298		continue;
299	<	thisphoton->r9 = it->r9;
300	<	thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
301	<	thisphoton->hadTowOverEm = it->hadTowOverEm;
302	<	thisphoton->pixelseed = it->nPixelSeeds;
303	<	photon.push_back( *thisphoton );
304	<	ht += thisphoton->pt;
299	>	thisphoton.eta = it->momentum.Eta();
300	>	thisphoton.phi = it->momentum.Phi();
301	>	thisphoton.r9 = it->r9;
302	>	thisphoton.sigmaIetaIeta = it->sigmaIetaIeta;
303	>	thisphoton.hadTowOverEm = it->hadTowOverEm;
304	>	thisphoton.pixelseed = it->nPixelSeeds;
305	>	thisphoton.conversionSafeVeto = it->passelectronveto;
306	>	photon.push_back( thisphoton );
307		if( loggingVerbosity > 2 )
308	<	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
308	>	std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
309		}
310
311		if( photon.size() == 0 && skim )
#	Line 248 | Line 323 | void TreeWriter::Loop() {
323		susy::PFJetCollection& jetColl = pfJets_it->second;
324
325		for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
326	<	it != jetColl.end(); it++) {
327	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
326	>	it != jetColl.end(); ++it) {
327	>	tree::Jet thisjet;
328	>	/*std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
329		if (s_it == it->jecScaleFactors.end()) {
330		std::cout << "JEC is not available for this jet!!!" << std::endl;
331		continue;
332		}
333		float scale = s_it->second;
334	+	*/
335	+	float scale = 1.;
336	+	if(it->jecScaleFactors.count("L2L3") == 0)
337	+	std::cout << "ERROR: JEC is not available for this jet" << std::endl;
338	+	else
339	+	scale = it->jecScaleFactors.find("L2L3")->second;
340		TLorentzVector corrP4 = scale * it->momentum;
341
342		if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
343	<	if(corrP4.Et() < 30 && skim ) continue;
344	<	thisjet->pt = corrP4.Et();
345	<	thisjet->eta = corrP4.Eta();
346	<	thisjet->phi = corrP4.Phi();
347	<	jet.push_back( *thisjet );
348	<	ht += thisjet->pt;
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;
359	>
360	>	if( loggingVerbosity > 2 )
361	>	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
362	>
363	>	jet.push_back( thisjet );
364	>	ht += thisjet.pt;
365		}// for jet
366		}// if, else
367		if( jet.size() < 2 && skim )
368		continue;
369		std::sort( jet.begin(), jet.end(), tree::EtGreater);
370	+	if( loggingVerbosity > 1 )
371	+	std::cout << "Found " << jet.size() << " jets" << std::endl;
372	+
373
374		// met
375		std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
376		susy::MET* metobj = &(met_it->second);
377		met = metobj->met();
378	+	met_phi = metobj->mEt.Phi();
379	+	if( loggingVerbosity > 2 )
380	+	std::cout << " met = " << met << std::endl;
381	+
382	+	std::map<TString, susy::MET>::iterator type1met_it = event->metMap.find("pfType1CorrectedMet");
383	+	susy::MET* type1metobj = &(type1met_it->second);
384	+	type1met = type1metobj->met();
385	+	type1met_phi = type1metobj->mEt.Phi();
386	+	if( loggingVerbosity > 2 )
387	+	std::cout << " type1met = " << type1met << std::endl;
388
389		// electrons
390	<	//std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
391	<	//nElectron = eVector.size();
390	>	std::vector<susy::Electron>* eVector;
391	>	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) {
398	>	tree::Particle thiselectron;
399	>	if( loggingVerbosity > 2 )
400	>	cout << " electron pt = " << it->momentum.Pt() << endl;
401	>	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
402	>	// use veto electrons
403	>	if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
404	>	continue; // spike rejection
405	>	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
406	>	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
407	>	) / it->momentum.Pt();
408	>	cout << iso << endl;
409	>	float d0 = d0correction( *it, *event );
410	>	cout << d0 << endl;
411	>	float dZ = std::abs( dZcorrection( *it, *event ) );
412	>	cout << dZ << endl;
413	>	if ( it->isEB() ){
414	>	cout << "electron in barrel" << endl;
415	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
416	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
417	>	|| it->sigmaIetaIeta > 0.01
418	>	|| it->hcalOverEcalBc > 0.15
419	>	|| d0 > 0.04
420	>	|| dZ > 0.2
421	>	|| iso > 0.15 ){
422	>	cout << " no real electron" << endl;
423	>	continue;
424	>	}
425	>	}
426	>	else if( it->isEE() ) {
427	>	cout << "electron in endcap" << endl;
428	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
429	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
430	>	|| it->sigmaIetaIeta > 0.03
431	>	|| d0 > 0.04
432	>	|| dZ > 0.2
433	>	|| iso > 0.15 ){
434	>	cout << "no real electron" << endl;
435	>	continue;
436	>	}
437	>	}
438	>	else{ // not in barrel nor in endcap
439	>	cout << " electron in gap" << endl;
440	>	continue;
441	>	// TODO: conversion rejection information not implemented yet, see twiki for more details
442	>	}
443	>	thiselectron.pt = it->momentum.Pt();
444	>	if( loggingVerbosity > 2 )
445	>	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
446	>	thiselectron.eta = it->momentum.Eta();
447	>	thiselectron.phi = it->momentum.Phi();
448	>	cout << " adde electron" << endl;
449	>	electron.push_back( thiselectron );
450	>	}
451	>	if( loggingVerbosity > 1 )
452	>	std::cout << "Found " << electron.size() << " electrons" << std::endl;
453	>
454	>	// muons
455	>	std::vector<susy::Muon> mVector = event->muons["muons"];
456	>	tree::Particle thismuon;
457	>	for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
458	>	if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
459	>	continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
460	>	thismuon.pt = it->momentum.Et();
461	>	if( thismuon.pt < 20 )
462	>	continue;
463	>	thismuon.eta = it->momentum.Eta();
464	>	thismuon.phi = it->momentum.Phi();
465	>	muon.push_back( thismuon );
466	>	}
467	>	if( loggingVerbosity > 1 )
468	>	std::cout << "Found " << muon.size() << " muons" << std::endl;
469
282	–	// vertices
470
471	+	// vertices
472		nVertex = event->vertices.size();
473
474		if( ht < 450 && skim)
475		continue;
476
477	+
478		tree->Fill();
479		} // for jentry
480
481
293	–
294	–	tree->Write();
482		outFile->cd();
483	+	eventNumbers->Write();
484	+	tree->Write();
485		outFile->Write();
486		outFile->Close();
487		}

Diff Legend

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