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

File Contents

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