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.13 by kiesel, Tue Apr 16 18:53:55 2013 UTC vs.
Revision 1.18 by kiesel, Tue Apr 23 11:32:05 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	+	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 TreeWriter::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 76 | 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 94 | 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 112 | 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 130 | Line 147 | float photonIso_corrected(susy::Photon g
147		return iso;
148		}
149
150	<	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
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());
155	>	return d0;
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	>
163	>	if(track.momentum.Pt() == 0.) return 1.e6;
164	>	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());
165	>	return dz;
166	>	}
167	>
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 142 | Line 177 | float TreeWriter::getPtFromMatchedJet( s
177		std::vector<susy::PFJet> nearJets;
178		nearJets.clear();
179
180	<	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
181	<	if(pfJets_it == myEvent.pfJets.end()){
182	<	if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
183	<	} else {
184	<	susy::PFJetCollection& jetColl = pfJets_it->second;
185	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
186	<	it != jetColl.end(); ++it) {
187	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
188	<	if (s_it == it->jecScaleFactors.end()) {
189	<	std::cout << "JEC is not available for this jet!!!" << std::endl;
190	<	continue;
191	<	}
192	<	float scale = s_it->second;
193	<	TLorentzVector corrP4 = scale * it->momentum;
194	<	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
195	<	if (deltaR_ > 0.3) continue;
196	<	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
180	>	for(std::vector<susy::PFJet>::const_iterator it = jetColl.begin();
181	>	it != jetColl.end(); ++it) {
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	>	}
190	>	TLorentzVector corrP4 = scale * it->momentum;
191	>	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
192	>	if (deltaR_ > 0.3) continue;
193	>	if( loggingVerbosity > 2 )
194	>	std::cout << " pT_jet / pT_gamma = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
195	>	nearJets.push_back( *it );
196	>	}// for jet
197
198		if ( nearJets.size() == 0 ) {
199		if( loggingVerbosity > 1 )
#	Line 190 | Line 221 | float TreeWriter::getPtFromMatchedJet( s
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;
234
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 );
239		if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;
240
241		if( loggingVerbosity > 0 )
242		std::cout << "Processing " << processNEvents << " ouf of "
243		<< nentries << " events. " << std::endl;
244
205	–	tree::Photon *thisphoton = new tree::Photon();
206	–	tree::Jet *thisjet = new tree::Jet();
207	–
245		tree->Branch("photon", &photon);
246		tree->Branch("jet", &jet);
247		tree->Branch("electron", &electron);
#	Line 218 | 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>0 && 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->GetEntry(jentry);
261	>	inputTree->LoadTree( jentry );
262	>	inputTree->GetEntry( jentry );
263
264		photon.clear();
265		jet.clear();
#	Line 242 | Line 280 | void TreeWriter::Loop() {
280		pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
281		}
282
283	+	// get ak5 jets
284	+	std::vector<susy::PFJet> jetVector = event->pfJets["ak5"];
285
286		// photons
287	<	if( loggingVerbosity > 1 )
288	<	std::cout << "Process photons" << std::endl;
289	<	std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
290	<	for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
251	<	it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
287	>	std::vector<susy::Photon> photonVector = event->photons["photons"];
288	>
289	>	for(std::vector<susy::Photon>::iterator it = photonVector.begin();
290	>	it != photonVector.end(); ++it ) {
291		if( !(it->isEB() || it->isEE()) && skim )
292		continue;
293	<	thisphoton->pt = getPtFromMatchedJet( *it, *event );
293	>	tree::Photon thisphoton;
294	>	thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
295
296	<	thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
297	<	thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
298	<	thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
296	>	thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
297	>	thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
298	>	thisphoton.photonIso = photonIso_corrected(*it, event->rho25);
299
300	<	bool loose_photon_barrel = thisphoton->pt>20
300	>	bool loose_photon_barrel = thisphoton.pt>20
301		&& it->isEB()
302		&& it->passelectronveto
303		&& it->hadTowOverEm<0.05
304		&& it->sigmaIetaIeta<0.012
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
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
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 )
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 )
316		continue;
317	<	thisphoton->eta = it->momentum.Eta();
318	<	thisphoton->phi = it->momentum.Phi();
319	<	thisphoton->r9 = it->r9;
320	<	thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
321	<	thisphoton->hadTowOverEm = it->hadTowOverEm;
322	<	thisphoton->pixelseed = it->nPixelSeeds;
323	<	thisphoton->conversionSafeVeto = it->passelectronveto;
324	<	photon.push_back( *thisphoton );
317	>	thisphoton.eta = it->momentum.Eta();
318	>	thisphoton.phi = it->momentum.Phi();
319	>	thisphoton.r9 = it->r9;
320	>	thisphoton.sigmaIetaIeta = it->sigmaIetaIeta;
321	>	thisphoton.hadTowOverEm = it->hadTowOverEm;
322	>	thisphoton.pixelseed = it->nPixelSeeds;
323	>	thisphoton.conversionSafeVeto = it->passelectronveto;
324	>	photon.push_back( thisphoton );
325		if( loggingVerbosity > 2 )
326	<	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
326	>	std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
327		}
328
329		if( photon.size() == 0 && skim )
#	Line 293 | Line 333 | void TreeWriter::Loop() {
333		std::cout << "Found " << photon.size() << " photons" << std::endl;
334
335		// jets
296	–	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
297	–	if(pfJets_it == event->pfJets.end()){
298	–	if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
299	–	} else {
336
301	–	susy::PFJetCollection& jetColl = pfJets_it->second;
337
338	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
339	<	it != jetColl.end(); ++it) {
340	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
341	<	if (s_it == it->jecScaleFactors.end()) {
342	<	std::cout << "JEC is not available for this jet!!!" << std::endl;
343	<	continue;
344	<	}
345	<	float scale = s_it->second;
346	<	TLorentzVector corrP4 = scale * it->momentum;
338	>	for(std::vector<susy::PFJet>::iterator it = jetVector.begin();
339	>	it != jetVector.end(); ++it) {
340	>	tree::Jet thisjet;
341	>
342	>	// scale with JEC
343	>	float scale = 1.;
344	>	if(it->jecScaleFactors.count("L2L3") == 0)
345	>	std::cout << "ERROR: JEC is not available for this jet" << std::endl;
346	>	else
347	>	scale = it->jecScaleFactors.find("L2L3")->second;
348	>	TLorentzVector corrP4 = scale * it->momentum;
349	>
350	>	// Calculate HT.
351	>	// The definiton differs from the saved jet, since trigger is described better
352	>	if( std::abs( corrP4.Eta() ) < 3 && corrP4.Pt > 40 )
353	>	ht += thisjet.pt;
354	>
355	>	if(std::abs(corrP4.Eta()) > 2.6 && skim ) continue;
356	>	if(corrP4.Pt() < 30 && skim ) continue;
357	>	thisjet.pt = corrP4.Pt();
358	>	thisjet.eta = corrP4.Eta();
359	>	thisjet.phi = corrP4.Phi();
360	>	thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
361	>	// jet composition
362	>	thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
363	>	thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
364	>	thisjet.photonEnergy = it->photonEnergy;
365	>	thisjet.electronEnergy = it->electronEnergy;
366	>	thisjet.muonEnergy = it->muonEnergy;
367	>	thisjet.HFHadronEnergy = it->HFHadronEnergy;
368	>	thisjet.HFEMEnergy = it->HFEMEnergy;
369	>	thisjet.chargedEmEnergy = it->chargedEmEnergy;
370	>	thisjet.chargedMuEnergy = it->chargedMuEnergy;
371	>	thisjet.neutralEmEnergy = it->neutralEmEnergy;
372
373	<	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
374	<	if(corrP4.Et() < 30 && skim ) continue;
375	<	thisjet->pt = corrP4.Et();
376	<	thisjet->eta = corrP4.Eta();
377	<	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	<
334	<	jet.push_back( *thisjet );
335	<	ht += thisjet->pt;
336	<	}// for jet
337	<	}// if, else
373	>	if( loggingVerbosity > 2 )
374	>	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
375	>
376	>	jet.push_back( thisjet );
377	>	}// for jet
378		if( jet.size() < 2 && skim )
379		continue;
380		std::sort( jet.begin(), jet.end(), tree::EtGreater);
#	Line 358 | Line 398 | void TreeWriter::Loop() {
398		std::cout << " type1met = " << type1met << std::endl;
399
400		// electrons
401	<	tree::Particle* thiselectron = new tree::Particle();
402	<	map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
403	<	if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
404	<	cout << "gsfElectrons not found!" << endl;
405	<	} else {
406	<	for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
407	<	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
408	<	if( loggingVerbosity > 0)
409	<	std::cout << "TODO: piAtVtx is not calculated correctly" << std::endl;
410	<	//float pAtVtx = it->trackMomentums.find("AtVtxWithConstraint")->second.P();
411	<	float pAtVtx = it->ecalEnergy;
412	<	if( it->momentum.Pt() < 1  || it->momentum.Pt() > 1e6 || pAtVtx == 0 )
401	>	std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
402	>	for(std::vector<susy::Electron>::iterator it = eVector.begin(); it < eVector.end(); ++it) {
403	>	tree::Particle thiselectron;
404	>	if( loggingVerbosity > 2 )
405	>	cout << " electron pt = " << it->momentum.Pt() << endl;
406	>	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
407	>	// use veto electrons
408	>	if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
409	>	continue; // spike rejection
410	>	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
411	>	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
412	>	) / it->momentum.Pt();
413	>	float d0 = d0correction( *it, *event );
414	>	float dZ = std::abs( dZcorrection( *it, *event ) );
415	>	if ( it->isEB() ){
416	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
417	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
418	>	|| it->sigmaIetaIeta > 0.01
419	>	|| it->hcalOverEcalBc > 0.15
420	>	|| d0 > 0.04
421	>	|| dZ > 0.2
422	>	|| iso > 0.15 )
423		continue;
424	<	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
425	<	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
426	<	) / it->momentum.Pt();
427	<	if ( it->isEE() ){
428	<	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
429	<	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.15
430	<	|| it->sigmaIetaIeta > 0.01
431	<	|| 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
424	>	}
425	>	else if( it->isEE() ) {
426	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
427	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
428	>	|| it->sigmaIetaIeta > 0.03
429	>	|| d0 > 0.04
430	>	|| dZ > 0.2
431	>	|| iso > 0.15 )
432		continue;
433	<	// TODO: conversion rejection information not implemented yet, see twiki for more details
433	>	}
434	>	else // not in barrel nor in endcap
435	>	continue;
436
437	<	thiselectron->pt = it->momentum.Pt();
438	<	if( thiselectron->pt < 20 )
439	<	continue;
440	<	if( loggingVerbosity > 2 )
441	<	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
442	<	thiselectron->eta = it->momentum.Eta();
409	<	thiselectron->phi = it->momentum.Phi();
410	<	electron.push_back( *thiselectron );
411	<	}
437	>	thiselectron.pt = it->momentum.Pt();
438	>	if( loggingVerbosity > 2 )
439	>	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
440	>	thiselectron.eta = it->momentum.Eta();
441	>	thiselectron.phi = it->momentum.Phi();
442	>	electron.push_back( thiselectron );
443		}
444		if( loggingVerbosity > 1 )
445		std::cout << "Found " << electron.size() << " electrons" << std::endl;
446
447		// muons
448	+	tree::Particle thismuon;
449		std::vector<susy::Muon> mVector = event->muons["muons"];
418	–	tree::Particle* thismuon = new tree::Particle();
450		for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
451		if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
452		continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
453	<	thismuon->pt = it->momentum.Et();
454	<	if( thismuon->pt < 20 )
453	>	thismuon.pt = it->momentum.Et();
454	>	if( thismuon.pt < 20 )
455		continue;
456	<	thismuon->eta = it->momentum.Eta();
457	<	thismuon->phi = it->momentum.Phi();
458	<	muon.push_back( *thismuon );
456	>	thismuon.eta = it->momentum.Eta();
457	>	thismuon.phi = it->momentum.Phi();
458	>	muon.push_back( thismuon );
459		}
460		if( loggingVerbosity > 1 )
461		std::cout << "Found " << muon.size() << " muons" << std::endl;
462
432	–
463		// vertices
464		nVertex = event->vertices.size();
465
466		if( ht < 450 && skim)
467		continue;
468
469	+
470		tree->Fill();
471		} // for jentry
472
473
474		outFile->cd();
475	+	eventNumbers->Write();
476		tree->Write();
477		outFile->Write();
478		outFile->Close();

Diff Legend

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