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 
465  double relIso = sumPtR03/smallestMuPt;
466 
467  if(relIso<0.1) return true;
468  else return false;
469 }
470 
471 bool
472 PFMuonAlgo::isTightMuonPOG(const reco::MuonRef& muonRef) {
473 
474  if(!muon::isGoodMuon(*muonRef,muon::GlobalMuonPromptTight)) return false;
475 
476  if(!muonRef->isTrackerMuon()) return false;
477 
478  if(muonRef->numberOfMatches()<2) return false;
479 
480  //const reco::TrackRef& combinedMuon = muonRef->combinedMuon();
481  const reco::TrackRef& combinedMuon = muonRef->globalTrack();
482 
483  if(combinedMuon->hitPattern().numberOfValidTrackerHits()<11) return false;
484 
485  if(combinedMuon->hitPattern().numberOfValidPixelHits()==0) return false;
486 
487  if(combinedMuon->hitPattern().numberOfValidMuonHits()==0) return false;
488 
489  return true;
490 
491 }
492 
493 void
494 PFMuonAlgo::printMuonProperties(const reco::MuonRef& muonRef){
495 
496  if ( !muonRef.isNonnull() ) return;
497 
498  bool isGL = muonRef->isGlobalMuon();
499  bool isTR = muonRef->isTrackerMuon();
500  bool isST = muonRef->isStandAloneMuon();
501 
502  std::cout<<" GL: "<<isGL<<" TR: "<<isTR<<" ST: "<<isST<<std::endl;
503  std::cout<<" nMatches "<<muonRef->numberOfMatches()<<std::endl;
504 
505  if ( muonRef->isGlobalMuon() ){
506  reco::TrackRef combinedMu = muonRef->combinedMuon();
507  std::cout<<" GL, pt: " << combinedMu->pt()
508  << " +/- " << combinedMu->ptError()/combinedMu->pt()
509  << " chi**2 GBL : " << combinedMu->normalizedChi2()<<std::endl;
510  std::cout<< " Total Muon Hits : " << combinedMu->hitPattern().numberOfValidMuonHits()
511  << "/" << combinedMu->hitPattern().numberOfLostMuonHits()
512  << " DT Hits : " << combinedMu->hitPattern().numberOfValidMuonDTHits()
513  << "/" << combinedMu->hitPattern().numberOfLostMuonDTHits()
514  << " CSC Hits : " << combinedMu->hitPattern().numberOfValidMuonCSCHits()
515  << "/" << combinedMu->hitPattern().numberOfLostMuonCSCHits()
516  << " RPC Hits : " << combinedMu->hitPattern().numberOfValidMuonRPCHits()
517  << "/" << combinedMu->hitPattern().numberOfLostMuonRPCHits()<<std::endl;
518 
519  std::cout<<" # of Valid Tracker Hits "<<combinedMu->hitPattern().numberOfValidTrackerHits()<<std::endl;
520  std::cout<<" # of Valid Pixel Hits "<<combinedMu->hitPattern().numberOfValidPixelHits()<<std::endl;
521  }
522  if ( muonRef->isStandAloneMuon() ){
523  reco::TrackRef standAloneMu = muonRef->standAloneMuon();
524  std::cout<<" ST, pt: " << standAloneMu->pt()
525  << " +/- " << standAloneMu->ptError()/standAloneMu->pt()
526  << " eta : " << standAloneMu->eta()
527  << " DT Hits : " << standAloneMu->hitPattern().numberOfValidMuonDTHits()
528  << "/" << standAloneMu->hitPattern().numberOfLostMuonDTHits()
529  << " CSC Hits : " << standAloneMu->hitPattern().numberOfValidMuonCSCHits()
530  << "/" << standAloneMu->hitPattern().numberOfLostMuonCSCHits()
531  << " RPC Hits : " << standAloneMu->hitPattern().numberOfValidMuonRPCHits()
532  << "/" << standAloneMu->hitPattern().numberOfLostMuonRPCHits()
533  << " chi**2 STA : " << standAloneMu->normalizedChi2()<<std::endl;
534  }
535 
536 
537  if ( muonRef->isTrackerMuon() ){
538  reco::TrackRef trackerMu = muonRef->track();
539  const reco::Track& track = *trackerMu;
540  std::cout<<" TR, pt: " << trackerMu->pt()
541  << " +/- " << trackerMu->ptError()/trackerMu->pt()
542  << " chi**2 TR : " << trackerMu->normalizedChi2()<<std::endl;
543  std::cout<<" nTrackerHits "<<track.hitPattern().numberOfValidTrackerHits()<<std::endl;
544  std::cout<< "TMLastStationAngLoose "
545  << muon::isGoodMuon(*muonRef,muon::TMLastStationAngLoose) << std::endl
546  << "TMLastStationAngTight "
547  << muon::isGoodMuon(*muonRef,muon::TMLastStationAngTight) << std::endl
548  << "TMLastStationLoose "
549  << muon::isGoodMuon(*muonRef,muon::TMLastStationLoose) << std::endl
550  << "TMLastStationTight "
551  << muon::isGoodMuon(*muonRef,muon::TMLastStationTight) << std::endl
552  << "TMOneStationLoose "
553  << muon::isGoodMuon(*muonRef,muon::TMOneStationLoose) << std::endl
554  << "TMOneStationTight "
555  << muon::isGoodMuon(*muonRef,muon::TMOneStationTight) << std::endl
556  << "TMLastStationOptimizedLowPtLoose "
557  << muon::isGoodMuon(*muonRef,muon::TMLastStationOptimizedLowPtLoose) << std::endl
558  << "TMLastStationOptimizedLowPtTight "
559  << muon::isGoodMuon(*muonRef,muon::TMLastStationOptimizedLowPtTight) << std::endl
560  << "TMLastStationOptimizedBarrelLowPtLoose "
561  << muon::isGoodMuon(*muonRef,muon::TMLastStationOptimizedBarrelLowPtLoose) << std::endl
562  << "TMLastStationOptimizedBarrelLowPtTight "
563  << muon::isGoodMuon(*muonRef,muon::TMLastStationOptimizedBarrelLowPtTight) << std::endl
564  << std::endl;
565 
566  }
567 
568  std::cout<< "TM2DCompatibilityLoose "
569  << muon::isGoodMuon(*muonRef,muon::TM2DCompatibilityLoose) << std::endl
570  << "TM2DCompatibilityTight "
571  << muon::isGoodMuon(*muonRef,muon::TM2DCompatibilityTight) << std::endl;
572 
573 
574 
575  if ( muonRef->isGlobalMuon() && muonRef->isTrackerMuon() && muonRef->isStandAloneMuon() ){
576  reco::TrackRef combinedMu = muonRef->combinedMuon();
577  reco::TrackRef trackerMu = muonRef->track();
578  reco::TrackRef standAloneMu = muonRef->standAloneMuon();
579 
580  double sigmaCombined = combinedMu->ptError()/(combinedMu->pt()*combinedMu->pt());
581  double sigmaTracker = trackerMu->ptError()/(trackerMu->pt()*trackerMu->pt());
582  double sigmaStandAlone = standAloneMu->ptError()/(standAloneMu->pt()*standAloneMu->pt());
583 
584  bool combined = combinedMu->ptError()/combinedMu->pt() < 0.20;
585  bool tracker = trackerMu->ptError()/trackerMu->pt() < 0.20;
586  bool standAlone = standAloneMu->ptError()/standAloneMu->pt() < 0.20;
587 
588  double delta1 = combined && tracker ?
589  fabs(1./combinedMu->pt() -1./trackerMu->pt())
590  /sqrt(sigmaCombined*sigmaCombined + sigmaTracker*sigmaTracker) : 100.;
591  double delta2 = combined && standAlone ?
592  fabs(1./combinedMu->pt() -1./standAloneMu->pt())
593  /sqrt(sigmaCombined*sigmaCombined + sigmaStandAlone*sigmaStandAlone) : 100.;
594  double delta3 = standAlone && tracker ?
595  fabs(1./standAloneMu->pt() -1./trackerMu->pt())
596  /sqrt(sigmaStandAlone*sigmaStandAlone + sigmaTracker*sigmaTracker) : 100.;
597 
598  double delta =
599  standAloneMu->hitPattern().numberOfValidMuonDTHits()+
600  standAloneMu->hitPattern().numberOfValidMuonCSCHits() > 0 ?
601  std::min(delta3,std::min(delta1,delta2)) : std::max(delta3,std::max(delta1,delta2));
602 
603  std::cout << "delta = " << delta << " delta1 "<<delta1<<" delta2 "<<delta2<<" delta3 "<<delta3<<std::endl;
604 
605  double ratio =
606  combinedMu->ptError()/combinedMu->pt()
607  / (trackerMu->ptError()/trackerMu->pt());
608  //if ( ratio > 2. && delta < 3. ) std::cout << "ALARM ! " << ratio << ", " << delta << std::endl;
609  std::cout<<" ratio "<<ratio<<" combined mu pt "<<combinedMu->pt()<<std::endl;
610  //bool quality3 = ( combinedMu->pt() < 50. || ratio < 2. ) && delta < 3.;
611 
612 
613  }
614 
615  double sumPtR03 = muonRef->isolationR03().sumPt;
616  double emEtR03 = muonRef->isolationR03().emEt;
617  double hadEtR03 = muonRef->isolationR03().hadEt;
618  double relIsoR03 = (sumPtR03 + emEtR03 + hadEtR03)/muonRef->pt();
619  double sumPtR05 = muonRef->isolationR05().sumPt;
620  double emEtR05 = muonRef->isolationR05().emEt;
621  double hadEtR05 = muonRef->isolationR05().hadEt;
622  double relIsoR05 = (sumPtR05 + emEtR05 + hadEtR05)/muonRef->pt();
623  std::cout<<" 0.3 Radion Rel Iso: "<<relIsoR03<<" sumPt "<<sumPtR03<<" emEt "<<emEtR03<<" hadEt "<<hadEtR03<<std::endl;
624  std::cout<<" 0.5 Radion Rel Iso: "<<relIsoR05<<" sumPt "<<sumPtR05<<" emEt "<<emEtR05<<" hadEt "<<hadEtR05<<std::endl;
625  return;
626 
627 }
628 
629 
630 
631 std::vector<reco::Muon::MuonTrackTypePair> PFMuonAlgo::goodMuonTracks(const reco::MuonRef& muon,bool includeSA) {
632  return muonTracks(muon,includeSA,maxDPtOPt_);
633 }
634 
635 
636 
637 std::vector<reco::Muon::MuonTrackTypePair> PFMuonAlgo::muonTracks(const reco::MuonRef& muon,bool includeSA,double dpt) {
638 
639 
640 
641  std::vector<reco::Muon::MuonTrackTypePair> out;
642 
643 
644  if(muon->globalTrack().isNonnull() && muon->globalTrack()->pt()>0)
645  if(muon->globalTrack()->ptError()/muon->globalTrack()->pt()<dpt)
646  out.push_back(std::make_pair(muon->globalTrack(),reco::Muon::CombinedTrack));
647 
648  if(muon->innerTrack().isNonnull() && muon->innerTrack()->pt()>0)
649  if(muon->innerTrack()->ptError()/muon->innerTrack()->pt()<dpt)//Here Loose!@
650  out.push_back(std::make_pair(muon->innerTrack(),reco::Muon::InnerTrack));
651 
652  bool pickyExists=false;
653  double pickyDpt=99999.;
654  if(muon->pickyTrack().isNonnull() && muon->pickyTrack()->pt()>0) {
655  pickyDpt = muon->pickyTrack()->ptError()/muon->pickyTrack()->pt();
656  if(pickyDpt<dpt)
657  out.push_back(std::make_pair(muon->pickyTrack(),reco::Muon::Picky));
658  pickyExists=true;
659  }
660 
661  bool dytExists=false;
662  double dytDpt=99999.;
663  if(muon->dytTrack().isNonnull() && muon->dytTrack()->pt()>0) {
664  dytDpt = muon->dytTrack()->ptError()/muon->dytTrack()->pt();
665  if(dytDpt<dpt)
666  out.push_back(std::make_pair(muon->dytTrack(),reco::Muon::DYT));
667  dytExists=true;
668  }
669 
670  //Magic: TPFMS is not a really good track especially under misalignment
671  //IT is kind of crap because if mu system is displaced it can make a change
672  //So allow TPFMS if there is no picky or the error of tpfms is better than picky
673  //AND if there is no DYT or the error of tpfms is better than DYT
674  if(muon->tpfmsTrack().isNonnull() && muon->tpfmsTrack()->pt()>0) {
675  double tpfmsDpt = muon->tpfmsTrack()->ptError()/muon->tpfmsTrack()->pt();
676  if( ( (pickyExists && tpfmsDpt<pickyDpt) || (!pickyExists) ) &&
677  ( (dytExists && tpfmsDpt<dytDpt) || (!dytExists) ) &&
678  tpfmsDpt<dpt )
679  out.push_back(std::make_pair(muon->tpfmsTrack(),reco::Muon::TPFMS));
680  }
681 
682  if(includeSA && muon->outerTrack().isNonnull())
683  if(muon->outerTrack()->ptError()/muon->outerTrack()->pt()<dpt)
684  out.push_back(std::make_pair(muon->outerTrack(),reco::Muon::OuterTrack));
685 
686  return out;
687 
688 }
689 
690 
691 
692 
693 
695 
696 
697 
698 
699 bool PFMuonAlgo::reconstructMuon(reco::PFCandidate& candidate, const reco::MuonRef& muon, bool allowLoose) {
700  using namespace std;
701  using namespace reco;
702 
703  if (!muon.isNonnull())
704  return false;
705 
706 
707 
708  bool isMu=false;
709 
710  if(allowLoose)
711  isMu = isMuon(muon) || isLooseMuon(muon);
712  else
713  isMu = isMuon(muon);
714 
715  if( !isMu)
716  return false;
717 
718  //get the valid tracks(without standalone except we allow loose muons)
719  //MIKE: Here we need to be careful. If we have a muon inside a dense
720  //jet environment often the track is badly measured. In this case
721  //we should not apply Dpt/Pt<1
722 
723  std::vector<reco::Muon::MuonTrackTypePair> validTracks = goodMuonTracks(muon);
724  if (!allowLoose)
725  validTracks = goodMuonTracks(muon);
726  else
727  validTracks = muonTracks(muon);
728 
729  if( validTracks.size() ==0)
730  return false;
731 
732 
733 
734  //check what is the track used.Rerun TuneP
735  reco::Muon::MuonTrackTypePair bestTrackPair = muon::tevOptimized(*muon);
736 
737  TrackRef bestTrack = bestTrackPair.first;
738  MuonTrackType trackType = bestTrackPair.second;
739 
740 
741 
742  MuonTrackTypePair trackPairWithSmallestError = getTrackWithSmallestError(validTracks);
743  TrackRef trackWithSmallestError = trackPairWithSmallestError.first;
744 
745  if( trackType == reco::Muon::InnerTrack &&
746  (!bestTrack->quality(trackQuality_) ||
747  bestTrack->ptError()/bestTrack->pt()> errorCompScale_*trackWithSmallestError->ptError()/trackWithSmallestError->pt() )) {
748  bestTrack = trackWithSmallestError;
749  trackType = trackPairWithSmallestError.second;
750  }
751  else if (trackType != reco::Muon::InnerTrack &&
752  bestTrack->ptError()/bestTrack->pt()> errorCompScale_*trackWithSmallestError->ptError()/trackWithSmallestError->pt()) {
753  bestTrack = trackWithSmallestError;
754  trackType = trackPairWithSmallestError.second;
755 
756  }
757 
758 
759  changeTrack(candidate,std::make_pair(bestTrack,trackType));
760  candidate.setMuonRef( muon );
761 
762  return true;
763 }
764 
765 
766 
767 
768 void PFMuonAlgo::changeTrack(reco::PFCandidate& candidate,const MuonTrackTypePair& track) {
769  using namespace reco;
770  reco::TrackRef bestTrack = track.first;
771  MuonTrackType trackType = track.second;
772  //OK Now redefine the canddiate with that track
773  double px = bestTrack->px();
774  double py = bestTrack->py();
775  double pz = bestTrack->pz();
776  double energy = sqrt(bestTrack->p()*bestTrack->p() + 0.1057*0.1057);
777 
778  candidate.setCharge(bestTrack->charge()>0 ? 1 : -1);
779  candidate.setP4(math::XYZTLorentzVector(px,py,pz,energy));
780  candidate.setParticleType(reco::PFCandidate::mu);
781  // candidate.setTrackRef( bestTrack );
782  candidate.setMuonTrackType(trackType);
783  if(trackType == reco::Muon::InnerTrack)
784  candidate.setVertexSource( PFCandidate::kTrkMuonVertex );
785  else if(trackType == reco::Muon::CombinedTrack)
786  candidate.setVertexSource( PFCandidate::kComMuonVertex );
787  else if(trackType == reco::Muon::TPFMS)
788  candidate.setVertexSource( PFCandidate::kTPFMSMuonVertex );
789  else if(trackType == reco::Muon::Picky)
790  candidate.setVertexSource( PFCandidate::kPickyMuonVertex );
791  else if(trackType == reco::Muon::DYT)
792  candidate.setVertexSource( PFCandidate::kDYTMuonVertex );
793  }
794 
795 
796 reco::Muon::MuonTrackTypePair
797 PFMuonAlgo::getTrackWithSmallestError(const std::vector<reco::Muon::MuonTrackTypePair>& tracks) {
798  TrackPtErrorSorter sorter;
799  return *std::min_element(tracks.begin(),tracks.end(),sorter);
800 }
801 
802 
803 
804 
805 
806 void PFMuonAlgo::estimateEventQuantities(const reco::PFCandidateCollection* pfc)
807 {
808  //SUM ET
809  sumetPU_ = 0.0;
810  METX_=0.;
811  METY_=0.;
812  sumet_=0.0;
813  for(reco::PFCandidateCollection::const_iterator i = pfc->begin();i!=pfc->end();++i) {
814  sumet_+=i->pt();
815  METX_+=i->px();
816  METY_+=i->py();
817  }
818 
819 }
820 
821 
822 
823 
824 
825 void PFMuonAlgo::postClean(reco::PFCandidateCollection* cands) {
826  using namespace std;
827  using namespace reco;
828  if (!postCleaning_)
829  return;
830 
831  //Initialize vectors
832 
833  if(pfCosmicsMuonCleanedCandidates_.get() )
834  pfCosmicsMuonCleanedCandidates_->clear();
835  else
836  pfCosmicsMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );
837 
838  if(pfCleanedTrackerAndGlobalMuonCandidates_.get() )
839  pfCleanedTrackerAndGlobalMuonCandidates_->clear();
840  else
841  pfCleanedTrackerAndGlobalMuonCandidates_.reset( new reco::PFCandidateCollection );
842 
843  if( pfFakeMuonCleanedCandidates_.get() )
844  pfFakeMuonCleanedCandidates_->clear();
845  else
846  pfFakeMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );
847 
848 
849  if( pfPunchThroughMuonCleanedCandidates_.get() )
850  pfPunchThroughMuonCleanedCandidates_->clear();
851  else
852  pfPunchThroughMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );
853 
854  if( pfPunchThroughHadronCleanedCandidates_.get() )
855  pfPunchThroughHadronCleanedCandidates_->clear();
856  else
857  pfPunchThroughHadronCleanedCandidates_.reset( new reco::PFCandidateCollection );
858 
859 
860 
861  pfPunchThroughHadronCleanedCandidates_->clear();
862 
863  maskedIndices_.clear();
864 
865  //Estimate MET and SumET
866 
867  estimateEventQuantities(cands);
868 
869 
870  std::vector<int> muons;
871  std::vector<int> cosmics;
872  //get the muons
873  for(unsigned int i=0;i<cands->size();++i) {
874  const reco::PFCandidate& cand = (*cands)[i];
875  if ( cand.particleId() == reco::PFCandidate::mu )
876  muons.push_back(i);
877  else if ( cand.particleId() == reco::PFCandidate::h && cand.muonRef().isNonnull()) {
878  //MET cleaning for muons that are not high purity and became charged hadrons
879  if (cand.pt()>100.0 && (!cand.trackRef()->quality(trackQuality_)) && cand.muonRef()->isGlobalMuon() && cand.rawHcalEnergy()<0.05*cand.p()) {
880  maskedIndices_.push_back(i);
881  pfPunchThroughMuonCleanedCandidates_->push_back(cand);
882  }
883  }
884 
885  }
886  //Then sort the muon indicess by decsending pt
887 
888  IndexPtComparator comparator(cands);
889  std::sort(muons.begin(),muons.end(),comparator);
890 
891 
892  //first kill cosmics
893  double METXCosmics=0;
894  double METYCosmics=0;
895  double SUMETCosmics=0.0;
896 
897  for(unsigned int i=0;i<muons.size();++i) {
898  const PFCandidate& pfc = cands->at(muons[i]);
899  double origin=0.0;
900  if(vertices_->size()>0&& vertices_->at(0).isValid() && ! vertices_->at(0).isFake())
901  origin = pfc.muonRef()->muonBestTrack()->dxy(vertices_->at(0).position());
902 
903  if( origin> cosmicRejDistance_) {
904  cosmics.push_back(muons[i]);
905  METXCosmics +=pfc.px();
906  METYCosmics +=pfc.py();
907  SUMETCosmics +=pfc.pt();
908  }
909  }
910  double MET2Cosmics = METXCosmics*METXCosmics+METYCosmics*METYCosmics;
911 
912  if ( SUMETCosmics > (sumet_-sumetPU_)/eventFactorCosmics_ && MET2Cosmics < METX_*METX_+ METY_*METY_)
913  for(unsigned int i=0;i<cosmics.size();++i) {
914  maskedIndices_.push_back(cosmics[i]);
915  pfCosmicsMuonCleanedCandidates_->push_back(cands->at(cosmics[i]));
916  }
917 
918  //Loop on the muons candidates and clean
919  for(unsigned int i=0;i<muons.size();++i) {
920  if( cleanMismeasured(cands->at(muons[i]),muons[i]))
921  continue;
922  cleanPunchThroughAndFakes(cands->at(muons[i]),cands,muons[i]);
923 
924  }
925 
926 
927 
928  //OK Now do the hard job ->remove the candidates that were cleaned
929  removeDeadCandidates(cands,maskedIndices_);
930 
931 
932 
933 
934 }
935 
936 void PFMuonAlgo::addMissingMuons(edm::Handle<reco::MuonCollection> muons, reco::PFCandidateCollection* cands) {
937  if(!postCleaning_)
938  return;
939 
940 
941  if( pfAddedMuonCandidates_.get() )
942  pfAddedMuonCandidates_->clear();
943  else
944  pfAddedMuonCandidates_.reset( new reco::PFCandidateCollection );
945 
946 
947 
948  for ( unsigned imu = 0; imu < muons->size(); ++imu ) {
949  reco::MuonRef muonRef( muons, imu );
950  bool used = false;
951  bool hadron = false;
952  for(unsigned i=0; i<cands->size(); i++) {
953  const PFCandidate& pfc = cands->at(i);
954  if ( !pfc.trackRef().isNonnull() ) continue;
955  if ( pfc.trackRef().isNonnull() && pfc.trackRef() == muonRef->track() )
956  hadron = true;
957  if ( !pfc.muonRef().isNonnull() ) continue;
958 
959  if ( pfc.muonRef()->innerTrack() == muonRef->innerTrack())
960  used = true;
961  else {
962  // Check if the stand-alone muon is not a spurious copy of an existing muon
963  // (Protection needed for HLT)
964  if ( pfc.muonRef()->isStandAloneMuon() && muonRef->isStandAloneMuon() ) {
965  double dEta = pfc.muonRef()->standAloneMuon()->eta() - muonRef->standAloneMuon()->eta();
966  double dPhi = pfc.muonRef()->standAloneMuon()->phi() - muonRef->standAloneMuon()->phi();
967  double dR = sqrt(dEta*dEta + dPhi*dPhi);
968  if ( dR < 0.005 ) {
969  used = true;
970  }
971  }
972  }
973 
974  if ( used ) break;
975  }
976 
977  if ( used ||hadron||(!muonRef.isNonnull()) ) continue;
978 
979 
980  TrackMETComparator comparator(METX_,METY_);
981  //Low pt dont need to be cleaned
982 
983  std::vector<reco::Muon::MuonTrackTypePair> tracks = goodMuonTracks(muonRef,true);
984  //If there is at least 1 track choice try to change the track
985  if(tracks.size()>0) {
986 
987  //Find tracks that change dramatically MET or Pt
988  std::vector<reco::Muon::MuonTrackTypePair> tracksThatChangeMET = tracksPointingAtMET(tracks);
989  //From those tracks get the one with smallest MET
990  if (tracksThatChangeMET.size()>0) {
991  reco::Muon::MuonTrackTypePair bestTrackType = *std::min_element(tracksThatChangeMET.begin(),tracksThatChangeMET.end(),comparator);
992 
993  //Make sure it is not cosmic
994  if((vertices_->size()==0) ||bestTrackType.first->dz(vertices_->at(0).position())<cosmicRejDistance_){
995 
996  //make a pfcandidate
997  int charge = bestTrackType.first->charge()>0 ? 1 : -1;
998  math::XYZTLorentzVector momentum(bestTrackType.first->px(),
999  bestTrackType.first->py(),
1000  bestTrackType.first->pz(),
1001  sqrt(bestTrackType.first->p()*bestTrackType.first->p()+0.1057*0.1057));
1002 
1003  cands->push_back( PFCandidate( charge,
1004  momentum,
1005  reco::PFCandidate::mu ) );
1006 
1007  changeTrack(cands->back(),bestTrackType);
1008 
1009  if (muonRef->track().isNonnull() )
1010  cands->back().setTrackRef( muonRef->track() );
1011 
1012  cands->back().setMuonRef(muonRef);
1013 
1014 
1015  pfAddedMuonCandidates_->push_back(cands->back());
1016 
1017  }
1018  }
1019  }
1020  }
1021 }
1022 
1023 std::pair<double,double>
1024 PFMuonAlgo::getMinMaxMET2(const reco::PFCandidate&pfc) {
1025  std::vector<reco::Muon::MuonTrackTypePair> tracks = goodMuonTracks((pfc.muonRef()),true);
1026 
1027  double METXNO = METX_-pfc.px();
1028  double METYNO = METY_-pfc.py();
1029  std::vector<double> met2;
1030  for (unsigned int i=0;i<tracks.size();++i) {
1031  met2.push_back(pow(METXNO+tracks.at(i).first->px(),2)+pow(METYNO+tracks.at(i).first->py(),2));
1032  }
1033 
1034  //PROTECT for cases of only one track. If there is only one track it will crash .
1035  //Has never happened but could likely happen!
1036 
1037  if(tracks.size()>1)
1038  return std::make_pair(*std::min_element(met2.begin(),met2.end()),*std::max_element(met2.begin(),met2.end()));
1039  else
1040  return std::make_pair(0,0);
1041 }
1042 
1043 
1044 bool PFMuonAlgo::cleanMismeasured(reco::PFCandidate& pfc,unsigned int i ){
1045  using namespace std;
1046  using namespace reco;
1047  bool cleaned=false;
1048 
1049  //First define the MET without this guy
1050  double METNOX = METX_ - pfc.px();
1051  double METNOY = METY_ - pfc.py();
1052  double SUMETNO = sumet_ -pfc.pt();
1053 
1054  TrackMETComparator comparator(METNOX,METNOY);
1055  //Low pt dont need to be cleaned
1056  if (pfc.pt()<minPostCleaningPt_)
1057  return false;
1058  std::vector<reco::Muon::MuonTrackTypePair> tracks = goodMuonTracks(pfc.muonRef(),false);
1059 
1060 
1061 
1062  //If there is more than 1 track choice try to change the track
1063  if(tracks.size()>1) {
1064  //Find tracks that change dramatically MET or Pt
1065  std::vector<reco::Muon::MuonTrackTypePair> tracksThatChangeMET = tracksWithBetterMET(tracks,pfc);
1066  //From those tracks get the one with smallest MET
1067  if (tracksThatChangeMET.size()>0) {
1068  reco::Muon::MuonTrackTypePair bestTrackType = *std::min_element(tracksThatChangeMET.begin(),tracksThatChangeMET.end(),comparator);
1069  changeTrack(pfc,bestTrackType);
1070 
1071  pfCleanedTrackerAndGlobalMuonCandidates_->push_back(pfc);
1072  //update eventquantities
1073  METX_ = METNOX+pfc.px();
1074  METY_ = METNOY+pfc.py();
1075  sumet_=SUMETNO+pfc.pt();
1076 
1077  }
1078  }
1079 
1080  //Now attempt to kill it
1081  if (!(pfc.muonRef()->isGlobalMuon() && pfc.muonRef()->isTrackerMuon())) {
1082  //define MET significance and SUM ET
1083  double MET2 = METX_*METX_+METY_*METY_;
1084  double newMET2 = METNOX*METNOX+METNOY*METNOY;
1085  double METSig = sqrt(MET2)/sqrt(sumet_-sumetPU_);
1086  if( METSig>metSigForRejection_)
1087  if((newMET2 < MET2/metFactorRejection_) &&
1088  ((SUMETNO-sumetPU_)/(sumet_-sumetPU_)<eventFractionRejection_)) {
1089  pfFakeMuonCleanedCandidates_->push_back(pfc);
1090  maskedIndices_.push_back(i);
1091  METX_ = METNOX;
1092  METY_ = METNOY;
1093  sumet_=SUMETNO;
1094  cleaned=true;
1095  }
1096 
1097  }
1098  return cleaned;
1099 
1100 }
1101 
1102 std::vector<reco::Muon::MuonTrackTypePair>
1103 PFMuonAlgo::tracksWithBetterMET(const std::vector<reco::Muon::MuonTrackTypePair>& tracks ,const reco::PFCandidate& pfc) {
1104  std::vector<reco::Muon::MuonTrackTypePair> outputTracks;
1105 
1106  double METNOX = METX_ - pfc.px();
1107  double METNOY = METY_ - pfc.py();
1108  double SUMETNO = sumet_ -pfc.pt();
1109  double MET2 = METX_*METX_+METY_*METY_;
1110  double newMET2=0.0;
1111  double newSUMET=0.0;
1112  double METSIG = sqrt(MET2)/sqrt(sumet_-sumetPU_);
1113 
1114 
1115  if(METSIG>metSigForCleaning_)
1116  for( unsigned int i=0;i<tracks.size();++i) {
1117  //calculate new SUM ET and MET2
1118  newSUMET = SUMETNO+tracks.at(i).first->pt()-sumetPU_;
1119  newMET2 = pow(METNOX+tracks.at(i).first->px(),2)+pow(METNOY+tracks.at(i).first->py(),2);
1120 
1121  if(newSUMET/(sumet_-sumetPU_)>eventFractionCleaning_ && newMET2<MET2/metFactorCleaning_)
1122  outputTracks.push_back(tracks.at(i));
1123  }
1124 
1125 
1126  return outputTracks;
1127 }
1128 
1129 
1130 std::vector<reco::Muon::MuonTrackTypePair>
1131 PFMuonAlgo::tracksPointingAtMET(const std::vector<reco::Muon::MuonTrackTypePair>& tracks) {
1132  std::vector<reco::Muon::MuonTrackTypePair> outputTracks;
1133 
1134 
1135  double newMET2=0.0;
1136 
1137  for( unsigned int i=0;i<tracks.size();++i) {
1138  //calculate new SUM ET and MET2
1139  newMET2 = pow(METX_+tracks.at(i).first->px(),2)+pow(METY_+tracks.at(i).first->py(),2);
1140 
1141  if(newMET2<(METX_*METX_+METY_*METY_)/metFactorCleaning_)
1142  outputTracks.push_back(tracks.at(i));
1143  }
1144 
1145 
1146  return outputTracks;
1147 }
1148 
1149 
1150 
1151 void PFMuonAlgo::setInputsForCleaning(const reco::VertexCollection* vertices) {
1152  vertices_ = vertices;
1153 }
1154 
1155 bool PFMuonAlgo::cleanPunchThroughAndFakes(reco::PFCandidate&pfc,reco::PFCandidateCollection* cands,unsigned int imu ){
1156  using namespace reco;
1157 
1158  bool cleaned=false;
1159 
1160  if (pfc.pt()<minPostCleaningPt_)
1161  return false;
1162 
1163 
1164  double METXNO = METX_-pfc.pt();
1165  double METYNO = METY_-pfc.pt();
1166  double MET2NO = METXNO*METXNO+METYNO*METYNO;
1167  double MET2 = METX_*METX_+METY_*METY_;
1168  bool fake1=false;
1169 
1170  std::pair<double,double> met2 = getMinMaxMET2(pfc);
1171 
1172  //Check for Fakes at high pseudorapidity
1173  if(pfc.muonRef()->standAloneMuon().isNonnull())
1174  fake1 =fabs ( pfc.eta() ) > 2.15 &&
1175  met2.first<met2.second/2 &&
1176  MET2NO < MET2/metFactorHighEta_ &&
1177  pfc.muonRef()->standAloneMuon()->pt() < pfc.pt()/ptFactorHighEta_;
1178 
1179  double factor = std::max(2.,2000./(sumet_-pfc.pt()-sumetPU_));
1180  bool fake2 = ( pfc.pt()/(sumet_-sumetPU_) < 0.25 && MET2NO < MET2/metFactorFake_ && met2.first<met2.second/factor );
1181 
1182  bool punchthrough =pfc.p() > minPunchThroughMomentum_ &&
1183  pfc.rawHcalEnergy() > minPunchThroughEnergy_ &&
1184  pfc.rawEcalEnergy()+pfc.rawHcalEnergy() > pfc.p()/punchThroughFactor_ &&
1185  !isIsolatedMuon(pfc.muonRef()) && MET2NO < MET2/punchThroughMETFactor_;
1186 
1187 
1188  if(fake1 || fake2||punchthrough) {
1189  // Find the block of the muon
1190  const PFCandidate::ElementsInBlocks& eleInBlocks = pfc.elementsInBlocks();
1191  if ( eleInBlocks.size() ) {
1192  PFBlockRef blockRefMuon = eleInBlocks[0].first;
1193  unsigned indexMuon = eleInBlocks[0].second;
1194  if (eleInBlocks.size()>1)
1195  indexMuon = eleInBlocks[1].second;
1196 
1197  // Check if the muon gave rise to a neutral hadron
1198  double iHad = 1E9;
1199  bool hadron = false;
1200  for ( unsigned i = imu+1; i < cands->size(); ++i ) {
1201  const PFCandidate& pfcn = cands->at(i);
1202  const PFCandidate::ElementsInBlocks& ele = pfcn.elementsInBlocks();
1203  if ( !ele.size() ) {
1204  continue;
1205  }
1206  PFBlockRef blockRefHadron = ele[0].first;
1207  unsigned indexHadron = ele[0].second;
1208  // We are out of the block -> exit the loop
1209  if ( blockRefHadron.key() != blockRefMuon.key() ) break;
1210  // Check that this particle is a neutral hadron
1211  if ( indexHadron == indexMuon &&
1212  pfcn.particleId() == reco::PFCandidate::h0 ) {
1213  iHad = i;
1214  hadron = true;
1215  }
1216  if ( hadron ) break;
1217  }
1218 
1219  if ( hadron ) {
1220 
1221  double rescaleFactor = cands->at(iHad).p()/cands->at(imu).p();
1222  METX_ -= cands->at(imu).px() + cands->at(iHad).px();
1223  METY_ -= cands->at(imu).py() + cands->at(iHad).py();
1224  sumet_ -=cands->at(imu).pt();
1225  cands->at(imu).rescaleMomentum(rescaleFactor);
1226  maskedIndices_.push_back(iHad);
1227  pfPunchThroughHadronCleanedCandidates_->push_back(cands->at(iHad));
1228  cands->at(imu).setParticleType(reco::PFCandidate::h);
1229  pfPunchThroughMuonCleanedCandidates_->push_back(cands->at(imu));
1230  METX_ += cands->at(imu).px();
1231  METY_ += cands->at(imu).py();
1232  sumet_ += cands->at(imu).pt();
1233 
1234  } else if ( fake1 || fake2 ) {
1235  METX_ -= cands->at(imu).px();
1236  METY_ -= cands->at(imu).py();
1237  sumet_ -= cands->at(imu).pt();
1238  maskedIndices_.push_back(imu);
1239  pfFakeMuonCleanedCandidates_->push_back(cands->at(imu));
1240  cleaned=true;
1241  }
1242  }
1243  }
1244  return cleaned;
1245 }
1246 
1247 
1248 void PFMuonAlgo::removeDeadCandidates(reco::PFCandidateCollection* obj, const std::vector<unsigned int>& indices)
1249 {
1250  size_t N = indices.size();
1251  size_t collSize = obj->size();
1252 
1253  for (size_t i = 0 ; i < N ; ++i)
1254  obj->at(indices.at(i)) = obj->at(collSize-i-1);
1255 
1256  obj->resize(collSize - indices.size());
1257 }
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:1044
reco::LeafCandidate::p
virtual double p() const
magnitude of momentum vector
Definition: LeafCandidate.h:108
PFMuonAlgo::isTightMuonPOG
static bool isTightMuonPOG(const reco::MuonRef &muonRef)
Definition: PFMuonAlgo.cc:472
digitizers_cfi.muon
tuple muon
Definition: digitizers_cfi.py:36
PFMuonAlgo::removeDeadCandidates
void removeDeadCandidates(reco::PFCandidateCollection *, const std::vector< unsigned int > &)
Definition: PFMuonAlgo.cc:1248
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:224
PFMuonAlgo::addMissingMuons
void addMissingMuons(edm::Handle< reco::MuonCollection >, reco::PFCandidateCollection *cands)
Definition: PFMuonAlgo.cc:936
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:93
PFMuonAlgo::setParameters
void setParameters(const edm::ParameterSet &)
Definition: PFMuonAlgo.cc:29
assert
assert(m_qm.get())
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:1155
reco::LeafCandidate::setP4
virtual void setP4(const LorentzVector &p4)
set 4-momentum
Definition: LeafCandidate.h:143
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:494
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:1131
PFMuonAlgo::isGlobalTightMuon
static bool isGlobalTightMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:185
reco::LeafCandidate::eta
virtual double eta() const
momentum pseudorapidity
Definition: LeafCandidate.h:137
reco::LeafCandidate::pt
virtual double pt() const
transverse momentum
Definition: LeafCandidate.h:131
reco::PFCandidate::ElementsInBlocks
std::vector< ElementInBlock > ElementsInBlocks
Definition: PFCandidate.h:386
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
muon::isLooseMuon
bool isLooseMuon(const reco::Muon &)
Definition: MuonSelectors.cc:762
reco::PFCandidate::trackRef
reco::TrackRef trackRef() const
Definition: PFCandidate.cc:433
PFMuonAlgo::getTrackWithSmallestError
MuonTrackTypePair getTrackWithSmallestError(const std::vector< MuonTrackTypePair > &)
Definition: PFMuonAlgo.cc:797
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:406
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:356
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:631
HLT_25ns14e33_v1_cff.vertices
tuple vertices
Definition: HLT_25ns14e33_v1_cff.py:4396
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:806
PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.dR
tuple dR
Definition: PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py:19
PFMuonAlgo::changeTrack
void changeTrack(reco::PFCandidate &, const MuonTrackTypePair &)
Definition: PFMuonAlgo.cc:768
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:125
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:1024
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
HLT_25ns14e33_v1_cff.quality
tuple quality
Definition: HLT_25ns14e33_v1_cff.py:33032
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:373
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:1103
PFMuonAlgo::postClean
void postClean(reco::PFCandidateCollection *)
Definition: PFMuonAlgo.cc:825
reco::PFCandidate::elementsInBlocks
const ElementsInBlocks & elementsInBlocks() const
Definition: PFCandidate.cc:682
PFMuonAlgo::setInputsForCleaning
void setInputsForCleaning(const reco::VertexCollection *)
Definition: PFMuonAlgo.cc:1151
PFMuonAlgo::muonTracks
std::vector< reco::Muon::MuonTrackTypePair > muonTracks(const reco::MuonRef &muon, bool includeSA=false, double dpt=1e+9)
Definition: PFMuonAlgo.cc:637
reco::LeafCandidate::py
virtual double py() const
y coordinate of momentum vector
Definition: LeafCandidate.h:127
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:234
PFMuonAlgo::reconstructMuon
bool reconstructMuon(reco::PFCandidate &, const reco::MuonRef &, bool allowLoose=false)
Definition: PFMuonAlgo.cc:699
Generated for CMSSW Reference Manual by doxygen 1.8.5