ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Anghel/macros/src/readTree.cc
Revision: 1.2
Committed: Fri Jan 29 23:22:41 2010 UTC (15 years, 3 months ago) by anghel
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +44 -11 lines
Error occurred while calculating annotation data.
Log Message:
Revision 29Jan2010

File Contents

# Content
1 #include "myGlobalMuon.h"
2 #include "myGlobalMuon_LinkDef.h"
3 #include "myJet.h"
4 #include "myJet_LinkDef.h"
5 #include "myParticle.h"
6 #include "myParticle_LinkDef.h"
7 #include "myGenParticle.h"
8 #include "myGenParticle_LinkDef.h"
9 #include "myEvent.h"
10 #include "myEvent_LinkDef.h"
11
12 #include "src/readTreeFunc.cc"
13
14 #include "TChain.h"
15 #include "TFile.h"
16 #include "TH1.h"
17 #include "TH1F.h"
18 #include "TMath.h"
19 #include <iostream>
20 #include <vector>
21 #include <string>
22 #include <fstream>
23
24 using namespace std;
25
26 int main(int argc, char* argv[]) {
27
28 string treeName;
29 string simType;
30 string filter;
31 string studyType;
32 string tagOpPoint;
33 double opP = -999.0;
34
35 //arguments used when running
36 for(int i = 1; i != argc; ++i) {
37 string arg = argv[i];
38 if (arg == "-s") {treeName = string(argv[i+1]); }
39 else if (arg == "-m") {simType = string(argv[i+1]);}
40 else if (arg == "-f") {filter = string(argv[i+1]); }
41 else if (arg == "-p") {studyType = string(argv[i+1]); }
42 else if (arg == "-t") {
43 if (string(argv[i+1]) == "medium") {opP = 2.05; tagOpPoint = argv[i+1];}
44 else if (string(argv[i+1]) == "tight") {opP = 4.07; tagOpPoint = argv[i+1];}
45 else if (string(argv[i+1]) == "none") {opP = -999.0; tagOpPoint = argv[i+1];}}
46
47 else if (arg == "-h") { WhichOption(); return 0; }
48 //else { std::cout << "\nUnknown parameter " << arg << std::endl; WhichOption(); return 0; }
49 }
50
51 cout << treeName << " ,\t" << simType << " ,\t" << filter << " ,\t" << tagOpPoint << endl;
52 //!!!!!!!!!!!!!!change it!!!!!! TFile *infile= TFile::Open(("/uscms_data/d2/omanghel/Analyses/PATtuples/"+realtreeName+"_"+filter+"_RPT.root").c_str());
53
54 string path = "/uscms_data/d2/omanghel/Analyses/PATtuples/";
55
56 //input file
57 TFile *infile;
58
59 //load all the input files
60 string treeRealName;
61 if (simType == "FastSim") treeRealName = treeName+simType;
62 else treeRealName = treeName;
63 TChain *chain = new TChain(treeRealName.c_str(), "");
64 chain->Add((path+treeName+simType+"/"+treeName+simType+"_"+filter+"_"+studyType+"/"+treeName+simType+"_"+filter+"_*_treeP.root").c_str());
65 //TChain *chain = new TChain((treeName+simType).c_str(), "");
66 //chain->Add((path+treeName+simType+"/"+treeName+simType+filter+"/"+"*.root").c_str());
67 //cout << (path+treeName+simType+"/"+treeName+simType+filter+"/"+"*.root").c_str() << endl;
68
69 //get the number of events
70 size_t nrEvents = chain->GetEntries();
71 cout <<"nrEvents = " << nrEvents <<endl;
72
73 //define the histograms that needs to be read directly from the input file
74 TH1F *histo_weight;
75
76 //define histograms put in the output file
77 TH1F histo_weightValue = TH1F("h_weightValue", "", 2, 0., 2.);
78 TH1F histo_jetFlavor = TH1F("h_jetFlavor", "", 2000, 0., 2000.);
79 TH1F histo_tagjetFlavor = TH1F("h_tagjetFlavor", "", 2000, 0., 2000.);
80 TH1F histo_tagBjetFlavor = TH1F("h_tagBjetFlavor", "", 2000, 0., 2000.);
81 TH1F histo_jetMultiplicity = TH1F("h_jetMultiplicity", "", 4, 0., 4.);
82 TH1F histo_jetMultiplicity2 = TH1F("h_jetMultiplicity_4binInclusive", "", 4, 0., 4.);
83 TH1F histo_tagjetMultiplicity = TH1F("h_tagjetMultiplicity", "", 4, 0., 4.);
84 TH1F histo_tagBjetMultiplicity = TH1F("h_tagBjetMultiplicity", "", 4, 0., 4.);
85 TH1F histo_jetPt = TH1F("h_jetPt", "", 500, 0., 500.);
86 TH1F histo_tagjetPt = TH1F("h_tagjetPt", "", 500, 0., 500.);
87 TH1F histo_tagBjetPt = TH1F("h_tagBjetPt", "", 500, 0., 500.);
88 TH1F histo_metEt = TH1F("h_metEt", "", 500, 0., 500.);
89 TH1F histo_bDiscriminator = TH1F("h_bDiscriminator", "", 20, -10., 10.);
90
91 histo_jetFlavor.Sumw2();
92 histo_jetMultiplicity.Sumw2();
93
94 //define the storing objects
95 vector<myJet> *jetPs_ = new vector<myJet>();
96 vector<myGlobalMuon> *muonPs_ = new vector<myGlobalMuon>();
97 vector<myParticle> *metPs_ = new vector<myParticle>();
98 vector<myGenParticle> *genPs_ = new vector<myGenParticle>();
99 vector<myEvent> *evtPs_ = new vector<myEvent>();
100
101 //set the branches
102 chain->SetBranchStatus("*", 0);
103 chain->SetBranchStatus("jetParams*", 1);
104 chain->SetBranchAddress("jetParams", &jetPs_);
105
106 chain->SetBranchStatus("muonParams*", 1);
107 chain->SetBranchAddress("muonParams", &muonPs_);
108
109 chain->SetBranchStatus("metParams*", 1);
110 chain->SetBranchAddress("metParams", &metPs_);
111
112 chain->SetBranchStatus("genParticleParams*", 1);
113 chain->SetBranchAddress("genParticleParams", &genPs_);
114
115 chain->SetBranchStatus("eventParams*", 1);
116 chain->SetBranchAddress("eventParams", &evtPs_);
117
118 typedef vector<myJet>::const_iterator myJetiter;
119 typedef vector<myGlobalMuon>::const_iterator myGlobalMuoniter;
120 typedef vector<myParticle>::const_iterator myParticleiter;
121 typedef vector<myGenParticle>::const_iterator myGenPiter;
122 typedef vector<myEvent>::const_iterator myEventiter;
123
124 //get the weight
125 double weight = 0.0;
126
127 int N1jet = 0;
128 int N2jet = 0;
129 int N3jet = 0;
130 int N4jet = 0;
131
132 int N1jet1tag = 0;
133 int N2jet1tag = 0;
134 int N3jet1tag = 0;
135 int N4jet1tag = 0;
136 int N2jet2tag = 0;
137 int N3jet2tag = 0;
138 int N4jet2tag = 0;
139
140 int NevVectBoson = 0;
141
142 int Nev0bjets[4];
143 int Nev1bjets[4];
144 int Nev2bjets[4];
145 int Nev0tagB[4];
146 int Nev1tagB[4];
147 int Nev2tagB[4];
148
149 for (int i=0; i<4; i++) {
150 Nev0bjets[i] = 0;
151 Nev1bjets[i] = 0;
152 Nev2bjets[i] = 0;
153 Nev0tagB[i] = 0;
154 Nev1tagB[i] = 0;
155 Nev2tagB[i] = 0;
156 }
157
158 //for testing only
159 int NjetsF[9];
160 for (int i=0; i!=9; i++) NjetsF[i] = 0;
161 int NtagJets = 0;
162 int isL;
163 int isT;
164
165 //loop over the number of events
166 //nrEvents = 50;
167 for(size_t n = 0; n != nrEvents; ++n) {
168 //cout <<"----------New Event---------------" << endl;
169
170 int isW = 0;
171 int isZ = 0;
172 isL = 0;
173 isT = 0;
174
175 chain->GetEntry(n);
176
177 if (n==0) {
178 //get input root file
179 infile = chain -> GetCurrentFile();
180 //get histograms already filled in the input file
181 infile->cd();
182 histo_weight = (TH1F*) infile->Get("h_weightFactor");
183 weight = histo_weight->GetBinContent(1);
184 cout << "weight factor = " << weight << endl;
185 }
186
187 /***************************************************************************/
188 /***************************************************************************/
189 /***************************************************************************/
190 /**************************************************************************/
191
192 for (myGenPiter igen = genPs_->begin(); igen != genPs_->end(); ++igen ) {
193 if (fabs(igen->id()) == 24) { isW = 1; break;}
194 if (fabs(igen->id()) == 23) { isZ = 1; break;}
195 if (fabs(igen->p4().Energy()) == fabs(igen->p4().Pz()))
196 cout << "GenParticles :::::: energy = " << igen->p4().Energy() << "\tPz= " << igen->p4().Pz() << "\ty = " << igen->p4().Rapidity() << endl;
197 }
198
199 for (myEventiter ev = evtPs_->begin(); ev != evtPs_->end(); ++ev) {
200 isL = ev->isLoose();
201 isT = ev->isTight();
202 }
203
204
205 //if ((isW == 0)&&(isZ == 0)) {
206 //if (isZ == 1) {
207 //if (isW == 1) {
208 if (isT == 1) {
209 NevVectBoson++;
210 if ((jetPs_->size() > 0) && (muonPs_->size() > 0)) {
211 int Nbjets = 0;
212 int Ntag = 0;
213 int NtagB = 0;
214
215 for (myGlobalMuoniter mu = muonPs_->begin(); mu != muonPs_->end(); ++mu) {
216 if (fabs(mu->p4().Energy()) == fabs(mu->p4().Pz()))
217 cout << "energy = " << mu->p4().Energy() << "\tPz= " << mu->p4().Pz() << "\ty = " << mu->p4().Rapidity() << endl;
218 }
219
220 for (myJetiter jp = jetPs_->begin(); jp != jetPs_->end(); ++jp) {
221 histo_jetFlavor.Fill(fabs(jp->flavor()),weight);
222 histo_jetPt.Fill(jp->p4().Pt(),weight);
223
224 //verify if the jet is a b jet
225 if (fabs(jp->flavor()) == 5) Nbjets++;
226
227 //store the flavor of the jet (tagged and untagged)
228 /*
229 int jetFlavor = fabs(jp->flavor());
230 if (fabs(jp->flavor()) < 7) NjetsF[jetFlavor]++;
231 if (fabs(jp->flavor()) == 21) NjetsF[7]++;
232 if ((fabs(jp->flavor()) > 7)&&(fabs(jp->flavor()) != 21)) NjetsF[8]++;*/
233
234 //verify if the jet is tagged
235 if (jp->bDiscrim() > opP) {
236 Ntag++;
237 histo_tagjetFlavor.Fill(fabs(jp->flavor()),weight);
238 histo_tagjetPt.Fill(jp->p4().Pt(),weight);
239
240 //store the flavor of the tagged jet
241 int jetFlavor = fabs(jp->flavor());
242 if (fabs(jp->flavor()) < 7) NjetsF[jetFlavor]++;
243 if (fabs(jp->flavor()) == 21) NjetsF[7]++;
244 if ((fabs(jp->flavor()) > 7)&&(fabs(jp->flavor()) != 21)) NjetsF[8]++;
245
246 //verify if the jet is a b jet
247 if (fabs(jp->flavor()) == 5) {
248 NtagB++;
249 histo_tagBjetFlavor.Fill(fabs(jp->flavor()),weight);
250 histo_tagBjetPt.Fill(jp->p4().Pt(),weight);
251 }
252 }
253 }
254
255 if (jetPs_->size() < 4) {
256 double content1_beforetagg = histo_jetMultiplicity.GetBinContent(jetPs_->size()) + weight;
257 histo_jetMultiplicity.SetBinContent(jetPs_->size(),content1_beforetagg);
258 if (Ntag > 0) {
259 double content1_aftertagg = histo_jetMultiplicity2.GetBinContent(jetPs_->size()) + weight;
260 histo_jetMultiplicity2.SetBinContent(jetPs_->size(),content1_aftertagg);
261 double content1tag = histo_tagjetMultiplicity.GetBinContent(Ntag) + weight;
262 histo_tagjetMultiplicity.SetBinContent(Ntag,content1tag);
263 }
264 if (NtagB > 0){
265 double content1tagB = histo_tagBjetMultiplicity.GetBinContent(Ntag) + weight;
266 histo_tagBjetMultiplicity.SetBinContent(NtagB,content1tagB);
267 }
268 }
269 else {
270 double content2_beforetagg = histo_jetMultiplicity.GetBinContent(4) + weight;
271 histo_jetMultiplicity.SetBinContent(4,content2_beforetagg);
272 if (Ntag > 0) {
273 double content2_aftertagg = histo_jetMultiplicity2.GetBinContent(4) + weight;
274 //cout << "# of jets = " << jetPs_->size() << " , content2 = " << content2 << endl;
275 histo_jetMultiplicity2.SetBinContent(4,content2_aftertagg);
276 //cout << "Ntag = " << Ntag << "\tNtagB = " << NtagB << endl;
277 if (Ntag >= 4) {
278 double content2tag = histo_tagjetMultiplicity.GetBinContent(4) + weight;
279 histo_tagjetMultiplicity.SetBinContent(4,content2tag);
280 }
281 }
282 if (NtagB == 2) {
283 double content2tagB = histo_tagBjetMultiplicity.GetBinContent(4) + weight;
284 histo_tagBjetMultiplicity.SetBinContent(4,content2tagB);
285 }
286 }
287
288 for (myParticleiter metp = metPs_->begin(); metp != metPs_->end(); ++metp)
289 histo_metEt.Fill(metp->p4().Et(),weight);
290
291 if (jetPs_->size() == 1) N1jet++;
292 if (jetPs_->size() == 2) N2jet++;
293 if (jetPs_->size() == 3) N3jet++;
294 if (jetPs_->size() >= 4) N4jet++;
295
296 if (Ntag==1) {
297 if (jetPs_->size() == 1) N1jet1tag++;
298 if (jetPs_->size() == 2) {
299 N2jet1tag++;
300 /*for (myJetiter jp = jetPs_->begin(); jp != jetPs_->end(); ++jp) {
301 if (jp->bDiscrim() > opP) {
302 int jetFlavor = fabs(jp->flavor());
303 if (fabs(jp->flavor()) < 7) NjetsF[jetFlavor]++;
304 if (fabs(jp->flavor()) == 21) NjetsF[7]++;
305 if ((fabs(jp->flavor()) > 7)&&(fabs(jp->flavor()) != 21)) NjetsF[8]++;
306 }
307 }*/
308 }
309 if (jetPs_->size() == 3) N3jet1tag++;
310 if (jetPs_->size() >= 4) N4jet1tag++;
311 }
312 if (Ntag==2) {
313 if (jetPs_->size() == 2) N2jet2tag++;
314 if (jetPs_->size() == 3) N3jet2tag++;
315 if (jetPs_->size() >= 4) N4jet2tag++;
316 }
317
318 int njets = (jetPs_->size()) - 1;
319 //cout << "Nbjets = " << Nbjets << "\tNtagB = " << NtagB << endl;
320 if (NtagB == 0) Nev0tagB[njets]++;
321 if (NtagB == 1) Nev1tagB[njets]++;
322 if (NtagB == 2) Nev2tagB[njets]++;
323
324 if (Nbjets == 0) Nev0bjets[njets]++;
325 if (Nbjets == 1) Nev1bjets[njets]++;
326 if (Nbjets == 2) Nev2bjets[njets]++;
327
328 }//end if empty collections
329 //}//end if W or Z
330 }//end if isL, isT
331 }//end loop # of events
332
333 cout <<"weight - before filling ---------- = " << weight << endl;
334 histo_weightValue.SetBinContent(1,weight);
335
336 //printed information - some will be removed later
337 cout << "Total # of events with vector bosons = " << NevVectBoson << endl;
338 cout << "The number of events with" << endl;
339 cout << "1jet -> " << N1jet << endl;
340 cout << "2jets -> " << N2jet << endl;
341 cout << "3jets -> " << N3jet << endl;
342 cout << ">= 4jets -> " << N4jet << endl;
343 cout << "-------------------------" << endl;
344
345 cout << "The number of events with" << endl;
346 cout << "1jet,1tag -> " << N1jet1tag << endl;
347 cout << "2jets,1tag -> " << N2jet1tag << endl;
348 cout << "3jets,1tag -> " << N3jet1tag << endl;
349 cout << ">= 4jets,1tag -> " << N4jet1tag << endl;
350 cout << "-------------------------" << endl;
351 cout << "The number of events with" << endl;
352 cout << "2jets,2tags -> " << N2jet2tag << endl;
353 cout << "3jets,2tags -> " << N3jet2tag << endl;
354 cout << ">= 4jets,2tags -> " << N4jet2tag << endl;
355 cout << "-------------------------" << endl;
356
357 for (int i=0; i<4; i++) {
358 cout << "For events with the # of jets = " << i+1 << endl;
359 cout << "Nev*(check) 0 bjets = " << Nev0bjets[i] << endl;
360 cout << "Nev*(check) 1 bjets = " << Nev1bjets[i] << endl;
361 cout << "Nev*(check) 2 bjets = " << Nev2bjets[i] << endl;
362 cout << "-------------------------" << endl;
363 cout << "Nev*(check) 0 tagged bjets = " << Nev0tagB[i] << endl;
364 cout << "Nev*(check) 1 tagged bjets = " << Nev1tagB[i] << endl;
365 cout << "Nev*(check) 2 tagged bjets = " << Nev2tagB[i] << endl;
366 cout << "-------------------------" << endl;
367 }
368
369 cout << "Counting tag jets and their flavors...." << endl;
370 cout << "Total # of tagged jets = " << NtagJets << endl;
371 cout << "Njets with flavor 0 is = " << NjetsF[0] << endl;
372 cout << "Njets with flavor 1 is = " << NjetsF[1] << endl;
373 cout << "Njets with flavor 2 is = " << NjetsF[2] << endl;
374 cout << "Njets with flavor 3 is = " << NjetsF[3] << endl;
375 cout << "Njets with flavor 4 is = " << NjetsF[4] << endl;
376 cout << "Njets with flavor 5 is = " << NjetsF[5] << endl;
377 cout << "Njets with flavor 6 is = " << NjetsF[6] << endl;
378 cout << "Njets with flavor 21 is = " << NjetsF[7] << endl;
379 cout << "Njets with other flavors is = " << NjetsF[8] << endl;
380 cout << "-------------------------" << endl;
381
382 cout << "=====Tagging efficiency====" << endl;
383 cout << "-> 1 jets , 1 tag = " << N1jet1tag*1.0/N1jet << endl;
384 cout << "-> 2 jets , 1 tag = " << N2jet1tag*1.0/N2jet << endl;
385 cout << "-> 3 jets , 1 tag = " << N3jet1tag*1.0/N3jet << endl;
386 cout << "-> 4 jets , 1 tag = " << N4jet1tag*1.0/N4jet << endl;
387 cout << "-> 2 jets , 2 tag = " << N2jet2tag*1.0/N2jet << endl;
388 cout << "-> 3 jets , 2 tag = " << N3jet2tag*1.0/N3jet << endl;
389 cout << "-> 4 jets , 2 tag = " << N4jet2tag*1.0/N4jet << endl;
390 cout << "-> 2 jets , any # of tagged jets = " << (N2jet1tag+N2jet2tag)*1.0/N2jet << endl;
391 cout << "-> 3 jets , any # of tagged jets = " << (N3jet1tag+N3jet2tag)*1.0/N3jet << endl;
392 cout << "-> 4 jets , any # of tagged jets = " << (N4jet1tag+N4jet2tag)*1.0/N4jet << endl;
393
394 //write the statistical information in a file
395 ofstream textFile;
396 string textFileName = treeName+simType+"_"+filter+"_TightSel_"+tagOpPoint+".txt";
397 textFile.open(textFileName.c_str());
398 textFile << N2jet1tag << "\t" << N2jet2tag << endl;
399 textFile << N3jet1tag << "\t" << N3jet2tag << endl;
400 textFile << N4jet1tag << "\t" << N4jet2tag << endl;
401 textFile.close();
402
403 //write the histograms into the root file
404 TFile *ofile;
405 ofile = new TFile(("histosFileTag_"+tagOpPoint+".root").c_str(), "UPDATE");
406 ofile->cd();
407 string jetFlavorName = histo_jetFlavor.GetName();
408 histo_jetFlavor.Write((jetFlavorName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
409 string tagjetFlavorName = histo_tagjetFlavor.GetName();
410 histo_tagjetFlavor.Write((tagjetFlavorName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
411 string jetMultiplicityName = histo_jetMultiplicity.GetName();
412 histo_jetMultiplicity.Write((jetMultiplicityName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
413 string jetMultiplicityName2 = histo_jetMultiplicity2.GetName();
414 histo_jetMultiplicity2.Write((jetMultiplicityName2+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
415 string tagjetMultiplicityName = histo_tagjetMultiplicity.GetName();
416 histo_tagjetMultiplicity.Write((tagjetMultiplicityName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
417 string tagBjetMultiplicityName = histo_tagBjetMultiplicity.GetName();
418 histo_tagBjetMultiplicity.Write((tagBjetMultiplicityName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
419 string jetPtName = histo_jetPt.GetName();
420 histo_jetPt.Write((jetPtName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
421 string tagjetPtName = histo_tagjetPt.GetName();
422 histo_tagjetPt.Write((tagjetPtName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
423 string tagBjetPtName = histo_tagBjetPt.GetName();
424 histo_tagBjetPt.Write((tagBjetPtName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
425 string metEtName = histo_metEt.GetName();
426 histo_metEt.Write((metEtName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
427 string weightName = histo_weightValue.GetName();
428 histo_weightValue.Write((weightName+"_"+treeName+"_"+filter).c_str(), TObject::kOverwrite);
429
430 return 0;
431 }