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

# 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     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     }