00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include <memory>
00015
00016
00017 #include "FWCore/Framework/interface/Frameworkfwd.h"
00018 #include "FWCore/Framework/interface/EDProducer.h"
00019
00020 #include "FWCore/Framework/interface/Event.h"
00021 #include "FWCore/Framework/interface/EventSetup.h"
00022
00023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00024
00025 #include "DataFormats/Common/interface/Handle.h"
00026 #include "DataFormats/TrackReco/interface/Track.h"
00027 #include "DataFormats/MuonReco/interface/Muon.h"
00028 #include "DataFormats/MuonReco/interface/CaloMuon.h"
00029 #include "DataFormats/MuonReco/interface/MuonCocktails.h"
00030 #include "DataFormats/MuonReco/interface/MuonTime.h"
00031 #include "DataFormats/MuonReco/interface/MuonTimeExtra.h"
00032 #include "DataFormats/MuonReco/interface/MuonTimeExtraMap.h"
00033 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
00034 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
00035
00036 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
00037 #include "Utilities/Timing/interface/TimerStack.h"
00038
00039 #include <boost/regex.hpp>
00040 #include "RecoMuon/MuonIdentification/plugins/MuonIdProducer.h"
00041 #include "RecoMuon/MuonIdentification/interface/MuonIdTruthInfo.h"
00042 #include "RecoMuon/MuonIdentification/interface/MuonArbitrationMethods.h"
00043 #include "RecoMuon/MuonIdentification/interface/MuonMesh.h"
00044
00045 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractorFactory.h"
00046 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00047
00048 #include <algorithm>
00049
00050 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00051 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
00052 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00053 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
00054
00055
00056 #include "DataFormats/RPCRecHit/interface/RPCRecHit.h"
00057 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
00058 #include "DataFormats/MuonReco/interface/MuonRPCHitMatch.h"
00059
00060 #include "DataFormats/Common/interface/ValueMap.h"
00061
00062 #include "RecoMuon/MuonIdentification/interface/MuonKinkFinder.h"
00063
00064 MuonIdProducer::MuonIdProducer(const edm::ParameterSet& iConfig):
00065 muIsoExtractorCalo_(0),muIsoExtractorTrack_(0),muIsoExtractorJet_(0)
00066 {
00067 produces<reco::MuonCollection>();
00068 produces<reco::CaloMuonCollection>();
00069 produces<reco::MuonTimeExtraMap>("combined");
00070 produces<reco::MuonTimeExtraMap>("dt");
00071 produces<reco::MuonTimeExtraMap>("csc");
00072
00073 minPt_ = iConfig.getParameter<double>("minPt");
00074 minP_ = iConfig.getParameter<double>("minP");
00075 minPCaloMuon_ = iConfig.getParameter<double>("minPCaloMuon");
00076 minNumberOfMatches_ = iConfig.getParameter<int>("minNumberOfMatches");
00077 addExtraSoftMuons_ = iConfig.getParameter<bool>("addExtraSoftMuons");
00078 maxAbsEta_ = iConfig.getParameter<double>("maxAbsEta");
00079 maxAbsDx_ = iConfig.getParameter<double>("maxAbsDx");
00080 maxAbsPullX_ = iConfig.getParameter<double>("maxAbsPullX");
00081 maxAbsDy_ = iConfig.getParameter<double>("maxAbsDy");
00082 maxAbsPullY_ = iConfig.getParameter<double>("maxAbsPullY");
00083 fillCaloCompatibility_ = iConfig.getParameter<bool>("fillCaloCompatibility");
00084 fillEnergy_ = iConfig.getParameter<bool>("fillEnergy");
00085 fillMatching_ = iConfig.getParameter<bool>("fillMatching");
00086 fillIsolation_ = iConfig.getParameter<bool>("fillIsolation");
00087 writeIsoDeposits_ = iConfig.getParameter<bool>("writeIsoDeposits");
00088 fillGlobalTrackQuality_ = iConfig.getParameter<bool>("fillGlobalTrackQuality");
00089 fillGlobalTrackRefits_ = iConfig.getParameter<bool>("fillGlobalTrackRefits");
00090
00091 fillTrackerKink_ = false;
00092 if (fillGlobalTrackQuality_) fillTrackerKink_ = iConfig.getParameter<bool>("fillTrackerKink");
00093
00094 ptThresholdToFillCandidateP4WithGlobalFit_ = iConfig.getParameter<double>("ptThresholdToFillCandidateP4WithGlobalFit");
00095 sigmaThresholdToFillCandidateP4WithGlobalFit_ = iConfig.getParameter<double>("sigmaThresholdToFillCandidateP4WithGlobalFit");
00096 caloCut_ = iConfig.getParameter<double>("minCaloCompatibility");
00097 arbClean_ = iConfig.getParameter<bool>("runArbitrationCleaner");
00098
00099
00100 edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
00101 parameters_.loadParameters( parameters );
00102
00103
00104 edm::ParameterSet timingParameters = iConfig.getParameter<edm::ParameterSet>("TimingFillerParameters");
00105 theTimingFiller_ = new MuonTimingFiller(timingParameters);
00106
00107 if (fillCaloCompatibility_){
00108
00109 parameters = iConfig.getParameter<edm::ParameterSet>("MuonCaloCompatibility");
00110 muonCaloCompatibility_.configure( parameters );
00111 }
00112
00113 if (fillIsolation_){
00114
00115 edm::ParameterSet caloExtractorPSet = iConfig.getParameter<edm::ParameterSet>("CaloExtractorPSet");
00116 std::string caloExtractorName = caloExtractorPSet.getParameter<std::string>("ComponentName");
00117 muIsoExtractorCalo_ = IsoDepositExtractorFactory::get()->create( caloExtractorName, caloExtractorPSet);
00118
00119 edm::ParameterSet trackExtractorPSet = iConfig.getParameter<edm::ParameterSet>("TrackExtractorPSet");
00120 std::string trackExtractorName = trackExtractorPSet.getParameter<std::string>("ComponentName");
00121 muIsoExtractorTrack_ = IsoDepositExtractorFactory::get()->create( trackExtractorName, trackExtractorPSet);
00122
00123 edm::ParameterSet jetExtractorPSet = iConfig.getParameter<edm::ParameterSet>("JetExtractorPSet");
00124 std::string jetExtractorName = jetExtractorPSet.getParameter<std::string>("ComponentName");
00125 muIsoExtractorJet_ = IsoDepositExtractorFactory::get()->create( jetExtractorName, jetExtractorPSet);
00126 }
00127 if (fillIsolation_ && writeIsoDeposits_){
00128 trackDepositName_ = iConfig.getParameter<std::string>("trackDepositName");
00129 produces<reco::IsoDepositMap>(trackDepositName_);
00130 ecalDepositName_ = iConfig.getParameter<std::string>("ecalDepositName");
00131 produces<reco::IsoDepositMap>(ecalDepositName_);
00132 hcalDepositName_ = iConfig.getParameter<std::string>("hcalDepositName");
00133 produces<reco::IsoDepositMap>(hcalDepositName_);
00134 hoDepositName_ = iConfig.getParameter<std::string>("hoDepositName");
00135 produces<reco::IsoDepositMap>(hoDepositName_);
00136 jetDepositName_ = iConfig.getParameter<std::string>("jetDepositName");
00137 produces<reco::IsoDepositMap>(jetDepositName_);
00138 }
00139
00140 inputCollectionLabels_ = iConfig.getParameter<std::vector<edm::InputTag> >("inputCollectionLabels");
00141 inputCollectionTypes_ = iConfig.getParameter<std::vector<std::string> >("inputCollectionTypes");
00142 if (inputCollectionLabels_.size() != inputCollectionTypes_.size())
00143 throw cms::Exception("ConfigurationError") << "Number of input collection labels is different from number of types. " <<
00144 "For each collection label there should be exactly one collection type specified.";
00145 if (inputCollectionLabels_.size()>7 ||inputCollectionLabels_.empty())
00146 throw cms::Exception("ConfigurationError") << "Number of input collections should be from 1 to 7.";
00147
00148 debugWithTruthMatching_ = iConfig.getParameter<bool>("debugWithTruthMatching");
00149 if (debugWithTruthMatching_) edm::LogWarning("MuonIdentification")
00150 << "========================================================================\n"
00151 << "Debugging mode with truth matching is turned on!!! Make sure you understand what you are doing!\n"
00152 << "========================================================================\n";
00153 if (fillGlobalTrackQuality_){
00154 globalTrackQualityInputTag_ = iConfig.getParameter<edm::InputTag>("globalTrackQualityInputTag");
00155 }
00156
00157 if (fillTrackerKink_) {
00158 trackerKinkFinder_.reset(new MuonKinkFinder(iConfig.getParameter<edm::ParameterSet>("TrackerKinkFinderParameters")));
00159 }
00160
00161
00162 meshAlgo_ = new MuonMesh(iConfig.getParameter<edm::ParameterSet>("arbitrationCleanerOptions"));
00163 }
00164
00165
00166 MuonIdProducer::~MuonIdProducer()
00167 {
00168 if (muIsoExtractorCalo_) delete muIsoExtractorCalo_;
00169 if (muIsoExtractorTrack_) delete muIsoExtractorTrack_;
00170 if (muIsoExtractorJet_) delete muIsoExtractorJet_;
00171 if (theTimingFiller_) delete theTimingFiller_;
00172 if (meshAlgo_) delete meshAlgo_;
00173
00174 }
00175
00176 void MuonIdProducer::init(edm::Event& iEvent, const edm::EventSetup& iSetup)
00177 {
00178
00179
00180
00181 innerTrackCollectionHandle_.clear();
00182 outerTrackCollectionHandle_.clear();
00183 linkCollectionHandle_.clear();
00184 muonCollectionHandle_.clear();
00185
00186 tpfmsCollectionHandle_.clear();
00187 pickyCollectionHandle_.clear();
00188 dytCollectionHandle_.clear();
00189
00190
00191
00192 edm::ESHandle<Propagator> propagator;
00193 iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
00194 trackAssociator_.setPropagator(propagator.product());
00195
00196 if (fillTrackerKink_) trackerKinkFinder_->init(iSetup);
00197
00198
00199 for ( unsigned int i = 0; i < inputCollectionLabels_.size(); ++i ) {
00200 if ( inputCollectionTypes_[i] == "inner tracks" ) {
00201 iEvent.getByLabel(inputCollectionLabels_[i], innerTrackCollectionHandle_);
00202 if (! innerTrackCollectionHandle_.isValid())
00203 throw cms::Exception("FatalError") << "Failed to get input track collection with label: " << inputCollectionLabels_[i];
00204 LogTrace("MuonIdentification") << "Number of input inner tracks: " << innerTrackCollectionHandle_->size();
00205 continue;
00206 }
00207 if ( inputCollectionTypes_[i] == "outer tracks" ) {
00208 iEvent.getByLabel(inputCollectionLabels_[i], outerTrackCollectionHandle_);
00209 if (! outerTrackCollectionHandle_.isValid())
00210 throw cms::Exception("FatalError") << "Failed to get input track collection with label: " << inputCollectionLabels_[i];
00211 LogTrace("MuonIdentification") << "Number of input outer tracks: " << outerTrackCollectionHandle_->size();
00212 continue;
00213 }
00214 if ( inputCollectionTypes_[i] == "links" ) {
00215 iEvent.getByLabel(inputCollectionLabels_[i], linkCollectionHandle_);
00216 if (! linkCollectionHandle_.isValid())
00217 throw cms::Exception("FatalError") << "Failed to get input link collection with label: " << inputCollectionLabels_[i];
00218 LogTrace("MuonIdentification") << "Number of input links: " << linkCollectionHandle_->size();
00219 continue;
00220 }
00221 if ( inputCollectionTypes_[i] == "muons" ) {
00222 iEvent.getByLabel(inputCollectionLabels_[i], muonCollectionHandle_);
00223 if (! muonCollectionHandle_.isValid())
00224 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputCollectionLabels_[i];
00225 LogTrace("MuonIdentification") << "Number of input muons: " << muonCollectionHandle_->size();
00226 continue;
00227 }
00228 if ( fillGlobalTrackRefits_ && inputCollectionTypes_[i] == "tev firstHit" ) {
00229 iEvent.getByLabel(inputCollectionLabels_[i], tpfmsCollectionHandle_);
00230 if (! tpfmsCollectionHandle_.isValid())
00231 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputCollectionLabels_[i];
00232 LogTrace("MuonIdentification") << "Number of input muons: " << tpfmsCollectionHandle_->size();
00233 continue;
00234 }
00235
00236 if ( fillGlobalTrackRefits_ && inputCollectionTypes_[i] == "tev picky" ) {
00237 iEvent.getByLabel(inputCollectionLabels_[i], pickyCollectionHandle_);
00238 if (! pickyCollectionHandle_.isValid())
00239 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputCollectionLabels_[i];
00240 LogTrace("MuonIdentification") << "Number of input muons: " << pickyCollectionHandle_->size();
00241 continue;
00242 }
00243
00244 if ( fillGlobalTrackRefits_ && inputCollectionTypes_[i] == "tev dyt" ) {
00245 iEvent.getByLabel(inputCollectionLabels_[i], dytCollectionHandle_);
00246 if (! dytCollectionHandle_.isValid())
00247 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputCollectionLabels_[i];
00248 LogTrace("MuonIdentification") << "Number of input muons: " << dytCollectionHandle_->size();
00249 continue;
00250 }
00251 throw cms::Exception("FatalError") << "Unknown input collection type: " << inputCollectionTypes_[i];
00252 }
00253 }
00254
00255 reco::Muon MuonIdProducer::makeMuon(edm::Event& iEvent, const edm::EventSetup& iSetup,
00256 const reco::TrackRef& track, MuonIdProducer::TrackType type)
00257 {
00258 LogTrace("MuonIdentification") << "Creating a muon from a track " << track.get()->pt() <<
00259 " Pt (GeV), eta: " << track.get()->eta();
00260 reco::Muon aMuon( makeMuon( *(track.get()) ) );
00261
00262 aMuon.setMuonTrack(type,track);
00263 aMuon.setBestTrack(type);
00264 aMuon.setTunePBestTrack(type);
00265
00266 return aMuon;
00267 }
00268
00269 reco::CaloMuon MuonIdProducer::makeCaloMuon( const reco::Muon& muon )
00270 {
00271 reco::CaloMuon aMuon;
00272 aMuon.setInnerTrack( muon.innerTrack() );
00273
00274 if (muon.isEnergyValid()) aMuon.setCalEnergy( muon.calEnergy() );
00275
00276 if (fillCaloCompatibility_) aMuon.setCaloCompatibility( muonCaloCompatibility_.evaluate(muon) );
00277 return aMuon;
00278 }
00279
00280
00281 reco::Muon MuonIdProducer::makeMuon( const reco::MuonTrackLinks& links )
00282 {
00283 LogTrace("MuonIdentification") << "Creating a muon from a link to tracks object";
00284
00285 reco::Muon aMuon;
00286 reco::Muon::MuonTrackTypePair chosenTrack;
00287 reco::TrackRef tpfmsRef;
00288 reco::TrackRef pickyRef;
00289 bool useSigmaSwitch = false;
00290
00291 if (tpfmsCollectionHandle_.isValid() && !tpfmsCollectionHandle_.failedToGet() &&
00292 pickyCollectionHandle_.isValid() && !pickyCollectionHandle_.failedToGet()) {
00293
00294 tpfmsRef = muon::getTevRefitTrack(links.globalTrack(), *tpfmsCollectionHandle_);
00295 pickyRef = muon::getTevRefitTrack(links.globalTrack(), *pickyCollectionHandle_);
00296
00297 if (tpfmsRef.isNull() && pickyRef.isNull()){
00298 edm::LogWarning("MakeMuonWithTEV")<<"Failed to get TEV refits, fall back to sigma switch.";
00299 useSigmaSwitch = true;
00300 }
00301 } else {
00302 useSigmaSwitch = true;
00303 }
00304
00305 if (useSigmaSwitch){
00306 chosenTrack = muon::sigmaSwitch( links.globalTrack(), links.trackerTrack(),
00307 sigmaThresholdToFillCandidateP4WithGlobalFit_,
00308 ptThresholdToFillCandidateP4WithGlobalFit_);
00309 } else {
00310 chosenTrack = muon::tevOptimized( links.globalTrack(), links.trackerTrack(),
00311 tpfmsRef, pickyRef,
00312 ptThresholdToFillCandidateP4WithGlobalFit_);
00313 }
00314 aMuon = makeMuon(*chosenTrack.first);
00315 aMuon.setInnerTrack( links.trackerTrack() );
00316 aMuon.setOuterTrack( links.standAloneTrack() );
00317 aMuon.setGlobalTrack( links.globalTrack() );
00318 aMuon.setBestTrack(chosenTrack.second);
00319 aMuon.setTunePBestTrack(chosenTrack.second);
00320
00321 if(fillGlobalTrackRefits_){
00322 if (tpfmsCollectionHandle_.isValid() && !tpfmsCollectionHandle_.failedToGet()) {
00323 reco::TrackToTrackMap::const_iterator it = tpfmsCollectionHandle_->find(links.globalTrack());
00324 if (it != tpfmsCollectionHandle_->end()) aMuon.setMuonTrack(reco::Muon::TPFMS, (it->val));
00325 }
00326 if (pickyCollectionHandle_.isValid() && !pickyCollectionHandle_.failedToGet()) {
00327 reco::TrackToTrackMap::const_iterator it = pickyCollectionHandle_->find(links.globalTrack());
00328 if (it != pickyCollectionHandle_->end()) aMuon.setMuonTrack(reco::Muon::Picky, (it->val));
00329 }
00330 if (dytCollectionHandle_.isValid() && !dytCollectionHandle_.failedToGet()) {
00331 reco::TrackToTrackMap::const_iterator it = dytCollectionHandle_->find(links.globalTrack());
00332 if (it != dytCollectionHandle_->end()) aMuon.setMuonTrack(reco::Muon::DYT, (it->val));
00333 }
00334 }
00335 return aMuon;
00336 }
00337
00338
00339 bool MuonIdProducer::isGoodTrack( const reco::Track& track )
00340 {
00341
00342 if (track.pt() < minPt_ || (track.p() < minP_ && track.p() < minPCaloMuon_)){
00343 LogTrace("MuonIdentification") << "Skipped low momentum track (Pt,P): " << track.pt() <<
00344 ", " << track.p() << " GeV";
00345 return false;
00346 }
00347
00348
00349 if ( fabs(track.eta()) > maxAbsEta_ ){
00350 LogTrace("MuonIdentification") << "Skipped track with large pseudo rapidity (Eta: " << track.eta() << " )";
00351 return false;
00352 }
00353 return true;
00354 }
00355
00356 unsigned int MuonIdProducer::chamberId( const DetId& id )
00357 {
00358 switch ( id.det() ) {
00359 case DetId::Muon:
00360 switch ( id.subdetId() ) {
00361 case MuonSubdetId::DT:
00362 {
00363 DTChamberId detId(id.rawId());
00364 return detId.rawId();
00365 }
00366 break;
00367 case MuonSubdetId::CSC:
00368 {
00369 CSCDetId detId(id.rawId());
00370 return detId.chamberId().rawId();
00371 }
00372 break;
00373 default:
00374 return 0;
00375 }
00376 default:
00377 return 0;
00378 }
00379 return 0;
00380 }
00381
00382
00383 int MuonIdProducer::overlap(const reco::Muon& muon, const reco::Track& track)
00384 {
00385 int numberOfCommonDetIds = 0;
00386 if ( ! muon.isMatchesValid() ||
00387 track.extra().isNull() ||
00388 track.extra()->recHits().isNull() ) return numberOfCommonDetIds;
00389 const std::vector<reco::MuonChamberMatch>& matches( muon.matches() );
00390 for ( std::vector<reco::MuonChamberMatch>::const_iterator match = matches.begin();
00391 match != matches.end(); ++match )
00392 {
00393 if ( match->segmentMatches.empty() ) continue;
00394 bool foundCommonDetId = false;
00395
00396 for ( TrackingRecHitRefVector::const_iterator hit = track.extra()->recHitsBegin();
00397 hit != track.extra()->recHitsEnd(); ++hit )
00398 {
00399
00400
00401
00402
00403 if ( chamberId(hit->get()->geographicalId()) == match->id.rawId() ) {
00404 foundCommonDetId = true;
00405 break;
00406 }
00407 }
00408 if ( foundCommonDetId ) {
00409 numberOfCommonDetIds++;
00410 break;
00411 }
00412 }
00413 return numberOfCommonDetIds;
00414 }
00415
00416 void MuonIdProducer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
00417 {
00418 edm::ESHandle<CSCGeometry> geomHandle;
00419 iSetup.get<MuonGeometryRecord>().get(geomHandle);
00420
00421 meshAlgo_->setCSCGeometry(geomHandle.product());
00422
00423 }
00424
00425 bool validateGlobalMuonPair( const reco::MuonTrackLinks& goodMuon,
00426 const reco::MuonTrackLinks& badMuon )
00427 {
00428 if ( std::min(goodMuon.globalTrack()->hitPattern().numberOfValidMuonHits(),
00429 badMuon.globalTrack()->hitPattern().numberOfValidMuonHits()) > 10 ){
00430 if ( goodMuon.globalTrack()->normalizedChi2() >
00431 badMuon.globalTrack()->normalizedChi2() )
00432 return false;
00433 else
00434 return true;
00435 }
00436 if ( goodMuon.globalTrack()->hitPattern().numberOfValidMuonHits() <
00437 badMuon.globalTrack()->hitPattern().numberOfValidMuonHits() ) return false;
00438 return true;
00439 }
00440
00441 void MuonIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00442 {
00443
00444
00445
00446 std::auto_ptr<reco::MuonCollection> outputMuons(new reco::MuonCollection);
00447 std::auto_ptr<reco::CaloMuonCollection> caloMuons( new reco::CaloMuonCollection );
00448
00449 init(iEvent, iSetup);
00450
00451 std::auto_ptr<reco::MuonTimeExtraMap> muonTimeMap(new reco::MuonTimeExtraMap());
00452 reco::MuonTimeExtraMap::Filler filler(*muonTimeMap);
00453 std::auto_ptr<reco::MuonTimeExtraMap> muonTimeMapDT(new reco::MuonTimeExtraMap());
00454 reco::MuonTimeExtraMap::Filler fillerDT(*muonTimeMapDT);
00455 std::auto_ptr<reco::MuonTimeExtraMap> muonTimeMapCSC(new reco::MuonTimeExtraMap());
00456 reco::MuonTimeExtraMap::Filler fillerCSC(*muonTimeMapCSC);
00457
00458 std::auto_ptr<reco::IsoDepositMap> trackDepMap(new reco::IsoDepositMap());
00459 reco::IsoDepositMap::Filler trackDepFiller(*trackDepMap);
00460 std::auto_ptr<reco::IsoDepositMap> ecalDepMap(new reco::IsoDepositMap());
00461 reco::IsoDepositMap::Filler ecalDepFiller(*ecalDepMap);
00462 std::auto_ptr<reco::IsoDepositMap> hcalDepMap(new reco::IsoDepositMap());
00463 reco::IsoDepositMap::Filler hcalDepFiller(*hcalDepMap);
00464 std::auto_ptr<reco::IsoDepositMap> hoDepMap(new reco::IsoDepositMap());
00465 reco::IsoDepositMap::Filler hoDepFiller(*hoDepMap);
00466 std::auto_ptr<reco::IsoDepositMap> jetDepMap(new reco::IsoDepositMap());
00467 reco::IsoDepositMap::Filler jetDepFiller(*jetDepMap);
00468
00469
00470
00471
00472 if ( muonCollectionHandle_.isValid() )
00473 for ( reco::MuonCollection::const_iterator muon = muonCollectionHandle_->begin();
00474 muon != muonCollectionHandle_->end(); ++muon )
00475 outputMuons->push_back(*muon);
00476
00477
00478 if ( linkCollectionHandle_.isValid() ){
00479 std::vector<bool> goodmuons(linkCollectionHandle_->size(),true);
00480 if ( goodmuons.size()>1 ){
00481
00482 for ( unsigned int i=0; i<linkCollectionHandle_->size()-1; ++i ){
00483 if (!checkLinks(&linkCollectionHandle_->at(i))) continue;
00484 for ( unsigned int j=i+1; j<linkCollectionHandle_->size(); ++j ){
00485 if (!checkLinks(&linkCollectionHandle_->at(j))) continue;
00486 if ( linkCollectionHandle_->at(i).trackerTrack().isNonnull() &&
00487 linkCollectionHandle_->at(i).trackerTrack() ==
00488 linkCollectionHandle_->at(j).trackerTrack() )
00489 {
00490
00491
00492
00493
00494 if ( validateGlobalMuonPair(linkCollectionHandle_->at(i),linkCollectionHandle_->at(j)) )
00495 goodmuons[j] = false;
00496 else
00497 goodmuons[i] = false;
00498 }
00499 }
00500 }
00501
00502 for ( unsigned int i=0; i<linkCollectionHandle_->size()-1; ++i ){
00503 if ( !goodmuons[i] ) continue;
00504 if (!checkLinks(&linkCollectionHandle_->at(i))) continue;
00505 for ( unsigned int j=i+1; j<linkCollectionHandle_->size(); ++j ){
00506 if ( !goodmuons[j] ) continue;
00507 if (!checkLinks(&linkCollectionHandle_->at(j))) continue;
00508 if ( linkCollectionHandle_->at(i).standAloneTrack().isNonnull() &&
00509 linkCollectionHandle_->at(i).standAloneTrack() ==
00510 linkCollectionHandle_->at(j).standAloneTrack() )
00511 {
00512 if ( validateGlobalMuonPair(linkCollectionHandle_->at(i),linkCollectionHandle_->at(j)) )
00513 goodmuons[j] = false;
00514 else
00515 goodmuons[i] = false;
00516 }
00517 }
00518 }
00519 }
00520 for ( unsigned int i=0; i<linkCollectionHandle_->size(); ++i ){
00521 if ( !goodmuons[i] ) continue;
00522 const reco::MuonTrackLinks* links = &linkCollectionHandle_->at(i);
00523 if ( ! checkLinks(links)) continue;
00524
00525 bool newMuon = true;
00526 for ( reco::MuonCollection::const_iterator muon = outputMuons->begin();
00527 muon != outputMuons->end(); ++muon )
00528 if ( muon->track() == links->trackerTrack() &&
00529 muon->standAloneMuon() == links->standAloneTrack() &&
00530 muon->combinedMuon() == links->globalTrack() )
00531 newMuon = false;
00532 if ( newMuon ) {
00533 outputMuons->push_back( makeMuon( *links ) );
00534 outputMuons->back().setType(reco::Muon::GlobalMuon | reco::Muon::StandAloneMuon);
00535 }
00536 }
00537 }
00538
00539
00540 if ( innerTrackCollectionHandle_.isValid() ) {
00541 LogTrace("MuonIdentification") << "Creating tracker muons";
00542 for ( unsigned int i = 0; i < innerTrackCollectionHandle_->size(); ++i )
00543 {
00544 const reco::Track& track = innerTrackCollectionHandle_->at(i);
00545 if ( ! isGoodTrack( track ) ) continue;
00546 bool splitTrack = false;
00547 if ( track.extra().isAvailable() &&
00548 TrackDetectorAssociator::crossedIP( track ) ) splitTrack = true;
00549 std::vector<TrackDetectorAssociator::Direction> directions;
00550 if ( splitTrack ) {
00551 directions.push_back(TrackDetectorAssociator::InsideOut);
00552 directions.push_back(TrackDetectorAssociator::OutsideIn);
00553 } else {
00554 directions.push_back(TrackDetectorAssociator::Any);
00555 }
00556 for ( std::vector<TrackDetectorAssociator::Direction>::const_iterator direction = directions.begin();
00557 direction != directions.end(); ++direction )
00558 {
00559
00560
00561 reco::Muon trackerMuon( makeMuon(iEvent, iSetup, reco::TrackRef( innerTrackCollectionHandle_, i ), reco::Muon::InnerTrack ) );
00562 trackerMuon.setType( reco::Muon::TrackerMuon | reco::Muon::RPCMuon );
00563 fillMuonId(iEvent, iSetup, trackerMuon, *direction);
00564
00565
00566 if ( debugWithTruthMatching_ ) {
00567
00568
00569 MuonIdTruthInfo::truthMatchMuon(iEvent, iSetup, trackerMuon);
00570 }
00571
00572
00573
00574
00575 bool newMuon = true;
00576 bool goodTrackerMuon = isGoodTrackerMuon( trackerMuon );
00577 bool goodRPCMuon = isGoodRPCMuon( trackerMuon );
00578 for ( reco::MuonCollection::iterator muon = outputMuons->begin();
00579 muon != outputMuons->end(); ++muon )
00580 {
00581 if ( muon->innerTrack().get() == trackerMuon.innerTrack().get() &&
00582 cos(phiOfMuonIneteractionRegion(*muon) -
00583 phiOfMuonIneteractionRegion(trackerMuon)) > 0 )
00584 {
00585 newMuon = false;
00586 muon->setMatches( trackerMuon.matches() );
00587 if (trackerMuon.isTimeValid()) muon->setTime( trackerMuon.time() );
00588 if (trackerMuon.isEnergyValid()) muon->setCalEnergy( trackerMuon.calEnergy() );
00589 if (goodTrackerMuon) muon->setType( muon->type() | reco::Muon::TrackerMuon );
00590 if (goodRPCMuon) muon->setType( muon->type() | reco::Muon::RPCMuon );
00591 LogTrace("MuonIdentification") << "Found a corresponding global muon. Set energy, matches and move on";
00592 break;
00593 }
00594 }
00595 if ( newMuon ) {
00596 if ( goodTrackerMuon ){
00597 outputMuons->push_back( trackerMuon );
00598 } else {
00599 LogTrace("MuonIdentification") << "track failed minimal number of muon matches requirement";
00600 const reco::CaloMuon& caloMuon = makeCaloMuon(trackerMuon);
00601 if ( ! caloMuon.isCaloCompatibilityValid() || caloMuon.caloCompatibility() < caloCut_ || caloMuon.p() < minPCaloMuon_) continue;
00602 caloMuons->push_back( caloMuon );
00603 }
00604 }
00605 }
00606 }
00607 }
00608
00609
00610 if ( outerTrackCollectionHandle_.isValid() ) {
00611 LogTrace("MuonIdentification") << "Looking for new muons among stand alone muon tracks";
00612 for ( unsigned int i = 0; i < outerTrackCollectionHandle_->size(); ++i )
00613 {
00614
00615 bool newMuon = true;
00616 for ( reco::MuonCollection::iterator muon = outputMuons->begin();
00617 muon != outputMuons->end(); ++muon )
00618 {
00619 if ( ! muon->standAloneMuon().isNull() ) {
00620
00621 if ( muon->standAloneMuon().get() == &(outerTrackCollectionHandle_->at(i)) ) {
00622 newMuon = false;
00623 break;
00624 }
00625 } else {
00626
00627
00628
00629
00630
00631
00632
00633 if ( overlap(*muon,outerTrackCollectionHandle_->at(i))>0 ) {
00634 LogTrace("MuonIdentification") << "Found associated tracker muon. Set a reference and move on";
00635 newMuon = false;
00636 muon->setOuterTrack( reco::TrackRef( outerTrackCollectionHandle_, i ) );
00637 muon->setType( muon->type() | reco::Muon::StandAloneMuon );
00638 break;
00639 }
00640 }
00641 }
00642 if ( newMuon ) {
00643 LogTrace("MuonIdentification") << "No associated stand alone track is found. Making a muon";
00644 outputMuons->push_back( makeMuon(iEvent, iSetup,
00645 reco::TrackRef( outerTrackCollectionHandle_, i ), reco::Muon::OuterTrack ) );
00646 outputMuons->back().setType( reco::Muon::StandAloneMuon );
00647 }
00648 }
00649 }
00650
00651 LogTrace("MuonIdentification") << "Dress up muons if it's necessary";
00652
00653 int nMuons=outputMuons->size();
00654
00655 std::vector<reco::MuonTimeExtra> dtTimeColl(nMuons);
00656 std::vector<reco::MuonTimeExtra> cscTimeColl(nMuons);
00657 std::vector<reco::MuonTimeExtra> combinedTimeColl(nMuons);
00658 std::vector<reco::IsoDeposit> trackDepColl(nMuons);
00659 std::vector<reco::IsoDeposit> ecalDepColl(nMuons);
00660 std::vector<reco::IsoDeposit> hcalDepColl(nMuons);
00661 std::vector<reco::IsoDeposit> hoDepColl(nMuons);
00662 std::vector<reco::IsoDeposit> jetDepColl(nMuons);
00663
00664
00665 unsigned int i=0;
00666 for ( reco::MuonCollection::iterator muon = outputMuons->begin(); muon != outputMuons->end(); ++muon )
00667 {
00668
00669
00670 if ( ( fillMatching_ && ! muon->isMatchesValid() ) ||
00671 ( fillEnergy_ && !muon->isEnergyValid() ) )
00672 {
00673
00674
00675 if ( muon->isStandAloneMuon() ) {
00676 if ( cos(phiOfMuonIneteractionRegion(*muon) - muon->phi()) > 0 )
00677 fillMuonId(iEvent, iSetup, *muon, TrackDetectorAssociator::InsideOut);
00678 else
00679 fillMuonId(iEvent, iSetup, *muon, TrackDetectorAssociator::OutsideIn);
00680 } else {
00681 LogTrace("MuonIdentification") << "THIS SHOULD NEVER HAPPEN";
00682 fillMuonId(iEvent, iSetup, *muon);
00683 }
00684 }
00685
00686 if (fillGlobalTrackQuality_){
00687
00688 fillGlbQuality(iEvent, iSetup, *muon);
00689 }
00690 LogDebug("MuonIdentification");
00691
00692 if (fillTrackerKink_) {
00693 fillTrackerKink(*muon);
00694 }
00695
00696
00697 if ( fillCaloCompatibility_ ) muon->setCaloCompatibility( muonCaloCompatibility_.evaluate(*muon) );
00698
00699
00700
00701 if ( fillIsolation_ ) fillMuonIsolation(iEvent, iSetup, *muon,
00702 trackDepColl[i], ecalDepColl[i], hcalDepColl[i], hoDepColl[i], jetDepColl[i]);
00703
00704
00705
00706 reco::MuonTime muonTime;
00707 reco::MuonTimeExtra dtTime;
00708 reco::MuonTimeExtra cscTime;
00709 reco::MuonTimeExtra combinedTime;
00710
00711 theTimingFiller_->fillTiming(*muon, dtTime, cscTime, combinedTime, iEvent, iSetup);
00712
00713 muonTime.nDof=combinedTime.nDof();
00714 muonTime.timeAtIpInOut=combinedTime.timeAtIpInOut();
00715 muonTime.timeAtIpInOutErr=combinedTime.timeAtIpInOutErr();
00716 muonTime.timeAtIpOutIn=combinedTime.timeAtIpOutIn();
00717 muonTime.timeAtIpOutInErr=combinedTime.timeAtIpOutInErr();
00718
00719 muon->setTime( muonTime);
00720 dtTimeColl[i] = dtTime;
00721 cscTimeColl[i] = cscTime;
00722 combinedTimeColl[i] = combinedTime;
00723
00724 i++;
00725
00726 }
00727
00728 LogTrace("MuonIdentification") << "number of muons produced: " << outputMuons->size();
00729
00730 if ( fillMatching_ ) fillArbitrationInfo( outputMuons.get() );
00731
00732 edm::OrphanHandle<reco::MuonCollection> muonHandle = iEvent.put(outputMuons);
00733
00734 filler.insert(muonHandle, combinedTimeColl.begin(), combinedTimeColl.end());
00735 filler.fill();
00736 fillerDT.insert(muonHandle, dtTimeColl.begin(), dtTimeColl.end());
00737 fillerDT.fill();
00738 fillerCSC.insert(muonHandle, cscTimeColl.begin(), cscTimeColl.end());
00739 fillerCSC.fill();
00740
00741 iEvent.put(muonTimeMap,"combined");
00742 iEvent.put(muonTimeMapDT,"dt");
00743 iEvent.put(muonTimeMapCSC,"csc");
00744
00745 if (writeIsoDeposits_ && fillIsolation_){
00746 trackDepFiller.insert(muonHandle, trackDepColl.begin(), trackDepColl.end());
00747 trackDepFiller.fill();
00748 iEvent.put(trackDepMap, trackDepositName_);
00749 ecalDepFiller.insert(muonHandle, ecalDepColl.begin(), ecalDepColl.end());
00750 ecalDepFiller.fill();
00751 iEvent.put(ecalDepMap, ecalDepositName_);
00752 hcalDepFiller.insert(muonHandle, hcalDepColl.begin(), hcalDepColl.end());
00753 hcalDepFiller.fill();
00754 iEvent.put(hcalDepMap, hcalDepositName_);
00755 hoDepFiller.insert(muonHandle, hoDepColl.begin(), hoDepColl.end());
00756 hoDepFiller.fill();
00757 iEvent.put(hoDepMap, hoDepositName_);
00758 jetDepFiller.insert(muonHandle, jetDepColl.begin(), jetDepColl.end());
00759 jetDepFiller.fill();
00760 iEvent.put(jetDepMap, jetDepositName_);
00761 }
00762
00763 iEvent.put(caloMuons);
00764 }
00765
00766
00767 bool MuonIdProducer::isGoodTrackerMuon( const reco::Muon& muon )
00768 {
00769 if(muon.track()->pt() < minPt_ || muon.track()->p() < minP_) return false;
00770 if ( addExtraSoftMuons_ &&
00771 muon.pt()<5 && fabs(muon.eta())<1.5 &&
00772 muon.numberOfMatches( reco::Muon::NoArbitration ) >= 1 ) return true;
00773 return ( muon.numberOfMatches( reco::Muon::NoArbitration ) >= minNumberOfMatches_ );
00774 }
00775
00776 bool MuonIdProducer::isGoodRPCMuon( const reco::Muon& muon )
00777 {
00778 if(muon.track()->pt() < minPt_ || muon.track()->p() < minP_) return false;
00779 if ( addExtraSoftMuons_ &&
00780 muon.pt()<5 && fabs(muon.eta())<1.5 &&
00781 muon.numberOfMatchedRPCLayers( reco::Muon::RPCHitAndTrackArbitration ) > 1 ) return true;
00782 return ( muon.numberOfMatchedRPCLayers( reco::Muon::RPCHitAndTrackArbitration ) > minNumberOfMatches_ );
00783 }
00784
00785 void MuonIdProducer::fillMuonId(edm::Event& iEvent, const edm::EventSetup& iSetup,
00786 reco::Muon& aMuon,
00787 TrackDetectorAssociator::Direction direction)
00788 {
00789
00790 const reco::Track* track = 0;
00791 if ( ! aMuon.track().isNull() )
00792 track = aMuon.track().get();
00793 else
00794 {
00795 if ( ! aMuon.standAloneMuon().isNull() )
00796 track = aMuon.standAloneMuon().get();
00797 else
00798 throw cms::Exception("FatalError") << "Failed to fill muon id information for a muon with undefined references to tracks";
00799 }
00800
00801 TrackDetMatchInfo info = trackAssociator_.associate(iEvent, iSetup, *track, parameters_, direction);
00802
00803 if ( fillEnergy_ ) {
00804 reco::MuonEnergy muonEnergy;
00805 muonEnergy.em = info.crossedEnergy(TrackDetMatchInfo::EcalRecHits);
00806 muonEnergy.had = info.crossedEnergy(TrackDetMatchInfo::HcalRecHits);
00807 muonEnergy.ho = info.crossedEnergy(TrackDetMatchInfo::HORecHits);
00808 muonEnergy.tower = info.crossedEnergy(TrackDetMatchInfo::TowerTotal);
00809 muonEnergy.emS9 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits,1);
00810 muonEnergy.emS25 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits,2);
00811 muonEnergy.hadS9 = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits,1);
00812 muonEnergy.hoS9 = info.nXnEnergy(TrackDetMatchInfo::HORecHits,1);
00813 muonEnergy.towerS9 = info.nXnEnergy(TrackDetMatchInfo::TowerTotal,1);
00814 muonEnergy.ecal_position = info.trkGlobPosAtEcal;
00815 muonEnergy.hcal_position = info.trkGlobPosAtHcal;
00816 if (! info.crossedEcalIds.empty() ) muonEnergy.ecal_id = info.crossedEcalIds.front();
00817 if (! info.crossedHcalIds.empty() ) muonEnergy.hcal_id = info.crossedHcalIds.front();
00818
00819 DetId emMaxId = info.findMaxDeposition(TrackDetMatchInfo::EcalRecHits,2);
00820 for(std::vector<const EcalRecHit*>::const_iterator hit=info.ecalRecHits.begin();
00821 hit!=info.ecalRecHits.end(); ++hit) {
00822 if ((*hit)->id() != emMaxId) continue;
00823 muonEnergy.emMax = (*hit)->energy();
00824 muonEnergy.ecal_time = (*hit)->time();
00825 }
00826 DetId hadMaxId = info.findMaxDeposition(TrackDetMatchInfo::HcalRecHits,1);
00827 for(std::vector<const HBHERecHit*>::const_iterator hit=info.hcalRecHits.begin();
00828 hit!=info.hcalRecHits.end(); ++hit) {
00829 if ((*hit)->id() != hadMaxId) continue;
00830 muonEnergy.hadMax = (*hit)->energy();
00831 muonEnergy.hcal_time = (*hit)->time();
00832 }
00833 aMuon.setCalEnergy( muonEnergy );
00834 }
00835 if ( ! fillMatching_ && ! aMuon.isTrackerMuon() && ! aMuon.isRPCMuon() ) return;
00836
00837 edm::Handle<RPCRecHitCollection> rpcRecHits;
00838 iEvent.getByLabel(edm::InputTag("rpcRecHits"), rpcRecHits);
00839
00840
00841 std::vector<reco::MuonChamberMatch> muonChamberMatches;
00842 unsigned int nubmerOfMatchesAccordingToTrackAssociator = 0;
00843 for( std::vector<TAMuonChamberMatch>::const_iterator chamber=info.chambers.begin();
00844 chamber!=info.chambers.end(); chamber++ )
00845 {
00846 if (chamber->id.subdetId() == 3 && rpcRecHits.isValid() ) continue;
00847 reco::MuonChamberMatch matchedChamber;
00848
00849 LocalError localError = chamber->tState.localError().positionError();
00850 matchedChamber.x = chamber->tState.localPosition().x();
00851 matchedChamber.y = chamber->tState.localPosition().y();
00852 matchedChamber.xErr = sqrt( localError.xx() );
00853 matchedChamber.yErr = sqrt( localError.yy() );
00854
00855 matchedChamber.dXdZ = chamber->tState.localDirection().z()!=0?chamber->tState.localDirection().x()/chamber->tState.localDirection().z():9999;
00856 matchedChamber.dYdZ = chamber->tState.localDirection().z()!=0?chamber->tState.localDirection().y()/chamber->tState.localDirection().z():9999;
00857
00858 AlgebraicSymMatrix55 trajectoryCovMatrix = chamber->tState.localError().matrix();
00859 matchedChamber.dXdZErr = trajectoryCovMatrix(1,1)>0?sqrt(trajectoryCovMatrix(1,1)):0;
00860 matchedChamber.dYdZErr = trajectoryCovMatrix(2,2)>0?sqrt(trajectoryCovMatrix(2,2)):0;
00861
00862 matchedChamber.edgeX = chamber->localDistanceX;
00863 matchedChamber.edgeY = chamber->localDistanceY;
00864
00865 matchedChamber.id = chamber->id;
00866 if ( ! chamber->segments.empty() ) ++nubmerOfMatchesAccordingToTrackAssociator;
00867
00868
00869 for( std::vector<TAMuonSegmentMatch>::const_iterator segment = chamber->segments.begin();
00870 segment != chamber->segments.end(); segment++ )
00871 {
00872 reco::MuonSegmentMatch matchedSegment;
00873 matchedSegment.x = segment->segmentLocalPosition.x();
00874 matchedSegment.y = segment->segmentLocalPosition.y();
00875 matchedSegment.dXdZ = segment->segmentLocalDirection.z()?segment->segmentLocalDirection.x()/segment->segmentLocalDirection.z():0;
00876 matchedSegment.dYdZ = segment->segmentLocalDirection.z()?segment->segmentLocalDirection.y()/segment->segmentLocalDirection.z():0;
00877 matchedSegment.xErr = segment->segmentLocalErrorXX>0?sqrt(segment->segmentLocalErrorXX):0;
00878 matchedSegment.yErr = segment->segmentLocalErrorYY>0?sqrt(segment->segmentLocalErrorYY):0;
00879 matchedSegment.dXdZErr = segment->segmentLocalErrorDxDz>0?sqrt(segment->segmentLocalErrorDxDz):0;
00880 matchedSegment.dYdZErr = segment->segmentLocalErrorDyDz>0?sqrt(segment->segmentLocalErrorDyDz):0;
00881 matchedSegment.t0 = segment->t0;
00882 matchedSegment.mask = 0;
00883 matchedSegment.dtSegmentRef = segment->dtSegmentRef;
00884 matchedSegment.cscSegmentRef = segment->cscSegmentRef;
00885 matchedSegment.hasZed_ = segment->hasZed;
00886 matchedSegment.hasPhi_ = segment->hasPhi;
00887
00888 bool matchedX = false;
00889 bool matchedY = false;
00890 LogTrace("MuonIdentification") << " matching local x, segment x: " << matchedSegment.x <<
00891 ", chamber x: " << matchedChamber.x << ", max: " << maxAbsDx_;
00892 LogTrace("MuonIdentification") << " matching local y, segment y: " << matchedSegment.y <<
00893 ", chamber y: " << matchedChamber.y << ", max: " << maxAbsDy_;
00894 if (matchedSegment.xErr>0 && matchedChamber.xErr>0 )
00895 LogTrace("MuonIdentification") << " xpull: " <<
00896 fabs(matchedSegment.x - matchedChamber.x)/sqrt(pow(matchedSegment.xErr,2) + pow(matchedChamber.xErr,2));
00897 if (matchedSegment.yErr>0 && matchedChamber.yErr>0 )
00898 LogTrace("MuonIdentification") << " ypull: " <<
00899 fabs(matchedSegment.y - matchedChamber.y)/sqrt(pow(matchedSegment.yErr,2) + pow(matchedChamber.yErr,2));
00900
00901 if (fabs(matchedSegment.x - matchedChamber.x) < maxAbsDx_) matchedX = true;
00902 if (fabs(matchedSegment.y - matchedChamber.y) < maxAbsDy_) matchedY = true;
00903 if (matchedSegment.xErr>0 && matchedChamber.xErr>0 &&
00904 fabs(matchedSegment.x - matchedChamber.x)/sqrt(pow(matchedSegment.xErr,2) + pow(matchedChamber.xErr,2)) < maxAbsPullX_) matchedX = true;
00905 if (matchedSegment.yErr>0 && matchedChamber.yErr>0 &&
00906 fabs(matchedSegment.y - matchedChamber.y)/sqrt(pow(matchedSegment.yErr,2) + pow(matchedChamber.yErr,2)) < maxAbsPullY_) matchedY = true;
00907 if (matchedX && matchedY) matchedChamber.segmentMatches.push_back(matchedSegment);
00908 }
00909 muonChamberMatches.push_back(matchedChamber);
00910 }
00911
00912
00913 if ( rpcRecHits.isValid() )
00914 {
00915
00916 for( std::vector<TAMuonChamberMatch>::const_iterator chamber=info.chambers.begin();
00917 chamber!=info.chambers.end(); chamber++ )
00918 {
00919
00920 if ( chamber->id.subdetId() != 3 ) continue;
00921
00922 reco::MuonChamberMatch matchedChamber;
00923
00924 LocalError localError = chamber->tState.localError().positionError();
00925 matchedChamber.x = chamber->tState.localPosition().x();
00926 matchedChamber.y = chamber->tState.localPosition().y();
00927 matchedChamber.xErr = sqrt( localError.xx() );
00928 matchedChamber.yErr = sqrt( localError.yy() );
00929
00930 matchedChamber.dXdZ = chamber->tState.localDirection().z()!=0?chamber->tState.localDirection().x()/chamber->tState.localDirection().z():9999;
00931 matchedChamber.dYdZ = chamber->tState.localDirection().z()!=0?chamber->tState.localDirection().y()/chamber->tState.localDirection().z():9999;
00932
00933 AlgebraicSymMatrix55 trajectoryCovMatrix = chamber->tState.localError().matrix();
00934 matchedChamber.dXdZErr = trajectoryCovMatrix(1,1)>0?sqrt(trajectoryCovMatrix(1,1)):0;
00935 matchedChamber.dYdZErr = trajectoryCovMatrix(2,2)>0?sqrt(trajectoryCovMatrix(2,2)):0;
00936
00937 matchedChamber.edgeX = chamber->localDistanceX;
00938 matchedChamber.edgeY = chamber->localDistanceY;
00939
00940 matchedChamber.id = chamber->id;
00941
00942 for ( RPCRecHitCollection::const_iterator rpcRecHit = rpcRecHits->begin();
00943 rpcRecHit != rpcRecHits->end(); ++rpcRecHit )
00944 {
00945 reco::MuonRPCHitMatch rpcHitMatch;
00946
00947 if ( rpcRecHit->rawId() != chamber->id.rawId() ) continue;
00948
00949 rpcHitMatch.x = rpcRecHit->localPosition().x();
00950 rpcHitMatch.mask = 0;
00951 rpcHitMatch.bx = rpcRecHit->BunchX();
00952
00953 const double AbsDx = fabs(rpcRecHit->localPosition().x()-chamber->tState.localPosition().x());
00954 if( AbsDx <= 20 or AbsDx/sqrt(localError.xx()) <= 4 ) matchedChamber.rpcMatches.push_back(rpcHitMatch);
00955 }
00956
00957 muonChamberMatches.push_back(matchedChamber);
00958 }
00959 }
00960
00961 aMuon.setMatches(muonChamberMatches);
00962
00963 LogTrace("MuonIdentification") << "number of muon chambers: " << aMuon.matches().size() << "\n"
00964 << "number of chambers with segments according to the associator requirements: " <<
00965 nubmerOfMatchesAccordingToTrackAssociator;
00966 LogTrace("MuonIdentification") << "number of segment matches with the producer requirements: " <<
00967 aMuon.numberOfMatches( reco::Muon::NoArbitration );
00968
00969
00970 }
00971
00972 void MuonIdProducer::fillArbitrationInfo( reco::MuonCollection* pOutputMuons )
00973 {
00974
00975
00976
00977 std::vector<std::pair<reco::MuonChamberMatch*,reco::MuonSegmentMatch*> > chamberPairs;
00978 std::vector<std::pair<reco::MuonChamberMatch*,reco::MuonSegmentMatch*> > stationPairs;
00979 std::vector<std::pair<reco::MuonChamberMatch*,reco::MuonSegmentMatch*> > arbitrationPairs;
00980
00981
00982 for( unsigned int muonIndex1 = 0; muonIndex1 < pOutputMuons->size(); ++muonIndex1 )
00983 {
00984
00985 for( std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = pOutputMuons->at(muonIndex1).matches().begin();
00986 chamberIter1 != pOutputMuons->at(muonIndex1).matches().end(); ++chamberIter1 )
00987 {
00988 if(chamberIter1->segmentMatches.empty()) continue;
00989 chamberPairs.clear();
00990
00991
00992 for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
00993 segmentIter1 != chamberIter1->segmentMatches.end(); ++segmentIter1 )
00994 {
00995 chamberPairs.push_back(std::make_pair(&(*chamberIter1), &(*segmentIter1)));
00996 if(!segmentIter1->isMask())
00997 {
00998 arbitrationPairs.clear();
00999 arbitrationPairs.push_back(std::make_pair(&(*chamberIter1), &(*segmentIter1)));
01000
01001
01002
01003
01004
01005 if (pOutputMuons->at(muonIndex1).isTrackerMuon()) {
01006
01007 for( unsigned int muonIndex2 = muonIndex1+1; muonIndex2 < pOutputMuons->size(); ++muonIndex2 )
01008 {
01009
01010 if (! pOutputMuons->at(muonIndex2).isTrackerMuon()) continue;
01011
01012 for( std::vector<reco::MuonChamberMatch>::iterator chamberIter2 = pOutputMuons->at(muonIndex2).matches().begin();
01013 chamberIter2 != pOutputMuons->at(muonIndex2).matches().end(); ++chamberIter2 )
01014 {
01015
01016 for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter2 = chamberIter2->segmentMatches.begin();
01017 segmentIter2 != chamberIter2->segmentMatches.end(); ++segmentIter2 )
01018 {
01019 if(segmentIter2->isMask()) continue;
01020 if(fabs(segmentIter2->x - segmentIter1->x ) < 1E-3 &&
01021 fabs(segmentIter2->y - segmentIter1->y ) < 1E-3 &&
01022 fabs(segmentIter2->dXdZ - segmentIter1->dXdZ ) < 1E-3 &&
01023 fabs(segmentIter2->dYdZ - segmentIter1->dYdZ ) < 1E-3 &&
01024 fabs(segmentIter2->xErr - segmentIter1->xErr ) < 1E-3 &&
01025 fabs(segmentIter2->yErr - segmentIter1->yErr ) < 1E-3 &&
01026 fabs(segmentIter2->dXdZErr - segmentIter1->dXdZErr) < 1E-3 &&
01027 fabs(segmentIter2->dYdZErr - segmentIter1->dYdZErr) < 1E-3)
01028 arbitrationPairs.push_back(std::make_pair(&(*chamberIter2), &(*segmentIter2)));
01029 }
01030 }
01031 }
01032 }
01033
01034
01035 if(arbitrationPairs.empty()) continue;
01036 if(arbitrationPairs.size()==1) {
01037 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDRSlope);
01038 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDXSlope);
01039 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDR);
01040 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDX);
01041 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::Arbitrated);
01042 } else {
01043 sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDRSlope));
01044 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDRSlope);
01045 sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDXSlope));
01046 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDXSlope);
01047 sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDR));
01048 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDR);
01049 sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDX));
01050 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDX);
01051 for( unsigned int it = 0; it < arbitrationPairs.size(); ++it )
01052 arbitrationPairs.at(it).second->setMask(reco::MuonSegmentMatch::Arbitrated);
01053 }
01054 }
01055
01056
01057 if( chamberIter1->id.subdetId() == MuonSubdetId::CSC && arbClean_ ) {
01058 CSCDetId cscid(chamberIter1->id);
01059 if(cscid.ring() == 4)
01060 for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter2 = chamberIter1->segmentMatches.begin();
01061 segmentIter2 != chamberIter1->segmentMatches.end(); ++segmentIter2 ) {
01062 if( segmentIter1->cscSegmentRef.isNonnull() && segmentIter2->cscSegmentRef.isNonnull() )
01063 if( meshAlgo_->isDuplicateOf(segmentIter1->cscSegmentRef,segmentIter2->cscSegmentRef) &&
01064 (segmentIter2->mask & 0x1e0000) &&
01065 (segmentIter1->mask & 0x1e0000) )
01066 segmentIter2->setMask(reco::MuonSegmentMatch::BelongsToTrackByME1aClean);
01067
01068 }
01069 }
01070
01071 }
01072
01073
01074 if(chamberPairs.empty()) continue;
01075 if(chamberPairs.size()==1) {
01076 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDRSlope);
01077 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDXSlope);
01078 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDR);
01079 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDX);
01080 } else {
01081 sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDRSlope));
01082 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDRSlope);
01083 sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDXSlope));
01084 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDXSlope);
01085 sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDR));
01086 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDR);
01087 sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDX));
01088 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDX);
01089 }
01090 }
01091
01092
01093 for( int stationIndex = 1; stationIndex < 5; ++stationIndex )
01094 for( int detectorIndex = 1; detectorIndex < 4; ++detectorIndex )
01095 {
01096 stationPairs.clear();
01097
01098
01099 for( std::vector<reco::MuonChamberMatch>::iterator chamberIter = pOutputMuons->at(muonIndex1).matches().begin();
01100 chamberIter != pOutputMuons->at(muonIndex1).matches().end(); ++chamberIter )
01101 {
01102 if(!(chamberIter->station()==stationIndex && chamberIter->detector()==detectorIndex)) continue;
01103 if(chamberIter->segmentMatches.empty()) continue;
01104
01105 for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter = chamberIter->segmentMatches.begin();
01106 segmentIter != chamberIter->segmentMatches.end(); ++segmentIter )
01107 stationPairs.push_back(std::make_pair(&(*chamberIter), &(*segmentIter)));
01108 }
01109
01110 if(stationPairs.empty()) continue;
01111 if(stationPairs.size()==1) {
01112 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDRSlope);
01113 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDXSlope);
01114 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDR);
01115 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDX);
01116 } else {
01117 sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDRSlope));
01118 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDRSlope);
01119 sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDXSlope));
01120 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDXSlope);
01121 sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDR));
01122 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDR);
01123 sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDX));
01124 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDX);
01125 }
01126 }
01127
01128 }
01129
01130 if(arbClean_) {
01131
01132 meshAlgo_->clearMesh();
01133 meshAlgo_->runMesh(pOutputMuons);
01134 }
01135 }
01136
01137 void MuonIdProducer::fillMuonIsolation(edm::Event& iEvent, const edm::EventSetup& iSetup, reco::Muon& aMuon,
01138 reco::IsoDeposit& trackDep, reco::IsoDeposit& ecalDep, reco::IsoDeposit& hcalDep, reco::IsoDeposit& hoDep,
01139 reco::IsoDeposit& jetDep)
01140 {
01141 reco::MuonIsolation isoR03, isoR05;
01142 const reco::Track* track = 0;
01143 if ( ! aMuon.track().isNull() )
01144 track = aMuon.track().get();
01145 else
01146 {
01147 if ( ! aMuon.standAloneMuon().isNull() )
01148 track = aMuon.standAloneMuon().get();
01149 else
01150 throw cms::Exception("FatalError") << "Failed to compute muon isolation information for a muon with undefined references to tracks";
01151 }
01152
01153
01154 reco::IsoDeposit depTrk = muIsoExtractorTrack_->deposit(iEvent, iSetup, *track );
01155 std::vector<reco::IsoDeposit> caloDeps = muIsoExtractorCalo_->deposits(iEvent, iSetup, *track);
01156 reco::IsoDeposit depJet = muIsoExtractorJet_->deposit(iEvent, iSetup, *track );
01157
01158 if(caloDeps.size()!=3) {
01159 LogTrace("MuonIdentification") << "Failed to fill vector of calorimeter isolation deposits!";
01160 return;
01161 }
01162
01163 reco::IsoDeposit depEcal = caloDeps.at(0);
01164 reco::IsoDeposit depHcal = caloDeps.at(1);
01165 reco::IsoDeposit depHo = caloDeps.at(2);
01166
01167 trackDep = depTrk;
01168 ecalDep = depEcal;
01169 hcalDep = depHcal;
01170 hoDep = depHo;
01171 jetDep = depJet;
01172
01173 isoR03.sumPt = depTrk.depositWithin(0.3);
01174 isoR03.emEt = depEcal.depositWithin(0.3);
01175 isoR03.hadEt = depHcal.depositWithin(0.3);
01176 isoR03.hoEt = depHo.depositWithin(0.3);
01177 isoR03.nTracks = depTrk.depositAndCountWithin(0.3).second;
01178 isoR03.nJets = depJet.depositAndCountWithin(0.3).second;
01179 isoR03.trackerVetoPt = depTrk.candEnergy();
01180 isoR03.emVetoEt = depEcal.candEnergy();
01181 isoR03.hadVetoEt = depHcal.candEnergy();
01182 isoR03.hoVetoEt = depHo.candEnergy();
01183
01184 isoR05.sumPt = depTrk.depositWithin(0.5);
01185 isoR05.emEt = depEcal.depositWithin(0.5);
01186 isoR05.hadEt = depHcal.depositWithin(0.5);
01187 isoR05.hoEt = depHo.depositWithin(0.5);
01188 isoR05.nTracks = depTrk.depositAndCountWithin(0.5).second;
01189 isoR05.nJets = depJet.depositAndCountWithin(0.5).second;
01190 isoR05.trackerVetoPt = depTrk.candEnergy();
01191 isoR05.emVetoEt = depEcal.candEnergy();
01192 isoR05.hadVetoEt = depHcal.candEnergy();
01193 isoR05.hoVetoEt = depHo.candEnergy();
01194
01195
01196 aMuon.setIsolation(isoR03, isoR05);
01197
01198 }
01199
01200 reco::Muon MuonIdProducer::makeMuon( const reco::Track& track )
01201 {
01202
01203 double energy = sqrt(track.p() * track.p() + 0.011163691);
01204 math::XYZTLorentzVector p4(track.px(),
01205 track.py(),
01206 track.pz(),
01207 energy);
01208 return reco::Muon( track.charge(), p4, track.vertex() );
01209 }
01210
01211 double MuonIdProducer::sectorPhi( const DetId& id )
01212 {
01213 double phi = 0;
01214 if( id.subdetId() == MuonSubdetId::DT ) {
01215 DTChamberId muonId(id.rawId());
01216 if ( muonId.sector() <= 12 )
01217 phi = (muonId.sector()-1)/6.*M_PI;
01218 if ( muonId.sector() == 13 ) phi = 3/6.*M_PI;
01219 if ( muonId.sector() == 14 ) phi = 9/6.*M_PI;
01220 }
01221 if( id.subdetId() == MuonSubdetId::CSC ) {
01222 CSCDetId muonId(id.rawId());
01223 phi = M_PI/4+(muonId.triggerSector()-1)/3.*M_PI;
01224 }
01225 if ( phi > M_PI ) phi -= 2*M_PI;
01226 return phi;
01227 }
01228
01229 double MuonIdProducer::phiOfMuonIneteractionRegion( const reco::Muon& muon ) const
01230 {
01231 if ( muon.isStandAloneMuon() ) return muon.standAloneMuon()->innerPosition().phi();
01232
01233 if ( muon.matches().empty() ){
01234 if ( muon.innerTrack().isAvailable() &&
01235 muon.innerTrack()->extra().isAvailable() )
01236 return muon.innerTrack()->outerPosition().phi();
01237 else
01238 return muon.phi();
01239 }
01240 return sectorPhi(muon.matches().at(0).id);
01241 }
01242
01243 void MuonIdProducer::fillGlbQuality(edm::Event& iEvent, const edm::EventSetup& iSetup, reco::Muon& aMuon)
01244 {
01245 edm::Handle<edm::ValueMap<reco::MuonQuality> > glbQualH;
01246 iEvent.getByLabel(globalTrackQualityInputTag_, glbQualH);
01247
01248 if(aMuon.isGlobalMuon() && glbQualH.isValid() && !glbQualH.failedToGet()) {
01249 aMuon.setCombinedQuality((*glbQualH)[aMuon.combinedMuon()]);
01250 }
01251
01252 LogDebug("MuonIdentification") << "tkChiVal " << aMuon.combinedQuality().trkRelChi2;
01253
01254 }
01255
01256 void MuonIdProducer::fillTrackerKink( reco::Muon& aMuon ) {
01257
01258 if (aMuon.innerTrack().isNull()) return;
01259
01260 reco::MuonQuality quality = (aMuon.isQualityValid() ? aMuon.combinedQuality() : reco::MuonQuality());
01261
01262 bool filled = trackerKinkFinder_->fillTrkKink(quality, *aMuon.innerTrack());
01263
01264 if (filled || aMuon.isQualityValid()) aMuon.setCombinedQuality(quality);
01265 }
01266
01267 bool MuonIdProducer::checkLinks(const reco::MuonTrackLinks* links) const {
01268 bool trackBAD = links->trackerTrack().isNull();
01269 bool staBAD = links->standAloneTrack().isNull();
01270 bool glbBAD = links->globalTrack().isNull();
01271 if (trackBAD || staBAD || glbBAD )
01272 {
01273 edm::LogWarning("muonIDbadLinks") << "Global muon links to constituent tracks are invalid: trkBad "
01274 <<trackBAD <<" standaloneBad "<<staBAD<<" globalBad "<<glbBAD
01275 <<". There should be no such object. Muon is skipped.";
01276 return false;
01277 }
01278 return true;
01279 }
1.7.3