ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.73
Committed: Thu May 3 10:22:10 2012 UTC (13 years ago) by ceballos
Content type: text/plain
Branch: MAIN
Changes since 1.72: +22 -21 lines
Log Message: 
fix a critical fix to avoid crashes

File Contents

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