ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/ControlSelection.cc
Revision: 1.2
Committed: Mon Feb 13 09:57:14 2012 UTC (13 years, 3 months ago) by khahn
Content type: text/plain
Branch: MAIN
CVS Tags: synced_FSR_2, synced_FSR, synched2, synched
Changes since 1.1: +376 -506 lines
Log Message:
*** empty log message ***

File Contents

# Content
1 #include "PassHLT.h"
2 #include "MuonSelection.h"
3 #include "HZZBDTElectronSelection.h"
4 #include "SelectionStatus.h"
5 #include "SelectionDefs.h"
6 #include "ExternData.h"
7 #include "ControlSelection.h"
8
9 #include <fstream>
10 #include <sstream>
11 #include "TMath.h"
12 #include "TFile.h"
13 #include "CEffUser2D.h"
14
15 int getSFetaBin( float eta, vector<float> &SFeta ) {
16 for( int i=0; i<SFeta.size()-1; i++ ) {
17 if ( eta > SFeta[i] && eta <= SFeta[i+1] )
18 return i;
19 }
20 if( eta > SFeta[SFeta.size()-1] && eta <= 2.5 )
21 return SFeta.size()-1;
22 return -1;
23 }
24
25 void getSF( char * sfname, vector<float> &SFeta, vector<float> &SFscale, vector<float> &SFres ) {
26
27 ifstream sf;
28 sf.open(sfname);
29
30 std::string line;
31 float eta, scale, scaleerr, res, reserr;
32 while( getline(sf,line) ) {
33 stringstream ss(line);
34 ss >> eta;
35 ss >> scale;
36 ss >> scaleerr;
37 ss >> res;
38 ss >> reserr;
39 SFeta.push_back(eta);
40 SFscale.push_back(scale);
41 SFres.push_back(res);
42 }
43
44 };
45
46 TFile *mufrfile,*elefrfile;
47 TH2D *hFRmu,*hFRele;
48 TH2D *hFRmuErrl,*hFReleErrl;
49 TH2D *hFRmuErrh,*hFReleErrh;
50 CEffUser2D mufr, elefr;
51
52
53 void initEMUFR(ControlFlags ctrl)
54 {
55 TString mufrfilename("data/FakeRate/muon-ksWWdenomNoIso-fuckedfinal-fr.root");
56 cout << "loading fake file: " << mufrfilename << endl;
57 mufrfile = TFile::Open(mufrfilename);
58 hFRmu = (TH2D*)mufrfile->Get("frEtaPt");
59 hFRmuErrl = (TH2D*)mufrfile->Get("errlEtaPt");
60 hFRmuErrh = (TH2D*)mufrfile->Get("errhEtaPt");
61 mufr.loadEff(hFRmu,hFRmuErrl,hFRmuErrh);
62
63 // TString elefrfilename("data/fr-ele-bdt-medium.root");
64 TString elefrfilename("data/fr-ele-bdt-medium-looserdenom.root");
65 cout << "loading fake file: " << elefrfilename << endl;
66 elefrfile = TFile::Open(elefrfilename);
67 hFRele = (TH2D*)elefrfile->Get("frEtaPt");
68 hFReleErrl = (TH2D*)elefrfile->Get("errlEtaPt");
69 hFReleErrh = (TH2D*)elefrfile->Get("errhEtaPt");
70 elefr.loadEff(hFRele,hFReleErrl,hFReleErrh);
71
72 }
73
74
75 EventData fails_Control_selection(ControlFlags &ctrl, // input control
76 mithep::TEventInfo *info, // input event inof
77 TClonesArray *electronArr, // input electrons
78 SelectionStatus (*ElectronPreSelector)( ControlFlags &, const mithep::TElectron*),
79 SelectionStatus (*ElectronIDSelector)( ControlFlags &, const mithep::TElectron*),
80 SelectionStatus (*ElectronIsoSelector)( ControlFlags &, const mithep::TElectron*),
81 TClonesArray *muonArr,
82 SelectionStatus (*MuonPreSelector)( ControlFlags &, const mithep::TMuon*),
83 SelectionStatus (*MuonIDSelector)( ControlFlags &, const mithep::TMuon*),
84 SelectionStatus (*MuonIsoSelector)( ControlFlags &, const mithep::TMuon*) )
85 {
86
87 EventData ret;
88 unsigned evtfail = 0x0;
89
90 vector<float> SFeta, SFscale, SFres;
91 // if( ctrl.do_ele_scale_res ) {
92 // getSF( "results.txt", SFeta, SFscale, SFres );
93 // }
94
95 if( ctrl.debug ) {
96 cout << "Run: " << info->runNum
97 << "\tEvt: " << info->evtNum
98 << "\tLumi: " << info->lumiSec
99 << endl;
100 }
101
102
103 //********************************************************
104 // Trigger
105 //********************************************************
106 if( !ctrl.mc ) {
107 // if( !(passHLT(info->triggerBits, info->runNum, channel) ) ) {
108 if( !(passHLT(info->triggerBits, info->runNum, 999) ) ) {
109 evtfail |= (1<<EVTFAIL_TRIGGER);
110 ret.status.setStatus(0);
111 return ret;
112 }
113 // not accounting for overlap atm
114 RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);
115 if( !(rlrm.HasRunLumi(rl)) ) {
116 evtfail |= (1<<EVTFAIL_JSON);
117 ret.status.setStatus(0);
118 return ret;
119 }
120 }
121
122 if( ctrl.debug ) {
123 cout << "presel nlep: " << muonArr->GetEntries() + electronArr->GetEntries()
124 << "\tnmuon: " << muonArr->GetEntries()
125 << "\tnelectron: " << electronArr->GetEntries()
126 << endl;
127 }
128
129 //********************************************************
130 // Lepton Selection
131 //********************************************************
132 vector<SimpleLepton> lepvec;
133
134 //
135 if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
136 //----------------------------------------------------
137 for(Int_t i=0; i<muonArr->GetEntries(); i++)
138 {
139 mithep::TMuon *mu = (mithep::TMuon*)((*muonArr)[i]);
140
141 bool isFO = isMuFO(mu);
142
143 SelectionStatus musel;
144 musel |= (*MuonPreSelector)(ctrl,mu);
145 musel |= (*MuonIDSelector)(ctrl,mu);
146 musel |= (*MuonIsoSelector)(ctrl,mu);
147
148 if( ctrl.debug ) {
149 cout << "muon:: pt: " << mu->pt
150 << "\teta: " << mu->eta
151 << "\tmask: 0x" << hex << musel.getStatus() << dec
152 << endl;
153 }
154
155 if ( isFO ) {
156 SimpleLepton tmplep;
157 tmplep.vec->SetPtEtaPhiM(mu->pt,
158 mu->eta,
159 mu->phi,
160 MUON_MASS);
161 tmplep.type = 13;
162 tmplep.index = i;
163 tmplep.charge = mu->q;
164 tmplep.isoTrk = mu->trkIso03;
165 tmplep.isoEcal = mu->emIso03;
166 tmplep.isoHcal = mu->hadIso03;
167 tmplep.isoPF03 = mu->pfIso03;
168 tmplep.isoPF04 = mu->pfIso04;
169 tmplep.ip3dSig = mu->ip3dSig;
170 tmplep.is4l = false;
171 tmplep.isEB = (fabs(mu->eta) < 1.479 ? 1 : 0 );
172 tmplep.isTight = musel.tight();
173 tmplep.isLoose = musel.loose();
174
175 float fakerate = 1;
176 double ptmax = 35.; double etamax = 2.4; double epsilon = 0.001;
177 double ptval = (mu->pt<ptmax) ? mu->pt : ptmax-epsilon;
178 double etaval = (fabs(mu->eta)<etamax) ? fabs(mu->eta) : etamax-epsilon;
179 double rate = mufr.getEff(etaval,ptval);
180 fakerate = rate/(1-rate);
181 tmplep.FR = fakerate;
182 lepvec.push_back(tmplep);
183 if( ctrl.debug ) { cout << "muon passes ... " << endl;}
184 }
185 }
186
187 if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
188 //---------------------------------------------------------------------------------
189 for(Int_t i=0; i<electronArr->GetEntries(); i++)
190 {
191 const mithep::TElectron *ele = (mithep::TElectron*)((*electronArr)[i]);
192
193 Bool_t isMuonOverlap = kFALSE;
194 for (int k=0; k<lepvec.size(); ++k) {
195 TVector3 tmplep;
196 tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
197 if ( lepvec[k].isLoose && lepvec[k].type == 13 && lepvec[k].vec->Vect().DrEtaPhi(tmplep) < 0.1 ) {
198 if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
199 isMuonOverlap = kTRUE;
200 break;
201 }
202 }
203
204 SelectionStatus elesel;
205 elesel |= (*ElectronPreSelector)(ctrl,ele);
206 elesel |= (*ElectronIDSelector)(ctrl,ele);
207 elesel |= (*ElectronIsoSelector)(ctrl,ele);
208
209 // why doe we need LooseFO???
210 if ( isLooseEleFO(ele) && !isMuonOverlap ) {
211 SimpleLepton tmplep;
212 tmplep.vec->SetPtEtaPhiM( ele->pt,
213 ele->eta,
214 ele->phi,
215 ELECTRON_MASS );
216
217 tmplep.type = 11;
218 tmplep.index = i;
219 tmplep.charge = ele->q;
220 tmplep.isoTrk = ele->trkIso03;
221 tmplep.isoEcal = ele->emIso03;
222 tmplep.isoHcal = ele->hadIso03;
223 tmplep.isoPF03 = ele->pfIso03;
224 tmplep.isoPF04 = ele->pfIso04;
225 tmplep.ip3dSig = ele->ip3dSig;
226 tmplep.is4l = false;
227 tmplep.isEB = ele->isEB;
228 tmplep.isTight = elesel.tight();
229 tmplep.isLoose = elesel.loose();
230
231 float fakerate = 1;
232 // if ( !isEleTight ) {
233 double ptmax = 35.; double etamax = 2.5; double epsilon = 0.001;
234 double ptval = (ele->pt<ptmax) ? ele->pt : ptmax-epsilon;
235 double etaval = (fabs(ele->eta)<etamax) ? fabs(ele->eta) : etamax-epsilon;
236 double rate = elefr.getEff(etaval,ptval);
237 fakerate = rate/(1-rate);
238 // }
239 tmplep.FR = fakerate;
240
241 lepvec.push_back(tmplep);
242 if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
243 }
244 }
245
246 //********************************************************
247 // Dump Stuff
248 //********************************************************
249 sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
250 int nmu=0, nele=0;
251 for( int i=0; i<lepvec.size(); i++ ) {
252 if(ctrl.debug) cout << "lepvec :: index: " << i
253 << "\tpt: " << lepvec[i].vec->Pt()
254 << "\ttype: " << lepvec[i].type
255 << endl;
256 if( abs(lepvec[i].type) == 11 ) nele++;
257 else nmu++;
258 }
259 if( ctrl.debug ) {
260 cout << "postsel nlep: " << lepvec.size()
261 << "\tnmuon: " << nmu
262 << "\tnelectron: " << nele
263 << endl;
264 }
265
266 //******************************************************************************
267 //Z1 Selection
268 //******************************************************************************
269 int Z1LeptonPlusIndex = -1;
270 int Z1LeptonMinusIndex = -1;
271 double BestZ1Mass = -999;
272 if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
273 for(int i = 0; i < lepvec.size(); i++) {
274 if ( !(lepvec[i].isLoose) ) continue;
275 for(int j = i+1; j < lepvec.size(); j++) {
276 if ( !(lepvec[j].isLoose) ) continue;
277 if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
278 if (!(lepvec[i].vec->Pt() > 20.0 || lepvec[j].vec->Pt() > 20.0)) continue;
279 if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
280 if (!(lepvec[i].vec->Pt() > 10.0 && lepvec[j].vec->Pt() > 10.0)) continue;
281 if( ctrl.debug ) { cout << "\tthe other is > 10 GeV" << endl; }
282 if (lepvec[i].charge == lepvec[j].charge) continue;
283 if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
284 if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;
285 if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
286
287
288 //Make Z1 hypothesis
289 TLorentzVector leptonPlus, leptonMinus;
290 if ( lepvec[i].charge > 0 ) {
291 leptonPlus = *(lepvec[i].vec);
292 leptonMinus = *(lepvec[j].vec);
293 } else {
294 leptonPlus = *(lepvec[j].vec);
295 leptonMinus = *(lepvec[i].vec);
296 }
297
298 float tmpZ1Mass = (leptonPlus+leptonMinus).M();
299 if( ctrl.debug ) cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
300 if( tmpZ1Mass > 60 ) {
301 if (fabs(tmpZ1Mass - 91.1876) < fabs(BestZ1Mass - 91.1876)) {
302 BestZ1Mass = tmpZ1Mass;
303 if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass
304 << "\tdM: " << fabs(BestZ1Mass - 91.1876)
305 << endl;
306 if (lepvec[i].charge > 0) {
307 Z1LeptonPlusIndex = i;
308 Z1LeptonMinusIndex = j;
309 } else {
310 Z1LeptonPlusIndex = j;
311 Z1LeptonMinusIndex = i;
312 }
313 }
314 }
315 }
316 }
317 // stop if no Z1 candidate is found
318 if( BestZ1Mass < 0 ) {
319 evtfail |= (1<<EVTFAIL_Z1);
320 ret.status.setStatus(0);
321 return ret;
322 }
323 if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum << endl;
324 if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
325 TLorentzVector Z1LeptonPlus = *(lepvec[Z1LeptonPlusIndex].vec);
326 TLorentzVector Z1LeptonMinus = *(lepvec[Z1LeptonMinusIndex].vec);
327 TLorentzVector Z1Candidate = Z1LeptonPlus + Z1LeptonMinus;
328
329
330 //******************************************************************************
331 // Z1 + l
332 //******************************************************************************
333 if( lepvec.size() < 3 ) {
334 ret.status.setStatus(0);
335 return ret;
336 }
337
338 //******************************************************************************
339 // 4l/Z2 Selection
340 //******************************************************************************
341 Int_t Z2LeptonPlusIndex = -1;
342 Int_t Z2LeptonMinusIndex = -1;
343 Double_t BestZ2Mass = -1;
344 if( ctrl.debug ) cout << "hey ... looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
345 for(int i = 0; i < lepvec.size(); ++i) {
346 for(int j = i+1; j < lepvec.size(); ++j) {
347 if( ctrl.debug ) cout << "i: " << i << "\tj: " << j << endl;
348 if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) {
349 if( ctrl.debug ) cout << "\ti matches a Z1 index, skipping ..." << endl;
350 continue; //skip Z1 leptons
351 }
352 if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) {
353 if( ctrl.debug ) cout << "\tj matches a Z1 index, skipping ..." << endl;
354 continue; //skip Z1 leptons
355 }
356
357 //Make Z2 hypothesis
358 TLorentzVector leptonPlus, leptonMinus;
359 if (lepvec[i].charge > 0 ) {
360 leptonPlus = *(lepvec[i].vec);
361 leptonMinus = *(lepvec[j].vec);
362 } else {
363 leptonPlus = *(lepvec[j].vec);
364 leptonMinus = *(lepvec[i].vec);
365 }
366 TLorentzVector dilepton = leptonPlus+leptonMinus;
367 TLorentzVector fourLepton = Z1Candidate + dilepton;
368
369
370 //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
371 if (Z2LeptonPlusIndex < 0) {
372 if (lepvec[i].charge > 0) {
373 Z2LeptonPlusIndex = i;
374 Z2LeptonMinusIndex = j;
375 } else {
376 Z2LeptonPlusIndex = j;
377 Z2LeptonMinusIndex = i;
378 }
379 } else {
380 Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec->Pt();
381 Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec->Pt();
382 if (lepvec[Z2LeptonMinusIndex].vec->Pt() > BestPairPtMax) {
383 BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec->Pt();
384 BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec->Pt();
385 }
386
387 Double_t CurrentPairPtMax = lepvec[i].vec->Pt();
388 Double_t CurrentPairPtMin = lepvec[j].vec->Pt();
389 if (lepvec[j].vec->Pt() > CurrentPairPtMax) {
390 CurrentPairPtMax = lepvec[j].vec->Pt();
391 CurrentPairPtMin = lepvec[i].vec->Pt();
392 }
393
394 if (CurrentPairPtMax > BestPairPtMax) {
395 if (lepvec[i].charge > 0) {
396 Z2LeptonPlusIndex = i;
397 Z2LeptonMinusIndex = j;
398 } else {
399 Z2LeptonPlusIndex = j;
400 Z2LeptonMinusIndex = i;
401 }
402 } else if (CurrentPairPtMax == BestPairPtMax) {
403 if (CurrentPairPtMin > BestPairPtMin) {
404 if (lepvec[i].charge > 0) {
405 Z2LeptonPlusIndex = i;
406 Z2LeptonMinusIndex = j;
407 } else {
408 Z2LeptonPlusIndex = j;
409 Z2LeptonMinusIndex = i;
410 }
411 }
412 }
413 }
414 }
415 }
416
417 // stop if no Z2 candidate is found
418 if (Z2LeptonPlusIndex == -1) {
419 evtfail |= ( 1<<EVTFAIL_4L );
420 ret.status.setStatus(0);
421 return ret;
422 }
423 if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
424 if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
425 << "\tminusindex: " << Z2LeptonMinusIndex << endl;
426 TLorentzVector Z2LeptonPlus = *(lepvec[Z2LeptonPlusIndex].vec);
427 TLorentzVector Z2LeptonMinus = *(lepvec[Z2LeptonMinusIndex].vec);
428 TLorentzVector Z2Candidate = Z2LeptonPlus+Z2LeptonMinus;
429 TLorentzVector ZZSystem = Z1Candidate + Z2Candidate;
430
431 lepvec[Z1LeptonPlusIndex].is4l = true;
432 lepvec[Z1LeptonMinusIndex].is4l = true;
433 lepvec[Z2LeptonPlusIndex].is4l = true;
434 lepvec[Z2LeptonMinusIndex].is4l = true;
435
436 //
437 // need a FO def for high IP above ...
438 //
439 //***************************************************************
440 // IP significance
441 //***************************************************************
442 /*
443 bool failip = false;
444 for( int i=0; i<lepvec.size(); i++ ) {
445 if( !(lepvec[i].is4l) ) continue;
446 if( lepvec[i].ip3dSig > 4 ) {
447 if( ( i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex ))
448 failip=true;
449 }
450 }
451 if( failip ) {
452 //evtfail |= (1<<EVTFAIL_IP );
453 //return evtfail;
454 //h_evtfail->Fill( evtfail, eventweight );
455 // h_evtfail->Fill( evtfail );
456 // cout << "evtfail: " << hex << evtfail << dec << endl;
457 // continue;
458 }
459 */
460
461 //***************************************************************
462 // remaining kinematic cuts
463 //***************************************************************
464 double Z2massCut=0;
465 if ( ctrl.kinematics == "loose" ) Z2massCut = 0;
466 else if ( ctrl.kinematics == "tight" ) Z2massCut = 0;
467 else { cout << "error! kinematic tightness not defined!" << endl; assert(0); }
468
469 if ( Z1Candidate.M() > 120 ||
470 Z2Candidate.M() < Z2massCut ||
471 Z2Candidate.M() > 120 ||
472 !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec->Pt() > 20.0) ||
473 !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec->Pt() > 10.0)
474 ) {
475 evtfail |= (1<<EVTFAIL_KINEMATICS );
476 ret.status.setStatus(0);
477 return ret;
478 }
479
480 int channel;
481 if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
482 if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
483 if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) ||
484 (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
485
486
487
488
489 //***************************************************************
490 // finish
491 //***************************************************************
492 if( !evtfail ) {
493 ret.status.setStatus(SelectionStatus::EVTPASS);
494 ret.Z1leptons.push_back(lepvec[Z1LeptonMinusIndex]);
495 ret.Z1leptons.push_back(lepvec[Z1LeptonPlusIndex]);
496 ret.Z2leptons.push_back(lepvec[Z2LeptonMinusIndex]);
497 ret.Z2leptons.push_back(lepvec[Z2LeptonPlusIndex]);
498 }
499 return ret;
500
501 }
502
503
504
505
506