ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/FSR.cc
Revision: 1.9
Committed: Mon Dec 17 17:12:27 2012 UTC (12 years, 5 months ago) by dkralph
Content type: text/plain
Branch: MAIN
CVS Tags: compiled
Changes since 1.8: +123 -369 lines
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 khahn 1.1 #include <vector>
2     #include "MathUtils.h"
3     #include "FSR.h"
4 khahn 1.3 #include "IsolationSelection.h"
5 khahn 1.1 #include "CommonDefs.h"
6     #include "TLorentzVector.h"
7 dkralph 1.8 #include "Various.h"
8 khahn 1.1
9     extern vector<bool> PFnoPUflag;
10    
11 dkralph 1.9 //----------------------------------------------------------------------------------------
12     void addPhotonToEventData(EventData &ret, TLorentzVector &pvec)
13     //
14     // add a photon to the list of fsr photons in the EventData
15     //
16 dkralph 1.8 {
17     assert(ret.fsrPhotons.size() < 3);
18     SimpleLepton photon;
19     photon.vec.SetPtEtaPhiM( pvec.Pt(), pvec.Eta(), pvec.Phi(), 0);
20    
21     // don't add the same photon twice...
22     if(ret.fsrPhotons.size() > 0)
23     if(dr(ret.fsrPhotons[0], photon) < 0.01) return;
24     if(ret.fsrPhotons.size() > 1)
25     if(dr(ret.fsrPhotons[1], photon) < 0.01) return;
26    
27     if(ret.fsrPhotons.size() == 2) {
28     // only want to save two, so if there's already two in the vector, store the highest two out of the three photons
29     sort( ret.fsrPhotons.begin(), ret.fsrPhotons.end(), SimpleLepton::lep_pt_sort );
30     assert(ret.fsrPhotons[0].vec.Pt() > ret.fsrPhotons[1].vec.Pt());
31     if(ret.fsrPhotons[1].vec.Pt() < photon.vec.Pt()) {
32     ret.fsrPhotons.erase(ret.fsrPhotons.begin() + 1);
33     ret.fsrPhotons.push_back(photon);
34     }
35     } else
36     ret.fsrPhotons.push_back(photon);
37     sort( ret.fsrPhotons.begin(), ret.fsrPhotons.end(), SimpleLepton::lep_pt_sort );
38     }
39 khahn 1.1 //--------------------------------------------------------------------------------------------------
40 dkralph 1.9 // typeI = PF IDed photons.
41     //----------------------------------------------------------------------------------------
42     pair<TLorentzVector,int> findFsrPhoton( ControlFlags & ctrl,
43     EventData &ret,
44     const ChargedParticle *lep, // the non-constant copy of lepvec[i] (or [j])
45     const int lepIndex,
46     vector<SimpleLepton> &lepvec, // really lepvec_i (or _j)
47     const Array<PFCandidate> * pfArr,
48     const Array<Electron> *eleArr,
49     TLorentzVector * Zvec)
50 khahn 1.2 //--------------------------------------------------------------------------------------------------
51 khahn 1.3 {
52 dkralph 1.9 if( lepvec[lepIndex].fsrRecoveryAttempted ) return pair<TLorentzVector,int>(TLorentzVector(0,0,0,0), -1);
53 khahn 1.4
54 dkralph 1.9 vector<PFCandidate> photons;
55     vector<int> photonIndices; // index in the PFCandidate Array (wrong pT for muons FSRs)
56 khahn 1.3 for( int i=0; i<pfArr->GetEntries(); i++ ) {
57 dkralph 1.8 if( !(PFnoPUflag[i])) continue;
58 dkralph 1.9 const PFCandidate *pfOrig = (PFCandidate*)((*pfArr)[i]);
59     PFCandidate pf(*pfOrig);
60    
61     bool isMuonFsr(false);
62     if(abs(pf.PFType())==PFCandidate::eMuon && pf.EECal()>0) {
63     TLorentzVector pfvec;
64     pfvec.SetPtEtaPhiM(pf.EECal()*pf.Pt()/pf.P(), pf.Eta(), pf.Phi(), 0);
65     pf.SetMom(pfvec.Px(), pfvec.Py(), pfvec.Pz(), pfvec.E());
66     pf.SetPFType(PFCandidate::eGamma);
67     isMuonFsr = true;
68     }
69    
70     if(abs(pf.PFType()) != PFCandidate::eGamma) continue;
71     if(pf.Pt() <= 2 ) continue;
72     if(fabs(pf.Eta()) >= 2.4 ) continue;
73    
74     float dR = MathUtils::DeltaR(pf.Phi(),pf.Eta(), lepvec[lepIndex].vec.Phi(), lepvec[lepIndex].vec.Eta());
75    
76     if( ctrl.debug ) {
77     cout << " --> pass pre, " << setprecision(5) << pf.Pt() << " " << dR;
78     if(isMuonFsr) cout << " (muon fsr)";
79     cout << endl;
80     }
81    
82     // veto if close to an electron SC
83     bool flagEleSC = false;
84     for( int j=0; j<lepvec.size(); j++ ) {
85     if( !(abs(lepvec[j].type) == 11 ) ) continue;
86     if( !(lepvec[j].status.looseIDAndPre()) ) continue;
87     double eeta=lepvec[j].vec.Eta(); double ephi=lepvec[j].vec.Phi();
88     float dPhi = fabs(MathUtils::DeltaPhi(pf.Phi(),ephi));
89     float dEta = fabs(pf.Eta()-eeta);
90     float sc_dR = MathUtils::DeltaR(pf.Phi(),pf.Eta(), ephi, eeta);
91     if(ctrl.debug) cout << " sc_dR:" << sc_dR << endl;
92     if( (dPhi<2.&& dEta<0.05) || sc_dR<0.15 ) {
93     flagEleSC = true;
94     break;
95 khahn 1.3 }
96 dkralph 1.9 if( flagEleSC ) break;
97     }
98     if( flagEleSC ) continue;
99     if( ctrl.debug ) cout << " no match SC" << endl;
100 dkralph 1.8
101 dkralph 1.9 // check that input lepton is the closest lepton to this photon
102     bool found_closer_lepton=false;
103     for( int j=0; j<lepvec.size(); j++ ) {
104     if( j == lepIndex ) continue;
105     if( !(lepvec[j].status.looseIDAndPre()) ) continue;
106     float tmp_dR = MathUtils::DeltaR(pf.Phi(),pf.Eta(),
107     lepvec[j].vec.Phi(), lepvec[j].vec.Eta());
108     if( tmp_dR < dR ) {
109     if(ctrl.debug) cout << " found closer lepton (j=" << j << " : " << tmp_dR << " vs " << dR << ") skipping..." << endl;
110 khahn 1.4 found_closer_lepton=true;
111     break;
112 khahn 1.2 }
113     }
114 dkralph 1.9 if( found_closer_lepton ) continue;
115 khahn 1.3
116 dkralph 1.9 // Z mass OK?
117     TLorentzVector pvec;
118     pvec.SetPtEtaPhiM( pf.Pt(), pf.Eta(), pf.Phi(), 0.);
119     float newMass = (pvec + *Zvec).M();
120     if( !( newMass > 4. &&
121     newMass < 100. &&
122     (fabs(newMass-Z_MASS) < fabs(Zvec->M()-Z_MASS))
123     ) ) continue;
124     if( ctrl.debug ) cout << " better mass " << Zvec->M() << " -> " << newMass << endl;
125    
126     // "keep all photons close to one of the 4L leptons ..."
127     bool use(false);
128     if(dR < 0.07) {
129     if( ctrl.debug ) cout << " push loose " << i << endl;
130     use = true;
131     } else if(dR<0.5 && pf.Pt()>4 && isoDr03ForFsr(ctrl, &pf, lep, pfArr, false) < 1) { // "need tighter cuts for other photons ..."
132     if( ctrl.debug ) cout << " push tight " << i << endl;
133     use = true;
134     }
135     if(use) {
136     photons.push_back(pf);
137     photonIndices.push_back(i); // note: will *not* have the right kinematics for muon fsr candidates
138 khahn 1.3 }
139 khahn 1.2 }
140    
141 dkralph 1.9 // choose the best one
142     float highest_pt(-1),smallest_dR(9999999);
143     int highest_pt_index(-1),smallest_dR_index(-1);
144     PFCandidate *hiPtPf(NULL),*smallDrPf(NULL);
145     for( int i=0; i<photons.size(); i++ ) {
146     const PFCandidate pfPhoton = photons[i];
147     if(pfPhoton.Pt() > highest_pt) {
148     highest_pt_index = photonIndices[i]; // index in the PFCandidate Array (wrong pT for muons FSRs)
149     highest_pt = pfPhoton.Pt();
150     hiPtPf = &photons[i];
151     }
152     float this_dR = MathUtils::DeltaR(pfPhoton.Phi(), pfPhoton.Eta(), lep->Phi(), lep->Eta());
153     if(this_dR < smallest_dR) {
154     smallest_dR_index = photonIndices[i]; // index in the PFCandidate Array (wrong pT for muons FSRs)
155     smallest_dR = this_dR;
156     smallDrPf = &photons[i];
157 khahn 1.1 }
158     }
159 khahn 1.3
160 dkralph 1.9 const PFCandidate *pfBest(NULL);
161     int iPfOrig(-1); // index of the final photon in the original PFCandidate array
162     if(highest_pt > 4) {
163     if(ctrl.debug) cout << " using hi pt " << highest_pt_index << endl;
164     pfBest = hiPtPf;
165     iPfOrig = highest_pt_index;
166     } else if(smallest_dR < 99999) {
167     if(ctrl.debug) cout << " using small dR " << smallest_dR_index << endl;
168     pfBest = smallDrPf;
169     iPfOrig = smallest_dR_index;
170 khahn 1.3 } else {
171 dkralph 1.9 return pair<TLorentzVector,int>(TLorentzVector(0,0,0,0), -1);
172 khahn 1.1 }
173    
174 dkralph 1.9 // add to the lepton
175     assert(pfBest);
176     TLorentzVector phvec;
177     phvec.SetPtEtaPhiM( pfBest->Pt(), pfBest->Eta(), pfBest->Phi(), 0);
178     return pair<TLorentzVector,int>(phvec, iPfOrig);
179 khahn 1.1 }