kiesel/TreeWriter/treeWriter.cc

#include<iostream>
#include<math.h>
#include<string>

#include "TSystem.h"

#include "treeWriter.h"

using namespace std;

TreeWriter::TreeWriter(TString inputName, TString outputName, int loggingVerbosity_ ) {
        // read the input file
        inputTree = new TChain("susyTree");
        if (loggingVerbosity_ > 0)
                std::cout << "Add files to chain" << std::endl;
        inputTree->Add( inputName );

        if (loggingVerbosity_ > 0)
                std::cout << "Set Branch Address of susy::Event" << std::endl;
        event = new susy::Event;
        inputTree->SetBranchAddress("susyEvent", &event);

        if (outputName == "") {
                // Here the output filename is generated automaticly
                try{
                        string fName = (string) inputName;
                        // Search for 'susyEvents_', and get the unique characters afterwards.
                        // eg /path/susyEvents_123_Myl.root -> susyTree_123_Mly.root
                        outputName = "susyTree_" + fName.substr(fName.find("susyEvents_")+11,-1);
                } catch (int e ) {
                        outputName = "susyTree.root";
                }
        }

        if (loggingVerbosity_>0)
                std::cout << "Open file " << outputName << " for writing." << std::endl;
        // open the output file
        outFile = new TFile( outputName, "recreate" );
        tree = new TTree("susyTree","Tree for single photon analysis");

        // set default parameter
        processNEvents = -1;
        reportEvery = 1000;
        loggingVerbosity = loggingVerbosity_;
        skim = true;

}

TreeWriter::~TreeWriter() {
        if (!inputTree) return;
        delete inputTree->GetCurrentFile();
}

float TreeWriter::deltaR( TLorentzVector v1, TLorentzVector v2 ) {
        return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) );
}

float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
        /**
         * \brief Takes jet p_T as photon p_T
         *
         * At first all jets with DeltaR < 0.3 (isolation cone) are searched.
         * If several jets are found, take the one with the minimal pt difference
         * compared to the photon. If no such jets are found, keep the photon_pt
         * TODO: remove photon matched jet from jet-selection?
         */
        std::vector<susy::PFJet> nearJets;
        nearJets.clear();

        std::map<TString,susy::PFJetCollection>::iterator pfJets_it = myEvent.pfJets.find("ak5");
        if(pfJets_it == myEvent.pfJets.end()){
                if(myEvent.pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
        } else {
                susy::PFJetCollection& jetColl = pfJets_it->second;
                for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
                                it != jetColl.end(); it++) {
                        std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
                        if (s_it == it->jecScaleFactors.end()) {
                                std::cout << "JEC is not available for this jet!!!" << std::endl;
                                continue;
                        }
                        float scale = s_it->second;
                        TLorentzVector corrP4 = scale * it->momentum;
                        float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
                        if (deltaR_ > 0.3) continue;
                        if( loggingVerbosity > 0 )
                                std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
                        nearJets.push_back( *it );
                }// for jet
        }// if, else

        if ( nearJets.size() == 0 ) {
                //std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
                return myPhoton.momentum.Et();
        }

        float pt = 0;
        float minPtDifferenz = 1E20; // should be very high
        for( std::vector<susy::PFJet>::iterator it = nearJets.begin(), jetEnd = nearJets.end();
                        it != jetEnd; it++ ) {
                float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
                if (  ptDiff < minPtDifferenz ) {
                        minPtDifferenz = ptDiff;
                        pt = it->momentum.Et();
                }
        }

        // testing
        if( nearJets.size() > 1 && loggingVerbosity > 0 )
                std::cout << "There are several jets matching to this photon. "
                                        << "Please check if jet-matching is correct." << std::endl;
        return pt;
}


