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