ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/Selection.cc
Revision: 1.24
Committed: Wed May 9 14:57:20 2012 UTC (13 years ago) by khahn
Content type: text/plain
Branch: MAIN
CVS Tags: synced_FSR_2, synced_FSR, synched2
Changes since 1.23: +5 -3 lines
Log Message:
*** empty log message ***

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