ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/ReferenceSelection.cc
Revision: 1.17
Committed: Fri Oct 12 12:26:12 2012 UTC (12 years, 7 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.16: +14 -12 lines
Log Message:
updated selection with regression

File Contents

# Content
1 //***************************************************************************************************
2 //
3 // selection sync'ed with https://twiki.cern.ch/twiki/bin/viewauth/CMS/HiggsZZ4l2012SummerSync
4 //
5 //***************************************************************************************************
6
7 // system headers
8 #include <map>
9 #include <utility>
10
11 // mit headers
12 #include "Vertex.h"
13
14 // 4L stuff
15 #include "SelectionStatus.h"
16 #include "EventData.h"
17 #include "SimpleLepton.h"
18 #include "EfficiencyWeightsInterface.h"
19 #include "ElectronSelection.h"
20 #include "MuonSelection.h"
21 #include "IsolationSelection.h"
22 #include "ReferenceSelection.h"
23 #include "Selection.h"
24 #include "CommonDefs.h"
25 #include "SelectionDefs.h"
26 #include "FSR.h"
27 #include "SelectionFuncs.h"
28 #include "ElectronMomentumCorrection.h"
29 #include "MuonMomentumCorrection.h"
30
31
32 extern vector<SimpleLepton> failingLeptons;
33 extern vector<SimpleLepton> passingLeptons;
34
35 extern vector<unsigned> cutvec;
36 extern vector<vector<unsigned> > zcutvec;
37 extern vector<vector<unsigned> > zzcutvec;
38 extern map<unsigned,float> evtrhoMap;
39 extern bool passes_HLT_MC;
40
41 //
42 // prototypes
43 //--------------------------------------------------------------------------------------------------
44 void updateSimpleLepton(SimpleLepton &tmplep);
45 void fillVetoArrays( ControlFlags & ctrl,
46 const mithep::Array<mithep::Muon> *muonArr,
47 vector< const mithep::Muon*> & muonsToVeto,
48 const mithep::Array<mithep::Electron> *electronArr,
49 vector< const mithep::Electron*> & electronsToVeto,
50 const mithep::Vertex * vtx ) ;
51 //--------------------------------------------------------------------------------------------------
52
53 //--------------------------------------------------------------------------------------------------
54 EventData apply_HZZ4L_reference_selection(ControlFlags &ctrl, // input control
55 const mithep::EventHeader *info, // input event info
56 const mithep::Array<mithep::Vertex> * vtxArr ,
57 const mithep::Array<mithep::PFCandidate> *pfCandidates,
58 const mithep::Array<mithep::PileupEnergyDensity> *puEnergyDensity,
59 const mithep::Array<mithep::Electron> *electronArr, // input electrons
60 SelectionStatus (*ElectronPreSelector)( ControlFlags &,
61 const mithep::Electron*,
62 const mithep::Vertex *),
63 SelectionStatus (*ElectronIDSelector)( ControlFlags &,
64 const mithep::Electron*,
65 const mithep::Vertex *),
66 SelectionStatus (*ElectronIsoSelector)( ControlFlags &,
67 const mithep::Electron*,
68 const mithep::Vertex *,
69 const mithep::Array<mithep::PFCandidate> *,
70 const mithep::Array<mithep::PileupEnergyDensity> *,
71 mithep::ElectronTools::EElectronEffectiveAreaTarget,
72 vector<const mithep::PFCandidate*>),
73 const mithep::Array<mithep::Muon> *muonArr, // input muons
74 SelectionStatus (*MuonPreSelector)( ControlFlags &,
75 const mithep::Muon*,
76 const mithep::Vertex *,
77 const mithep::Array<mithep::PFCandidate> *),
78 SelectionStatus (*MuonIDSelector)( ControlFlags &,
79 const mithep::Muon*,
80 // const mithep::Vertex &),
81 const mithep::Vertex *,
82 const mithep::Array<mithep::PFCandidate> *),
83 SelectionStatus (*MuonIsoSelector)( ControlFlags &,
84 const mithep::Muon*,
85 const mithep::Vertex *,
86 const mithep::Array<mithep::PFCandidate> *,
87 const mithep::Array<mithep::PileupEnergyDensity> *,
88 mithep::MuonTools::EMuonEffectiveAreaTarget,
89 vector<const mithep::PFCandidate*>)
90 )
91 //--------------------------------------------------------------------------------------------------
92 {
93
94 EventData ret;
95 unsigned evtfail = 0x0;
96 TRandom3 r;
97
98 failingLeptons.clear();
99 passingLeptons.clear();
100
101 mithep::MuonTools::EMuonEffectiveAreaTarget eraMu;
102 mithep::ElectronTools::EElectronEffectiveAreaTarget eraEle;
103 getEATargets(ctrl,eraMu,eraEle);
104
105 if( ctrl.debug ) {
106 cout << "presel nlep: " << muonArr->GetEntries() + electronArr->GetEntries()
107 << "\tnmuon: " << muonArr->GetEntries()
108 << "\tnelectron: " << electronArr->GetEntries()
109 << endl;
110 }
111
112 //correct muon momentum
113 if(ctrl.correct_muon_momentum){
114 for(int i=0; i<muonArr->GetEntries(); i++) {
115 const mithep::Muon *const_mu = (mithep::Muon*)((*muonArr)[i]);
116 mithep::Muon *mu = const_cast<mithep::Muon*>(const_mu);
117 correct_muon_momentum(mu,info->RunNum());
118 }
119 }
120
121 //********************************************************
122 // Skim 0 :
123 //********************************************************
124 int nlep_above_10=0;
125 int nlep_above_20=0;
126 for(int i=0; i<muonArr->GetEntries(); i++)
127 {
128 const mithep::Muon *mu = (mithep::Muon*)((*muonArr)[i]);
129
130
131 if( !(mu->IsTrackerMuon() || mu->IsGlobalMuon()) ) continue;
132 if( fabs(mu->Eta()) > 2.4 ) continue;
133 if( mu->Pt() > 10 ) nlep_above_10++;
134 if( mu->Pt() > 20 ) nlep_above_20++;
135 }
136 for(int i=0; i<electronArr->GetEntries(); i++)
137 {
138 const mithep::Electron *ele = (mithep::Electron*)((*electronArr)[i]);
139 if( fabs(ele->Eta()) > 2.5 ) continue;
140 if( ele->Pt() > 10 ) nlep_above_10++;
141 if( ele->Pt() > 20 ) nlep_above_20++;
142 }
143 if( (nlep_above_10 > 1 && nlep_above_20 > 0) ) {
144 ret.status.selectionBits.flip(PASS_SKIM0);
145 cutvec[PASS_SKIM0] +=1;
146 } else {
147 ret.status.setStatus(SelectionStatus::FAIL);
148 return ret;
149 }
150
151 //********************************************************
152 // Skim 0.1 : 1 SF pair with mLL > 40
153 //********************************************************
154 bool ossf_pair=false;
155 for(int i=0; i<muonArr->GetEntries(); i++)
156 {
157 const mithep::Muon *mu1 = (mithep::Muon*)((*muonArr)[i]);
158 if( !(mu1->IsTrackerMuon() || mu1->IsGlobalMuon()) ) continue;
159 if( fabs(mu1->Eta()) > 2.4 ) continue;
160 if( fabs(mu1->Pt()) < 3 ) continue;
161 for(int j=i+1; j<muonArr->GetEntries(); j++)
162 {
163 const mithep::Muon *mu2 = (mithep::Muon*)((*muonArr)[j]);
164 if( !(mu2->IsTrackerMuon() || mu2->IsGlobalMuon()) ) continue;
165 if( fabs(mu2->Eta()) > 2.4 ) continue;
166 if( fabs(mu2->Pt()) < 3 ) continue;
167 TLorentzVector v1,v2;
168 v1.SetPtEtaPhiM( mu1->Pt(), mu1->Eta(), mu1->Phi(), MUON_MASS);
169 v2.SetPtEtaPhiM( mu2->Pt(), mu2->Eta(), mu2->Phi(), MUON_MASS);
170 if( (v1+v2).M() >= 40) ossf_pair = true;
171 }
172 }
173 for(int i=0; i<electronArr->GetEntries(); i++)
174 {
175 const mithep::Electron *el1 = (mithep::Electron*)((*electronArr)[i]);
176 if( fabs(el1->Eta()) > 2.5 ) continue;
177 if( el1->Pt() < 5 ) continue;
178 for(int j=i+1; j<electronArr->GetEntries(); j++)
179 {
180 const mithep::Electron *el2 = (mithep::Electron*)((*electronArr)[j]);
181 if( fabs(el2->Eta()) > 2.5 ) continue;
182 if( el2->Pt() < 5 ) continue;
183 TLorentzVector v1,v2;
184 v1.SetPtEtaPhiM( el1->Pt(), el1->Eta(), el1->Phi(), ELECTRON_MASS);
185 v2.SetPtEtaPhiM( el2->Pt(), el2->Eta(), el2->Phi(), ELECTRON_MASS);
186 if( (v1+v2).M() >= 40 ) ossf_pair = true;
187 }
188 }
189 if( ossf_pair ) {
190 ret.status.selectionBits.flip(PASS_SKIM1);
191 cutvec[PASS_SKIM1] +=1;
192 } else {
193 ret.status.setStatus(SelectionStatus::FAIL);
194 return ret;
195 }
196
197 const mithep::Vertex * vtx;
198 bool goodVertex = setPV( ctrl, vtxArr, vtx );
199 if(goodVertex) {
200 ret.status.selectionBits.flip(PASS_SKIM2);
201 cutvec[PASS_SKIM2] +=1;
202 } else {
203 if(ctrl.debug) cout << "found bad vertex" << endl;
204 ret.status.setStatus(SelectionStatus::FAIL);
205 return ret;
206 }
207 if(ctrl.debug) {
208 cerr << "vtx :: ntrks: " << vtx->NTracksFit() << endl;
209 cerr.flush();
210 }
211
212 //***********************************************************
213 // Trigger Selection
214 //***********************************************************
215 if( ctrl.mc ) {
216 if( passes_HLT_MC ) {
217 ret.status.selectionBits.flip(PASS_TRIGGER);
218 cutvec[PASS_TRIGGER] +=1;
219 } else {
220 ret.status.setStatus(SelectionStatus::FAIL);
221 return ret;
222 }
223 } else {
224 ret.status.selectionBits.flip(PASS_TRIGGER);
225 cutvec[PASS_TRIGGER] +=1;
226 }
227
228 //***********************************************************
229 // Lepton Selection
230 //***********************************************************
231 vector<SimpleLepton> lepvec;
232
233 // empty, reference is not applying additional vetos
234 // vector<const mithep::Muon*> muonsToVeto;
235 // vector<const mithep::Electron*> electronsToVeto;
236 // fillVetoArrays( ctrl, muonArr, muonsToVeto, electronArr, electronsToVeto, vtx );
237 vector<const mithep::PFCandidate*> photonsToVeto;
238
239
240 //
241 if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
242 //----------------------------------------------------
243 for(int i=0; i<muonArr->GetEntries(); i++)
244 {
245 const mithep::Muon *mu = (mithep::Muon*)((*muonArr)[i]);
246
247 SelectionStatus musel;
248 if(ctrl.debug) cout << "musel.status before anything: " << musel.getStatus() << endl;
249
250 musel |= (*MuonPreSelector)(ctrl,mu,vtx,pfCandidates);
251 if(ctrl.debug) cout << "musel.status after presel: " << musel.getStatus() << endl;
252 if( !(musel.getStatus() & SelectionStatus::PRESELECTION) ) continue;
253 if( ctrl.debug ) cout << endl;
254
255 musel |= (*MuonIDSelector)(ctrl,mu,vtx,pfCandidates );
256 if(ctrl.debug) cout << "musel.status after ID: " << musel.getStatus() << endl;
257 if( ctrl.debug ) cout << endl;
258
259 if( !(ctrl.doFSR) ) {
260 musel |= (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto);
261 if(ctrl.debug) cout << "musel.status after iso: " << musel.getStatus() << endl;
262 if( ctrl.debug ) cout << "isomva : " << musel.isoMVA << endl;
263 if( ctrl.debug ) cout << endl;
264 }
265
266 if( ctrl.debug ) {
267 cout << "muon:: pt: " << mu->Pt()
268 << "\teta: " << mu->Eta()
269 << "\tstatus: " << hex << musel.getStatus() << dec
270 << endl;
271 }
272
273
274 SimpleLepton tmplep;
275 float pt = mu->Pt();
276 tmplep.vec.SetPtEtaPhiM(pt,
277 mu->Eta(),
278 mu->Phi(),
279 MUON_MASS);
280
281 tmplep.type = 13;
282 tmplep.index = i;
283 tmplep.charge = mu->Charge();
284 tmplep.isoTrk = mu->IsoR03SumPt();
285 tmplep.isoEcal = mu->IsoR03EmEt();
286 tmplep.isoHcal = mu->IsoR03HadEt();
287 tmplep.isoPF04 = musel.isoPF04;
288 tmplep.chisoPF04 = musel.chisoPF04;
289 tmplep.gaisoPF04 = musel.gaisoPF04;
290 tmplep.neisoPF04 = musel.neisoPF04;
291 // tmplep.isoPF03 = computePFMuonIso(mu,vtx,pfCandidates,0.3);
292 // tmplep.isoPF04 = computePFMuonIso(mu,vtx,pfCandidates,0.4);
293 tmplep.ip3dSig = mu->Ip3dPVSignificance();
294 tmplep.is4l = false;
295 tmplep.isEB = (fabs(mu->Eta()) < 1.479 ? 1 : 0 );
296 tmplep.isoMVA = musel.isoMVA;
297 tmplep.isTight = musel.tight();
298 tmplep.isLoose = musel.loose();
299 tmplep.status = musel;
300 tmplep.fsrRecoveryAttempted = false;
301 lepvec.push_back(tmplep);
302 if( ctrl.debug ) cout << endl;
303 }
304
305
306
307 //
308 if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
309 // --------------------------------------------------------------------------------
310 for(int i=0; i<electronArr->GetEntries(); i++)
311 {
312 // const mithep::Electron *ele = (mithep::Electron*)((*electronArr)[i]);
313 const mithep::Electron *const_ele = (mithep::Electron*)((*electronArr)[i]);
314 mithep::Electron *ele = const_cast<mithep::Electron*>(const_ele);
315 // if(ctrl.do_escale) correct_electron_energy(ele);
316 // if(ctrl.do_eregression) apply_electron_energy_regression(ele,puEnergyDensity->At(0)->RhoKt6PFJets(),vtxArr->GetEntries());
317
318 SelectionStatus elesel;
319 if( ctrl.debug ) cout << endl;
320 if( ctrl.debug ) cout << "--> status before anything: " << hex << elesel.getStatus() << dec << endl;
321
322 elesel |= (*ElectronIDSelector)(ctrl,ele,vtx);
323 if( ctrl.debug ) cout << "--> status after ID: " << hex << elesel.getStatus() << dec << endl;
324 if( ctrl.debug ) cout << endl;
325
326 ElectronMomentumCorrection electron_momentum_correction;
327 electron_momentum_correction.correct_electron_momentum(ctrl, ele,puEnergyDensity->At(0)->Rho(),vtxArr->GetEntries());
328
329 elesel |= (*ElectronPreSelector)(ctrl,ele,vtx);
330 if( ctrl.debug ) cout << "--> status after presel: " << hex << elesel.getStatus() << dec << endl;
331 if( !(elesel.getStatus() & SelectionStatus::PRESELECTION) ) continue;
332 if( ctrl.debug ) cout << endl;
333
334
335 if( !(ctrl.doFSR) ) {
336 elesel |= (*ElectronIsoSelector)(ctrl,ele,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto);
337 if( ctrl.debug ) cout << "--> status after iso: " << hex << elesel.getStatus() << dec << endl;
338 if( ctrl.debug ) cout << endl;
339 }
340
341 if( ctrl.debug ){
342 cout << "\tscEt: " << ele->SCluster()->Et()
343 << "\tscEta: " << ele->SCluster()->Eta()
344 << "\tstatus: " << hex << elesel.getStatus() << dec
345 << endl;
346 }
347
348
349 SimpleLepton tmplep;
350 float pt = ele->Pt();
351 tmplep.vec.SetPtEtaPhiM( pt,
352 ele->Eta(),
353 ele->Phi(),
354 ELECTRON_MASS );
355
356 tmplep.type = 11;
357 tmplep.index = i;
358 tmplep.charge = ele->Charge();
359 tmplep.isoTrk = ele->TrackIsolationDr03();
360 tmplep.isoEcal = ele->EcalRecHitIsoDr03();
361 tmplep.isoHcal = ele->HcalTowerSumEtDr03();
362 tmplep.isoPF04 = elesel.isoPF04;
363 tmplep.chisoPF04 = elesel.chisoPF04;
364 tmplep.gaisoPF04 = elesel.gaisoPF04;
365 tmplep.neisoPF04 = elesel.neisoPF04;
366 // tmplep.isoPF03 = computePFEleIso(ele,vtx,pfCandidates,0.3);
367 // tmplep.isoPF04 = computePFEleIso(ele,vtx,pfCandidates,0.4);
368 tmplep.ip3dSig = ele->Ip3dPVSignificance();
369 tmplep.is4l = false;
370 tmplep.isEB = ele->IsEB();
371 tmplep.scID = ele->SCluster()->GetUniqueID();
372 tmplep.isTight = elesel.tight();
373 tmplep.isLoose = elesel.loose();
374 tmplep.status = elesel;
375 tmplep.idMVA = elesel.idMVA;
376 tmplep.isoMVA = elesel.isoMVA;
377 tmplep.fsrRecoveryAttempted = false;
378 lepvec.push_back(tmplep);
379 if( ctrl.debug ) cout << endl;
380 }
381
382 //********************************************************
383 // Dump Stuff
384 //********************************************************
385 sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
386 // sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_reversept_sort );
387 int nmu=0, nele=0;
388 for( int i=0; i<lepvec.size(); i++ ) {
389 if(ctrl.debug) {
390 cout << "lepvec :: evt: " << info->EvtNum()
391 << "\tindex: " << i
392 << "\ttype: " << lepvec[i].type
393 << "\tpt: " << lepvec[i].vec.Pt()
394 << "\teta: " << lepvec[i].vec.Eta();
395 if( abs(lepvec[i].type) == 11 ) {
396 const mithep::Electron *tmpele = (mithep::Electron*)((*electronArr)[lepvec[i].index]);
397 cout << "\tSCeta: " << tmpele->SCluster()->Eta()
398 << "\tidMVA: " << lepvec[i].idMVA;
399 }
400 cout << "\tloose: " << lepvec[i].isLoose
401 << "\tpf: " << lepvec[i].isoPF04
402 << "\tch: " << lepvec[i].chisoPF04
403 << "\tga: " << lepvec[i].gaisoPF04
404 << "\tne: " << lepvec[i].neisoPF04 ;
405 cout << endl;
406 }
407 if( abs(lepvec[i].type) == 11 ) nele++;
408 else nmu++;
409 }
410 if( ctrl.debug ) {
411 cout << "postsel nlep: " << lepvec.size()
412 << "\tnmuon: " << nmu
413 << "\tnelectron: " << nele
414 << endl;
415 }
416
417
418
419
420 //********************************************************
421 // Step 2: Lepton Cleaning
422 //********************************************************
423 vector<vector<SimpleLepton>::iterator> electrons_to_erase;
424 for (vector<SimpleLepton>::iterator it1=lepvec.begin();
425 it1 != lepvec.end(); it1++ ) {
426 if ( abs(it1->type) != 11 ) continue;
427 TVector3 evec = it1->vec.Vect();
428
429 bool erase_this_electron=false;
430 for (vector<SimpleLepton>::iterator it2=lepvec.begin();
431 it2 != lepvec.end(); it2++ ) {
432 if ( it2 == it1 ) continue;
433 if ( abs(it2->type) != 13 ) continue;
434 if( !(it2->status.looseIDAndPre()) ) continue;
435 TVector3 mvec = it2->vec.Vect();
436
437 if ( evec.DrEtaPhi(mvec) < 0.05 ) {
438 erase_this_electron=true;
439 break;
440 }
441 }
442 if( erase_this_electron ) {
443 if( ctrl.debug ) cout << "erasing electron with pt " << it1->vec.Pt() << endl;
444 electrons_to_erase.push_back(it1);
445 }
446 }
447 for( int i=0; i<electrons_to_erase.size(); i++ ) {
448 lepvec.erase(electrons_to_erase[i]);
449 }
450 if( lepvec.size() >= 4 ) {
451 // ret.status.selectionBits.flip(3);
452 } else {
453 // ret.status.setStatus(SelectionStatus::FAIL);
454 // return ret;
455 }
456
457
458
459 //********************************************************
460 // Step 3: Good Leptons
461 //********************************************************
462 vector<double> pt_of_leptons_to_erase;
463 for (int i=0; i<lepvec.size(); i++ ) {
464 bool already_pushed=false;
465 if( !(lepvec[i].status.looseIDAndPre()) ) {
466 pt_of_leptons_to_erase.push_back(lepvec[i].vec.Pt());
467 already_pushed = true;
468 if(ctrl.debug)
469 cout << "pushing failed lepton type: " << lepvec[i].type
470 << "\tpt: " << lepvec[i].vec.Pt()
471 << "\teta: " << lepvec[i].vec.Eta()
472 << endl;
473 failingLeptons.push_back(lepvec[i]); // these should pass preselection
474 } else {
475 passingLeptons.push_back(lepvec[i]);
476 }
477 if( !already_pushed && fabs(lepvec[i].ip3dSig)>4 )
478 pt_of_leptons_to_erase.push_back(lepvec[i].vec.Pt());
479 }
480 for( int i=0; i<pt_of_leptons_to_erase.size(); i++ ) {
481 for( vector<SimpleLepton>::iterator it=lepvec.begin();
482 it != lepvec.end(); it++ ) {
483 SimpleLepton flep = *it;
484 if( flep.vec.Pt() != pt_of_leptons_to_erase[i] ) continue;
485 if(ctrl.debug) cout << "erasing lepton : "
486 << flep.vec.Pt() << "\t"
487 << flep.type << "\t"
488 << endl;
489 lepvec.erase(it);
490 break;
491 }
492 }
493 if( ctrl.debug ) cout << "good leptons : " << lepvec.size() << endl;
494
495 //********************************************************
496 // Step 4: Z candidate preselection
497 //********************************************************
498 std::vector<std::pair<int,int> > ZCandidates;
499 std::vector<std::pair<SimpleLepton,SimpleLepton> > ZCandidatesLeptons;
500 std::vector<std::pair<SelectionStatus,SelectionStatus> > ZCandidatesSelectionStatus;
501 for(int i = 0; i < lepvec.size(); ++i) {
502 for(int j = i+1; j < lepvec.size(); ++j) {
503 if( abs(lepvec[i].type) != abs(lepvec[j].type) ) continue;
504 if( lepvec[i].charge == lepvec[j].charge ) continue;
505
506 TLorentzVector Zvec = (lepvec[i].vec+lepvec[j].vec);
507
508 vector<SimpleLepton> lepvec_i = lepvec;
509 vector<SimpleLepton> lepvec_j = lepvec;
510
511 if( ctrl.doFSR ) {
512 cout << endl;
513 cout << "----------------> FSR ("<<i<<","<<j<<") <----------------------" << endl;
514 photonsToVeto.clear();
515 float old_pt_i = lepvec[i].vec.Pt();
516 float old_pt_j = lepvec[j].vec.Pt();
517 float old_M = Zvec.M();
518
519
520 cout << "i: " << i << endl;
521 if( abs(lepvec[i].type) == 13 ) {
522 const mithep::Muon *mu = (mithep::Muon*)((*muonArr)[lepvec[i].index]);
523 mithep::Muon * newmu = const_cast<mithep::Muon *>(mu);
524 if( recover_typeI_Photon( ctrl, newmu, i, lepvec_i, pfCandidates, electronArr, &Zvec, photonsToVeto ) ) {
525 if(ctrl.debug) cout << "FSR TYPEI :: oldpt: " << old_pt_i
526 << "\tnewpt: " << lepvec_i[i].vec.Pt()
527 << "\tindex: " << i
528 << endl;
529 //lepvec[i].fsrRecoveryAttempted=true;
530 }
531 if( recover_typeII_Photon( ctrl, newmu, i, lepvec_i, pfCandidates ) ) {
532 if(ctrl.debug) cout << "FSR TYPEII :: oldpt: " << old_pt_i
533 << "\tnewpt: " << lepvec_i[i].vec.Pt()
534 << "\tindex: " << i
535 << endl;
536 //lepvec[i].fsrRecoveryAttempted=true;
537 }
538 } else {
539 const mithep::Electron *el = (mithep::Electron*)((*electronArr)[lepvec[i].index]);
540 mithep::Electron* newel = const_cast<mithep::Electron*>(el);
541 if( recover_typeI_Photon( ctrl, newel, i, lepvec_i, pfCandidates, electronArr, &Zvec, photonsToVeto ) ) {
542 if(ctrl.debug) cout << "FSR TYPEI :: oldpt: " << old_pt_i
543 << "\tnewpt: " << lepvec_i[i].vec.Pt()
544 << "\tindex: " << i
545 << endl;
546 //lepvec[i].fsrRecoveryAttempted=true;
547 }
548 }
549
550
551 cout << "j: " << j << endl;
552 mithep::ChargedParticle *cp;
553 if( abs(lepvec[j].type) == 13 ) {
554 const mithep::Muon *mu = (mithep::Muon*)((*muonArr)[lepvec[j].index]);
555 mithep::Muon * newmu = const_cast<mithep::Muon *>(mu);
556 if( recover_typeI_Photon( ctrl, newmu, j, lepvec_j, pfCandidates, electronArr, &Zvec, photonsToVeto ) ) {
557 if(ctrl.debug) cout << "FSR TYPEI :: oldpt: " << old_pt_j
558 << "\tnewpt: " << lepvec_j[j].vec.Pt()
559 << "\tindex: " << j
560 << endl;
561 //lepvec[j].fsrRecoveryAttempted=true;
562 }
563 if( recover_typeII_Photon( ctrl, newmu, j, lepvec_j, pfCandidates ) ) {
564 if(ctrl.debug) cout << "FSR TYPEII :: oldpt: " << old_pt_j
565 << "\tnewpt: " << lepvec_j[j].vec.Pt()
566 << "\tindex: " << j
567 << endl;
568 //lepvec[j].fsrRecoveryAttempted=true;
569 }
570 } else {
571 const mithep::Electron *el = (mithep::Electron*)((*electronArr)[lepvec[j].index]);
572 mithep::Electron* newel = const_cast<mithep::Electron*>(el);
573 if( recover_typeI_Photon( ctrl, newel, j, lepvec_j, pfCandidates, electronArr, &Zvec, photonsToVeto ) ) {
574 if(ctrl.debug) cout << "FSR TYPEI :: oldpt: " << old_pt_j
575 << "\tnewpt: " << lepvec_j[j].vec.Pt()
576 << "\tindex: " << j
577 << endl;
578 // lepvec[j].fsrRecoveryAttempted=true;
579 }
580 }
581
582 //
583 // now fix isolation
584 //
585
586 if( abs(lepvec[i].type) == 11 ) {
587 const mithep::Electron *el = (mithep::Electron*)((*electronArr)[lepvec_i[i].index]);
588 lepvec_i[i].status |=
589 (*ElectronIsoSelector)(ctrl,el,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto);
590
591 } else {
592 const mithep::Muon *mu = (mithep::Muon*)((*muonArr)[lepvec_i[i].index]);
593 lepvec_i[i].status |=
594 (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto);
595 }
596 updateSimpleLepton(lepvec_i[i]);
597
598 if( abs(lepvec[j].type) == 11 ) {
599 const mithep::Electron *el = (mithep::Electron*)((*electronArr)[lepvec_j[j].index]);
600 lepvec_j[j].status |=
601 (*ElectronIsoSelector)(ctrl,el,vtx,pfCandidates,puEnergyDensity,eraEle,photonsToVeto);
602 } else {
603 const mithep::Muon *mu = (mithep::Muon*)((*muonArr)[lepvec_j[j].index]);
604 lepvec_j[j].status |=
605 (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,photonsToVeto);
606 }
607 updateSimpleLepton(lepvec_j[j]);
608
609 float new_M = (lepvec_i[i].vec+lepvec_j[j].vec).M();
610 float new_pt_i = lepvec_i[i].vec.Pt();
611 float new_pt_j = lepvec_j[j].vec.Pt();
612 if( ctrl.debug ) {
613 cout << "\toldM: " << old_M << "\tnewM:" << new_M << endl;
614 cout << "\told_pt_i: " << old_pt_i << "\tnew_pt_i:" << new_pt_i << endl;
615 cout << "\told_pt_j: " << old_pt_j << "\tnew_pt_j:" << new_pt_j << endl;
616 }
617
618 } // doFSR
619
620 if( !(lepvec_i[i].status.loose()) || !(lepvec_j[j].status.loose()) ) continue;
621 ZCandidates.push_back(std::pair<int,int> (i,j) );
622 ZCandidatesLeptons.push_back(std::pair<SimpleLepton,SimpleLepton> (lepvec_i[i],lepvec_j[j]) );
623 if( ctrl.debug ) cout << "Z candidate ("<< i << "," << j << ")"
624 << "\tmass: " << (lepvec_i[i].vec+lepvec_j[j].vec).M() << endl;
625 }
626 }
627 if( ZCandidates.size() > 0 ) {
628 ret.status.selectionBits.flip(PASS_ZCANDIDATE);
629 if( ctrl.debug ) cout << "event has >0 Z candidates" << endl;
630 cutvec[PASS_ZCANDIDATE] +=1;
631 } else {
632 ret.status.setStatus(SelectionStatus::FAIL);
633 return ret;
634 }
635
636
637 //
638 // !!!!!!!!!!!!!! Z1 SELECTED HERE
639 //
640 int best_Z1_index;
641 float best_Z1_mass = 9999.;
642 TLorentzVector Z1vec;
643 for( int i=0; i<ZCandidates.size(); i++ ) {
644 // TLorentzVector tmpZ1vec = (lepvec[ZCandidates[i].first].vec) +
645 // (lepvec[ZCandidates[i].second].vec);
646 TLorentzVector tmpZ1vec = (ZCandidatesLeptons[i].first.vec) +
647 (ZCandidatesLeptons[i].second.vec);
648 if( fabs(tmpZ1vec.M()-Z_MASS) < fabs(best_Z1_mass-Z_MASS) ) {
649 best_Z1_index=i;
650 best_Z1_mass=tmpZ1vec.M();
651 Z1vec = tmpZ1vec;
652 }
653 }
654 // update the leptons, damn FSR ...
655 // ret.Z1leptons.push_back(lepvec[ZCandidates[best_Z1_index].first]);
656 // ret.Z1leptons.push_back(lepvec[ZCandidates[best_Z1_index].second]);
657 ret.Z1leptons.push_back(ZCandidatesLeptons[best_Z1_index].first);
658 ret.Z1leptons.push_back(ZCandidatesLeptons[best_Z1_index].second);
659
660
661 cout << "best mZ1: " << best_Z1_mass << " from (" <<
662 ZCandidates[best_Z1_index].first << "," << ZCandidates[best_Z1_index].second << ")" << endl;
663 int Z1type;
664 if( abs(ret.Z1leptons[0].type) == 11 ) Z1type=0;
665 else Z1type=1;
666 zcutvec[Z1type][PASS_ZCANDIDATE] +=1;
667
668
669 //******************************************************************************
670 // Step 6.3 : require Z1 with 40<m<120
671 //******************************************************************************
672 if( Z1vec.M() > 40. && Z1vec.M() < 120. ) {
673 ret.status.selectionBits.flip(PASS_GOODZ1);
674 cutvec[PASS_GOODZ1] +=1;
675 zcutvec[Z1type][PASS_GOODZ1] +=1;
676 } else {
677 ret.status.setStatus(SelectionStatus::FAIL);
678 return ret;
679 }
680
681 //******************************************************************************
682 // Step 6.3 : 4 good leptons
683 //******************************************************************************
684 if( lepvec.size() >= 4 ) {
685 if( ctrl.debug) cout << "pass4L: " << info->EvtNum() << endl;
686 ret.status.selectionBits.flip(PASS_4L);
687 cutvec[PASS_4L] +=1;
688 zcutvec[Z1type][PASS_4L] +=1;
689 } else {
690 ret.status.setStatus(SelectionStatus::FAIL);
691 return ret;
692 }
693 int nEl=0, nMu=0;
694 for( int i=0; i<4; i++ ) {
695 if(abs(lepvec[i].type) == 11 ) nEl++;
696 if(abs(lepvec[i].type) == 13 ) nMu++;
697 }
698 if( nEl >= 4 ) zzcutvec[0][PASS_4L] +=1;
699 else if( nMu >= 4 ) zzcutvec[1][PASS_4L] +=1;
700 else zzcutvec[2][PASS_4L] +=1;
701
702
703 //********************************************************
704 // Step 5: ZZ candidates
705 //********************************************************
706 int nZZCandidates=0;
707 std::vector<std::pair<int,int> > ZZCandidates;
708 int Z2type;
709 for(int z2index=0; z2index<ZCandidates.size(); ++z2index) {
710 int z1index = best_Z1_index;
711 if ( z2index == z1index ) continue;
712 if( ZCandidates[z1index].first == ZCandidates[z2index].first ) continue;
713 if( ZCandidates[z1index].first == ZCandidates[z2index].second ) continue;
714 if( ZCandidates[z1index].second == ZCandidates[z2index].first ) continue;
715 if( ZCandidates[z1index].second == ZCandidates[z2index].second ) continue;
716 //ghost removal
717 if(mithep::MathUtils::DeltaR(ZCandidatesLeptons[z1index].first.vec.Phi(),ZCandidatesLeptons[z1index].first.vec.Eta(),ZCandidatesLeptons[z2index].first.vec.Phi(),ZCandidatesLeptons[z2index].first.vec.Eta()) < 0.02) continue;
718 if(mithep::MathUtils::DeltaR(ZCandidatesLeptons[z1index].first.vec.Phi(),ZCandidatesLeptons[z1index].first.vec.Eta(),ZCandidatesLeptons[z2index].second.vec.Phi(),ZCandidatesLeptons[z2index].second.vec.Eta()) < 0.02) continue;
719 if(mithep::MathUtils::DeltaR(ZCandidatesLeptons[z1index].second.vec.Phi(),ZCandidatesLeptons[z1index].second.vec.Eta(),ZCandidatesLeptons[z2index].first.vec.Phi(),ZCandidatesLeptons[z2index].first.vec.Eta()) < 0.02) continue;
720 if(mithep::MathUtils::DeltaR(ZCandidatesLeptons[z1index].second.vec.Phi(),ZCandidatesLeptons[z1index].second.vec.Eta(),ZCandidatesLeptons[z2index].second.vec.Phi(),ZCandidatesLeptons[z2index].second.vec.Eta()) < 0.02) continue;
721
722 ZZCandidates.push_back(std::pair<int,int> (z1index,z2index));
723 // Z2type = abs(lepvec[ZCandidates[z2index].first].type);
724 Z2type = abs(ZCandidatesLeptons[z2index].first.type);
725 }
726 if( ZZCandidates.size() > 0 ) {
727 if( ctrl.debug) cout << "passZZ: " << info->EvtNum() << endl;
728 ret.status.selectionBits.flip(PASS_ZZCANDIDATE);
729 cutvec[PASS_ZZCANDIDATE] +=1;
730 if( ZZCandidates.size() > 1 ) Z2type=999;
731 if( ctrl.debug ) cout << "nZZcandidates: " << ZZCandidates.size() << endl;
732 // if( ctrl.debug ) {
733 cout << "evt: " << info->EvtNum() << "\tnZZcandidates: " << ZZCandidates.size()
734 << "\tZ1f: " << abs(ret.Z1leptons[0].type) << "\tZ2f: " << Z2type << endl;
735 cout << "-------------------------------------------------------" << endl;
736 for( int l=0; l<lepvec.size(); l++ ) lepvec[l].print();
737 cout << "-------------------------------------------------------" << endl;
738 // }
739 } else {
740 ret.status.setStatus(SelectionStatus::FAIL);
741 return ret;
742 }
743
744 //
745 // !!!!!!!!!!!!!! Z2 SELECTED HERE
746 //
747 int best_Z2_index;
748 float best_Z2_pt = -1.;
749 for( int i=0; i<ZZCandidates.size(); i++ ) {
750 int z2index = ZZCandidates[i].second;
751 // TLorentzVector Z2 = (lepvec[ZCandidates[z2index].first].vec) +
752 // (lepvec[ZCandidates[z2index].second].vec);
753 TLorentzVector Z2 = (ZCandidatesLeptons[z2index].first.vec) +
754 (ZCandidatesLeptons[z2index].second.vec);
755 if( Z2.Pt() > best_Z2_pt ) {
756 best_Z2_index=z2index;
757 best_Z2_pt=Z2.Pt();
758 }
759 }
760 // update the leptons, damn FSR ...
761 // ret.Z2leptons.push_back(lepvec[ZCandidates[best_Z2_index].first]);
762 // ret.Z2leptons.push_back(lepvec[ZCandidates[best_Z2_index].second]);
763 ret.Z2leptons.push_back(ZCandidatesLeptons[best_Z2_index].first);
764 ret.Z2leptons.push_back(ZCandidatesLeptons[best_Z2_index].second);
765
766
767 cout << "best mZ2: " << (ret.Z2leptons[0].vec+ret.Z2leptons[1].vec).M() << endl;
768 int ZZtype;
769 if( Z1type == 0 && abs(ret.Z2leptons[0].type) == 11 ) ZZtype=0;
770 else if ( Z1type == 1 && abs(ret.Z2leptons[0].type) == 13 ) ZZtype=1;
771 else ZZtype=2;
772 zzcutvec[ZZtype][PASS_ZZCANDIDATE] += 1;
773
774 if(ctrl.debug)cout << "ZZ :: evt: " << info->EvtNum()
775 << "\tmZ1: " << (ret.Z1leptons[0].vec+ret.Z1leptons[1].vec).M()
776 << "\tmZ2: " << (ret.Z2leptons[0].vec+ret.Z2leptons[1].vec).M()
777 << endl;
778
779
780
781
782 //******************************************************************************
783 // Step 6.4 : require Z2 with 4<m<120
784 //******************************************************************************
785 TLorentzVector Z2vec = (ZCandidatesLeptons[best_Z2_index].first.vec) +
786 (ZCandidatesLeptons[best_Z2_index].second.vec);
787 if( Z2vec.M() > 4 && Z2vec.M() < 120 ) {
788 ret.status.selectionBits.flip(PASS_GOODZ2);
789 cutvec[PASS_GOODZ2] +=1;
790 zzcutvec[ZZtype][PASS_GOODZ2] += 1;
791 } else {
792 ret.status.setStatus(SelectionStatus::FAIL);
793 return ret;
794 }
795
796
797 //******************************************************************************
798 // Step 6.1 : any two leptons 20/10
799 //******************************************************************************
800 vector<SimpleLepton> zzleptons;
801 zzleptons.push_back( ZCandidatesLeptons[best_Z1_index].first );
802 zzleptons.push_back( ZCandidatesLeptons[best_Z1_index].second );
803 zzleptons.push_back( ZCandidatesLeptons[best_Z2_index].first );
804 zzleptons.push_back( ZCandidatesLeptons[best_Z2_index].second );
805 nlep_above_10=0; nlep_above_20=0;
806 for( int i=0; i<zzleptons.size(); i++ ) {
807 if( zzleptons[i].vec.Pt() > 10 ) nlep_above_10++;
808 if( zzleptons[i].vec.Pt() > 20 ) nlep_above_20++;
809 }
810 if( nlep_above_10 > 1 && nlep_above_20 > 0 ) {
811 ret.status.selectionBits.flip(PASS_ZZ_20_10);
812 cutvec[PASS_ZZ_20_10] +=1;
813 zzcutvec[ZZtype][PASS_ZZ_20_10] += 1;
814 if( ctrl.debug ) cout << "passess 20/10 ..." << endl;
815 } else {
816 ret.status.setStatus(SelectionStatus::FAIL);
817 return ret;
818 }
819
820
821
822
823 //******************************************************************************
824 // Step 6.5 : resonance killing (4/4)
825 //******************************************************************************
826 bool resonance = false;
827 for( int i=0; i<zzleptons.size(); i++ ) {
828 for( int j=i+1; j<zzleptons.size(); j++ ) {
829 if( zzleptons[i].charge == zzleptons[j].charge ) continue; // 4/4
830 if( (zzleptons[i].vec+zzleptons[j].vec).M() < 4. ) {
831 resonance = true;
832 break;
833 }
834 }
835 }
836 if( !resonance ) {
837 ret.status.selectionBits.flip(PASS_RESONANCE);
838 cutvec[PASS_RESONANCE] +=1;
839 zzcutvec[ZZtype][PASS_RESONANCE] += 1;
840 if( ctrl.debug ) cout << "\tpasses resonance killing ... " << endl;
841 } else {
842 ret.status.setStatus(SelectionStatus::FAIL);
843 return ret;
844 }
845
846
847
848 //******************************************************************************
849 // Step 6.6 : m(4l) > 70 , m(4l) > 100
850 //******************************************************************************
851 TLorentzVector zzvec = (ZCandidatesLeptons[best_Z1_index].first.vec) +
852 (ZCandidatesLeptons[best_Z1_index].second.vec) +
853 (ZCandidatesLeptons[best_Z2_index].first.vec) +
854 (ZCandidatesLeptons[best_Z2_index].second.vec);
855
856 if( zzvec.M() > 70. ) {
857 if(ctrl.debug) cout << "passes mzz > 70, mzz: " << zzvec.M() << endl;
858 ret.status.selectionBits.flip(PASS_m4l_GT_70);
859 cutvec[PASS_m4l_GT_70] +=1;
860 zzcutvec[ZZtype][PASS_m4l_GT_70] += 1;
861 } else {
862 ret.status.setStatus(SelectionStatus::FAIL);
863 return ret;
864 }
865
866 if( zzvec.M() > 100. && (ret.Z2leptons[0].vec+ret.Z2leptons[1].vec).M() > 12) {
867 if(ctrl.debug) cout << "passes mzz > 100, mzz: " << zzvec.M() << "\t mZ2 > 12" << endl;
868 ret.status.selectionBits.flip(PASS_m4l_GT_100);
869 cutvec[PASS_m4l_GT_100] +=1;
870 zzcutvec[ZZtype][PASS_m4l_GT_100] += 1;
871 } else {
872 ret.status.setStatus(SelectionStatus::FAIL);
873 return ret;
874 }
875
876 //***************************************************************
877 // finish
878 //***************************************************************
879
880 TLorentzVector theZ1 = (ZCandidatesLeptons[best_Z1_index].first.vec) +
881 (ZCandidatesLeptons[best_Z1_index].second.vec);
882 TLorentzVector theZ2 = (ZCandidatesLeptons[best_Z2_index].first.vec) +
883 (ZCandidatesLeptons[best_Z2_index].second.vec);
884 TLorentzVector theZZ = theZ1 + theZ2;
885 int theZ1type = ZCandidatesLeptons[best_Z1_index].first.type;
886 int theZ2type = ZCandidatesLeptons[best_Z2_index].first.type;
887
888 if( ctrl.debug ) cout << "run: " << info->RunNum()
889 << "\tevt: " << info->EvtNum()
890 << "\tZ1channel: " << theZ1type
891 << "\tZ2channel: " << theZ2type
892 << "\tmZ1: " << theZ1.M()
893 << "\tmZ2: " << theZ2.M()
894 << "\tm4l: " << theZZ.M()
895 << "\tevtfail: " << hex << evtfail << dec
896 // << "\ttrigbits: " << hex << info->triggerBits << dec
897 // << "\ttree: " << inputFiles[q][f]
898 << endl;
899
900 if( !evtfail ) {
901 ret.status.setStatus(SelectionStatus::EVTPASS);
902 // already done ..
903 // ret.Z1leptons.push_back(lepvec[ZCandidates[best_Z1_index].first]);
904 // ret.Z1leptons.push_back(lepvec[ZCandidates[best_Z1_index].second]);
905 // ret.Z2leptons.push_back(lepvec[ZCandidates[best_Z2_index].first]);
906 // ret.Z2leptons.push_back(lepvec[ZCandidates[best_Z2_index].second]);
907 }
908
909 return ret;
910 }
911 //--------------------------------------------------------------------------------------------------
912 void fillVetoArrays( ControlFlags & ctrl,
913 const mithep::Array<mithep::Muon> *muonArr,
914 vector< const mithep::Muon*> & muonsToVeto,
915 const mithep::Array<mithep::Electron> *electronArr,
916 vector< const mithep::Electron*> & electronsToVeto,
917 const mithep::Vertex * vtx )
918 //--------------------------------------------------------------------------------------------------
919 {
920
921 if( ctrl.debug ) cout << "looping for isolation ..." << endl;
922 /*
923 for(int i=0; i<muonArr->GetEntries(); i++)
924 {
925 const mithep::Muon *mu = (const mithep::Muon*)((*muonArr)[i]);
926 SelectionStatus musel;
927 // musel |= muonCutBasedVeto(ctrl,mu,vtx);
928 musel |= muonDummyVeto(ctrl,mu,vtx);
929 if( !(musel.getStatus() & SelectionStatus::PRESELECTION) ) continue;
930 if(ctrl.debug) cout << "pushing mu for isol veto ... " << endl;
931 muonsToVeto.push_back( mu );
932 }
933 */
934 for(int i=0; i<electronArr->GetEntries(); i++)
935 {
936 const mithep::Electron *ele = (const mithep::Electron*)((*electronArr)[i]);
937 SelectionStatus esel;
938 // esel |= electronCutBasedVeto(ctrl,ele,vtx);
939 esel |= electronDummyVeto(ctrl,ele,vtx);
940 if( !(esel.getStatus() & SelectionStatus::PRESELECTION) ) continue;
941 if(ctrl.debug) cout << "pushing ele for isol veto ... " << endl;
942 electronsToVeto.push_back( ele );
943 }
944 if( ctrl.debug ) cout << "done selecting for isolation veto ..." << endl << endl;;
945 }
946
947
948
949
950 //----------------------------------------------------------------------------
951 void updateSimpleLepton(SimpleLepton &tmplep)
952 //----------------------------------------------------------------------------
953 {
954 tmplep.isoPF04 = tmplep.status.isoPF04;
955 tmplep.chisoPF04 = tmplep.status.chisoPF04;
956 tmplep.gaisoPF04 = tmplep.status.gaisoPF04;
957 tmplep.neisoPF04 = tmplep.status.neisoPF04;
958 tmplep.isTight = tmplep.status.tight();
959 tmplep.isLoose = tmplep.status.loose();
960 // tmplep.fsrRecoveryAttempted = true;
961 }