ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/Selection.cc
Revision: 1.11
Committed: Mon Oct 17 16:23:43 2011 UTC (13 years, 7 months ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.10: +75 -16 lines
Log Message:
*** empty log message ***

File Contents

# Content
1 #include "Selection.h"
2 #include "PassHLT.h"
3
4 #include "SiMVAElectronSelection.h"
5
6 #include "HZZCiCElectronSelection.h"
7 #include "HZZLikelihoodElectronSelection.h"
8 #include "HZZBDTElectronSelection.h"
9 #include "RunLumiRangeMap.h"
10
11 RunLumiRangeMap rlrm;
12
13 void initRunLumiRangeMap() {
14 rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
15 // rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
16 rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
17 rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));
18 rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));
19 };
20
21 unsigned fails_HZZ4L_selection(ControlFlags &ctrl, // input control
22 mithep::TEventInfo *info, // input event inof
23 TClonesArray *electronArr, // input electrons
24 TClonesArray *muonArr, // input muons
25 double eventweight, // weight
26 TTree * passtuple ) {
27
28 fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, passtuple, NULL );
29
30 };
31
32 unsigned fails_HZZ4L_selection(ControlFlags &ctrl, // input control
33 mithep::TEventInfo *info, // input event inof
34 TClonesArray *electronArr, // input electrons
35 TClonesArray *muonArr, // input muons
36 double eventweight, // weight
37 LabVectors *l ) {
38
39 fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, NULL, l );
40
41 };
42
43
44 unsigned fails_HZZ4L_selection(ControlFlags &ctrl, // input control
45 mithep::TEventInfo *info, // input event inof
46 TClonesArray *electronArr, // input electrons
47 TClonesArray *muonArr, // input muons
48 double eventweight, // weight
49 TTree * passtuple,
50 LabVectors * l) { // output ntuple
51
52 unsigned evtfail = 0x0;
53
54
55 if( ctrl.debug ) {
56 cout << "Run: " << info->runNum
57 << "\tEvt: " << info->evtNum
58 << "\tLumi: " << info->lumiSec
59 << endl;
60 }
61
62 if( !ctrl.mc ) {
63 // not accounting for overlap atm
64 RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);
65 if( !(rlrm.HasRunLumi(rl)) ) {
66 if( ctrl.debug ) cout << "\tfails JSON" << endl;
67 evtfail |= (1<<EVTFAIL_JSON);
68 return evtfail;
69 }
70 }
71
72
73
74
75
76 //********************************************************
77 // Trigger
78 //********************************************************
79 if( !ctrl.mc ) {
80 // if( !(passHLT(info->triggerBits, info->runNum, channel) ) ) {
81 if( !(passHLT(info->triggerBits, info->runNum, 999) ) ) {
82 evtfail |= (1<<EVTFAIL_TRIGGER);
83 return evtfail;
84 }
85 } else {
86 if( !(passHLTMC(info->triggerBits)) ) {
87 evtfail |= (1<<EVTFAIL_TRIGGER);
88 return evtfail;
89 }
90 // cout << "MC trigger bits: " << hex << info->triggerBits << dec << endl;
91 }
92
93 if( ctrl.debug ) {
94 cout << "presel nlep: " << muonArr->GetEntries() + electronArr->GetEntries()
95 << "\tnmuon: " << muonArr->GetEntries()
96 << "\tnelectron: " << electronArr->GetEntries()
97 << endl;
98 }
99
100 //********************************************************
101 // Lepton Selection
102 //********************************************************
103 vector<SimpleLepton> lepvec;
104
105 //
106 if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
107 //----------------------------------------------------
108 for(Int_t i=0; i<muonArr->GetEntries(); i++) {
109 const mithep::TMuon *mu = (mithep::TMuon*)((*muonArr)[i]);
110 unsigned muonfail;
111 if( ctrl.muSele == "ksWW" )
112 muonfail = passMuonSelection(mu);
113 else
114 muonfail = passMuonSelectionZZ(mu);
115 if( ctrl.debug ) {
116 cout << "muon:: pt: " << mu->pt
117 << "\teta: " << mu->eta
118 << "\tmask: 0x" << hex << muonfail << dec
119 << endl;
120 }
121 if ( !muonfail ) {
122 SimpleLepton tmplep;
123 tmplep.vec.SetPtEtaPhiM(mu->pt,
124 mu->eta,
125 mu->phi,
126 105.658369e-3);
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 lepvec.push_back(tmplep);
139 if( ctrl.debug ) { cout << "muon passes ... " << endl;}
140 }
141 }
142
143 if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
144
145 //----------------------------------------------------
146 for(Int_t i=0; i<electronArr->GetEntries(); i++) {
147 const mithep::TElectron *ele = (mithep::TElectron*)((*electronArr)[i]);
148
149 Bool_t isMuonOverlap = kFALSE;
150 for (int k=0; k<lepvec.size(); ++k) {
151 TVector3 tmplep;
152 tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
153 if ( lepvec[k].type == 13 && lepvec[k].vec.Vect().DrEtaPhi(tmplep) < 0.1 ) {
154 if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
155 isMuonOverlap = kTRUE;
156 break;
157 }
158 }
159
160 unsigned FAIL=0, isEleTight=0;
161
162 unsigned failsCIC=0;
163 CICStruct ciccuts_tight, ciccuts_medium, ciccuts_loose;
164 if(ctrl.eleSele=="cic") {
165 if( ctrl.eleSeleScheme == "mediumloose" ) {
166 ciccuts_medium = getCiCCuts("medium");
167 unsigned failsCICMedium = failsCicSelection(ctrl, ele, ciccuts_medium, ctrl.kinematics);
168 ciccuts_loose = getCiCCuts("loose");
169 unsigned failsCICLoose = failsCicSelection(ctrl, ele, ciccuts_loose, ctrl.kinematics);
170 failsCIC = ( failsCICLoose | failsCICMedium );
171 if( !failsCICMedium ) isEleTight=1;
172 }
173 else {
174 ciccuts_tight = getCiCCuts(ctrl.eleSeleScheme);
175 failsCIC = failsCicSelection(ctrl, ele, ciccuts_tight, ctrl.kinematics);
176 if( !failsCIC ) isEleTight=1;
177 }
178 FAIL = failsCIC;
179 }
180
181 LikStruct likcuts;
182 unsigned failsLike=0;
183 if(ctrl.eleSele=="lik") {
184 likcuts = getLikCuts(ctrl.eleSeleScheme);
185 failsLike = failsLikelihoodSelection(ele, likcuts, ctrl.kinematics);
186 FAIL = failsLike;
187 }
188 unsigned failsBDT=0;
189 if(ctrl.eleSele=="bdt") {
190 if( ctrl.eleSeleScheme == "mediumloose" ) {
191 unsigned failsBDTMedium = failsBDTSelection(ctrl,"medium",ele);
192 unsigned failsBDTLoose = failsBDTSelection(ctrl,"loose",ele);
193 failsBDT = ( failsBDTLoose | failsBDTMedium );
194 if( !failsBDTMedium ) isEleTight=1;
195 } else {
196 failsBDT = failsBDTSelection(ctrl,"tight",ele);
197 if( !failsBDT ) isEleTight=1;
198 }
199 FAIL = failsBDT;
200 }
201 unsigned failsSi=0;
202 if(ctrl.eleSele=="si") {
203 failsSi = failsSiMVAElectronSelection(ctrl, ele, 0.95, ctrl.kinematics);
204 FAIL = failsSi;
205 }
206
207
208
209
210 if( ctrl.debug ){
211 cout << "CIC category: " << cicCategory(ele)
212 << "\tlikelihood: " << ele->likelihood
213 << "\tFAIL: 0x" << hex << FAIL << dec
214 << "\tfailsCIC: 0x" << hex << failsCIC << dec
215 << "\tfailsLike: 0x" << hex << failsLike << dec
216 << "\tfailsBDT: 0x" << hex << failsBDT << dec
217 << "\tscEt: " << ele->scEt
218 << "\tscEta: " << ele->scEta
219 << "\tncluster: " << ele->ncluster
220 << endl;
221 }
222 if ( !FAIL && !isMuonOverlap ) {
223 SimpleLepton tmplep;
224 tmplep.vec.SetPtEtaPhiM( ele->pt,
225 ele->eta,
226 ele->phi,
227 0.51099892e-3 );
228 tmplep.type = 11;
229 tmplep.index = i;
230 tmplep.charge = ele->q;
231 tmplep.isoTrk = ele->trkIso03;
232 tmplep.isoEcal = ele->emIso03;
233 tmplep.isoHcal = ele->hadIso03;
234 tmplep.isoPF03 = ele->pfIso03;
235 tmplep.isoPF04 = ele->pfIso04;
236 tmplep.ip3dSig = ele->ip3dSig;
237 tmplep.is4l = false;
238 tmplep.isTight = isEleTight;
239 tmplep.isEB = ele->isEB;
240 lepvec.push_back(tmplep);
241 if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
242 }
243 }
244
245 sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
246
247 int nmu=0, nele=0;
248 for( int i=0; i<lepvec.size(); i++ ) {
249 if( abs(lepvec[i].type) == 11 ) nele++;
250 else nmu++;
251 }
252 if( ctrl.debug ) {
253 cout << "postsel nlep: " << lepvec.size()
254 << "\tnmuon: " << nmu
255 << "\tnelectron: " << nele
256 << endl;
257 }
258
259 //******************************************************************************
260 //Z1 Selection
261 //******************************************************************************
262 int Z1LeptonPlusIndex = -1;
263 int Z1LeptonMinusIndex = -1;
264 double BestZ1Mass = -999;
265 if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
266 for(int i = 0; i < lepvec.size(); ++i) {
267 for(int j = i+1; j < lepvec.size(); ++j) {
268 if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
269 if (!(lepvec[i].vec.Pt() > 20.0 || lepvec[j].vec.Pt() > 20.0)) continue;
270 if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
271 if (!(lepvec[i].vec.Pt() > 10.0 && lepvec[j].vec.Pt() > 10.0)) continue;
272 if( ctrl.debug ) { cout << "\tthe other is > 10 GeV" << endl; }
273 if (lepvec[i].charge == lepvec[j].charge) continue;
274 if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
275 if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;
276 if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
277
278 //Make Z1 hypothesis
279 TLorentzVector leptonPlus, leptonMinus;
280 if ( lepvec[i].charge > 0 ) {
281 leptonPlus = lepvec[i].vec;
282 leptonMinus = lepvec[j].vec;
283 } else {
284 leptonPlus = lepvec[j].vec;
285 leptonMinus = lepvec[i].vec;
286 }
287
288 float tmpZ1Mass = (leptonPlus+leptonMinus).M();
289 if( ctrl.debug ) cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
290 if( tmpZ1Mass > 60 ) {
291 if (fabs(tmpZ1Mass - 91.1876) < fabs(BestZ1Mass - 91.1876)) {
292 BestZ1Mass = tmpZ1Mass;
293 if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass
294 << "\tdM: " << fabs(BestZ1Mass - 91.1876)
295 << endl;
296 if (lepvec[i].charge > 0) {
297 Z1LeptonPlusIndex = i;
298 Z1LeptonMinusIndex = j;
299 } else {
300 Z1LeptonPlusIndex = j;
301 Z1LeptonMinusIndex = i;
302 }
303 }
304 }
305 }
306 }
307 // stop if no Z1 candidate is found
308 if( BestZ1Mass < 0 ) {
309 evtfail |= (1<<EVTFAIL_Z1);
310 return evtfail;
311 }
312 if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum << endl;
313 if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
314 TLorentzVector Z1LeptonPlus = lepvec[Z1LeptonPlusIndex].vec;
315 TLorentzVector Z1LeptonMinus = lepvec[Z1LeptonMinusIndex].vec;
316 TLorentzVector Z1Candidate = Z1LeptonPlus + Z1LeptonMinus;
317 if( l != NULL ) {
318 l->vecz1 = Z1Candidate;
319 l->vecl1p = Z1LeptonPlus;
320 l->vecl1m = Z1LeptonMinus;
321 }
322
323 //******************************************************************************
324 // Z1 + l
325 //******************************************************************************
326 if( lepvec.size() < 3 ) {
327 evtfail |= (1<<EVTFAIL_Z1_PLUSL);
328 return evtfail;
329 }
330
331 //******************************************************************************
332 // 4l/Z2 Selection
333 //******************************************************************************
334 Int_t Z2LeptonPlusIndex = -1;
335 Int_t Z2LeptonMinusIndex = -1;
336 Double_t BestZ2Mass = -1;
337 if( ctrl.debug ) cout << "looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
338 for(int i = 0; i < lepvec.size(); ++i) {
339 if( abs(lepvec[i].type) == 11 &&
340 ctrl.eleSeleScheme == "mediumloose" &&
341 !(lepvec[i].isTight) ) continue;
342
343 for(int j = i+1; j < lepvec.size(); ++j) {
344 if( abs(lepvec[j].type) == 11 &&
345 ctrl.eleSeleScheme == "mediumloose" &&
346 !(lepvec[j].isTight) ) continue;
347
348 // cout << "i: " << i << "\tj: " << j << endl;
349 if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) {
350 // cout << "\ti matches a Z1 index, skipping ..." << endl;
351 continue; //skip Z1 leptons
352 }
353 if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) {
354 // cout << "\tj matches a Z1 index, skipping ..." << endl;
355 continue; //skip Z1 leptons
356 }
357 if (lepvec[i].charge == lepvec[j].charge) {
358 // cout << "\ti and j are same sign, skipping ..." << endl;
359 continue; //require opp sign
360 }
361 if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
362 // cout << "\ti and j are not same flavor, skipping ..." << endl;
363 continue; //require same flavor
364 }
365
366
367 //Make Z2 hypothesis
368 TLorentzVector leptonPlus, leptonMinus;
369
370 if (lepvec[i].charge > 0 ) {
371 leptonPlus = lepvec[i].vec;
372 leptonMinus = lepvec[j].vec;
373 } else {
374 leptonPlus = lepvec[j].vec;
375 leptonMinus = lepvec[i].vec;
376 }
377
378 TLorentzVector dilepton = leptonPlus+leptonMinus;
379 TLorentzVector fourLepton = Z1Candidate + dilepton;
380
381 if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
382 if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;
383
384 if (!(dilepton.M() > 12.0)) continue;
385 if (!(fourLepton.M() > 100.0)) continue;
386
387 /*
388 //for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
389 if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
390 TLorentzVector pair1 = Z1LeptonPlus+leptonMinus;
391 TLorentzVector pair2 = Z1LeptonMinus+leptonPlus;
392 if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
393 if (!(pair1.M() > 12 || pair2.M() > 12)) continue;
394 }
395 */
396
397 //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
398 if (Z2LeptonPlusIndex < 0) {
399 if (lepvec[i].charge > 0) {
400 Z2LeptonPlusIndex = i;
401 Z2LeptonMinusIndex = j;
402 } else {
403 Z2LeptonPlusIndex = j;
404 Z2LeptonMinusIndex = i;
405 }
406 } else {
407 Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec.Pt();
408 Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec.Pt();
409 if (lepvec[Z2LeptonMinusIndex].vec.Pt() > BestPairPtMax) {
410 BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec.Pt();
411 BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec.Pt();
412 }
413
414 Double_t CurrentPairPtMax = lepvec[i].vec.Pt();
415 Double_t CurrentPairPtMin = lepvec[j].vec.Pt();
416 if (lepvec[j].vec.Pt() > CurrentPairPtMax) {
417 CurrentPairPtMax = lepvec[j].vec.Pt();
418 CurrentPairPtMin = lepvec[i].vec.Pt();
419 }
420
421 if (CurrentPairPtMax > BestPairPtMax) {
422 if (lepvec[i].charge > 0) {
423 Z2LeptonPlusIndex = i;
424 Z2LeptonMinusIndex = j;
425 } else {
426 Z2LeptonPlusIndex = j;
427 Z2LeptonMinusIndex = i;
428 }
429 } else if (CurrentPairPtMax == BestPairPtMax) {
430 if (CurrentPairPtMin > BestPairPtMin) {
431 if (lepvec[i].charge > 0) {
432 Z2LeptonPlusIndex = i;
433 Z2LeptonMinusIndex = j;
434 } else {
435 Z2LeptonPlusIndex = j;
436 Z2LeptonMinusIndex = i;
437 }
438 }
439 }
440 }
441 }
442 }
443
444 // stop if no Z2 candidate is found
445 if (Z2LeptonPlusIndex == -1) {
446 evtfail |= ( 1<<EVTFAIL_4L );
447 return evtfail;
448 // h_evtfail->Fill( evtfail );
449 // cout << "evtfail: " << hex << evtfail << dec << endl;
450 // continue;
451 }
452 if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
453 if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
454 << "\tminusindex: " << Z2LeptonMinusIndex << endl;
455 TLorentzVector Z2LeptonPlus = lepvec[Z2LeptonPlusIndex].vec;
456 TLorentzVector Z2LeptonMinus = lepvec[Z2LeptonMinusIndex].vec;
457 TLorentzVector Z2Candidate = Z2LeptonPlus+Z2LeptonMinus;
458 TLorentzVector ZZSystem = Z1Candidate + Z2Candidate;
459 if( l != NULL ) {
460 l->vecz2 = Z2Candidate;
461 l->vecl2p = Z2LeptonPlus;
462 l->vecl2m = Z2LeptonMinus;
463 l->vec4l = ZZSystem;
464 }
465 lepvec[Z1LeptonPlusIndex].is4l = true;
466 lepvec[Z1LeptonMinusIndex].is4l = true;
467 lepvec[Z2LeptonPlusIndex].is4l = true;
468 lepvec[Z2LeptonMinusIndex].is4l = true;
469
470 //***************************************************************
471 // Isolation
472 //***************************************************************
473 bool failiso=false;
474
475 if( ctrl.isoScheme == "pf" ) {
476
477 for( int i=0; i<lepvec.size(); i++ ) {
478
479 if( !(lepvec[i].is4l) ) continue;
480
481 if( abs(lepvec[i].type) == 11 ) {
482 float reliso = lepvec[i].isoPF04/lepvec[i].vec.Pt();
483 if( lepvec[i].isEB && lepvec[i].vec.Pt() > 20 && reliso > PFISO_ELE_LOOSE_EB_HIGHPT ) {
484 failiso = true;
485 break;
486 }
487 if( lepvec[i].isEB && lepvec[i].vec.Pt() < 20 && reliso > PFISO_ELE_LOOSE_EB_LOWPT ) {
488 failiso = true;
489 break;
490 }
491 if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() > 20 && reliso > PFISO_ELE_LOOSE_EE_HIGHPT ) {
492 failiso = true;
493 break;
494 }
495 if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() < 20 && reliso > PFISO_ELE_LOOSE_EE_LOWPT ) {
496 failiso = true;
497 break;
498 }
499 }
500
501 if( abs(lepvec[i].type) == 13 ) {
502
503 float reliso = lepvec[i].isoPF03/lepvec[i].vec.Pt();
504 if( lepvec[i].isEB && lepvec[i].vec.Pt() > 20 && reliso > PFISO_MU_TIGHT_EB_HIGHPT ) { //0.13
505 failiso = true;
506 break;
507 }
508 if( lepvec[i].isEB && lepvec[i].vec.Pt() < 20 && reliso > PFISO_MU_TIGHT_EB_LOWPT ) { //0.06
509 failiso = true;
510 break;
511 }
512 if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() > 20 && reliso > PFISO_MU_TIGHT_EE_HIGHPT ) { //0.09
513 failiso = true;
514 break;
515 }
516 if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() < 20 && reliso > PFISO_MU_TIGHT_EE_LOWPT ) { //0.05
517 failiso = true;
518 break;
519 }
520 }
521 }
522 } else if ( ctrl.isoScheme == "pairwise" ) {
523 float rho = info->rho;
524 for( int i=0; i<lepvec.size(); i++ ) {
525 if( !(lepvec[i].is4l) ) continue;
526 float effArea_ecal_i, effArea_hcal_i;
527 if( lepvec[i].isEB ) {
528 if( lepvec[i].type == 11 ) {
529 effArea_ecal_i = 0.101;
530 effArea_hcal_i = 0.021;
531 } else {
532 effArea_ecal_i = 0.074;
533 effArea_hcal_i = 0.022;
534 }
535 } else {
536 if( lepvec[i].type == 11 ) {
537 effArea_ecal_i = 0.046;
538 effArea_hcal_i = 0.040;
539 } else {
540 effArea_ecal_i = 0.045;
541 effArea_hcal_i = 0.030;
542 }
543 }
544 float isoEcal_corr_i = lepvec[i].isoEcal - (effArea_ecal_i*rho);
545 float isoHcal_corr_i = lepvec[i].isoHcal - (effArea_hcal_i*rho);
546 for( int j=i+1; j<lepvec.size(); j++ ) {
547 if( !(lepvec[j].is4l) ) continue;
548 float effArea_ecal_j, effArea_hcal_j;
549 if( lepvec[j].isEB ) {
550 if( lepvec[j].type == 11 ) {
551 effArea_ecal_j = 0.101;
552 effArea_hcal_j = 0.021;
553 } else {
554 effArea_ecal_j = 0.074;
555 effArea_hcal_j = 0.022;
556 }
557 } else {
558 if( lepvec[j].type == 11 ) {
559 effArea_ecal_j = 0.046;
560 effArea_hcal_j = 0.040;
561 } else {
562 effArea_ecal_j = 0.045;
563 effArea_hcal_j = 0.030;
564 }
565 }
566 float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
567 float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
568 float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec.Pt();
569 float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec.Pt();
570 float comboIso = RIso_i + RIso_j;
571 if( info->evtNum == 1038911933 ) {
572 float tmpdR = lepvec[i].vec.DrEtaPhi(lepvec[j].vec);
573 cout << "i: " << i
574 << "\tdR: " << tmpdR
575 << "\trho: " << rho
576 << "\tRIso_i: " << RIso_i
577 << "\ttkrel: " << lepvec[i].isoTrk/lepvec[i].vec.Pt()
578 << "\tecalrel: " << lepvec[i].isoEcal/lepvec[i].vec.Pt()
579 << "\tecalrelcor: " << isoEcal_corr_i/lepvec[i].vec.Pt()
580 << "\thcalrel: " << lepvec[i].isoHcal/lepvec[i].vec.Pt()
581 << "\thcalrelcor: " << isoHcal_corr_i/lepvec[i].vec.Pt()
582 << "\tpt_i: " << lepvec[i].vec.Pt()
583 << "\tj: " << j
584 << "\tRIso_j: " << RIso_j
585 << "\ttkrel: " << lepvec[j].isoTrk/lepvec[j].vec.Pt()
586 << "\tecalrel: " << lepvec[j].isoEcal/lepvec[j].vec.Pt()
587 << "\tecalrelcor: " << isoEcal_corr_j/lepvec[j].vec.Pt()
588 << "\thcalrel: " << lepvec[j].isoHcal/lepvec[j].vec.Pt()
589 << "\thcalrelcor: " << isoHcal_corr_j/lepvec[j].vec.Pt()
590 << "\tpt_j: " << lepvec[j].vec.Pt()
591 << "\tcombo: " << comboIso
592 << endl;
593 cout.flush();
594 }
595 if( comboIso > 0.35 ) {
596 if( ctrl.debug ) cout << "combo failing for indices: " << i << "," << j << endl;
597 failiso = true;
598 // break;
599 }
600 }
601 }
602 }
603
604
605 if( failiso ) {
606 evtfail |= ( 1<<EVTFAIL_ISOLATION );
607 return evtfail;
608 //h_evtfail->Fill( evtfail, eventweight );
609 // h_evtfail->Fill( evtfail );
610 // cout << "evtfail: " << hex << evtfail << dec << endl;
611 // continue;
612 }
613
614 //***************************************************************
615 // IP significance
616 //***************************************************************
617 bool failip = false;
618 for( int i=0; i<lepvec.size(); i++ ) {
619 if( !(lepvec[i].is4l) ) continue;
620 if( lepvec[i].ip3dSig > 4 ) {
621 failip=true;
622 break;
623 }
624 }
625 if( failip ) {
626 evtfail |= (1<<EVTFAIL_IP );
627 return evtfail;
628 //h_evtfail->Fill( evtfail, eventweight );
629 // h_evtfail->Fill( evtfail );
630 // cout << "evtfail: " << hex << evtfail << dec << endl;
631 // continue;
632 }
633
634 //***************************************************************
635 // remaining kinematic cuts
636 //***************************************************************
637 double Z2massCut=0;
638 if ( ctrl.kinematics == "loose" ) Z2massCut = 12;
639 else if ( ctrl.kinematics == "tight" ) Z2massCut = 20;
640 else { cout << "error! kinematic tightness not defined!" << endl; assert(0); }
641
642 if ( Z1Candidate.M() > 120 ||
643 Z2Candidate.M() < Z2massCut ||
644 Z2Candidate.M() > 120 ||
645 !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec.Pt() > 20.0) ||
646 !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec.Pt() > 10.0)
647 ) {
648 evtfail |= (1<<EVTFAIL_KINEMATICS );
649 return evtfail;
650 //h_evtfail->Fill( evtfail, eventweight );
651 // h_evtfail->Fill( evtfail );
652 // cout << "evtfail: " << hex << evtfail << dec << endl;
653 // continue;
654 }
655
656 int channel;
657 if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
658 if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
659 if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) ||
660 (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
661
662
663
664 if( passtuple != NULL ) {
665 unsigned run = info->runNum;
666 unsigned evt = info->evtNum;
667 unsigned lumi = info->lumiSec;
668 unsigned chan = channel;
669 double w = eventweight;
670 float mZ1 = Z1Candidate.M() ;
671 float mZ2 = Z2Candidate.M() ;
672 float m4l = ZZSystem.M() ;
673 float pt4l = ZZSystem.Pt() ;
674 unsigned tZ1 = abs(lepvec[Z1LeptonMinusIndex].type);
675 unsigned tZ2 = abs(lepvec[Z2LeptonMinusIndex].type);
676 passtuple->SetBranchAddress("run", &run);
677 passtuple->SetBranchAddress("evt", &evt);
678 passtuple->SetBranchAddress("lumi", &lumi);
679 passtuple->SetBranchAddress("mZ1", &mZ1);
680 passtuple->SetBranchAddress("mZ2", &mZ2);
681 passtuple->SetBranchAddress("tZ1", &tZ1);
682 passtuple->SetBranchAddress("tZ2", &tZ2);
683 passtuple->SetBranchAddress("m4l", &m4l);
684 passtuple->SetBranchAddress("pt4l", &pt4l);
685 passtuple->SetBranchAddress("w", &w);
686 passtuple->Fill( );
687 }
688
689 if( ctrl.debug ) cout << "run: " << info->runNum
690 << "\tevt: " << info->evtNum
691 << "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
692 << "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
693 << "\tmZ1: " << Z1Candidate.M()
694 << "\tmZ2: " << Z2Candidate.M()
695 << "\tm4l: " << ZZSystem.M()
696 << "\tevtfail: " << hex << evtfail << dec
697 << "\ttrigbits: " << hex << info->triggerBits << dec
698 // << "\ttree: " << inputFiles[q][f]
699 << endl;
700
701 return evtfail;
702
703 }
704
705
706
707
708