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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines