ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.78
Committed: Sun May 6 10:32:30 2012 UTC (13 years ago) by ceballos
Content type: text/plain
Branch: MAIN
Changes since 1.77: +2 -2 lines
Log Message: 
adding WWMuIdV4

File Contents

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