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