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.57 by ceballos, Tue Oct 11 17:00:24 2011 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 13 | Line 16 | ClassImp(mithep::MuonIDMod)
16		BaseMod(name,title),
17		fMuonBranchName(Names::gkMuonBrn),
18		fCleanMuonsName(ModNames::gkCleanMuonsName),
19	<	fMuonIDType("Loose"),
20	<	fMuonIsoType("TrackCaloSliding"),
19	>	fNonIsolatedMuonsName("random"),
20	>	fNonIsolatedElectronsName("random"),
21	>	fVertexName(ModNames::gkGoodVertexesName),
22	>	fBeamSpotName(Names::gkBeamSpotBrn),
23	>	fTrackName(Names::gkTrackBrn),
24	>	fPFCandidatesName(Names::gkPFCandidatesBrn),
25	>	fMuonIDType("WWMuIdV3"),
26	>	fMuonIsoType("PFIso"),
27		fMuonClassType("Global"),
28		fTrackIsolationCut(3.0),
29		fCaloIsolationCut(3.0),
30	<	fCombIsolationCut(5.0),
30	>	fCombIsolationCut(0.15),
31	>	fCombRelativeIsolationCut(0.15),
32	>	fPFIsolationCut(-1.0),
33		fMuonPtMin(10),
34	<	fMuons(0)
34	>	fApplyD0Cut(kTRUE),
35	>	fApplyDZCut(kTRUE),
36	>	fD0Cut(0.020),
37	>	fDZCut(0.10),
38	>	fWhichVertex(-1),
39	>	fEtaCut(2.4),
40	>	fMuIDType(kIdUndef),
41	>	fMuIsoType(kIsoUndef),
42	>	fMuClassType(kClassUndef),
43	>	fMuons(0),
44	>	fVertices(0),
45	>	fBeamSpot(0),
46	>	fTracks(0),
47	>	fPFCandidates(0),
48	>	fIntRadius(0.0),
49	>	fNonIsolatedMuons(0),
50	>	fNonIsolatedElectrons(0),
51	>	fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
52	>	fPileupEnergyDensity(0)
53		{
54		// Constructor.
55		}
#	Line 30 | Line 59 | void MuonIDMod::Process()
59		{
60		// Process entries of the tree.
61
62	<	LoadBranch(fMuonBranchName);
63	<
62	>	if(fMuIsoType != kPFIsoNoL) {
63	>	LoadEventObject(fMuonBranchName, fMuons);
64	>	}
65	>	else {
66	>	fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
67	>	}
68	>	LoadEventObject(fBeamSpotName, fBeamSpot);
69	>	LoadEventObject(fTrackName, fTracks);
70	>	LoadEventObject(fPFCandidatesName, fPFCandidates);
71	>	if(fMuIsoType == kTrackCaloSliding || fMuIsoType == kCombinedRelativeConeAreaCorrected) {
72	>	LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
73	>	}
74		MuonOArr *CleanMuons = new MuonOArr;
75		CleanMuons->SetName(fCleanMuonsName);
76
77	+	fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
78	+
79		for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
80		const Muon *mu = fMuons->At(i);
81
82		Bool_t pass = kFALSE;
83	<	Double_t pt = -1; // make sure pt is taken from the correct track!
83	>	Double_t pt = 0;  // make sure pt is taken from the correct track!
84	>	Double_t eta = 0; // make sure eta is taken from the correct track!
85		switch (fMuClassType) {
86		case kAll:
87		pass = kTRUE;
88	<	pt = mu->Pt();
88	>	if (mu->HasTrk()) {
89	>	pt  = mu->Pt();
90	>	eta = TMath::Abs(mu->Eta());
91	>	}
92		break;
93		case kGlobal:
94	<	pass = (mu->GlobalTrk() != 0);
95	<	if (pass)
96	<	pt = mu->GlobalTrk()->Pt();
97	<	break;
94	>	pass = mu->HasGlobalTrk() && mu->IsTrackerMuon();
95	>	if (pass && mu->TrackerTrk()) {
96	>	pt  = mu->TrackerTrk()->Pt();
97	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
98	>	}
99	>	else {
100	>	pt  = mu->Pt();
101	>	eta = TMath::Abs(mu->Eta());
102	>	}
103	>	break;
104	>	case kGlobalTracker:
105	>	pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
106	>	(mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
107	>	(mu->IsTrackerMuon() &&
108	>	mu->Quality().Quality(MuonQuality::TMLastStationTight));
109	>	if (pass) {
110	>	pt  = mu->TrackerTrk()->Pt();
111	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
112	>	}
113	>	else {
114	>	pt  = mu->Pt();
115	>	eta = TMath::Abs(mu->Eta());
116	>	}
117	>	break;
118		case kSta:
119	<	pass = (mu->StandaloneTrk() != 0);
120	<	if (pass)
121	<	pt = mu->StandaloneTrk()->Pt();
122	<	break;
123	<	case kTrackerOnly:
124	<	pass = (mu->TrackerTrk() != 0);
125	<	if (pass)
126	<	pt = mu->TrackerTrk()->Pt();
119	>	pass = mu->HasStandaloneTrk();
120	>	if (pass) {
121	>	pt  = mu->StandaloneTrk()->Pt();
122	>	eta = TMath::Abs(mu->StandaloneTrk()->Eta());
123	>	}
124	>	break;
125	>	case kTrackerMuon:
126	>	pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
127	>	mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
128	>	if (pass) {
129	>	pt  = mu->TrackerTrk()->Pt();
130	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
131	>	}
132	>	break;
133	>	case kCaloMuon:
134	>	pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
135	>	if (pass) {
136	>	pt  = mu->TrackerTrk()->Pt();
137	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
138	>	}
139	>	break;
140	>	case kTrackerBased:
141	>	pass = mu->HasTrackerTrk();
142	>	if (pass) {
143	>	pt  = mu->TrackerTrk()->Pt();
144	>	eta = TMath::Abs(mu->TrackerTrk()->Eta());
145	>	}
146		break;
147		default:
148		break;
#	Line 70 | Line 154 | void MuonIDMod::Process()
154		if (pt <= fMuonPtMin)
155		continue;
156
157	+	if (eta >= fEtaCut)
158	+	continue;
159	+
160	+	Double_t RChi2 = 0.0;
161	+	if     (mu->HasGlobalTrk()) {
162	+	RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
163	+	}
164	+	else if(mu->BestTrk() != 0){
165	+	RChi2 = mu->BestTrk()->Chi2()/mu->BestTrk()->Ndof();
166	+	}
167		Bool_t idpass = kFALSE;
168		switch (fMuIDType) {
169	+	case kWMuId:
170	+	idpass = mu->BestTrk() != 0 &&
171	+	mu->BestTrk()->NHits() > 10 &&
172	+	RChi2 < 10.0 &&
173	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
174	+	mu->BestTrk()->NPixelHits() > 0 &&
175	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
176	+	break;
177	+	case kZMuId:
178	+	idpass = mu->BestTrk() != 0 &&
179	+	mu->BestTrk()->NHits() > 10 &&
180	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
181	+	mu->BestTrk()->NPixelHits() > 0 &&
182	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
183	+	break;
184		case kLoose:
185	<	idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationLoose) &&
186	<	fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityLoose);
185	>	idpass = mu->BestTrk() != 0 &&
186	>	mu->Quality().Quality(MuonQuality::TMOneStationLoose) &&
187	>	mu->Quality().Quality(MuonQuality::TM2DCompatibilityLoose) &&
188	>	mu->BestTrk()->NHits() > 10 &&
189	>	RChi2 < 10.0 &&
190	>	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
191		break;
192		case kTight:
193	<	idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationTight) &&
194	<	fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityTight);
193	>	idpass = mu->BestTrk() !=  0 &&
194	>	mu->Quality().Quality(MuonQuality::TMOneStationTight) &&
195	>	mu->Quality().Quality(MuonQuality::TM2DCompatibilityTight) &&
196	>	mu->BestTrk()->NHits() > 10 &&
197	>	RChi2 < 10.0 &&
198	>	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
199	>	break;
200	>	case kWWMuIdV1:
201	>	idpass = mu->BestTrk() != 0 &&
202	>	mu->BestTrk()->NHits() > 10 &&
203	>	mu->BestTrk()->NPixelHits() > 0 &&
204	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
205	>	RChi2 < 10.0 &&
206	>	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
207	>	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
208	>	break;
209	>	case kWWMuIdV2:
210	>	idpass = mu->BestTrk() != 0 &&
211	>	mu->BestTrk()->NHits() > 10 &&
212	>	mu->BestTrk()->NPixelHits() > 0 &&
213	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
214	>	break;
215	>	case kWWMuIdV3:
216	>	idpass = mu->BestTrk() != 0 &&
217	>	mu->BestTrk()->NHits() > 10 &&
218	>	mu->BestTrk()->NPixelHits() > 0 &&
219	>	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
220	>	mu->TrkKink() < 20.0;
221	>	break;
222	>	case kNoId:
223	>	idpass = kTRUE;
224		break;
225		default:
226		break;
#	Line 87 | Line 229 | void MuonIDMod::Process()
229		if (!idpass)
230		continue;
231
232	<	Bool_t isopass = kFALSE;
232	>	Bool_t isocut = kFALSE;
233		switch (fMuIsoType) {
234		case kTrackCalo:
235	<	isopass = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
235	>	isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
236		(mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
237		break;
238		case kTrackCaloCombined:
239	<	isopass = (1.0 * mu->IsoR03SumPt() + 1.0 * mu->IsoR03EmEt() +
240	<	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
239	>	isocut = (1.0 * mu->IsoR03SumPt() +
240	>	1.0 * mu->IsoR03EmEt()  +
241	>	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
242		break;
243		case kTrackCaloSliding:
244		{
245	<	Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
246	<	1.0 * mu->IsoR03EmEt() +
247	<	1.0 * mu->IsoR03HadEt();
248	<	if ((totalIso < (mu->Pt()-10.0)*5.0/15.0 && mu->Pt() <= 25) ||
249	<	(totalIso < 5.0 && mu->Pt() > 25) ||
250	<	totalIso <= 0)
251	<	isopass = kTRUE;
245	>	const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
246	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
247	>	// trick to change the signal region cut
248	>	double theIsoCut = fCombIsolationCut;
249	>	if(theIsoCut < 0.20){
250	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
251	>	else                 theIsoCut = 0.10;
252	>	}
253	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
254	>	}
255	>	break;
256	>	case kTrackCaloSlidingNoCorrection:
257	>	{
258	>	Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
259	>	1.0 * mu->IsoR03EmEt()  +
260	>	1.0 * mu->IsoR03HadEt();
261	>	// trick to change the signal region cut
262	>	double theIsoCut = fCombIsolationCut;
263	>	if(theIsoCut < 0.20){
264	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
265	>	else                 theIsoCut = 0.10;
266	>	}
267	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
268	>	}
269	>	break;
270	>	case kCombinedRelativeConeAreaCorrected:
271	>	{
272	>	const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
273	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
274	>	double theIsoCut = fCombRelativeIsolationCut;
275	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
276		}
277	+	break;
278	+	case kPFIso:
279	+	{
280	+	Double_t pfIsoCutValue = 9999;
281	+	if(fPFIsolationCut > 0){
282	+	pfIsoCutValue = fPFIsolationCut;
283	+	} else {
284	+	if (mu->AbsEta() < 1.479) {
285	+	if (mu->Pt() > 20) {
286	+	pfIsoCutValue = 0.13;
287	+	} else {
288	+	pfIsoCutValue = 0.06;
289	+	}
290	+	} else {
291	+	if (mu->Pt() > 20) {
292	+	pfIsoCutValue = 0.09;
293	+	} else {
294	+	pfIsoCutValue = 0.05;
295	+	}
296	+	}
297	+	}
298	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
299	+	if (totalIso < (mu->Pt()*pfIsoCutValue) )
300	+	isocut = kTRUE;
301	+	}
302	+	break;
303	+	case kPFIsoNoL:
304	+	{
305	+	fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
306	+	fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
307	+
308	+	Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
309	+	if(beta == 0) beta = 1.0;
310	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), fNonIsolatedMuons, fNonIsolatedElectrons, 0.2, 1.0, 0.4, 0.0, 3);
311	+	if (totalIso < (mu->Pt()*fPFIsolationCut) )
312	+	isocut = kTRUE;
313	+	}
314		break;
315		case kNoIso:
316	<	isopass = kTRUE;
316	>	isocut = kTRUE;
317		break;
318		case kCustomIso:
319		default:
320		break;
321		}
322
323	<	if (!isopass)
323	>	if (isocut == kFALSE)
324		continue;
325
326	+	// apply d0 cut
327	+	if (fApplyD0Cut) {
328	+	Bool_t passD0cut = kTRUE;
329	+	if(fD0Cut < 0.05) { // trick to change the signal region cut
330	+	if      (mu->Pt() >  20.0) fD0Cut = 0.02;
331	+	else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
332	+	}
333	+	if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
334	+	else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
335	+	if (!passD0cut)
336	+	continue;
337	+	}
338	+
339	+	// apply dz cut
340	+	if (fApplyDZCut) {
341	+	Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
342	+	if (!passDZcut)
343	+	continue;
344	+	}
345	+
346		// add good muon
347		CleanMuons->Add(mu);
348		}
#	Line 136 | Line 360 | void MuonIDMod::SlaveBegin()
360		// Run startup code on the computer (slave) doing the actual analysis. Here,
361		// we just request the muon collection branch.
362
363	<	ReqBranch(fMuonBranchName, fMuons);
364	<
365	<	fMuonTools = new MuonTools;
363	>	// In this case we cannot have a branch
364	>	if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
365	>	ReqEventObject(fMuonBranchName, fMuons, kTRUE);
366	>	}
367	>	ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
368	>	ReqEventObject(fTrackName, fTracks, kTRUE);
369	>	ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
370	>	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
371	>	|| fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0) {
372	>	ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
373	>	}
374
375	<	if (fMuonIDType.CompareTo("Tight") == 0)
375	>	if (fMuonIDType.CompareTo("WMuId") == 0)
376	>	fMuIDType = kWMuId;
377	>	else if (fMuonIDType.CompareTo("ZMuId") == 0)
378	>	fMuIDType = kZMuId;
379	>	else if (fMuonIDType.CompareTo("Tight") == 0)
380		fMuIDType = kTight;
381		else if (fMuonIDType.CompareTo("Loose") == 0)
382		fMuIDType = kLoose;
383	+	else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
384	+	fMuIDType = kWWMuIdV1;
385	+	else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
386	+	fMuIDType = kWWMuIdV2;
387	+	else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
388	+	fMuIDType = kWWMuIdV3;
389		else if (fMuonIDType.CompareTo("NoId") == 0)
390		fMuIDType = kNoId;
391		else if (fMuonIDType.CompareTo("Custom") == 0) {
#	Line 163 | Line 405 | void MuonIDMod::SlaveBegin()
405		fMuIsoType = kTrackCaloCombined;
406		else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
407		fMuIsoType = kTrackCaloSliding;
408	+	else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
409	+	fMuIsoType = kTrackCaloSlidingNoCorrection;
410	+	else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
411	+	fMuIsoType = kCombinedRelativeConeAreaCorrected;
412	+	else if (fMuonIsoType.CompareTo("PFIso") == 0)
413	+	fMuIsoType = kPFIso;
414	+	else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
415	+	fMuIsoType = kPFIsoNoL;
416		else if (fMuonIsoType.CompareTo("NoIso") == 0)
417		fMuIsoType = kNoIso;
418		else if (fMuonIsoType.CompareTo("Custom") == 0) {
#	Line 180 | Line 430 | void MuonIDMod::SlaveBegin()
430		fMuClassType = kAll;
431		else if (fMuonClassType.CompareTo("Global") == 0)
432		fMuClassType = kGlobal;
433	+	else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
434	+	fMuClassType = kGlobalTracker;
435		else if (fMuonClassType.CompareTo("Standalone") == 0)
436		fMuClassType = kSta;
437	<	else if (fMuonClassType.CompareTo("TrackerOnly") == 0)
438	<	fMuClassType = kTrackerOnly;
437	>	else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
438	>	fMuClassType = kTrackerMuon;
439	>	else if (fMuonClassType.CompareTo("CaloMuon") == 0)
440	>	fMuClassType = kCaloMuon;
441	>	else if (fMuonClassType.CompareTo("TrackerBased") == 0)
442	>	fMuClassType = kTrackerBased;
443		else {
444		SendError(kAbortAnalysis, "SlaveBegin",
445		"The specified muon class %s is not defined.",

Diff Legend

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