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.9 by ceballos, Fri Nov 28 10:31:38 2008 UTC vs.
Revision 1.46 by ceballos, Tue Apr 5 06:37:07 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 <  fMuonIDType("Loose"),
20 <  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("TrackCaloSliding"),
27    fMuonClassType("Global"),  
28    fTrackIsolationCut(3.0),
29    fCaloIsolationCut(3.0),
30 <  fCombIsolationCut(5.0),
30 >  fCombIsolationCut(0.15),
31    fMuonPtMin(10),
32 <  fMuons(0)
32 >  fApplyD0Cut(kTRUE),
33 >  fApplyDZCut(kTRUE),
34 >  fD0Cut(0.020),
35 >  fDZCut(0.20),
36 >  fWhichVertex(-1),
37 >  fEtaCut(2.4),
38 >  fMuIDType(kIdUndef),
39 >  fMuIsoType(kIsoUndef),
40 >  fMuClassType(kClassUndef),
41 >  fMuons(0),
42 >  fVertices(0),
43 >  fBeamSpot(0),
44 >  fTracks(0),
45 >  fPFCandidates(0),
46 >  fNonIsolatedMuons(0),
47 >  fNonIsolatedElectrons(0),
48 >  fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
49 >  fPileupEnergyDensity(0)
50   {
51    // Constructor.
52   }
53  
28
54   //--------------------------------------------------------------------------------------------------
55   void MuonIDMod::Process()
56   {
57    // Process entries of the tree.
58  
59 <  LoadBranch(fMuonBranchName);
60 <
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 = -1; // 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 <        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->GlobalTrk() != 0);
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->StandaloneTrk() != 0);
103 <        if (pass)
104 <          pt = mu->StandaloneTrk()->Pt();
105 <        break;
106 <      case kTrackerOnly:
107 <        pass = (mu->TrackerTrk() != 0);
108 <        if (pass)
109 <          pt = mu->TrackerTrk()->Pt();
102 >        pass = mu->HasStandaloneTrk();
103 >        if (pass) {
104 >          pt  = mu->StandaloneTrk()->Pt();
105 >          eta = TMath::Abs(mu->StandaloneTrk()->Eta());
106 >        }
107 >        break;
108 >      case kTrackerMuon:
109 >        pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
110 >               mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
111 >        if (pass) {
112 >          pt  = mu->TrackerTrk()->Pt();
113 >          eta = TMath::Abs(mu->TrackerTrk()->Eta());
114 >        }
115 >        break;
116 >      case kCaloMuon:
117 >        pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
118 >        if (pass) {
119 >          pt  = mu->TrackerTrk()->Pt();
120 >          eta = TMath::Abs(mu->TrackerTrk()->Eta());
121 >        }
122 >        break;
123 >      case kTrackerBased:
124 >        pass = mu->HasTrackerTrk();
125 >        if (pass) {
126 >          pt  = mu->TrackerTrk()->Pt();
127 >          eta = TMath::Abs(mu->TrackerTrk()->Eta());
128 >        }
129          break;
130        default:
131          break;
# Line 71 | 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;
194          break;
195        default:
196          break;
# Line 88 | 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 < (mu->Pt()-10.0)*5.0/15.0 && mu->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 beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
243 >          if(beta == 0) beta = 1.0;
244 >          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 0, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
245 >          if (totalIso < (mu->Pt()*fCombIsolationCut) )
246 >            isocut = kTRUE;
247 >        }
248 >        break;
249 >      case kPFIsoNoL:
250 >        {
251 >          fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
252 >          fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
253 >
254 >          Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
255 >          if(beta == 0) beta = 1.0;
256 >          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 3, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
257 >          if (totalIso < (mu->Pt()*fCombIsolationCut) )
258 >            isocut = kTRUE;
259 >        }
260 >        break;
261 >      case kNoIso:
262 >        isocut = kTRUE;
263          break;
112        case kNoIso:
113          isopass = kTRUE;
114          break;
264        case kCustomIso:
265        default:
266          break;
267      }
268  
269 <    if (!isopass)
269 >    if (isocut == kFALSE)
270        continue;
271  
272 +    // apply d0 cut
273 +    if (fApplyD0Cut) {
274 +      Bool_t passD0cut = kTRUE;
275 +      if(fD0Cut < 0.05) { // trick to change the signal region cut
276 +        if      (mu->Pt() >  20.0) fD0Cut = 0.02;
277 +        else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
278 +      }
279 +      if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
280 +      else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
281 +      if (!passD0cut)
282 +        continue;
283 +    }
284 +
285 +    // apply dz cut
286 +    if (fApplyDZCut) {
287 +      Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
288 +      if (!passDZcut)
289 +        continue;
290 +    }
291 +
292      // add good muon
293      CleanMuons->Add(mu);
294    }
295  
296 +  // sort according to pt
297 +  CleanMuons->Sort();
298 +
299    // add objects for other modules to use
300    AddObjThisEvt(CleanMuons);  
301   }
302  
131
303   //--------------------------------------------------------------------------------------------------
304   void MuonIDMod::SlaveBegin()
305   {
306    // Run startup code on the computer (slave) doing the actual analysis. Here,
307    // we just request the muon collection branch.
308  
309 <  ReqBranch(fMuonBranchName, fMuons);
310 <
311 <  fMuonTools = new MuonTools;
309 >   // In this case we cannot have a branch
310 >  if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
311 >    ReqEventObject(fMuonBranchName, fMuons, kTRUE);
312 >  }
313 >  ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
314 >  ReqEventObject(fTrackName, fTracks, kTRUE);
315 >  ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
316 >  if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0) {
317 >    ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
318 >  }
319  
320 <  if (fMuonIDType.CompareTo("Tight") == 0)
320 >  if (fMuonIDType.CompareTo("WMuId") == 0)
321 >    fMuIDType = kWMuId;
322 >  else if (fMuonIDType.CompareTo("ZMuId") == 0)
323 >    fMuIDType = kZMuId;
324 >  else if (fMuonIDType.CompareTo("Tight") == 0)
325      fMuIDType = kTight;
326    else if (fMuonIDType.CompareTo("Loose") == 0)
327      fMuIDType = kLoose;
328 +  else if (fMuonIDType.CompareTo("WWMuId") == 0)
329 +    fMuIDType = kWWMuId;
330 +  else if (fMuonIDType.CompareTo("NoId") == 0)
331 +    fMuIDType = kNoId;
332    else if (fMuonIDType.CompareTo("Custom") == 0) {
333      fMuIDType = kCustomId;
334      SendError(kWarning, "SlaveBegin",
# Line 160 | Line 346 | void MuonIDMod::SlaveBegin()
346      fMuIsoType = kTrackCaloCombined;
347    else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
348      fMuIsoType = kTrackCaloSliding;
349 +  else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
350 +    fMuIsoType = kTrackCaloSlidingNoCorrection;
351 +  else if (fMuonIsoType.CompareTo("PFIso") == 0)
352 +    fMuIsoType = kPFIso;
353 +  else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
354 +    fMuIsoType = kPFIsoNoL;
355    else if (fMuonIsoType.CompareTo("NoIso") == 0)
356      fMuIsoType = kNoIso;
357    else if (fMuonIsoType.CompareTo("Custom") == 0) {
# Line 179 | Line 371 | void MuonIDMod::SlaveBegin()
371      fMuClassType = kGlobal;
372    else if (fMuonClassType.CompareTo("Standalone") == 0)
373      fMuClassType = kSta;
374 <  else if (fMuonClassType.CompareTo("TrackerOnly") == 0)
375 <    fMuClassType = kTrackerOnly;
374 >  else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
375 >    fMuClassType = kTrackerMuon;
376 >  else if (fMuonClassType.CompareTo("CaloMuon") == 0)
377 >    fMuClassType = kCaloMuon;
378 >  else if (fMuonClassType.CompareTo("TrackerBased") == 0)
379 >    fMuClassType = kTrackerBased;
380    else {
381      SendError(kAbortAnalysis, "SlaveBegin",
382                "The specified muon class %s is not defined.",

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines