ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/csander/HEPTutorial/MyAnalysis.C
Revision: 1.9
Committed: Tue Feb 28 22:48:16 2012 UTC (13 years, 2 months ago) by csander
Content type: text/plain
Branch: MAIN
Changes since 1.8: +188 -87 lines
Log Message:
added more plots for bg/signal comparisons

File Contents

# User Rev Content
1 csander 1.1 #define MyAnalysis_cxx
2     // The class definition in MyAnalysis.h has been generated automatically
3     // by the ROOT utility TTree::MakeSelector(). This class is derived
4     // from the ROOT class TSelector. For more information on the TSelector
5     // framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
6    
7     // The following methods are defined in this file:
8     // Begin(): called every time a loop on the tree starts,
9     // a convenient place to create your histograms.
10     // SlaveBegin(): called after Begin(), when on PROOF called only on the
11     // slave servers.
12     // Process(): called for each event, in this function you decide what
13     // to read and fill your histograms.
14     // SlaveTerminate: called at the end of the loop on the tree, when on PROOF
15     // called only on the slave servers.
16     // Terminate(): called at the end of the loop on the tree,
17     // a convenient place to draw/fit your histograms.
18     //
19     // To use this file, try the following session on your Tree T:
20     //
21     // Root > T->Process("MyAnalysis.C")
22     // Root > T->Process("MyAnalysis.C","some options")
23     // Root > T->Process("MyAnalysis.C+")
24     //
25    
26     #include "MyAnalysis.h"
27 csander 1.3 #include <iostream>
28 csander 1.2 #include <TH1F.h>
29 csander 1.9 #include <TLatex.h>
30 csander 1.1
31 csander 1.3 using namespace std;
32    
33     void MyAnalysis::BuildEvent() {
34     Jets.clear();
35     for (int i = 0; i < NJet; ++i) {
36     MyJet jet(Jet_Px[i], Jet_Py[i], Jet_Pz[i], Jet_E[i]);
37     jet.SetBTagDiscriminator(Jet_btag[i]);
38     Jets.push_back(jet);
39     }
40     Muons.clear();
41     for (int i = 0; i < NMuon; ++i) {
42     MyMuon muon(Muon_Px[i], Muon_Py[i], Muon_Pz[i], Muon_E[i]);
43     muon.SetIsolation(Muon_Iso[i]);
44 csander 1.4 muon.SetCharge(Muon_Charge[i]);
45 csander 1.3 Muons.push_back(muon);
46     }
47 csander 1.4 Electrons.clear();
48     for (int i = 0; i < NElectron; ++i) {
49     MyElectron electron(Electron_Px[i], Electron_Py[i], Electron_Pz[i], Electron_E[i]);
50     electron.SetIsolation(Electron_Iso[i]);
51     electron.SetCharge(Electron_Charge[i]);
52     Electrons.push_back(electron);
53     }
54     Photons.clear();
55     for (int i = 0; i < NPhoton; ++i) {
56     MyPhoton photon(Photon_Px[i], Photon_Py[i], Photon_Pz[i], Photon_E[i]);
57     photon.SetIsolation(Photon_Iso[i]);
58     Photons.push_back(photon);
59     }
60 csander 1.6 hadB.SetXYZM(MChadronicBottom_px, MChadronicBottom_py, MChadronicBottom_pz, 4.8);
61     lepB.SetXYZM(MCleptonicBottom_px, MCleptonicBottom_py, MCleptonicBottom_pz, 4.8);
62     hadWq.SetXYZM(MChadronicWDecayQuark_px, MChadronicWDecayQuark_py, MChadronicWDecayQuark_pz, 0.0);
63     hadWqb.SetXYZM(MChadronicWDecayQuarkBar_px, MChadronicWDecayQuarkBar_py, MChadronicWDecayQuarkBar_pz, 0.0);
64     lepWl.SetXYZM(MClepton_px, MClepton_py, MClepton_pz, 0.0);
65     lepWn.SetXYZM(MCneutrino_px, MCneutrino_py, MCneutrino_pz, 0.0);
66 csander 1.8 met.SetXYZM(MET_px, MET_py, 0., 0.);
67 csander 1.9
68     if (EventWeight == 0.)
69     EventWeight = 1.;
70     EventWeight *= weight_factor;
71    
72 csander 1.3 }
73    
74 csander 1.2 void MyAnalysis::Begin(TTree * /*tree*/) {
75 csander 1.1 // The Begin() function is called at the start of the query.
76     // When running with PROOF Begin() is only called on the client.
77     // The tree argument is deprecated (on PROOF 0 is passed).
78    
79     TString option = GetOption();
80    
81     }
82    
83 csander 1.2 void MyAnalysis::SlaveBegin(TTree * /*tree*/) {
84 csander 1.1 // The SlaveBegin() function is called after the Begin() function.
85     // When running with PROOF SlaveBegin() is called on each slave server.
86     // The tree argument is deprecated (on PROOF 0 is passed).
87    
88     TString option = GetOption();
89    
90 csander 1.7 h_Mmumu = new TH1F("Mmumu", "Invariant di-muon mass", 50, 0, 150);
91 csander 1.4 h_Mmumu->SetXTitle("m_{#mu#mu}");
92     h_Mmumu->SetYTitle("a.u.");
93     h_Mmumu->Sumw2();
94    
95 csander 1.7 h_Mbqqb_mc = new TH1F("Mbqqb_mc", "Invariant hadronic top mass", 100, 170, 175);
96     h_Mbqqb_mc->SetXTitle("m_{hadronic top}");
97     h_Mbqqb_mc->SetYTitle("a.u.");
98     h_Mbqqb_mc->Sumw2();
99    
100     h_Mbln_mc = new TH1F("Mbln_mc", "Invariant leptonic top mass", 100, 170, 175);
101     h_Mbln_mc->SetXTitle("m_{leptonic top}");
102     h_Mbln_mc->SetYTitle("a.u.");
103     h_Mbln_mc->Sumw2();
104    
105     h_Mbqqb_reco = new TH1F("Mbqqb_reco", "Invariant hadronic top mass", 250, 0, 500);
106     h_Mbqqb_reco->SetXTitle("m_{hadronic top}");
107     h_Mbqqb_reco->SetYTitle("a.u.");
108     h_Mbqqb_reco->Sumw2();
109    
110 csander 1.8 h_Mbln_reco = new TH1F("Mbln_reco", "Invariant transverse leptonic top mass", 250, 0, 500);
111     h_Mbln_reco->SetXTitle("m_{T, leptonic top}");
112     h_Mbln_reco->SetYTitle("a.u.");
113     h_Mbln_reco->Sumw2();
114    
115 csander 1.9 h_NJet = new TH1F("NJet", "Number of jets", 7, 0, 7);
116     h_NJet->SetXTitle("NJet");
117     h_NJet->Sumw2();
118    
119     h_NBJet = new TH1F("NBJet", "Number of b-jets", 7, 0, 7);
120     h_NBJet->SetXTitle("NBJet");
121     h_NBJet->Sumw2();
122    
123     h_Jet1_Pt = new TH1F("Jet1Pt", "Pt of jet1", 100, 0, 250);
124     h_Jet1_Pt->SetXTitle("P_t [GeV]");
125     h_Jet1_Pt->Sumw2();
126    
127     h_Jet2_Pt = new TH1F("Jet2Pt", "Pt of jet2", 100, 0, 250);
128     h_Jet2_Pt->SetXTitle("P_t [GeV]");
129     h_Jet2_Pt->Sumw2();
130    
131     h_Jet3_Pt = new TH1F("Jet3Pt", "Pt of jet3", 100, 0, 250);
132     h_Jet3_Pt->SetXTitle("P_t [GeV]");
133     h_Jet3_Pt->Sumw2();
134    
135     h_Jet1_Eta = new TH1F("Jet1Eta", "#Eta of jet1", 50, -5, 5);
136     h_Jet1_Eta->SetXTitle("#eta");
137     h_Jet1_Eta->Sumw2();
138    
139     h_Jet2_Eta = new TH1F("Jet2Eta", "#Eta of jet2", 50, -5, 5);
140     h_Jet2_Eta->SetXTitle("#eta");
141     h_Jet2_Eta->Sumw2();
142    
143     h_Jet3_Eta = new TH1F("Jet3Eta", "#Eta of jet3", 50, -5, 5);
144     h_Jet3_Eta->SetXTitle("#eta");
145     h_Jet3_Eta->Sumw2();
146    
147     h_BJet1_Pt = new TH1F("BJet1Pt", "Pt of b-jet1", 100, 0, 250);
148     h_BJet1_Pt->SetXTitle("P_t [GeV]");
149     h_BJet1_Pt->Sumw2();
150    
151     h_BJet2_Pt = new TH1F("BJet2Pt", "Pt of b-jet2", 100, 0, 250);
152     h_BJet2_Pt->SetXTitle("P_t [GeV]");
153     h_BJet2_Pt->Sumw2();
154    
155     h_BJet1_Eta = new TH1F("BJet1Eta", "#Eta of b-jet1", 50, -5, 5);
156     h_BJet1_Eta->SetXTitle("#eta");
157     h_BJet1_Eta->Sumw2();
158    
159     h_BJet2_Eta = new TH1F("BJet2Eta", "#Eta of b-jet2", 50, -5, 5);
160     h_BJet2_Eta->SetXTitle("#eta");
161     h_BJet2_Eta->Sumw2();
162    
163     h_Muon1_Pt = new TH1F("Muon1Pt", "Pt of muon1", 100, 0, 250);
164     h_Muon1_Pt->SetXTitle("P_t [GeV]");
165     h_Muon1_Pt->Sumw2();
166    
167     h_Muon2_Pt = new TH1F("Muon2Pt", "Pt of muon2", 100, 0, 250);
168     h_Muon2_Pt->SetXTitle("P_t [GeV]");
169     h_Muon2_Pt->Sumw2();
170    
171     h_Muon1_Eta = new TH1F("Muon1Eta", "#Eta of muon1", 50, -5, 5);
172     h_Muon1_Eta->SetXTitle("#eta");
173     h_Muon1_Eta->Sumw2();
174    
175     h_Muon2_Eta = new TH1F("Muon2Eta", "#Eta of muon2", 50, -5, 5);
176     h_Muon2_Eta->SetXTitle("#eta");
177     h_Muon2_Eta->Sumw2();
178    
179     h_Muon1_Iso = new TH1F("Muon1Iso", "Isolation of muon1", 50, 0, 25);
180     h_Muon1_Iso->SetXTitle("[GeV]");
181     h_Muon1_Iso->Sumw2();
182    
183     h_Muon2_Iso = new TH1F("Muon2Iso", "Isolation of muon2", 50, 0, 25);
184     h_Muon2_Iso->SetXTitle("[GeV]");
185     h_Muon2_Iso->Sumw2();
186    
187     h_MET = new TH1F("MET", "MET", 30, 0, 300.);
188     h_MET->SetXTitle("MET");
189     h_MET->Sumw2();
190    
191     h_minDeltaPhi_MET_Muon = new TH1F("minDeltaPhi_MET_Muon", "min#Delta#phi(MET,muon)", 50, -TMath::Pi(), TMath::Pi());
192     h_minDeltaPhi_MET_Muon->SetXTitle("min#Delta#phi [rad]");
193     h_minDeltaPhi_MET_Muon->Sumw2();
194    
195     h_minDeltaPhi_MET_BJet = new TH1F("minDeltaPhi_MET_BJet", "min#Delta#phi(MET,b-jet)", 50, -TMath::Pi(), TMath::Pi());
196     h_minDeltaPhi_MET_BJet->SetXTitle("min#Delta#phi [rad]");
197     h_minDeltaPhi_MET_BJet->Sumw2();
198    
199     h_selectedEvents_Muon1_Pt = new TH1F("selectedEvents_Muon1_Pt", "Pt of muon1 (selected)", 100, 0, 250);
200     h_selectedEvents_Muon1_Pt->SetXTitle("P_t [GeV]");
201     h_selectedEvents_Muon1_Pt->Sumw2();
202    
203     h_selectedEvents_triggered_Muon1_Pt = new TH1F("selectedEvents_triggered_Muon1_Pt",
204     "Pt of muon1 (selected and triggered)", 100, 0, 250);
205     h_selectedEvents_triggered_Muon1_Pt->SetXTitle("P_t [GeV]");
206     h_selectedEvents_triggered_Muon1_Pt->Sumw2();
207    
208 csander 1.1 }
209    
210 csander 1.2 Bool_t MyAnalysis::Process(Long64_t entry) {
211 csander 1.1 // The Process() function is called for each entry in the tree (or possibly
212     // keyed object in the case of PROOF) to be processed. The entry argument
213     // specifies which entry in the currently loaded tree is to be processed.
214     // It can be passed to either MyAnalysis::GetEntry() or TBranch::GetEntry()
215     // to read either all or the required parts of the data. When processing
216     // keyed objects with PROOF, the object is already loaded and is available
217     // via the fObject pointer.
218     //
219     // This function should contain the "body" of the analysis. It can contain
220     // simple or elaborate selection criteria, run algorithms on the data
221     // of the event and typically fill histograms.
222     //
223     // The processing can be stopped by calling Abort().
224     //
225     // Use fStatus to set the return value of TTree::Process().
226     //
227     // The return value is currently not used.
228    
229 csander 1.3 ++TotalEvents;
230    
231     GetEntry(entry);
232    
233 csander 1.4 if (TotalEvents % 10000 == 0)
234     cout << "Next event -----> " << TotalEvents << endl;
235    
236 csander 1.3 BuildEvent();
237    
238 csander 1.7 // cout << "Jets: " << endl;
239     // for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
240     // cout << "pt, eta, phi, btag: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", " << it->IsBTagged()
241     // << endl;
242     // }
243     // cout << "Muons: " << endl;
244     // for (vector<MyMuon>::iterator it = Muons.begin(); it != Muons.end(); ++it) {
245     // cout << "pt, eta, phi, iso, charge: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", "
246     // << it->GetIsolation() << ", " << it->GetCharge() << endl;
247     // }
248     // cout << "Electrons: " << endl;
249     // for (vector<MyElectron>::iterator it = Electrons.begin(); it != Electrons.end(); ++it) {
250     // cout << "pt, eta, phi, iso, charge: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", "
251     // << it->GetIsolation() << ", " << it->GetCharge() << endl;
252     // }
253     // cout << "Photons: " << endl;
254     // for (vector<MyPhoton>::iterator it = Photons.begin(); it != Photons.end(); ++it) {
255     // cout << "pt, eta, phi, iso: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", " << it->GetIsolation()
256     // << endl;
257     // }
258    
259     if (NMuon > 1 && Muons.at(0).GetCharge() * Muons.at(1).GetCharge() < 1) {
260     h_Mmumu->Fill((Muons.at(0) + Muons.at(1)).M(), EventWeight);
261 csander 1.5 }
262 csander 1.4
263 csander 1.9 h_NJet->Fill(Jets.size(), EventWeight);
264 csander 1.7
265 csander 1.9 int N_BJet = 0;
266     int N_Jet = 0;
267 csander 1.7 for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
268 csander 1.9 ++N_Jet;
269     if (N_Jet == 1) {
270     h_Jet1_Pt->Fill(it->Pt(), EventWeight);
271     h_Jet1_Eta->Fill(it->Eta(), EventWeight);
272     } else if (N_Jet == 2) {
273     h_Jet2_Pt->Fill(it->Pt(), EventWeight);
274     h_Jet2_Eta->Fill(it->Eta(), EventWeight);
275     } else if (N_Jet == 3) {
276     h_Jet2_Pt->Fill(it->Pt(), EventWeight);
277     h_Jet2_Eta->Fill(it->Eta(), EventWeight);
278     }
279 csander 1.7 if (it->IsBTagged()) {
280 csander 1.9 ++N_BJet;
281     if (N_BJet == 1) {
282     h_BJet1_Pt->Fill(it->Pt(), EventWeight);
283     h_BJet1_Eta->Fill(it->Eta(), EventWeight);
284     } else if (N_BJet == 2) {
285     h_BJet2_Pt->Fill(it->Pt(), EventWeight);
286     h_BJet2_Eta->Fill(it->Eta(), EventWeight);
287 csander 1.7 }
288     }
289 csander 1.3 }
290 csander 1.9 h_NBJet->Fill(N_BJet, EventWeight);
291    
292     int N_Muons = 0;
293     for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
294     ++N_Muons;
295     if (N_Muons == 1) {
296     h_Muon1_Pt->Fill(jt->Pt(), EventWeight);
297     h_Muon1_Eta->Fill(jt->Eta(), EventWeight);
298     h_Muon1_Iso->Fill(jt->GetIsolation(), EventWeight);
299     } else if (N_Muons == 2) {
300     h_Muon2_Pt->Fill(jt->Pt(), EventWeight);
301     h_Muon2_Eta->Fill(jt->Eta(), EventWeight);
302     h_Muon2_Iso->Fill(jt->GetIsolation(), EventWeight);
303     }
304     }
305    
306     h_MET->Fill(met.Pt(), EventWeight);
307    
308     float dPhiMuon = 999;
309     for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
310     float tmp;
311     if (jt->Pt() > 24.) {
312     tmp = jt->DeltaPhi(met);
313     if (tmp < dPhiMuon)
314     dPhiMuon = tmp;
315     }
316     }
317     if (dPhiMuon < 10.)
318     h_minDeltaPhi_MET_Muon->Fill(dPhiMuon, EventWeight);
319 csander 1.1
320 csander 1.9 float dPhiBJet = 999;
321 csander 1.8 for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
322 csander 1.9 float tmp;
323 csander 1.8 if (it->IsBTagged()) {
324 csander 1.9 tmp = it->DeltaPhi(met);
325     if (tmp < dPhiBJet)
326     dPhiBJet = tmp;
327 csander 1.8 }
328     }
329 csander 1.9 if (dPhiBJet < 10.)
330     h_minDeltaPhi_MET_BJet->Fill(dPhiBJet, EventWeight);
331    
332     if (Muons.size() == 1) {
333     h_selectedEvents_Muon1_Pt->Fill(Muons.at(0).Pt(), EventWeight);
334     if (triggerIsoMu24)
335     h_selectedEvents_triggered_Muon1_Pt->Fill(Muons.at(0).Pt(), EventWeight);
336     }
337    
338     // h_Mbqqb_mc->Fill((hadB + hadWq + hadWqb).M(), EventWeight);
339     // h_Mbln_mc->Fill((lepB + lepWl + lepWn).M(), EventWeight);
340     //
341     // for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
342     // if (it->IsBTagged()) {
343     // for (vector<MyJet>::iterator jt = Jets.begin(); jt != Jets.end(); ++jt) {
344     // if (jt != it && !(jt->IsBTagged())) {
345     // for (vector<MyJet>::iterator kt = Jets.begin(); kt != Jets.end(); ++kt) {
346     // if (kt != it && kt != jt && !(kt->IsBTagged())) {
347     // h_Mbqqb_reco->Fill(((*it) + (*jt) + (*kt)).M(), EventWeight);
348     // }
349     // }
350     // }
351     // }
352     // }
353     // }
354     //
355     // for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
356     // if (it->IsBTagged()) {
357     // for (vector<MyMuon>::iterator jt = Muons.begin(); jt != Muons.end(); ++jt) {
358     // if ((jt->GetIsolation() / jt->Pt() < 0.1) && jt->Pt() > 24.)
359     // h_Mbln_reco->Fill(((*it) + (*jt) + met).M(), EventWeight);
360     // }
361     // }
362     // }
363 csander 1.8
364 csander 1.1 return kTRUE;
365     }
366    
367 csander 1.2 void MyAnalysis::SlaveTerminate() {
368 csander 1.1 // The SlaveTerminate() function is called after all entries or objects
369     // have been processed. When running with PROOF SlaveTerminate() is called
370     // on each slave server.
371    
372     }
373    
374 csander 1.2 void MyAnalysis::Terminate() {
375 csander 1.1 // The Terminate() function is the last function to be called during
376     // a query. It always runs on the client, it can be used to present
377     // the results graphically or save the results to file.
378    
379     }