00001
00016
00017 #include <memory>
00018 #include <string>
00019
00020
00021 #include "FWCore/Framework/interface/Event.h"
00022 #include "FWCore/Framework/interface/EventSetup.h"
00023 #include "FWCore/Framework/interface/Frameworkfwd.h"
00024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00025 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
00026
00027 #include "DataFormats/TrackReco/interface/Track.h"
00028 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00029 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00030 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00031 #include "DataFormats/DetId/interface/DetId.h"
00032 #include "DataFormats/VertexReco/interface/Vertex.h"
00033 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00034
00035 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00036 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
00037 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00038 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00039 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
00040
00041 #include "RecoMuon/MuonIdentification/interface/MuonCosmicCompatibilityFiller.h"
00042 #include "RecoMuon/MuonIdentification/interface/MuonCosmicsId.h"
00043
00044 #include "DataFormats/MuonReco/interface/MuonCosmicCompatibility.h"
00045 #include "DataFormats/MuonReco/interface/Muon.h"
00046 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00047 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00048
00049 #include "TMath.h"
00050
00051
00052 using namespace edm;
00053 using namespace std;
00054
00055 MuonCosmicCompatibilityFiller::MuonCosmicCompatibilityFiller(const edm::ParameterSet& iConfig):
00056 inputMuonCollections_(iConfig.getParameter<std::vector<edm::InputTag> >("InputMuonCollections")),
00057 inputTrackCollections_(iConfig.getParameter<std::vector<edm::InputTag> >("InputTrackCollections")),
00058 inputCosmicMuonCollection_(iConfig.getParameter<edm::InputTag>("InputCosmicMuonCollection")),
00059 inputVertexCollection_(iConfig.getParameter<edm::InputTag>("InputVertexCollection")),
00060 service_(0)
00061 {
00062
00063 edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
00064 service_ = new MuonServiceProxy(serviceParameters);
00065
00066
00067 angleThreshold_ = iConfig.getParameter<double>("angleCut");
00068 deltaPt_ = iConfig.getParameter<double>("deltaPt");
00069
00070 offTimePosTightMult_ = iConfig.getParameter<double>("offTimePosTightMult");
00071 offTimeNegTightMult_ = iConfig.getParameter<double>("offTimeNegTightMult");
00072 offTimePosTight_ = iConfig.getParameter<double>("offTimePosTight");
00073 offTimeNegTight_ = iConfig.getParameter<double>("offTimeNegTight");
00074 offTimePosLooseMult_ = iConfig.getParameter<double>("offTimePosLooseMult");
00075 offTimeNegLooseMult_ = iConfig.getParameter<double>("offTimeNegLooseMult");
00076 offTimePosLoose_ = iConfig.getParameter<double>("offTimePosLoose");
00077 offTimeNegLoose_ = iConfig.getParameter<double>("offTimeNegLoose");
00078 corrTimeNeg_ = iConfig.getParameter<double>("corrTimeNeg");
00079 corrTimePos_ = iConfig.getParameter<double>("corrTimePos");
00080
00081 sharedHits_ = iConfig.getParameter<int>("sharedHits");
00082 sharedFrac_ = iConfig.getParameter<double>("sharedFrac");
00083 ipThreshold_ = iConfig.getParameter<double>("ipCut");
00084
00085 nChamberMatches_ = iConfig.getParameter<int>("nChamberMatches");
00086 segmentComp_ = iConfig.getParameter<double>("segmentComp");
00087
00088 maxdzLooseMult_ = iConfig.getParameter<double>("maxdzLooseMult");
00089 maxdxyLooseMult_ = iConfig.getParameter<double>("maxdxyLooseMult");
00090 maxdzTightMult_ = iConfig.getParameter<double>("maxdzTightMult");
00091 maxdxyTightMult_ = iConfig.getParameter<double>("maxdxyTightMult");
00092 maxdzLoose_ = iConfig.getParameter<double>("maxdzLoose");
00093 maxdxyLoose_ = iConfig.getParameter<double>("maxdxyLoose");
00094 maxdzTight_ = iConfig.getParameter<double>("maxdzTight");
00095 maxdxyTight_ = iConfig.getParameter<double>("maxdxyTight");
00096 largedxyMult_ = iConfig.getParameter<double>("largedxyMult");
00097 largedxy_ = iConfig.getParameter<double>("largedxy");
00098 hIpTrdxy_ = iConfig.getParameter<double>("hIpTrdxy");
00099 hIpTrvProb_ = iConfig.getParameter<double>("hIpTrvProb");
00100 minvProb_ = iConfig.getParameter<double>("minvProb");
00101 maxvertZ_ = iConfig.getParameter<double>("maxvertZ");
00102 maxvertRho_ = iConfig.getParameter<double>("maxvertRho");
00103
00104 }
00105
00106 MuonCosmicCompatibilityFiller::~MuonCosmicCompatibilityFiller() {
00107 if (service_) delete service_;
00108 }
00109
00110 reco::MuonCosmicCompatibility
00111 MuonCosmicCompatibilityFiller::fillCompatibility( const reco::Muon& muon, edm::Event& iEvent, const edm::EventSetup& iSetup )
00112 {
00113 const std::string theCategory = "MuonCosmicCompatibilityFiller";
00114
00115 reco::MuonCosmicCompatibility returnComp;
00116
00117 service_->update(iSetup);
00118
00119 float timeCompatibility = muonTiming(iEvent, muon, false);
00120 float backToBackCompatibility = backToBack2LegCosmic(iEvent,muon);
00121 float overlapCompatibility = isOverlappingMuon(iEvent,iSetup,muon);
00122 float ipCompatibility = pvMatches(iEvent,muon,false);
00123 float vertexCompatibility = eventActivity(iEvent,muon);
00124 float combinedCompatibility = combinedCosmicID(iEvent,iSetup,muon,false,false);
00125
00126 returnComp.timeCompatibility = timeCompatibility;
00127 returnComp.backToBackCompatibility = backToBackCompatibility;
00128 returnComp.overlapCompatibility = overlapCompatibility;
00129 returnComp.cosmicCompatibility = combinedCompatibility;
00130 returnComp.ipCompatibility = ipCompatibility;
00131 returnComp.vertexCompatibility = vertexCompatibility;
00132
00133 return returnComp;
00134
00135 }
00136
00137
00138
00139
00140 float
00141 MuonCosmicCompatibilityFiller::muonTiming(const edm::Event& iEvent, const reco::Muon& muon, bool isLoose) const {
00142
00143 float offTimeNegMult, offTimePosMult, offTimeNeg, offTimePos;
00144
00145 if (isLoose) {
00146
00147 offTimeNegMult = offTimeNegLooseMult_;
00148 offTimePosMult = offTimePosLooseMult_;
00149 offTimeNeg = offTimeNegLoose_;
00150 offTimePos = offTimePosLoose_;
00151 } else {
00152 offTimeNegMult = offTimeNegTightMult_;
00153 offTimePosMult = offTimePosTightMult_;
00154 offTimeNeg = offTimeNegTight_;
00155 offTimePos = offTimePosTight_;
00156 }
00157
00158 float result = 0.0;
00159
00160 if( muon.isTimeValid() ) {
00161
00162 if (nMuons(iEvent) > 1) {
00163
00164 float positiveTime = 0;
00165 if ( muon.time().timeAtIpInOut < offTimeNegMult || muon.time().timeAtIpInOut > offTimePosMult) result = 1.;
00166 if ( muon.time().timeAtIpInOut > 0.) positiveTime = muon.time().timeAtIpInOut;
00167
00168
00169
00170 if (!isLoose && result == 0 && positiveTime > corrTimePos_) {
00171
00172
00173 bool isUp = false;
00174 reco::TrackRef outertrack = muon.outerTrack();
00175 if( outertrack.isNonnull() ) {
00176 if( outertrack->phi() > 0 ) isUp = true;
00177
00178
00179 edm::Handle<reco::MuonCollection> muonHandle;
00180 iEvent.getByLabel(inputMuonCollections_[1], muonHandle);
00181
00182 if( !muonHandle.failedToGet() ) {
00183 for ( reco::MuonCollection::const_iterator iMuon = muonHandle->begin(); iMuon != muonHandle->end(); ++iMuon ) {
00184 if (!iMuon->isGlobalMuon()) continue;
00185
00186 reco::TrackRef checkedTrack = iMuon->outerTrack();
00187 if( muon.isTimeValid() ) {
00188
00189
00190 if (checkedTrack->phi() < 0 && isUp) {
00191 if (iMuon->time().timeAtIpInOut < corrTimeNeg_) result = 1.0;
00192 break;
00193 } else if (checkedTrack->phi() > 0 && !isUp) {
00194
00195 if (iMuon->time().timeAtIpInOut < corrTimeNeg_) result = 1.0;
00196 break;
00197 }
00198 }
00199 }
00200 }
00201 }
00202 }
00203 } else {
00204
00205 if ( muon.time().timeAtIpInOut < offTimeNeg || muon.time().timeAtIpInOut > offTimePos) result = 1.;
00206 }
00207 }
00208
00209
00210 if (!isLoose && result > 0) {
00211
00212 if (pvMatches(iEvent, muon, true) == 0) result *= 2.;
00213 }
00214
00215 return result;
00216 }
00217
00218
00219
00220
00221 unsigned int
00222 MuonCosmicCompatibilityFiller::backToBack2LegCosmic(const edm::Event& iEvent, const reco::Muon& muon) const {
00223
00224 unsigned int result = 0;
00225 reco::TrackRef track;
00226 if ( muon.isGlobalMuon() ) track = muon.innerTrack();
00227 else if ( muon.isTrackerMuon() ) track = muon.track();
00228 else if ( muon.isStandAloneMuon() ) return false;
00229
00230 for (unsigned int iColl = 0; iColl<inputTrackCollections_.size(); ++iColl){
00231 edm::Handle<reco::TrackCollection> trackHandle;
00232 iEvent.getByLabel(inputTrackCollections_[iColl],trackHandle);
00233 if (muonid::findOppositeTrack(trackHandle, *track, angleThreshold_, deltaPt_).isNonnull()) {
00234 result++;
00235 }
00236 }
00237
00238 return result;
00239 }
00240
00241
00242
00243
00244 unsigned int
00245 MuonCosmicCompatibilityFiller::nMuons(const edm::Event& iEvent) const {
00246
00247 unsigned int nGlb = 0;
00248
00249 edm::Handle<reco::MuonCollection> muonHandle;
00250 iEvent.getByLabel(inputMuonCollections_[1], muonHandle);
00251
00252 if( !muonHandle.failedToGet() ) {
00253 for ( reco::MuonCollection::const_iterator iMuon = muonHandle->begin(); iMuon != muonHandle->end(); ++iMuon ) {
00254 if (!iMuon->isGlobalMuon()) continue;
00255 nGlb++;
00256 }
00257 }
00258
00259 return nGlb;
00260 }
00261
00262
00263
00264
00265
00266 bool
00267 MuonCosmicCompatibilityFiller::isOverlappingMuon(const edm::Event& iEvent, const edm::EventSetup& iSetup, const reco::Muon& muon) const {
00268
00269
00270
00271
00272
00273
00274
00275
00276 bool overlappingMuon = false;
00277 if( !muon.isGlobalMuon() ) return false;
00278
00279
00280 edm::Handle<reco::MuonCollection> muonHandle;
00281 iEvent.getByLabel(inputCosmicMuonCollection_, muonHandle);
00282
00283
00284 ESHandle<GlobalTrackingGeometry> trackingGeometry;
00285 iSetup.get<GlobalTrackingGeometryRecord>().get(trackingGeometry);
00286
00287
00288 math::XYZPoint RefVtx;
00289 RefVtx.SetXYZ(0, 0, 0);
00290
00291 edm::Handle<reco::VertexCollection> pvHandle;
00292 iEvent.getByLabel(inputVertexCollection_,pvHandle);
00293 const reco::VertexCollection & vertices = *pvHandle.product();
00294 for(reco::VertexCollection::const_iterator it=vertices.begin() ; it!=vertices.end() ; ++it) {
00295 RefVtx = it->position();
00296 }
00297
00298
00299 if( !muonHandle.failedToGet() ) {
00300 for ( reco::MuonCollection::const_iterator cosmicMuon = muonHandle->begin();cosmicMuon != muonHandle->end(); ++cosmicMuon ) {
00301 if ( cosmicMuon->innerTrack() == muon.innerTrack() || cosmicMuon->outerTrack() == muon.outerTrack()) return true;
00302
00303 reco::TrackRef outertrack = muon.outerTrack();
00304 reco::TrackRef costrack = cosmicMuon->outerTrack();
00305
00306 bool isUp = false;
00307 if( outertrack->phi() > 0 ) isUp = true;
00308
00309
00310 int RecHitsMuon = outertrack->numberOfValidHits();
00311 int RecHitsCosmicMuon = 0;
00312 int shared = 0;
00313
00314 if( costrack.isNonnull() ) {
00315 int nhitsUp = 0;
00316 int nhitsDown = 0;
00317 bool isCosmic1Leg = false;
00318 bool isCloseIP = false;
00319 bool isCloseRef = false;
00320
00321 for( trackingRecHit_iterator coshit = costrack->recHitsBegin(); coshit != costrack->recHitsEnd(); coshit++ ) {
00322 if( (*coshit)->isValid() ) {
00323 DetId id((*coshit)->geographicalId());
00324 double hity = trackingGeometry->idToDet(id)->position().y();
00325 if( hity > 0 ) nhitsUp++;
00326 if( hity < 0 ) nhitsDown++;
00327
00328 if( isUp && hity > 0 ) RecHitsCosmicMuon++;
00329 if( !isUp && hity < 0 ) RecHitsCosmicMuon++;
00330 }
00331 }
00332
00333 if( nhitsUp > 0 && nhitsDown > 0 ) isCosmic1Leg = true;
00334
00335
00336 if( outertrack.isNonnull() ) {
00337
00338 const double ipErr = (double)outertrack->d0Error();
00339 double ipThreshold = max(ipThreshold_, ipErr);
00340 if( fabs(outertrack->dxy(RefVtx) + costrack->dxy(RefVtx)) < ipThreshold ) isCloseIP = true;
00341
00342
00343
00344 GlobalPoint muonRefVtx( outertrack->vx(), outertrack->vy(), outertrack->vz() );
00345 GlobalPoint cosmicRefVtx( costrack->vx(), costrack->vy(), costrack->vz() );
00346 float dist = (muonRefVtx - cosmicRefVtx).mag();
00347 if( dist < 0.1 ) isCloseRef = true;
00348
00349
00350 for( trackingRecHit_iterator trkhit = outertrack->recHitsBegin(); trkhit != outertrack->recHitsEnd(); trkhit++ ) {
00351 if( (*trkhit)->isValid() ) {
00352 for( trackingRecHit_iterator coshit = costrack->recHitsBegin(); coshit != costrack->recHitsEnd(); coshit++ ) {
00353 if( (*coshit)->isValid() ) {
00354 if( (*trkhit)->geographicalId() == (*coshit)->geographicalId() ) {
00355 if( ((*trkhit)->localPosition() - (*coshit)->localPosition()).mag()< 10e-5 ) shared++;
00356 }
00357
00358 }
00359 }
00360 }
00361 }
00362 }
00363 }
00364
00365 double fraction = -1;
00366 if( RecHitsMuon != 0 ) fraction = shared/(double)RecHitsMuon;
00367
00368 if( shared > sharedHits_ && fraction > sharedFrac_ ) {
00369 overlappingMuon = true;
00370 break;
00371 }
00372 }
00373 }
00374
00375 return overlappingMuon;
00376 }
00377
00378
00379
00380
00381 unsigned int
00382 MuonCosmicCompatibilityFiller::pvMatches(const edm::Event& iEvent, const reco::Muon& muon, bool isLoose) const {
00383
00384 float maxdxyMult, maxdzMult, maxdxy, maxdz;
00385
00386 if (isLoose) {
00387
00388 maxdxyMult = maxdxyLooseMult_;
00389 maxdzMult = maxdzLooseMult_;
00390 maxdxy = maxdxyLoose_;
00391 maxdz = maxdzLoose_;
00392 } else {
00393 maxdxyMult = maxdxyTightMult_;
00394 maxdzMult = maxdzTightMult_;
00395 maxdxy = maxdxyTight_;
00396 maxdz = maxdzTight_;
00397 }
00398
00399 unsigned int result = 0;
00400
00401 reco::TrackRef track;
00402 if ( muon.isGlobalMuon() ) track = muon.innerTrack();
00403 else if ( muon.isTrackerMuon() ) track = muon.track();
00404 else if ( muon.isStandAloneMuon()) track = muon.standAloneMuon();
00405
00406 bool multipleMu = false;
00407 if (nMuons(iEvent) > 1) multipleMu = true;
00408
00409 math::XYZPoint RefVtx;
00410 RefVtx.SetXYZ(0, 0, 0);
00411
00412 edm::Handle<reco::VertexCollection> pvHandle;
00413 iEvent.getByLabel(inputVertexCollection_,pvHandle);
00414 const reco::VertexCollection & vertices = *pvHandle.product();
00415 for(reco::VertexCollection::const_iterator it=vertices.begin() ; it!=vertices.end() ; ++it){
00416 RefVtx = it->position();
00417
00418 if ( track.isNonnull() ) {
00419 if (multipleMu) {
00420
00421
00422 if ( fabs( (*track).dxy(RefVtx) ) < maxdxyMult || fabs( (*track).dz(RefVtx) ) < maxdzMult) {
00423 result++;
00424
00425
00426 if (!isLoose && fabs( (*track).dxy(RefVtx) ) > largedxyMult_) result -= 1;
00427
00428 }
00429 } else {
00430
00431
00432 if (fabs( (*track).dxy(RefVtx) ) < maxdxy || fabs( (*track).dz(RefVtx) ) < maxdz) {
00433 result++;
00434
00435
00436 if (!isLoose && fabs( (*track).dxy(RefVtx) ) > largedxy_) result -= 1;
00437
00438 }
00439 }
00440 }
00441 }
00442
00443
00444 if (result == 0 && multipleMu) {
00445
00446 edm::Handle<reco::MuonCollection> muonHandle;
00447 iEvent.getByLabel(inputMuonCollections_[1], muonHandle);
00448
00449
00450 edm::Handle<reco::VertexCollection> pvHandle;
00451 iEvent.getByLabel(inputVertexCollection_,pvHandle);
00452 const reco::VertexCollection & vertices = *pvHandle.product();
00453
00454
00455 if( !muonHandle.failedToGet() ) {
00456 for ( reco::MuonCollection::const_iterator muons = muonHandle->begin(); muons != muonHandle->end(); ++muons ) {
00457 if (!muons->isGlobalMuon()) continue;
00458
00459 if ( muons->innerTrack() == muon.innerTrack() && muons->outerTrack() == muon.outerTrack()) continue;
00460
00461 reco::TrackRef tracks;
00462 if (muons->isGlobalMuon()) tracks = muons->innerTrack();
00463 if (fabs((*tracks).dxy(RefVtx)) > hIpTrdxy_) continue;
00464
00465 if (vertices.begin() == vertices.end()) continue;
00466
00467
00468
00469 if (TMath::Prob(vertices.front().chi2(),(int)(vertices.front().ndof())) > hIpTrvProb_) result = 1;
00470
00471 }
00472 }
00473 }
00474
00475 return result;
00476
00477 }
00478
00479 float
00480 MuonCosmicCompatibilityFiller::combinedCosmicID(const edm::Event& iEvent,
00481 const edm::EventSetup& iSetup, const reco::Muon& muon, bool CheckMuonID, bool checkVertex) const {
00482
00483 float result = 0.0;
00484
00485
00486
00487 if( muon.isGlobalMuon() ) {
00488
00489 unsigned int cosmicVertex = eventActivity(iEvent, muon);
00490 bool isOverlapping = isOverlappingMuon(iEvent, iSetup, muon);
00491 unsigned int looseIp = pvMatches(iEvent, muon, true);
00492 unsigned int tightIp = pvMatches(iEvent, muon, false);
00493 float looseTime = muonTiming(iEvent, muon, true);
00494 float tightTime = muonTiming(iEvent, muon, false);
00495 unsigned int backToback = backToBack2LegCosmic(iEvent,muon);
00496
00497
00498
00499 if (checkVertex && cosmicVertex == 0) return 10.0;
00500
00501
00502
00503
00504
00505
00506 double weight_btob = 2.0;
00507 double weight_ip = 2.0;
00508 double weight_time = 1.0;
00509 double weight_overlap = 0.5;
00510
00511
00512
00513
00514
00515 if( backToback >= 1 ) {
00516
00517 result += weight_btob*2.;
00518 if( tightIp == 1 ) {
00519
00520 if( looseIp == 1 ) {
00521 if( backToback < 2 ) result -= weight_btob*0.5;
00522 }
00523 }
00524 }
00525
00526
00527 if( tightIp == 0 ) {
00528
00529 result += weight_ip*2.0;
00530 if( backToback == 0 ) {
00531
00532 if( tightTime == 0 ) {
00533 if( looseTime == 0 && !isOverlapping ) result -= weight_ip*1.0;
00534 }
00535 }
00536 }
00537 else if( tightIp >= 2 ) {
00538
00539
00540 if( backToback >= 1 ) result -= weight_ip*1.0;
00541 }
00542
00543
00544 if( tightTime > 0 ) {
00545
00546 if( looseTime > 0 ) {
00547 if( backToback >= 1 ) {
00548 if( tightIp == 0 ) result += weight_time*tightTime;
00549 else if( looseIp == 0 ) result += weight_time*0.25;
00550 }
00551 }
00552 else {
00553 if( backToback >= 1 && tightIp == 0 ) result += weight_time*0.25;
00554 }
00555 }
00556
00557
00558 if( backToback == 0 && isOverlapping ) {
00559
00560 if( tightIp == 0 && tightTime >= 1 ) {
00561 result += weight_overlap*1.0;
00562 }
00563 }
00564 }
00565
00566
00567
00568
00569 return result;
00570 }
00571
00572
00573
00574
00575
00576
00577 unsigned int
00578 MuonCosmicCompatibilityFiller::eventActivity(const edm::Event& iEvent, const reco::Muon& muon) const
00579 {
00580
00581 unsigned int result = 0;
00582
00583
00584 edm::Handle<reco::TrackCollection> tracks;
00585 iEvent.getByLabel(inputTrackCollections_[0],tracks);
00586 if (!tracks.failedToGet() && tracks->size() < 3) return 0;
00587
00588
00589 edm::Handle<reco::VertexCollection> pvHandle;
00590 if (!iEvent.getByLabel(inputVertexCollection_,pvHandle)) {return 0;} else {
00591 const reco::VertexCollection & vertices = *pvHandle.product();
00592
00593 if (vertices.begin() == vertices.end()) return 0;
00594 for(reco::VertexCollection::const_iterator it=vertices.begin() ; it!=vertices.end() ; ++it){
00595 if ((TMath::Prob(it->chi2(),(int)it->ndof()) > minvProb_) && (fabs(it->z()) <= maxvertZ_) && (fabs(it->position().rho()) <= maxvertRho_)) result++;
00596 }
00597 }
00598 return result;
00599 }
00600
00601
00602
00603
00604 bool MuonCosmicCompatibilityFiller::checkMuonID( const reco::Muon& imuon ) const
00605 {
00606 bool result = false;
00607
00608 if( muon::isGoodMuon(imuon, muon::GlobalMuonPromptTight) && muon::isGoodMuon(imuon, muon::TMOneStationLoose) ) result = true;
00609
00610 return result;
00611 }
00612
00613 bool MuonCosmicCompatibilityFiller::checkMuonSegments(const reco::Muon& imuon) const {
00614
00615 bool result = false;
00616 if (imuon.numberOfMatches() < nChamberMatches_ && muon::segmentCompatibility(imuon) < segmentComp_) result = true;
00617
00618 return result;
00619 }