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

Diff Legend

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