void TreeWriter::Loop() {

        // here the event loop is implemented and the tree is filled
        if (inputTree == 0) return;

        // get number of events to be proceeded
        Long64_t nentries = inputTree->GetEntries();
        if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries;

        if( loggingVerbosity > 0 )
                std::cout << "Processing " << processNEvents << " ouf of "
                        << nentries << " events. " << std::endl;

        tree::Photon *thisphoton = new tree::Photon();
        tree::Jet *thisjet = new tree::Jet();

        tree->Branch("photon", &photon);
        tree->Branch("jet", &jet);
        tree->Branch("met", &met, "met/F");
        tree->Branch("ht", &ht, "ht/F");
        tree->Branch("nVertex", &nVertex, "nVertex/I");
        tree->Branch("nElectron", &nElectron, "nElectron/I");


        for (long jentry=0; jentry < processNEvents; jentry++) {
                if ( loggingVerbosity>0 && jentry%reportEvery==0 )
                        std::cout << jentry << " / " << processNEvents << std :: endl;
                inputTree->GetEntry(jentry);

                photon.clear();
                jet.clear();
                ht = 0;

                // photons
                if( loggingVerbosity > 1 )
                        std::cout << "Process photons" << std::endl;
                std::map<TString, std::vector<susy::Photon> >::iterator phoMap = event->photons.find("photons");
                for(std::vector<susy::Photon>::iterator it = phoMap->second.begin();
                                it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) {
                        if( ! it->isEB() && skim )
                                continue;
                        thisphoton->pt = getPtFromMatchedJet( *it, *event );
                        if( thisphoton->pt < 80 && skim )
                                continue;
                        thisphoton->eta = it->momentum.Eta();
                        thisphoton->phi = it->momentum.Phi();
                        thisphoton->chargedIso = it->chargedHadronIso;
                        thisphoton->neutralIso = it->neutralHadronIso;
                        thisphoton->photonIso = it->photonIso;
                        if ( it->r9 > 1 && skim ) // if == 1 ?
                                continue;
                        thisphoton->r9 = it->r9;
                        thisphoton->sigmaIetaIeta = it->sigmaIetaIeta;
                        thisphoton->hadTowOverEm = it->hadTowOverEm;
                        thisphoton->pixelseed = it->nPixelSeeds;
                        photon.push_back( *thisphoton );
                        ht += thisphoton->pt;
                        if( loggingVerbosity > 2 )
                                std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
                }

                if( photon.size() == 0 && skim )
                        continue;
                std::sort( photon.begin(), photon.end(), tree::EtGreater);
                if( loggingVerbosity > 1 )
                        std::cout << "Found " << photon.size() << " photons" << std::endl;

                // jets
                std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
                if(pfJets_it == event->pfJets.end()){
                        if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
                } else {

                        susy::PFJetCollection& jetColl = pfJets_it->second;

                        for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
                                        it != jetColl.end(); it++) {
                                std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
                                if (s_it == it->jecScaleFactors.end()) {
                                        std::cout << "JEC is not available for this jet!!!" << std::endl;
                                        continue;
                                }
                                float scale = s_it->second;
                                TLorentzVector corrP4 = scale * it->momentum;

                                if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
                                thisjet->pt = corrP4.Et();
                                thisjet->eta = corrP4.Eta();
                                thisjet->phi = corrP4.Phi();
                                jet.push_back( *thisjet );
                                ht += thisjet->pt;
                        }// for jet
                }// if, else
                if( jet.size() < 2 && skim )
                        continue;
                std::sort( jet.begin(), jet.end(), tree::EtGreater);

                // met
                std::map<TString, susy::MET>::iterator met_it = event->metMap.find("pfMet");
                susy::MET* metobj = &(met_it->second);
                met = metobj->met();

                // electrons
                //std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
                //nElectron = eVector.size();

                // vertices

                nVertex = event->vertices.size();

                if( ht < 450 && skim)
                        continue;

                tree->Fill();
        } // for jentry


        tree->Write();
        outFile->cd();
        outFile->Write();
        outFile->Close();
}

