ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/SelectionZ4L.cc
Revision: 1.4
Committed: Wed Feb 29 09:43:13 2012 UTC (13 years, 2 months ago) by anlevin
Content type: text/plain
Branch: MAIN
CVS Tags: synced_FSR_2, synced_FSR, synched2, synched
Changes since 1.3: +21 -23 lines
Log Message:
changed Z4l selection code

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 SimpleLepton tmplep;
108 float pt = mu->pt;
109 tmplep.vecorig->SetPtEtaPhiM(pt,
110 mu->eta,
111 mu->phi,
112 MUON_MASS);
113
114 if( ctrl.do_escale_up ) {
115 pt=scale_smear_muon_Up(pt, 1, r);
116 }
117 if( ctrl.do_escale_down ) {
118 pt=scale_smear_muon_Down(pt, 1, r);
119 }
120
121 tmplep.vec->SetPtEtaPhiM(pt,
122 mu->eta,
123 mu->phi,
124 MUON_MASS);
125
126 tmplep.type = 13;
127 tmplep.index = i;
128 tmplep.charge = mu->q;
129 tmplep.isoTrk = mu->trkIso03;
130 tmplep.isoEcal = mu->emIso03;
131 tmplep.isoHcal = mu->hadIso03;
132 tmplep.isoPF03 = mu->pfIso03;
133 tmplep.isoPF04 = mu->pfIso04;
134 tmplep.ip3dSig = mu->ip3dSig;
135 tmplep.is4l = false;
136 tmplep.isEB = (fabs(mu->eta) < 1.479 ? 1 : 0 );
137 tmplep.isTight = musel.tightID();
138 tmplep.isLoose = musel.looseID();
139 lepvec.push_back(tmplep);
140 if( ctrl.debug ) { cout << "muon passes ... " << endl;}
141 }
142 // }
143 }
144
145
146
147 //
148 if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
149 // --------------------------------------------------------------------------------
150 for(Int_t i=0; i<electronArr->GetEntries(); i++)
151 {
152 const mithep::TElectron *ele = (mithep::TElectron*)((*electronArr)[i]);
153
154 Bool_t isMuonOverlap = kFALSE;
155 for (int k=0; k<lepvec.size(); ++k) {
156 TVector3 tmplep;
157 tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
158 if ( lepvec[k].isLoose && lepvec[k].type == 13 && lepvec[k].vec->Vect().DrEtaPhi(tmplep) < 0.05 ) {
159 if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
160 isMuonOverlap = kTRUE;
161 break;
162 }
163 }
164
165 SelectionStatus elesel;
166 if( ctrl.debug ) cout << "--> status before anything: " << hex << elesel.getStatus() << dec << endl;
167 elesel |= (*ElectronPreSelector)(ctrl,ele);
168 if( ctrl.debug ) cout << "--> status after presel: " << hex << elesel.getStatus() << dec << endl;
169 elesel |= (*ElectronIDSelector)(ctrl,ele,cicCuts);
170 if( ctrl.debug ) cout << "--> status after ID: " << hex << elesel.getStatus() << dec << endl;
171
172 if( ctrl.debug ){
173 cout << "\tscEt: " << ele->scEt
174 << "\tscEta: " << ele->scEta
175 << "\tncluster: " << ele->ncluster
176 << "\tisorel: " << ele->pfIso04/ele->pt
177 << "\tstatus: " << hex << elesel.getStatus() << dec
178 << endl;
179 }
180
181 if ( elesel.passIDAndPre() && !isMuonOverlap )
182 {
183 SimpleLepton tmplep;
184
185
186 float pt = ele->pt;
187
188 tmplep.vecorig->SetPtEtaPhiM( pt,
189 ele->eta,
190 ele->phi,
191 ELECTRON_MASS );
192
193 if( ctrl.do_escale ) {
194 pt=scale_smear_electron(pt, ele->isEB, r);
195 }
196 if( ctrl.do_escale_up ) {
197 pt=scale_smear_electron_Up(pt, ele->isEB, r);
198 }
199 if( ctrl.do_escale_down ) {
200 pt=scale_smear_electron_Down(pt, ele->isEB, r);
201 }
202
203
204 tmplep.vec->SetPtEtaPhiM( pt,
205 ele->eta,
206 ele->phi,
207 ELECTRON_MASS );
208
209 tmplep.type = 11;
210 tmplep.index = i;
211 tmplep.charge = ele->q;
212 tmplep.isoTrk = ele->trkIso03;
213 tmplep.isoEcal = ele->emIso03;
214 tmplep.isoHcal = ele->hadIso03;
215 tmplep.isoPF03 = ele->pfIso03;
216 tmplep.isoPF04 = ele->pfIso04;
217 tmplep.ip3dSig = ele->ip3dSig;
218 tmplep.is4l = false;
219 tmplep.isEB = ele->isEB;
220 tmplep.scID = ele->scID;
221 tmplep.isTight = elesel.tightID();
222 tmplep.isLoose = elesel.looseID();
223 lepvec.push_back(tmplep);
224 if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
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 // Z1 Selection
250 //******************************************************************************
251 int Z1LeptonPlusIndex = -1;
252 int Z1LeptonMinusIndex = -1;
253 double BestZ1Mass = -999;
254 if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
255 for(int i = 0; i < lepvec.size(); ++i) {
256 if( !(lepvec[i].isTight) ) continue;
257 for(int j = i+1; j < lepvec.size(); ++j) {
258 if( !(lepvec[j].isTight) ) continue;
259 if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
260 if (!(lepvec[i].vec->Pt() > 20.0 || lepvec[j].vec->Pt() > 20.0)) continue;
261 if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
262 if (!(lepvec[i].vec->Pt() > 10.0 && lepvec[j].vec->Pt() > 10.0)) continue;
263 if( ctrl.debug ) { cout << "\tthe other is > 10 GeV" << endl; }
264 if (lepvec[i].charge == lepvec[j].charge) continue;
265 if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
266 if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;
267 if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
268
269 //Make Z1 hypothesis
270 TLorentzVector *leptonPlus, *leptonMinus;
271 if ( lepvec[i].charge > 0 ) {
272 leptonPlus = lepvec[i].vec;
273 leptonMinus = lepvec[j].vec;
274 } else {
275 leptonPlus = lepvec[j].vec;
276 leptonMinus = lepvec[i].vec;
277 }
278
279 float tmpZ1Mass = (*leptonPlus + *leptonMinus).M();
280 if( ctrl.debug ) cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
281 if( tmpZ1Mass > 60 ) {
282 if (fabs(tmpZ1Mass - Z_MASS) < fabs(BestZ1Mass - Z_MASS)) {
283 BestZ1Mass = tmpZ1Mass;
284 if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass
285 << "\tdM: " << fabs(BestZ1Mass - Z_MASS)
286 << endl;
287 if (lepvec[i].charge > 0) {
288 Z1LeptonPlusIndex = i;
289 Z1LeptonMinusIndex = j;
290 } else {
291 Z1LeptonPlusIndex = j;
292 Z1LeptonMinusIndex = i;
293 }
294 }
295 }
296 }
297 }
298 // stop if no Z1 candidate is found
299 if( BestZ1Mass < 0 ) {
300 evtfail |= (1<<EVTFAIL_Z1);
301 //ret.status = evtfail;
302 ret.status.setStatus(0);
303 return ret;
304 }
305 if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum << endl;
306 if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
307 TLorentzVector Z1LeptonPlus = *(lepvec[Z1LeptonPlusIndex].vec);
308 TLorentzVector Z1LeptonMinus = *(lepvec[Z1LeptonMinusIndex].vec);
309 TLorentzVector Z1Candidate = Z1LeptonPlus + Z1LeptonMinus;
310
311 //******************************************************************************
312 // Z1 + l
313 //******************************************************************************
314 if( lepvec.size() < 3 ) {
315 evtfail |= (1<<EVTFAIL_Z1_PLUSL);
316 //ret.status = evtfail;
317 ret.status.setStatus(0);
318 return ret;
319 }
320
321 //******************************************************************************
322 // 4l/Z2 Selection
323 //******************************************************************************
324 Int_t Z2LeptonPlusIndex = -1;
325 Int_t Z2LeptonMinusIndex = -1;
326 Double_t BestZ2Mass = -1;
327 if( ctrl.debug ) cout << "looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
328 for(int i = 0; i < lepvec.size(); ++i) {
329
330 if( ctrl.debug) cout << "i: " << i
331 << "\tpt: " << lepvec[i].vec->Pt()
332 << "\ttype: " << lepvec[i].type
333 << endl;
334
335 if( !(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( !(lepvec[j].isTight) ) {
347 if( ctrl.debug) cout << "it's not tight, skipping ... " << endl;
348 continue;
349 }
350
351
352 if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) {
353 if( ctrl.debug) cout << "\ti matches a Z1 index, skipping ..." << endl;
354 continue; //skip Z1 leptons
355 }
356 if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) {
357 if( ctrl.debug) cout << "\tj matches a Z1 index, skipping ..." << endl;
358 continue; //skip Z1 leptons
359 }
360 if (lepvec[i].charge == lepvec[j].charge) {
361 if( ctrl.debug) cout << "\ti and j are same sign, skipping ..." << endl;
362 continue; //require opp sign
363 }
364 if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
365 if( ctrl.debug) cout << "\ti and j are not same flavor, skipping ..." << endl;
366 continue; //require same flavor
367 }
368
369 //Make Z2 hypothesis
370 TLorentzVector *leptonPlus, *leptonMinus;
371
372 if (lepvec[i].charge > 0 ) {
373 leptonPlus = lepvec[i].vec;
374 leptonMinus = lepvec[j].vec;
375 } else {
376 leptonPlus = lepvec[j].vec;
377 leptonMinus = lepvec[i].vec;
378 }
379
380 TLorentzVector dilepton = *leptonPlus + *leptonMinus;
381 TLorentzVector fourLepton = Z1Candidate + dilepton;
382
383 if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
384 if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;
385
386 if (!(dilepton.M() > 2.0)) continue;
387 if (!(fourLepton.M() > 60.0)) continue;
388
389 //for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
390 if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
391 TLorentzVector pair1 = Z1LeptonPlus + *leptonMinus;
392 TLorentzVector pair2 = Z1LeptonMinus + *leptonPlus;
393 if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
394 if (!(pair1.M() > 2 && pair2.M() > 2)) continue;
395 }
396
397
398
399 //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
400 if (Z2LeptonPlusIndex < 0) {
401 if (lepvec[i].charge > 0) {
402 Z2LeptonPlusIndex = i;
403 Z2LeptonMinusIndex = j;
404 } else {
405 Z2LeptonPlusIndex = j;
406 Z2LeptonMinusIndex = i;
407 }
408 } else {
409 Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec->Pt();
410 Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec->Pt();
411 if (lepvec[Z2LeptonMinusIndex].vec->Pt() > BestPairPtMax) {
412 BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec->Pt();
413 BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec->Pt();
414 }
415
416 Double_t CurrentPairPtMax = lepvec[i].vec->Pt();
417 Double_t CurrentPairPtMin = lepvec[j].vec->Pt();
418 if (lepvec[j].vec->Pt() > CurrentPairPtMax) {
419 CurrentPairPtMax = lepvec[j].vec->Pt();
420 CurrentPairPtMin = lepvec[i].vec->Pt();
421 }
422
423 if (CurrentPairPtMax > BestPairPtMax) {
424 if (lepvec[i].charge > 0) {
425 Z2LeptonPlusIndex = i;
426 Z2LeptonMinusIndex = j;
427 } else {
428 Z2LeptonPlusIndex = j;
429 Z2LeptonMinusIndex = i;
430 }
431 } else if (CurrentPairPtMax == BestPairPtMax) {
432 if (CurrentPairPtMin > BestPairPtMin) {
433 if (lepvec[i].charge > 0) {
434 Z2LeptonPlusIndex = i;
435 Z2LeptonMinusIndex = j;
436 } else {
437 Z2LeptonPlusIndex = j;
438 Z2LeptonMinusIndex = i;
439 }
440 }
441 }
442 }
443 }
444 }
445
446 // stop if no Z2 candidate is found
447 if (Z2LeptonPlusIndex == -1) {
448 evtfail |= ( 1<<EVTFAIL_4L );
449 // ret.status = evtfail;
450 ret.status.setStatus(0);
451 return ret;
452 }
453 if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
454 if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
455 << "\tminusindex: " << Z2LeptonMinusIndex << endl;
456 TLorentzVector Z2LeptonPlus = *(lepvec[Z2LeptonPlusIndex].vec);
457 TLorentzVector Z2LeptonMinus = *(lepvec[Z2LeptonMinusIndex].vec);
458 TLorentzVector Z2Candidate = Z2LeptonPlus+Z2LeptonMinus;
459 TLorentzVector ZZSystem = Z1Candidate + Z2Candidate;
460 lepvec[Z1LeptonPlusIndex].is4l = true;
461 lepvec[Z1LeptonMinusIndex].is4l = true;
462 lepvec[Z2LeptonPlusIndex].is4l = true;
463 lepvec[Z2LeptonMinusIndex].is4l = true;
464
465
466
467 //***************************************************************
468 // remaining cuts : kinematic & iso
469 //***************************************************************
470
471 if( !(PairwiseIsoSelector(ctrl, lepvec, info->rho)) ) {
472 ret.status.setStatus(0);
473 return ret;
474 }
475
476 double Z2massCut= 2;
477
478 if ( Z1Candidate.M() > 120 ||
479 Z2Candidate.M() < Z2massCut ||
480 Z2Candidate.M() > 120 ||
481 !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec->Pt() > 20.0) ||
482 !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec->Pt() > 10.0)
483 ) {
484 evtfail |= (1<<EVTFAIL_KINEMATICS );
485 // ret.status = evtfail;
486 ret.status.setStatus(0);
487 return ret;
488 }
489
490 unsigned channel;
491 if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
492 if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
493 if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) ||
494 (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
495
496
497
498
499 if( ctrl.debug ) cout << "run: " << info->runNum
500 << "\tevt: " << info->evtNum
501 << "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
502 << "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
503 << "\tmZ1: " << Z1Candidate.M()
504 << "\tmZ2: " << Z2Candidate.M()
505 << "\tm4l: " << ZZSystem.M()
506 << "\tevtfail: " << hex << evtfail << dec
507 << "\ttrigbits: " << hex << info->triggerBits << dec
508 // << "\ttree: " << inputFiles[q][f]
509 << endl;
510
511
512
513 //***************************************************************
514 // finish
515 //***************************************************************
516
517 if( !evtfail ) {
518 ret.status.setStatus(SelectionStatus::EVTPASS);
519 ret.Z1leptons.push_back(lepvec[Z1LeptonMinusIndex]);
520 ret.Z1leptons.push_back(lepvec[Z1LeptonPlusIndex]);
521 ret.Z2leptons.push_back(lepvec[Z2LeptonMinusIndex]);
522 ret.Z2leptons.push_back(lepvec[Z2LeptonPlusIndex]);
523 }
524
525
526
527 return ret;
528 }
529
530