ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/ReferenceSelection.cc
Revision: 1.14
Committed: Sun Jul 8 17:26:58 2012 UTC (12 years, 10 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.13: +6 -0 lines
Log Message:
added ghost removal to the reference selection

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