ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/ElectronFakeRateMod.cc
Revision: 1.2
Committed: Thu Jul 5 14:26:26 2012 UTC (12 years, 10 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.1: +7 -3 lines
Log Message:
updated tag and probe and fake rate code

File Contents

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