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

Diff Legend

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