ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/MuonFakeRateMod.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: +6 -3 lines
Log Message:
updated tag and probe and fake rate code

File Contents

# Content
1 #include "MuonFakeRateMod.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 passMuonDenominatorCuts(const mithep::Muon *mu, 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_muon_mass = 105.658369e-3;
29
30
31 void mithep::MuonFakeRateMod::SlaveBegin()
32 {
33 denominatorType.push_back(1);
34 denominatorType.push_back(2);
35
36 nProbes = 0;
37
38 ReqBranch("PFMet", fPFMet);
39
40 ReqBranch(mithep::Names::gkMvfConversionBrn, fConversions);
41
42 fMuonName = mithep::Names::gkMuonBrn;
43 ReqBranch(fMuonName, fMuons);
44
45 fElectronName = mithep::Names::gkElectronBrn;
46 ReqBranch(fElectronName, fElectrons);
47
48 fPFCandidateName = mithep::Names::gkPFCandidatesBrn;
49 ReqBranch(fPFCandidateName, fPFCandidates);
50
51 fPrimVtxName = mithep::Names::gkPVBrn;
52 ReqBranch(fPrimVtxName, fPrimVerts);
53
54 fPUEnergyDensityName = mithep::Names::gkPileupEnergyDensityBrn;
55 ReqBranch(fPUEnergyDensityName, fPileupEnergyDensity);
56
57 fTrigMaskName = mithep::Names::gkHltBitBrn;
58 ReqBranch(fTrigMaskName, fTrigMask);
59
60 fPileupName = mithep::Names::gkPileupInfoBrn;
61 ReqBranch(fPileupName, fPileup);
62
63 fOutputFile = new TFile(fOutputName, "RECREATE");
64 fOutputTrees = vector<TTree*>(denominatorType.size(),0);
65 for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
66 fOutputTrees[denominatorTypeIndex] = new TTree(TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])),TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])));
67 }
68
69 fake_muon_pt = vector<Float_t>(denominatorType.size(),0);
70 fake_muon_eta = vector<Float_t>(denominatorType.size(),0);
71 fake_muon_phi = vector<Float_t>(denominatorType.size(),0);
72 fake_muon_pass = vector<UInt_t>(denominatorType.size(),0);
73
74
75 for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
76 fOutputTrees[denominatorTypeIndex]->Branch("nvertices",&nvertices,"nvertices/i");
77 fOutputTrees[denominatorTypeIndex]->Branch("evt_num",&evt_num,"evt_num/i");
78 fOutputTrees[denominatorTypeIndex]->Branch("lumi_sec",&lumi_sec,"lumi_sec/i");
79 fOutputTrees[denominatorTypeIndex]->Branch("run_num",&run_num,"run_num/i");
80
81 fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_pt"),&fake_muon_pt[denominatorTypeIndex],TString("fake_muon_pt/F"));
82 fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_eta"),&fake_muon_eta[denominatorTypeIndex],TString("fake_muon_eta/F"));
83 fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_phi"),&fake_muon_phi[denominatorTypeIndex],TString("fake_muon_phi/F"));
84 fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_pass"),&fake_muon_pass[denominatorTypeIndex],TString("fake_muon_pass/i"));
85 }
86
87 rlrm.AddJSONFile("/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions12/8TeV/Prompt/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt");
88 }
89
90 void mithep::MuonFakeRateMod::SlaveTerminate()
91 {
92 fOutputFile->cd();
93 for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
94 fOutputTrees[denominatorTypeIndex]->Print();
95
96 }
97 fOutputFile->Write();
98 fOutputFile->Close();
99
100 }
101
102
103 void mithep::MuonFakeRateMod::Process()
104 {
105 vector<UInt_t> nLeptons(denominatorType.size(),0);
106 vector<UInt_t> nMuons(denominatorType.size(),0);
107
108
109 gDebugMask = mithep::Debug::kAnalysis; // debug message category
110 gDebugLevel = 1; // higher level allows more messages to print
111
112 LoadBranch(fMuonName);
113 LoadBranch(fElectronName);
114 LoadBranch(fPFCandidateName);
115 LoadBranch(fPrimVtxName);
116 LoadBranch(fPUEnergyDensityName);
117 LoadBranch(fTrigMaskName);
118 LoadBranch("PFMet");
119
120 LoadBranch(mithep::Names::gkMvfConversionBrn);
121 if(store_npu)
122 LoadBranch(fPileupName);
123
124 mithep::RunLumiRangeMap::RunLumiPairType rl(GetEventHeader()->RunNum(), GetEventHeader()->LumiSec());
125 if(useJSON && !rlrm.HasRunLumi(rl)) return;
126
127 fillTriggerBits( GetHLTTable(), fTrigMask, fTriggerBits );
128
129 makePFnoPUArray(fPFCandidates, PFnoPUflag, fPrimVerts );
130
131 if(fPFMet->At(0)->Et() >25) return;
132
133 vector<const Muon*> goodMuons;
134 vector<const Electron*> goodElectrons;
135
136 for(Int_t i=0; i<fMuons->GetEntries(); i++) {
137 const Muon *mu = fMuons->At(i);
138 for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
139 if (passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[denominatorTypeIndex]))
140 nMuons[denominatorTypeIndex]++;
141 }
142 ControlFlags ctrl;
143 ctrl.era = 2012;
144 if ( //maybe you should do a loop to make sure the muon passes all of the denominators? right now this is ok because 2 is tighter...
145 passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[1])
146 &&
147 muon2012CutBasedIDTight(ctrl,mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity,mithep::MuonTools::kMuEAData2012)
148 )
149 goodMuons.push_back(mu);
150
151 }
152 for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
153 const Electron *ele = fElectrons->At(i);
154 ControlFlags ctrl;
155 ctrl.era = 2012;
156 if(electron2012CutBasedIDMedium(ctrl,ele, fPrimVerts->At(0), fPFCandidates, fConversions, fPileupEnergyDensity,mithep::ElectronTools::kEleEAData2012)) goodElectrons.push_back(ele);
157
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(), 105.658369e-3 );
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(), 105.658369e-3 );
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(), 0.51099892e-3 );
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(), 0.51099892e-3 );
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 std::bitset<1024> zmuon1_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZMuon1->Pt(), ZMuon1->Eta(), ZMuon1->Phi());
210 std::bitset<1024> zmuon2_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZMuon2->Pt(), ZMuon2->Eta(), ZMuon2->Phi());
211
212 Bool_t pass_dimuon_trigger1 = (ZDecayType == 13)
213 && (fTriggerBits.test(kHLT_Mu17_Mu8))
214 &&
215 (
216 ((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj))
217 && (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj)))
218 ||
219 ((zmuon1_matchbits.test(kHLT_Mu17_Mu8_Mu2Obj))
220 && (zmuon2_matchbits.test(kHLT_Mu17_Mu8_Mu1Obj)))
221 );
222
223 Bool_t pass_dimuon_trigger2 = (ZDecayType == 13)
224 && (fTriggerBits.test(kHLT_Mu17_TkMu8))
225 &&
226 (
227 ((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj))
228 && (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj)))
229 ||
230 ((zmuon1_matchbits.test(kHLT_Mu17_TkMu8_Mu2Obj))
231 && (zmuon2_matchbits.test(kHLT_Mu17_TkMu8_Mu1Obj)))
232 );
233
234
235 if(!pass_dimuon_trigger1 && !pass_dimuon_trigger2) return;
236 if(!(die_or_dimu == e_dimuon_data)) return;
237 }
238 else if (ZDecayType == 11){
239
240 std::bitset<1024> zele1_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZEle1->Pt(), ZEle1->Eta(), ZEle1->Phi());
241 std::bitset<1024> zele2_matchbits = mithep::MuonFakeRateMod::GetHLTMatchBits(ZEle2->Pt(), ZEle2->Eta(), ZEle2->Phi());
242
243 Bool_t pass_dielectron_trigger =
244 (fTriggerBits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL))
245 &&
246 (
247 ((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj))
248 && (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj)))
249 ||
250 ((zele1_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj))
251 && (zele2_matchbits.test(kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj)))
252 );
253
254 if(!pass_dielectron_trigger) return;
255 if(!(die_or_dimu == e_dielectron_data)) return;
256
257 }
258 else assert(0);
259
260 nvertices = fPrimVerts->GetEntries();
261 run_num = GetEventHeader()->RunNum();
262 lumi_sec = GetEventHeader()->LumiSec();
263 evt_num = GetEventHeader()->EvtNum();
264
265
266 for(Int_t i=0; i<fMuons->GetEntries(); i++) {
267 const Muon *mu = fMuons->At(i);
268
269 if (ZDecayType == 13 && (mu == ZMuon1 || mu == ZMuon2)) continue;
270 if (ZEle1 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZEle1->Eta(), ZEle1->Phi()) < 0.3) continue;
271 if (ZEle2 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZEle2->Eta(), ZEle2->Phi()) < 0.3) continue;
272 if (ZMuon1 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZMuon1->Eta(), ZMuon1->Phi()) < 0.3) continue;
273 if (ZMuon2 && mithep::MathUtils::DeltaR(mu->Eta(), mu->Phi(), ZMuon2->Eta(), ZMuon2->Phi()) < 0.3) continue;
274
275 for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
276 if ((ZDecayType == 11 && nMuons[denominatorTypeIndex] > 1)
277 ||
278 (ZDecayType == 13 && nMuons[denominatorTypeIndex] > 3)
279 ) continue;
280
281 if (passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[denominatorTypeIndex])) {
282
283 fake_muon_pt[denominatorTypeIndex] = mu->Pt();
284 fake_muon_eta[denominatorTypeIndex] = mu->Eta();
285 fake_muon_phi[denominatorTypeIndex] = mu->Phi();
286
287 Bool_t passNumerator = kTRUE;
288 vector<const mithep::PFCandidate*> photonsToVeto;
289
290 ControlFlags ctrl;
291 mithep::MuonTools::EMuonEffectiveAreaTarget eraMu = mithep::MuonTools::kMuEAData2012;
292 if(!muonReferencePreSelection(ctrl, mu, fPrimVerts->At(0), fPFCandidates ).passPre()) passNumerator = kFALSE;
293 if(!muonIDPFSelection(ctrl, mu, fPrimVerts->At(0), fPFCandidates ).looseID()) passNumerator = kFALSE;
294 if(fabs(mu->Ip3dPVSignificance()) > 4) passNumerator = kFALSE;
295 if(!muonReferenceIsoSelection(ctrl,mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity,eraMu,photonsToVeto).passLooseIso()) passNumerator = kFALSE;
296
297 if (passNumerator)
298 fake_muon_pass[denominatorTypeIndex] = passNumerator ? 1 : 0;
299
300 fOutputTrees[denominatorTypeIndex]->Fill();
301
302 }//if pass denominator cut
303 } //loop over denominator types
304 } //loop over muons
305 }
306
307 unsigned makePFnoPUArray(const mithep::Array<mithep::PFCandidate> * fPFCandidates,
308 vector<bool> &pfNoPileUpflag,
309 const mithep::Array<mithep::Vertex> * vtxArr )
310 {
311 for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
312 const mithep::PFCandidate *pf = fPFCandidates->At(i);
313 assert(pf);
314
315 if(pf->PFType() == mithep::PFCandidate::eHadron) {
316 if(pf->HasTrackerTrk() &&
317 vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
318 vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
319
320 pfNoPileUpflag.push_back(1);
321
322 } else {
323
324 Bool_t vertexFound = kFALSE;
325 const mithep::Vertex *closestVtx = 0;
326 Double_t dzmin = 10000;
327
328 // loop over vertices
329 for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
330 const mithep::Vertex *vtx = vtxArr->At(j);
331 assert(vtx);
332
333 if(pf->HasTrackerTrk() &&
334 vtx->HasTrack(pf->TrackerTrk()) &&
335 vtx->TrackWeight(pf->TrackerTrk()) > 0) {
336 vertexFound = kTRUE;
337 closestVtx = vtx;
338 break;
339 }
340 Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
341 if(dz < dzmin) {
342 closestVtx = vtx;
343 dzmin = dz;
344 }
345 }
346
347 // if(fCheckClosestZVertex) {
348 if(1) {
349 // Fallback: if track is not associated with any vertex,
350 // associate it with the vertex closest in z
351 if(vertexFound || closestVtx != vtxArr->At(0)) {
352 // pfPileUp->Add(pf);
353 pfNoPileUpflag.push_back(0);
354 } else {
355 pfNoPileUpflag.push_back(1);
356 }
357 } else {
358 if(vertexFound && closestVtx != vtxArr->At(0)) {
359 // pfPileUp->Add(pf);
360 pfNoPileUpflag.push_back(0);
361 } else {
362 // PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
363 pfNoPileUpflag.push_back(1);
364 }
365 }
366 } //hadron & trk stuff
367 } else { // hadron
368 // PFCandidate * ptr = pfNoPileUp->AddNew();
369 pfNoPileUpflag.push_back(1);
370 }
371 } // Loop over PF candidates
372
373 return pfNoPileUpflag.size();
374 }
375
376 std::bitset<1024> mithep::MuonFakeRateMod::GetHLTMatchBits(const Double_t pt, const Double_t eta, const Double_t phi)
377 {
378 std::bitset<1024> object_bitset;
379 const Double_t hltMatchR = 0.2;
380
381 if(HasHLTInfo()) {
382 const TriggerTable *hltTable = GetHLTTable();
383 assert(hltTable);
384 for(UInt_t itrig=0; itrig<triggerNames.size(); itrig++) {
385 const TriggerName *trigname = hltTable->Get(triggerNames[itrig].c_str());
386 if(!trigname) continue;
387
388 const TList *list = GetHLTObjectsTable()->GetList(triggerNames[itrig].c_str());
389 if(!list) continue;
390 TIter iter(list->MakeIterator());
391 const TriggerObject *to = dynamic_cast<const TriggerObject*>(iter.Next());
392
393 while(to) {
394 if(to->IsHLT()) {
395
396 //assert(triggerObjNames1[itrig].length()>0);
397
398 if(triggerObjNames1[itrig].length()>0 && (triggerObjNames1[itrig] == string(to->ModuleName()))) {
399 if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
400 object_bitset.set(triggerObjIds1[itrig]);
401 }
402 }
403
404 if(triggerObjNames2[itrig].length()>0 && (triggerObjNames2[itrig] == string(to->ModuleName()))) {
405 if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
406 object_bitset.set(triggerObjIds2[itrig]);
407 }
408 }
409
410 }
411 to = dynamic_cast<const TriggerObject*>(iter.Next());
412 }
413 }
414 }
415
416 return object_bitset;
417 }
418
419 Bool_t passMuonDenominatorCuts(const mithep::Muon *mu, const mithep::Vertex * vtx, Int_t DenominatorType) {
420
421 Bool_t pass = kTRUE;
422
423 if (DenominatorType == 1) {
424
425 if (fabs(mu->Eta()) > 2.4) pass = kFALSE;
426
427 if (mu->Pt() < 5) pass = kFALSE;
428
429 if(!mu->IsTrackerMuon() && !mu->IsGlobalMuon()) pass = kFALSE;
430
431 if(fabs(mu->Ip3dPVSignificance()) > 100) pass = kFALSE;
432
433 if(!(mu->HasTrackerTrk() && fabs(mu->TrackerTrk()->D0Corrected(*vtx)) < 0.5 && fabs(mu->TrackerTrk()->DzCorrected(*vtx)) < 1.0)) pass = kFALSE;
434
435 }
436 else if (DenominatorType == 2) {
437
438 if (fabs(mu->Eta()) > 2.4) pass = kFALSE;
439
440 if (mu->Pt() < 5) pass = kFALSE;
441
442 if(!mu->IsTrackerMuon() && !mu->IsGlobalMuon())
443 pass = kFALSE;
444
445 if(fabs(mu->Ip3dPVSignificance()) > 4) pass = kFALSE;
446
447 if(!(mu->HasTrackerTrk() && fabs(mu->TrackerTrk()->D0Corrected(*vtx)) < 0.5 && fabs(mu->TrackerTrk()->DzCorrected(*vtx)) < 1.0)) pass = kFALSE;
448
449 }
450 else
451 assert(0);
452
453 return pass;
454 }
455
456