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.35 by ceballos, Wed Oct 20 02:44:16 2010 UTC vs.
Revision 1.76 by ceballos, Sat May 5 08:51:06 2012 UTC

#	Line 1 | Line 1
1	<	// $Id$
1	>	// $Id$
2
3		#include "MitPhysics/Mods/interface/MuonIDMod.h"
4		#include "MitCommon/MathTools/interface/MathUtils.h"
5	<	#include "MitAna/DataTree/interface/MuonCol.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
#	Line 13 | Line 14 | ClassImp(mithep::MuonIDMod)
14		//--------------------------------------------------------------------------------------------------
15		MuonIDMod::MuonIDMod(const char *name, const char *title) :
16		BaseMod(name,title),
17	+	fPrintMVADebugInfo(kFALSE),
18		fMuonBranchName(Names::gkMuonBrn),
19		fCleanMuonsName(ModNames::gkCleanMuonsName),
20	+	fNonIsolatedMuonsName("random"),
21	+	fNonIsolatedElectronsName("random"),
22		fVertexName(ModNames::gkGoodVertexesName),
23	<	fMuonIDType("Minimal"),
24	<	fMuonIsoType("TrackCaloSliding"),
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(5.0),
31	>	fCombIsolationCut(0.15),
32	>	fCombRelativeIsolationCut(0.15),
33	>	fPFIsolationCut(-1.0),
34		fMuonPtMin(10),
35		fApplyD0Cut(kTRUE),
36	+	fApplyDZCut(kTRUE),
37		fD0Cut(0.020),
38	+	fDZCut(0.10),
39	+	fWhichVertex(-1),
40		fEtaCut(2.4),
29	–	fReverseIsoCut(kFALSE),
30	–	fReverseD0Cut(kFALSE),
41		fMuIDType(kIdUndef),
42		fMuIsoType(kIsoUndef),
43		fMuClassType(kClassUndef),
44		fMuons(0),
45		fVertices(0),
46	<	fMuonTools(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	>	fTheRhoType(RhoUtilities::DEFAULT)
58		{
59		// Constructor.
60		}
#	Line 43 | Line 64 | void MuonIDMod::Process()
64		{
65		// Process entries of the tree.
66
67	<	LoadEventObject(fMuonBranchName, fMuons);
67	>	if(fMuIsoType != kPFIsoNoL) {
68	>	LoadEventObject(fMuonBranchName, fMuons);
69	>	}
70	>	else {
71	>	fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
72	>	}
73	>	LoadEventObject(fBeamSpotName, fBeamSpot);
74	>	LoadEventObject(fTrackName, fTracks);
75	>	LoadEventObject(fPFCandidatesName, fPFCandidates);
76	>	if(fMuIsoType == kTrackCaloSliding ||
77	>	fMuIsoType == kCombinedRelativeConeAreaCorrected ||
78	>	fMuIsoType == kPFIsoEffectiveAreaCorrected ||
79	>	fMuIsoType == kMVAIso_BDTG_IDIso ||
80	>	fMuIsoType == kIsoRingsV0_BDTG_Iso ||
81	>	fMuIsoType == kIsoDeltaR
82	>	) {
83	>	LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
84	>	}
85	>	if(fMuIsoType == kPFRadialIso || fMuIsoType == kIsoDeltaR){
86	>	// Name is hardcoded, can be changed if someone feels to do it
87	>	fPFNoPileUpCands = GetObjThisEvt<PFCandidateCol>("PFNoPileUp");
88	>	}
89
90		MuonOArr *CleanMuons = new MuonOArr;
91		CleanMuons->SetName(fCleanMuonsName);
92
93	<	for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
93	>	fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
94	>
95	>	for (UInt_t i=0; i<fMuons->GetEntries() && fVertices->GetEntries() > 0 ; ++i) {
96		const Muon *mu = fMuons->At(i);
97
98		Bool_t pass = kFALSE;
#	Line 73 | Line 117 | void MuonIDMod::Process()
117		eta = TMath::Abs(mu->Eta());
118		}
119		break;
120	+	case kGlobalTracker:
121	+	pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
122	+	(mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
123	+	(mu->IsTrackerMuon() &&
124	+	mu->Quality().Quality(MuonQuality::TMLastStationTight));
125	+	if (pass) {
126	+	pt  = mu->TrackerTrk()->Pt();
127	+	eta = TMath::Abs(mu->TrackerTrk()->Eta());
128	+	}
129	+	else {
130	+	pt  = mu->Pt();
131	+	eta = TMath::Abs(mu->Eta());
132	+	}
133	+	break;
134		case kSta:
135		pass = mu->HasStandaloneTrk();
136		if (pass) {
#	Line 115 | Line 173 | void MuonIDMod::Process()
173		if (eta >= fEtaCut)
174		continue;
175
176	+
177	+	//***********************************************************************************************
178	+	//Debug Info For Lepton MVA
179	+	//***********************************************************************************************
180	+	if( fPrintMVADebugInfo &&
181	+	(fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso)
182	+	) {
183	+	cout << "Event: " << GetEventHeader()->RunNum() << " " << GetEventHeader()->LumiSec() << " "
184	+	<< GetEventHeader()->EvtNum() << " : Rho = " << fPileupEnergyDensity->At(0)->Rho()
185	+	<< " : Muon " << i << " "
186	+	<< endl;
187	+	fMuonIDMVA->MVAValue(mu,fVertices->At(0),fMuonTools,fPFCandidates,fPileupEnergyDensity,kTRUE);
188	+	}
189	+	//***********************************************************************************************
190	+
191	+
192		Double_t RChi2 = 0.0;
193		if     (mu->HasGlobalTrk()) {
194		RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
#	Line 128 | Line 202 | void MuonIDMod::Process()
202		idpass = mu->BestTrk() != 0 &&
203		mu->BestTrk()->NHits() > 10 &&
204		RChi2 < 10.0 &&
205	<	mu->NSegments() > 1 &&
205	>	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
206		mu->BestTrk()->NPixelHits() > 0 &&
207		mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
208		break;
209		case kZMuId:
210		idpass = mu->BestTrk() != 0 &&
211		mu->BestTrk()->NHits() > 10 &&
212	<	mu->NSegments() > 1 &&
212	>	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
213		mu->BestTrk()->NPixelHits() > 0 &&
214		mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
215		break;
#	Line 155 | Line 229 | void MuonIDMod::Process()
229		RChi2 < 10.0 &&
230		mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
231		break;
232	<	case kMinimal:
232	>	case kWWMuIdV1:
233		idpass = mu->BestTrk() != 0 &&
234		mu->BestTrk()->NHits() > 10 &&
235	+	mu->BestTrk()->NPixelHits() > 0 &&
236	+	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
237		RChi2 < 10.0 &&
238	+	(mu->NSegments() > 1 || mu->NMatches() > 1) &&
239		mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
240		break;
241	+	case kWWMuIdV2:
242	+	idpass = mu->BestTrk() != 0 &&
243	+	mu->BestTrk()->NHits() > 10 &&
244	+	mu->BestTrk()->NPixelHits() > 0 &&
245	+	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
246	+	break;
247	+	case kWWMuIdV3:
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	+	mu->TrkKink() < 20.0;
253	+	break;
254	+	case kMVAID_BDTG_IDIso:
255	+	{
256	+	Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
257	+	mu->BestTrk()->NHits() > 10 &&
258	+	mu->BestTrk()->NPixelHits() > 0 &&
259	+	mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
260	+	MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
261	+	MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
262	+	mu->TrkKink() < 20.0
263	+	);
264	+	idpass =  passDenominatorM2;
265	+	//only evaluate MVA if muon passes M2 denominator to save time
266	+	if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
267	+	}
268	+	break;
269		case kNoId:
270		idpass = kTRUE;
271		break;
#	Line 171 | Line 276 | void MuonIDMod::Process()
276		if (!idpass)
277		continue;
278
279	+	Double_t Rho = 0.;
280	+	if(fPileupEnergyDensity){
281	+	const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
282	+
283	+	switch (fTheRhoType) {
284	+	case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
285	+	Rho = rho->RhoLowEta();
286	+	break;
287	+	case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
288	+	Rho = rho->Rho();
289	+	break;
290	+	case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
291	+	Rho = rho->RhoRandomLowEta();
292	+	break;
293	+	case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
294	+	Rho = rho->RhoRandom();
295	+	break;
296	+	default:
297	+	Rho = rho->Rho();
298	+	}
299	+
300	+	if ((TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho()))) Rho = 0.;
301	+	}
302	+
303		Bool_t isocut = kFALSE;
304		switch (fMuIsoType) {
305		case kTrackCalo:
#	Line 178 | Line 307 | void MuonIDMod::Process()
307		(mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
308		break;
309		case kTrackCaloCombined:
310	<	isocut = (1.0 * mu->IsoR03SumPt() + 1.0 * mu->IsoR03EmEt() +
311	<	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
310	>	isocut = (1.0 * mu->IsoR03SumPt() +
311	>	1.0 * mu->IsoR03EmEt()  +
312	>	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
313		break;
314		case kTrackCaloSliding:
315		{
316	<	Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
317	<	1.0 * mu->IsoR03EmEt() +
318	<	1.0 * mu->IsoR03HadEt();
319	<	if (totalIso < (mu->Pt()*0.15) )
316	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
317	>	// trick to change the signal region cut
318	>	double theIsoCut = fCombIsolationCut;
319	>	if(theIsoCut < 0.20){
320	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
321	>	else                 theIsoCut = 0.10;
322	>	}
323	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
324	>	}
325	>	break;
326	>	case kTrackCaloSlidingNoCorrection:
327	>	{
328	>	Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
329	>	1.0 * mu->IsoR03EmEt()  +
330	>	1.0 * mu->IsoR03HadEt();
331	>	// trick to change the signal region cut
332	>	double theIsoCut = fCombIsolationCut;
333	>	if(theIsoCut < 0.20){
334	>	if(mu->Pt() >  20.0) theIsoCut = 0.15;
335	>	else                 theIsoCut = 0.10;
336	>	}
337	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
338	>	}
339	>	break;
340	>	case kCombinedRelativeConeAreaCorrected:
341	>	{
342	>	//const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0); // Fabian: made Rho customable
343	>	Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
344	>	double theIsoCut = fCombRelativeIsolationCut;
345	>	if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
346	>	}
347	>	break;
348	>	case kCombinedRelativeEffectiveAreaCorrected:
349	>	{
350	>	Double_t tmpRho = Rho;   // Fabian: made the Rho type customable.
351	>	//if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))
352	>	//tmpRho = fPileupEnergyDensity->At(0)->Rho();
353	>
354	>	isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
355	>	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
356	>	-  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
357	>	) < (mu->Pt()* 0.40);
358	>	}
359	>	break;
360	>	case kPFIso:
361	>	{
362	>	Double_t pfIsoCutValue = 9999;
363	>	if(fPFIsolationCut > 0){
364	>	pfIsoCutValue = fPFIsolationCut;
365	>	} else {
366	>	if (mu->AbsEta() < 1.479) {
367	>	if (mu->Pt() > 20) {
368	>	pfIsoCutValue = 0.13;
369	>	} else {
370	>	pfIsoCutValue = 0.06;
371	>	}
372	>	} else {
373	>	if (mu->Pt() > 20) {
374	>	pfIsoCutValue = 0.09;
375	>	} else {
376	>	pfIsoCutValue = 0.05;
377	>	}
378	>	}
379	>	}
380	>	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
381	>	if (totalIso < (mu->Pt()*pfIsoCutValue) )
382	>	isocut = kTRUE;
383	>	}
384	>	break;
385	>	case kPFRadialIso:
386	>	{
387	>	Double_t pfIsoCutValue = 9999;
388	>	if(fPFIsolationCut > 0){
389	>	pfIsoCutValue = fPFIsolationCut;
390	>	} else {
391	>	if (mu->Pt() > 20) {
392	>	pfIsoCutValue = 0.10;
393	>	} else {
394	>	pfIsoCutValue = 0.05;
395	>	}
396	>	}
397	>	Double_t totalIso =  IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
398	>	if (totalIso < (mu->Pt()*pfIsoCutValue) )
399	>	isocut = kTRUE;
400	>	}
401	>	break;
402	>	case kPFIsoEffectiveAreaCorrected:
403	>	{
404	>	Double_t pfIsoCutValue = 9999;
405	>	if(fPFIsolationCut > 0){
406	>	pfIsoCutValue = fPFIsolationCut;
407	>	} else {
408	>	pfIsoCutValue = fPFIsolationCut; //leave it like this for now
409	>	}
410	>	Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
411	>	- Rho * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
412	>	//- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());  // Fabian: made Rho-type customable
413	>	isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
414	>	break;
415	>	}
416	>	case kPFIsoNoL:
417	>	{
418	>	fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
419	>	fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
420	>
421	>	Double_t pfIsoCutValue = 9999;
422	>	if(fPFIsolationCut > 0){
423	>	pfIsoCutValue = fPFIsolationCut;
424	>	} else {
425	>	if (mu->AbsEta() < 1.479) {
426	>	if (mu->Pt() > 20) {
427	>	pfIsoCutValue = 0.13;
428	>	} else {
429	>	pfIsoCutValue = 0.06;
430	>	}
431	>	} else {
432	>	if (mu->Pt() > 20) {
433	>	pfIsoCutValue = 0.09;
434	>	} else {
435	>	pfIsoCutValue = 0.05;
436	>	}
437	>	}
438	>	}
439	>	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
440	>	if (totalIso < (mu->Pt()*pfIsoCutValue) )
441		isocut = kTRUE;
191	–
192	–	if     (fReverseIsoCut == kTRUE &&
193	–	isocut == kFALSE && totalIso < 10)
194	–	isocut = kTRUE;
195	–	else if(fReverseIsoCut == kTRUE)
196	–	isocut = kFALSE;
442		}
443		break;
444	+	case kMVAIso_BDTG_IDIso:
445	+	{
446	+
447	+	// **************************************************************************
448	+	// Don't use effective area correction denominator. Instead use the old one.
449	+	// **************************************************************************
450	+
451	+	//         Double_t tmpRho = 0;
452	+	//         if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || isinf(fPileupEnergyDensity->At(0)->Rho())))
453	+	//           tmpRho = fPileupEnergyDensity->At(0)->Rho();
454	+
455	+	//         isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
456	+	//                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
457	+	//                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
458	+	//           ) < (mu->Pt()* 0.40);
459	+
460	+	Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
461	+	isocut = (totalIso < (mu->Pt()*0.4));
462	+
463	+	}
464	+	break;
465	+	case kIsoRingsV0_BDTG_Iso:
466	+	{
467	+
468	+	isocut = PassMuonIsoRingsV0_BDTG_Iso(mu, fVertices->At(0), fPileupEnergyDensity);
469	+
470	+	}
471	+	break;
472	+	case kIsoDeltaR:
473	+	{
474	+
475	+	isocut = PassMuonIsoDeltaR(mu, fVertices->At(0), fPileupEnergyDensity);
476	+
477	+	}
478	+	break;
479		case kNoIso:
480		isocut = kTRUE;
481		break;
#	Line 207 | Line 487 | void MuonIDMod::Process()
487		if (isocut == kFALSE)
488		continue;
489
490	+	// apply d0 cut
491		if (fApplyD0Cut) {
492	<	fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
493	<	Bool_t passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fReverseD0Cut);
492	>	Bool_t passD0cut = kTRUE;
493	>	if(fD0Cut < 0.05) { // trick to change the signal region cut
494	>	if      (mu->Pt() >  20.0) fD0Cut = 0.02;
495	>	else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
496	>	}
497	>	if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
498	>	else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
499		if (!passD0cut)
500		continue;
501		}
502
503	+	// apply dz cut
504	+	if (fApplyDZCut) {
505	+	Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
506	+	if (!passDZcut)
507	+	continue;
508	+	}
509	+
510		// add good muon
511		CleanMuons->Add(mu);
512		}
#	Line 231 | Line 524 | void MuonIDMod::SlaveBegin()
524		// Run startup code on the computer (slave) doing the actual analysis. Here,
525		// we just request the muon collection branch.
526
527	<	ReqEventObject(fMuonBranchName, fMuons, kTRUE);
527	>	// In this case we cannot have a branch
528	>	if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
529	>	ReqEventObject(fMuonBranchName, fMuons, kTRUE);
530	>	}
531	>	ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
532	>	ReqEventObject(fTrackName, fTracks, kTRUE);
533	>	ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
534	>	if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
535	>	|| fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
536	>	|| fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
537	>	|| fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
538	>	|| fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
539	>	|| fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0
540	>	|| fMuonIsoType.CompareTo("IsoDeltaR") == 0
541	>	) {
542	>	ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
543	>	}
544
236	–	fMuonTools = new MuonTools;
545
546		if (fMuonIDType.CompareTo("WMuId") == 0)
547		fMuIDType = kWMuId;
#	Line 243 | Line 551 | void MuonIDMod::SlaveBegin()
551		fMuIDType = kTight;
552		else if (fMuonIDType.CompareTo("Loose") == 0)
553		fMuIDType = kLoose;
554	<	else if (fMuonIDType.CompareTo("Minimal") == 0)
555	<	fMuIDType = kMinimal;
554	>	else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
555	>	fMuIDType = kWWMuIdV1;
556	>	else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
557	>	fMuIDType = kWWMuIdV2;
558	>	else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
559	>	fMuIDType = kWWMuIdV3;
560		else if (fMuonIDType.CompareTo("NoId") == 0)
561		fMuIDType = kNoId;
562		else if (fMuonIDType.CompareTo("Custom") == 0) {
563		fMuIDType = kCustomId;
564		SendError(kWarning, "SlaveBegin",
565		"Custom muon identification is not yet implemented.");
566	+	} else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
567	+	fMuIDType = kMVAID_BDTG_IDIso;
568		} else {
569		SendError(kAbortAnalysis, "SlaveBegin",
570		"The specified muon identification %s is not defined.",
#	Line 264 | Line 578 | void MuonIDMod::SlaveBegin()
578		fMuIsoType = kTrackCaloCombined;
579		else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
580		fMuIsoType = kTrackCaloSliding;
581	+	else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
582	+	fMuIsoType = kTrackCaloSlidingNoCorrection;
583	+	else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
584	+	fMuIsoType = kCombinedRelativeConeAreaCorrected;
585	+	else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
586	+	fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
587	+	else if (fMuonIsoType.CompareTo("PFIso") == 0)
588	+	fMuIsoType = kPFIso;
589	+	else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
590	+	fMuIsoType = kPFRadialIso;
591	+	else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
592	+	fMuIsoType = kPFIsoEffectiveAreaCorrected;
593	+	else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
594	+	fMuIsoType = kPFIsoNoL;
595		else if (fMuonIsoType.CompareTo("NoIso") == 0)
596		fMuIsoType = kNoIso;
597		else if (fMuonIsoType.CompareTo("Custom") == 0) {
598		fMuIsoType = kCustomIso;
599		SendError(kWarning, "SlaveBegin",
600		"Custom muon isolation is not yet implemented.");
601	+	} else if (fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0) {
602	+	fMuIsoType = kMVAIso_BDTG_IDIso;
603	+	} else if (fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0) {
604	+	fMuIsoType = kIsoRingsV0_BDTG_Iso;
605	+	} else if (fMuonIsoType.CompareTo("IsoDeltaR") == 0) {
606	+	fMuIsoType = kIsoDeltaR;
607		} else {
608		SendError(kAbortAnalysis, "SlaveBegin",
609		"The specified muon isolation %s is not defined.",
#	Line 281 | Line 615 | void MuonIDMod::SlaveBegin()
615		fMuClassType = kAll;
616		else if (fMuonClassType.CompareTo("Global") == 0)
617		fMuClassType = kGlobal;
618	+	else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
619	+	fMuClassType = kGlobalTracker;
620		else if (fMuonClassType.CompareTo("Standalone") == 0)
621		fMuClassType = kSta;
622		else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
#	Line 295 | Line 631 | void MuonIDMod::SlaveBegin()
631		fMuonClassType.Data());
632		return;
633		}
634	+
635	+
636	+	//If we use MVA ID, need to load MVA weights
637	+	if     (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
638	+	fMuonTools = new MuonTools();
639	+	fMuonIDMVA = new MuonIDMVA();
640	+	fMuonIDMVA->Initialize("BDTG method",
641	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin0_IDIsoCombined_BDTG.weights.xml"))),
642	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin0_IDIsoCombined_BDTG.weights.xml"))),
643	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin1_IDIsoCombined_BDTG.weights.xml"))),
644	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin1_IDIsoCombined_BDTG.weights.xml"))),
645	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin2_IDIsoCombined_BDTG.weights.xml"))),
646	+	string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin2_IDIsoCombined_BDTG.weights.xml"))),
647	+	MuonIDMVA::kIDIsoCombinedDetIso);
648	+	}
649	+	else if(fMuIsoType == kIsoRingsV0_BDTG_Iso) {
650	+	std::vector<std::string> muonidiso_weightfiles;
651	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml"))));
652	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml"))));
653	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml"))));
654	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml"))));
655	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml"))));
656	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml"))));
657	+	fMuonTools = new MuonTools();
658	+	fMuonIDMVA = new MuonIDMVA();
659	+	fMuonIDMVA->Initialize("MuonIso_BDTG_IsoRings",
660	+	MuonIDMVA::kIsoRingsV0,
661	+	kTRUE,
662	+	muonidiso_weightfiles);
663	+	}
664	+	else if(fMuIsoType == kIsoDeltaR) {
665	+	std::vector<std::string> muonidiso_weightfiles;
666	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LB_BDT.weights.xml"))));
667	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LE_BDT.weights.xml"))));
668	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HB_BDT.weights.xml"))));
669	+	muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HE_BDT.weights.xml"))));
670	+	fMuonTools = new MuonTools();
671	+	fMuonIDMVA = new MuonIDMVA();
672	+	fMuonIDMVA->Initialize("muonHZZ2012IsoDRMVA",
673	+	MuonIDMVA::kIsoDeltaR,
674	+	kTRUE,
675	+	muonidiso_weightfiles);
676	+	}
677	+
678	+	}
679	+
680	+
681	+	//--------------------------------------------------------------------------------------------------
682	+	Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
683	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
684	+	{
685	+
686	+	const Track *muTrk=0;
687	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
688	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
689	+
690	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
691	+
692	+	Int_t subdet = 0;
693	+	if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
694	+	else subdet = 1;
695	+	Int_t ptBin = 0;
696	+	if (muTrk->Pt() > 14.5) ptBin = 1;
697	+	if (muTrk->Pt() > 20.0) ptBin = 2;
698	+
699	+	Int_t MVABin = -1;
700	+	if      (subdet == 0 && ptBin == 0) MVABin = 0;
701	+	else if (subdet == 1 && ptBin == 0) MVABin = 1;
702	+	else if (subdet == 0 && ptBin == 1) MVABin = 2;
703	+	else if (subdet == 1 && ptBin == 1) MVABin = 3;
704	+	else if (subdet == 0 && ptBin == 2) MVABin = 4;
705	+	else if (subdet == 1 && ptBin == 2) MVABin = 5;
706	+
707	+	Double_t MVACut = -999;
708	+	if      (MVABin == 0) MVACut = -0.5618;
709	+	else if (MVABin == 1) MVACut = -0.3002;
710	+	else if (MVABin == 2) MVACut = -0.4642;
711	+	else if (MVABin == 3) MVACut = -0.2478;
712	+	else if (MVABin == 4) MVACut =  0.1706;
713	+	else if (MVABin == 5) MVACut =  0.8146;
714	+
715	+	if (MVAValue > MVACut) return kTRUE;
716	+	return kFALSE;
717	+	}
718	+
719	+	//--------------------------------------------------------------------------------------------------
720	+	Bool_t MuonIDMod::PassMuonIsoRingsV0_BDTG_Iso(const Muon *mu, const Vertex *vertex,
721	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
722	+	{
723	+
724	+	const Track *muTrk=0;
725	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
726	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
727	+
728	+	ElectronOArr *tempElectrons = new  ElectronOArr;
729	+	MuonOArr     *tempMuons     = new  MuonOArr;
730	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,
731	+	PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,kFALSE);
732	+	delete tempElectrons;
733	+	delete tempMuons;
734	+
735	+	Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
736	+
737	+	Double_t MVACut = -999;
738	+	if      (MVABin == 0) MVACut = -0.593;
739	+	else if (MVABin == 1) MVACut =  0.337;
740	+	else if (MVABin == 2) MVACut = -0.767;
741	+	else if (MVABin == 3) MVACut =  0.410;
742	+	else if (MVABin == 4) MVACut = -0.989;
743	+	else if (MVABin == 5) MVACut = -0.995;
744	+
745	+	if (MVAValue > MVACut) return kTRUE;
746	+	return kFALSE;
747	+	}
748	+
749	+	//--------------------------------------------------------------------------------------------------
750	+	Bool_t MuonIDMod::PassMuonIsoDeltaR(const Muon *mu, const Vertex *vertex,
751	+	const PileupEnergyDensityCol *PileupEnergyDensity) const
752	+	{
753	+
754	+	const Track *muTrk=0;
755	+	if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
756	+	else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
757	+
758	+	ElectronOArr *tempElectrons = new  ElectronOArr;
759	+	MuonOArr     *tempMuons     = new  MuonOArr;
760	+	Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFNoPileUpCands,
761	+	PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,kTRUE);
762	+	delete tempElectrons;
763	+	delete tempMuons;
764	+
765	+	Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
766	+
767	+	Double_t MVACut = -999;
768	+	if      (MVABin == 0) MVACut =  0.000;
769	+	else if (MVABin == 1) MVACut =  0.000;
770	+	else if (MVABin == 2) MVACut =  0.000;
771	+	else if (MVABin == 3) MVACut =  0.000;
772	+
773	+	if (MVAValue > MVACut) return kTRUE;
774	+	return kFALSE;
775	+	}
776	+
777	+	//--------------------------------------------------------------------------------------------------
778	+	void MuonIDMod::Terminate()
779	+	{
780	+	// Run finishing code on the computer (slave) that did the analysis
781	+	delete fMuonIDMVA;
782	+
783	+	delete fMuonTools;
784		}

Diff Legend

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