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.4 by kiesel, Mon Mar 25 12:00:06 2013 UTC vs.
Revision 1.14 by kiesel, Wed Apr 17 09:57:52 2013 UTC

#	Line 1 | Line 1
1		#include<iostream>
2		#include<math.h>
3	+	#include<string>
4
5	<	#include "TFile.h"
5	<	#include "TTree.h"
6	<	#include "TString.h"
7	<
8	<	#include "SusyEvent.h"
9	<	#include "TreeObjects.h"
10	<
11	<
12	<	class TreeWriter {
13	<	public :
14	<	TreeWriter(TString inputName, TString outputName );
15	<	virtual ~TreeWriter();
16	<	virtual void Loop();
17	<
18	<	void SetProcessNEvents(int nEvents) { processNEvents = nEvents; }
19	<	void SetReportEvents(int nEvents) { reportEvery = nEvents; }
20	<
21	<	TFile *inputFile;
22	<	TTree *inputTree;
23	<	susy::Event *event;
24	<
25	<	TFile *outFile;
26	<	TTree *tree;
27	<
28	<	float getPtFromMatchedJet( susy::Photon, susy::Event );
29	<	float deltaR( TLorentzVector, TLorentzVector );
30	<
31	<	private:
32	<	int processNEvents; // number of events to be processed
33	<	int reportEvery;
34	<	int loggingVerbosity;
35	<
36	<	// variables which will be stored in the tree
37	<	std::vector<tree::Photon> photon;
38	<	std::vector<tree::Jet> jet;
39	<	float met;
40	<	int nVertex;
41	<	int nElectron;
42	<	float weight;
43	<	};
5	>	#include "TSystem.h"
6
7	<	TreeWriter::TreeWriter(TString inputName, TString outputName) {
7	>	#include "treeWriter.h"
8	>
9	>	using namespace std;
10	>
11	>	TreeWriter::TreeWriter( TString inputName, TString outputName, int loggingVerbosity_ ) {
12		// read the input file
13	<	inputFile = new TFile( inputName, "read" );
14	<	inputTree = (TTree*) inputFile->Get("susyTree");
13	>	inputTree = new TChain("susyTree");
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;
30		event = new susy::Event;
31		inputTree->SetBranchAddress("susyEvent", &event);
32
33	+	// Here the number of proceeded events will be stored. For plotting, simply use L*sigma/eventNumber
34	+	eventNumbers = new TH1F("eventNumbers", "Histogram containing number of generated Events", 1, 0, 1);
35	+	eventNumbers->GetXaxis()->SetBinLabel(1,"Number of generated Events");
36	+
37		// open the output file
38	+	if (loggingVerbosity_>0)
39	+	std::cout << "Open file " << outputName << " for writing." << std::endl;
40		outFile = new TFile( outputName, "recreate" );
41		tree = new TTree("susyTree","Tree for single photon analysis");
42
43		// set default parameter
44		processNEvents = -1;
45		reportEvery = 1000;
46	<	loggingVerbosity = 0;
46	>	loggingVerbosity = loggingVerbosity_;
47	>	skim = true;
48	>	pileupHisto = 0;
49	>	}
50
51	+	void TreeWriter::PileUpWeightFile( string pileupFileName ) {
52	+	TFile *puFile = new TFile( pileupFileName.c_str() );
53	+	pileupHisto = (TH1F*) puFile->Get("pileup");
54		}
55
56		TreeWriter::~TreeWriter() {
57	<	if (!inputTree) return;
58	<	delete inputTree->GetCurrentFile();
57	>	if (pileupHisto != 0 )
58	>	delete pileupHisto;
59	>	inputTree->GetCurrentFile()->Close();
60		}
61
62	+	// useful functions
63		float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
64		return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
65		}
66
67	<	float TreeWriter::getPtFromMatchedJet( susy::Photon photon, susy::Event event ) {
67	>	float effectiveAreaElectron( float eta ) {
68	>	// see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection
69	>	// only for Delta R = 0.3 on 2012 Data
70	>	eta = fabs( eta );
71	>	float ea;
72	>	if( eta < 1.0 ) ea = 0.13;
73	>	else if( eta < 1.479 ) ea = 0.14;
74	>	else if( eta < 2.0 ) ea = 0.07;
75	>	else if( eta < 2.2 ) ea = 0.09;
76	>	else if( eta < 2.3 ) ea = 0.11;
77	>	else if( eta < 2.4 ) ea = 0.11;
78	>	else ea = 0.14;
79	>	return ea;
80	>	}
81	>
82	>	// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
83	>	float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
84	>	float eta = fabs(gamma.caloPosition.Eta());
85	>	float ea;
86	>
87	>	if(eta < 1.0) ea = 0.012;
88	>	else if(eta < 1.479) ea = 0.010;
89	>	else if(eta < 2.0) ea = 0.014;
90	>	else if(eta < 2.2) ea = 0.012;
91	>	else if(eta < 2.3) ea = 0.016;
92	>	else if(eta < 2.4) ea = 0.020;
93	>	else ea = 0.012;
94	>
95	>	float iso = gamma.chargedHadronIso;
96	>	iso = max(iso - rho*ea, (float)0.);
97	>
98	>	return iso;
99	>	}
100	>
101	>	float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
102	>	float eta = fabs(gamma.caloPosition.Eta());
103	>	float ea;
104	>
105	>	if(eta < 1.0) ea = 0.030;
106	>	else if(eta < 1.479) ea = 0.057;
107	>	else if(eta < 2.0) ea = 0.039;
108	>	else if(eta < 2.2) ea = 0.015;
109	>	else if(eta < 2.3) ea = 0.024;
110	>	else if(eta < 2.4) ea = 0.039;
111	>	else ea = 0.072;
112	>
113	>	float iso = gamma.neutralHadronIso;
114	>	iso = max(iso - rho*ea, (float)0.);
115	>
116	>	return iso;
117	>	}
118	>
119	>	float photonIso_corrected(susy::Photon gamma, float rho) {
120	>	float eta = fabs(gamma.caloPosition.Eta());
121	>	float ea;
122	>
123	>	if(eta < 1.0) ea = 0.148;
124	>	else if(eta < 1.479) ea = 0.130;
125	>	else if(eta < 2.0) ea = 0.112;
126	>	else if(eta < 2.2) ea = 0.216;
127	>	else if(eta < 2.3) ea = 0.262;
128	>	else if(eta < 2.4) ea = 0.260;
129	>	else ea = 0.266;
130	>
131	>	float iso = gamma.photonIso;
132	>	iso = max(iso - rho*ea, (float)0.);
133	>
134	>	return iso;
135	>	}
136	>
137	>	float d0correction( susy::Electron electron, susy::Event event ) {
138	>	TVector3 beamspot = event.vertices[0].position;
139	>	susy::Track track = event.tracks[electron.gsfTrackIndex];
140	>	float d0 = track.d0() - beamspot.X()*sin(track.phi()) + beamspot.Y()*cos(track.phi());
141	>	cout << "return " << d0 << endl;
142	>	return d0;
143	>	}
144	>
145	>	float dZcorrection( susy::Electron electron, susy::Event event ) {
146	>	TVector3 beamspot = event.vertices[0].position;
147	>	susy::Track track = event.tracks[electron.gsfTrackIndex];
148	>
149	>	if(track.momentum.Pt() == 0.) return 1.e6;
150	>	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());
151	>	return dz;
152	>	}
153	>
154	>	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
155		/**
156		* \brief Takes jet p_T as photon p_T
157		*
158		* At first all jets with DeltaR < 0.3 (isolation cone) are searched.
159		* If several jets are found, take the one with the minimal pt difference
160	<	* compared to the photon. If no such jets are found, return 0
160	>	* compared to the photon. If no such jets are found, keep the photon_pt
161		* TODO: remove photon matched jet from jet-selection?
162		*/
163		std::vector<susy::PFJet> nearJets;
164		nearJets.clear();
165
166	<	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event.pfJets.find("ak5");
167	<	if(pfJets_it == event.pfJets.end()){
168	<	if(event.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
166	>	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
167	>	if(pfJets_it == myEvent.pfJets.end()){
168	>	if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
169		} else {
170		susy::PFJetCollection& jetColl = pfJets_it->second;
171		for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
172	<	it != jetColl.end(); it++) {
172	>	it != jetColl.end(); ++it) {
173		std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
174		if (s_it == it->jecScaleFactors.end()) {
175		std::cout << "JEC is not available for this jet!!!" << std::endl;
#	Line 95 | Line 177 | float TreeWriter::getPtFromMatchedJet( s
177		}
178		float scale = s_it->second;
179		TLorentzVector corrP4 = scale * it->momentum;
180	<	float deltaR_ = deltaR(photon.momentum, corrP4 );
180	>	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
181		if (deltaR_ > 0.3) continue;
182	+	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
188		if ( nearJets.size() == 0 ) {
189	<	std::cout << "No jet with ΔR < .3 found, set p_T to 0" << std::endl;
190	<	return 0;
189	>	if( loggingVerbosity > 1 )
190	>	std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
191	>	return myPhoton.momentum.Et();
192		}
193
194		float pt = 0;
195		float minPtDifferenz = 1E20; // should be very high
196		for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
197	<	it != jetEnd; it++ ) {
198	<	float ptDiff = fabs(photon.momentum.Et() - it->momentum.Et());
197	>	it != jetEnd; ++it ) {
198	>	float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
199		if (  ptDiff < minPtDifferenz ) {
200		minPtDifferenz = ptDiff;
201		pt = it->momentum.Et();
#	Line 118 | Line 203 | float TreeWriter::getPtFromMatchedJet( s
203		}
204
205		// testing
206	<	if( nearJets.size() > 1 )
206	>	if( nearJets.size() > 1 && loggingVerbosity > 0 )
207		std::cout << "There are several jets matching to this photon. "
208	<	<< "Please check if the upper assumtion is reasonable" << std::endl;
208	>	<< "Please check if jet-matching is correct." << std::endl;
209		return pt;
210		}
211
212
213		void TreeWriter::Loop() {
129	–	// only for testing
130	–	bool skim = true;
214
215		// here the event loop is implemented and the tree is filled
216		if (inputTree == 0) return;
217
218		// get number of events to be proceeded
219		Long64_t nentries = inputTree->GetEntries();
220	+	// store them in histo
221	+	eventNumbers->Fill( "Number of generated Events", nentries );
222		if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;
223
224		if( loggingVerbosity > 0 )
#	Line 145 | Line 230 | void TreeWriter::Loop() {
230
231		tree->Branch("photon", &photon);
232		tree->Branch("jet", &jet);
233	+	tree->Branch("electron", &electron);
234	+	tree->Branch("muon", &muon);
235		tree->Branch("met", &met, "met/F");
236	+	tree->Branch("metPhi", &met_phi, "metPhi/F");
237	+	tree->Branch("type1met", &type1met, "type1met/F");
238	+	tree->Branch("type1metPhi", &type1met_phi, "type1metPhi/F");
239	+	tree->Branch("ht", &ht, "ht/F");
240		tree->Branch("nVertex", &nVertex, "nVertex/I");
241	<	tree->Branch("weigth", &weight, "weight/F");
151	<	tree->Branch("nElectron", &nElectron, "nElectron/I");
241	>	tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
242
243
244	<	for (long jentry=0; jentry < processNEvents; jentry++) {
244	>	for (long jentry=0; jentry < processNEvents; ++jentry) {
245		if ( loggingVerbosity>0 && jentry%reportEvery==0 )
246		std::cout << jentry << " / " << processNEvents << std :: endl;
247		inputTree->GetEntry(jentry);
248
249		photon.clear();
250		jet.clear();
251	+	electron.clear();
252	+	muon.clear();
253	+	ht = 0;
254	+
255	+	// weights
256	+	if (pileupHisto == 0) {
257	+	pu_weight = 1.;
258	+	} else {
259	+	float trueNumInteractions = -1;
260	+	for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
261	+	iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) {
262	+	if (iBX->BX == 0)
263	+	trueNumInteractions = iBX->trueNumInteractions;
264	+	}
265	+	pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
266	+	}
267	+
268
269		// photons
270	+	if( loggingVerbosity > 1 )
271	+	std::cout << "Process photons" << std::endl;
272		std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
273		for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
274	<	it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
275	<	if( ! it->isEB() && skim )
274	>	it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
275	>	if( !(it->isEB() || it->isEE()) && skim )
276		continue;
277		thisphoton->pt = getPtFromMatchedJet( *it, *event );
278	<	if( thisphoton->pt < 80 && skim )
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->eta = it->momentum.Eta();
301		thisphoton->phi = it->momentum.Phi();
173	–	thisphoton->chargedIso = it->chargedHadronIso;
174	–	thisphoton->neutralIso = it->neutralHadronIso;
175	–	thisphoton->photonIso = it->photonIso;
176	–	if ( it->r9 > 1 && skim ) // if == 1 ?
177	–	continue;
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	+
312		if( photon.size() == 0 && skim )
313		continue;
314		std::sort( photon.begin(), photon.end(), tree::EtGreater);
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 194 | 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++) {
327	>	it != jetColl.end(); ++it) {
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;
#	Line 203 | Line 333 | void TreeWriter::Loop() {
333		float scale = s_it->second;
334		TLorentzVector corrP4 = scale * it->momentum;
335
336	<	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
336	>	if(std::fabs(corrP4.Eta()) > 3.0 && skim ) continue;
337	>	if(corrP4.Et() < 30 && skim ) continue;
338		thisjet->pt = corrP4.Et();
339		thisjet->eta = corrP4.Eta();
340		thisjet->phi = corrP4.Phi();
341	+	thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
342	+	// jet composition
343	+	thisjet->chargedHadronEnergy = it->chargedHadronEnergy;
344	+	thisjet->neutralHadronEnergy = it->neutralHadronEnergy;
345	+	thisjet->photonEnergy = it->photonEnergy;
346	+	thisjet->electronEnergy = it->electronEnergy;
347	+	thisjet->muonEnergy = it->muonEnergy;
348	+	thisjet->HFHadronEnergy = it->HFHadronEnergy;
349	+	thisjet->HFEMEnergy = it->HFEMEnergy;
350	+	thisjet->chargedEmEnergy = it->chargedEmEnergy;
351	+	thisjet->chargedMuEnergy = it->chargedMuEnergy;
352	+	thisjet->neutralEmEnergy = it->neutralEmEnergy;
353	+
354	+	if( loggingVerbosity > 2 )
355	+	std::cout << " p_T, jet = " << thisjet->pt << std::endl;
356	+
357		jet.push_back( *thisjet );
358	+	ht += thisjet->pt;
359		}// for jet
360		}// if, else
361		if( jet.size() < 2 && skim )
362		continue;
363		std::sort( jet.begin(), jet.end(), tree::EtGreater);
364	+	if( loggingVerbosity > 1 )
365	+	std::cout << "Found " << jet.size() << " jets" << std::endl;
366	+
367
368		// met
369		std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
370		susy::MET* metobj = &(met_it->second);
371		met = metobj->met();
372	+	met_phi = metobj->mEt.Phi();
373	+	if( loggingVerbosity > 2 )
374	+	std::cout << " met = " << met << std::endl;
375	+
376	+	std::map<TString, susy::MET>::iterator type1met_it = event->metMap.find("pfType1CorrectedMet");
377	+	susy::MET* type1metobj = &(type1met_it->second);
378	+	type1met = type1metobj->met();
379	+	type1met_phi = type1metobj->mEt.Phi();
380	+	if( loggingVerbosity > 2 )
381	+	std::cout << " type1met = " << type1met << std::endl;
382
383		// electrons
384	<	std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
385	<	nElectron = eVector.size();
384	>	tree::Particle* thiselectron = new tree::Particle();
385	>	map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
386	>	if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
387	>	cout << "gsfElectrons not found!" << endl;
388	>	} else {
389	>	cout << "now begin electrons " << endl;
390	>	for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
391	>	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification for veto electrons
392	>	if( it->momentum.Pt() < 1  || it->momentum.Pt() > 1e6 )
393	>	continue; // spike rejection
394	>	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso
395	>	- effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. )
396	>	) / it->momentum.Pt();
397	>	cout << iso << endl;
398	>	cout << d0correction( *it, *event ) << endl;
399	>	float test = d0correction( *it, *event );
400	>	cout << test << endl;
401	>	float dZ = std::fabs( dZcorrection( *it, *event ) );
402	>	cout <<" find e" << endl;
403	>	float d0 = 0.;
404	>	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
424	>	continue;
425	>	// TODO: conversion rejection information not implemented yet, see twiki for more details
426
427	<	// vertices
427	>	thiselectron->pt = it->momentum.Pt();
428	>	if( thiselectron->pt < 20 )
429	>	continue;
430	>	if( loggingVerbosity > 2 )
431	>	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
432	>	thiselectron->eta = it->momentum.Eta();
433	>	thiselectron->phi = it->momentum.Phi();
434	>	electron.push_back( *thiselectron );
435	>	}
436	>	}
437	>	if( loggingVerbosity > 1 )
438	>	std::cout << "Found " << electron.size() << " electrons" << std::endl;
439
440	+	// muons
441	+	std::vector<susy::Muon> mVector = event->muons["muons"];
442	+	tree::Particle* thismuon = new tree::Particle();
443	+	for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
444	+	if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) )
445	+	continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon
446	+	thismuon->pt = it->momentum.Et();
447	+	if( thismuon->pt < 20 )
448	+	continue;
449	+	thismuon->eta = it->momentum.Eta();
450	+	thismuon->phi = it->momentum.Phi();
451	+	muon.push_back( *thismuon );
452	+	}
453	+	if( loggingVerbosity > 1 )
454	+	std::cout << "Found " << muon.size() << " muons" << std::endl;
455	+
456	+
457	+	// vertices
458		nVertex = event->vertices.size();
229	–	weight = 1;
459
460	<	// bjets
460	>	if( ht < 450 && skim)
461	>	continue;
462
463		tree->Fill();
464		} // for jentry
465
466
237	–
238	–	tree->Write();
467		outFile->cd();
468	+	eventNumbers->Write();
469	+	tree->Write();
470		outFile->Write();
471		outFile->Close();
472		}
473
244	–	int main(int argc, char** argv) {
245	–	TreeWriter *tw = new TreeWriter("../root-files/susyEvents_qcd_1000-inf_part.root", "myTree.root");
246	–	tw->SetProcessNEvents(-1);
247	–	tw->Loop();
248	–	}
249	–

Diff Legend

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