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 |
|
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 |
|
294 |
fake_muon_pass[denominatorTypeIndex] = passNumerator ? 1 : 0;
|
295 |
|
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 |
|