File indexing completed on 2024-09-07 04:37:46
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 #include <memory>
0016 #include <string>
0017
0018
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023
0024 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0025 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0026 #include "DataFormats/DetId/interface/DetId.h"
0027
0028 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0029 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0030 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0031
0032 #include "RecoMuon/MuonIdentification/interface/MuonCosmicCompatibilityFiller.h"
0033 #include "RecoMuon/MuonIdentification/interface/MuonCosmicsId.h"
0034
0035 #include "DataFormats/MuonReco/interface/MuonCosmicCompatibility.h"
0036 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0037
0038 #include "TMath.h"
0039
0040 using namespace edm;
0041 using namespace std;
0042
0043 MuonCosmicCompatibilityFiller::MuonCosmicCompatibilityFiller(const edm::ParameterSet& iConfig,
0044 edm::ConsumesCollector& iC)
0045 : inputMuonCollections_(iConfig.getParameter<std::vector<edm::InputTag> >("InputMuonCollections")),
0046 inputTrackCollections_(iConfig.getParameter<std::vector<edm::InputTag> >("InputTrackCollections")),
0047 inputCosmicMuonCollection_(iConfig.getParameter<edm::InputTag>("InputCosmicMuonCollection")),
0048 inputVertexCollection_(iConfig.getParameter<edm::InputTag>("InputVertexCollection")) {
0049
0050 angleThreshold_ = iConfig.getParameter<double>("angleCut");
0051 deltaPt_ = iConfig.getParameter<double>("deltaPt");
0052
0053 offTimePosTightMult_ = iConfig.getParameter<double>("offTimePosTightMult");
0054 offTimeNegTightMult_ = iConfig.getParameter<double>("offTimeNegTightMult");
0055 offTimePosTight_ = iConfig.getParameter<double>("offTimePosTight");
0056 offTimeNegTight_ = iConfig.getParameter<double>("offTimeNegTight");
0057 offTimePosLooseMult_ = iConfig.getParameter<double>("offTimePosLooseMult");
0058 offTimeNegLooseMult_ = iConfig.getParameter<double>("offTimeNegLooseMult");
0059 offTimePosLoose_ = iConfig.getParameter<double>("offTimePosLoose");
0060 offTimeNegLoose_ = iConfig.getParameter<double>("offTimeNegLoose");
0061 corrTimeNeg_ = iConfig.getParameter<double>("corrTimeNeg");
0062 corrTimePos_ = iConfig.getParameter<double>("corrTimePos");
0063
0064 sharedHits_ = iConfig.getParameter<int>("sharedHits");
0065 sharedFrac_ = iConfig.getParameter<double>("sharedFrac");
0066 ipThreshold_ = iConfig.getParameter<double>("ipCut");
0067
0068 nChamberMatches_ = iConfig.getParameter<int>("nChamberMatches");
0069 segmentComp_ = iConfig.getParameter<double>("segmentComp");
0070
0071 maxdzLooseMult_ = iConfig.getParameter<double>("maxdzLooseMult");
0072 maxdxyLooseMult_ = iConfig.getParameter<double>("maxdxyLooseMult");
0073 maxdzTightMult_ = iConfig.getParameter<double>("maxdzTightMult");
0074 maxdxyTightMult_ = iConfig.getParameter<double>("maxdxyTightMult");
0075 maxdzLoose_ = iConfig.getParameter<double>("maxdzLoose");
0076 maxdxyLoose_ = iConfig.getParameter<double>("maxdxyLoose");
0077 maxdzTight_ = iConfig.getParameter<double>("maxdzTight");
0078 maxdxyTight_ = iConfig.getParameter<double>("maxdxyTight");
0079 largedxyMult_ = iConfig.getParameter<double>("largedxyMult");
0080 largedxy_ = iConfig.getParameter<double>("largedxy");
0081 hIpTrdxy_ = iConfig.getParameter<double>("hIpTrdxy");
0082 hIpTrvProb_ = iConfig.getParameter<double>("hIpTrvProb");
0083 minvProb_ = iConfig.getParameter<double>("minvProb");
0084 maxvertZ_ = iConfig.getParameter<double>("maxvertZ");
0085 maxvertRho_ = iConfig.getParameter<double>("maxvertRho");
0086
0087
0088 for (unsigned int i = 0; i < inputMuonCollections_.size(); ++i)
0089 muonTokens_.push_back(iC.consumes<reco::MuonCollection>(inputMuonCollections_.at(i)));
0090 for (unsigned int i = 0; i < inputTrackCollections_.size(); ++i)
0091 trackTokens_.push_back(iC.consumes<reco::TrackCollection>(inputTrackCollections_.at(i)));
0092
0093 cosmicToken_ = iC.consumes<reco::MuonCollection>(inputCosmicMuonCollection_);
0094 vertexToken_ = iC.consumes<reco::VertexCollection>(inputVertexCollection_);
0095 geometryToken_ = iC.esConsumes<GlobalTrackingGeometry, GlobalTrackingGeometryRecord>();
0096 }
0097
0098 MuonCosmicCompatibilityFiller::~MuonCosmicCompatibilityFiller() {}
0099
0100 reco::MuonCosmicCompatibility MuonCosmicCompatibilityFiller::fillCompatibility(const reco::Muon& muon,
0101 edm::Event& iEvent,
0102 const edm::EventSetup& iSetup) {
0103 const std::string theCategory = "MuonCosmicCompatibilityFiller";
0104
0105 reco::MuonCosmicCompatibility returnComp;
0106
0107 float timeCompatibility = muonTiming(iEvent, muon, false);
0108 float backToBackCompatibility = backToBack2LegCosmic(iEvent, muon);
0109 float overlapCompatibility = isOverlappingMuon(iEvent, iSetup, muon);
0110 float ipCompatibility = pvMatches(iEvent, muon, false);
0111 float vertexCompatibility = eventActivity(iEvent, muon);
0112 float combinedCompatibility = combinedCosmicID(iEvent, iSetup, muon, false, false);
0113
0114 returnComp.timeCompatibility = timeCompatibility;
0115 returnComp.backToBackCompatibility = backToBackCompatibility;
0116 returnComp.overlapCompatibility = overlapCompatibility;
0117 returnComp.cosmicCompatibility = combinedCompatibility;
0118 returnComp.ipCompatibility = ipCompatibility;
0119 returnComp.vertexCompatibility = vertexCompatibility;
0120
0121 return returnComp;
0122 }
0123
0124
0125
0126
0127 float MuonCosmicCompatibilityFiller::muonTiming(const edm::Event& iEvent, const reco::Muon& muon, bool isLoose) const {
0128 float offTimeNegMult, offTimePosMult, offTimeNeg, offTimePos;
0129
0130 if (isLoose) {
0131
0132 offTimeNegMult = offTimeNegLooseMult_;
0133 offTimePosMult = offTimePosLooseMult_;
0134 offTimeNeg = offTimeNegLoose_;
0135 offTimePos = offTimePosLoose_;
0136 } else {
0137 offTimeNegMult = offTimeNegTightMult_;
0138 offTimePosMult = offTimePosTightMult_;
0139 offTimeNeg = offTimeNegTight_;
0140 offTimePos = offTimePosTight_;
0141 }
0142
0143 float result = 0.0;
0144
0145 if (muon.isTimeValid()) {
0146
0147 if (nMuons(iEvent) > 1) {
0148 float positiveTime = 0;
0149 if (muon.time().timeAtIpInOut < offTimeNegMult || muon.time().timeAtIpInOut > offTimePosMult)
0150 result = 1.;
0151 if (muon.time().timeAtIpInOut > 0.)
0152 positiveTime = muon.time().timeAtIpInOut;
0153
0154
0155
0156 if (!isLoose && result == 0 && positiveTime > corrTimePos_) {
0157
0158 bool isUp = false;
0159 reco::TrackRef outertrack = muon.outerTrack();
0160 if (outertrack.isNonnull()) {
0161 if (outertrack->phi() > 0)
0162 isUp = true;
0163
0164
0165 edm::Handle<reco::MuonCollection> muonHandle;
0166 iEvent.getByToken(muonTokens_[1], muonHandle);
0167
0168 if (!muonHandle.failedToGet()) {
0169 for (reco::MuonCollection::const_iterator iMuon = muonHandle->begin(); iMuon != muonHandle->end();
0170 ++iMuon) {
0171 if (!iMuon->isGlobalMuon())
0172 continue;
0173
0174 reco::TrackRef checkedTrack = iMuon->outerTrack();
0175 if (muon.isTimeValid()) {
0176
0177 if (checkedTrack->phi() < 0 && isUp) {
0178 if (iMuon->time().timeAtIpInOut < corrTimeNeg_)
0179 result = 1.0;
0180 break;
0181 } else if (checkedTrack->phi() > 0 && !isUp) {
0182
0183 if (iMuon->time().timeAtIpInOut < corrTimeNeg_)
0184 result = 1.0;
0185 break;
0186 }
0187 }
0188 }
0189 }
0190 }
0191 }
0192 } else {
0193
0194 if (muon.time().timeAtIpInOut < offTimeNeg || muon.time().timeAtIpInOut > offTimePos)
0195 result = 1.;
0196 }
0197 }
0198
0199 if (!isLoose && result > 0) {
0200
0201 if (pvMatches(iEvent, muon, true) == 0)
0202 result *= 2.;
0203 }
0204
0205 return result;
0206 }
0207
0208
0209
0210
0211 unsigned int MuonCosmicCompatibilityFiller::backToBack2LegCosmic(const edm::Event& iEvent,
0212 const reco::Muon& muon) const {
0213 unsigned int result = 0;
0214 reco::TrackRef track;
0215 if (muon.isGlobalMuon())
0216 track = muon.innerTrack();
0217 else if (muon.isTrackerMuon())
0218 track = muon.track();
0219 else if (muon.isStandAloneMuon() || muon.isRPCMuon() || muon.isGEMMuon() || muon.isME0Muon())
0220 return false;
0221
0222 for (unsigned int iColl = 0; iColl < trackTokens_.size(); ++iColl) {
0223 edm::Handle<reco::TrackCollection> trackHandle;
0224 iEvent.getByToken(trackTokens_[iColl], trackHandle);
0225 if (muonid::findOppositeTrack(trackHandle, *track, angleThreshold_, deltaPt_).isNonnull()) {
0226 result++;
0227 }
0228 }
0229
0230 return result;
0231 }
0232
0233
0234
0235
0236 unsigned int MuonCosmicCompatibilityFiller::nMuons(const edm::Event& iEvent) const {
0237 unsigned int nGlb = 0;
0238
0239 edm::Handle<reco::MuonCollection> muonHandle;
0240 iEvent.getByToken(muonTokens_[1], muonHandle);
0241
0242 if (!muonHandle.failedToGet()) {
0243 for (reco::MuonCollection::const_iterator iMuon = muonHandle->begin(); iMuon != muonHandle->end(); ++iMuon) {
0244 if (!iMuon->isGlobalMuon())
0245 continue;
0246 nGlb++;
0247 }
0248 }
0249
0250 return nGlb;
0251 }
0252
0253
0254
0255
0256 bool MuonCosmicCompatibilityFiller::isOverlappingMuon(const edm::Event& iEvent,
0257 const edm::EventSetup& iSetup,
0258 const reco::Muon& muon) const {
0259
0260
0261
0262
0263
0264
0265
0266 bool overlappingMuon = false;
0267 if (!muon.isGlobalMuon())
0268 return false;
0269
0270
0271 edm::Handle<reco::MuonCollection> muonHandle;
0272 iEvent.getByToken(cosmicToken_, muonHandle);
0273
0274
0275 ESHandle<GlobalTrackingGeometry> trackingGeometry = iSetup.getHandle(geometryToken_);
0276
0277
0278 math::XYZPoint RefVtx;
0279 RefVtx.SetXYZ(0, 0, 0);
0280
0281 edm::Handle<reco::VertexCollection> pvHandle;
0282 iEvent.getByToken(vertexToken_, pvHandle);
0283 const reco::VertexCollection& vertices = *pvHandle.product();
0284 for (reco::VertexCollection::const_iterator it = vertices.begin(); it != vertices.end(); ++it) {
0285 RefVtx = it->position();
0286 }
0287
0288 if (!muonHandle.failedToGet()) {
0289 for (reco::MuonCollection::const_iterator cosmicMuon = muonHandle->begin(); cosmicMuon != muonHandle->end();
0290 ++cosmicMuon) {
0291 if (cosmicMuon->innerTrack() == muon.innerTrack() || cosmicMuon->outerTrack() == muon.outerTrack())
0292 return true;
0293
0294 reco::TrackRef outertrack = muon.outerTrack();
0295 reco::TrackRef costrack = cosmicMuon->outerTrack();
0296
0297
0298 int RecHitsMuon = outertrack->numberOfValidHits();
0299 int shared = 0;
0300 if (costrack.isNonnull()) {
0301
0302
0303
0304
0305
0306 if (outertrack.isNonnull()) {
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320 for (trackingRecHit_iterator trkhit = outertrack->recHitsBegin(); trkhit != outertrack->recHitsEnd();
0321 trkhit++) {
0322 if ((*trkhit)->isValid()) {
0323 for (trackingRecHit_iterator coshit = costrack->recHitsBegin(); coshit != costrack->recHitsEnd();
0324 coshit++) {
0325 if ((*coshit)->isValid()) {
0326 if ((*trkhit)->geographicalId() == (*coshit)->geographicalId()) {
0327 if (((*trkhit)->localPosition() - (*coshit)->localPosition()).mag() < 10e-5)
0328 shared++;
0329 }
0330 }
0331 }
0332 }
0333 }
0334 }
0335 }
0336
0337 double fraction = -1;
0338 if (RecHitsMuon != 0)
0339 fraction = shared / (double)RecHitsMuon;
0340 if (shared > sharedHits_ && fraction > sharedFrac_) {
0341 overlappingMuon = true;
0342 break;
0343 }
0344 }
0345 }
0346
0347 return overlappingMuon;
0348 }
0349
0350
0351
0352
0353 unsigned int MuonCosmicCompatibilityFiller::pvMatches(const edm::Event& iEvent,
0354 const reco::Muon& muon,
0355 bool isLoose) const {
0356 float maxdxyMult, maxdzMult, maxdxy, maxdz;
0357
0358 if (isLoose) {
0359
0360 maxdxyMult = maxdxyLooseMult_;
0361 maxdzMult = maxdzLooseMult_;
0362 maxdxy = maxdxyLoose_;
0363 maxdz = maxdzLoose_;
0364 } else {
0365 maxdxyMult = maxdxyTightMult_;
0366 maxdzMult = maxdzTightMult_;
0367 maxdxy = maxdxyTight_;
0368 maxdz = maxdzTight_;
0369 }
0370
0371 unsigned int result = 0;
0372
0373 reco::TrackRef track;
0374 if (muon.isGlobalMuon())
0375 track = muon.innerTrack();
0376 else if (muon.isTrackerMuon() || muon.isRPCMuon())
0377 track = muon.track();
0378 else if (muon.isStandAloneMuon())
0379 track = muon.standAloneMuon();
0380
0381 bool multipleMu = false;
0382 if (nMuons(iEvent) > 1)
0383 multipleMu = true;
0384
0385 math::XYZPoint RefVtx;
0386 RefVtx.SetXYZ(0, 0, 0);
0387
0388 edm::Handle<reco::VertexCollection> pvHandle;
0389 iEvent.getByToken(vertexToken_, pvHandle);
0390 const reco::VertexCollection& vertices = *pvHandle.product();
0391 for (reco::VertexCollection::const_iterator it = vertices.begin(); it != vertices.end(); ++it) {
0392 RefVtx = it->position();
0393
0394 if (track.isNonnull()) {
0395 if (multipleMu) {
0396
0397
0398 if (fabs((*track).dxy(RefVtx)) < maxdxyMult || fabs((*track).dz(RefVtx)) < maxdzMult) {
0399 result++;
0400
0401
0402 if (!isLoose && fabs((*track).dxy(RefVtx)) > largedxyMult_)
0403 result -= 1;
0404 }
0405 } else {
0406
0407
0408 if (fabs((*track).dxy(RefVtx)) < maxdxy || fabs((*track).dz(RefVtx)) < maxdz) {
0409 result++;
0410
0411
0412 if (!isLoose && fabs((*track).dxy(RefVtx)) > largedxy_)
0413 result -= 1;
0414 }
0415 }
0416 }
0417 }
0418
0419
0420 if (result == 0 && multipleMu) {
0421
0422 edm::Handle<reco::MuonCollection> muonHandle;
0423 iEvent.getByToken(muonTokens_[1], muonHandle);
0424
0425
0426 edm::Handle<reco::VertexCollection> pvHandle;
0427 iEvent.getByToken(vertexToken_, pvHandle);
0428 const reco::VertexCollection& vertices = *pvHandle.product();
0429
0430
0431 if (!muonHandle.failedToGet()) {
0432 for (reco::MuonCollection::const_iterator muons = muonHandle->begin(); muons != muonHandle->end(); ++muons) {
0433 if (!muons->isGlobalMuon())
0434 continue;
0435
0436 if (muons->innerTrack() == muon.innerTrack() && muons->outerTrack() == muon.outerTrack())
0437 continue;
0438
0439 reco::TrackRef tracks;
0440 if (muons->isGlobalMuon())
0441 tracks = muons->innerTrack();
0442 if (fabs((*tracks).dxy(RefVtx)) > hIpTrdxy_)
0443 continue;
0444
0445 if (vertices.begin() == vertices.end())
0446 continue;
0447
0448
0449
0450 if (TMath::Prob(vertices.front().chi2(), (int)(vertices.front().ndof())) > hIpTrvProb_)
0451 result = 1;
0452
0453 }
0454 }
0455 }
0456
0457 return result;
0458 }
0459
0460 float MuonCosmicCompatibilityFiller::combinedCosmicID(const edm::Event& iEvent,
0461 const edm::EventSetup& iSetup,
0462 const reco::Muon& muon,
0463 bool CheckMuonID,
0464 bool checkVertex) const {
0465 float result = 0.0;
0466
0467
0468
0469 if (muon.isGlobalMuon()) {
0470 unsigned int cosmicVertex = eventActivity(iEvent, muon);
0471 bool isOverlapping = isOverlappingMuon(iEvent, iSetup, muon);
0472 unsigned int looseIp = pvMatches(iEvent, muon, true);
0473 unsigned int tightIp = pvMatches(iEvent, muon, false);
0474 float looseTime = muonTiming(iEvent, muon, true);
0475 float tightTime = muonTiming(iEvent, muon, false);
0476 unsigned int backToback = backToBack2LegCosmic(iEvent, muon);
0477
0478
0479
0480 if (checkVertex && cosmicVertex == 0)
0481 return 10.0;
0482
0483
0484
0485
0486
0487
0488 double weight_btob = 2.0;
0489 double weight_ip = 2.0;
0490 double weight_time = 1.0;
0491 double weight_overlap = 0.5;
0492
0493
0494
0495
0496
0497 if (backToback >= 1) {
0498
0499 result += weight_btob * 2.;
0500 if (tightIp == 1) {
0501
0502 if (looseIp == 1) {
0503 if (backToback < 2)
0504 result -= weight_btob * 0.5;
0505 }
0506 }
0507 }
0508
0509
0510 if (tightIp == 0) {
0511
0512 result += weight_ip * 2.0;
0513 if (backToback == 0) {
0514
0515 if (tightTime == 0) {
0516 if (looseTime == 0 && !isOverlapping)
0517 result -= weight_ip * 1.0;
0518 }
0519 }
0520 } else if (tightIp >= 2) {
0521
0522
0523 if (backToback >= 1)
0524 result -= weight_ip * 1.0;
0525 }
0526
0527
0528 if (tightTime > 0) {
0529
0530 if (looseTime > 0) {
0531 if (backToback >= 1) {
0532 if (tightIp == 0)
0533 result += weight_time * tightTime;
0534 else if (looseIp == 0)
0535 result += weight_time * 0.25;
0536 }
0537 } else {
0538 if (backToback >= 1 && tightIp == 0)
0539 result += weight_time * 0.25;
0540 }
0541 }
0542
0543
0544 if (backToback == 0 && isOverlapping) {
0545
0546 if (tightIp == 0 && tightTime >= 1) {
0547 result += weight_overlap * 1.0;
0548 }
0549 }
0550 }
0551
0552
0553
0554 return result;
0555 }
0556
0557
0558
0559
0560 unsigned int MuonCosmicCompatibilityFiller::eventActivity(const edm::Event& iEvent, const reco::Muon& muon) const {
0561 unsigned int result = 0;
0562
0563
0564 edm::Handle<reco::TrackCollection> tracks;
0565 iEvent.getByToken(trackTokens_[0], tracks);
0566 if (!tracks.failedToGet() && tracks->size() < 3)
0567 return 0;
0568
0569
0570 edm::Handle<reco::VertexCollection> pvHandle;
0571 if (!iEvent.getByToken(vertexToken_, pvHandle)) {
0572 return 0;
0573 } else {
0574 const reco::VertexCollection& vertices = *pvHandle.product();
0575
0576 if (vertices.begin() == vertices.end())
0577 return 0;
0578 for (reco::VertexCollection::const_iterator it = vertices.begin(); it != vertices.end(); ++it) {
0579 if ((TMath::Prob(it->chi2(), (int)it->ndof()) > minvProb_) && (fabs(it->z()) <= maxvertZ_) &&
0580 (fabs(it->position().rho()) <= maxvertRho_))
0581 result++;
0582 }
0583 }
0584 return result;
0585 }
0586
0587
0588
0589
0590 bool MuonCosmicCompatibilityFiller::checkMuonID(const reco::Muon& imuon) const {
0591 bool result = false;
0592
0593 if (muon::isGoodMuon(imuon, muon::GlobalMuonPromptTight) && muon::isGoodMuon(imuon, muon::TMOneStationLoose))
0594 result = true;
0595
0596 return result;
0597 }
0598
0599 bool MuonCosmicCompatibilityFiller::checkMuonSegments(const reco::Muon& imuon) const {
0600 bool result = false;
0601 if (imuon.numberOfMatches() < nChamberMatches_ && muon::segmentCompatibility(imuon) < segmentComp_)
0602 result = true;
0603
0604 return result;
0605 }