ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/FSR.cc
Revision: 1.8
Committed: Tue Oct 23 10:39:59 2012 UTC (12 years, 6 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Changes since 1.7: +119 -82 lines
Log Message:
*** empty log message ***

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