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.1 by loizides, Wed Oct 15 06:05:00 2008 UTC vs.
Revision 1.41 by ceballos, Fri Mar 11 15:13:09 2011 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"),
17	>	fMuonBranchName(Names::gkMuonBrn),
18	>	fCleanMuonsName(ModNames::gkCleanMuonsName),
19	>	fNonIsolatedMuonsName("random"),
20	>	fNonIsolatedElectronsName("random"),
21	>	fVertexName(ModNames::gkGoodVertexesName),
22	>	fTrackName(Names::gkTrackBrn),
23	>	fPFCandidatesName(Names::gkPFCandidatesBrn),
24	>	fMuonIDType("WWMuId"),
25	>	fMuonIsoType("TrackCaloSliding"),
26	>	fMuonClassType("Global"),
27	>	fTrackIsolationCut(3.0),
28	>	fCaloIsolationCut(3.0),
29	>	fCombIsolationCut(0.15),
30	>	fMuonPtMin(10),
31	>	fApplyD0Cut(kTRUE),
32	>	fD0Cut(0.020),
33	>	fEtaCut(2.4),
34	>	fReverseIsoCut(kFALSE),
35	>	fReverseD0Cut(kFALSE),
36	>	fMuIDType(kIdUndef),
37	>	fMuIsoType(kIsoUndef),
38	>	fMuClassType(kClassUndef),
39		fMuons(0),
40	<	fTrackIsolationCut(5.0),
41	<	fCaloIsolationCut(5.0),
42	<	fNEventsProcessed(0)
40	>	fVertices(0),
41	>	fTracks(0),
42	>	fPFCandidates(0),
43	>	fNonIsolatedMuons(0),
44	>	fNonIsolatedElectrons(0),
45	>	fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
46	>	fPileupEnergyDensity(0)
47		{
48		// Constructor.
49		}
50
51		//--------------------------------------------------------------------------------------------------
30	–	void MuonIDMod::Begin()
31	–	{
32	–	// Run startup code on the client machine. For this module, we dont do
33	–	// anything here.
34	–	}
35	–
36	–	//--------------------------------------------------------------------------------------------------
52		void MuonIDMod::Process()
53		{
54		// Process entries of the tree.
55
56	<	fNEventsProcessed++;
57	<
58	<	if (fNEventsProcessed % 1000 == 0 || fPrintDebug) {
59	<	time_t systime;
60	<	systime = time(NULL);
61	<
62	<	cerr << endl << "MuonIDMod : Process Event " << fNEventsProcessed << "  Time: " << ctime(&systime) << endl;
63	<	}
64	<
65	<	//Get Muons
66	<	LoadBranch(fMuonName);
67	<	ObjArray<Muon> *CleanMuons = new ObjArray<Muon>;
56	>	if(fMuIsoType != kPFIsoNoL) {
57	>	LoadEventObject(fMuonBranchName, fMuons);
58	>	}
59	>	else {
60	>	fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
61	>	}
62	>	LoadEventObject(fTrackName, fTracks);
63	>	LoadEventObject(fPFCandidatesName, fPFCandidates);
64	>	if(fMuIsoType == kTrackCaloSliding) {
65	>	LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
66	>	}
67	>	MuonOArr *CleanMuons = new MuonOArr;
68	>	CleanMuons->SetName(fCleanMuonsName);
69	>
70	>	fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
71	>
72		for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
73	<	Muon *mu = fMuons->At(i);
55	<
56	<	Double_t MuonClass = -1;
57	<	if (mu->GlobalTrk())
58	<	MuonClass = 0;
59	<	else if (mu->StandaloneTrk())
60	<	MuonClass = 1;
61	<	else if (mu->TrackerTrk())
62	<	MuonClass = 2;
63	<
64	<	//These cuts are from the 1.6.X analysis. I'm waiting for Phil to finalize his Muon ID class
65	<	const int nCuts = 4;
66	<	double cutValue[nCuts] = {0.1, 3.0, 3.0, 1.5 };
67	<	bool passCut[nCuts] = {false, false, false, false};
68	<	double muonD0 = fabs(mu->BestTrk()->D0());
69	<	if(muonD0 < cutValue[0] &&  MuonClass == 0 )
70	<	passCut[0] = true;
71	<	if(mu->IsoR03SumPt() < cutValue[1]) passCut[1] = true;
72	<	if(mu->IsoR03EmEt() +
73	<	mu->IsoR03HadEt() < cutValue[2]) passCut[2] = true;
74	<	if(mu->Pt() > 10)
75	<	passCut[3] = true;
76	<
77	<	// Final decision
78	<	bool allCuts = true;
79	<	for(int c=0; c<nCuts; c++) {
80	<	allCuts = allCuts & passCut[c];
81	<	}
82	<
83	<	if ( allCuts
84	<	&& abs(mu->Eta()) < 2.5
85	<	) {
86	<	CleanMuons->Add(mu);
87	<	}
73	>	const Muon *mu = fMuons->At(i);
74
75	+	Bool_t pass = kFALSE;
76	+	Double_t pt = 0;  // make sure pt is taken from the correct track!
77	+	Double_t eta = 0; // make sure eta is taken from the correct track!
78	+	switch (fMuClassType) {
79	+	case kAll:
80	+	pass = kTRUE;
81	+	if (mu->HasTrk()) {
82	+	pt  = mu->Pt();
83	+	eta = TMath::Abs(mu->Eta());
84	+	}
85	+	break;
86	+	case kGlobal:
87	+	pass = mu->HasGlobalTrk() && mu->IsTrackerMuon();
88	+	if (pass && mu->TrackerTrk()) {
89	+	pt  = mu->TrackerTrk()->Pt();
90	+	eta = TMath::Abs(mu->TrackerTrk()->Eta());
91	+	}
92	+	else {
93	+	pt  = mu->Pt();
94	+	eta = TMath::Abs(mu->Eta());
95	+	}
96	+	break;
97	+	case kSta:
98	+	pass = mu->HasStandaloneTrk();
99	+	if (pass) {
100	+	pt  = mu->StandaloneTrk()->Pt();
101	+	eta = TMath::Abs(mu->StandaloneTrk()->Eta());
102	+	}
103	+	break;
104	+	case kTrackerMuon:
105	+	pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
106	+	mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
107	+	if (pass) {
108	+	pt  = mu->TrackerTrk()->Pt();
109	+	eta = TMath::Abs(mu->TrackerTrk()->Eta());
110	+	}
111	+	break;
112	+	case kCaloMuon:
113	+	pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
114	+	if (pass) {
115	+	pt  = mu->TrackerTrk()->Pt();
116	+	eta = TMath::Abs(mu->TrackerTrk()->Eta());
117	+	}
118	+	break;
119	+	case kTrackerBased:
120	+	pass = mu->HasTrackerTrk();
121	+	if (pass) {
122	+	pt  = mu->TrackerTrk()->Pt();
123	+	eta = TMath::Abs(mu->TrackerTrk()->Eta());
124	+	}
125	+	break;
126	+	default:
127	+	break;
128	+	}
129	+
130	+	if (!pass)
131	+	continue;
132	+
133	+	if (pt <= fMuonPtMin)
134	+	continue;
135	+
136	+	if (eta >= fEtaCut)
137	+	continue;
138	+
139	+	Double_t RChi2 = 0.0;
140	+	if     (mu->HasGlobalTrk()) {
141	+	RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
142	+	}
143	+	else if(mu->BestTrk() != 0){
144	+	RChi2 = mu->BestTrk()->Chi2()/mu->BestTrk()->Ndof();
145	+	}
146	+	Bool_t idpass = kFALSE;
147	+	switch (fMuIDType) {
148	+	case kWMuId:
149	+	idpass = mu->BestTrk() != 0 &&
150	+	mu->BestTrk()->NHits() > 10 &&
151	+	RChi2 < 10.0 &&
152	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
153	+	mu->BestTrk()->NPixelHits() > 0 &&
154	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
155	+	break;
156	+	case kZMuId:
157	+	idpass = mu->BestTrk() != 0 &&
158	+	mu->BestTrk()->NHits() > 10 &&
159	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
160	+	mu->BestTrk()->NPixelHits() > 0 &&
161	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
162	+	break;
163	+	case kLoose:
164	+	idpass = mu->BestTrk() != 0 &&
165	+	mu->Quality().Quality(MuonQuality::TMOneStationLoose) &&
166	+	mu->Quality().Quality(MuonQuality::TM2DCompatibilityLoose) &&
167	+	mu->BestTrk()->NHits() > 10 &&
168	+	RChi2 < 10.0 &&
169	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
170	+	break;
171	+	case kTight:
172	+	idpass = mu->BestTrk() !=  0 &&
173	+	mu->Quality().Quality(MuonQuality::TMOneStationTight) &&
174	+	mu->Quality().Quality(MuonQuality::TM2DCompatibilityTight) &&
175	+	mu->BestTrk()->NHits() > 10 &&
176	+	RChi2 < 10.0 &&
177	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
178	+	break;
179	+	case kWWMuId:
180	+	idpass = mu->BestTrk() != 0 &&
181	+	mu->BestTrk()->NHits() > 10 &&
182	+	RChi2 < 10.0 &&
183	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
184	+	mu->BestTrk()->NPixelHits() > 0 &&
185	+	mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight) &&
186	+	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
187	+	break;
188	+	case kNoId:
189	+	idpass = kTRUE;
190	+	break;
191	+	default:
192	+	break;
193	+	}
194	+
195	+	if (!idpass)
196	+	continue;
197	+
198	+	Bool_t isocut = kFALSE;
199	+	switch (fMuIsoType) {
200	+	case kTrackCalo:
201	+	isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
202	+	(mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
203	+	break;
204	+	case kTrackCaloCombined:
205	+	isocut = (1.0 * mu->IsoR03SumPt() +
206	+	1.0 * mu->IsoR03EmEt()  +
207	+	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
208	+	break;
209	+	case kTrackCaloSliding:
210	+	{
211	+	//Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
212	+	//if(beta == 0) beta = 1.0;
213	+	const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
214	+	Double_t totalIso =  mu->IsoR03SumPt() + TMath::Max(mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3, 0.0);
215	+	if (totalIso < (mu->Pt()*fCombIsolationCut) )
216	+	isocut = kTRUE;
217	+
218	+	if     (fReverseIsoCut == kTRUE &&
219	+	isocut == kFALSE && totalIso < 10)
220	+	isocut = kTRUE;
221	+	else if(fReverseIsoCut == kTRUE)
222	+	isocut = kFALSE;
223	+	}
224	+	break;
225	+	case kTrackCaloSlidingNoCorrection:
226	+	{
227	+	Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
228	+	1.0 * mu->IsoR03EmEt()  +
229	+	1.0 * mu->IsoR03HadEt();
230	+	if (totalIso < (mu->Pt()*fCombIsolationCut) )
231	+	isocut = kTRUE;
232	+	}
233	+	break;
234	+	case kPFIso:
235	+	{
236	+	Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
237	+	if(beta == 0) beta = 1.0;
238	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 0, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
239	+	if (totalIso < (mu->Pt()*fCombIsolationCut) )
240	+	isocut = kTRUE;
241	+	}
242	+	break;
243	+	case kPFIsoNoL:
244	+	{
245	+	fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
246	+	fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
247	+
248	+	Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
249	+	if(beta == 0) beta = 1.0;
250	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 3, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
251	+	if (totalIso < (mu->Pt()*fCombIsolationCut) )
252	+	isocut = kTRUE;
253	+	}
254	+	break;
255	+	case kNoIso:
256	+	isocut = kTRUE;
257	+	break;
258	+	case kCustomIso:
259	+	default:
260	+	break;
261	+	}
262	+
263	+	if (isocut == kFALSE)
264	+	continue;
265	+
266	+	if (fApplyD0Cut) {
267	+	Bool_t passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut);
268	+	if (!passD0cut)
269	+	continue;
270	+	}
271
272	+	// add good muon
273	+	CleanMuons->Add(mu);
274		}
275
276	<	//Final Summary Debug Output
277	<	if ( fPrintDebug ) {
94	<	cerr << "Event Dump: " << fNEventsProcessed << endl;
95	<	cerr << "Muons" << endl;
96	<	for (UInt_t i = 0; i < CleanMuons->GetEntries(); i++) {
97	<	cerr << i << " " << CleanMuons->At(i)->Pt() << " " << CleanMuons->At(i)->Eta()
98	<	<< " " << CleanMuons->At(i)->Phi() << endl;
99	<	}
100	<	}
101	<
102	<	//Save Objects for Other Modules to use
103	<	AddObjThisEvt(CleanMuons, fCleanMuonsName.Data());
104	<	}
276	>	// sort according to pt
277	>	CleanMuons->Sort();
278
279	+	// add objects for other modules to use
280	+	AddObjThisEvt(CleanMuons);
281	+	}
282
283		//--------------------------------------------------------------------------------------------------
284		void MuonIDMod::SlaveBegin()
285		{
286		// Run startup code on the computer (slave) doing the actual analysis. Here,
287	<	// we typically initialize histograms and other analysis objects and request
112	<	// branches. For this module, we request a branch of the MitTree.
287	>	// we just request the muon collection branch.
288
289	<	ReqBranch(fMuonName,              fMuons);
290	<	}
289	>	// In this case we cannot have a branch
290	>	if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
291	>	ReqEventObject(fMuonBranchName, fMuons, kTRUE);
292	>	}
293	>	ReqEventObject(fTrackName, fTracks, kTRUE);
294	>	ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
295	>	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0) {
296	>	ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
297	>	}
298
299	<	//--------------------------------------------------------------------------------------------------
300	<	void MuonIDMod::SlaveTerminate()
301	<	{
302	<	// Run finishing code on the computer (slave) that did the analysis. For this
303	<	// module, we dont do anything here.
299	>	if (fMuonIDType.CompareTo("WMuId") == 0)
300	>	fMuIDType = kWMuId;
301	>	else if (fMuonIDType.CompareTo("ZMuId") == 0)
302	>	fMuIDType = kZMuId;
303	>	else if (fMuonIDType.CompareTo("Tight") == 0)
304	>	fMuIDType = kTight;
305	>	else if (fMuonIDType.CompareTo("Loose") == 0)
306	>	fMuIDType = kLoose;
307	>	else if (fMuonIDType.CompareTo("WWMuId") == 0)
308	>	fMuIDType = kWWMuId;
309	>	else if (fMuonIDType.CompareTo("NoId") == 0)
310	>	fMuIDType = kNoId;
311	>	else if (fMuonIDType.CompareTo("Custom") == 0) {
312	>	fMuIDType = kCustomId;
313	>	SendError(kWarning, "SlaveBegin",
314	>	"Custom muon identification is not yet implemented.");
315	>	} else {
316	>	SendError(kAbortAnalysis, "SlaveBegin",
317	>	"The specified muon identification %s is not defined.",
318	>	fMuonIDType.Data());
319	>	return;
320	>	}
321
322	<	}
322	>	if (fMuonIsoType.CompareTo("TrackCalo") == 0)
323	>	fMuIsoType = kTrackCalo;
324	>	else if (fMuonIsoType.CompareTo("TrackCaloCombined") == 0)
325	>	fMuIsoType = kTrackCaloCombined;
326	>	else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
327	>	fMuIsoType = kTrackCaloSliding;
328	>	else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
329	>	fMuIsoType = kTrackCaloSlidingNoCorrection;
330	>	else if (fMuonIsoType.CompareTo("PFIso") == 0)
331	>	fMuIsoType = kPFIso;
332	>	else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
333	>	fMuIsoType = kPFIsoNoL;
334	>	else if (fMuonIsoType.CompareTo("NoIso") == 0)
335	>	fMuIsoType = kNoIso;
336	>	else if (fMuonIsoType.CompareTo("Custom") == 0) {
337	>	fMuIsoType = kCustomIso;
338	>	SendError(kWarning, "SlaveBegin",
339	>	"Custom muon isolation is not yet implemented.");
340	>	} else {
341	>	SendError(kAbortAnalysis, "SlaveBegin",
342	>	"The specified muon isolation %s is not defined.",
343	>	fMuonIsoType.Data());
344	>	return;
345	>	}
346
347	<	//--------------------------------------------------------------------------------------------------
348	<	void MuonIDMod::Terminate()
349	<	{
350	<	// Run finishing code on the client computer. For this module, we dont do
351	<	// anything here.
347	>	if (fMuonClassType.CompareTo("All") == 0)
348	>	fMuClassType = kAll;
349	>	else if (fMuonClassType.CompareTo("Global") == 0)
350	>	fMuClassType = kGlobal;
351	>	else if (fMuonClassType.CompareTo("Standalone") == 0)
352	>	fMuClassType = kSta;
353	>	else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
354	>	fMuClassType = kTrackerMuon;
355	>	else if (fMuonClassType.CompareTo("CaloMuon") == 0)
356	>	fMuClassType = kCaloMuon;
357	>	else if (fMuonClassType.CompareTo("TrackerBased") == 0)
358	>	fMuClassType = kTrackerBased;
359	>	else {
360	>	SendError(kAbortAnalysis, "SlaveBegin",
361	>	"The specified muon class %s is not defined.",
362	>	fMuonClassType.Data());
363	>	return;
364	>	}
365		}

Diff Legend

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