ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/Selection.cc
(Generate patch)

Comparing UserCode/MitHzz4l/Selection/src/Selection.cc (file contents):
Revision 1.11 by khahn, Mon Oct 17 16:23:43 2011 UTC vs.
Revision 1.19 by khahn, Mon Feb 13 09:47:12 2012 UTC

# Line 1 | Line 1
1 < #include "Selection.h"
2 < #include "PassHLT.h"
3 <
4 < #include "SiMVAElectronSelection.h"
1 > #include "SelectionStatus.h"
2 > #include "EventData.h"
3 > #include "SimpleLepton.h"
4 > #include "EfficiencyWeightsInterface.h"
5  
6 #include "HZZCiCElectronSelection.h"
7 #include "HZZLikelihoodElectronSelection.h"
6   #include "HZZBDTElectronSelection.h"
7 < #include "RunLumiRangeMap.h"
8 <
9 < RunLumiRangeMap rlrm;
7 > #include "IsolationSelection.h"
8 > #include "PassHLT.h"
9 > #include "Selection.h"
10  
11 < void initRunLumiRangeMap() {
12 <  rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
15 <  //  rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
16 <  rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
17 <  rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));  
18 <  rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));  
19 < };
20 <
21 < unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control
22 <                               mithep::TEventInfo *info,     // input event inof
23 <                               TClonesArray *electronArr,    // input electrons
24 <                               TClonesArray *muonArr,        // input muons
25 <                               double eventweight,           // weight
26 <                               TTree * passtuple ) {
27 <
28 <  fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, passtuple, NULL );
29 <
30 < };
31 <
32 < unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control
33 <                               mithep::TEventInfo *info,     // input event inof
34 <                               TClonesArray *electronArr,    // input electrons
35 <                               TClonesArray *muonArr,        // input muons
36 <                               double eventweight,           // weight
37 <                               LabVectors  *l ) {
38 <
39 <  fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, NULL, l );
40 <
41 < };
42 <
43 <
44 < unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control
45 <                               mithep::TEventInfo *info,     // input event inof
46 <                               TClonesArray *electronArr,    // input electrons
47 <                               TClonesArray *muonArr,        // input muons
48 <                               double eventweight,           // weight
49 <                               TTree * passtuple,
50 <                               LabVectors * l) {       // output ntuple
11 > #include "ExternData.h"
12 > #include "SelectionDefs.h"
13  
52  unsigned evtfail = 0x0;
14  
15 + // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16 + // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
17 + // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
18 + EventData apply_HZZ4L_selection(ControlFlags &ctrl,           // input control
19 +                                mithep::TEventInfo *info,     // input event info
20 +                                TClonesArray *electronArr,    // input electrons
21 +                                SelectionStatus (*ElectronPreSelector)( ControlFlags &, const mithep::TElectron*),
22 +                                SelectionStatus (*ElectronIDSelector)( ControlFlags &, const mithep::TElectron*),
23 +                                SelectionStatus (*ElectronIsoSelector)( ControlFlags &, const mithep::TElectron*),
24 +                                TClonesArray *muonArr,         // input muons
25 +                                SelectionStatus (*MuonPreSelector)( ControlFlags &, const mithep::TMuon*),
26 +                                SelectionStatus (*MuonIDSelector)( ControlFlags &, const mithep::TMuon*),
27 +                                SelectionStatus (*MuonIsoSelector)( ControlFlags &, const mithep::TMuon*) )
28 +  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
29 + {      
30  
31 +  EventData ret;
32 +  unsigned evtfail = 0x0;
33 +  TRandom3 r;
34 +  
35    if( ctrl.debug ) {
36      cout << "Run: " << info->runNum
37           << "\tEvt: " << info->evtNum
# Line 64 | Line 44 | unsigned fails_HZZ4L_selection(ControlFl
44      RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);      
45      if( !(rlrm.HasRunLumi(rl)) )  {
46        if( ctrl.debug ) cout << "\tfails JSON" << endl;
47 <      evtfail |= (1<<EVTFAIL_JSON);
48 <      return evtfail;
47 >      ret.status.setStatus(0);
48 >      return ret;
49      }
50    }
51    
52    
73
74  
75  
53    //********************************************************
54    // Trigger
55    //********************************************************
56    if( !ctrl.mc ) {
57 <    //  if( !(passHLT(info->triggerBits, info->runNum, channel) )  ) {
57 >    //if( !(passHLT(info->triggerBits, info->runNum, channel) )  ) {
58      if( !(passHLT(info->triggerBits, info->runNum, 999) )  ) {
59 +      if( ctrl.debug ) cout << "\tfails trigger" << endl;
60        evtfail |= (1<<EVTFAIL_TRIGGER);
61 <      return evtfail;
61 >      ret.status.setStatus(0);
62 >      return ret;
63      }  
64 <  } else {
86 <    if( !(passHLTMC(info->triggerBits)) ) {
87 <      evtfail |= (1<<EVTFAIL_TRIGGER);
88 <      return evtfail;
89 <    }      
90 <    //    cout << "MC trigger bits: " << hex << info->triggerBits << dec << endl;
91 <  }
92 <
64 >  }
65    if( ctrl.debug ) {
66      cout << "presel nlep: " << muonArr->GetEntries() + electronArr->GetEntries()
67           << "\tnmuon: "    << muonArr->GetEntries()
# Line 105 | Line 77 | unsigned fails_HZZ4L_selection(ControlFl
77    //    
78    if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
79    //----------------------------------------------------
80 <  for(Int_t i=0; i<muonArr->GetEntries(); i++) {
81 <    const mithep::TMuon *mu = (mithep::TMuon*)((*muonArr)[i]);      
82 <    unsigned muonfail;
83 <    if( ctrl.muSele == "ksWW" )
84 <      muonfail = passMuonSelection(mu);
85 <    else
86 <      muonfail = passMuonSelectionZZ(mu);
87 <    if( ctrl.debug ) {
88 <      cout << "muon:: pt: " << mu->pt
89 <           << "\teta: " << mu->eta
90 <           << "\tmask: 0x" << hex << muonfail << dec
91 <           << endl;
92 <    }
93 <    if ( !muonfail ) {
94 <      SimpleLepton tmplep;
95 <      tmplep.vec.SetPtEtaPhiM(mu->pt,
96 <                              mu->eta,
97 <                              mu->phi,
98 <                              105.658369e-3);
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 <      lepvec.push_back(tmplep);
139 <      if( ctrl.debug ) { cout << "muon passes ... " << endl;}
140 <    }
141 <  }
142 <  
143 <  if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
144 <
145 <  //----------------------------------------------------
146 <  for(Int_t i=0; i<electronArr->GetEntries(); i++) {
147 <    const mithep::TElectron *ele = (mithep::TElectron*)((*electronArr)[i]);
148 <
149 <    Bool_t isMuonOverlap = kFALSE;
150 <    for (int k=0; k<lepvec.size(); ++k) {
151 <      TVector3 tmplep;
152 <      tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
153 <      if ( lepvec[k].type == 13 && lepvec[k].vec.Vect().DrEtaPhi(tmplep) < 0.1 ) {
154 <        if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
155 <        isMuonOverlap = kTRUE;
156 <        break;
80 >  for(Int_t i=0; i<muonArr->GetEntries(); i++)
81 >    {
82 >      const mithep::TMuon *mu = (mithep::TMuon*)((*muonArr)[i]);      
83 >      
84 >      SelectionStatus musel;
85 >      if(ctrl.debug) cout << "musel.status  before anything: " << musel.getStatus() << endl;
86 >      musel |= (*MuonPreSelector)(ctrl,mu);
87 >      if(ctrl.debug) cout << "musel.status  after presel: " << musel.getStatus() << endl;
88 >      musel |= (*MuonIDSelector)(ctrl,mu);
89 >      if(ctrl.debug) cout << "musel.status  after ID: " << musel.getStatus() << endl;
90 >      musel |= (*MuonIsoSelector)(ctrl,mu);
91 >      if(ctrl.debug) cout << "musel.status  after iso: " << musel.getStatus() << endl;
92 >
93 >      if( ctrl.debug ) {
94 >        cout << "muon:: pt: " << mu->pt
95 >             << "\teta: " << mu->eta
96 >             << "\tisorel: " <<  mu->pfIso03/mu->pt
97 >             << "\tstatus: " << hex << musel.getStatus() << dec
98 >             << endl;
99        }
100 <    }
101 <
160 <    unsigned FAIL=0, isEleTight=0;
100 >      
101 >      if ( musel.pass() ) {
102  
103 <    unsigned  failsCIC=0;
104 <    CICStruct ciccuts_tight, ciccuts_medium, ciccuts_loose;
105 <    if(ctrl.eleSele=="cic") {
106 <      if( ctrl.eleSeleScheme == "mediumloose" ) {
107 <        ciccuts_medium = getCiCCuts("medium");
108 <        unsigned failsCICMedium = failsCicSelection(ctrl, ele, ciccuts_medium, ctrl.kinematics);
109 <        ciccuts_loose = getCiCCuts("loose");
110 <        unsigned failsCICLoose = failsCicSelection(ctrl, ele, ciccuts_loose, ctrl.kinematics);
111 <        failsCIC = ( failsCICLoose | failsCICMedium );
112 <        if( !failsCICMedium ) isEleTight=1;
113 <      }
114 <      else {
115 <        ciccuts_tight = getCiCCuts(ctrl.eleSeleScheme);
116 <        failsCIC = failsCicSelection(ctrl, ele, ciccuts_tight, ctrl.kinematics);
117 <        if( !failsCIC ) isEleTight=1;
103 >        SimpleLepton tmplep;
104 >        float pt = mu->pt;
105 >        tmplep.vecorig->SetPtEtaPhiM(pt,
106 >                                    mu->eta,
107 >                                    mu->phi,
108 >                                    MUON_MASS);
109 >        
110 >        if( ctrl.do_escale_up ) {
111 >          pt=scale_smear_muon_Up(pt, 1,  r);
112 >        }
113 >        if( ctrl.do_escale_down ) {
114 >          pt=scale_smear_muon_Down(pt, 1,  r);
115 >        }
116 >        
117 >        tmplep.vec->SetPtEtaPhiM(pt,
118 >                                mu->eta,
119 >                                mu->phi,
120 >                                MUON_MASS);
121 >        
122 >        tmplep.type    = 13;
123 >        tmplep.index   = i;
124 >        tmplep.charge  = mu->q;
125 >        tmplep.isoTrk  = mu->trkIso03;
126 >        tmplep.isoEcal = mu->emIso03;
127 >        tmplep.isoHcal = mu->hadIso03;
128 >        tmplep.isoPF03 = mu->pfIso03;
129 >        tmplep.isoPF04 = mu->pfIso04;
130 >        tmplep.ip3dSig = mu->ip3dSig;
131 >        tmplep.is4l    = false;
132 >        tmplep.isEB    = (fabs(mu->eta) < 1.479 ? 1 : 0 );
133 >        tmplep.isTight = musel.tight();
134 >        tmplep.isLoose = musel.loose();
135 >        lepvec.push_back(tmplep);
136 >        if( ctrl.debug ) { cout << "muon passes ... " << endl;}
137        }
138 <      FAIL = failsCIC;
139 <    }
138 >      //  }
139 >    }    
140 >  
141 >  
142 >    
143 >    //
144 >    if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
145 >    // --------------------------------------------------------------------------------
146 >    for(Int_t i=0; i<electronArr->GetEntries(); i++)
147 >      {
148 >        const mithep::TElectron *ele = (mithep::TElectron*)((*electronArr)[i]);
149 >        
150 >        Bool_t isMuonOverlap = kFALSE;
151 >        for (int k=0; k<lepvec.size(); ++k) {
152 >          TVector3 tmplep;
153 >          tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
154 >          if ( lepvec[k].isLoose && lepvec[k].type == 13 && lepvec[k].vec->Vect().DrEtaPhi(tmplep) < 0.1 ) {
155 >            if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
156 >            isMuonOverlap = kTRUE;
157 >            break;
158 >          }
159 >        }
160 >        
161 >        SelectionStatus elesel;
162 >        if( ctrl.debug ) cout << "--> status before anything: " << hex << elesel.getStatus() << dec << endl;
163 >        elesel |= (*ElectronPreSelector)(ctrl,ele);
164 >        if( ctrl.debug ) cout << "--> status after presel: " << hex << elesel.getStatus() << dec << endl;
165 >        elesel |= (*ElectronIDSelector)(ctrl,ele);
166 >        if( ctrl.debug ) cout << "--> status after ID: " << hex << elesel.getStatus() << dec << endl;
167 >        elesel |= (*ElectronIsoSelector)(ctrl,ele);
168 >        if( ctrl.debug ) cout << "--> status after iso: " << hex << elesel.getStatus() << dec << endl;
169 >        
170 >        if( ctrl.debug ){
171 >          cout << "\tscEt: " << ele->scEt
172 >               << "\tscEta: " << ele->scEta
173 >               << "\tncluster: " << ele->ncluster
174 >               << "\tisorel: " <<  ele->pfIso04/ele->pt
175 >               << "\tstatus: " << hex << elesel.getStatus() << dec
176 >               << endl;
177 >        }
178  
179 <    LikStruct likcuts;
180 <    unsigned failsLike=0;
181 <    if(ctrl.eleSele=="lik") {
182 <      likcuts = getLikCuts(ctrl.eleSeleScheme);
183 <      failsLike = failsLikelihoodSelection(ele, likcuts, ctrl.kinematics);
184 <      FAIL = failsLike;
185 <    }
186 <    unsigned failsBDT=0;
187 <    if(ctrl.eleSele=="bdt") {
188 <      if( ctrl.eleSeleScheme == "mediumloose" ) {
189 <        unsigned failsBDTMedium = failsBDTSelection(ctrl,"medium",ele);
190 <        unsigned failsBDTLoose  = failsBDTSelection(ctrl,"loose",ele);
191 <        failsBDT = ( failsBDTLoose | failsBDTMedium );
192 <        if( !failsBDTMedium ) isEleTight=1;
193 <      } else {
194 <        failsBDT = failsBDTSelection(ctrl,"tight",ele);
195 <        if( !failsBDT ) isEleTight=1;
179 >        if ( elesel.pass() && !isMuonOverlap )
180 >          {
181 >            SimpleLepton tmplep;
182 >            
183 >            float pt = ele->pt;
184 >            tmplep.vecorig->SetPtEtaPhiM( pt,
185 >                                          ele->eta,
186 >                                          ele->phi,
187 >                                          ELECTRON_MASS );
188 >            
189 >            if( ctrl.do_escale ) {
190 >              pt=scale_smear_electron(pt, ele->isEB, r);
191 >            }
192 >            if( ctrl.do_escale_up ) {
193 >              pt=scale_smear_electron_Up(pt, ele->isEB,  r);
194 >            }
195 >            if( ctrl.do_escale_down ) {
196 >              pt=scale_smear_electron_Down(pt, ele->isEB,  r);
197 >            }
198 >
199 >            
200 >            tmplep.vec->SetPtEtaPhiM( pt,
201 >                                      ele->eta,
202 >                                      ele->phi,
203 >                                      ELECTRON_MASS );
204 >            
205 >            tmplep.type    = 11;
206 >            tmplep.index   = i;
207 >            tmplep.charge  = ele->q;
208 >            tmplep.isoTrk  = ele->trkIso03;
209 >            tmplep.isoEcal = ele->emIso03;
210 >            tmplep.isoHcal = ele->hadIso03;
211 >            tmplep.isoPF03 = ele->pfIso03;
212 >            tmplep.isoPF04 = ele->pfIso04;
213 >            tmplep.ip3dSig = ele->ip3dSig;
214 >            tmplep.is4l    = false;
215 >            tmplep.isEB    = ele->isEB;
216 >            tmplep.scID    = ele->scID;
217 >            tmplep.isTight = elesel.tight();
218 >            tmplep.isLoose = elesel.loose();
219 >            lepvec.push_back(tmplep);
220 >            if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
221 >          }
222        }
223 <      FAIL = failsBDT;
224 <    }
225 <    unsigned  failsSi=0;
226 <    if(ctrl.eleSele=="si") {
227 <      failsSi = failsSiMVAElectronSelection(ctrl, ele, 0.95, ctrl.kinematics);
228 <      FAIL = failsSi;
223 >    
224 >    
225 >    //********************************************************
226 >    // Dump Stuff
227 >    //********************************************************
228 >    sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
229 >    int nmu=0, nele=0;
230 >    for( int i=0; i<lepvec.size(); i++ ) {
231 >      if(ctrl.debug) cout << "lepvec :: index: " << i
232 >                          << "\tpt: " << lepvec[i].vec->Pt()
233 >                          << "\ttype: " << lepvec[i].type
234 >                          << endl;
235 >      if( abs(lepvec[i].type) == 11 ) nele++;
236 >      else nmu++;
237      }
238 <
239 <
240 <
241 <
210 <    if( ctrl.debug ){
211 <      cout << "CIC category: " << cicCategory(ele)
212 <           << "\tlikelihood: " << ele->likelihood
213 <           << "\tFAIL:  0x"     << hex << FAIL      << dec
214 <           << "\tfailsCIC:  0x" << hex << failsCIC  << dec
215 <           << "\tfailsLike: 0x" << hex << failsLike << dec
216 <           << "\tfailsBDT:  0x" << hex << failsBDT << dec
217 <           << "\tscEt: " << ele->scEt
218 <           << "\tscEta: " << ele->scEta
219 <           << "\tncluster: " << ele->ncluster
238 >    if( ctrl.debug ) {
239 >      cout << "postsel nlep: " << lepvec.size()
240 >           << "\tnmuon: " << nmu
241 >           << "\tnelectron: " << nele
242             << endl;
243      }
244 <    if ( !FAIL && !isMuonOverlap ) {
245 <      SimpleLepton tmplep;
246 <      tmplep.vec.SetPtEtaPhiM( ele->pt,
247 <                               ele->eta,
248 <                               ele->phi,
249 <                               0.51099892e-3 );
250 <      tmplep.type    = 11;
251 <      tmplep.index   = i;
252 <      tmplep.charge  = ele->q;
253 <      tmplep.isoTrk  = ele->trkIso03;
254 <      tmplep.isoEcal = ele->emIso03;
255 <      tmplep.isoHcal = ele->hadIso03;
256 <      tmplep.isoPF03 = ele->pfIso03;
257 <      tmplep.isoPF04 = ele->pfIso04;
258 <      tmplep.ip3dSig = ele->ip3dSig;
259 <      tmplep.is4l    = false;
260 <      tmplep.isTight = isEleTight;
261 <      tmplep.isEB    = ele->isEB;
262 <      lepvec.push_back(tmplep);
263 <      if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
244 >    
245 >    
246 >    //******************************************************************************
247 >    // Z1 Selection
248 >    //******************************************************************************
249 >    int Z1LeptonPlusIndex = -1;
250 >    int Z1LeptonMinusIndex = -1;
251 >    double BestZ1Mass = -999;
252 >    if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
253 >    for(int i = 0; i < lepvec.size(); ++i) {
254 >      if( !(lepvec[i].isLoose) ) continue;
255 >      for(int j = i+1; j < lepvec.size(); ++j) {
256 >        if( !(lepvec[j].isLoose) ) continue;
257 >        if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
258 >        if (!(lepvec[i].vec->Pt() > 20.0 || lepvec[j].vec->Pt() > 20.0)) continue;
259 >        if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
260 >        if (!(lepvec[i].vec->Pt() > 10.0 && lepvec[j].vec->Pt() > 10.0)) continue;
261 >        if( ctrl.debug ) { cout << "\tthe other  is > 10 GeV" << endl; }
262 >        if (lepvec[i].charge == lepvec[j].charge) continue;          
263 >        if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
264 >        if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;  
265 >        if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
266 >        
267 >        //Make Z1 hypothesis
268 >        TLorentzVector *leptonPlus, *leptonMinus;
269 >        if ( lepvec[i].charge > 0 ) {
270 >          leptonPlus  = lepvec[i].vec;
271 >          leptonMinus = lepvec[j].vec;
272 >        } else {
273 >          leptonPlus  = lepvec[j].vec;
274 >          leptonMinus = lepvec[i].vec;
275 >        }
276 >        
277 >        float tmpZ1Mass = (*leptonPlus + *leptonMinus).M();
278 >        if( ctrl.debug ) cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
279 >        if( tmpZ1Mass > 60 ) {
280 >          if (fabs(tmpZ1Mass - Z_MASS) < fabs(BestZ1Mass - Z_MASS)) {
281 >            BestZ1Mass = tmpZ1Mass;
282 >            if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass
283 >                                  << "\tdM: " << fabs(BestZ1Mass - Z_MASS)
284 >                                  << endl;
285 >            if (lepvec[i].charge > 0) {
286 >              Z1LeptonPlusIndex = i;
287 >              Z1LeptonMinusIndex = j;
288 >            } else {
289 >              Z1LeptonPlusIndex = j;
290 >              Z1LeptonMinusIndex = i;
291 >            }
292 >          }
293 >        }
294 >      }
295      }
296 <  }
297 <  
298 <  sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
299 <
300 <  int nmu=0, nele=0;
301 <  for( int i=0; i<lepvec.size(); i++ ) {
302 <    if( abs(lepvec[i].type) == 11 ) nele++;
303 <    else nmu++;
304 <  }
305 <  if( ctrl.debug ) {
306 <    cout << "postsel nlep: " << lepvec.size()
307 <         << "\tnmuon: " << nmu
308 <         << "\tnelectron: " << nele
309 <         << endl;
310 <  }
311 <
312 <  //******************************************************************************
313 <  //Z1 Selection
314 <  //******************************************************************************
315 <  int Z1LeptonPlusIndex = -1;
316 <  int Z1LeptonMinusIndex = -1;
317 <  double BestZ1Mass = -999;
318 <  if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
319 <  for(int i = 0; i < lepvec.size(); ++i) {
320 <    for(int j = i+1; j < lepvec.size(); ++j) {
321 <      if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
322 <      if (!(lepvec[i].vec.Pt() > 20.0 || lepvec[j].vec.Pt() > 20.0)) continue;
323 <      if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
324 <      if (!(lepvec[i].vec.Pt() > 10.0 && lepvec[j].vec.Pt() > 10.0)) continue;
325 <      if( ctrl.debug ) { cout << "\tthe other  is > 10 GeV" << endl; }
326 <      if (lepvec[i].charge == lepvec[j].charge) continue;          
327 <      if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
328 <      if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;  
329 <      if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
296 >    // stop if no Z1 candidate is found
297 >    if( BestZ1Mass < 0 ) {
298 >      evtfail |= (1<<EVTFAIL_Z1);
299 >      //ret.status = evtfail;
300 >      ret.status.setStatus(0);
301 >      return ret;
302 >    }
303 >    if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum  << endl;
304 >    if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
305 >    TLorentzVector Z1LeptonPlus  = *(lepvec[Z1LeptonPlusIndex].vec);
306 >    TLorentzVector Z1LeptonMinus = *(lepvec[Z1LeptonMinusIndex].vec);
307 >    TLorentzVector Z1Candidate   =  Z1LeptonPlus + Z1LeptonMinus;
308 >    
309 >    
310 >    //******************************************************************************
311 >    // Z1 + l
312 >    //******************************************************************************
313 >    if( lepvec.size() < 3 ) {
314 >      evtfail |= (1<<EVTFAIL_Z1_PLUSL);
315 >      //ret.status = evtfail;
316 >      ret.status.setStatus(0);
317 >      return ret;
318 >    }
319 >    
320 >    //******************************************************************************
321 >    // 4l/Z2 Selection
322 >    //******************************************************************************
323 >    Int_t Z2LeptonPlusIndex = -1;
324 >    Int_t Z2LeptonMinusIndex = -1;
325 >    Double_t BestZ2Mass = -1;
326 >    if( ctrl.debug ) cout << "looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
327 >    for(int i = 0; i < lepvec.size(); ++i) {
328 >      
329 >      if( ctrl.debug)  cout << "i: " << i
330 >                            << "\tpt: " << lepvec[i].vec->Pt()
331 >                            << "\ttype: " << lepvec[i].type
332 >                            << endl;
333        
334 <      //Make Z1 hypothesis
335 <      TLorentzVector leptonPlus, leptonMinus;
336 <      if ( lepvec[i].charge > 0 ) {
337 <        leptonPlus  = lepvec[i].vec;
282 <        leptonMinus = lepvec[j].vec;
283 <      } else {
284 <        leptonPlus  = lepvec[j].vec;
285 <        leptonMinus = lepvec[i].vec;
334 >      if( ctrl.eleSeleScheme == "mediumloose" &&
335 >          !(lepvec[i].isTight) ) {
336 >        if( ctrl.debug)  cout << "it's not tight, skipping ... " << endl;
337 >        continue;
338        }
339        
340 <      float tmpZ1Mass = (leptonPlus+leptonMinus).M();
341 <      if( ctrl.debug ) cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
342 <      if( tmpZ1Mass > 60 ) {
343 <        if (fabs(tmpZ1Mass - 91.1876) < fabs(BestZ1Mass - 91.1876)) {
344 <          BestZ1Mass = tmpZ1Mass;
345 <          if( ctrl.debug ) cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass
346 <               << "\tdM: " << fabs(BestZ1Mass - 91.1876)
347 <               << endl;
340 >      for(int j = i+1; j < lepvec.size(); ++j) {
341 >        if( ctrl.debug)  cout << "\t\tj: " << j
342 >                              << "\tpt: " << lepvec[j].vec->Pt()
343 >                              << "\ttype: " << lepvec[j].type
344 >                              << endl;
345 >        
346 >        if( ctrl.eleSeleScheme == "mediumloose" &&
347 >            !(lepvec[j].isTight) ) {
348 >          if( ctrl.debug)  cout << "it's not tight, skipping ... " << endl;
349 >          continue;
350 >        }
351 >        
352 >        
353 >        if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) {
354 >          if( ctrl.debug)  cout << "\ti matches a Z1 index, skipping ..." << endl;
355 >          continue; //skip Z1 leptons
356 >        }
357 >        if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) {
358 >          if( ctrl.debug)  cout << "\tj matches a Z1 index, skipping ..." << endl;
359 >          continue; //skip Z1 leptons
360 >        }
361 >        if (lepvec[i].charge == lepvec[j].charge) {
362 >          if( ctrl.debug)  cout << "\ti and j are same sign, skipping ..." << endl;
363 >          continue;         //require opp sign
364 >        }
365 >        if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
366 >          if( ctrl.debug) cout << "\ti and j are not same flavor, skipping ..." << endl;
367 >          continue; //require same flavor
368 >        }
369 >        
370 >        
371 >        //Make Z2 hypothesis
372 >        TLorentzVector *leptonPlus, *leptonMinus;
373 >        
374 >        if (lepvec[i].charge > 0 ) {
375 >          leptonPlus  = lepvec[i].vec;
376 >          leptonMinus = lepvec[j].vec;
377 >        } else {
378 >          leptonPlus  = lepvec[j].vec;
379 >          leptonMinus = lepvec[i].vec;
380 >        }
381 >        
382 >        TLorentzVector dilepton = *leptonPlus + *leptonMinus;
383 >        TLorentzVector fourLepton = Z1Candidate + dilepton;
384 >        
385 >        if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
386 >        if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;
387 >        
388 >        if (!(dilepton.M() > 12.0)) continue;
389 >        if (!(fourLepton.M() > 100.0)) continue;
390 >        
391 >        //for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
392 >        if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
393 >          TLorentzVector pair1 = Z1LeptonPlus + *leptonMinus;
394 >          TLorentzVector pair2 = Z1LeptonMinus + *leptonPlus;
395 >          if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
396 >          if (!(pair1.M() > 12 || pair2.M() > 12)) continue;
397 >        }
398 >      
399 >      
400 >        //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
401 >        if (Z2LeptonPlusIndex < 0) {
402            if (lepvec[i].charge > 0) {
403 <            Z1LeptonPlusIndex = i;
404 <            Z1LeptonMinusIndex = j;
403 >            Z2LeptonPlusIndex = i;
404 >            Z2LeptonMinusIndex = j;
405            } else {
406 <            Z1LeptonPlusIndex = j;
407 <            Z1LeptonMinusIndex = i;
406 >            Z2LeptonPlusIndex = j;
407 >            Z2LeptonMinusIndex = i;
408 >          }
409 >        } else {
410 >          Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec->Pt();              
411 >          Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec->Pt();
412 >          if (lepvec[Z2LeptonMinusIndex].vec->Pt() > BestPairPtMax) {
413 >            BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec->Pt();
414 >            BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec->Pt();
415            }
303        }
304      }
305    }
306  }
307  // stop if no Z1 candidate is found
308  if( BestZ1Mass < 0 ) {
309    evtfail |= (1<<EVTFAIL_Z1);
310    return evtfail;
311  }
312  if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum  << endl;
313  if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
314  TLorentzVector Z1LeptonPlus  = lepvec[Z1LeptonPlusIndex].vec;
315  TLorentzVector Z1LeptonMinus = lepvec[Z1LeptonMinusIndex].vec;
316  TLorentzVector Z1Candidate   =  Z1LeptonPlus + Z1LeptonMinus;
317  if( l != NULL ) {
318    l->vecz1 = Z1Candidate;
319    l->vecl1p = Z1LeptonPlus;
320    l->vecl1m = Z1LeptonMinus;
321  }
322  
323  //******************************************************************************
324  // Z1 + l
325  //******************************************************************************
326  if( lepvec.size() < 3 ) {
327    evtfail |= (1<<EVTFAIL_Z1_PLUSL);
328    return evtfail;
329  }
416          
417 <  //******************************************************************************
418 <  // 4l/Z2 Selection
419 <  //******************************************************************************
420 <  Int_t Z2LeptonPlusIndex = -1;
421 <  Int_t Z2LeptonMinusIndex = -1;
422 <  Double_t BestZ2Mass = -1;
423 <  if( ctrl.debug ) cout << "looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
424 <  for(int i = 0; i < lepvec.size(); ++i) {
425 <    if( abs(lepvec[i].type) == 11 &&
426 <        ctrl.eleSeleScheme == "mediumloose" &&
427 <        !(lepvec[i].isTight) ) continue;
342 <
343 <    for(int j = i+1; j < lepvec.size(); ++j) {
344 <      if( abs(lepvec[j].type) == 11 &&
345 <          ctrl.eleSeleScheme == "mediumloose" &&
346 <          !(lepvec[j].isTight) ) continue;
347 <
348 <      //            cout << "i: " << i << "\tj: " << j << endl;
349 <      if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) {
350 <        //            cout << "\ti matches a Z1 index, skipping ..." << endl;
351 <        continue; //skip Z1 leptons
352 <      }
353 <      if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) {
354 <        //            cout << "\tj matches a Z1 index, skipping ..." << endl;
355 <              continue; //skip Z1 leptons
356 <            }
357 <            if (lepvec[i].charge == lepvec[j].charge) {
358 <              //              cout << "\ti and j are same sign, skipping ..." << endl;
359 <              continue;         //require opp sign
360 <            }
361 <            if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
362 <              //              cout << "\ti and j are not same flavor, skipping ..." << endl;
363 <              continue; //require same flavor
364 <            }
365 <            
366 <            
367 <            //Make Z2 hypothesis
368 <            TLorentzVector leptonPlus, leptonMinus;
369 <            
370 <            if (lepvec[i].charge > 0 ) {
371 <              leptonPlus  = lepvec[i].vec;
372 <              leptonMinus = lepvec[j].vec;
373 <            } else {
374 <              leptonPlus  = lepvec[j].vec;
375 <              leptonMinus = lepvec[i].vec;
376 <            }
377 <
378 <            TLorentzVector dilepton = leptonPlus+leptonMinus;
379 <            TLorentzVector fourLepton = Z1Candidate + dilepton;
380 <
381 <            if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
382 <            if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;
383 <
384 <            if (!(dilepton.M() > 12.0)) continue;
385 <            if (!(fourLepton.M() > 100.0)) continue;
386 <
387 <            /*
388 <            //for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
389 <            if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
390 <              TLorentzVector pair1 = Z1LeptonPlus+leptonMinus;
391 <              TLorentzVector pair2 = Z1LeptonMinus+leptonPlus;
392 <              if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
393 <              if (!(pair1.M() > 12 || pair2.M() > 12)) continue;
394 <            }
395 <            */
396 <
397 <            //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
398 <            if (Z2LeptonPlusIndex < 0) {
399 <              if (lepvec[i].charge > 0) {
400 <                Z2LeptonPlusIndex = i;
401 <                Z2LeptonMinusIndex = j;
402 <              } else {
403 <                Z2LeptonPlusIndex = j;
404 <                Z2LeptonMinusIndex = i;
405 <              }
406 <            } else {
407 <              Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec.Pt();              
408 <              Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec.Pt();
409 <              if (lepvec[Z2LeptonMinusIndex].vec.Pt() > BestPairPtMax) {
410 <                BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec.Pt();
411 <                BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec.Pt();
412 <              }
413 <
414 <              Double_t CurrentPairPtMax = lepvec[i].vec.Pt();              
415 <              Double_t CurrentPairPtMin = lepvec[j].vec.Pt();
416 <              if (lepvec[j].vec.Pt() > CurrentPairPtMax) {
417 <                CurrentPairPtMax = lepvec[j].vec.Pt();
418 <                CurrentPairPtMin = lepvec[i].vec.Pt();
419 <              }
420 <
421 <              if (CurrentPairPtMax > BestPairPtMax) {
422 <                if (lepvec[i].charge > 0) {
423 <                  Z2LeptonPlusIndex = i;
424 <                  Z2LeptonMinusIndex = j;
425 <                } else {
426 <                  Z2LeptonPlusIndex = j;
427 <                  Z2LeptonMinusIndex = i;
428 <                }
429 <              } else if (CurrentPairPtMax  == BestPairPtMax) {
430 <                if (CurrentPairPtMin > BestPairPtMin) {
431 <                  if (lepvec[i].charge > 0) {
432 <                    Z2LeptonPlusIndex = i;
433 <                    Z2LeptonMinusIndex = j;
434 <                  } else {
435 <                    Z2LeptonPlusIndex = j;
436 <                    Z2LeptonMinusIndex = i;
437 <                  }                  
438 <                }
439 <              }
440 <            }            
441 <          }
442 <        }
443 <
444 <        // stop if no Z2 candidate is found
445 <        if (Z2LeptonPlusIndex == -1) {
446 <          evtfail |= ( 1<<EVTFAIL_4L );
447 <          return evtfail;
448 <          //      h_evtfail->Fill( evtfail );
449 <          //      cout << "evtfail: " << hex << evtfail << dec << endl;
450 <          //      continue;
451 <        }
452 <        if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
453 <        if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
454 <                              << "\tminusindex: " << Z2LeptonMinusIndex << endl;
455 <        TLorentzVector  Z2LeptonPlus  = lepvec[Z2LeptonPlusIndex].vec;
456 <        TLorentzVector  Z2LeptonMinus = lepvec[Z2LeptonMinusIndex].vec;
457 <        TLorentzVector  Z2Candidate   = Z2LeptonPlus+Z2LeptonMinus;
458 <        TLorentzVector  ZZSystem      = Z1Candidate + Z2Candidate;
459 <        if( l != NULL ) {
460 <          l->vecz2 = Z2Candidate;
461 <          l->vecl2p = Z2LeptonPlus;
462 <          l->vecl2m = Z2LeptonMinus;
463 <          l->vec4l = ZZSystem;
464 <        }
465 <        lepvec[Z1LeptonPlusIndex].is4l = true;
466 <        lepvec[Z1LeptonMinusIndex].is4l = true;
467 <        lepvec[Z2LeptonPlusIndex].is4l = true;
468 <        lepvec[Z2LeptonMinusIndex].is4l = true;
469 <
470 <        //***************************************************************
471 <        // Isolation
472 <        //***************************************************************
473 <        bool failiso=false;
474 <
475 <        if( ctrl.isoScheme == "pf" ) {
476 <
477 <          for( int i=0; i<lepvec.size(); i++ ) {
478 <
479 <            if( !(lepvec[i].is4l) ) continue;
480 <
481 <            if( abs(lepvec[i].type) == 11 ) {
482 <              float reliso = lepvec[i].isoPF04/lepvec[i].vec.Pt();
483 <              if( lepvec[i].isEB && lepvec[i].vec.Pt() > 20 && reliso > PFISO_ELE_LOOSE_EB_HIGHPT ) {
484 <                failiso = true;
485 <                break;
486 <              }
487 <              if( lepvec[i].isEB && lepvec[i].vec.Pt() < 20 && reliso > PFISO_ELE_LOOSE_EB_LOWPT ) {
488 <                failiso = true;
489 <                break;
490 <              }
491 <              if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() > 20 && reliso > PFISO_ELE_LOOSE_EE_HIGHPT ) {
492 <                failiso = true;
493 <                break;
494 <              }
495 <              if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() < 20 && reliso > PFISO_ELE_LOOSE_EE_LOWPT ) {
496 <                failiso = true;
497 <                break;
498 <              }
499 <            }
500 <
501 <            if( abs(lepvec[i].type) == 13 ) {
502 <
503 <              float reliso = lepvec[i].isoPF03/lepvec[i].vec.Pt();
504 <              if( lepvec[i].isEB && lepvec[i].vec.Pt() > 20 && reliso > PFISO_MU_TIGHT_EB_HIGHPT ) {  //0.13
505 <                failiso = true;
506 <                break;
507 <              }
508 <              if( lepvec[i].isEB && lepvec[i].vec.Pt() < 20 && reliso > PFISO_MU_TIGHT_EB_LOWPT ) { //0.06
509 <                failiso = true;
510 <                break;
511 <              }
512 <              if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() > 20 && reliso > PFISO_MU_TIGHT_EE_HIGHPT ) { //0.09
513 <                failiso = true;
514 <                break;
515 <              }
516 <              if( !(lepvec[i].isEB) && lepvec[i].vec.Pt() < 20 && reliso > PFISO_MU_TIGHT_EE_LOWPT ) { //0.05
517 <                failiso = true;
518 <                break;
519 <              }
520 <            }
521 <          }  
522 <        } else if ( ctrl.isoScheme == "pairwise" ) {
523 <          float rho = info->rho;
524 <          for( int i=0; i<lepvec.size(); i++ ) {
525 <            if( !(lepvec[i].is4l) ) continue;
526 <            float effArea_ecal_i, effArea_hcal_i;
527 <            if( lepvec[i].isEB ) {
528 <              if( lepvec[i].type == 11 ) {
529 <                effArea_ecal_i = 0.101;
530 <                effArea_hcal_i = 0.021;
531 <              } else {
532 <                effArea_ecal_i = 0.074;
533 <                effArea_hcal_i = 0.022;
534 <              }
417 >          Double_t CurrentPairPtMax = lepvec[i].vec->Pt();              
418 >          Double_t CurrentPairPtMin = lepvec[j].vec->Pt();
419 >          if (lepvec[j].vec->Pt() > CurrentPairPtMax) {
420 >            CurrentPairPtMax = lepvec[j].vec->Pt();
421 >            CurrentPairPtMin = lepvec[i].vec->Pt();
422 >          }
423 >        
424 >          if (CurrentPairPtMax > BestPairPtMax) {
425 >            if (lepvec[i].charge > 0) {
426 >              Z2LeptonPlusIndex = i;
427 >              Z2LeptonMinusIndex = j;
428              } else {
429 <              if( lepvec[i].type == 11 ) {
430 <                effArea_ecal_i = 0.046;
538 <                effArea_hcal_i = 0.040;
539 <              } else {
540 <                effArea_ecal_i = 0.045;
541 <                effArea_hcal_i = 0.030;
542 <              }
429 >              Z2LeptonPlusIndex = j;
430 >              Z2LeptonMinusIndex = i;
431              }
432 <            float isoEcal_corr_i = lepvec[i].isoEcal - (effArea_ecal_i*rho);
433 <            float isoHcal_corr_i = lepvec[i].isoHcal - (effArea_hcal_i*rho);
434 <            for( int j=i+1; j<lepvec.size(); j++ ) {
435 <              if( !(lepvec[j].is4l) ) continue;
436 <              float effArea_ecal_j, effArea_hcal_j;
549 <              if( lepvec[j].isEB ) {
550 <                if( lepvec[j].type == 11 ) {
551 <                  effArea_ecal_j = 0.101;
552 <                  effArea_hcal_j = 0.021;
553 <                } else {
554 <                  effArea_ecal_j = 0.074;
555 <                  effArea_hcal_j = 0.022;
556 <                }
432 >          } else if (CurrentPairPtMax  == BestPairPtMax) {
433 >            if (CurrentPairPtMin > BestPairPtMin) {
434 >              if (lepvec[i].charge > 0) {
435 >                Z2LeptonPlusIndex = i;
436 >                Z2LeptonMinusIndex = j;
437                } else {
438 <                if( lepvec[j].type == 11 ) {
439 <                  effArea_ecal_j = 0.046;
440 <                  effArea_hcal_j = 0.040;
561 <                } else {
562 <                  effArea_ecal_j = 0.045;
563 <                  effArea_hcal_j = 0.030;
564 <                }
565 <              }
566 <              float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
567 <              float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
568 <              float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec.Pt();
569 <              float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec.Pt();          
570 <              float comboIso = RIso_i + RIso_j;
571 <              if( info->evtNum == 1038911933 ) {
572 <                float tmpdR = lepvec[i].vec.DrEtaPhi(lepvec[j].vec);
573 <                cout << "i: " << i
574 <                     << "\tdR: " << tmpdR
575 <                     << "\trho: " << rho
576 <                     << "\tRIso_i: " << RIso_i  
577 <                     << "\ttkrel: " << lepvec[i].isoTrk/lepvec[i].vec.Pt()
578 <                     << "\tecalrel: " << lepvec[i].isoEcal/lepvec[i].vec.Pt()
579 <                     << "\tecalrelcor: " << isoEcal_corr_i/lepvec[i].vec.Pt()
580 <                     << "\thcalrel: " << lepvec[i].isoHcal/lepvec[i].vec.Pt()
581 <                     << "\thcalrelcor: " << isoHcal_corr_i/lepvec[i].vec.Pt()
582 <                     << "\tpt_i: " << lepvec[i].vec.Pt()
583 <                     << "\tj: " << j
584 <                     << "\tRIso_j: " << RIso_j  
585 <                     << "\ttkrel: "  << lepvec[j].isoTrk/lepvec[j].vec.Pt()
586 <                     << "\tecalrel: " << lepvec[j].isoEcal/lepvec[j].vec.Pt()
587 <                     << "\tecalrelcor: " << isoEcal_corr_j/lepvec[j].vec.Pt()
588 <                     << "\thcalrel: " << lepvec[j].isoHcal/lepvec[j].vec.Pt()
589 <                     << "\thcalrelcor: " << isoHcal_corr_j/lepvec[j].vec.Pt()
590 <                     << "\tpt_j: " << lepvec[j].vec.Pt()
591 <                     << "\tcombo: " << comboIso
592 <                     << endl;
593 <                cout.flush();
594 <              }
595 <              if( comboIso > 0.35 ) {
596 <                if( ctrl.debug ) cout << "combo failing for indices: " << i << "," << j << endl;
597 <                failiso = true;
598 <                //            break;
599 <              }
438 >                Z2LeptonPlusIndex = j;
439 >                Z2LeptonMinusIndex = i;
440 >              }                  
441              }
442            }
443 <        }
444 <
443 >        }            
444 >      }
445 >    }
446 >  
447 >    // stop if no Z2 candidate is found
448 >    if (Z2LeptonPlusIndex == -1) {
449 >      evtfail |= ( 1<<EVTFAIL_4L );
450 >      //      ret.status = evtfail;
451 >      ret.status.setStatus(0);
452 >      return ret;
453 >    }
454 >    if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
455 >    if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
456 >                          << "\tminusindex: " << Z2LeptonMinusIndex << endl;
457 >    TLorentzVector  Z2LeptonPlus  = *(lepvec[Z2LeptonPlusIndex].vec);
458 >    TLorentzVector  Z2LeptonMinus = *(lepvec[Z2LeptonMinusIndex].vec);
459 >    TLorentzVector  Z2Candidate   = Z2LeptonPlus+Z2LeptonMinus;
460 >    TLorentzVector  ZZSystem      = Z1Candidate + Z2Candidate;
461 >    lepvec[Z1LeptonPlusIndex].is4l = true;
462 >    lepvec[Z1LeptonMinusIndex].is4l = true;
463 >    lepvec[Z2LeptonPlusIndex].is4l = true;
464 >    lepvec[Z2LeptonMinusIndex].is4l = true;
465 >  
466  
605        if( failiso ) {
606          evtfail |= ( 1<<EVTFAIL_ISOLATION );
607          return evtfail;
608          //h_evtfail->Fill( evtfail, eventweight );
609          //      h_evtfail->Fill( evtfail );
610          //      cout << "evtfail: " << hex << evtfail << dec << endl;
611          //      continue;
612        }
613
614        //***************************************************************
615        // IP significance
616        //***************************************************************
617        bool failip = false;
618        for( int i=0; i<lepvec.size(); i++ ) {
619          if( !(lepvec[i].is4l) ) continue;
620          if( lepvec[i].ip3dSig > 4 ) {
621            failip=true;
622            break;
623          }
624        }
625        if( failip ) {
626          evtfail |= (1<<EVTFAIL_IP );
627          return evtfail;
628          //h_evtfail->Fill( evtfail, eventweight );
629          //      h_evtfail->Fill( evtfail );
630          //      cout << "evtfail: " << hex << evtfail << dec << endl;
631          //      continue;
632        }
633
634        //***************************************************************
635        // remaining kinematic cuts
636        //***************************************************************
637        double Z2massCut=0;
638        if      ( ctrl.kinematics == "loose" ) Z2massCut = 12;
639        else if ( ctrl.kinematics == "tight" ) Z2massCut = 20;
640        else { cout <<  "error! kinematic tightness not defined!" << endl; assert(0); }
641
642        if ( Z1Candidate.M() > 120        ||
643             Z2Candidate.M() < Z2massCut  ||
644             Z2Candidate.M() > 120        ||
645             !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec.Pt() > 20.0) ||
646             !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec.Pt() > 10.0)
647             ) {
648          evtfail |= (1<<EVTFAIL_KINEMATICS );
649          return evtfail;
650          //h_evtfail->Fill( evtfail, eventweight );
651          //      h_evtfail->Fill( evtfail );
652          //      cout << "evtfail: " << hex << evtfail << dec << endl;
653          //      continue;
654        }
655        
656        int channel;
657        if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
658        if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
659        if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) ||
660            (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
661        
662
663
664        if( passtuple != NULL ) {
665          unsigned run   = info->runNum;
666          unsigned evt   = info->evtNum;
667          unsigned lumi  = info->lumiSec;
668          unsigned chan  = channel;
669          double   w     = eventweight;
670          float mZ1      = Z1Candidate.M() ;
671          float mZ2      = Z2Candidate.M() ;
672          float m4l      = ZZSystem.M() ;
673          float pt4l     = ZZSystem.Pt() ;
674          unsigned tZ1   = abs(lepvec[Z1LeptonMinusIndex].type);
675          unsigned tZ2   = abs(lepvec[Z2LeptonMinusIndex].type);
676          passtuple->SetBranchAddress("run",  &run);
677          passtuple->SetBranchAddress("evt",  &evt);
678          passtuple->SetBranchAddress("lumi", &lumi);
679          passtuple->SetBranchAddress("mZ1",  &mZ1);
680          passtuple->SetBranchAddress("mZ2",  &mZ2);
681          passtuple->SetBranchAddress("tZ1",  &tZ1);
682          passtuple->SetBranchAddress("tZ2",  &tZ2);
683          passtuple->SetBranchAddress("m4l",  &m4l);
684          passtuple->SetBranchAddress("pt4l", &pt4l);
685          passtuple->SetBranchAddress("w",    &w);
686          passtuple->Fill( );
687        }
688
689        if( ctrl.debug ) cout  << "run: " << info->runNum  
690              << "\tevt: " << info->evtNum
691              << "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
692              << "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
693              << "\tmZ1: " << Z1Candidate.M()
694              << "\tmZ2: " << Z2Candidate.M()
695              << "\tm4l: " << ZZSystem.M()
696              << "\tevtfail: " << hex << evtfail << dec
697              << "\ttrigbits: " << hex << info->triggerBits << dec
698          //          << "\ttree: " << inputFiles[q][f]
699              << endl;
467  
468 <        return evtfail;
469 <          
470 < }
468 >  
469 >    //***************************************************************
470 >    // remaining kinematic cuts
471 >    //***************************************************************
472 >    double Z2massCut=0;
473 >    if      ( ctrl.kinematics == "loose" ) Z2massCut = 12;
474 >    else if ( ctrl.kinematics == "tight" ) Z2massCut = 20;
475 >    else { cout <<  "error! kinematic tightness not defined!" << endl; assert(0); }
476 >  
477 >    if ( Z1Candidate.M() > 120        ||
478 >         Z2Candidate.M() < Z2massCut  ||
479 >         Z2Candidate.M() > 120        ||
480 >         !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec->Pt() > 20.0) ||
481 >         !(lepvec[Z1LeptonPlusIndex].vec->Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec->Pt() > 10.0)
482 >         ) {
483 >      evtfail |= (1<<EVTFAIL_KINEMATICS );
484 >      //      ret.status = evtfail;
485 >      ret.status.setStatus(0);
486 >      return ret;
487 >    }
488 >  
489 >    unsigned channel;
490 >    if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
491 >    if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
492 >    if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) ||
493 >        (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
494 >  
495 >  
496 >  
497 >  
498 >    if( ctrl.debug ) cout  << "run: " << info->runNum  
499 >                           << "\tevt: " << info->evtNum
500 >                           << "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
501 >                           << "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
502 >                           << "\tmZ1: " << Z1Candidate.M()
503 >                           << "\tmZ2: " << Z2Candidate.M()
504 >                           << "\tm4l: " << ZZSystem.M()
505 >                           << "\tevtfail: " << hex << evtfail << dec
506 >                           << "\ttrigbits: " << hex << info->triggerBits << dec
507 >      //              << "\ttree: " << inputFiles[q][f]
508 >                           << endl;
509 >  
510  
511  
512 +    //***************************************************************
513 +    // finish
514 +    //***************************************************************
515 +
516 +    if( !evtfail ) {
517 +      ret.status.setStatus(SelectionStatus::EVTPASS);
518 +      ret.Z1leptons.push_back(lepvec[Z1LeptonMinusIndex]);
519 +      ret.Z1leptons.push_back(lepvec[Z1LeptonPlusIndex]);
520 +      ret.Z2leptons.push_back(lepvec[Z2LeptonMinusIndex]);
521 +      ret.Z2leptons.push_back(lepvec[Z2LeptonPlusIndex]);
522 +    }
523  
524 +    return ret;
525 + }
526  
527  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines