ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/ReferenceSelection.cc
Revision: 1.15
Committed: Wed Jul 25 13:29:14 2012 UTC (12 years, 9 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.14: +2 -0 lines
Log Message:
added electron energy regression stuff

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