CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/RecoMuon/MuonIdentification/src/MuonCosmicCompatibilityFiller.cc

Go to the documentation of this file.
00001 
00016 // system include files
00017 #include <memory>
00018 #include <string>
00019 
00020 // user include files
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   // service parameters
00063   edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
00064   service_ = new MuonServiceProxy(serviceParameters);     
00065   
00066   //kinematic vars
00067   angleThreshold_ = iConfig.getParameter<double>("angleCut");
00068   deltaPt_ = iConfig.getParameter<double>("deltaPt");
00069   //time
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   //rechits
00081   sharedHits_ = iConfig.getParameter<int>("sharedHits");
00082   sharedFrac_ = iConfig.getParameter<double>("sharedFrac");
00083   ipThreshold_ = iConfig.getParameter<double>("ipCut");
00084   //segment comp, matches
00085   nChamberMatches_ = iConfig.getParameter<int>("nChamberMatches");
00086   segmentComp_ = iConfig.getParameter<double>("segmentComp");
00087   //ip, vertex
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 //  nTrackThreshold_ = iConfig.getParameter<unsigned int>("nTrackThreshold");
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 //Timing: 0 - not cosmic-like
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      //use "loose" parameter set
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     //case of multiple muon event
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       //special case, looking for time-correlation 
00169       // between muons in opposite hemispheres 
00170       if (!isLoose && result == 0 && positiveTime > corrTimePos_) { 
00171 
00172           //check hemi of this muon
00173           bool isUp = false;
00174           reco::TrackRef outertrack = muon.outerTrack();
00175           if( outertrack.isNonnull() ) {
00176               if( outertrack->phi() > 0 ) isUp = true;
00177 
00178               //loop over muons in that event and find if there are any in the opposite hemi 
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                           // from bottom up
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                                // from top down 
00195                                if (iMuon->time().timeAtIpInOut < corrTimeNeg_) result = 1.0; 
00196                                break;
00197                              }
00198                          } //muon is time valid
00199                      } 
00200                  }
00201              } //track is nonnull
00202          } //double check timing
00203      } else {
00204         //case of a single muon event
00205         if ( muon.time().timeAtIpInOut < offTimeNeg || muon.time().timeAtIpInOut > offTimePos) result = 1.;
00206       }
00207     } //is time valid
00208 
00209   
00210    if (!isLoose && result > 0) {
00211      //check loose ip
00212      if (pvMatches(iEvent, muon, true) == 0) result *= 2.;
00213    } 
00214 
00215    return result;
00216 }
00217 
00218 //
00219 //Back-to-back selector 
00220 //
00221 unsigned int 
00222 MuonCosmicCompatibilityFiller::backToBack2LegCosmic(const edm::Event& iEvent, const reco::Muon& muon) const {
00223 
00224   unsigned int result = 0; //no partners - collision
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    } //loop over track collections
00237 
00238   return result;
00239 }
00240 
00241 //
00242 //Check the number of global muons in an event, return true if there are more than 1 muon
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 //Check overlap between collections, use shared hits info
00265 //
00266 bool 
00267 MuonCosmicCompatibilityFiller::isOverlappingMuon(const edm::Event& iEvent, const edm::EventSetup& iSetup, const reco::Muon& muon) const {
00268 
00269   // 4 steps in this module
00270   // step1 : check whether it's 1leg cosmic muon or not
00271   // step2 : both muons (muons and muonsFromCosmics1Leg) should have close IP
00272   // step3 : both muons should share very close reference point
00273   // step4 : check shared hits in both muon tracks
00274 
00275   // check if this muon is available in muonsFromCosmics collection
00276   bool overlappingMuon = false; //false - not cosmic-like
00277   if( !muon.isGlobalMuon() ) return false;
00278   
00279   // reco muons for cosmics
00280   edm::Handle<reco::MuonCollection> muonHandle;
00281   iEvent.getByLabel(inputCosmicMuonCollection_, muonHandle);
00282   
00283   // Global Tracking Geometry
00284   ESHandle<GlobalTrackingGeometry> trackingGeometry;
00285   iSetup.get<GlobalTrackingGeometryRecord>().get(trackingGeometry);
00286   
00287   // PV
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       // shared hits 
00310       int RecHitsMuon = outertrack->numberOfValidHits();
00311       int RecHitsCosmicMuon = 0;
00312       int shared = 0;
00313       // count hits for same hemisphere
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         // step1
00333         if( nhitsUp > 0 && nhitsDown > 0 ) isCosmic1Leg = true;
00334         //if( !isCosmic1Leg ) continue;
00335         
00336         if( outertrack.isNonnull() ) {
00337           // step2
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           //if( !isCloseIP ) continue;
00342 
00343           // step3
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           //if( !isCloseRef ) continue;
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       // step4
00365       double fraction = -1;
00366       if( RecHitsMuon != 0 ) fraction = shared/(double)RecHitsMuon;
00367     //       std::cout << "shared = " << shared << " " << fraction << " " << RecHitsMuon << " " << RecHitsCosmicMuon << std::endl;
00368       if( shared > sharedHits_ && fraction > sharedFrac_ ) {
00369         overlappingMuon = true;
00370         break;
00371       }
00372     }
00373   }
00374   
00375   return overlappingMuon;
00376 }
00377 
00378 //
00379 //pv matches
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      //use "loose" parameter set
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          //multiple muon event
00421 
00422          if ( fabs( (*track).dxy(RefVtx) ) < maxdxyMult || fabs( (*track).dz(RefVtx) ) < maxdzMult) {
00423            result++;
00424 
00425            //case of extra large dxy
00426            if (!isLoose && fabs( (*track).dxy(RefVtx) ) > largedxyMult_) result -= 1;  
00427 
00428           }
00429      } else {
00430         //single muon event
00431 
00432         if (fabs( (*track).dxy(RefVtx) ) < maxdxy || fabs( (*track).dz(RefVtx) ) < maxdz) {
00433           result++;
00434 
00435           //case of extra large dxy
00436           if (!isLoose && fabs( (*track).dxy(RefVtx) ) > largedxy_) result -= 1; 
00437 
00438          }
00439       }
00440    }//track is nonnull
00441 }//loop over vertices
00442 
00443    //special case for non-cosmic large ip muons
00444    if (result == 0 && multipleMu) {
00445       // consider all reco muons in an event
00446       edm::Handle<reco::MuonCollection> muonHandle;
00447       iEvent.getByLabel(inputMuonCollections_[1], muonHandle);
00448 
00449       //cosmic event should have zero good vertices
00450       edm::Handle<reco::VertexCollection> pvHandle;
00451       iEvent.getByLabel(inputVertexCollection_,pvHandle);
00452       const reco::VertexCollection & vertices = *pvHandle.product();
00453 
00454       //find the "other" one
00455       if( !muonHandle.failedToGet() ) {
00456          for ( reco::MuonCollection::const_iterator muons = muonHandle->begin(); muons !=  muonHandle->end(); ++muons ) {
00457             if (!muons->isGlobalMuon()) continue;
00458             //skip this track  
00459             if ( muons->innerTrack() == muon.innerTrack() && muons->outerTrack() == muon.outerTrack()) continue;
00460                 //check ip and vertex of the "other" muon
00461                 reco::TrackRef tracks;
00462                 if (muons->isGlobalMuon()) tracks = muons->innerTrack();
00463                 if (fabs((*tracks).dxy(RefVtx)) > hIpTrdxy_) continue; 
00464                 //check if vertex collection is empty
00465                 if (vertices.begin() == vertices.end()) continue;
00466                 //for(reco::VertexCollection::const_iterator it=vertices.begin() ; it!=vertices.end() ; ++it) {
00467                     //find matching vertex by position
00468                     //if (fabs(it->z() - tracks->vz()) > 0.01) continue; //means will not be untagged from cosmics 
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   // return >=1 = identify as cosmic muon (the more like cosmics, the higher is the number) 
00486   // return 0.0 = identify as collision muon
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     //bool cosmicSegment = checkMuonSegments(muon);
00497 
00498     //short cut to reject cosmic event
00499     if (checkVertex && cosmicVertex == 0) return 10.0;
00500 
00501     // compatibility (0 - 10)
00502     // weight is assigned by the performance of individual module
00503     // btob: ~90% eff / ~0% misid
00504     // ip: ~90% eff / ~0% misid
00505     // time: ~30% eff / ~0% misid
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     // collision muon should have compatibility < 4 (0 - 4)
00512     // cosmic muon should have compatibility >= 4 (4 - 10)
00513 
00514     // b-to-b (max comp.: 4.0)
00515     if( backToback >= 1 ) {
00516        //in this case it is cosmic for sure
00517        result += weight_btob*2.;
00518        if( tightIp == 1 ) {
00519            // check with other observables to reduce mis-id (subtract compatibilities)
00520            if( looseIp == 1 ) {
00521              if( backToback < 2 ) result -= weight_btob*0.5;
00522            }
00523        }
00524     }
00525 
00526     // ip (max comp.: 4.0)
00527     if( tightIp == 0 ) {
00528         //in this case it is cosmic for sure
00529         result += weight_ip*2.0;
00530         if( backToback == 0 ) {
00531            // check with other observables to reduce mis-id (subtract compatibilities)
00532            if( tightTime == 0 ) {
00533              if( looseTime == 0 && !isOverlapping ) result -= weight_ip*1.0;
00534            }
00535         }
00536     }
00537     else if( tightIp >= 2 ) {
00538         // in this case it is almost collision-like (reduce compatibility)
00539         // if multi pvs: comp = -2.0 
00540         if( backToback >= 1 ) result -= weight_ip*1.0;
00541     }
00542 
00543     // timing (max comp.: 2.0)
00544     if( tightTime > 0 ) {
00545         // bonus track
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     // overlapping
00558     if( backToback == 0 && isOverlapping ) {
00559         // bonus track
00560         if( tightIp == 0 && tightTime >= 1 ) {
00561             result += weight_overlap*1.0;
00562         }
00563     }
00564  }//is global muon
00565 
00566 
00567 //  if (CheckMuonID && cosmicSegment) result += 4;  
00568 
00569   return result;
00570 }
00571 
00572 
00573 
00574 //
00575 //Track activity/vertex quality, count good vertices
00576 //
00577 unsigned int 
00578 MuonCosmicCompatibilityFiller::eventActivity(const edm::Event& iEvent, const reco::Muon& muon) const
00579 {
00580 
00581   unsigned int result = 0; //no good vertices - cosmic-like
00582 
00583   //check track activity
00584   edm::Handle<reco::TrackCollection> tracks;
00585   iEvent.getByLabel(inputTrackCollections_[0],tracks);
00586   if (!tracks.failedToGet() && tracks->size() < 3) return 0; 
00587 
00588   //cosmic event should have zero good vertices
00589   edm::Handle<reco::VertexCollection> pvHandle;
00590   if (!iEvent.getByLabel(inputVertexCollection_,pvHandle)) {return 0;} else {
00591       const reco::VertexCollection & vertices = *pvHandle.product();
00592       //check if vertex collection is empty
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 //Muon iD variables
00603 //
00604 bool MuonCosmicCompatibilityFiller::checkMuonID( const reco::Muon& imuon ) const
00605 {
00606   bool result = false;
00607   // initial set up using Jordan's study: GlobalMuonPromptTight + TMOneStationLoose
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 }