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

Diff Legend

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