ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
(Generate patch)

Comparing UserCode/MitPhysics/Mods/src/MuonIDMod.cc (file contents):
Revision 1.20 by loizides, Mon May 11 08:02:43 2009 UTC vs.
Revision 1.49 by ceballos, Thu May 12 23:15:48 2011 UTC

# Line 2 | Line 2
2  
3   #include "MitPhysics/Mods/interface/MuonIDMod.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;
# Line 13 | Line 16 | ClassImp(mithep::MuonIDMod)
16    BaseMod(name,title),
17    fMuonBranchName(Names::gkMuonBrn),
18    fCleanMuonsName(ModNames::gkCleanMuonsName),  
19 <  fVertexName("PrimaryVertexesBeamSpot"),
20 <  fMuonIDType("Loose"),
21 <  fMuonIsoType("TrackCaloSliding"),  
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("PFIso"),
27    fMuonClassType("Global"),  
28    fTrackIsolationCut(3.0),
29    fCaloIsolationCut(3.0),
30 <  fCombIsolationCut(5.0),
30 >  fCombIsolationCut(-1.0),
31    fMuonPtMin(10),
32 <  fD0Cut(0.025),
33 <  fReverseIsoCut(kFALSE),
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 <  fMuonTools(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   }
# Line 38 | Line 56 | void MuonIDMod::Process()
56   {
57    // Process entries of the tree.
58  
59 <  LoadBranch(fMuonBranchName);
60 <  LoadBranch(fVertexName);
61 <
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      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!
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();
85 >        if (mu->HasTrk()) {
86 >          pt  = mu->Pt();
87 >          eta = TMath::Abs(mu->Eta());
88 >        }
89          break;
90        case kGlobal:
91 <        pass = mu->HasGlobalTrk();
92 <        if (pass)
93 <          pt = mu->GlobalTrk()->Pt();
94 <        break;
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();
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 kTrackerOnly:
123 >      case kTrackerBased:
124          pass = mu->HasTrackerTrk();
125 <        if (pass)
126 <          pt = mu->TrackerTrk()->Pt();
125 >        if (pass) {
126 >          pt  = mu->TrackerTrk()->Pt();
127 >          eta = TMath::Abs(mu->TrackerTrk()->Eta());
128 >        }
129          break;
130        default:
131          break;
# Line 80 | Line 137 | void MuonIDMod::Process()
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 = fMuonTools->IsGood(mu, MuonTools::kTMOneStationLoose) &&
169 <                 fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityLoose);
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 = fMuonTools->IsGood(mu, MuonTools::kTMOneStationTight) &&
177 <                 fMuonTools->IsGood(mu, MuonTools::kTM2DCompatibilityTight);
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;
# Line 100 | Line 199 | void MuonIDMod::Process()
199      if (!idpass)
200        continue;
201  
202 <    Bool_t isopass = kFALSE;
202 >    Bool_t isocut = kFALSE;
203      switch (fMuIsoType) {
204        case kTrackCalo:
205 <        isopass = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
205 >        isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
206            (mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
207          break;
208        case kTrackCaloCombined:
209 <        isopass = (1.0 * mu->IsoR03SumPt() + 1.0 * mu->IsoR03EmEt() +
210 <                   1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
209 >        isocut = (1.0 * mu->IsoR03SumPt() +
210 >                  1.0 * mu->IsoR03EmEt()  +
211 >                  1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
212          break;
213        case kTrackCaloSliding:
214          {
215 <          Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
216 <                              1.0 * mu->IsoR03EmEt() +
217 <                              1.0 * mu->IsoR03HadEt();
218 <          if ((totalIso < (pt-10.0)*5.0/15.0 && pt <= 25) ||
219 <              (totalIso < 5.0 && mu->Pt() > 25) ||
220 <               totalIso <= 0)
221 <            isopass = kTRUE;
222 <        }
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 pfIsoCutValue = 9999;
243 >          if(fCombIsolationCut > 0){
244 >            pfIsoCutValue = fCombIsolationCut;
245 >          } else {
246 >            if (mu->Pt() > 20) {
247 >              pfIsoCutValue = 0.17;
248 >            } else {
249 >              pfIsoCutValue = 0.13;
250 >            }
251 >          }
252 >          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0),0.2, 1.0, 0.4, 0.0);
253 >          if (totalIso < (mu->Pt()*pfIsoCutValue) )
254 >            isocut = kTRUE;
255 >        }
256 >        break;
257 >      case kPFIsoNoL:
258 >        {
259 >          fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
260 >          fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
261 >
262 >          Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
263 >          if(beta == 0) beta = 1.0;
264 >          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), fNonIsolatedMuons, fNonIsolatedElectrons, 0.2, 1.0, 0.4, 0.0, 3);
265 >          if (totalIso < (mu->Pt()*fCombIsolationCut) )
266 >            isocut = kTRUE;
267 >        }
268          break;
269        case kNoIso:
270 <        isopass = kTRUE;
270 >        isocut = kTRUE;
271          break;
272        case kCustomIso:
273        default:
274          break;
275      }
276  
277 <    if ((isopass == kFALSE && fReverseIsoCut == kFALSE) ||
133 <        (isopass == kTRUE  && fReverseIsoCut == kTRUE))
277 >    if (isocut == kFALSE)
278        continue;
279  
280 <    // d0 cut only for muons that have a track
281 <    const Track *mt = mu->BestTrk();
282 <    if (!mt)
283 <      continue;
284 <    Double_t d0_real = 1e30;
285 <    for(UInt_t i0 = 0; i0 < fVertices->GetEntries(); i0++) {
286 <      Double_t pD0 = mt->D0Corrected(*fVertices->At(i0));
287 <      if(TMath::Abs(pD0) < TMath::Abs(d0_real))
288 <        d0_real = TMath::Abs(pD0);
280 >    // apply d0 cut
281 >    if (fApplyD0Cut) {
282 >      Bool_t passD0cut = kTRUE;
283 >      if(fD0Cut < 0.05) { // trick to change the signal region cut
284 >        if      (mu->Pt() >  20.0) fD0Cut = 0.02;
285 >        else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
286 >      }
287 >      if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
288 >      else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
289 >      if (!passD0cut)
290 >        continue;
291 >    }
292 >
293 >    // apply dz cut
294 >    if (fApplyDZCut) {
295 >      Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
296 >      if (!passDZcut)
297 >        continue;
298      }
146    if(d0_real >= fD0Cut)
147      continue;
299  
300      // add good muon
301      CleanMuons->Add(mu);
# Line 163 | Line 314 | void MuonIDMod::SlaveBegin()
314    // Run startup code on the computer (slave) doing the actual analysis. Here,
315    // we just request the muon collection branch.
316  
317 <  ReqBranch(fMuonBranchName, fMuons);
318 <  ReqBranch(fVertexName, fVertices);
319 <
320 <  fMuonTools = new MuonTools;
317 >   // In this case we cannot have a branch
318 >  if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
319 >    ReqEventObject(fMuonBranchName, fMuons, kTRUE);
320 >  }
321 >  ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
322 >  ReqEventObject(fTrackName, fTracks, kTRUE);
323 >  ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
324 >  if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0) {
325 >    ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
326 >  }
327  
328 <  if (fMuonIDType.CompareTo("Tight") == 0)
328 >  if (fMuonIDType.CompareTo("WMuId") == 0)
329 >    fMuIDType = kWMuId;
330 >  else if (fMuonIDType.CompareTo("ZMuId") == 0)
331 >    fMuIDType = kZMuId;
332 >  else if (fMuonIDType.CompareTo("Tight") == 0)
333      fMuIDType = kTight;
334    else if (fMuonIDType.CompareTo("Loose") == 0)
335      fMuIDType = kLoose;
336 +  else if (fMuonIDType.CompareTo("WWMuId") == 0)
337 +    fMuIDType = kWWMuId;
338    else if (fMuonIDType.CompareTo("NoId") == 0)
339      fMuIDType = kNoId;
340    else if (fMuonIDType.CompareTo("Custom") == 0) {
# Line 191 | Line 354 | void MuonIDMod::SlaveBegin()
354      fMuIsoType = kTrackCaloCombined;
355    else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
356      fMuIsoType = kTrackCaloSliding;
357 +  else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
358 +    fMuIsoType = kTrackCaloSlidingNoCorrection;
359 +  else if (fMuonIsoType.CompareTo("PFIso") == 0)
360 +    fMuIsoType = kPFIso;
361 +  else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
362 +    fMuIsoType = kPFIsoNoL;
363    else if (fMuonIsoType.CompareTo("NoIso") == 0)
364      fMuIsoType = kNoIso;
365    else if (fMuonIsoType.CompareTo("Custom") == 0) {
# Line 210 | Line 379 | void MuonIDMod::SlaveBegin()
379      fMuClassType = kGlobal;
380    else if (fMuonClassType.CompareTo("Standalone") == 0)
381      fMuClassType = kSta;
382 <  else if (fMuonClassType.CompareTo("TrackerOnly") == 0)
383 <    fMuClassType = kTrackerOnly;
382 >  else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
383 >    fMuClassType = kTrackerMuon;
384 >  else if (fMuonClassType.CompareTo("CaloMuon") == 0)
385 >    fMuClassType = kCaloMuon;
386 >  else if (fMuonClassType.CompareTo("TrackerBased") == 0)
387 >    fMuClassType = kTrackerBased;
388    else {
389      SendError(kAbortAnalysis, "SlaveBegin",
390                "The specified muon class %s is not defined.",

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines