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.2 by dkralph, Wed Sep 14 12:11:57 2011 UTC vs.
Revision 1.22 by khahn, Mon Apr 30 23:54:34 2012 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines