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