ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/Selection.cc
Revision: 1.21
Committed: Mon Apr 30 21:42:16 2012 UTC (13 years ago) by khahn
Content type: text/plain
Branch: MAIN
CVS Tags: synched
Changes since 1.20: +22 -9 lines
Log Message:
synched ...

File Contents

# User Rev Content
1 khahn 1.19 #include "SelectionStatus.h"
2     #include "EventData.h"
3     #include "SimpleLepton.h"
4     #include "EfficiencyWeightsInterface.h"
5 khahn 1.7
6 khahn 1.20 #include "ElectronSelection.h"
7     #include "MuonSelection.h"
8 khahn 1.19 #include "IsolationSelection.h"
9 khahn 1.20 //#include "PassHLT.h"
10 khahn 1.19 #include "Selection.h"
11 khahn 1.15
12 khahn 1.19 #include "ExternData.h"
13     #include "SelectionDefs.h"
14 khahn 1.1
15 khahn 1.13
16 khahn 1.19 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
17     // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
18     // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
19     EventData apply_HZZ4L_selection(ControlFlags &ctrl, // input control
20 khahn 1.20 const mithep::EventHeader *info, // input event info
21     const mithep::Vertex & vtx,
22     const mithep::Array<mithep::PFCandidate> *pfCandidates,
23     const mithep::Array<mithep::PileupEnergyDensity> *puEnergyDensity,
24     const mithep::Array<mithep::Electron> *electronArr, // input electrons
25     SelectionStatus (*ElectronPreSelector)( ControlFlags &,
26     const mithep::Electron*,
27     const mithep::Vertex &),
28     SelectionStatus (*ElectronIDSelector)( ControlFlags &,
29     const mithep::Electron*,
30     const mithep::Vertex &),
31     SelectionStatus (*ElectronIsoSelector)( ControlFlags &,
32     const mithep::Electron*,
33     const mithep::Vertex &,
34     const mithep::Array<mithep::PFCandidate> *,
35     const mithep::Array<mithep::PileupEnergyDensity> *,
36     mithep::ElectronTools::EElectronEffectiveAreaTarget,
37     vector<const mithep::Muon*>,
38     vector<const mithep::Electron*> ),
39     const mithep::Array<mithep::Muon> *muonArr, // input muons
40     SelectionStatus (*MuonPreSelector)( ControlFlags &,
41     const mithep::Muon*,
42     const mithep::Vertex &,
43     const mithep::Array<mithep::PFCandidate> *),
44     SelectionStatus (*MuonIDSelector)( ControlFlags &,
45     const mithep::Muon*,
46     const mithep::Vertex &),
47     SelectionStatus (*MuonIsoSelector)( ControlFlags &,
48     const mithep::Muon*,
49     const mithep::Vertex &,
50     const mithep::Array<mithep::PFCandidate> *,
51     const mithep::Array<mithep::PileupEnergyDensity> *,
52     mithep::MuonTools::EMuonEffectiveAreaTarget,
53     vector<const mithep::Muon*>,
54     vector<const mithep::Electron*> )
55     )
56 khahn 1.19 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
57     {
58 khahn 1.12
59 khahn 1.19 EventData ret;
60 khahn 1.1 unsigned evtfail = 0x0;
61 khahn 1.19 TRandom3 r;
62    
63 khahn 1.5 if( ctrl.debug ) {
64 khahn 1.21 cout << "-----------------------------------------------------------------" << endl;
65     cout << "-----------------------------------------------------------------" << endl;
66 khahn 1.20 cout << "Run: " << info->RunNum()
67     << "\tEvt: " << info->EvtNum()
68     << "\tLumi: " << info->LumiSec()
69 khahn 1.5 << endl;
70 khahn 1.21 cout << "-----------------------------------------------------------------" << endl;
71 khahn 1.5 }
72    
73 khahn 1.1 if( !ctrl.mc ) {
74     // not accounting for overlap atm
75 khahn 1.20 RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
76 khahn 1.1 if( !(rlrm.HasRunLumi(rl)) ) {
77 khahn 1.5 if( ctrl.debug ) cout << "\tfails JSON" << endl;
78 khahn 1.19 ret.status.setStatus(0);
79     return ret;
80 khahn 1.1 }
81     }
82 khahn 1.20
83 khahn 1.21 mithep::MuonTools::EMuonEffectiveAreaTarget eraMu = mithep::MuonTools::kMuEAFall11MC;
84     mithep::ElectronTools::EElectronEffectiveAreaTarget eraEle = mithep::ElectronTools::kEleEAFall11MC;
85     if( !ctrl.mc ) {
86     eraMu = mithep::MuonTools::kMuEAData2011;
87     eraEle = mithep::ElectronTools::kEleEAData2011;
88 khahn 1.20 }
89 khahn 1.1
90    
91     //********************************************************
92     // Trigger
93     //********************************************************
94 khahn 1.20 //
95     // still have to port this part to bambu
96     //
97     /*
98 khahn 1.1 if( !ctrl.mc ) {
99 khahn 1.19 //if( !(passHLT(info->triggerBits, info->runNum, channel) ) ) {
100     if( !(passHLT(info->triggerBits, info->runNum, 999) ) ) {
101     if( ctrl.debug ) cout << "\tfails trigger" << endl;
102 khahn 1.1 evtfail |= (1<<EVTFAIL_TRIGGER);
103 khahn 1.19 ret.status.setStatus(0);
104     return ret;
105     }
106     }
107 khahn 1.21 */
108    
109 khahn 1.1 if( ctrl.debug ) {
110     cout << "presel nlep: " << muonArr->GetEntries() + electronArr->GetEntries()
111     << "\tnmuon: " << muonArr->GetEntries()
112     << "\tnelectron: " << electronArr->GetEntries()
113     << endl;
114     }
115 khahn 1.21
116 khahn 1.1
117     //********************************************************
118     // Lepton Selection
119     //********************************************************
120     vector<SimpleLepton> lepvec;
121 khahn 1.20
122     // do something hacky for vetos for now
123     vector<const mithep::Muon*> muonsToVeto;
124     vector<const mithep::Electron*> electronsToVeto;
125 khahn 1.21 if( ctrl.debug ) cout << "looping for isolation ..." << endl;
126 khahn 1.20 for(Int_t i=0; i<muonArr->GetEntries(); i++)
127     {
128     const mithep::Muon *mu = (mithep::Muon*)((*muonArr)[i]);
129     SelectionStatus musel;
130     musel |= (*MuonPreSelector)(ctrl,mu,vtx,pfCandidates);
131     if( !(musel.getStatus() & SelectionStatus::PRESELECTION) ) continue;
132     musel |= (*MuonIDSelector)(ctrl,mu,vtx );
133 khahn 1.21 if(ctrl.debug) cout << "status : " << musel.getStatus() << endl;
134     if( musel.tightIDAndPre() ) {
135     if(ctrl.debug) cout << "pushing mu for isol veto ... " << endl;
136 khahn 1.20 muonsToVeto.push_back( mu );
137 khahn 1.21 }
138 khahn 1.20 }
139     for(Int_t i=0; i<electronArr->GetEntries(); i++)
140     {
141     const mithep::Electron *ele = (mithep::Electron*)((*electronArr)[i]);
142     SelectionStatus esel;
143     esel |= (*ElectronPreSelector)(ctrl,ele,vtx);
144     if( !(esel.getStatus() & SelectionStatus::PRESELECTION) ) continue;
145     esel |= (*ElectronIDSelector)(ctrl,ele,vtx );
146 khahn 1.21 if(ctrl.debug) cout << "status : " << esel.getStatus() << endl;
147     if( esel.tightIDAndPre() ) {
148     if(ctrl.debug) cout << "pushing ele for isol veto ... " << endl;
149 khahn 1.20 electronsToVeto.push_back( ele );
150 khahn 1.21 }
151 khahn 1.20 }
152 khahn 1.21 if( ctrl.debug ) cout << "done looping for isolation ..." << endl << endl;;
153 khahn 1.20
154 khahn 1.1 //
155 dkralph 1.3 if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
156 khahn 1.1 //----------------------------------------------------
157 khahn 1.19 for(Int_t i=0; i<muonArr->GetEntries(); i++)
158     {
159 khahn 1.20 const mithep::Muon *mu = (mithep::Muon*)((*muonArr)[i]);
160 khahn 1.19
161     SelectionStatus musel;
162     if(ctrl.debug) cout << "musel.status before anything: " << musel.getStatus() << endl;
163 khahn 1.20 musel |= (*MuonPreSelector)(ctrl,mu,vtx,pfCandidates);
164 khahn 1.19 if(ctrl.debug) cout << "musel.status after presel: " << musel.getStatus() << endl;
165 khahn 1.20 if( !(musel.getStatus() & SelectionStatus::PRESELECTION) ) continue;
166     musel |= (*MuonIDSelector)(ctrl,mu,vtx );
167 khahn 1.19 if(ctrl.debug) cout << "musel.status after ID: " << musel.getStatus() << endl;
168 khahn 1.20 musel |= (*MuonIsoSelector)(ctrl,mu,vtx,pfCandidates,puEnergyDensity,eraMu,muonsToVeto,electronsToVeto);
169 khahn 1.19 if(ctrl.debug) cout << "musel.status after iso: " << musel.getStatus() << endl;
170    
171     if( ctrl.debug ) {
172 khahn 1.20 cout << "muon:: pt: " << mu->Pt()
173     << "\teta: " << mu->Eta()
174 khahn 1.19 << "\tstatus: " << hex << musel.getStatus() << dec
175     << endl;
176     }
177    
178     if ( musel.pass() ) {
179 khahn 1.17
180 khahn 1.19 SimpleLepton tmplep;
181 khahn 1.20 float pt = mu->Pt();
182 khahn 1.19 tmplep.vecorig->SetPtEtaPhiM(pt,
183 khahn 1.20 mu->Eta(),
184     mu->Phi(),
185     MUON_MASS);
186 khahn 1.19
187     if( ctrl.do_escale_up ) {
188     pt=scale_smear_muon_Up(pt, 1, r);
189     }
190     if( ctrl.do_escale_down ) {
191     pt=scale_smear_muon_Down(pt, 1, r);
192     }
193    
194     tmplep.vec->SetPtEtaPhiM(pt,
195 khahn 1.20 mu->Eta(),
196     mu->Phi(),
197     MUON_MASS);
198 khahn 1.19
199     tmplep.type = 13;
200     tmplep.index = i;
201 khahn 1.20 tmplep.charge = mu->Charge();
202     tmplep.isoTrk = mu->IsoR03SumPt();
203     tmplep.isoEcal = mu->IsoR03EmEt();
204     tmplep.isoHcal = mu->IsoR03HadEt();
205     tmplep.isoPF03 = computePFMuonIso(mu,vtx,pfCandidates,0.3);
206     tmplep.isoPF04 = computePFMuonIso(mu,vtx,pfCandidates,0.4);
207     tmplep.ip3dSig = mu->Ip3dPVSignificance();
208 khahn 1.19 tmplep.is4l = false;
209 khahn 1.20 tmplep.isEB = (fabs(mu->Eta()) < 1.479 ? 1 : 0 );
210 khahn 1.19 tmplep.isTight = musel.tight();
211     tmplep.isLoose = musel.loose();
212     lepvec.push_back(tmplep);
213     if( ctrl.debug ) { cout << "muon passes ... " << endl;}
214     }
215     // }
216     }
217    
218 khahn 1.1
219 khahn 1.19
220     //
221     if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
222     // --------------------------------------------------------------------------------
223     for(Int_t i=0; i<electronArr->GetEntries(); i++)
224     {
225 khahn 1.20 const mithep::Electron *ele = (mithep::Electron*)((*electronArr)[i]);
226 khahn 1.19
227     Bool_t isMuonOverlap = kFALSE;
228     for (int k=0; k<lepvec.size(); ++k) {
229     TVector3 tmplep;
230 khahn 1.20 tmplep.SetPtEtaPhi(ele->Pt(), ele->Eta(), ele->Phi());
231 khahn 1.19 if ( lepvec[k].isLoose && lepvec[k].type == 13 && lepvec[k].vec->Vect().DrEtaPhi(tmplep) < 0.1 ) {
232     if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
233     isMuonOverlap = kTRUE;
234     break;
235     }
236     }
237    
238     SelectionStatus elesel;
239     if( ctrl.debug ) cout << "--> status before anything: " << hex << elesel.getStatus() << dec << endl;
240 khahn 1.20 elesel |= (*ElectronPreSelector)(ctrl,ele,vtx);
241 khahn 1.19 if( ctrl.debug ) cout << "--> status after presel: " << hex << elesel.getStatus() << dec << endl;
242 khahn 1.20 elesel |= (*ElectronIDSelector)(ctrl,ele,vtx);
243 khahn 1.19 if( ctrl.debug ) cout << "--> status after ID: " << hex << elesel.getStatus() << dec << endl;
244 khahn 1.20 elesel |= (*ElectronIsoSelector)(ctrl,ele,vtx,pfCandidates,puEnergyDensity,eraEle,muonsToVeto,electronsToVeto);
245 khahn 1.19 if( ctrl.debug ) cout << "--> status after iso: " << hex << elesel.getStatus() << dec << endl;
246    
247     if( ctrl.debug ){
248 khahn 1.20 cout << "\tscEt: " << ele->SCluster()->Et()
249     << "\tscEta: " << ele->SCluster()->Eta()
250 khahn 1.19 << "\tstatus: " << hex << elesel.getStatus() << dec
251     << endl;
252     }
253 dkralph 1.3
254 khahn 1.19 if ( elesel.pass() && !isMuonOverlap )
255     {
256     SimpleLepton tmplep;
257    
258 khahn 1.20 float pt = ele->Pt();
259 khahn 1.19 tmplep.vecorig->SetPtEtaPhiM( pt,
260 khahn 1.20 ele->Eta(),
261     ele->Phi(),
262 khahn 1.19 ELECTRON_MASS );
263    
264     if( ctrl.do_escale ) {
265 khahn 1.20 pt=scale_smear_electron(pt, ele->IsEB(), r);
266 khahn 1.19 }
267     if( ctrl.do_escale_up ) {
268 khahn 1.20 pt=scale_smear_electron_Up(pt, ele->IsEB(), r);
269 khahn 1.19 }
270     if( ctrl.do_escale_down ) {
271 khahn 1.20 pt=scale_smear_electron_Down(pt, ele->IsEB(), r);
272 khahn 1.19 }
273 khahn 1.1
274 khahn 1.19
275     tmplep.vec->SetPtEtaPhiM( pt,
276 khahn 1.20 ele->Eta(),
277     ele->Phi(),
278 khahn 1.19 ELECTRON_MASS );
279    
280     tmplep.type = 11;
281     tmplep.index = i;
282 khahn 1.20 tmplep.charge = ele->Charge();
283     tmplep.isoTrk = ele->TrackIsolationDr03();
284     tmplep.isoEcal = ele->EcalRecHitIsoDr03();
285     tmplep.isoHcal = ele->HcalTowerSumEtDr03();
286     tmplep.isoPF03 = computePFEleIso(ele,vtx,pfCandidates,0.3);
287     tmplep.isoPF04 = computePFEleIso(ele,vtx,pfCandidates,0.4);
288     tmplep.ip3dSig = ele->Ip3dPVSignificance();
289 khahn 1.19 tmplep.is4l = false;
290 khahn 1.20 tmplep.isEB = ele->IsEB();
291     tmplep.scID = ele->SCluster()->GetUniqueID();
292 khahn 1.19 tmplep.isTight = elesel.tight();
293     tmplep.isLoose = elesel.loose();
294     lepvec.push_back(tmplep);
295     if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
296     }
297 khahn 1.1 }
298 khahn 1.19
299    
300     //********************************************************
301     // Dump Stuff
302     //********************************************************
303     sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
304     int nmu=0, nele=0;
305     for( int i=0; i<lepvec.size(); i++ ) {
306     if(ctrl.debug) cout << "lepvec :: index: " << i
307     << "\tpt: " << lepvec[i].vec->Pt()
308     << "\ttype: " << lepvec[i].type
309     << endl;
310     if( abs(lepvec[i].type) == 11 ) nele++;
311     else nmu++;
312 khahn 1.1 }
313 khahn 1.19 if( ctrl.debug ) {
314     cout << "postsel nlep: " << lepvec.size()
315     << "\tnmuon: " << nmu
316     << "\tnelectron: " << nele
317 khahn 1.1 << endl;
318     }
319 khahn 1.19
320    
321     //******************************************************************************
322     // Z1 Selection
323     //******************************************************************************
324     int Z1LeptonPlusIndex = -1;
325     int Z1LeptonMinusIndex = -1;
326     double BestZ1Mass = -999;
327     if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
328     for(int i = 0; i < lepvec.size(); ++i) {
329     if( !(lepvec[i].isLoose) ) continue;
330     for(int j = i+1; j < lepvec.size(); ++j) {
331     if( !(lepvec[j].isLoose) ) continue;
332     if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
333     if (!(lepvec[i].vec->Pt() > 20.0 || lepvec[j].vec->Pt() > 20.0)) continue;
334     if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
335     if (!(lepvec[i].vec->Pt() > 10.0 && lepvec[j].vec->Pt() > 10.0)) continue;
336     if( ctrl.debug ) { cout << "\tthe other is > 10 GeV" << endl; }
337     if (lepvec[i].charge == lepvec[j].charge) continue;
338     if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
339     if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;
340     if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
341    
342     //Make Z1 hypothesis
343     TLorentzVector *leptonPlus, *leptonMinus;
344     if ( lepvec[i].charge > 0 ) {
345     leptonPlus = lepvec[i].vec;
346     leptonMinus = lepvec[j].vec;
347     } else {
348     leptonPlus = lepvec[j].vec;
349     leptonMinus = lepvec[i].vec;
350     }
351    
352     float tmpZ1Mass = (*leptonPlus + *leptonMinus).M();
353     if( ctrl.debug ) cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
354 khahn 1.21 if( tmpZ1Mass > 50 ) {
355 khahn 1.19 if (fabs(tmpZ1Mass - Z_MASS) < fabs(BestZ1Mass - Z_MASS)) {
356     BestZ1Mass = tmpZ1Mass;
357     if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass
358     << "\tdM: " << fabs(BestZ1Mass - Z_MASS)
359     << endl;
360     if (lepvec[i].charge > 0) {
361     Z1LeptonPlusIndex = i;
362     Z1LeptonMinusIndex = j;
363     } else {
364     Z1LeptonPlusIndex = j;
365     Z1LeptonMinusIndex = i;
366     }
367     }
368 khahn 1.13 }
369     }
370     }
371 khahn 1.19 // stop if no Z1 candidate is found
372     if( BestZ1Mass < 0 ) {
373     evtfail |= (1<<EVTFAIL_Z1);
374     //ret.status = evtfail;
375     ret.status.setStatus(0);
376     return ret;
377     }
378 khahn 1.20 if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->RunNum() << "\tevt: " << info->EvtNum() << endl;
379 khahn 1.19 if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
380     TLorentzVector Z1LeptonPlus = *(lepvec[Z1LeptonPlusIndex].vec);
381     TLorentzVector Z1LeptonMinus = *(lepvec[Z1LeptonMinusIndex].vec);
382     TLorentzVector Z1Candidate = Z1LeptonPlus + Z1LeptonMinus;
383    
384    
385     //******************************************************************************
386     // Z1 + l
387     //******************************************************************************
388     if( lepvec.size() < 3 ) {
389     evtfail |= (1<<EVTFAIL_Z1_PLUSL);
390     //ret.status = evtfail;
391     ret.status.setStatus(0);
392     return ret;
393     }
394    
395     //******************************************************************************
396     // 4l/Z2 Selection
397     //******************************************************************************
398     Int_t Z2LeptonPlusIndex = -1;
399     Int_t Z2LeptonMinusIndex = -1;
400     Double_t BestZ2Mass = -1;
401     if( ctrl.debug ) cout << "looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
402     for(int i = 0; i < lepvec.size(); ++i) {
403    
404     if( ctrl.debug) cout << "i: " << i
405     << "\tpt: " << lepvec[i].vec->Pt()
406     << "\ttype: " << lepvec[i].type
407     << endl;
408 khahn 1.1
409 khahn 1.19 if( ctrl.eleSeleScheme == "mediumloose" &&
410     !(lepvec[i].isTight) ) {
411     if( ctrl.debug) cout << "it's not tight, skipping ... " << endl;
412     continue;
413 khahn 1.1 }
414    
415 khahn 1.19 for(int j = i+1; j < lepvec.size(); ++j) {
416     if( ctrl.debug) cout << "\t\tj: " << j
417     << "\tpt: " << lepvec[j].vec->Pt()
418     << "\ttype: " << lepvec[j].type
419     << endl;
420    
421     if( ctrl.eleSeleScheme == "mediumloose" &&
422     !(lepvec[j].isTight) ) {
423     if( ctrl.debug) cout << "it's not tight, skipping ... " << endl;
424     continue;
425     }
426    
427    
428     if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) {
429     if( ctrl.debug) cout << "\ti matches a Z1 index, skipping ..." << endl;
430     continue; //skip Z1 leptons
431     }
432     if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) {
433     if( ctrl.debug) cout << "\tj matches a Z1 index, skipping ..." << endl;
434     continue; //skip Z1 leptons
435     }
436     if (lepvec[i].charge == lepvec[j].charge) {
437     if( ctrl.debug) cout << "\ti and j are same sign, skipping ..." << endl;
438     continue; //require opp sign
439     }
440     if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
441     if( ctrl.debug) cout << "\ti and j are not same flavor, skipping ..." << endl;
442     continue; //require same flavor
443     }
444    
445    
446     //Make Z2 hypothesis
447     TLorentzVector *leptonPlus, *leptonMinus;
448    
449     if (lepvec[i].charge > 0 ) {
450     leptonPlus = lepvec[i].vec;
451     leptonMinus = lepvec[j].vec;
452     } else {
453     leptonPlus = lepvec[j].vec;
454     leptonMinus = lepvec[i].vec;
455     }
456    
457     TLorentzVector dilepton = *leptonPlus + *leptonMinus;
458     TLorentzVector fourLepton = Z1Candidate + dilepton;
459    
460     if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
461     if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;
462    
463     if (!(dilepton.M() > 12.0)) continue;
464     if (!(fourLepton.M() > 100.0)) continue;
465    
466     //for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
467     if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
468     TLorentzVector pair1 = Z1LeptonPlus + *leptonMinus;
469     TLorentzVector pair2 = Z1LeptonMinus + *leptonPlus;
470     if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
471     if (!(pair1.M() > 12 || pair2.M() > 12)) continue;
472     }
473    
474    
475     //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
476     if (Z2LeptonPlusIndex < 0) {
477 khahn 1.1 if (lepvec[i].charge > 0) {
478 khahn 1.19 Z2LeptonPlusIndex = i;
479     Z2LeptonMinusIndex = j;
480 khahn 1.1 } else {
481 khahn 1.19 Z2LeptonPlusIndex = j;
482     Z2LeptonMinusIndex = i;
483     }
484     } else {
485     Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec->Pt();
486     Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec->Pt();
487     if (lepvec[Z2LeptonMinusIndex].vec->Pt() > BestPairPtMax) {
488     BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec->Pt();
489     BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec->Pt();
490     }
491    
492     Double_t CurrentPairPtMax = lepvec[i].vec->Pt();
493     Double_t CurrentPairPtMin = lepvec[j].vec->Pt();
494     if (lepvec[j].vec->Pt() > CurrentPairPtMax) {
495     CurrentPairPtMax = lepvec[j].vec->Pt();
496     CurrentPairPtMin = lepvec[i].vec->Pt();
497 khahn 1.1 }
498    
499 khahn 1.19 if (CurrentPairPtMax > BestPairPtMax) {
500     if (lepvec[i].charge > 0) {
501     Z2LeptonPlusIndex = i;
502     Z2LeptonMinusIndex = j;
503 khahn 1.1 } else {
504 khahn 1.19 Z2LeptonPlusIndex = j;
505     Z2LeptonMinusIndex = i;
506 khahn 1.1 }
507 khahn 1.19 } else if (CurrentPairPtMax == BestPairPtMax) {
508     if (CurrentPairPtMin > BestPairPtMin) {
509     if (lepvec[i].charge > 0) {
510     Z2LeptonPlusIndex = i;
511     Z2LeptonMinusIndex = j;
512 khahn 1.1 } else {
513 khahn 1.19 Z2LeptonPlusIndex = j;
514     Z2LeptonMinusIndex = i;
515     }
516 khahn 1.1 }
517 khahn 1.9 }
518 khahn 1.19 }
519     }
520     }
521    
522     // stop if no Z2 candidate is found
523     if (Z2LeptonPlusIndex == -1) {
524     evtfail |= ( 1<<EVTFAIL_4L );
525     // ret.status = evtfail;
526     ret.status.setStatus(0);
527     return ret;
528     }
529     if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
530     if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
531     << "\tminusindex: " << Z2LeptonMinusIndex << endl;
532     TLorentzVector Z2LeptonPlus = *(lepvec[Z2LeptonPlusIndex].vec);
533     TLorentzVector Z2LeptonMinus = *(lepvec[Z2LeptonMinusIndex].vec);
534     TLorentzVector Z2Candidate = Z2LeptonPlus+Z2LeptonMinus;
535     TLorentzVector ZZSystem = Z1Candidate + Z2Candidate;
536     lepvec[Z1LeptonPlusIndex].is4l = true;
537     lepvec[Z1LeptonMinusIndex].is4l = true;
538     lepvec[Z2LeptonPlusIndex].is4l = true;
539     lepvec[Z2LeptonMinusIndex].is4l = true;
540    
541 khahn 1.8
542    
543 khahn 1.19
544     //***************************************************************
545     // remaining kinematic cuts
546     //***************************************************************
547     double Z2massCut=0;
548     if ( ctrl.kinematics == "loose" ) Z2massCut = 12;
549     else if ( ctrl.kinematics == "tight" ) Z2massCut = 20;
550     else { cout << "error! kinematic tightness not defined!" << endl; assert(0); }
551    
552     if ( Z1Candidate.M() > 120 ||
553     Z2Candidate.M() < Z2massCut ||
554     Z2Candidate.M() > 120 ||
555     !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec->Pt() > 20.0) ||
556     !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec->Pt() > 10.0)
557     ) {
558     evtfail |= (1<<EVTFAIL_KINEMATICS );
559     // ret.status = evtfail;
560     ret.status.setStatus(0);
561     return ret;
562     }
563    
564     unsigned channel;
565     if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
566     if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
567     if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) ||
568     (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
569    
570    
571    
572    
573 khahn 1.20 if( ctrl.debug ) cout << "run: " << info->RunNum()
574     << "\tevt: " << info->EvtNum()
575 khahn 1.19 << "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
576     << "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
577     << "\tmZ1: " << Z1Candidate.M()
578     << "\tmZ2: " << Z2Candidate.M()
579     << "\tm4l: " << ZZSystem.M()
580     << "\tevtfail: " << hex << evtfail << dec
581 khahn 1.20 // << "\ttrigbits: " << hex << info->triggerBits << dec
582 khahn 1.19 // << "\ttree: " << inputFiles[q][f]
583     << endl;
584    
585 khahn 1.14
586    
587 khahn 1.19 //***************************************************************
588     // finish
589     //***************************************************************
590    
591     if( !evtfail ) {
592     ret.status.setStatus(SelectionStatus::EVTPASS);
593     ret.Z1leptons.push_back(lepvec[Z1LeptonMinusIndex]);
594     ret.Z1leptons.push_back(lepvec[Z1LeptonPlusIndex]);
595     ret.Z2leptons.push_back(lepvec[Z2LeptonMinusIndex]);
596     ret.Z2leptons.push_back(lepvec[Z2LeptonPlusIndex]);
597     }
598 khahn 1.1
599 khahn 1.19 return ret;
600 khahn 1.1 }
601    
602