ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/MuonFakeRateMod.cc
Revision: 1.7
Committed: Wed Jul 18 09:05:38 2012 UTC (12 years, 10 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.6: +1 -0 lines
Log Message:
corrected pfnopileupbug with fake rates

File Contents

# User Rev Content
1 anlevin 1.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 anlevin 1.5
34 anlevin 1.1 denominatorType.push_back(1);
35     denominatorType.push_back(2);
36    
37     nProbes = 0;
38    
39     ReqBranch("PFMet", fPFMet);
40    
41     ReqBranch(mithep::Names::gkMvfConversionBrn, fConversions);
42    
43     fMuonName = mithep::Names::gkMuonBrn;
44     ReqBranch(fMuonName, fMuons);
45    
46     fElectronName = mithep::Names::gkElectronBrn;
47     ReqBranch(fElectronName, fElectrons);
48    
49     fPFCandidateName = mithep::Names::gkPFCandidatesBrn;
50     ReqBranch(fPFCandidateName, fPFCandidates);
51    
52     fPrimVtxName = mithep::Names::gkPVBrn;
53     ReqBranch(fPrimVtxName, fPrimVerts);
54    
55     fPUEnergyDensityName = mithep::Names::gkPileupEnergyDensityBrn;
56     ReqBranch(fPUEnergyDensityName, fPileupEnergyDensity);
57    
58     fTrigMaskName = mithep::Names::gkHltBitBrn;
59     ReqBranch(fTrigMaskName, fTrigMask);
60    
61     fPileupName = mithep::Names::gkPileupInfoBrn;
62     ReqBranch(fPileupName, fPileup);
63    
64     fOutputFile = new TFile(fOutputName, "RECREATE");
65     fOutputTrees = vector<TTree*>(denominatorType.size(),0);
66     for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
67     fOutputTrees[denominatorTypeIndex] = new TTree(TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])),TString("denominator_v"+IntToString(denominatorType[denominatorTypeIndex])));
68     }
69    
70     fake_muon_pt = vector<Float_t>(denominatorType.size(),0);
71     fake_muon_eta = vector<Float_t>(denominatorType.size(),0);
72     fake_muon_phi = vector<Float_t>(denominatorType.size(),0);
73     fake_muon_pass = vector<UInt_t>(denominatorType.size(),0);
74    
75    
76     for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
77     fOutputTrees[denominatorTypeIndex]->Branch("nvertices",&nvertices,"nvertices/i");
78     fOutputTrees[denominatorTypeIndex]->Branch("evt_num",&evt_num,"evt_num/i");
79     fOutputTrees[denominatorTypeIndex]->Branch("lumi_sec",&lumi_sec,"lumi_sec/i");
80     fOutputTrees[denominatorTypeIndex]->Branch("run_num",&run_num,"run_num/i");
81    
82     fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_pt"),&fake_muon_pt[denominatorTypeIndex],TString("fake_muon_pt/F"));
83     fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_eta"),&fake_muon_eta[denominatorTypeIndex],TString("fake_muon_eta/F"));
84     fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_phi"),&fake_muon_phi[denominatorTypeIndex],TString("fake_muon_phi/F"));
85     fOutputTrees[denominatorTypeIndex]->Branch(TString("fake_muon_pass"),&fake_muon_pass[denominatorTypeIndex],TString("fake_muon_pass/i"));
86     }
87    
88     rlrm.AddJSONFile("/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions12/8TeV/Prompt/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt");
89 anlevin 1.5
90 anlevin 1.1 }
91    
92     void mithep::MuonFakeRateMod::SlaveTerminate()
93     {
94     fOutputFile->cd();
95     for (UInt_t denominatorTypeIndex = 0; denominatorTypeIndex < denominatorType.size(); ++denominatorTypeIndex) {
96     fOutputTrees[denominatorTypeIndex]->Print();
97    
98     }
99     fOutputFile->Write();
100     fOutputFile->Close();
101    
102     }
103    
104    
105     void mithep::MuonFakeRateMod::Process()
106     {
107     vector<UInt_t> nLeptons(denominatorType.size(),0);
108     vector<UInt_t> nMuons(denominatorType.size(),0);
109    
110    
111     gDebugMask = mithep::Debug::kAnalysis; // debug message category
112     gDebugLevel = 1; // higher level allows more messages to print
113    
114     LoadBranch(fMuonName);
115     LoadBranch(fElectronName);
116     LoadBranch(fPFCandidateName);
117     LoadBranch(fPrimVtxName);
118     LoadBranch(fPUEnergyDensityName);
119     LoadBranch(fTrigMaskName);
120     LoadBranch("PFMet");
121    
122     LoadBranch(mithep::Names::gkMvfConversionBrn);
123     if(store_npu)
124     LoadBranch(fPileupName);
125    
126     mithep::RunLumiRangeMap::RunLumiPairType rl(GetEventHeader()->RunNum(), GetEventHeader()->LumiSec());
127     if(useJSON && !rlrm.HasRunLumi(rl)) return;
128    
129     fillTriggerBits( GetHLTTable(), fTrigMask, fTriggerBits );
130    
131 anlevin 1.7 PFnoPUflag.clear();
132 anlevin 1.1 makePFnoPUArray(fPFCandidates, PFnoPUflag, fPrimVerts );
133    
134     if(fPFMet->At(0)->Et() >25) return;
135    
136     vector<const Muon*> goodMuons;
137     vector<const Electron*> goodElectrons;
138    
139     for(Int_t i=0; i<fMuons->GetEntries(); i++) {
140     const Muon *mu = fMuons->At(i);
141     for ( UInt_t denominatorTypeIndex = 0 ; denominatorTypeIndex < denominatorType.size() ; ++denominatorTypeIndex ) {
142     if (passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[denominatorTypeIndex]))
143     nMuons[denominatorTypeIndex]++;
144     }
145 anlevin 1.3
146 anlevin 1.1 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...
147     passMuonDenominatorCuts(mu, fPrimVerts->At(0), denominatorType[1])
148     &&
149 anlevin 1.4 muonPOG2012CutBasedIDTight(mu,fPrimVerts->At(0),fPFCandidates,fPileupEnergyDensity)
150 anlevin 1.1 )
151     goodMuons.push_back(mu);
152    
153     }
154     for(Int_t i=0; i<fElectrons->GetEntries(); i++) {
155     const Electron *ele = fElectrons->At(i);
156 dkralph 1.6 if(electronPOG2012CutBasedIDMedium(ele, fPrimVerts->At(0), fPFCandidates, fConversions, fPileupEnergyDensity)) goodElectrons.push_back(ele);
157 anlevin 1.1
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 anlevin 1.5
297     fake_muon_pass[denominatorTypeIndex] = passNumerator ? 1 : 0;
298 anlevin 1.1
299     fOutputTrees[denominatorTypeIndex]->Fill();
300    
301     }//if pass denominator cut
302     } //loop over denominator types
303     } //loop over muons
304     }
305    
306     unsigned makePFnoPUArray(const mithep::Array<mithep::PFCandidate> * fPFCandidates,
307     vector<bool> &pfNoPileUpflag,
308     const mithep::Array<mithep::Vertex> * vtxArr )
309     {
310     for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
311     const mithep::PFCandidate *pf = fPFCandidates->At(i);
312     assert(pf);
313    
314     if(pf->PFType() == mithep::PFCandidate::eHadron) {
315     if(pf->HasTrackerTrk() &&
316     vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
317     vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
318    
319     pfNoPileUpflag.push_back(1);
320    
321     } else {
322    
323     Bool_t vertexFound = kFALSE;
324     const mithep::Vertex *closestVtx = 0;
325     Double_t dzmin = 10000;
326    
327     // loop over vertices
328     for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
329     const mithep::Vertex *vtx = vtxArr->At(j);
330     assert(vtx);
331    
332     if(pf->HasTrackerTrk() &&
333     vtx->HasTrack(pf->TrackerTrk()) &&
334     vtx->TrackWeight(pf->TrackerTrk()) > 0) {
335     vertexFound = kTRUE;
336     closestVtx = vtx;
337     break;
338     }
339     Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
340     if(dz < dzmin) {
341     closestVtx = vtx;
342     dzmin = dz;
343     }
344     }
345    
346     // if(fCheckClosestZVertex) {
347     if(1) {
348     // Fallback: if track is not associated with any vertex,
349     // associate it with the vertex closest in z
350     if(vertexFound || closestVtx != vtxArr->At(0)) {
351     // pfPileUp->Add(pf);
352     pfNoPileUpflag.push_back(0);
353     } else {
354     pfNoPileUpflag.push_back(1);
355     }
356     } else {
357     if(vertexFound && closestVtx != vtxArr->At(0)) {
358     // pfPileUp->Add(pf);
359     pfNoPileUpflag.push_back(0);
360     } else {
361     // PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
362     pfNoPileUpflag.push_back(1);
363     }
364     }
365     } //hadron & trk stuff
366     } else { // hadron
367     // PFCandidate * ptr = pfNoPileUp->AddNew();
368     pfNoPileUpflag.push_back(1);
369     }
370     } // Loop over PF candidates
371    
372     return pfNoPileUpflag.size();
373     }
374    
375     std::bitset<1024> mithep::MuonFakeRateMod::GetHLTMatchBits(const Double_t pt, const Double_t eta, const Double_t phi)
376     {
377     std::bitset<1024> object_bitset;
378     const Double_t hltMatchR = 0.2;
379    
380     if(HasHLTInfo()) {
381     const TriggerTable *hltTable = GetHLTTable();
382     assert(hltTable);
383     for(UInt_t itrig=0; itrig<triggerNames.size(); itrig++) {
384     const TriggerName *trigname = hltTable->Get(triggerNames[itrig].c_str());
385     if(!trigname) continue;
386    
387     const TList *list = GetHLTObjectsTable()->GetList(triggerNames[itrig].c_str());
388     if(!list) continue;
389     TIter iter(list->MakeIterator());
390     const TriggerObject *to = dynamic_cast<const TriggerObject*>(iter.Next());
391    
392     while(to) {
393     if(to->IsHLT()) {
394    
395     //assert(triggerObjNames1[itrig].length()>0);
396    
397     if(triggerObjNames1[itrig].length()>0 && (triggerObjNames1[itrig] == string(to->ModuleName()))) {
398     if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
399     object_bitset.set(triggerObjIds1[itrig]);
400     }
401     }
402    
403     if(triggerObjNames2[itrig].length()>0 && (triggerObjNames2[itrig] == string(to->ModuleName()))) {
404     if(MathUtils::DeltaR(phi,eta,to->Phi(),to->Eta()) < hltMatchR){
405     object_bitset.set(triggerObjIds2[itrig]);
406     }
407     }
408    
409     }
410     to = dynamic_cast<const TriggerObject*>(iter.Next());
411     }
412     }
413     }
414    
415     return object_bitset;
416     }
417    
418     Bool_t passMuonDenominatorCuts(const mithep::Muon *mu, const mithep::Vertex * vtx, Int_t DenominatorType) {
419    
420     Bool_t pass = kTRUE;
421    
422     if (DenominatorType == 1) {
423    
424     if (fabs(mu->Eta()) > 2.4) pass = kFALSE;
425    
426     if (mu->Pt() < 5) pass = kFALSE;
427    
428     if(!mu->IsTrackerMuon() && !mu->IsGlobalMuon()) pass = kFALSE;
429    
430     if(fabs(mu->Ip3dPVSignificance()) > 100) pass = kFALSE;
431    
432     if(!(mu->HasTrackerTrk() && fabs(mu->TrackerTrk()->D0Corrected(*vtx)) < 0.5 && fabs(mu->TrackerTrk()->DzCorrected(*vtx)) < 1.0)) pass = kFALSE;
433    
434     }
435     else if (DenominatorType == 2) {
436    
437     if (fabs(mu->Eta()) > 2.4) pass = kFALSE;
438    
439     if (mu->Pt() < 5) pass = kFALSE;
440    
441     if(!mu->IsTrackerMuon() && !mu->IsGlobalMuon())
442     pass = kFALSE;
443    
444     if(fabs(mu->Ip3dPVSignificance()) > 4) pass = kFALSE;
445    
446     if(!(mu->HasTrackerTrk() && fabs(mu->TrackerTrk()->D0Corrected(*vtx)) < 0.5 && fabs(mu->TrackerTrk()->DzCorrected(*vtx)) < 1.0)) pass = kFALSE;
447    
448     }
449     else
450     assert(0);
451    
452     return pass;
453     }
454    
455