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.1 by kiesel, Wed Mar 20 09:58:55 2013 UTC vs.
Revision 1.8 by kiesel, Wed Apr 10 20:46:33 2013 UTC

#	Line 1 | Line 1
1	<	#include "TFile.h"
2	<	#include "TTree.h"
3	<	#include "TString.h"
4	<	#include "TLorentzVector.h"
1	>	#include<iostream>
2	>	#include<math.h>
3	>	#include<string>
4
5	<	#include "SusyEvent.h"
5	>	#include "TSystem.h"
6
7	+	#include "treeWriter.h"
8
9	<	namespace tree {
10	<	// In this namespace classes for the trees are defined.
9	>	using namespace std;
10
11	<	class Particle {
13	<	public:
14	<	float pt, eta, phi;
15	<	};
16	<
17	<	class Photon : public Particle {
18	<	public:
19	<	float r9, sigmaIetaIeta, hadTowOverEm, pixelseed;
20	<	float chargedIso, neutralIso, photonIso;
21	<	};
22	<
23	<	class Jet : public Particle{
24	<	public:
25	<	float pt, eta;
26	<	};
27	<
28	<	bool EtGreater(const tree::Particle p1, const tree::Particle p2) {
29	<	return p1.pt > p2.pt;
30	<	}
31	<
32	<	} // end namespace definition
33	<
34	<	class TreeWriter {
35	<	public :
36	<	TreeWriter(TString inputName, TString outputName );
37	<	virtual ~TreeWriter();
38	<	virtual void Loop();
39	<
40	<	void SetProcessNEvents(int nEvents) { processNEvents = nEvents; }
41	<	void SetReportEvents(int nEvents) { reportEvery = nEvents; }
42	<
43	<	TFile *inputFile;
44	<	TTree *inputTree;
45	<	susy::Event *event;
46	<
47	<	TFile *outFile;
48	<	TTree *tree;
49	<
50	<	private:
51	<	int processNEvents; // number of events to be processed
52	<	int reportEvery;
53	<	int loggingVerbosity;
54	<
55	<	// variables which will be stored in the tree
56	<	std::vector<tree::Photon> photon;
57	<	std::vector<tree::Jet> jet;
58	<	float met;
59	<	int nVertex;
60	<	float weight;
61	<	};
62	<
63	<
64	<	TreeWriter::TreeWriter(TString inputName, TString outputName) {
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	+	if (outputName == "") {
24	+	// Here the output filename is generated automaticly
25	+	try{
26	+	string fName = (string) inputName;
27	+	// Search for 'susyEvents_', and get the unique characters afterwards.
28	+	// eg /path/susyEvents_123_Myl.root -> susyTree_123_Mly.root
29	+	outputName = "susyTree_" + fName.substr(fName.find("susyEvents_")+11,-1);
30	+	} catch (int e ) {
31	+	outputName = "susyTree.root";
32	+	}
33	+	}
34	+
35	+	if (loggingVerbosity_>0)
36	+	std::cout << "Open file " << outputName << " for writing." << std::endl;
37		// open the output file
38		outFile = new TFile( outputName, "recreate" );
39		tree = new TTree("susyTree","Tree for single photon analysis");
#	Line 75 | Line 41 | TreeWriter::TreeWriter(TString inputName
41		// set default parameter
42		processNEvents = -1;
43		reportEvery = 1000;
44	<	loggingVerbosity = 0;
44	>	loggingVerbosity = loggingVerbosity_;
45	>	skim = true;
46
47		}
48
#	Line 84 | Line 51 | TreeWriter::~TreeWriter() {
51		delete inputTree->GetCurrentFile();
52		}
53
54	+	// useful functions
55	+	float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
56	+	return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
57	+	}
58	+
59	+	// correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012
60	+	float chargedHadronIso_corrected(susy::Photon gamma, float rho) {
61	+	float eta = fabs(gamma.caloPosition.Eta());
62	+	float ea;
63	+
64	+	if(eta < 1.0) ea = 0.012;
65	+	else if(eta < 1.479) ea = 0.010;
66	+	else if(eta < 2.0) ea = 0.014;
67	+	else if(eta < 2.2) ea = 0.012;
68	+	else if(eta < 2.3) ea = 0.016;
69	+	else if(eta < 2.4) ea = 0.020;
70	+	else ea = 0.012;
71	+
72	+	float iso = gamma.chargedHadronIso;
73	+	iso = max(iso - rho*ea, (float)0.);
74	+
75	+	return iso;
76	+	}
77	+
78	+	float neutralHadronIso_corrected(susy::Photon gamma, float rho) {
79	+	float eta = fabs(gamma.caloPosition.Eta());
80	+	float ea;
81	+
82	+	if(eta < 1.0) ea = 0.030;
83	+	else if(eta < 1.479) ea = 0.057;
84	+	else if(eta < 2.0) ea = 0.039;
85	+	else if(eta < 2.2) ea = 0.015;
86	+	else if(eta < 2.3) ea = 0.024;
87	+	else if(eta < 2.4) ea = 0.039;
88	+	else ea = 0.072;
89	+
90	+	float iso = gamma.neutralHadronIso;
91	+	iso = max(iso - rho*ea, (float)0.);
92	+
93	+	return iso;
94	+	}
95	+
96	+	float photonIso_corrected(susy::Photon gamma, float rho) {
97	+	float eta = fabs(gamma.caloPosition.Eta());
98	+	float ea;
99	+
100	+	if(eta < 1.0) ea = 0.148;
101	+	else if(eta < 1.479) ea = 0.130;
102	+	else if(eta < 2.0) ea = 0.112;
103	+	else if(eta < 2.2) ea = 0.216;
104	+	else if(eta < 2.3) ea = 0.262;
105	+	else if(eta < 2.4) ea = 0.260;
106	+	else ea = 0.266;
107	+
108	+	float iso = gamma.photonIso;
109	+	iso = max(iso - rho*ea, (float)0.);
110	+
111	+	return iso;
112	+	}
113	+
114	+
115	+
116	+
117	+	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
118	+	/**
119	+	* \brief Takes jet p_T as photon p_T
120	+	*
121	+	* At first all jets with DeltaR < 0.3 (isolation cone) are searched.
122	+	* If several jets are found, take the one with the minimal pt difference
123	+	* compared to the photon. If no such jets are found, keep the photon_pt
124	+	* TODO: remove photon matched jet from jet-selection?
125	+	*/
126	+	std::vector<susy::PFJet> nearJets;
127	+	nearJets.clear();
128	+
129	+	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
130	+	if(pfJets_it == myEvent.pfJets.end()){
131	+	if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
132	+	} else {
133	+	susy::PFJetCollection& jetColl = pfJets_it->second;
134	+	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
135	+	it != jetColl.end(); it++) {
136	+	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
137	+	if (s_it == it->jecScaleFactors.end()) {
138	+	std::cout << "JEC is not available for this jet!!!" << std::endl;
139	+	continue;
140	+	}
141	+	float scale = s_it->second;
142	+	TLorentzVector corrP4 = scale * it->momentum;
143	+	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
144	+	if (deltaR_ > 0.3) continue;
145	+	if( loggingVerbosity > 0 )
146	+	std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
147	+	nearJets.push_back( *it );
148	+	}// for jet
149	+	}// if, else
150	+
151	+	if ( nearJets.size() == 0 ) {
152	+	//std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
153	+	return myPhoton.momentum.Et();
154	+	}
155	+
156	+	float pt = 0;
157	+	float minPtDifferenz = 1E20; // should be very high
158	+	for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
159	+	it != jetEnd; it++ ) {
160	+	float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
161	+	if (  ptDiff < minPtDifferenz ) {
162	+	minPtDifferenz = ptDiff;
163	+	pt = it->momentum.Et();
164	+	}
165	+	}
166	+
167	+	// testing
168	+	if( nearJets.size() > 1 && loggingVerbosity > 0 )
169	+	std::cout << "There are several jets matching to this photon. "
170	+	<< "Please check if jet-matching is correct." << std::endl;
171	+	return pt;
172	+	}
173	+
174	+
175		void TreeWriter::Loop() {
176	+
177		// here the event loop is implemented and the tree is filled
178		if (inputTree == 0) return;
179
#	Line 102 | Line 191 | void TreeWriter::Loop() {
191		tree->Branch("photon", &photon);
192		tree->Branch("jet", &jet);
193		tree->Branch("met", &met, "met/F");
194	+	tree->Branch("ht", &ht, "ht/F");
195		tree->Branch("nVertex", &nVertex, "nVertex/I");
196	<	tree->Branch("weigth", &weight, "weight/F");
196	>	tree->Branch("nElectron", &nElectron, "nElectron/I");
197
198
199		for (long jentry=0; jentry < processNEvents; jentry++) {
200	<	if (jentry%reportEvery==0)
200	>	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
201		std::cout << jentry << " / " << processNEvents << std :: endl;
202		inputTree->GetEntry(jentry);
203
204		photon.clear();
205		jet.clear();
206	+	ht = 0;
207
208		// photons
209	+	if( loggingVerbosity > 1 )
210	+	std::cout << "Process photons" << std::endl;
211		std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
212		for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
213		it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
214	<	if( ! it->isEB() )
214	>	if( ! it->isEB() && skim )
215		continue;
216	<	thisphoton->pt = it->momentum.Et();
217	<	if( thisphoton->pt < 80 )
216	>	thisphoton->pt = getPtFromMatchedJet( *it, *event );
217	>	if( thisphoton->pt < 80 && skim )
218		continue;
219		thisphoton->eta = it->momentum.Eta();
220	<	thisphoton->chargedIso = it->chargedHadronIso;
221	<	thisphoton->neutralIso = it->neutralHadronIso;
222	<	thisphoton->photonIso = it->photonIso;
223	<	if ( it->r9 > 1 ) // if == 1 ?
220	>	thisphoton->phi = it->momentum.Phi();
221	>	thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25);
222	>	thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25);
223	>	thisphoton->photonIso = photonIso_corrected(*it, event->rho25);
224	>	if ( it->r9 >= 1 && skim )
225		continue;
226		thisphoton->r9 = it->r9;
227		thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
228		thisphoton->hadTowOverEm = it->hadTowOverEm;
229		thisphoton->pixelseed = it->nPixelSeeds;
230		photon.push_back( *thisphoton );
231	+	ht += thisphoton->pt;
232	+	if( loggingVerbosity > 2 )
233	+	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
234		}
235	<	if( photon.size() == 0 )
235	>
236	>	if( photon.size() == 0 && skim )
237		continue;
238		std::sort( photon.begin(), photon.end(), tree::EtGreater);
239	<
239	>	if( loggingVerbosity > 1 )
240	>	std::cout << "Found " << photon.size() << " photons" << std::endl;
241
242		// jets
243		std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
#	Line 158 | Line 257 | void TreeWriter::Loop() {
257		float scale = s_it->second;
258		TLorentzVector corrP4 = scale * it->momentum;
259
260	<	if(std::abs(corrP4.Eta()) > 3.0) continue;
260	>	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
261	>	if(corrP4.Et() < 30 && skim ) continue;
262		thisjet->pt = corrP4.Et();
263		thisjet->eta = corrP4.Eta();
264	+	thisjet->phi = corrP4.Phi();
265		jet.push_back( *thisjet );
266	+	ht += thisjet->pt;
267		}// for jet
268		}// if, else
269	<	if( jet.size() == 0 )
270	<	std::cout << "error, no jets found " << std::endl;
271	<	else
170	<	std::sort( jet.begin(), jet.end(), tree::EtGreater);
269	>	if( jet.size() < 2 && skim )
270	>	continue;
271	>	std::sort( jet.begin(), jet.end(), tree::EtGreater);
272
273		// met
274		std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
275		susy::MET* metobj = &(met_it->second);
276		met = metobj->met();
277
278	+	// electrons
279	+	//std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
280	+	//nElectron = eVector.size();
281	+
282		// vertices
283
284		nVertex = event->vertices.size();
285	<	weight = 1;
285	>
286	>	if( ht < 450 && skim)
287	>	continue;
288
289		tree->Fill();
290		} // for jentry
#	Line 189 | Line 296 | void TreeWriter::Loop() {
296		outFile->Write();
297		outFile->Close();
298		}
299	<	/*
193	<	int main(int argc, char** argv) {
194	<	TreeWriter *tw = new TreeWriter("qcd-1000-nTuple-test.root", "myQCDTree.root");
195	<	tw->SetProcessNEvents(10);
196	<	tw->Loop();
197	<	}
198	<	*/
299	>

Diff Legend

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