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.57 by ceballos, Tue Oct 11 17:00:24 2011 UTC vs.
Revision 1.70 by ceballos, Sat Apr 28 19:10:00 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"),
#	Line 45 | Line 46 | ClassImp(mithep::MuonIDMod)
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	>	fMuonMVAWeights_Subdet0Pt10To14p5(""),
58	>	fMuonMVAWeights_Subdet1Pt10To14p5(""),
59	>	fMuonMVAWeights_Subdet0Pt14p5To20(""),
60	>	fMuonMVAWeights_Subdet1Pt14p5To20(""),
61	>	fMuonMVAWeights_Subdet0Pt20ToInf(""),
62	>	fMuonMVAWeights_Subdet1Pt20ToInf("")
63		{
64		// Constructor.
65		}
#	Line 68 | Line 78 | void MuonIDMod::Process()
78		LoadEventObject(fBeamSpotName, fBeamSpot);
79		LoadEventObject(fTrackName, fTracks);
80		LoadEventObject(fPFCandidatesName, fPFCandidates);
81	<	if(fMuIsoType == kTrackCaloSliding || fMuIsoType == kCombinedRelativeConeAreaCorrected) {
81	>	if(fMuIsoType == kTrackCaloSliding ||
82	>	fMuIsoType == kCombinedRelativeConeAreaCorrected ||
83	>	fMuIsoType == kPFIsoEffectiveAreaCorrected ||
84	>	fMuIsoType == kMVAIso_BDTG_IDIso
85	>	) {
86		LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
87		}
88	+	if(fMuIsoType == kPFRadialIso){
89	+	// Name is hardcoded, can be changed if someone feels to do it
90	+	fPFNoPileUpCands = GetObjThisEvt<PFCandidateCol>("PFNoPileUp");
91	+	}
92	+
93		MuonOArr *CleanMuons = new MuonOArr;
94		CleanMuons->SetName(fCleanMuonsName);
95
96		fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
97
98	<	for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
98	>	for (UInt_t i=0; i<fMuons->GetEntries() && fVertices->GetEntries() > 0 ; ++i) {
99		const Muon *mu = fMuons->At(i);
100
101		Bool_t pass = kFALSE;
#	Line 157 | Line 176 | void MuonIDMod::Process()
176		if (eta >= fEtaCut)
177		continue;
178
179	+
180	+	//***********************************************************************************************
181	+	//Debug Info For Lepton MVA
182	+	//***********************************************************************************************
183	+	if( fPrintMVADebugInfo &&
184	+	(fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso)
185	+	) {
186	+	cout << "Event: " << GetEventHeader()->RunNum() << " " << GetEventHeader()->LumiSec() << " "
187	+	<< GetEventHeader()->EvtNum() << " : Rho = " << fPileupEnergyDensity->At(0)->Rho()
188	+	<< " : Muon " << i << " "
189	+	<< endl;
190	+	fMuonIDMVA->MVAValue(mu,fVertices->At(0),fMuonTools,fPFCandidates,fPileupEnergyDensity,kTRUE);
191	+	}
192	+	//***********************************************************************************************
193	+
194	+
195		Double_t RChi2 = 0.0;
196		if     (mu->HasGlobalTrk()) {
197		RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
#	Line 219 | Line 254 | void MuonIDMod::Process()
254		mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
255		mu->TrkKink() < 20.0;
256		break;
257	+	case kMVAID_BDTG_IDIso:
258	+	{
259	+	Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
260	+	mu->BestTrk()->NHits() > 10 &&
261	+	mu->BestTrk()->NPixelHits() > 0 &&
262	+	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
263	+	MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
264	+	MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
265	+	mu->TrkKink() < 20.0
266	+	);
267	+	idpass =  passDenominatorM2;
268	+	//only evaluate MVA if muon passes M2 denominator to save time
269	+	if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
270	+	}
271	+	break;
272		case kNoId:
273		idpass = kTRUE;
274		break;
#	Line 275 | Line 325 | void MuonIDMod::Process()
325		if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
326		}
327		break;
328	+	case kCombinedRelativeEffectiveAreaCorrected:
329	+	{
330	+	Double_t tmpRho = 0;
331	+	if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))
332	+	tmpRho = fPileupEnergyDensity->At(0)->Rho();
333	+
334	+	isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
335	+	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
336	+	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
337	+	) < (mu->Pt()* 0.40);
338	+	}
339	+	break;
340		case kPFIso:
341	<	{
341	>	{
342		Double_t pfIsoCutValue = 9999;
343		if(fPFIsolationCut > 0){
344		pfIsoCutValue = fPFIsolationCut;
#	Line 300 | Line 362 | void MuonIDMod::Process()
362		isocut = kTRUE;
363		}
364		break;
365	+	case kPFRadialIso:
366	+	{
367	+	Double_t pfIsoCutValue = 9999;
368	+	if(fPFIsolationCut > 0){
369	+	pfIsoCutValue = fPFIsolationCut;
370	+	} else {
371	+	if (mu->Pt() > 20) {
372	+	pfIsoCutValue = 0.10;
373	+	} else {
374	+	pfIsoCutValue = 0.05;
375	+	}
376	+	}
377	+	Double_t totalIso =  IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
378	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
379	+	isocut = kTRUE;
380	+	}
381	+	break;
382	+	case kPFIsoEffectiveAreaCorrected:
383	+	{
384	+	Double_t pfIsoCutValue = 9999;
385	+	if(fPFIsolationCut > 0){
386	+	pfIsoCutValue = fPFIsolationCut;
387	+	} else {
388	+	pfIsoCutValue = fPFIsolationCut; //leave it like this for now
389	+	}
390	+	Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
391	+	- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
392	+	isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
393	+	break;
394	+	}
395		case kPFIsoNoL:
396		{
397		fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
398		fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
399
400	<	Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
401	<	if(beta == 0) beta = 1.0;
402	<	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), fNonIsolatedMuons, fNonIsolatedElectrons, 0.2, 1.0, 0.4, 0.0, 3);
403	<	if (totalIso < (mu->Pt()*fPFIsolationCut) )
400	>	Double_t pfIsoCutValue = 9999;
401	>	if(fPFIsolationCut > 0){
402	>	pfIsoCutValue = fPFIsolationCut;
403	>	} else {
404	>	if (mu->AbsEta() < 1.479) {
405	>	if (mu->Pt() > 20) {
406	>	pfIsoCutValue = 0.13;
407	>	} else {
408	>	pfIsoCutValue = 0.06;
409	>	}
410	>	} else {
411	>	if (mu->Pt() > 20) {
412	>	pfIsoCutValue = 0.09;
413	>	} else {
414	>	pfIsoCutValue = 0.05;
415	>	}
416	>	}
417	>	}
418	>	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
419	>	if (totalIso < (mu->Pt()*pfIsoCutValue) )
420		isocut = kTRUE;
421		}
422		break;
423	+	case kMVAIso_BDTG_IDIso:
424	+	{
425	+
426	+	// **************************************************************************
427	+	// Don't use effective area correction denominator. Instead use the old one.
428	+	// **************************************************************************
429	+
430	+	//         Double_t tmpRho = 0;
431	+	//         if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || isinf(fPileupEnergyDensity->At(0)->Rho())))
432	+	//           tmpRho = fPileupEnergyDensity->At(0)->Rho();
433	+
434	+	//         isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
435	+	//                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
436	+	//                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
437	+	//           ) < (mu->Pt()* 0.40);
438	+
439	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
440	+	isocut = (totalIso < (mu->Pt()*0.4));
441	+
442	+	}
443	+	break;
444		case kNoIso:
445		isocut = kTRUE;
446		break;
#	Line 368 | Line 497 | void MuonIDMod::SlaveBegin()
497		ReqEventObject(fTrackName, fTracks, kTRUE);
498		ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
499		if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
500	<	|| fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0) {
500	>	|| fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
501	>	|| fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
502	>	|| fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
503	>	|| fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
504	>	) {
505		ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
506		}
507
508	+
509		if (fMuonIDType.CompareTo("WMuId") == 0)
510		fMuIDType = kWMuId;
511		else if (fMuonIDType.CompareTo("ZMuId") == 0)
#	Line 392 | Line 526 | void MuonIDMod::SlaveBegin()
526		fMuIDType = kCustomId;
527		SendError(kWarning, "SlaveBegin",
528		"Custom muon identification is not yet implemented.");
529	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
530	+	fMuIDType = kMVAID_BDTG_IDIso;
531		} else {
532		SendError(kAbortAnalysis, "SlaveBegin",
533		"The specified muon identification %s is not defined.",
#	Line 409 | Line 545 | void MuonIDMod::SlaveBegin()
545		fMuIsoType = kTrackCaloSlidingNoCorrection;
546		else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
547		fMuIsoType = kCombinedRelativeConeAreaCorrected;
548	+	else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
549	+	fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
550		else if (fMuonIsoType.CompareTo("PFIso") == 0)
551		fMuIsoType = kPFIso;
552	+	else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
553	+	fMuIsoType = kPFRadialIso;
554	+	else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
555	+	fMuIsoType = kPFIsoEffectiveAreaCorrected;
556		else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
557		fMuIsoType = kPFIsoNoL;
558		else if (fMuonIsoType.CompareTo("NoIso") == 0)
#	Line 419 | Line 561 | void MuonIDMod::SlaveBegin()
561		fMuIsoType = kCustomIso;
562		SendError(kWarning, "SlaveBegin",
563		"Custom muon isolation is not yet implemented.");
564	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
565	+	fMuIsoType = kMVAIso_BDTG_IDIso;
566		} else {
567		SendError(kAbortAnalysis, "SlaveBegin",
568		"The specified muon isolation %s is not defined.",
#	Line 446 | Line 590 | void MuonIDMod::SlaveBegin()
590		fMuonClassType.Data());
591		return;
592		}
593	+
594	+
595	+	//If we use MVA ID, need to load MVA weights
596	+	if(fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
597	+	fMuonTools = new MuonTools();
598	+	fMuonIDMVA = new MuonIDMVA();
599	+	fMuonIDMVA->Initialize("BDTG method",
600	+	fMuonMVAWeights_Subdet0Pt10To14p5,
601	+	fMuonMVAWeights_Subdet1Pt10To14p5,
602	+	fMuonMVAWeights_Subdet0Pt14p5To20,
603	+	fMuonMVAWeights_Subdet1Pt14p5To20,
604	+	fMuonMVAWeights_Subdet0Pt20ToInf,
605	+	fMuonMVAWeights_Subdet1Pt20ToInf,
606	+	MuonIDMVA::kIDIsoCombinedDetIso);
607	+	}
608	+
609	+	}
610	+
611	+
612	+	//--------------------------------------------------------------------------------------------------
613	+	Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
614	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
615	+	{
616	+
617	+	const Track *muTrk=0;
618	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
619	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
620	+
621	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
622	+
623	+	Int_t subdet = 0;
624	+	if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
625	+	else subdet = 1;
626	+	Int_t ptBin = 0;
627	+	if (muTrk->Pt() > 14.5) ptBin = 1;
628	+	if (muTrk->Pt() > 20.0) ptBin = 2;
629	+
630	+	Int_t MVABin = -1;
631	+	if      (subdet == 0 && ptBin == 0) MVABin = 0;
632	+	else if (subdet == 1 && ptBin == 0) MVABin = 1;
633	+	else if (subdet == 0 && ptBin == 1) MVABin = 2;
634	+	else if (subdet == 1 && ptBin == 1) MVABin = 3;
635	+	else if (subdet == 0 && ptBin == 2) MVABin = 4;
636	+	else if (subdet == 1 && ptBin == 2) MVABin = 5;
637	+
638	+	Double_t MVACut = -999;
639	+	if      (MVABin == 0) MVACut = -0.5618;
640	+	else if (MVABin == 1) MVACut = -0.3002;
641	+	else if (MVABin == 2) MVACut = -0.4642;
642	+	else if (MVABin == 3) MVACut = -0.2478;
643	+	else if (MVABin == 4) MVACut =  0.1706;
644	+	else if (MVABin == 5) MVACut =  0.8146;
645	+
646	+	if (MVAValue > MVACut) return kTRUE;
647	+	return kFALSE;
648		}

Diff Legend

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