Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-30 23:36:34

0001 // -*- C++ -*- // // Package:  MuonSegmentMatcher // Class:  MuonSegmentMatcher // /**\class MuonSegmentMatcher MuonSegmentMatcher.cc
0002 // Description: <one line class summary>
0003 // Implementation:
0004 //     <Notes on implementation>
0005 //*/
0006 //
0007 // Original Author:  Alan Tua
0008 //         Created:  Wed Jul  9 21:40:17 CEST 2008
0009 //
0010 //
0011 
0012 #include "RecoMuon/TrackingTools/interface/MuonSegmentMatcher.h"
0013 
0014 // user include files
0015 #include "FWCore/Framework/interface/Event.h"
0016 
0017 #include "DataFormats/MuonReco/interface/Muon.h"
0018 
0019 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021 
0022 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0023 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0024 
0025 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0026 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0027 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0028 #include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h"
0029 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0030 
0031 // system include files
0032 #include <vector>
0033 #include <iostream>
0034 
0035 class MuonSegmentMatcher;
0036 
0037 using namespace std;
0038 
0039 // constructors and destructor
0040 
0041 MuonSegmentMatcher::MuonSegmentMatcher(const edm::ParameterSet& matchParameters, edm::ConsumesCollector& iC)
0042     : DTSegmentTags_(matchParameters.getParameter<edm::InputTag>("DTsegments")),
0043       CSCSegmentTags_(matchParameters.getParameter<edm::InputTag>("CSCsegments")),
0044       dtRadius_(matchParameters.getParameter<double>("DTradius")),
0045       dtTightMatch(matchParameters.getParameter<bool>("TightMatchDT")),
0046       cscTightMatch(matchParameters.getParameter<bool>("TightMatchCSC")) {
0047   if (matchParameters.existsAs<edm::InputTag>("RPChits")) {
0048     RPCHitTags_ = matchParameters.getParameter<edm::InputTag>("RPChits");
0049   } else
0050     RPCHitTags_ = edm::InputTag("hltRpcRecHits");
0051 
0052   dtRecHitsToken = iC.consumes<DTRecSegment4DCollection>(DTSegmentTags_);
0053   allSegmentsCSCToken = iC.consumes<CSCSegmentCollection>(CSCSegmentTags_);
0054   rpcRecHitsToken = iC.consumes<RPCRecHitCollection>(RPCHitTags_);
0055 }
0056 
0057 MuonSegmentMatcher::~MuonSegmentMatcher() {}
0058 
0059 // member functions
0060 vector<const DTRecSegment4D*> MuonSegmentMatcher::matchDT(const reco::Track& muon, const edm::Event& event) {
0061   using namespace edm;
0062 
0063   edm::Handle<DTRecSegment4DCollection> dtRecHits;
0064   event.getByToken(dtRecHitsToken, dtRecHits);
0065 
0066   vector<const DTRecSegment4D*> pointerTo4DSegments;
0067 
0068   std::vector<TrackingRecHit const*> dtHits;
0069 
0070   bool segments = false;
0071 
0072   // Loop and select DT recHits
0073   for (auto const& hit : muon.recHits()) {
0074     if (!hit->isValid())
0075       continue;
0076     if (hit->geographicalId().det() != DetId::Muon)
0077       continue;
0078     if (hit->geographicalId().subdetId() != MuonSubdetId::DT)
0079       continue;
0080     if (!hit->recHits().empty())
0081       if ((*hit->recHits().begin())->recHits().size() > 1)
0082         segments = true;
0083     dtHits.push_back(hit);
0084   }
0085 
0086   //  cout << "Muon DT hits found: " << dtHits.size() << " segments " << segments << endl;
0087 
0088   double PhiCutParameter = dtRadius_;
0089   double ZCutParameter = dtRadius_;
0090   double matchRatioZ = 0;
0091   double matchRatioPhi = 0;
0092 
0093   for (DTRecSegment4DCollection::const_iterator rechit = dtRecHits->begin(); rechit != dtRecHits->end(); ++rechit) {
0094     LocalPoint pointLocal = rechit->localPosition();
0095 
0096     if (segments) {
0097       // Loop over muon recHits
0098       for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0099         // Pick the one in the same DT Chamber as the muon
0100         DetId idT = (*hit)->geographicalId();
0101         if (!(rechit->geographicalId().rawId() == idT.rawId()))
0102           continue;
0103 
0104         // and compare the local positions
0105         LocalPoint segLocal = (*hit)->localPosition();
0106         if ((fabs(pointLocal.x() - segLocal.x()) < ZCutParameter) &&
0107             (fabs(pointLocal.y() - segLocal.y()) < ZCutParameter))
0108           pointerTo4DSegments.push_back(&(*rechit));
0109       }
0110       continue;
0111     }
0112 
0113     double nhitsPhi = 0;
0114     double nhitsZ = 0;
0115 
0116     if (rechit->hasZed()) {
0117       double countMuonDTHits = 0;
0118       double countAgreeingHits = 0;
0119 
0120       const DTRecSegment2D* segmZ;
0121       segmZ = dynamic_cast<const DTRecSegment2D*>(rechit->zSegment());
0122       nhitsZ = segmZ->recHits().size();
0123 
0124       const vector<DTRecHit1D> hits1d = segmZ->specificRecHits();
0125       DTChamberId chamberSegIdT((segmZ->geographicalId()).rawId());
0126 
0127       // Loop over muon recHits
0128       for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0129         if (!(*hit)->isValid())
0130           continue;
0131 
0132         DetId idT = (*hit)->geographicalId();
0133         DTChamberId dtDetIdHitT(idT.rawId());
0134         DTSuperLayerId dtDetLayerIdHitT(idT.rawId());
0135 
0136         LocalPoint pointLocal = (*hit)->localPosition();
0137 
0138         if ((chamberSegIdT == dtDetIdHitT) && (dtDetLayerIdHitT.superlayer() == 2))
0139           countMuonDTHits++;
0140 
0141         for (vector<DTRecHit1D>::const_iterator hiti = hits1d.begin(); hiti != hits1d.end(); hiti++) {
0142           if (!hiti->isValid())
0143             continue;
0144 
0145           // Pick the one in the same DT Layer as the 1D hit
0146           if (!(hiti->geographicalId().rawId() == idT.rawId()))
0147             continue;
0148 
0149           // and compare the local positions
0150           LocalPoint segLocal = hiti->localPosition();
0151           //      cout << "Zed Segment Point = "<<pointLocal<<"    Muon Point = "<<segLocal<<"  Dist:  "
0152           //           << (fabs(pointLocal.x()-segLocal.x()))+(fabs(pointLocal.y()-segLocal.y()))<< endl;
0153           if ((fabs(pointLocal.x() - segLocal.x()) < ZCutParameter) &&
0154               (fabs(pointLocal.y() - segLocal.y()) < ZCutParameter))
0155             countAgreeingHits++;
0156         }  //End Segment Hit Iteration
0157       }    //End Muon Hit Iteration
0158 
0159       matchRatioZ = countMuonDTHits == 0 ? 0 : countAgreeingHits / countMuonDTHits;
0160       if (nhitsZ)
0161         if (countAgreeingHits / nhitsZ > matchRatioZ)
0162           matchRatioZ = countAgreeingHits / nhitsZ;
0163     }  //End HasZed Check
0164 
0165     if (rechit->hasPhi()) {
0166       double countMuonDTHits = 0;
0167       double countAgreeingHits = 0;
0168 
0169       //PREPARE PARAMETERS FOR SEGMENT DETECTOR GEOMETRY
0170       const DTRecSegment2D* segmPhi;
0171       segmPhi = dynamic_cast<const DTRecSegment2D*>(rechit->phiSegment());
0172       nhitsPhi = segmPhi->recHits().size();
0173 
0174       const vector<DTRecHit1D> hits1d = segmPhi->specificRecHits();
0175       DTChamberId chamberSegIdT((segmPhi->geographicalId()).rawId());
0176 
0177       // Loop over muon recHits
0178       for (auto hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
0179         if (!(*hit)->isValid())
0180           continue;
0181 
0182         DetId idT = (*hit)->geographicalId();
0183         DTChamberId dtDetIdHitT(idT.rawId());
0184         DTSuperLayerId dtDetLayerIdHitT(idT.rawId());
0185 
0186         LocalPoint pointLocal =
0187             (*hit)
0188                 ->localPosition();  //Localposition is in DTLayer http://cmslxr.fnal.gov/lxr/source/DataFormats/DTRecHit/interface/DTRecHit1D.h
0189 
0190         if ((chamberSegIdT == dtDetIdHitT) &&
0191             ((dtDetLayerIdHitT.superlayer() == 1) || (dtDetLayerIdHitT.superlayer() == 3)))
0192           countMuonDTHits++;
0193 
0194         for (vector<DTRecHit1D>::const_iterator hiti = hits1d.begin(); hiti != hits1d.end(); hiti++) {
0195           if (!hiti->isValid())
0196             continue;
0197 
0198           // Pick the one in the same DT Layer as the 1D hit
0199           if (!(hiti->geographicalId().rawId() == idT.rawId()))
0200             continue;
0201 
0202           // and compare the local positions
0203           LocalPoint segLocal = hiti->localPosition();
0204           //      cout << "     Phi Segment Point = "<<pointLocal<<"    Muon Point = "<<segLocal<<"  Dist:   "
0205           //           << (fabs(pointLocal.x()-segLocal.x()))+(fabs(pointLocal.y()-segLocal.y()))<< endl;
0206 
0207           if ((fabs(pointLocal.x() - segLocal.x()) < PhiCutParameter) &&
0208               (fabs(pointLocal.y() - segLocal.y()) < PhiCutParameter))
0209             countAgreeingHits++;
0210         }  // End Segment Hit Iteration
0211       }    // End Muon Hit Iteration
0212 
0213       matchRatioPhi = countMuonDTHits != 0 ? countAgreeingHits / countMuonDTHits : 0;
0214       if (nhitsPhi)
0215         if (countAgreeingHits / nhitsPhi > matchRatioPhi)
0216           matchRatioPhi = countAgreeingHits / nhitsPhi;
0217     }  // End HasPhi Check
0218        //    DTChamberId chamberSegId2((rechit->geographicalId()).rawId());
0219     if (dtTightMatch && nhitsPhi && nhitsZ) {
0220       if ((matchRatioPhi > 0.9) && (matchRatioZ > 0.9)) {
0221         //  cout<<"Making a tight match in Chamber "<<chamberSegId2<<endl;
0222         pointerTo4DSegments.push_back(&(*rechit));
0223       }
0224     } else {
0225       if ((matchRatioPhi > 0.9 && nhitsPhi) || (matchRatioZ > 0.9 && nhitsZ)) {
0226         //  cout<<"Making a loose match in Chamber "<<chamberSegId2<<endl;
0227         pointerTo4DSegments.push_back(&(*rechit));
0228       }
0229     }
0230 
0231   }  //End DT Segment Iteration
0232 
0233   return pointerTo4DSegments;
0234 }
0235 
0236 vector<const CSCSegment*> MuonSegmentMatcher::matchCSC(const reco::Track& muon, const edm::Event& event) {
0237   using namespace edm;
0238 
0239   edm::Handle<CSCSegmentCollection> allSegmentsCSC;
0240   event.getByToken(allSegmentsCSCToken, allSegmentsCSC);
0241 
0242   vector<const CSCSegment*> pointerToCSCSegments;
0243 
0244   double matchRatioCSC = 0;
0245   int numCSC = 0;
0246   double CSCXCut = 0.001;
0247   double CSCYCut = 0.001;
0248   double countMuonCSCHits = 0;
0249 
0250   for (CSCSegmentCollection::const_iterator segmentCSC = allSegmentsCSC->begin(); segmentCSC != allSegmentsCSC->end();
0251        segmentCSC++) {
0252     double CSCcountAgreeingHits = 0;
0253 
0254     if (!segmentCSC->isValid())
0255       continue;
0256 
0257     numCSC++;
0258     const vector<CSCRecHit2D>& CSCRechits2D = segmentCSC->specificRecHits();
0259     countMuonCSCHits = 0;
0260     CSCDetId myChamber((*segmentCSC).geographicalId().rawId());
0261 
0262     bool segments = false;
0263 
0264     for (auto const& hitC : muon.recHits()) {
0265       if (!hitC->isValid())
0266         continue;
0267       if (hitC->geographicalId().det() != DetId::Muon)
0268         continue;
0269       if (hitC->geographicalId().subdetId() != MuonSubdetId::CSC)
0270         continue;
0271       if (!hitC->isValid())
0272         continue;
0273       if (hitC->recHits().size() > 1)
0274         segments = true;
0275 
0276       //DETECTOR CONSTRUCTION
0277       DetId id = hitC->geographicalId();
0278       CSCDetId cscDetIdHit(id.rawId());
0279 
0280       if (segments) {
0281         if (!(myChamber.rawId() == cscDetIdHit.rawId()))
0282           continue;
0283 
0284         // and compare the local positions
0285         LocalPoint positionLocalCSC = hitC->localPosition();
0286         LocalPoint segLocalCSC = segmentCSC->localPosition();
0287         if ((fabs(positionLocalCSC.x() - segLocalCSC.x()) < CSCXCut) &&
0288             (fabs(positionLocalCSC.y() - segLocalCSC.y()) < CSCYCut))
0289           pointerToCSCSegments.push_back(&(*segmentCSC));
0290         continue;
0291       }
0292 
0293       if (!(cscDetIdHit.ring() == myChamber.ring()))
0294         continue;
0295       if (!(cscDetIdHit.station() == myChamber.station()))
0296         continue;
0297       if (!(cscDetIdHit.endcap() == myChamber.endcap()))
0298         continue;
0299       if (!(cscDetIdHit.chamber() == myChamber.chamber()))
0300         continue;
0301 
0302       countMuonCSCHits++;
0303 
0304       LocalPoint positionLocalCSC = hitC->localPosition();
0305 
0306       for (vector<CSCRecHit2D>::const_iterator hiti = CSCRechits2D.begin(); hiti != CSCRechits2D.end(); hiti++) {
0307         if (!hiti->isValid())
0308           continue;
0309         CSCDetId cscDetId((hiti->geographicalId()).rawId());
0310 
0311         if (hitC->geographicalId().rawId() != (hiti->geographicalId()).rawId())
0312           continue;
0313 
0314         LocalPoint segLocalCSC = hiti->localPosition();
0315         //      cout<<"Layer Id (MuonHit) =  "<<cscDetIdHit<<" Muon Local Position (det frame) "<<positionLocalCSC <<endl;
0316         //      cout<<"Layer Id  (CSCHit) =  "<<cscDetId<<"  Hit Local Position (det frame) "<<segLocalCSC <<endl;
0317         if ((fabs(positionLocalCSC.x() - segLocalCSC.x()) < CSCXCut) &&
0318             (fabs(positionLocalCSC.y() - segLocalCSC.y()) < CSCYCut)) {
0319           CSCcountAgreeingHits++;
0320           //          cout << "   Matched." << endl;
0321         }
0322       }  //End 2D rechit iteration
0323     }    //End muon hit iteration
0324 
0325     matchRatioCSC = countMuonCSCHits == 0 ? 0 : CSCcountAgreeingHits / countMuonCSCHits;
0326 
0327     if ((matchRatioCSC > 0.9) && ((countMuonCSCHits > 1) || !cscTightMatch))
0328       pointerToCSCSegments.push_back(&(*segmentCSC));
0329 
0330   }  //End CSC Segment Iteration
0331 
0332   return pointerToCSCSegments;
0333 }
0334 
0335 vector<const RPCRecHit*> MuonSegmentMatcher::matchRPC(const reco::Track& muon, const edm::Event& event) {
0336   using namespace edm;
0337 
0338   edm::Handle<RPCRecHitCollection> rpcRecHits;
0339   event.getByToken(rpcRecHitsToken, rpcRecHits);
0340 
0341   vector<const RPCRecHit*> pointerToRPCRecHits;
0342   double RPCCut = 0.001;
0343 
0344   for (RPCRecHitCollection::const_iterator hitRPC = rpcRecHits->begin(); hitRPC != rpcRecHits->end(); hitRPC++) {
0345     if (!hitRPC->isValid())
0346       continue;
0347 
0348     RPCDetId myChamber((*hitRPC).geographicalId().rawId());
0349     LocalPoint posLocalRPC = hitRPC->localPosition();
0350     bool matched = false;
0351 
0352     for (auto const& hitC : muon.recHits()) {
0353       if (!hitC->isValid())
0354         continue;
0355       if (hitC->geographicalId().det() != DetId::Muon)
0356         continue;
0357       if (hitC->geographicalId().subdetId() != MuonSubdetId::RPC)
0358         continue;
0359       if (!hitC->isValid())
0360         continue;
0361 
0362       //DETECTOR CONSTRUCTION
0363       DetId id = hitC->geographicalId();
0364       RPCDetId rpcDetIdHit(id.rawId());
0365 
0366       if (rpcDetIdHit != myChamber)
0367         continue;
0368       LocalPoint posLocalMuon = hitC->localPosition();
0369 
0370       //        cout<<"Layer Id (MuonHit) =  "<<rpcDetIdHit<<" Muon Local Position (det frame) "<<posLocalMuon <<endl;
0371       //        cout<<"Layer Id  (RPCHit) =  "<<myChamber<<"  Hit Local Position (det frame) "<<posLocalRPC <<endl;
0372       if ((fabs(posLocalMuon.x() - posLocalRPC.x()) < RPCCut)) {
0373         matched = true;
0374         break;
0375       }
0376     }
0377 
0378     if (matched)
0379       pointerToRPCRecHits.push_back(&(*hitRPC));
0380   }
0381 
0382   return pointerToRPCRecHits;
0383 }
0384 
0385 //define this as a plug-in
0386 //DEFINE_FWK_MODULE(MuonSegmentMatcher);