ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FGolf/Tools/EfficiencyModelTools.cc
Revision: 1.1
Committed: Sat Mar 17 19:53:30 2012 UTC (13 years, 1 month ago) by fgolf
Content type: text/plain
Branch: MAIN
Log Message:
useful tools for deriving and testing efficiency models

File Contents

# Content
1 #include "EfficiencyModelTools.h"
2
3 #include "CORE/CMS2.h"
4
5 #include "CORE/ssSelections.h"
6 #include "CORE/mcSelections.h"
7
8 #include "Math/VectorUtil.h"
9
10 std::vector<unsigned int> efftools::getRecoJets () {
11
12 const float mupt = 20.;
13 const float elpt = 20.;
14
15 std::vector<unsigned int> goodJets_;
16 for (unsigned int jidx = 0; jidx < cms2.pfjets_p4().size(); jidx++) {
17
18 float jet_eta = cms2.pfjets_p4().at(jidx).eta();
19 float jet_pt = cms2.pfjets_p4().at(jidx).pt() * cms2.pfjets_corL1FastL2L3().at(jidx);
20
21 if (fabs(jet_eta) > 2.5)
22 continue;
23 if (jet_pt < 40.)
24 continue;
25
26 bool jetIsLep = false;
27 for (unsigned int eidx = 0; eidx < cms2.els_p4().size(); eidx++) {
28
29 if (cms2.els_p4().at(eidx).pt() < elpt)
30 continue;
31 if (fabs(cms2.els_p4().at(eidx).eta()) > 2.4)
32 continue;
33 if (!samesign::isNumeratorLepton(11, eidx, samesign::TIGHT_DET_ISO))
34 continue;
35
36 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), cms2.els_p4().at(eidx)) < 0.4) {
37 jetIsLep = true;
38 break;
39 }
40 }
41
42 if (jetIsLep)
43 continue;
44
45 for (unsigned int midx = 0; midx < cms2.mus_p4().size(); midx++) {
46
47 if (cms2.mus_p4().at(midx).pt() < mupt)
48 continue;
49 if (fabs(cms2.mus_p4().at(midx).eta()) > 2.4)
50 continue;
51 if (!isNumeratorLepton(13, midx, samesign::TIGHT_DET_ISO))
52 continue;
53
54 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), cms2.mus_p4().at(midx)) < 0.4) {
55 jetIsLep = true;
56 break;
57 }
58 }
59
60 if (jetIsLep)
61 continue;
62
63 goodJets_.push_back(jidx);
64 }
65
66 return goodJets_;
67 }
68
69 std::vector<unsigned int> efftools::getRecoJets (struct GenParticleStruct& struct1, struct GenParticleStruct& struct2) {
70
71 const float mupt = 20.;
72 const float elpt = 20.;
73
74 std::vector<unsigned int> goodJets_;
75 for (unsigned int jidx = 0; jidx < cms2.pfjets_p4().size(); jidx++) {
76
77 float jet_eta = cms2.pfjets_p4().at(jidx).eta();
78 float jet_pt = cms2.pfjets_p4().at(jidx).pt() * cms2.pfjets_corL1FastL2L3().at(jidx);
79
80 if (fabs(jet_eta) > 2.5)
81 continue;
82 if (jet_pt < 40.)
83 continue;
84
85 std::pair<LorentzVector, int> recoLep1 = getRecoLepton(struct1);
86 if (recoLep1.second > -1) {
87 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), recoLep1.first) < 0.4)
88 continue;
89 }
90
91 std::pair<LorentzVector, int> recoLep2 = getRecoLepton(struct2);
92 if (recoLep2.second > -1) {
93 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), recoLep2.first) < 0.4)
94 continue;
95 }
96
97 bool jetIsLep = false;
98 for (unsigned int eidx = 0; eidx < cms2.els_p4().size(); eidx++) {
99
100 if (cms2.els_p4().at(eidx).pt() < elpt)
101 continue;
102 if (fabs(cms2.els_p4().at(eidx).eta()) > 2.4)
103 continue;
104 if (!samesign::isNumeratorLepton(11, eidx, samesign::TIGHT_DET_ISO))
105 continue;
106
107 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), cms2.els_p4().at(eidx)) < 0.4) {
108 jetIsLep = true;
109 break;
110 }
111 }
112
113 if (jetIsLep)
114 continue;
115
116 for (unsigned int midx = 0; midx < cms2.mus_p4().size(); midx++) {
117
118 if (cms2.mus_p4().at(midx).pt() < mupt)
119 continue;
120 if (fabs(cms2.mus_p4().at(midx).eta()) > 2.4)
121 continue;
122 if (!isNumeratorLepton(13, midx, samesign::TIGHT_DET_ISO))
123 continue;
124
125 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), cms2.mus_p4().at(midx)) < 0.4) {
126 jetIsLep = true;
127 break;
128 }
129 }
130
131 if (jetIsLep)
132 continue;
133
134 goodJets_.push_back(jidx);
135 }
136
137 return goodJets_;
138 }
139
140 std::vector<unsigned int> efftools::getRecoJets (unsigned int hidx) {
141
142 const float mupt = 20.;
143 const float elpt = 20.;
144
145 std::vector<unsigned int> goodJets_;
146 for (unsigned int jidx = 0; jidx < cms2.pfjets_p4().size(); jidx++) {
147
148 float jet_eta = cms2.pfjets_p4().at(jidx).eta();
149 float jet_pt = cms2.pfjets_p4().at(jidx).pt() * cms2.pfjets_corL1FastL2L3().at(jidx);
150
151 if (fabs(jet_eta) > 2.5)
152 continue;
153 if (jet_pt < 40.)
154 continue;
155
156 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), cms2.hyp_lt_p4().at(hidx)) < 0.4)
157 continue;
158 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), cms2.hyp_ll_p4().at(hidx)) < 0.4)
159 continue;
160
161 bool jetIsLep = false;
162 for (unsigned int eidx = 0; eidx < cms2.els_p4().size(); eidx++) {
163
164 if (cms2.els_p4().at(eidx).pt() < elpt)
165 continue;
166 if (fabs(cms2.els_p4().at(eidx).eta()) > 2.4)
167 continue;
168 if (!isNumeratorLepton(11, eidx, samesign::TIGHT_DET_ISO))
169 continue;
170
171 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), cms2.els_p4().at(eidx)) < 0.4) {
172 jetIsLep = true;
173 break;
174 }
175 }
176
177 if (jetIsLep)
178 continue;
179
180 for (unsigned int midx = 0; midx < cms2.mus_p4().size(); midx++) {
181
182 if (cms2.mus_p4().at(midx).pt() < mupt)
183 continue;
184 if (fabs(cms2.mus_p4().at(midx).eta()) > 2.4)
185 continue;
186 if (!isNumeratorLepton(13, midx, samesign::TIGHT_DET_ISO))
187 continue;
188
189 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), cms2.mus_p4().at(midx)) < 0.4) {
190 jetIsLep = true;
191 break;
192 }
193 }
194
195 if (jetIsLep)
196 continue;
197
198 goodJets_.push_back(jidx);
199 }
200
201 return goodJets_;
202 }
203
204 float efftools::getHT (std::vector<unsigned int>& jets) {
205
206 float ht_ = 0;
207 for (unsigned int jidx = 0; jidx < jets.size(); jidx++)
208 ht_ += cms2.pfjets_p4().at(jets.at(jidx)).pt() * cms2.pfjets_corL1FastL2L3().at(jets.at(jidx));
209
210 return ht_;
211 }
212
213 std::pair<LorentzVector, int> efftools::getRecoLepton (struct GenParticleStruct& struct1) {
214
215 const float mupt = 20.;
216 const float elpt = 20.;
217
218 int id = abs(struct1.id_) == 15 ? struct1.did_ : struct1.id_;
219 unsigned int idx = abs(struct1.id_) == 15 ? struct1.didx_ : struct1.idx_;
220
221 if (abs(id) == 11) {
222 for (unsigned int eidx = 0; eidx < cms2.els_p4().size(); eidx++) {
223
224 if (fabs(cms2.els_p4().at(eidx).eta()) > 2.4)
225 continue;
226 if (cms2.els_p4().at(eidx).pt() < elpt)
227 continue;
228 if (leptonIsFromW(eidx, 11, true) < 1)
229 continue;
230 if (mc3idx_eormu(11, eidx, 0.3, elpt) == idx)
231 return std::make_pair(cms2.els_p4().at(eidx), eidx);
232 }
233 }
234
235 if (abs(id) == 13) {
236 for (unsigned int midx = 0; midx < cms2.mus_p4().size(); midx++) {
237
238 if (fabs(cms2.mus_p4().at(midx).eta()) > 2.4)
239 continue;
240 if (cms2.mus_p4().at(midx).pt() < mupt)
241 continue;
242 if (leptonIsFromW(midx, 13, true) < 1)
243 continue;
244 if (mc3idx_eormu(13, midx, 0.3, mupt) == idx)
245 return std::make_pair(cms2.mus_p4().at(midx), midx);
246 }
247 }
248
249 return std::make_pair(LorentzVector(0., 0., 0., 0.), -1);
250 }
251
252 std::vector<std::pair<GenParticleStruct, GenParticleStruct> > efftools::makeGenHyps (float eta_cut, bool removeLeptonsOverlappingWithPartons) {
253
254 //
255 // loop over gen particles
256 //
257 std::vector<GenParticleStruct> gleps;
258 for (unsigned int gidx = 0; gidx < cms2.genps_p4().size(); gidx++) {
259
260 if (cms2.genps_status().at(gidx) != 3)
261 continue;
262
263 int id = cms2.genps_id().at(gidx);
264 float pt = cms2.genps_p4().at(gidx).pt();
265 float eta = cms2.genps_p4().at(gidx).eta();
266
267 if (fabs(eta) > eta_cut)
268 continue;
269
270 if (abs(id) == 11 || abs(id) == 13) {
271 if (removeLeptonsOverlappingWithPartons && leptonOverlapsWithParton(cms2.genps_p4().at(gidx)))
272 continue;
273
274 gleps.push_back(GenParticleStruct(id, gidx, pt, -999, 999, -999.));
275 }
276
277 //
278 // need to add protection here to not get more than one lepton from the same tau
279 //
280
281 if (abs(id) == 15) {
282
283 GenParticleStruct tmp_gp = GenParticleStruct(0, 0, 0, 0, 0, 0);
284 for (unsigned int didx = 0; didx < cms2.genps_lepdaughter_id().at(gidx).size(); didx++) {
285
286 int did = cms2.genps_lepdaughter_id().at(gidx).at(didx);
287 if (abs(did) != 11 && abs(did) != 13)
288 continue;
289
290 if (fabs(cms2.genps_lepdaughter_p4().at(gidx).at(didx).eta()) > eta_cut)
291 continue;
292
293 if (removeLeptonsOverlappingWithPartons && leptonOverlapsWithParton(cms2.genps_lepdaughter_p4().at(gidx).at(didx)))
294 continue;
295
296 float dpt = cms2.genps_lepdaughter_p4().at(gidx).at(didx).pt();
297 if (dpt > tmp_gp.dpt_)
298 tmp_gp = GenParticleStruct(id, gidx, pt, did, didx, dpt);
299 }
300
301 if (tmp_gp.dpt_ > 0.1)
302 gleps.push_back(tmp_gp);
303
304 } // end tau block
305
306 } // end loop over gen particles
307
308
309 //
310 // now make gen hyps
311 //
312 std::vector<std::pair<GenParticleStruct, GenParticleStruct> > glepPairs;
313 for (unsigned int idx1 = 0; idx1 < gleps.size(); idx1++) {
314 for (unsigned int idx2 = idx1 + 1; idx2 < gleps.size(); idx2++) {
315 glepPairs.push_back(make_pair(gleps.at(idx1), gleps.at(idx2)));
316 }
317 } // loop to make gen pairs
318
319 return glepPairs;
320 }
321
322 std::pair<GenParticleStruct, GenParticleStruct> efftools::getGenHyp (float pt1_cut, float pt2_cut, enum HypType hypType) {
323
324 float min_pt = std::min(pt1_cut, pt2_cut);
325 float max_pt = std::max(pt1_cut, pt2_cut);
326
327 std::vector<std::pair<GenParticleStruct, GenParticleStruct> > glepPairs = makeGenHyps();
328 unsigned int npairs = glepPairs.size();
329
330 GenParticleStruct gp = GenParticleStruct(0, -1, 0., 0, -1, 0.);
331 std::pair<GenParticleStruct, GenParticleStruct> good_gen_hyp = std::make_pair(gp, gp);
332 enum DileptonHypType good_hyp_type = DILEPTON_ALL;
333 for (unsigned int idx = 0; idx < npairs; idx++) {
334
335 GenParticleStruct gp1 = glepPairs.at(idx).first;
336 GenParticleStruct gp2 = glepPairs.at(idx).second;
337
338 float gpt1 = (abs(gp1.id_) == 15) ? gp1.dpt_ : gp1.pt_;
339 int gid1 = (abs(gp1.id_) == 15) ? gp1.did_ : gp1.id_;
340 float gpt2 = (abs(gp2.id_) == 15) ? gp2.dpt_ : gp2.pt_;
341 int gid2 = (abs(gp2.id_) == 15) ? gp2.did_ : gp2.id_;
342
343 // require SS
344 if (gid1 * gid2 < 0 && hypType == DILEPTON_SS)
345 continue;
346 else if (gid1 * gid2 > 0 && hypType == DILEPTON_OS)
347 continue;
348
349 if (std::min(gpt1, gpt2) < min_pt)
350 continue;
351
352 if (std::max(gpt1, gpt2) < max_pt)
353 continue;
354
355 enum DileptonHypType tmp_type = getHypType(gid1, gid2);
356 if (tmp_type == DILEPTON_ALL)
357 continue;
358
359 if (good_hyp_type == DILEPTON_ALL) {
360 good_hyp_type = tmp_type;
361 good_gen_hyp = glepPairs.at(idx);
362 }
363 else if (tmp_type < good_hyp_type) {
364 good_hyp_type = tmp_type;
365 good_gen_hyp = glepPairs.at(idx);
366 }
367 else if (tmp_type == good_hyp_type) {
368
369 float ggpt1 = (abs(good_gen_hyp.first.id_) == 15) ? good_gen_hyp.first.dpt_ : good_gen_hyp.first.pt_;
370 int ggid1 = (abs(good_gen_hyp.first.id_) == 15) ? good_gen_hyp.first.did_ : good_gen_hyp.first.id_;
371 float ggpt2 = (abs(good_gen_hyp.second.id_) == 15) ? good_gen_hyp.second.dpt_ : good_gen_hyp.second.pt_;
372 int ggid2 = (abs(good_gen_hyp.second.id_) == 15) ? good_gen_hyp.second.did_ : good_gen_hyp.second.id_;
373
374
375 if ( (gpt1+gpt2) > (ggpt1+ggpt2) ) {
376 good_hyp_type = tmp_type;
377 good_gen_hyp = glepPairs.at(idx);
378 }
379 }
380 } // end loop over gen hyps
381
382 return good_gen_hyp;
383 }
384
385 int efftools::getRecoJet (LorentzVector p4) {
386
387 float tmp_dr = 99.;
388 int b_index = -1;
389 for (unsigned int jidx = 0; jidx < cms2.pfjets_p4().size(); jidx++) {
390
391 if (fabs(cms2.pfjets_p4().at(jidx).eta()) > 2.5)
392 continue;
393
394 float jet_pt = cms2.pfjets_p4().at(jidx).pt() * cms2.pfjets_corL1FastL2L3().at(jidx);
395 if (jet_pt < 40.)
396 continue;
397 if (!passesPFJetID(jidx))
398 continue;
399
400 if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), p4) < tmp_dr) {
401 tmp_dr = ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), p4);
402 b_index = jidx;
403 }
404 }
405
406 return b_index;
407 }
408
409 std::vector<LorentzVector> efftools::getGenBjets (float pt_cut, float eta_cut) {
410
411 std::vector<LorentzVector> tmp_vec;
412 tmp_vec.clear();
413 for (unsigned int gidx = 0; gidx < cms2.genps_status().size(); gidx++) {
414 if (cms2.genps_status().at(gidx) != 3)
415 continue;
416 if (abs(cms2.genps_id().at(gidx)) != 5)
417 continue;
418 if (fabs(cms2.genps_p4().at(gidx).eta()) > eta_cut)
419 continue;
420 if (cms2.genps_p4().at(gidx).pt() < pt_cut)
421 continue;
422
423 tmp_vec.push_back(cms2.genps_p4().at(gidx));
424 }
425
426 return tmp_vec;
427 }
428
429 std::vector<LorentzVector> efftools::getGenJets (float pt_cut, float eta_cut) {
430
431 std::vector<LorentzVector> tmp_vec;
432 tmp_vec.clear();
433 for (unsigned int gidx = 0; gidx < cms2.genps_status().size(); gidx++) {
434 if (cms2.genps_status().at(gidx) != 3)
435 continue;
436 if ((abs(cms2.genps_id().at(gidx)) < 1 || abs(cms2.genps_id().at(gidx)) > 5) && abs(cms2.genps_id().at(gidx)) != 21)
437 continue;
438 if (fabs(cms2.genps_p4().at(gidx).eta()) > eta_cut)
439 continue;
440 if (cms2.genps_p4().at(gidx).pt() < pt_cut)
441 continue;
442
443 tmp_vec.push_back(cms2.genps_p4().at(gidx));
444 }
445
446 return tmp_vec;
447 }
448
449 float efftools::getGenHT(float pt_cut, float eta_cut) {
450
451 std::vector<LorentzVector> tmp_vec = getGenJets(pt_cut, eta_cut);
452
453 float tmp_ht = 0.;
454 for(unsigned int idx = 0; idx < tmp_vec.size(); idx++) {
455 tmp_ht += tmp_vec.at(idx).pt();
456 }
457
458 return tmp_ht;
459 }
460
461 bool efftools::leptonOverlapsWithParton(LorentzVector p4, float parton_pt, float dr) {
462
463 for (unsigned int idx = 0; idx < cms2.genps_p4().size(); idx++) {
464
465 if (cms2.genps_status().at(idx) != 3)
466 continue;
467 if ((abs(cms2.genps_id().at(idx)) < 1 || abs(cms2.genps_id().at(idx)) > 5) && abs(cms2.genps_id().at(idx)) != 21) // require a quark or gluon
468 continue;
469 if (cms2.genps_p4().at(idx).pt() < parton_pt)
470 continue;
471 if (ROOT::Math::VectorUtil::DeltaR(cms2.genps_p4().at(idx), p4) < dr)
472 return true;
473 }
474
475 return false;
476 }
477
478 enum DileptonHypType efftools::getHypType (int id1, int id2) {
479
480 if (abs(id1) != 11 && abs(id1) != 13)
481 return DILEPTON_ALL;
482 if (abs(id2) != 11 && abs(id2) != 13)
483 return DILEPTON_ALL;
484
485 if (abs(id1) == 11 && abs(id2) == 11)
486 return DILEPTON_EE;
487 else if (abs(id1) == 13 && abs(id2) == 13)
488 return DILEPTON_MUMU;
489 else
490 return DILEPTON_EMU;
491 }
492
493
494 int efftools::getGenParton (LorentzVector p4) {
495
496 float tmp_dr = 99.;
497 int index = -1;
498 for (unsigned int jidx = 0; jidx < cms2.genps_p4().size(); jidx++) {
499
500 if (fabs(cms2.genps_p4().at(jidx).eta()) > 2.5)
501 continue;
502 if (cms2.genps_p4().at(jidx).pt() < 40.)
503 continue;
504
505 if (ROOT::Math::VectorUtil::DeltaR(cms2.genps_p4().at(jidx), p4) < tmp_dr) {
506 tmp_dr = ROOT::Math::VectorUtil::DeltaR(cms2.genps_p4().at(jidx), p4);
507 index = jidx;
508 }
509 }
510
511 return index;
512 }