CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PFMuonAlgo.cc
Go to the documentation of this file.
1 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
2 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
3 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
4 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementTrack.h"
5 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementGsfTrack.h"
6 #include "DataFormats/TrackReco/interface/Track.h"
7 #include "DataFormats/MuonReco/interface/Muon.h"
8 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
9 #include "DataFormats/MuonReco/interface/MuonCocktails.h"
10 #include <iostream>
11 
12 
13 using namespace std;
14 using namespace reco;
15 using namespace boost;
16 
17 PFMuonAlgo::PFMuonAlgo() {
18  pfCosmicsMuonCleanedCandidates_ = std::auto_ptr<reco::PFCandidateCollection>(new reco::PFCandidateCollection);
19  pfCleanedTrackerAndGlobalMuonCandidates_= std::auto_ptr<reco::PFCandidateCollection>(new reco::PFCandidateCollection);
20  pfFakeMuonCleanedCandidates_= std::auto_ptr<reco::PFCandidateCollection>(new reco::PFCandidateCollection);
21  pfPunchThroughMuonCleanedCandidates_= std::auto_ptr<reco::PFCandidateCollection>(new reco::PFCandidateCollection);
22  pfPunchThroughHadronCleanedCandidates_= std::auto_ptr<reco::PFCandidateCollection>(new reco::PFCandidateCollection);
23  pfAddedMuonCandidates_= std::auto_ptr<reco::PFCandidateCollection>(new reco::PFCandidateCollection);
24 
25 }
26 
27 
28 
29 void PFMuonAlgo::setParameters(const edm::ParameterSet& iConfig )
30 {
31 
32  if(iConfig.exists("maxDPtOPt"))
33  maxDPtOPt_ = iConfig.getParameter<double>("maxDPtOPt");
34  else
35  maxDPtOPt_=1.0;
36 
37  if(iConfig.exists("minTrackerHits"))
38  minTrackerHits_ = iConfig.getParameter<int>("minTrackerHits");
39  else
40  minTrackerHits_ = 8;
41 
42  if(iConfig.exists("minPixelHits"))
43  minPixelHits_ = iConfig.getParameter<int>("minPixelHits");
44  else
45  minPixelHits_ = 1;
46 
47  if(iConfig.exists("trackQuality"))
48  trackQuality_ = reco::TrackBase::qualityByName(iConfig.getParameter<std::string>("trackQuality"));
49  else
50  trackQuality_ = reco::TrackBase::qualityByName("highPurity");
51 
52  if(iConfig.exists("ptErrorScale"))
53  errorCompScale_ = iConfig.getParameter<double>("ptErrorScale");
54  else
55  errorCompScale_ = 4.;
56 
57  if(iConfig.exists("eventFractionForCleaning"))
58  eventFractionCleaning_ = iConfig.getParameter<double>("eventFractionForCleaning");
59  else
60  eventFractionCleaning_ = 0.75;
61 
62  if(iConfig.exists("dzPV"))
63  dzPV_ = iConfig.getParameter<double>("dzPV");
64  else
65  dzPV_ = 0.2;
66 
67  if(iConfig.exists("postMuonCleaning"))
68  postCleaning_ = iConfig.getParameter<bool>("postMuonCleaning");
69  else
70  postCleaning_ = false; //Disable by default (for HLT)
71 
72  if(iConfig.exists("minPtForPostCleaning"))
73  minPostCleaningPt_ = iConfig.getParameter<double>("minPtForPostCleaning");
74  else
75  minPostCleaningPt_ = 20.;
76 
77  if(iConfig.exists("eventFactorForCosmics"))
78  eventFactorCosmics_ = iConfig.getParameter<double>("eventFactorForCosmics");
79  else
80  eventFactorCosmics_ = 10.;
81 
82 
83  if(iConfig.exists("metSignificanceForCleaning"))
84  metSigForCleaning_ = iConfig.getParameter<double>("metSignificanceForCleaning");
85  else
86  metSigForCleaning_ = 3.;
87 
88  if(iConfig.exists("metSignificanceForRejection"))
89  metSigForRejection_ = iConfig.getParameter<double>("metSignificanceForRejection");
90  else
91  metSigForRejection_ = 4.;
92 
93  if(iConfig.exists("metFactorForCleaning"))
94  metFactorCleaning_ = iConfig.getParameter<double>("metFactorForCleaning");
95  else
96  metFactorCleaning_ = 4.;
97 
98  if(iConfig.exists("eventFractionForRejection"))
99  eventFractionRejection_ = iConfig.getParameter<double>("eventFractionForRejection");
100  else
101  eventFractionRejection_ = 0.75;
102 
103  if(iConfig.exists("metFactorForRejection"))
104  metFactorRejection_ = iConfig.getParameter<double>("metFactorForRejection");
105  else
106  metFactorRejection_ =4.;
107 
108  if(iConfig.exists("metFactorForHighEta"))
109  metFactorHighEta_ = iConfig.getParameter<double>("metFactorForHighEta");
110  else
111  metFactorHighEta_=4;
112 
113  if(iConfig.exists("ptFactorForHighEta"))
114  ptFactorHighEta_ = iConfig.getParameter<double>("ptFactorForHighEta");
115  else
116  ptFactorHighEta_ = 2.;
117 
118  if(iConfig.exists("metFactorForFakes"))
119  metFactorFake_ = iConfig.getParameter<double>("metFactorForFakes");
120  else
121  metFactorFake_ = 4.;
122 
123  if(iConfig.exists("minMomentumForPunchThrough"))
124  minPunchThroughMomentum_ = iConfig.getParameter<double>("minMomentumForPunchThrough");
125  else
126  minPunchThroughMomentum_=100.;
127 
128  if(iConfig.exists("minEnergyForPunchThrough"))
129  minPunchThroughEnergy_ = iConfig.getParameter<double>("minEnergyForPunchThrough");
130  else
131  minPunchThroughEnergy_ = 100.;
132 
133  if(iConfig.exists("punchThroughFactor"))
134  punchThroughFactor_ = iConfig.getParameter<double>("punchThroughFactor");
135  else
136  punchThroughFactor_ = 3.;
137 
138  if(iConfig.exists("punchThroughMETFactor"))
139  punchThroughMETFactor_ = iConfig.getParameter<double>("punchThroughMETFactor");
140  else
141  punchThroughMETFactor_ = 4.;
142 
143  if(iConfig.exists("cosmicRejectionDistance"))
144  cosmicRejDistance_ = iConfig.getParameter<double>("cosmicRejectionDistance");
145  else
146  cosmicRejDistance_ = 1.0;
147 }
148 
149 
150 
151 
152 bool
153 PFMuonAlgo::isMuon( const reco::PFBlockElement& elt ) {
154 
155  const reco::PFBlockElementTrack* eltTrack
156  = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
157 
158  assert ( eltTrack );
159  reco::MuonRef muonRef = eltTrack->muonRef();
160 
161  return isMuon(muonRef);
162 
163 }
164 
165 bool
166 PFMuonAlgo::isLooseMuon( const reco::PFBlockElement& elt ) {
167 
168 
169  const reco::PFBlockElementTrack* eltTrack
170  = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
171 
172 
173 
174  assert ( eltTrack );
175 
176 
177  reco::MuonRef muonRef = eltTrack->muonRef();
178 
179 
180  return isLooseMuon(muonRef);
181 
182 }
183 
184 bool
185 PFMuonAlgo::isGlobalTightMuon( const reco::PFBlockElement& elt ) {
186 
187  const reco::PFBlockElementTrack* eltTrack
188  = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
189 
190  assert ( eltTrack );
191  reco::MuonRef muonRef = eltTrack->muonRef();
192 
193  return isGlobalTightMuon(muonRef);
194 
195 }
196 
197 bool
198 PFMuonAlgo::isGlobalLooseMuon( const reco::PFBlockElement& elt ) {
199 
200  const reco::PFBlockElementTrack* eltTrack
201  = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
202 
203  assert ( eltTrack );
204  reco::MuonRef muonRef = eltTrack->muonRef();
205 
206  return isGlobalLooseMuon(muonRef);
207 
208 }
209 
210 bool
211 PFMuonAlgo::isTrackerTightMuon( const reco::PFBlockElement& elt ) {
212 
213  const reco::PFBlockElementTrack* eltTrack
214  = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
215 
216  assert ( eltTrack );
217  reco::MuonRef muonRef = eltTrack->muonRef();
218 
219  return isTrackerTightMuon(muonRef);
220 
221 }
222 
223 bool
224 PFMuonAlgo::isIsolatedMuon( const reco::PFBlockElement& elt ) {
225 
226  const reco::PFBlockElementTrack* eltTrack
227  = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
228 
229  assert ( eltTrack );
230  reco::MuonRef muonRef = eltTrack->muonRef();
231 
232  return isIsolatedMuon(muonRef);
233 
234 }
235 
236 bool
237 PFMuonAlgo::isMuon(const reco::MuonRef& muonRef ){
238 
239  return isGlobalTightMuon(muonRef) || isTrackerTightMuon(muonRef) || isIsolatedMuon(muonRef);
240 }
241 
242 bool
243 PFMuonAlgo::isLooseMuon(const reco::MuonRef& muonRef ){
244 
245  return isGlobalLooseMuon(muonRef) || isTrackerLooseMuon(muonRef);
246 
247 }
248 
249 bool
250 PFMuonAlgo::isGlobalTightMuon( const reco::MuonRef& muonRef ) {
251 
252  if ( !muonRef.isNonnull() ) return false;
253 
254  if ( !muonRef->isGlobalMuon() ) return false;
255  if ( !muonRef->isStandAloneMuon() ) return false;
256 
257 
258  if ( muonRef->isTrackerMuon() ) {
259 
260  bool result = muon::isGoodMuon(*muonRef,muon::GlobalMuonPromptTight);
261 
262  bool isTM2DCompatibilityTight = muon::isGoodMuon(*muonRef,muon::TM2DCompatibilityTight);
263  int nMatches = muonRef->numberOfMatches();
264  bool quality = nMatches > 2 || isTM2DCompatibilityTight;
265 
266  return result && quality;
267 
268  } else {
269 
270  reco::TrackRef standAloneMu = muonRef->standAloneMuon();
271 
272  // No tracker muon -> Request a perfect stand-alone muon, or an even better global muon
273  bool result = false;
274 
275  // Check the quality of the stand-alone muon :
276  // good chi**2 and large number of hits and good pt error
277  if ( ( standAloneMu->hitPattern().numberOfValidMuonDTHits() < 22 &&
278  standAloneMu->hitPattern().numberOfValidMuonCSCHits() < 15 ) ||
279  standAloneMu->normalizedChi2() > 10. ||
280  standAloneMu->ptError()/standAloneMu->pt() > 0.20 ) {
281  result = false;
282  } else {
283 
284  reco::TrackRef combinedMu = muonRef->combinedMuon();
285  reco::TrackRef trackerMu = muonRef->track();
286 
287  // If the stand-alone muon is good, check the global muon
288  if ( combinedMu->normalizedChi2() > standAloneMu->normalizedChi2() ) {
289  // If the combined muon is worse than the stand-alone, it
290  // means that either the corresponding tracker track was not
291  // reconstructed, or that the sta muon comes from a late
292  // pion decay (hence with a momentum smaller than the track)
293  // Take the stand-alone muon only if its momentum is larger
294  // than that of the track
295  result = standAloneMu->pt() > trackerMu->pt() ;
296  } else {
297  // If the combined muon is better (and good enough), take the
298  // global muon
299  result =
300  combinedMu->ptError()/combinedMu->pt() <
301  std::min(0.20,standAloneMu->ptError()/standAloneMu->pt());
302  }
303  }
304 
305  return result;
306  }
307 
308  return false;
309 
310 }
311 
312 bool
313 PFMuonAlgo::isTrackerTightMuon( const reco::MuonRef& muonRef ) {
314 
315  if ( !muonRef.isNonnull() ) return false;
316 
317  if(!muonRef->isTrackerMuon()) return false;
318 
319  reco::TrackRef trackerMu = muonRef->track();
320  const reco::Track& track = *trackerMu;
321 
322  unsigned nTrackerHits = track.hitPattern().numberOfValidTrackerHits();
323 
324  if(nTrackerHits<=12) return false;
325 
326  bool isAllArbitrated = muon::isGoodMuon(*muonRef,muon::AllArbitrated);
327 
328  bool isTM2DCompatibilityTight = muon::isGoodMuon(*muonRef,muon::TM2DCompatibilityTight);
329 
330  if(!isAllArbitrated || !isTM2DCompatibilityTight) return false;
331 
332  if((trackerMu->ptError()/trackerMu->pt() > 0.10)){
333  //std::cout<<" PT ERROR > 10 % "<< trackerMu->pt() <<std::endl;
334  return false;
335  }
336  return true;
337 
338 }
339 
340 bool
341 PFMuonAlgo::isGlobalLooseMuon( const reco::MuonRef& muonRef ) {
342 
343  if ( !muonRef.isNonnull() ) return false;
344  if ( !muonRef->isGlobalMuon() ) return false;
345  if ( !muonRef->isStandAloneMuon() ) return false;
346 
347  reco::TrackRef standAloneMu = muonRef->standAloneMuon();
348  reco::TrackRef combinedMu = muonRef->combinedMuon();
349  reco::TrackRef trackerMu = muonRef->track();
350 
351  unsigned nMuonHits =
352  standAloneMu->hitPattern().numberOfValidMuonDTHits() +
353  2*standAloneMu->hitPattern().numberOfValidMuonCSCHits();
354 
355  bool quality = false;
356 
357  if ( muonRef->isTrackerMuon() ){
358 
359  bool result = combinedMu->normalizedChi2() < 100.;
360 
361  bool laststation =
362  muon::isGoodMuon(*muonRef,muon::TMLastStationAngTight);
363 
364  int nMatches = muonRef->numberOfMatches();
365 
366  quality = laststation && nMuonHits > 12 && nMatches > 1;
367 
368  return result && quality;
369 
370  }
371  else{
372 
373  // Check the quality of the stand-alone muon :
374  // good chi**2 and large number of hits and good pt error
375  if ( nMuonHits <=15 ||
376  standAloneMu->normalizedChi2() > 10. ||
377  standAloneMu->ptError()/standAloneMu->pt() > 0.20 ) {
378  quality = false;
379  }
380  else {
381  // If the stand-alone muon is good, check the global muon
382  if ( combinedMu->normalizedChi2() > standAloneMu->normalizedChi2() ) {
383  // If the combined muon is worse than the stand-alone, it
384  // means that either the corresponding tracker track was not
385  // reconstructed, or that the sta muon comes from a late
386  // pion decay (hence with a momentum smaller than the track)
387  // Take the stand-alone muon only if its momentum is larger
388  // than that of the track
389 
390  // Note that here we even take the standAlone if it has a smaller pT, in contrast to GlobalTight
391  if(standAloneMu->pt() > trackerMu->pt() || combinedMu->normalizedChi2()<5.) quality = true;
392  }
393  else {
394  // If the combined muon is better (and good enough), take the
395  // global muon
396  if(combinedMu->ptError()/combinedMu->pt() < std::min(0.20,standAloneMu->ptError()/standAloneMu->pt()))
397  quality = true;
398 
399  }
400  }
401  }
402 
403 
404  return quality;
405 
406 }
407 
408 
409 bool
410 PFMuonAlgo::isTrackerLooseMuon( const reco::MuonRef& muonRef ) {
411 
412  if ( !muonRef.isNonnull() ) return false;
413  if(!muonRef->isTrackerMuon()) return false;
414 
415  reco::TrackRef trackerMu = muonRef->track();
416 
417  if(trackerMu->ptError()/trackerMu->pt() > 0.20) return false;
418 
419  // this doesn't seem to be necessary on the small samples looked at, but keep it around as insurance
420  if(trackerMu->pt()>20.) return false;
421 
422  bool isAllArbitrated = muon::isGoodMuon(*muonRef,muon::AllArbitrated);
423  bool isTMLastStationAngTight = muon::isGoodMuon(*muonRef,muon::TMLastStationAngTight);
424 
425  bool quality = isAllArbitrated && isTMLastStationAngTight;
426 
427  return quality;
428 
429 }
430 
431 bool
432 PFMuonAlgo::isIsolatedMuon( const reco::MuonRef& muonRef ){
433 
434 
435  if ( !muonRef.isNonnull() ) return false;
436  if ( !muonRef->isIsolationValid() ) return false;
437 
438  // Isolated Muons which are missed by standard cuts are nearly always global+tracker
439  if ( !muonRef->isGlobalMuon() ) return false;
440 
441  // If it's not a tracker muon, only take it if there are valid muon hits
442 
443  reco::TrackRef standAloneMu = muonRef->standAloneMuon();
444 
445  if ( !muonRef->isTrackerMuon() ){
446  if(standAloneMu->hitPattern().numberOfValidMuonDTHits() == 0 &&
447  standAloneMu->hitPattern().numberOfValidMuonCSCHits() ==0) return false;
448  }
449 
450  // for isolation, take the smallest pt available to reject fakes
451 
452  reco::TrackRef combinedMu = muonRef->combinedMuon();
453  double smallestMuPt = combinedMu->pt();
454 
455  if(standAloneMu->pt()<smallestMuPt) smallestMuPt = standAloneMu->pt();
456 
457  if(muonRef->isTrackerMuon())
458  {
459  reco::TrackRef trackerMu = muonRef->track();
460  if(trackerMu->pt() < smallestMuPt) smallestMuPt= trackerMu->pt();
461  }
462 
463  double sumPtR03 = muonRef->isolationR03().sumPt;
464  double emEtR03 = muonRef->isolationR03().emEt;
465  double hadEtR03 = muonRef->isolationR03().hadEt;
466 
467  double relIso = (sumPtR03 + emEtR03 + hadEtR03)/smallestMuPt;
468 
469  if(relIso<0.1) return true;
470  else return false;
471 }
472 
473 bool
474 PFMuonAlgo::isTightMuonPOG(const reco::MuonRef& muonRef) {
475 
476  if(!muon::isGoodMuon(*muonRef,muon::GlobalMuonPromptTight)) return false;
477 
478  if(!muonRef->isTrackerMuon()) return false;
479 
480  if(muonRef->numberOfMatches()<2) return false;
481 
482  //const reco::TrackRef& combinedMuon = muonRef->combinedMuon();
483  const reco::TrackRef& combinedMuon = muonRef->globalTrack();
484 
485  if(combinedMuon->hitPattern().numberOfValidTrackerHits()<11) return false;
486 
487  if(combinedMuon->hitPattern().numberOfValidPixelHits()==0) return false;
488 
489  if(combinedMuon->hitPattern().numberOfValidMuonHits()==0) return false;
490 
491  return true;
492 
493 }
494 
495 void
496 PFMuonAlgo::printMuonProperties(const reco::MuonRef& muonRef){
497 
498  if ( !muonRef.isNonnull() ) return;
499 
500  bool isGL = muonRef->isGlobalMuon();
501  bool isTR = muonRef->isTrackerMuon();
502  bool isST = muonRef->isStandAloneMuon();
503 
504  std::cout<<" GL: "<<isGL<<" TR: "<<isTR<<" ST: "<<isST<<std::endl;
505  std::cout<<" nMatches "<<muonRef->numberOfMatches()<<std::endl;
506 
507  if ( muonRef->isGlobalMuon() ){
508  reco::TrackRef combinedMu = muonRef->combinedMuon();
509  std::cout<<" GL, pt: " << combinedMu->pt()
510  << " +/- " << combinedMu->ptError()/combinedMu->pt()
511  << " chi**2 GBL : " << combinedMu->normalizedChi2()<<std::endl;
512  std::cout<< " Total Muon Hits : " << combinedMu->hitPattern().numberOfValidMuonHits()
513  << "/" << combinedMu->hitPattern().numberOfLostMuonHits()
514  << " DT Hits : " << combinedMu->hitPattern().numberOfValidMuonDTHits()
515  << "/" << combinedMu->hitPattern().numberOfLostMuonDTHits()
516  << " CSC Hits : " << combinedMu->hitPattern().numberOfValidMuonCSCHits()
517  << "/" << combinedMu->hitPattern().numberOfLostMuonCSCHits()
518  << " RPC Hits : " << combinedMu->hitPattern().numberOfValidMuonRPCHits()
519  << "/" << combinedMu->hitPattern().numberOfLostMuonRPCHits()<<std::endl;
520 
521  std::cout<<" # of Valid Tracker Hits "<<combinedMu->hitPattern().numberOfValidTrackerHits()<<std::endl;
522  std::cout<<" # of Valid Pixel Hits "<<combinedMu->hitPattern().numberOfValidPixelHits()<<std::endl;
523  }
524  if ( muonRef->isStandAloneMuon() ){
525  reco::TrackRef standAloneMu = muonRef->standAloneMuon();
526  std::cout<<" ST, pt: " << standAloneMu->pt()
527  << " +/- " << standAloneMu->ptError()/standAloneMu->pt()
528  << " eta : " << standAloneMu->eta()
529  << " DT Hits : " << standAloneMu->hitPattern().numberOfValidMuonDTHits()
530  << "/" << standAloneMu->hitPattern().numberOfLostMuonDTHits()
531  << " CSC Hits : " << standAloneMu->hitPattern().numberOfValidMuonCSCHits()
532  << "/" << standAloneMu->hitPattern().numberOfLostMuonCSCHits()
533  << " RPC Hits : " << standAloneMu->hitPattern().numberOfValidMuonRPCHits()
534  << "/" << standAloneMu->hitPattern().numberOfLostMuonRPCHits()
535  << " chi**2 STA : " << standAloneMu->normalizedChi2()<<std::endl;
536  }
537 
538 
539  if ( muonRef->isTrackerMuon() ){
540  reco::TrackRef trackerMu = muonRef->track();
541  const reco::Track& track = *trackerMu;
542  std::cout<<" TR, pt: " << trackerMu->pt()
543  << " +/- " << trackerMu->ptError()/trackerMu->pt()
544  << " chi**2 TR : " << trackerMu->normalizedChi2()<<std::endl;
545  std::cout<<" nTrackerHits "<<track.hitPattern().numberOfValidTrackerHits()<<std::endl;
546  std::cout<< "TMLastStationAngLoose "
547  << muon::isGoodMuon(*muonRef,muon::TMLastStationAngLoose) << std::endl
548  << "TMLastStationAngTight "
549  << muon::isGoodMuon(*muonRef,muon::TMLastStationAngTight) << std::endl
550  << "TMLastStationLoose "
551  << muon::isGoodMuon(*muonRef,muon::TMLastStationLoose) << std::endl
552  << "TMLastStationTight "
553  << muon::isGoodMuon(*muonRef,muon::TMLastStationTight) << std::endl
554  << "TMOneStationLoose "
555  << muon::isGoodMuon(*muonRef,muon::TMOneStationLoose) << std::endl
556  << "TMOneStationTight "
557  << muon::isGoodMuon(*muonRef,muon::TMOneStationTight) << std::endl
558  << "TMLastStationOptimizedLowPtLoose "
559  << muon::isGoodMuon(*muonRef,muon::TMLastStationOptimizedLowPtLoose) << std::endl
560  << "TMLastStationOptimizedLowPtTight "
561  << muon::isGoodMuon(*muonRef,muon::TMLastStationOptimizedLowPtTight) << std::endl
562  << "TMLastStationOptimizedBarrelLowPtLoose "
563  << muon::isGoodMuon(*muonRef,muon::TMLastStationOptimizedBarrelLowPtLoose) << std::endl
564  << "TMLastStationOptimizedBarrelLowPtTight "
565  << muon::isGoodMuon(*muonRef,muon::TMLastStationOptimizedBarrelLowPtTight) << std::endl
566  << std::endl;
567 
568  }
569 
570  std::cout<< "TM2DCompatibilityLoose "
571  << muon::isGoodMuon(*muonRef,muon::TM2DCompatibilityLoose) << std::endl
572  << "TM2DCompatibilityTight "
573  << muon::isGoodMuon(*muonRef,muon::TM2DCompatibilityTight) << std::endl;
574 
575 
576 
577  if ( muonRef->isGlobalMuon() && muonRef->isTrackerMuon() && muonRef->isStandAloneMuon() ){
578  reco::TrackRef combinedMu = muonRef->combinedMuon();
579  reco::TrackRef trackerMu = muonRef->track();
580  reco::TrackRef standAloneMu = muonRef->standAloneMuon();
581 
582  double sigmaCombined = combinedMu->ptError()/(combinedMu->pt()*combinedMu->pt());
583  double sigmaTracker = trackerMu->ptError()/(trackerMu->pt()*trackerMu->pt());
584  double sigmaStandAlone = standAloneMu->ptError()/(standAloneMu->pt()*standAloneMu->pt());
585 
586  bool combined = combinedMu->ptError()/combinedMu->pt() < 0.20;
587  bool tracker = trackerMu->ptError()/trackerMu->pt() < 0.20;
588  bool standAlone = standAloneMu->ptError()/standAloneMu->pt() < 0.20;
589 
590  double delta1 = combined && tracker ?
591  fabs(1./combinedMu->pt() -1./trackerMu->pt())
592  /sqrt(sigmaCombined*sigmaCombined + sigmaTracker*sigmaTracker) : 100.;
593  double delta2 = combined && standAlone ?
594  fabs(1./combinedMu->pt() -1./standAloneMu->pt())
595  /sqrt(sigmaCombined*sigmaCombined + sigmaStandAlone*sigmaStandAlone) : 100.;
596  double delta3 = standAlone && tracker ?
597  fabs(1./standAloneMu->pt() -1./trackerMu->pt())
598  /sqrt(sigmaStandAlone*sigmaStandAlone + sigmaTracker*sigmaTracker) : 100.;
599 
600  double delta =
601  standAloneMu->hitPattern().numberOfValidMuonDTHits()+
602  standAloneMu->hitPattern().numberOfValidMuonCSCHits() > 0 ?
603  std::min(delta3,std::min(delta1,delta2)) : std::max(delta3,std::max(delta1,delta2));
604 
605  std::cout << "delta = " << delta << " delta1 "<<delta1<<" delta2 "<<delta2<<" delta3 "<<delta3<<std::endl;
606 
607  double ratio =
608  combinedMu->ptError()/combinedMu->pt()
609  / (trackerMu->ptError()/trackerMu->pt());
610  //if ( ratio > 2. && delta < 3. ) std::cout << "ALARM ! " << ratio << ", " << delta << std::endl;
611  std::cout<<" ratio "<<ratio<<" combined mu pt "<<combinedMu->pt()<<std::endl;
612  //bool quality3 = ( combinedMu->pt() < 50. || ratio < 2. ) && delta < 3.;
613 
614 
615  }
616 
617  double sumPtR03 = muonRef->isolationR03().sumPt;
618  double emEtR03 = muonRef->isolationR03().emEt;
619  double hadEtR03 = muonRef->isolationR03().hadEt;
620  double relIsoR03 = (sumPtR03 + emEtR03 + hadEtR03)/muonRef->pt();
621  double sumPtR05 = muonRef->isolationR05().sumPt;
622  double emEtR05 = muonRef->isolationR05().emEt;
623  double hadEtR05 = muonRef->isolationR05().hadEt;
624  double relIsoR05 = (sumPtR05 + emEtR05 + hadEtR05)/muonRef->pt();
625  std::cout<<" 0.3 Radion Rel Iso: "<<relIsoR03<<" sumPt "<<sumPtR03<<" emEt "<<emEtR03<<" hadEt "<<hadEtR03<<std::endl;
626  std::cout<<" 0.5 Radion Rel Iso: "<<relIsoR05<<" sumPt "<<sumPtR05<<" emEt "<<emEtR05<<" hadEt "<<hadEtR05<<std::endl;
627  return;
628 
629 }
630 
631 
632 
633 std::vector<reco::Muon::MuonTrackTypePair> PFMuonAlgo::goodMuonTracks(const reco::MuonRef& muon,bool includeSA) {
634  return muonTracks(muon,includeSA,maxDPtOPt_);
635 }
636 
637 
638 
639 std::vector<reco::Muon::MuonTrackTypePair> PFMuonAlgo::muonTracks(const reco::MuonRef& muon,bool includeSA,double dpt) {
640 
641 
642 
643  std::vector<reco::Muon::MuonTrackTypePair> out;
644 
645 
646  if(muon->globalTrack().isNonnull() && muon->globalTrack()->pt()>0)
647  if(muon->globalTrack()->ptError()/muon->globalTrack()->pt()<dpt)
648  out.push_back(std::make_pair(muon->globalTrack(),reco::Muon::CombinedTrack));
649 
650  if(muon->innerTrack().isNonnull() && muon->innerTrack()->pt()>0)
651  if(muon->innerTrack()->ptError()/muon->innerTrack()->pt()<dpt)//Here Loose!@
652  out.push_back(std::make_pair(muon->innerTrack(),reco::Muon::InnerTrack));
653 
654  bool pickyExists=false;
655  double pickyDpt=99999.;
656  if(muon->pickyTrack().isNonnull() && muon->pickyTrack()->pt()>0) {
657  pickyDpt = muon->pickyTrack()->ptError()/muon->pickyTrack()->pt();
658  if(pickyDpt<dpt)
659  out.push_back(std::make_pair(muon->pickyTrack(),reco::Muon::Picky));
660  pickyExists=true;
661  }
662 
663  bool dytExists=false;
664  double dytDpt=99999.;
665  if(muon->dytTrack().isNonnull() && muon->dytTrack()->pt()>0) {
666  dytDpt = muon->dytTrack()->ptError()/muon->dytTrack()->pt();
667  if(dytDpt<dpt)
668  out.push_back(std::make_pair(muon->dytTrack(),reco::Muon::DYT));
669  dytExists=true;
670  }
671 
672  //Magic: TPFMS is not a really good track especially under misalignment
673  //IT is kind of crap because if mu system is displaced it can make a change
674  //So allow TPFMS if there is no picky or the error of tpfms is better than picky
675  //AND if there is no DYT or the error of tpfms is better than DYT
676  if(muon->tpfmsTrack().isNonnull() && muon->tpfmsTrack()->pt()>0) {
677  double tpfmsDpt = muon->tpfmsTrack()->ptError()/muon->tpfmsTrack()->pt();
678  if( ( (pickyExists && tpfmsDpt<pickyDpt) || (!pickyExists) ) &&
679  ( (dytExists && tpfmsDpt<dytDpt) || (!dytExists) ) &&
680  tpfmsDpt<dpt )
681  out.push_back(std::make_pair(muon->tpfmsTrack(),reco::Muon::TPFMS));
682  }
683 
684  if(includeSA && muon->outerTrack().isNonnull())
685  if(muon->outerTrack()->ptError()/muon->outerTrack()->pt()<dpt)
686  out.push_back(std::make_pair(muon->outerTrack(),reco::Muon::OuterTrack));
687 
688  return out;
689 
690 }
691 
692 
693 
694 
695 
697 
698 
699 
700 
701 bool PFMuonAlgo::reconstructMuon(reco::PFCandidate& candidate, const reco::MuonRef& muon, bool allowLoose) {
702  using namespace std;
703  using namespace reco;
704 
705  if (!muon.isNonnull())
706  return false;
707 
708 
709 
710  bool isMu=false;
711 
712  if(allowLoose)
713  isMu = isMuon(muon) || isLooseMuon(muon);
714  else
715  isMu = isMuon(muon);
716 
717  if( !isMu)
718  return false;
719 
720  //get the valid tracks(without standalone except we allow loose muons)
721  //MIKE: Here we need to be careful. If we have a muon inside a dense
722  //jet environment often the track is badly measured. In this case
723  //we should not apply Dpt/Pt<1
724 
725  std::vector<reco::Muon::MuonTrackTypePair> validTracks = goodMuonTracks(muon);
726  if (!allowLoose)
727  validTracks = goodMuonTracks(muon);
728  else
729  validTracks = muonTracks(muon);
730 
731  if( validTracks.size() ==0)
732  return false;
733 
734 
735 
736  //check what is the track used.Rerun TuneP
737  reco::Muon::MuonTrackTypePair bestTrackPair = muon::tevOptimized(*muon);
738 
739  TrackRef bestTrack = bestTrackPair.first;
740  MuonTrackType trackType = bestTrackPair.second;
741 
742 
743 
744  MuonTrackTypePair trackPairWithSmallestError = getTrackWithSmallestError(validTracks);
745  TrackRef trackWithSmallestError = trackPairWithSmallestError.first;
746 
747  if( trackType == reco::Muon::InnerTrack &&
748  (!bestTrack->quality(trackQuality_) ||
749  bestTrack->ptError()/bestTrack->pt()> errorCompScale_*trackWithSmallestError->ptError()/trackWithSmallestError->pt() )) {
750  bestTrack = trackWithSmallestError;
751  trackType = trackPairWithSmallestError.second;
752  }
753  else if (trackType != reco::Muon::InnerTrack &&
754  bestTrack->ptError()/bestTrack->pt()> errorCompScale_*trackWithSmallestError->ptError()/trackWithSmallestError->pt()) {
755  bestTrack = trackWithSmallestError;
756  trackType = trackPairWithSmallestError.second;
757 
758  }
759 
760 
761  changeTrack(candidate,std::make_pair(bestTrack,trackType));
762  candidate.setMuonRef( muon );
763 
764  return true;
765 }
766 
767 
768 
769 
770 void PFMuonAlgo::changeTrack(reco::PFCandidate& candidate,const MuonTrackTypePair& track) {
771  using namespace reco;
772  reco::TrackRef bestTrack = track.first;
773  MuonTrackType trackType = track.second;
774  //OK Now redefine the canddiate with that track
775  double px = bestTrack->px();
776  double py = bestTrack->py();
777  double pz = bestTrack->pz();
778  double energy = sqrt(bestTrack->p()*bestTrack->p() + 0.1057*0.1057);
779 
780  candidate.setCharge(bestTrack->charge()>0 ? 1 : -1);
781  candidate.setP4(math::XYZTLorentzVector(px,py,pz,energy));
782  candidate.setParticleType(reco::PFCandidate::mu);
783  // candidate.setTrackRef( bestTrack );
784  candidate.setMuonTrackType(trackType);
785  if(trackType == reco::Muon::InnerTrack)
786  candidate.setVertexSource( PFCandidate::kTrkMuonVertex );
787  else if(trackType == reco::Muon::CombinedTrack)
788  candidate.setVertexSource( PFCandidate::kComMuonVertex );
789  else if(trackType == reco::Muon::TPFMS)
790  candidate.setVertexSource( PFCandidate::kTPFMSMuonVertex );
791  else if(trackType == reco::Muon::Picky)
792  candidate.setVertexSource( PFCandidate::kPickyMuonVertex );
793  else if(trackType == reco::Muon::DYT)
794  candidate.setVertexSource( PFCandidate::kDYTMuonVertex );
795  }
796 
797 
798 reco::Muon::MuonTrackTypePair
799 PFMuonAlgo::getTrackWithSmallestError(const std::vector<reco::Muon::MuonTrackTypePair>& tracks) {
800  TrackPtErrorSorter sorter;
801  return *std::min_element(tracks.begin(),tracks.end(),sorter);
802 }
803 
804 
805 
806 
807 
808 void PFMuonAlgo::estimateEventQuantities(const reco::PFCandidateCollection* pfc)
809 {
810  //SUM ET
811  sumetPU_ = 0.0;
812  METX_=0.;
813  METY_=0.;
814  sumet_=0.0;
815  for(reco::PFCandidateCollection::const_iterator i = pfc->begin();i!=pfc->end();++i) {
816  sumet_+=i->pt();
817  METX_+=i->px();
818  METY_+=i->py();
819  }
820 
821 }
822 
823 
824 
825 
826 
827 void PFMuonAlgo::postClean(reco::PFCandidateCollection* cands) {
828  using namespace std;
829  using namespace reco;
830  if (!postCleaning_)
831  return;
832 
833  //Initialize vectors
834 
835  if(pfCosmicsMuonCleanedCandidates_.get() )
836  pfCosmicsMuonCleanedCandidates_->clear();
837  else
838  pfCosmicsMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );
839 
840  if(pfCleanedTrackerAndGlobalMuonCandidates_.get() )
841  pfCleanedTrackerAndGlobalMuonCandidates_->clear();
842  else
843  pfCleanedTrackerAndGlobalMuonCandidates_.reset( new reco::PFCandidateCollection );
844 
845  if( pfFakeMuonCleanedCandidates_.get() )
846  pfFakeMuonCleanedCandidates_->clear();
847  else
848  pfFakeMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );
849 
850 
851  if( pfPunchThroughMuonCleanedCandidates_.get() )
852  pfPunchThroughMuonCleanedCandidates_->clear();
853  else
854  pfPunchThroughMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );
855 
856  if( pfPunchThroughHadronCleanedCandidates_.get() )
857  pfPunchThroughHadronCleanedCandidates_->clear();
858  else
859  pfPunchThroughHadronCleanedCandidates_.reset( new reco::PFCandidateCollection );
860 
861 
862 
863  pfPunchThroughHadronCleanedCandidates_->clear();
864 
865  maskedIndices_.clear();
866 
867  //Estimate MET and SumET
868 
869  estimateEventQuantities(cands);
870 
871 
872  std::vector<int> muons;
873  std::vector<int> cosmics;
874  //get the muons
875  for(unsigned int i=0;i<cands->size();++i) {
876  const reco::PFCandidate& cand = (*cands)[i];
877  if ( cand.particleId() == reco::PFCandidate::mu )
878  muons.push_back(i);
879  else if ( cand.particleId() == reco::PFCandidate::h && cand.muonRef().isNonnull()) {
880  //MET cleaning for muons that are not high purity and became charged hadrons
881  if (cand.pt()>100.0 && (!cand.trackRef()->quality(trackQuality_)) && cand.muonRef()->isGlobalMuon() && cand.rawHcalEnergy()<0.05*cand.p()) {
882  maskedIndices_.push_back(i);
883  pfPunchThroughMuonCleanedCandidates_->push_back(cand);
884  }
885  }
886 
887  }
888  //Then sort the muon indicess by decsending pt
889 
890  IndexPtComparator comparator(cands);
891  std::sort(muons.begin(),muons.end(),comparator);
892 
893 
894  //first kill cosmics
895  double METXCosmics=0;
896  double METYCosmics=0;
897  double SUMETCosmics=0.0;
898 
899  for(unsigned int i=0;i<muons.size();++i) {
900  const PFCandidate& pfc = cands->at(muons[i]);
901  double origin=0.0;
902  if(vertices_->size()>0&& vertices_->at(0).isValid() && ! vertices_->at(0).isFake())
903  origin = pfc.muonRef()->muonBestTrack()->dxy(vertices_->at(0).position());
904 
905  if( origin> cosmicRejDistance_) {
906  cosmics.push_back(muons[i]);
907  METXCosmics +=pfc.px();
908  METYCosmics +=pfc.py();
909  SUMETCosmics +=pfc.pt();
910  }
911  }
912  double MET2Cosmics = METXCosmics*METXCosmics+METYCosmics*METYCosmics;
913 
914  if ( SUMETCosmics > (sumet_-sumetPU_)/eventFactorCosmics_ && MET2Cosmics < METX_*METX_+ METY_*METY_)
915  for(unsigned int i=0;i<cosmics.size();++i) {
916  maskedIndices_.push_back(cosmics[i]);
917  pfCosmicsMuonCleanedCandidates_->push_back(cands->at(cosmics[i]));
918  }
919 
920  //Loop on the muons candidates and clean
921  for(unsigned int i=0;i<muons.size();++i) {
922  if( cleanMismeasured(cands->at(muons[i]),muons[i]))
923  continue;
924  cleanPunchThroughAndFakes(cands->at(muons[i]),cands,muons[i]);
925 
926  }
927 
928 
929 
930  //OK Now do the hard job ->remove the candidates that were cleaned
931  removeDeadCandidates(cands,maskedIndices_);
932 
933 
934 
935 
936 }
937 
938 void PFMuonAlgo::addMissingMuons(edm::Handle<reco::MuonCollection> muons, reco::PFCandidateCollection* cands) {
939  if(!postCleaning_)
940  return;
941 
942 
943  if( pfAddedMuonCandidates_.get() )
944  pfAddedMuonCandidates_->clear();
945  else
946  pfAddedMuonCandidates_.reset( new reco::PFCandidateCollection );
947 
948 
949 
950  for ( unsigned imu = 0; imu < muons->size(); ++imu ) {
951  reco::MuonRef muonRef( muons, imu );
952  bool used = false;
953  bool hadron = false;
954  for(unsigned i=0; i<cands->size(); i++) {
955  const PFCandidate& pfc = cands->at(i);
956  if ( !pfc.trackRef().isNonnull() ) continue;
957  if ( pfc.trackRef().isNonnull() && pfc.trackRef() == muonRef->track() )
958  hadron = true;
959  if ( !pfc.muonRef().isNonnull() ) continue;
960 
961  if ( pfc.muonRef()->innerTrack() == muonRef->innerTrack())
962  used = true;
963  else {
964  // Check if the stand-alone muon is not a spurious copy of an existing muon
965  // (Protection needed for HLT)
966  if ( pfc.muonRef()->isStandAloneMuon() && muonRef->isStandAloneMuon() ) {
967  double dEta = pfc.muonRef()->standAloneMuon()->eta() - muonRef->standAloneMuon()->eta();
968  double dPhi = pfc.muonRef()->standAloneMuon()->phi() - muonRef->standAloneMuon()->phi();
969  double dR = sqrt(dEta*dEta + dPhi*dPhi);
970  if ( dR < 0.005 ) {
971  used = true;
972  }
973  }
974  }
975 
976  if ( used ) break;
977  }
978 
979  if ( used ||hadron||(!muonRef.isNonnull()) ) continue;
980 
981 
982  TrackMETComparator comparator(METX_,METY_);
983  //Low pt dont need to be cleaned
984 
985  std::vector<reco::Muon::MuonTrackTypePair> tracks = goodMuonTracks(muonRef,true);
986  //If there is at least 1 track choice try to change the track
987  if(tracks.size()>0) {
988 
989  //Find tracks that change dramatically MET or Pt
990  std::vector<reco::Muon::MuonTrackTypePair> tracksThatChangeMET = tracksPointingAtMET(tracks);
991  //From those tracks get the one with smallest MET
992  if (tracksThatChangeMET.size()>0) {
993  reco::Muon::MuonTrackTypePair bestTrackType = *std::min_element(tracksThatChangeMET.begin(),tracksThatChangeMET.end(),comparator);
994 
995  //Make sure it is not cosmic
996  if((vertices_->size()==0) ||bestTrackType.first->dz(vertices_->at(0).position())<cosmicRejDistance_){
997 
998  //make a pfcandidate
999  int charge = bestTrackType.first->charge()>0 ? 1 : -1;
1000  math::XYZTLorentzVector momentum(bestTrackType.first->px(),
1001  bestTrackType.first->py(),
1002  bestTrackType.first->pz(),
1003  sqrt(bestTrackType.first->p()*bestTrackType.first->p()+0.1057*0.1057));
1004 
1005  cands->push_back( PFCandidate( charge,
1006  momentum,
1007  reco::PFCandidate::mu ) );
1008 
1009  changeTrack(cands->back(),bestTrackType);
1010 
1011  if (muonRef->track().isNonnull() )
1012  cands->back().setTrackRef( muonRef->track() );
1013 
1014  cands->back().setMuonRef(muonRef);
1015 
1016 
1017  pfAddedMuonCandidates_->push_back(cands->back());
1018 
1019  }
1020  }
1021  }
1022  }
1023 }
1024 
1025 std::pair<double,double>
1026 PFMuonAlgo::getMinMaxMET2(const reco::PFCandidate&pfc) {
1027  std::vector<reco::Muon::MuonTrackTypePair> tracks = goodMuonTracks((pfc.muonRef()),true);
1028 
1029  double METXNO = METX_-pfc.px();
1030  double METYNO = METY_-pfc.py();
1031  std::vector<double> met2;
1032  for (unsigned int i=0;i<tracks.size();++i) {
1033  met2.push_back(pow(METXNO+tracks.at(i).first->px(),2)+pow(METYNO+tracks.at(i).first->py(),2));
1034  }
1035 
1036  //PROTECT for cases of only one track. If there is only one track it will crash .
1037  //Has never happened but could likely happen!
1038 
1039  if(tracks.size()>1)
1040  return std::make_pair(*std::min_element(met2.begin(),met2.end()),*std::max_element(met2.begin(),met2.end()));
1041  else
1042  return std::make_pair(0,0);
1043 }
1044 
1045 
1046 bool PFMuonAlgo::cleanMismeasured(reco::PFCandidate& pfc,unsigned int i ){
1047  using namespace std;
1048  using namespace reco;
1049  bool cleaned=false;
1050 
1051  //First define the MET without this guy
1052  double METNOX = METX_ - pfc.px();
1053  double METNOY = METY_ - pfc.py();
1054  double SUMETNO = sumet_ -pfc.pt();
1055 
1056  TrackMETComparator comparator(METNOX,METNOY);
1057  //Low pt dont need to be cleaned
1058  if (pfc.pt()<minPostCleaningPt_)
1059  return false;
1060  std::vector<reco::Muon::MuonTrackTypePair> tracks = goodMuonTracks(pfc.muonRef(),false);
1061 
1062 
1063 
1064  //If there is more than 1 track choice try to change the track
1065  if(tracks.size()>1) {
1066  //Find tracks that change dramatically MET or Pt
1067  std::vector<reco::Muon::MuonTrackTypePair> tracksThatChangeMET = tracksWithBetterMET(tracks,pfc);
1068  //From those tracks get the one with smallest MET
1069  if (tracksThatChangeMET.size()>0) {
1070  reco::Muon::MuonTrackTypePair bestTrackType = *std::min_element(tracksThatChangeMET.begin(),tracksThatChangeMET.end(),comparator);
1071  changeTrack(pfc,bestTrackType);
1072 
1073  pfCleanedTrackerAndGlobalMuonCandidates_->push_back(pfc);
1074  //update eventquantities
1075  METX_ = METNOX+pfc.px();
1076  METY_ = METNOY+pfc.py();
1077  sumet_=SUMETNO+pfc.pt();
1078 
1079  }
1080  }
1081 
1082  //Now attempt to kill it
1083  if (!(pfc.muonRef()->isGlobalMuon() && pfc.muonRef()->isTrackerMuon())) {
1084  //define MET significance and SUM ET
1085  double MET2 = METX_*METX_+METY_*METY_;
1086  double newMET2 = METNOX*METNOX+METNOY*METNOY;
1087  double METSig = sqrt(MET2)/sqrt(sumet_-sumetPU_);
1088  if( METSig>metSigForRejection_)
1089  if((newMET2 < MET2/metFactorRejection_) &&
1090  ((SUMETNO-sumetPU_)/(sumet_-sumetPU_)<eventFractionRejection_)) {
1091  pfFakeMuonCleanedCandidates_->push_back(pfc);
1092  maskedIndices_.push_back(i);
1093  METX_ = METNOX;
1094  METY_ = METNOY;
1095  sumet_=SUMETNO;
1096  cleaned=true;
1097  }
1098 
1099  }
1100  return cleaned;
1101 
1102 }
1103 
1104 std::vector<reco::Muon::MuonTrackTypePair>
1105 PFMuonAlgo::tracksWithBetterMET(const std::vector<reco::Muon::MuonTrackTypePair>& tracks ,const reco::PFCandidate& pfc) {
1106  std::vector<reco::Muon::MuonTrackTypePair> outputTracks;
1107 
1108  double METNOX = METX_ - pfc.px();
1109  double METNOY = METY_ - pfc.py();
1110  double SUMETNO = sumet_ -pfc.pt();
1111  double MET2 = METX_*METX_+METY_*METY_;
1112  double newMET2=0.0;
1113  double newSUMET=0.0;
1114  double METSIG = sqrt(MET2)/sqrt(sumet_-sumetPU_);
1115 
1116 
1117  if(METSIG>metSigForCleaning_)
1118  for( unsigned int i=0;i<tracks.size();++i) {
1119  //calculate new SUM ET and MET2
1120  newSUMET = SUMETNO+tracks.at(i).first->pt()-sumetPU_;
1121  newMET2 = pow(METNOX+tracks.at(i).first->px(),2)+pow(METNOY+tracks.at(i).first->py(),2);
1122 
1123  if(newSUMET/(sumet_-sumetPU_)>eventFractionCleaning_ && newMET2<MET2/metFactorCleaning_)
1124  outputTracks.push_back(tracks.at(i));
1125  }
1126 
1127 
1128  return outputTracks;
1129 }
1130 
1131 
1132 std::vector<reco::Muon::MuonTrackTypePair>
1133 PFMuonAlgo::tracksPointingAtMET(const std::vector<reco::Muon::MuonTrackTypePair>& tracks) {
1134  std::vector<reco::Muon::MuonTrackTypePair> outputTracks;
1135 
1136 
1137  double newMET2=0.0;
1138 
1139  for( unsigned int i=0;i<tracks.size();++i) {
1140  //calculate new SUM ET and MET2
1141  newMET2 = pow(METX_+tracks.at(i).first->px(),2)+pow(METY_+tracks.at(i).first->py(),2);
1142 
1143  if(newMET2<(METX_*METX_+METY_*METY_)/metFactorCleaning_)
1144  outputTracks.push_back(tracks.at(i));
1145  }
1146 
1147 
1148  return outputTracks;
1149 }
1150 
1151 
1152 
1153 void PFMuonAlgo::setInputsForCleaning(const reco::VertexCollection* vertices) {
1154  vertices_ = vertices;
1155 }
1156 
1157 bool PFMuonAlgo::cleanPunchThroughAndFakes(reco::PFCandidate&pfc,reco::PFCandidateCollection* cands,unsigned int imu ){
1158  using namespace reco;
1159 
1160  bool cleaned=false;
1161 
1162  if (pfc.pt()<minPostCleaningPt_)
1163  return false;
1164 
1165 
1166  double METXNO = METX_-pfc.pt();
1167  double METYNO = METY_-pfc.pt();
1168  double MET2NO = METXNO*METXNO+METYNO*METYNO;
1169  double MET2 = METX_*METX_+METY_*METY_;
1170  bool fake1=false;
1171 
1172  std::pair<double,double> met2 = getMinMaxMET2(pfc);
1173 
1174  //Check for Fakes at high pseudorapidity
1175  if(pfc.muonRef()->standAloneMuon().isNonnull())
1176  fake1 =fabs ( pfc.eta() ) > 2.15 &&
1177  met2.first<met2.second/2 &&
1178  MET2NO < MET2/metFactorHighEta_ &&
1179  pfc.muonRef()->standAloneMuon()->pt() < pfc.pt()/ptFactorHighEta_;
1180 
1181  double factor = std::max(2.,2000./(sumet_-pfc.pt()-sumetPU_));
1182  bool fake2 = ( pfc.pt()/(sumet_-sumetPU_) < 0.25 && MET2NO < MET2/metFactorFake_ && met2.first<met2.second/factor );
1183 
1184  bool punchthrough =pfc.p() > minPunchThroughMomentum_ &&
1185  pfc.rawHcalEnergy() > minPunchThroughEnergy_ &&
1186  pfc.rawEcalEnergy()+pfc.rawHcalEnergy() > pfc.p()/punchThroughFactor_ &&
1187  !isIsolatedMuon(pfc.muonRef()) && MET2NO < MET2/punchThroughMETFactor_;
1188 
1189 
1190  if(fake1 || fake2||punchthrough) {
1191  // Find the block of the muon
1192  const PFCandidate::ElementsInBlocks& eleInBlocks = pfc.elementsInBlocks();
1193  if ( eleInBlocks.size() ) {
1194  PFBlockRef blockRefMuon = eleInBlocks[0].first;
1195  unsigned indexMuon = eleInBlocks[0].second;
1196  if (eleInBlocks.size()>1)
1197  indexMuon = eleInBlocks[1].second;
1198 
1199  // Check if the muon gave rise to a neutral hadron
1200  double iHad = 1E9;
1201  bool hadron = false;
1202  for ( unsigned i = imu+1; i < cands->size(); ++i ) {
1203  const PFCandidate& pfcn = cands->at(i);
1204  const PFCandidate::ElementsInBlocks& ele = pfcn.elementsInBlocks();
1205  if ( !ele.size() ) {
1206  continue;
1207  }
1208  PFBlockRef blockRefHadron = ele[0].first;
1209  unsigned indexHadron = ele[0].second;
1210  // We are out of the block -> exit the loop
1211  if ( blockRefHadron.key() != blockRefMuon.key() ) break;
1212  // Check that this particle is a neutral hadron
1213  if ( indexHadron == indexMuon &&
1214  pfcn.particleId() == reco::PFCandidate::h0 ) {
1215  iHad = i;
1216  hadron = true;
1217  }
1218  if ( hadron ) break;
1219  }
1220 
1221  if ( hadron ) {
1222 
1223  double rescaleFactor = cands->at(iHad).p()/cands->at(imu).p();
1224  METX_ -= cands->at(imu).px() + cands->at(iHad).px();
1225  METY_ -= cands->at(imu).py() + cands->at(iHad).py();
1226  sumet_ -=cands->at(imu).pt();
1227  cands->at(imu).rescaleMomentum(rescaleFactor);
1228  maskedIndices_.push_back(iHad);
1229  pfPunchThroughHadronCleanedCandidates_->push_back(cands->at(iHad));
1230  cands->at(imu).setParticleType(reco::PFCandidate::h);
1231  pfPunchThroughMuonCleanedCandidates_->push_back(cands->at(imu));
1232  METX_ += cands->at(imu).px();
1233  METY_ += cands->at(imu).py();
1234  sumet_ += cands->at(imu).pt();
1235 
1236  } else if ( fake1 || fake2 ) {
1237  METX_ -= cands->at(imu).px();
1238  METY_ -= cands->at(imu).py();
1239  sumet_ -= cands->at(imu).pt();
1240  maskedIndices_.push_back(imu);
1241  pfFakeMuonCleanedCandidates_->push_back(cands->at(imu));
1242  cleaned=true;
1243  }
1244  }
1245  }
1246  return cleaned;
1247 }
1248 
1249 
1250 void PFMuonAlgo::removeDeadCandidates(reco::PFCandidateCollection* obj, const std::vector<unsigned int>& indices)
1251 {
1252  size_t N = indices.size();
1253  size_t collSize = obj->size();
1254 
1255  for (size_t i = 0 ; i < N ; ++i)
1256  obj->at(indices.at(i)) = obj->at(collSize-i-1);
1257 
1258  obj->resize(collSize - indices.size());
1259 }
delta
dbl * delta
Definition: mlp_gen.cc:36
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
reco::PFBlockElement
Abstract base class for a PFBlock element (track, cluster...)
Definition: PFBlockElement.h:28
PFMuonAlgo::isIsolatedMuon
static bool isIsolatedMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:224
i
int i
Definition: DBlmapReader.cc:9
PFMuonAlgo::cleanMismeasured
bool cleanMismeasured(reco::PFCandidate &, unsigned int)
Definition: PFMuonAlgo.cc:1046
reco::LeafCandidate::p
virtual double p() const
magnitude of momentum vector
Definition: LeafCandidate.h:157
PFMuonAlgo::isTightMuonPOG
static bool isTightMuonPOG(const reco::MuonRef &muonRef)
Definition: PFMuonAlgo.cc:474
PFMuonAlgo::removeDeadCandidates
void removeDeadCandidates(reco::PFCandidateCollection *, const std::vector< unsigned int > &)
Definition: PFMuonAlgo.cc:1250
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
reco::PFCandidate::rawEcalEnergy
double rawEcalEnergy() const
return corrected Ecal energy
Definition: PFCandidate.h:220
PFMuonAlgo::addMissingMuons
void addMissingMuons(edm::Handle< reco::MuonCollection >, reco::PFCandidateCollection *cands)
Definition: PFMuonAlgo.cc:938
reco::LeafCandidate::pt
virtual float pt() const
transverse momentum
Definition: LeafCandidate.h:178
reco::isMuon
bool isMuon(const Candidate &part)
Definition: pdgIdUtils.h:11
PFMuonAlgo::isMuon
static bool isMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:153
PFMuonAlgo::isTrackerLooseMuon
static bool isTrackerLooseMuon(const reco::PFBlockElement &elt)
reco::LeafCandidate::setCharge
virtual void setCharge(Charge q)
set electric charge
Definition: LeafCandidate.h:142
PFMuonAlgo::setParameters
void setParameters(const edm::ParameterSet &)
Definition: PFMuonAlgo.cc:29
PFMuonAlgo::isGlobalLooseMuon
static bool isGlobalLooseMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:198
edm::ParameterSet::exists
bool exists(std::string const &parameterName) const
checks if a parameter exists
Definition: ParameterSet.cc:758
PFMuonAlgo::cleanPunchThroughAndFakes
bool cleanPunchThroughAndFakes(reco::PFCandidate &, reco::PFCandidateCollection *, unsigned int)
Definition: PFMuonAlgo.cc:1157
reco::LeafCandidate::setP4
virtual void setP4(const LorentzVector &p4)
set 4-momentum
Definition: LeafCandidate.h:190
reco::VertexCollection
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
PFMuonAlgo::MuonTrackTypePair
reco::Muon::MuonTrackTypePair MuonTrackTypePair
Definition: PFMuonAlgo.h:14
edm::Ref::key
key_type key() const
Accessor for product key.
Definition: Ref.h:266
PFMuonAlgo::printMuonProperties
static void printMuonProperties(const reco::MuonRef &muonRef)
Definition: PFMuonAlgo.cc:496
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:107
PFMuonAlgo::tracksPointingAtMET
std::vector< MuonTrackTypePair > tracksPointingAtMET(const std::vector< MuonTrackTypePair > &)
Definition: PFMuonAlgo.cc:1133
PFMuonAlgo::isGlobalTightMuon
static bool isGlobalTightMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:185
reco::PFCandidate::ElementsInBlocks
std::vector< ElementInBlock > ElementsInBlocks
Definition: PFCandidate.h:382
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
muon::isLooseMuon
bool isLooseMuon(const reco::Muon &)
Definition: MuonSelectors.cc:762
trackingTruthProducerFastSim_cfi.muon
tuple muon
Definition: trackingTruthProducerFastSim_cfi.py:10
reco::PFCandidate::trackRef
reco::TrackRef trackRef() const
Definition: PFCandidate.cc:433
PFMuonAlgo::getTrackWithSmallestError
MuonTrackTypePair getTrackWithSmallestError(const std::vector< MuonTrackTypePair > &)
Definition: PFMuonAlgo.cc:799
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:248
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
reco::PFCandidate::setVertexSource
void setVertexSource(PFVertexType vt)
Definition: PFCandidate.h:402
muon::tevOptimized
reco::Muon::MuonTrackTypePair tevOptimized(const reco::TrackRef &combinedTrack, const reco::TrackRef &trackerTrack, const reco::TrackRef &tpfmsTrack, const reco::TrackRef &pickyTrack, const reco::TrackRef &dytTrack, const double ptThreshold=200., const double tune1=17., const double tune2=40., const double dptcut=0.25)
Definition: MuonCocktails.cc:9
dPhi
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
reco::PFCandidate::setMuonTrackType
void setMuonTrackType(const reco::Muon::MuonTrackType &type)
set the Best Muon Track Ref
Definition: PFCandidate.h:352
reco::LeafCandidate::eta
virtual float eta() const
momentum pseudorapidity
Definition: LeafCandidate.h:184
reco::PFCandidate::setParticleType
void setParticleType(ParticleType type)
set Particle Type
Definition: PFCandidate.cc:256
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
query.result
tuple result
Definition: query.py:137
PFMuonAlgo::goodMuonTracks
std::vector< reco::Muon::MuonTrackTypePair > goodMuonTracks(const reco::MuonRef &muon, bool includeSA=false)
Definition: PFMuonAlgo.cc:633
patCandidatesForDimuonsSequences_cff.tracker
tuple tracker
Definition: patCandidatesForDimuonsSequences_cff.py:69
min
T min(T a, T b)
Definition: MathUtil.h:58
PFMuonAlgo::estimateEventQuantities
void estimateEventQuantities(const reco::PFCandidateCollection *)
Definition: PFMuonAlgo.cc:808
PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.dR
tuple dR
Definition: PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py:19
PFMuonAlgo::changeTrack
void changeTrack(reco::PFCandidate &, const MuonTrackTypePair &)
Definition: PFMuonAlgo.cc:770
muon::isGoodMuon
bool isGoodMuon(const reco::Muon &muon, SelectionType type, reco::Muon::ArbitrationType arbitrationType=reco::Muon::SegmentAndTrackArbitration)
main GoodMuon wrapper call
Definition: MuonSelectors.cc:567
reco::PFCandidate::muonRef
reco::MuonRef muonRef() const
Definition: PFCandidate.cc:450
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
reco::PFCandidateCollection
std::vector< reco::PFCandidate > PFCandidateCollection
collection of PFCandidates
Definition: PFCandidateFwd.h:12
reco::TrackBase::qualityByName
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:108
N
#define N
Definition: blowfish.cc:9
reco::LeafCandidate::px
virtual double px() const
x coordinate of momentum vector
Definition: LeafCandidate.h:172
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:384
testEve_cfg.tracks
tuple tracks
Definition: testEve_cfg.py:39
PFMuonAlgo::isLooseMuon
static bool isLooseMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:166
PFMuonAlgo::getMinMaxMET2
std::pair< double, double > getMinMaxMET2(const reco::PFCandidate &)
Definition: PFMuonAlgo.cc:1026
MCScenario_CRAFT1_22X.sorter
def sorter
Definition: MCScenario_CRAFT1_22X.py:96
reco::HitPattern::numberOfValidTrackerHits
int numberOfValidTrackerHits() const
Definition: HitPattern.h:739
reco::PFCandidate
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:39
reco::PFCandidate::setMuonRef
void setMuonRef(const reco::MuonRef &ref)
set muon reference
Definition: PFCandidate.cc:436
patZpeak.muons
tuple muons
Definition: patZpeak.py:38
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
PFMuonAlgo::PFMuonAlgo
PFMuonAlgo()
constructor
Definition: PFMuonAlgo.cc:17
reco::Muon::MuonTrackType
MuonTrackType
map for Global Muon refitters
Definition: Muon.h:38
reco::PFBlockElementTrack::muonRef
const reco::MuonRef & muonRef() const
Definition: PFBlockElementTrack.h:95
reco::PFCandidate::particleId
virtual ParticleType particleId() const
Definition: PFCandidate.h:369
reco::Muon::MuonTrackTypePair
std::pair< TrackRef, Muon::MuonTrackType > MuonTrackTypePair
Definition: Muon.h:40
PFMuonAlgo::tracksWithBetterMET
std::vector< MuonTrackTypePair > tracksWithBetterMET(const std::vector< MuonTrackTypePair > &, const reco::PFCandidate &)
Definition: PFMuonAlgo.cc:1105
PFMuonAlgo::postClean
void postClean(reco::PFCandidateCollection *)
Definition: PFMuonAlgo.cc:827
reco::PFCandidate::elementsInBlocks
const ElementsInBlocks & elementsInBlocks() const
Definition: PFCandidate.cc:682
PFMuonAlgo::setInputsForCleaning
void setInputsForCleaning(const reco::VertexCollection *)
Definition: PFMuonAlgo.cc:1153
PFMuonAlgo::muonTracks
std::vector< reco::Muon::MuonTrackTypePair > muonTracks(const reco::MuonRef &muon, bool includeSA=false, double dpt=1e+9)
Definition: PFMuonAlgo.cc:639
reco::LeafCandidate::py
virtual double py() const
y coordinate of momentum vector
Definition: LeafCandidate.h:174
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
PFMuonAlgo::isTrackerTightMuon
static bool isTrackerTightMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:211
reco::PFCandidate::rawHcalEnergy
double rawHcalEnergy() const
return raw Hcal energy
Definition: PFCandidate.h:230
PFMuonAlgo::reconstructMuon
bool reconstructMuon(reco::PFCandidate &, const reco::MuonRef &, bool allowLoose=false)
Definition: PFMuonAlgo.cc:701
Generated for CMSSW Reference Manual by doxygen 1.8.5