Revision:	1.7
Committed:	Tue Apr 9 12:13:55 2013 UTC (12 years ago) by kiesel
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+15 -0 lines
Log Message:	passing filename as command line argument
#	User	Rev	Content
1	kiesel	1.2	#include<iostream>
2	kiesel	1.4	#include<math.h>
3	kiesel	1.7	#include<string>
4	kiesel	1.2
5	kiesel	1.5	#include "TSystem.h"
6
7			#include "treeWriter.h"
8
9			using namespace std;
10	kiesel	1.1
11	kiesel	1.6	TreeWriter::TreeWriter(TString inputName, TString outputName, int loggingVerbosity_ ) {
12	kiesel	1.1	// read the input file
13	kiesel	1.5	inputTree = new TChain("susyTree");
14	kiesel	1.6	if (loggingVerbosity_ > 0)
15			std::cout << "Add files to chain" << std::endl;
16	kiesel	1.5	inputTree->Add( inputName );
17
18	kiesel	1.6	if (loggingVerbosity_ > 0)
19			std::cout << "Set Branch Address of susy::Event" << std::endl;
20	kiesel	1.1	event = new susy::Event;
21			inputTree->SetBranchAddress("susyEvent", &event);
22
23	kiesel	1.7	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	kiesel	1.1	// open the output file
38			outFile = new TFile( outputName, "recreate" );
39			tree = new TTree("susyTree","Tree for single photon analysis");
40
41			// set default parameter
42			processNEvents = -1;
43			reportEvery = 1000;
44	kiesel	1.6	loggingVerbosity = loggingVerbosity_;
45	kiesel	1.5	skim = true;
46	kiesel	1.1
47			}
48
49			TreeWriter::~TreeWriter() {
50			if (!inputTree) return;
51			delete inputTree->GetCurrentFile();
52			}
53
54	kiesel	1.4	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	kiesel	1.5	float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) {
59	kiesel	1.4	/**
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	kiesel	1.5	* compared to the photon. If no such jets are found, keep the photon_pt
65	kiesel	1.4	* TODO: remove photon matched jet from jet-selection?
66			*/
67			std::vector<susy::PFJet> nearJets;
68			nearJets.clear();
69
70	kiesel	1.5	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	kiesel	1.4	} 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	kiesel	1.5	float deltaR_ = deltaR(myPhoton.momentum, corrP4 );
85	kiesel	1.4	if (deltaR_ > 0.3) continue;
86	kiesel	1.5	if( loggingVerbosity > 0 )
87			std::cout << "gamma pt jet matching factor = " << it->momentum.Et() / myPhoton.momentum.Et() << std::endl;
88	kiesel	1.4	nearJets.push_back( *it );
89			}// for jet
90			}// if, else
91
92			if ( nearJets.size() == 0 ) {
93	kiesel	1.5	//std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl;
94			return myPhoton.momentum.Et();
95	kiesel	1.4	}
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	kiesel	1.5	float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et());
102	kiesel	1.4	if ( ptDiff < minPtDifferenz ) {
103			minPtDifferenz = ptDiff;
104			pt = it->momentum.Et();
105			}
106			}
107
108			// testing
109	kiesel	1.6	if( nearJets.size() > 1 && loggingVerbosity > 0 )
110	kiesel	1.4	std::cout << "There are several jets matching to this photon. "
111	kiesel	1.6	<< "Please check if jet-matching is correct." << std::endl;
112	kiesel	1.4	return pt;
113			}
114
115
116	kiesel	1.1	void TreeWriter::Loop() {
117	kiesel	1.4
118	kiesel	1.1	// here the event loop is implemented and the tree is filled
119			if (inputTree == 0) return;
120
121			// get number of events to be proceeded
122			Long64_t nentries = inputTree->GetEntries();
123			if(processNEvents <= 0 \|\| processNEvents > nentries) processNEvents = nentries;
124
125			if( loggingVerbosity > 0 )
126			std::cout << "Processing " << processNEvents << " ouf of "
127			<< nentries << " events. " << std::endl;
128
129			tree::Photon *thisphoton = new tree::Photon();
130			tree::Jet *thisjet = new tree::Jet();
131
132			tree->Branch("photon", &photon);
133			tree->Branch("jet", &jet);
134			tree->Branch("met", &met, "met/F");
135	kiesel	1.5	tree->Branch("ht", &ht, "ht/F");
136	kiesel	1.1	tree->Branch("nVertex", &nVertex, "nVertex/I");
137	kiesel	1.4	tree->Branch("nElectron", &nElectron, "nElectron/I");
138	kiesel	1.1
139
140			for (long jentry=0; jentry < processNEvents; jentry++) {
141	kiesel	1.4	if ( loggingVerbosity>0 && jentry%reportEvery==0 )
142	kiesel	1.1	std::cout << jentry << " / " << processNEvents << std :: endl;
143			inputTree->GetEntry(jentry);
144
145			photon.clear();
146			jet.clear();
147	kiesel	1.5	ht = 0;
148	kiesel	1.1
149			// photons
150	kiesel	1.5	if( loggingVerbosity > 1 )
151	kiesel	1.6	std::cout << "Process photons" << std::endl;
152	kiesel	1.1	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	kiesel	1.4	if( ! it->isEB() && skim )
156	kiesel	1.1	continue;
157	kiesel	1.4	thisphoton->pt = getPtFromMatchedJet( it, event );
158			if( thisphoton->pt < 80 && skim )
159	kiesel	1.1	continue;
160			thisphoton->eta = it->momentum.Eta();
161	kiesel	1.4	thisphoton->phi = it->momentum.Phi();
162	kiesel	1.1	thisphoton->chargedIso = it->chargedHadronIso;
163			thisphoton->neutralIso = it->neutralHadronIso;
164			thisphoton->photonIso = it->photonIso;
165	kiesel	1.4	if ( it->r9 > 1 && skim ) // if == 1 ?
166	kiesel	1.1	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	kiesel	1.5	ht += thisphoton->pt;
173			if( loggingVerbosity > 2 )
174	kiesel	1.6	std::cout << " p_T, gamma = " << thisphoton->pt << std::endl;
175	kiesel	1.1	}
176	kiesel	1.5
177	kiesel	1.4	if( photon.size() == 0 && skim )
178	kiesel	1.1	continue;
179			std::sort( photon.begin(), photon.end(), tree::EtGreater);
180	kiesel	1.5	if( loggingVerbosity > 1 )
181			std::cout << "Found " << photon.size() << " photons" << std::endl;
182	kiesel	1.1
183			// jets
184			std::map<TString,susy::PFJetCollection>::iterator pfJets_it = event->pfJets.find("ak5");
185			if(pfJets_it == event->pfJets.end()){
186			if(event->pfJets.size() > 0) std::cout << "JetCollection is not available!!!" << std::endl;
187			} else {
188
189			susy::PFJetCollection& jetColl = pfJets_it->second;
190
191			for(std::vector<susy::PFJet>::iterator it = jetColl.begin();
192			it != jetColl.end(); it++) {
193			std::map<TString,Float_t>::iterator s_it = it->jecScaleFactors.find("L2L3");
194			if (s_it == it->jecScaleFactors.end()) {
195			std::cout << "JEC is not available for this jet!!!" << std::endl;
196			continue;
197			}
198			float scale = s_it->second;
199			TLorentzVector corrP4 = scale * it->momentum;
200
201	kiesel	1.4	if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue;
202	kiesel	1.1	thisjet->pt = corrP4.Et();
203			thisjet->eta = corrP4.Eta();
204	kiesel	1.4	thisjet->phi = corrP4.Phi();
205	kiesel	1.1	jet.push_back( *thisjet );
206	kiesel	1.5	ht += thisjet->pt;
207	kiesel	1.1	}// for jet
208			}// if, else
209	kiesel	1.4	if( jet.size() < 2 && skim )
210			continue;
211			std::sort( jet.begin(), jet.end(), tree::EtGreater);
212	kiesel	1.1
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	kiesel	1.4	// electrons
219	kiesel	1.5	//std::vector<susy::Electron> eVector = event->electrons["gsfElectrons"];
220			//nElectron = eVector.size();
221	kiesel	1.4
222	kiesel	1.1	// vertices
223
224			nVertex = event->vertices.size();
225
226	kiesel	1.5	if( ht < 450 && skim)
227			continue;
228	kiesel	1.4
229	kiesel	1.1	tree->Fill();
230			} // for jentry
231
232
233
234			tree->Write();
235			outFile->cd();
236			outFile->Write();
237			outFile->Close();
238			}
239	kiesel	1.3