ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/ElectronFakeRateMod.cc
Revision: 1.7
Committed: Wed Jul 18 12:35:27 2012 UTC (12 years, 10 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.6: +1 -74 lines
Log Message:
added electron energy stuff and miscellaneous other stuff

File Contents

# User Rev Content
1 anlevin 1.1 #include "ElectronFakeRateMod.h"
2 anlevin 1.7 #include "SelectionFuncs.h"
3 anlevin 1.1
4    
5     //needed because of externs
6     vector<bool> PFnoPUflag;
7     Float_t computePFMuonIso(const mithep::Muon *mu,
8     const mithep::Vertex * vtx,
9     const mithep::Array<mithep::PFCandidate> * pfCandidates,
10     const Double_t dRMax);
11    
12     Bool_t passElectronDenominatorCuts(const mithep::Electron *ele, const mithep::Vertex * vtx, Int_t DenominatorType);
13    
14     string IntToString(int i) {
15     char temp[100];
16     sprintf(temp, "%d", i);
17     string str = temp;
18     return str;
19     }
20    
21     class TTree;
22     class TFile;
23     class TString;
24    
25     const Float_t g_electron_mass = 0.51099892e-3;
26     const Float_t g_muon_mass = 105.658369e-3;
27    
28    
29     void mithep::ElectronFakeRateMod::SlaveBegin()
30     {
31    
32    
33     denominatorType.push_back(1);
34     denominatorType.push_back(2);
35     denominatorType.push_back(3);
36     denominatorType.push_back(4);
37    
38     initElectronIDMVAV1();
39    
40     ReqBranch("PFMet", fPFMet);
41    
42     ReqBranch(mithep::Names::gkMvfConversionBrn, fConversions);
43    
44     fMuonName = mithep::Names::gkMuonBrn;
45     ReqBranch(fMuonName, fMuons);
46    
47     fElectronName = mithep::Names::gkElectronBrn;
48     ReqBranch(fElectronName, fElectrons);
49    
50     fPFCandidateName = mithep::Names::gkPFCandidatesBrn;
51     ReqBranch(fPFCandidateName, fPFCandidates);
52    
53     fPrimVtxName = mithep::Names::gkPVBrn;
54     ReqBranch(fPrimVtxName, fPrimVerts);
55    
56     fPUEnergyDensityName = mithep::Names::gkPileupEnergyDensityBrn;
57     ReqBranch(fPUEnergyDensityName, fPileupEnergyDensity);
58    
59     fTrigMaskName = mithep::Names::gkHltBitBrn;
60     ReqBranch(fTrigMaskName, fTrigMask);
61    
62     fPileupName = mithep::Names::gkPileupInfoBrn;
63     ReqBranch(fPileupName, fPileup);
64    
65    
66     fOutputFile = new TFile(fOutputName, "RECREATE");
67     fOutputTrees = vector<TTree*>(denominatorType.size(),0);
68     for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
69     fOutputTrees[denominatorTypeIndex] = new TTree(TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])),TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])));
70     }
71    
72     rlrm.AddJSONFile("/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions12/8TeV/Prompt/Cert_190456-193557_8TeV_PromptReco_Collisions12_JSON.txt");
73    
74     fake_electron_pt = vector<Float_t>(denominatorType.size(),0);
75     fake_electron_eta = vector<Float_t>(denominatorType.size(),0);
76     fake_electron_phi = vector<Float_t>(denominatorType.size(),0);
77     fake_electron_pass = vector<UInt_t>(denominatorType.size(),0);
78    
79    
80     for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
81     fOutputTrees[denominatorTypeIndex]->Branch("nvertices",&nvertices,"nvertices/i");
82     fOutputTrees[denominatorTypeIndex]->Branch("evt_num",&evt_num,"evt_num/i");
83     fOutputTrees[denominatorTypeIndex]->Branch("lumi_sec",&lumi_sec,"lumi_sec/i");
84     fOutputTrees[denominatorTypeIndex]->Branch("run_num",&run_num,"run_num/i");
85    
86     fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_pt"),&fake_electron_pt[denominatorTypeIndex],TString("fake_electron_pt/F"));
87     fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_eta"),&fake_electron_eta[denominatorTypeIndex],TString("fake_electron_eta/F"));
88     fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_phi"),&fake_electron_phi[denominatorTypeIndex],TString("fake_electron_phi/F"));
89     fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_electron_pass"),&fake_electron_pass[denominatorTypeIndex],TString("fake_electron_pass/i"));
90     }
91    
92     }
93    
94     void mithep::ElectronFakeRateMod::SlaveTerminate()
95     {
96     fOutputFile->cd();
97     for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
98     fOutputTrees[denominatorTypeIndex]->Print();
99    
100     }
101     fOutputFile->Write();
102     fOutputFile->Close();
103    
104     }
105    
106    
107     void mithep::ElectronFakeRateMod::Process()
108     {
109    
110     vector<UInt_t> nElectrons(denominatorType.size(),0);
111    
112     gDebugMask = mithep::Debug::kAnalysis; // debug message category
113     gDebugLevel = 1; // higher level allows more messages to print
114    
115     LoadBranch(fMuonName);
116     LoadBranch(fElectronName);
117     LoadBranch(fPFCandidateName);
118     LoadBranch(fPrimVtxName);
119     LoadBranch(fPUEnergyDensityName);
120     LoadBranch(fTrigMaskName);
121     LoadBranch("PFMet");
122     LoadBranch(mithep::Names::gkMvfConversionBrn);
123    
124     if(fPFMet->At(0)->Et() >25) return;
125    
126     if(store_npu)
127     LoadBranch(fPileupName);
128    
129     mithep::RunLumiRangeMap::RunLumiPairType rl(GetEventHeader()->RunNum(), GetEventHeader()->LumiSec());
130     if(useJSON && !rlrm.HasRunLumi(rl)) return;
131    
132     fillTriggerBits( GetHLTTable(), fTrigMask, fTriggerBits );
133    
134 anlevin 1.6 PFnoPUflag.clear();
135 anlevin 1.1 makePFnoPUArray(fPFCandidates, PFnoPUflag, fPrimVerts );
136    
137     vector<const Muon*> goodMuons;
138     vector<const Electron*> goodElectrons;
139    
140     for(Int_t i=0; i<fMuons->GetEntries(); i++) {
141 anlevin 1.3 const Muon *mu = fMuons->At(i);
142 anlevin 1.4 if( muonPOG2012CutBasedIDTight(mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity) )
143 anlevin 1.1 goodMuons.push_back(mu);
144     }
145     for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
146     const Electron *ele = fElectrons->At(i);
147     Bool_t fails_some_denominator = kFALSE;
148     for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
149     if (passElectronDenominatorCuts(ele, fPrimVerts->At(0), denominatorType[denominatorTypeIndex]))
150     nElectrons[denominatorTypeIndex]++;
151     else
152     fails_some_denominator = kTRUE;
153     }
154    
155 dkralph 1.5 if(!fails_some_denominator && electronPOG2012CutBasedIDMedium(ele, fPrimVerts->At(0), fPFCandidates, fConversions, fPileupEnergyDensity))
156 anlevin 1.1 goodElectrons.push_back(ele);
157     }
158    
159     Int_t NZCandidates = 0;
160     const mithep::Muon* ZMuon1 = 0;
161     const mithep::Muon* ZMuon2 = 0;
162     const mithep::Electron *ZEle1 = 0;
163     const mithep::Electron *ZEle2 = 0;
164     Double_t ZPt = 0;
165     Double_t ZMass = 0;
166     for(UInt_t i=0; i<goodMuons.size(); i++) {
167     FourVectorM mu1;
168     mu1.SetCoordinates(goodMuons[i]->Pt(), goodMuons[i]->Eta(), goodMuons[i]->Phi(), g_muon_mass );
169     for(UInt_t j=i+1; j<goodMuons.size(); j++) {
170     FourVectorM mu2;
171     mu2.SetCoordinates(goodMuons[j]->Pt(), goodMuons[j]->Eta(), goodMuons[j]->Phi(), g_muon_mass );
172     FourVectorM dilepton = mu1+mu2;
173     if (dilepton.M() > massLo && dilepton.M() < massHi) {
174     ZMuon1 = goodMuons[i];
175     ZMuon2 = goodMuons[j];
176     ZPt = dilepton.Pt();
177     ZMass = dilepton.M();
178     NZCandidates++;
179     }
180     }
181     }
182    
183     for(UInt_t i=0; i<goodElectrons.size(); i++) {
184     FourVectorM ele1;
185     ele1.SetCoordinates(goodElectrons[i]->Pt(), goodElectrons[i]->Eta(), goodElectrons[i]->Phi(),g_electron_mass );
186     for(UInt_t j=i+1; j<goodElectrons.size(); j++) {
187     FourVectorM ele2;
188     ele2.SetCoordinates(goodElectrons[j]->Pt(), goodElectrons[j]->Eta(), goodElectrons[j]->Phi(), g_electron_mass );
189     FourVectorM dilepton = ele1+ele2;
190     if (dilepton.M() > massLo && dilepton.M() < massHi) {
191     ZEle1 = goodElectrons[i];
192     ZEle2 = goodElectrons[j];
193     ZPt = dilepton.Pt();
194     ZMass = dilepton.M();
195     NZCandidates++;
196     }
197     }
198     }
199    
200     if (NZCandidates != 1) return;
201     Int_t ZDecayType = 0;
202     if (ZMuon1 && ZMuon2) ZDecayType = 13;
203     else if (ZEle1 && ZEle2) ZDecayType = 11;
204     else cout << "Error. Z Leptons not properly found.\n";
205    
206     if(ZDecayType == 13){
207    
208     if(!(ZMuon1->Pt() > 20 && ZMuon2->Pt() > 10) && !(ZMuon2->Pt() > 20 && ZMuon1->Pt() > 10)) return;
209    
210     std::bitset<1024> zmuon1_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZMuon1->Pt(), ZMuon1->Eta(), ZMuon1->Phi());
211     std::bitset<1024> zmuon2_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZMuon2->Pt(), ZMuon2->Eta(), ZMuon2->Phi());
212    
213     Bool_t pass_dimuon_trigger1 = (ZDecayType == 13)
214     && (fTriggerBits.test(kHLT_Mu17_Mu8))
215     &&
216     (
217     ((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj))
218     && (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj)))
219     ||
220     ((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj))
221     && (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj)))
222     );
223    
224     Bool_t pass_dimuon_trigger2 = (ZDecayType == 13)
225     && (fTriggerBits.test(kHLT_Mu17_TkMu8))
226     &&
227     (
228     ((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj))
229     && (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj)))
230     ||
231     ((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj))
232     && (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj)))
233     );
234    
235    
236     if(!pass_dimuon_trigger1 && !pass_dimuon_trigger2) return;
237     if(!(die_or_dimu == e_dimuon_data)) return;
238     }
239     else if (ZDecayType == 11){
240    
241     if(!(ZEle1->Pt() > 20 && ZEle2->Pt() > 10) && !(ZEle2->Pt() > 20 && ZEle1->Pt() > 10)) return;
242    
243     std::bitset<1024> zele1_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZEle1->Pt(), ZEle1->Eta(), ZEle1->Phi());
244     std::bitset<1024> zele2_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZEle2->Pt(), ZEle2->Eta(), ZEle2->Phi());
245    
246     Bool_t pass_dielectron_trigger =
247     (fTriggerBits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL))
248     &&
249     (
250     ((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj))
251     && (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj)))
252     ||
253     ((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj))
254     && (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj)))
255     );
256    
257    
258     if(!pass_dielectron_trigger) return;
259     if(!(die_or_dimu == e_dielectron_data)) return;
260    
261     }
262     else assert(0);
263    
264     nvertices = fPrimVerts->GetEntries();
265     run_num = GetEventHeader()->RunNum();
266     lumi_sec = GetEventHeader()->LumiSec();
267     evt_num = GetEventHeader()->EvtNum();
268    
269    
270     for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
271     const Electron *ele = fElectrons->At(i);
272    
273     if (ZDecayType == 11 && (ele == ZEle1 || ele == ZEle2)) continue;
274     if (ZEle1 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZEle1->Eta(), ZEle1->Phi()) < 0.3) continue;
275     if (ZEle2 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZEle2->Eta(), ZEle2->Phi()) < 0.3) continue;
276     if (ZMuon1 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZMuon1->Eta(), ZMuon1->Phi()) < 0.3) continue;
277     if (ZMuon2 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZMuon2->Eta(), ZMuon2->Phi()) < 0.3) continue;
278    
279     for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
280    
281     if ((ZDecayType == 11 && nElectrons[denominatorTypeIndex] > 3)
282     ||
283     (ZDecayType == 13 && nElectrons[denominatorTypeIndex] > 1)
284     )
285     continue;
286    
287     if (passElectronDenominatorCuts(ele, fPrimVerts->At(0), denominatorType[denominatorTypeIndex])) {
288    
289     fake_electron_pt[denominatorTypeIndex] = ele->Pt();
290     fake_electron_eta[denominatorTypeIndex] = ele->Eta();
291     fake_electron_phi[denominatorTypeIndex] = ele->Phi();
292    
293     Bool_t passNumerator = kTRUE;
294    
295     ControlFlags ctrl;
296     mithep::ElectronTools::EElectronEffectiveAreaTarget eraEle = mithep::ElectronTools::kEleEAData2012;
297     vector<const mithep::PFCandidate*> photonsToVeto;
298    
299     if(!electronReferencePreSelection(ctrl, ele, fPrimVerts->At(0)).passPre()) passNumerator = kFALSE;
300     if(fabs(ele->Ip3dPVSignificance()) > 4) passNumerator = kFALSE;
301     if(!electronReferenceIDMVASelectionV1(ctrl, ele, fPrimVerts->At(0)).looseID()) passNumerator = kFALSE;
302     if(!electronReferenceIsoSelection(ctrl,ele,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity,eraEle,photonsToVeto).passLooseIso()) passNumerator = kFALSE;
303    
304     fake_electron_pass[denominatorTypeIndex] = passNumerator ? 1 : 0;
305    
306     fOutputTrees[denominatorTypeIndex]->Fill();
307     }//if pass denominator cut
308     } //loop over denominator types
309     } //loop over muons
310     }
311    
312    
313     std::bitset<1024> mithep::ElectronFakeRateMod::GetHLTMatchBits(const Double_t pt, const Double_t eta, const Double_t phi)
314     {
315     std::bitset<1024> object_bitset;
316     const Double_t hltMatchR = 0.2;
317    
318     if(HasHLTInfo()) {
319     const TriggerTable *hltTable = GetHLTTable();
320     assert(hltTable);
321     for(UInt_t itrig=0; itrig<triggerNames.size(); itrig++) {
322     const TriggerName *trigname = hltTable->Get(triggerNames[itrig].c_str());
323     if(!trigname) continue;
324    
325     const TList *list = GetHLTObjectsTable()->GetList(triggerNames[itrig].c_str());
326     if(!list) continue;
327     TIter iter(list->MakeIterator());
328     const TriggerObject *to = dynamic_cast<const TriggerObject*>(iter.Next());
329    
330     while(to) {
331     if(to->IsHLT()) {
332    
333     //assert(triggerObjNames1[itrig].length()>0);
334    
335     if(triggerObjNames1[itrig].length()>0 && (triggerObjNames1[itrig] == string(to->ModuleName()))) {
336     if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
337     object_bitset.set(triggerObjIds1[itrig]);
338     }
339     }
340    
341     if(triggerObjNames2[itrig].length()>0 && (triggerObjNames2[itrig] == string(to->ModuleName()))) {
342     if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
343     object_bitset.set(triggerObjIds2[itrig]);
344     }
345     }
346    
347     }
348     to = dynamic_cast<const TriggerObject*>(iter.Next());
349     }
350     }
351     }
352    
353     return object_bitset;
354     }
355    
356     Bool_t passElectronDenominatorCuts(const mithep::Electron *ele, const mithep::Vertex * vtx, Int_t DenominatorType) {
357    
358     Bool_t pass = kTRUE;
359    
360     if (DenominatorType == 1) {
361     if (fabs(ele->Ip3dPVSignificance()) > 100) pass = kFALSE;
362     if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
363     if (ele->Pt() < 7) pass = kFALSE;
364     if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
365     if (fabs(ele->BestTrk()->D0Corrected(*vtx)) > 0.5 ||
366     fabs(ele->BestTrk()->DzCorrected(*vtx)) > 1.0 ) pass = kFALSE;
367     }
368     else if (DenominatorType == 2) {
369     if (fabs(ele->Ip3dPVSignificance()) > 4) pass = kFALSE;
370     if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
371     if (ele->Pt() < 7) pass = kFALSE;
372     if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
373     if (fabs(ele->BestTrk()->D0Corrected(*vtx)) > 0.5 ||
374     fabs(ele->BestTrk()->DzCorrected(*vtx)) > 1.0 ) pass = kFALSE;
375     }
376     else if (DenominatorType == 3) {
377     if (fabs(ele->Ip3dPVSignificance()) > 100) pass = kFALSE;
378     if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
379     if (ele->Pt() < 7) pass = kFALSE;
380     if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
381     }
382     else if (DenominatorType == 4) {
383     if (fabs(ele->Ip3dPVSignificance()) > 4) pass = kFALSE;
384     if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
385     if (ele->Pt() < 7) pass = kFALSE;
386     if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
387     }
388     else
389     assert(0);
390    
391     return pass;
392     }
393    
394    
395