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.14 by kiesel, Wed Apr 17 09:57:52 2013 UTC vs.
Revision 1.26 by kiesel, Thu Apr 25 17:37:32 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;
#	Line 31 | Line 24 | void TreeWriter::Init( TString outputNam
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");
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 57 | 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	+	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 TreeWriter::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 80 | 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 98 | 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 116 | 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 134 | Line 154 | float photonIso_corrected(susy::Photon g
154		return iso;
155		}
156
157	<	float d0correction( susy::Electron electron, susy::Event event ) {
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());
141	–	cout << "return " << d0 << endl;
162		return d0;
163		}
164
165	<	float dZcorrection( susy::Electron electron, susy::Event event ) {
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 151 | Line 172 | float dZcorrection( susy::Electron elect
172		return dz;
173		}
174
175	<	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
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 163 | Line 184 | float TreeWriter::getPtFromMatchedJet( s
184		std::vector<susy::PFJet> nearJets;
185		nearJets.clear();
186
187	<	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
188	<	if(pfJets_it == myEvent.pfJets.end()){
189	<	if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
190	<	} else {
191	<	susy::PFJetCollection& jetColl = pfJets_it->second;
192	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
193	<	it != jetColl.end(); ++it) {
194	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
195	<	if (s_it == it->jecScaleFactors.end()) {
196	<	std::cout << "JEC is not available for this jet!!!" << std::endl;
197	<	continue;
198	<	}
199	<	float scale = s_it->second;
200	<	TLorentzVector corrP4 = scale * it->momentum;
201	<	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
202	<	if (deltaR_ > 0.3) continue;
203	<	if( loggingVerbosity > 0 )
183	<	std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
184	<	nearJets.push_back( *it );
185	<	}// for jet
186	<	}// if, else
187	>	for(std::vector<susy::PFJet>::const_iterator it = jetColl.begin();
188	>	it != jetColl.end(); ++it) {
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	>	}
197	>	TLorentzVector corrP4 = scale * it->momentum;
198	>	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
199	>	if (deltaR_ > 0.3) continue;
200	>	if( loggingVerbosity > 2 )
201	>	std::cout << " pT_jet / pT_gamma = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
202	>	nearJets.push_back( *it );
203	>	}// for jet
204
205		if ( nearJets.size() == 0 ) {
206		if( loggingVerbosity > 1 )
#	Line 211 | Line 228 | float TreeWriter::getPtFromMatchedJet( s
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 218 | Line 242 | void TreeWriter::Loop() {
242		// get number of events to be proceeded
243		Long64_t nentries = inputTree->GetEntries();
244		// store them in histo
245	<	eventNumbers->Fill( "Number of generated Events", nentries );
245	>	eventNumbers->Fill( "Number of generated events", nentries );
246		if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;
247
248		if( loggingVerbosity > 0 )
249		std::cout << "Processing " << processNEvents << " ouf of "
250		<< nentries << " events. " << std::endl;
251
228	–	tree::Photon *thisphoton = new tree::Photon();
229	–	tree::Jet *thisjet = new tree::Jet();
230	–
252		tree->Branch("photon", &photon);
253		tree->Branch("jet", &jet);
254		tree->Branch("electron", &electron);
#	Line 239 | Line 260 | void TreeWriter::Loop() {
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	+	tree->Branch("genElectron", &genElectron);
264	+	tree->Branch("genPhoton", &genPhoton);
265
266	<
267	<	for (long jentry=0; jentry < processNEvents; ++jentry) {
245	<	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
266	>	for (unsigned long jentry=0; jentry < processNEvents; ++jentry) {
267	>	if ( loggingVerbosity>1 || jentry%reportEvery==0 )
268		std::cout << jentry << " / " << processNEvents << std :: endl;
269	<	inputTree->GetEntry(jentry);
269	>	inputTree->LoadTree( jentry );
270	>	inputTree->GetEntry( jentry );
271
272		photon.clear();
273		jet.clear();
274		electron.clear();
275		muon.clear();
276	+	genElectron.clear();
277	+	genPhoton.clear();
278		ht = 0;
279
280		// weights
#	Line 265 | Line 290 | void TreeWriter::Loop() {
290		pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
291		}
292
293	+	// get ak5 jets
294	+	std::vector<susy::PFJet> jetVector = event->pfJets["ak5"];
295
296		// photons
297	<	if( loggingVerbosity > 1 )
298	<	std::cout << "Process photons" << std::endl;
299	<	std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
300	<	for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
301	<	it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
275	<	if( !(it->isEB() || it->isEE()) && skim )
297	>	std::vector<susy::Photon> photonVector = event->photons["photons"];
298	>
299	>	for(std::vector<susy::Photon>::iterator it = photonVector.begin();
300	>	it != photonVector.end(); ++it ) {
301	>	if( !(it->isEE() || it->isEB()) && it->isEBEtaGap() && it->isEBPhiGap() && it->isEERingGap() && it->isEEDeeGap() && it->isEBEEGap() && skim )
302		continue;
303	<	thisphoton->pt = getPtFromMatchedJet( *it, *event );
303	>	tree::Photon thisphoton;
304	>	thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
305
306	<	thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
307	<	thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
308	<	thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
306	>	thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
307	>	thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
308	>	thisphoton.photonIso = photonIso_corrected(*it, event->rho25);
309
310	<	bool loose_photon_barrel = thisphoton->pt>20
310	>	bool loose_photon_barrel = thisphoton.pt>20
311		&& it->isEB()
312		&& it->passelectronveto
313		&& it->hadTowOverEm<0.05
314		&& it->sigmaIetaIeta<0.012
315	<	&& thisphoton->chargedIso<2.6
316	<	&& thisphoton->neutralIso<3.5+0.04*thisphoton->pt
317	<	&& thisphoton->photonIso<1.3+0.005*thisphoton->pt;
318	<	bool loose_photon_endcap = thisphoton->pt > 20
315	>	&& thisphoton.chargedIso<2.6
316	>	&& thisphoton.neutralIso<3.5+0.04*thisphoton.pt
317	>	&& thisphoton.photonIso<1.3+0.005*thisphoton.pt;
318	>	bool loose_photon_endcap = thisphoton.pt > 20
319		&& it->isEE()
320		&& it->passelectronveto
321		&& it->hadTowOverEm<0.05
322		&& it->sigmaIetaIeta<0.034
323	<	&& thisphoton->chargedIso<2.3
324	<	&& thisphoton->neutralIso<2.9+0.04*thisphoton->pt;
325	<	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton->pt > 75 ) && skim )
323	>	&& thisphoton.chargedIso<2.3
324	>	&& thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
325	>
326	>	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim )
327		continue;
328	<	thisphoton->eta = it->momentum.Eta();
329	<	thisphoton->phi = it->momentum.Phi();
330	<	thisphoton->r9 = it->r9;
331	<	thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
332	<	thisphoton->hadTowOverEm = it->hadTowOverEm;
333	<	thisphoton->pixelseed = it->nPixelSeeds;
334	<	thisphoton->conversionSafeVeto = it->passelectronveto;
335	<	photon.push_back( *thisphoton );
328	>	thisphoton.eta = it->momentum.Eta();
329	>	thisphoton.phi = it->momentum.Phi();
330	>	thisphoton.r9 = it->r9;
331	>	thisphoton.sigmaIetaIeta = it->sigmaIetaIeta;
332	>	thisphoton.hadTowOverEm = it->hadTowOverEm;
333	>	thisphoton.pixelseed = it->nPixelSeeds;
334	>	thisphoton.conversionSafeVeto = it->passelectronveto;
335	>	photon.push_back( thisphoton );
336		if( loggingVerbosity > 2 )
337	<	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
337	>	std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
338		}
339
340		if( photon.size() == 0 && skim )
#	Line 316 | Line 344 | void TreeWriter::Loop() {
344		std::cout << "Found " << photon.size() << " photons" << std::endl;
345
346		// 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 {
347
324	–	susy::PFJetCollection& jetColl = pfJets_it->second;
348
349	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
350	<	it != jetColl.end(); ++it) {
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	<	TLorentzVector corrP4 = scale * it->momentum;
349	>	for(std::vector<susy::PFJet>::iterator it = jetVector.begin();
350	>	it != jetVector.end(); ++it) {
351	>	tree::Jet thisjet;
352	>
353	>	// scale with JEC
354	>	float scale = 1.;
355	>	if(it->jecScaleFactors.count("L2L3") == 0)
356	>	std::cout << "ERROR: JEC is not available for this jet" << std::endl;
357	>	else
358	>	scale = it->jecScaleFactors.find("L2L3")->second;
359	>	TLorentzVector corrP4 = scale * it->momentum;
360
361	<	if(std::fabs(corrP4.Eta()) > 3.0 && skim ) continue;
362	<	if(corrP4.Et() < 30 && skim ) continue;
363	<	thisjet->pt = corrP4.Et();
364	<	thisjet->eta = corrP4.Eta();
365	<	thisjet->phi = corrP4.Phi();
366	<	thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
367	<	// jet composition
368	<	thisjet->chargedHadronEnergy = it->chargedHadronEnergy;
369	<	thisjet->neutralHadronEnergy = it->neutralHadronEnergy;
370	<	thisjet->photonEnergy = it->photonEnergy;
371	<	thisjet->electronEnergy = it->electronEnergy;
372	<	thisjet->muonEnergy = it->muonEnergy;
373	<	thisjet->HFHadronEnergy = it->HFHadronEnergy;
374	<	thisjet->HFEMEnergy = it->HFEMEnergy;
375	<	thisjet->chargedEmEnergy = it->chargedEmEnergy;
376	<	thisjet->chargedMuEnergy = it->chargedMuEnergy;
377	<	thisjet->neutralEmEnergy = it->neutralEmEnergy;
378	<
379	<	if( loggingVerbosity > 2 )
380	<	std::cout << " p_T, jet = " << thisjet->pt << std::endl;
381	<
382	<	jet.push_back( *thisjet );
383	<	ht += thisjet->pt;
384	<	}// for jet
385	<	}// if, else
361	>	// Calculate HT.
362	>	// The definiton differs from the saved jet, since trigger is described better
363	>	if( std::abs( corrP4.Eta() ) < 3 && corrP4.Pt() > 40 )
364	>	ht += thisjet.pt;
365	>
366	>	if(std::abs(corrP4.Eta()) > 2.6 && skim ) continue;
367	>	if(corrP4.Pt() < 30 && skim ) continue;
368	>	thisjet.pt = corrP4.Pt();
369	>	thisjet.eta = corrP4.Eta();
370	>	thisjet.phi = corrP4.Phi();
371	>	thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
372	>	// jet composition
373	>	thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
374	>	thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
375	>	thisjet.photonEnergy = it->photonEnergy;
376	>	thisjet.electronEnergy = it->electronEnergy;
377	>	thisjet.muonEnergy = it->muonEnergy;
378	>	thisjet.HFHadronEnergy = it->HFHadronEnergy;
379	>	thisjet.HFEMEnergy = it->HFEMEnergy;
380	>	thisjet.chargedEmEnergy = it->chargedEmEnergy;
381	>	thisjet.chargedMuEnergy = it->chargedMuEnergy;
382	>	thisjet.neutralEmEnergy = it->neutralEmEnergy;
383	>
384	>	if( loggingVerbosity > 2 )
385	>	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
386	>
387	>	jet.push_back( thisjet );
388	>	}// for jet
389		if( jet.size() < 2 && skim )
390		continue;
391		std::sort( jet.begin(), jet.end(), tree::EtGreater);
392		if( loggingVerbosity > 1 )
393		std::cout << "Found " << jet.size() << " jets" << std::endl;
394
395	+	if( ht < 450 && skim)
396	+	continue;
397	+
398	+
399
400		// met
401		std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
#	Line 381 | Line 413 | void TreeWriter::Loop() {
413		std::cout << " type1met = " << type1met << std::endl;
414
415		// electrons
416	<	tree::Particle* thiselectron = new tree::Particle();
417	<	map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
418	<	if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
419	<	cout << "gsfElectrons not found!" << endl;
420	<	} else {
421	<	cout << "now begin electrons " << endl;
422	<	for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
423	<	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification for veto electrons
424	<	if( it->momentum.Pt() < 1  || it->momentum.Pt() > 1e6 )
425	<	continue; // spike rejection
426	<	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
427	<	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
428	<	) / it->momentum.Pt();
429	<	cout << iso << endl;
430	<	cout << d0correction( *it, *event ) << endl;
431	<	float test = d0correction( *it, *event );
432	<	cout << test << endl;
433	<	float dZ = std::fabs( dZcorrection( *it, *event ) );
434	<	cout <<" find e" << endl;
435	<	float d0 = 0.;
436	<	if ( it->isEB() ){
405	<	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
406	<	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
407	<	|| it->sigmaIetaIeta > 0.01
408	<	|| it->hcalOverEcalBc > 0.15
409	<	|| d0 > 0.04
410	<	|| dZ > 0.2
411	<	|| iso > 0.15 )
412	<	continue;
413	<	}
414	<	else if( it->isEE() ) {
415	<	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
416	<	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
417	<	|| it->sigmaIetaIeta > 0.03
418	<	|| d0 > 0.04
419	<	|| dZ > 0.2
420	<	|| iso > 0.15 )
421	<	continue;
422	<	}
423	<	else // not in barrel nor in endcap
416	>	std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
417	>	for(std::vector<susy::Electron>::iterator it = eVector.begin(); it < eVector.end(); ++it) {
418	>	tree::Particle thiselectron;
419	>	if( loggingVerbosity > 2 )
420	>	cout << " electron pt = " << it->momentum.Pt() << endl;
421	>	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
422	>	// use veto electrons
423	>	if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
424	>	continue; // spike rejection
425	>	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (float)0. )
426	>	) / it->momentum.Pt();
427	>	float d0 = d0correction( *it, *event );
428	>	float dZ = std::abs( dZcorrection( *it, *event ) );
429	>	if ( it->isEB() ){
430	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
431	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
432	>	|| it->sigmaIetaIeta > 0.01
433	>	|| it->hcalOverEcalBc > 0.15
434	>	|| d0 > 0.04
435	>	|| dZ > 0.2
436	>	|| iso > 0.15 )
437		continue;
438	<	// TODO: conversion rejection information not implemented yet, see twiki for more details
439	<
440	<	thiselectron->pt = it->momentum.Pt();
441	<	if( thiselectron->pt < 20 )
438	>	}
439	>	else if( it->isEE() ) {
440	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
441	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
442	>	|| it->sigmaIetaIeta > 0.03
443	>	|| d0 > 0.04
444	>	|| dZ > 0.2
445	>	|| iso > 0.15 )
446		continue;
447	<	if( loggingVerbosity > 2 )
448	<	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
449	<	thiselectron->eta = it->momentum.Eta();
450	<	thiselectron->phi = it->momentum.Phi();
451	<	electron.push_back( *thiselectron );
452	<	}
447	>	}
448	>	else // not in barrel nor in endcap
449	>	continue;
450	>
451	>	thiselectron.pt = it->momentum.Pt();
452	>	if( loggingVerbosity > 2 )
453	>	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
454	>	thiselectron.eta = it->momentum.Eta();
455	>	thiselectron.phi = it->momentum.Phi();
456	>	electron.push_back( thiselectron );
457		}
458		if( loggingVerbosity > 1 )
459		std::cout << "Found " << electron.size() << " electrons" << std::endl;
460
461		// muons
462	+	tree::Particle thismuon;
463		std::vector<susy::Muon> mVector = event->muons["muons"];
442	–	tree::Particle* thismuon = new tree::Particle();
464		for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
465		if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
466		continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
467	<	thismuon->pt = it->momentum.Et();
468	<	if( thismuon->pt < 20 )
467	>	thismuon.pt = it->momentum.Et();
468	>	if( thismuon.pt < 20 )
469		continue;
470	<	thismuon->eta = it->momentum.Eta();
471	<	thismuon->phi = it->momentum.Phi();
472	<	muon.push_back( *thismuon );
470	>	thismuon.eta = it->momentum.Eta();
471	>	thismuon.phi = it->momentum.Phi();
472	>	muon.push_back( thismuon );
473		}
474		if( loggingVerbosity > 1 )
475		std::cout << "Found " << muon.size() << " muons" << std::endl;
476
456	–
477		// vertices
478		nVertex = event->vertices.size();
479
480	<	if( ht < 450 && skim)
481	<	continue;
480	>	tree::Particle thisGenParticle;
481	>	for( std::vector<susy::Particle>::iterator it = event->genParticles.begin(); it != event->genParticles.end(); ++it ) {
482	>	if( it->momentum.Pt() < 20 ) continue;
483	>	thisGenParticle.pt = it->momentum.Pt();
484	>	thisGenParticle.eta = it->momentum.Eta();
485	>	thisGenParticle.phi = it->momentum.Phi();
486	>	switch( std::abs(it->pdgId) ) {
487	>	case 22: // photon
488	>	genPhoton.push_back( thisGenParticle );
489	>	break;
490	>	case 11: // electron
491	>	// Demand a W boson as mother particle of electron
492	>	if( abs(event->genParticles[it->motherIndex].pdgId) == 24 )
493	>	genElectron.push_back( thisGenParticle );
494	>	break;
495	>	}
496	>	}
497
498		tree->Fill();
499		} // for jentry

Diff Legend

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