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.2 by kiesel, Wed Mar 20 11:04:30 2013 UTC vs.
Revision 1.16 by kiesel, Thu Apr 18 16:26:44 2013 UTC

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

Diff Legend

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