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.3 by kiesel, Thu Mar 21 09:01:14 2013 UTC vs.
Revision 1.9 by kiesel, Sat Apr 13 09:53:02 2013 UTC

#	Line 1 | Line 1
1		#include<iostream>
2	+	#include<math.h>
3	+	#include<string>
4
5	<	#include "TFile.h"
4	<	#include "TTree.h"
5	<	#include "TString.h"
6	<
7	<	#include "SusyEvent.h"
8	<	#include "TreeObjects.h"
9	<
10	<
11	<	class TreeWriter {
12	<	public :
13	<	TreeWriter(TString inputName, TString outputName );
14	<	virtual ~TreeWriter();
15	<	virtual void Loop();
16	<
17	<	void SetProcessNEvents(int nEvents) { processNEvents = nEvents; }
18	<	void SetReportEvents(int nEvents) { reportEvery = nEvents; }
19	<
20	<	TFile *inputFile;
21	<	TTree *inputTree;
22	<	susy::Event *event;
23	<
24	<	TFile *outFile;
25	<	TTree *tree;
26	<
27	<	private:
28	<	int processNEvents; // number of events to be processed
29	<	int reportEvery;
30	<	int loggingVerbosity;
31	<
32	<	// variables which will be stored in the tree
33	<	std::vector<tree::Photon> photon;
34	<	std::vector<tree::Jet> jet;
35	<	float met;
36	<	int nVertex;
37	<	float weight;
38	<	float leadingPhotonPt;
39	<	};
5	>	#include "TSystem.h"
6
7	+	#include "treeWriter.h"
8
9	<	TreeWriter::TreeWriter(TString inputName, TString outputName) {
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	>
18	>	if (loggingVerbosity_ > 0)
19	>	std::cout << "Set Branch Address of susy::Event" << std::endl;
20		event = new susy::Event;
21		inputTree->SetBranchAddress("susyEvent", &event);
22
23		// open the output file
24	+	if (loggingVerbosity_>0)
25	+	std::cout << "Open file " << outputName << " for writing." << std::endl;
26		outFile = new TFile( outputName, "recreate" );
27		tree = new TTree("susyTree","Tree for single photon analysis");
28
29		// set default parameter
30		processNEvents = -1;
31		reportEvery = 1000;
32	<	loggingVerbosity = 0;
32	>	loggingVerbosity = loggingVerbosity_;
33	>	skim = true;
34	>	pileupHisto = 0;
35	>	}
36
37	+	void TreeWriter::PileUpWeightFile( string pileupFileName ) {
38	+	TFile *puFile = new TFile( pileupFileName.c_str() );
39	+	pileupHisto = (TH1F*) puFile->Get("pileup");
40		}
41
42	+
43	+
44		TreeWriter::~TreeWriter() {
45	<	if (!inputTree) return;
46	<	delete inputTree->GetCurrentFile();
45	>	if (pileupHisto != 0 )
46	>	delete pileupHisto;
47	>	inputTree->GetCurrentFile()->Close();
48	>	}
49	>
50	>	// useful functions
51	>	float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
52	>	return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
53	>	}
54	>
55	>	// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
56	>	float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
57	>	float eta = fabs(gamma.caloPosition.Eta());
58	>	float ea;
59	>
60	>	if(eta < 1.0) ea = 0.012;
61	>	else if(eta < 1.479) ea = 0.010;
62	>	else if(eta < 2.0) ea = 0.014;
63	>	else if(eta < 2.2) ea = 0.012;
64	>	else if(eta < 2.3) ea = 0.016;
65	>	else if(eta < 2.4) ea = 0.020;
66	>	else ea = 0.012;
67	>
68	>	float iso = gamma.chargedHadronIso;
69	>	iso = max(iso - rho*ea, (float)0.);
70	>
71	>	return iso;
72	>	}
73	>
74	>	float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
75	>	float eta = fabs(gamma.caloPosition.Eta());
76	>	float ea;
77	>
78	>	if(eta < 1.0) ea = 0.030;
79	>	else if(eta < 1.479) ea = 0.057;
80	>	else if(eta < 2.0) ea = 0.039;
81	>	else if(eta < 2.2) ea = 0.015;
82	>	else if(eta < 2.3) ea = 0.024;
83	>	else if(eta < 2.4) ea = 0.039;
84	>	else ea = 0.072;
85	>
86	>	float iso = gamma.neutralHadronIso;
87	>	iso = max(iso - rho*ea, (float)0.);
88	>
89	>	return iso;
90		}
91
92	+	float photonIso_corrected(susy::Photon gamma, float rho) {
93	+	float eta = fabs(gamma.caloPosition.Eta());
94	+	float ea;
95	+
96	+	if(eta < 1.0) ea = 0.148;
97	+	else if(eta < 1.479) ea = 0.130;
98	+	else if(eta < 2.0) ea = 0.112;
99	+	else if(eta < 2.2) ea = 0.216;
100	+	else if(eta < 2.3) ea = 0.262;
101	+	else if(eta < 2.4) ea = 0.260;
102	+	else ea = 0.266;
103	+
104	+	float iso = gamma.photonIso;
105	+	iso = max(iso - rho*ea, (float)0.);
106	+
107	+	return iso;
108	+	}
109	+
110	+	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
111	+	/**
112	+	* \brief Takes jet p_T as photon p_T
113	+	*
114	+	* At first all jets with DeltaR < 0.3 (isolation cone) are searched.
115	+	* If several jets are found, take the one with the minimal pt difference
116	+	* compared to the photon. If no such jets are found, keep the photon_pt
117	+	* TODO: remove photon matched jet from jet-selection?
118	+	*/
119	+	std::vector<susy::PFJet> nearJets;
120	+	nearJets.clear();
121	+
122	+	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
123	+	if(pfJets_it == myEvent.pfJets.end()){
124	+	if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
125	+	} else {
126	+	susy::PFJetCollection& jetColl = pfJets_it->second;
127	+	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
128	+	it != jetColl.end(); ++it) {
129	+	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
130	+	if (s_it == it->jecScaleFactors.end()) {
131	+	std::cout << "JEC is not available for this jet!!!" << std::endl;
132	+	continue;
133	+	}
134	+	float scale = s_it->second;
135	+	TLorentzVector corrP4 = scale * it->momentum;
136	+	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
137	+	if (deltaR_ > 0.3) continue;
138	+	if( loggingVerbosity > 0 )
139	+	std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
140	+	nearJets.push_back( *it );
141	+	}// for jet
142	+	}// if, else
143	+
144	+	if ( nearJets.size() == 0 ) {
145	+	//std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
146	+	return myPhoton.momentum.Et();
147	+	}
148	+
149	+	float pt = 0;
150	+	float minPtDifferenz = 1E20; // should be very high
151	+	for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
152	+	it != jetEnd; ++it ) {
153	+	float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
154	+	if (  ptDiff < minPtDifferenz ) {
155	+	minPtDifferenz = ptDiff;
156	+	pt = it->momentum.Et();
157	+	}
158	+	}
159	+
160	+	// testing
161	+	if( nearJets.size() > 1 && loggingVerbosity > 0 )
162	+	std::cout << "There are several jets matching to this photon. "
163	+	<< "Please check if jet-matching is correct." << std::endl;
164	+	return pt;
165	+	}
166	+
167	+
168		void TreeWriter::Loop() {
169	+
170		// here the event loop is implemented and the tree is filled
171		if (inputTree == 0) return;
172
#	Line 79 | Line 183 | void TreeWriter::Loop() {
183
184		tree->Branch("photon", &photon);
185		tree->Branch("jet", &jet);
186	+	tree->Branch("electron", &electron);
187	+	tree->Branch("muon", &muon);
188		tree->Branch("met", &met, "met/F");
189	+	tree->Branch("ht", &ht, "ht/F");
190		tree->Branch("nVertex", &nVertex, "nVertex/I");
191	<	tree->Branch("weigth", &weight, "weight/F");
85	<	tree->Branch("photonPt", &leadingPhotonPt, "photonPt/F");
191	>	tree->Branch("pu_weight", &pu_weight, "pu_weight/F");
192
193
194	<	for (long jentry=0; jentry < processNEvents; jentry++) {
195	<	if (jentry%reportEvery==0)
194	>	for (long jentry=0; jentry < processNEvents; ++jentry) {
195	>	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
196		std::cout << jentry << " / " << processNEvents << std :: endl;
197		inputTree->GetEntry(jentry);
198
199		photon.clear();
200		jet.clear();
201	+	electron.clear();
202	+	muon.clear();
203	+	ht = 0;
204	+
205	+	// weights
206	+	if (pileupHisto == 0) {
207	+	pu_weight = 1.;
208	+	} else {
209	+	float trueNumInteractions = -1;
210	+	for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin();
211	+	iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) {
212	+	if (iBX->BX == 0)
213	+	trueNumInteractions = iBX->trueNumInteractions;
214	+	}
215	+	pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) );
216	+	}
217	+
218
219		// photons
220	+	if( loggingVerbosity > 1 )
221	+	std::cout << "Process photons" << std::endl;
222		std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
223		for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
224	<	it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
225	<	if( ! it->isEB() )
224	>	it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) {
225	>	if( ! it->isEB() && skim )
226		continue;
227	<	thisphoton->pt = it->momentum.Et();
228	<	if( thisphoton->pt < 80 )
227	>	thisphoton->pt = getPtFromMatchedJet( *it, *event );
228	>	if( thisphoton->pt < 80 && skim )
229		continue;
230		thisphoton->eta = it->momentum.Eta();
231	<	thisphoton->chargedIso = it->chargedHadronIso;
232	<	thisphoton->neutralIso = it->neutralHadronIso;
233	<	thisphoton->photonIso = it->photonIso;
234	<	if ( it->r9 > 1 ) // if == 1 ?
231	>	thisphoton->phi = it->momentum.Phi();
232	>	thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
233	>	thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
234	>	thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
235	>	if ( it->r9 >= 1 && skim )
236		continue;
237		thisphoton->r9 = it->r9;
238		thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
239		thisphoton->hadTowOverEm = it->hadTowOverEm;
240		thisphoton->pixelseed = it->nPixelSeeds;
241		photon.push_back( *thisphoton );
242	+	if( loggingVerbosity > 2 )
243	+	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
244		}
245	<	if( photon.size() == 0 )
245	>
246	>	if( photon.size() == 0 && skim )
247		continue;
248		std::sort( photon.begin(), photon.end(), tree::EtGreater);
249	<	//leadingPhotonPt = photon.at(0)->pt;
250	<
249	>	if( loggingVerbosity > 1 )
250	>	std::cout << "Found " << photon.size() << " photons" << std::endl;
251
252		// jets
253		std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
#	Line 129 | Line 258 | void TreeWriter::Loop() {
258		susy::PFJetCollection& jetColl = pfJets_it->second;
259
260		for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
261	<	it != jetColl.end(); it++) {
261	>	it != jetColl.end(); ++it) {
262		std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
263		if (s_it == it->jecScaleFactors.end()) {
264		std::cout << "JEC is not available for this jet!!!" << std::endl;
#	Line 138 | Line 267 | void TreeWriter::Loop() {
267		float scale = s_it->second;
268		TLorentzVector corrP4 = scale * it->momentum;
269
270	<	if(std::abs(corrP4.Eta()) > 3.0) continue;
270	>	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
271	>	if(corrP4.Et() < 30 && skim ) continue;
272		thisjet->pt = corrP4.Et();
273		thisjet->eta = corrP4.Eta();
274	+	thisjet->phi = corrP4.Phi();
275	+	thisjet->bCSV = it->bTagDiscriminators[susy::kCSV];
276	+	if( loggingVerbosity > 2 )
277	+	std::cout << " p_T, jet = " << thisjet->pt << std::endl;
278	+
279		jet.push_back( *thisjet );
280	+	ht += thisjet->pt;
281		}// for jet
282		}// if, else
283	<	if( jet.size() == 0 )
284	<	std::cout << "error, no jets found " << std::endl;
285	<	else
286	<	std::sort( jet.begin(), jet.end(), tree::EtGreater);
283	>	if( jet.size() < 2 && skim )
284	>	continue;
285	>	std::sort( jet.begin(), jet.end(), tree::EtGreater);
286	>	if( loggingVerbosity > 1 )
287	>	std::cout << "Found " << jet.size() << " jets" << std::endl;
288	>
289
290		// met
291		std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
292		susy::MET* metobj = &(met_it->second);
293		met = metobj->met();
294	+	met_phi = metobj->mEt.Phi();
295	+	if( loggingVerbosity > 2 )
296	+	std::cout << " met = " << met << std::endl;
297	+
298	+	std::map<TString, susy::MET>::iterator type1met_it = event->metMap.find("pfType1CorrectedMet");
299	+	susy::MET* type1metobj = &(type1met_it->second);
300	+	type1met = type1metobj->met();
301	+	type1met_phi = type1metobj->mEt.Phi();
302	+	if( loggingVerbosity > 2 )
303	+	std::cout << " type1met = " << type1met << std::endl;
304	+
305	+	// electrons
306	+	tree::Particle* thiselectron = new tree::Particle();
307	+
308	+	map<TString, vector<susy::Electron> >::iterator eleMap = event->electrons.find("gsfElectrons");
309	+	if(eleMap == event->electrons.end() && loggingVerbosity > 0) {
310	+	cout << "gsfElectrons not found!" << endl;
311	+	} else {
312	+	for(vector<susy::Electron>::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) {
313	+	thiselectron->pt = it->momentum.Et();
314	+	if( thiselectron->pt < 20 )
315	+	continue;
316	+	if( loggingVerbosity > 2 )
317	+	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
318	+	thiselectron->eta = it->momentum.Eta();
319	+	thiselectron->phi = it->momentum.Phi();
320	+	electron.push_back( *thiselectron );
321	+	}
322	+	}
323	+	if( loggingVerbosity > 1 )
324	+	std::cout << "Found " << electron.size() << " electrons" << std::endl;
325
326	<	// vertices
326	>	// this seems not to work yet, where is the bug?
327	>	/*
328	>	std::vector<susy::Electron> eVector = event->electrons["gsfElectronsx"];
329	>	for( std::vector<susy::Electron>::iterator it = eVector.begin(); it != eVector.end(); ++it) {
330	>	thiselectron->pt = it->momentum.Et();
331	>	if( thiselectron->pt < 20 )
332	>	continue;
333	>	if( loggingVerbosity > 2 )
334	>	std::cout << " p_T, electron = " << it->momentum.Et() << std::endl;
335	>	thiselectron->eta = it->momentum.Eta();
336	>	thiselectron->phi = it->momentum.Phi();
337	>	electron.push_back( *thiselectron );
338	>	}
339	>	*/
340	>
341	>	// muons
342	>	std::vector<susy::Muon> mVector = event->muons["muons"];
343	>	tree::Particle* thismuon = new tree::Particle();
344	>	for( std::vector<susy::Muon>::iterator it = mVector.begin(); it != mVector.end(); ++it) {
345	>	thismuon->pt = it->momentum.Et();
346	>	if( thismuon->pt < 20 )
347	>	continue;
348	>	thismuon->eta = it->momentum.Eta();
349	>	thismuon->phi = it->momentum.Phi();
350	>	muon.push_back( *thismuon );
351	>	}
352	>	if( loggingVerbosity > 1 )
353	>	std::cout << "Found " << muon.size() << " muons" << std::endl;
354
355	+
356	+	// vertices
357		nVertex = event->vertices.size();
358	<	weight = 1;
358	>
359	>	if( ht < 450 && skim)
360	>	continue;
361
362		tree->Fill();
363		} // for jentry
364
365
166	–
167	–	tree->Write();
366		outFile->cd();
367	+	tree->Write();
368		outFile->Write();
369		outFile->Close();
370		}
371
173	–	int main(int argc, char** argv) {
174	–	TreeWriter *tw = new TreeWriter("qcd-1000-nTuple-test.root", "myQCDTree.root");
175	–	tw->SetProcessNEvents(10);
176	–	tw->Loop();
177	–	}
178	–

Diff Legend

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