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.62 by sixie, Wed Jan 4 16:31:00 2012 UTC

#	Line 19 | Line 19 | ClassImp(mithep::MuonIDMod)
19		fNonIsolatedMuonsName("random"),
20		fNonIsolatedElectronsName("random"),
21		fVertexName(ModNames::gkGoodVertexesName),
22	+	fBeamSpotName(Names::gkBeamSpotBrn),
23		fTrackName(Names::gkTrackBrn),
24		fPFCandidatesName(Names::gkPFCandidatesBrn),
25	<	fMuonIDType("WWMuId"),
26	<	fMuonIsoType("TrackCaloSliding"),
25	>	fMuonIDType("WWMuIdV3"),
26	>	fMuonIsoType("PFIso"),
27		fMuonClassType("Global"),
28		fTrackIsolationCut(3.0),
29		fCaloIsolationCut(3.0),
30		fCombIsolationCut(0.15),
31	+	fCombRelativeIsolationCut(0.15),
32	+	fPFIsolationCut(-1.0),
33		fMuonPtMin(10),
34		fApplyD0Cut(kTRUE),
35	+	fApplyDZCut(kTRUE),
36		fD0Cut(0.020),
37	+	fDZCut(0.10),
38	+	fWhichVertex(-1),
39		fEtaCut(2.4),
34	–	fReverseIsoCut(kFALSE),
35	–	fReverseD0Cut(kFALSE),
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)
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 59 | Line 73 | void MuonIDMod::Process()
73		else {
74		fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
75		}
76	+	LoadEventObject(fBeamSpotName, fBeamSpot);
77		LoadEventObject(fTrackName, fTracks);
78		LoadEventObject(fPFCandidatesName, fPFCandidates);
79	<	if(fMuIsoType == kTrackCaloSliding) {
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
#	Line 94 | Line 114 | void MuonIDMod::Process()
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->HasStandaloneTrk();
133		if (pass) {
#	Line 176 | Line 210 | void MuonIDMod::Process()
210		RChi2 < 10.0 &&
211		mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
212		break;
213	<	case kWWMuId:
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 &&
185	–	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight) &&
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	+	( IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
245	+	- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta())
246	+	) < (mu->Pt()* 0.40)
247	+	);
248	+
249	+	idpass = passDenominatorM2;
250	+	//only evaluate MVA if muon passes M2 denominator to save time
251	+	if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
252	+	}
253	+	break;
254		case kNoId:
255		idpass = kTRUE;
256		break;
#	Line 208 | Line 274 | void MuonIDMod::Process()
274		break;
275		case kTrackCaloSliding:
276		{
211	–	//Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
212	–	//if(beta == 0) beta = 1.0;
277		const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
278	<	Double_t totalIso =  mu->IsoR03SumPt() + TMath::Max(mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3, 0.0);
279	<	if (totalIso < (mu->Pt()*fCombIsolationCut) )
280	<	isocut = kTRUE;
281	<
282	<	if     (fReverseIsoCut == kTRUE &&
283	<	isocut == kFALSE && totalIso < 10)
284	<	isocut = kTRUE;
285	<	else if(fReverseIsoCut == kTRUE)
222	<	isocut = kFALSE;
278	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
279	>	// trick to change the signal region cut
280	>	double theIsoCut = fCombIsolationCut;
281	>	if(theIsoCut < 0.20){
282	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
283	>	else                 theIsoCut = 0.10;
284	>	}
285	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
286		}
287		break;
288		case kTrackCaloSlidingNoCorrection:
#	Line 227 | Line 290 | void MuonIDMod::Process()
290		Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
291		1.0 * mu->IsoR03EmEt()  +
292		1.0 * mu->IsoR03HadEt();
293	<	if (totalIso < (mu->Pt()*fCombIsolationCut) )
294	<	isocut = kTRUE;
293	>	// trick to change the signal region cut
294	>	double theIsoCut = fCombIsolationCut;
295	>	if(theIsoCut < 0.20){
296	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
297	>	else                 theIsoCut = 0.10;
298	>	}
299	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
300		}
301		break;
302	+	case kCombinedRelativeConeAreaCorrected:
303	+	{
304	+	const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
305	+	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
306	+	double theIsoCut = fCombRelativeIsolationCut;
307	+	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
308	+	}
309	+	break;
310		case kPFIso:
311		{
312	<	Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
313	<	if(beta == 0) beta = 1.0;
314	<	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 0, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
315	<	if (totalIso < (mu->Pt()*fCombIsolationCut) )
312	>	Double_t pfIsoCutValue = 9999;
313	>	if(fPFIsolationCut > 0){
314	>	pfIsoCutValue = fPFIsolationCut;
315	>	} else {
316	>	if (mu->AbsEta() < 1.479) {
317	>	if (mu->Pt() > 20) {
318	>	pfIsoCutValue = 0.13;
319	>	} else {
320	>	pfIsoCutValue = 0.06;
321	>	}
322	>	} else {
323	>	if (mu->Pt() > 20) {
324	>	pfIsoCutValue = 0.09;
325	>	} else {
326	>	pfIsoCutValue = 0.05;
327	>	}
328	>	}
329	>	}
330	>	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
331	>	if (totalIso < (mu->Pt()*pfIsoCutValue) )
332		isocut = kTRUE;
333		}
334		break;
335	+	case kPFIsoEffectiveAreaCorrected:
336	+	{
337	+	Double_t pfIsoCutValue = 9999;
338	+	if(fPFIsolationCut > 0){
339	+	pfIsoCutValue = fPFIsolationCut;
340	+	} else {
341	+	pfIsoCutValue = fPFIsolationCut; //leave it like this for now
342	+	}
343	+	Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
344	+	- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
345	+	isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
346	+	break;
347	+	}
348		case kPFIsoNoL:
349		{
350		fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
351		fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
352
353	<	Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
354	<	if(beta == 0) beta = 1.0;
355	<	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 3, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
356	<	if (totalIso < (mu->Pt()*fCombIsolationCut) )
353	>	Double_t pfIsoCutValue = 9999;
354	>	if(fPFIsolationCut > 0){
355	>	pfIsoCutValue = fPFIsolationCut;
356	>	} else {
357	>	if (mu->AbsEta() < 1.479) {
358	>	if (mu->Pt() > 20) {
359	>	pfIsoCutValue = 0.13;
360	>	} else {
361	>	pfIsoCutValue = 0.06;
362	>	}
363	>	} else {
364	>	if (mu->Pt() > 20) {
365	>	pfIsoCutValue = 0.09;
366	>	} else {
367	>	pfIsoCutValue = 0.05;
368	>	}
369	>	}
370	>	}
371	>	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
372	>	if (totalIso < (mu->Pt()*pfIsoCutValue) )
373		isocut = kTRUE;
374		}
375		break;
376	+	case kMVAIso_BDTG_IDIso:
377	+	isocut = ( IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
378	+	- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta())
379	+	) < (mu->Pt()* 0.40);
380	+	break;
381		case kNoIso:
382		isocut = kTRUE;
383		break;
#	Line 263 | Line 389 | void MuonIDMod::Process()
389		if (isocut == kFALSE)
390		continue;
391
392	+	// apply d0 cut
393		if (fApplyD0Cut) {
394	<	Bool_t passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut);
394	>	Bool_t passD0cut = kTRUE;
395	>	if(fD0Cut < 0.05) { // trick to change the signal region cut
396	>	if      (mu->Pt() >  20.0) fD0Cut = 0.02;
397	>	else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
398	>	}
399	>	if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
400	>	else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
401		if (!passD0cut)
402		continue;
403		}
404
405	+	// apply dz cut
406	+	if (fApplyDZCut) {
407	+	Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
408	+	if (!passDZcut)
409	+	continue;
410	+	}
411	+
412		// add good muon
413		CleanMuons->Add(mu);
414		}
#	Line 290 | Line 430 | void MuonIDMod::SlaveBegin()
430		if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
431		ReqEventObject(fMuonBranchName, fMuons, kTRUE);
432		}
433	+	ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
434		ReqEventObject(fTrackName, fTracks, kTRUE);
435		ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
436	<	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0) {
436	>	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
437	>	|| fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
438	>	|| fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
439	>	|| fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
440	>	) {
441		ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
442		}
443
444	+
445		if (fMuonIDType.CompareTo("WMuId") == 0)
446		fMuIDType = kWMuId;
447		else if (fMuonIDType.CompareTo("ZMuId") == 0)
#	Line 304 | Line 450 | void MuonIDMod::SlaveBegin()
450		fMuIDType = kTight;
451		else if (fMuonIDType.CompareTo("Loose") == 0)
452		fMuIDType = kLoose;
453	<	else if (fMuonIDType.CompareTo("WWMuId") == 0)
454	<	fMuIDType = kWWMuId;
453	>	else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
454	>	fMuIDType = kWWMuIdV1;
455	>	else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
456	>	fMuIDType = kWWMuIdV2;
457	>	else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
458	>	fMuIDType = kWWMuIdV3;
459		else if (fMuonIDType.CompareTo("NoId") == 0)
460		fMuIDType = kNoId;
461		else if (fMuonIDType.CompareTo("Custom") == 0) {
462		fMuIDType = kCustomId;
463		SendError(kWarning, "SlaveBegin",
464		"Custom muon identification is not yet implemented.");
465	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
466	+	fMuIDType = kMVAID_BDTG_IDIso;
467		} else {
468		SendError(kAbortAnalysis, "SlaveBegin",
469		"The specified muon identification %s is not defined.",
#	Line 327 | Line 479 | void MuonIDMod::SlaveBegin()
479		fMuIsoType = kTrackCaloSliding;
480		else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
481		fMuIsoType = kTrackCaloSlidingNoCorrection;
482	+	else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
483	+	fMuIsoType = kCombinedRelativeConeAreaCorrected;
484		else if (fMuonIsoType.CompareTo("PFIso") == 0)
485		fMuIsoType = kPFIso;
486	+	else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
487	+	fMuIsoType = kPFIsoEffectiveAreaCorrected;
488		else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
489		fMuIsoType = kPFIsoNoL;
490		else if (fMuonIsoType.CompareTo("NoIso") == 0)
#	Line 337 | Line 493 | void MuonIDMod::SlaveBegin()
493		fMuIsoType = kCustomIso;
494		SendError(kWarning, "SlaveBegin",
495		"Custom muon isolation is not yet implemented.");
496	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
497	+	fMuIsoType = kMVAIso_BDTG_IDIso;
498		} else {
499		SendError(kAbortAnalysis, "SlaveBegin",
500		"The specified muon isolation %s is not defined.",
#	Line 348 | Line 506 | void MuonIDMod::SlaveBegin()
506		fMuClassType = kAll;
507		else if (fMuonClassType.CompareTo("Global") == 0)
508		fMuClassType = kGlobal;
509	+	else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
510	+	fMuClassType = kGlobalTracker;
511		else if (fMuonClassType.CompareTo("Standalone") == 0)
512		fMuClassType = kSta;
513		else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
#	Line 362 | Line 522 | void MuonIDMod::SlaveBegin()
522		fMuonClassType.Data());
523		return;
524		}
525	+
526	+
527	+	//If we use MVA ID, need to load MVA weights
528	+	if(fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
529	+	fMuonTools = new MuonTools();
530	+	fMuonIDMVA = new MuonIDMVA();
531	+	fMuonIDMVA->Initialize("BDTG method",
532	+	fMuonMVAWeights_Subdet0Pt10To14p5,
533	+	fMuonMVAWeights_Subdet1Pt10To14p5,
534	+	fMuonMVAWeights_Subdet0Pt14p5To20,
535	+	fMuonMVAWeights_Subdet1Pt14p5To20,
536	+	fMuonMVAWeights_Subdet0Pt20ToInf,
537	+	fMuonMVAWeights_Subdet1Pt20ToInf,
538	+	MuonIDMVA::kV8);
539	+	}
540	+
541	+	}
542	+
543	+
544	+	//--------------------------------------------------------------------------------------------------
545	+	Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
546	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
547	+	{
548	+
549	+	const Track *muTrk=0;
550	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
551	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
552	+
553	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
554	+
555	+	Int_t subdet = 0;
556	+	if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
557	+	else subdet = 1;
558	+	Int_t ptBin = 0;
559	+	if (muTrk->Pt() > 14.5) ptBin = 1;
560	+	if (muTrk->Pt() > 20.0) ptBin = 2;
561	+
562	+	Int_t MVABin = -1;
563	+	if (subdet == 0 && ptBin == 0) MVABin = 0;
564	+	if (subdet == 1 && ptBin == 0) MVABin = 1;
565	+	if (subdet == 0 && ptBin == 1) MVABin = 2;
566	+	if (subdet == 1 && ptBin == 1) MVABin = 3;
567	+	if (subdet == 0 && ptBin == 2) MVABin = 4;
568	+	if (subdet == 1 && ptBin == 2) MVABin = 5;
569	+
570	+	Double_t MVACut = -999;
571	+	if (MVABin == 0) MVACut = -0.377;
572	+	if (MVABin == 1) MVACut = -0.0902;
573	+	if (MVABin == 2) MVACut = -0.221;
574	+	if (MVABin == 3) MVACut = -0.0154;
575	+	if (MVABin == 4) MVACut = 0.459;
576	+	if (MVABin == 5) MVACut = 0.9158;
577	+
578	+	if (MVAValue > MVACut) return kTRUE;
579	+	return kFALSE;
580		}

Diff Legend

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