ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/SFrameAnalysis/src/Cleaner.cxx
Revision: 1.27
Committed: Wed Jan 9 19:03:22 2013 UTC (12 years, 3 months ago) by bazterra
Content type: text/plain
Branch: MAIN
CVS Tags: Jan-17-2013-v1, Jan-16-2012-v1, Jan-10-2012-v1, Jan-09-2012-v1
Changes since 1.26: +42 -1 lines
Log Message:
Adding the new electron mvaid based on top reference seleciton for 53x.

File Contents

# User Rev Content
1 peiffer 1.1 #include "../include/Cleaner.h"
2    
3 bazterra 1.23 Cleaner::Cleaner( BaseCycleContainer* input)
4     {
5 peiffer 1.1
6 bazterra 1.23 bcc = input;
7 peiffer 1.1
8     }
9    
10 bazterra 1.23 Cleaner::Cleaner()
11     {
12 peiffer 1.9
13 bazterra 1.23 ObjectHandler* objs = ObjectHandler::Instance();
14     bcc = objs->GetBaseCycleContainer();
15 peiffer 1.9
16     }
17    
18 bazterra 1.23 void Cleaner::resetEventCalc()
19     {
20 peiffer 1.9
21 bazterra 1.23 //reset everything in EventCalc except the event weight
22     EventCalc* evc = EventCalc::Instance();
23     double weight = evc->GetWeight();
24     evc->Reset();
25     evc->ProduceWeight(weight);
26 peiffer 1.9
27     }
28    
29 bazterra 1.23 void Cleaner::JetEnergyResolutionShifter(E_SystShift syst_shift, bool sort)
30     {
31 peiffer 1.2
32 bazterra 1.23 LorentzVector met(0,0,0,0);
33     if(bcc->met) {
34     met = bcc->met->v4();
35     }
36    
37     for(unsigned int i=0; i<bcc->jets->size(); ++i) {
38 peiffer 1.2
39 bazterra 1.23 float genpt = bcc->jets->at(i).genjet_pt();
40     //ignore unmatched jets (which have zero vector) or jets with very low pt:
41     if(genpt < 15.0) {
42     // std::cout << "1.0 | " << bcc->jets->at(i).pt() << " | " << bcc->jets->at(i).eta() << " | " << genpt << std::endl;
43     continue;
44     }
45 peiffer 1.15
46 bazterra 1.23 LorentzVector jet_v4 = bcc->jets->at(i).v4();
47    
48     LorentzVector jet_v4_raw = jet_v4*bcc->jets->at(i).JEC_factor_raw();
49 peiffer 1.16
50 bazterra 1.23 float recopt = jet_v4.pt();
51     float factor = 0.0;
52     float abseta = fabs(jet_v4.eta());
53 peiffer 1.2
54 bazterra 1.23 //numbers taken from https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
55     if(syst_shift==e_Default) {
56     if(abseta < 0.5)
57     factor = 0.052;
58     else if(abseta >= 0.5 && abseta <1.1)
59     factor = 0.057;
60     else if(abseta >= 1.1 && abseta <1.7)
61     factor = 0.096;
62     else if(abseta >= 1.7 && abseta <2.3)
63     factor = 0.134;
64     else if(abseta >= 2.3)
65     factor = 0.288;
66     else
67     factor = 0.288;
68     } else if(syst_shift==e_Up) {
69     if(abseta < 0.5)
70     factor = 0.115;
71     else if(abseta >= 0.5 && abseta <1.1)
72     factor = 0.114;
73     else if(abseta >= 1.1 && abseta <1.7)
74     factor = 0.161;
75     else if(abseta >= 1.7 && abseta <2.3)
76     factor = 0.228;
77     else if(abseta >= 2.3)
78     factor = 0.488;
79     } else if(syst_shift==e_Down) {
80     if(abseta < 0.5)
81     factor = -0.01;
82     else if(abseta >= 0.5 && abseta <1.1)
83     factor = 0.00;
84     else if(abseta >= 1.1 && abseta <1.7)
85     factor = 0.032;
86     else if(abseta >= 1.7 && abseta <2.3)
87     factor = 0.042;
88     else if(abseta >= 2.3)
89     factor = 0.089;
90     }
91 peiffer 1.2
92 bazterra 1.23 float ptscale = std::max(0.0f, 1 + factor * (recopt - genpt) / recopt);
93 peiffer 1.16
94 bazterra 1.23 jet_v4*=ptscale;
95 peiffer 1.16
96 bazterra 1.23 bcc->jets->at(i).set_v4(jet_v4);
97 peiffer 1.2
98 bazterra 1.23 //propagate JER shifts to MET
99     //if(jet_v4.Pt()>25) {
100     met += jet_v4_raw;
101     jet_v4_raw*=ptscale;
102     met -= jet_v4_raw;
103     //}
104 peiffer 1.16
105 bazterra 1.23 //std::cout << ptscale << " | " << jet_v4.Pt() << " | "<< jet_v4.eta() << " | " << genpt << std::endl;
106 peiffer 1.16
107 bazterra 1.23 }
108 peiffer 1.6
109 bazterra 1.23 //store changed MET
110     if(bcc->met) {
111     bcc->met->set_pt(met.Pt());
112     bcc->met->set_phi(met.Phi());
113     }
114    
115     if(sort) std::sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
116     resetEventCalc();
117 peiffer 1.2 }
118    
119 peiffer 1.1
120 bazterra 1.23 void Cleaner::JetLeptonSubtractor(FactorizedJetCorrector *corrector, bool sort)
121     {
122    
123     //std::cout<< "event: " <<bcc->event <<std::endl;
124     //std::cout<< "ID| pt raw | PAT pt | area | rho | corr | old cor| JER ptscale | pt new | eta new | ptgen " <<std::endl;
125    
126 peiffer 1.7
127 bazterra 1.23 for(unsigned int i=0; i<bcc->jets->size(); ++i) {
128 peiffer 1.20
129 bazterra 1.23 LorentzVector jet_v4_raw = bcc->jets->at(i).v4()*bcc->jets->at(i).JEC_factor_raw();
130 peiffer 1.7
131 bazterra 1.23 //subtract lepton momenta from raw jet momentum
132 peiffer 1.7
133 bazterra 1.23 //std::cout << i+1 << " | " << jet_v4_raw.Pt() << std::endl;
134     double ele_energy = bcc->jets->at(i).chargedEmEnergyFraction()*jet_v4_raw.E();
135     double mu_energy = bcc->jets->at(i).muonEnergyFraction()*jet_v4_raw.E();
136 peiffer 1.11
137 bazterra 1.23 if(bcc->electrons) {
138     for(unsigned int j=0; j<bcc->electrons->size(); ++j) {
139     if(bcc->jets->at(i).deltaR(bcc->electrons->at(j))<0.5) {
140 peiffer 1.13 // if(jet_v4_raw.pt() >= bcc->electrons->at(j).pt()){
141 bazterra 1.23 jet_v4_raw -= bcc->electrons->at(j).v4();
142 peiffer 1.13 // bcc->jets->at(i).set_electronMultiplicity(bcc->jets->at(i).electronMultiplicity()-1);
143     // ele_energy -= bcc->electrons->at(j).energy();
144 bazterra 1.23 //}
145     //else{
146     //std::cout << bcc->jets->at(i).pt()<< " -ele pt:"<< bcc->electrons->at(j).pt() << std::endl;
147     //jet_v4_raw -= jet_v4_raw;
148     //}
149     }
150     }
151     }
152     if(bcc->muons) {
153     for(unsigned int j=0; j<bcc->muons->size(); ++j) {
154     if(bcc->jets->at(i).deltaR(bcc->muons->at(j))<0.5) {
155     //if(jet_v4_raw.pt()>= bcc->muons->at(j).pt()){
156     jet_v4_raw -= bcc->muons->at(j).v4();
157 peiffer 1.13 // bcc->jets->at(i).set_muonMultiplicity(bcc->jets->at(i).muonMultiplicity()-1);
158     // mu_energy -= bcc->muons->at(j).energy();
159 bazterra 1.23 //}
160     //else{
161     //std::cout << bcc->jets->at(i).pt() << " -mu pt:"<< bcc->muons->at(j).pt() << std::endl;
162     //jet_v4_raw -= jet_v4_raw;
163     //}
164     }
165     }
166     }
167     if(ele_energy<=jet_v4_raw.E())
168     bcc->jets->at(i).set_chargedEmEnergyFraction(ele_energy/jet_v4_raw.E());
169     if(mu_energy<=jet_v4_raw.E())
170     bcc->jets->at(i).set_muonEnergyFraction(mu_energy/jet_v4_raw.E());
171    
172     //apply jet energy corrections to modified raw momentum
173     corrector->setJetPt(jet_v4_raw.Pt());
174     corrector->setJetEta(jet_v4_raw.Eta());
175     corrector->setJetE(jet_v4_raw.E());
176     corrector->setJetA(bcc->jets->at(i).jetArea());
177     corrector->setRho(bcc->rho);
178     corrector->setNPV(bcc->pvs->size());
179    
180     float correctionfactor = corrector->getCorrection();
181    
182     LorentzVector jet_v4_corrected = jet_v4_raw *correctionfactor;
183    
184     //std::cout << i+1 << " | " << jet_v4_raw.Pt() << " | " << bcc->jets->at(i).pt() << " | " << bcc->jets->at(i).jetArea()<< " | " << bcc->rho << " | " << correctionfactor << " | " << 1./bcc->jets->at(i).JEC_factor_raw() << " | " << std::endl;
185    
186     bcc->jets->at(i).set_v4(jet_v4_corrected);
187     bcc->jets->at(i).set_JEC_factor_raw(1./correctionfactor);
188    
189     }
190    
191     if(sort) std::sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
192     resetEventCalc();
193     }
194    
195     void Cleaner::JetRecorrector(FactorizedJetCorrector *corrector, bool sort)
196     {
197    
198     for(unsigned int i=0; i<bcc->jets->size(); ++i) {
199    
200     LorentzVector jet_v4_raw = bcc->jets->at(i).v4()*bcc->jets->at(i).JEC_factor_raw();
201     corrector->setJetPt(jet_v4_raw.Pt());
202     corrector->setJetEta(jet_v4_raw.Eta());
203     corrector->setJetE(jet_v4_raw.E());
204     corrector->setJetA(bcc->jets->at(i).jetArea());
205     corrector->setRho(bcc->rho);
206     corrector->setNPV(bcc->pvs->size());
207    
208     float correctionfactor = corrector->getCorrection();
209    
210     LorentzVector jet_v4_corrected = jet_v4_raw *correctionfactor;
211    
212     bcc->jets->at(i).set_v4(jet_v4_corrected);
213     bcc->jets->at(i).set_JEC_factor_raw(1./correctionfactor);
214     }
215    
216     if(sort) std::sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
217     resetEventCalc();
218     }
219    
220 bazterra 1.27 // Old electron id
221     /*bool Cleaner::passElectronId(BaseCycleContainer * bcc, unsigned int index)
222 bazterra 1.26 {
223     Electron ele = bcc->electrons->at(index);
224     if(fabs(ele.supercluster_eta())<1.4442 || fabs(ele.supercluster_eta())>1.5660) {
225     if(bcc->pvs->size()>0) {
226     if(fabs(ele.gsfTrack_dxy_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y()))<0.02) {
227     if(fabs(ele.gsfTrack_dz_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y(), bcc->pvs->at(0).z()))<0.1) {
228     if(ele.passconversionveto()) {
229     if(ele.mvaTrigV0()>0.0) {
230     if(ele.eleID(Electron::e_Tight)) {
231     return true;
232     }
233     }
234     }
235     }
236     }
237     }
238     }
239     return false;
240 bazterra 1.27 }*/
241    
242     // Base on the new top simple cut recomendation for 53x
243     /*bool Cleaner::passElectronId(BaseCycleContainer * bcc, unsigned int index)
244     {
245     Electron ele = bcc->electrons->at(index);
246     if(fabs(ele.supercluster_eta())<1.4442 || fabs(ele.supercluster_eta())>1.5660) {
247     if(bcc->pvs->size()>0) {
248     if(fabs(ele.gsfTrack_dxy_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y()))<0.02) {
249     if(fabs(ele.gsfTrack_dz_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y(), bcc->pvs->at(0).z()))<0.1) {
250     if(ele.passconversionveto()) {
251     if(ele.eleID(Electron::e_Tight)) {
252     return true;
253     }
254     }
255     }
256     }
257     }
258     }
259     return false;
260     }*/
261    
262     // Base on the new top mva recomendation for 53x
263     bool Cleaner::passElectronId(BaseCycleContainer * bcc, unsigned int index)
264     {
265     Electron ele = bcc->electrons->at(index);
266     if(fabs(ele.supercluster_eta())<1.4442 || fabs(ele.supercluster_eta())>1.5660) {
267     if(bcc->pvs->size()>0) {
268     if(fabs(ele.gsfTrack_dxy_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y()))<0.02) {
269     if(ele.passconversionveto()) {
270     if(ele.mvaTrigV0()>0.9) {
271     if (ele.gsfTrack_trackerExpectedHitsInner_numberOfLostHits() <= 0) {
272     return true;
273     }
274     }
275     }
276     }
277     }
278     }
279     return false;
280 bazterra 1.26 }
281    
282     void Cleaner::ElectronCleaner_noID_noIso(double ptmin, double etamax)
283 bazterra 1.23 {
284     std::vector<Electron> good_eles;
285     for(unsigned int i=0; i<bcc->electrons->size(); ++i) {
286     Electron ele = bcc->electrons->at(i);
287     if(ele.pt()>ptmin) {
288     if(fabs(ele.eta())<etamax) {
289     good_eles.push_back(ele);
290     }
291     }
292     }
293    
294     bcc->electrons->clear();
295    
296     for(unsigned int i=0; i<good_eles.size(); ++i) {
297     bcc->electrons->push_back(good_eles[i]);
298     }
299     sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
300     resetEventCalc();
301     }
302    
303 bazterra 1.26 void Cleaner::ElectronCleaner_noIso(double ptmin, double etamax, bool reverseID)
304     {
305     ElectronCleaner_noID_noIso(ptmin, etamax);
306     std::vector<Electron> good_eles;
307     for(unsigned int i=0; i<bcc->electrons->size(); ++i) {
308     Electron ele = bcc->electrons->at(i);
309     if(!reverseID && passElectronId(bcc, i))
310     good_eles.push_back(ele);
311     else if (reverseID && !passElectronId(bcc, i))
312     good_eles.push_back(ele);
313     }
314     bcc->electrons->clear();
315    
316     for(unsigned int i=0; i<good_eles.size(); ++i) {
317     bcc->electrons->push_back(good_eles[i]);
318     }
319     sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
320     resetEventCalc();
321     }
322    
323    
324     void Cleaner::ElectronCleaner(double ptmin, double etamax, double relisomax, bool reverseID, bool reverseIso)
325 bazterra 1.23 {
326 bazterra 1.26 ElectronCleaner_noIso(ptmin, etamax, reverseID);
327 bazterra 1.23 std::vector<Electron> good_eles;
328     for(unsigned int i=0; i<bcc->electrons->size(); ++i) {
329     Electron ele = bcc->electrons->at(i);
330 bazterra 1.26 if(!reverseIso && ele.relIsorho(bcc->rho)<relisomax)
331     good_eles.push_back(ele);
332     else if (reverseIso && ele.relIsorho(bcc->rho)>=relisomax)
333 bazterra 1.23 good_eles.push_back(ele);
334     }
335     bcc->electrons->clear();
336    
337     for(unsigned int i=0; i<good_eles.size(); ++i) {
338     bcc->electrons->push_back(good_eles[i]);
339     }
340     sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
341     resetEventCalc();
342     }
343    
344 bazterra 1.26
345 bazterra 1.23 void Cleaner::MuonCleaner_noID_noIso(double ptmin, double etamax)
346     {
347     std::vector<Muon> good_mus;
348     for(unsigned int i=0; i<bcc->muons->size(); ++i) {
349     Muon mu = bcc->muons->at(i);
350     if(mu.pt()>ptmin) {
351     if(fabs(mu.eta())<etamax) {
352     good_mus.push_back(mu);
353     }
354     }
355     }
356     bcc->muons->clear();
357    
358     for(unsigned int i=0; i<good_mus.size(); ++i) {
359     bcc->muons->push_back(good_mus[i]);
360     }
361     sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
362     resetEventCalc();
363     }
364    
365     void Cleaner::MuonCleaner_noIso(double ptmin, double etamax)
366     {
367    
368     MuonCleaner_noID_noIso(ptmin, etamax);
369     std::vector<Muon> good_mus;
370     for(unsigned int i=0; i<bcc->muons->size(); ++i) {
371     Muon mu = bcc->muons->at(i);
372     if(mu.isGlobalMuon()) {
373     if(mu.isPFMuon()) {
374     if(mu.globalTrack_chi2()/mu.globalTrack_ndof()<10) {
375     if(mu.globalTrack_numberOfValidMuonHits()>0) {
376     if(mu.innerTrack_trackerLayersWithMeasurement()>5) {
377 peiffer 1.24 if(mu.dB()<0.2) {
378 bazterra 1.23 if(fabs(mu.vertex_z()-bcc->pvs->at(0).z())<0.5) {
379     if(mu.innerTrack_numberOfValidPixelHits()>0) {
380     if(mu.numberOfMatchedStations()>1) {
381     good_mus.push_back(mu);
382     }
383     }
384     }
385     }
386     }
387     }
388     }
389     }
390     }
391     }
392     bcc->muons->clear();
393    
394     for(unsigned int i=0; i<good_mus.size(); ++i) {
395     bcc->muons->push_back(good_mus[i]);
396     }
397     sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
398     resetEventCalc();
399     }
400    
401     void Cleaner::MuonCleaner(double ptmin, double etamax, double relisomax)
402     {
403    
404     MuonCleaner_noIso(ptmin, etamax);
405     std::vector<Muon> good_mus;
406     for(unsigned int i=0; i<bcc->muons->size(); ++i) {
407     Muon mu = bcc->muons->at(i);
408     if(mu.relIso()<relisomax) {
409     good_mus.push_back(mu);
410     }
411     }
412     bcc->muons->clear();
413    
414     for(unsigned int i=0; i<good_mus.size(); ++i) {
415     bcc->muons->push_back(good_mus[i]);
416     }
417     sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
418     resetEventCalc();
419     }
420    
421     void Cleaner::TauCleaner(double ptmin, double etamax)
422     {
423    
424     std::vector<Tau> good_taus;
425     for(unsigned int i=0; i<bcc->taus->size(); ++i) {
426     Tau tau = bcc->taus->at(i);
427     if(tau.pt()>ptmin) {
428     if(fabs(tau.eta())<etamax) {
429     if(bcc->taus->at(i).decayModeFinding()) {
430     if(bcc->taus->at(i).byMediumCombinedIsolationDeltaBetaCorr()) {
431     if(bcc->taus->at(i).againstElectronTight()) {
432     if(bcc->taus->at(i).againstMuonTight()) {
433     good_taus.push_back(tau);
434     }
435     }
436     }
437     }
438     }
439     }
440     }
441    
442     bcc->taus->clear();
443    
444     for(unsigned int i=0; i<good_taus.size(); ++i) {
445     bcc->taus->push_back(good_taus[i]);
446     }
447     sort(bcc->taus->begin(), bcc->taus->end(), HigherPt());
448     resetEventCalc();
449     }
450    
451    
452     void Cleaner::JetCleaner(double ptmin, double etamax, bool doPFID)
453     {
454    
455     std::vector<Jet> good_jets;
456     for(unsigned int i=0; i<bcc->jets->size(); ++i) {
457     Jet jet = bcc->jets->at(i);
458     if(jet.pt()>ptmin) {
459     if(fabs(jet.eta())<etamax) {
460     if(!doPFID || jet.pfID()) {
461     good_jets.push_back(jet);
462     }
463     }
464     }
465     }
466    
467     bcc->jets->clear();
468    
469     for(unsigned int i=0; i<good_jets.size(); ++i) {
470     bcc->jets->push_back(good_jets[i]);
471     }
472     sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
473     resetEventCalc();
474     }
475    
476     void Cleaner::TopJetCleaner(double ptmin, double etamax, bool doPFID)
477     {
478    
479     std::vector<TopJet> good_topjets;
480     for(unsigned int i=0; i<bcc->topjets->size(); ++i) {
481     TopJet topjet = bcc->topjets->at(i);
482     if(topjet.pt()>ptmin) {
483     if(fabs(topjet.eta())<etamax) {
484     if(!doPFID || topjet.pfID()) {
485     good_topjets.push_back(topjet);
486     }
487     }
488     }
489     }
490     bcc->topjets->clear();
491    
492     for(unsigned int i=0; i<good_topjets.size(); ++i) {
493     bcc->topjets->push_back(good_topjets[i]);
494     }
495     sort(bcc->topjets->begin(), bcc->topjets->end(), HigherPt());
496     resetEventCalc();
497     }
498    
499    
500     void Cleaner::PrimaryVertexCleaner(int ndofmax, double zmax, double rhomax)
501     {
502    
503     std::vector<PrimaryVertex> good_pvs;
504     for(unsigned int i=0; i<bcc->pvs->size(); ++i) {
505     PrimaryVertex pv = bcc->pvs->at(i);
506     if(pv.ndof()>=ndofmax && fabs(pv.z())<zmax && pv.rho()<=rhomax) {
507     good_pvs.push_back(pv);
508     }
509     }
510     bcc->pvs->clear();
511    
512     for(unsigned int i=0; i<good_pvs.size(); ++i) {
513     bcc->pvs->push_back(good_pvs[i]);
514     }
515     resetEventCalc();
516 peiffer 1.7 }