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.22 by kiesel, Wed Apr 24 13:24:30 2013 UTC

#	Line 1 | Line 1
1	<	#include<iostream>
2	<	#include<math.h>
1	>	#include "treeWriter.h"
2
3	<	#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	<	};
3	>	using namespace std;
4
5	<	TreeWriter::TreeWriter(TString inputName, TString outputName) {
5	>	TreeWriter::TreeWriter( std::string inputName, std::string outputName, int loggingVerbosity_ ) {
6		// read the input file
7	<	inputFile = new TFile( inputName, "read" );
8	<	inputTree = (TTree*) inputFile->Get("susyTree");
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	>	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	<	outFile = new TFile( outputName, "recreate" );
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		// set default parameter
37		processNEvents = -1;
38		reportEvery = 1000;
39	<	loggingVerbosity = 0;
39	>	loggingVerbosity = loggingVerbosity_;
40	>	skim = true;
41	>	pileupHisto = 0;
42	>	}
43
44	+	void TreeWriter::PileUpWeightFile( string pileupFileName ) {
45	+	TFile *puFile = new TFile( pileupFileName.c_str() );
46	+	pileupHisto = (TH1F*) puFile->Get("pileup");
47		}
48
49		TreeWriter::~TreeWriter() {
50	<	if (!inputTree) return;
51	<	delete inputTree->GetCurrentFile();
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	<	float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
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 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 TreeWriter::getPtFromMatchedJet( susy::Photon photon, susy::Event event ) {
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 );
84	>	float ea;
85	>	if( eta < 1.0 ) ea = 0.13;
86	>	else if( eta < 1.479 ) ea = 0.14;
87	>	else if( eta < 2.0 ) ea = 0.07;
88	>	else if( eta < 2.2 ) ea = 0.09;
89	>	else if( eta < 2.3 ) ea = 0.11;
90	>	else if( eta < 2.4 ) ea = 0.11;
91	>	else ea = 0.14;
92	>	return ea;
93	>	}
94	>
95	>	// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
96	>	float chargedHadronIso_corrected(const susy::Photon& gamma, float rho) {
97	>	float eta = fabs(gamma.caloPosition.Eta());
98	>	float ea;
99	>
100	>	if(eta < 1.0) ea = 0.012;
101	>	else if(eta < 1.479) ea = 0.010;
102	>	else if(eta < 2.0) ea = 0.014;
103	>	else if(eta < 2.2) ea = 0.012;
104	>	else if(eta < 2.3) ea = 0.016;
105	>	else if(eta < 2.4) ea = 0.020;
106	>	else ea = 0.012;
107	>
108	>	float iso = gamma.chargedHadronIso;
109	>	iso = max(iso - rho*ea, (float)0.);
110	>
111	>	return iso;
112	>	}
113	>
114	>	float neutralHadronIso_corrected(const susy::Photon& gamma, float rho) {
115	>	float eta = fabs(gamma.caloPosition.Eta());
116	>	float ea;
117	>
118	>	if(eta < 1.0) ea = 0.030;
119	>	else if(eta < 1.479) ea = 0.057;
120	>	else if(eta < 2.0) ea = 0.039;
121	>	else if(eta < 2.2) ea = 0.015;
122	>	else if(eta < 2.3) ea = 0.024;
123	>	else if(eta < 2.4) ea = 0.039;
124	>	else ea = 0.072;
125	>
126	>	float iso = gamma.neutralHadronIso;
127	>	iso = max(iso - rho*ea, (float)0.);
128	>
129	>	return iso;
130	>	}
131	>
132	>	float photonIso_corrected(const susy::Photon& gamma, float rho) {
133	>	float eta = fabs(gamma.caloPosition.Eta());
134	>	float ea;
135	>
136	>	if(eta < 1.0) ea = 0.148;
137	>	else if(eta < 1.479) ea = 0.130;
138	>	else if(eta < 2.0) ea = 0.112;
139	>	else if(eta < 2.2) ea = 0.216;
140	>	else if(eta < 2.3) ea = 0.262;
141	>	else if(eta < 2.4) ea = 0.260;
142	>	else ea = 0.266;
143	>
144	>	float iso = gamma.photonIso;
145	>	iso = max(iso - rho*ea, (float)0.);
146	>
147	>	return iso;
148	>	}
149	>
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		*
172		* At first all jets with DeltaR < 0.3 (isolation cone) are searched.
173		* If several jets are found, take the one with the minimal pt difference
174	<	* compared to the photon. If no such jets are found, return 0
174	>	* compared to the photon. If no such jets are found, keep the photon_pt
175		* TODO: remove photon matched jet from jet-selection?
176		*/
177		std::vector<susy::PFJet> nearJets;
178		nearJets.clear();
179
180	<	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event.pfJets.find("ak5");
181	<	if(pfJets_it == event.pfJets.end()){
182	<	if(event.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(photon.momentum, corrP4 );
195	<	if (deltaR_ > 0.3) continue;
196	<	nearJets.push_back( *it );
101	<	}// for jet
102	<	}// 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	<	std::cout << "No jet with ΔR < .3 found, set p_T to 0" << std::endl;
200	<	return 0;
199	>	if( loggingVerbosity > 1 )
200	>	std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
201	>	return myPhoton.momentum.Et();
202		}
203
204		float pt = 0;
205		float minPtDifferenz = 1E20; // should be very high
206		for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
207	<	it != jetEnd; it++ ) {
208	<	float ptDiff = fabs(photon.momentum.Et() - it->momentum.Et());
207	>	it != jetEnd; ++it ) {
208	>	float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
209		if (  ptDiff < minPtDifferenz ) {
210		minPtDifferenz = ptDiff;
211		pt = it->momentum.Et();
#	Line 118 | Line 213 | float TreeWriter::getPtFromMatchedJet( s
213		}
214
215		// testing
216	<	if( nearJets.size() > 1 )
216	>	if( nearJets.size() > 1 && loggingVerbosity > 0 )
217		std::cout << "There are several jets matching to this photon. "
218	<	<< "Please check if the upper assumtion is reasonable" << std::endl;
218	>	<< "Please check if jet-matching is correct." << std::endl;
219		return pt;
220		}
221
222
223		void TreeWriter::Loop() {
224	<	// only for testing
225	<	bool skim = true;
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
143	–	tree::Photon *thisphoton = new tree::Photon();
144	–	tree::Jet *thisjet = new tree::Jet();
145	–
245		tree->Branch("photon", &photon);
246		tree->Branch("jet", &jet);
247	+	tree->Branch("electron", &electron);
248	+	tree->Branch("muon", &muon);
249		tree->Branch("met", &met, "met/F");
250	+	tree->Branch("metPhi", &met_phi, "metPhi/F");
251	+	tree->Branch("type1met", &type1met, "type1met/F");
252	+	tree->Branch("type1metPhi", &type1met_phi, "type1metPhi/F");
253	+	tree->Branch("ht", &ht, "ht/F");
254		tree->Branch("nVertex", &nVertex, "nVertex/I");
255	<	tree->Branch("weigth", &weight, "weight/F");
256	<	tree->Branch("nElectron", &nElectron, "nElectron/I");
257	<
255	>	tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
256	>	tree->Branch("genElectron", &genElectron);
257	>	tree->Branch("genPhoton", &genPhoton);
258
259	<	for (long jentry=0; jentry < processNEvents; jentry++) {
260	<	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
259	>	for (unsigned long jentry=0; jentry < processNEvents; ++jentry) {
260	>	if ( loggingVerbosity>1 || jentry%reportEvery==0 )
261		std::cout << jentry << " / " << processNEvents << std :: endl;
262	<	inputTree->GetEntry(jentry);
262	>	inputTree->LoadTree( jentry );
263	>	inputTree->GetEntry( jentry );
264
265		photon.clear();
266		jet.clear();
267	+	electron.clear();
268	+	muon.clear();
269	+	genElectron.clear();
270	+	genPhoton.clear();
271	+	ht = 0;
272	+
273	+	// weights
274	+	if (pileupHisto == 0) {
275	+	pu_weight = 1.;
276	+	} else {
277	+	float trueNumInteractions = -1;
278	+	for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
279	+	iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) {
280	+	if (iBX->BX == 0)
281	+	trueNumInteractions = iBX->trueNumInteractions;
282	+	}
283	+	pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
284	+	}
285	+
286	+	// get ak5 jets
287	+	std::vector<susy::PFJet> jetVector = event->pfJets["ak5"];
288
289		// photons
290	<	std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
291	<	for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
292	<	it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
293	<	if( ! it->isEB() && skim )
294	<	continue;
168	<	thisphoton->pt = getPtFromMatchedJet( *it, *event );
169	<	if( thisphoton->pt < 80 && skim )
290	>	std::vector<susy::Photon> photonVector = event->photons["photons"];
291	>
292	>	for(std::vector<susy::Photon>::iterator it = photonVector.begin();
293	>	it != photonVector.end(); ++it ) {
294	>	if( !(it->isEE() || it->isEB()) && it->isEBEtaGap() && it->isEBPhiGap() && it->isEERingGap() && it->isEEDeeGap() && it->isEBEEGap() && skim )
295		continue;
296	<	thisphoton->eta = it->momentum.Eta();
297	<	thisphoton->phi = it->momentum.Phi();
298	<	thisphoton->chargedIso = it->chargedHadronIso;
299	<	thisphoton->neutralIso = it->neutralHadronIso;
300	<	thisphoton->photonIso = it->photonIso;
301	<	if ( it->r9 > 1 && skim ) // if == 1 ?
296	>	tree::Photon thisphoton;
297	>	thisphoton.pt = getPtFromMatchedJet( *it, jetVector, loggingVerbosity );
298	>
299	>	thisphoton.chargedIso = chargedHadronIso_corrected(*it, event->rho25);
300	>	thisphoton.neutralIso = neutralHadronIso_corrected(*it, event->rho25);
301	>	thisphoton.photonIso = photonIso_corrected(*it, event->rho25);
302	>
303	>	bool loose_photon_barrel = thisphoton.pt>20
304	>	&& it->isEB()
305	>	&& it->passelectronveto
306	>	&& it->hadTowOverEm<0.05
307	>	&& it->sigmaIetaIeta<0.012
308	>	&& thisphoton.chargedIso<2.6
309	>	&& thisphoton.neutralIso<3.5+0.04*thisphoton.pt
310	>	&& thisphoton.photonIso<1.3+0.005*thisphoton.pt;
311	>	bool loose_photon_endcap = thisphoton.pt > 20
312	>	&& it->isEE()
313	>	&& it->passelectronveto
314	>	&& it->hadTowOverEm<0.05
315	>	&& it->sigmaIetaIeta<0.034
316	>	&& thisphoton.chargedIso<2.3
317	>	&& thisphoton.neutralIso<2.9+0.04*thisphoton.pt;
318	>
319	>	if(!(loose_photon_endcap || loose_photon_barrel || thisphoton.pt > 75 ) && skim )
320		continue;
321	<	thisphoton->r9 = it->r9;
322	<	thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
323	<	thisphoton->hadTowOverEm = it->hadTowOverEm;
324	<	thisphoton->pixelseed = it->nPixelSeeds;
325	<	photon.push_back( *thisphoton );
321	>	thisphoton.eta = it->momentum.Eta();
322	>	thisphoton.phi = it->momentum.Phi();
323	>	thisphoton.r9 = it->r9;
324	>	thisphoton.sigmaIetaIeta = it->sigmaIetaIeta;
325	>	thisphoton.hadTowOverEm = it->hadTowOverEm;
326	>	thisphoton.pixelseed = it->nPixelSeeds;
327	>	thisphoton.conversionSafeVeto = it->passelectronveto;
328	>	photon.push_back( thisphoton );
329	>	if( loggingVerbosity > 2 )
330	>	std::cout << " p_T, gamma = " << thisphoton.pt << std::endl;
331		}
332	+
333		if( photon.size() == 0 && skim )
334		continue;
335		std::sort( photon.begin(), photon.end(), tree::EtGreater);
336	+	if( loggingVerbosity > 1 )
337	+	std::cout << "Found " << photon.size() << " photons" << std::endl;
338
339		// jets
189	–	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
190	–	if(pfJets_it == event->pfJets.end()){
191	–	if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
192	–	} else {
340
194	–	susy::PFJetCollection& jetColl = pfJets_it->second;
341
342	<	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
343	<	it != jetColl.end(); it++) {
344	<	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
345	<	if (s_it == it->jecScaleFactors.end()) {
346	<	std::cout << "JEC is not available for this jet!!!" << std::endl;
347	<	continue;
348	<	}
349	<	float scale = s_it->second;
350	<	TLorentzVector corrP4 = scale * it->momentum;
342	>	for(std::vector<susy::PFJet>::iterator it = jetVector.begin();
343	>	it != jetVector.end(); ++it) {
344	>	tree::Jet thisjet;
345	>
346	>	// scale with JEC
347	>	float scale = 1.;
348	>	if(it->jecScaleFactors.count("L2L3") == 0)
349	>	std::cout << "ERROR: JEC is not available for this jet" << std::endl;
350	>	else
351	>	scale = it->jecScaleFactors.find("L2L3")->second;
352	>	TLorentzVector corrP4 = scale * it->momentum;
353	>
354	>	// Calculate HT.
355	>	// The definiton differs from the saved jet, since trigger is described better
356	>	if( std::abs( corrP4.Eta() ) < 3 && corrP4.Pt() > 40 )
357	>	ht += thisjet.pt;
358	>
359	>	if(std::abs(corrP4.Eta()) > 2.6 && skim ) continue;
360	>	if(corrP4.Pt() < 30 && skim ) continue;
361	>	thisjet.pt = corrP4.Pt();
362	>	thisjet.eta = corrP4.Eta();
363	>	thisjet.phi = corrP4.Phi();
364	>	thisjet.bCSV = it->bTagDiscriminators[susy::kCSV];
365	>	// jet composition
366	>	thisjet.chargedHadronEnergy = it->chargedHadronEnergy;
367	>	thisjet.neutralHadronEnergy = it->neutralHadronEnergy;
368	>	thisjet.photonEnergy = it->photonEnergy;
369	>	thisjet.electronEnergy = it->electronEnergy;
370	>	thisjet.muonEnergy = it->muonEnergy;
371	>	thisjet.HFHadronEnergy = it->HFHadronEnergy;
372	>	thisjet.HFEMEnergy = it->HFEMEnergy;
373	>	thisjet.chargedEmEnergy = it->chargedEmEnergy;
374	>	thisjet.chargedMuEnergy = it->chargedMuEnergy;
375	>	thisjet.neutralEmEnergy = it->neutralEmEnergy;
376	>
377	>	if( loggingVerbosity > 2 )
378	>	std::cout << " p_T, jet = " << thisjet.pt << std::endl;
379
380	<	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
381	<	thisjet->pt = corrP4.Et();
208	<	thisjet->eta = corrP4.Eta();
209	<	thisjet->phi = corrP4.Phi();
210	<	jet.push_back( *thisjet );
211	<	}// for jet
212	<	}// if, else
380	>	jet.push_back( thisjet );
381	>	}// for jet
382		if( jet.size() < 2 && skim )
383		continue;
384		std::sort( jet.begin(), jet.end(), tree::EtGreater);
385	+	if( loggingVerbosity > 1 )
386	+	std::cout << "Found " << jet.size() << " jets" << std::endl;
387	+
388	+	if( ht < 450 && skim)
389	+	continue;
390	+
391	+
392
393		// met
394		std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
395		susy::MET* metobj = &(met_it->second);
396		met = metobj->met();
397	+	met_phi = metobj->mEt.Phi();
398	+	if( loggingVerbosity > 2 )
399	+	std::cout << " met = " << met << std::endl;
400	+
401	+	std::map<TString, susy::MET>::iterator type1met_it = event->metMap.find("pfType1CorrectedMet");
402	+	susy::MET* type1metobj = &(type1met_it->second);
403	+	type1met = type1metobj->met();
404	+	type1met_phi = type1metobj->mEt.Phi();
405	+	if( loggingVerbosity > 2 )
406	+	std::cout << " type1met = " << type1met << std::endl;
407
408		// electrons
409		std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
410	<	nElectron = eVector.size();
410	>	for(std::vector<susy::Electron>::iterator it = eVector.begin(); it < eVector.end(); ++it) {
411	>	tree::Particle thiselectron;
412	>	if( loggingVerbosity > 2 )
413	>	cout << " electron pt = " << it->momentum.Pt() << endl;
414	>	// for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification
415	>	// use veto electrons
416	>	if( it->momentum.Pt() < 20  || it->momentum.Pt() > 1e6 )
417	>	continue; // spike rejection
418	>	float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (float)0. )
419	>	) / it->momentum.Pt();
420	>	float d0 = d0correction( *it, *event );
421	>	float dZ = std::abs( dZcorrection( *it, *event ) );
422	>	if ( it->isEB() ){
423	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007
424	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8
425	>	|| it->sigmaIetaIeta > 0.01
426	>	|| it->hcalOverEcalBc > 0.15
427	>	|| d0 > 0.04
428	>	|| dZ > 0.2
429	>	|| iso > 0.15 )
430	>	continue;
431	>	}
432	>	else if( it->isEE() ) {
433	>	if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01
434	>	|| fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7
435	>	|| it->sigmaIetaIeta > 0.03
436	>	|| d0 > 0.04
437	>	|| dZ > 0.2
438	>	|| iso > 0.15 )
439	>	continue;
440	>	}
441	>	else // not in barrel nor in endcap
442	>	continue;
443
444	<	// vertices
444	>	thiselectron.pt = it->momentum.Pt();
445	>	if( loggingVerbosity > 2 )
446	>	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
447	>	thiselectron.eta = it->momentum.Eta();
448	>	thiselectron.phi = it->momentum.Phi();
449	>	electron.push_back( thiselectron );
450	>	}
451	>	if( loggingVerbosity > 1 )
452	>	std::cout << "Found " << electron.size() << " electrons" << std::endl;
453	>
454	>	// muons
455	>	tree::Particle thismuon;
456	>	std::vector<susy::Muon> mVector = event->muons["muons"];
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
470	+	// vertices
471		nVertex = event->vertices.size();
229	–	weight = 1;
472
473	<	// bjets
473	>	tree::Particle thisGenParticle;
474	>	for( std::vector<susy::Particle>::iterator it = event->genParticles.begin(); it != event->genParticles.end(); ++it ) {
475	>	if( it->status == 3 ) { // hard interaction
476	>	switch( std::abs(it->pdgId) ) {
477	>	thisGenParticle.pt = it->momentum.Pt();
478	>	thisGenParticle.eta = it->momentum.Eta();
479	>	thisGenParticle.phi = it->momentum.Phi();
480	>	case 22: // photon
481	>	if( thisGenParticle.pt > 75 )
482	>	genPhoton.push_back( thisGenParticle );
483	>	case 11: // electron
484	>	if( thisGenParticle.pt > 20 ) // pt cut is lower to estimate fake rate in all pt bins
485	>	genElectron.push_back( thisGenParticle );
486	>	}
487	>	}
488	>	}
489
490		tree->Fill();
491		} // for jentry
492
493
237	–
238	–	tree->Write();
494		outFile->cd();
495	+	eventNumbers->Write();
496	+	tree->Write();
497		outFile->Write();
498		outFile->Close();
499		}
500
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