ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/ReferenceSelection.cc
Revision: 1.4
Committed: Mon May 28 16:41:19 2012 UTC (12 years, 11 months ago) by khahn
Content type: text/plain
Branch: MAIN
CVS Tags: synched2
Changes since 1.3: +38 -32 lines
Log Message:
adding FSR recovery

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