ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/MuonFakeRateMod.cc
Revision: 1.8
Committed: Wed Jul 18 12:35:27 2012 UTC (12 years, 10 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.7: +2 -74 lines
Log Message:
added electron energy stuff and miscellaneous other stuff

File Contents

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