Mods/src/MuonIDMod.cc

// $Id: MuonIDMod.cc,v 1.23 2009/06/15 15:00:21 loizides Exp $

#include "MitPhysics/Mods/interface/MuonIDMod.h"
#include "MitCommon/MathTools/interface/MathUtils.h"
#include "MitAna/DataTree/interface/MuonCol.h"
#include "MitAna/DataTree/interface/VertexCol.h"
#include "MitPhysics/Init/interface/ModNames.h"

using namespace mithep;

ClassImp(mithep::MuonIDMod)

//--------------------------------------------------------------------------------------------------
  MuonIDMod::MuonIDMod(const char *name, const char *title) : 
  BaseMod(name,title),
  fMuonBranchName(Names::gkMuonBrn),
  fCleanMuonsName(ModNames::gkCleanMuonsName),  
  fVertexName("PrimaryVertexesBeamSpot"),
  fMuonIDType("Loose"),
  fMuonIsoType("TrackCaloSliding"),  
  fMuonClassType("Global"),  
  fTrackIsolationCut(3.0),
  fCaloIsolationCut(3.0),
  fCombIsolationCut(5.0),
  fMuonPtMin(10),
  fApplyD0Cut(kTRUE),
  fD0Cut(0.025),
  fReverseIsoCut(kFALSE),
  fReverseD0Cut(kFALSE),
  fMuIDType(kIdUndef),
  fMuIsoType(kIsoUndef),
  fMuClassType(kClassUndef),
  fMuons(0),
  fVertices(0),
  fMuonTools(0)
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
void MuonIDMod::Process()
{
  // Process entries of the tree. 

  LoadEventObject(fMuonBranchName, fMuons);
  LoadEventObject(fVertexName,     fVertices);

  MuonOArr *CleanMuons = new MuonOArr;
  CleanMuons->SetName(fCleanMuonsName);

  for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
    const Muon *mu = fMuons->At(i);

    Bool_t pass = kFALSE;
    Double_t pt = 0; // make sure pt is taken from the correct track!
    switch (fMuClassType) {
      case kAll:
        pass = kTRUE;
        if (mu->HasTrk())
          pt = mu->Pt();
        break;
      case kGlobal:
        pass = mu->HasGlobalTrk();
        if (pass) 
          pt = mu->GlobalTrk()->Pt();
        break;
      case kSta:
        pass = mu->HasStandaloneTrk();
        if (pass) 
          pt = mu->StandaloneTrk()->Pt();
        break;
      case kTrackerMuon:
        pass = mu->HasTrackerTrk() && mu->IsTrackerMuon();
        if (pass) 
          pt = mu->TrackerTrk()->Pt();
        break;
      case kCaloMuon:
        pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
        if (pass) 
          pt = mu->TrackerTrk()->Pt();
        break;
      case kTrackerBased:
        pass = mu->HasTrackerTrk();
        if (pass) 
          pt = mu->TrackerTrk()->Pt();
        break;
      default:
        break;
    }

    if (!pass)
      continue;

    if (pt <= fMuonPtMin) 
      continue;

    Bool_t idpass = kFALSE;
    switch (fMuIDType) {
      case kLoose:
        idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationLoose) &&
                 fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityLoose);
        break;
      case kTight:
        idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationTight) &&
                 fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityTight);
        break;
      case kNoId:
        idpass = kTRUE;
        break;
      default:
        break;
    }

    if (!idpass)
      continue;

    Bool_t isocut = kFALSE;
    switch (fMuIsoType) {
      case kTrackCalo:
        isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
          (mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
        break;
      case kTrackCaloCombined:
        isocut = (1.0 * mu->IsoR03SumPt() + 1.0 * mu->IsoR03EmEt() + 
                   1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
        break;
      case kTrackCaloSliding:
        { 
          Double_t totalIso = 1.0 * mu->IsoR03SumPt() + 
                              1.0 * mu->IsoR03EmEt() + 
                              1.0 * mu->IsoR03HadEt();
          if ((totalIso < (pt-10.0)*5.0/15.0 && pt <= 25) ||
              (totalIso < 5.0 && mu->Pt() > 25) ||
               totalIso <= 0) 
            isocut = kTRUE;

          if     (fReverseIsoCut == kTRUE &&
                  isocut == kFALSE && totalIso < 10)
            isocut = kTRUE;
          else if(fReverseIsoCut == kTRUE)
            isocut = kFALSE;
        }
        break;
      case kNoIso:
        isocut = kTRUE;
        break;
      case kCustomIso:
      default:
        break;
    }

    if (isocut == kFALSE)
      continue;

    if (fApplyD0Cut) {
      Bool_t d0cut = kFALSE;
      const Track *mt = mu->BestTrk();
      if (!mt)
        continue;
      Double_t d0_real = 1e30;
      for(UInt_t i0 = 0; i0 < fVertices->GetEntries(); i0++) {
        Double_t pD0 = mt->D0Corrected(*fVertices->At(i0));
        if(TMath::Abs(pD0) < TMath::Abs(d0_real)) 
          d0_real = TMath::Abs(pD0);
      }
      if(d0_real < fD0Cut) d0cut = kTRUE;

      if     (fReverseD0Cut == kTRUE &&
              d0cut == kFALSE && d0_real < 0.05)
        d0cut = kTRUE;
      else if(fReverseD0Cut == kTRUE)
        d0cut = kFALSE;

      if (d0cut == kFALSE)
        continue;
    }

    // add good muon
    CleanMuons->Add(mu);
  }

  // sort according to pt
  CleanMuons->Sort();

  // add objects for other modules to use
  AddObjThisEvt(CleanMuons);  
}

