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.72 by fabstoec, Thu May 3 08:45:28 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 11 | 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	<	fMuonIDType("Loose"),
21	<	fMuonIsoType("TrackCaloSliding"),
20	>	fNonIsolatedMuonsName("random"),
21	>	fNonIsolatedElectronsName("random"),
22	>	fVertexName(ModNames::gkGoodVertexesName),
23	>	fBeamSpotName(Names::gkBeamSpotBrn),
24	>	fTrackName(Names::gkTrackBrn),
25	>	fPFCandidatesName(Names::gkPFCandidatesBrn),
26	>	fMuonIDType("WWMuIdV3"),
27	>	fMuonIsoType("PFIso"),
28		fMuonClassType("Global"),
29		fTrackIsolationCut(3.0),
30		fCaloIsolationCut(3.0),
31	<	fCombIsolationCut(5.0),
31	>	fCombIsolationCut(0.15),
32	>	fCombRelativeIsolationCut(0.15),
33	>	fPFIsolationCut(-1.0),
34		fMuonPtMin(10),
35	<	fMuons(0)
35	>	fApplyD0Cut(kTRUE),
36	>	fApplyDZCut(kTRUE),
37	>	fD0Cut(0.020),
38	>	fDZCut(0.10),
39	>	fWhichVertex(-1),
40	>	fEtaCut(2.4),
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),
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	>	fTheRhoType(RhoUtilities::DEFAULT)
65		{
66		// Constructor.
67		}
#	Line 30 | Line 71 | void MuonIDMod::Process()
71		{
72		// Process entries of the tree.
73
74	<	LoadBranch(fMuonBranchName);
74	>	if(fMuIsoType != kPFIsoNoL) {
75	>	LoadEventObject(fMuonBranchName, fMuons);
76	>	}
77	>	else {
78	>	fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
79	>	}
80	>	LoadEventObject(fBeamSpotName, fBeamSpot);
81	>	LoadEventObject(fTrackName, fTracks);
82	>	LoadEventObject(fPFCandidatesName, fPFCandidates);
83	>	if(fMuIsoType == kTrackCaloSliding ||
84	>	fMuIsoType == kCombinedRelativeConeAreaCorrected ||
85	>	fMuIsoType == kPFIsoEffectiveAreaCorrected ||
86	>	fMuIsoType == kMVAIso_BDTG_IDIso
87	>	) {
88	>	LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
89	>	}
90	>	if(fMuIsoType == kPFRadialIso){
91	>	// Name is hardcoded, can be changed if someone feels to do it
92	>	fPFNoPileUpCands = GetObjThisEvt<PFCandidateCol>("PFNoPileUp");
93	>	}
94
95		MuonOArr *CleanMuons = new MuonOArr;
96		CleanMuons->SetName(fCleanMuonsName);
97
98	<	for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
98	>	fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
99	>
100	>	for (UInt_t i=0; i<fMuons->GetEntries() && fVertices->GetEntries() > 0 ; ++i) {
101		const Muon *mu = fMuons->At(i);
102
103		Bool_t pass = kFALSE;
104	<	Double_t pt = -1; // make sure pt is taken from the correct track!
104	>	Double_t pt = 0;  // make sure pt is taken from the correct track!
105	>	Double_t eta = 0; // make sure eta is taken from the correct track!
106		switch (fMuClassType) {
107		case kAll:
108		pass = kTRUE;
109	<	pt = mu->Pt();
109	>	if (mu->HasTrk()) {
110	>	pt  = mu->Pt();
111	>	eta = TMath::Abs(mu->Eta());
112	>	}
113		break;
114		case kGlobal:
115	<	pass = (mu->GlobalTrk() != 0);
116	<	if (pass)
117	<	pt = mu->GlobalTrk()->Pt();
118	<	break;
115	>	pass = mu->HasGlobalTrk() && mu->IsTrackerMuon();
116	>	if (pass && mu->TrackerTrk()) {
117	>	pt  = mu->TrackerTrk()->Pt();
118	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
119	>	}
120	>	else {
121	>	pt  = mu->Pt();
122	>	eta = TMath::Abs(mu->Eta());
123	>	}
124	>	break;
125	>	case kGlobalTracker:
126	>	pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
127	>	(mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
128	>	(mu->IsTrackerMuon() &&
129	>	mu->Quality().Quality(MuonQuality::TMLastStationTight));
130	>	if (pass) {
131	>	pt  = mu->TrackerTrk()->Pt();
132	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
133	>	}
134	>	else {
135	>	pt  = mu->Pt();
136	>	eta = TMath::Abs(mu->Eta());
137	>	}
138	>	break;
139		case kSta:
140	<	pass = (mu->StandaloneTrk() != 0);
141	<	if (pass)
142	<	pt = mu->StandaloneTrk()->Pt();
143	<	break;
144	<	case kTrackerOnly:
145	<	pass = (mu->TrackerTrk() != 0);
146	<	if (pass)
147	<	pt = mu->TrackerTrk()->Pt();
140	>	pass = mu->HasStandaloneTrk();
141	>	if (pass) {
142	>	pt  = mu->StandaloneTrk()->Pt();
143	>	eta = TMath::Abs(mu->StandaloneTrk()->Eta());
144	>	}
145	>	break;
146	>	case kTrackerMuon:
147	>	pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
148	>	mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
149	>	if (pass) {
150	>	pt  = mu->TrackerTrk()->Pt();
151	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
152	>	}
153	>	break;
154	>	case kCaloMuon:
155	>	pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
156	>	if (pass) {
157	>	pt  = mu->TrackerTrk()->Pt();
158	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
159	>	}
160	>	break;
161	>	case kTrackerBased:
162	>	pass = mu->HasTrackerTrk();
163	>	if (pass) {
164	>	pt  = mu->TrackerTrk()->Pt();
165	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
166	>	}
167		break;
168		default:
169		break;
#	Line 70 | Line 175 | void MuonIDMod::Process()
175		if (pt <= fMuonPtMin)
176		continue;
177
178	+	if (eta >= fEtaCut)
179	+	continue;
180	+
181	+
182	+	//***********************************************************************************************
183	+	//Debug Info For Lepton MVA
184	+	//***********************************************************************************************
185	+	if( fPrintMVADebugInfo &&
186	+	(fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso)
187	+	) {
188	+	cout << "Event: " << GetEventHeader()->RunNum() << " " << GetEventHeader()->LumiSec() << " "
189	+	<< GetEventHeader()->EvtNum() << " : Rho = " << fPileupEnergyDensity->At(0)->Rho()
190	+	<< " : Muon " << i << " "
191	+	<< endl;
192	+	fMuonIDMVA->MVAValue(mu,fVertices->At(0),fMuonTools,fPFCandidates,fPileupEnergyDensity,kTRUE);
193	+	}
194	+	//***********************************************************************************************
195	+
196	+
197	+	Double_t RChi2 = 0.0;
198	+	if     (mu->HasGlobalTrk()) {
199	+	RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
200	+	}
201	+	else if(mu->BestTrk() != 0){
202	+	RChi2 = mu->BestTrk()->Chi2()/mu->BestTrk()->Ndof();
203	+	}
204		Bool_t idpass = kFALSE;
205		switch (fMuIDType) {
206	+	case kWMuId:
207	+	idpass = mu->BestTrk() != 0 &&
208	+	mu->BestTrk()->NHits() > 10 &&
209	+	RChi2 < 10.0 &&
210	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
211	+	mu->BestTrk()->NPixelHits() > 0 &&
212	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
213	+	break;
214	+	case kZMuId:
215	+	idpass = mu->BestTrk() != 0 &&
216	+	mu->BestTrk()->NHits() > 10 &&
217	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
218	+	mu->BestTrk()->NPixelHits() > 0 &&
219	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
220	+	break;
221		case kLoose:
222	<	idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationLoose) &&
223	<	fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityLoose);
222	>	idpass = mu->BestTrk() != 0 &&
223	>	mu->Quality().Quality(MuonQuality::TMOneStationLoose) &&
224	>	mu->Quality().Quality(MuonQuality::TM2DCompatibilityLoose) &&
225	>	mu->BestTrk()->NHits() > 10 &&
226	>	RChi2 < 10.0 &&
227	>	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
228		break;
229		case kTight:
230	<	idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationTight) &&
231	<	fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityTight);
230	>	idpass = mu->BestTrk() !=  0 &&
231	>	mu->Quality().Quality(MuonQuality::TMOneStationTight) &&
232	>	mu->Quality().Quality(MuonQuality::TM2DCompatibilityTight) &&
233	>	mu->BestTrk()->NHits() > 10 &&
234	>	RChi2 < 10.0 &&
235	>	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
236	>	break;
237	>	case kWWMuIdV1:
238	>	idpass = mu->BestTrk() != 0 &&
239	>	mu->BestTrk()->NHits() > 10 &&
240	>	mu->BestTrk()->NPixelHits() > 0 &&
241	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
242	>	RChi2 < 10.0 &&
243	>	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
244	>	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
245	>	break;
246	>	case kWWMuIdV2:
247	>	idpass = mu->BestTrk() != 0 &&
248	>	mu->BestTrk()->NHits() > 10 &&
249	>	mu->BestTrk()->NPixelHits() > 0 &&
250	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
251	>	break;
252	>	case kWWMuIdV3:
253	>	idpass = mu->BestTrk() != 0 &&
254	>	mu->BestTrk()->NHits() > 10 &&
255	>	mu->BestTrk()->NPixelHits() > 0 &&
256	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
257	>	mu->TrkKink() < 20.0;
258	>	break;
259	>	case kMVAID_BDTG_IDIso:
260	>	{
261	>	Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
262	>	mu->BestTrk()->NHits() > 10 &&
263	>	mu->BestTrk()->NPixelHits() > 0 &&
264	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
265	>	MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
266	>	MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
267	>	mu->TrkKink() < 20.0
268	>	);
269	>	idpass =  passDenominatorM2;
270	>	//only evaluate MVA if muon passes M2 denominator to save time
271	>	if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
272	>	}
273	>	break;
274	>	case kNoId:
275	>	idpass = kTRUE;
276		break;
277		default:
278		break;
#	Line 87 | Line 281 | void MuonIDMod::Process()
281		if (!idpass)
282		continue;
283
284	<	Bool_t isopass = kFALSE;
284	>	Double_t Rho = 0.;
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	>	default:
301	>	Rho = rho->Rho();
302	>	}
303	>
304	>	if ((TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho()))) Rho = 0.;
305	>
306	>	Bool_t isocut = kFALSE;
307		switch (fMuIsoType) {
308		case kTrackCalo:
309	<	isopass = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
309	>	isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
310		(mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
311		break;
312		case kTrackCaloCombined:
313	<	isopass = (1.0 * mu->IsoR03SumPt() + 1.0 * mu->IsoR03EmEt() +
314	<	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
313	>	isocut = (1.0 * mu->IsoR03SumPt() +
314	>	1.0 * mu->IsoR03EmEt()  +
315	>	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
316		break;
317		case kTrackCaloSliding:
318		{
319	<	Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
320	<	1.0 * mu->IsoR03EmEt() +
321	<	1.0 * mu->IsoR03HadEt();
322	<	if ((totalIso < (mu->Pt()-10.0)*5.0/15.0 && mu->Pt() <= 25) ||
323	<	(totalIso < 5.0 && mu->Pt() > 25) ||
324	<	totalIso <= 0)
325	<	isopass = kTRUE;
319	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
320	>	// trick to change the signal region cut
321	>	double theIsoCut = fCombIsolationCut;
322	>	if(theIsoCut < 0.20){
323	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
324	>	else                 theIsoCut = 0.10;
325	>	}
326	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
327	>	}
328	>	break;
329	>	case kTrackCaloSlidingNoCorrection:
330	>	{
331	>	Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
332	>	1.0 * mu->IsoR03EmEt()  +
333	>	1.0 * mu->IsoR03HadEt();
334	>	// trick to change the signal region cut
335	>	double theIsoCut = fCombIsolationCut;
336	>	if(theIsoCut < 0.20){
337	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
338	>	else                 theIsoCut = 0.10;
339	>	}
340	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
341	>	}
342	>	break;
343	>	case kCombinedRelativeConeAreaCorrected:
344	>	{
345	>	//const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0); // Fabian: made Rho customable
346	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
347	>	double theIsoCut = fCombRelativeIsolationCut;
348	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
349		}
350	+	break;
351	+	case kCombinedRelativeEffectiveAreaCorrected:
352	+	{
353	+	Double_t tmpRho = Rho;   // Fabian: made the Rho type customable.
354	+	//if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))
355	+	//tmpRho = fPileupEnergyDensity->At(0)->Rho();
356	+
357	+	isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
358	+	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
359	+	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
360	+	) < (mu->Pt()* 0.40);
361	+	}
362	+	break;
363	+	case kPFIso:
364	+	{
365	+	Double_t pfIsoCutValue = 9999;
366	+	if(fPFIsolationCut > 0){
367	+	pfIsoCutValue = fPFIsolationCut;
368	+	} else {
369	+	if (mu->AbsEta() < 1.479) {
370	+	if (mu->Pt() > 20) {
371	+	pfIsoCutValue = 0.13;
372	+	} else {
373	+	pfIsoCutValue = 0.06;
374	+	}
375	+	} else {
376	+	if (mu->Pt() > 20) {
377	+	pfIsoCutValue = 0.09;
378	+	} else {
379	+	pfIsoCutValue = 0.05;
380	+	}
381	+	}
382	+	}
383	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
384	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
385	+	isocut = kTRUE;
386	+	}
387	+	break;
388	+	case kPFRadialIso:
389	+	{
390	+	Double_t pfIsoCutValue = 9999;
391	+	if(fPFIsolationCut > 0){
392	+	pfIsoCutValue = fPFIsolationCut;
393	+	} else {
394	+	if (mu->Pt() > 20) {
395	+	pfIsoCutValue = 0.10;
396	+	} else {
397	+	pfIsoCutValue = 0.05;
398	+	}
399	+	}
400	+	Double_t totalIso =  IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
401	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
402	+	isocut = kTRUE;
403	+	}
404	+	break;
405	+	case kPFIsoEffectiveAreaCorrected:
406	+	{
407	+	Double_t pfIsoCutValue = 9999;
408	+	if(fPFIsolationCut > 0){
409	+	pfIsoCutValue = fPFIsolationCut;
410	+	} else {
411	+	pfIsoCutValue = fPFIsolationCut; //leave it like this for now
412	+	}
413	+	Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
414	+	- Rho * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
415	+	//- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());  // Fabian: made Rho-type customable
416	+	isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
417	+	break;
418	+	}
419	+	case kPFIsoNoL:
420	+	{
421	+	fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
422	+	fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
423	+
424	+	Double_t pfIsoCutValue = 9999;
425	+	if(fPFIsolationCut > 0){
426	+	pfIsoCutValue = fPFIsolationCut;
427	+	} else {
428	+	if (mu->AbsEta() < 1.479) {
429	+	if (mu->Pt() > 20) {
430	+	pfIsoCutValue = 0.13;
431	+	} else {
432	+	pfIsoCutValue = 0.06;
433	+	}
434	+	} else {
435	+	if (mu->Pt() > 20) {
436	+	pfIsoCutValue = 0.09;
437	+	} else {
438	+	pfIsoCutValue = 0.05;
439	+	}
440	+	}
441	+	}
442	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
443	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
444	+	isocut = kTRUE;
445	+	}
446	+	break;
447	+	case kMVAIso_BDTG_IDIso:
448	+	{
449	+
450	+	// **************************************************************************
451	+	// Don't use effective area correction denominator. Instead use the old one.
452	+	// **************************************************************************
453	+
454	+	//         Double_t tmpRho = 0;
455	+	//         if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || isinf(fPileupEnergyDensity->At(0)->Rho())))
456	+	//           tmpRho = fPileupEnergyDensity->At(0)->Rho();
457	+
458	+	//         isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
459	+	//                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
460	+	//                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
461	+	//           ) < (mu->Pt()* 0.40);
462	+
463	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
464	+	isocut = (totalIso < (mu->Pt()*0.4));
465	+
466	+	}
467		break;
468		case kNoIso:
469	<	isopass = kTRUE;
469	>	isocut = kTRUE;
470		break;
471		case kCustomIso:
472		default:
473		break;
474		}
475
476	<	if (!isopass)
476	>	if (isocut == kFALSE)
477		continue;
478
479	+	// apply d0 cut
480	+	if (fApplyD0Cut) {
481	+	Bool_t passD0cut = kTRUE;
482	+	if(fD0Cut < 0.05) { // trick to change the signal region cut
483	+	if      (mu->Pt() >  20.0) fD0Cut = 0.02;
484	+	else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
485	+	}
486	+	if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
487	+	else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
488	+	if (!passD0cut)
489	+	continue;
490	+	}
491	+
492	+	// apply dz cut
493	+	if (fApplyDZCut) {
494	+	Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
495	+	if (!passDZcut)
496	+	continue;
497	+	}
498	+
499		// add good muon
500		CleanMuons->Add(mu);
501		}
#	Line 136 | Line 513 | void MuonIDMod::SlaveBegin()
513		// Run startup code on the computer (slave) doing the actual analysis. Here,
514		// we just request the muon collection branch.
515
516	<	ReqBranch(fMuonBranchName, fMuons);
516	>	// In this case we cannot have a branch
517	>	if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
518	>	ReqEventObject(fMuonBranchName, fMuons, kTRUE);
519	>	}
520	>	ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
521	>	ReqEventObject(fTrackName, fTracks, kTRUE);
522	>	ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
523	>	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
524	>	|| fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
525	>	|| fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
526	>	|| fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
527	>	|| fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
528	>	) {
529	>	ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
530	>	}
531
141	–	fMuonTools = new MuonTools;
532
533	<	if (fMuonIDType.CompareTo("Tight") == 0)
533	>	if (fMuonIDType.CompareTo("WMuId") == 0)
534	>	fMuIDType = kWMuId;
535	>	else if (fMuonIDType.CompareTo("ZMuId") == 0)
536	>	fMuIDType = kZMuId;
537	>	else if (fMuonIDType.CompareTo("Tight") == 0)
538		fMuIDType = kTight;
539		else if (fMuonIDType.CompareTo("Loose") == 0)
540		fMuIDType = kLoose;
541	+	else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
542	+	fMuIDType = kWWMuIdV1;
543	+	else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
544	+	fMuIDType = kWWMuIdV2;
545	+	else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
546	+	fMuIDType = kWWMuIdV3;
547		else if (fMuonIDType.CompareTo("NoId") == 0)
548		fMuIDType = kNoId;
549		else if (fMuonIDType.CompareTo("Custom") == 0) {
550		fMuIDType = kCustomId;
551		SendError(kWarning, "SlaveBegin",
552		"Custom muon identification is not yet implemented.");
553	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
554	+	fMuIDType = kMVAID_BDTG_IDIso;
555		} else {
556		SendError(kAbortAnalysis, "SlaveBegin",
557		"The specified muon identification %s is not defined.",
#	Line 163 | Line 565 | void MuonIDMod::SlaveBegin()
565		fMuIsoType = kTrackCaloCombined;
566		else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
567		fMuIsoType = kTrackCaloSliding;
568	+	else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
569	+	fMuIsoType = kTrackCaloSlidingNoCorrection;
570	+	else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
571	+	fMuIsoType = kCombinedRelativeConeAreaCorrected;
572	+	else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
573	+	fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
574	+	else if (fMuonIsoType.CompareTo("PFIso") == 0)
575	+	fMuIsoType = kPFIso;
576	+	else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
577	+	fMuIsoType = kPFRadialIso;
578	+	else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
579	+	fMuIsoType = kPFIsoEffectiveAreaCorrected;
580	+	else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
581	+	fMuIsoType = kPFIsoNoL;
582		else if (fMuonIsoType.CompareTo("NoIso") == 0)
583		fMuIsoType = kNoIso;
584		else if (fMuonIsoType.CompareTo("Custom") == 0) {
585		fMuIsoType = kCustomIso;
586		SendError(kWarning, "SlaveBegin",
587		"Custom muon isolation is not yet implemented.");
588	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
589	+	fMuIsoType = kMVAIso_BDTG_IDIso;
590		} else {
591		SendError(kAbortAnalysis, "SlaveBegin",
592		"The specified muon isolation %s is not defined.",
#	Line 180 | Line 598 | void MuonIDMod::SlaveBegin()
598		fMuClassType = kAll;
599		else if (fMuonClassType.CompareTo("Global") == 0)
600		fMuClassType = kGlobal;
601	+	else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
602	+	fMuClassType = kGlobalTracker;
603		else if (fMuonClassType.CompareTo("Standalone") == 0)
604		fMuClassType = kSta;
605	<	else if (fMuonClassType.CompareTo("TrackerOnly") == 0)
606	<	fMuClassType = kTrackerOnly;
605	>	else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
606	>	fMuClassType = kTrackerMuon;
607	>	else if (fMuonClassType.CompareTo("CaloMuon") == 0)
608	>	fMuClassType = kCaloMuon;
609	>	else if (fMuonClassType.CompareTo("TrackerBased") == 0)
610	>	fMuClassType = kTrackerBased;
611		else {
612		SendError(kAbortAnalysis, "SlaveBegin",
613		"The specified muon class %s is not defined.",
614		fMuonClassType.Data());
615		return;
616		}
617	+
618	+
619	+	//If we use MVA ID, need to load MVA weights
620	+	if(fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
621	+	fMuonTools = new MuonTools();
622	+	fMuonIDMVA = new MuonIDMVA();
623	+	fMuonIDMVA->Initialize("BDTG method",
624	+	fMuonMVAWeights_Subdet0Pt10To14p5,
625	+	fMuonMVAWeights_Subdet1Pt10To14p5,
626	+	fMuonMVAWeights_Subdet0Pt14p5To20,
627	+	fMuonMVAWeights_Subdet1Pt14p5To20,
628	+	fMuonMVAWeights_Subdet0Pt20ToInf,
629	+	fMuonMVAWeights_Subdet1Pt20ToInf,
630	+	MuonIDMVA::kIDIsoCombinedDetIso);
631	+	}
632	+
633	+	}
634	+
635	+
636	+	//--------------------------------------------------------------------------------------------------
637	+	Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
638	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
639	+	{
640	+
641	+	const Track *muTrk=0;
642	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
643	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
644	+
645	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
646	+
647	+	Int_t subdet = 0;
648	+	if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
649	+	else subdet = 1;
650	+	Int_t ptBin = 0;
651	+	if (muTrk->Pt() > 14.5) ptBin = 1;
652	+	if (muTrk->Pt() > 20.0) ptBin = 2;
653	+
654	+	Int_t MVABin = -1;
655	+	if      (subdet == 0 && ptBin == 0) MVABin = 0;
656	+	else if (subdet == 1 && ptBin == 0) MVABin = 1;
657	+	else if (subdet == 0 && ptBin == 1) MVABin = 2;
658	+	else if (subdet == 1 && ptBin == 1) MVABin = 3;
659	+	else if (subdet == 0 && ptBin == 2) MVABin = 4;
660	+	else if (subdet == 1 && ptBin == 2) MVABin = 5;
661	+
662	+	Double_t MVACut = -999;
663	+	if      (MVABin == 0) MVACut = -0.5618;
664	+	else if (MVABin == 1) MVACut = -0.3002;
665	+	else if (MVABin == 2) MVACut = -0.4642;
666	+	else if (MVABin == 3) MVACut = -0.2478;
667	+	else if (MVABin == 4) MVACut =  0.1706;
668	+	else if (MVABin == 5) MVACut =  0.8146;
669	+
670	+	if (MVAValue > MVACut) return kTRUE;
671	+	return kFALSE;
672		}

Diff Legend

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