ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/SelectionZ4L.cc
Revision: 1.2
Committed: Fri Feb 17 15:09:42 2012 UTC (13 years, 3 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.1: +1 -1 lines
Log Message:
*** empty log message ***

File Contents

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