//--------------------------------------------------------------------------------------------------
void MuonIDMod::SlaveBegin()
{
  // Run startup code on the computer (slave) doing the actual analysis. Here,
  // we just request the muon collection branch.

  ReqEventObject(fMuonBranchName, fMuons, kTRUE);

  if (fApplyD0Cut)
    ReqEventObject(fVertexName, fVertices, kTRUE);

  fMuonTools = new MuonTools;

  if (fMuonIDType.CompareTo("Tight") == 0) 
    fMuIDType = kTight;
  else if (fMuonIDType.CompareTo("Loose") == 0) 
    fMuIDType = kLoose;
  else if (fMuonIDType.CompareTo("NoId") == 0) 
    fMuIDType = kNoId;
  else if (fMuonIDType.CompareTo("Custom") == 0) {
    fMuIDType = kCustomId;
    SendError(kWarning, "SlaveBegin",
              "Custom muon identification is not yet implemented.");
  } else {
    SendError(kAbortAnalysis, "SlaveBegin",
              "The specified muon identification %s is not defined.",
              fMuonIDType.Data());
    return;
  }

  if (fMuonIsoType.CompareTo("TrackCalo") == 0)
    fMuIsoType = kTrackCalo;
  else if (fMuonIsoType.CompareTo("TrackCaloCombined") == 0)
    fMuIsoType = kTrackCaloCombined;
  else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
    fMuIsoType = kTrackCaloSliding;
  else if (fMuonIsoType.CompareTo("NoIso") == 0)
    fMuIsoType = kNoIso;
  else if (fMuonIsoType.CompareTo("Custom") == 0) {
    fMuIsoType = kCustomIso;
    SendError(kWarning, "SlaveBegin",
              "Custom muon isolation is not yet implemented.");
  } else {
    SendError(kAbortAnalysis, "SlaveBegin",
              "The specified muon isolation %s is not defined.",
              fMuonIsoType.Data());
    return;
  }

  if (fMuonClassType.CompareTo("All") == 0) 
    fMuClassType = kAll;
  else if (fMuonClassType.CompareTo("Global") == 0) 
    fMuClassType = kGlobal;
  else if (fMuonClassType.CompareTo("Standalone") == 0) 
    fMuClassType = kSta;
  else if (fMuonClassType.CompareTo("TrackerMuon") == 0) 
    fMuClassType = kTrackerMuon;
  else if (fMuonClassType.CompareTo("CaloMuon") == 0) 
    fMuClassType = kCaloMuon;
  else if (fMuonClassType.CompareTo("TrackerBased") == 0) 
    fMuClassType = kTrackerBased;
  else {
    SendError(kAbortAnalysis, "SlaveBegin",
              "The specified muon class %s is not defined.",
              fMuonClassType.Data());
    return;
  }
}
Revision:	1.24
Committed:	Fri Jul 17 12:11:36 2009 UTC (15 years, 9 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_011, Mit_010a, Mit_010
Changes since 1.23:	+18 -4 lines
Log Message:	Added kTrackerMuon kCaloMuon kTrackerBased.
#	User	Rev	Content
1	loizides	1.24	// $Id: MuonIDMod.cc,v 1.23 2009/06/15 15:00:21 loizides Exp $
2	loizides	1.1
3			#include "MitPhysics/Mods/interface/MuonIDMod.h"
4	loizides	1.6	#include "MitCommon/MathTools/interface/MathUtils.h"
5	loizides	1.23	#include "MitAna/DataTree/interface/MuonCol.h"
6			#include "MitAna/DataTree/interface/VertexCol.h"
7	loizides	1.6	#include "MitPhysics/Init/interface/ModNames.h"
8	loizides	1.1
9			using namespace mithep;
10
11			ClassImp(mithep::MuonIDMod)
12
13			//--------------------------------------------------------------------------------------------------
14			MuonIDMod::MuonIDMod(const char name, const char title) :
15			BaseMod(name,title),
16	loizides	1.6	fMuonBranchName(Names::gkMuonBrn),
17			fCleanMuonsName(ModNames::gkCleanMuonsName),
18	loizides	1.15	fVertexName("PrimaryVertexesBeamSpot"),
19	loizides	1.8	fMuonIDType("Loose"),
20			fMuonIsoType("TrackCaloSliding"),
21	loizides	1.6	fMuonClassType("Global"),
22	ceballos	1.2	fTrackIsolationCut(3.0),
23			fCaloIsolationCut(3.0),
24	loizides	1.8	fCombIsolationCut(5.0),
25	ceballos	1.3	fMuonPtMin(10),
26	loizides	1.21	fApplyD0Cut(kTRUE),
27	loizides	1.15	fD0Cut(0.025),
28	loizides	1.20	fReverseIsoCut(kFALSE),
29	ceballos	1.22	fReverseD0Cut(kFALSE),
30	loizides	1.15	fMuIDType(kIdUndef),
31			fMuIsoType(kIsoUndef),
32			fMuClassType(kClassUndef),
33	ceballos	1.14	fMuons(0),
34	loizides	1.15	fVertices(0),
35	loizides	1.20	fMuonTools(0)
36	loizides	1.1	{
37			// Constructor.
38			}
39
40			//--------------------------------------------------------------------------------------------------
41			void MuonIDMod::Process()
42			{
43			// Process entries of the tree.
44
45	loizides	1.21	LoadEventObject(fMuonBranchName, fMuons);
46			LoadEventObject(fVertexName, fVertices);
47	loizides	1.1
48	loizides	1.7	MuonOArr *CleanMuons = new MuonOArr;
49			CleanMuons->SetName(fCleanMuonsName);
50	loizides	1.1
51			for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
52	loizides	1.8	const Muon *mu = fMuons->At(i);
53	loizides	1.6
54			Bool_t pass = kFALSE;
55	loizides	1.19	Double_t pt = 0; // make sure pt is taken from the correct track!
56	loizides	1.6	switch (fMuClassType) {
57			case kAll:
58			pass = kTRUE;
59	loizides	1.19	if (mu->HasTrk())
60			pt = mu->Pt();
61	loizides	1.6	break;
62			case kGlobal:
63	loizides	1.19	pass = mu->HasGlobalTrk();
64	loizides	1.6	if (pass)
65			pt = mu->GlobalTrk()->Pt();
66			break;
67			case kSta:
68	loizides	1.19	pass = mu->HasStandaloneTrk();
69	loizides	1.6	if (pass)
70			pt = mu->StandaloneTrk()->Pt();
71			break;
72	loizides	1.24	case kTrackerMuon:
73			pass = mu->HasTrackerTrk() && mu->IsTrackerMuon();
74			if (pass)
75			pt = mu->TrackerTrk()->Pt();
76			break;
77			case kCaloMuon:
78			pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
79			if (pass)
80			pt = mu->TrackerTrk()->Pt();
81			break;
82			case kTrackerBased:
83	loizides	1.19	pass = mu->HasTrackerTrk();
84	loizides	1.6	if (pass)
85			pt = mu->TrackerTrk()->Pt();
86			break;
87			default:
88			break;
89	ceballos	1.5	}
90	loizides	1.6
91			if (!pass)
92			continue;
93
94			if (pt <= fMuonPtMin)
95			continue;
96
97			Bool_t idpass = kFALSE;
98			switch (fMuIDType) {
99			case kLoose:
100			idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationLoose) &&
101			fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityLoose);
102			break;
103			case kTight:
104			idpass = fMuonTools->IsGood(mu, MuonTools::kTMOneStationTight) &&
105			fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityTight);
106			break;
107	loizides	1.18	case kNoId:
108			idpass = kTRUE;
109			break;
110	loizides	1.6	default:
111			break;
112	ceballos	1.4	}
113	loizides	1.6
114			if (!idpass)
115			continue;
116
117	ceballos	1.22	Bool_t isocut = kFALSE;
118	loizides	1.6	switch (fMuIsoType) {
119			case kTrackCalo:
120	ceballos	1.22	isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
121	loizides	1.6	(mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
122			break;
123			case kTrackCaloCombined:
124	ceballos	1.22	isocut = (1.0 * mu->IsoR03SumPt() + 1.0 * mu->IsoR03EmEt() +
125	loizides	1.6	1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
126			break;
127			case kTrackCaloSliding:
128			{
129			Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
130			1.0 * mu->IsoR03EmEt() +
131			1.0 * mu->IsoR03HadEt();
132	loizides	1.19	if ((totalIso < (pt-10.0)*5.0/15.0 && pt <= 25) \|\|
133	ceballos	1.9	(totalIso < 5.0 && mu->Pt() > 25) \|\|
134			totalIso <= 0)
135	ceballos	1.22	isocut = kTRUE;
136
137			if (fReverseIsoCut == kTRUE &&
138			isocut == kFALSE && totalIso < 10)
139			isocut = kTRUE;
140			else if(fReverseIsoCut == kTRUE)
141			isocut = kFALSE;
142			}
143	loizides	1.6	break;
144	loizides	1.13	case kNoIso:
145	ceballos	1.22	isocut = kTRUE;
146	loizides	1.13	break;
147	loizides	1.6	case kCustomIso:
148			default:
149			break;
150	ceballos	1.4	}
151	ceballos	1.3
152	ceballos	1.22	if (isocut == kFALSE)
153	loizides	1.6	continue;
154
155	loizides	1.21	if (fApplyD0Cut) {
156	ceballos	1.22	Bool_t d0cut = kFALSE;
157	loizides	1.21	const Track *mt = mu->BestTrk();
158			if (!mt)
159			continue;
160			Double_t d0_real = 1e30;
161			for(UInt_t i0 = 0; i0 < fVertices->GetEntries(); i0++) {
162			Double_t pD0 = mt->D0Corrected(*fVertices->At(i0));
163			if(TMath::Abs(pD0) < TMath::Abs(d0_real))
164			d0_real = TMath::Abs(pD0);
165			}
166	ceballos	1.22	if(d0_real < fD0Cut) d0cut = kTRUE;
167
168			if (fReverseD0Cut == kTRUE &&
169			d0cut == kFALSE && d0_real < 0.05)
170			d0cut = kTRUE;
171			else if(fReverseD0Cut == kTRUE)
172			d0cut = kFALSE;
173
174			if (d0cut == kFALSE)
175	loizides	1.21	continue;
176	ceballos	1.14	}
177
178	loizides	1.6	// add good muon
179			CleanMuons->Add(mu);
180	loizides	1.1	}
181
182	loizides	1.10	// sort according to pt
183			CleanMuons->Sort();
184
185	loizides	1.6	// add objects for other modules to use
186	loizides	1.7	AddObjThisEvt(CleanMuons);
187	loizides	1.1	}
188
189			//--------------------------------------------------------------------------------------------------
190			void MuonIDMod::SlaveBegin()
191			{
192			// Run startup code on the computer (slave) doing the actual analysis. Here,
193	loizides	1.6	// we just request the muon collection branch.
194
195	loizides	1.21	ReqEventObject(fMuonBranchName, fMuons, kTRUE);
196
197			if (fApplyD0Cut)
198			ReqEventObject(fVertexName, fVertices, kTRUE);
199	loizides	1.1
200	loizides	1.6	fMuonTools = new MuonTools;
201	loizides	1.1
202	loizides	1.6	if (fMuonIDType.CompareTo("Tight") == 0)
203			fMuIDType = kTight;
204			else if (fMuonIDType.CompareTo("Loose") == 0)
205			fMuIDType = kLoose;
206	loizides	1.12	else if (fMuonIDType.CompareTo("NoId") == 0)
207			fMuIDType = kNoId;
208	loizides	1.6	else if (fMuonIDType.CompareTo("Custom") == 0) {
209			fMuIDType = kCustomId;
210			SendError(kWarning, "SlaveBegin",
211			"Custom muon identification is not yet implemented.");
212			} else {
213			SendError(kAbortAnalysis, "SlaveBegin",
214			"The specified muon identification %s is not defined.",
215			fMuonIDType.Data());
216			return;
217			}
218	loizides	1.1
219	loizides	1.6	if (fMuonIsoType.CompareTo("TrackCalo") == 0)
220			fMuIsoType = kTrackCalo;
221			else if (fMuonIsoType.CompareTo("TrackCaloCombined") == 0)
222			fMuIsoType = kTrackCaloCombined;
223			else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
224			fMuIsoType = kTrackCaloSliding;
225			else if (fMuonIsoType.CompareTo("NoIso") == 0)
226			fMuIsoType = kNoIso;
227			else if (fMuonIsoType.CompareTo("Custom") == 0) {
228			fMuIsoType = kCustomIso;
229			SendError(kWarning, "SlaveBegin",
230			"Custom muon isolation is not yet implemented.");
231			} else {
232			SendError(kAbortAnalysis, "SlaveBegin",
233			"The specified muon isolation %s is not defined.",
234			fMuonIsoType.Data());
235			return;
236			}
237	loizides	1.1
238	loizides	1.6	if (fMuonClassType.CompareTo("All") == 0)
239			fMuClassType = kAll;
240			else if (fMuonClassType.CompareTo("Global") == 0)
241			fMuClassType = kGlobal;
242			else if (fMuonClassType.CompareTo("Standalone") == 0)
243			fMuClassType = kSta;
244	loizides	1.24	else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
245			fMuClassType = kTrackerMuon;
246			else if (fMuonClassType.CompareTo("CaloMuon") == 0)
247			fMuClassType = kCaloMuon;
248			else if (fMuonClassType.CompareTo("TrackerBased") == 0)
249			fMuClassType = kTrackerBased;
250	loizides	1.6	else {
251			SendError(kAbortAnalysis, "SlaveBegin",
252			"The specified muon class %s is not defined.",
253			fMuonClassType.Data());
254			return;
255			}
256	loizides	1.1	}