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.7 by kiesel, Tue Apr 9 12:13:55 2013 UTC

#	Line 1 | Line 1
1	<	#include "TFile.h"
2	<	#include "TTree.h"
3	<	#include "TString.h"
4	<	#include "TLorentzVector.h"
5	<
6	<	#include "SusyEvent.h"
7	<
8	<
9	<	namespace tree {
10	<	// In this namespace classes for the trees are defined.
11	<
12	<	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	<	};
1	>	#include<iostream>
2	>	#include<math.h>
3	>	#include<string>
4
5	<	bool EtGreater(const tree::Particle p1, const tree::Particle p2) {
29	<	return p1.pt > p2.pt;
30	<	}
31	<
32	<	} // end namespace definition
5	>	#include "TSystem.h"
6
7	<	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	<	};
7	>	#include "treeWriter.h"
8
9	+	using namespace std;
10
11	<	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	+	float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
55	+	return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
56	+	}
57	+
58	+	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
59	+	/**
60	+	* \brief Takes jet p_T as photon p_T
61	+	*
62	+	* At first all jets with DeltaR < 0.3 (isolation cone) are searched.
63	+	* If several jets are found, take the one with the minimal pt difference
64	+	* compared to the photon. If no such jets are found, keep the photon_pt
65	+	* TODO: remove photon matched jet from jet-selection?
66	+	*/
67	+	std::vector<susy::PFJet> nearJets;
68	+	nearJets.clear();
69	+
70	+	std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
71	+	if(pfJets_it == myEvent.pfJets.end()){
72	+	if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
73	+	} else {
74	+	susy::PFJetCollection& jetColl = pfJets_it->second;
75	+	for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
76	+	it != jetColl.end(); it++) {
77	+	std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
78	+	if (s_it == it->jecScaleFactors.end()) {
79	+	std::cout << "JEC is not available for this jet!!!" << std::endl;
80	+	continue;
81	+	}
82	+	float scale = s_it->second;
83	+	TLorentzVector corrP4 = scale * it->momentum;
84	+	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
85	+	if (deltaR_ > 0.3) continue;
86	+	if( loggingVerbosity > 0 )
87	+	std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
88	+	nearJets.push_back( *it );
89	+	}// for jet
90	+	}// if, else
91	+
92	+	if ( nearJets.size() == 0 ) {
93	+	//std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
94	+	return myPhoton.momentum.Et();
95	+	}
96	+
97	+	float pt = 0;
98	+	float minPtDifferenz = 1E20; // should be very high
99	+	for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
100	+	it != jetEnd; it++ ) {
101	+	float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
102	+	if (  ptDiff < minPtDifferenz ) {
103	+	minPtDifferenz = ptDiff;
104	+	pt = it->momentum.Et();
105	+	}
106	+	}
107	+
108	+	// testing
109	+	if( nearJets.size() > 1 && loggingVerbosity > 0 )
110	+	std::cout << "There are several jets matching to this photon. "
111	+	<< "Please check if jet-matching is correct." << std::endl;
112	+	return pt;
113	+	}
114	+
115	+
116		void TreeWriter::Loop() {
117	+
118		// here the event loop is implemented and the tree is filled
119		if (inputTree == 0) return;
120
#	Line 102 | Line 132 | void TreeWriter::Loop() {
132		tree->Branch("photon", &photon);
133		tree->Branch("jet", &jet);
134		tree->Branch("met", &met, "met/F");
135	+	tree->Branch("ht", &ht, "ht/F");
136		tree->Branch("nVertex", &nVertex, "nVertex/I");
137	<	tree->Branch("weigth", &weight, "weight/F");
137	>	tree->Branch("nElectron", &nElectron, "nElectron/I");
138
139
140		for (long jentry=0; jentry < processNEvents; jentry++) {
141	<	if (jentry%reportEvery==0)
141	>	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
142		std::cout << jentry << " / " << processNEvents << std :: endl;
143		inputTree->GetEntry(jentry);
144
145		photon.clear();
146		jet.clear();
147	+	ht = 0;
148
149		// photons
150	+	if( loggingVerbosity > 1 )
151	+	std::cout << "Process photons" << std::endl;
152		std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
153		for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
154		it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
155	<	if( ! it->isEB() )
155	>	if( ! it->isEB() && skim )
156		continue;
157	<	thisphoton->pt = it->momentum.Et();
158	<	if( thisphoton->pt < 80 )
157	>	thisphoton->pt = getPtFromMatchedJet( *it, *event );
158	>	if( thisphoton->pt < 80 && skim )
159		continue;
160		thisphoton->eta = it->momentum.Eta();
161	+	thisphoton->phi = it->momentum.Phi();
162		thisphoton->chargedIso = it->chargedHadronIso;
163		thisphoton->neutralIso = it->neutralHadronIso;
164		thisphoton->photonIso = it->photonIso;
165	<	if ( it->r9 > 1 ) // if == 1 ?
165	>	if ( it->r9 > 1 && skim ) // if == 1 ?
166		continue;
167		thisphoton->r9 = it->r9;
168		thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
169		thisphoton->hadTowOverEm = it->hadTowOverEm;
170		thisphoton->pixelseed = it->nPixelSeeds;
171		photon.push_back( *thisphoton );
172	+	ht += thisphoton->pt;
173	+	if( loggingVerbosity > 2 )
174	+	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
175		}
176	<	if( photon.size() == 0 )
176	>
177	>	if( photon.size() == 0 && skim )
178		continue;
179		std::sort( photon.begin(), photon.end(), tree::EtGreater);
180	<
180	>	if( loggingVerbosity > 1 )
181	>	std::cout << "Found " << photon.size() << " photons" << std::endl;
182
183		// jets
184		std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
#	Line 158 | Line 198 | void TreeWriter::Loop() {
198		float scale = s_it->second;
199		TLorentzVector corrP4 = scale * it->momentum;
200
201	<	if(std::abs(corrP4.Eta()) > 3.0) continue;
201	>	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
202		thisjet->pt = corrP4.Et();
203		thisjet->eta = corrP4.Eta();
204	+	thisjet->phi = corrP4.Phi();
205		jet.push_back( *thisjet );
206	+	ht += thisjet->pt;
207		}// for jet
208		}// if, else
209	<	if( jet.size() == 0 )
210	<	std::cout << "error, no jets found " << std::endl;
211	<	else
170	<	std::sort( jet.begin(), jet.end(), tree::EtGreater);
209	>	if( jet.size() < 2 && skim )
210	>	continue;
211	>	std::sort( jet.begin(), jet.end(), tree::EtGreater);
212
213		// met
214		std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
215		susy::MET* metobj = &(met_it->second);
216		met = metobj->met();
217
218	+	// electrons
219	+	//std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
220	+	//nElectron = eVector.size();
221	+
222		// vertices
223
224		nVertex = event->vertices.size();
225	<	weight = 1;
225	>
226	>	if( ht < 450 && skim)
227	>	continue;
228
229		tree->Fill();
230		} // for jentry
#	Line 189 | Line 236 | void TreeWriter::Loop() {
236		outFile->Write();
237		outFile->Close();
238		}
239	<	/*
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	<	*/
239	>

Diff Legend

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