ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FGolf/Tools/EfficiencyModelTools.cc
Revision: 1.2
Committed: Mon May 28 22:11:53 2012 UTC (12 years, 11 months ago) by fgolf
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +23 -25 lines
Log Message:
update for new gen particle struct; change btag eta cut from 2.5 to 2.4; update for 2012 lepton selections

File Contents

# User Rev Content
1 fgolf 1.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 fgolf 1.2 if (fabs(jet_eta) > 2.4)
22 fgolf 1.1 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 fgolf 1.2 if (!samesign::isNumeratorLepton(11, eidx))
34 fgolf 1.1 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 fgolf 1.2 if (!samesign::isNumeratorLepton(13, midx))
52 fgolf 1.1 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 fgolf 1.2 if (fabs(jet_eta) > 2.4)
81 fgolf 1.1 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 fgolf 1.2 if (!samesign::isNumeratorLepton(11, eidx))
105 fgolf 1.1 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 fgolf 1.2 if (!samesign::isNumeratorLepton(13, midx))
123 fgolf 1.1 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 fgolf 1.2 if (fabs(jet_eta) > 2.4)
152 fgolf 1.1 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 fgolf 1.2 if (!samesign::isNumeratorLepton(11, eidx))
169 fgolf 1.1 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 fgolf 1.2 if (!samesign::isNumeratorLepton(13, midx))
187 fgolf 1.1 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 fgolf 1.2 int idx = abs(struct1.id_) == 15 ? struct1.didx_ : struct1.idx_;
220 fgolf 1.1
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 fgolf 1.2 gleps.push_back(GenParticleStruct(id, gidx, pt, eta, -999, 999, -999., -999.));
275 fgolf 1.1 }
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 fgolf 1.2 GenParticleStruct tmp_gp = GenParticleStruct(0, 0, 0, 0, 0, 0, 0, 0);
284 fgolf 1.1 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 fgolf 1.2 float dpt = cms2.genps_lepdaughter_p4().at(gidx).at(didx).pt();
297     float deta = cms2.genps_lepdaughter_p4().at(gidx).at(didx).eta();
298 fgolf 1.1 if (dpt > tmp_gp.dpt_)
299 fgolf 1.2 tmp_gp = GenParticleStruct(id, gidx, pt, eta, did, didx, dpt, deta);
300 fgolf 1.1 }
301    
302     if (tmp_gp.dpt_ > 0.1)
303     gleps.push_back(tmp_gp);
304    
305     } // end tau block
306    
307     } // end loop over gen particles
308    
309    
310     //
311     // now make gen hyps
312     //
313     std::vector<std::pair<GenParticleStruct, GenParticleStruct> > glepPairs;
314     for (unsigned int idx1 = 0; idx1 < gleps.size(); idx1++) {
315     for (unsigned int idx2 = idx1 + 1; idx2 < gleps.size(); idx2++) {
316     glepPairs.push_back(make_pair(gleps.at(idx1), gleps.at(idx2)));
317     }
318     } // loop to make gen pairs
319    
320     return glepPairs;
321     }
322    
323     std::pair<GenParticleStruct, GenParticleStruct> efftools::getGenHyp (float pt1_cut, float pt2_cut, enum HypType hypType) {
324    
325     float min_pt = std::min(pt1_cut, pt2_cut);
326     float max_pt = std::max(pt1_cut, pt2_cut);
327    
328     std::vector<std::pair<GenParticleStruct, GenParticleStruct> > glepPairs = makeGenHyps();
329     unsigned int npairs = glepPairs.size();
330    
331 fgolf 1.2 GenParticleStruct gp = GenParticleStruct(0, -1, 0., 0., 0, -1, 0., 0.);
332 fgolf 1.1 std::pair<GenParticleStruct, GenParticleStruct> good_gen_hyp = std::make_pair(gp, gp);
333     enum DileptonHypType good_hyp_type = DILEPTON_ALL;
334     for (unsigned int idx = 0; idx < npairs; idx++) {
335    
336     GenParticleStruct gp1 = glepPairs.at(idx).first;
337     GenParticleStruct gp2 = glepPairs.at(idx).second;
338    
339 fgolf 1.2 float gpt1 = (abs(gp1.id_) == 15) ? gp1.dpt_ : gp1.pt_;
340     int gid1 = (abs(gp1.id_) == 15) ? gp1.did_ : gp1.id_;
341     float gpt2 = (abs(gp2.id_) == 15) ? gp2.dpt_ : gp2.pt_;
342     int gid2 = (abs(gp2.id_) == 15) ? gp2.did_ : gp2.id_;
343 fgolf 1.1
344     // require SS
345     if (gid1 * gid2 < 0 && hypType == DILEPTON_SS)
346     continue;
347     else if (gid1 * gid2 > 0 && hypType == DILEPTON_OS)
348     continue;
349    
350     if (std::min(gpt1, gpt2) < min_pt)
351     continue;
352    
353     if (std::max(gpt1, gpt2) < max_pt)
354     continue;
355    
356     enum DileptonHypType tmp_type = getHypType(gid1, gid2);
357     if (tmp_type == DILEPTON_ALL)
358     continue;
359    
360     if (good_hyp_type == DILEPTON_ALL) {
361     good_hyp_type = tmp_type;
362     good_gen_hyp = glepPairs.at(idx);
363     }
364     else if (tmp_type < good_hyp_type) {
365     good_hyp_type = tmp_type;
366     good_gen_hyp = glepPairs.at(idx);
367     }
368     else if (tmp_type == good_hyp_type) {
369    
370     float ggpt1 = (abs(good_gen_hyp.first.id_) == 15) ? good_gen_hyp.first.dpt_ : good_gen_hyp.first.pt_;
371 fgolf 1.2 float ggpt2 = (abs(good_gen_hyp.second.id_) == 15) ? good_gen_hyp.second.dpt_ : good_gen_hyp.second.pt_;
372 fgolf 1.1
373     if ( (gpt1+gpt2) > (ggpt1+ggpt2) ) {
374     good_hyp_type = tmp_type;
375     good_gen_hyp = glepPairs.at(idx);
376     }
377     }
378     } // end loop over gen hyps
379    
380     return good_gen_hyp;
381     }
382    
383     int efftools::getRecoJet (LorentzVector p4) {
384    
385     float tmp_dr = 99.;
386     int b_index = -1;
387     for (unsigned int jidx = 0; jidx < cms2.pfjets_p4().size(); jidx++) {
388    
389 fgolf 1.2 if (fabs(cms2.pfjets_p4().at(jidx).eta()) > 2.4)
390 fgolf 1.1 continue;
391    
392     float jet_pt = cms2.pfjets_p4().at(jidx).pt() * cms2.pfjets_corL1FastL2L3().at(jidx);
393     if (jet_pt < 40.)
394     continue;
395     if (!passesPFJetID(jidx))
396     continue;
397    
398     if (ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), p4) < tmp_dr) {
399     tmp_dr = ROOT::Math::VectorUtil::DeltaR(cms2.pfjets_p4().at(jidx), p4);
400     b_index = jidx;
401     }
402     }
403    
404     return b_index;
405     }
406    
407     std::vector<LorentzVector> efftools::getGenBjets (float pt_cut, float eta_cut) {
408    
409     std::vector<LorentzVector> tmp_vec;
410     tmp_vec.clear();
411     for (unsigned int gidx = 0; gidx < cms2.genps_status().size(); gidx++) {
412     if (cms2.genps_status().at(gidx) != 3)
413     continue;
414     if (abs(cms2.genps_id().at(gidx)) != 5)
415     continue;
416     if (fabs(cms2.genps_p4().at(gidx).eta()) > eta_cut)
417     continue;
418     if (cms2.genps_p4().at(gidx).pt() < pt_cut)
419     continue;
420    
421     tmp_vec.push_back(cms2.genps_p4().at(gidx));
422     }
423    
424     return tmp_vec;
425     }
426    
427     std::vector<LorentzVector> efftools::getGenJets (float pt_cut, float eta_cut) {
428    
429     std::vector<LorentzVector> tmp_vec;
430     tmp_vec.clear();
431     for (unsigned int gidx = 0; gidx < cms2.genps_status().size(); gidx++) {
432     if (cms2.genps_status().at(gidx) != 3)
433     continue;
434     if ((abs(cms2.genps_id().at(gidx)) < 1 || abs(cms2.genps_id().at(gidx)) > 5) && abs(cms2.genps_id().at(gidx)) != 21)
435     continue;
436     if (fabs(cms2.genps_p4().at(gidx).eta()) > eta_cut)
437     continue;
438     if (cms2.genps_p4().at(gidx).pt() < pt_cut)
439     continue;
440    
441     tmp_vec.push_back(cms2.genps_p4().at(gidx));
442     }
443    
444     return tmp_vec;
445     }
446    
447     float efftools::getGenHT(float pt_cut, float eta_cut) {
448    
449     std::vector<LorentzVector> tmp_vec = getGenJets(pt_cut, eta_cut);
450    
451     float tmp_ht = 0.;
452     for(unsigned int idx = 0; idx < tmp_vec.size(); idx++) {
453     tmp_ht += tmp_vec.at(idx).pt();
454     }
455    
456     return tmp_ht;
457     }
458    
459     bool efftools::leptonOverlapsWithParton(LorentzVector p4, float parton_pt, float dr) {
460    
461     for (unsigned int idx = 0; idx < cms2.genps_p4().size(); idx++) {
462    
463     if (cms2.genps_status().at(idx) != 3)
464     continue;
465     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
466     continue;
467     if (cms2.genps_p4().at(idx).pt() < parton_pt)
468     continue;
469     if (ROOT::Math::VectorUtil::DeltaR(cms2.genps_p4().at(idx), p4) < dr)
470     return true;
471     }
472    
473     return false;
474     }
475    
476     enum DileptonHypType efftools::getHypType (int id1, int id2) {
477    
478     if (abs(id1) != 11 && abs(id1) != 13)
479     return DILEPTON_ALL;
480     if (abs(id2) != 11 && abs(id2) != 13)
481     return DILEPTON_ALL;
482    
483     if (abs(id1) == 11 && abs(id2) == 11)
484     return DILEPTON_EE;
485     else if (abs(id1) == 13 && abs(id2) == 13)
486     return DILEPTON_MUMU;
487     else
488     return DILEPTON_EMU;
489     }
490    
491    
492     int efftools::getGenParton (LorentzVector p4) {
493    
494     float tmp_dr = 99.;
495     int index = -1;
496     for (unsigned int jidx = 0; jidx < cms2.genps_p4().size(); jidx++) {
497    
498 fgolf 1.2 if (fabs(cms2.genps_p4().at(jidx).eta()) > 2.4)
499 fgolf 1.1 continue;
500     if (cms2.genps_p4().at(jidx).pt() < 40.)
501     continue;
502    
503     if (ROOT::Math::VectorUtil::DeltaR(cms2.genps_p4().at(jidx), p4) < tmp_dr) {
504     tmp_dr = ROOT::Math::VectorUtil::DeltaR(cms2.genps_p4().at(jidx), p4);
505     index = jidx;
506     }
507     }
508    
509     return index;
510     }