ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/SelectionZ4L.cc
Revision: 1.3
Committed: Fri Feb 17 15:19:16 2012 UTC (13 years, 3 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.2: +3 -5 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 "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_Z4L_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( !(lepvec[i].isTight) ) {
338 if( ctrl.debug) cout << "it's not tight, skipping ... " << endl;
339 continue;
340 }
341
342 for(int j = i+1; j < lepvec.size(); ++j) {
343 if( ctrl.debug) cout << "\t\tj: " << j
344 << "\tpt: " << lepvec[j].vec->Pt()
345 << "\ttype: " << lepvec[j].type
346 << endl;
347
348 if( !(lepvec[j].isTight) ) {
349 if( ctrl.debug) cout << "it's not tight, skipping ... " << endl;
350 continue;
351 }
352
353
354 if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) {
355 if( ctrl.debug) cout << "\ti matches a Z1 index, skipping ..." << endl;
356 continue; //skip Z1 leptons
357 }
358 if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) {
359 if( ctrl.debug) cout << "\tj matches a Z1 index, skipping ..." << endl;
360 continue; //skip Z1 leptons
361 }
362 if (lepvec[i].charge == lepvec[j].charge) {
363 if( ctrl.debug) cout << "\ti and j are same sign, skipping ..." << endl;
364 continue; //require opp sign
365 }
366 if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
367 if( ctrl.debug) cout << "\ti and j are not same flavor, skipping ..." << endl;
368 continue; //require same flavor
369 }
370
371
372 //Make Z2 hypothesis
373 TLorentzVector *leptonPlus, *leptonMinus;
374
375 if (lepvec[i].charge > 0 ) {
376 leptonPlus = lepvec[i].vec;
377 leptonMinus = lepvec[j].vec;
378 } else {
379 leptonPlus = lepvec[j].vec;
380 leptonMinus = lepvec[i].vec;
381 }
382
383 TLorentzVector dilepton = *leptonPlus + *leptonMinus;
384 TLorentzVector fourLepton = Z1Candidate + dilepton;
385
386 if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
387 if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;
388
389 if (!(dilepton.M() > 12.0)) continue;
390 if (!(fourLepton.M() > 100.0)) continue;
391
392 //for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
393 if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
394 TLorentzVector pair1 = Z1LeptonPlus + *leptonMinus;
395 TLorentzVector pair2 = Z1LeptonMinus + *leptonPlus;
396 if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
397 if (!(pair1.M() > 12 || pair2.M() > 12)) continue;
398 }
399
400
401 //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
402 if (Z2LeptonPlusIndex < 0) {
403 if (lepvec[i].charge > 0) {
404 Z2LeptonPlusIndex = i;
405 Z2LeptonMinusIndex = j;
406 } else {
407 Z2LeptonPlusIndex = j;
408 Z2LeptonMinusIndex = i;
409 }
410 } else {
411 Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec->Pt();
412 Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec->Pt();
413 if (lepvec[Z2LeptonMinusIndex].vec->Pt() > BestPairPtMax) {
414 BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec->Pt();
415 BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec->Pt();
416 }
417
418 Double_t CurrentPairPtMax = lepvec[i].vec->Pt();
419 Double_t CurrentPairPtMin = lepvec[j].vec->Pt();
420 if (lepvec[j].vec->Pt() > CurrentPairPtMax) {
421 CurrentPairPtMax = lepvec[j].vec->Pt();
422 CurrentPairPtMin = lepvec[i].vec->Pt();
423 }
424
425 if (CurrentPairPtMax > BestPairPtMax) {
426 if (lepvec[i].charge > 0) {
427 Z2LeptonPlusIndex = i;
428 Z2LeptonMinusIndex = j;
429 } else {
430 Z2LeptonPlusIndex = j;
431 Z2LeptonMinusIndex = i;
432 }
433 } else if (CurrentPairPtMax == BestPairPtMax) {
434 if (CurrentPairPtMin > BestPairPtMin) {
435 if (lepvec[i].charge > 0) {
436 Z2LeptonPlusIndex = i;
437 Z2LeptonMinusIndex = j;
438 } else {
439 Z2LeptonPlusIndex = j;
440 Z2LeptonMinusIndex = i;
441 }
442 }
443 }
444 }
445 }
446 }
447
448 // stop if no Z2 candidate is found
449 if (Z2LeptonPlusIndex == -1) {
450 evtfail |= ( 1<<EVTFAIL_4L );
451 // ret.status = evtfail;
452 ret.status.setStatus(0);
453 return ret;
454 }
455 if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
456 if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
457 << "\tminusindex: " << Z2LeptonMinusIndex << endl;
458 TLorentzVector Z2LeptonPlus = *(lepvec[Z2LeptonPlusIndex].vec);
459 TLorentzVector Z2LeptonMinus = *(lepvec[Z2LeptonMinusIndex].vec);
460 TLorentzVector Z2Candidate = Z2LeptonPlus+Z2LeptonMinus;
461 TLorentzVector ZZSystem = Z1Candidate + Z2Candidate;
462 lepvec[Z1LeptonPlusIndex].is4l = true;
463 lepvec[Z1LeptonMinusIndex].is4l = true;
464 lepvec[Z2LeptonPlusIndex].is4l = true;
465 lepvec[Z2LeptonMinusIndex].is4l = true;
466
467
468 //***************************************************************
469 // remaining cuts : kinematic & iso
470 //***************************************************************
471
472 if( !(PairwiseIsoSelector(ctrl, lepvec, info->rho)) ) {
473 ret.status.setStatus(0);
474 return ret;
475 }
476
477 double Z2massCut=0;
478 if ( ctrl.kinematics == "loose" ) Z2massCut = 12;
479 else if ( ctrl.kinematics == "tight" ) Z2massCut = 20;
480 else { cout << "error! kinematic tightness not defined!" << endl; assert(0); }
481
482 if ( Z1Candidate.M() > 120 ||
483 Z2Candidate.M() < Z2massCut ||
484 Z2Candidate.M() > 120 ||
485 !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec->Pt() > 20.0) ||
486 !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec->Pt() > 10.0)
487 ) {
488 evtfail |= (1<<EVTFAIL_KINEMATICS );
489 // ret.status = evtfail;
490 ret.status.setStatus(0);
491 return ret;
492 }
493
494 unsigned channel;
495 if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
496 if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
497 if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) ||
498 (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
499
500
501
502
503 if( ctrl.debug ) cout << "run: " << info->runNum
504 << "\tevt: " << info->evtNum
505 << "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
506 << "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
507 << "\tmZ1: " << Z1Candidate.M()
508 << "\tmZ2: " << Z2Candidate.M()
509 << "\tm4l: " << ZZSystem.M()
510 << "\tevtfail: " << hex << evtfail << dec
511 << "\ttrigbits: " << hex << info->triggerBits << dec
512 // << "\ttree: " << inputFiles[q][f]
513 << endl;
514
515
516
517 //***************************************************************
518 // finish
519 //***************************************************************
520
521 if( !evtfail ) {
522 ret.status.setStatus(SelectionStatus::EVTPASS);
523 ret.Z1leptons.push_back(lepvec[Z1LeptonMinusIndex]);
524 ret.Z1leptons.push_back(lepvec[Z1LeptonPlusIndex]);
525 ret.Z2leptons.push_back(lepvec[Z2LeptonMinusIndex]);
526 ret.Z2leptons.push_back(lepvec[Z2LeptonPlusIndex]);
527 }
528
529 return ret;
530 }
531
532