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())
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())
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  if(muon->pickyTrack().isNonnull()) {
656  if(muon->pickyTrack()->ptError()/muon->pickyTrack()->pt()<dpt)
657  out.push_back(std::make_pair(muon->pickyTrack(),reco::Muon::Picky));
658  pickyExists=true;
659  }
660 
661  //Magic: TPFMS is not a really good track especially under misalignment
662  //IT is kind of crap because if mu system is displaced it can make a change
663  //So allow TPFMS if there is no picky or the error of tpfms is better than picky
664  if(muon->tpfmsTrack().isNonnull() && ((pickyExists && muon->tpfmsTrack()->ptError()/muon->tpfmsTrack()->pt()<muon->pickyTrack()->ptError()/muon->pickyTrack()->pt())||(!pickyExists)) )
665  if(muon->tpfmsTrack()->ptError()/muon->tpfmsTrack()->pt()<dpt)
666  out.push_back(std::make_pair(muon->tpfmsTrack(),reco::Muon::TPFMS));
667 
668  if(includeSA && muon->outerTrack().isNonnull())
669  if(muon->outerTrack()->ptError()/muon->outerTrack()->pt()<dpt)
670  out.push_back(std::make_pair(muon->outerTrack(),reco::Muon::OuterTrack));
671 
672  return out;
673 
674 }
675 
676 
677 
678 
679 
681 
682 
683 
684 
685 bool PFMuonAlgo::reconstructMuon(reco::PFCandidate& candidate, const reco::MuonRef& muon, bool allowLoose) {
686  using namespace std;
687  using namespace reco;
688 
689  if (!muon.isNonnull())
690  return false;
691 
692 
693 
694  bool isMu=false;
695 
696  if(allowLoose)
697  isMu = isMuon(muon) || isLooseMuon(muon);
698  else
699  isMu = isMuon(muon);
700 
701  if( !isMu)
702  return false;
703 
704  //get the valid tracks(without standalone except we allow loose muons)
705  //MIKE: Here we need to be careful. If we have a muon inside a dense
706  //jet environment often the track is badly measured. In this case
707  //we should not apply Dpt/Pt<1
708 
709  std::vector<reco::Muon::MuonTrackTypePair> validTracks = goodMuonTracks(muon);
710  if (!allowLoose)
711  validTracks = goodMuonTracks(muon);
712  else
713  validTracks = muonTracks(muon);
714 
715  if( validTracks.size() ==0)
716  return false;
717 
718 
719 
720  //check what is the track used.Rerun TuneP
721  reco::Muon::MuonTrackTypePair bestTrackPair = muon::tevOptimized(*muon);
722 
723  TrackRef bestTrack = bestTrackPair.first;
724  MuonTrackType trackType = bestTrackPair.second;
725 
726 
727 
728  MuonTrackTypePair trackPairWithSmallestError = getTrackWithSmallestError(validTracks);
729  TrackRef trackWithSmallestError = trackPairWithSmallestError.first;
730 
731  if( trackType == reco::Muon::InnerTrack &&
732  (!bestTrack->quality(trackQuality_) ||
733  bestTrack->ptError()/bestTrack->pt()> errorCompScale_*trackWithSmallestError->ptError()/trackWithSmallestError->pt() )) {
734  bestTrack = trackWithSmallestError;
735  trackType = trackPairWithSmallestError.second;
736  }
737  else if (trackType != reco::Muon::InnerTrack &&
738  bestTrack->ptError()/bestTrack->pt()> errorCompScale_*trackWithSmallestError->ptError()/trackWithSmallestError->pt()) {
739  bestTrack = trackWithSmallestError;
740  trackType = trackPairWithSmallestError.second;
741 
742  }
743 
744 
745  changeTrack(candidate,std::make_pair(bestTrack,trackType));
746  candidate.setMuonRef( muon );
747 
748  return true;
749 }
750 
751 
752 
753 
754 void PFMuonAlgo::changeTrack(reco::PFCandidate& candidate,const MuonTrackTypePair& track) {
755  using namespace reco;
756  reco::TrackRef bestTrack = track.first;
757  MuonTrackType trackType = track.second;
758  //OK Now redefine the canddiate with that track
759  double px = bestTrack->px();
760  double py = bestTrack->py();
761  double pz = bestTrack->pz();
762  double energy = sqrt(bestTrack->p()*bestTrack->p() + 0.13957*0.13957);
763 
764  candidate.setCharge(bestTrack->charge()>0 ? 1 : -1);
765  candidate.setP4(math::XYZTLorentzVector(px,py,pz,energy));
766  candidate.setParticleType(reco::PFCandidate::mu);
767  // candidate.setTrackRef( bestTrack );
768  candidate.setMuonTrackType(trackType);
769  if(trackType == reco::Muon::InnerTrack)
770  candidate.setVertexSource( PFCandidate::kTrkMuonVertex );
771  else if(trackType == reco::Muon::CombinedTrack)
772  candidate.setVertexSource( PFCandidate::kComMuonVertex );
773  else if(trackType == reco::Muon::TPFMS)
774  candidate.setVertexSource( PFCandidate::kTPFMSMuonVertex );
775  else if(trackType == reco::Muon::Picky)
776  candidate.setVertexSource( PFCandidate::kPickyMuonVertex );
777  }
778 
779 
780 reco::Muon::MuonTrackTypePair
781 PFMuonAlgo::getTrackWithSmallestError(const std::vector<reco::Muon::MuonTrackTypePair>& tracks) {
782  TrackPtErrorSorter sorter;
783  return *std::min_element(tracks.begin(),tracks.end(),sorter);
784 }
785 
786 
787 
788 
789 
790 void PFMuonAlgo::estimateEventQuantities(const reco::PFCandidateCollection* pfc)
791 {
792  //SUM ET from PU
793  sumetPU_ = 0.0;
794  METX_=0.;
795  METY_=0.;
796  for (unsigned short i=1 ;i<vertices_->size();++i ) {
797  if ( !vertices_->at(i).isValid() || vertices_->at(i).isFake() ) continue;
798  vertices_->at(i);
799  for ( reco::Vertex::trackRef_iterator itr = vertices_->at(i).tracks_begin();
800  itr < vertices_->at(i).tracks_end(); ++itr ) {
801  sumetPU_ += (*itr)->pt();
802  }
803  }
804  sumetPU_ /= 0.65;
805  //SUM ET and MET
806  sumet_=0.0;
807  double METXCh=0.0;
808  double METYCh=0.0;
809  double METXNeut=0.0;
810  double METYNeut=0.0;
811 
812 
813  for(reco::PFCandidateCollection::const_iterator i = pfc->begin();i!=pfc->end();++i) {
814  sumet_+=i->pt();
815 
816  if (vertices_->size()>0 && vertices_->at(0).isValid()&& !vertices_->at(0).isFake()) {
817  //If charged and from PV or muon
818  if( (i->charge() !=0 && i->trackRef().isNonnull() && vertices_->size()>0&& i->trackRef()->dz(vertices_->at(0).position())<dzPV_)||(abs(i->pdgId())==13)) {
819  METXCh+=i->px();
820  METYCh+=i->py();
821  }
822  //If charged and not from PV(assume there is a neutral balancing it)
823  else if( i->charge() !=0 && i->trackRef().isNonnull() && i->trackRef()->dz(vertices_->at(0).position())>dzPV_) {
824  METXNeut-=i->px();
825  METYNeut-=i->py();
826  }
827  //Neutral
828  else if( !(i->charge() !=0 && i->trackRef().isNonnull())) {
829  METXNeut+=i->px();
830  METYNeut+=i->py();
831  }
832  } //else if we dont have a vertex make standard PFMET
833  else {
834  METXCh+=i->px();
835  METYCh+=i->py();
836  }
837  METX_ = (METXCh+METXNeut);
838  METY_ = (METYCh+METYNeut);
839  }
840 
841 }
842 
843 
844 
845 
846 void PFMuonAlgo::postClean(reco::PFCandidateCollection* cands) {
847  using namespace std;
848  using namespace reco;
849  if (!postCleaning_)
850  return;
851 
852  //Initialize vectors
853 
854  if(pfCosmicsMuonCleanedCandidates_.get() )
855  pfCosmicsMuonCleanedCandidates_->clear();
856  else
857  pfCosmicsMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );
858 
859  if(pfCleanedTrackerAndGlobalMuonCandidates_.get() )
860  pfCleanedTrackerAndGlobalMuonCandidates_->clear();
861  else
862  pfCleanedTrackerAndGlobalMuonCandidates_.reset( new reco::PFCandidateCollection );
863 
864  if( pfFakeMuonCleanedCandidates_.get() )
865  pfFakeMuonCleanedCandidates_->clear();
866  else
867  pfFakeMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );
868 
869 
870  if( pfPunchThroughMuonCleanedCandidates_.get() )
871  pfPunchThroughMuonCleanedCandidates_->clear();
872  else
873  pfPunchThroughMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );
874 
875  if( pfPunchThroughHadronCleanedCandidates_.get() )
876  pfPunchThroughHadronCleanedCandidates_->clear();
877  else
878  pfPunchThroughHadronCleanedCandidates_.reset( new reco::PFCandidateCollection );
879 
880 
881 
882  pfPunchThroughHadronCleanedCandidates_->clear();
883 
884  maskedIndices_.clear();
885 
886  //Estimate MET and SumET
887 
888  estimateEventQuantities(cands);
889 
890 
891  std::vector<int> muons;
892  std::vector<int> cosmics;
893  //get the muons
894  for(unsigned int i=0;i<cands->size();++i)
895  if ( cands->at(i).particleId() == reco::PFCandidate::mu )
896  muons.push_back(i);
897 
898  //Then sort the muon indicess by decsending pt
899 
900  IndexPtComparator comparator(cands);
901  std::sort(muons.begin(),muons.end(),comparator);
902 
903 
904  //first kill cosmics
905  double METXCosmics=0;
906  double METYCosmics=0;
907  double SUMETCosmics=0.0;
908 
909  for(unsigned int i=0;i<muons.size();++i) {
910  const PFCandidate& pfc = cands->at(muons[i]);
911  double origin=0.0;
912  if(vertices_->size()>0&& vertices_->at(0).isValid() && ! vertices_->at(0).isFake())
913  origin = pfc.muonRef()->muonBestTrack()->dxy(vertices_->at(0).position());
914 
915  if( origin> cosmicRejDistance_) {
916  cosmics.push_back(muons[i]);
917  METXCosmics +=pfc.px();
918  METYCosmics +=pfc.py();
919  SUMETCosmics +=pfc.pt();
920  }
921  }
922  double MET2Cosmics = METXCosmics*METXCosmics+METYCosmics*METYCosmics;
923 
924  if ( SUMETCosmics > (sumet_-sumetPU_)/eventFactorCosmics_ && MET2Cosmics < METX_*METX_+ METY_*METY_)
925  for(unsigned int i=0;i<cosmics.size();++i)
926  pfCosmicsMuonCleanedCandidates_->push_back(cands->at(muons[i]));
927 
928 
929  //Loop on the muons candidates and clean
930  for(unsigned int i=0;i<muons.size();++i) {
931 
932  if( cleanMismeasured(cands->at(muons[i]),muons[i]))
933  continue;
934  cleanPunchThroughAndFakes(cands->at(muons[i]),cands,muons[i]);
935 
936  }
937 
938 
939 
940  //OK Now do the hard job ->remove the candidates that were cleaned
941  removeDeadCandidates(cands,maskedIndices_);
942 
943 
944 
945 
946 }
947 
948 void PFMuonAlgo::addMissingMuons(edm::Handle<reco::MuonCollection> muons, reco::PFCandidateCollection* cands) {
949  if(!postCleaning_)
950  return;
951 
952 
953  if( pfAddedMuonCandidates_.get() )
954  pfAddedMuonCandidates_->clear();
955  else
956  pfAddedMuonCandidates_.reset( new reco::PFCandidateCollection );
957 
958 
959 
960  for ( unsigned imu = 0; imu < muons->size(); ++imu ) {
961  reco::MuonRef muonRef( muons, imu );
962  bool used = false;
963  bool hadron = false;
964  for(unsigned i=0; i<cands->size(); i++) {
965  const PFCandidate& pfc = cands->at(i);
966  if ( !pfc.trackRef().isNonnull() ) continue;
967  if ( pfc.trackRef().isNonnull() && pfc.trackRef() == muonRef->track() )
968  hadron = true;
969  if ( !pfc.muonRef().isNonnull() ) continue;
970 
971  if ( pfc.muonRef()->innerTrack() == muonRef->innerTrack())
972  used = true;
973  else {
974  // Check if the stand-alone muon is not a spurious copy of an existing muon
975  // (Protection needed for HLT)
976  if ( pfc.muonRef()->isStandAloneMuon() && muonRef->isStandAloneMuon() ) {
977  double dEta = pfc.muonRef()->standAloneMuon()->eta() - muonRef->standAloneMuon()->eta();
978  double dPhi = pfc.muonRef()->standAloneMuon()->phi() - muonRef->standAloneMuon()->phi();
979  double dR = sqrt(dEta*dEta + dPhi*dPhi);
980  if ( dR < 0.005 ) {
981  used = true;
982  }
983  }
984  }
985 
986  if ( used ) break;
987  }
988 
989  if ( used ||hadron||(!muonRef.isNonnull()) ) continue;
990 
991 
992  TrackMETComparator comparator(METX_,METY_);
993  //Low pt dont need to be cleaned
994 
995  std::vector<reco::Muon::MuonTrackTypePair> tracks = goodMuonTracks(muonRef,true);
996  //If there is at least 1 track choice try to change the track
997  if(tracks.size()>0) {
998 
999  //Find tracks that change dramatically MET or Pt
1000  std::vector<reco::Muon::MuonTrackTypePair> tracksThatChangeMET = tracksPointingAtMET(tracks);
1001  //From those tracks get the one with smallest MET
1002  if (tracksThatChangeMET.size()>0) {
1003  reco::Muon::MuonTrackTypePair bestTrackType = *std::min_element(tracksThatChangeMET.begin(),tracksThatChangeMET.end(),comparator);
1004 
1005  //Make sure it is not cosmic
1006  if((vertices_->size()==0) ||bestTrackType.first->dz(vertices_->at(0).position())<cosmicRejDistance_){
1007 
1008  //make a pfcandidate
1009  int charge = bestTrackType.first->charge()>0 ? 1 : -1;
1010  math::XYZTLorentzVector momentum(bestTrackType.first->px(),
1011  bestTrackType.first->py(),
1012  bestTrackType.first->pz(),
1013  sqrt(bestTrackType.first->p()*bestTrackType.first->p()+0.1057*0.1057));
1014 
1015  cands->push_back( PFCandidate( charge,
1016  momentum,
1017  reco::PFCandidate::mu ) );
1018 
1019  changeTrack(cands->back(),bestTrackType);
1020 
1021  if (muonRef->track().isNonnull() )
1022  cands->back().setTrackRef( muonRef->track() );
1023 
1024  cands->back().setMuonRef(muonRef);
1025 
1026 
1027  pfAddedMuonCandidates_->push_back(cands->back());
1028 
1029  }
1030  }
1031  }
1032  }
1033 }
1034 
1035 std::pair<double,double>
1036 PFMuonAlgo::getMinMaxMET2(const reco::PFCandidate&pfc) {
1037  std::vector<reco::Muon::MuonTrackTypePair> tracks = goodMuonTracks((pfc.muonRef()),true);
1038 
1039  double METXNO = METX_-pfc.px();
1040  double METYNO = METY_-pfc.py();
1041  std::vector<double> met2;
1042  for (unsigned int i=0;i<tracks.size();++i) {
1043  met2.push_back(pow(METXNO+tracks.at(i).first->px(),2)+pow(METYNO+tracks.at(i).first->py(),2));
1044  }
1045 
1046  //PROTECT for cases of only one track. If there is only one track it will crash .
1047  //Has never happened but could likely happen!
1048 
1049  if(tracks.size()>1)
1050  return std::make_pair(*std::min_element(met2.begin(),met2.end()),*std::max_element(met2.begin(),met2.end()));
1051  else
1052  return std::make_pair(0,0);
1053 }
1054 
1055 
1056 bool PFMuonAlgo::cleanMismeasured(reco::PFCandidate& pfc,unsigned int i ){
1057  using namespace std;
1058  using namespace reco;
1059  bool cleaned=false;
1060 
1061  //First define the MET without this guy
1062  double METNOX = METX_ - pfc.px();
1063  double METNOY = METY_ - pfc.py();
1064  double SUMETNO = sumet_ -pfc.pt();
1065 
1066  TrackMETComparator comparator(METNOX,METNOY);
1067  //Low pt dont need to be cleaned
1068  if (pfc.pt()<minPostCleaningPt_)
1069  return false;
1070  std::vector<reco::Muon::MuonTrackTypePair> tracks = goodMuonTracks(pfc.muonRef(),false);
1071 
1072 
1073 
1074  //If there is more than 1 track choice try to change the track
1075  if(tracks.size()>1) {
1076  //Find tracks that change dramatically MET or Pt
1077  std::vector<reco::Muon::MuonTrackTypePair> tracksThatChangeMET = tracksWithBetterMET(tracks,pfc);
1078  //From those tracks get the one with smallest MET
1079  if (tracksThatChangeMET.size()>0) {
1080  reco::Muon::MuonTrackTypePair bestTrackType = *std::min_element(tracksThatChangeMET.begin(),tracksThatChangeMET.end(),comparator);
1081  changeTrack(pfc,bestTrackType);
1082 
1083  pfCleanedTrackerAndGlobalMuonCandidates_->push_back(pfc);
1084  //update eventquantities
1085  METX_ = METNOX+pfc.px();
1086  METY_ = METNOY+pfc.py();
1087  sumet_=SUMETNO+pfc.pt();
1088 
1089  }
1090  }
1091 
1092  //Now attempt to kill it
1093  if (!(pfc.muonRef()->isGlobalMuon() && pfc.muonRef()->isTrackerMuon())) {
1094  //define MET significance and SUM ET
1095  double MET2 = METX_*METX_+METY_*METY_;
1096  double newMET2 = METNOX*METNOX+METNOY*METNOY;
1097  double METSig = sqrt(MET2)/sqrt(sumet_-sumetPU_);
1098  if( METSig>metSigForRejection_)
1099  if((newMET2 < MET2/metFactorRejection_) &&
1100  ((SUMETNO-sumetPU_)/(sumet_-sumetPU_)<eventFractionRejection_)) {
1101  pfFakeMuonCleanedCandidates_->push_back(pfc);
1102  maskedIndices_.push_back(i);
1103  METX_ = METNOX;
1104  METY_ = METNOY;
1105  sumet_=SUMETNO;
1106  cleaned=true;
1107  }
1108 
1109  }
1110  return cleaned;
1111 
1112 }
1113 
1114 std::vector<reco::Muon::MuonTrackTypePair>
1115 PFMuonAlgo::tracksWithBetterMET(const std::vector<reco::Muon::MuonTrackTypePair>& tracks ,const reco::PFCandidate& pfc) {
1116  std::vector<reco::Muon::MuonTrackTypePair> outputTracks;
1117 
1118  double METNOX = METX_ - pfc.px();
1119  double METNOY = METY_ - pfc.py();
1120  double SUMETNO = sumet_ -pfc.pt();
1121  double MET2 = METX_*METX_+METY_*METY_;
1122  double newMET2=0.0;
1123  double newSUMET=0.0;
1124  double METSIG = sqrt(MET2)/sqrt(sumet_-sumetPU_);
1125 
1126 
1127  if(METSIG>metSigForCleaning_)
1128  for( unsigned int i=0;i<tracks.size();++i) {
1129  //calculate new SUM ET and MET2
1130  newSUMET = SUMETNO+tracks.at(i).first->pt()-sumetPU_;
1131  newMET2 = pow(METNOX+tracks.at(i).first->px(),2)+pow(METNOY+tracks.at(i).first->py(),2);
1132 
1133  if(newSUMET/(sumet_-sumetPU_)>eventFractionCleaning_ && newMET2<MET2/metFactorCleaning_)
1134  outputTracks.push_back(tracks.at(i));
1135  }
1136 
1137 
1138  return outputTracks;
1139 }
1140 
1141 
1142 std::vector<reco::Muon::MuonTrackTypePair>
1143 PFMuonAlgo::tracksPointingAtMET(const std::vector<reco::Muon::MuonTrackTypePair>& tracks) {
1144  std::vector<reco::Muon::MuonTrackTypePair> outputTracks;
1145 
1146 
1147  double newMET2=0.0;
1148 
1149  for( unsigned int i=0;i<tracks.size();++i) {
1150  //calculate new SUM ET and MET2
1151  newMET2 = pow(METX_+tracks.at(i).first->px(),2)+pow(METY_+tracks.at(i).first->py(),2);
1152 
1153  if(newMET2<(METX_*METX_+METY_*METY_)/metFactorCleaning_)
1154  outputTracks.push_back(tracks.at(i));
1155  }
1156 
1157 
1158  return outputTracks;
1159 }
1160 
1161 
1162 
1163 void PFMuonAlgo::setInputsForCleaning(const reco::VertexCollection* vertices) {
1164  vertices_ = vertices;
1165 }
1166 
1167 bool PFMuonAlgo::cleanPunchThroughAndFakes(reco::PFCandidate&pfc,reco::PFCandidateCollection* cands,unsigned int imu ){
1168  using namespace reco;
1169 
1170  bool cleaned=false;
1171 
1172  if (pfc.pt()<minPostCleaningPt_)
1173  return false;
1174 
1175 
1176  double METXNO = METX_-pfc.pt();
1177  double METYNO = METY_-pfc.pt();
1178  double MET2NO = METXNO*METXNO+METYNO*METYNO;
1179  double MET2 = METX_*METX_+METY_*METY_;
1180  bool fake1=false;
1181 
1182  std::pair<double,double> met2 = getMinMaxMET2(pfc);
1183 
1184  //Check for Fakes at high pseudorapidity
1185  if(pfc.muonRef()->standAloneMuon().isNonnull())
1186  fake1 =fabs ( pfc.eta() ) > 2.15 &&
1187  met2.first<met2.second/2 &&
1188  MET2NO < MET2/metFactorHighEta_ &&
1189  pfc.muonRef()->standAloneMuon()->pt() < pfc.pt()/ptFactorHighEta_;
1190 
1191  double factor = std::max(2.,2000./(sumet_-pfc.pt()-sumetPU_));
1192  bool fake2 = ( pfc.pt()/(sumet_-sumetPU_) < 0.25 && MET2NO < MET2/metFactorFake_ && met2.first<met2.second/factor );
1193 
1194  bool punchthrough =pfc.p() > minPunchThroughMomentum_ &&
1195  pfc.rawHcalEnergy() > minPunchThroughEnergy_ &&
1196  pfc.rawEcalEnergy()+pfc.rawHcalEnergy() > pfc.p()/punchThroughFactor_ &&
1197  !isIsolatedMuon(pfc.muonRef()) && MET2NO < MET2/punchThroughMETFactor_;
1198 
1199 
1200  if(fake1 || fake2||punchthrough) {
1201  // Find the block of the muon
1202  const PFCandidate::ElementsInBlocks& eleInBlocks = pfc.elementsInBlocks();
1203  if ( eleInBlocks.size() ) {
1204  PFBlockRef blockRefMuon = eleInBlocks[0].first;
1205  unsigned indexMuon = eleInBlocks[0].second;
1206  for ( unsigned iele = 1; iele < eleInBlocks.size(); ++iele ) {
1207  indexMuon = eleInBlocks[iele].second;
1208  break;
1209  }
1210 
1211  // Check if the muon gave rise to a neutral hadron
1212  double iHad = 1E9;
1213  bool hadron = false;
1214  for ( unsigned i = imu+1; i < cands->size(); ++i ) {
1215  const PFCandidate& pfcn = cands->at(i);
1216  const PFCandidate::ElementsInBlocks& ele = pfcn.elementsInBlocks();
1217  if ( !ele.size() ) {
1218  continue;
1219  }
1220  PFBlockRef blockRefHadron = ele[0].first;
1221  unsigned indexHadron = ele[0].second;
1222  // We are out of the block -> exit the loop
1223  if ( blockRefHadron.key() != blockRefMuon.key() ) break;
1224  // Check that this particle is a neutral hadron
1225  if ( indexHadron == indexMuon &&
1226  pfcn.particleId() == reco::PFCandidate::h0 ) {
1227  iHad = i;
1228  hadron = true;
1229  }
1230  if ( hadron ) break;
1231  }
1232 
1233  if ( hadron ) {
1234 
1235  double rescaleFactor = cands->at(iHad).p()/cands->at(imu).p();
1236  METX_ -= cands->at(imu).px() + cands->at(iHad).px();
1237  METY_ -= cands->at(imu).py() + cands->at(iHad).py();
1238  sumet_ -=cands->at(imu).pt();
1239  cands->at(imu).rescaleMomentum(rescaleFactor);
1240  maskedIndices_.push_back(iHad);
1241  pfPunchThroughHadronCleanedCandidates_->push_back(cands->at(iHad));
1242  cands->at(imu).setParticleType(reco::PFCandidate::h);
1243  pfPunchThroughMuonCleanedCandidates_->push_back(cands->at(imu));
1244  METX_ += cands->at(imu).px();
1245  METY_ += cands->at(imu).py();
1246  sumet_ += cands->at(imu).pt();
1247 
1248  } else if ( fake1 || fake2 ) {
1249  METX_ -= cands->at(imu).px();
1250  METY_ -= cands->at(imu).py();
1251  sumet_ -= cands->at(imu).pt();
1252  maskedIndices_.push_back(imu);
1253  pfFakeMuonCleanedCandidates_->push_back(cands->at(imu));
1254  cleaned=true;
1255  }
1256  }
1257  }
1258  return cleaned;
1259 }
1260 
1261 
1262 void PFMuonAlgo::removeDeadCandidates(reco::PFCandidateCollection* obj, const std::vector<unsigned int>& indices)
1263 {
1264  size_t N = indices.size();
1265  size_t collSize = obj->size();
1266 
1267  for (size_t i = 0 ; i < N ; ++i)
1268  obj->at(indices.at(i)) = obj->at(collSize-i-1);
1269 
1270  obj->resize(collSize - indices.size());
1271 }
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:1056
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:1262
reco::PFCandidate::rawEcalEnergy
double rawEcalEnergy() const
return corrected Ecal energy
Definition: PFCandidate.h:207
PFMuonAlgo::addMissingMuons
void addMissingMuons(edm::Handle< reco::MuonCollection >, reco::PFCandidateCollection *cands)
Definition: PFMuonAlgo.cc:948
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:759
PFMuonAlgo::cleanPunchThroughAndFakes
bool cleanPunchThroughAndFakes(reco::PFCandidate &, reco::PFCandidateCollection *, unsigned int)
Definition: PFMuonAlgo.cc:1167
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
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:1143
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
PFMuonAlgo::isGlobalTightMuon
static bool isGlobalTightMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:185
reco::PFCandidate::ElementsInBlocks
std::vector< ElementInBlock > ElementsInBlocks
Definition: PFCandidate.h:368
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:429
PFMuonAlgo::getTrackWithSmallestError
MuonTrackTypePair getTrackWithSmallestError(const std::vector< MuonTrackTypePair > &)
Definition: PFMuonAlgo.cc:781
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
muon::tevOptimized
reco::Muon::MuonTrackTypePair tevOptimized(const reco::TrackRef &combinedTrack, const reco::TrackRef &trackerTrack, const reco::TrackRef &tpfmsTrack, const reco::TrackRef &pickyTrack, const double ptThreshold=200., const double tune1=17., const double tune2=40., const double dptcut=0.25)
Definition: MuonCocktails.cc:9
reco::PFCandidate::setVertexSource
void setVertexSource(PFVertexType vt)
Definition: PFCandidate.h:388
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:338
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:252
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
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
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:790
PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.dR
tuple dR
Definition: PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py:19
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:221
PFMuonAlgo::changeTrack
void changeTrack(reco::PFCandidate &, const MuonTrackTypePair &)
Definition: PFMuonAlgo.cc:754
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:446
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:46
N
#define N
Definition: blowfish.cc:9
reco::LeafCandidate::px
virtual double px() const
x coordinate of momentum vector
Definition: LeafCandidate.h:172
testEve_cfg.tracks
tuple tracks
Definition: testEve_cfg.py:39
PFMuonAlgo::isLooseMuon
static bool isLooseMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:166
edm::Ref::key
key_type key() const
Accessor for product key.
Definition: Ref.h:266
PFMuonAlgo::getMinMaxMET2
std::pair< double, double > getMinMaxMET2(const reco::PFCandidate &)
Definition: PFMuonAlgo.cc:1036
MCScenario_CRAFT1_22X.sorter
def sorter
Definition: MCScenario_CRAFT1_22X.py:96
reco::HitPattern::numberOfValidTrackerHits
int numberOfValidTrackerHits() const
Definition: HitPattern.h:593
reco::PFCandidate
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:38
reco::PFCandidate::setMuonRef
void setMuonRef(const reco::MuonRef &ref)
set muon reference
Definition: PFCandidate.cc:432
patZpeak.muons
tuple muons
Definition: patZpeak.py:38
reco::Vertex::trackRef_iterator
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:37
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:355
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:1115
PFMuonAlgo::postClean
void postClean(reco::PFCandidateCollection *)
Definition: PFMuonAlgo.cc:846
reco::PFCandidate::elementsInBlocks
const ElementsInBlocks & elementsInBlocks() const
Definition: PFCandidate.cc:675
PFMuonAlgo::setInputsForCleaning
void setInputsForCleaning(const reco::VertexCollection *)
Definition: PFMuonAlgo.cc:1163
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:217
PFMuonAlgo::reconstructMuon
bool reconstructMuon(reco::PFCandidate &, const reco::MuonRef &, bool allowLoose=false)
Definition: PFMuonAlgo.cc:685
Generated for CMSSW Reference Manual by doxygen 1.8.5