Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:25:57

0001 /**
0002  *  \package: MuonIdentification
0003  *  \class: MuonCosmicCompatibilityFiller
0004  *
0005  *  Description: class for cosmic muon identification
0006  *
0007  *
0008  *  \author: A. Everett, Purdue University
0009  *  \author: A. Svyatkovskiy, Purdue University
0010  *  \author: H.D. Yoo, Purdue University
0011  *
0012  **/
0013 
0014 // system include files
0015 #include <memory>
0016 #include <string>
0017 
0018 // user include files
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   //kinematic vars
0050   angleThreshold_ = iConfig.getParameter<double>("angleCut");
0051   deltaPt_ = iConfig.getParameter<double>("deltaPt");
0052   //time
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   //rechits
0064   sharedHits_ = iConfig.getParameter<int>("sharedHits");
0065   sharedFrac_ = iConfig.getParameter<double>("sharedFrac");
0066   ipThreshold_ = iConfig.getParameter<double>("ipCut");
0067   //segment comp, matches
0068   nChamberMatches_ = iConfig.getParameter<int>("nChamberMatches");
0069   segmentComp_ = iConfig.getParameter<double>("segmentComp");
0070   //ip, vertex
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   //  nTrackThreshold_ = iConfig.getParameter<unsigned int>("nTrackThreshold");
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 //Timing: 0 - not cosmic-like
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     //use "loose" parameter set
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     //case of multiple muon event
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       //special case, looking for time-correlation
0155       // between muons in opposite hemispheres
0156       if (!isLoose && result == 0 && positiveTime > corrTimePos_) {
0157         //check hemi of this muon
0158         bool isUp = false;
0159         reco::TrackRef outertrack = muon.outerTrack();
0160         if (outertrack.isNonnull()) {
0161           if (outertrack->phi() > 0)
0162             isUp = true;
0163 
0164           //loop over muons in that event and find if there are any in the opposite hemi
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                 // from bottom up
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                   // from top down
0183                   if (iMuon->time().timeAtIpInOut < corrTimeNeg_)
0184                     result = 1.0;
0185                   break;
0186                 }
0187               }  //muon is time valid
0188             }
0189           }
0190         }  //track is nonnull
0191       }    //double check timing
0192     } else {
0193       //case of a single muon event
0194       if (muon.time().timeAtIpInOut < offTimeNeg || muon.time().timeAtIpInOut > offTimePos)
0195         result = 1.;
0196     }
0197   }  //is time valid
0198 
0199   if (!isLoose && result > 0) {
0200     //check loose ip
0201     if (pvMatches(iEvent, muon, true) == 0)
0202       result *= 2.;
0203   }
0204 
0205   return result;
0206 }
0207 
0208 //
0209 //Back-to-back selector
0210 //
0211 unsigned int MuonCosmicCompatibilityFiller::backToBack2LegCosmic(const edm::Event& iEvent,
0212                                                                  const reco::Muon& muon) const {
0213   unsigned int result = 0;  //no partners - collision
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   }  //loop over track collections
0229 
0230   return result;
0231 }
0232 
0233 //
0234 //Check the number of global muons in an event, return true if there are more than 1 muon
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 //Check overlap between collections, use shared hits info
0255 //
0256 bool MuonCosmicCompatibilityFiller::isOverlappingMuon(const edm::Event& iEvent,
0257                                                       const edm::EventSetup& iSetup,
0258                                                       const reco::Muon& muon) const {
0259   // 4 steps in this module
0260   // step1 : check whether it's 1leg cosmic muon or not
0261   // step2 : both muons (muons and muonsFromCosmics1Leg) should have close IP
0262   // step3 : both muons should share very close reference point
0263   // step4 : check shared hits in both muon tracks
0264 
0265   // check if this muon is available in muonsFromCosmics collection
0266   bool overlappingMuon = false;  //false - not cosmic-like
0267   if (!muon.isGlobalMuon())
0268     return false;
0269 
0270   // reco muons for cosmics
0271   edm::Handle<reco::MuonCollection> muonHandle;
0272   iEvent.getByToken(cosmicToken_, muonHandle);
0273 
0274   // Global Tracking Geometry
0275   ESHandle<GlobalTrackingGeometry> trackingGeometry = iSetup.getHandle(geometryToken_);
0276 
0277   // PV
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       bool isUp = false;
0298       if (outertrack->phi() > 0)
0299         isUp = true;
0300 
0301       // shared hits
0302       int RecHitsMuon = outertrack->numberOfValidHits();
0303       int RecHitsCosmicMuon = 0;
0304       int shared = 0;
0305       // count hits for same hemisphere
0306       if (costrack.isNonnull()) {
0307         int nhitsUp = 0;
0308         int nhitsDown = 0;
0309         // unused
0310         //  bool isCosmic1Leg = false;
0311         //  bool isCloseIP = false;
0312         //  bool isCloseRef = false;
0313 
0314         for (trackingRecHit_iterator coshit = costrack->recHitsBegin(); coshit != costrack->recHitsEnd(); coshit++) {
0315           if ((*coshit)->isValid()) {
0316             DetId id((*coshit)->geographicalId());
0317             double hity = trackingGeometry->idToDet(id)->position().y();
0318             if (hity > 0)
0319               nhitsUp++;
0320             if (hity < 0)
0321               nhitsDown++;
0322 
0323             if (isUp && hity > 0)
0324               RecHitsCosmicMuon++;
0325             if (!isUp && hity < 0)
0326               RecHitsCosmicMuon++;
0327           }
0328         }
0329         // step1
0330         //UNUSED:   if( nhitsUp > 0 && nhitsDown > 0 ) isCosmic1Leg = true;
0331         //if( !isCosmic1Leg ) continue;
0332 
0333         if (outertrack.isNonnull()) {
0334           // step2
0335           //UNUSED:          const double ipErr = (double)outertrack->d0Error();
0336           //UNUSED:          double ipThreshold  = max(ipThreshold_, ipErr);
0337           //UNUSED:   if( fabs(outertrack->dxy(RefVtx) + costrack->dxy(RefVtx)) < ipThreshold ) isCloseIP = true;
0338           //if( !isCloseIP ) continue;
0339 
0340           // step3
0341           GlobalPoint muonRefVtx(outertrack->vx(), outertrack->vy(), outertrack->vz());
0342           GlobalPoint cosmicRefVtx(costrack->vx(), costrack->vy(), costrack->vz());
0343           //UNUSED:   float dist = (muonRefVtx - cosmicRefVtx).mag();
0344           //UNUSED:   if( dist < 0.1 ) isCloseRef = true;
0345           //if( !isCloseRef ) continue;
0346 
0347           for (trackingRecHit_iterator trkhit = outertrack->recHitsBegin(); trkhit != outertrack->recHitsEnd();
0348                trkhit++) {
0349             if ((*trkhit)->isValid()) {
0350               for (trackingRecHit_iterator coshit = costrack->recHitsBegin(); coshit != costrack->recHitsEnd();
0351                    coshit++) {
0352                 if ((*coshit)->isValid()) {
0353                   if ((*trkhit)->geographicalId() == (*coshit)->geographicalId()) {
0354                     if (((*trkhit)->localPosition() - (*coshit)->localPosition()).mag() < 10e-5)
0355                       shared++;
0356                   }
0357                 }
0358               }
0359             }
0360           }
0361         }
0362       }
0363       // step4
0364       double fraction = -1;
0365       if (RecHitsMuon != 0)
0366         fraction = shared / (double)RecHitsMuon;
0367       //         std::cout << "shared = " << shared << " " << fraction << " " << RecHitsMuon << " " << RecHitsCosmicMuon << std::endl;
0368       if (shared > sharedHits_ && fraction > sharedFrac_) {
0369         overlappingMuon = true;
0370         break;
0371       }
0372     }
0373   }
0374 
0375   return overlappingMuon;
0376 }
0377 
0378 //
0379 //pv matches
0380 //
0381 unsigned int MuonCosmicCompatibilityFiller::pvMatches(const edm::Event& iEvent,
0382                                                       const reco::Muon& muon,
0383                                                       bool isLoose) const {
0384   float maxdxyMult, maxdzMult, maxdxy, maxdz;
0385 
0386   if (isLoose) {
0387     //use "loose" parameter set
0388     maxdxyMult = maxdxyLooseMult_;
0389     maxdzMult = maxdzLooseMult_;
0390     maxdxy = maxdxyLoose_;
0391     maxdz = maxdzLoose_;
0392   } else {
0393     maxdxyMult = maxdxyTightMult_;
0394     maxdzMult = maxdzTightMult_;
0395     maxdxy = maxdxyTight_;
0396     maxdz = maxdzTight_;
0397   }
0398 
0399   unsigned int result = 0;
0400 
0401   reco::TrackRef track;
0402   if (muon.isGlobalMuon())
0403     track = muon.innerTrack();
0404   else if (muon.isTrackerMuon() || muon.isRPCMuon())
0405     track = muon.track();
0406   else if (muon.isStandAloneMuon())
0407     track = muon.standAloneMuon();
0408 
0409   bool multipleMu = false;
0410   if (nMuons(iEvent) > 1)
0411     multipleMu = true;
0412 
0413   math::XYZPoint RefVtx;
0414   RefVtx.SetXYZ(0, 0, 0);
0415 
0416   edm::Handle<reco::VertexCollection> pvHandle;
0417   iEvent.getByToken(vertexToken_, pvHandle);
0418   const reco::VertexCollection& vertices = *pvHandle.product();
0419   for (reco::VertexCollection::const_iterator it = vertices.begin(); it != vertices.end(); ++it) {
0420     RefVtx = it->position();
0421 
0422     if (track.isNonnull()) {
0423       if (multipleMu) {
0424         //multiple muon event
0425 
0426         if (fabs((*track).dxy(RefVtx)) < maxdxyMult || fabs((*track).dz(RefVtx)) < maxdzMult) {
0427           result++;
0428 
0429           //case of extra large dxy
0430           if (!isLoose && fabs((*track).dxy(RefVtx)) > largedxyMult_)
0431             result -= 1;
0432         }
0433       } else {
0434         //single muon event
0435 
0436         if (fabs((*track).dxy(RefVtx)) < maxdxy || fabs((*track).dz(RefVtx)) < maxdz) {
0437           result++;
0438 
0439           //case of extra large dxy
0440           if (!isLoose && fabs((*track).dxy(RefVtx)) > largedxy_)
0441             result -= 1;
0442         }
0443       }
0444     }  //track is nonnull
0445   }    //loop over vertices
0446 
0447   //special case for non-cosmic large ip muons
0448   if (result == 0 && multipleMu) {
0449     // consider all reco muons in an event
0450     edm::Handle<reco::MuonCollection> muonHandle;
0451     iEvent.getByToken(muonTokens_[1], muonHandle);
0452 
0453     //cosmic event should have zero good vertices
0454     edm::Handle<reco::VertexCollection> pvHandle;
0455     iEvent.getByToken(vertexToken_, pvHandle);
0456     const reco::VertexCollection& vertices = *pvHandle.product();
0457 
0458     //find the "other" one
0459     if (!muonHandle.failedToGet()) {
0460       for (reco::MuonCollection::const_iterator muons = muonHandle->begin(); muons != muonHandle->end(); ++muons) {
0461         if (!muons->isGlobalMuon())
0462           continue;
0463         //skip this track
0464         if (muons->innerTrack() == muon.innerTrack() && muons->outerTrack() == muon.outerTrack())
0465           continue;
0466         //check ip and vertex of the "other" muon
0467         reco::TrackRef tracks;
0468         if (muons->isGlobalMuon())
0469           tracks = muons->innerTrack();
0470         if (fabs((*tracks).dxy(RefVtx)) > hIpTrdxy_)
0471           continue;
0472         //check if vertex collection is empty
0473         if (vertices.begin() == vertices.end())
0474           continue;
0475         //for(reco::VertexCollection::const_iterator it=vertices.begin() ; it!=vertices.end() ; ++it) {
0476         //find matching vertex by position
0477         //if (fabs(it->z() - tracks->vz()) > 0.01) continue; //means will not be untagged from cosmics
0478         if (TMath::Prob(vertices.front().chi2(), (int)(vertices.front().ndof())) > hIpTrvProb_)
0479           result = 1;
0480         //}
0481       }
0482     }
0483   }
0484 
0485   return result;
0486 }
0487 
0488 float MuonCosmicCompatibilityFiller::combinedCosmicID(const edm::Event& iEvent,
0489                                                       const edm::EventSetup& iSetup,
0490                                                       const reco::Muon& muon,
0491                                                       bool CheckMuonID,
0492                                                       bool checkVertex) const {
0493   float result = 0.0;
0494 
0495   // return >=1 = identify as cosmic muon (the more like cosmics, the higher is the number)
0496   // return 0.0 = identify as collision muon
0497   if (muon.isGlobalMuon()) {
0498     unsigned int cosmicVertex = eventActivity(iEvent, muon);
0499     bool isOverlapping = isOverlappingMuon(iEvent, iSetup, muon);
0500     unsigned int looseIp = pvMatches(iEvent, muon, true);
0501     unsigned int tightIp = pvMatches(iEvent, muon, false);
0502     float looseTime = muonTiming(iEvent, muon, true);
0503     float tightTime = muonTiming(iEvent, muon, false);
0504     unsigned int backToback = backToBack2LegCosmic(iEvent, muon);
0505     //bool cosmicSegment = checkMuonSegments(muon);
0506 
0507     //short cut to reject cosmic event
0508     if (checkVertex && cosmicVertex == 0)
0509       return 10.0;
0510 
0511     // compatibility (0 - 10)
0512     // weight is assigned by the performance of individual module
0513     // btob: ~90% eff / ~0% misid
0514     // ip: ~90% eff / ~0% misid
0515     // time: ~30% eff / ~0% misid
0516     double weight_btob = 2.0;
0517     double weight_ip = 2.0;
0518     double weight_time = 1.0;
0519     double weight_overlap = 0.5;
0520 
0521     // collision muon should have compatibility < 4 (0 - 4)
0522     // cosmic muon should have compatibility >= 4 (4 - 10)
0523 
0524     // b-to-b (max comp.: 4.0)
0525     if (backToback >= 1) {
0526       //in this case it is cosmic for sure
0527       result += weight_btob * 2.;
0528       if (tightIp == 1) {
0529         // check with other observables to reduce mis-id (subtract compatibilities)
0530         if (looseIp == 1) {
0531           if (backToback < 2)
0532             result -= weight_btob * 0.5;
0533         }
0534       }
0535     }
0536 
0537     // ip (max comp.: 4.0)
0538     if (tightIp == 0) {
0539       //in this case it is cosmic for sure
0540       result += weight_ip * 2.0;
0541       if (backToback == 0) {
0542         // check with other observables to reduce mis-id (subtract compatibilities)
0543         if (tightTime == 0) {
0544           if (looseTime == 0 && !isOverlapping)
0545             result -= weight_ip * 1.0;
0546         }
0547       }
0548     } else if (tightIp >= 2) {
0549       // in this case it is almost collision-like (reduce compatibility)
0550       // if multi pvs: comp = -2.0
0551       if (backToback >= 1)
0552         result -= weight_ip * 1.0;
0553     }
0554 
0555     // timing (max comp.: 2.0)
0556     if (tightTime > 0) {
0557       // bonus track
0558       if (looseTime > 0) {
0559         if (backToback >= 1) {
0560           if (tightIp == 0)
0561             result += weight_time * tightTime;
0562           else if (looseIp == 0)
0563             result += weight_time * 0.25;
0564         }
0565       } else {
0566         if (backToback >= 1 && tightIp == 0)
0567           result += weight_time * 0.25;
0568       }
0569     }
0570 
0571     // overlapping
0572     if (backToback == 0 && isOverlapping) {
0573       // bonus track
0574       if (tightIp == 0 && tightTime >= 1) {
0575         result += weight_overlap * 1.0;
0576       }
0577     }
0578   }  //is global muon
0579 
0580   //  if (CheckMuonID && cosmicSegment) result += 4;
0581 
0582   return result;
0583 }
0584 
0585 //
0586 //Track activity/vertex quality, count good vertices
0587 //
0588 unsigned int MuonCosmicCompatibilityFiller::eventActivity(const edm::Event& iEvent, const reco::Muon& muon) const {
0589   unsigned int result = 0;  //no good vertices - cosmic-like
0590 
0591   //check track activity
0592   edm::Handle<reco::TrackCollection> tracks;
0593   iEvent.getByToken(trackTokens_[0], tracks);
0594   if (!tracks.failedToGet() && tracks->size() < 3)
0595     return 0;
0596 
0597   //cosmic event should have zero good vertices
0598   edm::Handle<reco::VertexCollection> pvHandle;
0599   if (!iEvent.getByToken(vertexToken_, pvHandle)) {
0600     return 0;
0601   } else {
0602     const reco::VertexCollection& vertices = *pvHandle.product();
0603     //check if vertex collection is empty
0604     if (vertices.begin() == vertices.end())
0605       return 0;
0606     for (reco::VertexCollection::const_iterator it = vertices.begin(); it != vertices.end(); ++it) {
0607       if ((TMath::Prob(it->chi2(), (int)it->ndof()) > minvProb_) && (fabs(it->z()) <= maxvertZ_) &&
0608           (fabs(it->position().rho()) <= maxvertRho_))
0609         result++;
0610     }
0611   }
0612   return result;
0613 }
0614 
0615 //
0616 //Muon iD variables
0617 //
0618 bool MuonCosmicCompatibilityFiller::checkMuonID(const reco::Muon& imuon) const {
0619   bool result = false;
0620   // initial set up using Jordan's study: GlobalMuonPromptTight + TMOneStationLoose
0621   if (muon::isGoodMuon(imuon, muon::GlobalMuonPromptTight) && muon::isGoodMuon(imuon, muon::TMOneStationLoose))
0622     result = true;
0623 
0624   return result;
0625 }
0626 
0627 bool MuonCosmicCompatibilityFiller::checkMuonSegments(const reco::Muon& imuon) const {
0628   bool result = false;
0629   if (imuon.numberOfMatches() < nChamberMatches_ && muon::segmentCompatibility(imuon) < segmentComp_)
0630     result = true;
0631 
0632   return result;
0633 }