ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/FSR.cc

Comparing UserCode/MitHzz4l/LeptonSelection/src/FSR.cc (file contents):
Revision 1.1 by khahn, Mon May 28 16:40:45 2012 UTC vs.
Revision 1.5 by khahn, Thu Jun 7 10:10:55 2012 UTC

#	Line 1 | Line 1
1		#include <vector>
2		#include "MathUtils.h"
3		#include "FSR.h"
4	+	#include "IsolationSelection.h"
5		#include "CommonDefs.h"
6		#include "TLorentzVector.h"
7
#	Line 12 | Line 13 | extern vector<bool> PFnoPUflag;
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( mithep::Muon * mu,
17	<	const mithep::Array<mithep::PFCandidate> * pfArr )
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
#	Line 23 | Line 32 | bool recover_typeI_Photon( mithep::Muon
32		if( abs(pf->PFType()) == mithep::PFCandidate::eGamma &&
33		pf->Pt() > 2.0 && fabs(pf->Eta()) < 2.4 ) {
34
35	<	float dR = mithep::MathUtils::DeltaR(pf->Phi(),pf->Eta(), mu->Phi(), mu->Eta());
36	<
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	<	if( dR < 0.07 ) photonIndices.push_back(i);
263	<
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	<	if( dR < 0.5 && pf->Pt() > 4. )  photonIndices.push_back(i);
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	<	int index=-1;
280	<	float highest_pt = -1;
281	<	const mithep::PFCandidate * thepf;
282	<	for( int i=0; i<photonIndices.size(); i++ ) {
283	<	const mithep::PFCandidate *pf = (mithep::PFCandidate*)(pfArr->At(photonIndices[i]));
284	<	if( pf->Pt() > highest_pt ) {
285	<	thepf      = pf;
286	<	highest_pt = pf->Pt();
287	<	index      = photonIndices[i];
288	<	}
289	<	}
290	<
291	<	if( thepf != NULL ) {
48	<	// "... remove it from lepton isolation ..."
49	<	PFnoPUflag[index] = 1;
50	<	// add to the muon
51	<	TLorentzVector muvec,phvec,newmuvec;
52	<	muvec.SetPtEtaPhiM( mu->Pt(), mu->Eta(), mu->Phi(), MUON_MASS);
53	<	phvec.SetPtEtaPhiM( thepf->Pt(), thepf->Eta(), thepf->Phi(), 0.);
54	<	newmuvec = muvec+phvec;
55	<	mu->SetPtEtaPhi        (newmuvec.Pt(),
56	<	newmuvec.Eta(),
57	<	newmuvec.Phi());
58	<	return true;
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( mithep::Muon * mu,
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( (pf->TrackerTrk() == mu->TrackerTrk()) && abs(pf->PFType()) == mithep::PFCandidate::eMuon ) {
356	<	foundPF = true;
357	<	thepf = pf;
358	<	break;
359	<	}
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	<	mu->SetPtEtaPhi        (mu->Pt()+phpt,
370	<	mu->Eta(),
371	<	mu->Phi());
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
95	–
382		return false;
383		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines