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