ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/ElectronFakeRateMod.cc
Revision: 1.9
Committed: Tue Oct 23 12:43:04 2012 UTC (12 years, 6 months ago) by dkralph
Content type: text/plain
Branch: MAIN
CVS Tags: compiled, HEAD
Changes since 1.8: +2 -1 lines
Log Message:
*** empty log message ***

File Contents

# Content
1 #include "ElectronFakeRateMod.h"
2 #include "SelectionFuncs.h"
3
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 PFnoPUflag.clear();
135 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 const Muon *mu = fMuons->At(i);
142 if( muonPOG2012CutBasedIDTight(mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity) )
143 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 // if(!fails_some_denominator && electronPOG2012CutBasedIDMedium(ele, fPrimVerts->At(0), fPFCandidates, fConversions, fPileupEnergyDensity->At(0)->RhoLowEta(), mithep::ElectronTools::kEleEAData2011))
156 assert(0); // doesn't compile
157 goodElectrons.push_back(ele);
158 }
159
160 Int_t NZCandidates = 0;
161 const mithep::Muon* ZMuon1 = 0;
162 const mithep::Muon* ZMuon2 = 0;
163 const mithep::Electron *ZEle1 = 0;
164 const mithep::Electron *ZEle2 = 0;
165 Double_t ZPt = 0;
166 Double_t ZMass = 0;
167 for(UInt_t i=0; i<goodMuons.size(); i++) {
168 FourVectorM mu1;
169 mu1.SetCoordinates(goodMuons[i]->Pt(), goodMuons[i]->Eta(), goodMuons[i]->Phi(), g_muon_mass );
170 for(UInt_t j=i+1; j<goodMuons.size(); j++) {
171 FourVectorM mu2;
172 mu2.SetCoordinates(goodMuons[j]->Pt(), goodMuons[j]->Eta(), goodMuons[j]->Phi(), g_muon_mass );
173 FourVectorM dilepton = mu1+mu2;
174 if (dilepton.M() > massLo && dilepton.M() < massHi) {
175 ZMuon1 = goodMuons[i];
176 ZMuon2 = goodMuons[j];
177 ZPt = dilepton.Pt();
178 ZMass = dilepton.M();
179 NZCandidates++;
180 }
181 }
182 }
183
184 for(UInt_t i=0; i<goodElectrons.size(); i++) {
185 FourVectorM ele1;
186 ele1.SetCoordinates(goodElectrons[i]->Pt(), goodElectrons[i]->Eta(), goodElectrons[i]->Phi(),g_electron_mass );
187 for(UInt_t j=i+1; j<goodElectrons.size(); j++) {
188 FourVectorM ele2;
189 ele2.SetCoordinates(goodElectrons[j]->Pt(), goodElectrons[j]->Eta(), goodElectrons[j]->Phi(), g_electron_mass );
190 FourVectorM dilepton = ele1+ele2;
191 if (dilepton.M() > massLo && dilepton.M() < massHi) {
192 ZEle1 = goodElectrons[i];
193 ZEle2 = goodElectrons[j];
194 ZPt = dilepton.Pt();
195 ZMass = dilepton.M();
196 NZCandidates++;
197 }
198 }
199 }
200
201 if (NZCandidates != 1) return;
202 Int_t ZDecayType = 0;
203 if (ZMuon1 && ZMuon2) ZDecayType = 13;
204 else if (ZEle1 && ZEle2) ZDecayType = 11;
205 else cout << "Error. Z Leptons not properly found.\n";
206
207 if(ZDecayType == 13){
208
209 if(!(ZMuon1->Pt() > 20 && ZMuon2->Pt() > 10) && !(ZMuon2->Pt() > 20 && ZMuon1->Pt() > 10)) return;
210
211 std::bitset<1024> zmuon1_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZMuon1->Pt(), ZMuon1->Eta(), ZMuon1->Phi());
212 std::bitset<1024> zmuon2_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZMuon2->Pt(), ZMuon2->Eta(), ZMuon2->Phi());
213
214 Bool_t pass_dimuon_trigger1 = (ZDecayType == 13)
215 && (fTriggerBits.test(kHLT_Mu17_Mu8))
216 &&
217 (
218 ((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj))
219 && (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj)))
220 ||
221 ((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj))
222 && (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj)))
223 );
224
225 Bool_t pass_dimuon_trigger2 = (ZDecayType == 13)
226 && (fTriggerBits.test(kHLT_Mu17_TkMu8))
227 &&
228 (
229 ((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj))
230 && (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj)))
231 ||
232 ((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj))
233 && (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj)))
234 );
235
236
237 if(!pass_dimuon_trigger1 && !pass_dimuon_trigger2) return;
238 if(!(die_or_dimu == e_dimuon_data)) return;
239 }
240 else if (ZDecayType == 11){
241
242 if(!(ZEle1->Pt() > 20 && ZEle2->Pt() > 10) && !(ZEle2->Pt() > 20 && ZEle1->Pt() > 10)) return;
243
244 std::bitset<1024> zele1_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZEle1->Pt(), ZEle1->Eta(), ZEle1->Phi());
245 std::bitset<1024> zele2_matchbits = mithep::ElectronFakeRateMod::GetHLTMatchBits(ZEle2->Pt(), ZEle2->Eta(), ZEle2->Phi());
246
247 Bool_t pass_dielectron_trigger =
248 (fTriggerBits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL))
249 &&
250 (
251 ((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj))
252 && (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj)))
253 ||
254 ((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj))
255 && (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj)))
256 );
257
258
259 if(!pass_dielectron_trigger) return;
260 if(!(die_or_dimu == e_dielectron_data)) return;
261
262 }
263 else assert(0);
264
265 nvertices = fPrimVerts->GetEntries();
266 run_num = GetEventHeader()->RunNum();
267 lumi_sec = GetEventHeader()->LumiSec();
268 evt_num = GetEventHeader()->EvtNum();
269
270
271 for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
272 const Electron *ele = fElectrons->At(i);
273
274 if (ZDecayType == 11 && (ele == ZEle1 || ele == ZEle2)) continue;
275 if (ZEle1 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZEle1->Eta(), ZEle1->Phi()) < 0.3) continue;
276 if (ZEle2 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZEle2->Eta(), ZEle2->Phi()) < 0.3) continue;
277 if (ZMuon1 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZMuon1->Eta(), ZMuon1->Phi()) < 0.3) continue;
278 if (ZMuon2 && mithep::MathUtils::DeltaR(ele->Eta(), ele->Phi(), ZMuon2->Eta(), ZMuon2->Phi()) < 0.3) continue;
279
280 for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
281
282 if ((ZDecayType == 11 && nElectrons[denominatorTypeIndex] > 3)
283 ||
284 (ZDecayType == 13 && nElectrons[denominatorTypeIndex] > 1)
285 )
286 continue;
287
288 if (passElectronDenominatorCuts(ele, fPrimVerts->At(0), denominatorType[denominatorTypeIndex])) {
289
290 fake_electron_pt[denominatorTypeIndex] = ele->Pt();
291 fake_electron_eta[denominatorTypeIndex] = ele->Eta();
292 fake_electron_phi[denominatorTypeIndex] = ele->Phi();
293
294 Bool_t passNumerator = kTRUE;
295
296 ControlFlags ctrl;
297 mithep::ElectronTools::EElectronEffectiveAreaTarget eraEle = mithep::ElectronTools::kEleEAData2012;
298 vector<const mithep::PFCandidate*> photonsToVeto;
299
300 if(!electronReferencePreSelection(ctrl, ele, fPrimVerts->At(0)).passPre()) passNumerator = kFALSE;
301 if(fabs(ele->Ip3dPVSignificance()) > 4) passNumerator = kFALSE;
302 if(!electronReferenceIDMVASelectionV1(ctrl, ele, fPrimVerts->At(0)).looseID()) passNumerator = kFALSE;
303 if(!electronReferenceIsoSelection(ctrl,ele,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity,eraEle,photonsToVeto).passLooseIso()) passNumerator = kFALSE;
304
305 fake_electron_pass[denominatorTypeIndex] = passNumerator ? 1 : 0;
306
307 fOutputTrees[denominatorTypeIndex]->Fill();
308 }//if pass denominator cut
309 } //loop over denominator types
310 } //loop over muons
311 }
312
313
314 std::bitset<1024> mithep::ElectronFakeRateMod::GetHLTMatchBits(const Double_t pt, const Double_t eta, const Double_t phi)
315 {
316 std::bitset<1024> object_bitset;
317 const Double_t hltMatchR = 0.2;
318
319 if(HasHLTInfo()) {
320 const TriggerTable *hltTable = GetHLTTable();
321 assert(hltTable);
322 for(UInt_t itrig=0; itrig<triggerNames.size(); itrig++) {
323 const TriggerName *trigname = hltTable->Get(triggerNames[itrig].c_str());
324 if(!trigname) continue;
325
326 const TList *list = GetHLTObjectsTable()->GetList(triggerNames[itrig].c_str());
327 if(!list) continue;
328 TIter iter(list->MakeIterator());
329 const TriggerObject *to = dynamic_cast<const TriggerObject*>(iter.Next());
330
331 while(to) {
332 if(to->IsHLT()) {
333
334 //assert(triggerObjNames1[itrig].length()>0);
335
336 if(triggerObjNames1[itrig].length()>0 && (triggerObjNames1[itrig] == string(to->ModuleName()))) {
337 if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
338 object_bitset.set(triggerObjIds1[itrig]);
339 }
340 }
341
342 if(triggerObjNames2[itrig].length()>0 && (triggerObjNames2[itrig] == string(to->ModuleName()))) {
343 if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
344 object_bitset.set(triggerObjIds2[itrig]);
345 }
346 }
347
348 }
349 to = dynamic_cast<const TriggerObject*>(iter.Next());
350 }
351 }
352 }
353
354 return object_bitset;
355 }
356
357 Bool_t passElectronDenominatorCuts(const mithep::Electron *ele, const mithep::Vertex * vtx, Int_t DenominatorType) {
358
359 Bool_t pass = kTRUE;
360
361 if (DenominatorType == 1) {
362 if (fabs(ele->Ip3dPVSignificance()) > 100) pass = kFALSE;
363 if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
364 if (ele->Pt() < 7) pass = kFALSE;
365 if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
366 if (fabs(ele->BestTrk()->D0Corrected(*vtx)) > 0.5 ||
367 fabs(ele->BestTrk()->DzCorrected(*vtx)) > 1.0 ) pass = kFALSE;
368 }
369 else if (DenominatorType == 2) {
370 if (fabs(ele->Ip3dPVSignificance()) > 4) pass = kFALSE;
371 if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
372 if (ele->Pt() < 7) pass = kFALSE;
373 if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
374 if (fabs(ele->BestTrk()->D0Corrected(*vtx)) > 0.5 ||
375 fabs(ele->BestTrk()->DzCorrected(*vtx)) > 1.0 ) pass = kFALSE;
376 }
377 else if (DenominatorType == 3) {
378 if (fabs(ele->Ip3dPVSignificance()) > 100) pass = kFALSE;
379 if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
380 if (ele->Pt() < 7) pass = kFALSE;
381 if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
382 }
383 else if (DenominatorType == 4) {
384 if (fabs(ele->Ip3dPVSignificance()) > 4) pass = kFALSE;
385 if (fabs(ele->Eta()) > 2.5) pass = kFALSE;
386 if (ele->Pt() < 7) pass = kFALSE;
387 if (ele->CorrectedNExpectedHitsInner() > 1) pass = kFALSE;
388 }
389 else
390 assert(0);
391
392 return pass;
393 }
394
395
396