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