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.1 by khahn, Thu Sep 8 13:33:18 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"
3 #include "HZZCiCElectronSelection.h"
4 #include "RunLumiRangeMap.h"
5
6 mithep::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                               TNtuple * passtuple ) {
20
21  fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, passtuple, NULL );
22
23 };
24
25 unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control
26                               mithep::TEventInfo *info,     // input event inof
27                               TClonesArray *electronArr,    // input electrons
28                               TClonesArray *muonArr,        // input muons
29                               LabVectors  *l ) {
30
31  fails_HZZ4L_selection( ctrl, info, electronArr, muonArr, NULL, l );
32
33 };
34
35
36 unsigned fails_HZZ4L_selection(ControlFlags &ctrl,           // input control
37                               mithep::TEventInfo *info,     // input event inof
38                               TClonesArray *electronArr,    // input electrons
39                               TClonesArray *muonArr,        // input muons
40                               TNtuple * passtuple,
41                               LabVectors * l ) {       // output ntuple
11  
12 + #include "ExternData.h"
13 + #include "SelectionDefs.h"
14 +
15 +
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 <  double eventweight = info->eventweight;
61 >  TRandom3 r;
62 >  
63 >  if( ctrl.debug ) {
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 <    mithep::RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);      
72 >    RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());      
73      if( !(rlrm.HasRunLumi(rl)) )  {
74 <      evtfail |= (1<<EVTFAIL_JSON);
75 <      return evtfail;
74 >      if( ctrl.debug ) cout << "\tfails JSON" << endl;
75 >      ret.status.setStatus(0);
76 >      return ret;
77      }
78    }
79 <  
80 <  
81 <  if( ctrl.debug ) {
82 <    cout << "Run: " << info->runNum
83 <         << "\tEvt: " << info->evtNum
84 <         << "\tLumi: " << info->lumiSec
60 <         << endl;
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 {
74 <    if( !(passHLTMC(info->triggerBits)) ) {
75 <      evtfail |= (1<<EVTFAIL_TRIGGER);
76 <      return evtfail;
77 <    }      
78 <    //    cout << "MC trigger bits: " << hex << info->triggerBits << dec << endl;
79 <  }
80 <
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;
92  
93  //    
94  cout << "\tnMuons: " << muonArr->GetEntries() << endl;
95  //----------------------------------------------------
96  for(Int_t i=0; i<muonArr->GetEntries(); i++) {
97    const mithep::TMuon *mu = (mithep::TMuon*)((*muonArr)[i]);      
98    unsigned muonfail = passMuonSelectionZZ(mu);
99    if( ctrl.debug ) {
100      cout << "muon:: pt: " << mu->pt
101           << "\teta: " << mu->eta
102           << "\tmask: 0x" << hex << muonfail << dec
103           << endl;
104    }
105    if ( !muonfail ) {
106      SimpleLepton tmplep;
107      tmplep.vec.SetPtEtaPhiM(mu->pt,
108                              mu->eta,
109                              mu->phi,
110                              105.658369e-3);
111      tmplep.type    = 13;
112      tmplep.index   = i;
113      tmplep.charge  = mu->q;
114      tmplep.isoTrk  = mu->trkIso03;
115      tmplep.isoEcal = mu->emIso03;
116      tmplep.isoHcal = mu->hadIso03;
117      tmplep.ip3dSig = mu->ip3dSig;
118      tmplep.is4l    = false;
119      tmplep.isEB    = (fabs(mu->eta) < 1.479 ? 1 : 0 );
120      lepvec.push_back(tmplep);
121      if( ctrl.debug ) { cout << "muon passes ... " << endl;}
122    }
123  }
124  
125  //
126  if( ctrl.debug ) { cout << "\tnElectron: " << electronArr->GetEntries() << endl; }
127  //----------------------------------------------------
128  for(Int_t i=0; i<electronArr->GetEntries(); i++) {
129    const mithep::TElectron *ele = (mithep::TElectron*)((*electronArr)[i]);
130
131    Bool_t isMuonOverlap = kFALSE;
132    for (int k=0; k<lepvec.size(); ++k) {
133      TVector3 tmplep;
134      tmplep.SetPtEtaPhi(ele->pt, ele->eta, ele->phi);
135      if ( lepvec[k].type == 13 && lepvec[k].vec.Vect().DrEtaPhi(tmplep) < 0.1 ) {
136        cout << "-----> isMuonOverlap! " << endl;
137        isMuonOverlap = kTRUE;
138        break;
139      }
140    }
116  
117 <    CICStruct tightcuts = getTightCuts();
118 <    unsigned failsCIC = failsCicSelection(ele, tightcuts);
119 <    if( ctrl.debug ){
120 <      cout << "CIC category: " << cicCategory(ele)
121 <           << "\tfailmask : 0x" << hex << failsCIC  << dec
122 <           << "\tscEt: " << ele->scEt
123 <           << "\tscEta: " << ele->scEta
124 <           << endl;
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 <    if ( !failsCIC && !isMuonOverlap ) {
131 <      SimpleLepton tmplep;
132 <      tmplep.vec.SetPtEtaPhiM( ele->pt,
133 <                               ele->eta,
134 <                               ele->phi,
135 <                               0.51099892e-3 );
136 <      tmplep.type    = 11;
137 <      tmplep.index   = i;
138 <      tmplep.charge  = ele->q;
160 <      tmplep.isoTrk  = ele->trkIso03;
161 <      tmplep.isoEcal = ele->emIso03;
162 <      tmplep.isoHcal = ele->hadIso03;
163 <      tmplep.ip3dSig = ele->ip3dSig;
164 <      tmplep.is4l    = false;
165 <      tmplep.isEB    = ele->isEB;
166 <      lepvec.push_back(tmplep);
167 <      if( ctrl.debug ) { cout << "\telectron passes ... " << endl; }
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      }
169  }
170  
171  sort( lepvec.begin(), lepvec.end(), SimpleLepton::lep_pt_sort );
172
173  int nmu=0, nele=0;
174  for( int i=0; i<lepvec.size(); i++ ) {
175    if( abs(lepvec[i].type) == 11 ) nele++;
176    else nmu++;
177  }
178  if( ctrl.debug ) {
179    cout << "postsel nlep: " << lepvec.size()
180         << "\tnmuon: " << nmu
181         << "\tnelectron: " << nele
182         << endl;
183  }
140  
141 <  //******************************************************************************
142 <  //Z1 Selection
143 <  //******************************************************************************
144 <  int Z1LeptonPlusIndex = -1;
145 <  int Z1LeptonMinusIndex = -1;
146 <  double BestZ1Mass = -999;
191 <  if( ctrl.debug ) { cout << "looking for a Z1 ..." << endl; }
192 <  for(int i = 0; i < lepvec.size(); ++i) {
193 <    for(int j = i+1; j < lepvec.size(); ++j) {
194 <      if( ctrl.debug ) { cout << "\tconsidering leptons " << i << " & " << j << endl; }
195 <      if (!(lepvec[i].vec.Pt() > 20.0 || lepvec[j].vec.Pt() > 20.0)) continue;
196 <      if( ctrl.debug ) { cout << "\tat least one is > 20 GeV" << endl; }
197 <      if (!(lepvec[i].vec.Pt() > 10.0 && lepvec[j].vec.Pt() > 10.0)) continue;
198 <      if( ctrl.debug ) { cout << "\tthe other  is > 10 GeV" << endl; }
199 <      if (lepvec[i].charge == lepvec[j].charge) continue;          
200 <      if( ctrl.debug ) { cout << "\tthey're opposite charge" << endl; }
201 <      if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) continue;  
202 <      if( ctrl.debug ) { cout << "\tthey're same flavor" << endl; }
141 >  //    
142 >  if( ctrl.debug ) cout << "\tnMuons: " << muonArr->GetEntries() << endl;
143 >  //----------------------------------------------------
144 >  for(Int_t i=0; i<muonArr->GetEntries(); i++)
145 >    {
146 >      const mithep::Muon *mu = (mithep::Muon*)((*muonArr)[i]);      
147        
148 <      //Make Z1 hypothesis
149 <      TLorentzVector leptonPlus, leptonMinus;
150 <      if ( lepvec[i].charge > 0 ) {
151 <        leptonPlus  = lepvec[i].vec;
152 <        leptonMinus = lepvec[j].vec;
153 <      } else {
154 <        leptonPlus  = lepvec[j].vec;
155 <        leptonMinus = lepvec[i].vec;
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 <      float tmpZ1Mass = (leptonPlus+leptonMinus).M();
166 <      cout << "Z1 selection, tmpZ1Mass: " << tmpZ1Mass << endl;
167 <      if( tmpZ1Mass > 60 ) {
168 <        if (fabs(tmpZ1Mass - 91.1876) < fabs(BestZ1Mass - 91.1876)) {
169 <          BestZ1Mass = tmpZ1Mass;
170 <          cout << "Z1 selection, new BestZ1Mass: " << BestZ1Mass
171 <               << "\tdM: " << fabs(BestZ1Mass - 91.1876)
172 <               << endl;
173 <          if (lepvec[i].charge > 0) {
174 <            Z1LeptonPlusIndex = i;
175 <            Z1LeptonMinusIndex = j;
225 <          } else {
226 <            Z1LeptonPlusIndex = j;
227 <            Z1LeptonMinusIndex = i;
228 <          }
165 >      if ( musel.pass() ) {
166 >
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 <    }
203 <  }
233 <  // stop if no Z1 candidate is found
234 <  if( BestZ1Mass < 0 ) {
235 <    evtfail |= (1<<EVTFAIL_Z1);
236 <    return evtfail;
237 <  }
238 <  cout << "\tgot a Z1 ... run: " << info->runNum << "\tevt: " << info->evtNum  << endl;
239 <  cout << "\tZ1 plusindex: " << Z1LeptonPlusIndex << "\tminusindex: " << Z1LeptonMinusIndex << endl;
240 <  TLorentzVector Z1LeptonPlus  = lepvec[Z1LeptonPlusIndex].vec;
241 <  TLorentzVector Z1LeptonMinus = lepvec[Z1LeptonMinusIndex].vec;
242 <  TLorentzVector Z1Candidate   =  Z1LeptonPlus + Z1LeptonMinus;
243 <  if( l != NULL ) {
244 <    l->vecz1 = Z1Candidate;
245 <    l->vecl1p = Z1LeptonPlus;
246 <    l->vecl1m = Z1LeptonMinus;
247 <  }
202 >      //  }
203 >    }    
204    
205 <  //******************************************************************************
206 <  // Z1 + l
207 <  //******************************************************************************
208 <  if( lepvec.size() < 3 ) {
209 <    evtfail |= (1<<EVTFAIL_Z1_PLUSL);
210 <    return evtfail;
211 <  }
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 <  //******************************************************************************
215 <  // 4l/Z2 Selection
216 <  //******************************************************************************
217 <  Int_t Z2LeptonPlusIndex = -1;
218 <  Int_t Z2LeptonMinusIndex = -1;
219 <  Double_t BestZ2Mass = -1;
220 <  cout << "looking for a Z2 ... out of " << lepvec.size() << " leptons" <<endl;
221 <  for(int i = 0; i < lepvec.size(); ++i) {
222 <    for(int j = i+1; j < lepvec.size(); ++j) {
223 <      //            cout << "i: " << i << "\tj: " << j << endl;
224 <      if (i == Z1LeptonPlusIndex || i == Z1LeptonMinusIndex) {
225 <        //            cout << "\ti matches a Z1 index, skipping ..." << endl;
226 <        continue; //skip Z1 leptons
227 <      }
228 <      if (j == Z1LeptonPlusIndex || j == Z1LeptonMinusIndex) {
229 <        //            cout << "\tj matches a Z1 index, skipping ..." << endl;
230 <              continue; //skip Z1 leptons
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 (lepvec[i].charge == lepvec[j].charge) {
255 <              //              cout << "\ti and j are same sign, skipping ..." << endl;
277 <              continue;         //require opp sign
254 >            if( ctrl.do_escale_up ) {
255 >              pt=scale_smear_electron_Up(pt, ele->IsEB(),  r);
256              }
257 <            if (fabs(lepvec[i].type) != fabs(lepvec[j].type)) {
258 <              //              cout << "\ti and j are not same flavor, skipping ..." << endl;
281 <              continue; //require same flavor
257 >            if( ctrl.do_escale_down ) {
258 >              pt=scale_smear_electron_Down(pt, ele->IsEB(),  r);
259              }
260 +
261              
262 +            tmplep.vec->SetPtEtaPhiM( pt,
263 +                                      ele->Eta(),
264 +                                      ele->Phi(),
265 +                                      ELECTRON_MASS );
266              
267 <            //Make Z2 hypothesis
268 <            TLorentzVector leptonPlus, leptonMinus;
269 <            
270 <            if (lepvec[i].charge > 0 ) {
271 <              leptonPlus  = lepvec[i].vec;
272 <              leptonMinus = lepvec[j].vec;
273 <            } else {
274 <              leptonPlus  = lepvec[j].vec;
275 <              leptonMinus = lepvec[i].vec;
276 <            }
277 <
278 <            TLorentzVector dilepton = leptonPlus+leptonMinus;
279 <            TLorentzVector fourLepton = Z1Candidate + dilepton;
280 <
281 <            cout << "dilepton.M() : " << dilepton.M() << endl;
282 <            cout << "fourLepton.M() : " << fourLepton.M() << endl;
283 <
284 <            if (!(dilepton.M() > 12.0)) continue;
285 <            if (!(fourLepton.M() > 100.0)) continue;
286 <
287 <            
288 <            //for 4e and 4mu, require at least 1 of the other opp sign lepton pairs have mass > 12
289 <            if (fabs(lepvec[i].type) == fabs(lepvec[Z1LeptonPlusIndex].type)) {
290 <              TLorentzVector pair1 = Z1LeptonPlus+leptonMinus;
291 <              TLorentzVector pair2 = Z1LeptonMinus+leptonPlus;
292 <              cout << "pair1: " << pair1.M() << "\tpair2: "<< pair2.M() << endl;
293 <              if (!(pair1.M() > 12 || pair2.M() > 12)) continue;
294 <            }
295 <            
296 <            //Disambiguiation is done by choosing the pair with the largest ptMax and largest ptMin
297 <            if (Z2LeptonPlusIndex < 0) {
298 <              if (lepvec[i].charge > 0) {
299 <                Z2LeptonPlusIndex = i;
300 <                Z2LeptonMinusIndex = j;
301 <              } else {
302 <                Z2LeptonPlusIndex = j;
303 <                Z2LeptonMinusIndex = i;
304 <              }
305 <            } else {
306 <              Double_t BestPairPtMax = lepvec[Z2LeptonPlusIndex].vec.Pt();              
307 <              Double_t BestPairPtMin = lepvec[Z2LeptonMinusIndex].vec.Pt();
308 <              if (lepvec[Z2LeptonMinusIndex].vec.Pt() > BestPairPtMax) {
309 <                BestPairPtMax = lepvec[Z2LeptonMinusIndex].vec.Pt();
310 <                BestPairPtMin = lepvec[Z2LeptonPlusIndex].vec.Pt();
311 <              }
312 <
313 <              Double_t CurrentPairPtMax = lepvec[i].vec.Pt();              
314 <              Double_t CurrentPairPtMin = lepvec[j].vec.Pt();
315 <              if (lepvec[j].vec.Pt() > CurrentPairPtMax) {
316 <                CurrentPairPtMax = lepvec[j].vec.Pt();
317 <                CurrentPairPtMin = lepvec[i].vec.Pt();
318 <              }
319 <
320 <              if (CurrentPairPtMax > BestPairPtMax) {
321 <                if (lepvec[i].charge > 0) {
322 <                  Z2LeptonPlusIndex = i;
323 <                  Z2LeptonMinusIndex = j;
324 <                } else {
325 <                  Z2LeptonPlusIndex = j;
326 <                  Z2LeptonMinusIndex = i;
327 <                }
328 <              } else if (CurrentPairPtMax  == BestPairPtMax) {
329 <                if (CurrentPairPtMin > BestPairPtMin) {
330 <                  if (lepvec[i].charge > 0) {
331 <                    Z2LeptonPlusIndex = i;
332 <                    Z2LeptonMinusIndex = j;
333 <                  } else {
334 <                    Z2LeptonPlusIndex = j;
335 <                    Z2LeptonMinusIndex = i;
336 <                  }                  
337 <                }
338 <              }
339 <            }            
340 <          }
341 <        }
342 <
343 <        // stop if no Z2 candidate is found
344 <        if (Z2LeptonPlusIndex == -1) {
345 <          evtfail |= ( 1<<EVTFAIL_4L );
346 <          return evtfail;
347 <          //      h_evtfail->Fill( evtfail );
348 <          //      cout << "evtfail: " << hex << evtfail << dec << endl;
349 <          //      continue;
368 <        }
369 <        cout << "\tgot a Z2 ..." << endl;
370 <        cout << "\tZ2 plusindex: " << Z2LeptonPlusIndex << "\tminusindex: " << Z2LeptonMinusIndex << endl;
371 <        TLorentzVector  Z2LeptonPlus  = lepvec[Z2LeptonPlusIndex].vec;
372 <        TLorentzVector  Z2LeptonMinus = lepvec[Z2LeptonMinusIndex].vec;
373 <        TLorentzVector  Z2Candidate   = Z2LeptonPlus+Z2LeptonMinus;
374 <        TLorentzVector  ZZSystem      = Z1Candidate + Z2Candidate;
375 <        if( l != NULL ) {
376 <          l->vecz2 = Z2Candidate;
377 <          l->vecl2p = Z2LeptonPlus;
378 <          l->vecl2m = Z2LeptonMinus;
379 <          l->vec4l = ZZSystem;
380 <        }
381 <        lepvec[Z1LeptonPlusIndex].is4l = true;
382 <        lepvec[Z1LeptonMinusIndex].is4l = true;
383 <        lepvec[Z2LeptonPlusIndex].is4l = true;
384 <        lepvec[Z2LeptonMinusIndex].is4l = true;
385 <
386 <        //***************************************************************
387 <        // Isolation
388 <        //***************************************************************
389 <        bool failiso=false;
390 <
391 <        /*
392 <        int i,j;
393 <        i=Z1LeptonPlusIndex;
394 <        j=Z1LeptonMinusIndex;
395 <        float RIso1 = (lepvec[i].isoTrk+lepvec[i].isoEcal+lepvec[i].isoHcal)/lepvec[i].vec.Pt();
396 <        float RIso2 = (lepvec[j].isoTrk+lepvec[j].isoEcal+lepvec[j].isoHcal)/lepvec[j].vec.Pt();
397 <        float comboIso12 = RIso1 + RIso2;
398 <        i=Z2LeptonPlusIndex;
399 <        j=Z2LeptonMinusIndex;
400 <        float RIso3 = (lepvec[i].isoTrk+lepvec[i].isoEcal+lepvec[i].isoHcal)/lepvec[i].vec.Pt();
401 <        float RIso4 = (lepvec[j].isoTrk+lepvec[j].isoEcal+lepvec[j].isoHcal)/lepvec[j].vec.Pt();
402 <        float comboIso34 = RIso3 + RIso4;
403 <        if( comboIso12 > 0.35 || comboIso34 > 0.35 ) {
404 <              failiso = true;
405 <        }
406 <        */
407 <
408 <        float rho = info->PileupEnergyDensity;
409 <        for( int i=0; i<lepvec.size(); i++ ) {
410 <          if( !(lepvec[i].is4l) ) continue;
411 <            float effArea_ecal_i, effArea_hcal_i;
412 <            if( lepvec[i].isEB ) {
413 <              if( lepvec[i].type == 11 ) {
414 <                effArea_ecal_i = 0.101;
415 <                effArea_hcal_i = 0.021;
416 <              } else {
417 <                effArea_ecal_i = 0.074;
418 <                effArea_hcal_i = 0.022;
419 <              }
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 >    
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 >    if( ctrl.debug ) {
301 >      cout << "postsel nlep: " << lepvec.size()
302 >           << "\tnmuon: " << nmu
303 >           << "\tnelectron: " << nele
304 >           << endl;
305 >    }
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 <              if( lepvec[i].type == 11 ) {
352 <                effArea_ecal_i = 0.046;
423 <                effArea_hcal_i = 0.040;
424 <              } else {
425 <                effArea_ecal_i = 0.045;
426 <                effArea_hcal_i = 0.030;
427 <              }
351 >              Z1LeptonPlusIndex = j;
352 >              Z1LeptonMinusIndex = i;
353              }
429            float isoEcal_corr_i = lepvec[i].isoEcal - (effArea_ecal_i*rho);
430            float isoHcal_corr_i = lepvec[i].isoHcal - (effArea_hcal_i*rho);
431            for( int j=i+1; j<lepvec.size(); j++ ) {
432              if( !(lepvec[j].is4l) ) continue;
433              float effArea_ecal_j, effArea_hcal_j;
434              if( lepvec[j].isEB ) {
435                if( lepvec[j].type == 11 ) {
436                  effArea_ecal_j = 0.101;
437                  effArea_hcal_j = 0.021;
438                } else {
439                  effArea_ecal_j = 0.074;
440                  effArea_hcal_j = 0.022;
441                }
442              } else {
443                if( lepvec[j].type == 11 ) {
444                  effArea_ecal_j = 0.046;
445                  effArea_hcal_j = 0.040;
446                } else {
447                  effArea_ecal_j = 0.045;
448                  effArea_hcal_j = 0.030;
449                }
450              }
451              float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
452              float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
453              float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec.Pt();
454              float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec.Pt();          
455              float comboIso = RIso_i + RIso_j;
456              if( info->evtNum == 1038911933 ) {
457                float tmpdR = lepvec[i].vec.DrEtaPhi(lepvec[j].vec);
458                cout << "i: " << i
459                     << "\tdR: " << tmpdR
460                     << "\trho: " << rho
461                     << "\tRIso_i: " << RIso_i  
462                     << "\ttkrel: " << lepvec[i].isoTrk/lepvec[i].vec.Pt()
463                     << "\tecalrel: " << lepvec[i].isoEcal/lepvec[i].vec.Pt()
464                     << "\tecalrelcor: " << isoEcal_corr_i/lepvec[i].vec.Pt()
465                     << "\thcalrel: " << lepvec[i].isoHcal/lepvec[i].vec.Pt()
466                     << "\thcalrelcor: " << isoHcal_corr_i/lepvec[i].vec.Pt()
467                     << "\tpt_i: " << lepvec[i].vec.Pt()
468                     << "\tj: " << j
469                     << "\tRIso_j: " << RIso_j  
470                     << "\ttkrel: "  << lepvec[j].isoTrk/lepvec[j].vec.Pt()
471                     << "\tecalrel: " << lepvec[j].isoEcal/lepvec[j].vec.Pt()
472                     << "\tecalrelcor: " << isoEcal_corr_j/lepvec[j].vec.Pt()
473                     << "\thcalrel: " << lepvec[j].isoHcal/lepvec[j].vec.Pt()
474                     << "\thcalrelcor: " << isoHcal_corr_j/lepvec[j].vec.Pt()
475                     << "\tpt_j: " << lepvec[j].vec.Pt()
476                     << "\tcombo: " << comboIso
477                     << endl;
478                cout.flush();
479              }
480              if( comboIso > 0.35 ) {
481                cout << "combo failing for indices: " << i << "," << j << endl;
482                failiso = true;
483                //            break;
484              }
485            }
486        }
487        if( failiso ) {
488          evtfail |= ( 1<<EVTFAIL_ISOLATION );
489          return evtfail;
490          //h_evtfail->Fill( evtfail, eventweight );
491          //      h_evtfail->Fill( evtfail );
492          //      cout << "evtfail: " << hex << evtfail << dec << endl;
493          //      continue;
494        }
495
496        //***************************************************************
497        // IP significance
498        //***************************************************************
499        bool failip = false;
500        for( int i=0; i<lepvec.size(); i++ ) {
501          if( !(lepvec[i].is4l) ) continue;
502          if( lepvec[i].ip3dSig > 4 ) {
503            failip=true;
504            break;
354            }
355          }
356 <        if( failip ) {
357 <          evtfail |= (1<<EVTFAIL_IP );
358 <          return evtfail;
359 <          //h_evtfail->Fill( evtfail, eventweight );
360 <          //      h_evtfail->Fill( evtfail );
361 <          //      cout << "evtfail: " << hex << evtfail << dec << endl;
362 <          //      continue;
363 <        }
364 <
365 <        //***************************************************************
366 <        // remaining kinematic cuts
367 <        //***************************************************************
368 <        if ( Z1Candidate.M() > 120 ||
369 <             Z2Candidate.M() < 20  ||
370 <             Z2Candidate.M() > 120 ||
371 <             !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 20.0 || lepvec[Z1LeptonMinusIndex].vec.Pt() > 20.0) ||
372 <             !(lepvec[Z1LeptonPlusIndex].vec.Pt() > 10.0 && lepvec[Z1LeptonMinusIndex].vec.Pt() > 10.0)
373 <             ) {
374 <          evtfail |= (1<<EVTFAIL_KINEMATICS );
375 <          return evtfail;
376 <          //h_evtfail->Fill( evtfail, eventweight );
377 <          //      h_evtfail->Fill( evtfail );
378 <          //      cout << "evtfail: " << hex << evtfail << dec << endl;
379 <          //      continue;
380 <        }
381 <        
382 <        int channel;
383 <        if( lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 11 ) channel=0;
384 <        if( lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 13 ) channel=1;
385 <        if( (lepvec[Z1LeptonMinusIndex].type == 11 && lepvec[Z2LeptonMinusIndex].type == 13) ||
386 <            (lepvec[Z1LeptonMinusIndex].type == 13 && lepvec[Z2LeptonMinusIndex].type == 11)) channel=2;
356 >      }
357 >    }
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 >      if( ctrl.eleSeleScheme == "mediumloose" &&
397 >          !(lepvec[i].isTight) ) {
398 >        if( ctrl.debug)  cout << "it's not tight, skipping ... " << endl;
399 >        continue;
400 >      }
401 >      
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 <
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( passtuple != NULL ) {
451 <          passtuple->Fill( info->runNum,
452 <                           info->evtNum,
453 <                           info->lumiSec,
454 <                           channel,
455 <                           Z1Candidate.M(),
456 <                           Z2Candidate.M(),
457 <                           ZZSystem.M(),
458 <                           ZZSystem.Pt(),
459 <                           eventweight);
460 <        }
461 <
462 <        cout  << "run: " << info->runNum  
463 <              << "\tevt: " << info->evtNum
464 <              << "\tZ1channel: " << lepvec[Z1LeptonMinusIndex].type
465 <              << "\tZ2channel: " << lepvec[Z2LeptonMinusIndex].type
466 <              << "\tmZ1: " << Z1Candidate.M()
467 <              << "\tmZ2: " << Z2Candidate.M()
468 <              << "\tm4l: " << ZZSystem.M()
469 <              << "\tevtfail: " << hex << evtfail << dec
470 <              << "\ttrigbits: " << hex << info->triggerBits << dec
471 <          //          << "\ttree: " << inputFiles[q][f]
472 <              << endl;
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 >            Z2LeptonPlusIndex = i;
466 >            Z2LeptonMinusIndex = j;
467 >          } else {
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 >          }
478 >        
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 >              Z2LeptonPlusIndex = j;
492 >              Z2LeptonMinusIndex = i;
493 >            }
494 >          } else if (CurrentPairPtMax  == BestPairPtMax) {
495 >            if (CurrentPairPtMin > BestPairPtMin) {
496 >              if (lepvec[i].charge > 0) {
497 >                Z2LeptonPlusIndex = i;
498 >                Z2LeptonMinusIndex = j;
499 >              } else {
500 >                Z2LeptonPlusIndex = j;
501 >                Z2LeptonMinusIndex = i;
502 >              }                  
503 >            }
504 >          }
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  
529 <        return evtfail;
530 <          
531 < }
529 >
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