ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/Selection.cc
Revision: 1.19
Committed: Mon Feb 13 09:47:12 2012 UTC (13 years, 3 months ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.18: +453 -827 lines
Log Message:
*** empty log message ***

File Contents

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