Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*- // // Package:  MuonSegmentMatcher // Class:  MuonSegmentMatcher // /**\class MuonSegmentMatcher MuonSegmentMatcher.cc
00002 // Description: <one line class summary>
00003 // Implementation:
00004 //     <Notes on implementation>
00005 //*/
00006 //
00007 // Original Author:  Alan Tua
00008 //         Created:  Wed Jul  9 21:40:17 CEST 2008
00009 // $Id: MuonSegmentMatcher.cc,v 1.12 2011/11/09 12:19:25 ptraczyk Exp $
00010 //
00011 //
00012 
00013 #include "RecoMuon/TrackingTools/interface/MuonSegmentMatcher.h"
00014 
00015 
00016 // user include files
00017 #include "FWCore/Framework/interface/EDAnalyzer.h"
00018 #include "FWCore/Framework/interface/Event.h"
00019 
00020 #include "DataFormats/MuonReco/interface/Muon.h"
00021 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
00022 
00023 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00024 #include "DataFormats/TrackReco/interface/Track.h"
00025 
00026 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00027 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00028 
00029 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
00030 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
00031 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
00032 #include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h"
00033 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00034 
00035 // system include files
00036 #include <vector>
00037 #include <iostream>
00038 
00039 class MuonServiceProxy;
00040 class MuonSegmentMatcher;
00041 
00042 using namespace std;
00043 
00044 // constructors and destructor
00045 
00046 MuonSegmentMatcher::MuonSegmentMatcher(const edm::ParameterSet& matchParameters, MuonServiceProxy* service)
00047   :
00048   theService(service),
00049   DTSegmentTags_(matchParameters.getParameter<edm::InputTag>("DTsegments")),
00050   CSCSegmentTags_(matchParameters.getParameter<edm::InputTag>("CSCsegments")),
00051   dtRadius_(matchParameters.getParameter<double>("DTradius")),
00052   dtTightMatch(matchParameters.getParameter<bool>("TightMatchDT")),
00053   cscTightMatch(matchParameters.getParameter<bool>("TightMatchCSC"))
00054 {
00055 }
00056 
00057 MuonSegmentMatcher::~MuonSegmentMatcher()
00058 {
00059 }
00060 
00061 // member functions
00062 vector<const DTRecSegment4D*> MuonSegmentMatcher::matchDT(const reco::Track &muon, const edm::Event& event)
00063 {
00064   using namespace edm;
00065 
00066   edm::Handle<DTRecSegment4DCollection> dtRecHits;
00067   event.getByLabel(DTSegmentTags_, dtRecHits);  
00068   
00069   vector<const DTRecSegment4D*> pointerTo4DSegments;
00070 
00071   TrackingRecHitRefVector dtHits;
00072 
00073   bool segments = false;
00074 
00075   // Loop and select DT recHits
00076   for(trackingRecHit_iterator hit = muon.recHitsBegin(); hit != muon.recHitsEnd(); ++hit) {
00077     if ( !(*hit)->isValid()) continue; 
00078     if ( (*hit)->geographicalId().det() != DetId::Muon ) continue; 
00079     if ( (*hit)->geographicalId().subdetId() != MuonSubdetId::DT ) continue; 
00080     if (!(*hit)->isValid()) continue; 
00081     if ((*hit)->recHits().size()>1) segments = true;
00082     dtHits.push_back(*hit);
00083   }
00084   
00085   double PhiCutParameter=dtRadius_;
00086   double ZCutParameter=dtRadius_;
00087   double matchRatioZ=0;
00088   double matchRatioPhi=0;
00089 
00090   for (DTRecSegment4DCollection::const_iterator rechit = dtRecHits->begin(); rechit!=dtRecHits->end();++rechit) {
00091   
00092     if ( !rechit->isValid()) continue; 
00093     LocalPoint pointLocal = rechit->localPosition();
00094 
00095     if (segments) {
00096       // Loop over muon recHits
00097       for(trackingRecHit_iterator hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
00098         if ( !(*hit)->isValid()) continue; 
00099                                         
00100         // Pick the one in the same DT Chamber as the muon
00101         DetId idT = (*hit)->geographicalId();
00102         if(!(rechit->geographicalId().rawId()==idT.rawId())) continue; 
00103 
00104         // and compare the local positions
00105         LocalPoint segLocal = (*hit)->localPosition();
00106         if ((fabs(pointLocal.x()-segLocal.x())<ZCutParameter) && 
00107             (fabs(pointLocal.y()-segLocal.y())<ZCutParameter)) 
00108           pointerTo4DSegments.push_back(&(*rechit));
00109       }
00110       continue;
00111     }
00112 
00113     double nhitsPhi = 0;
00114     double nhitsZ = 0;
00115           
00116     if(rechit->hasZed()) {
00117       double countMuonDTHits = 0;
00118       double countAgreeingHits=0;
00119 
00120       const DTRecSegment2D* segmZ;
00121       segmZ = dynamic_cast<const DTRecSegment2D*>(rechit->zSegment());
00122       nhitsZ = segmZ->recHits().size(); 
00123                 
00124       const vector<DTRecHit1D> hits1d = segmZ->specificRecHits();
00125       DTChamberId chamberSegIdT((segmZ->geographicalId()).rawId());
00126                 
00127       // Loop over muon recHits
00128       for(trackingRecHit_iterator hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
00129 
00130         if ( !(*hit)->isValid()) continue; 
00131         
00132         DetId idT = (*hit)->geographicalId();
00133         DTChamberId dtDetIdHitT(idT.rawId());
00134         DTSuperLayerId dtDetLayerIdHitT(idT.rawId());
00135 
00136         LocalPoint  pointLocal = (*hit)->localPosition();
00137                                         
00138         if ((chamberSegIdT==dtDetIdHitT) && (dtDetLayerIdHitT.superlayer()==2)) countMuonDTHits++;
00139 
00140         for (vector<DTRecHit1D>::const_iterator hiti=hits1d.begin(); hiti!=hits1d.end(); hiti++) {
00141 
00142           if ( !hiti->isValid()) continue; 
00143 
00144           // Pick the one in the same DT Layer as the 1D hit
00145           if(!(hiti->geographicalId().rawId()==idT.rawId())) continue; 
00146 
00147           // and compare the local positions
00148           LocalPoint segLocal = hiti->localPosition();
00149 //        cout << "Zed Segment Point = "<<pointLocal<<"    Muon Point = "<<segLocal<<"  Dist:  "
00150 //             << (fabs(pointLocal.x()-segLocal.x()))+(fabs(pointLocal.y()-segLocal.y()))<< endl;
00151           if ((fabs(pointLocal.x()-segLocal.x())<ZCutParameter) && 
00152               (fabs(pointLocal.y()-segLocal.y())<ZCutParameter)) 
00153             countAgreeingHits++;
00154         } //End Segment Hit Iteration
00155       } //End Muon Hit Iteration
00156                 
00157       matchRatioZ = countMuonDTHits == 0 ? 0 : countAgreeingHits/countMuonDTHits;
00158       if (nhitsZ)
00159         if (countAgreeingHits/nhitsZ>matchRatioZ) matchRatioZ=countAgreeingHits/nhitsZ;
00160     } //End HasZed Check
00161                         
00162     if(rechit->hasPhi()) {
00163       double countMuonDTHits = 0;
00164       double countAgreeingHits=0;
00165 
00166       //PREPARE PARAMETERS FOR SEGMENT DETECTOR GEOMETRY
00167       const DTRecSegment2D* segmPhi;
00168       segmPhi = dynamic_cast<const DTRecSegment2D*>(rechit->phiSegment());
00169       nhitsPhi = segmPhi->recHits().size();
00170                 
00171       const vector<DTRecHit1D> hits1d = segmPhi->specificRecHits();
00172       DTChamberId chamberSegIdT((segmPhi->geographicalId()).rawId());
00173                 
00174       // Loop over muon recHits
00175       for(trackingRecHit_iterator hit = dtHits.begin(); hit != dtHits.end(); ++hit) {
00176 
00177         if ( !(*hit)->isValid()) continue; 
00178 
00179         DetId idT = (*hit)->geographicalId();
00180         DTChamberId dtDetIdHitT(idT.rawId());
00181         DTSuperLayerId dtDetLayerIdHitT(idT.rawId());
00182 
00183         LocalPoint pointLocal = (*hit)->localPosition(); //Localposition is in DTLayer http://cmslxr.fnal.gov/lxr/source/DataFormats/DTRecHit/interface/DTRecHit1D.h
00184 
00185         if ((chamberSegIdT==dtDetIdHitT)&&((dtDetLayerIdHitT.superlayer()==1)||(dtDetLayerIdHitT.superlayer()==3))) 
00186           countMuonDTHits++;
00187 
00188         for (vector<DTRecHit1D>::const_iterator hiti=hits1d.begin(); hiti!=hits1d.end(); hiti++) {
00189 
00190           if ( !hiti->isValid()) continue; 
00191 
00192           // Pick the one in the same DT Layer as the 1D hit
00193           if(!(hiti->geographicalId().rawId()==idT.rawId())) continue; 
00194                                          
00195           // and compare the local positions
00196           LocalPoint segLocal = hiti->localPosition();
00197 //        cout << "     Phi Segment Point = "<<pointLocal<<"    Muon Point = "<<segLocal<<"  Dist:   " 
00198 //             << (fabs(pointLocal.x()-segLocal.x()))+(fabs(pointLocal.y()-segLocal.y()))<< endl;
00199 
00200           if ((fabs(pointLocal.x()-segLocal.x())<PhiCutParameter) && 
00201               (fabs(pointLocal.y()-segLocal.y())<PhiCutParameter))
00202             countAgreeingHits++; 
00203         } // End Segment Hit Iteration
00204       } // End Muon Hit Iteration
00205 
00206       matchRatioPhi = countMuonDTHits != 0 ? countAgreeingHits/countMuonDTHits : 0;
00207       if (nhitsPhi)
00208         if (countAgreeingHits/nhitsPhi>matchRatioPhi) matchRatioPhi=countAgreeingHits/nhitsPhi;
00209     } // End HasPhi Check
00210 //    DTChamberId chamberSegId2((rechit->geographicalId()).rawId());
00211     if (dtTightMatch && nhitsPhi && nhitsZ) {
00212       if((matchRatioPhi>0.9)&&(matchRatioZ>0.9)) {
00213 //      cout<<"Making a tight match in Chamber "<<chamberSegId2<<endl;
00214         pointerTo4DSegments.push_back(&(*rechit));
00215       }
00216     } else {
00217       if((matchRatioPhi>0.9 && nhitsPhi)||(matchRatioZ>0.9 && nhitsZ)) {
00218 //      cout<<"Making a loose match in Chamber "<<chamberSegId2<<endl;
00219         pointerTo4DSegments.push_back(&(*rechit));
00220       }
00221     }
00222     
00223   } //End DT Segment Iteration
00224 
00225   return pointerTo4DSegments;
00226 }
00227 
00228 
00229 
00230 vector<const CSCSegment*> MuonSegmentMatcher::matchCSC(const reco::Track& muon, const edm::Event& event)
00231 {
00232 
00233   using namespace edm;
00234 
00235   edm::Handle<CSCSegmentCollection> allSegmentsCSC;
00236   event.getByLabel(CSCSegmentTags_, allSegmentsCSC);
00237 
00238   vector<const CSCSegment*> pointerToCSCSegments;
00239 
00240   double matchRatioCSC=0;
00241   int numCSC = 0;
00242   double CSCXCut = 0.001;
00243   double CSCYCut = 0.001;
00244   double countMuonCSCHits = 0;
00245 
00246   for(CSCSegmentCollection::const_iterator segmentCSC = allSegmentsCSC->begin(); segmentCSC != allSegmentsCSC->end(); segmentCSC++) {
00247     double CSCcountAgreeingHits=0;
00248 
00249     if ( !segmentCSC->isValid()) continue; 
00250 
00251     numCSC++;
00252     const vector<CSCRecHit2D>& CSCRechits2D = segmentCSC->specificRecHits();
00253     countMuonCSCHits = 0;
00254     CSCDetId myChamber((*segmentCSC).geographicalId().rawId());
00255 
00256     bool segments = false;
00257 
00258     for(trackingRecHit_iterator hitC = muon.recHitsBegin(); hitC != muon.recHitsEnd(); ++hitC) {
00259       if (!(*hitC)->isValid()) continue; 
00260       if ( (*hitC)->geographicalId().det() != DetId::Muon ) continue; 
00261       if ( (*hitC)->geographicalId().subdetId() != MuonSubdetId::CSC ) continue;
00262       if (!(*hitC)->isValid()) continue;
00263       if ( (*hitC)->recHits().size()>1) segments = true;
00264 
00265       //DETECTOR CONSTRUCTION
00266       DetId id = (*hitC)->geographicalId();
00267       CSCDetId cscDetIdHit(id.rawId());
00268 
00269       if (segments) {
00270         if(!(myChamber.rawId()==cscDetIdHit.rawId())) continue; 
00271 
00272         // and compare the local positions
00273         LocalPoint positionLocalCSC = (*hitC)->localPosition();
00274         LocalPoint segLocalCSC = segmentCSC->localPosition();
00275         if ((fabs(positionLocalCSC.x()-segLocalCSC.x())<CSCXCut) && 
00276             (fabs(positionLocalCSC.y()-segLocalCSC.y())<CSCYCut)) 
00277           pointerToCSCSegments.push_back(&(*segmentCSC)); 
00278         continue;
00279       }
00280 
00281       if(!(cscDetIdHit.ring()==myChamber.ring())) continue;
00282       if(!(cscDetIdHit.station()==myChamber.station())) continue;
00283       if(!(cscDetIdHit.endcap()==myChamber.endcap())) continue;
00284       if(!(cscDetIdHit.chamber()==myChamber.chamber())) continue;
00285 
00286       countMuonCSCHits++;
00287 
00288       LocalPoint positionLocalCSC = (*hitC)->localPosition();
00289         
00290       for (vector<CSCRecHit2D>::const_iterator hiti=CSCRechits2D.begin(); hiti!=CSCRechits2D.end(); hiti++) {
00291 
00292         if ( !hiti->isValid()) continue; 
00293         CSCDetId cscDetId((hiti->geographicalId()).rawId());
00294                 
00295         if ((*hitC)->geographicalId().rawId()!=(hiti->geographicalId()).rawId()) continue;
00296 
00297         LocalPoint segLocalCSC = hiti->localPosition();
00298         //              cout<<"Layer Id (MuonHit) =  "<<cscDetIdHit<<" Muon Local Position (det frame) "<<positionLocalCSC <<endl;
00299         //              cout<<"Layer Id  (CSCHit) =  "<<cscDetId<<"  Hit Local Position (det frame) "<<segLocalCSC <<endl;
00300         if((fabs(positionLocalCSC.x()-segLocalCSC.x())<CSCXCut) && 
00301            (fabs(positionLocalCSC.y()-segLocalCSC.y())<CSCYCut)) {
00302           CSCcountAgreeingHits++;
00303           //              cout << "   Matched." << endl;
00304         }  
00305       }//End 2D rechit iteration
00306     }//End muon hit iteration
00307     
00308     matchRatioCSC = countMuonCSCHits == 0 ? 0 : CSCcountAgreeingHits/countMuonCSCHits;
00309                 
00310     if ((matchRatioCSC>0.9) && ((countMuonCSCHits>1) || !cscTightMatch)) pointerToCSCSegments.push_back(&(*segmentCSC));
00311 
00312   } //End CSC Segment Iteration 
00313 
00314   return pointerToCSCSegments;
00315 
00316 }
00317 
00318 //define this as a plug-in
00319 //DEFINE_FWK_MODULE(MuonSegmentMatcher);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoMuon/TrackingTools/src/MuonSegmentMatcher.cc
