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