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.3 by ceballos, Wed Nov 5 14:06:09 2008 UTC vs.
Revision 1.87 by mingyang, Mon Jan 7 22:14:56 2013 UTC

#	Line 1 | Line 1
1		// $Id$
2
3		#include "MitPhysics/Mods/interface/MuonIDMod.h"
4	–	#include "MitAna/DataTree/interface/Names.h"
5	–	#include "MitAna/DataCont/interface/ObjArray.h"
6	–	#include "MitPhysics/Utils/interface/IsolationTools.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;
11
#	Line 13 | Line 14 | ClassImp(mithep::MuonIDMod)
14		//--------------------------------------------------------------------------------------------------
15		MuonIDMod::MuonIDMod(const char *name, const char *title) :
16		BaseMod(name,title),
17	<	fPrintDebug(false),
18	<	fMuonName(Names::gkMuonBrn),
19	<	fCleanMuonsName(Names::gkCleanMuonsName),
20	<	fMuonIDType("Tight"),
21	<	fMuonIsoType("TrackCalo"),
22	<	fMuons(0),
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	>	fPFNoPileUpName("PFNoPileUp"),
27	>	fPFPileUpName("PFPileUp"),
28	>	fMuonIDType("WWMuIdV3"),
29	>	fMuonIsoType("PFIso"),
30	>	fMuonClassType("Global"),
31		fTrackIsolationCut(3.0),
32		fCaloIsolationCut(3.0),
33	+	fCombIsolationCut(0.15),
34	+	fCombRelativeIsolationCut(0.15),
35	+	fPFIsolationCut(-999.0),
36		fMuonPtMin(10),
37	<	fNEventsProcessed(0)
37	>	fApplyD0Cut(kTRUE),
38	>	fApplyDZCut(kTRUE),
39	>	fD0Cut(0.020),
40	>	fDZCut(0.10),
41	>	fWhichVertex(-1),
42	>	fEtaCut(2.4),
43	>	fMuIDType(kIdUndef),
44	>	fMuIsoType(kIsoUndef),
45	>	fMuClassType(kClassUndef),
46	>	fMuons(0),
47	>	fVertices(0),
48	>	fBeamSpot(0),
49	>	fTracks(0),
50	>	fPFCandidates(0),
51	>	fPFNoPileUpCands(0),
52	>	fPFPileUpCands(0),
53	>	fIntRadius(0.0),
54	>	fNonIsolatedMuons(0),
55	>	fNonIsolatedElectrons(0),
56	>	fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
57	>	fPileupEnergyDensity(0),
58	>	fMuonTools(0),
59	>	fMuonIDMVA(0),
60	>	fTheRhoType(RhoUtilities::DEFAULT),
61	>	fPVName(Names::gkPVBeamSpotBrn)
62		{
63		// Constructor.
64		}
65
66		//--------------------------------------------------------------------------------------------------
31	–	void MuonIDMod::Begin()
32	–	{
33	–	// Run startup code on the client machine. For this module, we dont do
34	–	// anything here.
35	–	}
36	–
37	–	//--------------------------------------------------------------------------------------------------
67		void MuonIDMod::Process()
68		{
69		// Process entries of the tree.
70
71	<	fNEventsProcessed++;
72	<
73	<	if (fNEventsProcessed % 1000 == 0 || fPrintDebug) {
74	<	time_t systime;
75	<	systime = time(NULL);
76	<
77	<	cerr << endl << "MuonIDMod : Process Event " << fNEventsProcessed << "  Time: " << ctime(&systime) << endl;
78	<	}
79	<
80	<	//Get Muons
81	<	LoadBranch(fMuonName);
82	<	ObjArray<Muon> *CleanMuons = new ObjArray<Muon>;
83	<	for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
84	<	Muon *mu = fMuons->At(i);
85	<
86	<	Double_t MuonClass = -1;
87	<	if (mu->GlobalTrk())
88	<	MuonClass = 0;
89	<	else if (mu->StandaloneTrk())
90	<	MuonClass = 1;
91	<	else if (mu->TrackerTrk())
92	<	MuonClass = 2;
93	<
94	<	bool allCuts = false;
95	<
96	<	if(MuonClass == 0) allCuts = true;
97	<
98	<	if(mu->IsoR03SumPt() >= fTrackIsolationCut) allCuts = false;
99	<
100	<	if(mu->IsoR03EmEt() +
101	<	mu->IsoR03HadEt() >= fCaloIsolationCut) allCuts = false;
102	<
103	<	if(mu->Pt() <= fMuonPtMin) allCuts = false;
104	<
105	<	if(allCuts) {
106	<	CleanMuons->Add(mu);
71	>	if(fCleanMuonsName.CompareTo("HggLeptonTagMuons") == 0 ){
72	>	LoadEventObject(fPVName,fVertices);
73	>	}
74	>	else{
75	>	fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
76	>	}
77	>
78	>	if(fMuIsoType != kPFIsoNoL) {
79	>	LoadEventObject(fMuonBranchName, fMuons);
80	>	}
81	>	else {
82	>	fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
83	>	}
84	>	LoadEventObject(fBeamSpotName, fBeamSpot);
85	>	LoadEventObject(fTrackName, fTracks);
86	>	LoadEventObject(fPFCandidatesName, fPFCandidates);
87	>	if(fMuIsoType == kTrackCaloSliding ||
88	>	fMuIsoType == kCombinedRelativeConeAreaCorrected ||
89	>	fMuIsoType == kPFIsoEffectiveAreaCorrected ||
90	>	fMuIsoType == kMVAIso_BDTG_IDIso ||
91	>	fMuIsoType == kIsoRingsV0_BDTG_Iso ||
92	>	fMuIsoType == kIsoDeltaR
93	>	) {
94	>	LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
95	>	}
96	>	if(fMuIsoType == kPFRadialIso || fMuIsoType == kIsoDeltaR || fMuIsoType == kPFIsoBetaPUCorrected){
97	>	// Name is hardcoded, can be changed if someone feels to do it *** did it--Heng
98	>	fPFNoPileUpCands = GetObjThisEvt<PFCandidateCol>(fPFNoPileUpName);
99	>	fPFPileUpCands = GetObjThisEvt<PFCandidateCol>(fPFPileUpName);
100	>	}
101	>
102	>	MuonOArr *CleanMuons = new MuonOArr;
103	>	CleanMuons->SetName(fCleanMuonsName);
104	>
105	>	for (UInt_t i=0; i<fMuons->GetEntries() && fVertices->GetEntries() > 0 ; ++i) {
106	>	const Muon *mu = fMuons->At(i);
107	>
108	>	Bool_t pass = kFALSE;
109	>	Double_t pt = 0;  // make sure pt is taken from the correct track!
110	>	Double_t eta = 0; // make sure eta is taken from the correct track!
111	>	switch (fMuClassType) {
112	>	case kAll:
113	>	pass = kTRUE;
114	>	if (mu->HasTrk()) {
115	>	pt  = mu->Pt();
116	>	eta = TMath::Abs(mu->Eta());
117	>	}
118	>	break;
119	>	case kGlobal:
120	>	pass = mu->HasGlobalTrk() && mu->IsTrackerMuon();
121	>	if (pass && mu->TrackerTrk()) {
122	>	pt  = mu->TrackerTrk()->Pt();
123	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
124	>	}
125	>	else {
126	>	pt  = mu->Pt();
127	>	eta = TMath::Abs(mu->Eta());
128	>	}
129	>	break;
130	>	case kGlobalTracker:
131	>	pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
132	>	(mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
133	>	(mu->IsTrackerMuon() &&
134	>	mu->Quality().Quality(MuonQuality::TMLastStationTight));
135	>	if (pass) {
136	>	pt  = mu->TrackerTrk()->Pt();
137	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
138	>	}
139	>	else {
140	>	pt  = mu->Pt();
141	>	eta = TMath::Abs(mu->Eta());
142	>	}
143	>	break;
144	>	case kSta:
145	>	pass = mu->HasStandaloneTrk();
146	>	if (pass) {
147	>	pt  = mu->StandaloneTrk()->Pt();
148	>	eta = TMath::Abs(mu->StandaloneTrk()->Eta());
149	>	}
150	>	break;
151	>	case kTrackerMuon:
152	>	pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
153	>	mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
154	>	if (pass) {
155	>	pt  = mu->TrackerTrk()->Pt();
156	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
157	>	}
158	>	break;
159	>	case kCaloMuon:
160	>	pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
161	>	if (pass) {
162	>	pt  = mu->TrackerTrk()->Pt();
163	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
164	>	}
165	>	break;
166	>	case kTrackerBased:
167	>	pass = mu->HasTrackerTrk();
168	>	if (pass) {
169	>	pt  = mu->TrackerTrk()->Pt();
170	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
171	>	}
172	>	break;
173	>	case kGlobalOnly:
174	>	pass = mu->HasGlobalTrk();
175	>	if (pass && mu->TrackerTrk()) {
176	>	pt  = mu->TrackerTrk()->Pt();
177	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
178	>	}
179	>	else {
180	>	pt  = mu->Pt();
181	>	eta = TMath::Abs(mu->Eta());
182	>	}
183	>	break;
184	>	default:
185	>	break;
186	>	}
187	>
188	>	if (!pass)
189	>	continue;
190	>
191	>	if (pt <= fMuonPtMin)
192	>	continue;
193	>
194	>	if (eta >= fEtaCut)
195	>	continue;
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 = 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 = mu->BestTrk() != 0 &&
231	>	mu->NTrkLayersHit() > 5 &&
232	>	mu->IsPFMuon() == kTRUE &&
233	>	mu->BestTrk()->NPixelHits() > 0 &&
234	>	RChi2 < 10.0;
235	>	break;
236	>	// 2012 WW analysis for 42x (there is no PFMuon link)
237	>	case kWWMuIdV1:
238	>	idpass = mu->BestTrk() != 0 &&
239	>	mu->NTrkLayersHit() > 5 &&
240	>	mu->BestTrk()->NPixelHits() > 0 &&
241	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
242	>	mu->TrkKink() < 20.0;
243	>	break;
244	>	// 2010 WW analysis
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	>	// 2011 WW analysis
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	>	// 2012 WW analysis
260	>	case kWWMuIdV4:
261	>	idpass = mu->BestTrk() != 0 &&
262	>	mu->NTrkLayersHit() > 5 &&
263	>	mu->IsPFMuon() == kTRUE &&
264	>	mu->BestTrk()->NPixelHits() > 0 &&
265	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
266	>	mu->TrkKink() < 20.0;
267	>	break;
268	>	case kMVAID_BDTG_IDIso:
269	>	{
270	>	Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
271	>	mu->BestTrk()->NHits() > 10 &&
272	>	mu->BestTrk()->NPixelHits() > 0 &&
273	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
274	>	MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
275	>	MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
276	>	mu->TrkKink() < 20.0
277	>	);
278	>	idpass =  passDenominatorM2;
279	>	//only evaluate MVA if muon passes M2 denominator to save time
280	>	if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
281	>	}
282	>	break;
283	>	case kNoId:
284	>	idpass = kTRUE;
285	>	break;
286	>	default:
287	>	break;
288		}
289	+
290	+	if (!idpass)
291	+	continue;
292	+
293	+	Double_t Rho = 0.;
294	+	if(fPileupEnergyDensity){
295	+	const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
296	+
297	+	switch (fTheRhoType) {
298	+	case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
299	+	Rho = rho->RhoLowEta();
300	+	break;
301	+	case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
302	+	Rho = rho->Rho();
303	+	break;
304	+	case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
305	+	Rho = rho->RhoRandomLowEta();
306	+	break;
307	+	case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
308	+	Rho = rho->RhoRandom();
309	+	break;
310	+	case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
311	+	Rho = rho->RhoKt6PFJets();
312	+	break;
313	+	default:
314	+	Rho = rho->Rho();
315	+	}
316	+
317	+	if ((TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho()))) Rho = 0.;
318	+	}
319	+
320	+	Bool_t isocut = kFALSE;
321	+	switch (fMuIsoType) {
322	+	case kTrackCalo:
323	+	isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
324	+	(mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
325	+	break;
326	+	case kTrackCaloCombined:
327	+	isocut = (1.0 * mu->IsoR03SumPt() +
328	+	1.0 * mu->IsoR03EmEt()  +
329	+	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
330	+	break;
331	+	case kTrackCaloSliding:
332	+	{
333	+	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
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 kTrackCaloSlidingNoCorrection:
344	+	{
345	+	Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
346	+	1.0 * mu->IsoR03EmEt()  +
347	+	1.0 * mu->IsoR03HadEt();
348	+	// trick to change the signal region cut
349	+	double theIsoCut = fCombIsolationCut;
350	+	if(theIsoCut < 0.20){
351	+	if(mu->Pt() >  20.0) theIsoCut = 0.15;
352	+	else                 theIsoCut = 0.10;
353	+	}
354	+	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
355	+	}
356	+	break;
357	+	case kCombinedRelativeConeAreaCorrected:
358	+	{
359	+	//const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0); // Fabian: made Rho customable
360	+	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
361	+	double theIsoCut = fCombRelativeIsolationCut;
362	+	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
363	+	}
364	+	break;
365	+	case kCombinedRelativeEffectiveAreaCorrected:
366	+	{
367	+	Double_t tmpRho = Rho;   // Fabian: made the Rho type customable.
368	+	//if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))
369	+	//tmpRho = fPileupEnergyDensity->At(0)->Rho();
370	+
371	+	isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
372	+	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
373	+	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
374	+	) < (mu->Pt()* 0.40);
375	+	}
376	+	break;
377	+	case kPFIso:
378	+	{
379	+	Double_t pfIsoCutValue = 9999;
380	+	if(fPFIsolationCut > 0){
381	+	pfIsoCutValue = fPFIsolationCut;
382	+	} else {
383	+	if (mu->AbsEta() < 1.479) {
384	+	if (mu->Pt() > 20) {
385	+	pfIsoCutValue = 0.13;
386	+	} else {
387	+	pfIsoCutValue = 0.06;
388	+	}
389	+	} else {
390	+	if (mu->Pt() > 20) {
391	+	pfIsoCutValue = 0.09;
392	+	} else {
393	+	pfIsoCutValue = 0.05;
394	+	}
395	+	}
396	+	}
397	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
398	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
399	+	isocut = kTRUE;
400	+	}
401	+	break;
402	+	case kPFRadialIso:
403	+	{
404	+	Double_t pfIsoCutValue = 9999;
405	+	if(fPFIsolationCut > 0){
406	+	pfIsoCutValue = fPFIsolationCut;
407	+	} else {
408	+	if (mu->Pt() > 20) {
409	+	pfIsoCutValue = 0.10;
410	+	} else {
411	+	pfIsoCutValue = 0.05;
412	+	}
413	+	}
414	+	Double_t totalIso =  IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
415	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
416	+	isocut = kTRUE;
417	+	}
418	+	break;
419	+	case kPFIsoBetaPUCorrected:
420	+	{
421	+	Double_t pfIsoCutValue = 9999;
422	+	if(fPFIsolationCut > 0){
423	+	pfIsoCutValue = fPFIsolationCut;
424	+	} else {
425	+	if (mu->Pt() > 20) {
426	+	pfIsoCutValue = 0.2;
427	+	} else {
428	+	pfIsoCutValue = 0.2;
429	+	}
430	+	}
431	+	Double_t totalIso =  IsolationTools::BetaMwithPUCorrection(fPFNoPileUpCands, fPFPileUpCands, mu, 0.4);
432	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
433	+	isocut = kTRUE;
434	+	}
435	+	break;
436	+	case kPFIsoEffectiveAreaCorrected:
437	+	{
438	+	Double_t pfIsoCutValue = 9999;
439	+	if(fPFIsolationCut > 0){
440	+	pfIsoCutValue = fPFIsolationCut;
441	+	} else {
442	+	pfIsoCutValue = fPFIsolationCut; //leave it like this for now
443	+	}
444	+	Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
445	+	- Rho * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
446	+	//- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());  // Fabian: made Rho-type customable
447	+	isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
448	+	break;
449	+	}
450	+
451	+
452	+	case kPFIsoNoL:
453	+	{
454	+	fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
455	+	fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
456	+
457	+	Double_t pfIsoCutValue = 9999;
458	+	if(fPFIsolationCut > 0){
459	+	pfIsoCutValue = fPFIsolationCut;
460	+	} else {
461	+	if (mu->AbsEta() < 1.479) {
462	+	if (mu->Pt() > 20) {
463	+	pfIsoCutValue = 0.13;
464	+	} else {
465	+	pfIsoCutValue = 0.06;
466	+	}
467	+	} else {
468	+	if (mu->Pt() > 20) {
469	+	pfIsoCutValue = 0.09;
470	+	} else {
471	+	pfIsoCutValue = 0.05;
472	+	}
473	+	}
474	+	}
475	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
476	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
477	+	isocut = kTRUE;
478	+	}
479	+	break;
480	+	case kMVAIso_BDTG_IDIso:
481	+	{
482	+
483	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
484	+	isocut = (totalIso < (mu->Pt()*0.4));
485	+
486	+	}
487	+	break;
488	+	case kIsoRingsV0_BDTG_Iso:
489	+	{
490	+
491	+	isocut = PassMuonIsoRingsV0_BDTG_Iso(mu, fVertices->At(0), fPileupEnergyDensity);
492	+
493	+	}
494	+	break;
495	+	case kIsoDeltaR:
496	+	{
497	+
498	+	isocut = PassMuonIsoDeltaR(mu, fVertices->At(0), fPileupEnergyDensity);
499	+
500	+	}
501	+	break;
502	+	case kNoIso:
503	+	isocut = kTRUE;
504	+	break;
505	+	case kCustomIso:
506	+	default:
507	+	break;
508	+	}
509	+
510	+	if (isocut == kFALSE)
511	+	continue;
512	+
513	+	// apply d0 cut
514	+	if (fApplyD0Cut) {
515	+	Bool_t passD0cut = kTRUE;
516	+	if(fD0Cut < 0.05) { // trick to change the signal region cut
517	+	if      (mu->Pt() >  20.0) fD0Cut = 0.02;
518	+	else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
519	+	}
520	+	if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
521	+	else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
522	+	if (!passD0cut)
523	+	continue;
524	+	}
525	+
526	+	// apply dz cut
527	+	if (fApplyDZCut) {
528	+	Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
529	+	if (!passDZcut)
530	+	continue;
531	+	}
532	+
533	+	// add good muon
534	+	CleanMuons->Add(mu);
535		}
536
537	<	//Final Summary Debug Output
538	<	if ( fPrintDebug ) {
83	<	cerr << "Event Dump: " << fNEventsProcessed << endl;
84	<	cerr << "Muons" << endl;
85	<	for (UInt_t i = 0; i < CleanMuons->GetEntries(); i++) {
86	<	cerr << i << " " << CleanMuons->At(i)->Pt() << " " << CleanMuons->At(i)->Eta()
87	<	<< " " << CleanMuons->At(i)->Phi() << endl;
88	<	}
89	<	}
90	<
91	<	//Save Objects for Other Modules to use
92	<	AddObjThisEvt(CleanMuons, fCleanMuonsName.Data());
93	<	}
537	>	// sort according to pt
538	>	CleanMuons->Sort();
539
540	+	// add objects for other modules to use
541	+	AddObjThisEvt(CleanMuons);
542	+	}
543
544		//--------------------------------------------------------------------------------------------------
545		void MuonIDMod::SlaveBegin()
546		{
547		// Run startup code on the computer (slave) doing the actual analysis. Here,
548	<	// we typically initialize histograms and other analysis objects and request
549	<	// branches. For this module, we request a branch of the MitTree.
548	>	// we just request the muon collection branch.
549	>
550	>	if(fCleanMuonsName.CompareTo("HggLeptonTagMuons") == 0 ){
551	>	ReqEventObject(fPVName,fVertices,true);
552	>	}
553	>
554	>	// In this case we cannot have a branch
555	>	if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
556	>	ReqEventObject(fMuonBranchName, fMuons, kTRUE);
557	>	}
558	>	ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
559	>	ReqEventObject(fTrackName, fTracks, kTRUE);
560	>	ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
561	>	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
562	>	|| fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
563	>	|| fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
564	>	|| fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
565	>	|| fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
566	>	|| fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0
567	>	|| fMuonIsoType.CompareTo("IsoDeltaR") == 0
568	>	) {
569	>	ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
570	>	}
571	>
572	>
573	>	if (fMuonIDType.CompareTo("WMuId") == 0)
574	>	fMuIDType = kWMuId;
575	>	else if (fMuonIDType.CompareTo("ZMuId") == 0)
576	>	fMuIDType = kZMuId;
577	>	else if (fMuonIDType.CompareTo("Tight") == 0)
578	>	fMuIDType = kTight;
579	>	else if (fMuonIDType.CompareTo("Loose") == 0)
580	>	fMuIDType = kLoose;
581	>	else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
582	>	fMuIDType = kWWMuIdV1;
583	>	else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
584	>	fMuIDType = kWWMuIdV2;
585	>	else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
586	>	fMuIDType = kWWMuIdV3;
587	>	else if (fMuonIDType.CompareTo("WWMuIdV4") == 0)
588	>	fMuIDType = kWWMuIdV4;
589	>	else if (fMuonIDType.CompareTo("NoId") == 0)
590	>	fMuIDType = kNoId;
591	>	else if (fMuonIDType.CompareTo("Custom") == 0) {
592	>	fMuIDType = kCustomId;
593	>	SendError(kWarning, "SlaveBegin",
594	>	"Custom muon identification is not yet implemented.");
595	>	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
596	>	fMuIDType = kMVAID_BDTG_IDIso;
597	>	} else {
598	>	SendError(kAbortAnalysis, "SlaveBegin",
599	>	"The specified muon identification %s is not defined.",
600	>	fMuonIDType.Data());
601	>	return;
602	>	}
603	>
604	>	if (fMuonIsoType.CompareTo("TrackCalo") == 0)
605	>	fMuIsoType = kTrackCalo;
606	>	else if (fMuonIsoType.CompareTo("TrackCaloCombined") == 0)
607	>	fMuIsoType = kTrackCaloCombined;
608	>	else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
609	>	fMuIsoType = kTrackCaloSliding;
610	>	else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
611	>	fMuIsoType = kTrackCaloSlidingNoCorrection;
612	>	else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
613	>	fMuIsoType = kCombinedRelativeConeAreaCorrected;
614	>	else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
615	>	fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
616	>	else if (fMuonIsoType.CompareTo("PFIso") == 0)
617	>	fMuIsoType = kPFIso;
618	>	else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
619	>	fMuIsoType = kPFRadialIso;
620	>	else if (fMuonIsoType.CompareTo("PFIsoBetaPUCorrected") == 0)
621	>	fMuIsoType = kPFIsoBetaPUCorrected;
622	>	else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
623	>	fMuIsoType = kPFIsoEffectiveAreaCorrected;
624	>	else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
625	>	fMuIsoType = kPFIsoNoL;
626	>	else if (fMuonIsoType.CompareTo("NoIso") == 0)
627	>	fMuIsoType = kNoIso;
628	>	else if (fMuonIsoType.CompareTo("Custom") == 0) {
629	>	fMuIsoType = kCustomIso;
630	>	SendError(kWarning, "SlaveBegin",
631	>	"Custom muon isolation is not yet implemented.");
632	>	} else if (fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0) {
633	>	fMuIsoType = kMVAIso_BDTG_IDIso;
634	>	} else if (fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0) {
635	>	fMuIsoType = kIsoRingsV0_BDTG_Iso;
636	>	} else if (fMuonIsoType.CompareTo("IsoDeltaR") == 0) {
637	>	fMuIsoType = kIsoDeltaR;
638	>	} else {
639	>	SendError(kAbortAnalysis, "SlaveBegin",
640	>	"The specified muon isolation %s is not defined.",
641	>	fMuonIsoType.Data());
642	>	return;
643	>	}
644	>
645	>	if (fMuonClassType.CompareTo("All") == 0)
646	>	fMuClassType = kAll;
647	>	else if (fMuonClassType.CompareTo("Global") == 0)
648	>	fMuClassType = kGlobal;
649	>	else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
650	>	fMuClassType = kGlobalTracker;
651	>	else if (fMuonClassType.CompareTo("Standalone") == 0)
652	>	fMuClassType = kSta;
653	>	else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
654	>	fMuClassType = kTrackerMuon;
655	>	else if (fMuonClassType.CompareTo("CaloMuon") == 0)
656	>	fMuClassType = kCaloMuon;
657	>	else if (fMuonClassType.CompareTo("TrackerBased") == 0)
658	>	fMuClassType = kTrackerBased;
659	>	else if (fMuonClassType.CompareTo("GlobalOnly") == 0)
660	>	fMuClassType = kGlobalOnly;
661	>	else {
662	>	SendError(kAbortAnalysis, "SlaveBegin",
663	>	"The specified muon class %s is not defined.",
664	>	fMuonClassType.Data());
665	>	return;
666	>	}
667	>
668	>
669	>	//If we use MVA ID, need to load MVA weights
670	>	if     (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
671	>	fMuonTools = new MuonTools();
672	>	fMuonIDMVA = new MuonIDMVA();
673	>	fMuonIDMVA->Initialize("BDTG method",
674	>	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin0_IDIsoCombined_BDTG.weights.xml"))),
675	>	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin0_IDIsoCombined_BDTG.weights.xml"))),
676	>	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin1_IDIsoCombined_BDTG.weights.xml"))),
677	>	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin1_IDIsoCombined_BDTG.weights.xml"))),
678	>	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin2_IDIsoCombined_BDTG.weights.xml"))),
679	>	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin2_IDIsoCombined_BDTG.weights.xml"))),
680	>	MuonIDMVA::kIDIsoCombinedDetIso,
681	>	fTheRhoType);
682	>	}
683	>	else if(fMuIsoType == kIsoRingsV0_BDTG_Iso) {
684	>	std::vector<std::string> muonidiso_weightfiles;
685	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml"))));
686	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml"))));
687	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml"))));
688	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml"))));
689	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml"))));
690	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml"))));
691	>	fMuonTools = new MuonTools();
692	>	fMuonIDMVA = new MuonIDMVA();
693	>	fMuonIDMVA->Initialize("MuonIso_BDTG_IsoRings",
694	>	MuonIDMVA::kIsoRingsV0,
695	>	kTRUE,
696	>	muonidiso_weightfiles,
697	>	fTheRhoType);
698	>	}
699	>	else if(fMuIsoType == kIsoDeltaR) {
700	>	std::vector<std::string> muonidiso_weightfiles;
701	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LB_BDT.weights.xml"))));
702	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LE_BDT.weights.xml"))));
703	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HB_BDT.weights.xml"))));
704	>	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HE_BDT.weights.xml"))));
705	>	fMuonTools = new MuonTools();
706	>	fMuonIDMVA = new MuonIDMVA();
707	>	fMuonIDMVA->Initialize("muonHZZ2012IsoDRMVA",
708	>	MuonIDMVA::kIsoDeltaR,
709	>	kTRUE,
710	>	muonidiso_weightfiles,
711	>	fTheRhoType);
712	>	}
713	>
714	>	}
715	>
716	>
717	>	//--------------------------------------------------------------------------------------------------
718	>	Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
719	>	const PileupEnergyDensityCol *PileupEnergyDensity) const
720	>	{
721
722	<	ReqBranch(fMuonName,              fMuons);
722	>	const Track *muTrk=0;
723	>	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
724	>	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
725	>
726	>	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
727	>
728	>	Int_t subdet = 0;
729	>	if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
730	>	else subdet = 1;
731	>	Int_t ptBin = 0;
732	>	if (muTrk->Pt() > 14.5) ptBin = 1;
733	>	if (muTrk->Pt() > 20.0) ptBin = 2;
734	>
735	>	Int_t MVABin = -1;
736	>	if      (subdet == 0 && ptBin == 0) MVABin = 0;
737	>	else if (subdet == 1 && ptBin == 0) MVABin = 1;
738	>	else if (subdet == 0 && ptBin == 1) MVABin = 2;
739	>	else if (subdet == 1 && ptBin == 1) MVABin = 3;
740	>	else if (subdet == 0 && ptBin == 2) MVABin = 4;
741	>	else if (subdet == 1 && ptBin == 2) MVABin = 5;
742	>
743	>	Double_t MVACut = -999;
744	>	if      (MVABin == 0) MVACut = -0.5618;
745	>	else if (MVABin == 1) MVACut = -0.3002;
746	>	else if (MVABin == 2) MVACut = -0.4642;
747	>	else if (MVABin == 3) MVACut = -0.2478;
748	>	else if (MVABin == 4) MVACut =  0.1706;
749	>	else if (MVABin == 5) MVACut =  0.8146;
750	>
751	>	if (MVAValue > MVACut) return kTRUE;
752	>	return kFALSE;
753		}
754
755		//--------------------------------------------------------------------------------------------------
756	<	void MuonIDMod::SlaveTerminate()
756	>	Bool_t MuonIDMod::PassMuonIsoRingsV0_BDTG_Iso(const Muon *mu, const Vertex *vertex,
757	>	const PileupEnergyDensityCol *PileupEnergyDensity) const
758		{
109	–	// Run finishing code on the computer (slave) that did the analysis. For this
110	–	// module, we dont do anything here.
759
760	+	Bool_t isDebug = kFALSE;
761	+	const Track *muTrk=0;
762	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
763	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
764	+
765	+	ElectronOArr *tempElectrons = new  ElectronOArr;
766	+	MuonOArr     *tempMuons     = new  MuonOArr;
767	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,
768	+	PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,isDebug);
769	+	delete tempElectrons;
770	+	delete tempMuons;
771	+
772	+	Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
773	+
774	+	Double_t MVACut = -1.0;
775	+	Double_t eta = mu->AbsEta();
776	+	if     (mu->Pt() <  20 && eta <  1.479) MVACut = 0.86;
777	+	else if(mu->Pt() <  20 && eta >= 1.479) MVACut = 0.82;
778	+	else if(mu->Pt() >= 20 && eta <  1.479) MVACut = 0.82;
779	+	else if(mu->Pt() >= 20 && eta >= 1.479) MVACut = 0.86;
780	+
781	+	if(fPFIsolationCut > -1.0) MVACut = fPFIsolationCut;
782	+
783	+	if(isDebug == kTRUE){
784	+	printf("PassMuonIsoRingsV0_BDTG_IsoDebug: %d, pt, eta = %f, %f, rho = %f(%f) : RingsMVA = %f, bin: %d\n",
785	+	GetEventHeader()->EvtNum(),mu->Pt(), mu->Eta(),
786	+	fPileupEnergyDensity->At(0)->Rho(),fPileupEnergyDensity->At(0)->RhoKt6PFJets(),MVAValue,MVABin);
787	+	}
788	+
789	+	if (MVAValue > MVACut) return kTRUE;
790	+	return kFALSE;
791	+	}
792	+
793	+	//--------------------------------------------------------------------------------------------------
794	+	Bool_t MuonIDMod::PassMuonIsoDeltaR(const Muon *mu, const Vertex *vertex,
795	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
796	+	{
797	+
798	+	const Track *muTrk=0;
799	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
800	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
801	+
802	+	ElectronOArr *tempElectrons = new  ElectronOArr;
803	+	MuonOArr     *tempMuons     = new  MuonOArr;
804	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFNoPileUpCands,
805	+	PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,kFALSE);
806	+	delete tempElectrons;
807	+	delete tempMuons;
808	+
809	+	Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
810	+
811	+	Double_t MVACut = -999;
812	+	if      (MVABin == 0) MVACut =  0.000;
813	+	else if (MVABin == 1) MVACut =  0.000;
814	+	else if (MVABin == 2) MVACut =  0.000;
815	+	else if (MVABin == 3) MVACut =  0.000;
816	+
817	+	if (MVAValue > MVACut) return kTRUE;
818	+	return kFALSE;
819		}
820
821		//--------------------------------------------------------------------------------------------------
822		void MuonIDMod::Terminate()
823		{
824	<	// Run finishing code on the client computer. For this module, we dont do
825	<	// anything here.
824	>	// Run finishing code on the computer (slave) that did the analysis
825	>	delete fMuonIDMVA;
826	>
827	>	delete fMuonTools;
828		}

Diff Legend

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