ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/csander/HEPTutorial/MyAnalysis.C
Revision: 1.15
Committed: Wed May 23 14:48:20 2012 UTC (12 years, 11 months ago) by csander
Content type: text/plain
Branch: MAIN
Changes since 1.14: +2 -0 lines
Log Message:
added EventSelectorByList

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 csander 1.14 jet.SetJetID(Jet_ID[i]);
39     jet.SetMatchedGenJet(Jet_Px_gen[i], Jet_Py_gen[i], Jet_Pz_gen[i], Jet_E_gen[i]);
40     jet.SetFlavor(Jet_Flavor_gen[i]);
41 csander 1.3 Jets.push_back(jet);
42     }
43     Muons.clear();
44     for (int i = 0; i < NMuon; ++i) {
45     MyMuon muon(Muon_Px[i], Muon_Py[i], Muon_Pz[i], Muon_E[i]);
46     muon.SetIsolation(Muon_Iso[i]);
47 csander 1.4 muon.SetCharge(Muon_Charge[i]);
48 csander 1.14 muon.SetChargeGen(Muon_Charge_gen[i]);
49     muon.SetMatchedGenMuon(Muon_Px_gen[i], Muon_Py_gen[i], Muon_Pz_gen[i], Muon_E_gen[i]);
50 csander 1.3 Muons.push_back(muon);
51     }
52 csander 1.4 Electrons.clear();
53     for (int i = 0; i < NElectron; ++i) {
54     MyElectron electron(Electron_Px[i], Electron_Py[i], Electron_Pz[i], Electron_E[i]);
55     electron.SetIsolation(Electron_Iso[i]);
56     electron.SetCharge(Electron_Charge[i]);
57 csander 1.14 electron.SetChargeGen(Electron_Charge_gen[i]);
58     electron.SetMatchedGenElectron(Electron_Px_gen[i], Electron_Py_gen[i], Electron_Pz_gen[i], Electron_E_gen[i]);
59 csander 1.4 Electrons.push_back(electron);
60     }
61     Photons.clear();
62     for (int i = 0; i < NPhoton; ++i) {
63     MyPhoton photon(Photon_Px[i], Photon_Py[i], Photon_Pz[i], Photon_E[i]);
64     photon.SetIsolation(Photon_Iso[i]);
65 csander 1.14 photon.SetMatchedGenPhoton(Photon_Px_gen[i], Photon_Py_gen[i], Photon_Pz_gen[i], Photon_E_gen[i]);
66 csander 1.4 Photons.push_back(photon);
67     }
68 csander 1.14 Taus.clear();
69     for (int i = 0; i < NTau; ++i) {
70     MyTau tau(Tau_Px[i], Tau_Py[i], Tau_Pz[i], Tau_E[i]);
71     tau.SetTauID(Jet_ID[i]);
72     tau.SetMatchedGenJet(Tau_Px_gen[i], Tau_Py_gen[i], Tau_Pz_gen[i], Tau_E_gen[i]);
73     Taus.push_back(tau);
74     }
75 csander 1.8 met.SetXYZM(MET_px, MET_py, 0., 0.);
76 csander 1.14 genmet.SetXYZM(MET_px_gen, MET_py_gen, 0., 0.);
77 csander 1.9
78 csander 1.10 if (EventWeight < 1.e-5)
79 csander 1.9 EventWeight = 1.;
80     EventWeight *= weight_factor;
81    
82 csander 1.3 }
83    
84 csander 1.2 void MyAnalysis::Begin(TTree * /*tree*/) {
85 csander 1.1 // The Begin() function is called at the start of the query.
86     // When running with PROOF Begin() is only called on the client.
87     // The tree argument is deprecated (on PROOF 0 is passed).
88    
89     TString option = GetOption();
90    
91     }
92    
93 csander 1.2 void MyAnalysis::SlaveBegin(TTree * /*tree*/) {
94 csander 1.1 // The SlaveBegin() function is called after the Begin() function.
95     // When running with PROOF SlaveBegin() is called on each slave server.
96     // The tree argument is deprecated (on PROOF 0 is passed).
97    
98     TString option = GetOption();
99    
100 csander 1.13 h_Mmumu = new TH1F("Mmumu", "Invariant di-muon mass", 70, 0, 700);
101 csander 1.4 h_Mmumu->SetXTitle("m_{#mu#mu}");
102     h_Mmumu->Sumw2();
103    
104 csander 1.9 h_NJet = new TH1F("NJet", "Number of jets", 7, 0, 7);
105     h_NJet->SetXTitle("NJet");
106     h_NJet->Sumw2();
107    
108     h_NBJet = new TH1F("NBJet", "Number of b-jets", 7, 0, 7);
109     h_NBJet->SetXTitle("NBJet");
110     h_NBJet->Sumw2();
111    
112 csander 1.10 h_Jet1_Pt = new TH1F("Jet1Pt", "Pt of jet1", 50, 0, 250);
113     h_Jet1_Pt->SetXTitle("p_{T} [GeV]");
114 csander 1.9 h_Jet1_Pt->Sumw2();
115    
116 csander 1.10 h_Jet2_Pt = new TH1F("Jet2Pt", "Pt of jet2", 50, 0, 250);
117     h_Jet2_Pt->SetXTitle("p_{T} [GeV]");
118 csander 1.9 h_Jet2_Pt->Sumw2();
119    
120 csander 1.10 h_Jet3_Pt = new TH1F("Jet3Pt", "Pt of jet3", 50, 0, 250);
121     h_Jet3_Pt->SetXTitle("p_{T} [GeV]");
122 csander 1.9 h_Jet3_Pt->Sumw2();
123    
124 csander 1.10 h_Jet1_Eta = new TH1F("Jet1Eta", "#eta of jet1", 50, -4, 4);
125 csander 1.9 h_Jet1_Eta->SetXTitle("#eta");
126     h_Jet1_Eta->Sumw2();
127    
128 csander 1.10 h_Jet2_Eta = new TH1F("Jet2Eta", "#eta of jet2", 50, -4, 4);
129 csander 1.9 h_Jet2_Eta->SetXTitle("#eta");
130     h_Jet2_Eta->Sumw2();
131    
132 csander 1.10 h_Jet3_Eta = new TH1F("Jet3Eta", "#eta of jet3", 50, -4, 4);
133 csander 1.9 h_Jet3_Eta->SetXTitle("#eta");
134     h_Jet3_Eta->Sumw2();
135    
136 csander 1.10 h_BJet1_Pt = new TH1F("BJet1Pt", "Pt of b-jet1", 50, 0, 250);
137     h_BJet1_Pt->SetXTitle("p_{T} [GeV]");
138 csander 1.9 h_BJet1_Pt->Sumw2();
139    
140 csander 1.10 h_BJet2_Pt = new TH1F("BJet2Pt", "Pt of b-jet2", 50, 0, 250);
141     h_BJet2_Pt->SetXTitle("p_{T} [GeV]");
142 csander 1.9 h_BJet2_Pt->Sumw2();
143    
144 csander 1.10 h_BJet1_Eta = new TH1F("BJet1Eta", "#eta of b-jet1", 50, -4, 4);
145 csander 1.9 h_BJet1_Eta->SetXTitle("#eta");
146     h_BJet1_Eta->Sumw2();
147    
148 csander 1.10 h_BJet2_Eta = new TH1F("BJet2Eta", "#eta of b-jet2", 50, -4, 4);
149 csander 1.9 h_BJet2_Eta->SetXTitle("#eta");
150     h_BJet2_Eta->Sumw2();
151    
152     h_MET = new TH1F("MET", "MET", 30, 0, 300.);
153     h_MET->SetXTitle("MET");
154     h_MET->Sumw2();
155    
156 csander 1.1 }
157    
158 csander 1.2 Bool_t MyAnalysis::Process(Long64_t entry) {
159 csander 1.1 // The Process() function is called for each entry in the tree (or possibly
160     // keyed object in the case of PROOF) to be processed. The entry argument
161     // specifies which entry in the currently loaded tree is to be processed.
162     // It can be passed to either MyAnalysis::GetEntry() or TBranch::GetEntry()
163     // to read either all or the required parts of the data. When processing
164     // keyed objects with PROOF, the object is already loaded and is available
165     // via the fObject pointer.
166     //
167     // This function should contain the "body" of the analysis. It can contain
168     // simple or elaborate selection criteria, run algorithms on the data
169     // of the event and typically fill histograms.
170     //
171     // The processing can be stopped by calling Abort().
172     //
173     // Use fStatus to set the return value of TTree::Process().
174     //
175     // The return value is currently not used.
176    
177 csander 1.15 EventSelectorByList singleTau("CandidatesSingleTau.txt");
178    
179 csander 1.3 ++TotalEvents;
180    
181     GetEntry(entry);
182    
183 csander 1.4 if (TotalEvents % 10000 == 0)
184     cout << "Next event -----> " << TotalEvents << endl;
185    
186 csander 1.3 BuildEvent();
187    
188 csander 1.7 // cout << "Jets: " << endl;
189     // for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
190     // cout << "pt, eta, phi, btag: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", " << it->IsBTagged()
191     // << endl;
192     // }
193     // cout << "Muons: " << endl;
194     // for (vector<MyMuon>::iterator it = Muons.begin(); it != Muons.end(); ++it) {
195     // cout << "pt, eta, phi, iso, charge: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", "
196     // << it->GetIsolation() << ", " << it->GetCharge() << endl;
197     // }
198     // cout << "Electrons: " << endl;
199     // for (vector<MyElectron>::iterator it = Electrons.begin(); it != Electrons.end(); ++it) {
200     // cout << "pt, eta, phi, iso, charge: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", "
201     // << it->GetIsolation() << ", " << it->GetCharge() << endl;
202     // }
203     // cout << "Photons: " << endl;
204     // for (vector<MyPhoton>::iterator it = Photons.begin(); it != Photons.end(); ++it) {
205     // cout << "pt, eta, phi, iso: " << it->Pt() << ", " << it->Eta() << ", " << it->Phi() << ", " << it->GetIsolation()
206     // << endl;
207     // }
208    
209 csander 1.11 //////////////////////////////
210 csander 1.10 if (NMuon > 1) {
211     if (Muons.at(0).GetCharge() * Muons.at(1).GetCharge() < 1 && Muons.at(0).IsIsolated() && Muons.at(1).IsIsolated()) {
212     h_Mmumu->Fill((Muons.at(0) + Muons.at(1)).M(), EventWeight);
213     }
214 csander 1.5 }
215 csander 1.11 //////////////////////////////
216    
217     //////////////////////////////
218 csander 1.14 int N_BJet = 0;
219     int N_Jet = 0;
220     for (vector<MyJet>::iterator it = Jets.begin(); it != Jets.end(); ++it) {
221     ++N_Jet;
222     if (N_Jet == 1) {
223     h_Jet1_Pt->Fill(it->Pt(), EventWeight);
224     h_Jet1_Eta->Fill(it->Eta(), EventWeight);
225     } else if (N_Jet == 2) {
226     h_Jet2_Pt->Fill(it->Pt(), EventWeight);
227     h_Jet2_Eta->Fill(it->Eta(), EventWeight);
228     } else if (N_Jet == 3) {
229     h_Jet3_Pt->Fill(it->Pt(), EventWeight);
230     h_Jet3_Eta->Fill(it->Eta(), EventWeight);
231 csander 1.11 }
232 csander 1.14 if (it->IsBTagged()) {
233     ++N_BJet;
234     if (N_BJet == 1) {
235     h_BJet1_Pt->Fill(it->Pt(), EventWeight);
236     h_BJet1_Eta->Fill(it->Eta(), EventWeight);
237     } else if (N_BJet == 2) {
238     h_BJet2_Pt->Fill(it->Pt(), EventWeight);
239     h_BJet2_Eta->Fill(it->Eta(), EventWeight);
240 csander 1.11 }
241     }
242 csander 1.9 }
243 csander 1.14 h_NBJet->Fill(N_BJet, EventWeight);
244     h_NJet->Fill(N_Jet, EventWeight);
245 csander 1.9
246 csander 1.14 h_MET->Fill(met.Pt(), EventWeight);
247 csander 1.11 //////////////////////////////
248 csander 1.8
249 csander 1.1 return kTRUE;
250     }
251    
252 csander 1.2 void MyAnalysis::SlaveTerminate() {
253 csander 1.1 // The SlaveTerminate() function is called after all entries or objects
254     // have been processed. When running with PROOF SlaveTerminate() is called
255     // on each slave server.
256    
257     }
258    
259 csander 1.2 void MyAnalysis::Terminate() {
260 csander 1.1 // The Terminate() function is the last function to be called during
261     // a query. It always runs on the client, it can be used to present
262     // the results graphically or save the results to file.
263    
264     }