ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/SFrameAnalysis/src/Cleaner.cxx
Revision: 1.23
Committed: Mon Dec 3 20:08:29 2012 UTC (12 years, 5 months ago) by bazterra
Content type: text/plain
Branch: MAIN
CVS Tags: Dec-05-2012-v1, Dec-04-2012-v1, Dec-03-2012-v4
Changes since 1.22: +467 -453 lines
Log Message:
Adding inverse qcd selection for electron channel.

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    
221     void Cleaner::ElectronCleaner_noID_noIso(double ptmin, double etamax)
222     {
223    
224     std::vector<Electron> good_eles;
225     for(unsigned int i=0; i<bcc->electrons->size(); ++i) {
226     Electron ele = bcc->electrons->at(i);
227     if(ele.pt()>ptmin) {
228     if(fabs(ele.eta())<etamax) {
229     good_eles.push_back(ele);
230     }
231     }
232     }
233    
234     bcc->electrons->clear();
235    
236     for(unsigned int i=0; i<good_eles.size(); ++i) {
237     bcc->electrons->push_back(good_eles[i]);
238     }
239     sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
240     resetEventCalc();
241     }
242    
243     void Cleaner::ElectronCleaner_noIso(double ptmin, double etamax)
244     {
245     ElectronCleaner_noID_noIso(ptmin, etamax);
246     std::vector<Electron> good_eles;
247     for(unsigned int i=0; i<bcc->electrons->size(); ++i) {
248     Electron ele = bcc->electrons->at(i);
249     if(fabs(ele.supercluster_eta())<1.4442 || fabs(ele.supercluster_eta())>1.5660) {
250     if(bcc->pvs->size()>0) {
251     if(fabs(ele.gsfTrack_dxy_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y()))<0.02) {
252     if(fabs(ele.gsfTrack_dz_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y(), bcc->pvs->at(0).z()))<0.1) {
253     if(ele.passconversionveto()) {
254     if(ele.mvaTrigV0()>0.0) {
255     if(ele.eleID(Electron::e_Tight)) {
256     good_eles.push_back(ele);
257     }
258     }
259     }
260     }
261 pturner 1.22 }
262     }
263 bazterra 1.23 }
264     }
265    
266     bcc->electrons->clear();
267    
268     for(unsigned int i=0; i<good_eles.size(); ++i) {
269     bcc->electrons->push_back(good_eles[i]);
270     }
271     sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
272     resetEventCalc();
273     }
274    
275    
276     void Cleaner::ElectronCleaner_noIso_reverse(double ptmin, double etamax)
277     {
278    
279     ElectronCleaner_noID_noIso(ptmin, etamax);
280     std::vector<Electron> good_eles;
281     for(unsigned int i=0; i<bcc->electrons->size(); ++i) {
282     int pass = 0;
283     Electron ele = bcc->electrons->at(i);
284     if(fabs(ele.supercluster_eta())<1.4442 || fabs(ele.supercluster_eta())>1.5660) {
285     if(bcc->pvs->size()>0) {
286     if(fabs(ele.gsfTrack_dxy_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y()))<0.02) {
287     if(fabs(ele.gsfTrack_dz_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y(), bcc->pvs->at(0).z()))<0.1) {
288     if(ele.passconversionveto()) {
289     if(ele.mvaTrigV0()>0.0) {
290     if(ele.eleID(Electron::e_Tight)) {
291     // good_eles.push_back(ele);
292     pass = 1;
293     // cout << "Failed EleID req." << endl;
294     }
295     }
296     }
297     }
298     }
299     }
300     }
301    
302     if(pass == 0) {
303     good_eles.push_back(ele);
304     }
305    
306     }
307    
308     bcc->electrons->clear();
309    
310     for(unsigned int i=0; i<good_eles.size(); ++i) {
311     bcc->electrons->push_back(good_eles[i]);
312     }
313     sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
314     resetEventCalc();
315     }
316    
317     void Cleaner::ElectronCleaner(double ptmin, double etamax, double relisomax)
318     {
319    
320     ElectronCleaner_noIso(ptmin, etamax);
321     std::vector<Electron> good_eles;
322     for(unsigned int i=0; i<bcc->electrons->size(); ++i) {
323     Electron ele = bcc->electrons->at(i);
324     if(ele.relIsorho(bcc->rho)<relisomax) {
325     good_eles.push_back(ele);
326     }
327     }
328     bcc->electrons->clear();
329    
330     for(unsigned int i=0; i<good_eles.size(); ++i) {
331     bcc->electrons->push_back(good_eles[i]);
332     }
333     sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
334     resetEventCalc();
335     }
336    
337     void Cleaner::MuonCleaner_noID_noIso(double ptmin, double etamax)
338     {
339    
340     std::vector<Muon> good_mus;
341     for(unsigned int i=0; i<bcc->muons->size(); ++i) {
342     Muon mu = bcc->muons->at(i);
343     if(mu.pt()>ptmin) {
344     if(fabs(mu.eta())<etamax) {
345     good_mus.push_back(mu);
346     }
347     }
348     }
349     bcc->muons->clear();
350    
351     for(unsigned int i=0; i<good_mus.size(); ++i) {
352     bcc->muons->push_back(good_mus[i]);
353     }
354     sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
355     resetEventCalc();
356     }
357    
358     void Cleaner::MuonCleaner_noIso(double ptmin, double etamax)
359     {
360    
361     MuonCleaner_noID_noIso(ptmin, etamax);
362     std::vector<Muon> good_mus;
363     for(unsigned int i=0; i<bcc->muons->size(); ++i) {
364     Muon mu = bcc->muons->at(i);
365     if(mu.isGlobalMuon()) {
366     if(mu.isPFMuon()) {
367     if(mu.globalTrack_chi2()/mu.globalTrack_ndof()<10) {
368     if(mu.globalTrack_numberOfValidMuonHits()>0) {
369     if(mu.innerTrack_trackerLayersWithMeasurement()>5) {
370     if(mu.dB()<0.02) { //0.2 ?
371     if(fabs(mu.vertex_z()-bcc->pvs->at(0).z())<0.5) {
372     if(mu.innerTrack_numberOfValidPixelHits()>0) {
373     if(mu.numberOfMatchedStations()>1) {
374     good_mus.push_back(mu);
375     }
376     }
377     }
378     }
379     }
380     }
381     }
382     }
383     }
384     }
385     bcc->muons->clear();
386    
387     for(unsigned int i=0; i<good_mus.size(); ++i) {
388     bcc->muons->push_back(good_mus[i]);
389     }
390     sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
391     resetEventCalc();
392     }
393    
394     void Cleaner::MuonCleaner(double ptmin, double etamax, double relisomax)
395     {
396    
397     MuonCleaner_noIso(ptmin, etamax);
398     std::vector<Muon> good_mus;
399     for(unsigned int i=0; i<bcc->muons->size(); ++i) {
400     Muon mu = bcc->muons->at(i);
401     if(mu.relIso()<relisomax) {
402     good_mus.push_back(mu);
403     }
404     }
405     bcc->muons->clear();
406    
407     for(unsigned int i=0; i<good_mus.size(); ++i) {
408     bcc->muons->push_back(good_mus[i]);
409     }
410     sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
411     resetEventCalc();
412     }
413    
414     void Cleaner::TauCleaner(double ptmin, double etamax)
415     {
416    
417     std::vector<Tau> good_taus;
418     for(unsigned int i=0; i<bcc->taus->size(); ++i) {
419     Tau tau = bcc->taus->at(i);
420     if(tau.pt()>ptmin) {
421     if(fabs(tau.eta())<etamax) {
422     if(bcc->taus->at(i).decayModeFinding()) {
423     if(bcc->taus->at(i).byMediumCombinedIsolationDeltaBetaCorr()) {
424     if(bcc->taus->at(i).againstElectronTight()) {
425     if(bcc->taus->at(i).againstMuonTight()) {
426     good_taus.push_back(tau);
427     }
428     }
429     }
430     }
431     }
432     }
433     }
434    
435     bcc->taus->clear();
436    
437     for(unsigned int i=0; i<good_taus.size(); ++i) {
438     bcc->taus->push_back(good_taus[i]);
439     }
440     sort(bcc->taus->begin(), bcc->taus->end(), HigherPt());
441     resetEventCalc();
442     }
443    
444    
445     void Cleaner::JetCleaner(double ptmin, double etamax, bool doPFID)
446     {
447    
448     std::vector<Jet> good_jets;
449     for(unsigned int i=0; i<bcc->jets->size(); ++i) {
450     Jet jet = bcc->jets->at(i);
451     if(jet.pt()>ptmin) {
452     if(fabs(jet.eta())<etamax) {
453     if(!doPFID || jet.pfID()) {
454     good_jets.push_back(jet);
455     }
456     }
457     }
458     }
459    
460     bcc->jets->clear();
461    
462     for(unsigned int i=0; i<good_jets.size(); ++i) {
463     bcc->jets->push_back(good_jets[i]);
464     }
465     sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
466     resetEventCalc();
467     }
468    
469     void Cleaner::TopJetCleaner(double ptmin, double etamax, bool doPFID)
470     {
471    
472     std::vector<TopJet> good_topjets;
473     for(unsigned int i=0; i<bcc->topjets->size(); ++i) {
474     TopJet topjet = bcc->topjets->at(i);
475     if(topjet.pt()>ptmin) {
476     if(fabs(topjet.eta())<etamax) {
477     if(!doPFID || topjet.pfID()) {
478     good_topjets.push_back(topjet);
479     }
480     }
481     }
482     }
483     bcc->topjets->clear();
484    
485     for(unsigned int i=0; i<good_topjets.size(); ++i) {
486     bcc->topjets->push_back(good_topjets[i]);
487     }
488     sort(bcc->topjets->begin(), bcc->topjets->end(), HigherPt());
489     resetEventCalc();
490     }
491    
492    
493     void Cleaner::PrimaryVertexCleaner(int ndofmax, double zmax, double rhomax)
494     {
495    
496     std::vector<PrimaryVertex> good_pvs;
497     for(unsigned int i=0; i<bcc->pvs->size(); ++i) {
498     PrimaryVertex pv = bcc->pvs->at(i);
499     if(pv.ndof()>=ndofmax && fabs(pv.z())<zmax && pv.rho()<=rhomax) {
500     good_pvs.push_back(pv);
501     }
502     }
503     bcc->pvs->clear();
504    
505     for(unsigned int i=0; i<good_pvs.size(); ++i) {
506     bcc->pvs->push_back(good_pvs[i]);
507     }
508     resetEventCalc();
509 peiffer 1.7 }