ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UNL/macros/MyHbbAnalyzer.C
Revision: 1.3
Committed: Mon Sep 5 18:18:12 2011 UTC (13 years, 8 months ago) by malbouis
Content type: text/plain
Branch: MAIN
Changes since 1.2: +343 -87 lines
Log Message:
update analysis code

File Contents

# User Rev Content
1 malbouis 1.1 #define MyHbbAnalyzer_cxx
2     #include "MyHbbAnalyzer.h"
3     #include <TH2.h>
4     #include <TStyle.h>
5     #include <TCanvas.h>
6 malbouis 1.3 #include <iostream>
7     #include <fstream>
8 malbouis 1.1
9     void MyHbbAnalyzer::Loop()
10     {
11     // In a ROOT session, you can do:
12     // Root > .L MyHbbAnalyzer.C
13     // Root > MyHbbAnalyzer t
14     // Root > t.GetEntry(12); // Fill t data members with entry number 12
15     // Root > t.Show(); // Show values of entry 12
16     // Root > t.Show(16); // Read and show values of entry 16
17     // Root > t.Loop(); // Loop on all entries
18     //
19    
20     // This is the loop skeleton where:
21     // jentry is the global entry number in the chain
22     // ientry is the entry number in the current Tree
23     // Note that the argument to GetEntry must be:
24     // jentry for TChain::GetEntry
25     // ientry for TTree::GetEntry and TBranch::GetEntry
26     //
27     // To read only selected branches, Insert statements like:
28     // METHOD1:
29     // fChain->SetBranchStatus("*",0); // disable all branches
30     // fChain->SetBranchStatus("branchname",1); // activate branchname
31     // METHOD2: replace line
32     // fChain->GetEntry(jentry); //read all branches
33     //by b_branchname->GetEntry(ientry); //read only this branch
34     if (fChain == 0) return;
35    
36     Long64_t nentries = fChain->GetEntriesFast();
37    
38     Long64_t nbytes = 0, nb = 0;
39    
40     ////////
41 malbouis 1.3 // open output text file
42     ////////
43     ofstream myfile;
44     myfile.open ("sync10K.txt");
45    
46     ////////
47     // get samples weight:
48     ////////
49     double weight = 1.0;
50     double ZH_weight = SetTTreeWeight("ZH");
51     double ZZ_weight = SetTTreeWeight("ZZ");
52     double DY_weight = SetTTreeWeight("DY");
53    
54     ////////
55     // initialize variables for significance
56     ////////
57    
58     // CSV vars:
59     // variables for significance:
60     // variables for significance:
61     double cut_CSV = 0, cut2_CSV = 0;
62     double cut_interval_CSV = 0.025;
63     int const max_cut_iter_CSV = 200;
64    
65     int evt_count_CSV_zh[max_cut_iter_CSV][2];
66     //initializeArr(evt_count_CSV_zh, max_cut_iter_CSV);
67     int evt_count_CSV_zz[max_cut_iter_CSV][2];
68     //initializeArr(evt_count_CSV_zz, max_cut_iter_CSV);
69     int evt_count_CSV_z[max_cut_iter_CSV][2];
70     //initializeArr(evt_count_CSV_z, max_cut_iter_CSV);
71     // 2D significance:
72     int evt_count_CSV2D_zh[max_cut_iter_CSV][max_cut_iter_CSV];
73     //initializeArr2D(evt_count_CSV2D_zh, max_cut_iter_CSV);
74     int evt_count_CSV2D_zz[max_cut_iter_CSV][max_cut_iter_CSV];
75     //initializeArr2D(evt_count_CSV2D_zz, max_cut_iter_CSV);
76     int evt_count_CSV2D_z[max_cut_iter_CSV][max_cut_iter_CSV];
77     //initializeArr2D(evt_count_CSV2D_z, max_cut_iter_CSV);
78    
79     for (int i = 0; i != max_cut_iter_CSV; i++){
80     for (int j = 0; j != max_cut_iter_CSV; j++){
81     evt_count_CSV2D_zh[i][j] = 0;
82     evt_count_CSV2D_zz[i][j] = 0;
83     evt_count_CSV2D_z[i][j] = 0;
84     if (j < 2){
85     evt_count_CSV_zh[i][j] = 0;
86     evt_count_CSV_zz[i][j] = 0;
87     evt_count_CSV_z[i][j] = 0;
88     }// j < 2
89     }// cuts
90     }// cuts
91    
92     ////////
93 malbouis 1.1 // cut and count variables
94     ////////
95     int allEvts = 0;
96     int preSelection = 0;
97     int pT_jj = 0, pT_Z = 0, CSV1 = 0, CSV2 = 0, DPhi_ZH = 0, N_aj = 0, M_jj = 0;
98 malbouis 1.3 int allEvts_zz = 0;
99     int preSelection_zz = 0;
100     int pT_jj_zz = 0, pT_Z_zz = 0, CSV1_zz = 0, CSV2_zz = 0, DPhi_ZH_zz = 0, N_aj_zz = 0, M_jj_zz = 0;
101     int allEvts_z = 0;
102     int preSelection_z = 0;
103     int pT_jj_z = 0, pT_Z_z = 0, CSV1_z = 0, CSV2_z = 0, DPhi_ZH_z = 0, N_aj_z = 0, M_jj_z = 0;
104    
105 malbouis 1.1 ////////
106     // loop on entries
107     ////////
108     pct = 0;
109     total_entries = fChain->GetEntries();
110    
111     //nentries = 1000;
112     for (Long64_t jentry=0; jentry<nentries;jentry++) {
113     Long64_t ientry = LoadTree(jentry);
114     if (ientry < 0) break;
115     nb = fChain->GetEntry(jentry); nbytes += nb;
116     // if (Cut(ientry) < 0) continue;
117    
118     // jobProgress:
119     jobProgress(jentry);
120    
121 malbouis 1.3 // set weights
122     if (isZH()) weight = ZH_weight;
123     if (isZZ()) weight = ZZ_weight;
124     if (isDY()) weight = DY_weight;
125    
126     //cout << "sample ZH: " << isZH() << " DY: " << isDY() << " ZZ: " << isZZ() << endl;
127     //cout << "weight = " << weight << endl;
128    
129 malbouis 1.1 // count events:
130 malbouis 1.3 if (isZH()) allEvts++;
131     if (isZZ()) allEvts_zz++;
132     if (isDY()) allEvts_z++;
133     basicMuonDistributions("allEvts", weight);
134     basicJetDistributions("allEvts", weight);
135 malbouis 1.1
136     // pre-selection cuts:
137     if (!preSelect()) continue;
138 malbouis 1.3 if (isZH()) preSelection++;
139     if (isZZ()) preSelection_zz++;
140     if (isDY()) preSelection_z++;
141     basicMuonDistributions("preSel", weight);
142     basicJetDistributions("preSel", weight);
143     myfile << "passed pre-selection:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << "\n";
144 malbouis 1.1
145     // pt_jj cut:
146     if (!Select_pTjj()) continue;
147 malbouis 1.3 if (isZH()) pT_jj++;
148     if (isZZ()) pT_jj_zz++;
149     if (isDY()) pT_jj_z++;
150     basicMuonDistributions("pT_jj", weight);
151     basicJetDistributions("pT_jj", weight);
152     myfile << "passed pT(jj):" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
153    
154 malbouis 1.1
155     // pT_Z cut:
156     if (!Select_pTZ()) continue;
157 malbouis 1.3 if (isZH()) pT_Z++;
158     if (isZZ()) pT_Z_zz++;
159     if (isDY()) pT_Z_z++;
160     basicMuonDistributions("pT_Z", weight);
161     basicJetDistributions("pT_Z", weight);
162     myfile << "passed pTZ:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
163    
164     ////////
165     // look at CSV1 optimization:
166     // fill arrays with discriminators after cuts
167     ////////
168     int Jet1Idx, Jet2Idx;
169     jetIndices(Jet1Idx, Jet2Idx, bDisc_CSV);
170    
171     for(int icut = 0; icut != max_cut_iter_CSV; icut++){
172     cut_CSV = icut*cut_interval_CSV;
173     if (bDisc_CSV[Jet1Idx] > cut_CSV){
174     if (isZH()) evt_count_CSV_zh[icut][0] = evt_count_CSV_zh[icut][0] + 1;
175     if (isZZ()) evt_count_CSV_zz[icut][0] = evt_count_CSV_zz[icut][0] + 1;
176     if (isDY()) evt_count_CSV_z[icut][0] = evt_count_CSV_z[icut][0] + 1;
177     }
178     }// cuts
179    
180 malbouis 1.1 // CSV1 cut:
181     if (!Select_CSV1()) continue;
182 malbouis 1.3 if (isZH()) CSV1++;
183     if (isZZ()) CSV1_zz++;
184     if (isDY()) CSV1_z++;
185     basicMuonDistributions("CSV1", weight);
186     basicJetDistributions("CSV1", weight);
187     myfile << "passed CSV1:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
188    
189     ////////
190     // look at CSV2 optimization, AFTER CSV1 cut!
191     // fill arrays with discriminators after cuts
192     ////////
193    
194     for(int icut = 0; icut != max_cut_iter_CSV; icut++){
195     cut_CSV = icut*cut_interval_CSV;
196     if (bDisc_CSV[Jet2Idx] > cut_CSV){
197     if (isZH()) evt_count_CSV_zh[icut][1] = evt_count_CSV_zh[icut][1] + 1;
198     if (isZZ()) evt_count_CSV_zz[icut][1] = evt_count_CSV_zz[icut][1] + 1;
199     if (isDY()) evt_count_CSV_z[icut][1] = evt_count_CSV_z[icut][1] + 1;
200     }
201     }// cuts
202    
203 malbouis 1.1
204     // CSV2 cut:
205     if (!Select_CSV2()) continue;
206 malbouis 1.3 if (isZH()) CSV2++;
207     if (isZZ()) CSV2_zz++;
208     if (isDY()) CSV2_z++;
209     basicMuonDistributions("CSV2", weight);
210     basicJetDistributions("CSV2", weight);
211     myfile << "passed CSV2:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
212    
213 malbouis 1.1 // DPhi cut:
214     if (!Select_DPhi()) continue;
215 malbouis 1.3 if (isZH()) DPhi_ZH++;
216     if (isZZ()) DPhi_ZH_zz++;
217     if (isDY()) DPhi_ZH_z++;
218     basicMuonDistributions("DPhi", weight);
219     basicJetDistributions("DPhi", weight);
220     myfile << "passed DPhi:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
221    
222 malbouis 1.1
223     // Naj cut:
224     if (!Select_Naj()) continue;
225 malbouis 1.3 if (isZH()) N_aj++;
226     if (isZZ()) N_aj_zz++;
227     if (isDY()) N_aj_z++;
228     basicMuonDistributions("Naj", weight);
229     basicJetDistributions("Naj", weight);
230     myfile << "passed Naj:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
231    
232 malbouis 1.1
233     // Mjj cut:
234     if (!Select_Mjj()) continue;
235 malbouis 1.3 if (isZH()) M_jj++;
236     if (isZZ()) M_jj_zz++;
237     if (isDY()) M_jj_z++;
238     basicMuonDistributions("Mjj", weight);
239     basicJetDistributions("Mjj", weight);
240     myfile << "passed Mjj:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
241    
242 malbouis 1.1
243     }// loop on entries
244    
245 malbouis 1.3 //////
246     // write significance histos
247     ////////
248     int nbins = max_cut_iter_CSV;
249     double binMax = (max_cut_iter_CSV*cut_interval_CSV)+cut_interval_CSV;
250     double binMin = cut_interval_CSV;
251    
252     // fill cut histos:
253     for (int i = 0; i!=max_cut_iter_CSV; i++){
254     for (int ijet = 0; ijet != 2; ijet++){
255     string zh_csv_histo = Form("ZH_CSV_jet%i", ijet);
256     fillhisto(zh_csv_histo, i*cut_interval_CSV, ZH_weight*evt_count_CSV_zh[i][ijet], "", nbins, binMin, binMax);
257     string zz_csv_histo = Form("ZZ_CSV_jet%i", ijet);
258     fillhisto(zz_csv_histo, i*cut_interval_CSV, ZZ_weight*evt_count_CSV_zz[i][ijet], "", nbins, binMin, binMax);
259     string z_csv_histo = Form("DY_CSV_jet%i", ijet);
260     fillhisto(z_csv_histo, i*cut_interval_CSV, DY_weight*evt_count_CSV_z[i][ijet], "", nbins, binMin, binMax);
261     string sqrtBackground_CSV = Form("sqrtBackground_CSV_jet%i", ijet);
262     fillhisto(sqrtBackground_CSV, i*cut_interval_CSV, sqrt((DY_weight*evt_count_CSV_z[i][ijet])+(ZZ_weight*evt_count_CSV_zz[i][ijet])), "", nbins, binMin, binMax);
263     }// jets
264     }// cuts
265    
266    
267 malbouis 1.1 ////////
268     // write histos to file:
269     ////////
270 malbouis 1.3 TFile *f = new TFile("hbb_CSV.root","RECREATE");
271 malbouis 1.1 for (map<std::string,TH1*>::iterator it=histmap.begin(); it!=histmap.end();it++) {
272     (*it).second->Write();
273     }
274 malbouis 1.3 for (map<string,TH2*>::iterator it2=bidimhistmap.begin(); it2!=bidimhistmap.end();it2++) {
275     (*it2).second->Write();
276     delete (*it2).second;
277     }
278    
279 malbouis 1.1 f->Write();
280     f->Close();
281 malbouis 1.3
282     // close output text file
283     myfile.close();
284 malbouis 1.1
285     ////////
286     // print stats
287     ////////
288 malbouis 1.3 cout << endl;
289     cout << endl;
290     cout << "**** stats (no normalization) ****" << endl;
291     cout << "allEvts: " << allEvts << endl;
292 malbouis 1.1 cout << "preSelection: " << preSelection << endl;
293     cout << "pT_jj: " << pT_jj << endl;
294     cout << "pT_Z: " << pT_Z << endl;
295     cout << "CSV1: " << CSV1 << endl;
296     cout << "CSV2: " << CSV2 << endl;
297     cout << "DPhi: " << DPhi_ZH << endl;
298     cout << "Naj: " << N_aj << endl;
299     cout << "Mjj: " << M_jj << endl;
300 malbouis 1.3 cout << endl;
301     cout << endl;
302     cout << "**** stats normalized to 10 fb-1 ****" << endl;
303     cout << "allEvts: " << allEvts*ZH_weight << ", ZZ allEvts: " << allEvts_zz*ZZ_weight << ", DY allEvts: " << allEvts_z*DY_weight << endl;
304     cout << "preSelection: " << preSelection*ZH_weight << ", ZZ preSelection: " << preSelection_zz*ZZ_weight << ", DY preSelection: " << preSelection_z*DY_weight << endl;
305     cout << "pT_jj: " << pT_jj*ZH_weight << ", ZZ pT_jj: " << pT_jj_zz*ZZ_weight << ", DY pT_jj: " << pT_jj_z*DY_weight << endl;
306     cout << "pT_Z: " << pT_Z*ZH_weight << ", ZZ pT_Z: " << pT_Z_zz*ZZ_weight << ", DY pT_Z: " << pT_Z_z*DY_weight << endl;
307     cout << "CSV1: " << CSV1*ZH_weight << ", ZZ CSV1: " << CSV1_zz*ZZ_weight << ", DY CSV1: " << CSV1_z*DY_weight << endl;
308     cout << "CSV2: " << CSV2*ZH_weight << ", ZZ CSV2: " << CSV2_zz*ZZ_weight << ", DY CSV2: " << CSV2_z*DY_weight << endl;
309     cout << "DPhi: " << DPhi_ZH*ZH_weight << ", ZZ DPhi: " << DPhi_ZH_zz*ZZ_weight << ", DY DPhi: " << DPhi_ZH_z*DY_weight << endl;
310     cout << "Naj: " << N_aj*ZH_weight << ", ZZ Naj: " << N_aj_zz*ZZ_weight << ", DY Naj: " << N_aj_z*DY_weight << endl;
311     cout << "Mjj: " << M_jj*ZH_weight << ", ZZ Mjj: " << M_jj_zz*ZZ_weight << ", DY Mjj: " << M_jj_z*DY_weight << endl;
312     cout << endl;
313     cout << "signal eff. PRE-SELECT: " << (double)preSelection/allEvts << ", ZZ eff. PRE-SELECT: " << (double)preSelection_zz/allEvts_zz << ", DY eff. PRE-SELECT: " << (double)preSelection_z/allEvts_z << endl;
314     cout << "signal eff. pT_jj: " << (double)pT_jj/preSelection << ", ZZ eff. pT_jj: " << (double)pT_jj_zz/preSelection_zz << ", DY eff. pT_jj: " << (double)pT_jj_z/preSelection_z << endl;
315     cout << "signal eff. pT_Z: " << (double)pT_Z/pT_jj << ", ZZ signal eff. pT_Z: " << (double)pT_Z_zz/pT_jj_zz << ", DY signal eff. pT_Z: " << (double)pT_Z_z/pT_jj_z << endl;
316     cout << "signal eff. CSV1: " << (double)CSV1/pT_Z << ", ZZ eff. CSV1: " << (double)CSV1_zz/pT_Z_zz << ", DY eff. CSV1: " << (double)CSV1_z/pT_Z_z << endl;
317     cout << "signal eff. CSV2: " << (double)CSV2/CSV1 << ", ZZ eff. CSV2: " << (double)CSV2_zz/CSV1_zz << ", DY eff. CSV2: " << (double)CSV2_z/CSV1_z << endl;
318     cout << "signal eff. DPhi: " << (double)DPhi_ZH/CSV2 << ", ZZ eff. DPhi: " << (double)DPhi_ZH_zz/CSV2_zz << ", DY eff. DPhi: " << (double)DPhi_ZH_z/CSV2_z << endl;
319     cout << "signal eff. Naj: " << (double)N_aj/DPhi_ZH << ", ZZ eff. Naj: " << (double)N_aj_zz/DPhi_ZH_zz << ", DY eff. Naj: " << (double)N_aj_z/DPhi_ZH_z << endl;
320     cout << "signal eff. Mjj: " << (double)M_jj/N_aj << ", ZZ eff. Mjj: " << (double)M_jj_zz/N_aj_zz << ", DY eff. Mjj: " << (double)M_jj_z/N_aj_z << endl;
321     cout << endl;
322     cout << "**** signal only ****" << endl;
323     cout << "signal eff. PRE-SELECT: " << (double)preSelection*100/allEvts << " +- " << relativeError((double)preSelection, (double)allEvts)*100 << endl;
324     cout << "signal eff. pT_jj: " << (double)pT_jj*100/preSelection << " +- " << relativeError((double)pT_jj, (double)preSelection)*100 << endl;
325     cout << "signal eff. pT_Z: " << (double)pT_Z*100/pT_jj << " +- " << relativeError((double)pT_Z, (double)pT_jj)*100 << endl;
326     cout << "signal eff. CSV1: " << (double)CSV1*100/pT_Z << " +- " << relativeError((double)CSV1, (double)pT_Z)*100 << endl;
327     cout << "signal eff. CSV2: " << (double)CSV2*100/CSV1 << " +- " << relativeError((double)CSV2, (double)CSV1)*100 << endl;
328     cout << "signal eff. DPhi: " << (double)DPhi_ZH*100/CSV2 << " +- " << relativeError((double)DPhi_ZH, (double)CSV2)*100 << endl;
329     cout << "signal eff. Naj: " << (double)N_aj*100/DPhi_ZH << " +- " << relativeError((double)N_aj, (double)DPhi_ZH)*100 << endl;
330     cout << "signal eff. Mjj: " << (double)M_jj*100/N_aj << " +- " << relativeError((double)M_jj, (double)N_aj)*100 << endl;
331     cout << "Total signal eff.: " << (double)M_jj*100/allEvts << " +- " << relativeError((double)M_jj, (double)allEvts)*100 << endl;
332     cout << endl;
333     cout << "**** TOTAL SIGNIFICANCE, 95 < M_jj < 125 ****" << endl;
334     cout << "S = S/sqrt{B} = " << ((double)M_jj*ZH_weight)/sqrt((M_jj_zz*ZZ_weight) + (M_jj_z*DY_weight)) << endl;
335     cout << endl;
336     cout << "**** TOTAL SIGNIFICANCE before Mjj cut ****" << endl;
337     cout << "S = S/sqrt{B} = " << ((double)N_aj*ZH_weight)/sqrt((N_aj_zz*ZZ_weight) + (N_aj_z*DY_weight)) << endl;
338    
339 malbouis 1.1
340     }// Loop
341     bool MyHbbAnalyzer::preSelect(){
342     bool passedPreSel = true;
343 malbouis 1.3 // trigger requirement:
344     if (!HLT_IsoMu17_v5_trig) return false;
345 malbouis 1.1 if (nMuons < 2) passedPreSel = false;
346     if (nJets < 2) passedPreSel = false;
347     return passedPreSel;
348     }// preSelect
349    
350     bool MyHbbAnalyzer::Select_pTjj(){
351     bool passedSel_pTjj = true;
352 malbouis 1.3 double pTjj = computePtjj();
353     if (pTjj <= 100) passedSel_pTjj = false;
354 malbouis 1.1 return passedSel_pTjj;
355     }// pT_jj cut
356    
357     bool MyHbbAnalyzer::Select_pTZ(){
358     bool passed_pTZ = true;
359 malbouis 1.3 double pTZ = computePtZ();
360     if (pTZ <= 100) passed_pTZ = false;
361 malbouis 1.1 return passed_pTZ;
362     }// pT_Z cut
363    
364     bool MyHbbAnalyzer::Select_CSV1(){
365 malbouis 1.2 // max CSV tagged jet should satisfy the tight Operating Point (CSVT > 0.898)
366 malbouis 1.1 bool passedCSV1 = true;
367 malbouis 1.2 int firstJetIdx, secondJetIdx;
368 malbouis 1.3 jetIndices(firstJetIdx, secondJetIdx, bDisc_CSV);
369 malbouis 1.1 double bDiscMax;
370 malbouis 1.2 bDiscMax = bDisc_CSV[firstJetIdx] > bDisc_CSV[secondJetIdx] ? bDisc_CSV[firstJetIdx]:bDisc_CSV[secondJetIdx];
371 malbouis 1.1 if (bDiscMax <= 0.898) passedCSV1 = false;
372 malbouis 1.3 //if (bDiscMax <= 0.866) passedCSV1 = false;
373 malbouis 1.1 return passedCSV1;
374     }// CSV1
375    
376     bool MyHbbAnalyzer::Select_CSV2(){
377 malbouis 1.2 // min CSV tagged jet should satisfy CSV > 0.5
378 malbouis 1.1 bool passedCSV2 = true;
379 malbouis 1.2 int firstJetIdx, secondJetIdx;
380 malbouis 1.3 jetIndices(firstJetIdx, secondJetIdx, bDisc_CSV);
381 malbouis 1.1 double bDiscMin;
382 malbouis 1.2 bDiscMin = bDisc_CSV[firstJetIdx] < bDisc_CSV[secondJetIdx] ? bDisc_CSV[firstJetIdx]:bDisc_CSV[secondJetIdx];
383 malbouis 1.1 if (bDiscMin <= 0.5) passedCSV2 = false;
384     return passedCSV2;
385     }// CSV2
386    
387     bool MyHbbAnalyzer::Select_DPhi(){
388     bool passedDPhi = true;
389 malbouis 1.3 double DPhi = computeDPhi();
390     if (fabs(DPhi) <= 2.90) passedDPhi = false;
391 malbouis 1.1 return passedDPhi;
392     }// DPhi
393    
394     bool MyHbbAnalyzer::Select_Naj(){
395     bool passedNaj = true;
396     if (nJets > 3) passedNaj = false;
397     return passedNaj;
398     }// Naj
399    
400     bool MyHbbAnalyzer::Select_Mjj(){
401     bool passedMjj = true;
402 malbouis 1.3 double Mjj = computeMjj();
403     if (Mjj < 95 || Mjj > 125) passedMjj = false;
404 malbouis 1.1 return passedMjj;
405     }
406 malbouis 1.3 void MyHbbAnalyzer::jetIndices(int &index1, int &index2, double bDisc[5]){
407 malbouis 1.2 int fJetIdx = -1, sJetIdx = -1;
408     double highestBdiscSum = 0;
409     for (int i = 0; i != nJets; ++i){
410 malbouis 1.3 if (bDisc[i] < 0) continue;
411 malbouis 1.2 for (int j = i+1; j != nJets; ++j){
412 malbouis 1.3 if (bDisc[j] < 0) continue;
413     if (bDisc[i] + bDisc[j] > highestBdiscSum){
414     highestBdiscSum = bDisc[i] + bDisc[j];
415 malbouis 1.2 fJetIdx = i;
416     sJetIdx = j;
417     }// if greater sum
418     }// second loop
419     }// 1st loop
420     swap(fJetIdx, index1);
421     swap(sJetIdx, index2);
422     }// jetIndices
423 malbouis 1.1
424 malbouis 1.3 void MyHbbAnalyzer::initializeArr(int evt_count[10000][2], int const max_cut_iter){
425     for (int i = 0; i != max_cut_iter; i++){
426     for (int ijet = 0; ijet != 2; ijet++){
427     evt_count[i][ijet] = 0;
428     }// jets
429     }// cuts
430     }// initializeArr: end function
431    
432     void MyHbbAnalyzer::initializeArr2D(int evt_count[100][100], int const max_cut_iter){
433     for (int i = 0; i != max_cut_iter; i++){
434     for (int j = 0; j != max_cut_iter; j++){
435     evt_count[i][j] = 0;
436     }// cuts
437     }// cuts
438    
439     }// initializeArr2D: end function
440    
441     void MyHbbAnalyzer::basicMuonDistributions(string cut, double ph_weight){
442 malbouis 1.1 // var definition
443     string ph_process;
444    
445     if (isZH()) ph_process = "ZH";
446     if (isZZ()) ph_process = "ZZ";
447     if (isDY()) ph_process = "DY";
448    
449     for (int i=0; i != nMuons; ++i) {
450     string muonpT = "muonPt_" + ph_process + cut;
451     string muoneta = "muonEta_" + ph_process + cut;
452     string muonphi = "muonPhi_" + ph_process + cut;
453    
454     //cout << muonpT << endl;
455    
456 malbouis 1.3 fillhisto(muonpT, muonPt[i], ph_weight, "#mu pT", 400, 0, 400);
457     fillhisto(muoneta, muonEta[i], ph_weight, "#mu #eta", 100, -5, 5);
458     fillhisto(muonphi, muonPhi[i], ph_weight, "#mu #phi", 80, -4, 4);
459 malbouis 1.1 }
460     }// basicMuonDistributions
461    
462 malbouis 1.3 void MyHbbAnalyzer::basicJetDistributions(string cut, double ph_weight){
463 malbouis 1.1 // var definition
464     string ph_process;
465    
466 malbouis 1.3 if (isZH()){ ph_process = "ZH";}
467     if (isZZ()){ ph_process = "ZZ";}
468     if (isDY()){ ph_process = "DY";}
469 malbouis 1.1
470     for (int i=0; i != nJets; ++i) {
471     string jetpT = "jetPt_" + ph_process + cut;
472     string jeteta = "jetEta_" + ph_process + cut;
473     string jetphi = "jetPhi_" + ph_process + cut;
474    
475     //cout << muonpT << endl;
476    
477 malbouis 1.3 fillhisto(jetpT, jetPt[i], ph_weight, "jet pT", 400, 0, 400);
478     fillhisto(jeteta, jetEta[i], ph_weight, "jet #eta", 100, -5, 5);
479     fillhisto(jetphi, jetPhi[i], ph_weight, "jet #phi", 80, -4, 4);
480 malbouis 1.1 }
481 malbouis 1.3
482     // Mjj
483     double Mjj = computeMjj();
484     string mjj = "Mjj_" + ph_process + cut;
485     fillhisto(mjj, Mjj, ph_weight, "", 1000, 0, 1000);
486 malbouis 1.1 }// basicJetDistributions
487 malbouis 1.3
488     double MyHbbAnalyzer::computeMjj(){
489     int firstJetIdx, secondJetIdx;
490     jetIndices(firstJetIdx, secondJetIdx, bDisc_CSV);
491    
492     TLorentzVector firstJet_p4 = TLorentzVector(jetPx[firstJetIdx], jetPy[firstJetIdx], jetPz[firstJetIdx], jetP[firstJetIdx]);
493     TLorentzVector secondJet_p4 = TLorentzVector(jetPx[secondJetIdx], jetPy[secondJetIdx], jetPz[secondJetIdx], jetP[secondJetIdx]);
494    
495     TLorentzVector jj = firstJet_p4 + secondJet_p4;
496    
497     return jj.M();
498     }// Mjj
499    
500     double MyHbbAnalyzer::computePtjj(){
501     double ptjj = 0.;
502     // get jet indices that maximize the b-tag discriminant
503     int firstJetIdx, secondJetIdx;
504     jetIndices(firstJetIdx, secondJetIdx, bDisc_CSV);
505     // reject event with less than 2 CSV tags
506     if (firstJetIdx == -1 || secondJetIdx == -1) return false;
507     // get the dijet system
508     TLorentzVector firstJet_p4 = TLorentzVector(jetPx[firstJetIdx], jetPy[firstJetIdx], jetPz[firstJetIdx], jetP[firstJetIdx]);
509     TLorentzVector secondJet_p4 = TLorentzVector(jetPx[secondJetIdx], jetPy[secondJetIdx], jetPz[secondJetIdx], jetP[secondJetIdx]);
510     TLorentzVector jj = firstJet_p4 + secondJet_p4;
511     ptjj = jj.Pt();
512    
513     return ptjj;
514     }// Ptjj
515    
516     double MyHbbAnalyzer::computePtZ(){
517     double ptz = 0.;
518     TLorentzVector firstMu_p4 = TLorentzVector(muonPx[0], muonPy[0], muonPz[0], muonP[0]);
519     TLorentzVector secondMu_p4 = TLorentzVector(muonPx[1], muonPy[1], muonPz[1], muonP[1]);
520     TLorentzVector Z = firstMu_p4 + secondMu_p4;
521     ptz = Z.Pt();
522    
523     return ptz;
524     }// PtZ
525    
526     double MyHbbAnalyzer::computeDPhi(){
527    
528     double dphi = 99.;
529     TLorentzVector firstMu_p4 = TLorentzVector(muonPx[0], muonPy[0], muonPz[0], muonP[0]);
530     TLorentzVector secondMu_p4 = TLorentzVector(muonPx[1], muonPy[1], muonPz[1], muonP[1]);
531     TLorentzVector Z = firstMu_p4 + secondMu_p4;
532     int firstJetIdx, secondJetIdx;
533     jetIndices(firstJetIdx, secondJetIdx, bDisc_CSV);
534     TLorentzVector firstJet_p4 = TLorentzVector(jetPx[firstJetIdx], jetPy[firstJetIdx], jetPz[firstJetIdx], jetP[firstJetIdx]);
535     TLorentzVector secondJet_p4 = TLorentzVector(jetPx[secondJetIdx], jetPy[secondJetIdx], jetPz[secondJetIdx], jetP[secondJetIdx]);
536     TLorentzVector jj = firstJet_p4 + secondJet_p4;
537     dphi = Z.DeltaPhi(jj);
538    
539     return dphi;
540    
541     }// dPhi
542    
543    
544    
545