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.7 by khahn, Fri Oct 14 11:20:25 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-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));  
17 <  rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));  
18 < };
19 <
20 < unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control
21 <                               mithep::TEventInfo *info,     // input event inof
22 <                               TClonesArray *electronArr,    // input electrons
23 <                               TClonesArray *muonArr,        // input muons
24 <                               double eventweight,           // weight
25 <                               TTree * passtuple ) {
26 <
27 <  fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, passtuple, NULL );
28 <
29 < };
30 <
31 < unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control
32 <                               mithep::TEventInfo *info,     // input event inof
33 <                               TClonesArray *electronArr,    // input electrons
34 <                               TClonesArray *muonArr,        // input muons
35 <                               double eventweight,           // weight
36 <                               LabVectors  *l ) {
37 <
38 <  fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, eventweight, NULL, l );
39 <
40 < };
41 <
42 <
43 < unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control
44 <                               mithep::TEventInfo *info,     // input event inof
45 <                               TClonesArray *electronArr,    // input electrons
46 <                               TClonesArray *muonArr,        // input muons
47 <                               double eventweight,           // weight
48 <                               TTree * passtuple,
49 <                               LabVectors * l) {       // output ntuple
11 > #include "ExternData.h"
12 > #include "SelectionDefs.h"
13  
51  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 63 | 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    
72
73  
74  
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 {
85 <    if( !(passHLTMC(info->triggerBits)) ) {
86 <      evtfail |= (1<<EVTFAIL_TRIGGER);
87 <      return evtfail;
88 <    }      
89 <    //    cout << "MC trigger bits: " << hex << info->triggerBits << dec << endl;
90 <  }
91 <
64 >  }
65    if( ctrl.debug ) {
66      cout << "presel nlep: " << muonArr->GetEntries() + electronArr->GetEntries()
67           << "\tnmuon: "    << muonArr->GetEntries()
# Line 104 | 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 = passMuonSelectionZZ(mu);
83 <    if( ctrl.debug ) {
84 <      cout << "muon:: pt: " << mu->pt
85 <           << "\teta: " << mu->eta
86 <           << "\tmask: 0x" << hex << muonfail << dec
87 <           << endl;
88 <    }
89 <    if ( !muonfail ) {
90 <      SimpleLepton tmplep;
91 <      tmplep.vec.SetPtEtaPhiM(mu->pt,
92 <                              mu->eta,
93 <                              mu->phi,
94 <                              105.658369e-3);
95 <      tmplep.type    = 13;
96 <      tmplep.index   = i;
97 <      tmplep.charge  = mu->q;
98 <      tmplep.isoTrk  = mu->trkIso03;
99 <      tmplep.isoEcal = mu->emIso03;
100 <      tmplep.isoHcal = mu->hadIso03;
101 <      tmplep.ip3dSig = mu->ip3dSig;
102 <      tmplep.is4l    = false;
103 <      tmplep.isEB    = (fabs(mu->eta) < 1.479 ? 1 : 0 );
104 <      lepvec.push_back(tmplep);
105 <      if( ctrl.debug ) { cout << "muon passes ... " << endl;}
106 <    }
107 <  }
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 >      if ( musel.pass() ) {
102 >
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 >      //  }
139 >    }    
140    
141 <  if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
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 <  //----------------------------------------------------
180 <  for(Int_t i=0; i<electronArr->GetEntries(); i++) {
181 <    const mithep::TElectron *ele = (mithep::TElectron*)((*electronArr)[i]);
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 <    Bool_t isMuonOverlap = kFALSE;
200 <    for (int k=0; k<lepvec.size(); ++k) {
201 <      TVector3 tmplep;
202 <      tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
203 <      if ( lepvec[k].type == 13 && lepvec[k].vec.Vect().DrEtaPhi(tmplep) < 0.1 ) {
204 <        if( ctrl.debug ) cout << "-----> isMuonOverlap! " << endl;
205 <        isMuonOverlap = kTRUE;
206 <        break;
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 +    
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 <    unsigned FAIL=0;
240 <    CICStruct ciccuts = getCiCCuts(ctrl.eleSeleScheme);
241 <    unsigned  failsCIC=0;
156 <    if(ctrl.eleSele=="cic") {
157 <      failsCIC = failsCicSelection(ctrl, ele, ciccuts, ctrl.kinematics);
158 <      FAIL = failsCIC;
159 <    }
160 <    LikStruct likcuts;
161 <    unsigned failsLike=0;
162 <    if(ctrl.eleSele=="lik") {
163 <      likcuts = getLikCuts(ctrl.eleSeleScheme);
164 <      failsLike = failsLikelihoodSelection(ele, likcuts, ctrl.kinematics);
165 <      FAIL = failsLike;
166 <    }
167 <    unsigned  failsBDT=0;
168 <    if(ctrl.eleSele=="bdt") {
169 <      failsBDT = failsBDTSelection(ctrl,ele);
170 <      FAIL = failsBDT;
171 <    }
172 <    unsigned  failsSi=0;
173 <    if(ctrl.eleSele=="si") {
174 <      failsSi = failsSiMVAElectronSelection(ctrl, ele, 0.95, ctrl.kinematics);
175 <      FAIL = failsSi;
176 <    }
177 <
178 <
179 <
180 <
181 <    if( ctrl.debug ){
182 <      cout << "CIC category: " << cicCategory(ele)
183 <           << "\tlikelihood: " << ele->likelihood
184 <           << "\tFAIL:  0x"     << hex << FAIL      << dec
185 <           << "\tfailsCIC:  0x" << hex << failsCIC  << dec
186 <           << "\tfailsLike: 0x" << hex << failsLike << dec
187 <           << "\tfailsBDT:  0x" << hex << failsBDT << dec
188 <           << "\tscEt: " << ele->scEt
189 <           << "\tscEta: " << ele->scEta
190 <           << "\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.ip3dSig = ele->ip3dSig;
257 <      tmplep.is4l    = false;
258 <      tmplep.isEB    = ele->isEB;
259 <      lepvec.push_back(tmplep);
260 <      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; }
243 <      if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;  
244 <      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 <      //Make Z1 hypothesis
330 <      TLorentzVector leptonPlus, leptonMinus;
331 <      if ( lepvec[i].charge > 0 ) {
332 <        leptonPlus  = lepvec[i].vec;
333 <        leptonMinus = lepvec[j].vec;
334 <      } else {
335 <        leptonPlus  = lepvec[j].vec;
336 <        leptonMinus = lepvec[i].vec;
329 >      if( ctrl.debug)  cout << "i: " << i
330 >                            << "\tpt: " << lepvec[i].vec->Pt()
331 >                            << "\ttype: " << lepvec[i].type
332 >                            << endl;
333 >      
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            }
271        }
272      }
273    }
274  }
275  // stop if no Z1 candidate is found
276  if( BestZ1Mass < 0 ) {
277    evtfail |= (1<<EVTFAIL_Z1);
278    return evtfail;
279  }
280  if( ctrl.debug ) cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum  << endl;
281  if( ctrl.debug ) cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
282  TLorentzVector Z1LeptonPlus  = lepvec[Z1LeptonPlusIndex].vec;
283  TLorentzVector Z1LeptonMinus = lepvec[Z1LeptonMinusIndex].vec;
284  TLorentzVector Z1Candidate   =  Z1LeptonPlus + Z1LeptonMinus;
285  if( l != NULL ) {
286    l->vecz1 = Z1Candidate;
287    l->vecl1p = Z1LeptonPlus;
288    l->vecl1m = Z1LeptonMinus;
289  }
290  
291  //******************************************************************************
292  // Z1 + l
293  //******************************************************************************
294  if( lepvec.size() < 3 ) {
295    evtfail |= (1<<EVTFAIL_Z1_PLUSL);
296    return evtfail;
297  }
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 <    for(int j = i+1; j < lepvec.size(); ++j) {
426 <      //            cout << "i: " << i << "\tj: " << j << endl;
427 <      if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) {
310 <        //            cout << "\ti matches a Z1 index, skipping ..." << endl;
311 <        continue; //skip Z1 leptons
312 <      }
313 <      if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) {
314 <        //            cout << "\tj matches a Z1 index, skipping ..." << endl;
315 <              continue; //skip Z1 leptons
316 <            }
317 <            if (lepvec[i].charge == lepvec[j].charge) {
318 <              //              cout << "\ti and j are same sign, skipping ..." << endl;
319 <              continue;         //require opp sign
320 <            }
321 <            if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
322 <              //              cout << "\ti and j are not same flavor, skipping ..." << endl;
323 <              continue; //require same flavor
324 <            }
325 <            
326 <            
327 <            //Make Z2 hypothesis
328 <            TLorentzVector leptonPlus, leptonMinus;
329 <            
330 <            if (lepvec[i].charge > 0 ) {
331 <              leptonPlus  = lepvec[i].vec;
332 <              leptonMinus = lepvec[j].vec;
333 <            } else {
334 <              leptonPlus  = lepvec[j].vec;
335 <              leptonMinus = lepvec[i].vec;
336 <            }
337 <
338 <            TLorentzVector dilepton = leptonPlus+leptonMinus;
339 <            TLorentzVector fourLepton = Z1Candidate + dilepton;
340 <
341 <            if( ctrl.debug ) cout << "dilepton.M() : " << dilepton.M() << endl;
342 <            if( ctrl.debug ) cout << "fourLepton.M() : " << fourLepton.M() << endl;
343 <
344 <            if (!(dilepton.M() > 12.0)) continue;
345 <            if (!(fourLepton.M() > 100.0)) continue;
346 <
347 <            
348 <            //for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
349 <            if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
350 <              TLorentzVector pair1 = Z1LeptonPlus+leptonMinus;
351 <              TLorentzVector pair2 = Z1LeptonMinus+leptonPlus;
352 <              if( ctrl.debug ) cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
353 <              if (!(pair1.M() > 12 || pair2.M() > 12)) continue;
354 <            }
355 <            
356 <            //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
357 <            if (Z2LeptonPlusIndex < 0) {
358 <              if (lepvec[i].charge > 0) {
359 <                Z2LeptonPlusIndex = i;
360 <                Z2LeptonMinusIndex = j;
361 <              } else {
362 <                Z2LeptonPlusIndex = j;
363 <                Z2LeptonMinusIndex = i;
364 <              }
365 <            } else {
366 <              Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec.Pt();              
367 <              Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec.Pt();
368 <              if (lepvec[Z2LeptonMinusIndex].vec.Pt() > BestPairPtMax) {
369 <                BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec.Pt();
370 <                BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec.Pt();
371 <              }
372 <
373 <              Double_t CurrentPairPtMax = lepvec[i].vec.Pt();              
374 <              Double_t CurrentPairPtMin = lepvec[j].vec.Pt();
375 <              if (lepvec[j].vec.Pt() > CurrentPairPtMax) {
376 <                CurrentPairPtMax = lepvec[j].vec.Pt();
377 <                CurrentPairPtMin = lepvec[i].vec.Pt();
378 <              }
379 <
380 <              if (CurrentPairPtMax > BestPairPtMax) {
381 <                if (lepvec[i].charge > 0) {
382 <                  Z2LeptonPlusIndex = i;
383 <                  Z2LeptonMinusIndex = j;
384 <                } else {
385 <                  Z2LeptonPlusIndex = j;
386 <                  Z2LeptonMinusIndex = i;
387 <                }
388 <              } else if (CurrentPairPtMax  == BestPairPtMax) {
389 <                if (CurrentPairPtMin > BestPairPtMin) {
390 <                  if (lepvec[i].charge > 0) {
391 <                    Z2LeptonPlusIndex = i;
392 <                    Z2LeptonMinusIndex = j;
393 <                  } else {
394 <                    Z2LeptonPlusIndex = j;
395 <                    Z2LeptonMinusIndex = i;
396 <                  }                  
397 <                }
398 <              }
399 <            }            
400 <          }
401 <        }
402 <
403 <        // stop if no Z2 candidate is found
404 <        if (Z2LeptonPlusIndex == -1) {
405 <          evtfail |= ( 1<<EVTFAIL_4L );
406 <          return evtfail;
407 <          //      h_evtfail->Fill( evtfail );
408 <          //      cout << "evtfail: " << hex << evtfail << dec << endl;
409 <          //      continue;
410 <        }
411 <        if( ctrl.debug ) cout << "\tgot a Z2 ..." << endl;
412 <        if( ctrl.debug ) cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex
413 <                              << "\tminusindex: " << Z2LeptonMinusIndex << endl;
414 <        TLorentzVector  Z2LeptonPlus  = lepvec[Z2LeptonPlusIndex].vec;
415 <        TLorentzVector  Z2LeptonMinus = lepvec[Z2LeptonMinusIndex].vec;
416 <        TLorentzVector  Z2Candidate   = Z2LeptonPlus+Z2LeptonMinus;
417 <        TLorentzVector  ZZSystem      = Z1Candidate + Z2Candidate;
418 <        if( l != NULL ) {
419 <          l->vecz2 = Z2Candidate;
420 <          l->vecl2p = Z2LeptonPlus;
421 <          l->vecl2m = Z2LeptonMinus;
422 <          l->vec4l = ZZSystem;
423 <        }
424 <        lepvec[Z1LeptonPlusIndex].is4l = true;
425 <        lepvec[Z1LeptonMinusIndex].is4l = true;
426 <        lepvec[Z2LeptonPlusIndex].is4l = true;
427 <        lepvec[Z2LeptonMinusIndex].is4l = true;
428 <
429 <        //***************************************************************
430 <        // Isolation
431 <        //***************************************************************
432 <        bool failiso=false;
433 <
434 <        /*
435 <        int i,j;
436 <        i=Z1LeptonPlusIndex;
437 <        j=Z1LeptonMinusIndex;
438 <        float RIso1 = (lepvec[i].isoTrk+lepvec[i].isoEcal+lepvec[i].isoHcal)/lepvec[i].vec.Pt();
439 <        float RIso2 = (lepvec[j].isoTrk+lepvec[j].isoEcal+lepvec[j].isoHcal)/lepvec[j].vec.Pt();
440 <        float comboIso12 = RIso1 + RIso2;
441 <        i=Z2LeptonPlusIndex;
442 <        j=Z2LeptonMinusIndex;
443 <        float RIso3 = (lepvec[i].isoTrk+lepvec[i].isoEcal+lepvec[i].isoHcal)/lepvec[i].vec.Pt();
444 <        float RIso4 = (lepvec[j].isoTrk+lepvec[j].isoEcal+lepvec[j].isoHcal)/lepvec[j].vec.Pt();
445 <        float comboIso34 = RIso3 + RIso4;
446 <        if( comboIso12 > 0.35 || comboIso34 > 0.35 ) {
447 <              failiso = true;
448 <        }
449 <        */
450 <
451 <        float rho = info->rho;
452 <        for( int i=0; i<lepvec.size(); i++ ) {
453 <          if( !(lepvec[i].is4l) ) continue;
454 <            float effArea_ecal_i, effArea_hcal_i;
455 <            if( lepvec[i].isEB ) {
456 <              if( lepvec[i].type == 11 ) {
457 <                effArea_ecal_i = 0.101;
458 <                effArea_hcal_i = 0.021;
459 <              } else {
460 <                effArea_ecal_i = 0.074;
461 <                effArea_hcal_i = 0.022;
462 <              }
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;
466 <                effArea_hcal_i = 0.040;
467 <              } else {
468 <                effArea_ecal_i = 0.045;
469 <                effArea_hcal_i = 0.030;
470 <              }
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;
477 <              if( lepvec[j].isEB ) {
478 <                if( lepvec[j].type == 11 ) {
479 <                  effArea_ecal_j = 0.101;
480 <                  effArea_hcal_j = 0.021;
481 <                } else {
482 <                  effArea_ecal_j = 0.074;
483 <                  effArea_hcal_j = 0.022;
484 <                }
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;
489 <                } else {
490 <                  effArea_ecal_j = 0.045;
491 <                  effArea_hcal_j = 0.030;
492 <                }
493 <              }
494 <              float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
495 <              float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
496 <              float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec.Pt();
497 <              float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec.Pt();          
498 <              float comboIso = RIso_i + RIso_j;
499 <              if( info->evtNum == 1038911933 ) {
500 <                float tmpdR = lepvec[i].vec.DrEtaPhi(lepvec[j].vec);
501 <                cout << "i: " << i
502 <                     << "\tdR: " << tmpdR
503 <                     << "\trho: " << rho
504 <                     << "\tRIso_i: " << RIso_i  
505 <                     << "\ttkrel: " << lepvec[i].isoTrk/lepvec[i].vec.Pt()
506 <                     << "\tecalrel: " << lepvec[i].isoEcal/lepvec[i].vec.Pt()
507 <                     << "\tecalrelcor: " << isoEcal_corr_i/lepvec[i].vec.Pt()
508 <                     << "\thcalrel: " << lepvec[i].isoHcal/lepvec[i].vec.Pt()
509 <                     << "\thcalrelcor: " << isoHcal_corr_i/lepvec[i].vec.Pt()
510 <                     << "\tpt_i: " << lepvec[i].vec.Pt()
511 <                     << "\tj: " << j
512 <                     << "\tRIso_j: " << RIso_j  
513 <                     << "\ttkrel: "  << lepvec[j].isoTrk/lepvec[j].vec.Pt()
514 <                     << "\tecalrel: " << lepvec[j].isoEcal/lepvec[j].vec.Pt()
515 <                     << "\tecalrelcor: " << isoEcal_corr_j/lepvec[j].vec.Pt()
516 <                     << "\thcalrel: " << lepvec[j].isoHcal/lepvec[j].vec.Pt()
517 <                     << "\thcalrelcor: " << isoHcal_corr_j/lepvec[j].vec.Pt()
518 <                     << "\tpt_j: " << lepvec[j].vec.Pt()
519 <                     << "\tcombo: " << comboIso
520 <                     << endl;
521 <                cout.flush();
522 <              }
523 <              if( comboIso > 0.35 ) {
524 <                if( ctrl.debug ) cout << "combo failing for indices: " << i << "," << j << endl;
525 <                failiso = true;
526 <                //            break;
527 <              }
438 >                Z2LeptonPlusIndex = j;
439 >                Z2LeptonMinusIndex = i;
440 >              }                  
441              }
529        }
530        if( failiso ) {
531          evtfail |= ( 1<<EVTFAIL_ISOLATION );
532          return evtfail;
533          //h_evtfail->Fill( evtfail, eventweight );
534          //      h_evtfail->Fill( evtfail );
535          //      cout << "evtfail: " << hex << evtfail << dec << endl;
536          //      continue;
537        }
538
539        //***************************************************************
540        // IP significance
541        //***************************************************************
542        bool failip = false;
543        for( int i=0; i<lepvec.size(); i++ ) {
544          if( !(lepvec[i].is4l) ) continue;
545          if( lepvec[i].ip3dSig > 4 ) {
546            failip=true;
547            break;
442            }
443 <        }
444 <        if( failip ) {
445 <          evtfail |= (1<<EVTFAIL_IP );
446 <          return evtfail;
447 <          //h_evtfail->Fill( evtfail, eventweight );
448 <          //      h_evtfail->Fill( evtfail );
449 <          //      cout << "evtfail: " << hex << evtfail << dec << endl;
450 <          //      continue;
451 <        }
452 <
453 <        //***************************************************************
454 <        // remaining kinematic cuts
455 <        //***************************************************************
456 <        double Z2massCut=0;
457 <        if      ( ctrl.kinematics == "loose" ) Z2massCut = 12;
458 <        else if ( ctrl.kinematics == "tight" ) Z2massCut = 20;
459 <        else { cout <<  "error! kinematic tightness not defined!" << endl; assert(0); }
460 <
461 <        if ( Z1Candidate.M() > 120        ||
462 <             Z2Candidate.M() < Z2massCut  ||
463 <             Z2Candidate.M() > 120        ||
464 <             !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec.Pt() > 20.0) ||
465 <             !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec.Pt() > 10.0)
572 <             ) {
573 <          evtfail |= (1<<EVTFAIL_KINEMATICS );
574 <          return evtfail;
575 <          //h_evtfail->Fill( evtfail, eventweight );
576 <          //      h_evtfail->Fill( evtfail );
577 <          //      cout << "evtfail: " << hex << evtfail << dec << endl;
578 <          //      continue;
579 <        }
580 <        
581 <        int channel;
582 <        if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
583 <        if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
584 <        if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) ||
585 <            (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
586 <        
587 <
588 <
589 <        if( passtuple != NULL ) {
590 <          unsigned run   = info->runNum;
591 <          unsigned evt   = info->evtNum;
592 <          unsigned lumi  = info->lumiSec;
593 <          unsigned chan  = channel;
594 <          double   w     = eventweight;
595 <          float mZ1      = Z1Candidate.M() ;
596 <          float mZ2      = Z2Candidate.M() ;
597 <          float m4l      = ZZSystem.M() ;
598 <          float pt4l     = ZZSystem.Pt() ;
599 <          passtuple->SetBranchAddress("run",  &run);
600 <          passtuple->SetBranchAddress("evt",  &evt);
601 <          passtuple->SetBranchAddress("lumi", &lumi);
602 <          passtuple->SetBranchAddress("mZ1",  &mZ1);
603 <          passtuple->SetBranchAddress("mZ2",  &mZ2);
604 <          passtuple->SetBranchAddress("m4l",  &m4l);
605 <          passtuple->SetBranchAddress("pt4l", &pt4l);
606 <          passtuple->SetBranchAddress("w",    &w);
607 <          passtuple->Fill( );
608 <        }
609 <
610 <        if( ctrl.debug ) cout  << "run: " << info->runNum  
611 <              << "\tevt: " << info->evtNum
612 <              << "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
613 <              << "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
614 <              << "\tmZ1: " << Z1Candidate.M()
615 <              << "\tmZ2: " << Z2Candidate.M()
616 <              << "\tm4l: " << ZZSystem.M()
617 <              << "\tevtfail: " << hex << evtfail << dec
618 <              << "\ttrigbits: " << hex << info->triggerBits << dec
619 <          //          << "\ttree: " << inputFiles[q][f]
620 <              << endl;
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  
467 <        return evtfail;
468 <          
469 < }
467 >
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