ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/FSR.cc
Revision: 1.4
Committed: Tue Jun 5 19:35:47 2012 UTC (12 years, 11 months ago) by khahn
Content type: text/plain
Branch: MAIN
CVS Tags: synced_FSR
Changes since 1.3: +74 -60 lines
Log Message:
final sync for FSR

File Contents

# Content
1 #include <vector>
2 #include "MathUtils.h"
3 #include "FSR.h"
4 #include "IsolationSelection.h"
5 #include "CommonDefs.h"
6 #include "TLorentzVector.h"
7
8 using namespace std;
9
10 extern vector<bool> PFnoPUflag;
11
12 //--------------------------------------------------------------------------------------------------
13 // typeI = PF IDed photons. NB : repurpose PFnoPUflag, flip for recovered photons
14 // so that they are skipped in the isolation calculation
15 //--------------------------------------------------------------------------------------------------
16 bool recover_typeI_Photon( ControlFlags & ctrl,
17 mithep::Electron * el,
18 const int electronIndex,
19 vector<SimpleLepton> &lepvec,
20 const mithep::Array<mithep::PFCandidate> * pfArr,
21 const mithep::Array<mithep::Electron> *eleArr,
22 TLorentzVector * Zvec,
23 vector<const mithep::PFCandidate*> &photonsToVeto )
24 //--------------------------------------------------------------------------------------------------
25 {
26 if( lepvec[electronIndex].fsrRecoveryAttempted ) return false;
27
28 vector<int> photonIndices;
29 for( int i=0; i<pfArr->GetEntries(); i++ ) {
30 if( !(PFnoPUflag[i])) continue; // my PF no PU hack
31 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfArr)[i]);
32 if( abs(pf->PFType()) == mithep::PFCandidate::eGamma &&
33 pf->Pt() > 2.0 && fabs(pf->Eta()) < 2.4 ) {
34
35 if( ctrl.debug ) std::cerr << "FSR :: pass preselection ... pt: "<< pf->Pt() << std::endl;
36 // float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), el->Phi(), el->Eta());
37 float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(),
38 lepvec[electronIndex].vec.Phi(), lepvec[electronIndex].vec.Eta());
39 if( ctrl.debug ) std::cerr << "FSR :: dR = " << dR << std::endl;
40
41 //
42 // veto if close to an electron SC
43 //
44 bool flagEleSC = false;
45 for( int j=0; j<lepvec.size(); j++ ) {
46 if( !(abs(lepvec[j].type) == 11 ) ) continue;
47 if( !(lepvec[j].status.looseIDAndPre()) ) continue;
48 double eeta=lepvec[j].vec.Eta(); double ephi=lepvec[j].vec.Phi();
49 float dPhi = fabs(mithep::MathUtils::DeltaPhi(pf->Phi(),ephi));
50 float dEta = fabs(pf->Eta()-eeta);
51 float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), ephi, eeta);
52 if(ctrl.debug) cout << "FSR :: comparing to ele, dPhi: " << dPhi
53 << "\tdEta: " << dEta
54 << "\tetaPH: " << pf->Eta()
55 << "\tetaELH: " << eeta
56 << "\tdR:" << dR << endl;
57 if( (dPhi<2.&& dEta<0.05) || dR<0.15 ) {
58 flagEleSC = true;
59 break;
60 }
61 if( flagEleSC ) break;
62 }
63 if( flagEleSC ) continue;
64 if( ctrl.debug ) std::cerr << "FSR :: not matched to an ele SC ... " << std::endl;
65
66
67 //
68 // check that input electron is the closest lepton to this photon
69 //
70 bool found_closer_lepton=false;
71 for( int j=0; j<lepvec.size(); j++ ) {
72 if( j == electronIndex ) continue;
73 if( !(lepvec[j].status.looseIDAndPre()) ) continue;
74 float tmp_dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(),
75 lepvec[j].vec.Phi(), lepvec[j].vec.Eta());
76 if( tmp_dR < dR ) {
77 if(ctrl.debug) cout << "FSR :: found closer lepton (j="<<j<<" : "
78 <<tmp_dR<<" vs "<<dR<<") skipping..." << endl;
79 found_closer_lepton=true;
80 break;
81 }
82 }
83 if( found_closer_lepton ) continue;
84
85
86 //
87 // Z mass OK?
88 //
89 TLorentzVector pvec;
90 pvec.SetPtEtaPhiM( pf->Pt(), pf->Eta(), pf->Phi(), 0.);
91 float newMass = (pvec + *Zvec).M();
92 if( !( newMass > 4. &&
93 newMass < 100. &&
94 (fabs(newMass-Z_MASS) < fabs(Zvec->M()-Z_MASS))
95 ) ) continue;
96 if( ctrl.debug ) std::cerr << "FSR :: improved Zmass ... " <<
97 Zvec->M() << " -> " << newMass << std::endl;
98
99
100 //
101 // "keep all photons close to one of the 4L electrons ..."
102 //
103 if( dR < 0.07 ) {
104 if( ctrl.debug ) std::cerr << "FSR :: dR < 0.07, pushing ... " << std::endl;
105 photonIndices.push_back(i);
106 }
107
108 //
109 // "need tighter cuts for other photons ..."
110 //
111 if( dR < 0.5 && pf->Pt() > 4. && dbetaCorrectedIsoDr03(ctrl, pf, el, pfArr) < 1.0) {
112 // if( dR < 0.5 && pf->Pt() > 4. && nonCorrectedIsoDr03(ctrl, pf, el, pfArr) < 1.0) {
113 if( ctrl.debug ) std::cerr << "FSR :: tighter cuts, pushing ... " << std::endl;
114 photonIndices.push_back(i);
115 }
116 }
117 }
118
119 float highest_pt = -1; int highest_pt_index=-1;
120 float smallest_dR = 999.; int smallest_dR_index=-1;
121 for( int i=0; i<photonIndices.size(); i++ ) {
122 const mithep::PFCandidate *pf = (mithep::PFCandidate*)(pfArr->At(photonIndices[i]));
123 float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), el->Phi(), el->Eta());
124 if( pf->Pt() > highest_pt ) {
125 highest_pt_index = photonIndices[i];
126 highest_pt = pf->Pt();
127 }
128 if( dR < smallest_dR ) {
129 smallest_dR_index = photonIndices[i];
130 smallest_dR = dR;
131 }
132 }
133
134 const mithep::PFCandidate * thepf;
135 if( highest_pt > 4. ) {
136 thepf = (const mithep::PFCandidate*)(pfArr->At(highest_pt_index));
137 // "... remove it from lepton isolation ..."
138 // PFnoPUflag[highest_pt_index] = 0;
139 // TMP, commented flip above for FSR study
140 // gammaMatches[highest_pt_index].push_back(lepvec[electronIndex].index);
141 photonsToVeto.push_back(thepf);
142 } else if( smallest_dR != 999. ) {
143 thepf = (const mithep::PFCandidate*)(pfArr->At(smallest_dR_index));
144 // "... remove it from lepton isolation ..."
145 // PFnoPUflag[smallest_dR_index] = 0;
146 // TMP, commented flip above for FSR study
147 // gammaMatches[smallest_dR_index].push_back(lepvec[electronIndex].index);
148 photonsToVeto.push_back(thepf);
149 } else {
150 return false;
151 }
152
153 if( thepf != NULL ) {
154 // add to the electron
155 TLorentzVector elvec,phvec,newelvec;
156 elvec.SetPtEtaPhiM( el->Pt(), el->Eta(), el->Phi(), ELECTRON_MASS);
157 phvec.SetPtEtaPhiM( thepf->Pt(), thepf->Eta(), thepf->Phi(), 0.);
158 newelvec = elvec+phvec;
159 // don't update the electron object, just simplelepton
160 // el->SetPtEtaPhi (newelvec.Pt(),
161 // newelvec.Eta(),
162 // newelvec.Phi());
163 lepvec[electronIndex].vec += phvec;
164 lepvec[electronIndex].fsrRecoveryAttempted = true;
165 return true;
166 }
167 return false;
168 }
169
170
171
172 //--------------------------------------------------------------------------------------------------
173 // typeI = PF IDed photons. NB : repurpose PFnoPUflag, flip for recovered photons
174 // so that they are skipped in the isolation calculation
175 //--------------------------------------------------------------------------------------------------
176 bool recover_typeI_Photon( ControlFlags & ctrl,
177 mithep::Muon * mu,
178 const int muonIndex,
179 vector<SimpleLepton> &lepvec,
180 const mithep::Array<mithep::PFCandidate> * pfArr,
181 const mithep::Array<mithep::Electron> *eleArr,
182 TLorentzVector * Zvec,
183 vector<const mithep::PFCandidate*> &photonsToVeto )
184 //--------------------------------------------------------------------------------------------------
185 {
186 if( lepvec[muonIndex].fsrRecoveryAttempted ) return false;
187
188 vector<int> photonIndices;
189 for( int i=0; i<pfArr->GetEntries(); i++ ) {
190 if( !(PFnoPUflag[i])) continue; // my PF no PU hack
191 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfArr)[i]);
192 if( abs(pf->PFType()) == mithep::PFCandidate::eGamma &&
193 pf->Pt() > 2.0 && fabs(pf->Eta()) < 2.4 ) {
194
195 if( ctrl.debug ) std::cerr << "FSR :: pass preselection ... pt: "<< pf->Pt() << std::endl;
196 // float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), mu->Phi(), mu->Eta());
197 float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(),
198 lepvec[muonIndex].vec.Phi(), lepvec[muonIndex].vec.Eta());
199 if( ctrl.debug ) std::cerr << "FSR :: dR = " << dR << std::endl;
200
201 //
202 // veto if close to an electron SC
203 //
204 bool flagEleSC = false;
205 for( int j=0; j<lepvec.size(); j++ ) {
206 if( !(abs(lepvec[j].type) == 11 ) ) continue;
207 if( !(lepvec[j].status.looseIDAndPre()) ) continue;
208 double eeta=lepvec[j].vec.Eta(); double ephi=lepvec[j].vec.Phi();
209 float dPhi = fabs(mithep::MathUtils::DeltaPhi(pf->Phi(),ephi));
210 float dEta = fabs(pf->Eta()-eeta);
211 float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), ephi, eeta);
212 if(ctrl.debug) cout << "FSR :: comparing to ele, dPhi: " << dPhi
213 << "\tdEta: " << dEta
214 << "\tetaPH: " << pf->Eta()
215 << "\tetaELH: " << eeta
216 << "\tdR:" << dR << endl;
217 if( (dPhi<2.&& dEta<0.05) || dR<0.15 ) {
218 flagEleSC = true;
219 break;
220 }
221 if( flagEleSC ) break;
222 }
223 if( flagEleSC ) continue;
224 if( ctrl.debug ) std::cerr << "FSR :: not matched to an ele SC ... " << std::endl;
225
226
227 //
228 // check that input muon is the closest lepton to this photon
229 //
230 bool found_closer_lepton=false;
231 for( int j=0; j<lepvec.size(); j++ ) {
232 if( j == muonIndex ) continue;
233 if( !(lepvec[j].status.looseIDAndPre()) ) continue;
234 float tmp_dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(),
235 lepvec[j].vec.Phi(), lepvec[j].vec.Eta());
236 if( tmp_dR < dR ) {
237 if(ctrl.debug) cout << "FSR :: found closer lepton (j="<<j<<" : "
238 <<tmp_dR<<" vs "<<dR<<") skipping..." << endl;
239 found_closer_lepton=true;
240 break;
241 }
242 }
243 if( found_closer_lepton ) continue;
244
245
246
247 //
248 // Z mass OK?
249 //
250 TLorentzVector pvec;
251 pvec.SetPtEtaPhiM( pf->Pt(), pf->Eta(), pf->Phi(), 0.);
252 float newMass = (pvec + *Zvec).M();
253 if( !( newMass > 4. &&
254 newMass < 100. &&
255 (fabs(newMass-Z_MASS) < fabs(Zvec->M()-Z_MASS))
256 ) ) continue;
257 if( ctrl.debug ) std::cerr << "FSR :: improved Zmass ... " <<
258 Zvec->M() << " -> " << newMass << std::endl;
259
260 //
261 // "keep all photons close to one of the 4L muons ..."
262 //
263 if( dR < 0.07 ) {
264 if( ctrl.debug ) std::cerr << "FSR :: dR < 0.07, pushing ... " << std::endl;
265 photonIndices.push_back(i);
266 }
267
268 //
269 // "need tighter cuts for other photons ..."
270 //
271 if( ctrl.debug ) std::cerr << "FSR :: pass tighter?, pT: " << pf->Pt() << std::endl;
272 if( dR < 0.5 && pf->Pt() > 4. && dbetaCorrectedIsoDr03(ctrl, pf, mu, pfArr) < 1.0) {
273 // if( dR < 0.5 && pf->Pt() > 4. && nonCorrectedIsoDr03(ctrl, pf, mu, pfArr) < 1.0) {
274 if( ctrl.debug ) std::cerr << "FSR :: tighter cuts, pushing index= " << i << std::endl;
275 photonIndices.push_back(i);
276 }
277 }
278 }
279
280 float highest_pt = -1; int highest_pt_index=-1;
281 float smallest_dR = 999.; int smallest_dR_index=-1;
282 for( int i=0; i<photonIndices.size(); i++ ) {
283 const mithep::PFCandidate *pf = (mithep::PFCandidate*)(pfArr->At(photonIndices[i]));
284 float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), mu->Phi(), mu->Eta());
285 if( pf->Pt() > highest_pt ) {
286 highest_pt_index = photonIndices[i];
287 highest_pt = pf->Pt();
288 }
289 if( dR < smallest_dR ) {
290 smallest_dR_index = photonIndices[i];
291 smallest_dR = dR;
292 }
293 }
294
295 const mithep::PFCandidate * thepf;
296 if( highest_pt > 4. ) {
297 if(ctrl.debug) std::cerr << "FSR :: taking highest pt gamma, index = " << highest_pt_index << endl;
298 thepf = (const mithep::PFCandidate*)(pfArr->At(highest_pt_index));
299 // "... remove it from lepton isolation ..."
300 // PFnoPUflag[highest_pt_index] = 0;
301 // TMP, commented flip above for FSR study
302 // gammaMatches[highest_pt_index].push_back(lepvec[muonIndex].index);
303 photonsToVeto.push_back(thepf);
304 } else if( smallest_dR != 999. ) {
305 if(ctrl.debug) std::cerr << "FSR :: taking smallest dR gamma, index = " << highest_pt_index << endl;
306 thepf = (const mithep::PFCandidate*)(pfArr->At(smallest_dR_index));
307 // "... remove it from lepton isolation ..."
308 // PFnoPUflag[smallest_dR_index] = 0;
309 // TMP, commented flip above for FSR study
310 //gammaMatches[smallest_dR_index].push_back(lepvec[muonIndex].index);
311 photonsToVeto.push_back(thepf);
312 } else {
313 return false;
314 }
315
316 TLorentzVector pvec;
317 if( thepf != NULL ) {
318 // add to the muon
319 if( ctrl.debug ) cerr << "FSR :: before return, oldpT=" << mu->Pt() << endl;
320 TLorentzVector muvec,phvec,newmuvec;
321 muvec.SetPtEtaPhiM( mu->Pt(), mu->Eta(), mu->Phi(), MUON_MASS);
322 phvec.SetPtEtaPhiM( thepf->Pt(), thepf->Eta(), thepf->Phi(), 0.);
323 pvec = phvec;
324 newmuvec = muvec+phvec;
325 // don't update the muon object, just simplelepton
326 // mu->SetPtEtaPhi (newmuvec.Pt(),
327 // newmuvec.Eta(),
328 // newmuvec.Phi());
329 lepvec[muonIndex].vec += phvec;
330 lepvec[muonIndex].fsrRecoveryAttempted = true;
331 return true;
332 }
333 return false;
334 }
335
336
337
338 //--------------------------------------------------------------------------------------------------
339 // typeII = "PFClusters linked to muons"
340 //--------------------------------------------------------------------------------------------------
341 bool recover_typeII_Photon( ControlFlags & ctrl,
342 mithep::Muon * mu,
343 const int muonIndex,
344 vector<SimpleLepton> &lepvec,
345 const mithep::Array<mithep::PFCandidate> * pfArr )
346 //--------------------------------------------------------------------------------------------------
347 {
348 if( lepvec[muonIndex].fsrRecoveryAttempted ) return false;
349
350 bool foundPF=false;
351 const mithep::PFCandidate * thepf;
352 for( int i=0; i<pfArr->GetEntries(); i++ ) {
353 if( !(PFnoPUflag[i]) ) continue; // my PF no PU hack
354 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfArr)[i]);
355 if( abs(pf->PFType()) == mithep::PFCandidate::eMuon
356 && (mu->TrackerTrk()==pf->TrackerTrk()))
357 {
358 if(ctrl.debug) cout << "FSR :: t2, found pf muon, pt " << mu->Pt() << endl;
359 foundPF = true;
360 thepf = pf;
361 break;
362 }
363 }
364
365 if( foundPF ) {
366 double sintet = thepf->Pt()/thepf->E();
367 double phpt = thepf->EECal() * sintet;
368 if ( thepf->EECal() >= 2.0 && phpt >= 2.0 ) {
369 // don't update the muon object, just simplelepton
370 // mu->SetPtEtaPhi (mu->Pt()+phpt,
371 // mu->Eta(),
372 // mu->Phi());
373 TLorentzVector pvec;
374 pvec.SetPtEtaPhiM( phpt,mu->Eta(),mu->Phi(),0.);
375 lepvec[muonIndex].vec += pvec;
376 if(ctrl.debug) cout << "FSR :: t2, new pt " << lepvec[muonIndex].vec.Pt() << endl;
377 lepvec[muonIndex].fsrRecoveryAttempted = true;
378 return true;
379 }
380 }
381
382 return false;
383 }