ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc

Comparing UserCode/MitPhysics/Mods/src/MuonIDMod.cc (file contents):
Revision 1.13 by loizides, Thu Dec 11 15:53:03 2008 UTC vs.
Revision 1.60 by sixie, Wed Jan 4 10:33:27 2012 UTC

#	Line 2 | Line 2
2
3		#include "MitPhysics/Mods/interface/MuonIDMod.h"
4		#include "MitCommon/MathTools/interface/MathUtils.h"
5	+	#include "MitAna/DataTree/interface/MuonFwd.h"
6	+	#include "MitAna/DataTree/interface/ElectronFwd.h"
7	+	#include "MitAna/DataTree/interface/VertexCol.h"
8		#include "MitPhysics/Init/interface/ModNames.h"
9
10		using namespace mithep;
#	Line 13 | Line 16 | ClassImp(mithep::MuonIDMod)
16		BaseMod(name,title),
17		fMuonBranchName(Names::gkMuonBrn),
18		fCleanMuonsName(ModNames::gkCleanMuonsName),
19	<	fMuonIDType("Loose"),
20	<	fMuonIsoType("TrackCaloSliding"),
19	>	fNonIsolatedMuonsName("random"),
20	>	fNonIsolatedElectronsName("random"),
21	>	fVertexName(ModNames::gkGoodVertexesName),
22	>	fBeamSpotName(Names::gkBeamSpotBrn),
23	>	fTrackName(Names::gkTrackBrn),
24	>	fPFCandidatesName(Names::gkPFCandidatesBrn),
25	>	fMuonIDType("WWMuIdV3"),
26	>	fMuonIsoType("PFIso"),
27		fMuonClassType("Global"),
28		fTrackIsolationCut(3.0),
29		fCaloIsolationCut(3.0),
30	<	fCombIsolationCut(5.0),
30	>	fCombIsolationCut(0.15),
31	>	fCombRelativeIsolationCut(0.15),
32	>	fPFIsolationCut(-1.0),
33		fMuonPtMin(10),
34	<	fMuons(0)
34	>	fApplyD0Cut(kTRUE),
35	>	fApplyDZCut(kTRUE),
36	>	fD0Cut(0.020),
37	>	fDZCut(0.10),
38	>	fWhichVertex(-1),
39	>	fEtaCut(2.4),
40	>	fMuIDType(kIdUndef),
41	>	fMuIsoType(kIsoUndef),
42	>	fMuClassType(kClassUndef),
43	>	fMuons(0),
44	>	fVertices(0),
45	>	fBeamSpot(0),
46	>	fTracks(0),
47	>	fPFCandidates(0),
48	>	fIntRadius(0.0),
49	>	fNonIsolatedMuons(0),
50	>	fNonIsolatedElectrons(0),
51	>	fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
52	>	fPileupEnergyDensity(0),
53	>	fMuonTools(0),
54	>	fMuonIDMVA(0),
55	>	fMuonMVAWeights_Subdet0Pt10To14p5(""),
56	>	fMuonMVAWeights_Subdet1Pt10To14p5(""),
57	>	fMuonMVAWeights_Subdet0Pt14p5To20(""),
58	>	fMuonMVAWeights_Subdet1Pt14p5To20(""),
59	>	fMuonMVAWeights_Subdet0Pt20ToInf(""),
60	>	fMuonMVAWeights_Subdet1Pt20ToInf("")
61		{
62		// Constructor.
63		}
#	Line 30 | Line 67 | void MuonIDMod::Process()
67		{
68		// Process entries of the tree.
69
70	<	LoadBranch(fMuonBranchName);
70	>	if(fMuIsoType != kPFIsoNoL) {
71	>	LoadEventObject(fMuonBranchName, fMuons);
72	>	}
73	>	else {
74	>	fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
75	>	}
76	>	LoadEventObject(fBeamSpotName, fBeamSpot);
77	>	LoadEventObject(fTrackName, fTracks);
78	>	LoadEventObject(fPFCandidatesName, fPFCandidates);
79	>	if(fMuIsoType == kTrackCaloSliding ||
80	>	fMuIsoType == kCombinedRelativeConeAreaCorrected ||
81	>	fMuIsoType == kPFIsoEffectiveAreaCorrected ||
82	>	fMuIsoType == kMVAIso_BDTG_IDIso
83	>	) {
84	>	LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
85	>	}
86
87		MuonOArr *CleanMuons = new MuonOArr;
88		CleanMuons->SetName(fCleanMuonsName);
89
90	+	fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
91	+
92		for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
93		const Muon *mu = fMuons->At(i);
94
95		Bool_t pass = kFALSE;
96	<	Double_t pt = -1; // make sure pt is taken from the correct track!
96	>	Double_t pt = 0;  // make sure pt is taken from the correct track!
97	>	Double_t eta = 0; // make sure eta is taken from the correct track!
98		switch (fMuClassType) {
99		case kAll:
100		pass = kTRUE;
101	<	pt = mu->Pt();
101	>	if (mu->HasTrk()) {
102	>	pt  = mu->Pt();
103	>	eta = TMath::Abs(mu->Eta());
104	>	}
105		break;
106		case kGlobal:
107	<	pass = (mu->GlobalTrk() != 0);
108	<	if (pass)
109	<	pt = mu->GlobalTrk()->Pt();
110	<	break;
107	>	pass = mu->HasGlobalTrk() && mu->IsTrackerMuon();
108	>	if (pass && mu->TrackerTrk()) {
109	>	pt  = mu->TrackerTrk()->Pt();
110	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
111	>	}
112	>	else {
113	>	pt  = mu->Pt();
114	>	eta = TMath::Abs(mu->Eta());
115	>	}
116	>	break;
117	>	case kGlobalTracker:
118	>	pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
119	>	(mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
120	>	(mu->IsTrackerMuon() &&
121	>	mu->Quality().Quality(MuonQuality::TMLastStationTight));
122	>	if (pass) {
123	>	pt  = mu->TrackerTrk()->Pt();
124	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
125	>	}
126	>	else {
127	>	pt  = mu->Pt();
128	>	eta = TMath::Abs(mu->Eta());
129	>	}
130	>	break;
131		case kSta:
132	<	pass = (mu->StandaloneTrk() != 0);
133	<	if (pass)
134	<	pt = mu->StandaloneTrk()->Pt();
135	<	break;
136	<	case kTrackerOnly:
137	<	pass = (mu->TrackerTrk() != 0);
138	<	if (pass)
139	<	pt = mu->TrackerTrk()->Pt();
132	>	pass = mu->HasStandaloneTrk();
133	>	if (pass) {
134	>	pt  = mu->StandaloneTrk()->Pt();
135	>	eta = TMath::Abs(mu->StandaloneTrk()->Eta());
136	>	}
137	>	break;
138	>	case kTrackerMuon:
139	>	pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
140	>	mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
141	>	if (pass) {
142	>	pt  = mu->TrackerTrk()->Pt();
143	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
144	>	}
145	>	break;
146	>	case kCaloMuon:
147	>	pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
148	>	if (pass) {
149	>	pt  = mu->TrackerTrk()->Pt();
150	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
151	>	}
152	>	break;
153	>	case kTrackerBased:
154	>	pass = mu->HasTrackerTrk();
155	>	if (pass) {
156	>	pt  = mu->TrackerTrk()->Pt();
157	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
158	>	}
159		break;
160		default:
161		break;
#	Line 70 | Line 167 | void MuonIDMod::Process()
167		if (pt <= fMuonPtMin)
168		continue;
169
170	+	if (eta >= fEtaCut)
171	+	continue;
172	+
173	+	Double_t RChi2 = 0.0;
174	+	if     (mu->HasGlobalTrk()) {
175	+	RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
176	+	}
177	+	else if(mu->BestTrk() != 0){
178	+	RChi2 = mu->BestTrk()->Chi2()/mu->BestTrk()->Ndof();
179	+	}
180		Bool_t idpass = kFALSE;
181		switch (fMuIDType) {
182	+	case kWMuId:
183	+	idpass = mu->BestTrk() != 0 &&
184	+	mu->BestTrk()->NHits() > 10 &&
185	+	RChi2 < 10.0 &&
186	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
187	+	mu->BestTrk()->NPixelHits() > 0 &&
188	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
189	+	break;
190	+	case kZMuId:
191	+	idpass = mu->BestTrk() != 0 &&
192	+	mu->BestTrk()->NHits() > 10 &&
193	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
194	+	mu->BestTrk()->NPixelHits() > 0 &&
195	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
196	+	break;
197		case kLoose:
198	<	idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationLoose) &&
199	<	fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityLoose);
198	>	idpass = mu->BestTrk() != 0 &&
199	>	mu->Quality().Quality(MuonQuality::TMOneStationLoose) &&
200	>	mu->Quality().Quality(MuonQuality::TM2DCompatibilityLoose) &&
201	>	mu->BestTrk()->NHits() > 10 &&
202	>	RChi2 < 10.0 &&
203	>	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
204		break;
205		case kTight:
206	<	idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationTight) &&
207	<	fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityTight);
206	>	idpass = mu->BestTrk() !=  0 &&
207	>	mu->Quality().Quality(MuonQuality::TMOneStationTight) &&
208	>	mu->Quality().Quality(MuonQuality::TM2DCompatibilityTight) &&
209	>	mu->BestTrk()->NHits() > 10 &&
210	>	RChi2 < 10.0 &&
211	>	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
212	>	break;
213	>	case kWWMuIdV1:
214	>	idpass = mu->BestTrk() != 0 &&
215	>	mu->BestTrk()->NHits() > 10 &&
216	>	mu->BestTrk()->NPixelHits() > 0 &&
217	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
218	>	RChi2 < 10.0 &&
219	>	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
220	>	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
221	>	break;
222	>	case kWWMuIdV2:
223	>	idpass = mu->BestTrk() != 0 &&
224	>	mu->BestTrk()->NHits() > 10 &&
225	>	mu->BestTrk()->NPixelHits() > 0 &&
226	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
227	>	break;
228	>	case kWWMuIdV3:
229	>	idpass = mu->BestTrk() != 0 &&
230	>	mu->BestTrk()->NHits() > 10 &&
231	>	mu->BestTrk()->NPixelHits() > 0 &&
232	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
233	>	mu->TrkKink() < 20.0;
234	>	break;
235	>	case kMVAID_BDTG_IDIso:
236	>	{
237	>	Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
238	>	mu->BestTrk()->NHits() > 10 &&
239	>	mu->BestTrk()->NPixelHits() > 0 &&
240	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
241	>	MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
242	>	MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
243	>	mu->TrkKink() < 20.0
244	>	);
245	>	idpass =  ( passDenominatorM2 &&
246	>	PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity->At(0)->Rho()) );
247	>	}
248	>	break;
249	>	case kNoId:
250	>	idpass = kTRUE;
251		break;
252		default:
253		break;
#	Line 87 | Line 256 | void MuonIDMod::Process()
256		if (!idpass)
257		continue;
258
259	<	Bool_t isopass = kFALSE;
259	>	Bool_t isocut = kFALSE;
260		switch (fMuIsoType) {
261		case kTrackCalo:
262	<	isopass = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
262	>	isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
263		(mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
264		break;
265		case kTrackCaloCombined:
266	<	isopass = (1.0 * mu->IsoR03SumPt() + 1.0 * mu->IsoR03EmEt() +
267	<	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
266	>	isocut = (1.0 * mu->IsoR03SumPt() +
267	>	1.0 * mu->IsoR03EmEt()  +
268	>	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
269		break;
270		case kTrackCaloSliding:
271		{
272	<	Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
273	<	1.0 * mu->IsoR03EmEt() +
274	<	1.0 * mu->IsoR03HadEt();
275	<	if ((totalIso < (mu->Pt()-10.0)*5.0/15.0 && mu->Pt() <= 25) ||
276	<	(totalIso < 5.0 && mu->Pt() > 25) ||
277	<	totalIso <= 0)
278	<	isopass = kTRUE;
272	>	const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
273	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
274	>	// trick to change the signal region cut
275	>	double theIsoCut = fCombIsolationCut;
276	>	if(theIsoCut < 0.20){
277	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
278	>	else                 theIsoCut = 0.10;
279	>	}
280	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
281	>	}
282	>	break;
283	>	case kTrackCaloSlidingNoCorrection:
284	>	{
285	>	Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
286	>	1.0 * mu->IsoR03EmEt()  +
287	>	1.0 * mu->IsoR03HadEt();
288	>	// trick to change the signal region cut
289	>	double theIsoCut = fCombIsolationCut;
290	>	if(theIsoCut < 0.20){
291	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
292	>	else                 theIsoCut = 0.10;
293	>	}
294	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
295	>	}
296	>	break;
297	>	case kCombinedRelativeConeAreaCorrected:
298	>	{
299	>	const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
300	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
301	>	double theIsoCut = fCombRelativeIsolationCut;
302	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
303		}
304	+	break;
305	+	case kPFIso:
306	+	{
307	+	Double_t pfIsoCutValue = 9999;
308	+	if(fPFIsolationCut > 0){
309	+	pfIsoCutValue = fPFIsolationCut;
310	+	} else {
311	+	if (mu->AbsEta() < 1.479) {
312	+	if (mu->Pt() > 20) {
313	+	pfIsoCutValue = 0.13;
314	+	} else {
315	+	pfIsoCutValue = 0.06;
316	+	}
317	+	} else {
318	+	if (mu->Pt() > 20) {
319	+	pfIsoCutValue = 0.09;
320	+	} else {
321	+	pfIsoCutValue = 0.05;
322	+	}
323	+	}
324	+	}
325	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
326	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
327	+	isocut = kTRUE;
328	+	}
329	+	break;
330	+	case kPFIsoEffectiveAreaCorrected:
331	+	{
332	+	Double_t pfIsoCutValue = 9999;
333	+	if(fPFIsolationCut > 0){
334	+	pfIsoCutValue = fPFIsolationCut;
335	+	} else {
336	+	pfIsoCutValue = fPFIsolationCut; //leave it like this for now
337	+	}
338	+	Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
339	+	- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
340	+	isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
341	+	break;
342	+	}
343	+	case kPFIsoNoL:
344	+	{
345	+	fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
346	+	fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
347	+
348	+	Double_t pfIsoCutValue = 9999;
349	+	if(fPFIsolationCut > 0){
350	+	pfIsoCutValue = fPFIsolationCut;
351	+	} else {
352	+	if (mu->AbsEta() < 1.479) {
353	+	if (mu->Pt() > 20) {
354	+	pfIsoCutValue = 0.13;
355	+	} else {
356	+	pfIsoCutValue = 0.06;
357	+	}
358	+	} else {
359	+	if (mu->Pt() > 20) {
360	+	pfIsoCutValue = 0.09;
361	+	} else {
362	+	pfIsoCutValue = 0.05;
363	+	}
364	+	}
365	+	}
366	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
367	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
368	+	isocut = kTRUE;
369	+	}
370	+	break;
371	+	case kMVAIso_BDTG_IDIso:
372	+	isocut = ( IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
373	+	- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta())
374	+	) < (mu->Pt()* 0.40);
375		break;
376		case kNoIso:
377	<	isopass = kTRUE;
377	>	isocut = kTRUE;
378		break;
379		case kCustomIso:
380		default:
381		break;
382		}
383
384	<	if (!isopass)
384	>	if (isocut == kFALSE)
385		continue;
386
387	+	// apply d0 cut
388	+	if (fApplyD0Cut) {
389	+	Bool_t passD0cut = kTRUE;
390	+	if(fD0Cut < 0.05) { // trick to change the signal region cut
391	+	if      (mu->Pt() >  20.0) fD0Cut = 0.02;
392	+	else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
393	+	}
394	+	if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
395	+	else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
396	+	if (!passD0cut)
397	+	continue;
398	+	}
399	+
400	+	// apply dz cut
401	+	if (fApplyDZCut) {
402	+	Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
403	+	if (!passDZcut)
404	+	continue;
405	+	}
406	+
407		// add good muon
408		CleanMuons->Add(mu);
409		}
#	Line 136 | Line 421 | void MuonIDMod::SlaveBegin()
421		// Run startup code on the computer (slave) doing the actual analysis. Here,
422		// we just request the muon collection branch.
423
424	<	ReqBranch(fMuonBranchName, fMuons);
424	>	// In this case we cannot have a branch
425	>	if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
426	>	ReqEventObject(fMuonBranchName, fMuons, kTRUE);
427	>	}
428	>	ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
429	>	ReqEventObject(fTrackName, fTracks, kTRUE);
430	>	ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
431	>	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
432	>	|| fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
433	>	|| fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
434	>	|| fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
435	>	) {
436	>	ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
437	>	}
438
141	–	fMuonTools = new MuonTools;
439
440	<	if (fMuonIDType.CompareTo("Tight") == 0)
440	>	if (fMuonIDType.CompareTo("WMuId") == 0)
441	>	fMuIDType = kWMuId;
442	>	else if (fMuonIDType.CompareTo("ZMuId") == 0)
443	>	fMuIDType = kZMuId;
444	>	else if (fMuonIDType.CompareTo("Tight") == 0)
445		fMuIDType = kTight;
446		else if (fMuonIDType.CompareTo("Loose") == 0)
447		fMuIDType = kLoose;
448	+	else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
449	+	fMuIDType = kWWMuIdV1;
450	+	else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
451	+	fMuIDType = kWWMuIdV2;
452	+	else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
453	+	fMuIDType = kWWMuIdV3;
454		else if (fMuonIDType.CompareTo("NoId") == 0)
455		fMuIDType = kNoId;
456		else if (fMuonIDType.CompareTo("Custom") == 0) {
457		fMuIDType = kCustomId;
458		SendError(kWarning, "SlaveBegin",
459		"Custom muon identification is not yet implemented.");
460	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
461	+	fMuIDType = kMVAID_BDTG_IDIso;
462		} else {
463		SendError(kAbortAnalysis, "SlaveBegin",
464		"The specified muon identification %s is not defined.",
#	Line 163 | Line 472 | void MuonIDMod::SlaveBegin()
472		fMuIsoType = kTrackCaloCombined;
473		else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
474		fMuIsoType = kTrackCaloSliding;
475	+	else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
476	+	fMuIsoType = kTrackCaloSlidingNoCorrection;
477	+	else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
478	+	fMuIsoType = kCombinedRelativeConeAreaCorrected;
479	+	else if (fMuonIsoType.CompareTo("PFIso") == 0)
480	+	fMuIsoType = kPFIso;
481	+	else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
482	+	fMuIsoType = kPFIsoEffectiveAreaCorrected;
483	+	else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
484	+	fMuIsoType = kPFIsoNoL;
485		else if (fMuonIsoType.CompareTo("NoIso") == 0)
486		fMuIsoType = kNoIso;
487		else if (fMuonIsoType.CompareTo("Custom") == 0) {
488		fMuIsoType = kCustomIso;
489		SendError(kWarning, "SlaveBegin",
490		"Custom muon isolation is not yet implemented.");
491	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
492	+	fMuIsoType = kMVAIso_BDTG_IDIso;
493		} else {
494		SendError(kAbortAnalysis, "SlaveBegin",
495		"The specified muon isolation %s is not defined.",
#	Line 180 | Line 501 | void MuonIDMod::SlaveBegin()
501		fMuClassType = kAll;
502		else if (fMuonClassType.CompareTo("Global") == 0)
503		fMuClassType = kGlobal;
504	+	else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
505	+	fMuClassType = kGlobalTracker;
506		else if (fMuonClassType.CompareTo("Standalone") == 0)
507		fMuClassType = kSta;
508	<	else if (fMuonClassType.CompareTo("TrackerOnly") == 0)
509	<	fMuClassType = kTrackerOnly;
508	>	else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
509	>	fMuClassType = kTrackerMuon;
510	>	else if (fMuonClassType.CompareTo("CaloMuon") == 0)
511	>	fMuClassType = kCaloMuon;
512	>	else if (fMuonClassType.CompareTo("TrackerBased") == 0)
513	>	fMuClassType = kTrackerBased;
514		else {
515		SendError(kAbortAnalysis, "SlaveBegin",
516		"The specified muon class %s is not defined.",
517		fMuonClassType.Data());
518		return;
519		}
520	+
521	+
522	+	//If we use MVA ID, need to load MVA weights
523	+	if(fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
524	+	fMuonTools = new MuonTools();
525	+	fMuonIDMVA = new MuonIDMVA();
526	+	fMuonIDMVA->Initialize("BDTG method",
527	+	fMuonMVAWeights_Subdet0Pt10To14p5,
528	+	fMuonMVAWeights_Subdet1Pt10To14p5,
529	+	fMuonMVAWeights_Subdet0Pt14p5To20,
530	+	fMuonMVAWeights_Subdet1Pt14p5To20,
531	+	fMuonMVAWeights_Subdet0Pt20ToInf,
532	+	fMuonMVAWeights_Subdet1Pt20ToInf,
533	+	MuonIDMVA::kV8);
534	+	}
535	+
536	+	}
537	+
538	+
539	+	//--------------------------------------------------------------------------------------------------
540	+	Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex, Double_t Rho) const
541	+	{
542	+
543	+	const Track *muTrk=0;
544	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
545	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
546	+
547	+	Double_t muNchi2 = 0.0;
548	+	if(mu->HasGlobalTrk())          { muNchi2 = mu->GlobalTrk()->RChi2();     }
549	+	else if(mu->HasStandaloneTrk()) { muNchi2 = mu->StandaloneTrk()->RChi2(); }
550	+	else if(mu->HasTrackerTrk())    { muNchi2 = mu->TrackerTrk()->RChi2();    }
551	+
552	+	Double_t ChargedIso03 = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 99999, 0.3, 0.0, 0.0);
553	+	Double_t NeutralIso03_05Threshold = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.0, 0.5, 0.3, 0.0, 0.0);
554	+	Double_t ChargedIso04 = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 99999, 0.4, 0.0, 0.0);
555	+	Double_t NeutralIso04_05Threshold = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.0, 0.5, 0.4, 0.0, 0.0);
556	+
557	+	Double_t MVAValue = fMuonIDMVA->MVAValue(
558	+	muTrk->Pt(), muTrk->Eta(), muTrk->RChi2(),muNchi2,
559	+	mu->NValidHits(),
560	+	muTrk->NHits(),
561	+	muTrk->NPixelHits(),
562	+	mu->NMatches(),
563	+	muTrk->D0Corrected(*vertex),
564	+	mu->Ip3dPV(),
565	+	mu->Ip3dPVSignificance(),
566	+	mu->TrkKink(),
567	+	fMuonTools->GetSegmentCompatability(mu),
568	+	fMuonTools->GetCaloCompatability(mu, kTRUE, kTRUE),
569	+	(mu->HadEnergy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadEnergy,muTrk->Eta()))/muTrk->Pt(),
570	+	(mu->HoEnergy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHoEnergy,muTrk->Eta()))/muTrk->Pt(),
571	+	(mu->EmEnergy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEmEnergy,muTrk->Eta()))/muTrk->Pt(),
572	+	(mu->HadS9Energy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadS9Energy,muTrk->Eta()))/muTrk->Pt(),
573	+	(mu->HoS9Energy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHoS9Energy,muTrk->Eta()))/muTrk->Pt(),
574	+	(mu->EmS9Energy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEmS9Energy,muTrk->Eta()))/muTrk->Pt(),
575	+	(ChargedIso03 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuChargedIso03,muTrk->Eta()))/muTrk->Pt(),
576	+	(NeutralIso03_05Threshold - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03,muTrk->Eta()))/muTrk->Pt(),
577	+	(ChargedIso04 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuChargedIso04,muTrk->Eta()))/muTrk->Pt(),
578	+	(NeutralIso04_05Threshold - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso04,muTrk->Eta()))/muTrk->Pt()
579	+	);
580	+
581	+	Int_t subdet = 0;
582	+	if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
583	+	else subdet = 1;
584	+	Int_t ptBin = 0;
585	+	if (muTrk->Pt() > 14.5) ptBin = 1;
586	+	if (muTrk->Pt() > 20.0) ptBin = 2;
587	+
588	+	Int_t MVABin = -1;
589	+	if (subdet == 0 && ptBin == 0) MVABin = 0;
590	+	if (subdet == 1 && ptBin == 0) MVABin = 1;
591	+	if (subdet == 0 && ptBin == 1) MVABin = 2;
592	+	if (subdet == 1 && ptBin == 1) MVABin = 3;
593	+	if (subdet == 0 && ptBin == 2) MVABin = 4;
594	+	if (subdet == 1 && ptBin == 2) MVABin = 5;
595	+
596	+	Double_t MVACut = -999;
597	+	if (MVABin == 0) MVACut = -0.377;
598	+	if (MVABin == 1) MVACut = -0.0902;
599	+	if (MVABin == 2) MVACut = -0.221;
600	+	if (MVABin == 3) MVACut = -0.0154;
601	+	if (MVABin == 4) MVACut = 0.459;
602	+	if (MVABin == 5) MVACut = 0.9158;
603	+
604	+	if (MVAValue > MVACut) return kTRUE;
605	+	return kFALSE;
606		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines