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.41 by ceballos, Fri Mar 11 15:13:09 2011 UTC vs.
Revision 1.84 by ceballos, Fri Jun 15 11:58:54 2012 UTC

#	Line 14 | Line 14 | ClassImp(mithep::MuonIDMod)
14		//--------------------------------------------------------------------------------------------------
15		MuonIDMod::MuonIDMod(const char *name, const char *title) :
16		BaseMod(name,title),
17	+	fPrintMVADebugInfo(kFALSE),
18		fMuonBranchName(Names::gkMuonBrn),
19		fCleanMuonsName(ModNames::gkCleanMuonsName),
20		fNonIsolatedMuonsName("random"),
21		fNonIsolatedElectronsName("random"),
22		fVertexName(ModNames::gkGoodVertexesName),
23	+	fBeamSpotName(Names::gkBeamSpotBrn),
24		fTrackName(Names::gkTrackBrn),
25		fPFCandidatesName(Names::gkPFCandidatesBrn),
26	<	fMuonIDType("WWMuId"),
27	<	fMuonIsoType("TrackCaloSliding"),
26	>	fMuonIDType("WWMuIdV3"),
27	>	fMuonIsoType("PFIso"),
28		fMuonClassType("Global"),
29		fTrackIsolationCut(3.0),
30		fCaloIsolationCut(3.0),
31		fCombIsolationCut(0.15),
32	+	fCombRelativeIsolationCut(0.15),
33	+	fPFIsolationCut(-999.0),
34		fMuonPtMin(10),
35		fApplyD0Cut(kTRUE),
36	+	fApplyDZCut(kTRUE),
37		fD0Cut(0.020),
38	+	fDZCut(0.10),
39	+	fWhichVertex(-1),
40		fEtaCut(2.4),
34	–	fReverseIsoCut(kFALSE),
35	–	fReverseD0Cut(kFALSE),
41		fMuIDType(kIdUndef),
42		fMuIsoType(kIsoUndef),
43		fMuClassType(kClassUndef),
44		fMuons(0),
45		fVertices(0),
46	+	fBeamSpot(0),
47		fTracks(0),
48		fPFCandidates(0),
49	+	fPFNoPileUpCands(0),
50	+	fIntRadius(0.0),
51		fNonIsolatedMuons(0),
52		fNonIsolatedElectrons(0),
53		fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
54	<	fPileupEnergyDensity(0)
54	>	fPileupEnergyDensity(0),
55	>	fMuonTools(0),
56	>	fMuonIDMVA(0),
57	>	fTheRhoType(RhoUtilities::DEFAULT)
58		{
59		// Constructor.
60		}
#	Line 59 | Line 70 | void MuonIDMod::Process()
70		else {
71		fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
72		}
73	+	LoadEventObject(fBeamSpotName, fBeamSpot);
74		LoadEventObject(fTrackName, fTracks);
75		LoadEventObject(fPFCandidatesName, fPFCandidates);
76	<	if(fMuIsoType == kTrackCaloSliding) {
76	>	if(fMuIsoType == kTrackCaloSliding ||
77	>	fMuIsoType == kCombinedRelativeConeAreaCorrected ||
78	>	fMuIsoType == kPFIsoEffectiveAreaCorrected ||
79	>	fMuIsoType == kMVAIso_BDTG_IDIso ||
80	>	fMuIsoType == kIsoRingsV0_BDTG_Iso ||
81	>	fMuIsoType == kIsoDeltaR
82	>	) {
83		LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
84		}
85	+	if(fMuIsoType == kPFRadialIso || fMuIsoType == kIsoDeltaR){
86	+	// Name is hardcoded, can be changed if someone feels to do it
87	+	fPFNoPileUpCands = GetObjThisEvt<PFCandidateCol>("PFNoPileUp");
88	+	}
89	+
90		MuonOArr *CleanMuons = new MuonOArr;
91		CleanMuons->SetName(fCleanMuonsName);
92
93		fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
94
95	<	for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
95	>	for (UInt_t i=0; i<fMuons->GetEntries() && fVertices->GetEntries() > 0 ; ++i) {
96		const Muon *mu = fMuons->At(i);
97
98		Bool_t pass = kFALSE;
#	Line 94 | Line 117 | void MuonIDMod::Process()
117		eta = TMath::Abs(mu->Eta());
118		}
119		break;
120	+	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		case kSta:
135		pass = mu->HasStandaloneTrk();
136		if (pass) {
#	Line 123 | Line 160 | void MuonIDMod::Process()
160		eta = TMath::Abs(mu->TrackerTrk()->Eta());
161		}
162		break;
163	+	case kGlobalOnly:
164	+	pass = mu->HasGlobalTrk();
165	+	if (pass && mu->TrackerTrk()) {
166	+	pt  = mu->TrackerTrk()->Pt();
167	+	eta = TMath::Abs(mu->TrackerTrk()->Eta());
168	+	}
169	+	else {
170	+	pt  = mu->Pt();
171	+	eta = TMath::Abs(mu->Eta());
172	+	}
173	+	break;
174		default:
175		break;
176		}
#	Line 169 | Line 217 | void MuonIDMod::Process()
217		mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
218		break;
219		case kTight:
220	<	idpass = mu->BestTrk() !=  0 &&
221	<	mu->Quality().Quality(MuonQuality::TMOneStationTight) &&
222	<	mu->Quality().Quality(MuonQuality::TM2DCompatibilityTight) &&
223	<	mu->BestTrk()->NHits() > 10 &&
224	<	RChi2 < 10.0 &&
225	<	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
220	>	idpass = mu->BestTrk() != 0 &&
221	>	mu->NTrkLayersHit() > 5 &&
222	>	mu->IsPFMuon() == kTRUE &&
223	>	mu->BestTrk()->NPixelHits() > 0 &&
224	>	RChi2 < 10.0;
225	>	break;
226	>	// 2012 WW analysis for 42x (there is no PFMuon link)
227	>	case kWWMuIdV1:
228	>	idpass = mu->BestTrk() != 0 &&
229	>	mu->NTrkLayersHit() > 5 &&
230	>	mu->BestTrk()->NPixelHits() > 0 &&
231	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
232	>	mu->TrkKink() < 20.0;
233		break;
234	<	case kWWMuId:
234	>	// 2010 WW analysis
235	>	case kWWMuIdV2:
236		idpass = mu->BestTrk() != 0 &&
237		mu->BestTrk()->NHits() > 10 &&
182	–	RChi2 < 10.0 &&
183	–	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
238		mu->BestTrk()->NPixelHits() > 0 &&
185	–	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight) &&
239		mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
240		break;
241	+	// 2011 WW analysis
242	+	case kWWMuIdV3:
243	+	idpass = mu->BestTrk() != 0 &&
244	+	mu->BestTrk()->NHits() > 10 &&
245	+	mu->BestTrk()->NPixelHits() > 0 &&
246	+	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
247	+	mu->TrkKink() < 20.0;
248	+	break;
249	+	// 2012 WW analysis
250	+	case kWWMuIdV4:
251	+	idpass = mu->BestTrk() != 0 &&
252	+	mu->NTrkLayersHit() > 5 &&
253	+	mu->IsPFMuon() == kTRUE &&
254	+	mu->BestTrk()->NPixelHits() > 0 &&
255	+	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
256	+	mu->TrkKink() < 20.0;
257	+	break;
258	+	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	+	mu->TrkKink() < 20.0
267	+	);
268	+	idpass =  passDenominatorM2;
269	+	//only evaluate MVA if muon passes M2 denominator to save time
270	+	if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
271	+	}
272	+	break;
273		case kNoId:
274		idpass = kTRUE;
275		break;
#	Line 195 | Line 280 | void MuonIDMod::Process()
280		if (!idpass)
281		continue;
282
283	+	Double_t Rho = 0.;
284	+	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	+	case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
301	+	Rho = rho->RhoKt6PFJets();
302	+	break;
303	+	default:
304	+	Rho = rho->Rho();
305	+	}
306	+
307	+	if ((TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho()))) Rho = 0.;
308	+	}
309	+
310		Bool_t isocut = kFALSE;
311		switch (fMuIsoType) {
312		case kTrackCalo:
#	Line 208 | Line 320 | void MuonIDMod::Process()
320		break;
321		case kTrackCaloSliding:
322		{
323	<	//Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
324	<	//if(beta == 0) beta = 1.0;
325	<	const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
326	<	Double_t totalIso =  mu->IsoR03SumPt() + TMath::Max(mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3, 0.0);
327	<	if (totalIso < (mu->Pt()*fCombIsolationCut) )
328	<	isocut = kTRUE;
329	<
330	<	if     (fReverseIsoCut == kTRUE &&
219	<	isocut == kFALSE && totalIso < 10)
220	<	isocut = kTRUE;
221	<	else if(fReverseIsoCut == kTRUE)
222	<	isocut = kFALSE;
323	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
324	>	// trick to change the signal region cut
325	>	double theIsoCut = fCombIsolationCut;
326	>	if(theIsoCut < 0.20){
327	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
328	>	else                 theIsoCut = 0.10;
329	>	}
330	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
331		}
332		break;
333	<	case kTrackCaloSlidingNoCorrection:
333	>	case kTrackCaloSlidingNoCorrection:
334		{
335		Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
336		1.0 * mu->IsoR03EmEt()  +
337		1.0 * mu->IsoR03HadEt();
338	<	if (totalIso < (mu->Pt()*fCombIsolationCut) )
339	<	isocut = kTRUE;
338	>	// trick to change the signal region cut
339	>	double theIsoCut = fCombIsolationCut;
340	>	if(theIsoCut < 0.20){
341	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
342	>	else                 theIsoCut = 0.10;
343	>	}
344	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
345		}
346		break;
347	<	case kPFIso:
348	<	{
349	<	Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
350	<	if(beta == 0) beta = 1.0;
351	<	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 0, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
352	<	if (totalIso < (mu->Pt()*fCombIsolationCut) )
347	>	case kCombinedRelativeConeAreaCorrected:
348	>	{
349	>	//const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0); // Fabian: made Rho customable
350	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
351	>	double theIsoCut = fCombRelativeIsolationCut;
352	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
353	>	}
354	>	break;
355	>	case kCombinedRelativeEffectiveAreaCorrected:
356	>	{
357	>	Double_t tmpRho = Rho;   // Fabian: made the Rho type customable.
358	>	//if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))
359	>	//tmpRho = fPileupEnergyDensity->At(0)->Rho();
360	>
361	>	isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
362	>	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
363	>	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
364	>	) < (mu->Pt()* 0.40);
365	>	}
366	>	break;
367	>	case kPFIso:
368	>	{
369	>	Double_t pfIsoCutValue = 9999;
370	>	if(fPFIsolationCut > 0){
371	>	pfIsoCutValue = fPFIsolationCut;
372	>	} else {
373	>	if (mu->AbsEta() < 1.479) {
374	>	if (mu->Pt() > 20) {
375	>	pfIsoCutValue = 0.13;
376	>	} else {
377	>	pfIsoCutValue = 0.06;
378	>	}
379	>	} else {
380	>	if (mu->Pt() > 20) {
381	>	pfIsoCutValue = 0.09;
382	>	} else {
383	>	pfIsoCutValue = 0.05;
384	>	}
385	>	}
386	>	}
387	>	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
388	>	if (totalIso < (mu->Pt()*pfIsoCutValue) )
389	>	isocut = kTRUE;
390	>	}
391	>	break;
392	>	case kPFRadialIso:
393	>	{
394	>	Double_t pfIsoCutValue = 9999;
395	>	if(fPFIsolationCut > 0){
396	>	pfIsoCutValue = fPFIsolationCut;
397	>	} else {
398	>	if (mu->Pt() > 20) {
399	>	pfIsoCutValue = 0.10;
400	>	} else {
401	>	pfIsoCutValue = 0.05;
402	>	}
403	>	}
404	>	Double_t totalIso =  IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
405	>	if (totalIso < (mu->Pt()*pfIsoCutValue) )
406		isocut = kTRUE;
407		}
408	<	break;
409	<	case kPFIsoNoL:
408	>	break;
409	>	case kPFIsoEffectiveAreaCorrected:
410	>	{
411	>	Double_t pfIsoCutValue = 9999;
412	>	if(fPFIsolationCut > 0){
413	>	pfIsoCutValue = fPFIsolationCut;
414	>	} else {
415	>	pfIsoCutValue = fPFIsolationCut; //leave it like this for now
416	>	}
417	>	Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
418	>	- Rho * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
419	>	//- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());  // Fabian: made Rho-type customable
420	>	isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
421	>	break;
422	>	}
423	>
424	>
425	>	case kPFIsoNoL:
426		{
427		fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
428		fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
429
430	<	Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
431	<	if(beta == 0) beta = 1.0;
432	<	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 3, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
433	<	if (totalIso < (mu->Pt()*fCombIsolationCut) )
430	>	Double_t pfIsoCutValue = 9999;
431	>	if(fPFIsolationCut > 0){
432	>	pfIsoCutValue = fPFIsolationCut;
433	>	} else {
434	>	if (mu->AbsEta() < 1.479) {
435	>	if (mu->Pt() > 20) {
436	>	pfIsoCutValue = 0.13;
437	>	} else {
438	>	pfIsoCutValue = 0.06;
439	>	}
440	>	} else {
441	>	if (mu->Pt() > 20) {
442	>	pfIsoCutValue = 0.09;
443	>	} else {
444	>	pfIsoCutValue = 0.05;
445	>	}
446	>	}
447	>	}
448	>	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
449	>	if (totalIso < (mu->Pt()*pfIsoCutValue) )
450		isocut = kTRUE;
451		}
452		break;
453	+	case kMVAIso_BDTG_IDIso:
454	+	{
455	+
456	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
457	+	isocut = (totalIso < (mu->Pt()*0.4));
458	+
459	+	}
460	+	break;
461	+	case kIsoRingsV0_BDTG_Iso:
462	+	{
463	+
464	+	isocut = PassMuonIsoRingsV0_BDTG_Iso(mu, fVertices->At(0), fPileupEnergyDensity);
465	+
466	+	}
467	+	break;
468	+	case kIsoDeltaR:
469	+	{
470	+
471	+	isocut = PassMuonIsoDeltaR(mu, fVertices->At(0), fPileupEnergyDensity);
472	+
473	+	}
474	+	break;
475		case kNoIso:
476		isocut = kTRUE;
477		break;
#	Line 263 | Line 483 | void MuonIDMod::Process()
483		if (isocut == kFALSE)
484		continue;
485
486	+	// apply d0 cut
487		if (fApplyD0Cut) {
488	<	Bool_t passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut);
488	>	Bool_t passD0cut = kTRUE;
489	>	if(fD0Cut < 0.05) { // trick to change the signal region cut
490	>	if      (mu->Pt() >  20.0) fD0Cut = 0.02;
491	>	else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
492	>	}
493	>	if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
494	>	else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
495		if (!passD0cut)
496		continue;
497		}
498
499	+	// apply dz cut
500	+	if (fApplyDZCut) {
501	+	Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
502	+	if (!passDZcut)
503	+	continue;
504	+	}
505	+
506		// add good muon
507		CleanMuons->Add(mu);
508		}
#	Line 290 | Line 524 | void MuonIDMod::SlaveBegin()
524		if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
525		ReqEventObject(fMuonBranchName, fMuons, kTRUE);
526		}
527	+	ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
528		ReqEventObject(fTrackName, fTracks, kTRUE);
529		ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
530	<	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0) {
530	>	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
531	>	|| fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
532	>	|| fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
533	>	|| fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
534	>	|| fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
535	>	|| fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0
536	>	|| fMuonIsoType.CompareTo("IsoDeltaR") == 0
537	>	) {
538		ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
539		}
540
541	+
542		if (fMuonIDType.CompareTo("WMuId") == 0)
543		fMuIDType = kWMuId;
544		else if (fMuonIDType.CompareTo("ZMuId") == 0)
#	Line 304 | Line 547 | void MuonIDMod::SlaveBegin()
547		fMuIDType = kTight;
548		else if (fMuonIDType.CompareTo("Loose") == 0)
549		fMuIDType = kLoose;
550	<	else if (fMuonIDType.CompareTo("WWMuId") == 0)
551	<	fMuIDType = kWWMuId;
550	>	else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
551	>	fMuIDType = kWWMuIdV1;
552	>	else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
553	>	fMuIDType = kWWMuIdV2;
554	>	else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
555	>	fMuIDType = kWWMuIdV3;
556	>	else if (fMuonIDType.CompareTo("WWMuIdV4") == 0)
557	>	fMuIDType = kWWMuIdV4;
558		else if (fMuonIDType.CompareTo("NoId") == 0)
559		fMuIDType = kNoId;
560		else if (fMuonIDType.CompareTo("Custom") == 0) {
561		fMuIDType = kCustomId;
562		SendError(kWarning, "SlaveBegin",
563		"Custom muon identification is not yet implemented.");
564	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
565	+	fMuIDType = kMVAID_BDTG_IDIso;
566		} else {
567		SendError(kAbortAnalysis, "SlaveBegin",
568		"The specified muon identification %s is not defined.",
#	Line 327 | Line 578 | void MuonIDMod::SlaveBegin()
578		fMuIsoType = kTrackCaloSliding;
579		else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
580		fMuIsoType = kTrackCaloSlidingNoCorrection;
581	+	else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
582	+	fMuIsoType = kCombinedRelativeConeAreaCorrected;
583	+	else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
584	+	fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
585		else if (fMuonIsoType.CompareTo("PFIso") == 0)
586		fMuIsoType = kPFIso;
587	+	else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
588	+	fMuIsoType = kPFRadialIso;
589	+	else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
590	+	fMuIsoType = kPFIsoEffectiveAreaCorrected;
591		else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
592		fMuIsoType = kPFIsoNoL;
593		else if (fMuonIsoType.CompareTo("NoIso") == 0)
#	Line 337 | Line 596 | void MuonIDMod::SlaveBegin()
596		fMuIsoType = kCustomIso;
597		SendError(kWarning, "SlaveBegin",
598		"Custom muon isolation is not yet implemented.");
599	+	} else if (fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0) {
600	+	fMuIsoType = kMVAIso_BDTG_IDIso;
601	+	} else if (fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0) {
602	+	fMuIsoType = kIsoRingsV0_BDTG_Iso;
603	+	} else if (fMuonIsoType.CompareTo("IsoDeltaR") == 0) {
604	+	fMuIsoType = kIsoDeltaR;
605		} else {
606		SendError(kAbortAnalysis, "SlaveBegin",
607		"The specified muon isolation %s is not defined.",
#	Line 348 | Line 613 | void MuonIDMod::SlaveBegin()
613		fMuClassType = kAll;
614		else if (fMuonClassType.CompareTo("Global") == 0)
615		fMuClassType = kGlobal;
616	+	else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
617	+	fMuClassType = kGlobalTracker;
618		else if (fMuonClassType.CompareTo("Standalone") == 0)
619		fMuClassType = kSta;
620		else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
#	Line 356 | Line 623 | void MuonIDMod::SlaveBegin()
623		fMuClassType = kCaloMuon;
624		else if (fMuonClassType.CompareTo("TrackerBased") == 0)
625		fMuClassType = kTrackerBased;
626	+	else if (fMuonClassType.CompareTo("GlobalOnly") == 0)
627	+	fMuClassType = kGlobalOnly;
628		else {
629		SendError(kAbortAnalysis, "SlaveBegin",
630		"The specified muon class %s is not defined.",
631		fMuonClassType.Data());
632		return;
633		}
634	+
635	+
636	+	//If we use MVA ID, need to load MVA weights
637	+	if     (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
638	+	fMuonTools = new MuonTools();
639	+	fMuonIDMVA = new MuonIDMVA();
640	+	fMuonIDMVA->Initialize("BDTG method",
641	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin0_IDIsoCombined_BDTG.weights.xml"))),
642	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin0_IDIsoCombined_BDTG.weights.xml"))),
643	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin1_IDIsoCombined_BDTG.weights.xml"))),
644	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin1_IDIsoCombined_BDTG.weights.xml"))),
645	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin2_IDIsoCombined_BDTG.weights.xml"))),
646	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin2_IDIsoCombined_BDTG.weights.xml"))),
647	+	MuonIDMVA::kIDIsoCombinedDetIso,
648	+	fTheRhoType);
649	+	}
650	+	else if(fMuIsoType == kIsoRingsV0_BDTG_Iso) {
651	+	std::vector<std::string> muonidiso_weightfiles;
652	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml"))));
653	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml"))));
654	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml"))));
655	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml"))));
656	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml"))));
657	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml"))));
658	+	fMuonTools = new MuonTools();
659	+	fMuonIDMVA = new MuonIDMVA();
660	+	fMuonIDMVA->Initialize("MuonIso_BDTG_IsoRings",
661	+	MuonIDMVA::kIsoRingsV0,
662	+	kTRUE,
663	+	muonidiso_weightfiles,
664	+	fTheRhoType);
665	+	}
666	+	else if(fMuIsoType == kIsoDeltaR) {
667	+	std::vector<std::string> muonidiso_weightfiles;
668	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LB_BDT.weights.xml"))));
669	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LE_BDT.weights.xml"))));
670	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HB_BDT.weights.xml"))));
671	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HE_BDT.weights.xml"))));
672	+	fMuonTools = new MuonTools();
673	+	fMuonIDMVA = new MuonIDMVA();
674	+	fMuonIDMVA->Initialize("muonHZZ2012IsoDRMVA",
675	+	MuonIDMVA::kIsoDeltaR,
676	+	kTRUE,
677	+	muonidiso_weightfiles,
678	+	fTheRhoType);
679	+	}
680	+
681	+	}
682	+
683	+
684	+	//--------------------------------------------------------------------------------------------------
685	+	Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
686	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
687	+	{
688	+
689	+	const Track *muTrk=0;
690	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
691	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
692	+
693	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
694	+
695	+	Int_t subdet = 0;
696	+	if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
697	+	else subdet = 1;
698	+	Int_t ptBin = 0;
699	+	if (muTrk->Pt() > 14.5) ptBin = 1;
700	+	if (muTrk->Pt() > 20.0) ptBin = 2;
701	+
702	+	Int_t MVABin = -1;
703	+	if      (subdet == 0 && ptBin == 0) MVABin = 0;
704	+	else if (subdet == 1 && ptBin == 0) MVABin = 1;
705	+	else if (subdet == 0 && ptBin == 1) MVABin = 2;
706	+	else if (subdet == 1 && ptBin == 1) MVABin = 3;
707	+	else if (subdet == 0 && ptBin == 2) MVABin = 4;
708	+	else if (subdet == 1 && ptBin == 2) MVABin = 5;
709	+
710	+	Double_t MVACut = -999;
711	+	if      (MVABin == 0) MVACut = -0.5618;
712	+	else if (MVABin == 1) MVACut = -0.3002;
713	+	else if (MVABin == 2) MVACut = -0.4642;
714	+	else if (MVABin == 3) MVACut = -0.2478;
715	+	else if (MVABin == 4) MVACut =  0.1706;
716	+	else if (MVABin == 5) MVACut =  0.8146;
717	+
718	+	if (MVAValue > MVACut) return kTRUE;
719	+	return kFALSE;
720	+	}
721	+
722	+	//--------------------------------------------------------------------------------------------------
723	+	Bool_t MuonIDMod::PassMuonIsoRingsV0_BDTG_Iso(const Muon *mu, const Vertex *vertex,
724	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
725	+	{
726	+
727	+	Bool_t isDebug = kFALSE;
728	+	const Track *muTrk=0;
729	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
730	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
731	+
732	+	ElectronOArr *tempElectrons = new  ElectronOArr;
733	+	MuonOArr     *tempMuons     = new  MuonOArr;
734	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,
735	+	PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,isDebug);
736	+	delete tempElectrons;
737	+	delete tempMuons;
738	+
739	+	Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
740	+
741	+	Double_t MVACut = -1.0;
742	+	Double_t eta = mu->AbsEta();
743	+	if     (mu->Pt() <  20 && eta <  1.479) MVACut = 0.86;
744	+	else if(mu->Pt() <  20 && eta >= 1.479) MVACut = 0.82;
745	+	else if(mu->Pt() >= 20 && eta <  1.479) MVACut = 0.82;
746	+	else if(mu->Pt() >= 20 && eta >= 1.479) MVACut = 0.86;
747	+
748	+	if(fPFIsolationCut > -1.0) MVACut = fPFIsolationCut;
749	+
750	+	if(isDebug == kTRUE){
751	+	printf("PassMuonIsoRingsV0_BDTG_IsoDebug: %d, pt, eta = %f, %f, rho = %f(%f) : RingsMVA = %f, bin: %d\n",
752	+	GetEventHeader()->EvtNum(),mu->Pt(), mu->Eta(),
753	+	fPileupEnergyDensity->At(0)->Rho(),fPileupEnergyDensity->At(0)->RhoKt6PFJets(),MVAValue,MVABin);
754	+	}
755	+
756	+	if (MVAValue > MVACut) return kTRUE;
757	+	return kFALSE;
758	+	}
759	+
760	+	//--------------------------------------------------------------------------------------------------
761	+	Bool_t MuonIDMod::PassMuonIsoDeltaR(const Muon *mu, const Vertex *vertex,
762	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
763	+	{
764	+
765	+	const Track *muTrk=0;
766	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
767	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
768	+
769	+	ElectronOArr *tempElectrons = new  ElectronOArr;
770	+	MuonOArr     *tempMuons     = new  MuonOArr;
771	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFNoPileUpCands,
772	+	PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,kFALSE);
773	+	delete tempElectrons;
774	+	delete tempMuons;
775	+
776	+	Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
777	+
778	+	Double_t MVACut = -999;
779	+	if      (MVABin == 0) MVACut =  0.000;
780	+	else if (MVABin == 1) MVACut =  0.000;
781	+	else if (MVABin == 2) MVACut =  0.000;
782	+	else if (MVABin == 3) MVACut =  0.000;
783	+
784	+	if (MVAValue > MVACut) return kTRUE;
785	+	return kFALSE;
786	+	}
787	+
788	+	//--------------------------------------------------------------------------------------------------
789	+	void MuonIDMod::Terminate()
790	+	{
791	+	// Run finishing code on the computer (slave) that did the analysis
792	+	delete fMuonIDMVA;
793	+
794	+	delete fMuonTools;
795		}

Diff Legend

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