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

Diff Legend

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