ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.93
Committed: Fri Oct 18 14:10:13 2013 UTC (11 years, 6 months ago) by ceballos
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.92: +2 -2 lines
Log Message: 
minor fix

File Contents

# User Rev Content
1 ceballos 1.93 // $Id: MuonIDMod.cc,v 1.92 2013/10/18 14:09:34 ceballos Exp $
2 loizides 1.1
3     #include "MitPhysics/Mods/interface/MuonIDMod.h"
4 loizides 1.6 #include "MitCommon/MathTools/interface/MathUtils.h"
5 ceballos 1.38 #include "MitAna/DataTree/interface/MuonFwd.h"
6     #include "MitAna/DataTree/interface/ElectronFwd.h"
7 loizides 1.23 #include "MitAna/DataTree/interface/VertexCol.h"
8 loizides 1.6 #include "MitPhysics/Init/interface/ModNames.h"
9 loizides 1.1
10     using namespace mithep;
11    
12     ClassImp(mithep::MuonIDMod)
13    
14     //--------------------------------------------------------------------------------------------------
15     MuonIDMod::MuonIDMod(const char *name, const char *title) :
16     BaseMod(name,title),
17 sixie 1.66 fPrintMVADebugInfo(kFALSE),
18 loizides 1.6 fMuonBranchName(Names::gkMuonBrn),
19     fCleanMuonsName(ModNames::gkCleanMuonsName),
20 ceballos 1.39 fNonIsolatedMuonsName("random"),
21     fNonIsolatedElectronsName("random"),
22 ceballos 1.35 fVertexName(ModNames::gkGoodVertexesName),
23 ceballos 1.42 fBeamSpotName(Names::gkBeamSpotBrn),
24 ceballos 1.38 fTrackName(Names::gkTrackBrn),
25     fPFCandidatesName(Names::gkPFCandidatesBrn),
26 fabstoec 1.85 fPFNoPileUpName("PFNoPileUp"),
27     fPFPileUpName("PFPileUp"),
28 mdecross 1.91 fMuonIDType("NoId"),
29 ceballos 1.49 fMuonIsoType("PFIso"),
30 mdecross 1.91 fMuonClassType("GlobalorTracker"),
31 ceballos 1.2 fTrackIsolationCut(3.0),
32     fCaloIsolationCut(3.0),
33 sixie 1.54 fCombIsolationCut(0.15),
34     fCombRelativeIsolationCut(0.15),
35 ceballos 1.82 fPFIsolationCut(-999.0),
36 mdecross 1.91 fMuonPtMin(10.),
37 loizides 1.21 fApplyD0Cut(kTRUE),
38 ceballos 1.42 fApplyDZCut(kTRUE),
39 ceballos 1.30 fD0Cut(0.020),
40 ceballos 1.50 fDZCut(0.10),
41 ceballos 1.42 fWhichVertex(-1),
42 ceballos 1.30 fEtaCut(2.4),
43 loizides 1.15 fMuIDType(kIdUndef),
44     fMuIsoType(kIsoUndef),
45     fMuClassType(kClassUndef),
46 ceballos 1.14 fMuons(0),
47 loizides 1.15 fVertices(0),
48 ceballos 1.42 fBeamSpot(0),
49 ceballos 1.38 fTracks(0),
50     fPFCandidates(0),
51 ceballos 1.70 fPFNoPileUpCands(0),
52 fabstoec 1.85 fPFPileUpCands(0),
53 ceballos 1.56 fIntRadius(0.0),
54 ceballos 1.39 fNonIsolatedMuons(0),
55 ceballos 1.41 fNonIsolatedElectrons(0),
56     fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
57 sixie 1.59 fPileupEnergyDensity(0),
58     fMuonTools(0),
59     fMuonIDMVA(0),
60 mdecross 1.91 fPVName(Names::gkPVBeamSpotBrn),
61     fTheRhoType(RhoUtilities::DEFAULT)
62 loizides 1.1 {
63     // Constructor.
64     }
65    
66     //--------------------------------------------------------------------------------------------------
67     void MuonIDMod::Process()
68     {
69     // Process entries of the tree.
70    
71 mingyang 1.87 if(fCleanMuonsName.CompareTo("HggLeptonTagMuons") == 0 ){
72     LoadEventObject(fPVName,fVertices);
73     }
74     else{
75     fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
76     }
77    
78 ceballos 1.38 if(fMuIsoType != kPFIsoNoL) {
79     LoadEventObject(fMuonBranchName, fMuons);
80     }
81     else {
82     fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
83     }
84 ceballos 1.42 LoadEventObject(fBeamSpotName, fBeamSpot);
85 ceballos 1.38 LoadEventObject(fTrackName, fTracks);
86     LoadEventObject(fPFCandidatesName, fPFCandidates);
87 sixie 1.59 if(fMuIsoType == kTrackCaloSliding ||
88 fabstoec 1.81 fMuIsoType == kCombinedRelativeConeAreaCorrected ||
89     fMuIsoType == kPFIsoEffectiveAreaCorrected ||
90 ceballos 1.74 fMuIsoType == kMVAIso_BDTG_IDIso ||
91 ceballos 1.75 fMuIsoType == kIsoRingsV0_BDTG_Iso ||
92     fMuIsoType == kIsoDeltaR
93 sixie 1.59 ) {
94 ceballos 1.41 LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
95     }
96 fabstoec 1.85 if(fMuIsoType == kPFRadialIso || fMuIsoType == kIsoDeltaR || fMuIsoType == kPFIsoBetaPUCorrected){
97     // Name is hardcoded, can be changed if someone feels to do it *** did it--Heng
98     fPFNoPileUpCands = GetObjThisEvt<PFCandidateCol>(fPFNoPileUpName);
99     fPFPileUpCands = GetObjThisEvt<PFCandidateCol>(fPFPileUpName);
100 ceballos 1.70 }
101 sixie 1.59
102 loizides 1.7 MuonOArr *CleanMuons = new MuonOArr;
103     CleanMuons->SetName(fCleanMuonsName);
104 loizides 1.1
105 fabstoec 1.69 for (UInt_t i=0; i<fMuons->GetEntries() && fVertices->GetEntries() > 0 ; ++i) {
106 loizides 1.8 const Muon *mu = fMuons->At(i);
107 loizides 1.6
108     Bool_t pass = kFALSE;
109 ceballos 1.30 Double_t pt = 0; // make sure pt is taken from the correct track!
110     Double_t eta = 0; // make sure eta is taken from the correct track!
111 loizides 1.6 switch (fMuClassType) {
112     case kAll:
113     pass = kTRUE;
114 ceballos 1.30 if (mu->HasTrk()) {
115     pt = mu->Pt();
116     eta = TMath::Abs(mu->Eta());
117     }
118 loizides 1.6 break;
119     case kGlobal:
120 ceballos 1.31 pass = mu->HasGlobalTrk() && mu->IsTrackerMuon();
121 ceballos 1.32 if (pass && mu->TrackerTrk()) {
122 ceballos 1.30 pt = mu->TrackerTrk()->Pt();
123     eta = TMath::Abs(mu->TrackerTrk()->Eta());
124     }
125 ceballos 1.32 else {
126     pt = mu->Pt();
127     eta = TMath::Abs(mu->Eta());
128     }
129 ceballos 1.30 break;
130 mdecross 1.91 case kGlobalorTracker:
131     pass = mu->HasGlobalTrk() || mu->IsTrackerMuon();
132     if (pass && mu->TrackerTrk()) {
133     pt = mu->TrackerTrk()->Pt();
134     eta = TMath::Abs(mu->TrackerTrk()->Eta());
135     }
136     else{
137     pt = mu->Pt();
138     eta = TMath::Abs(mu->Eta());
139     }
140 ceballos 1.53 case kGlobalTracker:
141     pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
142     (mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
143     (mu->IsTrackerMuon() &&
144     mu->Quality().Quality(MuonQuality::TMLastStationTight));
145     if (pass) {
146     pt = mu->TrackerTrk()->Pt();
147     eta = TMath::Abs(mu->TrackerTrk()->Eta());
148     }
149     else {
150     pt = mu->Pt();
151     eta = TMath::Abs(mu->Eta());
152     }
153     break;
154 loizides 1.6 case kSta:
155 loizides 1.19 pass = mu->HasStandaloneTrk();
156 ceballos 1.30 if (pass) {
157     pt = mu->StandaloneTrk()->Pt();
158     eta = TMath::Abs(mu->StandaloneTrk()->Eta());
159     }
160 loizides 1.6 break;
161 loizides 1.24 case kTrackerMuon:
162 ceballos 1.28 pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
163     mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
164 ceballos 1.30 if (pass) {
165     pt = mu->TrackerTrk()->Pt();
166     eta = TMath::Abs(mu->TrackerTrk()->Eta());
167     }
168 loizides 1.24 break;
169     case kCaloMuon:
170     pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
171 ceballos 1.30 if (pass) {
172     pt = mu->TrackerTrk()->Pt();
173     eta = TMath::Abs(mu->TrackerTrk()->Eta());
174     }
175 loizides 1.24 break;
176     case kTrackerBased:
177 loizides 1.19 pass = mu->HasTrackerTrk();
178 ceballos 1.30 if (pass) {
179     pt = mu->TrackerTrk()->Pt();
180     eta = TMath::Abs(mu->TrackerTrk()->Eta());
181     }
182 loizides 1.6 break;
183 ceballos 1.84 case kGlobalOnly:
184     pass = mu->HasGlobalTrk();
185     if (pass && mu->TrackerTrk()) {
186     pt = mu->TrackerTrk()->Pt();
187     eta = TMath::Abs(mu->TrackerTrk()->Eta());
188     }
189     else {
190     pt = mu->Pt();
191     eta = TMath::Abs(mu->Eta());
192     }
193     break;
194 loizides 1.6 default:
195     break;
196 ceballos 1.5 }
197 loizides 1.6
198     if (!pass)
199     continue;
200    
201 mdecross 1.91 if (pt <= fMuonPtMin) continue;
202 loizides 1.6
203 mdecross 1.91 if (eta >= fEtaCut) continue;
204 ceballos 1.30
205 ceballos 1.34 Double_t RChi2 = 0.0;
206     if (mu->HasGlobalTrk()) {
207     RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
208     }
209     else if(mu->BestTrk() != 0){
210     RChi2 = mu->BestTrk()->Chi2()/mu->BestTrk()->Ndof();
211     }
212 loizides 1.6 Bool_t idpass = kFALSE;
213 mdecross 1.91
214    
215    
216 loizides 1.6 switch (fMuIDType) {
217 mdecross 1.91
218 ceballos 1.34 case kWMuId:
219     idpass = mu->BestTrk() != 0 &&
220     mu->BestTrk()->NHits() > 10 &&
221     RChi2 < 10.0 &&
222 ceballos 1.36 (mu->NSegments() > 1 || mu->NMatches() > 1) &&
223 ceballos 1.34 mu->BestTrk()->NPixelHits() > 0 &&
224     mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
225     break;
226     case kZMuId:
227     idpass = mu->BestTrk() != 0 &&
228     mu->BestTrk()->NHits() > 10 &&
229 ceballos 1.36 (mu->NSegments() > 1 || mu->NMatches() > 1) &&
230 ceballos 1.34 mu->BestTrk()->NPixelHits() > 0 &&
231     mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
232     break;
233 loizides 1.6 case kLoose:
234 ceballos 1.32 idpass = mu->BestTrk() != 0 &&
235     mu->Quality().Quality(MuonQuality::TMOneStationLoose) &&
236 ceballos 1.28 mu->Quality().Quality(MuonQuality::TM2DCompatibilityLoose) &&
237 ceballos 1.29 mu->BestTrk()->NHits() > 10 &&
238 ceballos 1.34 RChi2 < 10.0 &&
239 ceballos 1.30 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
240 loizides 1.6 break;
241     case kTight:
242 ceballos 1.84 idpass = mu->BestTrk() != 0 &&
243     mu->NTrkLayersHit() > 5 &&
244     mu->IsPFMuon() == kTRUE &&
245     mu->BestTrk()->NPixelHits() > 0 &&
246     RChi2 < 10.0;
247 ceballos 1.28 break;
248 veverka 1.90 case kmuonPOG2012CutBasedIDTight:
249     idpass = mu->IsGlobalMuon() &&
250     mu->IsPFMuon() &&
251     mu->GlobalTrk()->RChi2() < 10 &&
252     mu->NValidHits() != 0 &&
253     mu->NMatches() > 1 &&
254     mu->BestTrk()->NPixelHits() != 0 &&
255     mu->NTrkLayersHit() > 5;
256     break;
257 ceballos 1.83 // 2012 WW analysis for 42x (there is no PFMuon link)
258 ceballos 1.53 case kWWMuIdV1:
259 ceballos 1.32 idpass = mu->BestTrk() != 0 &&
260 ceballos 1.83 mu->NTrkLayersHit() > 5 &&
261 ceballos 1.53 mu->BestTrk()->NPixelHits() > 0 &&
262     mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
263 ceballos 1.83 mu->TrkKink() < 20.0;
264 ceballos 1.53 break;
265 ceballos 1.83 // 2010 WW analysis
266 ceballos 1.53 case kWWMuIdV2:
267     idpass = mu->BestTrk() != 0 &&
268     mu->BestTrk()->NHits() > 10 &&
269 ceballos 1.36 mu->BestTrk()->NPixelHits() > 0 &&
270     mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
271 loizides 1.6 break;
272 ceballos 1.83 // 2011 WW analysis
273 ceballos 1.57 case kWWMuIdV3:
274     idpass = mu->BestTrk() != 0 &&
275     mu->BestTrk()->NHits() > 10 &&
276     mu->BestTrk()->NPixelHits() > 0 &&
277     mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
278     mu->TrkKink() < 20.0;
279     break;
280 ceballos 1.83 // 2012 WW analysis
281 ceballos 1.77 case kWWMuIdV4:
282     idpass = mu->BestTrk() != 0 &&
283     mu->NTrkLayersHit() > 5 &&
284 ceballos 1.78 mu->IsPFMuon() == kTRUE &&
285 ceballos 1.77 mu->BestTrk()->NPixelHits() > 0 &&
286     mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
287     mu->TrkKink() < 20.0;
288     break;
289 sixie 1.59 case kMVAID_BDTG_IDIso:
290     {
291     Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
292     mu->BestTrk()->NHits() > 10 &&
293     mu->BestTrk()->NPixelHits() > 0 &&
294     mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
295     MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
296     MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
297 sixie 1.64 mu->TrkKink() < 20.0
298     );
299     idpass = passDenominatorM2;
300 sixie 1.62 //only evaluate MVA if muon passes M2 denominator to save time
301     if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
302 sixie 1.59 }
303     break;
304 loizides 1.18 case kNoId:
305 mdecross 1.91 {
306 loizides 1.18 idpass = kTRUE;
307 mdecross 1.91 }
308 loizides 1.18 break;
309 loizides 1.6 default:
310     break;
311 ceballos 1.4 }
312 loizides 1.6
313     if (!idpass)
314     continue;
315    
316 fabstoec 1.71 Double_t Rho = 0.;
317 ceballos 1.73 if(fPileupEnergyDensity){
318     const PileupEnergyDensity *rho = fPileupEnergyDensity->At(0);
319    
320     switch (fTheRhoType) {
321     case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
322     Rho = rho->RhoLowEta();
323     break;
324     case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
325     Rho = rho->Rho();
326     break;
327     case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
328     Rho = rho->RhoRandomLowEta();
329     break;
330     case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
331     Rho = rho->RhoRandom();
332     break;
333 ceballos 1.80 case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
334     Rho = rho->RhoKt6PFJets();
335     break;
336 ceballos 1.73 default:
337     Rho = rho->Rho();
338     }
339    
340     if ((TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho()))) Rho = 0.;
341 fabstoec 1.71 }
342    
343 ceballos 1.22 Bool_t isocut = kFALSE;
344 loizides 1.6 switch (fMuIsoType) {
345     case kTrackCalo:
346 ceballos 1.22 isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
347 loizides 1.6 (mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
348     break;
349     case kTrackCaloCombined:
350 ceballos 1.38 isocut = (1.0 * mu->IsoR03SumPt() +
351     1.0 * mu->IsoR03EmEt() +
352     1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
353 loizides 1.6 break;
354     case kTrackCaloSliding:
355     {
356 ceballos 1.93 Double_t totalIso = mu->IsoR03SumPt() + TMath::Max(mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3, 0.0);
357 ceballos 1.46 // trick to change the signal region cut
358     double theIsoCut = fCombIsolationCut;
359     if(theIsoCut < 0.20){
360     if(mu->Pt() > 20.0) theIsoCut = 0.15;
361     else theIsoCut = 0.10;
362     }
363     if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
364 ceballos 1.22 }
365 loizides 1.6 break;
366 fabstoec 1.81 case kTrackCaloSlidingNoCorrection:
367 ceballos 1.40 {
368     Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
369 ceballos 1.41 1.0 * mu->IsoR03EmEt() +
370     1.0 * mu->IsoR03HadEt();
371 ceballos 1.46 // trick to change the signal region cut
372     double theIsoCut = fCombIsolationCut;
373     if(theIsoCut < 0.20){
374     if(mu->Pt() > 20.0) theIsoCut = 0.15;
375     else theIsoCut = 0.10;
376     }
377     if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
378 ceballos 1.40 }
379     break;
380 fabstoec 1.81 case kCombinedRelativeConeAreaCorrected:
381     {
382     //const PileupEnergyDensity *rho = fPileupEnergyDensity->At(0); // Fabian: made Rho customable
383 ceballos 1.92 Double_t totalIso = mu->IsoR03SumPt() + TMath::Max(mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3, 0.0);
384 fabstoec 1.81 double theIsoCut = fCombRelativeIsolationCut;
385     if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
386     }
387     break;
388     case kCombinedRelativeEffectiveAreaCorrected:
389     {
390     Double_t tmpRho = Rho; // Fabian: made the Rho type customable.
391     //if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))
392     //tmpRho = fPileupEnergyDensity->At(0)->Rho();
393    
394     isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
395     - tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
396     - tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
397     ) < (mu->Pt()* 0.40);
398     }
399     break;
400     case kPFIso:
401     {
402     Double_t pfIsoCutValue = 9999;
403     if(fPFIsolationCut > 0){
404     pfIsoCutValue = fPFIsolationCut;
405     } else {
406     if (mu->AbsEta() < 1.479) {
407     if (mu->Pt() > 20) {
408     pfIsoCutValue = 0.13;
409     } else {
410     pfIsoCutValue = 0.06;
411     }
412     } else {
413     if (mu->Pt() > 20) {
414     pfIsoCutValue = 0.09;
415     } else {
416     pfIsoCutValue = 0.05;
417 ceballos 1.51 }
418 fabstoec 1.81 }
419 ceballos 1.38 }
420 ceballos 1.86 Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
421 fabstoec 1.81 if (totalIso < (mu->Pt()*pfIsoCutValue) )
422     isocut = kTRUE;
423     }
424     break;
425     case kPFRadialIso:
426     {
427     Double_t pfIsoCutValue = 9999;
428     if(fPFIsolationCut > 0){
429     pfIsoCutValue = fPFIsolationCut;
430     } else {
431     if (mu->Pt() > 20) {
432     pfIsoCutValue = 0.10;
433     } else {
434 ceballos 1.70 pfIsoCutValue = 0.05;
435     }
436     }
437     Double_t totalIso = IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
438     if (totalIso < (mu->Pt()*pfIsoCutValue) )
439     isocut = kTRUE;
440 fabstoec 1.85 }
441     break;
442     case kPFIsoBetaPUCorrected:
443     {
444     Double_t pfIsoCutValue = 9999;
445     if(fPFIsolationCut > 0){
446     pfIsoCutValue = fPFIsolationCut;
447     } else {
448     if (mu->Pt() > 20) {
449     pfIsoCutValue = 0.2;
450     } else {
451     pfIsoCutValue = 0.2;
452     }
453 ceballos 1.70 }
454 fabstoec 1.85 Double_t totalIso = IsolationTools::BetaMwithPUCorrection(fPFNoPileUpCands, fPFPileUpCands, mu, 0.4);
455 veverka 1.90
456 fabstoec 1.85 if (totalIso < (mu->Pt()*pfIsoCutValue) )
457     isocut = kTRUE;
458     }
459 fabstoec 1.81 break;
460     case kPFIsoEffectiveAreaCorrected:
461     {
462     Double_t pfIsoCutValue = 9999;
463     if(fPFIsolationCut > 0){
464     pfIsoCutValue = fPFIsolationCut;
465     } else {
466     pfIsoCutValue = fPFIsolationCut; //leave it like this for now
467     }
468 fabstoec 1.85 Double_t EffectiveAreaCorrectedPFIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
469     - Rho * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
470 fabstoec 1.81 //- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta()); // Fabian: made Rho-type customable
471     isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
472     break;
473     }
474    
475    
476     case kPFIsoNoL:
477 ceballos 1.38 {
478 ceballos 1.39 fNonIsolatedMuons = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
479     fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
480 ceballos 1.38
481 ceballos 1.58 Double_t pfIsoCutValue = 9999;
482     if(fPFIsolationCut > 0){
483     pfIsoCutValue = fPFIsolationCut;
484     } else {
485     if (mu->AbsEta() < 1.479) {
486     if (mu->Pt() > 20) {
487     pfIsoCutValue = 0.13;
488     } else {
489     pfIsoCutValue = 0.06;
490     }
491     } else {
492     if (mu->Pt() > 20) {
493     pfIsoCutValue = 0.09;
494     } else {
495     pfIsoCutValue = 0.05;
496     }
497     }
498     }
499     Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
500     if (totalIso < (mu->Pt()*pfIsoCutValue) )
501 ceballos 1.38 isocut = kTRUE;
502     }
503     break;
504 sixie 1.59 case kMVAIso_BDTG_IDIso:
505 sixie 1.64 {
506 sixie 1.65
507     Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
508     isocut = (totalIso < (mu->Pt()*0.4));
509    
510 sixie 1.64 }
511 sixie 1.59 break;
512 ceballos 1.74 case kIsoRingsV0_BDTG_Iso:
513     {
514    
515     isocut = PassMuonIsoRingsV0_BDTG_Iso(mu, fVertices->At(0), fPileupEnergyDensity);
516    
517     }
518     break;
519 ceballos 1.75 case kIsoDeltaR:
520     {
521    
522     isocut = PassMuonIsoDeltaR(mu, fVertices->At(0), fPileupEnergyDensity);
523    
524     }
525     break;
526 loizides 1.13 case kNoIso:
527 ceballos 1.22 isocut = kTRUE;
528 loizides 1.13 break;
529 loizides 1.6 case kCustomIso:
530     default:
531     break;
532 ceballos 1.4 }
533 ceballos 1.3
534 ceballos 1.22 if (isocut == kFALSE)
535 loizides 1.6 continue;
536    
537 ceballos 1.42 // apply d0 cut
538 loizides 1.21 if (fApplyD0Cut) {
539 ceballos 1.42 Bool_t passD0cut = kTRUE;
540 ceballos 1.46 if(fD0Cut < 0.05) { // trick to change the signal region cut
541     if (mu->Pt() > 20.0) fD0Cut = 0.02;
542     else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
543     }
544 ceballos 1.42 if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
545     else passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
546 ceballos 1.31 if (!passD0cut)
547 loizides 1.21 continue;
548 ceballos 1.14 }
549    
550 ceballos 1.42 // apply dz cut
551     if (fApplyDZCut) {
552 ceballos 1.45 Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
553 ceballos 1.42 if (!passDZcut)
554     continue;
555     }
556    
557 loizides 1.6 // add good muon
558     CleanMuons->Add(mu);
559 loizides 1.1 }
560    
561 loizides 1.10 // sort according to pt
562     CleanMuons->Sort();
563    
564 loizides 1.6 // add objects for other modules to use
565 loizides 1.7 AddObjThisEvt(CleanMuons);
566 loizides 1.1 }
567    
568     //--------------------------------------------------------------------------------------------------
569     void MuonIDMod::SlaveBegin()
570     {
571     // Run startup code on the computer (slave) doing the actual analysis. Here,
572 loizides 1.6 // we just request the muon collection branch.
573    
574 mingyang 1.87 if(fCleanMuonsName.CompareTo("HggLeptonTagMuons") == 0 ){
575     ReqEventObject(fPVName,fVertices,true);
576     }
577    
578 ceballos 1.38 // In this case we cannot have a branch
579     if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
580     ReqEventObject(fMuonBranchName, fMuons, kTRUE);
581     }
582 ceballos 1.42 ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
583 ceballos 1.38 ReqEventObject(fTrackName, fTracks, kTRUE);
584     ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
585 sixie 1.55 if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
586 fabstoec 1.81 || fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
587     || fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
588     || fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
589 sixie 1.59 || fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
590 ceballos 1.74 || fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0
591 ceballos 1.75 || fMuonIsoType.CompareTo("IsoDeltaR") == 0
592 mingyang 1.87 ) {
593 ceballos 1.41 ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
594     }
595 loizides 1.1
596 sixie 1.59
597 ceballos 1.34 if (fMuonIDType.CompareTo("WMuId") == 0)
598     fMuIDType = kWMuId;
599     else if (fMuonIDType.CompareTo("ZMuId") == 0)
600     fMuIDType = kZMuId;
601     else if (fMuonIDType.CompareTo("Tight") == 0)
602 loizides 1.6 fMuIDType = kTight;
603 veverka 1.90 else if (fMuonIDType.CompareTo("muonPOG2012CutBasedIDTight") == 0)
604     fMuIDType = kmuonPOG2012CutBasedIDTight;
605 loizides 1.6 else if (fMuonIDType.CompareTo("Loose") == 0)
606     fMuIDType = kLoose;
607 ceballos 1.53 else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
608     fMuIDType = kWWMuIdV1;
609     else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
610     fMuIDType = kWWMuIdV2;
611 ceballos 1.57 else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
612     fMuIDType = kWWMuIdV3;
613 ceballos 1.80 else if (fMuonIDType.CompareTo("WWMuIdV4") == 0)
614     fMuIDType = kWWMuIdV4;
615 loizides 1.12 else if (fMuonIDType.CompareTo("NoId") == 0)
616     fMuIDType = kNoId;
617 loizides 1.6 else if (fMuonIDType.CompareTo("Custom") == 0) {
618     fMuIDType = kCustomId;
619     SendError(kWarning, "SlaveBegin",
620     "Custom muon identification is not yet implemented.");
621 sixie 1.59 } else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
622     fMuIDType = kMVAID_BDTG_IDIso;
623 loizides 1.6 } else {
624     SendError(kAbortAnalysis, "SlaveBegin",
625     "The specified muon identification %s is not defined.",
626     fMuonIDType.Data());
627     return;
628     }
629 loizides 1.1
630 loizides 1.6 if (fMuonIsoType.CompareTo("TrackCalo") == 0)
631     fMuIsoType = kTrackCalo;
632     else if (fMuonIsoType.CompareTo("TrackCaloCombined") == 0)
633     fMuIsoType = kTrackCaloCombined;
634     else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
635     fMuIsoType = kTrackCaloSliding;
636 ceballos 1.41 else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
637     fMuIsoType = kTrackCaloSlidingNoCorrection;
638 fabstoec 1.81 else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
639     fMuIsoType = kCombinedRelativeConeAreaCorrected;
640     else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
641     fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
642 ceballos 1.38 else if (fMuonIsoType.CompareTo("PFIso") == 0)
643     fMuIsoType = kPFIso;
644 ceballos 1.70 else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
645     fMuIsoType = kPFRadialIso;
646 fabstoec 1.85 else if (fMuonIsoType.CompareTo("PFIsoBetaPUCorrected") == 0)
647     fMuIsoType = kPFIsoBetaPUCorrected;
648 fabstoec 1.81 else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
649     fMuIsoType = kPFIsoEffectiveAreaCorrected;
650 ceballos 1.38 else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
651     fMuIsoType = kPFIsoNoL;
652 loizides 1.6 else if (fMuonIsoType.CompareTo("NoIso") == 0)
653     fMuIsoType = kNoIso;
654     else if (fMuonIsoType.CompareTo("Custom") == 0) {
655     fMuIsoType = kCustomIso;
656     SendError(kWarning, "SlaveBegin",
657     "Custom muon isolation is not yet implemented.");
658 ceballos 1.74 } else if (fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0) {
659 sixie 1.59 fMuIsoType = kMVAIso_BDTG_IDIso;
660 ceballos 1.74 } else if (fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0) {
661     fMuIsoType = kIsoRingsV0_BDTG_Iso;
662 ceballos 1.75 } else if (fMuonIsoType.CompareTo("IsoDeltaR") == 0) {
663     fMuIsoType = kIsoDeltaR;
664 loizides 1.6 } else {
665     SendError(kAbortAnalysis, "SlaveBegin",
666     "The specified muon isolation %s is not defined.",
667     fMuonIsoType.Data());
668     return;
669     }
670 mdecross 1.91
671 loizides 1.6 if (fMuonClassType.CompareTo("All") == 0)
672     fMuClassType = kAll;
673     else if (fMuonClassType.CompareTo("Global") == 0)
674     fMuClassType = kGlobal;
675 ceballos 1.53 else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
676     fMuClassType = kGlobalTracker;
677 loizides 1.6 else if (fMuonClassType.CompareTo("Standalone") == 0)
678     fMuClassType = kSta;
679 loizides 1.24 else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
680     fMuClassType = kTrackerMuon;
681     else if (fMuonClassType.CompareTo("CaloMuon") == 0)
682     fMuClassType = kCaloMuon;
683     else if (fMuonClassType.CompareTo("TrackerBased") == 0)
684     fMuClassType = kTrackerBased;
685 ceballos 1.84 else if (fMuonClassType.CompareTo("GlobalOnly") == 0)
686     fMuClassType = kGlobalOnly;
687 mdecross 1.91 else if (fMuonClassType.CompareTo("GlobalorTracker") == 0)
688     fMuClassType = kGlobalorTracker;
689 loizides 1.6 else {
690     SendError(kAbortAnalysis, "SlaveBegin",
691     "The specified muon class %s is not defined.",
692     fMuonClassType.Data());
693     return;
694     }
695 sixie 1.59
696    
697     //If we use MVA ID, need to load MVA weights
698 ceballos 1.74 if (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
699 sixie 1.59 fMuonTools = new MuonTools();
700     fMuonIDMVA = new MuonIDMVA();
701     fMuonIDMVA->Initialize("BDTG method",
702 ceballos 1.76 string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin0_IDIsoCombined_BDTG.weights.xml"))),
703     string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin0_IDIsoCombined_BDTG.weights.xml"))),
704     string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin1_IDIsoCombined_BDTG.weights.xml"))),
705     string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin1_IDIsoCombined_BDTG.weights.xml"))),
706     string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin2_IDIsoCombined_BDTG.weights.xml"))),
707     string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin2_IDIsoCombined_BDTG.weights.xml"))),
708 ceballos 1.80 MuonIDMVA::kIDIsoCombinedDetIso,
709     fTheRhoType);
710 sixie 1.59 }
711 ceballos 1.74 else if(fMuIsoType == kIsoRingsV0_BDTG_Iso) {
712     std::vector<std::string> muonidiso_weightfiles;
713     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml"))));
714     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml"))));
715     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml"))));
716     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml"))));
717     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml"))));
718     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml"))));
719     fMuonTools = new MuonTools();
720     fMuonIDMVA = new MuonIDMVA();
721     fMuonIDMVA->Initialize("MuonIso_BDTG_IsoRings",
722     MuonIDMVA::kIsoRingsV0,
723     kTRUE,
724 ceballos 1.80 muonidiso_weightfiles,
725     fTheRhoType);
726 ceballos 1.74 }
727 ceballos 1.75 else if(fMuIsoType == kIsoDeltaR) {
728     std::vector<std::string> muonidiso_weightfiles;
729     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LB_BDT.weights.xml"))));
730     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LE_BDT.weights.xml"))));
731     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HB_BDT.weights.xml"))));
732     muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HE_BDT.weights.xml"))));
733     fMuonTools = new MuonTools();
734     fMuonIDMVA = new MuonIDMVA();
735     fMuonIDMVA->Initialize("muonHZZ2012IsoDRMVA",
736     MuonIDMVA::kIsoDeltaR,
737     kTRUE,
738 ceballos 1.80 muonidiso_weightfiles,
739     fTheRhoType);
740 ceballos 1.75 }
741 sixie 1.59
742     }
743    
744    
745     //--------------------------------------------------------------------------------------------------
746 sixie 1.61 Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
747     const PileupEnergyDensityCol *PileupEnergyDensity) const
748 sixie 1.59 {
749    
750     const Track *muTrk=0;
751     if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
752     else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
753 sixie 1.64
754 sixie 1.61 Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
755 sixie 1.59
756     Int_t subdet = 0;
757     if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
758     else subdet = 1;
759     Int_t ptBin = 0;
760     if (muTrk->Pt() > 14.5) ptBin = 1;
761     if (muTrk->Pt() > 20.0) ptBin = 2;
762    
763     Int_t MVABin = -1;
764 ceballos 1.67 if (subdet == 0 && ptBin == 0) MVABin = 0;
765     else if (subdet == 1 && ptBin == 0) MVABin = 1;
766     else if (subdet == 0 && ptBin == 1) MVABin = 2;
767     else if (subdet == 1 && ptBin == 1) MVABin = 3;
768     else if (subdet == 0 && ptBin == 2) MVABin = 4;
769     else if (subdet == 1 && ptBin == 2) MVABin = 5;
770 sixie 1.59
771     Double_t MVACut = -999;
772 ceballos 1.67 if (MVABin == 0) MVACut = -0.5618;
773     else if (MVABin == 1) MVACut = -0.3002;
774     else if (MVABin == 2) MVACut = -0.4642;
775     else if (MVABin == 3) MVACut = -0.2478;
776     else if (MVABin == 4) MVACut = 0.1706;
777     else if (MVABin == 5) MVACut = 0.8146;
778 sixie 1.59
779     if (MVAValue > MVACut) return kTRUE;
780     return kFALSE;
781 loizides 1.1 }
782 ceballos 1.74
783     //--------------------------------------------------------------------------------------------------
784     Bool_t MuonIDMod::PassMuonIsoRingsV0_BDTG_Iso(const Muon *mu, const Vertex *vertex,
785     const PileupEnergyDensityCol *PileupEnergyDensity) const
786     {
787    
788 ceballos 1.80 Bool_t isDebug = kFALSE;
789 ceballos 1.74 const Track *muTrk=0;
790     if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
791     else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
792    
793     ElectronOArr *tempElectrons = new ElectronOArr;
794     MuonOArr *tempMuons = new MuonOArr;
795     Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,
796 ceballos 1.80 PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,isDebug);
797 ceballos 1.74 delete tempElectrons;
798     delete tempMuons;
799    
800     Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
801    
802 ceballos 1.80 Double_t MVACut = -1.0;
803     Double_t eta = mu->AbsEta();
804     if (mu->Pt() < 20 && eta < 1.479) MVACut = 0.86;
805     else if(mu->Pt() < 20 && eta >= 1.479) MVACut = 0.82;
806     else if(mu->Pt() >= 20 && eta < 1.479) MVACut = 0.82;
807     else if(mu->Pt() >= 20 && eta >= 1.479) MVACut = 0.86;
808    
809 ceballos 1.82 if(fPFIsolationCut > -1.0) MVACut = fPFIsolationCut;
810    
811 ceballos 1.80 if(isDebug == kTRUE){
812     printf("PassMuonIsoRingsV0_BDTG_IsoDebug: %d, pt, eta = %f, %f, rho = %f(%f) : RingsMVA = %f, bin: %d\n",
813     GetEventHeader()->EvtNum(),mu->Pt(), mu->Eta(),
814     fPileupEnergyDensity->At(0)->Rho(),fPileupEnergyDensity->At(0)->RhoKt6PFJets(),MVAValue,MVABin);
815     }
816 ceballos 1.74
817     if (MVAValue > MVACut) return kTRUE;
818     return kFALSE;
819     }
820 ceballos 1.75
821     //--------------------------------------------------------------------------------------------------
822     Bool_t MuonIDMod::PassMuonIsoDeltaR(const Muon *mu, const Vertex *vertex,
823     const PileupEnergyDensityCol *PileupEnergyDensity) const
824     {
825    
826     const Track *muTrk=0;
827     if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
828     else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
829    
830     ElectronOArr *tempElectrons = new ElectronOArr;
831     MuonOArr *tempMuons = new MuonOArr;
832     Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFNoPileUpCands,
833 ceballos 1.80 PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,kFALSE);
834 ceballos 1.75 delete tempElectrons;
835     delete tempMuons;
836    
837     Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
838    
839     Double_t MVACut = -999;
840     if (MVABin == 0) MVACut = 0.000;
841     else if (MVABin == 1) MVACut = 0.000;
842     else if (MVABin == 2) MVACut = 0.000;
843     else if (MVABin == 3) MVACut = 0.000;
844    
845     if (MVAValue > MVACut) return kTRUE;
846     return kFALSE;
847     }
848    
849     //--------------------------------------------------------------------------------------------------
850     void MuonIDMod::Terminate()
851     {
852     // Run finishing code on the computer (slave) that did the analysis
853     delete fMuonIDMVA;
854    
855     delete fMuonTools;
856     }