ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.88
Committed: Sat Feb 23 14:51:45 2013 UTC (12 years, 2 months ago) by mingyang
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_029c, Mit_029b, Mit_029a
Changes since 1.87: +13 -1 lines
Log Message: 
*** empty log message ***

File Contents

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