ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.66
Committed: Fri Jan 27 11:48:26 2012 UTC (13 years, 3 months ago) by sixie
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025d
Changes since 1.65: +18 -1 lines
Log Message: 
add option flag to print out MVA debug information

File Contents

# User Rev Content
1 sixie 1.66 // $Id: MuonIDMod.cc,v 1.65 2012/01/26 15:29:13 sixie 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 ceballos 1.57 fMuonIDType("WWMuIdV3"),
27 ceballos 1.49 fMuonIsoType("PFIso"),
28 loizides 1.6 fMuonClassType("Global"),
29 ceballos 1.2 fTrackIsolationCut(3.0),
30     fCaloIsolationCut(3.0),
31 sixie 1.54 fCombIsolationCut(0.15),
32     fCombRelativeIsolationCut(0.15),
33     fPFIsolationCut(-1.0),
34 ceballos 1.3 fMuonPtMin(10),
35 loizides 1.21 fApplyD0Cut(kTRUE),
36 ceballos 1.42 fApplyDZCut(kTRUE),
37 ceballos 1.30 fD0Cut(0.020),
38 ceballos 1.50 fDZCut(0.10),
39 ceballos 1.42 fWhichVertex(-1),
40 ceballos 1.30 fEtaCut(2.4),
41 loizides 1.15 fMuIDType(kIdUndef),
42     fMuIsoType(kIsoUndef),
43     fMuClassType(kClassUndef),
44 ceballos 1.14 fMuons(0),
45 loizides 1.15 fVertices(0),
46 ceballos 1.42 fBeamSpot(0),
47 ceballos 1.38 fTracks(0),
48     fPFCandidates(0),
49 ceballos 1.56 fIntRadius(0.0),
50 ceballos 1.39 fNonIsolatedMuons(0),
51 ceballos 1.41 fNonIsolatedElectrons(0),
52     fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
53 sixie 1.59 fPileupEnergyDensity(0),
54     fMuonTools(0),
55     fMuonIDMVA(0),
56     fMuonMVAWeights_Subdet0Pt10To14p5(""),
57     fMuonMVAWeights_Subdet1Pt10To14p5(""),
58     fMuonMVAWeights_Subdet0Pt14p5To20(""),
59     fMuonMVAWeights_Subdet1Pt14p5To20(""),
60     fMuonMVAWeights_Subdet0Pt20ToInf(""),
61     fMuonMVAWeights_Subdet1Pt20ToInf("")
62 loizides 1.1 {
63     // Constructor.
64     }
65    
66     //--------------------------------------------------------------------------------------------------
67     void MuonIDMod::Process()
68     {
69     // Process entries of the tree.
70    
71 ceballos 1.38 if(fMuIsoType != kPFIsoNoL) {
72     LoadEventObject(fMuonBranchName, fMuons);
73     }
74     else {
75     fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
76     }
77 ceballos 1.42 LoadEventObject(fBeamSpotName, fBeamSpot);
78 ceballos 1.38 LoadEventObject(fTrackName, fTracks);
79     LoadEventObject(fPFCandidatesName, fPFCandidates);
80 sixie 1.59 if(fMuIsoType == kTrackCaloSliding ||
81     fMuIsoType == kCombinedRelativeConeAreaCorrected ||
82     fMuIsoType == kPFIsoEffectiveAreaCorrected ||
83     fMuIsoType == kMVAIso_BDTG_IDIso
84     ) {
85 ceballos 1.41 LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
86     }
87 sixie 1.59
88 loizides 1.7 MuonOArr *CleanMuons = new MuonOArr;
89     CleanMuons->SetName(fCleanMuonsName);
90 loizides 1.1
91 ceballos 1.38 fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
92    
93 loizides 1.1 for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
94 loizides 1.8 const Muon *mu = fMuons->At(i);
95 loizides 1.6
96     Bool_t pass = kFALSE;
97 ceballos 1.30 Double_t pt = 0; // make sure pt is taken from the correct track!
98     Double_t eta = 0; // make sure eta is taken from the correct track!
99 loizides 1.6 switch (fMuClassType) {
100     case kAll:
101     pass = kTRUE;
102 ceballos 1.30 if (mu->HasTrk()) {
103     pt = mu->Pt();
104     eta = TMath::Abs(mu->Eta());
105     }
106 loizides 1.6 break;
107     case kGlobal:
108 ceballos 1.31 pass = mu->HasGlobalTrk() && mu->IsTrackerMuon();
109 ceballos 1.32 if (pass && mu->TrackerTrk()) {
110 ceballos 1.30 pt = mu->TrackerTrk()->Pt();
111     eta = TMath::Abs(mu->TrackerTrk()->Eta());
112     }
113 ceballos 1.32 else {
114     pt = mu->Pt();
115     eta = TMath::Abs(mu->Eta());
116     }
117 ceballos 1.30 break;
118 ceballos 1.53 case kGlobalTracker:
119     pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
120     (mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
121     (mu->IsTrackerMuon() &&
122     mu->Quality().Quality(MuonQuality::TMLastStationTight));
123     if (pass) {
124     pt = mu->TrackerTrk()->Pt();
125     eta = TMath::Abs(mu->TrackerTrk()->Eta());
126     }
127     else {
128     pt = mu->Pt();
129     eta = TMath::Abs(mu->Eta());
130     }
131     break;
132 loizides 1.6 case kSta:
133 loizides 1.19 pass = mu->HasStandaloneTrk();
134 ceballos 1.30 if (pass) {
135     pt = mu->StandaloneTrk()->Pt();
136     eta = TMath::Abs(mu->StandaloneTrk()->Eta());
137     }
138 loizides 1.6 break;
139 loizides 1.24 case kTrackerMuon:
140 ceballos 1.28 pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
141     mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
142 ceballos 1.30 if (pass) {
143     pt = mu->TrackerTrk()->Pt();
144     eta = TMath::Abs(mu->TrackerTrk()->Eta());
145     }
146 loizides 1.24 break;
147     case kCaloMuon:
148     pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
149 ceballos 1.30 if (pass) {
150     pt = mu->TrackerTrk()->Pt();
151     eta = TMath::Abs(mu->TrackerTrk()->Eta());
152     }
153 loizides 1.24 break;
154     case kTrackerBased:
155 loizides 1.19 pass = mu->HasTrackerTrk();
156 ceballos 1.30 if (pass) {
157     pt = mu->TrackerTrk()->Pt();
158     eta = TMath::Abs(mu->TrackerTrk()->Eta());
159     }
160 loizides 1.6 break;
161     default:
162     break;
163 ceballos 1.5 }
164 loizides 1.6
165     if (!pass)
166     continue;
167    
168     if (pt <= fMuonPtMin)
169     continue;
170    
171 ceballos 1.30 if (eta >= fEtaCut)
172     continue;
173    
174 sixie 1.66
175     //***********************************************************************************************
176     //Debug Info For Lepton MVA
177     //***********************************************************************************************
178     if( fPrintMVADebugInfo &&
179     (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso)
180     ) {
181     cout << "Event: " << GetEventHeader()->RunNum() << " " << GetEventHeader()->LumiSec() << " "
182     << GetEventHeader()->EvtNum() << " : Rho = " << fPileupEnergyDensity->At(0)->Rho()
183     << " : Muon " << i << " "
184     << endl;
185     fMuonIDMVA->MVAValue(mu,fVertices->At(0),fMuonTools,fPFCandidates,fPileupEnergyDensity,kTRUE);
186     }
187     //***********************************************************************************************
188    
189    
190 ceballos 1.34 Double_t RChi2 = 0.0;
191     if (mu->HasGlobalTrk()) {
192     RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
193     }
194     else if(mu->BestTrk() != 0){
195     RChi2 = mu->BestTrk()->Chi2()/mu->BestTrk()->Ndof();
196     }
197 loizides 1.6 Bool_t idpass = kFALSE;
198     switch (fMuIDType) {
199 ceballos 1.34 case kWMuId:
200     idpass = mu->BestTrk() != 0 &&
201     mu->BestTrk()->NHits() > 10 &&
202     RChi2 < 10.0 &&
203 ceballos 1.36 (mu->NSegments() > 1 || mu->NMatches() > 1) &&
204 ceballos 1.34 mu->BestTrk()->NPixelHits() > 0 &&
205     mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
206     break;
207     case kZMuId:
208     idpass = mu->BestTrk() != 0 &&
209     mu->BestTrk()->NHits() > 10 &&
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 loizides 1.6 case kLoose:
215 ceballos 1.32 idpass = mu->BestTrk() != 0 &&
216     mu->Quality().Quality(MuonQuality::TMOneStationLoose) &&
217 ceballos 1.28 mu->Quality().Quality(MuonQuality::TM2DCompatibilityLoose) &&
218 ceballos 1.29 mu->BestTrk()->NHits() > 10 &&
219 ceballos 1.34 RChi2 < 10.0 &&
220 ceballos 1.30 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
221 loizides 1.6 break;
222     case kTight:
223 ceballos 1.32 idpass = mu->BestTrk() != 0 &&
224     mu->Quality().Quality(MuonQuality::TMOneStationTight) &&
225 ceballos 1.28 mu->Quality().Quality(MuonQuality::TM2DCompatibilityTight) &&
226 ceballos 1.29 mu->BestTrk()->NHits() > 10 &&
227 ceballos 1.34 RChi2 < 10.0 &&
228 ceballos 1.30 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
229 ceballos 1.28 break;
230 ceballos 1.53 case kWWMuIdV1:
231 ceballos 1.32 idpass = mu->BestTrk() != 0 &&
232     mu->BestTrk()->NHits() > 10 &&
233 ceballos 1.53 mu->BestTrk()->NPixelHits() > 0 &&
234     mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
235 ceballos 1.34 RChi2 < 10.0 &&
236 ceballos 1.36 (mu->NSegments() > 1 || mu->NMatches() > 1) &&
237 ceballos 1.53 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
238     break;
239     case kWWMuIdV2:
240     idpass = mu->BestTrk() != 0 &&
241     mu->BestTrk()->NHits() > 10 &&
242 ceballos 1.36 mu->BestTrk()->NPixelHits() > 0 &&
243     mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
244 loizides 1.6 break;
245 ceballos 1.57 case kWWMuIdV3:
246     idpass = mu->BestTrk() != 0 &&
247     mu->BestTrk()->NHits() > 10 &&
248     mu->BestTrk()->NPixelHits() > 0 &&
249     mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
250     mu->TrkKink() < 20.0;
251     break;
252 sixie 1.59 case kMVAID_BDTG_IDIso:
253     {
254     Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
255     mu->BestTrk()->NHits() > 10 &&
256     mu->BestTrk()->NPixelHits() > 0 &&
257     mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
258     MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
259     MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
260 sixie 1.64 mu->TrkKink() < 20.0
261     );
262     idpass = passDenominatorM2;
263 sixie 1.62 //only evaluate MVA if muon passes M2 denominator to save time
264     if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
265 sixie 1.59 }
266     break;
267 loizides 1.18 case kNoId:
268     idpass = kTRUE;
269     break;
270 loizides 1.6 default:
271     break;
272 ceballos 1.4 }
273 loizides 1.6
274     if (!idpass)
275     continue;
276    
277 ceballos 1.22 Bool_t isocut = kFALSE;
278 loizides 1.6 switch (fMuIsoType) {
279     case kTrackCalo:
280 ceballos 1.22 isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
281 loizides 1.6 (mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
282     break;
283     case kTrackCaloCombined:
284 ceballos 1.38 isocut = (1.0 * mu->IsoR03SumPt() +
285     1.0 * mu->IsoR03EmEt() +
286     1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
287 loizides 1.6 break;
288     case kTrackCaloSliding:
289     {
290 ceballos 1.41 const PileupEnergyDensity *rho = fPileupEnergyDensity->At(0);
291 sixie 1.54 Double_t totalIso = mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
292 ceballos 1.46 // trick to change the signal region cut
293     double theIsoCut = fCombIsolationCut;
294     if(theIsoCut < 0.20){
295     if(mu->Pt() > 20.0) theIsoCut = 0.15;
296     else theIsoCut = 0.10;
297     }
298     if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
299 ceballos 1.22 }
300 loizides 1.6 break;
301 ceballos 1.41 case kTrackCaloSlidingNoCorrection:
302 ceballos 1.40 {
303     Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
304 ceballos 1.41 1.0 * mu->IsoR03EmEt() +
305     1.0 * mu->IsoR03HadEt();
306 ceballos 1.46 // trick to change the signal region cut
307     double theIsoCut = fCombIsolationCut;
308     if(theIsoCut < 0.20){
309     if(mu->Pt() > 20.0) theIsoCut = 0.15;
310     else theIsoCut = 0.10;
311     }
312     if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
313 ceballos 1.40 }
314     break;
315 sixie 1.54 case kCombinedRelativeConeAreaCorrected:
316     {
317     const PileupEnergyDensity *rho = fPileupEnergyDensity->At(0);
318     Double_t totalIso = mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
319     double theIsoCut = fCombRelativeIsolationCut;
320     if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
321     }
322     break;
323 sixie 1.64 case kCombinedRelativeEffectiveAreaCorrected:
324     {
325     Double_t tmpRho = 0;
326     if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || isinf(fPileupEnergyDensity->At(0)->Rho())))
327     tmpRho = fPileupEnergyDensity->At(0)->Rho();
328    
329     isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
330     - tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
331     - tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
332     ) < (mu->Pt()* 0.40);
333     }
334     break;
335 ceballos 1.38 case kPFIso:
336 sixie 1.59 {
337 ceballos 1.49 Double_t pfIsoCutValue = 9999;
338 sixie 1.54 if(fPFIsolationCut > 0){
339     pfIsoCutValue = fPFIsolationCut;
340 ceballos 1.49 } else {
341 ceballos 1.51 if (mu->AbsEta() < 1.479) {
342     if (mu->Pt() > 20) {
343 ceballos 1.52 pfIsoCutValue = 0.13;
344 ceballos 1.51 } else {
345 ceballos 1.52 pfIsoCutValue = 0.06;
346 ceballos 1.51 }
347 ceballos 1.49 } else {
348 ceballos 1.51 if (mu->Pt() > 20) {
349 ceballos 1.52 pfIsoCutValue = 0.09;
350 ceballos 1.51 } else {
351 ceballos 1.52 pfIsoCutValue = 0.05;
352 ceballos 1.51 }
353     }
354 ceballos 1.49 }
355 ceballos 1.56 Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
356 ceballos 1.49 if (totalIso < (mu->Pt()*pfIsoCutValue) )
357 ceballos 1.38 isocut = kTRUE;
358     }
359     break;
360 sixie 1.59 case kPFIsoEffectiveAreaCorrected:
361     {
362     Double_t pfIsoCutValue = 9999;
363     if(fPFIsolationCut > 0){
364     pfIsoCutValue = fPFIsolationCut;
365     } else {
366     pfIsoCutValue = fPFIsolationCut; //leave it like this for now
367     }
368     Double_t EffectiveAreaCorrectedPFIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
369     - fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
370     isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
371     break;
372     }
373 ceballos 1.38 case kPFIsoNoL:
374     {
375 ceballos 1.39 fNonIsolatedMuons = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
376     fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
377 ceballos 1.38
378 ceballos 1.58 Double_t pfIsoCutValue = 9999;
379     if(fPFIsolationCut > 0){
380     pfIsoCutValue = fPFIsolationCut;
381     } else {
382     if (mu->AbsEta() < 1.479) {
383     if (mu->Pt() > 20) {
384     pfIsoCutValue = 0.13;
385     } else {
386     pfIsoCutValue = 0.06;
387     }
388     } else {
389     if (mu->Pt() > 20) {
390     pfIsoCutValue = 0.09;
391     } else {
392     pfIsoCutValue = 0.05;
393     }
394     }
395     }
396     Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
397     if (totalIso < (mu->Pt()*pfIsoCutValue) )
398 ceballos 1.38 isocut = kTRUE;
399     }
400     break;
401 sixie 1.59 case kMVAIso_BDTG_IDIso:
402 sixie 1.64 {
403 sixie 1.65
404     // **************************************************************************
405     // Don't use effective area correction denominator. Instead use the old one.
406     // **************************************************************************
407    
408     // Double_t tmpRho = 0;
409     // if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || isinf(fPileupEnergyDensity->At(0)->Rho())))
410     // tmpRho = fPileupEnergyDensity->At(0)->Rho();
411    
412     // isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
413     // - tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
414     // - tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
415     // ) < (mu->Pt()* 0.40);
416    
417     Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
418     isocut = (totalIso < (mu->Pt()*0.4));
419    
420 sixie 1.64 }
421 sixie 1.59 break;
422 loizides 1.13 case kNoIso:
423 ceballos 1.22 isocut = kTRUE;
424 loizides 1.13 break;
425 loizides 1.6 case kCustomIso:
426     default:
427     break;
428 ceballos 1.4 }
429 ceballos 1.3
430 ceballos 1.22 if (isocut == kFALSE)
431 loizides 1.6 continue;
432    
433 ceballos 1.42 // apply d0 cut
434 loizides 1.21 if (fApplyD0Cut) {
435 ceballos 1.42 Bool_t passD0cut = kTRUE;
436 ceballos 1.46 if(fD0Cut < 0.05) { // trick to change the signal region cut
437     if (mu->Pt() > 20.0) fD0Cut = 0.02;
438     else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
439     }
440 ceballos 1.42 if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
441     else passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
442 ceballos 1.31 if (!passD0cut)
443 loizides 1.21 continue;
444 ceballos 1.14 }
445    
446 ceballos 1.42 // apply dz cut
447     if (fApplyDZCut) {
448 ceballos 1.45 Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
449 ceballos 1.42 if (!passDZcut)
450     continue;
451     }
452    
453 loizides 1.6 // add good muon
454     CleanMuons->Add(mu);
455 loizides 1.1 }
456    
457 loizides 1.10 // sort according to pt
458     CleanMuons->Sort();
459    
460 loizides 1.6 // add objects for other modules to use
461 loizides 1.7 AddObjThisEvt(CleanMuons);
462 loizides 1.1 }
463    
464     //--------------------------------------------------------------------------------------------------
465     void MuonIDMod::SlaveBegin()
466     {
467     // Run startup code on the computer (slave) doing the actual analysis. Here,
468 loizides 1.6 // we just request the muon collection branch.
469    
470 ceballos 1.38 // In this case we cannot have a branch
471     if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
472     ReqEventObject(fMuonBranchName, fMuons, kTRUE);
473     }
474 ceballos 1.42 ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
475 ceballos 1.38 ReqEventObject(fTrackName, fTracks, kTRUE);
476     ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
477 sixie 1.55 if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
478 sixie 1.59 || fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
479 sixie 1.64 || fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
480     || fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
481 sixie 1.59 || fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
482     ) {
483 ceballos 1.41 ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
484     }
485 loizides 1.1
486 sixie 1.59
487 ceballos 1.34 if (fMuonIDType.CompareTo("WMuId") == 0)
488     fMuIDType = kWMuId;
489     else if (fMuonIDType.CompareTo("ZMuId") == 0)
490     fMuIDType = kZMuId;
491     else if (fMuonIDType.CompareTo("Tight") == 0)
492 loizides 1.6 fMuIDType = kTight;
493     else if (fMuonIDType.CompareTo("Loose") == 0)
494     fMuIDType = kLoose;
495 ceballos 1.53 else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
496     fMuIDType = kWWMuIdV1;
497     else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
498     fMuIDType = kWWMuIdV2;
499 ceballos 1.57 else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
500     fMuIDType = kWWMuIdV3;
501 loizides 1.12 else if (fMuonIDType.CompareTo("NoId") == 0)
502     fMuIDType = kNoId;
503 loizides 1.6 else if (fMuonIDType.CompareTo("Custom") == 0) {
504     fMuIDType = kCustomId;
505     SendError(kWarning, "SlaveBegin",
506     "Custom muon identification is not yet implemented.");
507 sixie 1.59 } else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
508     fMuIDType = kMVAID_BDTG_IDIso;
509 loizides 1.6 } else {
510     SendError(kAbortAnalysis, "SlaveBegin",
511     "The specified muon identification %s is not defined.",
512     fMuonIDType.Data());
513     return;
514     }
515 loizides 1.1
516 loizides 1.6 if (fMuonIsoType.CompareTo("TrackCalo") == 0)
517     fMuIsoType = kTrackCalo;
518     else if (fMuonIsoType.CompareTo("TrackCaloCombined") == 0)
519     fMuIsoType = kTrackCaloCombined;
520     else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
521     fMuIsoType = kTrackCaloSliding;
522 ceballos 1.41 else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
523     fMuIsoType = kTrackCaloSlidingNoCorrection;
524 sixie 1.54 else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
525     fMuIsoType = kCombinedRelativeConeAreaCorrected;
526 sixie 1.64 else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
527     fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
528 ceballos 1.38 else if (fMuonIsoType.CompareTo("PFIso") == 0)
529     fMuIsoType = kPFIso;
530 sixie 1.59 else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
531     fMuIsoType = kPFIsoEffectiveAreaCorrected;
532 ceballos 1.38 else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
533     fMuIsoType = kPFIsoNoL;
534 loizides 1.6 else if (fMuonIsoType.CompareTo("NoIso") == 0)
535     fMuIsoType = kNoIso;
536     else if (fMuonIsoType.CompareTo("Custom") == 0) {
537     fMuIsoType = kCustomIso;
538     SendError(kWarning, "SlaveBegin",
539     "Custom muon isolation is not yet implemented.");
540 sixie 1.59 } else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
541     fMuIsoType = kMVAIso_BDTG_IDIso;
542 loizides 1.6 } else {
543     SendError(kAbortAnalysis, "SlaveBegin",
544     "The specified muon isolation %s is not defined.",
545     fMuonIsoType.Data());
546     return;
547     }
548 loizides 1.1
549 loizides 1.6 if (fMuonClassType.CompareTo("All") == 0)
550     fMuClassType = kAll;
551     else if (fMuonClassType.CompareTo("Global") == 0)
552     fMuClassType = kGlobal;
553 ceballos 1.53 else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
554     fMuClassType = kGlobalTracker;
555 loizides 1.6 else if (fMuonClassType.CompareTo("Standalone") == 0)
556     fMuClassType = kSta;
557 loizides 1.24 else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
558     fMuClassType = kTrackerMuon;
559     else if (fMuonClassType.CompareTo("CaloMuon") == 0)
560     fMuClassType = kCaloMuon;
561     else if (fMuonClassType.CompareTo("TrackerBased") == 0)
562     fMuClassType = kTrackerBased;
563 loizides 1.6 else {
564     SendError(kAbortAnalysis, "SlaveBegin",
565     "The specified muon class %s is not defined.",
566     fMuonClassType.Data());
567     return;
568     }
569 sixie 1.59
570    
571     //If we use MVA ID, need to load MVA weights
572     if(fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
573     fMuonTools = new MuonTools();
574     fMuonIDMVA = new MuonIDMVA();
575     fMuonIDMVA->Initialize("BDTG method",
576     fMuonMVAWeights_Subdet0Pt10To14p5,
577     fMuonMVAWeights_Subdet1Pt10To14p5,
578     fMuonMVAWeights_Subdet0Pt14p5To20,
579     fMuonMVAWeights_Subdet1Pt14p5To20,
580     fMuonMVAWeights_Subdet0Pt20ToInf,
581     fMuonMVAWeights_Subdet1Pt20ToInf,
582 sixie 1.64 MuonIDMVA::kIDIsoCombinedDetIso);
583 sixie 1.59 }
584    
585     }
586    
587    
588     //--------------------------------------------------------------------------------------------------
589 sixie 1.61 Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
590     const PileupEnergyDensityCol *PileupEnergyDensity) const
591 sixie 1.59 {
592    
593     const Track *muTrk=0;
594     if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
595     else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
596 sixie 1.64
597 sixie 1.61 Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
598 sixie 1.59
599     Int_t subdet = 0;
600     if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
601     else subdet = 1;
602     Int_t ptBin = 0;
603     if (muTrk->Pt() > 14.5) ptBin = 1;
604     if (muTrk->Pt() > 20.0) ptBin = 2;
605    
606     Int_t MVABin = -1;
607     if (subdet == 0 && ptBin == 0) MVABin = 0;
608     if (subdet == 1 && ptBin == 0) MVABin = 1;
609     if (subdet == 0 && ptBin == 1) MVABin = 2;
610     if (subdet == 1 && ptBin == 1) MVABin = 3;
611     if (subdet == 0 && ptBin == 2) MVABin = 4;
612     if (subdet == 1 && ptBin == 2) MVABin = 5;
613    
614     Double_t MVACut = -999;
615 sixie 1.64 if (MVABin == 0) MVACut = -0.5618;
616     if (MVABin == 1) MVACut = -0.3002;
617     if (MVABin == 2) MVACut = -0.4642;
618     if (MVABin == 3) MVACut = -0.2478;
619     if (MVABin == 4) MVACut = 0.1706;
620     if (MVABin == 5) MVACut = 0.8146;
621 sixie 1.59
622     if (MVAValue > MVACut) return kTRUE;
623     return kFALSE;
624 loizides 1.1 }