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